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

翟春婕，曹兆楼(1.南京森林警察学院森林消防系，南京,1003；. 南京信息工程大学光电工程系，南京,10044)

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

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

火焰锋面的追踪是建立林火蔓延模型及预测锋面位置的关键环节。基于水平集法追踪锋面具有求解灵活稳定、易于处理锋面拓扑结构变化的特点。针对传统水平集法重新初始化过程复杂、计算量较大的问题，根据距离函数定义，提出通过形态学膨胀的方法以锋面点为圆心向外膨胀，获得锋面附近窄带内每个点距离锋面的最近距离，重建窄带内的水平集函数，物理意义明确，有效减小了计算量，在此基础上基于林火速度场经验模型模拟了单火源及双火源的锋面蔓延情况。结果表明，该方法能够有效模拟坡度、风速等因素对火焰锋面的影响，为其在林火蔓延模型中的应用提供了基础。

锋面追踪；林火蔓延；水平集；形态学

0 引言

林火蔓延模型通过有效预测林火行为从而帮助防治火灾，提高扑火效率，指导林火管理工作，有巨大的社会经济效益，受到了广泛的关注。由于实际林火过程影响因素众多，如坡度、风向及可燃物含量等，尽管近年来计算机的计算能力增长迅速，对林火动力学理解也愈加深刻，但从基本物理原理出发同时建模求解由大尺度地形决定的气候影响、小尺度地形决定的局部燃烧动力学影响以及瞬态湍流等因素仍非常困难，难以实现林火蔓延的实时预测。因此目前林火蔓延软件一般使用经验模型，通过火蔓延速率与方向总结规律，如澳大利亚的McArthur模型[1]、美国的Rothermel模型[2]、加拿大的国家林火蔓延模型[3]及我国的王正非模型[4]等。

林火蔓延模型在速度场建模之外还包括给定速度场时火焰锋面演进追踪，针对火焰锋面的演进追踪，很多学者提出了有效方法。传统方法一般使用拉格朗日粒子法，比如FARSITE[5]，通过连接相邻的粒子可获得锋面，但该方法计算量与火焰尺寸相关且当锋面拓扑结构发生相交变化时需要特殊处理。相比之下，水平集法通过求解双曲型偏微分方程实现锋面的演进，求解结果经过多年发展具有灵活稳定，且易于处理拓扑结构变化的特点。近年来已有学者将其应用于林火蔓延模型，如Mallet[6], Mandel[7], Kim[8], Rehm[9]等，实现了火蔓延锋面随时间演进的追踪。Lautenberger[10]基于数据融合的思想，通过水平集法追踪锋面成功获得了速度模型中的各个参数。Liu等[11]将其应用于预混火焰中追踪锋面，实现了预期的效果。但水平集法是通过设定距离函数将锋面信息隐藏在高维空间之中，且在追踪过程中需不断重新初始化距离函数以保证其满足定义，运算量较大。由于林火蔓延模型中仅关注火蔓延锋面的位置，其它区域水平集函数的变化并不影响结果，因此有学者提出窄带水平集法，仅计算锋面附近而非整个火焰区域距离函数的变化，显著减小了计算量[12]。虽然这些工作为水平集法更广泛的应用提供了基础，但距离函数的重新初始化仍需求解时变偏微分方程达到稳定状态，时间复杂度较大，影响了模型的实时性。

本文根据目前研究现状，提出结合速度场经验模型及水平集法预测火蔓延锋面位置，针对距离函数重新初始化计算量较大的问题，提出根据距离函数定义使用形态学膨胀的方法重建锋面附近窄带内的距离函数，减小计算量，并数值模拟了单点火源及双点火源在不同外部条件下火焰锋面的蔓延情况，验证了本文提出的方法。

1 基本原理与方法

林火蔓延模型包含两个部分：速度场模型及火蔓延锋面追踪。前者根据锋面位置确定局域的环境如可燃物浓度，风向与锋面法向夹角等进而获得蔓延的速度场，决定了林火蔓延的规则，后者在给定锋面处速度场的条件下实现锋面位置的移动，两者不断迭代，实现不同时刻锋面位置的预测。

1.1 速度场模型

目前已有Rothermel，McArthur，王正非等多种经验模型，各有不同的适用场合，但本质均是根据坡度、风速、局域可燃物浓度等参数决定锋面处蔓延速度。本文侧重于火焰锋面追踪，选择Fendell[13]模型，该模型主要适用于平坦地区可燃物均匀且较为稀疏的场合，来进行风速对林火蔓延的影响研究。一般情况下，当给定风速时，顺风火头处的蔓延较快，逆风火尾处的蔓延较慢，使得锋面形成半圆加抛物线形状，速度表达式如式(1)。

(1)

为了便于调节参数，可进一步化简式(1)，得式(2)：

(2)

其中α为火尾处速度(θ= π)与火焰侧边(θ= π/2)速度之比。简化后火头处(θ= 0)与式(1)一致，但|θ|>π/2时的速度分布将不受风速直接影响，这主要是因为其与风速的关系难以准确建模。式(2)并未改变锋面速度场整体分布情况，但显著缩小了调节范围，便于根据实际测量值进行修正。

1.2 水平集法

与传统的拉格朗日粒子法跟踪界面不同，水平集法将界面信息隐藏在高维空间中，无需显式跟踪界面，因此易于处理锋面拓扑结构变化，基本思想为人工构造光滑水平集函数φ(t,x,y)，使得任意t0时刻锋面Γ(t0)处满足φ(t0,x0,y0) = 0，因此尽管φ(t,x,y)没有明确的物理意义，其不同时刻零值的位置仍反映锋面的变化，为了简化计算锋面法向及曲率等参数，一般选择距离函数作为水平集函数，即φ(t,x,y)的取值为该点至锋面Γ(t)的距离，如式(3)：

(3)

(4)

根据式(4)获得锋面法向后，可由几何关系进一步计算法向与风向之间的夹角θ，代入速度场经验模型式(2)后，可得锋面任意位置的蔓延速度场。

水平集函数随时间的变化可由双曲微分方程描述：

(5)

其中F为整个区域内的速度场分布。由于式(2)仅在锋面处有意义，因此为了将锋面处速度场推广至整个区域，本文设置求解区域内每个点的速度为距离其最近锋面点处的速度。式(5)描述了给定速度场时水平集函数的变化，该方程无需拉格朗日粒子法显式追踪锋面，而是将整个求解域中水平集函数进行更新，后将值为0的点设置为新的锋面。

水平集法通过设定距离函数将锋面信息隐藏在高维空间之中，传统水平集法求解整个区域中水平集函数的变化，由于锋面位置连续变化且水平集函数为光滑距离函数，取值不会突变，其余区域的水平集函数值并不影响锋面位置的确定，因此追踪火焰锋面只需锋面附近窄带中水平集函数值。同时使用1.3节中形态学膨胀方法重新初始化距离函数时，时间复杂度正比于窄带宽度的平方，使得对整个区域重新初始化的计算量较大，因此本文使用窄带水平集法，仅计算锋面附近而非整个火焰区域距离函数的变化，重新初始化时可显著减小计算量。

1.3 数值求解

本文在窄带水平集法的基础上根据距离函数定义使用形态学膨胀的方法重建锋面附近窄带内的距离函数，避免了使用水平集法重新初始化求解变偏微分方程时稳定态的问题，改善了速度场模型模拟的实时性。数值求解林火蔓延模型的思路如图1所示，首先初始化求解域内的网格，根据初始条件设定水平集函数，若为圆形或矩形等理想火源，可使用解析法直接计算距离函数进行设定，若为火场图像，则需提取火蔓延锋面并按照重新初始化水平集函数中步骤进行设定。重新初始化基于形态学中膨胀的概念，过程如图1所示。

图1 数值模拟流程Fig.1 Flow chart of numerical simulation

(1)设定窄带宽度W，根据t0时刻水平集函数φ(x0,y0) = 0提取火蔓延锋面位置(x0,y0)，记录于数组P[N]中。由于已将求解域离散化，求解过程中φ(x,y)不会严格等于0，因此本文设定当一个网格点上下或者左右两个网格的水平集函数发生异号即认为其为锋面点。求解域内其余点φ(x,y)值均赋为Wsgn[φ(x,y)]，其中sgn(φ)为符号函数，φ>0时sgn(φ) = 1，φ< 0时sgn(φ) = -1；

由于本文根据距离函数定义重建窄带内水平集函数，因此重建后可确保窄带内任一点的水平集函数值为其距锋面最近点的距离，满足了距离函数的要求。同时重建的时间复杂度为O(N)，相比传统通过求解偏微分方程重新初始化水平集函数的方式计算量更小，减小了模拟的时间，对于林火现场的实时建模较为有利。

重建水平集函数后，根据式(2)可计算获得火蔓延锋面处的速度场分布，计算中θ由式(4)描述的锋面法向与设定的风向组成。如前文所述，水平集法需求解整个求解域中的水平集函数，意味着速度场在锋面以外位置仍然有意义，即需将锋面处速度延拓至整个求解域。一个合理的想法是保证每个点的速度均与距离其最近的锋面点速度相同，使得锋面仍然会按照设定的速度场移动，由于重建过程中本文已保存了每个点距离最近的锋面点位置数组Px及Py，因此可直接使用这两个数组，无需遍历比较，减小了计算量，这也是此重建方法的另一优点。

完成当前t0时刻的水平集函数及速度场计算后，需求解式(5)获得下一(t0+t)时刻的水平集函数，基本思路为使用有限差分法离散化式(5)。

(6)

将式(6)代入式(5)，可按照nt时刻水平集函数及速度场求解(n+1)t时刻的水平集函数：

(7)

(8)

为了确保数值求解的收敛性，在按照图1所示流程，不断循环求解过程中，需要在每一步求解前均对锋面处每一点验证式(8)，若不满足要求，则需减小时间间隔t。当达到指定时间后，保存数据，结束数值模拟。

2 结果分析与讨论

本文首先使用理想火蔓延对林火蔓延模型的正确性进行了验证。初始火源为圆形，无风，火焰均匀向外蔓延。模型的基本参数如表1所示，其中参数均为无量纲数，设定窄带的宽度W时需保证每一时刻锋面移动后仍然位于窄带范围内，本文中设为0.01。初始时刻锋面半径为0.2，蔓延速度为1，因此每一时刻的锋面半径R= 0.2 +t，与模拟结果对比可验证模型的正确性。

表1 验证模型正确性的基本参数

模拟结果如图2(a)所示，其中每一时刻的锋面均为水平集函数为0的位置。不同时刻的锋面构成同心圆，对其使用最小二乘法拟合后获得t=0, 0.04, 0.08, 0.12, 0.16及0.2时刻圆半径分别为0.2000, 0.2399, 0.2790, 0.3181, 0.3573及0.3966，与理论解误差小于1%，意味着提出模型的模拟结果可信，由于数值求解，误差主要来源于离散时引入的量化误差。图2(b)为t=0.08时刻水平集函数重建结果，可见窄带的设定与预期相符，只在锋面附近才有意义。我们进一步研究了窄带宽度与计算时间的关系，根据形态学膨胀法重新初始化的原理可知，重新初始化的时间TN×D2，其中N为锋面点数，D为窄带宽度，因此计算时间对于窄带宽度非常敏感。模拟使用的CPU为Intel Core i5-4200U @ 1.60 GHz，D= 10, 20, 80时，T= 203 ms, 516 ms, 5797 ms，可见减小窄带宽度能够显著提高计算速度，但D受到火蔓延速度影响，需确保下一时刻火焰锋面仍位于窄带范围内。

图2 (a)不同时刻的锋面位置，(b) t=0.08时刻重建窄带水平集函数分布Fig.2 (a)Fire fronts at different time moments, (b) Distribution of narrow-band level set function at t=0.08

进一步模拟了风向对火蔓延的影响，使用参数如表2所示，其中U,n,a,α及0均为式(2)中使用。

表2 验证风速影响的基本参数

模拟结果如图3所示，风向为水平向右，与预期结果相符，在顺风处蔓延速度较快，逆风处蔓延较慢，整体构成半圆形+抛物线形结构。

图3 风向水平向右时不同时刻的锋面位置Fig.3 Fire fronts at different time moments with wind along the right direction

当求解域存在斜坡时，蔓延的速度场发生变化，向上蔓延时速度增加，向下蔓延时速度减小，修改后速度场为：

Fslope=F×e2s

(9)

其中s为斜坡的角度，单位为弧度。图4为求解结果，其中x[0.65, 0.7]时s=π/8，x(0.7, 0.75]时s= -π/8，显然上坡处蔓延变快而在下坡处减慢。

图4 坡度影响时不同时刻的锋面位置Fig.4 Fire fronts at different time moments under the influence of slope

图5 (a)双火源不同时刻的锋面位置，(b) t=0.24时刻重建窄带水平集函数分布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

水平集法的一个优点是可以方便处理锋面拓扑结构的变化，本文对此进行了验证，模型中使用两个火源，初始时半径均为0.1，圆心分别位于(0.35, 0.45)及(0.65, 0.55)处，其余参数与表2相同。图5(a)为模拟结果，模拟初期两个火源为分离状态，蔓延互不影响，当t=0.16时，锋面相交，融合处锋面逐渐消失，与实际情况相符，较好地描述了锋面拓扑结构的改变。图5(b)为t=0.24时刻的水平集函数，可见其并未受到锋面融合影响，仍然满足距离函数的定义，进一步验证了本文提出的重新初始化方法。

3 结论

本文针对传统水平集法重新初始化过程复杂、计算量较大的问题，提出根据定义使用形态学膨胀的方法重建锋面附近窄带内的距离函数，并基于林火速度场经验模型数值模拟了单火源及双火源的锋面蔓延情况，结果表明，文中方法可显著节约计算量，有效模拟坡度、风速等因素对火焰锋面的影响，提高了实时性，并降低设备成本，为其在林火蔓延模型中应用提供了基础，在森林火灾发生时，通过与GIS系统联用，可方便灭火指挥员现场及时预测，为指挥灭火提供依据。

[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] 王正非. 通用森林火险等级系统[J]. 自然灾害学报, 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

中央高校基本科研业务费专项资金项目(LGYB201615)；江苏省自然科学基金(BK20150929)；国家自然科学基金(61605081)

翟春婕(1988-)，女，汉族，讲师，研究方向为燃烧诊断及林火蔓延机理。

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

1004-5309(2017)-0037-06

10.3969/j.issn.1004-5309.2017.01.05

X932

A