邵泉杰，歐陽繕

(桂林電子科技大學 信息與通信學院，廣西 桂林541004)

探地雷達(ground penetrating radar，簡稱GPR)是一種地下目標無損檢測技術(shù)，具有探測速度快、分辨率高、操作靈活方便、探測范圍廣及探測費用低等優(yōu)點，被廣泛應用于公路檢測［1］、橋梁的無損檢測［2］、地質(zhì)工程探測和考古探測等［3］。GPR的基本原理為:發(fā)射天線在地面以寬頻帶窄脈沖的形式向地下發(fā)射高頻電磁波，入射波在介電性差異的2種介質(zhì)分界面(地層界面或目的體)產(chǎn)生反射波，接收天線接收記錄反射波的波形、振幅及到達時刻，并以雷達圖象方式顯示探測結(jié)果;根據(jù)測量的雙程走時和波速計算出目標體深度，連續(xù)測量剖面各點的反射波，形成雷達圖像。對探地雷達回波數(shù)據(jù)進行處理，最終目的是解釋地下異常的情況，其中目標成像是最直觀有效的數(shù)據(jù)處理方式。已有的成像算法研究主要包括反向投影(back projection，簡稱BP)算法［4］、時間反轉(zhuǎn)(time reversal，簡稱TR)成像算法［5］、壓縮感知(compressive sensing，簡稱CS)算法［6-7］、反演(inversion)成像［8］等。其中BP算法最為經(jīng)典，該算法對成像域進行像素點劃分，通過計算回波時延進行能量累加成像，實現(xiàn)簡單，但運算量大。TR算法充分利用電磁波傳播的物理過程，適合復雜多散射環(huán)境，具有較高的成像分辨率。反演成像一般基于Maxwell方程，能計算出目標的位置、性質(zhì)和尺寸等，但其過程需迭代等運算，復雜度高。上述方法大都計算量大，且不適合實時處理。CS算法雖可降低數(shù)據(jù)存儲和運算復雜度，以較少數(shù)據(jù)對信號進行重構(gòu)，但對測量矩陣和重構(gòu)算法要求嚴格。

橢圓包絡(luò)線算法［9］是一種目標邊界成像算法，該算法基于收發(fā)分置天線提出，利用天線與目標邊界所構(gòu)成的幾何相切關(guān)系成像，實時性好，計算效率高。穿墻雷達只需考慮墻體的影響，對墻體進行補償后，電磁波相當于在單一的空氣介質(zhì)中傳播。但探地雷達對目標體進行探測過程中，由于地下介質(zhì)的復雜度高，地下情況未知，可能存在多層介質(zhì)，且多層介質(zhì)的形狀不像墻體那樣水平，補償?shù)姆椒黠@不適用。為此，以一種近似構(gòu)造橢圓的方法將包絡(luò)線算法擴展到探地雷達地下目標成像，使其適用于多層介質(zhì)的地下環(huán)境。近似后的橢圓包絡(luò)線算法對多層介質(zhì)地下目標成像深度有微小的誤差，通過分析誤差影響因素，可發(fā)現(xiàn)算法的最佳適用環(huán)境。

1 橢圓包絡(luò)線目標成像算法

1.1 穿墻雷達中的橢圓包絡(luò)線算法

包絡(luò)線算法用于穿墻雷達時，墻體是對成像產(chǎn)生誤差的主要外界因素。由于墻體比較均勻，厚度一定，可采用補償方法消除墻體影響。基于收發(fā)分置的橢圓包絡(luò)線算法用于穿墻雷達的幾何原理如圖1所示，圖中墻體已進行補償。收發(fā)天線沿著側(cè)線即x軸自左向右移動，目標物體位于空氣介質(zhì)中，假設(shè)在一次天線位置，發(fā)射天線坐標為F1(XT，0)，接收天線坐標為F2(XR，0)，電磁波照射目標體的C(x，y)點，天線與目標體之間的回波往返距離為Y=Y1+Y2。

圖1 穿墻雷達幾何原理圖Fig.1 Geometry relationship of through-wall radar

幾何學上橢圓的定義為:平面上到兩點距離之和為常數(shù)的點的軌跡。為此，在固定天線位置處，收發(fā)天線坐標視為2個定點，電磁波在空氣中傳播的距離為常數(shù)，便可得天線、回波路徑以及目標邊界之間的橢圓關(guān)系。將發(fā)送天線與接收天線坐標看成橢圓的2個焦點坐標F1(XT，0)和F2(XR，0)，則半焦距為。天線與目標點之間的回波路徑長度Y=c t=2a。其中:c為電磁波在空氣中的傳播速度;t為電磁波從發(fā)射天線到達目標邊界點后返回接收天線所用的時間。由此計算橢圓長軸a=v t/2，根據(jù)橢圓的性質(zhì)，可得短軸。如此，回波路徑與天線坐標可建立幾何關(guān)系，得到橢圓方程表達式:

收發(fā)天線均以某一固定間隔距離ΔX沿橫軸自左向右同時移動，發(fā)送天線坐標集合XT={X1，X2，…，XN}，對應的接收天線坐標集合XR={X1－ΔX，X2－ΔX，…，XN－ΔX}。每個天線位置對應一個橢圓，并與目標邊界的一個點相切。眾多橢圓包絡(luò)簇實現(xiàn)回波雙曲線的能量聚焦，基本無誤差地實現(xiàn)對目標上邊界的還原。

1.2 探地雷達中的橢圓包絡(luò)線算法

探地雷達針對地下目標體進行探測，假設(shè)目標體為圓形，收發(fā)分置天線位于地表，電磁波進入地下環(huán)境介質(zhì)中傳播，設(shè)地下介電常數(shù)為μ1，目標體的介電常數(shù)為μ2，目標體的上下界面均為曲面，電磁波穿透目標體，上下表面均會產(chǎn)生回波。上表面的回波為單一介質(zhì)情況，與墻體補償后的穿墻雷達成像相似，唯一區(qū)別為介質(zhì)中的波速。下表面的回波經(jīng)過2種介質(zhì)，即地下環(huán)境介質(zhì)及目標體介質(zhì)，而且2種介質(zhì)的分界面均為曲面，這種情況下對下表面進行成像，即多層介質(zhì)成像。構(gòu)造橢圓時需要考慮2種介質(zhì)的情況，其幾何關(guān)系如圖2所示。對于多層介質(zhì)情況，橢圓焦點坐標仍然為發(fā)射天線和接收天線的坐標。由于接收到的回波數(shù)據(jù)為電磁波的時延，長軸的計算采用不同介質(zhì)內(nèi)電磁波運行距離分段相加近似計算的方法，

其中:Y1+Y4為電磁波在介質(zhì)環(huán)境中往返路程;Y2+Y3為目標體內(nèi)往返路程;t1為目標上表面回波往返時間;t2為目標下表面的回波往返時間;v1為電磁波在環(huán)境介質(zhì)中的傳播速度;v2為電磁波在目標介質(zhì)中的傳播速度。橢圓長軸a=Y/2，通過b=計算短軸，構(gòu)造橢圓包絡(luò)線，從而實現(xiàn)下界面的邊界重建。對于實際中的天線離地較遠情況，相當于增加一層空氣介質(zhì)，分段相加的方法仍適用。

圖2 地下多層介質(zhì)幾何關(guān)系Fig.2 Geometric relationship of underground multi-medias

所有天線位置構(gòu)造的包絡(luò)線簇呈現(xiàn)目標體的邊界形狀。為了更直觀地觀察，提取反映目標區(qū)域的邊界，這需確定邊界的起始點和終止點以及2個端點之間的值。假設(shè)天線移動過程中總計產(chǎn)生M個橢圓，天線離目標體越遠，電磁波照射到的目標體的位置越邊緣，故以第1個和第2個橢圓的交點橫坐標為起始點xmin，第M－1個和第M個橢圓的交點為終止點xmax。目標邊界有凹和凸2種情況，與橢圓的對應關(guān)系分別為內(nèi)切和外切，圖2所示的目標體下界面為凹形，構(gòu)造的橢圓與下界面外切，按照文獻［9］確定形狀與橢圓的相切關(guān)系。在劃定成像區(qū)域范圍［xmin，xmax］且y≥0中，橢圓的最小個數(shù)為N，橢圓與坐標軸圍成的橢圓區(qū)域集合為?E，即?E=。?E1與?EN交點為P，目標邊界為?B。設(shè)?Eα為?E邊界點的集合，?Eβ為?E內(nèi)部點的集合，滿足并集關(guān)系?E=?Eα∪?Eβ。根據(jù)P所在的集合建立邊界提取的數(shù)學描述，其表達式為:

其中:?Eout和?Ein分別為橢圓與目標外切和內(nèi)切時的成像邊界集合;ai、bi、ci和xi分別為橢圓?Ei的半長軸、半短軸、半焦距和發(fā)射天線橫坐標。

1.3 算法性能及可行性分析

穿透多層介質(zhì)時，以一種分段相加近似計算的關(guān)系構(gòu)造橢圓長軸。由Snell折射定律知，電磁波在介電特性不同的2種介質(zhì)中，入射角和折射角的大小與2種介質(zhì)的介電特性有關(guān)。例如，圖2中以Y1+Y2+Y3+Y4近似計算長軸，而實際長軸計算應當為2點間直線距離F1C+CF2，由三角形關(guān)系知，三角形任意二邊之和大于第三邊，所以穿透多層介質(zhì)時計算的深度略微有所增加。

誤差產(chǎn)生的原因主要是多層介質(zhì)介電特性的差異，反映到距離上，即為電磁波在不同介質(zhì)中的運行距離。與電磁波運行距離最直接相關(guān)的因素是目標體埋深。驗證目標體埋深對誤差影響，設(shè)置不同深度的模型，如圖3所示。電磁波從介質(zhì)1入射，經(jīng)過交界面進入介質(zhì)2，介質(zhì)2內(nèi)存在目標體，假設(shè)介質(zhì)1和介質(zhì)2的介電常數(shù)不變，目標體深度不同。實際中的長軸為x1，近似構(gòu)造的長軸為a1+a2，由三角形性質(zhì)可得如下關(guān)系:

在介電常數(shù)一定情況下，角度θ是固定的，誤差Δ只與天線在介質(zhì)1的高度以及目標體的埋深有關(guān)。深度a1越小，天線離地距離a2越小，誤差越小。

圖3 不同深度模型Fig.3 Different depth model

2 算法仿真與實測數(shù)據(jù)分析

2.1 仿真結(jié)果與分析

采用GprMaxV2.0進行仿真，驗證橢圓包絡(luò)線算法對目標邊界成像的有效性，并探討探地雷達成像的性能以及深度對多層介質(zhì)成像的影響。仿真場景為:地下區(qū)域長2 m、寬1.2 m，環(huán)境的介電常數(shù)為3，地下埋一個圓柱體，圓心位于地下0.5 m，半徑分別設(shè)置為0.2、0.1 m，介電常數(shù)為10;天線采用收發(fā)分置間距0.1 m，每次移動0.03 m，總計移動60步，頻率900 MHz，2種介質(zhì)的介電特性已知，求得環(huán)境介質(zhì)和目標體內(nèi)電磁波傳播速度，對圓柱上下界面進行包絡(luò)線成像，仿真結(jié)果如圖4所示。

圖4 包絡(luò)線成像的仿真結(jié)果Fig.4 Simulation results of envelope imaging

由仿真結(jié)果可知，上表面只涉及單一介質(zhì)，在已知介質(zhì)波速情況下，不管深度還是形狀，其還原均十分準確。由于分段相加近似計算的長軸偏長，在半徑為0.2 m情況下，圓的下表面成像產(chǎn)生了接近0.02 m的深度誤差。但相同仿真條件下，半徑為0.1 m的圓下表面的誤差基本可忽略。上述結(jié)果驗證了深度越小，近似計算的包絡(luò)線算法誤差越小?？梢姡瑢τ跍\層目標成像，近似計算的包絡(luò)線算法具有較好的性能。

2.2 實測數(shù)據(jù)成像與分析

實測數(shù)據(jù)實驗場景描述:沙坑中埋有一個圓弧形的鐵片，長度為0.18 m，每次移動0.02 m，總計移動41步。使用的實驗設(shè)備為脈沖源、喇叭天線、數(shù)據(jù)采集示波器、低噪聲放大器等。脈沖源為二階高斯脈沖，6 db帶寬0.32～1.56 GHz;發(fā)射接收天線型號ATB1G-12-003，均為極化天線;低噪聲放大器型號ADL5523;示波器型號安捷倫ium 80000，采樣頻率40 GHz，數(shù)據(jù)經(jīng)示波器顯示并保存。

波速通過已知目標深度換算的方法［10］精確求取。提取回波延時，根據(jù)天線移動情況，應用多層介質(zhì)情況下的包絡(luò)線成像方法，成像結(jié)果如圖5所示。由于物體埋深很淺，原始數(shù)據(jù)采集較為準確，與物體上表面的實線進行對比，成像后誤差較小，上表面的形狀還原精確。

圖5 包絡(luò)線成像結(jié)果Fig.5 Envelope imaging result

與CBP成像算法［11］比較，原始數(shù)據(jù)經(jīng)CBP成像后的結(jié)果如圖6所示。從圖6可見，改進的BP算法只能粗略地對能量進行聚焦，雖然物體的大致位置信息等基本正確，但根據(jù)成像結(jié)果很難分辨出物體的形狀，圖中目標物體位置出現(xiàn)2個類似目標體的情況，這會導致地下情況的錯誤解讀。BP成像算法需要對眾多像素點進行能量疊加，導致計算復雜度較高。在奔騰E5400處理器2 GB內(nèi)存情況下，Matlab進行BP成像所需時間為49.54 s，包絡(luò)線成像所需時間為2.58 s?？梢?，包絡(luò)線算法的處理速度明顯優(yōu)于BP算法，更適合實時處理。

圖6 CBP成像結(jié)果Fig.6 CBP imaging result

3 結(jié)束語

將包絡(luò)線成像算法應用于探地雷達地下目標成像，根據(jù)多層介質(zhì)內(nèi)電磁波照射物體往返路徑分段相加的近似方法構(gòu)造橢圓的長軸，通過理論分析及實測數(shù)據(jù)驗證淺層目標成像，近似相加方法的誤差在可允許范圍之內(nèi)。與BP成像算法進行對比，包絡(luò)線算法對邊界的構(gòu)造準確，運算速度快，更適合實時處理。

［1］盧成明，秦臻，朱海龍，等.探地雷達檢測公路結(jié)構(gòu)層隱含裂縫實用方法研究［J］.地球物理學報，2007，50(5):1558-1568.

［2］Wang Z W，Zhou M，Slabaugh G G，et al.Automatic detection of bridge deck condition from ground penetrating radar images［J］.IEEE Transactions on Automation Science and Engineering，2011，8(3):633-640.

［3］曾昭發(fā)，劉四新，王者江.探地雷達方法原理及應用［M］.北京:科學出版社，2006:213-223.

［4］Zhou L，Huang C，Su Y.A fast back-projection algorithm based on cross correlation for GPR imaging［J］.Geoscience and Remote Sensing Letters，2012，9(2):228-232.

［5］鄭文軍，趙志欽，張薇，等.超寬帶探地雷達中TRM-SAR成像技術(shù)研究［J］.電子科技大學學報，2011，40(3):363-366.

［6］Tuncer M A C，Gurbuz A C.Ground reflection removal in compressive sensing ground penetrating radars［J］.Geoscience and Remote Sensing Letters，2012，9(1):23-27.

［7］Gurbuz A C，Mcclellan J H，JR Scott W R.Compressive sensing of underground structures using GPR［J］.Digital Signal Processing，2012，22(1):66-73.

［8］丁亮，韓波，劉潤澤，等.基于探地雷達的混凝土無損檢測反演成像方法［J］.地球物理學報，2012，55(1):317-326.

［9］李育暉，歐陽繕，晉良念，等.超寬帶穿墻雷達橢圓包絡(luò)線目標邊界成像算法［J］.電子與信息學報，2014，36(7):1532-1537.

［10］鄧小燕，王通.探地雷達探測中對媒質(zhì)相對介電常數(shù)的測定［J］.物探與化探，2009，33(1):43-45.

［11］周琳，粟毅.基于互相關(guān)的探地雷達反向投影成像算法［J］.電子與信息學報，2011，33(11):2714-2719.