一种基于网格的卫星访问区域目标的时间窗口快速计算方法
2018-10-17王守斌程承旗李海东
张 玮，王守斌，程承旗，陈 波，李海东
（1. 北京大学 地球与空间科学学院，北京 100871；2. 61646部队 北京 100192；3. 北京大学 工学院，北京 100871）
0 引 言
时间窗口是指卫星经过地面目标的上空时，可对该目标进行观测的时间范围。时间窗口的计算具有重要意义，在观测卫星调度中有着重要作用[1-2]。
相对于点目标而言，区域目标的时间窗口计算更加困难：区域目标形状可能不规则，时间窗口计算通常只能利用空间采样点或区域边界线进行近似求解，而这种求解的瓶颈在于时间效率。
1 研究现状
当前，卫星对区域目标的时间窗口计算主要有如下两类方法：①基于空间采样点的区域目标时间窗口计算方法。较为典型的有：跟踪传播法[3]、星下点大圆近似迭代法[4]、视函数法[5]、庞加莱映射解析法[6]等。当采样点划分尺度和时间步长足够细时，跟踪传播法的时间窗口计算结果非常精确，但缺点是时耗很高。星下点大圆近似迭代法的算法精确度较低，同时存在无法收敛的可能。视函数法的精度很高，但区域目标空间范围较大时，采样点数量剧增，时间窗口的实时计算时效较差。庞加莱映射法基于数值计算，能精确获得时间窗口，但计算过程涉及大量积分运算，算法复杂度高，对区域目标进行计算时效率较低。②基于边界线段的区域目标时间窗口求解方法。较为典型的方法是星下点轨迹-边界求交法[1,2,7]，其优势是计算卫星对大范围的区域目标的时间窗口时相对基于空间采样点的计算方法的算法效率更高，但缺点在于计算时耗与区域目标边界的数量成正比，如果多边形区域本身范围不大但边界数量很短碎时，反而时耗会比基于空间采样点的方法更多；同时，该方法不适用于地球椭球体的情况。
上述两种方法均无法存储卫星地面覆盖信息，因此卫星访问目标计算的过程是相互独立的，如图1所示，其中存在大量实时、重复的覆盖计算，导致耗时严重。
图1 现有时间窗口计算方法：实时覆盖和独立过程计算处理Fig.1 The existing methods in time window computing access time:the real-time and independent coverage computing process
本文拟提出一种适用于卫星访问区域目标的时间窗口计算方法。这种方法将卫星覆盖信息提前存储下来，卫星访问区域目标的时间窗口实时进行复杂计算转化为覆盖信息的查询、整合，因此既能保持较高的计算精度，还具备更高的时间效率。同时，多星并行计算过程转变为针对区域目标的一次覆盖查询过程，计算耗时相对稳定，因此更具优势。
2 方法与模型
2.1 基本思路
传统卫星时间窗口计算通常采用等经纬度、等面积或等距离等方式去生成网格，再用网格点集合来近似拟合区域目标。本文同样也需要借助网格，但立足点不同：①传统方法使用的是分布均匀的离散点，网格是生成这些点的方法，而本文方法需要的是网格本身，因为网格代表了固定的小尺度地表区域；②传统方法用网格点来表征目标，本文方法用网格来存储卫星覆盖的时空信息。本文选择北京大学程承旗教授团队提出的全球等经纬度剖分网格体系（Geographical coordinate global Subdivision grid with One - dimension - integer on Two to n-th power，GeoSOT）[8]作为卫星时空覆盖网格的空间基准，剖分编码原理如图2所示。GeoSOT的优势在于：①全球尺度无缝无叠；②网格粒度适宜，且多层级嵌套；③网格编码唯一且便于计算机存储、计算。
图2 卫星地面覆盖区域剖分成同一层级网格Fig.2 The satellite coverage area subdivided into the same layer of GeoSOT grids
本文采用时间剖分编码[9]作为卫星覆盖时间信息存储载体。时间剖分编码借鉴了四叉树倒排聚合的做法，将时间剖分形成同时支持单尺度和多尺度表达的64位二进制1维编码形式，既方便数据库存储和索引，又能基于编码进行时间聚合和跨度运算。同时，这种编码还具有十进制表达数值越小、时间越早，数值越大、时间越晚的特点，因此适合于时刻排序过程。
该方法需要计算瞬时时刻卫星地面覆盖区域。目前，卫星地面覆盖计算方法研究较多，主要有球面三角形法和矢量法两大类。本文主要采用矢量法[10]进行卫星地面覆盖计算。矢量法的优势在于计算精度较高，同时兼容地球球体和地球椭球体的情形，且还兼容不同视场形状的卫星传感器。
主要思路是：①提前计算卫星地面覆盖，剖分化后用网格记录卫星的时空覆盖信息，存储至数据库中；②区域目标访问分析计算转化为按边界网格遍历查询数据库的过程，寻找卫星覆盖区域目标边界的所有时刻，生成时间窗口；③多星分析区域目标时，只需一次遍历查询即可获得所有卫星访问区域目标的所有时间窗口。技术路线流程图如图3所示。
图3 技术路线流程图Fig.3 Technical route fl ow chart
2.2 卫星地面覆盖的预存储
卫星地面覆盖预存储是指将仿真时域按等间隔划分成一系列时间区间，依次计算时间区间的边界时刻下卫星地面覆盖区域，并将覆盖区域网格化，将该时刻的所有覆盖网格的网格编码附上该时刻生成的时间剖分编码，存入数据库，最终形成卫星覆盖数据库的过程。
2.2.1 仿真时域的划分
假设仿真时域的起始时刻为Tbegin，终止时刻为Tend，设定时间步长为Δt，那么该时域可以分割成数量为N的小时域。从第0段开始到第N-1段，时间段的起始和终止时刻分别为：[Tbegin,Tbegin+Δt],[Tbegin+Δt,Tbegin+2Δt],…,[Tend-Δt,Tend]。其中，前一个小时域的终止时刻是后一个小时域的起始时刻，那么将起始时刻为Tbegin、终止时刻为Tend以及所有分割时刻集合，可得
此时，T被称为仿真时域的时刻集。
2.2.2 卫星瞬时地面覆盖区域计算与剖分化
遍历T中的时刻：
Step1：求解该时刻下卫星的地面覆盖区域，结果用区域边界点的经纬度序列CoverLine _pt_seq表示。
Step2：选择符合卫星应用精度需求的最高网格层级（即基础层级，一般可选择第14层级，该层级网格大小在赤道附近约为4 km左右），将计算得到的卫星传感器地面覆盖区域剖分至该网格层级，如图2所示，得到一个由网格组成的集合Cs,t-＞e；用该集合来表征卫星传感器此时刻的地面覆盖，即：
Step3：为了尽可能减少数据量，将Cs,t-＞e中的网格聚合，形成多尺度网格集合Cs,t-＞e,multi，如图4所示。
图4 地面覆盖网格聚合Fig.4 Satellite coverage grids merged into multi-layer coverage grids
2.2.3 卫星覆盖时空网格编码入库
遍历T中的时刻Ti，都能得到与之对应的Cs,t-＞e,multi。将Cs,t-＞e,multi编码化可以得到多层级的编码集合Gs,t→e：
同时，将Ti编码化为Tcode，那么结合Gs,t→e中的每个空间网格编码元素，都能形成一个
卫星覆盖数据库的表结构字段应当包括：
1）卫星编号（satID）：INT32类型；
2）传感器编号（senID）：INT32类型；
3）卫星传感器名称（Name）：String类型；
4）空间网格编码（gcode）：UNSIGNED LONG LONG或 __INT64类型；
5）时间剖分编码（tcode）：UNSIGNED LONG LONG或 __INT64类型；
6）其他属性码。
其中，数据库中每个数据条目代表的含义是：某颗卫星的某个传感器在某时刻（或分割时域，取决于时间精度）能够覆盖某个网格。此时，卫星传感器地面覆盖网格的映射关系全部打散、存储至数据库的一张表中，这张表的条目数量就是所有卫星在所有时刻的所有网格编码的总和。在此基础上，基于空间网格编码建立数据库正序普通索引。
2.3 区域目标的时间窗口计算
区域目标的时间窗口分析的基本做法是，先将区域目标的边界点序列转换成基础层级的网格，再将网格代入数据库中查询，返回与网格相关的所有时间编码，最后将时间编码整合生成时间窗口。
2.3.1 区域目标边界剖分化
以点序列输入的区域目标边界后，将点序列内相邻点两两连接，形成若干条边界线段，依次遍历所有线段，将线段转换成基础层级下的线段网格集合，如图5所示；最后将所有集合合并，得到区域目标的边界网格编码集合。
图5 区域目标边界剖分化Fig.5 The area target boundary differentiation
2.3.2 数据库查询
将边界网格编码集合中的网格编码依次代入卫星覆盖数据库中查询，获得覆盖所有边界网格时的时间编码。由于数据库存储的是多层级网格时空覆盖信息，因此如果该网格能被卫星覆盖，可能存在两种情况：
1）该网格编码存在于数据库中；
2）该网格编码不存在于数据库，但该网格的多层级父网格存在于数据库中。
因此，应当使用代入网格编码，如下两个条件满足其一即返回的查询条件进行查询：
条件1. 等值判断，空间编码与查询输入相等，此时该输入网格存在于数据库中；
条件2. 异或求交判断，查询输入与空间编码的经向编码与纬向编码进行一次异或运算，并对两次异或运算结果进行求交运算，结果值为1。
按2.4.1节得到的边界网格集合是顺序排列的。考虑到卫星在仿真时间内进出区域目标的次数（假设为N次）是有限的，卫星对区域边界的覆盖集中在2N个边界部分，在这些部分之外的区域边界不被覆盖。因此，在卫星覆盖数据库查询中，可以采用跳跃查询的方式缩减查询次数。具体做法是：
Step 1：设置合适的跳跃步长step，step=2n，设当前查询网格的集合索引值为record，设回退状态布尔值为isStep为false；
Step 2：从第0个网格开始遍历，即record初始化为0：
Step ①：将record指向的网格编码代入数据库中查询；
Step ②：如果返回为空，则判断record+step是否大于边界网格集合的元素个数，如果不是，则record+=step，进入Step①，如果是，则record指向最后一个元素，进入Step①；如果返回不为空，则判断isStep是否为true，如果为true，则进入回退搜索临近最远的非空返回的索引值，如图6所示，如果为false，则记录该网格的查询返回时间编码，如果此时record指向最后一个元素，则进入Step3，否则record+1，进入Step①；
Step 3：遍历结束，将时间编码集合并去重。
图6 回退搜索过程示意图Fig.6 The processing of the rollback search method
2.3.3 生成时间窗口
获得所有覆盖边界网格的时间编码后，按值大小排序。此时，相邻元素之间的时间差值存在3种情况：①差值等于时间步长（如1 s）；②差值略大于时间步长（如10 s）；③差值远大于时间步长（如30 min）。情况①说明相邻元素均属于同一个时间窗口，较早者可能为时间窗口的起始时刻；情况②说明相邻元素属于同一个时间窗口，较早者为卫星进入区域时覆盖边界的最晚时刻，较晚者为卫星离开区域时覆盖边界的最早时刻；情况③说明相邻元素不属于同一个时间窗口，较早者为前一个时间窗口的终止时刻，较晚者为后一个时间窗口的起始时刻。
基于上述情况，生成时间窗口的主要做法是利用区域边界的经纬度极值，生成区域目标的最小外包矩形，计算外包矩形的较长对角大圆弧长度，再除以卫星星下点轨迹运动平均速度，以此预估卫星一次覆盖区域目标可能出现的时间窗口最大值△T；以该预估值为临界值，将时间编码集合分簇，然后获得各簇的最小和最大编码，作为时间窗口的起始时刻和终止时刻，如图7所示。
图7 时刻分簇生成时间窗口示意图Fig.7 Time clustering generation time windows
3 仿真实验
基于单星的地面覆盖预存储结果，分析计算两个假想区域目标的时间窗口，并以跟踪传播法作为比较对象，对比分析本文方法的结果精度和时间效率；同时，对多星地面覆盖进行存储，计算一个区域目标的时间窗口，对比分析本文方法和传统方法多星计算的时间效率。
3.1 实验环境
实验的硬件环境和软件环境分别见表1和表2。
表1 硬件环境Tab.1 Hardware Environment of the simulation experiment
表2 软件环境Tab.2 Software Environment of the simulation experiment
3.2 实验数据
使用20颗卫星进行实验，其中单星分析所用的卫星轨道根数见表3。仿真区域目标为两个，目标1为中国大陆区域（不含中国台湾、海南等岛屿、海域），边界点数量为4 626，目标2为中国南海部分区域，边界点经纬度坐标见表4，不考虑卫星的姿态角变化，假定视场中心点投影到地表处与星下点始终重合，同时设定水平半视场角和垂直半视场角均为15°。仿真时域的起始时刻为2018年5月21日9点整（UT0），终止时刻为2018年5月22日9点整（UT0），时域跨度为1 d（即86 400 s）。
表3 卫星轨道根数Tab.3 Orbital elements of the satellite in the single-satellite experiment
表4 仿真区域目标边界点经纬度Tab.4 Geographical coordinates of the boundary points of the area targets
3.3 实验结果与分析
3.3.1 单星计算实验结果与分析
在卫星覆盖预存储阶段，选定时间步长为1 s，覆盖网格基础层级为GeoSOT第14层级（2′网格，网格大小对应赤道附近4 km左右），计算卫星覆盖时空网格并完成入库。输入目标1和目标2的边界点序列，分析两者的时间窗口，计算结果见表5。
表5 目标1和目标2的时间窗口计算结果Tab.5 Time window calculation results of target 1 and target 2
将目标2的起始时刻和终止时刻转换为时间编码，代入数据库中查询相关网格编码，再将网格绘制成图，与STK二维标准输出窗口结果进行比对可知，本文采用的数据库查询方法是可靠的，如图8所示。
图8 本文方法结果与传统方法结果的对比Fig.8 Comparison of the results of this method with the traditional method
以STK分析结果为标准值，采用绝对误差Δ和相对误差δ[1]作为分析结果精度的指标。
绝对误差的计算公式为：
式(4)、(5)中，N为时间窗口的个数，tSi,true和tEi,true分别是第i个时间窗口真实值的起始时刻和终止时刻，tSi和tEi分别是通过本文方法得到的第i个时间窗口的起始时刻和终止时刻。
绝对误差反映了计算所得的时间窗口与实际时间窗口在数值上的绝对差别。相对误差则反映了误差相比于时间窗口的大小。
取跟踪传播法的时间步长为1 s，计算时间窗口，将所得结果与本文方法进行对比，计算时耗与结果精度见表6和表7。
表6 本文方法与跟踪传播法的误差对比Tab.6 Time error comparison between this method and the tracking propagation method
表7 本文方法与跟踪传播法的计算效率对比Tab.7 Comparison of computational efficiency between this method and the tracking propagation method
由表6可知，与跟踪传播法相比，本文方法的绝对误差为秒级误差，与其他快速计算方法持平[11]。误差主要取决于如下3个方面：
1）时间步长。时间步长越小，仿真时域的划分就越细，相邻时刻计算得到的卫星地面覆盖区域重合度越高，单个覆盖网格包含的时间覆盖信息就越丰富（原理类似于视频帧数）；反之，则单个覆盖网格包含的时间覆盖信息就越少。当时间步长取值为1 s时，卫星在1 s时间内的覆盖动态变化信息就无法获取，那么网格只能记录1 s时间能否被卫星覆盖，至于这1 s内是否都能被覆盖，无从得知。因此，时间步长直接决定了时间窗口分析精度的分辨率，即基于1 s步长的预存储覆盖数据进行区域目标查询分析的结果绝对误差不可能小于1 s。
2）网格基础层级。本文提出的方法本质上是将卫星覆盖能力提前计算、存储至剖分网格中，剖分网格是卫星时空覆盖信息的载体。网格划分越细，对卫星瞬时覆盖区域的拟合就越好，基于区域边界网格的时间窗口查询分析结果也越收敛、逼近真实时间窗口。
3）卫星瞬时覆盖的计算精度。卫星瞬时覆盖的结果精度也会影响本文方法的结果精度。瞬时覆盖精度越高，数据库存储的时空覆盖信息就越接近真实情况，最终经过查询检索得到的时间窗口也越精确。
表7可知，本文提出的方法计算效率相对于传统的跟踪传播法更高，主要原因是该方法将复杂覆盖计算变成了预存储过程，面向具体区域目标时，只需要判断网格与网格之间的包含关系，而GeoSOT采用了64位二进制Z序网格编码的形式，网格包含关系都可用编码位运算解决，因此区域目标的时间窗口计算过程被分离、简化，得以提升时间效率。
3.3.2 多星计算实验结果与分析
将20颗卫星地面覆盖存储，对目标2进行多星计算。以跟踪传播法为对比，计算效率结果如图3-2所示，横坐标为卫星数量，纵坐标为传统方法计算耗时与本文方法耗时的比率。
图9 多星分析效率比较Fig.9 Computational efficiency comparison of the multi-satellite analysis
由图9可知，随着卫星数量的增加，本文方法的时间效率比传统方法越来越高，其原因在于多星地面覆盖都用网格进行存储，而基于边界网格进行覆盖时刻查询时，可以做到一次遍历、同时获取所有卫星的时刻，进而得到所有卫星的时间窗口，因此用本文方法进行多星计算时，时间损耗是较为稳定的；而传统方法进行多星计算时，需要分别对每颗卫星计算访问目标的时间窗口，因此时间损耗随着卫星数量的增加而增加。因此，本文方法更适合于多星计算的情形。
4 结束语
本文提出了一种基于网格的卫星访问区域目标的时间窗口计算方法，即将卫星地面覆盖提前计算、剖分化形成时空覆盖网格信息，存储至数据库中，在计算区域目标的时间窗口时，将区域边界剖分得到的网格遍历代入数据库查询，寻找卫星覆盖区域目标边界的所有时刻，生成时间窗口。该方法具有卫星覆盖信息“一次存储、反复查询”的技术优势，仿真试验证明，该方法比传统方法更适用于多星时间窗口计算，因此可为多星联合快速调度分析提供技术支撑。该方法在卫星覆盖计算及网格存储方面，未考虑敏捷卫星姿态角变化、传感器侧摆等因素，同时需要消耗一定的计算机存储空间。因此，下一步工作需要改进时空覆盖网格存储的方法，以便兼容卫星姿态角、传感器侧摆等复杂情况下的卫星覆盖信息存储，同时还要进一步压缩数据库大小，以便节省更多的空间开支。
