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

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

0　引　言

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

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

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

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

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

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

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

1.1　坐標(biāo)系定義

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

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

δ=arcsinzZ

(1)

(2)

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

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

(3)

(4)

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

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

2　投影面積Monte Carlo模擬

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

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

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

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

(5)

2.2　邊界矩形的確定

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

z

(6)

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

z=zmin-1

(7)

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

(8)

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

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

(9)

式中，R1、R2均為(0,1)之間均勻分布的隨機數(shù)。由于投影方向為oz方向，試驗粒子的運動方向為

v=(0,0,1)

(10)

2.4　試驗粒子的追蹤

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

(11)

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

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

(12)

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

(13)

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

(14)

若同時滿足

(15)

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

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

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

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

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

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

(16)

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

(17)

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

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

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

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

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

Mi=M/P

(18)

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

(19)

顯然，上述分配滿足

(20)

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

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

(21)

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

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

4　方法校驗

4.1　算例描述

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

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

4.2　計算結(jié)果

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

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

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

(22)

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

圖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個太陽帆，太陽帆的展開和折疊、投影方向的改變導(dǎo)致其投影面積跨度較大。SAMSON衛(wèi)星復(fù)雜的工程外形和投影面積跨度較大的特性，使得其既可以用來考察方法對復(fù)雜工程外形的適用性，也有助于分析投影面積隨投影方向的變化規(guī)律。

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

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

方向角/(°)投影面積/m2相對誤差δα本文結(jié)果文獻[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　對比校驗

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

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

圖 5是SAMSON 衛(wèi)星沿不同方向的投影面積，(a)為投影面積三維分布圖，(b)為對應(yīng)的二維等值線圖，赤緯的變化范圍是δ∈[-90°,90°]，赤經(jīng)的變化范圍為α∈[0°,360°]。赤緯δ=-90°時，投影方向為Z軸反方向；赤緯δ=0°時，投影方向平行于XOY平面；赤緯δ=90°時，投影方向為Z軸正方向。顯然，SAMSON 衛(wèi)星沿不同方向的投影面積變化范圍跨度較大，從最小值A(chǔ)Pmin=0.03 m2變化到最大值A(chǔ)Pmax=0.194 m2。此外，當(dāng)太陽帆展開時，赤緯，即投影方向與XOY平面的夾角，對投影面積影響更為明顯。相對而言，赤經(jīng)的影響更小一些，這是SAMSON衛(wèi)星的幾何外形所決定的。在赤緯δ=0°附近，投影面積最小，此時赤經(jīng)對投影面積的影響最大；赤緯δ=±90°時，投影面積最大,此時赤經(jīng)對投影面積沒有影響。圖 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)用，采用該方法計算GOCE衛(wèi)星沿不同方向的投影面積。

6.1　算例描述

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

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

6.2　計算結(jié)果

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

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

7　并行性能

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

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

表2　并行性能Table 2　The parallel performance

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

8　結(jié)　論

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

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

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

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

(4)本文發(fā)展的Monte Carlo方法能夠給出衛(wèi)星姿態(tài)與其投影面積的關(guān)系，是太陽翼發(fā)電量計算與監(jiān)測、大氣阻力和太陽光壓預(yù)測的輸入?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é)非平衡流動和氣動加熱相似律研究[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]周偉勇, 張育林, 劉昆. 超低軌航天器氣動力分析與減阻設(shè)計 [J]. 宇航學(xué)報, 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)外流一體化航天器氣動特性分析與減阻設(shè)計. 宇航學(xué)報, 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é)報, 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ù)理統(tǒng)計教程 [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.