基于多邊形布爾運算的衛(wèi)星區(qū)域覆蓋分析算法

汪榮峰

(裝備學(xué)院 航天指揮系, 北京 101416)

?

基于多邊形布爾運算的衛(wèi)星區(qū)域覆蓋分析算法

汪榮峰

(裝備學(xué)院 航天指揮系, 北京 101416)

摘要傳統(tǒng)網(wǎng)格點法計算衛(wèi)星區(qū)域覆蓋性能的精度受網(wǎng)格大小影響且效率低，為此提出了一種基于多邊形布爾運算的新算法。算法在計算衛(wèi)星覆蓋帶與待分析區(qū)域相交多邊形的基礎(chǔ)上，基于多邊形交、差運算，將覆蓋多邊形分解為具有單一覆蓋屬性的組成部分；將分解結(jié)果三角化后利用球面三角形面積公式計算面積；最后統(tǒng)計面積以計算覆蓋率，并設(shè)計實現(xiàn)了分類統(tǒng)計方式和可視化表現(xiàn)方法。與傳統(tǒng)網(wǎng)格點法相比，該算法覆蓋率計算結(jié)果穩(wěn)定，不受類似網(wǎng)格大小之類因素影響，在接近精度情況下效率比網(wǎng)格點法提高約20倍。

關(guān)鍵詞衛(wèi)星；區(qū)域覆蓋；布爾運算；三角化；覆蓋率

衛(wèi)星區(qū)域覆蓋分析對于衛(wèi)星任務(wù)規(guī)劃[1]、星座設(shè)計[2]等都具有非常重要的意義，其主要任務(wù)是計算單顆或多顆衛(wèi)星(星座)在給定時間范圍對待分析區(qū)域的覆蓋率、覆蓋次數(shù)、總覆蓋時間、平均覆蓋時間、最大覆蓋間隔、平均覆蓋間隔等指標(biāo)[3]，其中覆蓋率是核心指標(biāo)。

衛(wèi)星區(qū)域覆蓋分析方法主要有解析法、網(wǎng)格點法和基于幾何運算的方法。解析法[4]是基于衛(wèi)星與地球的幾何關(guān)系直接得到覆蓋面積計算的解析公式，這種方法只適于單顆衛(wèi)星覆蓋性能分析，且待分析區(qū)域必須包含衛(wèi)星覆蓋范圍。

目前最常用的是網(wǎng)格點法[5-6]，待分析區(qū)域劃分為一系列網(wǎng)格(可按經(jīng)緯度、距離和面積劃分網(wǎng)格點)，對于每顆衛(wèi)星按一定步長計算其覆蓋網(wǎng)格，根據(jù)覆蓋網(wǎng)格數(shù)與總網(wǎng)格數(shù)的關(guān)系得到覆蓋率等指標(biāo)。這種方法易于實現(xiàn)、應(yīng)用廣泛，且可以避免重合覆蓋區(qū)域的多次統(tǒng)計，但計算量大、重復(fù)計算多、計算結(jié)果受網(wǎng)格大小影響，基于網(wǎng)格交點和按一定步長的計算方式都有可能導(dǎo)致誤判。秦睿杰等[7]提出了網(wǎng)格中抽樣的快速算法，但精度會受到一定影響。

白萌等[8]將衛(wèi)星瞬時覆蓋投影到二維平面上形成覆蓋多邊形，通過多邊形并運算獲得總的覆蓋多邊形，從而得到總瞬時覆蓋率。該方法效率高，且計算得到的面積較為精確，但只適于瞬時覆蓋分析，且只能得到總覆蓋率，無法得到覆蓋次數(shù)等其他信息。

本文算法通過計算出衛(wèi)星在經(jīng)緯度平面上覆蓋帶與待分析區(qū)域的相交多邊形，一方面減少了網(wǎng)格點法中各步長覆蓋之間的重復(fù)計算，另一方面消除了步長取值過大所導(dǎo)致的漏判。通過多邊形交、差運算將區(qū)域覆蓋多邊形劃分為具有單一覆蓋屬性的子多邊形，算法不但可用于瞬時覆蓋分析，也可用于一段時間的覆蓋分析；不但可計算總覆蓋率，也可計算各衛(wèi)星覆蓋率、覆蓋次數(shù)等其他關(guān)鍵指標(biāo)。算法的時間復(fù)雜度與分解后多邊形數(shù)有關(guān)，而不像網(wǎng)格點法取決于劃分的網(wǎng)格數(shù)目，效率更高。在覆蓋面積和覆蓋率的計算上，采用多邊形三角剖分之后運用球面三角形面積公式計算的方法，計算結(jié)果穩(wěn)定，不像網(wǎng)格法會受到網(wǎng)格大小之類因素影響。

1待分析區(qū)域內(nèi)衛(wèi)星覆蓋多邊形的計算

待分析區(qū)域往往以經(jīng)緯度坐標(biāo)給定，在經(jīng)緯度平面上既可能是規(guī)則矩形，也可能是復(fù)雜多邊形，如某個國家或地區(qū)。衛(wèi)星飛行過程中只有有限時間經(jīng)過該區(qū)域，首先計算這段時間形成的覆蓋多邊形。

衛(wèi)星傳感器覆蓋其星下點周圍一定范圍，沿衛(wèi)星飛行方向在地球表面形成覆蓋帶，如圖1a)所示。根據(jù)衛(wèi)星位置，計算覆蓋帶兩側(cè)點的經(jīng)緯度坐標(biāo)，如圖1a)中a、b、c、d點。根據(jù)衛(wèi)星位置和覆蓋角，構(gòu)造三維空間中射線，計算該射線與地球表面的交點，根據(jù)交點坐標(biāo)反算其經(jīng)緯度；衛(wèi)星軌道高、覆蓋角大時，射線可能與地球表面無交，此時構(gòu)造經(jīng)過該射線和地球中心的平面，以該平面與地球表面的切點作為覆蓋帶邊界點。得到覆蓋帶在經(jīng)緯度平面的投影，如圖1b)所示。

首先利用多邊形求交算法計算相交覆蓋帶，對于衛(wèi)星軌道上每個采樣點，算法為：

Step1利用該采樣點的2個覆蓋邊界點和下一采樣點的覆蓋邊界點構(gòu)造一四邊形；

Step2如果四邊形與待分析區(qū)域有交，轉(zhuǎn)Step3，否則轉(zhuǎn)Step5；

Step3如果沒有已形成的相交區(qū)域，構(gòu)造2個空的點表leftpts和rightpts；

Step4將邊界點分別輸出到leftpts和rightpts中，轉(zhuǎn)Step1；

Step5如果leftpts和rightpts不為空，合并leftpts和rightpts為多邊形并輸出，轉(zhuǎn)Step1。

以圖1b)為例，開始AabB與區(qū)域12345無交，繼續(xù)下一個；BbcC與12345有交，則將BC輸出到leftpts中，bc輸出到rightpts中；依次處理，分別將DEFG加入到leftpts，defg加入到rightpts中；GghH不再與12345有交，則將leftpts和rightpts合并為bcdefgGFEDCB，即為相交的覆蓋帶。

上述算法中，對軌道上第1個和最后1個采樣點，要根據(jù)傳感器類型進(jìn)行特殊處理。如為圓錐形傳感器，則需在對應(yīng)方向加入半圓(離散為多邊形)；如為四棱錐形傳感器，則需在對應(yīng)方向擴(kuò)展出2個點。

a) 衛(wèi)星覆蓋帶

b) 衛(wèi)星覆蓋經(jīng)緯度投影與待分析區(qū)域圖1　與待分析區(qū)域有交的衛(wèi)星覆蓋帶計算

上述算法僅得到與待分析區(qū)域有交的覆蓋帶，還非實際的區(qū)域內(nèi)覆蓋，如圖2a)所示，還須將得到的多邊形與待分析區(qū)域求交，得到區(qū)域內(nèi)覆蓋多邊形，如圖2b)所示。如果待分析區(qū)域跨越180°經(jīng)度，將其以180°經(jīng)度為界分解為2個或多個區(qū)域，然后再應(yīng)用本文算法求解。

a) 相交覆蓋帶

b) 區(qū)域覆蓋多邊形圖2　區(qū)域內(nèi)覆蓋多邊形計算

2區(qū)域覆蓋多邊形分解與面積計算

經(jīng)過第1節(jié)處理，得到的相交多邊形并不能直接計算覆蓋率，因為存在區(qū)域被多顆衛(wèi)星覆蓋，或者被同一衛(wèi)星多次覆蓋等情況，即多邊形存在復(fù)雜的相交情況。如果僅僅計算總覆蓋率1個指標(biāo)，可以通過計算所有多邊形的并,然后計算其面積來實現(xiàn)。為了支持其他指標(biāo)的計算，將這些相交多邊形分解為獨立的具有單一覆蓋屬性的小多邊形，即分解后的每個小多邊形覆蓋次數(shù)、覆蓋衛(wèi)星等情況一致。

如圖3所示，ABCD和abcd為一顆衛(wèi)星的覆蓋，EFGH為另一顆衛(wèi)星的覆蓋，三者之間都存在相交關(guān)系，需分解為單一覆蓋性質(zhì)的小多邊形。AE21、1bB、5cF、a47d為僅被第1顆衛(wèi)星覆蓋的區(qū)域；DH43、CG76為僅被第2顆衛(wèi)星覆蓋的區(qū)域；ED32、FC65為被第1顆和第2顆衛(wèi)星同時覆蓋的部分；125cB為被第1顆衛(wèi)星2次覆蓋的區(qū)域；2365則是被第1顆衛(wèi)星覆蓋2次、第2顆衛(wèi)星覆蓋1次的區(qū)域。

圖3　相交多邊形分解示意圖

基于多邊形交、差運算進(jìn)行覆蓋多邊形分解，任意2多邊形之間關(guān)系分為4種：

1) 2多邊形無交，不必特殊處理，當(dāng)多邊形與所有其他多邊形都無交的時候，該多邊形已是分解最終結(jié)果，輸出。

2) 2多邊形有交，且交與2多邊形均不重合，如圖4a)所示，ABCD和abcd2個多邊形。首先計算得到交1234，其覆蓋屬性為2多邊形屬性的和；然后計算ABCD與abcd的差，得到A14D和BC32；最后計算abcd與ABCD的差，得到bc21和da43，差的覆蓋屬性仍維持原來的值。

3) 2多邊形的交為第1個多邊形，如圖4b)所示，多邊形AabD與ABCD求交，結(jié)果為AabD。此時首先將AabD的覆蓋屬性改為2多邊形覆蓋屬性之和；然后計算第2個多邊形與第1個多邊形之差，得到BCba，其覆蓋屬性維持不變。

4) 2多邊形的交為第2個多邊形，與第3)種情況類似處理。多邊形的覆蓋屬性定義為一個數(shù)組，每一項包括覆蓋的衛(wèi)星、分辨率等各種指標(biāo)，覆蓋屬性合并即為該數(shù)組的合并。

a) 交不等于輸入多邊形　　　　b) 交等于第1個多邊形圖4　2多邊形相交關(guān)系

通過一個先進(jìn)先出隊列來實現(xiàn)覆蓋多邊形分解，首先將第1節(jié)計算得到的覆蓋多邊形加入隊列，然后以隊列是否為空來控制循環(huán)，執(zhí)行如下邏輯：取出隊列中第1個多邊形；判斷該多邊形與隊列中所有其他多邊形的關(guān)系，如果均無交，則輸出；否則按上1段的描述生成新多邊形并加入到隊尾。

以上算法中，需使用多邊形交、差布爾運算，由于各多邊形相交情況復(fù)雜，布爾運算組合情況多，因此在多次運算后會出現(xiàn)多邊形有公共點、公共邊的情形，對算法穩(wěn)定性要求高，因此使用開源幾何引擎(Geometry Engine Open Source，GEOS)作為多邊形布爾運算工具。而在第1部分的相交覆蓋帶計算算法中，不存在上述奇異情況，但對衛(wèi)星軌道上每個采樣點都需應(yīng)用多邊形求交運算，效率要求高，使用效率更高的算法[9]。

不同于網(wǎng)格點法以網(wǎng)格數(shù)比值來表示覆蓋率，本文算法通過計算每個分解多邊形的面積來實現(xiàn)。采用的方法是：首先將多邊形三角化，然后將其反算到地球表面，使用球面三角形面積公式進(jìn)行計算。

使用OpenGL的GLU輔助庫中的GLUtesselator來實現(xiàn)多邊形三角化。GLUtesselator對多邊形剖分的結(jié)果是OpenGL中的基本圖元，而并非都是三角形。因此，通過其回調(diào)函數(shù)記錄剖分結(jié)果之后，再根據(jù)圖元類型(包括三角形扇、三角形帶、矩形、凸多邊形等)，將其轉(zhuǎn)化為三角形。

雖然也有橢球面三角形面積的計算算法[10]，但是其計算復(fù)雜、計算量大，對地球采用球形近似對精度影響很小，因此采用球面三角形面積計算公式[11]

(1)

式中：S為三角形球面面積；A、B、C分別為球面三角形的3個內(nèi)角(大圓弧在交點處切線的夾角)；R為地球半徑。內(nèi)角計算方法為：根據(jù)頂點經(jīng)緯度坐標(biāo)計算地心坐標(biāo)，過地心到頂點構(gòu)造單位矢量；每2個矢量計算矢量積得到大圓平面的法向量；通過單位法向量的數(shù)量積，反余弦計算得到內(nèi)角。

3精度效率分析及可視化表達(dá)方式設(shè)計

通過統(tǒng)計分解后的多邊形信息，得到區(qū)域覆蓋率等指標(biāo)。累加所有分解多邊形的面積，得到總覆蓋面積；將待分析區(qū)域三角化后，計算得到區(qū)域總面積；以總覆蓋面積除以總面積，得到覆蓋率。

首先對本文算法和網(wǎng)格點法的效率進(jìn)行簡單的理論分析。極端情況下，任意2個衛(wèi)星覆蓋條帶都有交，相交多邊形為n2量級，此時多邊形布爾運算次數(shù)為n3量級；對于網(wǎng)格點法，衛(wèi)星個數(shù)為n，網(wǎng)格劃分為m，采用判斷覆蓋條帶與網(wǎng)格關(guān)系的方法(原始網(wǎng)格點法是針對衛(wèi)星軌道上每個采樣點，對所有網(wǎng)格進(jìn)行判斷，非常耗時，但可以記錄每個網(wǎng)格的時間信息尤其是持續(xù)信息，這是本文算法和采用覆蓋條帶與網(wǎng)格相交關(guān)系方法無法支持的)，其運算量為nm。有2個條件決定了本文算法的效率優(yōu)勢：一是衛(wèi)星區(qū)域覆蓋分析一般針對星座或有限多顆衛(wèi)星進(jìn)行(一方面n值有限(20～30)，另一方面這些衛(wèi)星軌道具有一定的相關(guān)性，任意2個覆蓋帶都相交的極端情況出現(xiàn)的可能性非常小)；二是m的值遠(yuǎn)遠(yuǎn)大于n，如對于經(jīng)緯度范圍1°左右的區(qū)域，按1 km劃分網(wǎng)格，其m值約為10 000，即nm遠(yuǎn)大于n3。

對東經(jīng)110°～130°，北緯35°～45°之間的區(qū)域，分析12顆衛(wèi)星(兩行軌道根數(shù)略，覆蓋角均為5°)在20140501T080000時刻開始24 h的覆蓋情況，網(wǎng)格點法和本文算法得到的覆蓋率及用時如表1所示。

表1　網(wǎng)格點法和本文算法覆蓋率計算結(jié)果及用時對比

可以看出，網(wǎng)格點法計算精度受到網(wǎng)格大小的影響。在本算例中，網(wǎng)格大小約為1 km時，其精度才接近本文算法。在算法效率方面，本文算法用時與網(wǎng)格大小為10 km左右時接近，而與精度接近的1 km網(wǎng)格相比，用時不到網(wǎng)格點法的1/20。需要說明的是，本文算法對比所用的網(wǎng)格點法實現(xiàn)中，已經(jīng)采用了本文算法第1部分的技術(shù)對衛(wèi)星軌道采樣點進(jìn)行了快速排除，否則其效率更低。

除了覆蓋面積，總覆蓋率之外，還可以根據(jù)多邊形覆蓋屬性得到分類統(tǒng)計信息：按覆蓋次數(shù)統(tǒng)計，根據(jù)多邊形由幾個原始多邊形相交得到，可知其覆蓋次數(shù)，將相同覆蓋次數(shù)的所有多邊形面積相加，得到統(tǒng)計結(jié)果；按衛(wèi)星統(tǒng)計，首先遍歷所有多邊形，確定其覆蓋衛(wèi)星組合，然后就每種組合累加相應(yīng)多邊形面積。此外，可以支持按分辨率、頻段等各種統(tǒng)計方式。

設(shè)計實現(xiàn)了柱形圖和餅圖2種統(tǒng)計結(jié)果可視化表現(xiàn)形式，顯示總面積、覆蓋面積、圖例等信息，并以柱形圖和餅圖形式顯示各種衛(wèi)星組合、覆蓋次數(shù)的覆蓋面積與百分比等。

4結(jié) 束 語

與網(wǎng)格點法相比，本文提出算法在覆蓋面積和覆蓋率計算的精度和效率方面有明顯優(yōu)勢，存在的主要不足有3點：

1) 無法支持總覆蓋時間、平均覆蓋時間等與時間有關(guān)指標(biāo)的分析，這是由于網(wǎng)格點法可以將時間作為網(wǎng)格的屬性，但是本文算法使用的多邊形大小不一，無法做此處理，應(yīng)用中相應(yīng)指標(biāo)可以通過對地面目標(biāo)的覆蓋時間窗口分析技術(shù)來獲得；

2) 采用的面積計算方法，假定地球為球形，與應(yīng)用橢球模型或地理投影方式相比，存在可忽略的微小誤差；

3) 當(dāng)衛(wèi)星軌道傾角很大，地表覆蓋帶很寬，且待分析區(qū)域也在高緯度地區(qū)時，星下線轉(zhuǎn)折相對陡峭，算法第1步中，星下線兩側(cè)點位置關(guān)系可能錯亂，如圖1中Aa、Bb、Cc可能相交，此時生成的覆蓋多邊形有時不再是簡單多邊形，結(jié)果不再準(zhǔn)確。

參考文獻(xiàn)(References)

[1]王慧林,邱滌珊,黃小軍，等.面向區(qū)域覆蓋的電子偵察衛(wèi)星規(guī)劃方法研究[J].兵工學(xué)報,2011,32(11):1365-1372.

[2]曾德林.快速響應(yīng)小衛(wèi)星星座設(shè)計及覆蓋性能仿真分析[J].計算機(jī)仿真,2014,31(6):73-77.

[3]劉文,張育林.區(qū)域覆蓋低軌衛(wèi)星移動通信系統(tǒng)星座優(yōu)化設(shè)計[J].上海航天,2007(4):43-47.

[4]張潤.基于重訪周期的對地偵察小衛(wèi)星星座設(shè)計[D].西安:西安電子科技大學(xué),2011:23-32.

[5]馬吉康.通信衛(wèi)星組網(wǎng)仿真系統(tǒng)的設(shè)計與實現(xiàn)[D].北京: 北京郵電大學(xué),2008:18-23.

[6]沈欣.光學(xué)遙感衛(wèi)星軌道設(shè)計若干關(guān)鍵技術(shù)研究[D]. 武漢:武漢大學(xué), 2012:67-68.

[7]秦睿杰,戴光明,王茂才，等.一種計算星座區(qū)域覆蓋率的高效抽樣網(wǎng)格點法[J].計算機(jī)應(yīng)用研究,2015,32(4):65-68.

[8]白萌,李大林,陳夢云.衛(wèi)星對地覆蓋區(qū)域的融合算法研究[C]//中國空間科學(xué)學(xué)會.第二十三屆全國空間探測學(xué)術(shù)交流會論文集.廈門:中國空間科學(xué)學(xué)會,2010：341-346.

[9]汪榮峰,廖學(xué)軍.格網(wǎng)劃分的雙策略跟蹤多邊形裁剪算法[J].圖學(xué)學(xué)報,2012,33(6):45-49.

[10]施一民,朱紫陽.利用測地坐標(biāo)計算橢球面上凸多邊形面積的算法研究[J].同濟(jì)大學(xué)學(xué)報(自然科學(xué)版),2006,36(4):504-507.

[11]張楚賓.球面三角學(xué)[M].北京:人民教育出版社,1978:52-53.

(編輯：李江濤)

Analysis Algorithm for Satellite Regional Coverage Based on Polygonal Boolean Operation

WANG Rongfeng

(Department of Space Command, Equipment Academy, Beijing 101416, China)

AbstractTraditional mesh point calculation of satellite regional coverage is affected by the grid size in precision and inefficient, so the author proposes a new algorithm based on polygonal Boolean operation. With this algorithm, based on the calculation of polygons formed by the intersection of coverage area and the zone to be analyzed, by using polygonal intersection and subtraction operation, the author decomposes the coverage polygon into many components with single attribute of coverage; and then, the author triangulates the decomposition result and calculates the area with spherical triangle area formula; in the end, the author calculates the area to obtain the coverage and realizes the categorization-based statistic method and visualized presentation method through deliberate design. Comparing with traditional mesh point method, this algorithm can draw stable result, free of influence of factors like size of grid and its efficiency is approx. 20 times of that of mesh point method when similar accuracy .

Keywordssatellite; regional coverage; Boolean operation; triangulation; coverage ratio

文獻(xiàn)標(biāo)志碼A DOI10.3783/j.issn.2095-3828.2016.02.019

文章編號2095-3828(2016)02-0083-05

中圖分類號TP391.9

作者簡介汪榮峰(1973-),男，副教授，主要研究方向為空間態(tài)勢與分析。wrflly@163.com

收稿日期2015-07-21