CN100349018C - 一种星敏感器内外方元素校准方法 - Google Patents

Abstract

本发明航天测量技术，涉及对星敏感器校准方法的改进。它基于一个由气垫平台、单星星光模拟器、星敏感器、2维轴向转台以及数据处理计算机组成的校准系统，其步骤是：校准系统建模；采集数据；数据处理。本发明采用了整体建模的方法，避免了因外部参数偏差引入到内部参数的估计过程中而影响参数校准精度。该方法减少了校准过程中复杂的星敏感器安装对准过程，使得整个校准过程变得相对简单。

Description

一种星敏感器内外方元素校准方法

技术领域

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

背景技术

星敏感器是一种利用恒星观测，为空间飞行器提供高精度姿态信息的航天测量仪器。其工作原理为：星敏感器前端摄像头单元利用CCD(或CMOS)图像传感器拍摄得到星图图像，经过图像处理程序得到恒星像点的质心坐标和亮度的信息，然后星图识别程序利用这些信息在导航星库中找到对应的恒星，最后计算出星敏感器的三轴姿态。在星敏感器投入使用前，其主点、焦距以及畸变系数等内部参数必须进行精确的测量，称为星敏感器校准。常见的星敏感器校准有几种方法，第一种是在星光实验室利用星光模拟器配合高精度2轴转台进行数据采集和校准；第二种是在大气环境良好的地点对晴朗夜空进行拍摄以获得数据并进行校准；第三种是当星敏感器在轨工作时进行校准。其中在星光实验室的校准是最基础和精度最高的校准。如图1所示，该校准系统主要有气垫平台，单星星光模拟器，星敏感器，2维轴向转台以及数据处理计算机组成。目前的校准方法为：

(1)安装星敏感器到2维轴向转台上，使得星敏感器的视轴垂直于2维轴向转台的两个旋转轴。

(2)调整单星星光模拟器，使得其发生的星光方向和星敏感器的视轴对齐。

(3)然后根据图2所示，连续转动转台采集星点数据，并记录当时的转台转动坐标，最终使得星点数据遍布整个星敏感器图像传感器靶面。

(4)根据小孔成像原理，建立星敏感器星点成像模型，模型参数包括主点、焦距以及畸变系数等，利用采集得到数据计算星敏感器模型参数。

以上方法存在的问题是：

(1)星敏感器在转台上的安装存在偏差，使得星敏感器的视轴不能够和2维轴向转台的两个旋转轴恰好垂直。

(2)调整星光模拟器和转台坐标系两个旋转轴组成的平面不垂直。

由于这些安装对齐等外部因素影响，使得仅仅对星敏感器内部参数建立成像模型的校准方法存在误差，从而最终影响星敏感器内部参数的估计精度。

发明内容

本发明的目的是：针对上述校准方法存在的问题，提出一种对星敏感器校准系统内部参数和外部参数联合建模的方法，利用非线性最小二乘法和共线性公式，迭代计算出星敏感器校准系统内部参数和外部参数。

本发明的技术方案是：一种星敏感器内外方元素校准方法，基于一个由气垫平台、单星星光模拟器、星敏感器、2维轴向转台以及数据处理计算机组成的校准系统，单星星光模拟器、星敏感器、2维轴向转台安装在气垫平台上，星敏感器安装在2维轴向转台的内框上，星敏感器的视轴垂直于2维轴向转台的两个旋转轴，单星星光模拟器发生的星光方向和星敏感器的视轴对齐，其特征在于，校准的步骤如下：

1、校准系统建模；

1.1、外部参数建模；影响星敏感器采集得到的星点坐标位置的外部因素有：星光模拟器发生的星光方向，转台内外两个转动框偏离转台初始位置的角度，星敏感器在转台内框上的安装偏差，以及星敏感器图像传感器成像靶面和机壳之间的安装偏差；通过建立多个坐标系，将上述外部参数联系起来，以便分析星敏感器靶面上的星点成像位置；外部参数建模的步骤如下：

1.1.1、根据转台的初始位置建立坐标系N；

坐标系N以内框的转动轴为Xn轴，外框转动轴为Yn轴，Xn轴和Yn轴交于星敏感器镜头中心点O点，O为坐标系原点，过O点沿着星敏感器视轴方向为坐标系的Zn轴；星光模拟器发生的星光方向矢量Vn在坐标系N中的表达式为：

$v_n=\left[\begin{array}{c}{n}_1\\ n_2\\ n_3\end{array}\right]=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[1\right]$

式中，n₁、n₂、n₃为星光矢量Vn在坐标系N中的3个方向分量，α和β分别为该矢量在坐标系N中的偏航和俯仰角；

1.1.2、建立坐标系B；

当外框转过一个角度θ₁后，根据转动后的转台状态建立坐标系B，坐标系B的原点和坐标系N相同，其Yb轴和Yn轴重合，Xb轴、Zb轴分别位于Xn轴、Zn轴转动后的位置，当外框转过θ₁角后，星光矢量在坐标系B中的表达式为：

$v_b=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]=R_{\mathrm{bn}}{v}_n=\left[\begin{array}{ccc}\cos {\mathrm{\θ}}_1& 0& -{\sin \mathrm{\θ}}_1\\ 0& 1& 0\\ \sin {\mathrm{\θ}}_1& 0& \cos {\mathrm{\θ}}_1\end{array}\right]\left[\begin{array}{c}{n}_1\\ n_2\\ n_3\end{array}\right]---\left[2\right]$

式中，b₁、b₂、b₃为星光矢量在坐标系B中的3个方向分量，R_bn为转换矩阵；

1.1.3、建立坐标系C；

当内框转过一个角度θ₂后，根据转动后的转台状态建立坐标系C；坐标系C的原点和坐标系B相同，其Xc轴和Xb轴重合，Yc轴、Zc轴分别位于Yb轴、Zb轴转动后的位置；当外框转过θ₂角后，星光矢量在坐标系C中的表达式为：

$v_c=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=R_{\mathrm{cb}}{v}_b=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\mathrm{\θ}}_2& \sin {\mathrm{\θ}}_2\\ 0& -\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]---\left[3\right]$

式中，c₁、c₂、c₃为星光矢量在坐标系C中的3个方向分量，R_cb为转换矩阵；

1.1.4、建立坐标系D和E；

星敏感器在转台内框上的安装偏角为₁和₂，₁为偏航方向安装误差，₂为俯仰方向安装误差；对应于星敏感器这两个安装位置，分别建立坐标系D和E；

坐标系D的原点和坐标系C相同，其Xd轴和Xc轴重合，Yd轴、Zd轴分别位于Yc轴、Zc轴转动后的位置，当外框转过₁角后，星光矢量在坐标系D中的表达式为：

式中，d₁、d₂、d₃为星光矢量在坐标系D中的3个方向分量，R_dc为转换矩阵；

坐标系E的原点和坐标系D相同，其Xe轴和Xd轴重合，Ye轴、Ze轴分别位于Yd轴、Zd轴转动后的位置，当外框转过₂角后，星光矢量在坐标系E中的表达式为：

式中，e₁、e₂、e₃为星光矢量在坐标系E中的3个方向分量，R_ed为转换矩阵；坐标系E即为星敏感器机壳坐标系；

1.1.5、建立坐标系F；

图像传感器靶面的X、Y轴和机壳坐标系E的Xe，Ye轴之间的安装偏差造成的偏置角为₃，根据图像传感器靶面建立星敏感器的针孔成像坐标系F；坐标系F和坐标系E原点相同，Zf轴和Ze轴相重合，Xf轴和靶面图像传感器靶面的X方向一致，Yf轴和靶面Y方向一致；靶面坐标系与图像传感器采集得到的图像一致，原点位于图像左下角，X、Y方向即为图像的横、纵方向，星光矢量从坐标系E到坐标系F的转换过程为：

式中，f₁，f₂，f₃为星光矢量在坐标系F中的3个方向分量，R_fe为转换矩阵；

根据以上坐标系的转换过程，得到从转台坐标系N到星敏感器坐标系F的转换关系为：

$v_f=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=R_{\mathrm{fe}}{R}_{\mathrm{ed}}{R}_{\mathrm{dc}}{R}_{\mathrm{cb}}{R}_{\mathrm{bn}}{v}_n---\left[7\right]$

展开式[7]得到：

f₁＝cos₃cos₂cosθ₁cosβcosα+sin₃sin(θ₂+₁)sinθ₁cosβcosα

-cos₃sin₂cos(θ₂+₁)sinθ₁cosβcosα

+sin₃cos(θ₂+₁)cosβsinα+cos₃sin₂sin(θ₂+₁)cosβsinα

-cos₃cos₂sinθ₁sinβ+sin₃sin(θ₂+₁)cosθ₁sinβ

-cos₃sin₂cos(θ₂+₁)cosθ₁sinβ

f₂＝-sin₃cos₂cosθ₁cosβcosα+cos₃sin(θ₂+₁)sinθ₁cosβcosα

+sin₃sin₂cos(θ₂+₁)sinθ₁cosβcosα

+cos₃cos(θ₂+₁)cosβsinα-sin₃sin₃sin(θ₂+₁)cosβsinα

+sin₃cos(₂)sinθ₁sinβ+cos₃sin(θ₂+₁)cosθ₁sinβ

+sin₃sin₂cos(θ₂+₁)cosθ₁sinβ

f₃＝sin₂cosθ₁cosβcosα+cos₂cos(θ₂+₁)sinθ₁cosβcosα

-cosθ₂sin(θ₂+₁)cosβsinα-sin₂sinθ₁sinβ+cos₂cos(θ₂+₁)cosθ₁sinβ

上述外部参数中，转台的偏航角和俯仰角由转台自身提供，需要求取的外部参数有5个，即为α、β、φ₁、φ₂、φ₃；

1.2、针孔成像建模；

根据针孔成像原理，在坐标系F中，Zf为星敏感器视轴方向，Xf，Yf和星敏感器靶面坐标系的X，Y轴方向一致，这里的靶面坐标系为2维平面坐标系，以主点O’为原点，其X，Y轴对应于采集得到图像的横纵轴，星光经过镜头中心点O后，投影成像在星敏感器靶面P(x，y)点；假设星敏感器的焦距为f_c，主点O’的图像坐标为(x₀，y₀)，根据共线性公式有：

$x=f_c\frac{f_1}{f_3}+x_0,$ $y=f_c\frac{f_2}{f_3}+y_0---\left[8\right]$

式中，主点(x₀，y₀)和焦距f_c为3个未知参数；

1.3、镜头畸变建模；

dx，dy为星敏感器在x方向和y方向的畸变偏差，有：

$\mathrm{dy}=\stackrel{\‾}{y}[q_1{r}^2+q_2{r}^4+q_3{r}^6]+[p_2(r^2+2{\stackrel{\‾}{y}}^2)+2p_1\stackrel{\‾}{\mathrm{xy}}][1+p_3{r}^2]---\left[9\right]$

式中，

$\left\{\begin{array}{c}\stackrel{\‾}{x}=x-x_0\\ \stackrel{\‾}{y}=y-y_0\end{array}\right.;$

$r^2={\stackrel{\‾}{x}}^2+{\stackrel{\‾}{y}}^2;$

q₁、q₂、q₃为径向畸变系数；

p₁、p₂、p₃为偏心畸变系数；

于是内部参数总共有9个，即为x₀、y₀、f_c、q₁、q₂、q₃、p₁、p₂、p₃，根据公式[8]有：

$x=f_c\frac{f_1}{f_3}+x_0+\mathrm{dx}$

$y=f_c\frac{f_2}{f_3}+y_0+\mathrm{dy}---\left[10\right]$

上述模型中，总共有外部参数和内部参数14个，即：α、β、φ₁、φ₂、φ₃、x₀、y₀、f_c、q₁、q₂、q₃、p₁、p₂、p₃；

2、采集数据；

转动转台，从-6度到+6度，每间隔1度为一个采集位置，在一个位置采集n次数据，n可取100～1000，并记录当时的转台转动坐标，最终使得星点成像遍布靶面。公式为：

$\tilde{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

$\stackrel{\‾}{y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i$

这里x_i、y_i为每次采集得到星点成像坐标，

为n次采集后的平均坐标；

3、数据处理；

主点位置采用自准直方法，星光垂直投射到靶面上，求出此时的图像坐标；此时，需要校准的参数有12个，用一个向量 来表示这些向量有：

根据公式进一步可以得到：

$\stackrel{\‾}{x}=x-x_0=f_c\frac{f_1}{f_3}+\mathrm{dx}=f_x\left(\stackrel{\→}{x}\right)---\left[11\right]$

$\stackrel{\‾}{y}=y-y_0=f_c\frac{f_2}{f_3}+\mathrm{dy}=f_y\left(\stackrel{\→}{x}\right)$

由于f_x和f_y均为非线性函数，因此采用非线性最小二乘迭代方法来估计参数向量

为实际测量得到的

的估计值，

为向量估计偏差；则有：

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

$\mathrm{\Δy}=\stackrel{\‾}{y}-\hat{y}\≈\mathrm{B\Δ}\stackrel{\→}{x}$

这里A、B为敏感矩阵；

星点采集个数为m，联合x和y方向的偏差和敏感矩阵，

$p=\left(\begin{array}{c}{\mathrm{\Δx}}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{x}_m\\ \mathrm{\Δ}{y}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{y}_m\end{array}\right),$ $M=\left(\begin{array}{c}{A}_1\\ .\\ .\\ .\\ A_m\\ B_1\\ .\\ .\\ .\\ B_m\end{array}\right)$

这里p是由x和y方向残留偏差组成的向量，M为A和B两个敏感矩阵组成的整体敏感矩阵；

于是有迭代方程为：

$\mathrm{\Δ}{\stackrel{\→}{x}}^{(k+1)}=\mathrm{\Δ}{\stackrel{\→}{x}}^{\left(k\right)}-{\left(M_k^T{M}_k\right)}^{-1}{M}_k^T{p}^{\left(k\right)}---\left[12\right]$

这里k为迭代序号，k取5～20，迭代结束后得到稳定的数据值，这时的参数即为最后的校准结果。

本发明的优点是：第一，由于采用了整体建模的方法，将星敏感器校准系统的内部参数和外部参数同时进行估计，外部参数对应于世界坐标系，求出了外部参数的估计偏差后，可以避免外部参数偏差引入到内部参数的估计过程中，影响参数校准精度。第二，该方法减少了校准过程中复杂的星敏感器安装对准过程，使得整个校准过程变得相对简单。第三，建立了世界坐标系后，可以很方便的将星敏感器坐标系引出到和机壳固联的光学立方体坐标系上，为将来星敏感器安装到空间飞行器上时的方位测量提供了基础。

附图说明

图1是现有的星敏感器校准装置构成示意图。

图2是现有的星敏感器校准装置在校准时星点采集过程示意图。

图3是建立转台初始位置坐标系N示意图。

图4是转台坐标系N中的星向量Vn示意图。

图5是由坐标系N到坐标系B的转换示意图。

图6是由坐标系B到坐标系C的转换示意图。

图7是由坐标系C到坐标系D的转换示意图。

图8是由坐标系D到坐标系E的转换示意图。

图9是由坐标系E到坐标系F的转换示意图。

图10是星点针孔成像前视示意图。

具体实施方式

下面对本发明做进一步详细说明。本发明提出一种对星敏感器校准系统内部参数和外部参数联合建模的方法，并利用非线性最小二乘法和共线性公式，迭代计算出星敏感器校准系统内部参数和外部参数。这种方法能够将星敏感器校准系统的内部参数和外部参数分离开来，使得最终获得的星敏感器内部参数精度更高。同时，该校准过程中建立了星敏感器校准系统的外部世界坐标系，从而很方便地将星敏感器坐标系引出到和星敏感器固联的光学立方体坐标系上，为将来星敏感器安装到空间飞行器上时的方位测量提供了基础。下面详细说明本发明方法的步骤。

1、系统建模；

1.1、外部参数建模；

星敏感器测试系统在进行数据采集过程中，影响星敏感器采集得到的星点坐标位置的外部因素有：星光模拟器发生的星光方向，转台内外两个转动框偏离转台初始位置的角度，星敏感器在转台内框上的安装偏差，以及星敏感器图像传感器成像靶面和机壳之间的安装偏差。于是需要建立多个坐标系，将这些外部参数联系起来，以便分析星敏感器靶面上的星点成像位置。

1.1.1、首先根据转台的初始位置建立坐标系N，如图3所示。坐标系N以内框的转动轴为Xn轴，外框转动轴为Yn轴，Xn轴和Yn轴交星敏感器镜头中心点O点，O为坐标系原点。过O点沿着星敏感器视轴方向为坐标系的Zn轴。星光模拟器发生的星光方向矢量Vn在坐标系N中的表达式为：

$v_n=\left[\begin{array}{c}{n}_1\\ n_2\\ n_3\end{array}\right]=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[1\right]$

其中，n₁，n₂，n₃为星光矢量Vn在坐标系N中的3个方向分量，α和β分别为该矢量在坐标系N中的偏航和俯仰角，如图4所示。

1.1.2、当外框转过一个角度θ₁后，根据转动后的转台状态建立坐标系B，如图5所示。坐标系B的原点和坐标系N相同，其Yb轴和Yn轴重合。Xb轴、Zb轴分别位于Xn轴、Zn轴转动后的位置。当外框转过θ₁角后，星光矢量在坐标系B中的表达式为：

$v_b=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]=R_{\mathrm{bn}}{v}_n=\left[\begin{array}{ccc}\cos {\mathrm{\θ}}_1& 0& -\sin {\mathrm{\θ}}_1\\ 0& 1& 0\\ \sin {\mathrm{\θ}}_1& 0& \cos {\mathrm{\θ}}_1\end{array}\right]\left[\begin{array}{c}{n}_1\\ n_2\\ n_3\end{array}\right]---\left[2\right]$

这里，b₁，b₂，b₃为星光矢量在坐标系B中的3个方向分量，R_bn为转换矩阵。

1.1.3、当内框转过一个角度θ₂后，根据转动后的转台状态建立坐标系C，如图6所示。坐标系C的原点和坐标系B相同，其Xc轴和Xb轴重合。Yc轴、Zc轴分别位于Yb轴、Zb轴转动后的位置。当外框转过θ₂角后，星光矢量在坐标系C中的表达式为：

$v_c=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=R_{\mathrm{cb}}{v}_b=\left[\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_2& \sin {\mathrm{\θ}}_2\\ 0& -\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]---\left[3\right]$

这里c₁，c₂，c₃为星光矢量在坐标系C中的3个方向分量，R_cb为转换矩阵。

1.1.4、星敏感器安装在转台内框上，这里假设星敏感器的安装偏角为₁和₂，分别对应于偏航方向误差和俯仰方向安装误差。分别建立坐标系D和E，对应于星敏感器这两个安装位置，分别如图7和图8所示。这里坐标系E即为星敏感器机壳坐标系。坐标系D的原点和坐标系C相同，其Xd轴和Xc轴重合。Yd轴、Zd轴分别位于Yc轴、Zc轴转动后的位置。当外框转过₁角后，星光矢量在坐标系D中的表达式为：

这里d₁，d₂，d₃为星光矢量在坐标系D中的3个方向分量，R_dc为转换矩阵。

坐标系E的原点和坐标系D相同，其Xe轴和Xd轴重合。Ye轴、Ze轴分别位于Yd轴、Zd轴转动后的位置。当外框转过₂角后，星光矢量在坐标系E中的表达式为：

这里e₁， e₂，e₃为星光矢量在坐标系E中的3个方向分量，R_ed为转换矩阵。

1.1.5、在星敏感器装配过程中，图像传感器靶面的X，Y轴和机壳坐标系E的Xe，Ye轴之间可能出现偏差，设此偏置角为₃。于是，根据图像传感器靶面建立星敏感器的针孔成像坐标系F，如图9所示。坐标系F和坐标系E原点相同，Zf轴和Ze轴相重合。Xf轴和靶面图像传感器靶面的X方向一致，Yf轴和靶面Y方向一致。靶面坐标系与图像传感器采集得到的图像一致，原点位于图像左下角，X，Y方向即为图像的横，纵方向。于是星光矢量从坐标系E到坐标系F的转换过程为：

这里f₁，f₂，f₃为星光矢量在坐标系F中的3个方向分量，R_fe为转换矩阵。

总结以上坐标系的转换过程，可以得到从转台坐标系N到星敏感器坐标系F的转换关系，即：

$v_f=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=R_{\mathrm{fe}}{R}_{\mathrm{ed}}{R}_{\mathrm{dc}}{R}_{\mathrm{cb}}{R}_{\mathrm{bn}}{v}_n---\left[7\right]$

展开该公式可以得到：

f₁＝cos₃cos₂cosθ₁cosβcosα+sin₃sin(θ₂+₁)sinθ₁cosβcosα

一cos₃sin₂cos(θ₂+₁)sinθ₁cosβcosα

+sin₃cos(θ₂+₁)cosβsinα+cos₃sin₂sin(θ₂+₁)cosβsinα

-cos₃cos₂sinθ₁sinβ+sin₃sin(θ₂+₁)cosθ₁sinβ

-cos₃sin₂cos(θ₂+₁)cosθ₁sinβ

f₂＝-sin₃cos₂cosθ₁cosβcosα+cos₃sin(θ₂+₁)sinθ₁cosβcosα

+sin₃sin₂cos(θ₂+₁)sinθ₁cosβcosα

+cos₃cos(θ₂+₁)cosβsinα-sin₃sin₂sin(θ₂+₁)cosβsinα

+sin₃cos(₂)sinθ₁sinβ+cos₃sin(θ₂+₁)cosθ₁sinβ

+sin₃sin₂cos(θ₂+₁)cosθ₁sinβ

f₃＝sin₂cosθ₁cosβcosα+cos₂cos(θ₂+₁)sinθ₁cosβcosα

-cosθ₂sin(θ₂+₁)cosβsinα-sin₂sinθ₁sinβ+cos₂cos(θ₂+₁)cosθ₁sinβ

考虑到转台的偏航角和俯仰角为转台自身提供，所以这里需要求取的外部参数实际有5个，即为(α，β，φ₁，φ₂，φ₃)。

1.2、针孔成像建模；

根据针孔成像原理，坐标系F中，Zf为星敏感器视轴方向，Xf，Yf和星敏感器靶面坐标系的X，Y轴方向一致。这里的靶面坐标系为2维平面坐标系，以主点O’为原点，其X，Y轴对应于采集得到图像的横纵轴。如图10所示，星光经过镜头中心点O后，投影成像在星敏感器靶面P(x，y)点。假设星敏感器的焦距为f_c，主点O’的图像坐标为(x₀，y₀)。根据共线性公式有：

$x=f_c\frac{f_1}{f_3}+x_0,$ $y=f_c\frac{f_2}{f_3}+y_0---\left[8\right]$

这里，主点(x₀，y₀)和焦距f_c为3个未知参数。

1.3、镜头畸变建模；

由于制造、安装的工艺限制，使得星敏感器镜头实际工作状态并不完全符合针孔成像透视原理，而总是会存在一点的畸变，下面对星敏感器的畸变建立数学模型。假设dx，dy为x方向和y方向的畸变偏差，有：

$\mathrm{dy}=\stackrel{\‾}{y}[q_1{r}^2+q_2{r}^4+q_3{r}^6]+[p_2(r^2+2{\stackrel{\‾}{y}}^2)+2p_1\stackrel{\‾}{\mathrm{xy}}][1+p_3{r}^2]---\left[9\right]$

这里，

$\left\{\begin{array}{c}\stackrel{\‾}{x}=x-x_0\\ \stackrel{\‾}{y}=y-y_0\end{array}\right.;$

$r^2={\stackrel{\‾}{x}}^2+{\stackrel{\‾}{y}}^2;$

q₁，q₂，q₃为径向畸变系数；

p₁，p₂，p₃为偏心畸变系数。

于是内部参数总共有9个，即为(x₀，y₀，f_c，q₁，q₂，q₃，p₁，p₂，p₃)，根据公式[8]有：

$x=f_c\frac{f_1}{f_3}+x_0+\mathrm{dx}---\left[10\right]$

$y=f_c\frac{f_2}{f_3}+y_0+\mathrm{dy}$

考虑外部参数和内部参数，上述模型中总共有1 4个参数，即：(α，β，φ₁，φ₂，φ₃，x₀，y₀，f_c，q₁，q₂，q₃，p₁，p₂，p₃)。

2、数据采集；

转动转台，从-6度到+6度，每间隔1度为一个采集位置，在一个位置采集n次数据，n可取100～1000，并记录当时的转台转动坐标，最终使得星点成像遍布靶面。公式为：

$\tilde{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

$\tilde{y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i$

这里x_i，y_i为每次采集得到星点成像坐标，

为n次采集后的平均坐标；如果n＝100，则质心噪声水平可以下降一个数量级。

3、数据处理；

设主点位置采用自准直方法，星光垂直投射到靶面上，求出此时的图像坐标。之后，需要校准的参数还有12个，用一个向量

来表示这些向量有：

根据公式进一步可以得到：

$\stackrel{\‾}{x}=x-x_0=f_c\frac{f_1}{f_3}+\mathrm{dx}=f_x\left(\stackrel{\→}{x}\right)---\left[11\right]$

$\stackrel{\‾}{y}=y-y_0=f_c\frac{f_2}{f_3}+\mathrm{dy}=f_y\left(\stackrel{\→}{x}\right)$

由于f_x和f_y均为非线性函数，因此采用非线性最小二乘迭代方法来估计参数向量

假设

为实际测量得到的

的估计值，

为向量估计偏差。则有：

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

$\mathrm{\Δy}=\stackrel{\‾}{y}-\hat{y}\≈\mathrm{B\Δ}\stackrel{\→}{x}$

这里A、B为敏感矩阵。

假设星点采集个数为m，联合x和y方向的偏差和敏感矩阵，假设：

$p=\left(\begin{array}{c}\mathrm{\Δ}{x}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{x}_m\\ \mathrm{\Δ}{y}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{y}_m\end{array}\right),$ $M=\left(\begin{array}{c}{A}_1\\ .\\ .\\ .\\ A_m\\ B_1\\ .\\ .\\ .\\ B_m\end{array}\right)$

这里p是由x和y方向残留偏差组成的向量，M为A和B两个敏感矩阵组成的整体敏感矩阵。

于是有迭代方程为：

$\mathrm{\Δ}{\stackrel{\→}{x}}^{(k+1)}=\mathrm{\Δ}{\stackrel{\→}{x}}^{\left(k\right)}-{\left(M_k^T{M}_k\right)}^{-1}{M}_k^T{p}^{\left(k\right)}---\left[12\right]$

这里k为迭代序号。k取5～20，例如k取10，迭代结束后得到稳定的数据值，这时的参数即为最后的校准结果。

仿真及结果分析。

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

视场：12×12；

像素阵列：1024×1024；

像素尺寸：0.015mm×0.015mm；

焦距：73.6059mm。

假设星点质心噪声为0均值，标准差为0.05像素的高斯噪声，每个采集位置采集次数n＝100，于是噪声水平的期望值会下降一个数量级，即0.005像素。这时主要的误差来源将会是来自转台，转台0.5角秒的精度大约相当于0.01像素。假设图像中心即为主点位置，即：

x₀＝512×0.015mm，

y₀＝512×O.015mm，

单星星光模拟器发生的星光方向在转台初始坐标系N中的坐标为：

α＝45度，

β＝89度，

安装偏转角φ₁＝-1度，

安装俯仰角φ₂＝1度，

机壳和靶面安装滚转角φ₃＝2度

焦距偏差为0.2mm，

径向畸变参数为：

q1＝2e-4，q2＝-4e-7，q3＝1e-8，

p1＝2e-4，p2＝2e-4，p3＝4e-6，

进行参数计算前，首先给出参数的初始估计值，我们分别假设初始值为：

α＝0度，

β＝90度，

安装偏转角φ₁＝0度，

安装俯仰角φ₂＝0度，

机壳和靶面安装滚转角φ₃＝0度

初始焦距为73.8059mm，

径向畸变参数为：

q1＝0，q2＝0，q3＝0，

p1＝0，p2＝0，p3＝0，

在不添加噪声的情况下，根据上述的步骤，经过10次迭代计算可以得到：

	α	β	α0	β0	φ0	Δf
							1	0	0	0	0	0	0
2	0	-0.78134	-0.27722	1.0802	2.0207	-0.28309
							3	52.265	-0.70166	-1.0056	0.99466	2	-0.20342
4	43.379	-0.83839	-0.88902	0.89255	1.9977	-0.13412
							5	44.591	-1.008	-0.99813	1.0134	1.9998	-0.19795
6	44.834	-0.996 32	-0.99536	0.99943	2	-0.19999
							7	45.007	-1.0001	-1.0002	1	2	-0.2
8	44.999	-0.99998	-0.99998	1	2	-0.2
							9	45	-1	-1	1	2	-0.2
10	45	-1	-1	1	2	-0.2
								q1	q2	q3	p1	p2	p3
1	0	0	0	0	0	0
							2	0.00019752	-3.6089e-007	9.8124e-009	4.8357e-005	6.4017e-005	0
3	0.00020039	-3.9649e-007	1.0049e-008	0.0002042	0.00019477	-0.00014556
							4	0.00016925	2.1958e-007	6.5892e-009	0.00026416	0.00020973	0.0029887
5	0.0001983	-3.7644e-007	9.8208e-009	0.00020072	0.00020132	0.00052435
							6	0.00019994	-3.9905e-007	9.9935e-009	0.00020109	0.00019873	1.2493e-005
7	0.0002	-3.9994e-007	9.9996e-009	0.00019999	0.00020001	4.169e-006
							8	0.0002	-3.9999e-007	1e-008	0.0002	0.0002	4.0596e-006
9	0.0002	-4e-007	1e-008	0.0002	0.0002	4.0006e-006
							10	0.0002	-4e-007	1e-008	0.0002	0.0002	4e-006

从上表可以看出，在没有噪声时候，经过5次迭代，可以得到校准参数精确的估计值。

进而考虑实际校准过程中噪声的引入，假设经过多次采样取平均值后的质心残留噪声水平为0均值，均方差为0.005像素的高斯白噪声。仿真结果如下：

	α	β	α0	β0	φ0	Δf
							1	0	0	0	0	0	0
2	0	-0.78034	-0.27723	1.0792	2.0207	-0.28443

3	52.712	-0.70021	-1.0108	0.99319	2	-0.20414
							4	43.696	-0.83972	-0.89394	0.88975	1.9977	-0.13338
5	44.864	-1.0112	-1.0039	1.0122	1.9998	-0.19851
							6	45.092	-0.99869	-1.0002	0.99791	1.9999	-0.20065
7	45.257	-1.0025	-1.0049	0.99855	2	-0.20065
							8	45.25	-1.0023	-1.0047	0.99853	2	-0.20066
9	45.25	-1.0023	-1.0048	0.99853	2	-0.20066
							10	45.25	-1.0023	-1.0048	0.99853	2	-0.20066
	q1	q2	q3	p1	p2	p3
							1	0	0	0	0	0	0
2	0.00019872	-3.8045e-007	9.9013e-009	4.8446e-005	6.4127e-005	0
							3	0.00020102	-4.0623e-007	1.0093e-008	0.00020398	0.00019447	-8.8312e-005
4	0.00016913	2.2429e-007	6.5547e-009	0.00026478	0.00021029	0.0030273
							5	0.00019881	-3.8412e-007	9.8517e-009	0.00020039	0.00020137	0.0005469
6	0.00020054	-4.0829e-007	1.0034e-008	0.00020081	0.00019856	2.6804e-005
							7	0.00020059	-4.0912e-007	1.004e-008	0.00019976	0.00019979	1.7888e-005
8	0.0002006	-4.092e-007	1.004e-008	0.00019977	0.00019978	1.792e-005
							9	0.0002006	-4.092e-007	1.004e-008	0.00019976	0.00019978	1.7862e-005
10	0.0002006	-4.092e-007	1.004e-008	0.00019976	0.00019978	1.7862e-005

以上结果可以看出，由于噪声的影响，最小二乘法求出的参数估计值和真实值之间有所误差。由于该优化方法是全局优化的方法，参数的偏差并不会影响最终的校准精度。为了验证所获得的参数的精度，随机采集100个点，每个点的噪声水平为0.05像素。通过对比仿真的真实位置和校准后的星点位置，可以得到误差统计平方根为：x方向0.052像素，y方向为0.049像素，和质心算法噪声水平相一致，因此本发明校准精度完全满足要求。

Claims (1)

1、一种星敏感器内外方元素校准方法，基于一个由气垫平台、单星星光模拟器、星敏感器、2维轴向转台以及数据处理计算机组成的校准系统，单星星光模拟器、星敏感器、2维轴向转台安装在气垫平台上，星敏感器安装在2维轴向转台的内框上，星敏感器的视轴垂直于2维轴向转台的两个旋转轴，单星星光模拟器发生的星光方向和星敏感器的视轴对齐，其特征在于，校准的步骤如下：

1.1、校准系统建模；

1.1.1、外部参数建模；影响星敏感器采集得到的星点坐标位置的外部因素有：星光模拟器发生的星光方向，转台内外两个转动框偏离转台初始位置的角度，星敏感器在转台内框上的安装偏差，以及星敏感器图像传感器成像靶面和机壳之间的安装偏差；通过建立多个坐标系，将上述外部参数联系起来，以便分析星敏感器靶面上的星点成像位置；外部参数建模的步骤如下：

1.1.1.1、根据转台的初始位置建立坐标系N；

坐标系N以内框的转动轴为Xn轴，外框转动轴为Yn轴，Xn轴和Yn轴交于星敏感器镜头中心点0点，0为坐标系原点，过0点沿着星敏感器视轴方向为坐标系的Zn轴；星光模拟器发生的星光方向矢量Vn在坐标系N中的表达式为：

$v_n=\left[\begin{array}{c}{n}_1\\ n_2\\ n_3\end{array}\right]=\left(\begin{array}{c}\cos \mathrm{\β}\cos \mathrm{\α}\\ \cos \mathrm{\β}\sin \mathrm{\α}\\ \sin \mathrm{\β}\end{array}\right)---\left[1\right]$

式中，n₁、n₂、n₃为星光矢量Vn在坐标系N中的3个方向分量，α和β分别为该矢量在坐标系N中的偏航和俯仰角；

1.1.1.2、建立坐标系B；

当外框转过一个角度θ₁后，根据转动后的转台状态建立坐标系B，坐标系B的原点和坐标系N相同，其Yb轴和Yn轴重合，Xb轴、Zb轴分别位于Xn轴、Zn轴转动后的位置，当外框转过θ₁角后，星光矢量在坐标系B中的表达式为：

$v_b=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]=R_{\mathrm{bn}}{v}_n=\left[\begin{array}{ccc}\cos {\mathrm{\θ}}_1& 0& -\sin {\mathrm{\θ}}_1\\ 0& 1& 0\\ \sin {\mathrm{\θ}}_1& 0& \cos {\mathrm{\θ}}_1\end{array}\right]\left[\begin{array}{c}{n}_1\\ n_2\\ n_3\end{array}\right]---\left[2\right]$

式中，b₁，b₂，b₃为星光矢量在坐标系B中的3个方向分量，R_bn为转换矩阵；

1.1.1.3、建立坐标系C；

当内框转过一个角度θ₂后，根据转动后的转台状态建立坐标系C；坐标系C的原点和坐标系B相同，其Xc轴和Xb轴重合，Yc轴、Zc轴分别位于Yb轴、Zb轴转动后的位置；当外框转过θ₂角后，星光矢量在坐标系C中的表达式为：

$v_c=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]R_{\mathrm{cb}}{v}_b=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\mathrm{\θ}}_2& \sin {\mathrm{\θ}}_2\\ 0& -\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]---\left[3\right]$

式中，c₁，c₂，c₃为星光矢量在坐标系C中的3个方向分量，R_cb为转换矩阵；

1.1.1.4、建立坐标系D和E；

假设星敏感器在转台内框上的安装偏角为₁和₂，₁为偏航方向安装误差，₂为俯仰方向安装误差；对应于星敏感器这两个安装位置，分别建立坐标系D和E；

坐标系D的原点和坐标系C相同，其Xd轴和Xc轴重合，Yd轴、Zd轴分别位于Yc轴、Zc轴转动后的位置，当外框转过₁角后，星光矢量在坐标系D中的表达式为：

式中，d₁，d₂，d₃为星光矢量在坐标系D中的3个方向分量，R_dc为转换矩阵；

坐标系E的原点和坐标系D相同，其Xe轴和Xd轴重合，Ye轴、Ze轴分别位于Yd轴、Zd轴转动后的位置，当外框转过₂角后，星光矢量在坐标系E中的表达式为：

式中，e₁、e₂、e₃为星光矢量在坐标系E中的3个方向分量，R_ed为转换矩阵；坐标系E即为星敏感器机壳坐标系；

1.1.1.5、建立坐标系F；

图像传感器靶面的X，Y轴和机壳坐标系E的Xe，Ye轴之间的安装偏差造成的偏置角为₃，根据图像传感器靶面建立星敏感器的针孔成像坐标系F；坐标系F和坐标系E原点相同，Zf轴和Ze轴相重合，Xf轴和靶面图像传感器靶面的X方向一致，Yf轴和靶面Y方向一致；靶面坐标系与图像传感器采集得到的图像一致，原点位于图像左下角，X，Y方向即为图像的横、纵方向，星光矢量从坐标系E到坐标系F的转换过程为：

式中，f₁、f₂、f₃为星光矢量在坐标系F中的3个方向分量，R_fe为转换矩阵；

根据以上坐标系的转换过程，得到从转台坐标系N到星敏感器坐标系F的转换关系为：

$v_f=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=R_{\mathrm{fe}}{R}_{\mathrm{ed}}{R}_{\mathrm{dc}}{R}_{\mathrm{cb}}{R}_{\mathrm{bn}}{v}_n---\left[7\right]$

展开式[7]得到：

f₁＝cos₃ cos₂cosθ₁cosβcosα+sin₃sin(θ₂+₁)sinθ₁cosβcosα

-cos₃sin₂cos(θ₂+₁)sinθ₁cosβcosα

+sin₃cos(θ₂+₁)cosβsinα+cos₃sin₂sin(θ₂+₁)cosβsinα

-cos₃cos₂sinθ₁sinβ+sin₃sin(θ₂+₁)cosθ₁sinβ

-cos₃sin₂cos(θ₂+₁)cosθ₁sinβ

f₂＝-sin₃cos₂cosθ₁cosβcosα+cos₃sin(θ₂+₁)sinθ₁cosβcosα

+sin₃sin₂cos(θ₂+₁)sinθ₁cosβcosα

+cos₃cos(θ₂+₁)cosβsinα-sin₃sin₂sin(θ₂+₁)cosβsinα

+sin₃cos(₂)sinθ₁sinβ+cos₃sin(θ₂+₁)cosθ₁sinβ

+sin₃sin₂cos(θ₂+₁)cosθ₁sinβ

f₃＝sin₂cosθ₁cosβcosα+cos₂cos(θ₂+₁)sinθ₁cosβcosα

-cosθ₂sin(θ₂+₁)cosβsinα-sin₂sinθ₁sinβ+cos₂cos(θ₂+₁)cosθ₁sinβ

上述外部参数中，转台的偏航角和俯仰角由转台自身提供，需要求取的外部参数有5个，即为α、β、φ₁、φ₂、φ₃；

取的外部参数有5个，即为α、β、φ₁、φ₂、φ₃；

1.1.2、针孔成像建模；

根据针孔成像原理，在坐标系F中，Zf为星敏感器视轴方向，Xf，Yf和星敏感器靶面坐标系的X，Y轴方向一致，这里的靶面坐标系为2维平面坐标系，以主点0’为原点，其X，Y轴对应于采集得到图像的横纵轴，星光经过镜头中心点0后，投影成像在星敏感器靶面P(x，y)点；假设星敏感器的焦距为f_c，主点0’的图像坐标为(x₀，y₀)，根据共线性公式有：

$x=f_c\frac{f_1}{f_3}+x_0,$ $y=f_c\frac{f_2}{f_3}+y_0---\left[8\right]$

式中，主点(x₀，y₀)和焦距f_c为3个未知参数；

1.1.3、镜头畸变建模；

dx，dy为星敏感器在x方向和y方向的畸变偏差，有：

$\mathrm{dy}=\stackrel{\‾}{y}[q_1{r}^2+q_2{r}^4+q_3{r}^6]+[p_2(r^2+2{\stackrel{\‾}{y}}^2)+2p_1\stackrel{\‾}{\mathrm{xy}}][1+p_3{r}^2]$

式中，

$\left\{\begin{array}{c}\stackrel{\‾}{x}=x-x_0\\ \stackrel{\‾}{y}=y-y_0\end{array}\right.;$

$r^2={\stackrel{\‾}{x}}^2+{\stackrel{\‾}{y}}^2;$

q₁、q₂、q₃为径向畸变系数；

p₁、p₂、p₃为偏心畸变系数；

于是内部参数总共有9个，即为x₀、y₀、f_c、q₁、q₂、q₃、p₁、p₂、p₃，根据公式[8]有：

$x=f_c\frac{f_1}{f_3}+x_0+\mathrm{dx}---\left[10\right]$

$y=f_c\frac{f_2}{f_3}+y_0+\mathrm{dy}$

上述模型中，总共有外部参数和内部参数14个，即：α、β、φ₁、φ₂、φ₃、x₀、y₀、f_c、q₁、q₂、q₃、p₁、p₂、p₃；

1.2、采集数据；

转动转台，从-6度到+6度，每间隔1度为一个采集位置，在一个位置采集n次数据，n可取100～1000，并记录当时的转台转动坐标，最终使得星点成像遍布靶面，公式为：

$\tilde{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

$\tilde{y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i$

这里x_i，y_i为每次采集得到星点成像坐标，

为n次采集后的平均坐标。

1.3、数据处理；

主点位置采用自准直方法，星光垂直投射到靶面上，求出此时的图像坐标；此时，需要校准的参数有12个，用一个向量 来表示这些向量有： 根据公式进一步可以得到：

$\stackrel{\‾}{x}=x-x_0=f_c\frac{f_1}{f_3}+\mathrm{dx}=f_x\left(\stackrel{\→}{x}\right)---\left[11\right]$

$\stackrel{\‾}{y}=y-y_0=f_c\frac{f_2}{f_3}+\mathrm{dy}=f_y\left(\stackrel{\→}{x}\right)$

由于f_x和f_y均为非线性函数，因此采用非线性最小二乘迭代方法来估计参数向量

假设

为实际测量得到的

的估计值，

为向量估计偏差；则有：

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

$\mathrm{\Δy}=\stackrel{\‾}{y}-\hat{y}\≈\mathrm{B\Δ}\stackrel{\→}{x}$

这里A、B为敏感矩阵；

星点采集个数为m，联合x和y方向的偏差和敏感矩阵，

$\begin{array}{cc}p=\left(\begin{array}{c}\mathrm{\Δ}{x}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{x}_m\\ \mathrm{\Δ}{y}_1\\ .\\ .\\ .\\ \mathrm{\Δ}{y}_m\end{array}\right),& M=\left(\begin{array}{c}{A}_1\\ .\\ .\\ .\\ A_m\\ B_1\\ .\\ .\\ .\\ B_m\end{array}\right)\end{array}$

这里p是由x和y方向残留偏差组成的向量，M为A和B两个敏感矩阵组成的整体敏感矩阵；

于是有迭代方程为：

${\mathrm{\Δ}\stackrel{\→}{x}}^{(k+1)}=\mathrm{\Δ}{\stackrel{\→}{x}}^{\left(k\right)}-{\left(M_k^T{M}_k\right)}^{-1}{M}_k^T{p}^{\left(k\right)}---\left[12\right]$

这里k为迭代序号，k取5～20，迭代结束后得到稳定的数据值，这时的参数即为最后的校准结果。

Images