CN116907503A - 一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统 - Google Patents

一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统 Download PDF

Info

Publication number
CN116907503A
CN116907503A CN202310840461.9A CN202310840461A CN116907503A CN 116907503 A CN116907503 A CN 116907503A CN 202310840461 A CN202310840461 A CN 202310840461A CN 116907503 A CN116907503 A CN 116907503A
Authority
CN
China
Prior art keywords
distribution
remote sensing
noise
state
measurement
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.)
Pending
Application number
CN202310840461.9A
Other languages
English (en)
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.)
Harbin Institute of Technology Weihai
Original Assignee
Harbin Institute of Technology Weihai
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 Harbin Institute of Technology Weihai filed Critical Harbin Institute of Technology Weihai
Priority to CN202310840461.9A priority Critical patent/CN116907503A/zh
Publication of CN116907503A publication Critical patent/CN116907503A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20—Instruments for performing navigational calculations
    • 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
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统，涉及遥感卫星编队定位技术。本发明要解决的技术问题为：在面对各种无法精确确定的扰动力、强扰动、复杂的宇宙环境、导航传感器不足等问题时，导航系统的系统噪声与量测噪声出现野值的问题，进而导致编队卫星之间的定位精度受到严重影响。技术要点：建立一个具有随机状态空间模型的非线性系统；进行时间更新解算出状态误差协方差矩阵；进行量测更新中的迭代初始化；进行量测更新，利用加权求和的方法解算出量测预测矩阵；参数迭代更新，将厚尾噪声分层高斯模型的联合后验分布通过变分贝叶斯方法找到联合后验分布的联合近似分布；计算状态误差和量测噪声协方差修正值，由此减小非高斯厚尾噪声对系统定位精度的影响；更新状态变量Xk+1的后验近似分布，从而完成遥感编队卫星定位。

Description

一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及 系统
技术领域
本发明涉及遥感卫星编队定位技术，具体一种基于抗野值鲁棒定位算法的遥感卫星编队定位方法。
背景技术
针对近地环境下运行的遥感卫星编队，在面对各种无法精确确定的扰动力、强扰动、复杂的宇宙环境、导航传感器不足等问题时，导航系统的系统噪声与量测噪声出现野值的问题，进而编队卫星之间的定位精度受到严重影响。在这种情况下如何提高编队卫星之间的相对定位精度和鲁棒性，改进抗野值鲁棒定位现有技术中没有人提及和解决此技术问题。
文献号为CN111017264B现有技术，公开了一种高效立体遥感卫星编队方法，将立体成像所需的多个空间相机分别安装在多颗卫星上，每颗卫星可以彼此独立实施姿态控制，通过姿态控制的相互配合，使多颗卫星从目标的不同方向实施观察。每颗卫星都可以进行长时间大条带扫描成像，在立体成像过程中无需进行大角度姿态机动。但是针对遥感编队卫星定位没有提及。
发明内容
本发明要解决的技术问题是：
本发明的目的是提供一种基于抗野值鲁棒定位算法的遥感卫星编队定位方法及系统，以解决如下技术问题：在面对各种无法精确确定的扰动力、强扰动、复杂的宇宙环境、导航传感器不足等问题时，导航系统的系统噪声与量测噪声出现野值的问题，进而导致编队卫星之间的定位精度受到严重影响。
本发明为解决上述技术问题所采用的技术方案为：
一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，所述的方法的实现过程为：
步骤一：建立一个具有随机状态空间模型的非线性系统，所述非线性系统体现系统函数、量测函数、系统状态变量以及量测变量，用学生t分布描述非线性系统噪声具有非高斯厚尾现象；
步骤二：进行时间更新，解算出状态预测矩阵和状态误差协方差矩阵；
步骤三：进行量测更新中的迭代初始化；
步骤四：进行量测更新，利用加权求和的方法解算出量测预测矩阵；
步骤五：参数迭代更新，将厚尾噪声分层高斯模型的联合后验分布通过变分贝叶斯方法求解出联合后验分布的联合近似分布；
步骤六：计算状态误差和量测噪声协方差修正值，由此减小非高斯厚尾噪声对系统定位精度的影响；
步骤七：利用修正后的状态量进行可变参数更新，可变参数包括卡尔曼滤波增益、状态变量、误差协方差矩阵，并给出自相关协方差矩阵和互相关协方差矩阵，根据更新的可变参数从而完成遥感编队卫星定位。
进一步地，步骤一的具体实现过程为：
具有随机状态空间模型的非线性系统表达式为：
其中k+1表示系统所处时间，fk(·)和hk()分别表示系统的系统函数和量测函数，Xk+1表示系统状态在k+1时刻的值，系统的状态变量为 为X轴相对速度，x为X轴相对位置，/>为Y轴相对速度，y为Y轴相对位置，/>为Z轴相对速度，z为Z轴相对位置；Zk+1表示系统量测在k+1时刻的值，系统的量测变量为[ραβ]，其中，ρ表示参考卫星与跟随卫星的相对位置，α代表遥感卫星编队中无线电所测得的参考卫星到跟随卫星的方位角，β代表遥感卫星编队中无线电所测得的参考卫星到跟随卫星的俯仰角；wk表示关系k时刻的系统噪声，vk+1表示k+1时刻的量测噪声；
通过学生t分布来描述非线性系统中噪声具有非高斯厚尾现象，将系统的量测噪声的概率分布写成学生t分布的形式：
p(vk+1)＝St(vk+1；0,Rk+1,υk+1) (2)
式(2)中，学生t分布的均值等于0，Rk+1和υk+1分别指的是其尺度矩阵与自由度参数；为获取非线性系统的后验分布概率近似解，可将从前一时刻到下一时刻的一步预测概率p(Xk+1|Z1:k)与似然概率p(Zk+1|Xk+1)可分别表示成学生t分布的形式：
p(Zk+1|Xk+1)＝St(Zk+1；h(Xk+1),Rk+1,υk+1) (4)
式(3)和(4)中，μk+1，Σk+1和δ分别指学生t分布下的均值、尺度矩阵与自由度参数；
因为学生t分布相当于是无穷个高斯分布的混合，则可将式(3)与式(4)改写成：
式(5)和(6)中，ζk+1与λk+1分别指的是p(Xk+1|Z1:k)和p(Zk+1|Xk+1)的辅助随机参数，G(a,b,c)指的是满足形状参数为b且尺度参数为c的Gamma分布；
则p(Xk+1|Z1:k)与p(Zk+1|Xk+1)的分层高斯形式可表示成：
通过上述过程可完成学生t分布下的厚尾噪声分层高斯模型的建立，将具有非高斯厚尾特性的随机状态空间模型下的状态估计问题转变成基于学生t分布的厚尾噪声分层高斯模型下的状态估计问题。
进一步地，步骤二中进行时间更新的具体实现过程为：
初始化：
对协方差矩阵Pk|k进行Cholesky分解，分解为上三角矩阵Uk|k的乘积：
通过因式分解求得容积点：
其中ζj为第j个容积点，其表达式为：
其中，n为状态变量的维度，[1]为单位矩阵；
利用系统函数对解算得的容积点进行传播：
利用加权求和的方法解算出状态预测值：
解算出状态误差协方差矩阵：
进一步地，步骤三中进行量测更新中的迭代初始化的具体为：
初始状态预测值：
初始状态协方差矩阵：
进一步地，步骤四中进行量测更新的具体实现过程为：
对时间更新中所解算的状态预测值进行Cholesky分解，分解为上三角矩阵Uk+1|k的乘积：
利用因式分解的结果解算新生成的容积点：
使用量测函数对解算得的容积点进行传播：
Zj,k+1|k＝h(Xj,k+1|k) (18)
利用加权求和的方法解算出k+1时刻的量测预测矩阵：
进一步地，步骤五中参数迭代更新的具体实现过程为：
设定所求厚尾噪声分层高斯模型的联合后验分布为P(θk+1|Z1:k+1)，通过变分贝叶斯方法可以找到联合后验分布的联合近似分布：
P(θk+1|Z1:k+1)≈q(Xk+1)q(ζk+1)q(λk+1)q(Σk+1) (20)
式(20)中，θk+1指的是由所需参数Xk+1、ζk+1、λk+1、Σk+1组成的集合，q(·)指的是各变分参数的近似后验分布；通过KL散度最小化的方法，可以分别得到变分参数的近似后验分布；利用定点迭代方法，在更新一个所需参数的同时其他变分参数维持不变，以此求得各变分参数的局部最优解；
将q(i+1)(ζk+1)以Gamma分布的形式进行更新：
其中q(i+1)(·)指的是i+1次迭代时所求变分参数的近似后验分布，其形状参数为尺度参数为/>Σk+1和/>分别指学生t分布下的尺度矩阵与自由度参数；E(i)[·]指的是i次迭代所得到的其余变分参数期望，tr(·)指的是求迹计算，/>是一个参数，可由下式求出:
将q(i+1)(λk+1)以Gamma分布的形式进行更新，则：
其中表示i+1次迭代时所求变分参数的近似后验分布的形状参数，其尺度参数为/>Rk+1和/>分别指学生t分布下的尺度矩阵与自由度参数，参数/>可以通过下式得到：
将q(i+1)(Σk+1)以Inv-Wishart分布的形式进行更新，则：
其中表示Σk+1的自由度参数，/>为Σk+1的逆尺度矩阵；
将q(i+1)(Xk+1)近似成高斯分布的形式：
进一步地，步骤六中计算状态误差和量测噪声协方差修正值的具体实现过程为：
状态误差协方差修正值：
量测噪声协方差修正值：
式(27)与(28)中，各期望值可由下式计算：
式(29)中，E(i+1)(·)指的是各个变分参数的期望值，ψ(·)指的是Gamma函数的对数的导数。
进一步地，步骤七中更新状态变量Xk+1的后验近似分布的具体实现过程为：
卡尔曼滤波增益：
状态变量：
误差协方差矩阵：
其中：
一种基于抗野值鲁棒定位算法的遥感编队卫星定位系统，该系统具有与上述技术方案的步骤对应的程序模块，运行时执行上述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法中的步骤。
一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序配置为由处理器调用时实现上述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法的步骤。
本发明具有以下有益技术效果：
本发明以变分贝叶斯定位算法为基础，在学生t分布的基础上将后验分布近似为高斯分布，以改进抗野值鲁棒定位算法，提高提高编队卫星之间的相对定位精度和鲁棒性。
经实验证明：当系统噪声和量测噪声均存在野值时，本发明所提出的SIVBCKF定位算法的定位误差为0.841m，与HMCKF定位算法的0.9581m相比，具有更高的定位精度。其相对定位精度比HMCKF定位算法提升了12.2％，具有较好的鲁棒性，这是因为SIVBCKF定位算法在学生t分布的基础上将后验分布近似为高斯分布，然后通过变分贝叶斯的方法对定位误差进行估计，具有较高的定位精度。从二者的平均运行时间来看，SIVBCKF定位算法的运行时间略长于HMCKF定位算法，这是由于SIVBCKF定位算法中有固定点迭代过程，但其获得了更高的定位精度。
附图说明
图1为本发明所述方法的流程框图(所提算法的流程图)；图2为遥感卫星空间圆编队仿真图；图3为跟随卫星F1的三轴定位误差图；图4为跟随卫星F1的相对定位误差图；图5为跟随卫星F1的三轴速度误差图。
具体实施方式
结合图1至图5，针对本发明所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法的实现是进行如下阐述：
步骤一：建立一个具有随机状态空间模型的非线性系统：
其中k+1表示系统所处时间，fk(·)和hk(·)分别表示系统的系统函数和量测函数，Xk+1表示系统状态在k+1时刻的值，系统的状态变量为 为X轴相对速度，x为X轴相对位置，/>为Y轴相对速度，y为Y轴相对位置，/>为Z轴相对速度，z为Z轴相对位置。Zk+1表示系统量测在k+1时刻的值，系统的量测变量为[ρ α β]，其中，ρ表示参考卫星与跟随卫星的相对位置，α代表遥感卫星编队中无线电所测得的参考卫星到跟随卫星的方位角，β代表遥感卫星编队中无线电所测得的参考卫星到跟随卫星的俯仰角。wk表示关系k时刻的系统噪声，vk+1表示k+1时刻的量测噪声。
在理论上，通过学生t分布来描述非线性系统中噪声具有非高斯厚尾现象的问题较为适合，因此将系统的量测噪声的概率分布写成学生t分布的形式：
p(vk+1)＝St(vk+1；0,Rk+1,υk+1) (2)
式(2)中，学生t分布的均值等于0，Rk+1和υk+1分别指的是其尺度矩阵与自由度参数。为获取非线性系统的后验分布概率近似解，可将一步预测概率p(xk+1|Z1:k)与似然概率p(Zk+1|Xk+1)可分别表示成学生t分布的形式：
式(3)和(4)中，μk+1，Σk+1和δ分别指学生t分布下的均值、尺度矩阵与自由度参数。
因为学生t分布相当于是无穷个高斯分布的混合，则可将式(3)与式(4)改写成：
式(5)和(6)中，ζk+1与λk+1分别指的是p(Xk+1|Z1:k)和p(Zk+1|Xk+1)的辅助随机参数，G(a,b,c)指的是满足形状参数为b且尺度参数为c的Gamma分布。
则p(Xk+1|Z1:k)与p(Zk+1|Xk+1)的分层高斯形式可表示成：
通过上述过程可完成学生t分布下的厚尾噪声分层高斯模型的建立。此时将具有非高斯厚尾特性的随机状态空间模型下的状态估计问题转变成基于学生t分布的厚尾噪声分层高斯模型下的状态估计问题。
步骤二：进行时间更新：
对协方差矩阵进行Cholesky分解，分解为上三角矩阵Uk|k的乘积：
通过因式分解求得容积点：
其中ζj为第j个容积点，其表达式为：
其中，n为状态变量的维度，[1]为单位矩阵。
利用系统函数对解算得的容积点进行传播：
利用加权求和的方法解算出状态预测值：
解算出状态误差协方差矩阵：
步骤三：进行量测更新中的迭代初始化：
初始状态预测值：
初始状态协方差矩阵：
步骤四：进行量测更新：
对时间更新中所解算的状态预测值进行Cholesky分解，分解为上三角矩阵Uk+1|k的乘积：
利用因式分解的结果解算新生成的容积点：
使用量测函数对解算得的容积点进行传播：
Zj,k+1|k＝h(Xj,k+1|k) (18)
利用加权求和的方法解算出k+1时刻的量测预测矩阵：
步骤五：参数迭代更新
设定所求厚尾噪声分层高斯模型的联合后验分布为P(θk+1|Z1:k+1)，通过变分贝叶斯方法可以找到联合后验分布的联合近似分布：
P(θk+1|Z1:k+1)≈q(Xk+1)q(ζk+1)q(λk+1)q(Σk+1) (20)
式(20)中，θk+1指的是由所需参数Xk+1、ζk+1、λk+1、Σk+1组成的集合，q(·)指的是各变分参数的近似后验分布。通过KL散度最小化的方法，可以分别得到变分参数的近似后验分布。由于变分参数之间存在耦合问题，所以利用定点迭代方法，在更新一个所需参数的同时其他变分参数维持不变，以此求得各变分参数的局部最优解。
将q(i+1)(ζk+1)以Gamma分布的形式进行更新：
其中q(i+1)(·)指的是i+1次迭代时所求变分参数的近似后验分布，其形状参数为尺度参数为/>Σk+1和/>分别指学生t分布下的尺度矩阵与自由度参数。E(i)[·]指的是i次迭代所得到的其余变分参数期望，tr(·)指的是求迹计算，/>是一个参数，可由下式求出:
将q(i+1)(λk+1)以Gamma分布的形式进行更新，则：
其中表示i+1次迭代时所求变分参数的近似后验分布的形状参数，其尺度参数为/>Rk+1和/>分别指学生t分布下的尺度矩阵与自由度参数，参数/>可以通过下式得到：
将q(i+1)(Σk+1)以Inv-Wishart分布的形式进行更新，则：
其中表示Σk+1的自由度参数，/>为Σk+1的逆尺度矩阵。
将q(i+1)(Xk+1)近似成高斯分布的形式：
步骤六：计算状态误差和量测噪声协方差修正值：
状态误差协方差修正值：
量测噪声协方差修正值：
式(27)与(28)中，各期望值可由下式计算：
式(29)中，E(i+1)(g)指的是各个变分参数的期望值，ψ(g)指的是Gamma函数的对数的导数。
步骤七：更新状态变量Xk+1的后验近似分布：
卡尔曼滤波增益：
状态变量：
误差协方差矩阵：
其中：
本发明的仿真分析：
(1)遥感卫星编队设计：
设置遥感卫星编队中的参考卫星的轨道参数如下：轨道高度为520km，偏心率为0°，轨道倾角为97°，升交点赤经为0°，近地点幅角为0°，真近点角为0°，周期为5701.76s，轨道类型为圆形太阳同步轨道。其余跟随卫星以参考卫星为圆心，初始相位分别为45°、135°、225°、315°，构成一个参考卫星与跟随卫星星间距离为100m的“一主四从”空间圆编队。跟随卫星的具体轨道参数如表1所示：
表1遥感卫星编队中各跟随卫星轨道参数
把表1的参数代入到STK轨道参数模块中，得到所设计的遥感卫星编队仿真图如图2所示：
(2)算法仿真：
由于四颗跟随卫星相对于参考卫星的运动轨迹都是半径为100m的空间圆，具有一致的运动特性，所以只以跟随卫星F1为例，设置系统状态更新频率为1s，系统运行时间为6000s，相对测距精度为1m，方位角与俯仰角的相对测角精度均为0.000175rad。在遥感卫星编队中，跟随卫星F1的系统真实状态初值X0、初始系统噪声协方差Q0与初始量测噪声协方差R0可分别表示为：
X0＝[35.3943 70.6327 61.30470.0389-0.07800.0674]T (36)
Q0＝diag[10-10 10-10 10-10 10-9 10-9 10-9] (37)
R0＝diag[1 1 1] (38)
设置产生野值噪声的条件：
现对改进变分贝叶斯抗野值定位算法(SIVBCKF)与现有的Huber-M抗野值定位算法(HMCKF)进行性能对比：
图3展示了使用HMCKF和SIVBCKF两种算法进行相对定位导航后，X、Y、Z三个轴的定位误差。如图3所示，当系统噪声与量测噪声呈非高斯厚尾分布时，相较于HMCKF定位算法，SIVBCKF定位算法具有更好的跟踪效果。具体而言，SIVBCKF定位算法在X轴上的最大误差约为6m，Y轴上的最大误差约为2m，Z轴上的最大误差约为3m.HMCKF定位算法的X轴最大误差约7.6m，Y轴最大误差约6.6m，Z轴最大误差约为3.8m；SIVBCKF定位算法具有较好的鲁棒性。
跟随卫星F1分别应用HMCKF定位算法与SIVBCKF定位算法进行相对定位导航的相对定位误差图如图4所示。从图中可以看出SIVBCKF定位算法的定位误差小于HMCKF定位算法的定位误差，SIVBCKF定位算法的曲线数值波动较小，SIVBCKF定位算法的最大定位误差在3.3m之内且最终定位精度大致为0.3m，HMCKF定位算法的曲线数值波动较大，其最大定位误差在3.8m之内且最终定位精度大致为0.4m，说明在该场景下，SIVBCKF的定位精度较高。
跟随卫星F1分别应用HMCKF定位算法与SIVBCKF定位算法进行相对定位导航的X，Y，Z三个方向的速度误差图如图5所示。从图中可知，在系统运行结束时，HMCKF定位算法与SIVBCKF定位算法的速度误差都是收敛的，且SIVBCKF定位算法的曲线整体在HMCKF定位算法之下，各轴的速度误差突变幅度较小，说明其具有较好的鲁棒性能。
下面对该场景下的两种定位算法的定位性能进行定量分析，在系统运行时间内，具体均方根误差数据和系统平均运行时间如表2所示：
表2系统噪声与量测噪声非高斯厚尾时两种定位算法性能比较
根据表2，当系统噪声和量测噪声均存在野值时，本发明所提出的SIVBCKF定位算法的定位误差为0.841m，与HMCKF定位算法的0.9581m相比，具有更高的定位精度。其相对定位精度比HMCKF定位算法提升了12.2％，具有较好的鲁棒性，这是因为SIVBCKF定位算法在学生t分布的基础上将后验分布近似为高斯分布，然后通过变分贝叶斯的方法对定位误差进行估计，具有较高的定位精度。从二者的平均运行时间来看，SIVBCKF定位算法的运行时间略长于HMCKF定位算法，这是由于SIVBCKF定位算法中有固定点迭代过程，但其获得了更高的定位精度。

Claims (10)

1.一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，所述的方法的实现过程为：
步骤一：建立一个具有随机状态空间模型的非线性系统，所述非线性系统体现系统函数、量测函数、系统状态变量以及量测变量，用学生t分布描述非线性系统噪声具有非高斯厚尾现象；
步骤二：进行时间更新，解算出状态预测矩阵和状态误差协方差矩阵；
步骤三：进行量测更新中的迭代初始化；
步骤四：进行量测更新，利用加权求和的方法解算出量测预测矩阵；
步骤五：参数迭代更新，将厚尾噪声分层高斯模型的联合后验分布通过变分贝叶斯方法求解出联合后验分布的联合近似分布；
步骤六：计算状态误差和量测噪声协方差修正值，由此减小非高斯厚尾噪声对系统定位精度的影响；
步骤七：利用修正后的状态量进行可变参数更新，可变参数包括卡尔曼滤波增益、状态变量、误差协方差矩阵，并给出自相关协方差矩阵和互相关协方差矩阵，根据更新的可变参数从而完成遥感编队卫星定位。
2.根据权利要求1所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤一的具体实现过程为：
具有随机状态空间模型的非线性系统表达式为：
其中k+1表示系统所处时间，fk(·)和hk(·)分别表示系统的系统函数和量测函数，Xk+1表示系统状态在k+1时刻的值，系统的状态变量为 为X轴相对速度，x为X轴相对位置，/>为Y轴相对速度，y为Y轴相对位置，/>为Z轴相对速度，z为Z轴相对位置；Zk+1表示系统量测在k+1时刻的值，系统的量测变量为[ρ α β]，其中，ρ表示参考卫星与跟随卫星的相对位置，α代表遥感卫星编队中无线电所测得的参考卫星到跟随卫星的方位角，β代表遥感卫星编队中无线电所测得的参考卫星到跟随卫星的俯仰角；wk表示关系k时刻的系统噪声，vk+1表示k+1时刻的量测噪声；
通过学生t分布来描述非线性系统中噪声具有非高斯厚尾现象，将系统的量测噪声的概率分布写成学生t分布的形式：
p(vk+1)＝St(vk+1；0,Rk+1,υk+1) (2)
式(2)中，学生t分布的均值等于0，Rk+1和υk+1分别指的是其尺度矩阵与自由度参数；为获取非线性系统的后验分布概率近似解，可将从前一时刻到下一时刻的一步预测概率p(Xk+1|Z1:k)与似然概率p(Zk+1|Xk+1)可分别表示成学生t分布的形式：
p(Zk+1|Xk+1)＝St(Zk+1；h(Xk+1),Rk+1,υk+1) (4)
式(3)和(4)中，μk+1，Σk+1和δ分别指学生t分布下的均值、尺度矩阵与自由度参数；
因为学生t分布相当于是无穷个高斯分布的混合，则可将式(3)与式(4)改写成：
式(5)和(6)中，ζk+1与λk+1分别指的是p(Xk+1|Z1:k)和p(Zk+1|Xk+1)的辅助随机参数，G(a,b,c)指的是满足形状参数为b且尺度参数为c的Gamma分布；
则p(Xk+1|Z1:k)与p(Zk+1|Xk+1)的分层高斯形式可表示成：
通过上述过程可完成学生t分布下的厚尾噪声分层高斯模型的建立，将具有非高斯厚尾特性的随机状态空间模型下的状态估计问题转变成基于学生t分布的厚尾噪声分层高斯模型下的状态估计问题。
3.根据权利要求1或2所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤二中进行时间更新的具体实现过程为：
初始化：
对协方差矩阵Pk|k进行Cholesky分解，分解为上三角矩阵Uk|k的乘积：
通过因式分解求得容积点：
其中ζj为第j个容积点，其表达式为：
其中，n为状态变量的维度，[1]为单位矩阵；
利用系统函数对解算得的容积点进行传播：
利用加权求和的方法解算出状态预测值：
解算出状态误差协方差矩阵：
4.根据权利要求3所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤三中进行量测更新中的迭代初始化的具体为：
初始状态预测值：
初始状态协方差矩阵：
5.根据权利要求4所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤四中进行量测更新的具体实现过程为：
对时间更新中所解算的状态预测值进行Cholesky分解，分解为上三角矩阵Uk+1|k的乘积：
利用因式分解的结果解算新生成的容积点：
使用量测函数对解算得的容积点进行传播：
Zj,k+1|k＝h(Xj,k+1|k) (18)
利用加权求和的方法解算出k+1时刻的量测预测矩阵：
6.根据权利要求5所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤五中参数迭代更新的具体实现过程为：
设定所求厚尾噪声分层高斯模型的联合后验分布为P(θk+1|Z1:k+1)，通过变分贝叶斯方法可以找到联合后验分布的联合近似分布：
P(θk+1|Z1:k+1)≈q(Xk+1)q(ζk+1)q(λk+1)q(Σk+1) (20)
式(20)中，θk+1指的是由所需参数Xk+1、ζk+1、λk+1、Σk+1组成的集合，q()指的是各变分参数的近似后验分布；通过KL散度最小化的方法，可以分别得到变分参数的近似后验分布；利用定点迭代方法，在更新一个所需参数的同时其他变分参数维持不变，以此求得各变分参数的局部最优解；
将q(i+1)(ζk+1)以Gamma分布的形式进行更新：
其中q(i+1)(·)指的是i+1次迭代时所求变分参数的近似后验分布，其形状参数为尺度参数为/>Σk+1和/>分别指学生t分布下的尺度矩阵与自由度参数；E(i)[·]指的是i次迭代所得到的其余变分参数期望，tr(·)指的是求迹计算，/>是一个参数，可由下式求出:
将q(i+1)(λk+1)以Gamma分布的形式进行更新，则：
其中表示i+1次迭代时所求变分参数的近似后验分布的形状参数，其尺度参数为Rk+1和/>分别指学生t分布下的尺度矩阵与自由度参数，参数/>可以通过下式得到：
将q(i+1)(Σk+1)以Inv-Wishart分布的形式进行更新，则：
其中表示Σk+1的自由度参数，/>为Σk+1的逆尺度矩阵；
将q(i+1)(Xk+1)近似成高斯分布的形式：
7.根据权利要求6所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤六中计算状态误差和量测噪声协方差修正值的具体实现过程为：
状态误差协方差修正值：
量测噪声协方差修正值：
式(27)与(28)中，各期望值可由下式计算：
式(29)中，E(i+1)(·)指的是各个变分参数的期望值，ψ(·)指的是Gamma函数的对数的导数。
8.根据权利要求7所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法，其特征在于，步骤七中更新状态变量Xk+1的后验近似分布的具体实现过程为：
卡尔曼滤波增益：
状态变量：
误差协方差矩阵：
其中：
9.一种基于抗野值鲁棒定位算法的遥感编队卫星定位系统，其特征在于：该系统具有与上述权利要求1-8任一项权利要求的步骤对应的程序模块，运行时执行上述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法中的步骤。
10.一种计算机可读存储介质，其特征在于：所述计算机可读存储介质存储有计算机程序，所述计算机程序配置为由处理器调用时实现权利要求1-8中任一项所述的一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法的步骤。
CN202310840461.9A 2023-07-10 2023-07-10 一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统 Pending CN116907503A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310840461.9A CN116907503A (zh) 2023-07-10 2023-07-10 一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310840461.9A CN116907503A (zh) 2023-07-10 2023-07-10 一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统

Publications (1)

Publication Number Publication Date
CN116907503A true CN116907503A (zh) 2023-10-20

Family

ID=88352421

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310840461.9A Pending CN116907503A (zh) 2023-07-10 2023-07-10 一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统

Country Status (1)

Country Link
CN (1) CN116907503A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117367431A (zh) * 2023-10-23 2024-01-09 哈尔滨工业大学（威海） 带未知量测偏置的mems与uwb紧组合定位方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117367431A (zh) * 2023-10-23 2024-01-09 哈尔滨工业大学（威海） 带未知量测偏置的mems与uwb紧组合定位方法及系统
CN117367431B (zh) * 2023-10-23 2024-05-14 哈尔滨工业大学（威海） 带未知量测偏置的mems与uwb紧组合定位方法及系统

Similar Documents

Publication Publication Date Title
CN109211276B (zh) 基于gpr与改进的srckf的sins初始对准方法
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN111156987B (zh) 基于残差补偿多速率ckf的惯性/天文组合导航方法
CN109724599B (zh) 一种抗野值的鲁棒卡尔曼滤波sins/dvl组合导航方法
CN111238467B (zh) 一种仿生偏振光辅助的无人作战飞行器自主导航方法
CN111783307A (zh) 一种高超声速飞行器状态估计方法
CN110567455B (zh) 一种求积更新容积卡尔曼滤波的紧组合导航方法
CN114777812B (zh) 一种水下组合导航系统行进间对准与姿态估计方法
Yoo et al. Improvement of terrain referenced navigation using a point mass filter with grid adaptation
CN116907503A (zh) 一种基于抗野值鲁棒定位算法的遥感编队卫星定位方法及系统
CN113916226A (zh) 一种基于最小方差的组合导航系统抗扰滤波方法
CN110912535B (zh) 一种新型无先导卡尔曼滤波方法
CN111426322B (zh) 同时估计状态和输入的自适应目标跟踪滤波方法及系统
CN110567490B (zh) 一种大失准角下sins初始对准方法
CN111982126A (zh) 一种全源BeiDou/SINS弹性状态观测器模型设计方法
CN113503891B (zh) 一种sinsdvl对准校正方法、系统、介质及设备
CN114296069B (zh) 一种基于视觉雷达的小天体探测器多模型导航方法
CN113985732B (zh) 针对飞行器系统的自适应神经网络控制方法及装置
CN114923483A (zh) 一种基于状态变换卡尔曼滤波的dvl/sins组合导航方法
CN110595470A (zh) 一种基于外定界椭球集员估计的纯方位目标跟踪方法
CN114323007A (zh) 一种载体运动状态估计方法及装置
CN109581356B (zh) 一种常值机动空间目标的约束滤波追踪方法
Yuan et al. Multisensor Integrated Autonomous Navigation Based on Intelligent Information Fusion
CN117784114B (zh) 异常噪声下基于混合熵的不规则扩展目标跟踪方法
CN111796271B (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