卫星遥感影像多尺度融合变化检测方法、介质及电子装置与流程

1.本发明涉及遥感图像处理领域，涉及一种卫星遥感影像多尺度融合变化检测方法、介质及电子装置。


背景技术：

2.遥感变化检测目前分为有监督和无监督两种方式。
3.有监督变化检测需要较多的样本，其检测对象局限于某一类地物，例如建筑物的变化检测。基于ai的变化检测经过训练后表现较好，且能生成规则的范围边界，但是对输入的影像较为严格。如果模型训练时采用的rgb波段，卫星影像也只能舍弃近红外等其他波段，训练时若只针对建筑物，则只能对建筑物进行检测，无法探测其他变化。有监督变化检测还存在分类后检测的方法，分为面向对象和面向像元的，该方法极度依赖分类准确度。
4.无监督变化检测是指直接比较两幅影像的差异，将影像分为变化与非变化两类。现有方法可分为面向像元和面向对象的分类，基于像元的方法目前已较为成熟，且应用广泛。然而，随着遥感影像分辨率提升，基于像元的方法在处理高分辨率影像时存在噪声过多的问题，检测结果很难转换为理想的矢量结果。面向对象的方法很好的解决了这一问题。目前面向对象变化检测比较流行的方式采用图像分割的方法得到分割对象，该对象可包含特定地物的光谱信息、纹理信息、形状信息等。利用这些信息计算得到特征图，再采用传统的基于像元的方法进行分类。
5.面对多源遥感数据时，为了充分利用影像的光谱信息和纹理信息，可以采用面对对象的无监督变化检测方法。在实际应用中发现，目前的面向对象无监督变化检测存在受阴影干扰、纹理计算复杂、缺少大图像处理策略等问题。


技术实现要素：

6.为了解决上述任一技术问题，本技术实施例提供了一种卫星遥感影像多尺度融合变化检测方法、介质及电子装置。
7.为了达到本技术实施例目的，本技术实施例提供了一种卫星遥感影像多尺度融合变化检测方法，包括：将变化前的遥感影像和变化后的遥感影像均划分为多个图像块；按照每个图像块中地物的大小，确定每个地物的分割尺度；按照每个地物的分割尺度，从每个图像块提取地物，得到每个图像块中的对象；生成在变化前的遥感影像和变化后的遥感影像中每个对象的变化强度图；对不同分割尺度对应的变化强度图进行融合；确定融合后的变化强度图的变化检测结果。
8.一种存储介质，所述存储介质中存储有计算机程序，其中，所述计算机程序被设置为运行时执行上文所述的方法。
9.一种电子装置，包括存储器和处理器，所述存储器中存储有计算机程序，所述处理
器被设置为运行所述计算机程序以执行上文所述的方法。
10.本技术的有益效果为：根据地物的大小，确定图像块中地物的分割尺度，并采用多尺度分割进行对象的提取，并按照分割尺度，获取每个分割尺度的变化强度图，再将全部分割尺度的变化强度图进行融合，得到变化前的遥感影像和变化后的遥感影像对应的变化强度图，实现对变化信息的识别，提高图像变化信息的识别准确度。此外，通过基于直方图的阴影检测降低了阴影影像元的干扰，以对象为单位生成灰度共生矩阵能有效降低计算量，本技术所用方法能充分利用计算机的多核性能，并行处理大区域的高分影像。
11.将最大类间方差阈值法与神经网络模型中sigmoid激活函数结合,使得特征在归一化的同时，进行预分类，减弱了极大异常点的影像和极小异常点的影像，使得后续变化强度图可分性更强。
12.利用具有一定重叠度的滑窗策略使得算法适配大区域遥感影像，在分块出的重叠部分使得由于分块造成的分割差异被降到最小，使得分块策略可发挥计算机的并行性能，算法的运行时间更短。
13.本技术实施例的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本技术实施例而了解。本技术实施例的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。
附图说明
14.附图用来提供对本技术实施例技术方案的进一步理解，并且构成说明书的一部分，与本技术实施例的实施例一起用于解释本技术实施例的技术方案，并不构成对本技术实施例技术方案的限制。
15.图1为本技术实施例提供的卫星遥感影像多尺度融合变化检测方法的流程图；图2为本技术实施例提供的图像分块的示意图；图3为阴影像元在直方图中的分布示意图；图4为本技术实施例提供的曲线拟合结果的示意图；图5为本技术实施例提供的阴影检测结果；图6为变化前影像和变化后影像的图像和人工绘制标签图；图7为本发明提出方法的标签图和置信图。
具体实施方式
16.为使本技术实施例的目的、技术方案和优点更加清楚明白，下文中将结合附图对本技术实施例的实施例进行详细说明。需要说明的是，在不冲突的情况下，本技术实施例中的实施例及实施例中的特征可以相互任意组合。
17.图1为本技术实施例提供的卫星遥感影像多尺度融合变化检测方法的流程图。如图1所示，该方法包括：步骤101、将变化前的遥感影像和变化后的遥感影像均划分为多个图像块；其中，变化前的遥感影像和变化后的遥感影像的图片尺寸和显示比例相同，可以按照相同的分割策略将上述两个遥感影像分割为多个图像块。
18.步骤102、按照每个图像块中地物的大小，确定每个地物的分割尺度；地物指的是地表面的固定性物体，包括自然形成和人工建造的。在遥感图像中，地物可以指地球表面上相对固定的物体，例如，山脉、树木、楼宇等。
19.按照图像块中地物的大小，可以将地物划分为多个等级，每个等级的地物有各自对应的分割尺寸，其中，对于小地物，利用较细的分割尺度能从图像块中得到体现该地物轮廓和细节的分割对象；对于大地物，则利用较粗的分割尺度从图像块中即可得到体现该地图轮廓和特征的分割对象。
20.由于遥感影像上地物复杂多样，按照图像块中地物的大小，确定地物的分割尺度，在实现对地物分割的前提下，使得不同地物能够得到用于体现地物图像信息的最佳分割尺度，为从图像块中分割地物提供处理依据。
21.步骤103、按照每个地物的分割尺度，从每个图像块提取地物，得到每个图像块中的对象；具体的，识别每个图像块中的地物，按照地物的分割尺度，从图像块中提取该地物的图像信息，得到图像块中的对象。通过采用由粗到细的多尺度分割，得到多分割尺度的对象，实现对象的精细化分类，为后续图像处理提供支持。
22.步骤104、生成在变化前的遥感影像和变化后的遥感影像中每个对象的变化强度图；具体的，利用同一对象在变化前的遥感影像中图像信息以及在变化后的遥感影像中的图像信息，可以得到该对象的变化强度图。
23.步骤105、根据每个对象的变化强度图，生成每个分割尺度对应的变化强度图；由于同一分割尺度的对象在大小相似度高，因此，可以利用同一分割尺度的对象生成该分割尺度对应的变化强度图，达到利用该分割尺度对应的变化强度图记录同一分割尺度的多个对象的图像变化信息的目的。
24.步骤106、对不同分割尺度对应的变化强度图进行融合；由于每个分割尺寸的变化强度图均记录有该分割尺度对象的变化信息，因此，通过将全部分割尺度对应的变化强度图进行融合，可以得到变化前的遥感影像和变化后的遥感影像中全部对象的变化信息。
25.步骤107、确定融合后的变化强度图的变化检测结果。
26.由于该融合后的变化强度图记录了全部对象的变化信息，可以根据该融合变化强度图得到发生变化的图像信息。
27.本技术实施例提供的方法，根据地物的大小，确定图像块中地物的分割尺度，并采用多尺度分割进行对象的提取，并按照分割尺度，获取每个分割尺度的变化强度图，再将全部分割尺度的变化强度图进行融合，得到变化前的遥感影像和变化后的遥感影像对应的变化强度图，实现对变化信息的识别，提高图像变化信息的识别准确度。
28.在实现本技术过程中，发现相关技术存在如下问题，包括：计算对象的光谱属性值通常采用均值法，对象内包含阴影像元会使得对象的光谱偏离。
29.针对纹理计算复杂且对分割尺度敏感，通过基于图像生成灰度共生矩阵，再采用滑窗计算。由于滑窗计算是以像元为步长，当分割尺度发生变化时，对象的纹理值将会发生
变化。除此之外，滑窗计算方式还存在效率较低和耗时长的问题。
30.另外，由于遥感影像的数据量较大，现有技术中的图像处理手段不适用于遥感影像的处理。
31.针对上述技术问题，本技术提出对应的解决方案，包括：针对光谱受阴影干扰的问题，本技术实施例提出基于直方图的阴影检测，将阴影像元检测出来后，再计算对象的光谱属性值，从而达到降低阴影像元干扰的面对。
32.针对纹理计算复杂且对分割尺度敏感的问题，本技术实施例提出以对象为单位生成灰度共生矩阵，针对对象计算纹理，该方法对分割尺度不敏感，且分割后对象数远比像元的数量少，从而有效降低计算量。
33.针对缺少大图像处理策略的问题，本技术实施例提出充分利用计算机的多核性能，并行处理大区域的高分影像。
34.下面对本技术实施例提供的方法进行说明：本技术实施例提供的方法，将变化前后的影像划分为多个图像块，对每个图像块依据多个分割尺度将遥感影像分割成对象，依次计算每个对象的光谱均值和纹理值，计算每个分割尺度的变化强度图，通过融合多尺度的变化强度图得到最终的变化检测结果。
35.以图像分割采用简单的线性迭代聚类（simple linear iterative clustering，slic）为例，通过对图像块做高斯平滑，并依据光谱和空间距离做k均值聚类，将k均值的期望类别数作为分割尺度的系数。
36.进一步的，通过对影像进行分开处理，充分利用计算机的多核性能，并行处理大区域的高分影像。
37.在对变化前后的影像进行分块操作时，相邻的图像分块中在横向和/或竖向具有重叠区域。
38.图2为本技术实施例系统的图像分块的示意图。如图2所示，在图像分块时，相邻分块具有重叠区域，其中圆点纹理为重叠区。
39.采用具有横向和竖向具有重叠区域的分块方式，这样能减弱边缘地区由于分块导致同一地物被分割为两个的情况。
40.在针对每个图像块进行处理时，由于该重叠区的图像信息会被多次处理，因此采用计数器对每个重叠区的处理次数进行计数，其中重叠区的最终结果为同一项特征值的累加和除去计数器的计数得到的均值作为该项特征值的大小。
41.进一步的，可以获取同一对象的至少两种类型的特征图来得到每个对象对应的变化强度图，但是由于类型的特征图的量纲不统一，通过对每种类型的特征图的归一化结果进行加权计算，以提升变化强度图的准确性。
42.具体的，可以执行如下步骤，包括：步骤a1、获取每种类型的特征图的最大类间方差的阈值；其中，表示特征图的最大类间方差的阈值；步骤a2、计算每种类型的特征图的最大类间方差与最大类间方差的阈值之间的差值，得到每种类型的特征图对应的差值；其中，表示特征图的最大类间方差；步骤a3、将每种类型的特征图差值与该至少两种类型特征图的标准差进行除法运
算，得到每种类型的特征图的商；其中， 表示特征图的标准差；步骤a4、采用归一化函数对每种类型的特征图的商进行处理，得到每种类型的特征图的归一化结果；具体的，可以sigmoid函数的归一化方式来统一不同特征图的量纲；具体计算表达式如下：其中，。
43.经过归一化处理后，不同类型的特征图的量纲实现统一，使得最大类间方差的阈值附近的差异被放大，极端值被抑制。
44.步骤a5、对每种类型的特征图的归一化结果进行加权计算，得到每个对象对应的变化强度图。
45.采用加权的方法融合光谱特征和纹理特征得到各个尺度的变化强度图。根据经验取值，光谱特征的权重为0.7，纹理特征的权重为0.3时最佳。
46.本技术实施例基于最大类间方差otsu和sigmoid函数的归一化方式来统一不同特征图的量纲。
47.其中该至少两种类型包括光谱特征图和纹理特征图，其中光谱特征图可以为每个波段的光谱特征图，纹理特征图可以为不同纹理指数的特征图。
48.可选的，如果一种类型的特征图包括至少两张特征图，则对该至少两张特征图中每个像元的特征值进行比较，确定数值最大的像元，利用数值最大的像元生成最终的特征图。
49.例如，对某一具有4个波段的影像，计算得到5个光谱特征图，对5个特征图中每个像元作比较，保留最大的值作为最终的光谱特征值。
50.通过如下方式得到光谱特征图，包括：分别计算同一对象在变化前的遥感影像和变化后的遥感影像的光谱属性；计算同一对象在变化前的遥感影像和变化后的遥感影像的光谱属性之间的差值的绝对值作为该对象的光谱特征图；具体的，计算每个分割对象的光谱均值作为该对象的光谱属性，在计算光谱均值时，利用阴影检测结果去除阴影像元对地物对象的干扰。
51.进一步的，若阴影像元数小于对象像元数的一半，则认为该对象内受到阴影干扰，计算该对象的光谱均值时将不考虑阴影像元；若阴影像元数大于对象像元数的一半，则该对象被判定为纯阴影对象，则将计算得到的光谱均值乘以预设的衰减系数后作为对象的光谱属性，减弱其对变化强度图的贡献。
52.通过如下方式得到纹理特征图，包括：分别计算同一对象在变化前的遥感影像和变化后的遥感影像的纹理指数；计算同一对象在变化前的遥感影像和变化后的遥感影像的纹理指数之间的差值的绝对值作为该
对象的纹理特征图；其中，每个对象的纹理指数是通过如下方式得到的，包括：提取每个对象的外接矩形，利用该外接矩阵计算灰度共生矩阵，再利用灰度共生矩阵计算纹理指数。
53.具体的，设置每个对象的外接矩形，其中该外接矩形中不属于对象的像元使用对象像元光谱均值代替。用这个外接矩阵计算不同角度(如，0
°
、 90
°
、 180
°
以及 270
°
)的灰度共生矩阵，再计算对应的纹理指数。
54.其中上述纹理指数可以为非相似性（dissimilarity）和能量（energy）。
55.进一步的，计算两幅影像特征图差值的绝对值作为最终特征图，再加入光谱的差值的l2范数作为补充。
56.对于多光谱影像，可根据应用需求加入植被指数和基于形态学的建筑物指数mbi等补充特征，其中植被指数归为光谱特征，形态学指数归为纹理特征。
57.进一步的，在执行光谱属性值的计算时，对象中的阴影像元是根据预设的阴影像元的阈值进行阴影检测得到的；通过如下方式得到阴影像元的阈值，包括：步骤b1、提取遥感影像的红绿蓝波段，将红绿蓝波段的反射率值转换到色彩表示法lab颜色空间中，得到亮度通道的直方图；由于亮度通道对阴影非常敏感，可统计整幅影像亮度通道的直方图。
58.步骤b2、确定该亮度通道的直方图中阴影像元的阈值；相关技术可以采用计算直方图的势函数，具体参见计算如下表达式：其中，，；在计算表达式（1）中， 表示在=时，的个数值求和，其中为大于或等于1的整数；其中，表示亮度值；表示势函数的结果值；表示中的最大值；表示直方图的取值；表示调节系数，选择合适的值使得势函数和直方图的各峰特性接近。
59.图3为阴影像元在直方图中的分布示意图。如图3所示，阴影像元的阈值在亮度通道的直方图的第一个波谷结束的地方。
60.根据势函数求得第一个波谷的位置即为阈值。但这种方法在高分影像中存在漏检的情况，因此，本技术实施例提供的方法为求得这一位置，将直方图的像元值和频率值作为x、y值带入十次多项式中，具体参见如下计算表达式：
其中，至均为实数。
61.利用最小二乘法求解。求解后得到多项式的表达式。对该表达式求一阶导和二阶导。利用一阶导为0处为极值和二阶导为0处为拐点的性质找到在该十次多项式曲线第一个极小值后的拐点位置，该位置即为阴影阈值。求解的结果如图4所示，曲线拟合结果及拐点位置，其中，黑点为拐点位置，该点的x值为阈值。
62.步骤b3、采用该阴影像元的阈值对遥感影像进行检测；如果像元的像元值小于该阈值，则确定该像元为阴影像元；否则，确定该像元不是阴影像元。
63.图5为本技术实施例提供的阴影检测结果。如图5所示，左侧图为实际遥感影像图，右侧图为阴影检测结果的示意图。其中，图像中阴影像元主要由于地形、树木和建筑等造成。
64.进一步的，在对不同分割尺度对应的变化强度图进行融合时，需要确定确每个分割尺度的变化强度图的权重，再利用每个分割尺度的权重，对不同分割尺度对应的变化强度图进行融合；其中，采用如下表达式计算第i个分割尺度的特征图的权重，包括：其中，为第i个分割尺度的特征图，std()为该特征图的标准差， 为第j个分割尺度的特征图。
65.用上述方法融合变化强度图后，最终再次采用最大类间方差法otsu计算变化强度图的二分阈值（0.4-0.6），大于或等于阈值的判断为变化，小于阈值的判断为非变化。最终输出结果为一张标签图和置信图，标签图指示了变化区域和非变化区域，置信图描述了每个像元属于变化类别的置信度，置信图即为变化强度图。图6为变化前影像和变化后影像的图像和人工绘制标签图。图7为本发明提出方法的标签图和置信图。
66.与现有技术相比，本技术实施例提供的解决方案具有如下优势，包括：1. 将阴影检测纳入到变化检测中，排除了不同传感器导致的阴影差异，使结果更加准确。
67.2. 阴影检测为现有直方图检测的改进方案，使得阴影阈值描述更加准确，阴影漏检率大幅降低。
68.3. 将最大类间方差阈值法与神经网络模型中sigmoid激活函数结合,使得特征在归一化的同时，进行预分类，减弱了极大异常点的影像和极小异常点的影像，使得后续变化强度图可分性更强。
69.4.利用具有一定重叠度的滑窗策略使得算法适配大区域遥感影像，在分块出的重叠部分使得由于分块造成的分割差异被降到最小，使得分块策略可发挥计算机的并行性
能，算法的运行时间更短。
70.本技术实施例提供一种存储介质，该存储介质中存储有计算机程序，其中，该计算机程序被设置为运行时执行上文任一项中的方法。
71.本技术实施例提供一种电子装置，包括存储器和处理器，该存储器中存储有计算机程序，该处理器被设置为运行该计算机程序以执行上文任一项中的方法。
72.本领域普通技术人员可以理解，上文中所公开方法中的全部或某些步骤、系统、装置中的功能模块/单元可以被实施为软件、固件、硬件及其适当的组合。在硬件实施方式中，在以上描述中提及的功能模块/单元之间的划分不一定对应于物理组件的划分；例如，一个物理组件可以具有多个功能，或者一个功能或步骤可以由若干物理组件合作执行。某些组件或所有组件可以被实施为由处理器，如数字信号处理器或微处理器执行的软件，或者被实施为硬件，或者被实施为集成电路，如专用集成电路。这样的软件可以分布在计算机可读介质上，计算机可读介质可以包括计算机存储介质（或非暂时性介质）和通信介质（或暂时性介质）。如本领域普通技术人员公知的，术语计算机存储介质包括在用于存储信息（诸如计算机可读指令、数据结构、程序模块或其他数据）的任何方法或技术中实施的易失性和非易失性、可移除和不可移除介质。计算机存储介质包括但不限于 ram、rom、eeprom、闪存或其他存储器技术、cd-rom、数字多功能盘（dvd）或其他光盘存储、磁盒、磁带、磁盘存储或其他磁存储装置、或者可以用于存储期望的信息并且可以被计算机访问的任何其他的介质。此外，本领域普通技术人员公知的是，通信介质通常包含计算机可读指令、数据结构、程序模块或者诸如载波或其他传输机制之类的调制数据信号中的其他数据，并且可包括任何信息递送介质。