基于加權(quán)相關(guān)的探地雷達(dá)后向投影成像算法*

陳鑫澎 李靜霞 * 劉 麗徐 航王冰潔張建國(guó)

(1.太原理工大學(xué)新型傳感器與智能控制教育部重點(diǎn)實(shí)驗(yàn)室，山西 晉中 030600；2.太原理工大學(xué)物理與光電工程學(xué)院，山西 晉中 030600)

地下管線(xiàn)在城市中扮演著能量輸送、物質(zhì)傳輸和信息傳遞的角色，被譽(yù)為城市的“生命線(xiàn)”[1－2]。然而，當(dāng)今許多城市由于缺乏對(duì)地下管線(xiàn)的有效管理，造成對(duì)地下管線(xiàn)分布和埋設(shè)深度信息等的誤判，導(dǎo)致在工程施工過(guò)程中事故頻發(fā)，甚至造成爆炸等嚴(yán)重事故，因此地下管線(xiàn)的精準(zhǔn)探測(cè)顯得尤其重要[3－5]。在各種管線(xiàn)探測(cè)方法中，探地雷達(dá)因具有測(cè)量精度高、操作簡(jiǎn)單等優(yōu)點(diǎn)被廣泛應(yīng)用。通常探地雷達(dá)使用一維或二維回波來(lái)完成目標(biāo)識(shí)別，這種工作方式對(duì)操作人員的專(zhuān)業(yè)技能和先驗(yàn)知識(shí)提出了較高層次的要求，且目標(biāo)檢測(cè)存在處理速度慢、誤判和漏判率高等問(wèn)題。為了降低目標(biāo)檢測(cè)和識(shí)別的難度，探地雷達(dá)成像技術(shù)成為研究熱點(diǎn)。其中，較為常用的成像算法有后向投影(back projection，BP)[6－8]，逆時(shí)偏移(reverse time migration，RTM)[9－10]，距離遷移(range migration，RM)[11－12]，衍射層析(diffraction tomography，DT)[13－15]。利用DT 算法對(duì)探測(cè)目標(biāo)進(jìn)行反演時(shí)，反演過(guò)程過(guò)度依賴(lài)對(duì)探測(cè)環(huán)境的精準(zhǔn)建模以及苛刻的近似條件，使其在實(shí)際使用中受到較大限制。采用傳統(tǒng)的RTM 算法和RM 算法對(duì)分層介質(zhì)進(jìn)行探測(cè)時(shí)，目標(biāo)的距離信息估計(jì)不準(zhǔn)確，容易出現(xiàn)欠聚焦現(xiàn)象。而使用BP 算法對(duì)分層介質(zhì)進(jìn)行探測(cè)時(shí)，BP算法會(huì)對(duì)電磁波在不同介質(zhì)表面發(fā)生的折射現(xiàn)象進(jìn)行精確補(bǔ)償，成像效果好，且該算法計(jì)算過(guò)程簡(jiǎn)單，易于工程應(yīng)用，因此BP 算法在雷達(dá)成像技術(shù)領(lǐng)域得到了廣泛應(yīng)用。

經(jīng)典的BP 算法通過(guò)計(jì)算各通道沿傳播時(shí)延曲線(xiàn)的能量累加值，從而得到目標(biāo)區(qū)域的能量聚焦。然而，這種“延時(shí)－求和”的方法會(huì)在成像結(jié)果中引入較強(qiáng)的偽影。在對(duì)目標(biāo)進(jìn)行探測(cè)時(shí)，該偽影易被誤認(rèn)為目標(biāo)，從而造成誤判。為了抑制偽影，國(guó)內(nèi)外學(xué)者提出了不同的改進(jìn)算法。Foo 和Zetik 等[16－17]提出了基于參考信號(hào)與回波信號(hào)互相關(guān)的BP 成像算法。該算法雖然在一定程度上抑制了偽影，但是對(duì)參考信號(hào)的選擇有要求，因而增加了系統(tǒng)復(fù)雜度。周琳等[18－19]借鑒“延時(shí)－相乘－求和”的算法思想，提出了基于數(shù)據(jù)間互相關(guān)后向投影算法(cross-correlated back projection，CBP)，該算法能抑制部分偽影，且無(wú)需引入額外參考信號(hào)。此外，張海如等[20]人提出了基于自相關(guān)的BP 算法來(lái)抑制偽影，但這種方法需要預(yù)先根據(jù)探測(cè)背景設(shè)定相關(guān)系數(shù)的判別閾值，不便于實(shí)際應(yīng)用。

本文提出了一種基于加權(quán)相關(guān)的后向投影算法(correlation weighted back projection，CWBP)來(lái)抑制偽影。將該算法應(yīng)用于頻率步進(jìn)探地雷達(dá)系統(tǒng)，用于對(duì)地下管線(xiàn)進(jìn)行探測(cè)。進(jìn)而，與經(jīng)典BP 算法和CBP 算法相比較，分析了CWBP 算法的偽影抑制能力。

1 頻率步進(jìn)探地雷達(dá)測(cè)量原理

頻率步進(jìn)探地雷達(dá)采用頻率步進(jìn)信號(hào)作為探測(cè)信號(hào)，該信號(hào)是一種大時(shí)寬帶寬積信號(hào)，基于此信號(hào)的探地雷達(dá)能同時(shí)兼顧分辨率和探測(cè)距離兩方面的要求，因此在探地雷達(dá)領(lǐng)域中得到廣泛運(yùn)用[21－22]。

頻率步進(jìn)雷達(dá)的發(fā)射信號(hào)表示為:

式中:T為脈沖重復(fù)周期，τ為脈沖寬度，f0為初始頻率，Δf為頻率步進(jìn)量，N為頻率步進(jìn)數(shù)。

本振信號(hào)表示為:

若目標(biāo)管線(xiàn)與天線(xiàn)的徑向距離為R，則回波信號(hào)可表示為:

式中:c為光速，ε為介質(zhì)的介電常數(shù)。

回波信號(hào)與本振信號(hào)進(jìn)行混頻后，其輸出信號(hào)為:

對(duì)式(4)進(jìn)行歸一化，之后對(duì)歸一化后的方程進(jìn)行逆傅里葉變換并求模可得到:

式中:S(n)為逆傅里葉變換后的信號(hào)。

根據(jù)式(5)可以得到目標(biāo)的回波時(shí)延，進(jìn)一步可將其轉(zhuǎn)化為目標(biāo)的一維距離信息。

2 加權(quán)相關(guān)后向投影算法

然后，采用CWBP 成像算法實(shí)現(xiàn)探測(cè)區(qū)域的二維成像。

經(jīng)典BP 成像算法應(yīng)用于探地雷達(dá)時(shí)，若合成孔徑位置數(shù)目為M，對(duì)于任一成像點(diǎn)A，其到第i個(gè)合成孔徑位置的時(shí)延τA，i，則得到A點(diǎn)到第i個(gè)孔徑位置處的回波響應(yīng):

式中:N為采用點(diǎn)個(gè)數(shù)；Si，N(t)為預(yù)處理后第i個(gè)孔徑位置處的A-scan 數(shù)據(jù)。

通過(guò)式(7)進(jìn)行疊加，得到A點(diǎn)的成像結(jié)果:

經(jīng)典BP 算法會(huì)在成像結(jié)果中引入較強(qiáng)的偽影。為了抑制偽影，我們提出CWBP 成像算法。與經(jīng)典BP 算法不同的是，任一成像點(diǎn)A 的回波響應(yīng)在進(jìn)行累加前，首先通過(guò)加權(quán)和相關(guān)處理進(jìn)行重構(gòu)。其目的是:通過(guò)加權(quán)處理，增強(qiáng)目標(biāo)區(qū)域的反射能量；之后通過(guò)相關(guān)處理，削弱了非目標(biāo)區(qū)域反射能量。因此，采用CWBP 算法，成像結(jié)果中的偽影得到了極大抑制。CWBP 成像算法具體的步驟為:

(1)對(duì)于任一成像點(diǎn)A，計(jì)算第i個(gè)合成孔徑位置的時(shí)延τA，i，進(jìn)而得到A點(diǎn)到第i個(gè)孔徑位置處的回波響應(yīng)XA，i。

(2)對(duì)成像點(diǎn)A到第i個(gè)合成孔徑位置的回波響應(yīng)XA，i進(jìn)行加權(quán)處理。這里需要兩組長(zhǎng)度相等的數(shù)據(jù)Si，n(t)和SM/2，n(t)。其中，Si，n(t)為Si，N(t)中以t＝τA，i為中心的一組數(shù)據(jù)；SM/2，n(t)為Si，N(t)中以t＝τA，M/2為中心的一組數(shù)據(jù)。根據(jù)Si，n(t)和SM/2，n(t)得到皮爾遜加權(quán)系數(shù):

式中:Cov 函數(shù)表示兩個(gè)向量的協(xié)方差。

利用該加權(quán)系數(shù)，對(duì)A點(diǎn)到第i個(gè)孔徑位置的回波響應(yīng)進(jìn)行加權(quán)，得到加權(quán)后的回波響應(yīng):

重復(fù)上述步驟，得到A點(diǎn)到各個(gè)合成孔徑位置加權(quán)的回波響應(yīng):{HA，1，HA，2，…，HA，M}。

(3)對(duì)A點(diǎn)到各個(gè)合成孔徑位置加權(quán)的回波響應(yīng)進(jìn)行相關(guān)處理。具體地，將A點(diǎn)到各個(gè)合成孔徑位置加權(quán)的回波響應(yīng)兩兩相乘，如表1 所示，然后取上三角部分，再按照式(10)，將各列相加求和，重新構(gòu)造得到A點(diǎn)到各個(gè)合成孔徑位置的一次相關(guān)回波響應(yīng):。一次相關(guān)后，虛假目標(biāo)和雜波干擾得到一定程度的消除。

表1 相乘矩陣示意圖

為了進(jìn)一步消除雜波干擾，將一次相關(guān)的回波響應(yīng)再次兩兩相乘，如表2 所示，然后取上三角部分，再按照式(11)將各列相加求和，最終構(gòu)造得到A點(diǎn)到各個(gè)合成孔徑位置的二次相關(guān)回波響應(yīng):。

表2 相乘矩陣示意圖

通過(guò)上述兩步，得到A點(diǎn)的成像結(jié)果。通過(guò)遍歷所有成像點(diǎn)，完成整個(gè)成像過(guò)程。

3 管線(xiàn)探測(cè)實(shí)驗(yàn)結(jié)果

3.1 實(shí)驗(yàn)裝置

管線(xiàn)探測(cè)的實(shí)驗(yàn)裝置如圖1 所示，將管線(xiàn)埋設(shè)于規(guī)格為2.0 m×1.5 m×1.5 m 的沙池內(nèi)，池內(nèi)裝有干燥河沙，干沙的相對(duì)介電常數(shù)約為3.3。實(shí)驗(yàn)中采用矢量網(wǎng)絡(luò)分析儀(Rohde ＆ Schwarz ZNB40)作為頻率步進(jìn)探地雷達(dá)，其產(chǎn)生的步進(jìn)頻率連續(xù)波的起始頻率為1.8 GHz，終止頻率為5 GHz，中心頻率為3.4 GHz。收/發(fā)天線(xiàn)間距2 cm 并沿同一測(cè)線(xiàn)以相同的方向同時(shí)移動(dòng)，每移動(dòng)5 cm 進(jìn)行一次測(cè)量，每個(gè)探測(cè)位置采集401 個(gè)頻率點(diǎn)數(shù)據(jù)。

圖1 實(shí)驗(yàn)場(chǎng)景圖

3.2 單目標(biāo)探測(cè)結(jié)果

首先對(duì)單根金屬或PVC 管線(xiàn)進(jìn)行了探測(cè)，分析了本文所提出的CWBP 成像算法的偽影抑制效果。實(shí)驗(yàn)中，PVC 管線(xiàn)的直徑為0.20 m，長(zhǎng)度為0.52 m，管壁厚度為3 mm，埋深度為0.60 m；金屬管線(xiàn)的直徑為0.10 m，長(zhǎng)度為0.60 m，管壁厚度為2 mm，埋深度為0.32 m。管線(xiàn)探測(cè)的實(shí)驗(yàn)結(jié)果如圖2 所示。對(duì)于PVC 管線(xiàn)，直接對(duì)回波信號(hào)進(jìn)行經(jīng)典BP 成像算法處理，從圖2(a)中可以看到，成像結(jié)果中存在嚴(yán)重的直達(dá)波干擾。為此，在成像前首先對(duì)回波信號(hào)進(jìn)行均值濾波，從圖2(b)中看到地面直達(dá)波得到了有效的抑制。然而，經(jīng)典BP 成像算法因采用“延時(shí)－求和”的方法實(shí)現(xiàn)對(duì)目標(biāo)的定位，其成像結(jié)果中存在很強(qiáng)的偽影，嚴(yán)重影響了目標(biāo)的識(shí)別。為了消除偽影，我們首先采用文獻(xiàn)[18]提出的CBP 成像算法來(lái)代替經(jīng)典BP 算法，可以看到:成像結(jié)果(圖2(c))中的偽影得到了部分的抑制，殘留偽影的存在仍然會(huì)影響目標(biāo)的識(shí)別。為此，采用本文所提出的CWBP 成像算法，從成像結(jié)果(圖2(d))觀察到偽影幾乎被消除。CWBP成像算法的偽影抑制效果顯著優(yōu)于經(jīng)典BP 成像算法和CBP 成像算法。

圖2 PVC 管線(xiàn)和金屬管線(xiàn)成像結(jié)果

對(duì)于金屬管線(xiàn)的探測(cè)結(jié)果如圖2(e)～(h)所示，結(jié)果與PVC 管線(xiàn)的探測(cè)結(jié)果類(lèi)似。通過(guò)使用均值濾波有效地去除了地表直達(dá)波。此外，對(duì)比經(jīng)典的BP 算法、CBP 算法、CWBP 算法的成像結(jié)果，可以看出CWBP 成像算法能顯著抑制偽影，有利于對(duì)目標(biāo)的識(shí)別。

此外，從分辨率的角度，分別給出了三種算法成像結(jié)果中峰值點(diǎn)位置在水平方向的剖面圖。從圖3和圖4 中看到:不論是探測(cè)金屬管線(xiàn)還是PVC 管線(xiàn)，CWBP 成像算法獲得的剖面圖中，峰值的主瓣寬度均有明顯的縮減。該結(jié)果證明了CWBP 算法相較于經(jīng)典BP 算法和CBP 算法，目標(biāo)周?chē)园旮蓴_和背景雜波得到很好的抑制，目標(biāo)的聚焦能力得到增強(qiáng)，進(jìn)而提高了成像分辨率。

圖3 PVC 管線(xiàn)峰值處水平剖面圖比較

圖4 金屬管線(xiàn)峰值處水平剖面圖比較

進(jìn)而采用積分旁瓣比(integrated side lobe ratio，ISLR)、峰值旁瓣比(peak side lobe ratio，PSLR)、目標(biāo)雜波比(target clutter ratio，TCR)定量分析了三種成像算法對(duì)雜波抑制的能力。

ISRL 定義為[23]:

式中:Etotal和Emain分別表示圖像的總能量和目標(biāo)的主瓣能量。

PSLR 定義為[24]:

式中:Pside和Pmain分別表示旁瓣峰值能量和主瓣峰值能量。

TCR 數(shù)學(xué)表達(dá)式為[25]:

式中:Iq為圖像中像素點(diǎn)幅值，Nt和Nc分別為目標(biāo)區(qū)域At和雜波區(qū)域Ac的像素點(diǎn)數(shù)目。

表3 給出了上述三種成像算法的定量比較結(jié)果，其中ISLR值和PSLR值越小，說(shuō)明旁瓣能量和旁瓣峰值相對(duì)目標(biāo)主瓣能量和目標(biāo)峰值越小，圖像成像效果越好。TCR值越大，說(shuō)明圖像中包含雜波干擾越少，圖像成像效果越好。從表3 中可以看出，無(wú)論探測(cè)管線(xiàn)類(lèi)型如何，CWBP算法的上述三個(gè)指標(biāo)值均遠(yuǎn)優(yōu)于經(jīng)典BP算法和CBP算法。定量分析結(jié)果證明CWBP算法具有優(yōu)秀的偽影抑制效果。

表3 成像算法衡量指標(biāo)

3.3 雙目標(biāo)探測(cè)結(jié)果

在實(shí)際工程應(yīng)用中，地下通常埋設(shè)有多根不同材質(zhì)的管線(xiàn)。因此，本文將金屬管線(xiàn)和PVC 管線(xiàn)同時(shí)埋設(shè)于干沙中進(jìn)行探測(cè)，其示意圖如圖5 所示。其中，金屬管線(xiàn)直徑為0.13 m，長(zhǎng)度為0.60 m，厚度2 mm，埋深0.38 m，水平位置距離0.38 m。塑料管線(xiàn)直徑0.20 m，長(zhǎng)0.60 m，厚度4 mm，埋深0.24 m，水平位置1.28 m，兩者水平相距為0.90 m。

圖5 金屬和PVC 管線(xiàn)埋設(shè)示意圖

成像的結(jié)果如圖6 所示。未經(jīng)預(yù)處理的經(jīng)典BP 成像結(jié)果中，存在嚴(yán)重的直達(dá)波干擾(圖6(a))。對(duì)數(shù)據(jù)進(jìn)行均值濾波后，從6(b)中可以看出，直達(dá)波被有效抑制。然而，目標(biāo)周?chē)€存在較多的偽影。且兩目標(biāo)中間因偽影產(chǎn)生了疊影，造成了虛假目標(biāo)。采用了CBP 算法，其成像結(jié)果(6(c))相較于經(jīng)典BP 成像結(jié)果，偽影得到了部分的抑制。殘留偽影以及疊影仍然會(huì)對(duì)目標(biāo)識(shí)別造成干擾。采用CWBP算法，從圖6(d)中可以看到，偽影及疊影得到了顯著抑制，更加利于目標(biāo)的識(shí)別。

圖6 雙目標(biāo)成像結(jié)果圖

4 結(jié)論

本文提出了一種CWBP 成像算法，從實(shí)驗(yàn)數(shù)據(jù)得到的成像結(jié)果來(lái)看，該算法能顯著抑制偽影。此外，對(duì)于雙目標(biāo)探測(cè)，該算法還能抑制雙目標(biāo)偽影產(chǎn)生的疊影。進(jìn)而，分別采用了ISLR，PSLR，TCR 對(duì)結(jié)果進(jìn)行量化分析，結(jié)果進(jìn)一步證明CWBP 算法在目標(biāo)成像方面，其偽影抑制能力均顯著優(yōu)于經(jīng)典BP算法和CBP 算法。