基于MUSIC譜估計的探地雷達回波處理方法

(長安大學 信息工程學院，陜西 西安 710064)

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

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

1 探地雷達數據處理方法

1.1 直流偏移量處理

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

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

(1)

式中，X(n)為未經處理的原始信號;n為采樣點數量;X′(n)為去除直流偏移量后得到的回波信號。

1.2 均值法濾波

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

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

(2)

依據一元函數求極值的原理可以得到

(3)

在隨機噪聲干擾的環境下，每次采樣接收的信號均可以表示為

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

(4)

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

探地雷達對同一目標多次采樣時，N次采樣值的算術平均值為

sk(n)+ek(n)]

(5)

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

(6)

(7)

假定采樣過程中，噪聲環境不發生變化，則噪聲ei(n)為相互獨立且概率密度相同的隨機變量。這樣根據噪聲ei(n)的性質，經過平均值濾波后的噪聲大小為

(8)

1.3 MUSIC功率譜估計

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

(9)

此時，探地雷達信號的濾波問題變為由已測得到的回波信號X(n)恢復出回波信號中的未知量Ai、ωi。下面采用多重信號分類的MUSIC算法來進行ωi估算。

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

依據MUSIC算法的工作原理，探地雷達信號X(n)的自相關函數可表示為

(10)

假設接收到的雷達信號可表示為N個采樣值的組合，則信號矢量可表示為

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

(11)

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

(12)

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

Rx=Rs+Rw

(13)

將矩陣Rs特征分解得

(14)

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

(15)

對于噪聲相關矩陣也可以表示為

(16)

式(10)可以表示為

(17)

(18)

實際做法是構造函數，構造的函數為

(19)

設當ω=ωi時功率譜為無窮大，此時相關矩陣估計值為ωi，為PMUSIC(ω)取到最大值時所對應的值。由于相關矩陣是估計值，必然存在誤差，因此PMUSIC(ωi)為有限值，但會呈現出很尖的峰值，此時峰值對應的頻率就是復正弦信號的頻率。

1.4 計算幅值與相角

信號的最佳頻率ωi和正弦信號個數M可以由MUSIC譜方法準確地估算出來，但是通過該方法無法檢測出正弦波的幅值和相角，因此需要結合其他方法進行信號參數計算。主要采用求解過程簡單、計算量較小的最小二乘法[12]來估算振幅值和相角。

假設已知信號記為X(t)，振幅和初相分別為Ai和φi，最佳組成頻率成分為fi，其中fi=2π/ωi,噪聲為w(t)。此時信號的數學關系表達式為

(20)

由式(20)可知，雷達信號經過時域n點采樣后可以表達為

(21)

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

A=

當噪聲較小時，式(21)可寫成如下形式：

B=AX

(22)

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

X=A-1B

(23)

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

X=A+B

(24)

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

(25)

(26)

2 仿真與實測數據分析

實驗仿真主要分為兩個部分，一部分用于驗證MUSIC算法的可行性，另一部分為實測探地雷達回波處理過程。

2.1 基于模擬數據的MUSIC算法驗證

假設信號的采樣率為Fs=1000 Hz，信號長度為N=2048，表達式為

利用MUSIC譜估計方法， 加入0.01倍信號噪聲后的測試結果如表1所示。

表1 仿真信號的參數估計值

由表1中列出的測試結果可以發現，原始信號與加入噪聲后的仿真結果基本一致，只有幅值和相位結果存在少許差別，這主要是因為噪聲的存在使得結果產生了一些誤差。而在處理實際測量數據時，會采用第1節中的均值法和去直流偏移量進行回波信號預處理，這樣可以有效地降低噪聲，而這種近似帶來的誤差是可以接受的。

2.2 探地雷達實測數據分析

如圖1所示，以探地雷達探測層狀結構層為例來說明雷達回波數據的處理過程。先采用消除水平同軸法與均值濾波法對回波信號進行數據預處理，再采用MUSIC譜估計與反褶積方法進行去噪聲處理，以期達到更好的效果。

圖1 探地雷達回波信號處理過程

采用的雷達設備為WGPR系列無線雷達，雷達中心頻率為1 GHz。為了盡可能降低外部環境干擾和人為干擾對測量結果的影響，實驗中的沙箱尺寸為150 cm×50 cm×50 cm，將探地雷達懸掛在支架上，探地雷達距離沙子的距離為1.45 m，沙箱中沙子厚度為10 cm。圖2為探地雷達測量層狀結構示意圖。

圖2 探地雷達測量層狀結構示意圖

雷達采樣后得到的實時圖形如圖3所示，該段數據為10 cm砂層的雷達剖面圖，由450條采樣數據構成，采樣點為 8192。其中，圖3中的右側波形即為采樣的單道回波信號，在實際數據處理時將提取每道回波的采樣值進行數據預處理。

圖3 厚度為10 cm沙層的采樣數據

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

去直流波后利用均值法濾波將450條雷達回波進行數據處理后，得到結果如圖4所示。

圖4 均值法濾波

(2) MUSIC譜估計去噪聲。

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

(3) 回波信號處理有效性驗證。

利用重構波形計算沙層的厚度來驗證本文中回波信號處理過程的有效性。

圖5 MUSIC方法進行譜估計的功率譜圖

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

已知結構層厚度的通用公式為

(27)

式中，c為電磁波在真空中的速度；Δt為探地雷達電磁波在該層的傳播時間;d為結構層厚度；ε為該結構層材料的介電常數值。

本實驗中通過介電常數測量儀來確定沙層的實際介電常數。取多次測量平均值為實際沙層的介電常數，大小為εsand=2.96。

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

3 結束語

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

圖7 均值法濾波與重構回波時間計算