類天宮航天器迎風(fēng)面積建模及變化特性分析

楊 成，李 勰，孫 軍，李志輝

（1.航天飛行動力學(xué)技術(shù)重點實驗室，北京100094；2.北京航天飛行控制中心，北京100094；3.中國空氣動力研究與發(fā)展中心超高速空氣動力研究所，綿陽621000；4.空氣動力前沿技術(shù)研究中心國家計算流力學(xué)實驗室，北京100191）

1 引言

大氣阻力是航天器在軌運行的一種重要非保守力，對軌道處理特別是中長期軌道預(yù)報精度具有重大影響［1-3］。航天器迎風(fēng)面積是計算大氣阻力的重要因素，定義為沿航天器相對于當?shù)卮髿馑俣确较蚝教炱鲙缀谓Y(jié)構(gòu)的投影面積。

大型航天器一般由多功能艙段組成，且具有較大的太陽帆板，構(gòu)成一個復(fù)雜的幾何外形。航天器在軌運行中姿態(tài)的變化導(dǎo)致迎風(fēng)面也在不斷變化，同時太陽帆板會根據(jù)太陽方位不斷調(diào)整角度，對航天器的外形產(chǎn)生影響。因此，在航天器運行過程中，其迎風(fēng)面積將發(fā)生非常大的變化。

工程計算中通常將迎風(fēng)面積作為待估計參數(shù)，使用大氣阻力系數(shù)與面質(zhì)比乘積的等效彈道系數(shù)法［3-4］，結(jié)合歷史軌道數(shù)據(jù)，選取一定時間范圍（一周或更長時間），根據(jù)該時段大型航天器軌道半長軸衰減情況，估算航天器等效彈道系數(shù)，進行定軌與軌道外推壽命預(yù)報。根據(jù)前期數(shù)據(jù)將迎風(fēng)面積作為參數(shù)進行擬合，在航天器飛行預(yù)報階段會引入較大誤差，從而影響軌道預(yù)報精度。因此，根據(jù)航天器的實際幾何外形、運行姿態(tài)等信息，實時地計算航天器迎風(fēng)面積，對提高軌道處理精度有重要意義。

航天器一般包括本體和太陽帆板兩部分，為便于處理，本體簡化為圓柱形，太陽帆板簡化為長方體。唐歌實等［3］提出了一種航天器俯仰滾動偏航角均為0°的正飛姿態(tài)情況下的等效迎風(fēng)面積計算方法。首先計算太陽帆板在繞旋轉(zhuǎn)軸一周的面積變化，求取平均值，加上本體的截面積即為迎風(fēng)面積，實際所得是航天器正飛狀態(tài)的平均迎風(fēng)面積。這種方法是不同狀態(tài)下的數(shù)值平均，并不能代替實際情況下的迎風(fēng)面積變化；另一類為近似外形幾何投影方法［5-6］。首先將航天器外形簡化，分解為方形、柱形等基礎(chǔ)形狀的組合。分別對這些簡單形狀在航天器速度方向進行投影處理，獲得各自在投影面上的投影區(qū)域，得到航天器的迎風(fēng)面積。這種方法采用簡單幾何體來對航天器外形進行簡化，在計算精度上難以滿足精密定軌的要求，并且計算過程較為復(fù)雜，缺乏通用性。

楊成等［7］提出了一種應(yīng)用陰影圖計算迎風(fēng)面積的方法，通過對航天器幾何外形進行精確網(wǎng)格建模，同時考慮航天器姿態(tài)變化和帆板轉(zhuǎn)動情況，利用幾何投影方式對航天器進行投影；根據(jù)投影面積，確定航天器的迎風(fēng)面積，并以天宮一號航天器為例，對其理論上的迎風(fēng)面積進行了計算。朱戰(zhàn)霞等［8］基于微元劃分思想提出了射擊線掃描法以求解有效迎風(fēng)面積。這類方法在外形描述和投影計算方面具有較好的通用性，但計算過程需要的航天器狀態(tài)參數(shù)較多，在實際工程應(yīng)用中受測控條件影響［9］，在測控弧段外航天器狀態(tài)數(shù)據(jù)缺失，此時如何建立合適的計算模型并解決模型參數(shù)獲取是個問題。

本文提出了一種類天宮航天器迎風(fēng)面積建模方法。該模型按照地球經(jīng)緯度定義方式，在航天器本體系建立經(jīng)緯度網(wǎng)格，將計算迎風(fēng)面積的投影方向轉(zhuǎn)換為經(jīng)度和緯度，與太陽帆板的轉(zhuǎn)動角度一起構(gòu)成迎風(fēng)面積模型的參數(shù)；在模型參數(shù)取值范圍內(nèi)按照一定的步長，利用幾何投影［7］計算所有狀態(tài)的迎風(fēng)面積數(shù)據(jù)，形成航天器的迎風(fēng)面積模型；在軌道處理過程中，按照航天器長期飛行期間的控制模式，計算航天器各時刻的姿態(tài)和帆板轉(zhuǎn)動角度，確定對應(yīng)的模型參數(shù)，通過插值獲得迎風(fēng)面積。

2 迎風(fēng)面積幾何投影算法

由于航天器在軌飛行過程中姿態(tài)變化，帆板等大尺寸部件發(fā)生相對運動，迎風(fēng)面積將發(fā)生較大的變化。精確迎風(fēng)面積計算方法［7］通過構(gòu)建航天器精確的幾何外形，根據(jù)實際航天器姿態(tài)、活動部件運動變形情況，在相對速度方向進行幾何投影，利用產(chǎn)生的投影面積即可確定對應(yīng)迎風(fēng)面積的大小。基本原理包括模型表示和幾何投影兩部分，如圖1所示。利用非結(jié)構(gòu)化三角網(wǎng)格表征航天器的幾何外形，使用樹結(jié)構(gòu)表征航天器部件之間的關(guān)系，根據(jù)航天器姿態(tài)、帆板轉(zhuǎn)動角度，在航天器相對大氣速度方向進行幾何投影，獲得航天器該狀態(tài)下的迎風(fēng)面積。

2.1 模型表示

使用樹結(jié)構(gòu)表示航天器活動部件之間的層次連接關(guān)系和運動狀態(tài)。樹結(jié)構(gòu)的節(jié)點對應(yīng)一個航天器部件，使用非結(jié)構(gòu)三角網(wǎng)格表征該部件實體的幾何外形信息，可逼近任何復(fù)雜外觀和結(jié)構(gòu)的航天器，提高模型表示精度［10-11］，同時降低網(wǎng)格數(shù)量。利用非結(jié)構(gòu)網(wǎng)格表示航天器如圖2所示，幾何形狀信息里包括頂點和面片兩部分信息，其中頂點信息包括所有頂點的三維坐標，面片信息包括組成幾何形狀外表面所有面片的頂點序號。按照航天器的結(jié)構(gòu)，對不同部分使用單獨的三角網(wǎng)格進行表示，平緩表面網(wǎng)格稀疏，彎曲表面使用更稠密網(wǎng)格擬合。

圖1 類天宮航天器迎風(fēng)面積計算方法Fig.1 Calculation method of cross section area of TG-like spacecraft

圖2 航天器部件幾何外形的非結(jié)構(gòu)網(wǎng)格Fig.2 Unstructured triangular surface elementmesh of spacecraft com ponent

2.2 幾何投影

航天器迎風(fēng)面積幾何投影原理如圖3所示。以航天器相對大氣的速度矢量方向為投影方向，利用計算機圖形學(xué)方法［11］，對航天器幾何網(wǎng)格進行投影，在垂直于投影方向的投影面產(chǎn)生航天器投影的影像。通過計算影像所占面積大小，即可確定航天器的投影面積。

3 迎風(fēng)面積建模方法

軌道計算一般包括定軌和預(yù)報兩部分，需要對一定弧段內(nèi)若干時刻的迎風(fēng)面積進行計算。在軌道處理過程中，直接采用2.1節(jié)的方法，若頻繁調(diào)用迎風(fēng)面積計算模塊，以圖形學(xué)幾何投影原理來獲取某個時刻的迎風(fēng)面積數(shù)值，將影響軌道處理速度。

圖3 航天器幾何投影示意圖Fig.3 Geometric projection diagram of spacecraft

根據(jù)類天宮航天器外形和在軌飛行特點，本文使用投影經(jīng)度、投影緯度、帆板角度3個參數(shù)對計算狀態(tài)進行表示，根據(jù)各參數(shù)取值范圍，以一定步長計算狀態(tài)空間中所有可能狀態(tài)的迎風(fēng)面積，離線建立迎風(fēng)面積模型。在軌道處理過程中，根據(jù)飛行狀態(tài)和飛行模式，計算航天器相對大氣的速度、飛行姿態(tài)［12］、帆板轉(zhuǎn)角，并轉(zhuǎn)化為投影經(jīng)度、投影緯度、帆板角度3個模型參數(shù)，使用線性插值的方式從模型中獲得迎風(fēng)面積，如圖4所示。

圖4 類天宮航天器迎風(fēng)面積建模方法Fig.4 M odeling on cross section area for Tiangongtype spacecraft

3.1 離線建模

類天宮航天器的主要結(jié)構(gòu)包括本體和太陽帆板，本體近似圓柱體，太陽帆板連接在本體兩側(cè)，其迎風(fēng)面積可認為與投影方向和帆板轉(zhuǎn)角有關(guān)［12］。引入經(jīng)緯度球的概念，將投影方向矢量從三維轉(zhuǎn)化為經(jīng)度、緯度二維，降低模型參數(shù)個數(shù)。參照地球地理信息系統(tǒng)中的經(jīng)緯度定義［3］，在航天器本體坐標系建立經(jīng)緯度球，其中：x軸和z軸構(gòu)成的平面指向x軸正方向?qū)?yīng)經(jīng)度為0°，y軸和z軸構(gòu)成的平面指向y軸正方向?qū)?yīng)經(jīng)度為90°，經(jīng)度范圍為［0°，360°）；x軸和y軸構(gòu)成的平面對應(yīng)緯度為0°，沿z軸正方向緯度增加，北極為90°，沿z軸負方向緯度降低，南極為-90°，緯度范圍為［-90°，90°］，如圖5所示。

圖5 航天器本體系下的相對速度Fig.5 Relative velocity of spacecraft in body fixed coordinate system

航天器相對大氣的速度歸一化為單位向量［Vx，Vy，Vz］，即為計算迎風(fēng)面積的投影方向，對應(yīng)航天器經(jīng)緯度球面的一個點，可以用該點的經(jīng)度Pj和緯度Pw進行表示。［Vx，Vy，Vz］和［Pj，Pw］之間的轉(zhuǎn)化關(guān)系如下：

已知投影方向?qū)?yīng)的經(jīng)度、緯度和帆板轉(zhuǎn)角即可確定迎風(fēng)面積計算狀態(tài)，本文的迎風(fēng)面積模型使用投影經(jīng)度、投影緯度、帆板轉(zhuǎn)角作為模型參數(shù)，記為［Pj，Pw，Rf］，其中Pj∈ ［0，360），Pw∈［-90，90］，Rf∈ ［0，180）。

按照模型參數(shù)Pj、Pw、Rf的取值范圍，對每個參數(shù)以一定步長取值，產(chǎn)生所有參數(shù)對，根據(jù)公式（1）計算投影方向，然后利用投影方向和帆板轉(zhuǎn)角，應(yīng)用幾何投影法計算該參數(shù)對應(yīng)狀態(tài)的迎風(fēng)面積，形成迎風(fēng)面積模型。

3.2 實時計算

在對軌道數(shù)據(jù)進行定軌和利用軌道根數(shù)進行軌道預(yù)報過程中，都需要軌道處理弧段內(nèi)的航天器迎風(fēng)面積數(shù)據(jù)。航天器在軌長期飛行期間，為保證太陽帆板獲取太陽能的效率，姿態(tài)控制模式在正飛和偏航2種模式之間切換，如圖6所示。定義地心到太陽的方向為太陽矢量方向，當太陽矢量與航天器軌道面之間夾角在閾值以內(nèi)時，采取正飛模式，偏航、俯仰、滾動3個姿態(tài)角為0°，帆板通過調(diào)整角度，使得太陽與航天器連線位于帆板轉(zhuǎn)軸與帆板法線組成的平面內(nèi)，以便在保持正飛姿態(tài)前提下太陽帆板獲得最大太陽能量；當太陽矢量與航天器軌道面之間夾角大于閾值時，采取偏航模式，俯仰、滾動保持0°，同時調(diào)整偏航角度和帆板轉(zhuǎn)動角度，使得太陽與航天器連線垂直于帆板平面，太陽帆板獲得最大太陽能量。

圖6 航天器飛行模式與太陽方位關(guān)系示意圖Fig.6 Diagram of relationship between spacecraftflightmode and solar orientation

在地球慣性系下，記航天器位置為P=［x，y，z］T，速度為下面給出按照飛行模式確定迎風(fēng)面積模型參數(shù)的方法。

3.2.1 正飛模式

正飛模式下，航天器x軸指向速度V方向，z軸指向地心方向，與x軸垂直，y軸與z、x軸形成右手坐標系。航天器相對大氣速度轉(zhuǎn)換到航天器本體坐標系下，按照公式（2）可計算出投影經(jīng)緯度［Pj，Pw］；將太陽方向矢量投影到航天器z、x軸構(gòu)成的平面上，并歸一化為單位向量，在航天器本體系下的坐標記為S=［Sx，0，Sz］T，則太陽帆板轉(zhuǎn)動角度Rf計算如式（3）所示：

3.2.2 偏航模式

偏航模式下，航天器將根據(jù)太陽方向進行偏航，使得航天器的x軸沿z軸轉(zhuǎn)動一個偏航角到達x'位置，與太陽方向矢量S和z軸共面，從而帆板轉(zhuǎn)動一定角度就能使得帆板法線與太陽矢量S重合，獲得最大太陽能轉(zhuǎn)化效率，如圖7所示。x'為太陽方向矢量S在x-y平面的投影，通過S與X矢量的內(nèi)積LSX=S·X和S與Y矢量的內(nèi)積LSY=S·Y，可計算出偏航角度ψ為式（4）所示：

圖7 偏航模式下偏航角與太陽方位示意圖Fig.7 Diagram of yaw angle and sun orientation in yaw mode

航天器偏航ψ角度后，其本體系坐標軸轉(zhuǎn)動到x'、y'和z位置，y'為帆板轉(zhuǎn)動軸方向，由于太陽方向矢量S在x'和z軸構(gòu)成的平面內(nèi)，y'與S是垂直關(guān)系，通過轉(zhuǎn)動帆板即可使得帆板法向量與太陽方向矢量S重合。同理，根據(jù)公式（3）可計算太陽帆板轉(zhuǎn)動角度。

利用航天器飛行模式確定該時刻迎風(fēng)面積模型的參數(shù)［Pj，Pw，Rf］，通過線性插值，得到對應(yīng)迎風(fēng)面積。

4 模型方法驗證

天宮一號為中國目前最大的復(fù)雜航天器，主要由本體、太陽帆板組成，其中太陽帆板可繞本體系的y軸旋轉(zhuǎn)。本體由資源艙和實驗艙組成，如圖8所示，總長約為10 m，直徑約為3 m，太陽帆板展開后長約為20 m。由于該航天器尺寸較大，大氣阻力對軌道的影響較為顯著。本文以天宮一號為研究對象，驗證類天宮航天器迎風(fēng)面積模型，分析其在軌飛行期間迎風(fēng)面積的變化特性。

圖8 類天宮航天器幾何外形Fig.8 Geometric surface of Tiangong-type spacecraft

4.1 模型建立

天宮一號本體部分的幾何外形使用8342個三角面片表示，帆板部分幾何構(gòu)型比較規(guī)則，包含648個三角面片。迎風(fēng)面積模型使用投影經(jīng)度、投影緯度、帆板轉(zhuǎn)角作為模型參數(shù)，記為［Pj，Pw，Rf］，取值范圍Pj∈［0，360），Pw∈［-90，90］，Rf∈ ［0，180）。參數(shù)步長設(shè)置為1°，這樣Pj存在360種取值可能性，Pw存在181種取值可能性，Rf存在180種取值可能性，組合起來對應(yīng)360×181×180種狀態(tài)。幾何投影功能基于OpenGL圖形庫實現(xiàn)［10］，處理電腦配置為CPU i3、內(nèi)存4G、顯存1G，計算所有狀態(tài)的迎風(fēng)面積耗時2 h。

4.2 在軌運行期間迎風(fēng)面積變化情況

利用2015年天宮一號在軌飛行數(shù)據(jù)，驗證本文迎風(fēng)面積計算方法。數(shù)據(jù)為153天的連續(xù)觀測弧段，期間航天器的飛行姿態(tài)模式分為正飛和偏航兩種。整個實驗過程中，隨著太陽矢量方向與軌道面夾角的變化，天宮一號進行了多次正飛和偏航模式的切換。

首先根據(jù)航天器飛行模式，對其偏航角度和帆板轉(zhuǎn)動角度的計算方法進行實驗。受航天器機械控制系統(tǒng)的影響，實際角度與本文計算角度會存在一定誤差。2015年11月1日天宮一號處于偏航模式，在16時左右存在兩段約8 min的測控弧段，偏航角度和帆板轉(zhuǎn)動角度數(shù)據(jù)通過遙測進行下傳，使用該數(shù)據(jù)與本文方法計算的偏航角度和帆板轉(zhuǎn)動角度數(shù)據(jù)進行對比驗證。

圖9為帆板轉(zhuǎn)動角度對比情況。天宮一號帆板安裝在資源艙左右兩側(cè)，將左右帆板轉(zhuǎn)動當作相同的，實際由于受2套獨立機構(gòu)驅(qū)動，左右帆板的實際轉(zhuǎn)動角度會有不同。圖9中虛線為本文計算的帆板轉(zhuǎn)動角度，實線為左右帆板通過遙測下傳的實際轉(zhuǎn)動角度。經(jīng)統(tǒng)計，帆板轉(zhuǎn)動角度的平均誤差為2.65°。圖10為偏航計算角度和遙測下傳實際角度的對比，偏航角度平均誤差為0.53°。從圖9、圖10可以看出，按照本文方法計算角度和遙測下傳實際角度在角度變化上趨勢一致，角度值誤差較小。

圖9 迎風(fēng)面積方法計算的帆板轉(zhuǎn)動角度和遙測實際數(shù)據(jù)比較Fig.9 Comparison between telemetry data and calculation of the sailboard rotation angle

圖10 迎風(fēng)面積模型偏航角度和遙測實際數(shù)據(jù)比較Fig.10 Com parison between telemetry data and calcu lation of the yaw angle

每種模式飛行期間內(nèi)，可以根據(jù)該模式計算該時刻航天器的姿態(tài)和帆板轉(zhuǎn)動角度，得到模型參數(shù)，然后采用插值方式計算迎風(fēng)面積。在153天實驗弧段內(nèi)，每隔30 s進行一次時間采樣，共計算440 640個時刻的迎風(fēng)面積，結(jié)果如圖11所示。圖中可以看出，正飛過程迎風(fēng)面積變化范圍更大，可以達到理論上的極值，該模式下近似為本體前端面和太陽帆板產(chǎn)生迎風(fēng)面積，本體前端面為圓柱形的截面積保持不變，太陽帆板循環(huán)轉(zhuǎn)動360°，當帆板轉(zhuǎn)動到水平和垂直狀態(tài)，迎風(fēng)面積就會接近最小值和最大值；偏航模式中，由于本體存在偏航轉(zhuǎn)動，本體的前端面和側(cè)面都會產(chǎn)生迎風(fēng)面積，所以迎風(fēng)面積的最小值比正飛模式大，由于偏航時本體和太陽帆板之間的遮擋抵消，迎風(fēng)面積的最大值比正飛模式小。

圖11 在軌飛行過程中迎風(fēng)面積變化情況Fig.11 Variation of cross section area during orbit flight

圖12為每隔3 h統(tǒng)計一次的迎風(fēng)面積平均值，由于天宮一號的軌道周期約1.5 h，該平均值反映的是軌道周期時間尺度上迎風(fēng)面積的變化情況。可以看出，正飛模式下，航天器姿態(tài)固定，在每個軌道周期內(nèi)太陽帆板旋轉(zhuǎn)360°，以軌道周期時間范圍進行平均后，迎風(fēng)面積平均值約為36.5～37.5 m2，變化幅度很小；偏航模式下，軌道面與太陽位置關(guān)系確定了航天器的偏航角度，迎風(fēng)面積平均值約為32.6～40.3 m2，在一個偏航模式范圍內(nèi)先變小再變大，變化范圍比正飛模式大。分析可以看出，正飛模式下天宮一號迎風(fēng)面積受帆板轉(zhuǎn)動影響變化幅度很大，但軌道周期尺度的均值幾乎恒定；偏航模式下迎風(fēng)面積的下限更高，軌道周期尺度的均值具有較大的變化幅度。因此，天宮一號迎風(fēng)面積在不同飛行模式下具有不同的變化情況，軌道計算中使用固定等效面積或者彈道系數(shù)的方法不能反映迎風(fēng)面積的這種變化，帶來了較大誤差，需要根據(jù)實際狀態(tài)實時計算。

圖12 每3 h的平均迎風(fēng)面積變化情況Fig.12 Variation of average cross section area every 3 hours

5 結(jié)論

1）本文提出了一種航天器迎風(fēng)面積建模方法，解決了可變復(fù)雜外形航天器迎風(fēng)面積的精確計算問題。該方法能處理航天器復(fù)雜幾何外形、部件運動變形、飛行姿態(tài)變化對迎風(fēng)面積的影響，實驗結(jié)果表明了方法的有效性；

2）對天宮一號在軌飛行期間不同時間尺度的迎風(fēng)面積變化情況進行了分析，正飛模式和偏航模式下迎風(fēng)面積具有不同的極值和變化規(guī)律，表明了姿態(tài)變化、帆板轉(zhuǎn)動因素對類天宮航天器迎風(fēng)面積具有不可忽略的影響。