基于相移偏移和匹配追蹤的探地雷達成像方法

張建中 黃月琴

（1.中國海洋大學海洋地球科學學院，山東 青島266100；2.廈門大學物理與機電工程學院，福建 廈門361005）

引 言

偏移成像是探地雷達（GPR）目標識別和定位的主要技術之一。相移偏移方法因簡單實用被廣泛用于GPR成像［1－4］，這類方法中的一個重要步驟是用反傅里葉變換把波數域的相移數據變回到空間域。但是，由于均勻采樣頻率與波數之間的非線性關系，使得波數空間的波場被非均勻采樣，常規的快速傅里葉變換（FFT）難以直接使用。現有的做法是先對波數空間非均勻采樣的波場進行插值運算，獲得均勻采樣的波數空間波場，再直接使用FFT，如線性插值－FFT 算法［5］，Stolt插 值－FFT 算 法［6－7］等。雖然FFT運行很快，但插值的計算量卻較大，而且，插值后的波場將不可避免地含有一定的誤差，從而影響成像精度。Liu等人［8］提出了精確的非均勻采樣快速傅里葉變換算法（NUFFT）。把NUFFT引入到相移算法中，可以避免插值運算，使成像精度和計算效率得到進一步的提高。

實測GPR數據通常被噪聲所污染，致使偏移成像質量和目標定位精度受到嚴重影響。匹配追蹤（MP）算法［9］是目前最流行的信號分析工具之一，被用于地震和雷達信號處理中［10－12］。MP算法選擇一個超完備展開函數集合，根據信號特征，自適應地選擇展開函數，通過逐步逼近，求得信號的表示，達到分解信號和去除噪聲的效果。但是，MP算法在每一次求解的分解過程中，需要在超完備展開函數集合或原子庫中尋找最優原子，這是一個非常耗時的過程，導致MP算法具有很大的計算量，從而限制了該算法的實際應用。為此，將最小二乘原理引入到MP算法中，直接計算最優原子參數，避免了在龐大的原子庫中搜尋最優原子的耗時過程，形成了具有較高計算速度的最小二乘匹配追蹤算法（LSMP）。

根據相移偏移方法和MP算法的上述問題和改進途徑，對它們進行了改進，并把它們綜合起來，對GPR資料進行成像處理，獲得了偏移結果的振幅譜和相位譜，顯示出了良好的成像效果。

1 改進的相移偏移方法

相移偏移算法屬于頻率－波數域方法，是目前常用的GPR偏移成像方法之一。該方法所要解決的問題可描述為：利用地表z＝0處接收的GPR信號s（x，z＝0，t）來重建波場s（x，z，t＝0）.其方法原理如下：

1）對GPR回波信號s（x，z＝0，t）進行二維傅里葉變換，得

2）在頻域－波數空間，任意深度z處的波場S（kx，z，ω）可由z＝0處的波場表示為

式中kz是z方向波數。kx，kz和ω滿足

式中ν是電磁波在地下的平均傳播速度。

3）對式（2）進行二維反傅里葉變換，可得時域空間內任意深度z處的波場為

4）將式（2）代入在式（4），并令t＝0，就得到偏移后的波場s（x，z，t＝0）為

由于GPR系統通常在時間t方向和測線x方向上等間距采樣，式（1）中傅里葉變換后的頻率ω和kx波數都為均勻分布，由式（3）的非線性關系式得到的波數kz則是非均勻分布的。這樣，式（5）右端實際上是一個非均勻采樣的傅里葉逆變換式。為了避免因使用FFT而需要把波數空間非均勻采樣數據插值成均勻采樣的問題，我們采用Liu等人提出的NUFFT［8］進行快速計算。式（5）的離散形式為［13］

式中：Δkx和Δω為頻域－波數域的采樣率；M和N分別是x和t方向上的樣點數；xp和zq是樣點坐標。令

則有

2 最小二乘匹配追蹤算法

假定GPR某道回波信號u（t）可以表示成子波w（t）經過不同的時移和尺度變換后的子波的線性組合，

式中：aj、tj、fj和φj分別表示第j個子波的振幅、中心時間、主頻和相位，j＝1，2，…，N；n是噪聲。由于Ricker子波可以很好地表征GPR脈沖信號，采用Ricker子波，其表達式為

為了能夠得到回波信號的振幅譜和相位譜，采用解析信號進行處理。利用Hilbert變換，GPR信號u（t）和小波w（t）的解析信號分別為

式中H表示Hilbert變換。這樣，式（9）的GPR回波解析道信號可表示為

式中Aj為復振幅，且

每一次匹配追蹤后的解析道信號殘差R（t）為［14］

根據最小二乘原理，解析道信號殘差能量定義為

式中，ti為回波信號的時間樣點，i＝1，2，… ，M.為使殘差能量E達到最小，解式（16）極小值問題得到

I是單位矩陣，ε是阻尼因子，T表示矩陣轉置。

利用式（17）就可求出每一道回波信號的解析子波復振幅。具體計算步驟如下：

1）計算中心時間tj，即輸入信號U（t）振幅包絡的峰值所對應的時間。

2）計算主頻fj.采取窗函數截取中心時間tj的局部波形，并利用FFT計算瞬時頻譜。該頻譜振幅包絡的峰值所對應的頻率即為所需要的主頻。

3）將中心時間tj和主頻fj代入式（10），并利用式（12）計算當前解析子波。

4）利用式（17）計算出復數振幅Aj.

5）由式（15）計算信號殘差R（t），且令

U（t）＝R（t）.

6）計算殘差能量E，并判斷是否滿足終止條件，若是，結束；否則，轉至步驟1）。

執行上述方法步驟，就得到了組成GPR某一道隨時間變化信號的各個子波的復振幅Aj.這些子波的線性疊加就近似為有效信號，而殘差則被認為是噪聲。去噪后的解析道信號用獲得的解析子波表示為

該信號的幅度譜和相位譜可分別表示為

式中Re和Im分別表示取復函數的實部和虛部。

3 實驗結果

分別使用數值模擬獲得的GPR信號及實測的GPR資料進行實驗，主要測試本文算法的去噪和聚焦性能。為了將相移偏移結果與本文的相移偏移和LSMP綜合算法結果進行一一對應的比較，對相移偏移結果也通過Hilbert變換形成對應的解析信號后，再利用與式（19）和（20）類似的方法計算出相應的振幅譜和相位譜。

3.1 合成數據實驗

圖1（a）所示為一GPR理論模型，長度為1.4m，厚度為0.6m.背景媒質為干燥土壤，深度從0.05m到0.60m，介電常數εr為10；土壤上層為空氣；土壤中分布有一個混凝土板，該板在中間部位垂直斷裂而錯位，板的總長度為1.4m，厚度為0.1m，介電常數εr為6.利用軟件 Gpr MaxV2.0［15］模擬了該模型的 GPR記錄，其中，雷達脈沖主頻為900MHz.圖1（b）是消除了直達波和地面反射波的模擬記錄，在界面突變處有明顯的繞射波，可以用來檢驗偏移方法的歸位效果。圖1（c）和（d）是在圖1（b）模擬記錄中加入隨機噪聲后的含噪記錄，其中圖1（c）的信噪比SNR為0dB，回波信號還比較明顯，以下稱為含弱噪聲信號；圖1（d）的信噪比SNR降為－10dB，回波信號已完全淹沒在隨機噪聲中，以下稱為含強噪聲信號。為了與模型進行對比，所有圖中的時間坐標被換算成了深度坐標。

分別使用改進的相移偏移方法及改進的相移偏移與LSMP綜合方法對合成的無噪信號和含噪聲信號進行了成像處理。圖2（a）和（b）是無噪信號偏移結果的幅度譜和相位譜，斷裂處的繞射波被較好的歸位，說明了使用NUFFT相移偏移方法的良好效果。圖2（c）和（d）是對信噪比SNR為0dB的含弱噪聲信號，用相移偏移方法的幅度譜和相位譜，而圖2（e）和（f）是對該含弱噪聲信號用相移偏移與LSMP綜合方法得到的幅度譜和相位譜。可以看出，只用相移偏移方法的結果雖然尚能反映目標，但背景噪聲較大；而相移偏移與LSMP綜合方法能很好的壓制噪聲，成像效果得到了明顯改善。圖2（g）和（h）是對信噪比SNR為－10dB的含強噪聲信號，用相移偏移方法得到的幅度譜和相位譜，而圖2（i）和（j）是對該含強噪聲信號用相移偏移與LSMP綜合方法得到的幅度譜和相位譜。由于目標信號幾乎被強噪聲淹沒，在相移偏移圖上噪聲依然很強，而目標信號非常微弱，幾乎不能識別出來；而在相移偏移與LSMP綜合成像圖上，噪聲得到壓制，目標信號反映清晰。可見，相移偏移與LSMP綜合方法明顯優于相移偏移方法，不僅能夠對同相軸正確歸位，更具有很強的壓制噪聲的能力，從而極大地提高了低信噪比資料的成像質量，較好地刻畫目標體的位置和形狀。另外，除了振幅譜外，相位譜也能較好地反映目標信息。由于相位譜反映信號相位變化的連續性，與反射信號能量無關，對于同相軸的連續展布以及突變，相位譜的特征更明顯。特別是對于含噪聲資料，相位譜比振幅譜能更明顯地突出目標信號，因此，將幅度譜和相位譜進行綜合分析，可以更好地識別和圈定目標。

3.2 實測數據實驗

圖3（a）是我們用pulseEKKO PRO型探地雷達在英國東約克郡的Beverley高爾夫球場采集到的雷達數據剖面圖。該雷達系統采用Ricker脈沖，工作頻率為500MHz，脈沖電壓為180 V.共采集540道數據，道間距為2cm，每道記錄250個采樣點，采樣時窗為50ns.通過試驗得到電磁波在該球場土壤中的傳播速度約為0.09m／ns.在進行偏移前，對該雷達數據進行了去除背景信號及相關噪聲處理，結果如圖3（b）所示。用相移偏移算法得到的幅度譜和相位譜如圖3（c）和（d）所示；用相移偏移與LSMP綜合方法進行成像，得到的幅度譜和相位譜分別如圖3（e）和（f）所示。可以看出，后一種方法成像結果的背景噪聲更弱，同相軸聚焦更好，成像質量明顯優于前一種方法。

4 結 論

最小二乘匹配追蹤算法具有很強的壓制隨機噪聲的能力，引入NUFFT的相移方法具有較高的精度和計算效率。綜合最小二乘匹配追蹤算法與相移偏移算法，形成了一種壓噪能力很強的快速GPR偏移成像方法，極大地提高了含噪聲資料的成像質量。另外，成像結果的相位譜也能較好地反映目標信息。特別是對于含噪聲資料，有時相位譜比振幅譜對目標信息的反映更明顯。因此，將幅度譜和相位譜進行綜合分析，可以更好地識別和圈定目標，這對于低信噪比資料的解釋和弱小目標的識別具有重要意義。

［1］CAFFORIO C，PRATI C，ROCCA F.SAR data focusing using seismic migration techniques［J］.IEEE Transactions on Aerospace and Electronic Systems，1991，27（2）：194－207.

［2］LAMBOT S，SLOB E，VAN DEN BOSCH I，et al.Modeling of ground－penetrating radar for accurate characterization of the subsurface dielectric properties［J］.IEEE Transactions on Geoscience and Remote Sensing，2004，42（11）：2555－2568.

［3］ODEN C P，POWERS M H，WRIGHT D L，et al.Improving GPR image resolution in lossy ground using dispersive migration［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（8）：2492－2500.

［4］修志杰，陳 潔，方廣有，等.基于F－K偏移及最小熵技術的探地雷達成像法［J］.電子與信息學報，2007，29（4）：827－830.XIU Zhijie，CHEN Jie，FANG Guangyou，et al.Ground penetrating radar imaging based on F－K migration and minimum entropy method［J］.Journal of Electronics ＆Information Technology，2007，29（4）：827－830.（in Chinese）

［5］GAZDAG J，SGUAZZERO P.Migration of seismic data by phase shift plus interpolation［J］.Geophysics，1984，49（2）：124－131.

［6］STOLT R H.Migration by Fourier transform［J］.Geophysics，1978，42（1）：23－48.

［7］張春城，周正歐.基于Stolt偏移的探地雷達合成孔徑成像研究［J］.電波科學學報，2004，19（3）：316－320.ZHANG Chuncheng，ZHOU Zhengou.Ground penetrating radar synthetic aperture imaging based on Stolt migration［J］.Chinese Journal of Radar Science，2004，19（3）：316－320.（in Chinese）

［8］LIU Q H，NGUYEN N.An accurate algorithm for nonuniform fast Fourier transform （NUFFT’s）［J］.IEEE Microwave and Guided Wave Letters，1998，8（1）：18－20.

［9］MALLAT S，ZHANG Z.Matching pursuits with time－frequency dictionaries［J］.IEEE Transaction on signal processing，1993，41（12）：3377－3415.

［10］WANG Y.Seismic time－frequency spectral decomposition by matching pursuit［J］.Geophysics，2007，72（1）：V13－V20.

［11］朱 明，普運偉，金煒東，等.基于時頻原子方法的雷達輻射源信號特征提取［J］.電波科學學報，2007，22（3）：458－462.ZHU Ming，PU Yunwei，JIN Weidong，et al.Feature extraction of radar emitter signals based on timefrequency atoms［J］.Chinese Journal of Radar Science，2007，22（3）：458－462.（in Chinese）

［12］孫閩紅，唐 斌.基于原子分解理論的雷達欺騙式干擾信號特征提取.電波科學學報，2008，23（3）：550－554.SUN Minhong，TANG Bin.Teature extraction of radar deceptive－jamming signal based on atomic decomposition［J］.Chinese Journal of Radar Science，2008，23（3）：550－554.（in Chinese）

［13］HUANG Y Q，LIU Y H，Liu Q H，ZHANG J Z.Improved 3－D GPR detection by NUFFT combined with MPD method［J］.Progress in Electromagnetic Research，2010，103：185－199.

［14］LIU J，MARFURT K J.Instantaneous spectral attributes to detect channels［J］.Geophysics，2007，72（2）：23－31.

［15］GIANNOPOULOS A.GprMax2DV2.0 （Electromagnetic simulator for ground probing radar）［CP／OL］. 2005－10－31 ［2010－06－08］. http：／／www.gprmax.org／download.html.