短时间序列的复杂程度分析方法研究*

陈强强 成建波 张 刚 张 婕

（1.海军研究院 上海 200436）（2.海军航空大学 烟台 264000）

1 引言

随着航空工业的飞速发展，航空装备的地位作用及使用效率都得到了极大的提升。同时，随着航空科技的不断进步，航空装备的复杂性及综合性也越来越高，保证航空装备具有良好的工作状态也显得尤为重要，这对航空装备的可靠性提出了更高的要求。航空发动机健康状态监视技术是航空装备视情维修的前提及基础，同时也是提高航空发动机使用效率并保证航空发动机发挥其稳定效能的重要保障。

针对航空发动机气路参数进行准确有效的分析，是实现航空发动机健康状态监视的重要步骤及方法。文献[1]通过对涡桨发动机转子系统振动信号进行时域特征、能量特征分析，完成了涡桨发动机转子的故障诊断；文献[2]以民航发动机排气温度裕度作为研究参数，通过对排气温度裕度序列进行分析，从而判断水洗对拆发间隔的影响。文献[3]通过对飞参系统记录的发动机状态参数时间序列进行分析，得到反映发动机状态参数的特征信息。对发动机气路参数时间序列进行分析，能够为发动机状态监测提供决策依据[4]。

在对时间序列进行复杂程度分析方面，排列熵算法具有计算简单、效率高、抗噪声能力强、适合在线监测等优点，在时间序列分析、故障特征提取、故障诊断等方面具有重要应用[5]。文献[6]采用排列熵方法对转子振动信号进行检测并实现特征提取、故障诊断；文献[7]对航空器陀螺仪序列进行排列熵算法检测，提取其故障特征并实现故障诊断。文献[8]对液压泵参数信号进行排列熵分析，验证了其有效性。排列熵算法在机械系统时间序列分析中，取得了重要应用。但在解算过程中，受限于排列熵算法的特征，所采集的时间序列长度均设置在较长的合理范围内。然而航空发动机气路参数具有获取较为困难、数据量较少、样本长度有限等局限性，所采集到的数据往往具有序列长度短的特征[9]。为了弥补排列熵算法在短时间序列复杂程度计算方面的不足，本文通过引入插值思想，实现了对原始数据的延拓，并将原始序列重构至排列熵计算的合理范围。通过对排气温度裕度（Delta Exhaust Gas Temperature，DEGT）时间序列及随机序列进行实验分析，结果表明所提出的基于插值重构的排列熵算法可以有效地表达短时间序列的复杂程度，是一种有效的复杂性分析方法。

2 基于排列熵的复杂性分析

排列熵算法通过衡量时间序列内部排序过程，从而实现对时间序列复杂程度的衡量，排列熵算法的基本原理如下。

1）首先，针对长度为N的原始时间序列{x(i)，i=1，2，...，N} ，对x(i)进行相空间重构分析，得到矩阵如下[10]：

式中，m为嵌入维数，τ为延迟时间，k=N-(m-1)τ为重构向量的个数。

2）将Xj中的数据元素按照升序的方式进行排列，得：

式中，j1，j2，…jm为各元素在排序之前位于相空间所在列的索引。

3）如果Xj中存在两元素相等，则按原始顺序进行排列。通过对相空间重构后的时间序列分析，对于任意一个Xj，均能得到相应的符号序列Sl={j1，j2，…，jm} ，其中，l=1，2，...，k，且k≤m!。可定义排列熵的计算公式为[11]

式中，Pj为符号序列的概率，且。

此时，Hp的取值为[0， 1]，可反映出时间序列的复杂程度。Hp越大，则序列复杂程度越高。

排列熵作为衡量时间序列内部复杂度的算法，计算过程简单方便，能够有效地分析非线性、非平稳信号，适用于发动机气路参数数据。在排列熵的计算过程中主要包括3个参数：时间序列长度N、嵌入维数m、延迟时间τ。排列熵算法的提出者Bandt建议选择嵌入维数m为3～7，一般随着时间序列长度增加进行相应调整[12]。因为如果嵌入维数m过小，相空间重构的向量中包含状态太少，算法失去意义和有效性；而嵌入维数m过大，则相空间重构步骤会均匀化时间序列，无法反映序列的微小变化。延迟时间τ对时间序列的计算影响较小，一般选择τ=1。文献[13]以随机信号为例对排列熵计算过程中时间序列长度的影响进行了分析，证明了时间序列长度在1024或2048左右较为合适，此时可对应的选择嵌入维数m=6，延迟时间τ=1。通过文献中对排列熵的分析可知，如果时间序列长度过短，则嵌入维数的取值应进行相应的减小，而过小的嵌入维数又会导致排列熵算法失去其有效性，从而影响复杂度的衡量[14]。因此，有必要对短时间序列进行进一步处理，使其适用于排列熵算法。

3 三次样条插值

为了对采集到的原始DEGT短序列进行处理，使其时间序列长度适用于排列熵的解算，借鉴经验模态分解（Empirical Mode Decomposition，EMD）算法中构造上下包络曲线的原理[15]，引入三次样条插值（Cubic Splines Interpolation）思想[16]对原始短时间序列进行重构。三次样条插值方法得到的曲线具有良好的平滑性，逼近效果好，能够很好地实现数据的插值。

在经验模态分解的过程中，三次样条插值的运用主要是为了求取原始时间序列的固有模态函数（Intrinsic Mode Function，IMF）。其中心思想为对原始时间序列x(i)进行分析，分别求取其极大值点和极小值点，然后对确定的极大值点和极小值点，采用三次样条曲线进行拟合，得到时间序列x(i)的上包络线U(i)及下包络线L(i)，并求取U(i)与L(i)的均值，得到中值包络线M(i)；以此类推，得到相应的固有模态函数，从而完成对整个时间序列的分解。

借鉴经验模态分解过程中对三次样条插值的应用，本文提出的短时间序列复杂程度分析方法，通过对原始短时间序列进行三次样条插值，从而得到符合排列熵解算过程所需要长度的新时间序列。

三次样条插值法是一种特殊多项式进行插值的方法，可以减小低阶多项式的插值误差[17]。线性数学分析中常用的线性插值（Linear）在估计两个数据点之间的数值时较为常用也较为简单，但线性插值在使用过程中，对数值的模拟较为粗糙，在曲线之间的连接处不够平滑。相比之下，三次样条插值使用一条平滑曲线对各数值点进行模拟，经过实验验证，其结果比线性插值更接近工程实际情况，所得结果也更具有优势。

对样条的定义是用来表述光滑曲线的一种方法，在数学模型中，为了实现光滑的闭合曲线，通常需要在样点上固定样条。非样点区域可以任意弯曲，这样连接样点的曲线被定义为样条曲线[18]。

三次样条插值的数学定义为

假设对y=f(x)在区间[a，b]上给定一组节点a=x0＜x1＜x2…＜xn=b，并且给出相应的函数值y0，y1，y2…yn，如果s(x)具有以下特征：

1）每个子区间[xi-1，xi]，其中i=1，2，3，…，n上s(x)是不高于三次的多项式；

2）其中s(x)，s′(x) ，s″(x)在区间[a，b]上连续，s(x)为三次样条函数；

3）s(xi)=yi，其中其中i=1，2，3，…，n，s(x)是s(x)=y的三次插值函数。

4 短时间序列插值重构

基于三次样条插值理论，针对长度为N的原始时间序列{x(i)，i=1，2，...，N} ，原始时间序列的数据点间隔可默认为1，即每隔1个横坐标数据点，存在与其对应的1个纵坐标数据点。构造短时间序列插值重构算法模型如下。

1）设原始时间序列{x(i)，i=1，2，...，N} ，原始时间序列可以看做x(i)与i之间的函数，每个x(i)对应的坐标位置分别为i=1，2，3，…，n，为了直观展示横坐标i与纵坐标x(i)之间的关系，也可表述为i=1，1+1，2*1+1，…，(n-1)*1+1；

2）参考小波分析（Wavelet Transform，WT）中的尺度因子变化方法，引入重构维数k，在保证每个纵坐标x(i)数值不变的前提下，对原始数据序列的横坐标位置进行重构，得到重构后的横坐标位置对应为i′=1，k+1，2k+1，…，(n-1)k+1，此时构建得到x(i)与i′之间的函数，即为重构时间序列；

3）以原始时间序列x(i)的值为纵坐标，以重构后的坐标位置i′为横坐标，补全后的时间序列间隔仍为1，进行三次样条插值，最终得到重构后的时间序列x′(i)长度为(n-1)k+1。

针对x(i) 与x′(i)进行分析，本文提出的基于三次样条插值方法的时间序列重构模型遵循原始时间序列结构，通过改变时间序列对应的坐标，对时间序列进行延拓，提取出原始时间序列的内部信息，并按照原有的时间间隔重新分布索引位置。该时间序列重构方法未改变原始数据结构，遵循了原始时间序列规律。

基于此，确定短时间序列的复杂程度分析方法如下：

1）针对原始短时间序列，确定其数据长度及需要重构的维数；

2）对原始序列进行重构，将其程度延拓至排列熵解算所需要的合理长度；

3）计算排列熵值，完成对重构后的时间序列复杂程度分析。

为了验证本文所提出的基于三次样条插值重构方法的合理性，利用正弦信号进行分析实验，在[0 ，6]区间内（约为[0 ，2π]区间，即一个周期，实验中2π的值不为整数，为了对实验结果进行直观展示，取整选择[0 ，6]区间），以步长x等于0.1，构建正弦信号y=sinx。此时y与x之间的映射为y=f(x)，即每一个以0.1为横坐标间隔的点，对应一个纵坐标的值。

以重构维数取4为例，原始时间序列长度为（61-1）*1+1=61，则重构后的时间序列长度为（61-1）*4+1=241。

针对原始时间序列y=f(x)，其横坐标x分别对应为0，0.1，0.2…，64，重构后的横坐标x′分别对应为0，0.4，0.8，…，241；保持纵坐标的取值不变，得到映射y′=g(x′)。

保证时间序列的间隔与原始时间序列间隔相等，即为0.1，对y′=g(x′)进行三次样条插值并补全时间序列。得到时间序列间隔为0.1的y′=g(x′)。

原始时间序列y=f(x)=sinx与重构后时间序列y′=g(x′)之间的对比如图1所示。

图1 原始时间序列与重构后时间序列

如图1所示，通过原始时间序列y=f(x)与重构后的时间序列y′=g(x′)之间的对比可知，本文所提出的基于三次样条插值的时间序列重构算法仅通过插值重构增加了时间序列长度，将原始时间序列的长度从61变为241，没有改变原始时间序列的表现形式，仍然保持原始时间序列的内部规律。

验证1：相关性验证。为了从数值形式上直观表示原始时间序列与重构后时间序列的关系，在[0 ， 24]区间内，以步长0.1对横坐标进行取值，构建正弦信号，并与重构后的时间序列进行对比。采用互相关系数ρy′y″表示重构后时间序列y′与在[0 ， 24]区间内的正弦信号y″之间的相关性；ρx3t的取值越接近于1，则表示两种信号之间的相关性越强[19]。

计算得到在[0 ，24]区间内的正弦信号y″与重构后时间序列y′之间的互相关系数为0.9956，表明两者呈现强烈的相关性，重构后的时间序列未改变原始时间序列的物理特性，仍然遵循了原始时间序列的内部规律。

验证2：误差验证。为了验证重构后时间序列y′与 [0 ，24]区间内的正弦信号y″之间的误差，采用均方根误差（Root Mean Squared Error，RMSE）为指标对重构时间序列的误差进行衡量[20]。RMSE的表达式为

根据式（6），计算得到正弦信号y″与重构后时间序列y′之间的RMSE为0.0680。表明三次样条插值具有很好的拟合效果。

5 实验分析

以国航某发动机的巡航数据DEGT观测值为实验样本，采样间隔约为15个飞行循环，可近似为等间隔采样。从该型发动机第1个飞行循环到第12501个飞行循环，共计101个数据采样点。原始DEGT数据序列如图2所示。

图2 原始DEGT时间序列

如图2所示，原始时间序列长度为101，当时间序列长度过短时，由于嵌入维数的选择问题，排列熵算法失去了其对时间序列复杂程度衡量的意义，无法准确衡量时间序列的复杂度。根据文献[13]中对排列熵算法的分析，选取数据长度约为1024较为理想，此时选择嵌入维数m=6，延迟时间τ=1。根据本文提出的短时间序列插值重构算法对其进行重构，确定重构维数k=10，则重构后的时间序列长度为(101-1)*10+1=1001。原始序列与重构后序列的对比如图3所示。

图3 原始序列与重构后序列

如图3所示，通过重构前后的时间序列对比可知，本文所提出的基于三次样条插值的时间序列重构算法仅通过插值重构增加了时间序列长度，没有改变原始时间序列的内部规律。对重构后的序列进行排列熵解算，得到其熵值为0.2964。

为验证本文所提出的短时间序列插值重构方法的合理性，在延迟时间τ=1时，验证在不同嵌入维数下，原始时间序列和重构后时间序列的排列熵变化情况。由于原始DEGT序列长度过短，因此在计算时嵌入维数的选择不宜过大，排列熵计算结果如表1所示。

表1 DEGT序列排列熵计算结果

根据排列熵的定义可知，排列熵的求解取值范围为[0， 1]，数值越小，表明时间序列复杂性越低。如表1所示，随着嵌入维数的变化，短时间序列的排列熵变化范围过小，无法准确反映排列熵结果随嵌入维数的变化过程，不能准确衡量时间序列的复杂程度。

为了进一步验证插值重构算法在排列熵计算过程中的有效性，选取长度为101的随机白噪声序列，并与重构后的时间序列进行对比验证，所得结果如表2所示。

表2 随机序列排列熵计算结果

如表2所示，直接对短时间序列进行排列熵求解，仍然无法准确反映排列熵结果随嵌入维数的变化过程，同时，对比表1及表2，也验证了在排列熵求解过程中，嵌入维数的取值不宜过小（如嵌入维数m=2时，均无法有效衡量时间序列复杂程度），否则相空间重构的向量中包含状态太少，算法失去意义和有效性。

6 结语

排列熵算法能够准确有效地衡量时间序列复杂程度，但在工程应用中，受限于数据样本获取困难的不足，时间序列呈现长度较短的特征。本文通过对排列熵原理进行分析，针对短时间序列无法准确使用排列熵衡量其时间序列复杂程度的问题，提出了采用基于三次样条插值的时间序列重构方法对原始数据进行延拓。以随机序列及采集到的航空发动机气路参数序列为实验对象，对插值重构方法进行了实验验证，实验结果表明：

1）基于三次样条插值的短时间序列重构模型，能够将短时间序列延拓至合理长度，为排列熵求解过程奠定基础，有效解决了短时间序列复杂程度无法准确衡量的问题；同时，时间序列长度的确定，也避免了排列熵求解中的参数选择问题；

2）采用插值重构算法，未改变原始时间序列结构，遵循了原始序列规律；

3）针对不同序列，均可采用插值重构算法进行延拓，在工程中具有一定的普适性。