基于折射波靜校正的未知參數穿墻成像合成延遲時算法*

李家強**1a,，劉 然，徐才秀2，陳金立1a,，陳焱博，朱艷萍

(1.南京信息工程大學 a.氣象災害預報預警與評估協同創新中心;b.電子與信息工程學院，南京210044；2.江蘇北斗衛星導航檢測中心公司，南京210044)

0 引 言

穿墻成像雷達(Though-the-Wall Radar Imaging，TWRI)作為一種無損檢測技術，能夠探測視覺非透明障礙物后的目標，在反恐、救災等方面發揮了重要作用，廣泛應用于軍事、民用等領域[1]。

近年來，在穿墻成像雷達方面有許多研究致力于對未知參數的單層墻體后目標進行成像。其中文獻[2]討論了相對介電常數和墻壁厚度對目標點成像的影響，通過設置兩組與墻壁不同距離的天線陣列，對預先假設的幾組墻壁厚度和介電常數的場景進行成像，兩組陣列所成圖像都會隨著墻體厚度和相對介電常數的誤差而產生偏移，其中兩者交叉的點即為真實目標位置。由于該方法需要對設置的不同墻體參數進行成像，此過程計算復雜，且墻體參數的估計結果誤差較大。另外，支持向量機方法(Support Vector Machine,SVM)[3]根據散射信號提取出信號作為輸入、墻體參數作為輸出得到的函數關系，對墻后目標進行成像。此方法所估計出的墻體參數值易受目標位置、大小以及形狀的變化影響。文獻[4]運用兩列垂直天線陣列對穿墻場景進行成像，利用中心關聯目標的位置關系對穿墻場景進行重聚焦，但該方法易受目標的種類和位置的影響。對于多層墻體成像的研究，文獻[5]主要針對成像聚焦時延的計算問題，其每層墻體的參數都是假設已知的；文獻[6]通過構造濾波器去除回波域中墻體產生的影響以獲得清晰的成像結果，該方法每次更新值時都要進行重新計算并成像，故計算效率低，且成像效果易受實際墻體參數的影響;文獻[7]通過對墻體厚度和相對介電常數的二維搜索，得到延時補償值的取值范圍，進而對該取值范圍進行遍歷，找到最佳聚焦時延近似補償由于墻體存在產生的異常傳播時延。在對墻體參數進行搜索時，該方法不能精確地計算聚焦時延的取值范圍，導致在遍歷時迭代次數過多。

在實際應用中，墻體一般并非單層均勻介質，且墻體參數未知。由于電磁波在墻體內的傳播速度比在空氣中傳播速度慢，導致其回波時延增長，成像過程中目標點位置會產生偏移或出現散焦現象，這給成像工作帶來巨大挑戰。聚焦時延的計算是影響穿墻成像質量至關重要的因素，因此，解決未知參數多層墻體的時延計算問題是目前穿墻成像面臨的挑戰之一。為解決未知參數的多層墻體介質，如常見的雙層墻體成像問題，本文提出了基于折射波靜校正的合成延遲時法對墻后目標進行精準地聚焦成像，通過計算出每一個通道內的聚焦時延，進而依次對每個通道中由于墻體存在產生的異常傳播時延進行一一校正。

1 穿墻成像雷達場景模型

本文運用多發多收的天線陣列探測墻后目標。如圖1所示，墻體由兩層不同介質組成，內層墻體為石膏材料，外層墻體為混凝土材料，兩者厚度分別為d1和d2，相對介電常數分別為ε1和ε2。天線陣列緊貼墻體依次擺放，發射信號為超寬帶窄脈沖信號[7]，信號經過墻體的透射和反射至目標點處，再由接收天線所接收。

圖1 電磁波發射接收傳播時延

本文運用基于時域有限差分法(Finite Difference Time Domain,FDTD)的GprMax2D/3D[8]軟件對穿墻場景進行仿真并獲取回波數據。設發射信號為s(t)為負一階高斯脈沖信號，其表達式為

(1)

設第m(m=1,2,…,M)根發射天線坐標為(xtm,0)，第n(n=1,2,…,N)根接收天線坐標為(xrn,0)，目標點q坐標為(xq,yq)，在不考慮噪聲的影響下，則接收天線接收的回波信號為aqs(t-τmqn)，其中aq為目標的回波幅值，τmqn為電磁波從第m根發射天線到第n根接收天線的傳播時延。由圖1可知,

(2)

式中：c為電磁波在空氣中的傳播速度,τtm1與τtm2分別為電磁波在外層和內層墻體中發射天線延遲時，τrn1與τrn2分別為電磁波在兩層墻體中接收天線的延遲時,l1與l2為電磁波在空氣中的傳播路徑。由后向投影(Back Projection)算法知，成像場景區域某一像素點p(xp,yp)的像素值為

(3)

式中：τmpn為成像區域任意點p從第m根發射天線到第n根接收天線的傳播時延。

若兩層墻體參數已知，則可通過折射定理計算折射點坐標，進而計算出電磁波在墻體內的傳播時延，完成對墻后目標的成像與精準定位。而在實際應用中，墻體厚度以及相對介電常數這些先驗信息通常是未知的，這給后續的成像工作帶來影響，導致目標位置偏移和成像散焦。本文提出基于折射波靜校正的合成延遲時法，通過計算每一個通道的聚焦時延，使得目標更接近真實位置，圖像聚焦程度更高。

2 基于折射波靜校正的合成延遲時算法

2.1 折射波靜校正原理

靜校正技術屬于地震勘探領域，是一項處理地震資料的技術。假設存在兩層不同介質，電磁波在下層介質中的傳播速度大于在上層介質的傳播速度時，則透射波會沿水平界面滑動，滑行波在上層介質中激發新的波動即折射波。基于基本折射方程[9]，折射波靜校正通過讀取折射波初至時間[10]計算介質的厚度以及速度信息，從而得到地震波在介質中的傳播時間，完成對由于地表條件起伏引起的時差的校正。故應用折射波靜校正必須具備兩個基本條件：一是要存在水平的折射面；二是下層速度要大于上層速度[11]。在本文設置的穿墻成像場景模型中，由于墻體設定為平行，故其擁有相對穩定的折射面；另外，由于外層墻體的相對介電常數比內層墻體相對介電常數大，電磁波在介質內的傳播速度與其相對介電常數成反比，故電磁波在外層墻體內傳播速度比在內層傳播速度小。因此，折射波靜校正理論運用于穿墻成像領域能夠計算電磁波在墻體內的傳播時延，以校正每一通道中由于墻體存在引起的異常傳播時延誤差。

在本文設置的穿墻場景中，墻體由兩層不同介質組成。當電磁波從發射天線發出后，一部分透射波會沿著兩層墻體的界面滑動，另一部分透射波在內層墻體與空氣的界面上滑動，引起墻體振動，產生折射波。折射波的折射路徑如圖2所示，在本文模型中折射面有兩個：內層墻體與空氣接觸的界面設為折射面1，兩層相對介電常數不同的墻體接觸的界面設為折射面2。由基本折射方程知：

圖2 折射波延遲時路徑

(4)

(5)

2.2 外層墻體電磁波傳播延遲時

(6)

令

(7)

可得到一系列延遲時超定方程組，則式(5)可寫為

Gx=b。

(8)

式中：G=[g1,g2,…,gm]T為該超定方程組的系數矩陣，x為發射天線以及接收天線延遲時構成的未知數向量。對式(8)作差可得

(9)

對式(9)構建最小二乘：

(10)

對x求導得

(11)

式中：tr[xxTGTG]可寫為tr[xTGTGx]，整理可得

(12)

其中：

(13)

將式(13)代入式(12)中得

(14)

所以有

(15)

則

(16)

則有

x=(GTG)-1b。

(17)

求得的向量x為該超定方程組的解，即電磁波在外層墻體內的發射點傳播延遲時τm,2和接收點傳播延遲時τn,2。

2.3 內層墻體電磁波傳播時延

(18)

(19)

對于相鄰接收天線陣元n與n-1，其截距時間差[17]可表示為

(20)

則求得的相鄰發射天線截距時間差即相鄰發射天線延遲時間增量，而相鄰接收天線截距時間差即為相鄰接收天線延遲時間增量。

τ=τm,1+τn,1+τm,2+τn,2+τ″ 。

(21)

其中：

(22)

3 實驗與分析

采用GprMax2D/3D建立圖1所示的穿墻仿真場景。天線陣列緊貼外側墻面擺放，采用3發13收天線裝置，目標為球形理想導體，位于坐標(1.1,1.2)處。若不考慮墻體存在，使用后向投影算法對目標場景進行直接成像，其成像結果如圖3所示，可見由于墻體存在引起的異常傳播時延，成像目標位置會嚴重偏離真實坐標。

圖3 原始成像效果圖

圖4 接收天線1的回波信號波形

(a)回波時刻離散點

如上文所述，首先通過解式(8)超定方程組求得外層墻體的發射天線與接收天線延遲時，再計算每個通道內的發射天線與接收天線延遲時增量來確定內層墻體的延遲時，最后令其兩者相加則可得電磁波在兩層墻體內的傳播延遲時。圖6給出了本文算法的流程圖。

圖6 合成延遲時算法流程圖

為驗證本文算法在不同墻體兩層介質參數下的有效性，設置以下三組不同情況墻體參數值，每組墻體由介質1和介質2組成，相對介電常數和厚度分別為(ε1,d1)和(ε2,d2)，如表1所示。運用基于折射波靜校正的合成延遲時算法對三種不同墻體參數情況時穿墻場景的成像結果如圖7所示，同時對比了運用回波域濾波方法[6]以及自聚焦方法[7]在不同墻體參數值下的成像情況。與回波域濾波算法對比，本文所提算法對墻后目標位置的定位更準確，且更不易受墻體參數的影響。自聚焦算法成像遍歷電磁波的傳播時延的取值范圍，通過圖像評價標準選取成像結果最佳圖像確定電磁波傳播時延的值，故當墻體參數改變時，其成像目標位置不易受墻體參數影響，而本文算法與自聚焦算法對比，其成像目標位置更精準。

表1 不同墻體參數取值

圖7 本文算法與兩種算法成像結果對比

為進一步分析算法的性能，對本文算法成像結果的圖像熵進行對比。圖像熵能夠衡量墻后目標物體成像的聚焦程度，其值越小則聚焦程度越好。若像素點個數為X×Y，p(i,j)表示第i行第j列的像素值，則該圖像的圖像熵ξ可表示為[14]

(23)

根據式(23) 求得在不同墻體參數情況下本文所提算法與回波域濾波算法和自聚焦算法成像結果的圖像熵值，其對比如表2所示。

表2 三種算法成像結果圖像熵對比

由表2可知，在不同墻體參數環境下，本文算法與兩種算法相比，其圖像熵值較小，即本文算法成像結果在聚焦程度方面比上述兩種算法效果更優。

分別對回波域濾波算法、自聚焦算法以及本文所提算法進行計算復雜度對比分析。對于回波域濾波算法，其計算復雜度可表示為

C0=mnlgp+mN+3N。

(24)

式中：P為采樣點數，m、n分別為穿墻成像區域網格橫向和縱向的劃分個數。對于自聚焦算法，其計算復雜度可表示為

C1=LmnP+mN。

(25)

式中：L為自聚焦算法延時補償值取值范圍的遍歷次數。而本文所提基于折射波靜校正的合成延遲時算法，其計算復雜度如式(26)所示：

C2=mnlg 2P。

(26)

圖8給出了三種算法在不同采樣點數情況下計算復雜度的曲線圖。由圖8可知，在相同采樣點數條件下，本文所提算法相較于其他兩種算法，其計算復雜度相對較低。

圖8 三種算法的計算復雜度對比曲線

為驗證在噪聲影響下本文所提算法的魯棒性，在設置同一墻體參數下運用本文算法以及上述兩種算法對穿墻場景進行成像，并給出了不同信噪比條件下該算法的輸入輸出目標雜波比(Target Clutter Ratio,TCR)。由圖9可知，相較于自聚焦算法，基于折射波靜校正的合成延遲時算法在不同信噪比的噪聲情況下能夠輸出更高的目標雜波比；相較于回波域濾波方法，本文所提算法在噪聲信噪比為-3 dB之后能輸出更高的目標雜波比。

圖9 噪聲影響下TCR的變化

4 結 論

針對未知參數的多層墻體穿墻雷達成像目標散焦或位置偏移問題，本文提出基于折射波合成延遲時靜校正算法對電磁波的異常傳播時延進行校正并成像。該算法首先讀取回波信息建立超定方程組計算外層墻體內電磁波的傳播時延，再通過合成延遲時增量估計內層墻體每一通道的電磁波傳播時延，最后對每一通道電磁波的異常傳播時延進行一一校正，使得成像目標被校正到準確位置。仿真結果表明，本文所提算法能夠有效校正墻體存在聚焦時延的偏差問題。與傳統方法相比，本文所提方法能夠更加精確地將目標點校正至正確位置，且圖像聚焦效果優于傳統算法。理論推導以及程序仿真結果表明，本文方法在未知墻體參數穿墻成像的實際應用中具有一定的參考和指導意義。