集合卡曼滤波算法对FARSITE模型林火蔓延预测的修正效果研究

钱 兰，张启兴，张永明

(中国科学技术大学火灾科学国家重点实验室，合肥，230026)

0 引言

森林火灾是一种全球范围内普遍存在的灾害现象，具有突发性和不确定性。火灾不仅会破坏生态系统、污染大气环境、还有可能对人类生命和财产安全造成威胁。近年来，伴随着气候变暖，火灾形势更加严峻，森林火灾频繁发生[1]。

利用林火数值模型进行计算机仿真模拟是指导火灾预防和扑救、辅助林火管理的有效手段。但计算机模拟也存在一定的局限性，模型构建过程中对复杂火灾演化机理的不完全理解以及输入参数的缺失、不准确，都会降低模型预测结果的可信度。动态数据驱动的仿真模拟技术能在传统模拟方法的基础上融入观测数据对模型预测轨迹进行动态调整，从而有效提高模型预测精度，近年来这一新方法被越来越多地应用到森林火灾领域。

杨广斌[2]提出了基于人工神经网络的模型适宜性选择技术和林火模拟误差修正参数自动生成技术；周宇飞[3]提出了利用真实火场数据和专家干预来调整模型参数的方法；Rochoux等[4]验证了扩展卡曼滤波(EKF)算法对Rothermel模型参数进行修正的可行性；一些研究利用遗传算法(GA)对林火蔓延预测模型中的输入参数进行组合和迭代，搜索输出与观测数据匹配度最高的近似最优解[5-7]；粒子滤波(PF)算法也被广泛应用于林火模型的状态修正[8-12]。近年来，基于集合卡曼滤波(EnKF)算法的森林火灾蔓延预测修正方法被提出并得到了发展，被用于林火模型中可燃物、气象、火线位置等参数的修正[13-16]。

以上方法在实际应用过程中：遗传算法中问题的编码过程比较复杂；PF算法需要的粒子数众多，计算量过大；EKF算法则要对系统做线性近似处理，而基于集合思想的EnKF算法不仅可以很好地适用于林火这类非线性系统，还适合开展并行计算，运算效率高。EnKF算法对林火模型的修正研究方面：与文献[13,14]使用FIREFLY模拟简单工况下的小尺度(百米级别)林火不同，本文选用FARSITE作为林火蔓延预测模型，模拟了多可燃物、不均匀地形、变化气候条件下的大尺度(千米级别)森林地表火灾，更加接近真实火灾场景；同时区别于文献[15,16]，本文采用了不同的火灾案例和火源集合生成方式，并综合分析了集合元素个数、观测数据标准差、同化频率对修正效果的影响。研究可为算法的参数分析和同化方案设计提供借鉴，并为林火精准扑救、火灾应急决策提供理论参考。

1 FARSITE模型

FARSITE模型是一个二维矢量模拟系统，由美国林务局(US-FS)Missoula火灾科学实验室开发。它能够模拟不同环境条件下的地表火、树冠火、飞火等多种类型火灾的蔓延过程，接近真实林火场景，实用性较强。由于地表火是最常见的一种林火类型，因此本文模拟的也是地表火。Rothermel模型和Huygens原理是FARSITE地表火子模型的两大重要理论基础[17]：FARSITE将地表火的蔓延视为位于燃烧区和未燃区之间的极薄火焰区(称为火线)垂直于自身向前扩散的过程，火线的形状由其上一系列离散标记点决定，Rothermel模型则用于求解这些标记点在风和坡度矢量合成方向上的最大稳定传播速率[18]。Rothermel模型是一个半经验模型，模型的速度求解方程是依据能量守恒定律建立的，方程中的一些参数需要通过试验来获取。模型从宏观尺度描述火灾蔓延，对不同的传热模式不做区分，因此对一些复杂的燃烧流动现象不做考虑。FARSITE借助Huygens波动理论将火灾的空间增长模拟为椭圆波的扩散，在单位时间步长内，火线上的每个标记点都被看作独立火源并以椭圆形式向前蔓延，如图1。椭圆的形状(即a/b，a/c的值)和方向(即火头前进方向)由风矢量和坡矢量合成矢量决定，而椭圆大小则由Rothermel模型根据由风速、坡向、可燃物等因子计算得到的火头蔓延速度和模拟时间步长决定。所有标记点在前一时刻火线基础上扩散形成的波动椭圆所构成的包络面就是新时刻的火线位置。其中周线分辨率用于定义火线上相邻标记点间的距离，距离分辨率是标记点的最大投影传播距离，在该距离内环境参数保持不变，原理如图2[19]。

图1 火线上标记点的椭圆蔓延示意图Fig.1 Diagram of the elliptical spread of a marker on the fire perimeter

图2 FARSITE模型中的惠更斯原理说明图Fig.2 Illustration of the Huygens' principle in FARSITE simulator

FARSITE能综合考虑地形(海拔，坡度，坡向)、可燃物、气候(风况，降水，湿度)等因素对林火蔓延的影响，但不对火灾的复杂物理化学过程做精细求解，因此可以对大空间尺度的森林火灾进行快速模拟预测，能满足实时预报的要求，同时可以二维矢量图的方式直观展现林火行为状态，在林火管理中具有重要应用价值。

为满足编程需要，本研究使用的是FARSITE的Linux1.0版本，可以通过在终端输入TestFARSITE+路径参数的方式来实现模型调用。模型运行必需的输入文件有：景观要素文件、模拟参数文件、火源位置文件和路径文件，可选的为障碍物位置文件。火源位置和障碍物位置文件是shape格式。景观要素文件的扩展名为.LCP，是FARSITE特有的文件格式，该文件至少包含5类主题数据：海拔、坡度、坡向、可燃物类型和树冠覆盖度。模拟参数文件则用于存储自定义可燃物的特征参数、可燃物湿度、气象等数据，并定义一些全局性变量：如模拟开始和结束时间、模拟时间步长、距离分辨率、周线分辨率等。路径文件记录了所有输入文件的绝对路径，路径文件和模拟参数文件为文本格式。Linux版本的FARSITE共有26个输出文件，这些文件记录了模拟过程中每个时间步长下的火线位置、火焰长度、火线强度、火焰蔓延速度等关键林火行为数据。

2 FARSITE-EnKF火线预测动态修正方法

FARSITE是一个确定性模型，它以一组特定参数为输入，得到的输出结果也是一定的。实际模拟过程中，FARSITE输入参数众多，难免存在数据缺失和误差，Rothermel速度计算模型的应用也具有一定局限性(如要求地形、可燃物空间分布连续，风等动态环境参数不能变化太快[20])，这些都会在一定程度上影响预测结果的准确性。而EnKF算法能在模型预测过程中同化反映系统真实状态的观测数据，有效提高模型预测能力。EnKF[21]由海洋学者Geir Evensen于1994年提出并应用推广到同化领域，它是一种顺序数据同化算法，主要应用于非线性系统。EnKF算法的状态估计主要包括预测和更新两个过程：首先根据k时刻状态参量的统计特征，产生一组蒙特卡罗集合用于表征状态误差，然后利用模型向前积分直到有新的观测值输入，得到k+1时刻的状态预测集合，完成预测。再依据k+1时刻观测数据的统计特征，同样地产生一组观测集合，将所有集合元素的状态预测值和观测值进行相应的加权得到分析值，加权系数则根据各自误差大小确定。所有集合元素分析值的平均就是状态变量的最终分析值，集合的统计误差也即状态变量的分析误差协方差，更新结束。重复上述过程，直到完成所有观测数据的同化。

为提高FARSITE模型的预测准确性，本文提出了FARSITE-EnKF动态修正方法。考虑到火线位置相比于其他林火行为指标对火灾管控和扑救更具有直观指导意义，因此选取火线位置作为待修正的模型状态参量，并忽略模型本身及其他输入参数的误差，FARSITE-EnKF方法流程图如图3所示。

图3 FARSITE-EnKF方法流程图Fig.3 Flowchart of the proposed FARSITE-EnKF method

结合上文思路，可以将FARSITE模型抽象为一个非线性离散系统，系统的状态方程和观测方程可以写成：

Xk+1=M(Xk,θk)

(1)

(2)

状态方程(1)中X为状态参量，X可以表示成一个由火线上标记点坐标组成的2m维向量X=(e1,n1,e2,n2…ej,nj…em,nm)T(m为火线上标记点总个数，ej,nj为火线上第j个点的坐标)。Xk可以理解为k时刻的火线位置输入，θk是该时刻模拟所需的其他输入参量，如地形、可燃物、气象数据等。M是一个非线性状态转换方程，实现从k时刻火线位置Xk到k+1时刻火线位置Xk+1的状态预测，对应于FARSITE的正向模拟过程。观测方程(2)中Yo是观测量，Yo∈R2s,(s

综合EnKF计算公式[22]，FARSITE-EnKF方法的计算步骤如下：

(3)

Xi,k=(ei,1,ni,1,ei,2,ni,2…ei,j,ni,j…ei,mk,ni,mk)T

(4)

(2)利用FARSITE进行集合预测。

(5)

集合预测结果归为模拟组。

(6)

(7)

(4)产生观测集合。

(8)

(9)

(10)

(5)计算卡曼增益矩阵。

(11)

(12)

(13)

(14)

(15)

(16)

对每个集合元素计算分析值，所有集合元素分析值的平均即状态参量的最终分析值，也就是火线预测的最终修正结果，此时的结果归为EnKF组。

3 算例设计

本文借助了数据同化研究中广泛使用的观测系统模拟试验方法(OSSE)[23]：首先做一组无参数误差的模拟(真实组)，将模拟结果作为真实值记录下来，“观测数据”从中产生。另外做一组有参数误差的模拟(模拟组)，产生一组不准确的预测值。然后用EnKF算法同化“观测数据”，对模拟组的模型预测结果做修正，得到最终分析值(EnKF组)。

由于模拟所需输入参数众多，实际获取困难，而本文的主要研究目的是从理论上验证修正方法的可行性，因此同Lin和Liu[24]一样，本文选用了FARSITE自带的Ashley案例进行算法验证。Ashley案例中模拟区域在东西方向的跨度为8 km，南北方向跨度为11 km。火场环境因素复杂：地形随空间位置发生变化、可燃物成分复杂且非均匀分布、气候条件也是时变的，接近真实火灾场景。图4用不同颜色给出了Ashley案例的可燃物分布并标出了本实验中真实组和误差组的火源位置。其他参数设置见表1。

图4 Ashley案例的可燃物分布及两组火源位置Fig.4 Fuel distribution and the locations of two designed fire sources in Ashley case

表1 基本输入参数设置Table 1 Arrangement of the basic input parameters

4 结果分析与讨论

4.1 EnKF算法对FARSITE火线预测的修正效果

我们输入Ashely的基础数据和表1中的自定义参数进行模拟，验证FARSITE-EnKF方法的可行性。图5给出了真实火源、猜测火源和stdc=100 m时对火源猜测值施加扰动所生成的20个集合元素的分布情况。

图5 真实火源、猜测火源及stdc=100 m时的火源扰动集合Fig.5 Real,guessed fire source and initial fire source ensemble with a centroid standard deviation of 100 m

图6、图7分别为T=9 h—T=18 h真实火源和猜测火源对应的火线预测结果。

图6 基于真实火源的火线预测Fig.6 Fire perimeters prediction based on real fire source

图7 基于猜测火源的火线预测Fig.7 Fire perimeters prediction based on guessed fire source

可以看到，由于火源的输入偏差，FARSITE给出了不同的火线输出结果。基于火源猜测值模拟得到的火线形成的过火面积相比于同一时刻的真实面积明显偏大。两组结果的预测偏差可以用同化算法中常用的均方根误差(RMSE)来衡量，RMSE值越小，偏差越小。RMSE的计算公式为：

(17)

图8为根据公式17计算得到的模拟组的火线预测偏差随时间的变化图。从图8中可以看出：未使用EnKF算法同化观测数据时，从T=9 h开始RMSE呈现快速增长趋势，虽然在T=16 h～18 h时段内有小幅度减小，但整体仍维持在一个较高的数值水平。

图8 真实组和模拟组预测火线间的RMSEFig.8 RMSE of fire perimeters between real prediction and the error one

相反地，如果从T=10 h开始每小时输入一次来自真实组的观测数据，并用EnKF算法同化观测数据对预测结果进行修正，模型预测结果就能够得到很好改善。由于篇幅有限，本文仅选取了T=11 h、T=13 h、T=15 h和T=17 h四个时刻修正后的效果图作为展示，如图9～图12所示。

图9 T=11 h时真实组、模拟组和EnKF组的火线预测结果比较Fig.9 Comparison of fire perimeter prediction at T=11 h from real group,simulated group and the group using EnKF algorithm

图10 T=13 h时真实组、模拟组和EnKF组的火线预测结果比较Fig.10 Comparison of fire perimeter prediction at T=13 h from real group,simulated group and the group using EnKF algorithm

图11 T=15 h时真实组、模拟组和EnKF组的火线预测结果比较Fig.11 Comparison of fire perimeter prediction at T=15 h from real group,simulated group and the group using EnKF algorithm

图12 T=17 h时真实组、模拟组和EnKF组的火线预测结果比较Fig.12 Comparison of fire perimeter prediction at T=17 h from real group,simulated group and the group using EnKF algorithm

结果表明，同化观测数据后，所有时刻的火线预测偏差都得到了很好的修正。修正后的火线与真实火线不仅在空间位置上十分接近，而且在火线形状上也呈现出高度一致性。为了定性分析算法的修正效果，本文计算了EnKF算法逐时同化后的RMSE值，如图13。

图13 真实组和EnKF组火线间的RMSEFig.13 RMSE of fire perimeters between real group and the group using EnKF algorithm

图13中RMSE值代表火线位置这一状态参量的修正值与真实值间的偏差，RMSE在第一次同化后骤降，修正效果非常明显。此后RMSE缓慢减小并维持在20 m左右，模型预测误差的增长能够得到持续控制。在本案例条件下只需经过前3次的连续同化，就能够使RMSE降低到一个较小的稳定值。

4.2 集合元素个数N对修正效果的影响

EnKF算法中集合元素个数N越多，集合就越能表征状态变量的真实分布，修正效果就会越好[25]。但是集合元素个数越多，模拟次数也会变多，相应地模型计算所花的时间也越多，因此N的值并不是越大越好。本节在4.1节的基础上具体研究了N=5、10、20、30、40五种情况下EnKF算法修正后RMSE的变化，结果如图14所示。

图14 不同集合元素个数下真实组和EnKF组火线间的RMSEFig.14 RMSE of fire perimeters between real group and the group using EnKF algorithm under different ensemble members

图14表明，N=5时，RMSE值并不能像图13中那样不断减小，预测偏差即使在上一时刻得到了修正，下一时刻仍然可能会继续增大。这在很大程度上归因于集合元素产生的随机性，这时样本过少，统计失真，不仅无法反映总体的分布情况，还可能加大统计结果的计算误差。除N=5外，其他RMSE曲线的变化趋势基本一致。当N从5增大到20时，随着集合元素个数增加，同一时刻的RMSE值不断减小，算法对火线的修正效果越来越好。但当N继续增大时，不同N值间RMSE的差异不再明显。N=30时，已经能够获得较为满意的修正效果，继续增大N值对修正效果的改善意义不大。由此可见，N值太小，修正效果会不理想，N值太大计算资源会被浪费，还需根据研究需要，选取合适的集合元素个数。

4.3 观测数据标准差stdobs对修正效果的影响

图15给出了其他参数不变，观测数据标准差stdobs不同时RMSE随时间的变化图。

从图15可以看出，stdobs为60 m、120 m、180 m时，RMSE的变化趋势一致：都是先急剧减小，然后逐步递减，最后缓慢下降。这表现为所有曲线在T=10 h时刻下降幅度最大，但在T=10 h之后曲线下降得越来越缓慢。这是因为初始时刻火源输入的不确定性造成了较大的预测偏差，此时存在较大的算法修正空间。但T=10 h之后，所有时刻的预测结果都以上一时刻的状态修正值作为输入模拟得到的，经过一次同化后，修正值已经进一步接近真实值，之后预测结果的修正空间就会变小。但不同stdobs下算法具备的修正能力是不同的，这表现在不同stdobs下，同一时刻的RMSE值会随stdobs的减小而增大，从整体上来看，RMSE大小与stdobs值呈正相关。stdobs大小代表观测数据的可信度，stdobs越小，观测数据的误差越小，越接近真实值。EnKF算法计算得到的修正值是状态变量的预测值和观测值的加权平均，用标准差更小的观测数据对模型预测结果做修正意味着得到的修正结果的误差也会越小。

图15 观测数据标准差为60 m、120 m、180 m时真实组和EnKF组火线间的RMSEFig.15 RMSE of fire perimeters between real group and the group using EnKF algorithm when the standard deviation of the observation is 60 m,120 m and 180 m

4.4 同化频率对修正效果的影响

前文都是每小时输入一次观测数据，本节主要讨论将同化频率由每小时1次降低为每两小时1次后RMSE的变化，因此仅在T=11 h、T=13 h、T=15 h、T=17 h时刻加入观测数据，总共做了四次修正，结果如图16所示。对比图13可以发现：同化频率减小后，只在同化时刻RMSE的值会减小，其他时刻由于没有进行算法修正，RMSE值仍会增大，但相比于图8完全不做任何修正而言，RMSE的增长还是得到了很大程度的抑制。整个模拟时长内，图13的RMSE值都比对应时刻图8中的值要小。

图16 每两小时同化一次的RMSE曲线Fig.16 RMSE curve with every two hours' assimilation

5 结论

本文论证了复杂森林火灾场景下，利用EnKF算法对因火源位置输入不准确导致的FARSITE火线预测偏差进行修正的可行性。并综合分析了集合元素个数、观测数据标准差和同化频率对修正效果的影响。该研究可以为林火精准扑救提供理论参考，研究提出的FARSITE-EnKF框架也可以推广到其他非线性仿真模型的应用领域。研究初步结论如下：

(1)经EnKF算法逐时同化后，FARSITE模型的火线预测偏差会不断减小并达到一个较低的稳定值，修正结果在位置和形状上都更加接近真实火线。

(2)N过小时，模型预测偏差不能得到很好修正。随着N继续增大，修正效果会越来越好，但N增大到一定值后，对修正效果的影响就不再明显。

(3)逐时同化时，stdobs越小，算法在同一时刻的修正效果越好，RMSE大小与stdobs呈正相关。

(4)同化频率减小后，同化时刻的RMSE会明显减小，未同化时刻的RMSE仍有所增长，但相比于不做任何修正的情况，整体上，误差的增长还是得到了较好抑制。