目標紅外熱像的等效輻射-匯聚法數值模擬

孫海鋒, 艾 青, 夏新林

(哈爾濱工業大學 能源科學與工程學院, 哈爾濱 150001)

?

目標紅外熱像的等效輻射-匯聚法數值模擬

孫海鋒, 艾 青, 夏新林

(哈爾濱工業大學 能源科學與工程學院, 哈爾濱 150001)

為得到紅外成像陣列面上的目標輻照強度分布，提出等效輻射-匯聚法.將目標紅外成像過程分解為目標表面輻射場的產生和空間傳輸與紅外系統中的匯聚成像兩個階段.采用蒙特卡羅法模擬目標表面的等效輻射形成及其向紅外系統的傳輸過程；利用光學成像位置和光強變換關系確定目標任意點的成像方位及匯聚到成像點的輻照強度.以SDM標準飛機模型為對象，采用該方法模擬獲得了不同探測角度下的紅外熱像，并與從紅外成像陣列面發射光線的反向蒙特卡羅方法進行了計算效率的對比.計算結果表明：對于單個探測方位，等效輻射法與反向蒙特卡羅的計算效率相當；而對于多個探測方位，由于目標等效輻射(或輻射傳遞因子)的獨立性和可復用性，相對于反向蒙特卡羅法，等效輻射-匯聚法能夠提高計算效率，更便于工程應用.

紅外成像; 等效輻射; 光學匯聚; 蒙特卡羅法; 大氣衰減; 紅外傳輸; 區域分割

高速飛機等飛行器的紅外性能是其綜合性能的重要組成部分.由于探測角度和波段的差異，相對于尾噴焰紅外輻射，高速氣動熱導致的目標表面的紅外輻射不可忽視[1].相對于飛行試驗，數值模擬方法由于其成本優勢和復雜問題的易實現性，在獲取高速飛行器的紅外熱像的研究中得到了廣泛應用.一些學者針對目標的紅外輻射形成和分布特征開展了理論分析和數值模擬.如杜勝華等[2]分析了某空間飛行器的紅外光譜特性；單勇等[3]計算了亞音速和超音速條件下導彈蒙皮的紅外光譜輻射強度；劉立等[4]分析了光學頭罩在超音速條件下的熱紅外特性.這類研究通常不考慮探測系統的成像傳輸過程.而對目標在光學系統中的成像，現有的研究大多不考慮探測器口徑和焦距，采用單視點近似或小孔成像近似，通過反向蒙特卡羅(RMC)法模擬，跟蹤從像面或視點發射的光線,分析表面的方向輻射強度.如Coiro[5]采用RMC法分析了B737客機的紅外光譜成像特性；Surzhikov[6]對熱氣團的方向光譜成像特性進行了RMC模擬計算；Ross等[7]利用RMC法分析了再入飛行器的紅外熱成像.本文作者曾采用RMC方法對目標輻射到探測陣列的紅外熱像傳輸過程進行了一體化數值模擬，獲得了目標在紅外成像陣列平面的能流密度[8].該方法將目標有效輻射計算與其在紅外光學系統中的成像過程模擬緊密地聯系在一起，仿真效果好，但如果紅外探測器的參數(焦距、孔徑和探測方位等)發生變化，那么需要重新進行模擬計算，導致不必要的計算消耗.

為實現對目標表面紅外熱像特征和探測系統成像的快速預測評估，考慮到實際中的大多數紅外目標表面在探測波段內可近似為漫射面，本文將目標紅外成像過程分解為目標表面輻射的產生、空間傳輸與紅外系統中的傳輸匯聚兩個階段，提出了計算目標紅外熱像的等效輻射-匯聚法.在導出相應計算關系式的基礎上，給出了該方法的模擬過程；并以SDM飛機模型為對象，得到并分析了不同探測角度下的紅外熱像；通過與RMC法的目標紅外成像傳輸過程一體化數值模擬相比較，驗證了該方法的可靠性和計算的高效率.

1 目標紅外熱像的等效輻射-匯聚法計算原理

1.1 目標表面等效輻射及其大氣衰減的計算

如圖1所示，目標在紅外光學系統中的成像包含其自身等效輻射的形成、大氣衰減及其在光學系統中的聚集成像兩個物理過程.目標表面的等效輻射，包含其自身輻射和對投入輻射的反射或散射.對于飛機等高速飛行器而言，由于強烈的氣動熱效應，其表面的自身紅外輻射強度遠大于環境大氣和地表的紅外輻射；而太陽光主要集中在可見光波段[9]，照射到飛行器表面的紅外波段能量遠小于目標自身的紅外輻射.為分析方便，對高速飛行目標表面某處的有效輻射主要考慮該處的自身輻射和對其它部位紅外輻射的反射.

圖1 目標紅外成像物理過程

由于溫度差異導致目標自身輻射與反射輻射的黑體光譜特性不一，考慮到大氣吸收的強烈光譜性，所以二者之間以及反射輻射之間的大氣衰減率不同.為此，本研究將表面的自身輻射形成與其大氣衰減過程作為一個整體進行考量.

真實氣體的紅外傳輸特性的計算主要包含譜帶模型法[10]和逐線計算的方法[11].譜帶模型方法采用了一定程度的近似，所以精度相對較低，但是逐線計算的方法的計算量相對較大.為應用逐線計算方法，引入以下假設以降低其計算量：考慮到目標輻射的大氣傳輸距離遠大于目標的尺寸，因此可以忽略輻射在單元間相互反射中的衰減；同時光學系統入瞳面上不同點與目標的距離相差很小，所以目標上任意點輻射的大氣衰減率對于固定的探測方位可認為是確定的數值.

考慮目標形狀及溫度、輻射物性的非均勻性，將目標表面劃分為大量離散單元.對某一離散單元i，其經過大氣衰減的等效輻射強度為

(1)

(2)

式中:Rji為譜帶內單元j對單元i的輻射傳遞因子，Ai和Aj分別為單元i和單元j的面積.蒙特卡羅法采用光線隨機抽樣的方法計算單元間的輻射傳遞因子.令某一單元i發射的隨機光線的數量為Ni，跟蹤每一條隨機光線，統計系統內各個單元j吸收的份額Nij，那么單元i對單元j的輻射傳遞因子近似表示為

(3)

光線的隨機發射點、發射方向和反射方向按參考文獻[12]中的方法進行確定.在傳統方法中，光線與目標求交的計算最為耗時，因為上述計算要在光線與所有離散單元間進行.對于實際目標，其表面的幾何、物理量(溫度、材料輻射物性等)的非均勻性都需要巨大的網格數量來分辨，所以上述時間消耗更為顯著.針對于此，本研究將目標的離散單元進行分組，而求交只在光線與特定分組里的單元進行，從而減少計算量.具體采用八叉樹自適應分區方法，八叉樹的節點為長方體空間區域，每個節點可能包含8個子節點或不包含子節點，而節點是否繼續分割取決于其包含單元的數量.以圖2所示的二維情況為例進行說明，分區的過程包括創建根節點ABCD；一次分區將ABCD分為AGEI、GBIF、EICH和IFHD；二次分區將AGEI和IFHD進一步劃分.對于圖中紅色箭頭所示的光線(垂直向內)，其求交計算的過程為：光線依次與節點ABCD、AGEI、LMEQ的求交計算和其與葉節點LMEQ包含單元的求交計算.

圖2 目標八叉樹自適應分區過程

1.2 目標在紅外光學系統中匯聚與成像計算

對于如圖1所示的經過大氣衰減的表面等效輻射在紅外光學系統中匯聚成像的計算，按以下步驟進行：確定目標上任意點的成像位置；計算成像位置匯聚的等效輻射能量.

實際的成像系統可以等效為一個像方焦距與物方焦距不同的透鏡，如圖3所示.因此可以利用透鏡的成像法則計算目標的成像位置.對目標上的一點，如果其位置矢量為rp，經推導，其成像點的位置矢量rm可以表示為

式中：F和f分別為物方和像方焦距，n為透鏡平面的單位內法向量，r0為等效透鏡的中心的位置矢量.

圖3 目標與成像系統

利用光線穿越成像系統時的強度和方向變化規律可以得出像點匯聚的輻射照度.經推導得出，從目標表面上的一點rp發出，被入瞳面S匯聚到點rm的單位面積的紅外輻射能流密度Ei可以表示為

2 算例的幾何與物理模型

以圖4所示的SDM標準飛機模型為研究對象，該模型是一種三代機風動試驗模型[13]，其幾何特征為：模型弦向總長0.3 m，共包含機身、上下凸起部分、機翼、后翼和頂翼幾個部分，同時為減小尾部回流的影響，增加了一個錐形體結構.本研究僅關注模型表面的紅外特性，沒有考慮尾噴焰對其整體紅外成像特性的貢獻.

圖4 SDM模型的幾何結構

蒙皮溫度是飛機紅外圖像計算的輸入條件[14],其分布受到外部與機艙內熱環境的共同作用.本文關注點在目標的紅外成像，所以僅考慮對溫度具有顯著影響的氣動熱.以CFD軟件ANSYS CFX作為工具對氣動熱進行模擬，計算的來流條件見表1.

表1 SDM模型的來流條件

大氣中的H2O蒸汽和CO2對目標的紅外輻射具有強烈的衰減作用，大氣的光譜吸收參數利用HITRAN2008[15]光譜數據庫進行計算.大氣的溫度和壓力對譜線的線形具有很大的影響，而HITRAN數據庫僅提供了標準條件下的參數，所以非標準條件下的數值需要進行換算，具體可參見文獻[16].

為表示目標和紅外成像系統的位置關系，建立如圖5所示的全局坐標系.成像系統的方位以飛機模型中心到入瞳面中心的距離和從目標中心到入瞳面中心的向量的分量形式表示.

圖5 目標表面笛卡爾坐標和成像系統方位

Fig.5 The Cartesian coordinates system of target and the azimuth of the optical image system

3 計算結果與討論

圖6為SDM標準模型的CFD模擬計算得到的溫度分布云圖.從圖6中可以看出，機身上部和下部凸起部分的前方，由于氣動壓縮作用較為強烈，所以其表面溫度顯著地高于其他部分.除此之外，在機身后部等回流較為強烈的部分，也存在高溫區域.

以SDM模型溫度分布作為基礎，計算了8～14 μm波段模型的紅外熱像.目標表面的離散單元數為243 578個，發射率為0.3，每個單元發射的光線數為104根.在每個探測方位，探測器入瞳面中心點到目標幾何中心點的距離均為103m，入瞳面的內法向經過飛機模型的幾何中心.紅外成像系統的其他參數如表2所示.

(a) 上部

(b)下部

表2 紅外成像系統的參數

圖7為利用等效輻射-匯聚法計算的SDM飛機模型紅外熱像.從圖7中可以看出，焦平面上的能流分布的變化較為連續.在氣動熱和回流較強的區域，如上下部凸起物的前方和機身尾部，紅外信號的強度相對較大.而在機翼-機身相連的部分，由于輻射的多次反射，其紅外信號也相對較強. 在光線數為104的條件下，同時計算了發射率為0.6和0.9時的紅外熱像.圖8為不同表面發射率，(0,-1,0)方向、距離為1 km，SDM模型在像面豎直中間線上的入射能流，其中的相對值以橫坐標原點的數值為基準計算.從圖8中可以看出，目標的紅外信號隨其發射率的降低而減弱.此外，在橫軸兩端附近，信號的強度近似與發射率成正比；而在機身-機翼等連接部位，隨著發射率的降低，紅外信號的相對強度顯著增強. 為進一步評價本文提出的方法，采用RMC[8]法計算了SDM模型在表面發射率為0.3時的紅外熱像，并將兩種方法進行了對比.計算中，探測陣列的像素為500×500，每個像素發射104根光線.圖9為像面的入射輻射能流，圖10(a)為像面豎直中線上兩種方法的入射能流的對比.無論是圖7(a)和圖9間的定性還是圖10(a)中的定量對比都表明，本文提出的等效輻射-匯聚法的計算結果可信.

圖7 等效輻射-匯聚法SDM模型成像面入射能流分布(ε=0.3)

Fig.7 The radiative flux on imaging arrays of the SDM model by the equivalent radiation-concentrating method(ε=0.3)

現對等效輻射-匯聚法和RMC法的計算效率進行對比.在表面發射率為0.3、單元(或像素)發射光線數為100的條件下，經過100次重復計算，得出了兩種方法像面入射能流的統計標準差.此計算條件下二者的計算時間不相同，為此，利用上述標準差與計算時間的平方根成反比這樣的關系，以等效輻射-匯聚法為基準將RMC的標準差進行了換算.

(a)絕對值

(b)相對值

圖8 等效輻射-匯聚法目標(0,-1,0)方位、距離1 000 m時圖6(a)成像面水平中線能流密度

Fig.8 The radiative flux on the horizontal middle cross section of the imaging arrays from the (0,-1, 0) azimuth by the equivalent radiation-concentrating method

圖9 (0,-1,0) 方位RMC法SDM模型成像面入射能流密度分布

Fig.9 The radiative flux on imaging arrays of the SDM model from the (0,-1, 0) azimuth by RMC method

圖10(b)為經過換算的在相同計算時間條件下，圖7(a)中像面豎直中線上的上述標準差的分布.從圖10(b)中可以看出，二者的差別很小.二者的標準差最大值分別為1.2 W·m-2左右.由于標準差的數值與單元抽樣光線數的平方根成反比，因此可以認為光線數為104根時的標準差在0.12 W·m-2左右，即光線數為104的計算結果可認為與抽樣光線數無關.由于相同計算時間條件下兩種方法計算的能流的標準差相當，所以可以認為：對于如圖7(a)所示的單個探測方位，等效輻射-匯聚法與RMC方法的計算效率相當.同時，對于如圖7所示多個探測方位，RMC法在每個探測方位的計算時間都與圖7(a)所示方位相當，總的計算時間將為圖7(a)所示方位的4倍左右，因為RMC法將等效輻射的形成及在光學系統中的匯聚成像進行一體化計算；而對于等效輻射-匯聚法，由于等效輻射或輻射傳遞因子的可復用性，總的計算時間接近于圖7(a)所示方位的1倍，因為在該法中，等效輻射的計算時間遠大于其在紅外光學系統中成像的計算時間.

(a)能流

(b)能流標準差

4 結 論

針對漫射表面目標，從其在光學系統中成像的物理過程出發，提出了其紅外成像模擬計算的等效輻射-匯聚法.通過算例的計算與對比驗證了計算方法的準確性；同時，與反向蒙特卡洛方法相比，在多探測方位紅外成像的模擬計算中，等效輻射-匯聚法由于其中目標的等效輻射(或輻射傳遞因子)具有獨立性和可復用性，所以具有更高的計算效率.

本文的主要創新點歸納如下：

1)考慮到目標等效輻射(或輻射傳遞因子)的相對獨立性，將目標紅外成像過程分解為目標表面輻射場的產生、空間傳輸與其紅外系統中的匯聚成像兩個階段分別進行計算.

2)基于等效輻射的各個部分之間(直接發射和反射)由于發射源的溫度差異導致的大氣衰減率的差異，通過蒙特卡羅法，將等效輻射場的產生與其大氣傳輸耦合在一起進行模擬計算.

3)目標到光學系統間的距離通常遠大于目標尺寸，并且光學系統入瞳面上各點到目標的距離相差很小.基于上述假設，忽略表面輻射在相互反射過程中的大氣衰減，并利用目標點到入瞳面中心的距離采用逐線計算法計算單元直接輻射的衰減率.

4)為減少蒙特卡羅方法中光線與單元求交的計算量，利用了八叉樹區域分割方法對目標表面離散單元進行了分組.

5)對于目標表面等效輻射在紅外光學系統中匯聚成像的計算，利用光學成像變換關系確定了目標任意點的成像方位；根據成像系統中的光強傳輸變化得出了匯聚到成像點的紅外輻照強度.

[1] BARANOWSKI L C, DOUGLAS M.Surface infrared signature computer program[C]//AIAA/ASME 3rd joint thermophysics, fluids, plasma and heat Transfer Conference.Louis: AIAA, 1982: 1-8.DOI: 10.2514/6.1982-914.

[2] 杜勝華,龔加明,夏新林.某飛行器紅外輻射特性研究[J].紅外與激光工程,2008,37(增刊):432-436.

DU Shenghua, GONG Jiaming, XIA Xinlin.Infrared characteristics of a spacecraft[J].Infrared and laser engineering, 2008, 37(s):432-436.

[3] 單勇,張靖周,郭榮偉.導彈蒙皮紅外輻射特性的數值計算與分析[J].航空動力學報,2008,23(2):251-255.

SHAN Yong, ZHANG Jingzhou, GUO Rongwei.Numerical computation and analysis of the infrared radiation characteristic of missile scarfskin[J].Journal of aerospace power, 2008, 23(2): 251-255.

[4] 劉立,孟衛華,潘國慶.超音速飛行環境中光學頭罩熱輻射建模與分析[J].紅外與激光工程, 2011,40(7): 1193-1198.

LIU Li, MENG Weihua, PAN Guoqing.Modeling and analysis of infrared radiation from the dome flying at supersonic speed[J].Infrared and laser engineering, 2011, 40(7): 1193-1198.

[5] COIRO E.Global illumination technique for aircraft infrared signature calculations [J].Journal of Aircraft, 2013, 50(1):103-113.DOI: 10.2514/1.C031787.

[6] SURZHIKOV S T.Hybrid Monte-Carlo/random model of molecular lines algorithm for signature prediction[C]//44th AIAA Aerospace Sciences Meeting and Exhibit.Reno: AIAA, 2006: 1-18.DOI:10.2514/6.2006-1187.

[7] ROSS M, WERNER M, MAZUK S, et al.Infrared imagery of the space shuttle at hypersonic entry conditions[C]//46th AIAA Aerospace Sciences Meeting and Exhibit.Reno: AIAA, 2008:1-16.DOI:10.2514/6.2008-636.

[8] SUN Haifeng, XIA Xinlin, SUN Chuang, et al.Spectral backward Monte Carlo method for surface infrared image simulation[C]//International Symposium on Optoelectronic Technology and Application: Infrared Technology and Applications.Beijing: SPIE, 2014: 9300L1-9300L9.DOI: 10.1117/12.2070699.

[9] 周彥平,盧春蓮,楊莉莉,等.空間目標可見光反射特性的研究[J].哈爾濱工業大學學報,2010,42(11): 1717-1719.

ZHOU Yanping, LU Chunlian, YANG Lili, et al.Visible light reflection characteristics of space target[J].Journal of Harbin Institute of Technology, 2010, 42(11):1717-1719.

[10]王雁鳴,董士奎,談和平,等.窄譜帶模型數值研究高溫噴流動態紅外特性[J].哈爾濱工業大學學報,2009,41(11):1771-1774.

WANG Yanming, DONG Shikui, TAN Heping, et al.Numerical study on dynamic infrared properties of high temperature jet flow using narrow-band model.Journal of Harbin Institute of Technology,2009,41(11):1771-1774.

[11]ROTHMAN L S, GORDON I E, BARBE A, et al.The HITRAN 2008 molecular spectroscopic database [J].Journal of quantitive spectroscopy and radiative transfer, 2009, 110(9/10):533-572.DOI:10.1016/j.jqsrt.2009.02.013.

[12]談和平,夏新林,劉林華,等.紅外輻射與傳輸的數值計算: 計算輻射學[M].哈爾濱: 哈爾濱工業大學出版社,2006:18-19.

TAN Heping, XIA Xinlin, LIU Linhua, et al.Numerical simulation on Infrared Radiation and its Transfer:Computational Radiation[M].Harbin: Harbin Institude of Technology Press, 2006 :18-19.

[13]COULTER S M, MARQUART E J.Cross and cross-coupling derivative measurements on the standard dynamic at AEDC[C]//12th Aerodynamic Testing Conference.Williamsburg: AIAA,1982:202-214.DOI: 10.2514/6.1982-596.

[14]夏新林,艾青,任德鵬.飛機蒙皮紅外輻射的瞬態溫度場分析[J].紅外與毫米波學報,2007,26(3): 174-177.DOI: 10.3321/j.issn:1001-9014.2007.03.004.

XIA Xinlin, AI Qing, REN Depeng.Analysis on the transient temperature fields for infrared radiation of aircraft skin[J].J Infrared Millim.Waves,2007, 26(3): 174-177.DOI: 10.3321/j.issn:1001-9014.2007.03.004.

[15]劉林華,董士奎,余其錚,等.紅外1～14μm波長間隔0.1μm 上大氣平均透過率: (Ⅰ)二氧化碳的透過率[J].哈爾濱工業大學學報,1998,30(5):8-12.

LIU Linhua, DONG Shikui, YU Qizheng, et al.Atmospheric mean transmittance in wavelength interval 0.1 μ m from infrared 1 to 14 μm: (Ⅰ) transmittance of carbon dioxide[J].Journal of Harbin Institute of Technology, 1998,30(5):8-12.

[16]梅飛,江勇,陳世國,等.一種氣體吸收的逐線計算模型及其實驗驗證[J].光學學報,2012,32(3):21-27.

MEI Fei, JIANG Yong, CHEN Shiguo, et al.Experimental verification for line by line prediction model of gas absorption[J].Acta Optica Sinica, 2012, 32(3):21-27.

(編輯 楊 波)

Simulating the infrared image of targets by an equivalent radiation-concentrating method

SUN Haifeng, AI qing, XIA Xinlin

(School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001，China)

An equivalent radiation-concentrating method was developed for calculating the radiative flux of targets into sensor pixels array in an infrared (IR) detecting system.This method contains division of the imaging process: the formation and transfer of the radiation field of targets and the concentration of it in the IR optical system.In the first process, the Monte Carlo (MC) method is adopted for simulating the formation of the equivalent radiation and for simulating the transfer of it to the IR system.Then in the later process, the optical relations are employed for locating the image and for calculating the light intensity transformation.The developed method was applied for calculating the infrared image of the SDM airplane at different detecting azimuth, and the results were compared in terms of efficiency with those from the reverse Monte Carlo method that traced rays from the sensor pixels array.It is shown that for one single detecting azimuth, the developed method is closed to the reverse Monte Carlo method with respect to efficiency; for multiple detecting azimuths, compared to the reverse Monte Carlo method that coupled the generation and concentration of the effective radiation of targets, the developed method results in computing cost savings and so is more suitable in engineering applications.

infrared image; equivalent radiation; optical concentration; Monte Carlo method; air attenuation; radiation transfer; region grouping

10.11918/j.issn.0367-6234.2017.01.011

2016-02-04

國家自然科學基金(51176038, 51106036)

孫海鋒(1981—)，男，博士研究生; 夏新林(1966—)，男，教授，博士生導師

夏新林，xiaxl@hit.edu.cn

TK124

A

0367-6234(2017)01-0080-07