基于图像处理的航空发动机自适应频谱报警技术研究

2019-09-25

计算机测量与控制 2019年9期

(1.航空工业西飞西飞设计院，西安 710089; 2.中国飞行试验研究院测试技术研究所，陕西西安 710089)

0 引言

航空发动机结构复杂、工作环境多变且恶劣，经常出现各种故障，美国航空运输协会调查了上世纪80年代至21世纪初期间的航空事故，发现最常见的机械失效就是起落架和航空发动机故障[1-3]。而航空发动机所带来的故障往往直接影响整套飞行系统的运行，从而造成巨大经济损失，甚至导致机毁人亡的严重后果。因此，对这些航空发动机的运行状态进行监测，对可能发生的故障给予报警并准确确定故障原因是十分必要的。而振动信号能够迅速直接地反映出发动机的运行状态，不需要对发动机停机或者拆卸，便于实现在线和离线监测，同时振动信号便于测量，振动信号的分析理论丰富而完备，包括时域分析方法、频谱分析方法、时频分析方法、图形分析法等，因此广泛地应用于旋转机械、往复机械、流体机械、转轴、轴承、齿轮等的状态监测[4-5]。

报警技术是振动监测系统的核心技术之一，报警信号是对监测对象状态的直接评估，通过划分不同报警等级，可表征当前状态的严重程度。当监测对象在设定条件下发生报警时，意味着发动机运行状态恶化并且可能发生某种故障。尽管报警信息在多数情况下不能直接指出故障源，只是指示一种现象，但维修人员可根据所指现象进行分析，缩小检测范围，进行有目的地检查或维修，从而防止故障的发生，并降低维修成本[4-5]。因此从预防事故的角度来讲，发动机运行状态的报警比故障诊断更为重要和有效，状态监测的误报、漏报等情况很大程度上与报警设置有直接关系，报警的准确性一直是故障诊断可靠性中的一个关键问题。可靠的报警技术不仅可以确保航空发动机运行的安全，还可以降低发动机状态监测诊断的成本，对安全生产和国民经济都有重要意义[6-9]。

目前，常用的报警技术有参数化报警方法[10-11]和频谱报警方法[12-14]两种。其中，参数化报警方法利用历史数据得到的统计指标作为报警阈值，并通过制定各种振动标准进行报警。而频谱报警则是通过判断各频率成分的幅值相位是否发生异常来报警的。运行中的发动机在各种内外因素综合影响下，状态处于动态迁移之中[15]，对监测对象的各种测量数据也在不断变化之中。而传统的参数化报警技术和频谱报警技术侧重设置固定的报警阈值，对报警对象的变化适应性差，无法满足发动机不同工况下的报警准确性。

针对上述问题，本文通过研究基于图像处理的航空发动机自适应频谱报警方法，将新型包络生成方法和频段自动划分方法应用到均化频谱的报警中，以提高航空发动机故障报警的准确性。

1 频谱报警算法

频谱各点依次连接形成频谱曲线，其幅值构成一组一维数据序列，如果将该序列连续复制后组成二维数据序列，则二维序列在一个方向上是频谱曲线的幅值，另一个方向上所有数据相同。灰度图像是RGB色彩模式下，R(红)、G(绿)、B(蓝)3个颜色通道数据均相等时的图像，每个通道有0～255共256级强度值。将二维序列归一化，归一化范围取[0, 255]，则各点数值对应某灰度级别，由各点灰度级别构成二维图像，便得到频谱曲线对应的灰度图像，图1给出了一个简单示例。

图1 频谱曲线及其对应的灰度图像

在图1中，频谱中各点归一化为0～255，数值越小则能量越低，扩展到二维灰度图像后，高能量区(亮度高)和低能量区(亮度低)区别明显。频谱曲线可以看作二维图像的一种压缩，因此频谱报警可以借鉴图像处理的相关理论。

1.1 基于平滑滤波的包络线生成算法

发动机振动信号的频谱在截止频率范围内的所有频率和对应幅值构成一条曲线C，离散频谱情况下，假设C共有n个点，在频率轴上的坐标为{Fi,i=1,2,…,n}，在幅值轴上的坐标为{Ai,i=1,2,…,n}。考虑到频率轴上的坐标是等间隔的(FFT频谱)，因此在进行包络计算时，不需要关心具体频率数值，可将{Fi}简化为一个计数序列{i,i=1,2,…,n}。包络生成只是对{Ai}进行计算，其结果和{Fi}仍有对应关系。

下面分析频谱曲线某个点Ai和相邻点之间的关系。图2给出了全频带中的一段频谱曲线。

图2 频谱曲线示意图(含5个点)

一个频率点可能含有来自多个振动源的振动特性对应能量，在图2中，Ai点的振动能量较大。由于傅里叶变换时的窗函数会导致能量泄露，其相邻的两个点Ai-1和Ai+1，可能含有一些与Ai相同的振动特性对应能量。当频谱的频率分辨率较高时，Ai-2和Ai+2也与Ai具有类似联系。振动源的振动特性是变化的，振动能量可能会在某几个频率点中重新分配。即使Ai比相邻点幅值偏小，或者差别不大，其能量贡献的意义相同，只是接受相邻点的能量贡献程度不同。某个频率点的包络计算需要考虑相邻频率点的能量贡献。

图像平滑滤波用于去除图像中的噪声，高斯平滑将每个像素的灰度值用其邻域的加权平均值代替，以此为基础设计包络生成算法。算法核心思想是，相邻点按照一定加权比例对中心点做出幅值贡献，比例分配依据高斯正态分布设定。正态分布曲线关于均值中心对称，曲线的高能量区集中(在中心左右3个标准差范围内集中了99.73%的能量)。虽然中心点与各相邻点的幅值并不具有正态分布规律，但相邻点对中心点的幅值贡献比较集中，可认为该贡献比例近似服从正态分布。

1)幅值贡献比例计算。

2)包络线上各点生成。

包络线上各点(包括中心点)的幅值与幅值贡献比例bi中对应的比例系数相乘，得到幅值贡献分配，显然中心点的幅值贡献是自身(bj+1=1)。包络线E上的点由下式计算确定，

(1)

对于全频带频谱，其频率下限是0 Hz，上限是截止频率F。0 Hz以下频率点没有任何实际意义，可认为其对应幅值为零；而F以上的频率点一般是信号中的高频干扰成分，幅值贡献很小，也可认为其对应幅值为零。在计算接近0 Hz或者F的频率点包络时，中心点的某一侧可能在有效频率范围之外，根据上面分析，通过补充零幅值的点来解决。按照所选定的扩展点数j，在频谱两端各加入j个零幅值点，根据式(1)代入{Ai}计算便得到完整的包络线E。这种处理方法使程序编写十分方便。控制包络线细化程度的参数只有j，能够方便调节。

包络报警只需要将待检频谱和新型包络线比较即可，根据超出包络部分的能量累计值进行报警。由于在包络生成时，考虑了相邻点的能量贡献，因此获得的均化包络敏感度适中，报警准确性高。

1.2 基于边缘检测的频段划分算法

频段报警[16-17]的主要目的是简化频谱监测，各频段单独设定段内总值报警，原监测点从几百个可缩减至几个。频段的自动划分可提高报警信息处理效率。

具有监测价值的频率点及其相邻点对应幅值一般较大，这些点所构成的区域为高能量区，而各高能量区之间是监测意义不大的低能量区。频谱曲线高低能量区的划分与灰度图像中亮度差异较大区域的分割类似，需要找寻分割位置，即图像的边缘。图像边缘检测是指找寻其周围像素灰度发生阶跃变化或屋顶状变化的像素集合。在图3中给出了两种典型边缘的示意，图(b)中的边缘其实可以看作图(a)中边缘的特例，因为其边缘两侧均发生较大阶跃。

图3 两种典型的边缘示意图

边缘检测[18-19]的主要工具是基于各种算子的边缘检测模板。经典边缘算子考察原始图像各像素的领域内灰度阶跃变化，利用邻近一阶或二阶方向导数变化规律检测边缘。最优算子根据信噪比求得检测边缘的最优滤波器，主要算子有LoG算子、Canny算子。

频谱曲线的形状“参差不齐”，对应灰度图像属于高噪声图像。一般的梯度算子不适用于高噪声图像，可借鉴Canny算子的原理进行频谱曲线的边缘检测。Canny算子利用高斯滤波器平滑图像，然后采用一阶偏导数的有限差分计算梯度的幅值和方向，算法沿屋脊带的顶部进行跟踪，最后采用两个阈值来连接边缘。频谱曲线的边缘只须检测一个方向，因此自动划分方法有3个步骤：曲线平滑、梯度计算和边缘搜索。

1.2.1 频谱曲线平滑算法

频谱曲线从形状角度来看，具有很强的噪声，为防止检测时产生大量“虚假边缘”，因此首先需要平滑处理。频段划分要求分割出高能量区(简称H区)和低能量区(简称L区)，两类区域的相邻处存在边缘，因此一小段数据的平滑要求如下：数据整体差异小，则平滑程度低；数据差异大，各点构成阶跃段，则保留该阶跃特性；含有明显的局部极值点，则平滑程度高。从上述要求出发，同时为了提高算法效率，采用平均平滑的方法。

曲线的纵坐标是幅值序列{Ai,i=1,2,…,n}，横坐标是简化的计数序列{i,i=1,2,…,n}(后续计算都是按照序列顺序处理幅值，因此不需要考虑频率值)。平均平滑与前面的包络求取类似，选取当前中心点Ai和两侧各j(j∈N)个对称点，求出所有点的均值作为平滑后的点。

(2)

与包络计算类似，对于曲线两端的平滑采用添加零幅值点的方法保证算法完成。经过平滑处理得到新序列{Si,i=1,2,…,n}。

1.2.2 基于梯度的边缘选取

{Si}在H或者L区内的数值差异较小，在区域交界处具有明显阶跃变化，通过梯度值可找到具有边缘特性的阶跃变化点。对{Si}作前向差分，得到梯度序列{DSi}。

Dsi=Si+1-Si,i=1,2,...,n-1

(3)

由于是前向差分，故{DSi}分别对应{Si}中的前n-1个点。

梯度值较大的点，往往是边缘所在位置。设定梯度阈值ε(ε>0)，如果某个梯度点的绝对值超过阈值，则该梯度点对应幅值点可作为参考边缘点(下文为叙述方便，以梯度点指代原幅值点)，边缘从参考点中按照一定规则确定。从梯度序列DSi求出边缘的过程称为边缘检测。

下面介绍边缘检测中的几个概念。图4是一个梯度序列进行边缘检测的示意图。

图4 梯度序列的边缘搜索

在图4中，实线相连接的各点组成梯度序列，关于横坐标轴对称的两条点划线对应梯度阈值，所有在区间(-ε, +ε)之外的梯度点都可以作为边缘选取的参考点。

从低能量区向高能量区过渡的阶跃变化点称为上升边缘，反之称为下降边缘，一个上升边缘和其右侧(序列递增方向。下文提到的左右方向均指按照序列递减或者递增的方向)紧邻的下降边缘中间是高能量区H，而一个下降边缘和其右侧紧邻的上升边缘中间是低能量区L。

在图4中，U1、U2、U3、U4可以作为上升边缘的参考点；D1虽然超出了阈值，但其对应的原幅值因前向差分关系，会偏左一个横坐标单位，因此下降边缘参考点是D1的右侧紧邻D2，这样处理也是为了和上升边缘偏左特性相适应，使H区的范围适当扩大(对报警监测有利)；对于D3，选取下降边缘参考点则取右侧紧邻U3，但U3又是上升边缘参考点，U3对应原幅值是两个高幅值点之间的一个较低点，一般无须从该点划分频段，因此U3不作为下降边缘参考点。

当梯度阈值过大时，可能无法检测到边缘，为保证结果的一致性，加入极限边缘的概念。全频带中，i=1对应的频率为0 Hz，这个频率处的振动停止，因此幅值为零；在上限截止频率处，待分析的主要频段在全频带内，更高的频率范围可认为振动能量很小，因此上限截止频率点处的幅值也为零。当无法检测到边缘时，这两个特殊零幅值点之间视为L区，0 Hz处有一个极限下降边缘DE，截止频率处有一个极限上升边缘UE。极限边缘保证了边缘检测结果不为空，同时能够有效处理边界边缘问题。在图3中，DE的位置在i=1处，D0的边缘属性根据后续算法调整，可能是DE，也可能演化为上升边缘；UE的位置在i=n处，相当于新添加了一个梯度点，如果i-1对应点成为下降边缘参考点，i=n点也可能演化为下降边缘。

1.2.3 边缘搜索算法

采用搜索方式对梯度序列扫描，依据规则记录边缘点位置。搜索算法的基础是，偏向选取规则和同类型边缘合并规则。从选取规则可知，上升边缘偏左、下降边缘偏右，梯度序列按递增方向搜索上升边缘和按递减方向搜索下降边缘的方法一致，只是下降边缘位置需要右移一位。为了减少搜索中的各种判断，简化算法，采用双向搜索。

算法中需要用到的变量有，梯度序列Ds，上升边缘位置集合Ue，下降边缘位置集合De，梯度阈值e=ε，上升边缘类型标签TagN，下降边缘类型标签TagI。边缘类型标签取值，当Ds≥e，TagN=TagI=+1；当Ds≤-e，TagN=TagI=-1。

1)初始化：

给定Ds共n-1个元素(相比原频谱曲线的点数n)；Ue和De为空；给定梯度阈值e；i=1处对应梯度值不会小于零，在该点只可能检测到上升边缘，因此TagN=-1；i=n-1处梯度值不会大于零，在该点只可能检测到下降边缘，因此TagI=+1。

2)主循环：

Fork=1,k

正向搜索

IfDs(k)≥e&&TagN==-1

Ue中记录位置k；TagN=1;

IfDs(k)≤-e

TagN=-1；

反向搜索

l=n-k;

IfDs(l)≤-e&&TagI==+1

(1)实验所使用的样品有4种，分别为马尾松、北美短叶松、杉木和毛白杨。4种木材均是在180℃恒温条件下进行热处理。但实验结果显示4种不一样的木材在热处理后粗糙度出现了3种变化趋势，其中北美短叶松的粗糙度在1～4 h过程中逐渐变大，即并不是所有的木材都是随热处理时间增加，其粗糙度变小。

De中记录位置l+1；TagI=-1;

IfDs(l)≥e

TagI=+1;

k++;

3)边缘位置重合处理：

反向搜索得到的下降边缘需要进行右移一位处理(l+1)，可能和某个上升边缘重合，因此将Ue和De中的相同元素删除。

根据Se完成频段自动划分后，各频段内的报警对象是段内振动能量，一个简单的能量计算方法是将段内各点振动幅值相加。均化频谱的自动频段划分，使窗滑移过程中的所有报警计算由软件算法实现，因此提高了监测系统的自动化程度。

2 航空发动机自适应频谱报警方法

随着发动机规模增加和监测复杂度的提高，传统报警技术显露出诸多不足，无法满足现代振动监测系统的报警需求。

传统的参数化报警方法是通过制定各种振动标准进行报警，重点是计算历史数据的统计指标，但这些标准的不确定性较强。简单依据发动机种类和操作人员经验设置报警标准，并未考虑发动机运行的实际情况和飞行中的工作状态，这些标准只是从共性基础上提出具有参考意义的判断标准。在飞机飞行中，不同型号的发动机、同类发动机不同的工作负载、发动机运行的历史过程等对发动机的影响都不尽相同，这时根据相同标准对发动机状态进行评估显然是不合理的。报警标准常将发动机当前状态划分为几个等级，但等级无法从物理意义上明确说明发动机处于何种状态，应该对发动机进行什么样的操作。当前的发动机状态报警技术仍然建立在静态报警的基础上，用这些报警标准对发动机状态进行评估，无法适应发动机的个性化。在进行报警的阈值设定和等级划分时，一般都是人工进行报警计算有关的设置。在有效报警的前提下，操作人员需要不断去检查和更新各类设置，随着监测系统的数据处理量的增大，类似的人工环节会大大影响系统的自动化程度，进而影响到飞行的连续性。

另一方面，参数化报警没有考虑发动机的状态变化历程和个体特性，报警阈值来源于各种统计指标，具有静态特点，不能满足系统长期监测的报警准确性。智能化报警要求建立具有自适应特性的报警信息处理模型。运行中的发动机在各种内外因素综合影响下，状态处于动态迁移之中，各种监测对象的测量数据也在不断变化之中。具有自适应特性的报警模型能够跟踪这种变化，通过对当前测量数据的学习更新模型，使报警阈值或报警指标也随之动态调整。具有自适应特性增强了报警的动态性和灵活性，也提高了报警的准确性。频谱图直观反映了不同频率的振动幅值分布，负载波动和故障发生等诸多因素都会引起频谱形状的变化，相应的报警线需要能够在这种变化中判断幅值异常。包络报警中的报警线是对多个频谱进行统计后得到的幅值上限，智能化的包络报警可以对形状变化后的频谱曲线进行敏感度适中的报警监测。全频带的报警如果从振动幅值的总量考虑，不够精细，不能发现局部振动异常，一种较好的办法是将全频带划分为几个频段，每个频段内单独监测振动幅值。这种划分与频谱曲线形状有关，智能化的频谱报警可以自动划分多个频段进行监测，在频谱形状小幅变化后将主要频率分量仍然分割在同一频段中；而在频谱形状大幅变化后可自动重新分割全频带，有效监测各个频段。频段报警是一种精细化的频谱报警技术，加入智能化手段可增强报警的适应性，提高了报警计算的自动化程度。

包络报警和频段报警有利监测细化，应用广泛，前者敏感性强，后者实际意义明确。传统包络生成方法是，取幅值的若干局部极大点相连，或对所有幅值放大一定倍数，具有“点对点”的特点。传统频段报警，可将高能量区划分出来，使监测对象集中，但分割位置选择十分重要。由于负载波动，振动分量的随机性，发动机状态迁移等因素都会引起频谱数据的变化，因此，相应的报警设置可能失效。

综上所述，基于图像处理的航空发动机自适应频谱报警方法具体流程图见图5，具体步骤如下：

1)设运行数据的时刻为[t0,t0+T/a]，取[t0-T,t0]时刻发动机振动信号作为历史数据，将该段信号平均分为a段，得到每段信号长度为T/a的多组发动机历史振动信号，其中通过设置T和a的数值来改变报警的灵敏度。

2)对平均分为a段的发动机历史振动信号分别进行FFT变换，并且将作为报警参考的多组频谱进行均化处理，即求取各频率点对应的幅值均值，得到均化频谱。均化频谱是分析时长内幅值分布的统计结果，可作为频谱报警计算的样本模板。

3)对均化频谱使用基于平滑滤波的包络线生成算法获取均化频谱包络线。

4)对均化频谱使用基于边缘检测的频段划分算法得到均化分割点和均化频段能量值。

5)对运行数据频谱使用基于平滑滤波的包络线生成算法获取运行数据频谱包络线，将运行数据频谱包络线与均化频谱包络线进行差值，得出包络报警结论。

6)按照均化频谱的均化分割点对运行数据频谱进行分割，比较均化频谱与运行数据频谱各个频段能量值，得出频段报警结论。

7)对包络报警和频段报警结论进行综合得到最终报警结论，其中包络报警或频段报警其中任一出现报警时，均视为发出报警。

本方法将作为报警参考的多组频谱进行均化处理，即求取各频率点对应的幅值均值，得到均化频谱。均化频谱是分析时长内幅值分布的统计结果，可作为频谱报警计算的样本模板。对均化频谱，使用包络报警方法获取均化包络线，使用频段报警方法得到均化分割点和均化频段能量值，这些结果反映整个分析时长内的振动特性。做当前时间窗外下一组频谱的报警分析，该组频谱与均化包络线比较得到包络报警结论，按照均化分割点划分该组频谱，将各段能量值与均化频段能量值比较得到频段报警结论。

在发动机运行过程中，具有统计意义的均化频谱不断改变，包络报警和频段报警通过对均化频谱的监测跟踪，实现了自适应的频谱报警。

图5 自适应频谱报警原理

3 实验结果与分析

采用以上的自适应频谱报警方法应用于某型发动机试验中，对其振动情况进行分析，测点编号R1，取2016年8月到2017年1月所采集的振动速度数据，经FFT变换得到40组频谱数据。求取均化频谱，作为报警分析的主要对象。

对均化频谱进行包络计算，由式(1)计算包络线，参数j=3(最高谱峰两侧各约3个点幅值较大)，其包络计算结果如图6所示。

图6 包络计算结果

在图6中，包络线显示为虚线，由于临近点的能量贡献，保证了包络线是频谱曲线的上限。该包络生成过程对每个点依次求解，所以包络细化水平较高。参数j调节包络高度，j越小，包络越低，对频谱形状保留越多。

选取两组频谱，将其和包络线绘制在图7中(对图进行了适当放大)。

图7 两组频谱与包络线比较

在图7中可以看到，包络线对两组频谱具有良好的上限覆盖能力，同时对形状差异较大的位置也能够实现覆盖(例如40 Hz和60 Hz之间)。

将2017年3月3日的测量数据对应频谱与包络线作对比，其结果如图8所示。

图8 一组频谱的包络报警

从8图中可看出，25 Hz处的幅值远大于包络线对应值，发生报警。这与测量当日，由于发动机的叶片磨损导致了振动烈度显著增加的情况一致。

按照式(3)和式(5)求取梯度序列和边缘复杂度，得CI=0.158 8，梯度的绝对值序列最大值为0.613 6，给出梯度阈值参考值ε=0.613 6×0.158 8=0.097 4。由双向边缘搜索算法得到上升边缘和下降边缘的序列位置，去除其中的重合边缘，其结果如图9所示(对图进行了适当放大)。

图9 梯度序列的边缘选取结果

图9中标出了梯度序列进行边缘检测的阈值上下限，依据双向边缘搜索算法得到上升边缘和下降边缘集合，按照坐标序列位置(此处不是频率)得到边缘位置，U={6，17，27，37}，D={16，25，35，45}。再根据边界处理规则，得到完整边缘序列，其中邻域大小δ=7(参考曲线平滑参数j)，则U={1，17，27，37，400}，D={16，25，35，45}。将两个序列合并，去除较窄的低能量区，得到最终的分割点序列Se={1，16，27，35，45，400}。分割点序列对应的均化频谱点如图10所示。

图10 频段划分结果

在图10中，划分出四个明显谱峰和一个较宽的低能量区，可见基于边缘检测的频谱自动划分方法能够有效分割高能量和低能量区。分割点序列对应的频率点为{0 Hz，40 Hz，67.5 Hz，87.5 Hz，112.5 Hz，1 000 Hz}，以此划分出五个频段，统计各段能量累积值，将中间各分割点对应幅值计入左侧频段，结果见表1。