CN109682383B - 一种使用深空三向测量距离和数据的实时滤波定位方法 - Google Patents

一种使用深空三向测量距离和数据的实时滤波定位方法 Download PDF

Info

Publication number
CN109682383B
CN109682383B CN201811407741.6A CN201811407741A CN109682383B CN 109682383 B CN109682383 B CN 109682383B CN 201811407741 A CN201811407741 A CN 201811407741A CN 109682383 B CN109682383 B CN 109682383B
Authority
CN
China
Prior art keywords
time
measurement
data
distance
way
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
CN201811407741.6A
Other languages
English (en)
Other versions
CN109682383A (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.)
China Xian Satellite Control Center
Original Assignee
China Xian Satellite Control Center
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 China Xian Satellite Control Center filed Critical China Xian Satellite Control Center
Priority to CN201811407741.6A priority Critical patent/CN109682383B/zh
Publication of CN109682383A publication Critical patent/CN109682383A/zh
Application granted granted Critical
Publication of CN109682383B publication Critical patent/CN109682383B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种使用深空三向测量距离和数据的实时滤波定位方法，简化了数据预处理过程；原理简单，可操作性强，易推广和使用，通过采用适于非线性系统的免导数UKF滤波算法，省去了计算雅克比矩阵带来的复杂数学推导过程，且外推模型可使用基于动力学的建模方法，也可使用非动力学的统计建模方法，具有较广泛的适用性；在滤波过程中可实时估计测量误差，提高了计算结果的可靠性和稳定性，为其推广提供了条件；本发明的实时滤波定位方法，适用于双向测量及单程测量外测数据定位计算，对单程测量外测数据，只需要在步骤一的外测数据规整过程中，将单程测距*2变为双程测距，发端站与收端站信息给定为相同站即可使用该方法。

Description

一种使用深空三向测量距离和数据的实时滤波定位方法
技术领域
本发明属于航天导航方法技术领域，具体涉及一种使用深空三向 测量距离和数据的实时滤波定位方法。
背景技术
在深空探测器定位数据源中，除了常用的甚长基线干涉测量 (VLBI)数据外，基于扩频及USB等统一测控设备的外测测量也是 一类重要的观测数据。但由于深空探测任务电磁波传播时延较长，以 及地球自转等因素影响，测站有时无法接收到对应的发射信号，使得 传统的双向测距体制应用受限，这时就需要通过不同站的收发进行外 测距离及距离变化率测量。另外，深空目标的USB外测测角数据精 度较差，在定位计算中较少使用，此时就需要通过三个测站的测距观 测来几何确定探测器位置，此种情况一般采用主站发射上行信号，经 应答机相干转发后，两个副站及主站同时接收下行信号的方式进行距 离测量(如图1所示)。在深空任务中，将这种通过一个站发送上行， 另一个站(或多个站)接收下行信号进行距离测量的体制称作三向测 量，而将一个站自发自收测量的情况称作双向测量。
三向测量数据的测量元素包括三向测距(距离和)及三向测速(距 离和变化率)，是深空探测任务的一类重要的观测数据。
在深空探测任务中，三向测量数据是进行实时定位计算的重要数 据源之一。
传统的利用三向测量数据定位计算方法有以下几类：1)几何方 法，先将三向测量数据解算为各站的单程测距、测距变化率，然后利 用三站数据进行几何定位计算；2)滤波方法，也是先将三向测量数 据解算为各站的单程测距、测距变化率，然后使用单程测量数据进行 滤波计算(类似地球附近航天器的外测观测数据的使用)。
这些方法存在的不足主要是：1)首先需要进行单程测距、单程 测距变化率的解算，而不能直接将原始测量数据作为观测值；2)单 程测量量的解算需要较复杂的数据预处理计算，且可能带来一定的精 度损失；3)几何方法需要同时将三站数据插值到同一时间点，且其 计算结果受测量数据的质量影响显著，计算结果曲线光滑性有时较差。
无味卡尔曼滤波(UKF)是一种使用采样策略近似函数非线性分 布的递归贝叶斯估计算法，它以采样变换(UT变换)为基础。滤波 性能与卡尔曼滤波相当，但对于非线性系统，该方法的滤波性能明显 好于EKF。其最大特点是不需要直接对函数进行线性化处理，因而具 有实现过程简单、避免了繁琐的雅可比矩阵求导计算等优点，成为了 实时轨道及定位计算中使用较多的一种方法。
采用适用于非线性动力学过程的免导数无味卡尔曼滤波框架，且 直接以距离和、距离和变化率为观测量进行深空探测器实时定位计算， 其优点是：1)避免了复杂的单程测量量解算的预处理过程，使数据 预处理软件简单可靠；2)减少了数据预处理可能带来的精度损失；3) 借助了UKF滤波算法在非线性状态求解中的优越性，当模型噪声和 测量噪声矩阵取值合适时，采样变换能保证收敛。
发明内容
本发明的目的在于提供一种使用深空三向测量距离和数据的实 时滤波定位方法，不再需要预先解算单程测距、单程测距变化率，能 够简化数据预处理过程。
本发明采用的技术方案为，一种使用深空三向测量距离和数据的 实时滤波定位方法，具体按照以下步骤实施：
步骤1、外测数据的规整及完善，包括测量数据信息及各个测量 信息的有效性标志、收端站站址坐标、发端站站址坐标；
步骤2、建立系统状态方程；
步骤3、通过迭代算法解算数据对应探测器时间和发射端时间；
步骤4、建立观测模型，在滤波计算中直接将观测模型建立在观 测坐标系下，且以测量数据信息作为观测量；
步骤5、采用UKF算法完成位置估算。
本发明的特点还在于：
步骤1测量数据信息包括距离和测量值、距离和变化率测量值。
步骤2具体过程为：在J2000地心惯性系下建立系统状态向量为：
Figure BDA0001877778250000031
则轨道机动过程增广的动力学模型为：
Figure BDA0001877778250000032
式(1)中
Figure BDA0001877778250000033
gE、gM、gS分别为地球引力、 月球引力、太阳光压摄动力产生的加速度，aF为推力加速度矢量，其 大小为a，b为单位质量的质量变化率。
推力加速度矢量aF的计算，由探测器惯性系下的姿态四元数求得 本体系到惯性系转换矩阵，记为M，以及已知的推力方向在探测器本 体系下的方向矢量，记为Fb，即得到惯性系下推力矢量方向为MFb， 然后乘上推力加速度大小a，得到推力加速度矢量。
步骤2另一种实现方法具体过程为：以当前统计模型为积分模型， 并将该积分模型扩充到三维空间，定义系统状态矢量为
Figure BDA0001877778250000041
其中
Figure BDA0001877778250000042
为加速度矢量，则k时刻到k+1时刻的系统状态外推方程为：
Figure BDA0001877778250000043
式(2)中，F(k)为状态转移矩阵，U(k)为输入的状态矩阵，
Figure BDA0001877778250000044
为当前加速度均值，w(k)为离散白噪声序列；
状态转移矩阵的一种自适应表达式为：
Figure BDA0001877778250000045
式中，T为数据采样的时间间隔。
步骤3具体过程为：
第一步、赋初值：
Figure BDA0001877778250000046
式(3)中，tn为探测器时间，c为光速，τ1为上行信号传播时间， τ2为下行信号传播时间；
第二步、计算探测器在tn时间的位置
Figure BDA0001877778250000047
速度
Figure BDA0001877778250000048
计算发端站时间表示为：te-(τ1+τ2)，在该时间点的位置速度 RA(te-(τ1+τ2))，
Figure BDA0001877778250000049
计算收端站在te时刻的位置RB(te)、速度
Figure BDA00018777782500000410
第三步、计算双站测距和：
Figure BDA0001877778250000051
第四步、判断
Figure BDA0001877778250000052
时或满足收敛条件：迭代次数超过 最大门限时退出迭代，其中ε为给定的判断门限；
第五步、若收敛条件不成立，则置
Figure BDA0001877778250000053
Figure BDA0001877778250000054
然后重复上述迭代步骤，直到满足收敛条件，退出，输出相应的 探测器时间、发端站时间。
步骤4建立观测模型具体过程为：根据探测器时间和发射端时间 和发端站站址坐标、收端站站址坐标，计算出两站分别在TS、TE时刻 的地心惯性系位置rS、rE，以及速度矢量
Figure BDA0001877778250000055
由滤波器状态外推 计算出TN时刻探测器的位置rN和速度
Figure BDA0001877778250000056
得到TE时刻的距离和预测 值为：
r＝|rN-rE|+|rN-rS|
距离和变化率预测值为：
Figure BDA0001877778250000057
TE时刻的距离和预测值、距离和变化率预测值即为观测模型。
步骤5具体过程为：
步骤5.1、滤波初始化：tk-1时刻系统位置估计值
Figure BDA0001877778250000058
估计值方 差矩阵Pk-1，系统状态模型噪声矩阵Qk-1，观测模型噪声矩阵Rk-1；
步骤5.2、状态模型外推；
步骤5.3、计算Sigma采样点及权系数、状态先验均值
Figure BDA0001877778250000061
与方差
Figure BDA0001877778250000062
步骤5.4、计算观测预测值、观测先验均值
Figure BDA0001877778250000063
与方差
Figure BDA0001877778250000064
步骤5.5、计算滤波增益矩阵K；
步骤5.6、计算状态及协方差矩阵更新
Figure BDA0001877778250000065
本发明一种使用深空三向测量距离和数据的实时滤波定位方法 有益效果是：
1)本发明的实时滤波定位方法，不再需要预先解算单程测距、 单程测距变化率，简化了数据预处理过程；
2)本发明的实时滤波定位方法，原理简单，可操作性强，易推 广和使用，通过采用适于非线性系统的免导数UKF滤波算法，省去 了计算雅克比矩阵带来的复杂数学推导过程，且外推模型可使用基于 动力学的建模方法(如李恒年的αβ变质量机动过程动力学模型等)， 也可使用非动力学的统计建模方法(如多项式模型、周宏仁的“当前 统计模型”等)，具有较广泛的适用性；在滤波过程中可实时估计测量 误差，提高了计算结果的可靠性和稳定性，为其推广提供了条件；
3)本发明的实时滤波定位方法，适用于双向测量及单程测量外 测数据定位计算，对单程测量外测数据，只需要在步骤一的外测数据 规整过程中，将单程测距*2变为双程测距，发端站与收端站信息给 定为相同站即可使用该方法。
附图说明
图1是本发明一种使用深空三向测量距离和数据的实时滤波定 位方法流程图；
图2是本发明中距离和解算示意图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明一种使用深空三向测量距离和数据的实时滤波定位方法， 如图1所示，具体按照以下步骤实施：
步骤1、外测数据的规整及完善，包括测量数据信息及各个测量 信息的有效性标志、收端站站址坐标、发端站站址坐标；
测量数据信息包括距离和测量值、距离和变化率测量值。
步骤2、建立系统状态方程；
设探测器t时刻的位置矢量为r＝(x,y,z)，速度矢量为
Figure BDA0001877778250000071
探测器在轨飞行过程中，主要受到地月引力、太阳光压摄动力影响， 以及机动过程中发送机推力等其它力的影响。
由于深空探测器飞行过程可能存在在轨运行、动力下降等不同动 力学阶段，有些阶段较难进行动力学建模(或模型复杂)，为了适应 不同动力学运行方式，本发明采用两种较成熟方法进行系统状态建模：
方法一、为了对连续推力作用下的轨道机动过程进行适应，以αβ 变质量机动过程动力学模型为积分模型。
具体过程为：在J2000地心惯性系下建立系统状态向量为：
Figure BDA0001877778250000072
则轨道机动过程增广的动力学模型为：
Figure BDA0001877778250000073
式(1)中
Figure BDA0001877778250000074
gE、gM、gS分别为地球引力、 月球引力、太阳光压摄动力产生的加速度，aF为推力加速度矢量，其 大小为a，b为单位质量的质量变化率。
推力加速度矢量aF的计算，由探测器惯性系下的姿态四元数求得 本体系到惯性系转换矩阵，记为M，以及已知的推力方向在探测器本 体系下的方向矢量，记为Fb，即得到惯性系下推力矢量方向为MFb， 然后乘上推力加速度大小a，得到推力加速度矢量。
方法二、以当前统计模型为积分模型，并将该积分模型扩充到三 维空间。
具体过程为：定义系统状态矢量为
Figure BDA0001877778250000081
其中
Figure BDA0001877778250000082
为加速度 矢量，则k时刻到k+1时刻的系统状态外推方程为：
Figure BDA0001877778250000083
式(2)中，F(k)为状态转移矩阵，U(k)为输入的状态矩阵，
Figure BDA0001877778250000084
为当前加速度均值，w(k)为离散白噪声序列；
状态转移矩阵的一种自适应表达式为：
Figure BDA0001877778250000085
式中，T为数据采样的时间间隔。
步骤3、通过迭代算法解算数据对应探测器时间和发射端时间；
由于三项测量数据的距离和测量值、距离和变化率测量值的时标 (记为TE)均打在信号收端。考虑到电波传播时延，其对应的探测器 时间(记为TN)、发端站时间(记为TS)必不相同，为此，滤波计算 中需要由探测器飞行的当前状态矢量X及滤波状态外推模型进行外 推及迭代计算出探测器时间和发端站时间，如图2所示。
具体过程为：
第一步、赋初值：
Figure BDA0001877778250000091
式(3)中，tn为探测器时间，c为光速，τ1为上行信号传播时间， τ2为下行信号传播时间；
第二步、计算探测器在tn时间的位置
Figure BDA0001877778250000092
速度
Figure BDA0001877778250000093
计算发端站时间表示为：te-(τ1+τ2)，在该时间点的位置速度 RA(te-(τ1+τ2))，
Figure BDA0001877778250000094
计算收端站在te时刻的位置RB(te)、速度
Figure BDA0001877778250000095
第三步、计算双站测距和：
Figure BDA0001877778250000096
第四步、判断
Figure BDA0001877778250000097
时或满足收敛条件：迭代次数超过 最大门限时退出迭代，其中ε为给定的判断门限；
第五步、若收敛条件不成立，则置
Figure BDA0001877778250000098
Figure BDA0001877778250000099
然后重复上述迭代步骤，直到满足收敛条件，退出，输出相应的 探测器时间、发端站时间。
步骤4、建立观测模型，在滤波计算中直接将观测模型建立在观 测坐标系下，且以测量数据信息作为观测量；
三向测量的观测数据主要有距离和测量值(记为ρ)、距离和变 化率测量值(记为
Figure BDA0001877778250000101
)等，为了更好的利用三向测量数据，减少数据 预处理复杂度及精度损失，本发明在滤波计算中直接将观测模型建立 在观测坐标系(如东北天测量坐标系)下，且距离和测量值、距离和 变化率测量值作为观测量，而非传统的单程测距测速。这也只是本发 明区别于其它方法的最显著特征。
具体过程为：根据探测器时间和发射端时间和发端站站址坐标、 收端站站址坐标，计算出两站分别在TS、TE时刻的地心惯性系位置rS、 rE，以及速度矢量
Figure BDA0001877778250000102
由滤波器状态外推计算出TN时刻探测器的 位置rN和速度
Figure BDA0001877778250000103
得到TE时刻的距离和预测值为：
r＝|rN-rE|+|rN-rS|
距离和变化率预测值为：
Figure BDA0001877778250000104
TE时刻的距离和预测值、距离和变化率预测值即为观测模型。
步骤5、采用UKF算法完成位置估算；
考虑到用固定数目的参数逼近状态的条件分布比逼近它的非线 性函数更容易，UKF滤波以确定性的随机量采样策略(Sigma点采样) 和变换Sigma点来表示随机状态的先验和后验分布，分别计算Sigma 点和权值。算法能够保证随机过程在经过非线性函数传递后，仍能表 征随机函数的主要概率分布特征。在UKF滤波计算中不对非线性方 程进行线性化处理，省去了计算雅克比矩阵带来的复杂数学推导过程； 采用UT变换进行Sigma采样，能够保证滤波计算的收敛和精度。
UKF滤波的顺序递推具体过程为：
步骤5.1、滤波初始化：tk-1时刻系统位置估计值
Figure BDA0001877778250000111
估计值方 差矩阵Pk-1，系统状态模型噪声矩阵Qk-1，观测模型噪声矩阵Rk-1；
步骤5.2、状态模型外推；
步骤5.3、计算Sigma采样点及权系数、状态先验均值
Figure BDA0001877778250000112
与方差
Figure BDA0001877778250000113
步骤5.4、计算观测预测值、观测先验均值
Figure BDA0001877778250000114
与方差
Figure BDA0001877778250000115
步骤5.5、计算滤波增益矩阵K；
步骤5.6、计算状态及协方差矩阵更新
Figure BDA0001877778250000116
实施例1，以嫦娥三号近月制动变轨为例，如图1所示，本发明 一种直接使用深空三向测量距离和数据的实时滤波定位方法，具体步 骤如下：
1、对外测数据帧进行规整和补齐：
首先对每帧外测数据补齐发端站站址坐标、收端站站址坐标，便 于后续计算。
2、建立系统状态方程：
设探测器t时刻的位置矢量为r＝(x,y,z)，速度矢量为
Figure BDA0001877778250000117
建立系统变质量机动过程状态方程为
Figure BDA0001877778250000118
其中，g为引力及各类摄动力产生的加速度；a为推力矢量aF的 大小，设由姿态四元数计算的本体系到惯性系转换矩阵为M，已知推 力矢量在探测器本体系下的方向矢量为Fb(由安装情况得到)，可得 到惯性系下推力矢量方向为MFb，推力矢量为aMFb；b为单位质量的 质量变化率(初值可给定为0)。状态外推通过该式积分得到。
3、迭代解算数据对应探测器时间及发端站时间
设当前外测数据帧数据时标为TE(收端)，通过迭代计算发端站 时标TS、对应探测器时标TN，记距离和为ρ2_way，光速为c，上行信 号传播时间为τ1，下行信号传播时间为τ2。
迭代计算步骤为：
(a)赋初值：
Figure BDA0001877778250000121
式(3)中，tn为探测器时间，c为光速，τ1为上行信号传播时间， τ2为下行信号传播时间；
(b)计算探测器在tn时间的位置
Figure BDA0001877778250000122
速度
Figure BDA0001877778250000123
计算发端站时间表示为：te-(τ1+τ2)，在该时间点的位置速度 RA(te-(τ1+τ2))，
Figure BDA0001877778250000124
计算收端站在te时刻的位置RB(te)、速 度
Figure BDA0001877778250000125
(c)计算双站测距和：
Figure BDA0001877778250000126
(d)判断
Figure BDA0001877778250000127
时或满足收敛条件：迭代次数超过最 大门限时退出迭代，其中ε为给定的判断门限；
(e)若收敛条件不成立，则置
Figure BDA0001877778250000128
Figure BDA0001877778250000131
然后重复上述迭代步骤，直到满足收敛条件，退出迭代，输出相 应的探测器时间、发端站时间。
4、建立量测方程
由发端站、收端站大地坐标可计算出两站分别在TS、TE时刻的地 心惯性系位置rS、rE，以及速度矢量
Figure BDA0001877778250000132
又由滤波器状态外推可 计算出TN时刻探测器的位置rN和速度
Figure BDA0001877778250000133
则距离和预测值为
r＝|rN-rE|+|rN-rS|
距离和变化率的预测值为
Figure BDA0001877778250000134
5、采用UKF算法完成位置估计
步骤5.1、滤波初始化：tk-1时刻系统位置估计值
Figure BDA0001877778250000135
估计值方 差矩阵Pk-1，系统状态模型噪声矩阵Qk-1，观测模型噪声矩阵Rk-1；
步骤5.2、状态模型外推；
状态量
Figure BDA0001877778250000136
中r、
Figure BDA0001877778250000137
的初值可由定轨结果外推得到，推力 加速度大小a初始可给定为0，单位质量质量变化率初始可给定为0。
步骤5.3、计算Sigma采样点及权系数、状态先验均值
Figure BDA0001877778250000138
与方差
Figure BDA0001877778250000139
Sigma采样点计算公式为：
Figure BDA00018777782500001310
Figure BDA00018777782500001311
Figure BDA00018777782500001312
式中：λ＝α2(n+κ)-n，n为状态向量维数
时间更新：
状态预测均值：
Figure BDA0001877778250000141
状态协方差：
Figure BDA0001877778250000142
WX为状态均值变换权系数，WP为协方差变换权系数，Q为状态 协方差矩阵。
步骤5.4、计算观测预测值、观测先验均值
Figure BDA0001877778250000143
与方差
Figure BDA0001877778250000144
观测量预测:
Figure BDA0001877778250000145
观测预测均值：
Figure BDA0001877778250000146
观测预测协方差：
Figure BDA0001877778250000147
其中，R观测噪声协方差矩阵；
步骤5.5、计算滤波增益矩阵K；
Figure BDA0001877778250000148
步骤5.6、计算状态及协方差矩阵更新
Figure BDA0001877778250000149
Figure BDA00018777782500001410
Figure BDA00018777782500001411
实施例2，以嫦娥三号月面动力下降和月面软着陆过程为例。
与例1处理过程主要部分相同，不同点在于因动力下降过程动力 学模型较难建模，故采用当前统计模型建立。
即例1第二步建立系统状态方程中，建立系统状态方程修改为
Figure BDA00018777782500001412
其中，F(k)为状态转移矩阵，U(k)为输入的状态矩阵，
Figure BDA00018777782500001413
为当 前加速度均值，w(k)为离散白噪声序列。
其余步骤同例1。
通过上述方式，本发明的实时滤波定位方法，不再需要预先解算 单程测距、单程测距变化率，简化了数据预处理过程；本发明的实时 滤波定位方法，原理简单，可操作性强，易推广和使用，通过采用适 于非线性系统的免导数UKF滤波算法，省去了计算雅克比矩阵带来 的复杂数学推导过程，且外推模型可使用基于动力学的建模方法(如 李恒年的αβ变质量机动过程动力学模型等)，也可使用非动力学的统 计建模方法(如多项式模型、周宏仁的“当前统计模型”等)，具有较 广泛的适用性；在滤波过程中可实时估计测量误差，提高了计算结果 的可靠性和稳定性，为其推广提供了条件；本发明的实时滤波定位方 法，适用于双向测量及单程测量外测数据定位计算，对单程测量外测 数据，只需要在步骤一的外测数据规整过程中，将单程测距*2变为 双程测距，发端站与收端站信息给定为相同站即可使用该方法。

Claims (6)

1.一种使用深空三向测量距离和数据的实时滤波定位方法，其特征在于，具体按照以下步骤实施：
步骤1、外测数据的规整及完善，包括测量数据信息及各个测量信息的有效性标志、收端站站址坐标、发端站站址坐标；
步骤2、建立系统状态方程；
步骤3、通过迭代算法解算数据对应探测器时间和发端站时间；具体过程为：
第一步、赋初值：
Figure FDA0003843703010000011
式(3)中，tn为探测器时间，c为光速，τ1为上行信号传播时间，τ2为下行信号传播时间；
第二步、计算探测器在tn时间的位置
Figure FDA0003843703010000012
速度
Figure FDA0003843703010000013
计算发端站时间表示为：te-(τ1+τ2)，在该时间点的位置RA(te-(τ1+τ2))、速度
Figure FDA0003843703010000014
计算收端站在te时刻的位置RB(te)、速度
Figure FDA0003843703010000015
第三步、计算双站测距和：
Figure FDA0003843703010000016
第四步、判断
Figure FDA0003843703010000017
时，其中ε为给定的判断门限，或满足收敛条件：迭代次数超过最大门限时，退出迭代；
第五步、若收敛条件不成立，则置
Figure FDA0003843703010000021
Figure FDA0003843703010000022
然后重复迭代步骤，直到满足收敛条件，退出，输出相应的探测器时间、发端站时间；
步骤4、建立观测模型，在滤波计算中直接将观测模型建立在观测坐标系下，且以测量数据信息作为观测量；
步骤5、采用UKF算法完成位置估算；具体过程为：
步骤5.1、滤波初始化：tk-1时刻系统位置估计值
Figure FDA0003843703010000023
估计值方差矩阵Pk-1，系统状态模型噪声矩阵Qk-1，观测模型噪声矩阵Rk-1；
步骤5.2、状态模型外推；
步骤5.3、计算Sigma采样点及权系数、状态先验均值
Figure FDA0003843703010000024
与方差
Figure FDA0003843703010000025
步骤5.4、计算观测预测值、观测先验均值
Figure FDA0003843703010000026
与方差
Figure FDA0003843703010000027
步骤5.5、计算滤波增益矩阵K；
步骤5.6、计算状态及协方差矩阵更新
Figure FDA0003843703010000028
2.根据权利要求1所述一种使用深空三向测量距离和数据的实时滤波定位方法，其特征在于，步骤1所述测量数据信息包括距离和测量值、距离和变化率测量值。
3.根据权利要求1所述一种使用深空三向测量距离和数据的实时滤波定位方法，其特征在于，步骤2具体过程为：探测器t时刻的位置矢量为r＝(x,y,z)，速度矢量为
Figure FDA0003843703010000031
在J2000地心惯性系下建立系统状态向量为：
Figure FDA0003843703010000032
则轨道机动过程增广的动力学模型为：
Figure FDA0003843703010000033
式(1)中
Figure FDA0003843703010000034
gE、gM、gS分别为地球引力、月球引力、太阳光压摄动力产生的加速度，aF为推力加速度矢量，其大小为a，b为单位质量的质量变化率。
4.根据权利要求3所述一种使用深空三向测量距离和数据的实时滤波定位方法，其特征在于，所述推力加速度矢量aF的计算，由探测器惯性系下的姿态四元数求得本体系到惯性系转换矩阵，记为M，以及已知的推力方向在探测器本体系下的方向矢量，记为Fb，即得到惯性系下推力矢量方向为MFb，然后乘上推力加速度大小a，得到推力加速度矢量。
5.根据权利要求1所述一种使用深空三向测量距离和数据的实时滤波定位方法，其特征在于，步骤2具体过程为：以当前统计模型为积分模型，并将该积分模型扩充到三维空间，探测器t时刻的位置矢量为r＝(x,y,z)，速度矢量为
Figure FDA0003843703010000035
定义系统状态矢量为
Figure FDA0003843703010000036
其中
Figure FDA0003843703010000037
为加速度矢量，则k时刻到k+1时刻的系统状态外推方程为：
Figure FDA0003843703010000038
式(2)中，F(k)为状态转移矩阵，U(k)为输入的状态矩阵，
Figure FDA0003843703010000041
为当前加速度均值，w(k)为离散白噪声序列；
状态转移矩阵的一种自适应表达式为：
Figure FDA0003843703010000042
式中，T为数据采样的时间间隔。
6.根据权利要求1所述一种使用深空三向测量距离和数据的实时滤波定位方法，其特征在于，步骤4建立观测模型具体过程为：根据探测器时间和发端站时间和发端站站址坐标、收端站站址坐标，计算出两站分别在TS、TE时刻的地心惯性系位置rS、rE，以及速度矢量
Figure FDA0003843703010000043
由滤波器状态外推计算出TN时刻探测器的位置rN和速度
Figure FDA0003843703010000044
得到TE时刻的距离和预测值为：
A＝|rN-rE|+|rN-rS|
距离和变化率预测值为：
Figure FDA0003843703010000045
TE时刻的距离和预测值、距离和变化率预测值即为观测模型。
CN201811407741.6A 2018-11-23 2018-11-23 一种使用深空三向测量距离和数据的实时滤波定位方法 Active CN109682383B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811407741.6A CN109682383B (zh) 2018-11-23 2018-11-23 一种使用深空三向测量距离和数据的实时滤波定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811407741.6A CN109682383B (zh) 2018-11-23 2018-11-23 一种使用深空三向测量距离和数据的实时滤波定位方法

Publications (2)

Publication Number Publication Date
CN109682383A CN109682383A (zh) 2019-04-26
CN109682383B true CN109682383B (zh) 2022-11-04

Family

ID=66185446

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811407741.6A Active CN109682383B (zh) 2018-11-23 2018-11-23 一种使用深空三向测量距离和数据的实时滤波定位方法

Country Status (1)

Country Link
CN (1) CN109682383B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110466790A (zh) * 2019-08-30 2019-11-19 上海电机学院 一种基于机器视觉的无人机目标跟踪系统
CN113569430B (zh) * 2021-08-31 2023-07-04 中国西安卫星测控中心 一种仅外测观测下再入飞行转弯方向辨识方法
CN114061622B (zh) * 2021-11-11 2023-08-22 中国西安卫星测控中心 一种深空三向测距系统误差标定方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104864875A (zh) * 2015-04-03 2015-08-26 北京控制工程研究所 一种基于非线性h∞滤波的航天器自主定位方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0175967B1 (en) * 1984-09-07 1990-10-24 Deutsches Zentrum für Luft- und Raumfahrt e.V. Navigation, communication, and surveillance system based on dme
US6847587B2 (en) * 2002-08-07 2005-01-25 Frank K. Patterson System and method for identifying and locating an acoustic event
US8914038B2 (en) * 2009-02-11 2014-12-16 Telefonaktiebolaget L M Ericsson (Publ) Method and arrangement for determining terminal position
CN104219620B (zh) * 2013-06-05 2018-11-06 华为技术有限公司 一种终端的定位方法和装置
CN103424116B (zh) * 2013-07-23 2015-09-23 中国西安卫星测控中心 一种适应轨道机动的地球同步卫星精密定轨方法
CN105158817B (zh) * 2015-08-04 2017-08-25 中国科学院上海天文台 深空探测器多普勒频率被动式测量方法
CN107421550B (zh) * 2017-07-25 2020-08-28 北京航空航天大学 一种基于星间测距的地球-Lagrange联合星座自主定轨方法
CN108844536B (zh) * 2018-04-04 2020-07-03 中国科学院国家空间科学中心 一种基于量测噪声协方差矩阵估计的地磁导航方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104864875A (zh) * 2015-04-03 2015-08-26 北京控制工程研究所 一种基于非线性h∞滤波的航天器自主定位方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Jiang Liu ; Mingquan Lu.An Adaptive UKF Filtering Algorithm for GPS Position Estimation.《2009 5th International Conference on Wireless Communications, Networking and Mobile Computing》.2009,全文. *
一种火箭及上面级外弹道实时滤波算法;淡鹏等;《雷达科学与技术》;20161031;第14卷(第05期);全文 *
三向测量技术在深空探测中的应用研究;黄磊，王宏，樊敏;《飞行器测控学报》;20120630;第31卷(第3期);全文 *

Also Published As

Publication number Publication date
CN109682383A (zh) 2019-04-26

Similar Documents

Publication Publication Date Title
Liu et al. Maximum correntropy square-root cubature Kalman filter with application to SINS/GPS integrated systems
CN109682383B (zh) 一种使用深空三向测量距离和数据的实时滤波定位方法
US8639476B2 (en) Process for estimation of ballistic missile boost state
CN109323698B (zh) 一种空间目标陨落多模型跟踪引导方法
CN108255791B (zh) 基于分布式传感器一致性的机动目标跟踪方法
Li et al. A machine learning-based approach for improved orbit predictions of LEO space debris with sparse tracking data from a single station
CN107270933B (zh) 一种基于多星协同的空间碎片运动状态联合确定方法
US7725259B2 (en) Trajectory estimation system for an orbiting satellite
CN110779518A (zh) 一种具有全局收敛性的水下航行器单信标定位方法
CN105203101A (zh) 一种基于目标天体星历修正的深空探测器捕获段天文导航方法
Biswas et al. A novel a priori state computation strategy for the unscented Kalman filter to improve computational efficiency
CN110779519A (zh) 一种具有全局收敛性的水下航行器单信标定位方法
CN104048676A (zh) 基于改进粒子滤波的mems陀螺随机误差补偿方法
CN108490973B (zh) 航天器编队相对轨道确定方法及装置
CA2699137A1 (en) Hybrid inertial system with non-linear behaviour and associated method of hybridization by multi-hypothesis filtering
CN114510076A (zh) 基于无迹变换的目标协同探测与制导一体化方法及系统
CN104932266B (zh) 一种基于前馈补偿的着陆器进入段精确控制方法
CN110146092B (zh) 基于导航信息评价的双体小行星探测轨迹优化方法
Ascorti An application of the extended Kalman filter to the attitude control of a quadrotor
CN107483131B (zh) 高速飞行器双卫星联合信道马尔科夫状态序列产生方法
CN112346010A (zh) 基于尺度差和时差的双机无源定位方法
Hough Improved performance of recursive tracking filters using batch initialization and process noise adaptation
Karlgaard et al. Mars Science Laboratory entry, descent, and landing trajectory and atmosphere reconstruction
CN115859839A (zh) 基于rcs序列阶跃效应的抛物面天线载荷指向估计方法
CN110132283A (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