8.1.背景

由于发动机试车过程具有较高的危险性,测控室往往与试车台相距十几米甚至更远,通过测控线路将试验控制和数据采集设备与试车台上的传感器和控制器相连。测控线路带来的信号污染加上充斥着各种电子设备的复杂电磁环境的影响,致使测试数据中往往不可避免地存在一定量的噪声,这些噪声不但使试验数据进一步偏离实际值,还会改变其数学形态,如光滑性甚至是连续性,对一些包含梯度和导数求取运算的算法实施造成了障碍。因此,本系统中实现了滤波和去噪的功能,能够为后续计算分析的开展奠定良好基础。图1为某试验发动机推力-时间曲线,该曲线存在多处明显振荡,可见噪声的影响较为显著。

在发动机研制过程中的方案阶段,需要设计、制造厚壁试验发动机或全尺寸发动机,并进行静止试验检验内弹道性能。在随后的初样阶段,也即技术设计阶段,要试制产品发动机,完成发动机各项地面综合试验及点火试车,全面检验发动机结构及各项技术指标,修改设计。每次试车试验结束后,都需根据数据处理规范对所得数据进行处理,得到发动机内弹道性能参数。其中广泛采用的标准为QJ 1047-92“固体火箭发动机压强-时间推力-时间数据处理规范”。该标准制定推出时,基于计算机的数据处理方式尚未在各院所推广开来,标准中一些特征点的选取通过人工作图确定,最终处理结果受人工误差影响较大;另外,每次都需要重复进行的数据处理工作不仅耗费人力,更影响设计工作的效率。本系统能够依据国军标中规定的数据处理规范,对发动机静止试车数据进行标准处理,计算得出表征发动机内弹道性能的各参数,从而为下一步与设计预期进行对比分析,改进设计提供数据支持。

在获取试验发动机内弹道性能参数之后,需要将其与设计目标进行对比分析,根据初样发动机在试制、试验中出现的问题,修改设计并重新进行有关分析计算和试验。从经济角度考虑,大型发动机应尽可能减少全尺寸发动机的点火试验次数,适当多做一些单项试验,采用计算机进行发动机工作过程仿真。但精细仿真如RocStar等往往需耗费巨大的计算资源和较长的计算时间。鉴于性能偏差或故障原因的最终确定对分析过程高度依赖,分析过程的较长耗时无疑会导致下一步设计修改环节的停滞,影响项目总体进度,为项目的顺利进行引入了不确定因素。因此,如何初步发现问题可能存在的位置,提出合理假设并进行精细分析对设计研发工作的顺利进行有着重要意义。

图 8-1-1含噪声的发动机推力-时间曲线

固体火箭发动机结构简单,部组件数量极少,具有较高的可靠性,其故障模式也不甚复杂。同时,由于发动机的工作过程是推进剂燃烧、燃气膨胀加速并由喷管排出的过程,故而发动机的主要性能偏差和故障绝大部分与装药、燃烧以及流动过程密切相关,如装药的包覆和裂纹、燃面局部变化、侵蚀燃烧、燃速异常以及燃烧不稳定等。

从装药和燃烧的角度分析,推进剂燃速特性与预期值之间的偏差是导致偏差的重要原因之一,这主要体现在实际发动机试车条件下推进剂燃速公式中的重要参数如燃速系数a和压强指数n与在燃烧器、标准发动机以及实际发动机中往往不同,而这些参数在发动机内弹道性能预示中扮演着极为重要的角色,且最终预测结果对其敏感性较强。如果通过参数辨识方法能够获得实际发动机工作条件下的燃速性能参数,并在这些参数条件下进行性能预示,能够与测试数据获得更好的一致性,则一方面可能完全确定了偏差的起因,另一方面可能有助于在因素影响显著性分析中排除一些非显著因素。参数辨识功能将是该系统的一项重要功能,对数据分析工作的开展有重要指导意义。

从另一个角度分析发动机的工作过程,燃烧室平衡压强的维持直接依赖于推进剂燃烧——燃气的生成和加速排出两个过程。内弹道的异常与这两个过程紧密联系。如果通过一些分析方法,能够根据现有试验数据,反推得出与实际试验内弹道曲线相匹配的燃面-时间数据,并与设计的理论燃面-时间曲线进行对比,定位异常发生时刻的燃面情况,将有利于问题的精确定位。因此,燃面的反算是异常诊断系统的核心功能。

总体而言,需要一个系统的固体火箭发动机内弹道性能分析与诊断环境,对发动机实际试车数据进行标准化、精细化处理,依据相关标准处理规范得出发动机各项性能参数并生成符合标准的报告文件;同时,针对试车结果与设计预期之间的差距,深入挖掘数据所隐含的信息,通过参数辨识从参数角度探寻差异产生的原因,利用燃面反推获取与发动机实际工作过程相对应的燃面-时间变化过程,定量分析判断异常原因,定位异常位置,为异常原因的精确定位提供量化依据。同时,该系统应提供全面、灵活的功能接口,方便用户调用。

2. 算法原理

2.1滤波原理

2.1.1研究背景

固体火箭发动机点火过程是一个十分复杂的物理化学过程,包括点火器启动,点火药各组分间的化学反应并产生高温燃气,点火器产生的高温燃气向推进剂表面的热传导、对流传热和热辐射,推进剂表面点燃,火焰在推进剂表面的传播,燃烧室内流场的建立,燃烧室压力增大到准稳态工作状态等。点火过渡过程包含一系列复杂的相互关联的过程或事件,诸如:点火器启动信号;点火器各组分间化学反应引起的热生成;从点火器产物向推进剂表面的热传导、对流传热和热辐射;火焰在整个推进剂表面的扩展;燃烧室内流场的建立;燃烧室压力增大到准稳态工作状态。有时还伴随有许多不规则现象,诸如超压,滞火(被延迟的点火),破坏性冲击波(爆轰),燃烧震荡,间歇和碎灭。

固体火箭发动机的点火是一个复杂的、瞬态的物理-化学变化过程,研究人员对固体推进剂燃烧机理、发动机点火过程进行了大量的理论和试验研究,并建立了各种点火模型,为固体火箭发动机点火装置设计提供了基础。但在工程设计中,由于影响点火性能的因素及其影响程度各不相同,因此已建立的各种模型只能定性地说明各种影响因素对点火性能的影响趋势。研制经验表明,一种新型点火装置的设计,往往要通过发动机点火试验才能确定,在许多情况下可能要经过若干次反复,这便存在着研制费用高、研制周期长的问题。

由于点火过程的复杂性,以及测试水平的局限性,目前点火过程的试验数据处理是个比较关键的问题。首先点火器的发火系统的启动将对整个发动机燃烧室产生初始的压强冲击,另外由于点火过程的特殊性,这一过程中敏感数据较多。但是试验时间很短,由试验数据提取有关点火启动的各项压力数据具有很大的工程应用价值,这将为理论分析提供强有力的数据支持。

点火装置是固体火箭发动机的一个重要部件,其可靠性直接影响到固体火箭发动机工作的可靠性,而且它含有爆炸性危险敏感部件,经常处于各种复杂的环境条件之下,要在地面环境条件下运输、贮存和勤务处理;又要经受导弹或卫星发射时的振动、冲击和加速度的影响;甚至还要经受外层空间的电子和质子辐射的影响。因此必须对点火装置进行各种模拟试验,以检验其安全性、可靠性和适应性,以满足固体发动机总体技术要求。

另外,发动机点火过程中参数变化剧烈,加之各种干扰的存在,试验所采集到的数据中往往在在噪声的作用下呈现出明显的振荡,偏离其所反映的真实过程特性,为后续分析工作引入了不小的误差,也为一些对数据形态有要求的算法的实施设置了障碍。因此,有必要对数据进行高效的预处理,在保证信号不失真的前提下最大限度地滤除噪声干扰。

点火器或是点火发动机通常在模拟发动机内做点火试验,检验其设计的正确性和合理性。本模块利用内弹道曲线数据,通过对其进行预处理,滤除噪声干扰,并依据国军标和所标的相关规定,结合专业理论,处理出点火阶段的点火延迟时间、最大压强值、点火时间和点火上升速率。

鉴于数据预处理的通用性,以及噪声干扰在发动机试验数据中存在的普遍性,本系统将点火过程试验数据处理功能进行了划分,分为数据预处理部分和点火特征量计算两部分。其中数据预处理作为一个单独模块,经过其处理过的数据可用于通用试验数据处理和典型发动机试验数据处理;点火特征量的计算集成到了典型发动机试验数据处理模块中,最终处理结果连同各参数值一并显示到结果列表中,并在生成报表时插入其中。因此,本部分文档主要对数据预处理模块及其使用进行详细介绍。

试验数据本质上是反映试验研究对象属性和特征的物理量随时间变化的离散时间数字信号。数据处理的本质是信号处理。从信号与系统的角度出发,发动机的设计与研究过程实质上是对系统进行建模,并通过试验对比检验模型和信号之间的差异,不断完善数学模型,改进设计模型与方法,提高设计水平。具体而言,在建立数学模型之后,便能求解出在预定工作条件下的发动机内弹道性能,即系统模型在典型信号激励下的输出响应,并将相同工作条件下的试验数据与预示性能也即系统预期响应与相同激励下的实际系统输出响应作比较。若两个响应之间的吻合程度在允许范围之内,则可以用该模型对系统进行分析设计。否则,需要进一步对模型进行改进。真实系统的数学分析过程如图1-1-2所示。

图 8-1-2真实系统数学分析过程图

实际采集到的试验数据由有用信息和噪声干扰两部分组成,其中信息部分是描述试验发动机各种物理特性的物理量变化过程的信号,噪声则是干扰信号传输或测量系统中不应存在却无时无刻无处不在的信号。噪声自身携带了有关噪声源和其所传播于的环境的信息,使所测信号产生一定程度的失真。噪声和失真的种类和来源多种多样,主要包括如下几个大类:

  • 电子噪声,如热噪声与散弹噪声;
  • 由机械运动、振动和碰撞引起的声学噪声,如旋转机械引发的噪声;
  • 干扰电磁信号的传播与接收的电磁噪声;
  • 由于电压的存在而引发的的静电噪声;
  • 通信信道失真与衰减;
  • 数模转换噪声以及由于测试网络阻塞造成的数据包丢失。

就一般的发动机试车试验而言,噪声和失真的主要来源包括上述列出的各个类别。失真使试验数据偏离真实值,无法反映发动机实际工作情况,为后续的性能分析引入了不必要的误差。另外,由于物理过程本身的连续性,真实数据应该具有一定的光滑性,但噪声的引入和失真的产生对数据的形态产生不利影响,如导致局部奇异值和毛刺的出现,影响梯度依赖算法的实施等。同时,数据采集的起始和结束时刻分别提前和滞后于试验过程,最终数据中含有一定量的冗余数据,从中提取有用的数据将减少后续数据分析处理过程的运算量,节省计算时间,提高程序运行效率。因此,有必要在进行性能分析、参数辨识和异常诊断等一系列后续数据操作处理之前,对试验数据进行一系列预处理,包括数据提取和噪声消除。图3为数据预处理模块的工作流程。

图 3数据预处理模块工作流程图

2.1.2信号变

数字信号处理的基本方法主要包括四种,分别是:基于变换的信号处理、基于源-滤波器模型的信号处理、基于贝叶斯统计模型的信号处理和基于神经网络的信号处理。基于变换的信号处理利用一组基础信号(如正弦信号、特征向量和小波)对系统或信号进行表达,以方便对信号进行进一步的分析、诠释和操作。变换类信号处理方法主要包括傅里叶变换、拉普拉斯变换、z变换和小波变换,其中傅里叶变换应用最为广泛,小波因其具有多分辨率的特点近年来得以快速发展,应用的领域和范围也迅速扩张。基于源-滤波器模型的信号处理方法能够利用更多的信息,但对偏离模型特征的信号偏差较为敏感。基于贝叶斯统计模型的信号处理方法主要处理随机过程。神经网络信号处理方法在信号空间的非线性划分、特征提取和模式识别以及决策系统中具有明显优势。

发动机试车数据是离散时间信号,即信号是时间的函数。通常对数据的定性描述和定量分析主要针对时域,这些处理方法直观、物理概念强,易于理解。除时域以外,信号在其他域上也有相应的物理意义。域变换能够从另一个角度诠释信号,使时域中信号的隐藏特征在变换域中更为显著,从而有利于进一步的分析和操作。常见的变换域有频域、尺度域、时频域和角度域。根据域变换揭示出的信号特征对信号进行进一步分析处理能够最大程度利用信号所承载的信息,避免直接处理方法带来的失真以及其他不良影响。本系统主要采用基于变换的数据处理方法。以下分别讨论两种应用广泛且极具代表性的域变换方法——傅里叶变换和小波变换。

(1)傅里叶变换

对于周期信号,可用傅里叶级数展开的方法将时域信号变换为频域信号,以幅值谱、相位谱和功率谱来表示信号的频域特征。对于非周期信号,采用傅里叶变换方法,变换后的信号相应地用幅值谱密度、相位谱密度以及功率谱密度来表示。傅里叶变换的具体数学定义如下:

设\(x\left( t \right)\)是一个连续时间信号,如果\(x\left( t \right)\)属于L2空间,即

\(\int_{-\infty }^{\infty }{{{\left| x\left( t \right) \right|}^{2}}dt}<\infty \)

x(t)的傅里叶变换(FT)存在,并定义为

\(FT\left[ x\left( t \right) \right]=X\left( i\Omega  \right)=\int_{-\infty }^{\infty }{x\left( t \right)}{{e}^{-i\Omega t}}dt\)

其逆变换为

\(F{{T}^{-1}}\left[ X\left( i\Omega  \right) \right]=x\left( t \right)=\frac{1}{2\pi }\int_{-\infty }^{\infty }{X\left( i\Omega  \right)}{{e}^{i\Omega t}}d\Omega \)

式中\(\Omega =2\pi f\),单位为rad/s;f为频率,单位为HZ。将\(X\left( i\Omega  \right)\)表示为\(\left| X\left( i\Omega  \right) \right|{{e}^{i\varphi \left( \Omega  \right)}}\)的形式,可得到\(\left| X\left( i\Omega  \right) \right|\)和\(\varphi \left( \Omega  \right)\)随Ω变化的曲线,分别称为x(t)的幅频特性和相频特性曲线。由此,通过式(2)和(3),傅里叶变换完成了将信号从时域到频域的变换。

对离散的数字信号,需要借助离散傅里叶变换(DFT)将信号变换到频域。对N点序列x(n),其DFT变换定义为

\(DFT\left[ x\left( n \right) \right]=X\left( k \right)=\sum\limits_{n=0}^{N-1}{x\left( n \right)W_{N}^{nk}}    k=0,\ 1,\ \cdots ,\ N-1,\quad {{W}_{N}}={{e}^{-i\frac{2\pi }{N}}}\)

相应逆变换为

\(DF{{T}^{-1}}\left[ X\left( k \right) \right]=x\left( n \right)=\frac{1}{N}\sum\limits_{k=0}^{N-1}{X\left( k \right)W_{N}^{-nk}}   n=0,\ 1,\ \cdots ,\ N-1\)

从式可以看出,由于WNx(n)都可能是复数,对N点序列x(n),需要进行次N2复数乘法和N(N-1)次复数加法运算,而一次复数乘法运算相当于四次实数乘法运算,一次复数加法运算相当于两次实数加法运算。显然,当序列长度N增大时,DFT的计算量将以N2的速度增加。实际应用中,多采用快速傅里叶变换(FFT)算法,该算法能有效减少计算次数,显著提高计算效率。

傅里叶变换经过百年来的发展,已逐渐发展成为一项重要的数学分支,以及信号分析和处理领域内的重要工具,但在实际应用过程中,人们发现了傅里叶变换的不足之处,主要表现为三个方面:

  • 傅里叶变换缺乏时间和频率的定位能力,即对于给定信号,难以通过傅里叶变换得出某一特定时刻对应的频率值或反过来找出特定频率分量对应的具体时刻。
  • 傅里叶变换对非平稳信号即频率随时间变化的时变信号的分析存在局限性,而现实中绝大多数信号频率都随时间变化。
  • 傅里叶变换在分辨率上存在不足。

(2)小波变换

以上讨论了傅里叶变换存在的不足,而正是这些不足成为了建立新的信号分析和处理方法的直接推动力,短时傅里叶变换(STFT)、Wigner-Ville变换、Hilbert变换以及80年代后期开始兴起的小波变换等方法陆续出现并实际应用于信号分析和处理领域,相应的多分辨率分析和时频联合分析方法也飞速发展。近10多年来,小波分析在量子物理、信号及图像处理、语音处理、数值计算、模式识别、故障诊断等学科和领域得到了广泛应用,并取得了具有科学意义和应用价值的丰硕成果,被认为是近年来在工具和方法上的重大突破,为非平稳信号分析提供了便利。

小波变换是一种时频分析方法,并具有多尺度的特点,它的窗口面积固定,但时间窗和频率窗都可以改变。因此,小波变换在时域和频域都具有表征信号局部特征的能力,以及多分辨率的特点。具体而言,小波变换在低频部分具有较低的时间分辨率和较高的频率分辨率,在高频部分具有较高的时间分辨率和较低的频率分辨率,这一特性是它与傅里叶变换相比,更适合于非平稳信号的处理、分析以及特征提取。连续小波变换和离散小波变换以及各自相应逆变换的定义分别如下。

对于平方可积函数\(\psi \left( t \right)\in {{L}^{2}}\left( R \right)\),如果满足小波的容许条件

\({{C}_{\psi }}=\int_{-\infty }^{+\infty }{{{\left| \hat{\psi }\left( \omega  \right) \right|}^{2}}{{\left| \omega  \right|}^{-1}}d\omega }<0\)

则称其为一个基本小波或母小波(Mother Wavelet)。对母小波进行伸缩和平移后得到

\({{\psi }_{a,b}}\left( t \right)=\frac{1}{\sqrt{a}}\psi \left( \frac{t-b}{a} \right)  a>0  a,b\in R\)

式中a被称为尺度因子,b被称为平移因子,\({{\psi }_{a,b}}\left( t \right)\)称为小波基函数或简称为小波基。给定平方可积信号x(t),即\(x\left( t \right)\in {{L}^{2}}\left( R \right)\),则x(t)的小波变换(Wavelet

Transform, WT)定义为

\(W{{T}_{x}}\left( a,b \right)=\frac{1}{\sqrt{a}}\int{x\left( t \right){{\psi }^{*}}\left( \frac{t-b}{a} \right)}dt=\int_{-\infty }^{+\infty }{x\left( t \right)\psi _{a,b}^{*}\left( t \right)dt}=\left\langle x\left( t \right),{{\psi }_{a,b}}\left( t \right) \right\rangle \)

式中a、b和t都是连续变量,因此该变换又名连续小波变换(CWT)。母小波可以是实函数,也可以是复函数。若是实信号,\(\psi \left( t \right)\)是实函数,\(W{{T}_{x}}\left( a,b \right)\)也是实函数,反之\(W{{T}_{x}}\left( a,b \right)\)为复函数。连续小波变换的逆变换为

\(x\left( t \right)=W{{T}^{-1}}\int_{-\infty }^{+\infty }{\int_{-\infty }^{+\infty }{W{{T}_{x}}\left( a,b \right)\tilde{\psi }\left( \frac{t-b}{a} \right){{a}^{-2}}dadb}}\)

其中,\(\tilde{\psi }\left( t \right)\)是小波函数的对偶。对于正交小波来说,其对偶就是它本身。

以上讨论了连续小波变换,该变换针对连续时间信号。为了对发动机试车数据这样的离散时间信号进行处理,同时鉴于连续小波变换是冗余的变换,不仅计算量大,所占的存储空间也大,故而采用离散小波变换,该变换可以在不损失信息的前提下,大大减少计算量和存储空间。通常,离散小波变换针对尺度即式(8)中的 和 进行二进离散化处理,取a=2jj为整数,如果\(\psi \left( t \right)\)满足下列容许条件

\(\sum\limits_{j\in Z}{{{\left| \psi \left( {{2}^{j}}\omega  \right) \right|}^{2}}=1}\)

则称\(\psi \left( t \right)\)为基本二进小波,其对应小波基

\({{\psi }_{{{2}^{j}}}}\left( t-x \right)={{2}^{-{j}/{2}\;}}\psi \left( \frac{t-x}{{{2}^{j}}} \right)\)

称为二进小波。

二进小波变换定义为

\({{W}_{{{2}^{j}}}}x\left( b \right)=\left\langle x\left( t \right),{{\psi }_{{{2}^{j}}}}\left( t \right) \right\rangle ={{2}^{-{j}/{2}\;}}\int_{-\infty }^{+\infty }{x\left( t \right)}\overline{\psi \left( \frac{t-b}{{{2}^{j}}} \right)}dt\)

二进小波有如下重构公式

\(x\left( t \right)=\sum\limits_{j\in Z}{{{2}^{-j}}\int_{-\infty }^{+\infty }{{{W}_{{{2}^{j}}}}x\left( b \right)\overline{{{\psi }_{{{2}^{j}}}}\left( b-t \right)}}db}\)

(1)问题描述与方法讨论

设长度为N的信号xn被噪声en所污染,则所测得的实际含噪声信号为

\({{X}_{n}}={{x}_{n}}+{{e}_{n}}\)

因此,去噪的目标是通过一系列信号处理方法,得到信号xn的一个逼近信号\({{\hat{x}}_{n}}\),并使其在某种误差估计下是信号xn的最优逼近,既去除噪声又保持边缘和局部特性。去噪的基本思想是利用一些途径,使信号和噪声之间的差别更加明显,从而将两者从所测信号中有效分离出来,保留真实信号,去除噪声信号,从而达到去噪的目的。日常出现的噪声大都服从或近似服从高斯分布,因此后续去噪效果对比分析中采用高斯噪声为主。

现有的去噪方法主要有移动平均、基于傅里叶变换的去噪和小波去噪。移动平均是一种特定的低通滤波器,可去除信号中的高频噪声,该方法会导致低频部分的幅值存在一定程度上的衰减,衰减程度随平均长度的增加而增大。基于傅里叶变换的去噪方法是一种较为传统的方法,它相对于移动平均法,可灵活调整滤波器的频响,在去除高频噪声的同时又尽可能保持低频成分的幅度,具有可控的过渡带和低频衰减特性。但总体而言,移动平均和基于傅里叶变换的低通滤波方法不适用于宽频、时变、多尺度信号的去噪。相比之下,小波去噪方法能够有效地去除高频噪声,低通并保留高频中幅度较大的部分,对非平稳信号去噪效果良好。鉴于傅里叶变换去噪方法较为经典且应用广泛,因此在介绍小波去噪方法的同时叙述傅里叶变换去噪的基本原理和方法,并在后续的去噪实例分析中对两种方法的去噪能力和效果进行对比分析。

2.1.3傅里叶变换去噪

对于时变信号的分析,通常在时域和频域分别展开,时域描述信号随时间的强度变化,频域描述信号各个组成频率的幅值分布。信号变化率大的部分,对应于高频分量,信号变化率小的部分则相应于低频分量。傅里叶变换去噪方法在这一前提下通过频域变换,对频域中的高频信号成分进行抑制,对低频成分进行保留,从而去除噪声,尽可能还原真实的信号。

低通滤波的频域表达式为

\(G\left( u,v \right)=H\left( u,v \right)\cdot F\left( u,v \right)\)

式中,\(G\left( u,v \right)\)是去噪后信号的傅里叶变换,\(H\left( u,v \right)\)是滤波器的传递函数,\(F\left( u,v \right)\)是含噪声信号的傅里叶变换。在传递函数\(H\left( u,v \right)\)的作用下,原含噪声信号的高频分量被衰减,平稳的低频分量得以保留,从而得到最终的去噪后信号\(G\left( u,v \right)\),将其进行傅里叶逆变换即可重构出降噪后的时域信号。低通滤波器有多种,其中巴特沃斯(Butterworth)滤波器和切比雪夫(Chebyshev)滤波器最为常用。低通滤波器在去噪过程中发挥了比较重要的作用,其特性对去噪效果的影响也较为明显:若截止频率太低,则滤波后仍有大量噪声残留;反之,若截止频率太高,则一部分真实信号将被当做噪声滤除。

作为经典的信号去噪方法,傅里叶变换去噪建立在真实信号频谱和噪声信号频谱没有重叠的前提下。在实际情况中两者的频谱往往有一定的重叠,重叠范围的大小和重叠程度的强弱都与信号自身特点密切相关。如前所述,傅里叶变换缺乏时间和频率的定位能力,无法给出信号在具体时刻上的频率特征,使得信号在时域中的任何突变都会影响到整个频域,主要适用于平稳周期性信号。因此,在傅里叶变换的基础上进行去噪,难以在不同时刻、不同尺度上对信号和噪声进行有效区分,从而最终难以获得理想的去噪效果。

巴特沃斯滤波器是较简单而常用的模拟低通滤波器,具有在通带内最大平坦,在正频率范围内随频率升高而单调下降的幅频特性,阶次越高特性曲线越接近矩形等特点,又被称为通带最平坦幅度特性滤波器。切比雪夫滤波器采用使通带内或阻带内幅度呈等幅波动的等波纹设计来逼近理想特性,可以提高选择性,其系统设计的方法与巴特沃斯型相同,在相同指标情况下,切比雪夫滤波器的阶数比巴特沃斯滤波器要低,但相频特性要差,且设计相对要复杂。

巴特沃斯滤波器和切比雪夫滤波器参数过程设置基本相同。首先进行FFT频谱分析,根据频率特性选择滤波器种类(低通、高通、带通、带阻),低通和高通只需设置低截止频率,带通和带阻需要设置低截止频率和高截止频率,再选定滤波器阶数即可完成滤波器的设置。根据实际经验和频谱分析,需要的压强数据集中在低频信号区,因此选择低通滤波器,具体的低截止频率则应根据实际频谱图设置。

2.1.4卡尔曼滤波器

卡尔曼滤波器用于估计离散时间过程的状态变量\(x\in {{R}^{n}}\),由以下离散随机差分方程描述:

\({{x}_{k}}=A{{x}_{k-1}}+B{{u}_{k-1}}+{{\omega }_{k-1}}\)                                                        (16)

定义观测变量\(z\in {{R}^{m}}\),得到量测方程:

\({{z}_{k}}=H{{x}_{k}}+{{\upsilon }_{k}}\)     (17)

式(17)中矩阵H表示状态变量xk对测量变量zk的增益。实际中H可能随时间变化,但在这儿假设为常数。实际系统中,过程激励噪声协方差矩阵Q,观测噪声协方差矩阵R也可能会随每次迭代计算而变化,在这儿我们也假设它们是常数。

定义\(\hat{x}_{k}^{-}\in {{R}^{n}}\)(-代表先验,^代表估计)为在已知第k步以前状态情况下第k步的先验状态估计。定义\(\hat{x}_{k}^{{}}\in {{R}^{n}}\)为已知测量变量zkk步的后验状态估计。

先验估计\(\hat{x}_{k}^{-}\)和加权的测量变量zk及其预测\(H\hat{x}_{k}^{-}\)之差的线性组合构成了后验状态估计\(\hat{x}_{k}^{{}}\),即卡尔曼滤波器表达式:

\({{\hat{x}}_{k}}=\hat{x}_{k}^{-}+K({{z}_{k}}-H\hat{x}_{k}^{-})\)                               (18)

式(18)中K为卡尔曼增益,K的一种表示形式为:

\({{K}_{k}}=P_{k}^{-}{{H}^{T}}{{(HP_{k}^{-}{{H}^{T}}+R)}^{-1}}\)

\(=\frac{P_{k}^{-}{{H}^{T}}}{HP_{k}^{-}{{H}^{T}}+R} \)                                   (19)

由(19)式可知,观测噪声协方差R越小,K越大。另一方面,先验估计误差协方差\(P_{k}^{-}\)越小,K越小。K越大,则测量变量zk的权重越大。

卡尔曼滤波器可分为两个部分:时间更新方程和测量更新方程。时间更新方程负责及时向前推算当前状态变量和误差协方差估计的值,以便为下一个时间状态构造先验估计。测量更新方程负责反馈——也就是说,它将先验估计和新的测量变量结合以构造改进的后验估计。

离散卡尔曼滤波器时间更新方程:

\(\hat{x}_{k}^{-}=A{{\hat{x}}_{k-1}}+B{{u}_{k-1}}\)                                                                 (20)

\(P_{k}^{-}=A{{P}_{k-1}}+Q{{A}^{T}}+Q\)       (21)

离散卡尔曼滤波器状态更新方程:

\({{K}_{k}}=P_{k}^{-}{{H}^{T}}{{(HP_{k}^{-}{{H}^{T}}+R)}^{-1}}\)                                                       (22)

\({{\hat{x}}_{k}}=\hat{x}_{k}^{-}+{{K}_{k}}({{z}_{k}}-H\hat{x}_{k}^{-})\)                               (23)

\({{P}_{k}}=(1-{{K}_{k}}H)P_{k}^{-}\)   (24)

卡尔曼滤波器整个计算过程不断递归推算,完成一轮计算后,将上一次的计算结果带入下一次的计算式中,不断向前推算。图4显示了滤波器的整个操作流程。

图 4卡尔曼滤波流程示意图

对于压强传感器,过程的状态不随时间变化,所以A=1;没有控制输入,所以u=0;包含噪声的观测值是状态变量的直接体现,所以H=1。因此,式(16)和式(17)改写为:

\({{x}_{k}}={{x}_{k-1}}+{{\omega }_{k-1}}\)                                                                 (25)

\({{z}_{k}}={{x}_{k}}+{{\upsilon }_{k}}\)      (26)

则滤波器的时间更新方程和测量更新方程分别为:

时间更新方程:

\(\hat{x}_{k}^{-}={{\hat{x}}_{k-1}}\)      (27)

\(P_{k}^{-}={{P}_{k-1}}+Q\)             (28)

测量更新方程:

\( {{K}_{k}}=P_{k}^{-}{{(P_{k}^{-}+R)}^{-1}} \)

\( =\frac{P_{k}^{-}}{P_{k}^{-}+R} \)

(29)

\({{\hat{x}}_{k}}=\hat{x}_{k}^{-}+{{K}_{k}}({{z}_{k}}-\hat{x}_{k}^{-})\)                                (30)

\({{P}_{k}}=(1-{{K}_{k}})P_{k}^{-}\)     (31)

该滤波器需要输入两个初始值\({{\hat{x}}_{0}}\)和P0,两个系数RQ。在理想情况下,压强传感器没有信号输入时理论值为0,因此令\({{\hat{x}}_{0}}=0\)。P0的选择并不关键,几乎任何\({{P}_{0}}\ne 0\)都会使滤波器最终收敛,可先令P0=1。测量噪声方差R可以通过观测得到,是滤波器的已知条件。在LabVIEW的程序实现中通过截取点火前一段压强数据,在控件选板中选择数学→概率与统计→标准差和方差控件,使用该控件可以很方便地计算出一段压强数据的方差,即可作为R值。由于无法直接观测到过程信号xk,所以难以确定过程激励噪声协方差Q的值,可通过Q的选择来产生可以接受的结果,一般Q取一个小的非零常数可以方便地调整滤波器参数(如Q=10-5)。无论是否有一个合理的标准来选择系数,通常都可以通过调整滤波器的系数来获得更好的性能。初步选定系数值后再进行调整直到获得一个理想的结果。

2.1.5小波去噪

用传统的傅里叶分析方法进行去噪时具有很多的条件限制,只是通过一个低通或带通滤波器,对于过程信号并不平稳或含有宽带噪声的信号并不适用,同时傅里叶分析去除噪声不能给出该信号在某一时刻的突然变化的情况,分析的过程只是完全在频域里进行信号分析。而小波变换是一种时频局部化分析方法,优点在于:分辨率多,窗口面积不变、时域窗口和频窗口都可以改变,对于去除非稳态信号的噪声,不仅不会损坏有用信号的突变部分,而且又能很好地去除噪声。目前流行地一维小波去噪的方法主要有三大类:阈值法、模极大值法、平移不变量法。

小波系数与信号是对应的,其中包含有要获得信号的重要信息,信号的小波系数特点是:数量较少,但幅值一般比较大,而与噪声对应的小波系数的特点是:数量较多,分布一致,但是幅值一般比较小。小波变换阈值去噪法的阈值选取有软阈值、硬阈值可以选择,这种方法已经相对全面成熟了,小波变换阈值去噪法也是最经典的小波去噪方法,主要步骤分为三步:

  • 首先进行含噪声信号的小波变换。选择合适的小波分解层数k和合适的小波基,然后进行小波分解,分解带有噪声的信号,得到与之对应的小波分解系数。
  • 其次阈值量化分解后的高频系数。对于从1到k的每一层,选择一个适当的阈值和合适的阈值函数,将分解得到的高频系数进行阈值量化,得到估计小波系数。
  • 进行小波反变换。根据小波分解后的第k层的低频系数,即尺度系数和经过阈值量化处理的各层高频系数,即小波系数。最后利用重构算法进行小波重构,得到去噪后信号。

小波变换在各尺度上的传播特性对于信号和噪声来说是不同,根据这种不同,噪声产生的模极大值点将被剔除掉,信号对应的模极大值点会被保留,最后利用剩下的模极大值点重新生成小波系数,信号由此将被恢复。

在小波变换阈值去噪法基础上改进得到更优的平移不变量小波去噪法。虽然阈值法已经可以取得很好的去除噪声的效果,但在有些信号不连续的邻域内,阈值法去噪会表现出极为显著的非正常信号,比如在不连续的信号点的附近可能会产生在一个特定的数值水平上下不断变化的现象,即伪吉布斯现象。使用平移不变量去噪的方法,可以有效地抑制伪吉布斯现象。其原理是:“平移-去噪-平均”,即将有噪声的信号循环平移n次,再用阈值法对循环平移后变化的信号进行去噪,最后再对去噪的结果求平均值,这就是所谓的小波平移不变量去噪法。

小波变换阈值去噪法具有能近似估计原始信号最优值、若用软阈值的方法去噪则能够使被估计信号最小化最大均方误差、得到的估计信号肯定会比原始信号光滑且不会产生附加振荡、具有广泛的适应性以及计算速度快等优点,是小波去噪方法中应用最广泛的一种。一般情况下,均可选用该方法去噪。

现有去噪方法多具有低通性,因此在对噪声进行分离抑制的同时,也平滑了非平稳信号的突变点。高通滤波器虽然可以使边缘更加突出,但背景噪声也同时被放大。由此可以看出,基于傅里叶变换的去噪方法存在着保护信号局部性和抑制噪声之间的矛盾。相比傅里叶变换去噪,小波去噪能力更加强大,更适合于对非平稳信号进行去噪,这主要是由小波的以下重要特性所决定的:

  • 低熵性。小波系数的稀疏分布,使得信号变换后的熵值降低;
  • 多分辨率特性。由于采用了多分辨率的方法,可以良好地刻画信号的非平稳特征,如边缘、尖峰、断点等,以便于提取特征以及保护真实信号;
  • 去相关性。因为小波变换可以对信号进行去相关,且噪声在变换后有白化趋势,所以在小波域相比在时域更利于去噪;
  • 小波基选择的多样性。由于小波变换可以灵活地选择变换基,所以可以针对不同应用场合选用不同的小波函数,以取得最佳的处理效果。

小波分解与重构算法本质上相当于一个多通道的带通滤波器,适用于信号和噪声频带相互分离的确定性噪声,该方法复杂程度低,计算量较小。小波阈值收缩法具有能够得到原始信号最有近似估计、计算速度快和适应性广的特点,是小波去噪方法中应用最为广泛的一种。平移不变量小波法是在阈值法基础上的改进,可以有效去除信号不练续处所产生的伪吉布斯现象,但在计算速度上存在不足。小波变换模极大值法在信号中含有较多奇异点时去噪性能较为优越,但计算量太大,速度太慢。综上分析,本系统采用阈值收缩法作为去噪算法。

小波阈值收缩去噪法利用信号经小波变换后信号的能量主要集中于较大的小波系数,噪声的能量仍分散于整个小波域内的特性,通过设定阈值,保留信号部分,去除噪声部分,然后对其做逆小波变换重构信号,从而实现信号去噪。该算法的具体实施步骤为:

\(Y=W(X)\)                         (32)

\(Z=D(Y,T)\)                           (33)

\(\hat{X}={{W}^{-1}}(Z)\)              (34)

式中,W()和W-1()分别表示小波分解和重构算法,D(Y,T)表示给定阈值T的小波分解系数的收缩操作,X为含噪信号,Y为X的小波分解系数,Z表示经过收缩操作后的小波系数,\(\hat{X}\)为去噪后的信号,该去噪算法的实施流程如图5所示。

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注

You cannot copy content of this page