李志强 张香燕 田华东

(北京空间飞行器总体设计部,北京 100094)

卫星遥测数据包含了星上各部件在轨运行的丰富信息,地面通常根据下传的遥测数据对其运行情况进行监视和评估。研究分析在轨卫星的遥测数据变化规律,在此基础上对遥测数据进行长期和短期预测,不仅有助于发现各部件和性能的老化衰退程度,还能根据预测数据进行提前预警和故障预示,以便提前发现非预期的异常变化并及时有效地处置,对降低卫星在轨运行风险、确保长期安全可靠运行具有重要意义[1-2]。

由于卫星面临复杂的空间环境,影响遥测参数变化的空间因素较多,如辐射粒子、中性大气、地磁等空间要素,还有运行轨道、姿态等影响,这些不确定或者周期性因素与卫星各部件相互作用产生效应,造成遥测数据呈现不同规律的波动,给遥测数据预测带来了困难。遥测数据的变化规律是时间序列的各种统计特征的叠加,既有趋势性,又有季节性,还有周期性的和不规则的变动,这些都反映了遥测数据变动规律的复杂性。研究时间序列数据的统计规律,并假定这种规律在未来时刻仍将发挥作用,就可以建立合适的时间序列数学模型来描述这种规律,并利用它作出合理的预测。

传统上,对遥测数据进行预测一般采用几种常见的预测方法,如曲线拟合、返向传播(BP)神经网络、小波分析、灰色系统预测等。曲线拟合方法一般适用于单一或者固定变化的简单预测,且在长期预测时容易偏离原有的趋势而发散[3];BP神经网络预测依赖于大量样本训练和算法的有效性,且一般仅用于短期预测,对长期预测误差较大[4];小波分析方法在进行剔除噪声和提取周期处理时,容易遗漏数据突变和引起信息失真[5];灰色系统预测一般仅对单调变化具有较好的预测精度,对叠加了周期性和随机性的数据预测精度不够,且预测有效期较短[6-7]。HP(Hodrick-Prescott)滤波是一种常用的时间序列分析工具[8-9]。基于时间序列分析原理,采用HP滤波将数据分解成相应时间序列的趋势成分、波动成分,可弱化序列的趋势性、波动性等因素之间的相互影响,用来分析不同序列的规律,同时保留原序列的完整数据特性,为进一步的数据预测提供了基础。HP 滤波已在工业能源分析、经济预测等领域应用,将HP滤波应用在卫星遥测数据分析中尚未看到相关报道。

本文将HP 滤波引入到遥测数据分析与预测中,提出应用HP滤波进行时间序列分解的遥测参数预测方法,即应用HP滤波对数据序列进行趋势成分和波动成分的分解,分开建模、组合预测,最终通过叠加得到参数的预测值,具有较高的预测精度,可为遥测参数的分析和预测、故障诊断和预警提供一种新的技术方法。

1 预测方法

卫星遥测数据的变动规律一般既包含有趋势性的长期规律,又包含有季节性、周期性和不规则变动的短期波动规律。基于此,本文提出的预测方法是:首先,应用HP滤波将卫星遥测数据Y 分解成趋势成分数据G 和波动成分数据C;其次,根据趋势成分和波动成分的不同特点,对两者分别建模进行预测;最后,将2个预测结果进行叠加,得到原始遥测数据的组合预测结果。预测方法过程见图1。对遥测数据应用HP滤波,并没有消除原始遥测数据本身的趋势性、周期性和季节性,这些性质被保留在分解出的趋势成分G 与波动成分C 中。观察两类数据成分的特征,发现趋势成分G 带有显著的线性趋势性,因此考虑建立多元自回归模型进行预测;对波动成分C进行检验,发现C 明显带有周期性和季节性,考虑建立季节型单整自回归移动平均(ARIMA)模型进行预测。

图1 预测方法过程Fig.1 Process of prediction method

1.1 方法步骤

预测方法的具体步骤如下。

(1)根据HP 滤波,将原始遥测数据Y＝{y1,y2,…,yn}进行分解,得到带长期趋势性的趋势成分G＝{g1,g2,…,gn}与带短期波动性的波动成分C＝{c1,c2,…,cn}[8-9]。它们与原始遥测数据的关系为

式中:t＝1,2,…,n,其中,n 为采样的数据个数。

HP滤波目的在于将数据的趋势成分从不平滑的原始遥测数据中分解出来,分离过程必须满足损失函数最小原则,即

式中:平滑参数λ＝σ12/σ22,σ12和σ22分别为趋势成分G 和波动成分C 的方差。

(2)将得到的趋势成分G 视作t的函数,G 具有明显的线性趋势,定义各个数据点到自回归曲线的误差项为ε1,ε2,…,εn。根据min{ε12＋ε22＋…＋εn2}确定拟合多元自回归多项式。

式中:G(t－i)表示滞后项的函数,即gt同时与gt－i的线性组合有关。

(3)根据拟合多项式,对趋势成分G 进行预测,记预测值为。

(4)检验波动成分C 的平稳性:若C 平稳,则进行相关性分析;否则,对C 进行差分运算直至平稳,并记为差分数据{wC}。

(5)对差分数据{wC}进行相关性分析,检验数据的季节性,确定模型的阶数,建立季节型ARIMA模型。ARIMA 模型是一种精度较高的时间序列预测方法,其应用于遥测数据yt的表达式为[10]

式中:εt为白噪声;L 为滞后算子;通过差分得到的平稳时间序列Δdyt＝(1－L)dyt;Φ(L)＝1－φ1Lφ2L2－…－φpLp和Θ(L)＝1＋θ1L＋θ2L2＋…＋θqLq分别为自回归算子和移动平均算子多项式,其中,φ 为自回归系数,θ 为移动平均系数,p 和q 分别为非季节自回归和移动平均算子的最大滞后阶数。

针对遥测数据的季节性变动规律,考虑应用季节型ARIMA 模型[11-12],即

式中:φp(L)为非季节自回归算子多项式;ΦP(Ls)为季节自回归算子多项式;θq(L)为非季节移动平均算子多项式;ΘQ(Ls)为季节移动平均算子多项式;D 为季节差分次数;P,Q 分别为季节自回归和移动平均算子的最大滞后阶数;s 为周期。

式(5)即季节型ARIMA(p,d,q)(P,D,Q)s模型,依据这一模型可得到波动成分C 的预测公式为

(6)检验残差数据{εt}是否为白噪声,通过对残差进行单位根扩展的迪克富勒(ADF)检验,当相应的t统计量小于一定置信水平的临界值时,表明残差序列不存在单位根,序列是平稳的,即为白噪声过程;判断季节型ARIMA 模型的合理性。

(7)根据季节型ARIMA 模型对波动成分C 进行预测,预测值记为。

(8)根据HP滤波原理,得到原始遥测数据的组合预测结果为

(9)对模型预测结果进行误差分析,检验预测数据的相对误差为

1.2 关键问题

本文预测方法需要解决的关键问题如下。

(1)HP滤波平滑参数λ 取值。λ 用来调节趋势成分对原始遥测数据的跟踪程度与趋势成分的光滑程度。根据不同周期的数据和计算需求,λ 取值不同,随着λ 值的增加,估计的趋势更加平滑;当λ 无穷大时,估计的趋势将接近线性函数。一般根据经验和试算结果,综合考虑趋势的平滑程度并有利于波动成分ARIM A 建模来确定λ。

(2)季节型ARIMA 模型的模型识别及确定。通过分析波动成分的自相关函数和偏自相关函数的截尾性和拖尾性来识别模型的阶数,并检验模型参数估计后的残差,通过比较不同模型的调整拟合优度R2、赤池信息准则(AIC)值和施瓦茨准则(SC)值,筛选确定最优预测模型。其中:调整的拟合优度也称为修正样本决定系数,它表示总离差平方和中由回归方程可以解释的部分所占的比例,这一比例越大,表明回归方程可以解释的部分越多,模型越精确,回归的效果越显著;它一般介于0～1,越接近1说明回归拟合效果越好。AIC 和SC 均是基于对数似然函数的统计量,它们的值越小越好,越小意味着模型的简洁性和精确性都较好。选择调整拟合优度R2较大、同时其AIC 值和SC 值较小的模型,即确定了最优预测模型。

本文预测方法根据遥测数据的趋势成分G 和波动成分C 的不同特点分别建立模型,既有效地拟合了原数据的趋势性,也降低了季节性对趋势性的影响,从而建立了较为精确的组合预测模型。

2 预测方法实现及实例验证

依据本文方法,设计遥测数据预测和预警系统方案,并进行实例验证。

2.1 预测和预警系统方案

系统方案基于预测方法,采用分布式模块化设计,由数据接收分发子系统、数据预测子系统、预警管理子系统、任务配置和中心控制子系统、故障诊断子系统及数据库管理子系统组成,如图2所示。用户终端可接收任务配置和中心控制子系统发送的预测结果信息,用于在轨管理日常监视和预警。

图2 系统组成Fig.2 System composition

数据预测子系统中是系统方案中的核心子系统,构成了预测方法的实现过程,主要包括时间序列分析模块、预测算法配置模块、实时预测模块和历史数据诊断模块。其中:时间序列分析模块配置HP滤波方法,对遥测数据进行趋势成分和波动成分的分解;预测算法模块接收来自时间序列分析模块的分解数据,配置多项式自回归算法模型和ARIMA模型,分别进行建模预测,并叠加二者预测值,最终得到组合预测结果;实时在轨遥测数据和历史数据分别调用时间序列分析模块和预测算法模块,分别进行实时预测和历史数据诊断,输出预测和诊断信息。各模块的结果发送给任务配置和中心控制子系统。

预警管理子系统获取预测数据和提前预警阈值配置文件,并进行预测数据和预警阈值的比对,当预测结果超出预警阈值时,输出相应的预警信息,并将信息发送给任务配置和中心控制子系统。故障诊断管理子系统获取预测数据和历史数据,配置诊断知识,并完成预测数据和历史数据的比对,当预测结果与历史数据规律不一致时,输出相应的诊断信息,并根据诊断提示故障信息,同时将信息发送给任务配置和中心控制子系统。

任务配置和中心控制子系统完成整个系统的控制和运行管理,主要完成预测参数的选择、预测任务的配置、平滑参数设定、ARIMA 模型识别、输出结果管理等,尤其要解决预测方法实现中的两个难点:一是时间序列分析模块中的HP 滤波参数λ 的选择,需要根据遥测数据分解结果和数据特点,反复迭代、综合考虑来确定;二是ARIMA 模型的识别,在预测算法模块计算数据自相关函数和偏自相关函数,通过分析其截尾性和拖尾性识别模型的阶数,完成模型的确定和筛选。

2.2 实例验证

为验证预测方法的有效性,本文选取某卫星行波管阳压(阳极电压)2012年1月－2021年4月的在轨遥测数据用来验证。其中:2012年1月－2020年10月的数据作为模型识别和建模的参考数据;再用此模型预测2020年11月－2021年4月的数据,并与同期实际数据比较并计算预测误差。本文利用Matlab和Eviews软件对数据进行建模分析。

2012－2021年阳压随时间变化趋势如图3 所示。阳压随时间一直有缓慢上升趋势,且每年阳压有季节性变动。这是由于随着行波管工作时间增加,其内部电极有老化衰退现象,表现为阳压有逐渐缓慢上升的趋势;同时又受到运行轨道和光照影响,阳压又与行波管的温度环境有关,每年表现出一定的波峰波动性,即季节性变动。针对此遥测数据,首先采用HP滤波分解出趋势成分和波动成分。滤波处理的关键是参数λ 的选择,经过反复试算拟合,综合考虑趋势的光滑程度并有利于波动成分建模,最终确定λ 取值为14 400。

原始遥测数据经HP 滤波分解后如图3所示,趋势成分G 曲线代表长期趋势,且增长率在寿命初期较为明显,中后期较为平滑。波动成分C 曲线代表包含季节性、周期性的波动成分,每年波动情况不同。

图3 遥测数据HP滤波分解结果Fig.3 Telemetry data separation results using HP filter

2.2.1 趋势成分预测

趋势成分G 近似一条平滑的曲线,记t(t＝1,2,3,…)为累积时间,依据式(3),将G 对t进行回归拟合,根据多项式曲线拟合的结果调整次数和滞后项。选择2012年1月－2020年10月数据拟合模型,用来预测2020年11月起未来6个月数据。依据点到拟合曲线的误差平方和最小原则,得到模型估计系数计入表1,即建立了趋势成分的预测模型。由表1可知:多项式自回归曲线对G 的拟合效果非常好,R2和调整后的R2达到了1.000;AIC值和SC值分别为－9.714和－9.561,模型的精确性较好。各参数均通过了显著性检验。可见,趋势gt不仅与时序t有关,还与2阶滞后项有关。

表1 趋势模型估计结果Table 1 Estimated results by trend model

根据模型估计结果和式(3),对2020年11月－2021年4月的趋势成分G 进行预测,得到预测曲线如图4所示。从图4中可见:置信区间很窄,采用多项式自回归曲线对趋势成分进行预测时,对整体平滑走势预测效果较好。

图4 趋势成分预测值及置信区间Fig.4 Prediction values and confidence interval of trend term

2.2.2 波动成分预测

通过建立季节性ARIMA 模型,对波动成分C进行分析与预测。经分析,HP 滤波后的波动成分C 自相关曲线没有呈指数快速衰减,表明数据的趋势性没有完全消除。为进一步消除趋势性,对C 进行1阶差分,差分数据是平稳的,故d＝1。分析1阶差分数据的自相关函数和偏自相关函数,发现C 也具有季节变化规律,周期为12;为了消除季节性,进行滞后12期的季节差分处理,故D＝1,s＝12。分析季节差分数据的偏自相关函数和自相关函数[10],可知p＝3,q＝3或1,P＝1,Q＝1,考虑建立模型ARIMA(3,1,1)(1,1,1)12或ARIMA(3,1,3)(1,1,1)12。在确定模型阶数后,要筛选出最优模型。一般,调整的R2越大、AIC值和SC值越小,相应的模型越好。调整的R2较大的模型是ARIMA(3,1,1)(1,1,1)12,同时其AIC值和SC值也相对较小,因此较合适的模型是ARIMA(3,1,1)(1,1,1)12。模型估计结果见表2,其R2和调整后的R2分别为0.595和0.560,AIC值和SC值分别为2.960和3.171,相对其他模型较好,对波动能进行合适的刻画。

表2 波动模型估计结果Table 2 Estimated results by fluctuation model

在模型估计之后,需要对模型的残差是否为白噪声进行检验,通过白噪声检验才能进行预测。对残差进行单位根ADF 检验后,可确定残差为白噪声。因此,ARIMA(3,1,1)(1,1,1)12模型为理想预测模型,可以用于预测。依据这一模型和式(6),对波动成分2020年11月－2021年4月的数据进行预测,得到预测曲线见图5。

图5 波动成分预测值及置信区间Fig.5 Prediction values and confidence interval of fluctuation term

图5显示模型可以准确地预测随时间变化的波动性,并随着预测期的延长,预测置信区间也变大,表明越往后模型的预测精度越差。这是由于动态模型中的AR 项和MA 项依赖于新的信息,而新信息的不确定性逐渐积累,从而导致新的预测精度下降。

2.2.3 组合预测及验证

根据HP滤波原理yt＝gt＋ct,利用前述模型进行组合预测,得到2020年11月－2021年4月的阳压预测结果,并与实际值比较,进行可靠性验证。灰色理论预测是一种应用广泛的数据预测方法,一般可应用于对趋势性数据的预测。为验证本文预测方法的有效性,利用灰色模型GM(1,1)对数据进行同期预测,结果与本文方法对比,如表3 所示。从表3中看出:利用组合模型进行预测,得到的预测结果非常精确,预测值与实际值差距非常小,绝对误差在－1.6～－0.1,相对误差都在0.04%之内,最小相对误差仅为－0.003%,如图6所示,在半年时间的预测区间内,达到了很高的精度。作为对比,灰色理论预测值绝对误差明显偏大。同时,本文预测方法获得的预测值也准确模拟跟踪了实际曲线的波动,跟踪实时性较强,这主要是由于ARIMA 模型对波动细节和随机性的刻画程度高。相比较而言,灰色模型GM(1,1)模型的预测值仅仅反映了趋势性增长,而对波动性预测不足,导致预测结果不能紧密跟踪实际值变化,其预测绝对误差也远大于组合预测误差。因此,本文预测方法的组合模型识别及验证结果理想,用近10年的数据预测半年的数据,计算时间约5 min,在行波管状态不发生突变的情况下,可利用此模型对行波管在轨阳压数据进行较为快速准确的预测。

图6 组合模型预测、灰色模型预测及相对误差Fig.6 Combined model prediction,GM mode prediction and relative errors

表3 2020年11月-2021年4月阳压预测结果Table 3 Prediction results of anode voltage from November 2020 to April 2021

事实上,阳压遥测数据是由遥测原码通过数据处理方法计算得到的,其原码变化的最小单位是一个分层值,而一个分层值对应计算得到的阳压数值约为2.00 V;因此该组合模型的预测精度(约1.65 V)已经达到预测一个分层值的微小变化,甚至先于遥测值更本质地反映其变化趋势。由于阳压受空间温度环境、行波管自身性能衰退等诸多外界和自身因素的影响,趋势成分和波动成分相互叠加,使得基于物理机制的预测得到的预测精度往往不佳。而本文预测方法的验证结果表明:该方法不需要考虑多种复杂的物理和外界因素,仅仅依据数据本身的变化规律,在较长时间的跨度内能够实现优于单个分层值的高精度预测,可以有效提升对遥测参数的预测能力。

本文预测方法还可以应用于在轨故障诊断和寿命预测中。卫星数据长期表现包含有各种因素共同作用的影响,评估性能时分离出各种成分的数据是必要的,可以更精准地发现数据变化的规律。在实际卫星的在轨监视中,利用本文预测方法对卫星历史遥测数据进行回溯建模,可以发现参数理论预测值与实际值的差别,当实际值偏离预测值较大时,就可能诊断出相应设备的潜在故障状态;同时,针对老化部件和性能衰退部件,利用此方法可以对设备长期性能进行诊断和预警,设置动态预警监视范围,提前预警异常问题,进而实现对运行寿命的有效预测。

3 结束语

本文将HP滤波应用到卫星遥测数据领域,提出了一种基于HP滤波分解的遥测参数预测方法。将原始遥测数据应用HP滤波按照趋势成分、波动成分分解,弱化了数据的趋势性、波动性之间的相互影响,从而更精确地建立组合预测模型。通过在轨阳压数据实例分析,验证了该组合预测方法的有效性。该方法具有拟合精度高、预测精度高的特点,同时也能保证在长时间跨度内多步预测的精度。基于该方法提出的应用验证方案,可在卫星故障诊断预警、寿命预测等方面实际应用。同时,将HP滤波思想应用到遥测数据,引入时间序列分析的成熟工具和模型,可为卫星状态监视、健康评估与故障预警提供新的思路和技术方法,为在轨管理提供高效有力的分析手段,具有重要的工程应用价值。