基于MUSIC譜估計(jì)的探地雷達(dá)回波處理方法

(長(zhǎng)安大學(xué) 信息工程學(xué)院，陜西 西安 710064)

探地雷達(dá)(Ground Penetrating Radar，GPR，也稱地質(zhì)雷達(dá))系統(tǒng)是一種通過(guò)頻段位于106～109Hz的無(wú)線電波來(lái)探測(cè)地下結(jié)構(gòu)的無(wú)損探測(cè)儀器。探地雷達(dá)可以探測(cè)金屬及非金屬埋藏物，現(xiàn)在已經(jīng)廣泛應(yīng)用于考古、礦產(chǎn)資源勘探、災(zāi)害地質(zhì)勘察、路面結(jié)構(gòu)檢測(cè)等眾多領(lǐng)域。在實(shí)際應(yīng)用中，探地雷達(dá)接收到的波形常常會(huì)產(chǎn)生畸變且信號(hào)的信噪比較弱[1]。而由于被測(cè)物體具有強(qiáng)反射性、被測(cè)媒質(zhì)內(nèi)的大尺寸顆粒反射及不均勻的物體電導(dǎo)率、雷達(dá)設(shè)備天線間的耦合等因素導(dǎo)致探地雷達(dá)接收到的信號(hào)中含有大量的噪聲、直達(dá)波等各種散射體的雜波分量。這些信號(hào)由于傳播路徑短且衰減少，因此在到達(dá)時(shí)間上相互重疊，導(dǎo)致實(shí)際所測(cè)的目標(biāo)信號(hào)常常被淹沒(méi)，因此需要對(duì)采集到的探地雷達(dá)回波信號(hào)進(jìn)行濾波處理，提高所測(cè)數(shù)據(jù)的精度[2-5]。

探地雷達(dá)回波信號(hào)的處理通常從時(shí)域和頻域兩個(gè)方面進(jìn)行：① 時(shí)域處理方法，對(duì)接收到的信號(hào)直接提取并分離出有效的速度、振幅等信號(hào)值,例如數(shù)字濾器的處理方法；② 頻域處理方法，對(duì)接收到的雷達(dá)回波信號(hào)進(jìn)行頻譜分析。經(jīng)典常用的算法有小波分析法、FFT算法等。其中，F(xiàn)FT算法對(duì)數(shù)據(jù)樣本的長(zhǎng)度要求較高，而小波分析在特定情況下(已知噪聲的頻率范圍且信號(hào)和噪聲的頻帶相互分離時(shí))非常有效。但在實(shí)際應(yīng)用中對(duì)廣泛存在的白噪聲處理效果較差[5-8]。因此，筆者在雷達(dá)信號(hào)預(yù)處理后引入譜分析中信號(hào)子空間與噪聲子空間的概念，將MUSIC算法與最小二乘法相結(jié)合來(lái)進(jìn)行探地雷達(dá)雜波及噪聲處理。

1 探地雷達(dá)數(shù)據(jù)處理方法

1.1 直流偏移量處理

在探地雷達(dá)系統(tǒng)中，有時(shí)雷達(dá)剖面上的數(shù)據(jù)含有直流漂移量，這時(shí)的數(shù)據(jù)全是正的或全是負(fù)的，亦或正負(fù)半周出現(xiàn)不對(duì)稱情況。此時(shí)需要對(duì)數(shù)據(jù)進(jìn)行處理，若直接利用接收到的回波信號(hào)，肯定會(huì)產(chǎn)生誤差，這對(duì)后續(xù)波形的利用產(chǎn)生影響，因此需要消除或壓制[9]直流成分。通常直流波去除采用消除水平同相軸的方法，該方法假設(shè)直流偏移量是一個(gè)常量，因此在探地雷達(dá)測(cè)量過(guò)程中,將每道回波數(shù)據(jù)求平均值即可得到直流偏移量，再將每條原始回波數(shù)據(jù)依次減去直流偏移量即可得到移除直流偏移量的回波數(shù)據(jù)。

假設(shè)回波信號(hào)為X(n)，消除水平同相軸去除直流偏移的方法可以表示為

(1)

式中，X(n)為未經(jīng)處理的原始信號(hào);n為采樣點(diǎn)數(shù)量;X′(n)為去除直流偏移量后得到的回波信號(hào)。

1.2 均值法濾波

均值濾波是一種線性處理技術(shù)，是一種最簡(jiǎn)單且最常用的數(shù)字濾波方法，對(duì)抑制高斯噪聲有較好的效果。該方法是在某時(shí)刻連續(xù)多次進(jìn)行信號(hào)采樣，對(duì)得到的采樣值進(jìn)行算術(shù)平均處理，得到的值將作為此刻的最終信號(hào)值，其中連續(xù)采樣的次數(shù)視具體情況而定。

假設(shè)在測(cè)量過(guò)程中探地雷達(dá)獲得N段采樣數(shù)據(jù)，在此記為xi(i=1,2,…,n)。根據(jù)最小二乘法的原理可知：存在一個(gè)數(shù)值y，使得各采樣值間的偏差平方和為最小值，則y就是最接近真實(shí)值的數(shù)據(jù)，則有

(2)

依據(jù)一元函數(shù)求極值的原理可以得到

(3)

在隨機(jī)噪聲干擾的環(huán)境下，每次采樣接收的信號(hào)均可以表示為

xi(n)=si(n)+ei(n)

(4)

式中，si(n)為目標(biāo)的回波信號(hào)；ei(n)為噪聲。

探地雷達(dá)對(duì)同一目標(biāo)多次采樣時(shí)，N次采樣值的算術(shù)平均值為

sk(n)+ek(n)]

(5)

當(dāng)采樣時(shí)間很短時(shí)，采樣信號(hào)可被看作近似相等，因此在N次采樣后，得到的信號(hào)s(n)的分量可表示為

(6)

(7)

假定采樣過(guò)程中，噪聲環(huán)境不發(fā)生變化，則噪聲ei(n)為相互獨(dú)立且概率密度相同的隨機(jī)變量。這樣根據(jù)噪聲ei(n)的性質(zhì)，經(jīng)過(guò)平均值濾波后的噪聲大小為

(8)

1.3 MUSIC功率譜估計(jì)

探地雷達(dá)接收到的信號(hào)包括被測(cè)物體回波、噪聲和強(qiáng)雜波信號(hào)。通常雷達(dá)回波可以用時(shí)間的實(shí)數(shù)振蕩函數(shù)來(lái)模擬，而雷達(dá)系統(tǒng)中的所接收回波中的強(qiáng)雜波可以近似看成服從高斯分布的有色噪聲，在本文中將雜波信號(hào)假定為經(jīng)過(guò)濾波后的高斯白噪聲[10]。因此，雷達(dá)回波信號(hào)X(n)可看作K個(gè)復(fù)正弦信號(hào)與高斯白噪聲信號(hào)的組合，即

(9)

此時(shí)，探地雷達(dá)信號(hào)的濾波問(wèn)題變?yōu)橛梢褱y(cè)得到的回波信號(hào)X(n)恢復(fù)出回波信號(hào)中的未知量Ai、ωi。下面采用多重信號(hào)分類的MUSIC算法來(lái)進(jìn)行ωi估算。

MUSIC(Multiple Signal Classification)譜估計(jì)法，又稱為多重信號(hào)分類算法。該算法是Schmidt于1979年提出的一種高分辨率空間譜頻率估計(jì)方法[11]。此算法依據(jù)矩陣特征值分解理論，將信號(hào)分為信號(hào)子空間和噪聲子空間，利用信號(hào)與噪聲的正交性，從自相關(guān)矩陣求出最小的特征值及相對(duì)應(yīng)的特征向量并構(gòu)建空間譜函數(shù)，再通過(guò)譜峰搜索估計(jì)出未知的頻率。

依據(jù)MUSIC算法的工作原理，探地雷達(dá)信號(hào)X(n)的自相關(guān)函數(shù)可表示為

(10)

假設(shè)接收到的雷達(dá)信號(hào)可表示為N個(gè)采樣值的組合，則信號(hào)矢量可表示為

ei=[1,ejωi,…，ej(N-1)ωi]T

(11)

此時(shí)式(10)就可以寫成如下形式：

(12)

式中，I為N×N的單位矩陣。

Rx=Rs+Rw

(13)

將矩陣Rs特征分解得

(14)

式中，λi為特征向量υi所對(duì)應(yīng)的特征值。由于特征向量υi之間是正交的，矩陣Rs的秩為M,所以其特征值中必有N-M個(gè)為0。此時(shí)式(14)可以寫成如下形式:

(15)

對(duì)于噪聲相關(guān)矩陣也可以表示為

(16)

式(10)可以表示為

(17)

(18)

實(shí)際做法是構(gòu)造函數(shù)，構(gòu)造的函數(shù)為

(19)

設(shè)當(dāng)ω=ωi時(shí)功率譜為無(wú)窮大，此時(shí)相關(guān)矩陣估計(jì)值為ωi，為PMUSIC(ω)取到最大值時(shí)所對(duì)應(yīng)的值。由于相關(guān)矩陣是估計(jì)值，必然存在誤差，因此PMUSIC(ωi)為有限值，但會(huì)呈現(xiàn)出很尖的峰值，此時(shí)峰值對(duì)應(yīng)的頻率就是復(fù)正弦信號(hào)的頻率。

1.4 計(jì)算幅值與相角

信號(hào)的最佳頻率ωi和正弦信號(hào)個(gè)數(shù)M可以由MUSIC譜方法準(zhǔn)確地估算出來(lái)，但是通過(guò)該方法無(wú)法檢測(cè)出正弦波的幅值和相角，因此需要結(jié)合其他方法進(jìn)行信號(hào)參數(shù)計(jì)算。主要采用求解過(guò)程簡(jiǎn)單、計(jì)算量較小的最小二乘法[12]來(lái)估算振幅值和相角。

假設(shè)已知信號(hào)記為X(t)，振幅和初相分別為Ai和φi，最佳組成頻率成分為fi，其中fi=2π/ωi,噪聲為w(t)。此時(shí)信號(hào)的數(shù)學(xué)關(guān)系表達(dá)式為

(20)

由式(20)可知，雷達(dá)信號(hào)經(jīng)過(guò)時(shí)域n點(diǎn)采樣后可以表達(dá)為

(21)

通常情況下，當(dāng)n>>2m時(shí)，式(19)為可表示為一組線性超定方程組，將此方程寫為矩陣形式：

A=

當(dāng)噪聲較小時(shí)，式(21)可寫成如下形式：

B=AX

(22)

取n=2m,矩陣A即為常數(shù)陣，將式(22)兩邊同乘以A-1,則可以求解出X，即

X=A-1B

(23)

但在實(shí)際應(yīng)用中常取n>2m,此時(shí)矩陣A不為方陣不存在逆矩陣，可以用A的偽逆矩陣A+來(lái)求得未知數(shù)X的解，式(22)的解可表示成如下形式：

X=A+B

(24)

求出X中的所有元素后，可以用式(25)和式(26)求出振幅和相位：

(25)

(26)

2 仿真與實(shí)測(cè)數(shù)據(jù)分析

實(shí)驗(yàn)仿真主要分為兩個(gè)部分，一部分用于驗(yàn)證MUSIC算法的可行性，另一部分為實(shí)測(cè)探地雷達(dá)回波處理過(guò)程。

2.1 基于模擬數(shù)據(jù)的MUSIC算法驗(yàn)證

假設(shè)信號(hào)的采樣率為Fs=1000 Hz，信號(hào)長(zhǎng)度為N=2048，表達(dá)式為

利用MUSIC譜估計(jì)方法， 加入0.01倍信號(hào)噪聲后的測(cè)試結(jié)果如表1所示。

表1 仿真信號(hào)的參數(shù)估計(jì)值

由表1中列出的測(cè)試結(jié)果可以發(fā)現(xiàn)，原始信號(hào)與加入噪聲后的仿真結(jié)果基本一致，只有幅值和相位結(jié)果存在少許差別，這主要是因?yàn)樵肼暤拇嬖谑沟媒Y(jié)果產(chǎn)生了一些誤差。而在處理實(shí)際測(cè)量數(shù)據(jù)時(shí)，會(huì)采用第1節(jié)中的均值法和去直流偏移量進(jìn)行回波信號(hào)預(yù)處理，這樣可以有效地降低噪聲，而這種近似帶來(lái)的誤差是可以接受的。

2.2 探地雷達(dá)實(shí)測(cè)數(shù)據(jù)分析

如圖1所示，以探地雷達(dá)探測(cè)層狀結(jié)構(gòu)層為例來(lái)說(shuō)明雷達(dá)回波數(shù)據(jù)的處理過(guò)程。先采用消除水平同軸法與均值濾波法對(duì)回波信號(hào)進(jìn)行數(shù)據(jù)預(yù)處理，再采用MUSIC譜估計(jì)與反褶積方法進(jìn)行去噪聲處理，以期達(dá)到更好的效果。

圖1 探地雷達(dá)回波信號(hào)處理過(guò)程

采用的雷達(dá)設(shè)備為WGPR系列無(wú)線雷達(dá)，雷達(dá)中心頻率為1 GHz。為了盡可能降低外部環(huán)境干擾和人為干擾對(duì)測(cè)量結(jié)果的影響，實(shí)驗(yàn)中的沙箱尺寸為150 cm×50 cm×50 cm，將探地雷達(dá)懸掛在支架上，探地雷達(dá)距離沙子的距離為1.45 m，沙箱中沙子厚度為10 cm。圖2為探地雷達(dá)測(cè)量層狀結(jié)構(gòu)示意圖。

圖2 探地雷達(dá)測(cè)量層狀結(jié)構(gòu)示意圖

雷達(dá)采樣后得到的實(shí)時(shí)圖形如圖3所示，該段數(shù)據(jù)為10 cm砂層的雷達(dá)剖面圖，由450條采樣數(shù)據(jù)構(gòu)成，采樣點(diǎn)為 8192。其中，圖3中的右側(cè)波形即為采樣的單道回波信號(hào)，在實(shí)際數(shù)據(jù)處理時(shí)將提取每道回波的采樣值進(jìn)行數(shù)據(jù)預(yù)處理。

圖3 厚度為10 cm沙層的采樣數(shù)據(jù)

(1) 直流波去除和均值法濾波。

去直流波后利用均值法濾波將450條雷達(dá)回波進(jìn)行數(shù)據(jù)處理后，得到結(jié)果如圖4所示。

圖4 均值法濾波

(2) MUSIC譜估計(jì)去噪聲。

對(duì)圖4中進(jìn)行濾波后的信號(hào)運(yùn)用MUSIC法進(jìn)行功率譜估計(jì)，得到圖5中的功率譜曲線及圖6中的濾波結(jié)果，其中紅線代表0.2×max(20lg|PMUSIC(ω)|)的閾值線。

(3) 回波信號(hào)處理有效性驗(yàn)證。

利用重構(gòu)波形計(jì)算沙層的厚度來(lái)驗(yàn)證本文中回波信號(hào)處理過(guò)程的有效性。

圖5 MUSIC方法進(jìn)行譜估計(jì)的功率譜圖

圖6 MUSIC與均值法濾波比較圖

已知結(jié)構(gòu)層厚度的通用公式為

(27)

式中，c為電磁波在真空中的速度；Δt為探地雷達(dá)電磁波在該層的傳播時(shí)間;d為結(jié)構(gòu)層厚度；ε為該結(jié)構(gòu)層材料的介電常數(shù)值。

本實(shí)驗(yàn)中通過(guò)介電常數(shù)測(cè)量?jī)x來(lái)確定沙層的實(shí)際介電常數(shù)。取多次測(cè)量平均值為實(shí)際沙層的介電常數(shù)，大小為εsand=2.96。

分別對(duì)均值法濾波與高斯重構(gòu)回波數(shù)據(jù)進(jìn)行時(shí)間計(jì)算，計(jì)算得出沙層的傳播時(shí)間Δt約為108 ns，代入式(27)中，通過(guò)計(jì)算得到砂層厚度為9.493 cm，與沙層厚度實(shí)際值的絕對(duì)誤差為0.507，相對(duì)誤差為5.07%。而利用本文中方法進(jìn)行高斯信號(hào)回波重構(gòu)，計(jì)算得出雷達(dá)波在沙層中的傳播時(shí)間Δt為113 ns，通過(guò)計(jì)算得到沙層厚度為9.865 cm，與沙層厚度實(shí)際值的絕對(duì)誤差僅為0.135，相對(duì)誤差為1.35 %，誤差在5%之內(nèi)，相較于均值法濾波計(jì)算結(jié)果，精度得到了顯著提高，說(shuō)明本文中的回波處理方法有效，適用于實(shí)際測(cè)量。

3 結(jié)束語(yǔ)

針對(duì)探地雷達(dá)提出一種簡(jiǎn)單有效的處理回波信號(hào)方法，對(duì)雷達(dá)接收機(jī)接收到的回波信號(hào)進(jìn)行了時(shí)域和頻域聯(lián)合處理，先將回波信號(hào)進(jìn)行去直流波、均值濾波等回波后預(yù)處理，利用MUSIC譜估計(jì)方法估計(jì)各信號(hào)分量的頻率信息，再利用最小二乘法估計(jì)其幅值和相位，最終通過(guò)反褶積計(jì)算進(jìn)行波形重建，從而達(dá)到去除噪聲精確估計(jì)實(shí)用參數(shù)的目的。通過(guò)理論仿真和實(shí)測(cè)探地雷達(dá)數(shù)據(jù)處理，表明了該方法能很好地消除雷達(dá)回波數(shù)據(jù)的噪聲和雜波，能有效提高測(cè)量數(shù)據(jù)的精度，適合實(shí)際應(yīng)用。

圖7 均值法濾波與重構(gòu)回波時(shí)間計(jì)算