一种基于对象的高分辨率遥感影像变化检测方法及系统

resolution triplet network，hrtnet)框架，作用于高分辨率遥感图像。其使用一个动态初始模块来增强表示多尺度特征的能力，使整体模型对不同大小的变化区域更加敏感。利用差分图像中包含的时间特征，确保了模型对表面变化的关注，增强了其对伪变化识别的鲁棒性。该框架的可以识别感兴趣的变化信息，过滤掉不相关的变化信息的干扰。
7.面向对象的变化检测(object-oriented change detection,oocd)算法主要分为三类：基于类对象的、基于图像对象的和基于多时相对象的方法。基于类对象的方法高度依赖于图像分类的准确性，而图像分类的高精度则依赖于训练样本的质量和数量。通过设置阈值，基于图像对象的方法可以直接比较不同的时间图像。然而，多时相图像的独立分割结果呈现出不同的大小和形状，这使得对象难以被找到。
8.现有技术10利用合并的多时相图像分割结果构建图像目标，计算归一化植被指数(normalized difference vegetation index，ndvi)、平均值和标准差，并将其作为变化指标。其将统计异常值挖掘方法与变化检测相结合。方法假设样本数据服从高斯分布，采用卡方迭代裁剪方法寻找变化对象。但是，该方法主要用于森林变化检测，而且仅使用ndvi作为变化指数，总体准确率为90％，总体kappa系数为0.8。现有技术11对技术10的方法进行了改进，引入马氏距离作为相似度量方法，放弃使用特征差作为测量反射率变化程度的指标。
9.现有技术12采用卡方和迭代修剪的变化检测方法来寻找被困在样本数据中的变化对象，该方法并不严格服从高斯分布，导致检测效果不理想。现有技术13为了修正马氏距离的不稳定性和对小变量的敏感性等缺点，提出了一种基于余弦定律和箱线图的变化检测方法。现有技术14使用对象作为最小分类单位，克服了基于像素的分类方法导致的椒盐效应问题。现有技术15提出了一个使用基于图像对象的方法来寻找给定特征类型的最优分割的框架。现有技术16公开了一种基于目标/邻域相关性的图像分析和图像分割技术的变化检测方法。现有技术17基于对均匀区域计算特征的比较，提出了一些无监督的方法。
10.通过上述分析，现有技术存在的问题及缺陷为：现有技术采用单一尺度进行检测，不确定性大，与真实变化误差大，检测结果存在不确定性，不能应用于高分辨率遥感图像检测中；同时现有技术受噪声影响大，无法有效去除遥感影像的高斯、脉冲噪声。


技术实现要素：

11.针对现有技术存在的问题，本发明提供了一种基于对象的高分辨率遥感影像变化检测方法及系统。
12.本发明是这样实现的，一种基于对象的高分辨率遥感影像变化检测方法，所述基于对象的高分辨率遥感影像变化检测方法包括：
13.首先，对高分辨率遥感数据进行辐射校正、几何校正；对校正后的高分辨率遥感数据依次进行对等组滤波、颜色量化、j值计算、空间分割处理；
14.其次，使用ssim与psnr两种相似性测度分别对对象的相似性进行定量评价；利用多尺度融合策略进行变化检测，使用景观格局指数中的分维数过滤检测结果中的疑似误检区域得到变化检测结果。
15.进一步，所述基于对象的高分辨率遥感影像变化检测方法包括以下步骤：
16.步骤一，对高分辨率遥感数据进行辐射校正、几何校正；通过对等组滤波对两幅经过校正的遥感图像进行平滑与去噪；
17.步骤二，通过使用修改的快速k均值聚类，合并颜色相近的聚类得到颜色量化图像；
18.步骤三，在不同的尺度上对量化图像进行jseg图像分割得到多尺度地物对象；使用ssim与psnr两种相似性测度分别对对象的相似性进行定量评价；
19.步骤四，利用多尺度融合策略进行变化检测；使用景观格局指数中的分维数过滤疑似误检区域得到变化检测结果。
20.进一步，所述通过对等组滤波对两幅经过校正的遥感图像进行平滑与去噪包括：
21.(1)计算di(n)的一阶差分fi(n)：
22.fi(n)＝d
i+1
(n)-di(n)；
23.其中，di(n)表示窗口内所有像素到中心像素x0(n)的欧式距离；
24.(2)利用下式对xi(n)的前后m个点进行测试，判断各个点否属于脉冲噪声：
25.fi(n)≤α；
26.其中，m＝w/2，窗口大小的一半；α对于高度损坏的图像设置为较大值，对于轻微损坏的图像设置为较小值；
27.(3)若fi(n)不满足α，则令j≤i；或j》i的终点xj(n)视为脉冲噪声并去除；剩下的dj(n)用于估计真正地对等组；
28.(4)将像素x0(n)替换为对等组成员的加权平均值：
[0029][0030]
其中，wi表示标准高斯滤波权值；
[0031]
(5)若x0(n)属于脉冲噪声，并且没有对等组，则所述位置的真实对等组通过窗口中的其他对等像素进行估计。
[0032]
进一步，所述通过使用修改的快速k均值聚类，合并颜色相近的聚类得到颜色量化图像包括：
[0033]
1)获取每个对等组最大距离t(n)：
[0034]
t(n)＝d
m(n)-1
(n)；
[0035]
其中，di(i＝0，m(n)-1)表示对等组成员像素到中心像素的距离；m(n)表示对等组的大小，即对等组成员数量；
[0036]
2)得到对等组的最大距离t(n)的值之后，以像素x(n)为中心表示局部区域的平滑度，利用下式计算每个像素x(n)的权重v(n)：
[0037]
v(n)＝exp(-t(n))；
[0038]
3)利用下式确定颜色量化中聚类的初始数量k：
[0039]
k＝βt
avg
；
[0040]
其中，t
avg
表示t(n)的平均值，β表示参数；
[0041]
4)设置包含像素权重质心的聚类中心的更新规则；通过下式计算一个颜色类ci的质心ci：
[0042]
[0043]
5)使用结合权重的k均值确定初始聚类的中心，利用下式进行加权失真测量值，确定待拆分集群，直到达到初始的集群数量k：
[0044]di
＝∑v(n)||x(n)-ci||2,x(n)∈ci；
[0045]
6)对相近的聚类进行合并，令两个质心之间的最小距离满足预设的阈值，分配每一个像素到最近的聚类质心，得到颜色量化图像。
[0046]
进一步，所述在不同的尺度上对量化图像进行jseg图像分割得到多尺度地物对象包括：
[0047]
(1)确定颜色量化图像的所有n个点的集合得到点集z，将z按颜色类分成c部分，z＝{zi，1≤i≤c}，计算点集zi的中心mi：
[0048][0049]
其中，z＝(x，y)，z∈z，是点在图像中的坐标；
[0050]
(2)确定区域的j值：
[0051]
j＝(s
t-sw)/sw；
[0052]
其中：
[0053][0054][0055]
(3)把颜色量化图像的全域u先划分为k个区域ui(1≤i≤k)，每个区域内有mi个点，分别计算区域ui的j值ji，且计算区域u的平均j值j
avg
：
[0056][0057]
(4)构造一个图像，将所述图像像素值对应于以像素为中心的小窗口上计算的j值，得到j图像，并将相应的像素值作为局部j值；
[0058]
(5)计算j图像的平均值和标准差，分别用μj和σj表示；并计算阈值变量tj：
[0059]
tj＝μj+aσj；
[0060]
其中，a表示几个不同的预设值；
[0061]
(6)筛选初始种子：将局部j值小于tj的像素作为候选种子点：基于4连通性连接候选种子点，得到候选种子区域；利用不同的a计算不同的阈值，得到不同的种子区及数量，取数量最多的作为最后结果；若候选种子区域的尺寸大于相应最小尺寸，则确定为种子；
[0062]
(7)填充种子中的空洞；平均区域剩余的未分配部分的局部j值，并连接低于平均值的像素，形成增长区域；如果一个生长区域毗邻一颗且只有一颗种子，被分配给种子；
[0063]
(8)计算下一个较小尺度下剩余像素的局部j值，重复步骤(7)，定位边界；以最小的尺度逐个增长剩余的像素。将j值最小的像素分配给相邻的种子，并更新用于存放种子边界处的未分类像素的缓冲区，直到对所有像素进行分类。
[0064]
进一步，所述利用多尺度融合策略进行变化检测；使用景观格局指数中的分维数过滤检测结果中疑似误检区域得到变化检测结果包括：
[0065]
1)利用下式计算每个对象rj的多尺度相似性sj：
[0066][0067]
2)判断每个对象rj的多尺度相似性sj与变化检测的相似性测定阈值ts的大小关系，若sj≥ts，则rj保持不变，标记为0；否则，判断rj发生了变化，标记为1；
[0068]
3)重复步骤1)至步骤2)，直至完成分割结果中的所有对象；
[0069]
4)对所有检测得到的变化区域对象利用下式计算分维数：
[0070][0071]
其中，c表示对象的周长；s表示对象的面积；
[0072]
5)对比计算得到的变化区域对象的分维数index是否大于t
idx
，若大于则判断非建筑，并在变化检测结果中排除，遍历所有变化区域对象得到变化检测结果。
[0073]
本发明的另一目的在于提供一种实施所述基于对象的高分辨率遥感影像变化检测方法的基于对象的高分辨率遥感影像变化检测系统，所述基于对象的高分辨率遥感影像变化检测系统包括：
[0074]
数据处理模块，用于对高分辨率遥感数据进行辐射校正、几何校正；
[0075]
对象提取模块，用于对校正后的高分辨率遥感数据依次进行对等组滤波、颜色量化、j值计算、空间分割处理，得到多个对象；
[0076]
对象分析模块，用于使用ssim与psnr两种相似性测度分别对对象的相似性进行定量评价，得到检测结果；
[0077]
多尺度融合和分维数过滤模块，用于利用多尺度融合策略进行变化检测，使用景观格局指数中的分维数过滤检测结果中的疑似误检区域得到变化检测结果。
[0078]
本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述基于对象的高分辨率遥感影像变化检测方法的步骤。
[0079]
本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行所述基于对象的高分辨率遥感影像变化检测方法的步骤。
[0080]
本发明的另一目的在于提供一种信息数据处理终端，所述信息数据处理终端用于实现所述基于对象的高分辨率遥感影像变化检测系统。
[0081]
结合上述的技术方案和解决的技术问题，本发明所要保护的技术方案所具备的优点及积极效果为：
[0082]
本发明使用对等组滤波方法对遥感影像进行平滑与去噪，能有效地去除高斯、脉冲噪声并能在较好地保留边缘信息的情况下进行平滑处理。
[0083]
本发明提供了一种改进的k均值聚类取代普通的非监督聚类方法。使用启发式方法找到聚类的初始质心，改进了k均值算法的运行时间和最终解的质量。通过在初始质心的选取与聚类质心迭代更新的过程增加权重以在纹理区域生成更多的聚类，再对质心相近的
聚类加以合并，以更精确得到检测边缘。
[0084]
本发明提供的对边缘感知比较敏感、对像素分布均匀的平滑区域感知较弱的分割算法，有效地提取对象区域与边界，实现了图像的分割。在多个不同的尺度下对遥感图像进行分割，得到了不同尺度的对象。
[0085]
本发明通过计算两种评价图像相似性的相似性测度ssim与psnr，在每一个尺度下计算每一个对象相似度以定量地对对象间的相似程度进行分析。
[0086]
本发明通过构建多尺度融合策略以有效地结合不同尺度下图像对象的特征，对图像进行变化检测。通过使用景观格局指数的分维数作为指标衡量变化区域形状的规整程度，过滤检测结果中的疑似误检的区域。
[0087]
本发明在数据集whu building dataset上进行实验取得较好结果，总体精度达到96.78％，kappa系数达到0.8402。证明本发明方法的有效性。
[0088]
本发明在高分辨率遥感图像的城市变化检测中是有效和可靠的。使用jseg算法不仅可以实现对场景中目标的准确提取，还可以利用j图像序列中包含的多特征进行变化检测，进一步应用融合策略可以得到最终结果。实验证明，该方法克服了单尺度检测的不确定性，产生了更接近真实变化的检测结果。
[0089]
本发明即使在单一尺度上的性能也有着不错的效果，从而证明了面向对象的变化检测算法研究可以满足对不断发展的高分辨率遥感图像的需求。
[0090]
本发明中的两种相似性测度都有各自的优劣。ssim测度虽在本实验所用的数据集中没能取得较好的效果，但它的有界特性使得它能在其他环境下取得比psnr更优的结果。psnr测度在本实验中有较高的精度，较低的误报率与漏检率；相反它的无界的特性会给阈值的设定造成困难。
[0091]
本发明的景观格局指数过滤操作能提高检测结果的精度，使变化区域更加清晰明了。
[0092]
本发明的变化检测对区域的划分可以为实地工作提供有价值的参考信息，从而减少工作量，节约资源。
[0093]
本发明不易受噪声影响，应用图像纹理特征进行变化检测的oocd变化检测方法，构建了一个变化检测框架，并在数据集whu building dataset上取得较好结果，且能够快速准确地检测出具体的变化区域。
附图说明
[0094]
图1是本发明实施例提供的基于对象的高分辨率遥感影像变化检测方法原理图；
[0095]
图2是本发明实施例提供的基于对象的高分辨率遥感影像变化检测方法流程图；
[0096]
图3是本发明实施例提供的j＝1.720的class-map图；
[0097]
图4是本发明实施例提供的j＝0的class-map图；
[0098]
图5是本发明实施例提供的j＝0.855的class-map图；
[0099]
图6是本发明实施例提供的j
+
＝0，j
*
＝0，jo＝0，j
avg
＝0最优分割图示；
[0100]
图7是本发明实施例提供的j
+
＝0，j
{*,o}
＝0.011，j
avg
＝0.05最优分割图；
[0101]
图8是本发明实施例提供的计算局部j值的基本窗口示意图；图中只有“+”点用于计算局部j值；
[0102]
图9是本发明实施例提供的2级窗口的降采样示意图；图中只有“+”点用于计算局部j值；
[0103]
图10是本发明实施例提供的2012年某地区遥感图；
[0104]
图11是本发明实施例提供的2016年某地区遥感图；
[0105]
图12是本发明实施例提供的2012年某地区多时相图像的尺度3的j图像；
[0106]
图13是本发明实施例提供的2016年某地区多时相图像的尺度3的j图像；
[0107]
图14是本发明实施例提供的2012年某地区多时相图像的尺度3下的分割结果示意图；
[0108]
图15是本发明实施例提供的2012年某地区多时相图像的尺度3下的分割结果示意图。
具体实施方式
[0109]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0110]
如图1－图2所示，本发明实施例提供的基于对象的高分辨率遥感影像变化检测方法包括以下步骤：
[0111]
s101，对高分辨率遥感数据进行辐射校正、几何校正；通过对等组滤波对两幅经过校正的遥感图像进行平滑与去噪；
[0112]
s102，通过使用修改的快速k均值聚类，合并颜色相近的聚类得到颜色量化图像；
[0113]
s103，在不同的尺度上对量化图像进行jseg图像分割得到多尺度地物对象；使用ssim与psnr两种相似性测度分别对对象的相似性进行定量评价；
[0114]
s104，利用多尺度融合策略进行变化检测；使用景观格局指数中的分维数过滤检测结果中的疑似误检区域得到变化检测结果。
[0115]
下面对本发明的具体方案进行详细说明：
[0116]
1、变化检测
[0117]
1.1变化检测一般流程
[0118]
传统的直接比较法由于其简单、易于实现，在实践中得到了广泛的应用。其他的方法在一定程度上都是借鉴了直接比较法的部分思想发展而来。首先是对多时相遥感图像进行预处理，该步骤一般包含辐射校正、几何校正、图像配准等操作。然后是计算多时相图像之间的差分图像，差分图像可以是绝对差值、各种距离和植被指数差分等。最后一般是设置一个阈值对差分图像进行变化检测，阈值可以是根据研究人员的经验设置一个或多个固定值，也可以通过算法求取全局阈值或动态阈值。
[0119]
1.2基于ndvi植被指数差分的变化检测
[0120]
基于归一化植被指数(normalized difference vegetation index，ndvi)差分的卫星图像变化检测方法是一种简单有效的基于代数的变化检测，它已经被许多学者用于研究各地区的植被覆盖变化。例如gandhi等利用tm图像研究了印度vellore地区在不同的ndvi差分图像阈值下的变化检测与其他土地植被特征。
[0121]
基于植被指数dnvi差分的变化检测方法简单概括为以下几步：(1)通过已经配准
两个时相的多光谱反射率图像计算ndvi图像。(2)计算两幅ndvi图像的差分图像。(3)差分图像与设定的阈值进行比较，小于该阈值的则认为是不变区域，否则为变化区域。
[0122]
方法的第一步是计算两个时相的多光谱反射率图像，然后计算ndvi图像。ndvi已被广泛用于研究光谱变异性与植被生长速率变化之间的关系，它也有助于确定绿色植被的生长，以及检测植被的变化。考虑两幅在同一区域不同时相t1和t2的遥感图像x1与x2，大小都是h
×
w像素，每个时相的图像都有b个波段，设图像已经经过校正与配准等预处理。设x
i，b
为原图像xi的第b个波段中大小为h
×
w像素的单波段图像(i＝1，2；b＝1，b)。时相i的ndvi图像ndvii由下式计算。
[0123][0124]
其中，niri为时相i的近红外反射率图像，在tm图像中为x
i，4
；ri为时相i的可见光红色光反射率图像，在tm图像中为x
i，3
。
[0125]
方法的第二步是计算ndvi差分图像，计算公式如下。
[0126]
ndvid＝|ndvi
1-ndvi2|
ꢀꢀ
(1-2)
[0127]
方法的第三步是与阈值比较，从而实现变化检测。检测结果一般是二值图像，计算公式如下。
[0128]
c＝ndvid≥t
ndvi
ꢀꢀ
(1-3)
[0129]
其中t
ndvi
是变化检测的阈值，可以由研究人员通过经验设置，也可以通过其他自适应方法得到。c是变化检测的结果图像，通过逻辑表达式计算，值为1表示变化像素，值为0表示不变像素。
[0130]
基于植被指数差分的阈值变化检测简单而有效，适用于多地的林草地变化检测，也可以根据检测目标的不同更换不同的指数。方法阈值的选取依赖研究人员的经验与统计，自动化较差。
[0131]
1.3基于遥感图像无监督变化的动态阈值检测方法
[0132]
为了区分变化像素和不变像素，最常见的解决方案是使用阈值算法来选择全局阈值对差分图像进行二值化。尽管其简单且效果良好，但最佳阈值的选择通常应与有关场景或视觉解释的先验知识相关联以赋予意义。然而，大多数的阈值决策方法都没有考虑到先验信息。此外，由于差分图中两个簇的普遍重叠问题，全局阈值通常不适合所有像素。因此，有时一个全局阈值会产生很差的结果。
[0133]
现有技术将fcm对每个像素的隶属度值作为先验信息，指导每个像素的阈值决策。方法如下：首先，采用变化向量分析(change vector analysis，cva)技术生成差分图像。然后利用高斯混合模型(gaussian mixture model，gmm)对差图像的变化和无变化像素集进行建模，采用期望最大值(expectation maximum，em)算法估计差分图像的统计参数。因此，基于贝叶斯决策理论，可以确定一个全局的初始阈值tg。最后，通过利用fcm算法生成的每个像素u(y，x)的隶属度值来修改初始全局阈值，可以得到一个动态阈值td。
[0134]
方法的第一步是生成差分图像。考虑两幅在同一区域不同时相t1和t2的遥感图像x1与x2，大小都是h
×
w像素，每个时相的图像都有b个波段，设图像已经经过校正与配准等预处理。设x
i，b
为原图像xi的第b个波段中大小为h
×
w像素的单波段图像(i＝1，2；b＝1，b)。双时相多波段图像的二维差分图像xd可以由下式定义：
[0135][0136]
方法的第二步是初始全局阈值识别。差分图像的变化像素和不变像素，分别标注为cc、cn，首先通过em算法和贝叶斯决策理论(em-bayes)识别的全局阈值进行区分。在简单的假设下，差分图像中像素的统计分布可以用cc和cn的混合高斯模型进行建模，它可以表示为：
[0137][0138]
其中ni(x；μi，σ
i2
)为正态分布，期望为μi，方差为σ
i2
；αi是每个分布的混合权值；m是高斯分布的个数，这里m为2。
[0139]
在此基础上，利用em算法通过迭代过程估计gmm的参数。然后，根据贝叶斯最小误差规则，将差分图像中的每个像素分配给后验条件概率最大的类。最后，通过求解以下二次方程，可以得到初始的全局阈值。
[0140][0141]
式中(μn、σ
n2
，αn)和(μc，σ
c2
，αc)分别为变化和无变化像素集的期望、方差和混合权重。
[0142]
方法的第三步是结合fcm的隶属度计算动态阈值进行变化检测。fcm算法是一种迭代聚类方法，同一聚类中的像素具有较高的相似度，而不同聚类之间的相似度较低，并且每个像素的聚类通过隶属度矩阵的模糊方式进行识别。fcm试图通过迭代最小化目标函数jm来寻找给定数据集的模糊划分：
[0143][0144]
其中x＝{x1，x2，xn}是m维向量空间中的数据集，c(2≤c≤n)是聚类数，u
ji
是xi在第j个聚类中的隶属度，m是每个模糊隶属度的加权指数，cj是聚类j的中心，d2(xi，cj)是xi和中心cj之间的距离度量。
[0145]
得到fcm隶属度矩阵之后，每个像素的动态阈值修改为：
[0146][0147]
其中，tg为上述em算法和贝叶斯决策理论得到的最优初始阈值，un(y，x)和uc(y，x)分别为位置(y，x)处像素属于无变化像素和有变化像素的隶属度值。如果0《un(y，x)/uc(y，x)《1，则表示像素有更大的概率被识别为变化像素，阈值将被公式(1-8)降低。如果1《un(y，x)/uc(y，x)《e-1，像素的阈值也将降低，以解决fcm算法的高误检率。如果un(y，x)/uc(y，x)》e-1，这表示像素有更大的概率被识别为不变像素，阈值被升高。
[0148]
最后，使用动态阈值对差分图像xd进行变化检测：
[0149]
c＝xd≥tdꢀꢀ
(1-9)”表示向下取整。该集合的平均向量被定义为：
[0159][0160]
每个向量与平均向量求差：
[0161][0162]
主成分分析寻找n个标准正交向量及其最能描述数据分布的相关标量λs并应用于差异向量集δ
p
。向量es和相应的标量λs分别为协方差矩阵c的特征向量和特征值。
[0163][0164]
其中，上标符号“t”对应于向量的转置。矩阵c大小为h2×
h2，它决定了h2个特征向量和特征值。假设c生成的特征向量根据它们的特征值按降序排序，即λs≥λ
s+1
。
[0165]
特征向量空间是通过将xd(i，j)投影到空间位置(i，j)处的每个像素的特征向量空间上来创建的，即：
[0166]
v(i,j)＝[v1,v2,
…
,vs]
t
ꢀꢀ
(1-15)
[0167]
其中
[0168][0169]
其中1≤s≤h2，1≤s≤s。参数s决定了特征向量v(i，j)在空间位置(i，j)处得维数，并且只是用于将xd(i，j)投影到特征向量空间上的特征向量的数量。向量xd(i，j)与公式(1-11)得到的方法相同。
[0170]
该方法的下一阶段是利用k均值聚类算法对特征向量空间进行聚类，生成两个聚类。设v
wu
和v
wc
分别为wu和wc类的聚类均值特征向量。为了对k均值聚类算法生成的聚类进行标记，利用k均值聚类算法中的标记像素来寻找差异图像上的两个平均值。当特定区域的两幅图像之间发生变化时，预计该区域的差分像素值高于没有变化区域的像素值。利用此假设，将差分图像中像素平均值较低的聚类划分为wu类，将另一个聚类划分为wc类。
[0171]
通过v
wc
和v
wu
，可以创建一个二进制变化检测图cm＝{cm(i，j)|1≤i≤h，1≤j≤w}，其中“1”表示相应的像素位置发生变化(即属于wc类)，而“0”表示不发生任何变化(即属于wu类)。这个过程可以看作是无监督的阈值。
[0172][0173]
其中，符号“d()”表示欧氏距离。
[0174]
该算法计算简单，适用于实时检测应用，在对抗均值为零的高斯噪声和散斑噪声方面均具有较好的性能。
[0175]
以上详细地介绍了变化检测的一般流程与一些变化检测的相关理论方法。这些方法中包括基于ndvi植被指数差分的变化检测、基于遥感图像无监督变化的动态阈值检测和主成分分析结合k均值聚类的无监督变化检测。详细介绍了差分图像、植被指数差分图像与多波段图像的二维差分图像的生成。详细介绍了全局阈值的经验确定与基于em算法和贝叶斯决策理论的全局阈值确定。介绍了fcm基本概念与全局阈值与fcm结合的动态阈值确定方
法。详细介绍了非重叠图像块的生成、特征向量空间的创建以及如何利用k均值聚类进行变化检测。
[0176]
2、滤波与颜色量化
[0177]
2.1对等组滤波
[0178]
噪声去除和图像平滑是许多图像处理应用的重要内容。例如，在彩色图像量化、运动估计和图像分割中，经常使用高斯滤波和中值滤波作为预处理步骤。对于彩色图像，去除脉冲噪声的一种常见方法是通过向量中值滤波(vector median filtering，vmf)。其他方法包括矢量方向滤波(vector directional filtering，vdf)和方向距离滤波(directional-distance filtering，ddf)。后者是vmf方法和vdf方法的结合。这些方法的缺点之一是，它们通常倾向于修改未被噪声破坏的像素。在矢量方向滤波中，首先使用一个类teager的算子来检测异常值，以便只替换有噪声的像素。但检测过程是对每个单独的颜色成分进行的，这可能会导致最终结果的错误。针对混合高斯噪声和脉冲噪声的情况，在方向距离滤波中提出了一种自适应非线性多元滤波方法。但是，由于整个局部窗口的平均值被用来估计原始像素值，它可能会模糊边缘和细节。
[0179]
为了解决这些缺点，现有技术公开了一种非线性的彩色图像噪声去除算法对等组滤波(peer group filtering，pgf)，这个方法不会模糊边缘。
[0180]
设x0(n)表示一个图像像素向量，取以其为中心的w
×
w大小的窗口，计算窗口内所有像素到x0(n)的欧氏距离并按升序排序，并将其表示为xi(n)(i＝0、k＝w
2-1)。欧氏距离di(n)：
[0181]di
(n)＝||x0(n)-xi(n)||,i＝0,
…
,k
ꢀꢀ
(2-1)
[0182]
d0(n)≤d1(n)≤
…
≤dk(n)
ꢀꢀ
(2-2)
[0183]
中心像素x0(n)地对等组p(n)定义为：
[0184]
p(n)＝{xi(n),i＝0,
…
,m(n)-1}
ꢀꢀ
(2-3)
[0185]
其中m(n)为p(n)的大小，即成员像素数量。
[0186]
对等组p(n)由x0(n)本身及其具有相似颜色的相邻点组成。如何根据局部统计数据为每个对等组选择合适的大小m(n)对算法的成功具有重要意义。一种简单的方法可以是设置一个阈值t(n)，使m(n)满足：
[0187]dm(n)-1
(n)≤t(n)且d
m(n)
(n)≥t(n)
ꢀꢀ
(2-4)
[0188]
然而，由于信号和噪声统计可以在不同的图像甚至在同一图像中发生变化，所以很难找到一个固定的t(n)最优值。
[0189]
如果窗口中有两个颜色聚类，则可以用使类间散射与类内散射的比值最大化的fisher线性判别来分离这两个类。但是，对于两个以上的类，该方法将无法分离包含中心像素x0(n)的聚类。此外，在三维空间中，这个方法的计算复杂度也较高。
[0190]
要规避这些问题，一种简单的方法是只使用一维距离di(n)的fisher的判别估计。使之最大化的判别标准是：
[0191][0192]
其中：
[0193][0194][0195][0196][0197]
该算法为每个i计算j(i)，并找到j(i)为最大值的截止位置，即：
[0198]
m(n)＝argmax(j(i)),i＝0,
…
,k
ꢀꢀ
(2-10)
[0199]
为了消除脉冲噪声的影响，在进行对等组分类前，计算了距离di(n)的一阶差分fi(n)：
[0200]fi
(n)＝d
i+1
(n)-di(n)
ꢀꢀ
(2-11)
[0201]
对xi(n)的前后m个点进行以下测试，检查它们是否属于脉冲噪声：
[0202]fi
(n)≤α
ꢀꢀ
(2-12)
[0203]
其中，m＝w/2，窗口大小的一半，α对于高度损坏的图像设置为较大值，对于轻微损坏的图像设置为较小值，实验中设置α＝12。如果fi(n)不满足该阈值，则将j≤i或j》i的终点xj(n)视为脉冲噪声并去除。剩下的dj(n)用于估计真正地对等组。
[0204]
在脉冲噪声去除和对等组分类后，将像素x0(n)替换为其对等组成员的加权平均值：
[0205][0206]
其中，wi是标准高斯滤波权值，这取决于pi(n)相对于x0(n)的相对位置。如果确定x0(n)属于脉冲噪声，并且没有对等组，则该位置的真实对等组通过窗口中的其他对等像素进行估计。
[0207]
2.2颜色量化
[0208]
在上一个步骤已经使用对等组滤波对图像进行平滑和噪声去除，在对等组滤波的过程中，获取每个对等组最大距离t(n)：
[0209]
t(n)＝d
m(n)-1
(n)
ꢀꢀ
(2-14)
[0210]
其中di(i＝0，m(n)-1)是对等组成员像素到中心像素的距离，m(n)是对等组的大小，即对等组成员数量。得到了地对等组的最大距离t(n)的值以像素x(n)为中心的表示局部区域的平滑度，且每个像素x(n)的权重v(n)计算方法由下式给出：
[0211]
v(n)＝exp(-t(n))
ꢀꢀ
(2-15)
[0212]
噪声区域中的像素加权小于平滑区域中的像素。
[0213]
t(n)的平均值t
avg
能表示整个图像的平滑度，一般来说，t
avg
越高，图像就越不光滑，需要更多的聚类来量化图像中的颜色。通过下式确定颜色量化中聚类的初始数量k：
[0214]
k＝βt
avg
ꢀꢀ
(2-16)
[0215]
其中，在实验中，β被设为2。
[0216]
在颜色量化中使用加权修改的k均值算法。聚类中心的更新规则被修改为包含像素权重质心。通过下式计算一个颜色类ci的质心ci：
[0217][0218]
由于像素向量被加上了权重，故质心向权值较高的点移动。
[0219]
初始聚类的中心使用结合权重的k均值加加确定。加权失真测量值通过下式定义，用于确定要拆分哪些集群，直到达到初始的集群数量k。
[0220]di
＝∑v(n)||x(n)-ci||2,x(n)∈ciꢀꢀ
(2-18)
[0221]
因此，权重较小的点将被分配更少的聚类。所以在详细区域中的颜色簇的数量就被抑制了。
[0222]
因为初始聚类中心的选择方法中，权重较小的点将被分配更少的聚类，纹理区域中的聚类的数量就被抑制了。因此，在颜色量化之后，具有相同颜色的大量像素将有多个聚类，因为k均值聚类的目标是最小化全局失真。因此，需要对相近的聚类进行进一步合并，使两个质心之间的最小距离满足预设的阈值，实验中使用的阈值是12。分配每一个像素到其最近的聚类质心，得到最终的量化图像。
[0223]
2.3 jseg图像分割
[0224]
基于对象的变化检测目标是要将两个时相遥感图像中的地物对象准确地提取出来。jseg分割算法是最常用的彩色图像分割方法之一。jseg对边缘比较敏感，而对于像素分布均匀的对象内部感知较弱，因此可以很好地提取对象的边界，进而实现图像分割。所以将使用jseg主要部分算法提取两个时相遥感图像中的地物对象，以便于完成oocd变化检测。
[0225]
2.3.1 j值与j图像
[0226]
在之前的步骤中图像中的颜色被粗量化，而不会显著降低颜色质量，其目的是只提取几种具有代表性的颜色，可以区分图像中的邻近区域。通常，在自然场景的图像中只需要10到20种颜色。颜色被量化后，量化的颜色被分成多个类，在这里这种以标签作为像素值的图像被称为class-map类图。
[0227]
设点集z是一个class-map中所有n个点的集合，u是class-map的全域。z＝(x，y)，z∈z，是点在图像中的坐标。m是点集的坐标中心。
[0228][0229]
将z按颜色类分成c部分，即z＝{zi，1≤i≤c}，则mi为点集zi的中心。
[0230][0231]
定义区域的j值：
[0232]
j＝(s
t-sw)/swꢀꢀ
(2-21)
[0233]
其中：
[0234]
[0235][0236]
根据图3、图4、图5的class-map图及对应j值可以看到，j值表征了颜色类在区域内分布的均匀情况。图3几个类分别扎堆，分布并不均匀；图4每个类在区域u中都是均匀分布；图5是前面二者的综合情况。
[0237]
若把区域u先划分成k个区域ui(1≤i≤k)，每个区域内有mi个点，分别计算区域ui的j值ji，且定义区域u的平均j值j
avg
：
[0238][0239]
以j
avg
作为标准，提出所有可能的分割图像的方法。对于固定数量的区域，一个较好的分割往往有一个较低的j
avg
，所以对于图3与图5的最优分割分别如图6、图7所示。
[0240]
图6、图7最优分割图示。(a)j
+
＝0，j
*
＝0，jo＝0，j
avg
＝0；(b)j
+
＝0，j
{*,o}
＝0.011，j
avg
＝0.05。
[0241]
整个图像的j
ave
全局最小化是不实用的，因为有数百万种方法来分割图像。如果j值应用于类图的局部区域，可以很好地很好地指示该区域是在区域内部还是靠近区域边界。所以，可以考虑构造一个图像，其像素值对应于在以像素为中心的小窗口上计算的这些j值。将这些图像称为j图像，并将相应的像素值称为局部j值。局部j值越高，对应的像素越有可能靠近区域边界。j图像就像一个包含三维地形图的山谷和丘陵，它们实际上分别代表了区域内部和区域边界。
[0242]
小尺寸的窗口在定位强度/颜色边缘时很有用，而大尺寸的窗口对于检测纹理边界也很有用。通常，需要多个尺度来分割一个图像。在本发明实施例中，在最小尺度上的基本窗口是一个没有角的9
×
9呈圆形的窗口，如图8、图9所示。
[0243]
图8、图9是计算局部j值的窗口示意图。图8是计算局部j值的基本窗口；图9是2级窗口的降采样示意图，只有“+”点用于计算局部j值。
[0244]
2.3.2种子测定
[0245]
j图像的特点使得可以使用区域增长的方法来分割图像。该算法在一个粗糙的初始尺度上开始对图像进行分割。然后，在下一个更细的尺度上对新分割的区域重复相同的过程。区域生长包括确定种子点和从这些种子位置生长。区域增长之后是一个区域合并操作，得到最终的分割图像。
[0246]
一组初始种子区被确定为区域生长的基础之后，这些区域对应于局部j值的最小值。一般来说，找到最好的初始种子区域非常重要。以下简单的启发式方法在实验中提供了良好的结果：
[0247]
(1)计算该j图像的平均值和标准差，分别用μj和σj表示。
[0248]
(2)设置一个阈值变量tj：
[0249]
tj＝μj+aσjꢀꢀ
(2-25)
[0250]
其中，a是几个不同的预设值。将局部j值小于tj的像素作为候选种子点之后，基于4连通性连接候选种子点，获得候选种子区域。不同的阈值将得到不同的种子区及其数量，取数量最多的作为最终结果。
[0251]
(3)如果候选种子区域的尺寸大于相应尺度上表1中列出的最小尺寸，则确定为种子。
[0252]
表1不同尺度下的窗口尺寸
[0253]
尺度窗口(像素)采样(1/像素)最小种子大小(像素)19
×
91/(1
×
1)32217
×
171/(2
×
2)128333
×
331/(4
×
4)512465
×
651/(8
×
8)2048
[0254]
2.3.3种子增长
[0255]
确定了初始种子之后，从种子中生长出新的区域。种子逐个像素的慢慢生长。我们在实现中使用了一种更快的方法：
[0256]
(1)填充种子中的空洞。
[0257]
(2)平均该区域剩余未分配部分的局部j值，并连接低于平均值的像素，形成增长区域。如果一个生长区域毗邻一颗且只有一颗种子，它就被分配给该种子。
[0258]
(3)计算下一个较小尺度下剩余像素的局部j值，重复步骤2，以更准确地定位边界。
[0259]
(4)以最小的尺度逐个增长剩余的像素。种子边界处的未分类像素存储在一个缓冲区中。每次，将j值最小的像素分配给其相邻的种子，并更新缓冲区，直到对所有像素进行分类为止。
[0260]
2.4结果与分析
[0261]
实验用到的数据覆盖了2011年2月发生6.3级地震的区域，并在接下来的几年中重建。该数据集由2012年4月获得的航空图像组成，其中包含20.5km2范围内的12796栋建筑(2016年数据集为同一地区的16077栋建筑)。通过在地面上手动选择30个地面控制点(gcp)，进行地理校正，得到精度为1.6像素，分辨率为0.2米的航空数据集。
[0262]
从数据集选择分块序号为032的图像1和图像2进行实验，如图10和11所示。图像大小为512
×
512像素，为rgb 24位真彩色图像。
[0263]
j值的窗口大小集被设置为9
×
9、17
×
17、33
×
33和65
×
65像素，尺度数量m＝4。图12和13分别为图像1和图像2的33
×
33像素窗口的j图像。
[0264]
通过在不同尺度下分别对图像1和2进行图像分割，得到了图像1和图像2的分割对象序列，这些对象将用于变化检测。图14与图15为两幅图像在尺度3下的分割对象结果。通过视觉比较图10、图11与图14、图15，在尺度3下，方法可以较好地将地物建筑分割开来，分割边缘基本符合地物实际边缘。
[0265]
3、基于对象变化检测
[0266]
3.1相似性测度
[0267]
针对上述j图像的特征，需要根据分割结果分别对多尺度j图像中的每个对象进行分析和比较。此时，选择一个合适的相似度量来描述某一对象在不同时间内的相似性是至关重要的。常见的测量方法包括各种“距离”，如欧氏距离和马氏距离、直方图匹配、协方差等。
[0268]
3.1.1 ssim
[0269]
结构相似性(structural similarity，ssim)考虑了向量的平均值、方差和协方差，因此可以很好地表示向量之间的相似性。向量x和向量y之间的ssim由下式定义：
[0270]
s(x,y)＝[l(x,y)]
α
·
[c(x,y)]
β
·
[s(x,y)]
γ
ꢀꢀ
(3-1)
[0271]
其中：
[0272][0273][0274][0275]
式中，μ
x
、μy、σ
x
、σy、σ
2x
、σ
2y
分别为x和y的平均值、标准差和方差。σ
xy
是x和y之间的协方差。α、β和γ是三个向量的权值，c1、c2和c3是公式中的常数，以防止分母接近零时的不稳定。
[0276]
常取α＝β＝γ＝1，c3＝c2/2，则式(3-1)可以化简为：
[0277][0278]
s(x，y)越大，多时相图像之间的对象变化越小，相似度越高。此外，根据其定义，ssim具有以下特征：
[0279]
它有界s(x，y)∈[0，1]；它是对称的：s(x，y)＝s(y，x)；它具有唯一的最大值，当x＝y时，s(x，y)＝1。
[0280]
3.1.2psnr
[0281]
峰值信噪比(peak signal to noise ratio，psnr)，这个指标以分贝为单位计算信噪比。它通常用于对原始图像进行细化后的质量评估。psnr的值决定了细化图像的质量。较高的值意味着高质量，反之亦然。psnr由下式计算：
[0282][0283]
其中，max表示图像的最大像素数值，如果图像是8位灰度图像，那么max等于255；mse表示实际图像和预测/退化图像之间的均方误差，在两幅大小都是m
×
n的图像i和k中，mse由下式计算：
[0284][0285]
mse越小，表示两幅图像越相似。而mse在式(3-6)的分母位置，所以psnr的值越大，两幅图像越相似。与ssim相比，psnr不满足“有界”的特征，所以在进行变化检测之前要对这个相似性测定进行归一化。
[0286]
3.2多尺度融合与变化检测
[0287]
与单尺度检测相比，多尺度融合策略可以有效地提高检测精度，产生更可靠的结
果。融合策略采用了加权数据融合。wi(i＝1，m)表示每个尺度上检测结果的权重值。多尺度融合变化检测决策规则可以解释如下：
[0288]
(1)首先，对于每个对象rj，多尺度相似性sj通过由下式计算：
[0289][0290]
(2)然后，给定变化检测的相似性测定阈值ts，如果sj≥ts，则rj保持不变，标记为0。
[0291]
(3)否则，rj就被认为发生了变化，标记为1。
[0292]
(4)重复步骤(1)到(3)，直到完成分割结果中的所有对象。
[0293]
3.3景观格局指数过滤
[0294]
建筑的形状通常较为规整，多为矩形或矩形的组合。若能检测出形状较为规整的变化区域，则能过滤掉一些形状不规则地误检的非建筑区域，所以本发明引入景观格局指数。
[0295]
景观格局是景观的空间结构特征。景观格局指数反映了景观格局的结构组成和空间配置等方面的特征，通常包括面积、周长面积比、分维数等。本发明采用的分维数是反映景观物体形状复杂度的一个指标，其计算公式如下：
[0296][0297]
其中c为对象的周长、s为对象的面积。
[0298]
当分维数接近1时，代表其形状非常简单，如正方形；分维数越大代表形状越复杂。因此可以利用分维数分辨地面的建筑物、道路等人为的建筑。对每一个检测出来的变化区域对象，设定一个阈值t
idx
，当对象的分维数index大于该阈值时，则认为是非建筑。
[0299]
在融合策略中，设权值为w1＝0.1，w2＝0.2，w3＝0.4和w4＝0.3。在ssim相似性测度下，变化检测的阈值设置为0.5；在psnr测度下，归一化变化检测阈值设置为0.27。
[0300]
基于图10、图11的视觉观察和比较，得到两个测度下变化检测结果的分析标注图。分析如下：(1)ssim测度的变化检测对于新建建筑的检测不完整；psnr测度对于这部分则几乎完整地检测了出来。(2)对于两幅图像的b区域，ssim测度误把树木当成了变化区域，而psnr测度则没有误判。(3)ssim测度把停车场误认为是变化区域；psnr测度对于停车场的阴影部分比较敏感。(4)ssim测度把空地误判为变化区域。(5)ssim测度的变化检测图像出现条纹，并不能把整块建筑检测出来。(6)psnr测度把阴影区域检测为了变化区域。
[0301]
因此可以总结出：ssim的变化检测会出现不能将实际地物完整检测出来的现象，有时还会出现将不变区域误判为变化区域的情况。psnr的变化检测能够完整地检测出不变区域或者变化区域，但容易受到阴影的影响。
[0302]
在此基础上，基于建筑形状较为规整的特点，使用景观格局指数中的分维数对psnr的变化检测结果进行过滤，记为psnr+jg。分维数阈值t
idx
设为1.1。与psnr的变化检测相比，psnr+jg将变化结果小且不规则的区域和细长区域进行了去除，并进行了平滑。因此psnr+jg方法检测的变化区域更符合实际情况。
[0303]
通过计算总体精度、误报率、漏检率和kappa系数来定量评价两种测度评价方法和psnr+jg方法的性能，结果如表2所示。
[0304]
表2两种测度的多尺度融合变化检测与psnr+jg的检测精度
[0305][0306]
结合精度评价指标分析，在本实验中，ssim测度作为变化检测的相似性测度显然不合适。psnr相似性测度的变化检测效果则比ssim测度好很多。而psnr+jg的变化检测具有最高的总体精度和kappa系数，具有更低的误报率。
[0307]
为了分析变化检测对尺度的依赖性，进一步对检测结果及其精度参数进行比较。表3为上述检测结果的定量评价指标的比较。
[0308]
表3 psnr与psnr+jg的单一尺度变化检测结果检测精度
[0309][0310]
结合表3的精度评价，分析得出如下结论：(1)在一定阈值ts下，检测出的变化区域随着尺度的增大而减少。出现了尺度1多检而尺度4几乎检测不出变化区域的情况。(2)尺度3在所有单尺度变化检测中取得的效果最好。(3)无论是从视觉效果还是从精度评价来看，多尺度融合的效果最优。(4)增加jg操作能使尺度3或多尺度融合的总体精度更高。但jg操作有效的前提是检测结果区域较为完整，有较好的区域划分。
[0311]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所做的做的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。