韓 意，陳 明，張 宇，孔 靜，段成林

1 引言

在空間實驗室任務(wù)階段，各航天器有效迎風(fēng)面積參數(shù)是否準確，直接影響航天器大氣阻力建模、精密定軌以及軌道預(yù)報精度[1]。2016年9月15日天宮二號及其搭載的伴星發(fā)射入軌，10月19日神舟十一號在393 km軌道與天宮二號交會對接構(gòu)成組合體，以三軸對地定向的模式飛行(下文簡稱正飛)。10月23日伴星從天宮二號上成功釋放[2]。在10月23日至11月15日共計約24天的時間里，伴星與組合體共同運行在平均高度約388 km的軌道上。與歷次載人航天工程任務(wù)相比，這是首次有2個航天器長時間運行在相同軌道高度，對各目標精密定軌中大氣阻力系數(shù)求解提供了一致性對比的條件。天宮二號組合體等大型航天器外形復(fù)雜，且姿軌運動形式多樣，不同部件之間存在遮擋關(guān)系，難以準確計算其迎風(fēng)面積。我國即將建造的空間站也飛行在相近軌道高度上，由于尺寸結(jié)構(gòu)以及構(gòu)型更加復(fù)雜，其迎風(fēng)面積計算難度更大。研究該高度大氣阻力特點及航天器迎風(fēng)面積計算方法，對于提高空間站定軌預(yù)報精度具有重要意義。

國內(nèi)外研究人員在航天器迎風(fēng)面積與阻力系數(shù)CD分析方面開展了相關(guān)研究[3]，盧偉[4]提出了迎風(fēng)面積包絡(luò)微元計算法，解決了具有規(guī)則表面形狀衛(wèi)星的迎風(fēng)面積精確計算問題；楊成等[5]采用了基于陰影圖計算迎風(fēng)面積的方法；朱戰(zhàn)霞等[6]提出了以微元劃分思想建立航天器三維網(wǎng)格模型、使用射擊線掃描法求解迎風(fēng)面積的方法，計算結(jié)果與參考值誤差在10%以內(nèi)；劉衛(wèi)等[7]從小波分解的角度分析了衛(wèi)星阻力系數(shù)與空間環(huán)境參量變化的關(guān)系，為CD預(yù)報提供了方法；Mehta等[8]采用了一種響應(yīng)面模型方法來研究航天器的阻力系數(shù)，分析了球體、GRACE衛(wèi)星和哈勃望遠鏡以及國際空間站阻力系數(shù)的建模情況。

本文根據(jù)組合體和伴星的精密定軌需求，采用一種基于開放式圖形庫OpenGL和3ds格式三維模型的迎風(fēng)面積計算方法，開發(fā)工具軟件，給出正飛與連續(xù)偏航2種狀態(tài)下的迎風(fēng)面積計算方式，解決復(fù)雜結(jié)構(gòu)航天器有效迎風(fēng)面積計算困難的問題。

2 大氣阻力模型分析

除地球非球形攝動力之外，大氣阻力是影響低軌航天器軌道確定與預(yù)報精度的最主要的攝動力[9-10]。根據(jù)自由分子流理論，大氣分子與航天器表面撞擊后會發(fā)生阻留、漫散射、鏡反射等運動，大氣阻力的一般模型為式(1)[1]。

式中，F(xiàn)D為大氣阻力，ρ為航天器所處位置的大氣密度，v為航天器相對大氣的速度大小，uv單位速度矢量，CD為大氣阻力系數(shù)，Sa為有效迎風(fēng)面積，微積分計算公式如式(2)所示。

式中，d s為航天器表面在速度方向上的有效微平面元，θ為該微面元法向與速度方向的夾角。

由式(1)可知，ρ、Sa和CD這3個不確定性因素是大氣阻力難以精確建模的主要原因。大氣密度隨機性很強，各種影響因素比較復(fù)雜。目前使用的各種大氣密度模型有指數(shù)模型、改進的Harris-Priester模型、DTM模型等，大都屬于半經(jīng)驗?zāi)Ｐ?，誤差約為10%~20%。大氣密度不僅與高度有關(guān)，還與周日項、太陽活動、地磁活動、季節(jié)、經(jīng)緯度等因素密切相關(guān)[11-12]?？臻g環(huán)境參數(shù)主要用于計算大氣密度，基本都是通過觀測獲得并基于相應(yīng)的模型進行預(yù)報得到的。迎風(fēng)面積分析是從航天器自身角度出發(fā)，計算其在飛行方向上的有效截面積。通過準確獲取空間環(huán)境參數(shù)和計算航天器有效面積，有助于提高測定軌精度以及軌道預(yù)報的精度。航天器飛行過程中姿態(tài)在不斷變化，迎風(fēng)面積也隨之變化。在定軌中對迎風(fēng)面積的處理通常有2種方法，一是將航天器表面積的1/4作為有效迎風(fēng)面積值，二是將迎風(fēng)面積作為未知參數(shù)與軌道狀態(tài)參數(shù)一起進行估計。這2種方法都是近似處理，精度不高。實際應(yīng)用中為便于計算和提高估計的準確度，通常將航天器設(shè)計為對稱體或姿態(tài)容易掌握的小面質(zhì)比的構(gòu)型[4]。

鑒于CD的準確數(shù)值難以確定，因此在精密定軌中可先計算出目標迎風(fēng)面積，然后把CD作為未知量與衛(wèi)星運動狀態(tài)矢量一起解算，求解的CD一定程度上補償了由于大氣密度的模型誤差所導(dǎo)致的大氣阻力的誤差，能夠使動力學(xué)模型和觀測數(shù)據(jù)更好地擬合，從而得到更高的內(nèi)符合精度[13]。

3 有效迎風(fēng)面積計算方法

3.1 基于OpenGL的計算方法

天宮二號、神舟十一號組合體和伴星的三維模型分別如圖1所示，圖中oxyz為各自的本體坐標系，正飛時oz軸指向地心。從圖中可以看出，伴星結(jié)構(gòu)較為簡單，而天宮二號與神舟十一號結(jié)構(gòu)復(fù)雜，迎風(fēng)方向上各部件之間存在相互遮擋關(guān)系。

圖1 組合體和伴星三維模型Fig.1 3D models of complex and accompanying satellite

采用微元劃分思想計算目標迎風(fēng)面積時，需要3個步驟：目標三維建模與細分、讀取和處理模型、面積計算。首先需要將目標細分為許多個微小平面面元，并計算每個小面元的頂點坐標、法向矢量和材質(zhì)類型，判斷出有效迎風(fēng)的小面元，還要判斷面元的相互遮擋關(guān)系，將被遮擋的面元剔除掉，之后計算有效面元的面積，最后通過累加方式得到整個目標的迎風(fēng)面積。

上述過程中，所需目標模型可以通過三維建模獲得，每個面元對應(yīng)的θi也可通過數(shù)值計算獲得，但面元之間相互遮擋關(guān)系的判斷則比較復(fù)雜，需要多次循環(huán)計算和比較。通過軟件編程執(zhí)行循環(huán)計算和判斷，比較費時且精度較低，尤其對于復(fù)雜目標，面元劃分越細則面元數(shù)量越多，計算時間呈指數(shù)級增加。但若減少面元數(shù)量、增大面元面積，相互遮擋關(guān)系的判斷誤差則大幅增加，計算精度降低。

OpenGL提供了強大的圖形函數(shù)將大量數(shù)據(jù)轉(zhuǎn)換成圖形或圖像在屏幕上顯示出來并可進行交互處理，即快速可視化功能。此外OpenGL在深度緩存中為每個像素存儲了該像素與視點之間的深度值，可以通過啟用深度緩沖區(qū)和深度測試，利用計算機圖形加速卡在硬件層級迅速實現(xiàn)目標消隱。

為解決計算效率及精度與面元細分和遮擋判斷之間的矛盾，采用一種基于OpenGL的有效迎風(fēng)面積計算方法，利用OpenGL的可視化與消隱功能開發(fā)了工具軟件，能夠計算航天器在任意姿態(tài)下的有效迎風(fēng)面積，每次用時小于0.05 s，速度快、效率高，可用于飛行試驗任務(wù)中面積實時計算的場合，其基本流程如圖2所示。

1)根據(jù)工業(yè)部門提供的航天器模型文件，使用三維建模軟件3DSMAX進行模型格式轉(zhuǎn)換、編輯、加工和處理，用不同的材質(zhì)表示和區(qū)別航天器不同部件和載荷(艙體、太陽能電池帆板和天線等)，設(shè)定材質(zhì)對光源的散射特性是漫反射。

2)通過軟件編程讀取3ds格式目標模型文件。面積計算需要用到目標表面各點位置、法線方向、面元面積、面元所在的材質(zhì)種類等數(shù)據(jù)。3ds文件格式是一個由塊組成的層級結(jié)構(gòu)，基本塊即主塊包括版本、編輯信息和關(guān)鍵幀信息等，每個基本塊中的一級塊又包含了多種信息。需要讀取的模型數(shù)據(jù)包括目標對象材質(zhì)信息、頂點信息和面元信息等。本文將其讀入到自定義數(shù)據(jù)結(jié)構(gòu)之中，然后利用OpenGL繪制目標模型。

圖2 迎風(fēng)面積計算流程示意圖Fig.2 Flowchart of windward area calculation

3)根據(jù)當(dāng)前時刻目標姿態(tài)對三維模型進行旋轉(zhuǎn)變換，應(yīng)用OpenGL的可視化與消隱功能，采用正交投影模式將三維目標投影到二維屏幕并啟用深度緩沖區(qū)進行深度測試，并完成復(fù)雜目標的消隱和可見部位的顯示，提高了計算的準確度和速度。

4)采用Phong光照模型，通過設(shè)置合適的參數(shù)，讀取光照渲染處理后的幀緩存中的像素面元顏色值，得到各個面元的法向矢量、各面元與速度方向的夾角以及面元所代表的部件種類等參數(shù)信息，從而計算各面元迎風(fēng)面積，最后通過累加計算出整個目標的迎風(fēng)面積。

設(shè)置光源為方向漫反射白光源，強度為1。設(shè)第(i，j)個像素面元的單位法矢為n，面元材質(zhì)屬性為Cm=(R，G，B)，速度方向單位矢量為l，設(shè)θij為n和l的夾角，則顯示區(qū)域上該像素的顏色值為式(3)。

通過讀取目標模型可見部分每個像素的顏色值Iij及其組合方式，可計算出面元法向與速度方向的夾角θij以及面元材質(zhì)種類。設(shè)該像素對應(yīng)的目標面元的實際面積為s0，其有效迎風(fēng)面積為sij為式(4)。

根據(jù)式(2)，通過計算各面元迎風(fēng)面積sij的代數(shù)和，可得到整個目標的迎風(fēng)面積Sa，即式(5)。

為檢驗上述方法的準確性，以長方體為目標開展比對驗證實驗。如圖3所示，建立一個長方體模型，其長度為l，寬度為w，高度為h，長方體繞坐標系o-xyz的oy軸做順時針旋轉(zhuǎn)，旋轉(zhuǎn)角度用θ表示(對應(yīng)于航天器本體的偏航轉(zhuǎn)動)，速度方向與ox軸相同，圖中是θ=0°時的姿態(tài)情況。

圖3 長方體及其迎風(fēng)面積計算結(jié)果Fig.3 Calculation results of the cuboid and effective windward area

對于任一旋轉(zhuǎn)角θ，長方體有效迎風(fēng)面積Sa的理論值為式(6)。

當(dāng) l=1.5 m，w=1.0 m，h=2.0 m，θ從 0°以0.5°的間隔增加至180°時，根據(jù)式(6)計算的Sa理論值、本文方法計算的迎風(fēng)面積與理論值之間的偏差如圖3中曲線圖所示。

該長方體Sa的理論值范圍為2.0~3.547 m2，本文的計算結(jié)果與理論值的偏差小于0.005 m2，相對偏差小于0.2%，表明了本文方法的正確性。

3.2 正飛與偏航時迎風(fēng)面積計算方式

大型航天器正飛時以三軸穩(wěn)定的方式在軌運行，其偏航、俯仰、滾轉(zhuǎn)角基本為0，只有太陽電池板根據(jù)太陽與軌道面的夾角轉(zhuǎn)動。為補償大氣密度建模誤差，在精密定軌時使用3圈以上弧段，弧長6 h、9 h、12 h，則該時間段有效迎風(fēng)面積S0可認為是太陽翼轉(zhuǎn)動一圈的平均迎風(fēng)面積與本體最大截面積之和。航天器連續(xù)偏航時，其迎風(fēng)面積則難以用數(shù)學(xué)表達式來描述。

如圖4所示，設(shè)太陽翼的寬為w，高為h，面積A0=wh，初始狀態(tài)時太陽翼與速度風(fēng)向(迎風(fēng)方向)垂直，繞旋轉(zhuǎn)軸轉(zhuǎn)動一周總共需要轉(zhuǎn)動N次，轉(zhuǎn)到第i次(i=1，… N)時與初始狀態(tài)的夾角為θi，這時的有效迎風(fēng)面積Ai如式(7)所示。

圖4 太陽翼電池板轉(zhuǎn)動示意圖Fig.4 Schematic diagram of solar panel rotation

則太陽翼轉(zhuǎn)動N次的平均迎風(fēng)面積為式(8)。

當(dāng)N趨向于無窮大時，這時的S(N)即可認為是太陽翼轉(zhuǎn)動一周的迎風(fēng)面積S，如式(9)所示。

使用積分函數(shù)來表示上式，設(shè)x=2πi/N，則d x=2π/N，因此可得式(10)。

則航天器正飛狀態(tài)下的有效迎風(fēng)面積S0為式(11)。

式中，A0為航天器太陽翼(包括電池片和支架)的面積，S2為航天器艙體及外部載荷的截面積。

對于天宮二號與神舟十一號，單體正飛時可用式(11)直接估算其迎風(fēng)面積，而二者對接而成組合體后，正飛時由于兩個艙體尺寸結(jié)構(gòu)不同，且各自太陽翼存在前后遮擋關(guān)系，因此組合體的迎風(fēng)面積無法直接根據(jù)式(11)計算，需要使用本文的工具軟件采取不同的計算方式，具體如下：

1)正飛狀態(tài)時，保持航天器模型本體不動，太陽翼繞旋轉(zhuǎn)軸等間隔轉(zhuǎn)動0.05°或更小，累計轉(zhuǎn)動7200次，將各次計算的迎風(fēng)面積求平均值，即可認為是正飛狀態(tài)下的迎風(fēng)面積；

2)連續(xù)偏航時，需要讀取定軌所使用的測量弧段內(nèi)航天器偏航角和太陽翼轉(zhuǎn)動角度等姿態(tài)參數(shù)文件，計算各個時刻對應(yīng)的迎風(fēng)面積，然后取平均值作為該偏航弧段內(nèi)的有效迎風(fēng)面積。

4 實際定軌解算情況

在天宮二號與神舟十一號飛行任務(wù)中，利用各航天器的地基、天基測量數(shù)據(jù)并考慮軌控影響進行精密軌道確定，計算和預(yù)報各目標在過去、當(dāng)前和未來一段時間的運動狀態(tài)(包括星歷和姿態(tài)參數(shù))[11]，軌道確定位置誤差優(yōu)于10 m量級。

任務(wù)中測定軌采用MSIS90半經(jīng)驗大氣密度模型，主要使用當(dāng)日地磁指數(shù)日均值A(chǔ)p、前1天10.7 cm太陽輻射流量日均值F10.7和81天10.7 cm輻射流量平均值F10.7p等空間環(huán)境參數(shù)[14]。參數(shù)來源為每日國家空間科學(xué)中心提供的空間環(huán)境參數(shù)值及3天預(yù)報值。該模型在空間實驗室任務(wù)階段應(yīng)用效果良好[15]。

從10月23日07：31至 11月 15日 12：00約24天的時間里，伴星與組合體均以正飛的方式運行在高度約為388 km的軌道上，期間兩目標每天早8點的相對距離變化情況如圖5所示，相對距離范圍為0~280 km，但平均軌道高度相差在2 km以內(nèi)。

圖5 組合體與伴星相對距離變化情況Fig.5 Changes of relative distance between the complex and accompanying satellite

24天里的空間環(huán)境參數(shù)變化情況如圖6所示。除10月27日發(fā)生一次較為明顯的磁暴外，整個過程中空間環(huán)境比較平靜，太陽活動穩(wěn)定，F(xiàn)10.7變化不大，均值在77左右。

圖6 空間環(huán)境參數(shù)變化情況Fig.6 Changes of space environment parameters with time

在迎風(fēng)面積參數(shù)方面，最初計算時只考慮了艙體最大直徑和太陽翼電池板的尺寸，與實際情況有一定偏差。本文根據(jù)各目標的精細三維模型，綜合考慮了艙體及艙外載荷和零部件、太陽翼厚度及支架的尺寸等來計算迎風(fēng)面積。伴星和組合體精密定軌時使用了3整圈以上的相同的數(shù)據(jù)弧段，其空間環(huán)境參數(shù)一致。設(shè)使用原先方式計算的組合體正飛時的迎風(fēng)面積為M1、伴星面積為N1，使用本文方法計算得到組合體迎風(fēng)面積為M2、伴星的面積為N2，在10月23日至11月2日期間，每天1~2次的兩目標相同弧段精密定軌中，分別使用M1、N1求解得到的大氣阻力系數(shù)CDm1、CDn1及二者之差，分別使用 M2、N2求解得到的大氣阻力系數(shù)CDm2、CDn2及二者之差，具體如表1所示，表中序號表示不同的定軌歷元時間。

從表1中可以看出，使用原先面積參數(shù)時兩目標分別解CD值相差較大，差值基本在0.2~0.35，均值為0.26；使用本文方法計算出的面積參數(shù)后，兩目標CD相差約為0.01~0.07，均值為0.03。

從11月2日至11月15日期間，每天1~2次的兩目標同弧段精密定軌時，分別使用M2、N2求解得到的大氣阻力系數(shù)CDm2、CDn2及二者之差如表2所示，伴星求解CD均值為2.033，組合體解CD均值為2.026，兩目標CD之差小于0.01。

表1 伴星和組合體定軌求解C D結(jié)果對比(10月23日—11月2日)Table 1 The first group of C Dcomparisons of the complex and accompanying satellite(Oct.23-Nov.2)

表2 伴星和組合體定軌求解C D結(jié)果對比(11月2日—11月15日)Table 2 The second group of C D comparisons of the complex and accompanying satellite(Nov.2—Nov.15)

理論上相同軌道、相同表面材料、相同形狀的兩目標精密定軌求解的大氣阻力系數(shù)結(jié)果應(yīng)當(dāng)相同。由于組合體與伴星軌道高度接近、二者有效迎風(fēng)面為平面，且迎風(fēng)面材質(zhì)大部分是太陽能電池片，因此二者求解的大氣阻力系數(shù)結(jié)果應(yīng)當(dāng)相近。根據(jù)表1和表2可以看出，使用本文方法計算出的迎風(fēng)面積參數(shù)，兩目標精密定軌求解大氣阻力系數(shù)結(jié)果一致，與理論分析相符合。

5 結(jié)論

本文采用一種基于OpenGL和3ds格式三維模型的大型復(fù)雜結(jié)構(gòu)航天器迎風(fēng)面積計算方法，給出了正飛與連續(xù)偏航2種狀態(tài)下的迎風(fēng)面積計算方式；根據(jù)組合體與伴星共同正飛的24天里的精密軌道確定結(jié)果，求解得到的兩目標大氣阻力系數(shù)十分接近，二者之差約0.01~0.07。此外，兩目標定軌結(jié)果表明388 km軌道高度的CD平均值約為2.0，與傳統(tǒng)使用的CD值2.2有一定的差異，進行軌道預(yù)報時CD值使用2.0比2.2的精度更高。

與其他方法相比，本文方法和軟件具有自主知識產(chǎn)權(quán)并經(jīng)過了實際航天任務(wù)的數(shù)據(jù)驗證，重點針對大型航天器，目標三維模型與材質(zhì)模型的精細化程度更高，具有一定的工程實用價值，并且便于開展大量仿真分析實驗。后續(xù)將根據(jù)實際需求進行功能擴展和升級，為載人空間站任務(wù)航天器有效面積計算、所在軌道高度大氣阻力特性研究提供參考和支持。