彭 建，鐘景陽，白 潔，孟 旭，王洪華，劉 海

(1. 中國建筑第八工程局有限公司，上海 200135；2. 廣州大學(xué) 土木工程學(xué)院，廣東 廣州 510006；3. 桂林理工大學(xué) 地球科學(xué)學(xué)院，廣西 桂林 541006)

0 引言

探地雷達(Ground penetrating radar，GPR)是一種高效、高分辨率的淺層地球物理勘探方法，其原理是利用高頻寬帶電磁波(10 MHz～10 GHz)對地下目標(biāo)或者物體內(nèi)部進行探測。通過發(fā)射天線向待探測目標(biāo)發(fā)射高頻電磁波，之后利用接收天線采集回波信號，并分析其雙程走時、振幅和相位等信息來解析目標(biāo)體的空間分布及物理性質(zhì)。探地雷達方法相比其他地球物理勘探方法具有操作簡單、分辨率高、實時成像、適應(yīng)性強等優(yōu)點，目前已被廣泛應(yīng)用于土木工程檢測[1-3]、深空探測[4]、地質(zhì)調(diào)查、考古勘探[5]等領(lǐng)域。

在實際探測中，探地雷達接收到的信號不僅包含地下目標(biāo)產(chǎn)生的回波信號，還存在著各種雜波，如直達波、耦合波、環(huán)境干擾和隨機噪聲等。目標(biāo)信號往往會被雜波信號掩蓋，增加了對目標(biāo)體的辨識難度。雜波信號的形成有內(nèi)部因素和外部因素，內(nèi)部因素是雷達采集系統(tǒng)自身引起的，而外部因素是由于外界環(huán)境引起的，比如，探地雷達發(fā)射的探測信號從空氣進入地下時，由地表面引起的雜波，由于地下介質(zhì)非均勻隨機分布引起的地下介質(zhì)雜波，同時發(fā)射天線和接收天線之間還存在天線互耦雜波，雷達系統(tǒng)在接收和存儲數(shù)據(jù)時也會引入一定的噪聲等[6]。

雜波抑制是探地雷達數(shù)據(jù)處理的必要步驟，目的是去除雷達數(shù)據(jù)中的非目標(biāo)成分，使目標(biāo)信號得到增強。目前常見的探地雷達數(shù)據(jù)雜波抑制方法大致可分為基于雜波模型、基于信號統(tǒng)計和基于子空調(diào)投影三類方法。其中，基于雜波模型的減平均道法[7]和Scale and Shift法[8]較為常用，但由于計算方法簡單，雜波抑制效果不理想。基于信號統(tǒng)計的最大似然估計法[9]、卡爾曼濾波法[10]，以及Dahan Liao[11]提出的統(tǒng)計極化法、Fawzy Abujarad[12]提出的獨立特征分析法等算法僅對抑制地表不平整時產(chǎn)生的雜波有較好的效果，而且這些算法的穩(wěn)定性不高，對所使用的參數(shù)比較敏感?；谧涌臻g投影的奇異分解法[13]和主成分分析法[14]較其他算法具有較強的雜波抑制能力，但算法參數(shù)的選取尤為關(guān)鍵。近年來有學(xué)者提出使用背景矩陣減法對視頻和圖片中的噪聲進行處理，獲得了較好的結(jié)果[15]。本文將背景矩陣減法應(yīng)用于探地雷達信號雜波抑制，并將減平均道法、奇異分解法、背景矩陣減法以上三種方法分別對鋼筋混凝土結(jié)構(gòu)探測和地下管線探測兩種情形的仿真數(shù)據(jù)和一組實測數(shù)據(jù)進行雜波抑制處理，并對結(jié)果展開對比分析。

1 雜波抑制方法

1.1 減平均道法

減平均道法是探地雷達數(shù)據(jù)處理中最常用的雜波抑制算法。這種算法假設(shè)雜波信號沿測線方向基本不變，即具有大致相同的波形、振幅和時延。通過對構(gòu)成 B-scan 數(shù)據(jù)的每道 A-scan雷達數(shù)據(jù)沿橫向進行平均，得到平均道，之后將每道 A-scan 數(shù)據(jù)減去平均道來消除雜波信號。假設(shè)探地雷達的 B-scan 數(shù)據(jù)矩陣為A(A為M行N列的矩陣)，S表示雜波抑制后得到的目標(biāo)信號矩陣估計，減平均道法見數(shù)學(xué)公式(1)：

i=1，2，…，M；j=1，2，…，N

(1)

1.2 奇異分解法

基于子空間投影的奇異分解法是將矩陣分解為3個矩陣相乘的形式，其中分解的矩陣中包括2個正交矩陣和1個對角陣。

對于M×N的雷達數(shù)據(jù)，M代表采樣點數(shù)，N代表道數(shù)，雷達數(shù)據(jù)的矩陣表示形式見公式(2)：

A=(xi，j)(i=1，2，…，M；j=1，2，…，N)

(2)

式中：xi，j代表在第i行，第j列的值。

矩陣A奇異分解表示可寫成

A=UDVT=[u1…ur…um]·

(3)

式中：U是M×M的正交矩陣，其每列為協(xié)方差矩陣AAT特征向量ui，V是N×N正交矩陣，其每列為協(xié)方差矩陣ATA的特征向量vi，D是由r(r是A的秩)個元素組成的M×N的對角矩陣，其元素從大到小排列，即σ1≥σ2≥…≥σr≥0，σi稱為矩陣A的奇異值。將式(3)展開得，

(4)

(5)

式中：k是一個小于維度M和N的正整數(shù)。將式(5)代表的雜波信號去除即得到處理后的雷達數(shù)據(jù)。

1.3 背景矩陣減法

背景矩陣減法將探地雷達數(shù)據(jù)作為矩陣處理，矩陣的列數(shù)對應(yīng)雷達數(shù)據(jù)的道數(shù)，行數(shù)對應(yīng)于雷達數(shù)據(jù)的采樣點數(shù)。背景矩陣減法通過計算僅由雜波組成與雷達數(shù)據(jù)矩陣相同尺寸的背景矩陣，然后將其從原始雷達信號矩陣中減去，達到抑制雜波信號的目的。

背景矩陣減法的實現(xiàn)過程包括以下步驟：

1)在GPR剖面數(shù)據(jù)矩陣每個采樣點布置一維水平窗口，窗口長度取決于需要處理的雷達數(shù)據(jù)道數(shù)；

2)計算窗口內(nèi)雷達數(shù)據(jù)的平均值，見公式(6)、(7)：

(6)

L=N×α

(7)

式中：N為一維水平窗口的寬度；α為裁剪系數(shù)，不大于0.5；L為需要排除數(shù)據(jù)點數(shù)；ai為第i個數(shù)據(jù)點對應(yīng)雷達信號值。

3)排除窗口內(nèi)與Aα極性不同的樣本，具體表達式見公式(8)：

(8)

4)根據(jù)步驟3)中剩余的數(shù)據(jù)值與經(jīng)過α修剪后的數(shù)據(jù)均值接近程度進行加權(quán)，其中越接近的數(shù)據(jù)值將獲得更高的權(quán)重，權(quán)重系數(shù)記為s，見計算公式(9)、(10)：

(9)

(10)

5)對樣本權(quán)重進行歸一化，見公式(11)：

(11)

6)每個樣本乘以相應(yīng)的歸一化權(quán)重并求和，將結(jié)果值作為選取的一維窗口數(shù)據(jù)點的背景噪聲，見公式(12)：

(12)

7)沿水平和垂直方向滑動窗口，并將元素值分配到相應(yīng)的位置，創(chuàng)建與雷達數(shù)據(jù)矩陣具有相同維數(shù)的背景矩陣。

背景矩陣減法計算過程中涉及的參數(shù)有窗口長度、裁剪系數(shù)α、加權(quán)系數(shù)s，以上參數(shù)與具體雷達數(shù)據(jù)相關(guān)，需要根據(jù)要處理的GPR數(shù)據(jù)來決定。但是，經(jīng)驗表明，選擇錯誤的參數(shù)只會增加獲得滿意結(jié)果所需的迭代次數(shù)，并不會影響最終結(jié)果。

2 數(shù)值仿真試驗

為了對比以上三種算法的雜波抑制效果，首先進行了混凝土中鋼筋探測和地下管線探測的數(shù)值仿真試驗，并對仿真數(shù)據(jù)進行處理。

2.1 鋼筋混凝土結(jié)構(gòu)探測

本小節(jié)對探地雷達在鋼筋混凝土結(jié)構(gòu)鋼筋探測進行仿真，建立了仿真模型(圖1)，模型大小為3 m×0.5 m，其中混凝土的相對介電常數(shù)為6，鋼筋埋深為5 cm。離散網(wǎng)格數(shù)為1500×250，網(wǎng)格大小為0.002 m。激勵源采用中心頻率為2 GHz的雷克子波，采樣時窗為8 ns。利用時域有限差分法進行仿真，天線收發(fā)間距為0.05 m，獲得575道信號，道間距為0.005 m。

仿真得到的探地雷達剖面(B-scan)見圖1b，從圖1b中能夠看到較強的直達波信號、地表反射以及鋼筋產(chǎn)生的雙曲線反射。為了對比分析雜波抑制效果，使用減平均道法、奇異分解法、背景矩陣減法三種算法分別對數(shù)據(jù)進行處理，結(jié)果見圖2。背景矩陣減法去除直達波效果最好，處理后的圖像幾乎無虛假反射產(chǎn)生。對于減平均道法和奇異分解法，由于算法存在缺陷，計算雜波信號存在偏差，導(dǎo)致在鋼筋位置產(chǎn)生虛假信號。選取圖2紅色虛線位置繪制單道信號(圖3)，讀取目標(biāo)振幅后計算減平均道法、奇異分解法、背景矩陣減法與原始數(shù)據(jù)中振幅差異分別為-7.2%，-9.2%，8.2%。可以看出三種方法中背景矩陣減法處理后目標(biāo)信號保真度較好，而另兩種算法會對目標(biāo)信號造成衰減。

圖1 鋼筋混凝土仿真模型(a)和鋼筋反射雷達剖面(b)Fig.1 Reinforced concrete simulation model (a) and rebar reflection radar profile (b)

圖2 雜波信號抑制后鋼筋反射雷達剖面Fig.2 Rebar reflection radar profile after clutter suppression(a)減平均道法(average trace subtraction) (b)奇異分解法(singular value decomposition) (c)背景矩陣減法(BMS)

圖3 雜波信號抑制后鋼筋單道反射信號Fig.3 Rebar single channel reflection signal after clutter suppression

2.2 地下管線探測

圖4a為具有分層結(jié)構(gòu)的空心塑料管道數(shù)值模型。模型的大小為3 m×1.8 m，雙層背景介質(zhì)的相對介電常數(shù)由上至下分別為6和10，管的直徑為20 cm。同樣采用時域有限差分(Finite difference time domain， FDTD)法對該模型進行正演數(shù)值模擬，離散網(wǎng)格數(shù)量為600×360，網(wǎng)格大小為0.005 m×0.005 m。激勵源采用中心頻率為900 MHz的雷克子波，采樣時窗為40 ns。天線收發(fā)間距為0.04 m，每0.02 m采集一道雷達數(shù)據(jù)，仿真得到的雷達剖面共有124道數(shù)據(jù)。

圖4 管線仿真模型(a)和反射雷達剖面(b)Fig.4 Pipeline simulation model (a) and pipeline reflection radar profile (b)

仿真得到的GPR剖面圖見圖4b，由圖4b可見較強的直達波信號、地面反射、管線的雙曲線反射和地下介質(zhì)分層處的反射。由于直達波信號較強，雷達剖面圖中地下管線的反射信號較弱，不利于我們對地下目標(biāo)特征分析，需要對數(shù)據(jù)進行去雜波信號處理。為了對比雜波抑制效果，分別使用減平均道法、奇異值分解法和背景矩陣減法三種算法對圖4b數(shù)據(jù)進行處理，處理后的結(jié)果見圖5。圖5表明不同方法均能較好地去除直達波，但奇異分解法和減平均道法對最大振幅處雜波值計算存在偏差，進行雜波抑制時過度去除信號，導(dǎo)致雷達剖面中出現(xiàn)虛假信號(約6 ns和10 ns處)。但背景矩陣減法處理后效果最好，表現(xiàn)為其抑制結(jié)果更加干凈，虛假信號(artefact)能量較低，尤其是6 ns左右的虛假信號。

圖5 雜波信號抑制后管線反射雷達剖面Fig.5 Pipeline reflection radar profile after clutter suppression(a)減平均道法(average trace subtraction) (b)奇異分解法(singular value decomposition) (c)背景矩陣減法(BMS)

為進一步對比不同方法的雜波抑制效果，從圖5中提取出目標(biāo)反射振幅，并通過計算處理前后振幅變化百分比來評估不同算法的優(yōu)劣。選取圖5紅色虛線位置處目標(biāo)反射單道信號，繪制目標(biāo)反射信號的A-scan對比圖(圖6)。從圖6中讀取出目標(biāo)反射振幅后發(fā)現(xiàn)三種算法處理后導(dǎo)致目標(biāo)信號有所擾動，但對于背景矩陣減法處理后的目標(biāo)信號振幅的衰減率0.5%遠小于其他兩種算法的衰減率5.4%。

圖6 雜波抑制后管線單道反射信號Fig.6 Single channel reflectedion signal of pipeline after clutter suppression

3 現(xiàn)場試驗

本節(jié)進一步對比分析了三種雜波去除方法在探地雷達實地數(shù)據(jù)處理中的效果。

3.1 鋼筋混凝土結(jié)構(gòu)探測

本小節(jié)對混凝土中鋼筋探測數(shù)據(jù)進行了處理，數(shù)據(jù)采集使用意大利IDS公司生產(chǎn)的TRHF型號探地雷達，天線中心頻率為2 GHz。測線水平距離約為0.7 m，共采集353道數(shù)據(jù)。采樣時窗長度為8 ns，共512個采樣點。

圖7為混凝土中鋼筋的雷達反射剖面(B-scan)，可以清晰地看到兩根鋼筋的雙曲線反射信號以及地面雜波信號。采用前文介紹的三種雜波抑制算法進行處理，結(jié)果見圖8。從圖8中可以看到三種算法對雜波處理的效果相似，均能將直達波很好的去除，可以更加清晰地識別鋼筋的雙曲線反射。然而減平均道法(圖8a)和奇異分解法(圖8b)在處理后雙曲線頂點附近產(chǎn)生虛假信號，這是因為算法本身的缺陷導(dǎo)致的。同樣的，提取圖8紅色虛線位置的單道A-Scan，結(jié)果見圖9。由圖9可見，三種算法對直達波信號能夠很好的去除，讀取目標(biāo)振幅后發(fā)現(xiàn)均存在一定衰減。但相較于其他兩種算法，背景矩陣減法處理后對目標(biāo)信號衰減值為原信號的31%，而減平均道和奇異分解法為原信號的39%，所以背景矩陣減法仍優(yōu)于其他兩種方法。

圖7 鋼筋反射雷達剖面Fig.7 Radar profile of rebar reflection

圖9 雜波信號抑制后鋼筋單道反射信號Fig.9 Single channel reflection signal of rebar after clutter suppression

3.2 地下管線探測

本小節(jié)對某地下管線的雷達探測數(shù)據(jù)進行了處理，試驗采用意大利IDS探地雷達，選用中心頻率400 MHz天線。測線長度約2.2 m，共采集92道數(shù)據(jù)，采樣時窗長度為50 ns，采樣點數(shù)為512個。

圖10為實際采集的地下管線探地雷達剖面。由于有較強的直達波及不均勻背景介質(zhì)引起的雜波，地下管線產(chǎn)生的雙曲線型反射信號不能清楚呈現(xiàn)。采用前文介紹的三種雜波抑制算法分別對實測數(shù)據(jù)進行處理，處理后結(jié)果見圖11。由圖11可見三種算法對雜波處理的效果接近。直達波被很好地去除，直流成分也在一定程度上得到了去除，從而使管線的雙曲線反射變得明顯。但減平均道法(圖11a)和奇異分解法(圖11b)處理后的結(jié)果在目標(biāo)雙曲線頂點位置處存在虛假反射，這與前文仿真數(shù)據(jù)雜波抑制后的結(jié)果相似，而背景矩陣減法(圖11c)沒有出現(xiàn)虛假反射。同樣繪制圖11紅色虛線位置處不同算法處理后管線目標(biāo)單道信號，見圖12，讀取目標(biāo)振幅后，計算出減平均道法、奇異分解法、背景矩陣減法三種算法處理后相對原始數(shù)據(jù)振幅百分比分別為120%，123%，139%。結(jié)果表明背景矩陣減法優(yōu)于其余兩種方法。

圖10 管線反射雷達剖面Fig.10 Pipeline reflection radar profil

圖11 雜波信號抑制后管線反射雷達剖面Fig.11 Radar profile of pipeline reflection after clutter suppression(a)減平均道法(average trace subtraction) (b)奇異分解法(singular value decomposition) (c)背景矩陣減法(BMS)

圖12 雜波信號抑制后管線單道反射信號Fig.12 Single channel reflection signal of pipeline after clutter suppression

4 結(jié)論

本文提出使用背景矩陣減法對探地雷達數(shù)據(jù)中雜波信號進行抑制，并對比研究了減平均道法、奇異分解法和背景矩陣減法三種算法在探地雷達仿真和實測數(shù)據(jù)處理效果。結(jié)果表明：

1)減平均道法雖然算法簡單，但容易受強目標(biāo)信號影響，導(dǎo)致計算得到噪聲信號部分由背景噪聲組成，部分由目標(biāo)的特征組成，處理后的效果不理想，甚至產(chǎn)生較明顯的虛假信號。

2)奇異分解法需要選擇適當(dāng)?shù)钠娈愔祦磉_到抑制雜波目的，但奇異值的選取往往不唯一，不同雷達數(shù)據(jù)可能不同，適應(yīng)性較差。與減平均道法類似，由于算法存在缺陷，處理后的雷達剖面中易產(chǎn)生虛假反射，不能達到令人滿意的雜波抑制效果。

3)背景矩陣減法是用較復(fù)雜過程創(chuàng)建一個完整的雜波噪聲矩陣，因此它是由純雜波組成的，而且選擇錯誤參數(shù)只會影響迭代計算次數(shù)，相比于減平均道法和奇異分解法其在探地雷達數(shù)據(jù)雜波信號抑制中有更好的穩(wěn)健性和適應(yīng)性。