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

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

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

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

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

1 頻率步進探地雷達測量原理

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

頻率步進雷達的發(fā)射信號表示為:

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

本振信號表示為:

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

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

回波信號與本振信號進行混頻后，其輸出信號為:

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

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

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

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

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

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

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

通過式(7)進行疊加，得到A點的成像結(jié)果:

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

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

(2)對成像點A到第i個合成孔徑位置的回波響應(yīng)XA，i進行加權(quán)處理。這里需要兩組長度相等的數(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ù)表示兩個向量的協(xié)方差。

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

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

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

表1 相乘矩陣示意圖

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

表2 相乘矩陣示意圖

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

3 管線探測實驗結(jié)果

3.1 實驗裝置

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

圖1 實驗場景圖

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

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

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

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

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

圖3 PVC 管線峰值處水平剖面圖比較

圖4 金屬管線峰值處水平剖面圖比較

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

ISRL 定義為[23]:

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

PSLR 定義為[24]:

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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