遥感卫星对区域目标的可见窗口是指其搭载的视场负载能够观测到区域目标的时间范围。该可见窗口计算是卫星观测任务的基础，是区域目标覆盖计算、成像任务规划等多种任务的前提。Voigt等^[1]认为，随着在轨卫星增多及卫星相关技术发展，合理调配卫星资源将成为新研究趋势。多星多目标调度问题逐渐成为研究的难点和热点问题^[2-6]。由于调度问题的解空间巨大，而可见窗口计算是其必要的求解前提。若能进一步提高可见窗口的计算效率，将显著缩短观测任务规划等后续任务计算耗时，具有重要的研究价值与工程意义。

可见性问题的研究最早聚焦在卫星对点目标的可见性，经典的算法是传播算法^[7]，即通过固定步长遍历离散时间点依次判断可见性，该算法计算准确，但计算效率低，常在文献中作为对比算法。后续研究目标均着力在减少计算耗时。文[7]对卫星-卫星、卫星-点目标的可见性第一次提出使用视函数方法解析进行预测，大大提升了计算效率，但只适用于小偏心率轨道；文[8]改进了视函数的适应性，使其可以推广至任意轨道；文[9]使用大圆弧近似一个轨道周期下的星下点轨迹，其模型简单、计算效率高，但只适用于低轨道卫星；文[10]提出用粗搜索-精细估计的两阶段方法进行可见窗口的快速求解，显著缩短了计算用时。后续一部分研究^[11-13]是基于以上方法的改进。

上述研究均为卫星对点目标可见窗口的计算方法，而对于区域目标可见性，现有研究方法可分为3类：1) 轨道周期过滤，通过卫星轨道及区域目标信息筛除大部分的不可见时段，如文[14]针对圆锥与矩形视场，用大圆近似星下点轨迹，降低了计算用时，但只适用于低轨卫星；2) 通过对区域内部或区域边界减小离散采样点，如文[15]考虑圆锥视场，将问题转化为区域每条边界的时间窗口计算可见窗口的并集；3) 考虑使用并行、两阶段搜索、动态步长等计算手段减小计算量，如文[16]在计算过程中考虑了大数据并行计算，文[17-18]在计算过程中考虑了动态步长的情况，文[19]提出用二分法进行精确搜索以增大初始搜索步长，同时只对区域目标的边界进行离散采样以减小计算规模。

上述区域目标可见窗口计算方法，仿真1 d的耗时最快约为0.1 s量级^{[14-17, 19]}(相比于传播算法近1 000倍加速比)。在面对如大规模卫星调度问题时，由于可见窗口被反复调用，现有求解效率仍需提高。

本文针对这一问题，将区域边界合理近似，进而使用空间矢量提出了一种基于卫星圆锥、矩形视场的区域目标可见窗口半解析快速计算方法。在保证计算精度仍满足工程目标的基础下，该方法大幅提高原有计算效率。主要步骤如下：首先，利用空间几何关系，将所有区域的边界描述为大圆弧或普通圆弧；其次，根据该问题数学模型得到快速判断可见性解析方法，与原有算法相比，该方法不再对区域目标进行离散采样计算；再利用判断算法解析地判断卫星在每个离散时间步长下对区域目标的可见性；最后，引入二分法精确搜索可见起始结束时刻，并通过Lagrange插值得到卫星的位置和速度信息，以最终求解区域目标的可见窗口。

1 区域目标可见性问题建模

卫星在运行过程中，其搭载的视场传感器将对地面特定目标进行间断的观测。本文所研究的区域目标就是该特定目标的一类，具体是指覆盖地球表面的封闭凸多边形。若卫星视场能覆盖区域目标任意部分时，就视为卫星对该区域目标可见。在具体求解之前，为明确求解过程中的数学问题以及相应的问题假设，同时为了便于描述后续提出的解析算法，本文将在此给出区域目标可见性的几何判定条件。

1.1 卫星对区域目标可见性的数学描述及假设

本文在计算卫星对区域目标可见性的过程中，假设地球为球体，忽略地球表面海拔和地球扁率的影响。本文算例表明结果误差已满足工程要求。若考虑精确模型，可在本文方法的基础上以可见窗口的起止点作为初值进行精确迭代求解。

在求解区域目标的可见窗口时，程序需要判断：在某一时刻下，卫星对区域目标是否可见。假定已经知道卫星在该时刻下的位置、速度以及搭载的传感器类型，则卫星在该时刻视场投影到地球表面的区域是确定的。设该区域为A，同时设地球表面的区域目标为B。那么上述实际问题的数学表述即为：判断球面上的区域A与区域B是否有重叠部分。

1.2 判断区域目标可见性的几何判定条件

区域可以通过离散点采样方式得到对区域形状及大小的估计，但该做法会显著增大计算量。当需要求解的区域面积较大时，为获得同等精度结果，离散点数增加而导致计算耗时激增，将大大降低计算效率。

为了解决该问题，本文利用球面几何数学方法，针对卫星视场投影区域边界与区域目标边界的特殊性，得到更高效的判断方式。

在求解过程中使用到的数学概念以及可见性的空间几何计算基础算法如图 1所示。其中，蓝色圆为一大圆，单位法向量为v′，与它平行的黄色平面为一普通圆，该黄色普通圆的圆周对球心的张角为θ。另一黄色圆为一大圆，其单位法向量为v。该圆周上有2个点，其位矢分别为d₁和d₂，而Q是蓝色大圆与黄色大圆的交点。p与p′是球面上2个位矢，其中p′是p沿着v′方向旋转一定角度得到的。

图 1中大圆是指圆心与球心重合，圆周在球表面上的圆。大圆又称测地线，因为球体上2点在球上沿测地线距离最短。普通圆是圆周在球上的任意圆，为了方便区分，若无特殊说明，本文后续所述的普通圆均不包含大圆。

算法1：求球面一位矢沿某一方向旋转角度后的矢量

如图 1所示，若要求p沿着v′(单位矢量)方向旋转α角后得到的位矢p′。则

$ \boldsymbol{p}^{\prime}=\cos (\alpha) \boldsymbol{p}+\sin (\alpha) \boldsymbol{v}^{*}. $

(1)

其中v^*=norm(v′－(v′·p)v′)，norm(·)表示矢量的单位化。

算法2：判断球面上一点是否在大圆的某一侧

由于大圆及普通圆会将球体表面划分为两个封闭区域，为判断球体表面一点是否在其中一个封闭区域上，本文将此判断方法称为算法2、3。为了描述简洁，在后续的表述中，将该点位于被圆的法向量指向一侧的封闭区域上称为在圆的“内侧”；反之，则称之为“外侧”。

如图 1所示，为判断p是否在黄色大圆的“内侧”，需计算p·v的大小。若p·v≥0，则p在该黄色大圆的“内侧”；反之，则在该大圆的“外侧”。

算法3：判断球面一点是否在普通圆的某一侧

同样地，若要判断p点是否在黄色普通圆的“内侧”，需计算p·v′的大小。

若p·v≥cosθ，则说明p在黄色普通圆的“内侧”，反之，则在该圆的“外侧”。

算法4：判断两段大圆弧是否有交点

判断2段大圆弧是否有交点的做法分为2个阶段：首先计算大圆弧所对应的两个大圆的交点；再判断该交点是否在两个圆弧上。

以图 1为例，首先计算黄色大圆与蓝色大圆的交点Q，易知Q=±norm(v×v′)，需要注意的是该交点有2个。

接下来，判断Q是否在大圆弧$\overset\frown{{\mathit{\boldsymbol{d}}_{1}}{\mathit{\boldsymbol{d}}_{2}}}$内部，则需要判断是否满足条件(d₁×Q)·(Q×d₂)≥0与Q·(d₁+d₂)≥0，若以上两个条件都成立，则说明Q在大圆弧$\overset\frown{{\mathit{\boldsymbol{d}}_{1}}{\mathit{\boldsymbol{d}}_{2}}}$内部。

算法5：判断一段大圆弧和一段普通圆弧是否有交点

算法5与算法4思路相似，判断过程仍分为两阶段，但需考虑普通圆圆心与球心不重合的问题，计算流程如表 1所述。下面用图 1以垂直于v与v′方向的侧视图图 2来作详细说明。

其中，黄色普通圆与黄色大圆的交点为Q′，并设v与v′的夹角为β。

与算法4一致，首先要求交点Q′，由于普通圆与大圆的交点可能存在多个情况，故要先讨论是否存在交点。

若β+θ=90°，则两条曲线相切，只有1个交点；若β+θ > 90°，则两条曲线有2个交点；若β+θ < 90°，则两条曲线无交点。

其中，若考虑大圆的法向单位向量有2种方向可能，β计算如下：

$ \beta=\arccos (|\boldsymbol{c} \cdot \boldsymbol{n}|) . $

(2)

若判断出有交点，则需要计算交点的位置。

设交点Q′与当前投影平面的夹角为γ，则有如下关系式成立：

$ \cos \gamma=\cos \theta / \sin \beta. $

(3)

构造平面内方向向量为v_‖=norm(v′－(v′·v)v)以及垂直面内的单位矢量v_⊥=norm(v′×v)。

则交点Q′可按照下式求出：

$ \boldsymbol{Q}^{\prime}=\cos (\gamma) \boldsymbol{v}_{\|} \pm \sin (\gamma) \boldsymbol{v}_{\perp}. $

(4)

至此，求出了交点Q′，最后需要判断交点是否在2个圆弧内。判断交点是否在普通圆弧内需要首先将点的位矢投影到普通圆平面上，之后的判断方法与算法4保持一致，在此不再赘述。

2 任一时刻区域目标可见性快速判断条件

根据节1对该问题的数学描述可以发现，当已知卫星在该时刻下的轨道信息后，问题转化为：判断球面上的区域A与区域B是否有重叠部分。其中，卫星视场投影到地球表面的区域为A，区域目标为B。

区域A的形状与卫星所搭载的传感器视场类型直接相关，且影响最终计算结果，本文研究的瞬时视场类型为：圆锥、矩形视场。区域B的定义如下：给定地球表面区域目标的边界点位置，将这些点顺次用大圆弧连接，则由该大圆弧围成的凸多边形区域即为区域B。

针对此数学问题，本文将问题拆分如下：当区域A与区域B有重叠，则只可能是以下3种情况之一：1) A在B内；2) A与B边界有交点；3) B在A内。

因此在某一时刻下判断该区域目标是否可见，只需要判断上述3种情况，若有一种情况满足，则该时刻可见；若均不满足，则该时刻不可见。

下面，将针对圆锥、矩形两种常见瞬时视场类型^{[14, 17, 19]}，分别说明其详细计算流程。

2.1 圆锥视场

圆锥视场投影到平面上为一个圆，用圆锥半角表示其视场大小。整体成像示意图如图 3所示。注意到，圆锥视场投影到地球表面后形成的视场边界是个普通圆。

根据圆锥视场按照如上所述的3种情况进行详细说明，其实际与地球成像情况如图 4所示。

针对情况1：计算过程见表 2。

其中，“极限点”是按照每条区域目标的边界定义的，如图 5所示。每次计算的“极限点”即为图中紫色点，表示的含义是：若在当前步骤，正在判断投影点在边界线“1”(图中红颜色数字)的“内侧”，则需要计算视场内与边界线“1”的距离最近与最远点，即图中的紫颜色点。

对于每条边界的判断，需要计算一次“极限点”，当判断的边界由“1”转为“2”时，其“极限点”也发生变化，如图 5所示。

针对情况2：遍历判断每条区域目标的边(大圆弧)，是否与视场投影的普通圆有交点，若有交点，则满足情况2，返回“是”；若否，则返回“否”。

针对情况3：遍历判断每个区域坐标位矢是否在视场投影的普通圆“内侧”，若全部位矢在普通圆“内侧”，则满足情况3，返回“是”；若有一位矢不在，则返回“否”。

2.2 矩形视场

矩形视场投影到平面上为矩形，整体可视空间为一个四棱锥，通过垂直半角和水平半角可以表示该四棱锥的形状和大小，如图 6所示。值得注意的是，当矩形视场投影到地球表面时，其投影区域的边界都是矩形视场的可视空间四棱锥中的各锥面与地球表面形成的交线，故该边界将变成四段普通圆弧。

根据矩形视场按照如上所述的3种情况进行详细说明，实际与地球成像的情况如图 7所示。

针对情况1：遍历判断矩形视场投影区域的边界点4个位矢是否均在区域目标区域边界的“内侧”，若全部符合，则说明满足情况1，返回“是”；若有一个点不满足，则返回“否”。

针对情况2：遍历判断每条区域目标的大圆弧边界是否与视场投影的4条普通圆弧有交点，若存在一个交点，则满足情况2，返回“是”；若否，则返回“否”。

针对情况3：遍历判断每个区域坐标位矢是否在视场投影的普通圆“外侧”，若全部位矢在普通圆“外侧”，则满足情况3，返回“是”；若有一位矢不在，则返回“否”。需要注意的是，由于矩形视场区域边界为普通圆弧，故在视场内部是普通圆弧的“外侧”。

3 可见窗口半解析快速计算算法

根据上述的任一时刻区域目标可见性快速判断条件，本文提出了一种针对区域目标快速计算可见窗口的半解析算法。该方法采用“粗略-精确变步长搜索”策略，其求解过程如图 8所示。

根据预设固定步长依次判断区域目标的可见性，并将判断结果按照时间顺序依次排序，该过程称为“粗略搜索”：每一个时刻的区域目标可见性判断方法，均使用节2所提出的快速判断条件，即判断区域A与B区域是否有重叠的3种情况。

针对前后判断结果不一致的时刻区间进行“精确搜索”，得到精确可见窗口起止时刻。该搜索过程采用文[19]提出的二分搜索方法，其中所需要的该时刻下卫星位置、速度信息，通过预处理后用6阶Lagrange插值得到。算法总流程图如图 9所示。

4 数值模拟

通过一组仿真算例来验证本文提出方法的精确性以及快速性。在精确性的验证上，本文以商用软件STK的结果为基准，对圆锥、矩形、复杂圆锥视场进行可见窗口的求解并计算其误差；在快速性的验证上，本文与经典传播算法、文[19]提出的快速算法进行对比。

对于可见窗口的精确性验证，本文选用的仿真各算例如下所述：仿真时长为1 d，仿真起始时刻及卫星轨道参数如表 3所示，轨道参数按照J2000坐标系定义，其搭载的圆锥视场载荷为圆锥半角30°、矩形视场载荷的水平及垂直半角均为30°，区域目标顶点参数如表 4所示。

计算结果如表 5所示，为了便于对比其与STK的误差结果，计算绝对误差和相对误差的箱线图如图 10和11所示。

由图 10和11可见，本文方法与STK计算结果绝对误差在0.6 s以内，相对误差在0.3%以内，平均相对误差在0.1%以内，误差水平符合工程要求。在误差的来源上，本文采用地球球体的假设，这表明误差主要是系统误差：忽略地球扁率的误差、忽略地球极移岁差影响产生的误差、卫星位置速度信息误差、STK本身误差。

在快速性的验证上，为了便于比较以及保证结果精度的展示，仿真时长在1 d或1个月(31 d)，仿真开始时刻以及轨道参数、区域参数与上述算例一致，卫星搭载的传感器载荷使用圆锥视场。在计算的过程中，使用了文[19]提出的二分法搜索策略。仿真环境为：CPU：i7-7700HQ 2.80 GHZ，内存：32.0 GB，运行平台C++。计算结果如表 6所示，其中，表格中所述的原有快速算法为文[19]提出的计算方法。

由表 6可见，本文算法的计算效率相较于现有算法有明显提升。相较于传统的传播算法，本算法在单一离散时间步长下的计算效率提高数千倍。本文提出的算法结合二分法搜索最终可以在计算效率上超过传播算法10⁵倍、超过现有快速算法近千倍。该算例的计算用时为原有快速算法的0.098%，效率提高1 024倍。这种快速性是通过将每次判断可见性时的离散方法(即将区域目标边界离散数千个点)改为解析几何方式判断直接达到的。在获得卫星位置、速度信息上，通过预处理后用6阶Lagrange插值得到。由于可见性判断模块与轨道数据计算相互独立，故本文提出算法可应用于各类轨道模型。

5 结论

本文提出了一种针对区域目标快速计算可见窗口的半解析方法，该算法将可见窗口计算耗时大幅降低。在地球球体的假设下，利用空间几何关系，将所有区域目标的边界描述为大圆弧或普通圆弧，并针对圆锥、矩形视场，得出快速判断可见性解析条件。而后引入二分搜索快速预测区域目标的可见窗口。由于每次判断可见性引入了解析判断条件，大大加快了对大区域目标的判断速度，且精度更高、适用范围更广。数值算例表明本文提出方法的具有高精度与高效率，该算法相较于商用软件STK，可见窗口持续时间相对误差满足工程要求；在计算速度上较现有快速计算算法平均速度提升显著。显著提升的计算速度将便于遥感卫星后续任务的快速处理。由于可见性判断模块和轨道数据计算相互独立，故各类轨道模型均可适用。