土地利用和土地覆盖变化(LUCC)是全球环境变化研究和可持续发展研究的核心问题^{[1, 2, 3]}。它不仅客观记录了地球表面景观的空间格局，也集中体现了自然因素和人类活动改变地球表面景观的时空变化过程^[4]。遥感技术因其具有快速、实时、宏观、周期性等特点而成为土地覆盖变化检测主要的技术手段^[5]。变化检测就是从不同时期的遥感影像中，定量的分析和确定地表变化的特征与过程，是土地覆盖遥感监测的关键环节^[6]。

遥感影像变化检测按最小处理单元可以分为基于像元和面向对象两种。随着影像空间分辨率的提高，基于像元的变化检测存在信息单一、易产生离散像元等不足而受到一定限制^[7]。相比而言，面向对象的检测方法具有提高影像信息利用程度，避免“椒盐效应”等特点而受到广泛的重视与应用^{[8, 9, 10]}。其核心思想是在变化检测之前将影像分割为若干个具有相近特征的“同质”对象，提供丰富的形状、纹理、拓扑等语义信息，并能结合人工智能进行分析，从而达到提高检测精度的目的。在相关方法研究中，Lambin等^[11]利用变化矢量分析对Sudanian地区的土地覆盖变化进行了检测，该类方法易于实现但对变化阈值的确定具有主观性。盛辉等^[12]在相关分析变化检测中采用基于贝叶斯最小误差理论的方法确定变化阈值，避免人为确定阈值的不足。Im等^[8]和张雨霁等^[13]分别将人工智能和机器学习引入变化检测中，建立了基于多尺度分割和决策树的变化检测方法，有效提高了检测精度。为了更加充分利用面向对象所提供的丰富信息，Gamanya等^[14]构建TM/ETM+影像对象的多种光谱、形状和纹理特征，利用决策树提取分类规则对津巴布韦首都Harare的城市扩张进行了检测，精度达到90%以上。以上研究为发展复杂地形条件下适用的变化检测方法提供了参考与借鉴。

相比平原地区，山区遥感影像变化检测由于地形效应的影响而面临以下问题：①几何畸变严重，增加影像配准难度，发生偏移的像元易被检测为虚假变化；②由于山体阴影的遮挡，阴坡地物的光谱信号十分微弱，不同成像条件下易被误检为变化信息；③山区地形起伏改变传感器-地物间的几何成像关系，导致辐射信号失真；④山区地形起伏导致土地覆盖格局破碎，细微的变化信息在检测过程中容易被忽略。同时，山区植被类型丰富且同类间光谱异质性更强；土地覆盖变化方式多样，植被物候特征变化迅速，增加变化检测难度。因此，构建复杂条件下适用的变化检测方法是山区大面积土地覆盖遥感监测的重要技术。然而，目前针对山区的遥感影像变化检测研究较少，且从已有研究来看还存在以下不足：未充分考虑地形效应对变化检测的影响,以及地形特征和形状特征等在山区变化检测中的独特作用；检测方法单一，易将山区地形阴影等“噪声”误检为变化信息；未充分考虑山区空间格局破碎、覆盖类型丰富、变化方式多样等对变化检测的影响。

环境保护部和中国科学院于2011年底联合启动了全国生态环境十年(2000—2010年)变化遥感调查与评估项目，土地覆盖制图和动态变化监测是该项目的主要研究内容和重要基础。中国科学院水利部成都山地灾害与环境研究所承担了西南片区(以山区为主，包括四川省、重庆市、西藏自治区、云南省、贵州省)土地覆盖制图及十年动态变化监测任务。论文在充分总结项目成果和已有变化检测研究的基础上，构建了一种山区复杂条件下面向对象的遥感影像变化检测方法。并以四川攀西地区的多期TM影像为例，对方法中涉及的影像预处理、多特征构建、训练样本选择、变化信息规则提取和精度评价等内容进行了详细的阐述与探讨，以期为复杂山区遥感影像变化检测提供一种适用方法和借鉴。

选择Landsat-TM轨道号为P130R041(WRS2)对应的四川省南部典型山地为研究区(图 1)，主要包含凉山州和攀枝花市的部分区域，简称攀西地区。该区位于横断山脉构造带，海拔介于488—4409m之间，地形起伏明显。区内水系发达，包括金沙江、雅砻江、安宁河等长江干支流，地貌以深切河谷为主。结合遥感影像和野外调查可知，研究区土地覆盖类型丰富，空间格局破碎且异质性强。作为四川省重要的金属冶炼基地、水电基地和经济作物种植区，攀西地区的土地覆盖格局在过去20年间发生了剧烈的变化。此外，安宁河谷平原的耕地在两景影像中分别处于休耕期和生长期，光谱差异明显。综上，攀西地区在地形条件、覆盖类型以及变化方式等方面都具有代表性和典型性，适于作为复杂山区遥感影像变化检测方法研究的试验区。

图 1 研究区概况 Fig. 1 Overview of Panxi region

图选项

变化检测所用TM影像来源于马里兰大学GLCF(The Global Land Cover Facility)，分辨率为30m。成像时间分别为1989-05-11(基准影像，记为Image1)和2009-04-16(检测影像，记为Image2)，其中含云量均小于1%。地形数据为最新版本的SRTM 3产品数据，分辨率为90m，其全球绝对高程精度为±16m，绝对平面精度为±20m^[15]，能够体现出地形的细微起伏变化。

相比平原地区，山区地形起伏使影像几何畸变更加严重，通过同名点分析发现两景影像间存在约2个像元的位置偏移，不能满足变化检测叠加分析的精度要求。为此，采用美国NASA下设LEDAPS课题组开发的自动配准与正射纠正程序包AROP对两景影像进行正射纠正和几何配准^[16]，已有研究结果表明AROP对于山区TM影像的配准误差小于0.5个像元^[17]。同时为了消除山区气溶胶对辐射值的影响，增强影像对比度和地物间的区分能力，同样采用LEDAPS课题组开发的大气校正模型对影像进行校正^[18]，完成从DN值到反射率(扩大10000倍，0—9999)的转变。在此基础上，通过叠加分析确定重叠区域并进行掩膜和裁剪处理，得到重叠区影像大小为7122列×6569行，有效总面积约2.97万km²。通过影像判读，攀西地区土地覆盖典型的变化方式如表 1所示。

表 1 攀西地区1989—2009年土地覆盖典型变化方式 Table 1 The typical way of land cover change in Panxi region during 1989 to 2009

影像 Image	典型变化方式 The typical changes
	矿山开采 Mining	退耕还林 Grain for green project	城镇扩张 Urban sprawl	兴建水库 Reservoir construction	耕地扩张 Farmland expansion	森林损毁 Forest damage

表选项

基于对象和规则自动提取的山区遥感影像变化检测的主要研究内容包括：设定合适的参数对影像进行多尺度分割以获得同质性高、大小适中的斑块对象，构建对象的光谱、地形和形状特征，结合地面调查数据选择“变化”和“未变化”对象的典型训练样本，利用决策树算法提取检测规则并对变化信息进行自动提取，最后对检测结果进行精度验证以评价方法的合理性,并提出改进措施。。

影像多尺度分割是面向对象遥感分析的重要先决条件，其基本思想是综合考虑多光谱影像的颜色(光谱)特征和形状特征等因素，采用自下而上的迭代合并算法将影像分割为高度同质性的斑块对象^[19]。对象的同质性可以通过斑块内像元的标准差来衡量，异质性则由对象的光谱异质性和形状异质性共同确定^[20]。对于遥感影像分析而言，光谱特征是决定对象的首要条件，形状特征的引入则会提高分割结果的质量，光谱和形状两个因子的权重之和为1。而形状异质性又通过光滑度和紧凑度来衡量^[21]，这两个因子的权重之和也是1。影像多尺度分割的主要参数包括：①波段及其权重，信息载量大的波段可以赋予较高的权重。②分割尺度，决定生成斑块的大小及对象层的破碎程度。③均质性因子，包括光谱与形状因子(光滑度和紧凑度)的权重，通常情况下光谱特征最为重要，其权重值大于0.6。

文中采用德国Definiens Imaging公司开发的面向对象遥感处理系统eCognition Developer 8.7对影像进行多尺度分割。该平台能够充分利用对象特征和类间关系等信息，是当前应用较为广泛的面向对象分析软件。

无论是分类还是变化检测，影像特征都是起决定作用的关键因素。面向对象的遥感分析不仅可以利用影像的光谱特征，还可以利用影像的形状特征和地形特征等，为山区遥感影像变化检测提供更加丰富的描述信息。文中构建的光谱特征包括反射率均值、标准差、变化矢量强度(CVI)^[22]、向量相似度(VS)^[23]、归一化植被指数(NDVI)、改进型归一化水体指数(MNDWI)等；形状特征包括形状指数(SI)、长宽比、面积(像元数量)等；地形特征包括海拔、坡度、坡向等。共构建特征24个，部分特征的计算公式见表 2。

表 2 部分特征的计算公式 Table 2 Formulas of partial features

特征名称 Feature name	计算公式 Formula	说明 Explanation
变化矢量强度 Change Vector Intensity,CVI		xi，yi(i=1，2，…，n)分别表示对象在基准影像和检测影像中第i个波段的地表反射率
向量相似度 Vector Similarity,VS
归一化植被指数 Normalized Difference Vegetation Index,NDVI	NDVI=(b4－b3)/(b4+b3)	bi(i=2，3，4，5)表示对象在TM影像第i个波段的地表反射率
改进型归一化水体指数 Modified Normalized Difference Water Index,MNDWI	MNDWI=(b2－b5)/(b2+b5)
NDVI差值 Difference of NDVI,Diff_NDVI	Diff_NDVI = NDVI_1989-NDVI_2009
形状指数 Shape Index,SI		A为斑块面积，P为斑块周长

表选项

文中构建的部分特征的释义如下，光谱特征中变化矢量强度CVI表示对象在不同时期光谱矢量间的欧式距离，值越大表示该对象发生变化的可能性越大。向量相似度VS表示对象在不同时期影像中光谱向量的相似性，其值介于0—1之间。VS越接近于1，说明对象在不同时期的光谱向量相似性较高，其发生变化的概率较低，反之则说明对象发生变化的概率较高。NDVI、MNDWI分别对植被和水体有较好的表征作用。而形状指数SI和长宽比对山区河流、道路等线状对象有较好的表征作用。地形特征中，坡度体现斑块所处位置的地形起伏度，坡向对阴坡和阳坡有较好的区分作用。

变化检测可以认为是只有“变化”和“未变化”两个类别的遥感影像分类过程。在众多的分类规则提取算法中，决策树作为一种非参数算法具有对训练样本分布没有特殊要求、可以同时处理连续和离散数据、规则易于理解、分类速度快、精度高等优点而被广泛采用^[24]。其基本思想是：按照一定的规则把数据集逐级往下细分以得到具有不同属性的各个子类别。分类决策树由一个根节点、一系列内部节点和终极节点组成，每一个内部节点只有一个父节点和两个或多个子节点。在每一个内部节点(包括根节点)处根据一系列规则将该节点处的数据集划分为两个或多个子集，如此往复直至所有的数据被分为预期设定的各个子集为止^[25]。目前在分类规则构建中使用广泛的决策树算法为C5.0决策树^[26]。

决策树作为一种机器学习算法受训练样本的影响较大，尤其是在复杂山区，训练样本的代表性和全面性是提高变化检测精度的关键。经过多次试验总结得出训练样本的选取原则如下：①训练样本在研究区内大致均匀分布，且在不同土地覆盖类型、变化方式、地形条件下都要有“变化”和“未变化”样本的分布；②选择一定数量的虚假变化作为训练样本，如将地形遮挡、云及云阴影、植被物候等引起的虚假变化均当作“未变化”样本。

通常情况下，精度评价的参照样本由地面实测获取。然而，历史影像的地面调查数据获取是精度评价的难点，多年来一直未得到满意的解决。尤其是实例中基准影像获取于1989年，对应时刻的地面调查数据稀少。作者所在课题组从1995年开始从事西南地区土地覆盖制图及野外调查工作，至今已积累大量全省范围内不同时期的地面调查数据，攀西地区也有一定数量的样点分布，主要集中在安宁河谷平原地区。

除了地面调查数据，研究还通过分层随机采样选取部分参照样本作为补充。具体步骤为：将影像划分为50×50个大小一致的矩形区域，在每个矩形区域内随机选取某个对象并结合多期历史遥感影像、Google Earth等进行综合判读，避免人为选择参照样本的主观性和随意性。在此基础上，通过构建误差矩阵来对检测结果的精度作出评价。

经过多次对比试验，确定攀西山区遥感影像多尺度分割的最优参数设置如下：波段为两景影像除热红外(TM6)以外的12个波段，各个波段的权重均设为1；分割尺度设为25(相对于大气校正后的反射率数据)，光谱权重和形状权重分别设为0.9和0.1，光滑度和紧凑度权重分别设为0.3和0.7。攀西地区部分典型地物多尺度分割的结果对比见图 2。

图 2 多尺度分割结果 Fig. 2 Results of multiscale segmentation

图选项

从图 2中典型地物的分割细节可以看出，在上述参数设定条件下，多尺度分割获得了较为满意的结果。主要体现为：分割结果不至于过于破碎但能够体现出斑块间的差异信息，对象内部均质性高且与地物实际边界吻合较好。相对而言，变化区域的分割结果更加细致，如图 2(a)、图 2(b)中的矿场和水库；未发生变化且影像光谱差异不大的区域分割斑块则相对较大，如图 2(d)中的水库。变化检测必须选用两景影像参与多尺度分割，其原因在于若只采用其中一景影像参与分割，则该景影像中的均质斑块在另外一景影像中可能是一个包含多种地物的混合对象。如图 2(b)中的水库，在只有Image 2参与分割的情形下易被分割为面积更大的斑块，不能体现出Image 1中的河流信息，并由此降低对象的均质性，给规则提取及检测结果带来影响。

结合课题组积累的历史地面调查数据和Google Earth高分辨率影像，对试验影像共选取“变化”和“未变化”训练样本393个，其中“变化”类型的训练样本数为152个，“未变化”类型的训练样本数为241个。在eCognition中完成2.2节中描述特征的构建并将其导入See5.0决策树软件中，自动提取遥感影像变化检测规则见图 3。

图 3 决策树 Fig. 3 Decision tree NDVI_1989和NDVI_2009分别表示对象在两景影像中的归一化植被指数NDVI，B3_1989表示TM波段3在1989年影像中的反射率，B2_2009表示TM波段2在2009年影像中的反射率，其余变量的释义见表 2; 终极节点括号中的数字x/y释义：x表示该节点处的样本数，y表示分类错误的样本数

图选项

从图 3可知，C5.0算法自动构建的决策树包含8个特征、6个层级和12个终极节点。对训练样本的分类正确数为382个，精度为97.20%，其中“变化”和“未变化”训练样本的分类精度分别为98.02%和96.68%。在根节点处，利用Diff_NDVI=0.1343可以将训练样本一分为二，其中Diff_NDVI>0.1343的子集主要包括“变化”样本，其比例为79.10%；Diff_NDVI≤0.1343的子集则主要包括“未变化”样本，其比例为82.24%，表明本试验中NDVI差值对于变化信息有较高的区分能力。在Diff_NDVI>0.1343的子集中，主要的“变化”样本检测规则为：(NDVI_2009≤0.2194)和(B5_1989>631.05)，该部分“变化”主要由植被向非植被转变，其NDVI在2009年的影像中较小。该子集中的“未变化”样本主要包括地形阴影、物候等因素导致的虚假变化，其NDVI在两景影像中的差值同样较大，但其属性并未发生变化。而在Diff_NDVI≤0.1343的子集中，主要的“未变化”检测规则为：(CVI≤5289.4)和(Diff_NDVI>-0.1775)，该类样本在两景影像中的NDVI差异不大(-0.1775＜Diff_NDVI≤0.1343)，光谱向量差值CVI也被限制在一定范围内，故其能将大部分“未变化”样本提取出来。该子集中的“变化”样本主要为退耕还林，其NDVI_2009明显大于NDVI_1989，Diff_NDVI一般小于-0.2。如规则(CVI≤5289.4)和(Diff_NDVI≤-0.1775)和(B2_2009≤1228.3)可以将大部分退耕还林检测出来，B2_2009≤1228.3表示植被在TM2中的反射率较小。此外，在以CVI和VS为划分依据的节点处，CVI较大(或VS较小)的样本主要为“变化”，反之则为“未变化”，这与CVI(VS)的定义与作用相符。

C5.0决策树算法对于只有“变化”和“未变化”两类且同类训练样本间差异明显的分类仍能取得满意的结果，表明该算法具有较强的适应能力和鲁棒性。总体而言，光谱特征仍然是山区TM影像遥感分析最重要的信息源，其中利用频率最高的特征为Diff_NDVI，对其设定合适的阈值可以将约80%的“变化”和“未变化”样本区分开。这是由于攀西地区植被覆盖度高，且主要的变化方式如退耕还林、矿山开采等均发生在植被覆盖区域，因此NDVI差值能够作为上述变化的有效表征。此外，NDVI_2009、CVI等特征在决策树构建中也得到了充分的利用。

在eCognition 8.7平台中构建图 3中的分类决策树对试验影像进行变化检测。同时为了消除结果中的碎小斑块，使结果更加美观和便于后续应用，将变化区域中像元个数小于6的斑块进行合并或者去除(6个像元为项目要求最小制图单元)，得到攀西地区1989—2009年遥感影像和典型变化的检测结果如图 4所示。可以看出，文中构建的方法对各类典型变化均取得了理想的检测结果，提取的变化信息与实际变化区域在空间上能够较好的吻合。此外，图 4b-B中分别处于休耕期和生长期的耕地并未被检测为变化信息，表明该方法对物候导致的虚假变化具有较强的抗干扰能力。

图 4 变化检测结果 Fig. 4 The results of change detection

图选项

从检测结果来看，攀西地区20a间土地覆盖格局变化明显，共检测出变化面积约740.2km²，占影像总面积的2.49%。其中典型的变化方式包括退耕还林、水利工程、矿山开采和城镇扩张等(图 4b)，人类活动是攀西地区土地覆盖格局变化最主要的驱动因子，如退耕还林(草)(1999年实施)、二滩水电站(2000年建成)等重大工程迅速改变了该区的土地覆盖格局。典型变化的空间分布特征为：退耕还林主要分布在具有一定地形起伏的坡耕地上；耕地扩张主要分布在盐源平原北部的经济作物种植区，以及城镇和水利工程附近人地矛盾相对突出的地区。雅砻江河谷两侧由于二滩水电站的修建，大量林地、耕地被淹没。除了人为因素的影响，自然条件下主要的变化方式为低覆盖度的草地、灌木向高覆盖度的灌木、林地转变。

利用课题组积累的历史地面调查数据和分层随机采样共获取参照样本249个，其中“变化”样本113个，“未变化”样本136个，计算得到的误差矩阵和相关统计量见表 3。

检测结果精度评价的相关统计量如下：总体精度为93.57%，Kappa系数为0.8706，表明方法的总体检测精度较高。其中“变化”对象的漏检率和误检率分别为6.19%和7.83%，容易被漏检的变化信息主要为退耕时间过于短暂的耕地，其还未发育为典型的森林，光谱特征与农作物相似，检测过程中容易被忽略；误检率(7.83%)稍高表明该方法对变化信息较为敏感，其中易被误检为“变化”的覆盖类型包括部分落叶林、季节性变化的草地以及部分不同生长季的耕地。“未变化”对象的漏检率和误检率分别为6.62%和5.22%，误检率(5.22%)较低表明该方法不易将山区众多“噪声”因素检测为变化信息。

论文在对影像进行多尺度分割的基础上，构建对象的一系列光谱、形状和地形等特征，将地形阴影、物候等“噪声”引起的虚假变化当作“未变化”训练样本，利用决策树算法自动提取检测规则，建立复杂条件下适用的遥感影像变化检测方法。该方法对攀西复杂山区的遥感影像取得了较好的检测结果，总体验证精度达93.57%，为大面积山区变化检测提供了参考和借鉴。主要结论包括：光谱特征仍是TM影像遥感分析主要的信息源，选入决策树的特征以地表反射率及NDVI、CVI等衍生光谱特征为主，其中NDVI差值对于植被覆盖区域土地覆盖格局的变化具有较好的表征作用。C5.0决策树算法具有自动选择少数最佳特征及其阈值的能力，对于只有“变化”和“未变化”两类且同类训练样本间差异明显的影像分类仍能取得满意的结果。攀西地区1989—2009年间土地覆盖格局变化明显，典型变化方式包括退耕还林、水利工程、矿山开采等。共检测出变化面积740.2km²，占影像总面积的比例为2.49%，人类活动是土地覆盖变化的主导因素。

虽然文中构建的变化检测方法取得了较为满意的结果，但在复杂山区开展遥感影像变化检测仍然面临许多挑战。

正射校正和影像配准是山区遥感影像变化检测的关键预处理过程，且随着空间分辨率的提高，影像的配准难度增大。当像元位置存在偏差时，发生偏移的像元可能被检测为变化信息。传统基于像元的分析方法对几何配准的精度要求较高，有研究认为当配准精度达到0.5个像元时可以忽略配准误差对变化检测的影响^[27]。相比之下，面向对象的分析方法用近似“同质”的对象代替单个像元作为最小处理单元，对象内部的光谱差异得到平滑，在一定程度上减弱配准误差的影响^[28]。文中采用AROP程序包对山区遥感影像进行正射校正与几何配准，精度优于0.5个像元，由配准误差导致的虚假变化几乎可以忽略。除了几何畸变，山区起伏地形还会导致阴坡地物的光谱信号失真，此时对影像进行地形辐射校正^[29]有利于减少地形造成的虚假变化。

山区空间格局破碎，设定合适的输入波段、分割尺度和均质性因子，获取合理的分割结果是变化检测的重要基础。在没有先验矢量数据的情况下应选择两景影像参与分割，以保证信息的充分利用和分割结果的均质性。此外，还可以将遥感影像与土地利用图等先验数据进行套合，以先验数据的矢量边界为准进行分割，使对象边界与现有数据保持一致，便于叠加分析^[30]。分割尺度决定对象的同质性、差异性以及破碎程度。若分割结果过于破碎，则会增加同类地物间的光谱异质性，改变地物的形状特征并增加数据处理的计算量。若分割结果过于粗糙，则会使不同类的光谱差异变小，同时降低对象内部的均质性，掩盖细微变化信息。均质性因子表示光谱和形状在分割中的权重，文中的研究结果表明光谱特征仍是TM影像主要的信息源，其权重应赋予较大的数值。通过试验还发现，在其它参数设定一致的情况下，增加形状因子的权重会使分割结果趋于破碎。

与多数影像分类和变化检测不同，文中的训练样本只有“变化”和“未变化”两类，且同类训练样本间的差异明显，增加规则提取的难度。为此，在选择训练样本时需注意样本的代表性、全面性及其在空间分布的均匀性。此外，通过将地形阴影、物候等因素造成的虚假变化作为“未变化”样本参与规则提取，有利于提高虚假变化的识别精度。针对历史地面调查数据难以获取的难点，可以利用多源遥感数据降低其影响，如在样本的选择过程中，对研究区多期、不同时相的遥感影像开展综合判读来提高样本属性辨识的精度，降低成像条件、植被物候等“噪声”的影响。

从检测规则来看，形状特征并未充分发挥作用，究其原因可能在于：影像空间分辨率不足，变化检测所用的TM影像分辨率为30m，形状特征对地物边界的刻画存在不足。目前在形状特征能给对象提供准确描述信息的研究中，多是基于QuickBird、IKONOS和WordView等分辨率达到米级的遥感影像开展^{[31, 32]}。除了分辨率的影响，山区多为自然植被且空间格局破碎，大多数覆盖类型的形状特征无明显规律可循，形状特征难以发挥作用。如何将景观尺度(如30m分辨率)的形状特征有效的引入到遥感分析中值得进一步研究。

除了形状特征，地形特征在检测规则中也未得到体现，原因可能在于遥感影像的分辨率为30m，而地形特征是基于90m分辨率的SRTM计算得出，两者在空间尺度上的差异是否对检测结果带来影响将在后续的研究中开展。然而，地形特征在土地覆盖变化检测中仍有用途，如研究区中矿山开采、兴建水库等大型工程必然导致其地形发生变化，在检测过程中加入不同时期的数字高程模型结合某种算法即可将变化区域的位置及结构信息提取出来^[33]。此外，本研究只是确定变化信息的空间位置，在变化信息的分类过程中增加坡度信息可以提高“退耕还林”类别的识别精度，因为该区退耕还林(草)工程多发生在具有一定地形起伏的坡耕地上面。可以预见，随着影像分辨率的提高以及相关研究的深入，形状特征和地形特征将在面向对象遥感分析中发挥独特的作用。