CN104848860A - 一种敏捷卫星成像过程姿态机动规划方法 - Google Patents

Abstract

一种敏捷卫星成像过程姿态机动规划方法，首先根据指定的地面成像条带的起始与结束地理位置以及两个端点对应的成像时间，计算每个时间点对应的成像点位置；然后根据卫星的轨道参数以及对应时刻的成像点位置，计算卫星与成像点在地心赤道惯性坐标系中的相对位置矢量；再根据地心赤道惯性坐标系到卫星轨道坐标系的一系列转换矩阵，计算相对位置矢量在卫星本体坐标系中的分量；随后根据相对位置矢量在卫星轨道坐标系中的分量，计算卫星成像时的滚转角和俯仰角以及滚转角速度和俯仰角速度；最后根据得到的滚转角、滚转角速度、俯仰角和俯仰角速度，轨道参数以及对应的成像点位置计算卫星的偏航角与偏航角速度，得到卫星在成像过程中规划的姿态信息。

Description

一种敏捷卫星成像过程姿态机动规划方法

技术领域

本发明涉及一种适用于具有敏捷姿态机动能力的低轨道遥感卫星对地成像过程中姿态机动规划算法，特别涉及一种敏捷卫星成像过程姿态机动规划方法。

背景技术

随着遥感卫星的不断发展，遥感卫星不断丰富其自身功能并不断提高自身性能，其中敏捷机动成像能力是目前比较关注的技术。敏捷机动成像能力主要是指卫星平台具备大范围快速姿态机动能力，能够在最大姿态角和姿态机动速度的限制范围内实施灵活的对地观测，从而实现对观测目标的快速响应和多种模式的成像能力。具备敏捷机动能力的遥感卫星的成像模式除了传统卫星被动推扫成像模式外，增加了一般轨迹主动推扫成像能力。

与传统卫星相比，敏捷遥感卫星不再是卫星姿态摆到哪里就对哪里成像，而是具有根据任意位置的成像条带，卫星进行姿态机动，且能够边姿态机动边进行成像。随之而来的问题是，面对各种复杂的成像情况，如果依然采用成像任务前地面上注指令与相关信息的方法，增加了地面系统的工作量，且姿态机动成像所需上注的指令内容复杂、数据量大。面对这个问题，需要星上增加自主计算能力，实现根据成像任务的起始结束时间、起始结束位置，自行计算成像过程中卫星的姿态信息，以保证卫星相机指向正确的成像点。因此需要建立一套面向星上应用的敏捷遥感卫星成像过程中姿态机动规划的算法，实现星上自主计算姿态机动规划结果。

为了实现敏捷卫星的姿态机动规划，需要精确计算卫星的滚转角度和俯仰角度，及其这两个角度的角速度。同时为了保证卫星的成像质量，还需要精确计算卫星考虑偏流修正后的偏航角角度和角速度，滚转角度、俯仰角度和偏航角度，滚转角速度、俯仰角速度和偏航角速度分别构成了姿态信息中的姿态角度信息和姿态角速度信息。

从以上分析可知，有必要提供一种面向敏捷遥感卫星成像过程姿态机动规划算法，在满足星上使用约束的前提下，保证计算的卫星姿态准确性与精度。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提出一种敏捷卫星成像过程姿态机动规划方法，该方法针对实现星上应用，利用空间解析几何与坐标系转换，以及适当简化的地球自转模型，采用实时计算卫星成像过程中姿态信息的方法，在满足容量大小与计算速度等卫星实际工程约束的前提下，最大限度地保证卫星姿态的计算准确性与精度。

本发明的技术方案：一种敏捷卫星成像过程姿态机动规划方法，步骤如下：

(1)根据指定的地面成像条带的起始点地理位置与结束点地理位置，以及成像的起始时间与结束时间，按照星上计算步长计算获得每个时间点对应的成像点位置；所述的地理位置包括地理经度、地理纬度和海拔高度；所述的成像点位置表示为成像点在地固坐标系下的三个位置分量；

(2)获取每个时间点对应的卫星轨道参数，根据步骤(1)获得的每个时间点对应时刻的成像点位置，计算获得卫星与成像点在地心赤道惯性坐标系中的相对位置矢量；

(21)根据卫星的轨道六根数，计算卫星在地心赤道惯性坐标系中的位置；

(22)计算成像点在地心赤道惯性坐标系中的位置；

(23)根据卫星和成像点在地心赤道惯性坐标系中的位置，采用矢量计算方法计算得到地心赤道惯性系中卫星与成像点的相对位置矢量；

(3)计算步骤(2)获得的相对位置矢量在卫星本体坐标系中的分量；

(4)根据相对位置矢量在卫星轨道坐标系中的分量，计算卫星成像时的滚转角、俯仰角、滚转角的欧拉角速度、俯仰角的欧拉角速度；其中计算滚转角的欧拉角速度、俯仰角的欧拉角速度时，使用二次多项式插值的数值方法得到；

(5)根据得到的滚转角、滚转角的欧拉角速度、俯仰角、俯仰角的欧拉角速度、每个时间点对应的卫星轨道参数以及每个时间点对应的成像点位置，计算卫星的偏航角与偏航角的欧拉角速度；其中计算偏航角的欧拉角速度时，使用二次多项式插值的数值方法得到；

(6)根据得到的滚转角的欧拉角速度、俯仰角的欧拉角速度和偏航角欧拉角速度，计算得到滚转角的角速度、俯仰角的角速度和偏航角的角速度；

(7)根据得到的卫星成像时的滚转角、俯仰角、偏航角、滚转角的角速度、俯仰角的角速度、偏航角的角速度，利用星上姿态控制设备对卫星实施姿态规划。

步骤(1)中获得每个时间点对应的成像点位置的具体方法为：设地面成像条带的起始点为A，其地理经纬度为Lon_A,d和Lat_A,d，对应的起始时间为t_A；设地面成像条带的结束点为B，其地理经纬度为Lon_B,d和Lat_B,d，对应的结束时间为t_B；星上计算时间的步长为Δt；地面各个成像点位置的具体计算步骤如下：

(11)计算A点和B点分别在地固坐标系ECF中的位置坐标 $A_{\mathrm{ECF}}=\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}\end{array}\right]$ 和 $B_{\mathrm{ECF}}=\left[\begin{array}{c}{x}_{B,\mathrm{ECF}}\\ y_{B,\mathrm{ECF}}\\ z_{B,\mathrm{ECF}}\end{array}\right];$

(12)根据A点、B点的地固坐标系中的位置坐标和地心原点O的位置坐标，确定OAB平面在地固坐标系中的平面方程和该平面的法线方程；

所述OAB平面的方程具体为：

$\begin{array}{c}\left|\begin{array}{cccc}x& y& z& 1\\ x_{A,\mathrm{ECF}}& y_{A,\mathrm{ECF}}& z_{A,\mathrm{ECF}}& 1\\ x_{B,\mathrm{ECF}}& y_{B,\mathrm{ECF}}& z_{B,\mathrm{ECF}}& 1\\ 0& 0& 0& 1\end{array}\right|=(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})x\\ -(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})y\\ +(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})z\end{array}$

OAB平面的法线方向数分别为(y_A,ECFz_B,ECF-z_A,ECFy_B,ECF)、-(x_A,ECFz_B,ECF-z_A,ECFx_B,ECF)和(x_A,ECFy_B,ECF-y_A,ECFx_B,ECF)；

平面的法线方程为：

$\frac{x}{y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}}=\frac{y}{z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}}=\frac{z}{x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}$

根据法线方程，过平面法线上的一点的坐标为 $\left(\begin{array}{c}{y}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}\\ x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\end{array}\right);$

(13)在地固坐标系中求解法线与地固坐标系的z轴的夹角i；

平面OAB法线与地固坐标系的z轴两直线的夹角计算公式为：

$\cos i=\frac{(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}{\sqrt{{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2+{(z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}})}^2+{(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2}}$

，且夹角i的范围为0°～180°；

(14)求解平面OAB与地球在地固坐标系中xy平面上的交线；

设平面OAB与地固坐标系xy平面的交线矢量为地固坐标系的z轴矢量为则有：

$\begin{array}{c}\stackrel{\→}{\mathrm{ON}}=\stackrel{\→}{\mathrm{OZ}}\×\stackrel{\→}{\mathrm{oh}}=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\×\left[\begin{array}{c}{y}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}\\ x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\end{array}\right]\\ =\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ 0\end{array}\right]\end{array}$

计算得到地固坐标系x轴到矢量的旋转角度Ω：

$\cos \mathrm{\Ω}=\frac{{\mathrm{ON}}_{x,\mathrm{ECF}}}{\left|\stackrel{\→}{\mathrm{ON}}\right|}=\frac{x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}{\sqrt{{(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2+{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2}}$

如果在地固坐标系下y轴的分量y_A,ECFz_B,ECF-z_A,ECFy_B,ECF小于0，则有：

如果在地固坐标系下y轴的分量y_A,ECFz_B,ECF-z_A,ECFy_B,ECF大于等于0，则有：

$\mathrm{\Ω}=\arccos \left(\frac{x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}{\sqrt{{(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2+{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2}}\right)$

，且夹角Ω的范围为0°～360°；

其中ON_x,ECF表示为在地固坐标系下x轴的分量；

(15)建立中间坐标系S1；

建立方法为首先绕地固坐标系的z轴旋转角度Ω，然后绕此时地固坐标系的x轴旋转角度i，即可得到中间坐标系S1；两次坐标系旋转均符合右手螺旋定则；在中间坐标系S1中A点和B点的坐标为：

$\left[\begin{array}{c}{x}_{A,S1}\\ y_{A,S1}\\ z_{A,S1}\end{array}\right]=L_{S1,\mathrm{ECF}}\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}\end{array}\right],\left[\begin{array}{c}{x}_{B,S1}\\ y_{B,S1}\\ z_{B,S1}\end{array}\right]=L_{S1,\mathrm{ECF}}\left[\begin{array}{c}{x}_{B,\mathrm{ECF}}\\ y_{B,\mathrm{ECF}}\\ z_{B,\mathrm{ECF}}\end{array}\right]$

其中L_S1,ECF为：

$\begin{array}{c}{L}_{S1,\mathrm{ECF}}=L_x\left(i\right)L_z\left(\mathrm{\Ω}\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]\\ =\left[\begin{array}{ccc}\cos \mathrm{\Ω}& \sin \mathrm{\Ω}& 0\\ -\cos i\sin \mathrm{\Ω}& \cos i\cos \mathrm{\Ω}& \sin i\\ \sin i\sin \mathrm{\Ω}& -\sin i\cos \mathrm{\Ω}& \cos i\end{array}\right]\end{array};$

(16)建立坐标系S2；

计算坐标系S1下x轴到A点转动角度ω；具体方法为：

设A点在中间坐标系S1的坐标系为(x_A,S1,y_A,S1,z_A,S1)′，有：

$\cos \mathrm{\ω}=\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}$

如果z_A,ECF<0，则 $\mathrm{\&omega;}=360-\arccos \left(\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}\right),$ 如果 $z_{A,\mathrm{ECF}}\&GreaterEqual;0\mathrm{\&omega;}=\arccos \left(\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}\right);$ 且夹角ω的范围为0°～360°；

绕坐标系S1的z轴转动角度ω后，得到坐标系S2；坐标系S2中A点和B点的坐标计算公式为：

$\left[\begin{array}{c}{x}_{A,S2}\\ y_{A,S2}\\ z_{A,S2}\end{array}\right]=L_{S2,S1}\left[\begin{array}{c}{x}_{A,S1}\\ y_{A,S1}\\ z_{A,S1}\end{array}\right],\left[\begin{array}{c}{x}_{B,S2}\\ y_{B,S2}\\ z_{B,S2}\end{array}\right]=L_{S2,S1}\left[\begin{array}{c}{x}_{B,S1}\\ y_{B,S1}\\ z_{B,S1}\end{array}\right]$

其中L_S2,S1为：

$L_{S2,S1}=L_z\left(\mathrm{\&omega;}\right)=\left[\begin{array}{ccc}\cos \mathrm{\&omega;}& \sin \mathrm{\&omega;}& 0\\ -\sin \mathrm{\&omega;}& \cos \mathrm{\&omega;}& 0\\ 0& 0& 1\end{array}\right]$

(17)计算成像观测条带中A、B两点之间，卫星每隔Δt时间观测的地面点；

根据在S2坐标系中A、B两点的坐标，计算A、B两点之间的夹角α：

$\cos \mathrm{\&alpha;}=\frac{x_{A,S2}{x}_{B,S2}+y_{A,S2}{y}_{B,S2}+z_{A,S2}{z}_{B,S2}}{\sqrt{x_{A,S2}^2+y_{A,S2}^2+z_{A,S2}^2}\sqrt{x_{B,S2}^2+y_{B,S2}^2+z_{B,S2}^2}}$

得到A、B两点之间的夹角α后，根据卫星星上计算的时间间隔Δt，和成像起始时间t_A和结束时间t_B，确定观测的地面点个数：

$\mathrm{sum}=\frac{t_B-t_A}{\mathrm{\&Delta;t}}$

则卫星成像过程中指向的每个成像点在球体模型上大圆的角度间隔为

$\mathrm{\&Delta;\&alpha;}=\frac{\mathrm{\&alpha;}}{\mathrm{sum}};$

设A点的地心距为则A点在坐标系S2中的单位向量为 $\left[\begin{array}{ccc}{\tilde{x}}_{A,S2}& {\tilde{y}}_{A,S2}& {\tilde{z}}_{A,S2}\end{array}\right]={\left[\begin{array}{ccc}\frac{x_{A,S2}}{r_A}& \frac{y_{A,S2}}{r_A}& \frac{z_{A,S2}}{r_A}\end{array}\right]}^{\&prime;},$ 继而有A点在坐标系S2中的单位向量的模为 ${\tilde{r}}_A=\sqrt{{\tilde{x}}_{A,S2}^2+{\tilde{y}}_{A,S2}^2+{\tilde{z}}_{A,S2}^2}=1;$

设成像观测条带中的成像点为A，P1，P2，……，Pn，……，B，在坐标系S2中的每个成像点位置的单位向量有：

$\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},S2}\\ {\tilde{y}}_{\mathrm{Pn},S2}\\ {\tilde{z}}_{\mathrm{Pn},S2}\end{array}\right]=\left[\begin{array}{c}{\tilde{r}}_A\cos (n\&times;\mathrm{\&beta;})\\ {\tilde{r}}_A\sin (n\&times;\mathrm{\&beta;})\\ 0\end{array}\right],(n=\mathrm{0,1},...,\mathrm{sum}),\left\{\begin{array}{cc}\mathrm{\&beta;}=\mathrm{\&Delta;\&alpha;},& y_{B,S2}>0\\ \mathrm{\&beta;}=-\mathrm{\&Delta;\&alpha;},& y_{B,S2}<0\end{array}\right.$

得到Pn在坐标系S2中单位向量的位置坐标后，通过坐标系转换，计算其在地固坐标系中的位置坐标：

$\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}\\ {\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}\\ {\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}\end{array}\right]=L_z(-\mathrm{\&Omega;})L_x(-i)L_z(-\mathrm{\&omega;})\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},S2}\\ {\tilde{y}}_{\mathrm{Pn},S2}\\ {\tilde{z}}_{\mathrm{Pn},S2}\end{array}\right]$

其中 $L_z(-\mathrm{\&Omega;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&Omega;})& \sin (-\mathrm{\&Omega;})& 0\\ -\sin (-\mathrm{\&Omega;})& \cos (-\mathrm{\&Omega;})& 0\\ 0& 0& 1\end{array}\right],L_x(-i)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos (-i)& \sin (-i)\\ 0& -\sin (-i)& \cos (-i)\end{array}\right],$

$L_z(-\mathrm{\&omega;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&omega;})& \sin (-\mathrm{\&omega;})& 0\\ -\sin (-\mathrm{\&omega;})& \cos (-\mathrm{\&omega;})& 0\\ 0& 0& 1\end{array}\right];$

得到在地固坐标系中的分量后，利用空间解析几何，求解直线与空间椭球的交点：

设过原点和的空间直线为设地球椭球模型为其中a＝Re＝6378140m为地球赤道半径，c＝Rp＝6356755m为地球极半径；

根据空间解析几何，得到空间直线的两点式方程为：

$\frac{x-0}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}-0}=\frac{y-0}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}-0}=\frac{z-0}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}-0}\&DoubleRightArrow;\frac{x}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{y}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{z}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}}$

将的空间直线方程与地球椭球方程联立得到： $\left\{\begin{array}{c}\frac{x}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{y}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{z}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}}\\ \frac{x^2+y^2}{a^2}+\frac{z^2}{c^2}=1\end{array}\right.,$ 求解得到的两组解[x₁ y₁ z₁]′和[x₂ y₂ z₂]′，分别求两个解与成像起始点A之间的夹角 $\cos {\mathrm{\&gamma;}}_1=\frac{x_{A,\mathrm{ECF}}{x}_1+y_{A,\mathrm{ECF}}{y}_1+z_{A,\mathrm{ECF}}{z}_1}{\sqrt{x_{A,\mathrm{ECF}}^2+y_{A,\mathrm{ECF}}^2+z_{A,\mathrm{ECF}}^2}\sqrt{x_1^2+y_1^2+z_1^2}}$ 和 $\cos {\mathrm{\&gamma;}}_2=\frac{x_{A,\mathrm{ECF}}{x}_2+y_{A,\mathrm{ECF}}{y}_2+z_{A,\mathrm{ECF}}{z}_2}{\sqrt{x_{A,\mathrm{ECF}}^2+y_{A,\mathrm{ECF}}^2+z_{A,\mathrm{ECF}}^2}\sqrt{x_2^2+y_2^2+z_2^2}},$ 与成像起始点A夹角小的即为成像点P_n在地固坐标系中的位置。

步骤22)中计算时利用地固坐标系到地心赤道惯性坐标系的转换矩阵，且转换矩阵中只考虑岁差修正以及章动的前五项修正。

步骤(5)中计算偏航角的具体方法为：

(51)定义卫星本体系Z轴指向地心，X轴在轨道平面内垂直z轴指向卫星飞行方向，Y轴满足右手定则；设卫星本体的Z轴在卫星本体坐标系中的矢量为 ${\stackrel{\&RightArrow;}{z}}_b=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],$ 根据已经得到的滚转角与俯仰角θ，计算卫星本体Z轴在地固坐标系中的矢量表示形式 ${\stackrel{\&RightArrow;}{z}}_{\mathrm{ECF}}=\left[\begin{array}{c}{z}_{x,\mathrm{ECF}}\\ z_{y,\mathrm{ECF}}\\ z_{z,\mathrm{ECF}}\end{array}\right],$ 的计算公式为：

${\stackrel{\&RightArrow;}{z}}_{\mathrm{ECF}}=L_{\mathrm{ECF},\mathrm{ECI}}\&CenterDot;L_{\mathrm{io}}\&CenterDot;L_{\mathrm{ob}}\&CenterDot;{\stackrel{\&RightArrow;}{z}}_b$

其中L_io＝L_oi ^T，上标T表示对矩阵进行转置；L_ECF,ECI为地心赤道坐标系到地固坐标系的转换矩阵；L_ob为卫星本体坐标系到卫星轨道坐标系的转换矩阵；

(52)根据成像点在地固坐标系中的位置和卫星在地固坐标系中的位置计算卫星到成像点的距离矢量h，其数值为其中地心角为r是卫星地心距，R是成像点的地心距；

(53)计算在地心赤道惯性坐标系中，成像点相对于卫星的相对速度矢量，即地速矢量具体的计算公式为：

$\stackrel{\&RightArrow;}{v}={\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_e\&times;\stackrel{\&RightArrow;}{R}-[{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n\&times;\stackrel{\&RightArrow;}{r}+({\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n+{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_s)\&times;\stackrel{\&RightArrow;}{h}+{\stackrel{\&RightArrow;}{v}}_r]={\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_e\&times;\stackrel{\&RightArrow;}{R}-{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n\&times;\stackrel{\&RightArrow;}{R}-{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_s\&times;\stackrel{\&RightArrow;}{h}-{\stackrel{\&RightArrow;}{v}}_r$

其中为地球自转角速度矢量，是地心到成像点的矢量，是地心指向卫星的矢量，是轨道角速度矢量，其数值为 是卫星绝对速度的径向分量，其大小为 $v_r=\sqrt{\frac{\mathrm{\&mu;}}{p}}e\sin f,$ 为卫星的姿态角速度；

(54)根据上一步得到的地速以及地心赤道惯性坐标系到卫星轨道坐标系的转换矩阵L_oi，转换得到该地速在卫星轨道坐标系中的矢量其具体计算公式为：

${\stackrel{\&RightArrow;}{v}}_o=L_{\mathrm{oi}}\&CenterDot;\stackrel{\&RightArrow;}{v}$

(55)根据卫星滚转角、俯仰角和卫星轨道坐标系到卫星本体坐标系的转换矩阵L_bo，计算地速在卫星本体坐标系中的矢量具体的计算公式为：

${\stackrel{\&RightArrow;}{v}}_b=L_{\mathrm{bo}}\&CenterDot;{\stackrel{\&RightArrow;}{v}}_o$

(56)根据地速在卫星本体坐标系中的矢量 ${\stackrel{\&RightArrow;}{v}}_b=\left[\begin{array}{c}{v}_{\mathrm{xb}}\\ v_{\mathrm{yb}}\\ v_{\mathrm{zb}}\end{array}\right],$ 计算获得偏航角：

$\mathrm{\&psi;}=\mathrm{arctg}\frac{v_{\mathrm{yb}}}{v_{\mathrm{xb}}}.$

本发明的有益效果：

(1)本发明方法采用空间解析几何和坐标转换的方法计算成像条带中每个成像点的位置，并对模型进行了适当简化，在满足精度要求的前提下减少计算量。

(2)本发明采用了实时计算成像点在地心赤道坐标系中的位置，计算过程中考虑岁差、章动前5项对地球自转的影响，可以实时保证数传天线指向角度计算的准确性，并能满足星上计算机对计算速度的要求。

(3)本发明采用的算法中考虑了偏流修正的偏航角算法，能够适用于任意卫星姿态的偏流角计算。

附图说明

图1为本发明姿态机动规划算法流程框图；

图2为地面成像观测条带的示意图；

图3为地面成像点位置计算流程框图；

图4为卫星与成像点相对位置的示意图；

图5为滚转角、滚转角速度，俯仰角、俯仰角速度4个参数的计算流程框图；

图6为卫星与成像点的几何关系示意图；

图7为偏航角与偏航角速度计算流程框图。

具体实施方式

如图1所示，为本发明姿态机动规划算法的流程框图，具体步骤如下：

(1)根据指定的地面成像条带的起始与结束点地理位置以及成像的起始时间与结束时间，根据星上的计算步长，计算每个时间点对应的成像点位置。

设成像观测条带的起始点为A，其地理经纬度为Lon_A,d和Lat_A,d，对应时间为t_A；设成像观测条带的结束点为B，其地理经纬度为Lon_B,d和Lat_B,d，对应时间为t_B。星上计算时间的步长为Δt。

地面成像观测条带示意如图2所示。在本方法中，对起始点与结束点之间形成的观测条带按照切地球大圆的弧段进行计算，这样进行适当简化，能够在满足精度要求的同时，尽可能减少计算量。

地面各个成像点位置的具体计算步骤如下：

a.计算A点和B点分别在地固坐标系ECF中的位置坐标 $A_{\mathrm{ECF}}=\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}\end{array}\right]$ 和 $B_{\mathrm{ECF}}=\left[\begin{array}{c}{x}_{B,\mathrm{ECF}}\\ y_{B,\mathrm{ECF}}\\ z_{B,\mathrm{ECF}}\end{array}\right],$ 地理经纬度到地固坐标系位置坐标的转换方法可以参见肖业伦著的《航天器飞行动力学原理》的第五章。

b.由于卫星对观测条带成像时，在地面上按照大圆进行推扫，因此可以根据A点、B点的地固坐标系中的位置坐标和地心原点O的位置坐标，确定OAB平面在地固坐标系中的平面方程和该平面的法线方程。

其中OAB平面的平面方程可以利用空间解析几何的三点式平面方程得到，OAB平面的方程具体为：

$\begin{array}{c}\left|\begin{array}{cccc}x& y& z& 1\\ x_{A,\mathrm{ECF}}& y_{A,\mathrm{ECF}}& z_{A,\mathrm{ECF}}& 1\\ x_{B,\mathrm{ECF}}& y_{B,\mathrm{ECF}}& z_{B,\mathrm{ECF}}& 1\\ 0& 0& 0& 1\end{array}\right|=(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})x\\ -(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})y\\ +(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})z\end{array}$

根据空间解析几何中平面方程的一般形式，确定OAB平面的法线方向数分别为(y_A,ECFz_B,ECF-z_A,ECFy_B,ECF)、_-(x_A,ECFz_B,ECF-z_A,ECFx_B,ECF)和(x_A,ECFy_B,ECF-y_A,ECFx_B,ECF)。得到平面法线的方向数后，可以继续计算得到OAB平面的法线方程，该平面的法线过原点O，且具有方向数(y_A,ECFz_B,ECF-z_A,ECFy_B,ECF)、-(x_A,ECFz_B,ECF-z_A,ECFx_B,ECF)和(x_A,ECFy_B,ECF-y_A,ECFx_B,ECF)，则该平面法线的直线方程为：

$\frac{x}{y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}}=\frac{y}{z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}}=\frac{z}{x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}$

根据法线方程，过平面法线上的一点的坐标为 $\left(\begin{array}{c}{y}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}\\ x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\end{array}\right).$

c.在地固坐标系中求解法线与地固坐标系的z轴的夹角i。

具体解算方法如下：

在步骤b中已经得到了平面OAB的法线方程，而地固坐标系的z轴方程的方向数为(0,0,1)。则平面OAB法线与z轴两直线的夹角计算公式为：

$\cos i=\frac{(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}{\sqrt{{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2+{(z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}})}^2+{(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2}}$

夹角i的范围为0°～180°。

d.求解平面OAB与地球在地固坐标系中xy平面上的交线ON。

设平面与地固坐标系xy平面的交线ON的矢量为地固坐标系的z轴矢量为则有：

$\begin{array}{c}\stackrel{\&RightArrow;}{\mathrm{ON}}=\stackrel{\&RightArrow;}{\mathrm{OZ}}\&times;\stackrel{\&RightArrow;}{\mathrm{oh}}=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\&times;\left[\begin{array}{c}{y}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}\\ x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\end{array}\right]\\ =\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ 0\end{array}\right]\end{array}$

得到矢量后，可以计算得到地固坐标系x轴到矢量的旋转角度Ω，Ω可利用和在地固坐标系下x轴的分量ON_x,ECF获得：

$\cos \mathrm{\&Omega;}=\frac{{\mathrm{ON}}_{x,\mathrm{ECF}}}{\left|\stackrel{\&RightArrow;}{\mathrm{ON}}\right|}=\frac{x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}{\sqrt{{(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2+{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2}}$

如果在地固坐标系下y轴的分量y_A,ECFz_B,ECF-z_A,ECFy_B,ECF小于0，则有：

如果在地固坐标系下y轴的分量y_A,ECFz_B,ECF-z_A,ECFy_B,ECF大于等于0，则有：

$\mathrm{\&Omega;}=\arccos \left(\frac{x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}{\sqrt{{(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2+{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2}}\right)$

夹角Ω的范围为0°～360°。

e.建立中间坐标系S1。

建立方法为首先绕地固坐标系的z轴旋转角度Ω，然后绕此时的x轴旋转角度i，即可得到中间坐标系S1。两次坐标系旋转均符合右手螺旋定则。在新建立的中间坐标系中，S1的x轴与矢量重合，方向相同，z轴与OAB平面的法线重合且方向相同，y轴满足右手螺旋定则。在中间坐标系S1中A点和B点的坐标为：

$\left[\begin{array}{c}{x}_{A,S1}\\ y_{A,S1}\\ z_{A,S1}\end{array}\right]=L_{S1,\mathrm{ECF}}\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}\end{array}\right],\left[\begin{array}{c}{x}_{B,S1}\\ y_{B,S1}\\ z_{B,S1}\end{array}\right]=L_{S1,\mathrm{ECF}}\left[\begin{array}{c}{x}_{B,\mathrm{ECF}}\\ y_{B,\mathrm{ECF}}\\ z_{B,\mathrm{ECF}}\end{array}\right]$

其中L_S1,ECF为：

$\begin{array}{c}{L}_{S1,\mathrm{ECF}}=L_x\left(i\right)L_z\left(\mathrm{\&Omega;}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos i& \sin i\\ 0& -\sin i& \cos i\end{array}\right]\&CenterDot;\left[\begin{array}{ccc}\cos \mathrm{\&Omega;}& \sin \mathrm{\&Omega;}& 0\\ -\sin \mathrm{\&Omega;}& \cos \mathrm{\&Omega;}& 0\\ 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{ccc}\cos \mathrm{\&Omega;}& \sin \mathrm{\&Omega;}& 0\\ -\cos i\sin \mathrm{\&Omega;}& \cos i\cos \mathrm{\&Omega;}& \sin i\\ \sin i\sin \mathrm{\&Omega;}& -\sin i\cos \mathrm{\&Omega;}& \cos i\end{array}\right]\end{array}$

f.建立坐标系S2。

在中间坐标系S1中，计算x轴到A点转动角度ω。具体方法为：

设A点在中间坐标系S1的坐标系为(x_A,S1,y_A,S1,z_A,S1)′，有：

$\cos \mathrm{\&omega;}=\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}$

如果z_A,ECF<0，则 $\mathrm{\&omega;}=360-\arccos \left(\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}\right),$ 如果z_A,ECF≥0，则有 $\mathrm{\&omega;}=\arccos \left(\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}\right),$ 夹角ω的范围为0°～360°。

得到ω后，在中间坐标系S1的基础上，绕z轴转动角度ω后，得到坐标系S2。在坐标系S2中，A点位于x轴上，B点位于xy平面上。S2中A点和B点的坐标计算公式为：

$\left[\begin{array}{c}{x}_{A,S2}\\ y_{A,S2}\\ z_{A,S2}\end{array}\right]=L_{S2,S1}\left[\begin{array}{c}{x}_{A,S1}\\ y_{A,S1}\\ z_{A,S1}\end{array}\right],\left[\begin{array}{c}{x}_{B,S2}\\ y_{B,S2}\\ z_{B,S2}\end{array}\right]=L_{S2,S1}\left[\begin{array}{c}{x}_{B,S1}\\ y_{B,S1}\\ z_{B,S1}\end{array}\right]$

其中L_S2,S1为：

$L_{S2,S1}=L_z\left(\mathrm{\&omega;}\right)=\left[\begin{array}{ccc}\cos \mathrm{\&omega;}& \sin \mathrm{\&omega;}& 0\\ -\sin \mathrm{\&omega;}& \cos \mathrm{\&omega;}& 0\\ 0& 0& 1\end{array}\right]$

g.计算成像观测条带中AB两点之间，卫星每隔Δt时间观测的地面点。

根据在S2坐标系中AB两点的坐标，计算AB两点之间的夹角α。α的计算公式为：

$\cos \mathrm{\&alpha;}=\frac{x_{A,S2}{x}_{B,S2}+y_{A,S2}{y}_{B,S2}+z_{A,S2}{z}_{B,S2}}{\sqrt{x_{A,S2}^2+y_{A,S2}^2+z_{A,S2}^2}\sqrt{x_{B,S2}^2+y_{B,S2}^2+z_{B,S2}^2}}$

得到A、B两点之间的夹角α后，根据卫星星上计算的时间间隔Δt，和成像起始时间t_A和结束时间t_B，确定观测的地面点个数sum：

$\mathrm{sum}=\frac{t_B-t_A}{\mathrm{\&Delta;t}}$

则卫星成像过程中指向的地面成像点共sum+1个，每个成像点在球体模型上大圆的角度间隔Δα为

设A点的地心距为则A点在坐标系S2中的单位向量为 $\left[\begin{array}{ccc}{\tilde{x}}_{A,S2}& {\tilde{y}}_{A,S2}& {\tilde{z}}_{A,S2}\end{array}\right]={\left[\begin{array}{ccc}\frac{x_{A,S2}}{r_A}& \frac{y_{A,S2}}{r_A}& \frac{z_{A,S2}}{r_A}\end{array}\right]}^{\&prime;},$ 继而有A点在坐标系S2中的单位向量的模为 ${\tilde{r}}_A=\sqrt{{\tilde{x}}_{A,S2}^2+{\tilde{y}}_{A,S2}^2+{\tilde{z}}_{A,S2}^2}=1$

设成像观测条带中的成像点为A，P1，P2，……，Pn，……，B，在坐标系S2中的每个成像点位置的单位向量有：

$\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},S2}\\ {\tilde{y}}_{\mathrm{Pn},S2}\\ {\tilde{z}}_{\mathrm{Pn},S2}\end{array}\right]=\left[\begin{array}{c}{\tilde{r}}_A\cos (n\&times;\mathrm{\&beta;})\\ {\tilde{r}}_A\sin (n\&times;\mathrm{\&beta;})\\ 0\end{array}\right],(n=\mathrm{0,1},...,\mathrm{sum}),\left\{\begin{array}{cc}\mathrm{\&beta;}=\mathrm{\&Delta;\&alpha;},& y_{B,S2}>0\\ \mathrm{\&beta;}=-\mathrm{\&Delta;\&alpha;},& y_{B,S2}<0\end{array}\right.$

得到Pn在坐标系S2中单位向量的位置坐标后，通过坐标系转换，计算其在地固坐标系中的位置坐标：

$\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}\\ {\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}\\ {\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}\end{array}\right]=L_z(-\mathrm{\&Omega;})L_x(-i)L_z(-\mathrm{\&omega;})\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},S2}\\ {\tilde{y}}_{\mathrm{Pn},S2}\\ {\tilde{z}}_{\mathrm{Pn},S2}\end{array}\right]$

其中 $L_z(-\mathrm{\&Omega;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&Omega;})& \sin (-\mathrm{\&Omega;})& 0\\ -\sin (-\mathrm{\&Omega;})& \cos (-\mathrm{\&Omega;})& 0\\ 0& 0& 1\end{array}\right],L_x(-i)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos (-i)& \sin (-i)\\ 0& -\sin (-i)& \cos (-i)\end{array}\right],$

$L_z(-\mathrm{\&omega;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&omega;})& \sin (-\mathrm{\&omega;})& 0\\ -\sin (-\mathrm{\&omega;})& \cos (-\mathrm{\&omega;})& 0\\ 0& 0& 1\end{array}\right].$

得到在地固坐标系中的分量后，利用空间解析几何，求解直线与空间椭球的交点：

设过原点和的空间直线为设地球椭球模型为其中a＝Re＝6378140m为地球赤道半径，c＝Rp＝6356755m为地球极半径。

根据空间解析几何，可以得到空间直线的两点式方程为：

$\frac{x-0}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}-0}=\frac{y-0}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}-0}=\frac{z-0}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}-0}\&DoubleRightArrow;\frac{x}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{y}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{z}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}}$

将的空间直线方程与地球椭球方程联立得到： $\left\{\begin{array}{c}\frac{x}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{y}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{z}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}}\\ \frac{x^2+y^2}{a^2}+\frac{z^2}{c^2}=1\end{array}\right.,$ 求解得到的两组解[x₁ y₁ z₁]′和[x₂ y₂ z₂]′，分别求两个解与成像起始点A之间的夹角 $\cos {\mathrm{\&gamma;}}_1=\frac{x_{A,\mathrm{ECF}}{x}_1+y_{A,\mathrm{ECF}}{y}_1+z_{A,\mathrm{ECF}}{z}_1}{\sqrt{x_{A,\mathrm{ECF}}^2+y_{A,\mathrm{ECF}}^2+z_{A,\mathrm{ECF}}^2}\sqrt{x_1^2+y_1^2+z_1^2}}$ 和 $\cos {\mathrm{\&gamma;}}_2=\frac{x_{A,\mathrm{ECF}}{x}_2+y_{A,\mathrm{ECF}}{y}_2+z_{A,\mathrm{ECF}}{z}_2}{\sqrt{x_{A,\mathrm{ECF}}^2+y_{A,\mathrm{ECF}}^2+z_{A,\mathrm{ECF}}^2}\sqrt{x_2^2+y_2^2+z_2^2}},$ 与成像起始点A夹角小的即为成像点P_n在地固坐标系中的位置。

地面成像点位置计算流程如图3所示。

(2)获取每个时间点对应的卫星轨道参数，根据步骤(1)获得的每个时间点对应时刻的成像点位置，计算获得卫星与成像点在地心赤道惯性坐标系中的相对位置矢量。

a.根据卫星的轨道六根数(a、e、i、Ω、ω、M)，计算卫星在地心赤道惯性坐标系中的位置具体计算公式可参见肖业伦著的《航天器飞行动力学原理》的第二章第6节。

b.计算成像点在地心赤道惯性坐标系中的位置计算时需要利用地固坐标系到地心赤道惯性坐标系的转换矩阵。

${\stackrel{\&RightArrow;}{r}}_{P,\mathrm{ECI}}=L_{\mathrm{ECI},\mathrm{ECF}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_{P,\mathrm{ECF}}$

其中为地面成像点在地固坐标系中的位置矢量。

从地固坐标系到地心赤道坐标系要分别经过岁差修正、章动修正、地球自转以及极移修正四次坐标转换，转换过程中涉及的模型较为复杂，但均已非常成熟。由于算法将用于星上使用，对运算速度和存储大小有要求，故在满足计算精度的前提下，进行适当的简化，只考虑了地球岁差修正、地球章动前5项影响和地球自转三次坐标转换。

在地固坐标系到地心赤道惯性坐标系的转换中需要由标准历元到计算历元的赤道坐标系之间转换的三个赤道岁差参数ζ_A,z_A,θ_A以及格林尼治恒星时S_G四个参数。三个赤道岁差参数计算公式为：

$\left\{\begin{array}{c}{\mathrm{\&zeta;}}_A={0306}^{\&prime;\&prime;}.2181t+0^{\&prime;\&prime;}.30188t^2+0^{\&prime;\&prime;}.017998t^3\\ z_A={2306}^{\&prime;\&prime;}.2181t+1^{\&prime;\&prime;}.09468t^2+0^{\&prime;\&prime;}.018203t^3\\ {\mathrm{\&theta;}}_A={2004}^{\&prime;\&prime;}.3109t-0^{\&prime;\&prime;}.42665t^2-0^{\&prime;\&prime;}.041833t^3\end{array}\right.$

岁差转换矩阵为(PR)＝R_z(-z_A)R_y(θ_A)R_z(-ζ_A)。

在从地固坐标系到地心赤道坐标系中除了岁差影响以外，还需考虑章动的影响。

从地心赤道坐标系到地固坐标系的章动矩阵为：

(NR)＝R_x(-Δε)R_y(Δθ)R_z(-Δμ)

其中Δμ＝Δφcosε，Δθ＝Δφsinε是赤经和赤纬章动，Δε和Δφ为交角章动和黄经章动，ε为考虑岁差影响的黄赤交角。交角章动和黄经章动以及黄赤交角的计算方法如下：

考虑岁差影响的黄赤交角计算公式为：

ε＝23°26'21”.448-46”.8150t-0”.00059t²+0”.001813t³

根据IAU1980章动序列，黄经章动Δφ和交角章动Δε的计算公式如下：

$\left\{\begin{array}{c}\mathrm{\&Delta;\&phi;}=\underset{j=1}{\overset{106}{\mathrm{\&Sigma;}}}(A_{0j}+A_{1j}t)\sin \left(\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ji}}{\mathrm{\&alpha;}}_i\left(t\right)\right)\\ \mathrm{\&Delta;\&epsiv;}=\underset{j=1}{\overset{106}{\mathrm{\&Sigma;}}}(B_{0j}+B_{1j}t)\cos \left(\underset{i=1}{\overset{5}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ji}}{\mathrm{\&alpha;}}_i\left(t\right)\right)\end{array}\right.$

其中A_0j、A_1j、B_0j、B_1j和k_ji为相应的系数，5个基本变量α_i(α₁～α₅)分别是月球的平近点角l，太阳的平近点角l'，月球平升角距F，日、月平角距D和月球轨道升交点平黄经Ω，有：

其中1^r＝360°＝1296000”，t是从标准历元J2000.0起算的世纪数。下表列出了章动参数计算公式中所用系数的最大的前5项。

在考虑章动对转换矩阵的影响后，则在考虑地球自转的格林威治恒星时时，使用真恒星时S_G而不是平恒星时

平恒星时的计算公式为：

${\stackrel{\&OverBar;}{S}}_G={18}^h.6973746+{879000}^h.0513367t+0^s.093104t^2-6^s.2\&times;{10}^{-6}{t}^3$

真恒星时S_G的计算公式为：

其中Δφcosε为赤经章动，是黄经章动Δφ在赤道上的分量。

上述公式中t为从J2000.0起算的儒略日的时间间隔，又称儒略世纪数，计算公式为：

$t=\frac{\mathrm{JD}\left(t\right)-\mathrm{JD}\left(J2000.0\right)}{36525.0}$

其中JD(t)是计算时刻对应的儒略日，JD(J2000.0)＝2451545.0是历元J2000.0对应的儒略日。

地球自转矩阵为(ER)。

在考虑岁差、章动前5项影响，从地心赤道惯性坐标系到地固坐标系的转换矩阵L_ECF,ECI为：

L_ECF,ECI＝(ER)(NR)(PR)

＝R_z(S_G)R_x(-Δε)R_y(Δθ)R_z(-Δμ)R_z(-z_A)R_y(θ_A)R_z(-ζ_A)

从地固坐标系到地心赤道惯性坐标系的转换矩阵L_ECI,ECF为：

L_ECI,ECF＝L_ECF,ECI′

c.根据卫星和成像点在地心赤道惯性坐标系中的位置，计算惯性系中卫星与成像点的相对位置矢量。图4为卫星与成像点相对位置的示意图，计算相对位置的具体方法为：

表示卫星在地心赤道惯性坐标系下的位置矢量，表示成像点在地心赤道惯性坐标系下的位置矢量。表示卫星与地面站之间的相对位置矢量，有如下关系：

$\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{r}}_{\mathrm{ECI}}={\stackrel{\&RightArrow;}{r}}_{\mathrm{earth}}-{\stackrel{\&RightArrow;}{r}}_{\mathrm{SAT}}$

即：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{r}_x=r_{x,\mathrm{earth}}-r_{x,\mathrm{SAT}}\\ \mathrm{\&Delta;}{r}_y=r_{y,\mathrm{earth}}-r_{y,\mathrm{SAT}}\\ \mathrm{\&Delta;}{r}_z=r_{z,\mathrm{earth}}-r_{z,\mathrm{SAT}}\end{array}\right.$

(3)根据地心赤道惯性坐标系到卫星轨道坐标系的一系列转换矩阵，计算相对位置矢量在卫星本体坐标系中的分量。计算卫星与地面成像点相对位置矢量在卫星轨道坐标系中的计算公式为：

$\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{r}}_{\mathrm{orbit}}=L_{\mathrm{oi}}\&CenterDot;\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{r}}_{\mathrm{ECI}}$

其中 $\begin{array}{c}{L}_{\mathrm{oi}}=L_y(-(\mathrm{\&pi;}/2+u))L_x(-(\mathrm{\&pi;}/2-i))L_z\left(\mathrm{\&Omega;}\right)\\ =\left[\begin{array}{ccc}-\cos \mathrm{\&Omega;}\sin u-\sin \mathrm{\&Omega;}\cos i\cos u& \cos \mathrm{\&Omega;}\cos i\cos u-\sin \mathrm{\&Omega;}\sin u& \cos u\sin i\\ -\sin \mathrm{\&Omega;}\sin i& \cos \mathrm{\&Omega;}\sin i& -\cos i\\ \sin \mathrm{\&Omega;}\cos i\sin u-\cos \mathrm{\&Omega;}\cos u& -\sin \mathrm{\&Omega;}\cos u-\cos \mathrm{\&Omega;}\cos i\sin u& -\sin i\sin u\end{array}\right]\end{array}.$

(4)根据相对位置矢量在卫星轨道坐标系中的分量，计算卫星成像时的滚转角、俯仰角、滚转角的欧拉角速度、俯仰角的欧拉角速度。

a.计算卫星成像时的滚转角和俯仰角以及滚转角速度和俯仰角速度。设卫星指向地面成像点所需转动的滚转角为俯仰角为θ。计算公式为：

$\mathrm{\&theta;}=\arcsin \frac{r_{x,\mathrm{orbit}}}{r_{\mathrm{orbit}}}$

其中 $r_{\mathrm{orbit}}=\sqrt{r_{x,\mathrm{orbit}}^2+r_{y,\mathrm{orbit}}^2+r_{z,\mathrm{orbit}}^2}.$

b.计算卫星的滚转角的欧拉角速度和俯仰角的欧拉角速度。

计算卫星滚转角的欧拉角速度w_x和俯仰角的欧拉角速度w_y。卫星姿态机动过程中，卫星的姿态角度变化是一个连续变化过程，因此，在计算滚转角的欧拉角速度和俯仰角的欧拉角速度时可以采用二次多项式插值的数值逼近方法进行计算。

二次多项式插值方法，首先需要获得三个不同时刻的角度数据，分别记作Angle₁、Angle₂、Angle₃，其对应的时刻分别为t₁、t₂、t₃。二次多项式插值得到的多项式公式如下：

设函数y＝f(x)，则二次多项式插值函数为：

$L_n\left(x\right)=\underset{i=0}{\overset{2}{\mathrm{\&Sigma;}}}{y}_i\underset{k\&NotEqual;i}{\overset{2}{\underset{k=0}{\mathrm{\&Pi;}}}}\left(\frac{x-x_k}{x_i-x_k}\right)$

针对滚转角，已知三个不同时刻的滚转角数据和其对应的时刻t₁、t₂、t₃，则任意的t_t时刻的滚转角的计算公式为：

上式为二次多项式，其对于求关于时间t的导数即为滚转角对应的滚转角的欧拉角速度其导数公式为：

在得到三个滚转角后，可根据上述二次插值方法计算各时刻滚转角对应的滚转角的欧拉角速度。

针对俯仰角，已知三个不同时刻的俯仰角数据θ₁、θ₂、θ₃和其对应的时刻t₁、t₂、t₃，则任意的t_t时刻的滚转角θ_t的计算公式为：

${\mathrm{\&theta;}}_t={\mathrm{\&theta;}}_1\left(\frac{t_t-t_2}{t_1-t_2}\right)\left(\frac{t_t-t_3}{t_1-t_3}\right)+{\mathrm{\&theta;}}_2\left(\frac{t_t-t_1}{t_2-t_1}\right)\left(\frac{t_t-t_3}{t_2-t_3}\right)+{\mathrm{\&theta;}}_3\left(\frac{t_t-t_1}{t_3-t_1}\right)\left(\frac{t_t-t_2}{t_3-t_2}\right)$

上式为二次多项式，其对于θ_t求关于时间t的导数即为俯仰角θ_t对应的俯仰角的欧拉角速度其导数公式为：

$\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_t=[\frac{2{\mathrm{\&theta;}}_1}{(t_1-t_2)(t_1-t_3)}+\frac{2{\mathrm{\&theta;}}_2}{(t_2-t_1)(t_2-t_3)}+\frac{2{\mathrm{\&theta;}}_3}{(t_3-t_1)(t_3-t_2)}]t_t\\ -[\frac{(t_2+t_3){\mathrm{\&theta;}}_1}{(t_1-t_2)(t_1-t_3)}+\frac{(t_1+t_3){\mathrm{\&theta;}}_2}{(t_2-t_1)(t_2-t_3)}+\frac{(t_1+t_2){\mathrm{\&theta;}}_3}{(t_3-t_1)(t_3-t_2)}]\end{array}$

在得到三个俯仰角后，可根据上述二次插值方法计算各时刻俯仰角对应的滚转角的欧拉角速度。

在得到滚转角的欧拉角速度和俯仰角的欧拉角速度，可以计算得到滚转角速度和俯仰角速度。滚转角速度和俯仰角速度将用于下一步偏航角和偏航角的欧拉角速度的计算。计算公式和卫星姿态转序有密切关系，具体计算方法可以参考《卫星姿态动力学与控制(1)》一书。在此卫星以1-2-3转序为例，给出欧拉角速率与姿态角速度的转换关系。

此时令偏航角为0，偏航角的欧拉角速度为0，则滚转角速度w_x、俯仰角速度w_y为：

滚转角、滚转角的欧拉角速度、俯仰角、俯仰角的欧拉角速度的计算流程如图5所示。

(5)根据得到的滚转角、滚转角的欧拉角速度、俯仰角和俯仰角的欧拉角速度，轨道参数以及对应的成像点位置计算卫星的偏航角与偏航角的欧拉角速度。

a.设卫星本体的Z轴在卫星本体坐标系中的矢量为 ${\stackrel{\&RightArrow;}{z}}_b=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],$ 根据已经得到的滚转角与俯仰角θ，计算卫星本体Z轴在地固坐标系中的矢量表示形式 ${\stackrel{\&RightArrow;}{z}}_{\mathrm{ECF}}=\left[\begin{array}{c}{z}_{x,\mathrm{ECF}}\\ z_{y,\mathrm{ECF}}\\ z_{z,\mathrm{ECF}}\end{array}\right],$ 的计算公式为：

${\stackrel{\&RightArrow;}{z}}_{\mathrm{ECF}}=L_{\mathrm{ECF},\mathrm{ECI}}\&CenterDot;L_{\mathrm{io}}\&CenterDot;L_{\mathrm{ob}}\&CenterDot;{\stackrel{\&RightArrow;}{z}}_b$

其中L_io＝L_oi ^T，上标T表示对矩阵进行转置；L_ECF,ECI为地心赤道坐标系到地固坐标系的转换矩阵；L_ob为卫星本体坐标系到卫星轨道坐标系的转换矩阵：

其中为滚动角，θ为俯仰角，ψ为偏航角，需要指出的是L_ob中偏航角ψ的数值为0，这里卫星的姿态转序为1-2-3。

b.根据成像点的地固坐标系中的位置和卫星在地固坐标系中的位置计算卫星到成像点的距离矢量h。

成像点为随地球运动的地球表面一点，卫星与成像点的几何关系如图6所示，点S表示星上相机，星下点摄影时成像点为T1，像平面垂直于ST1；在卫星一般姿态机动后，星地连线和相机视轴的夹角为φ，摄影点为T，像平面垂直于ST，ST方向就是本体系Z_b轴的方向。地心角为其中r是卫星地心距，R是成像点的地心距。h的数值为

c.计算在地心赤道惯性坐标系中，成像点相对于卫星的相对速度矢量，即地速矢量具体的计算公式为：

$\stackrel{\&RightArrow;}{v}={\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_e\&times;\stackrel{\&RightArrow;}{R}-[{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n\&times;\stackrel{\&RightArrow;}{r}+({\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n+{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_s)\&times;\stackrel{\&RightArrow;}{h}+{\stackrel{\&RightArrow;}{v}}_r]={\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_e\&times;\stackrel{\&RightArrow;}{R}-{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n\&times;\stackrel{\&RightArrow;}{R}-{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_s\&times;\stackrel{\&RightArrow;}{h}-{\stackrel{\&RightArrow;}{v}}_r$

其中为地球自转角速度矢量，是地心到成像点的矢量，是地心指向卫星的矢量，是轨道角速度矢量，其数值为 是卫星绝对速度的径向分量，其大小为 $v_r=\sqrt{\frac{\mathrm{\&mu;}}{p}}e\sin f,$ $\stackrel{\&RightArrow;}{\mathrm{\&omega;}}=\left[\begin{array}{c}{w}_x\\ w_y\\ 0\end{array}\right]$ 为卫星的姿态角速度。

d.根据上一步得到的地速以及地心赤道惯性坐标系到卫星轨道坐标系的转换矩阵L_oi，转换得到该地速在卫星轨道坐标系中的矢量其具体计算公式为：

${\stackrel{\&RightArrow;}{v}}_o=L_{\mathrm{oi}}\&CenterDot;\stackrel{\&RightArrow;}{v}$

e.根据卫星滚转角、俯仰角和卫星轨道坐标系到卫星本体坐标系的转换矩阵L_bo，计算地速在卫星本体坐标系中的矢量具体的计算公式为：

${\stackrel{\&RightArrow;}{v}}_b=L_{\mathrm{bo}}\&CenterDot;{\stackrel{\&RightArrow;}{v}}_o$

f.根据地速在卫星本体坐标系中的矢量，计算偏航角。由矢量 ${\stackrel{\&RightArrow;}{v}}_b=\left[\begin{array}{c}{v}_{\mathrm{xb}}\\ v_{\mathrm{yb}}\\ v_{\mathrm{zb}}\end{array}\right]$ 计算偏航角，具体计算公式为：

$\mathrm{\&psi;}=\mathrm{arctg}\frac{v_{\mathrm{yb}}}{v_{\mathrm{xb}}}$

g.计算偏航角速度，利用二次多项式插值的数值逼近方法进行计算。具体的计算方法同滚转角的欧拉角速度、俯仰角的欧拉角速度的计算公式为：

针对偏航角，已知三个不同时刻的偏航角数据ψ₁、ψ₂、ψ₃和其对应的时刻t₁、t₂、t₃，则任意的t_t时刻的偏航角ψ_t的计算公式为：

${\mathrm{\&psi;}}_t={\mathrm{\&psi;}}_1\left(\frac{t_t-t_2}{t_1-t_2}\right)\left(\frac{t_t-t_3}{t_1-t_3}\right)+{\mathrm{\&psi;}}_2\left(\frac{t_t-t_1}{t_2-t_1}\right)\left(\frac{t_t-t_3}{t_2-t_3}\right)+{\mathrm{\&psi;}}_3\left(\frac{t_t-t_1}{t_3-t_1}\right)\left(\frac{t_t-t_2}{t_3-t_2}\right)$

上式为二次多项式，其对于ψ_t求关于时间t的导数即为偏航角ψ_t对应的偏航角欧拉角速度其导数公式为：

$\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&psi;}}}_t=[\frac{2{\mathrm{\&psi;}}_1}{(t_1-t_2)(t_1-t_3)}+\frac{2{\mathrm{\&psi;}}_2}{(t_2-t_1)(t_2-t_3)}+\frac{2{\mathrm{\&psi;}}_3}{(t_3-t_1)(t_3-t_2)}]t_t\\ -[\frac{(t_2+t_3){\mathrm{\&psi;}}_1}{(t_1-t_2)(t_1-t_3)}+\frac{(t_1+t_3){\mathrm{\&psi;}}_2}{(t_2-t_1)(t_2-t_3)}+\frac{(t_1+t_2){\mathrm{\&psi;}}_3}{(t_3-t_1)(t_3-t_2)}]\end{array}$

在得到三个偏航角后，可计算各时刻偏航角对应的偏航角欧拉角速度。偏航角与偏航角速度计算流程如图7所示。

(6)根据得到的滚转角的欧拉角速度、俯仰角的欧拉角速度和偏航角欧拉角速度，计算得到滚转角的角速度、俯仰角的角速度和偏航角的角速度：

(7)根据得到的卫星成像时的滚转角、俯仰角、偏航角、滚转角的角速度、俯仰角的角速度、偏航角的角速度，利用星上姿态控制设备对卫星实时姿态规划。

Claims (4)

1.一种敏捷卫星成像过程姿态机动规划方法，其特征在于步骤如下：

(1)根据指定的地面成像条带的起始点地理位置与结束点地理位置，以及成像的起始时间与结束时间，按照星上计算步长计算获得每个时间点对应的成像点位置；所述的地理位置包括地理经度、地理纬度和海拔高度；所述的成像点位置表示为成像点在地固坐标系下的三个位置分量；

(2)获取每个时间点对应的卫星轨道参数，根据步骤(1)获得的每个时间点对应时刻的成像点位置，计算获得卫星与成像点在地心赤道惯性坐标系中的相对位置矢量；

(21)根据卫星的轨道六根数，计算卫星在地心赤道惯性坐标系中的位置；

(22)计算成像点在地心赤道惯性坐标系中的位置；

(23)根据卫星和成像点在地心赤道惯性坐标系中的位置，采用矢量计算方法计算得到地心赤道惯性系中卫星与成像点的相对位置矢量；

(3)计算步骤(2)获得的相对位置矢量在卫星本体坐标系中的分量；

(4)根据相对位置矢量在卫星轨道坐标系中的分量，计算卫星成像时的滚转角、俯仰角、滚转角的欧拉角速度、俯仰角的欧拉角速度；其中计算滚转角的欧拉角速度、俯仰角的欧拉角速度时，使用二次多项式插值的数值方法得到；

(5)根据得到的滚转角、滚转角的欧拉角速度、俯仰角、俯仰角的欧拉角速度、每个时间点对应的卫星轨道参数以及每个时间点对应的成像点位置，计算卫星的偏航角与偏航角的欧拉角速度；其中计算偏航角的欧拉角速度时，使用二次多项式插值的数值方法得到；

(6)根据得到的滚转角的欧拉角速度、俯仰角的欧拉角速度和偏航角欧拉角速度，计算得到滚转角的角速度、俯仰角的角速度和偏航角的角速度；

(7)根据得到的卫星成像时的滚转角、俯仰角、偏航角、滚转角的角速度、俯仰角的角速度、偏航角的角速度，利用星上姿态控制设备对卫星实施姿态规划。

2.根据权利要求1所述的一种敏捷卫星成像过程姿态机动规划方法，其特征在于：步骤(1)中获得每个时间点对应的成像点位置的具体方法为：设地面成像条带的起始点为A，其地理经纬度为Lon_A,d和Lat_A,d，对应的起始时间为t_A；设地面成像条带的结束点为B，其地理经纬度为Lon_B,d和Lat_B,d，对应的结束时间为t_B；星上计算时间的步长为Δt；地面各个成像点位置的具体计算步骤如下：

(11)计算A点和B点分别在地固坐标系ECF中的位置坐标 $A_{\mathrm{ECF}}=\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}\end{array}\right]$ 和 $B_{\mathrm{ECF}}=\left[\begin{array}{c}{x}_{B,\mathrm{ECF}}\\ y_{B,\mathrm{ECF}}\\ z_{B,\mathrm{ECF}}\end{array}\right];$

(12)根据A点、B点的地固坐标系中的位置坐标和地心原点O的位置坐标，确定OAB平面在地固坐标系中的平面方程和该平面的法线方程；

所述OAB平面的方程具体为：

$\begin{array}{c}\left|\left|\begin{array}{cccc}x& y& z& 1\\ x_{A,\mathrm{ECF}}& y_{A,\mathrm{ECF}}& z_{A,\mathrm{ECF}}& 1\\ x_{B,\mathrm{ECF}}& y_{B,\mathrm{ECF}}& z_{E,\mathrm{CF}}& 1\\ 0& 0& 0& 1\end{array}\right.\right|=(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})x\\ -(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})y\\ +(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})z\end{array}$

OAB平面的法线方向数分别为(y_A,ECFz_B,ECF-z_A,ECFy_B,ECF)、-(x_A,ECFz_B,ECF-z_A,ECFx_B,ECF)和(x_A,ECFy_B,ECF-y_A,ECFx_B,ECF)；

平面的法线方程为：

$\frac{x}{y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}}=\frac{y}{z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}}=\frac{z}{x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}$

根据法线方程，过平面法线上的一点的坐标为 $\left(\begin{array}{c}{y}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}\\ x_{A,\mathrm{ECF}}{t}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\end{array}\right);$

(13)在地固坐标系中求解法线与地固坐标系的z轴的夹角i；

平面OAB法线与地固坐标系的z轴两直线的夹角计算公式为：

$\cos i=\frac{(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}{\sqrt{{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2+{(z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}})}^2+{(x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2}}$

，且夹角i的范围为0°～180°；

(14)求解平面OAB与地球在地固坐标系中xy平面上的交线；

设平面OAB与地固坐标系xy平面的交线矢量为地固坐标系的z轴矢量为则有：

$\begin{array}{c}\stackrel{\&RightArrow;}{\mathrm{ON}}=\stackrel{\&RightArrow;}{\mathrm{OZ}}\&times;\stackrel{\&RightArrow;}{\mathrm{oh}}=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\&times;\left[\begin{array}{c}{y}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}-x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}\\ x_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}-y_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\end{array}\right]\\ =\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}}\\ 0\end{array}\right]\end{array}$

计算得到地固坐标系x轴到矢量的旋转角度Ω：

$\cos \mathrm{\&Omega;}=\frac{{\mathrm{ON}}_{x,\mathrm{ECF}}}{\left|\stackrel{\&RightArrow;}{\mathrm{ON}}\right|}=\frac{x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}{\sqrt{{(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2+{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2}}$

如果在地固坐标系下y轴的分量y_A,ECFz_B,ECF-z_A,ECFy_B,ECF小于0，则有：

如果在地固坐标系下y轴的分量y_A,ECFz_B,ECF-z_A,ECFy_B,ECF大于等于0，则有：

$\mathrm{\&Omega;}=\arccos \left(\frac{x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}}}{\sqrt{{(x_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{x}_{B,\mathrm{ECF}})}^2+{(y_{A,\mathrm{ECF}}{z}_{B,\mathrm{ECF}}-z_{A,\mathrm{ECF}}{y}_{B,\mathrm{ECF}})}^2}}\right)$

，且夹角Ω的范围为0°～360°；

其中ON_x,ECF表示为在地固坐标系下x轴的分量；

(15)建立中间坐标系S1；

建立方法为首先绕地固坐标系的z轴旋转角度Ω，然后绕此时地固坐标系的x轴旋转角度i，即可得到中间坐标系S1；两次坐标系旋转均符合右手螺旋定则；在中间坐标系S1中A点和B点的坐标为：

$\left[\begin{array}{c}{x}_{A,S1}\\ y_{A,S1}\\ z_{A,S1}\end{array}\right]=L_{S1,\mathrm{ECF}}\left[\begin{array}{c}{x}_{A,\mathrm{ECF}}\\ y_{A,\mathrm{ECF}}\\ z_{A,\mathrm{ECF}}\end{array}\right],\left[\begin{array}{c}{x}_{B,S1}\\ y_{B,S1}\\ z_{B,S1}\end{array}\right]=L_{S1,\mathrm{ECF}}\left[\begin{array}{c}{x}_{B,\mathrm{ECF}}\\ y_{B,\mathrm{ECF}}\\ z_{B,\mathrm{ECF}}\end{array}\right]$

其中L_S1,ECF为：

$\begin{array}{c}{L}_{S1,\mathrm{ECF}}=L_x\left(i\right)L_z\left(\mathrm{\&Omega;}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos i& \sin i\\ 0& -\sin i& \cos i\end{array}\right]\&CenterDot;\left[\begin{array}{ccc}\cos \mathrm{\&Omega;}& \sin \mathrm{\&Omega;}& 0\\ -\sin \mathrm{\&Omega;}& \cos \mathrm{\&Omega;}& 0\\ 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{ccc}\cos \mathrm{\&Omega;}& \sin \mathrm{\&Omega;}& 0\\ -\cos i\sin \mathrm{\&Omega;}& \cos i\cos \mathrm{\&Omega;}& \sin i\\ \sin i\sin \mathrm{\&Omega;}& -\sin i\cos \mathrm{\&Omega;}& \cos i\end{array}\right]\end{array};$

(16)建立坐标系S2；

计算坐标系S1下x轴到A点转动角度ω；具体方法为：

设A点在中间坐标系S1的坐标系为(x_A,S1,y_A,S1,z_A,S1)′，有：

$\cos \mathrm{\&omega;}=\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}$

如果z_A,ECF<0，则 $\mathrm{\&omega;}=360-\arccos \left(\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}\right),$ 如果 $z_{A,\mathrm{ECF}}\&GreaterEqual;0\mathrm{\&omega;}=\arccos \left(\frac{x_{A,S1}}{\sqrt{x_{A,S1}^2+y_{A,S1}^2+z_{A,S1}^2}}\right);$ 且夹角ω的范围为0°～360°；

绕坐标系S1的z轴转动角度ω后，得到坐标系S2；坐标系S2中A点和B点的坐标计算公式为：

$\left[\begin{array}{c}{x}_{A,S2}\\ y_{A,S2}\\ z_{A,S2}\end{array}\right]=L_{S2,S1}\left[\begin{array}{c}{x}_{A,S1}\\ y_{A,S1}\\ z_{A,S1}\end{array}\right],\left[\begin{array}{c}{x}_{B,S2}\\ y_{B,S2}\\ z_{B,S2}\end{array}\right]=L_{S2,S1}\left[\begin{array}{c}{x}_{B,S1}\\ y_{B,S1}\\ z_{B,S1}\end{array}\right]$

其中L_S2,S1为：

$L_{S2,S1}=L_z\left(\mathrm{\&omega;}\right)=\left[\begin{array}{ccc}\cos \mathrm{\&omega;}& \sin \mathrm{\&omega;}& 0\\ -\sin \mathrm{\&omega;}& \cos \mathrm{\&omega;}& 0\\ 0& 0& 1\end{array}\right]$

(17)计算成像观测条带中A、B两点之间，卫星每隔Δt时间观测的地面点；

根据在S2坐标系中A、B两点的坐标，计算A、B两点之间的夹角α：

$\cos \mathrm{\&alpha;}=\frac{x_{A,S2}{x}_{B,S2}+y_{A,S2}{y}_{B,S2}+z_{A,S2}{z}_{B,S2}}{\sqrt{x_{A,S2}^2+y_{A,S2}^2+z_{A,S2}^2}\sqrt{x_{B,S2}^2+y_{B,S2}^2+z_{B,S2}^2}}$

得到A、B两点之间的夹角α后，根据卫星星上计算的时间间隔Δt，和成像起始时间t_A和结束时间t_B，确定观测的地面点个数：

$\mathrm{sum}=\frac{t_B-t_A}{\mathrm{\&Delta;t}}$

则卫星成像过程中指向的每个成像点在球体模型上大圆的角度间隔为 $\mathrm{\&Delta;\&alpha;}=\frac{\mathrm{\&alpha;}}{\mathrm{sum}};$

设A点的地心距为则A点在坐标系S2中的单位向量为 $\left[\begin{array}{ccc}{\tilde{x}}_{A,S2}& {\tilde{y}}_{A,S2}& {\tilde{z}}_{A,S2}\end{array}\right]={\left[\begin{array}{ccc}\frac{x_{A,S2}}{r_A}& \frac{y_{A,S2}}{r_A}& \frac{z_{A,S2}}{r_A}\end{array}\right]}^{\&prime;},$ 继而有A点在坐标系S2中的单位向量的模为 ${\tilde{r}}_A=\sqrt{{\tilde{x}}_{A,S2}^2+{\tilde{y}}_{A,S2}^2+{\tilde{z}}_{A,S2}^2}=1;$

设成像观测条带中的成像点为A，P1，P2，……，Pn，……，B，在坐标系S2中的每个成像点位置的单位向量有：

$\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},S2}\\ {\tilde{y}}_{\mathrm{Pn},S2}\\ {\tilde{z}}_{\mathrm{Pn},S2}\end{array}\right]=\left[\begin{array}{c}{\tilde{r}}_A\cos (n\&times;\mathrm{\&beta;})\\ {\tilde{r}}_A\sin (n\&times;\mathrm{\&beta;})\\ 0\end{array}\right],(n=\mathrm{0,1},...,\mathrm{sum}),\left\{\begin{array}{cc}\mathrm{\&beta;}=\mathrm{\&Delta;\&alpha;},& y_{B,S2}>0\\ \mathrm{\&beta;}=-\mathrm{\&Delta;\&alpha;},& y_{B,S2}<0\end{array}\right.$

得到Pn在坐标系S2中单位向量的位置坐标后，通过坐标系转换，计算其在地固坐标系中的位置坐标：

$\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}\\ {\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}\\ {\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}\end{array}\right]=L_z(-\mathrm{\&Omega;})L_x(-i)L_z(-\mathrm{\&omega;})\left[\begin{array}{c}{\tilde{x}}_{\mathrm{Pn},S2}\\ {\tilde{y}}_{\mathrm{Pn},S2}\\ {\tilde{z}}_{\mathrm{Pn},S2}\end{array}\right]$

其中 $L_z(-\mathrm{\&Omega;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&Omega;})& \sin (-\mathrm{\&Omega;})& 0\\ -\sin (-\mathrm{\&Omega;})& \cos (-\mathrm{\&Omega;})& 0\\ 0& 0& 1\end{array}\right],L_x(-i)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos (-i)& \sin (-i)\\ 0& -\sin (-i)& \cos (-i)\end{array}\right],$ $L_z(-\mathrm{\&omega;})=\left[\begin{array}{ccc}\cos (-\mathrm{\&omega;})& \sin (-\mathrm{\&omega;})& 0\\ -\sin (-\mathrm{\&omega;})& \cos (-\mathrm{\&omega;})& 0\\ 0& 0& 1\end{array}\right];$

得到在地固坐标系中的分量后，利用空间解析几何，求解直线与空间椭球的交点：

设过原点和的空间直线为设地球椭球模型为其中a＝Re＝6378140m为地球赤道半径，c＝Rp＝6356755m为地球极半径；

根据空间解析几何，得到空间直线的两点式方程为：

$\frac{x-0}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}-0}=\frac{y-0}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}-0}=\frac{z-0}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}-0}\&DoubleRightArrow;\frac{x}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{y}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{z}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}}$

将的空间直线方程与地球椭球方程联立得到： $\left\{\begin{array}{c}\frac{x}{{\tilde{x}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{y}{{\tilde{y}}_{\mathrm{Pn},\mathrm{ECF}}}=\frac{z}{{\tilde{z}}_{\mathrm{Pn},\mathrm{ECF}}}\\ \frac{x^2+y^2}{a^2}+\frac{z^2}{c^2}=1\end{array}\right.,$ 求解得到的两组解[x₁ y₁ z₁]′和[x₂ y₂ z₂]′，分别求两个解与成像起始点A之间的夹角 $\cos {\mathrm{\&gamma;}}_1=\frac{x_{A,\mathrm{ECF}}{x}_1+y_{A,\mathrm{ECF}}{y}_1+z_{A,\mathrm{ECF}}{z}_1}{\sqrt{x_{A,\mathrm{ECF}}^2+y_{A,\mathrm{ECF}}^2+z_{A,\mathrm{ECF}}^2}\sqrt{x_1^2+y_1^2+z_1^2}}$ 和 $\cos {\mathrm{\&gamma;}}_2=\frac{x_{A,\mathrm{ECF}}{x}_2+y_{A,\mathrm{ECF}}{y}_2+z_{A,\mathrm{ECF}}{z}_2}{\sqrt{x_{A,\mathrm{ECF}}^2+y_{A,\mathrm{ECF}}^2+z_{A,\mathrm{ECF}}^2}\sqrt{x_2^2+y_2^2+z_2^2}}$ 与成像起始点A夹角小的即为成像点P_n在地固坐标系中的位置。

3.根据权利要求1所述的一种敏捷卫星成像过程姿态机动规划方法，其特征在于：步骤22)中计算时利用地固坐标系到地心赤道惯性坐标系的转换矩阵，且转换矩阵中只考虑岁差修正以及章动的前五项修正。

4.根据权利要求1所述的一种敏捷卫星成像过程姿态机动规划方法，其特征在于：步骤(5)中计算偏航角的具体方法为：

(51)定义卫星本体系Z轴指向地心，X轴在轨道平面内垂直z轴指向卫星飞行方向，Y轴满足右手定则；设卫星本体的Z轴在卫星本体坐标系中的矢量为 ${\stackrel{\&RightArrow;}{z}}_b=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right],$ 根据已经得到的滚转角与俯仰角θ，计算卫星本体Z轴在地固坐标系中的矢量表示形式 ${\stackrel{\&RightArrow;}{z}}_{\mathrm{ECF}}=\left[\begin{array}{c}{z}_{x,\mathrm{ECF}}\\ z_{y,\mathrm{ECF}}\\ z_{z,\mathrm{ECF}}\end{array}\right],$ 的计算公式为：

${\stackrel{\&RightArrow;}{z}}_{\mathrm{ECF}}=L_{\mathrm{ECF},\mathrm{ECI}}\&CenterDot;L_{\mathrm{io}}\&CenterDot;L_{\mathrm{ob}}\&CenterDot;{\stackrel{\&RightArrow;}{z}}_b$

其中L_io＝L_oi ^T，上标T表示对矩阵进行转置；L_ECF,ECI为地心赤道坐标系到地固坐标系的转换矩阵；L_ob为卫星本体坐标系到卫星轨道坐标系的转换矩阵；

(52)根据成像点在地固坐标系中的位置和卫星在地固坐标系中的位置计算卫星到成像点的距离矢量h，其数值为其中地心角为r是卫星地心距，R是成像点的地心距；

(53)计算在地心赤道惯性坐标系中，成像点相对于卫星的相对速度矢量，即地速矢量具体的计算公式为：

$\stackrel{\&RightArrow;}{v}={\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_e\&times;\stackrel{\&RightArrow;}{R}-[{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n\&times;\stackrel{\&RightArrow;}{r}+({\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n+{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_s)\&times;\stackrel{\&RightArrow;}{h}+{\stackrel{\&RightArrow;}{v}}_r]={\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_e\&times;\stackrel{\&RightArrow;}{R}-{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_n\&times;\stackrel{\&RightArrow;}{R}-{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_s\&times;\stackrel{\&RightArrow;}{h}-{\stackrel{\&RightArrow;}{v}}_r$

其中为地球自转角速度矢量，是地心到成像点的矢量，是地心指向卫星的矢量，是轨道角速度矢量，其数值为 是卫星绝对速度的径向分量，其大小为 为卫星的姿态角速度；

(54)根据上一步得到的地速以及地心赤道惯性坐标系到卫星轨道坐标系的转换矩阵L_oi，转换得到该地速在卫星轨道坐标系中的矢量其具体计算公式为：

${\stackrel{\&RightArrow;}{v}}_o=L_{\mathrm{oi}}\&CenterDot;\stackrel{\&RightArrow;}{v}$

(55)根据卫星滚转角、俯仰角和卫星轨道坐标系到卫星本体坐标系的转换矩阵L_bo，计算地速在卫星本体坐标系中的矢量具体的计算公式为：

${\stackrel{\&RightArrow;}{v}}_b=L_{\mathrm{bo}}\&CenterDot;{\stackrel{\&RightArrow;}{v}}_o$

(56)根据地速在卫星本体坐标系中的矢量 ${\stackrel{\&RightArrow;}{v}}_b=\left[\begin{array}{c}{v}_{\mathrm{xb}}\\ v_{\mathrm{yb}}\\ v_{\mathrm{zb}}\end{array}\right],$ 计算获得偏航角： $\mathrm{\&psi;}=\mathrm{arctg}\frac{v_{\mathrm{yb}}}{v_{\mathrm{xb}}}.$