CN106643744A - 一种基于四程中继跟踪模式的远月面着陆器精密定位方法 - Google Patents

Abstract

本发明提供一种基于四程中继跟踪模式的远月面着陆器精密定位方法，基于观测模型实现远月面着陆器精密定位，所述四程中继跟踪测量模式，将月球轨道器作为中继星，使用月球着陆器和月球轨道器之间的链路进行测量。本发明克服传统直接跟踪测量模式无法对远月面着陆器进行定位的缺陷，提供一种基于四程中继跟踪模式的远月面着陆器精密定位方法，借助于中继星的中继跟踪测量，可以消除月球自身的遮挡，实现对远月面着陆器的精密定位。由于该方式对着陆器以及中继星有很强的几何约束，将使精密定轨、定位的精度得到极大地提高。

Description

一种基于四程中继跟踪模式的远月面着陆器精密定位方法

技术领域

本发明属于月球着陆器精密定位与深空探测领域，特别是涉及一种基于四程中继跟踪模式的远月面着陆器精密定位方法。

背景技术

目前，深空探测是世界各国竞相争逐的热点，进行深空探测，可以进一步解答地球如何起源与演变、行星和太阳系形成和演化的过程、人类是不是宇宙中唯一的智慧生命、地球的未来将如何等一系列问题，同时有利于人类积极开发和利用空间资源。在深空探测中，探测器的定轨定位是任务成败的关键，也是各种科学任务顺利进行的前提。精密的探测器轨道是进行地形地貌测绘的基础数据，也可以用于行星重力场的解算，进而反演行星的内部构造。

月球背面的着陆探测一直是国际上的热点和难点。月球背面具有不同于月球正面的独特地质构造，是研究地月起源的重要突破口；月球背面没有地球电磁波干扰，是进行低频射电天文观测的天然理想场所。但是由于月球自转与公转同步，月球背面着陆器无法与地球测站进行直接通信，任务难度和风险比较大，目前对月球背面的就位探测仍旧是一项国际空白。因此，开展月球背面的就位探测具有重要的科学意义和工程意义。

对于近月面着陆器的定位问题，一般采用运动学统计定位方法，采用测量模式包括：双程/三程测距、测速，VLBI时延/时延率模型等等。在“嫦娥三号”着陆器的精密定位中使用了运动学统计定位方法，综合了测距、测速数据和VLBI时延、时延率数据，着陆器定位的绝对精度在10m左右(曹建峰,张宇,胡松杰,等.2016.嫦娥三号着陆器精确定位与精度分析.武汉大学学报(信息科学版),41(2):274-278.doi:10.13203/j.whugis20140123.)，“玉兔号”漫游器与着陆器的相对定位采用了同波束VLBI技术，相对位置精度可达到米级(黄勇,昌胜骐,李培佳,等.2014.“嫦娥三号”月球探测器的轨道确定和月面定位.科学通报,59(23):2268-2277.)。

但是对于月球背面着陆器的精密定位问题，如图1所示，图1a为近月面示意图，图1b为远月面示意图，由于月球自身的遮挡，远月面着陆器与地球深空测控站无法通视，传统的直接测量模式，如“嫦娥三号”中的双程、三程测距/测速，VLBI时延/时延率测量模式将不再适用，因此有必要采取新的跟踪测量模式对远月面着陆器进行定位。

在对月球着陆器定位的过程中，涉及到坐标系的转换。地面测站位于地球，一般用地球参考框架描述；月球着陆器需要用月球主轴坐标系来表示；而信号在太空中的传输要考虑大天体引力时延以及相对论效应的影响，最终均需要在太阳系质心坐标系下(BCRS)表示。为了得到高精度的定位结果，需要不同坐标系之间进行高精度的转换。各坐标系之间的转换如图2所示。

地面测站坐标由协议地固系(ITRS)转为J2000地球惯性坐标系下的坐标，此转换需要岁差旋转、章动旋转、极移旋转以及地球自转旋转。关于岁差、章动模型，随着观测值精度的不断提高，也在不断地修正。目前IERS2010规范推荐采用IAU2006决议的岁差章动模型。J2000地球惯性系坐标经过平移、洛伦兹变换可以转到太阳系质心坐标系下。

同地固坐标系统的转换中用到岁差章动和极移一样，从月球主轴坐标系(PA)到月心天球坐标系，要考虑月球天平动(Libration)的影响。目前用于月球物理天平动研究的主要手段仍然是对月激光测距(LLR)，由月球自转引起测距的变化，通过数值积分来求得月球自转相应的三个欧拉角，通过旋转三个欧拉角即可转换到月心天球坐标系下，进一步经过平移、洛伦兹变换可以转到太阳系质心坐标系(BCRS)下。

目前对月球背面的着陆一直是国际空白，月球背面的着陆最大的难题是克服信号被月球遮挡的影响，同时面临无法实时测控跟踪着陆器等难题。

发明内容

本发明的目的是克服传统直接跟踪测量模式无法对远月面着陆器进行定位的缺陷，提供一种基于四程中继跟踪模式的远月面着陆器精密定位方法，该方法借助于中继星的中继跟踪测量，可以消除月球自身的遮挡，实现对远月面着陆器的精密定位。

为达到上述发明目的，本发明的技术方案提供一种基于四程中继跟踪模式的远月面着陆器精密定位方法，建立月球着陆器-轨道器的四程中继跟踪测量模式的观测模型，基于观测模型实现远月面着陆器精密定位，

所述四程中继跟踪测量模式，将月球轨道器作为中继星，使用月球着陆器和月球轨道器之间的链路进行测量，设地球跟踪站T_i在时刻i发送一个上行的信号给月球轨道器，月球轨道器位置为S_j，时刻为j；经过转发，信号传送给月面着陆器，相应位置为L_k，时刻为k，L_k位于月球背面；月面着陆器再次转发给在轨的月球轨道器，此时轨道器已经从S_j的位置运动至S_m，时刻为m；之后该信号下行至地球跟踪站T_n，时刻为n；设时刻n、m、k、j、i分别相应的参与者记为T_n,S_m,L_k,S_j,T_i，则有以下4个几何距离，

R₁＝|X(S_m)-X(T_n)|

R₂＝|X(S_m)-X(L_k)|

R₃＝|X(S_j)-X(L_k)|

R₄＝|X(S_j)-X(T_i)|

其中，X(j)为参与者j在太阳系质心坐标系下下的位置矢量；

四程中继测距的观测模型建立如下，

R＝(R₁+c·RLT_nm)+R₂+R₃+(R₄+c·RLT_ij)+

c·[TDB(i)-UTC(i)]-c·[TDB(n)-UTC(n)]

＝c·[UTC(n)-UTC(i)]

其中，RLT_nm为由参与者T_n发射信号到S_m的广义相对论时延，RLT_ij为由参与者T_i发射信号到S_j的广义相对论时延，UTC(n)和UTC(i)分别为时刻n和时刻i对应的协调世界时时标，R为四程中继测距值，TDB(n)和TDB(i)分别为时刻n和时刻i对应的太阳系质心力学时，c为光速；

设在一个多普勒积分周期内，起始时刻T_s和终止时刻T_e各自对应的四程中继测距值为R_s、R_e，建立四程中继测速的观测模型如下，

其中，T_c为一个多普勒积分周期，RR为四程中继测速值。

而且，根据四程中继跟踪测量模式，将月球轨道器的轨道和远月面着陆器的位置同时视为待估参数，采用精密定轨的方式进行解算，最后得到两者的最佳估值。

而且，采用精密定轨的方式进行解算时，

计算四程中继测距值R对远月面着陆器坐标X(L_k)的偏导数H_R公式如下，

按照上式计算积分起始时刻T_s和积分终止时刻T_e分别对应的四程中继测距值R对远月面着陆器坐标X(L_k)的偏导数H_s和H_e，则四程中继测速值RR对远月面着陆器坐标X(L_k)的偏导数H_RR为

将结果作为观测时刻的观测值对状态向量的偏导数。

而且，将探测器轨道跟踪数据文件分为多个观测弧段，将初始参考轨道x₀ ^ref＝[r₀v₀ p₀ q₀]^T中的中继星的初始位置r₀、初始速度v₀以及力模型参数p₀作为局部参数处理，记作x₁，维数为m₁；参数q₀是远月面着陆器位置坐标，将其作为全局参数处理，记作x₂，维数为3；采用精密定轨的方式进行解算时，融合多个观测弧段相应法方程，将非线性的方程线性化，使用迭代法求解。

本发明的有益效果：本发明基于四程中继跟踪模式，利用一颗中继星，实现对远月面着陆器的事后精密定位；由于基于四程中继跟踪模式建立了中继星和远月面着陆器间的联系，解决了地球深空站无法直接跟踪测量远月面着陆器的问题。由于该方式对着陆器以及中继星有很强的几何约束，将使精密定轨、定位的精度得到极大地提高。此外，随着此模式观测数据的增加，可将用于支持对月球自转欧拉角，月球K2勒夫数以及月球重力场等反演，得到更为精确的月球物理模型，进而为了解月球的起源与演化，深入研究月球内部构造提供重要的科学数据。本发明还提出将远月球着陆器坐标作为全局参数进行解算，理论上更为严密，可得到高精度的远月面着陆器坐标。因此，本发明具有重要的科研价值和市场价值，是我国领先世界的重要研究项目。

附图说明

图1为传统的双程/三程测距、测速及VLBI模型示意图，其中图1a为近月面示意图，图1b为远月面示意图；

图2为本发明涉及到的各类坐标系之间的转换示意图；

图3为本发明实施例的基于月球着陆器-轨道器的四程中继跟踪测量模型示意图。

图4为本发明实施例的解算流程示意图。

具体实施方式

以下结合附图和实施例对本发明技术方案进行具体描述。

为便于实施参考起见，首先分别介绍本发明涉及的精密定轨原理和基于月球着陆器-轨道器的四程中继跟踪测量模式(4W L-O)的基本原理：

(1)精密定轨原理

本发明中涉及到的精密定轨理论不考虑探测器搭乘火箭被送入指定轨道的过程，而是进入指定轨道后，通过探测器自身发射无线电信号，建立与地面(或空间站和卫星)跟踪测控系统的联系，从而对原始轨道进行调控和微分改进。具体阐述如下：

第一，根据牛顿运动定律，探测器(本文指月球轨道器)进入指定轨道(初始轨道)后，受到宇宙中各种力(天体引力、太阳光压力等)的影响将会产生加速度，从而改变原有的初始运动状态(初始轨道)。该过程可用一个常微分方程组描述，称之为力模型。但由于太空环境复杂且人类目前的认知程度有限，不可能用该方程准确无误地描述出探测器的实际受力情况，所以该方程组存在误差即力模型存在误差，具体体现在两方面：第一，力模型中的已知参数不够准确；第二，实际受力中还存在未建模的部分。对于第二点，需要整个科学的推动才能对受力有进一步的认识，所以力模型的误差主要指第一点。

第二，如果不能对探测器保持联系，任何探测器都会很快失效。所以，探测器将会主动发射无线电信号，由地面(或空间站、卫星)跟踪站的跟踪测量设备进行接收。根据收发信号链路的几何关系，可以建立相应的观测模型。信号在实际传播过程中，会受到复杂空间环境、相对论效应等的影响，该影响通常作为改正项融入观测模型中，可以视为观测模型的模型参数。同样，由于对这些影响的认知有限等原因，这些模型参数也存在误差。同时，跟踪站的观测设备本身存在一定的误差，所以获取的观测资料也是存在误差的。

结合上述两点，精密定轨实际上可以这样描述：利用带有大量带有误差的观测数据和并非精确的运动学微分方程，使用统计学原理解算在某种意义下的探测器初始轨道和相关模型(力模型、观测模型)参数的最优估值的过程。

由于一般将观测误差视为正态分布的随机误差，所以精密定轨实质上为基于最小二乘的参数估计问题，即：寻找初始轨道和相关模型(力模型、观测模型)参数的估计值，使得利用估计值计算的观测残差平方和最小。

其中，第一点中所述力模型的形式如下：

式(1)是一个描述力模型的常微分方程组，其中：x(t)为随时间t变化的m维状态向量，包含探测器的位置r(t)、速度v(t)、力模型相关参数p以及与观测模型相关的参数q，m表示以上这些参数的个数的总和；为随时间变化的m维状态向量的一阶导数，f(t,x(t))为的具体函数形式，x₀为历元初始时刻t₀不精确的初值x(t₀)，包含探测器的初始位置r₀(t₀)、初始速度v₀(t₀)、(以下简写为r₀和v₀)力模型相关初始参数p₀以及与观测模型相关的参数q₀，如式(2)所示。

其中，

令z表示n维观测向量，即：

z_i(t_i)表示第i个观测量，即在观测时间t_i获取的观测量，i＝1,2,…n，n表示观测弧段内的观测次数；z_i(t_i)用于表示第i个观测量与探测器状态的关系如下：

z_i(t_i)＝g_i(t_i,x(t_i))+ε_i＝h_i(t_i,x₀)+ε_i (4)

式(4)中：g_i表示观测量与观测时间t_i及该时刻的探测器状态x(t_i)的函数关系；h_i()表示观测量与观测时间t_i及历元初始时刻t₀的航天器状态的函数关系；ε_i为在观测时间t_i由于误差造成的理论观测量与实际观测量的偏差。式(5)写为向量表达式如下：

z＝h(x₀)+ε (5)

其中，z,h(x₀),ε分别为式(4)中所述的z_i(t_i),h_i(t_i,x₀),ε_i的向量形式；一方面，h(x₀)是一个关于未知矢量x₀的非线性函数，利用最小二乘方法求解轨道确定问题非常复杂；另一方面，真实的矢量x₀虽然未知，但是实际给出的初始探测器状态可以视为探测器真实初始状态x₀的近似值，称之为初始参考轨道，记为x₀ ^ref。在该点进行泰勒展开，并忽略高阶项，称之为线性化，可得：

写为误差方程的形式，有：

ε＝Δz-HΔx₀ (7)

Δx₀＝x₀-x₀ ^ref (8)

Δz＝z-h(x₀ ^ref) (9)

式(6)至式(10)中：△x₀为真实初始轨道x₀与不精确的初始参考轨道x₀ ^ref的偏差；△z为实际观测量z与利用该不精确的参考轨道x₀ ^ref计算的观测量的残差；H为利用历元t₀时刻参考轨道初值(即初始参考轨道x₀ ^ref)计算的观测量相对于历元t₀时刻状态向量的偏导数，x(t_i)为历元t_i时刻的参考轨道，H可以利用链式法则进行展开成(10)式所示的形式。是t_i时刻的观测量对状态向量x(t_i)的偏导数。

通过线性化，精密定轨问题可以简化为线性最小二乘问题。获取状态改正量△x₀的最优估值△x₀ ^lsq，使得经状态改正△x₀后的目标函数J(Δx₀)取最小值。

J(Δx₀)＝ε^Tε＝(Δz-HΔx₀)^T(Δz-HΔx₀) (11)

根据线性方程理论，如果H满秩，则目标函数取极小值的条件为

计算上述偏导数，可得法方程形式如下：

(H^TH)Δx₀＝H^TΔz (13)

解得线性最小二乘问题解Δx₀ ^lsq的一般形式：

Δx₀ ^lsq＝(H^TH)^-1(H^TΔz) (14)

矩阵H^TH是一个m维对称矩阵，通常称为法方程矩阵。

如果在轨道估计之前，在已有待估参数△x₀初始值的前提下，同时具备待估参数△x₀的验前信息

式(15)中，E(△x₀)表示对△x₀取期望，cov(△x₀-x₀)表示对某个随机变量序列取方差-协方差矩阵；给出的先验信息即为待估参数△x₀的期望x₀和方差-协方差矩阵P₀。引入变量如下，

这样式(11)所示的目标函数就可以写为

相应的法方程可以写为

其最小二乘解Δx₀ ^lsq为

则修正后的航天器状态x₀ ^lsq为

x₀ ^lsq＝x₀ ^ref+Δx₀ ^lsq (20)

对于不同观测类型，测量数据精度各有差异。而上述式中将观测数据按相同观测精度处理，这在实际应用不尽合理，故需要对数据进行加权，相应的权矩阵可以写为如下形式：

W＝diag(σ₁ ^-2,…,σ_n ^-2) (21)

σ_i为互不相关的测量误差，i＝1,2,…n。当误差相互关联时，权矩阵为非对角矩阵。相应的法方程的形式为：

式(18)的最小二乘解为如下形式：

当不考虑先验信息时，P₀ ^-1＝0；相应的法方程变为：

(H^TWH)Δx₀＝H^TWΔz (24)

相应的最小二乘解为：

Δx₀ ^lsq＝(H^TWH)^-1(H^TWΔz) (25)

其中，式中的H按照式(10)进行计算；式中的△z按照式(9)计算。

由于对非线性方程进行了线性化，为了保证最终解的准确性，根据精密定轨理论，整个过程需要迭代进行。即将解算得到的x₀ ^lsq重新作为x₀ ^ref，重复式(1)至式(25)，直至与上一次得到的x₀ ^lsq之差的绝对值小于设定的限差δ即可。

故本发明可概括为：结合一种新的观测模型——四程中继跟踪测量模型，将月球轨道器(中继星)的轨道和远月面着陆器的位置同时视为待估参数，按照精密定轨的流程进行解算，最后得到两者的最佳估值的过程。

(2)月球着陆器-轨道器的四程中继跟踪测量模式(4W L-O)原理

和传统双程模式相比，基于月球着陆器-轨道器的四程中继跟踪模式使用了着陆器和轨道器之间的链路。如图3所示，地球跟踪站T_i(时刻i)发送一个上行信号给月球轨道器(中继星，下文轨道器均指中继星)(位置为S_j，时刻j)，经过转发，该信号传送给月面着陆器(位置为L_k，时刻k，L_k位于月球背面)，着陆器再次转发给在轨的月球轨道器(此时轨道器已经从S_j的位置运动至S_m，时刻为m)，之后该信号下行至地球跟踪站T_n(时刻n)。上述过程就是基于月球着陆器-轨道器的四程中继跟踪模式信号收发的一个完整过程。在实际测量中，地球跟踪站T_n可以通过积分多普勒的方式得到4W L-O的距离变率，和双程测速的模式类似。对于4W L-O测距和测速的理论值C的计算，和传统的双程测距测速类似，采用“回溯法”，即按照T_n-＞S_m-＞L_k-＞S_j-＞T_i的路径，为便于计算，可视为有五个参与者T_n,S_m,L_k,S_j,T_i。地球跟踪站T_n在n时刻接收到下行信号并记录，但m，k，j，i时刻的相关参与者的状态向量是未知的，需要进行光行时解算依次获取。其中，光行时解算公式如下：

式(26)和式(27)中，c为光速；下标a为该支路的信号发射方，下标b为该支路的信号接收方；t_b表示该支路的信号接收方的接收时刻，t_a表示该支路的信号发射方的发射时刻；对于4W L-O模式的四条链路，a＝i时b＝j，a＝j时b＝k，a＝k时b＝m，a＝m时b＝n。式(26)中，r_ab是BCRS(Barycentric Celestial Reference System)坐标系下该支路信号发射方和接收方的几何距离，是该支路信号传播过程中的牛顿光行时，RLT_ab是该支路信号传播过程中的广义相对论时延。式(27)中，μ_s是太阳引力常数，μ_B是行星、外行星系统或月球的引力常数。γ为后牛顿参数，和代表该支路信号发射方和接收方分别在对应时刻的离太阳质心的几何距离，和代表该支路信号发射方和接收方在对应时刻的离行星、外行星系统或月球质心的几何距离。其中和代表该支路信号接收方和发射方在对应时刻的日心坐标向量，和代表该支路信号接收方和发射方在对应时刻的行星质心坐标向量。代表在日心坐标系下该支路信号接收方和发射方的几何距离，代表在行星质心坐标系下该支路信号接收方和发射方的几何距离。

本发明利用提出的月球着陆器-轨道器的四程中继跟踪观测模型，同时结合传统的双程、三程测距测速观测模式，实现中继星的精密定轨和远月面着陆器的精密定位。具体步骤如下：

根据式(2)，设与轨道器相关的初始状态的近似值为：x₀ ^ref＝[r₀ v₀ p₀ q₀]^T。其中，x₀ ^ref为随时间变化的m维状态向量，包含轨道器的初始位置r₀、速度v₀、力模型相关参数p₀以及与观测模型相关的参数q₀；其中q₀为远月面着陆器的初始近似位置坐标同时，设置迭代收敛的限差δ，在设置迭代收敛的限差时，分别设置局部参数和全局参数迭代收敛的限差δ₁和δ₂。

①根据轨道器的实际受力情况，同时在给定的观测站的观测弧段内，选定x₀ ^ref，并将式(1)所示的微分方程组进行数值积分，解算出一条以x₀ ^ref为初始条件且并不精确的参考轨道x^ref(t_i)，其中x^ref(t_i)代表观测弧段内某时刻t_i的参考轨道x^ref，i＝1,2,…n。同时，式(10)所示的也可以由式(1)所示的微分方程组数值积分得来。具体实施时，初始值、微分方程组的具体细节以及限差可由本领域技术人员自行预先设定。

②建立月球着陆器-轨道器的四程中继跟踪测量模式(4W L-O)的观测模型

本步骤相当于计算式(6)中的h(x₀ ^ref)。

四程中继测距(4W L-O R)：

将下行信号到达T_n的时刻n记录为协调世界时时标UTC(n)，协调世界时通常简称UTC，首先需要将该UTC时标转换为太阳系质心力学时标TDB，可得到下行信号到达T_n的时刻的太阳系质心力学时标TDB(n)。之后通过光行时解依次得到S_m,L_k,S_j,T_i处对应的太阳系质心力学时TDB(m)，TDB(k)，TDB(j)，TDB(i)，最后再次通过时间转换得到T_i处的协调世界时时标UTC(i)。在进行各个链路的光行时解算时，如式(26)和(27)所示，轨道器与地面之间的链路考虑广义相对论时延项，轨道器与着陆器之间的链路由于距离很近，只考虑牛顿光行时，不考虑广义相对论时延项。通过此“回溯”，可以得到相应时刻各个参与者的状态向量。定义X(J)为参与者J在BCRS下的位置矢量，J＝T_n,S_m,L_k,S_j,T_i，那么4个几何距离依次为：

R₁＝|X(S_m)-X(T_n)| (28)

R₂＝|X(S_m)-X(L_k)| (29)

R₃＝|X(S_j)-X(L_k)| (30)

R₄＝|X(S_j)-X(T_i)| (31)

以上可作为四程测距测量方程。

至此，四程中继测距(4W L-O R)的观测模型建立如下：

其中，RLT_nm和RLT_ij为式(26)和式(27)中的广义相对论时延，即RLT_nm为由参与者T_n发射信号到S_m的广义相对论时延，RLT_ij为由参与者T_i发射信号到S_j的广义相对论时延；UTC(n)和UTC(i)分别为上述T_n相应时刻n和T_i相应时刻i对应的协调世界时时标，TDB(n)和TDB(i)分别为时刻n和时刻i对应的太阳系质心力学时，R为四程中继测距值，即基于月球着陆器-轨道器的四程中继测距理论计算值。

四程中继测速(4W L-O RR)：

基于月球着陆器-轨道器的四程中继测速理论值按照积分多普勒的形式可表达为：在一个多普勒积分周期内，起始时刻T_s和终止时刻T_e各自对应的四程中继测距值R_s、R_e之差相对于时间的变化率。首先将积分起始时刻T_s视为下行信号的到达时刻T_n，按照式(32)所示的观测模型计算四程中继测距值R_s，然后将积分终止时刻T_e同样视为下行信号的到达时刻T_n，按照式(32)所示的观测模型计算四程中继测距值R_e。建立四程中继测速(4W L-ORR)的观测模型为：

上式中T_c为一个多普勒积分周期，等于T_e-T_s，RR为四程中继测速值，即基于月球着陆器-轨道器的四程中继测速理论值。

③计算观测模型所得四程中继测距观测量(包括四程中继测距值R、四程中继测速值RR)对远月面着陆器坐标的偏导数

本步骤相当于计算式(10)中的

四程中继测距值R对远月面着陆器坐标X(L_k)的偏导数H_R：

按照式(34)计算积分起始时刻T_s和积分终止时刻T_e分别对应的四程中继测距值R对远月面着陆器坐标X(L_k)的偏导数H_s和H_e，则四程中继测速值RR对远月面着陆器坐标X(L_k)的偏导数H_RR为：

计算所得偏导数H_RR即偏导数

④输入实际观测值序列(即n维观测向量z)，按照式(9)计算Δz，并按照式(10)计算H矩阵。

⑤结合传统的双程/三程测距、测速等测量模型，结合实际测量精度，按照式(21)确定权矩阵W；其中，W也可以由本领域技术人员根据不同的算法(如：赫尔模特方差分量估计法)预先确定。

⑥生成如式(24)所示的法方程，根据式(25)解算出改正量Δx₀ ^lsq，根据式(20)计算出最佳估值，即修正后的航天器状态x₀ ^lsq。其中，x₀ ^lsq向量中包含有远月面着陆器的位置坐标。

⑦将解算出的x₀ ^lsq作为①中的x₀ ^ref，重复①～⑥，直至两次得到的Δx₀ ^lsq之差的绝对值小于预设的限差δ时停止迭代，得到最终的最佳估值x₀ ^lsq。

本发明进一步提出，在实际情况中，由于地球测站可视范围有限，探测器变轨等原因，不能保证永远有连续且有效的观测跟踪数据，所以探测器轨道跟踪数据文件可分为多个观测弧段；同时，根据待估参数本身的性质，一般将x₀ ^ref＝[r₀ v₀ p₀ q₀]^T中的r₀v₀p₀作为局部参数处理，记作x₁，维数为m₁，其中局部参数的特点是已知其每个弧段的初始近似值x₁ ^ref(i)，并解算每个弧段的参数的改正值Δx₁ ^lsq(i)，最终得到每个弧段的参数的最佳估值x₁ ^lsq(i)，其中，i＝n,n-1,…,1，代表第i个弧段；对于参数q₀，本发明中指远月面着陆器位置坐标X(L_k)，将其作为全局参数处理，记作x₂，维数为3，其中全局参数的特点是与弧段无关，解算时融合多个弧段获得最终的全局参数的改正值Δx₂ ^lsq以及最佳估值x₂ ^lsq。具体阐述如下：

对于每一个弧段，一般取1-2天的弧长。输入某个弧段的初始近似值x₀ ^ref，综合大量的四程中继跟踪测量值和传统的双程/三程测距、测速等测量值，对中继星(轨道器)进行精密定轨，在不考虑参数先验信息的前提下，按照步骤①～⑥，根据四程中继跟踪测量模式，针对单弧段中继星精密定轨生成对着陆器坐标的偏导数，生成如式(24)所示的单弧段法方程，为了区分局部参数和全局参数，将法方程矩阵按照局部参数和全局参数的维数分块，可得如下形式：

Δx₁和Δx₂分别代表待估的局部参数(包括中继星初轨、力模型参数，中继星初轨此处为位置和速度)的改正值和全局参数(即远月面着陆器的位置坐标)的改正值，A₁₁对应为局部参数的法方程矩阵，具体形式对应式(24)，为m₁×m₁维的矩阵A₂₂对应为全局参数的法方程矩阵，具体形式对应式(24)，为3×3的矩阵A₂₁和A₁₂为矩阵分块后的左下部分和右上部分；b₁和b₂按照式(24)所示对应为m₁×1维的矩阵和3×1维的矩阵

为了更有效的解算待估参数，联合多个弧段按照步骤①～⑥形成法方程，设有N个弧段可得到类似式(36)的一个序列的弧段法方程：

其中上标(1)～(N)表示弧段编号1～N。式(37)中各符号的具体含义与式(36)相同。

融合多个弧段的法方程，形成如下融合法方程：

对上式中的联合矩阵每个弧段部分采取矩阵行变换，可得到如下形式：

式(38)和式(39)中的各项含义与式(36)相同，其中上标(1)～(N)表示弧段编号1～N。I代表单位矩阵，从下至上解算式(39)，首先解得全局参数向量的改正值Δx₂ ^lsq，即远月面着陆器坐标的改正值，按照式(20)加到远月面着陆器坐标的初始近似值x₂ ^ref得到最优估值x₂ ^lsq。解算出Δx₂ ^lsq后，对式(39)从下至上依次回代，即可按顺序解算出局部参数向量按照式(20)加到局部参数初始近似值x₁ ^ref(n),x₁ ^ref(n-1),…,x₁ ^ref(2),x₁ ^ref(1)，即可得各个弧段局部参数对应的最优估值x₁ ^lsq(n),x₁ ^lsq(n-1),…,x₁ ^lsq(2),x₁ ^lsq(1)。

由于上述过程基于将非线性的方程线性化，为更准确的解算参数，数学上通常将线性化后的非线性方程使用迭代法求解。参见图4，实施例提供的具体流程阐述如下：

第一，对于第一个弧段，按照①～⑥形成如式(36)所示的法方程。

第二，对于第2～N个弧段，按照①～⑥与第一个弧段按照上一步骤形成的法方程一起形成如式(37)所示的法方程。

第三，按照式(37)和式(38)融合多弧段法方程，形成如式(39)所示的形式；即生成弧段(1)～(N)对中继星初轨、力模型参数(局部参数)的偏导数，生成对着陆器坐标(全局参数)的偏导数，可生成弧段(1)～(N)的法方程，开始融合弧段(1)～(N)的法方程进行着陆器坐标(全局参数)解算。

第四，解算出全局参数向量的第一次的改正值Δx₂ ^lsq，按照式(20)加到远月面着陆器坐标的初始近似值x₂ ^ref得到第一次的最优估值x₂ ^lsq。

第五，将当前的最优估值x₂ ^lsq重新作为初始近似值x₂ ^ref，并按照①～⑥的步骤更新第三步中的法方程中的矩阵形成新的多弧段融合的法方程。其中，i＝N,N-1,…1，表示第i个弧段。第一次执行第五步时，当前的最优估值即采用第四步中第一次得到的全局参数的最优估值x₂ ^lsq。

第六，根据新的法方程解算出新的全局参数向量的改正值Δx₂ ^lsq，按照式(20)加到远月面着陆器坐标的初始近似值x₂ ^ref得到新的最优估值x₂ ^lsq。

第七，重复第五步和第六步，直至当次迭代和上一次迭代得到的Δx₂ ^lsq之差的绝对值小于预设的相应限差δ₂停止迭代，得到最终的全局参数的最优估值x₂ ^lsq。具体实施时，本领域技术人员可自行根据精度需要设置限差δ₂，可将流程设计为判断是否满足着陆器坐标收敛条件，是则完成着陆器坐标(全局参数)的最后一次迭代，获取着陆器坐标(全局参数)精密坐标，否则返回进行全局参数迭代。

第八，解算出最新的Δx₂ ^lsq后，按照从下到上的次序回代到由以上步骤形成的最新的多弧段融合的法方程中，形式如式(39)所示，依次解算出第一次的局部参数的改正值按照式(20)加到局部参数初始近似值x₁ ^ref(n),x₁ ^ref(n-1),…,x₁ ^ref(2),x₁ ^ref(1)，即可得各个弧段局部参数对应的第一次的最优估值x₁ ^lsq(n),x₁ ^lsq(n-1),…,x₁ ^lsq(2),x₁ ^lsq(1)。

第九，将第一次得到的局部参数的最优估值x₁ ^lsq(i)重新作为初始近似值x₁ ^ref(i)，按照①～⑥的步骤更新法方程中的矩阵形成新的多弧段融合的法方程。其中，i＝N,N-1,…1，表示第i个弧段。

第十，根据新的法方程解算出新的局部参数向量的改正值按照式(20)加到各个弧段局部参数的初始近似值x₁ ^ref(i)得到新的最优估值x₁ ^lsq(i)，其中，i＝N,N-1,…1，表示第i个弧段。

第十一，重复第九和第十步，直至当次迭代和上次迭代得到的之差的绝对值小于预设的相应限差δ₁停止迭代，得到最终局部参数的最优估值x₁ ^lsq(i)，其中，i＝N,N-1,…1，表示第i个弧段。具体实施时，本领域技术人员可自行根据精度需要设置限差δ₁，可将流程设计为判断是否满足局部参数收敛条件，是则完成中继星初轨、力模型参数(局部参数)的最后一次迭代，获取中继星初轨、力模型参数估值，否则返回进行局部参数迭代。

至此，根据以上流程，得到最终稳定的中继星精密轨道和远月面着陆器的精密坐标。具体实施时，可采用计算机软件技术实现自动流程运行。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。

Claims (4)

1.一种基于四程中继跟踪模式的远月面着陆器精密定位方法，其特征在于：建立月球着陆器-轨道器的四程中继跟踪测量模式的观测模型，基于观测模型实现远月面着陆器精密定位，

所述四程中继跟踪测量模式，将月球轨道器作为中继星，使用月球着陆器和月球轨道器之间的链路进行测量，设地球跟踪站T_i在时刻i发送一个上行的信号给月球轨道器，月球轨道器位置为S_j，时刻为j；经过转发，信号传送给月面着陆器，相应位置为L_k，时刻为k，L_k位于月球背面；月面着陆器再次转发给在轨的月球轨道器，此时轨道器已经从S_j的位置运动至S_m，时刻为m；之后该信号下行至地球跟踪站T_n，时刻为n；设时刻n、m、k、j、i分别相应的参与者记为T_n,S_m,L_k,S_j,T_i，则有以下4个几何距离，

R₁＝|X(S_m)-X(T_n)|

R₂＝|X(S_m)-X(L_k)|

R₃＝|X(S_j)-X(L_k)|

R₄＝|X(S_j)-X(T_i)|

其中，X(j)为参与者j在太阳系质心坐标系下下的位置矢量；

四程中继测距的观测模型建立如下，

R＝(R₁+c·RLT_nm)+R₂+R₃+(R₄+c·RLT_ij)+

c·[TDB(i)-UTC(i)]-c·[TDB(n)-UTC(n)]

＝c·[UTC(n)-UTC(i)]

其中，RLT_nm为由参与者T_n发射信号到S_m的广义相对论时延，RLT_ij为由参与者T_i发射信号到S_j的广义相对论时延，UTC(n)和UTC(i)分别为时刻n和时刻i对应的协调世界时时标，R为四程中继测距值，TDB(n)和TDB(i)分别为时刻n和时刻i对应的太阳系质心力学时，c为光速；

设在一个多普勒积分周期内，起始时刻T_s和终止时刻T_e各自对应的四程中继测距值为R_s、R_e，建立四程中继测速的观测模型如下，

$RR=\frac{(R_e-R_s)}{T_c}$

其中，T_c为一个多普勒积分周期，RR为四程中继测速值。

2.根据权利要求1所述基于四程中继跟踪模式的远月面着陆器精密定位方法，其特征在于：根据四程中继跟踪测量模式，将月球轨道器的轨道和远月面着陆器的位置同时视为待估参数，采用精密定轨的方式进行解算，最后得到两者的最佳估值。

3.根据权利要求1所述基于四程中继跟踪模式的远月面着陆器精密定位方法，其特征在于：采用精密定轨的方式进行解算时，

计算四程中继测距值R对远月面着陆器坐标X(L_k)的偏导数H_R公式如下，

$H_R=\frac{X\left(L_k\right)-X\left(S_m\right)}{|X\left(L_k\right)-X\left(S_m\right)|}+\frac{X\left(L_k\right)-X\left(S_j\right)}{|X\left(L_k\right)-X\left(S_j\right)|}$

按照上式计算积分起始时刻T_s和积分终止时刻T_e分别对应的四程中继测距值R对远月面着陆器坐标X(L_k)的偏导数H_s和H_e，则四程中继测速值RR对远月面着陆器坐标X(L_k)的偏导数H_RR为

将结果作为观测时刻的观测值对状态向量的偏导数。

4.根据权利要求2或3所述基于四程中继跟踪模式的远月面着陆器精密定位方法，其特征在于：将探测器轨道跟踪数据文件分为多个观测弧段，将初始参考轨道x₀ ^ref＝[r₀ v₀ p₀q₀]^T中的中继星的初始位置r₀、初始速度v₀以及力模型参数p₀作为局部参数处理，记作x₁，维数为m₁；参数q₀是远月面着陆器位置坐标，将其作为全局参数处理，记作x₂，维数为3；采用精密定轨的方式进行解算时，融合多个观测弧段相应法方程，将非线性的方程线性化，使用迭代法求解。