柴油燃燒煙塵干擾的紅外輻射特性研究

田 益，馬向超，蔣佳麗

（西安電子科技大學，陜西西安 710071）

0 引 言

在現代軍事戰爭中，柴油燃燒所形成的煙塵是非常普遍的。例如：坦克、裝甲車等重型武器運行過程中，柴油經過柴油機燃燒后會釋放出大量的煙塵；加油站和油井等能源設施也經常被作為軍事打擊的目標；在某些特殊情況下，為了隱蔽目標，有時也會故意點燃油井、油桶等形成大量燃燒煙塵，以減少目標被命中的幾率，如海灣戰爭中伊拉克故意點燃油井躲避美軍紅外武器的攻擊。

隨著紅外成像導引頭技術［1］的快速發展，雖然采用第三代凝視紅外焦平面探測器的導引頭系統綜合性能大幅提高，但是燃燒煙塵等擴展干擾物依然是制約紅外制導導彈作戰效果的主要措施之一。當導引頭和目標之間存在燃燒煙塵時，由于其消光遮蔽作用，可破壞導引頭對目標的發現、識別和跟蹤，使得處于搜索狀態的導彈不能轉入跟蹤狀態；處于跟蹤狀態的導彈則會丟失目標，由跟蹤轉入搜索，或者導致其跟蹤誤差增大。因此，柴油燃燒所形成的煙塵不僅是己方裝備對抗紅外成像制導武器的潛在手段之一，也是己方裝備的導彈紅外導引頭所面臨的潛在干擾之一。

目前對柴油燃燒煙塵的研究主要在于對煙塵的成分進行物化分析，如柴油完全燃燒時的主要產物是二氧化碳和水，不完全燃燒時還產生一氧化碳甚至碳微小顆粒。

由于柴油的燃燒煙塵干擾普遍存在且容易釋放，本文擬通過對這種煙塵的物理化學特性、時空運動特性及其紅外輻射的空間、時間和光譜傳輸特性等進行研究，構建燃燒煙塵的擴散運動和紅外輻射傳輸特性模型，為紅外制導武器抗燃燒煙塵干擾能力的評估提供基礎。

1 柴油燃燒煙塵的收集和微觀形貌分析

柴油燃燒煙塵的組成成分較為復雜，受燃燒條件的影響很大。我們使用煤油燈對柴油進行直接可控燃燒，然后采集其燃燒煙塵。

直接采集煙塵所使用的方法如圖1 所示，其目的是為了獲得用于微觀形貌分析的顆粒物。首先，將導電膠帶粘貼到電鏡樣品臺上，靠近燃燒煙塵排放通道，使煙塵顆粒直接流經導電膠帶被沾附于樣品臺，從而獲得燃燒煙塵樣品。

圖1 燃燒煙塵采集示意圖Fig.1 Schematic of smoke collection from diesel combustion

圖2為所收集的燃燒煙塵樣品的掃描電子顯微鏡（scanning electron microscope，SEM）照片，顯示了500 nm 標尺下的表面形貌。從圖2 可以看出，燃燒煙塵主要是由小尺寸的顆粒組成，小顆粒物的形狀近似為球形，尺寸在20～40 nm之間。

圖2 柴油燃燒煙塵的SEM圖（標尺500 nm）Fig.2 SEM image of smoke from diesel combustion（scale ruler 500 nm）

2 紅外輻射傳輸模型

2.1 單個粒子消光特性

在建立煙塵的紅外輻射傳輸模型［2］前，需要知道單個煙塵粒子的消光特性參量，如散射截面σsca、消光截面σext、不對稱因子g等。為了簡化，根據實驗分析結果，我們將煙塵粒子視為球形粒子。對于球形粒子，這些參量可由Mie散射理論求解得到，也可采用離散偶極近似（discrete dipole approximation，DDA）算法計算。考慮到通用性，本文采用DDA 算法來計算單個煙塵粒子在特定波長下的消光特性。

當波長不同時，粒子的復折射率會發生改變，因此粒子對輻射的消光特性會隨著波長變化。Chang等［3］對煙塵顆粒提出了一組對數擬合函數，復折射率的表達式為n′＝n＋ik，n和k的擬合函數如下所示。

式中：λ為入射光波長，單位μm。

根據對燃燒煙塵微觀形貌的分析，球形粒子半徑取20 nm，依據Chang 等提出的復折射率模型，在紅外波段上等間隔選取m個采樣點來計算對應波長下單個煙塵顆粒的散射截面、消光截面以及不對稱因子，最終得出煙塵粒子在采樣點處的光學參量。

2.2 粒子系消光特性

燃燒煙塵是由多個粒子組成的彌散粒子系。根據粒子復折射率及尺度參數計算出單個粒子消光截面及不對稱因子后，再結合粒子系的物理參數（如粒子濃度、形狀及粒徑）可求出粒子系的消光及散射特性。

為了簡化模型，本文將粒子系中的粒子視為相同粒徑的球體，則單個粒子在某波長處消光截面σext乘以網格區域單位體積內粒子數N的結果即為粒子系在該波長下的消光系數βsys：

由于在單個顆粒的消光中，有一部分是前向散射，依然會進入視場，因此，需要考慮前向散射的貢獻。Henyey 和Greenstein 根據多年研究，提出了一個與已知實驗數據較為符合的散射相函數經驗公式，簡稱為HG相函數［4］，如式（4）所示。

式中：g是一個無量綱的控制散射相函數分布的常數，為單個煙塵顆粒的不對稱因子。當g＜1時，散射是各項異性的，g越大，分布越集中，前向散射占據的比率越大；當g＝1時，散射方向完全集中在前向。

單個網格內粒子系在各角度散射的能量比率為

式中：∑σsca是流場計算中在某波長下單個網格內粒子的累加散射截面；Δs是垂直入射方向的單個網格橫截面積；PHG（Θ）為該波長處單個煙塵顆粒HG 的散射相函數。當Θ 取0 時，F（Θ）即為單個網格內粒子系的前向散射比率。

2.3 粒子系光譜消光特性

燃燒煙塵顆粒的光學特性具有強烈的光譜依賴性，粒子系在某個波段內的消光特性通常是對光譜平均獲得的，其中常用的平均方法是普朗克平均法，即對之前計算的粒子系在該特定波段內的光譜消光特性結果進行普朗克（Planck）平均，普朗克平均消光系數為

式中：f1為在紅外波段λ1～λm上分別等間隔選取m個采樣點后計算采樣點對應的消光系數βsys，并將βsys數值進行擬合而得的4 階多項式函數；λ1和λm分別為波段積分下限與上限；Ebλ為特定波長范圍內等間距所取的各個波長對應的黑體的光譜輻射照度。

同理，普朗克平均前向散射比率可表達為

式中：f2為利用在紅外波段λ1～λm上等間隔選取m個采樣點后計算采樣點對應的前向散射比率F（Θ），并將F（Θ）數值進行擬合而得的4階多項式函數。

上述計算過程可用圖3所示流程圖表達。

圖3 粒子系在不同波段的消光特性及前向散射計算流程圖Fig.3 Flow chart for calculating the extinction properties and forward scattering of particles in different wavebands

2.4 出射強度

沿厚度方向將煙塵平均分割為Ω份，每份的厚度為Δx，如圖4所示。

圖4 煙塵平均分割為Ω份的示意圖Fig.4 Schematic of dividing the smoke into N parts

出射強度為

其中：I0為入射紅外輻射強度。

3 時空擴散模型

3.1 物理模型

燃燒煙塵的擴散運動符合流體的運動規律，因此可以基于流體力學理論對煙塵的時空擴散運動進行建模。基于物理過程的煙塵擴散運動仿真的理論基礎是求解其對應的Navier-Stokes（NS）方程，可使用半拉格朗日法與壓力泊松方程來高效地求解方程，并使用渦量約束法［5］來降低由半拉格朗日法引起的數值耗散問題，從而準確地模擬出煙塵流場的湍流現象。

一般而言，煙塵的黏性可忽略不計。當煙塵的擴散速度遠小于聲速時，可壓縮性的影響也可以忽略不計，因此視煙塵為無黏性、不可壓縮流體［6］。求解煙塵流場的方程可表達為［6］

式中：ρ為密度；p為壓強；a為體積力（如重力，慣性力）的加速度；u為速度。

流場中煙塵顆粒受到的浮力與溫度有關，溫度較高的煙塵顆粒受到浮力的作用上升，通常用如下經驗公式描述［7］。

式中：α表示溫度對浮力的影響系數；T為煙塵溫度；Tair為環境溫度；y表示浮力方向上的單位向量。

溫度在空間分布中不是均勻的，導致浮力的空間分布也不均勻，所以該浮力對速度場的影響不會被壓力項完全消除，會形成對流。

通過對溫度的計算可進而求出煙塵顆粒所受的浮力大小與煙塵的自發輻射強度，為此需要引入一個溫度標量場。溫度場不僅自身具有擴散效應，速度場也會對其產生對流影響。溫度場可表示為

式中：cT為溫度冷卻系數；Tmax為流場最高溫度；Tair為環境溫度。

風與大氣條件和海拔高度有關，流場中用于描述風場的公式為

式中：u0是處于海拔為z0時測量出的風速；z是地面海拔；r用于描述大氣中風速隨海拔高度的變化情況。在實際仿真計算中取z0為10 m，并設r的默認值為0.3。

3.2 基于網格的煙塵流體仿真方法

不可壓縮流體控制方程可化為拉格朗日形式和歐拉形式，兩種形式在理論上是等價的。在CFD 中，由于計算精度有限，且兩種形式的離散化方式不同，會形成兩種截然不同的仿真方法：基于粒子的方法、基于網格的方法。

本文采用的流體計算方法是歐拉法，將流場計算區域劃分為離散網格。一個三維空間的規則網格的數據結構就是一個三維數組，只需要記錄數組的大小，通過下標就能完成對數組元素的訪問。本文的流場計算采用規則的靜態結構化網格來描述。依據流場維度將計算區域劃分為X×Y×Z的三維立方體網格。

網格將空間均勻地劃分，仿真過程中速度、溫度、壓強、密度這些物理量都記錄在網格單元中心。為了簡化流場邊界處理，在實際仿真過程中需要在流場周圍額外加上一層網格［8］。因此在程序中，對每個物理量分配的數組大小為（X＋2）×（Y＋2）×（Z＋2）。

3.3 紅外煙塵渲染

紅外煙塵，特別是軍事上用于遮擋目標而使用的紅外煙幕，發煙器都會具有比環境溫度高的初始溫度，故需要考慮煙塵的紅外自發輻射。渲染紅外煙塵主要需要考慮的因素是透過率和對輻射的散射。煙塵粒子對光子的相互作用概率決定了其透過率，從視線方向穿透煙塵體的總透過率就是最終被渲染出來的α通道，α通道即在OpenGL 渲染引擎內用于存儲透明度值的通道。煙塵的消光系數與消光截面和密度相關，輻射一部分會被外圍的煙塵吸收和散射，另一部分會透射出來并進入相機［9］。本文只考慮在視線積分路徑上的前向散射。

由于煙塵運動模型采用網格方法構建，這種網格無法通過常規的光柵化方法進行渲染，故采用的是體繪制算法［10］，該算法能夠高效地處理煙塵體在視線方向上的透過率和前向散射，以達到實時性的需求。

本文采用了體繪制中最常用的光線投射算法，如圖5 所示，展示了從相機視口發出的射線穿過體紋理并進行采樣的過程。首先需要計算煙塵占據區域，當從相機引出的射線命中煙塵前表面時開始進行體繪制，光線傳輸方向與射線方向相反：射線不斷地在煙塵體積內步進，不斷累計亮度和透過率，當射線離開長方體后表面，體繪制結束，其過程中累計的亮度和透過率就是煙塵的亮度和透過率。

將圖5 射線穿越體紋理的過程作為采樣合成過程，則合成公式可表達為

圖5 射線穿過體紋理［11］Fig.5 Rays passing through the volume texture

式中：Ci和Ai分別是在體紋理上采樣所得到的顏色值和不透明度，其實也就是體素中蘊含的數據；和表示累加的顏色值和不透明度。

基于本文所建立的紅外輻射傳輸模型，對體紋理的顏色值Ci與不透明度Ai有

式中：Iself為網格內煙塵的自發輻射強度；Itrans為透射進來的輻射強度；對于煙塵自身紅外輻射，將其簡化為灰體輻射模型，且使其輻射強度與濃度成正比，該網格單元發出的自身紅外輻射強度為

式中：Ibλ為T溫度下λ波長的黑體輻射強度；λ1～λ2為所計算輻射的波段；ε為由煙塵粒子的等效發射率；N為網格內的粒子數。

3.4 具體實現方法

針對煙塵的擴散運動模型，基于NS 流體力學控制方程，采用歐拉網格方法求解流場的擴散運動，主要包括：求解壓力項、計算對流、降低流場渦量和其他內外力項的解算。將流場計算區域分為40×40×40個網格，對每個網格內的煙塵運動狀態進行計算迭代，最終可得出網格內的煙塵的速度、濃度、溫度等物理量。

針對煙塵的自發紅外輻射，首先建立煙塵的等效發射率模型，然后基于灰體模型建立煙塵粒子系的自身紅外輻射模型。針對煙塵紅外輻射的時間、空間、波段紅外輻射傳輸特性，采用DDA 算法，且分別在3～5 μm 以及8～12 μm 紅外波段等間隔選取40 個采樣點，計算單個煙塵粒子在不同波長處的光譜和波段光學參量，如散射與吸收截面、散射相函數等；繼而可計算粒子系在各波長采樣點處的消光系數和前向散射比例；然后，基于普朗克平均方法獲得煙塵粒子系在3～5 μm、8～12 μm 紅外波段的消光和散射特性；最后，采用對煙塵仿真區域進行等間隔采樣的光線投射算法，結合單個網格內粒子系的消光、散射和自發輻射模型，從成像平面的每個像元引出一條射線，計算途經每一條射線到達視點處的紅外輻射，如圖6所示。

圖6 煙塵輻射傳輸計算示意圖［12］Fig.6 Schematic for calculating the radiation transfer of smoke

依據以上一系列理論和方法，我們使用OpenGL圖形API 與CUDA 并行程序搭建整個煙塵的計算和渲染的框架，實現了基于歐拉網格方法的煙塵擴散運動和紅外輻射傳輸特性的聯合解算。

4 結果與分析

煙團的運動過程主要由煙塵源（煙塵發生器）類型和邊界條件決定，運動導致煙塵的溫度和密度隨時間變化，所以在某一時刻的溫度和密度的分布決定了煙塵的亮度和透過率。下面將結合仿真結果分析擴散時間、風場邊界條件以及波段對煙塵運動的影響，并檢驗其對恒定溫度黑體的遮擋效果。

4.1 風場對仿真結果的影響

圖7 顯示的是在5 m×5 m×5 m 的空間內，無風條件下不同時刻煙的擴散圖。發射源持續釋放初始速度為1 m/s 且速度方向垂直向上的煙塵，煙塵初始濃度為0.1 kg/m3。發射源溫度恒定為350 K，環境溫度設為280 K，并且在流場區域放置一個溫度與發射源相同的黑體，在每個時刻對煙塵的擴散情況進行計算仿真。

圖7 無風情況下不同時刻煙的擴散圖Fig.7 Spread of smoke at different times without wind

將煙塵源位置設置為側面，煙塵初始釋放速度設為3 m/s 且方向設為水平方向，其它初始條件保持不變。在有風情況下煙塵的時空運動如圖8所示。

圖8 有風情況下不同時刻煙的擴散圖Fig.8 Spread of smoke at different times in the wind

由仿真結果可以看出，在風力作用下，煙塵的整體形態發生了變化，向著風場方向偏移。一般而言，風速越大，偏移越明顯。經與無風情況進行對比還可看出，在有風與無風的條件下，煙塵對350 K 方形黑體遮蔽能力不同：在有風的情況下，煙塵迅速擴散，這導致對黑體的遮蔽時間和效果都明顯下降。

4.2 不同波段下的仿真結果

設煙塵初始釋放速度為1 m/s，釋放濃度為0.1 kg/m3。量化范圍上限設置為350 K所對應的黑體輻射亮度，下限設置為280 K。波段分別取3～5 μm 和8～12 μm，并在場景中加入與煙塵初始溫度（350 K）相同的黑體方塊，用于觀察煙塵的遮蔽特性，得到仿真結果如圖9所示，左側是8～12 μm波段，右側是3～5 μm。

圖9 不同波段紅外煙塵Fig.9 Smoke in different infrared bands

將3～5 μm 和8～12 μm 波段中煙塵對黑體的遮蔽情況進行對比，可知：在8～12 μm 波段內輻射透過率較高，煙塵的遮蔽能力較弱。紅外圖像中煙塵底端比較暗，這是因為底部煙塵濃度較高，輻射衰減很快；而在紅外煙塵頂端，煙塵濃度與溫度較低，輻射衰減較慢。

5 結束語

由于柴油燃燒煙塵在戰場中對目標的遮蔽特性，本文建立了柴油燃燒煙塵的紅外仿真模型。通過獲取煙塵形貌特征、建立紅外輻射傳輸模型、建立煙塵時空擴散模型等，完成了對紅外煙塵的實時仿真渲染，研究了不同風場下的柴油燃燒煙塵干擾的紅外輻射特性。研究結果可為紅外制導武器抗燃燒煙塵干擾能力的評估提供基礎。