何臘梅 楊帆

摘 要：以 LMS自適應(yīng)濾波理論為基礎(chǔ)，借助Matlab編程軟件成功濾除了航空瞬變電磁信號(hào)中的噪聲，為航空瞬變電磁信號(hào)的后期正反演研究提供了可靠的數(shù)據(jù)。通過(guò)實(shí)驗(yàn)證明，LMS自適應(yīng)濾波器能夠合理、有效地去除天電噪聲和地質(zhì)噪聲。

關(guān)鍵詞：航空瞬變電磁信號(hào)；數(shù)據(jù)去噪；濾波方法；LMS自適應(yīng)濾波

中圖分類號(hào)：TP39；P631.3+26 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：2095-1302（2018）01-00-03

0 引 言

航空瞬變電磁法（Airborne Transient Electromagnetic Method，ATEM）是20世紀(jì)中期問(wèn)世的一種快速普查良導(dǎo)電金屬礦的航空物探方法，其具有速度快，勘探成本低，探測(cè)范圍廣，可有效抑制復(fù)雜地形影響等優(yōu)點(diǎn)，已成為國(guó)內(nèi)外廣泛使用的一種地質(zhì)勘探方法[1]。但航空瞬變電磁勘探法主要觀測(cè)的是二次場(chǎng)信號(hào)，其有效信號(hào)幅值弱，頻帶寬[2]，因而測(cè)得的二次場(chǎng)信號(hào)的真實(shí)性或準(zhǔn)確性難以保證。此外，航空瞬變電磁法在實(shí)際應(yīng)用中受噪聲干擾的影響更為嚴(yán)重，甚至得到的觀測(cè)信號(hào)是已被噪聲信號(hào)掩蓋的無(wú)用信號(hào)[3]，導(dǎo)致后期處理進(jìn)入錯(cuò)誤的方向。為保證實(shí)測(cè)信號(hào)的可靠性，對(duì)測(cè)得的數(shù)據(jù)必須先做濾波、去噪等預(yù)處理。

航空瞬變電磁信號(hào)的噪聲類型分為天電噪聲、儀器振動(dòng)噪聲、地質(zhì)噪聲以及人文噪聲等[4]。由于平均濾波、中值濾波等濾波方法對(duì)天電噪聲的濾除不夠理想，本文提出運(yùn)用LMS自適應(yīng)濾波方式濾除其高頻范圍內(nèi)的天電干擾和地質(zhì)噪聲。

1 LMS自適應(yīng)濾波原理

自適應(yīng)濾波（adaptive filtering）是信號(hào)處理領(lǐng)域一個(gè)非常重要的分支。自1959年Widrow提出自適應(yīng)的概念以來(lái)，自適應(yīng)濾波理論一直受到普遍關(guān)注，并得到了不斷發(fā)展與完善。信號(hào)處理理論和應(yīng)用的發(fā)展為自適應(yīng)濾波理論提供了必要的理論基礎(chǔ)，其已在通信、雷達(dá)、自動(dòng)控制、圖像與語(yǔ)音處理等領(lǐng)域得到了廣泛應(yīng)用[5]。

顧名思義，自適應(yīng)濾波器是一種能夠根據(jù)輸入信號(hào)自動(dòng)調(diào)整自身性能并進(jìn)行數(shù)字信號(hào)處理的數(shù)字濾波器，其最本質(zhì)的特點(diǎn)在于具有自學(xué)習(xí)和自調(diào)整的能力，即自適應(yīng)能力。與固定濾波器相比，自適應(yīng)濾波器能夠根據(jù)當(dāng)前自身的狀態(tài)和環(huán)境自動(dòng)調(diào)整和校正當(dāng)前濾波器的參數(shù)。由于信號(hào)與噪聲的時(shí)變特性未知，因而其統(tǒng)計(jì)特性不確定，故可認(rèn)為其是最優(yōu)濾波方法。自適應(yīng)濾波器具有可調(diào)整系數(shù)的濾波結(jié)構(gòu)與可調(diào)整和校正濾波器系數(shù)的自適應(yīng)算法[6]。

Widrow等人提出的最小均方算法（Least Mean Square，LMS）是一種以期望響應(yīng)和濾波器輸出信號(hào)之間誤差的均方值最小為準(zhǔn)則，依據(jù)輸入信號(hào)在迭代過(guò)程中估計(jì)梯度質(zhì)量，并更新權(quán)系數(shù)以達(dá)到最優(yōu)的自適應(yīng)迭代算法。LMS算法是一種梯度最速下降算法，其顯著特點(diǎn)在于其簡(jiǎn)單性。這種算法無(wú)需計(jì)算相關(guān)矩陣及矩陣運(yùn)算[5]。

LMS算法是一種線性自適應(yīng)濾波算法，包括兩個(gè)基本過(guò)程，即濾波過(guò)程與自適應(yīng)過(guò)程。在濾波過(guò)程中，自適應(yīng)濾波器計(jì)算其對(duì)輸入的響應(yīng)，并通過(guò)與期望響應(yīng)比較，得到估計(jì)的誤差信號(hào)。在自適應(yīng)過(guò)程中，系統(tǒng)估計(jì)誤差自動(dòng)調(diào)整濾波器自身的參數(shù)。這兩個(gè)過(guò)程共同組成一個(gè)反饋環(huán)，LMS自適應(yīng)濾波原理如圖1所示。

其中，s （k）是輸入信號(hào)，d（k）是期望信號(hào)，y（k）是濾波后的信號(hào)，e（k）是誤差信號(hào)。

該濾波器根據(jù)e（k）和s（k），通過(guò)自適應(yīng)算法找到E[e2（k）]最小時(shí)濾波器的權(quán)值，從而找到濾波系數(shù)權(quán)值，實(shí)現(xiàn)自適應(yīng)濾波。

LMS算法系數(shù)的更新規(guī)則是通過(guò)迭代尋找極值，即最速下降法，該方法后一時(shí)刻的系數(shù)由前一時(shí)刻的系數(shù)通過(guò)一定的迭代運(yùn)算得出：

式中，W為濾波系數(shù)，k為當(dāng)前迭代值，Wk和Wk+1分別代表k時(shí)刻和k+1時(shí)刻的自適應(yīng)系數(shù)，X為待濾波數(shù)據(jù)，μ為放大系數(shù)，稱為收斂因子，用于調(diào)整自適應(yīng)迭代的步長(zhǎng)，ε是誤差信號(hào)（期望值與預(yù)測(cè)輸出值之差），為誤差信號(hào)平方的梯度。假設(shè)當(dāng)前為k時(shí)刻，其后一時(shí)刻為k+1時(shí)刻，則k+1時(shí)刻的系數(shù)為k時(shí)刻系數(shù)加上歸一化函數(shù)的負(fù)梯度。系數(shù)函數(shù)也稱為性能函數(shù)，等同于期望信號(hào)的均方誤差。從（3）式可得出，LMS算法與平方、平均、矩陣逆變換無(wú)關(guān)，但其收斂判斷標(biāo)準(zhǔn)與最速下降方法確定的值相等。對(duì)μ的標(biāo)準(zhǔn)由下式得出：

n為濾波系數(shù)的階數(shù)，R為輸入的相關(guān)矩陣E[XkXkT]，tr為矩陣跡。濾波系數(shù)的階數(shù)n和濾波增益μ是失調(diào)M和自適應(yīng)處理時(shí)間常數(shù)τmse的折中。

2 LMS自適應(yīng)濾波實(shí)現(xiàn)流程

圖2是運(yùn)用LMS自適應(yīng)濾波算法對(duì)航空瞬變電磁去噪的實(shí)現(xiàn)流程[7]。由圖可知，讀入數(shù)據(jù)后，首先計(jì)算信號(hào)長(zhǎng)度和功率，然后計(jì)算出步長(zhǎng)值并將該值作為初始步長(zhǎng)存入算法中，設(shè)置好前兩個(gè)抽頭系數(shù)值后進(jìn)入算法的自適應(yīng)計(jì)算環(huán)節(jié)，首先計(jì)算前兩個(gè)值，之后再計(jì)算誤差與下一個(gè)抽頭系數(shù)，依次循環(huán)，直到最后一個(gè)值為止，從而得到濾波后的信號(hào)。

3 LMS自適應(yīng)濾波應(yīng)用

1988年，Spies運(yùn)用預(yù)測(cè)（非自適應(yīng)）技術(shù)（或稱局部噪聲預(yù)測(cè)濾波）來(lái)估計(jì)垂直磁場(chǎng)的噪聲[8]。固定系數(shù)的消噪方法在去除時(shí)不變（平穩(wěn)）噪聲時(shí)效果較好，但當(dāng)前的電磁數(shù)據(jù)噪聲幾乎都是非平穩(wěn)噪聲，因此這些濾波方法力不從心。本文結(jié)合天然磁場(chǎng)的性質(zhì)和LMS自適應(yīng)濾波算法的特性，運(yùn)用LMS自適應(yīng)濾波算法對(duì)由模擬勘探模型計(jì)算得出的數(shù)據(jù)和實(shí)測(cè)的航空電磁數(shù)據(jù)消噪[9]。

3.1 LMS自適應(yīng)濾波算法對(duì)模擬信號(hào)去噪

讀入模擬數(shù)據(jù)后就可運(yùn)用該濾波算法濾除噪聲，為了更好地濾波，將信號(hào)分為早期和后期兩部分，分別進(jìn)行LMS自適應(yīng)濾波，程序如下：

//加高斯白噪聲

snr=5 ；

sigpower=measured；

x=awgn（x0，snr，sigpower）； //x0為模擬的原始信號(hào)，x為加入高斯白噪聲后的信號(hào)

N=length（x）；

n=randn（1，N）；

u=0.03；

w（1：N）=0；

w（1）=2；

for i=1：20//將信號(hào)分為兩部分，分別運(yùn)用LMS自適應(yīng)算法

y（i）=n（i）*w（i）；

e（i）=x（i）-y（i）；

w（i+1）=w（i）+u*e（i）*n（i）；

end；

u=0.01；

for i=21：N

y（i）=n（i）*w（i）；

e（i）=x（i）-y（i）；

w（i+1）=w（i）+u*e（i）*n（i）；

end；

圖3所示為L(zhǎng)MS自適應(yīng)濾波方法對(duì)含噪模擬數(shù)據(jù)的濾波結(jié)果。

由圖可知，該方法早期效果較好，但在早期向中期過(guò)渡時(shí)，由于幅度變化太快，使得自適應(yīng)濾波的效果不太理想，到了中后期濾波后曲線比較光滑。LMS濾波方法可以較好地去除干擾信號(hào)，擬合曲線，可有效保證信號(hào)幅值不被削弱。

3.2 LMS自適應(yīng)濾波算法對(duì)實(shí)測(cè)信號(hào)去噪

分兩部分運(yùn)用LMS自適應(yīng)算法，濾波程序如下：

N=length（x0）； //信號(hào)長(zhǎng)度

t=1：N；

Nw=21；

fs=50；

nfft=N；

ff=（0：N/2-1）*fs/（N-1）；

x1=x0；

snr=5 ；

sigpower=measured；

n（1：N1）=0.5；

n=awgn（n，snr，sigpower，linear）；

u=0.005；

w（1：N1）=0；

w（1）=5.5；

for i=1：50

y（i）=n（i）*w（i）；

e（i）=x1（i）-y（i）；

w（i+1）=w（i）+u*e（i）*n（i）；

end；

u=0.33；

for i=51：N1

y（i）=n（i）*w（i）；

e（i）=x1（i）-y（i）；

w（i+1）=w（i）+u*e（i）*n（i）；

end；

圖4所示為L(zhǎng)MS自適應(yīng)濾波方法對(duì)實(shí)測(cè)數(shù)據(jù)濾波的結(jié)果。從圖中可知，50～100點(diǎn)間發(fā)生了較大跳變，無(wú)法較好地反映原始數(shù)據(jù)特性，但中后期處理效果較好，雖然曲線仍有毛刺，但與原始數(shù)據(jù)相比曲線光滑很多，可認(rèn)為已基本實(shí)現(xiàn)了去噪效果。濾波后早期信號(hào)與原始早期信號(hào)擬合較好，但毛刺較多，處理效果并不理想，但這由LMS自適應(yīng)算法本身的原因?qū)е?。LMS自適應(yīng)方法的關(guān)鍵步驟在于自適應(yīng)算子的計(jì)算。但從整體消噪效果和保幅性方面看，已基本去除干擾成分，較好地保留了有用信號(hào)。

4 結(jié) 語(yǔ)

綜上所述，LMS自適應(yīng)算法能夠快速進(jìn)行自學(xué)習(xí)和自調(diào)整，適合處理非平穩(wěn)的干擾信號(hào)，且算法實(shí)現(xiàn)簡(jiǎn)單。航空瞬變電磁信號(hào)中最主要的干擾是天電噪聲，通過(guò)運(yùn)用LMS自適應(yīng)濾波算法能夠較好地濾除這一主要噪聲信號(hào)，因而能夠做好航空瞬變電磁信號(hào)前期的預(yù)處理工作，為后期分析提供可靠的數(shù)據(jù)。

參考文獻(xiàn)

[1]雷棟，胡祥云，張素芳.航空電磁法的發(fā)展現(xiàn)狀[J].地質(zhì)找礦論叢，2006，21（1）：40-44.

[2]嵇艷鞠，林君，程德福，等.瞬變電磁法中數(shù)據(jù)取樣處理方法的研究[J].物探與化探，2003， 27（2）： 142-145.

[3]何臘梅.航空瞬變電磁數(shù)據(jù)時(shí)頻分析及去噪方法研究[D].成都： 成都理工大學(xué)，2013.

[4] Jewell T R，Ward S H. The influence of conductivity inhomogeneities upon audio-frequency magnetic fields[J].Geophysics， 1963，28（2）：201-221.

[5]邱天爽，郭瑩.信號(hào)處理與數(shù)據(jù)分析[M].北京：清華大學(xué)出版社，2015.

[6]何臘梅，王宇航.STFT與FIR在航空瞬變電磁數(shù)據(jù)處理中的應(yīng)用[J].工程地球物理學(xué)報(bào)，2013，10（1）： 15-20.

[7] Buselli G， Hwang H S， Pik J P. AEM noise reduction with remote referencing[J].Exploration Geophysics， 1998， 29（2）： 71-76.

[8] Spies B R. Local noise prediction filtering for central induction transient electromagnetic sounding[J]. Geophysics， 1988， 53（8）： 1068-1079.

[9] Jeng Y， Li Y W， Chen C S，et al. Adaptive filtering of random noise in near-surface seismic and ground-penetrating radar data[J].Journal of applied geophysics，2009，68：36-46.