大兴安岭呼中林区林火蔓延三维模拟

柳生吉 (铁道第三勘察设计院集团有限公司，天津 300251)

1 研究综述

森林火灾是一种突发性强、破坏性大的自然灾害。林火烧毁森林及地表物，破坏了森林生态系统，造成生态系统内生物因子、生态因子的混乱，大大降低了森林保持水土、涵养水源、调节气候的作用。森林防火工作事关森林资源和生态安全，研究林火蔓延规律和模拟林火蔓延现象有助于森林防火。

林火燃烧是指森林中的可燃物在一定温度下快速与氧气结合并发光发热的化学物理反应［1］。林火发生后，会沿着不同的方向蔓延。林火蔓延是一种林火行为，是森林可燃物从点燃开始直至熄灭的整个过程中所表现出的特性。林火专家将林火蔓延的影响因子主要归结为可燃物、地形和气象3个方面［2］。林火蔓延可视化是根据研究区的实际情况选取合适的林火蔓延模型，采用高效的计算机模拟算法和可视化技术展现林火发生发展的过程。

目前常见的林火蔓延算法按照对空间要素的处理方法可分为基于矢量数据形式和栅格数据形式2类算法［3］。前者即为惠更斯原理，后者中较流行的有边界插值算法、迷宫算法以及自动元胞机算法等。矢量或惠更斯原理的特点是林火蔓延的过火区域是一个随时间变化的连续扩展的多边形。该多边形的形状由一系列过火区边缘上的二维顶点决定［4］。为了保证一定的精度要求，随着多边形的增大，二维顶点数会不断增加。在林火蔓延的过程中，每个顶点被认为是一个独立的火点，其蔓延形状也被认为是椭圆形，蔓延方向由风速矢量和坡度矢量叠加决定，蔓延速率则由风速、坡向、可燃物等因子通过计算林火蔓延速率模型得到。这种算法精度较高，但是计算效率较低。目前应用较多的是模拟精度略低，而计算效率较高的栅格形式算法，如美国农业部开发的FARSITE林火管理系统［5－6］。国内学者对林火蔓延可视化的研究也集中于栅格形式的算法［7］。刘月文等［8］、陈喆等［9］采用王正非模型，分别实现了基于二维和三维元胞自动机模型的林火蔓延模拟，而黄华国等则采用Rothermel模型实现了基于三维曲面元胞自动机模型的林火蔓延模拟［10］。王海军等在元胞自动机林火蔓延模型的基础上提出元胞作用域的概念，综合了元胞自身属性和邻居元胞属性对林火蔓延的影响［11］。杨广斌等在动态数据驱动应用系统(DDDAS)技术方面开展的研究在北京市森林防火系统中得到了很好的应用，它体现了一种基于林火蔓延模型适宜性自动选择的林火蔓延模拟的全新结构体系和技术框架［12］。DDDAS是一种功能强大的仿真应用的全新模式，它将是林火蔓延模拟发展的一个重要趋势［13］。

随着多维可视化技术及虚拟现实技术的发展，林火蔓延三维模拟及森林景观虚拟有所发展。虚拟的林火蔓延包括三维环境林火造型模拟、地形上林火蔓延过程模拟及林火与其他景观交互模拟等方面内容［14］。虽然这些技术在模拟林火蔓延时还存在诸多不足，但代表了现代模拟林火蔓延的发展方向。三维可视化开发分为基于底层3DAPI开发和基于二次开发平台的开发。底层3DAPI如DirectX、OpenScene-Graph等都是免费的，开发出的程序也可以脱离其运行环境，就是开发的过程比较复杂。Erdas公司的Virtual GIS，ESRI公司的3D Analyst都是较典型的二次开发三维可视化平台，它们都脱离不了其运行环境，并且价格较昂贵。目前常见的林火造型模拟技术有纹理合成、粒子系统［15］以及基于物理的模拟等3种［16］。粒子系统是应用最广泛的三维林火造型技术，它具有实时性好、复杂度低、真实感强、适合于大范围等优点。大兴安岭呼中林区林火发生频繁，对当地的林业资源造成严重破坏。笔者在遵循林火蔓延规律的基础上开发出一款林火蔓延三维模拟软件，并针对呼中林区进行初步验证。

2 数据来源与研究方法

2．1 研究区概况 呼中林区位于黑龙江省北部，其位置如图1所示。

图1 呼中林区地理位置

呼中区地理坐标为 122°39'30″～124°21'00″E、51°14'40″～52°25'00″N，总面积约741 999 hm2。呼中春秋2 季受蒙古干旱风影响，天气条件变化剧烈，常出现高温、低湿和大风天气，因而春季和秋季是林火的高发期。呼中地貌属石质中低山地，坡度平缓，一般小于15°。呼中植被在植物区系上属泛北极植物区东西伯利亚植物区系，以西伯利亚植物区系成分为主，混有东北植物区系和蒙古植物区系成分。地带性植被类型为寒温性针叶林，以兴安落叶松(Larix gmelinii)为单优势种。包括5种植被类型:以兴安落叶松为主或混生有樟子松(Pinus sylvestris var．mongolica)、白桦(Betula ptyphylla)的兴安落叶松林;樟子松林或伴生兴安落叶松、白桦、山杨(Populus davidiana)的樟子松林;伴生兴安落叶松和白桦的红皮云杉(Picea koratensis)林;沿河甜杨林和沿河钻天柳(Chosenia arbutifolia)林;白桦林和山杨林。另外，在海拔高于800 m的地方分布有可燃性很好的林下灌木偃松。

2．2 数据来源 该研究中用到的DEM是SRTM的分辨率为30 m的数据，林区地形可视化选择Landsat5的TM影像作纹理。可燃物分类数据利用呼中林区2000年林相图获取。林火发生当日的气象数据从相关部门的气象观测数据获得。

2．3 研究方法

2．3．1 林区地形三维可视化。OSG中三维场景组织是通过场景树方式实现［17］，场景树的基本单位是节点，只有一个根节点(Root)，根节点下面包括组节点(Group)和叶节点(Node)。叶节点中管理一个或者多个可绘制对象(OSG固有的和用户自定义的)的信息，并可以通过其接口函数对可绘制体的信息进行查询［18］。该研究利用海量地形数据处理工具VPB(Virtual PlanetBuilder)对获取整个呼中林区的高程数据和影像数据格式化处理，将转换后的地形数据作为OSG所组织的三维场景的一个叶子节点进行三维显示。3DVFS中的OSG场景树结构如图2所示。

2．3．2 林区坡度分类。利用ArcGIS软件Spatial Analysis模块的Slope工具处理DEM数据，获取研究区用于林火蔓延的坡度数据。

图23DVFS中的OSG场景树结构

2．3．3 林区可燃物分类。根据吴志伟等的研究［19］，应用ArcGIS将呼中林区的可燃物主要划分为以下4类:偃松灌丛;以白桦和杨树为主的阔叶林;分布在阴坡，以杜香和越橘为主要林下灌木的针叶林;分布在阳坡，以兴安杜鹃为主要林下灌木的针叶林。在BehavePlus［20］软件中输入吴志伟等的研究结果中的相关参数，计算出不同可燃物模型的地表火蔓延速率作为3DVFS的林火蔓延初始化速率。

2．3．4 3DVFS 设计

2．3．4．1 林火蔓延模型。Rothermel林火蔓延模型以能量守恒定律等物理学理论为基础，以林火实验为依据，是一个半经验的模型［21－22］，适用范围较为广泛，但是该模型需要的输入森林可燃物相关的参数较多。Rothermel模型的林火蔓延公式如下:

式中，R是蔓延速度(m/min);IR是反应强度(kJ/min·m2);ζ是林火蔓延率;φW是风速修正系数;φS是坡度修正系数;ρb这是可燃物体积密度(kg/m3);ε为有效热参数(无因次量);Qt指点燃单位质量可燃物所需要的热量(kJ/kg)。

呼中林区是大兴安岭地区林火发生比较多的地方高发区，近些年来大量林火研究者皆在呼中展开研究，积累了大量关于该林区可燃物方面的珍贵资料。这些资料可以引用到美国Rothermel模型中做林火蔓延模拟相关的研究［23］。因此，该研究选择Rothermel模型做为林火蔓延模型。

2．3．4．2 林火蔓延算法。该研究应用的算法是Mark A．Finney介绍的一种方法:minimum travel time(最短传播时间算法)［24］。最短传播时间算法是在一个与地形、可燃物等具有相同空间位置、相同空间分辨率的二维矩形栅格中进行计算。最短传播时间指的是从火源栅格单元点直线传播到与它相邻的栅格单元点所需要的最短时间。

图3 计算距着火点dY和dX处林火传播最短时间示意

如图3所示，对于最短传播时间算法，林火行为特征根据主林火传播方向、最大传播速度和强度，以及椭圆形的火形状尺寸等计算得到。火形状参数和最大传播速度决定了一场平行于地面的椭圆形状的林火在笛卡尔坐标系下的传播速度［25－26］:

式中，θ是指定栅格点和椭圆中心的连线与最大传播速度方向的夹角，θ∈［－π，π］;b+c是向前最大传播速度(θ=0);a是火翼传播速度，着火点不一定是在上述椭圆的后交点处。当横切线段(dX，dY)距离给定之后，就可以用公式(2)、(3)计算出θ值:

式中，α是最大传播速度的方向;β是横切线段与最大传播速度方向的夹角。根据以上公式，推算出任意方向上垂直于火线边缘的火传播速度计算公式如下:

上面提到计算方法是基于二维栅格，而林火蔓延椭圆计算是平行于地表面的。因此，当地形具有坡度时，必须用横切线段长度除以其移动方向坡度的余弦值来计算对应的林火传播时间。

林火蔓延计算FireEngine模块由4个自定义类构成:FireCell、ParameterSet、FireStart以及 FireSpread。FireCell类派生自Icomparable接口，其主要属性有行列号、林火到该栅格的最小传播时间，主要方法有SetValue和CompareTo，分别提供栅格单元属性的赋值功能和与其他栅格相互比较最短传播时间的功能。ParameterSet提供在系统运行过程中，存储模拟参数的功能，它包含了所有的模拟参数。FireStart类主要包含 CreateWorkSpace、InitilizeFireMap、Readfuelmap、Read-SlopeMap、SetParamSpace等方法，分别提供创建工作空间、给模拟所需数组变量分配内存、初始化林火蔓延状态栅格图、从坡度文件、可燃物文件读取相应数据，将保存在电脑的模拟参数列表文件中的参数值赋给ParameterSet。FireSpread类是该模块的核心算法部分。该研究中用来计算林火传播的栅格有3种状态:未燃烧、点燃、熄灭，分别用数字0、1、2表示。随着时间的增加，林火传播栅格的状态将按照最短传播时间算法的规则不断转换。如图4所示，该研究中最短传播时间算法的基本步骤分为:

(1)将所有的栅格单元的最短到达时间赋值为∞。

(2)初始化起火点并计算林火传播到紧邻起火点的8个栅格单元所需要的时间Ti。

(3)找出所用传播时间最短的栅格单元，将其状态修改为“燃烧”。

(4)计算林火从起火点传播到新引燃栅格单元再传播到与新引燃点紧邻的8个栅格单元(处于点燃或熄灭状态的除外)所需的累积时间Ti。

(5)将步骤(2)算得的时间和步骤(4)算得的时间进行比较，找出传播用时最短的栅格单元，将其状态赋为“燃烧”，而将步骤(3)中的栅格单元的状态赋为“熄灭”。

(6)从步骤(4)开始不断地重复计算，直至该过程满足模拟控制条件，运算结束。

该研究将模拟过程中不同时间段产生的林火蔓延状态栅格存储到NetCDF(Network Common Data Form)网络通用数据格式文件中。

2．3．4．3 开发平台。基于3DVFS系统所预期的运行操作系统以及具体的可视化编程方面的要求，该研究选择微软推出的功能强大的开发环境Microsoft Visual Studio 2010作为编程工具，选择开源的、跨平台的OSG(OpenSceneGraph)作为三维渲染平台，从底层构建该模拟系统。

3 结果与分析

3．1 3DVFS功能实现 3DVFS功能设计采用C#编程实现，界面主要划分为菜单和工具栏操作区、模拟参数输入和模拟结果输出区、三维场景可视化显示区3部分。菜单和工具栏操作区提供一些用户常用的三维视图操作命令，如三维可视化模型数据的导入，三维场景蔓延控制、场景漫游速度设置以及退出系统等。模拟参数输入和模拟结果输出区根据模型具体的运行需求，提供林火模拟输入参数设置、启动停止模拟功能以及模拟结束后相关统计结果的输出。三维场景可视化显示区是一个展示林火蔓延的动态过程的视口，它是一个自定义的封装OSG相关功能的控件。3DVFS系统界面如图5所示。

图5 3DVFS界面

图6 林火蔓延模拟各时间段火场边界变化

3．2 呼中林区林火模拟 该研究选择2010年6月10日发生在呼中国家级自然保护区的一场林火进行模拟。这场火的起火点位置坐标为x=493 951 m，y=574 839 0 m，实测过火面积约为580．14 hm2。模拟过程中，设置风向为45°方向(西南)，风速为12．5m/s，模拟总步长为106步，并加载可燃物类型及坡度等其他参数。模拟过程中各个时间阶段火场边界变化如图6所示。模拟所得的火场形状近似为椭圆，火场边缘呈现锯齿形状，标示为黄色的栅格单元是着火点，粒子系统火焰表示正在燃烧的栅格单元，标示为黑色的栅格单元是熄灭点。

模拟结束后，对模拟火场面积精度进行评价。该研究选择的评价指标为周宇飞等的计算方法［27］:

式中，SD是模拟和实测过火面积不重叠的面积，SR是实测过火面积，是相对误差。将本实验中实测过火面积和模拟过火面积叠合(图7)，经计算ε为40．46%，二者叠合程度为中等。

图7 实测面积和模拟面积的叠合效果

4 结论与讨论

(1林火燃烧是一个及其复杂的现象，使得林火蔓延模拟十分困难。林火蔓延模型是对林火蔓延现象的数学或物理上的近似，模型的简化效果使得要非常精确的预测林火行为基本上不可能。该研究只是初步完成了一个林火蔓延模拟软件的开发，对于该软件的模拟结果，只是比较了实测过火面积和模拟过火面积的叠合程度，暂时无法对其预测精度进行十分准确地评价。

(2)在考虑地形坡度对林火蔓延的影响时，为了提高运算效率，研究中简化了一些几何方面的计算，势必会对模拟结果的精度产生影响。但是，目前还没有找到更为有效的解决方法。

(3)该研究中可燃物类型是参考一些研究人员在实验林区的相关研究结果来粗略划分的。实际上可燃物的状况又是随时间发生变化的，这里采用的是一种对可燃物状况相对静态的描述。所以，可燃物类型划分对林火蔓延模拟结果影响的不确定性是比较大的，今后要探索更为精确的可燃物类型划分方法。

(4)对于该研究开发的3DVFS系统，在整个一场林火蔓延模拟的过程中，假设风速以及风向是不发生变化的，但实际中风速和风向是因时因地发生变化的，这是将来需要改进的地方。当然，目前在采集实时变化的风速和风向数据方面也存在一些问题，需要在今后的研究中加以重视。

(5)该研究采用了OpenSceneGraph的粒子系统技术来实现对真实火焰的模拟［28］。在一个比较大的尺度上模拟林火蔓延时采用粒子系统将会面临模拟计算用时过长对粒子系统动态更新的影响。针对这个问题，该研究采用将不同模拟时段林火蔓延的状态存储在自定义的NetCDF文件中，等模拟计算结束之后再从NetCDF读入到计算机内存中用于动态显示，取得了较好的效果。

(6)目前比较常见的林火蔓延模拟很多是基于OpenGL或者像ArcGIS等一些平台，而该研究采用的是OpenScene-Graph三维渲染引擎。实践证明，采用OpenSceneGraph开发具有很大的优势，相比OpenGL开发复杂度明显减小，工作量减少。该研究后期的目标是完成3DVFS的模拟精度验证以及将SpeedTree树木建模软件［29］和OpenSceneGraph相结合，进而完成呼中林区整个森林景观的虚拟可视化，为林区森林资源管理提供服务。

［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］．世界林业研究，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］唐晓燕，孟宪宇，葛宏立，等．基于栅格结构的林火蔓延模拟研究及其实现［J］．北京林业大学学报，2003，25(1):53 －57．

［8］刘月文，杨宏业，王硕，等．一种基于CA的林火蔓延模型的设计与实现 － －以内蒙古地区为例［J］．灾害学，2009，24(3):98 －102．

［9］陈喆，孙涛，张凌寒，等．三维元胞自动机各向异性林火蔓延快速模型［J］．北京林业大学学报，2012，34(1):86 －91．

［10］黄华国，张晓丽，王蕾．基于三维曲面元胞自动机模型的林火蔓延模拟［J］．北京林业大学学报，2005，27(3):94 －97．

［11］王海军，张文婷，陈莹莹，等．利用元胞自动机作用域构建林火蔓延模型［J］．武汉大学学报:信息科学版，2011，36(5):575 －578．

［12］杨广斌，刘鹏举，唐小明．动态数据驱动的林火蔓延模型适宜性选择［J］．林业科学，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］．林业调查规划，2012，37(2):53 －58．

［15］赵春霞，张艳，战守义．基于粒子系统方法的三维火焰模拟［J］．计算机工程与应用，2004，28(3):35 －37．

［16］国志兴，钟兴春，方伟华，等．野火蔓延灾害风险评估研究进展［J］．地理科学进展，2010，29(7):778 －788．

［17］王锐，钱学雷．OpenSceneGraph三维渲染引擎设计与实践［M］．北京:清华大学出版社，2009．

［18］许列，韦群，王珏．基于OSG的三维场景管理及实时绘制技术研究与实现［J］．装备指挥技术学院学报，2011，22(3):100 －104．

［19］吴志伟，贺红士，胡远满，等．FARSITE火行为模型的原理，结构及其应用［J］．生态学杂志，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，Forest Service，US Department of Agriculture，1972．

［22］单延龙，舒立福，王洪伟，等．Rothermel火蔓延模型特征参数的解析［J］．森林防火，2003(1):22 －25．

［23］单彦龙，李华，其其格．黑龙江大兴安岭主要树种燃烧性及理化性质的实验分析［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］．北京林业大学学报，2010，32(2):21 －26．

［28］辛佳佳．基于OSG的粒子特效仿真研究［J］．计算机光盘软件与应用，2011(13):65－66．

［29］刘颖，罗岱，黄心渊．基于OSG的Speedtree植物模型绘制研究［J］．计算机工程与设计，2012，33(6):2406 －2409．