基于遥感影像的林地变化检测方法

李 恒，臧 卓，唐 宪

(1.中南林业科技大学，长沙 410004；2.三亚市林业科学研究院，海南 三亚 572023)

引言

森林作为一种重要的生物资源，在改善环境、调节气候、维持生物多样性和促进经济可持续发展等方面有着不可替代的作用[1]。为了及时掌握森林资源的增减动态，自2017年开始，我国每年都要开展森林督查，暨森林资源管理“一张图”年度更新(原林地变更调查)工作，不断完善森林资源“一张图”数据库的建设。当前森林资源“一张图”更新工作的一般流程为首先判读区划变化图斑，再结合森林经营档案和外业实地调查核实图斑变化原因，最后更新森林资源数据库。目前生成变化图斑的主要方式是以高分辨率影像为底图，通过人工目视解译的方法来区划变化图斑[2]，在此过程中不仅需要耗费大量的人力物力，同时也对判读人员有着较高的专业和经验要求。因此，如何利用计算机图像处理技术和计算机视觉技术实现变化图斑的自动提取，以减少人力物力的投入，提高工作效率是当前需要解决的问题之一。

遥感图像变化检测是综合使用图像处理技术和计算机视觉技术，对比同一区域不同时期的影像，判断发生变化的区域及其空间分布的一门技术。最早将变化检测算法应用于光学遥感影像的研究可以追溯到上世纪70年代[3]，威斯米勒等[4]以德克萨斯州海岸的马塔哥尔达湾河口为研究区，分别用影像分类提取法、影像差值法、光谱—时间分析法、分层光谱—时间分析法四种方法对多期Landsat卫星数据进行处理来提取河口的变化区域。经过40多年的发展，国内外学者针对不同的类型提出了许多变化检测算法，大多数学者[5-6]将这些算法分为代数运算法、图像变换法、模型法、视觉分析法、影像分类提取法、GIS集成法和其他方法等7类。近年来，这些方法在森林资源遥感动态监测研究方面得到了广泛的应用。

代数运算法是指两期影像经过配准后，对前后期同一位置的像元值进行代数运算，然后得到差异图像，最后选取合适的阈值来提取变化区域。冯林艳等[7]利用归一化植被指数(NDVI)差值法进行林地变化检测，该方法在快速提取林地变化信息的同时具有较高的精度，但是无法有效地克服噪声和不同传感器之间辐射差异的影响。

对于多光谱图像来说，各波段之间在一定程度上存在相关性，造成图像信息冗余。因此，在变化检测过程中可以对多光谱图像进行线性变换，消除波段间的相关性，使图像信息集中到较少的分量。经过变化后的图形更有利于变化信息的提取[8]。常见的变换方法有主成分变换[9](PCA变换)和缨帽变换[10](K-T变换)。莫德林等[11]对比分析了基于主成分变换的主成分差异法、差异主成分法和多波段主成分分析法等在影像上的应用效果，最后得出主成分差异法能够全面的显示变化信息、差异主成分法提取变化信息细节的效果较好和多波段主成分分析能够检测出更符合实际的轮廓形状等结论。范应龙等[12]对遥感影像进行缨帽变换后得到研究区的亮度、绿度和湿度指数，然后根据植被覆盖变化与三种指数之间的关系，利用决策树分类来提取变化信息。虽然图像变换后能够消除波段间的相关性，实现数据降维和压缩，但是通常来说，图像变换的计算量大、过程复杂，并且变换后会丢失一部分影像的原始信息。

视觉分析法是指通过对影像进行假彩色合成，突出显示变化区域并结合目视解译的方法来勾画变化区域。该方法是目前森林资源年度更新中常用的方法[2]，对判读人员的专业知识和经验要求较高，且工作量较大。同时，假彩色合成法应用于资源卫星影像会有较大的色差，不利于目视解译[13]。

GIS集成法使用已知的地理信息数据辅助遥感影像来提取特征，通过特征对比来提取林地变化信息。赵珍珍等[14]提出了一种矢量数据与遥感影像相结合的变化检测算法，该算法利用历史矢量数据对遥感影像进行分割，获取前后期遥感影像的像斑，并根据矢量数据的属性信息对遥感影像进行分类，最后通过两期像斑分类结果的叠加分析提取变化信息。这种方法充分利用了GIS数据中丰富的历史知识，使得遥感影像分析变为基于对象的一种模式，有助于提高变化检测精度。但是在该方法中，变化检测结果的准确性依赖于地理信息数据的类型和准确性，且地理信息数据获取较为困难。

模型法是通过建立影像特征与实际林地变化之间的关系模型，定量地表达林地变化信息的方法。如龚鑫烨等[15]应用统计分布分析法和基于Zone模型的方法提取林地变化信息，试验结果表明基于Zone模型的方法的检测精度要高于统计分布分析法。使用生物模型进行林地变化信息提取的方法虽然能够简化流程并且有效地克服噪声和光照的影响，但是一个精度高且适用范围广的模型需要用大量的样本数据来训练调整。因此，相对于其他几种方法，模型法的鲁棒性较低。

影像分类提取法是对多期影像进行分类，通过比较分类结果来提取林地变化信息。如庞博等[16]利决策树分类法对两期TM影像进行分类，通过分类结果对比提取了1999年至2010年间伊春市林地面积的变化趋势；任冲等[17]以天水市为研究区，分别利用随机森林分类法和参数优化支持向量机分类法对多期Landsat影像进行了分类，然后逐期对比分类结果提取了1990年至2015年间森林资源的变化情况。影像分类提取法能够避免不同传感器和物候的影响，同时能够获得变化的方向信息。变化检测的精度主要取决于分类体系是否合理以及分类结果的精度。

在上述方法中，最常用的变化检测算法主要是图像代数运算法和影像分类提取法。因此，本文以广西省梧州市苍梧县木双镇为研究区，分别采用植被指数差值法、波段替换差值法、非监督分类提取法、监督分类提取法进行林地变化检测，其中非监督分类采用ISODATA分类法、监督分类采用最大似然分类法和支持向量机分类法。最后，对每种方法的结果进行了对比分析。

1 研究区和数据来源

1.1 研究区概况

木双镇位于广西省中东部，属梧州市苍梧县管辖，具体位置为23°41′14″N，111°34′14″ E，行政区总面积14 975 hm2，其中林地面积12 896 hm2，占比为86.1%；森林面积12 156 hm2。全镇森林覆盖率为81.2%，总蓄积量为737 632 m3。

1.2 数据源

选用木双镇2018年和2019年两个时期的高分二号卫星影像为基础数据。高分二号卫星发射于2014年08月19日，轨高度为631 km，轨道倾角为97.908 0°，其搭载了空间分辨率1 m的全色相机和空间分辨率4 m的多光谱相机，幅宽为45 km，重访周期为5 d。高分二号卫星在星下点空间分辨率可达到0.8 m，是我国首颗实现了亚米级分辨率的遥感卫星，它的成功发射标志着我国开始进入亚米级“高分时代”。

2 方法与步骤

2.1 影像预处理

影像的预处理主要包括辐射校正、几何校正和波段组合等过程。对于多时相遥感影像变化检测来说，前后两期影像的辐射差异和位置配准对检测结果的影响较大。

1)几何校正：本试验以前期(2018年)影像为参考，对后期(2019年)影像进行校准。经校准后，将两幅影像的精度控制在0.5个像素以内[18]。

2)辐射校正：由于天气、传感器角度等因素的影响，前后期遥感影像存在一定的辐射差异。本试验采用直方图匹配法，以前期影像为参考，对后期影像的灰度分布进行变换，使两幅影像的灰度分布相近，缩小影像间的辐射差异。

3)波段组合：许多学者研究发现，可见光—近红外波段在林业遥感研究中具有重要意义[19]，因此，本文选择的波段组合为近红外—红—绿波段。图1为经过几何校正和辐射校正后的双时相遥感影像，在重新组合的遥感图像中，红色区域为植被覆盖区域，绿色区域为无植被覆盖区域，深蓝色区域为水域。

图1 研究区前后期遥感影像

2.2 研究方法

2.2.1 植被指数法

植被指数是植被多光谱遥感数据中的多个波段之间经过一定的数学运算得到的数值。该数值能够定性和定量地表示植被覆盖度、生长状况和生物量的多少。植被指数是监测地面植物状态行之有效的方法之一。在众多的植被指数中最常用，应用范围最广的是归一化植被指数[20](NDVI)，它不但对地面植被敏感度高，而且能够消除大气和土壤的影响。

NDVI=(NIR-R)/(NIR+R)

(1)

式中，NIR红外波段，R是红光波段。

试验首先计算分别得到两期影像的NDVI图，然后用前期NDVI减去后期NDVI得到差值影像，最后根据经验阈值法确定一个合适的阈值来区分变化区域和未变化区域。

2.2.2 波段替换法

在遥感动态监测方法中，图像差值法是最快捷简单的方法，但是由于天气、传感器角度等因素的影响，使得图像中未发生变化的区域被误检出来。在本实验中，为了消除未变化区域的影响，用前一时相的R波段和B波段以及后一时相的G波段重新组合作为新的后期影像。然后用加权平均法将前期影像和新合成的后期影像转为灰度图像。

IM=0.298 9×R+0.578×G+0.114×B

(2)

式中，IM是计算所得的灰度图像；R，G，B是影像三个通道的DN值。

对两期灰度图像进行减法运算得到差值影像，最后选择合适的阈值提取变化区域。

2.2.3 监督分类与非监督分类提取算法

利用监督分类和非监督分类的算法，对两个时相的遥感影像进行分类，提取两期遥感影像的林地范围，最后将两期提取结果作差值运算，提取林地变化范围。该方法不仅能够提取变化区域，还能够提供变化矩阵信息。这种方法的检测精度主要受分类精度的影响。本试验分别使用ISODATA分类法、最大似然分类法和支持向量机分类法对遥感影像进行分类并评价了每种方法的精度。其中，ISODATA分类法属于非监督分类，最大似然分类法和支持向量机分类法属于监督分类。得到每种分类方法的结果后，通过前后期分类结果比较来提取变化区域。

2.2.4 精度评价

通过构建变化检测结果与地表真实变化数据之间的误差矩阵[21]，计算总体精度、漏检率和错检率等评价指标，用以评价变化检测结果的精度。其中，地表真实变化数据由人工结合历史高分辨率影像目视判读区划生成。

3 结果与分析

3.1 变化检测结果

本试验的目的主要是检测林地上发生变化的区域，因此所有的检测结果都用研究区的林地范围进行掩膜处理。

3.1.1 植被指数法

两期NDVI影像(图2)相减得到植被指数差值影像后，根据经验阈值法选取的阈值逐像元判断，大于阈值的像元认定为变化区域，其值赋值为1，其余未变化像元赋值为0，最后得到二值化的检测结果(图3)。图中白色区域为林地发生变化的区域。

图2 植被指数图

图3 植被指数法变化检测图

经过植被指数法的检测，木双镇在2018—2019年间共有1 483 hm2的林地发生了变化，主要分布于研究区的西部和东北部。

3.1.2 波段替换法

用前期影像的R波段、B波段分量替换后期影像对应分量后的效果如图4所示。经过替换后，后期影像中未变化区域的灰度值与前期影像更为接近，同时变化区域的灰度值依然明显区别于前期影像，使得变化区域与未变化区域更容易区分开。

(a)2018年遥感影像

前期影像与经过替换的后期影像进行减法运算，然后使用经验阈值法确定一个合适的阈值来逐像元判断差值影像，大于设定阈值即为变化区域，将其赋值为1，其余像元赋值为0，最后的变化检测结果如图5所示，白色为林地发生变化的区域。

图5 波段替换法变化检测图

经过植被指数法的检测，木双镇在2018—2019年间共有513 hm2林地发生了变化，主要分布于研究区的西南部和东北部。

3.1.3 监督分类与非监督分类的提取结果

本试验的目的主要是提取林地上发生变化的区域，因此选择将研究区的地物分为植被，非植被和水体三种类型。对于ISODATA分类法，得到初步分类结果后根据影像确定每个类别的归属。对于监督分类，首先在研究区范围内均匀地区划训练样本，本试验中森林样本区划37块、非森林样本36块、水体样本28块，且所有样本在前后期的影像中均为同一种地类。最终各方法的分类结果如图6所示。

(a)ISODATA分类

在基于影像分类提取变化信息的算法中，遥感影像分类结果的精度直接影像变化检测结果的准确程度，因此在本次试验中采用混淆矩阵对上述遥感影像分类结果进行精度评价。用于评价分类精度的验证样本数据由人工结合历史高分辨率影像目视判读区划，其中，森林36块、非森林45块、水体28块，共计109块，并且各地类在前后期均为一致。具体评价结果如表1所示。

表1 分类结果精度评价%年度ISODATA分类法最大似然分类法支持向量机分类法201878.1699.1399.54201984.5099.6199.74

结果表明，监督分类的总体精度要优于非监督分类，支持向量机分类法总体精度稍高于最大似然分类法。此外，在本次分类由于只区分植被、非植被和水体三种类型，且三种类型之间差异很明显，因此三种分类方法的精度都比较高。

同样地，得到两期分类图之后进行代数计算，如果前后两期分类结果相同，那么相减之后结果为0，因此差值影像中像元值不为0的区域即为变化区域。三种分类方法得到的变化检测结果如图7所示，绿色区域代表林地正向变化，红色区域代表林地负向变化。

(a)ISODATA分类法结果

3.2 精度评价

为了评价每种变化检测方法的精度，在研究区林地范围内通过目视解译的方法人工区划变化区域，并将其作为研究区内林地实际发生变化的区域，如图7所示，图中白色为林地发生变化的区域，黑色为林地未变化区域。

然后通过构建混淆矩阵的方法来评价每种方法的精度，主要的评价指标有总体精度、漏检率、错检率和Kappa系数。总体精度是指检测结果与实际数据相符的样本占总体的比例；漏检率是指实际发生变化而检测结果为非变化的样本数占实际变化总样本数的比例；错检率是指实际未发生变化而检测结果为变化的样本数占检测结果为变化总样本数的比例。每种检测方法的混淆矩阵如表2所示，评价结果见表3。

表2 变化检测混淆矩阵变化检测结果实际变化像元实际未变化像元合计NDVI差值法变化像元1 434 3762 273 7743 708 150NDVI差值法未变化像元1 041 06930 482 13031 523 199波段替换法变化像元984 606298 0101 282 616波段替换法未变化像元1 490 83932 457 89433 948 733ISODATA分类法变化像元1 764 43512 917 44114 681 876

表3 变化检测精度统计分类法总体精度/%漏检率/%错检率/%Kappa系数NDVI差值法90.5942.0661.320.41波段替换法94.9260.2323.230.50ISODATA分类法61.3228.7287.980.10最大似然分类法92.4920.3452.080.56支持向量机分类法93.4120.7047.960.60

由表3可知，ISODATA分类法的精度最低，为0.1；其余四种方法的Kappa系数均为0.4以上，支持向量机分类法的最高，为0.6。在漏检率检查中，波段替换法的误差最高，其次为NDVI差值法，其原因可能是根据阈值提取变化区域，由于前后期差异较小，导致部分区域被漏检；ISODATA分类法的错检率最高，原因是ISODATA分类法的精度低，且错分类现象严重。综合来看，支持向量机分类法的提取结果最优。

4 结论与讨论

4.1 结论

以梧州市木双镇2018年和2019年两期高分二号影像为基础数据，利用植被指数差值法、波段替换法、ISODATA分类法、最大似然分类法、持向量机分类法来提取林地变化信息，并对每种方法进行了精度评价。经过综合对比分析每种方法的总体精度、漏检率和错检率等评价指标，得出以下结论：

1)基于监督分类的变化检测算法提取林地变化结果的精度比植被指数差值法，波段替换法和基于非监督分类的变化检测算法要高。

2)对于基于影像分类的变化检测算法来说，分类结果精度对变化检测结果影响很大，分类精度越高，变化检测的精度越高。

4.2 讨论

本试验中变化检测方法都是基于像素来进行的，因此在变化检测结果中有很多噪声，降低了检测精度。在实际的地理环境中，每种地理实体在空间上都具有相关性，因此，可以将地理实体认为是一个对象。下一步可通过分别对两期影像进行面向对象分类，以对象而非像元为单位来进行变化检测的方法来改善此类情况[22-23]。同时，近年来越来越多的学者致力于研究深度学习理论在遥感影像变化检测中的应用[24-27]，并且取得了较高的精度，所以未来基于学习型算法的林地变化检测方法，也值得进一步去探索研究的。