CN106156499A - 一种基于Fourier级数模型的时序数据预测方法 - Google Patents

Abstract

本发明公开了一种基于Fourier级数模型的时序数据预测方法。使用本发明能够对在轨卫星的遥测数据进行有效预测，且预测精度较高。本发明首先对遥测采样数据进行预处理，剔除粗大误差、剔野，然后建立Fourier级数模型，并针对遥测采样数据特点，综合应用了快速傅里叶变换、矩阵的QR分解算法以及最小二乘法，确定Fourier级数模型的各个参数，最后利用Fourier级数模型对卫星遥测数据进行预测，获得较好的预测结果，与实际的遥测数据误差小，预测精度高。

Description

一种基于Fourier级数模型的时序数据预测方法

技术领域

本发明涉及卫星遥测数据处理技术领域，具体涉及一种基于Fourier级数模型的时序数据预测方法。

背景技术

卫星在轨工作后，长期运行在空间环境中，其性能与功能随着时间、测控事件及故障等因素发生变化，这些变化能够在遥测数据上体现出来，因此分析在轨卫星的遥测数据(时序数据)变化规律，对遥测数据未来变化趋势进行预测，在此基础上进行提前预警，可以在早期及时发现卫星性能的异常变化，并采取有效措施，避免可能发生的重大故障，降低卫星在轨运行的风险。

在航天领域，由于卫星设备结构、运行环境复杂、设备不可维修和回收分析等众多因素，造成对卫星在轨运行状态的预测，尤其是卫星长期状态的预测较为困难。现阶段，基于历史数据的卫星在轨运行状态的预测一般都是短期的，而且对于变化比较复杂的情况预测准确性较低。

发明内容

有鉴于此，本发明提供了一种基于Fourier级数模型的时序数据预测方法，能够对在轨卫星的遥测数据进行有效预测，且预测精度较高。

本发明的基于Fourier级数模型的时序数据预测方法，包括如下步骤：

步骤1，对遥测采样数据进行预处理，剔除粗大误差、剔野；

步骤2，建立Fourier级数模型；

其中，Fourier级数模型为：

y＝Φ(x,t)＝a₀+a₁cos(ωt)+b₁sin(ωt)+a₂cos(2ωt)+b₂sin(2ωt)+a₃cos(3ωt)+b₃sin(3ωt)+a₄cos(4ωt)+b₄sin(4ωt)+…+a_Kcos(Kωt)+b_Ksin(Kωt)

其中，x＝(a₀,a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄,…,a_K,b_K,ω)，为待求模型参数；ω为基频，(a₀,a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄,…,a_K,b_K)为模型系数；t为采样时间；y为采样时间对应的遥测采样数据；K为Fourier级数模型的阶数；

步骤3，求取Fourier级数模型的参数x；

步骤3.1，利用FFT变换求取步骤1预处理后的采样数据的最大频率ω_peak：

F_Y＝FFT(Y)

${\mathrm{\ω}}_{peak}=2\mathrm{\π}\frac{max(0.5,\max \_loc-1)}{t_m-t_1}$

其中，Y＝[y₁,y₂,…，y_m]^T，y₁,y₂,…，y_m为采样时刻t₁,t₂,…，t_m对应的遥测数据采样值；max_loc为FFT变换序列F_Y前一半的最大值索引号；

步骤3.2，令ω＝ω_peak/k，k＝1,2,…,K；针对K个ω，均进行如下处理：

(1)构造矩阵A：

$A={\left[\begin{array}{cccccccc}1& \cos \left(t_1\mathrm{\ω}\right)& \sin \left(t_1\mathrm{\ω}\right)& \cos \left(2t_1\mathrm{\ω}\right)& \sin \left(2t_1\mathrm{\ω}\right)& \mathrm{...}& \cos \left(4t_1\mathrm{\ω}\right)& \sin \left(4t_1\mathrm{\ω}\right)\\ 1& \cos \left(t_2\mathrm{\ω}\right)& \sin \left(t_2\mathrm{\ω}\right)& \cos \left(2t_2\mathrm{\ω}\right)& \sin \left(2t_2\mathrm{\ω}\right)& \mathrm{...}& \cos \left(4t_2\mathrm{\ω}\right)& \sin \left(4t_2\mathrm{\ω}\right)\\ .& & & & & & & \\ .& & & & & & & \\ .& & & & & & & \\ 1& \cos \left(t_m\mathrm{\ω}\right)& \sin \left(t_m\mathrm{\ω}\right)& \cos \left(2t_m\mathrm{\ω}\right)& \sin \left(2t_m\mathrm{\ω}\right)& \mathrm{...}& \cos \left(4t_m\mathrm{\ω}\right)& \sin \left(4t_m\mathrm{\ω}\right)\end{array}\right]}_{m\×9}$

(2)对矩阵A进行QR分解：

[Q,R,E]＝qr(A,0)

(3)求取矩阵R的秩rand(R)，令p＝rand(R)，取矩阵Q的前p列，即

Q′＝Q(:,1:p)

(4)求取拟合序列Y_fit＝Q′×(Q′^T×Y)；

(5)计算||Y-Y_fit||；

步骤3.3，针对步骤3.2获得的K个||Y-Y_fit||值，选取其中最小值||Y-Y_fit||_min对应的频率为最佳基频ω_bset；

步骤3.4，根据步骤3.3获得的最佳基频ω_bset，利用最小二乘求取其他模型系数：[a₀ a₁ b₁…a₄ b₄]^T＝(A^T A)^-1A^TY，其中，

$A={\left[\begin{array}{cccccc}1& \cos \left(t_1{\mathrm{\ω}}_{best}\right)& \sin \left(t_1{\mathrm{\ω}}_{best}\right)& \mathrm{...}& \cos \left(4t_1{\mathrm{\ω}}_{best}\right)& \sin \left(4t_1{\mathrm{\ω}}_{best}\right)\\ 1& \cos \left(t_2{\mathrm{\ω}}_{best}\right)& \sin \left(t_2{\mathrm{\ω}}_{best}\right)& \mathrm{...}& \cos \left(4t_2{\mathrm{\ω}}_{best}\right)& \sin \left(4t_2{\mathrm{\ω}}_{best}\right)\\ .& & & & & \\ .& & & & & \\ .& & & & & \\ 1& \cos \left(t_m{\mathrm{\ω}}_{best}\right)& \sin \left(t_m{\mathrm{\ω}}_{best}\right)& \mathrm{...}& \cos \left(4t_m{\mathrm{\ω}}_{best}\right)& \sin \left(4t_m{\mathrm{\ω}}_{best}\right)\end{array}\right]}_{m\×9}$

从而获得确定的Fourier级数模型；

步骤4，利用步骤3确定的Fourier级数模型对遥测数据进行预测。

较优的，Fourier级数模型的阶数为4。

有益效果：

本发明采用Fourier级数模型，对具有一定规则周期的遥测数据进行建模，并结合遥测数据特点，结合QR算法和最小二乘法确定模型的参数，最后利用确定的Fourier级数模型对遥测数据进行预测，获得较好的预测结果，与实际的遥测数据误差小，预测精度高。

附图说明

图1为本发明流程图。

图2为采用本发明方法的预测效果图。其中，图2(a)为预测信息显示，图2(b)为误差信息显示。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提供了一种基于Fourier级数模型的时序数据预测方法，采用Fourier级数模型，对具有一定规则周期的遥测数据进行建模并确定模型的参数，进而采用该模型对遥测数据进行预测。

具体包括如下步骤：

步骤1，对遥测采样数据进行预处理，剔除粗大误差、剔野，获得能真正反映卫星在轨运行状态的数据集合。

步骤2，建立Fourier级数模型。

Fourier级数模型能够较好地建模拟合各种有周期规律的数据类型，也是一种基于谐波频率的数据模型。本具体实施方式以4阶Fourier级数模型为例进行说明。该4阶Fourier级数模型如下：

y＝Φ(x,t)＝a₀+a₁cos(ωt)+b₁sin(ωt)+a₂cos(2ωt)+b₂sin(2ωt)+a₃cos(3ωt)+b₃sin(3ωt)+a₄cos(4ωt)+b₄sin(4ωt)

其中，x＝(a₀,a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄,ω)为待求模型参数，ω为基频，模型参数个数n＝10，t为采样时间，y为采样时间对应的遥测采样数据。

步骤3，求取Fourier级数模型的参数。

以4阶Fourier级数模型为例，说明模型参数x的求取思路：先利用FFT变换求取步骤1预处理后的采用数据的最大主频，将最大主频为最大基频，以最大基频依次倍除1,2,3,…,k，其中k为Fourier级数模型的阶数，获得的频率作为基频；然后利用QR分解测试出最佳基频；最后结合最佳基频，通过最小二乘求取其他模型系数，从而获得模型参数的初始值。具体包括如下子步骤。

步骤3.1，利用FFT变换求取步骤1预处理后的采样数据的最大频率ω_peak：

F_Y＝FFT(Y)

${\mathrm{\ω}}_{peak}=2\mathrm{\π}\frac{max(0.5,\max \_loc-1)}{t_m-t_1}$

其中，Y＝[y₁,y₂,…，y_m]^T，y₁,y₂,…，y_m为采样时刻t₁,t₂,…，t_m对应的遥测数据采样值；max_loc为FFT变换序列F_Y前一半的最大值索引号，记最佳基频ω_best为ω_best＝ω_peak。

步骤3.2，令k＝1,2,…,4，ω＝ω_peak/k，初始状态normr＝Inf，Inf为无穷大；循环完成以下测试：

(1)构造如下矩阵

$A={\left[\begin{array}{cccccccc}1& \cos \left(t_1\mathrm{\ω}\right)& \sin \left(t_1\mathrm{\ω}\right)& \cos \left(2t_1\mathrm{\ω}\right)& \sin \left(2t_1\mathrm{\ω}\right)& \mathrm{...}& \cos \left(4t_1\mathrm{\ω}\right)& \sin \left(4t_1\mathrm{\ω}\right)\\ 1& \cos \left(t_2\mathrm{\ω}\right)& \sin \left(t_2\mathrm{\ω}\right)& \cos \left(2t_2\mathrm{\ω}\right)& \sin \left(2t_2\mathrm{\ω}\right)& \mathrm{...}& \cos \left(4t_2\mathrm{\ω}\right)& \sin \left(4t_2\mathrm{\ω}\right)\\ .& & & & & & & \\ .& & & & & & & \\ .& & & & & & & \\ 1& \cos \left(t_m\mathrm{\ω}\right)& \sin \left(t_m\mathrm{\ω}\right)& \cos \left(2t_m\mathrm{\ω}\right)& \sin \left(2t_m\mathrm{\ω}\right)& \mathrm{...}& \cos \left(4t_m\mathrm{\ω}\right)& \sin \left(4t_m\mathrm{\ω}\right)\end{array}\right]}_{m\×9}$

其中，{t₁,t₂,…,t_m}为采样时刻序列。

(2)对A矩阵进行QR分解

[Q,R,E]＝qr(A,0)

(3)求取矩阵R的秩rand(R)，令p＝rand(R)，取矩阵Q的前p列，即

Q＝Q(:,1:p)

(4)求取拟合序列y_fit＝Q×(Q^T×y)；

(5)如果||y-y_fit||＜normr，则normr＝||y-y_fit||，且最佳基频ω_bset＝ω。

步骤3.3，根据最佳基频ω_bset再次利用最小二乘求取其他模型系数：

$A={\left[\begin{array}{cccccc}1& \cos \left(t_1{\mathrm{\ω}}_{best}\right)& \sin \left(t_1{\mathrm{\ω}}_{best}\right)& \mathrm{...}& \cos \left(4t_1{\mathrm{\ω}}_{best}\right)& \sin \left(4t_1{\mathrm{\ω}}_{best}\right)\\ 1& \cos \left(t_2{\mathrm{\ω}}_{best}\right)& \sin \left(t_2{\mathrm{\ω}}_{best}\right)& \mathrm{...}& \cos \left(4t_2{\mathrm{\ω}}_{best}\right)& \sin \left(4t_2{\mathrm{\ω}}_{best}\right)\\ .& & & & & \\ .& & & & & \\ .& & & & & \\ 1& \cos \left(t_m{\mathrm{\ω}}_{best}\right)& \sin \left(t_m{\mathrm{\ω}}_{best}\right)& \mathrm{...}& \cos \left(4t_m{\mathrm{\ω}}_{best}\right)& \sin \left(4t_m{\mathrm{\ω}}_{best}\right)\end{array}\right]}_{m\×9}$

A×[a₀ a₁ b₁…a₄ b₄]^T＝Y

则[a₀ a₁ b₁…a₄ b₄]^T＝(A^T A)^-1A^TY

步骤4，利用确定的Fourier级数模型对遥测数据进行预测。

其中，采用Fourier基数模型能够较好地对具有一定周期性规律和趋势性的数据进行建模，能够反应出数据的周期与趋势。

本发明的Fourier模型初始化算法综合应用了快速傅里叶变换与矩阵的QR分解算法，由于QR分解需要的矩阵空间非常大，因此考虑到实用效果，处理的速度受数据量大小的影响，步骤1中的遥测采样数据的数据量不要超过50000为宜。

下面结合具体数据进行说明。

根据用户选择的模型阶数，调用对应的数据预测方法对历史数据进行预测，并统计误差信息。

输入项

本方法软件实现后需要的输入数据如表1所示：

表1输入项

输出项

本方法用软件计算后的输出数据如下：

表2输出项

预测效果图如图2所示。其中，图2(a)为预测信息显示，其中，深灰色为历史的遥测采样数据；浅灰色为采用本发明方法预测得到的数据；图2(b)为误差信息显示，其中，浅灰色为本发明预测得到的数据相对于历史遥测数据的相对误差，深灰色为本发明预测得到的数据的绝对误差。

针对周期性遥测数据和周期性不明显的平稳遥测数据，采用本发明预测方法进行预测，结果如表3所示，可以看出，本发明预测方法的预测精度较高。

表3预测信息误差分析表

预测数据类型	平稳遥测数据	周期性遥测数据
			预测数据长度	10000	10000
最大绝对误差	7.99361XE-15	0.455594
			平均绝对误差	7.99361XE-15	0.145056
均方差	7.99361XE-15	0.099329
			平均误差	7.99361XE-15	0.178464
平均相对误差	1.70077XE-13％	2.31％

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于Fourier级数模型的时序数据预测方法，其特征在于，包括如下步骤：

步骤1，对遥测采样数据进行预处理，剔除粗大误差、剔野；

步骤2，建立Fourier级数模型；

其中，Fourier级数模型为：

y＝Φ(x,t)＝a₀+a₁cos(ωt)+b₁sin(ωt)+a₂cos(2ωt)+b₂sin(2ωt)+a₃cos(3ωt)+b₃sin(3ωt)

+a₄cos(4ωt)+b₄sin(4ωt)+…+a_Kcos(Kωt)+b_Ksin(Kωt)

其中，x＝(a₀,a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄,…,a_K,b_K,ω)，为待求模型参数；ω为基频，(a₀,a₁,b₁,a₂,b₂,a₃,b₃,a₄,b₄,…,a_K,b_K)为模型系数；t为采样时间；y为采样时间对应的遥测采样数据；K为Fourier级数模型的阶数；

步骤3，求取Fourier级数模型的参数x；

步骤3.1，利用FFT变换求取步骤1预处理后的采样数据的最大频率ω_peak：

F_Y＝FFT(Y)

${\mathrm{\ω}}_{peak}=2\mathrm{\π}\frac{max(0.5,\max \_loc-1)}{t_m-t_1}$

其中，Y＝[y₁,y₂,…，y_m]^T，y₁,y₂,…，y_m为采样时刻t₁,t₂,…，t_m对应的遥测数据采样值；max_loc为FFT变换序列F_Y前一半的最大值索引号；

步骤3.2，令ω＝ω_peak/k，k＝1,2,…,K；针对K个ω，均进行如下处理：

(1)构造矩阵A：

$A={\left[\begin{array}{cccccccc}1& \cos \left(t_1\mathrm{\ω}\right)& \sin \left(t_1\mathrm{\ω}\right)& \cos \left(2t_1\mathrm{\ω}\right)& \sin \left(2t_1\mathrm{\ω}\right)& ...& \cos \left(4t_1\mathrm{\ω}\right)& \sin \left(4t_1\mathrm{\ω}\right)\\ 1& \cos \left(t_2\mathrm{\ω}\right)& \sin \left(t_2\mathrm{\ω}\right)& \cos \left(2t_2\mathrm{\ω}\right)& \sin \left(2t_2\mathrm{\ω}\right)& ...& \cos \left(4t_2\mathrm{\ω}\right)& \sin \left(4t_2\mathrm{\ω}\right)\\ .& & & & & & & \\ .& & & & & & & \\ .& & & & & & & \\ 1& \cos \left(t_m\mathrm{\ω}\right)& \sin \left(t_m\mathrm{\ω}\right)& \cos \left(2t_m\mathrm{\ω}\right)& \sin \left(2t_m\mathrm{\ω}\right)& ...& \cos \left(4t_m\mathrm{\ω}\right)& \sin \left(4t_m\mathrm{\ω}\right)\end{array}\right]}_{m\×9}$

(2)对矩阵A进行QR分解：

[Q,R,E]＝qr(A,0)

(3)求取矩阵R的秩rand(R)，令p＝rand(R)，取矩阵Q的前p列，即

Q′＝Q(:,1:p)

(4)求取拟合序列Y_fit＝Q′×(Q′^T×Y)；

(5)计算||Y-Y_fit||；

步骤3.3，针对步骤3.2获得的K个||Y-Y_fit||值，选取其中最小值||Y-Y_fit||_min对应的频率为最佳基频ω_bset；

步骤3.4，根据步骤3.3获得的最佳基频ω_bset，利用最小二乘求取其他模型系数：[a₀ a₁b₁ … a₄ b₄]^T＝(A^TA)^-1A^TY，其中，

$A={\left[\begin{array}{cccccc}1& \cos \left(t_1{\mathrm{\ω}}_{best}\right)& \sin \left(t_1{\mathrm{\ω}}_{best}\right)& ...& \cos \left(4t_1{\mathrm{\ω}}_{best}\right)& \sin \left(4t_1{\mathrm{\ω}}_{best}\right)\\ 1& \cos \left(t_2{\mathrm{\ω}}_{best}\right)& \sin \left(t_2{\mathrm{\ω}}_{best}\right)& ...& \cos \left(4t_2{\mathrm{\ω}}_{best}\right)& \sin \left(4t_2{\mathrm{\ω}}_{best}\right)\\ .& & & & & \\ .& & & & & \\ .& & & & & \\ 1& \cos \left(t_m{\mathrm{\ω}}_{best}\right)& \sin \left(t_m{\mathrm{\ω}}_{best}\right)& ...& \cos \left(4t_m{\mathrm{\ω}}_{best}\right)& \sin \left(4t_m{\mathrm{\ω}}_{best}\right)\end{array}\right]}_{m\×9}$

从而获得确定的Fourier级数模型；

步骤4，利用步骤3确定的Fourier级数模型对遥测数据进行预测。

2.如权利要求1所述的基于Fourier级数模型的时序数据预测方法，其特征在于，Fourier级数模型的阶数为4。