CN100348947C - 基于weng模型的星敏感器在轨校准方法 - Google Patents

Abstract

本发明属于航天测量技术，涉及对星敏感器在轨校准方法的改进。校准的步骤如下：建立星敏感器姿态转换矩阵；建立星敏感器weng成像模型；数据采集；通过参数计算进行校准。本发明降低了外部参数和内部参数的耦合关联影响；畸变模型完整，有良好的适应性；以线性闭式求解参数为主，计算量适中。

Description

基于weng模型的星敏感器在轨校准方法

技术领域

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

背景技术

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

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

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

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

发明内容

本发明的目的是：针对目前星敏感器在轨校准存在的问题，提出一种基于weng模型的在轨校准方法。该方法利用机器视觉中的weng模型，分线性参数和非线性参数两步来估计内方元素。该方法畸变模型完整，可以处理复杂的星敏感器镜头畸变情况。数据处理分两步完成并且避免了复杂的参数优化，数据收敛快速稳定。

本发明的技术方案是：

基于weng模型的星敏感器在轨校准方法，其特征在于，校准的步骤如下：

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

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

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

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

R可简化表示为： $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],$ r₁～r₉表示R矩阵9个元素；

2、建立星敏感器weng成像模型；

设星光方向矢量为： $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₃表示星光向量在星敏感器坐标系下的x，y，z坐标值；α，β为恒星的赤经、赤纬坐标；

2.1、使用星敏感器对星光进行径向畸变成像；O-XYZ为星敏感器坐标系，o-xy为2维靶面坐标系，∏为靶面，Oo之间距离为焦距；星光向量w理想透视成像点为P′，由于发生径向畸变，实际成像点为P；

2.2、建立线性成像模型；由于x和y方向之间存在不确定比例因子，用f_u和f_v分别表示x和y方向的焦距，建立线性成像模型为：

$\frac{x^{\′}}{f_u}=\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}$ [3]

$\frac{y^{\′}}{f_v}=\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}$

式中，x′，y′为透视投影的星点位置，定义线性参数向量d₁＝(f_u f_v α₀ β₀ ₀)；

2.3、考虑到畸变影响，建立weng非线性畸变模型公式为：

$\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}=\hat{u}+{(g_1+g_3)\hat{u}}^2+g_4\hat{u}\hat{v}+g_1{\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$ [4]

$\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}=\hat{v}+g_2{\hat{u}}^2+g_3\hat{u}\hat{v}+(g_2+g_4){\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

式中， $\hat{u}=\frac{x}{f_u},$ $\hat{v}=\frac{y}{f_v},$ x，y为星点实际测量坐标，k₁为径向畸变系数，g₁～g₄为偏心畸变和薄棱镜畸变的综合公式系数，于是非线性的畸变参数有5个，定义参数向量d₂＝(k₁ g₁ g₂ g₃ g₄)；

3、数据采集；根据星敏感器采集得到的一帧星图，定义采集得到的全部星点数据集合为Ω₀，星点个数为n；以主点为中心，直径为图像尺寸1/2做圆，记圆内的数据集合为Ω₁，星点个数为n₁；

4、通过参数计算进行校准；

校准过程分为两步完成，根据中心数据Ω₁受畸变影响较小的原理，第一步利用Ω₁构造线性计算模型，计算得到参数向量d₁的初始值；第二步对d₁的参数进行优化处理；第三步，根据weng模型，构造畸变系数线性方程，在全局数据Ω₀下计算参数向量d₂；最后，判断数据结果是否满足要求，如果不满足则返回第二步进行迭代计算，如果满足则输出参数向量的最终解d₁和d₂；具体计算过程如下：

4.1、进行线性参数估计；根据公式[3]有：

x′w₁r₇+x′w₂r₈+x′w₃r₉＝w₁r₁f_u+w₂r₂f_u+w₃r₃f_u

                                               [5]

y′w₁r₇+y′w₂r₈+y′w₃r₉＝w₁r₄f_v+w₂r₅f_v+w₃r₆f_v

整理有：

$w_1\frac{r_1{f}_u}{r_9}+w_2\frac{r_2{f}_u}{r_9}+w_3\frac{r_3{f}_u}{r_9}-{\mathrm{xw}}_1\frac{r_7}{r_9}-{\mathrm{xw}}_2\frac{r_8}{r_9}={\mathrm{xw}}_3$ [6]

$w_1\frac{r_4{f}_v}{r_9}+w_2\frac{r_5{f}_v}{r_9}+w_3\frac{r_6{f}_v}{r_9}-{\mathrm{yw}}_1\frac{r_8}{r_9}-{\mathrm{yw}}_2\frac{r_8}{r_9}={\mathrm{yw}}_3$

假设 $r_1^{\′}=\frac{r_1{f}_u}{r_9},$ $r_2^{\′}=\frac{r_2{f}_u}{r_9},$ $r_3^{\′}=\frac{r_3{f}_u}{r_9},$ $r_4^{\′}=\frac{r_4{f}_v}{r_9},$ $r_5^{\′}=\frac{r_5{f}_u}{r_9},$ $r_6^{\′}=\frac{r_6{f}_v}{r_9},$ $r_7^{\′}=\frac{r_7}{r_9},$ $r_8^{\′}=\frac{r_8}{r_9},$ 共8个未知数，定义参数向量m＝[r₁′r₂′r₃′r₄′r₅′r₆′r₇′r₈′]^T，对于每个星点数据可以得到方程：

[w₁ w₂ w₃  0 0 0-xw₁  -xw₂]m＝xw₃

                                  [7]

[0 0 0 w₁ w₂ w₃-yw₁ -yw₂]m＝yw₃

对于数据集合Ω₁内的星点，可以得到2n₁个方程，构成“Am＝b”形式；A为2n₁个方程左式构成的矩阵，b为2n₁个方程右式构成的向量；可用最小二乘法使得 $\underset{m}{\min }\left|\right|A\hat{m}-b\left|\right|,$ 解得8个未知数r₁′～r₈′；然后根据R为正交矩阵特点，有：

$r_9=\sqrt{1+r_7^{\′}+r_8^{\′}}\text{,}$ r₇＝r₇′·r₉，  r₈＝r₈′·r₉

$f_u=r_9\·\sqrt{{r_1^{\′}}^2+{r_2^{\′}}^2+{r_3^{\′}}^2}$

r₁＝r₁′·r₉/r_u，  r₂＝r₂′·r₉/f_u，  r₃＝r₃′·r₉/f_u    [8]

$f_v=r_9\·\sqrt{{r_4^{\′}}^2+{r_5^{\′}}^2+{r_6^{\′}}^2}$

r₄＝r₄′·r₉/f_v，  r₅＝r₅′·r₉/f_v，  r₆＝r₆′·r₉/f_v

先假设r9为正，得到姿态估计矩阵，然后选择视场内的某一颗距离中心较远的恒星来判断r9是否真的为正，如果星光向量经姿态估计矩阵成像后坐标和测量坐标的符合一致，则说明r9为正，否则为负；

4.2、优化参数向量d₁；

设线性参数向量和非线性参数向量的真实值为d₁，d₂，估计值为

最优的估计值应该满足：

$\underset{d}{\min }\left|\right|P-f({\hat{d}}_1,{\hat{d}}_2)\left|\right|$

式中，P为点集Ω₁的星点坐标；以公式[8]得到的解为初值，采用非线性优化方法对d₁进行优化；设：

$x^{\′}=f_x(d_1,{\stackrel{\‾}{d}}_2)$

                                             [9]

$y^{\′}=f_y(d_1,{\stackrel{\‾}{d}}_2)$

式中，

表示参数向量d₂固定；由于f_x和f_y均为非线性函数，因此采用非线性线性化方法来估计参数向量d₁；假设

为估计值，

为测量值，Δd₁为向量估计偏差；

$\mathrm{\Δx}=\stackrel{\‾}{x}-\hat{x}\≈{\mathrm{A\Δd}}_1$

                             [10]

$\mathrm{\Δy}=\stackrel{\‾}{y}-\hat{y}\≈{\mathrm{B\Δd}}_1$

式中，A、B为敏感矩阵，求各项偏导数即可获得；

于是得到参数向量更新值d₁＝d₁+Δd₁；

4.3、计算畸变系数向量d₂；

得到线性参数后，下一步求解非线性参数d₂；根据公式[4]，构造线性估计模型：

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& {\hat{u}}^2+{\hat{v}}^2& 0& {\hat{u}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=\mathrm{du}$

                                         [11]

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& 0& {\hat{u}}^2+{\hat{v}}^2& {\hat{v}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=\mathrm{dv}$

$\mathrm{du}=\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}-\hat{u}$

式中：

$\mathrm{dv}=\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}-\hat{v}$

构成“Bd₂＝c”形式，B为2n个方程左式构成的矩阵，c为2n个方程右式构成的向量；利用最小二乘法使得 $\underset{d_2}{\min }\left|\right|B{d}_2-c\left|\right|,$ 求出畸变参数向量d₂的估计值，这里的数据为全部点集Ω₀；

得到模型的全部参数d₁和d₂后，通过对比理想坐标值和测量坐标值进行判断；如果误差不满足要求，则回到步骤1.4.2进一步优化d₁，进行迭代计算；如果误差满足要求，则停止迭代，输出结果。

本发明的优点是：

1、降低了外部参数和内部参数的耦合关联影响；

2、畸变模型完整，有良好的适应性；

3、以线性闭式求解参数为主，计算量适中。

附图说明

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

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

图3是本发明中的参数计算过程示意图。

图4是模拟星场示意图。

具体实施方式

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

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

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

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

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

R可简化表示为： $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],$ r₁～r₉表示R矩阵9个元素；

2、建立星敏感器weng成像模型；

设星光方向矢量为： $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₃表示星光向量在星敏感器坐标系下的x，y，z坐标值；α，β为恒星的赤经、赤纬坐标；

2.1、使用星敏感器对星光进行径向畸变成像；如图2所示为星敏感器径向畸变成像示意图。O-XYZ为星敏感器坐标系，o-xy为2维靶面坐标系，∏为靶面，Oo之间距离为焦距；星光向量w理想透视成像点为P′，由于发生径向畸变，实际成像点为P；

2.2、建立线性成像模型；由于x和y方向之间存在不确定比例因子，用f_u和f_v分别表示x和y方向的焦距，建立线性成像模型为：

$\frac{x^{\′}}{f_u}=\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}$

                                             [3]

$\frac{y^{\′}}{f_v}=\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}$

式中，x′，y′为透视投影的星点位置，定义线性参数向量d₁＝(f_u f_v α₀ β₀ ₀)；

2.3、考虑到畸变影响，建立weng非线性畸变模型公式为：

$\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}=\hat{u}+(g_1+g_3){\hat{u}}^2+g_4\hat{u}\hat{v}+g_1{\hat{v}}^2\text{+}{k}_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

                                            [4]

$\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}=\hat{v}+g_2{\hat{u}}^2+g_3\hat{u}\hat{v}+(g_2+g_4){\hat{v}}^2\text{+}{k}_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

式中， $\hat{u}=\frac{x}{f_u},$ $\hat{v}=\frac{y}{f_v},$ x，y为星点实际测量坐标，k₁为径向畸变系数，g₁～g₄为偏心畸变和薄棱镜畸变的综合公式系数，于是非线性的畸变参数有5个，定义参数向量d₂＝(k₁ g₁ g₂ g₃ g₄)；

3、数据采集；根据星敏感器采集得到的一帧星图，定义采集得到的全部星点数据集合为Ω₀，星点个数为n；以主点为中心，直径为图像尺寸1/2做圆，记圆内的数据集合为Ω₁，星点个数为n₁；

4、通过参数计算进行校准；

图3为参数计算过程示意图。校准过程分为两步完成，根据中心数据Ω₁受畸变影响较小的原理，第一步利用Ω₁构造线性计算模型，计算得到参数向量d₁的初始值；第二步对d₁的参数进行优化处理；第三步，根据weng模型，构造畸变系数线性方程，在全局数据Ω₀下计算参数向量d₂；最后，判断数据结果是否满足要求，如果不满足则返回第二步进行迭代计算，如果满足则输出参数向量的最终解d₁和d₂；具体计算过程如下：

4.1、进行线性参数估计；根据公式[3]有：

x′w₁r₇+x′w₂r₈+x′w₃r₉＝w₁r₁f_u+w₂r₂f_u+w₃r₃f_u

                                              [5]

y′w₁r₇+y′w₂r₈+y′w₃r₉＝w₁r₄f_v+w₂r₅f_v+w₃r₆f_v整理有：

$w_1\frac{r_1{f}_u}{r_9}+w_2\frac{r_2{f}_u}{r_9}+w_3\frac{r_3{f}_u}{r_9}-x{w}_1\frac{r_7}{r_9}-x{w}_2\frac{r_8}{r_9}=x{w}_3$

                                      [6]

$w_1\frac{r_4{f}_v}{r_9}+w_2\frac{r_5{f}_v}{r_9}+w_3\frac{r_6{f}_v}{r_9}-{\mathrm{yw}}_1\frac{r_8}{r_9}-{\mathrm{rw}}_2\frac{r_8}{r_9}=y{w}_3$

假设 $r_1^{\′}=\frac{r_1{f}_u}{r_9},$ $r_2^{\′}=\frac{r_2{f}_u}{r_9},$ $r_3^{\′}=\frac{r_3{f}_u}{r_9},$ $r_4^{\′}=\frac{r_4{f}_v}{r_9},$ $r_5^{\′}=\frac{r_5{f}_u}{r_9},$ $r_6^{\′}=\frac{r_6{f}_v}{r_9},$ $r_7^{\′}=\frac{r_7}{r_9},$ $r_8^{\′}=\frac{r_8}{r_9},$ 共8个未知数，定义参数向量m＝[r₁′r₂′r₃′r₄′r₅′r₆′r₇′r₈′]^T，对于每个星点数据可以得到方程：

[w₁ w₂ w₃ 0 0 0-xw₁-xw₂]m＝xw₃

                               [7]

[0 0 0 w₁ w₂ w₃-yw₁-yw₂]m＝yw₃

对于数据集合Ω₁内的星点，可以得到2n₁个方程，构成“Am＝b”形式；A为2n₁个方程左式构成的矩阵，c为2n₁个方程右式构成的向量；可用最小二乘法使得 $\underset{m}{\min }\left|\right|A\hat{m}-b\left|\right|,$ 解得8个未知数r₁′～r₈′；然后根据R为正交矩阵特点，有：

$r_9=\sqrt{1+r_7^{\′}+r_8^{\′}},$ r₇＝r₇′·r₉，  r₈＝r₈′·r₉

$f_u=r_9\·\sqrt{{r_1^{\′}}^2+{r_2^{\′}}^2+{r_3^{\′}}^2}$

r₁＝r₁′·r₉/f_u，  r₂＝r₂′·r₉/f_u，  r₃＝r₃′·r₉/f_u    [8]

$f_v=r_9\·\sqrt{{r_4^{\′}}^2+{r_5^{\′}}^2+{r_6^{\′}}^2}$

r₄＝r₄′·r₉/f_v，  r₅＝r₅′·r₉/f_v，  r₆＝r₆′·r₉/f_v

先假设r9为正，得到姿态估计矩阵，然后选择视场内的某一颗距离中心较远的恒星来判断r9是否真的为正，如果星光向量经姿态估计矩阵成像后坐标和测量坐标的符合一致，则说明r9为正，否则为负；

4.2、优化参数向量d₁；

设线性参数向量和非线性参数向量的真实值为d₁，d₂，估计值为

最优的估计值应该满足：

$\underset{d}{\min }\left|\right|P-f({\hat{d}}_1,{\hat{d}}_2)\left|\right|$

式中，P为点集Ω₁的星点坐标；以公式[8]得到的解为初值，采用非线性优化方法对d₁进行优化；设：

$x^{\′}=f_x(d_1,{\stackrel{\‾}{d}}_2)$

                                  [9]

$y^{\′}=f_y(d_1,{\stackrel{\‾}{d}}_2)$

式中，

表示参数向量d₂固定；由于f_x和f_y均为非线性函数，因此采用非线性线性化方法来估计参数向量d₁；假设

为估计值，

为测量值，Δd₁为向量估计偏差；

$\mathrm{\Δx}=\stackrel{\‾}{x}-\hat{x}\≈\mathrm{A\Δ}{d}_1$

                                      [10]

$\mathrm{\Δy}=\stackrel{\‾}{y}-\hat{y}\≈\mathrm{B\Δ}{d}_1$

式中，A、B为敏感矩阵，求各项偏导数即可获得；

于是得到参数向量更新值d₁＝d₁+Δd₁；

4.3、计算畸变系数向量d₂；

得到线性参数后，下一步求解非线性参数d₂；根据公式[4]，构造线性估计模型：

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& {\hat{u}}^2+{\hat{v}}^2& 0& {\hat{u}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=\mathrm{du}$

                                    [11]

$\left(\begin{array}{ccccc}\hat{u}({\hat{u}}^2+{\hat{v}}^2)& 0& {\hat{u}}^2+{\hat{v}}^2& {\hat{v}}^2& \hat{u}\hat{v}\end{array}\right){\hat{d}}_2=\mathrm{dv}$

$\mathrm{du}=\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}-\hat{u}$

式中：

$\mathrm{dv}=\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}-\hat{v}$

构成“Bd₂＝c”形式，B为2n个方程左式构成的矩阵，c为2n个方程右式构成的向量；利用最小二乘法使得 $\underset{d_2}{\min }\left|\right|B{d}_2-c\left|\right|,$ 求出畸变参数向量d₂的估计值，这里的数据为全部点集Ω₀；

得到模型的全部参数d₁和d₂后，通过对比理想坐标值和测量坐标值进行判断；如果误差不满足要求，则回到步骤1.4.2进一步优化d₁，进行迭代计算；如果误差满足要求，则停止迭代，输出结果。

仿真和结果分析

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

视场：12×12；

像素阵列：1024×1024；

像素尺寸：0.015mm×0.015mm；

焦距：73.0703mm

基于weng畸变模型的在轨校准方法，对视场内的星数要求比较多。可以采用二维随机分布点来模拟恒星在视场的分布，这样做的好处是很方便产生一定星数要求的星场，并可以研究星数和校准参数收敛的稳定性。这里假设星场内有100颗星。

图4为模拟星场示意图，图中圆圈直径为9.2mm，即靶面面积的28.3％，圆内有28颗星。假设姿态角为(0，0，0)，靶面星点坐标均匀随机分布。

仿真参数为：

f_u＝76.7238mm，f_v＝73.0703mm。

k₁＝-0.3，g₁＝0.02，g₂＝-0.009，g₃＝-0.02，g₄＝0.009。

星点噪声水平为0.05像素方差的高斯白噪声。

第一步的闭式计算结果为：

fu/mm	fv/mm	α0	β0	φ0
					75.957	72.223	0.015	0.140	-0.00116

第二，三步的迭代计算过程为：

首先式参数向量d₁的计算过程见下表：

迭代次数	fu/mm	fv/mm	α0/arcsec	β0/arcsec	φ0/arcsec
						1	77.373	74.495	-67.2277	939.3451	239.1954
2	76.917	73.206	-15.8128	201.4925	473.6030
						3	76.754	73.081	-2.4443	41.4535	103.5382
4	76.727	73.061	-0.7597	8.0169	20.8955
						5	76.724	73.058	-0.5817	1.0852	3.6631
6	76.723	73.058	-0.5772	-0.3536	0.0849
						7	76.723	73.058	-0.5828	-0.6531	-0.6587

从上表中的姿态欧拉角偏差可以看出，得到的外部姿态角误差已经小于1角秒，完全可以满足工作要求。

参数向量d₂的计算过程见下表：

迭代次数	k1	g1	g2	g3	g4
						1	-1.1	0.035399	0.43	-0.12328	0.024027
2	-0.40897	0.027119	0.079602	-0.045578	0.016899
						3	-0.31141	0.022141	0.0089006	-0.025147	0.011051
4	-0.29845	0.020565	-0.005874	-0.020933	0.0098816

5	-0.29733	0.020158	-0.008952	-0.020081	0.0096507
						6	-0.29744	0.020059	-0.009594	-0.019909	0.0096042
7	-0.29753	0.020035	-0.009728	-0.019874	0.0095947

最终的星点坐标估计值和测量值的偏差均方根为x方向0.0571，y方向0.0575像素。由于噪声水平为0.05像素，因此最终的参数可以满足校准精度要求。

Claims (1)

1、基于weng模型的星敏感器在轨校准方法，其特征在于，校准的步骤如下：

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

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

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

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

式中，s表示正弦函数，c表示余弦函数；R可简化表示为： $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],$ r₁～r₉表示R矩阵9个元素；

1.2、建立星敏感器weng成像模型；

设星光方向矢量为： $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₃表示星光向量在星敏感器坐标系下的x，y，z坐标值；α，β为恒星的赤经、赤纬坐标；

1.2.1、使用星敏感器对星光进行径向畸变成像；0-XYZ为星敏感器坐标系，o-xy为2维靶面坐标系，П为靶面，0o之间距离为焦距；星光向量w理想透视成像点为P′，由于发生径向畸变，实际成像点为P；

1.2.2、建立线性成像模型；由于x和y方向之间存在不确定比例因子，用f_u和f_v分别表示x和y方向的焦距，建立线性成像模型为：

$\frac{x^{\′}}{f_u}=\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]$

$\frac{y^{\′}}{f_v}=\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}$

式中，x′，y′为透视投影的星点位置，定义线性参数向量d₁＝(f_u f_v α₀ β₀ ₀)；

1.2.3、考虑到畸变影响，建立weng非线性畸变模型公式为：

$\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}=\hat{u}+(g_1+g_3){\hat{u}}^2+g_4\hat{u}\hat{v}+g_1{\hat{v}}^2+k_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)---\left[4\right]$

$\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}=\hat{v}+g_2{\hat{u}}^2+g_3\hat{u}\hat{v}+(g_2+g_4){\hat{v}}^2{+k}_1\hat{u}({\hat{u}}^2+{\hat{v}}^2)$

式中， $\hat{u}=\frac{x}{f_u},$ $\hat{v}=\frac{y}{f_v},$ x，y为星点实际测量坐标，k₁为径向畸变系数，g₁～g₄为偏心畸变和薄棱镜畸变的综合公式系数，于是非线性的畸变参数有5个，定义参数向量d₂＝(k₁ g₁ g₂ g₃ g₄)；

1.3、数据采集；根据星敏感器采集得到的一帧星图，定义采集得到的全部星点数据集合为Ω₀，星点个数为n；以主点为中心，直径为图像尺寸1/2做圆，记圆内的数据集合为Ω₁，星点个数为n₁；

1.4、通过参数计算进行校准；

校准过程分为两步完成，根据中心数据Ω₁受畸变影响较小的原理，第一步利用Ω₁构造线性计算模型，计算得到参数向量d₁的初始值；第二步对d₁的参数进行优化处理；第三步，根据weng模型，构造畸变系数线性方程，在全局数据Ω₀下计算参数向量d₂；最后，判断数据结果是否满足要求，如果不满足则返回第二步进行迭代计算，如果满足则输出参数向量的最终解d₁和d₂；具体计算过程如下：

1.4.1、进行线性参数估计；根据公式[3]有：

x′w₁r₇+x′w₂r₈+x′w₃r₉＝w₁r₁f_u+w₂r₂f_u+w₃r₃r_u

y′w₁r₇+y′w₂r₈+y′w₃r₉＝w₁r₄f_v+w₂r₅f_v+w₃r₆f_v     [5]

整理有：

$w_1\frac{r_1{f}_u}{r_9}+w_2\frac{r_2{f}_u}{r_9}{w}_3\frac{r_3{f}_u}{r_9}-x{w}_1\frac{r_7}{r_9}-x{w}_2\frac{r_8}{r_9}=x{w}_3---\left[6\right]$

$w_1\frac{r_4{f}_v}{r_9}+w_2\frac{r_5{f}_v}{r_9}+w_3\frac{r_6{f}_v}{r_9}-y{w}_1\frac{r_8}{r_9}-y{w}_2\frac{r_8}{r_9}=y{w}_3$

假设 $r_1^{\′}=\frac{r_1{f}_u}{r_9},$ $r_2^{\′}=\frac{r_2{f}_u}{r_9},$ $r_3^{\′}=\frac{r_3{f}_u}{r_9},$ $r_4^{\′}=\frac{r_4{f}_v}{r_9},$ $r_5^{\′}=\frac{r_5{f}_u}{r_9},$ $r_6^{\′}=\frac{r_6{f}_v}{r_9},$ $r_7^{\′}=\frac{r_7}{r_9},$ $r_8^{\′}=\frac{r_8}{r_9},$ 共8个未知数，定义参数向量m＝[r₁′ r₂′ r₃′ r₄′ r₅′ r₆′ r₇′r₈′]^T，对于每个星点数据可以得到方程：

[w₁ w₂ w₃ 0 0 0-xw₁-xw₂]m＝xw₃

[0 0 0 w₁ w₂ w₃-yw₁-yw₂]m＝yw₃   [7]

对于数据集合Ω₁内的星点，可以得到2n₁个方程，构成“Am＝b”形式；A为2n₁个方程左式构成的矩阵，b为2n₁个方程右式构成的向量；可用最小二乘法使得 $\underset{m}{\min }\left|\right|A\hat{m}-b\left|\right|,$ 解得8个未知数r₁′～r₈′；然后根据R为正交矩阵特点，有：

$r_9=\sqrt{1+r_7^{\′}+r_8^{\′}},$ r₇＝r₇′·r₉，r₈＝r₈′·r₉

$f_u=r_9\·\sqrt{{r_1^{\′}}^2+{r_2^{\′}}^2+{r_3^{\′}}^2}$

r₁＝r₁′·r₉/f_u，r₂＝r₂′·r₉/f_u，r₃＝r₃′·r₉/f_u

$f_v=r_9\·\sqrt{{r_4^{\′}}^2+{r_5^{\′}}^2+{r_6^{\′}}^2}$

r₄＝r₄′·r₉/f_v，r₅＝r₅′·r₉/f_v ,r₆＝r₆′·r₉/f_v

先假设r₉为正，得到姿态估计矩阵，然后选择视场内的某一颗距离中心较远的恒星来判断r₉是否真的为正，如果星光向量经姿态估计矩阵成像后坐标和测量坐标的符合一致，则说明r₉为正，否则为负；

1.4.2、优化参数向量d₁；

设线性参数向量和非线性参数向量的真实值为d₁，d₂，估计值为

最优的估计值应该满足：

$\underset{d}{\min }\left|\right|P-f({\hat{d}}_1,{\hat{d}}_2)\left|\right|$

式中，P为点集Ω₁的星点坐标；以公式[8]得到的解为初值，采用非线性优化方法对d₁进行优化；设：

$x^{\′}=f_x(d_1,{\stackrel{\‾}{d}}_2)---\left[9\right]$

$y^{\′}=f_y(d_1,{\stackrel{\‾}{d}}_2)$

式中，

表示参数向量d₂固定；由于f_x和f_y均为非线性函数，因此采用非线性线性化方法来估计参数向量d₁；假设

为估计值，

为测量值，Δd₁为向量估计偏差；

$\mathrm{\Δx}=\stackrel{\‾}{x}-\hat{x}\≈\mathrm{A\Δ}{d}_1---\left[10\right]$

$\mathrm{\Δy}=\stackrel{\‾}{y}-\hat{y}\≈\mathrm{B\Δ}{d}_1$

式中，A、B为敏感矩阵，求各项偏导数即可获得；

于是得到参数向量更新值d₁＝d₁+Δd₁；

1.4.3、计算畸变系数向量d₂；

得到线性参数后，下一步求解非线性参数d₂；根据公式[4]，构造线性估计模型：

$(\hat{u}({\hat{u}}^2+{\hat{v}}^2){\hat{u}}^2+{\hat{v}}^{2}0{\hat{u}}^2\hat{u}\hat{v}){\hat{d}}_2=\mathrm{du}---\left[11\right]$

$(\hat{u}({\hat{u}}^2+{\hat{v}}^2)0{\hat{u}}^2+{\hat{v}}^2{\hat{v}}^2\hat{u}\hat{v}){\hat{d}}_2=\mathrm{dv}$

式中： $\mathrm{du}=\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}-\hat{u}$

$\mathrm{dv}=\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}-\hat{v}$

构成“Bd₂＝c”形式，B为2n个方程左式构成的矩阵，b为2n个方程右式构成的向量；利用最小二乘法使得 $\underset{d_2}{\min }\left|\right|{\mathrm{Bd}}_2-c\left|\right|,$ 求出畸变参数向量d₂的估计值，这里的数据为全部点集Ω₀；

得到模型的全部参数d₁和d₂后，通过对比理想坐标值和测量坐标值进行判断；如果误差不满足要求，则回到步骤1.4.2进一步优化d₁，进行迭代计算；如果误差满足要求，则停止迭代，输出结果。

Images