一种基于星场的星敏感器校准方法

专利名称：一种基于星场的星敏感器校准方法
技术领域：
本发明属于航天测量技术，涉及对星敏感器校准方法的改进。
背景技术：
星敏感器是一种利用恒星观测，为空间飞行器提供高精度姿态信息的航天测量仪器。其工作原理为星敏感器前端摄像头单元利用CCD(或CMOS)图像传感器拍摄得到星像，经过图像处理程序得到恒星像点的质心坐标和亮度的信息，然后星图识别程序利用这些信息在导航星库中找到对应的恒星，最后计算出星敏感器的三轴姿态。
星敏感器校准一般分为地面校准和在轨校准两种方式。在飞行器发射前，星敏感器要首先进行地面校准和测试，在实验室的条件下标定出星敏感器的内部参数，如主点、焦距和畸变系数等。在轨校准是在飞行器发射升空后，星敏感器根据发射时的冲击和工作环境的变化情况，如重力、大气和温度等都会不同于地面情况，需要在地面校准参数的基础上适时地对这些参数进行修正。
如图1所示，П表示星敏感器靶面，星敏感器主点为O，图像坐标为(x0，y0)，焦距为f，表示镜头中心到靶面的距离。以主点为原点建立星敏感器靶面坐标系，Z轴为主光轴，X轴和Y轴分别对应星敏感器图像采集的行方向和列方向。设星点i和j的成像点为Pi(xi，yi)，Pj(xj，yj)。星库星向量夹角成为星内角，设为θij，同时测量成像点向量夹角，称为测量星内角，设为θ′ij。根据星向量夹角的正交变换不变性，理想小孔成像模型时，θij＝θ′ij。
主点和焦距的校准。最常见的一种主点和焦距的在轨校准方法是利用星库星向量间夹角和测量星向量间夹角相等原理(如果测量参数出现偏差，则该夹角也将会出现偏差)，采用非线性参数优化的方法，求解出主点和焦距偏差，使得最终的星向量夹角误差在最小二乘意义下最小。
畸变系数的校准。在轨校准畸变系数的方法常见的有两种，一种是利用姿态信息的校准方法，另一种则是不依赖姿态信息的校准方法。利用姿态信息的校准方法是跟据主点和焦距建立理想的小孔成像模型，然后利用姿态信息计算星场星向量在靶面上的理想投影位置，进而计算它和实际测量的星点质心位置偏差。利用2维多项式模型来拟合这一偏差，完成对畸变的校准。这种方法的一个明显问题是由于对姿态信息的依赖而引入误差。由于星敏感器本身是空间飞行器上精度最高的姿态测量仪器，飞行器在星敏感器校准前只能利用陀螺或者其它姿态测量仪器获得当前姿态信息，这些外部姿态信息带有明显的误差。于是该校准方法将会把这些姿态误差引入星敏感器的参数估计过程中，引起星敏感器的最终工作误差。
不依赖姿态信息的校准方法是和校准主点和焦距的方法相似，利用星内角不变原理，计算测量星内角和星库星内角的偏差。利用2维多项式建立测量星向量模型，由于不依赖姿态信息，需要将两个测量星向量相减来表达星向量之间的相对值，然后利用2维多项式模型来拟合该相对值，求出最小二乘解。这种方法的第一个缺点是由于星向量相减而损失了0阶项，使得误差拟合不完全；第二个缺点是算法对2维多项式的建模方法有很强的依赖性，模型不好会导致数值计算不稳定，误差会很大。

发明内容
本发明的目的是针对上述星敏感器在轨校准方法存在的问题，提出一种基于星场的星敏感器校准方法。该方法对星敏感器姿态参数和内部参数同时进行估计运算，并利用非线性最小二乘法和共线性公式，迭代计算出星敏感器姿态参数和内部参数。该方法不仅可以用在星敏感器的在轨校准过程中，也可以用在夜空拍摄和实验室校准，因此该方法有着广泛是适用性。
本发明的技术方案是一种基于星场的星敏感器校准方法，其特征在于，1、建立星敏感器姿态转换矩阵；1.1、建立天球坐标系；以地球中心O′为坐标原点，过春分点的轴线为Xn轴，过北极点的轴线为Zn轴，Yn轴则定义为垂直于XnZn平面的直线；1.2、建立星敏感器坐标系O-XYZ，以主点为原点建立星敏感器靶面坐标系，Z轴为主光轴，X轴和Y轴分别对应星敏感器图像采集的行方向和列方向；星敏感器的姿态角由赤经α0、赤纬β0、和滚转角φ0组成，α0为Z轴在XnYn面上的投影与Xn轴的夹角，从Xn轴起逆时针计算；β0为Z轴与它在XnYn面上的投影之间的夹角，从投影逆时针计算；φ0为Zn轴在XY面上的投影与Yn轴的夹角，从投影起顺时针计算；φ0为天球坐标系子午面和像平面的交线和像面Y轴之间的夹角；1.3、建立姿态转换矩阵；天球坐标系经过三次旋转即可到达星敏感器坐标系，旋转过程为第一次绕Zn轴旋转 使得Xn轴和子午面相垂直；第二次绕旋转后的Xn轴旋转 使得Zn轴与Z轴重合；第三次绕两次旋转后的Zn轴旋转φ0，则天球坐标系O′-XnYnZn和星敏感器坐标系O-XYZ重合；设姿态转换矩阵为M，则有M=a1a2a3b1b2b3c1c2c3=cos(α0-π2)-sin(α0-π2)0sin(α0-π2)cos(α0-π2)0001]]> [1]这里a1＝sinα0cos0-cosα0sinβ0sin0；a2＝-sinα0sin0-cosα0sinβ0cos0；a3＝-cosα0cos0；b1＝-cosα0cos0-sinα0sinβ0sin0；b2＝cosα0sin0-sinα0sinβ0cos0；b3＝-sinα0cos0；c1＝cosβ0sin0；c2＝cosβ0cos0；c3＝-sinβ0；设星场中有n颗星，第i颗星的惯性坐标为(αi，βi)，i＝1，...，n，则第i颗星的方向矢量为vi=n1n2n3=cosβicosαicosβisinαisinβi---[2]]]>这里n1，n2，n3为星光向量在天球坐标系3个轴的投影；假设焦距为fc，星敏感器成像的共线性公式为x‾=-fca1n1+b1n2+c1n3a3n1+b3n2+c3n3---[3]]]>y‾=-fca2n1+b2n2+c2n3a3n1+b3n2+c3n3---[4]]]>这里x‾=x-x0y‾=y-y0,]]>(x0，y0)为主点坐标，(x，y)为像平面坐标；2、建立星敏感器畸变模型；假设dx，dy为x方向和y方向的畸变偏差，有 dy=y‾[q1r2+q2r4+q3r6]+[p2(r2+2y‾2)+2p1xy‾][1+p3r2]]]>[5]这里，x‾=x-x0y‾=y-y0;]]>r2＝x2+y2；q1，q2，q3为径向畸变系数；p1，p2，p3为偏心畸变系数；这里内部参数有9个，即为x0，y0，fc，q1，q2，q3，p1，p2，p3，整体成像的共线性公式为x‾=-fca1n1+b1n2+c1n3a3n1+b3n2+c3n3+dx]]>y‾=-fca2n1+b2n2+c2n3a3n1+b3n2+c3n3+dy]]>[6]3、最小二乘参数估计；主点偏差采用地面自准直校准结果；在星场校准过程中，假设主点位置为已知，这样总共需要考虑的校准参数有10个，即为α0，β0，φ0，fc，q1，q2，q3，p1，P2，p3，用参数向量来表示为
根据共线性公式有x‾=fx(x→)]]>y‾=fy(x→)]]>[8]由于fx和fy均为非线性函数，因此采用非线性最小二乘迭代方法来估计参数向量 假设 为实际测量得到的x，y的估计值， 为向量估计偏差，则有Δx=x‾-x^≈AΔx→]]>Δy=y‾-y^≈BΔx→]]>[9] 这里A、B为敏感矩阵，求出 的各项偏导数即可获得；假设星点采集个数为m，联合x和y方向的偏差和敏感矩阵，假设p=Δx1···ΔxmΔy1···Δym,M=A1···AmB1···Bm]]>这里p是由x和y方向残留偏差组成的向量，M为A和B两个敏感矩阵组成的整体敏感矩阵；于是有迭代方程为Δx→(k+1)=Δx→(k)-(MkTMk)-1MkTp(k)---[10]]]>这里k为迭代序号，k取5～20，迭代结束后得到稳定的数据值，这时的参数即为最后的校准结果。
本发明的优点是第一，该基于星场的方法同时将姿态估计过程和星敏感器参数校准过程结合起来，从而消除了姿态误差引入对星敏感器参数的影响；第二，该方法适用广泛，不仅可以应用于在轨校准，也可以应用于地面的夜空观测校准和实验室内部的星场校准。第三，该方法有着很好数值稳定性。


图1是星敏感器星场成像示意图。
图2是本发明方法中星敏感器姿态角示意图。
图3是模拟星场示意图。
具体实施例方式
下面对本发明做进一步详细说明。星场校准的方法在星敏感器的在轨校准、夜空拍摄校准和实验室利用星场模拟器进行校准的场合都是需要的，因此研究一种通用的基于星场的校准方法有着现实意义。本发明方法星敏感器姿态参数和内部参数同时进行估计运算，并利用非线性最小二乘法和共线性公式，迭代计算出星敏感器姿态参数和内部参数。它不仅可以用在星敏感器的在轨校准过程中，也可以用在夜空拍摄和实验室校准，因此本发明方法有着广泛是适用性。同时，由于该方法同时将星敏感器的姿态参数和内部参数同时进行校准，有效消除了姿态误差引入带来的星敏感器内部参数估计误差。本发明方法的步骤如下1、建立星敏感器姿态转换矩阵。
1.1、建立天球坐标系。在球面天文学中，为了与人们的直观感觉相适应，把天空假想成一个巨大的球面，即为天球。天球坐标系以地球中心O′为坐标原点，过春分点的轴线为Xn轴，过北极点的轴线为Zn轴，Yn轴则定义为垂直于XnZn平面的直线，如图2所示。
1.2、建立星敏感器坐标系。图2中O-XYZ为星敏感器坐标系。星敏感器的姿态角由赤经α0、赤纬β0、和滚转角φ0组成。这里α0为Z轴在XnYn面上的投影与Xn轴的夹角，从Xn轴起逆时针计算；β0为Z轴与它在XnYn面上的投影之间的夹角，从投影逆时针计算；φ0为Zn轴在XY面上的投影与Yn轴的夹角，从投影起顺时针计算。
1.3、建立姿态转换矩阵。φ0为天球坐标系子午面和像平面的交线和像面Y轴之间的夹角。天球坐标系经过三次旋转即可到达星敏感器坐标系。旋转过程为第一次绕Zn轴旋转 使得Xn轴和子午面相垂直；第二次绕旋转后的Xn轴旋转 使得Zn轴与Z轴重合；第三次绕两次旋转后的Zn轴旋转φ0，则天球坐标系O′-XnYnZn和星敏感器坐标系O-XYZ重合。设姿态转换矩阵为M，则有
M=a1a2a3b1b2b3c1c2c3=cos(α0-π2)-sin(α0-π2)0sin(α0-π2)cos(α0-π2)0001]]> [1]这里a1＝sinα0cos0-cosα0sinβ0sin0；a2＝-sinα0sin0-cosα0sinβ0cos0；a3＝-cosα0cos0；b1＝-cosα0cos0-sinα0sinβ0sin0；b2＝cosα0sin0-sinα0sinβ0cos0；b3＝-sinα0cos0；c1＝cosβ0sin0；c2＝cosβ0cos0；c3＝-sinβ0。
设星场中有n颗星，第i颗星的惯性坐标为(αi，βi)，i＝1，...，n，则第i颗星的方向矢量为vi=n1n2n3=cosβicosαicosβisinαisinβi---[2]]]>这里n1，n2，n3为星光向量在天球坐标系3个轴的投影。
假设焦距为fc，星敏感器成像的共线性公式为x‾=-fca1n1+b1n2+c1n3a3n1+b3n2+c3n3---[3]]]>y‾=-fca2n1+b2n2+c2n3a3n1+b3n2+c3n3---[4]]]>这里x‾=x-x0y‾=y-y0,]]>(x0，y0)为主点坐标，(x，y)为像平面坐标。
2、建立星敏感器畸变模型。
假设dx，dy为x方向和y方向的畸变偏差，有
dy=y‾[q1r2+q2r4+q3r6]+[p2(r2+2y‾2)+2p1xy‾][1+p3r2]]]>[5]这里，x‾=x-x0y‾=y-y0;]]>r2＝x2+y2；q1，q2，q3为径向畸变系数；p1，p2，p3为偏心畸变系数。
这里内部参数有9个，即为(x0，y0，fc，q1，q2，q3，p1，p2，p3)，整体成像的共线性公式为x‾=-fca1n1+b1n2+c1n3a3n1+b3n2+c3n3+dx]]>y‾=-fca2n1+b2n2+c2n3a3n1+b3n2+c3n3+dy]]>[6]3、最小二乘参数估计。
主点偏差对于星敏感器的工作精度影响不大，可以采用地面自准直校准结果。在该星场校准过程中，假设主点位置为已知，这样总共需要考虑的校准参数有10个，即为(α0，β0，φ0，fc，q1，q2，q3，p1，p2，p3)，用参数向量来表示为 根据共线性公式有x‾=fx(x→)]]>y‾=fy(x→)]]>[8]由于fx和fy均为非线性函数，因此采用非线性最小二乘迭代方法来估计参数向量 假设 为实际测量得到的x，y的估计值， 为向量估计偏差。则有Δx=x‾-x^≈AΔx→]]>Δy=y‾-y^≈BΔx→]]>[9] 这里A、B为敏感矩阵，求出 的各项偏导数即可获得。
假设星点采集个数为m，联合x和y方向的偏差和敏感矩阵，假设
p=Δx1···ΔxmΔy1···Δym,M=A1···AmB1···Bm]]>这里p是由x和y方向残留偏差组成的向量，M为A和B两个敏感矩阵组成的整体敏感矩阵。
于是有迭代方程为Δx→(k+1)=Δx→(k)-(MkTMk)-1MkTp(k)---[10]]]>这里k为迭代序号。k取5～20，例如取10，迭代结束后得到稳定的数据值，这时的参数即为最后的校准结果。
仿真及结果分析。
仿真的星敏感器基本参数为视场12度×12度；像素阵列1024×1024；像素尺寸0.015mm×0.015mm；焦距73.6059mm。
假设星点质心噪声为0均值，标准差为0.05像素的高斯噪声。在实验室内部利用星场模拟器，可以进行多次测量求平均值的方法降低质心算法噪声水平。但是如果数据来源是在轨校准或者夜空拍摄，那么数据只可能是一次性的，噪声的影响无法通过这种方法来降低。
假设图像中心即为主点位置，即x0＝512×0.015mm，y0＝512×0.015mm，星敏感器的姿态角为赤经ra＝0度，估计误差为0.5度；赤纬dec＝0度，估计误差为0.4度；滚转roll＝0度，估计误差为0.3度；这时模拟星场如图3所示设焦距偏差为0.2mm，径向畸变参数为q1＝2e-4，q2＝-4e-7，q3＝1e-8，
p1＝2e-4，p2＝2e-4，p3＝4e-6，畸变参数的估计值均为0。
为了验证算法，首先不添加噪声，根据上面提出的算法，经过10次迭代计算可以得到

从上表可以看出经过4次迭代，参数已经收敛于当初的假设值。
下面添加均值为0，均方差为0.05像素的质心高斯噪声。

同样噪声水平下，根据估计参数产生一组数据，对应于测量获得的真实数据，计算统计偏差得到x方向0.0737像素，y方向为0.0744像素。姿态估计误差角的统计值为2.4角秒，由于仿真星图中星点分布不均匀，同时星点个数不多，导致该误差角偏大，实际的姿态误差估计精度在1个多角秒。对比现有的在轨校准方法，本校准方法的精度提高了约1倍。
权利要求
1.一种基于星场的星敏感器校准方法，其特征在于，1.1、建立星敏感器姿态转换矩阵；1.1.1、建立天球坐标系；以地球中心O′为坐标原点，过春分点的轴线为Xn轴，过北极点的轴线为Zn轴，Yn轴则定义为垂直于XnZn平面的直线；1.1.2、建立星敏感器坐标系O-XYZ，以主点为原点建立星敏感器靶面坐标系，Z轴为主光轴，X轴和Y轴分别对应星敏感器图像采集的行方向和列方向；星敏感器的姿态角由赤经α0、赤纬β0、和滚转角φ0组成，α0为Z轴在XnYn面上的投影与Xn轴的夹角，从Xn轴起逆时针计算；β0为Z轴与它在XnYn面上的投影之间的夹角，从投影逆时针计算；φ0为Zn轴在XY面上的投影与Yn轴的夹角，从投影起顺时针计算；φ0为天球坐标系子午面和像平面的交线和像面Y轴之间的夹角；1.1.3、建立姿态转换矩阵；天球坐标系经过三次旋转即可到达星敏感器坐标系，旋转过程为第一次绕Zn轴旋转 使得Xn轴和子午面相垂直；第二次绕旋转后的Xn轴旋转 使得Zn轴与Z轴重合；第三次绕两次旋转后的Zn轴旋转φ0，则天球坐标系O′-XnYnZn和星敏感器坐标系O-XYZ重合；设姿态转换矩阵为M，则有M=a1a2a3b1b2b3c1c2c3=cos(α0-π2)-sin(α0-π2)0sin(α0-π2)cos(α0-π2)0001]]> [1]这里a1＝sinα0cos0-cosα0sinβ0sin0；a2＝-sinα0sin0-cosα0sinβ0cos0；a3＝-cosα0cos0；b1＝-cosα0cos0-sinα0sinβ0sin0；b2＝cosα0sin0-sinα0sinβ0cos0；b3＝-sinα0cos0；c1＝cosβ0sin0；c2＝cosβ0cos0；c3＝-sinβ0；设星场中有n颗星，第i颗星的惯性坐标为(αi，βi)，i＝1，...，n，则第i颗星的方向矢量为vi=n1n2n3=cosβicosαicosβisinαisinβi...[2]]]>这里n1，n2，n3为星光向量在天球坐标系3个轴的投影；假设焦距为fc，星敏感器成像的共线性公式为x‾=-fca1n1+b1n2+c1n3a3n1+b3n2+c3n3...[3]]]>y‾=-fca2n1+b2n2+c2n3a3n1+b3n2+c3n3...[4]]]>这里x‾=x-x0y‾=y-y0,]]>(x0，y0)为主点坐标，(x，y)为像平面坐标；1.2、建立星敏感器畸变模型；假设dx，dy为x方向和y方向的畸变偏差，有 dy=y‾[q1r2+q2r4+q3r6]+[p2(r2+2y‾2)+2p1xy‾I1+p3r2]]]>[5]这里，x‾=x-x0y‾=y-y0;]]>r2＝x2+y2；q1，q2，q3为径向畸变系数；p1，p2，p3为偏心畸变系数；这里内部参数有9个，即为x0，y0，fc，q1，q2，q3，p1，p2，p3，整体成像的共线性公式为x‾=-fca1n1+b1n2+c1n3a3n1+b3n2+c3n3+dx]]>y‾=-fca2n1+b2n2+c2n3a3n1+b3n2+c3n3+dy]]>[6]1.3、最小二乘参数估计；主点偏差采用地面自准直校准结果；在星场校准过程中，假设主点位置为已知，这样总共需要考虑的校准参数有10个，即为α0，β0，φ0，fc，q1，q2，q3，p1，p2，p3，用参数向量来表示为 根据共线性公式有x‾=fx(x→)]]>y‾=fy(x→)]]>[8]由于fx和fy均为非线性函数，因此采用非线性最小二乘迭代方法来估计参数向量 假设 为实际测量得到的x，y的估计值， 为向量估计偏差，则有Δx=x‾-x^≈AΔx→]]>Δy=y‾-y^≈BΔx→]]>[9] 这里A、B为敏感矩阵，求出 的各项偏导数即可获得；假设星点采集个数为m，联合x和y方向的偏差和敏感矩阵，假设p=Δx1···ΔxmΔy1···Δym,M=A1···AmB1···Bm]]>这里p是由x和y方向残留偏差组成的向量，M为A和B两个敏感矩阵组成的整体敏感矩阵；于是有迭代方程为Δx→(k+1)=Δx→(k)-(MkTMk)-1MkTp(k)---[10]]]>这里k为迭代序号，k取5～20，迭代结束后得到稳定的数据值，这时的参数即为最后的校准结果。
全文摘要
本发明属于航天测量技术，涉及对星敏感器校准方法的改进。本发明的步骤是建立星敏感器姿态转换矩阵；建立星敏感器畸变模型；进行最小二乘参数估计。本发明方法同时将姿态估计过程和星敏感器参数校准过程结合起来，从而消除了姿态误差引入对星敏感器参数的影响；适用广泛，不仅可以应用于在轨校准，也可以应用于地面的夜空观测校准和实验室内部的星场校准；有着很好数值稳定性。
文档编号B64G1/66GK1948085SQ20051011255
公开日2007年4月18日 申请日期2005年10月12日 优先权日2005年10月12日
发明者张广军, 郝雪涛, 江洁 申请人:北京航空航天大学