一種航天器返回覆蓋區(qū)面積快速計(jì)算方法

淡 鵬， 何雨帆， 陳俊收

(1.宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 西安 710043； 2.西安衛(wèi)星測控中心， 西安 710043)

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 球面緯線分割

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)圖3 所示的地面覆蓋區(qū)域緯度范圍進(jìn)行n等份分割，依次掃描每一個(gè)緯度分割帶。 統(tǒng)計(jì)第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)計(jì)時(shí)，可從覆蓋區(qū)及國境線球面多邊形左側(cè)向右掃描目標(biāo)緯度線(中間緯度線)，計(jì)算相應(yīng)部分的圓心張角之和。

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

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

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

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

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

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

3)若弧段2 左端點(diǎn)在弧段1 右側(cè)，則兩者可能的位置關(guān)系如圖5(c)所示。

從圖5 可見:

1)若α21≤α30，弧段1 與弧段2 無共有弧段；

2)若α20≥α31，弧段1 與弧段2 無共有弧段；

3)其他情況時(shí)，弧段1 與弧段2 的共有弧段的圓心張角為min(α21，α31)－max(α20，α30) 。 這樣就計(jì)算出了這2 個(gè)弧段的共有部分的緯度線圓心張角。 通過累加這些圓心張角即可得到交集部分的緯度線圓心張角之和μi，再由公式(2)～(4)就可得到完整的緯線掃描算法。

3 仿真分析

為了驗(yàn)證本文算法的可行性，進(jìn)行仿真試算。算例1:采用緯線掃描算法的公式(2)計(jì)算德國國境線(大陸部分，不含任何島嶼)內(nèi)區(qū)域面積(圖6 中連線繪制的曲線部分)，該國境線為極不規(guī)則的多邊形。 在緯線掃描計(jì)算時(shí)，對(duì)緯度方向進(jìn)行1000 等份分割，計(jì)算結(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中散點(diǎn)繪制的曲線部分)，計(jì)算其進(jìn)入德國大陸邊境線內(nèi)區(qū)域面積以及該面積與整個(gè)覆蓋區(qū)面積的比值。 同時(shí)采用網(wǎng)格分割算法[13]進(jìn)行計(jì)算比較，網(wǎng)格分割算法方法將球面多邊形在經(jīng)度、緯度方向進(jìn)行等間隔分割，并計(jì)算各網(wǎng)格的面積，然后累加出覆蓋區(qū)及其與國境線交集部分面積。

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

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

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

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

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

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

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

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

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

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

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

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

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

從算例2、3 可見，覆蓋區(qū)邊界點(diǎn)從234 減少到117 時(shí)，覆蓋區(qū)面積發(fā)生變化，但入境面積占比變化不明顯，這是因?yàn)榇烁采w區(qū)圖形相對(duì)較平滑，點(diǎn)數(shù)減少一半后形狀變化不明顯。

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

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

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

4 結(jié)論

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

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