一種復(fù)雜外形航天器投影面積計(jì)算的Monte Carlo方法

靳旭紅，黃　飛，程曉麗，丁延衛(wèi)，王　強(qiáng)

(1. 中國(guó)航天空氣動(dòng)力技術(shù)研究院， 北京 100074；2. 清華大學(xué)航天航空學(xué)院， 北京 100084; 3. 航天東方紅衛(wèi)星有限公司， 北京 100094)

0　引　言

近年來(lái)，由于重力場(chǎng)和穩(wěn)態(tài)海洋環(huán)流精確測(cè)量的需要，超低地球軌道(200～500 km)航天器逐漸展現(xiàn)出廣闊的應(yīng)用前景[1-3]，然而，實(shí)現(xiàn)重力梯度精確測(cè)量的先決條件是給低軌衛(wèi)星提供一套無(wú)拖拽姿態(tài)控制系統(tǒng) (Drag-Free and Altitude Control, DFAC)，該系統(tǒng)旨在采用一套離子發(fā)動(dòng)機(jī)和磁力矩器來(lái)補(bǔ)償大氣阻力等非重力的影響，從而進(jìn)行姿態(tài)和軌道控制[4]。DFAC正常工作的先決條件是低軌衛(wèi)星氣動(dòng)力，尤其是氣動(dòng)阻力的快速計(jì)算。對(duì)于大多數(shù)低軌衛(wèi)星外形，其氣動(dòng)阻力主要由當(dāng)?shù)卮髿饷芏取⒆枇ο禂?shù)和衛(wèi)星沿運(yùn)動(dòng)方向的投影面積共同決定。當(dāng)?shù)卮髿饷芏纫话阌纱髿饽Ｐ蚚5-8]給出，阻力系數(shù)的研究屬于稀薄氣體動(dòng)力學(xué)自由分子流理論的范疇，前人已經(jīng)進(jìn)行過(guò)大量的研究[9-13]：對(duì)于簡(jiǎn)單外形的低軌衛(wèi)星(衛(wèi)星各部件之間無(wú)相互遮擋)，一般采用面元積分(Panel Integral, PI)法[14]或者射線跟蹤面元積分 (Ray-Tracing Panel Integral, RTPI) 法[15]；更精確的阻力系數(shù)計(jì)算可以采用試驗(yàn)粒子Monte Carlo (Test Particle Monte Carlo, TPMC) 方法[16]，該方法能準(zhǔn)確模擬復(fù)雜外形部件之間的相互遮擋和多次反射效應(yīng)。因此，低軌航天器氣動(dòng)阻力預(yù)測(cè)需要解決計(jì)算其沿運(yùn)動(dòng)方向的投影面積 (Projected Cross-Sectional Area, PCSA)[17]的問(wèn)題。

2010年5月21日，日本航天局設(shè)計(jì)的太陽(yáng)帆探測(cè)器IKAROS (Interplanetary Kite-craft Accelerated by Radiation of the Sun) 由H-IIA運(yùn)載火箭發(fā)射升空[18]。太陽(yáng)帆探測(cè)器利用太陽(yáng)光壓力獲取推動(dòng)力，不需要攜帶大量推進(jìn)劑，雖然太陽(yáng)光壓很小，但持續(xù)的加速可使軌道器達(dá)到一個(gè)可觀的速度，最終可達(dá)到傳統(tǒng)航天器5～10倍的速度。太陽(yáng)帆探測(cè)器的姿態(tài)和軌道控制系統(tǒng)需要太陽(yáng)光壓力的輸入，太陽(yáng)光壓力則與太陽(yáng)通量、太陽(yáng)帆表面反射率和太陽(yáng)帆沿太陽(yáng)光線方向的投影面積有關(guān)。太陽(yáng)通量為物理常數(shù)，反射率為材料屬性，可通過(guò)查閱相關(guān)手冊(cè)或?qū)嶒?yàn)測(cè)量獲得，太陽(yáng)光壓力的計(jì)算就歸結(jié)為太陽(yáng)帆沿太陽(yáng)光線方向投影面積的計(jì)算[20]。

此外，航天器防護(hù)系統(tǒng)設(shè)計(jì)和風(fēng)險(xiǎn)評(píng)估中重要的一環(huán)是微流星體或空間碎片與航天器部件撞擊概率的預(yù)估[19]。進(jìn)行航天器易損性評(píng)估的前提是提供航天器及其各個(gè)部件在給定威脅攻擊方向上的投影面積以及存在重疊區(qū)域的重疊部分的面積，同時(shí)應(yīng)提供各個(gè)部件的位置遮擋關(guān)系[21]。2009年2月份發(fā)生的銥 (Iridium) 通信衛(wèi)星與一枚廢棄的俄羅斯衛(wèi)星碰撞事件說(shuō)明航天器防護(hù)系統(tǒng)設(shè)計(jì)和風(fēng)險(xiǎn)評(píng)估的重要性和挑戰(zhàn)性。同時(shí)，衛(wèi)星沿某一方向的投影面積是更直接姿態(tài)控制的必要輸入?yún)?shù)，比如太陽(yáng)瞄準(zhǔn) (Sun-pointing) 姿態(tài)模式需要輸入衛(wèi)星沿太陽(yáng)光線方向的投影面積，最小或最大阻力 (Minimum-drag or Maximum-drag) 姿態(tài)模式則需要識(shí)別并計(jì)算出衛(wèi)星的最小或最大投影面積[22]。

生成幾何實(shí)體沿某方向的投影，在計(jì)算機(jī)圖形學(xué)領(lǐng)域稱為遮擋算法，是計(jì)算幾何實(shí)體投影面積的一種直觀數(shù)值方法，至今仍在研究和完善之中。Appel[23]通過(guò)多邊形邊界的投影形成遮擋區(qū)域的邊界，發(fā)展了一種多邊形投影算法，應(yīng)用于由平面構(gòu)成的幾何實(shí)體投影圖像的生成。Atherton等[24]把上述多邊形投影算法加以推廣，以顯示多個(gè)光源下實(shí)體的投影圖像。目前，計(jì)算機(jī)圖形學(xué)廣泛使用的投影算法是光線投射和光線追蹤方法[25]：光線投射方法計(jì)算從光源發(fā)射出的光線與幾何實(shí)體表面的交點(diǎn)，通過(guò)最小距離識(shí)別出可見(jiàn)面(非可見(jiàn)面為背光面)，然后采用多邊形投影算法對(duì)可見(jiàn)面進(jìn)行投影即可生成幾何實(shí)體的投影圖像；光線追蹤方法在光線投射方法的基礎(chǔ)上，考慮了光線散射和反射效應(yīng)。

無(wú)論是低軌衛(wèi)星大氣阻力、太陽(yáng)帆探測(cè)器太陽(yáng)光壓的計(jì)算，還是航天器防護(hù)設(shè)計(jì)與風(fēng)險(xiǎn)評(píng)估，以及航天器姿態(tài)控制，都只需要航天器(部件)投影面積的輸入，通過(guò)計(jì)算機(jī)圖形學(xué)的方法生成投影圖像再計(jì)算投影面積顯得耗時(shí)耗力。最近，Ben等[26]基于凸多邊形理論，提出一種計(jì)算航天器沿任意方向投影面積的解析方法。該方法首先對(duì)航天器(部件)外表面進(jìn)行多邊形單元離散；然后分別計(jì)算每個(gè)多邊形單元的投影面積；再通過(guò)投影交點(diǎn)和內(nèi)部點(diǎn)識(shí)別和計(jì)算出投影多邊形相交部分的面積；最后從總的投影面積中減去相交部分的面積得到航天器的實(shí)際投影面積。由于每個(gè)過(guò)程均采用解析解的形式來(lái)實(shí)現(xiàn)，此方法的計(jì)算效率很高。然而，復(fù)雜外形的航天器進(jìn)行多邊形離散后，其投影多邊形難免存在兩兩相交甚至三個(gè)或三個(gè)以上多邊形同時(shí)相交的情形，導(dǎo)致識(shí)別計(jì)算投影多邊形相交部分的面積很困難，甚至不可能。解析的特性導(dǎo)致其魯棒性和通用性較差。極端情況下，計(jì)算機(jī)固有的舍入誤差會(huì)導(dǎo)致計(jì)算結(jié)果的錯(cuò)誤。航天器外形越復(fù)雜，出現(xiàn)極端情況的可能性越大。因此，上述解析方法雖然高效，卻很難應(yīng)用于航天工程實(shí)踐中。

借鑒統(tǒng)計(jì)學(xué)中Monte Carlo技術(shù)的思想，結(jié)合計(jì)算幾何學(xué)中基本的射線與三角形相交理論，提出一種計(jì)算復(fù)雜外形航天器沿任意方向投影面積的Monte Carlo方法。方法繼承了Monte Carlo類方法的優(yōu)點(diǎn)，即魯棒性和通用性好，計(jì)算時(shí)間隨著計(jì)算規(guī)模、問(wèn)題復(fù)雜程度和空間維數(shù)的增大呈線性關(guān)系。首先給出方法實(shí)現(xiàn)過(guò)程的詳細(xì)描述，并且以投影面積存在解析解的單面圓形平板驗(yàn)證了方法的正確性。然后，計(jì)算歐洲空間局SAMSON衛(wèi)星在典型姿態(tài)下的投影面積，與文獻(xiàn)報(bào)道的數(shù)據(jù)對(duì)比，進(jìn)一步驗(yàn)證該方法的可靠性以及對(duì)復(fù)雜工程外形的適用性。此外，采用共享存儲(chǔ)模式實(shí)現(xiàn)了該方法的并行化，以充分利用當(dāng)前高性能計(jì)算領(lǐng)域主流的多核處理器技術(shù)。最后，將該方法應(yīng)用于歐空局SAMSON小衛(wèi)星和GOCE低軌衛(wèi)星大量姿態(tài)下投影面積的計(jì)算，分析赤緯和赤經(jīng)對(duì)投影面積的影響，為衛(wèi)星的姿態(tài)和軌道控制提供輸入?yún)?shù)。

1　坐標(biāo)系統(tǒng)

1.1　坐標(biāo)系定義

本文使用以下2套坐標(biāo)系：機(jī)體坐標(biāo)系B和投影坐標(biāo)系P。機(jī)體坐標(biāo)系B的坐標(biāo)原點(diǎn)位于航天器質(zhì)心O，X、Y和Z軸與航天器慣性主軸重合。投影坐標(biāo)系P的原點(diǎn)同樣位于航天器質(zhì)心O，其z軸與投影方向重合，即投影方向?yàn)镺z軸，令Oz在平面的投影為OzP，則Oy位于XOY平面內(nèi)且垂直于OzP，Ox通過(guò)右手定則確定，如圖 1所示。

采用航天器軌道力學(xué)或天文學(xué)術(shù)語(yǔ)[27]，Oz與XOY平面的傾角定義為赤緯δ，其范圍為δ∈[-π/2,π/2]，OzP與OX的夾角定義為赤經(jīng)α，其范圍為α∈[0,2π)。令投影方向(即Oz軸)在機(jī)體坐標(biāo)系中歸一化的坐標(biāo)為(zX,zY,zZ)，可得

δ=arcsinzZ

(1)

(2)

1.2　坐標(biāo)轉(zhuǎn)換

根據(jù)坐標(biāo)系的定義，機(jī)體坐標(biāo)系OXYZ與投影坐標(biāo)系Oxyz之間的轉(zhuǎn)換關(guān)系為：

(3)

(4)

其中，(X,Y,Z)為某個(gè)點(diǎn)或矢量在機(jī)體坐標(biāo)系B中的坐標(biāo)，(x,y,z)為該點(diǎn)或矢量在投影坐標(biāo)系P中的坐標(biāo)。若非特殊說(shuō)明，下文中的坐標(biāo)或者矢量均為投影坐標(biāo)系P中的量。

圖1　機(jī)體坐標(biāo)系和投影坐標(biāo)系Fig.1　The body-axis and projected reference system

2　投影面積Monte Carlo模擬

Monte Carlo方法，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。該方法源于美國(guó)在第二次世界大戰(zhàn)期間研制原子彈的“曼哈頓計(jì)劃”。烏拉姆 (S. M. Ulam) 和該計(jì)劃的主持人之一、數(shù)學(xué)家馮·諾伊曼(J. von Neumann)用馳名世界的賭城—摩納哥 (Monaco) 的Monte Carlo—來(lái)命名這種方法。其實(shí)，Monte Carlo的基本思想最早可以追溯到18世紀(jì)初期的蒲豐 (Buffon) 投針問(wèn)題，法國(guó)科學(xué)家蒲豐用投針試驗(yàn)的方法來(lái)計(jì)算圓周率[28]。

采用Monte Carlo方法計(jì)算航天器沿任意方向投影面積的思想非常簡(jiǎn)單：在垂直于投影方向足夠大的矩形平面隨機(jī)產(chǎn)生大量試驗(yàn)粒子，所有試驗(yàn)粒子沿投影方向直線運(yùn)動(dòng)，計(jì)算出與航天器表面相交(碰撞)的比例，則投影面積即為該比例與上述矩形面積的乘積。下面給出詳細(xì)步驟。

2.1　包圍盒子的構(gòu)建

首先，構(gòu)建一個(gè)包圍航天器外表面的長(zhǎng)方體盒子 (Axis Aligned Bounding Box)，該長(zhǎng)方體的所有棱邊均平行于投影坐標(biāo)系P的3條坐標(biāo)軸。通常，復(fù)雜外形航天器表面由三角形單元離散，設(shè)表面三角形單元數(shù)量為N，單元節(jié)點(diǎn)在投影坐標(biāo)系中的坐標(biāo)為(xi,yi,zi),i∈[1,N]，則該長(zhǎng)方體盒子的6個(gè)面由下式確定

(5)

2.2　邊界矩形的確定

完成長(zhǎng)方體包圍盒子的構(gòu)建后，需要確定產(chǎn)生試驗(yàn)粒子的“邊界面”。邊界面垂直于投影方向，且位于包圍盒子的上游(相對(duì)于投影方向而言)，即邊界面上的坐標(biāo)z滿足

z

(6)

顯然，滿足上述條件的邊界面有無(wú)窮多個(gè)，不失一般性，本文取邊界面為

z=zmin-1

(7)

邊界矩形為包圍盒子沿投影反方向在邊界面上的投影，故其坐標(biāo)(x,y,z)滿足

(8)

2.3　試驗(yàn)粒子的產(chǎn)生

確定邊界矩形之后，即可產(chǎn)生試驗(yàn)粒子。試驗(yàn)粒子在邊界矩形隨機(jī)均勻分布，其位置由下式確定

(9)

式中，R1、R2均為(0,1)之間均勻分布的隨機(jī)數(shù)。由于投影方向?yàn)閛z方向，試驗(yàn)粒子的運(yùn)動(dòng)方向?yàn)?/p>

v=(0,0,1)

(10)

2.4　試驗(yàn)粒子的追蹤

試驗(yàn)粒子的跟蹤需要依次判斷試驗(yàn)粒子的軌跡是否與每一個(gè)三角形單元相交，涉及到計(jì)算幾何學(xué)中基本的射線與三角形的交點(diǎn)問(wèn)題。令試驗(yàn)粒子的坐標(biāo)為(x0,y0,z0)，三角形ABC三個(gè)頂點(diǎn)的坐標(biāo)分別為(x1,y1,z1)、(x2,y2,z2)和(x3,y3,z3)，表面外法向矢量為n=(n1,n2,n3)，試驗(yàn)粒子射線軌跡與三角形ABC所在平面的交點(diǎn)為P(x,y,z)，則P滿足射線的參數(shù)方程

(11)

式中，t為參數(shù)，同時(shí)，P滿足三角形ABC所在平面的方程，即

n·(x-x1,y-y1,z-z1)=0

(12)

若n3=0，則平面平行于射線，兩者顯然不相交，否則容易得

(13)

若t≤0，則射線與平面不相交，否則需要進(jìn)一步判斷交點(diǎn)是否在三角形ABC內(nèi)部。令

(14)

若同時(shí)滿足

(15)

則交點(diǎn)P在三角形ABC內(nèi)部，可知試驗(yàn)粒子的軌跡射線與該三角形單元相交；若式(15)不成立，即交點(diǎn)P在三角形ABC外部，則試驗(yàn)粒子的軌跡與該三角形單元不相交。

對(duì)于復(fù)雜外形而言，在某個(gè)試驗(yàn)粒子的跟蹤過(guò)程中，試驗(yàn)粒子的軌跡很可能與不止一個(gè)三角形單元相交，即出現(xiàn)了單元之間的相互遮擋現(xiàn)象，需要區(qū)別出“可見(jiàn)單元”和“被遮擋單元”。假定與某試驗(yàn)粒子軌跡相交的表面單元數(shù)為l，則對(duì)于每一個(gè)相交三角形單元，根據(jù)式(13)可以求出對(duì)應(yīng)的參數(shù)t，則最小參數(shù)值對(duì)應(yīng)的三角形單元為“可見(jiàn)”單元，其他(l-1)個(gè)單元為“被遮擋”單元。

通常，為了保證空間分辨率，復(fù)雜外形的航天器需要成千上萬(wàn)個(gè)三角形單元進(jìn)行離散，導(dǎo)致計(jì)算量比較大。為了減少計(jì)算時(shí)間，可以先排除“背光面”的表面三角形單元，即排除表面外法向向量沿z軸的分量為正的三角形單元(“背光”單元)，這可以降低將近一半的工作量而不影響精度。同時(shí)，對(duì)于表面三角形單元數(shù)量較多的復(fù)雜航天器外形，可以采用計(jì)算幾何學(xué)領(lǐng)域的加速搜索技術(shù)，比如交錯(cuò)數(shù)字樹(shù) (Alternating Digital Tree, ADT) 技術(shù)[29]，該算法可以將計(jì)算量從O(N)降低到O(lnN)。

2.5　投影面積的統(tǒng)計(jì)

設(shè)產(chǎn)生并跟蹤的試驗(yàn)粒子數(shù)為M，其中，運(yùn)動(dòng)軌跡與航天器表面相交的試驗(yàn)粒子數(shù)為Q，則航天器沿z軸的投影面積AP為

AP=(xmax-xmin)(ymax-ymin)(Q/M)

(16)

Monte Carlo方法屬于統(tǒng)計(jì)方法，難免存在統(tǒng)計(jì)誤差。統(tǒng)計(jì)誤差ε定義為多次計(jì)算結(jié)果誤差的均方根，根據(jù)統(tǒng)計(jì)學(xué)的中心極限定理[30]，其值為單一誤差分布(Single Error Distribution) 標(biāo)準(zhǔn)差σ與試驗(yàn)粒子數(shù)平方根的商，即

(17)

3　共享存儲(chǔ)并行模式實(shí)現(xiàn)

對(duì)于本文發(fā)展的計(jì)算投影面積的Monte Carlo方法，計(jì)算量大致與航天器表面網(wǎng)格單元數(shù)和試驗(yàn)粒子數(shù)的乘積成正比。一方面，為保證空間分辨率，復(fù)雜外形的航天器需要成千上萬(wàn)個(gè)三角形單元進(jìn)行離散；另一方面，為了減小統(tǒng)計(jì)誤差，需要采用上千萬(wàn)甚至幾十億個(gè)試驗(yàn)粒子：上述兩個(gè)方面都會(huì)導(dǎo)致計(jì)算量劇增。實(shí)際上，Monte Carlo方法雖然通用性和魯棒性好，能夠很好地解決“維數(shù)災(zāi)難”問(wèn)題，但是巨大的計(jì)算量一直是一個(gè)缺點(diǎn)。然而，隨著計(jì)算資源的穩(wěn)步增加、計(jì)算能力的逐步提高和并行算法的改善，這個(gè)缺點(diǎn)將會(huì)弱化。

采用共享存儲(chǔ)編程模式Open Multi-Processing (OpenMP) 實(shí)現(xiàn)投影面積Monte Carlo方法的并行化，這種并行模式能夠充分利用當(dāng)前主流的單節(jié)點(diǎn)多核 (Multiple-Core-Single-Node) 系統(tǒng)。隨著傳統(tǒng)單核處理器達(dá)到其運(yùn)行速度和集成復(fù)雜度的物理極限，多核處理器，即單節(jié)點(diǎn)多核系統(tǒng)成為實(shí)現(xiàn)高性能計(jì)算的有效途徑。比如，目前超級(jí)計(jì)算機(jī)世界500強(qiáng)中，有426個(gè)采用了四核處理器，59個(gè)采用了雙核處理器，只有4個(gè)仍在使用單核處理器[32]?？梢?jiàn)，多核處理器已經(jīng)成為高性能計(jì)算的一個(gè)趨勢(shì)。

對(duì)于投影面積Monte Carlo模擬的5個(gè)步驟中，包圍盒子的構(gòu)建和邊界矩形的確定計(jì)算量極小，沒(méi)有必要并行處理；試驗(yàn)粒子的產(chǎn)生和跟蹤最為耗時(shí)，尤其是試驗(yàn)粒子的跟蹤。因此，需要對(duì)這兩個(gè)步驟進(jìn)行并行處理。由于每個(gè)試驗(yàn)粒子之間沒(méi)有相互作用，這兩步的并行化只需要執(zhí)行以下步驟：

首先，把M個(gè)試驗(yàn)粒子大致平均分配到P個(gè)線程，具體的分配方式為：若M能被P整除，則

Mi=M/P

(18)

式中，i∈[1,P]且i∈N。若無(wú)特殊說(shuō)明，本節(jié)的i均為此范圍。若M不能被P整除，令余數(shù)為mod(M,P)，則

(19)

顯然，上述分配滿足

(20)

然后，分別跟蹤每個(gè)線程內(nèi)部的試驗(yàn)粒子，并計(jì)算出與航天器表面相交的試驗(yàn)粒子數(shù)Qi。

投影面積的統(tǒng)計(jì)需要對(duì)P個(gè)線程中與航天器表面相交的試驗(yàn)粒子數(shù)Qi進(jìn)行歸約求和，求出所有線程與航天器表面相交的試驗(yàn)粒子數(shù)Q，即

(21)

最終，可通過(guò)公式(16)計(jì)算出投影面積。

需要說(shuō)明的是，上述試驗(yàn)粒子向各個(gè)線程的分配以及對(duì)所有線程的歸約求和都直接通過(guò)OpenMP提供的指令“!$omp parallel do”完成，并將Q聲明為歸約求和 (reduction( + :Q) )變量即可。并行處理只需要在對(duì)應(yīng)串行程序的基礎(chǔ)上添加幾行OpenMP指令，實(shí)現(xiàn)過(guò)程容易，這也是OpenMP的一大優(yōu)勢(shì)。

4　方法校驗(yàn)

4.1　算例描述

考慮外形簡(jiǎn)單、投影面積存在解析解的單面圓形平板，便于對(duì)比分析，以驗(yàn)證方法的正確性。單面圓形平板示意圖、機(jī)體坐標(biāo)系和投影方向如圖 2所示，軸向、展向和法向分別為X、Y和Z軸，投影軸為z軸。圓形平板半徑為1 m，由于軸對(duì)稱特性，只需要考慮赤經(jīng)α=0°的情形。因此，z軸位于XOZ平面，與X軸的夾角為赤緯δ，與Z軸的夾角為赤緯的余角。采用608個(gè)三角形單元進(jìn)行離散，如圖 2所示，試驗(yàn)粒子數(shù)為2×107，此時(shí)統(tǒng)計(jì)誤差為O(10-4)。

圖2　單面圓形平板其三角形面網(wǎng)格離散Fig.2　The one-side circular plate and computational grid

4.2　計(jì)算結(jié)果

圖3給出了該單面平板投影面積隨赤緯的變化，赤緯的變化范圍是δ∈[0°,90°]，即投影方向從平行于平板表面變化到垂直于平板表面。投影方向平行于平板表面時(shí)(δ=0°)，投影面積為0；投影方向垂直于平板表面時(shí)(δ=90°)，投影面積為平板面積π；隨著投影方向從平行于平板表面變化到垂直于平板表面時(shí)，投影面積呈正弦函數(shù)增大。

對(duì)于該簡(jiǎn)單外形，容易得到投影面積的理論解為平板面積與赤緯正弦的乘積，即

AP=πcos(π/2-δ)=πsinδ

(22)

顯然，本文的Monte Carlo統(tǒng)計(jì)結(jié)果與理論解完全一致，驗(yàn)證了方法的可靠性。

圖3　單面圓形平板投影面積隨赤緯的變化Fig.3　The projected area of the circular plate vs. declination

5　SAMSON衛(wèi)星投影面積

SAMSON 衛(wèi)星，全稱為Space Autonomous Mission for Swarming and Geo-locating Nanosatellites 衛(wèi)星，屬于歐空局的小衛(wèi)星[31]。SAMSON衛(wèi)星裝載了2個(gè)太陽(yáng)帆，太陽(yáng)帆的展開(kāi)和折疊、投影方向的改變導(dǎo)致其投影面積跨度較大。SAMSON衛(wèi)星復(fù)雜的工程外形和投影面積跨度較大的特性，使得其既可以用來(lái)考察方法對(duì)復(fù)雜工程外形的適用性，也有助于分析投影面積隨投影方向的變化規(guī)律。

SAMSON的真實(shí)模型如圖 4a所示，計(jì)算采用的三維模型為真實(shí)模型的主要輪廓，略去了天線等對(duì)投影面積影響微小的部件，如圖 4b所示，圖中的坐標(biāo)系為機(jī)體坐標(biāo)系。

圖4　SAMSON衛(wèi)星Fig.4　The SAMSON satellite

方向角/(°)投影面積/m2相對(duì)誤差δα本文結(jié)果文獻(xiàn)[26]結(jié)果ER/(%)000．0300040．03000．013330300．0559170．05590．030410600．0669080．06690．011960900．0599940．0600-0．010003000．1075920．1076-0．0074330300．1337890．1338-0．0082230600．1348140．13480．0103930900．1181900．1182-0．008466000．1680380．16800．0226260300．1784870．1785-0．0072860600．1779830．1780-0．0095560900．1680790．16800．0470290-0．1940010．19400．00052

5.1　對(duì)比校驗(yàn)

表1列出了SAMSON衛(wèi)星沿幾個(gè)典型方向(赤緯δ和赤經(jīng)α)的投影面積，分別給出本文Monte Carlo方法和文獻(xiàn)[26]報(bào)道的結(jié)果，ER為兩者的相對(duì)誤差。赤緯δ=90°時(shí)，投影方向與機(jī)體坐標(biāo)系Z軸重合，故投影面積與赤經(jīng)無(wú)關(guān)。SAMSON衛(wèi)星主要由平面構(gòu)成，采用56個(gè)三角形單元對(duì)其表面進(jìn)行離散即可獲得準(zhǔn)確的空間分辨率，試驗(yàn)粒子總數(shù)為2×107，對(duì)應(yīng)的統(tǒng)計(jì)誤差為O(10-4)，滿足工程需求。從表中可以看到，本文計(jì)算結(jié)果和文獻(xiàn)報(bào)道的數(shù)據(jù)一致，且相對(duì)誤差的量級(jí)為O(10-4)，與統(tǒng)計(jì)誤差量級(jí)相同。因此，本文發(fā)展的Monte Carlo技術(shù)適用于復(fù)雜工程應(yīng)用外形，具備準(zhǔn)確計(jì)算復(fù)雜外形航天器沿任意方向投影面積的能力，且相對(duì)誤差滿足統(tǒng)計(jì)分布規(guī)律。

5.2　投影面積變化規(guī)律

圖 5是SAMSON 衛(wèi)星沿不同方向的投影面積，(a)為投影面積三維分布圖，(b)為對(duì)應(yīng)的二維等值線圖，赤緯的變化范圍是δ∈[-90°,90°]，赤經(jīng)的變化范圍為α∈[0°,360°]。赤緯δ=-90°時(shí)，投影方向?yàn)閆軸反方向；赤緯δ=0°時(shí)，投影方向平行于XOY平面；赤緯δ=90°時(shí)，投影方向?yàn)閆軸正方向。顯然，SAMSON 衛(wèi)星沿不同方向的投影面積變化范圍跨度較大，從最小值A(chǔ)Pmin=0.03 m2變化到最大值A(chǔ)Pmax=0.194 m2。此外，當(dāng)太陽(yáng)帆展開(kāi)時(shí)，赤緯，即投影方向與XOY平面的夾角，對(duì)投影面積影響更為明顯。相對(duì)而言，赤經(jīng)的影響更小一些，這是SAMSON衛(wèi)星的幾何外形所決定的。在赤緯δ=0°附近，投影面積最小，此時(shí)赤經(jīng)對(duì)投影面積的影響最大；赤緯δ=±90°時(shí)，投影面積最大,此時(shí)赤經(jīng)對(duì)投影面積沒(méi)有影響。圖 5揭示了衛(wèi)星姿態(tài)與其投影面積的關(guān)系,是衛(wèi)星姿態(tài)控制的必要輸入條件。

圖5　SAMSON衛(wèi)星沿不同方向的投影面積Fig.5　Projected areas of the SAMSON along different projected directions

6　GOCE衛(wèi)星投影面積

作為投影面積Monte Carlo模擬方法在航天工程的應(yīng)用，采用該方法計(jì)算GOCE衛(wèi)星沿不同方向的投影面積。

6.1　算例描述

GOCE衛(wèi)星，全稱為重力場(chǎng)和穩(wěn)態(tài)海洋環(huán)流探測(cè)者(Gravity Field and Steady State Ocean Circulation Explorer)[3]。GOCE衛(wèi)星主要由衛(wèi)星本體、太陽(yáng)翼和穩(wěn)定翼構(gòu)成，本體呈八棱柱體，橫截面是八邊形，太陽(yáng)翼面積較大，穩(wěn)定翼面積較小，如圖 6所示，圖中的坐標(biāo)系為機(jī)體坐標(biāo)系。GOCE衛(wèi)星雖然外形比較復(fù)雜，但主要由平面構(gòu)成，采用190個(gè)三角形單元對(duì)其表面進(jìn)行離散即可獲得準(zhǔn)確的空間分辨率，試驗(yàn)粒子總數(shù)為2×107，對(duì)應(yīng)的統(tǒng)計(jì)誤差為O(10-4)，滿足工程需求。

圖6　GOCE衛(wèi)星計(jì)算模型Fig.6　The computational model of GOCE satellite

6.2　計(jì)算結(jié)果

圖7是GOCE衛(wèi)星沿不同方向的投影面積，(a)為投影面積三維分布圖，(b)為對(duì)應(yīng)的二維等值線圖，赤緯的變化范圍是δ∈[-90°,90°]，赤經(jīng)的變化范圍為α∈[0°,360°]。赤緯δ=-90°時(shí)，投影方向?yàn)閆軸反方向；赤緯δ=0°時(shí)，投影方向平行于XOY平面；赤緯δ=90°時(shí)，投影方向?yàn)閆軸正方向。與SAMON衛(wèi)星類似，在赤緯δ=0°附近，GOCE投影面積最小，此時(shí)赤經(jīng)對(duì)投影面積的影響最大；赤緯δ=±90°時(shí)，投影面積最大，此時(shí)赤經(jīng)對(duì)投影面積沒(méi)有影響。然而，GOCE 衛(wèi)星細(xì)長(zhǎng)體的外形，外加太陽(yáng)翼較大的面積，導(dǎo)致其沿不同方向的投影面積跨度比SAMSON衛(wèi)星更大，從最小值A(chǔ)Pmin=0.817 m2變化到最大值A(chǔ)Pmax=10.273 m2。隨著其姿態(tài)(即赤緯和赤經(jīng))的改變，GOCE衛(wèi)星投影面積如此之大的變化范圍導(dǎo)致其所受大氣阻力和太陽(yáng)光壓變化巨大，同時(shí)也會(huì)引起太陽(yáng)翼發(fā)電效率的劇烈變化，對(duì)其穩(wěn)定性、姿態(tài)和軌道控制提出了巨大的挑戰(zhàn)。圖 7揭示了衛(wèi)星姿態(tài)與其投影面積的關(guān)系，是太陽(yáng)翼發(fā)電量計(jì)算與監(jiān)測(cè)、大氣阻力和太陽(yáng)光壓預(yù)測(cè)的輸入?yún)?shù)，也是衛(wèi)星姿態(tài)和軌道控制的必要條件。

圖7　GOCE衛(wèi)星沿不同方向的投影面積Fig.7　Projected areas of the GOCE along different projected directions

7　并行性能

并行性能最重要的2個(gè)指標(biāo)是加速比和并行效率?？疾祛惡教祜w機(jī)外形，采用8746個(gè)三角形單元進(jìn)行表面離散，試驗(yàn)粒子數(shù)均為1.6×107。表 2列出了采用典型計(jì)算狀態(tài)的加速比和并行效率，圖 8則給出了加速比隨CPU數(shù)的變化曲線?？紤]到并行的額外開(kāi)銷，比如線程的分叉 (fork)和合并 (join)、歸約之前的時(shí)間同步，以及歸約求和不能完全并行化等因素，本文的加速比和并行效率比較理想。而且，隨著試驗(yàn)粒子數(shù)和表面三角形單元的增加，相對(duì)于并行額外開(kāi)銷，有效并行計(jì)算時(shí)間所占比例增大，同時(shí)，不能完全并行處理所占時(shí)間比例減小，加速比和并行效率會(huì)進(jìn)一步得到改善。

需要注意的是，大多數(shù)衛(wèi)星雖然裝備太陽(yáng)帆導(dǎo)致其外形較復(fù)雜(比如SAMSON和GOCE衛(wèi)星)，但是衛(wèi)星本體和太陽(yáng)帆一般都由平板構(gòu)成，采用數(shù)百甚至數(shù)十個(gè)三角形單元即可對(duì)其表面進(jìn)行準(zhǔn)確離散，單個(gè)CPU計(jì)算也僅僅需要數(shù)十秒的時(shí)間，完全沒(méi)有必要并行計(jì)算。對(duì)于需要采用成千上萬(wàn)個(gè)三角形單元進(jìn)行表面離散的航天飛機(jī)類外形，并行計(jì)算必不可少，此時(shí)的并行加速比和并行效率都較高。

表2　并行性能Table 2　The parallel performance

圖8　多核處理器的OpenMP并行加速比Fig.8　OpenMP parallel speedup for multiple-core processors

8　結(jié)　論

借鑒統(tǒng)計(jì)學(xué)中Monte Carlo技術(shù)的思想，結(jié)合計(jì)算幾何學(xué)中基本的射線與三角形相交理論，發(fā)展了一種魯棒性和通用性較好的計(jì)算復(fù)雜外形航天器投影面積的Monte Carlo方法，實(shí)現(xiàn)了該方法的并行化。將該方法應(yīng)用于單面圓形平板模型、SAMSON小衛(wèi)星和GOCE低軌衛(wèi)星投影面積的計(jì)算。結(jié)論如下：

(1)單面圓形平板投影面積Monte Carlo模擬結(jié)果和理論解的一致以及SAMSON小衛(wèi)星投影面積與文獻(xiàn)報(bào)道結(jié)果的一致性則進(jìn)一步驗(yàn)證了該方法的可靠性和對(duì)復(fù)雜工程外形的適用性。

(2)本文發(fā)展的投影面積的Monte Carlo方法的計(jì)算量大致與試驗(yàn)粒子總數(shù)與表面三角形單元數(shù)的乘積呈正比，采用共享存儲(chǔ)模式的OpenMP很容易實(shí)現(xiàn)并行化，且并行加速比和并行效率都較高。

(3)SAMSON和GOCE衛(wèi)星在不同姿態(tài)下的投影面積變化范圍跨度較大，赤緯對(duì)投影面積的影響明顯，赤經(jīng)的影響小一些，在赤緯δ=0°附近，投影面積最小，此時(shí)赤經(jīng)對(duì)投影面積的影響最大，赤緯δ=±90°時(shí)，投影面積最大，此時(shí)赤經(jīng)對(duì)投影面積沒(méi)有影響。

(4)本文發(fā)展的Monte Carlo方法能夠給出衛(wèi)星姿態(tài)與其投影面積的關(guān)系，是太陽(yáng)翼發(fā)電量計(jì)算與監(jiān)測(cè)、大氣阻力和太陽(yáng)光壓預(yù)測(cè)的輸入?yún)?shù)，也是衛(wèi)星姿態(tài)和軌道控制的必要條件。

[1]Tapley B D, Bettadpur S, Ries J C, et al. GRACE measurements of mass variability in the Earth system [J]. Science 2004, 305 (5683): 503-505.

[2]Reigber C, Luehr H, Schwintzer P. CHAMP mission status [J]. Advances in Space Research, 2002, 30(2): 129-134

[3]Rummel R. Yi W, Stummer C. GOCE gravitational gradiometry [J].Journal of Geodesy, 2011, 85(11): 777-790.

[4]Canuto E. Drag-free and alltitude control for the GOCE satellite [J]. Automatica, 2008, 44: 1766-1780.

[5]Jacchia L G. Thermospheric temperature, density, and composition: New models [R]. Special Report No. 375. Smithsonian Astrophysical Observatory, Cambridge, MA, 1977.

[6]Picone J M, Hedin A E, and Drob D P. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues [J]. Journal of Geophysical Research, 2002, 107(A12): 1468-1483.

[7]Bowmana B R, Tobiskab W K, Marcosc F A, et al. The JB2006 empirical thermospheric density model [J]. Journal of Atmospheric and Solar-Terrestrial Physics, 2008, 70: 774-793.

[8]Bruinsma S. The DTM-2013 thermosphere model [J]. Journal of Space Weather and Space Climate, 2015, 5(A1): 1-8.

[9]Cook G E. Satellite drag coefficients [J]. Planetary and Space Science, 1965, 13: 929-946.

[10]McLaughlin C A, Mance S, and Lechtenberg T. Drag coefficient estimation in orbit determination [J]. The Journal of Astronautical Sciences, 2011, 58(3): 513-530.

[11]Pilinski M D,Argrow B M, and Palo S E. Drag coefficients of satellites with concave geometries: comparing models and observations [J]. Journal of Spacecraft and Rocket, 2011, 48(2): 312-325.

[12]Mehta M M, McLaughlin C A, Sutton E K. Drag coefficient modeling for grace using Direct Simulation Monte Carlo [J]. Advances in Space Research, 2013, 52: 2035-2051.

[13]王智慧, 鮑麟, 童秉綱. 尖化前緣的稀薄氣體化學(xué)非平衡流動(dòng)和氣動(dòng)加熱相似律研究[J]. 氣體物理, 2016, 1(1): 5-12. [Wang Zhi-hui, Bao Lin, Tong Bing-gang. Similarity law of aero-heating to sharpened noses in rarefied chemical nonequilibrium flows [J]. Physics of Gases, 2016, 1(1): 5-12.]

[14]周偉勇, 張育林, 劉昆. 超低軌航天器氣動(dòng)力分析與減阻設(shè)計(jì) [J]. 宇航學(xué)報(bào), 2010, 31(2): 342-348. [Zhou Wei-yong, Zhang Yu-lin, Liu Kun. Aerodynamics analysis and reduced drag design for the lower LEO spacecraft [J]. Journal of Astronautics, 2010, 31(2): 342-348.]

[15]Fuller J D, Tolson R H. Improved method for estimation of spacecraft free-molecular aerodynamic properties [J]. Journal of Spacecraft and Rocket, 2009, 46(5): 938-947.

[16]靳旭紅, 黃飛, 程曉麗, 等. 內(nèi)外流一體化航天器氣動(dòng)特性分析與減阻設(shè)計(jì). 宇航學(xué)報(bào), 2017, 38(1): 10-17. [Jin Xu-hong, Huang Fei, Cheng Xiao-li, et al. Analysis of aerodynamic properties and drag-reduction design for spacecraft with an open orifice [J]. Journal of Astronautics, 2017, 38(1): 10-17.]

[17]Stone W C, Witzgall C. Evaluation of aerodynamic drag and torque for external tanks in low earth orbit.Journal of Research of the National Institute of Standards and Technology, 2006, 111 (2): 143-159.

[18]Ha J, Mimasu Y, Tsuda Y, et al. Solar and Thermal radiation models and flight evaluation for IKAROS solar sail [J]. Journal of Spacecraft and Rocket, 2015, 52(3): 958-967.

[19]白顯宗, 陳磊, 張翼, 等. 空間目標(biāo)碰撞預(yù)警技術(shù)研究綜述. 宇航學(xué)報(bào), 2013, 34(8): 1027-1039. [Bai Xian-zong, Chen Lei, Zhang Yi, et al. Survey on collision assessment and warning techniques for space object [J]. Journal of Astronautics, 2013, 34(8): 1027-1039.]

[20]Vallado D A. Fundamentals of astrodynamics and applications[M]. fourth edition. Hawthorne: Microcosm, 2013.

[21]Ashish T. Adaptive Vectored thrust deorbiting of space debris [J]. Journal of Spacecraft and Rocket, 2013, 50(2): 394-401.

[22]Markley F L, Crassidis J L. Fundamentals of spacecraft attitude determination and control [M]. New York: Springer, 2014.

[23]Appel A. Some techniques for shading machine renderings of solids[J].ACM Siggraph Comput. Graph, 1968, 2(1): 37-45.

[24]Atherton P, Weiler K, Greenberg D. Polygon shadow generation[J]. ACM Siggraph Comput. Graph, 1978, 12(3): 275-281.

[25]Woo A, Poulin P, Fournier A. A survey of shadow algorithms[J]. IEEE Computer Graphics & Applications, 1990, 10(6): 13-32.

[26]Ben Y O, Edlerman E, Curfil P. Analytical technique for satellite projected cross-sectional area calculation[J]. Advances in Space Research, 2015, 56: 205-217.

[27]Curtis H D. Orbital mechanics for engineering students[M]. Third edition, Hawthorne: Butterworth-Heinemann, 2014.

[28]Zio E. The Monte Carlo Simulation Method for System Reliability and Risk Analysis, Springer Series in Reliability Engineering[M]. London: Springer-Verlag, 2013.

[29]Bonet J, Peraire J. An Alternating Digital Tree (ADT) algorithm for geometric searching and intersection problems[J]. International Journal for Numerical Methods in Engineering, 1991, 31: 1-17.

[30]茆詩(shī)松, 程依明, 濮小龍. 概率論與數(shù)理統(tǒng)計(jì)教程 [M]. 北京: 高等教育出版社, 2011.

[31]Ben Y O, Gurfil P. Stability and performance of orbital elements feedback for cluster keeping using differential drag [J]. Journal of the Astronautical Sciences, 2014, 61(1): 1-29.

[32]Gao D, Schwartzentruber T E. Optimizations and OpenMP implementation for the direct simulation Monte Carlo method [J]. Computers and Fluids, 2011, 42(1): 73-81.