自適應(yīng)噪聲抑制器的設(shè)計和Simulink驗證

雷文英 ，黨瑞榮 ，牟明洋 ，戴乾軍

(1西安石油大學(xué)光電油氣測井與檢測教育部重點實驗室 陜西 西安 710065;2 中國石油集團測井有限公司長慶事業(yè)部 陜西 高陵 710020；3 中國石油集團測井有限公司吐哈事業(yè)部 新疆 鄯善 838202)

0 引言

測井儀器的傳感器采集的數(shù)據(jù)中往往含有大量的噪聲信息，這些疊加在有用信號中的噪聲信息使得人們很難直接從觀測信號中提取出反應(yīng)真實物理量的有用信息。為了從觀測信號中提取人們感興趣的有用數(shù)據(jù)，必須對傳感器采集的數(shù)據(jù)進行濾波處理濾除無用的噪聲信號，還原出被噪聲污染之前的原始信號。在測井儀器的接收傳感器信號的濾波處理的實踐中，由于油氣井下環(huán)境的復(fù)雜性，人們難以事先獲得噪聲的先驗知識，且噪聲可能是非平穩(wěn)的隨時間變化的，與此同時當(dāng)周圍的介質(zhì)發(fā)生變化時接收探頭中的感應(yīng)電動勢是未知的時變信號，被噪聲污染的原始信號往往微弱并且不穩(wěn)定，采用經(jīng)典的固定系數(shù)的FIR和IIR數(shù)字濾波器對接收的數(shù)據(jù)進行濾波處理效果都不理想，不能獲得令人滿意的濾波效果。針對此種情況，本文設(shè)計了一個自適應(yīng)濾波噪聲抑制系統(tǒng)，該系統(tǒng)實際上是一個基于ENLMS算法調(diào)整FIR濾波器的系數(shù)的自適應(yīng)濾波器，令濾波器的系數(shù)隨信號的變化和噪聲的變化而動態(tài)地進行調(diào)整，使其能夠從減輕傳感器接收信號中疊加的背景噪聲的影響，提高觀測信號的信噪比。

1 自適應(yīng)噪聲抑制的基本原理

自適應(yīng)噪聲抑制的基本原理是將被噪聲污染的信號與參考信號進行抵消運算，經(jīng)過自適應(yīng)濾波器系數(shù)的自動調(diào)整，使期望信號與濾波器輸出信號之間的均方誤差達到最小。實現(xiàn)自適應(yīng)濾波的算法很多[1]，由于應(yīng)用LMS算法比較簡便而且方便于FPGA采用FIR濾波器進行硬件實現(xiàn)，但同時由于最小均方LMS算法本身存在一些不足，故本文采用修正的歸一化ENLMS算法進行自適應(yīng)噪聲抑制。ENLMS算法的基礎(chǔ)是LMS算法，因此首先對LMS算法進行介紹[2]，然后推出ENLMS算法。對于M抽頭的FIR濾波器，設(shè)在第n時刻濾波器的輸入信號向量為X(n)，F(xiàn)IR濾波器的權(quán)系數(shù)向量為W(n)，濾波器的輸出信號為Y(n)，期望的輸出信號為d(n)則該自適應(yīng)濾波器的輸出信號y(n)，輸出信號y(n)與期望輸出信號d(n)的誤差e(n)，以及權(quán)系數(shù)更新之間的關(guān)系可以表示為[3]：

由于式(1)中的LMS算法中FIR濾波器的權(quán)值的修正量正比于輸入信號X(n)，當(dāng)輸入向量X(n)較大時，容易使修正量過大而造成噪聲放大，因此采用如(2)式所示的修正的ENLMS算法[4]對權(quán)值進行修正，其中ε是一個小的正數(shù)其目的是為了防止當(dāng)XT(n)X(n)過小時權(quán)值的修正量過大，給e(n)X(n)除以XT(n)X(n)+ε的目的是對調(diào)整量進行歸一化，μ 在這里代表歸一化權(quán)值。

圖1 FIR濾波器實現(xiàn)LMS算法的框圖[4]

一般的用于干擾消除的自適應(yīng)濾波器需要輸入兩路信號[5-6]，一路信號是含有加性噪聲的觀測信號，另一路是由于估計加性噪聲的參考信號。在測井中用于估計加性噪聲的參考信號往往難以獲得，這里采取一種無參考信號的自適應(yīng)FIR濾波器實現(xiàn)噪聲抑制。圖2示出了該噪聲抑制的原理框圖[7]，其基本原理是通過對圖2中的觀測信號d(n)進行一定的延時來產(chǎn)生一個參考信號x(n) 。圖2 中s(n) 是事先未知的原始信號，v(n) 是加性噪聲干擾信號，d(n) 是含有干擾信號的s(n) 的觀測值，e(n) 是觀測信號d(n) 與FIR 濾波器的輸出y(n) 的差值。均方誤差可以表示為：

實際中原始的信號s(n)一般是一個窄帶信號，噪聲信號v(n)是一個寬帶信號，且通常情況原始信號s(n)和噪聲信號v(n)之間互不相關(guān)。另外自適應(yīng)濾波器的輸入x(n)是d(n-L)，合理選擇L值可以使x(n)僅與s(n)相關(guān)而與噪聲v(n)不相關(guān)如(4)式所示：

故可以通過x(n)作為參考信號用于估計s(n)的值。x(n)經(jīng)過系數(shù)可變的FIR濾波器后其輸出y(n)可以表示為：在式(3)中由于v(n)和s(n)不相關(guān)故，EE{ {vv(n(n)s)(sn(n)}) }＝ = 0，0對于噪聲序列可以根據(jù)噪聲信號的自相關(guān)合理選取 L這個延遲量使E[v(n)y(n)]= 0，這時可以推出y(n)與v(n)不相關(guān)E[v(n)v(n?k?L)]= 0。此種情況下可以得出 式(3)E[e2(n)]=E[v2(n)]+E[s(n) ?y(n)]2，最小化均方誤差 等效于最小E[e2(n)]化E[s(n)-y(n)]2， 所以FIR濾波器的輸出y(n)就是s(n)的估計值，相應(yīng)的e(n)是噪聲信號v(n)的估計值。

根據(jù)經(jīng)驗，在接收的測井信號中夾雜的噪聲信號是有色噪聲而不是白噪聲，因此為了簡便起見用一個高斯白噪聲通過v(n) = 0.8v(n? 1 )+ σ(n)的AR(1)系統(tǒng)模擬產(chǎn)生該有色噪聲以便能反映真實的環(huán)境噪聲。

圖2 無參考信號的自適應(yīng)FIR噪聲抑制原理

2 自適應(yīng)FIR噪聲抑制器的Simulink驗證

要將上述自適應(yīng)噪聲抑制系統(tǒng)的第一步是編程仿真驗證該模型的正確性，本文將上述自適應(yīng)FIR噪聲抑制器在Simulink環(huán)境下進行建模仿真。一方面可以對上述噪聲抑制器的性能進行評估，另一方面可以對后續(xù)的硬件實現(xiàn)做好準(zhǔn)備。整個自適應(yīng)噪聲抑制仿真模型的各個模塊連接后的結(jié)構(gòu)如圖3所示，下面介紹具體的建模和分析過程。

圖3 自適應(yīng)噪聲抑制器的Simulink模型

2.1 自適應(yīng)噪聲抑制器Simulink模型的建立

首先啟動Simulink，新建一個名為adaptive.mdl的模型文件；然后給該模型文件中添加Sine Wave模塊用來表示原始輸入信號s(n)，Sine Wave模塊的參數(shù)設(shè)置為Aplitude=1V，F(xiàn)requency=25Hz，sample time=0.001s；接下來添加Random Source模塊用來產(chǎn)生均值為零的高斯白噪聲，Random Source模塊的參數(shù)設(shè)置為Source Type=Gaussian，Mean=0，

Variance=var， 具 體 的var的 值 由MATLAB Workspace中確定。在建立了高斯白噪聲Random Source后需要讓它通過一個AR(1)過程以產(chǎn)生非平穩(wěn)的有色噪聲[8]，該AR(1)過程用一個嵌入MATLAB函數(shù)實現(xiàn)。從Simulink Library Browser中找到Embedded MATLAB Function模塊添加到該模型文件中，雙擊該模塊編寫產(chǎn)生AR(1)過程的MATLAB代碼。Embedded MATLAB Function模塊實際上是Simulink中的一個MATLAB函數(shù)文件，該MATLAB函數(shù)的輸入?yún)?shù)自動作為Embedded MATLAB Function模塊的輸入，輸出參數(shù)自動作為Embedded MATLAB Function的輸出，產(chǎn)生有色噪聲v(n)的Embedded MATLAB Function1 模塊的MATLAB代碼如下：

function y = fColorNoise(u)

persistent y_old;

if(isempty(y_old))

y_old=0;

end

y = y_old*.8+u;

y_old=y;

下一步，一方面對有色噪聲進行數(shù)字特征分析，主要分析其均值和自相關(guān)，另一方面將該有色噪聲信號v(n)與原始信號s(n)疊加產(chǎn)生觀測信號d(n)。對有色噪聲進行數(shù)字特這分析的第一步需要將產(chǎn)生的有色噪聲通過一個output buffer size = 128，buffer overlap =0 初始條件為0的Buffer模塊，產(chǎn)生一個128點長度的噪聲序列，將該Buffer模塊的輸出作為Autocorrelation自相關(guān)模塊的輸入和均值Mean模塊的輸入。由于輸入到自相關(guān)模塊的有色噪聲是128個點的隨機序列，下一時刻的128點自相關(guān)勢必與當(dāng)前的自相關(guān)的結(jié)果不同但趨勢相同，為了能夠反映出自相關(guān)的變化趨勢，對該自相關(guān)序列的輸出加上均值模塊對前幾次的自相關(guān)結(jié)果取平均獲得自相關(guān)的變化趨勢。將fColorNoise模塊的輸出v(n)和Sine Wave模塊的輸出s(n)相加的結(jié)果d(n)作為觀測信號，同時將d(n)送入Integer Delay 模塊通過設(shè)置number of delays =L完成L點采樣的延遲作為自適應(yīng)FIR濾波器的輸入。

最后，完成整個自適應(yīng)噪聲抑制模型的核心部分---濾波器系數(shù)可變的FIR濾波器模塊的設(shè)計。該模塊同樣采用一段嵌入式MATLAB代碼來實現(xiàn)，主要實現(xiàn)ENLMS算法。在Simulink Library Browser 中拖動一個Embedded MATLAB function模塊到該Simulink文件中，在圖3中為Embedded MATLAB function2,雙擊該模塊編寫如下MATLAB代碼：

function [y,w,e]= fadaptive(x, x_delayed, mu, N)

persistent u;

persistent t;

persistent w_prev;

if(isempty(u))%初始化

u=zeros(N,1);

t=zeros(N,1);

w_prev=zeros(N,1);

end

u=[x_delayed;t(1:end-1)];

t=u;

y = w_prev'*u; % ENLMS自適應(yīng)算法

e = x-y;

w = w_prev+mu*e*u/(u'*u+.00001);

w_prev=w;

在上述fadaptive函數(shù)的輸入?yún)?shù)中，mu步長，N是FIR濾波器的階數(shù)，為了方便調(diào)試這兩個參數(shù)在仿真時通過MATLAB Workspace給出。Fadaptive函數(shù)的輸出值y對應(yīng)于圖2中的y(n)，輸出值w對應(yīng)于圖2中的權(quán)系數(shù)W，輸出值e對應(yīng)于圖2中的誤差e(n)。

2.2 自適應(yīng)噪聲抑制器Simulink模型的仿真

在連接好整個Simulink設(shè)計模型后，對上述Simulink模型設(shè)置適當(dāng)?shù)姆抡鎱?shù)后就可以對該模型進行仿真分析，進一步對該自適應(yīng)FIR濾波器的性能進行評估。在MATLAB Workspace中設(shè)置如表1所示的參數(shù)及參數(shù)值，然后選擇Simulation菜單下的Simulation Parameters選項設(shè)置仿真的結(jié)束時間為inf，讓模型開始仿真后持續(xù)運行直到用戶手動停止仿真，其他參數(shù)采用默認(rèn)設(shè)置。接下來選擇Simulating Start啟動模型仿真，雙擊圖3中的示波器模塊，波形顯示如圖4所示。

表1 自適應(yīng)噪聲抑制器Simulink模型中的參數(shù)

圖4中的第一個波形是含有噪聲的觀察信號，從該信號中可以看出經(jīng)過噪聲污染后的原始信號失真非常嚴(yán)重，已經(jīng)很難分辨出原始信號的特征；圖4中第二個波形是經(jīng)過自適應(yīng)FIR濾波器的輸出結(jié)果，可以看出該自適應(yīng)FIR濾波器的輸出結(jié)果與原始信號已經(jīng)非常接近，雖然在波形上有些失真，但是在不知道任何原始信號的先驗信息的情況下其效果已經(jīng)比較理想；圖4中第三個波形是圖2中d(n)與y(n)相減產(chǎn)生的誤差信號e(n)，可以看出該誤差信號是個噪聲信號。從圖4所示的仿真結(jié)果不但可以驗證圖2 所示的無參考信號的自適應(yīng)FIR噪聲抑制原理的正確性而且可以為下一步的不需要參考信號的自適應(yīng)噪聲抑制器的硬件實現(xiàn)提供理論依據(jù)。

圖4 自適應(yīng)噪聲抑制器仿真結(jié)果波形

噪聲序列的均值可以由圖3中的Vector Scope1可以觀測，雙擊圖3中的Vector Scope1彈出圖5上邊部分的圖像，該圖的橫軸是時間，縱軸是噪聲序列的均值，可以看出該噪聲信號的均值是與時間無關(guān)的量。噪聲序列的自相關(guān)的平均值變化趨勢可以從圖3中的Vector Scope中得到，如圖5下邊部分的圖像所示，從該圖像中可以看出噪聲信號當(dāng)時間間隔即圖2中的延遲量L約大于20后v(n)和v(n-L)不相關(guān)，即E[v(n)v(n-L)]=0，應(yīng)該注意的是該AR(1)噪聲是一個非平穩(wěn)隨時間變化的隨機信號，其自相關(guān)的結(jié)果同樣是實時波動變化的量，這里取其均值是為了看出自相關(guān)的變化趨勢。

圖5 噪聲信號的均值和自相關(guān)的均值

實際中按以上方案設(shè)計自適應(yīng)噪聲抑制器時，應(yīng)該注意兩點：第一應(yīng)該確保噪聲信號的自相關(guān)的變化趨勢在間隔一定的時間后變?yōu)?，記這個間隔為 ，這個要求一般的噪聲信號都能滿足；第二應(yīng)該確保圖2中的對觀測信號的延遲量L大于 ，這一點可以通過不斷的嘗試來選擇合適的值。

3 結(jié)束語

針對實際中電磁測井儀器測量信號中含有的噪聲采用固定系數(shù)的FIR或IIR濾波器濾除噪聲效果不理想的問題，本文提出了一種濾波器系數(shù)可變的自適應(yīng)FIR濾波器對噪聲進行抑制的方法，首先給出了該自適應(yīng)噪聲抑制器的基本原理和理論，然后給出了自適應(yīng)噪聲抑制器在Simulink中的建模和仿真過程，最后分析了仿真結(jié)果，對該自適應(yīng)噪聲抑制器的可行性進行了驗證。該方法相對于常見的自適應(yīng)濾波器進行噪聲抑制的方案不需要采集噪聲參考信號，僅僅利用觀測信號就可以完成對噪聲的抑制，而且噪聲抑制效果比較明顯，該模型的仿真結(jié)果不但可以為下一步用FPGA進行硬件實現(xiàn)提供理論依據(jù)，而且可以通過DSP Builder 結(jié)合該自適應(yīng)噪聲抑制器的Simulink仿真模型，將該自適應(yīng)FIR濾波器用DSP Builder提供的模塊進行實現(xiàn)[9]，將Simulink作為建模和實現(xiàn)環(huán)境，將由此DSP Builder模塊實現(xiàn)的自適應(yīng)濾波器的Simulink模型用Signal Compiler工具直接編譯為VHDL代碼和工具命令語言腳本TCL，對該VHDL語言代碼進行綜合[10]，最終可以完成FPGA中的硬件實現(xiàn)。

[1]黃建亮.基于FPGA的自適應(yīng)濾波器設(shè)計與實現(xiàn)[D].西安:西安電子科技大學(xué), 2006:9-26.

[2]張賢達. 現(xiàn)代信號處理[M]. 2版.北京:清華大學(xué)出版社,2002:190-192.

[3]高清運，李學(xué)初. 自適應(yīng)濾波器的FPGA實現(xiàn)[J].電子測量與儀器學(xué)報,2005,19(1):25.

[4]Alexander D. Poularikas, Zayed M. Ramadan. Adaptive Filtering Primer with MATLAB[M].CRC,2006:103,166.

[5]石艷麗.基于DSP的自適應(yīng)噪聲抵消系統(tǒng)研究[D].長春:長春理工大學(xué),2008:14-16.

[6]齊海兵.自適應(yīng)濾波器算法設(shè)計及其FPGA實現(xiàn)的研究與應(yīng)用[D].長沙：中南大學(xué),2006:14.

[7]楊綠溪.現(xiàn)代數(shù)字信號處理[M].北京:科學(xué)出版社,2007:321-324.

[8]王宏禹,邱天爽,陳喆. 非平穩(wěn)隨機信號分析與處理[M].2版.北京:機械工業(yè)出版社,2008:167-169.

[9]羅韓君,劉明偉,王成.基于DSPBuilder的FIR濾波器設(shè)計與實現(xiàn)[J].微計算機信息,2009,25(1).

[10]雷文英,黨瑞榮.多頻帶數(shù)字濾波器的設(shè)計和FPGA實現(xiàn)[J].電子測試,2009(11):29-34.