本发明涉及深空测控领域，具体涉及一种计算双探测器在地外天体表面相对位置的方法。
背景技术：
在月球及深空地外天体着陆巡视探测任务中，“着陆器+巡视器”组合探测模式十分常见。嫦娥三号任务就采用该模式实现了我国对月球的初次着陆勘测，后续发射的嫦娥四号任务也将类似地对月球背面开展探测，而将于2020年发射的我国首次火星探测任务也包含了着陆器和火星车巡视器。在该类任务中，获取着陆器和巡视器精确的相对位置对于工程实施和科学探测至关重要。
在国外月球和火星着陆探测以及我国嫦娥三号任务中，通常采用“视觉+惯导”的组合方法对两个探测器进行相对定位，但视觉方法作用距离有限，当两器之间的视线存在遮挡或者距离太远不可见时，视觉定位方法就无法开展；而单纯的惯导方法则会随时间积累产生巨大的误差(参见Liu，Z.Q.，Di，K.C.，Peng，M.，et al.，2015.High precision landing site mapping and rover localization for Chang’e-3 mission.Sci.China Phys.Mech.Astron.58，019601.)。另一种方法是通过地面无线电跟踪测量来计算两者的相对位置，但通常需要测距、测速和干涉测量三种数据综合才能进行定位，测距、测速数据用来约束视向位置关系，而干涉测量数据用来约束横向关系。
技术实现要素：
本发明要解决的技术问题是：提供一种不需要测距和测速数据，不受视觉可见性约束，不存在累积误差的计算双探测器在地外天体表面相对位置的方法。
为了解决上述技术问题，本发明的技术方案如下：
一种计算双探测器在地外天体表面相对位置的方法，包括如下步骤：
步骤一、计算着陆器在地心J2000坐标系中的赤经αL、赤纬δL和着陆器至地心的距离dL；
步骤二、计算巡视器在地心J2000坐标系中的赤经αR和赤纬δR，并给巡视器至地心的距离dR赋初值dL；
步骤三、计算巡视器在地外天体中心固连坐标系中的高程值偏差ΔH；
3.1根据巡视器在地心J2000坐标系中的赤经αR、赤纬δR，和巡视器至地心的距离dR，计算巡视器在地外天体表面的经度LR，纬度BR和高程值HR；
3.2根据步骤3.1得到的经度LR和纬度BR从地外天体数字高程模型中确定巡视器的数字高程值
3.3根据公式计算巡视器在地外天体中心固连坐标系中的高程值偏差ΔH；
步骤四、判断巡视器高程值偏差ΔH是否小于10-5m；
若ΔH＜10-5m，跳至步骤五；若ΔH≥10-5m，则按照公式dR_new＝dR+K·ΔH，K为调整系数，将巡视器在地心J2000坐标系中的距离值dR赋值为dR_new后，跳转至步骤三；
步骤五、计算巡视器相对着陆器的地外天体表面北东地位置。
进一步的，所述步骤五中，巡视器相对着陆器的地外天体表面北东地位置计算过程如下：
1)根据下式计算地外天体中心固连坐标系到着陆器地外天体表面北东地坐标系的转换矩阵M3；
式中，LL、BL分别为着陆器在地外天体表面的经度、纬度。
2)根据下式计算并输出巡视器相对着陆器的地外天体表面北东地位置(xR2L，yR2L，zR2L)：
式中，(xR，yR，zR)为巡视器在地外天体中心固连坐标系中的三维坐标，(xL，yL，zL)为着陆器在地外天体中心固连坐标系中的三维坐标，M3为地外天体中心固连坐标系到着陆器地外天体表面北东地坐标系的转换矩阵。
进一步的，所述步骤一还包括：
1)利用着陆器在地外天体表面的经度LL、纬度BL和高程HL，根据下式计算着陆器在地外天体中心固连坐标系中的三维坐标(xL，yL，zL)；
式中，N为地外天体的参考半径，LL、BL和HL分别为着陆器在地外天体表面的经度、纬度和高程；
2)从行星历表获得地外天体中心固连坐标系与地心J2000坐标系的转换矩阵M1及逆转换矩阵M2，和地外天体中心在地心J2000坐标系中的三维坐标(xO，yO，zO)；
3)根据下式计算着陆器在地心J2000坐标系中的三维坐标
式中，(xO，yO，zO)为地外天体中心在地心J2000坐标系中的三维坐标，M1为地外天体中心固连坐标系与地心J2000坐标系的转换矩阵，(xL，yL，zL)为着陆器在地外天体中心固连坐标系中的三维坐标。
4)根据下式计算着陆器在地心J2000坐标系中的赤经αL、赤纬δL和着陆器至地心的距离dL；
式中，为着陆器在地心J2000坐标系中的三维坐标。
进一步的，所述步骤三中，巡视器在地外天体表面的经度LR，纬度BR和高程值HR计算过程如下：
1)根据下式计算巡视器在地心J2000坐标系中的三维坐标
式中，dR为巡视器至地心的距离，αR、δR分别为巡视器在地心J2000坐标系中赤经和赤纬。
2)根据下式计算巡视器在地外天体中心固连坐标系中的三维坐标(xR，yR，zR)；
式中，为巡视器在地心J2000坐标系中的三维坐标，(xO，yO，zO)为地外天体中心在地心J2000坐标系中的三维坐标，M2为地外天体中心固连坐标系与地心J2000坐标系的逆转换矩阵。
3)根据下式计算巡视器在地外天体表面的经度LR，纬度BR和高程HR；
式中，(xR，yR，zR)为巡视器在地外天体中心固连坐标系中的三维坐标。
进一步的，所述的步骤二中，按照下式计算巡视器在地心J2000坐标系中的赤经αR和赤纬δR：
式中，Δα，Δδ为地面无线电干涉测量获得的某时刻巡视器相对着陆器的赤经差、赤纬差；αL、δL分别为着陆器在地心J2000坐标系中的赤经、赤纬。
本发明所达到的有益效果：
本发明通过循环迭代的方式快速计算出巡视器相对着陆器的地外天体表面北东地位置，收敛效果好，计算精度高，适用于太阳系内各个地外天体，不需要测距和测速数据，不受视觉可见性约束，不存在累积误差。
附图说明
图1为本发明技术方案的流程图。
具体实施方式
以下结合附图和嫦娥三号实测数据对本发明的具体实施方式作进一步详细说明。
本发明的原理：利用地面无线电干涉测量得到的两个探测器赤经差、赤纬差和地外天体数字高程模型，通过循环迭代修正巡视器高程值，计算两个探测器在地外天体表面的北东地相对位置。
假设两个探测器分别命名为着陆器和巡视器，其中着陆器在地外天体表面的经度LL、纬度BL和高程HL已经确定，巡视器在地外天体表面自主或半自主运动。由地面无线电干涉测量获得的某时刻巡视器相对着陆器的赤经差、赤纬差为(Δα，Δδ)。视地外天体为均匀球体，参考半径为N。图1示出了本发明实现的流程，具体如下：
步骤一、计算着陆器在地心J2000坐标系中的赤经αL、赤纬δL和着陆器至地心的距离dL；
1)利用着陆器在地外天体表面的经度LL、纬度BL和高程HL，根据式(1)计算着陆器在地外天体中心固连坐标系中的三维坐标(xL，yL，zL)；
式中，N为地外天体的参考半径，LL、BL和HL分别为着陆器在地外天体表面的经度、纬度和高程。
地外天体包括，但不限于月球、火星、金星、小行星。
2)从喷气推进实验室(JPL)行星历表获得地外天体中心固连坐标系与地心J2000坐标系的转换矩阵M1及逆转换矩阵M2，和地外天体中心在地心J2000坐标系中的三维坐标(xO，yO，zO)；
3)根据式(2)计算着陆器在地心J2000坐标系中的三维坐标
式中，(xO，yO，zO)为地外天体中心在地心J2000坐标系中的三维坐标，M1为地外天体中心固连坐标系与地心J2000坐标系的转换矩阵，(xL，yL，zL)为着陆器在地外天体中心固连坐标系中的三维坐标。
4)根据式(3)计算着陆器在地心J2000坐标系中的赤经αL、赤纬δL和着陆器至地心的距离dL；
式中，为着陆器在地心J2000坐标系中的三维坐标。
步骤二、计算巡视器在地心J2000坐标系中的赤经αR和赤纬δR，并给巡视器至地心的距离dR赋初值dL；
根据式(4)计算巡视器在地心J2000坐标系中赤经αR和赤纬δR，考虑到巡视器和着陆器相对距离远小于巡视器到地心的距离，给巡视器至地心的距离dR赋初值dL；
式中，Δα，Δδ为地面无线电干涉测量获得的某时刻巡视器相对着陆器的赤经差、赤纬差；
步骤三、计算巡视器在地外天体中心固连坐标系中的高程值偏差ΔH；
3.1计算巡视器在地外天体表面的的经度LR，纬度BR和高程值HR；
1)根据式(5)计算巡视器在地心J2000坐标系中的三维坐标
式中，dR为巡视器至地心的距离，αR、δR分别为巡视器在地心J2000坐标系中赤经和赤纬。
2)根据式(6)计算巡视器在地外天体中心固连坐标系中的三维坐标(xR，yR，zR)；
式中，为巡视器在地心J2000坐标系中的三维坐标，(xO，yO，zO)为地外天体中心在地心J2000坐标系中的三维坐标，M2为地外天体中心固连坐标系与地心J2000坐标系的逆转换矩阵。
3)根据式(7)计算巡视器在地外天体表面的经度LR，纬度BR和高程HR；
式中，(xR，yR，zR)为巡视器在地外天体中心固连坐标系中的三维坐标。
3.2根据步骤3.1中得到的经度LR和纬度BR从地外天体数字高程模型中确定巡视器的数字高程值
3.3根据式(8)计算巡视器在地外天体中心固连坐标系中的高程值偏差ΔH；
步骤四、判断巡视器高程值偏差ΔH是否小于10-5m；
若ΔH＜10-5m，认为巡视器的计算位置与其在地外天体数字高程模型中的位置完全吻合，跳至步骤五；如果ΔH≥10-5m，则按式(9)将巡视器在地心J2000坐标系中的距离值dR赋值为dR_new后，跳转至步骤三，开始循环迭代计算；
dR_new＝dR+K·ΔH (9)
式中，K为调整系数，可以根据迭代收敛情况适当选择。
步骤五、计算巡视器相对着陆器的地外天体表面北东地位置；
1)根据式(10)计算地外天体中心固连坐标系到着陆器地外天体表面北东地坐标系的转换矩阵M3；
式中，LL、BL分别为着陆器在地外天体表面的经度、纬度。
2)根据式(11)计算并输出巡视器相对着陆器的地外天体表面北东地位置(xR2L，yR2L，zR2L)；
式中，(xR，yR，zR)为巡视器在地外天体中心固连坐标系中的三维坐标，(xL，yL，zL)为着陆器在地外天体中心固连坐标系中的三维坐标，M3为地外天体中心固连坐标系到着陆器地外天体表面北东地坐标系的转换矩阵。
以下结合嫦娥三号实测数据对本发明的具体实施方式作进一步详细说明。
嫦娥三号任务是我国第一个月球表面着陆探测任务，包含着陆器和巡视器两部分。组合体着陆后，巡视器从着陆器上分离，开展月面半自主探测，着陆器则在原位开展定点观测。通过多种手段测得着陆器月面位置为44.1206°N，-19.5124°E，高程-2.632km(月球参考半径1737.4km)。四个地面测站(上海天马，云南昆明，北京密云和新疆乌鲁木齐)对两个探测器开展了无线电干涉测量，在2013年12月15日16:10测得巡视器相对着陆器的赤经、赤纬差为(-3.29，3.36)mas(毫角秒)。
Step 1：利用着陆器在月球表面的经度LL、纬度BL和高程HL计算得到着陆器在月球中心固连坐标系(平轴系)中的三维坐标为(1173.82707，-415.95989，1207.69508)km；
Step 2：从JPL行星历表获得该时刻月球中心固连坐标系与地心J2000坐标系的转换矩阵M1及逆转换矩阵M2，和月球中心在地心J2000坐标系中的三维坐标(xO，yO，zO)分别为：
(xO，yO，zO)＝(171364.65484，338150.65016，128316.59429)km；
Step 3：计算得到着陆器在地心J2000坐标系中的三维坐标为(170442.82374，336897.75265，129084.64069)km；
Step 4：计算得到着陆器在地心J2000坐标系中的赤经为63.1643050544°、赤纬为18.8751770102°，距离为399015.90992km；
Step 5：计算得到巡视器在地心J2000坐标系中赤经为63.1643041405°，赤纬为18.8751779435°，并给距离赋初值为399015.90992km；
Step 6：计算得到巡视器在地心J2000坐标系中的三维坐标为(170442.82837，336897.74848，129084.64701)km；
Step 7：计算得到巡视器在月球中心固连坐标系中的三维坐标为(1173.82584，-415.95538，1207.70263)km；
Step 8：计算得到巡视器在月球表面的经度为-19.51224°E，纬度为44.12084°N，高程为-2.62866km。由于着陆器和巡视器相隔很近，着陆器所在周边区域地形平坦，从月球数字高程模型中可以近似认为巡视器和着陆器处于同一水平面，再考虑两器信号发射天线的安装位置高度差为1.52m，可以得到巡视器的数字高程值为-2.63352km。
Step 9：计算得到巡视器高程值偏差为4.86m；
Step 10：巡视器高程值偏差大于10-5m，系数K设为0.0005，调整巡视器在地心J2000坐标系中的距离值dR为399015.91235km，跳转至Step 6；
经过32次循环迭代后，巡视器高程偏差小于10-5m，跳转至Step 11，此时dR为399015.91743km；
Step 11：计算得到月球中心固连坐标系到着陆器月球表面北东地坐标系的转换矩阵为
Step 12：计算并输出巡视器相对着陆器的月球表面北东地位置为(11.47，1.02，1.52)m。
根据视觉定位方法得到的两器相对位置为(10.87，1.11，1.52)m。
由于巡视器在月球表面运动，其相对着陆器的位置在不断变化，在2013年12月16日12:00、2013年12月20日17:00、2013年12月20日22:30、2013年12月21日19:10等观测时间，利用本发明方法与视觉定位方法得到的两器相对位置对比如下表。本发明方法与视觉定位方法的定位偏差均方根为(0.42，0.18，0.22)m，即位置偏差小于1m，充分验证本发明方法的有效性和高精度。
表1不同时刻得到的两器相对位置结果对比
通过上述实施例，完全有效地实现了本发明的目的。该领域的技术人员可以理解本发明包括但不限于附图和以上具体实施方式中描述的内容。虽然本发明就目前认为最为实用且优选的实施例进行说明，但应知道，本发明并不限于所公开的实施例，任何不偏离本发明的功能和结构原理的修改都将包括在权利要求书的范围中。