來五星 唐文靜 孫少山 鐘升

摘 要： 對脈沖探地雷達數(shù)據(jù)進行成像處理有利于地下目標的定位和探測。詳細闡述了探地雷達后向投影成像方法的實現(xiàn)原理與過程，并針對該方法存在計算效率低的問題，提出基于子孔徑劃分的改進后向投影成像方法，最后用頻率?波數(shù)成像、有限差分成像、標準后向投影成像以及改進的后向投影成像方法分別對仿真和實測數(shù)據(jù)進行成像處理并比較各成像方法的性能，驗證了改進算法的有效性。

關鍵詞： 探地雷達； 成像方法； 后向投影； 子孔徑劃分； 合成孔徑； 成像效率

中圖分類號： TN958?34 文獻標識碼： A 文章編號： 1004?373X（2018）15?0061?04

Research on back?projection imaging method of impulse ground penetrating radar

LAI Wuxing， TANG Wenjing， SUN Shaoshan， ZHONG Sheng

（School of Mechanical Science and Engineering， Huazhong University of Science and Technology， Wuhan 430074， China）

Abstract： The imaging processing of impulse ground penetrating radar （GPR） data is helpful for the location and detection of underground targets. The implementation principle and procedure of the standard back?projection imaging method of GPR are described in detail， and an improved back?projection imaging algorithm based on sub?aperture division is proposed to improve the computational efficiency of the standard imaging method. The frequency and wave number imaging， finite difference imaging， standard back?projection imaging and improved back?projection imaging methods are used to process the simulation data and actual measured data， their performances are compared， and the effectiveness of the improved algorithm is verified.

Keywords： ground penetrating radar； imaging method； back?projection； sub?aperture division； synthetic aperture； imaging efficiency

0 引 言

探地雷達（GPR）也稱地質(zhì)雷達，是對埋地目標進行定位與探測的有效工具，它集探測速度快、分辨率高、穿透力強、無損檢測等優(yōu)點于一身，被廣泛應用于市政管線檢測、道路與橋梁災害檢查、考古發(fā)掘探測、地雷與未爆彈藥的排查、礦物資源勘探等諸多方面[1?2]。脈沖探地雷達的原理是發(fā)射天線向地下輻射脈寬為納秒級的高頻電磁波，通過對接收到的回波信號的幅度、頻率以及相位等特征進行分析處理，得到被測物體的位置和大小信息。探地雷達成像就是對接收到的回波信號進行合成孔徑處理，從而準確地反映地下目標的位置。由于探地雷達回波數(shù)據(jù)與地震波數(shù)據(jù)具有相似性，地震數(shù)據(jù)處理使用的偏移技術可以應用在探地雷達的合成孔徑成像處理中，偏移可以將從地面接收到的來自地下的繞射波或反射波歸位到實際位置[3?4]，常用的偏移成像技術包括頻率?波數(shù)（[ω?k]）成像和后向投影（BP）成像。

[ω?k]成像基于波動方程以及“爆炸反射模型”的波場反向外推 [5?6]。該方法要求地下介質(zhì)均勻，即電磁波在地下介質(zhì)中的傳播速度要保持恒定，因此適用于淺地層目標探測或非層狀的均勻介質(zhì)。[ω?k]算法的優(yōu)勢是采用FFT，成像效率高。因此，[ω?k]算法在實際工程中廣泛使用。BP成像算法對地下分層介質(zhì)的適應能力強、成像精度高，在探地雷達成像領域的應用也較多，但該方法計算量大、實時性差。本文基于子孔徑的改進BP成像算法在保證成像質(zhì)量的前提下，可大幅減少BP算法的計算量，提高成像效率。

1 標準BP成像算法

后向投影（Back?Projection，BP）起源于醫(yī)學領域的CT技術，是McCorkle根據(jù)CT成像的投影切片理論推導出來的一種合成孔徑雷達時域成像方法[7]。其原理是對成像區(qū)域內(nèi)的任意一個待成像點先計算出該點距離各天線孔徑位置的時延，然后在各道接收回波上搜索對應該時延點的幅值，最后將各道回波中相應位置處的幅值進行時域相干疊加，從而完成成像過程[8?9]。因此，BP成像算法的核心思想就是“時延?求和”。

建立如圖1所示的探地雷達BP成像算法的場景模型。方位向沿[x]軸向右，距離向沿[z]軸向下，[k]個探地雷達收發(fā)天線與[x]軸平行，且與地面相距[h]，不同時刻的合成孔徑天線位置用倒三角來表示，其中圖示深色倒三角形表示當前發(fā)射和接收的第[i]個天線位置，其坐標記為[xi，-h]，淺色倒三角形為其余天線位置。成像場景被直線[z=0]分為兩部分，直線以上為空氣，直線以下假定為各向同性的均勻地下介質(zhì)。由于空氣與地下介質(zhì)的介電常數(shù)不相等，電磁波在空氣與地面分界處會發(fā)生折射，因此電磁波從空氣到地下介質(zhì)的實際傳播路徑為一條經(jīng)過折射點[（xr，0）]的折線，而不再是經(jīng)過點[（xl，0）]的直線。

令[θi]和[θr]分別為入射角和折射角，介質(zhì)相對介電常數(shù)為[εr]，成像場景中待成像點[P]的坐標為[xP，zP]，由Snell折射定律，有如下關系式：

根據(jù)三角幾何關系，有：

[sin θi=xr-xi（xr-xi）2+h2] （2）

[sin θr=xP-xr（xP-xi）2+z2P] （3）

聯(lián)立以上三式，可得：

[（xr-xi）2（xr-xi）2+h2?（xP-xr）2+z2P（xP-xr）2=εr] （4）

式（4）為關于折射點[xr]的一元四次方程，直接求解過程繁瑣。根據(jù)Mast[10]提出的方法，折射點[xr]的近似求解公式如下：

[xr≈x0+（xl-x0）εr] （5）

在折射點位置[xr]求出解后，可以求得各道回波上對應于待成像點[P]的時延為：

[τP，i=2（xr-xi）2+h2c+2（xP-xr）2+z2Pv] （6）

式中：[τP，i ]表示第[i]個合成孔徑位置到點[P]的雙程時延，[c]為電磁波在自由空間中的波速；[v=cεr]表示電磁波在介質(zhì)中的波速。

若第[i]道回波數(shù)據(jù)可表示為[Si（t）]，便可得到[P]點在第[i]道回波上的散射響應幅值：

[xP， i=Si（t=τP， i）] （7）

將[P]點在各道回波數(shù)據(jù)上的散射響應幅值進行疊加，便完成對[P]點的成像：

[EP=xP， i] （8）

式中[EP]為待成像點[P]的最終成像結果。

通過遍歷成像區(qū)域中全部待成像點，不斷重復上述“時延?求和”操作，即可實現(xiàn)對整個探地雷達圖像的BP成像。

2 改進的BP成像算法

假定待成像點[P]的坐標為[x，z]，由標準BP成像算法的過程可知，[P]點成像結果可表述為：

式中：[K]為總的天線數(shù)目；[Si（t）]為第[i]個天線孔徑位置處的回波數(shù)據(jù)；[τx，z，i]為[P]點坐標[x，z]到第[i]個天線位置的雙程時延。

由前述分析可知，標準BP成像算法存在一定的缺陷：一方面，在整個BP成像過程中，對于每個待成像點都需要遍歷全部[K]個天線孔徑位置，對大小為[M×N]的回波矩陣，算法的時間復雜度為[Θ（MNK）]，當探測區(qū)域較大時，算法運算速度慢、效率低；另一方面，對于成像區(qū)域內(nèi)的某個待成像點，并不是在所有天線位置都會產(chǎn)生散射響應，這些冗余計算進一步增大了BP算法的運算量，因此散射響應小或者沒有散射響應的位置在成像時可以不遍歷。準確地判斷埋地目標回波在接收的回波矩陣中的分布情況后，BP成像只利用目標回波矩陣進行運算，從而實現(xiàn)快速成像。基于此，改進后的BP成像算法將[K]個合成孔徑天線按照一定的準則劃分為[Nsub]個子孔徑，這樣式（9）中的求和操作無需遍歷全部[K]個合成孔徑位置，而是在天線孔徑數(shù)目為[ Ksub]的天線子孔徑內(nèi)進行相干疊加，得到若干個子孔徑圖像，最后合并全部子孔徑圖像即完成整個成像過程。因此，改進的BP算法復雜度為[ΘMNKsub]，運算時間為標準BP算法的[KsubK]。

根據(jù)改進的BP成像算法，第[j]個子孔徑對應的子孔徑圖像成像結果可表述為：

式中：[Ejx，z]為第[j]個子孔徑圖像中坐標點為[x，z]的最終成像結果；[Sj，it]為第[j]個子孔徑內(nèi)第[i]個天線孔徑位置的回波數(shù)據(jù)；[Ksub]為每個子孔徑內(nèi)所需遍歷的天線孔徑數(shù)；[Nsub=KKsub]表示子孔徑數(shù)。

3 算法性能對比

探地雷達成像算法的性能可以用方位向分辨率、積分旁瓣比（Integrated Side Lobe Ratio，ISLR））、運算時間等指標來衡量[11]。其中，方位向分辨率用歸一化掃描能量圖來表征，而ISLR是雷達信號處理中常用的一個指標，用來衡量算法對旁瓣和雜波的抑制程度，ISLR值越大則抑制效果越好，成像精度越高[12]。ISLR定義為目標所有旁瓣能量與主瓣能量之比，其表達式如下：

[ISLR=-10lgEtotal-EmainEmain] （11）

式中：[Etotal]為目標總能量；[Emain]為主瓣能量。

下面分別用仿真和實測數(shù)據(jù)進行實驗，并定量地用上述性能指標加以衡量和對比。

3.1 仿真數(shù)據(jù)實驗

GprMax2D為一款基于FDTD算法的探地雷達二維正演仿真軟件，常用于探地雷達成像方法的研究[13]。圖2為利用GprMax2D軟件模擬“鋼筋?混泥土”的場景并用不同方法進行成像。

仿真數(shù)據(jù)成像中各性能參數(shù)的大小見表1。

從表1可以看出：方位向分辨率上，標準BP算法和改進的BP算法最高；由ISLR可知，改進的BP算法ISLR值最大，表明其旁瓣和雜波抑制效果最好，成像精度最高，而[ω?k]成像的ISLR值最小，這在圖像上表現(xiàn)為寄生能量較BP成像多；從運行時間上看，[ω?k ]算法的時間開銷最少，標準BP成像算法的時間最長、效率最低；相對于標準BP成像算法而言，改進后的BP成像算法的運行時間僅為其1/5，不僅在時間消耗上大幅減少，而且在成像質(zhì)量上亦有所提升，證明了改進算法的有效性。

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

實測數(shù)據(jù)成像如圖3所示，圖3a）為國外某反雷技術研究中心實測的探地雷達原始B?Scan數(shù)據(jù)，圖像大小為512×98，反映的是對埋藏在地下的PMN?2地雷進行探測的結果。

實測數(shù)據(jù)成像中各算法的性能參數(shù)大小見表2。

從表2可以看出：實測數(shù)據(jù)實驗的處理結果與仿真實驗的結論不謀而合，即三種成像算法中，[ω?k]偏移成像算法耗時最短，實時性最好；標準BP成像算法耗時最長；改進的BP算法成像精度最高，耗時較標準BP算法明顯減少，但仍要大于[ω?k]算法。因此，當系統(tǒng)對成像算法的實時性要求較高時，應首選[ω?k]成像算法；而當注重成像質(zhì)量時，可采用改進的BP成像算法。

4 結 語

[ω?k ]成像基于波動方程和“爆炸反射模型”，成像質(zhì)量弱于BP算法，但實時性好。標準BP算法實現(xiàn)原理簡單，成像精度高，但運算效率低，實時性較差。仿真和實測數(shù)據(jù)實驗結果表明，本文改進的BP算法相對于標準BP算法，不僅成像質(zhì)量有所提升，而且成像重建的運算速度加快，成像效率得到明顯提升。

Fig. 3 Imaging results of actual measured data

參考文獻

[1] TREES H L V. Ground penetrating radar theory and applications [M]. US： Elsevier Science， 2009.

[2] 張春城.淺地層探地雷達中的信號處理技術研究[D].成都：電子科技大學，2005.

ZHANG Chuncheng. Study on signal processing technology in ground penetrating radar of shallow layer [D]. Chengdu： UESTC， 2005.

[3] 張金花.車載探地雷達天線特性分析及其成像處理研究[D].成都：西南交通大學，2014.

ZHANG Jinhua. Analysis of antenna characteristic of vehicle ground penetrating radar and study on imaging processing [D]. Chengdu： SWJTY， 2014.

[4] 蔚建斌，陳自力，江濤.基于偏移技術的探地雷達SAR成像方法[J].信號處理，2010，26（5）：778?782.

WEI J B， CHEN Z L， JIANG T. The SAR imaging method of GPR based on migration [J]. Signal processing， 2010， 26（5）： 778?782.

[5] CHEN G Y， BUI T D. Multiwavelets denoising using neighbo?ring coefficients [J]. IEEE signal processing， 2003， 10（7）： 211?214.

[6] ENGIN E， CIFTCIOGLU B， ?ZCAN M， et al. High resolution ultrawideband wall penetrating radar [J]. Microwave & optical technology letters， 2007， 49（2）： 320?325.

[7] 鄭文軍.復雜多散射環(huán)境下TRM成像技術研究[D].成都：電子科技大學，2010.

ZHENG W J. Study on TRM imaging technology under complex multi?scattering environment [D]. Chengdu： UESTC， 2010.

[8] 雷文太.脈沖GPR高分辨率成像算法研究[D].長沙：國防科技大學，2006.

LEI W T. Study on pulse GPR high?resolution imaging algorithm [D]. Changsha： NUDT， 2006.

[9] 孔令講.淺地層探地雷達信號處理算法的研究[D].成都：電子科技大學，2003.

KONG L J. Study on signal processing algorithm of ground penetrating radar of shallow layer [D]. Chengdu： UESTC， 2003.

[10] MAST J E， JOHANSSONE M. Three?dimensional ground pe?netrating radar imaging using synthetic aperture time?domain focusing [J]. Proceedings of SPIE： the international society for optical engineering， 1994， 2275： 205?214.

[11] OZDEMIR C， DEMIRCI S， YIGIT E， et al. A review on migration methods in B?scan ground penetrating radar imaging [J]. Mathematical problems in engineering， 2014 （1）： 1?16.

[12] YIGIT E， DEMIRCI S， OZDEMIR C， et al. Short?range ground?based synthetic aperture radar imaging： performance comparison between frequency?wavenumber migration and back?projection algorithms [J]. Journal of applied remote sen?sing， 2013， 7（1）： 382?385.

[13] 周奇才，李炳杰，鄭宇軒，等.基于GPRMax2D的探地雷達圖像正演模擬[J].工程地球物理學報，2008（4）：396?399.

ZHOU Q C， LI B J， ZHENG Y X， et al. Forward simulation of GPR image based on GPRMax2D [J]. Chinese journal of engineering geophysics， 2008（4）： 396?399.

[14] 李廣場.有限差分法探地雷達波動方程偏移[D].南京：河海大學，2004.

LI Guangchang. Wave equation offset of ground penetrating radar with finite difference method [D]. Nanjing： HHU， 2004.