李淑君，郑柯，唐娉，霍连志，袁媛
1.中国科学院空天信息创新研究院,北京 100094;
2.中国科学院大学,北京 100049;
3.南京邮电大学,南京 210042
森林火灾是一种突发性强、破坏性高、处置困难的自然灾害（Cao等，2020）。近年来，受全球气候变化与人类活动的影响，森林防火行业面临着愈发严峻的考验，对森林火灾识别技术提出了新的需求与挑战（魏书精等，2020；
王振师等，2019；
孙红斌等，2018）。火烧迹地是描述森林燃烧最重要的信息之一，是经火灾烧毁后尚未长成新林的土地，能够提供森林火灾发生时间、频度、位置、面积以及空间分布等重要信息（Ruiz等，2014），是火灾形成机制不可或缺的研究对象（Chuvieco等，2016），也是影响森林资源保护、植被恢复、碳排放等至关重要的变量（Chuvieco等，2018），以及碳循环扰动和全球变化研究的重要参数（孙桂芬等，2019）。遥感技术具有观测覆盖范围大、可动态观测等优势，使其成为最为有效的森林监测手段，如何利用时间序列遥感图像提取森林变化信息是当前研究的热点和难点问题（杨辰，2013）。
根据数学统计方法划分，传统的时间序列变化检测算法可分为6大类，包括阈值法、差分法、分段法、轨迹分类法、统计边界法和回归法（Zhu，2017）。阈值法采用了一个预定的阈值确定时间序列中的森林，当与阈值存在显著偏差时检测到变化，Hilker等（2009）基于Landsat和Modis数据的NDVI指数利用融合模型对森林实现变化检测，Huang等（2009）利用森林综合指数完成率对美国西部森林区域动态变化的评定（IFZ）；
差分法通过比较不同时间获取的影像，并且将变化很大的差异定义为变化区域，Bolton等（2015）使用差分法对火灾后的残留结构和森林恢复进行分析；
分段法将完整的历史时间序列数据根据残差或角度标准分为一系列直线段，检测突变和渐变信息（Chance等，2016）；
轨迹分类方法从时间序列中提取某种变化（训练目的）信息，用该信息对图像中的每个时间序列进行分类，Kennedy等（2007）基于轨迹对森林动态特征进行自动匹配实现变化检测；
统计边界法将明显偏离边界的检测为变化，Zhu等（2016）利用Landsat数据对绿度趋势分析土地覆盖变化；
回归法假定预测变量和响应变量之间存在线性关系，使用回归获得结果（Markham等，2012）。此外，还有一些其他常用的经典时间序列变化检测算法，如干扰和趋势监测（LandTrendr）方法（Kennedy等，2010）、植被追踪算法（VCT）（Huang等，2010）、bfast算法（Verbesselt等，2012）、连续变化检测和分类算法（CCDC）（Zhu和Woodcock，2014）等。
近年来，深度学习得以迅速发展，并越来越多地用于时间序列分析（Längkvist等，2014），在时间序列森林变化监测方面的应用也越来越广泛，Reddy和Prasad（2018）将LSTM用于海岛植被动态监测，预测岛上植被未来时间序列以获取植被变化信息；
Kong等（2018）基于LSTM模型对长时间序列森林火灾进行监测，检测率高于传统的bfast算法；
Yuan等（2020）进一步在LSTM模型中加入Attention机制并用于森林火灾监测，检测精度相比LSTM有很大提升。然而，LSTM模型主要是对图像中单个像元的时间序列进行分析，没有考虑像素的空间邻域信息，而空间邻域的相关性对于遥感图像分析具有重要意义（Huo等，2015）。
为在时间序列预测中同时考虑空间信息，Shi等（2015）提出了ConvLSTM，将原来LSTM中的全连接改进为卷积连接，可以较好地刻画时间序列的时空特征，并用于降雨预测。Stacked ConvLSTM堆叠多个ConvLSTM增加网络的表达性能，提升了对时间序列数据时空结构信息的提取能力，因此，Kim等（2017）用两层ConvLSTM基于雷达数据进行降雨的多步预测，Zhao等（2019）基于三层的ConvLSTM对超声波时间序列影像进行单幅影像预测。
由于Stacked ConvLSTM尚未用于森林火烧迹地检测，其在森林火烧迹地检测的有效性和效果有待验证，本文提出将时空检测方法Stacked ConvLSTM用于时间序列森林火烧迹地的检测，实现端到端提取森林火烧迹地信息。由于火烧迹地在空间分布上具有连续性（武晋雯等，2020），因而在时间序列变化检测中加入邻域空间信息具有重要意义。一方面，能够在一定程度上降低因数据缺失（云、阴影、地形遮挡等所致）造成的影响（Khan等，2017）；
另一方面，相比现有方法大多采用先分类再后处理抑制零散虚警信息的两阶段策略（张兆明等，2020），端到端数据驱动的方法能够避免后处理方法选择的主观性。
2.1 研究区域及数据
两个研究区域分布如图1所示，研究区域Ⅰ是黑龙江沾河林业局伊南河林场，地理位置在48°39＇N，128°14＇E附近区域，于2009年4月27日发生草甸森林火灾，过火总面积1557.562 km2，其中有林地1072.746 km2。
图1 研究区域（研究区域Ⅰ：伊南河林场火灾区域，研究区域Ⅱ：北大河林场火灾区域）Fig.1 Study area（study areaⅠ：Yinanhe Forest Farm fire area，study areaⅡ：Beidahe Forest Farm fire area）
研究区域Ⅱ在内蒙古自治区毕拉河林业局北大河林场，地理位置在49°31＇N，123°06＇E附近，2017年5月2日发生森林火灾，过火面积115 km2，有林地占60%，受害森林面积达82.816 km2。
选择这两个研究区域的原因有以下几点：（1）两个研究区域的火灾面积较大，在遥感影像上有明显的过火区域轮廓，便于研究实验的进行；
（2）两个区域的历史时间序列都比较平稳，除了已知的火灾外没有受到其他明显的森林干扰发生，便于模型算法的验证；
（3）两个研究区分别包含草甸森林和森林区域，为在不同的森林覆盖类型验证模型的有效性提供了条件。
本文采用的研究数据是中分辨率成像仪（MODIS）陆地产品MOD13Q1。MOD13Q1空间分辨率250 m，是16 d合成产品，每年有23幅影像。分别获取了研究区域Ⅰ的2001年—2009年以及研究区域Ⅱ的2001年—2017年间的时间序列数据。实验剔除了序列中受云、雪影响大的春冬季影像，采用一年中季相变化比较明显的14幅影像作为研究的时间窗口（DOY97—DOY305）。
MODIS数据的MOD13Q1产品中增强型植被指数（EVI）在植被监测中具有时相多、覆盖面广、且不易产生过饱和等优势（Huete等，2002；
王正兴等，2003），被广泛应用于大区域尺度的植被物候变化监测等研究领域（Peng等，2017；
Zhang等，2003）。因此本文选用增强型植被指数EVI来研究其时间序列上的动态变化，EVI的计算公式可表示为
式中，ρR、ρB、ρNIR分别表示MODIS数据的红、蓝和近红外波段的反射率。
选择MODIS数据的原因主要有两点：（1）MODIS数据的时间跨度长，可以提供足够长的时间序列；
（2）MODIS 16 d合成产品，时间间隔均匀，方便用于时间序列的统计建模。
2.2 Stacked ConvLSTM
时间序列数据预测是指学习过去的时间序列并预测未来的变化。传统的神经网络无法解决随时间轴变化的问题，随之诞生了RNN（Recurrent Neural Network）（Jordan，1997）；
但由于经典的RNN在提取较长时间序列信息的表现较差，提取的时序信息有限，Hochreiter和Schmidhuber（1997）提出了LSTM网络模型，在经典RNN中加入门（gates）的结构来选择性的增加和删除过去的时序信息，即加入输入门、输出门和遗忘门，分别控制本单元（一个LSTM结构为一个基本单元）数据的输入、输出以及上一步单元的输出的信息增减。LSTM公式表示如下（Hochreiter和Schmidhuber，1997）：
式中，I、F、O分别表示输入门、遗忘门和输出门，C和H分别表示细胞状态（经过门控输出的信息）和隐藏状态（每个时间点的输出值），W表示对应数据的权重，X表示输入数据，b表示偏置值，σ表示激活函数，o表示哈达玛积，下标t表示t时刻；
下标i、f、out分别表示3个控制门对应的权重和偏置值，下标c表示细胞状态C对应的权重和偏置；
W的下标x表示对应输入数据X的权重，h表示对应隐藏状态H的权重。
ConvLSTM是在LSTM基础上提出来的一种LSTM变体，将LSTM的输入层到隐藏层和隐藏层到隐藏层之间的全连接状态替换为卷积连接，对LSTM无法充分利用空间信息进行了改进。LSTM在处理图像数据时需要将图像数据转为一维向量，无法处理原图像数据的空间结构信息。相比LSTM模型，ConvLSTM能够更好地提取时间序列图像中的时空结构信息。ConvLSTM模型公式表示如下：
式（6）—（9）中字母表达含义和式（2）—（5）
相同，输入变量的全连接换成了卷积操作，*表示卷积操作，由图2的ConvLSTM内部结构可以看出，输入门、输出门和遗忘门对于输入和隐藏层都进行了卷积操作；
其中WcioCt-1、WcfoCt-1和WcoutoCt-1表示输入门、遗忘门和输出门与上一细胞状态的peephole（Gers和Schmidhuber，2000）连接。如图2所示，peephole连接在每个门都加入细胞状态的信息，由于网络单元可能会有门状态为0的情况，导致缺乏重要信息，添加peephole操作可以改进这个缺点。
图2 ConvLSTM内部结构Fig.2 Internal structure of ConvLSTM
本文的Stacked ConvLSTM网络结构如图3所示，Stacked ConvLSTM在ConvLSTM基础上进行改进，由多个ConvLSTM堆叠而成。假设网络输入时间序列数据为(X1，X2，…，Xn)，每次输入一个Xi经过n次数据输入提取输入的时序特征，经过k层ConvLSTM网络单元前向传递网络学习到的时空信息，最后一层三维卷积层输出最后结果(Xn+1，Xn+2，…，Xn+m)，将输出结果与目标序列计算损失函数，并反向传递进行网络训练。
图3 堆叠ConvLSTMFig.3 Stacked ConvLSTM
式中，Xi表示第i时间点目标序列值，表示第i时间点网络模型的预测值。
2.3 火烧迹地检测精度评价指标
混淆矩阵也称作误差矩阵，是表示分类精度的一个n×n矩阵（Liu等，2007）。在本文中火烧迹地和非火烧迹地属于二分类问题，混淆矩阵如表1所示。
表1 火烧迹地检测混淆矩阵Table 1 Confusion matrix of fire detection classification
基于混淆矩阵，进一步计算以下常用的结果评价指标，包括：精确度（P）、召回率（R）、准确度（Acc）、F1-score值（F1）。以上4个评估指标公式如下：
式中，TP表示正确分为火烧迹地的像素个数，FP表示错误分为火烧迹地的像素个数，TN表示正确分为非火烧迹地的像素个数，FN表示错误分为非火烧迹地的像素个数；
P＇是分为火烧迹地的总像元个数，N＇是分为非火烧迹地的总像元个数；
P是真实火烧迹地像元个数，N是真实非火烧迹地像元个数。
3.1 数据处理
基于深度学习框架Keras完成了Stacked ConvLSTM的构建（Chollet等，2015），实验设备环境是NVIDIA GeForce GTX1080 GPU。研究区域Ⅰ的MODIS时间序列数据集图像像素区域大小为512×512，研究区域Ⅱ图像像素区域大小为400×400。在这两个研究区域，MODIS时间序列中观测到森林火烧迹地信息的变化时间点分别为2009年和2017年的第126天，对应本研究的时间序列中第3个时间点影像。如图4是在两个研究区分别取像素区域大小为10×10的平均值在时间轴的变化，发生火灾的时间点在最后一年（图4中红色显示），是需要进行预测的年份。
图4 研究区域EVI在时间轴上变化Fig.4 EVI varies along the time axis in the study area
本文使用滑动窗口法（Cheng等，2016）将研究所用MODIS数据划分为监督训练的输入输出集。在时间上，输入序列长度为1年时序数据14幅影像，输出序列长度为后1年时序数据14幅像，滑动的步长为一年即14幅影像。例如，2001年数据作为输入，2002年数据作为输出，下一个滑动窗口的输入为2002年数据，输出为2003年数据，依次类推至时序最后一年。最后一次滑动作为测试数据，研究区域Ⅰ输入2008年数据预测2009年，研究区域Ⅱ输入2016年数据预测2017年。在空间上，将研究区域Ⅰ的像素区域大小512×512切分成像素区域大小64×64的子块，研究区域Ⅱ的像素区域大小400×400分成像素区域大小50×50的子块。对处理后的数据集，除了最后一年的数据作为测试数据，剩下的数据按照8∶2划分为训练集和验证集。
3.2 时序预测与火烧迹地检测
本文对Stacked ConvLSTM模型中的超参数设置了不同的值进行实验，得到的最优参数分别为卷积核大小3×3，Batch Size大小设置为16，采用的优化器为RMSprop（Zhong等，2012），学习率大小设置为0.001，梯度衰减值为0.9，网络的循环训练次数为100次。
Stacked ConvLSTM的Stacked Layers及每层的单元个数对于网络模型的学习能力影响很大，因此对不同的网络结构进行测试。网络所测试的层数以及单元数是在Stacked ConvLSTM的应用以及实验中得出的经验值，根据预测精度取在本文数据表现最佳的网络层数和单元数，作为的森林火烧迹地检测网络结构。
研究区域Ⅰ测试结果如表2所示，由平均绝对误差MAE（Mean Absolute Error）和均方误差MSE（Mean Squared Error）指标可以看出，研究区域Ⅰ的预测效果最好的是32-32-32-batch16结构；
研究区域Ⅱ的结构如表3所示，效果最好的结构是64-32-batch16。
表2 研究区域Ⅰ预测精度Table 2 Predict precision of study areaⅠ
表3 研究区域Ⅱ预测精度Table 3 Predict precision of study areaⅡ
将两个研究区域的Stacked ConvLSTM最佳网络结构用于MODIS数据的森林火烧迹地检测。取预测的时间序列以及验证的目标序列中像素区域大小10×10求平均值作图结果如图5所示，可以直观地看出，未变化区域的预测值和目标序列值很接近，而变化区域在第3个时间点（MODIS时间序列的第126天）变化之后一段时间EVI值下降很大，与预测值相距悬殊，植被在后半年有所恢复，EVI值回升至和预测值相近。
图5 未变化和变化区域时序拟合Fig.5 Unchanged and changed area fitting
对测试序列中每一个像元采用经验阈值的方法确定是否为变化像元（Kong等，2018）。根据检测精度分别调整两个研究区域的阈值，获得最佳阈值为-0.05。如图6所示，阈值分割线表示预测序列每个像素的EVI值减去0.05。并提取目标序列首次低于阈值分割线的像元作为变化像元，在测试时间序列的第3个时间点检测到变化。
图6 变化区域阈值Fig.6 Threshold of changed area
本文利用Stacked ConvLSTM、Stacked LSTM以及bfast算法在两个区域的MODIS时间序列中提取森林火烧迹地，并将火烧迹地检测结果与ESA发布的Fire_CCI 5.1火烧迹地产品进行了对比分析。地面真值数据分别由2009年5月23日Landsat-5 TM影像和2017年5月25日MODIS影像目视解译生成。
火烧迹地检测结果如图7所示，从目视效果来看，在研究区域Ⅰ，Stacked ConvLSTM检测的结果比Stacked LSTM和bfast算法错误检测点少，并且在空间分布也保持较高连续性；
在研究区域Ⅱ，Stacked ConvLSTM检测到了较完整的火烧迹地区域。在研究区域Ⅰ，Stacked ConvLSTM检测的结果相对Stacked LSTM以及bfast算法更少细碎的像素，错误检测点少并且在空间分布也保持相对高的连续性，检测的火烧迹地精度相对于另外两种方法也更高，Fire_CCI 5.1的错检区域比Stacked ConvLSTM更大。在研究区域Ⅱ，Stacked ConvLSTM检测到了较完整的火烧迹地区域，Stacked LSTM以及bfast算法的错误检测像元相对更多，并且在空间连续性上表现更差，Fire_CCI 5.1相比Stacked ConvLSTM漏检像元更多。
图7 森林火烧迹地检测结果图Fig.7 Forest burned areas detection results
由表4和表5可以看出，在两个研究区域Stacked ConvLSTM的检测精度相比Stacked LSTM和bfast算法更好。在研究区域Ⅰ，Stacked ConvLSTM的精确度比Stacked LSTM和bfast算法分别高出0.120和0.405，并且召回率、准确度和F1-score也更高，Fire_CCI 5.1召回率虽更高，由于错检区域较大，其他精度指标远低于Stacked ConvLSTM；
在研究区域Ⅱ，Stacked ConvLSTM精确度达0.924，召回率、准确度和F1-score相比Stacked LSTM和bfast算法以及Fire_CCI 5.1更高。
表4 研究区域Ⅰ火烧迹地检测精度表Table 4 Precision table of burned areas detection in study areaⅠ
表5 研究区域Ⅱ火烧迹地检测精度表Table 5 Precision table of burned areas detection in study areaⅡ
将两个区域的精度表绘制成柱状统计图，见图8所示，可以更直观地看出，在研究区域Ⅰ，Stacked ConvLSTM比Stacked LSTM的检测精度要好，并远高于bfast算法；
在研究区域Ⅱ，Stacked ConvLSTM的召回率以及F1-score比另两种方法高很多。总体来看，Stacked ConvLSTM的检测表现比Stacked LSTM和bfast算法表现更好。
图8 森林火烧迹地检测精度柱状统计图Fig.8 Precision of forest burned areas column chart
Stacked ConvLSTM在两个区域检测的火灾发生的时间点分别在MODIS时间序列中的2009年和2017年的第126天，而实际发生的火灾时间点分别在2009年的4月27日和2017年的5月2日，即第116天和第122天。以上3种方法在两个区域的检测时间点如表6和表7所示，Stacked ConvLSTM和Stacked LSTM比实际发生的时间点分别迟10天和4天，bfast算法的检测时间点相对另两种方法要迟一步。检测时间点的推迟和MOD13Q1数据的时间分辨率有关，MOD13Q1是16天合成产品，发生的火灾点和获取的数据点正好相差了一个时间周期，导致了检测时间点的延迟。
表6 研究区域Ⅰ变化时间点Table 6 Change detected time of study areaⅠ
表7 研究区域Ⅱ变化时间点Table 7 Change detected time of study areaⅡ
本文提出将Stacked ConvLSTM用于时间序列森林火烧迹地检测，利用Stacked ConvLSTM学习时空特征的优势，在保持结果具有较好空间连续性的基础上避免了具有主观性的后处理操作，实现端到端提取森林火烧迹地信息，提升了森林火烧迹地的提取精度。基于两个发生过特大火灾区域的MODIS时间序列数据，定性和定量分析Stacked ConvLSTM、Stacked LSTM以及经典的bfast算法火烧迹地提取结果，并将3种方法的火烧迹地提取结果与Fire_CCI 5.1火烧迹地产品进行了对比。研究结果表明，在检测精度和空间分布的连续性上，Stacked ConvLSTM检测的火烧迹地结果都有很大提高；
此外，在火灾发生时间点检测效果上，由于MODIS时间分辨率（MOD13Q1数据为16 d合成数据）为16 d，Stacked ConvLSTM和Stacked LSTM均落后真实火灾时间点一步，但均比bfast算法的在时间上更准确一步。
Stacked ConvLSTM在森林火烧迹地检测存在一些不足之处。首先，本文基于单个EVI指数开展Stacked ConvLSTM在森林火烧迹地检测的有效性研究，更多植被相关指数及其组合的检测效果值得进一步研究；
其次，确定火烧迹地区域的方法是基于经验阈值，有待研究自动的阈值计算方法；
再次，由于MODIS数据的分辨率属于中低分辨率，Stacked ConvLSTM模型的卷积核设置太大反而会模糊空间信息，因此提取的空间信息有限，在更高空间分辨率和时间分辨率的影像数据上，Stacked ConvLSTM模型的时空预测表现可能会更好；
最后，网络模型训练中是固定时间长度的输入输出，只能学习有限的时空序列信息，这大大限制了网络在时空预测上的表现。因此，在后续的研究中，考虑在Stacked ConvLSTM基础上加入一些可以输入输出不定长时序数据和更高效的措施如Encoder-Decoder结构以及Attention机制等，提高网络在时空序列上的预测效果。
志谢此次实验过程中，在袁媛老师的帮助下获取了本文的实验数据，在此表示由衷感谢！
猜你喜欢 迹地火烧火灾 封二、封三：寒葱岭战迹地 （伏击战） 新考兰台内外(2021年33期)2021-12-17火烧毒品南方周末(2019-06-27)2019-06-27我国科学家在国际上首发30 m分辨率全球火烧迹地数据产品军民两用技术与产品(2018年11期)2018-07-16我国科学家在国际上首发30米分辨率全球火烧迹地数据产品河南科技(2018年16期)2018-07-13驴肉火烧小学生作文(低年级适用)(2018年12期)2018-03-23起源滇池(2017年5期)2017-05-19不怕火烧的手帕莫愁(2017年15期)2017-05-11掌握火灾逃生知识文理导航·趣味课堂(2016年6期)2016-09-09离奇的火灾故事作文·高年级(2009年7期)2009-08-20