劉 雷, 譚方玉

(1.重慶市勘測院，重慶 401121； 2.貴州省有色金屬和核工業(yè)地質(zhì)勘查局 物化探總隊，都勻 558000)

0 前言

探地雷達(Ground penetrating radar, GPR)是利用高頻寬帶電磁(1Mhz-10GHz)波，對地下結(jié)構(gòu)和特性或者物體內(nèi)部不可見目標(biāo)體進行探測定位的一種淺地表地球物理方法。憑借其具有操作簡單、分辨率高、抗干擾能力強、適應(yīng)性強、解釋成果準(zhǔn)確及效率高等優(yōu)點，在交通建設(shè)、地質(zhì)勘探、水利工程、城市建設(shè)、地質(zhì)災(zāi)害、農(nóng)業(yè)地質(zhì)探測、考古、隧道、軍事等許多領(lǐng)域都表現(xiàn)出廣闊的應(yīng)用前景[1]。在實際探測過程中，淺部地下介質(zhì)復(fù)雜多變，具有頻散特性，電磁波在其中傳播衰減速度快，散射非常強烈，且存在各種人為設(shè)施的雜亂回波以及儀器本身或系統(tǒng)噪聲等多種因素的影響，使得得到的GPR剖面中含有各種各樣的雜波和噪聲，降低了雷達信號的質(zhì)量。因此，研究高信噪比的探地雷達數(shù)據(jù)去噪方法顯得尤為重要[2]。實際噪聲在實踐中是不可避免的，必須在預(yù)處理數(shù)據(jù)的同時對噪聲處理。對GPR數(shù)據(jù)進行去噪已經(jīng)得到了廣泛地研究，這些方法涉及將數(shù)據(jù)轉(zhuǎn)換為不同的域，然后利用波場在域內(nèi)的特性進行去噪(如傅立葉變換、加窗傅立葉變換[3]、小波變換[4]、小波-KL變換[5]、S變換[6]、希爾伯特-黃變換[7]、經(jīng)驗?zāi)B(tài)分解[8]以及獨立成份分析[9]等)。這些方法也存在著各自難以克服的缺點。以小波變換為例，近年來小波變換被廣泛應(yīng)用于信號處理領(lǐng)域，其中以小波閾值去噪應(yīng)用最為廣泛，雖然這種方法具有多分辨率的特點，但其不足的地方在于具體小波閾值函數(shù)、軟硬閾值、尺度系數(shù)的選擇對去噪結(jié)果影響較大，同時對于非連續(xù)信號，采用小波閾值去噪后容易產(chǎn)生偽吉布斯現(xiàn)象[8]，使得小波變換在探地雷達數(shù)據(jù)處理的應(yīng)用和推廣中受到了限制。因此，有必要引入一種快速便捷的方法來去除探地雷達信號中的噪聲干擾，以提高雷達數(shù)據(jù)的信噪比。

地下介質(zhì)通常具有一定規(guī)律，實測的GPR數(shù)據(jù)中由異常體或者界面所引起的數(shù)據(jù)信號一般具有相關(guān)性，而有噪聲所對應(yīng)的成分無明顯相關(guān)性[10]。因此可以找到一個域使得其相關(guān)性降序排列，達到信噪分離的目的。主成分分析[11](Principal component analysis, PCA)是一種常用的多元統(tǒng)計數(shù)據(jù)分析方法，該方法能夠?qū)⒃紨?shù)據(jù)中相關(guān)性較強的多個變量轉(zhuǎn)化成含有原有變量最大信息的彼此獨立的少數(shù)幾個變量。PCA 方法可以將數(shù)據(jù)中最“主要”的成分和結(jié)構(gòu)找出來，去除噪聲和冗余成分，使原始復(fù)雜數(shù)據(jù)的維數(shù)降低，揭示出在復(fù)雜數(shù)據(jù)下所隱藏的簡單結(jié)構(gòu)，采用這些主要成分重構(gòu)信號，可以達到去噪的效果。

筆者將主成分分析用于探地雷達信號的去噪處理中，為探地雷達信號去噪提供一種快速便捷的方法，以提高實際生產(chǎn)中GPR數(shù)據(jù)的信噪比。采用PCA確定二維剖面數(shù)據(jù)中最大變化的模式，將PCA分析之后的特征值降序排列，采用L曲線信噪分離準(zhǔn)則選取K值，并采用前K個主成份將信號重構(gòu)，在保留原始信號的同時消除大部分噪聲。

1 主成分分析(PCA)

主成分分析是一種統(tǒng)計方法，通過正交變換將一組可能存在相關(guān)性的變量轉(zhuǎn)換為一組線性不相關(guān)的變量，轉(zhuǎn)換后的這組變量叫主成分，其常用于降維和特征提[12-14]。

(1)

其中

(2)

(3)

令Sx為(xi)的協(xié)方差矩陣。因此tr(Sy) =tr(ASxAT)通過使用拉格朗日乘數(shù)并取導(dǎo)數(shù)，可以得到

Sxuk=λkuk

(4)

這意味著uk是Sx的特征向量。此時xi可以表示為：

(5)

因此xi可以近似為：

(6)

PCA的算法步驟：

1) 將原始數(shù)據(jù)按列組成N行M列矩陣x。

2) 將x的每一行進行零均值化，即減去這一行的均值。

3) 求出協(xié)方差矩陣Sx。

4) 求出協(xié)方差矩陣的特征值λk及對應(yīng)的特征向量uk。

5) 將特征向量按對應(yīng)特征值大小從上到下按行排列成矩陣，取前K行組成矩陣A。

6)y=Ax即為降維到K維后的數(shù)據(jù)。

2 GRP數(shù)據(jù)PCA去噪流程

實際測量中，噪聲的類型是多種多樣的，除開地層中不同介質(zhì)產(chǎn)生的背景干擾噪聲，還有儀器本底噪聲，直流漂移的干擾和各種外部因素干擾[15]。從野外采集到的探地雷達原始數(shù)據(jù)，其中既包含有用信息，也包含各種干擾噪聲，有時有用信息會被噪聲所掩蓋，此時需要經(jīng)過數(shù)據(jù)處理，才能得到有助于解釋的數(shù)據(jù)和圖像。需要壓制噪聲，增強信號，提高資料信噪比，以便從數(shù)據(jù)中提取速度、振幅、頻率、相位等特征信息，幫助解釋人員對資料進行地質(zhì)解釋。

地下介質(zhì)中會包含多種的天然或人為的雜質(zhì)，導(dǎo)致接收的回波中包含很多干擾信號，這些干擾統(tǒng)稱為背景隨機噪聲。探地雷達勘探中,記錄雷達波時, 為了保持更多的波的特征, 通常采用寬頻帶進行記錄。因此, 在寬頻帶范圍內(nèi)記錄了各種反射波的同時，也記錄了各種噪聲。假設(shè)探地雷達信號是長度為N的時間域離散信號s(n)，被噪聲e(n)污染，所測得的含噪數(shù)據(jù)表示為[16]式(7)。

f(n)=s(n)+e(n)

(7)

去噪的主要任務(wù)是盡可能地將實際信號與噪聲信號分離開。通常出現(xiàn)的噪聲都服從或近似服從高斯分布，因此假設(shè)e(n)為高斯噪聲。一般利用數(shù)學(xué)變換解決信號去噪問題，才能夠時域轉(zhuǎn)換到頻域加以解決，如傅里葉數(shù)字濾波和小波變換。筆者研究了PCA在時間域和頻率域下去除噪聲的有效性：在時間域直接將數(shù)據(jù)進行PCA去噪；對于頻率域而言，首先將每道單道數(shù)據(jù)進行的離散傅立葉變換，對非負頻率數(shù)據(jù)進行PCA去噪，負頻率數(shù)據(jù)可以通過恢復(fù)的非負頻率片的復(fù)共軛得到，然后使用逆離散傅立葉變換將完整頻率數(shù)據(jù)變換到時域與原始信號進行比較。時間域與頻率域的實驗流程如下。

1)時域處理流程:①在模擬數(shù)據(jù)上加入高斯白噪聲;②在時間域數(shù)據(jù)上進行PCA去噪;③對降噪質(zhì)量進行評估對比。

2)頻域處理流程:①在模擬數(shù)據(jù)上加入高斯白噪聲;②對時間維度進行離散傅里葉變換;③在非負頻率數(shù)據(jù)上執(zhí)行PCA;④對時間維度進行逆離散傅里葉變換;⑤對降噪質(zhì)量進行評估對比。

將二維時間域和頻率域數(shù)據(jù)進行PCA去噪，其特征值表示為每個相應(yīng)的特征向量的方差，使用前K個特征向量重構(gòu)數(shù)據(jù)，重構(gòu)數(shù)據(jù)中只包含最大變化模式的數(shù)據(jù)，達到保留原始信號的同時消除大部分噪聲。而不同K值的選取會對去噪效果產(chǎn)生影響。為了計算最佳模式數(shù)K，筆者采用均方根誤差(RMSE)以及信噪比(SNR)兩個性能指標(biāo)對降噪質(zhì)量進行評估對比。

均方根誤差(RMSE)：去噪信號與原信號的均方誤差為式(8)。

(8)

信噪比(SNR)：原始信號能量與噪聲能量的比值為式(9)。

(9)

信噪比越高，則去噪效果越好。

3 合成GRP數(shù)據(jù)PCA去噪實驗

為了研究PCA去噪效果，分析最佳K值的取值方法，首先對合成的理論GPR記錄進行處理。此處建立了兩個數(shù)值模擬模型數(shù)據(jù)。

建立圖1所示的層狀模型。模擬區(qū)域為1.0 m×1.0 m，模型相對介電常數(shù)從上至下依次為3、5、8、5，電導(dǎo)率依次為3 mS/m、5 mS/m、8 mS/m、5 mS/m。模型離散網(wǎng)格為201×201，網(wǎng)格間距為0.005 m，CPML吸收層為10層。將主頻為900 MHz布萊克曼-哈里斯脈沖置于地表，模擬時間間隔為0.01 ns，模擬時窗為25 ns，采樣時間間隔為0.04 ns。我們采用gprMax軟件[17]對該模型進行正演數(shù)值模擬，該GPR剖面共181道，道間距為0.005 m，收發(fā)距為0.1 m。

圖2顯示了模擬的GPR剖面和加入高斯白噪聲的GPR剖面。采用PCA去噪方法，分別在時間域和頻率域?qū)υ肼晹?shù)據(jù)進行去噪處理，為了計算PCA去噪中的最佳數(shù)量K，在去噪過程中保留了重構(gòu)信號和原始信號之間的均方誤差(RMSE)和信噪比(SNR)用于衡量恢復(fù)信號的準(zhǔn)確性(圖3)。

圖1 層狀模型Fig.1 Layered model(a)相對介電常數(shù);(b)導(dǎo)電率

圖2 高斯白噪聲GPR剖面圖Fig.2 GPR profile of white gaussian noise(a)層狀模型的原始剖面圖;(b)含噪剖面

圖3 不同域去噪處理后均方誤差(RMSE)和信噪比(SNR)對比圖Fig.3 Comparison diagram of mean square error (RMSE) and signal to noise ratio (SNR) after denoising in different domains(a)RMSE隨 k值的變化圖；(b)SNR隨 k值的變化圖

由圖3可以發(fā)現(xiàn)，兩種不同域去噪處理均存在一個合適的K值使得均方根誤差最小且信噪比最大。選擇低于最佳值的K會使得信號失真，而選擇高于最佳值的K會使得恢復(fù)信號中含有更多的噪音。由于合成噪聲具有隨機頻率成分，在時間和頻率范圍內(nèi)的去噪效果都相似的。時域和頻域中均方根誤差與信噪比與K值的關(guān)系滿足相同的規(guī)律，RMSE分別在K=5和K=4取得最小值。

主成分分析去噪可以看做是將觀測到的信號劃分成信號子空間和噪聲子空間兩部分。信號子空間是其中的主要成分，可以由前n個最大的特征值對應(yīng)的特征向量計算得到。相應(yīng)的，其余的特征值對應(yīng)的是噪聲子空間，也可以由特征值對應(yīng)的特征向量計算得到。這樣得到的信號空間中可以看成是不含噪聲的。主成分分析的關(guān)鍵問題是如何確定需要保留的主成分的個數(shù)K。在主成分分析中，確定合理的K值是一個比較難解決的問題，其將信號分成了兩個子空間，一般來說，信號的信噪比越低，合理的K值越難確定。

然而對于實際數(shù)據(jù)而言，我們并不能通過均方根誤差曲線以及信噪比曲線來選取最佳的K值。通常在主成分分析的過程中，選取主要成分時使用的方法是計算各成分貢獻率以及前個成分的累計貢獻率[13]。只要大致規(guī)定一個累計貢獻率的百分比(一般為 85%以上即可)便可以決定選擇幾個主成分。

這里我們采用L曲線法[18]來進行最佳K值的選取。在基于L曲線的信噪分離準(zhǔn)則[19]中，對特征值進行降序排列，以特征值索引作為橫坐標(biāo)，特征值作為縱坐標(biāo)。那么所獲得的特征值曲線可以近似為L曲線。一般來講，含噪數(shù)據(jù)的L曲線可明顯分為兩個部分：①對應(yīng)有效信號成分的特征值幅度較大，衰減較快；②而對應(yīng)噪聲成分的特征值幅度要小得多，衰減較慢且數(shù)值變化比較平穩(wěn)。根據(jù)這個特性，我們可以大致確定 L曲線的拐點位置。將拐點位置對應(yīng)的K值作為保留成分的依據(jù)，該值表示作為PCA分析中的信噪分離的基準(zhǔn)點。

圖4 不同域特征值的變化曲線對比圖Fig.4 Comparison graph of change curves of eigenvalues in different domains(a) 時間域特征值變化曲線；(b)頻率域特征值變化曲線

圖5 RMSE最小值點與L曲線拐點對應(yīng)K值的去噪結(jié)果Fig.5 Denoising results of the minimum value point of RMSE and the corresponding value of K at the inflection point of L curve(a) 時間域K=5時PCA去噪結(jié)果；(b) 時間域K=8(L曲線法)時PCA去噪結(jié)果；(c)頻率域K=4時PCA去噪結(jié)果；(d) 頻率域K=7(L曲線法)時PCA去噪結(jié)果

圖4為時間域和頻率域中特征值的變化曲線，其拐點分別為K=8和K=7。圖5為RMSE最小值點以及L曲線拐點對應(yīng)K值的去噪結(jié)果，由圖5可以發(fā)現(xiàn)，L曲線拐點處K值的去噪效果與RMSE最小點K值的去噪效果基本相當(dāng)，僅在細微地方有些差別，均能在保留原始信號的基礎(chǔ)上去除大部分噪聲。

為了進一觀察基于 L 曲線的信噪分離準(zhǔn)則的PCA去噪效果，圖6顯示了第90道的原始記錄、噪聲記錄以及不同域下不同K值的PCA去噪結(jié)果?？梢园l(fā)現(xiàn)RMSE對應(yīng)的K值在3個反射波處(10 ns ～15 ns)出現(xiàn)了一些波形失真，而基于L曲線信噪分離準(zhǔn)則選取K值PCA去噪曲線此處波形與原始信號基本相同

從上述分析可以看到，基于L曲線的信噪分離準(zhǔn)則PCA去噪方法可以在保持有用數(shù)據(jù)不丟失，原幅值不失真的基礎(chǔ)上，使噪聲得到了較好地抑制，提高數(shù)據(jù)的信噪比。

4 實測數(shù)據(jù)處理

圖7(a)為探測某隧道地下管線的實測探地雷達剖面，檢測過程中采用美國GSSI公司生產(chǎn)的SIR-3000型探地雷達儀進行檢測，按照探測目的及要求，選用900 MHz天線。

圖6 第90道的原始記錄、噪聲記錄以及不同域下不同K值的PCA去噪結(jié)果Fig.6 The 90th data of original records and noise records and PCA denoising results with different K values in different domains(a)原始單道記錄；(b)含噪單道記錄；(c)時間域K=5時PCA去噪結(jié)果；(d)時間域K=8(L曲線法)時PCA去噪結(jié)果；(e)頻率域K=4時PCA去噪結(jié)果；(f)頻率域K=7(L曲線法)時PCA去噪結(jié)果

圖7 實測探地雷達剖面Fig.7 Measured GPR profile of actual data(a)實測探地雷達信號剖面; (b)PCA去噪剖面

圖8 實測數(shù)據(jù)PCA變換特征值變化曲線Fig.8 Actual data PCA eigenvalues value change curve

圖9 實測數(shù)據(jù)第100道數(shù)據(jù)與PCA其消噪圖Fig.9 The 100th data before and after PCA denoising of actual data

共采集400道數(shù)據(jù)，每道512個采樣點，剖面水平距離為4 m，時窗長度為20 ns。采用之前討論的L曲線信噪分離準(zhǔn)則的PCA去噪方法在時間域?qū)υ搶崪y數(shù)據(jù)進行去噪處理。圖8為原始數(shù)據(jù)PCA的歸一化特征值曲線，根據(jù)L曲線信噪分離原則，其拐點位于K=27的位置。圖7(b)為取前27個主成份重構(gòu)得到的PCA去噪數(shù)據(jù)剖面圖，圖中去除了深部(采樣點250～500)的雪花狀的隨機噪聲并使得有效信息得到保留，深部數(shù)據(jù)剖面圖變得清晰，特別是位于第125道400采樣點的位置的雙曲線反射異常以及280道250采樣點位置的弱雙曲線反射異常，進一步證了該方法可行性。圖9為第125道的單道數(shù)據(jù)去噪前后的對比，從圖9可以看出，該方法在有效抑制噪聲的同時，能大幅地保留原始信號信息。本文方法能夠很好地去除高頻噪聲，并且可以很好應(yīng)用于實際生產(chǎn)中。

5 結(jié)論

筆者主要研究了主成份分析算法在探地雷達信號去噪中地應(yīng)用，對主成份分析方法基本理論和PCA算法流程進行了介紹，并應(yīng)用PCA算法對正演模擬的GPR剖面數(shù)據(jù)和實測數(shù)據(jù)作為實驗數(shù)據(jù)，對主成分分析法的去噪性能進行了測試與驗證。結(jié)果表明：

1) PCA算法對含有高斯白噪聲探地雷達數(shù)據(jù)的去噪處理有很好的適應(yīng)性和穩(wěn)健性。在時間和頻率范圍內(nèi)的去噪效果基本相同，兩個域內(nèi)的PCA去噪過程中的均方根誤差與信噪比與K值的關(guān)系滿足相同的規(guī)律，均能較好地對GPR數(shù)據(jù)進行噪聲抑制，達到提高數(shù)據(jù)信噪比目的，有助于突出探地雷達剖面中異常體特征。

2)使用基于L曲線的信噪分離準(zhǔn)則可以較為準(zhǔn)確的將有效信號成分提取出來，基于L曲線的信噪分離準(zhǔn)則的PCA去噪方法可以直接應(yīng)用與GPR實測數(shù)據(jù)的去噪中，能在保留原始信號信息得的基礎(chǔ)上有效地抑制噪聲，為主成分分析去噪中的最佳重構(gòu)成分K值的選取提取了新思路。