張昊春，曲博巖，金 亮

(哈爾濱工業(yè)大學(xué) 能源科學(xué)與工程學(xué)院，黑龍江 哈爾濱 150001)

引言

隨著光電子技術(shù)的發(fā)展，紅外成像、光電制導(dǎo)及跟蹤、激光測距、激光雷達、聲吶等技術(shù)開始廣泛應(yīng)用于陸?？斩喾N作戰(zhàn)平臺，從而大大改變了現(xiàn)代戰(zhàn)場的攻防形態(tài)。在很長一段時間里，針對戰(zhàn)斗機等飛行目標的探測手段主要分為雷達探測和非雷達探測兩種。而鑒于目前所開發(fā)的新一代飛行器多采用隱形材料，使得雷達的探測有效性明顯降低，而紅外探測技術(shù)以其不受無線電干擾，機動性較好，探測精確性高等優(yōu)勢，顯著地彌補了傳統(tǒng)雷達對飛行器目標探測能力的不足[1]，進而得到了發(fā)展。

飛行器在高速飛行時機身與周圍空氣摩擦將產(chǎn)生大量的熱，加上機體表面接收到的來自太陽、地表和大氣背景輻射的影響，因此機身對外會發(fā)出強烈的紅外輻射。哈爾濱工業(yè)大學(xué)的夏新林等[2]通過研究飛機蒙皮溫度場中的傳熱過程，建立起飛機內(nèi)部熱傳導(dǎo)、飛機蒙皮表面氣動對流及輻射的瞬態(tài)傳熱模型。北京航空航天大學(xué)的呂建偉等[3]通過反向的蒙特卡洛法(RMC)，引入表面多重遮擋的算法，計算具有復(fù)雜幾何外形的F22機身紅外輻射特性。國防科技大學(xué)的盧建等[4]通過3DMAX軟件建立了飛行器目標的三維坐標模型以及傳感器光學(xué)成像系統(tǒng)仿真模型。王科偉等[5]介紹了紅外成像波門形心跟蹤算法的原理，通過建立跟蹤誤差模型，討論分析了紅外成像波門形心跟蹤算法的誤差。成剛等[6]對多光譜成像探測技術(shù)進行了詳細研究，針對迷彩偽裝目標偵察設(shè)計了無人機載多光譜攝像機，重點分析了多光譜偵察探測、識別的計算方法。韓軍等[7]提出一種提高現(xiàn)有光電觀瞄系統(tǒng)成像空間分辨率的新方法。在不改變陣列探測器件像元尺寸和不移動探測器的前提下，利用雙光楔的較大移動使像平面發(fā)生微小位移，實現(xiàn)對探測器各相鄰像元和不感光間隔目標進行微位移采樣，提高目標信息的采樣率，通過超分辨重建技術(shù)來提高系統(tǒng)的分辨能力，達到提高光電觀瞄系統(tǒng)空間分辨率的目的。Nicolas Douchin[8]論述了如何利用SE-Workbench-EO軟件模擬飛機紅外輻射特征以及獲得實時渲染畫面，并采用寬帶、窄帶模型及Curtis-Godson逼近法計算出飛機羽流區(qū)內(nèi)的輻射亮度。Haynes[9]介紹了一種紫外至紅外波段的場景生成軟件CAMEO-SIM的各組件及其相應(yīng)的功能，說明了該軟件中加入紅外痕跡及陰影、尾焰、傳感器效應(yīng)等因素的方法。Gretzmacher Floris M等[10]采用一種改進的歐幾里得距離測量的視覺圖像來進行偽裝評估，在圖像中形狀相似的選擇區(qū)域可以可視化。Bastiaansen M C等[11]采用丹麥國防研究機構(gòu)(DDRE)開發(fā)的計算機偽裝評估方法CAMEVA，輸入由高分辨率目標和適當數(shù)量的背景組成的單個圖像。Jacobs P.A.M等[12]為了確定偽裝材料跟蹤背景溫度變化的能力，在各種氣象條件下進行偽裝和背景溫度的測量，以確定偽裝材料有效降低目標特征的情況。

截至目前，在紅外成像領(lǐng)域，較多學(xué)者針對飛行器自身的紅外輻射特性進行了成像仿真研究，而較少考慮了太陽輻射與天空背景對成像的影響，導(dǎo)致成像結(jié)果有一定誤差。而當探測時刻、地點或季節(jié)改變時，太陽相對飛機的位置變化將對飛機紅外成像產(chǎn)生較大的影響，綜合考慮太陽輻射等環(huán)境因素，針對不同的飛行姿態(tài)及速度，不同的探測時刻、地點對飛機紅外成像的影響進行研究。

1 飛機紅外輻射特性分析

1.1 飛機模型

采用COMSOL軟件建立J-16飛機簡化后的三維幾何模型，并對機體外表面劃分非結(jié)構(gòu)網(wǎng)格。機體表面網(wǎng)格沿飛機中軸線呈左右對稱分布，并將機頭、發(fā)動機兩處溫度變化梯度較大的位置上的網(wǎng)格進行加密處理，幾何模型與劃分結(jié)果分別如圖1和圖2所示。

圖1 飛機幾何模型示意圖Fig.1 Diagram of aircraft geometric model

圖2 網(wǎng)格劃分結(jié)果示意圖Fig.2 Diagram of grating division result

1.2 紅外輻射特性

氣動加熱的形成條件為：飛機處于高速運動過程中時，機體周圍的氣流將有一部分動能不可逆地轉(zhuǎn)化為熱能，從而在蒙皮表面產(chǎn)生一層熱層，與機身蒙皮進行換熱[4]。目前對于飛機蒙皮溫度分布的計算主要有數(shù)值仿真與經(jīng)驗公式兩種方法。第一種是通過ANSYS等軟件進行仿真模擬，進而得出蒙皮的溫度場分布。在工程應(yīng)用中，為了簡化計算，常常采用經(jīng)驗公式((1)式)計算得到駐點溫度Tr，以之代替蒙皮溫度。

(1)

式中：Ta為處于飛機飛行高度上的大氣溫度，單位K，可以查表獲得；k即氣體絕熱指數(shù)，文中取為1.4；γ即恢復(fù)系數(shù)，其中層流一般取0.85，湍流取為0.83；Ma即飛機飛行速度馬赫數(shù)。由普朗克定律可知，蒙皮區(qū)域單個面元上的光譜輻射出射度等于：

(2)

式中：λ為波長，單位μm；T為絕對溫度，單位K；c1、c2分別為第一、第二輻射常數(shù)。另外，為計算飛機蒙皮的定向紅外輻射亮度，在t時刻將沿視線方向能夠觀察到的若干面元的輻射強度作積分處理，即得蒙皮區(qū)域的總輻射強度：

(3)

式中：ε即面元的發(fā)射率，文中均取為0.9；Ai為面元i的面積；ωi為面元i的外法線相對于視線的交角。因為目標與其所處背景的紅外輻射強度分布不同，隨著探測器成像單元接收到的輻射強度增大，對應(yīng)的像素點灰度值增加，根據(jù)這一點，在捕捉到的紅外圖像中應(yīng)用灰度等級表示目標區(qū)域的輻射能量分布?；叶鹊燃売成潢P(guān)系由飛機蒙皮區(qū)域紅外輻射亮度的計算結(jié)果建立，計算公式如下[15]：

GL=[r+(1-r)·(L-Lmin)/

(Lmax-Lmin)]×255

(4)

式中：GL表示某一像素點的灰度值，取值范圍為0～255；r表示環(huán)境泛光的紅外等效值，取值范圍為0～1；Lmin和Lmax表示成像域輻射亮度的下限與上限，L即為該像素點處紅外輻射亮度的初始值。最后根據(jù)像素點輻射溫度和灰度值的轉(zhuǎn)化關(guān)系得出飛機蒙皮紅外仿真圖像，由紅外測溫原理知，飛機蒙皮任一面元輻射溫度Ts為[16]

(5)

式中：τa為探測器到達飛機蒙皮位置處的大氣透過率；ε(χ)為λ1～λ2波段內(nèi)飛機蒙皮面元法線方向的平均發(fā)射率；χ為入射角，即探測方向與面元的夾角；Tr為蒙皮溫度，根據(jù)(1)式進行計算；n的取值與波段相關(guān)，在2 μm～5 μm波段，n=8.68，在8 μm～13 μm波段，n=4.09。

2 太陽輻射模型

2.1 太陽高度角及方位角

太陽高度角Hs及方位角As的計算公式如下：

sinHs=sinφsinδ+cosφcosδcost

(6)

cosAs=(sinHs·sinφ-sinδ)/(cosHs·cosφ)

(7)

式中：φ為飛機所處的緯度；δ為太陽赤緯；t為太陽時角。太陽赤緯的計算公式如下：

δ=0.006 918-0.399 912cos(b)+

0.070 257sin(b)-0.007 758cos(2b)+

0.000 907sin(2b)-0.002 697cos(3b)+0.001 48sin(3b)

(8)

式中：b=2π(N-1)/365，N從每年的1月1號算起，直到仿真時所取日期為止的天數(shù)(1月1日則N=1，1月2日N=2，以此類推)。時角的計算公式如下：

t=ts·15°

(9)

式中：ts為視太陽時，其值為平太陽時(即地方時)減去12。太陽時角在正午時為 0 (即位于天空的正中間)，上午為負，下午為正，日出時為-90°，日落時為+90°，平均每小時時角變換為15°。以中國東北某城市(E125°，N44°，海拔170 m)為例，零方位角為S方向，計算6:00～18:00期間內(nèi)太陽高度角與方位角的變化曲線，結(jié)果如圖3和圖4。

圖3 某市太陽高度角Fig.3 Solar height angle in a certain city

圖4 某市太陽方位角Fig.4 Solar azimuth angle in a certain city

2.2 太陽輻射面積比例因子

由于所取飛機模型的結(jié)構(gòu)比較復(fù)雜，各個面元之間彼此存在著遮擋的情況，故需定義一個與面元遮擋關(guān)系相關(guān)的參量，即太陽輻射的面積比例因子。假設(shè)面元a由于b的遮擋，只有總面積中的一部分能夠接收到太陽的輻射能量，于是面積比例因子D定義為

(10)

式中：Ai即某面元的總面積大??；Ai,sun為面元能夠接受到太陽輻射的面積大小。遮擋可分為完全遮擋和不完全遮擋兩種情況，如果此面元完全沒有被其他面元遮擋，則D=1；而當此面元完全被其他面元遮擋時，D=0，由此可知D的取值范圍在0到1之間。

基于反向蒙特卡洛法確定每個面元的面積比例因子。首先，在一個面元中取總數(shù)為N的隨機分布的點，然后根據(jù)這個面元所接收的太陽光的入射方向，令所取的每個點沿著其反方向射出一束光線，再通過空間幾何關(guān)系推斷出每束光線與其余的面元是否相交。只要有一束光線與其余面元相交，即說明該隨機點被遮擋(所有的光束都相交即為完全遮擋)；如每一束都不相交，說明該點可以接受到全部的太陽輻射。對所有光線的判斷循環(huán)完成后，統(tǒng)計面元被遮擋的點的數(shù)目，從而得到可被照射到的點的數(shù)目，最終求出每個面元的D值。根據(jù)以上方法，假定太陽入射方向向量為(-1,-1,-1)時，D值的計算結(jié)果如圖5所示。

圖5 D值分布情況示意圖Fig.5 Distribution schematic diagram of solar radiation area scale factor

2.3 太陽輻射值

太陽常數(shù)定義為太陽在大氣層外產(chǎn)生的總輻照度Eall=1 353 W/m2，太陽的平均半徑約6.363 8×105km，日地平均距離約為1.5×108km[14]。當計算到達飛機所在位置的太陽輻照度時，需要考慮的影響因素諸如太陽的高度角、觀察者的海平面高度、天空中云霧、塵埃等，情況比較復(fù)雜。在工程計算中，太陽輻射值Es(W/m2)一般使用如下經(jīng)驗公式來進行估算：

Es=[1-A(U*,β)]·(0.349Eall)sinβ+

(11)

表1 不同地表狀態(tài)反射率Table 1 Reflectivity of different surface states

3 紅外成像仿真

3.1 坐標系轉(zhuǎn)換

在整個成像仿真過程中，涉及到各坐標系之間的轉(zhuǎn)換，其中包括用于標定飛機三維坐標的大地坐標系，描述飛機尺寸和飛行姿態(tài)的飛機坐標系，表示探測器三維坐標的探測器坐標系，確定成像位置的成像面坐標系以及位于焦平面的像元坐標系。

為表明各個坐標系之間的相對位置關(guān)系，畫出成像仿真坐標系統(tǒng)示意圖，如圖6所示。其中主要涉及的坐標系包括大地坐標系Xd-Yd-Zd，飛機模型坐標系X-Y-Z，成像面坐標系Xp-Yp和探測系統(tǒng)坐標系Xs-Ys-Zs。

圖6 系統(tǒng)坐標示意圖Fig.6 System coordinate diagram

進行三維空間的坐標旋轉(zhuǎn)變換要求已知旋轉(zhuǎn)角度和旋轉(zhuǎn)軸，假設(shè)在右手系中，全部坐標點繞x軸、y軸、z軸旋轉(zhuǎn)的角度分別為α、β和γ，旋轉(zhuǎn)變換的變換矩陣R如下式[9]：

(12)

假設(shè)坐標系中的全部坐標點相對于x軸、y軸、z軸的平移量分別為tx、ty、tz，平移變換的變換矩陣T如下式：

(13)

由(10)式和(11)式，得出飛機模型坐標系到大地坐標系的轉(zhuǎn)換關(guān)系式：

(14)

探測器坐標系到成像坐標系的轉(zhuǎn)換關(guān)系式如下：

(15)

式中f為焦距。接下來在成像坐標系中確定像元的位置坐標，假設(shè)像元坐標系下某個像元在x和y方向上的寬度分別為Δm和Δn，于是根據(jù)下式可將成像坐標系下的投影點映射到像元坐標系。

(16)

式中[x]代表對x向下取整。飛機機身蒙皮離散為很多個三角面元，如果三角面元至少有2個點位于某成像單元的有效成像范圍內(nèi)，則該面元在這個成像單元內(nèi)成像[15]。將機身蒙皮上各個面元逐一投影到像元坐標系上，即可得到飛機蒙皮面元在該坐標系下的3個頂點坐標，由這3個頂點確定的三角面之和即為目標的成像區(qū)域，而探測器陣列上其余位置則為背景區(qū)域。

3.2 仿真結(jié)果及分析

假定在大地坐標系下，飛機的三維坐標為(0,0,1 000)，探測器坐標為(0,0,0)，飛行速度1馬赫，大氣模式為中緯度夏季模式，探測地點為E120°，N30°，地表狀態(tài)為干燥裸地，探測時間為北京時間2019年6月1日星期六，8:00AM，探測波段取8 μm～12 μm，仿真結(jié)果如圖7所示，其中圖7(a)為灰度等級由(5)式轉(zhuǎn)化而得的紅外輻射溫度(單位K)，圖7(b)為處理后相應(yīng)的RGB圖像。

圖7 飛行速度1 Ma時仿真結(jié)果Fig.7 Simulation results at flight speed of 1 Ma

根據(jù)計算結(jié)果，結(jié)合圖像進行分析知，此時飛機處于高速運行狀態(tài)，氣動加熱效果顯著，機身的發(fā)射輻射中來自太陽輻射的份額占總輻射強度的比例較小，所以蒙皮的輻射亮度分布主要由氣動加熱層定向輻射強度決定；另一方面由于太陽光傾斜照射，飛機機頭被直接照射的區(qū)域亮度明顯高于其他區(qū)域。經(jīng)計算，圖7中框區(qū)域內(nèi)的輻射亮度比其他區(qū)域輻射亮度高18.2%。

飛機在飛行過程中速度的變化將會影響氣動加熱效果，進而影響仿真結(jié)果，另外，飛行高度的變化將會使得探測器光學(xué)路徑長度改變，也會對結(jié)果產(chǎn)生影響?，F(xiàn)保持上述計算條件不變的情況下，改變飛機的運行速度和高度，分別研究二者帶來的影響，其中飛行速度分別取為0.5 Ma、1 Ma(圖7)、3 Ma和6 Ma，將得到的仿真結(jié)果和圖7統(tǒng)一轉(zhuǎn)化為輻射溫度，結(jié)果如圖8至圖10所示。飛行高度分別取1 000 m(圖7)、4 000 m和8 000 m，最終得到的結(jié)果如圖11所示。

圖8 飛行速度0.5 Ma時的仿真結(jié)果Fig.8 Simulation results at flight speed of 0.5 Ma

圖9 飛行速度3 Ma時的仿真結(jié)果Fig.9 Simulation results at flight speed of 3 Ma

圖10 飛行速度6 Ma時的仿真結(jié)果Fig.10 Simulation results at flight speed of 6 Ma

圖11 改變飛行高度時的仿真結(jié)果Fig.11 Simulation results when flying altitude is changed

從圖中可以看出，隨著飛行速度的提高，飛機蒙皮輻射亮度逐漸增大，這一現(xiàn)象在翼面上更加明顯(圖8至圖10)，這是由于當速度提高時，飛機與表面氣流摩擦生成的熱量更高，換熱效果更加明顯，從而使蒙皮區(qū)域溫度和輻射亮度提高。圖8至圖10中框區(qū)域內(nèi)的輻射亮度比其他區(qū)域輻射亮度分別高21.4%，17.9%和16.1%。另一方面，隨著飛行高度的增加，目標在像平面中所占面積逐漸減小，當飛機從1 000 m高度上升到4 000 m時，蒙皮區(qū)域紅外輻射亮度整體上變化很小，分析認為，這是由兩方面因素導(dǎo)致的，首先隨著飛行高度的增加，目標與探測器之間大氣散射路徑增大，大氣介質(zhì)中輻射能量的吸收和散射衰減增強，探測器接收到的輻射能量減小，所以仿真得到的目標紅外輻射亮度減?。煌瑫r，隨著飛行高度增加，對應(yīng)于單個成像單元的飛機表面劃分的面元數(shù)量增加，從而使像平面接收的每個成像單元紅外輻射能量增加。所以，由于投影面元疊加和大氣衰減的雙重作用，從1 000 m到4 000 m高空的紅外仿真圖像整體變化很小。而當飛機從4 000 m爬升至8 000 m時，蒙皮區(qū)域整體的輻射亮度明顯下降，此時大氣衰減成為主導(dǎo)因素。

圖12 飛機偏航，俯仰及旋轉(zhuǎn)示意圖Fig.12 Diagram of aircraft yaw, pitch and rotation

接下來研究探測時間以及飛行姿態(tài)對成像效果的影響。假定飛行速度為0.5馬赫，飛行高度為1 000 m，大氣模式為中緯度夏季模式，探測地點為E120°，N30°，地表狀態(tài)為干燥裸地。以上條件均控制不變，現(xiàn)設(shè)定日期為6月1日，探測波段仍取8 μm～12 μm，探測時間及飛行姿態(tài)控制參數(shù)見表2，得到的仿真結(jié)果如圖13至圖16所示。規(guī)定飛機左偏航時偏航角為正，抬頭時俯仰角為正，右滾轉(zhuǎn)時旋轉(zhuǎn)角為正，如下圖所示。

圖13 12:00AM的仿真結(jié)果Fig.13 Simulation result at 12:00AM

圖14 14:00PM的仿真結(jié)果Fig.14 Simulation result at 14:00PM

圖15 16:00PM的仿真結(jié)果Fig.15 Simulation result at 16:00PM

圖16 18:00PM的仿真結(jié)果Fig.16 Simulation result at 18:00PM

表2 仿真計算參數(shù)設(shè)置Table 2

由計算結(jié)果分析知，因飛行速度較低，氣動加熱效果不明顯，機身的發(fā)射輻射中來自太陽輻射的份額占總輻射強度的比例較大，所以蒙皮的輻射亮度分布主要由接收的太陽輻射條件(高度角，面積比例因子等)決定。當探測時間為中午12點時，機身被太陽照射的面積上輻射亮度呈對稱分布，同時因機頭處D值較小，所以輻射亮度較小，經(jīng)計算得圖8中框區(qū)域的輻射亮度比其他區(qū)域輻射亮度低7.9%。

4 結(jié)論

建立了飛機目標紅外成像仿真模型，在此基礎(chǔ)上，著重考慮不同探測時間下，太陽輻射條件的變化對目標紅外成像特性的影響。仿真結(jié)果表明。

1) 飛機低速度運行時，飛機蒙皮各個區(qū)域溫差相對較小，由太陽輻射所造成的蒙皮直射區(qū)域的溫升較為明顯，太陽輻射亮度占總輻射亮度比例較大，蒙皮的輻射亮度分布主要由氣動加熱層定向輻射強度決定。

2) 飛機高速運行時，由太陽輻射量所引起的蒙皮直射區(qū)域的溫升相對較小，太陽輻射亮度占總輻射亮度比例較小，蒙皮的輻射亮度分布主要由接收的太陽輻射條件決定。

3) 調(diào)整探測時間，將會影響飛機位置所接收到的太陽輻射量，進而影響蒙皮區(qū)域的輻射亮度分布。當太陽光傾斜入射時，飛機機頭被直接照射的區(qū)域比其他區(qū)域的輻射亮度高18.2%；調(diào)整飛行姿態(tài)對蒙皮區(qū)域紅外成像特性局部特征產(chǎn)生影響。