航天器太陽光壓面積影響因素仿真分析

韓 意，陳 明，段成林，李翠蘭

（北京航天飛行控制中心，北京 100094）

引 言

航天器在軌過程中會(huì)受到太陽光壓的影響[1]，太陽光壓力是影響深空探測航天器軌道確定與預(yù)報(bào)精度的最主要的攝動(dòng)力，軌道動(dòng)力學(xué)模型誤差主要來源于太陽光壓模型。太陽光壓受到多種因素的影響，包括航天器自身特性、姿態(tài)變化及控制誤差和太陽活動(dòng)等，光壓面積是表征太陽光壓力的重要指標(biāo)，是航天器精密軌道確定和預(yù)報(bào)中所需的關(guān)鍵參數(shù)，直接影響定軌和軌道預(yù)報(bào)的精度[2-3]。國內(nèi)外學(xué)者針對太陽光壓問題開展了諸多研究，建立了多種光壓模型，主要包括分析型/物理型模型、經(jīng)驗(yàn)型模型和半經(jīng)驗(yàn)型模型3類。

分析型模型主要適用于新發(fā)射航天器，根據(jù)航天器的幾何形狀、大小、表面材料光學(xué)特性等物理特性進(jìn)行建模分析，具有清晰的物理內(nèi)涵。但航天器長期在軌運(yùn)行后光壓會(huì)發(fā)生變化，繼續(xù)使用先前模型會(huì)帶來精度不足的問題。最早提出的分析型模型是球模型，后來又發(fā)展了ROCK系列模型、T30模型和G2A模型等[4]。經(jīng)驗(yàn)型模型是根據(jù)航天器在軌運(yùn)行后的精密星歷等大量歷史數(shù)據(jù)統(tǒng)計(jì)計(jì)算建立的，精度較高，主要有Colombo模型、ECOM模型、擴(kuò)展CODE模型等，但此類方法需要依靠長期大量觀測數(shù)據(jù)做支撐，并且吸收了多種攝動(dòng)力的影響，物理內(nèi)涵缺乏，不利于單獨(dú)分析光壓的變化。半經(jīng)驗(yàn)?zāi)Ｐ徒Y(jié)合了分析型與經(jīng)驗(yàn)?zāi)Ｐ偷膬?yōu)點(diǎn)，在綜合考慮航天器發(fā)射前的自身基礎(chǔ)信息和屬性的同時(shí)，兼顧了發(fā)射后航天器在軌運(yùn)行狀態(tài)的運(yùn)動(dòng)信息，主要有JPL模型、Ad box-wing模型等[5-7]。

航天器的復(fù)雜性，主要體現(xiàn)在幾何結(jié)構(gòu)復(fù)雜、目標(biāo)表面材料種類多樣、材料光散射特性不同等方面[7]。在我國深空探測任務(wù)中，各航天器目前還沒有專用的太陽光壓模型，通常做法是簡化航天器的形狀和結(jié)構(gòu)組成、估算其在光照方向上的投影面積（橫截面積）及光壓面積。此種方法較為簡便，但未考慮部件之間相互遮擋關(guān)系以及表面不同材質(zhì)對太陽光的反射特性，與實(shí)際情況有較大差距。而要提高光壓模型精度，不僅需要考慮探測器不同部件在光照方向上的相互遮擋關(guān)系，還需要考慮表面不同材質(zhì)對太陽光的不同反射特性，遮擋判斷和處理計(jì)算過程相對比較復(fù)雜[8-9]。

本文針對實(shí)際任務(wù)特點(diǎn)進(jìn)行分析型光壓模型建模，了一種基于目標(biāo)特性的光壓面積快速計(jì)算方法，以仿真分析計(jì)算方式，側(cè)重于對航天器太陽光壓面積的影響因素（不同表面材料、不同尺寸、不同形狀）進(jìn)行初步分析和探討，嘗試分析目標(biāo)表面材料、尺寸和外形對目標(biāo)整體光壓面積的影響程度和一般性規(guī)律特點(diǎn)，有助于深空探測航天器太陽光壓的建模與分析，并對探測器研制和設(shè)計(jì)工作提供參考。

1 光壓面積計(jì)算方法

采用微面元的思想，將受照面劃分為各微小平面面元，然后以平面受照情況下的力學(xué)模型為基礎(chǔ)，分析整個(gè)航天器的太陽光壓[10]。設(shè)某一微小平面的法向單位矢量為n，太陽光線照在該平面上，入射角為θ，太陽光產(chǎn)生的力及方向如圖1所示。

圖1 平面受太陽光照壓力示意圖Fig. 1 Sketch map of lighting pressure on the tiny flat

在考慮熱輻射對受照平面表面力影響的情況下，該微平面受到的入射光產(chǎn)生的壓力分解到法線方向，經(jīng)過合成可以得到法向分力為[11]

其中：E為地球表面太陽輻射常數(shù)；c為光速；R1為面元到太陽的距離；R0為地球到太陽的距離；Sc為面元在入射方向上的投影面積，其值為Acosθ，A為面元的實(shí)際受照面積；γ為反射系數(shù)（無量綱）；β為鏡面系數(shù)（無量綱），漫反射系數(shù)等于γ（1–β）[12]。

用S表示該微平面的太陽光壓面積參數(shù)（下文簡稱光壓面積），即

用光壓系數(shù)因子k表示面元的光壓面積S與投影面積Sc的比值，即

顯然k與面元的光學(xué)特性以及面元的朝向有關(guān)。對于表面材質(zhì)有多種、三維立體的而非平面的航天器，設(shè)第i種材質(zhì)面元的反射系數(shù)為γi，鏡面系數(shù)為βi，系數(shù)因子k用k(γi,βi,θ)表示，根據(jù)微積分思想，可得到用積分形式表示的投影面積Sc為

光壓面積Sp可表示為

實(shí)際情況中航天器模型的表面難以用解析表達(dá)式來描述，且星體與帆板之間存在復(fù)雜的遮擋關(guān)系。為解決航天器光壓面積計(jì)算精度及效率與面元遮擋判斷之間的矛盾，本文采用了一種基于目標(biāo)特性的光壓面積計(jì)算方法，利用開放式圖形庫OpenGL開發(fā)了計(jì)算軟件，通過簡單幾何體光壓面積理論值與仿真值比對，驗(yàn)證了該方法的準(zhǔn)確性。本方法對目標(biāo)模型的消隱處理速度和幀緩存中像素信息讀取速度都很快，進(jìn)行某一姿態(tài)情況下的單次光壓面積計(jì)算用時(shí)小于0.05 s，可用于航天器飛行試驗(yàn)任務(wù)中光壓面積實(shí)時(shí)計(jì)算的場合。主要計(jì)算流程包括：

1）目標(biāo)三維建模：利用三維建模軟件3DMAX對復(fù)雜結(jié)構(gòu)航天器進(jìn)行幾何建模，通過修改和編輯含有材質(zhì)信息的目標(biāo)三維模型，輸出3ds文件格式；

2）讀取3ds模型：使用基于OpenGL開發(fā)的計(jì)算軟件，讀取處理好的3ds格式模型文件，主要讀取信息包括目標(biāo)表面各點(diǎn)的位置、外法線方向、面元面積、面元代表的材質(zhì)種類等，用以實(shí)現(xiàn)對模型的編輯、控制和繪制；

3）OpenGL消隱處理和計(jì)算：采用正交投影模式、根據(jù)目標(biāo)姿態(tài)進(jìn)行旋轉(zhuǎn)變換，啟用深度緩存進(jìn)行深度測試實(shí)現(xiàn)目標(biāo)消隱，讀取渲染處理后幀緩存中的像素面元信息，計(jì)算得到各個(gè)面元法向矢量、面元與入射光夾角以及面元材質(zhì)種類，進(jìn)而計(jì)算各面元光壓面積，最后計(jì)算出整個(gè)目標(biāo)的光壓面積。通過改變目標(biāo)姿態(tài)和光源位置即光照方向，可計(jì)算復(fù)雜目標(biāo)在不同條件下的光壓面積。

根據(jù)該計(jì)算方法，只需更換3ds格式航天器三維模型和相應(yīng)材質(zhì)的光學(xué)特性參數(shù)，即可快速分析不同航天器的光壓面積隨各種影響因素的變化情況，便于開展相關(guān)仿真實(shí)驗(yàn)。

2 目標(biāo)光壓面積影響因素分析

由于航天器通常由星體、一對或單個(gè)的較大的太陽能電池帆板構(gòu)成，星體通常為立方體或圓柱體，另外少量的衛(wèi)星是棱錐體，還有通訊天線和雷達(dá)等。因此，為分析光壓面積影響因素的一般規(guī)律，本文建立了立方體、圓柱體、圓錐體、球體四個(gè)凸面體目標(biāo)進(jìn)行光壓建模，分別計(jì)算各目標(biāo)在光照方向上的光壓面積Sp以及投影面積Sc的變化情況。

2.1 初始仿真條件

4個(gè)目標(biāo)體積都為1.0 m3，尺寸參數(shù)如表1所示。

定義計(jì)算坐標(biāo)系o-xyz及3個(gè)目標(biāo)的初始姿態(tài)如圖2所示，太陽光照方向與oz軸負(fù)向相同，設(shè)目標(biāo)繞x軸順時(shí)針旋轉(zhuǎn)的角度a為目標(biāo)姿態(tài)角。

表1 目標(biāo)尺寸參數(shù)Table 1 Targets' size parameters

圖2 目標(biāo)坐標(biāo)系及太陽光入射方向示意圖Fig. 2 Targets’ body coordinate systems and the lighting direction

分析中不考慮衛(wèi)星表面材料反射率的影響因素，如波長、表面粗糙度等，而是直接分析材質(zhì)反射率對目標(biāo)整體光壓面積的影響。假設(shè)目標(biāo)表面均為同一種材質(zhì)，反射系數(shù)用γ表示，鏡面系數(shù)用β表示，選取理想的、具有典型意義的完全漫反射、完全鏡面反射和綜合反射情況的3種材質(zhì)，如表2所示。

表2 表面材質(zhì)光學(xué)特性參數(shù)Table 2 Optical characteristic parameters of surface materials

假設(shè)材料的光學(xué)特性為完全吸收時(shí)，γ=β= 0，根據(jù)式（2）和式（3）可知，光壓面積與材料特性為完全漫反射時(shí)的值相同，因此本文不再單獨(dú)分析材料為完全吸收的情況。

設(shè)球體半徑為R，由于球中心對稱，故建立以球心為原點(diǎn)的坐標(biāo)系，如圖3所示。以球心到太陽的單位向量d的方向?yàn)閦軸正方向，取積分面元為dA，面元法向?yàn)閚，設(shè)n與d的夾角為θ，n與xoz平面的夾角為，n在xoz平面的投影n’與d的夾角為，定義n與+y軸同向時(shí)為90°，n’與+x軸同向時(shí)為90°，、取值范圍均為-90°～90°，則dA=R2cosφd?dφ，根據(jù)面元幾何關(guān)系為

圖3 球體光壓面積計(jì)算示意圖Fig. 3 Sketch map of calculating the sphere's SRP area

根據(jù)式（5）進(jìn)行積分，可得球體光壓面積S的與投影面積Sc= πR2的關(guān)系式為

從式（7）可知，理論上球體的光壓面積與表面材質(zhì)的反射率無關(guān)，比例系數(shù)k的值約為1.33。

使用工具軟件可計(jì)算出當(dāng)表面材質(zhì)不同時(shí)，半徑為0.620 3 m體積為1 m3的球體的光壓面積均為、投影面積和比值k如表3所示。從表3可知，使用工具軟件計(jì)算得到的球體光壓面積結(jié)果與理論分析結(jié)論相同。

表3 球體光壓面積計(jì)算結(jié)果Table 3 Calculation results of sphere SRP area

下文分別分析計(jì)算各目標(biāo)在不同表面材料、不同尺寸和不同形狀條件下的光壓面積變化情況。

2.2 表面材料對光壓面積的影響

在圖2中當(dāng)目標(biāo)繞ox軸旋轉(zhuǎn)一周即目標(biāo)姿態(tài)角a從0～360°時(shí)，立方體、圓柱、圓錐和球體4個(gè)目標(biāo)在垂直于太陽入射方向的平面上的投影面積（單位為：m2）變化曲線如圖4所示，從圖4中可以看出，圓錐體投影面積整體上最大，立方體的投影面積最小，圓柱體的居中。

設(shè)立方體、圓錐和圓柱體的光壓面積分別為Scube、Scone、Scder，設(shè)表面材質(zhì)分別為表2所示的完全漫反射、完全鏡面反射和綜合反射時(shí)目標(biāo)的光壓面積分別為Sdiff、Sspec、Smix。

當(dāng)姿態(tài)角a從0～360°時(shí)，3個(gè)目標(biāo)的光壓面積變化情況如圖5所示。從圖5可以得到以下結(jié)論：①對于立方體和圓柱體，在表面材質(zhì)分別為漫反射、鏡面反射和綜合情況時(shí)，整體上目標(biāo)光壓面積S的大小關(guān)系為Sspec>Smix>Sdiff，即鏡面反射材質(zhì)時(shí)目標(biāo)光壓面積最大，漫反射材質(zhì)時(shí)光壓面積較?。虎趫A錐體表面材質(zhì)為鏡面反射時(shí)，光壓面積的變化幅度大于漫反射和綜合情況時(shí)的變化幅度。

圖4 四種目標(biāo)的投影面積隨姿態(tài)角變化曲線Fig. 4 Four targets' projection areas changing with attitude angle a

圖5 目標(biāo)光壓面積隨姿態(tài)角a變化情況Fig. 5 Targets' SRP areas changing with attitude angle a

3個(gè)目標(biāo)的光壓系數(shù)k隨姿態(tài)角a從0～360°時(shí)的變化情況如圖6所示。3個(gè)目標(biāo)的光壓系數(shù)k的均值如表4所示。

從圖6和表4可知，鏡面反射材質(zhì)時(shí)，各目標(biāo)的光壓系數(shù)均值大于漫反射和綜合情況；整體上立方體的光壓系數(shù)最大，圓錐體的最小。

2.3 尺寸對光壓面積的影響

為分析目標(biāo)形狀相同但尺寸不同時(shí)的光壓面積變化情況，本文按照等比例關(guān)系新構(gòu)建了3個(gè)大尺寸目標(biāo)：立方體邊長2 m；圓柱體底面半徑1.128 4 m、高2 m；圓錐體底面半徑1.382 m、高4 m。這3個(gè)大目標(biāo)的尺寸是原來小目標(biāo)的2倍，體積是小目標(biāo)的8倍，都為8 m3。目標(biāo)初始姿態(tài)以及表面材質(zhì)類型都與前述情況相同。

圖6 目標(biāo)光壓系數(shù)k隨姿態(tài)角a變化情況Fig. 6 Targets' SRP coefficient k changing with attitude angle a

表4 不同目標(biāo)光壓系數(shù)計(jì)算結(jié)果Table 4 The mean value of targets' SRP coefficient k

當(dāng)姿態(tài)角a從0～360°、材質(zhì)為完全漫反射時(shí)，3種形狀的大尺寸目標(biāo)與相應(yīng)小尺寸目標(biāo)的光壓面積的比值變化曲線如圖7所示。

從圖7中可以看出，盡管形狀不同，但是大尺寸目標(biāo)的光壓面積與相應(yīng)小尺寸目標(biāo)的比值都是一致的，約為4。通過計(jì)算可知，材質(zhì)分別為綜合反射和鏡面反射時(shí)，大尺寸目標(biāo)與相應(yīng)小尺寸目標(biāo)的光壓面積比值也是4。

2.4 形狀對光壓面積的影響

當(dāng)目標(biāo)姿態(tài)角a從0～360°時(shí)，對于同一種材質(zhì)，不同目標(biāo)的光壓面積變化曲線如圖8所示。

圖9顯示了當(dāng)姿態(tài)角a從0～360°時(shí)，3種材質(zhì)條件下不同形狀目標(biāo)的光壓面積之差的變化情況，分別是①立方體與圓柱的光壓面積之差（Scube–Scder）；②圓錐與圓柱的光壓面積之差（Scone–Scder）；③圓錐與立方體的光壓面積之差（Scone–Scube）。

從圖8和圖9可以看出：①整體上在體積相同的情況下，對于同一種材質(zhì)，立方體和圓柱體相比，立方體的光壓面積最大，圓柱體的光壓面積較??；②圓錐體頂部錐尖朝向太陽附近時(shí)（對應(yīng)于姿態(tài)角60°～120°）光壓面積較小，其他情況下光壓面積相比立方體和圓柱體都要大。

3 結(jié) 論

圖8 目標(biāo)光壓面積隨姿態(tài)角a變化情況Fig. 8 Targets' SRP areas changing with attitude angle a

圖9 不同形狀光壓面積之差隨姿態(tài)角a變化情況Fig. 9 The difference of different shape targets’ SRP areas changing with a

開展航天器太陽光壓建模研究，對于提高深空探測任務(wù)的定軌與預(yù)報(bào)精度具有重要價(jià)值。本文采用一種基于分析型太陽光壓面積的計(jì)算方法，對表面材料光學(xué)特性、尺寸和形狀等影響因素進(jìn)行了仿真計(jì)算和分析，通過歸納概括，可以得到以下一般性結(jié)論：①對于同一個(gè)目標(biāo)，當(dāng)其表面材質(zhì)鏡面反射特性較強(qiáng)時(shí)，目標(biāo)光壓面積較大，而漫反射特性較強(qiáng)時(shí)光壓面積較??；②在目標(biāo)體積相同但形狀不同的情況下，立方體的光壓面積較大、投影面積較小、光壓系數(shù)較大，圓柱體的光壓面積較小，而圓錐體的投影面積較大；③平面的光壓系數(shù)整Scube–Scder體上高于曲面的光壓系數(shù)；④在形狀相同的情況下，若大目標(biāo)的尺寸是小目標(biāo)尺寸的k倍，則相應(yīng)地大目標(biāo)的光壓面積值是小目標(biāo)的k2倍，這一比值與材質(zhì)的反射特性無關(guān)。此外，經(jīng)分析可知目標(biāo)的光壓面積與其投影面積的比值k并非是一個(gè)定值，之前研究人員在工程應(yīng)用中計(jì)算目標(biāo)光壓面積時(shí)所用的簡便方法是用投影面積乘以一個(gè)系數(shù)，這種方法與實(shí)際情況有著一定的差距，結(jié)果不太準(zhǔn)確。

總體來說，目標(biāo)表面材料種類、目標(biāo)形狀和尺寸等是光壓面積的重要影響因素。在實(shí)際應(yīng)用中，如果要降低目標(biāo)的整體光壓面積、減小光壓攝動(dòng)影響，根據(jù)本文的分析結(jié)果可知，一方面可以選擇漫反射特性較強(qiáng)、反射率較低的表面材料；另一方面可以縮小尺寸，將目標(biāo)設(shè)計(jì)為圓柱體。本文只得到了各種影響因素對目標(biāo)整體光壓面積的影響的一般性的、定性的分析結(jié)論，后續(xù)將進(jìn)一步開展各種影響因素的量化分析研究，以期獲得對工程應(yīng)用更具有指導(dǎo)意義的有益結(jié)論和啟示。