窄帶水平集法下林火蔓延仿真研究

翟春婕，曹兆樓(1.南京森林警察學(xué)院森林消防系，南京,1003；. 南京信息工程大學(xué)光電工程系，南京,10044)

窄帶水平集法下林火蔓延仿真研究

翟春婕1*，曹兆樓2
(1.南京森林警察學(xué)院森林消防系，南京,210023；2. 南京信息工程大學(xué)光電工程系，南京,210044)

火焰鋒面的追蹤是建立林火蔓延模型及預(yù)測鋒面位置的關(guān)鍵環(huán)節(jié)。基于水平集法追蹤鋒面具有求解靈活穩(wěn)定、易于處理鋒面拓?fù)浣Y(jié)構(gòu)變化的特點(diǎn)。針對傳統(tǒng)水平集法重新初始化過程復(fù)雜、計(jì)算量較大的問題，根據(jù)距離函數(shù)定義，提出通過形態(tài)學(xué)膨脹的方法以鋒面點(diǎn)為圓心向外膨脹，獲得鋒面附近窄帶內(nèi)每個(gè)點(diǎn)距離鋒面的最近距離，重建窄帶內(nèi)的水平集函數(shù)，物理意義明確，有效減小了計(jì)算量，在此基礎(chǔ)上基于林火速度場經(jīng)驗(yàn)?zāi)Ｐ湍M了單火源及雙火源的鋒面蔓延情況。結(jié)果表明，該方法能夠有效模擬坡度、風(fēng)速等因素對火焰鋒面的影響，為其在林火蔓延模型中的應(yīng)用提供了基礎(chǔ)。

鋒面追蹤；林火蔓延；水平集；形態(tài)學(xué)

0 引言

林火蔓延模型通過有效預(yù)測林火行為從而幫助防治火災(zāi)，提高撲火效率，指導(dǎo)林火管理工作，有巨大的社會經(jīng)濟(jì)效益，受到了廣泛的關(guān)注。由于實(shí)際林火過程影響因素眾多，如坡度、風(fēng)向及可燃物含量等，盡管近年來計(jì)算機(jī)的計(jì)算能力增長迅速，對林火動力學(xué)理解也愈加深刻，但從基本物理原理出發(fā)同時(shí)建模求解由大尺度地形決定的氣候影響、小尺度地形決定的局部燃燒動力學(xué)影響以及瞬態(tài)湍流等因素仍非常困難，難以實(shí)現(xiàn)林火蔓延的實(shí)時(shí)預(yù)測。因此目前林火蔓延軟件一般使用經(jīng)驗(yàn)?zāi)Ｐ?，通過火蔓延速率與方向總結(jié)規(guī)律，如澳大利亞的McArthur模型[1]、美國的Rothermel模型[2]、加拿大的國家林火蔓延模型[3]及我國的王正非模型[4]等。

林火蔓延模型在速度場建模之外還包括給定速度場時(shí)火焰鋒面演進(jìn)追蹤，針對火焰鋒面的演進(jìn)追蹤，很多學(xué)者提出了有效方法。傳統(tǒng)方法一般使用拉格朗日粒子法，比如FARSITE[5]，通過連接相鄰的粒子可獲得鋒面，但該方法計(jì)算量與火焰尺寸相關(guān)且當(dāng)鋒面拓?fù)浣Y(jié)構(gòu)發(fā)生相交變化時(shí)需要特殊處理。相比之下，水平集法通過求解雙曲型偏微分方程實(shí)現(xiàn)鋒面的演進(jìn)，求解結(jié)果經(jīng)過多年發(fā)展具有靈活穩(wěn)定，且易于處理拓?fù)浣Y(jié)構(gòu)變化的特點(diǎn)。近年來已有學(xué)者將其應(yīng)用于林火蔓延模型，如Mallet[6], Mandel[7], Kim[8], Rehm[9]等，實(shí)現(xiàn)了火蔓延鋒面隨時(shí)間演進(jìn)的追蹤。Lautenberger[10]基于數(shù)據(jù)融合的思想，通過水平集法追蹤鋒面成功獲得了速度模型中的各個(gè)參數(shù)。Liu等[11]將其應(yīng)用于預(yù)混火焰中追蹤鋒面，實(shí)現(xiàn)了預(yù)期的效果。但水平集法是通過設(shè)定距離函數(shù)將鋒面信息隱藏在高維空間之中，且在追蹤過程中需不斷重新初始化距離函數(shù)以保證其滿足定義，運(yùn)算量較大。由于林火蔓延模型中僅關(guān)注火蔓延鋒面的位置，其它區(qū)域水平集函數(shù)的變化并不影響結(jié)果，因此有學(xué)者提出窄帶水平集法，僅計(jì)算鋒面附近而非整個(gè)火焰區(qū)域距離函數(shù)的變化，顯著減小了計(jì)算量[12]。雖然這些工作為水平集法更廣泛的應(yīng)用提供了基礎(chǔ)，但距離函數(shù)的重新初始化仍需求解時(shí)變偏微分方程達(dá)到穩(wěn)定狀態(tài)，時(shí)間復(fù)雜度較大，影響了模型的實(shí)時(shí)性。

本文根據(jù)目前研究現(xiàn)狀，提出結(jié)合速度場經(jīng)驗(yàn)?zāi)Ｐ图八郊A(yù)測火蔓延鋒面位置，針對距離函數(shù)重新初始化計(jì)算量較大的問題，提出根據(jù)距離函數(shù)定義使用形態(tài)學(xué)膨脹的方法重建鋒面附近窄帶內(nèi)的距離函數(shù)，減小計(jì)算量，并數(shù)值模擬了單點(diǎn)火源及雙點(diǎn)火源在不同外部條件下火焰鋒面的蔓延情況，驗(yàn)證了本文提出的方法。

1 基本原理與方法

林火蔓延模型包含兩個(gè)部分：速度場模型及火蔓延鋒面追蹤。前者根據(jù)鋒面位置確定局域的環(huán)境如可燃物濃度，風(fēng)向與鋒面法向夾角等進(jìn)而獲得蔓延的速度場，決定了林火蔓延的規(guī)則，后者在給定鋒面處速度場的條件下實(shí)現(xiàn)鋒面位置的移動，兩者不斷迭代，實(shí)現(xiàn)不同時(shí)刻鋒面位置的預(yù)測。

1.1 速度場模型

目前已有Rothermel，McArthur，王正非等多種經(jīng)驗(yàn)?zāi)Ｐ?，各有不同的適用場合，但本質(zhì)均是根據(jù)坡度、風(fēng)速、局域可燃物濃度等參數(shù)決定鋒面處蔓延速度。本文側(cè)重于火焰鋒面追蹤，選擇Fendell[13]模型，該模型主要適用于平坦地區(qū)可燃物均勻且較為稀疏的場合，來進(jìn)行風(fēng)速對林火蔓延的影響研究。一般情況下，當(dāng)給定風(fēng)速時(shí)，順風(fēng)火頭處的蔓延較快，逆風(fēng)火尾處的蔓延較慢，使得鋒面形成半圓加拋物線形狀，速度表達(dá)式如式(1)。

(1)

為了便于調(diào)節(jié)參數(shù)，可進(jìn)一步化簡式(1)，得式(2)：

(2)

其中α為火尾處速度(θ= π)與火焰?zhèn)冗?θ= π/2)速度之比。簡化后火頭處(θ= 0)與式(1)一致，但|θ|>π/2時(shí)的速度分布將不受風(fēng)速直接影響，這主要是因?yàn)槠渑c風(fēng)速的關(guān)系難以準(zhǔn)確建模。式(2)并未改變鋒面速度場整體分布情況，但顯著縮小了調(diào)節(jié)范圍，便于根據(jù)實(shí)際測量值進(jìn)行修正。

1.2 水平集法

與傳統(tǒng)的拉格朗日粒子法跟蹤界面不同，水平集法將界面信息隱藏在高維空間中，無需顯式跟蹤界面，因此易于處理鋒面拓?fù)浣Y(jié)構(gòu)變化，基本思想為人工構(gòu)造光滑水平集函數(shù)φ(t,x,y)，使得任意t0時(shí)刻鋒面Γ(t0)處滿足φ(t0,x0,y0) = 0，因此盡管φ(t,x,y)沒有明確的物理意義，其不同時(shí)刻零值的位置仍反映鋒面的變化，為了簡化計(jì)算鋒面法向及曲率等參數(shù)，一般選擇距離函數(shù)作為水平集函數(shù)，即φ(t,x,y)的取值為該點(diǎn)至鋒面Γ(t)的距離，如式(3)：

(3)

(4)

根據(jù)式(4)獲得鋒面法向后，可由幾何關(guān)系進(jìn)一步計(jì)算法向與風(fēng)向之間的夾角θ，代入速度場經(jīng)驗(yàn)?zāi)Ｐ褪?2)后，可得鋒面任意位置的蔓延速度場。

水平集函數(shù)隨時(shí)間的變化可由雙曲微分方程描述：

(5)

其中F為整個(gè)區(qū)域內(nèi)的速度場分布。由于式(2)僅在鋒面處有意義，因此為了將鋒面處速度場推廣至整個(gè)區(qū)域，本文設(shè)置求解區(qū)域內(nèi)每個(gè)點(diǎn)的速度為距離其最近鋒面點(diǎn)處的速度。式(5)描述了給定速度場時(shí)水平集函數(shù)的變化，該方程無需拉格朗日粒子法顯式追蹤鋒面，而是將整個(gè)求解域中水平集函數(shù)進(jìn)行更新，后將值為0的點(diǎn)設(shè)置為新的鋒面。

水平集法通過設(shè)定距離函數(shù)將鋒面信息隱藏在高維空間之中，傳統(tǒng)水平集法求解整個(gè)區(qū)域中水平集函數(shù)的變化，由于鋒面位置連續(xù)變化且水平集函數(shù)為光滑距離函數(shù)，取值不會突變，其余區(qū)域的水平集函數(shù)值并不影響鋒面位置的確定，因此追蹤火焰鋒面只需鋒面附近窄帶中水平集函數(shù)值。同時(shí)使用1.3節(jié)中形態(tài)學(xué)膨脹方法重新初始化距離函數(shù)時(shí)，時(shí)間復(fù)雜度正比于窄帶寬度的平方，使得對整個(gè)區(qū)域重新初始化的計(jì)算量較大，因此本文使用窄帶水平集法，僅計(jì)算鋒面附近而非整個(gè)火焰區(qū)域距離函數(shù)的變化，重新初始化時(shí)可顯著減小計(jì)算量。

1.3 數(shù)值求解

本文在窄帶水平集法的基礎(chǔ)上根據(jù)距離函數(shù)定義使用形態(tài)學(xué)膨脹的方法重建鋒面附近窄帶內(nèi)的距離函數(shù)，避免了使用水平集法重新初始化求解變偏微分方程時(shí)穩(wěn)定態(tài)的問題，改善了速度場模型模擬的實(shí)時(shí)性。數(shù)值求解林火蔓延模型的思路如圖1所示，首先初始化求解域內(nèi)的網(wǎng)格，根據(jù)初始條件設(shè)定水平集函數(shù)，若為圓形或矩形等理想火源，可使用解析法直接計(jì)算距離函數(shù)進(jìn)行設(shè)定，若為火場圖像，則需提取火蔓延鋒面并按照重新初始化水平集函數(shù)中步驟進(jìn)行設(shè)定。重新初始化基于形態(tài)學(xué)中膨脹的概念，過程如圖1所示。

圖1 數(shù)值模擬流程Fig.1 Flow chart of numerical simulation

(1)設(shè)定窄帶寬度W，根據(jù)t0時(shí)刻水平集函數(shù)φ(x0,y0) = 0提取火蔓延鋒面位置(x0,y0)，記錄于數(shù)組P[N]中。由于已將求解域離散化，求解過程中φ(x,y)不會嚴(yán)格等于0，因此本文設(shè)定當(dāng)一個(gè)網(wǎng)格點(diǎn)上下或者左右兩個(gè)網(wǎng)格的水平集函數(shù)發(fā)生異號即認(rèn)為其為鋒面點(diǎn)。求解域內(nèi)其余點(diǎn)φ(x,y)值均賦為Wsgn[φ(x,y)]，其中sgn(φ)為符號函數(shù)，φ>0時(shí)sgn(φ) = 1，φ< 0時(shí)sgn(φ) = -1；

由于本文根據(jù)距離函數(shù)定義重建窄帶內(nèi)水平集函數(shù)，因此重建后可確保窄帶內(nèi)任一點(diǎn)的水平集函數(shù)值為其距鋒面最近點(diǎn)的距離，滿足了距離函數(shù)的要求。同時(shí)重建的時(shí)間復(fù)雜度為O(N)，相比傳統(tǒng)通過求解偏微分方程重新初始化水平集函數(shù)的方式計(jì)算量更小，減小了模擬的時(shí)間，對于林火現(xiàn)場的實(shí)時(shí)建模較為有利。

重建水平集函數(shù)后，根據(jù)式(2)可計(jì)算獲得火蔓延鋒面處的速度場分布，計(jì)算中θ由式(4)描述的鋒面法向與設(shè)定的風(fēng)向組成。如前文所述，水平集法需求解整個(gè)求解域中的水平集函數(shù)，意味著速度場在鋒面以外位置仍然有意義，即需將鋒面處速度延拓至整個(gè)求解域。一個(gè)合理的想法是保證每個(gè)點(diǎn)的速度均與距離其最近的鋒面點(diǎn)速度相同，使得鋒面仍然會按照設(shè)定的速度場移動，由于重建過程中本文已保存了每個(gè)點(diǎn)距離最近的鋒面點(diǎn)位置數(shù)組Px及Py，因此可直接使用這兩個(gè)數(shù)組，無需遍歷比較，減小了計(jì)算量，這也是此重建方法的另一優(yōu)點(diǎn)。

完成當(dāng)前t0時(shí)刻的水平集函數(shù)及速度場計(jì)算后，需求解式(5)獲得下一(t0+t)時(shí)刻的水平集函數(shù)，基本思路為使用有限差分法離散化式(5)。

(6)

將式(6)代入式(5)，可按照nt時(shí)刻水平集函數(shù)及速度場求解(n+1)t時(shí)刻的水平集函數(shù)：

(7)

(8)

為了確保數(shù)值求解的收斂性，在按照圖1所示流程，不斷循環(huán)求解過程中，需要在每一步求解前均對鋒面處每一點(diǎn)驗(yàn)證式(8)，若不滿足要求，則需減小時(shí)間間隔t。當(dāng)達(dá)到指定時(shí)間后，保存數(shù)據(jù)，結(jié)束數(shù)值模擬。

2 結(jié)果分析與討論

本文首先使用理想火蔓延對林火蔓延模型的正確性進(jìn)行了驗(yàn)證。初始火源為圓形，無風(fēng)，火焰均勻向外蔓延。模型的基本參數(shù)如表1所示，其中參數(shù)均為無量綱數(shù)，設(shè)定窄帶的寬度W時(shí)需保證每一時(shí)刻鋒面移動后仍然位于窄帶范圍內(nèi)，本文中設(shè)為0.01。初始時(shí)刻鋒面半徑為0.2，蔓延速度為1，因此每一時(shí)刻的鋒面半徑R= 0.2 +t，與模擬結(jié)果對比可驗(yàn)證模型的正確性。

表1 驗(yàn)證模型正確性的基本參數(shù)

模擬結(jié)果如圖2(a)所示，其中每一時(shí)刻的鋒面均為水平集函數(shù)為0的位置。不同時(shí)刻的鋒面構(gòu)成同心圓，對其使用最小二乘法擬合后獲得t=0, 0.04, 0.08, 0.12, 0.16及0.2時(shí)刻圓半徑分別為0.2000, 0.2399, 0.2790, 0.3181, 0.3573及0.3966，與理論解誤差小于1%，意味著提出模型的模擬結(jié)果可信，由于數(shù)值求解，誤差主要來源于離散時(shí)引入的量化誤差。圖2(b)為t=0.08時(shí)刻水平集函數(shù)重建結(jié)果，可見窄帶的設(shè)定與預(yù)期相符，只在鋒面附近才有意義。我們進(jìn)一步研究了窄帶寬度與計(jì)算時(shí)間的關(guān)系，根據(jù)形態(tài)學(xué)膨脹法重新初始化的原理可知，重新初始化的時(shí)間TN×D2，其中N為鋒面點(diǎn)數(shù)，D為窄帶寬度，因此計(jì)算時(shí)間對于窄帶寬度非常敏感。模擬使用的CPU為Intel Core i5-4200U @ 1.60 GHz，D= 10, 20, 80時(shí)，T= 203 ms, 516 ms, 5797 ms，可見減小窄帶寬度能夠顯著提高計(jì)算速度，但D受到火蔓延速度影響，需確保下一時(shí)刻火焰鋒面仍位于窄帶范圍內(nèi)。

圖2 (a)不同時(shí)刻的鋒面位置，(b) t=0.08時(shí)刻重建窄帶水平集函數(shù)分布Fig.2 (a)Fire fronts at different time moments, (b) Distribution of narrow-band level set function at t=0.08

進(jìn)一步模擬了風(fēng)向?qū)鹇拥挠绊?，使用參?shù)如表2所示，其中U,n,a,α及0均為式(2)中使用。

表2 驗(yàn)證風(fēng)速影響的基本參數(shù)

模擬結(jié)果如圖3所示，風(fēng)向?yàn)樗较蛴遥c預(yù)期結(jié)果相符，在順風(fēng)處蔓延速度較快，逆風(fēng)處蔓延較慢，整體構(gòu)成半圓形+拋物線形結(jié)構(gòu)。

圖3 風(fēng)向水平向右時(shí)不同時(shí)刻的鋒面位置Fig.3 Fire fronts at different time moments with wind along the right direction

當(dāng)求解域存在斜坡時(shí)，蔓延的速度場發(fā)生變化，向上蔓延時(shí)速度增加，向下蔓延時(shí)速度減小，修改后速度場為：

Fslope=F×e2s

(9)

其中s為斜坡的角度，單位為弧度。圖4為求解結(jié)果，其中x[0.65, 0.7]時(shí)s=π/8，x(0.7, 0.75]時(shí)s= -π/8，顯然上坡處蔓延變快而在下坡處減慢。

圖4 坡度影響時(shí)不同時(shí)刻的鋒面位置Fig.4 Fire fronts at different time moments under the influence of slope

圖5 (a)雙火源不同時(shí)刻的鋒面位置，(b) t=0.24時(shí)刻重建窄帶水平集函數(shù)分布Fig.5 (a)Fire fronts at different time moments with two point sources, (b) Distribution of the narrow-band level set function at t= 0.24

水平集法的一個(gè)優(yōu)點(diǎn)是可以方便處理鋒面拓?fù)浣Y(jié)構(gòu)的變化，本文對此進(jìn)行了驗(yàn)證，模型中使用兩個(gè)火源，初始時(shí)半徑均為0.1，圓心分別位于(0.35, 0.45)及(0.65, 0.55)處，其余參數(shù)與表2相同。圖5(a)為模擬結(jié)果，模擬初期兩個(gè)火源為分離狀態(tài)，蔓延互不影響，當(dāng)t=0.16時(shí)，鋒面相交，融合處鋒面逐漸消失，與實(shí)際情況相符，較好地描述了鋒面拓?fù)浣Y(jié)構(gòu)的改變。圖5(b)為t=0.24時(shí)刻的水平集函數(shù)，可見其并未受到鋒面融合影響，仍然滿足距離函數(shù)的定義，進(jìn)一步驗(yàn)證了本文提出的重新初始化方法。

3 結(jié)論

本文針對傳統(tǒng)水平集法重新初始化過程復(fù)雜、計(jì)算量較大的問題，提出根據(jù)定義使用形態(tài)學(xué)膨脹的方法重建鋒面附近窄帶內(nèi)的距離函數(shù)，并基于林火速度場經(jīng)驗(yàn)?zāi)Ｐ蛿?shù)值模擬了單火源及雙火源的鋒面蔓延情況，結(jié)果表明，文中方法可顯著節(jié)約計(jì)算量，有效模擬坡度、風(fēng)速等因素對火焰鋒面的影響，提高了實(shí)時(shí)性，并降低設(shè)備成本，為其在林火蔓延模型中應(yīng)用提供了基礎(chǔ)，在森林火災(zāi)發(fā)生時(shí)，通過與GIS系統(tǒng)聯(lián)用，可方便滅火指揮員現(xiàn)場及時(shí)預(yù)測，為指揮滅火提供依據(jù)。

[1] McArthur AG. Weather and grassland fire behavior[M]. Forestryand Timber Bureau,Department of national Development，Commonwealth of Australia, 1966

[2] Rothermel RC. A mathematical model for predicting fire spread in wildland fuels[R], USDA Forest Service General Technical Report, 1972, 1-115.

[3] Stocks BJ, et al. Overview of the international crown fire modelling experiment[J]. Canadian Journal of Forest Research, 2004, 34(8):1543-1547.

[4] 王正非. 通用森林火險(xiǎn)等級系統(tǒng)[J]. 自然災(zāi)害學(xué)報(bào), 1992, 1(3): 39-44.

[5] Finney MA. FARSITE: Fire area simulator-model development and evaluation[M], United States Department of Agriculture Forest Service- Research Paper RMRS, 1998, 18: RP-4.

[6] Mallet V, et al. Modeling wildland fire propagation with level set methods[J]. Computers and Mathematics with Applications, 2009, 57(7):1089-1101.

[7] Mandel J, et al. Coupled atmosphere-wildland fire modeling with WRF 3.3 and SFIRE 2011[J]. Geoscientific Model Development, 2011, 4(3): 591-610.

[8] Kim M. Reaction diffusion equations and numerical wildland fire models [D]. University of Colorado Denver, Department of Applied Mathematics, Colorado, USA, 2011.

[9] Rehm RG, McDermott RJ. Fire-front propagation using the level set method[M], US National Institute of Standards and Technology, NIST Technical Note 1611, 2009.

[10] Lautenberger C. Wildland fire modeling with an Eulerian level set method and automated calibration[J].Fire Safety Journal, 2013, 62(6): 289-298.

[11] Liu Q, Chan CK. Simulation of a premixed turbulent flame by a conservative level set method[J], International Journal of Computer Mathematics, 2011, 88(10): 2154-2166.

[12] Chopp DL. Computing minimal surfaces via level set curvature flow[J], Journal of Computational Physics, 1993, 106(1): 77-91.

[13] Fendell FE, Wolff MF. Forest fires: behavior and ecological effects[M], San Diego, Academic Press, 2001: 171-223.

Numerical simulation on wildland fire spread with narrow-band level set method

ZHAI Chunjie1， CAO Zhaolou2
(1. Department of Forest Fire Protection，Nanjing Forestpolice College，Nanjing 210023，China；2. Department of Optical Engineering，Nanjing University of Information Science and Technology，Nanjing 210044，China)

Fire front tracking is an important step for modeling wildland fire spread and predicting the position of fire front. Level set method is known for its versatile advantage of handling topological change of front. In this study, we proposed to reconstruct the narrow band near fire front using image inflation based on the definition of sign function to simplify the complicated calculation in classic level set method, which has definite physical meaning. Single-point source and two-point source are simulated. Results show that the method here can effectively simulate the effects of slope and wind speed, which provides a solid foundation for its further use in practical fire spread model.

Fire front tracking; Wildland fire spreading; Level set method; Morphological image processing

2016-09-07；修改日期：2016-10-18

中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金項(xiàng)目(LGYB201615)；江蘇省自然科學(xué)基金(BK20150929)；國家自然科學(xué)基金(61605081)

翟春婕(1988-)，女，漢族，講師，研究方向?yàn)槿紵\斷及林火蔓延機(jī)理。

翟春婕，E-mail：zhaichunjie1988@163.com

1004-5309(2017)-0037-06

10.3969/j.issn.1004-5309.2017.01.05

X932

A