低轨光学卫星同轨立体成像姿态规划与控制方法

随着航天技术的发展，光学卫星的总体设计、姿态确定方法、控制以及成像参数优化等技术与方法得到深入的研究与长足的进步，并使得利用光学卫星获得遥感影像成为大范围地理测绘的重要手段。随着遥感与测绘技术的发展，高分辨率立体影像得到广泛应用^［1-4］。

卫星多次对同一地面目标进行不同角度观测的成像方式可实现对目标的立体测绘，获取三维影像。专用立体成像测绘卫星采用多线阵相机进行连续长条带立体成像^［5］，经济成本高、研制难度大。微小卫星体积小、重量轻，可采用单相机通过卫星俯仰轴的快速姿态机动实现同轨立体成像或通过多次过境进行异轨立体成像。其中，同轨立体成像缩短了立体像对的获取时间间隔，影像环境差异小，方便测图处理^［6-7］，成为获取立体像对的趋势。王中果等推导出基于圆形地球模型星下点轨迹下的立体成像条带长度通用公式和最小机动能力需求，为简化分析并不考虑地球椭率、自转以及滚动方向的姿态机动^［5］。张新伟等采用平面地球模型，对同轨立体成像时不同重叠条带长度对所需姿态机动能力进行了分析^［8］，未考虑地球曲率的影响，导致卫星成像姿态求解结果存在偏差。田原等从立体像对后续处理角度出发，利用基高比和几何分辨率差异两方面构建同轨立体成像规划目标函数，提出了一种基于步进搜索策略的成像规划方法，决策多次拍摄的时间窗口^［6］，为已确定成像参数的同轨立体成像时间窗口计算提供研究思路。

在卫星幅宽确定的情况下，条带长度决定了成像覆盖能力，因此要求卫星具有三轴大角度姿态快速机动能力，现有研究主要面向姿态规划与控制算法两个方向。针对侧摆方向大角度机动的姿态规划研究，资源三号与高分三号卫星分别设计了三段梯形和五段梯形角速度的姿态规划^［9-10］，但角速度存在突变。申晓宁等针对加速度突变问题提出了基于S型角速度曲线的机动姿态规划方法^［11］。Ha和Femiano等介绍了角加速度为连续曲线的姿态机动规划方法^［12］。面向控制算法的研究如：Li和Bainum等学者认为经典的欧拉定轴转动是最优控制或接近最优控制^［13］；华冰等设计了基于欧拉轴转动的递阶饱和姿态控制律^［14］，但机动过程存在快速性与稳定性的矛盾。汪礼成等的研究表明引入加速度前馈能够改善控制性能^［15］。韩京清在自抗扰控制中详细说明了通过姿态规划可以解决传统控制存在的快速性与超调性矛盾问题^［16］。

本文考虑地球椭率、自转以及滚动方向的姿态机动，进行同轨立体成像姿态规划和成像时间窗口确定，提高规划的精度；进一步考虑微小卫星执行机构输出力矩和角动量约束，设计了路径规划快速机动控制（Path planning Fast Maneuver Control，PFMC）算法，基于旋转轴不变约束设计了三轴机动的最短路径和角加速度连续的姿态引导路径规划器，同时设计了结合角加速度前馈与内外环控制的快速机动算法。最后通过仿真与在轨试验验证了同轨双视立体成像的姿态规划与PFMC算法的可行性和有效性。

低轨光学微小卫星进行同轨立体成像时，在同一轨道圈次内通过卫星三轴快速机动可对同一目标进行不同角度的多次观测，如图1所示。双视立体成像通过前视和后视两次成像实现，三视立体成像则须增加一次正视成像，但都要求重复覆盖，才可生产出立体影像。

图1  低轨光学卫星同轨立体成像示意图

Fig.1  Schematic diagram of along-track stereoscopic imaging for low orbit optical satellite

立体成像方案分为成像姿态规划和PFMC两部分。其中成像姿态规划用于高精度的成像姿态与成像时间计算，PFMC分为三轴路径规划设计与快速机动控制，卫星同轨立体成像控制流图如图2所示。

图2  同轨立体成像控制流图

Fig.2  Control flow diagram of along-track stereoscopic imaging

成像姿态规划中，为了提高立体成像姿态求解的精度，采用地球椭球模型确定成像点的位置，并考虑地球的自转确定成像点与卫星的相对位置与速度信息，进行含有侧摆机动的立体成像姿态求解。通过对沿卫星飞行方向的特性分析，获取成像窗口时间及相应的成像姿态四元数，详细计算过程如下。

2.1　卫星成像姿态求解

采用椭球形地球模型，考虑地球椭率的影响，由成像点的地理经度θlon、地理纬度θlat和高度h信息，求得在地固系下的位置向量P=(Px,Py,Pz)T，示意图如图3所示。

图3  卫星位置与成像点关系示意图

Fig.3  Schematic diagram of relationship between satellite position and imaging points

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪Px=(r1−e2sin2 (θlat)√+h)cos (θlat)cos (θlon)Py=(r1−e2sin2 (θlat)√+h)cos (θlat)sin(θlon)Pz=(r(1−e2)1−e2sin2 (θlat)√+h)sin (θlat)     $$$$

（1）

其中：r为赤道半径，e为偏心率。

通过卫星导航数据和轨道递推公式得到地固系下卫星的实时位置向量S，对应的沿轨道坐标系(OoXoYoZo)中Zo轴反方向的单位向量为：

eS=S∥S∥$$$$

（2）

其中，∥S∥向量S的模。

地固系下地球自转角速度为：

ωEarth=[00ωEe]T$$$$

（3）

其中，ωEe=7.292 115×10−5 rad/s。

考虑地球自转的影响，若卫星速度向量为V，卫星相对成像点的相对速度为：

Vs=V+ωEarth×S$$$$

（4）

eVs=Vs∥Vs∥$$$$

（5）

以卫星的位置、速度和成像点的经度、纬度、高度以及地球自转角速度作为输入，进行卫星的姿态求解，求解流程如图4所示。

图4  卫星成像姿态计算流程图

Fig.4  Flow chart of satellite image attitude calculation

Ce，P，S三点共面，卫星成像时相机指向的单位向量ek为：

ek=S−P∥S−P∥$$$$

（6）

平面CePS的法线向量en在地固系的表示为：

en=eS×ek∥eS×ek∥$$$$

（7）

由欧拉旋转定理知，向量eS以向量en为旋转轴旋转角度φ得到向量ek，旋转角度为：

φ=arcsin(∥eS×ek∥)$$$$

（8）

由S，k，n，Vs计算轨道坐标系三轴在地固系表示的单位向量eox，eoy，eoz分别为：

⎧⎩⎨⎪⎪⎪⎪⎪⎪eoy=eVs×eS∥eVs×eS∥eoz=−eS  eox=eoy×eoz∥eoy×eoz∥$$$$

（9）

旋转轴en在轨道系下表示enox，enoy，enoz为：

⎧⎩⎨⎪⎪enox=∥eoxen∥enoy=∥eoyen∥enoz=∥eozen∥$$$$

（10）

旋转角φ在轨道系下表示φnox，φnoy，φnoz为：

⎧⎩⎨⎪⎪⎪⎪φnox=enoxφφnoy=enoyφφnoz=enozφ$$$$

（11）

成像时卫星本体姿态相对轨道系旋转四元数：

qA=[cos (φ2)sin (φ2)eno]T$$$$

（12）

卫星成像时本体姿态相对惯性系旋转四元数：

qf=qOI⊗qA$$$$

（13）

其中，qOI为轨道系相对惯性系旋转四元数。

2.2　成像时间窗口确定

以同轨立体双视和三视为例，通过对沿卫星飞行方向的特性进行分析，获取成像时间窗口。沿轨方向的旋转角度为φnoy，立体成像前后视角度为θ，即前视成像要求φnoy=θ，后视成像要求φnoy=−θ。

2.2.1　双视模式

双视模式下卫星成像过程分为前视与后视两部分，如图5所示。

图5  双视立体成像机动时间序列示意图

Fig.5  Time series diagram of imaging and maneuver for double perspective stereoscopic imaging

双视成像的成像与机动时间关系表如表1。当观测角度θ确定后，时间Ttrans+Timage为定值，通过提高卫星的机动能力，Ttrans减小为Ttrans−ΔT，则成像时间Timage增加为Timage+ΔT。

表1  双视立体成像与机动时间关系表

Tab.1  Relationship between imaging and maneuver time for double perspective stereoscopic imaging

Item	Time quantum	Imaging time	Maneuver time
Front view （θ$\theta$）	A1~A2$A_1~A_2$	Timage$T_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}$	/
Transform	A2~B1$A_2~B_1$	/	Ttrans$T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$
Rear view（−θ$-\theta$）	B1~B2$B_1~B_2$	Timage$T_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}$	/

2.2.2　三视模式

三视模式可对目标区域进行三次成像，以获得三个视角的影像。三视成像过程分为前视、正视与后视三部分，其中前视与后视成像角度大小通常相等，如图6所示。

图6  三视立体成像机动时间序列示意图

Fig.6  Time series diagram of imaging and maneuver for triple perspective stereoscopic imaging

三视成像的成像与机动时间关系如表2，对比双视成像可知三视成像要求卫星具有更强的机动能力。

表2  三视成像与机动时间关系表

Tab.2  Relationship between imaging and maneuver time for triple perspective stereoscopic imaging

Item	Time quantum	Imaging time	Maneuver time
Front view （θ$\theta$）	A1~A2$A_1~A_2$	Timage$T_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}$	/
Transform1	A2~B1$A_2~B_1$	/	Ttrans1$T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}1}$
Nadir view	B1~B2$B_1~B_2$	Timage$T_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}$	/
Transform2	B2~C1$B_2~C_1$	/	Ttrans2$T_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}2}$
Rear view （−θ$-\theta$）	C1~C2$C_1~C_2$	Timage$T_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{e}}$	/

由于成像时前视和后视角度的存在，卫星的成像时间并非成像点的过境时间窗口，需要通过轨道递推并根据实时旋转角度的分量与前后视角度的关系进行反向求解获取。三视成像时间窗口确定流程图如图7所示。

图7  三视成像时间窗口计算流程图

Fig.7  Flow chart to calculate the imaging time window for triple perspective stereoscopic imaging

具体获取过程为：

Step 1.设卫星前、正、后视成像时间初值T1=T2=T3=1020，初始时间T=T0。

Step 2.轨道递推得到时间T对应的卫星位置、速度信息。

Step 3.由卫星位置和速度信息、成像点的经纬高以及地球自转角速度，经过卫星成像姿态求解，获得时间T对应的旋转角度在卫星飞行方向的分量φnoy。

Step 4.判定T，若T<T1，转至Step 5判前视成像；若T1<T<T2，转至Step 6判正视成像；若T2<T<T3转至Step 7判后视成像。

Step 5.判断前视角度φnoy=θ，如果成立，获取对应的前视成像时间T01与成像姿态qf1，修改成像时间值T1=T01，转至Step 6。否则直接转至Step 8。

Step 6.判断正视角度φnoy=0，如果成立，获取对应的正视成像时间T02与成像姿态qf2，修改成像时间值T2=T02，转至Step 5。否则直接转至Step 8。

Step 7.判断后视角度φnoy=−θ，如果成立，获取对应的后视成像时间T03与成像姿态qf3，转至Step 9。否则转至Step 8。

Step 8.轨道进行递推，递推时间增加ΔT，T=T+ΔT，转至Step 2。

Step 9.姿态求解结束，输出结果。

对于双视成像，前视成像时间为T04，成像姿态为qf4，后视成像时间为T05，成像姿态为qf5。

同轨立体成像要求卫星快速进行三轴大角度姿态机动，同时通过提高卫星的机动性能可延长同轨立体成像的成像时长，提高成像覆盖能力。而微小卫星的执行机构力矩和角动量的大小是有限制的，导致卫星机动角速度受限。首先，考虑执行机构的约束设计一种基于旋转轴不变约束的最短路径方法，然后设计一种角加速度连续的八段式路径规划器，最终，设计了角加速度前馈结合内外环控制的快速机动算法，提高姿态跟踪的动态性能。

3.1　旋转轴/角计算

目标坐标系 Ox'y'z'相对初始坐标系 Oxyz的姿态转换，如图8所示，用四元数形式表示为：

q=[q0q¯ ]T=[cos φ2 sin φ2enT ]T$$$$

（14）

其中：q20+q¯2=1，en=[ij k]T，i  2+j2  +   k2=1。

图8  两坐标系间的旋转示意图

Fig.8  Schematic of rotation between two coordinate systems

卫星的目标姿态相对初始姿态的期望四元数：

qQ=q−1C⊗qF$${\mathit{q}}_{\mathrm{Q}}={\mathit{q}}_{\mathrm{C}}^{-\mathrm{1}}\otimes {\mathit{q}}_{\mathrm{F}}$$ ，

（15）

其中：初始四元数qC为卫星的初始姿态相对惯性系的旋转四元数，目标四元数qF为卫星的目标姿态相对惯性系的旋转四元数。

由qQ反向求得旋转角θQ和旋转轴en：

{θA=2arccos(qQ0)en=q¯Q/sin θQ2    ,θQ≠0$$$$

（16）

当θQ=0时，对应的四元数为qQ=[1000]T，目标姿态与初始姿态重合，不需进行姿态转换。

3.2　旋转轴不变约束求解

在卫星姿态机动的过程中，角加速度、角速度以及控制力矩受执行机构性能的限制。考虑卫星的性能约束通过路径规划器对卫星的机动过程进行实时规划，提高卫星的机动能力。

在满足卫星转动惯量I、反作用飞轮力矩T和角动量H约束下，为了实现初始四元数qC至目标四元数qF的最短路径旋转，即整个期望四元数qQ的机动过程中均保证旋转轴en不变，在进行路径规划时，需要对三轴角加速度、角速度进行限幅和约束设置。

路径规划器的输入为期望旋转角度θQ${\theta }_{\mathrm{Q}}$，角加速度限值αLG${\alpha }_{\mathrm{L}\mathrm{G}}$，角速度限值ωLG${\omega }_{\mathrm{L}\mathrm{G}}$，输出为实时角度，实时角速度，PFMC三轴路径规划如图9所示。

图9  PMFC三轴路径规划图

Fig.9  Schematic diagram of PMFC triaxial path planning

角加速度限值αLG${\alpha }_{\mathrm{L}\mathrm{G}}$的计算如下：

α¯¯i={(I−1T)i/ei ,ei≠0 Mmax     ,ei=0$$$$

（17）

其中，Mmax=1020$M_{\mathrm{m}\mathrm{a}\mathrm{x}}={\mathrm{10}}^{\mathrm{20}}$，为一个较大的数。

αLG=[α¯¯xα¯¯yα¯¯z]Tminen$$$$

（18）

αLG=∥αLG∥$$$$

（19）

其中，[α¯¯xα¯¯yα¯¯z]Tmin${\left[\begin{array}{ccc}{\overline{\alpha }}_{\mathit{x}}& {\overline{\alpha }}_{\mathit{y}}& {\overline{\alpha }}_{\mathit{z}}\end{array}\right]}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{T}}$为求向量[α¯¯xα¯¯yα¯¯z]T${\left[\begin{array}{ccc}{\overline{\alpha }}_{\mathit{x}}& {\overline{\alpha }}_{\mathit{y}}& {\overline{\alpha }}_{\mathit{z}}\end{array}\right]}_{}^{\mathrm{T}}$的最小值。

角速度限值ωLG${\omega }_{\mathrm{L}\mathrm{G}}$的计算如下：

ω¯¯¯i={[(I−1H)i,(ωLim)i]minei,Mmax,ei≠0ei=0$$$$

（20）

其中，ωLim=[ωLimxωLimyωLimz]T${\mathit{\omega }}_{\mathrm{L}\mathrm{i}\mathrm{m}}={\left[\begin{array}{ccc}{\omega }_{\mathrm{L}\mathrm{i}\mathrm{m}x}& {\omega }_{\mathrm{L}\mathrm{i}\mathrm{m}y}& {\omega }_{\mathrm{L}\mathrm{i}\mathrm{m}z}\end{array}\right]}^{\mathrm{T}}$。

ωLG=[ω¯¯xω¯¯yω¯¯z]Tminen$$$$

（21）

ωLG=∥ωLG∥$$$$

（22）

其中，飞轮角动量H=[HxHyHz]T$\mathit{H}={\left[\begin{array}{ccc}{H}_x& H_y& H_z\end{array}\right]}^{\mathrm{T}}$。