再入覆蓋區域面積計算的一種球面矢量方法

淡 鵬，石 峰，王 丹，王 奧

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

0 引言

航天器升力式再入返回地球過程中，其制導和控制能力受到動力學系統的強擾動性、大氣模型參數不確定性等因素制約，因而再入滑翔技術已經開始在一些航天器返回過程中進行應用，這使得此類航天器能夠到達的地表位置區域呈現出較大的范圍。將航天器再入過程能夠到達的著陸或交班區域范圍稱作再入覆蓋區域。

隨著再入過程位置速度的變化，再入覆蓋區域也在不斷的變化中。為了評估實際飛行能力或必要時對航天器進行安控操作，常常需要計算再入覆蓋區域的面積，以及該區域進入某個指定區域內的面積。

航天器再入覆蓋區域一般可用一系列點圍成的不規則的多邊形來進行區域邊界描述，這樣覆蓋區面積實際上就轉化為這一多邊形區域面積。而覆蓋區進入指定區域內部分面積可轉化為兩個不規則多邊形相交部分的面積計算。張海濤等提出了一種不規則區域面積四等分割計算法，通過逐次將目標區域所在的矩形四等分割，直到滿足精度閾值要求為止，將區域內所有有效值相加得出區域總面積。葛偉華等提出了一種基于邊界跟蹤的區域面積計算方法，根據圖像邊界跟蹤時下一次和上一次跟蹤方向，確定圖像左右邊界，利用邊界像素的橫坐標進行加權求和計算，得到圖像區域面積。但這些方法主要針對平面圖形，應用到球面上時還需要一些改進。

而對于相交區域面積計算問題，魏許青利用Weiler-Atherton多邊形裁剪算法實現了平面多邊形的交集計算,但其也主要針對平面多邊形。

對于球面上的計算問題，施一民等給出了一種用三頂點測地坐標計算橢球面上三角形面積方法來計算凸多邊形面積；齊澄宇等提出了一種基于地球網格剖分的區域面積計算方法；劉洋等提出了一種基于拋物線擬合的橢球面區域面積計算方法。

文中從航天器再入覆蓋區域面積及進入特定區域部分面積計算的實際情況出發，基于一般矢量運算及面積矢量概念，結合球面三角形運算法則，實現了一種再入覆蓋區面積及相交面積的計算方法。

1 覆蓋區面積計算

1.1 面積計算方法

以逆時針順序排列的一系列經緯度點來表示再入覆蓋區的邊界。考慮到飛行器高度較高時，再入覆蓋區范圍可能較大，因而使用傳統的平面上區域面積計算方法將會存在較大偏差，故計算時需要考慮球面形狀。

在計算球面多邊形面積時，一種有效的方法就是將多邊形分割成多個三角形，然后利用球面三角形面積公式進行計算。

設如圖1所示的球面上,,三個點的經緯度坐標分別為[,]、[,]、[,]，地球半徑為，邊,,所對應的地心張角分別為,,。

圖1 球面三角形示意圖

由球面上兩點距離公式可計算出：

(1)

再由球面三角余弦公式可計算出點處三角形內角為:

=arccos[(cos-coscos)(sinsin)]

(2)

同樣可計算出點,處的三角形內角,。

則由球面三角形面積公式可得到三角形的面積為:

=(++-π)

(3)

將表示再入覆蓋區的多邊形分割成多個三角形，即可通過累加這些三角形面積得到可達區面積。但是由于再入覆蓋區是不規則的，且常常是凹多邊形，直接分割算法有時不便使用。

為此，借鑒面積矢量的定義，在三角形分割及面積計算時，將計算的面積定義為有向面積。

具體計算方法為:對如圖2所示的由個點組成的球面多邊形再入覆蓋區，將逆時針順序排列的各個頂點依次標記為,,,…,-1，以第0個點為基準點(為便于后續表達，同時將點標記為)，逆時針方向遍歷其它點，每兩個相鄰的點與基準點組成三角形，計算其有向面積，并進行累加。

圖2 球面多邊形三角劃分示意圖

設,,組成的球面三角的有向面積為;,,球面三角的有向面積為，依次類推，則整個球面多邊形面積為:

=|++…+0-(-2)-(-1)|

(4)

1.2 球面上面積正負判斷

使用式(4)計算時的重要一點是進行面積的正負符號判斷。根據有向面積的定義，當三角形第三個點出現在逆時針方向時為正，反之為負。

為了簡化順逆的判斷，一種思路是采用平面圖上矢量運算的判斷方法。以經度方向為方向，緯度方向為方向，建立各點經緯度坐標的平面圖。對平面圖中某個三角形Δ+1，則由矢量運算定律，在平面圖上，若點到矢量0叉乘點到+1的矢量0(+1)，所得結果大于0，則由右手法則，點+1在矢量0的逆時針方向，反之在其順時針方向。這樣即可判斷出各三角形的正負號。

但該方法對球面三角形有誤判的可能，為此改用球面上矢量關系運算。

如圖1所示，對點序列,,，判斷點在點序列,逆時針方向方法為:首先計算出逆時針遍歷的前兩點與球心所在平面的法向量,然后計算向量在該法向量上投影值,若其大于0，則在逆時針方向；若小于0，則在順時針方向，等于0時不構成球面三角形。

2 相交區域面積計算

對再入可達區評估計算中，需要計算可達區與一個指定的目標區域相交部分的面積，也就是兩個球面多邊形相交部分的面積計算。

2.1 相交部分計算

對球面上相交部分的計算，最簡單的是當兩個球面多邊形區域都是三角形的情況。

如圖3所示為兩個球面三角形，計算其相交區域，多邊形邊界點集可采用切割算法。

具體實現方法為：以逆時針點序遍歷第一個三角形每兩個相鄰頂點構成的連線，循環判斷第二個多邊形各個頂點位置，若該頂點在連線逆時針方向，則保留；若在順時針方向，則舍棄；若一個頂點與上一個頂點在不同方向，則計算當前連線與這兩個頂點連線的交點，保留該交點；對第二個多邊形遍歷完后，即可得到對用當前連線切割第二個多邊形一次后的剩余多邊形包絡點集。循環進行此切割,即可最終得到兩個三角形相交部分包絡點集。

圖3 三角形相交部分計算示意圖

切割算法也可擴充到兩個凸多邊形的情況，但是當多邊形為凹多邊形時就不再適用。為此，繼續使用有向面積算法，先將兩個多邊形以各自的一個基準點進行三角劃分(如圖4所示)；然后遍歷第一個多邊形各三角面片，分別計算其與第二個多邊形各三角面片的相交部分的面積，然后累加這些面積即可得到兩個多邊形相交部分面積。

圖4 球面多邊形相交部分計算示意圖

具體地，對第一個多邊形中一個三角面片Δ+1，首先判定該面片的正負號(符號記作Sign)，若為負值，則交換與+1坐標，形成臨時三角形1；若為非負數，則直接作為臨時三角形1。同樣對第二個多邊形中一個三角面片Δ+1(其符號記為Sign)，使用上述方法形成第二個臨時三角形2。然后計算臨時三角形1與2相交部分面積(記作，為正數)，則兩個三角面片相交部分有向面積為：

(Δ+1,Δ+1)=Sign·Sign·

(5)

進而得到兩個多邊形相交部分面積為：

(6)

式中:,分別為,兩個多邊形頂點數。

2.2 兩條弧線交點計算

設球面上一點的經緯度分別為,，則其在地球固連系下的直角坐標為:

·[coscos,sincos,sin]

(7)

式中為地球半徑。

則由點，，及球心可確定一個平面，設該平面方程為:

++=0

(8)

系數,,計算方法參見解析幾何相關內容，此處不再贅述。

由,及球心形成的平面方程為：

++=0

(9)

同樣，對半徑為的球體，有

++=

(10)

聯立式(7)～式(10)可計算出交點的直角坐標，進而計算出其經緯度值。

但該方法編程實現較為復雜，此處仍采用矢量運算的方法。

圖5 兩條弧線交點計算示意圖

此方法解出的點有兩個解。解選擇方法為：分別計算出兩個解離，,,四個點的球面距離和，取距離和最小的點作為所選擇的點。

2.3 僅一個相交區域時的簡化處理

上面的相交區域計算方法實現較為復雜。當覆蓋區與目標區域的相交部分僅為一個封閉區域，不可能有兩個或多個不相連接的封閉區域時，可以進一步作簡化處理。

如圖6所示，兩個多邊形只有一塊相交區域時，簡化處理方法為：逆時針順序遍歷第一個多邊形各邊，若判斷出該邊與第二個多邊形某邊有交點時，記錄下該交點，將其作為兩個多邊形交集的第一個點；然后繼續遍歷第一個多邊形的下一個邊，若與第二個沒交點，將該邊第一個頂點加入交集定點序列；若有交點，計算出交點加入交集序列；然后將第二個多邊形該交點后的所有頂點加入交集序列，直到到達第一個交點處時結束，將多邊形頂點序列作首尾相接處理，這樣即可得到相交部分多邊形頂點序列。進而通過此多邊形面積計算即可得到交集部分面積。

3 算法驗證

為了驗證該球面上覆蓋區面積算法的可行性，以多邊形面積計算中常用的網格分割算法為比較對象，將球面多邊形在經度、緯度方向進行等間隔分割并計算各網格的面積，進而累加出可達區及其與目標區域交集部分面積，然后與文中方法進行比較。

考慮到網格分割算法在邊界處可能存在面積多算的問題，故對結果的比較只要求近似即可。

計算時為了簡化處理，將地球半徑設置為1。對比如表1所示。由表中數據可見，兩種算法計算結果較為接近。

表1 計算結果對比(球半徑設為1)

4 結論

通過將平面上多邊形計算常用的有向面積概念擴充到球面上，給出了一種再入覆蓋區面積及其與目標區域交集部分面積的計算方法，計算結果表明該方法是可行的。

需要注意的是，因該算法使用了球面三角形面積公式，而球面三角形的邊是通過球心的大圓上一段弧線，但通過兩點的弧線可以存在很多條，因而使用時多邊形頂點個數不能過少，應能控制各邊走勢為宜。同時為了提高精度，應適當增加控制頂點的數量及減小間隔。

應該看到，此算法實現還是稍顯復雜，下一步將重點在算法的簡化處理及性能改進上進行研究。