CN109470267A - 一种卫星姿态滤波方法 - Google Patents
一种卫星姿态滤波方法 Download PDFInfo
- Publication number
- CN109470267A CN109470267A CN201811300220.0A CN201811300220A CN109470267A CN 109470267 A CN109470267 A CN 109470267A CN 201811300220 A CN201811300220 A CN 201811300220A CN 109470267 A CN109470267 A CN 109470267A
- Authority
- CN
- China
- Prior art keywords
- satellite
- attitude
- state
- value
- filtering
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C25/00—Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
Landscapes
- Engineering & Computer Science (AREA)
- Manufacturing & Machinery (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Navigation (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
本发明公开了一种卫星姿态滤波方法，包括根据星敏感器的偏转乘性误差，构建具有偏转乘性噪声的卫星姿态确定系统模型；对卫星姿态确定系统模型进行第一阶段状态滤波，得到第一状态滤波值以及状态的一步平滑值；对卫星姿态确定系统模型进行第二阶段状态滤波，得到第二状态滤波值；将第二状态滤波值的姿态四元数转换成姿态欧拉角。本发明基于星敏感器偏转乘性噪声以及卫星姿态运动学方程，构建具有成型偏转误差的卫星姿态确定系统模型，通过对非线性卫星姿态确定系统模型进行二次线性展开，实现对星敏感器偏转乘性噪声的抑制和卫星在轨运动姿态的估计，可以有效降低乘性噪声对卫星姿态确定精度的影响。
Description
技术领域
本发明涉及卫星分析领域，尤其涉及一种提高星敏感器偏转乘性噪声处理能力的卫星姿态滤波方法。
背景技术
卫星姿态确定的主要任务是根据带有噪声的姿态敏感器的测量值估计卫星相对于某个参考坐标系的姿态参数。姿态确定系统主要由姿态敏感器和相应的信息处理算法即姿态确定算法组成，姿态确定精度取决于姿态敏感器测量精度和姿态确定算法精度。
卫星姿态确定系统是卫星姿态控制系统中的重要组成部分，高精度的姿态确定是高分辨率成像、高精度测绘等卫星应用的重要基础和关键技术。
星敏感器和陀螺的组合姿态确定方案以陀螺为基准，星敏感器对陀螺的累积误差和各种漂移进行校正，成为现代高精度卫星姿态确定的主要手段，在卫星的姿态测量方案中得到了广泛的实际应用。
利用星敏感器和陀螺两种姿态敏感器进行组合来确定卫星的在轨运动姿态，需要构建并求解由卫星姿态运动学模型和星敏感器姿态测量模型所共同构成的非线性系统模型。常见的姿态确定系统模型求解方法包括传统的扩展Kalman滤波(EKF)，但是EKF滤波算法仅适用于仅含加型噪声的非线性系统模型。
如果星敏感器相对于卫星本体的安装矩阵在卫星的运动过程中发生了偏转，或者由于抖动，振动等在轨环境的影响，使得星敏感器测量方程中引入另外一个未知的坐标变换。因此，星敏感器的测量随机噪声在卫星姿态确定系统的观测方程中往往体现为乘性噪声。
然而，由于卫星在轨运行环境复杂，使得星敏感器测量模型不仅包括热影响、像差影响、探测器噪声等引起的内部误差，而且还包括如相对安装误差、基准偏差以及形变误差等影响的外部误差，而这些外部误差的引入，使得卫星本体的安装矩阵在卫星的运动过程中发生了偏转，进而使得星敏感器测量模型中引入另外一个未知的坐标变换。因此，星敏感器的测量随机噪声在卫星姿态确定系统的观测方程中不仅体现为加性噪声，而且还体现为乘性噪声。因此，采用传统的EKF滤波算法，对乘性噪声的处理很难得到保证。
针对带乘性噪声的动态系统，目前也存在一些理论研究，如针对一维乘性噪声系统，P.K.Rajasekaram提出了一种处理乘性噪声的最优线性估计方法，给出了最优滤波算法，有效地解决了存在乘性噪声情况下的线性问题。褚东升等在此基础上，分别将线性系统的乘性噪声推广为对角阵和随机阵，F.Wang,V.Balakrishnan等人对乘性噪声的线性系统设计了鲁棒Kalman滤波算法，通过凸优化的方法确定滤波参数，并采用线性矩阵不等式技术求解凸优化问题，还给出了另外一种用于处理乘性噪声的鲁棒滤波算法，适合在线应用需求。
但针对带随机乘性噪声的非动态线性系统还没有人做过相关研究，且目前的研究也未涉及具体的应用层面，尤其是由星敏感器/陀螺组合的卫星姿态确定系统，这在某种程度上提供了星敏感器观测模型优化和卫星姿态确定精度提升的空间。
综上所述，星敏感器和陀螺组合的卫星姿态确定系统中，由于星敏感器观测模型中同时存在的加性噪声和乘性噪声，设计一种提高星敏感器偏转乘性噪声处理能力的卫星姿态确定算法，降低偏转乘性噪声对卫星姿态确定精度的影响，是提高星敏感器/陀螺组合姿态确定系统精度的技术难点之一。
发明内容
本发明要解决的技术问题是：如何提高对星敏感器偏转乘性噪声的处理能力，从而使卫星姿态确定精度提高。
本发明解决其技术问题的解决方案是：
一种卫星姿态滤波方法，包括以下步骤：
步骤1，根据星敏感器的偏转乘性误差，构建具有偏转乘性噪声的卫星姿态确定系统模型；
步骤2，对所述卫星姿态确定系统模型进行第一阶段状态滤波，得到第一状态滤波值以及状态的一步平滑值；
步骤3，对所述卫星姿态确定系统模型进行第二阶段状态滤波，得到第二状态滤波值；
步骤4，将第二状态滤波值的姿态四元数转换成姿态欧拉角，得到卫星的在轨姿态，即卫星姿态的偏航角、俯仰角以及滚动角。
作为上述技术方案的进一步改进，步骤1具体包括以下步骤：
步骤1.1，根据星敏感器的偏转乘性误差，得到具有偏转乘性噪声的星敏感器的观测方程；
步骤1.2，根据卫星姿态的四元数状态方程和观测方程，得到具有偏转乘性噪声的卫星姿态确定系统模型。
作为上述技术方案的进一步改进，步骤1.1中所述观测方程如式1所示：
z(t)＝m(t)h[x(t)]+v(t) 式1
其中，为偏转乘性噪声，δC(t)是星敏感器的乘性偏转随机阵，偏转的欧拉角误差分别定义为δψ、δφ和δθ，x(t)为状态变量， 是由惯性系到卫星本体系的四元数姿态矩阵，为姿态四元数， μI1＝[1 0 0]T，μI2＝[0 1 0]T，μI3＝[0 0 1]T分别为星敏感器不同的安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。
作为上述技术方案的进一步改进，步骤1.2中卫星姿态的四元数状态方程如式2所示：
其中，x(t)为状态变量，b(t)和d(t)分别为陀螺的常值漂移和相关漂移，ω(t)＝[ωx(t) ωy(t)ωz(t)]T表示星体坐标系相对于惯性坐标系的旋转角速度，矩阵矩阵Dτ为由相关时间常数τ构成的对角阵， 表示系统过程噪声，ng表示陀螺的测量噪声，nd和nb分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，w(t)的协方差为Q(t)。
作为上述技术方案的进一步改进，步骤1.2中具有偏转乘性噪声的卫星姿态确定系统模型如式3所示：
其中k表示离散化时刻，v(k)表示星敏感器的加性噪声，加性噪声v(k)的统计特性满足E{v(k)}＝0以及E{v(k)v(j)T}＝Rδkj，且{v(k)}、{w(k)}以及x(0)之间不相关；m(k)表示星敏感器的乘性噪声，乘性噪声m(k)的统计特性满足m(k)＝(mij(k))9×9，记则有将该等式记为Nij,gl(k)，i，j，g，l＝1，2......9，且{m(k)}与{w(k),v(k),x(0)}统计独立。
作为上述技术方案的进一步改进，步骤2具体包括以下步骤：
步骤2.1，将卫星姿态确定系统模型线性化，得到第一阶段线性化离散模型；
步骤2.2，对第一阶段线性化离散模型进行状态滤波，构建第一阶段状态滤波算法，得到第一状态滤波值以及状态的一步平滑值。
作为上述技术方案的进一步改进，步骤2.1中第一阶段线性化离散模型如式4所示：
其中B1(k)＝Γ(x(k|k),k)， x(k|k)以及x(k|k-1)分别为k时刻的状态滤波值和状态预测值。
作为上述技术方案的进一步改进，步骤2.2中第一阶段状态滤波算法包括以下步骤：
步骤2.2.1，计算状态变量的预测值
步骤2.2.2，计算状态变量的预测值的一步预测协方差阵，
步骤2.2.3，设置Kalman增益K1(k)，其中 p1ij为矩阵第i行第j列的元素；
， q1ij为矩阵第i行第j列的元素；
步骤2.2.4，计算第一状态滤波值 其中为测量的一步预测误差；
步骤2.2.5，计算一步平滑值
作为上述技术方案的进一步改进，步骤3具体包括以下步骤：
步骤3.1，将卫星姿态确定系统模型线性化，得到第二阶段线性化离散模型；
步骤3.2，对第二阶段线性化离散模型进行状态滤波，构建第二阶段状态滤波算法，得到第二状态滤波值。
作为上述技术方案的进一步改进，步骤3.1中第二阶段线性化离散模块如式5所示：
其中，B(k)＝Γ(x(k|k+1),k)， x(k|k+1)与x(k|k)分别为k时刻的状态平滑值和状态滤波值。
作为上述技术方案的进一步改进，步骤3.2具体包括以下步骤：
步骤3.2.1，计算状态预测值x(k/k-1)，x(k/k-1)＝f(x(k-1/k),k-1)；
步骤3.2.2，预测误差协方差阵P(k/k-1)，P(k/k-1)＝A(k,k-1)P(k-1)AT(k,k-1)+B(k-1)Q(k-1)BT(k-1)；
步骤3.2.3，设置Kalman增益K(k)，其中，
步骤3.2.4，计算第二状态滤波值 其中，为测量的一步预测误差，
步骤3.2.5，计算滤波误差协方差阵P(k)，
作为上述技术方案的进一步改进，步骤4中将第二状态滤波值的姿态四元数转换成姿态欧拉角，其中偏航角俯仰角滚动角其中而表示姿态四元数滤波值。
本发明的有益效果是：本发明基于星敏感器偏转乘性噪声以及卫星姿态运动学方程，构建具有成型偏转误差的星敏感器/陀螺组合姿态确定系统模型，通过对非线性卫星姿态确定系统模型进行二次线性展开，并基于线性最小方差准则设计卫星姿态确定滤波算法，实现对星敏感器偏转乘性噪声的抑制和卫星在轨运动姿态的估计，可以有效降低乘性噪声对卫星姿态确定精度的影响。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单说明。显然，所描述的附图只是本发明的一部分实施例，而不是全部实施例，本领域的技术人员在不付出创造性劳动的前提下，还可以根据这些附图获得其他设计方案和附图。
图1是本发明的方法流程图。
具体实施方式
以下将结合实施例和附图对本发明的构思、具体结构及产生的技术效果进行清楚、完整的描述，以充分地理解本发明的目的、特征和效果。显然，所描述的实施例只是本发明的一部分实施例，而不是全部实施例，基于本发明的实施例，本领域的技术人员在不付出创造性劳动的前提下所获得的其他实施例，均属于本发明保护的范围。
参照图1，本发明申请公开了一种卫星姿态滤波方法，包括以下步骤：
步骤1，根据星敏感器的偏转乘性误差，构建具有偏转乘性噪声的卫星姿态确定系统模型；
步骤2，对所述卫星姿态确定系统模型进行第一阶段状态滤波，得到第一状态滤波值以及状态的一步平滑值；
步骤3，对所述卫星姿态确定系统模型进行第二阶段状态滤波，得到第二状态滤波值；
步骤4，将第二状态滤波值的姿态四元数转换成姿态欧拉角，得到卫星的在轨姿态，即卫星姿态的偏航角、俯仰角以及滚动角。
具体地，本发明基于星敏感器偏转乘性噪声以及卫星姿态运动学方程，构建具有成型偏转误差的星敏感器/陀螺组合姿态确定系统模型，通过对非线性卫星姿态确定系统模型进行二次线性展开，并基于线性最小方差准则设计卫星姿态确定滤波算法，实现对星敏感器偏转乘性噪声的抑制和卫星在轨运动姿态的估计，可以有效降低乘性噪声对卫星姿态确定精度的影响。
进一步作为优选的实施方式，本发明创造具体实施方式中，步骤1具体包括以下步骤：
步骤1.1，根据星敏感器的偏转乘性误差，得到具有偏转乘性噪声的星敏感器的观测方程；
步骤1.2，根据卫星姿态的四元数状态方程和观测方程，得到具有偏转乘性噪声的卫星姿态确定系统模型。
进一步作为优选的实施方式，本发明创造具体实施方式中，步骤1.1中所述观测方程如式1所示：
z(t)＝m(t)h[x(t)]+v(t) 式1
其中，为偏转乘性噪声，δC(t)是星敏感器的乘性偏转随机阵，偏转的欧拉角误差分别定义为δψ、δφ和δθ，x(t)为状态变量， 是由惯性系到卫星本体系的四元数姿态矩阵，为姿态四元数， μI1＝[1 0 0]T，μI2＝[0 1 0]T，μI3＝[0 0 1]T分别为星敏感器不同的安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。
更进一步地，步骤1.2中卫星姿态的四元数状态方程如式2所示：
其中，x(t)为状态变量，b(t)和d(t)分别为陀螺的常值漂移和相关漂移，ω(t)＝[ωx(t) ωy(t)ωz(t)]T表示星体坐标系相对于惯性坐标系的旋转角速度，矩阵矩阵Dτ为由相关时间常数τ构成的对角阵， 表示系统过程噪声，ng表示陀螺的测量噪声，nd和nb分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，w(t)的协方差为Q(t)。
更进一步地，步骤1.2中具有偏转乘性噪声的卫星姿态确定系统模型如式3所示：
其中k表示离散化时刻，v(k)表示星敏感器的加性噪声，加性噪声v(k)的统计特性满足E{v(k)}＝0以及E{v(k)v(j)T}＝Rδkj，且{v(k)}、{w(k)}以及x(0)之间不相关；m(k)表示星敏感器的乘性噪声，乘性噪声m(k)的统计特性满足m(k)＝(mij(k))9×9，记则有将该等式记为Nij,gl(k)，i，j，g，l＝1，2......9，且{m(k)}与{w(k),v(k),x(0)}统计独立。
进一步作为优选的实施方式，本发明创造具体实施方式中，步骤2具体包括以下步骤：
步骤2.1，将卫星姿态确定系统模型线性化，得到第一阶段线性化离散模型；
步骤2.2，对第一阶段线性化离散模型进行状态滤波，构建第一阶段状态滤波算法，得到第一状态滤波值以及状态的一步平滑值。
作为上述技术方案的进一步改进，步骤2.1中第一阶段线性化离散模型如式4所示：
其中B1(k)＝Γ(x(k|k),k)， x(k|k)以及x(k|k-1)分别为k时刻的状态滤波值和状态预测值。
更进一步地，步骤2.2中第一阶段状态滤波算法包括以下步骤：
步骤2.2.1，计算状态变量的预测值
步骤2.2.2，计算状态变量的预测值的一步预测协方差阵，
步骤2.2.3，设置Kalman增益K1(k)，其中 p1ij为矩阵第i行第j列的元素；
， q1ij为矩阵第i行第j列的元素；
步骤2.2.4，计算第一状态滤波值 其中为测量的一步预测误差；
步骤2.2.5，计算一步平滑值
进一步作为优选的实施方式，本发明创造具体实施方式中，步骤3具体包括以下步骤：
步骤3.1，将卫星姿态确定系统模型线性化，得到第二阶段线性化离散模型；
步骤3.2，对第二阶段线性化离散模型进行状态滤波，构建第二阶段状态滤波算法，得到第二状态滤波值。
更进一步地，步骤3.1中第二阶段线性化离散模块如式5所示：
其中，B(k)＝Γ(x(k|k+1),k)， x(k|k+1)与x(k|k)分别为k时刻的状态平滑值和状态滤波值。
更进一步地，步骤3.2具体包括以下步骤：
步骤3.2.1，计算状态预测值x(k/k-1)，x(k/k-1)＝f(x(k-1/k),k-1)；
步骤3.2.2，预测误差协方差阵P(k/k-1)，P(k/k-1)＝A(k,k-1)P(k-1)AT(k,k-1)+B(k-1)Q(k-1)BT(k-1)；
步骤3.2.3，设置Kalman增益K(k)，其中，
步骤3.2.4，计算第二状态滤波值 其中，为测量的一步预测误差，
步骤3.2.5，计算滤波误差协方差阵P(k)，
进一步作为优选的实施方式，本发明创造具体实施方式中，步骤4中将第二状态滤波值的姿态四元数转换成姿态欧拉角，其中偏航角俯仰角滚动角其中而表示姿态四元数滤波值。
为更好地说明本技术方案的处理流程，下面结合应用实例对本发明申请进行详细说明。
步骤一，根据星敏感器的偏转乘性误差，构建带有偏转乘性噪声的卫星姿态确定系统模型，取转序为zxy，偏转的欧拉角误差分别为δψ(t)、δφ(t)和δθ(t)，则星敏感器的乘性偏转随机阵δC(t)为：
由于偏转的欧拉角误差δψ(t)、δφ(t)和δθ(t)均为小角度，可以认为cosδψ(t)＝cosδφ(t)＝cosδθ(t)＝1，sinδψ(t)＝δψ(t)，sinδφ(t)＝δφ(t)，sinδθ(t)＝δθ(t)，同时忽略二阶小量，可得：
δC(t)即为有偏转乘性误差引起的坐标变换，对三个星敏感器作此变换，得观测方程中偏转乘性噪声为结合星敏感器的测量原理，得到星敏感器的观测方程为z(t)＝m(t)h[x(t)]+v(t)，其中状态变量 是由惯性系到卫星本体系的四元数姿态矩阵，为姿态四元数， μI1＝[1 0 0]T，μI2＝[0 1 0]T，μI3＝[0 0 1]T分别为星敏感器不同的安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。
而卫星姿态的四元数状态方程为其中x(t)为状态变量，b(t)和d(t)分别为陀螺的常值漂移和相关漂移，ω(t)＝[ωx(t) ωy(t) ωz(t)]T表示星体坐标系相对于惯性坐标系的旋转角速度，矩阵矩阵Dτ为由相关时间常数τ构成的对角阵，表示系统过程噪声，ng表示陀螺的测量噪声，nd和nb分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，w(t)的协方差为Q(t)。
具有偏转乘性噪声的卫星姿态确定系统模型如式3所示：
其中k表示离散化时刻，v(k)表示星敏感器的加性噪声，加性噪声v(k)的统计特性满足E{v(k)}＝0以及E{v(k)v(j)T}＝Rδkj，且{v(k)}、{w(k)}以及x(0)之间不相关；m(k)表示星敏感器的乘性噪声，乘性噪声m(k)的统计特性满足m(k)＝(mij(k))9×9，记则有将该等式记为Nij,gl(k)，i，j，g，l＝1，2......9，且{m(k)}与{w(k),v(k),x(0)}统计独立。
步骤二，将卫星姿态确定系统模型线性化，得到第一阶段线性化离散模型；对第一阶段线性化离散模型进行状态滤波，构建第一阶段状态滤波算法，得到第一状态滤波值以及状态的一步平滑值。
第一阶段线性化离散模型如下所示：
其中B1(k)＝Γ(x(k|k),k)， x(k|k)以及x(k|k-1)分别为k时刻的状态滤波值和状态预测值。
步骤三，将卫星姿态确定系统模型线性化，得到第二阶段线性化离散模型；对第二阶段线性化离散模型进行状态滤波，构建第二阶段状态滤波算法，得到第二状态滤波值。
第二阶段线性化离散模块如下所示：
其中，B(k)＝Γ(x(k|k+1),k)， x(k|k+1)与x(k|k)分别为k时刻的状态平滑值和状态滤波值。
步骤四，在获得k时刻状态滤波值后，将第二状态滤波值的姿态四元数转换成姿态欧拉角，其中偏航角俯仰角滚动角
下面进行仿真计算以及结果分析。
步骤一中，取星敏感器偏转乘性噪声的均值为0，协方差阵为：
选取不同星敏感器的安装光轴矢量μI1＝[1 0 0]T，μI2＝[0 1 0]T，μI3＝[0 0 1]T，星敏感器的加性噪声均值为0，均方差为15角秒(3σ)，陀螺测量噪声ng＝0.0010/h，相关漂移噪声nd＝0.050/h，常值漂移噪声nb＝0.010/h。
步骤二中取采样间隔T＝1秒，个各分量的初始真值选取如下：b(0)＝[1 -1 1]0/h，d(0)＝[0.1 0.1 0.1]0/h；姿态四元数初始估值取陀螺的常值漂移和相关漂移的初始估值均为0，初始误差协方差阵为：
根据k时刻状态滤波值将第二状态滤波值的姿态四元数转换成姿态欧拉角。
表1给出了应用本技术方案与未应用本技术方案时欧拉角估计值的均方差。
| 均方差(3σ) | 滚动角(单位：角秒) | 俯仰角(单位：角秒) | 偏航角(单位：角秒) |
| 未应用本技术方案 | 15.781 | 17.156 | 16.735 |
| 应用本技术方案 | 9.363 | 8.269 | 8.754 |
表1
从表1中可以看出，本发明创造在星敏感器观测方程中引入偏转乘性噪声，通过二次滤波抑制其对卫星姿态确定精度的影响，在星敏感器测量精度为15角秒时，加入标准差为5角秒的乘性噪声后，采用本发明后可使欧拉角估计值的均方差提高一倍左右。
以上对本发明的较佳实施方式进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可作出种种的等同变型或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。
Claims (12)
1.一种卫星姿态滤波方法，其特征在于，包括以下步骤：
步骤1，根据星敏感器的偏转乘性误差，构建具有偏转乘性噪声的卫星姿态确定系统模型；
步骤2，对所述卫星姿态确定系统模型进行第一阶段状态滤波，得到第一状态滤波值以及状态的一步平滑值；
步骤3，对所述卫星姿态确定系统模型进行第二阶段状态滤波，得到第二状态滤波值；
步骤4，将第二状态滤波值的姿态四元数转换成姿态欧拉角，得到卫星的在轨姿态。
2.根据权利要求1所述的一种卫星姿态滤波方法，其特征在于，步骤1具体包括以下步骤：
步骤1.1，根据星敏感器的偏转乘性误差，得到具有偏转乘性噪声的星敏感器的观测方程；
步骤1.2，根据卫星姿态的四元数状态方程和观测方程，得到具有偏转乘性噪声的卫星姿态确定系统模型。
3.根据权利要求2所述的一种卫星姿态滤波方法，其特征在于，步骤1.1中所述观测方程如式1所示：
z(t)＝m(t)h[x(t)]+v(t) 式1
其中，为偏转乘性噪声，δC(t)是星敏感器的乘性偏转随机阵，偏转的欧拉角误差分别定义为δψ、δφ和δθ，x(t)为状态变量， 是由惯性系到卫星本体系的四元数姿态矩阵，为姿态四元数， μI1＝[1 0 0]T，μI2＝[0 1 0]T，μI3＝[0 0 1]T分别为不同的星敏感器安装光轴矢量，为星敏感器的加性测量噪声，加性测量噪声的协方差为R(t)，t为时间参数。
4.根据权利要求3所述的一种卫星姿态滤波方法，其特征在于，步骤1.2中卫星姿态的四元数状态方程如式2所示：
其中，x(t)为状态变量，b(t)和d(t)分别为陀螺的常值漂移和相关漂移，ω(t)＝[ωx(t)ωy(t)ωz(t)]T表示星体坐标系相对于惯性坐标系的旋转角速度，矩阵矩阵Dτ为由相关时间常数τ构成的对角阵， 表示系统过程噪声，ng表示陀螺的测量噪声，nd和nb分别为陀螺的相关漂移噪声以及常值漂移噪声，I和0分别表示相应维数的单位矩阵和零矩阵，w(t)的协方差为Q(t)。
5.根据权利要求4所述的一种卫星姿态滤波方法，其特征在于，步骤1.2中具有偏转乘性噪声的卫星姿态确定系统模型如式3所示：
其中k表示离散化时刻，v(k)表示星敏感器的加性噪声，加性噪声v(k)的统计特性满足E{v(k)}＝0以及E{v(k)v(j)T}＝Rδkj，且{v(k)}、{w(k)}以及x(0)之间不相关；m(k)表示星敏感器的乘性噪声，乘性噪声m(k)的统计特性满足m(k)＝(mij(k))9×9，记则有将该等式记为Nij,gl(k)，i，j，g，l＝1，2......9，且{m(k)}与{w(k),v(k),x(0)}统计独立。
6.根据权利要求5所述的一种卫星姿态滤波方法，其特征在于，步骤2具体包括以下步骤：
步骤2.1，将卫星姿态确定系统模型线性化，得到第一阶段线性化离散模型；
步骤2.2，对第一阶段线性化离散模型进行状态滤波，构建第一阶段状态滤波算法，得到第一状态滤波值以及状态的一步平滑值。
7.根据权利要求6所述的一种卫星姿态滤波方法，其特征在于，步骤2.1中第一阶段线性化离散模型如式4所示：
其中B1(k)＝Γ(x(k|k),k)， x(k|k)以及x(k|k-1)分别为k时刻的状态滤波值和状态预测值。
8.根据权利要求7所述的一种卫星姿态滤波方法，其特征在于，步骤2.2中第一阶段状态滤波算法包括以下步骤：
步骤2.2.1，计算状态变量的预测值
步骤2.2.2，计算状态变量的预测值的一步预测协方差阵，
步骤2.2.3，设置Kalman增益K1(k)，其中 p1ij为矩阵第i行第j列的元素；
， q1ij为矩阵第i行第j列的元素；
步骤2.2.4，计算第一状态滤波值 其中为测量的一步预测误差；
步骤2.2.5，计算一步平滑值
9.根据权利要求8所述的一种卫星姿态滤波方法，其特征在于，步骤3具体包括以下步骤：
步骤3.1，将卫星姿态确定系统模型线性化，得到第二阶段线性化离散模型；
步骤3.2，对第二阶段线性化离散模型进行状态滤波，构建第二阶段状态滤波算法，得到第二状态滤波值。
10.根据权利要求9所述的一种卫星姿态滤波方法，其特征在于，步骤3.1中第二阶段线性化离散模块如式5所示：
其中，B(k)＝Γ(x(k|k+1),k)， x(k|k+1)与x(k|k)分别为k时刻的状态平滑值和状态滤波值。
11.根据权利要求10所述的一种卫星姿态滤波方法，其特征在于，步骤3.2具体包括以下步骤：
步骤3.2.1，计算状态预测值x(k/k-1)，x(k/k-1)＝f(x(k-1/k),k-1)；
步骤3.2.2，预测误差协方差阵P(k/k-1)，P(k/k-1)＝A(k,k-1)P(k-1)AT(k,k-1)+B(k-1)Q(k-1)BT(k-1)；
步骤3.2.3，设置Kalman增益K(k)，其中，
步骤3.2.4，计算第二状态滤波值 其中，为测量的一步预测误差，
步骤3.2.5，计算滤波误差协方差阵P(k)，
12.根据权利要求11所述的一种卫星姿态滤波方法，其特征在于，步骤4中将第二状态滤波值的姿态四元数转换成姿态欧拉角，其中偏航角俯仰角滚动角其中而表示姿态四元数滤波值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811300220.0A CN109470267A (zh) | 2018-11-02 | 2018-11-02 | 一种卫星姿态滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811300220.0A CN109470267A (zh) | 2018-11-02 | 2018-11-02 | 一种卫星姿态滤波方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109470267A true CN109470267A (zh) | 2019-03-15 |
Family
ID=65666893
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811300220.0A Pending CN109470267A (zh) | 2018-11-02 | 2018-11-02 | 一种卫星姿态滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109470267A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110411477A (zh) * | 2019-08-06 | 2019-11-05 | 广州泾渭信息科技有限公司 | 基于序列机动的星敏安装误差在轨标定方法 |
| CN111444475A (zh) * | 2020-03-24 | 2020-07-24 | 广东海洋大学深圳研究院 | 一种应用在飞行测试数据分析的容错ckf滤波融合方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2012240433A (ja) * | 2011-05-16 | 2012-12-10 | Mitsubishi Electric Corp | 人工衛星の姿勢決定装置および人工衛星の姿勢決定方法 |
| CN105300384A (zh) * | 2015-04-03 | 2016-02-03 | 东南大学 | 一种用于卫星姿态确定的交互式滤波方法 |
| CN108333922A (zh) * | 2017-12-26 | 2018-07-27 | 佛山科学技术学院 | 一种基于智能优化与约束推理的单星自主任务规划方法 |
-
2018
- 2018-11-02 CN CN201811300220.0A patent/CN109470267A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2012240433A (ja) * | 2011-05-16 | 2012-12-10 | Mitsubishi Electric Corp | 人工衛星の姿勢決定装置および人工衛星の姿勢決定方法 |
| CN105300384A (zh) * | 2015-04-03 | 2016-02-03 | 东南大学 | 一种用于卫星姿态确定的交互式滤波方法 |
| CN108333922A (zh) * | 2017-12-26 | 2018-07-27 | 佛山科学技术学院 | 一种基于智能优化与约束推理的单星自主任务规划方法 |
Non-Patent Citations (5)
| Title |
|---|
| 杨喆等: "《弹目姿轨复合交会精准起爆控制》", 30 June 2016, 国防工业出版社 * |
| 江丹等: "几何扩展卡尔曼滤波算法在卫星姿态确定系统中的应用", 《传感器与微系统》 * |
| 王炯琦等: "适合处理乘性噪声估计卫星姿态的非线性迭代滤波算法", 《电子学报》 * |
| 王飞: "基于乘性扩展卡尔曼滤波的微小卫星姿态测量算", 《电子设计工程》 * |
| 矫媛媛: "基于 MEKF 的卫星姿态确定精度影响因素分析", 《系统工程与电子》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110411477A (zh) * | 2019-08-06 | 2019-11-05 | 广州泾渭信息科技有限公司 | 基于序列机动的星敏安装误差在轨标定方法 |
| CN111444475A (zh) * | 2020-03-24 | 2020-07-24 | 广东海洋大学深圳研究院 | 一种应用在飞行测试数据分析的容错ckf滤波融合方法 |
| CN111444475B (zh) * | 2020-03-24 | 2023-07-14 | 广东海洋大学深圳研究院 | 一种应用在飞行测试数据分析的容错ckf滤波融合方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109001787B (zh) | 一种姿态角解算与定位的方法及其融合传感器 | |
| Wu et al. | Generalized linear quaternion complementary filter for attitude estimation from multisensor observations: An optimization approach | |
| US9417091B2 (en) | System and method for determining and correcting field sensors errors | |
| Vasconcelos et al. | Geometric approach to strapdown magnetometer calibration in sensor frame | |
| Pylvänäinen | Automatic and adaptive calibration of 3D field sensors | |
| CN109470266A (zh) | 一种处理乘性噪声的星敏感器陀螺组合定姿方法 | |
| CN110109470A (zh) | 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制系统 | |
| Schopp et al. | Self-calibration of accelerometer arrays | |
| KR20130143576A (ko) | 중력 기준계에서의 장치에 부착된 움직임 센서 및 자력계의 측정치를 사용하여 장치의 요각을 추정하는 장치 및 방법 | |
| WO2004029549A2 (en) | Method and system for processing pulse signals within an inertial navigation system | |
| CN113465599B (zh) | 定位定向方法、装置及系统 | |
| CN108344409B (zh) | 提高卫星姿态确定精度的方法 | |
| Ludwig et al. | Comparison of attitude and heading reference systems using foot mounted MIMU sensor data: Basic, Madgwick, and Mahony | |
| CN108332775A (zh) | 提高星敏感器姿态测量精度的方法 | |
| CN114754798A (zh) | 一种陀螺误差特性参数在轨辨识与标定方法 | |
| CN111189474A (zh) | 基于mems的marg传感器的自主校准方法 | |
| CN109470267A (zh) | 一种卫星姿态滤波方法 | |
| CN107860382B (zh) | 一种在地磁异常情况下应用ahrs测量姿态的方法 | |
| Garcia et al. | Unscented Kalman filter for spacecraft attitude estimation using quaternions and euler angles | |
| CN103954288A (zh) | 一种卫星姿态确定系统精度响应关系确定方法 | |
| CN110375773B (zh) | Mems惯导系统姿态初始化方法 | |
| CN110702142A (zh) | 一种采用三轴加速度计辅助的三轴磁强计全参数外场标定方法 | |
| CN114858166B (zh) | 基于最大相关熵卡尔曼滤波器的imu姿态解算方法 | |
| CN110672127A (zh) | 阵列式mems磁传感器实时标定方法 | |
| CN114964259A (zh) | 基于红外焦平面阵列时域信号与红外图像的姿态解算方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190315 |