应急任务响应时间最优的多星成像规划方法

航天工程大学，北京 101400

成像卫星是一类利用星载传感器获取特定区域图像信息的对地观测卫星，主要包括可见光、红外、多光谱、合成孔径雷达(Synthetic Aperture Radar, SAR)等成像卫星[1]，光学成像卫星因其高分辨率和信息的丰富性仍然是成像卫星的主流。目前，中国已初步建立天基对地观测系统，用户需求复杂、卫星资源统一调度困难等问题进一步凸显。当前自然灾害、工业事故等应急事件频发，用户的需求不再仅仅局限于是否能够完成成像任务，而对响应时间也提出了要求。卫星管控部门在进行任务规划时既需要考虑完成尽可能多的任务，又要考虑部分任务的优先性和响应时间，如何合理地调度卫星资源成为多星成像规划领域亟待解决的问题。

多星成像规划是NP-Hard问题，通常要在简化约束条件的基础上构建模型，国内外在建立模型和求解方面开展了大量研究。Globus等[2-3]建立了多星任务规划的约束满足问题(Constraint Satisfaction Problem，CSP)模型，考虑了卫星的多设备约束条件，并用遗传算法得到了较好的结果；Gabrel等[4]建立了图论模型，采用有向无环图表示卫星在观测任务之间的转换约束，使用最短路径算法进行求解；Bianchessi等[5]建立了混合整数规划模型，验证了禁忌搜索算法求解该模型的可能性。国内研究者对以上的经典规划模型持续改进并应用于不同背景。王智勇等[6]建立了考虑卫星载荷连续侧摆的多星联合任务规划模型；林晓辉等[7]针对敏捷光学成像卫星密集区域推扫成像任务规划问题进行了建模和求解；赵萍等[8]设计了一种改进遗传算法，验证了此算法应用在卫星自主任务调度问题上的有效性，但以上研究并没有考虑应急条件。随着实际中应急需求的突出，一些研究者针对应急任务规划也展开了研究。Verfaillie等[9]研究了新到达任务在不影响原方案基础上的插入方法；王建江等[10]提出了综合考虑任务合成、修复和向后移位的多星动态应急调度算法，但这些研究集中在应急任务对原调度方案的扰动上。杨正磊等[11]首次以应急需求的响应时间为目标进行综合规划研究，但局限于单个应急任务的规划。

目前的研究多是单独优先规划应急任务，虽然保证了应急任务的响应时间但牺牲了常规任务的完成度，使整体调度收益受到较大影响。针对此问题，本文建立综合考虑应急任务响应时间和任务总收益的多星成像规划模型，将规划问题分解为任务时间窗选择和单轨动态规划两部分，分别设计自适应免疫算法(Adaptive Immune Algorithm, AIA)和前向动态规划算法(Forward Dynamic Programming Algorithm, FDPA)进行优化求解，并进行了仿真校验。

1 多星成像规划模型

本文研究的多星成像规划问题可描述为：在某一规划时间段内，一批成像任务T在卫星资源集合S下具有多个时间窗。在满足卫星成像一系列约束的前提下为每个成像任务选择合适的时间窗，并确定每颗卫星的观测路径，使所有应急成像任务尽可能被完成和完成时间尽可能早，在此基础上再完成尽可能多的常规成像任务。依据问题特征，建立该问题的约束满足模型，从模型要素的定义、目标函数和约束条件3方面描述模型。

1.1 模型要素的定义

模型要素定义如表1所示。

(1)

式中：i为任务编号；j为卫星编号；k为卫星的轨道圈次编号；b和e分别为成像开始时间和结束时间；a为侧摆角。

1.2 目标函数

首先给出模型的决策变量xijk：

(2)

模型具有两级优化目标，一级目标是应急任务完成的平均响应时间最短，即下式中的f1最小。

(3)

二级目标是任务总收益最大，即下式中f2最大。

(4)

1.3 约束条件

1)本文假设成像任务均为单次执行即可完成，所以任务执行时最多占有一个卫星资源：

(5)

2)任务所选取卫星的最低分辨率要优于任务要求的分辨率：

(6)

3)任务的成像时间需在要求的执行时间范围内：

(7)

4)任一合成任务的侧摆角绝对值不能大于所用卫星的最大侧摆角度：

(8)

5)任一合成任务的成像持续时间不能大于所用卫星的遥感器单次最长开机时间：

(9)

6)同轨道圈次两个连续合成任务间的姿态转换时间要小于两合成任务的间隔时间：

(10)

7)卫星在任一轨道圈次上成像姿态机动次数不大于单轨任务数量和卫星单轨最大侧摆次数的较小值：

(11)

多星成像规划问题可以分解为选择任务的时间窗和确定卫星单轨道圈次的最优观测路径两个子问题，因此本文设计了时间窗选择和单轨动态规划两个模块耦合求解，并在两个模块分别使用自适应免疫算法和前向动态规划算法进行优化。单轨动态规划的结果反馈至时间窗选择模块计算目标函数，根据目标函数的性能不断调整部分任务的时间窗选择，迭代计算至最优解，模型求解流程如图1所示。

图1 模型求解流程Fig.1 Solving process of the model

2 基于自适应免疫算法的时间窗选择

本文依据免疫算法的基本原理，在亲和力函数中加入判断解的收益是否达到期望收益的自适应阈值λ，用抗体的浓度体现解的多样性，用期望繁殖率来寻找亲和力高的且浓度低的抗体，在抗体种群更新的同时选取部分亲和力较高的和部分期望繁殖概率较高的抗体进入记忆细胞(对优良解的保留)和下一代抗体种群。每次更新种群后根据记忆细胞中优良解的表现，逐步调整自适应阈值λ，迭代计算至一定次数，将当前记忆细胞中的最优解输出。

2.1 设计亲和力函数

(12)

2.2 抗体相似度和浓度

在浓度计算之前，需要计算抗体种群中抗体的相似度。抗体长度均为lz，抗体每一个编码位置对应一个任务的时间窗选择情况，若两抗体有Ng个任务选择了同一时间窗，则两抗体的相似度X为：

(13)

当X达到设定的相似度阈值υ时，则认为此两抗体为相似抗体。抗体浓度计算公式为：

(14)

式中：Cz为抗体浓度；M为抗体总数；NX为抗体种群中与抗体z相似的抗体数。

2.3 抗体期望繁殖概率

期望繁殖概率Pz由抗体的亲和力和浓度共同决定，其计算公式如下：

(15)

式中：α∈(0,1)为设定的多样性评价参数。由式(15)可见，Pz在一定程度上既鼓励高亲和力的抗体，又抑制高浓度的抗体，能够保持抗体种群的多样性，避免陷入局部最优解。

2.4 抗体种群、记忆细胞的产生和更新

初始抗体种群由随机贪婪算法[15]产生。记忆细胞实质是优良解的记忆库，设记忆细胞的容量为B。更新种群时，为防止亲和力高的抗体因为浓度高而被淘汰，先选取抗体种群中B/2个Az最大的抗体进入记忆细胞和下一代抗体，再在剩下抗体种群中分别选取B/2个和(M-3B)/2个Pz最大的抗体进入记忆细胞和下一代抗体，即淘汰B个Pz最小的抗体。

保留的抗体通过交叉和变异产生新抗体，再加入记忆细胞中的上一代优良抗体更新种群。考虑到应急任务和常规任务的不同，为更好地产生新个体，本文设计了两段交叉和变异方法，如图2所示，虚线左边为应急任务段，右边为常规任务段，抗体的两段交叉和变异是同时进行的。

图2 抗体交叉与变异Fig.2 Cross and mutation of antibodies

每次更新记忆细胞后，判断其中B/2的抗体收益IB/2是否均达到期望收益λIa，若均达到，说明λ的取值偏小，令λ增加1%，期望收益将在下一轮迭代中提高。自适应免疫算法实现流程如图3所示。

图3 自适应免疫算法实现流程Fig.3 Flowchat of adaptive immune algorithm

3 单轨动态规划

经过每轮的时间窗选择，应急任务的响应时间已经确定，每颗卫星的不同轨道圈次也已形成元任务序列，单轨规划的目的是在最大限度完成应急任务的同时使任务总收益最大。考虑到单轨规划结果要返回至时间窗选择模块，遍历寻找最优解耗时过高，且不利于求解，本文借鉴文献[16]的动态规划思想，对单轨任务动态合成方法进行改进，提出了一种面向应急任务的前向动态规划算法(FDPA)，将任务合成分为多个前后关联的阶段，按照顺序解法递推，直到达到卫星单轨最大侧摆次数或者规划至单轨最后一个任务。

3.1 面向应急任务的最优合成算法

设某卫星单轨道圈次上的元任务按照成像开始时间顺序排列后为{m1,m2,…,ml}，l为单轨任务数量。卫星每观测一次形成一个合成任务c，必包含成像时间最早的元任务mj和最晚的元任务mi，则定义cj,i为以j为起点和以i为终点的合成任务，Aj,i={aj,aj+1…,ai}为cj,i包含的元任务侧摆角集合。cj,i观测角度的浮动范围和覆盖范围分别如式(16)、式(17)所示：

[maxAj,i-θj/2,minAj,i+θj/2]

(16)

[maxAj,i-θj,minAj,i+θj]

(17)

cj,i对应一次观测活动，即决策过程中的一个阶段，只需要确定cj,i的最优观测角度，使其覆盖最多的应急任务和获得最大收益，即可认为此合成任务为此阶段以j为起点、i为终点的最优合成任务，设其最优合成收益为rj,i，最优观测角度为aj,i。面向应急任务的最优合成算法步骤为：

步骤1：判断mj、mi是否满足任务合成约束条件且i≥j，若不满足，令rj,i=0、aj,i=0，合成终止；若满足，转入步骤2。

步骤4：返还cm对应的am和rm，令aj,i=am，rj,i=rm，合成结束。

3.2 前向动态规划算法

定义动态规划第k阶段cj,i的前接最优合成任务集为：第k-1阶段终点在j之前，且与cj,i满足姿态转换时间约束的最优合成任务集合。第k阶段以i为终点的最优合成任务Ck(i)由{cj,i}(k≤j≤i)及其前接最优合成任务集共同确定。

设fk(i)表示Ck(i)的收益，sk(i)表示其起点。设Fk为第k阶段不同终节点的最优合成任务收益向量，Sk为对应的起点向量，并有：

(18)

则第k+1阶段以i为终点的最优合成任务的收益Fk+1(i)计算公式为：

(k+1≤j≤i)

(19)

图4 前接最优合成任务不存在的情况Fig.4 Non-existent situation of attachable optimal composite task

若第k+1阶段cj,i的前接最优合成任务不存在，遍历第k阶段终点在j之前的所有满足姿态转换时间约束的合成任务cx,y(k≤x≤y≤j-1)，保留在包含最多应急任务基础上收益最大的合成任务，用其收益max {rx,y}代替Fk(m)代入式(19)计算。

前向动态规划算法(FDPA)步骤如下：

步骤1：将卫星sj不同轨道圈次的元任务按照成像开始时间进行排序，形成每个轨道圈次元任务序列{m1,m2,…,ml}，每个轨道圈次任务序列均转入步骤2进行并行计算。

(20)

(21)

步骤3：计算动态规划第1阶段的收益向量F1和对应的起点向量S1，令k=1，转入步骤4。

4 解的修复

若调度方案确定后出现应急任务没有全部完成的情况，应进行适当修复以尝试完成更多的应急任务，为避免计算复杂度进一步升级，本文设计了较为简洁高效的修复方法。

步骤1:将剩下未安排的应急任务按照优先级从高到低排序，选择当前优先级最高的应急任务ts，转入步骤2。

步骤2:遍历ts的每个时间窗，计算是否有其他应急任务冲突，若均有应急任务冲突，则ts安排失败，用规划总时长代替其响应时间，结束修复；否则转入步骤3。

步骤3:对于只有常规任务冲突的时间窗，计算ts插入需删除的常规任务优先级之和并定义其为冲突损失。将ts插入至冲突损失最小的时间窗，若有多个时间窗冲突损失同时最小，选择开始时间最早的时间窗插入ts，结束修复。

5 仿真校验

本文通过大量仿真试验测试AIA+FDPA的性能，并与启发式算法和智能优化算法进行了比较。启发式算法选择文献[17]的动态合成启发式(DTMH)算法和文献[10]的综合考虑任务合成、修复和向后移位的多星动态应急调度(TMRBS-DES)算法的结合形式，常规任务使用DTMH算法，应急任务使用TMRBS-DES算法。智能优化算法选择文献[18]中的混合遗传禁忌搜索(HGTS)算法，HGTS算法先使用遗传算法在解空间进行全局搜索，迭代至一定次数获得多个局部最优解，再使用禁忌搜索算法在局部最优解邻域进行小规模搜索。

5.1 测试方法与参数

本文参考国内外民用光学遥感卫星，选取了8颗光学遥感卫星平台，其轨道根数如表2所示，平台能力的设定如表3所示。

表2 试验选取的卫星轨道六根数Table 2 Orbital elements of satellites selected in experiment

表3 试验选取的卫星平台能力Table 3 Capability of satellites selected in experiment

为加强选取任务的真实性，本文从中国地震台网获取了一年内全球和我国四川省发生地震的地理位置和震级信息，选取全球震级大于4级的100/200/300/400个地点和四川震级大于3级的25/50/75/100个地点作为观测目标，如图5和图6所示，分别从全球背景和特定区域背景下进行测试。

设某测试实例下任务集T对应的震级集合为ET，任务ti的震级为Ei，则本文定义ti的优先级wi计算公式如下：

(22)

本文设定规划起始时刻(UTC)为2019年3月24日00:00:00，结束时刻(UTC)为2019年3月25日00:00:00。假设用户提出需求的时刻均为规划起始时刻，目标可见时间窗采用STK软件计算，并设置了卫星拍摄地面目标时的太阳高度角不小于10°时为有效拍摄。所有算法均采用Matlab实现，统一在配置为Intel(R) Xeon(R) CPU E5-2630的计算机上运行，其中两种智能优化算法的相关参数如表4和表5所示。

5.2 任务规模对算法效果的影响

本节主要考察任务总数量对算法的影响，应急任务数量设定为任务总数的10%(向上取整)。测试结果如图7所示。

图5 全球观测目标分布情况Fig.5 Distribution of global ground targets

图6 特定区域(四川)观测目标分布情况Fig.6 Distribution of ground targets in specific area (Sichuan)

表4 自适应免疫算法参数Table 4 Parameters of AIA

表5 混合遗传禁忌搜索算法参数Table 5 Parameters of HGTS

从图7可以看出，随着任务总数量的增加，所有算法得出的应急任务平均响应时间和任务总收益都随之增大，应急任务均完成，故没有比较展示。当任务总数较小时，任务冲突较少，3种算法的结果差距很小。随着任务总数增加，DTMH+TMRBS-DES的解的质量下降较快，且应急任务平均响应时间指标下降程度较任务总收益更明显，这是由于其没有考虑已插入任务的退出规则，造成后续任务的时间窗后移甚至无法插入。HGTS算法解中的任务总收益指标略优于AIA+FDPA，但应急任务平均响应时间指标与AIA+FDPA相比较差，这是因为HGTS设计的目标函数中任务总收益占比较大，且在染色体交叉和变异时并没有针对应急任务进行处理。

图7 任务规模对算法效果的影响Fig.7 Influence of the number of tasks on algorithm performance

从目标分布来看，特定区域背景下各算法的表现都优于全球背景，一方面是任务数量较少减少了冲突，且目标的密集分布使任务合成的优势得到体现。 DTMH+TMRBS-DES算法中的任务合成启发式规则表现较好，其解的质量与两种智能优化算法得到的最优解差距变小。

从各算法的耗时来看，DTMH+TMRBS-DES作为启发式算法具有巨大的优势，在所有测试实例中耗时都小于1 s，而AIA+FDPA和HGTS算法耗时是其数百倍，并随着任务总数的增加有呈指数增长的趋势。且这两种智能优化算法在面对特定区域背景下目标密集分布时，计算耗时明显大于全球背景下目标较稀疏分布的情况，这是因为目标密集分布时，分配到部分卫星资源的任务变多，使单轨动态规划和约束检查计算耗时都变长。

5.3 应急任务比例对算法效果的影响

5.2节的测试中，由于应急任务数量较少，所有测试结果中应急任务均完成，本节分别在全球400个目标和特定区域100个目标背景下，使应急任务比例从10%到40%递增，测试应急任务数量占总任务数量的比例对各算法的影响，结果如图8所示。

图8 应急任务比例对算法效果的影响Fig.8 Influence of percentage of emergency tasks on algorithm performance

从图8可以看出，应急任务比例的增加对应急任务平均响应时间、未完成应急任务数具有明显的影响。DTMH+TMRBS-DES算法表现最差，这是因为其确定性的规则加剧了应急任务之间的冲突。而AIA+FDPA算法在保证应急任务平均响应时间最优的同时在应急任务完成率方面同样表现较好，这和此算法具备解的修复过程且在应急任务未完成时具有惩罚机制有关。

此外，应急任务平均响应时间随应急任务比例增大而增加，这是由于应急任务数量的增多导致部分应急任务的执行时间不断后移，且一旦错过了早期的时间窗将可能经历长期的夜晚不可观测时段。

任务总收益受应急任务比例的影响较小，即使应急任务比例从10%增长至40%，各算法得到的任务总收益下降仅在1.4%～3.6%之间。这是由于总任务数量是恒定的，只是部分任务从常规任务变为应急任务，因为追求响应时间而导致此类任务执行时间前移，影响了少数常规任务的执行。

6 结束语

本文为在成像卫星任务规划中可以综合考虑应急任务和常规任务，建立两级目标优化模型，提出自适应免疫算法(AIA)和前向动态规划算法(FDPA)相结合的分解优化策略进行求解，并通过大量仿真试验对AIA+FDPA算法进行了测试，并与DTMH+TMRBS-DES算法和HGTS算法进行了对比。仿真结果表明，本文的AIA+FDPA算法得到的应急成像任务的响应时间优于其他算法，获得的任务总收益和HGTS算法接近，且耗时小于HGTS算法，这为需同时考虑应急任务和常规任务的大规模多星成像规划问题建模和求解提供了很好的途径。同时，AIA+FDPA算法的计算耗时仍然较长，不利于对应急任务的规划，后续工作可在寻优策略和算法效率方面进一步开展。