CN109470266A - 一种处理乘性噪声的星敏感器陀螺组合定姿方法 - Google Patents

Abstract

本发明公开了一种处理乘性噪声的星敏感器陀螺组合定姿方法，包括确定乘性安装误差的表达式以及统计特性；建立具有星敏感器乘性安装误差的卫星姿态确定系统模型；将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型；对所述线性离散系统模型的系统状态进行时间更新，得到状态滤波参数；将状态滤波参数中的误差四元数转换为姿态欧拉角。本发明基于误差四元数建立具有乘性安装误差的星敏感器/陀螺组合的卫星姿态确定系统模型，利用线性最小方差估计和投影定理，设计基于线性最小方差的最优姿态滤波算法，利用此算法实现对星敏感器乘性测量噪声的抑制和卫星三轴运动姿态的估计，有效降低乘性测量噪声对姿态确定精度的影响。

Description

一种处理乘性噪声的星敏感器陀螺组合定姿方法

技术领域

本发明涉及卫星分析技术领域，更具体地说涉及一种用于处理乘性噪声的星敏感器/陀螺组合定姿方法。

背景技术

卫星姿态确定的主要任务是根据带有噪声的姿态敏感器的测量值估计卫星相对于某个参考坐标系的姿态参数。姿态确定系统主要由姿态敏感器和相应的信息处理算法即姿态确定算法组成，姿态确定精度取决于姿态敏感器测量精度和姿态确定算法精度。

卫星姿态确定系统是卫星姿态控制系统中的重要组成部分，高精度的姿态确定是高分辨率成像、高精度测绘等卫星应用的重要基础和关键技术。

星敏感器的姿态测量精度最高，但受其采样频率的限制，仅利用角度敏感器无法精确获得卫星姿态的动态变化信息，一般还需引入高精度陀螺，连续测量星体姿态角速度，作为星体姿态的基准。因此，星敏感器与陀螺组成的组合定姿系统成为现代高精度卫星姿态确定的主要手段，在卫星的姿态测量系统中得到了广泛的实际应用。

为提高星敏感器/陀螺组合的定姿精度，国内外许多学者对此做过大量的研究。如针对陀螺/星敏感器组合方案，利用状态估计法设计相应的Kalman或其改进的Kalman姿态滤波算法；针对星敏感器/陀螺在轨的常值误差、陀螺漂移误差以及测量噪声等，设计相应的星敏感器在轨标定模型和算法。但这些研究未考虑星敏感器的安装误差，或者只把安装误差作为常值误差进行处理，而未考虑安装误差为随机噪声的情况。

如果星敏感器相对于卫星本体的安装矩阵在卫星的运动过程中发生了偏转，或者由于抖动，振动等在轨环境的影响，会使得星敏感器测量方程中引入另外一个未知的坐标变换。因此，星敏感器的测量随机噪声在卫星姿态确定系统的观测方程中往往体现为乘性噪声。

与传统Kalman或改进的Kalman滤波模型中的加性噪声，乘性噪声的方差依赖于系统的状态变量，而状态变量的真值通常是未知的，这给滤波器的设计带来了困难。

针对带乘性噪声的动态系统，目前也存在一些理论研究，如针对一维乘性噪声系统，P.K.Rajasekaram提出了一种处理乘性噪声的最优线性估计方法，给出了最优滤波算法，有效地解决了存在乘性噪声情况下的线性问题。褚东升等在此基础上，分别将线性系统的乘性噪声推广为对角阵和随机阵。F.Wang,V.Balakrishnan等人对乘性噪声的线性系统设计了鲁棒Kalman滤波算法，通过凸优化的方法确定滤波参数，并采用线性矩阵不等式技术求解凸优化问题。还给出了另外一种用于处理乘性噪声的鲁棒滤波算法，适合在线应用需求。

但针对带随机乘性噪声的非动态线性系统还没有人做过相关研究，且目前的研究也未涉及具体的应用层面，尤其是由星敏感器/陀螺组合的卫星姿态确定系统。

综上所述，星敏感器/陀螺组合定姿系统中，由于星敏感器在轨安装过程中会产生安装误差，而这种安装误差在星敏感器的测量方程中体现为乘性噪声形式，设计处理乘性噪声的星敏感器/陀螺组合定姿方法，降低这种误差对姿态确定精度的影响，是提高星敏感器/陀螺组合姿态确定系统精度的技术难点之一。

发明内容

本发明要解决的技术问题是：如何提高卫星姿态确定精度。

本发明解决其技术问题的解决方案是：

一种处理乘性噪声的星敏感器陀螺组合定姿方法，包括以下步骤：

步骤1，确定乘性安装误差的表达式以及统计特性；

步骤2，根据卫星姿态误差四元数的状态方程和乘性安装误差，建立具有星敏感器乘性安装误差的卫星姿态确定系统模型；

步骤3，将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型；

步骤4，对所述线性离散系统模型的系统状态进行时间更新，并基于带有乘性安装误差的星敏感器测量数据对系统状态进行测量更新，得到状态滤波参数；

步骤5，将状态滤波参数中的误差四元数转换为姿态四元数，再将姿态四元数进一步转换为姿态欧拉角，得到卫星在轨姿态确定结果。

作为上述技术方案的进一步改进，采用单个星敏感器进行测量时，所述乘性安装误差的表达式如式1所示：

其中δψ、δφ以及δθ分别表示偏航、滚动以及俯仰三个方向的安装误差角，式1的近似形式如式2所示：

对于利用三个星敏感器进行测量时，所述乘性安装误差的表达式的近似形式如式3所示：

其中，0表示零矩阵。

作为上述技术方案的进一步改进，所述乘性安装误差的统计特性为：乘性安装误差U的均值为M，协方差阵为N，协方差阵N＝(N_ij)_9×9，其中N_ij表示乘性安装误差U第i行与第j列元素的方差。

作为上述技术方案的进一步改进，所述步骤2具体包括以下步骤：

步骤2.1，根据乘性安装误差，获取具有星敏感器乘性安装误差的测量方程；

步骤2.2，根据卫星姿态误差四元数的状态方程以及测量方程，得到具有星敏感器乘性安装误差的卫星姿态确定系统模型。

作为上述技术方案的进一步改进，卫星姿态误差四元数的状态方程如式4所示：

其中系统状态变量X_9×1＝(Δq₁₃ ^T Δd^T Δb^T)，Δq＝[Δq₀,(Δq₁₃)^T]^T为误差四元数，Δq₁₃＝[Δq₁ Δq₂ Δq₃]^T表示误差四元数的矢量部分，Δd＝[Δd₁ Δd₂ Δd₃]^T表示陀螺的相关漂移误差，Δb＝[Δb₁ Δb₂ Δb₃]^T表示陀螺的常值漂移误差，t为时间参数；

表示过程噪声，为通过陀螺测量输出的卫星角速度的估计值；矩阵D_τ是由相关时间常数τ_i(i＝x,y,z)构成的对角阵，n_g表示陀螺的测量噪声，n_d和n_b分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，过程噪声w(t)的协方差为Q(t)。

作为上述技术方案的进一步改进，具有星敏感器乘性安装误差的测量方程如式5所示：

Z(t)＝h(X,t)+v(t) 式5

其中δC为星敏感器的乘性安装误差，A(q)为四元数姿态矩阵，q＝[q₀ q₁ q₂ q₃]^T为姿态四元数，μ_I1＝[1 0 0]^T，μ_I2＝[0 1 0]^T，μ_I3＝[0 0 1]^T分别为不同的星敏感器的安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。

作为上述技术方案的进一步改进，具有星敏感器乘性安装误差的卫星姿态确定系统模型为

作为上述技术方案的进一步改进，步骤3中将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型如式6所示：

其中，状态转移矩阵 T为滤波器周期，为姿态四元数的预报值，为状态预测值，k为离散化的时间参数。

作为上述技术方案的进一步改进，步骤4具体包括：设置Kalman增益K(k)，其中P(k/k-1)＝Φ(k,k-1)P(k-1/k-1)Φ^T(k,k-1)+Q(k-1)，R_L(k)＝Y(k)+G(k)+M(k)H(k)P(k/k-1)H^T(k)M^T(k)+R(k)，Y(k)＝E[(U(k)-M(k))H(k)X(k)X^T(k)H^T(k)(U(k)-M(k))]＝(y_tl(k))_r×r，a_ij为矩阵H(k)X(k)X^T(k)H^T(k)第i行第j列的元素，Q(k-1)为过程噪声的方差，R(k)为星敏感器测量方程中加性噪声的方差，S(k)＝Φ(k,k-1)S(k-1)Φ^T(k,k-1)+Q(k-1)，G(k)＝E[(U(k)-M(k))L(k)L^T(k)(U(k)-M(k))]＝(g_tl(k))_r×r，b_ij为矩阵L(k)L^T(k)第i行第j列的元素，状态滤波值的测量更新为X(k/k)＝X(k/k-1)+K(k)(Z(k)-M(k)H(k)X(k/k-1))。

作为上述技术方案的进一步改进，将状态滤波值转换成姿态四元数过程如下：令 则状态滤波值为姿态四元数的滤波值，通过姿态四元数可得：偏航角俯仰角滚动角

本发明的有益效果是：本发明基于误差四元数建立具有乘性安装误差的星敏感器/陀螺组合的卫星姿态确定系统模型，利用线性最小方差估计和投影定理，设计基于线性最小方差的最优姿态滤波算法，利用此算法实现对星敏感器乘性测量噪声的抑制和卫星三轴运动姿态的估计，有效降低乘性测量噪声对姿态确定精度的影响。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单说明。显然，所描述的附图只是本发明的一部分实施例，而不是全部实施例，本领域的技术人员在不付出创造性劳动的前提下，还可以根据这些附图获得其他设计方案和附图。

图1是本发明的方法流程图；

图2是根据现有的技术方案得到的卫星姿态确定结果；

图3是根据本技术方案得到的卫星姿态确定结果。

具体实施方式

以下将结合实施例和附图对本发明的构思、具体结构及产生的技术效果进行清楚、完整的描述，以充分地理解本发明的目的、特征和效果。显然，所描述的实施例只是本发明的一部分实施例，而不是全部实施例，基于本发明的实施例，本领域的技术人员在不付出创造性劳动的前提下所获得的其他实施例，均属于本发明保护的范围。

参照图1，本发明申请公开了一种处理乘性噪声的星敏感器陀螺组合定姿方法，包括以下步骤：

步骤1，确定乘性安装误差的表达式以及统计特性；

步骤2，根据卫星姿态误差四元数的状态方程和乘性安装误差，建立具有星敏感器乘性安装误差的卫星姿态确定系统模型；

步骤3，将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型；

步骤4，对所述线性离散系统模型的系统状态进行时间更新，并基于带有乘性安装误差的星敏感器测量数据对系统状态进行测量更新，得到状态滤波参数；

步骤5，将状态滤波参数中的误差四元数转换为姿态四元数，再将姿态四元数进一步转换为姿态欧拉角，得到卫星在轨姿态确定结果。

具体地，本发明基于误差四元数建立具有乘性安装误差的星敏感器/陀螺组合的卫星姿态确定系统模型，利用线性最小方差估计和投影定理，设计基于线性最小方差的最优姿态滤波算法，利用此算法实现对星敏感器乘性测量噪声的抑制和卫星三轴运动姿态的估计，有效降低乘性测量噪声对姿态确定精度的影响。

进一步作为优选的实施方式，本发明创造具体实施方式中，采用单个星敏感器进行测量时，所述乘性安装误差的表达式如式1所示：

其中δψ、δφ以及δθ分别表示偏航、滚动以及俯仰三个方向的安装误差角，由于δψ，δφ，δθ均为小角度，可以认为cosδψ＝cosδφ＝cosδθ＝1，sinδψ＝δψ，sinδφ＝δφ，sinδθ＝δθ，同时忽略二阶小量，可得单个星敏感器的乘性安装误差如式2所示：

对于利用三个星敏感器进行测量时，所述乘性安装误差的表达式的近似形式如式3所示：

其中，0表示零矩阵。

进一步作为优选的实施方式，本发明创造具体实施方式中，所述乘性安装误差的统计特性为：乘性安装误差U的均值为M，协方差阵为N，协方差阵N＝(N_ij)_9×9，其中N_ij表示乘性安装误差U第i行与第j列元素的方差。

进一步作为优选的实施方式，本发明创造具体实施方式中，所述步骤2具体包括以下步骤：

步骤2.1，根据乘性安装误差，获取具有星敏感器乘性安装误差的测量方程；

步骤2.2，根据卫星姿态误差四元数的状态方程以及测量方程，得到具有星敏感器乘性安装误差的卫星姿态确定系统模型。

进一步作为优选的实施方式，本发明创造具体实施方式中，卫星姿态误差四元数的状态方程如式4所示：

其中系统状态变量X_9×1＝(Δq₁₃ ^T Δd^T Δb^T)，Δq＝[Δq₀,(Δq₁₃)^T]^T为误差四元数，Δq₁₃＝[Δq₁ Δq₂ Δq₃]^T表示误差四元数的矢量部分，Δd＝[Δd₁ Δd₂ Δd₃]^T表示陀螺的相关漂移误差，Δb＝[Δb₁ Δb₂ Δb₃]^T表示陀螺的常值漂移误差，t为时间参数；

表示过程噪声，为通过陀螺测量输出的卫星角速度的估计值；矩阵D_τ是由相关时间常数τ_i(i＝x,y,z)构成的对角阵，n_g表示陀螺的测量噪声，n_d和n_b分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，过程噪声w(t)的协方差为Q(t)。

进一步作为优选的实施方式，本发明创造具体实施方式中，具有星敏感器乘性安装误差的测量方程如式5所示：

Z(t)＝h(X,t)+v(t) 式5

其中δC为星敏感器的乘性安装误差，A(q)为四元数姿态矩阵，q＝[q₀ q₁ q₂ q₃]^T为姿态四元数，μ_I1＝[1 0 0]^T，μ_I2＝[0 1 0]^T，μ_I3＝[0 0 1]^T分别为不同的星敏感器的安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。

进一步作为优选的实施方式，本发明创造具体实施方式中，具有星敏感器乘性安装误差的卫星姿态确定系统模型为

进一步作为优选的实施方式，本发明创造具体实施方式中，步骤3中将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型如式6所示：

其中，状态转移矩阵 T为滤波器周期，为姿态四元数的预报值，为状态预测值，k为离散化的时间参数。

进一步作为优选的实施方式，本发明创造具体实施方式中，步骤4具体包括：设置Kalman增益K(k)，其中P(k/k-1)＝Φ(k,k-1)P(k-1/k-1)Φ^T(k,k-1)+Q(k-1)，R_L(k)＝Y(k)+G(k)+M(k)H(k)P(k/k-1)H^T(k)M^T(k)+R(k)，Y(k)＝E[(U(k)-M(k))H(k)X(k)X^T(k)H^T(k)(U(k)-M(k))]＝(y_tl(k))_r×r，a_ij为矩阵H(k)X(k)X^T(k)H^T(k)第i行第j列的元素，Q(k-1)为过程噪声的方差，R(k)为星敏感器测量方程中加性噪声的方差，S(k)＝Φ(k,k-1)S(k-1)Φ^T(k,k-1)+Q(k-1)，G(k)＝E[(U(k)-M(k))L(k)L^T(k)(U(k)-M(k))]＝(g_tl(k))_r×r，b_ij为矩阵L(k)L^T(k)第i行第j列的元素，状态滤波值的测量更新为X(k/k)＝X(k/k-1)+K(k)(Z(k)-M(k)H(k)X(k/k-1))。

进一步作为优选的实施方式，本发明创造具体实施方式中，将状态滤波值转换成姿态四元数过程如下：令 则状态滤波值 为姿态四元数的滤波值，通过姿态四元数可得：偏航角俯仰角滚动角

为更好地说明本技术方案的处理流程，下面结合应用实例对本发明申请进行详细说明。

步骤1中取E[δψ δφ δθ]＝[0 0 0]，则有算得均值M＝I_9×9，协方差阵N亦可通过上式得出；

步骤3中取滤波器周期T＝1秒，τ＝3600，得到线性离散系统模型，算出矩阵Φ(k,k-1)以及H(k)；

步骤4中取星敏感器测量精度为10角秒(3σ)，根据步骤一得出的协方差阵N，求出Kalman增益K(k)，根据星敏感器的测量值Z(k)计算测量更新后的状态滤波值X(k/k)。

步骤5中根据状态滤波值X(k/k)修正预报值以及得到以及进而把姿态四元数转换为姿态欧拉角。

参照图2和图3，分别给出了利用现有技术和本技术方案处理后得到的姿态欧拉角误差图，从上而下依次为滚动角、俯仰角以及偏航角，单位为角秒，表1给出了利用现有技术和本技术方案得到的欧拉角估计值的均方差。

均方差(3σ)	滚动角(单位：角秒)	俯仰角(单位：角秒)	偏航角(单位：角秒)
				应用本发明前	4.821	4.502	4.438
应用本发明后	4.629	4.321	4.198

表1

本发明创造中在星敏感器的测量方程中引入乘性安装误差，通过在滤波过程中改变增益矩阵来减少乘性安装误差带来的影响，在星敏感器测量精度为10角秒时，加入标准差为8角秒的乘性安装误差，采用本技术方案后得到的欧拉角估计值的均方差提高0.2角秒左右。

以上对本发明的较佳实施方式进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可作出种种的等同变型或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。

Claims (10)

1.一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，包括以下步骤：

步骤1，确定乘性安装误差的表达式以及统计特性；

步骤2，根据卫星姿态误差四元数的状态方程和乘性安装误差，建立具有星敏感器乘性安装误差的卫星姿态确定系统模型；

步骤3，将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型；

步骤4，对所述线性离散系统模型的系统状态进行时间更新，并基于带有乘性安装误差的星敏感器测量数据对系统状态进行测量更新，得到状态滤波参数；

步骤5，将状态滤波参数中的误差四元数转换为姿态四元数，再将姿态四元数进一步转换为姿态欧拉角，得到卫星在轨姿态确定结果。

2.根据权利要求1所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，采用单个星敏感器进行测量时，所述乘性安装误差的表达式如式1所示：

其中δψ、δφ以及δθ分别表示偏航、滚动以及俯仰三个方向的安装误差角，式1的近似形式如式2所示：

对于利用三个星敏感器进行测量时，所述乘性安装误差的表达式的近似形式如式3所示：

其中，0表示零矩阵。

3.根据权利要求2所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，所述乘性安装误差的统计特性为：乘性安装误差U的均值为M，协方差阵为N，协方差阵N＝(N_ij)_9×9，其中N_ij表示乘性安装误差U第i行与第j列元素的方差。

4.根据权利要求3所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，所述步骤2具体包括以下步骤：

步骤2.1，根据乘性安装误差，获取具有星敏感器乘性安装误差的测量方程；

步骤2.2，根据卫星姿态误差四元数的状态方程以及测量方程，得到具有星敏感器乘性安装误差的卫星姿态确定系统模型。

5.根据权利要求4所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，卫星姿态误差四元数的状态方程如式4所示：

其中系统状态变量X_9×1＝(Δq₁₃ ^T Δd^T Δb^T)，Δq＝[Δq₀,(Δq₁₃)^T]^T为误差四元数，Δq₁₃＝[Δq₁ Δq₂ Δq₃]^T表示误差四元数的矢量部分，Δd＝[Δd₁ Δd₂ Δd₃]^T表示陀螺的相关漂移误差，Δb＝[Δb₁ Δb₂ Δb₃]^T表示陀螺的常值漂移误差，t为时间参数；

表示过程噪声，为通过陀螺测量输出的卫星角速度的估计值；矩阵D_τ是由相关时间常数τ_i(i＝x,y,z)构成的对角阵，n_g表示陀螺的测量噪声，n_d和n_b分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，过程噪声w(t)的协方差为Q(t)。

6.根据权利要求5所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，具有星敏感器乘性安装误差的测量方程如式5所示：

Z(t)＝h(X,t)+v(t) 式5

其中δC为星敏感器的乘性安装误差，A(q)为四元数姿态矩阵，q＝[q₀ q₁ q₂ q₃]^T为姿态四元数，μ_I1＝[1 0 0]^T，μ_I2＝[0 1 0]^T，μ_I3＝[0 0 1]^T分别为不同的星敏感器的安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。

7.根据权利要求6所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，具有星敏感器乘性安装误差的卫星姿态确定系统模型为

8.根据权利要求7所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，步骤3中将所述卫星姿态确定系统模型进行离散化和线性化，得到线性离散系统模型如式6所示：

其中，状态转移矩阵 T为滤波器周期，为姿态四元数的预报值，为状态预测值，k为离散化的时间参数。

9.根据权利要求8所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，步骤4具体包括：设置Kalman增益K(k)，其中P(k/k-1)＝Φ(k,k-1)P(k-1/k-1)Φ^T(k,k-1)+Q(k-1)，R_L(k)＝Y(k)+G(k)+M(k)H(k)P(k/k-1)H^T(k)M^T(k)+R(k)，Y(k)＝E[(U(k)-M(k))H(k)X(k)X^T(k)H^T(k)(U(k)-M(k))]＝(y_tl(k))_r×r，

a_ij为矩阵H(k)X(k)X^T(k)H^T(k)第i行第j列的元素，Q(k-1)为过程噪声的方差，R(k)为星敏感器测量方程中加性噪声的方差，S(k)＝Φ(k,k-1)S(k-1)Φ^T(k,k-1)+Q(k-1)，G(k)＝E[(U(k)-M(k))L(k)L^T(k)(U(k)-M(k))]＝(g_tl(k))_r×r，

b_ij为矩阵L(k)L^T(k)第i行第j列的元素，状态滤波值的测量更新为X(k/k)＝X(k/k-1)+K(k)(Z(k)-M(k)H(k)X(k/k-1))。

10.根据权利要求9所述的一种处理乘性噪声的星敏感器陀螺组合定姿方法，其特征在于，将状态滤波值转换成姿态四元数过程如下：令 则状态滤波值 为姿态四元数的滤波值，通过姿态四元数可得：偏航角俯仰角滚动角