CN114255397A - 基于空间相关性的无监督光谱水域变化检测方法及装置 - Google Patents

Abstract

本发明公开了一种基于空间相关性的无监督光谱水域变化检测方法和装置，采用对双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；对引导滤波后的双时相一维特征图像进行自抑制；基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。本方法具有不需要提前得到图像先验信息，具有一定的普适性。

Description

基于空间相关性的无监督光谱水域变化检测方法及装置

技术领域

本发明涉及遥感卫星地面/星载图像处理领域，具体涉及一种基于空间相关性的无监督光谱水域变化检测方法及装置。

背景技术

近年来，全球变暖导致气候变化异常，地震，海啸等自然灾害时有发生，这些灾害具有突发性且破坏范围大，对人们的生命和财产安全造成的损失难以估量。常规的灾情获取方法如实地调查研究一般存在成本高、周期长、调查范围有限等缺陷。而遥感技术作为一种直观非接触的侦测技术，可以对灾区进行大范围的快速检测，获取灾情信息，具有成本低，周期短，能在短时间内实现对灾难的定量评估。利用图像变化检测，提取同一成像场景中不同时期的图像之间发生的变化信息，是获取灾害信息、检测灾区变化的有效手段。

目前针对光谱图像的变化检测算法基本上可以分为四类：基于图像代数的方法、基于图像变换的方法、基于图像分类的方法和基于深度学习的方法。其中基于图像分类的方法和深度学习的方法属于有监督学习方法，需要提前得到部分图像的先验信息，用于对分类器进行训练，或者采用阈值对图像进行分类。如果采用阈值法，阈值设定往往是由研究人员多次试验确定，方法的可迁移性差。传统基于图像代数的方法处理过程比较简单，广泛应用的变化检测模型是基于图像代数的F.Bovolo等人提出的变化向量分析(CVA)方法，但如果图像的波段数量增大，这种方法的检测计算量大幅增加，检测效能下降。

发明内容

有鉴于此，本发明提供了一种基于空间相关性的无监督光谱水域变化检测方法及装置，能够解决现有利用遥感光谱数据变化进行检测的可迁移性差、计算量大、检测效能低的技术问题。

为了解决上述技术问题，本发明是这样实现的。

一种基于空间相关性的无监督光谱水域变化检测方法，包括：

步骤S1：获取双时相遥感光谱图像，所述双时相遥感光谱图像是对同一区域在两个时间点获取的遥感光谱图像；对所述双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；

步骤S2：对所述双时相地表反射率图像进行主成分分析(Principal ComponentAnalysis，PCA)，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；

步骤S3：基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；

步骤S4：对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；

步骤S5：对引导滤波后的双时相一维特征图像进行自抑制；

步骤S6：基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；

步骤S7：对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。

优选地，所述步骤S2：对所述双时相地表反射率图像进行主成分分析，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像，其中，对所述双时相地表反射率图像的每一个图像，执行以下操作：

步骤S201：获取地表反射率图像信息，记为

，其中地表反射率图像，表示一个大小为M×N×l的矩阵，M，N为地表反射率图像中光谱数据的空间维度，l为光谱数据包含的波段的维度；

步骤S202：将所述地表反射率图像转化为维度为l×d的观测样本数据矩阵X＝[x₁,x₂,x₃,...,x_l]，其中d＝M×N是一维像元的总个数；

步骤S203：计算观测样本数据集X各波段的均值μ和协方差矩阵Σ：

步骤S204：通过求解特征方程来获得协方差矩阵Σ的特征值λ和特征向量U：

(λI-Σ)＝0，其中I为单位矩阵；

步骤S205：对所有的特征值进行降序排列，即λ₁≥λ₂...≥λ_l。λ₁，λ₂...，λ_l对应的特征向量分别为U₁,U₂...,U_l，得到PCA变换矩阵W＝[U₁,U₂...,U_l]，其中W为向量矩阵，满足条件WW^T＝W^TW＝I；

步骤S206：计算Y＝[Y₁,Y₂,...Y_l]＝W^TX，选取Y的第一个分量Y₁，作为地表反射率图像的一维特征图像。

优选地，所述步骤S3：基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作，其中：

确定属性滤波器为面积属性；设置阈值τ，所述形态学属性滤波开操作用于去除所述双时相一维特征图像中面积小于阈值τ的连通分量的亮度大于周边像素的部分；定义保留的连通分量为：

γ^T(f)(x)＝max{k:x∈Γ^T[Th_k(f)]}

其中，f为待操作的图像，Th_k(f)为通过灰度级k对图像进行阈值操作而获得的二值图像，Γ^T[Th_k(f)]为二值图像中值为1的像素点的集合，x为集合Γ^T[Th_k(f)]中的一个，γ^T(f)(x)为图像保留的连通分量。

优选地，所述步骤S4：对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波，包括：

对所述经过形态学属性滤波开操作的所述双时相一维特征图像中的每一个图像，执行以下操作：

确定引导图像，所述引导图像为原始的一维特征图像本身；

引导滤波的过程为：

其中q_m为引导滤波后的输出图像、I为引导图像，m，n分别表示所述的引导滤波后的输出图像和输入的一维特征图像中像素点的索引，p_n为所述经过形态学属性滤波开操作的作为输入的一维特征图像。

优选地，所述步骤S5：对引导滤波后的双时相一维特征图像进行自抑制，其中，对所述引导滤波后的双时相一维特征图像中的每一个图像，执行以下操作：

步骤S501：记所述引导滤波后的一维特征图像为z；设置迭代次数Num，抑制系数

，设置对应于z大小的全1矩阵G；初始化当前迭代次数num1为0；

步骤S502：计算z＝z.*G，

其中，β为抑制系数，z.表示对矩阵中的每个值进行与G矩阵中对应位置的值进行相乘的运算；将当前迭代次数num1赋值为num1加1；若当前迭代次数num1<Num，进入步骤S502；否则，进入步骤S503；

步骤S503：将最终得到的矩阵z作为自抑制后的图像。

优选地，所述步骤S6：基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理，其中，对所述自抑制后的双时相一维特征图像中的每一个图像，执行以下操作：假设阈值为t，w₀和w₁分别是一个像素属于前景、背景的概率，σ₁和σ₂分别表示前景类和后景类的方差；在给定的阈值t的情况下，w₀和w₁分别定义为：

其中，p(i)为不同像素等级下的直方图统计值，直方图有L个等级；

将阈值t在自抑制后的一维特征图像的灰度值最小值和最大值之间进行遍历操作；若满足σ_w(t)得到最大值，则此时t为最佳分割阈值；

基于最佳分割阈值，对自抑制后的一维特征图像进行二值化处理。

本发明所提供的一种基于空间相关性的无监督光谱水域变化检测装置，包括预处理模块：配置为获取双时相遥感光谱图像，所述双时相遥感光谱图像是对同一区域在两个时间点获取的遥感光谱图像；对所述双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；

降维模块：配置为对所述双时相地表反射率图像进行主成分分析(PrincipalComponent Analysis，PCA)，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；

形态滤波模块：配置为基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；

引导滤波模块：配置为对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；

自抑制模块：配置为对引导滤波后的双时相一维特征图像进行自抑制；

二值化模块：配置为基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；

差分模块：配置为对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。

有益效果：

本发明利用遥感影像光谱和空间相关性信息，提取两个时相遥感影像的变化信息，属于一种无监督方法，不需要提前得到图像先验信息，也不需要通过试验得到相关阈值，具有一定的普适性。

本发明先对光谱图像进行降维处理，不仅可用于多光谱数据，也适用于超/高光谱类高维数据的处理。

本发明利用形态学属性滤波有效较低变化检测结果中的噪声干扰，同时利用引导滤波边缘保持和自抑制算法提升检测精度。

基于总体精度的评价指标表明本发明方法优于传统的变化向量分析(CVA)、主成分分析-变化向量分析(PCA-CVA)、支持向量机(SVM)、极限学习机(ELM)。

附图说明

图1为基于空间相关性的无监督光谱水域变化检测方法流程示意图；

图2为遥感光谱影像预处理流程示意图；

图3为基于空间相关性的无监督光谱水域变化检测算法流程示意图；

图4为形态学属性开操作原理示意图；

图5为抑制函数示意图；

图6(A)为T1时刻高光谱图像；图6(B)为T2时刻高光谱图像；图6(C)为地面真值图；

图7(A)为T1时刻高光谱图经过主成分分析(PCA)后得到的一维图像；图7(B)为T2时刻高光谱图经过主成分分析(PCA)后得到的一维图像；图7(C)为PCA输出的特征图经过属性滤波的空间相关性优化之后的结果图；图7(D)为PCA输出的特征图经过属性滤波的空间相关性优化之后的结果图；图7(E)、7(F)为空间优化后的两时相图像先经过自抑制算法加强变化区域和非变化区域的对比度；图7(G)为最终检测结果；

图8(A)为CVA算法的检测结果；图8(B)为PCA-CVA算法的检测结果；图8(C)为SVM算法的检测结果；图8(D)为ELM算法的检测结果；图8(E)为本发明实施例的检测结果；图8(F)为地面真值数据；

图9为基于空间相关性的无监督光谱水域变化检测装置结构示意图。

具体实施方式

下面结合附图和实施例，对本发明进行详细描述。

本发明一种基于空间相关性的无监督光谱水域变化检测方法，该方法解决现有利用遥感光谱数据变化检测技术不足，通过将光谱特征和空间特征进行结合，有效提取遥感图像中的水体变化信息。利用主成分分析方法解决超/高光谱数据的数据冗余问题；利用形态学属性滤波有效较低变化检测结果中的噪声干扰，同时利用引导滤波边缘保持和自抑制算法提升检测精度。

如图1、3所示，本发明提供的基于空间相关性的无监督光谱水域变化检测方法，包括如下步骤：

步骤S1：获取双时相遥感光谱图像，所述双时相遥感光谱图像是对同一区域在两个时间点获取的遥感光谱图像；对所述双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；

步骤S2：对所述双时相地表反射率图像进行主成分分析(Principal ComponentAnalysis，PCA)，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；

步骤S3：基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；

步骤S4：对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；

步骤S5：对引导滤波后的双时相一维特征图像进行自抑制；

步骤S6：基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；

步骤S7：对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。

本实施例将双时相地表反射率图像分别输入主成分分析方法，降至一维，接着根据空间相关性对其进行形态学属性滤波开操作，然后使用引导滤波强化变化区域的边缘特征，再采用自抑制算法自抑制算法使得图像中明亮的部分与黑暗的部分对比更加明显，即提升目标与背景的对比度。然后利用OTSU阈值算法进行图像的二值化，最后，将两个时刻的图像差分，得到最终的变化检测结果。

所述步骤S1，所述双时相遥感光谱图像是对同一区域在两个时间点获取的感光谱图像，该图像是航天器通过光谱相机获得的，经过底噪去除、相对辐射校正、坏像元恢复，几何校正等处理得到的一级(L1)光谱图像数据。在一级光谱图像数据的基础上进行后续的应用处理工作，需要对所述一级光谱图像数据进行预处理。

所述预处理包括去除所述双时相遥感光谱图像中未校准谱段、噪声谱段，对所述双时相遥感光谱图像进行绝对辐射校正和大气校正；预处理流程如图2所示。超/高光谱影像，有些频段可能未经过校准或噪声比较大，例如受水汽吸收影响比较严重的波段，这些波段将不参与后续的应用处理过程，以免对最终应用的效果带来影响。需要将这些频段从光谱影像中去除。利用ENVI的FLAASH对光谱影像进行大气校正，得到地表反射率影像。

本实施例中，大气校正的方法可以使用ENVI自带的FLAASH或者其他软件代码，目的是去除大气对光谱信息的影响，得到地表的固有光谱特性。

所述步骤S2：对所述双时相地表反射率图像进行主成分分析，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像，其中，对所述双时相地表反射率图像的每一个图像，执行以下操作：

步骤S201：获取地表反射率图像信息，记为

其中地表反射率图像，表示一个大小为M×N×l的矩阵，M，N为地表反射率图像中光谱数据的空间维度，l为光谱数据包含的波段的维度；

步骤S202：将所述地表反射率图像转化为维度为l×d的观测样本数据矩阵X＝[x₁,x₂,x₃,...,x_l]，其中d＝M×N是一维像元的总个数；

步骤S203：计算观测样本数据集X各波段的均值μ和协方差矩阵Σ：

步骤S204：通过求解特征方程来获得协方差矩阵Σ的特征值λ和特征向量U：(λI-Σ)＝0，其中I为单位矩阵；

步骤S205：对所有的特征值进行降序排列，即λ₁≥λ₂...≥λ_l。λ₁，λ₂...，λ_l对应的特征向量分别为U₁,U₂...,U_l，得到PCA变换矩阵W＝[U₁,U₂...,U_l]，其中W为向量矩阵，满足条件WW^T＝W^TW＝I；

步骤S206：计算Y＝[Y₁,Y₂,...Y_l]＝W^TX，选取Y的第一个分量Y₁，作为地表反射率图像的一维特征图像。

本实施例中，利用主成分分析方法将多光谱特征图像降低为一维有效特征图像，以便后续经典滤波器对特征图像的空间滤波处理得到最优效果。PCA变换可以表示为Y＝[Y₁,Y₂,...Y_l]＝W^TX，并且Y的第一个分量的方差是最大的，第二个分量的方差次之，排在前面的结果主要成分包含了数据的大部分信息。选取Y的第一个分量，即将光谱特征信息压缩到一维特征图像，有助于开展后续的工作。

所述步骤S3：基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作，其中：

确定属性滤波器为面积属性；设置阈值τ，所述形态学属性滤波开操作用于去除所述双时相一维特征图像中面积小于阈值τ的连通分量的亮度大于周边像素的部分；定义保留的连通分量为：

γ^T(f)(x)＝max{k:x∈Γ^T[Th_k(f)]}

其中，f为待操作的图像，Th_k(f)为通过灰度级k对图像进行阈值操作而获得的二值图像，Γ^T[Th_k(f)]为二值图像中值为1的像素点的集合，x为集合Γ^T[Th_k(f)]中的一个，γ^T(f)(x)为图像保留的连通分量。

形态学属性滤波的原理如图4所示，如图4，例如假设属性滤波器设置的阈值τ为17，针对图像中某一中心点及周围像素点总数不超过17的区域，假设某一点的像素值为32，周围点的像素值为11，通过设置门限分别为0～255，大于门限值的置为1，小于门限值的置为0，可以得到256幅图像，可以发现，在门限为12时，像素值为11的点已经全部置为0，而像素值为32的点仍为1，此时就将仍为1的点的原像素置更新为此时的门限值，通过这种更新将原来图像中明亮的部分，即亮度大于周边像素的部分去除。

本实施例中，采用形态学属性滤波开操作去除图像中可能为噪声的小面积区域。根据图像中目标区域大小的不同，预先设定阈值。

所述步骤S4：对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波，包括：

对所述经过形态学属性滤波开操作的所述双时相一维特征图像中的每一个图像，执行以下操作：

确定引导图像，所述引导图像为原始的一维特征图像本身；

引导滤波的过程为：

其中q_m为引导滤波后的输出图像、I为引导图像，m，n分别表示所述的引导滤波后的输出图像和输入的一维特征图像中像素点的索引，p_n为所述经过形态学属性滤波开操作的作为输入的一维特征图像。

本实施例中，利用引导滤波对图像边缘内部的特征进行平滑，同时保持边界，有利于后面检测精度的提升。将引导图像选择为输入的初始图像，实现边缘保持滤波，以实现图像的边缘保持。

所述步骤S5：对引导滤波后的双时相一维特征图像进行自抑制，其中，对所述引导滤波后的双时相一维特征图像中的每一个图像，执行以下操作：

步骤S501：记所述引导滤波后的一维特征图像为z；设置迭代次数Num，抑制系数

，设置对应于z大小的全1矩阵G；初始化当前迭代次数num1为0；

步骤S502：计算z＝z.*G，

其中，β为抑制系数，z.表示对矩阵中的每个值进行与G矩阵中对应位置的值进行相乘的运算；将当前迭代次数num1赋值为num1加1；若当前迭代次数num1<Num，进入步骤S502；否则，进入步骤S503；

步骤S503：将最终得到的矩阵z作为自抑制后的图像。

本实施例中，设计自抑制算法的目的是使得图像中明亮的部分与黑暗的部分对比更加明显，即提升目标与背景的对比度，将有利于提取图像中具有区分性的变化特征。自抑制算法的公式

的函数图像如图5所示。输入的灰度图像的像素值越接近于1的像素点，在经过多次迭代后像素值会变小，但变小速率由斜率来看较慢，而对应于像素值越接近于0的像素点经过相同次迭代后像素值也会变小，但速率较快，使得图像中明暗对比更加明显。通过步骤S5的处理，实现图像对比度的提升。

所述步骤S6：基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理，OTSU阈值方法通过统计学方法自动选取二值化阈值度图像进行二值化处理。整个图像的灰度值被划分为前景色和背景色，阈值的选取使得两者尽可能的分开。最优阈值的判据为最大组间方差或最小组内方差，可以证明两者等价。

对所述自抑制后的双时相一维特征图像中的每一个图像，执行以下操作：假设阈值为t，w₀和w₁分别是一个像素属于前景、背景的概率，σ₁和σ₂分别表示前景类和后景类的方差；在给定的阈值t的情况下，w₀和w₁分别定义为：

其中，p(i)为不同像素等级下的直方图统计值，图像的直方图有L个等级。

本实施例中，自适应阈值方法使用OTSU阈值方法，通过统计学方法自动选取二值化阈值。

将阈值t在S5步骤得到的自抑制后的图像的灰度值最小值和最大值之间进行遍历操作；若满足σ_w(t)得到最大值，则此时t为最佳分割阈值；

基于最佳分割阈值，对自抑制后的一维特征图像进行二值化处理。

本实施例中，一般L＝256。

进一步地，本实施例中，采用总体精度(OA)对检测结果进行评价。采用总体精度(OA)作为评价指标。评价指标包括以下指标：

(1)真正例(TP)，表示正确检测出的发生变化的像素个数；

(2)真负例(TN)，表示正确检测出的未发生变化的像素个数；

(3)假正例(FP)，表示未发生变化的像素被错误检测为发生变化的像素个数；

(4)假负例(FN)，表示发生变化的像素被错误检测为未发生变化的像素个数。

总体精度(OA)是指被正确分类的类别的像素个数与总的类别个数的比值，在变化检测中发生变化和未发生变化可以看作二分类检测，因此，OA的定义为：

以下结合另一实施例说明本发明的方法。

采用一个真实的应用于变化检测的高光谱数据集River。这个数据集是由高光谱传感器采集的，该高光谱传感器包括从400nm至2500nm范围内的242个光谱波段，它的光谱分辨率大概为10nm，空间分辨率为30m。选取的高光谱元数据包含3幅图像，一幅为T1时刻的高光谱图像；一幅为T2时刻的高光谱图像，与T1时刻为同一拍摄场景；一幅为地面真值图，反映了同一地点两个不同时刻的变化信息，地面实况图为二值图像，白色部分代表发生过变化的区域，黑色部分代表未发生变化的区域。该数据集一共有198个波段，每个波段包括463*241个像素值。从地面实况图来看，图像中变化的部分主要集中在河流部分。三幅图如图6所示，图6(A)为T1时刻高光谱图像，图6(B)为T2时刻高光谱图像，图6(C)为地面真值图像。

对上述两个时刻的数据进行本文的基于空间相关性的无监督水域变化检测，检测结果如图7所示。图7(A)和图7(B)分别为T1和T2两个时刻高光谱图像经过主成分分析(PCA)后得到的一维图像；图7(C)和图7(D)分别为PCA输出的特征图经过属性滤波和引导滤波的空间相关性优化之后的结果图；

图7(E)和图7(F)是上一步空间优化后的两时相图像先经过自抑制算法加强变化区域和非变化区域的对比度，然后通过自适应阈值算法二值化后的结果，图7(G)是最终检测结果。

将本方法与四个传统变化检测算法进行对比，传统方法包括两个无监督算法：变化向量分析(CVA)算法、主成分分析-变化向量分析(PCA-CVA)算法；两个有监督算法：支持向量机(SVM)算法、极限学习机(ELM)算法。图8为利用这些对比方法变化检测的结果。图8(A)是CVA算法的检测结果，图8(B)是PCA-CVA算法的检测结果，图8(C)是SVM算法的检测结果，图8(D)是ELM算法的检测结果，图8(E)是本变化检测方法的检测结果，图8(F)为地面真值数据。从对比图中可以看出利用本实施例的方法，可以有效降低错检率。

定量化的对比结果如下表1所示：

表1变化检测算法检测性能对比结果

有表1可见，本实施例的方法在OA指标上获得了最优的效果。

本发明还提供了一种基于空间相关性的无监督光谱水域变化检测装置，如图9所示，该装置包括预处理模块：配置为获取双时相遥感光谱图像，所述双时相遥感光谱图像是对同一区域在两个时间点获取的遥感光谱图像；对所述双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；

降维模块：配置为对所述双时相地表反射率图像进行主成分分析(PrincipalComponent Analysis，PCA)，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；

形态滤波模块：配置为基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；

引导滤波模块：配置为对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；

自抑制模块：配置为对引导滤波后的双时相一维特征图像进行自抑制；

二值化模块：配置为基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；

差分模块：配置为对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。

以上的具体实施例仅描述了本发明的设计原理，该描述中的部件形状，名称可以不同，不受限制。所以，本发明领域的技术人员可以对前述实施例记载的技术方案进行修改或等同替换；而这些修改和替换未脱离本发明创造宗旨和技术方案，均应属于本发明的保护范围。

Claims (7)

1.一种基于空间相关性的无监督光谱水域变化检测方法，其特征在于，包括：

步骤S1：获取双时相遥感光谱图像，所述双时相遥感光谱图像是对同一区域在两个时间点获取的遥感光谱图像；对所述双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；

步骤S2：对所述双时相地表反射率图像进行主成分分析(Principal ComponentAnalysis，PCA)，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；

步骤S3：基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；

步骤S4：对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；

步骤S5：对引导滤波后的双时相一维特征图像进行自抑制；

步骤S6：基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；

步骤S7：对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。

2.如权利要求1所述的方法，其特征在于，所述步骤S2：对所述双时相地表反射率图像进行主成分分析，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像，其中，对所述双时相地表反射率图像的每一个图像，执行以下操作：

步骤S201：获取地表反射率图像信息，记为

其中地表反射率图像，表示一个大小为M×N×l的矩阵，M，N为地表反射率图像中光谱数据的空间维度，l为光谱数据包含的波段的维度；

步骤S202：将所述地表反射率图像转化为维度为l×d的观测样本数据矩阵X＝[x₁,x₂,x₃,...,x_l]，其中d＝M×N是一维像元的总个数；

步骤S203：计算观测样本数据集X各波段的均值μ和协方差矩阵Σ：

步骤S204：通过求解特征方程来获得协方差矩阵Σ的特征值λ和特征向量U：(λI-Σ)＝0，其中I为单位矩阵；

步骤S205：对所有的特征值进行降序排列，即λ₁≥λ₂...≥λ_l；λ₁，λ₂...，λ_l对应的特征向量分别为U₁,U₂...,U_l，得到PCA变换矩阵W＝[U₁,U₂...,U_l]，其中W为向量矩阵，满足条件WW^T＝W^TW＝I；

步骤S206：计算Y＝[Y₁,Y₂,...Y_l]＝W^TX，选取Y的第一个分量Y₁，作为地表反射率图像的一维特征图像。

3.如权利要求1所述的方法，其特征在于，所述步骤S3：基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作，其中：

确定属性滤波器为面积属性；设置阈值τ，所述形态学属性滤波开操作用于去除所述双时相一维特征图像中面积小于阈值τ的连通分量的亮度大于周边像素的部分；定义保留的连通分量为：

γ^T(f)(x)＝max{k:x∈Γ^T[Th_k(f)]}

其中，f为待操作的图像，Th_k(f)为通过灰度级k对图像进行阈值操作而获得的二值图像，Γ^T[Th_k(f)]为二值图像中值为1的像素点的集合，x为集合Γ^T[Th_k(f)]中的一个，γ^T(f)(x)为图像保留的连通分量。

4.如权利要求1所述的方法，其特征在于，所述步骤S4：对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波，包括：

对所述经过形态学属性滤波开操作的所述双时相一维特征图像中的每一个图像，执行以下操作：

确定引导图像，所述引导图像为原始的一维特征图像本身；

引导滤波的过程为：

其中q_m为引导滤波后的输出图像、I为引导图像，m，n分别表示所述的引导滤波后的输出图像和输入的一维特征图像中像素点的索引，p_n为所述经过形态学属性滤波开操作的作为输入的一维特征图像。

5.如权利要求1所述的方法，其特征在于，所述步骤S5：对引导滤波后的双时相一维特征图像进行自抑制，其中，对所述引导滤波后的双时相一维特征图像中的每一个图像，执行以下操作：

步骤S501：记引导滤波后的一维特征图像为z；设置迭代次数Num，抑制系数

设置对应于z大小的全1矩阵G；初始化当前迭代次数num1为0；

步骤S502：计算z＝z.*G，

其中，β为抑制系数，z.表示对矩阵中的每个值进行与G矩阵中对应位置的值进行相乘的运算；将当前迭代次数num1赋值为num1加1；若当前迭代次数num1<Num，进入步骤S502；否则，进入步骤S503；

步骤S503：将最终得到的矩阵z作为自抑制后的图像。

6.如权利要求1所述的方法，其特征在于，所述步骤S6：基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理，其中，对所述自抑制后的双时相一维特征图像中的每一个图像，执行以下操作：假设阈值为t，w₀和w₁分别是一个像素属于前景、背景的概率，σ₁和σ₂分别表示前景类和后景类的方差；在给定的阈值t的情况下，w₀和w₁分别定义为：

其中，p(i)为不同像素等级下的直方图统计值，直方图有L个等级；

将阈值t在自抑制后的一维特征图像的灰度值最小值和最大值之间进行遍历操作；若满足σ_w(t)得到最大值，则此时t为最佳分割阈值；

基于最佳分割阈值，对自抑制后的一维特征图像进行二值化处理。

7.一种基于空间相关性的无监督光谱水域变化检测装置，其特征在于，包括：

预处理模块：配置为获取双时相遥感光谱图像，所述双时相遥感光谱图像是对同一区域在两个时间点获取的遥感光谱图像；对所述双时相遥感光谱图像进行预处理，得到双时相地表反射率图像；

降维模块：配置为对所述双时相地表反射率图像进行主成分分析(PrincipalComponent Analysis，PCA)，将所述双时相地表反射率图像的光谱维度降为一维，得到双时相一维特征图像；

形态滤波模块：配置为基于空间相关性对所述双时相一维特征图像进行形态学属性滤波开操作；

引导滤波模块：配置为对经过形态学属性滤波开操作的所述双时相一维特征图像进行引导滤波；

自抑制模块：配置为对引导滤波后的双时相一维特征图像进行自抑制；

二值化模块：配置为基于自适应阈值方法对自抑制后的双时相一维特征图像进行二值化处理；

差分模块：配置为对二值化处理后的双时相一维特征图像进行差分，得到在获取双时相遥感光谱图像的两个时间点之间，水域发生变化的区域。

Images