石 习，冯占林，赵会朋，秦晓珊

(中国电子科技集团公司 电子科学研究院，北京 100026)

0 引言

低轨卫星凭借重量小、成本低、通信时延低等优势，逐渐成为各国空间领域研究的热潮。成百上千颗的卫星按照预先设定的构型，以组网的形式分布在不同轨道平面以及相位上，协同完成通信、导航、侦察遥感等任务。其中最具代表性的当属“星链”低轨卫星星座。“星链”卫星轨道高度低、周期短、重访率高，携带传感器后可对全球主要地区进行24小时侦察监视。轨道高度为550 km的“星链”卫星能够以最大44.85°的半张角对地探测，覆盖半径可达573.5 km。覆盖分析是评估卫星对地侦察监视能力的基础，计算方法主要有网格点法和基于几何运算的方法。

常用来分析卫星覆盖性能的指标主要有：覆盖时间、重访时间、覆盖面积以及N重覆盖等等。覆盖分析的经典算法网格点法由Morrison，J.J于1973年首次提出[1]。该方法基于特定的规则将目标区域划分为一系列网格，通过对网格点的覆盖性能统计表述卫星对区域的覆盖性能。划分规则决定着计算结果精度和计算量。文献[2-4]分别对网格点法做出了改进。基于几何运算的方法就是在经纬度平面上利用多边形交、并运算来计算覆盖性能[5]。文献[6-7]也通过多边形计算覆盖相关性能。

本文关注卫星对地覆盖性能中的覆盖时间这一指标，也即卫星对地面目标的可见时间窗口，需要计算卫星在未来特定时间内对于该目标的入境、出境时间。传统的时间窗口计算方法主要是跟踪传播模型，该模型需要对卫星位置连续采样，即对卫星进行轨道预报，并判断当前位置卫星与地面目标的可见关系，如果可见则记录当前时刻，最后形成的可见时间段即为卫星对地可见时间窗口。该模型采取固定步长的方式，当预报周期增加，此方法耗时严重。

针对覆盖时间窗口传统模型计算量大、所耗时间长的问题，解决方法大致可分以下几种：

1)通过轨道过滤、数学模型简化等近似方法。I.Ali等利用球面几何计算了圆形低地球轨道卫星对地面终端的可见性时间函数，其迭代计算方法速度快，但算法精度较低且只能针对点目标区域[8]。李冬等通过大圆近似星下点轨迹，通过求解关于偏近点角的超越方程从而计算时间窗口[9]。该方法误差主要来自于岁差，章动以及地球形状等。Feuerstein P等首先过滤掉不可见的轨道周期，减少了星地链路可见性时间窗口的计算量[10]。宋志明等首先在不考虑地球自转的情况下计算卫星对地可见时间窗口，然后根据地球自转特征迭代修正，从而得到整个周期内的时间窗口[11]。该模型未加入地球非球形引力等其他摄动因素，且无法对区域目标进行计算。

2)通过分布式，实时服务、查询方法。张玮等将区域目标时间窗口计算分离成提前计算并编码的预存储阶段和遍历查询的检索获取阶段[12]。卢万杰等提出了星地可见时间窗口计算的实时服务方法[13]，首先通过地面目标区域与星下点的空间关系，从而决定是否计算当前时刻覆盖。在此基础上构建时间窗口计算流程，根据系统CPU核数设置计算并行的任务数，并基于Storm实现实时数据处理，保证了计算与服务的实时性。

3)通过动态设置步长的方法，即通过由粗略到精细的搜索方式。张众等通过空间几何关系把区域目标的边界描述为圆弧，判定边界是否相交；针对卫星视场得出可见性的解析判断条件，最后用二分搜索方法计算卫星对区域目标可见窗口[14]。鄂智博等首先使用较大的步长，得到卫星对地面点目标的低精度可见时间窗口，然后使用二分搜索方法得到精确窗口计算结果。对于区域目标而言，则通过统计边界点的可见时间窗口的上下界得到[15]。汪荣峰等定义“预测距离”的概念[16]，分别将卫星与地下点距离和距离的一半作为预测距离，计算动态步长下卫星对点目标区域的可见时间窗口，大大减少了采样点数目，提高了效率。最后针对地面区域目标，给出预测距离的计算方法。

本文提出了一种基于经纬度平面上多边形距离的覆盖时间窗口快速计算方法。算法流程是：以观测矢量和地球的位置情况(包括相交、相切、不相交三种情况)计算卫星对地球覆盖边界坐标，将其转化为二维平面内的点或多边形，通过与目标区域多边形位置关系判断从而计算卫星对目标区域的出、入境时间。在此基础上，重新定义“预测距离”的概念，将多边形的距离作为动态调整算法步长的依据，避免了传统跟踪传播模型采样点多、计算慢的问题，实现了对地覆盖时间窗口的快速计算。

1 卫星对地覆盖建模

卫星携带传感器后会对地形成一个探测区域，如图1所示。传统的基于球面三角形的覆盖模型是一种探测区域的近似求解方法。如果要考虑传感器的形状、姿态以及半张角等因素，则可以使用观测矢量模型[17]。使用观测矢量模型进行覆盖范围边界计算的步骤是：首先根据传感器形状设置传感器边界点个数，建立卫星对地观测矢量，在地心固定坐标系下联立观测矢量和地球椭球方程，求解观测矢量在地球表面的映射位置，从而实现覆盖边界计算。

图1 卫星对地覆盖示意图

对于观测矢量和地球不相交的情况，需要计算相切时的观测矢量和地球的交点位置，作为不相交情况下的卫星对地覆盖边界，具体推导参考文献[18]。在卫星工具软件STK中生成一颗卫星，设置传感器半张角、姿态角等参数后生成会对地球表面的覆盖区域。利用卫星对地覆盖模型计算出的边界点坐标生成STK中的AreaTarget对象(记作MyAreaTarget)，从图2中可以看出，两区域几乎重合。

图2 对地覆盖仿真二维图

图3 射线法示意图

如果模型中传感器边界点个数较大，那么生成的覆盖区域边界点数目也会较多，手动添加边界点经纬度坐标至STK中的步骤就会愈加繁琐。为了避免这种重复性软件设置操作，接下来介绍将覆盖模型得到的一系列点坐标生成STK中AreaTarget对象的简单方法，具体通过STK提供的MATLAB接口发出控制指令实现。具体代码如下。

代码1：建立卫星覆盖区域

ui=actxGetRunningServer(‘STK11.application’);

root = ui.Personality2;

sc = root.CurrentScenario;

a=sc.Children.New(‘eAreaTarget’,‘MyAreaTarget’);

root.BeginUpdate();

boundary = {26.58 124.26

28.63 126.34

... ...

26.29 128.01

24.18 126.67;}

a.CommonTasks.SetAreaTypePattern(boundary);

root.EndUpdate();

2 卫星对地时间窗口计算方法

求解空间目标在一定周期内对目标区域覆盖时间窗口的集合，是地面目标预警和干预的基础。预报结果基于卫星对地覆盖范围求解结果，计算卫星对地面目标区域的可见性弧段，求解卫星在一定时间内经过目标区域的时间，可以为空间预警和航天侦察规避提供数据信息。

卫星是否携带传感器以及传感器的形状和参数都会直接影响时间窗口的计算结果。计算过程中，对未携带传感器或传感器半张角为0的情况，卫星对地的覆盖为星下点，卫星是否经过目标区域，可以理解为经纬度平面上星下点是否在目标区域之内，即判断点和多边形的平面位置关系。其他情况下，卫星对地球表面形成一个覆盖区域，此时需要判断其与目标区域之间是否存在重叠。

2.1 位置关系判断

2.1.1 点和多边形位置关系判断

点与多边形的位置关系判断有多种方法，本文使用射线法。

射线法的核心思想是：沿目标点任意方向做一条射线，通过射线与多边形边的交点数量判断点是否包含在多边形内部。交点数量为奇数表示在多边形内，偶数表示不在多边形内部。对于几种特殊情况的规定与处理，可以参考文献[19]。这里为简化计算，取射线方向为Y轴负方向，此时射线法示意图如下。

本文将多边形边界点作为内部点处理。

此时的计算流程如图4所示。

图4 射线法计算流程图

2.1.2 多边形位置关系判断

若覆盖区域和目标区域有重合部分，则认为卫星对目标区域进行了覆盖，重合的情况应该包括相交和包含关系。对于相交情况，可以通过MATLAB的POLYXPOLY函数实现，该函数返回平面笛卡尔坐标系中两个多边形的交点(注意这里的多边形也可以是线段，说明可以满足线阵式传感器的情况)。这里利用该函数返回值是否为空作为两个多边形是否有交点的判断依据。

将卫星覆盖区域抽象为二维坐标系下多边形时，需要对以下两种特殊情况进行处理：

当卫星对地覆盖区域包含极点时，需要添加极点作为多边形顶点辅助生成覆盖区域多边形；当卫星对地覆盖区域穿过180°经线时，需要添加180°经线上的一系列点将覆盖区域划分为两块，分别生成两个多边形。此时需要分别对这两个多边形与目标区域多边形进行位置关系判断，可参考文献[20]。

2.2 时间窗口计算

依据卫星轨道预报结果和传感器对地覆盖计算结果，可以实现对卫星是否经过目标区域的出、入境时间。算法流程如图5所示。

图5 时间窗口计算流程图

分析模型的误差主要来自于卫星对地覆盖模型的计算结果，对于简单圆锥形传感器而言，由于实际覆盖区域是一个圆形区域，因此可以通过增加多边形的顶点数来减小覆盖区域实际区域的误差，从而降低最后算法时间窗口结果的误差。

由观测矢量模型可知，多边形的顶点可以通过传感器边界点的个数进行修改。将模型中传感器圆边界点个数增加，此时对覆盖区域的拟合效果更好，算法计算的时间窗口精度越高，但算法耗时也随之增加。同时，固定步长1 s下，对卫星位置的采样点数目与计算周期成正比，算法时间耗时也随之增加。因此，需要对计算时间进行优化。

3 动态步长的快速计算方法

跟踪传播模型以一定步长连续计算预报周期内卫星对地覆盖区域边界点。当步长过大，卫星位置采样数量减少，相邻采样点对应覆盖多边形间隔越大，可能导致错过某些时刻对于卫星是否过境的判断，造成卫星过境预报的“漏警”，计算得到的出入境时间精度越低。当步长过小，容易造成整个预报周期内很多“无意义”的计算和判断。例如，一天内“星链”对乌克兰地区的单个覆盖时间窗口通常在几分钟内，其他大部分时间内，卫星覆盖区域与目标区域没有重合，甚至相距较远。对一天时间的预报周期，如果按以1秒为步长的采样方法，无疑会进行多次不必要的操作，导致算法计算量大、效率低。

卫星对地时间窗口算法耗时与卫星采样点的个数呈正相关。因此，实际有意义的计算只发生在相交时间及相交边界的领域。如何避免无意义的计算，使采样时间尽快“收敛”到相交时刻附近，是本文研究的重点。

以编号为44238的“星链”卫星为例，下载该卫星某天的TLE数据，解析出卫星历元时刻为7 Oct 2021 11：40：04.473(UTC)，计算一天内该卫星对地面某目标区域的可见时间窗口。

将一段时间内卫星对地覆盖区域的多边形投影到二维平面上，截取部分时刻的覆盖区域和目标区域多边形位置关系如图6和图7。当覆盖区域接近目标区域或者覆盖区域与目标区域相交时，此时步长应该设置较小，避免漏掉其他相交时刻的判断，导致过境时间精度低。

图6 覆盖区域与目标区域位置图

图7 覆盖区域与目标区域位置图(相距较远)

某时刻卫星覆盖区域和目标区域相距较远，位置关系如图7所示。此时步长应该设置较大，避免很多不必要的计算与判断，导致计算耗时严重。

至此，本文定义二维平面内覆盖区域和目标区域的预测距离为动态选取步长的因素，并做出如下定义：当两区域有重合时，预测距离为0；当两区域没有重合时，预测距离为平面内两多边形的最短距离。计算二维平面内覆盖区域和目标区域的瞬时预测距离的函数如下：

算法:多边形预测距离

输入:二维平面内覆盖多边形、目标区域多边形

输出:覆盖区域和目标区域的瞬时预测距离Pd

if(X与Y相交)

Pd=0

else

对X的第i条边,计算各顶点Yj与该边的最小距离Dj

对Y的第j条边,计算X各顶点Xi与该边的最小距离DjPd=min(Di,Dj)

end

计算上文两种不相交情况下的多边形预测距离，最短距离如图8和图9所示。

图8 不相交情况最小距离示意图(一)

图9 不相交情况最小距离示意图(二)

按照定义，此时预测距离分别为0.168 0和6.488 7。

选取预报周期内一段时间(包含覆盖时间)，计算预测距离随时间变化，结果如图10所示，其中横坐标表示距离预报起始时间(卫星历元时刻)过去的秒数。

图10 预测距离随时间变化图(部分)

在该时间段内，随着卫星位置的变化，预测距离先由70下降至0后再次上升。由前文预测距离的定义可知，预测距离为0即表示覆盖区域多边形和目标区域多边形有重合，也即卫星对目标区域可见。

于是，计算卫星出入境时间的问题转化为预测距离的函数等于最小值的优化问题。

由图可见，在选取的这段时间内，预测距离下降呈现一定的规律性，梯度较为平稳。

依据预测距离的变化率动态设置步长，其中步长最小取1秒，即当预测距离接近0时，卫星采样步长为1秒；预测距离较大时，卫星采样步长可达几至几十分钟。

计算预报周期内卫星对地可见时间窗口，对比固定步长为1秒的算法，结果如表1所示。

表1 两种方法计算结果对比

对比结果可知，两种方法计算卫星对地可见性时间窗口完全一致，动态设置步长后，参与覆盖边界计算、与目标区域相交判断的采样点数目由86 500个减少为457个。算法耗时与采样点数目成正比，故动态设置步长后，相比于固定步长方法效率提升约99.5%。

卫星对地可见时间窗口结果同时取决于卫星轨道、传感器参数和目标区域位置等因素。选取“星链”卫星中轨道倾角更高的极地轨道卫星，编号48 879，计算固定步长和动态设置步长两种方法下卫星对目标区域的可见时间窗口相关结果如表2所示。

表2 两种方法计算结果对比(44 879)

对于极地轨道卫星48 879而言，预报周期内对地覆盖时间较短，动态设置步长方法所采样点数目更少，算法效率提升更高。

4 计算结果对比

STK是个功能强大的仿真分析软件，而且兼具场景、模型逼真的特点、在卫星、导弹等领域应用广泛。但软件逻辑能力较差，无法通过编程实现循环计算、嵌套等复杂过程。通常情况下，航天应用场景规模庞大，需要进行手动重复性操作时便不再方便。为此，STK提供了编程接口便于其他应用程序调用STK，借助接口可以编程实现复杂任务的自动化实现。

为了验证算法的准确性，可以和STK结果对比。利用STK创建仿真场景，导入100 颗“星链”卫星的TLE数据，并设置仿真周期为一天。

卫星对地时间窗口算法的相关参数设置如下：传感器边界点数目设置为100，半锥角为44.85°，姿态角均为0。在STK中需要对每颗卫星添加和算法相同参数的传感器，这里依然使用STK提供的MATLAB接口编程实现。

本文在向STK导入卫星对象时获取添加路径，通过创建constellation集合(记作MyConst)存放所有的传感器，利用卫星路径和集合元素的简单映射关系，循环为每个卫星添加传感器。以下是MATLAB实现的关键代码。

代码2：自动为卫星添加传感器

cs=root.CurrentScenario.Children.New(‘eConstellation’,‘MyConst’);

%创建星座集合

%循环给卫星添加传感器

for i=1:count

sat=root.GetObjectFromPath(char(satPaths(i)));%获取卫星路径

Sen=sat.Children.New(‘eSensor’,‘Sen’);

Sen.CommonTasks.SetPatternSimpleConic(44.85,0.1);

cs.Objects.Add(Sensor.Path);%添加传感器路径

end

分析“星链”星座对区域的覆盖性能需要使用STK的覆盖模块(STK/Coverage)，分析对象主要是Coverage Definition，利用该模块可以设置“覆盖资源”(比如卫星或者传感器)、覆盖区域(指定经纬度区域或者全球)、覆盖网格(Grid)以及各种约束条件。

覆盖对象设置具体步骤如下：(1)在场景中添加Coverage Definition对象covDef；(2)设置Grid：本文使用Custom Regions自定义区域类型，关联场景中建立好的AreaTarget，并设定Point Granularity(网格点间隔)为Lat/Lon=0.5°；(3)设置Assets(覆盖资源)：选择预先建立的MyConst星座，添加至覆盖资源。以下是MATLAB设置覆盖对象的部分代码。

代码3：覆盖对象设置

covDef=sc.Children.New(‘eCoverageDefinition’,‘CovDef’);

covDef.Grid.BoundsType=‘eBoundsCustomRegions’;

covGrid=covDef.Grid;

bounds=covGrid.Bounds;

bounds.AreaTargets.Add(‘AreaTarget/’MyAreaTarget);

Res=covGrid.Resolution;

Res.LatLon=.5;%网格点经纬度大小

covDef.AssetList.Add(constellation.Path);

至此，100颗携带传感器的“星链”卫星覆盖场景建立完毕。为了验证本文卫星对地时间窗口算法的准确性，将STK分析结果中的覆盖起始时间和结束时间与算法结果中的入境时间和出境时间取差的绝对值作为误差，如果某颗卫星对目标区域不止有一个覆盖时间窗口，则取多个窗口误差的均值作为最终误差，入境时间和出境时间的误差结果如图11和图12所示。

图11 入境时间误差图

图12 出境时间误差图

从计算结果和误差比较图可以看出，本文算法获得的时间窗口准确性较高，结果与STK仿真结果之间的误差多在1 s以内。

5 结束语

本章首先构建了卫星对地覆盖模型，借助该模型可精确求解卫星对地瞬时覆盖区域的边界点，并利用STK软件验证了模型的正确性。本文模型考虑了传感器形状、视场角和姿态角因素，适用于观测矢量与地球表面相交、相切、不相交三种情况。

卫星对地覆盖模型计算卫星当前位置下，传感器对地球表面探测区域的覆盖边界点坐标。计算一段时间内的覆盖区域和卫星对地可见时间窗口，需要对卫星位置进行轨道预报，且最后结果误差受轨道预报误差的影响。“星链”卫星轨道周期低于225分钟，属于低轨目标，与其TLE数据相适应的轨道预报模型是SGP4模型。对于不同轨道特征的卫星应选择相适应的轨道预报模型。对不同轨道高度卫星的轨道预报模型以及模型误差分析是本文后续计划开展的工作。

其次，提出了一套计算卫星对地目标区域可见性分析的方法，可以求解卫星对目标区域的过境时间窗口，结果基于二维平面上点与多边形或者多边形与多边形之间的位置关系，并且可通过改变传感器边界点的个数得到了不同精度的时间窗口计算结果。随着预报时间的增长，按一定步长逐步采样计算的方法耗时严重。最后，本文在相关研究的基础上，将平面内多边形之间的距离作为对步长动态设置的依据，在保证计算精度的同时提升了卫星对地可见时间窗口算法的效率。对于步长设置过程的优化方法，也是后续研究的方向和重点。