基于时空融合技术的森林火灾遥感动态监测

黄武彪，栾海军，李大成

（1.厦门理工学院计算机与信息工程学院，福建厦门 361024；2.长安大学地质工程与测绘学院，陕西西安 710064；3.厦门理工学院数字福建自然灾害监测大数据研究所，福建厦门 361024；4.太原理工大学矿业工程学院，山西太原 030024）

引言

森林火灾是一种突发性强、破坏性大、救援困难的自然灾害，全球每年平均发生森林火灾数十万次，森林受灾面积达数百公顷［1］，对生态环境和经济发展造成了极大威胁。截止2019年，全国有记录森林火灾2 345次，其中重大火灾为8次，特别重大火灾为1次，火场总面积为39 705 hm2，受害森林面积为13 505 hm2，伤亡人数为76人，其它损失折款16 219.9万元［2］。森林火灾监测，要求及时发现着火点的位置及其变化、过火面积，并准确评估出火灾损失及影响，因此需要具有高时间分辨率和高空间分辨率的遥感影像来进行分析判读。但在现有传感器硬件条件限制下，卫星遥感数据无法同时满足高空间分辨率和高时间分辨率的要求［3-4］。因此，很多学者提出时空融合的技术方法来解决遥感传感器时间分辨率和空间分辨率的矛盾［5-8］，将该方法应用于森林火灾场景中，便于更加准确、快速得获得受灾区域高空间和高时间分辨率的遥感影像，为森林火灾演化的监测和灾损评估提供有力支撑。Gao等［9］提出的时空自适应反射率融合模型（Spatial and Temporal Adaptive Reflectance Fusion Model，STARFM），用于融合Landsat影像和MODIS影像得到高时空分辨率的数据，取得了较好的效果，该模型适用于同质地表季节变化情形。Hilker等［10］在提出了一种针对反射率变化的时空自适应融合模型算法（Spatial TemporalAdaptive Algorithm for mapping Reflectance Change，STA⁃ARCH），从低分辨率影像的密集时间序列中检测出变化点，以提高土地覆盖类型变化时的STARFM性能，适用于地表反射率突发扰动事件情形。邬明权等［11］提出基于混合像元分解的方法（Spatial and Temporal Data Fusion Model，STDFM）来融合MODIS和Landsat影像数据，进一步，Zhang等［12］对STDFM方法进行了改进，提出增强型基于混合像元分解的方法（Enhanced spatial and temporal data fusion model，ESTDFM），该模型适用于异质地表季节变化情形。Zhu等［13］基于STARFM提出了一种增强型时空自适应反射率融合模型（En⁃hanced Spatial and Temporal Adaptive Reflectance Fusion Model，ESTARFM），引入了一个转换系数，可以更好地预测异质性地表的反射率的变化。Huang等［14］提出的基于稀疏表示的时空反射率融合模型（SParse⁃repre⁃sentation⁃based SpatioTem⁃poral reflectance FusionModel，SPSTFM），将稀疏表达理论引入时空融合算法，该模型适用于地表季节和类别变化情形。除地表反射率参数外，学者们［15-17］对地表温度参数也进行了时空融合研究。学者们致力于研发出通用性与鲁棒性更优的时空融合算法。如，Cheng等［18］提出了一种时空非局部滤波融合模型（Spatial and Temporal Nonlocal Filter⁃Based Data FusionMethod，STNLFFM），对于异质性地表区域具有更高的预测精度。Zhao等［19］提出了一种针对复杂地表变化的鲁棒性自适应时空数据融合模型（Ro⁃bust Adaptive Spatial and Temporal Fusion Model，RASTFM），在捕捉地表变化现象时具有更高的准确度和鲁棒性。黄波和姜晓璐［20］提出一种增强型空间像元分解时空遥感影像融合算法（Unmixing Enhanced model for Spatial and temporal image fusion，EUSTFM），能够实现对季节性变化及复杂的地物类型变化的稳定预测，生成具有更高精度的融合影像。

鉴于森林属于地表相对均质的场景，故本次研究将选用经典的STARFM算法与具有组合优势的基于地物内组分时相变化模型的地表反射率时空融合算法，并综合使用多种空间分辨率更优（≤30m）的国内外传感器影像，对两种算法所得时空融合结果优化组合，进而应用于2019年3月30日四川凉山州木里县森林火灾遥感动态监测中。

1 研究区与研究数据

1.1 研究区概况

本次研究区域（如图1所示）为四川凉山木里县雅砻江镇立尔村附近的雅砻江边海拔约3 800 m的一处森林中，区域经纬度介于28°32′10″N～28°33′27″N，101°15′02″E～101°16′49″E之间。该区域位于典型的高山峡谷区，气候特点为冷热两季交替、干湿分明，天气干燥，日照强，大面积被植被覆盖，有少部分水及岩石。地形复杂、坡陡谷深，交通、通讯不便。火灾发生时间为2019年3月30日-4月4日。此次火灾造成30名消防员牺牲，据中国新闻网2019年4月5日报道，本次火灾起火点为一颗位于山脊上的树龄约80年的云南松，起火原因为雷击火［21-22］。

图1 研究区位置Fig.1 Location of the study area

1.2 研究数据

本次研究所用实验数据源为MOD09GA、Landsat8 OLI、Sentinel-2及GF-1 WFV影像，具体实验影像说明见表1，表中所列影像均为本次实验中所用到的数据。各遥感数据的特征与使用波段如表2所示。

表1 实验数据说明Table 1 Experimental data

表2 各遥感数据的特征及使用波段Table 2 The characteristics and use bands of various types of remote sensing data

Landsat8 OLI数据和Sentinel-2数据可从USGS（https：//earthexplorer.usgs.gov/）免费获取，MOD09GA数据可从NASA（https：//ladsweb.modaps.eosdis.nasa.gov/search/）免费获取，GF-1 WFV数据可从中国资源卫星应用中心陆地观测卫星数据服务平台（http：//218.247.138.119：7777/DSSPlatform/productSearch.html）免费获取。

在进行时空融合之前，需要对获取的影像数据使用遥感图像处理软件进行预处理。使用ENVI对Land⁃sat8 OLI影像进行辐射定标和大气校正，生成反射率数值影像。Sentinel-2影像使用Sen2Cor和SNAP软件进行大气校正并转换成ENVI所能打开的img格式。GF-1WFV影像除进行与Landsat8 OLI相同的预处理外，还需先使用ENVI遥感图像处理软件进行几何纠正，其辐射定标时绝对辐射定标系数在中国资源卫星应用中心网站上下载。对于MOD09GA影像需先使用MODIS Reprojection Tool（MRT）重投影至UTM/WGS84坐标系、GeoTIFF格式，并重采样为10 m、16 m和30 m分辨率。除上述预处理外，还应使用ENVI遥感图像处理软件进行各传感器影像裁剪，确保其研究区域范围相同。在本次研究中，均针对各影像的单波段进行研究。

2 研究方法

具体研究技术路线图如图2所示。主要技术流程阐述如下：

图2 本研究的技术路线图Fig.2 Flowchart of the research

（1）对获取到的MOD09GA、Landsat8 OLI、Sentinel-2、GF-1WFV影像数据进行预处理；

（2）使用STARFM算法与基于地物内组分时相变化模型的地表反射率时空融合算法分别对MOD09GA与Landsat8 OLI、Sentinel-2、GF-1 WFV影像进行融合，以达到结合多源影像的目的，从而生成待预测时刻的中等空间分辨率影像；

（3）基于预测影像计算火灾指标因子，进行火灾演化趋势分析。

2.1 面向多源数据的时空自适应反射率融合模型

时空自适应反射率融合模型（STARFM）在使用前需要对来自不同平台的观测数据进行几何配准，并进行大气校正转化至地表反射率。STARFM算法的总体流程描述如下［9］。STARFM算法的基本原理是：假定在tk日低分辨率影像数据地表反射率值M（xi，yj，t k）与中高分辨率影像数据地表反射率值L（xi，yj，tk）之间的关系可以表示为

则在t0日的低分辨率影像数据地表反射率值M（xi，yj，t0）与在t0日的中高分辨率影像数据地表反射率预测值L（xi，yj，t0）的关系可以表示为

εk与ε0分别表示观测到的低分辨率影像数据和中高分辨率影像数据地表反射率之间由于不同的波段宽度和太阳几何造成的差异。假设在预测日期t0和日期tk时像素（xi，y j）的地面覆盖类型和系统误差不变，则ε0=εk，因此

但这仅仅是一种理想情况，它们之间的关系受到以下3方面的影响：1）观测的低分辨率影像数据在用与高分辨率影像数据相同的空间分辨率考虑时可能包括混合的土地覆盖类型；2）在预测期间，土地覆盖可能从一种类型变为另一种类型；3）土地覆盖状态和太阳几何双向反射率分布函数的变化将改变从预测日期t0到日期tk的反射率。

因此，通过引入邻近像素的附加信息，使用加权函数来计算日期t0时中心像素的地表反射率：

其中ω是搜索窗口的大小是该移动窗口的中心像素。为确保使用来自邻近像素的正确信息，仅使用来自同一光谱类别的和来自移动窗口内中高分辨率影像数据地表反射率的无云像素来计算反射率。权重Wijk代表每个邻近像素对中心像素的预测反射率的贡献程度，从光谱、时间、距离3个方面来确定每个光谱相似像素的最终权重。

本次研究所做改进：上述内容通常针对单一类型中等空间分辨率传感器，当面向多源数据时，需要对STARFM时空融合模型进行优化组合应用。具体优化组合策略设计如下：在经典的STARFM算法基础上，在30 m Landsat卫星16天或更长的预测周期内利用多种国内外传感器影像（空间分辨率≤30 m）进行分段，分段原则为：时间最邻近及空间分辨率优先。上述原则以案例描述为：（1）若引入10 m Sentinel-2可见光影像，该影像在Landsat基期与末期时间段之间（如第5天），则将预测区间分为两段（Landsat基期-第5天）、（第6天-第16天）；（2）对于每一时间段（包含新的基期和末期两天），若天数为偶数，对其均分，前半段影像预测由本时间段基期影像完成，后半段影像预测由本时间段末期影像完成，若天数为奇数，前半段（如Landsat基期-第2天）的影像预测由“最邻近”的Landsat基期影像完成，后半段（如第4天-第5天）的影像预测由“最邻近”的Sentinel-2影像完成，中间一天（如Landsat基期-第5天时间段，中间一天为第3天）的影像预测依照空间分辨率优先原则，由空间分辨率高的传感器影像（也就是Sentinel-2影像）完成；（3）当引入的传感器影像更多时，参照上述案例进行扩展即可。各时间分段内，当新的基期或末期影像中高空间分辨率优于30 m，需要通过现有商业遥感软件进行空间重采样至30 m，然后进行预测影像生成。

2.2 面向多源数据的基于地物内组分时相变化模型的地表反射率时空融合算法

对于相同的土地覆盖类型的均质像素都具有相同的反射率，并且这些像素的季节性和双向反射率变化对于不同的像元大小也应该相同［23］。对于不同分辨率的影像，每种土地覆盖类型在获取日期和预测日期之间的变化关系大致相同，即在不同组分之间具有空间尺度不变性。因此，可以将基于获取日期的低分辨率影像数据和中高分辨率影像数据建立的变化关系映射到中高分辨率遥感影像中，从而预测出待预测时刻的中高分辨率影像。

基于上述前提，提出了基于地物内组分时相变化模型的地表反射率时空融合算法，具体流程如图3所示。先对中高空间分辨率遥感影像采用迭代自组织数据分析算法（ISODATA）进行非监督分类，将其分为多个地物类别，得到不同的地物组分。然后分别对各个地物组分的变化，建立相应时相变化模型进行预测。

图3 基于地物内组分时相变化模型的地表反射率时空融合算法流程图Fig.3 Flow chart of Spatiotemporal fusion algorithm of land surface reflectance basedon temporal phase change model of components in surface features

寻找“纯像元”的基础为将非监督分类得到的结果聚合到低分辨率影像像元大小，其聚合过程可理解为通过计算聚合前后的分辨率像元的比例n，将非监督分类影像自左上开始，组成若干n*n大小的矩阵，并记录聚合后的行列号。分别计算在每一个n*n矩阵中各组分所占比例，取比例最大的组分作为该像元的类别，并判断比例值是否大于20%（20%为经验值，针对不同的分辨率设置不同的阈值），若大于该值，则将该聚合后的像元称为该地类的一个“纯像元”。依次得到各类地物的“纯像元”，并找出两个时间的低分辨率影像中对应组分内部的无云像元是如何变化的，使用差值、比值或变化率建立不同地物组分的时相变化模型。整个过程基于MATLAB软件实现。

基于地物内组分时相变化模型的地表反射率时空融合算法三类模型依次描述如下。

差值模型：

比值模型：

变化率模型：

式中，M OD09G A_T0表示在T0时刻低分辨率遥感影像的地表反射率值，M O D09 G A_T1表示在T1时刻低分辨率遥感影像的地表反射率值，S entinel_T0表示在T0时刻的中高分辨率遥感影像地表反射率值，pre_Sen⁃tinel_T1表示预测得到的T1时刻的中高分辨率遥感影像地表反射率值，mean表示对所有的值取平均。

本次研究所做改进：上文所述算法通常针对单一类型中等空间分辨率传感器，对于多源数据的实施原则可参照2.1节“优化组合策略”，其实施过程相似，但无需“重采样至30 m”步骤。

2.3 火灾指标因子介绍

在森林火灾监测中应用遥感影像，需要对使用两种时空融合算法得到的影像数据计算燃烧面积指数和归一化燃烧指数。

2.3.1 燃烧面积指数

燃烧面积指数（Burn Area Index（BAI））是采用影像的红色波段（Red）和近红外波段（N I R）来增强火烧以后的地表信息，也就是增强火灾过火以后图像上的木炭信号，燃烧区域的BAI值比较大。计算公式［24］如下：

在ENVI5.3中利用光谱指数计算工具选择“Burn Area Index”即可计算研究区燃烧面积指数。

2.3.2 归一化燃烧指数

归一化燃烧指数（Normalized Burn Ratio（N B R））是基于近红外波段（N IR）和短波红外波段2（SWI R2）来增强较大范围的火灾区域，燃烧区域的NB R值比较小，计算公式［25-26］如下：

在ENVI5.3中利用光谱指数计算工具选择“Normalized Burn Ratio”即可计算研究区归一化燃烧指数。

3 分析与讨论

3.1 融合结果分析

利用STARFM算法和基于地物内组分时相变化模型的地表反射率时空融合算法分别对MOD09GA与Landsat8 OLI、Sentinel-2及GF-1 WFV进行时空融合。根据融合结果（见表3），将预测后的影像与可从其他类型传感器获得的同时刻的影像进行目视对比，分析本次实验所采用的时空融合算法的精度，并对火灾发生前后的火灾指标因子进行提取分析。如图4是3月22日GF-1影像与MOD09GA影像预测3月21日GF-1影像，并与Sentinel-2影像对比，图5是3月18日Landsat8 OLI影像与MOD09GA影像预测3月26日Landsat8 OLI影像，并与Sentinel-2影像对比，图6是3月26日Sentinel-2影像与MOD09GA影像预测3月30日Senti⁃nel-2影像，并与GF-1影像对比。

图4 使用3月22日GF-1与MODIS影像预测3月21日GF-1影像，并与Sentinel-2影像对比Fig.4 Use GF-1 and MODIS images on March 22 to predict GF-1 image on March 21 and compare with Sentinel-2 images

图5 3月18日Landsat与MODIS影像预测3月26日Landsat影像，并与Sentinel-2影像对比Fig.5 Use Landsat and MODIS images on March 18 to predict Landsat image on March 26 and compare with Sentinel-2 images

图6 3月26日Sentinel-2与MODIS影像预测3月30日Sentinel-2影像，并与GF-1影像对比Fig.6 Use Sentinel-2 and MODIS images on March 26 to predict Sentinel-2 image on March 30 and compare with GF-1 images

表3 使用两种算法的融合结果Table 3 Fusion results based on the two algorithms

两种算法在不同类型遥感数据融合应用中各有其局限性。将时空融合后的单波段影像进行叠加，得到RGB彩色图像，如图4～6所示。基于图4～6分析：使用STARFM算法对地表变化明显区域预测精度差，从图4（e）可发现影像左侧存在黑色区域，与实际地表存在明显差距，说明STARFM算法对GF-1 WFV的预测周边裸露岩石的效果不是很好；从图5（e）、图6（e）与图5（d）、图6（d）的山脉、植被、岩石等地区进行对比可以发现STARFM对于Landsat8 OLI和Sentinel-2的预测效果较好，但Landsat8 OLI的结果更优。基于地物内组分时相变化模型的地表反射率时空融合算法包括3种变化模型，将图4（f）、图5（f）和图6（f）与对应的（d）影像选取山脉、植被、岩石等地区对比，发现使用地物组分时相变化差值模型来进行预测时，对Landsat8 OLI、Sentinel-2和GF-1预测效果都比较好；将图4（g）、图5（g）和图6（g）与对应的（d）影像对比，发现使用地物组分时相变化比值模型进行预测时，图4（g）、图5（g）影像图左侧岩石地区呈现黄色，与真实地表存在明显差异，因此地物组分时相变化比值模型对GF-1和Landsat8 OLI的周边预测效果较差，对Sentinel-2的预测效果较好；将图4（h）、图5（h）和图6（h）与对应的（d）影像对比，发现使用地物组分时相变化的变化率模型进行预测时，图4（h）影像图左侧裸露岩石呈现明显差异，因此地物组分时相变化的变化率模型对GF-1预测效果较差，Sentinel-2和Landsat8 OLI预测效果较好。综合分析，对于Landsat8 OLI的预测可以选用STARFM和地物组分时相变化差值模型，对于Sentinel-2的预测，地物组分时相变化差值模型和变化率模型要好于STARFM和比值模型，对于GF-1的预测，尽量选择地物组分时相变化差值模型。因此，通过影像对比与分析，基于地物内组分时相变化的差值模型对于上述3种影像均适用，而使用比值模型进行预测的效果都不是很好。

对不同传感器中等空间分辨率影像与相同时间不同算法预测所得结果进行比对分析，计算相同波段上的影像相关系数，结果如表4所示。表4中，相关系数越大表明算法在该类型传感器该波段上的预测效果越好。可以看出，对于GF-1 WFV和Sentinel-2，在预测的每个波段，基于地物组分时相变化的差值模型的相关系数最高；对于Landsat8 OLI，在近红外和短波红外1两个波段上，基于地物组分时相变化的变化率模型更优。同时，不同的算法在不同传感器影像的同一波段上的预测效果也不同，具体为：1）对于蓝光波段，STARFM、基于地物组分时相变化的比值模型和变化率模型对Sentinel-2的预测效果最好，差值模型对GF-1 WFV的预测效果最优；2）对于绿光波段，STARFM和基于地物组分时相变化的比值模型对Sentinel-2的预测效果最好，差值模型对GF-1 WFV的预测效果最优，变化率模型对Landsat8 OLI的预测效果更好；3）对于红光波段，STARFM、基于地物组分时相变化的差值模型、比值模型和变化率模型分别对Landsat8 OLI、GF-1 WFV、Sentinel-2和Landsat8 OLI的预测效果最好；4）对于近红外波段，四种算法的预测效果最好的传感器影像分别为Landsat8 OLI、GF-1 WFV、Landsat8 OLI、Landsat8 OLI；5）对于短波红外1和2，4种算法在Landsat OLI上的预测效果均最好。对于提取2种火灾指标因子所用波段，可选用差值模型进行预测，其对几个波段均适用，其次是变化率模型，最后是比值模型和STARFM模型。

表4 各波段不同时空融合方法的相关系数Table 4 Correlation coefficients of different fusion methods in each band

3.2 火灾变化时序分析

对2种时空融合算法所得结果进行优化组合，获取质量更优的逐日中等空间分辨率预测影像，进而对预测影像计算火灾指标因子，由于GF-1 WFV影像无短波红外波段，故无法计算归一化燃烧指数，以“时间最邻近及空间分辨率优先”为原则选取符合条件的Landsat8 OLI/Sentinel-2预测影像作为替代计算该指数。经过计算，如图7所示是火灾发生前期、中期、后期3个时间的燃烧面积指数变化，图8所示是火灾发生前期、中期、后期3个时间的归一化燃烧指数对比。

图7 火灾发生前期、中期、后期的燃烧面积指数Fig.7 Burn Area Index before，during and after the fire

图8 火灾发生前期、中期、后期的归一化燃烧指数Fig.8 Normalized Burn Ratio before，during and after the fire

图7（a）中植被覆盖区域呈现黑色或暗灰色，图7（c）中植被覆盖区域（未发生火灾区域）呈现黑色，而火灾发生区域呈现白色，显著的对比度差异利于确定火灾发生以后的位置和面积；但图7（b）研究区呈现黑色，无法看到燃烧区域，这是由于在火灾发生时有植被燃烧产生的浓烟，而燃烧面积指数计算所使用的红波段（Red）和近红外波段（NIR）均无法穿透烟雾，所以对于燃烧面积指数会产生一定的影响。图8（a）中植被覆盖区域呈现白色或浅灰色，图8（c）中植被覆盖区域（未发生火灾区域）呈现白色或浅灰色，而火灾发生区域呈现黑色，明显的黑白对比可确定火灾发生以后的位置和面积；图8（b）中可以看到火灾发生时零星的燃烧区域，这是由于归一化燃烧指数计算所使用的短波红外波段（SWIR）可以穿透烟雾，GF-1 WFV影像没有短波红外波段，这也造成了GF-1在火灾指标因子分析中的局限性。从图7整体来看，未发生火灾区域仍然保持黑色或暗灰色，而发生火灾区域变为白色；图8则与图7相反，从整体来看未发生火灾区域仍然保持白色或浅灰色，而发生火灾区域变为黑色。且总体而言，火灾演化进程中NBR的变化更加显著，BAI效果略差。通过两种火灾指标因子的对比研究，结果表明：基于时空融合影像分析火灾演化态势时，归一化燃烧指数计算结果更敏感、更有效。

火灾发生时间为2019年3月30日-2019年4月4日。基于预测得到的影像并计算NBR，统计过火面积，我们发现：由于2019年3月30日影像获取时火灾还未发生，故NBR结果中没有火灾区域；2019年3月31日过火面积约为72 087.96m2；2019年4月1日过火面积约为187 616.86m2；2019年4月2日过火面积约为249 256.77m2；2019年4月3-5日MODIS数据云量太大无法使用，故无预测结果；2019年4月6日过火面积约为335 834.01m2。可以看到从2019年3月31日开始，火灾面积逐步扩大，最终过火面积为335 834.01m2，这与官方所公布的结果总体一致。

3.3 讨论

通过分析发现，结合多类型中高空间分辨率影像（如Landsat8 OLI、Sentinel-2和GF-1影像）、MODIS影像，并针对两种时空融合算法结果优化组合，可弥补单一时空融合方法或使用单一中等空间分辨率影像（如Landsat影像）和MODIS影像时空融合的不足，可获取更精确、更好的时间和空间分辨率预测影像。但对于实验存在的问题进行讨论：在应用时空融合时，我们发现在使用较差的预处理后的影像进行融合时，效果比较差，因此，时空融合的2种算法对数据的预处理工作有很高的要求，包括辐射定标、大气校正等，对于高分系列卫星还受到辐射校正和几何校正精度的影响，因此时空融合结果的精度一定程度上受影像预处理精度的影响；其次，该算法对于输入的影像数据有很高的要求，要求获取的影像必须是无云的纯净像素，MODIS影像的时间分辨率虽然比较好，但不可避免的也会出现研究区被云层覆盖的情况，如果是薄雾，可使用去云处理，但如果云层较厚，关于去云研究目前还未成熟。因此，在后续研究中，可以考虑加入雷达影像，利用雷达影像能够穿透云层，可以看到火灾的边界等信息的特点，探索雷达强度图像与光学数据的之间的关系，实现两者的融合，以期得到更好的应用。

4 结论

本文综合利用经典的STARFM算法与具有组合优势的基于地物内组分时相变化模型的地表反射率时空融合算法，联合使用多种空间分辨率更优（≤30m）的传感器影像（Landsat8 OLI、Sentinel-2、GF-1 WFV）与MODIS影像，以“时间最邻近及空间分辨率优先”为原则对传统单一中等空间分辨率影像预测周期（如Land⁃sat影像为16天）进行分段独立预测，并优化组合两种预测方法的预测结果，进而对四川凉山木里县3·30森林大火预测结果计算BAI和NBR指数，提取森林火灾区域动态变化信息。研究结果表明：本文所设计的时空融合策略可弥补单一时空融合方法或使用单一中等空间分辨率影像（如Landsat影像）和MODIS影像时空融合的不足，同时，使用NBR提取火灾变化信息效果更优。本文所提出的方法在森林火灾遥感动态监测场景中具有可行性，具有进一步深入研究的价值与意义。