基于高分二号的东川蒋家沟地表覆被变化检测分析

2021-10-19李新澳沈聪颖

地质灾害与环境保护 2021年3期

李新澳，甘淑,2，沈聪颖

(1.昆明理工大学国土资源工程学院，昆明 650093；2.云南省高校高原山区空间信息测绘技术应用工程研究中心，昆明 650093)

1 引言

随着各类传感器的不断发展，海量的遥感数据详尽地记录着各种地表地物变化过程，极大地推动了多时相高分辨率遥感数据在地表覆被变化、城市发展和各类环境问题检测中的应用。变化检测就是通过分析和解译同一区域多时相遥感影像，定性或定量地分析和确定地表变化特征和过程的技术[1]。变化检测主要有两个思路，第一，先对影像进行直接比较，后提取信息。基本思想是先对多时相遥感影像进行比较，再对比较得出的差异图像进行进一步的分类处理，最后得到变化检测结果。第二，先对影像提取信息，后比较。基本思想是先对多时相遥感影像进行分类标记，然后根据相应分类的差异来提取发生变化的区域，最后得到变化检测结果[2]。

多元变化检测(Multivariate Alteration Detection, MAD)算法的概念与方法由丹麦学者Nielsen等提出，这是一种基于典型相关分析(Canonical Correction Analysis, CCA)的变化检测方法，且在多时相遥感影像中具有明显的优势[3]。随后Nielsen在原始MAD方法的基础上引入最大期望值(Expectation Maximization, EM)算法，利用先验概率计算变化阈值，改进得到迭代加权多元变化检测算法，对多时相遥感影像进行辐射归一化，得到了较好的不变背景反衬变化区域，提高了变化检测结果的精度[4]。

东川独特的地形和地质构造，造成当地地质酥松，水土流失严重，加之长期以来铜矿开采和人为砍伐，森林面积减少，控制自然灾害的能力减弱，成为泥石流爆发的重灾区，严重威胁到当地人民的生命财产安全。自2012年以来，东川区将“生态修复”作为东川战略发展思路之一，因此检测其地表覆被的变化，为优化东川区生态修复方案提供有效的数据显得十分有必要。考虑到研究区地表覆被类型分布细碎，所采用的两期影像分别来自GF-2的PMS1和PMS2不同传感器等特点，本文采取了先比较后提取信息的方法进行检测。检测主要包括3个部分，确定地区是否发生变化；分析变化区域所在的位置；确定变化类型和变化特性[5]。

2 研究区和研究数据

2.1 研究区概况

本文选择了云南东川蒋家沟地区为研究区域，其位置分布如图1所示。研究区位于东经103°6′～103°12′、北纬26°13′～26°17′之间。东川区地处金沙江和小江深大断裂带，地质结构复杂、山体高、山谷深、沟壑纵横。气候主要为亚热带高原季风气候，降雨以夜雨和点暴雨为主且集中，复杂地形，加上人类活动，导致植被覆盖具有很大的差异[6]。

2.2 数据与预处理

GF-2是我国自主研制的首颗空间分辨率优于1m的民用光学遥感卫星，搭载有两台高分辨率1 m全色、4 m多光谱相机，具有亚米级空间分辨率、高定位精度和快速姿态机动能力等特点，有效地提升了卫星综合观测效能，达到了国际先进水平。本文采用的数据为覆盖东川区蒋家沟的两景无云GF-2影像，实验数据如图2所示。

图2 研究区遥感影像

借助ENVI5.3.1平台对遥感影像进行预处理，具体流程如图3所示。辐射定标是为了将影像的亮度转换为绝对的辐射亮度；大气校正是为了消除遥感影像中由于大气散射和吸收引起的辐射误差，本文运用了ENVI5.3.1软件中的FLAASH大气校正模块对影像进行校正；正射校正是对影像进行几何畸变纠正的一个过程，利用分辨率为5 m的DEM数据结合GF-2影像自带的RPC文件对影像进行正射校正[7]；图像融合既可以保留多光谱的真彩色信息，又保留全色波段的高空间分辨率信息，提高影像质量；几何校正是以前期影像为基准影像，配准后期影像；最后裁剪出两期影像在研究区域重叠部分。

图3 遥感影像预处理流程

3 方法

3.1 变化检测方法

为了检验IR-MAD算法的有效性，选取了PCA算法和差值算法进行对比实验。IR-MAD算法和PCA算法得到的结果通过密度分割，得到只包括变化部分和未变化部分信息的二值图像，与差值算法得到的二值图像进行对比，其中黑色区域为变化部分，白色区域为未发生变化部分，如图4所示。相比IR-MAD算法，被PCA算法确认为变化区域的面积较大，识别出的变化区域呈大面积块状分布，将水体全部归并为变化区域，对细节信息的把控能力较差，图像左下角的道路部分，PCA算法将其全部识别为变化区域；相比IR-MAD算法，差值算法将大部分未发生变化的水域部分确认为变化区域，对植被变化的识别能力较低。从以上结果的分析可以发现，PCA算法和差值算法在对蒋家沟地区进行变化检测时，存在大量错分、误分的现象，未能有效的进行检测；IR-MAD算法，根据其原理，将每个像元的权值与最终的阈值进行比较，确定像元是否发生变化，能有效的区分变化与未变化区域，达到较好的检测效果。

图4 各类算法变化检测结果图

3.2 方法运用

采用IR-MAD算法提取两期影像变化的区域。IR-MAD算法的原理是，每个像元的初始权重为1，每次迭代都赋予两幅影像中每个像元新的权重，通过计算，未发生变化的像元具有较大的权重。最终得到的权重是决定各个像元是否发生变化的唯一依据，在若干次迭代后，每个像元的权重会趋于稳定，直到变化小于设定的阈值则停止迭代[8]。设定不同的最大迭代次数和变化收敛阈值对研究区的处理结果如图5所示。在3幅灰度影像中，值越大，说明该处地物发生变化的可能性越大。对于变化较为明显的地区，几乎不受迭代次数和收敛阈值的影响，而在图5实线框内，可以明显的看出变化较为细碎的地区，受最大迭代次数和变化收敛阈值的影响较大。通过实验，将最大迭代次数确认为30，变化收敛阈值确认为0.001。

图5 不同最大迭代次数和变化收敛阈值处理结果图

在使用IR-MAD算法进行变化检测后，对该算法进行精度评价。主要是根据检测结果与实际数据建立混淆矩阵，通过Kappa系数精度指标，对变化检测结果进行评价。以 Google Earth影像作为参考影像，在影像中选择360个样本点进行精度的评定，样本点均匀分布于研究区域内，使得到的影像变化精度检测更具说服力，根据统计数据， Kappa系数为0.948。

通过IR-MAD算法，可以将研究区域划分为变化部分和未变化部分，根据研究区的具体情况，发生变化区域的地表覆盖类型划分为以下6类：建筑用地(有楼房、简易房、道路)、荒地(有岩石、裸地)、林地(有林地、园地)、耕地(有水田、旱地)、草地、水系(有河流、水沟、坑塘)。

由于研究区地处泥石流灾害区，泥石流是由粒径不等的砂粒、雨水和大量粘性土等不同的混合体组合而成的，在影像上主要表现为影像的混合像元较多。由于混合像元的影响，泥石流地区在影像上会出现较多“异物同谱”和“同物异谱”的现象[9]。全色影像与多光谱影像融合以后的GF-2影像，相较于中低分辨率影像，影像中的地物外部轮廓更加清晰，光谱异质性更强，且光谱波段数较少。基于研究区域和数据源的特点，采用面向对象分类的方法，对变化区域进行信息提取。

多尺度分割是面向对象分类变化检测的关键环节，多尺度分割是从任一个像元开始，采用自下而上的异质性最小的区域合并方法，生成影像对象[10]。由于不同的地物有不同的最优分割尺度，考虑研究区域的特点，选择恰当的参数，对于实现影像准确分类有较大的改善[11]。通过多次对比实验，GF-2的4个波段都参与分割且权重为1，分割尺度为150时，分割效果较好。同时影像进行多尺度分割时通过调节形状指数与紧凑度两个参数的取值可有效改善分割效果，两参数取值均在 0～1以内。通过参数评价可以得出规律，形状指数的值偏小，紧凑度的值偏大时，图像分割效果趋于最佳状态，如图6所示。相较于形状指数和紧凑度均为0.5的分割效果，形状指数为0.9，紧凑度为0.1，对细节的分割效果不佳；形状指数为0.1，紧凑度为0.9对细节的分割效果较好。通过对比实验，最终将形状指数确认为0.1，紧凑度为0.9。

图6 分割效果图

根据以上分割结果，应用eCognition软件，基于地物的本质特征，建立分类体系。影像对象的本质特征是指各类地物在地理空间的实际情况和成像时的传感器所决定的光谱、形状、纹理特征等物理属性特征。选择具有代表性、明显的特征，并将特征组合、建立分类规则是信息分类提取的关键。在对影像特征进行选取时，选取的特征过多会导致不同特征之间相互影响，数据冗余量过大，信息分类效率降低，应利用目标地物最为明显的本质特征进行描述[12]。林地、草地、部分农田属于植被类型，在近红外波段反射率较强，可以充分利用其光谱信息；房屋与道路具有明显的几何规则，以线形、矩形居多，应当充分利用其形状信息；水体在近红外波段的反射率最低，色调均一，连续性较强，水陆界限明显，可以充分利用光谱及形状特征。通过对研究区地物类型特征的分析，本文选取的分类特征包括：波段均值、亮度、长宽比、形状指数、归一化植被指数(normalized differential vegetation index,INDV)、归一化水体指数(normalized difference water index，INDW)，计算如下：

(1)

式中,bNIR和bR分别为GF-2近红外波段、红光波段的反射率，突出影像中的植被信息。

(2)

式中,bG为 GF-2绿光波段的反射率，抑制植被信息，突出影像中的水体。

基于上述选取的特征，选取样本进行训练并应用。发生变化区域的地表覆盖类型划分为建筑用地、荒地、林地、耕地、草地、水系6类，共选取训练样本500个，其中建筑用地样本79个，荒地样本91个，林地样本85个，耕地样本82个，草地样本93个，水系样本70个。在选取的训练样本上利用最近邻分类方法基于所选取的特征进行分类，选出最佳特征组合，然后应用选出的最佳特征组合，对影像进行分类。

4 结果与分析

本文基于GF-2多光谱数据与全色数据融合以后的两期影像，利用IR-MAD算法，提取两期影像的变化区域，通过对比实验，选取形状指数为0.1，紧凑度为0.9，分割尺度为150，进行多尺度分割，结合变化区域地物的波段均值、亮度、长宽比、形状指数、INDV、INDW等特征，对影像进行分类，通过QGIS对分类结果统计，得到2015年和2017年变化区域的地表覆被类型示意图，如图7所示。该研究区域内，大部分地区的地表覆被未发生变化，发生变化的区域在研究区内呈现碎块状分布，主要集中于北部和东部地区；在研究区的中北部靠近居民点分布地区的大量草地转变为了耕地；东北部有大部分的荒地转变为了草地和林地，该区域也是林地与草地相互转换分布较多的地区，草地与林地交错分布，相比林地，草地的分布更为细碎；在研究区的中部地区，沿着沟谷堤岸分布着大量荒地，此区域内的部分荒地转变成了草地和林地；靠近居民点的荒地，主要转变为了耕地；在建筑区附近多分布着荒地与草地相互转换的区域。

图7 2015年、2017年变化区域地表覆被类型图

对数据进行统计，如表1所示。2015年，用地类型发生变化的荒地占变化总区域的36.82%，草地所占比次之，为27.16%，林地和耕地所占的比例分别为16.25%、10.18%，建筑用地和水系的占比较小。到2017年，草地、耕地、林地、荒地占变化总区域的比例分别为32.16%、25.69%、23.02%、10.02%。2015～2017年间，在变化区域内，荒地的占比是在下降的，草地、林地、耕地的占比是依次提高的。