點在球面多邊形內(nèi)外的判定方法及應(yīng)用*

黃 躍，陳 黎，盧 宇

(1.北京師范大學(xué)天文學(xué)系，北京 100875;2.中國科學(xué)院高能物理研究所，北京 100039)

南大西洋輻射異常區(qū)(SAA)［1］，是指位于南美洲東邊的南大西洋的地磁異常區(qū)域，與相鄰區(qū)域相比，SAA區(qū)域的磁場強度僅為正常區(qū)域磁場強度的一半左右，屬于負(fù)磁異常區(qū)。這使得該區(qū)高能帶電粒子的通量非常大，對過該區(qū)的衛(wèi)星會造成嚴(yán)重危害。因此，很多衛(wèi)星在經(jīng)過SAA區(qū)時，都要關(guān)機保護(hù)重要儀器設(shè)備，延長衛(wèi)星運行壽命。在衛(wèi)星軌道已知的情況下，衛(wèi)星的開關(guān)機計劃必須知道衛(wèi)星是否進(jìn)入了SAA區(qū)。進(jìn)入SAA區(qū)時，如果關(guān)機時間提前，將浪費衛(wèi)星有效運行時間;如果關(guān)機時間延后，可能對衛(wèi)星的儀器造成損壞。離開SAA區(qū)時，如果開機時間提前，可能對衛(wèi)星的儀器造成損壞;如果開機時間延后，將浪費衛(wèi)星有效運行時間。因此，精確判斷衛(wèi)星進(jìn)出SAA區(qū)，對于儀器保護(hù)和觀測計劃制定是很重要的，這就是點定位問題 (P點是否位于多邊形區(qū)域R中)。

在計算幾何、計算機圖形系統(tǒng)、地理信息系統(tǒng)、CAD/CAM等研究領(lǐng)域，判斷點在多邊形內(nèi)外是一個基本問題。迄今為止，不少學(xué)者對判斷平面內(nèi)點在多邊形內(nèi)外進(jìn)行了深入研究。目前常用的算法有射線法和角度累加法等［2－4］。射線法是通過在待測點P上任作一條射線L，并統(tǒng)計出射線穿過多邊形R的次數(shù)，奇數(shù)在內(nèi)，偶數(shù)在外(圖1)。如果射線與多邊形的某一條邊共線或經(jīng)過多邊形的頂點時，會出現(xiàn)奇異，需要甄別和判定。

實際的SAA區(qū)域為球面多邊形區(qū)域，用傳統(tǒng)的平面內(nèi)的方法得到的誤差較大。本文根據(jù)球面的特點，對射線法加以改進(jìn)，不再是從待測點P發(fā)出一條射線，而是從P沿著過該點的經(jīng)圈向北延伸到北極點，形成一條大圓弧，統(tǒng)計該弧段穿過球面多邊形的次數(shù)，以此判斷點在球面多邊形內(nèi)外，并對可能出現(xiàn)的奇異情況進(jìn)行處理，實現(xiàn)了對衛(wèi)星進(jìn)出SAA區(qū)的精確判斷。

圖1 射線法圖示Fig.1 Illustration of the“Ray Method”

1 算法框架

輸入:球面多邊形R(P1，P2，…，Pn－1，Pn)，點P;輸出:點在多邊形內(nèi)、外或邊界上的判定結(jié)果。為方便，分別記點Q點的經(jīng)度和緯度為L(Q)和A(Q)。

判定待測點P在球面多邊形R(P1，P2，…，Pn－1，Pn)內(nèi)外的算法框架如下:

(1)頂點過濾

待測點P與多邊形頂點的重合與否會影響后續(xù)的判斷邏輯，需要首先予以甄別。如待測點P就是球面多邊形的某個頂點Pi，則直接認(rèn)定點屬于球面多邊形區(qū)域，判定結(jié)束。

(2)經(jīng)緯度過濾

約定球面緯度范圍為(－90°，90°)，經(jīng)度范圍為(－180°，180°)。經(jīng)緯度過濾是通過對球面區(qū)域R和待測點的相對位置的判斷，過濾掉不必考慮的邊界。如圖2，虛線部分的邊界都是屬于可以過濾掉的邊界。

由于R跨越東經(jīng)180°經(jīng)線時，會導(dǎo)致數(shù)據(jù)不連續(xù)，這時需作相應(yīng)的算法調(diào)整。判斷球面多邊形的某條邊是否跨過東經(jīng)180°的方法是:逐條判斷球面多邊形的每條邊PiPi+1，如果:

圖2 經(jīng)緯度過濾Fig.2 Filtering based on latitude-longitude ranges

則確認(rèn)該邊過東經(jīng)180°經(jīng)線。此時須對Pi、Pi+1的經(jīng)度進(jìn)行處理:把經(jīng)度小于0°的加上360°，使兩個頂點的經(jīng)度變成連續(xù)值。如果 L(P)＜0°，也要將P的經(jīng)度加上360°。

現(xiàn)在考慮經(jīng)度過濾問題，逐條判斷球面多邊形的每條邊，如果待測點P的經(jīng)度值不在邊的兩個頂點經(jīng)度之間，則該邊不予考慮。如果所有邊界都被過濾掉，則說明待測點P不屬于球面多邊形R。

當(dāng)R處于高緯度時，球面多邊形的形變會導(dǎo)致誤判，所以不做緯度過濾。

(3)交點數(shù)計算

經(jīng)過經(jīng)度過濾，只剩下圖示情況的邊需要考慮。圖3中P0為北極點，P為待測點，PiPi+1，PjPj+1為球面多邊形的邊。PP0與邊的交點經(jīng)度就是P點的經(jīng)度，只要計算交點緯度值。如果交點緯度值等于P點緯度值，則直接認(rèn)定點P在球面多邊形上，判定完成;如果交點緯度值大于P點緯度值，則交點在邊上，該邊與PP0相交，這時交點數(shù)增加1。

所有的邊都完成判定后，計算交點的數(shù)量。如果是奇數(shù)，則點P在球面多邊形內(nèi);如果是偶數(shù)，則點P在球面多邊形外。

(4)判斷球面多邊形區(qū)域是否包含北極點P0

將球面多邊形和北極點都投影到赤道面上，看北極的投影點是否在投影區(qū)域中。如果不在，則球面多邊形不包含北極點，之前得到的點在球面多邊形內(nèi)外的結(jié)果不變。如果在，還需要再做一個處理，截取在北半球部分的球面多邊形區(qū)域，用上面的方法判斷北極點是否在新的球面多邊形區(qū)域中，如果在，之前得到的點在球面多邊形內(nèi)外的結(jié)果相反。

圖3 交點緯度計算Fig.3 Illustration of the calculation of the latitude of the intersection point

(5)奇異情況處理，待測點P的經(jīng)度與球面多邊形某個頂點經(jīng)度相同，分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)}時，交點數(shù)加1

iii.待測點P的緯度在上述兩個頂點的緯度之間時，直接認(rèn)為待測點在球面多邊形上。

b.與待測點P經(jīng)度相同的頂點兩側(cè)的邊在待測點的兩側(cè)。這種情況交點數(shù)加1。

c.與待測點P經(jīng)度相同的頂點兩側(cè)的邊在待測點的一側(cè)。這種情況交點數(shù)不變。

圖4 奇異情況Fig.4 Singular cases

2 算法應(yīng)用

這里選取硬X射線調(diào)制望遠(yuǎn)鏡(HXMT)［5］做為考察對象。HXMT是基于直接解調(diào)方法，采用簡單可靠的成熟技術(shù)、非位置靈敏探測器建造的高靈敏度硬X射線望遠(yuǎn)鏡，致力于實現(xiàn)硬X射線波段(20～250 keV)的高靈敏度、高分辨率的巡天和定點觀測，對感興趣天體進(jìn)行長時標(biāo)能譜分析和時變分析，研究致密天體、黑洞強引力場中的動力學(xué)和高能輻射過程，研究早期宇宙的高能過程。探測器在SAA區(qū)要關(guān)機保護(hù)。SAA區(qū)以一組頂點形式給定。用上面的算法與傳統(tǒng)的平面內(nèi)的算法相比較。

當(dāng)HXMT位置采樣間隔為1 min時，兩種算法基本能得到相同的結(jié)果。

圖5 HXMT位置采樣間隔為1 min(圖中虛線為陸地邊界，實線為SAA區(qū)邊界。點為HXMT的軌跡。in:表示計算結(jié)果在SAA區(qū)內(nèi);out:表示計算結(jié)果在SAA區(qū)外)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)

但當(dāng)采樣間隔為1 s時，兩者的差別就顯現(xiàn)出來。圖6分別為用本文提出的算法(左圖)和用傳統(tǒng)的平面內(nèi)的方法(右圖)得到的結(jié)果。

圖6 HXMT位置采樣間隔為1s(圖中虛線為陸地邊界，實線為SAA區(qū)邊界。點為HXMT的軌跡。in:表示計算結(jié)果在SAA區(qū)內(nèi);out:表示計算結(jié)果在SAA區(qū)外)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可以看出用傳統(tǒng)的平面內(nèi)的方法得到的結(jié)果相對于本文的方法是不精確的，HXMT已經(jīng)進(jìn)入了SAA區(qū)，可用傳統(tǒng)的平面內(nèi)的方法得到的卻是在SAA區(qū)外(out)的結(jié)果?？紤]到HXMT的運行速度，誤差達(dá)到23 s。統(tǒng)計一個月的結(jié)果，用傳統(tǒng)的平面內(nèi)的方法產(chǎn)生的誤差，HXMT在SAA區(qū)外，得到的結(jié)果在SAA區(qū)內(nèi)，累計達(dá)到1727 s，浪費了一定的有效運行時間;HXMT在SAA區(qū)內(nèi)，得到的結(jié)果在SAA區(qū)外，累計達(dá)到2316 s，這可能對HXMT的儀器性能造成一定影響，從而影響觀測數(shù)據(jù)的質(zhì)量。

3 結(jié)語

本文提出的算法與傳統(tǒng)的射線法相比，能準(zhǔn)確地判斷點與球面多邊形區(qū)域的關(guān)系。同時算法的幾何意義明顯，便于理解和實現(xiàn)，不僅可用于簡單多邊形的判斷，還能用于環(huán)狀區(qū)域、自相交區(qū)域等復(fù)雜區(qū)域的判斷。并通過應(yīng)用實例證實本文的研究工作所具有的重要的實際應(yīng)用價值，可提高一定的衛(wèi)星有效運行時間和保護(hù)衛(wèi)星儀器設(shè)備。本文提出的算法用Matlab實現(xiàn)，實際使用效果良好。

［1］謝倫，濮祖蔭，焦維新，等.南大西洋異常區(qū)內(nèi)輻射帶高能質(zhì)子輻射環(huán)境長期變化的研究［J］.中國科學(xué)(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］周培德.計算幾何——算法設(shè)計與分析 ［M］.北京:清華大學(xué)出版社，2005.

［4］鄒有建，肖龍鑫，陳鼎.判斷某點是否在任意多邊形內(nèi)兩種算法的比較 ［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.