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

張建中 黃月琴

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

引 言

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

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

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

1 改進的相移偏移方法

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

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

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

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

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

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

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

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

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

則有

2 最小二乘匹配追蹤算法

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

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

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

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

式中Aj為復(fù)振幅，且

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

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

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

I是單位矩陣，ε是阻尼因子，T表示矩陣轉(zhuǎn)置。

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

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

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

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

4）利用式（17）計算出復(fù)數(shù)振幅Aj.

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

U（t）＝R（t）.

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

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

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

式中Re和Im分別表示取復(fù)函數(shù)的實部和虛部。

3 實驗結(jié)果

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

3.1 合成數(shù)據(jù)實驗

圖1（a）所示為一GPR理論模型，長度為1.4m，厚度為0.6m.背景媒質(zhì)為干燥土壤，深度從0.05m到0.60m，介電常數(shù)εr為10；土壤上層為空氣；土壤中分布有一個混凝土板，該板在中間部位垂直斷裂而錯位，板的總長度為1.4m，厚度為0.1m，介電常數(shù)ε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）是無噪信號偏移結(jié)果的幅度譜和相位譜，斷裂處的繞射波被較好的歸位，說明了使用NUFFT相移偏移方法的良好效果。圖2（c）和（d）是對信噪比SNR為0dB的含弱噪聲信號，用相移偏移方法的幅度譜和相位譜，而圖2（e）和（f）是對該含弱噪聲信號用相移偏移與LSMP綜合方法得到的幅度譜和相位譜。可以看出，只用相移偏移方法的結(jié)果雖然尚能反映目標，但背景噪聲較大；而相移偏移與LSMP綜合方法能很好的壓制噪聲，成像效果得到了明顯改善。圖2（g）和（h）是對信噪比SNR為－10dB的含強噪聲信號，用相移偏移方法得到的幅度譜和相位譜，而圖2（i）和（j）是對該含強噪聲信號用相移偏移與LSMP綜合方法得到的幅度譜和相位譜。由于目標信號幾乎被強噪聲淹沒，在相移偏移圖上噪聲依然很強，而目標信號非常微弱，幾乎不能識別出來；而在相移偏移與LSMP綜合成像圖上，噪聲得到壓制，目標信號反映清晰?？梢?，相移偏移與LSMP綜合方法明顯優(yōu)于相移偏移方法，不僅能夠?qū)ν噍S正確歸位，更具有很強的壓制噪聲的能力，從而極大地提高了低信噪比資料的成像質(zhì)量，較好地刻畫目標體的位置和形狀。另外，除了振幅譜外，相位譜也能較好地反映目標信息。由于相位譜反映信號相位變化的連續(xù)性，與反射信號能量無關(guān)，對于同相軸的連續(xù)展布以及突變，相位譜的特征更明顯。特別是對于含噪聲資料，相位譜比振幅譜能更明顯地突出目標信號，因此，將幅度譜和相位譜進行綜合分析，可以更好地識別和圈定目標。

3.2 實測數(shù)據(jù)實驗

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

4 結(jié) 論

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

［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偏移及最小熵技術(shù)的探地雷達成像法［J］.電子與信息學(xué)報，2007，29（4）：827－830.XIU Zhijie，CHEN Jie，F(xiàn)ANG 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］.電波科學(xué)學(xué)報，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］朱 明，普運偉，金煒東，等.基于時頻原子方法的雷達輻射源信號特征提?。跩］.電波科學(xué)學(xué)報，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］孫閩紅，唐 斌.基于原子分解理論的雷達欺騙式干擾信號特征提取.電波科學(xué)學(xué)報，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.