孫建平,齊 宏,王申領(lǐng),趙松慶,郝燕云,吳根水

(1.哈爾濱工業(yè)大學(xué)能源與科學(xué)工程學(xué)院,黑龍江 哈爾濱150001； 2.中國空空導(dǎo)彈研究院,河南 洛陽 471009；3.航空制導(dǎo)武器航空重點實驗室,河南 洛陽 471009)

1 引 言

隨著國防科技的發(fā)展,高速隱身戰(zhàn)機在大國軍事戰(zhàn)略競爭中的地位變得舉足輕重,而相對的制導(dǎo)則成為其中關(guān)鍵技術(shù)。目前的制導(dǎo)技術(shù)主要包括雷達制導(dǎo)、紅外制導(dǎo)和地圖定位制導(dǎo)等方法[1-4]。目前,隨著計算機和紅外傳感器性能的提高及雷達反射面的逐漸減小,紅外探測技術(shù)逐漸成為研究主流。

對于飛行器紅外目標(biāo)探測,國外研究較早且開發(fā)了相應(yīng)計算軟件。例如:美國RIT在美國DSC公司開發(fā)THERM模型基礎(chǔ)上開發(fā)了DIRSIG紅外仿真計算軟件[5-6]；美、英等國聯(lián)合開發(fā)了可用于飛機紅外信號仿真的NIRTRAN平臺[7]。國內(nèi)方面近年來發(fā)展迅速。例如:黃偉等[8-9]各輻射源及各種飛機飛行環(huán)境下的飛機蒙皮的紅外反射特性進行了系統(tǒng)的分析研究,并認為蒙皮發(fā)射率對紅外信號由很大影響；王霄等[10]和韓平麗等[11]分別對不同背景下的飛機進行了紅外目標(biāo)成像研究。

飛行器的紅外目標(biāo)探測以其表面的紅外反射特性分布規(guī)律為基礎(chǔ),而飛行狀態(tài)與飛行背景對飛行器表面反射特性具有重大影響,因此,本文采用有限元與商業(yè)軟件結(jié)合的方法對不同飛行狀態(tài)復(fù)雜背景下飛行器的表面紅外反射特性的分布規(guī)律進行研究,以為飛行器的紅外目標(biāo)探測提供參考。

2 理論計算模型

2.1 Schlick模型

1982年,Schlick提出一種計算材料表面雙向反射分布函數(shù)(BRDF)的數(shù)學(xué)模型,以表征不同材料的表面反射特性,具體可由式(1)表達[12]:

f(t,u,v,l,w)=Sλ(u)·D(t,v,l,w)

(1)

式中,Sλ表示入射方向函數(shù),具體如式(2):

Sλ(u)=Cλ+(1-Cλ)(1-u)5

(2)

其中,Cλ表示波長反射因子。

式(1)中,D表示BRDF模型的方向因子,由式(3)表示:

(3)

式中,Z和A分別為關(guān)于變量r和變量p的極坐標(biāo)函數(shù);G(v)和G(l)分別為反射和入射的輻射遮擋因子,各函數(shù)具體表達式如式(4)所示:

(4)

式(1)～(3)中,變量t、u、v、l和w為Schlick模型中相應(yīng)角度的余弦值,各變量具體表達式及相應(yīng)角度的定義見文獻[12],在此不再贅述。

2.2 坐標(biāo)系轉(zhuǎn)換及輻射遮擋判斷因子計算

隨著飛機的運動,飛機與探測器的相對位置一直在變化,同時對于飛機上單個面元來說,入射輻射方向與反射輻射方向也是不斷變化的。因此,需要建立當(dāng)?shù)刈鴺?biāo)系-飛機坐標(biāo)系和飛機坐標(biāo)系-面元坐標(biāo)系之間的坐標(biāo)轉(zhuǎn)換關(guān)系,以此得到探測器接收到的輻射信號與飛機面元之間直接的對應(yīng)關(guān)系。

飛機表面形狀的復(fù)雜性與運動軌跡和姿態(tài)的多變性導(dǎo)致部分飛機表面入射輻射為0,且探測器接收到的輻射信號也會存在遮擋現(xiàn)象,因此本文引入了兩個輻射遮擋判斷因子來表征這一現(xiàn)象。

2.2.1 當(dāng)?shù)?地球坐標(biāo)系

飛機所在地的太陽輻射初始入射方向取決于當(dāng)?shù)貢r間、當(dāng)?shù)鼐暥群图竟?jié)等參數(shù)。飛機緯度信息則取決于飛機與地球的相對位置,因此,需要建立當(dāng)?shù)?地球坐標(biāo)系描述這一關(guān)系。

建立以地球為參照的當(dāng)?shù)刂苯亲鴺?biāo)系i-j-k,規(guī)定i軸方向為正面向西,k軸方向為地面法向,根據(jù)右手法則確定j軸方向,坐標(biāo)原點為飛機所在位置P,如圖1所示。

圖1 當(dāng)?shù)?地球直角坐標(biāo)系

2.2.2 飛機坐標(biāo)系

飛機在飛行過程中存在旋轉(zhuǎn)和仰角變化,且其飛行高度和飛行速度同樣變化劇烈,其運動姿態(tài)和運動狀態(tài)的多變性導(dǎo)致飛機表面獲得的投射輻射方向不斷變化。同時,由于飛機與探測器相對位置不斷改變,探測器的探測方向也不斷改變。因此,需要建立飛機坐標(biāo)系,并明確當(dāng)?shù)刈鴺?biāo)系-飛機坐標(biāo)系之間的坐標(biāo)轉(zhuǎn)換關(guān)系。

建立飛機坐標(biāo)系x-y-z,定義x軸方向為中軸線向后,y軸方向為水平面向右,右手定則確定z軸方向,如圖2所示。

圖2 飛機直角坐標(biāo)系

則當(dāng)?shù)刈鴺?biāo)系中坐標(biāo)(i,j,k)與飛機坐標(biāo)系中坐標(biāo)(x,y,z)之間的轉(zhuǎn)換關(guān)系如式(5)所示:

(5)

當(dāng)飛機旋轉(zhuǎn)時,其繞x、y和z軸的右旋角度分別為θx、θy和θz,于是式(5)中的轉(zhuǎn)換矩陣可由式(6)表示:

(6c)

2.2.3 面元坐標(biāo)系

BRDF的計算是在每一個飛機表面的面元上進行的,即需要得到每一個面元的入射輻射和反射輻射的方向參數(shù),因此需要將飛機表面的入射輻射方向轉(zhuǎn)化為每個面元的入射輻射方向參數(shù),因而需要建立面元坐標(biāo)系,并明確當(dāng)?shù)?飛機-面元之間坐標(biāo)位置關(guān)系。

圖3 面元直角坐標(biāo)系

面元坐標(biāo)系如圖3所示。圖中P1、P2、P3和P4組成了一個四面體網(wǎng)格,P1、P2和P3為飛機表面面元的三個頂點,P4為飛機內(nèi)部節(jié)點,以此為基礎(chǔ),建立面元坐標(biāo)系u-w-v,規(guī)則如下:

(1)u軸:P1為原點,向量P1P2方向為u軸方向,其方向余弦可由式(7)表示:

(7)

(2)w軸:面元外法向方向為w軸方向,可由面元面方程Ax+By+Cz+D=0與飛機內(nèi)部點P4確定,如式(8):

(8)

式中,A、B和C為面元面方程系數(shù),可根據(jù)P1、P2和P3坐標(biāo)由式(9)確定:

(9)

(3)v軸:根據(jù)u軸與w軸,通過右手法則確定方向,其方向余弦如式(10):

(10)

三個坐標(biāo)系確定后,結(jié)合式(5)可以確定飛機表面上某點在當(dāng)?shù)刈鴺?biāo)系中的坐標(biāo)(i,j,k)與其所在的面元坐標(biāo)系的坐標(biāo)(u,w,v)之間的轉(zhuǎn)換關(guān)系,即式(11):

(11)

2.2.4 輻射遮擋判斷因子計算

為判斷單個面元是否有投射輻射和能否被探測器探測,本文引入投射輻射遮擋因子和探測遮擋因子,并利用反向蒙特卡洛法對這兩遮擋因子進行計算。

在無其他飛機面元遮擋時,記投射輻射直接遮擋因子為Rs1和探測直接遮擋因子為Rd1,記由面元重心出發(fā)沿投射反方向(探測方向反方向)的追蹤光束的方向向量為sis和面元法向量為ni,則Rs1和Rd1可由式(12)計算:

(12a)

(12b)

在考慮其他面元的遮擋時,記投射輻射間接遮擋因子為Rs2和探測間接遮擋因子Rd2。考慮面元j對面元i的遮擋情況,如圖4所示。

圖4 面元遮擋關(guān)系示意圖

記面元i重心為G,沿投射輻射方向向量sis方向反向追蹤光束交面元j于點R,則根據(jù)面元j面方程計算可得G到面元j距離lCR和面元j法向量與輻射方向夾角b,并由此得到GR距離lRG,進而求得R點坐標(biāo)(xR,yR,zR)。最后根據(jù)R點與面元j其余三點組成的三角形面積是否與面元j本身面積相等判斷R是否位于面元j內(nèi)。若R位于面元j內(nèi),則Rs2和Rd2為0,否則為1。

記投射輻射遮擋因子和探測遮擋因子分別為Rs和Rd,綜合直接與間接遮擋因子,則兩遮擋判斷因子可由式(13)計算:

(13)

2.2.5 紅外探測模型

在近距離紅外探測時,對于相同探測方向,飛機上不同面元相對于探測器的探測方向并不相同,因此需要建立探測模型,將以當(dāng)?shù)刈鴺?biāo)系為基準(zhǔn)的探測器探測方向轉(zhuǎn)換為飛機坐標(biāo)系為基礎(chǔ)的探測方向?？紤]100 m探測距離模型,如圖5所示。

圖5 紅外探測模型

探測方向天頂角為θd,圓周角為φd,當(dāng)?shù)刈鴺?biāo)系中探測點位置Pd(id,jd,kd)經(jīng)過坐標(biāo)轉(zhuǎn)換可得飛機坐標(biāo)系探測點位置Pdp(xd,yd,zd)和所在面元重心位置PGp(xGp,yGp,zGp),進而可以得到以飛機坐標(biāo)系為基準(zhǔn)的每個面元的探測方向,如式(14)所示:

(14)

式中,dp代表飛機坐標(biāo)系中探測器位置;Gp代表飛機坐標(biāo)系中面元重心位置。

最后,通過飛機坐標(biāo)系與面元坐標(biāo)系之間的轉(zhuǎn)換,即有計算BRDF所需探測方向的天頂角θr和圓周角φr。

2.3 飛機表面投射輻射模型

2.3.1 太陽輻射

太陽赤緯δ0由飛行時的日期決定,可由式(15)計算:

式中,n為從1月1日起已經(jīng)度過的天數(shù)。

太陽時角t0由飛行時刻決定,可由式(16)計算:

(16)

式中,t為一天中所處的時刻。

考慮飛行緯度對日出日落的影響,對飛機所處時刻位于日出前或日落后進行判斷,如式(17):

(17)

式中,ω為飛機所在緯度。

當(dāng)?shù)刈鴺?biāo)系中太陽入射輻射的天頂角θ和圓周角φ可由式(18)計算:

(18)

面元i獲得的太陽投射輻射力可由式(19)計算:

(19)

式中,Rsun為太陽投射方向的遮擋因子;Esun為太陽輻射到地球表面的輻射照度(取635W/m2);Si為單個面元i的面積;θi為太陽輻射方向與面元i面法向向量夾角;Pir為太陽輻射直射比例(取0.75);BN,Tsun為N個譜帶內(nèi)太陽的黑體輻射力(本文太陽表面溫度取5762 K)。

面元i反射的太陽直射輻射強度由式(20)計算:

(20)

面元i獲得的太陽散射輻射力由式(21)計算:

(21)

式中,Psc為太陽輻射中的散射比例,取0.25。

面元i反射的太陽散射輻射可由式(22)計算:

(22)

2.3.2 海面反射的太陽輻射

Cox與Munk[13-14]在1954年給出了海面反射率的計算模型,該模型較為復(fù)雜,在此只給出結(jié)論性公式(23),各變量含義及模型的具體內(nèi)容見文獻[13]和文獻[14]。

(23)

式中,r為海面反射率,計算式(23)所需的大氣折射率n1取1,海水折射率n2取1.334。

面元i獲得的海面反射的太陽輻射由式(24)計算:

(24)

式中,re代表反射;τλ代表大氣光譜透過率。

面元i對海面反射的太陽輻射的反射輻射強度由式(25)計算:

(25)

2.3.3 海面自身發(fā)射輻射

海面發(fā)生的輻射強度在本文的研究中并不占主導(dǎo)作用,因此可對海面自身發(fā)射的輻射模型進行一定簡化[15]。認為海面自身發(fā)射為漫發(fā)射體,發(fā)射率為0.98,則面元i獲得的海面自身發(fā)射輻射可由式(26)計算:

(26)

式中,sea代表海面自身發(fā)射的輻射。

面元i反射的海面自身發(fā)射輻射的輻射強度由式(27)計算:

(27)

2.3.4 大氣輻射

大氣本身存在一定自身發(fā)射輻射,可用式(28)計算:

(28)

式中,a、b為待定系數(shù)(分別取0.625、0.074[15]),γ為日期決定的季節(jié)角,可由式(29)計算:

γ=(n-1)·2π/365

(29)

式中,n與式(15)中含義相同。

則面元i反射的大氣輻射可由式(30)計算:

式中,Tair代表大氣溫度。

綜合式(15)～(30),可得面元i在探測方向的總反射輻射強度,如式(31)所示:

(31)

大氣中的水和CO2對紅外輻射具有較強的吸收作用,同時大氣中的氣溶膠顆粒等也會對傳輸中的輻射存在一定的散射作用,因此,對于遠距離的大氣輻射傳輸(例如:經(jīng)由海面作用產(chǎn)生的輻射投射至高空的飛機表面過程)需要考慮計算大氣透過率。

本文利用MODTRAN商業(yè)軟件對2～5 μm波段和8～14 μm波段的大氣透過率進行計算,計算條件為:飛行高度6 km,海拔高度0,探測天頂角150°,計算結(jié)果如圖(6)所示。

圖6 不同波段大氣透過率分布

3 飛機模型建立及計算網(wǎng)格劃分

3.1 飛機幾何模型

本文使用Catia商業(yè)軟件中的Imagine&Shape模塊和創(chuàng)成式設(shè)計模塊,根據(jù)實際飛機尺寸繪制飛機幾何模型,繪制結(jié)果如圖7所示,比例為1∶1。

圖7 繪制飛機模型與實際飛機模型

3.2 飛機計算網(wǎng)格

本文使用GID商業(yè)軟件對飛機幾何模型進行計算網(wǎng)格劃分,在劃分前,由于軟件兼容性的問題需要對Catia軟件中得到的飛機幾何模型進行修補,而后進行網(wǎng)格劃分,得到Triangle面元網(wǎng)格和Tetrahedra體元網(wǎng)格,劃分結(jié)果用于反射輻射強度程序的調(diào)用計算。本文反射輻射強度計算結(jié)果的后處理也采用GID商業(yè)軟件。

飛機幾何模型及計算網(wǎng)格數(shù)據(jù)見表1,修補后的飛機幾何模型及網(wǎng)格劃分結(jié)果見圖8。

表1 飛機幾何模型與計算網(wǎng)格數(shù)據(jù)表

圖8 修補后的幾何模型及網(wǎng)格劃分結(jié)果

4 結(jié)果與討論

4.1 探測角度對飛機表面紅外反射特性影響

本文首先對探測圓周角φr=270°不同探測天頂角條件下的飛機反射輻射分布特性進行了計算,計算條件如表2所示。

表2 不同探測角度下的計算條件參數(shù)

計算結(jié)果如圖9所示,圖9(d)中右側(cè)箭頭代表太陽輻射投射方向。綜合分析圖9(a)～(c)可以看到,隨著探測天頂角度的變化,飛機的主要反射輻射源由飛機右翼、兩翼和左翼逐次變化,出現(xiàn)這種現(xiàn)象是因為飛機表面幾何結(jié)構(gòu)復(fù)雜,采用BRDF模型計算時,飛機的主要輻射源出現(xiàn)在投射輻射與探測方向成鏡面反射關(guān)系的位置處。

圖9 不同探測天頂角反射輻射分布

圖10為探測天頂角θr=30°不同探測圓周角情況下的飛機表面紅外反射輻射分布計算結(jié)果。計算條件參數(shù)同表2,圖10(d)中下方箭頭代表太陽投射輻射方向。分析圖10可知,隨著探測圓周角的變化,飛機表面輻射源分布規(guī)律基本不變,集中在飛機兩翼和中部機身位置(其余位置存在遮擋效應(yīng)),但是其反射的強度由于遠離鏡面反射角度而逐漸降低。

圖10 不同探測圓周角反射輻射分布

4.2 飛機整體反射特性分析

針對特定探測平面對飛機紅外反射特性進行了計算,計算條件參數(shù)與表2相同。圖11為yoz平面不同波段不同探測天頂角下的飛機紅外反射特性計算結(jié)果,其所得輻射強度為面元輻射強度與其面積作積并在特定方向上累加的結(jié)果。

從圖11中可以看出,兩波段上下表面反射強度分布相反,這是由于在2～5 μm波段太陽輻射強于海面輻射,而在8～14 μm波段海面輻射強于太陽輻射。在圖11(a)中可以看到靠近與太陽入射方向成鏡反射方向的方向附近輻射強度明顯高于其他方向,凸顯4.1小節(jié)描述的BRDF模型中的鏡反射作用,圖11(b)中的分布則是由于海面輻射正對向飛機下表面且飛機下表面的復(fù)雜性綜合作用產(chǎn)生的。

圖11 yoz平面不同波段反射輻射強度分布

圖12為2～5 μm探測波段下不同平面不同探測角度的飛機表面紅外反射輻射特性分布計算結(jié)果,計算條件參數(shù)同表2。

圖12 2～5 μm波段不同平面反射輻射強度分布

從圖12(a)中可以看出飛機的表面反射輻射主要集中與飛機主體及兩翼,前端和后端反射很小,而從圖12(b)中可以看到飛機前端與后端的反射基本相等,且由于太陽輻射的存在導(dǎo)致上端輻射比下端輻射更強,而這中分布建立在飛機水平飛行的基礎(chǔ)上,所改變飛機姿態(tài),則分布規(guī)律會發(fā)生相應(yīng)改變。

5 結(jié) 論

本文采用有限元與商業(yè)軟件結(jié)合的方法計算了海面背景下飛機表面紅外反射特性分布情況,通過對計算結(jié)果的分析可以了解到,在BRDF模型下飛機表面鏡面反射作用對飛機的表面紅外反射特性分布影響較大,且太陽輻射與海面背景輻射在不同的探測波段分別起主導(dǎo)作用。

本文提出的計算方法并不局限于海面背景,對于不同的背景(如:沙漠、荒原或雪地等),只要獲得相應(yīng)背景的輻射模型,本文模型均有效,體現(xiàn)了本文所提出的計算方法的廣泛適用性和有效性。