淡 鵬， 何雨帆， 陳俊收

(1.宇航動力學(xué)國家重點實驗室， 西安 710043； 2.西安衛(wèi)星測控中心， 西安 710043)

1 引言

航天器升力式再入返回[1-2]地球過程中可通過主動控制傾側(cè)角和速度攻角來調(diào)節(jié)飛行過程中的升力和阻力，使其獲得較大的橫向和縱向機動能力[3]。 為了衡量航天器再入飛行的機動性能或?qū)ふ液线m的安控時機，常常需要實時計算飛行過程中的覆蓋范圍(例如飛行過程的可達區(qū)、濺落區(qū)等)及其相關(guān)幾何參數(shù)[4]。

航天器再入飛行的覆蓋范圍定義了航天器潛在的覆蓋區(qū)域[4-6]，可為航天器再入飛行任務(wù)規(guī)劃提供關(guān)鍵信息，同時可幫助航天器在正常的再入飛行和突發(fā)故障時變更或選擇潛在的著陸位置。 在覆蓋范圍分析時需要實時計算該覆蓋范圍與特定的區(qū)域(如期望的落點限制區(qū)域、國境線、海域等)相交部分的面積或面積占比等參數(shù)，以便對當前飛行能力及飛行的安全性進行實時評估。 此類計算實質(zhì)上可轉(zhuǎn)化為計算球面上不規(guī)則多邊形區(qū)域的面積以及多個多邊形相交部分求解等問題。 最常見的多是二維平面上的算法。 例如葛偉華等[7]基于邊界跟蹤方法實現(xiàn)了一種平面上的區(qū)域面積計算方法；齊澄宇等[8]在GeoST 全球參考網(wǎng)格基礎(chǔ)上，設(shè)計了一種可以擺脫經(jīng)緯度和積分運算的面積計算方法。 這些方法主要針對平面圖形，應(yīng)用到球面時誤差擴大。 對球面上的計算問題，王會然等[9]利用圖斑橢球和高斯投影方法減小了面積計算的誤差；Kumar 等[10]給出了一種采用近似展開公式計算橢球面上面積的方法。 但以上方法還稍顯復(fù)雜，編程實現(xiàn)有部分難度。

多邊形相交區(qū)域面積計算是該問題更深層次的計算，它需要首先找出不規(guī)則區(qū)域的相交部分，然后再計算出面積等參數(shù)。 針對該問題，朱磊等[11]提出了一種平面上相交面積計算的雙向鏈表算法；汪榮峰等[12]基于多邊形的交、差運算，將覆蓋多邊形進行分解和三角化，采用統(tǒng)計算法實現(xiàn)了一種基于多邊形布爾運算的相交面積計算方法。 但以上方法實現(xiàn)時還不夠簡便，對于覆蓋區(qū)面積占比計算時的實時性有待進一步提高。

考慮到航天器再入飛行過程中，為了及時尋找合適的安控時機，覆蓋區(qū)面積及占比計算常常需要較高的實時性和對各種不規(guī)則形狀的適應(yīng)能力，需要快速對航天器進行飛行狀態(tài)的判斷。 因而需要研究一些誤差較小、實現(xiàn)過程簡潔、實時性高、可靠性強，且支持復(fù)雜幾何邊界和易于工程應(yīng)用的計算方法。 本文提出的方法與常用的經(jīng)度、緯度二維方向網(wǎng)格劃分算法[13]不同，僅在緯度方向進行一維分割的區(qū)域面積及面積占比計算。 該方法直接建立在球面而非平面形狀之上，利用球冠面積公式和弧段張角的統(tǒng)計思路，實現(xiàn)覆蓋區(qū)相交部分面積及面積占比的快速估計。

2 覆蓋區(qū)相交面積及占比快速估計

在航天器返回覆蓋區(qū)與特定區(qū)域(此處以極不規(guī)則的國境線為例)相交部分的面積及面積占比快速估算中，傳統(tǒng)的網(wǎng)格劃分算法[13]因為要進行經(jīng)度、緯度二維方向上的劃分，使得分割塊數(shù)較多時，計算的復(fù)雜度成二次方增長，計算速度明顯下降，因而較難滿足實時性的要求。 本文對其進行改進，只在緯度方向上進行分割，使計算復(fù)雜度降低了1 個維數(shù)，相較于網(wǎng)格算法的計算速度有較大提升。

計算覆蓋區(qū)相交面積的算法流程如圖1所示。

圖1 覆蓋區(qū)相交面積計算流圖Fig.1 Flow chart for calculating the intersection area of landing footprint

2.1 界線數(shù)據(jù)的準備

航天器返回覆蓋區(qū)進入國境內(nèi)部分面積及面積占比計算中需要用到2 個球面多邊形，分別是地面覆蓋區(qū)域邊界多邊形和國境線多邊形。

設(shè)2 個多邊形分別用一系列逆時針順序排列的多個點表示。 其中，地面覆蓋區(qū)域邊界多邊形頂點個數(shù)為n0，國境線多邊形頂點個數(shù)為n1。 設(shè)地面覆蓋區(qū)域邊界多邊形頂點序列為p0，p1…pi…pn0( 0 ≤i

為方便后續(xù)計算，可事先統(tǒng)計出國境線多邊形點序列χ1的經(jīng)度及緯度的值域(即最大值和最小值)，設(shè)國境線多邊形經(jīng)度取值范圍為[α1min，α1max] ，緯度取值范圍為[β1min，β1max] 。

2.2 覆蓋區(qū)與國境線關(guān)系粗判斷

對某個時刻航天器地面覆蓋區(qū)多邊形的點序列χ0，統(tǒng)計其邊界點經(jīng)度及緯度取值范圍(最大值和最小值)。 設(shè)統(tǒng)計出的地面覆蓋區(qū)域邊界點經(jīng)度范圍為[α0min，α0max] ，緯度范圍為[β0min，β0max] 。

據(jù)此，可先對地面覆蓋區(qū)域與國境線相交的可能性進行粗判斷。 若β0min≥β1max，則無相交區(qū)域；若β0max≤β1min，則無相交區(qū)域；若α0max≤α1min，則無相交區(qū)域；若α0min≥α1max，則無相交區(qū)域。 無相交區(qū)域時，航天器返回覆蓋區(qū)進入國境線內(nèi)部分的面積占比為0，面積也為0。

2.3 球面緯線分割

對航天器返回覆蓋區(qū)的球面多邊形進行緯線分割，將覆蓋區(qū)緯度取值范圍[β0min，β0max] 等分為n份(n為事先給定的等分份數(shù)，例如取1000)。 如圖2 所示(圖中多邊形為一個覆蓋區(qū)多邊形示例)，這樣就形成了n個環(huán)狀緯度分割帶。

圖2 球面緯線等分示意圖Fig.2 Diagram of dividing latitude line on sphere

設(shè)地面覆蓋區(qū)域緯度取值范圍中各條等分線上的緯度值依次為β0、β1…βn。 其中第i(0 ≤i

考慮到此處對緯度進行的是n等分，則任意相鄰2 條緯度等分線間的緯度差是相同的，即有βi＋1－βi＝(β0max－β0min)/n。 同時，由球冠面積公式可推導(dǎo)出i和i＋1 條緯度線之間的整個環(huán)狀緯度帶區(qū)域的面積為Si＝2πR2×(sinβi＋1－sinβi) 。其中，R為地球半徑。

圖3 所示為當n足夠大時，對某一個緯度環(huán)狀帶，2 條緯度線長度足夠接近，緯度圈單位圓心張角上梯形部分的面積可近似為式(1)所示。

圖3 緯度環(huán)狀帶示意圖Fig.3 Schematic diagram of the latitude ring

即當n足夠大時，若能統(tǒng)計出該緯度圈在某個區(qū)域內(nèi)的圓心張角之和，即可快速估算出該區(qū)域內(nèi)部分的面積近似值。

由此可見，對第i個緯度環(huán)狀分割帶，可通過統(tǒng)計其在覆蓋區(qū)內(nèi)部分的緯度圈圓心張角之和(記為κi)即可得到該部分的面積，同樣通過統(tǒng)計在覆蓋區(qū)與國境線共有部分的緯度圈圓心張角之和(記為μi)就可得到相交部分的面積。

通過對各個緯度分割帶的累加即可得到整個地面覆蓋區(qū)面積、覆蓋區(qū)進入國境內(nèi)部分的面積。據(jù)此可得出地面覆蓋區(qū)面積為式(2)所示。

覆蓋區(qū)在國境線內(nèi)部分區(qū)域的面積為式(3)所示。

這樣就可得到地面覆蓋區(qū)進入國境內(nèi)部分區(qū)域面積占比的快速計算公式為式(4)所示。

本文方法僅僅進行了緯線方向的分割，而沒有進行經(jīng)線方向分割，因而相對于傳統(tǒng)的網(wǎng)格式分割算法，計算速度大幅度得到提升。

2.4 緯線掃描比較計算相交區(qū)域

對圖3 所示的地面覆蓋區(qū)域緯度范圍進行n等份分割，依次掃描每一個緯度分割帶。 統(tǒng)計第i和第i＋1 條緯度等分線的中間線在地面覆蓋區(qū)內(nèi)、覆蓋區(qū)與國境線相交部分內(nèi)的緯度線圓心張角和(即經(jīng)度跨度)。

圖4 給出了第i和第i＋1 條緯度等分線的中間線與球面多邊形某部分邊的相交關(guān)系示意圖。

圖4 緯線分割帶與多邊形邊相交示意圖Fig.4 Schematic diagram of the intersection of latitude dividing zone and the polygon

在緯度線圓心張角之和統(tǒng)計時，可從覆蓋區(qū)及國境線球面多邊形左側(cè)向右掃描目標緯度線(中間緯度線)，計算相應(yīng)部分的圓心張角之和。

由覆蓋區(qū)經(jīng)度最小值α0min、國境線經(jīng)度最小值α1min，可得兩者的經(jīng)度最小值為α′min＝min(α0min，α1min) ，則此中間緯度線上左側(cè)點的經(jīng)度可取為α′min－τ(τ為事先給定的大于零的值，如取1.0，以保證該點在2 個多邊形之外)，緯度值為β′i。

計算該左側(cè)點向右發(fā)出的中間緯度線與覆蓋區(qū)的交點坐標，因該左側(cè)點在覆蓋區(qū)及國境線多邊形的外部，則中間緯度線與覆蓋區(qū)多邊形交點個數(shù)必為偶數(shù)(或無交點)[14-15]。 計算時可逆時針方向遍歷覆蓋區(qū)邊界多邊形的各條邊，然后計算與各邊的交點(有時2 條邊的交點可能相同)，形成1 個交點序列。 設(shè)該中間緯度線與覆蓋區(qū)多邊形交點序列(記為序列χ2)為p20，p21，…，p2i，…(交點的總個數(shù)取決于在該緯度線上的邊數(shù))，計算每2 個交點(如p20與p21，p22與p23，…)之間的弧線對應(yīng)的緯度圈圓心張角(即這2 個點的經(jīng)度之差，因為緯度線上各點的緯度值是相同的)，然后對每2 個點的經(jīng)度差進行累加，即得到κi(該緯度線在覆蓋區(qū)內(nèi)的弧段的緯度圈圓心張角之和)。

如圖4 所示，某個中間緯度線與多邊形部分邊形成了A－F共6 個交點，圖中紅色弧線部分為被邊切割后在多邊形內(nèi)的弧段。 將這幾個累加就可得到κi。 同理，可計算出該向右的中間緯度線與國境線多邊形的交點序列(記為序列χ3)，將各交點標記為p30，p31，…，p3i，… (交點的總個數(shù)取決于在該緯度線上的邊數(shù))。

下面需要計算該向右的緯度線在覆蓋區(qū)與國境線共有部分的弧段的緯度圈圓心張角之和對應(yīng)的緯度圈張角度數(shù)和μi。 具體方法為:通過計算序列χ2中每個交線弧段(向右緯度線與覆蓋區(qū)每2 個交點間的弧段)與序列χ3中每個交線弧段(向右緯度線與國土邊界線每2 個交點間的弧段)共有部分的弧線的緯度圈圓心張角后，進行累加即可。 如對χ2中弧段1[p20，p21]，設(shè)2 個點經(jīng)度分別為α20、α21，與χ3中弧段2[p30，p31]，設(shè)經(jīng)度分別為α30、α31，因各點的緯度相同，只需要對經(jīng)度值作簡單地大小比較，即可得到共有部分弧線的張角。 具體地，當弧段1 給定時，弧段2 與其有以下幾種關(guān)系。 圖5 所示即為2 個弧段可能的關(guān)系圖，需要說明的是弧段1 與弧段2 實際在同一緯度，此處為了表述清晰對其在上下方向進行了平移。

1)若弧段2 左端點在弧段1 左側(cè)，則根據(jù)其右端點的位置不同，兩者可能的位置關(guān)系如圖5(a)所示；

2)若弧段2 左端點在弧段1 內(nèi)部(包含兩端點)，則兩者可能的位置關(guān)系如圖5(b)所示；

3)若弧段2 左端點在弧段1 右側(cè)，則兩者可能的位置關(guān)系如圖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)計算德國國境線(大陸部分，不含任何島嶼)內(nèi)區(qū)域面積(圖6 中連線繪制的曲線部分)，該國境線為極不規(guī)則的多邊形。 在緯線掃描計算時，對緯度方向進行1000 等份分割，計算結(jié)果為35.32 萬km2，這與公開的德國大陸部分面積(約34.93 萬km2，數(shù)據(jù)來源:維基百科)相近，誤差約為1.1%，說明算法是可行的。

圖6 覆蓋區(qū)與德國邊境線相交示意圖Fig.6 Schematic diagram of the intersection of landing footprint and German border

算例2:仿真給定一覆蓋區(qū)多邊形區(qū)域(圖6中散點繪制的曲線部分)，計算其進入德國大陸邊境線內(nèi)區(qū)域面積以及該面積與整個覆蓋區(qū)面積的比值。 同時采用網(wǎng)格分割算法[13]進行計算比較，網(wǎng)格分割算法方法將球面多邊形在經(jīng)度、緯度方向進行等間隔分割，并計算各網(wǎng)格的面積，然后累加出覆蓋區(qū)及其與國境線交集部分面積。

計算所用計算機為賽揚雙核2.2 GHz 的CPU，4GB 內(nèi)存，32 位的Win7 操作系統(tǒng)。

因覆蓋區(qū)及入境部分無準確的參考面積，故此處只對2 種算法計算情況進行對比，算法對比結(jié)果見表1。 表1 給出了各算例計算的覆蓋區(qū)面積、入境部分面積占比(百分比)以及計算耗時。

表1 各算例計算結(jié)果Table 1 Calculation results of each example

本算例中，覆蓋區(qū)邊界由234 個點組成，德國大陸界線由1612 個點組成；本文緯線掃描法計算時對覆蓋區(qū)進行1000 等份分割，網(wǎng)格算法在經(jīng)緯度方向各進行1000 等份分割；

為了考察本文算法精度影響因素，算例2 基礎(chǔ)上繼續(xù)給定以下算例。

算例3:在算例2 基礎(chǔ)上將覆蓋區(qū)邊界點數(shù)減少一半(實際使用117 個點)，其他條件保持不變。

算例4:在算例2 基礎(chǔ)上將德國大陸界線點數(shù)減少一半(實際使用806 點)，其他條件保持不變。

算例5:在算例2 基礎(chǔ)上對覆蓋區(qū)進行900等份分割，其他條件保持不變，同時使用經(jīng)、緯度方向各900 等份劃分的網(wǎng)格算法。

算例6:在算例2 基礎(chǔ)上對覆蓋區(qū)進行500等份分割，其他條件保持不變，同時使用經(jīng)、緯度方向各500 等份劃分的網(wǎng)格算法。

算例7:在算例2 基礎(chǔ)上對覆蓋區(qū)進行100等份分割，其他條件保持不變，同時使用經(jīng)、緯度方向各100 等份劃分的網(wǎng)格算法。

算例8:在算例2 基礎(chǔ)上對覆蓋區(qū)進行50 等份分割，其他條件保持不變，同時使用經(jīng)、緯度方向各50 等份劃分的網(wǎng)格算法。

算例9:在算例2 基礎(chǔ)上對覆蓋區(qū)進行10 等份分割，其他條件保持不變，同時使用經(jīng)、緯度方向各10 等份劃分的網(wǎng)格算法。

從算例2 可見，當進行1000 等份分割時，本文算法耗時只有網(wǎng)格方法的0.22%，入境部分百分比約相差0.1%。

從算例2、3 可見，覆蓋區(qū)邊界點從234 減少到117 時，覆蓋區(qū)面積發(fā)生變化，但入境面積占比變化不明顯，這是因為此覆蓋區(qū)圖形相對較平滑，點數(shù)減少一半后形狀變化不明顯。

算例4 中德國大陸界線點從1612 減少到806，覆蓋區(qū)點不變時，本文算法耗時約減少1/3，可見對于點數(shù)過多的邊界線數(shù)據(jù)，適當減少點數(shù)可提高計算效率。

從算例5～9 可見，當覆蓋區(qū)、邊界點數(shù)保持不變時，通過改變緯度劃分段數(shù)，可有效提高計算效率，但面積精度變化較大，段數(shù)太少時面積精度過低。

在上述算例中，本文算法耗時都在0.4 s 以下，耗時明顯比網(wǎng)格分割算法少，滿足實時性要求。 各算例中，相對于網(wǎng)格算法，本文算法計算的入境部分面積占比的誤差均在0.2%以內(nèi)，算法可行。

4 結(jié)論

1)本文通過將傳統(tǒng)的經(jīng)緯度二維方向劃分算法降維處理，僅在緯線方向上進行分割，并擴展到球面上，形成了一種再入覆蓋區(qū)相交面積及面積占比的快速計算方法，仿真結(jié)果表明本文方法可行，計算速度滿足實時性要求。

2)本文算法在具體實現(xiàn)中進行了部分近似，需滿足精度要求時緯線分割數(shù)目不能太少，下一步將重點在保證實時性前提下對提高算法精度進行研究。