李家強, 陳德昌, 陳金立, 朱艷萍
(1. 南京信息工程大學(xué)電子與信息工程學(xué)院, 江蘇南京 210044;2. 南京信息工程大學(xué)氣象災(zāi)害預(yù)報預(yù)警與評估協(xié)同創(chuàng)新中心, 江蘇南京 210044)
穿墻成像雷達(through the wall radar imaging)通過發(fā)射電磁波對墻體后靜止或運動的目標進行探測與成像,廣泛地應(yīng)用于救援任務(wù)、領(lǐng)土安全、城市監(jiān)控、軍事偵察等方面[1-2]。在穿墻雷達成像過程中,很多學(xué)者對前墻的強雜波進行了深入的研究,并對其進行了有效抑制。Solimene等提出一種基于熵值的算法[3],能夠有效地抑制前墻雜波,該算法是利用墻體與目標的差異來設(shè)置適宜的門限值,濾除前墻雜波。然而由于實際中側(cè)墻的存在,接收的雷達信號不僅包括目標與側(cè)墻間多徑傳播產(chǎn)生的回波信號,也包括側(cè)墻與后墻之間反射產(chǎn)生的強雜波信號,這必將導(dǎo)致多徑虛像與墻體雜波的產(chǎn)生,使得目標定位失準、真假目標難辨,降低穿墻雷達檢測與成像性能,因此如何降低多徑虛像與強雜波的干擾是很有意義的。
Setlur等[4]采用高斯加權(quán)函數(shù)的方法抑制多徑虛像,該方法首先利用最小二乘法來定位虛像的位置,然后對于目標與虛像的位置分別設(shè)定不同的加權(quán)函數(shù),使目標的像素值保持不變的同時將虛像的像素值置零。文獻[5]提出用子孔徑成像的方法來抑制多徑虛像,分別找出不能接收左右兩面?zhèn)葔w多徑回波的子孔徑,通過后向投影算法形成兩幅子圖像,最后再與全孔徑成像圖連乘,根據(jù)這三幅圖像多徑虛像的差異性來去掉左右墻的兩個多徑虛像。上述方法在虛像抑制之前都通過基于先驗的背景相消[6]的方法消除了前后墻體雜波以及側(cè)墻與后墻反射形成的強雜波,然而在實際工程應(yīng)用中是無法直接測量墻體后的背景環(huán)境,所以對封閉墻體室內(nèi)目標成像前預(yù)先濾除墻體雜波仍是穿墻雷達一個值得研究的方向。其中經(jīng)典的去除墻體雜波的方法是奇異值分解法(Singular Value Decomposition,SVD)[7],其主要通過墻體雜波與目標回波信號能量之間的大小差異來確定墻體子空間,進而通過墻體子空間投影方法抑制墻體雜波,然而當側(cè)墻存在時,側(cè)墻與后墻反射產(chǎn)生的墻體雜波能量與目標接近且分布不均勻,因此難以濾除墻體雜波。
為了實現(xiàn)更為實用的方法對墻體后的目標進行檢測,獲得精度更高的成像,本文提出一種基于伽馬校正的子孔徑圖像相乘的方法,該方法無需對墻體場景進行背景相消也能濾除虛像與墻體雜波。本文通過仿真實驗?zāi)M一個無目標的室內(nèi)成像結(jié)果圖分析得出兩面?zhèn)葔εc后墻兩個強雜波產(chǎn)生與目標多徑無關(guān),然后基于子孔徑成像相乘原理,抑制兩面?zhèn)葔Φ膬蓚€虛像與前后墻體雜波。同時采用伽馬校正的方法,根據(jù)每幅圖像對比度因子的比例關(guān)系,得到每幅子圖像指數(shù)權(quán)值,并對指數(shù)加權(quán)后的子圖像進行連乘,以消除側(cè)墻與后墻反射產(chǎn)生的強雜波,最后用閾值法濾除后墻產(chǎn)生的虛像殘余。
本文基于FDTD的電磁仿真軟件GprMax2D/3D[8],對穿墻雷達場景進行建模,探測模型如圖1所示,設(shè)墻體由單層均勻介質(zhì)組成,厚度為d,相對介電常數(shù)為εr,電導(dǎo)率為σ。目標為一理想的球狀電導(dǎo)體,它的半徑為r,圓心與前墻體的垂直距離為δ,全向天線與前墻體之間的垂直距離為h,沿著平行于墻體的測量線等間距掃描N次,發(fā)射天線發(fā)射超寬帶窄脈沖信號,接收天線在相同的位置接收回波信號。發(fā)射信號為Ricker子波,即一階高斯脈沖信號取負,其表達式為
(1)
天線接收到的回波信號模型表示為
e(t)=et(t)+ea(t)+ew(t)+em(t)
(2)
式中,et(t)為目標回波信號,ea(t)為天線失配噪聲,ew(t)為前后墻體雜波信號與側(cè)墻和后墻之間反射產(chǎn)生的強雜波,em(t)為由目標多徑產(chǎn)生的多徑回波信號。在本文中,假定全向天線為理想天線,因此不存在失配噪聲,而式(2)中后兩項是需要濾除和抑制掉的墻體間雜波和多徑虛像信號,即為了得到干凈的目標信號et(t),必須通過算法濾除墻體雜波信號ew(t)與目標多徑回波信號em(t)。
本文將成像區(qū)域約束在封閉的墻體中,如圖2所示,側(cè)墻與前墻的長度分別為D1和D2,墻體的厚度與介電常數(shù)分別為d和εr,目標位于P(x,y),每根天線都是單發(fā)單收型,且沿著前墻平行的直線排列,每相鄰兩根天線的距離相等,共N根天線。第n根天線的位置是(xn,yn),n=1,2,…,N。
圖2 多徑傳播模型
(3)
式中,c為電磁波的傳播速度。設(shè)s(t)為發(fā)射信號,則得到目標回波信號與目標多徑回波信號分別為
(4)
(5)
式中,Γpn為與傳播路徑path-p相關(guān)的目標散射參數(shù)。式(4)代表發(fā)射信號由第n根天線沿路徑A到目標,再由目標沿此路徑返回到天線,式(5)代表發(fā)射信號沿路徑A到目標,再由路徑B,C,D返回到天線。
本文提出的算法由3個步驟完成:第一步,找到兩個合適的子孔徑,即既要找到盡可能多的天線保證成像分辨率,又要確保只能接收單邊側(cè)墻的虛像;第二步,用兩個子孔徑圖像與全孔徑圖像連乘來抑制側(cè)墻的兩個虛像與前后墻體雜波;第三步,根據(jù)伽馬校正的方法對每幅圖像設(shè)置一個指數(shù)權(quán)值并連乘以消除側(cè)墻與后墻反射產(chǎn)生的強雜波。
如圖2所示,以左墻為例:當作P點與左墻鏡像對稱點時,此對稱點與前墻左端連線的延長線與天線有一個交點,這個交點距離天線左端有一段距離,根據(jù)電磁波反射的幾何關(guān)系可知在這段距離中排放天線無法接收左墻的目標多徑回波,然而在實際測量中無法預(yù)先探知目標的位置。選取子孔徑的方法如下:天線從左端第1根開始排放天線,并進行BP(back projection)后向投影算法,然后依次累加一根天線,當?shù)趉根天線開始能接收到左墻虛像時,則第1根到第k-1根天線為子孔徑1,同理可以得到子孔徑2。
假設(shè)有N幅子圖像,都是由BP后向投影成像算法產(chǎn)生,每個子圖像的行列像素都相同,第i幅子圖像表示為Pi,對其進行歸一化處理:
Ii(x,y)=Pi(x,y)/max?(x,y){Pi(x,y)}
(6)
式中,Ii(x,y)為圖像Ii的第x行第y列的歸一化像素值。將N幅圖像進行相乘融合[9],即為歸一化幅度圖像的逐像素乘積,設(shè)融合后的圖像像素值為If(x,y)。
If(x,y)=I1(x,y)*I2(x,y)…*IN(x,y)
(7)
鑒于本文多徑回波模型取N=3,其中I1是全孔徑成像圖。I2無法接收左墻虛像,即子孔徑1成像圖,I3無法接收右墻虛像,即子孔徑2成像圖。在每個子圖像中,雜波分布是不同的,左右兩面?zhèn)葔δ繕颂撓竦奈恢靡彩窍喈惖模悄繕说奈恢檬窍嗤蛔兊?,所以用子孔徑圖像相乘可以去除側(cè)墻的兩個虛像與前后墻體雜波。由于兩個側(cè)墻與后墻產(chǎn)生的強雜波能量與目標接近,所以上述的方法處理后仍然有雜波殘余,可進一步采用伽馬校正[10]的方法進行消除。
在圖像處理中,對圖像進行伽馬校正(又稱冪律變換)的修訂,能夠使得圖像信息的對比效果更加明顯,伽馬校正的基本公式為
s=c(r+ε)g
(8)
式中:s表示圖像的輸出數(shù)據(jù);c表示常量系數(shù),通??梢詫⑵湓O(shè)為1;r表示輸入圖像的輸入數(shù)據(jù);ε表示為了有一個可適度的輸出值而加入的一個偏移量。然而在一般情況下不在計算范圍內(nèi),因此可以省略偏移量ε,于是公式可以簡化為
s=crg
(9)
其中,g決定了像素變換的方式,g=1時代表線性變換,不會影響圖像對比度;g≠1時代表非線性變換,會導(dǎo)致圖像對比度發(fā)生改變。將常量系數(shù)c設(shè)為1,根據(jù)指數(shù)與對數(shù)的關(guān)系可知:
(10)
所以當指數(shù)權(quán)值g越大時,根據(jù)對數(shù)函數(shù)的性質(zhì),輸入像素值r越大的區(qū)域有更寬的輸出值范圍,而像素值越小的區(qū)域則有更窄的輸出值范圍,所以在穿墻雷達圖像中,一般來說目標區(qū)域的像素值最大,而其他雜波背景像素值較小,所以用伽馬變換可以讓穿墻雷達圖像中的目標區(qū)域更加鮮明。
設(shè)有N幅子圖像,每個子圖像賦予值為gi(i=1,2,3,…,N)的指數(shù)權(quán)值,并用一個圖像評估標準(NSSI)來約束它,把這個標準稱為圖像對比度[11],圖像對比度與指數(shù)權(quán)值的關(guān)系為
(11)
gn=NSSIn/NSSImin
(12)
在所有子圖像中對比度最低的gmin=1,其他子圖像gi值通過比例關(guān)系算出,不同的子圖像像素被賦予不同的指數(shù)權(quán)值然后相乘如下:
(13)
Iλ(x,y)=
(14)
式中,λth表示一個預(yù)定義的閾值。
為了方便比較虛像與雜波抑制的效果,現(xiàn)定義目標雜波比(Target-to-Clutter Ratio,TCR)如下[12]:
(15)
式中,I(·)為成像點所對應(yīng)的歸一化幅度值,At和Ac分別為目標區(qū)域和雜波區(qū)域(虛像也作雜波處理),Nt和Nc為各區(qū)域所對應(yīng)的成像點數(shù)目。
根據(jù)圖1所示的系統(tǒng)模型,建立仿真場景:全向天線陣列沿平行于墻體的測量線等間距進行26次掃描,且每根天線與前墻存在一段距離,天線的橫向掃描范圍為0.1~2.1 m。墻體是均勻介質(zhì),其厚度d=0.2 m,相對介電常數(shù)εr=6.4。目標為球狀的理想電導(dǎo)體,其球心距離前墻體δ=1 m,半徑r=0.1 m。通過GprMax2D/3D軟件獲取各個測量位置的回波信號,使用BP后向投影成像算法進行全孔徑成像,結(jié)果如圖3所示。
圖3 全孔徑成像圖
圖3中存在8個重點區(qū)域,如箭頭所指。除了目標外,其余分別為虛像與墻體雜波。根據(jù)虛像定位的基本理論可知M1,M2,M3分別由左墻、右墻、后墻的目標多徑回波產(chǎn)生的虛像。通過空場景成像圖4所示,C1與C2是由側(cè)墻與后墻的反射產(chǎn)生的雜波。
基于SVD去除四面墻體模型的墻體雜波時,如圖5所示,該算法失效,而圖6的基于背景相消法能夠有效地濾除墻體雜波,但是該方法需要事先獲取室內(nèi)空場景回波信息,實際情況很難滿足。
結(jié)合本文設(shè)定的模型,通過子孔徑選取計算得到子孔徑1為第1根到第8根天線,子孔徑2為第19根到第26根天線,子孔徑1,2的成像結(jié)果如圖7、圖8所示,圖7與圖8中分別抑制了虛像M1與M2,同時仿真實驗表明目標圖像始終存在,且位置不發(fā)生改變。
圖4 空目標成像圖
圖5 SVD去墻體雜波
圖6 背景相消去墻體雜波
圖7 子孔徑1成像
圖8 子孔徑2成像
將兩個子孔徑圖像與全孔徑圖像直接相乘,得到結(jié)果如圖9所示,虛像M1與M2及前后墻體雜波得到了有效抑制,但是仍然有大量的強雜波殘余C1和C2。
圖9 子孔徑圖像相乘
為了解決上述雜波殘余的問題,采用了伽馬校正的方法,即對每個子圖像設(shè)置一個指數(shù)權(quán)值gi,gi的大小計算方法如上文介紹,根據(jù)得到的3個子圖像的像素矩陣,通過比例關(guān)系得出g1=1,g2=4.11,g3=4.03,將三幅圖像設(shè)置對應(yīng)指數(shù)權(quán)值后連乘,并設(shè)置閾值λth=0.5得到如圖10,從圖中可以看出C1,C2及M3得到完全抑制。
圖10 伽馬校正后用閾值法處理
為了驗證本文算法的性能,給出目標位于不同位置時,傳統(tǒng)子孔徑圖像相乘方法與伽馬校正之后的信雜比,如表1所示,當目標位置不同時,信雜比增量最小為24 dB。
表1 不同橫坐標位置時輸出目標圖像TCR比較
由于在實際系統(tǒng)中噪聲是不可避免的,因此需要在噪聲的影響下評估提出算法的性能。圖11給出在回波信號中加入不同信噪比時,原始成像、子孔徑圖像相乘、伽馬校正之后目標區(qū)域TCR的變化,從圖11可知通過伽馬校正能提高TCR。圖12給出伽馬校正算法中3個指數(shù)權(quán)值隨加入噪聲的變化,可知g1值恒定不變?yōu)?,而g2與g3的值隨著SNR增加而增大,這是由于SNR增加,圖像對比度也隨著增大,由伽馬校正的公式可知每個指數(shù)權(quán)值也會相應(yīng)變大,當SNR很小時,3個指數(shù)權(quán)值都接近1。
圖11 加入噪聲時信雜比的變化
圖12 加入噪聲時指數(shù)權(quán)值的變化
為了解決穿墻雷達中目標與側(cè)墻、側(cè)墻與后墻之間的多徑虛像和雜波問題,首先對其產(chǎn)生原因進行了回波仿真實驗分析。根據(jù)試驗分析結(jié)果,提出基于伽馬校正的子孔徑圖像相乘的方法,即通過子孔徑圖像相乘的方法消除由側(cè)墻產(chǎn)生的虛像與前后墻體雜波,并以此為基礎(chǔ)采用伽馬校正的方法,計算出每個子圖像的圖像對比度。然后根據(jù)三幅圖像對比度比例關(guān)系,將每個子圖像設(shè)置相應(yīng)的指數(shù)權(quán)值,并進行連乘,有效地去除了側(cè)墻與后墻之間的強雜波。從仿真結(jié)果可以看出信雜比有了明顯的提高,充分說明了本文所提方法對于實際工程中強雜波背景下穿墻雷達多徑虛像抑制具有一定的理論指導(dǎo)意義。