柳生吉 (鐵道第三勘察設計院集團有限公司，天津 300251)

1 研究綜述

森林火災是一種突發(fā)性強、破壞性大的自然災害。林火燒毀森林及地表物，破壞了森林生態(tài)系統(tǒng)，造成生態(tài)系統(tǒng)內(nèi)生物因子、生態(tài)因子的混亂，大大降低了森林保持水土、涵養(yǎng)水源、調(diào)節(jié)氣候的作用。森林防火工作事關森林資源和生態(tài)安全，研究林火蔓延規(guī)律和模擬林火蔓延現(xiàn)象有助于森林防火。

林火燃燒是指森林中的可燃物在一定溫度下快速與氧氣結合并發(fā)光發(fā)熱的化學物理反應［1］。林火發(fā)生后，會沿著不同的方向蔓延。林火蔓延是一種林火行為，是森林可燃物從點燃開始直至熄滅的整個過程中所表現(xiàn)出的特性。林火專家將林火蔓延的影響因子主要歸結為可燃物、地形和氣象3個方面［2］。林火蔓延可視化是根據(jù)研究區(qū)的實際情況選取合適的林火蔓延模型，采用高效的計算機模擬算法和可視化技術展現(xiàn)林火發(fā)生發(fā)展的過程。

目前常見的林火蔓延算法按照對空間要素的處理方法可分為基于矢量數(shù)據(jù)形式和柵格數(shù)據(jù)形式2類算法［3］。前者即為惠更斯原理，后者中較流行的有邊界插值算法、迷宮算法以及自動元胞機算法等。矢量或惠更斯原理的特點是林火蔓延的過火區(qū)域是一個隨時間變化的連續(xù)擴展的多邊形。該多邊形的形狀由一系列過火區(qū)邊緣上的二維頂點決定［4］。為了保證一定的精度要求，隨著多邊形的增大，二維頂點數(shù)會不斷增加。在林火蔓延的過程中，每個頂點被認為是一個獨立的火點，其蔓延形狀也被認為是橢圓形，蔓延方向由風速矢量和坡度矢量疊加決定，蔓延速率則由風速、坡向、可燃物等因子通過計算林火蔓延速率模型得到。這種算法精度較高，但是計算效率較低。目前應用較多的是模擬精度略低，而計算效率較高的柵格形式算法，如美國農(nóng)業(yè)部開發(fā)的FARSITE林火管理系統(tǒng)［5－6］。國內(nèi)學者對林火蔓延可視化的研究也集中于柵格形式的算法［7］。劉月文等［8］、陳喆等［9］采用王正非模型，分別實現(xiàn)了基于二維和三維元胞自動機模型的林火蔓延模擬，而黃華國等則采用Rothermel模型實現(xiàn)了基于三維曲面元胞自動機模型的林火蔓延模擬［10］。王海軍等在元胞自動機林火蔓延模型的基礎上提出元胞作用域的概念，綜合了元胞自身屬性和鄰居元胞屬性對林火蔓延的影響［11］。楊廣斌等在動態(tài)數(shù)據(jù)驅(qū)動應用系統(tǒng)(DDDAS)技術方面開展的研究在北京市森林防火系統(tǒng)中得到了很好的應用，它體現(xiàn)了一種基于林火蔓延模型適宜性自動選擇的林火蔓延模擬的全新結構體系和技術框架［12］。DDDAS是一種功能強大的仿真應用的全新模式，它將是林火蔓延模擬發(fā)展的一個重要趨勢［13］。

隨著多維可視化技術及虛擬現(xiàn)實技術的發(fā)展，林火蔓延三維模擬及森林景觀虛擬有所發(fā)展。虛擬的林火蔓延包括三維環(huán)境林火造型模擬、地形上林火蔓延過程模擬及林火與其他景觀交互模擬等方面內(nèi)容［14］。雖然這些技術在模擬林火蔓延時還存在諸多不足，但代表了現(xiàn)代模擬林火蔓延的發(fā)展方向。三維可視化開發(fā)分為基于底層3DAPI開發(fā)和基于二次開發(fā)平臺的開發(fā)。底層3DAPI如DirectX、OpenScene-Graph等都是免費的，開發(fā)出的程序也可以脫離其運行環(huán)境，就是開發(fā)的過程比較復雜。Erdas公司的Virtual GIS，ESRI公司的3D Analyst都是較典型的二次開發(fā)三維可視化平臺，它們都脫離不了其運行環(huán)境，并且價格較昂貴。目前常見的林火造型模擬技術有紋理合成、粒子系統(tǒng)［15］以及基于物理的模擬等3種［16］。粒子系統(tǒng)是應用最廣泛的三維林火造型技術，它具有實時性好、復雜度低、真實感強、適合于大范圍等優(yōu)點。大興安嶺呼中林區(qū)林火發(fā)生頻繁，對當?shù)氐牧謽I(yè)資源造成嚴重破壞。筆者在遵循林火蔓延規(guī)律的基礎上開發(fā)出一款林火蔓延三維模擬軟件，并針對呼中林區(qū)進行初步驗證。

2 數(shù)據(jù)來源與研究方法

2．1 研究區(qū)概況 呼中林區(qū)位于黑龍江省北部，其位置如圖1所示。

圖1 呼中林區(qū)地理位置

呼中區(qū)地理坐標為 122°39'30″～124°21'00″E、51°14'40″～52°25'00″N，總面積約741 999 hm2。呼中春秋2 季受蒙古干旱風影響，天氣條件變化劇烈，常出現(xiàn)高溫、低濕和大風天氣，因而春季和秋季是林火的高發(fā)期。呼中地貌屬石質(zhì)中低山地，坡度平緩，一般小于15°。呼中植被在植物區(qū)系上屬泛北極植物區(qū)東西伯利亞植物區(qū)系，以西伯利亞植物區(qū)系成分為主，混有東北植物區(qū)系和蒙古植物區(qū)系成分。地帶性植被類型為寒溫性針葉林，以興安落葉松(Larix gmelinii)為單優(yōu)勢種。包括5種植被類型:以興安落葉松為主或混生有樟子松(Pinus sylvestris var．mongolica)、白樺(Betula ptyphylla)的興安落葉松林;樟子松林或伴生興安落葉松、白樺、山楊(Populus davidiana)的樟子松林;伴生興安落葉松和白樺的紅皮云杉(Picea koratensis)林;沿河甜楊林和沿河鉆天柳(Chosenia arbutifolia)林;白樺林和山楊林。另外，在海拔高于800 m的地方分布有可燃性很好的林下灌木偃松。

2．2 數(shù)據(jù)來源 該研究中用到的DEM是SRTM的分辨率為30 m的數(shù)據(jù)，林區(qū)地形可視化選擇Landsat5的TM影像作紋理?？扇嘉锓诸悢?shù)據(jù)利用呼中林區(qū)2000年林相圖獲取。林火發(fā)生當日的氣象數(shù)據(jù)從相關部門的氣象觀測數(shù)據(jù)獲得。

2．3 研究方法

2．3．1 林區(qū)地形三維可視化。OSG中三維場景組織是通過場景樹方式實現(xiàn)［17］，場景樹的基本單位是節(jié)點，只有一個根節(jié)點(Root)，根節(jié)點下面包括組節(jié)點(Group)和葉節(jié)點(Node)。葉節(jié)點中管理一個或者多個可繪制對象(OSG固有的和用戶自定義的)的信息，并可以通過其接口函數(shù)對可繪制體的信息進行查詢［18］。該研究利用海量地形數(shù)據(jù)處理工具VPB(Virtual PlanetBuilder)對獲取整個呼中林區(qū)的高程數(shù)據(jù)和影像數(shù)據(jù)格式化處理，將轉(zhuǎn)換后的地形數(shù)據(jù)作為OSG所組織的三維場景的一個葉子節(jié)點進行三維顯示。3DVFS中的OSG場景樹結構如圖2所示。

2．3．2 林區(qū)坡度分類。利用ArcGIS軟件Spatial Analysis模塊的Slope工具處理DEM數(shù)據(jù)，獲取研究區(qū)用于林火蔓延的坡度數(shù)據(jù)。

圖23DVFS中的OSG場景樹結構

2．3．3 林區(qū)可燃物分類。根據(jù)吳志偉等的研究［19］，應用ArcGIS將呼中林區(qū)的可燃物主要劃分為以下4類:偃松灌叢;以白樺和楊樹為主的闊葉林;分布在陰坡，以杜香和越橘為主要林下灌木的針葉林;分布在陽坡，以興安杜鵑為主要林下灌木的針葉林。在BehavePlus［20］軟件中輸入?yún)侵緜サ鹊难芯拷Y果中的相關參數(shù)，計算出不同可燃物模型的地表火蔓延速率作為3DVFS的林火蔓延初始化速率。

2．3．4 3DVFS 設計

2．3．4．1 林火蔓延模型。Rothermel林火蔓延模型以能量守恒定律等物理學理論為基礎，以林火實驗為依據(jù)，是一個半經(jīng)驗的模型［21－22］，適用范圍較為廣泛，但是該模型需要的輸入森林可燃物相關的參數(shù)較多。Rothermel模型的林火蔓延公式如下:

式中，R是蔓延速度(m/min);IR是反應強度(kJ/min·m2);ζ是林火蔓延率;φW是風速修正系數(shù);φS是坡度修正系數(shù);ρb這是可燃物體積密度(kg/m3);ε為有效熱參數(shù)(無因次量);Qt指點燃單位質(zhì)量可燃物所需要的熱量(kJ/kg)。

呼中林區(qū)是大興安嶺地區(qū)林火發(fā)生比較多的地方高發(fā)區(qū)，近些年來大量林火研究者皆在呼中展開研究，積累了大量關于該林區(qū)可燃物方面的珍貴資料。這些資料可以引用到美國Rothermel模型中做林火蔓延模擬相關的研究［23］。因此，該研究選擇Rothermel模型做為林火蔓延模型。

2．3．4．2 林火蔓延算法。該研究應用的算法是Mark A．Finney介紹的一種方法:minimum travel time(最短傳播時間算法)［24］。最短傳播時間算法是在一個與地形、可燃物等具有相同空間位置、相同空間分辨率的二維矩形柵格中進行計算。最短傳播時間指的是從火源柵格單元點直線傳播到與它相鄰的柵格單元點所需要的最短時間。

圖3 計算距著火點dY和dX處林火傳播最短時間示意

如圖3所示，對于最短傳播時間算法，林火行為特征根據(jù)主林火傳播方向、最大傳播速度和強度，以及橢圓形的火形狀尺寸等計算得到。火形狀參數(shù)和最大傳播速度決定了一場平行于地面的橢圓形狀的林火在笛卡爾坐標系下的傳播速度［25－26］:

式中，θ是指定柵格點和橢圓中心的連線與最大傳播速度方向的夾角，θ∈［－π，π］;b+c是向前最大傳播速度(θ=0);a是火翼傳播速度，著火點不一定是在上述橢圓的后交點處。當橫切線段(dX，dY)距離給定之后，就可以用公式(2)、(3)計算出θ值:

式中，α是最大傳播速度的方向;β是橫切線段與最大傳播速度方向的夾角。根據(jù)以上公式，推算出任意方向上垂直于火線邊緣的火傳播速度計算公式如下:

上面提到計算方法是基于二維柵格，而林火蔓延橢圓計算是平行于地表面的。因此，當?shù)匦尉哂衅露葧r，必須用橫切線段長度除以其移動方向坡度的余弦值來計算對應的林火傳播時間。

林火蔓延計算FireEngine模塊由4個自定義類構成:FireCell、ParameterSet、FireStart以及 FireSpread。FireCell類派生自Icomparable接口，其主要屬性有行列號、林火到該柵格的最小傳播時間，主要方法有SetValue和CompareTo，分別提供柵格單元屬性的賦值功能和與其他柵格相互比較最短傳播時間的功能。ParameterSet提供在系統(tǒng)運行過程中，存儲模擬參數(shù)的功能，它包含了所有的模擬參數(shù)。FireStart類主要包含 CreateWorkSpace、InitilizeFireMap、Readfuelmap、Read-SlopeMap、SetParamSpace等方法，分別提供創(chuàng)建工作空間、給模擬所需數(shù)組變量分配內(nèi)存、初始化林火蔓延狀態(tài)柵格圖、從坡度文件、可燃物文件讀取相應數(shù)據(jù)，將保存在電腦的模擬參數(shù)列表文件中的參數(shù)值賦給ParameterSet。FireSpread類是該模塊的核心算法部分。該研究中用來計算林火傳播的柵格有3種狀態(tài):未燃燒、點燃、熄滅，分別用數(shù)字0、1、2表示。隨著時間的增加，林火傳播柵格的狀態(tài)將按照最短傳播時間算法的規(guī)則不斷轉(zhuǎn)換。如圖4所示，該研究中最短傳播時間算法的基本步驟分為:

(1)將所有的柵格單元的最短到達時間賦值為∞。

(2)初始化起火點并計算林火傳播到緊鄰起火點的8個柵格單元所需要的時間Ti。

(3)找出所用傳播時間最短的柵格單元，將其狀態(tài)修改為“燃燒”。

(4)計算林火從起火點傳播到新引燃柵格單元再傳播到與新引燃點緊鄰的8個柵格單元(處于點燃或熄滅狀態(tài)的除外)所需的累積時間Ti。

(5)將步驟(2)算得的時間和步驟(4)算得的時間進行比較，找出傳播用時最短的柵格單元，將其狀態(tài)賦為“燃燒”，而將步驟(3)中的柵格單元的狀態(tài)賦為“熄滅”。

(6)從步驟(4)開始不斷地重復計算，直至該過程滿足模擬控制條件，運算結束。

該研究將模擬過程中不同時間段產(chǎn)生的林火蔓延狀態(tài)柵格存儲到NetCDF(Network Common Data Form)網(wǎng)絡通用數(shù)據(jù)格式文件中。

2．3．4．3 開發(fā)平臺?；?DVFS系統(tǒng)所預期的運行操作系統(tǒng)以及具體的可視化編程方面的要求，該研究選擇微軟推出的功能強大的開發(fā)環(huán)境Microsoft Visual Studio 2010作為編程工具，選擇開源的、跨平臺的OSG(OpenSceneGraph)作為三維渲染平臺，從底層構建該模擬系統(tǒng)。

3 結果與分析

3．1 3DVFS功能實現(xiàn) 3DVFS功能設計采用C#編程實現(xiàn)，界面主要劃分為菜單和工具欄操作區(qū)、模擬參數(shù)輸入和模擬結果輸出區(qū)、三維場景可視化顯示區(qū)3部分。菜單和工具欄操作區(qū)提供一些用戶常用的三維視圖操作命令，如三維可視化模型數(shù)據(jù)的導入，三維場景蔓延控制、場景漫游速度設置以及退出系統(tǒng)等。模擬參數(shù)輸入和模擬結果輸出區(qū)根據(jù)模型具體的運行需求，提供林火模擬輸入?yún)?shù)設置、啟動停止模擬功能以及模擬結束后相關統(tǒng)計結果的輸出。三維場景可視化顯示區(qū)是一個展示林火蔓延的動態(tài)過程的視口，它是一個自定義的封裝OSG相關功能的控件。3DVFS系統(tǒng)界面如圖5所示。

圖5 3DVFS界面

圖6 林火蔓延模擬各時間段火場邊界變化

3．2 呼中林區(qū)林火模擬 該研究選擇2010年6月10日發(fā)生在呼中國家級自然保護區(qū)的一場林火進行模擬。這場火的起火點位置坐標為x=493 951 m，y=574 839 0 m，實測過火面積約為580．14 hm2。模擬過程中，設置風向為45°方向(西南)，風速為12．5m/s，模擬總步長為106步，并加載可燃物類型及坡度等其他參數(shù)。模擬過程中各個時間階段火場邊界變化如圖6所示。模擬所得的火場形狀近似為橢圓，火場邊緣呈現(xiàn)鋸齒形狀，標示為黃色的柵格單元是著火點，粒子系統(tǒng)火焰表示正在燃燒的柵格單元，標示為黑色的柵格單元是熄滅點。

模擬結束后，對模擬火場面積精度進行評價。該研究選擇的評價指標為周宇飛等的計算方法［27］:

式中，SD是模擬和實測過火面積不重疊的面積，SR是實測過火面積，是相對誤差。將本實驗中實測過火面積和模擬過火面積疊合(圖7)，經(jīng)計算ε為40．46%，二者疊合程度為中等。

圖7 實測面積和模擬面積的疊合效果

4 結論與討論

(1林火燃燒是一個及其復雜的現(xiàn)象，使得林火蔓延模擬十分困難。林火蔓延模型是對林火蔓延現(xiàn)象的數(shù)學或物理上的近似，模型的簡化效果使得要非常精確的預測林火行為基本上不可能。該研究只是初步完成了一個林火蔓延模擬軟件的開發(fā)，對于該軟件的模擬結果，只是比較了實測過火面積和模擬過火面積的疊合程度，暫時無法對其預測精度進行十分準確地評價。

(2)在考慮地形坡度對林火蔓延的影響時，為了提高運算效率，研究中簡化了一些幾何方面的計算，勢必會對模擬結果的精度產(chǎn)生影響。但是，目前還沒有找到更為有效的解決方法。

(3)該研究中可燃物類型是參考一些研究人員在實驗林區(qū)的相關研究結果來粗略劃分的。實際上可燃物的狀況又是隨時間發(fā)生變化的，這里采用的是一種對可燃物狀況相對靜態(tài)的描述。所以，可燃物類型劃分對林火蔓延模擬結果影響的不確定性是比較大的，今后要探索更為精確的可燃物類型劃分方法。

(4)對于該研究開發(fā)的3DVFS系統(tǒng)，在整個一場林火蔓延模擬的過程中，假設風速以及風向是不發(fā)生變化的，但實際中風速和風向是因時因地發(fā)生變化的，這是將來需要改進的地方。當然，目前在采集實時變化的風速和風向數(shù)據(jù)方面也存在一些問題，需要在今后的研究中加以重視。

(5)該研究采用了OpenSceneGraph的粒子系統(tǒng)技術來實現(xiàn)對真實火焰的模擬［28］。在一個比較大的尺度上模擬林火蔓延時采用粒子系統(tǒng)將會面臨模擬計算用時過長對粒子系統(tǒng)動態(tài)更新的影響。針對這個問題，該研究采用將不同模擬時段林火蔓延的狀態(tài)存儲在自定義的NetCDF文件中，等模擬計算結束之后再從NetCDF讀入到計算機內(nèi)存中用于動態(tài)顯示，取得了較好的效果。

(6)目前比較常見的林火蔓延模擬很多是基于OpenGL或者像ArcGIS等一些平臺，而該研究采用的是OpenScene-Graph三維渲染引擎。實踐證明，采用OpenSceneGraph開發(fā)具有很大的優(yōu)勢，相比OpenGL開發(fā)復雜度明顯減小，工作量減少。該研究后期的目標是完成3DVFS的模擬精度驗證以及將SpeedTree樹木建模軟件［29］和OpenSceneGraph相結合，進而完成呼中林區(qū)整個森林景觀的虛擬可視化，為林區(qū)森林資源管理提供服務。

［1］鄭煥能．森林火災分類研究［J］．森林防火，1989(2):15 －17．

［2］毛賢敏．風和地形對林火蔓延速度的作用［J］．應用氣象學報，1993，4(1):100 －104．

［3］秦向東，林其釗．兩類常用森林火災蔓延模型的比較［J］．自然災害學報，2005(5):113－118．

［4］ANDERSON H E．Aids to Determining Fuel Models for Estimating Fire Behavior［M］．National Wildfire Coordinating Group，1982．

［5］田曉瑞，舒立福，王明玉．林火增長模型及應用軟件［J］．世界林業(yè)研究，2012，25(1):25 －29．

［6］WU ZW，HE HS，CHANG Y，et al．Development of customized fire behavior fuel models for boreal forests of Northeastern China［J］．Environmental Management，2011，48(6):1148 －1157．

［7］唐曉燕，孟憲宇，葛宏立，等．基于柵格結構的林火蔓延模擬研究及其實現(xiàn)［J］．北京林業(yè)大學學報，2003，25(1):53 －57．

［8］劉月文，楊宏業(yè)，王碩，等．一種基于CA的林火蔓延模型的設計與實現(xiàn) － －以內(nèi)蒙古地區(qū)為例［J］．災害學，2009，24(3):98 －102．

［9］陳喆，孫濤，張凌寒，等．三維元胞自動機各向異性林火蔓延快速模型［J］．北京林業(yè)大學學報，2012，34(1):86 －91．

［10］黃華國，張曉麗，王蕾．基于三維曲面元胞自動機模型的林火蔓延模擬［J］．北京林業(yè)大學學報，2005，27(3):94 －97．

［11］王海軍，張文婷，陳瑩瑩，等．利用元胞自動機作用域構建林火蔓延模型［J］．武漢大學學報:信息科學版，2011，36(5):575 －578．

［12］楊廣斌，劉鵬舉，唐小明．動態(tài)數(shù)據(jù)驅(qū)動的林火蔓延模型適宜性選擇［J］．林業(yè)科學，2011，47(001):107 －112．

［13］RODRIGUEZ R，CORTES A，MARGALEF T，et al．An adaptive system for forest fire behavior prediction［C］//Computational Science and Engineering，2008．CSE＇08．11th IEEE International Conference on．IEEE，2008:275 －282．

［14］葉影，齊永峰，游先祥，等．虛擬景觀下林火蔓延三維可視化研究進展［J］．林業(yè)調(diào)查規(guī)劃，2012，37(2):53 －58．

［15］趙春霞，張艷，戰(zhàn)守義．基于粒子系統(tǒng)方法的三維火焰模擬［J］．計算機工程與應用，2004，28(3):35 －37．

［16］國志興，鐘興春，方偉華，等．野火蔓延災害風險評估研究進展［J］．地理科學進展，2010，29(7):778 －788．

［17］王銳，錢學雷．OpenSceneGraph三維渲染引擎設計與實踐［M］．北京:清華大學出版社，2009．

［18］許列，韋群，王玨．基于OSG的三維場景管理及實時繪制技術研究與實現(xiàn)［J］．裝備指揮技術學院學報，2011，22(3):100 －104．

［19］吳志偉，賀紅士，胡遠滿，等．FARSITE火行為模型的原理，結構及其應用［J］．生態(tài)學雜志，2013，31(2):494 －500．

［20］ANDREWS P L．BEHAVE:fire behavior prediction and fuel modeling system －BURN subsystem，Part 1［R］．1986．

［21］ROTHERMEL R C．A mathematical model for predicting fire spread in wildland fuels［M］．Intermountain Forest＆ Range Experiment Station，F(xiàn)orest Service，US Department of Agriculture，1972．

［22］單延龍，舒立福，王洪偉，等．Rothermel火蔓延模型特征參數(shù)的解析［J］．森林防火，2003(1):22 －25．

［23］單彥龍，李華，其其格．黑龍江大興安嶺主要樹種燃燒性及理化性質(zhì)的實驗分析［J］．火災科學，2003，12(2):74 －78．

［24］FINNEY M A．Fire growth using minimum travel time methods［J］．Canadian Journal of Forest Research，2002，32(8):1420 －1424．

［25］ANDERSON D H，CATCHPOLE E A，DE MESTRE N J，et al．Modelling the spread of grass fires［J］．Journal of the Australian Mathematical Society，1982，23:451 －466．

［26］CATCHPOLE E A，DE MESTRE N J，GILL A M．Intensity of fire at its perimeter［J］．Australian Forest Research，1982，12．

［27］周宇飛，劉鵬舉，唐小明．林火蔓延模型模擬空間精度評價研究［J］．北京林業(yè)大學學報，2010，32(2):21 －26．

［28］辛佳佳．基于OSG的粒子特效仿真研究［J］．計算機光盤軟件與應用，2011(13):65－66．

［29］劉穎，羅岱，黃心淵．基于OSG的Speedtree植物模型繪制研究［J］．計算機工程與設計，2012，33(6):2406 －2409．