CN1923621A - 基于径向排列约束的星敏感器在轨校准方法 - Google Patents

Abstract

本发明属于航天测量技术，涉及对星敏感器在轨校准方法的改进。校准的步骤如下：建立星敏感器姿态转换矩阵；星敏感器径向排列约束成像；进行图像采集和星图识别；处理单帧星图；进行多帧星图的整体优化；本发明利用径向排列约束，闭式求解外部参数，使得数据估计有良好稳定性；需迭代优化的参数少，计算速度快；径向畸变代表了主要畸变，计算精度较高。

Description

基于径向排列约束的星敏感器在轨校准方法

技术领域

本发明属于航天测量技术，涉及对星敏感器在轨校准方法的改进。

背景技术

星敏感器是一种利用恒星观测，为空间飞行器提供高精度姿态信息的航天测量仪器。星敏感器在轨校准方法是星敏感器研究和应用中的一项关键技术。通常星敏感器发射前，其内部参数如主点，焦距和畸变系数等会在地面进行标定。标定方法有夜空标定和星光实验室标定等。但是在飞行器发射后，由于发射时的冲击和工作环境的变化情况，如重力、大气和温度等都会不同于地面情况，使得星敏感器的内部参数出现偏差，另外长期使用后老化也会影响工作精度。

现有的在轨校准方法主要分两类，一类是依赖于姿态信息的校准；一类是根据星内角距不变原理(不依赖于姿态信息)的校准。依赖于姿态信息的校准是由外部提供一个已知的姿态，根据此姿态可以计算出当前视场恒星在星敏感器靶面上的理想位置，如星敏感器内部参数，如主点，焦距和畸变系数发生变化时，恒星测量位置将会和理想位置发生偏差。根据获得的一系列偏差值，可以重新校准内部参数。该方法一个明显的缺陷是由于需要外部提供一个姿态，从而外部姿态如果有误差，则该误差会引入校准过程中，影响校准精度。

不依赖于姿态信息的在轨校准方法是根据星内角距正交变换不变原理，检测在轨飞行过程中星内角距测量值和真实值之间的偏差，利用基函数来拟合角距偏差。该方法存在的缺点是，由于拟合过程采用的是星内角相对值，误差的0阶项被消去；同时，基函数选择对拟合精度影响较大，基函数选择不合适，结果可能出现数值不稳定的情况。

另外美国JPL实验室有采用径向基神经网络来进行星敏感器在轨校准。该神经网络基本原理也是一种依赖姿态的方法，只是利用神经网络的并行运算能力来处理估计星向量和测量星向量之间的偏差。该方法缺点是精度不高，运算资源消耗比较大。

发明内容

本发明的目的是：针对目前星敏感器在轨校准存在的问题，提出一种基于径向排列约束的、不依赖于外部姿态信息的在轨校准方法。该方法利用镜头畸变通常以径向畸变为主的情况，根据径向排列约束来解决于星敏感器外部姿态和内部参数的耦合问题。该校准方法采用闭式推导，计算过程稳定性好，能够快速完成，并有良好的校准精度。

本发明的技术方案是：基于径向排列约束的星敏感器在轨校准方法，其特征在于，校准的步骤如下：

1、建立星敏感器姿态转换矩阵；

1.1、建立天球坐标系和星敏感器坐标系；设O’-XnYnZn为天球坐标系，O-XYZ为星敏感器坐标系，星敏感器的姿态角由赤经α₀、赤纬β₀、和滚转角φ₀组成，这里α₀为Z轴在XnYn面上的投影与Xn轴的夹角；β₀为Z轴与它在XnYn面上的投影之间的夹角；φ₀为子午面和XY面交线与Y轴的夹角；

1.2、建立星敏感器姿态转换矩阵；从天球坐标系到星敏感器坐标系的转换过程为：第一步，绕Zn轴转动 角，使得Xn轴和子午面垂直；第二步，绕新的Xn轴转动 角，使得Zn轴和Z轴一致；第三步，绕新的Zn轴转动φ₀角，使得天球坐标系和星敏感器坐标系重合。则该转动过程的方向余弦矩阵R即为姿态转换矩阵，表示为：

式中，s表示正弦函数，c表示余弦函数；简化表示为： $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$

2、星敏感器径向排列约束成像；

2.1、设星光方向矢量为： $w=\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[2\right]$

其中w₁，w₂，w₃分别为星光矢量在天球坐标系Xn，Yn，Zn轴的投影，α，β为恒星的赤经、赤纬坐标；

2.2、星光在星敏感器的图像传感器靶面上投影成像，O-XYZ为星敏感器坐标系，o-xy为2维靶面坐标系，∏为靶面，Oo之间距离为焦距，星光向量w理想透视成像点为P’，由于发生径向畸变，实际成像点为P；

2.3、建立星敏感器成像模型如下：

$x=f_c(1+k_1{r}^2)\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}---\left[3\right]$

$y=f_c(1+k_1{r}^2)\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

式中，f_c为焦距，k1为径向畸变系数，x，y为靶面坐标系坐标，当采集硬件获得图像后有：

                X＝S_x·x/DX+X₀    [4]

                Y＝y/DY+Y₀

式中，Sx为采样比例因子，DX，DY表示像素中心间距离，X0，Y0为图像中心，摄像机的主点坐标采用地面校准得到的主点坐标；

3、进行图像采集和星图识别；星敏感器在轨采集图像时，间隔1秒采集一帧图像，连续采集10～20帧图像；利用鲁棒性强的星图识别方法，对星图中的恒星进行识别，获得每颗恒星的赤经、赤纬坐标；

4、处理单帧星图；

4.1、计算姿态矩阵R和比例因子Sx；

4.1.1、估计姿态矩阵R；将图像坐标转换到以图像中心为原点的坐标系中，设：

        x_d＝(X-X₀)·DX

        y_d＝(Y-Y₀)·DY             [5]

这里，X，Y为以像素为单位的恒星图像坐标，X₀，Y₀为以像素为单位的星敏感器主点坐标，Dx，Dy为图像传感器像元在x，y方向的尺寸，单位毫米；

根据径向排比不变原理，有：

$\frac{x_d}{s_x{y}_d}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_4{w}_1+r_5{w}_2+r_6{w}_3}---\left[6\right]$

s_xr₁w₁+s_xr₂w₂+s_xr₃w₃＝x_dr₄w₁+x_dr₅w₂+x_dr₆w₃    [7]

利用奇异值分解方法求解齐次方程，

由上式得到：

y_dw₁s_xr₁+y_dw₂s_xr₂+y_dw₃s_xr₃-x_dw₁r₄-x_dw₂r₅-x_dw₃r₆＝0  [8]

设 $\stackrel{\→}{h}=\left(\begin{array}{cccccc}{h}_1& h_2& h_3& h_4& h_5& h_6\end{array}\right),$ 其中h₁＝s_xr₁，h₂＝s_xr₂，h₃＝s_xr₃，h₄＝r₄，h₅＝r₅，h₆＝r₆，有：

$\left[\begin{array}{cccccc}{y}_d{w}_1& y_d{w}_2& y_d{w}_3& -x_d{w}_1& -x_d{w}_2& -x_d{w}_3\end{array}\right]\stackrel{\→}{h}=0---\left[9\right]$

根据每一个采样点都可以得到上式，于是它们可以组成齐次方程：

$A\stackrel{\→}{h}=0---\left[10\right]$

该齐次方程的解为A矩阵的最小奇异值对应的右奇异值向量，该解向量和

有一个比例关系，记为

4.1.2、计算比例因子Sx；

姿态矩阵R为正交矩阵，因此有约束：

$r_1^2+r_2^2+r_3^2=1---\left[11\right]$

$r_4^2+r_5^2+r_6^2=1$

于是，设 $c=\sqrt{{\tilde{h}}_4^2+{\tilde{h}}_5^2+{\tilde{h}}_6^2}---\left[12\right]$

为齐次方程解和真实解的比例系数，即 $\stackrel{\→}{h}=\tilde{h}/c,$ 得到： $s_x=\frac{1}{c}\sqrt{({\left({\tilde{h}}_1\right)}^2+{\left({\tilde{h}}_2\right)}^2+{\left({\tilde{h}}_3\right)}^2)}---\left[13\right]$

4.1.3、计算r_1～9；从齐次解可以得到：

$r_4=\frac{{\tilde{h}}_4}{c},r_5=\frac{{\tilde{h}}_5}{c},r_6=\frac{{\tilde{h}}_6}{c},r_1=\frac{{\tilde{h}}_1}{s_{x}c},r_2=\frac{{\tilde{h}}_2}{s_{x}c},r_3=\frac{{\tilde{h}}_3}{s_{x}c}---\left[14\right]$

r₇，r₈和r₉可以根据R阵前两行的叉积获得，即：

r₇＝r₂r₆-r₃r₅，r₈＝r₃r₄-r₁r₆，r₉＝r₁r₅-r₂r₄        [15]

根据[12]可知，此处默认比例系数c为正，实际上c有可能为负，因此需要进一步判断c的符号，以当前视场内距离中心最远的一个星作为判断星点，设此星点序号为i，其天球坐标系内的向量为：

           ν_i＝[n_i1 n_i2  n_i3]^T

记： ${\stackrel{\‾}{x}}_i=\frac{r_1{n}_{i1}+r_2{n}_{i2}+r_3{n}_{i3}}{r_7{n}_{i1}+r_8{n}_{i2}+r_9{n}_{i3}}---\left[16\right]$

${\stackrel{\‾}{y}}_i=\frac{r_4{n}_{i1}+r_5{n}_{i2}+r_6{n}_{i3}}{r_7{n}_{i1}+r_8{n}_{i2}+r_9{n}_{i3}}$

设x_id，y_id为以主点为原点的图像测量坐标，如果sign( x_i)＝sign(x_id)且sign( y_i)＝sign(y_id)，则表示估计投影位置和实际投影位置在同一象限，c符号为正；否则c符号为负，姿态矩阵R的参数r_1～6取反；

4.1.4、R矩阵正交化；

由于参数r_1～6是通过解齐次方程求出，因此需要对R矩阵进行正交化处理，设 为前面计算得到姿态矩阵，正交化的目标是寻找正交矩阵R使得：

$\underset{R}{\min }{\left|\right|R-\tilde{R}\left|\right|}^2---\left[17\right]$

设

阵的奇异值分解为USV^T，于是可以求出正交矩阵R＝UV^T；

4.2、计算焦距f_c和径向畸变系数k₁；

根据公式[3]，利用最速下降法来优化计算焦距f_c和径向畸变系数k1，焦距f_c的初始值可以通过地面校准获得，并假设畸变k1初始值设为0，设：

$F_x(f_c,k_1)=f_c(1+k_1{r}^2)\frac{r_1{w}_1+r_2{w}_3+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}---\left[18\right]$

$F_y(f_c,k_1)=f_c(1+k_1{r}^2)\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

求其偏导数：

$\mathrm{dx}=\tilde{x}-x=\mathrm{Ax}\left[\begin{array}{c}{\mathrm{\Δf}}_c\\ {\mathrm{\Δk}}_1\end{array}\right]---\left[19\right]$

$\mathrm{dy}=\tilde{y}-y=\mathrm{Ay}\left[\begin{array}{c}\mathrm{\Δ}{f}_c\\ \mathrm{\Δ}{k}_1\end{array}\right]$

这里， $\mathrm{Ax}=\left(\begin{array}{cc}\frac{{\∂F}_x}{\∂f_c}& \frac{{\∂F}_x}{{\∂k}_1}\end{array}\right),\mathrm{Ay}=\left(\begin{array}{cc}\frac{{\∂F}_y}{{\∂f}_c}& \frac{{\∂F}_y}{{\∂k}_1}\end{array}\right)$ 称为敏感向量，所有该星图星点组成敏感矩阵，将估计点和测量点带入，即可求出焦距估计偏差和径向畸变系数估计偏差，通过迭代计算1～3次即可获得稳定的数值解；

5、进行多帧星图的整体优化；多帧星图的总星点数积累到80～120后，进行参数整体优化，获得稳定的最优解。优化方法可以采用Marquardt最小二乘法。

本发明的优点是：

(1)利用径向排列约束，闭式求解外部参数，使得数据估计有良好稳定性；

(2)需迭代优化的参数少，计算速度快；

(3)径向畸变代表了主要畸变，计算精度较高。

附图说明

图1是本发明中的星敏感器姿态角示意图。

图2是发明中的星敏感器星光成像示意图。

图3是发明中基于径向排列约束的校准方法主要步骤示意图。

具体实施方式

下面对本发明做进一步详细说明。本发明的在轨校准方法，其特征在于，校准的步骤如下：

1、建立星敏感器姿态转换矩阵；

1.1、建立天球坐标系和星敏感器坐标系；如图1所示，设O’-XnYnZn为天球坐标系，O-XYZ为星敏感器坐标系，星敏感器的姿态角由赤经α₀、赤纬β₀、和滚转角φ₀组成，这里α₀为Z轴在XnYn面上的投影与Xn轴的夹角；β₀为Z轴与它在XnYn面上的投影之间的夹角；φ₀为子午面和XY面交线与Y轴的夹角；

1.2、建立星敏感器姿态转换矩阵；从天球坐标系到星敏感器坐标系的转换过程为：第一步，绕Zn轴转动

角，使得Xn轴和子午面垂直；第二步，绕新的Xn轴转动

角，使得Zn轴和Z轴一致；第三步，绕新的Zn轴转动φ₀角，使得天球坐标系和星敏感器坐标系重合。则该转动过程的方向余弦矩阵R即为姿态转换矩阵，表示为：

式中，s表示正弦函数，c表示余弦函数；

简化表示为： $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$

2、星敏感器径向排列约束成像；

2.1、设星光方向矢量为： $w=\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[2\right]$

这里，X，Y为以像素为单位的恒星图像坐标，X₀，Y₀为以像素为单位的星敏感器主点坐标，Dx，Dy为图像传感器像元在x，y方向的尺寸，单位毫米。

2.2、星光在星敏感器的图像传感器靶面上投影成像，如图2所示为星敏感器径向畸变成像示意图。O-XYZ为星敏感器坐标系，o-xy为2维靶面坐标系，∏为靶面，Oo之间距离为焦距，星光向量w理想透视成像点为P’，由于发生径向畸变，实际成像点为P；

2.3、建立星敏感器成像模型如下：

$x=f_c(1+k_1{r}^2)\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}---\left[3\right]$

$y=f_c(1+k_1{r}^2)\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

式中，f_c为焦距，k1为径向畸变系数，x，y为靶面坐标系坐标，当采集硬件获得图像后有：

       X＝S_x·x/DX+X₀             [4]

       Y＝y/DY+Y₀

式中，Sx为采样比例因子，DX，DY表示像素中心间距离，X0，Y0为图像中心，摄像机的主点坐标采用地面校准得到的主点坐标；

根据径向排比约束(RAC，Radial Alignment Constraint)的特点，外部姿态参数独立于径向畸变，姿态和径向畸变没有耦合关系，因此可以将外部姿态和畸变分开进行校准。如图3所示，本发明方法在经过了图像采集和星图识别后，分两个阶段完成校准，首先是单帧星图处理，然后是多帧星图的整体优化。单帧星图处理分两步，第一步为姿态矩阵R和比例因子Sx的计算；第二步为焦距和畸变系数计算。多帧星图的总星点数积累到一定数值后，进行参数整体优化，以获得稳定的最优解。下面详细阐述该方法。

3、进行图像采集和星图识别；星敏感器在轨采集图像时，间隔1秒采集一帧图像，连续采集10～20帧图像；利用鲁棒性强的星图识别方法，对星图中的恒星进行识别，获得每颗恒星的赤经、赤纬坐标。

4、处理单帧星图；

4.1、计算姿态矩阵R和比例因子Sx；

4.1.1、估计姿态矩阵R；将图像坐标转换到以图像中心为原点的坐标系中，设：

            x_d＝(X-X₀)·DX             [5]

            y_d＝(Y-Y₀)·DY

这里，X，Y为以像素为单位的恒星图像坐标，X₀，Y₀为以像素为单位的星敏感器主点坐标，Dx，Dy为图像传感器像元在x，y方向的尺寸，单位毫米。

根据径向排比不变原理，有：

$\frac{x_d}{s_x{y}_d}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_4{w}_1+r_5{w}_2+r_6{w}_3}---\left[6\right]$

s_xr₁w₁+s_xr₂w₂+s_xr₃w₃＝x_dr₄w₁+x_dr₅w₂+x_dr₆w₃          [7]

利用奇异值分解方法求解齐次方程，

由上式得到：

y_dw₁s_xr₁+y_dw₂s_xr₂+y_dw₃s_xr₃-x_dw₁r₄-x_dw₂r₅-x_dw₃r₆＝0      [8]

设 $\stackrel{\→}{h}=\left(\begin{array}{cccccc}{h}_1& h_2& h_3& h_4& h_5& h_6\end{array}\right),$ 其中h₁＝s_xr₁，h₂＝s_xr₂，h₃＝s_xr₃，h₄＝r₄，h₅＝r₅，h₆＝r₆，有：

$\left[\begin{array}{cccccc}{y}_d{w}_1& y_d{w}_2& y_d{w}_3& -x_d{w}_1& -x_d{w}_2& -x_d{w}_3\end{array}\right]\stackrel{\→}{h}=0---\left[9\right]$

根据每一个采样点都可以得到上式，于是它们可以组成齐次方程：

$A\stackrel{\→}{h}=0---\left[10\right]$

该齐次方程的解为A矩阵的最小奇异值对应的右奇异值向量，该解向量和

有一个比例关系，记为

4.1.2、计算比例因子Sx；

姿态矩阵R为正交矩阵，因此有约束：

$r_1^2+r_2^2+r_3^2=1---\left[11\right]$

$r_4^2+r_5^2+r_6^2=1$

于是，设 $c=\sqrt{{\tilde{h}}_4^2+{\tilde{h}}_5^2+{\tilde{h}}_6^2}---\left[12\right]$

为齐次方程解和真实解的比例系数，即 $\stackrel{\→}{h}=\tilde{h}/c,$

得到： $s_x=\frac{1}{c}\sqrt{({\left({\tilde{h}}_1\right)}^2+{\left({\tilde{h}}_2\right)}^2+{\left({\tilde{h}}_3\right)}^2)}---\left[13\right]$

4.1.3、计算r_1～9；从齐次解可以得到：

$r_4=\frac{{\tilde{h}}_4}{c},r_5=\frac{{\tilde{h}}_5}{c},r_6=\frac{{\tilde{h}}_6}{c},r_1=\frac{{\tilde{h}}_1}{s_{x}c},r_2=\frac{{\tilde{h}}_2}{s_{x}c},r_3=\frac{{\tilde{h}}_3}{s_{x}c}---\left[14\right]$

r₇，r₈和r₉可以根据R阵前两行的叉积获得，即：

r₇＝r₂r₆-r₃r₅，r₈＝r₃r₄-r₁r₆，r₉＝r₁r₅-r₂r₄       [15]

根据[12]可知，此处默认比例系数c为正，实际上c有可能为负，因此需要进一步判断c的符号，以当前视场内距离中心最远的一个星作为判断星点，设此星点序号为i，其天球坐标系内的向量为：

         ν_i＝[n_i1  n_i2  n_i3]^T

记： ${\stackrel{\‾}{x}}_i=\frac{r_1{n}_{i1}+r_2{n}_{i2}+r_3{n}_{i3}}{r_7{n}_{i1}+r_8{n}_{i2}+r_9{n}_{i3}}---\left[16\right]$

${\stackrel{\‾}{y}}_i=\frac{r_4{n}_{i1}+r_5{n}_{i2}+r_6{n}_{i3}}{r_7{n}_{i1}+r_8{n}_{i2}+r_9{n}_{i3}}$

设x_id，y_id为以主点为原点的图像测量坐标，如果sign( x_i)＝sign(x_id)且sign( y_i)＝sign(y_id)，则表示估计投影位置和实际投影位置在同一象限，c符号为正；否则c符号为负，姿态矩阵R的参数r_1-6取反；

4.1.4、R矩阵正交化；

由于参数r_1～6是通过解齐次方程求出，因此需要对R矩阵进行正交化处理，设

为前面计算得到姿态矩阵，正交化的目标是寻找正交矩阵R使得：

$\underset{R}{\min }{\left|\right|R-\tilde{R}\left|\right|}^2---\left[17\right]$

设

阵的奇异值分解为USV^T，于是可以求出正交矩阵R＝UV^T；

4.2、计算焦距f_c和径向畸变系数k₁；

根据公式[3]，利用最速下降法来优化计算焦距f_c和径向畸变系数k1，焦距f_c的初始值可以通过地面校准获得，并假设畸变k1初始值设为0，设：

$F_x(f_c,k_1)=f_c(1+k_1{r}^2)\frac{r_1{w}_1+r_2{w}_3+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}---\left[18\right]$

$F_y(f_c,k_1)=f_c(1+k_1{r}^2)\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

求其偏导数：

$\mathrm{dx}=\tilde{x}-x=\mathrm{Ax}\left[\begin{array}{c}{\mathrm{\Δf}}_c\\ {\mathrm{\Δk}}_1\end{array}\right]---\left[19\right]$

$\mathrm{dy}=\tilde{y}-y=\mathrm{Ay}\left[\begin{array}{c}\mathrm{\Δ}{f}_c\\ \mathrm{\Δ}{k}_1\end{array}\right]$

这里， $\mathrm{Ax}=\left(\begin{array}{cc}\frac{{\∂F}_x}{\∂f_c}& \frac{{\∂F}_x}{{\∂k}_1}\end{array}\right),\mathrm{Ay}=\left(\begin{array}{cc}\frac{{\∂F}_y}{{\∂f}_c}& \frac{{\∂F}_y}{{\∂k}_1}\end{array}\right)$ 称为敏感向量，所有该星图星点组成敏感矩阵，将估计点和测量点带入，即可求出焦距估计偏差和径向畸变系数估计偏差，通过迭代计算1～3次即可获得稳定的数值解；

5、进行多帧星图的整体优化；多帧星图的总星点数积累到80～120后，进行参数整体优化，获得稳定的最优解。优化方法可以采用Marquardt最小二乘法。

仿真和结果分析

仿真的星敏感器基本参数为：

视场：12×12；

像素阵列：1024×1024；

像素尺寸：0.015mm×0.015mm；

比例因子：1.05；

焦距：73.0703mm；

径向畸变系数：-0.0005；

仿真采用了10帧星图，姿态分别是：

	1	2	3	4	5	6	7	8	9	10
											赤经	-45	-35	-25	-15	-5	5	15	25	35	45
赤纬	-35	-25	-15	-5	5	15	25	35	45	55
											滚转	20	20	20	20	20	20	20	20	20	20
星数	17	14	9	15	17	11	19	13	17	17

这里假设中心像素没有偏差，为(512，512)。

在没有噪声的情况下，仿真结果表明，计算得到的姿态以及内部参数和设定的内部参数没有偏差。说明该方法是正确的。

在每个星点上添加位置噪声，噪声设定为0均值，0.05像素误差的高斯噪声。仿真结果如下：

1	2	3	4	5	6	7	8	9	10
										-45.00	-35.001	-25.000	-15.0	-5.000	5.000	15.000	25.00	35.00	44.999
-35.00	-24.999	-15.001	-5.00	5.000	15.00	25.000	34.999	45.00	55.00
										19.999	19.996	19.997	19.998	20.001	20.004	20.004	20.001	20.002	20.004

赤经偏差均值为1.1角秒，赤纬偏差均值为1.0角秒，滚转偏差均值为9.7角秒。该结果符合该噪声水平下的姿态精度，说明计算精度可以接受。

比例因子采用10帧图像计算的平均值，Sx＝1.04995。

仿真假设初始焦距为72mm，径向畸变系数：k1＝0。最速下降法迭代两次后，数据稳定，如下表所示。

0

1

2

ff	72	73.072	73.072
				k1	0	-5.2225*10-4	-4.998*10-4

最后将校准后的畸变来修正最初的数据，进行参数的整体优化。整体优化后得到：

Sx＝1.0500，ff＝73.0723，k₁＝-5.0031×10^-4。代入一组验证数据得到残留像素误差均方根为x方向0.063像素，y方向0.053像素，对比0.05像素的噪声水平说明校准结果可以满意。

Claims (1)

1、基于径向排列约束的星敏感器在轨校准方法，其特征在于，校准的步骤如下：

1.1、建立星敏感器姿态转换矩阵；

1.1.1、建立天球坐标系和星敏感器坐标系；设O’-XnYnZn为天球坐标系，O-XYZ为星敏感器坐标系，星敏感器的姿态角由赤经α₀、赤纬β₀、和滚转角φ₀组成，这里α₀为Z轴在XnYn面上的投影与Xn轴的夹角；β₀为Z轴与它在XnYn面上的投影之间的夹角；φ₀为子午面和XY面交线与Y轴的夹角；

1.1.2、建立星敏感器姿态转换矩阵；从天球坐标系到星敏感器坐标系的转换过程为：第一步，绕Zn轴转动 角，使得Xn轴和子午面垂直；第二步，绕新的Xn轴转动 角，使得Zn轴和Z轴一致；第三步，绕新的Zn轴转动φ₀角，使得天球坐标系和星敏感器坐标系重合。则该转动过程的方向余弦矩阵R即为姿态转换矩阵，表示为：

式中，s表示正弦函数，c表示余弦函数；

简化表示为： $R=\left[\begin{array}{ccc}{r}_1& r_2& r_3\\ r_4& r_5& r_6\\ r_7& r_8& r_9\end{array}\right]$

1.2、星敏感器径向排列约束成像；

1.2.1、设星光方向矢量为： $\left(\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[2\right]$

其中w₁，w₂，w₃分别为星光矢量在天球坐标系Xn，Yn，Zn轴的投影，α，β为恒星的赤经、赤纬坐标；

1.2.2、星光在星敏感器的图像传感器靶面上投影成像，O-XYZ为星敏感器坐标系，o-xy为2维靶面坐标系，∏为靶面，Oo之间距离为焦距，星光向量w理想透视成像点为P’，由于发生径向畸变，实际成像点为P；

1.2.3、建立星敏感器成像模型如下：

$x=f_c(1+k_1{r}^2)\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

$y=f_c(1+k_1{r}^2)\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$ [3]

式中，f_c为焦距，k1为径向畸变系数，x，y为靶面坐标系坐标，当采集硬件获得图像后有：

                      X＝S_x·x/DX+X₀

                      Y＝y/DY+Y₀               [4]

式中，Sx为采样比例因子，DX，DY表示像素中心间距离，X0，Y0为图像中心，摄像机的主点坐标采用地面校准得到的主点坐标；

1.3、进行图像采集和星图识别；星敏感器在轨采集图像时，间隔1秒采集一帧图像，连续采集10～20帧图像；利用鲁棒性强的星图识别方法，对星图中的恒星进行识别，获得每颗恒星的赤经、赤纬坐标；

1.4、处理单帧星图；

1.4.1、计算姿态矩阵R和比例因子Sx；

1.4.1.1、估计姿态矩阵R；将图像坐标转换到以图像中心为原点的坐标系中，设：

                     x_d＝(X-X₀)·Dx

                     y_d＝(Y-Y₀)·Dy                [5]

这里，X，Y为以像素为单位的恒星图像坐标，X₀，Y₀为以像素为单位的星敏感器主点坐标，Dx，Dy为图像传感器像元在x，y方向的尺寸，单位毫米；

根据径向排比不变原理，有：

$\frac{x_d}{s_x{y}_d}=\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_4{w}_1+r_5{w}_2+r_6{w}_3}---\left[6\right]$

             s_xr₁w₁+s_xr₂w₂+s_xr₃w₃＝x_dr₄w₁+x_dr₅w₂+x_dr₆w₃    [7]

利用奇异值分解方法求解齐次方程，

由上式得到：

         y_dw₁s_xr₁+y_dw₂s_xr₂+y_dw₃s_xr₃-x_dw₁r₄-x_dw₂r₅-x_dw₃r₆＝0    [8]

设 $\stackrel{\→}{h}=\left(\begin{array}{cccccc}{h}_1& h_2& h_3& h_4& h_5& h_6\end{array}\right),$ 其中h₁＝s_xr₁，h₂＝s_xr₂，h₃＝s_xr₃，h₄＝r₄，h₅＝r₅，h₆＝r₆，有：

$\left[\begin{array}{cccccc}{y}_d{w}_1& y_d{w}_2& y_d{w}_3& -x_d{w}_1& -x_d{w}_2& -x_d{w}_3\end{array}\right]\stackrel{\→}{h}=0---\left[9\right]$

根据每一个采样点都可以得到上式，于是它们可以组成齐次方程：

$A\stackrel{\→}{h}=0---\left[10\right]$

该齐次方程的解为A矩阵的最小奇异值对应的右奇异值向量，该解向量和 有一个比例关系，记为

1.4.1.2、计算比例因子Sx；

姿态矩阵R为正交矩阵，因此有约束：

$r_1^2+r_2^2+r_3^2=1$

$r_4^2+r_5^2+r_6^2=1$ [11]

于是，设 $c=\sqrt{{\tilde{h}}_4^2+{\tilde{h}}_5^2+{\tilde{h}}_6^2}---\left[12\right]$

为齐次方程解和真实解的比例系数，即 $\stackrel{\→}{h}=\tilde{h}/c,$

得到： $s_x=\frac{1}{c}\sqrt{({\left({\tilde{h}}_1\right)}^2+{\left({\tilde{h}}_2\right)}^2+{\left({\tilde{h}}_3\right)}^2)}---\left[13\right]$

1.4.1.3、计算r_1～9；从齐次解可以得到：

$r_4=\frac{{\tilde{h}}_4}{c},r_5=\frac{{\tilde{h}}_5}{c},r_6=\frac{{\tilde{h}}_6}{c},r_1=\frac{{\tilde{h}}_1}{s_{x}c},r_2=\frac{{\tilde{h}}_2}{s_{x}c},r_3=\frac{{\tilde{h}}_3}{s_{x}c}---\left[14\right]$

r₇，r₈和r₉可以根据R阵前两行的叉积获得，即：

            r₇＝r₂r₆-r₃r₅，r₈＝r₃r₄-r₁r₆，r₉＝r₁r₅-r₂r₄    [15]

根据[12]可知，此处默认比例系数c为正，实际上c有可能为负，因此需要进一步判断c的符号。以当前视场内距离中心最远的一个星作为判断星点，设此星点序号为i，其天球坐标系内的向量为：

                     v_i＝[n_i1 n_i2 n_i3]^T

记：

${\stackrel{\‾}{x}}_i=\frac{r_1{n}_{i1}+r_2{n}_{i2}+r_3{n}_{i3}}{r_7{n}_{i1}+r_8{n}_{i2}+r_9{n}_{i3}}$

${\stackrel{\‾}{y}}_i=\frac{r_4{n}_{i1}+r_5{n}_{i2}+r_6{n}_{i3}}{r_7{n}_{i1}+r_8{n}_{i2}+r_9{n}_{i3}}$ [16]

设x_id，y_id为以主点为原点的图像测量坐标，如果sign( x_i)＝sign(x_id)且sign( y_i)＝sign(y_id)，则表示估计投影位置和实际投影位置在同一象限，c符号为正；否则c符号为负，姿态矩阵R的参数r_1～6取反；

1.4.1.4、R矩阵正交化；

由于参数r_1～6是通过解齐次方程求出，因此需要对R矩阵进行正交化处理，设 为前面计算得到姿态矩阵，正交化的目标是寻找正交矩阵R使得：

$m\underset{R}{i}n\left|\right|R-\tilde{R}|{|}^2---[17]$

设

阵的奇异值分解为USV^T，于是可以求出正交矩阵R＝UV^T；

1.4.2、计算焦距f_c和径向畸变系数k₁；

根据公式[3]，利用最速下降法来优化计算焦距f_c和径向畸变系数k1，焦距f_c的初始值可以通过地面校准获得，并假设畸变k₁初始值设为0，设：

$F_x(f_c,k_1)=f_c(1+k_1{r}^2)\frac{r_1{w}_1+r_2{w}_2+r_3{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$

$F_y(f_c,k_1)=f_c(1+k_1{r}^2)\frac{r_4{w}_1+r_5{w}_2+r_6{w}_3}{r_7{w}_1+r_8{w}_2+r_9{w}_3}$ [18]

求其偏导数：

$\mathrm{dx}=\tilde{x}-x=\mathrm{Ax}\left[\begin{array}{c}{\mathrm{\Δf}}_c\\ \mathrm{\Δ}{k}_1\end{array}\right]$

$\mathrm{dy}=\tilde{y}-y=\mathrm{Ay}\left[\begin{array}{c}{\mathrm{\Δf}}_c\\ \mathrm{\Δ}{k}_1\end{array}\right]$ [19]

这里， $\mathrm{Ax}=\left(\begin{array}{cc}\frac{{\∂F}_x}{\∂f_c}& \frac{{\∂F}_x}{{\∂k}_1}\end{array}\right),\mathrm{Ay}=\begin{array}{}\left(\begin{array}{cc}\frac{{\∂F}_y}{\∂f_c}& \frac{{\∂F}_y}{\∂k_1}\end{array}\right)\end{array}$ 称为敏感向量，所有该星图星点组成敏感矩阵，将估计点和测量点带入，即可求出焦距估计偏差和径向畸变系数估计偏差，通过迭代计算1～3次即可获得稳定的数值解；

1.5、进行多帧星图的整体优化；多帧星图的总星点数积累到80～120后，进行参数整体优化，获得稳定的最优解，优化方法可以采用Marquardt最小二乘法。

Images