空間目標軌道外熱流計算及輻射特性研究

鄭鴻儒，馬 巖 ，張 帥，王建超，曲友陽

（1.北京跟蹤與通信技術研究所,北京 100094；2.長光衛星技術股份有限公司,吉林 長春 130000）

1 引言

目前各國對空間資源的爭奪愈演愈烈，空間態勢感知技術對國防的戰略意義日益凸顯，如何對空間目標進行有效監視成為亟需解決的問題之一?？臻g目標由于距離太遠，往往在探測器上成像只有幾個像素，可提取的目標信息十分有限[1]。為了維持載荷正常工作，空間目標不可避免地會向外輻射多余熱量，因此，通過目標紅外信息獲得目標的工作狀態是空間態勢感知系統的重要組成部分[2]。

然而，空間目標紅外探測試驗成本較高，通過建立理論模型，開展仿真分析已成為研究目標紅外特性的重要手段。其中，目標軌道外熱流計算是重要一環，主要包括太陽輻射、地球紅外輻射和地球反照輻射。太陽輻射及遮擋的影響相對容易計算，而地球紅外輻射和反照輻射則需要關聯目標姿態和受照射情況，計算過程十分復雜。國內外的相關學者開展了大量研究，主要分為針對簡單結構的解析法或積分法，以及處理復雜結構的蒙特卡洛（Monte Carlo）法。前者在處理常見六面體結構時比較方便，如易樺[3]等提出了一種針對圓軌道航天器外熱流的計算方法，對偏航姿態瞬態外熱流進行了分析。GARZóN[4]等人針對3U 衛星建立了簡化的熱流和溫度模型，分析了貝塔角和導熱率的影響。李志松[5]等人針對微納衛星開展了軌道外熱流計算及在軌溫度分析。李世俊[6]、吳愉華[7]等人將相機簡化為六面體結構，對太陽同步軌道和地球同步軌道外熱流分別進行了分析，得到變姿態條件下的熱流值。

積分方法雖然能夠快速獲得初步結果，但對于處理復雜結構及遮擋問題時則比較困難，且精度不夠。蒙特卡洛方法是商業仿真軟件常用的計算方法，在處理復雜結構航天器及多次反射等情況具有明顯優勢。Atra 等人[8]對衛星熱環境的建模和分析進行了綜述，并對比分析了幾款商業軟件的計算效果。韓玉閣等人[9]利用蒙特卡洛方法計算衛星紅外輻射特征，結果表明衛星的散熱面是判斷衛星是否失效的依據。潘晴[10]等使用反向蒙特卡洛方法對帶天線結構的立方體衛星進行了軌道外熱流的計算，并與商業軟件進行了對比分析，以驗證其精度。劉巨[11]使用STK 軟件獲得太陽矢量關系，使用IDEAS/TMG 模塊獲得了軌道周期內的熱流變化。到目前為止，國內采用自主編程對空間目標外熱流的精細化計算研究的較少。

本文采用蒙特卡洛方法，基于非結構化四面體網格，編寫空間目標軌道外熱流仿真軟件，采用OpenMP 并行加速光線計算，獲得目標在任意時間、任意軌道、任意姿態下的軌道外熱流，并進行了對比驗證。進一步地，對存在遮擋情況下的目標表面熱流情況進行了仿真分析，并結合表面材料屬性，對表面溫度和輻射特性開展了研究。

2 計算方法

本文衛星目標軌道外熱流的計算軟件主要采用矢量坐標變換法，計算順序如下：首先讀取衛星軌道參數，根據衛星軌道信息，計算得到太陽和衛星在J2000 坐標系下的位置坐標，獲得轉換矩陣。然后，讀取網格文件，遍歷每一個表面網格，計算每個表面受到的太陽輻射、地球輻射和地球反照輻射。判斷衛星所在位置是否處于地影區，如果在地影區則將太陽輻射和地球反照太陽輻射值設置為0，不在地影區則總輻射為3 種輻射產生的熱流之和。最后判斷表面是否被遮擋，如果被遮擋，則熱流值設置為0。

軟件運行時衛星軌道信息輸入為6 個（半長軸、偏心率、傾角、升交點赤經、近地點幅角、真近點角），并在內部完成遞推。在不進行姿態調整時，本體系與軌道系重合。軌道系定義為，+X指向飛行方向，+Z指向地心，+Y遵循右手定則。

衛星軌道外熱流計算中，暫不考慮目標自身溫度產生的輻射。下面將逐一介紹太陽直接輻射、地球紅外輻射和地球反照輻射的計算方法。

2.1 太陽輻射計算

太陽輻射及其遮擋情況計算與地球輻射和反照輻射計算不同，不需要對大量光線進行隨機計算。太陽在J2000 系的位置計算過程參見文獻[7]。在獲得太陽矢量Ssun后，目標表面接收到的太陽輻射的熱流密度可表示為

式中，Sc為太陽輻射照度，φ為表面網格法向與太陽光向量的夾角，ρ為面元反射率。太陽輻射照度隨季節變化，通常取平均太陽輻射照度=1 367W/m2進行計算：

式中Rs為日地距離，為平均日地距離。目標在光照區受到太陽輻射，在地影區則設置Sc為0。地影區通過計算當前時刻衛星視角下的地球張角βe判斷，當衛星與地球連線和衛星與太陽連線的夾角大于 βe時，衛星受到光照。βe的計算方法如下：

式中Re為地球半徑，Rsat為衛星與地心的距離。

2.2 地球紅外輻射計算

地球紅外輻射及反照外熱流示意圖如圖1 所示。計算地球紅外輻射時認為地球處于熱平衡狀態，任意表面的輻射強度均勻且相等，則地球輻射熱流的表達式[3,12]為：

圖1 地球紅外輻射及反照外熱流示意圖Fig.1 Schematic diagram of earth’s infrared radiation and albedo external heat flow

其中，ρe為地球反射率，本文設定為0.3，ε為面元紅外發射率；α1和 α2分別為地球表面微元 ds與目標表面微元 dA的連線與二者法向之間的夾角，L為目標表面微元與地球表面之間的距離。令，為地球紅外輻射角系數，是主要求解對象。

通過積分方法計算時，主要判斷面元法向與目標-地球連線的夾角 β 與 βe的余角之間的關系，令k=Re/Rsat，則 ?e可由下式得到：

當(π-arccosk)≤β≤ 時，?e=0。

地球紅外角系數公式的推導詳見文獻[13]。在使用蒙特卡洛方法時，則通過統計光線實現地球紅外角系數的計算，此時地球輻射熱流表達式變換為：

式中，N是面元發出的總光線數量，Nabs是被吸收的光線數量，θi是光線發射位置的天頂角。

2.3 地球反照輻射計算

地球反照輻射計算中將地球設定為漫反射，依據蘭貝特余弦定理，則目標表面接收到的地球反照輻射熱流[3,12]可表示為：

令 ?er為地球反照角系數，(cosα2,0)·max(cosη,0)/πL2ds。由于涉及光線與地球表面交點，以及和太陽的位置關系，計算較為復雜。在工程上，地球反照輻射熱流可以由地球輻射熱流得到:

其中，?為衛星-地球連線與太陽光的夾角。在使用蒙特卡洛方法計算時，地球反照輻射熱流變換為：

2.4 溫度場計算

詳細計算空間目標溫度場的分布十分復雜，本文對計算條件進行了一定簡化。目標在大氣層外時，忽略熱對流；除帆板外，忽略星體各表面之間的導熱和輻射；除太陽和地球外，忽略其他天體輻射的影響。在這種條件下，空間目標接收到的外熱流即為太陽輻射、地球輻射和地球反照輻射。星體表面熱平衡方程如下：

對于太陽能帆板電池片表面，平衡方程如下：

對于太陽能帆板背板表面，平衡方程如下：

其中，αs、αe、αa分別為目標表面對太陽輻射、地球輻射、地球反照的吸收率，太陽輻射和地球反照輻射吸收率取值為材料吸收率，地球輻射吸收系數為表面發射率 ε；qs,i、qa,i、qe,i分別為第i個面元接收到的太陽輻射、地球反照、地球輻射的熱流；Ai是第i個面元的面積，Ns是該表面上面元的總數目，A是該表面所有面元的總面積，σ為斯忒芬－玻耳茲曼常數，σ=5.67×10-8W/(m2·K4)；T為目標表面溫度；T1為帆板電池的表面溫度；T2為帆板背板的表面溫度；ηs為太陽能帆板的光電轉換率，設為0.2；qin為目標表面的內熱源，此處等效為面熱源，對于本體表面，其值為20 W/m2。太陽能帆板蜂窩材料厚度h設為20 mm，熱導率ζ設為1.7 W/(m·K)[14]，表面材料屬性如表1 所示。

表1 表面材料熱參數[15]Tab.1 Thermal parameters of surface material

2.5 計算結果驗證

為了驗證仿真軟件計算結果的準確性，針對800 km 太陽同步軌道和地球同步軌道開展軌道外熱流計算，并與文獻[5]中給出的結果進行對比分析。仿真中采用的光線數為10 000，軌道遞推間隔為90 s，仿真日期選擇為春分日，對比結果如圖2～圖3（彩圖見期刊電子版）所示。其中大寫字母（線條）代表文獻中給出的各表面結果，小寫字母（符號）代表本文仿真結果。其中q為包括太陽輻射、地球輻射和地球反照輻射的軌道總熱流。由圖2～圖3 可以看出，兩種軌道高度下軌道外熱流的計算值與文獻中的計算值誤差小于5%，具有較高的一致性。

圖2 800 km 太陽同步軌道外熱流對比Fig.2 Comparison of external heat flow of the 800 km sun-synchronous orbit

圖3 地球同步軌道外熱流對比Fig.3 Comparison of external heat flow of the geosynchronous orbit

此外，對地球輻射角系數也進行了驗證，如表2 所示，給出了與文獻[16]的對比結果?？梢钥闯?，兩者一致性較好，可以認為本軟件計算精度較高。

表2 典型位置地球紅外角系數隨平板俯仰角變化對比Tab.2 Comparison of earth infrared angular coefficients varying with plate pitch angle at typical positions

2.6 計算域

進一步地，建立帶帆板的模擬衛星模型網格，計算域如圖4 所示。衛星本體系與軌道系重合（本體系原點為衛星質心）。帆板為固定式帆板，分布在+Y和-Y側，衛星本體尺寸為30 cm×20 cm×30 cm，帆板尺寸為30 cm×40 cm×2 cm，帆板邊緣與星體±Y側平面4 cm，網格數量為28 萬。

圖4 計算域Fig.4 Computational domain

2.7 遮擋計算

在計算目標本體結構對太陽輻射的遮擋時，由網格中心點沿著衛星指向太陽的矢量方向發射一條光線。通過給定足夠大的步長，保證矢量終點在計算域外。然后遍歷目標所有面元，如果光線矢量穿插了其他本體表面，則將該光線所在面元太陽輻射熱流值設置為0。

在計算對地球紅外輻射和地球反照輻射的遮擋時，由表面面元向 2π 空間內隨機發射N條光線，跟蹤每條光線的運動。在計算光線運動路徑時，采用臨近網格檢索方法，通過計算光線路徑與四面體網格的4 個三角面元相交情況，當光線穿插本體其他表面時刪除光線，以此循環，直到光線到達計算域邊界。此時，根據光線、目標-地球連線、地球-太陽連線等的關系，統計地球紅外輻射和地球反照輻射量。

3 分析與討論

3.1 不同模式下軌道外熱流計算

空間目標的輻射特性除與所處位置、結構形狀、表面屬性因素相關外，還受姿態變化的影響。對于天基目標觀測尤為明顯。通常，三軸穩定衛星在軌長期模式可分為三軸對地、三軸對日等。本節對固定帆板的衛星在兩種模式下的春分日軌道外熱流進行了計算。對日坐標系定義為，-Z軸指向太陽，+Y軸指向黃北極，X軸遵循右手定則。軌道參數為535 km 太陽同步軌道，降交點地方時為10:30，計算起始時刻為春分日UTC 時間12:30:00，迭代時間間隔為120 s，軌道周期為5 720 s，其中2 838 s 至4 901 s 為地影區。

圖5～圖10（彩圖見期刊電子版）給出了衛星本體各表面在一個軌道周期內的3 種熱流變化情況，其中mode 1 代表三軸對地模式，圖中以實線表示，mode 2 代表三軸對日模式，圖中以點劃線表示?？梢钥闯?，在對地模式下，除+Y面以外，本體其余各表面均受到太陽輻射，由于太陽輻射在總輻射中占比較大，因此，在衛星設計上一般選擇+Y面為散熱面。在各面中，-Y面受到太陽帆板遮擋，太陽輻射在0～550 W/m2范圍內波動，相較于其他受曬表面幅值較小。+Z面長期對地，只有進出地影區時短時間受曬，峰值在510 W/m2左右。除-Z面以外，其他表面均受到地球輻射及地球反照輻射影響。對地模式下，各表面受到的地球輻射恒定，反照輻射隨時間變化。

圖5 +X 面外熱流曲線Fig.5 External heat flow of +X surface

圖6 -X 面外熱流曲線Fig.6 External heat flow of -X surface

圖7 +Y 面外熱流曲線Fig.7 External heat flow of +Y surface

圖8 -Y 面外熱流曲線Fig.8 External heat flow of -Y surface

圖10 -Z 面外熱流曲線Fig.10 External heat flow of -Z surface

在對日模式下，各表面熱流情況與對地模式有顯著不同。太陽輻射僅存在于-Z面，在光照區恒定，約為1 378 W/m2，其他表面只受地球輻射和反照輻射，與太陽輻射相比量級較小，因此受熱相對穩定，在衛星熱設計上，±X、±Y、+Z面均可設計為散熱面。對于地球輻射和反照輻射來說，對日模式和對地模式相比，主要不同是由于姿態變化產生的在時間和幅值上的差別。對于載荷在+Z面的衛星來說，從圖9 可以看出，衛星在兩種模式下受到的3 種輻射值均變化不大，熱環境比較溫和。

3.2 遮擋影響分析

圖11（彩圖見期刊電子版）給出了計算起始時刻對地模式下星體各表面受到的總外熱流云圖。從圖中可以看出，由于衛星降交點地方時選擇為10:30 am，太陽從-Y側照射星體，星體和帆板的-Z、+X，-Y側受光照，+X面熱流值在1 000 W/m2以上，其他面在600 W/m2左右。受帆板遮擋影響，-Y側部分位置不被太陽直接照射，+Y側帆板靠近星體邊緣被遮擋，說明本文編寫的軟件可以很好地處理復雜結構體遮擋情況。

圖11 外熱流分布及遮擋影響示意圖Fig.11 Schematic diagram of external heat flow distribution and occluding effects

表3 給出了一個軌道周期內各表面的平均熱流情況。對于三軸對地模式（mode 1），除+Y面平均熱流在90 W/m2左右外，各表面熱流分布相對均勻，在340～400 W/m2左右變化。而對于對日模式，-Z側的星體表面和帆板表面熱流值在972 W/m2左右，其他表面則小于150 W/m2。針對結構遮擋的影響開展了仿真分析，無遮擋情況下不考慮光線被星體結構或帆板遮擋。計算結果表明，遮擋的影響主要表現在對地模式下，由于帆板的遮擋，-Y表面遮擋后熱流值降低了53.79 W/m2，由于星本體的遮擋，+Y-Z側帆板表面遮擋后熱流值降低了32.05 W/m2。對本體其他表面影響微弱，對于對日模式則無影響。

表3 各平面一個軌道周期內的平均外熱流Tab.3 Average external heat fluxes of each surface in one orbital period W/m2

3.3 表面溫度

結合各表面熱流變化及表1 中給出的表面材料熱參數，根據式（9）～式（11）計算得到各表面溫度。圖12～圖13（彩圖見期刊電子版）分別給出了對地和對日模式下的各表面溫度。由于各表面的質量和比熱容未知，本節計算結果僅考慮穩態解，相較于實際情況波動較大?？梢钥闯?，對地模式下各表面溫度在不同時刻變化較大，其中，星本體-Z面溫度變化最大，約為180 K，+Y側溫度變化最小，約為22 K。帆板由于存在導熱，電池片側（-Z）表面溫度在光照區和地影區的差值相較于本體-Z面較小，約為120 K；隨著時間變化，+X和-X側表面受到陽光照射，溫度變化較大，幅值約在140 K 波動。-Y面在光照區和地影區溫度都比較恒定，分別為283 K 和200 K，+Z表面溫度由于在進出地影區時受到太陽照射，存在兩個峰值，其余時間波動較小。

圖12 對地模式下的各表面溫度Fig.12 Temperature of each surface in the earth-pointing mode

圖13 對日模式下的各表面溫度Fig.13 Temperature of each surface in the sun-pointing mode

在對日模式下，除帆板和本體-Z面在光照區和地影區變化較大外，其余表面變化幅度較小，其中波動最大的是本體+Z面，波動區間約為132 K。對于帆板和本體-Z面，光照區溫度在340 K 左右波動，地影區在220 K 左右波動，整個區間相對恒定。

3.4 表面溫度對比驗證

為了驗證軌道外熱流計算的準確性，對“吉林一號”某衛星帆板溫度進行仿真與在軌實測值驗證分析。衛星運行在535 km 太陽同步軌道，降交點地方時為12:00，帆板位于-Y側，帆板運行模式為對日模式。在實際應用中帆板處于非平衡狀態，因此表面溫度使用以下簡化公式進行求解：

式中，c為比熱容，取值為350 J/(kg·K)，m是帆板質量，取值為1.5 kg，As為帆板面積，此處取值為0.3，Q為帆板表面總的吸收熱流，表面吸收率設為0.9。采用蒙特卡洛方法獲得表面外熱流后，代入公式迭代求解帆板溫度，并與衛星遙測數據做對比，結果如圖14 所示。

圖14 太陽帆板溫度對比Fig.14 Temperature comparison of solar panel

可以看出，計算結果在各個軌道周期內的變化幅值和規律與遙測數據符合較好。最高溫度約為95 °C，最低值約為-76 °C，在光照區升高，在地影區下降，不斷循環。本文計算結果可以為太陽帆板溫度場分析提供可信參考。同時也可以看出，計算值和遙測值在部分區域存在一定偏差，可能是參數設置與實際情況的差異影響了計算結果，需要開展進一步研究。

3.5 紅外輻射強度分布

目標自身紅外輻射由表面溫度和表面材料的發射率決定，在獲得溫度變化曲線后，目標自身紅外光譜輻射出射度可由式（13）表示：

其中，c1為第一輻射常量，值為3.742×10-16W ·m2，c2為第二輻射常量，其值為1.438 8×10-2m·K，λ為波長，單位為μm，Ti是面元表面溫度。εi(λ)是面元光譜發射率。

結合圖12 和圖13 給出的兩種模式下的各表面溫度，認為各表面為漫反射，根據朗伯余弦定律計算得到兩種模式下一個軌道周期內的各方向上的輻射強度，如圖15～圖18（彩圖見期刊電子版）所示。

圖15 對地模式下的各方向輻射強度（3～5 μm）Fig.15 Radiation intensity in each direction in the earthpointing mode (3～5 μm)

圖16 對日模式下的各方向輻射強度（3～5 μm）Fig.16 Radiation intensity in each direction in the sunpointing mode (3～5 μm)

圖17 對地模式下的各方向輻射強度（8～14 μm）Fig.17 Radiation intensity in each direction in the earthpointing mode (8～14 μm)

圖18 對日模式下的各方向輻射強度（8～14 μm）Fig.18 Radiation intensity in each direction in the sunpointing mode (8～14 μm)

可以看出，在常用的3～5 μm 和8～14 μm 兩個探測波段中，兩種模式下+Z和-Z方向的輻射強度較高，3～5 μm 波段最大值在1.5 W·sr-1左右。8～14 μm波段最大值在22 W·sr-1左右。主要原因是帆板溫度較高且面積較大，是紅外信號的主要來源。8～14 μm 波段輻射強度明顯強于3～5 μm 波段。同一譜段內，對地模式下，不同時間內，+X、-X和-Y方向也存在紅外信號較強時刻。但對日模式下，除+Z和-Z方向的輻射強度較高外，其他方向紅外輻射值較小，目標探測存在困難。

4 結論

本文針對繞地衛星軌道外熱流，采用蒙特卡洛法、基于非結構四面體網格和OpenMP 并行算法編寫了仿真軟件。針對文獻中的典型工況進行了對比分析，驗證了軟件計算結果的準確性。進一步，針對535km太陽同步軌道，考慮衛星結構遮擋的影響，對不同姿態控制策略下的衛星軌道外熱流進行了仿真分析，并將太陽帆板表面溫度與在軌遙測數據進行了對比分析。得出如下結論：

（1）不同姿態模式下目標的軌道外熱流區別較大。對于535km太陽同步，10:30地方時軌道來說，對地模式下除+Y面外，各表面熱流值隨時間變化波動較大，而對日模式下除-Z面本體及帆板表面波動較大外，其他表面變化較小。

（2）蒙特卡洛算法對空間目標復雜結構及結構遮擋具有很好的適應性，對地模式下，考慮遮擋后，-Y表面熱流值降低了53.79 W/m2，+Y-Z側帆板表面熱流值降低了32.05W/m2。

（3）不同模式下目標各表面溫度特性不同。對地模式下各表面溫度隨時間波動較大，使紅外觀測窗口規劃提高了難度。在兩種模式下，帆板在光照區溫度較高，具有明顯的紅外特征，便于開展空間目標紅外觀測。