煙霧擴散與氣體污染的動態(tài)仿真研究

唐 勇，甄志華，汪新宇，孫旭東

(1.燕山大學 信息科學與工程學院，河北 秦皇島 066004；2.河北省計算機虛擬技術(shù)與系統(tǒng)集成重點實驗室，河北 秦皇島 066004)

0 引言

氣體污染是生活中常見的污染之一，包括燃料燃燒、尾氣排放、工業(yè)污染等，其中，煙霧擴散導致的氣體污染對日常生活造成了極大的影響，因此，模擬氣體污染具有重要的現(xiàn)實意義和廣泛的應用前景。一方面，由于煙霧隨時間變化擴散較快，煙霧細節(jié)難以捕獲；另一方面，煙霧擴散同氣體污染難以建立聯(lián)系，值得深入研究。

近年來，煙霧模擬一直都是國內(nèi)外研究的焦點話題，自2003年Stam[1]和Fedkiw等人[2]提出模擬煙霧運動的方法以來，煙霧的模擬應用越來越廣泛。采用物理模型的方法通常能真實地展現(xiàn)煙霧細節(jié)，Xie等[3]通過將低分辨率煙霧進行神經(jīng)網(wǎng)絡訓練，合成具有真實細節(jié)的高分辨率煙霧，并將其應用于不同的物理模型，但是其神經(jīng)網(wǎng)絡訓練時間較長，無法達到實時；唐勇等[4-5]提出一種改進的空間自適應漩渦限制方法，生成清晰漩渦煙霧細節(jié)；在煙霧路徑方面，通過采用改進的有限差分法求解N-S方程，加快求解速度，并引入吸引力和驅(qū)動力，實現(xiàn)大規(guī)模煙霧路徑模擬。而基于經(jīng)驗模型的模擬煙霧以及氣體污染，雖然能迅速達到實時狀態(tài)，但犧牲煙霧細節(jié)，真實感不足。陸薇等[6]提出了球形霧化渲染模型，結(jié)合傳統(tǒng)Perlin噪聲，通過HDR實時繪制層次霧；Guo等[7]采用Perlin噪聲生成異質(zhì)密度分布紋理，再利用MRF模型結(jié)合大氣散射透視圖評估渲染，繪制異質(zhì)霧圖像。

因此，本文提出一種混合經(jīng)驗模型、物理模型的方法，繪制由煙霧擴散引起的動態(tài)氣體污染。首先，針對物理模型優(yōu)化煙霧運動軌跡，并采用基于物理渲染的方式對煙霧粒繪制，提升擴散細節(jié)；其次，根據(jù)大氣污染理論建立煙霧擴散同氣體污染之間的聯(lián)系，并根據(jù)優(yōu)化后的經(jīng)驗模型，模擬具有真實感、動態(tài)可變的氣體污染。

1 基于物理模型的煙霧模擬

為模擬真實的煙霧擴散，改進半拉格朗日法中的驅(qū)動力，并優(yōu)化漩渦計算過程，結(jié)合基于物理的光照模型生成真實煙霧效果。

1.1 N-S方程構(gòu)建物理模型

采用半拉格朗日法模擬煙霧粒子能更加真實地計算煙霧擴散運動，其動量守恒方程和質(zhì)量守恒方程為

(1)

?·u=0，

(2)

其中，u為不可壓縮流體速度場，ρ為流體密度，p為壓強，v為黏度系數(shù)，f表示流體合力項，?表示梯度算子。

在流體合力項中，通過引入漩渦約束力增加煙霧模擬細節(jié)，漩渦場的計算公式為

ω=?×u，

(3)

其中，×表示卷積操作。將生成的漩渦場代入式(1)得到渦度守恒方程

(4)

根據(jù)Biot-Savart公式，通過漩渦場帶動速度場更新，實現(xiàn)真實煙霧粒子物理運動。

1.2 引入K-D樹提升計算效率

由于引入漩渦后，速度場迭代會受多渦影響，計算量顯著增大。通過引入K-D樹來降低計算量。相比八叉樹，K-D樹在空間劃分方面具有明顯的優(yōu)越性，模擬效率顯著提升。

根據(jù)K-D樹分割點設置距離閾值，如果粒子到分割點距離小于閾值，計算距離閾值內(nèi)每個漩渦對粒子的速度場影響。單漩渦對環(huán)境影響的計算公式為

(5)

為簡化計算，將漩渦視為點，由積分轉(zhuǎn)化為求和方式，速度場可近似表示為

(6)

在計算過程中，若粒子到漩渦的距離大于閾值，則視渦旋為單個渦旋結(jié)構(gòu)，而小于閾值的多個渦旋簇，采用矢量和計算渦旋位置L，

(7)

(8)

1.3 改進外力項提升運動細節(jié)

計算外力項過程中，僅靠漩渦力并不能提供更多的細節(jié)，本文通過改善風力提升煙霧擴散軌跡。

空間中的風通常由順風向Vm，橫風向Vw以及豎風向Vh組成。采用Kaimal譜作為風速譜表達式，比與高度無關(guān)的Davenport功率譜能更好展現(xiàn)風場的自相關(guān)特性。因此風速計算中脈動風通過高斯過程求解，得到v(t)表達式

(9)

其中，p為AR模型的階數(shù)，Δt為模擬風速時程的時間步長，ψk為AR模型的自回歸矩陣系數(shù)，N(t)為獨立隨機過程向量。根據(jù)隨機振動理論，利用維納-辛欽公式求得相關(guān)函數(shù)后，利用期望操作得到RN，具體計算公式為

(10)

對RN進行Cholesky矩陣分解，可求得具有時間間隔的隨機風速向量，將其帶入風力項改善煙霧運動軌跡。

1.4 煙霧粒子真實感渲染繪制

煙霧作為非均勻介質(zhì)媒體，光線穿越其中會產(chǎn)生散射和反射現(xiàn)象，如圖1所示?？紤]到實時性問題，將真實煙霧紋理結(jié)合基于物理的雙向反射分布函數(shù)渲染更真實的煙霧粒子。

圖1 光線的散射和反射Fig.1 The scattering and reflection of light

該模型中，出射輻射率L0(v)等于所有入射方向的輻射率積分和BRDF值，以及余弦值的乘積，出射輻射率的計算公式為

(11)

反射項通常由次表面散射和反射組成。次表面散射的計算公式為

(12)

其中，F(xiàn)=0.5+2r(hI)2。

采用Torrance-Sparrow微面元模型計算高光項

(13)

其中，F(xiàn)(I,h)為菲涅爾反射模型，用于處理反射光和入射光比率，G(I,v,h)為陰影遮掩函數(shù)，D(h)為法線分布函數(shù)。

2 氣體污染動態(tài)模擬

考慮到實時性，采用基于經(jīng)驗模型的方法模擬氣體污染，并建立煙霧擴散同氣體污染之間的聯(lián)系，繪制真實動態(tài)的污染。

2.1 構(gòu)建煙霧與氣體污染關(guān)系模型

污染物高度由煙云抬升高度同顆粒物高度組成，其中煙云抬升高度根據(jù)煙霧粒子高度賦值，顆粒物高度采用線性高度場，為搭建煙霧擴散與氣體污染之間的關(guān)系，將煙霧建模同大氣污染理論模型結(jié)合。本文引入高斯煙羽污染模型

(14)

其中，X(x,y,z,t,H)為時間t下，豎風向x米，橫風向y米，距離地面z米的濃度，Q代表污染源強度，H為粉塵的有效高度，u為脈動風速，σy、σz為顆粒物的水平和垂直擴散系數(shù)。增加污染因子k并優(yōu)化為指數(shù)函數(shù)，實現(xiàn)污染濃度漸變效果。針對不同污染變化，渲染成不同程度的氣體污染。

2.2 優(yōu)化污染衰減公式增加污染細節(jié)

考慮到實時性，利用基于經(jīng)驗模型的屏幕技術(shù)進行繪制。首先，根據(jù)污染濃度衰減公式獲得濃度變化梯度，計算公式為

(15)

由于氣體污染擴散的不規(guī)則性，式(15)渲染產(chǎn)生的細節(jié)較少，為此采用基于3D Perlin噪聲優(yōu)化后的分型噪聲，并在其中引入風速因子，生成動態(tài)變化的氣體污染效果，具體計算公式為

(16)

其中，m為倍頻數(shù)，f為頻率，p為振幅，vp(x,y,z,t)代表P處的脈動風場，生成真實的污染衰減公式。

2.3 動態(tài)計算環(huán)境光

優(yōu)化生成的污染濃度衰減公式，雖然細節(jié)得以提升，但顏色無法隨時間動態(tài)變化。為解決該問題，引入優(yōu)化后的時間軸算法，生成動態(tài)可變的煙霧，不同時刻環(huán)境光的計算公式為

(17)

其中，Cinc表示環(huán)境光的插值顏色，T0表示初始時間(即凌晨初始時間)，Tcurrent表示當前時間，Ttotal表示總體時間。根據(jù)負余弦函數(shù)平滑過渡，計算不同時刻的環(huán)境光CTcurrent。

另外，引入透明度系數(shù)λ，并將污染顏色同光照顏色作點乘操作，獲得最終的真實光照，如公式(18)所示，進而渲染出真實的動態(tài)氣體污染場景。

C=f(CTcurrent+Clight)+λ(1-f)Cf

(18)

3 煙霧擴散與氣體污染仿真實現(xiàn)

3.1 本文整體程序框架

本文整體程序流程圖如圖2所示。

圖2 整體程序流程圖Fig.2 Overall program flow chart

3.2 實驗結(jié)果與分析

本文實驗采用基于Windows系統(tǒng)，Unity3d平臺開發(fā)的動態(tài)氣體污染仿真系統(tǒng)，硬件環(huán)境為：Intel Core i7-4790 CPU 3.60 GHz，16 G RAM，顯卡為NVIDIA GeForce GTX 750Ti。

圖3(a)、3(b)為文獻[3]與本文方法煙霧擴散對比。圖3(a)為文獻[3]中使用時間相關(guān)的生成模型來解決流體流動問題的實驗效果圖，網(wǎng)格分辨率為256×180×180，雖然僅使用了單個時間步長，但幀率僅有0.000 8 fps，無法實時渲染。本文方法能在實時基礎上，生成大量擴散細節(jié)。圖3(c)、圖3(d)為文獻[4]與本文方法煙霧對比。圖3(c)為文獻[4]采用改進的空間自適應漩渦限制方法模擬煙霧，網(wǎng)格分辨率為64×86×64，煙霧擴散細節(jié)不夠明顯。本文利用紋理結(jié)合物理渲染模型的方法模擬的煙霧色彩、光照更加逼真。圖3(e)、圖3(f)為真實圖片與本文方法應用至大規(guī)模場景的對比，保證實時的狀態(tài)下，模擬煙霧真實自然。

圖3 與文獻[3]、[4]以及真實煙霧對比實驗Fig.3 Comparison with literature [3],[4] and real smoke

圖4展示為導彈煙霧尾跡擴散對比實驗。其中圖4(a)為真實圖片，圖4(b)為文獻[5]方法，圖4(c)為采用本文方法應用到煙霧路徑方向。通過對比看出，本文利用優(yōu)化的物理模型方法比文獻[5]生成的更自然靈活，尾跡擴散效果更明顯，更接近真實圖片效果。

圖4 導彈煙霧尾跡對比實驗Fig.4 Missile smoke wake comparison experiment

圖5為優(yōu)化后的經(jīng)驗模型模擬的氣體污染，通過將不同污染因子代入氣體污染模擬效果進行對比。圖中分別為污染因子為k=0、0.4、1下的氣體污染模擬，可以明顯比較出不同污染因子產(chǎn)生不同梯度的氣體污染，效果明顯。

圖5 污染因子對氣體污染影響實驗Fig.5 Experiment on the influence of pollution factors on smoke pollution

圖6展示了清晨T=6、正午T=13、夜晚T=20時刻下，煙霧擴散以及氣體污染隨著時間軸的動態(tài)變化，不同時刻光照插值計算顏色所展現(xiàn)出的氣體污染明暗效果。

圖6 優(yōu)化的時間軸算法實驗Fig.6 Optimized timeline algorithm experiment

為了檢測動態(tài)污染模擬效率，表1列出了本文實驗以及部分文獻數(shù)據(jù)的對比狀況。表中粒子和網(wǎng)格指文獻中模擬煙霧效果使用的粒子數(shù)或網(wǎng)格大小。其中，可以明顯看出本文算法在顯著提升模擬的細節(jié)同時，保證了實時性。

表1 不同實驗場景幀率統(tǒng)計Tab.1 Frame rate statistics for different experimental scenes

4 結(jié)論

本文提出一種煙霧擴散動態(tài)污染模擬的方法。首先，采用半拉格朗日方法計算粒子運動軌跡，并使用K-D樹提升計算效率，引入基于Kaimal譜的脈動風模型并結(jié)合真實物理光照，在避免粒子顆粒感的同時，改善煙霧擴散細節(jié)；此外，將煙霧擴散同優(yōu)化的高斯煙羽模型結(jié)合，利用改進的經(jīng)驗模型以及Perlin分型噪聲生成更加真實的氣體污染；采用優(yōu)化的時間軸算法，解決了污染顏色無法動態(tài)變化的問題，大幅度提升污染真實感。實驗數(shù)據(jù)表明，本文方法能實現(xiàn)煙霧擴散下的氣體污染實時模擬。在未來工作中，需要進一步對氣體污染與環(huán)境的交互進行研究。