火灾、地震等自然灾害和人为因素干扰常常导致森林面积减少, 类型和健康质量变化。常规的地面调查手段, 难以监测大区域、不易达到区域的森林景观变化。遥感数据具有观测区域大、可以重复访问等优势, 在快速识别、监测森林景观类型变化方面发挥着重要作用, 对于生态恢复、景观重建等具有重要意义。

利用遥感数据进行变化检测，一直是遥感界重究; 覃先林等(2006)利用MODIS数据，分别采用红光-近红外法、共生纹理矩阵法和基于相似度的变化检测方法对我国东北林区的森林覆盖变化进行了研究，探明了大尺度森林火灾干扰引起的变化要的研究热点(Jensen et al., 1993; Ridd et al., 1998; Lunetta et al., 2000)。迄今为止，国内外专家学者已经开发出了众多的遥感变化检测识别算法，包括图像差值、缨帽变换、分类后比较、光谱变化向量分析等(Jensen，2005)，专家系统和神经网络也已经应用到变化检测领域(Lu et al., 2004)。在森林景观变化研究方面，杨存建等(2000)利用植被指数插值法提取了云南省泸水县的森林植被变化信息; Sunil等(2004)用分类后比较方法进行地表覆盖分类，获得了不同时期(1940s, 1960s, 1990s)爱荷华州国家公园管理区的森林、草地、农田等景观变化信息; Sean等(2005)对Landsat数据进行缨帽变换来监测植被的变化，获得了区域尺度上森林采伐干扰引起的变化; Jin等(2005)对缨帽变化湿度和归一化水分指数探测森林干扰并进行了比较分析; Okhandiar等(2008)利用印度的IRS卫星数据，采用邻近相关图像分析技术进行了森林景观变化研

邻近相关图像(neighborhood correlation image，NCI)分析技术是一种新的变化检测方法(Im et al., 2005; Im，2006)，利用2个时相遥感图像中像元及其周围像元之间的相关性进行变化探测。这项技术由Im等(2005)提出，最初用于高分辨率图像城市信息变化的检测，检测精度较高。本文主要探讨基于邻近相关图像和决策树分类的森林景观变化检测方法，利用2个时相遥感影像探测干扰条件下森林景观的变化，包括森林景观类型和范围的变化，以期为受灾林区、特别是自然保护区的生态恢复与重建提供科学依据。

1 研究区与数据 1.1 研究区概况

以5·12汶川特大地震中森林遭受严重破坏的龙溪—虹口国家级自然保护区为研究区。该保护区位于四川省都江堰市北部，地理位置31°04'—31°22'N，103°32'—103°43' E，总面积达310 km²; 属四川盆地亚热带湿润气候区，年平均气温10 ℃，雨量丰沛，年降水量约1 600~1 900 mm，有“华西雨屏”之称。整体海拔由南向北递增，相对高差大，从低到高呈现亚热带季风气候到寒温带的气候变化，植被也呈现明显的垂直分布，由低到高依次为常绿阔叶林、常绿落叶阔叶混交林、常绿针叶与落叶阔叶混交林、亚高山针叶林、高山灌木林、高山草甸、高山流石滩植被。

1.2 数据获取与数据处理 1.2.1 数据源

选取研究区域地震前后2个时相的Landsat5 TM多光谱卫星遥感数据，如图 1所示。地震前数据的过境时间为2007年9月18日，震后数据的过境时间为2008年7月18日。采用蓝波段(TM1)、绿波段(TM2)、红波段(TM3)、近红外波段(TM4) 和2个短波红外波段(TM5，TM7) 进行变化检测研究。

选择相同传感器的遥感数据保证了遥感观测条件的基本一致，即卫星在同一地方、同一方向通过同一地点，这样更利于图像的对比。获取时间均为夏季，减少了植被在不同物候条件下的光谱响应差异。以当地保护区管理局提供的植被图与2010年8月进行的野外实地调查灾害点的类型、位置数据作为研究的辅助参考数据。由于当地的植被有垂直地带性分布的特点，获取了龙溪—虹口国家级自然保护区范围的30 m空间分辨率ASTER全球数字高程数据(GDEM)数据作为辅助参考数据。

1.2.2 数据预处理

数据预处理包括辐射定标、大气校正和几何校正。首先对2个时相影像进行辐射定标，将其DN值转化为辐射亮度值，再使用ENVI软件中的Flaash模型对其进行大气校正，得到地面真实反射率值; 尔后对2个时相影像进行几何纠正，配准误差小于0.5个像素; 对于影像区域中有云、阴影的地方进行掩膜处理。

2 原理与方法 2.1 邻近相关图像分析原理

该方法是利用2个时相遥感图像中一定区域的像元值来计算该区域2个时相影像光谱值的相关系数，并采用最小二乘法得到该区域2个时相像元值回归方程的斜率和截距，并将值赋给该区域的中心像元，最终得到整幅邻近相关图像:相关系数图、斜率图、截距图。

邻近相关图像能够反映出2个时相影像在相同区域内光谱值变化的数量和方向。如果在2个时相图像中的特定区域内像元的光谱值发生变化，则该区域内的相关系数就会降低，变化越大相关系数越低，此时斜率和截距的增加或减少可以反映出光谱变化的大小和方向。假设没有发生任何变化，特定区域内同名点的2个时相亮度值点对(x，y)应该在其二维散点图中处于y=x的直线上，此时斜率为1，而截距为0;但如果发生了变化，点对就远离y=x这条直线，斜率不为1，截距不为0，而且变化越大，斜率越远离1，截距越远离0，斜率和截距远离的方向代表着光谱变化的方向。所以，斜率图像和截距图像可以辅助相关系数图像来进行变化检测。利用相关系数、斜率、截距，结合专家系统分类方法可以提取详细的变化类型信息。

2.2 基于邻近相关图像和决策树的变化检测方法

基于邻近相关图像和决策树分类的变化检测流程如图 2所示。该流程主要包括5个步骤: 1) 数据预处理，包括辐射定标、大气校正、几何纠正; 2) 创建由不同窗口大小生成的相关系数图像、斜率图像、截距图像; 3) 将2个时相图像分别与不同窗口的邻近相关图像叠加，结合辅助数据和调查数据，建立训练点和检验点; 4) 利用训练点数据建立分类规则，并生成变化检测图; 5) 利用检验点数据对变化检测结果进行精度评价。在该流程中，每一步都至关重要，因为任何一步的误差都会影响到最终的变化检测精度。

2.2.1 创建邻近相关图像

本研究利用IDL(interactive data language)编程，并选取不同窗口大小(3×3，5×5，7×7，9×9) 创建邻近相关图像(NCIs)，以探索变化检测模型的最佳配置窗口大小。将最高变化检测精度所对应的窗口大小认为是该区域的最佳窗口。相关系数、斜率、截距的计算公式如式(1)~(4) 所示:

(1)

(2)

(3)

(4)

式中: r是皮尔逊相关系数，cov₁₂是在一定区域(如3×3窗口)内2个时相图像所有波段亮度值的协方差，s₁，s₂是前、后时相图像的特定区域内所有波段亮度值的标准差，BV_i1是前一时相图像的一定区域内所有波段的第i个像元亮度值，BV_i2是后一时相图像的一定区域内所有波段的第i个像元亮度值，n是特定区域内所有波段的像元总数，μ₁和μ₂分别是2个时相图像特定区域内所有波段像元亮度值的平均值，a和b分别是斜率和截距(Im et al., 2005)。

2.2.2 提取训练点/检验点数据

在试验区域随机生成700个点，根据研究区域在森林干扰条件下的实际破坏情况和森林分布情况定义了7个变化类和7个不变化类。变化类别包括:阔叶林-裸地(BF-BL)、针叶林-裸地(CF-BL)、混交林-裸地(MF-BL)、灌木林-裸地(SF-BL)、草地-裸地(GL-BL)、高山流石滩植被-裸地(FV-BL)、植被-水体(堰塞湖，V-W); 不变化类别包括:阔叶林(BF)、针叶林(CF)、混交林(MF)、灌木林(SF)、草地(GL)、高山流石滩植被(FV)、裸地(BL)。参考保护区植被类型图，并结合2010年8月野外调查样点GPS数据，确定700个点的所属类别，其中400个作为训练点，300个作为检验点，所有点的类别信息见表 1。将2个时相图像的12个波段，分别与不同窗口大小的邻近相关图像(NCI-3，NCI-5，NCI-7，NCI-9) 和该区域DEM进行图层叠加处理，形成包括多个数据图层的输入数据。不同变化检测方法的输入数据见表 2。

2.2.3 训练分类规则，生成变化检测图

研究使用Quinlan的机器学习程序(C5.0) 来开发决策树分类规则，采用属性信息增益率(gain ratio)来选择属性，最大信息增益率是属性选择及样本分类的准则。作为非参数算法，C5.0算法对训练样本点的分布没有特殊要求，可以同时处理连续和离散数据，生成的分类规则易于理解，生成规则和利用规则分类速度较快。机器学习的决策树分类，其输入数据不但可以包括原始波段数据和波段运算生成的各种指数数据，还可以包括各种GIS数据，如DEM、坡度、坡向信息等。本研究的输入数据由2个时相影像数据、DEM数据和不同窗口大小的邻近相关图像组成。为了检测结果，不同变化检测方法使用相同的训练点和检验点，提取所有图层数据，导入C5.0开发分类规则。由表 2可知:输入数据差异主要表现为是否有NCI以及不同窗口大小的NCI数据，这样可以有效验证NCI对变化检测精度的影响。

2.2.4 精度评价

在分类结果精度评价中，主要是将分类结果与地表真实信息进行比较，最常采用的方法是基于误差矩阵的统计方法。通过计算得到基本的精度估计量，如分类总体精度、制图精度、用户精度，但这些指标的客观性依赖于采样样本以及方法。Kappa分析采用离散的多元技术来克服以上的缺点，是一种测定2幅图像之间吻合度或精度的指标，其计算公式如下:

(5)

式中: r是误差矩阵中总行数; x_ii是误差矩阵中第i行、第i列上样本数; x_i+和x_+i分别是第i行和第i列的总样本数; N是用于精度评价的总样本数。

3 结果与讨论

本研究采用方形窗口作为滤波窗口，5×5窗口生成的龙溪—虹口保护区地震前后邻近相关图像局部图如图 3所示。相关系数图中灰度值低的像元表示变化区域(图 3a); 而斜率和截距图中中等灰度值像元表示非变化区域，较高和较低灰度值的像元则表示变化区域(图 3b，c)。

图 4分别列出了不使用NCI和分别使用NCI-3，NCI-5，NCI-7，NCI-9，并结合决策树分类生成的森林景观变化检测图，图 4中红色区域即为地震所导致的森林破坏区域。表 3为不同分类方法得到的地类变化面积，其中是否结合NCI图像，分类结果中地类变化面积有较大差异。

本研究对选取的300个检测点分别计算了5种分类结果的误差矩阵，相关的精度信息见表 4。

通过表 4可以看出:

1) 不使用NCI的变化检测分类总体精度为75.00%，Kappa系数为0.728 8，利用NCI-3到NCI-9的变化检测分类总体精度和Kappa系数分别为:80.67%和0.790 4，82.33%和0.808 5，79.00%和0.772 6，78.66%和0.769 0。通过比较，可以发现使用邻近相关分析图像的分类精度会高一些，NCI-5的变化检测精度最高。NCI提高变化检测精度主要体现在对变化类别的识别上，因为邻近相关图像可以有效反映影像中的变化信息，且生成的决策树规则60%以上都与邻近相关变量有关。图 5为NCI-5的变化类别分布详图。

2) 利用NCI的变化检测精度从高到低为NCI-5＞NCI-3＞NCI-7＞NCI-9，5×5窗口被认为是创建该研究区域NCI的最佳窗口大小。由于30 m分辨率图像中的高频信息已很少，选取越大的窗口NCI图像越平滑，越不利于突出变化信息，因此变化监测的精度有所降低。

使用该变化检测模型应该注意以下几点: ① 多时相图像的高精度几何校正是正确探测森林景观变化的前提; ② 多时相影像的不同视角问题也会带来误差，所以应该选择相同传感器、相同时期和视角的遥感数据; ③ 样本的精度决定了变化检测结果，选取高精度、高质量的样本非常重要; ④ 由于森林生态系统的复杂性，反映到影像上常会出现同物异谱、同谱异物的现象，因而如何有效地减少森林自身变化的影响，还需进一步探讨。

4 结论

本文主要探讨在地震干扰的条件下利用邻近相关图像和决策树分类技术相结合的遥感变化检测方法来识别森林景观变化状况。研究表明: 1) 该变化检测方法可有效监测森林景观变化，自动化程度高并且节省大量成本; 2) 邻近相关图像包含变化信息，有助于提高变化检测的精度，其中5×5邻近窗口大小监测森林景观变化精度最高，为82.3%。