孫浩然 蘭 天 楊小鵬

（1.北京理工大學(xué)信息與電子學(xué)院雷達(dá)技術(shù)研究所，北京 100081；2.北京理工大學(xué)重慶創(chuàng)新中心，重慶 401120）

1 引言

探地雷達(dá)（Ground Penetrating Radar，GPR）是一種用于探測地下場景的雷達(dá)系統(tǒng)，其利用電磁波在介質(zhì)電磁特性變化處的反射和散射現(xiàn)象來獲取地下場景的信息［1］。GPR系統(tǒng)相比于傳統(tǒng)對地探測方式，具有無損性、分辨率高、易于操作等優(yōu)勢，因此在武器探測、冰層凍土探測、管線定位以及考古探測等工程及環(huán)境領(lǐng)域得到了廣泛應(yīng)用［2-5］。

近年來，隨著探地雷達(dá)應(yīng)用場景的不斷增多，傳統(tǒng)GPR 系統(tǒng)收發(fā)一體、多點(diǎn)采集的數(shù)據(jù)獲取方式越來越不能滿足應(yīng)用中實(shí)時數(shù)據(jù)采集和處理的需求，因此多輸入多輸出（Multiple Input Multiple Output，MIMO）雷達(dá)體制被引入GPR 系統(tǒng)中，其利用空間中分布的多個發(fā)射正交波形的收發(fā)天線，從多角度照射目標(biāo)［6］，既明顯提高了數(shù)據(jù)采集速度，也豐富了接收的目標(biāo)散射數(shù)據(jù)，從而提升了最終目標(biāo)的分辨率［7］。但隨之而來的是MIMO-GPR下回波數(shù)據(jù)處理算法的改進(jìn)工作，其中偏移成像技術(shù)作為GPR 數(shù)據(jù)解釋的基礎(chǔ)步驟，其處理質(zhì)量將直接影響后續(xù)對地下場景中目標(biāo)的探測及分類等的準(zhǔn)確性。然而近年來對MIMO-GPR的研究多是著眼于MIMO陣列拓?fù)湓O(shè)計［8］或利用MIMO 雷達(dá)體制進(jìn)行高精度地質(zhì)探測［9］及目標(biāo)定位［10］等，而針對MIMO-GPR 系統(tǒng)的成像算法研究相對較少，因此探索基于MIMO-GPR的成像技術(shù)很有必要。

傳統(tǒng)偏移成像算法的研究可追溯到20 世紀(jì)60年代末，Schneider等人首次引入了衍射疊加偏移技術(shù)來重建成像區(qū)域的反射率圖［11］，即通過反向傳播每個采集到的時域距離剖面來將觀測數(shù)據(jù)從采集位置移回其真實(shí)的空間位置，這種方法需要對所有接收信號進(jìn)行后向傳播并積分，從而完成偏移區(qū)域的衍射疊加計算。由于這種方法實(shí)現(xiàn)簡單，靈活性強(qiáng)，因此被廣泛用于不同場景下的偏移成像處理。然而這種算法也存在一些缺點(diǎn)，例如當(dāng)回波中只有一個尖峰時，這種算法也會使得其出現(xiàn)在所有可能的反射點(diǎn)的位置，從而帶來較大的旁瓣［12］，這主要是因?yàn)樵摲椒ú⒉皇腔趪?yán)格的波動方程理論。為改善衍射疊加方法的缺陷，提高成像質(zhì)量，Schneider 提出了嚴(yán)格基于波動方程的基爾霍夫偏移方法，而后X.Liu 等人對其做了進(jìn)一步完善并應(yīng)用于諸多成像問題［13-14］，這種方法代表了惠更斯原理的定量描述，根據(jù)Kirchhoff 積分定理，利用接收波場值及其導(dǎo)數(shù)表示任意點(diǎn)的場值，表現(xiàn)出了比衍射疊加方法更高的成像質(zhì)量和更低的旁瓣偽影。

由于最初提出該傳統(tǒng)Kirchhoff 成像方法時主要用于地震數(shù)據(jù)處理，所以其充分利用了地震數(shù)據(jù)采集中逐點(diǎn)收發(fā)的特點(diǎn)基于爆炸源模型得到了該算法，因而其也僅適用于收發(fā)一體的雷達(dá)體制，而MIMO 探地雷達(dá)的收發(fā)天線分布于空間中不同位置，這也就直接導(dǎo)致傳統(tǒng)成像算法無法直接應(yīng)用于處理MIMO-GPR 采集的收發(fā)天線位置不同時的回波數(shù)據(jù)。為解決這一問題，Zhuge X 等人在傳統(tǒng)SISO 的成像方法的基礎(chǔ)上根據(jù)MIMO 場景對算法進(jìn)行了直接拓展，提出了適用于MIMO 情況下的Kirchhoff 成像方法［15］，但未給出嚴(yán)格證明，也沒有考慮遠(yuǎn)點(diǎn)衰減以及介質(zhì)分層的情況，因此并不能真正應(yīng)用于實(shí)際探地雷達(dá)成像中；而后Zhongmin Wang 等人在基于有效相位中心原理和標(biāo)量衍射理論提出了一種用于近場MIMO 成像的距離偏移方法［16］，實(shí)現(xiàn)了頻域MIMO成像，但此方法也僅適用于陣元間隔相同的情況。為了彌補(bǔ)之前工作中的不足，本文提出了一種基于MIMO-GPR 的時域偏移成像算法，并給出了頻域的嚴(yán)格證明，同時引入遠(yuǎn)點(diǎn)補(bǔ)償和折射點(diǎn)近似計算方法，實(shí)現(xiàn)了對地下目標(biāo)的精確成像。最后，本文對所提算法進(jìn)行了仿真驗(yàn)證，仿真結(jié)果顯示，利用本文方法所得的成像結(jié)果相較于利用傳統(tǒng)成像方法的結(jié)果有更低的旁瓣偽影和更高的圖像分辨率。

2 MIMO-GPR成像算法推導(dǎo)

2.1 衍射疊加原理

衍射疊加原理是時域偏移成像方法的基礎(chǔ)，在連續(xù)均勻介質(zhì)中，未經(jīng)過聚焦處理的點(diǎn)狀反射物的回波數(shù)據(jù)所表現(xiàn)出的形狀稱為衍射雙曲線，因此在偏移成像時，對于回波中的每一個波峰數(shù)據(jù)，都要在成像平面中構(gòu)造一條曲線，然后取每條曲線在該點(diǎn)的強(qiáng)度值并疊加，從而得到該點(diǎn)的像素值。如果所用GPR 雷達(dá)為收發(fā)一體式，且接收到的回波數(shù)據(jù)為u(k，t)，則由爆炸源模型，根據(jù)逆時偏移算法原理［11］，可將待測場景中某一點(diǎn)(x，z)的像素值表示為：

其中，v代表介質(zhì)中電磁波的傳播速度，f(x，z)表示(x，z)位置處的像素值，u(k，t)表示第k道數(shù)據(jù)中t時刻的回波場值，R1和R2分別代表電磁波從發(fā)射天線傳播到成像點(diǎn)和從成像點(diǎn)傳播回接收天線的距離，w(x，z)表示加權(quán)系數(shù)，主要用于補(bǔ)償收發(fā)機(jī)的空間分布不均勻所帶來的成像誤差。

2.2 Kirchhoff算法的時域推廣

從傳統(tǒng)的SISO-GPR 體制到MIMO 體制，所需解決的問題主要是如何充分有效利用多個接收機(jī)所采集的全部數(shù)據(jù)。對于Kirchhoff 這類衍射疊加算法而言，可以從算法原理出發(fā)進(jìn)行考慮，衍射疊加算法的本質(zhì)就是對接收信號進(jìn)行反向傳播并積分，而從這個角度，MIMO 體制下的Kirchhoff 成像可以認(rèn)為是對每一個收發(fā)對的回波信號進(jìn)行反向傳播并疊加，故對MIMO-Kirchhoff 算法成像表達(dá)式的推導(dǎo)可以通過對常規(guī)Kirchhoff 算法進(jìn)行一般化推廣而得到。因此，本文先對常規(guī)Kirchhoff 算法進(jìn)行推導(dǎo)，再根據(jù)MIMO 場景的不同特點(diǎn)對成像公式進(jìn)行一般化處理。

在均勻各向同性介質(zhì)中，矢量波動方程可以簡化為標(biāo)量波動方程，利用其Kirchhoff 積分計算式并進(jìn)行時間反向后，可以將波場表示為：

此時，可以選定Q面為半球面和水平地面形成的閉合曲面。由于在半球面上積分項對場值的貢獻(xiàn)為0，可將上式化簡，而后令=-cosθ，θ表示陣列與成像點(diǎn)的連線和豎直方向的夾角，可得常規(guī)Kirchhoff成像表達(dá)式為：

其中，S表示傳統(tǒng)SISO-GPR 系統(tǒng)數(shù)據(jù)采集平面上的所有采集點(diǎn)，v表示電磁波在介質(zhì)中的傳播速度，r表示收發(fā)陣元到成像點(diǎn)(x，z) 的距離，若用(x0，z0) 表示收發(fā)陣元的位置，則距離r=，τ表示接收信號的反向傳播時延，即τ=。

但為使得Kirchhoff 偏移成像算法可以滿足MIMO 陣列配置的要求，需要對上式進(jìn)行修正。首先由于此時雷達(dá)系統(tǒng)的天線不再收發(fā)一體，所以電磁波發(fā)射到目標(biāo)的路徑和從目標(biāo)反射到接收機(jī)的路徑不再相同，因此需要用θ1和r1來分別表示成像點(diǎn)和發(fā)射天線之間的夾角和距離，用θ2和r2來分別表示成像點(diǎn)和接收天線之間的夾角和距離，具體幾何關(guān)系見下圖1。

此時上述公式可以修正為

其中S0和S1分別表示雷達(dá)系統(tǒng)中的發(fā)射孔徑和接收孔徑。

2.3 MIMO-Kirchhoff算法的頻域證明

上述表達(dá)式是在常規(guī)Kirchhoff 算法的基礎(chǔ)上針對MIMO 雷達(dá)的不同之處而做出的普適性外推，本節(jié)中，我們將在頻率-波數(shù)域中對上述成像表達(dá)式展開推導(dǎo)證明。

首先假設(shè)所有發(fā)射機(jī)的發(fā)射信號強(qiáng)度都相同，并通過傅里葉變換將發(fā)射信號和目標(biāo)散射到接收天線的場變換到波數(shù)-頻率域，而后根據(jù)波場延拓關(guān)系［9］，可以將發(fā)射波場從發(fā)射天線延拓至散射點(diǎn)處，則可得到散射前的波場為：

其中Δz′=z1-z′0表示波場延拓距離，表示發(fā)射天線輻射出的源波場，H(kx，Δz，w)函數(shù)表示電磁波場在介質(zhì)中的外推延拓運(yùn)算，w表示電磁波信號的角頻率。同樣的可以將接收波場逆向延拓至散射點(diǎn)處，即得散射后的波場為：

其中，波場延拓距離Δz=z1-z0，表示接收天線收到的波場。由于當(dāng)前散射場中介質(zhì)分界面相對明顯，可以進(jìn)一步假設(shè)Born 近似在本散射過程中仍然適用，則可以得到散射后的波場和散射前波場的關(guān)系式為：

式中，fr表示散射體的散射系數(shù)。

一般情況下，成像處理對分辨率要求較高且待成像區(qū)域普遍較淺，因此GPR 所用發(fā)射信號通常為高頻脈沖信號，有時還會采用超寬帶（UWB）信號［17］，所以此處可以將發(fā)射信號理想化為一沖激信號，則發(fā)射波場為一常數(shù)，同時根據(jù)上述H(kx，Δz，w)函數(shù)的特點(diǎn)可以得到：

而后將上式逆傅里葉變換回空間-時間域，可得：

根據(jù)參考文獻(xiàn)［13］可得式中，

而后將式（10）中的卷積運(yùn)算寫成二維積分形式，并對其進(jìn)一步離散、化簡，即可得到：

其中，r1和r2分別表示發(fā)射天線和接收天線到成像點(diǎn)的距離，τ和τ′分別表示發(fā)射天線和接收天線到成像點(diǎn)的時延。此時，反射系數(shù)在空間中的分布反映的就是最終GPR 逆時偏移處理所得的成像結(jié)果。另外，按照圖1 中所示的幾何關(guān)系，可得上式和式（4）所表示的MIMO-Kirchhoff 算法成像表達(dá)式是一致的。

2.4 遠(yuǎn)點(diǎn)補(bǔ)償處理

根據(jù)電磁波在介質(zhì)中的傳播規(guī)律可知，當(dāng)介質(zhì)為非線性時，電磁波在介質(zhì)中傳播時除了相位會以傳播常數(shù)隨距離變化外，其幅值也會以衰減常數(shù)隨距離指數(shù)衰減。因此，距離收發(fā)天線越遠(yuǎn)的目標(biāo)，其回波的波場幅值越低。而上述成像表達(dá)式（12）中，由于被積分項的系數(shù)亦會隨著成像點(diǎn)離收發(fā)天線的距離增大而減小，因此利用上式對雷達(dá)回波數(shù)據(jù)進(jìn)行成像處理，必定會使得本就微弱的深處目標(biāo)的場值變得更小，從而帶來較為嚴(yán)重的有限孔徑問題，故需要對上式進(jìn)行修正以補(bǔ)償遠(yuǎn)點(diǎn)的場值。此處，為了給遠(yuǎn)點(diǎn)目標(biāo)以更大的權(quán)值，可以將成像表達(dá)式中的系數(shù)修正為r1r2，即可得到最終的成像表達(dá)式：

3 空耦下的MIMO-GPR成像處理

3.1 折射點(diǎn)計算

對于實(shí)際應(yīng)用中的一般情況，GPR 天線通常都無法布設(shè)于待探測場景所在的介質(zhì)中，而往往是位于空氣中，且距待探測場景的界面有一定距離，此時的成像場景示意圖如下圖2。

從上圖可以看出，由于空氣和介質(zhì)分界面的作用，電磁波的傳輸路徑不再是一條直線，而是一條被折射點(diǎn)(xr，0)分為兩段的折線，此處為了幾何關(guān)系表示清晰，只繪制了發(fā)射天線發(fā)出的電磁波到達(dá)成像點(diǎn)的路徑，而根據(jù)光路可逆原理，電磁波從成像點(diǎn)到接收天線的路徑與此相似，故此處只需研究發(fā)射路徑即可。根據(jù)Snell 折射定律以及圖中所示的幾何關(guān)系，可得

其中，θi和θr分別表示入射角和折射角。結(jié)合上式（14）～（15），可得

可以發(fā)現(xiàn)，這是一個關(guān)于xr的一元四次方程，而由一元多次方程的一般理論可知此方程沒有一般通用的求解方法，在實(shí)際中一般采用牛頓法等數(shù)值解法進(jìn)行求解。但由于牛頓迭代法需要較多的迭代次數(shù)，計算量很大，會使得算法的運(yùn)算速度大大下降，所以Mast等人提出了一種折射點(diǎn)近似求解方法［18］，可將復(fù)雜的求解過程簡單化，其給出的計算公式為：

其中xl=，而后Zhou 對其進(jìn)行了進(jìn)一步優(yōu)化［19］，得到了當(dāng)成像點(diǎn)和收發(fā)天線方位向距離較大時仍適用的近似求解折射點(diǎn)的方法：

3.2 空耦下的MIMO-Kirchhoff算法

在空耦合的情況下，利用波場延拓關(guān)系，重新表示散射前和散射后的波場。由于此時的介質(zhì)不連續(xù)，且波場在介質(zhì)分界面處折射前后只有幅值的變化，而不發(fā)生相位的移動，因此可以將上述二式重新表示為：

其中H1(kx，Δz，w)表示電磁波場在空氣層中的外推延拓運(yùn)算。

而后將上述兩式代入式（8）中，并對其進(jìn)行化簡、逆傅里葉變換和離散后，可得空耦合情況下MIMO-Kirchhoff的成像表達(dá)式為：

其中ra和rb分別表示發(fā)射天線到折射點(diǎn)和折射點(diǎn)到散射點(diǎn)的距離，分別表示散射點(diǎn)到折射點(diǎn)和折射點(diǎn)到接收天線的距離，εr1和εr2分別表示上下兩層介質(zhì)的介電常數(shù)。

至此，本文成功利用波場延拓關(guān)系，將MIMOKirchhoff 算法拓展到空氣耦合的應(yīng)用場景中，并得到了空耦適用的成像表達(dá)式。

4 數(shù)值仿真

為了充分驗(yàn)證算法的性能，本文分別建立了地耦合多點(diǎn)目標(biāo)模型和空耦合W 形目標(biāo)模型，并利用基于FDTD 的仿真工具gprMax［20］對上述模型進(jìn)行仿真，從而獲得了雷達(dá)回波數(shù)據(jù)。在下文中，將對上述兩種模型及相應(yīng)的仿真結(jié)果分別進(jìn)行介紹和分析。

4.1 地耦合多點(diǎn)目標(biāo)模型

為了說明文中MIMO-Kirchhoff算法的有效性以及遠(yuǎn)點(diǎn)補(bǔ)償處理的必要性，建立了如下圖3 所示的地耦合多點(diǎn)目標(biāo)模型。

如下圖所示，本仿真場景中Z=0 處為介質(zhì)分界面，收發(fā)天線陣均為線陣，且緊貼分界面，其中天線均為赫茲偶極子天線，發(fā)射和接收陣元均從x=0.2 m 處開始布置，發(fā)射陣元間隔為0.16 m，接收陣元間隔為0.032 m；介質(zhì)的介電常數(shù)εr=6，電導(dǎo)率為0，其中有3 個半徑均為0.03 m 的小球，具體位置如圖3所示。

由于小球半徑相對較小，同時也為獲得較高的距離向分辨率，本模型所用的發(fā)射信號頻率為6 GHz，波形為Ricker 子波。利用gprMax 對模型進(jìn)行仿真，可以得到該模型的B-Scan 回波如下圖4（a）所示。

從圖4（a）中可以發(fā)現(xiàn)，回波數(shù)據(jù)中最強(qiáng)波為界面直耦波，利用空場景對消法將其去除后，得到的回波如圖4（b）所示。從圖中可以較明顯的看到三條雙曲線，這也就代表地下的三個點(diǎn)目標(biāo)。而后用與MIMO 陣列中接收陣元數(shù)相同道數(shù)的傳統(tǒng)SISO體制的Kirchhoff 算法［21］作為對比，分別繪制傳統(tǒng)算法的成像結(jié)果和本文中未經(jīng)遠(yuǎn)點(diǎn)補(bǔ)償?shù)腗IMOKirchhoff算法的成像結(jié)果如下圖5所示。

從圖5（a）中可以看到，傳統(tǒng)算法得到的像中包含很多旁瓣偽影，且由于采集到的目標(biāo)散射信息不夠豐富，距離維分辨率也略差一些，而由圖5（b）可以看出未經(jīng)補(bǔ)償?shù)腗IMO-Kirchhoff算法雖然可以實(shí)現(xiàn)對較近目標(biāo)的精確成像，但深層目標(biāo)在成像結(jié)果中基本完全不可見。

而后繪制補(bǔ)償后的MIMO-Kirchhoff算法成像結(jié)果，如下圖6所示，可以看出成像結(jié)果中不僅沒有明顯的旁瓣偽影，還實(shí)現(xiàn)了深層目標(biāo)的精確成像。但對于淺層目標(biāo)，補(bǔ)償后的成像質(zhì)量不如補(bǔ)償前好，這主要是因?yàn)檫h(yuǎn)點(diǎn)補(bǔ)償處理在增強(qiáng)深層目標(biāo)的同時也增強(qiáng)了近目標(biāo)附近的噪聲，略微降低了成像質(zhì)量。

為定量描述算法的成像質(zhì)量，本文引入積分旁瓣比（Integrates Side Lobe Ratio，ISLR）來衡量幾種成像算法對圖像中旁瓣能量的抑制程度。ISLR 的定義為目標(biāo)所有旁瓣能量與主瓣能量的比值：

其中Etotal和Emain分別表示成像結(jié)果的總能量和目標(biāo)主瓣能量。從上式可知，圖像的ISLR 值越低，說明算法對圖像中旁瓣和雜波的抑制能力越強(qiáng)，成像質(zhì)量越好。

而后本文分別測量上述三種算法所成圖像中目標(biāo)點(diǎn)的位置坐標(biāo)，并計算其各自在方位維和距離維的ISLR 值，最終將所得結(jié)果整理如下表1所示。

由上表1可知，根據(jù)成像結(jié)果中的目標(biāo)位置，傳統(tǒng)Kirchhoff 算法所成像的位置偏差最大，平均偏差距離為0.0064 m，而補(bǔ)償后的MIMO-Kirchhoff 算法成像目標(biāo)位置偏差最小，平均偏差距離為0.0037 m；同時，根據(jù)成像結(jié)果的ISLR 對比，傳統(tǒng)Kirchhoff 成像算法對圖像旁瓣和雜波的抑制能力最差，而補(bǔ)償前MIMO-Kirchhoff 算法對雜波的抑制能力最強(qiáng)，因此當(dāng)目標(biāo)處于較近距離時使用此種算法最佳，當(dāng)目標(biāo)處于較遠(yuǎn)距離時則需要對算法進(jìn)行補(bǔ)償，以顯示出深層目標(biāo)。

表1 三種算法的成像指標(biāo)對比Tab.1 Comparison of imaging indexes of three algorithms

4.2 空耦合W形目標(biāo)模型

為了說明空耦合情況下MIMO-Kirchhoff算法的有效性，本文建立了如下圖7（a）所示的空耦合W 形目標(biāo)模型。其中，天線陣列配置和介質(zhì)參數(shù)設(shè)置均與圖3 中相同，但天線陣列距介質(zhì)分界面的豎直距離為0.1 m，介質(zhì)中目標(biāo)的形狀為W 形，其各頂點(diǎn)的位置如圖7（a）所示，發(fā)射信號仍采用6 GHz 的Ricker波形。

利用gprMax 對上述目標(biāo)模型進(jìn)行建模仿真，并對回波數(shù)據(jù)進(jìn)行空場景對消處理，而后利用式（21）所示的成像表達(dá)式對處理后的數(shù)據(jù)進(jìn)行偏移成像，最終的成像結(jié)果如圖7（b），從中可以看出本文所提的空耦MIMO-Kirchhoff算法在地下場景相對復(fù)雜時也可以達(dá)到不錯的成像效果，充分說明了本算法的有效性。

5 結(jié)論

本文通過對傳統(tǒng)探地雷達(dá)偏移方法進(jìn)行推廣，提出了一種基于Kirchhoff的MIMO 探地雷達(dá)成像方法，解決了傳統(tǒng)探地雷達(dá)不能實(shí)時采集數(shù)據(jù)并成像的缺陷。通過遠(yuǎn)點(diǎn)補(bǔ)償修正成像表達(dá)式中的加權(quán)系數(shù)，實(shí)現(xiàn)了對較深層目標(biāo)的探測和成像，而后引入折射點(diǎn)近似計算方法并對成像表達(dá)式進(jìn)行修正，得到了空耦場景下的MIMO-Kirchhoff 算法，最后本文利用gprMax 對地耦合多點(diǎn)目標(biāo)模型和空耦合W 形目標(biāo)模型進(jìn)行建模和仿真，結(jié)果表明該方法在地耦和空耦場景下都能達(dá)到較好的成像效果，具有較高的分辨率，可以實(shí)現(xiàn)對地下場景或目標(biāo)的精確成像。

由于算法權(quán)值的作用，在增強(qiáng)遠(yuǎn)點(diǎn)的目標(biāo)的同時也增強(qiáng)了部分噪聲，使得成像結(jié)果的ISLR 降低，在后續(xù)工作中可以考慮利用Sigmoid 等非線性加權(quán)函數(shù)對成像表達(dá)式的權(quán)值進(jìn)行修正。