CN104090280B - 一种基于区域cors的电离层延迟修正预报方法 - Google Patents

Abstract

本发明公开了一种基于区域CORS的电离层延迟修正预报方法，包括以下步骤：电离层延迟初值解算，电离层延迟修正值解算，电离层延迟的解算，电离层延迟预报，电离层延迟修正预报。本发明可根据区域CORS的观测数据，快速地获取电离层延迟，并对电离层延迟进行精化，提高电离层延迟解算的精度。利用Holt-winter模型进行初步预报可以较好地描述电离层延迟的周日性变化，利用Markov模型进行电离层延迟修正预报，可以很好地描述电离层延迟的扰动，组合两种模型可以达到提高电离层延迟预报精度的目的，增加预报结果的可靠性和稳定性。

Description

一种基于区域CORS的电离层延迟修正预报方法

技术领域

本发明属于卫星导航定位技术领域，具体涉及一种基于区域CORS的电离层延迟修正预报方法。

背景技术

以全球定位系统(Global Positioning System，GPS)技术为代表的卫星空间定位技术因其全天候、自动化、高精度等特点日益成为导航定位、精密工程测量、地壳变形等领域的主要手段。电离层延迟已成为GPS定位中最大的误差源，不仅严重降低了GPS的定位精度，也限制了其应用领域和定位方法的可选择性。由于电离层的结构与变化规律极为复杂，目前的理论与方法更多地表现在从定性的角度认识和研究电离层的结构与活动特性，所得出的用于修正无线电波传播中电离层延迟影响的各种理论模型，不仅改正效果有限，而且数学描述复杂，难以满足准实时、实时的定位需求。随着研究和应用的不断发展，对电离层延迟模型的精度和准确性的要求也不断提高，如何进一步提高模型的精度，成为当前研究的热点。建立和发展大规模的连续运行参考站(Continuously Operating ReferenceStations,CORS)，已成为GPS应用技术新的发展趋势。CORS和电离层的研究是相辅相成的，一方面，不同尺度、密度的CORS参考站网为研究全球性或区域特性的电离层提供了技术和数据平台，可利用CORS进行精密解算、分析全球或区域电离层的结构及其分布特性；另一方面，利用电离层的研究成果建立更高精度的电离层延迟模型，有助于提高CORS的定位精度。

多项式函数建立电离层延迟模型是由函数的局部性质决定它的整体性质，其分时间段建模存在不连续问题，拟合曲面的光滑度和逼近精确度之间存在矛盾。此外，由于多项式函数模型本身是由短时间序列的观测数据拟合得到的，模型在时域上的可扩展性较差，无法实现对电离层延迟的精确预报。

发明内容

本发明的目的是为了消减电离层延迟对卫星导航定位精度的影响，提供一种基于区域CORS的电离层延迟建模和修正预报方法。该方法在利用B样条函数模型精化电离层延迟的基础上，提出了含有更新过程的电离层延迟预报模型，进一步提高了电离层延迟建模和预报的精度。

为实现上述目的，本发明提供了一种基于区域CORS的电离层延迟修正预报方法，包括：

(1)电离层延迟概略初值解算：

利用电离层单层模型和载波相位解算出天顶方向的总电子含量(VerticalTotal Electron Content,VTEC)：

$\mathrm{VTEC}=\frac{\cos z^{\′}}{40.28}\·\frac{{f_1}^2{f_2}^2}{{f_1}^2-{f_2}^2}[({\mathrm{\λ}}_2{\mathrm{\φ}}_2-{\mathrm{\λ}}_1{\mathrm{\φ}}_1)+({\mathrm{\λ}}_2{N}_2-{\mathrm{\λ}}_1{N}_1)-(b_2^S-b_1^S)+(b_2^R-b_1^R)]$

其中，z′为穿刺点处的天顶距，即地心至穿刺点方向与接收机至卫星方向的夹角，和分别为载波L₁和L₂相位观测的卫星硬件延迟，和分别载波L₁和L₂的接收机硬件延迟，f₁和f₂分别为载波L₁和L₂相位的频率，φ₁和φ₂分别为载波L₁和L₂的相位观测量，λ₁和λ₂分别为载波L₁和L₂相位的波长，N₁和N₂分别为载波L₁和L₂的整周模糊度，以周为单位。

多项式函数模型将VTEC看作是纬度差和太阳时角差S-S₀的函数，用一个规则的曲面来拟合各穿刺点处的VTEC值，多项式函数模型解算VTEC的表达式如下：

S-S₀＝(λ-λ₀)+(t-t₀)

其中，n_max、m_max为多项式函数模型的泰勒展开式中的最大维数，通常n_max取值为1-2，m_max取值为2-4，VTEC多项式函数模型对单层模型的高度H在计算中一般取值为300-500km。E_nm为多项式函数模型的系数，为测区中心点的地理纬度，为穿刺点或者星下点的地理纬度，S₀为测区中心点在该时段中央时刻t₀时的太阳时角，S为穿刺点或者星下点的太阳时角，t为观测时刻(世界时)，λ为穿刺点或者星下点的地理经度，λ₀为测区中心点的地理经度；

结合单层模型和载波相位解算的VTEC和多项式函数模型解算的VTEC，可以得到多项式函数模型的观测方程的最终形式：

利用多项式函数模型的观测方程的最终形式，解算电离层延迟概略初值；

(2)电离层延迟修正值解算：

根据多项式函数模型解算的误差，利用B样条函数解算电离层延迟的修正值其模型表达式为：

其中，mj_i＝2^ji+2，j₁、j₂为阶数，通常取值为1-2；为模型系数；为二维尺度函数。

(3)电离层延迟的解算：

根据已解算的电离层延迟的初值和利用B样条函数获取的电离层延迟修正值的模型，拟合出更为精确的电离层延迟：

(4)电离层延迟概略预报：

根据步骤(3)获得的电离层延迟，利用Holt-winter(HW)模型进行预报，HW模型的基本公式为：

稳定项： $S_t=\mathrm{\α}\left(\frac{X_t}{I_{t-L}}\right)+(1-\mathrm{\α})(S_{t-1}-b_{t-1})$

趋势项：b_t＝γ(S_t-S_t-1)+(1-γ)b_t-1

季节项： $I_t=\mathrm{\β}\left(\frac{X_t}{S_t}\right)+(1-\mathrm{\β})I_{t-L}$

预报：X_t+τ＝(S_t+τb_t)I_t-L+τ

其中，X_t为t时刻的观测值，S_t为t时刻的稳定项，I_t为t时刻的季节项，b_t为t时刻的趋势项，X_t+τ为τ期的预测值，τ为预测的期数，L为季节长度，α,β,γ∈[0,1]为平滑参数；

(5)电离层延迟修正预报：

利用Markov模型对HW模型的预报值进行更新修正，修正模型表达式为：

${\hat{X}}_t=X_t+{\hat{F}}_t$

其中，为修正后的电离层延迟预报值，X_t为HW模型的预报值，为残差预报值；

残差预报值的计算公式为：

${\hat{F}}_{t+k}=X_t\×E_t\×P^{\left(k\right)}$

其中，X_t为t时刻的电离层延迟，E_t为t时刻的Markov链的数据序列状态(以E₁,E₂,…,E_n表示，对应转移时间t₁,t₂,…,t_n)，P^(k)为k步状态转移概率矩阵，其计算公式为：

$P^{\left(k\right)}=\left[\begin{array}{ccccc}{P}_{11}^k& P_{12}^k& ...& P_{1(j-1)}^k& P_{1j}^k\\ P_{21}^k& P_{22}^k& ...& P_{2(j-1)}^k& P_{2j}^k\\ P_{31}^k& P_{32}^k& ...& P_{3(j-1)}^k& P_{3j}^k\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ P_{j1}^k& P_{j2}^k& ...& P_{j(j-1)}^k& P_{\mathrm{jj}}^k\end{array}\right]$

其中，表示数列由状态E_i经过k步变为E_j的概率，即：

$P_{\mathrm{ij}}^{\left(k\right)}=n_{\mathrm{ij}}^{\left(k\right)}/N_i$

其中，表示状态E_i经过k步变为E_j的次数，N_i表示状态E_i出现的总次数。

本发明根据区域CORS的观测数据，可快速地获取电离层延迟，并对电离层延迟进行精化，提高了电离层延迟解算的精度。HW模型进行初步预报可以较好地描述电离层延迟的周日性变化，Markov模型进行电离层延迟修正预报，可以很好地描述电离层延迟的扰动，组合两种模型可以达到提高电离层延迟预报精度的目的，增加预报结果的可靠性和稳定性。

附图说明

图1为本发明方法的流程图。

图2为电离层单层模型。

图3为电离层延迟概略初值解算结果(多项式函数模型)。

图4为电离层延迟修正值解算结果(B样条函数模型)。

图5为电离层延迟的解算结果。

图6为BTJR站电离层延迟多站解数据序列。

图7为BTJR站电离层延迟多站解与HW模型预报值的残差。

图8为HW模型残差与经Markov模型修正残差(DOY074)。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

一种基于区域CORS的电离层延迟修正预报方法，利用电离层延迟修正预报方法对区域CORS的电离层延迟进行预报，如图1所示，包括以下步骤：

(1)电离层延迟概略初值解算

利用电离层单层模型和载波相位解算出天顶方向的总电子含量(VerticalTotal Electron Content,VTEC)：

$\mathrm{VTEC}=\frac{\cos z^{\′}}{40.28}\·\frac{{f_1}^2{f_2}^2}{{f_1}^2-{f_2}^2}[({\mathrm{\λ}}_2{\mathrm{\φ}}_2-{\mathrm{\λ}}_1{\mathrm{\φ}}_1)+({\mathrm{\λ}}_2{N}_2-{\mathrm{\λ}}_1{N}_1)-(b_2^S-b_1^S)+(b_2^R-b_1^R)]$

其中，z′为穿刺点处的天顶距，即地心至穿刺点方向与接收机至卫星方向的夹角，和分别为载波L₁和L₂相位观测的卫星硬件延迟，和分别载波L₁和L₂的接收机硬件延迟，f₁和f₂分别为载波L₁和L₂相位的频率，φ₁和φ₂分别为载波L₁和L₂的相位观测量，λ₁和λ₂分别为载波L₁和L₂相位的波长，N₁和N₂分别为载波L₁和L₂的整周模糊度，以周为单位。

多项式函数模型将VTEC看作是纬度差和太阳时角差S-S₀的函数，用一个规则的曲面来拟合各穿刺点处的VTEC值，多项式函数模型解算VTEC的表达式如下：

其中，n_max、m_max为多项式函数模型的泰勒展开式中的最大维数，E_nm为多项式函数模型的系数，为测区中心点的地理纬度，为穿刺点或者星下点的地理纬度，S₀为测区中心点在该时段中央时刻t₀时的太阳时角，S为穿刺点或者星下点的太阳时角，t为观测时刻(世界时)，λ为穿刺点或者星下点的地理经度，λ₀为测区中心点的地理经度。

结合单层模型和载波相位解算的VTEC和多项式函数模型解算的VTEC，可以得到多项式函数模型的观测方程的最终形式：

9.52437[(λ₂φ₂-λ₁φ₁)+(λ₂N₂-λ₁N₁)]cosz′

利用多项式函数模型的观测方程的最终形式，解算电离层延迟概略初值。

电离层单层模型如图2所示。采用江苏省连续运行参考站网络(JiangsuContinuously Operating Reference Stations,JSCORS)的NJKF站和NJLS站，2009年年积日001的数据(15s采样间隔)，利用多项式函数模型的观测方程的最终形式以单站模式进行2小时时间间隔的电离层延迟计算，其中取n_max＝2，m_max＝3，H＝350m，解算结果如图3所示。由图3可知，利用多项式函数模型解算的电离层延迟具有明显的时段性，表现为2h的计算时段内电离层延迟比较平滑，而时段间存在明显的不连续性，整体建模精度有待提高。

(2)电离层延迟修正值解算

根据多项式函数模型解算的误差，利用B样条函数解算电离层延迟的修正值其模型表达式为：

其中，mj_i＝2^ji+2，j₁、j₂为阶数；为模型系数；为二维尺度函数。选择j₁＝j₂＝2，方程有36个待估参数，通过2小时的观测数据即可建立对应的模型。

数据选择同步骤(1)，利用B样条函数模型以单站模式进行2小时时间间隔的电离层延迟计算，解算结果如图4所示。利用B样条函数对电离层延迟的误差进行修正，解决了由于电离层延迟分时段解算的方式造成部分时刻误差较大和分布不均的问题，可以提高电离层高延迟解算的整体精度。

(3)电离层延迟的解算

根据已解算的电离层延迟的初值和电离层延迟的修正值得到电离层延迟：

数据选择同步骤(1)，结合多项式函数模型和B样条函数模型，解算结果如图5所示。经过B样条函数修正后的电离层延迟可以使各时段间的不连续性得到改善，电离层延迟的整体精度有所提高。这是由于B样条函数不仅保持了多项式函数模型的简单性和逼近的可行性，而且在各段之间保持了相对独立的局部性质，可以解决拟合曲面的光滑度和逼近精确度之间存在着矛盾。

(4)电离层延迟概略预报

根据步骤(3)获得的电离层延迟，利用Holt-winter(HW)模型进行预报，HW模型的基本公式为：

稳定项： $S_t=\mathrm{\α}\left(\frac{X_t}{I_{t-L}}\right)+(1-\mathrm{\α})(S_{t-1}-b_{t-1})$

趋势项：b_t＝γ(S_t-S_t-1)+(1-γ)b_t-1

季节项： $I_t=\mathrm{\β}\left(\frac{X_t}{S_t}\right)+(1-\mathrm{\β})I_{t-L}$

预报：X_t+τ＝(S_t+τb_t)I_t-L+τ

其中，X_t为t时刻的观测值，S_t为t时刻的稳定项，I_t为t时刻的季节项，b_t为t时刻的趋势项，X_t+τ为τ期的预测值，τ为预测的期数，L为季节长度，α,β,γ∈[0,1]为平滑参数，α,β,γ取值的原则是预测值和实算值之间的均方差最小。

采用经步骤(1)、(2)、(3)解算的JSCORS的BTJR站电离层延迟，参与预报的数据为8天(2009年年积日066～073)用以预报2009年年积日074的电离层延迟，采用电离层延迟IDPB模型计算的2小时时间间隔的中央时刻的多站解，其实算电离层延迟见图6。采用HW模型进行预报，其实算值与HW模型的预报值之差见图7。HW模型在用于短期的预报时各时段上的实算值与预报值之差的绝对值小于1TECU的约占45％，1-2TECU的约占30％。HW模型在以单日(12个电离层延迟)为周期进行预报时能较好地描述电离层延迟的周日性，可以用于短期的电离层延迟预报。但电离层延迟多站解与HW模型预报值的残差仍存在较大的波动，尤其是拟合精度较差的时段，HW模型未能很好地描述，需要对残差进行修正。

(5)电离层延迟修正预报

利用Markov模型对HW模型的预报值进行更新修正，修正模型的表达式为：

${\hat{X}}_t=X_t+{\hat{F}}_t$

其中，为修正后的电离层延迟预报值，X_t为HW模型的预报值，为Markov模型的残差预报值；

残差预报值的计算公式为：

${\hat{F}}_{t+k}=X_t\×E_t\×P^{\left(k\right)}$

其中，X_t为t时刻的电离层延迟，E_t为t时刻的Markov链的数据序列状态(以E₁,E₂,…,E_n表示，对应转移时间t₁,t₂,…,t_n)，P^(k)为k步状态转移概率矩阵，其计算公式为：

$P^{\left(k\right)}=\left[\begin{array}{ccccc}{P}_{11}^k& P_{12}^k& ...& P_{1(j-1)}^k& P_{1j}^k\\ P_{21}^k& P_{22}^k& ...& P_{2(j-1)}^k& P_{2j}^k\\ P_{31}^k& P_{32}^k& ...& P_{3(j-1)}^k& P_{3j}^k\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ P_{j1}^k& P_{j2}^k& ...& P_{j(j-1)}^k& P_{\mathrm{jj}}^k\end{array}\right]$

其中，表示数列由状态E_i经过k步变为E_j的概率，即：

$P_{\mathrm{ij}}^{\left(k\right)}=n_{\mathrm{ij}}^{\left(k\right)}/N_i$

其中，表示状态E_i经过k步变为E_j的次数，N_i表示状E_i出现的总次数。

采用Markov模型对步骤(4)的电离层延迟预报值残差进行修正，HW模型残差与经Markov模型修正的残差预报量见图8。利用Markov模型对HW模型的预报值进行更新修正，预报时各时段上的实算值与预报值之差的绝对值小于1TECU的约占50％，1-2TECU的约占38％，可以有效地对电离层延迟的预报值进行更新修正。

根据以上具体实施例的综合分析可知，本发明将电离层延迟分为概略初值和修正值分项解算的方法，可以提高电离层延迟解算的精度，为电离层延迟的预报提供更为精确的基础数据。利用Markov模型对HW模型预报的电离层延迟进行更新修正，分别利用了HW模型用于时序分析能够突出表达电离层延迟的周期性，Markov模型具有适合描述随机波动性较大的预测问题的优势。将两种预测模型组合运用，具有优势互补性，可以达到提高预测精度的目的，同时增加预测结果的可靠性和稳定性。

Claims (1)

1.一种基于区域CORS的电离层延迟修正预报方法，其特征在于：包括以下步骤：

(1)电离层延迟概略初值解算：

利用电离层单层模型和载波相位解算出天顶方向的总电子含量(Vertical TotalElectron Content,VTEC)：

$\mathrm{VTEC}=\frac{\cos z^{\′}}{40.28}\·\frac{{f_1}^2{f}_2^2}{{f_1}^2{f}_2^2}[({\mathrm{\λ}}_2{\mathrm{\φ}}_2-{\mathrm{\λ}}_1{\mathrm{\φ}}_1)+({\mathrm{\λ}}_2{N}_2-{\mathrm{\λ}}_1{N}_1)-(b_2^S-b_1^S)+(b_2^R-b_1^R)]$

其中，z′为穿刺点处的天顶距，即地心至穿刺点方向与接收机至卫星方向的夹角，和分别为载波L₁和L₂相位观测的卫星硬件延迟，和分别载波L₁和L₂的接收机硬件延迟，f₁和f₂分别为载波L₁和L₂相位的频率，φ₁和φ₂分别为载波L₁和L₂的相位观测量，λ₁和λ₂分别为载波L₁和L₂相位的波长，N₁和N₂分别为载波L₁和L₂的整周模糊度，以周为单位；

多项式函数模型将VTEC看作是纬度差和太阳时角差S-S₀的函数，用一个规则的曲面来拟合各穿刺点处的VTEC值，多项式函数模型解算VTEC的表达式如下：

S-S₀＝(λ-λ₀)+(t-t₀)

其中，n_max、m_max为多项式函数模型的泰勒展开式中的最大维数，E_nm为多项式函数模型的系数，为测区中心点的地理纬度，为穿刺点或者星下点的地理纬度，S₀为测区中心点在该时段中央时刻t₀时的太阳时角，S为穿刺点或者星下点的太阳时角，t为世界时的观测时刻，λ为穿刺点或者星下点的地理经度，λ₀为测区中心点的地理经度；

结合单层模型和载波相位解算的VTEC和多项式函数模型解算的VTEC，可以得到多项式函数模型的观测方程的最终形式：

利用多项式函数模型的观测方程的最终形式，解算电离层延迟概略初值；

(2)电离层延迟修正值解算：

根据多项式函数模型解算的误差，利用B样条函数解算电离层延迟的修正值其模型表达式为：

其中，j₁、j₂为阶数；为模型系数；为二维尺度函数；

(3)电离层延迟的解算：

根据已解算的电离层延迟的初值和利用B样条函数获取的电离层延迟修正值的模型，拟合出更为精确的电离层延迟：

(4)电离层延迟概略预报：

根据步骤(3)获得的电离层延迟，利用Holt-winter(HW)模型进行预报，HW模型的基本公式为：

稳定项： $S_t=\mathrm{\α}\left(\frac{X_t}{I_{t-L}}\right)+(1-\mathrm{\α})(S_{t-1}-b_{t-1})$

趋势项：b_t＝γ(S_t-S_t-1)+(1-γ)b_t-1

季节项： $I_t=\mathrm{\β}\left(\frac{X_t}{S_t}\right)+(1-\mathrm{\β})I_{t-L}$

预报：X_t+τ＝(S_t+τb_t)I_t-L+τ

其中，X_t为t时刻的观测值，S_t为t时刻的稳定项，I_t为t时刻的季节项，b_t为t时刻的趋势项，X_t+τ为τ期的预测值，τ为预测的期数，L为季节长度，α,β,γ∈[0,1]为平滑参数；

(5)电离层延迟修正预报：

利用Markov模型对HW模型的预报值进行更新修正，修正模型表达式为：

${\hat{X}}_t=X_t+{\hat{F}}_t$

其中，为修正后的电离层延迟预报值，X_t为HW模型的预报值，为残差预报值；

残差预报值的计算公式为：

${\hat{F}}_{t+k}=X_t\×E_t\×P^{\left(k\right)}$

其中，X_t为t时刻的电离层延迟，E_t为t时刻的Markov链的数据序列状态，其中以E₁,E₂,…,E_n表示，对应转移时间t₁,t₂,…,t_n，P^(k)为k步状态转移概率矩阵，其计算公式为：

$P^{\left(k\right)}=\left[\begin{array}{ccccc}{P}_{11}^k& P_{12}^k& ...& P_{1(j-1)}^k& P_{1j}^k\\ P_{21}^k& P_{22}^k& ...& P_{2(j-1)}^k& P_{2j}^k\\ P_{31}^k& P_{32}^k& ...& P_{3(j-1)}^k& P_{3j}^k\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ P_{j1}^k& P_{j2}^k& ...& P_{j(j-1)}^k& P_{\mathrm{jj}}^k\end{array}\right]$

其中，表示数列由状态E_i经过k步变为E_j的概率，即：

$P_{\mathrm{ij}}^{\left(k\right)}=n_{\mathrm{ij}}^{\left(k\right)}/N_i$

其中，表示状态E_i经过k步变为E_j的次数，N_i表示状态E_i出现的总次数。