一种重力卫星星间指向的在轨定标方法及系统与流程

1.本发明涉及地球物理测量技术领域，特别是涉及一种重力卫星星间指向的在轨定标方法及系统。


背景技术：

2.重力卫星的星间指向标定是为检核双星之间的对准情况，定义为卫星质心到天线相位中心矢量相对于视线连线方向的夹角，即
[0003][0004]
其中，r
sf
→
irf
为卫星参考坐标系到惯性系的转动矩阵，图1为卫星参考系(sf)，图2为天线相位中心参考坐标系(kf)和视线参考坐标系(los)。kf：其坐标原点位于卫星质心，x轴方向定义为从卫星质心指向天线相位中心，z轴方向定义为沿卫星飞行方向垂直向下，y轴方向通过右手定则确定。
[0005]
los：其坐标原点位于卫星质心，x轴方向为两个卫星质心视线方向的连线，z轴方向指向地心，y轴方向根据右手定则确定。
[0006]
视线参考坐标系定义为：
[0007][0008][0009][0010]
其中，上面式子中j＝1,2，表示卫星j，分别是卫星1和卫星2质心位置矢量，表示卫星j(1或2)在los坐标系下的x分量，表示卫星j(1或2)在los坐标系下的y分量，表示卫星j(1或2)在los坐标系下的z分量。
[0011]
天线相位中心参考坐标系定义为：
[0012][0013][0014][0015]
表示卫星j(1或者2)在kf坐标系下的x分量，表示卫星j(1或者2在kf坐标系下的y分量，表示卫星j(1或者2)在kf坐标系下的z分量，r
sf
→
irf
表示卫星参考系到惯
性系下的转换矩阵，为卫星质心到天线相位中心矢量。卫星参考系坐标轴与天线相位中心参考系一致。
[0016][0017]
星间指向要求角向偏差较小，其原因有两个：一是指向抖动影响星间测距观测量；二是由于微波信号在卫星表面进行多次反射形成的多径效应。
[0018]
现有的星间指向定标方法大多采用卫星稳定运行时的事后精密轨道及精密姿态数据对卫星的双星指向进行评价，此方法无法反映真实的星间指向情况，造成定标结果不准确。


技术实现要素：

[0019]
本发明的目的是提供一种重力卫星星间指向的在轨定标方法及系统，可对卫星实际在轨的双星指向变化进行标定，提高在轨定标结果的准确度。
[0020]
为实现上述目的，本发明提供了如下方案：
[0021]
一种重力卫星星间指向的在轨定标方法，包括：
[0022]
获取多个时刻下的观测数据；所述观测数据包括卫星的角度、精密定轨数据和双星的kbr测距数据；所述卫星的角度根据机动方案确定；所述角度为横滚角或者俯仰角；
[0023]
对于任意一个时刻，根据所述时刻下的kbr测距数据得到所述时刻下的双星中各卫星的质心到天线相位中心矢量的模值；
[0024]
根据所述时刻下的观测数据构建星间测距理论模型；
[0025]
根据所述星间测距理论模型计算所述时刻下的双星的星间测距的残差；
[0026]
将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的角度变化量；所述角度变化量为横滚角变化量或者俯仰角变化量。
[0027]
一种重力卫星星间指向的在轨定标系统，包括：
[0028]
获取模块，用于获取多个时刻下的观测数据；所述观测数据包括卫星的角度、精密定轨数据和双星的kbr测距数据；所述卫星的角度根据机动方案确定；所述角度为横滚角或者俯仰角；
[0029]
双星模值计算模块，用于对于任意一个时刻，根据所述时刻下的kbr测距数据得到所述时刻下的双星中各卫星的质心到天线相位中心矢量的模值；
[0030]
星间测距计算模块，用于根据所述时刻下的观测数据构建星间测距理论模型；
[0031]
残差计算模块，用于根据所述星间测距理论模型计算所述时刻下的双星的星间测距的残差；
[0032]
变化量计算模块，用于将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的角度变化量；所述角度变化量为横滚角变化量或者俯仰角变化量。
[0033]
根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明根据kbr测距数据得到双星中各卫星的质心到天线相位中心矢量的模值；根据所述时刻下的观测数据构
建星间测距理论模型；计算双星的星间测距的残差；将双星间的星间测距、残差、双星的kbr测距数据和双星中各卫星的质心到天线相位中心矢量的模值代入观测方程进行求解得到双星的角度变化量；根据双星中各卫星的质心到天线相位中心矢量的模值和角度变化量完成重力卫星星间指向的在轨定标，可对卫星实际在轨的双星指向变化进行标定，提高定标结果的准确度。
附图说明
[0034]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0035]
图1为sf示意图；
[0036]
图2为kf和los的示意图；
[0037]
图3为本发明实施例提供的一种重力卫星星间指向的在轨定标方法的流程图。
具体实施方式
[0038]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0039]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0040]
本发明所采用的技术方案是通过卫星指向机动周期搜索方法实现卫星双星指向定标，原则上，理想的星间指向使得两颗卫星质心和两个天线相位中心在两颗卫星质心连线上，记从kf转到los的转换矩阵为r
kf
→
los
，其可定义为：r
kf
→
los
＝r
x
(ψ)ry(θ)rz(φ)，其中，绕x，y和z轴的旋转矩阵，分别用rx(ψ)，ry(θ)，rz(φ)表示。r
kf
→
los
和r
sf
→
los
一样。根据和得到星间指向计算为：γ＝arccos(cosθcosφ)。
[0041]
设理想情况下，卫星坐标系到los坐标系预期的转换矩阵为r
sf,d
→
los
，真实情况下，卫星坐标系到los坐标系估计的转换矩阵为r
sf,e
→
los
，可写为r
sf,e
→
los
＝r
sf,d
→
los
·rsf,e
→
sf,d
，那么星间指向标定的内容是对转换矩阵r
sf,e
→
sf,d
进行标定。
[0042]
设r
sfj,e
→
sfj,d
＝r
x
(δψj)ry(δθj)rz(δφj),j＝1,2,则欧拉角δψj，δθj和δφj即为需要标定的量。r
sfj,d
→
los,11
变量的下标11代表矩阵1行1列；r
sfj,d
→
losfj,12
变量的下标12代表矩阵1行2列，r
sfj,d
→
losfj,13
变量的下标13代表矩阵1行3列。
[0043]
星间指向此时可以计算为cosγj＝r
sfj,d
→
los,11
cosδθjcosδφj+r
sfj,d
→
losfj,12
(-cosδψjsinδφj+sinδψjsinδθjcosδφj)+r
sfj,d
→
losfj,13
(sinδψjsinδφj+cosδψjsinδθjcosδφj)
[0044]
当欧拉角δψj，δθj和δφj为小量时，且r
sfj,d
→
losfj,11
＝1，r
sfj,d
→
losfj,12
＝0，rsfj,d
→
losfj,13
＝0，因此有
[0045]
下面研究利用机动开展星间指向标定：kbr数据给出的是两个相位中心的距离，可建模为两个相位中心的距离即双星距离
[0046]
其中，和分别为卫星1和卫星2天线相位中心在惯性系中的位置矢量，r
br
为系统差，r
nr
为测距噪声。天线相位中心在惯性系中的位置矢量可以表示为：
[0047]
和这里和分别为卫星1和卫星2质心在惯性系中的位置矢量，为卫星j上kbr天线相位中心相对于卫星j质心的位置矢量在卫星j本体坐标系中的表示，r
sfj
→
irf
(j＝1,2)为从卫星j坐标系转换到惯性系中的转换矩阵。记换到惯性系中的转换矩阵。记在计算l
12
时，可将同时转换到卫星1的los坐标系下，记为losf1中，r
irf
→
losf1
表示卫星惯性坐标系到卫星1los坐标系的转换矩阵；是卫星2天线相位中心到质心的距离；是卫星1天线相位中心到质心的距离，离，离，记记精确到1微米量级，可得到其中e为3维单位矩阵，其中e为3维单位矩阵，其中e为3维单位矩阵，是星间距离矢量的误差因，r
12
是双星间的质心距离；是卫星1在los坐标系下的x分量；δ是测量误差。因在losf1坐标系中，故有：
[0048][0049][0050][0051][0051]
指卫星j在惯性系下的速度矢量；卫星j在惯性系下的位置矢量。
[0052]
设计机动方案，改变卫星姿态时，可以得到：
[0053]rsfj
→
losf1
＝r
sfj,d
→
losf1rsfj,e
→
sfj,drsfj,t
→
sfj,d
，上式中，r
sfj,t
→
sfj,d
为机动引起的卫星姿态改变，利用欧拉角ψj，θj和φj，定义为：
[0054]rsfj,t
→
sfj,d
＝r
x
(ψj)ry(θj)rz(φj),j＝1,2。
[0055]
理想情况下，r
sf1,d
→
losf1
＝e,此时，若定义指代关系且对r
sfj,e
→
sfj,d
线性近似得到：
[0056][0057]
得到，
[0058][0059]
令φ为卫星横滚角、θ为卫星俯仰角、ψ为卫星偏航角，φ1表示卫星a的横滚角，φ2表示卫星b的横滚角，θ1表示卫星a的俯仰角，θ2表示卫星b的俯仰角，ψ1表示卫星a的偏航角，ψ2表示卫星b的偏航角，根据grace任务经验，周期性机动方案的设计含以下四个基本机动方案：
[0060]
①
通过卫星上的姿态控制系统产生特定的力矩，使得卫星a在偏航方向上，按以下方式机动：φ
10
是卫星机动初始偏置角；
[0061]
是卫星机动偏值振幅。
[0062]
此种机动方案下相位观测对卫星a沿滚动轴和俯仰轴的天线中心分量比较敏感，此种方案记为ma，在方案ma下，其它欧拉角
[0063]
ψ1＝θ1＝ψ2＝θ2＝φ2＝0。
[0064]
②
与
①
类似，不同的是，此机动控制卫星a在俯仰方向上，按以下方式机动：
[0065]
此种机动方案下相位观测对卫星a沿滚动轴和偏航轴的天线中心分量比较敏感，此种方案记为mb，在方案mb下，其它欧拉角
[0066]
ψ1＝φ1＝ψ2＝θ2＝φ2＝0。
[0067]
③
与
①
类似，不同的是，此机动控制卫星b在偏航方向上，按以下方式机动：
[0068]
此种机动方案下相位观测对卫星b沿滚动轴和俯仰轴的天线中心分量比较敏感，此种方案记为mc，在方案mc下，其它欧拉角
[0069]
ψ1＝θ1＝φ1＝ψ2＝θ2＝0。
[0070]
④
与
③
类似，不同的是，此机动控制卫星b在俯仰方向上，按以下方式机动：
[0071]
此种机动方案下相位观测对卫星b沿滚动轴和偏航轴的天线中心分量比较敏感，此种方案记为md，在方案md下，其它欧拉角
[0072]
ψ1＝θ1＝φ1＝ψ2＝φ2＝0。
[0073]
在机动方案ma中，有
[0074][0075][0076][0077][0078]
计算偏导数，有
[0079][0080][0081][0082][0083]
在机动方案mb中，有
[0084][0085]
[0086][0087][0088]
偏导数为
[0089][0090][0091][0092][0093]
在机动方案mc中，有
[0094][0095][0096][0097][0098]
偏导数为
[0099][0100]
[0101][0102][0103]
在机动方案md中，有
[0104][0105][0106][0107][0108]
偏导数为
[0109][0110][0111][0112][0113]
从以上偏导数计算可知，各方案可对不同参数进行估计，如表1所示。
[0114]
表1不同参数不同方案中是否可估计表(可：√，不可：
×
)
[0115]
[0116][0117]
综上所述，本发明实施例提供了一种重力卫星星间指向的在轨定标方法，如图3所示，包括：
[0118]
步骤101：获取多个时刻下的观测数据；所述观测数据包括卫星的角度、精密定轨数据和双星的kbr测距数据；所述卫星的角度根据机动方案确定；所述角度为横滚角或者俯仰角。
[0119]
步骤102：对于任意一个时刻，根据所述时刻下的kbr测距数据得到所述时刻下的双星中各卫星的质心到天线相位中心矢量的模值。
[0120]
步骤103：根据所述时刻下的观测数据构建星间测距理论模型。
[0121]
步骤104：根据所述星间测距理论模型计算所述时刻下的双星的星间测距的残差，在卫星进行ma机动后可以得到和利用就可以得到星间距离残差。
[0122]
星间测距理论模型的表达形式是根据公式星间测距理论模型的表达形式是根据公式推到得到的，利用泰勒公式根据星间测距理论模型就可以得到此方法为公知的。
[0123]
步骤105：将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的角度变化量；所述角度变化量为横滚角变化量或者俯仰角变化量，可以根据所述双星中各卫星的质心到天线相位中心矢量的模值和角度变化量完成重力卫星星间指向的在轨定标。
[0124]
在实际应用中，所述获取多个时刻下的观测数据，具体包括：
[0125]
当所述机动方案为使双星中的第一个卫星在偏航方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第一个卫星的横滚角。
[0126]
当所述机动方案为使双星中的第一个卫星在俯仰方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第一个卫星的俯仰角。
[0127]
当所述机动方案为使双星中的第二个卫星在偏航方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第二个卫星的横滚角。
[0128]
当所述机动方案为使双星中的第二个卫星在俯仰方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第二个卫星的俯仰角。
[0129]
在实际应用中，所述根据所述时刻下的kbr测距数据得到所述时刻下的双星中各卫星的质心到天线相位中心矢量的模值，具体包括：
[0130]
计算所述双星中第一个卫星的质心到天线相位中心矢量的模值的参考值d
pc10
与所述第一个卫星的质心到天线相位中心矢量的模值的变化量δd
pc1
之和得到第一个卫星的质心到天线相位中心矢量的模值d
pc1
。
[0131]
计算所述双星中第二个卫星的质心到天线相位中心矢量的模值的参考值d
pc20
与
所述第二个卫星的质心到天线相位中心矢量的模值的变化量δd
pc2
之和得到第二个卫星的质心到天线相位中心矢量的模值d
pc2
。
[0132]
在实际应用中，所述星间测距理论模型具体为：
[0133]
当所述机动方案为使双星中的第一个卫星在偏航方向上运行时，则所述星间测距理论模型为其中，l
12
为双星间的星间测距，r
12
为精密定轨数据，d
pc20
为所述第二个卫星的质心到天线相位中心矢量的模值的参考值，d
pc10
为所述第一个卫星的质心到天线相位中心矢量的模值的参考值，φ1为所述第一个卫星的横滚角。
[0134]
当所述机动方案为使双星中的第一个卫星在俯仰方向上运行时，则所述星间测距理论模型为其中，θ1为第一个卫星的俯仰角。
[0135]
当所述机动方案为使双星中的第二个卫星在偏航方向上运行时，则所述星间测距理论模型为其中，φ2为第二个卫星的横滚角。
[0136]
当所述机动方案为使双星中的第二个卫星在俯仰方向上运行时，则所述星间测距理论模型为其中θ2为第二个卫星的俯仰角。
[0137]
在实际应用中，所述将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的角度变化量，具体包括：
[0138]
当所述机动方案为使双星中的第一个卫星在偏航方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的横滚角变化量。
[0139]
当所述机动方案为使双星中的第一个卫星在俯仰方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的俯仰角变化量。
[0140]
当所述机动方案为使双星中的第二个卫星在偏航方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的横滚角变化量。
[0141]
当所述机动方案为使双星中的第二个卫星在俯仰方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的俯
仰角变化量。
[0142]
针对上述方法，本发明实施例提供了一种重力卫星星间指向的在轨定标系统，包括：
[0143]
获取模块，用于获取多个时刻下的观测数据；所述观测数据包括卫星的角度、精密定轨数据和双星的kbr测距数据；所述卫星的角度根据机动方案确定；所述角度为横滚角或者俯仰角。
[0144]
双星模值计算模块，用于对于任意一个时刻，根据所述时刻下的kbr测距数据得到所述时刻下的双星中各卫星的质心到天线相位中心矢量的模值。
[0145]
星间测距计算模块，用于根据所述时刻下的观测数据构建星间测距理论模型。
[0146]
残差计算模块，用于根据所述星间测距理论模型计算所述时刻下的双星的星间测距的残差。
[0147]
变化量计算模块，用于将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的角度变化量；所述角度变化量为横滚角变化量或者俯仰角变化量。
[0148]
作为一种可选的实施方式，所述获取模块，具体包括：
[0149]
ma获取单元，用于当所述机动方案为使双星中的第一个卫星在偏航方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第一个卫星的横滚角。
[0150]
mb获取单元，用于当所述机动方案为使双星中的第一个卫星在俯仰方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第一个卫星的俯仰角。
[0151]
mc获取单元，用于当所述机动方案为使双星中的第二个卫星在偏航方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第二个卫星的横滚角。
[0152]
md获取单元，用于当所述机动方案为使双星中的第二个卫星在俯仰方向上运行时，则获取多个时刻下的所述精密定轨数据、所述双星的kbr测距数据和所述第二个卫星的俯仰角。
[0153]
作为一种可选的实施方式，所述双星模值计算模块，具体包括：
[0154]
第一个卫星的模值计算单元，用于计算所述双星中第一个卫星的质心到天线相位中心矢量的模值的参考值与所述第一个卫星的质心到天线相位中心矢量的模值的变化量之和得到第一个卫星的质心到天线相位中心矢量的模值。
[0155]
第二个卫星的模值计算单元，用于计算所述双星中第二个卫星的质心到天线相位中心矢量的模值的参考值与所述第二个卫星的质心到天线相位中心矢量的模值的变化量之和得到第二个卫星的质心到天线相位中心矢量的模值。
[0156]
作为一种可选的实施方式，所述星间测距理论模型具体为：
[0157]
当所述机动方案为使双星中的第一个卫星在偏航方向上运行时，则所述星间测距理论模型为其中，l
12
为双星间
的星间测距，r
12
为精密定轨数据，d
pc20
为所述第二个卫星的质心到天线相位中心矢量的模值的参考值，d
pc10
为所述第一个卫星的质心到天线相位中心矢量的模值的参考值，φ1为所述第一个卫星的横滚角。
[0158]
当所述机动方案为使双星中的第一个卫星在俯仰方向上运行时，则所述星间测距理论模型为其中，θ1为第一个卫星的俯仰角。
[0159]
当所述机动方案为使双星中的第二个卫星在偏航方向上运行时，则所述星间测距理论模型为其中，φ2为第二个卫星的横滚角。
[0160]
当所述机动方案为使双星中的第二个卫星在俯仰方向上运行时，则所述星间测距理论模型为其中θ2为第二个卫星的俯仰角。
[0161]
作为一种可选的实施方式，所述变化量计算模块，具体包括：
[0162]
ma变化量计算单元，用于当所述机动方案为使双星中的第一个卫星在偏航方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的横滚角变化量。
[0163]
mb变化量计算单元，用于当所述机动方案为使双星中的第一个卫星在俯仰方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的俯仰角变化量。
[0164]
mc变化量计算单元，用于当所述机动方案为使双星中的第二个卫星在偏航方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的横滚角变化量。
[0165]
md变化量计算单元，用于当所述机动方案为使双星中的第二个卫星在俯仰方向上运行时，则将所有时刻下的所述双星间的星间测距、所述残差、所述双星的kbr测距数据和所述双星中各卫星的质心到天线相位中心矢量的模值代入观测方程，采用带权最小二乘估计进行求解得到双星的俯仰角变化量。
[0166]
本发明实施例提供了一种更加具体的重力卫星星间指向的在轨定标方法。
[0167]
第一步，进行ma的机动方案指向计算；
[0168]
机动方案ma可估计参数d
pc1
，δφ1，d
pc2
，δφ2。记d
pc1
＝d
pc10
+δd
pc1
，d
pc2
＝d
pc20
+δd
pc2
，其中d
pc10
和d
pc20
分别为d
pc1
和d
pc2
的参考值。则估计算法为：
[0169]
s1.读入kbr测距数据，获取d
pc20
、d
pc10
。kbr测距数据文件中包含d
pc20
、d
pc10
。
[0170]
s2.读入精密定轨数据，获取r
12
，构造星间测距理论模型
[0171][0172]
其中
[0173][0174]
s3.计算星间测距l
12
的残差y(t)。
[0175]
s4.确定估计模型，建立观测方程。利用4阶多项式拟合精密定轨产生的残差。建立观测方程
[0176][0177]
其中
[0178]
p4(t)＝a0+a1t+a2t2+a3t3+a4t4[0179]
令待估参数矢量为
[0180][0181]
观测方程可以写为
[0182][0183]
其中
[0184]
h(t)＝(h
poly
(t)h
vip
(t))
[0185]hpoly
(t)＝(1,t,t2,t3,t4[0186][0187]
上式中偏导数和
[0188]
将所有观测值集中，假定有n次观测量，定义
[0189][0190]
观测方程可写为
[0191][0192]
s5.利用带权最小二乘估计对上式进行求解计算，可得参数。
[0193]
[0194]
估计协方差矩阵为(h
t
wh)-1
。
[0195]
此处是最小二乘估计最终的参数，可实现双星横滚角变化量δφ1，δφ2的估计。
[0196]
第二步，进行mb的机动方案指向计算；
[0197]
机动方案mb可估计参数d
pc1
，δθ1，d
pc2
，δθ2。估计算法为：
[0198]
s1.读入kbr测距数据，获取d
pc20
、d
pc10
。
[0199]
s2.读入精密定轨数据，获取r
12
，构造星间测距理论模型
[0200][0201]
其中
[0202][0203]
s3.计算星间测距l
12
的残差y(t)。
[0204]
s4.确定估计模型，建立观测方程。利用4阶多项式拟合精密定轨产生的残差。建立观测方程
[0205][0206]
其中
[0207]
p4(t)＝a0+a1t+a2t2+a3t3+a4t4[0208]
令待估参数矢量为
[0209][0210]
观测方程可以写为
[0211][0212]
其中
[0213]
h(t)＝(h
poly
(t)h
vip
(t))
[0214]hpoly
(t)＝(1,t,t2,t3,t4)
[0215][0216]
上式中偏导数和将所有观测值集中，假定有n次观测量，定义
[0217]
[0218]
观测方程可写为
[0219][0220]
s5.利用带权最小二乘估计对上式进行求解计算，可得参数。
[0221][0222]
估计协方差矩阵为(h
t
wh)-1
。
[0223]
此处是最小二乘估计最终的参数，可实现双星俯仰角变化量δθ1，δθ2的估计。
[0224]
第三步，进行mc的机动方案指向计算；
[0225]
机动方案mc可估计参数d
pc1
，δφ1，d
pc2
，δφ2。估计算法为：
[0226]
s1.读入kbr测距数据，获取d
pc20
、d
pc10
。
[0227]
s2.读入精密定轨数据，获取r
12
，构造星间测距理论模型
[0228][0229]
其中
[0230][0231]
s3.计算星间测距l
12
的残差y(t)。
[0232]
s4.确定估计模型，建立观测方程。利用4阶多项式拟合精密定轨产生的残差。建立观测方程
[0233][0234]
其中
[0235]
p4(t)＝a0+a1t+a2t2+a3t3+a4t4[0236]
令待估参数矢量为
[0237][0238]
观测方程可以写为
[0239][0240]
其中
[0241]
h(t)＝(h
poly
(t)h
vip
(t))
[0242]hpoly
(t)＝(1,t,t2,t3,t4[0243]
[0244]
上式中偏导数和将所有观测值集中，假定有n次观测量，定义
[0245][0246]
观测方程可写为
[0247][0248]
s5.利用带权最小二乘估计对上式进行求解计算，可得参数。
[0249][0250]
估计协方差矩阵为(h
t
wh)-1
。
[0251]
此处是最小二乘估计最终的参数，可实现双星横滚角变化量δφ1，δφ2的估计。
[0252]
第四步，进行md的机动方案指向计算；
[0253]
机动方案md可估计参数d
pc1
，δθ1，d
pc2
，δθ2。估计算法为：
[0254]
s1.读入kbr测距数据，获取d
pc20
、d
pc10
。
[0255]
s2.读入精密定轨数据，获取r
12
，构造星间测距理论模型
[0256][0257]
其中
[0258][0259]
s3.计算星间测距l
12
的残差y(t)。
[0260]
s4.确定估计模型，建立观测方程。利用4阶多项式拟合精密定轨产生的残差。建立观测方程
[0261][0262]
其中
[0263]
p4(t)＝a0+a1t+a2t2+a3t3+a4t4[0264]
令待估参数矢量为
[0265][0266]
观测方程可以写为
[0267]
[0268]
其中
[0269]
h(t)＝(h
poly
(t)h
vip
(t))
[0270]hpoly
(t)＝(1,t,t2,t3,t4)
[0271][0272]
上式中偏导数和
[0273]
将所有观测值集中，假定有n次观测量，定义
[0274][0275]
观测方程可写为
[0276][0277]
s5.利用带权最小二乘估计对上式进行求解计算，可得参数。
[0278][0279]
估计协方差矩阵为(h
t
wh)-1
。
[0280]
此处是最小二乘估计最终的参数，可实现双星横滚角变化量δθ1，δθ2的估计。
[0281]
本实施例将上述方法应用于grace follow-on任务提供的实测数据作为测量误差的实测数据进行仿真，本仿真采用的数据包含2018年6月1日得到的1b级精密定轨数据gni1b、1b级相位中心间的测距数据kbr1b、及天线相位中心1b级标定数据vkb1b数据。
[0282]
在kbr1b数据当中，含有偏测距、测距变化率、测距加速度、飞行时间修正和天线相位中心偏差修正等，2018年6月1日当天的kbr1b数据，时间间隔为5秒。grace follow-on任务中2018年6月1获得的精密定轨数据gni1b给出的星间距离和星间距离变率，gni1b数据当中，含有地心惯性系下卫星位置和速度，时间间隔1秒。
[0283]
(1)时长1000s情况下三种机动周期仿真
[0284]
首先开展在机动时长为1000s时三种机动周期情况下仿真。卫星a偏航(yaw)方向偏差角真值为1mrad。进行了25次重复实验，计算得到的卫星a偏航(yaw)方向的估计值与理论值的绝对误差，从实验结果可以得，13次251s的机动方案的绝对误差百分比要小于其它两种情形，8次250s机动周期绝对值误差小于其它两种情形，4次249s的机动方案的绝对误差百分比要小于其它两种情形。因此，不同机动周期会对结果产生影响，相较而言，249s机动周期最差，251s机动周期最好。因此，后续，将针对250s和251s机动周期开展进一步仿真，通过仿真研究指向标定能力。
[0285]
表2仿真机动方案ma机动参数设置
[0286] 机动偏值(
°
)机动幅值(
°
)机动周期(s)机动时长(s)情形1212491000情形2212501000
情形3212511000
[0287]
(2)1000s机动时长下249s，250s和251s机动周期情况下解算不同指向偏差角的能力，并且同时解算相位中心矢量幅值。
[0288]
考虑到数据质量限制，指向偏差标定具有一定的能力。下面针对指向偏差角真值分别为1mrad，0.5mrad，0.1mrad，0.05mrad和0.01mrad时开展解算质量评估。针对每一个指向偏差真值，考虑机动时长为1000s，机动周期分别为249s，250s和251s三种情形下不同大小指向偏差角解算能力。进行了6次仿真结果。从这些仿真结果可得，真值为1mrad时，三种机动周期情形下，解算误差均小于10％。在真值为0.5mrad时，机动周期为249s，估计值与真值偏差在11.7％与29.2％之间。机动周期为250s，估计值与真值偏差在1.7％与18.7％之间，且有四种情况下小于10％。机动周期为251s，估计值与真值偏差在0.3％与12.7％之间，且有四种情况下小于10％。在真值为0.1mrad时，估计结果较差，其中机动周期为251s时，结果最好。机动周期为251s时，估计值与真值偏差在6.2％与60.2％之间，且有五种情况下大于10％。当偏差角真值为0.05mrad和0.01mrad时，解算精度更差。综上，在真值为0.5mrad时，利用机动周期为251s方案，勉强可实现估算值与偏差值的绝对误差接近于10％，说明绝对标定下指向角偏值不能过小。总的来说，在指向偏差真值小于1mrad时候解算精度较差，其中251s机动周期方案估计精度最好。
[0289]
(3)1000s机动时长下249s，250s和251s机动周期情况下解算不同指向偏差角的能力，不解算相位中心矢量幅值。
[0290]
下面针对指向偏差角真值分别为1mrad，0.5mrad,0.1mrad,0.05mrad和0.01mrad时开展解算质量评估。针对每一个指向偏差真值，考虑机动时长为1000s，机动周期分别为249s，250s和251s三种情形下不同大小指向偏差角解算能力。进行了6次仿真结果。从这些仿真结果可以看到，真值为1mrad时，三种机动周期情形下，解算误差均小于7％。在真值为0.5mrad时，机动周期为249s和250s情况下，估计值与真值偏差在7％以内。机动周期为251s，估计值与真值偏差在10％左右。在真值为0.1mrad，0.05mrad和0.01mrad时，解算精度变差，超过了10％。综上，在真值为0.5mrad时，利用机动周期为249s和250s方案，勉强可实现估算值与偏差值的绝对误差接近于10％，说明绝对标定下指向角偏值真值不能过小。总的来说，在指向偏差真值小于0.5mrad时候解算精度较差，其中249s机动周期方案估计精度最好，250s机动周期方案估计精度也可小于7％。
[0291]
从上述实验结果可得，当在同时解算相位中心幅值与指向偏差角时，指向偏差角小于1mrad，解算结果较差。如果为了能够解算更小指向偏差角，例如0.5mrad，采用3000s机动时长为宜。但是，总的来说，解算1mrad指向偏差角时能够保证解算精度。
[0292]
本发明的有益效果是通过卫星机动实现指向周期摆动，获得星间指向观测数据，估计出卫星星间指向，可以避免传统方法无法真实准确反应双星指向情况，从而获得卫星精确的双星指向控制情况参数，提高在轨定标结果的准确度。
[0293]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0294]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说
明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。