再入覆蓋區(qū)域面積計(jì)算的一種球面矢量方法

淡 鵬，石 峰，王 丹，王 奧

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

0 引言

航天器升力式再入返回地球過程中，其制導(dǎo)和控制能力受到動(dòng)力學(xué)系統(tǒng)的強(qiáng)擾動(dòng)性、大氣模型參數(shù)不確定性等因素制約，因而再入滑翔技術(shù)已經(jīng)開始在一些航天器返回過程中進(jìn)行應(yīng)用，這使得此類航天器能夠到達(dá)的地表位置區(qū)域呈現(xiàn)出較大的范圍。將航天器再入過程能夠到達(dá)的著陸或交班區(qū)域范圍稱作再入覆蓋區(qū)域。

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

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

而對(duì)于相交區(qū)域面積計(jì)算問題，魏許青利用Weiler-Atherton多邊形裁剪算法實(shí)現(xiàn)了平面多邊形的交集計(jì)算,但其也主要針對(duì)平面多邊形。

對(duì)于球面上的計(jì)算問題，施一民等給出了一種用三頂點(diǎn)測(cè)地坐標(biāo)計(jì)算橢球面上三角形面積方法來計(jì)算凸多邊形面積；齊澄宇等提出了一種基于地球網(wǎng)格剖分的區(qū)域面積計(jì)算方法；劉洋等提出了一種基于拋物線擬合的橢球面區(qū)域面積計(jì)算方法。

文中從航天器再入覆蓋區(qū)域面積及進(jìn)入特定區(qū)域部分面積計(jì)算的實(shí)際情況出發(fā)，基于一般矢量運(yùn)算及面積矢量概念，結(jié)合球面三角形運(yùn)算法則，實(shí)現(xiàn)了一種再入覆蓋區(qū)面積及相交面積的計(jì)算方法。

1 覆蓋區(qū)面積計(jì)算

1.1 面積計(jì)算方法

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

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

設(shè)如圖1所示的球面上,,三個(gè)點(diǎn)的經(jīng)緯度坐標(biāo)分別為[,]、[,]、[,]，地球半徑為，邊,,所對(duì)應(yīng)的地心張角分別為,,。

圖1 球面三角形示意圖

由球面上兩點(diǎn)距離公式可計(jì)算出：

(1)

再由球面三角余弦公式可計(jì)算出點(diǎn)處三角形內(nèi)角為:

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

(2)

同樣可計(jì)算出點(diǎn),處的三角形內(nèi)角,。

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

=(++-π)

(3)

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

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

具體計(jì)算方法為:對(duì)如圖2所示的由個(gè)點(diǎn)組成的球面多邊形再入覆蓋區(qū)，將逆時(shí)針順序排列的各個(gè)頂點(diǎn)依次標(biāo)記為,,,…,-1，以第0個(gè)點(diǎn)為基準(zhǔn)點(diǎn)(為便于后續(xù)表達(dá)，同時(shí)將點(diǎn)標(biāo)記為)，逆時(shí)針方向遍歷其它點(diǎn)，每兩個(gè)相鄰的點(diǎn)與基準(zhǔn)點(diǎn)組成三角形，計(jì)算其有向面積，并進(jìn)行累加。

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

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

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

(4)

1.2 球面上面積正負(fù)判斷

使用式(4)計(jì)算時(shí)的重要一點(diǎn)是進(jìn)行面積的正負(fù)符號(hào)判斷。根據(jù)有向面積的定義，當(dāng)三角形第三個(gè)點(diǎn)出現(xiàn)在逆時(shí)針方向時(shí)為正，反之為負(fù)。

為了簡化順逆的判斷，一種思路是采用平面圖上矢量運(yùn)算的判斷方法。以經(jīng)度方向?yàn)榉较颍暥确较驗(yàn)榉较颍⒏鼽c(diǎn)經(jīng)緯度坐標(biāo)的平面圖。對(duì)平面圖中某個(gè)三角形Δ+1，則由矢量運(yùn)算定律，在平面圖上，若點(diǎn)到矢量0叉乘點(diǎn)到+1的矢量0(+1)，所得結(jié)果大于0，則由右手法則，點(diǎn)+1在矢量0的逆時(shí)針方向，反之在其順時(shí)針方向。這樣即可判斷出各三角形的正負(fù)號(hào)。

但該方法對(duì)球面三角形有誤判的可能，為此改用球面上矢量關(guān)系運(yùn)算。

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

2 相交區(qū)域面積計(jì)算

對(duì)再入可達(dá)區(qū)評(píng)估計(jì)算中，需要計(jì)算可達(dá)區(qū)與一個(gè)指定的目標(biāo)區(qū)域相交部分的面積，也就是兩個(gè)球面多邊形相交部分的面積計(jì)算。

2.1 相交部分計(jì)算

對(duì)球面上相交部分的計(jì)算，最簡單的是當(dāng)兩個(gè)球面多邊形區(qū)域都是三角形的情況。

如圖3所示為兩個(gè)球面三角形，計(jì)算其相交區(qū)域，多邊形邊界點(diǎn)集可采用切割算法。

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

圖3 三角形相交部分計(jì)算示意圖

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

圖4 球面多邊形相交部分計(jì)算示意圖

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

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

(5)

進(jìn)而得到兩個(gè)多邊形相交部分面積為：

(6)

式中:,分別為,兩個(gè)多邊形頂點(diǎn)數(shù)。

2.2 兩條弧線交點(diǎn)計(jì)算

設(shè)球面上一點(diǎn)的經(jīng)緯度分別為,，則其在地球固連系下的直角坐標(biāo)為:

·[coscos,sincos,sin]

(7)

式中為地球半徑。

則由點(diǎn)，，及球心可確定一個(gè)平面，設(shè)該平面方程為:

++=0

(8)

系數(shù),,計(jì)算方法參見解析幾何相關(guān)內(nèi)容，此處不再贅述。

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

++=0

(9)

同樣，對(duì)半徑為的球體，有

++=

(10)

聯(lián)立式(7)～式(10)可計(jì)算出交點(diǎn)的直角坐標(biāo)，進(jìn)而計(jì)算出其經(jīng)緯度值。

但該方法編程實(shí)現(xiàn)較為復(fù)雜，此處仍采用矢量運(yùn)算的方法。

圖5 兩條弧線交點(diǎn)計(jì)算示意圖

此方法解出的點(diǎn)有兩個(gè)解。解選擇方法為：分別計(jì)算出兩個(gè)解離，,,四個(gè)點(diǎn)的球面距離和，取距離和最小的點(diǎn)作為所選擇的點(diǎn)。

2.3 僅一個(gè)相交區(qū)域時(shí)的簡化處理

上面的相交區(qū)域計(jì)算方法實(shí)現(xiàn)較為復(fù)雜。當(dāng)覆蓋區(qū)與目標(biāo)區(qū)域的相交部分僅為一個(gè)封閉區(qū)域，不可能有兩個(gè)或多個(gè)不相連接的封閉區(qū)域時(shí)，可以進(jìn)一步作簡化處理。

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

3 算法驗(yàn)證

為了驗(yàn)證該球面上覆蓋區(qū)面積算法的可行性，以多邊形面積計(jì)算中常用的網(wǎng)格分割算法為比較對(duì)象，將球面多邊形在經(jīng)度、緯度方向進(jìn)行等間隔分割并計(jì)算各網(wǎng)格的面積，進(jìn)而累加出可達(dá)區(qū)及其與目標(biāo)區(qū)域交集部分面積，然后與文中方法進(jìn)行比較。

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

計(jì)算時(shí)為了簡化處理，將地球半徑設(shè)置為1。對(duì)比如表1所示。由表中數(shù)據(jù)可見，兩種算法計(jì)算結(jié)果較為接近。

表1 計(jì)算結(jié)果對(duì)比(球半徑設(shè)為1)

4 結(jié)論

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

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

應(yīng)該看到，此算法實(shí)現(xiàn)還是稍顯復(fù)雜，下一步將重點(diǎn)在算法的簡化處理及性能改進(jìn)上進(jìn)行研究。