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

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

1 引 言

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

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

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

2 理論計(jì)算模型

2.1 Schlick模型

1982年,Schlick提出一種計(jì)算材料表面雙向反射分布函數(shù)(BRDF)的數(shù)學(xué)模型,以表征不同材料的表面反射特性,具體可由式(1)表達(dá)[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λ表示波長(zhǎng)反射因子。

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

(3)

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

(4)

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

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

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

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

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

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

建立以地球?yàn)閰⒄盏漠?dāng)?shù)刂苯亲鴺?biāo)系i-j-k,規(guī)定i軸方向?yàn)檎嫦蛭?k軸方向?yàn)榈孛娣ㄏ?根據(jù)右手法則確定j軸方向,坐標(biāo)原點(diǎn)為飛機(jī)所在位置P,如圖1所示。

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

2.2.2 飛機(jī)坐標(biāo)系

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

建立飛機(jī)坐標(biāo)系x-y-z,定義x軸方向?yàn)橹休S線向后,y軸方向?yàn)樗矫嫦蛴?右手定則確定z軸方向,如圖2所示。

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

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

(5)

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

(6c)

2.2.3 面元坐標(biāo)系

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

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

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

(1)u軸:P1為原點(diǎn),向量P1P2方向?yàn)閡軸方向,其方向余弦可由式(7)表示:

(7)

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

(8)

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

(9)

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

(10)

三個(gè)坐標(biāo)系確定后,結(jié)合式(5)可以確定飛機(jī)表面上某點(diǎn)在當(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 輻射遮擋判斷因子計(jì)算

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

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

(12a)

(12b)

在考慮其他面元的遮擋時(shí),記投射輻射間接遮擋因子為Rs2和探測(cè)間接遮擋因子Rd2?？紤]面元j對(duì)面元i的遮擋情況,如圖4所示。

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

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

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

(13)

2.2.5 紅外探測(cè)模型

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

圖5 紅外探測(cè)模型

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

(14)

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

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

2.3 飛機(jī)表面投射輻射模型

2.3.1 太陽(yáng)輻射

太陽(yáng)赤緯δ0由飛行時(shí)的日期決定,可由式(15)計(jì)算:

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

太陽(yáng)時(shí)角t0由飛行時(shí)刻決定,可由式(16)計(jì)算:

(16)

式中,t為一天中所處的時(shí)刻。

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

(17)

式中,ω為飛機(jī)所在緯度。

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

(18)

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

(19)

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

面元i反射的太陽(yáng)直射輻射強(qiáng)度由式(20)計(jì)算:

(20)

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

(21)

式中,Psc為太陽(yáng)輻射中的散射比例,取0.25。

面元i反射的太陽(yáng)散射輻射可由式(22)計(jì)算:

(22)

2.3.2 海面反射的太陽(yáng)輻射

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

(23)

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

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

(24)

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

面元i對(duì)海面反射的太陽(yáng)輻射的反射輻射強(qiáng)度由式(25)計(jì)算:

(25)

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

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

(26)

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

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

(27)

2.3.4 大氣輻射

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

(28)

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

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

(29)

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

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

式中,Tair代表大氣溫度。

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

(31)

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

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

圖6 不同波段大氣透過(guò)率分布

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

3.1 飛機(jī)幾何模型

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

圖7 繪制飛機(jī)模型與實(shí)際飛機(jī)模型

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

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

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

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

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

4 結(jié)果與討論

4.1 探測(cè)角度對(duì)飛機(jī)表面紅外反射特性影響

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

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

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

圖9 不同探測(cè)天頂角反射輻射分布

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

圖10 不同探測(cè)圓周角反射輻射分布

4.2 飛機(jī)整體反射特性分析

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

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

圖11 yoz平面不同波段反射輻射強(qiáng)度分布

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

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

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

5 結(jié) 論

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

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