对地观测卫星在众多领域都有着广泛的应用。在搭载可见光载荷的卫星中，太阳同步轨道是常被采用的卫星轨道类型。单颗对地观测卫星往往通过姿态机动的方法提高观测能力，可以通过定量的方法对卫星姿态机动的敏捷性进行评估^[1]。随着对地观测任务向着不断复杂化的方向发展，单颗卫星已经不能满足任务需求。卫星星座是指一个由多颗卫星构成的用于完成某项任务的空间系统^[2]。越来越多的对地观测任务需要卫星星座来完成。同时，为了增强星座的覆盖能力，卫星的轨道机动能力也被充分利用。为实现快速响应、多星协同的目标，对于由多颗具备快速机动能力的太阳同步轨道卫星组成的卫星星座，研究其成像任务规划的重要性日益凸显。

关于对地成像任务规划，国内外学者进行了大量研究^[3]。可见性计算方法方面，鄂智博等^[4]提出了一种快速计算方法，具备很高的效率和精度。任务规划算法方面，姜维等^[5]将对地成像任务规划算法分为确定性算法、启发式算法和智能搜索算法等。其中，启发式算法和智能搜索算法在解决大规模问题上更具优势，被多数学者采用。Chu等^[6]在双星任务规划问题中，提出了分支限界(branch and bound)算法。Tangpattanakul等^[7]针对单颗敏捷卫星，提出了一种基于指标的多目标局部搜索(indicator-based multi-objective local search, IBMOLS)。其优化目标为最大化观测任务收益，通过多目标局部搜索，在不断迭代过程中逐步缩小解空间，并配合遗传算法使得优化结果快速收敛。Li等^[8]在模拟退火算法的基础上进行改进，通过结合模拟退火算法与遗传算法提出了混合遗传算法，证明了其优于单纯的模拟退火算法。何磊等^[9]考虑了光学成像任务中的云层遮挡问题，并利用蚁群算法进行求解。Niu等^[10]针对自然灾害观测问题，建立了多目标整数规划模型，并提出了一种遗传算法(non-dominated sorting genetic algorithm-Ⅱ, NSGA-Ⅱ)，对平均覆盖率、成像时间等多指标进行优化。宋彦杰等^[11]提出了一种改进遗传算法，能够根据种群情况切换全局、局部优化，从而实现任务利益的最大化。Lee等^[12]在其研究中针对多星的任务规划也使用了遗传算法。Zheng等^[13]对卫星自主任务重规划方法进行了研究。在复杂的观测任务中，往往还需使卫星具备更强的姿态机动和轨道机动能力。部分学者考虑了卫星的姿态机动能力，研究了敏捷卫星的成像任务规划问题^[14-16]。Zhu等^[17]以最大观测覆盖时间及最少燃料消耗为优化指标，研究了应急观测任务中的卫星轨道机动优化问题，对不同的轨道机动方式进行了比较。盛婧等^[18]提出了脉冲轨道机动情况下的一种星下点调整方法。Xue等^[19]则研究了航天器在进行单脉冲轨道机动情况下的可达域问题。可见，现有研究集中于对任务规划算法的讨论，主要采用不同的启发式算法和智能搜索算法，或在原有算法基础上进行改进、结合、创新；然而较少考虑卫星的轨道机动能力。

针对机动卫星星座对多个地面点目标的成像任务规划问题，本文考虑了卫星的快速轨道机动能力，首先，通过结合解析和数值方法，提出单星单目标观测方法，并给出变轨策略；其次，利用遗传算法提高任务规划问题的求解效率；最后，通过仿真算例验证算法的可行性和效果。

1 相关约束及指标

本研究针对特定的星座构型及成像条件，相关约束包括：

1) 卫星星座由多颗太阳同步圆轨道卫星组成；

2) 每颗卫星的单星总速度增量约束为定值；

3) 卫星进行轨道机动后需在规定时间内返回原有构型；

4) 卫星上搭载载荷为光学相机，仅在光照情况下完成对地成像；

5) 卫星通过姿态机动实现相机侧摆，姿态机动须满足最大姿态机动角限制；

6) 不考虑太阳同步轨道卫星无法观测的两极区域；

7) 仅考虑地球J₂摄动。

工程实际中，尽管星座能够在长时间内覆盖全球，短时间(如24 h)内仍存在无法观测到的部分目标。在卫星星座对多目标成像任务规划过程中，通过轨道机动使卫星具备观测到尽可能多目标的能力，同时使所有目标的观测次数分布更加均匀。

由于星座已经具备对全球绝大部分区域的成像能力，考虑到工程实际，只考虑每颗卫星最多进行一次相位机动(不包括返回构型)的情况。

2 单星单目标观测方法

结合解析和数值2种方法，完成单颗卫星对单个地面点目标可见性分析。对于不可见的目标，给出变轨策略。

2.1 解析方法

如图 1所示，建立如下地心惯性坐标系：X轴由地心指向轨道升交点，Z轴沿地轴指向北极，Y轴构成右手正交系。其中S、T分别为卫星、目标初始时刻的位置。

假设卫星观测目标的一种理想情况：当目标处于卫星所在轨道平面内时，恰好与卫星星下点重合，即某时刻卫星和目标恰好分别运动到S^′、T^′。此时，卫星不需要进行变轨或调整视场即可对目标进行直接观测。

对于一般的情况，目标和卫星运行总会存在一个时间差，即当卫星运行到S^′时，目标并不能恰好处于T^′。T^″、$ T'''$ 分别为考虑了卫星的姿态机动角后，目标可以被卫星观测到的边界点。当卫星运行到成像点S^′点时，如果目标位于边界点T^″、$ T'''$之间，卫星能够通过调整姿态使目标进入相机视场，则认为目标可见；否则目标不可见，需要对卫星进行变轨，使卫星到达成像点S^′的时间提前或者推迟，从而满足成像条件。

首先考虑二体模型情况。根据所建立的坐标系，容易求得轨道倾角i满足如下关系：

$ \tan i = \frac{z}{y}. $

(1)

由于地心、X轴与赤道交点、升交点均在轨道平面上，代入轨道平面x+By+Cz+D=0中，有如下轨道平面方程：

$ y\tan i - z = 0. $

(2)

目标在该坐标系中的坐标为 $ T\left(R_{\mathrm{e}} \cos \varphi \cos \theta,\right.$$ \left.R_{\mathrm{e}} \cos \varphi \sin \theta, R_{\mathrm{e}} \sin \varphi\right)$ ，其中：θ为目标和轨道升交点的经度之差；φ为目标的纬度；R_e为地球半径。则目标所在的纬线圈的圆方程为

$ {x^2} + {y^2} = {({R_{\rm{e}}}\cos \varphi )^2}. $

(3)

该圆上，z值恒定不变，即

$ z = {R_{\rm{e}}}\sin \varphi . $

(4)

将式(3)、(4)联立求解，得到成像点的星下点T^′的坐标。进而可以求得目标通过地球自转到达该点的时间为

$ {t_{{T^\prime }}} = {\theta _{T{T^\prime }}}/{\omega _{\rm{e}}}. $

(5)

其中：θ_TT^′为目标由T运动到T^′绕Z轴转过的角度；ω_e为地球自转角速度，其值约为7.292×10^－5 rad/s。

接着，在二体模型基础上考虑地球J₂摄动的影响。由于J₂摄动，卫星轨道根数中的升交点赤经Ω、近地点幅角ω随时间长期变化，其平均变化率分别为

$ {\dot \varOmega = - \left[ {\frac{3}{2}\frac{{\sqrt \mu {J_2}R_{\rm{e}}^2}}{{{{\left( {1 - {e^2}} \right)}^2}{a^{\frac{7}{2}}}}}} \right]\cos i,} $

(6)

$ {\dot \omega = - \left[ {\frac{3}{2}\frac{{\sqrt \mu {J_2}R_{\rm{e}}^2}}{{{{\left( {1 - {e^2}} \right)}^2}{a^{\frac{7}{2}}}}}} \right]\left( {\frac{5}{2}{{\sin }^2}i - 2} \right).} $

(7)

其中：a为卫星轨道半长轴；e为卫星轨道偏心率；μ为地心引力常数，其值为3.986 004 418×10¹⁴ m³/s²。

注意到太阳同步轨道倾角i大于90°，则 $\mathit{ \dot{\Omega}}$ 为正，即逆行轨道的交线将向东移动。据此对t_T^′进行如下修正：

$ {t_{{T^\prime }}} = {\theta _{T{T^\prime }}}/\left( {{\omega _{\rm{e}}} - \dot \varOmega } \right). $

(8)

进一步地，通过卫星的姿态机动角，可以求得当卫星处于成像点S^′时，目标可见的边界点T^″、$T''' $ 。下面以T^″为例，如图 2所示。

卫星到达成像点S^′的时间为

$ {t_{S'}} = {t_1} + n{T_1}. $

(9)

其中：n为自然数，代表卫星转过的圈数；t₁为卫星从S第一次运动到成像点S^′的时间；T₁为卫星运行轨道周期。

在考虑J₂摄动的情况下，加入卫星近地点幅角平均变化率 $ \dot{\omega}$ ，T₁计算表达式为

$ {T_1} = \frac{{2{\rm{ \mathsf{ π} }}{a^{\frac{3}{2}}}}}{{\sqrt \mu }} - \dot \omega . $

(10)

卫星和目标的时间差为

$ \Delta t = {t_{{S^\prime }}} - {t_{{T^\prime }}}. $

(11)

那么，判断目标是否可见的条件为

$ |\Delta t| \le {t_{{T^\prime }{T^{\prime \prime }}}}. $

(12)

目标的可见性与Δt相关，具体分为4种情况，如表 1所示。

2.2 变轨策略

从表 1可以看出，若|Δt|>t_T^′T^″，目标不具备可见条件，需要卫星通过相位机动，提前或推迟到达成像点S^′。若Δt<0且|Δt|>t_T^′T^″，则需要向后相位机动；若Δt>0且|Δt|>t_T^′T^″，则需要向前相位机动。为实现卫星对目标的快速响应，卫星通过相位机动，仅在过渡轨道上运行一圈后即再次回到原轨道，最终到达成像点S^′。以向后相位机动为例，卫星的相位机动方式如图 3所示。

图中轨道1为原轨道；轨道2为相位机动的过渡轨道；T₂为轨道2的周期；Δv₁₂为由轨道1机动到轨道2施加的脉冲；Δv₂₁为由轨道2返回轨道1施加的脉冲。

卫星需要调整的时间差为

$ \delta = {t_{S'}} - {t_{{T^{\prime \prime }}}}. $

(13)

为使卫星通过变轨后具备成像条件，变轨时间最迟不应超过

$ {t_{{\rm{man}}}} = {t_{{T^\prime }}} - {t_{{T^\prime }{T^{\prime \prime }}}} - {T_2}. $

(14)

接下来考虑构型返回时的燃料限制问题。构型的返回需要进行与之前相反的相位机动，返回构型如图 4所示。

图中轨道3为返回构型时的过渡轨道；T₃为轨道3的周期；Δv₁₃为由轨道1机动到轨道3施加的脉冲；Δv₃₁为由轨道3返回轨道1施加的脉冲。

为了满足最大返回时间T_return，并使燃料消耗尽可能小，应该使卫星在轨道3上运行圈次N尽可能多，则N应满足如下条件：

$ {N{T_3} \le {T_{{\rm{return}}}}.} $

(15)

需要调整的时间差为δ，则又有如下条件：

$ {{T_3} = {T_1} - \frac{\delta }{N}.} $

(16)

N为自然数，则有最大圈数为

$ {N_{\max }} = \left\lfloor {\frac{{{T_{{\rm{return}}}} + \delta }}{{{T_1}}}} \right\rfloor . $

(17)

根据相位机动的相关理论，可以得到相位机动、返回构型过程中的相关变轨参数。为满足任务约束，卫星在进行相位机动时，应当保留能够提供返回构型所需速度增量的燃料。

2.3 数值方法

本方法在解析法的基础上进行精确计算，通过在卫星运动到成像点时刻t_S^′附近搜索目标计算可见性，可以精确获得目标的可见时段。

该方法建立在J₂摄动模型数值积分的轨道动力学模型基础上。在J2000惯性系中，考虑J₂摄动的动力学方程为

$ \left\{ {\begin{array}{*{20}{l}} {\dot x = {v_x},}\\ {\dot y = {v_y},}\\ {\dot z = {v_z}.} \end{array}} \right. $

(18)

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot v}_x} = - \frac{{\mu x}}{{{r^3}}}\left( {1 + \frac{3}{2}{J_2}{{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^2}\left( {1 - 5\frac{{{z^2}}}{{{r^2}}}} \right)} \right),}\\ {{{\dot v}_y} = - \frac{{\mu y}}{{{r^3}}}\left( {1 + \frac{3}{2}{J_2}{{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^2}\left( {1 - 5\frac{{{z^2}}}{{{r^2}}}} \right)} \right),}\\ {{{\dot v}_z} = - \frac{{\mu z}}{{{r^3}}}\left( {1 + \frac{3}{2}{J_2}{{\left( {\frac{{{R_{\rm{e}}}}}{r}} \right)}^2}\left( {3 - 5\frac{{{z^2}}}{{{r^2}}}} \right)} \right).} \end{array}} \right. $

(19)

其中：x、y、z分别为卫星在J2000惯性系下的位置分量；v_x、v_y、v_z分别为卫星在J2000惯性系下的速度分量；r为卫星与地心的距离。

根据卫星的初始位置和速度，利用数值积分，外推得到卫星任意时刻的位置和速度。在此基础上，结合解析法给出的卫星运动到成像点时刻t_s^′，在其附近分析目标的可见时段。

在可见性分析过程中，应当满足太阳高度角、最大地心角和卫星姿态机动角的限制。光学成像卫星只有在目标处于光照区时才能对目标进行成像，则需要计算太阳高度角。当太阳高度角达到一定角度时，目标可见；否则不可见。太阳高度角的计算公式^[20]为

$ \sin H = \sin \varphi \sin \delta + \cos \varphi \cos \delta \cos t. $

(20)

其中：H为太阳高度角；φ为当地的地理纬度；δ为当日的太阳赤纬，单位为度(°)；t为当时的太阳时角，单位为度(°)。

太阳赤纬计算公式为

$ \begin{array}{l} \begin{array}{*{20}{l}} {\delta = 0.006{\kern 1pt} {\kern 1pt} 918 - 0.399{\kern 1pt} {\kern 1pt} 912\cos b + }\\ {0.070{\kern 1pt} {\kern 1pt} 257\sin b - 0.006{\kern 1pt} {\kern 1pt} 758\cos 2b + } \end{array}\\ \begin{array}{*{20}{c}} {0.000{\kern 1pt} {\kern 1pt} 907\sin 2b - 0.002{\kern 1pt} {\kern 1pt} 697\cos 3b + }\\ {\quad 0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 48\sin 3b,} \end{array} \end{array} $

(21)

$ {b = 2{\rm{ \mathsf{ π} }}(N - 1)/365.} $

(22)

其中：N为每年从1月1日起，距计算日的天数；b为中间变量。

太阳时角计算公式为

$ {t = 15({t_s} - 12),} $

(23)

$ {{t_s} = h + m/60 + {\rm{lon}}/15.} $

(24)

其中：t_s为真太阳时；h为当地时间小时数；m为当地时间分钟数；lon为当地经度。

如图 5所示，卫星距地面的距离为H，卫星可视最大范围所对应的地心角度θ_max为

$ {\theta _{{\rm{max}}}} = {\cos ^{ - 1}}\left( {\frac{{{R_{\rm{e}}}}}{{H + {R_{\rm{e}}}}}} \right). $

(25)

假设给出了地面点T在地固系中的位置R，卫星于某时刻t相对地固系的位置r(t)和速度v(t)，若r(t)与R的夹角小于θ_max，则满足成像条件，目标可见；否则不可见。

v_ST为卫星指向目标的向量，v_SO为卫星指向地心的向量，卫星观测角为其夹角，即

$ {\alpha _{{\rm{obs}}}} = {\cos ^{ - 1}}\left( {\frac{{{\mathit{\boldsymbol{v}}_{ST}} \cdot {\mathit{\boldsymbol{v}}_{SO}}}}{{\left| {{\mathit{\boldsymbol{v}}_{ST}}} \right|\left| {{\mathit{\boldsymbol{v}}_{SO}}} \right|}}} \right). $

(26)

若目标处于卫星姿态机动范围内，即α_obs小于最大卫星姿态机动角，则目标可见；否则不可见。

在成像任务中，对于需要进行轨道机动的情况，根据解析方法给出的变轨策略完成变轨。变轨结束后，按照数值方法对目标可见时段进行分析。

3 多目标任务规划方法 3.1 整体可见性分析

以单星单目标观测方法为基础，分析卫星在不变轨情况下对全球随机分布的多个点目标的可见性，形成可见性分析表，如表 2所示。

3.2 变轨策略分析

在可见性分析表基础上，分别对每颗卫星进行变轨策略分析。对于每个不可见的目标，卫星通过变轨使其可见，若变轨及返回构型消耗的燃料满足最大燃料限制，则将该变轨策略加入策略集合。通过该方法，得到该卫星的变轨策略集合。

将所有卫星的变轨策略进行组合，即可得到任务规划的总策略集。当任务规模达到一定程度时，总策略集将会异常庞大，普通的搜索方法已经不再适用。

3.3 遗传算法搜索最终策略

为了解决上述问题，采用遗传算法进行任务规划。遗传算法(genetic algorithm, GA)^[21]是一种模拟自然界中自然选择、遗传、进化的算法，具有较快的收敛速度和较强的鲁棒性。本文采用的遗传算法具体思路如图 6所示。

为使卫星通过变轨覆盖更多目标，同时使目标观测次数分布更加均匀，需要对遗传算法中的适应性计算方法进行规定：若在总策略集中选定一个变轨策略，根据变轨后的可见性分析表得到每个目标的被观测次数；对每个目标而言，若其被观测次数为λ，则其对应的适应性为

$ {\rm{adaptability}} = \sum\limits_{k = 1}^\lambda 0 {.5^{k - 1}}. $

(27)

其中：当λ为0时，代表该目标未被观测，其适应性为0。

最后将所有目标对应的适应性求和，得到该策略对应的适应性。适应性计算流程如图 7所示。

遗传算法中的自然选择采用锦标赛选择方法；交叉遗传方法为PMX(partial-mapped crossover)；变异方法为随机变异。经过多代自然选择，种群逐渐显现出高适应性的优良特性，最终产生任务所需的个体即最终变轨策略。星座中各卫星根据此变轨策略进行变轨，从而完成对多目标的成像任务。

4 算例分析

本文分别对单星单目标观测方法和多目标任务规划方法选取算例进行仿真。仿真起始时刻为2018-11-07 04∶00∶00 (UTCG)，仿真结束时刻为2018-11-08 04∶00∶00 (UTCG)，仿真时间为24 h。仿真环境为3.6 GHz CPU的计算机，编程为MATLAB语言。卫星最大姿态机动角为30°，返回构型最大时长限制为12 h。数值方法积分输出步长设置为1 s，如有更高精度要求可以将步长设置为更短。

4.1 单星单目标可见性分析

将单星单目标观测方法与STK仿真方法(HPOP高精度模型)，在相同轨道和目标参数条件下进行对比，分析单星单目标观测方法的可行性和精度。选取一颗卫星和两个目标，参数如表 3和表 4所示。

利用本方法和STK方法，分析该卫星对以上2个目标的观测情况，结果对比如表 5所示。

通过表 5可以发现，本方法与STK方法得到的可见性判断一致，证明解析法有效可行。本方法开始和结束时刻均与STK方法接近，可见时长误差约为5.82%。同时，解析、数值方法的结合，由解析法给出目标大致可见区间供数值法搜索，与纯数值方法盲目搜索可见性相比，极大缩小了搜索范围，提高了目标可见性的分析效率。

针对表 5中目标2不可见的情况，利用本方法给出变轨策略，得到卫星进行相位机动时刻为13 359.3 s，总速度增量为234.45 m/s。利用STK方法对该变轨策略进行可行性验证，卫星实现了对目标2的观测，可见时长为31.5 s，证明该变轨策略确实可行。

4.2 多目标任务规划

星座共9颗卫星，均匀分布在3个太阳同步圆轨道上。轨道半长轴均为6 878 137.0 m，偏心率为0，倾角为97.406 5°，近地点幅角为0，其余轨道根数如表 6所示。目标集由随机生成的100个目标组成。

在不变轨的情况下，分析星座中各卫星分别对所有目标的观测情况，统计结果如表 7所示。在此基础上分析总策略数达到46 904 014 080个。使用遗传算法对总策略集进行搜索，初始种群规模为100个，遗传代数为500。经过200代自然选择的种群适应性变化情况如图 8所示。

通过观察图 8发现，随着种群不断遗传迭代，种群平均适应性和最高适应性均呈现显著提高；最高适应性在100代内即出现，其对应的个体即为任务规划的最终策略。将该策略应用到星座对目标的成像任务中，最终观测结果如表 7所示。

由表 7可知，进行任务规划后，单个目标最少观测次数提升至2次，完全覆盖100%的目标，提高了星座观测能力；观测次数提升约7.45%，提高了卫星的利用率；观测次数标准差降低，适应性提高，分布更加均匀。

在计算规模上，遗传算法每代只需要对100个个体(策略)进行操作，即使进行500代遗传，相对总策略数46 904 014 080而言规模极小；与此同时，由于自然选择率为0.5，本代种群的50%都会保留遗传到下一代，从而更加减小了计算规模。

5 结论

本文提出了一种机动卫星星座对多目标成像任务规划算法。利用太阳同步圆轨道性质提出解析法，对目标可见性进行快速估计，计算速度快，可见性判断准确；对于不可见目标，提出了通过相位机动使其可见的变轨策略，可行性高；数值法通过数值积分获得目标的准确可见时长，验证了解析法，且搜索时间短。对多目标进行成像任务规划时，在单星单目标观测方法的基础上，分析了多目标可见性；针对不可见目标提出变轨策略，形成变轨策略总集合；利用遗传算法，以适应性最高为优化指标，进行多代自然选择，最终获得星座的最优变轨策略，大大提高了计算效率。数值结果表明，规划后平均目标观测次数提高约7.45%，星座完全覆盖100%的目标。