CN101750067B - 一种成像式地球敏感器地球扁率修正方法 - Google Patents

Abstract

一种成像式地球敏感器地球扁率修正方法，首先提取地球轮廓点在敏感器光电探测器上的像点坐标，计算像点在航天器本体系中的单位矢量；然后假设地球为标准球形，计算地心矢量V_EB，并根据地心矢量V_EB得到航天器的滚动角

和俯仰角θ；根据前一时刻航天器的姿态预估当前时刻航天器的偏航角ψ，同时从航天器姿态轨道控制系统获取航天器的当前位置以及坐标系转换矩阵，计算地球轮廓点的方位角σ；再根据地球轮廓点的方位角σ计算其对应的视角半径η；最后利用地球轮廓点的视角半径η消除地球扁率的影响，获得修正后的地心矢量。本发明方法操作简单，可以提供更为精确的地心矢量。

Description

一种成像式地球敏感器地球扁率修正方法

技术领域

本发明涉及一种星载地球敏感器成像时对地球扁率进行修正的方法，可以用于提高敏感器的测量精度。

背景技术

地球敏感器通过对地心矢量的测量，可以确定航天器的滚动和俯仰姿态；地心矢量和其它天体矢量(如恒星矢量)结合，还可以确定航天器的速度和位置，因此地球敏感器是未来航天器姿态轨道控制系统所不可缺少的。

地球敏感器由光学头部、传感器和信号处理等部分组成，有些还包括机械扫描部件。按照是否含有机械扫描部件，可以分为动态扫描式和静态成像式两种。动态扫描地球敏感器利用运动机械部件，带动一个或者多个视场扫过地平圆，由探测器信号确定地球的宽度和相位，计算地平圆的位置，从而获得地心矢量。静态成像式地球敏感器采用探测器，通过探测投影在其焦平面上的地球图像来计算地心矢量，和动态扫描地球敏感器相比具有质量轻、功耗低、可靠性强等优点，并可以通过适当的算法对大气模型的误差进行修正，从而提高地心方向矢量的精度。

地球敏感器的测量模型中地球是标准的球体，而地球的实际形状近似为小扁率的椭球体，地球扁率的存在使地球敏感器引入系统误差。目前，针对动态扫描式地球敏感器的地球扁率修正方法较多，但对于静态成像式地球敏感器的地球扁率修正方法研究较少，而且静态成像式地球敏感器与动态扫描式地球敏感器的修正方法有很大的不同，不能替用。

北京控制工程研究所的郝云彩在航天控制第23卷第1期上发表的《紫外月球敏感器的几个关键问题》中首次提出了采用矢量法确定月心矢量，该方法具有实现简单、计算量小、鲁棒性强等特点。但由于月球是一个标准的球形，不存在扁率问题，因此将该文中的矢量法应用于地心矢量计算时，虽然可以得到地心矢量的计算值，但是其精度受到较大的影响。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种操作简单、能够有效克服地球扁率对成像式地球敏感器的影响，从而提高地心矢量测量精度的方法。

本发明的技术解决方案是：一种成像式地球敏感器地球扁率修正方法，步骤如下：

(1)提取地球轮廓点在敏感器光电探测器上的像点坐标，计算像点在航天器本体系中的单位矢量；

(2)假设地球为标准球形，计算地心矢量V_EB，并根据地心矢量V_EB得到航天器的滚动角

和俯仰角θ，

(3)根据前一时刻航天器的姿态预估当前时刻航天器的偏航角ψ，同时从航天器姿态轨道控制系统获取航天器的当前位置以及坐标系转换矩阵，计算地球轮廓点的方位角σ，计算公式为：

E_iG＝[x y z]^T＝C_GλC_λIC_IOC_OBE_iB

其中E_iB为地球轮廓点在航天器本体系中的单位矢量，E_iG为地球轮廓点在地理坐标系中的单位矢量，C_OB为航天器本体系到轨道坐标系的转换矩阵，且 $C_{\mathrm{OB}}=C_{\mathrm{BO}}^T,$ C_IO为轨道坐标系到惯性坐标系的转换矩阵，且 $C_{\mathrm{IO}}=C_{\mathrm{OI}}^T,$ C_λI为惯性坐标系到零经度辅助坐标系的转换矩阵，且 $C_{\mathrm{\λI}}=C_{\mathrm{I\λ}}^T,$ C_Gλ为零经度辅助坐标系到地理坐标系的转换矩阵，且 $C_{\mathrm{G\λ}}=C_{\mathrm{\λG}}^T;$

i为航天器的轨道倾角、Ω为升交点赤径、ω为近地点幅角、f为真近点角，且U＝ω+f，

ζ为地心到春分点的矢量与地心到零经度方向矢量的夹角，

λ为航天器所在位置的经度，δ为地平平面法线方向与赤道平面的夹角，

$\mathrm{\σ}=\arctan \left(\frac{x}{y}\right),x\≥0,y\≥0$

(4)根据地球轮廓点的方位角σ计算其对应的视角半径η，计算公式为：

cosη＝cosη₀[1+e·F(φσ)]

e为地球的扁率，R_e为地球的赤道半径，r为航天器到地心的距离，φ为航天器所在位置的纬度；

(5)利用地球轮廓点的视角半径η消除地球扁率的影响，获得修正后的地心矢量。

所述敏感器坐标系的原点位于等效光学镜头的中心，x轴沿探测器像素的行方向，y轴沿探测器像素的列方向，z轴构成正交坐标系。

所述步骤(2)中计算地心矢量V_EB时的方法为，以X＝[E_x E_y E_z]作为观测量，H_i＝[e_ix e_iy e_iz]，Z＝[1 1 …1]，采用最小二乘法，得到X的估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^{T}Z$

其中，

所述步骤(5)中计算修正后的地心矢量V_EB时的方法为，以X＝[E_x E_y E_z]作为观测量，H_i＝[e_ix e_iy e_iz]，Z＝[cosη₁ cosη₂…cosη_n]，应用最小二乘法，得到X的估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^{T}Z$

其中V_EB＝[E_x E_y E_z]^T，E_iB＝[e_ix e_iy e_iz]^T为地球轮廓点的像点在航天器本体系中的表示，η_i为各地球轮廓点对应的视角半径。

本发明与现有技术相比的优点在于：本发明采用二次计算的方法来修正地球扁率对地心矢量测量的影响，将地心矢量的测量分为两个步骤，第一步不需要航天器姿态轨道控制系统提供任何当前时刻的姿态估计值，直接获得地心矢量的测量值，或者是航天器的滚动角和俯仰角，可以作为航天器的初始姿态敏感器；第二步根据第一步获得的航天器滚动角和俯仰角，结合航天器姿态轨道控制系统提供的当前时刻偏航角估计值，以及坐标系之间的转换矩阵，计算地球轮廓点的视角半径，从而修正地球扁率对地心矢量测量的影响，操作简单，可以提供更为精确的地心矢量。

附图说明

图1为将地球表面取为旋转椭球时的地球轮廓和地平面示意图；

图2为将地球表面取为旋转椭球时的地平轮廓圆示意图；

图3为本发明方法的流程框图；

图4为本发明实施例中采用传统矢量法时航天器滚动角和俯仰角误差随纬度幅角变化曲线；

图5为本发明实施例中采用本发明方法，航天器姿态轨道控制系统提供先验偏航信息条件下修正后航天器滚动角和俯仰角误差随纬度幅角变化曲线；

图6为本发明实施例中采用本发明方法，航天器姿态轨道控制系统未提供先验偏航信息条件下修正后航天器滚动角和俯仰角误差随纬度幅角变化曲线。

具体实施方式

本发明方法适用于静态成像式地球敏感器。

地球的形状不是理想球体，而是近似球体的不规则几何形状。本发明以旋转椭球作为成像式地球敏感器的成像模型，其它成像模型可以根据本发明提供的思路作相应的更改。

本发明的主要思路是基于旋转椭球模型的成像式地球敏感器信息处理算法，提出了一种地心矢量二次确定以消除地球扁率带来的误差，提高地心矢量确定精度的方法，第一次确定不考虑地球扁率，以球型模型计算卫星的粗姿态和粗略轨道位置信息；再根据卫星的粗略位置和粗姿态信息，结合所建立的旋转椭球模型，消除扁率影响，实现更高精度的地心矢量确定。

一、地球边缘线的数学模型

如图1所示，将地球表面取为旋转椭球，记作E，R_e为赤道半径，R_p为极半径，以地心O为原点，在零经度辅助坐标系F_λ0(原点位于地心，Z_λ0轴指向北极，X_λ0轴和Y_λ0轴在赤道平面内，X_λ0轴指向零经度点，Y_λ0轴完成正交)中，考虑地球扁率情况下E的球面方程为：

$x^2+y^2+\frac{z^2}{{(1-e)}^2}=R_e^2$

其中，e为地球扁率，且e＝(R_e-R_p)/R_e。

卫星所在位置为S，经度为λ，纬度为φ，与地心的距离为r。在卫星上看到的地球边缘称为地平轮廓线，地平轮廓线上的称为地球轮廓点，地平轮廓线所在的平面称为地平平面∏，卫星与地心的连线与椭球面E相交于S′，与地平平面相交于S^*。设S^*与地心的距离为D，由于O、S′和S^*共线，则

S^*O＝D·(n_nad)

其中，n_nad为SO的单位矢量。

由于S^*的坐标满足地平平面方程，可求得

从S^*向地平轮廓线上任意一点P引矢量ρ＝S^*P，ρ与Yλ的夹角为方位角σ。根据坐标关系，在零经度辅助坐标系F_λ0下

${\mathrm{OP}|}_{F_{\mathrm{\λ}0}}={\mathrm{OS}}^{*}+{\mathrm{\ρ}|}_{F_{\mathrm{\λ}0}}$

如图2所示，以S^*为原点，在地平平面上建立地理坐标系F_G，X_θ和Y_λ沿地平平面，X_θ沿子午线方向向北，Y_λ沿纬线方向向东，从S^*向地平轮廓线上任意一点P引矢量ρ＝S^*P，ρ与Y_λ的夹角为方位角σ。在F_G中，

F_λ0绕Z_λ0轴旋转λ得到F₁，再将F₁绕Y₁轴旋转-(90°+δ)得到F_G，F_G到F_λ0的转换矩阵为：

其中：δ为地平平面法线方向与赤道平面的夹角，且

可得OP在零经度辅助坐标系F_λ0中的坐标

由于P在椭球表面上，所以其坐标满足椭球面方程E，可得：

对三角形SPS^*应用余弦定理可得：

$\cos \mathrm{\η}=\frac{{\mathrm{SP}}^2+{\mathrm{SS}}^{*2}-S^{*}{P}^2}{2\·\mathrm{SP}\·{\mathrm{SS}}^{*}}$

其中，η为地球相对卫星的半张角或天底角，也称为从卫星上所看到的视角半径。将S、P和S^*的坐标带入公式，并忽略地球扁率的高阶影响，计算结果可简化为：

cosη＝cosη₀[1+e·F(φσ)]

式中η₀为忽略地球扁率时的视角半径值，且

$\cos {\mathrm{\η}}_0=\frac{\sqrt{r^2-R_e^2}}{r}$

F(φσ)称为地球扁率误差修正函数，并且，

具体推导过程见2002年西北工业大学钱勇的博士学位论文《高精度三轴稳定卫星姿态确定和控制系统研究》第三章。

二、球形模型下地心矢量估计算法

静态成像式地球敏感器，其原理是通过地球轮廓点在敏感器光电探测器上像点的集合，来确定地心矢量VE在航天器本体坐标系中的表达V_EB。

航天器按照3-1-2的旋转顺序分别旋转ψ、

θ，则轨道坐标系到航天器本体坐标系的转换矩阵C_BO为：

其中ψ为航天器的偏航角，

为航天器的滚动角，θ为航天器的俯仰角。轨道坐标系下，地心矢量可以表示为：

V_EO＝[0 0 1]^T

航天器本体坐标系下，地心矢量可以表示为：

V_EB＝C_BO·V_EO则：

在不考虑地球扁率的条件下，每个地平轮廓点对应的视角半径相同，每个点到透镜中心的矢量集合所围成的面是一个圆锥面，圆锥的中轴矢量正是地心矢量，这一性质不随姿态角的变化而变化。因此，地心矢量的求解可描述为已知部分圆锥母线矢量来拟合圆锥中轴的问题。

假设某地平轮廓点在探测器上的像点坐标为(x y)，根据理想光学系统成像模型(敏感器坐标系的原点位于等效光学镜头的中心，x轴沿探测器像素的行方向，y轴沿探测器像素的列方向)，该像点在敏感器坐标系下的坐标为(x y-focus)，它与敏感器坐标系原点的连线定义为边缘点矢量，其单位矢量表示为：

其中：focus为光学系统的等效焦距；

地心矢量在敏感器坐标中的三轴坐标设定为 $V_{\mathrm{EB}}=\left[\begin{array}{c}{E}_x\\ E_y\\ E_z\end{array}\right]$ ，地心矢量和边缘点矢量的点乘积为：

E_iB·V_EB＝e_ix·E_x+e_iy·E_y+e_iz·E_z圆锥的半顶角即为视角半径η，即

E_iB·V_EB＝cosη在探测器选取边缘点矢量，组成如下的方程组：

选择非单位的地心矢量表达模式：

e_ix·E_x+e_iy·E_y+e_iz·E_z＝1则：

以X＝[E_x E_y E_z]作为观测量，H_i＝[e_ix e_iy e_iz]，Z＝[1 1…1]，应用最小二乘法，得到X的估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^{T}Z$

三、椭球模型下的地心矢量补偿算法由公式

可知，球形模型下，地平轮廓点的视角半径只与地球半径和航天器的轨道高度有关，而椭球模型下，每个地平轮廓点的视角半径并不相同，其大小由地球扁率、地球半径、航天器的轨道高度、所处的位置，以及这个点在地平平面的位置决定。因此，地球廓点视角半径的测量，是消除地球扁率影响的关键。

地球的半径非常大，而扁率又非常小，在以星下点S^*为中心的地理坐标系中，视角半径η随方位角σ的增大变化非常缓慢。计算表明，卫星在赤经0°，赤纬45°，地心距离6878km(等效轨道高度500km)时，地平轮廓线上点，当方位角σ从0°增大到180°，其视角半径η的变化值Δη＝0.18°。由地心矢量计算航天器的滚动角

和俯仰角θ，以及航天器姿态轨道控制系统提供的偏航角ψ，可以获得航天器完整的姿态信息，即轨道坐标系到本体系的转换矩阵C_BO：

航天器的轨道六要素分别为：轨道倾角i、近地点幅角ω、升交点赤径Ω、轨道半长轴a、偏心率e和过近地点时刻t_p。惯性坐标系绕OZ轴旋转Ω，得到坐标系OX₁Y₁Z₁；坐标系OX₁Y₁Z₁绕OX₁轴旋转i，得到坐标系OX₂Y₂Z₂；坐标系OX₂Y₂Z₂绕OZ₂轴旋转U(U＝ω+f)，得到轨道坐标系，则惯性坐标系到轨道坐标系的转换矩阵C_OI：

假设某一轨道时刻，地心到春分点的矢量与地心到零经度方向矢量的夹角为ζ，则零经度辅助坐标系绕OZ_λ0旋转ζ得到惯性坐标系，则零经度辅助坐标系到惯性坐标系的转换矩阵C_Iλ：

而地理坐标系到零经度辅助坐标系的转换矩阵C_λG如前面所示：

由前面所述，对探测器进行图像处理，可以获得地球轮廓点在探测器平面上的像点坐标，再根据敏感器光学系统的焦距，以及敏感器在航天器本体的安装矩阵，可以获得该像点在航天器本体系中的单位矢量E_iB，E_iB经过如下矩阵变化，即为E_iB在地理坐标系中的单位矢量E_iG：

E_iG＝C_GλC_λIC_IOC_OBE_iB

其中：C_OB为航天器本体系到轨道坐标系的转换矩阵，且 $C_{\mathrm{OB}}=C_{\mathrm{BO}}^T;$

C_IO为轨道坐标系到惯性坐标系的转换矩阵，且 $C_{\mathrm{IO}}=C_{\mathrm{OI}}^T;$

C_λI为惯性坐标系到零经度辅助坐标系的转换矩阵，且 $C_{\mathrm{\λI}}=C_{\mathrm{I\λ}}^T;$

C_Gλ为零经度辅助坐标系到地理坐标系的转换矩阵，且 $C_{\mathrm{G\λ}}=C_{\mathrm{\λG}}^T.$

由单位矢量E_iG＝[x y z]^T，计算该像点对应的地球轮廓线上点的方位角σ：

$\mathrm{\σ}=\arctan \left(\frac{x}{y}\right),x\≥0,y\≥0$

由方位角σ，以及航天器提供的位置信息(距地心距离r，赤纬φ)，以及地球半径R_e、扁率e，由公式

cosη＝cosη₀[1+e·F(φσ)]

可以确定地球轮廓点对应的视角半径η。

根据以上计算获得的视角半径，重新构建最小二乘法：

E_iB·V_EB＝e_ix·E_x+e_iy·E_y+e_iz·E_z＝cosη_i

以X＝[E_x E_y E_z]作为观测量，H_i＝[e_ix e_iy e_iz]，Z＝[cosη₁ cosη₂…cosη_n]，应用最小二乘法，得到X的估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^{T}Z$

上述流程如图3所示。

实施例

为了更清楚的表明本方法的优点，建立地球的椭球仿真模型，将扁率修正后的结果与矢量法的结果进行对比。

设某卫星的轨道高度为500km，轨道倾角80°，地球半径R_e＝6378.14km，地球扁率e＝3.292×10^-3。地球敏感器的安装矩阵为单位阵，安装角ξ_x＝ξ_y＝ξ_z＝0°。

轨道周期的某一时刻，假定航天器姿态，按照理想成像方式，将地球轮廓圆上的点投影在探测器平面上，确定其理想的平面坐标，以这些平面坐标作为输入，首先以圆球模型按照矢量法计算航天器的滚动和俯仰姿态，并与真实姿态对比；再按照椭球模型，对圆球模型的姿态进行修正，并与真实姿态对比。

图4是不考虑地球扁率条件下，采用矢量法确定的航天器滚动角和俯仰角误差随纬度幅角变化曲线，可以看出，地球扁率引起的滚动角误差可以达到0.16°。

图5是采用本发明方法，在航天器姿态轨道控制系统提供先验偏航信息条件下，航天器滚动角和俯仰角误差随纬度幅角的变化曲线，可以看出，采用本方法后地球扁率引起的误差减小到10^-5°。

图6是采用本发明方法，在航天器姿态轨道控制系统不能提供先验偏航信息的条件下，航天器滚动角和俯仰角误差随纬度幅角的变化曲线，可以看出，采用本发明方法后地球扁率引起的误差减小到10^-3°以下。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (2)

1.一种成像式地球敏感器地球扁率修正方法，其特征在于步骤如下：

(1)提取地球轮廓点在敏感器光电探测器上的像点坐标，计算像点在航天器本体系中的单位矢量；

(2)假设地球为标准球形，计算地心矢量V_EB，并根据地心矢量V_EB得到航天器的滚动角

和俯仰角θ，计算地心矢量V_EB时的方法为，以X＝[E_x E_y E_z]作为观测量，H_i＝[e_ix e_iy e_iz]，Z＝[1 1 … 1]，采用最小二乘法，得到X的估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^{T}Z$

其中，

(3)根据前一时刻航天器的姿态预估当前时刻航天器的偏航角ψ，同时从航天器姿态轨道控制系统获取航天器的当前位置以及坐标系转换矩阵，计算地球轮廓点的方位角σ，计算公式为：

E_iG＝[x y z]^T＝C_GλC_λIC_IOC_OBE_iB

其中E_iB为地球轮廓点在航天器本体系中的单位矢量，E_iG为地球轮廓点在地理坐标系中的单位矢量，C_OB为航天器本体系到轨道坐标系的转换矩阵，且

C_IO为轨道坐标系到惯性坐标系的转换矩阵，且C_λI为惯性坐标系到零经度辅助坐标系的转换矩阵，且

C_Gλ为零经度辅助坐标系到地理坐标系的转换矩阵，且

$C_{\mathrm{OI}}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& -1\\ -1& 0& 0\end{array}\right]\left[\begin{array}{ccc}\cos U& \sin U& 0\\ -\sin U& \cos U& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos i& \sin i\\ 0& -\sin i& \cos i\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\Ω}& \sin \mathrm{\Ω}& 0\\ -\sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right]$

i为航天器的轨道倾角、Ω为升交点赤径、ω为近地点幅角、f为真近点角，且U＝ω+f，

$C_{\mathrm{I\λ}}=\left[\begin{array}{ccc}\cos \mathrm{\ζ}& \sin \mathrm{\ζ}& 0\\ -\sin \mathrm{\ζ}& \cos \mathrm{\ζ}& 0\\ 0& 0& 1\end{array}\right]$

ζ为地心到春分点的矢量与地心到零经度方向矢量的夹角，

$C_{\mathrm{\λG}}=\left[\begin{array}{ccc}-\cos \mathrm{\λ}\sin \mathrm{\δ}& -\sin \mathrm{\λ}& -\cos \mathrm{\λ}\cos \mathrm{\δ}\\ -\sin \mathrm{\λ}\sin \mathrm{\δ}& \cos \mathrm{\λ}& -\sin \mathrm{\λ}\cos \mathrm{\δ}\\ \cos \mathrm{\δ}& 0& -\sin \mathrm{\δ}\end{array}\right]$

λ为航天器所在位置的经度，δ为地平平面法线方向与赤道平面的夹角，

$\mathrm{\σ}=\arctan \left(\frac{x}{y}\right),x\≥0,y\≥0$

(4)根据地球轮廓点的方位角σ计算其对应的视角半径η，计算公式为：

cosη＝cosη₀[1+e·F(φ σ)]

$F\left(\begin{array}{cc}\mathrm{\φ}& \mathrm{\σ}\end{array}\right)=\frac{R_e^4}{r^2(r^2-R_e^2)}[{\sin }^2\mathrm{\φ}+\frac{r^2-R_e^2}{R_e^2}{\cos }^2\mathrm{\φ}{\sin }^2\mathrm{\σ}+2\·\frac{\sqrt{r^2-R_e^2}}{R_e}\cos \mathrm{\φ}\sin \mathrm{\φ}\sin \mathrm{\σ}];$

e为地球的扁率，R_e为地球的赤道半径，r为航天器到地心的距离，φ为航天器所在位置的纬度；

(5)利用地球轮廓点的视角半径η消除地球扁率的影响，获得修正后的地心矢量，计算方法为，以X＝[E_x E_y E_z]作为观测量，H_i＝[e_ix e_iy e_iz]，Z＝[cosη₁ cosη₂ … cosη_n]，应用最小二乘法，得到X的估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^{T}Z$

其中V_EB＝[E_x E_y E_z]^T，E_iB＝[e_ix e_iy e_iz]^T为地球轮廓点的像点在航天器本体系中的表示，η_i为各地球轮廓点对应的视角半径。

2.根据权利要求1所述的一种成像式地球敏感器地球扁率修正方法，其特征在于：所述敏感器坐标系的原点位于等效光学镜头的中心，x轴沿探测器像素的行方向，y轴沿探测器像素的列方向，z轴构成正交坐标系。

Images