鄭鴻儒，馬 巖 ，張 帥，陳亞濤

（1.北京跟蹤與通信技術(shù)研究所,北京 100094；2.長光衛(wèi)星技術(shù)股份有限公司,吉林 長春 130000；3.北京航空航天大學(xué) 宇航學(xué)院,北京 100191）

1 引言

近年來，隨著衛(wèi)星技術(shù)的不斷發(fā)展，各國對(duì)空間資源的爭奪越來越激烈。空間目標(biāo)的探測、識(shí)別和監(jiān)視對(duì)于國防的意義日益凸顯[1,2]。開展空間目標(biāo)的紅外輻射特性研究對(duì)于敵方目標(biāo)探測和我方紅外隱身都具有重要意義。

國內(nèi)外空間目標(biāo)紅外輻射特性研究主要包括目標(biāo)表面溫度場計(jì)算、目標(biāo)紅外輻射計(jì)算，以及通過地基或天基平臺(tái)的目標(biāo)觀測等。目標(biāo)表面溫度計(jì)算主要是通過積分法或蒙特卡洛法等獲得表面熱流，然后結(jié)合材料屬性迭代求解表面溫度。獲得表面溫度后，結(jié)合表面發(fā)射和反射特性即可計(jì)算出目標(biāo)紅外輻射值。進(jìn)而可以得到地基或天基平臺(tái)的可觀測性以及目標(biāo)光譜輻射特性。其可為分析目標(biāo)工作狀態(tài)，反演目標(biāo)表面材料屬性提供參考。

劉建斌等人[3]將目標(biāo)分別視為朗伯表面和隨機(jī)粗糙表面，計(jì)算得到了風(fēng)云一號(hào)衛(wèi)星在地面觀測位置處的光照度。申文濤等人[4]采用隨機(jī)起伏表面算法模擬材料表面，并結(jié)合雙向反射分布函數(shù)（BRDF）計(jì)算了衛(wèi)星在可見光和紅外譜段的輻射特性，并與地面衛(wèi)星模型試驗(yàn)進(jìn)行了對(duì)比。楊帆等人[5]計(jì)算了探測器接收到的紅外輻射光譜，采用最小二乘算法反演了目標(biāo)各組件的有效面積和表面溫度。李文豪等人[6]分析了散熱面在不同功率下在距離5 km 處入瞳的中波和長波輻射照度。谷牧等人[7]考慮了地基探測器的可見范圍及大氣衰減的影響，分析了地基探測器入瞳處影響紅外光譜輻射照度的主要因素。孫成明等人[8]對(duì)風(fēng)云一號(hào)衛(wèi)星通過有限元法進(jìn)行了建模，并與縮比模型散射輻照度的實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析，結(jié)果基本一致。Skinner 等人[9-10]對(duì)地球同步軌道衛(wèi)星的紅外光譜進(jìn)行了長期觀測和分析，獲得了等效面積和等效溫度等信息。

然而，以上研究要么對(duì)目標(biāo)的建模過程介紹較少，要么采用結(jié)構(gòu)化面網(wǎng)格進(jìn)行計(jì)算，當(dāng)空間目標(biāo)結(jié)構(gòu)比較復(fù)雜時(shí)建模及處理遮擋的難度較大。其次，在計(jì)算軌道外熱流時(shí)，通常將太陽光線與表面夾角的余弦值作為太陽輻射角系數(shù)，但該方法僅適用于凸表面。最后，以往研究中對(duì)于目標(biāo)紅外輻射研究較多，對(duì)考慮大氣衰減的地基探測場景研究較少。

基于此，本文基于非結(jié)構(gòu)化四面體網(wǎng)格開展有限元建模，編寫空間目標(biāo)紅外輻射特性仿真程序，采用蒙特卡洛算法追蹤光線，計(jì)算軌道外熱流。結(jié)合表面材料屬性計(jì)算得到目標(biāo)表面溫度和輻射特性。進(jìn)一步，考慮大氣衰減效應(yīng)，對(duì)目標(biāo)在地基探測器入瞳處接收到的紅外光譜輻射特性進(jìn)行研究，給出升軌和降軌過程中典型位置的紅外輻射光譜特性。

2 計(jì)算方法

本文采用有限元方法開展了空間目標(biāo)紅外光譜特性研究，探測器接收到的紅外輻射主要有目標(biāo)自身輻射、目標(biāo)反射太陽輻射、目標(biāo)反射地球輻射和目標(biāo)反射地球反照輻射[11]。

2.1 仿真算法流程

首先，需要對(duì)目標(biāo)模型進(jìn)行有限元建模，采用非結(jié)構(gòu)四面體網(wǎng)格對(duì)空間目標(biāo)進(jìn)行離散，光線的步入和追蹤均以網(wǎng)格為單位進(jìn)行，網(wǎng)格結(jié)構(gòu)也是信息存儲(chǔ)結(jié)構(gòu)的基礎(chǔ)。程序中網(wǎng)格元、面元和節(jié)點(diǎn)分別定義為相應(yīng)的結(jié)構(gòu)體，并通過指針或序號(hào)相互映射。軟件通過讀取CGNS（CFD General Notation System）格式網(wǎng)格文件，自主建立網(wǎng)格-面元的映射關(guān)系，判斷相鄰網(wǎng)格和邊界網(wǎng)格等，并將空間目標(biāo)表面溫度和材料屬性等參數(shù)離散到各表面中心點(diǎn)。

然后，通過矢量坐標(biāo)變換得到各面元中心點(diǎn)坐標(biāo)在目標(biāo)本體坐標(biāo)系下與地球、太陽的位置關(guān)系。通過蒙特卡洛方法，由目標(biāo)各表面面元發(fā)出給定數(shù)目的光線，并通過計(jì)算光線是否被遮擋，以及光線方向和太陽矢量、地球矢量的關(guān)系，匯總所有光線計(jì)算結(jié)果，得到目標(biāo)接收到的軌道外熱流。此處的外熱流包括太陽輻射、地球輻射及地球反照輻射，忽略其他天體的影響。

詳細(xì)計(jì)算表面溫度特性十分復(fù)雜，本文在計(jì)算時(shí)忽略了星體各表面面元之間的導(dǎo)熱和輻射換熱，僅考慮太陽帆板正反面的導(dǎo)熱。使用牛頓迭代法求解表面熱平衡方程，并根據(jù)溫度特性計(jì)算紅外輻射量。在計(jì)算軌道外熱流、各表面對(duì)給定探測器入瞳方向的輻射強(qiáng)度的過程中使用Open-MP 并行加速算法進(jìn)行加速。

最后，匯總各表面面元的參數(shù)計(jì)算結(jié)果，得到目標(biāo)溫度及紅外輻射特性。

2.2 自身輻射

表面自身輻射由表面溫度和表面材料發(fā)射率決定，在獲得各表面溫度后，目標(biāo)自身紅外光譜輻射出射度可由下式表示：

其中：c1為第一輻射常量，其值為3.742×10-16W·m2；c2為第二輻射常量，其值為1.438 8×10-2m·K；λ為波長，單位為 μm；Ti是面元表面溫度；εi(λ)是面元光譜發(fā)射率，根據(jù)文獻(xiàn)[5]知，在紅外譜段發(fā)射率基本不變，本文設(shè)置為常值。

2.3 反射太陽輻射

計(jì)算目標(biāo)表面反射太陽輻射時(shí)，通常將太陽視為一個(gè)溫度為5 900 K 的黑體，由普朗克公式得到的光譜輻射出射度為：

其中，Ts是太陽溫度。則太陽在面元上的光譜輻射強(qiáng)度為：

其中：Rs是太陽半徑；Ds是太陽-地球之間的距離；?s為太陽輻射角系數(shù)。面元反射太陽輻射光譜輻射出射度為：

其中，ρs為表面材料對(duì)太陽輻射的反射率。

2.4 反射地球輻射

面元反射地球輻射包括地球紅外輻射和地球反照太陽輻射。地球紅外輻射可以將地球視為一個(gè)溫度為255 K 的黑體，由普朗克公式得到光譜輻射出射度為：

其中，Te是地球溫度。則地球自身輻射在面元上的光譜輻射強(qiáng)度為：

式中Rsat為衛(wèi)星與地心的距離，Re是地球半徑，?e是地球紅外輻射角系數(shù)。地球反射輻射在面元上的光譜輻射強(qiáng)度為：

式中 ρa(bǔ)為地球?qū)μ栞椛涞钠骄瓷渎剩?a是地球反照輻射角系數(shù)。

面元反射地球輻射光譜輻射出射度為：

其中，ρe為表面材料對(duì)地球紅外輻射的反射率。

3 表面屬性及模型參數(shù)

3.1 表面反射特性

對(duì)于衛(wèi)星表面反射參數(shù)設(shè)置，本文將太陽帆板背陽面設(shè)置為漫反射，根據(jù)朗伯余弦定律計(jì)算。對(duì)于鏡面屬性較強(qiáng)的太陽帆板電池片表面和多層包覆的衛(wèi)星本體表面采用雙向反射分布函數(shù)（BRDF）進(jìn)行計(jì)算。由于BRDF 的模型較多，計(jì)算過程也比較復(fù)雜，本文選用文獻(xiàn)[12]中介紹的改進(jìn)Phong 模型進(jìn)行計(jì)算，選用GaSa 材料的模型參數(shù)計(jì)算帆板正面電池片反射特性，選用聚酰亞胺薄膜的模型參數(shù)計(jì)算衛(wèi)星本體反射。計(jì)算方法如下式所示：

其中，fr(θin;θr,φ) 為表面的BRDF，θin、θr、φ分別為入射天頂角、觀測天頂角、觀測方位角；ρd為材料的漫反射系數(shù)；ρsp為鏡面反射系數(shù)；α為鏡向指數(shù)；β為觀測方向與鏡反射方向的夾角，具體定義可見參考文獻(xiàn)[12]。目標(biāo)表面與探測器入瞳的位置關(guān)系示意如圖1 所示。

圖1 目標(biāo)表面與探測器入瞳的位置關(guān)系示意圖Fig.1 Schematic diagram of the position relationship between the target surface and the probe's entrance pupil

對(duì)于地球紅外輻射和地球反照輻射，由于每個(gè)方向基本都有能量貢獻(xiàn)，因此本文只對(duì)太陽直射采用BRDF 計(jì)算，其他輻射采用漫反射計(jì)算面元的反射紅外輻射強(qiáng)度[13]。則目標(biāo)在探測器方向上的紅外光譜輻射強(qiáng)度對(duì)于漫反射表面表示為：

對(duì)于采用BRDF 計(jì)算的表面表示為：

其中，i是目標(biāo)表面離散的面元序號(hào)，N是面元總個(gè)數(shù)。

3.2 模型及參數(shù)設(shè)置

本文計(jì)算空間目標(biāo)輻射特性采用的模型主要參考風(fēng)云一號(hào)C 衛(wèi)星（FY-1C），衛(wèi)星本體尺寸為：1.42 m×1.42 m×1.2 m，太陽帆板位于兩側(cè)，尺寸為4.5 m×0.2 m×1.2 m。衛(wèi)星在太陽同步軌道運(yùn)行。太陽帆板沿飛行方向前后兩側(cè)固定安裝。帆板所在平面與軌道面平行，以保證在整個(gè)軌道周期內(nèi)太陽光照角度恒定[14]。軌道坐標(biāo)系定義為，+X指向飛行方向，+Z指向地心，+Y符合右手定則。仿真軟件運(yùn)行時(shí)需要輸入衛(wèi)星軌道信息（半長軸、偏心率、傾角、升交點(diǎn)赤經(jīng)、近地點(diǎn)幅角、真近點(diǎn)角）。軟件內(nèi)部完成遞推后。認(rèn)為本體系與軌道系重合。仿真中使用的6 個(gè)軌道參數(shù)如表1 所示。表面參數(shù)設(shè)置如表2 所示。衛(wèi)星本體設(shè)置為多層材料。

表1 軌道參數(shù)Tab.1 Orbit parameters

表2 表面材料熱參數(shù)Tab.2 Thermal parameters of surface material

衛(wèi)星運(yùn)行軌道降交點(diǎn)地方時(shí)約為8:34，太陽帆板電池片朝向?yàn)?Y平面。計(jì)算使用的模型如圖2 所示，使用的網(wǎng)格如圖3 所示，計(jì)算域?yàn)榘ㄐl(wèi)星本體在內(nèi)的13 m×3 m×3m 的區(qū)域，總四面體網(wǎng)格數(shù)約為13.5 萬，衛(wèi)星目標(biāo)表面面元網(wǎng)格數(shù)約1.4 萬。

圖2 仿真算例中使用的簡化模型Fig.2 The simplified model used in simulation

圖3 仿真算例中的網(wǎng)格劃分Fig.3 Grid division in simulation cases

4 分析與討論

4.1 表面溫度計(jì)算結(jié)果

首先，為了驗(yàn)證算法的準(zhǔn)確性，根據(jù)前述的仿真條件，計(jì)算得到不考慮遮擋影響的正六面體目標(biāo)在一個(gè)軌道周期內(nèi)各表面受到的軌道外熱流，并將其與商業(yè)軟件的計(jì)算結(jié)果進(jìn)行對(duì)比，如圖4（彩圖見期刊電子版）所示。其中，線條為本文軟件計(jì)算結(jié)果，符號(hào)為對(duì)應(yīng)表面的商業(yè)軟件計(jì)算結(jié)果，計(jì)算起始時(shí)刻為2000-03-21 04:58:00?？梢钥闯?，本軟件的軌道外熱流計(jì)算結(jié)果與商業(yè)軟件相近，誤差較小。

圖4 不考慮遮擋效應(yīng)時(shí)本文計(jì)算得到的軌道外熱流與商用軟件的結(jié)果對(duì)比Fig.4 Comparison of calculation results of orbit external heat flow by proposed method and commercial software without taking into account the occulusion effect

采用圖3 給出的網(wǎng)格模型計(jì)算考慮遮擋效應(yīng)的軌道外熱流分布，如圖5（彩圖見期刊電子版）所示。可以看出，由于目標(biāo)三軸穩(wěn)定對(duì)地，且帆板位置固定不變，故本體-Y面與帆板電池片受到的熱流最大，為1 200 W/m2左右，且在陽照區(qū)和地影區(qū)都比較穩(wěn)定。其中，+X和-X側(cè)帆板表面在陽照區(qū)相對(duì)于圖4 存在熱流波動(dòng)，這是由于衛(wèi)星本體對(duì)陽光的遮擋造成的。-Z表面在陽照區(qū)、+Z表面在衛(wèi)星星下點(diǎn)進(jìn)出地影區(qū)時(shí)存在被陽光直射區(qū)間，因此，熱流隨時(shí)間變化較大，其他表面受到的熱流變化幅度較小。

圖5 考慮遮擋效應(yīng)時(shí)軌道外熱流計(jì)算結(jié)果Fig.5 Calculation results of orbit external heat flow with taking into account the occlusion effect

結(jié)合表面材料參數(shù)，計(jì)算得到的表面溫度，如圖6（彩圖見期刊電子版）所示。此處的計(jì)算結(jié)果僅考慮穩(wěn)態(tài)解。假設(shè)本體各表面間相互絕熱，內(nèi)熱源等效為各本體表面受到的熱流，設(shè)為50 W/m2。帆板導(dǎo)熱系數(shù)設(shè)置為 1.7 W/(m·K)?？梢钥闯?，除+Y表面外，其余各表面溫度在一個(gè)軌道周期內(nèi)變化較大。溫度最高的位置在帆板上，溫度在311～330 K 范圍內(nèi)變化，本體-Y表面陽照區(qū)溫度約為316 K。帆板在地影區(qū)和陽照區(qū)溫差最大，變化幅度約為150 K。其他表面溫度變化趨勢與熱流變化基本相同。

圖6 各表面溫度Fig.6 Temperature of each surface

4.2 紅外輻射強(qiáng)度計(jì)算結(jié)果

進(jìn)一步地，計(jì)算目標(biāo)表面在3～5 μm 和8～14 μm的紅外輻射強(qiáng)度，如圖7～8（彩圖見期刊電子版）所示。此處選擇的位置點(diǎn)為陽照區(qū)過赤道時(shí)刻，方位角定義為繞Z軸逆時(shí)針旋轉(zhuǎn)方向，0 點(diǎn)位于+X軸，天頂角定義為繞Y軸逆時(shí)針旋轉(zhuǎn)方向，0 點(diǎn)位于+Z軸，探測器入瞳法線指向目標(biāo)，計(jì)算步長為10°。

圖7 3～5 μm 波段紅外輻射強(qiáng)度Fig.7 Infrared radiation intensity in 3～5 μm band

圖8 8～14 μm 波段紅外輻射強(qiáng)度Fig.8 Infrared radiation intensity in 8～14 μm band

從圖中可以看出，在3～5 μm 中波和8～14 μm長波波段的輻射強(qiáng)度分布都存在天頂角90°、方位角90°和270°兩個(gè)峰值，與文獻(xiàn)[4]分布趨勢一致。這主要是因?yàn)榉宸?、正兩面的溫度均較高，而空間目標(biāo)反射的太陽輻射主要在可見光波段，地球輻射相對(duì)較小，因此，目標(biāo)自身溫度產(chǎn)生的紅外輻射占比較高。

對(duì)于中波波段，輻射強(qiáng)度最大值在45 W/sr左右。對(duì)于長波波段，輻射強(qiáng)度最大值在770 W/sr 左右。兩個(gè)波段在天頂角90°、方位角90°，即本體+Y軸方向強(qiáng)度稍高。這是因?yàn)榉逭疵鏈囟认嘟那闆r下，太陽帆板背板的發(fā)射率設(shè)置稍高于正面電池片。

4.3 地基探測紅外計(jì)算結(jié)果

地基探測器位置選擇為北緯43.8°，東經(jīng)125.3°，高232 m 位置，設(shè)定最小俯仰角為20°后得到升軌過站弧段（UTC 時(shí)間）：2000-03-21 11:26:16 至11:34:13，降軌過站弧段：2000-03-21 23:48:26 至23:54:58。分別選取弧段的起始、中點(diǎn)和結(jié)束位置為觀測點(diǎn)，兩種情況下的各位置點(diǎn)方位、俯仰角以及距離如表3 所示。表3 中，序號(hào)1、2、3 分別代表升軌弧段的進(jìn)站、中點(diǎn)和出站；序號(hào)4、5、6 代表降軌弧段的進(jìn)站、中點(diǎn)及出站。地面站采用北東地坐標(biāo)系。

表3 不同位置點(diǎn)處的方位角、俯仰角和距離Tab.3 Azimuth,elevation,and range at different positions

計(jì)算得到在3～5 μm 和8～14 μm 波段的紅外光譜輻射強(qiáng)度分布，如圖9～圖10（彩圖見期刊電子版）所示，光譜分辨率為0.02 μm?？梢钥闯鰞蓚€(gè)弧段都是中間位置處紅外光譜輻射強(qiáng)度最大，升軌進(jìn)站時(shí)的紅外光譜輻射強(qiáng)度略低于出站，降軌進(jìn)站時(shí)的紅外光譜輻射強(qiáng)度高于出站。在升軌弧段對(duì)于觀測貢獻(xiàn)較大的+Z表面溫度逐漸升高，其他可探測表面溫度變化較小。而降軌弧段僅-Z表面溫度變化較大，但其對(duì)于地基觀測沒有貢獻(xiàn)，進(jìn)、出站時(shí)刻的紅外光譜輻射強(qiáng)度差別主要由地基探測位置可見的有效面積決定。

圖9 升軌觀測紅外光譜輻射強(qiáng)度Fig.9 Infrared radiation spectral intensity of ascending detection

圖10 降軌觀測紅外光譜輻射強(qiáng)度Fig.10 Infrared radiation spectral intensity of descending detection

從表3 可以看出，升軌弧段俯仰角較大，接近過頂，對(duì)于溫度較高的-Y面帆板和衛(wèi)星本體在探測器入瞳處的有效面積較小。光譜輻射強(qiáng)度最大值在中點(diǎn)位置處的9～10 μm 波段內(nèi)產(chǎn)生，約為23 W/(sr·μm)，進(jìn)、出站時(shí)刻最大值分別約為18 W/(sr·μm)和19 W/(sr·μm)。降軌弧段俯仰角較小，且探測器方向在衛(wèi)星本體的-Y側(cè)，探測器入瞳處的有效面積較大，因此光譜輻射強(qiáng)度較大。中間位置處最大值約為113 W/(sr·μm)，進(jìn)、出站時(shí)刻最大值分別約為85 W/(sr·μm)和76 W/(sr·μm)。

進(jìn)一步地，考慮到大氣吸收及觀測路徑上大氣的自身輻射，使用Modtran 軟件計(jì)算大氣透過率和沿程輻射值。計(jì)算條件設(shè)置為中緯度地區(qū)冬季大氣，得到探測器入瞳處接收到的疊加大氣衰減效應(yīng)和沿程輻射的空間目標(biāo)紅外光譜輻射強(qiáng)度如圖11～圖12（彩圖見期刊電子版）所示。

圖11 考慮大氣影響的升軌觀測紅外光譜輻射強(qiáng)度Fig.11 Infrared radiation spectral intensity of ascending detection with the atmosphere impact

圖12 考慮大氣影響的降軌觀測紅外光譜輻射強(qiáng)度Fig.12 Infrared radiation spectral intensity of descending detection with the atmosphere impact

可以看出，兩種觀測情況下大氣影響均較大。對(duì)3～5 μm 中波波段影響比8～14 μm 長波波段大，光譜強(qiáng)度降低比例較多。對(duì)于兩種觀測條件，由于進(jìn)站和出站時(shí)的俯仰角較小，沿程損失比中間時(shí)刻大。在4.3 μm 和14.7 μm 的CO2吸收帶附近，以及9.5 μm 的O3吸收帶附近損失較大，在進(jìn)行地基紅外探測時(shí)應(yīng)避免使用附近波段。同時(shí)也看到，地基觀測無法避免地會(huì)受到大氣衰減和背景輻射干擾，導(dǎo)致很難準(zhǔn)確反演目標(biāo)信息，因此，聯(lián)合天基平臺(tái)開展空間目標(biāo)聯(lián)合探測很有必要。

5 結(jié)論

本文針對(duì)空間目標(biāo)的紅外輻射特性，基于有限元方法編寫了仿真計(jì)算軟件，計(jì)算得到了目標(biāo)在太陽同步軌道上接受的外熱流，結(jié)合表面材料屬性，計(jì)算得到了軌道周期內(nèi)的表面溫度變化。進(jìn)一步地，結(jié)合改進(jìn)的Phong 模型和漫反射模型，計(jì)算了目標(biāo)在中波波段和長波波段的紅外輻射強(qiáng)度空間分布。結(jié)果表明，對(duì)地三軸穩(wěn)定太陽同步衛(wèi)星采用沿飛行方向前后布置固定太陽帆板時(shí)，陽照區(qū)帆板表面溫度較高，約在311～330 K 范圍內(nèi)變化，是紅外觀測的主要輻射來源。采用8～14 μm 波段進(jìn)行探測時(shí)，紅外輻射強(qiáng)度最大值在770 W/sr 左右，遠(yuǎn)大于3～5 μm 波段。最后，基于地基探測場景，使用Modtran 軟件計(jì)算得到了在不用觀測路徑上的大氣衰減和沿程輻射，計(jì)算得到空間目標(biāo)在陽照區(qū)和地影區(qū)觀測場景中的紅外光譜輻射特性。本文研究結(jié)果可為空間目標(biāo)信息反演及探測波段優(yōu)選提供助力。