淡 鹏， 何雨帆， 陈俊收

(1.宇航动力学国家重点实验室， 西安 710043； 2.西安卫星测控中心， 西安 710043)

1 引言

航天器升力式再入返回[1-2]地球过程中可通过主动控制倾侧角和速度攻角来调节飞行过程中的升力和阻力，使其获得较大的横向和纵向机动能力[3]。 为了衡量航天器再入飞行的机动性能或寻找合适的安控时机，常常需要实时计算飞行过程中的覆盖范围(例如飞行过程的可达区、溅落区等)及其相关几何参数[4]。

航天器再入飞行的覆盖范围定义了航天器潜在的覆盖区域[4-6]，可为航天器再入飞行任务规划提供关键信息，同时可帮助航天器在正常的再入飞行和突发故障时变更或选择潜在的着陆位置。 在覆盖范围分析时需要实时计算该覆盖范围与特定的区域(如期望的落点限制区域、国境线、海域等)相交部分的面积或面积占比等参数，以便对当前飞行能力及飞行的安全性进行实时评估。 此类计算实质上可转化为计算球面上不规则多边形区域的面积以及多个多边形相交部分求解等问题。 最常见的多是二维平面上的算法。 例如葛伟华等[7]基于边界跟踪方法实现了一种平面上的区域面积计算方法；齐澄宇等[8]在GeoST 全球参考网格基础上，设计了一种可以摆脱经纬度和积分运算的面积计算方法。 这些方法主要针对平面图形，应用到球面时误差扩大。 对球面上的计算问题，王会然等[9]利用图斑椭球和高斯投影方法减小了面积计算的误差；Kumar 等[10]给出了一种采用近似展开公式计算椭球面上面积的方法。 但以上方法还稍显复杂，编程实现有部分难度。

多边形相交区域面积计算是该问题更深层次的计算，它需要首先找出不规则区域的相交部分，然后再计算出面积等参数。 针对该问题，朱磊等[11]提出了一种平面上相交面积计算的双向链表算法；汪荣峰等[12]基于多边形的交、差运算，将覆盖多边形进行分解和三角化，采用统计算法实现了一种基于多边形布尔运算的相交面积计算方法。 但以上方法实现时还不够简便，对于覆盖区面积占比计算时的实时性有待进一步提高。

考虑到航天器再入飞行过程中，为了及时寻找合适的安控时机，覆盖区面积及占比计算常常需要较高的实时性和对各种不规则形状的适应能力，需要快速对航天器进行飞行状态的判断。 因而需要研究一些误差较小、实现过程简洁、实时性高、可靠性强，且支持复杂几何边界和易于工程应用的计算方法。 本文提出的方法与常用的经度、纬度二维方向网格划分算法[13]不同，仅在纬度方向进行一维分割的区域面积及面积占比计算。 该方法直接建立在球面而非平面形状之上，利用球冠面积公式和弧段张角的统计思路，实现覆盖区相交部分面积及面积占比的快速估计。

2 覆盖区相交面积及占比快速估计

在航天器返回覆盖区与特定区域(此处以极不规则的国境线为例)相交部分的面积及面积占比快速估算中，传统的网格划分算法[13]因为要进行经度、纬度二维方向上的划分，使得分割块数较多时，计算的复杂度成二次方增长，计算速度明显下降，因而较难满足实时性的要求。 本文对其进行改进，只在纬度方向上进行分割，使计算复杂度降低了1 个维数，相较于网格算法的计算速度有较大提升。

计算覆盖区相交面积的算法流程如图1所示。

图1 覆盖区相交面积计算流图Fig.1 Flow chart for calculating the intersection area of landing footprint

2.1 界线数据的准备

航天器返回覆盖区进入国境内部分面积及面积占比计算中需要用到2 个球面多边形，分别是地面覆盖区域边界多边形和国境线多边形。

设2 个多边形分别用一系列逆时针顺序排列的多个点表示。 其中，地面覆盖区域边界多边形顶点个数为n0，国境线多边形顶点个数为n1。 设地面覆盖区域边界多边形顶点序列为p0，p1…pi…pn0( 0 ≤i

为方便后续计算，可事先统计出国境线多边形点序列χ1的经度及纬度的值域(即最大值和最小值)，设国境线多边形经度取值范围为[α1min，α1max] ，纬度取值范围为[β1min，β1max] 。

2.2 覆盖区与国境线关系粗判断

对某个时刻航天器地面覆盖区多边形的点序列χ0，统计其边界点经度及纬度取值范围(最大值和最小值)。 设统计出的地面覆盖区域边界点经度范围为[α0min，α0max] ，纬度范围为[β0min，β0max] 。

据此，可先对地面覆盖区域与国境线相交的可能性进行粗判断。 若β0min≥β1max，则无相交区域；若β0max≤β1min，则无相交区域；若α0max≤α1min，则无相交区域；若α0min≥α1max，则无相交区域。 无相交区域时，航天器返回覆盖区进入国境线内部分的面积占比为0，面积也为0。

2.3 球面纬线分割

对航天器返回覆盖区的球面多边形进行纬线分割，将覆盖区纬度取值范围[β0min，β0max] 等分为n份(n为事先给定的等分份数，例如取1000)。 如图2 所示(图中多边形为一个覆盖区多边形示例)，这样就形成了n个环状纬度分割带。

图2 球面纬线等分示意图Fig.2 Diagram of dividing latitude line on sphere

设地面覆盖区域纬度取值范围中各条等分线上的纬度值依次为β0、β1…βn。 其中第i(0 ≤i

考虑到此处对纬度进行的是n等分，则任意相邻2 条纬度等分线间的纬度差是相同的，即有βi＋1－βi＝(β0max－β0min)/n。 同时，由球冠面积公式可推导出i和i＋1 条纬度线之间的整个环状纬度带区域的面积为Si＝2πR2×(sinβi＋1－sinβi) 。其中，R为地球半径。

图3 所示为当n足够大时，对某一个纬度环状带，2 条纬度线长度足够接近，纬度圈单位圆心张角上梯形部分的面积可近似为式(1)所示。

图3 纬度环状带示意图Fig.3 Schematic diagram of the latitude ring

即当n足够大时，若能统计出该纬度圈在某个区域内的圆心张角之和，即可快速估算出该区域内部分的面积近似值。

由此可见，对第i个纬度环状分割带，可通过统计其在覆盖区内部分的纬度圈圆心张角之和(记为κi)即可得到该部分的面积，同样通过统计在覆盖区与国境线共有部分的纬度圈圆心张角之和(记为μi)就可得到相交部分的面积。

通过对各个纬度分割带的累加即可得到整个地面覆盖区面积、覆盖区进入国境内部分的面积。据此可得出地面覆盖区面积为式(2)所示。

覆盖区在国境线内部分区域的面积为式(3)所示。

这样就可得到地面覆盖区进入国境内部分区域面积占比的快速计算公式为式(4)所示。

本文方法仅仅进行了纬线方向的分割，而没有进行经线方向分割，因而相对于传统的网格式分割算法，计算速度大幅度得到提升。

2.4 纬线扫描比较计算相交区域

对图3 所示的地面覆盖区域纬度范围进行n等份分割，依次扫描每一个纬度分割带。 统计第i和第i＋1 条纬度等分线的中间线在地面覆盖区内、覆盖区与国境线相交部分内的纬度线圆心张角和(即经度跨度)。

图4 给出了第i和第i＋1 条纬度等分线的中间线与球面多边形某部分边的相交关系示意图。

图4 纬线分割带与多边形边相交示意图Fig.4 Schematic diagram of the intersection of latitude dividing zone and the polygon

在纬度线圆心张角之和统计时，可从覆盖区及国境线球面多边形左侧向右扫描目标纬度线(中间纬度线)，计算相应部分的圆心张角之和。

由覆盖区经度最小值α0min、国境线经度最小值α1min，可得两者的经度最小值为α′min＝min(α0min，α1min) ，则此中间纬度线上左侧点的经度可取为α′min－τ(τ为事先给定的大于零的值，如取1.0，以保证该点在2 个多边形之外)，纬度值为β′i。

计算该左侧点向右发出的中间纬度线与覆盖区的交点坐标，因该左侧点在覆盖区及国境线多边形的外部，则中间纬度线与覆盖区多边形交点个数必为偶数(或无交点)[14-15]。 计算时可逆时针方向遍历覆盖区边界多边形的各条边，然后计算与各边的交点(有时2 条边的交点可能相同)，形成1 个交点序列。 设该中间纬度线与覆盖区多边形交点序列(记为序列χ2)为p20，p21，…，p2i，…(交点的总个数取决于在该纬度线上的边数)，计算每2 个交点(如p20与p21，p22与p23，…)之间的弧线对应的纬度圈圆心张角(即这2 个点的经度之差，因为纬度线上各点的纬度值是相同的)，然后对每2 个点的经度差进行累加，即得到κi(该纬度线在覆盖区内的弧段的纬度圈圆心张角之和)。

如图4 所示，某个中间纬度线与多边形部分边形成了A－F共6 个交点，图中红色弧线部分为被边切割后在多边形内的弧段。 将这几个累加就可得到κi。 同理，可计算出该向右的中间纬度线与国境线多边形的交点序列(记为序列χ3)，将各交点标记为p30，p31，…，p3i，… (交点的总个数取决于在该纬度线上的边数)。

下面需要计算该向右的纬度线在覆盖区与国境线共有部分的弧段的纬度圈圆心张角之和对应的纬度圈张角度数和μi。 具体方法为:通过计算序列χ2中每个交线弧段(向右纬度线与覆盖区每2 个交点间的弧段)与序列χ3中每个交线弧段(向右纬度线与国土边界线每2 个交点间的弧段)共有部分的弧线的纬度圈圆心张角后，进行累加即可。 如对χ2中弧段1[p20，p21]，设2 个点经度分别为α20、α21，与χ3中弧段2[p30，p31]，设经度分别为α30、α31，因各点的纬度相同，只需要对经度值作简单地大小比较，即可得到共有部分弧线的张角。 具体地，当弧段1 给定时，弧段2 与其有以下几种关系。 图5 所示即为2 个弧段可能的关系图，需要说明的是弧段1 与弧段2 实际在同一纬度，此处为了表述清晰对其在上下方向进行了平移。

1)若弧段2 左端点在弧段1 左侧，则根据其右端点的位置不同，两者可能的位置关系如图5(a)所示；

2)若弧段2 左端点在弧段1 内部(包含两端点)，则两者可能的位置关系如图5(b)所示；

3)若弧段2 左端点在弧段1 右侧，则两者可能的位置关系如图5(c)所示。

从图5 可见:

1)若α21≤α30，弧段1 与弧段2 无共有弧段；

2)若α20≥α31，弧段1 与弧段2 无共有弧段；

3)其他情况时，弧段1 与弧段2 的共有弧段的圆心张角为min(α21，α31)－max(α20，α30) 。 这样就计算出了这2 个弧段的共有部分的纬度线圆心张角。 通过累加这些圆心张角即可得到交集部分的纬度线圆心张角之和μi，再由公式(2)～(4)就可得到完整的纬线扫描算法。

3 仿真分析

为了验证本文算法的可行性，进行仿真试算。算例1:采用纬线扫描算法的公式(2)计算德国国境线(大陆部分，不含任何岛屿)内区域面积(图6 中连线绘制的曲线部分)，该国境线为极不规则的多边形。 在纬线扫描计算时，对纬度方向进行1000 等份分割，计算结果为35.32 万km2，这与公开的德国大陆部分面积(约34.93 万km2，数据来源:维基百科)相近，误差约为1.1%，说明算法是可行的。

图6 覆盖区与德国边境线相交示意图Fig.6 Schematic diagram of the intersection of landing footprint and German border

算例2:仿真给定一覆盖区多边形区域(图6中散点绘制的曲线部分)，计算其进入德国大陆边境线内区域面积以及该面积与整个覆盖区面积的比值。 同时采用网格分割算法[13]进行计算比较，网格分割算法方法将球面多边形在经度、纬度方向进行等间隔分割，并计算各网格的面积，然后累加出覆盖区及其与国境线交集部分面积。

计算所用计算机为赛扬双核2.2 GHz 的CPU，4GB 内存，32 位的Win7 操作系统。

因覆盖区及入境部分无准确的参考面积，故此处只对2 种算法计算情况进行对比，算法对比结果见表1。 表1 给出了各算例计算的覆盖区面积、入境部分面积占比(百分比)以及计算耗时。

表1 各算例计算结果Table 1 Calculation results of each example

本算例中，覆盖区边界由234 个点组成，德国大陆界线由1612 个点组成；本文纬线扫描法计算时对覆盖区进行1000 等份分割，网格算法在经纬度方向各进行1000 等份分割；

为了考察本文算法精度影响因素，算例2 基础上继续给定以下算例。

算例3:在算例2 基础上将覆盖区边界点数减少一半(实际使用117 个点)，其他条件保持不变。

算例4:在算例2 基础上将德国大陆界线点数减少一半(实际使用806 点)，其他条件保持不变。

算例5:在算例2 基础上对覆盖区进行900等份分割，其他条件保持不变，同时使用经、纬度方向各900 等份划分的网格算法。

算例6:在算例2 基础上对覆盖区进行500等份分割，其他条件保持不变，同时使用经、纬度方向各500 等份划分的网格算法。

算例7:在算例2 基础上对覆盖区进行100等份分割，其他条件保持不变，同时使用经、纬度方向各100 等份划分的网格算法。

算例8:在算例2 基础上对覆盖区进行50 等份分割，其他条件保持不变，同时使用经、纬度方向各50 等份划分的网格算法。

算例9:在算例2 基础上对覆盖区进行10 等份分割，其他条件保持不变，同时使用经、纬度方向各10 等份划分的网格算法。

从算例2 可见，当进行1000 等份分割时，本文算法耗时只有网格方法的0.22%，入境部分百分比约相差0.1%。

从算例2、3 可见，覆盖区边界点从234 减少到117 时，覆盖区面积发生变化，但入境面积占比变化不明显，这是因为此覆盖区图形相对较平滑，点数减少一半后形状变化不明显。

算例4 中德国大陆界线点从1612 减少到806，覆盖区点不变时，本文算法耗时约减少1/3，可见对于点数过多的边界线数据，适当减少点数可提高计算效率。

从算例5～9 可见，当覆盖区、边界点数保持不变时，通过改变纬度划分段数，可有效提高计算效率，但面积精度变化较大，段数太少时面积精度过低。

在上述算例中，本文算法耗时都在0.4 s 以下，耗时明显比网格分割算法少，满足实时性要求。 各算例中，相对于网格算法，本文算法计算的入境部分面积占比的误差均在0.2%以内，算法可行。

4 结论

1)本文通过将传统的经纬度二维方向划分算法降维处理，仅在纬线方向上进行分割，并扩展到球面上，形成了一种再入覆盖区相交面积及面积占比的快速计算方法，仿真结果表明本文方法可行，计算速度满足实时性要求。

2)本文算法在具体实现中进行了部分近似，需满足精度要求时纬线分割数目不能太少，下一步将重点在保证实时性前提下对提高算法精度进行研究。