空間目標(biāo)地基觀測紅外輻射特性研究

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

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

1 引言

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

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

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

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

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

2 計算方法

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

2.1 仿真算法流程

首先，需要對目標(biāo)模型進(jìn)行有限元建模，采用非結(jié)構(gòu)四面體網(wǎng)格對空間目標(biāo)進(jìn)行離散，光線的步入和追蹤均以網(wǎng)格為單位進(jìn)行，網(wǎng)格結(jié)構(gòu)也是信息存儲結(jié)構(gòu)的基礎(chǔ)。程序中網(wǎng)格元、面元和節(jié)點(diǎn)分別定義為相應(yīng)的結(jié)構(gòu)體，并通過指針或序號相互映射。軟件通過讀取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ù)目的光線，并通過計算光線是否被遮擋，以及光線方向和太陽矢量、地球矢量的關(guān)系，匯總所有光線計算結(jié)果，得到目標(biāo)接收到的軌道外熱流。此處的外熱流包括太陽輻射、地球輻射及地球反照輻射，忽略其他天體的影響。

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

最后，匯總各表面面元的參數(shù)計算結(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 反射太陽輻射

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

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

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

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

2.4 反射地球輻射

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

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

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

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

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

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

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

3.1 表面反射特性

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

其中，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

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

對于采用BRDF 計算的表面表示為：

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

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

本文計算空間目標(biāo)輻射特性采用的模型主要參考風(fēng)云一號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è)固定安裝。帆板所在平面與軌道面平行，以保證在整個軌道周期內(nèi)太陽光照角度恒定[14]。軌道坐標(biāo)系定義為，+X指向飛行方向，+Z指向地心，+Y符合右手定則。仿真軟件運(yùn)行時需要輸入衛(wèi)星軌道信息（半長軸、偏心率、傾角、升交點(diǎn)赤經(jīng)、近地點(diǎn)幅角、真近點(diǎn)角）。軟件內(nèi)部完成遞推后。認(rèn)為本體系與軌道系重合。仿真中使用的6 個軌道參數(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)地方時約為8:34，太陽帆板電池片朝向?yàn)?Y平面。計算使用的模型如圖2 所示，使用的網(wǎng)格如圖3 所示，計算域?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 表面溫度計算結(jié)果

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

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

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

結(jié)合表面材料參數(shù)，計算得到的表面溫度，如圖6（彩圖見期刊電子版）所示。此處的計算結(jié)果僅考慮穩(wěn)態(tài)解。假設(shè)本體各表面間相互絕熱，內(nèi)熱源等效為各本體表面受到的熱流，設(shè)為50 W/m2。帆板導(dǎo)熱系數(shù)設(shè)置為 1.7 W/(m·K)。可以看出，除+Y表面外，其余各表面溫度在一個軌道周期內(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)度計算結(jié)果

進(jìn)一步地，計算目標(biāo)表面在3～5 μm 和8～14 μm的紅外輻射強(qiáng)度，如圖7～8（彩圖見期刊電子版）所示。此處選擇的位置點(diǎn)為陽照區(qū)過赤道時刻，方位角定義為繞Z軸逆時針旋轉(zhuǎn)方向，0 點(diǎn)位于+X軸，天頂角定義為繞Y軸逆時針旋轉(zhuǎn)方向，0 點(diǎn)位于+Z軸，探測器入瞳法線指向目標(biāo)，計算步長為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°兩個峰值，與文獻(xiàn)[4]分布趨勢一致。這主要是因?yàn)榉宸础⒄齼擅娴臏囟染^高，而空間目標(biāo)反射的太陽輻射主要在可見光波段，地球輻射相對較小，因此，目標(biāo)自身溫度產(chǎn)生的紅外輻射占比較高。

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

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

地基探測器位置選擇為北緯43.8°，東經(jīng)125.3°，高232 m 位置，設(shè)定最小俯仰角為20°后得到升軌過站弧段（UTC 時間）：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 中，序號1、2、3 分別代表升軌弧段的進(jìn)站、中點(diǎn)和出站；序號4、5、6 代表降軌弧段的進(jìn)站、中點(diǎn)及出站。地面站采用北東地坐標(biāo)系。

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

計算得到在3～5 μm 和8～14 μm 波段的紅外光譜輻射強(qiáng)度分布，如圖9～圖10（彩圖見期刊電子版）所示，光譜分辨率為0.02 μm。可以看出兩個弧段都是中間位置處紅外光譜輻射強(qiáng)度最大，升軌進(jìn)站時的紅外光譜輻射強(qiáng)度略低于出站，降軌進(jìn)站時的紅外光譜輻射強(qiáng)度高于出站。在升軌弧段對于觀測貢獻(xiàn)較大的+Z表面溫度逐漸升高，其他可探測表面溫度變化較小。而降軌弧段僅-Z表面溫度變化較大，但其對于地基觀測沒有貢獻(xiàn)，進(jìn)、出站時刻的紅外光譜輻射強(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 可以看出，升軌弧段俯仰角較大，接近過頂，對于溫度較高的-Y面帆板和衛(wèi)星本體在探測器入瞳處的有效面積較小。光譜輻射強(qiáng)度最大值在中點(diǎn)位置處的9～10 μm 波段內(nèi)產(chǎn)生，約為23 W/(sr·μm)，進(jìn)、出站時刻最大值分別約為18 W/(sr·μm)和19 W/(sr·μm)。降軌弧段俯仰角較小，且探測器方向在衛(wèi)星本體的-Y側(cè)，探測器入瞳處的有效面積較大，因此光譜輻射強(qiáng)度較大。中間位置處最大值約為113 W/(sr·μm)，進(jìn)、出站時刻最大值分別約為85 W/(sr·μm)和76 W/(sr·μm)。

進(jìn)一步地，考慮到大氣吸收及觀測路徑上大氣的自身輻射，使用Modtran 軟件計算大氣透過率和沿程輻射值。計算條件設(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

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

5 結(jié)論

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