协同光学与雷达遥感数据的面向对象土地覆盖变化检测

翟鹏飞，李世华,2※，胡月明,3,4

（1. 电子科技大学资源与环境学院，成都 611731；2. 电子科技大学长三角研究院（湖州），湖州 313001；3. 海南大学热带作物学院，海口 570228；4. 广州市华南自然资源科学技术研究院，广州 510630）

0 引 言

社会的发展和人口的急剧增加导致人类对土地资源的需求也在不断增长，土地利用/覆盖（Land use and Land cover）变化成为全球或区域变化领域的研究热点[1-5]。遥感技术的快速进步使得大面积数据获取手段更加丰富，这对土地覆盖变化检测提出了更高的要求和挑战。目前，国内外很多学者利用光学数据提出了不同的方法并将这些方法应用在像素级的变化检测任务[3-4]。基于像元的方法是对前、后时相的同一位置的像元进行比较，通常采用差值或比值等代数运算获取差异影像，然后使用不同的阈值分割方法来确定该像元是否发生变化[6-7]，这对前、后时相遥感影像的配准、辐射校正等预处理有着较高的要求，同时会使得检测结果出现“椒盐现象”。随着影像空间分辨率逐渐提高，包括米级、亚米级数据的问世，同一地物在影像中往往由多个相邻像元组成，单个像元的光谱特性不足以反映相应地物的光谱特征，而相邻地物的像元多呈现一定的关联性，同时，“同物异谱”或“同谱异物”的现象极大地影响了高分辨率变化检测的精度。随着基于地理对象的遥感影像分析方法的发展，针对对象的变化检测方法以影像分割对象为基础，以具有相似光谱和空间特征的影像对象作为基本处理单元，考虑了像元及其邻域的光谱、几何、纹理及空间上下文特征，能够减少图像配准结果的影响，极大地提高了变化检测的精度，成为高分辨率光学遥感影像变化检测的研究重点[8-12]。

结合合成孔径雷达（Synthetic Aperture Radar，SAR）数据和面向对象方法在土地覆盖变化检测研究中取得了良好的效果。尤红建[13]采用基于方差的区域增长算法实现多时相SAR 影像联合分割，并计算斑块的差异指数获得变化检测结果。张明哲等[14]使用简单线性迭代分割方法对单极化Radarsat-2 数据进行分割，并融合像素和对象结果，通过众数投票的方式获得最终结果。Zhao 等[15]提出了利用图像回归的方式获取变化差异图，在超像素分割的基础上利用C 均值聚类获得检测结果。以上研究中将影像分割算法直接应用在SAR 数据上，然而由于SAR传感器与光学传感器不同，SAR 成像是主动相干成像，成像机理复杂，且受斑点噪声影响严重，各种图像分割的方法在SAR 影像难以获得良好的边界，进而影响后续的检测结果，这使得不同时相的SAR 影像的变化检测更具有挑战性。同时利用SAR 数据进行变化检测多使用单极化或双极化数据，也在一定程度上延缓了全极化SAR数据在土地覆盖变化检测应用的发展。西南地区常年多云多雨，高质量的光学影像数据通常难以获取，而全极化SAR 传感器重访周期较长，二者单独应用于土地覆盖变化检测任务难以满足获取高精度变化结果的要求。针对上述不足，提出一种协同光学与雷达遥感数据的面向对象土地覆盖变化检测的方法，该方法在一定程度上避免了单一使用两种数据的缺陷，充分发挥了二者的优势，为土地覆盖变化检测提供了一种新的范式。2 个样区包括多种不同地物类型并且涵盖了该地区多数的变化类型。所有数据都进行了相关的预处理并基于光学数据进行配准。为了验证本文方法的可行性，选择两个试验区进行抽样检验，样区大小均为488×488，如图2和图3 所示。试验区数据如表1 所示。

表1 试验区数据详情Table 1 Details of the experimental area data

1 材料与方法

1.1 研究区概况与数据

研究区位于四川省眉山市（ 30°08′25.31″～30°18′42.39″N，103°42′0.92″～103°53′50.92″E）。眉山市地处四川盆地西南边缘，成都平原的西南部，是成都平原连通川南、川西的咽喉。眉山市辖二区四县，境内河网密集，主要有岷江和青衣江两支水系。眉山市地处中低纬度，总体地势西高东低，南高北低，属中亚热带湿润季风气候，夏季温润多雨，多种植水稻等粮食作物。随着经济的发展，城市化建设的快速推进，眉山市大量水稻种植区改种席草、果林等经济作物；部分裸地也随着时间的推移施工完成，形成建设用地。

分别使用2019年Sentinel-2 光学数据和2016年、2019年全极化Radarsat-2 雷达数据，影像大小为2 000×2 000像元。图1 中红色和蓝色框图分别表示样区一和样区二。

1.2 研究方法

研究路线如图4 所示，主要包括以下步骤。首先使用分型网络演化分割算法对2019 年Sentinel-2 影像进行分割，得到初始分割结果，然后对前、后时相Radarsat-2数据利用对数比值法获得差异影像。利用初始分割过程中获取的矢量边界作为基准，再次使用分型网络演化分割算法对雷达差异影像进行引导分割获取最终分割结果；选取合适的变化与未变化样本并提取样本特征空间，将样本间距离度量可分性特征空间优化方法（Feature Space Optimization, FSO）结果与随机森林（Random Forest，RF）分类器有机的结合构成变化检测框架（以下称为FSO-RF），利用该框架获取最终的检测结果。

1.2.1 遥感影像分割

根据分割原理的不同，图像分割方法分为自顶向下的区域分裂和自底向上的区域生长法，本文采用的分割方法为Baatz 等[16]在2000 年提出的分型网络演化分割算法（Fractal Net Evolution Approach, FNEA），该算法根据区域生长的思想，基于图像对象的光谱和几何特征，由小对象逐步合并为较大的对象，并保证在不同对象中保持较大的异质性。对象的异质性由光谱异质性和形状异质性确定，而紧致度因子和光滑度因子共同确定形状异质性。异质性由以下公式定义：

对象的光滑度用来表示对象分割边界的破碎程度，利用对象边界的周长l与近似边界的长度b的比值来表示。

式中hcolor_m、hcompact_m、hsmooth_m分别表示对象合并后的光谱异质性、紧致度异质性和光滑度异质性，nobj1、nobj2、lobj1、lobj2、bobj1和bobj2分别表示合并前2 个影像对象的像元个数、对象边界周长以及对象的近似边界长度；nmerge、lmerge、bmerge分别表示影像合并后的像元个数、合并后对象边界的周长以及合并后对象的近似边界的长度。

分割尺度也是分割算法的重要参数，尺度是一个抽象的概念，它控制了分割过程中形成对象的大小和最终分割获得的对象个数，尺度参数不同，对象的大小和个数也不同。随着尺度参数的增加，分割获得的对象越大，最终获得的对象个数相应越少。本文目的在于分割获得的对象内部应仅包含变化或未变化像元，既要保证良好的分割边界，又不要使得分割结果过于破碎。

遥感影像分割的理想结果是分割后，对象内部具有较高的同质性，而相邻的不同对象间呈现良好的异质性[17]，并且已有研究表明，基于对象的方法能更好的抑制SAR 影像相干斑噪声对图像信息的破坏[18]。本文使用对数比值法获得雷达差异影像，该方法相较差值法、比值法能够更好的抑制噪声[19-21]。利用分型网络演化分割算法对Sentinel-2进行分割，分割的结果由波段权重、尺度参数、形状参数和紧致度构成。由于所选的研究区土地覆盖类型复杂多样，地物十分破碎，经过多次试验将Sentinel-2 影像的各波段权重设置为1，分割尺度设置为50，形状参数设置为0.2，紧致度因子设置为0.3，并以初始分割结果最为基准引导雷达差异影像进行二次分割获得最终的分割结果。

1.2.2 距离度量可分性特征空间优化

遥感影像经过分割形成对象后能够提供更加丰富的特征，包括光谱、纹理、形状、空间上下文等特征，可以充分利用这些特征来提高变化检测的精度。由于特征信息的冗余等因素，当特征的数量增加到一定程度时，变化检测的精度不仅不会继续增加，甚至有可能降低[10]，因此有必要对提取的特征空间进行优化，以获取最佳的特征集合进行变化检测任务。距离度量可分性特征空间优化方法的步骤和基本思想为选取可靠的变化、未变化样本并提取相应特征构建特征空间，其次在各个特征子空间中计算变化、未变化的样本可分距离，理论上可分距离越大则表明二者差异越大，可分性越好。对每一个维度所有样本的距离进行比较，获得每个维度的最大值，代表该维度空间下的最佳可分距离，然后比较每个维度的最大值，这些值中的最大可分距离，其所代表的特征组合即为最优的特征组合。变化与未变化类的可分距离由下式定义：

式中dmax表示在最优特征空间下的最大可分距离。

使用距离度量可分性特征空间优化的方法来寻找变化检测任务的最优特征空间。本文依据HH、HV、VV 三种极化方式，Freeman-Durden 分解获得的双次散射（double-bounce scattering，dbl）、体散射（volume scattering，vol）和表面散射（surface scattering，surf）参数以及Pauli 分解获得的Pauli_r、Pauli_g 和Pauli_b 参数来计算后续所要提取的特征。在选取置信度较高的变化与未变化样本后，提取样本的特征，包括HH、HV、VV后向散射均值和标准差，Pauli 分解参数的对象均值和标准差以及Freeman-Durden 分解参数的对象均值和对象标准差，同时根据Haralick 提出的灰度共生矩阵（Gray-level Co-occurrence Matrix,GLCM）提取纹理特征[22]，本文利用同质性（Homogeneity）、对比度（Contrast）、熵（Entropy）、均值（Mean）、标准差（Standard Deviation）和相关性（Correlation）6 种纹理特征。利用灰度共生矩阵可以按照不同方向进行扫描求得对应参数的值，包括0°、45°、90°以及135°，本文按照所有方向扫描并求取均值的方式求得对应纹理特征的值，共提取72 维特征如表2 所示。

表2 特征参数变量Table 2 Characteristic parameter variables

根据式（8）～（10）计算各维度特征子空间的变化、未变化样本的可分距离，并在各个维度中获取该维度的最大值，这些值中的最大可分距离所代表的特征组合即代表最佳特征组合。

1.2.3 随机森林面向对象变化检测方法

随机森林（Random Forest）是近年发展起来的一种机器学习模型，由Breiman[23]首先提出。由于其可以适用于高维特征样本和大的数据集上，因此已经广泛应用于遥感领域，并且在分类和土地覆盖变化检测的任务上获得了良好的效果[10,11,24]。随机森林的本质属于集成学习，该模型的理论基础是决策树，它由不同的决策树进行组合，最终的结果根据每棵决策树生成的结果进行投票产生最优的分类结果。在随机森林中，误差估计是在训练过程中进行的，训练数据根据Bootstrap 重采样获得不同的训练集，利用这些不同的训练集分别构造决策树；在采样过程中使用Bagging 的原理，同一个样本可以被多次选择，而其他样本可能没有被选择，大约由三分之一的样本没有被用于每个决策树的构建，这些样本用来测试以获得袋外误差（Out of Bag ,OOB）估计。

在初始样本选择时，规定以下样本选择原则：假设在分割对象内部共有Nn个像元，其中属于变化的像元个数为Nc，属于未变化的像元为Nu，显然Nn=Nc+Nu。定义对象内部变化率μ如下：

当μ=0，表示对象内部均为未变化像元，则该对象选择为未变化样本，当μ=1，则表示对象内部均为变化样本，相应的选择为变化样本。当μ>80%或μ<20%时，在该置信区间选择置信度较高的变化与未变化样本。图5 为置信度较高的变化与未变化样本示意图。

1.2.4 检验方法

为了验证所提方法的有效性，有必要对变化检测结果进行定量的精度验证，本文在引入混淆矩阵后采用准确率（Accuracy）、召回率（Recall）、F1、错分误差（Commission Errors, CE）和漏分误差（Omission Errors,OE）进行精度评价。

准确率表示正、负类别均被正确分类的概率。

式中TP 表示将正类预测为正类数量；TN 表示将负类预测为负类数量；FP 表示将负类预测为正类数量；FN 表示将正类预测为负类的数量。

2 结果与分析

2.1 面向对象分割

以样区一为例进行对比试验，将分割算法直接应用于雷达差异影像，并将分割结果套合于光学影像用于可视化对比分析如图6 所示。图6a、6b 分别表示对光学数据分割和利用光学数据引导SAR 影像分割；图6c 是直接将分割算法应用于雷达差异图，将分割结果套合于光学影像如图6d。

将图6a 中红色框图内区域的分割细节作为展示如图7，重点关注图7 中白色框图内的地物，其内部地物类型由草地变化为建设用地。根据分割算法获得的分割结果来看，无论是直接对SAR 影像分割还是利用光学影像引导SAR 影像分割都可以正确提取变化的范围，但是后者可以在雷达差异图中获得良好的地物边界，获取更加精细的变化地理实体，并且针对变化剧烈的区域分割对象内部获取了尽可能多的变化或未变化像元。

2.2 FSO-RF 变化检测框架结果分析

在选取置信度较高的变化和未变化样本后，依据选的样本计算样本各维度特征的可分距离如表3 所示，表中显示了各维度样本可分距离大小和具体的特征参数，以一维、二维特征为例，样本特征一维空间中GLCM Con_surf 特征使得72 个一维特征的变化、未变化样本的可分距离最大，因此它的距离代表了一维特征空间的可分距离。样本特征二维空间中GLCM Con_surf 与GLCM Mea_dbl 特征在各个二维特征组合中获得了最大的可分距离，因此它代表了二维特征空间的可分距离，其余同理。最终在第45 维空间中获得最大可分距离，理论上该距离所代表的特征空间即为本次任务的最佳特征组合。

表3 各维度特征可分距离Table 3 Separable distances for each dimensional feature

图8 直观地展示了样本特征距离可分性大小随着特征数量增加的变化趋势，随着特征的增加，样本间的可分距离先增加，并在第45 维特征空间获得了最大的可分距离，随后可分距离减小，在一定程度上表明了更多的特征不一定全部有益于变化检测，也凸显了特征优选的必要性。

为了验证距离度量可分性特征空间优化方法的有效性并实现FSO-RF 变化检测框架，本文将提取的特征按照维度由小到大的方式根据样本间可分距离的计算结果正序的将各维度的特征输入到随机森林分类器中，获取各维度特征下变化检测结果的准确率如图9 所示。准确率的总体变化趋势与不同维度的可分距离的计算结果基本吻合，同时在第 41 维特征获得了最高的准确率为92.90%，此时的样本类间可分距离为1.653，在一定程度上证明距离度量可分性特征空间优化方法的有效性。

2.3 精度评价与讨论

本文采用像素的对数比值法获得雷达影像差异图，然后根据分型网络演化分割方法对光学影像进行分割，并利用该分割结果引导雷达差异图进行分割获得最终分割结果。选取置信度较高的变化与未变化样本并提取样本特征，利用本文提出的FSO-RF 框架获得最终的变化检测结果。

将本文的方法与基于雷达差异影像分割的方法、变化向量分析（Change Vector Analysis, CVA）、主成分分析（Principal Component Analysis, PCA）、多元变化检测（Multivariate Alteration Detection, MAD）、迭代多元变化检测（ Iteration Multivariate Alteration Detection,IR-MAD）、以及基于像元的支持向量机（Support Vector Machines, SVM）和随机森林（Random Forest, RF）的方法进行对比分析。其中基于雷达差异影像分割的方法将分型网络演化分割方法直接应用于雷达差异影像，尺度参数设置为7，形状因子和紧致度因子分别设置为0.1 和0.4，根据1.2.3 节样本选择原则选择置信样本并提取1.2.2节所示特征后，根据本文提出的FSO-RF 变化检测流程框架获得变化检测结果；CVA 方法是根据计算前、后时相SAR 影像的特征的欧式距离获得变化结果；PCA 方法是在对特征进行降维，选取贡献度前75%的特征，然后利用CVA 获得变化检测结果；IR-MAD 是在MAD 的基础上改进，根据卡方距离进行加权迭代，最大迭代次数Max Iteration 设置为100，变化收敛阈值Delta Criterion 设置为0.001；随机森林和支持向量机方法是在选择像素级的变化与未变化样本进行分类获得变化检测结果。

同时，为了更直观地显示变化检测任务的结果，对样区一和样区二进行定性和定量的分析。图10 和图11分别展示了2 个样区的变化检测结果。各种方法的精度评价结果如表4 所示。图10a 和11a 是根据高分辨率遥感影像进行目视解译结合实地考察绘制的参考变化结果。

表4 不同方法下变化检测精度评价Table 4 Accuracy evaluation of change detection between different methods %

从整体以及2 个样区的变化检测结果和精度可以发现：

1）无论是传统经典的基于像元的变化检测方法还是有监督的SVM 和RF 方法，在研究区数据的变化检测结果精度都相对较低。这些方法在检测过程中只使用了前、后时相影像的后向散射特征和两种极化分解特征；更重要的原因在于即使采用了滤波的处理手段，SAR 影像的相干噪声也不能完全去除，这些斑点噪声在前、后时相的差异计算中以相应特征的形式干扰算法的计算，导致针对SAR 影像基于像元的方法会造成大量的椒盐噪声，对于变化剧烈的区域也难以获得良好的边界。同时也发现，改进后的IR-MAD 相较MAD 或者其他两种无监督的变化检测方法在检测结果上会减少大量的孤点；总体而言，有监督的SVM 和RF 对无监督的经典变化检测方法在精度上有不同幅度的提升。

2）基于雷达差异影像分割的方法在各项精度评价指标下大幅度提高了变化检测的精度，同时，利用对象分割可以极大程度的改善基于像元方法的椒盐噪声现象。然而，由于雷达成像原理与光学传感器成像有着本质的区别，分型网络演化分割方法不能分割出地物的边界，进而难以获得良好的变化检测结果。根据图10h、图11h也可以发现，检测结果仅能粗略地提取变化的大致范围，而难以获取细致、真实地理实体的变化结果。

3）2 个样区结果在精度上有一定差异，样区一和样区二的准确率分别为95.08%和88.16%，原因在于初始分割仅使用了后时相光学数据。针对样区一，前后时相的土地类型相差较大，多为建筑裸地或荒草地向建设用地转变，因此分割形成了不同的土地覆盖类型对象，进而分割对象内包含了变化或未变化像元；针对样区二，精度相对较低的原因在于不同类型的农用地的转变过程，样区二后时相农用地多种植为果林或成片的席草等经济作物，在使用后时相光学数据分割时，尽管已经考虑到目的是为了获取的对象内部应仅包含变化或未变化像元，但是为保证对象内部的同质性以及尽量不将地物分割的过于破碎，将成片状的经济作物分割为同一对象，而对于前一时相，该对象内部存在其他地物类型，如水稻等粮食作物，因此分割对象内部未获得尽可能多的变化或未变化像元导致在样区二精度低于样区一。

4）本文使用距离度量可分性特征空间优化方法对提取的72 维特征进行优化并在第45 维特征中获得最大可分距离，在验证的过程中表明第41 维特征中获得了最佳的变化检测结果。但是利用高维度特征在实际的变化检测应用中复杂且繁琐，而且本文方法可以计算出具体的可分距离的大小，并且在文中验证了可分距离与准确率的变化曲线基本吻合，因此可以利用可分距离大小的变化曲线，选取可分距离大小排名靠前的对应的特征或者设定当可分距离的增长率第一次小于某阈值时，选择当前维度特征进行变化检测来指导实际的变化检测应用。

5）本文将面向对象的思想与随机森林算法相结合，利用分型网络演化分割的方法对光学影像分割并引导分割雷达差异影像以获取拥有良好、真实的地理边界分割结果；同时，针对提取的多维特征根据各维度特征样本的特征距离进行特征空间优化，验证了距离度量可分性特征空间优化的方法进而获得良好的变化检测结果，各项精度指标均高于对比方法，准确率为92.90%，召回率为96.61%，F1 为96.07，错分误差和漏分误差分别为4.47%和3.39%，进一步证明了该方法的有效性和可行性。

综上所述，在多云多雨地区应用有限的光学数据引导雷达影像分割进而采用随机森林算法进行面向对象的土地覆盖变化检测研究是可行的，这种方法不仅能有效地发挥光学数据丰富光谱信息的优势，在引导雷达影像分割时能够获得良好的真实地物边界，同时能将SAR 影像所提供丰富的极化信息有机的相结合来获取良好的变化检测结果。

3 结 论

本文提出了一种基于随机森林算法的FSO-RF 土地覆盖变化检测算法框架，该框架以光学影像引导SAR 数据分割为基础获取良好的真实地物边界，在选择置信样本提取多维特征空间后，验证了变化检测结果精度与不同维度的特征之间的关系，并将距离度量可分性特征空间优化方法和随机森林算法相结合构建FSO-RF 变化检测算法框架，该算法框架能够有效、准确地提取变化区域，在实验区的总体精度达到92.90%，错分误差和漏分误差分别为4.47%和3.39%，相较于现有变化检测算法的精度有一定的提升。