CN104318119B - 一种高动态下星点质心误差补偿方法 - Google Patents

Abstract

本发明公布了一种高动态下星点质心误差补偿方法，包括：A.三轴角速度的计算，对三轴角速度分别进行最小二乘法拟合，并对拟合结果进行卡尔曼滤波；B.将曝光时间内三轴角速度近似用一阶多项式表示，将其积分得到三轴旋转的角度信息，代入星敏感器的旋转矩阵；C.对旋转矩阵进行二阶泰勒展开近似，经透视投影变换，得到星点光斑在图像传感器上的移动轨迹坐标表达式；D.对星点轨迹坐标表达式进行二阶泰勒展开近似，参考质心法公式求解星点质心坐标，计算出星点质心误差补偿量。本发明提高了星敏感器在高动态下星点定位精度，处理算法较为简单，补偿量计算精确，在高动态下有较好的应用前景。

Description

一种高动态下星点质心误差补偿方法

技术领域

本发明涉及航空航天领域星敏感器技术，尤其涉及一种高动态下星点质心误差补偿方法。

背景技术

星敏感器已经是当今航天飞行器中广泛采用的一种高精度、高可靠性的姿态测量装置。其工作原理为：首先利用高清镜头拍摄周天恒星图像，然后提取图像中全部恒星像星点并利用质心法得到其精确位置信息，之后进行星图识别得到星点在周天中对应的实际恒星，最后参考导航星表计算出得到星敏感器的三轴姿态。星点提取时星点的质心定位精度是星敏感器的重要精度指标，将直接影响三轴姿态的计算精度。

星敏感器一般工作在静态情况下，即星敏感器载体运动较为缓慢，在曝光时恒星和星敏感器几乎是相对静止，恒星只在图像传感器的固定位置感光形成一个理想的弥散点状光斑，这种情况下利用传统质心法可以得到很高的星点定位精度。但随着星敏感器的应用领域的不断扩展，尤其是在高动态载体上(如敏捷卫星、弹道导弹等)的应用，对星敏感器动态性能指标提出了更高的要求。

何谓高动态，即星敏感器载体角速度较大，且变化比较复杂，使得星敏感器在曝光时有明显的旋转运动，从而使成像星点是一条长的光斑移动轨迹，这会造成星点质心定位的精度明显降低。

目前提高高动态下星点定位精度的方法主要有三种：(1)通过改变CCD图像传感器内部时序设计和数据处理算法，使得恒星产生的电荷在CCD成像区跟随恒星同步移动，达到恒星始终在固定像元上感光成像的效果，以避免光斑轨迹的产生。但这种算法需要合理设计电荷转移次数，对于积分时间过长、星点运动轨迹不规则的情况，它的设计实现比较困难。(2)通过对高动态下星敏感器的主要运动形式建立星图退化模型，将多重模糊星图进行分步复原，以还原星点移动图像。但这种去模糊算法都需要已知星敏感器的运动函数，而在通常实际应用中这些信息都是未知的，且星图复原计算量很大。(3)利用高数据更新率的惯性器件(如陀螺)作为主姿态敏感器提供输出姿态，在考虑由于星像移动引起的运动误差的基础上，采用星敏感器作为附姿态敏感器进行姿态修正，以满足高动态下的姿态精度要求。该方法依赖于惯性器件，数学模型复杂，不利于参数求解。

可以看出，在提高高动态星点的定位精度方面目前亟需一种能适应星点各种运动轨迹，无需已知星敏感器的运动函数，且不依赖惯性器件，计算简单，效果较好的方法。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种高动态下星点定位误差补偿方法，解决高动态下星点定位精度较低的问题，其不仅能适应各种星点运动轨迹，且无需已知星敏感器运动函数，计算简单，不依赖惯性器件，能在高动态下有较好的应用前景。

为达到上述目的，本发明的技术方案是这样实现的：

一种高动态下星点质心误差补偿方法，包括以下步骤：

步骤201，三轴角速度的计算，即对三轴角速度分别进行最小二乘法拟合，并对拟合结果进行卡尔曼滤波；

步骤202，将曝光时间内三轴角速度近似用一阶多项式表示，将一阶多项式积分得到三轴旋转的角度信息，代入星敏感器的旋转矩阵；

步骤203，对星敏感器的旋转矩阵进行二阶泰勒展开近似，经透视投影变换，得到星点光斑在图像传感器上的移动轨迹坐标表达式；

步骤204，对所述移动轨迹坐标表达式进行二阶泰勒展开近似，参考质心法公式求解星点质心坐标，计算出星点质心误差补偿量。

所述步骤201中的卡尔曼滤波具体为：设角速度的拟合时间多项式的系数为系统状态向量，以角速度的最小二乘拟合公式建立量测模型，以多项式时间零点的转移建立系统预测模型，对三轴角速度拟合分别进行滤波。

所述步骤202中代入星敏感器的旋转矩阵具体为：

三轴角速度用一阶多项式表示后，将其积分可得到曝光时间内三轴旋转的角度信息，将星点光斑坐标表示成星敏感器坐标系下的单位方向矢量，将所述角度信息代入旋转矩阵，从而建立星点光斑方向矢量随时间的运动方程。

所述步骤203中对旋转矩阵的二阶泰勒展开的具体为：

由于旋转矩阵的复杂和非线性，对旋转矩阵每个元素进行二阶泰勒展开近似，将旋转矩阵每个元素化为角速度参数的二阶时间多项式，这样就将星点光斑的方向矢量用时间二阶多项式表示。

所述步骤204中对星点轨迹坐标表达式进行二阶泰勒展开近似的具体为：

经透视投影变换后，星点轨迹坐标表达式分子分母均为二阶时间多项式，对其进行二阶泰勒展开近似，化为简单的二阶时间多项式，系数均由角速度参数表示，计算出星点质心和星点理想中心，将星点质心和星点理想中心相减得出星点质心误差补偿量。

本发明与现有技术相比的优点在于：本发明提高了星敏感器在高动态下星点定位精度，处理算法较为简单，补偿量计算精确，其不仅能适应各种星点运动轨迹，且无需已知星敏感器运动函数，在高动态下有较好的应用前景。

附图说明

图1为本发明高动态下典型的光斑移动轨迹；

图2为本发明高动态下星点定位误差产生示意图；

图3为本发明方法实现流程图；

图4为本发明小角度机动下姿态四元数和角速度的变化曲线；

图5为本发明小角度机动下未补偿与补偿后误差比较图；

图6为本发明大角度机动下姿态四元数和角速度的变化曲线；

图7为本发明大角度机动下未补偿与补偿后误差比较图；

图8为本发明两种机动情况下用起始光斑坐标和质心坐标补偿效果对比图。

具体实施方式

为了更好的理解本发明，首先介绍一下高动态星点质心误差的产生原因。如图1所示，为高动态下典型的星点光斑移动轨迹。星敏感器在静态情况下其星点是近似符合二维高斯分布的理想弥散光斑。而在高动态下，由于星敏感器在曝光时有明显运动，从而使得恒星星点是一条长的光斑移动轨迹。根据透视投影模型，光斑移动轨迹一般不会是均匀的完美直线，特别是当角速度变化复杂时，光斑移动轨迹会有较大的弯曲度。

如图2所示，为高动态下星点质心误差产生示意图。高动态情况下，通常定义曝光时间中间时刻的光斑位置为整个星点的中心，如图2中A点。设曝光时间内光斑的移动曲线坐标为(X(t)，Y(t))，根据质心法公式可得质心

$\begin{array}{c}\stackrel{\‾}{x}=\frac{\underset{C}{\∫\∫}\mathrm{xS}(x,y)\mathrm{dxdy}}{\underset{C}{\∫\∫}S(x,y)\mathrm{dxdy}}=\frac{{\∫}_0^{T_e}X\left(t\right)\mathrm{dt}}{T_e}\\ \stackrel{\‾}{y}=\frac{\underset{C}{\∫\∫}\mathrm{yS}(x,y)\mathrm{dxdy}}{\underset{C}{\∫\∫}S(x,y)\mathrm{dxdy}}=\frac{{\∫}_0^{T_e}Y\left(t\right)\mathrm{dt}}{T_e}\end{array}---\left(1\right)$

其中，T_e为曝光时间，S(x，y)为高动态星点能量分布函数，C为整个高动态星点的质心法积分窗口。

得到的质心是光斑移动轨迹在曝光时间内的积和曝光时间的比值，并不是(X(T_e/2)，Y(T_e/2))，如图2中B点，两者之间存在一定的误差，这种误差主要由星点光斑移动曲线坐标(X(t)，Y(t))决定。而(X(t)，Y(t))是由于星敏感器在曝光时间内相对于恒星的旋转运动形成的，所以称其为星点质心的运动误差，其主要由星敏感器三轴旋转角速度决定。

下面结合附图和具体仿真实例对本发明的技术方案进一步详细阐述。如图3所示，

步骤201，三轴角速度的计算，对所述三轴角速度分别进行最小二乘法拟合，并对拟合结果进行卡尔曼滤波；

三轴角速度的计算方法可参见Crassidis.J.L的“Angular Velocity Determination DirectlyFrom Star Tracker Measurements.Journal of Guidance，Control，and Dynamics，2002，25(6)：1165～1168”。

补偿量计算中需要当前帧曝光时刻的角速度信息，而角速度计算只能计算出上一帧的曝光中间时刻的瞬时角速度，故当前帧角速度信息需要进行预测。三轴角速度得分别进行拟合，以下以x轴角速度为例。设定角速度拟合的多项式为三阶，表示为ω_x＝a+bt+ct²+dt³，拟合的时间零点为当前帧曝光起始时刻，设卡尔曼滤波的系统状态向量为X(n)＝[a(n)  b(n)  c(n)  d(n)]。每一帧系统状态的预测可由上一帧状态向量转移时间零点得到，由此推导的系统预测模型为：

$X(n+1)=A\left(n\right)X\left(n\right)+W\left(n\right)=\left[\begin{array}{cccc}1& \mathrm{\Δt}& {\mathrm{\Δt}}^2& {\mathrm{\Δt}}^3\\ 0& 1& 2\mathrm{\Δt}& 3\mathrm{\Δ}{t}^2\\ 0& 0& 1& 3\mathrm{\Δt}\\ 0& 0& 0& 1\end{array}\right]X\left(n\right)+W\left(n\right)---\left(2\right)$

其中，W(n)为系统过程噪声，Δt为星敏感器每一帧的周期。

设观测向量为Z(n)＝[ω_x(n)  ω_x(n-1)  ω_x(n-2)  ω_x(n-3)]，由多项式拟合公式可得量测模型：

$F=\left[\begin{array}{ccc}1& -\mathrm{\Δt}+T_e/2& {(-\mathrm{\Δt}+T_e/2)}^2{(-\mathrm{\Δt}+T_e/2)}^3\\ 1& -2\mathrm{\Δt}+T_e/2& {(-2\mathrm{\Δt}+T_e/2)}^2{(-2\mathrm{\Δt}+T_e/2)}^3\\ 1& -3\mathrm{\Δt}+T_e/2& {(-3\mathrm{\Δt}+T_e/2)}^2{(-3\mathrm{\Δt}+T_e/2)}^3\\ 1& -4\mathrm{\Δt}+T_e/2& {(-4\mathrm{\Δt}+T_e/2)}^2{(-4\mathrm{\Δt}+T_e/2)}^3\end{array}\right]---\left(3\right)$

Z(n)＝H(n)X(n)+V(n)＝(F^T)^-1(F^TF)X(n)+V(n)

其中，T_e为每帧的曝光时间，V(n)为观测噪声。系统过程噪声W(n)和测量噪声V(n)均为相互独立的零均值的高斯白噪声。

这样就把一个最小二乘法拟合问题转化成卡尔曼滤波问题，可得到最小均方误差意义下的最优拟合参数。由于补偿量计算只需将角速度表示成一阶多项式，那么只需要拟合中a，b参数，其分别为当前帧的角速度和角加速度。

步骤202，将曝光时间内三轴角速度近似用一阶多项式表示，将其积分得到三轴旋转的总角度，代入星敏感器的旋转矩阵；

星点质心的运动误差大小主要由星敏感器在曝光时三轴旋转角速度决定，但由于星敏感器载体运动的复杂多样，这角速度并不能用一简单函数表示。采用近似方法，假设曝光时星敏感器的三轴角加速度固定，那么角速度的时间表达式可以用简单一阶多项式表示，一般高动态下曝光时间都比较短，这种一阶多项式近似可以达到较高精度。

星敏感器三轴各在曝光时旋转的角度信息可用各自角速度随时间t积分表示：

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_x\left(t\right)\&ap;a_x+b_{x}t\\ {\mathrm{\&omega;}}_y\left(t\right)\&ap;a_y+b_{y}t\\ {\mathrm{\&omega;}}_z\left(t\right)\&ap;a_z+b_{z}t\end{array}\right.\&RightArrow;\left\{\begin{array}{c}\mathrm{\&Delta;x}\left(t\right)\&ap;{\&Integral;}_0^t{\mathrm{\&omega;}}_x\left(t\right)=a_{x}t+b_x{t}^2/2\\ \mathrm{\&Delta;y}\left(t\right)\&ap;{\&Integral;}_0^t{\mathrm{\&omega;}}_y\left(t\right)=a_{y}t+b_y{t}^2/2(0<t<T_e)\\ \mathrm{\&Delta;z}\left(t\right)\&ap;{\&Integral;}_0^t{\mathrm{\&omega;}}_z\left(t\right)=a_{z}t+b_z{t}^2/2\end{array}\right.---\left(4\right)$

其中，(Δx(t)，Δy(t)，Δz(t))为星敏感器三轴各在曝光时旋转的角度信息，(ω_x(t)，ω_y(t)，ω_z(t))为曝光时间内的三轴角速度，(a_x，a_y，a_z)为曝光起始时刻的三个轴的角速度，(b_x，b_y，b_z)为曝光时间内近似固定的三个轴的角加速度。

设曝光起始时刻的星点光斑在星敏感器坐标系上的单位方向矢量为(x₀，y₀，z₀)。由于星点光斑在曝光时间内是移动的，设其曝光到t时星点光斑的单位方向矢量为(x(t)，y(t)，z(t))。故依据3-1-2坐标轴旋转矩阵可得：

$\left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ z\left(t\right)\end{array}\right]=\left[\begin{array}{ccc}\cos \mathrm{\&Delta;y}\left(t\right)\cos \mathrm{\&Delta;z}\left(t\right)-\sin \mathrm{\&Delta;x}\left(t\right)\sin \mathrm{\&Delta;y}\left(t\right)\sin \mathrm{\&Delta;z}\left(t\right)& \cos \mathrm{\&Delta;y}\left(t\right)\sin \mathrm{\&Delta;z}\left(t\right)+\sin \mathrm{\&Delta;x}\left(t\right)\sin \mathrm{\&Delta;y}\left(t\right)\cos \mathrm{\&Delta;z}\left(t\right)& -\cos \mathrm{\&Delta;x}\left(t\right)\sin \mathrm{\&Delta;y}\left(t\right)\\ -\cos \mathrm{\&Delta;x}\left(t\right)\sin \mathrm{\&Delta;z}\left(t\right)& \cos \mathrm{\&Delta;x}\left(t\right)\cos \mathrm{\&Delta;z}\left(t\right)& \sin \mathrm{\&Delta;x}\left(t\right)\\ \sin \mathrm{\&Delta;y}\left(t\right)\cos \mathrm{\&Delta;z}\left(t\right)+\sin \mathrm{\&Delta;x}\left(t\right)\cos \mathrm{\&Delta;y}\left(t\right)\sin \mathrm{\&Delta;z}\left(t\right)& \sin \mathrm{\&Delta;y}\left(t\right)\sin \mathrm{\&Delta;z}\left(t\right)-\sin \mathrm{\&Delta;x}\left(t\right)\cos \mathrm{\&Delta;y}\left(t\right)\cos \mathrm{\&Delta;z}\left(t\right)& \cos \mathrm{\&Delta;x}\left(t\right)\cos \mathrm{\&Delta;y}\left(t\right)\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]---\left(5\right)$

$\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]=\frac{1}{\sqrt{X_0^2+Y_0^2+f^2}}\left[\begin{array}{c}{X}_0\\ Y_0\\ f\end{array}\right]---\left(6\right)$

其中，(Δx(t)，Δy(t)，Δz(t))为上述星敏感器三轴各在曝光时旋转的角度信息，(X₀，Y₀)为起始时刻星点光斑坐标，f为镜头焦距。

步骤203，对所述旋转矩阵进行二阶泰勒展开近似，经透视投影变换，得到星点在图像传感器上的光斑移动轨迹坐标表达式；

对旋转矩阵每个元素进行二阶的泰勒展开，近似可得曝光到t时星点光斑的单位方向矢量(x(t)，y(t)，z(t))：

$\begin{array}{c}\left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ z\left(t\right)\end{array}\right]\&ap;\left[\begin{array}{ccc}1-(a_y^2+a_z^2)t^2/2& a_{z}t+(a_x{a}_y+b_z/2)t^2& -a_{y}t-b_y{t}^2/2\\ -a_{z}t-b_z{t}^2/2& 1-(a_x^2+a_z^2)t^2/2& a_{x}t+b_x{t}^2/2\\ a_{y}t+(a_x{a}_z+b_y/2)t^2& -a_{x}t+(a_y{a}_z-b_x/2)t^2& 1-(a_x^2+a_y^2)t^2/2\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]\\ \left[\begin{array}{c}x\left(t\right)\\ y\left(t\right)\\ z\left(t\right)\end{array}\right]\&ap;\left[\begin{array}{ccc}-(a_y^2+a_z^2)/2& a_x{a}_y+b_z/2& -b_y/2\\ -b_z/2& -(a_x^2+a_z^2)/2& b_x/2\\ a_x{a}_z+b_y/2& a_y{a}_z-b_x/2& -(a_y^2+a_y^2)/2\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]t^2+\left[\begin{array}{ccc}1& a_z& -a_y\\ -a_z& 1& a_x\\ a_y& -a_x& 1\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]t+\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right]\end{array}---\left(7\right)$

其中，(a_x，a_y，a_z)为曝光起始时刻的三个轴的角速度，(b_x，b_y，b_z)为曝光时间内近似固定的三个轴的角加速度，(x₀，y₀，z₀)为曝光起始时刻的星点光斑在星敏感器坐标系上的单位方向矢量。

根据星敏感器透视投影模型，则星点光斑在曝光时间内移动轨迹坐标(X(t)，Y(t))为：

$X\left(t\right)=f\frac{x\left(t\right)}{z\left(t\right)},Y\left(t\right)=f\frac{y\left(t\right)}{z\left(t\right)}---\left(8\right)$

令

$X\left(t\right)=f\frac{x\left(t\right)}{x\left(t\right)}\&ap;f\frac{k_1{t}^2+k_{2}t+x_0}{k_5{t}^2+k_{6}t+f},Y\left(t\right)\&ap;f\frac{y\left(t\right)}{z\left(t\right)}=f\frac{k_3{t}^2+k_{4}t+y_0}{k_5{t}^2+k_{6}t+f}---\left(9\right)$

则

$\begin{array}{c}\left[\begin{array}{c}{k}_1\\ k_3\\ k_5\end{array}\right]=\left[\begin{array}{ccc}-(a_y^2+a_z^2)/2& a_x{a}_y+b_z/2& -b_y/2\\ -b_z/2& -(a_x^2+a_z^2)/2& b_x/2\\ a_x{a}_z\text{+}{b}_y/2& a_y{a}_z-b_x/2& -(a_x^2+a_y^2)/2\end{array}\right]\left[\begin{array}{c}{X}_0\\ Y_0\\ f\end{array}\right]\\ \left[\begin{array}{c}{k}_2\\ k_4\\ k_6\end{array}\right]=\left[\begin{array}{ccc}1& a_z& -a_y\\ -a_z& 1& a_x\\ a_y& -a_x& 1\end{array}\right]\left[\begin{array}{c}{X}_0\\ Y_0\\ f\end{array}\right]\end{array}---\left(10\right)$

其中，(k₁，k₂)为x(t)的二阶系数和一阶系数，(k₃，k₄)为y(t)的二阶系数和一阶系数，(k₅，k₆)为z(t)的二阶系数和一阶系数，(a_x，a_y，a_z)为上述曝光起始时刻的角速度，(b_x，b_y，b_z)为上述曝光时间内近似固定角加速度，(X₀，Y₀)为起始时刻星点光斑坐标，f为镜头焦距。

步骤204，对所述星点光斑轨迹坐标表达式(X(t)，Y(t))进行二阶泰勒展开近似，参考质心法公式求解星点质心坐标，计算出星点质心运动误差补偿量。

对星点光斑轨迹坐标表达式X(t)，Y(t)继续进行二阶泰勒展开近似得到：

$\begin{array}{c}X\left(t\right)\&ap;[(k_1-\frac{k_2{k}_6+k_5{X}_0}{f}+\frac{k_6^2{X}_0}{f^2})t^2+(k_2-\frac{k_6{X}_0}{f})t+X_0]\\ Y\left(t\right)\&ap;[(k_3-\frac{k_{4}g+k_5{Y}_0}{f}+\frac{k_6^2{Y}_0}{f^2})t^2+(k_4-\frac{k_6{Y}_0}{f})t+Y_0]\end{array}---\left(10\right)$

依据质心法公式可得：

$\begin{array}{c}\stackrel{\&OverBar;}{x}=\frac{{\&Integral;}_0^{T_e}X\left(t\right)\mathrm{dt}}{T_e}\&ap;[(k_1-\frac{k_2{k}_6+k_5{X}_0}{f}+\frac{k_6^2{X}_0}{f^2})\frac{T_e^2}{3}+(k_2-\frac{k_6{X}_0}{f})\frac{T_e}{2}+X_0]\\ \stackrel{\&OverBar;}{y}=\frac{{\&Integral;}_0^{T_e}Y\left(t\right)\mathrm{dt}}{T_e}\&ap;[(k_3-\frac{k_{4}g+k_5{Y}_0}{f}+\frac{k_6^2{Y}_0}{f^2})\frac{T_e^2}{3}+(k_4-\frac{k_6{Y}_0}{f})\frac{T_e}{2}{Y}_0]\end{array}---\left(11\right)$

其中，为高动态星点的质心，T_e为曝光时间，(k₁，k₂)为x(t)的二阶系数和一阶系数，(k₃，k₄)为y(t)的二阶系数和一阶系数，(k₅，k₆)为z(t)的二阶系数和一阶系数，(X₀，Y₀)为起始时刻星点光斑坐标，f为镜头焦距。

高动态下星点的真实中心是曝光中间时刻的星点光斑轨迹坐标，即(X(T_e/2)，Y(T_e/2))，故星点质心的运动误差补偿量(δ_x，δ_y)为：

$\begin{array}{c}{\mathrm{\&delta;}}_x=X(T_e/2)-\stackrel{\&OverBar;}{x}\&ap;-\frac{T_e^2}{12}(k_1-\frac{k_2{k}_6+k_5{X}_0}{f}+\frac{k_6^2{X}_0}{f^2})\\ {\mathrm{\&delta;}}_y=Y(T_e/2)-\stackrel{\&OverBar;}{y}\&ap;-\frac{T_e^2}{12}(k_3-\frac{k_{4}g+k_5{Y}_0}{f}+\frac{k_6^2{Y}_0}{f^2})\end{array}---\left(12\right)$

其中，(k₁，k₂)为x(t)的二阶系数和一阶系数，(k₃，k₄)为y(t)的二阶系数和一阶系数，(k₅，k₆)为z(t)的二阶系数和一阶系数，(X₀，Y₀)为起始时刻星点光斑坐标，f为镜头焦距，T_e为曝光时间。

将角速度拟合得到的a，b及曝光起始时刻光斑坐标(X₀，Y₀)代入上式就可求得星点质心运动误差的补偿量。

为了说明本发明的补偿方法，以一具体仿真实例进行说明。卫星是星敏感器的典型载体，其运行情况一般分为两种：(1)在空间中保持高精度定向(姿态稳定)；(2)在空间中按预定要求进行机动(姿态机动)。这种姿态机动就是控制卫星在大范围内进行快速姿态角调整。而在这个姿态调控过程中，姿态角速度会较大，且变化速度较快，会使得星敏感器在曝光时间内的成像星点在像平面内发生明显移动现象。以敏捷卫星的快速机动为例，仿真运动误差补偿在高动态下的实现效果，以下仿真分为两种机动情况：(1)小角度机动(小于3°/s)；(2)大角度机动(大于等于3°/s)。

星敏感器仿真参数：焦距f：46.16mm，像素阵列：2048×2048，视场：14.5°×14.5°，像素大小：5.5μm×5.5μm，主点：1024×1024。角速度的协方差计算公式计算得到三轴角速度的方差大致分别为(1.0*10^-6，1.0*10^-6，1.8*10^-6)。由于预测过程中角速度拟合阶次较高(3阶)，Q阵可以设定为较小，拟定为1.0*10^-6I_4*4。

(1)小角度机动

运动参数：初始姿态四元数：(0，0，0，1)，初始角速度设定：(0，0.015，0)rad/s，目标姿态四元数：(0.242，0.242，0，0.94)，目标角速度：(0，0，0.015)rad/s，机动角度：50°，机动时间为：22s，曝光时间：80ms，更新率：100ms。图3为小角度机动下姿态四元数和角速度的变化曲线。图4为小角度下未补偿与补偿后误差比较图。

(2)大角度机动

运动参数：初始姿态四元数：(0.42，0.42，0.63，0.5)，初始角速度和目标角速度均为：(0，0，0)，目标姿态四元数：(0，0，0，1)，机动角度：120°，机动时间：30s，曝光时间：80ms，更新率：100ms。图6为大角度机动下姿态四元数和角速度的变化曲线。图7为大角度机动下未补偿与补偿后误差比较图。

从图5和图7实线可以看出，两种机动情况下未补偿时运动误差都较大，且随着角速度的变化会有明显变化，这是因为此误差是由三轴角速度和星点位置决定，且主要影响因素为三轴角速度的变化。当三轴角速度变化速度较快时，如仿真中机动开始或机动结束，运动误差会明显变大，由此可见运动误差的大小主要与其角速度的变化速度相关，与其角速度本身大小关系不大。经计算第一种情况下未补偿下运动误差的平均值为0.0720pixel，第二种情况下其平均值为0.1572pixel。反观两图虚线，经过补偿后，误差的变化幅度变小，且误差数值也变小。另外可以看出，当角速度变化较快时，虽然运动误差较大，但补偿效果比角速度变化慢时要好。这主要是因为角速度变化快时，其变化规律近似线性，角速度的一阶多项式表示精度较高，从而补偿量计算较准确。而当角速度变化曲线曲率较大时，角速度一阶多项式表示精度较差，从而也使得补偿效果较差。总体来看，两种情况下补偿后误差平均值为0.0169pixel、0.0547pixel。两次的补偿效果分别减少了76.5％、65.2％的运动误差，补偿效果较好。

对于星图预处理，曝光时间起始时刻的星点光斑坐标比较难获取，为了将系统简化，直接用质心坐标代替。图8为两种机动情况下用起始光斑坐标和质心坐标补偿效果对比图。如图8所示，用起始光斑坐标和用质心两者的补偿效果基本差别不大。由此可见，将星点质心代替曝光时间起始光斑坐标对补偿量的影响很小，这是因为星点质心坐标与星点起始时刻光斑坐标本身相差不大，两者的方向矢量相差很小。另外，补偿量主要由三轴角速度决定，起始光斑方向矢量只是次要影响因素。最终导致将星点质心代替曝光起始时刻光斑坐标对补偿量的影响很小，可以省去提取星点的起始光斑坐标的麻烦，直接用星点质心坐标代替。

以上所述，仅为本发明的较佳仿真实例而已，并非用于限定本发明的保护规范。

Claims (5)

1.一种高动态下星点质心误差补偿方法，其特征在于，所述方法包括以下步骤：

步骤201，三轴角速度的计算，即对三轴角速度分别进行最小二乘法拟合，并对拟合结果进行卡尔曼滤波；

步骤202，将曝光时间内三轴角速度近似用一阶多项式表示，将一阶多项式积分得到三轴旋转的角度信息，代入星敏感器的旋转矩阵；

步骤203，对星敏感器的旋转矩阵进行二阶泰勒展开近似，经透视投影变换，得到星点光斑在图像传感器上的移动轨迹坐标表达式；

步骤204，对所述移动轨迹坐标表达式进行二阶泰勒展开近似，参考质心法公式求解星点质心坐标，计算出星点质心误差补偿量；所述质心法公式求解星点质心坐标的公式为：

$\begin{array}{c}\stackrel{\&OverBar;}{x}=\frac{{\&Integral;}_0^{T_e}X\left(t\right)dt}{T_e}\&ap;\&lsqb;(k_1-\frac{k_2{k}_6+k_5{X}_0}{f}+\frac{k_6^2{X}_0}{f^2})\frac{T_e^2}{3}+(k_2-\frac{k_6{X}_0}{f})\frac{T_e}{2}+X_0\&rsqb;\\ \stackrel{\&OverBar;}{y}=\frac{{\&Integral;}_0^{T_e}Y\left(t\right)dt}{T_e}\&ap;\&lsqb;(k_3-\frac{k_{4}g+k_5{Y}_0}{f}+\frac{k_6^2{Y}_0}{f^2})\frac{T_e^2}{3}+(k_4-\frac{k_6{Y}_0}{f})\frac{T_e}{2}+Y_0\&rsqb;\end{array}---\left(11\right)$

其中，为高动态星点的质心，T_e为曝光时间，(k₁，k₂)为x(t)的二阶系数和一阶系数，(k₃，k₄)为y(t)的二阶系数和一阶系数，(k₅，k₆)为z(t)的二阶系数和一阶系数，(X₀，Y₀)为起始时刻星点光斑坐标，f为镜头焦距。

2.根据权利要求1所述的高动态下星点质心误差补偿方法，其特征在于，所述步骤201中的卡尔曼滤波具体为：设角速度的拟合时间多项式的系数为系统状态向量，以角速度的最小二乘拟合公式建立量测模型，以多项式时间零点的转移建立系统预测模型，对三轴角速度拟合分别进行滤波；

所述角速度的最小二乘拟合公式建立量测模型为：设观测向量为Z(n)＝[ω_x(n) ω_x(n-1) ω_x(n-2) ω_x(n-3)]，由多项式拟合公式得量测模型：

$F=\left[\begin{array}{ccc}1& -\mathrm{\&Delta;}t+T_e/2& {(-\mathrm{\&Delta;}t+T_e/2)}^2{(-\mathrm{\&Delta;}t+T_e/2)}^3\\ 1& -2\mathrm{\&Delta;}t+T_e/2& {(-2\mathrm{\&Delta;}t+T_e/2)}^2{(-2\mathrm{\&Delta;}t+T_e/2)}^3\\ 1& -3\mathrm{\&Delta;}t+T_e/2& {(-3\mathrm{\&Delta;}t+T_e/2)}^2{(-3\mathrm{\&Delta;}t+T_e/2)}^3\\ 1& -4\mathrm{\&Delta;}t+T_e/2& {(-4\mathrm{\&Delta;}t+T_e/2)}^2{(-4\mathrm{\&Delta;}t+T_e/2)}^4\end{array}\right]--\left(3\right)$

Z(n)＝H(n)X(n)+V(n)＝(F^T)^-1(F^TF)X(n)+V(n)

其中，T_e为每帧的曝光时间，V(n)为观测噪声，系统过程噪声W(n)和测量噪声V(n)均为相互独立的零均值的高斯白噪声。

3.根据权利要求1所述的高动态下星点质心误差补偿方法，其特征在于，所述步骤202中代入星敏感器的旋转矩阵具体为：

三轴角速度用一阶多项式表示后，将其积分可得到曝光时间内三轴旋转的角度信息，将星点光斑坐标表示成星敏感器坐标系下的单位方向矢量，将所述角度信息代入旋转矩阵，从而建立星点光斑方向矢量随时间的运动方程。

4.根据权利要求1所述的高动态下星点质心误差补偿方法，其特征在于，所述步骤203中对旋转矩阵的二阶泰勒展开的具体为：

由于旋转矩阵的复杂和非线性，对旋转矩阵每个元素进行二阶泰勒展开近似，将旋转矩阵每个元素化为角速度参数的二阶时间多项式，这样就将星点光斑的方向矢量用时间二阶多项式表示。

5.根据权利要求1所述的高动态下星点质心误差补偿方法，其特征在于，所述步骤204中对移动轨迹坐标表达式进行二阶泰勒展开近似的具体为：

经透视投影变换后，移动轨迹坐标表达式分子分母均为二阶时间多项式，对其进行二阶泰勒展开近似，化为简单的二阶时间多项式，系数均由角速度参数表示，计算出星点质心和星点理想中心，将星点质心和星点理想中心相减得出星点质心误差补偿量。