自适应噪声抑制器的设计和Simulink验证

雷文英 ，党瑞荣 ，牟明洋 ，戴乾军

(1西安石油大学光电油气测井与检测教育部重点实验室 陕西 西安 710065;2 中国石油集团测井有限公司长庆事业部 陕西 高陵 710020；3 中国石油集团测井有限公司吐哈事业部 新疆 鄯善 838202)

0 引言

测井仪器的传感器采集的数据中往往含有大量的噪声信息，这些叠加在有用信号中的噪声信息使得人们很难直接从观测信号中提取出反应真实物理量的有用信息。为了从观测信号中提取人们感兴趣的有用数据，必须对传感器采集的数据进行滤波处理滤除无用的噪声信号，还原出被噪声污染之前的原始信号。在测井仪器的接收传感器信号的滤波处理的实践中，由于油气井下环境的复杂性，人们难以事先获得噪声的先验知识，且噪声可能是非平稳的随时间变化的，与此同时当周围的介质发生变化时接收探头中的感应电动势是未知的时变信号，被噪声污染的原始信号往往微弱并且不稳定，采用经典的固定系数的FIR和IIR数字滤波器对接收的数据进行滤波处理效果都不理想，不能获得令人满意的滤波效果。针对此种情况，本文设计了一个自适应滤波噪声抑制系统，该系统实际上是一个基于ENLMS算法调整FIR滤波器的系数的自适应滤波器，令滤波器的系数随信号的变化和噪声的变化而动态地进行调整，使其能够从减轻传感器接收信号中叠加的背景噪声的影响，提高观测信号的信噪比。

1 自适应噪声抑制的基本原理

自适应噪声抑制的基本原理是将被噪声污染的信号与参考信号进行抵消运算，经过自适应滤波器系数的自动调整，使期望信号与滤波器输出信号之间的均方误差达到最小。实现自适应滤波的算法很多[1]，由于应用LMS算法比较简便而且方便于FPGA采用FIR滤波器进行硬件实现，但同时由于最小均方LMS算法本身存在一些不足，故本文采用修正的归一化ENLMS算法进行自适应噪声抑制。ENLMS算法的基础是LMS算法，因此首先对LMS算法进行介绍[2]，然后推出ENLMS算法。对于M抽头的FIR滤波器，设在第n时刻滤波器的输入信号向量为X(n)，FIR滤波器的权系数向量为W(n)，滤波器的输出信号为Y(n)，期望的输出信号为d(n)则该自适应滤波器的输出信号y(n)，输出信号y(n)与期望输出信号d(n)的误差e(n)，以及权系数更新之间的关系可以表示为[3]：

由于式(1)中的LMS算法中FIR滤波器的权值的修正量正比于输入信号X(n)，当输入向量X(n)较大时，容易使修正量过大而造成噪声放大，因此采用如(2)式所示的修正的ENLMS算法[4]对权值进行修正，其中ε是一个小的正数其目的是为了防止当XT(n)X(n)过小时权值的修正量过大，给e(n)X(n)除以XT(n)X(n)+ε的目的是对调整量进行归一化，μ 在这里代表归一化权值。

图1 FIR滤波器实现LMS算法的框图[4]

一般的用于干扰消除的自适应滤波器需要输入两路信号[5-6]，一路信号是含有加性噪声的观测信号，另一路是由于估计加性噪声的参考信号。在测井中用于估计加性噪声的参考信号往往难以获得，这里采取一种无参考信号的自适应FIR滤波器实现噪声抑制。图2示出了该噪声抑制的原理框图[7]，其基本原理是通过对图2中的观测信号d(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)之间互不相关。另外自适应滤波器的输入x(n)是d(n-L)，合理选择L值可以使x(n)仅与s(n)相关而与噪声v(n)不相关如(4)式所示：

故可以通过x(n)作为参考信号用于估计s(n)的值。x(n)经过系数可变的FIR滤波器后其输出y(n)可以表示为：在式(3)中由于v(n)和s(n)不相关故，EE{ {vv(n(n)s)(sn(n)}) }＝ = 0，0对于噪声序列可以根据噪声信号的自相关合理选取 L这个延迟量使E[v(n)y(n)]= 0，这时可以推出y(n)与v(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)的估计值，相应的e(n)是噪声信号v(n)的估计值。

根据经验，在接收的测井信号中夹杂的噪声信号是有色噪声而不是白噪声，因此为了简便起见用一个高斯白噪声通过v(n) = 0.8v(n− 1 )+ σ(n)的AR(1)系统模拟产生该有色噪声以便能反映真实的环境噪声。

图2 无参考信号的自适应FIR噪声抑制原理

2 自适应FIR噪声抑制器的Simulink验证

要将上述自适应噪声抑制系统的第一步是编程仿真验证该模型的正确性，本文将上述自适应FIR噪声抑制器在Simulink环境下进行建模仿真。一方面可以对上述噪声抑制器的性能进行评估，另一方面可以对后续的硬件实现做好准备。整个自适应噪声抑制仿真模型的各个模块连接后的结构如图3所示，下面介绍具体的建模和分析过程。

图3 自适应噪声抑制器的Simulink模型

2.1 自适应噪声抑制器Simulink模型的建立

首先启动Simulink，新建一个名为adaptive.mdl的模型文件；然后给该模型文件中添加Sine Wave模块用来表示原始输入信号s(n)，Sine Wave模块的参数设置为Aplitude=1V，Frequency=25Hz，sample time=0.001s；接下来添加Random Source模块用来产生均值为零的高斯白噪声，Random Source模块的参数设置为Source Type=Gaussian，Mean=0，

Variance=var， 具 体 的var的 值 由MATLAB Workspace中确定。在建立了高斯白噪声Random Source后需要让它通过一个AR(1)过程以产生非平稳的有色噪声[8]，该AR(1)过程用一个嵌入MATLAB函数实现。从Simulink Library Browser中找到Embedded MATLAB Function模块添加到该模型文件中，双击该模块编写产生AR(1)过程的MATLAB代码。Embedded MATLAB Function模块实际上是Simulink中的一个MATLAB函数文件，该MATLAB函数的输入参数自动作为Embedded MATLAB Function模块的输入，输出参数自动作为Embedded MATLAB Function的输出，产生有色噪声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;

下一步，一方面对有色噪声进行数字特征分析，主要分析其均值和自相关，另一方面将该有色噪声信号v(n)与原始信号s(n)叠加产生观测信号d(n)。对有色噪声进行数字特这分析的第一步需要将产生的有色噪声通过一个output buffer size = 128，buffer overlap =0 初始条件为0的Buffer模块，产生一个128点长度的噪声序列，将该Buffer模块的输出作为Autocorrelation自相关模块的输入和均值Mean模块的输入。由于输入到自相关模块的有色噪声是128个点的随机序列，下一时刻的128点自相关势必与当前的自相关的结果不同但趋势相同，为了能够反映出自相关的变化趋势，对该自相关序列的输出加上均值模块对前几次的自相关结果取平均获得自相关的变化趋势。将fColorNoise模块的输出v(n)和Sine Wave模块的输出s(n)相加的结果d(n)作为观测信号，同时将d(n)送入Integer Delay 模块通过设置number of delays =L完成L点采样的延迟作为自适应FIR滤波器的输入。

最后，完成整个自适应噪声抑制模型的核心部分---滤波器系数可变的FIR滤波器模块的设计。该模块同样采用一段嵌入式MATLAB代码来实现，主要实现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自适应算法

e = x-y;

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

w_prev=w;

在上述fadaptive函数的输入参数中，mu步长，N是FIR滤波器的阶数，为了方便调试这两个参数在仿真时通过MATLAB Workspace给出。Fadaptive函数的输出值y对应于图2中的y(n)，输出值w对应于图2中的权系数W，输出值e对应于图2中的误差e(n)。

2.2 自适应噪声抑制器Simulink模型的仿真

在连接好整个Simulink设计模型后，对上述Simulink模型设置适当的仿真参数后就可以对该模型进行仿真分析，进一步对该自适应FIR滤波器的性能进行评估。在MATLAB Workspace中设置如表1所示的参数及参数值，然后选择Simulation菜单下的Simulation Parameters选项设置仿真的结束时间为inf，让模型开始仿真后持续运行直到用户手动停止仿真，其他参数采用默认设置。接下来选择Simulating Start启动模型仿真，双击图3中的示波器模块，波形显示如图4所示。

表1 自适应噪声抑制器Simulink模型中的参数

图4中的第一个波形是含有噪声的观察信号，从该信号中可以看出经过噪声污染后的原始信号失真非常严重，已经很难分辨出原始信号的特征；图4中第二个波形是经过自适应FIR滤波器的输出结果，可以看出该自适应FIR滤波器的输出结果与原始信号已经非常接近，虽然在波形上有些失真，但是在不知道任何原始信号的先验信息的情况下其效果已经比较理想；图4中第三个波形是图2中d(n)与y(n)相减产生的误差信号e(n)，可以看出该误差信号是个噪声信号。从图4所示的仿真结果不但可以验证图2 所示的无参考信号的自适应FIR噪声抑制原理的正确性而且可以为下一步的不需要参考信号的自适应噪声抑制器的硬件实现提供理论依据。

图4 自适应噪声抑制器仿真结果波形

噪声序列的均值可以由图3中的Vector Scope1可以观测，双击图3中的Vector Scope1弹出图5上边部分的图像，该图的横轴是时间，纵轴是噪声序列的均值，可以看出该噪声信号的均值是与时间无关的量。噪声序列的自相关的平均值变化趋势可以从图3中的Vector Scope中得到，如图5下边部分的图像所示，从该图像中可以看出噪声信号当时间间隔即图2中的延迟量L约大于20后v(n)和v(n-L)不相关，即E[v(n)v(n-L)]=0，应该注意的是该AR(1)噪声是一个非平稳随时间变化的随机信号，其自相关的结果同样是实时波动变化的量，这里取其均值是为了看出自相关的变化趋势。

图5 噪声信号的均值和自相关的均值

实际中按以上方案设计自适应噪声抑制器时，应该注意两点：第一应该确保噪声信号的自相关的变化趋势在间隔一定的时间后变为0，记这个间隔为 ，这个要求一般的噪声信号都能满足；第二应该确保图2中的对观测信号的延迟量L大于 ，这一点可以通过不断的尝试来选择合适的值。

3 结束语

针对实际中电磁测井仪器测量信号中含有的噪声采用固定系数的FIR或IIR滤波器滤除噪声效果不理想的问题，本文提出了一种滤波器系数可变的自适应FIR滤波器对噪声进行抑制的方法，首先给出了该自适应噪声抑制器的基本原理和理论，然后给出了自适应噪声抑制器在Simulink中的建模和仿真过程，最后分析了仿真结果，对该自适应噪声抑制器的可行性进行了验证。该方法相对于常见的自适应滤波器进行噪声抑制的方案不需要采集噪声参考信号，仅仅利用观测信号就可以完成对噪声的抑制，而且噪声抑制效果比较明显，该模型的仿真结果不但可以为下一步用FPGA进行硬件实现提供理论依据，而且可以通过DSP Builder 结合该自适应噪声抑制器的Simulink仿真模型，将该自适应FIR滤波器用DSP Builder提供的模块进行实现[9]，将Simulink作为建模和实现环境，将由此DSP Builder模块实现的自适应滤波器的Simulink模型用Signal Compiler工具直接编译为VHDL代码和工具命令语言脚本TCL，对该VHDL语言代码进行综合[10]，最终可以完成FPGA中的硬件实现。

[1]黄建亮.基于FPGA的自适应滤波器设计与实现[D].西安:西安电子科技大学, 2006:9-26.

[2]张贤达. 现代信号处理[M]. 2版.北京:清华大学出版社,2002:190-192.

[3]高清运，李学初. 自适应滤波器的FPGA实现[J].电子测量与仪器学报,2005,19(1):25.

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

[5]石艳丽.基于DSP的自适应噪声抵消系统研究[D].长春:长春理工大学,2008:14-16.

[6]齐海兵.自适应滤波器算法设计及其FPGA实现的研究与应用[D].长沙：中南大学,2006:14.

[7]杨绿溪.现代数字信号处理[M].北京:科学出版社,2007:321-324.

[8]王宏禹,邱天爽,陈喆. 非平稳随机信号分析与处理[M].2版.北京:机械工业出版社,2008:167-169.

[9]罗韩君,刘明伟,王成.基于DSPBuilder的FIR滤波器设计与实现[J].微计算机信息,2009,25(1).

[10]雷文英,党瑞荣.多频带数字滤波器的设计和FPGA实现[J].电子测试,2009(11):29-34.