來五星 唐文靜 孫少山 鐘升

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

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

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

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

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

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

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

令[θi]和[θr]分別為入射角和折射角，介質(zhì)相對(duì)介電常數(shù)為[εr]，成像場(chǎng)景中待成像點(diǎn)[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)于折射點(diǎn)[xr]的一元四次方程，直接求解過程繁瑣。根據(jù)Mast[10]提出的方法，折射點(diǎn)[xr]的近似求解公式如下：

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

在折射點(diǎn)位置[xr]求出解后，可以求得各道回波上對(duì)應(yīng)于待成像點(diǎn)[P]的時(shí)延為：

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

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

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

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

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

[EP=xP， i] （8）

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

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

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

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

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

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

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

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

3 算法性能對(duì)比

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

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

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

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

3.1 仿真數(shù)據(jù)實(shí)驗(yàn)

GprMax2D為一款基于FDTD算法的探地雷達(dá)二維正演仿真軟件，常用于探地雷達(dá)成像方法的研究[13]。圖2為利用GprMax2D軟件模擬“鋼筋?混泥土”的場(chǎng)景并用不同方法進(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成像多；從運(yùn)行時(shí)間上看，[ω?k ]算法的時(shí)間開銷最少，標(biāo)準(zhǔn)BP成像算法的時(shí)間最長(zhǎng)、效率最低；相對(duì)于標(biāo)準(zhǔn)BP成像算法而言，改進(jìn)后的BP成像算法的運(yùn)行時(shí)間僅為其1/5，不僅在時(shí)間消耗上大幅減少，而且在成像質(zhì)量上亦有所提升，證明了改進(jìn)算法的有效性。

3.2 實(shí)測(cè)數(shù)據(jù)實(shí)驗(yàn)

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

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

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

4 結(jié) 語

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

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á)中的信號(hào)處理技術(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].信號(hào)處理，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].長(zhǎng)沙：國防科技大學(xué)，2006.

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

[9] 孔令講.淺地層探地雷達(dá)信號(hào)處理算法的研究[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é)報(bào)，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] 李廣場(chǎng).有限差分法探地雷達(dá)波動(dòng)方程偏移[D].南京：河海大學(xué)，2004.

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