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

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

0 引言

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

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

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

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

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

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

實(shí)際中原始的信號(hào)s(n)一般是一個(gè)窄帶信號(hào)，噪聲信號(hào)v(n)是一個(gè)寬帶信號(hào)，且通常情況原始信號(hào)s(n)和噪聲信號(hào)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)作為參考信號(hào)用于估計(jì)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對(duì)于噪聲序列可以根據(jù)噪聲信號(hào)的自相關(guān)合理選取 L這個(gè)延遲量使E[v(n)y(n)]= 0，這時(shí)可以推出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)的估計(jì)值，相應(yīng)的e(n)是噪聲信號(hào)v(n)的估計(jì)值。

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

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

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

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

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

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

首先啟動(dòng)Simulink，新建一個(gè)名為adaptive.mdl的模型文件；然后給該模型文件中添加Sine Wave模塊用來表示原始輸入信號(hào)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后需要讓它通過一個(gè)AR(1)過程以產(chǎn)生非平穩(wěn)的有色噪聲[8]，該AR(1)過程用一個(gè)嵌入MATLAB函數(shù)實(shí)現(xiàn)。從Simulink Library Browser中找到Embedded MATLAB Function模塊添加到該模型文件中，雙擊該模塊編寫產(chǎn)生AR(1)過程的MATLAB代碼。Embedded MATLAB Function模塊實(shí)際上是Simulink中的一個(gè)MATLAB函數(shù)文件，該MATLAB函數(shù)的輸入?yún)?shù)自動(dòng)作為Embedded MATLAB Function模塊的輸入，輸出參數(shù)自動(dòng)作為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;

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

最后，完成整個(gè)自適應(yīng)噪聲抑制模型的核心部分---濾波器系數(shù)可變的FIR濾波器模塊的設(shè)計(jì)。該模塊同樣采用一段嵌入式MATLAB代碼來實(shí)現(xiàn)，主要實(shí)現(xiàn)ENLMS算法。在Simulink Library Browser 中拖動(dòng)一個(gè)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)試這兩個(gè)參數(shù)在仿真時(shí)通過MATLAB Workspace給出。Fadaptive函數(shù)的輸出值y對(duì)應(yīng)于圖2中的y(n)，輸出值w對(duì)應(yīng)于圖2中的權(quán)系數(shù)W，輸出值e對(duì)應(yīng)于圖2中的誤差e(n)。

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

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

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

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

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

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

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

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

3 結(jié)束語

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

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

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

[3]高清運(yùn)，李學(xué)初. 自適應(yīng)濾波器的FPGA實(shí)現(xiàn)[J].電子測(cè)量與儀器學(xué)報(bào),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è)計(jì)及其FPGA實(shí)現(xiàn)的研究與應(yīng)用[D].長沙：中南大學(xué),2006:14.

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

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

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

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