來五星 唐文靜 孫少山 鐘升

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

關(guān)鍵詞： 探地雷達(dá)； 成像方法； 后向投影； 子孔徑劃分； 合成孔徑； 成像效率

中圖分類號： TN958?34 文獻(xiàn)標(biāo)識碼： 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 引 言

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

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

1 標(biāo)準(zhǔn)BP成像算法

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

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

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

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

[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）為關(guān)于折射點[xr]的一元四次方程，直接求解過程繁瑣。根據(jù)Mast[10]提出的方法，折射點[xr]的近似求解公式如下：

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

在折射點位置[xr]求出解后，可以求得各道回波上對應(yīng)于待成像點[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]道回波上的散射響應(yīng)幅值：

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

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

[EP=xP， i] （8）

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

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

2 改進(jìn)的BP成像算法

假定待成像點[P]的坐標(biāo)為[x，z]，由標(biāo)準(zhǔn)BP成像算法的過程可知，[P]點成像結(jié)果可表述為：

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

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

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

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

3 算法性能對比

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

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

式中：[Etotal]為目標(biāo)總能量；[Emain]為主瓣能量。

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

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

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

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

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

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

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

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

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

4 結(jié) 語

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

Fig. 3 Imaging results of actual measured data

參考文獻(xiàn)

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

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

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

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

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

[4] 蔚建斌，陳自力，江濤.基于偏移技術(shù)的探地雷達(dá)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] 鄭文軍.復(fù)雜多散射環(huán)境下TRM成像技術(shù)研究[D].成都：電子科技大學(xué)，2010.

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

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

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

[9] 孔令講.淺地層探地雷達(dá)信號處理算法的研究[D].成都：電子科技大學(xué)，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的探地雷達(dá)圖像正演模擬[J].工程地球物理學(xué)報，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á)波動方程偏移[D].南京：河海大學(xué)，2004.

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