《武汉工程大学学报》 2014年02期
73-78
出版日期：2014-02-28
ISSN:1674-2869
CN:42-1779/TQ
马特拉算法在遥测数据短期预测中的应用
0引言遥测数据直接反映了卫星在轨运行的状态.卫星的遥测数据具有非平稳变化的特点，而具有非平稳性的遥测数据的各阶常用统计量（如均值、自相关函数等）经常随时间的变化而改变［1］，这给遥测数据的预测带来了很大的困难.传统的预测方法，如自回归模型（AutoRegressive model，简称AR）、自动回归移动平均模型（AutoRegressive and Moving Average model，简称ARMA）等常用预测模型，适合处理平稳的数据［2］，对于非平稳的数据不能满足预测需要.笔者通过对卫星历史遥测数据分析研究建立适用于非平稳时间序列的预测模型，对遥测参数的未来趋势进行短期预测，通过预测遥测数据序列潜在的故障趋势，为指导故障发生前的正确决策提供重要支持［3］，降低卫星潜在的风险，提高卫星运行可靠性. 笔者以卫星运行期间太阳翼的输出功率数据为研究对象.太阳翼的输出功率序列具有随机性和一定的周期性，这是不同频率分量的叠加结果.遥测数据的不同频率分量可分两部分，即慢变和快变，其中慢变部分反映了序列主体，快变部分体现了序列的细节［4］.对遥测数据进行短期的预测要求能够捕捉到在时域上表现瞬时、随机的分量.由于遥测数据的趋势中心随时间变化，传统的时域分析方法不能准确分析瞬时、随机的分量变化规律，而小波分析可对遥测数据信号中的频率分量进行粗细分离.1小波变换相关工作采用多分辨率分析思想用小波变换对遥测数据序列进行分解，得到反映序列主体变化的低频分量和代表序列细节的高频部分［5］.遥测数据原始序列根据所选择的小波基函数和适宜的分解尺度被分解为频率不同的分量，各分量在长度上与原序列保持一致［6］.低频分量代表原始序列中基本不变的主体部分［7］，即太阳翼基本输出功率；高频分量表示原序列中的瞬时、随机的分量，即太阳翼随机输出功率.1.1小波基选择在工程应用中小波函数除了要满足容许条件和正则性条件，还要满足以下3个条件：良好的紧支撑性；Ψ(t)具有消失矩；满足正交性.常用的小波函数，如Morlet、Mexican小波函数不存在尺度函数，即不满足正交性；Haar小波在时域上不是连续的，不适合做小波基.笔者选择dbN（Daubechies小波）小波函数做小波基，dbN小波是具有高阶消失矩的紧支集正交小波函数，阶数N的具体取值通过小波变换对实际序列分解的结果判定.图1是某同步卫星太阳翼输出功率序列采用dbN小波在N取1、2、3时进行分解后，近似部分与原序列的比较结果.N为3时，较好的体现了序列的变化趋势，且趋势中的突出点被保留下来；N为1、2时，序列的平滑度不利于分量预测模型的建立.所以，选择db3小波作为小波基对原序列进行分解.图1dbN小波分解结果Fig.1Decomposition result based on dbN第2期任国恒，等：马特拉算法在遥测数据短期预测中的应用武汉工程大学学报第36卷1.2分解尺度研究遥测数据的变化是缓变和快变结合，与之对应的是长度周期嵌套.这是时间序列的能量集中分布在一些频率带上的结果［8］.因此，将不同的频率分量分开，使其变化规律更加直观.对遥测序列进行分解结果是将原序列分解为几个细节部分和一个近似部分.分解尺度偏大，序列的采样密度会变稀，序列的主体曲线会越来越平滑，导致获得的近似序列失真；分解尺度若偏小，序列的主体变化趋势又不太明显，不易观察到各分量的变化规律［9］.图2是经db3小波分解后的近似部分aN在不同尺度下的分解结果.分析图2发现分解尺度图2不同分解尺度下近似部分比较Fig.2Comparation of aN with different scale为2时，近似部分的曲线a2已经足够光滑，同时保持了原曲线的形状；而a3、a4随着分解次数的增加，采样点减少，所得曲线过分平滑，序列的变化趋势已失真，因此本文选择的分解尺度为2.2基于小波变换的预测模型构建2.1周期自回归模型若有一时间序列X，表达式为Xt=a0t+a1tXt-1+…+aptXt-p+εt(1)满足以下条件：（1）εt为独立序列，且其期望值Eεt=0，方差Eε2t=σ2t；（2）对任意i=0,1,…p，有ait=ait+T，σ2t=σ2it+T，t=0,±1,…，其中T为一正数，则上述模型为周期自回归模型(Periodic Autoregression Model，简称：PAR模型)，T为PAR模型的周期长度，t为PAR模型的相位.此时遥测数据的预测模型为XKT+T=a0T+a1TXKT+T-1+…+apTXKT+T-p+εKT+T(2)设遥测数据每分钟样本为X1,X2,…,Xn，未来时刻h的值Xt(h)是Xt+h在时刻t的条件期望，并以此作为它的预测值，记为X^n(1)，则有X^n(1)=E(Xn+1|Xn-1,…,X0)=a01+a11Xn+…+ap1Xn-p-1(3)第k步预测值X^n(k)为X^n(k)=E(Xn+k|Xn,Xn-1,…,X0)=a01+a1kX^n+…+apKX^n-p+1(4)对于遥测数据分解序列建立一小时预测模型，选取周期的长度为60，即p=T=60.负荷序列的PAR模型如式(5)：X^n(k)=akT,k+akT+1,kX^t(k-1)+…+akT+60,kX^(k-60)(5)其中akT+i,k=ai,k(i,k=1,2,…,60),X^t(k-j)=Xt-k+i-1(k-j<0).遥测数据短期预测模型共有未知参数(p+1)T个，记为a(i)=(a0,i,a1,i,…,a60,i)(6)式(5)中的未知参数可由式(6)确定.已知ai拟合式(5)的残差平方和定义为St=\[Xn+k-a0,kXt(k-1)-…-a60,kXt(k-60)\]2(7)使得St达到最小值的ai为ai的最小二乘估计，其满足条件Staj,i=0 ，即Staj,i=aj,i∑ktk=0XKT+1-a0,i-a1,iXt(k-1)-…-a24,kXt(k-24)=0j=0,1,2,…,p(8)2.2遥测数据短期预测模型利用马特拉（Mallat）小波变换算法和db3小波基对太阳翼输出功率数据序列进行2尺度的分解，图3是对4小时内的输出功率负荷的原始序列和进行分解后的小波分量结果.图3太阳翼输出功率数据分解结果Fig.3Decomposition result of solar panel output power从图3可以看出对数据进行2尺度分解后，得到的低频数据a2具有较强的周期性，这是因为低频数据反映原序列的主体信息，变化规律较强；高频数据d1、d2体现原序列的细节，是待分解序列中变化较快的部分.由于分量d1、d2和a2的曲线变化特征不同，应根据不同分量的特点构建预测模型.笔者对各个子序列作如下区别对待：(1)主体信息a2周期性显著，变化相对平缓，对分解得到的15个点进行周期自回归预测；(2)高频分量d1、d2随机性较强，为提高预测的实时性，高频部分采用二次指数平滑法进行预测.指数平滑法由R.G.Brown提出，该算法的优点是预测模型构建过程中只需少量的历史数据，计算量小，便于实时预测.模型建立过程中是利用算法对原始时间序列的不规则性进行平滑，获得原序列的变化规律和趋势，对未来某时刻的数据进行预测，更多的考虑数据的更新.指数平滑法计算公式如式(9)：Yt=yt=axt-1+(1-a)yt-1(9)式(9)中:xt-1、yt-1分别是t-1时刻的实际值和预测值；yt为t时刻的平滑值；Yt为t时刻的预测值；a为平滑系数，取值范围为0