類天宮航天器迎風面積建模及變化特性分析

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

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

1 引言

大氣阻力是航天器在軌運行的一種重要非保守力，對軌道處理特別是中長期軌道預報精度具有重大影響［1-3］。航天器迎風面積是計算大氣阻力的重要因素，定義為沿航天器相對于當地大氣速度方向航天器幾何結構的投影面積。

大型航天器一般由多功能艙段組成，且具有較大的太陽帆板，構成一個復雜的幾何外形。航天器在軌運行中姿態的變化導致迎風面也在不斷變化，同時太陽帆板會根據太陽方位不斷調整角度，對航天器的外形產生影響。因此，在航天器運行過程中，其迎風面積將發生非常大的變化。

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

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

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

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

2 迎風面積幾何投影算法

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

2.1 模型表示

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

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

圖2 航天器部件幾何外形的非結構網格Fig.2 Unstructured triangular surface elementmesh of spacecraft com ponent

2.2 幾何投影

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

3 迎風面積建模方法

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

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

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

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

3.1 離線建模

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

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

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

已知投影方向對應的經度、緯度和帆板轉角即可確定迎風面積計算狀態，本文的迎風面積模型使用投影經度、投影緯度、帆板轉角作為模型參數，記為［Pj，Pw，Rf］，其中Pj∈ ［0，360），Pw∈［-90，90］，Rf∈ ［0，180）。

按照模型參數Pj、Pw、Rf的取值范圍，對每個參數以一定步長取值，產生所有參數對，根據公式（1）計算投影方向，然后利用投影方向和帆板轉角，應用幾何投影法計算該參數對應狀態的迎風面積，形成迎風面積模型。

3.2 實時計算

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

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

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

3.2.1 正飛模式

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

3.2.2 偏航模式

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

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

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

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

4 模型方法驗證

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

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

4.1 模型建立

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

4.2 在軌運行期間迎風面積變化情況

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

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

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

圖9 迎風面積方法計算的帆板轉動角度和遙測實際數據比較Fig.9 Comparison between telemetry data and calculation of the sailboard rotation angle

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

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

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

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

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

5 結論

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

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