CN108763811A - 动态数据驱动林火蔓延预测方法 - Google Patents

Abstract

本发明公开了一种动态数据驱动林火蔓延预测方法，通过对火线位置施加扰动产生集合，FARSITE软件实现集合的正向预测，集合卡曼滤波算法则通过融入观测数据对预测结果进行动态修正，提高了FARSITE林火蔓延软件模拟预测的准确性。同时集合卡曼滤波算法是通过Monte‑Carlo方法计算状态的预报误差协方差，用集合的思想解决了实际应用中背景误差协方差矩阵估计和预报困难的问题，易于程序代码编写及调试与维护，可并行计算，实施起来简单高效。

Description

动态数据驱动林火蔓延预测方法

技术领域

本发明涉及林火蔓延模拟预测领域，尤其涉及一种动态数据驱动林火蔓延预测方法。

背景技术

森林火灾是一种失去人为控制，在林地内自由蔓延和扩展的林火行为，具有过火面积大，扑救困难，生态环境影响严重的特点。火势时空蔓延趋势的准确预测是对森林火灾进行及时控制，准确扑救和减少灾害损失的前提和关键。

森林火灾的蔓延扩散是一个多环境因素控制与跨区域尺度的复杂火灾演化过程。目前多依靠计算机仿真模拟的方法对火灾的发展状况进行预测。现有的林火蔓延模型可以划分为三大类：物理/类物理模型；经验模型；半经验模型。其中最广泛使用的是美国Richard C.Rothermel提出的Rothermel半经验模型，火场面积模拟软件(FARSITE)是基于该模型开发的代表性软件，在各国林业及研究部门中得到广泛应用。作为一个确定性模型，FARSITE输入的准确性会显著影响预测结果的可靠性。实际情况下，由于输入数据的缺失和误差，数据处理过程中的近似，模型模拟结果与真实情况会存在一定偏差，从而使得预测结果不完全可靠，因此有必要对模型的预测结果进行修正。

数据同化的概念起源于数值天气预报，它包括动力模型、状态量的直接或间接观测数据、数据同化算法、驱动模型运行的基础参量数据四个基本要素。数据同化算法利用观测数据不断降低模型模拟过程中的不确定性，使得模型模拟预测结果向真实轨迹靠拢。基于此国内外开展了一系列基于林火蔓延模拟的数据同化研究。其中，以遗传算法为代表的数据同化算法通过让模型输入参数(即变量)在一定范围内随机取值来构造个体、利用模拟结果和真实结果的偏差构造适应度函数，模仿自然界的选择与遗传的机理来寻找最优解，即使得模拟结果和真实情况最接近的变量值组合。但是遗传算法编码过程困难，局部搜索能力较差。

发明内容

本发明的目的是提供一种动态数据驱动林火蔓延预测方法，可以用于森林火灾中火线位置的动态预测，能够降低因初始火源位置误差导致的软件预测结果与真实观测值的偏差，使预测的准确性大大提高了，并且，易于程序代码编写及调试与维护，可并行计算，且局部搜索能力较强。

本发明的目的是通过以下技术方案实现的：

一种动态数据驱动林火蔓延预测方法，包括：

步骤P1、获取带有时间和空间位置信息的森林火灾热红外图像；

步骤P2、对获取到的森林火灾热红外图像进行处理得到火线位置观测数据；

步骤P3、对输入的火线位置施加扰动，利用FARSITE模拟生成直至t’时刻的时间序列的火场火线位置；

步骤P4、利用集合卡曼滤波数据同化算法并结合火线位置观测数据对火场火线位置进行修正；

步骤P5、若在t’+Δt时刻还有火线位置观测数据，则重复步骤P3～步骤P4，其中的Δt为时间步长的整数倍；若否，则根据火场火线位置修正结果进行FARSITE模拟预测；

步骤P6、对预测结果进行可视化显示。

由上述本发明提供的技术方案可以看出，通过对火线位置施加扰动产生集合，FARSITE软件实现集合的正向预测，集合卡曼滤波算法则通过融入观测数据对预测结果进行动态修正，提高了FARSITE林火蔓延软件模拟预测的准确性。同时集合卡曼滤波算法是通过Monte-Carlo方法计算状态的预报误差协方差，用集合的思想解决了实际应用中背景误差协方差矩阵估计和预报困难的问题，易于程序代码编写及调试与维护，可并行计算，实施起来简单高效。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种动态数据驱动林火蔓延预测方法的流程图；

图2为本发明实施例提供的森林火灾热红外图像的处理流程图；

图3为本发明实施例提供的FARSITE软件输入输出示意图；

图4为本发明实施例提供的基于FARSITE软件和集合卡曼滤波算法的数据同化过程示意图；

图5为本发明实施例提供的集合卡曼滤波算法计算流程图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

本发明实施例提供一种动态数据驱动林火蔓延预测方法，如图1所示，其主要包括如下步骤：

步骤P1、获取带有时间和空间位置信息的森林火灾热红外图像。

本发明实施例中，可以通过无人机勘察的方式，利用机载热红外成像仪，拍摄整个火场范围的热红外视频图像，同时搜集火场的真实地理位置信息；并对热红外视频图像进行处理，为其附加时间信息和空间地理位置信息，从而获得带有时间和空间位置信息的森林火灾热红外图像。

步骤P2、对获取到的森林火灾热红外图像进行处理得到火线位置观测数据。

本发明实施例中，对步骤P1获取的热红外视频图像中选取质量较好的帧进行边缘检测，边缘细化，和点的离散化处理。边缘检测和边缘细化处理是为了得到由单像素构成的火线轮廓，点的离散化处理则是为了匹配FARSITE输出文件的格式。

如图2所示，为针对森林火灾热红外图像进行处理的流程图。

1)利用Canny边缘检测算子对森林火灾热红外图像中图像质量符合要求的帧进行多个步骤的复核检测：

a、图像平滑：首先采用二维高斯函数的任意方向上的一阶方向导数为噪声滤波器，与图像卷积进行滤波。

b、梯度计算：用一阶偏导的有限差分来计算梯度的幅值，以便进一步确定图像的边缘。

c、非极大值抑制：为确定边缘，在全局梯度的基础上保留局部梯度最大的点，而抑制非极值。

d、双阈值检测：选取高阈值Th和低阈值Tl，然后开始扫描图像，从候选边缘图像中确定边缘点，非边缘点和疑似边缘点，对疑似边缘点依据边缘的连通性对其进行判断，最终检测和连接出最终的边缘。

本领域技术人员可以理解，图像质量符合要求的帧可以从各个指标来确定，例如，图像清晰度等。

2)通过Hilditch算法对检测到的森林火灾热红外图像边缘进行细化，经过边缘检测之后的得到的火线轮廓是具有一定像素宽度的线条状区域，通过Hilditch算法设定目标像素为背景值的条件，每次扫描删除图像上满足该条件的目标轮廓像素，直到图像上不存在可删除的轮廓像素为止，即将二值图像中具有一定宽度的线条状区域变成一条薄线，形成一个像素的宽度，得到由单像素组成的火线轮廓。

3)对火线上的坐标点进行离散：森林火灾热红外图像经上述处理后得到由单像素构成的火线轮廓，确定火线轮廓上每个像素点对应空间位置中心的x,y坐标，得到表征火线轮廓的全部点的坐标；为实现集合卡曼滤波算法实施过程中模型输出结果和观测数据的匹配，需将所获取的点进一步离散，可以处理成点间距相等的固定点数S的点集，其中S的个数可以自定义，S的大小表征了观测数据的密集程度，也会影响同化效果的优劣，点数越多，吸收的观测信息越多。最后仿照FARSITE火线位置输出的格式，按照逆时针对离散点进行数字标记，得到火线位置观测数据。

步骤P3、对输入的火线位置施加扰动，利用FARSITE模拟生成直至t’时刻的时间序列的火场火线位置。

本步骤中，首先定义一个初始点火位置，假定真实火源位置并不清楚，对点火位置做一个近似猜测，以此来模拟火源位置的初始输入误差。

所输入的火线位置为一组矢量shp文件，也即基于该初始点火位置的一组扰动集合；扰动通过状态参量(x₁,y₁,x₂,y₂···x_m,y_m)的均值和协方差两个特征值施加，施加方式为：

其中，状态参量为随机变量，m为t时刻火线边界上离散火点的个数，(x _j,y_j)为火线边界上第j个离散火点的坐标，j＝1,2,...,m；表示状态参量在t时刻的分析值的平均，表示状态参量的分析误差协方差矩阵，x_i,t是从以为均值，为方差的多维正态分布中抽样产生的第i个集合元素，i＝1,2,..,n，n为取样个数，也就是集合数。初始时刻t＝0时，即代表初始点火位置，此时和参数的大小都可以根据要求自定义。此后t不为0时刻，和的大小根据计算求取。

然后对所有集合利用FARSITE软件进行正向模拟，得到下一时刻的火线位置输出，此时FARSITE软件模拟过程可以概括为一个状态转换过程表示为：

其中，表示第t’时刻第i个集合元素的预测值，也即对t时刻基于火线位置状态分布产生的初始集合中的每一个元素做模型预测得到的对应的t’时刻的火线位置预测结果；M_t,t’表示从t到t’时刻的状态变化关系，即对应于FARSITE模型将火线位置坐标从t时刻的值向t’时刻的值进行求解计算的过程，也即状态参量的预测过程，上标f表示模型预测值，下标t、t’分别对应t时刻、t’时刻(二者的关系可以理解为上一时刻与当前时刻)。

如图3所示，为FARSITE软件输入输出示意图：为了满足编程需要，本发明中使用的FARSITE软件为Linux版本。Linux版本的FARSITE软件是一个经过编译链接后的可执行程序，通过在终端命令行输入指令调用程序完成林火蔓延的模拟过程。Linux版本的输入文件类型包括：txt模型运行文件，lcp地形植被文件，input输入参数文件，火线位置的shp文件，障碍物位置的shp文件。输出文件有26个，输出参量包括：火线位置，火线强度，火焰到达时间，蔓延速度等。其中前四个输入文件是必须的，本发明中输出参量只关注火线位置，其输出文件格式也是shp格式。其中input输入参数文件中可以定义模拟所需的各种时空参数，包括时间步长ΔT和火线上离散点的点间距ΔP。其中火线位置的shp文件是一个包含离散点的数据集，这些点组成的拓扑结构即火源的形状。根据输入的火线位置shp文件，FARSITE利用Rothermel模型计算火线上各点蔓延速度，各点以椭圆形状，按照Huygens波动方式向四周蔓延，给定时间步长ΔT内所有波动椭圆的包络面就构成了下一时刻的火线蔓延边界。

步骤P4、利用集合卡曼滤波数据同化算法并结合火线位置观测数据对火场火线位置进行修正。

火场火线位置进行修正的过程主要包括：

预报误差的分析：

上式中，为t’时刻状态量的n个集合元素的预测值的均值；H为状态转移矩阵，实现预测值向观测数据空间的映射；表示第t’时刻第i个集合元素的预测值的误差协方差矩阵，其含义可以用计算公式来理解：表示预测和观测值之间的协方差矩阵，这里通常将当成一个整体理解，因为H^T通常不是线性矩阵；表示将预测值映射到观测空间的过程，当于对进行和步骤P2相同的离散过程，从其中获取同样的数目为S的等间距点集。

卡曼增益的计算：

上式中，K_t'为t’时刻卡曼增益矩阵；

火场火线位置及状态参量修正：

上式中，为第i个集合的修正结果，从计算公式可以看出，是通过在第t’时刻(即当前时刻)预测的基础上对第t时刻(即上一时刻)的观测数据和预测的偏差)进行系数修正K_t'得到；为状态参量的修正结果；和一样，是通过均值和方差的形式得到的，经过步骤P2离散后的t时刻的火线位置观测数据作为通过数据采集过程中的误差分析，得到观测值的方差均值通过正态分布得到的t时刻第i个火线位置观测数据集合元素

协方差矩阵的计算：

如图4所示，为本发明基于FARSITE软件和集合卡曼滤波算法的数据同化过程示意图，数据同化算法作为数据同化的重要组成部分，是连接观测数据和模型模拟预测的关键。如图4所示，首先在t＝0时刻输入真实火场的地形植被条件、气候条件文件，设置控制参数和初始火源的火线位置，然后FARSITE按照时间步长ΔT进行逐步积分，完成模型的正向模拟过程，输出指定时刻t’的火线位置的预测值。同时当t’时刻可以获取到较高质量热红外图像帧时，则可以通过P2步骤处理，得到观测数据。利用集合卡曼滤波数据同化算法，将观测数据和模型输出进行比较，根据观测量和输出量的协方差矩阵所定义的可信度的大小分别赋予它们一定的权重，对模型输出进行修正，完成单步的同化过程。因为FARSITE的火线位置的输入输出文件都是shp格式，因此可将修正后的输出同时作为t’时刻输入，继续此后的预测和同化过程。其中是否有观测数据输入的判定条件是能否从热红外图像中处理得到高质量的火线位置的离散点数据。可根据需要自行设定观测数据输入的时间。

图5为本发明集合卡曼滤波算法计算流程图，t＝0时刻将自定义的初始点火位置和协方差矩阵作为多维正态分布的均值和方差进行抽样产生n个初始火线位置的集合(x_1,0,x_2,0,x_3,0,···,x_n,0)，然后FARSITE将所有集合向前积分，得到t’时刻的n个预报值的集合为实现预报值的修正更新，利用和计算得到t’时刻卡曼增益矩阵K_t'，其中表示对所有输出进行和步骤P2一样的点的离散处理过程。最后利用K_t'，利用相应公式计算得到n个集合分析值，最后对分析值的集合取平均得到状态参量的最终修正结果。其中(i～(1-n))表示n个观测值的集合，其构造方法同x_i,0。最后为继续下一步的同化过程，需计算即分析误差协方差矩阵。

步骤P5、若在t’+Δt时刻还有火线位置观测数据，则重复步骤P3～步骤P4，其中的Δt为时间步长的整数倍；若否，则根据火场火线位置修正结果进行FARSITE模拟预测。

得到后，若在t’+Δt时刻还有火线位置观测数据则将作为输入，即将作为t＝0时刻的对继续施加扰动，重复P3～P4过程；若无观测值，则继续FARSITE的预测过程。

步骤P6、对预测结果进行可视化显示。

将预测结果统一处理并导入到地理信息处理软件，与地理信息数据叠加显示；具体的可以将FARSITE输出的预测结果，提取的观测数据和预测结果的修正值统一处理成一定文件格式，将其导入到Windows版本FARSITE，GIS，QGIS，Google Earth等软件，与地形，高程，坡度，植被分布等地理信息数据叠加显示，直观地展示森林火灾的火势蔓延情况，同时还可以对同化效果进行定性比较。例如可以将上述数据处理成shp格式文件，通过Windows版本FARSITE软件中的矢量导入功能，将结果叠加显示在lcp文件上

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.一种动态数据驱动林火蔓延预测方法，其特征在于，包括：

步骤P1、获取带有时间和空间位置信息的森林火灾热红外图像；

步骤P2、对获取到的森林火灾热红外图像进行处理得到火线位置观测数据；

步骤P3、对输入的火线位置施加扰动，利用FARSITE模拟生成直至t’时刻的时间序列的火场火线位置；

步骤P4、利用集合卡曼滤波数据同化算法并结合火线位置观测数据对火场火线位置进行修正；

步骤P5、若在t’+Δt时刻还有火线位置观测数据，则重复步骤P3～步骤P4，其中的Δt为时间步长的整数倍；若否，则根据火场火线位置修正结果进行FARSITE模拟预测；

步骤P6、对预测结果进行可视化显示。

2.根据权利要求1所述的一种动态数据驱动林火蔓延预测方法，其特征在于，所述获取带有时间和空间位置信息的森林火灾热红外图像包括：

拍摄整个火场范围的热红外视频图像，并对热红外视频图像进行处理，为其附加时间信息和空间地理位置信息，从而获得带有时间和空间位置信息的森林火灾热红外图像。

3.根据权利要求1所述的一种动态数据驱动林火蔓延预测方法，其特征在于，所述对获取到的森林火灾热红外图像进行处理得到火线位置观测数据的步骤包括：

利用Canny边缘检测算子对森林火灾热红外图像中图像质量符合要求的帧进行多个步骤的复核检测：图像平滑、梯度计算、非极大值抑制与双阈值检测；再通过Hilditch算法对检测到的森林火灾热红外图像边缘进行细化，即将二值图像中具有一定宽度的线条状区域变成一条薄线，形成一个像素的宽度；

森林火灾热红外图像经上述处理后得到由单像素构成的火线轮廓，确定火线轮廓上每个像素点对应空间位置中心的x,y坐标，得到表征火线轮廓的全部点的坐标；将所获取的点进一步离散，处理成点间距相等的固定点数S的点集，并按照逆时针进行数字标记，得到火线位置观测数据。

4.根据权利要求1所述的一种动态数据驱动林火蔓延预测方法，其特征在于，所述对输入的火线位置施加扰动包括：

首先定义一个初始点火位置，所输入的火线位置为一组矢量shp文件，也即基于该初始点火位置的一组扰动集合；扰动通过状态参量(x₁,y₁,x₂,y₂···x_m,y_m)的均值和协方差两个特征值施加，施加方式为：

其中，状态参量为随机变量，m为t时刻火线边界上离散火点的个数，(x_j,y_j)为火线边界上第j个离散火点的坐标，j＝1,2,...,m；表示状态参量在t时刻的分析值的平均，表示状态参量的分析误差协方差矩阵，x_i,t是从以为均值，为方差的多维正态分布中抽样产生的第i个集合。

5.根据权利要求4所述的一种动态数据驱动林火蔓延预测方法，其特征在于，利用FARSITE模拟的过程即为一个状态转换过程，表示为：

其中，表示第t’时刻第i个集合元素的预测值，M_t,t'表示从t到t’时刻的状态变化关系，即对应于FARSITE模型将火线位置坐标从t时刻的值向t’时刻的值进行求解计算的过程，也即状态参量的预测过程，上标f表示模型预测值。

6.根据权利要求5所述的一种动态数据驱动林火蔓延预测方法，其特征在于，所述利用集合卡曼滤波数据同化算法并结合火线位置观测数据对火场火线位置进行修正的过程包括：

预报误差的分析：

上式中，为t’时刻状态量的n个集合元素的预测值的均值；H为状态转移矩阵，实现预测值向观测数据空间的映射；表示将预测值映射到观测空间的过程；表示第t’时刻第i个集合元素的预测值的误差协方差矩阵；

卡曼增益的计算：

上式中，K_t'为t’时刻卡曼增益矩阵；

火场火线位置及状态参量修正：

上式中，为通过正态分布得到的t时刻第i个火线位置观测数据集合元素，为第i个集合的修正结果，为状态参量的修正结果；

协方差矩阵的计算：

7.根据权利要求6所述的一种动态数据驱动林火蔓延预测方法，其特征在于，若在t’+Δt时刻还有火线位置观测数据则将作为输入，即将作为t＝0时刻的

8.根据权利要求1所述的一种动态数据驱动林火蔓延预测方法，其特征在于，所述对预测结果进行可视化显示包括：将预测结果统一处理并导入到地理信息处理软件，与地理信息数据叠加显示。