基于曲波域凸集投影算法的缺失GPR 信號高精度重建

吳祺銘，王洪華，席宇何，王欲成

(桂林理工大學 地球科學學院，廣西 桂林 541004)

探地雷達(Ground Penetrating Radar,GPR)是一種利用高頻電磁波來探測地下淺部地質結構的地球物理方法[1-3]，具有操作便捷、分辨高、效率高、無損等優(yōu)點，被廣泛應用于路基病害、塌陷空洞、隧道襯砌等淺層目標體的精細探測[4-7]。高質量GPR 信號是利用成像方法精確分析淺層物性和幾何結構的前提。然而，受探測現場河溝、施工機器等障礙物、起伏地形引起的天線耦合不佳、儀器性能等影響，實測GPR 剖面中不可避免會出現信號缺失和壞道，易造成目標體產生反射波和繞射波同相軸斷開現象[8]。常規(guī)處理方法如線性插值法根據相鄰道之間的信號連續(xù)進行重建，當信號大規(guī)模連續(xù)缺失時，重建誤差較大，嚴重影響后續(xù)偏移與反演成像效果。因此，深入研究缺失信號的高精度重建方法，是GPR 研究的熱點之一。

從信號處理角度看，利用低于奈奎斯特(Nyquist)采樣定理的未缺失信號去重建缺失信號是一個欠定問題，若先驗信息未知，通常難以高精度重建缺失信號[9]。近年來，壓縮感知理論為解決缺失信號的高精度重建問題提供了一種有效方法[10-12]。壓縮感知理論指出，信號在某個變換域內可稀疏表示條件下，通過采用特定的采樣方式和合適的稀疏變換系數迭代方式求解，可有效實現缺失信號的高精度重建[13-16]，其不需背景介質和目標體的先驗信息，且兼具計算精度高和效率高的優(yōu)點[17]。高精度重建缺失GPR 信號的關鍵是對其稀疏表示，通常GPR 信號呈時空域分布，不具備稀疏性。為此，傅里葉變換[18-21]、小波變換[22-24]、曲波變換[25-26]等稀疏變換被引入到GPR 信號處理領域并得到廣泛應用[17]。上述稀疏變換方法由于基函數的不同，稀疏表示GPR 信號的性能各有優(yōu)劣[27-28]，如傅里葉變換僅適合表征線性GPR 信號，對彎曲GPR 信號同相軸表征能力較弱。曲波變換(Curvelet Transform)具有各向異性、局部性、多尺度等特性，能更好地表征GPR 信號的非線性和多尺度特征，被認為是當前GPR 信號的最優(yōu)稀疏表示方法之一[29-30]。

凸集投影(Projection Onto Convex Sets,POCS)算法是壓縮感知理論下基于稀疏變換的重建方法之一[31-32]，其基本原理是在根據信號屬性構造不同凸集約束集合的基礎上，對目標函數進行迭代求解，以獲取滿足所有約束集合的解，從而消除因信號丟失導致的不相干噪聲，高精度重建缺失信號[33]。相比于迭代閾值法、Bregman迭代法等重建方法，POCS 算法因其具有理論簡單，易實現，計算效率高的優(yōu)點而在圖像處理領域[34-37]和地球物理信號重建領域[38-42]得到廣泛應用。與此同時，M.Sato 等[43]應用頻率域POCS 算法實現了缺失GPR信號的重建；Li Yi 等[44]利用該算法分析了不同GPR信號缺失比例和閾值迭代模型對缺失信號重建精度的影響。然而，頻率域POCS 算法采用的傅里葉變換僅能展現信號的全局特征，未考慮地球物理信號的非線性和多尺度特征，重建精度較低[45]。為此，一些學者將具有各向異性、多方向和多尺度特征的曲波變換代替傅里葉變換，提出了曲波域POCS 算法并應用于缺失地震信號的高精度重建，數值試驗結果表明，曲波變換能更好地稀疏表示非線性、多尺度、多方向的地震信號，更有利于提高缺失信號重建精度[46-50]。

為此，筆者在上述研究基礎上，為克服傅里葉變換稀疏表示的不足，將曲波變換與POCS 算法相結合，提出一種基于曲波域POCS 算法的缺失GPR 信號高精度重建方法。通過模擬和實測缺失GPR 信號的頻率域和曲波域POCS 重建試驗，以分析曲波域POCS 重建方法的有效性及其在提高重建精度方面的優(yōu)勢。

1 基于曲波域凸集投影算法的缺失GPR 信號重建方法

根據壓縮感知理論，實測部分缺失GPR 信號dobs與完整信號d的稀疏采樣關系[51-52]表示為：

考慮到曲波變換的非線性、多方向、多尺度特征，且曲波基類似于GPR 信號的曲線特征，是非線性信號最優(yōu)稀疏表示的方法之一[53-54]，本文采用曲波變換對GPR 信號進行稀疏表示。

曲波變換通過計算曲波基與信號的內積來實現信號的稀疏[55]表示：

因此，利用曲波變換作為稀疏變換，完整信號d可表示為：

根據曲波變換和最小二乘原理，可構建缺失GPR信號重建的目標函數為：

利用信號的曲波變換和曲波系數迭代閾值模型分別構建的凸集約束集合[45]為：

利用dobs的曲波變換系數作為初始解在2 個凸集約束集合內進行多次迭代投影，獲得的滿足約束條件的解[56]可表示為：

dobs在曲波域內的迭代更新可看作凸集Q1的不斷迭代投影過程[57]，即：

集合Q2的投影過程是利用閾值算子將小于閾值的曲波系數置零，保留大于閾值的曲波系數，并作為第k次迭代的初始解：

本文分別采用線性閾值和指數閾值模型重建缺失的GPR 信號[38]。

線性閾值模型可表示為：

指數閾值模型可表示為：

將式(8)和式(9)以及閾值算子Tk合并運算[45]，可得到：

由式(13)可知：在每次迭代時減小閾值Tk可得到具有信號局部特征的曲波變換系數，并不斷更新變換系數xk-1，即可獲得第k次迭代后的曲波變換系數。

將缺失位置處的重建信號加載到待重建信號dobs，則可獲得用于下次迭代的時域GPR 信號：

將式(14)代入式(15)，則第k+1次的曲波域POCS重建的時域迭代公式為：

為定量評價缺失信號的POCS 重建精度，本文計算缺失位置處每次迭代后的重建信號與完整信號的相對誤差公式[39]為：

此外，應用重建后的GPR 剖面與完整GPR 剖面之間的平均絕對誤差EMA(Mean Absolute Error,MAE)、信噪比RS/N(Signal to Noise Ratio,SNR)和峰值信噪比RPS/N(Peak Signal to Noise Ratio,PSNR)來進一步評價GPR 剖面的重建質量。平均絕對誤差、信噪比、峰值信噪比的計算公式[58-59]分別為：

根據上述POCS 重建理論，構建的缺失GPR 信號曲波域POCS 重建流程如圖1 所示。

圖1 缺失GPR 信號的曲波域POCS 重建流程Fig.1 Reconstruction process of missing GPR signals based on curvelet-domain POCS

2 數值算例

2.1 曲波域POCS 重建方法的優(yōu)勢分析

為驗證曲波域POCS 算法應用于缺失GPR 信號重建的有效性及其優(yōu)勢，利用GprMax 軟件對大小為3.0 m×2.0 m 的空洞模型進行計算，獲得的模擬GPR 剖面如圖2a 所示。模型背景介質的相對介電常數為4，電導率為0.005 S/m，模型正中心埋有一個半徑為0.03 m 的圓形空洞。由圖可見，空洞產生的雙曲線繞射波清晰可見，同相軸光滑、連續(xù)性較好。對圖2a 隨機缺失40%信號獲得的結果如圖2b 所示，在水平位置1.0 m 附近出現多道連續(xù)缺失，部分信號缺失后的直達波和繞射波同相軸連續(xù)性較差，雙曲線繞射波信號不完整，嚴重干擾其識別與解釋。

圖2 空洞模型的GPR 模擬剖面Fig.2 The simulated GPR profile of cavity model

圖3 是應用不同閾值模型的POCS 算法對圖2b 所示的GPR 信號進行重建獲得的GPR 剖面。其中，圖3a和圖3b 分別為線性和指數閾值模型的頻率域POCS 重建剖面；圖3c 和圖3d 分別為線性和指數閾值模型的曲波域POCS 重建剖面。由圖可見：頻率域和曲波域POCS 重建方法均能較好地重建缺失位置處的直達波；但對比4 個分圖可知，不同閾值模型和稀疏域的POCS算法重建缺失處繞射波的能力有顯著不同，圖3a、圖3b 中水平位置約1.0 m 處出現能量較強的縱向偽影，如黑色箭頭所示。與圖3a、圖3b 相比，圖3c 展示的線性閾值模型的曲波域POCS 重建剖面中雖然縱向偽影能量較弱，但繞射波同相軸的兩側出現明顯斷開現象，如黑色箭頭所示；圖3d 所示的指數閾值模型的曲波域POCS 重建剖面與圖2b 所示的原始完整剖面非常吻合，重建后的雙曲線繞射波同相軸完整、連續(xù)性較好且未見明顯的縱向偽影。

圖3 頻率域和曲波域POCS 重建剖面Fig.3 Reconstructed profiles by using POCS in frequency domain and curvelet domain

圖4 為圖3 中各重建GPR 剖面的第103 道(水平位置1.02 m)的波形對比。由圖可見，線性閾值模型、指數閾值模型的頻率域POCS 和線性閾值模型的曲波域POCS 重建信號與原始信號相比出現不同程度的噪聲，重建精度較差，而指數閾值模型的曲波域POCS 重建信號與原始信號擬合較好，這個現象也可從15.6 ns附近處的波形放大圖中觀察到。由此可知，指數閾值模型的曲波域POCS 重建精度更高。

圖4 圖3 中各剖面的第103 道(水平位置1.02 m)波形對比Fig.4 Waveform comparison of the channel 103st (horizontal position 1.02 m) in each profile in Fig.3

為進一步分析曲波域POCS 重建精度和質量，將圖2a 所示的完整剖面與圖3 所示的重建剖面做差，獲得的結果如圖5 所示。圖5a 和圖5b 中僅在水平位置1.0 m 附近出現明顯的重建誤差；而線性閾值模型的曲波域POCS 重建誤差主要集中于雙曲線繞射波位置，如圖5c 所示。與圖5a-圖5c 相比，圖5d 中的重建誤差非常微弱，幾乎不可見。由此可見，與其他重建方法相比，指數閾值模型的曲波域POCS 重建誤差更小，精度更高。

圖5 空洞模型的GPR 模擬剖面與頻率域和曲波域POCS 重建剖面(圖3)誤差Fig.5 The error between the GPR simulation profiles of the cavity model and reconstructed profiles (Fig.3) by using POCS in frequency domain and curvelet domain

利用式(17)計算的4 種POCS 重建方法的迭代相對誤差曲線如圖6 所示。由圖可知，在前25 次重建迭代過程中指數閾值模型的頻率域POCS 重建誤差下降最快，經40 次迭代后相對誤差基本穩(wěn)定在2.81×10-2%；而線性閾值模型的頻率域POCS 重建迭代誤差緩慢下降，100 次迭代后相對誤差基本穩(wěn)定在6.18×10-2%。線性閾值模型的曲波域POCS 重建誤差下降最慢，100 次迭代后誤差最高，如圖中綠線所示。指數閾值模型的曲波域POCS 重建誤差下降速度較快，經80 次迭代后相對誤差基本穩(wěn)定在5.21×10-4%，如圖中紅線所示。由此可見，相比于其他POCS 重建方法，指數閾值模型的曲波域POCS 算法迭代誤差最小，重建精度更高。

圖6 不同重建方法的迭代相對誤差Fig.6 Iterative relative errors of different reconstruction methods

利用式(18)-式(20)計算圖3 各分圖與圖2a 的平均絕對誤差EMA、信噪比RS/N、峰值信噪比RPS/N，見表1。與其他3 種POCS 重建方法的相比，指數閾值模型的曲波域POCS 重建后的平均絕對誤差下降86%～99%，平均絕對誤差最低；信噪比和峰值信噪比提高9～20 dB，這充分說明了指數閾值模型的曲波域POCS 算法的重建精度更高，平均絕對誤差更小，信噪比與峰值信噪比更高。為此，后續(xù)算例均采用指數閾值模型的曲波域POCS 算法進行計算。

表1 不同POCS 重建方法的重建誤差對比Table 1 Error comparison of different POCS-based reconstruction methods

2.2 復雜模型

為進一步分析本文提出的曲波域POCS 算法用于重建復雜地質結構GPR 缺失信號的有效性，利用GprMax 軟件對如圖7a 所示的復雜地電模型進行計算，獲得的模擬GPR 剖面如圖7b 所示。模型被起伏界面和水平界面分成上、中、下3 層，其相對介電常數分別為6、12、15，電導率分別為0.001、0.020、0.050 S/m，在第2 層的左右兩側分別埋有一個圓形空洞與長方形空洞。由圖7b 可見，起伏界面與水平界面產生的反射波和圓形空洞、長方形空洞產生的繞射波同相軸光滑連續(xù)。

圖7 GPR 模型及其正演剖面Fig.7 GPR model diagram and forward modeling profile

對圖7b 所示的完整GPR 剖面隨機缺失40%，50%，60%，70%的結果如圖8 所示。由圖可見，GPR 信號隨機缺失比例越大，繞射波與反射波同相軸連續(xù)性越差；缺失部分信號后，空洞、起伏界面與水平界面產生的繞射波和反射波同相軸不完整，有效信息被嚴重缺失，影響后續(xù)偏移與反演成像效果。

圖8 圖7b 中隨機缺失不同占比后的GPR 剖面Fig.8 GPR profiles with random signal missing proportions (derived from Fig.7b)

圖9 為采用頻率域POCS 重建方法對圖8 中缺失GPR 信號進行重建獲得的GPR 剖面。由圖9 可見，對不同比例隨機缺失GPR 信號重建時，起伏界面與水平界面產生的反射波與空洞產生的雙曲線繞射波大部分信號都能得到較好的重建，重建后同相軸較為連續(xù)，僅在多段連續(xù)缺失處存在斷開現象。對比圖9 中4 個分圖可見，信號隨機缺失比例越大，重建后的GPR 剖面中雜波和縱向偽影能量越強，連續(xù)多道缺失處重建效果越差，重建精度越低。

圖9 圖8 經頻率域POCS 重建后的GPR 剖面Fig.9 GPR reconstructed profiles based on frequency-domain POCS (derived from Fig.8)

圖10 是采用曲波域POCS 重建方法對圖8 中的缺失GPR 剖面進行重建獲得的GPR 剖面。由圖10 可見，隨機缺失40%和50%的重建剖面中直達波、反射波和繞射波同相軸均非常連續(xù)，未出現明顯斷開現象，缺失位置處的繞射波與反射波能量均得到有效恢復，重建精度較高；當隨機缺失比例為60%和70%時，重建剖面中直達波和界面反射波同相軸也較為光滑，僅在連續(xù)多道缺失處出現未重建完全的現象。與圖9 所示的頻率域POCS 重建結果相比，曲波域POCS 重建后GPR 剖面中直達波、反射波和繞射波同相軸連續(xù)性更好，雜波能量更小。

利用圖9 和圖10 所示GPR 信號重建剖面計算的迭代相對誤差曲線如圖11 所示。由圖可見，4 種隨機缺失比例條件下頻率域POCS 重建誤差均大于0.3%，且重建誤差隨缺失比例的增大而增大；而曲波域POCS重建誤差都小于0.2%，當缺失比例為70%，重建誤差僅為0.14%。由此可見，但本文所提出的POCS 算法的重建誤差均相對較小，可有效用于大比例隨機缺失的GPR 信號的高精度重建。

圖11 GPR 信號缺失的重建迭代誤差Fig.11 Iterative reconstruction errors of missing GPR signals

將圖7b 所示的完整剖面與圖9 中重建剖面做差獲得的誤差剖面如圖12 所示。由圖可見，頻率域POCS重建精度不高，重建誤差較大，且缺失比例越大，整個剖面的誤差越大，如圖12 所示。將圖7b 所示的完整剖面與圖10 中重建剖面做差獲得的誤差剖面如圖13 所示。相較于頻率域POCS 算法，曲波域POCS 算法在相同缺失比例下重建誤差更小，重建誤差主要分布在信號連續(xù)多道缺失位置處，且未見明顯的縱向偽影，可見曲波域POCS 重建精度更高，重建誤差更小。

圖12 頻率域POCS 重建后的誤差剖面Fig.12 Reconstruction error profiles based on frequency-domain POCS

圖13 曲波域POCS 重建后的誤差剖面Fig.13 Reconstruction error profiles based on curvelet-domain POCS

計算圖9 和圖10 中各分圖與圖7b 的平均絕對誤差EMA、信噪比RS/N、峰值信噪比RPS/N見表2。缺失比例越大，頻率域POCS 與曲波域POCS 重建后的平均絕對誤差越大，信噪比與峰值信噪比越低；與頻率域POCS 重建相比，相同缺失比例下的曲波域POCS 重建后的平均絕對誤差下降82%～96%，且信噪比和峰值信噪比提高6～17 dB?？梢?，本文提出的曲波域POCS 算法可高精度重建復雜地質結構的缺失GPR 信號。

表2 不同缺失比例和不同POCS 重建方法的重建誤差對比Table 2 Error comparison of different POCS-based reconstruction methods with different signal missing proportions

3 實測算例

為分析曲波域POCS 重建算法的實際效果，對桂林某道路GPR 信號剖面進行重建測試。圖14 為經過零時校正、濾波、增益等處理的實測剖面。由圖可見，直達波、反射波和繞射波能量較強，且其同相軸連續(xù)對圖14 所示的完整GPR 剖面隨機缺失40%和70%，獲得的結果如圖15 所示。與圖14 相比，隨機缺失后的直達波與反射波同相軸出現明顯斷開現象，嚴重破壞了其波形特征。

圖14 某道路結構探測的GPR 剖面Fig.14 GPR profile of a road structure

圖15 圖14 中隨機缺失信號后的GPR 剖面Fig.15 GPR profiles with random signal missing proportions (derived from Fig.14)

圖16 是應用不同稀疏域對圖15 中不同隨機缺失比例的GPR 信號重建獲得的GPR 剖面圖。由圖16a、圖16c 可見，當缺失比例為40%時，2 種重建方法都能有效地重建反射波和繞射波能量，重建精度較高。但當GPR 信號大比例缺失時，頻率域POCS 重建剖面中出現了多處縱向偽影，如圖16b 所示。而由圖16d 可見，在70%隨機缺失比例下曲波域POCS 重建效果較好，未見明顯的縱向偽影。

圖16 圖15 經頻率域和曲波域重建后的GPR 剖面Fig.16 GPR reconstructed profiles in frequency and curvelet domains (derived from Fig.15)

將圖14 所示的完整剖面與圖16 所示的重建剖面做差，獲得的誤差剖面如圖17 所示。由圖可見，當隨機缺失比例為40%時，頻率域POCS 與曲波域POCS 重建誤差較小，僅出現在連續(xù)多道缺失處，且曲波域POCS 重建誤差幾乎不可見；當隨機缺失比例為70%時，頻率域POCS 重建誤差廣泛分布在整個剖面中，而曲波域POCS 重建誤差僅在連續(xù)多道缺失處存在較弱的能量。計算圖16 各分圖與圖14 獲得的平均絕對誤差EMA、信噪比RS/N、峰值信噪比RPS/N見表3，曲波域POCS 重建后的平均絕對誤差相比于頻率域POCS 重建誤差下降45%～79%；信噪比和峰值信噪比提高1～9dB，且缺失GPR 信號比例越多，其信噪比提升越大。由圖可見，本文提出的基于曲波域POCS 算法的GPR信號重建方法可有效用于實測GPR 信號缺失的高精度重建，且對信號大比例缺失具有較好的適應性。

圖17 實測剖面(圖14)與頻率域和曲波域POCS 重建剖面(圖16)的誤差Fig.17 Errors between the measured profile (Fig.14) and GPR reconstructed profiles based on POCS in frequency and curvelet domains(Fig.16)

4 結論

a.針對連續(xù)多道缺失GPR 信號的高精度重建問題，本文將非線性、多方向、多尺度的曲波變換與POCS 算法相結合，基于壓縮感知理論推導了缺失GPR 信號曲波域POCS 重建迭代的相關公式，提出了一種基于曲波域POCS 算法的缺失GPR 信號高精度重建方法。

b.模擬與實測GPR 信號試驗結果表明：POCS 算法可有效重建GPR 剖面中的缺失信號；與線性和指數閾值模型的頻率域POCS 算法、線性閾值模型的曲波域POCS 算法相比，指數閾值模型的曲波域POCS 重建方法的重建精度更高、平均絕對誤差下降45%～99%，信噪比和峰值信噪比提高1～20 dB，其重建結果可為后續(xù)處理與解釋提供高質量GPR 信號。

符號注釋：

a1為所在行數；a2為所在列數；b為離散采樣間隔；C為曲波正變換；C-1為曲波反變換；d為完整信號；dobs為缺失信號；drec為重建后信號；EMA為平均絕對誤差；f為時間域信號；為頻率域信號；I為單位矩陣；j為尺度參數；k為迭代次數；l為角度參數；m1為信號采樣點數；m2為信號道數；n為最大采樣頻率；N為最大迭代次數；p為頻移參數；p1為最小頻移參數；p2為最大頻移參數；Q1為曲波變換凸集集合；Q2為閾值集合；R是由對角元素為1 和0 的稀疏采樣矩陣，1 和0 分別代表未缺失信號與缺失信號；RS/N為信噪比；RPS/N為峰值信噪比；Sθl為剪切矩陣；Tk為閾值算子；t為時間；t1為時域離散信號中最小時刻；t2為時域離散信號中最大時刻；W為連續(xù)曲波系數；WD為離散曲波系數；x為完整信號估計值；為信號更新解；λ為閾值門限；λmax為曲波變換系數的最大值；β為凸集Q2中閾值；ω為角頻率；為離散曲波；的復共軛；Φ為目標函數；ε為迭代誤差；δ為極小值；()*為伴隨矩陣。