点在球面多边形内外的判定方法及应用*

黄 跃，陈 黎，卢 宇

(1.北京师范大学天文学系，北京 100875;2.中国科学院高能物理研究所，北京 100039)

南大西洋辐射异常区(SAA)［1］，是指位于南美洲东边的南大西洋的地磁异常区域，与相邻区域相比，SAA区域的磁场强度仅为正常区域磁场强度的一半左右，属于负磁异常区。这使得该区高能带电粒子的通量非常大，对过该区的卫星会造成严重危害。因此，很多卫星在经过SAA区时，都要关机保护重要仪器设备，延长卫星运行寿命。在卫星轨道已知的情况下，卫星的开关机计划必须知道卫星是否进入了SAA区。进入SAA区时，如果关机时间提前，将浪费卫星有效运行时间;如果关机时间延后，可能对卫星的仪器造成损坏。离开SAA区时，如果开机时间提前，可能对卫星的仪器造成损坏;如果开机时间延后，将浪费卫星有效运行时间。因此，精确判断卫星进出SAA区，对于仪器保护和观测计划制定是很重要的，这就是点定位问题 (P点是否位于多边形区域R中)。

在计算几何、计算机图形系统、地理信息系统、CAD/CAM等研究领域，判断点在多边形内外是一个基本问题。迄今为止，不少学者对判断平面内点在多边形内外进行了深入研究。目前常用的算法有射线法和角度累加法等［2－4］。射线法是通过在待测点P上任作一条射线L，并统计出射线穿过多边形R的次数，奇数在内，偶数在外(图1)。如果射线与多边形的某一条边共线或经过多边形的顶点时，会出现奇异，需要甄别和判定。

实际的SAA区域为球面多边形区域，用传统的平面内的方法得到的误差较大。本文根据球面的特点，对射线法加以改进，不再是从待测点P发出一条射线，而是从P沿着过该点的经圈向北延伸到北极点，形成一条大圆弧，统计该弧段穿过球面多边形的次数，以此判断点在球面多边形内外，并对可能出现的奇异情况进行处理，实现了对卫星进出SAA区的精确判断。

图1 射线法图示Fig.1 Illustration of the“Ray Method”

1 算法框架

输入:球面多边形R(P1，P2，…，Pn－1，Pn)，点P;输出:点在多边形内、外或边界上的判定结果。为方便，分别记点Q点的经度和纬度为L(Q)和A(Q)。

判定待测点P在球面多边形R(P1，P2，…，Pn－1，Pn)内外的算法框架如下:

(1)顶点过滤

待测点P与多边形顶点的重合与否会影响后续的判断逻辑，需要首先予以甄别。如待测点P就是球面多边形的某个顶点Pi，则直接认定点属于球面多边形区域，判定结束。

(2)经纬度过滤

约定球面纬度范围为(－90°，90°)，经度范围为(－180°，180°)。经纬度过滤是通过对球面区域R和待测点的相对位置的判断，过滤掉不必考虑的边界。如图2，虚线部分的边界都是属于可以过滤掉的边界。

由于R跨越东经180°经线时，会导致数据不连续，这时需作相应的算法调整。判断球面多边形的某条边是否跨过东经180°的方法是:逐条判断球面多边形的每条边PiPi+1，如果:

图2 经纬度过滤Fig.2 Filtering based on latitude-longitude ranges

则确认该边过东经180°经线。此时须对Pi、Pi+1的经度进行处理:把经度小于0°的加上360°，使两个顶点的经度变成连续值。如果 L(P)＜0°，也要将P的经度加上360°。

现在考虑经度过滤问题，逐条判断球面多边形的每条边，如果待测点P的经度值不在边的两个顶点经度之间，则该边不予考虑。如果所有边界都被过滤掉，则说明待测点P不属于球面多边形R。

当R处于高纬度时，球面多边形的形变会导致误判，所以不做纬度过滤。

(3)交点数计算

经过经度过滤，只剩下图示情况的边需要考虑。图3中P0为北极点，P为待测点，PiPi+1，PjPj+1为球面多边形的边。PP0与边的交点经度就是P点的经度，只要计算交点纬度值。如果交点纬度值等于P点纬度值，则直接认定点P在球面多边形上，判定完成;如果交点纬度值大于P点纬度值，则交点在边上，该边与PP0相交，这时交点数增加1。

所有的边都完成判定后，计算交点的数量。如果是奇数，则点P在球面多边形内;如果是偶数，则点P在球面多边形外。

(4)判断球面多边形区域是否包含北极点P0

将球面多边形和北极点都投影到赤道面上，看北极的投影点是否在投影区域中。如果不在，则球面多边形不包含北极点，之前得到的点在球面多边形内外的结果不变。如果在，还需要再做一个处理，截取在北半球部分的球面多边形区域，用上面的方法判断北极点是否在新的球面多边形区域中，如果在，之前得到的点在球面多边形内外的结果相反。

图3 交点纬度计算Fig.3 Illustration of the calculation of the latitude of the intersection point

(5)奇异情况处理，待测点P的经度与球面多边形某个顶点经度相同，分3种情况，如图4。

a.L(P)=L(Pi)=L(Pi+1)

i.A(P)＞max{ A(Pi)，A(Pi+1)}时，不计算交点

ii.A(P)＜min{ A(Pi)，A(Pi+1)}时，交点数加1

iii.待测点P的纬度在上述两个顶点的纬度之间时，直接认为待测点在球面多边形上。

b.与待测点P经度相同的顶点两侧的边在待测点的两侧。这种情况交点数加1。

c.与待测点P经度相同的顶点两侧的边在待测点的一侧。这种情况交点数不变。

图4 奇异情况Fig.4 Singular cases

2 算法应用

这里选取硬X射线调制望远镜(HXMT)［5］做为考察对象。HXMT是基于直接解调方法，采用简单可靠的成熟技术、非位置灵敏探测器建造的高灵敏度硬X射线望远镜，致力于实现硬X射线波段(20～250 keV)的高灵敏度、高分辨率的巡天和定点观测，对感兴趣天体进行长时标能谱分析和时变分析，研究致密天体、黑洞强引力场中的动力学和高能辐射过程，研究早期宇宙的高能过程。探测器在SAA区要关机保护。SAA区以一组顶点形式给定。用上面的算法与传统的平面内的算法相比较。

当HXMT位置采样间隔为1 min时，两种算法基本能得到相同的结果。

图5 HXMT位置采样间隔为1 min(图中虚线为陆地边界，实线为SAA区边界。点为HXMT的轨迹。in:表示计算结果在SAA区内;out:表示计算结果在SAA区外)Fig.5 Application of the algorithm to the positions of the HXMT satellite.The sampling of the HXMT positions has a time interval of 1 min.The coastal lines are in dashed lines.The boundaries of the SAA region are in solid lines.The points are for the projected HXMT positions(“in:”within the SAA according to the algorithm; “out:”out of the SAA according to the algorithm)

但当采样间隔为1 s时，两者的差别就显现出来。图6分别为用本文提出的算法(左图)和用传统的平面内的方法(右图)得到的结果。

图6 HXMT位置采样间隔为1s(图中虚线为陆地边界，实线为SAA区边界。点为HXMT的轨迹。in:表示计算结果在SAA区内;out:表示计算结果在SAA区外)Fig.6 Application of the algorithm to the positions of the HXMT satellite.The sampling of the HXMT positions has a time interval of 1s.The coastal lines are in dashed lines.The boundaries of the SAA region are in solid lines.The points are for the projected HXMT positions(“in:”within the SAA according to the algorithm; “out:”out of the SAA according to the algorithm)

从图6可以看出用传统的平面内的方法得到的结果相对于本文的方法是不精确的，HXMT已经进入了SAA区，可用传统的平面内的方法得到的却是在SAA区外(out)的结果。考虑到HXMT的运行速度，误差达到23 s。统计一个月的结果，用传统的平面内的方法产生的误差，HXMT在SAA区外，得到的结果在SAA区内，累计达到1727 s，浪费了一定的有效运行时间;HXMT在SAA区内，得到的结果在SAA区外，累计达到2316 s，这可能对HXMT的仪器性能造成一定影响，从而影响观测数据的质量。

3 结语

本文提出的算法与传统的射线法相比，能准确地判断点与球面多边形区域的关系。同时算法的几何意义明显，便于理解和实现，不仅可用于简单多边形的判断，还能用于环状区域、自相交区域等复杂区域的判断。并通过应用实例证实本文的研究工作所具有的重要的实际应用价值，可提高一定的卫星有效运行时间和保护卫星仪器设备。本文提出的算法用Matlab实现，实际使用效果良好。

［1］谢伦，濮祖荫，焦维新，等.南大西洋异常区内辐射带高能质子辐射环境长期变化的研究［J］.中国科学(D 辑)，2005，35(7):658－663.

［2］Kai Hormann，Alexander Agathos.The Point in Polygon Problem for Arbitrary Polygons ［J］.Computational Geometry，2001，20(3):131－144.

［3］周培德.计算几何——算法设计与分析 ［M］.北京:清华大学出版社，2005.

［4］邹有建，肖龙鑫，陈鼎.判断某点是否在任意多边形内两种算法的比较 ［J］.地矿测绘，2009，25(3):28－30，36.Zou Youjian，Xiao Longxin，Chen Ding.A Contrast between Two Approaches to Find Whether the Point Being inside a Polygon ［J］.Surveying and Mapping of Geology and Mineral Resources，2009，25(3):28－30，36.

［5］Ti-Pei Li，We Mei，Zhang Shuangnan，et al.The Hard X-ray Modulation Telescope HXMT［C］//International Cosmic Ray Conference 28th.Universal Academy Press，2003:2775－2778.