復雜環境背景下典型飛行器表面紅外反射特性研究

孫建平,齊 宏,王申領,趙松慶,郝燕云,吳根水

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

1 引 言

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

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

飛行器的紅外目標探測以其表面的紅外反射特性分布規律為基礎,而飛行狀態與飛行背景對飛行器表面反射特性具有重大影響,因此,本文采用有限元與商業軟件結合的方法對不同飛行狀態復雜背景下飛行器的表面紅外反射特性的分布規律進行研究,以為飛行器的紅外目標探測提供參考。

2 理論計算模型

2.1 Schlick模型

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

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

(1)

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

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

(2)

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

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

(3)

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

(4)

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

2.2 坐標系轉換及輻射遮擋判斷因子計算

隨著飛機的運動,飛機與探測器的相對位置一直在變化,同時對于飛機上單個面元來說,入射輻射方向與反射輻射方向也是不斷變化的。因此,需要建立當地坐標系-飛機坐標系和飛機坐標系-面元坐標系之間的坐標轉換關系,以此得到探測器接收到的輻射信號與飛機面元之間直接的對應關系。

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

2.2.1 當地-地球坐標系

飛機所在地的太陽輻射初始入射方向取決于當地時間、當地緯度和季節等參數。飛機緯度信息則取決于飛機與地球的相對位置,因此,需要建立當地-地球坐標系描述這一關系。

建立以地球為參照的當地直角坐標系i-j-k,規定i軸方向為正面向西,k軸方向為地面法向,根據右手法則確定j軸方向,坐標原點為飛機所在位置P,如圖1所示。

圖1 當地-地球直角坐標系

2.2.2 飛機坐標系

飛機在飛行過程中存在旋轉和仰角變化,且其飛行高度和飛行速度同樣變化劇烈,其運動姿態和運動狀態的多變性導致飛機表面獲得的投射輻射方向不斷變化。同時,由于飛機與探測器相對位置不斷改變,探測器的探測方向也不斷改變。因此,需要建立飛機坐標系,并明確當地坐標系-飛機坐標系之間的坐標轉換關系。

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

圖2 飛機直角坐標系

則當地坐標系中坐標(i,j,k)與飛機坐標系中坐標(x,y,z)之間的轉換關系如式(5)所示:

(5)

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

(6c)

2.2.3 面元坐標系

BRDF的計算是在每一個飛機表面的面元上進行的,即需要得到每一個面元的入射輻射和反射輻射的方向參數,因此需要將飛機表面的入射輻射方向轉化為每個面元的入射輻射方向參數,因而需要建立面元坐標系,并明確當地-飛機-面元之間坐標位置關系。

圖3 面元直角坐標系

面元坐標系如圖3所示。圖中P1、P2、P3和P4組成了一個四面體網格,P1、P2和P3為飛機表面面元的三個頂點,P4為飛機內部節點,以此為基礎,建立面元坐標系u-w-v,規則如下:

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

(7)

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

(8)

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

(9)

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

(10)

三個坐標系確定后,結合式(5)可以確定飛機表面上某點在當地坐標系中的坐標(i,j,k)與其所在的面元坐標系的坐標(u,w,v)之間的轉換關系,即式(11):

(11)

2.2.4 輻射遮擋判斷因子計算

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

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

(12a)

(12b)

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

圖4 面元遮擋關系示意圖

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

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

(13)

2.2.5 紅外探測模型

在近距離紅外探測時,對于相同探測方向,飛機上不同面元相對于探測器的探測方向并不相同,因此需要建立探測模型,將以當地坐標系為基準的探測器探測方向轉換為飛機坐標系為基礎的探測方向。考慮100 m探測距離模型,如圖5所示。

圖5 紅外探測模型

探測方向天頂角為θd,圓周角為φd,當地坐標系中探測點位置Pd(id,jd,kd)經過坐標轉換可得飛機坐標系探測點位置Pdp(xd,yd,zd)和所在面元重心位置PGp(xGp,yGp,zGp),進而可以得到以飛機坐標系為基準的每個面元的探測方向,如式(14)所示:

(14)

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

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

2.3 飛機表面投射輻射模型

2.3.1 太陽輻射

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

式中,n為從1月1日起已經度過的天數。

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

(16)

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

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

(17)

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

當地坐標系中太陽入射輻射的天頂角θ和圓周角φ可由式(18)計算:

(18)

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

(19)

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

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

(20)

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

(21)

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

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

(22)

2.3.2 海面反射的太陽輻射

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

(23)

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

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

(24)

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

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

(25)

2.3.3 海面自身發射輻射

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

(26)

式中,sea代表海面自身發射的輻射。

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

(27)

2.3.4 大氣輻射

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

(28)

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

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

(29)

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

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

式中,Tair代表大氣溫度。

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

(31)

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

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

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

3 飛機模型建立及計算網格劃分

3.1 飛機幾何模型

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

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

3.2 飛機計算網格

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

飛機幾何模型及計算網格數據見表1,修補后的飛機幾何模型及網格劃分結果見圖8。

表1 飛機幾何模型與計算網格數據表

圖8 修補后的幾何模型及網格劃分結果

4 結果與討論

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

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

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

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

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

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

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

4.2 飛機整體反射特性分析

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

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

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

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

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

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

5 結 論

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

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