马长刚,李 青,陈 明,孙大林
(1.空军勤务学院航空弹药系,江苏 徐州 221000;2.解放军94676 部队,上海 202150)
通过某装备检测参数可以判断出该装备的质量状态,但测试参数种类众多,数量庞大,其中包含了众多表征该装备质量状态的冗余测试参数。在进行其质量状态评估、故障预测、检测周期优化等内容时,无论是准确度还是计算的效率,均受到了数量庞大的测试参数群的影响。完全可从大量参数群中筛选出具有代表性的,同时与该装备质量状态相关性较强的一部分测试参数,去除大部分冗余检测参数,作为评估装备质量状态、预测装备故障等的状态特征量主成分,方便后续基于状态检测参数的问题研究。
假设某装备在k 个测试时刻其质量状态测试的测试参数进行收集,共有m 种状态测试参数,用Zt表示t 时刻检测参数值向量,zit表示在时刻t 质量状态测试参数的测试值,则Zt=(z1t,z2t,…,zmk)是一个k 维的随机列向量,而且假设每种质量状态检测参数所有的历史测试数据个数为n。由于该装备检测参数数据量比较大,如果考虑到当前t 时刻之前该装备质量状态检测所有历史检测时刻的检测参数,庞大的数据量会给计算带来很大麻烦,而且由于叠加效应,有可能会导致计算结果偏离实际值的情况。所以选取该装备质量状态检测特征量参数的关键是确定延迟期p,解决到底将历史检测时刻往后推到哪个时刻比较合理这一问题。利用时序分析法中多维自回归(AR)[1]确定该装备状态检测信息在时序上的延迟期p,具体步骤如下。
要利用自回归模型确定延迟期p,首先需要保证不同测试参数项在同一时间序列中数据值的均值为零,对均值不为零的原始质量状态检测数据做如下处理:
一阶差分之后,利用自协方差式(3)判断是否满足平稳序列要求:
按照一阶跟二阶的差分方法继续对非平稳序列进行差分处理,直至其满足平稳化条件。平稳化处理的目的就是消除质量状态检测参数在时间序列上的规律性和线性特点[3]。
由式(1)变换计算后得到该装备质量状态检测参数的平稳时间序列,建立此时k 维质量状态参数变量的AR 模型[4]。
根据以上定义,建立该装备质量状态检测参数的AR(p)模型:
式(6)可以表示为:
式(8)可以用下面的式子来表示:
要求解式(9)中的最小值,借助Q 关于Φ 的偏导来求解,令偏导为0,即:
即:
通过该装备质量状态检测参数的零均值处理、平稳化处理确定了时间序列延迟期p,也确定了其AR(p)模型,并对模型中的参数进行了估计,但确定的时间序列延迟期p 是否满足要求,是否精确,估计值是否合理,需要对其进行检验。在这里,根据估计值,取显著性水平α 进行检验。
可通过分析残差序列来检验模型的合适与否,采用Ljung 和Box 改进并提出的LB 统计量检验法进行检验:
假定原假设H0:模型延迟期p 符合要求;
H1:模型延迟期p 不符合要求;
主成分分析法(PCA)利用降维的思想,把多指标转化为少数几个综合指标,即主成分,其中每个主成分都能够反映原始变量的大部分信息,且所含信息互不重复[6]。PCA 方法中仅仅考虑到了状态监测信息间的互相关关系,而动态主成分DPCA 不仅考虑到了状态检测信息间的互相关而且考虑到同一状态检测信息不同时序时刻间的自相关[7]。
在求得模型阶次p 的基础上,考虑到该装备状态检测信息间既互相关又在不同时序上自相关,采用动态主成分分析法提取状态检测信息的主成分。步骤如下:
1)计算协方差矩阵
通过时序模型的计算检验,确定模型的阶次为p,状态检测信息间的协方差矩阵由互相关协方差矩阵C(0),以及自相关协方差矩阵C(1),C(2),…,C(p)组成,此时的协方差矩阵为一个由(p+1)×(p+1)个m×m 矩阵块组成的协方差矩阵,即C。每一矩阵块表示为:
i,j=1,2,…,p+1,Zit表示t 时刻状态检测信息i的检测值,表示其检测值的平均值。
2)确定相关矩阵
3)计算相关矩阵的特征值
4)提取主成分
主成分的提取,通过主成分的综合贡献率进行确定,根据步骤3)中的特征值的大小排序,前k 个主成分的综合贡献率[8]为:
根据综合贡献率要求就可确定出k 个主成分。
本章前面部分已经详细地介绍了该装备质量状态特征量选取的方法,方法合适与否,都需要根据真实的数据进行验证。从一线保障部队收集到了部分该装备寿命周期内完整的测试数据。用本文提出来的将时间序列和主成分分析法结合起来的动态主成分分析法来提取这几枚该装备质量状态特征量。选择收集到的某型号该装备完整历史测试数据,进行特征量的选取,该装备总共8 次历史测试次数的原始测试数据如表1 所示,表中C1~C90分别对应该装备各项参数。
表1 测试参数原始数据
根据表1 中该装备原始的8 次历史状态检测数据,下面利用前面一节中提到的该装备质量状态特征量参数选取方法提取该装备的状态特征量。具体实现步骤如下。
表1 中收集了该装备从出厂后到最近一次测试的数据,数据链比较完整,首先可以选择其中的一部分参数,从数据层面来了解一下随着检测次数的增加每一项测试参数数据值的变化趋势,在这里以表1 中的部分测试参数在8 次测试中数据值的变化情况为例,画出其变化曲线图,如图1 所示:
图1 原始参数数据变化趋势
图1 中选择了90 多项测试参数中的3 项,以它们为例,分析其原始数据的变化趋势,从图中可以看出,每项参数总共8 次测试的均值都不为零,数据波动起伏比较大,并没有严格围绕0 值线规律地上下波动,因此,数据值不满足模型求解的零均值要求,需要对其进行零均值处理。
收集到的该装备的测试参数项目总共90 项,如果在后面故障预测以及质量状态评估的过程中,考虑所有测试参数,将这么多的数据代入进行计算的话,计算起来比较困难,而且计算的准确度也受到很大的影响,前面已经分析过了,这90 项测试参数中,有很大一部分参数之间是相互关联的,它们的变化具有同步性。由图1 中每一项参数变化情况的分析中,其变化不具备零均值和平稳化的要求。根据特征量选取模型计算的要求,首先,需要对原始数据零均值、平稳化地处理。
1)数据零均值处理
2)数据平稳化处理
在满足零均值的情况下,根据平稳化处理的方法,结合式(2),对1)中零均值后的数据进行平稳化处理,并利用式(3)对其平稳化效果进行检验,最终可以得到满足后面质量状态特征量选取模型的零均值、平稳化后的数据。经过编写MATLAB 程序进行计算,平稳化处理后的数据满足平稳化要求。图1中3 项参数在8 次测试中的参数值经平稳后处理后的变化趋势图如图2 所示。
图2 平稳化后参数数据变化趋势
根据建立的该装备质量状态检测参数AR(p)的模型,要确定AR 模型的阶数p,需要求出式(7)、式(10)中的所有测试参数在P 个时刻内完整的测试数据矩阵X 和当前时刻质量状态参数测试数据矩阵Z。算法实现步骤如下:
Step1:首先假定延迟时序P=1;
Step2:根据p 的值,确定出矩阵X、Z;
Step6:如果通过检验则算法结束,认为p 值符合要求,若检验不通过,则令p=p+1,则跳转至Step2进入循环,直至p 的值通过检验,算法结束。
根据上面的算法步骤,编写MATLAB 代码进行计算,经计算,当P=2 时,通过检验,因此,最终确定p 的值为2。
确定出了时间序列模型的延迟时序P,并检验了P 值的合理性,运用主分量分析的方法提取出该装备的质量状态特征量,提取步骤及计算过程如下。
首先,根据前面所确定的延迟时序P 为2,所以在考虑同一参数不同检测时刻测试值间相互关系时考虑到往前推移2 个时刻,则本例中,根据式(14)计算协方差矩阵,这里计算出的协方差矩阵是由9 个90×90 维的矩阵组成的高维矩阵,即:
C11,…,C33分别为上面提到的9 个90×90 维的协方差矩阵,然后分别求出协方差矩阵所对应的相关矩阵R,并求解出相关矩阵R 的特征值和特征向量。由于协方差阵和相关阵维数较高,具体矩阵在这里不作展示。本例中,最终求解出的大部分特征值极小,这里仅列出部分特征值,如表2 所示。
表2 部分特征值
对所有特征值进行排序,并选取前100 项作出特征值散点图进行分析特征值分布规律,如3图所示。
图3 特征值分布规律
从图3 中可以看出来,特征值从第7 项开始已经趋近于零,直观判断,可以将前6 项特征值所对应的质量状态测试参数作为我们要选取的状态特征量。但根据前面建模过程中的方法,应从其特征值综合贡献率出发来计算,前7 项特征值综合贡献率的和为96.12%,满足对该装备质量状态特征量选取综合贡献率大于95%的要求。因此,最终选取出来的该装备质量状态特征量为6 项,分别是特征值排序中前6 项所对应的参数项。这6 项参数分别对应表1 中的C76、C77、C78、C79、C80、C81。
针对某装备测试参数众多,冗余参数较多,在利用参数进行该装备故障预测或者状态评估时计算很不方便,预测和评估的准确度受该装备状态参数维度高的影响较大。因此,借助时间序列分析确定该装备质量状态不同时序条件下的自相关、互相关关系,又利用动态主分量分析法从众多测试参数中选取出决定该装备质量状态的部分主分量,选取效果非常明显。