CN111291615A - 一种多时相遥感影像变化监测方法 - Google Patents

Abstract

本发明公开了一种多时相遥感影像变化监测方法，通过将像元级变化检测方法和对象级变化检测方法进行集成，提取变化区域，通过顾及跳跃度的分级特征选择方法和跨时相的样本迁移学习方法识别变化区域的变化转移类型。本发明一种多时相遥感影像变化监测方法能够提高变化提取的自动化水平，得到的分类精度达90%以上，能够满足精细化变化监测的要求。

Description

一种多时相遥感影像变化监测方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种多时相遥感影像变化监测方法，可用于土地利用/土地覆被监测、地理信息更新、生态环境监测、城市规划研究。

背景技术

遥感影像变化检测技术通过数理模型分析能够得到不同时相上地表所发生的变化，是遥感影像变化监测的一种有效技术手段。目前已发展了较多的遥感影像变化检测方法，如影像代数法、影像变换法等方法，然而这些方法只能提取变化区域，无法识别变化区域具体的变化类型。分类后比较法不仅能够提取变化区域，而且能够获得变化区域具体的变化类型。采用分类后比较法识别变化类型过程中，为了减少对未变化区域地物的重复分类，按照“先变化提取后变化类型识别”的思路，首先提取变化区域，然后仅对变化区域的地物类型进行识别，是实际应用中一种有效的变化监测方法。

遥感变化检测方法从分析单元的角度可分为像元级变化检测方法和对象级变化检测方法。像元级变化检测方法是常采用的变化检测方法，具有无需人工干预、简单易行等优点，但像元级方法存在对影像配准、辐射校正误差敏感，检测结果椒盐现象严重等问题；对象级变化检测方法具有对影像配准、辐射校正要求低，改善椒盐现象等优点，但对象级方法受影像分割尺度等参数影响较大、实现过程繁琐。

另外，变化识别中利用高维特征能够提高不同类型的分离度，但将高维特征全部输入至分类器，容易因数据过大造成内存溢出、“维数灾难”问题，同时多时相影像分类时，单独选择训练样本进行分类，自动化程度低。

发明内容

针对上述问题，本发明提供了一种多时相遥感影像变化监测方法，通过将像元级变化检测方法和对象级变化检测方法进行集成，提取变化区域，通过顾及跳跃度的分级特征选择方法和跨时相的样本迁移学习方法识别变化区域的变化转移类型。该方法能够提高变化提取的自动化水平，得到的分类精度达90％以上，能够满足精细化变化监测的要求。

为实现上述目的，本发明提供了一种多时相遥感影像变化监测方法，具体包括如下步骤：

步骤1：获取T₁和T₂时相的遥感影像，包括光学影像和SAR影像，然后对影像进行预处理，具体是：

对光学影像进行图像配准、图像裁剪和大气校正；

对SAR影像进行图像配准、图像裁剪和去噪处理；

步骤2：在步骤1预处理后的光学影像和SAR影像上，分别提取T₁和T₂时相的如光谱特征、纹理特征、形态学特征、基于图像变换的特征以及植被指数等影像特征；

步骤3：根据步骤2提取的T₁和T₂时相影像特征，进行显著图引导的结合像元级与对象级变化检测过程，提取出变化区域，所述的结合像元级与对象级变化检测方法具体过程如下：

31、采用差值法对步骤2提取的T₁和T₂时相影像特征进行差值运算，得到特征差异影像；

32、对步骤31的特征差异影像进行显著图引导的像元级变化检测过程，具体方法如下：

321、利用最大对称环绕显著性检测(MSS)算法对特征差异影像进行显著性检测，生成显著图；

322、之后采用结合模糊C均值和马尔可夫随机场方法对显著区域进行分割得到初始像元级变化检测结果；

33、将步骤31特征差异影像叠加，进行面向对象的影像分割；

34、根据步骤33获得的影像对象，进行结合像元级与对象级变化检测过程，具体方法如下：341、通过构建对象的不确定性指数自适应选择训练样本；

342、之后基于训练样本采用随机森林分类器进行分类得到最终变化检测结果；

步骤4：引入跳跃度概念，采用顾及跳跃度的分级特征选择方法，依据特征重要性值对步骤2提取的T₁时相影像特征进行优化，得到最优特征子集；

步骤5：在步骤1获得的T₁时相影像上采集训练样本，利用随机森林分类器进行分类得到T₁时相的分类结果，将步骤3提取的变化区域结果与T₁时相的分类结果进行叠加，得到T₁时相变化区域的分类结果；

步骤6：根据步骤4获得的最优特征子集，得到T₂时相的最优特征子集，采用跨时相的样本迁移学习方法对T₂时相的变化区域进行分类，得到T₂时相的变化区域分类结果；

步骤7：对比T₁，T₂时相的变化区域分类结果，得到T₁、T₂时相间的变化转移矩阵。

进一步的，所述步骤32所述的显著图引导的像元级变化检测过程，其具体方法如下：采用MSS算法对特征差异影像进行显著性检测，得到显著图S(x_i)，然后采用Otsu算法对显著图进行分割，根据公式(1)提取显著性区域，

其中，S_E表示显著度；τ表示阈值；

当S(x_i)大于阈值τ时，表明像元x_i位于显著区域，即为可能发生变化的区域，显著度S_E用1表示；否则像元x_i位于非显著区域，即为未发生变化区域，显著度S_E用0表示；利用结合模糊C均值和马尔可夫随机场方法仅对显著区域进行处理，得到初始的变化检测结果。

进一步的，步骤3所述的结合像元级与对象级变化检测过程，其具体方法如下：将T₁和T₂时相的差异特征影像进行叠加，采用多尺度分割算法对叠加影像进行面向对象的分割，根据分割对象R_i中变化和未变化像元个数占整个对象R_i像元个数的比值，利用公式(2)计算分割对象的不确定性指数T，

其中，n_c，n_u和n分别为对象R_i中像元级变化检测得到的变化、未变化像元个数和对象总像元个数，设置不确定性阈值T_m，根据公式(3)判定对象R_i的属性l_i，

其中，l_i＝1,2,3表示对象R_i的属性为未变化对象、不确定是否发生变化的对象和变化对象，

当不确定性指数T>T_m时，分割对象为变化对象，被选为变化训练样本；

当不确定性指数T<-T_m时，分割对象为未变化对象，被选为未变化训练样本；

当不确定性指数位于区间-T_m≤T≤T_m时，分割对象为不确定是否发生变化的对象，不被选为训练样本；

之后利用选择的变化和未变化训练样本，采用随机森林分类方法对不确定是否发生变化的对象进行分类，得到最终的变化提取结果。

进一步的，步骤4中跳跃度的定义为；假设数据序列{X₁,X₂,X₃,…,X_n}是服从总体分布为F(x,θ)，个数为n的统计量，

为仅依赖数据序列{X₁,X₂,X₃,…,X_n}期望

的点估计，计算如公式(4)，那么

为

在k点的跳跃度，其中，

其中，x_i为序列第i个统计值，x_k为序列第k个统计值。

进一步的，步骤4中分级特征选择过程，其具体方法如下：

步骤41：利用随机森林算法，根据公式(5)计算影像特征的重要性，

其中，E_OOB1为决策树的袋外误差；

E_OOB2为将特征F_i替换后决策树的袋外误差；

N_tree为决策树个数；

步骤42：根据步骤41获得的特征重要性值，按照由低到高进行升序排序s₁<s₂<s₃<…<s_n，根据公式(6)和(7)计算每个特征的跳跃度t_k，

其中，s_i为第i个特征的重要性值，s_k为第k个特征的重要性值，n为特征总数；

步骤43：若t_k(k≥2)大于前面(k-1)个特征的跳跃度，则将s_k看作不同特征层的分界点，此时将前面(k-1)个特征作为特征层f₁，去掉f₁后的剩余特征重复进行以上操作，直到所有特征分级为止，从而得到最优特征子空间。

进一步的，步骤6所述的跨时相的样本迁移学习方法，具体过程如下：对T₁时相影像分类时，已采集训练样本，那么将T₁时相采集的训练样本中未发生变化的样本保留下来，转移作为T₂时相的训练样本；根据步骤3提取的变化区域，将步骤4获得的T₂时相最优特征子集输入至随机森林分类器，利用转移的训练样本进行训练，根据训练得到的分类模型对T₂时相上变化区域进行分类，得到T₂时相变化区域的分类结果。

与现有技术相比，本发明的有益效果在于：

本发明所提供的一种多时相遥感影像变化监测方法，具有三方面的优点：

(1)采用显著图引导的结合像元级与对象级变化检测方法进行变化提取，缓解了像元级变化检测结果的椒盐现象；由于将像元级方法与对象级方法有效串行，自适应选择训练样本，提高了变化提取的自动化水平；

(2)采用分级特征选择方法对提取的高维特征进行优化，为最优特征子集的选择提供了定量参考；

(3)采用跨时相的样本迁移学习方法利用转移的训练样本进行分类，得到的分类总体精度达90％以上，能够满足精细化变化监测的要求。

附图说明

图1：本发明基于多时相遥感影像变化监测方法的流程图；

图2：本发明实施例图像；

其中，2a为2013年ZY-3影像；2b为2018年ZY-3影像；2c为2013年TerraSAR-X影像；2d为2018年高分三号影像；

图3：本发明实例中本发明方法和传统方法的变化提取结果；

其中，3a为像元级方法；3b为对象级方法；3c为SG-PCAK方法；3d为SG-RCVA-RF方法；3e为结合像元级与对象级方法；3f为真实变化参考；

图4：本发明实例中本发明方法和传统方法的变化提取精度对比图；

其中，4a为ROC曲线图；4b为不同变化提取方法的精度对比图；

图5：本发明实施例中特征分级结果；

图6：本发明实施例中不同特征组合的分类精度对比图；

具体实施方式

以下结合附图和实施例详细说明本发明具体实施过程，需要说明的是，此处所描述的具体实施例仅用以解释本发明，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

实施例用2013年和2018年获取的江苏盐城滨海湿地TerraSAR-X、高分三号(GF3)SAR影像和资源三号(ZY-3)光学影像做实验，ZY-3影像包含4个光谱波段，TerraSAR-X影像和GF3影像仅使用了HH极化影像。数据基本信息如表1所示，影像数据如图2所示。

由图2a和图2b标注可知，实施例发生的变化是：海域转为盐沼植被(即米草盐沼向海域扩散，矩形框A所示)、海域转为养殖塘(即海域围垦为养殖塘，矩形框B所示)、农田转为养殖塘(矩形框C所示)、盐沼植被转为农田(矩形框D所示)和盐沼植被转为养殖塘(矩形框E所示)。

表1实验数据

图1是本发明提供的一种多时相遥感影像变化监测方法流程图，包括如下步骤：

步骤1：对2013年和2018年的SAR影像和光学影像进行预处理。对光学影像进行预处理，首先是图像配准，在两时相影像上分别选取一定数量的同名点进行自动配准，配准误差控制在1个像元之内，其次是图像裁剪，将两时相影像进行叠加，裁剪公共区域，最后利用FLAASH大气校正模块进行大气校正。对SAR影像进行预处理，除图像配准和图像裁剪外，采用非局域均值滤波算法对SAR影像进行去噪处理。

步骤2：在步骤1预处理的光学影像和SAR影像上，分别提取两个时相的4个光谱波段，1个亮度图像，1个SAR强度影像，240个基于灰度共生矩阵(GLCM)的纹理特征，100个形态学特征，15个基于非下采样Shearlet变换(NSST)的特征，10个植被指数，5个基于Sobel算子的边缘特征，共376个特征。具体如下：

a)纹理特征：窗口对GLCM纹理特征的影响最大，且不同窗口大小对纹理特征的表达存在差异，因而计算了3×3，7×7，11×11，15×15，19×19和23×23六种不同窗口大小的GLCM纹理特征。对每一个时相4个光谱波段和1幅SAR影像，分别计算6种不同窗口大小的8个纹理特征，共提取SAR与光学影像的GLCM纹理特征240(5×6×8)个。

b)形态学特征：提取形态学特征时，采用不同尺度的结构元素对不同地物的形态学特征进行描述，这里尺度取值为[1,2,3,4,5]。对每一个时相4个光谱波段和1幅SAR影像，分别计算5种尺度的形态学谱开重建和闭重建、差分形态学谱开重建和闭重建，共提取形态学特征100(5×5×4)个。

c)采用NSST变换对两个时相4个光谱波段和1幅SAR强度影像进行多尺度分解，每幅图像得到3幅不同尺度的特征图像，共提取15(3×5＝15)个图像变换特征。

d)植被指数：提取归一化差异植被指数(NDVI)、归一化差异水体指数(NDWI)、绿度比值、差值植被指数、比值植被指数、土壤调节植被指数、优化土壤调节植被指数、修正土壤调节植被指数、垂直植被指数、增强型植被指数，共10个植被指数。

e)边缘特征：4个光谱波段和1个SAR强度影像各提取1个Sobel边缘特征。

步骤3：发明人利用ZY-3影像的亮度、NDVI、NDWI三个特征，能很好地将实验例的变化信息体现出来，所以采用差值法对亮度、NDVI、NDWI特征进行差值运算，得到特征差异影像。不同特征的取值范围不同，为了使不同特征在后续处理中具有相同的地位，需进行归一化。假设F_i(i＝1,2)为提取的T_i(i＝1,2)时刻影像特征，根据公式(8)进行归一化，

根据公式(9)对归一化特征F_i ⁿ进行运算，生成特征差值图DF，表达式为

其中，

为归一化处理后的变化后影像特征，F₁ ⁿ为归一化处理后的变化前影像特征。

步骤3的结合模糊C均值和马尔可夫随机场方法对显著区域进行分类，具体方法如下：

传统马尔可夫随机场模型利用“0”和“1”表示像元x_i与其邻域N(i)的空间关系，1表示中心像元x_i的邻域像元x_j发生变化，0表示邻域像元x_j未发生变化，通过计算邻域像元x_j中发生变化和未变化的数量确定当前中心像元x_i的变化类型。当邻域像元x_j发生变化的个数大于未发生变化的个数时，当前中心像元x_i确定为发生变化像元，否则中心像元确定为未发生变化像元。由于地表覆盖的复杂性和影像分辨率的限制，图像上通常存在大量的混合像元，此时利用“0”和“1”不再能准确表达空间邻域关系，而且还容易造成空间邻域信息的过度利用。为此，利用空间引力模型改进的Potts模型来定义像元x_i的空间邻域关系，公式如(10)，(11)，

其中，x_j,j∈N(i){j＝1,2,…,8}表示中心像元x_i的3×3邻域窗口内的像元，z(p_i)和z(p_j)分别表示中心像元x_i与邻域像元x_j的类别标记，D_ij表示中心像元x_i与邻域像元x_j的欧式距离，p_i和p_j分别表示中心像元x_i和邻域像元x_j属于类别k,k∈{未变化,变化}的隶属度，是利用FCM算法计算得到。

图3和图4对比了结合像元级与对象级方法和传统方法的变化提取结果和变化提取精度。由图3和图4可知，发明方法的AUC值(0.86174)最大，像元级方法的AUC值(0.81123)次之，SG-PCAK方法的AUC值(0.75266)最低。

由图4b可知，结合像元级与对象级方法的变化提取总体精度(OA)、精度比(AR)最高，总体精度为93.51％，高于SG-PCAK方法5.95％，同时虚检率(FDR)最低，与SG-PCAK方法相比，虚检率降低了35.96％，精度比提高了29.24％。

与SG-RCVA-RF方法相比，精度比提高了14.23％，与像元级方法和对象级方法相比，虚检率分别降低了29.04％和22.78％。AUC、OA、AR和FDR四个指标分析结果均表明，结合像元级与对象级方法的变化提取结果优于像元级、对象级、SG-PCAK方法和SG-RCVA-RF方法。

SG-PCAK方法具体参见文献“Zheng Y G,Jiao L C,Liu H Y,et al.Unsupervisedsaliency-guided SAR image change detection.Pattern Recognition,2017,61:309-326”。

SG-RCVA-RF方法具体参见文献“Feng W Q,Sui H G,Chen X.Saliency-guidedchange detection of remotely sensed images using random forest.TheInternational Archives of the Photogrammetry,Remote Sensing and SpatialInformation Sciences,2018,XLII-3:341-348”。

步骤4：利用跳跃度对2013年提取的376个特征进行特征分级，得到最优特征子集空间。

采用特征分级方法将提取的376个特征分为10级，如图5所示，F＝{f₁,f₂,f₃,f₄,f₅,f₆,f₇,f₈,f₉,f₁₀}，其中，特征子集的特征重要性由低到高排序分别如下：f₁＝[1,102]，f₂＝[103,136]，f₃＝[137,158]，f₄＝[159,177]，f₅＝[178,194]，f₆＝[195,254]，f₇＝[255,292]，f₈＝[293,314]，f₉＝[315,348]，f₁₀＝[349,376]。通过分析不同的特征组合与分类精度的关系，见图6所示，可见特征组合F₂＝{f₃,f₄,f₅,f₆,f₇,f₈,f₉,f₁₀}的分类精度最高，总体精度为97.39％。可见特征组合F₂＝{f₃,f₄,f₅,f₆,f₇,f₈,f₉,f₁₀}为分级特征选择方法得到的最优特征子空间。

步骤5：在2013年影像上采集训练样本，利用随机森林分类器进行分类得到2013年的分类结果。

步骤6：依据步骤3提取出的变化区域，将2013年采集的训练样本中未发生变化的样本保留下来，转移作为2018年的训练样本。将2018年影像上最优特征子集输入至随机森林分类器，利用转移的训练样本进行训练，仅对2018年变化区域进行分类，得到2018年时相变化区域的分类结果；

步骤7：对比2013年和2018年变化区域的分类结果，得到2013-2018年间变化转移矩阵，如表2所示。

表2 2013-2018年变化转移矩阵(单位:km²)

Claims (7)

1.一种多时相遥感影像变化监测方法，其特征在于，包括以下步骤：

步骤1：获取T₁和T₂时相的遥感影像，包括光学影像和SAR影像，然后对影像进行预处理，具体是：

对光学影像进行图像配准、图像裁剪和大气校正；

对SAR影像进行图像配准、图像裁剪和去噪处理；

步骤2：在步骤1预处理后的光学影像和SAR影像上，分别提取T₁和T₂时相的影像特征；

步骤3：根据步骤2提取的T₁和T₂时相影像特征，进行显著图引导的结合像元级与对象级变化检测过程，提取出变化区域，所述的结合像元级与对象级变化检测方法具体过程如下：

31、采用差值法对步骤2提取的T₁和T₂时相影像特征进行差值运算，得到特征差异影像；

32、对步骤31的特征差异影像进行显著图引导的像元级变化检测过程，具体方法如下：

321、利用最大对称环绕显著性检测(MSS)算法对特征差异影像进行显著性检测，生成显著图；

322、之后采用结合模糊C均值和马尔可夫随机场方法对显著区域进行分割得到初始像元级变化检测结果；

33、将步骤31特征差异影像叠加，进行面向对象的影像分割；

34、根据步骤33获得的影像对象，进行结合像元级与对象级变化检测过程，

具体方法如下：341、通过构建对象的不确定性指数自适应选择训练样本；

342、之后基于训练样本采用随机森林分类器进行分类得到最终变化检测结果；

步骤4：引入跳跃度概念，采用顾及跳跃度的分级特征选择方法，依据特征重要性值对步骤2提取的T₁时相影像特征进行优化，得到最优特征子集；

步骤5：在步骤1获得的T₁时相影像上采集训练样本，利用随机森林分类器进行分类得到T₁时相的分类结果，将步骤3提取的变化区域结果与T₁时相的分类结果进行叠加，得到T₁时相变化区域的分类结果；

步骤6：根据步骤4获得的最优特征子集，得到T₂时相的最优特征子集，将步骤3提取的变化区域结果与T₂时相影像进行叠加得到T₂时相的变化区域，采用跨时相的样本迁移学习方法对T₂时相的变化区域进行分类，得到T₂时相的变化区域分类结果；

步骤7：对比T₁，T₂时相的变化区域分类结果，得到T₁、T₂时相间的变化转移矩阵。

2.如权利要求1所述的多时相遥感影像变化监测方法，其特征在于：步骤2中的所述影像特征包括光谱特征、纹理特征、形态学特征、基于图像变换的特征以及植被指数。

3.如权利要求1所述的多时相遥感影像变化监测方法，其特征在于：

所述步骤32所述的显著图引导的像元级变化检测过程，其具体方法如下：采用MSS算法对特征差异影像进行显著性检测，得到显著图S(x_i)，然后采用Otsu算法对显著图进行分割，根据公式(1)提取显著性区域，

其中，S_E表示显著度；τ表示阈值；

当S(x_i)大于阈值τ时，表明像元x_i位于显著区域，即为可能发生变化的区域，显著度S_E用1表示；否则像元x_i位于非显著区域，即为未发生变化区域，显著度S_E用0表示；利用结合模糊C均值和马尔可夫随机场方法仅对显著区域进行处理，得到初始的变化检测结果。

4.如权利要求1所述的多时相遥感影像变化监测方法，其特征在于：步骤3所述的结合像元级与对象级变化检测过程，其具体方法如下：将T₁和T₂时相的差异特征影像进行叠加，采用多尺度分割算法对叠加影像进行面向对象的分割，根据分割对象R_i中变化和未变化像元个数占整个对象R_i像元个数的比值，利用公式(2)计算分割对象的不确定性指数T，

其中，n_c，n_u和n分别为对象R_i中像元级变化检测得到的变化、未变化像元个数和对象总像元个数，设置不确定性阈值T_m，根据公式(3)判定对象R_i的属性l_i，

其中，l_i＝1,2,3表示对象R_i的属性为未变化对象、不确定是否发生变化的对象和变化对象，

当不确定性指数T>T_m时，分割对象为变化对象，被选为变化训练样本；

当不确定性指数T<-T_m时，分割对象为未变化对象，被选为未变化训练样本；

当不确定性指数位于区间-T_m≤T≤T_m时，分割对象为不确定是否发生变化的对象，不被选为训练样本；

之后利用选择的变化和未变化训练样本，采用随机森林分类方法对不确定是否发生变化的对象进行分类，得到最终的变化提取结果。

5.如权利要求1所述的多时相遥感影像变化监测方法，其特征在于：步骤4中跳跃度的定义为；假设数据序列{X₁,X₂,X₃,…,X_n}是服从总体分布为F(x,θ)，个数为n的统计量，

为仅依赖数据序列{X₁,X₂,X₃,…,X_n}期望

的点估计，计算如公式(4)，那么

为

在k点的跳跃度，其中，

其中，x_i为序列第i个统计值，x_k为序列第k个统计值。

6.如权利要求1所述的多时相遥感影像变化监测方法，其特征在于：步骤4中分级特征选择过程，其具体方法如下：

步骤41：利用随机森林算法，根据公式(5)计算影像特征的重要性VI，

其中，E_OOB1为决策树的袋外误差；

E_OOB2为将特征F_i替换后决策树的袋外误差；

N_tree为决策树个数；

步骤42：根据步骤41获得的特征重要性值，按照由低到高进行升序排序s₁<s₂<s₃<…<s_n，根据公式(6)和(7)计算每个特征的跳跃度t_k，

其中，s_i为第i个特征的重要性值，s_k为第k个特征的重要性值，n为特征总数；

步骤43：若t_k(k≥2)大于前面(k-1)个特征的跳跃度，则将s_k看作不同特征层的分界点，此时将前面(k-1)个特征作为特征层f₁，去掉f₁后的剩余特征重复进行以上操作，直到所有特征分级为止，从而得到最优特征子空间。

7.如权利要求1所述的多时相遥感影像变化监测方法，其特征在于：步骤6所述的跨时相的样本迁移学习方法，具体过程如下：对T₁时相影像分类时，已采集训练样本，那么将T₁时相采集的训练样本中未发生变化的样本保留下来，转移作为T₂时相的训练样本；根据步骤3提取的变化区域，将步骤4获得的T₂时相最优特征子集输入至随机森林分类器，利用转移的训练样本进行训练，根据训练得到的分类模型对T₂时相上变化区域进行分类，得到T₂时相变化区域的分类结果。

Images