一种自适应匹配的地球自转参数预报方法

摘要

本发明公开了一种自适应匹配的地球自转参数预报方法。使用本发明能够实现ERP预报所需的几类关键参数的自适应匹配，提高ERP预报精度。本发明采用LS+AR方法进行ERP预报，并给出了ERP预报时间序列训练长度经验设置准则；采用差分方法提高ERP残差的平稳性，并引入卡方检验方法，通过计算卡方值，定量评估PMX、PMY和UT1‑UTC残差的分布特性，并给出了卡方值与差分阶次的对应关系经验准则；最后在ERP残差的AR预报中，给出了AR模型的模型阶数确定方法，从而实现自适应匹配的ERP预报，并提高了预报精度。

说明书

技术领域

本发明涉及航天测控技术领域，具体涉及一种自适应匹配的地球自转参数预报方法。

背景技术

地球自转参数(ERP)是实现国际天球参考系(ICRS)与国际地球参考系(ITRS)转换的必须参数，是深空航天器跟踪导航、地球空间航天器测控等航天任务中不可缺少的重要状态参量。由于航天器定轨一般是在国际天球参考系中进行，而导航定位系统的最终目的是确定航天器在地球参考系中的位置和速度，发布的星历数据均处于国际地球参考系中，因此，需要进行ICRS与ITRS的相互转换。ERP用于实现ICRS与ITRS相互转换，其精度直接影响两参考系转换精度，进而对航天器的高精度导航、定位结果直接产生影响。此外，ERP作为大地测量学、天体测量学和地球物理学研究最主要的参数之一，与地球的核、幔、地壳、海洋、大气各圈层的物质运动及相互作用有密切的关系。

高精度实时导航定位任务离不开ERP数据支持，但由于从ERP观测过程到数据传输、处理过程需要花费大量时间，且数据处理过程较复杂，致使ERP产品发布具有一定的时间滞后。一方面为了准实时与实时导航任务要求，必须采用ERP预报产品用于支持任务实施，同时涉及航天器测控轨道超前任务规划等也需ERP预报产品支持；另一方面，当航天器进入自主定轨模式后，地面控制系统不能再上传最新的ERP数据来完成ITRS与ICRS的转换，而只能通过ERP预报产品来保证两者之间的转换。因此，开展ERP的预报方法研究对于航天器测控具有十分重要的作用。此外，ERP预报在地球物理现象的动力学分析方面具有重要应用价值。

目前为止，国内外不同学者与学术机构对ERP预报进行了深入研究，有多种预报方法应用于ERP的预报中，如最小二乘(LS)外推法、最小二乘联合自回归分析法(LS+AR)、神经网络预报法、最小二乘卡尔曼滤波预报方法、小波分解与自回归分析方法等。在众多ERP预报方法中，基于LS+AR预报方法被国际上公认为最有效的ERP预报方法之一，目前国际上ERP预报产品的主要提供机构，国际地球自转与参考架服务组织(IERS)与美国海军天文台(USNO)，均主要基于LS+AR方法进行ERP预报。

在传统的LS+AR的ERP预报方法中，还有以下问题：(1)ERP时间序列的平稳特性评估，本发明提出采用差分处理提高ERP时间序列平稳性，但是当前并没有明确的差分阶次的确定方法；(2)虽然AR模型对平稳时间序列具有很好的预报特性，但是当前并没有明确的AR模型的最优阶数确定方法，且不同ERP特性对应不同的AR模型阶数，而AR模型的最优阶数确定在AR模型中尤为关键；(3)ERP预报时需要利用一定长度的ERP时间序列用于训练，从而实现ERP外推预报，用于训练的ERP长度与ERP预报精度直接相关，也是ERP预报的关键问题之一，但当前并没有相关的ERP时间序列训练长度的选取方法。

发明内容

有鉴于此，本发明提供了一种自适应匹配的地球自转参数预报方法，能够实现ERP预报所需的几类关键参数的自适应匹配，提高ERP预报精度。

本发明的自适应匹配的地球自转参数预报方法，包括如下步骤：

步骤1，分别选取给定训练长度的PMX、PMY和UT1-UTC数据；其中，选取的PMX、PMY和UT1-UTC统称为原始ERP数据；

步骤2，将步骤1获得的三种原始ERP数据分别进行周期函数与二次多项式的最小二乘拟合，获得ERP拟合值；用步骤1的原始ERP数据减去ERP拟合值，获取ERP残差；

步骤3，利用原始ERP数据的周期函数与二次多项式函数对步骤2的ERP拟合值进行外推，其中，ERP外推长度与期望的预报长度一致；

步骤4，对步骤2获得的ERP残差进行卡方值计算，根据ERP残差的卡方值χ²和式(1)确定差分阶次：

步骤5，根据步骤4确定的差分阶次，对ERP残差进行差分运算，获得差分后的ERP残差；

步骤6，针对步骤5获得的差分后的ERP残差，采用AR模型进行ERP残差的AR预报，获得ERP残差的预报值；

步骤7，将步骤6获得ERP残差的预报值进行逆差分处理，并将逆差分处理结果与步骤3获得的外推结果进行叠加，即获取ERP的预报初值；

步骤8，步骤7获得的PMX预报初值和PMY预报初值即分别为PMX和PMY的最终预报值；UT1-UTC预报初值经过带谐固体潮汐项预报添加与跳秒恢复后获得最终的UT1-UTC预报值。

进一步地，所述步骤1中，依照下式选取PMX、PMY和UT1-UTC数据的训练长度：

$\left(\begin{array}{c}3\le length(PMX\_N)\le 4\\ 3\le length(PMY\_N)\le 4\\ 19\le length(UT1-UTC\_N)\le 22\end{array}\right)---\left(2\right)$

式中的单位为年，length(PMX_N)为PMX的训练长度，length(PMY_N)为PMY的训练长度，length(UT1-UTC_N)为UT1-UTC的训练长度。

进一步地，所述步骤6中，所述AR模型的模型阶数确定方法如下：

步骤6.1，AR模型为：

其中，为AR模型参数；a_t为白噪声；p为AR模型阶数；x_t,x_t-1,x_t-2,…,x_t-p为t,t-1,t-2,…,t-p时刻的平稳时间序列的状态量；

设置AR模型阶数p的搜索空间为[1,P]；将模型阶数p的P个取值分别代入式(4)中，获得P个AR模型；

步骤6.2，针对步骤6.1获得的P个AR模型，采用解Yule-Walker方程法估计每一个模型阶数p对应的模型参数然后将估计出的模型参数代回式(4)中，计算获得平稳时间序列的模型值，然后用平稳时间序列的状态值减去模型值得到噪声残差序列；

步骤6.3，采用终预误差准则确定每一个AR阶数p所对应的判断准则系数ξ_p：

${\xi }_p=FPE\left(p\right)=\frac{N+p}{N-p}{\sigma }_n^2---\left(5\right)$

其中，N为样本长度，为噪声残差序列的方差；

步骤6.4，对P个AR模型阶数p对应的判断准则系数进行判断，最小判断准则系数对应的模型阶数即为确定的AR模型的AR模型阶数。

有益效果：

本发明提出ERP预报差分经验准则、自适应匹配的AR阶数搜索确定方法和预报时间序列训练长度经验准则，联合多项式与周期函数的最小二乘平滑方法、卡方检验方法与AR预报方法，综合确定ERP预报结果，并评定ERP预报精度，实现自适应匹配的ERP预报。

附图说明

图1为自适应匹配的ERP预报数据处理流程图。

图2为用于预报训练的PMX原始时域波形。

图3为PMX的最小二乘拟合与残差图。

图4为用于预报训练的PMY原始时域波形。

图5为PMY的最小二乘拟合与残差图。

图6为UT1-UTC的跳秒剔除与待谐固体潮汐项消除。

图7为UT1-UTC的最小二乘拟合与残差图。

图8为PMX的预报结果与国际比对。

图9为PMY的预报结果与国际比对。

图10为UT1-UTC的预报结果与国际比对。

图11为ERP预报误差统计图。

具体实施方式

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

本发明提供了一种自适应匹配的地球自转参数预报方法，给出了LS+AR的ERP预报方法中关键参数的确定方法：

(1)为了能利用AR模型方法实现ERP预报，待处理的ERP时间序列需满足平稳性要求。对于ERP时间序列，原始的PMX(极移X分量)、PMY(极移Y分量)和UT1-UTC(世界时与世界协调时差值)时间序列均不满足平稳性要求，一方面PMX、PMY和UT1-UTC时间序列中存在已知的周期分量，这些周期分量在进行AR预报之前必须提前予以扣除，得到PMX、PMY和UT1-UTC残差；然后，可通过差分的方式，提供PMX、PMY和UT1-UTC残差的平稳性。为了选用合理差分阶次，同时避免差分阶次欠缺或者过度差分，需根据不同的待分析的PMX、PMY和UT1-UTC残差的平稳性，选用合适的差分阶次。本发明引入卡方检验方法，通过计算卡方值，定量评估PMX、PMY和UT1-UTC残差的分布特性，并提出了卡方值与差分阶次的对应关系准则(式(1))，从而实现PMX、PMY和UT1-UTC的差分阶次自适应匹配。

(2)为了实现AR模型阶数的自适应匹配，本发明提出自适应匹配的AR阶数搜索确定方法。针对AR模型，在一个指定的AR阶数搜索空间内，采用解Yule-Walker方程法用于估计每一个模型阶数p对应的AR模型参数，然后利用终预误差准则(FPE)进行判断，选取最优的AR模型阶数。

(3)为了选用合理的ERP预报时间序列训练长度，实现ERP的有效预报，本发明在大量试验与分析基础上，给出了ERP预报时间训练长度经验准则如下：

$\left(\begin{array}{c}3\le length(PMX\_N)\le 4\\ 3\le length(PMY\_N)\le 4\\ 19\le length(UT1-UTC\_N)\le 22\end{array}\right)---\left(2\right)$

式中的单位为年(year)，PMX_N为PMX的训练长度，PMY_N为PMY的训练长度，UT1-UTC_N为UT1-UTC的训练长度。

下面结合具体实施例进行说明。

本发明用于预报的ERP原始数据来自IERS网站的EOP 08C04加上1个月的USNO的ERP解算值。本发明以从2016年2月19日(对应的MJD为57437)开始预报为例，往后预报90天的ERP值。

本发明的自适应匹配的ERP预报方法的具体数据处理过程如下：

步骤1，依据本发明中提出的ERP预报时间序列训练长度经验准则(式(2))，确定PMX/PMY/UT1-UTC的训练长度。

按照本发明提出的ERP预报时间序列训练长度经验准则设置PMX、PMY的训练长度为1300天(3.56年)，UT1-UTC的训练长度为7700天(21.08年)。

步骤2，进行ERP的最小二乘拟合，获取ERP残差。

将步骤1选取的PMX、PMY依据其已知周期函数与二次多项式直接进行最小二乘拟合获得拟合值；UT1-UTC经过跳秒检测、跳秒消除与带谐固体潮汐项剔除后依据其已知周期函数与二次多项式直接进行最小二乘拟合，获得拟合值。用原始ERP数据减去ERP拟合值，获取ERP残差。

PMX的原始波形如图2所示，PMX经最小二乘拟合后的残差图如图3所示。PMY的原始波形如图4所示，PMY经最小二乘拟合后的残差图如图5所示。UT1-UTC的原始波形经过跳秒检测与剔除、带谐固体潮汐项剔除后的图形如图6所示，UT1-UTC经最小二乘拟合后的残差图如图7所示。

步骤3，利用已知的周期函数与二次多项式函数进行ERP外推。

利用PMX、PMY、UT1-UTC已知的周期函数与二次多项式函数分别对PMX、PMY及剔除跳秒与带谐固体潮汐项的UT1-UTC进行外推预报。PMX，PMY的已知周期项为钱德勒周期项(435天)、1年周期项、1/2年周期项和1/3年周期项；UT1-UTC的已知周期项为18.6年、9.3年、1年和1/2年；ERP的二次多项式系数通过最小二乘拟合方法确定。ERP外推长度与期望的预报长度一致，可事先设定，本发明中的外推长度为90天，即三个月。

步骤4，引入卡方检验方法计算ERP残差的卡方值。

为了能利用AR模型方法实现ERP预报，待处理的ERP时间序列需满足平稳性要求。对于ERP时间序列，原始的PMX、PMY和UT1-UTC时间序列均不满足平稳性要求，一方面PMX、PMY和UT1-UTC时间序列中存在已知的周期分量，这些周期分量在进行AR预报之前必须提前予以扣除，得到PMX、PMY和UT1-UTC残差；另一方面，可通过差分的方式，提供PMX、PMY和UT1-UTC残差的平稳性，为了选用合理差分阶次，同时避免差分阶次欠缺或者过度，需根据不同待分析的PMX、PMY和UT1-UTC残差的平稳性，选用合适的差分阶次。本发明中通过引入卡方检验方法，通过计算卡方值，定量评估PMX、PMY和UT1-UTC残差的分布特性。其中卡方值的计算公式如式(3)所示。

${\chi }^2=\Sigma \frac{{(A-E)}^2}{E}=\underset{i=1}{\overset{k}{\Sigma }}\frac{{(A_i-E_i)}^2}{E_i}=\underset{i=1}{\overset{k}{\Sigma }}\frac{{(A_i-{\mathrm{np}}_i)}^2}{{\mathrm{np}}_i}---\left(3\right)$

其中，A_i为i水平的观测频数，E_i为i水平的期望频数，n为总频数，p_i为i水平的期望频率。i水平的期望频数E_i等于总频数n×i水平的期望概率p_i，k为单元格数。当n比较大时，χ²统计量近似服从k-1(计算E_i时用到的参数个数)个自由度的卡方分布。

基于以上方法，分别对PMX、PMY、UT1-UTC的拟合残差进行卡方值计算，其结果为χ²_pmx＝0.41，χ²_pmy＝0.36，χ²_ut1-utc＝0.83。

步骤6，提出ERP差分经验准则，确定ERP残差的差分阶次，并进行ERP残差的差分运算。

本发明提出卡方值与差分阶次的对应关系经验准则，从而实现PMX、PMY和UT1-UTC的差分阶次自适应匹配。本发明提出ERP预报差分经验准则，用于确定对ERP序列进行差分的阶次。ERP预报差分经验准则如式(1)所示：

因此，依据本发明提出的ERP差分经验准则，对PMX残差经过1次差分，对PMY残差经过1次差分，对UT1-UTC残差经过2次差分。

步骤7，提出AR模型阶数自适应匹配搜索确定方法，获取AR模型阶数，进行ERP残差的AR预报。

为了实现AR模型阶数的自适应匹配，本发明提出自适应匹配的AR阶数搜索确定方法。其基本思想是，在一个指定的AR阶数区间内对AR阶数进行最优搜索，并按照判断准则找出最优的AR模型阶数。具体实现方法如下：

首先介绍AR模型的数学表达式如式(4)所示。

式中为AR模型参数；a_t为白噪声；p为AR模型阶数；x_t,x_t-1,x_t-2,…,x_t-p为t,t-1,t-2,…,t-p时刻的平稳时间序列(即差分后的ERP残差)的状态量；称式(4)为p阶自回归模型，简记为AR(p)；a_t满足正态分布，a_t～N(0,σ²)，σ²为白噪声的方差。

从式(4)中可以看出AR模型阶数直接影响AR预报结果。

然后介绍本发明的AR模型阶数确定的方法步骤为：(a)设置AR模型阶数的搜索空间为1至200，即p＝1:200，这时将p分别带入式(4)中，获得200个AR模型；(b)为了解算式(4)的AR模型参数这里采用解Yule-Walker方程法估计每1个模型阶数p对应的模型参数然后将估计出的模型参数带入式(4)中，计算获得平稳时间序列的模型值，然后用平稳时间序列的状态值减去平稳时间序列的模型值，得到噪声残差序列；(c)采用终预误差准则(FPE)，如式(5)所示，用于确定每一个AR阶数p所对应的判断准则系数ξ_p；(d)比较ξ_p,p＝1:200，当ξ(min)＝ξ_p的p即为最优AR模型阶数。

${\xi }_p=FPE\left(p\right)=\frac{N+p}{N-p}{\sigma }_n^2---\left(5\right)$

式中N为样本长度，p为模型阶次，为噪声残差序列的方差。

对ERP残差的差分值进行AR预报，其中AR模型阶数按照本发明提出的自适应匹配AR阶数搜索确定方法予以获取。PMX残差预报的AR阶数为54，PMY残差预报的AR阶数为54，UT1-UTC残差预报的AR阶数为82。

步骤8，ERP残差预报的逆差分处理。

对ERP残差的预报值进行逆差分处理，逆差分的阶次等于上面步骤5中的差分阶次。

步骤9，联合最小二乘外推与AR预报预报值，获取ERP预报初值。

ERP残差逆差分预报结果结合步骤3中的ERP周期函数与二次多项式函数外推结果，两者相加，获取ERP的预报初值。

步骤10，ERP预报最终值确定。

PMX、PMY此时的预报结果即为其最终预报结果，UT1-UTC预报初值经过带谐固体潮汐项预报添加与跳秒恢复后获得UT1-UTC预报值，基于此最终获得ERP预报值。

为了直观说明本发明的ERP预报效果，将本发明获取的ERP预报值与国际地球自转与参考架服务组织(IERS)、美国海军天文台(USNO)的ERP预报值进行国际比对，其中IERS的预报跨度为180天，USNO的预报跨度为90天，本发明的预报天数为90天，在图中显示的代号为BACC(北京航天飞行控制中心的缩写)。图8为PMX预报结果与国际比对结果，图9为PMY预报结果与国际比对结果，图10为UT1-UTC预报结果与国际比对结果。

从图8～图10显示的ERP(包含PMX、PMY和UT1-UTC)的预报结果直观可以看出，本发明的PMX、PMY预报结果与IERS、USNO的一致性较好，其中在图10中，由于本发明中所用于ERP预报的最近一个月的UT1-UTC值是利用USNO的解算值，所以UT1-UTC预报结果与USNO的一致性较好。

步骤11，利用平均绝对误差(MAE)准则评估ERP预报精度。

为了定量化描述本发明的ERP预报精度，利用MAE准则，如式(6)所示，这里统计了利用本发明方法计算的2015年全年(共365天)的ERP预报误差情况，统计的预报跨度为30天。其预报误差如图11所示，具体预报误差数值如表1所示，从中可以看出本发明的ERP预报具有很好的精度，这个ERP预报精度达到国际先进水平。

${\mathrm{MAE}}_i=\frac{1}{n}\underset{j=1}{\overset{n}{\Sigma }}\left(\right|{p_j}^i-{o_j}^i\left|\right)---\left(6\right)$

式中o为实际观测值；p为预报值；i为预报跨度；n为预报期数。

表1ERP预报误差统计表

Pre DayPMX Pre Error(mas)PMY Pre Error(mas)UT1-UTC Pre Error(ms)10.2050.1790.04651.710.9670.291102.281.520.818202.352.472.13303.143.183.58

本发明针对以上三个ERP预报中的问题，提出问题解决方法，并综合数据处理方法，实现ERP的高精度预报。

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