CN110203424B - 利用测速数据估计航天器自旋运动的方法和设备 - Google Patents

Abstract

本发明提出利用测速数据估计航天器自旋运动的方法和设备，适用于航天器姿态发生自旋故障，地面无法获得稳定遥测数据的情况下，进行航天器自旋运动的估计，从而辅助地面工程技术人员及时判断和处理姿态故障。本发明选定航天器的自旋角速度大小和方向为参数，建立测速数据观测模型。本发明通过航天器测速数据估计航天器自旋运动，适应性强，收敛迅速。该方案不依赖航天器遥测信息，可为航天器姿态异常情况下的姿态故障诊断和处置提供有效依据。

Description

利用测速数据估计航天器自旋运动的方法和设备

技术领域

本发明一般地涉及航天器测控技术领域。更具体地，本发明涉及利用测速数据估计航天器自旋运动的方法和设备。

背景技术

航天器姿态的测量、机动、指向等主要由其姿态控制系统负责，是航天器在轨供电、遥测遥控、科学目标实现等的基础。目前，地面主要是通过航天器载姿态敏感器的下行遥测数据，获得航天器姿态信息。一旦在轨航天器的姿态控制系统发生故障，航天器在内外力矩的作用下将逐渐形成稳定的自旋运动。在自旋状态下，航天器测控天线将失去稳定指向，在月球和深空探测任务等测控链路较为紧张的情况下，会导致地面接收航天器遥测“闪锁”甚至“失锁”。在这种情况下，地面无法通过遥测及时了解航天器自旋状态，更无法实施有效的姿态控制。若长期自旋，航天器将面临电源耗尽、结构解体等重大风险。因此，在航天器姿态故障导致自旋运动，地面无稳定遥测数据的情况下，必须及时估计航天器的自旋运动状态，这是故障诊断和故障处置的基本前提。

自旋运动导致航天器测控天线相对于航天器质心产生周期性的运动，会引起航天器(测控天线)对地面测控设备的附加多普勒，并表现在测速数据中；而且，通常情况下，地面测控设备锁定探测器下行载波比下行遥测的链路余量更大(对深空测控一般约大10dB)。换言之，在航天器自旋运动时，地面测控设备通常仍能锁定其下行载波，获取有效的测速数据，可望用于航天器自旋运动估计。

发明内容

基于前述背景，本发明针对航天器自旋故障导致地面无稳定遥测数据的情况下的自旋运动估计问题，选定航天器的自旋角速度大小和方向为参数建立测速数据观测模型，提出了一种基于谱估计确定自旋角速度大小，基于间接平差估计自旋角速度方向的航天器自旋运动估计方法。

为此，在一个方面中，本发明提供一种利用测速数据估计航天器自旋运动的方法，包括：

步骤1)：根据航天器精密轨道，确定航天器质心相对于多个地面测控设备在观测弧段内相对于航天器的速度并确定由于所述自旋运动引起的附加速度；

步骤2)：针对每个地面测控设备得到的航天器自旋运动速度进行变换，以获得谱估计；

步骤3)：针对每个地面测控设备的谱估计，确定其最强频谱分量对应的频率，以确定航天器自旋角速度大小的估计值；

步骤4)：基于所述估计值建立航天器自旋运动速度的观测方程；

步骤5)：选择状态向量X，并且基于观测方程和状态向量推导Jacobi矩阵；

步骤6)：基于间接平差方法确定所述状态向量X的改正量

步骤7)：基于迭代操作更新状态向量X的估计值，其中所述估计值

X⁰表示状态向量X的初始估计值；以及

步骤8)：将迭代操作结束后的估计值的第一至第三元素确定为自旋角速度的方向在航天器本体坐标系三个坐标轴上投影的估计值。

在一个实施例中，其中步骤1)包括：

根据航天器精密轨道，计算航天器质心相对于地面测控设备n(n＝1，2，...，N)在时刻idt对航天器的速度v_orbit(n，i)，然后根据式(1)计算由于自旋运动引起的附加速度v_obs.spin(n，i)。

v_obs.spin(n，i)＝v_obs(n，i)-v_orbit(n，i) (1)

在另一个实施例中，其中步骤2)包括：

根据式(2)对每个地面测控设备得到的航天器自旋运动速度进行傅里叶变换，以完成谱估计，

P(n，f)＝FT(v_obs.spin(n，i)) (2)

其中FT表示基于傅里叶变换进行谱估计。

在一个实施例中，其中步骤3)包括：

对每一地面测控设备的测量数据所得频谱，找到其最强频谱分量对应的频率f_n，并且根据下式(3)确定航天器自旋角速度大小的估计值，

在又一个实施例中，其中步骤4)包括：

根据式(4)建立航天器自旋运动速度的观测方程，

其中I₃为3×3的单位矩阵，e_r(n，i)为地面测控设备n在时刻idt至航天器质心方向的单位向量，H₀＝[H₀₁ H₀₂ H₀₃]^T是与时间、地面测控设备无关的常数未知向量，E_ω根据下式(5)确定，

其中e_ωx、e_ωy、e_ωz应满足如下式(6)的约束条件，

在一个实施例中，其中步骤5)包括：

选择如下的待定参数作为状态向量X，

X＝[e_ωx e_ωy e_ωz H₀₁ H₀₂ H₀₃]^T (7)

确定所述观测方程对状态向量的Jacobi矩阵B(n，i)为下式(8)

在又一个实施例中，其中步骤6)包括：

利用如下式(9)，求解状态向量X的改正量

(向量)，

其中[·]_NI表示将地面测控设备1～N、时刻0～(I-1)的全部数据组成NI×1的列向量，并且组成列向量的顺序不限，并且式(9)中所有[·]_NI的顺序一致，[·]|_X0表示计算[·]在X₀处的结果，P为地面测控设备1～N、时刻0～(I-1)的全部数据的权矩阵(NI×NI)方阵，元素顺序与[·]_NI的元素排列顺序一致，X₀为状态向量X的初始估计值，在缺少足够先验信息的情况下，按照下式(10)取值，

其中，A₀为初始时刻i＝0时，航天器本体坐标系O_b-X_bY_bZ_b至地球赤道惯性坐标系O-XYZ的坐标转换矩阵，并且若无先验信息，则取为3×3单位矩阵，ρ_antenna为航天器天线在本体坐标系O_b-X_bY_bZ_b下的安装位置。

在另外的实施例中，其中步骤7)包括：

更新参数估计值，迭代求解直至收敛，包括：

判断参数改正量

的每个元素的绝对值与收敛门限Tol的大小：

当

退出迭代；

否则，取

作为新的初始化值X⁰，根据式(9)计算对应新的初始化值的参数改正量

当迭代结束，根据下式(11)得到参数的估计值：

其中

的第1、2、3元素表示自旋角速度的方向在航天器本体坐标系三个坐标轴上的投影e_ωx、e_ωy、e_ωz的估计值。

在另一个方面中，本发明提供一种利用测速数据估计航天器自旋运动的设备，包括：

处理器；

存储器，其包括计算机指令，当该计算机指令由所述处理器执行时，使得所述设备执行以下步骤：

步骤1)：根据航天器精密轨道，确定航天器质心相对于多个地面测控设备在观测弧段内相对于航天器的速度并确定由于所述自旋运动引起的附加速度；

步骤2)：针对每个地面测控设备得到的航天器自旋运动速度进行变换，以获得谱估计；

步骤3)：针对每个地面测控设备的谱估计，确定其最强频谱分量对应的频率，以确定航天器自旋角速度大小的估计值；

步骤4)：基于所述估计值建立航天器自旋运动速度的观测方程；

步骤5)：选择状态向量X，并且基于观测方程和状态向量推导Jacobi矩阵；

步骤6)：基于间接平差方法确定所述状态向量X的改正量

步骤7)：基于迭代操作更新状态向量X的估计值，其中所述估计值

X⁰表示状态向量X的初始估计值；以及

步骤8)：将迭代操作结束后的估计值的第一至第三元素确定为自旋角速度的方向在航天器本体坐标系三个坐标轴上投影的估计值。

在又一个方面中，本发明提供一种计算机可读存储介质，其包括由处理器执行的用于利用测速数据估计航天器自旋运动的程序，当所述程序由处理器执行时，执行以下的操作步骤：

步骤1)：根据航天器精密轨道，确定航天器质心相对于多个地面测控设备在观测弧段内相对于航天器的速度并确定由于所述自旋运动引起的附加速度；

步骤2)：针对每个地面测控设备得到的航天器自旋运动速度进行变换，以获得谱估计；

步骤3)：针对每个地面测控设备的谱估计，确定其最强频谱分量对应的频率，以确定航天器自旋角速度大小的估计值；

步骤4)：基于所述估计值建立航天器自旋运动速度的观测方程；

步骤5)：选择状态向量X，并且基于观测方程和状态向量推导Jacobi矩阵；

步骤6)：基于间接平差方法确定所述状态向量X的改正量

步骤7)：基于迭代操作更新状态向量X的估计值，其中所述估计值

X⁰表示状态向量X的初始估计值；以及

步骤8)：将迭代操作结束后的估计值的第一至第三元素确定为自旋角速度的方向在航天器本体坐标系三个坐标轴上投影的估计值。

本发明的技术方案通过航天器测速数据估计航天器自旋运动状态，适应性强，收敛迅速。该方法不需要航天器遥测信息，可为航天器自旋故障情况下的故障诊断和处置提供有效依据。

附图说明

通过阅读仅作为示例提供并且参考附图进行的以下描述，将更好地理解本发明及其优点，其中：

图1示出根据本发明实施例的利用测速数据估计航天器自旋运动的方法的流程图；

图2示出根据本发明实施例的航天器自旋角速度在本体坐标系下的投影图示；

图3示出根据本发明实施例的地面测控设备对航天器测速的几何示意图；

图4示出根据本发明的地面测控设备得到的航天器自旋运动速度及频谱图示；以及

图5示出利用测速数据估计航天器自旋角速度方向的收敛过程。

具体实施方式

本发明基于前述背景，对航天器自旋故障导致地面无稳定遥测数据的情况下的自旋运动估计问题，选定航天器的自旋角速度大小和方向为参数建立测速数据观测模型，提出了一种基于谱估计确定自旋角速度大小，基于间接平差估计自旋角速度方向的航天器自旋运动估计方法。

为了便于对本发明的深入理解，首先对坐标系的建立进行介绍。在一个实施例中，本发明的方案选取地球赤道惯性坐标系O-XYZ(原点位于地心，X轴在地球赤道平面内指向春分点，Z轴与地球自转轴重合，Y轴与X轴、Z轴构成右手系)和航天器本体坐标系O_b-X_bY_bZ_b。航天器在本体坐标系O_b-X_bY_bZ_b下的自旋角速度为ω，其大小为ω，在本体坐标系三个坐标轴上的投影分别为ω_x、ω_y、ω_z；自旋角速度的方向为单位向量e_ω，在本体坐标系三个坐标轴上的投影分别为e_ωx、e_ωy、e_ωz，如图2中所示出的。

图1示出根据本发明实施例的利用测速数据估计航天器自旋运动的方法100的流程图。本发明实现的流程如下：

在航天器自旋运动期间，地面测控设备n(n＝1，2，...，N)在时刻i×dt(dt为采样间隔，i＝0，1，...，I-1)对航天器的测速结果为v_obs(n，i)。

在步骤1：根据航天器精密轨道，计算航天器质心相对于地面测控设备n(n＝1，2，...，N)在时刻idt对航天器的速度v_orbit(n，i)，然后根据式(1)计算由于自旋运动引起的附加速度v_obs.spin(n，i)。

v_obs.spin(n，i)＝v_obs(n，i)-v_orbit(n，i) (1)

步骤2：对每个地面测控设备得到的航天器自旋运动速度(时间序列)进行傅里叶变换(Fourier Transform)，完成谱估计，如式(2)所示。

P(n，f)＝FT(v_obs.spin(n，i)) (2)

在式(2)中，FT表示基于傅里叶变换进行谱估计(如周期图法等)。

步骤3：对每一地面测控设备的测量数据所得频谱，找到其最强频谱分量对应的频率f_n，得到航天器自旋角速度大小的估计值如式(3)所示。

步骤4：建立航天器自旋运动速度的观测方程如下：

在式(4)中，I₃为3×3的单位矩阵，e_r(n，i)为地面测控设备n在时刻idt至航天器质心方向的单位向量(如图3所示，在已知航天器精密轨道的情况下为已知量)，H₀＝[H₀₁ H₀₂H₀₃]^T是与时间、地面测控设备无关的常数未知向量，E_ω根据如下(5)确定，

此外，e_ωx、e_ωy、e_ωz应满足如下的约束条件，

步骤5：选择如下的待定参数作为状态向量X，

X＝[e_ωx e_ωy e_ωz H₀₁ H₀₂ H₀₃]^T (7)

推导步骤4中的观测方程对状态向量的Jacobi矩阵B(n，i)，即：

步骤6：根据间接平差方法(例如“附有约束条件的间接平差”)的相关模型和理论，利用如下公式，求解状态向量X的改正量

(向量)，

在式(9)中，[·]_NI表示将地面测控设备1～N、时刻0～(I-1)的全部数据组成NI×1的列向量(组成列向量的顺序不限，但式(9)中所有[·]_NI的顺序应一致)，[·]|_X0表示计算[·]在X₀处的结果，P为地面测控设备1～N、时刻0～(I-1)的全部数据的权矩阵(NI×NI方阵，元素顺序与[·]_NI的元素排列顺序一致)，X₀为状态向量X的初始估计值，在缺少足够先验信息的情况下，可按照下式取值，

在式(10)中，A₀为初始时刻(i＝0)航天器本体坐标系O_b-X_bY_bZ_b至地球赤道惯性坐标系O-XYZ的坐标转换矩阵(若无先验信息，则取为3×3单位矩阵)，ρ_antenna为航天器天线在本体坐标系O_b-X_bY_bZ_b下的安装位置。

步骤7：更新参数估计值，迭代求解直至收敛。判断参数改正量

的每个元素的绝对值与收敛门限Tol(根据实际需要设定，参考值1×10^-8)的大小：

1、若

退出迭代；

2、否则，取

作为新的初始化值X⁰，根据式(9)计算计算对应新的初始化值的参数改正量

迭代结束，得到参数的估计值：

该参数(向量)的第1、2、3元素即自旋角速度的方向在航天器本体坐标系三个坐标轴上的投影e_ωx、e_ωy、e_ωz的估计值。

上面结合附图描述了本发明的技术方案及其多个实施例，以下将结合我国探月任务的典型地月转移轨道，模拟产生航天器发生自旋故障情况下两个测站的测速数据，对本发明的具体实施方式作进一步说明。

假设航天器的轨道根数为：半长轴196478km，偏心率0.9665，轨道倾角28.5°，升交点赤经205.5°，近地点幅角173.8°，真近点角150.0°(初始时刻)。天线安装位置ρ_antenna＝[-1.0 2.0 -1.5]^T(m)。自旋运动状态为：ω_x＝60°/s、ω_y＝-20°/s、ω_z＝10°/s，即ω＝1.11756rad/s、e_ωx＝0.9370、e_ωy＝-0.3123、e_ωz＝0.1562。可获得的测速数据为地面测控设备1、2(即N＝2)的测量结果，地面测控设备1、2在地球固连坐标系下的坐标分别为(-2872，3331，4603)km、(1150，4870，3943)km。基于以上数据，模拟生成idt(dt＝0.2s，i＝0，1，2，3，...，600，即I＝601)时刻地面测控设备1、2对航天器的测速数据v_obs(n，i)，并在测速数据中添加零均值白噪声(1σ＝5×10^-3m/s)。

基于以上模型及数据，按照如下流程进行航天器自旋运动估计。

步骤2-1：根据航天器轨道根数，计算航天器质心相对于地面测控设备1、2在时刻idt对航天器的速度v_orbit(n，i)，然后根据式(12)计算由于自旋运动引起的附加速度v_obs.spin(n，i)。

v_obs.spin(n，i)＝v_obs(n，i)-v_orbit(n，i) (12)

为了便于理解，图4上图所示为地面测控设备1得到的由于航天器自旋运动引起的附加速度。

步骤2-1：对地面测控设备1、2得到的航天器自旋运动速度(时间序列)进行傅里叶变换(Fourier Transform)，采用Welch方法完成谱估计，如式(13)所示。

P(n，f)＝FT(v_obs.spin(n，i)) (13)

在式(13)中，FT表示基于傅里叶变换进行谱估计。为了便于理解，图4下图所示为对地面测控设备1得到的航天器自旋运动速度进行谱估计得到的频谱。

步骤2-3：根据地面测控设备1、2的测量数据所得频谱，其最强频谱分量对应的频率分别为f₁＝0.17790Hz、f₂＝0.17785Hz，从而航天器自旋角速度大小的估计值为：

步骤2-4：建立航天器自旋运动速度的观测方程如下：

在式(15)中，I₃为3×3的单位矩阵，e_r(n，i)为地面测控设备n在时刻idt至航天器质心方向的单位向量(如图4所示，其为已知量)，H₀＝[H₀₁ H₀₂ H₀₃]^T是与时间、地面测控设备无关的常数未知向量，E_ω可以如下表达：

此外，e_ωx、e_ωy、e_ωz应满足如下的约束条件，

步骤2-5：选择如下的待定参数作为状态向量X，

X＝[e_ωx e_ωy e_ωz H₀₁ H₀₂ H₀₃]^T (18)

推导步骤2-4中的观测方程对状态向量的Jacobi矩阵B(n，i)，即：

步骤2-6：根据“附有约束条件的间接平差”的相关模型和理论，利用如下公式，求解状态向量X的改正量

(向量)，

在式(20)中，[·]₁₂₀₂表示将地面测控设备1～2、时刻0～600的全部数据组成1202×1的列向量(第1-601行是地面测控设备1的测量数据，第602-1202行是地面测控设备2的测量数据)，[·]|_X0表示计算[·]在X₀处的结果，P为全部测量数据的权矩阵(1202×1202方阵，第1-601行/列表示地面测控设备1的测量数据，第602-1202行/列表示地面测控设备2的测量数据)，X₀为状态向量X的初始估计值，根据当前的先验信息，取值如下，

步骤2-7：更新参数估计值，迭代求解直至收敛。根据式(20)，经过一次求解后得到：

设定收敛门限Tol为1×10^-8，因

则取

作为新的初始化值X⁰，根据式(20)计算计算对应新的初始化值的参数改正量

重复上述判断和计算步骤，直至第6次迭代计算达到收敛条件，迭代结束，得到参数的估计值：

即自旋角速度的方向在航天器本体坐标系三个坐标轴上的投影e_ωx、e_ωy、e_ωz的估计值分别为0.9365、-0.3117、0.1604。

图5示出利用测速数据估计航天器自旋角速度方向的收敛过程。具体地，图5所示为e_ωx、e_ωy、e_ωz估计值的迭代收敛过程，图中横轴为迭代次数，纵轴为航天器自旋角速度方向在本体系下的投影；同时，图中以黑色虚线显示了自旋角速度方向的真实值。从图中直观可见，航天器自旋角速度方向的参数估计迅速收敛至真实值附近。

本例中，航天器的自旋角速度大小和方向的估计值和真实值的比对情况见下表1所示。从表1中看出，航天器自旋角速度大小和方向的估计较为准确，能够满足航天器自旋故障情况下的故障诊断和处置需要。

表1 航天器的自旋角速度大小和方向的估计值和真实值

在一些实施例中，本发明的方案还可以通过在计算机可读的记录介质以计算机可读代码来体现。计算机可读记录介质包括存储可通过计算机系统解读的数据的所有种类的记录介质。该记录介质例如可以包括但不限于只读存储器(ROM，“Read Only Memory”)、随机存取存储器(RAM，“Random Access Memory”)、磁盘、磁盘、光盘、闪存等。进一步，这些计算机可读的记录介质可以通过通信网络(包括计算机通信网络、蜂窝通信网络或局域域通信网络)在各个通信实体之间传播或扩散，从而也可以通过任意的方式来运行存储在计算机可读存储介质上的计算机可读指令或计算机可执行代码。

虽然本发明所实施的方式如上，但所述内容只是为便于理解本发明而采用的实施例，并非用以限定本发明的范围和应用场景。任何本发明所述技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种利用测速数据估计航天器自旋运动的方法，包括：

步骤1)：根据航天器精密轨道，确定航天器质心相对于多个地面测控设备在观测弧段内相对于航天器的速度并确定由于所述自旋运动引起的附加速度；

步骤2)：针对每个地面测控设备得到的航天器自旋运动速度进行变换，以获得谱估计；

步骤3)：针对每个地面测控设备的谱估计，确定其最强频谱分量对应的频率，以确定航天器自旋角速度大小的估计值；

步骤4)：基于所述估计值建立航天器自旋运动速度的观测方程；

步骤5)：选择状态向量X，并且基于观测方程和状态向量推导Jacobi矩阵；

步骤6)：基于间接平差方法确定所述状态向量X的改正量

步骤7)：基于迭代操作更新状态向量X的估计值，其中所述估计值

X⁰表示状态向量X的初始估计值；以及

步骤8)：将迭代操作结束后的估计值的第一至第三元素确定为自旋角速度的方向在航天器本体坐标系三个坐标轴上投影的估计值。

2.根据权利要求1所述的方法，其中步骤1)包括：

根据航天器精密轨道，计算航天器质心相对于地面测控设备n(n＝1，2，...，N)在时刻idt对航天器的速度v_orbit(n，i)，然后根据式(1)计算由于自旋运动引起的附加速度v_obs.spin(n，i)：

v_obs.spin(n，i)＝v_obs(n，i)-v_orbit(n，i) (1)

其中，v_obs(n，i)为地面测控设备n在时刻idt对航天器的测速结果。

3.根据权利要求1所述的方法，其中步骤2)包括：

根据式(2)对每个地面测控设备得到的航天器自旋运动速度进行傅里叶变换，以完成谱估计，

P(n，f)＝FT(v_obs.spin(n，i)) (2)

其中，v_obs.spin(n，i)表示由于自旋运动引起的附加速度，FT表示基于傅里叶变换进行谱估计。

4.根据权利要求1所述的方法，其中步骤3)包括：

对每一地面测控设备的测量数据所得频谱，找到其最强频谱分量对应的频率f_n，并且根据下式(3)确定航天器自旋角速度大小的估计值，

其中，N为地面测控设备的总数。

5.根据权利要求1所述的方法，其中步骤4)包括：

根据式(4)建立航天器自旋运动速度的观测方程，

其中I₃为3×3的单位矩阵，e_r(n，i)为地面测控设备n在时刻idt至航天器质心方向的单位向量，H₀＝[H₀₁ H₀₂ H₀₃]^T是与时间、地面测控设备无关的常数未知向量，

为航天器自旋角速度大小的估计值，E_ω根据下式(5)确定，

其中，e_ωx、e_ωy、e_ωz分别为自旋角速度方向在航天器本体坐标系三个坐标轴上的投影，应满足如下式(6)的约束条件，

6.根据权利要求1所述的方法，其中步骤5)包括：

选择如下的待定参数作为状态向量X，

X＝[e_ωx e_ωy e_ωz H₀₁ H₀₂ H₀₃]^T (7)

其中，e_ωx、e_ωy、e_ωz分别为自旋角速度方向在航天器本体坐标系三个坐标轴上的投影，H₀＝[H₀₁ H₀₂ H₀₃]^T是与时间、地面测控设备无关的常数未知向量；

确定所述观测方程对状态向量的Jacobi矩阵B(n，i)为下式(8)

其中，v_spin(n，i)表示航天器自旋运动速度的观测方程。

7.根据权利要求1所述的方法，其中步骤6)包括：

利用如下式(9)，求解状态向量X的改正量

(向量)，

其中，B(n，i)为所述观测方程对状态向量的Jacobi矩阵；e_ωx、e_ωy、e_ωz分别为自旋角速度方向在航天器本体坐标系三个坐标轴上的投影，[·]_NI表示将地面测控设备1～N、时刻0～(I-1)的全部数据组成NI×1的列向量，并且组成列向量的顺序不限，并且式(9)中所有[·]_NI的顺序一致，[·]|_X0表示计算[·]在X₀处的结果，P为地面测控设备1～N、时刻0～(I-1)的全部数据的权矩阵(NI×NI)方阵，元素顺序与[·]_NI的元素排列顺序一致，X₀为状态向量X的初始估计值，在缺少足够先验信息的情况下，按照下式(10)取值，

其中，

为航天器自旋角速度大小的估计值，v_obs.spin(n，i)表示由于自旋运动引起的附加速度，A₀为初始时刻i＝0时，航天器本体坐标系O_b-X_bY_bZ_b至地球赤道惯性坐标系O-XYZ的坐标转换矩阵，并且若无先验信息，则取为3×3单位矩阵，ρ_antenna为航天器天线在本体坐标系O_b-X_bY_bZ_b下的安装位置；v_spin(n，i)表示航天器自旋运动速度的观测方程。

8.根据权利要求7所述的方法，其中步骤7)包括：

更新参数估计值，迭代求解直至收敛，包括：

判断参数改正量

的每个元素的绝对值与收敛门限Tol的大小：

当

退出迭代；

否则，取

作为新的初始化值X⁰，根据所述式(9)：

N_BB＝B^TPB

W＝B^TPl

计算对应新的初始化值的参数改正量

当迭代结束，根据下式(11)得到参数的估计值：

其中

的第1、2、3元素表示自旋角速度的方向在航天器本体坐标系三个坐标轴上的投影e_ωx、e_ωy、e_ωz的估计值。

9.一种利用测速数据估计航天器自旋运动的设备，包括：

处理器；

存储器，其包括计算机指令，当该计算机指令由所述处理器执行时，使得所述设备执行以下步骤：

步骤1)：根据航天器精密轨道，确定航天器质心相对于多个地面测控设备在观测弧段内相对于航天器的速度并确定由于所述自旋运动引起的附加速度；

步骤2)：针对每个地面测控设备得到的航天器自旋运动速度进行变换，以获得谱估计；

步骤3)：针对每个地面测控设备的谱估计，确定其最强频谱分量对应的频率，以确定航天器自旋角速度大小的估计值；

步骤4)：基于所述估计值建立航天器自旋运动速度的观测方程；

步骤5)：选择状态向量X，并且基于观测方程和状态向量推导Jacobi矩阵；

步骤6)：基于间接平差方法确定所述状态向量X的改正量

步骤7)：基于迭代操作更新状态向量X的估计值，其中所述估计值

X⁰表示状态向量X的初始估计值；以及

步骤8)：将迭代操作结束后的估计值的第一至第三元素确定为自旋角速度的方向在航天器本体坐标系三个坐标轴上投影的估计值。

10.一种计算机可读存储介质，其包括由处理器执行的用于利用测速数据估计航天器自旋运动的程序，当所述程序由处理器执行时，执行以下的操作步骤：

步骤1)：根据航天器精密轨道，确定航天器质心相对于多个地面测控设备在观测弧段内相对于航天器的速度并确定由于所述自旋运动引起的附加速度；

步骤2)：针对每个地面测控设备得到的航天器自旋运动速度进行变换，以获得谱估计；

步骤3)：针对每个地面测控设备的谱估计，确定其最强频谱分量对应的频率，以确定航天器自旋角速度大小的估计值；

步骤4)：基于所述估计值建立航天器自旋运动速度的观测方程；

步骤5)：选择状态向量X，并且基于观测方程和状态向量推导Jacobi矩阵；

步骤6)：基于间接平差方法确定所述状态向量X的改正量

步骤7)：基于迭代操作更新状态向量X的估计值，其中所述估计值

X⁰表示状态向量X的初始估计值；以及

步骤8)：将迭代操作结束后的估计值的第一至第三元素确定为自旋角速度的方向在航天器本体坐标系三个坐标轴上投影的估计值。

Images