基于主成份分析的探地雷達(dá)信號(hào)去噪研究

劉 雷, 譚方玉

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

0 前言

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

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

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

1 主成分分析(PCA)

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

(1)

其中

(2)

(3)

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

Sxuk=λkuk

(4)

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

(5)

因此xi可以近似為：

(6)

PCA的算法步驟：

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

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

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

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

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

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

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

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

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

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

(7)

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

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

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

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

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

(8)

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

(9)

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

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

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

建立圖1所示的層狀模型。模擬區(qū)域?yàn)?.0 m×1.0 m，模型相對(duì)介電常數(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布萊克曼-哈里斯脈沖置于地表，模擬時(shí)間間隔為0.01 ns，模擬時(shí)窗為25 ns，采樣時(shí)間間隔為0.04 ns。我們采用gprMax軟件[17]對(duì)該模型進(jìn)行正演數(shù)值模擬，該GPR剖面共181道，道間距為0.005 m，收發(fā)距為0.1 m。

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

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

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

圖3 不同域去噪處理后均方誤差(RMSE)和信噪比(SNR)對(duì)比圖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)，兩種不同域去噪處理均存在一個(gè)合適的K值使得均方根誤差最小且信噪比最大。選擇低于最佳值的K會(huì)使得信號(hào)失真，而選擇高于最佳值的K會(huì)使得恢復(fù)信號(hào)中含有更多的噪音。由于合成噪聲具有隨機(jī)頻率成分，在時(shí)間和頻率范圍內(nèi)的去噪效果都相似的。時(shí)域和頻域中均方根誤差與信噪比與K值的關(guān)系滿足相同的規(guī)律，RMSE分別在K=5和K=4取得最小值。

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

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

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

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

圖5 RMSE最小值點(diǎn)與L曲線拐點(diǎn)對(duì)應(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) 時(shí)間域K=5時(shí)PCA去噪結(jié)果；(b) 時(shí)間域K=8(L曲線法)時(shí)PCA去噪結(jié)果；(c)頻率域K=4時(shí)PCA去噪結(jié)果；(d) 頻率域K=7(L曲線法)時(shí)PCA去噪結(jié)果

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

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

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

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

圖7(a)為探測(cè)某隧道地下管線的實(shí)測(cè)探地雷達(dá)剖面，檢測(cè)過(guò)程中采用美國(guó)GSSI公司生產(chǎn)的SIR-3000型探地雷達(dá)儀進(jìn)行檢測(cè)，按照探測(cè)目的及要求，選用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)時(shí)間域K=5時(shí)PCA去噪結(jié)果；(d)時(shí)間域K=8(L曲線法)時(shí)PCA去噪結(jié)果；(e)頻率域K=4時(shí)PCA去噪結(jié)果；(f)頻率域K=7(L曲線法)時(shí)PCA去噪結(jié)果

圖7 實(shí)測(cè)探地雷達(dá)剖面Fig.7 Measured GPR profile of actual data(a)實(shí)測(cè)探地雷達(dá)信號(hào)剖面; (b)PCA去噪剖面

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

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

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

5 結(jié)論

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

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

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