CN110543668B - 一种电离层建模中Kalman滤波状态误差协方差阵的确定方法 - Google Patents

一种电离层建模中Kalman滤波状态误差协方差阵的确定方法 Download PDF

Info

Publication number
CN110543668B
CN110543668B CN201910671565.5A CN201910671565A CN110543668B CN 110543668 B CN110543668 B CN 110543668B CN 201910671565 A CN201910671565 A CN 201910671565A CN 110543668 B CN110543668 B CN 110543668B
Authority
CN
China
Prior art keywords
ionospheric
spherical harmonic
ionosphere
kalman filtering
covariance matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910671565.5A
Other languages
English (en)
Other versions
CN110543668A (zh
Inventor
陈鹏
李政
刘航
马永超
郑乃铨
刘丽霞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Science and Technology
Original Assignee
Xian University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian University of Science and Technology filed Critical Xian University of Science and Technology
Priority to CN201910671565.5A priority Critical patent/CN110543668B/zh
Publication of CN110543668A publication Critical patent/CN110543668A/zh
Application granted granted Critical
Publication of CN110543668B publication Critical patent/CN110543668B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明公开了一种电离层建模中Kalman滤波状态误差协方差阵的确定方法，包括以下步骤：对每一个电离层球谐系数Ci进行历元间差分；计算Kalman滤波历元间隔内电离层球谐系数的变化率；计算电离层球谐系数变化率vCi(k)序列的标准差；得到电离层建模中Kalman滤波状态误差协方差阵Q，Q为对角阵，对角线元素分别为每个电离层球谐系数变化率vCi(k)标准差的三倍的平方。本发明提出的方法有别于现有纯经验确定的方法，该方法根据电离层球谐系数计算Q，可以精确给出电离层球谐系数中每个元素的方差。

Description

一种电离层建模中Kalman滤波状态误差协方差阵的确定方法
技术领域
本发明涉及属于电离层建模领域，涉及一种状态误差协方差阵的确定方法，特别是一种电离层建模中Kalman滤波状态误差协方差阵Q的确定方法。
背景技术
实时电离层建模是当前国际电离层研究领域的热点，不仅在导航定位方面有重大作用，可以显著提高单频实时单点定位用户的定位精度，在科学研究领域同样意义重大，实时获得全球电离层信息，进而实时监测和研究电离层的异常变化，可以对海啸、地震等重大自然灾害作出预报和预警。
目前，全球电离层建模中广泛使用单层电离层假设，即假设电离层电子全部集中在高度为350米～450米的一个无限薄层中。利用如(1)式所示的投影函数MSLM(modifiedsingle-layer mapping)将信号传播路径上的STEC(总电子含量)投影到天顶方向VTEC(垂直总电子含量)。
Figure BDA0002141896000000011
其中，m(z)为电离层投影函数，αm＝0.9782,平均地球半径Re＝6371km，H＝506.7km，角度z是卫星在接收机位置的天顶角，角度z'是电离层穿刺点的天顶角。
球谐函数对球面上数据的拟合具有独特的优势,将得到的穿刺点处的VTEC值，利用合适的模型进行拟合，便可得到全球电离层模型。目前,欧洲定轨中心和武汉大学等国际电离层分析中心均采用球谐函数模型。全球电离层VTEC球谐函数公式如式(2)所示：
Figure BDA0002141896000000021
式中，β为穿刺点的地磁纬度，s＝λ-λ0为穿刺点的日固经度，λ为穿刺点的地磁经度，λ0为太阳的地磁经度，N为球谐函数展开式的最高阶数，
Figure BDA0002141896000000022
为完全正规化的n阶m次勒让德函数，
Figure BDA0002141896000000023
和
Figure BDA0002141896000000024
为未知球谐系数。
Kalman滤波算法是一种基于状态空间模型的线性最小方差估计算法，它考虑了系统的过程噪声和测量噪声的统计特性，在算法上采用递推的形式，能够处理多维和非平稳的随机过程。自从引入Kalman滤波技术以来，就普遍应用于动态参数估计领域，它的主要思想是：当获得新观测值时，状态估计和状态估计的误差被更新。线性Kalman滤波的状态方程和观测方程为：
xk＝Φxk-1+wk-1   (3)
zk＝Hkxk+vk   (4)
式中，xk为k时刻的状态参数向量，在电离层建模中，xk为k时刻的球谐函数模型系数(即式(2)中的
Figure BDA0002141896000000025
和
Figure BDA0002141896000000026
)；Φ为k与k-1时刻的状态转移矩阵，假定状态变量的变化服从随机游走过程，则Φ为单位矩阵；zk为k时刻观测值向量，Hk为设计矩阵，wk-1为状态方程噪声向量，vk为观测噪声向量。wk-1和vk被假定为数学期望为零的白噪声时间序列，两个噪声的协方差假定随着时间的推移是平稳的，并由以下公式给出：
Figure BDA0002141896000000027
Figure BDA0002141896000000028
Kalman滤波的特点是需要精确给出初始信息，如状态向量及其方差协方差阵初始值，状态噪声方差协方差阵Q以及观测噪声方差协方差阵R。
在初始信息全部给定后，即可实现状态向量的预报和修正。在预报阶段，利用状态转移矩阵Φ对k-1时刻的状态向量及方差协方差阵进行预报，获得k时刻状态向量的预报值
Figure BDA0002141896000000031
及方差协方差阵的预报值P′k：
Figure BDA0002141896000000032
Pk'＝ΦPk-1ΦT+Q   (8)
然后，利用k时刻的观测值向量对k时刻的预报值进行修正，得到k时刻状态向量的修正值
Figure BDA0002141896000000033
及其方差协方差阵的修正值Pk：
Figure BDA0002141896000000034
Pk＝(I-KkH)P′k   (10)
式中，Kk称作增益矩阵:
Kk＝P′kHT(HP′kHT+R)-1   (11)
随着时间的不断推进，预报和修正过程反复执行，实现状态向量的更新，直到最后一个时刻。
目前，在使用Kalman滤波求解电离层模型参数时，尚未有精确确定Kalman滤波状态误差协方差阵Q的方法。假设状态向量(球谐系数)随时间的变化服从随机游走过程，则状态转移矩阵为单位阵，而状态方程的不确定性则由Kalman滤波状态误差协方差阵Q反映。Q的精确确定是滤波得到精确结果的关键，如果不能准确给出，则滤波不能得到高精度的结果，甚至可能出现滤波发散。
发明内容
针对现有技术的不足，本发明的目的是提出了一种电离层建模中Kalman滤波状态误差协方差阵的确定方法，根据电离层球谐系数计算电离层建模中Kalman滤波状态误差协方差阵Q的方法，可以精确给出电离层球谐系数中每个元素的方差。
本发明的技术解决方案是：一种电离层建模中Kalman滤波状态误差协方差阵的确定方法，包括以下步骤：
(1)对每一个电离层球谐系数Ci进行历元间差分，即dCi(k)＝Ci(k+1)-Ci(k)，i为电离层球谐系数的序号，dCi(k)表示历元间差分，k表示历元；
(2)假设短时间内电离层球谐系数为线性变化，计算Kalman滤波历元间隔内(如30s或1min)电离层球谐系数的变化率，即
Figure BDA0002141896000000041
vCi(k)表示电离层球谐系数的变化率，Δt表示历元间隔；
(3)计算电离层球谐系数变化率vCi(k)序列的标准差；
(4)得到电离层建模中Kalman滤波状态误差协方差阵Q，Q为对角阵，对角线元素分别为每个电离层球谐系数变化率vCi(k)标准差的3倍的平方。
附图说明
图1是2018年第一个电离层球谐系数C00的时间序列。
图2是电离层球谐系数C00在Kalman滤波历元间隔(1分钟)的变化率。
图3是2018年一年内电离层球谐系数C00每分钟变化率的直方图。
图4是利用2018年全年全部的电离层球谐系数计算出的全球电离层建模中Kalman滤波状态误差协方差阵的对角线元素序列图。
图5是本发明的方法示意图。
具体实施方式
下面结合附图，以2018年1小时间隔的电离层球谐系数为例，说明电离层建模中Kalman滤波状态误差协方差阵的计算过程，并予以详细描述。
如图5所示，一种电离层建模中Kalman滤波状态误差协方差阵的确定方法，包括以下步骤：
(1)对每一个电离层球谐系数Ci进行历元间差分，即dCi(k)＝Ci(k+1)-Ci(k)，i为电离层球谐系数的序号(i＝1,2,...,256)，dCi(k)表示历元间差分，k表示历元。如图1所示，给出了2018年第一个电离层球谐系数C00的时间序列,可以看出C00在一年内呈现出明显的周期性变化。
(2)假设短时间内电离层球谐系数为线性变化，计算Kalman滤波历元间隔内(如30s或1min)电离层球谐系数的变化率，即
Figure BDA0002141896000000051
vCi(k)表示电离层球谐系数的变化率，Δt表示历元间隔。如图2所示，对电离层球谐系数C00求差，并假设电离层球谐系数C00在短时间内为线性变化，得到电离层球谐系数C00在Kalman滤波历元间隔(1分钟)的变化率，系数的变化率无明显的趋势性变化，以随机变化为主。
(3)计算电离层球谐系数变化率vCi(k)序列的标准差。如图3所示，给出了2018年一年内电离层球谐系数C00每分钟变化率的直方图，变化率服从均值为零，标准差为2.7×10-3的正态分布，则电离层建模中Kalman滤波状态误差协方差阵对角线上第一个元素即为2.7×10-3×3×2.7×10-3×3＝6.561×10-5。
(4)电离层建模中Kalman滤波状态误差协方差阵中的其他对角线元素均按照上述方法计算，得到全球电离层建模中Kalman滤波状态误差协方差阵Q，Q为对角阵，对角线元素分别为每个电离层球谐系数变化率vCi(k)标准差的3倍的平方。如图4所示，利用2018年全年全部的电离层球谐系数计算出的Q矩阵对角线元素序列图。
本发明提出了一种电离层建模中Kalman滤波状态误差协方差阵的确定方法，可以精确给出状态向量中每个元素的方差。在使用Kalman滤波求解电离层模型参数时，通过本发明可以得到精确确定的Kalman滤波状态误差协方差阵Q，避免可能出现的滤波发散情况，得到高精度的滤波结果，从而得到高精度的全球电离层模型。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种电离层建模中Kalman滤波状态误差协方差阵的确定方法，其特征在于，根据电离层球谐系数计算Kalman滤波状态误差协方差阵Q，精确给出电离层球谐系数中每个元素的方差；
包括以下步骤：
(1)对每一个电离层球谐系数Ci进行历元间差分，即dCi(k)＝Ci(k+1)-Ci(k)，i为电离层球谐系数的序号(i＝1,2,...,256)，dCi(k)表示历元间差分，k表示历元；
(2)假设短时间内电离层球谐系数为线性变化，计算Kalman滤波历元间隔内电离层球谐系数的变化率，即
Figure FDA0002141895990000011
vCi(k)表示电离层球谐系数的变化率，Δt表示历元间隔；
(3)计算电离层球谐系数变化率vCi(k)序列的标准差；
(4)得到电离层建模中Kalman滤波状态误差协方差阵Q，Q为对角阵，对角线元素分别为每个电离层球谐系数变化率vCi(k)标准差的三倍的平方。
CN201910671565.5A 2019-07-24 2019-07-24 一种电离层建模中Kalman滤波状态误差协方差阵的确定方法 Active CN110543668B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910671565.5A CN110543668B (zh) 2019-07-24 2019-07-24 一种电离层建模中Kalman滤波状态误差协方差阵的确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910671565.5A CN110543668B (zh) 2019-07-24 2019-07-24 一种电离层建模中Kalman滤波状态误差协方差阵的确定方法

Publications (2)

Publication Number Publication Date
CN110543668A CN110543668A (zh) 2019-12-06
CN110543668B true CN110543668B (zh) 2023-04-07

Family

ID=68709831

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910671565.5A Active CN110543668B (zh) 2019-07-24 2019-07-24 一种电离层建模中Kalman滤波状态误差协方差阵的确定方法

Country Status (1)

Country Link
CN (1) CN110543668B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013061768A (ja) * 2011-09-13 2013-04-04 Nippon Telegr & Teleph Corp <Ntt> 最適モデル推定装置、方法、及びプログラム
CN109828288A (zh) * 2019-01-23 2019-05-31 东南大学 一种基于区域cors的实时电离层建模与监测方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2013061768A (ja) * 2011-09-13 2013-04-04 Nippon Telegr & Teleph Corp <Ntt> 最適モデル推定装置、方法、及びプログラム
CN109828288A (zh) * 2019-01-23 2019-05-31 东南大学 一种基于区域cors的实时电离层建模与监测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
几种地基GPS区域电离层TEC建模方法的比较及其一致性研究;柳景斌等;《武汉大学学报(信息科学版)》(第05期);全文 *
球谐函数电离层模型研究进展;刘继州;《洛阳师范学院学报》(第02期);全文 *

Also Published As

Publication number Publication date
CN110543668A (zh) 2019-12-06

Similar Documents

Publication Publication Date Title
US11237276B2 (en) System and method for gaussian process enhanced GNSS corrections generation
US20240036216A1 (en) Systems and methods for high-integrity satellite positioning
US6505122B1 (en) Method and apparatus for providing accurate position estimates in instances of severe dilution of precision
CN110568459B (zh) 基于igs和cors站的区域电离层tec实时监测方法
US20240142637A1 (en) System and method for gaussian process enhanced gnss corrections generation
CN110455287A (zh) 自适应无迹卡尔曼粒子滤波方法
JP2011221035A (ja) ローカルｒｔｋシステムと、地域、広域、またはグローバル搬送波位相測位システムを組み合わせて利用する方法
JP2007529010A (ja) ２周波数の一方で測定データが利用できない場合に短期間バックアップ２周波数ナビゲーションを行なう方法
EP2330382B1 (en) Method and system for latitude adaptive navigation quality estimation
CN114966760B (zh) 一种电离层加权的非差非组合ppp-rtk技术实现方法
EP3430437A1 (en) Navigation satellite wide-lane bias determination system and method
CN115373007B (zh) 基于手机gnss模糊度相对变化估计的里程计定位方法
CN109597105A (zh) 一种顾及载波系统间偏差的gps/glonass紧组合定位方法
CN113504557B (zh) 面向实时应用的gps频间钟差新预报方法
CN111123345A (zh) 一种基于gnss测量的经验电离层模型数据驱动方法
CN112799110A (zh) 一种顾及Doppler的北斗修正伪距单点定位方法、系统与设备
CN110543668B (zh) 一种电离层建模中Kalman滤波状态误差协方差阵的确定方法
Chen et al. Near real-time global ionospheric modeling based on an adaptive Kalman filter state error covariance matrix determination method
CN110716219A (zh) 一种提高定位解算精度的方法
CN113671551B (zh) Rtk定位解算方法
CN112987051B (zh) 一种提高卫星导航定位性能的方法
CN115980803B (zh) 基于双频码伪距和载波相位观测量进行伪距平滑的方法
Shenglong et al. A multidimensional gross error positioning method for GNSS navigation satellite based on UKF algorithm
CN115061166B (zh) 一种载波相位重构方法、装置、电子设备及介质
CN112444187B (zh) 形变监测方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant