尹子翔 楊小鵬 蘭 天,2
(1.北京理工大學信息與電子學院雷達技術研究所,北京 100081;2.北京理工大學重慶創(chuàng)新中心,重慶 401120)
近些年來,城市成為了災害救援、反恐行動的主戰(zhàn)場。這些作戰(zhàn)往往需要進入密閉建筑物內(nèi)部,如果能夠在進入建筑物之前就能夠獲取到建筑物內(nèi)部的結構信息,可以減少己方損失、提高作戰(zhàn)行動成功率。因此,研究建筑布局成像,具有重要的現(xiàn)實意義和研究價值[1-2]。
對建筑布局成像,國內(nèi)外許多科研機構都開展了相關研究。文獻[3-4]利用邊緣檢測和Hough 變換在穿墻成像結果中突出顯示建筑物墻壁,但是僅能顯示墻體大致位置,沒有對成像結果進行進一步處理。文獻[5]提出了一種多視角建筑布局成像算法,使用M-N-K檢測器融合成像,能夠獲取建筑物內(nèi)部結構信息。文獻[6]提出了一種基于多方位多尺度的建筑布局成像方法,使用非下采樣Contourlet變換,將變換后的高頻分量和低頻分量采用相應的融合規(guī)則進行融合,可以增強建筑物墻體的輪廓和細節(jié)信息。以上兩種算法完全基于圖像處理的手段,并不考慮墻體參數(shù),導致最終成像結果中墻體位置偏移、厚度展寬。文獻[7]提出了一種適用于穿墻雷達的建筑布局成像算法,該算法通過Radon 變換提取墻體位置來進行墻體補償,但是增強了多徑鬼影等雜波的幅度,最終得到的建筑布局相對粗糙。
針對以上問題,本文提出了一種基于墻體補償與多視角融合的建筑布局成像方法。在圖像融合前進行墻體位置矯正,解決圖像融合中出現(xiàn)的墻體展寬、位置偏移問題;在圖像融合時采用最小值融合,能夠在保留墻體的同時,抑制多徑鬼影等雜波。首先,采用后向投影方法對多個視角的觀測回波進行成像,對成像結果應用Radon變換,通過恒虛警率(CFAR)檢測獲取Radon平面峰值,提取墻體的位置信息;然后,以墻體的前后表面中心為分界線,對圖像進行分割;進一步,對包含偏移墻體的子圖像進行補償,矯正墻體位置;最后,對多個視角墻體補償后的圖像進行最小值融合,抑制多徑鬼影等雜波,得到正確清晰的建筑布局全景圖像。
對于一個穿墻雷達建筑布局成像場景,第k個合成孔徑單元接收到的回波Yk(t)(k=1,2,…,L)由第i(i=1,2,…,N)堵墻的前表面反射波Yki,1(t)和后表面反射波Yki,2(t)組成。圖1為三墻場景下電磁波傳播示意圖。在本文的成像場景中,成像目標為墻體而不是點目標,在成像中起主要作用的是垂直入射波的回波。因此可以忽略電磁波的折射,僅關注電磁波的垂直透射與反射。假設發(fā)射信號為S(t),那么Yk(t)可以寫為
其中,Ri為第i堵墻和第i-1堵墻之間的距離,R1為第一堵墻和收發(fā)天線之間的距離,Di為第i堵墻的厚度,c為光速,為電磁波在第i堵墻中的傳播速度,其值為
其中,εi是第i堵墻的相對介電常數(shù)。
后向投影法(Back Projection,BP)基本思想是首先對成像區(qū)域進行網(wǎng)格化,計算區(qū)域中像素點到收發(fā)天線的距離,進而得到傳輸時延,根據(jù)該時延尋址各個收發(fā)天線對應的回波信號并進行疊加處理。完成整個成像區(qū)域的遍歷運算,得到最終的成像結果[8]。
以成像平面I(X,Y) 中任意一像素點Pq=(xq,yq)為例,該像素點對應的成像延時Tqk為
式中PT,k和PR,k(k=1,2,…,L)分別表示發(fā)射天線和接收天線的位置。將發(fā)射天線坐標PT,k=(xtk,ytk),接收天線坐標PR,k=(xrk,yrk)代入公式(3),像素點Pq的成像聚焦時延可以寫為
將該像素點Pq對應的所有回波信號疊加,得到像素點Pq的成像結果
最后遍歷整個成像區(qū)域,得到最終的成像結果I(x,y)。
實際應用中,墻體大多為平直的,利用這一特性采用直線檢測方法檢測原始圖像I(X,Y)中的直線,獲取墻體的位置信息,對偏移墻體進行矯正。本文采用Radon變換檢測直線來獲取墻體位置。
Radon變換定義如下
其中,I為成像平面,I(x,y)為BP 成像結果的像素值,ρ為直線到坐標原點的距離,θ為點(x,y)與x軸的夾角,δ為Dirac delta函數(shù),其值為
這個函數(shù)使I(x,y)沿著直線
進行積分,如圖2所示。
Radon 變換將圖像投影到(ρ,θ)平面,(ρ,θ)平面中每一點對應(x,y)平面中一條直線的積分。因此,(x,y)中的墻體所在直線會在(ρ,θ)平面形成峰值,這就將原直線的檢測轉化為對峰值的檢測。使用CFAR 檢測得到峰值位置,提取峰值對應的(ρ,θ)值,得到(x,y)平面中該直線的斜率和與圖像中心的斜距,從而確定該直線即墻體位置。
由公式(4)計算得到的聚焦時延未考慮電磁波在墻體和在空氣中傳播的速度的差異,因此需在不同墻體及所在區(qū)域上矯正上述差異帶來的位置偏移。
設有圖1所示場景,墻體位置矯正的步驟如下:
1)對回波數(shù)據(jù)進行BP 成像得到原始圖像I1(x,y)。
2)對I1(x,y)應用Radon 變換得到R1(ρ,θ),通過CFAR 檢測,得到的第一峰值和第二峰值分別對應I1(x,y)中第一堵墻的前表面和后表面。二者位置相減得到第一堵墻的成像厚度d1,其與真實厚度D1的關系為
3)以第一堵墻中心為分界線將圖像I1(x,y)分割為兩個子圖像I11(x,y)和I12(x,y)。子圖像I11(x,y)包含第一堵墻前表面及之前的部分,子圖像I12(x,y)包含I1(x,y)中剩下的部分。其中,I11(x,y)中墻體位置正確,I12(x,y)中墻體位置偏移。
4)對第一堵墻造成的時延進行補償,這個補償時延為
補償后I12(x,y)實際的成像結果應為
式(11)中對子圖像每個像素點補償一個常數(shù)時延Δτ1,等效于將該子圖像向天線方向平移Δτ1·c=d1-D1的距離。因此無需對該子圖像重新BP 成像,直接將該子圖像向天線方向平移d1-D1,得到補償后的圖像I2(x,y),即能達到與重新BP 成像相同的效果。對于公式(10)中d1可由步驟2)得到。D1為第一堵墻的真實厚度,可在真實場景中由相鄰墻體厚度估計得到。通過這種方式,我們無需提前得知墻體的相對介電常數(shù)ε1便可對該墻體進行補償。
5)對I2(x,y)應用Radon變換,得到第二堵墻的成像厚度d2,以第二堵墻中心為分界線將I2(x,y)分割為兩個子圖像I21(x,y) 和I22(x,y)。子圖像I21(x,y)包含第一堵墻的后表面、第二堵墻的前表面以及它們之前的部分,子圖像I22(x,y) 包含I2(x,y)剩下的部分。其中,I21(x,y)中墻體位置正確,I22(x,y)中墻體位置偏移。
6)對第二堵墻造成的時延進行補償,這個補償時延為
將I22(x,y)向天線方向平移(d1-D1)+(d2-D2),得到補償后的圖像I3(x,y)。
7)對I3(x,y)應用Radon變換,以第三堵墻中心為分界線分割I3(x,y),并對靠后的子圖像進行補償,得到I31(x,y)和補償后的圖像I4(x,y)。
8)將所有包含正確墻體位置的子圖像I11(x,y),I21(x,y),I31(x,y),I4(x,y)分別進行歸一化并拼接,得到完整的圖像I'(x,y)。
穿墻雷達探測的場景大多是封閉的室內(nèi)環(huán)境,電磁波除了產(chǎn)生墻體前后表面回波外,還會產(chǎn)生墻體間多徑回波、墻體內(nèi)多徑回波等。多徑回波會導致成像中出現(xiàn)鬼影,嚴重影響對建筑布局的判斷。
Q.Tan等研究發(fā)現(xiàn)多徑鬼影對雷達探測位置具有方向依賴性(Aspect Dependence),即多徑鬼影位置會隨著雷達觀測位置的變化而改變[9]。利用這一特性,從多個視角對建筑進行觀測,通過對多個觀測結果取最小值來抑制多徑鬼影。
基于多視角融合的多徑鬼影抑制算法原理如下:
(1)多視角觀測下墻體成像位置不變
如圖3 所示,從視角1 和視角2 對該墻體的左墻面進行成像。由于視角1 發(fā)射的電磁波僅在空氣中傳播,經(jīng)過一次反射,接收到信號,因此墻體左表面成像位置就是真實位置。從視角2 觀測,由于電磁波要穿透墻體,電磁波的傳播速度發(fā)生改變,成像位置會變?yōu)閳D中左側虛線位置。但通過第4 節(jié)墻體補償后,墻體位置被矯正到圖中右側虛線位置,即真實位置。因此,對同一個墻面,從視角1 與視角2 進行觀測得到的墻體成像結果一致且均為真實位置。此為多視角觀測下目標墻體成像位置不變。
(2)多視角觀測下墻體內(nèi)多徑鬼影成像位置不同
從視角1 與視角2 對墻體進行觀測得到的墻體內(nèi)多徑鬼影如圖4 所示。視角1 發(fā)射的電磁波在墻體內(nèi)進行了三次反射,同時電磁波在墻體內(nèi)的傳播速度也會發(fā)生改變,因此該路徑下多徑回波產(chǎn)生的多徑鬼影位置如圖中右側虛線所示。同理,由視角2 觀測得到的墻體內(nèi)多徑鬼影如圖中左側虛線所示。綜上可以發(fā)現(xiàn),對同一個墻體由視角1 與視角2觀測得到的墻體內(nèi)多徑鬼影位于墻體兩側。此為多視角觀測下墻體內(nèi)多徑鬼影成像位置不同。
(3)多視角觀測下墻體間多徑鬼影成像位置不同
從視角1 與視角2 對墻體進行觀測得到的墻體間多徑鬼影如圖5 所示。視角1 發(fā)射的電磁波在墻體間進行了三次反射,該路徑下多徑回波產(chǎn)生的多徑鬼影位置如圖中右側虛線所示。同理,由視角2觀測得到的墻體內(nèi)多徑鬼影如圖中左側虛線所示。因此,對同一組墻體,由視角1 與視角2 觀測得到的墻體間多徑鬼影位于該組墻體的兩側。此為多視角觀測下墻體間多徑鬼影成像位置不同。
(4)多徑鬼影的方位依賴性
由前三點可以得出結論,在應用第4 節(jié)墻體補償?shù)那疤嵯?,多視角觀測下墻體的成像位置不變,墻體內(nèi)多徑鬼影、墻體間多徑鬼影成像位置不同,即多徑鬼影的方位依賴性。基于這一特性,可以從多個視角對建筑進行觀測,對多個觀測結果取最小值來抑制多徑鬼影,同時墻體信息得到保留。多視角觀測示意圖如圖6所示。
假設有四個觀測視角,每個觀測視角得到的墻體補償完畢后的圖像為(xk,yk)(k=1,2,3,4)。由于每一個視角的成像結果都有自己對應的坐標系,因此在圖像融合時需要統(tǒng)一各個視角的坐標系。假設全局坐標系為圖6 所示的XOY坐標系,利用下式將局部成像結果轉換為全局成像結果。
其中,Ik(x,y)為轉換后的全局坐標系成像結果,rot為轉換操作符,ψk=[90° 0° 270° 180°]為轉換角度。
將多個視角的觀測信息進行融合,得到
對每個像素點取所有圖像對應位置的最小值,融合后成像結果僅在墻壁處有較高幅值。因為同一個多徑回波不是所有觀測視角都存在,所以取最小值后多徑回波得到抑制,而且該方法不需要場景的先驗信息即可有效的抑制多徑鬼影。最終的I(x,y)為矯正墻體位置并抑制多徑鬼影的建筑布局圖像。
以三層墻體為例,墻體補償與多視角融合方法流程圖如圖7所示。
下面通過仿真驗證本文所提墻體補償與多視角融合方法成像效果。使用gprMax[10-11]對一個房間模型進行仿真,房間結構布局如圖8所示。建筑的墻體為0.24 m厚的磚砌成,相對介電常數(shù)為3.8。此外,模型中具有5 cm 厚的混凝土地板與天花板,相對介電常數(shù)為6.8。整個房間的大小為3 m×3 m×1 m。
發(fā)射信號為gaussiandot波形,中心頻率1.0 GHz,通過偶極子天線發(fā)射。依次從四個視角對建筑物進行觀測,每次觀測時天線距離墻壁0.4 m,天線步進0.1 m,步進次數(shù)27 次,發(fā)射天線與接收天線距離0.2 m。
墻體補償?shù)姆抡娼Y果如圖9 所示。圖9 中,圖9(a)為沿觀測視角1 去掉直達波后的BP 成像。在此圖中,第一堵墻前表面位置正確,其他墻體位置偏移。圖9(b)為對圖9(a)旋轉90°后進行Radon變換后的結果,運用CFAR 檢測得到第一堵墻的位置與厚度信息。圖9(b)中,θ=90°列中的峰值對應圖9(a)中的墻面。圖9(c)和圖9(d)分別為補償?shù)谝欢聣εc第二堵墻后的結果。在圖9(c)中,第一堵墻后表面與第二堵墻前表面位置正確;在圖9(d)中,第二堵墻后表面與第三堵墻前表面位置正確,其余墻體位置偏移。圖9(e)為補償?shù)谌聣蟮膱D像,在此圖中,第三堵墻后表面得到了矯正。圖9(f)為所有墻體位置正確的子圖像歸一化后拼接而成的完整圖像,在此圖中,所有墻體的位置均正確。
多視角融合的仿真結果如圖10 所示。圖10中,圖10(a)為對觀測視角1 進行墻體補償后的圖像。圖10(b)為對觀測視角3 進行墻體補償后的圖像。圖10(c)為圖10(b)與圖10(a)取最小值后的結果,即視角1 與視角3 融合后的結果。圖10(d)為對圖10(c)進行灰度線性變換增加對比度后的圖像,在此圖中,能夠明顯辨認三堵墻的位置,并且多徑鬼影相比圖10(a)、(b)顯著減少。對觀測視角2、4 進行相同的操作,得到圖11 所示結果。
圖12 給出了不同建筑布局成像方法的比較,圖12(a)為對四個視角觀測回波BP 成像后直接相加的結果,圖12(b)為文獻[5]所提方法得到的結果,圖12(c)為文獻[7]所提方法得到的結果,圖12(d)為本文算法得到的結果??梢钥闯?,與算數(shù)融合方法和文獻[5]所提的M-N-K 方法相比,本文方法可以還原墻體的真實厚度與位置。文獻[7]的方法雖然能還原墻體位置與厚度,但是增強了多徑鬼影等雜波的強度,本文算法與之相比抑制了多徑鬼影,改善了成像效果。
本文針對穿墻雷達建筑布局成像中的墻體展寬、位置偏移和多徑鬼影問題,提出了一種基于墻體補償與多視角融合的穿墻雷達建筑布局成像方法。首先,建立了建筑布局成像的回波模型,推導了BP 成像算法;其次,針對墻體位置偏移問題,采用墻體位置檢測并補償聚焦時延的方法矯正墻體位置;之后,針對多徑鬼影問題,采用多視角觀測融合的方法抑制多徑鬼影。仿真結果表明,與傳統(tǒng)的建筑布局成像方法相比,本文所提方法能解決傳統(tǒng)方法中的墻體展寬和位置偏移問題,并且能在圖像融合的同時抑制多徑鬼影,得到質量更高的建筑布局全景圖。