CN101539629A - 基于多特征证据融合与结构相似度的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于多特征证据融合与结构相似度的遥感图像变化检测方法。其步骤：A.数据准备：选取需要进行变化检测的不同时相、相同空间范围的遥感图像。B.特征提取：分别提取出两时相遥感图像的纹理、梯度和边缘特征。C.结构相似度计算：计算两时相遥感图像在工作窗口内的纹理结构相似度、梯度结构相似度和边缘结构相似度。D.确定基本概率赋值函数：以纹理、梯度和边缘结构相似度为基础，构建证据理论的基本概率赋值函数(BPAF)。E.证据合成并获得变化区域，从而实现了变化检测。本发明方法基于图像边缘、纹理、梯度的区域统计特征进行变化检测，可有效避免仅使用单一特征进行变化检测造成的漏检与误检，算法具稳健性。

Description

基于多特征证据融合与结构相似度的遥感图像变化检测方法

技术领域

本发明涉及一种遥感图像变化内容的检测方法，具体说是一种基于多特征证据融合与结构相似度的遥感图像变化检测方法，属于遥感图像处理与信息提取领域。

背景技术

遥感图像变化检测是利用对同一地区在不同时期拍摄的两幅遥感数字图像(以下简称遥感图像)，采用数字图像处理和模式识别等手段，检测出该地区的地物变化信息的一项技术，它在环境、农业、水利和军事等诸多领域都有着非常广泛的应用。近年来，随着航天技术、传感器技术、计算机技术及其相关科学的迅猛发展，遥感技术得到了飞速进步，如GeoEye-1遥感图像的空间分辨率已达到0.41米。高空间分辨率遥感图像呈现地物的更多细节，光谱混淆现象更为严重，这给遥感图像变化检测带来了新的困难。

遥感图像变化检测主要依据的图像特征包括灰度、边缘、结构、形状、纹理等，其中最基本与最常用的是图像灰度特征，即每个像素上记录一组灰度值，它反映了地物的波谱反射特性。遥感技术应用的一个基本前提就是不同类型的地物具有不同的波谱反射特性，对于遥感图像变化检测来说就是地表覆盖类型的变化在不同时相的遥感图像上反映为波谱反射特性发生了变化，即导致了波谱反射值的变化。正是基于这种思路，传统的遥感图像变化检测主要依靠图像的灰度特征进行。它对于中低分辨率的遥感图像，效果较好，但对于光谱混淆现象严重的高空间分辨率遥感图像，其效果就较差。正是基于这种原因，以相对稳定的图像特征为依据进行遥感图像变化检测，日渐成为研究的热点。

相对于灰度特征，图像边缘、纹理、梯度等特征更为稳定，不易受遥感图像时相变化的影响。此外，不同的图像特征是反映图像不同侧面并互为补充的图像表征，如纹理描述了图像灰度分布的空间结构状况，梯度描述了区域内相邻像素的灰度变化程度，边缘则往往代表了目标边界。在变化检测时，单独使用某一特征进行变化检测都有可能造成漏检或误检，为此，应该考虑综合利用多特征实现变化检测。

曹广真等人在《电子与信息学报》2008年第30卷第8期提出了一种基于改进证据理论的多时相微波遥感图像融合与变化检测的方法。他们首先通过综合考虑证据本身的确信度和证据对辨别框架中子集的平均支持度进行证据间的加权合成，改进证据理论对证据的合成，提高其可靠性；然后提取不同时相图像间的散射幅度的对比度和概率密度分布函数在皮尔逊图中的距离，两种特征参数代表了像素级和区域级不同空间尺度下微波遥感图像中关于地表变化的信息；最后将改进的证据理论用于两特征之间的融合处理，得到地表的变化信息。曹等人的变化检测方法是针对微波遥感图像设计的，且没有使用包括图像边缘、纹理、梯度在内的其他重要图像特征。光学遥感图像的成像机理和微波遥感图像有很大不同，地物变化往往反映为图像边缘、纹理、形状等多特征的综合变化。为此，有必要发展综合利用图像更多特征，适用于光学遥感图像的新的变化检测方法。

发明内容

本发明的目的在于克服现有技术存在的缺陷，提出一种基于多特征证据融合与结构相似度的遥感图像变化检测方法。本方法对图像边缘变化、纹理变化、梯度变化等多变化证据进行证据融合(也称作证据合成)，实现了多种变化证据的综合利用，从而能有效地避免仅使用单一特征进行变化检测所造成的漏检或误检。相比现有技术，区别与特点在于：由于综合使用了图像的边缘、纹理、梯度等特征，所以方法更适合于高空间分辨率光学图像的遥感图像变化检测。

本发明方法原理与基本步骤

本发明涉及到的技术包括：1)图像特征(纹理、梯度、边缘)的提取；2)以上特征的直方图统计获得对应的特征向量；3)两时相遥感图像的特征向量的结构相似度计算；4)依托多特征、利用D-S证据理论，进行证据融合提取遥感图像变化区域。

在本方法中，图像边缘特征采用Canny边缘检测方法获取，纹理特征采用灰度共生矩阵方法(GLCM)获取，梯度特征采用Sobel算子获取。

Canny边缘检测方法步骤如下：

(1)用高斯滤波器对遥感图像进行平滑处理，得到平滑图像S(x，y)。在此需要设置高斯滤波器的标准差参数，以构造高斯滤波模板。

(2)用一阶偏导数的有限差分计算图像梯度的幅度和方向，得到图像f(x，y)上(i，j)点处的边缘强度A(i，j)和法向矢量

所有的A(i，j)构成边缘强度图像。

(3)对梯度幅值进行非极大值抑制：由于以图像f(x，y)上(i，j)点处的边缘强度A(i，j)值的大小并不能确定该点是否为边缘点，为此需要采用非极大值抑制方法确定边缘点。非极大值抑制的基本过程是：如果图像f(x，y)上(i，j)像素点的边缘强度A(i，j)小于沿着梯度线方向上的两个相邻像素点的边缘强度，则认为该像素点为非边缘点，将A(i，j)置为0，否则保留原值。

(4)用双阈值方法检测和连接边缘。其过程是首先设定两个阈值t1(低阈值)和t2(高阈值)(t2＞t1)，对非极大值抑制处理后的边缘强度图像进行双阈值化，可得到一个高阈值检测结果和一个低阈值检测结果。然后在高阈值检测结果中连接边缘轮廓，连接到端点时，到低阈值检测结果中寻找边缘点，直到高阈值检测结果所有间隙连接起来为止。在此需设置小于高阈值像素数占边缘强度图像的总像素数的比例，从而计算得到高阈值，并通过设置低阈值和高阈值之比来计算得到低阈值。

创建灰度共生矩阵的基本步骤是：在图像上任意取一点(x，y)，以及偏离它的另一点(x+a，y+b)，形成一个点对，设该点对的灰度值为(i，j)，即点(x，y)的灰度为i，点(x+a，y+b)的灰度为j。固定a和b，令点(x，y)在整幅图像上移动，则会得到各种(i，j)值。设灰度值的级数为L，则i和j的组合有L²种。对于整幅图像，统计出每一种(i，j)值的出现次数，再将它们归一化为出现的概率P_ij，则称方阵[P_ij]_L×L为灰度联合概率矩阵，也称为灰度共生矩阵。

在得到灰度共生矩阵之后，就可以利用它来计算一系列的纹理特征统计量，本方法选用灰度共生矩阵对比度特征来表征纹理，其计算公式如下：

$W=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{L-1}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}{(i-j)}^2---\left(1\right)$

其中W代表灰度共生矩阵对比度特征，L为灰度级数。

采用Sobel算子计算图像梯度。所采用的模板如下：

$\left|\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right|,$ $\left|\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right|.$

获取各种特征后，利用滑动窗口技术(我们称该滑动窗口为工作窗口，以下“两时相遥感图像工作窗口”均指在两个时相遥感图像的相同位置上同时开设的、相同大小的、并进行同步滑动的子窗口)，同时对两时相遥感图像进行特征统计。对于纹理、梯度特征，进行特征归一化以及分段(也称为特征量化)，构成频数直方图，以直方图的每个分段(量化等级)的频数形成特征向量；对于边缘特征，我们对其边缘分布模式进行分析，统计各种边缘模式的出现频率，构成频数直方图，形成特征向量(具体见方法详细步骤的说明)。

获取了各个特征的特征向量后，分别计算两个时相遥感图像工作窗口内的边缘、纹理、梯度的结构相似度作为多个特征的各自相似性指标，然后依此构建D-S证据理论所需的基本概率赋值函数(BPAF)。选用结构相似度作为特征向量相似性度量的理由如下。

为保证特征向量间相似性度量指标的通用性，我们要求其满足如下条件(以下X，Y代表特征向量)：

①有界性：即0≤|SIM(X，Y)|≤1，越接近于0，表示两向量相似性越弱；越

接近于1，表示两向量相似性越强；

②对称性：即SIM(X，Y)＝SIM(Y，X)；

③唯一最大值：即SIM(X，Y)＝1，当且仅当X＝Y；

在向量空间中定义的各种距离，例如欧式距离、马氏距离等，它们不满足有界性条件。对于相关系数和夹角余弦等指标，它们不满足唯一最大值条件，当两组向量在特征空间中“平行”时，其值都为1。对于直方图相交方法，它不满足对称性的条件，稳定性较差。

结构相似度从两个向量的均值、方差和协方差三个方面对两个向量的相似性进行度量，向量X与Y的结构相似度SSIM(X，Y)定义为：

SSIM(X，Y)＝[l(X，Y)]^α·[c(X，Y)]^β·[s(X，Y)]^γ    (2)

其中

$l(X,Y)=\frac{2{\mathrm{\μ}}_X{\mathrm{\μ}}_Y+C1}{{\mathrm{\μ}}_X^2+{\mathrm{\μ}}_Y^2+C2}$

$c(X,Y)=\frac{2{\mathrm{\σ}}_X{\mathrm{\σ}}_Y+C2}{{\mathrm{\σ}}_X^2{\mathrm{\σ}}_Y^2+C2}$

$s(X,Y)=\frac{{\mathrm{\σ}}_{\mathrm{XY}}+C3}{{\mathrm{\σ}}_X{\mathrm{\σ}}_Y+C3}$

其中μ_X、μ_Y、σ_X、σ_Y、σ_X ²、σ_Y ²、σ_XY分别是X与Y的均值、标准差、方差和协方差。α，β，γ＞0，它们是三个分量的权重，C1、C2、C3是为了防止当分母接近零时产生不稳定现象所添加的常数。当α＝β＝γ＝1，C3＝C2/2时，公式(2)可以简化为：

$\mathrm{SSIM}(X,Y)=\frac{(2{\mathrm{\μ}}_X{\mathrm{\μ}}_Y+C1)(2{\mathrm{\σ}}_{\mathrm{XY}}+C2)}{({\mathrm{\μ}}_X^2+{\mathrm{\μ}}_Y^2+C1)({\mathrm{\σ}}_X^2+{\mathrm{\σ}}_Y^2+C2)}---\left(3\right)$

结构相似度满足所需有界、对称、唯一性条件。而且相对于相关系数，结构相似度综合考虑了两变量的均值、方差和协方差，能更好地表征两变量之间的相似性。例如，对于向量X＝{53，17，59，37，4，69，27，53，29，49，11，67}，Y＝{123，87，129，107，74，139，97，123，99，119，81，137}，图1是两个向量的分布曲线图，经计算X和Y的相关系数为1.0，这违背了我们要求只有两个向量完全相等时其相似性才能为1.0的要求。经计算X和Y的结构相似度为0.64。由此可以看出结构相似度在度量特征向量相似性方面性能较好。

D-S证据理论作为一种数学工具，以其在不确定性的表示、量测和组合方面的优势而广受重视。它允许人们对不精确、不确定性问题进行建模，并进行推理，这为不确定信息的处理提供了新的思路。D-S证据理论的基本原理如下：称一个非空集合U为辨别框架，它由一系列互斥且穷举的元素组成。对于问题域中任意命题A，都应属于幂集2^U。在2^U上定义基本概率赋值函数m：2^U→[0，1]，使m满足

m(Φ)＝0                    (4)

$\underset{A\⊆U}{\mathrm{\Σ}}m\left(A\right)=1---\left(5\right)$

基本概率赋值函数m(A)表示证据对U的子集A的一种信任度量。m(A)的意义是：

(1)若 $A\⋐U,$ 且A≠U，则m(A)表示对A的确定信任度；

(2)若A＝U，则m(A)表示这个数不知如何分配；

(3)若 $A\⊆U,$ 且m(A)＞0，则称A是m的一个焦元。

对于不同的证据源，D-S证据理论采用正交和来组合这些证据源。设m₁，m₂，……，m_n为2^U上的n个基本概率赋值函数，它们的正交和表示为：

$m=m_1\⊕m_2\⊕...\⊕m_n---\left(6\right)$

且定义：

m(φ)＝0                        (7)

$m\left(A\right)=\frac{\underset{\∩A_i=A}{\mathrm{\Σ}}\underset{1\≤j\≤n}{\mathrm{\Π}}{m}_j\left(A_i\right)}{1-k},(\∀A\⋐U)---\left(8\right)$

其中

$k=\underset{\∩A_i=\mathrm{\φ}}{\mathrm{\Σ}}\underset{1\≤j\≤n}{\mathrm{\Π}}{m}_j\left(A_i\right)---\left(9\right)$

k反映了证据的冲突程度。式(7)、(8)(9)也称为Dempster合成法则。

本发明基于多特征证据融合与结构相似度的遥感图像变化检测方法，其步骤如下：

(1)数据准备

选取需要进行变化检测的两个不同时相、相同空间范围的遥感图像。

(2)特征提取

分别提取出两时相遥感图像的纹理、梯度和边缘特征。分别计算两时相遥感图像的灰度共生矩阵的对比度特征来描述两时相遥感图像的纹理特征，使用Canny边缘检测方法计算两时相遥感图像的边缘特征，使用Sobel算子计算两时相遥感图像的梯度特征。

(3)结构相似度计算

在两时相遥感图像的相同位置上同时开设相同大小的工作窗口，并进行逐行列同步滑动，计算两时相遥感图像在工作窗口内的纹理结构相似度、梯度结构相似度和边缘结构相似度。

对于纹理结构相似度，若直接用工作窗口内的原始特征值构成的特征向量进行结构相似度计算，它们对图像的配准精度要求较高，且易受噪声影响。为此采用直方图统计的方式对特征进行量化，量化公式为：

$M=\frac{\mathrm{Max}-\mathrm{Min}}{L}---\left(10\right)$

$h=\frac{x-\mathrm{Min}}{M}$

其中：x为当前要量化的值，Max和Min是分别是两时相工作窗口内纹理特征的最大值和最小值，L为量化等级的个数，M为每个量化等级的宽度(均相同)，h为纹理特征值经过量化后的等级归属。将两时相工作窗口内的纹理特征分别进行量化，并统计其各个量化等级的相应频数，生成各自的直方图，构成各自的纹理特征向量，然后用公式(3)计算两时相遥感图像在该工作窗口内的纹理特征向量的纹理结构相似度S_t，此时公式(3)中，X，Y代表的是两时相工作窗口内的纹理特征向量。

对于梯度相似度，采用和纹理相似度相同的方式计算。同样利用公式(10)进行量化。此时公式(10)中的x代表当前要量化的梯度特征值，Max和Min是分别是两时相工作窗口内梯度特征的最大值和最小值，L为量化等级，M为每个量化等级的宽度(均相同)，h为梯度特征值经过量化后的等级归属。将两时相工作窗口内的梯度特征分别进行量化，并统计其各个量化等级的频数，生成各自的直方图，构成各自的梯度特征向量，然后同样用公式(3)计算两时相遥感图像在该工作窗口内的梯度特征向量的梯度结构相似度S_g。此时公式(3)中的X，Y代表的是两时相工作窗口内的梯度特征向量。

由于边缘特征已经二值化而不能沿用上面的方法。所以我们设计了边缘模式分布直方图来计算边缘特征的结构相似度。边缘模式分布直方图反映了在某个区域内，图像边缘点分布的各种组合模式与其出现频率间的统计关系。在这里，在工作窗口内再开设一个2×2的子窗口用于特征统计。由于边缘是二值化的，如用黑白两种颜色代表的话(我们用黑色代表背景，白色代表边缘像素)，在2×2的窗口内，根据排列组合知识可得到如图3所示的14种情形。以这14种情形作为模板，对两个时相遥感图像的对应工作窗口内的边缘图像分别进行逐行列遍历，并记录每个模板匹配的次数，就得到了两个时相遥感图像的工作窗口内的边缘模式分布直方图，构成各自的边缘特征向量。然后用式(3)计算两个边缘特征向量的结构相似度，作为边缘结构相似度S_e。此时公式中的X，Y代表的是两时相工作窗口内的边缘特征向量。

(4)确定基本概率赋值函数

以纹理结构相似度、梯度结构相似度和边缘结构相似度为基础，构建证据理论的基本概率赋值函数(BPAF)。

首先确定证据理论的辨别框架U。将遥感图像的所有像素作为检测对象，定义辨别框架U为：

U＝{Y，N}

其中Y表示变化类，N表示未变化类，因此2^U的非空子集为：{Y}，{N}，{Y，N}。

纹理特征的基本概率赋值函数(BPAF1)为：

m₁({Y})＝(1.0-S_t)*α₁                (11)

m₁({N})＝S_t*α₁                      (12)

m₁({Y，N})＝1.0-α₁                  (13)

其中α₁为该证据(纹理)对判别的信任度，该值依经验给出。

梯度特征的基本概率赋值函数(BPAF2)为：

m₂({Y})＝(1.0-S_g)*α₂                (14)

m₂({N})＝S_g*α₂                      (15)

m₂({Y，N})＝1.0-α₂                  (16)

其中α₂为该证据(梯度)对判别的信任度，该值依经验给出。

边缘特征的基本概率赋值函数(BPAF3)为：

m₃({Y})＝(1.0-S_e)*α₃                (17)

m₃({N})＝S_e*α₃                      (18)

m₃({Y，N})＝1.0-α₃                  (19)

其中α₃为该证据(边缘)对判别的信任度，该值依经验给出。

(5)证据合成并获得变化区域

获得两时相遥感图像某工作窗口内的边缘、纹理、梯度的各自BPAF后，采用公式(8)对BPAF进行证据合成，得到合成后的BPAF(BPAF4)。最后对未变化类(N)的BPAF4和变化类(Y)的BPAF4分别设置如下阈值，判断某个像素是否发生变化：支持变化类的BPAF4的值大于0.25或者支持未变化类的BPAF4的值小于0.7。如满足条件，该像素是变化像素，否则为非变化像素。在整幅遥感图像上滑动工作窗口，对整幅图像的所有像素均进行以上处理，则得到所有变化区域，从而实现了变化检测。

本发明方法具有以下优点：

(1)方法鲁棒性高。传统的变化检测方法，对配准精度和辐射校正的精度要求较高，该算法利用窗口特征统计技术计算两时相图像的相似性，为此对预处理误差具有较强的免疫力；此外，由于该算法不是直接基于图像灰度，而是基于较为稳定的边缘、纹理、梯度的区域统计特征进行变化检测，算法具稳健性。

(2)从边缘变化、纹理变化和梯度变化多个证据源入手进行变化检测，有效避免了仅使用单一特征进行变化检测的漏检与误检。

(3)采用结构相似度作为相似性匹配准则，它综合了向量的均值、方差和协方差等三个因素，对相似性进行度量，能较好的表达相似性。

附图说明

图1：向量分布曲线

图2：基于多特征证据融合与结构相似度的遥感图像变化检测方法流程图

图3：2×2窗口内14种边缘像素的分布模式(黑色代表背景，白色代表边缘像素)

图4：实施区遥感图像及变化检测结果：(a)实施区03年ALOS图像，(b)实施区05年ALOS图像，(c)使用边缘、纹理、梯度三种特征的检测结果，(d)使用梯度与边缘特征的检测结果，(e)单独使用梯度特征的检测结果，(f)单独使用边缘特征的检测结果，(g)单独使用纹理特征的检测结果

图5：纹理直方图示意：(a)第一时相某工作窗口纹理直方图，(b)第二时相某工作窗口纹理直方图

图6：梯度直方图示意：(a)第一时相某工作窗口梯度直方图，(b)第二时相某工作窗口梯度直方图

图7：边缘模式直方图示意：(a)第一时相某工作窗口边缘模式直方图，(b)第二时相某工作窗口边缘模式直方图

具体实施方式

下面结合附图和实施例，对本发明作详细说明。

实施例

(1)数据准备

实验图像为南京江宁的某区域的ALOS全色图像，空间分辨率2.5米。成像时间分别为2003年11月和2005年12月。图像的坐标空间范围是31°39′49.84″N~31°40′30.01″N，119°2′9.66″E~119°3′44.74″E。图像大小是1012×477像素，如图4(a)、4(b)所示。

(2)特征提取

分别计算两时相遥感图像灰度共生矩阵的对比度指标来获得两时相遥感图像的纹理特征，其中灰度级数为32，点对间的偏离方向取0°，偏离步长是1个像元。使用Canny边缘检测方法计算两时相遥感图像的边缘特征，高斯滤波使用的标准差为0.6，双阈值法中的低阈值和高阈值之比为0.7，小于高阈值像素数占像素总数的比例为0.8。使用Sobel算子计算两时相图像的梯度特征，采用的模板见方法原理与步骤小节。

(3)结构相似度计算

在两时相遥感图像上同时开设工作窗口，窗口大小为11×11像素，并进行逐行列移动，并计算两时相图像在工作窗口内的纹理、梯度和边缘结构相似度。

将两时相工作窗口内的纹理特征通过公式(10)进行量化，量化等级L为14。统计其各个量化等级的频数，生成各自的直方图(见图5示意)，构成各自的纹理特征向量，然后用公式(3)计算该工作窗口内两时相遥感图像的纹理特征向量的纹理结构相似度S_t。在计算纹理结构相似度时，C1和C2分别为0.3和0.9。

将工作窗口内的两时相梯度特征通过公式(10)进行量化，量化等级L为14。将两个工作窗口内的梯度特征分别进行量化，并统计其各个量化等级的频数，生成各自的直方图(见图6示意)，构成各自的梯度特征向量，然后同样用公式(3)计算该工作窗口内两时相遥感图像的梯度特征向量的梯度结构相似度S_g。在计算梯度结构相似度时，C1和C2分别为0.3和0.9。

采用边缘模式分布直方图来计算边缘特征的结构相似度。边缘模式分布直方图反映了在某个区域内，图像边缘点分布的各种组合模式与其出现频率间的统计关系。在这里，在工作窗口内再开设一个2×2的子窗口用于特征统计。由于边缘是二值化的，如用黑白两种颜色代表的话(我们用黑色代表背景，白色代表边缘像素)，在2×2的窗口内，根据排列组合知识可得到如图3所示的14种情形。以这14种情形作为模板，对两个时相工作窗口内的边缘图像分别进行遍历，并记录每个模板匹配的个数，就得到了相应工作窗口的边缘模式分布直方图(见图7示意)，构成各自的边缘特征向量。然后用公式(3)计算两个边缘特征向量的结构相似度，作为边缘结构相似度S_e，在计算边缘结构相似度时，C1和C2分别为0.3和0.9。

(4)确定基本概率赋值函数

以纹理、梯度和边缘结构相似度为基础，构建证据理论的纹理BPAF、梯度BPAF、边缘BPAF。对纹理、梯度和边缘三个BPAF进行证据融合，然后进行规则判断，得到变化检测结果。具体步骤如下。

首先确定证据理论的辨别框架U。将图像的所有像素作为检测对象，定义辨别框架U为：

U＝{Y，N}

其中Y表示变化类，N表示未变化类，因此2^U的非空子集为：{Y}，{N}，{Y，N}。纹理BPAF(BPAF1)为式(11)、(12)和(13)，其中α₁为0.9。梯度BPAF(BPAF2)为(14)、(15)和(16)，其中α₂为0.9。边缘BPAF(BPAF3)为(17)、(18)和(19)，其中α₃为0.95。

(5)证据合成并获得变化区域

获得两时相遥感图像某工作窗口内的边缘、纹理、梯度的各自BPAF后，采用公式(8)对BPAF进行证据合成，得到合成后的BPAF(BPAF4)。最后对未变化类(N)的BPAF4和变化类(Y)的BPAF4分别设置如下阈值，判断某个像素是否发生变化：支持变化类的BPAF4的值大于0.25或者支持未变化类的BPAF4的值小于0.7。如满足条件，该像素是变化像素，否则为非变化像素。在两时相遥感图像上同步滑动工作窗口，对整幅图像的所有像素均进行以上处理，则得到所有变化区域，从而实现了变化检测。

图4是实施区的变化检测结果。其中4(c)是使用纹理、梯度和边缘三个特征进行变化检测的结果，4(d)是使用梯度和边缘二个特征进行变化检测的结果，4(e)是单独使用梯度特征的变化检测结果，4(f)是单独使用边缘特征的变化检测结果，4(g)是单独使用纹理特征的变化检测结果。由此可以看出，采用多特征证据融合的方式进行变化检测，相对于使用单一特征的变化检测，检测精度得到了有效提高。此外，梯度、纹理和边缘特征三者进行证据融合的检测方法，相对梯度、边缘两个特征进行证据融合的方法，对变化相对敏感，可用来检测一些相对细微的图像变化。

Claims (1)

1、基于多特征证据融合与结构相似度的遥感图像变化检测方法，其步骤如下：

A、数据准备

选取需要进行变化检测的两个不同时相、相同空间范围的遥感图像；

B、特征提取

分别提取出两时相遥感图像的纹理特征、梯度特征和边缘特征；即：分别计算两时相遥感图像的灰度共生矩阵的对比度特征来描述两时相遥感图像的纹理特征，使用Canny边缘检测方法计算两时相遥感图像的边缘特征，使用Sobel算子计算两时相遥感图像的梯度特征；

C、结构相似度计算

在两时相遥感图像的相同位置上同时开设相同大小的工作窗口，并进行逐行列同步滑动，计算两时相遥感图像在工作窗口内的纹理结构相似度、梯度结构相似度和边缘结构相似度；

对于纹理结构相似度，采用直方图统计的方式对特征进行量化，量化公式为：

$M=\frac{\mathrm{Max}-\mathrm{Min}}{L}---\left(a\right)$

$h=\frac{x-\mathrm{Min}}{M}$

其中：x为当前要量化的值，Max和Min是分别是两时相工作窗口内纹理特征的最大值和最小值，L为量化等级的个数，M为每个量化等级的宽度(均相同)，h为纹理特征值经过量化后的等级归属；

将两时相工作窗口内的纹理特征分别进行量化，并统计其各个量化等级的相应频数，生成各自的直方图，构成各自的纹理特征向量，然后用下列公式计算两时相遥感图像在该工作窗口内的纹理特征向量的纹理结构相似度S_t：

$\mathrm{SSIM}(X,Y)=\frac{(2{\mathrm{\μ}}_X{\mathrm{\μ}}_Y+C1)(2{\mathrm{\σ}}_{\mathrm{XY}}+C2)}{({\mathrm{\μ}}_X^2+{\mathrm{\μ}}_Y^2+C1)({\mathrm{\σ}}_X^2+{\mathrm{\σ}}_Y^2+C2)}---\left(b\right)$

其中X，Y代表特征向量。μ_X、μ_Y、σ_X、σ_Y、σ_X ²、σ_Y ²、σ_XY分别是X与Y的均值、标准差、方差和协方差，C1、C2是为了防止当分母接近零时产生不稳定现象所添加的常数；在这里公式(b)中的X，Y代表的是两时相工作窗口内的纹理特征向量；

对于梯度相似度，采用和纹理相似度相同的方式计算；即：利用公式(a)进行量化，此时公式(a)中的x代表当前要量化的梯度特征值，Max和Min是分别是两时相工作窗口内梯度特征的最大值和最小值，L为量化等级，M为每个量化等级的宽度，h为梯度特征值经过量化后的等级归属；将两时相工作窗口内的梯度特征分别进行量化，并统计其各个量化等级的频数，生成各自的直方图，构成各自的梯度特征向量，然后同样用公式(b)计算两时相遥感图像在该工作窗口内的梯度特征向量的梯度结构相似度S_g；此时公式(b)中的X，Y代表的是两时相工作窗口内的梯度特征向量；

采用边缘模式分布直方图来计算边缘特征的结构相似度；即：在工作窗口内再开设一个2×2的子窗口用于特征统计；由于边缘是二值化的，采用黑色代表背景，白色代表边缘像素，在2×2的窗口内，根据排列组合知识可得到14种排列组合情形；以这14种情形作为模板，对两个时相遥感图像的对应工作窗口内的边缘图像分别进行逐行列遍历，并记录每个模板匹配的次数，就得到了两个时相遥感图像的工作窗口内的边缘模式分布直方图，构成各自的边缘特征向量；然后用公式(b)计算两个边缘特征向量的结构相似度，作为边缘结构相似度S_e；此时公式(b)中的X，Y代表的是两时相工作窗口内的边缘特征向量；

D、确定基本概率赋值函数

以纹理结构相似度、梯度结构相似度和边缘结构相似度为基础，构建证据理论的基本概率赋值函数(BPAF)：

首先确定证据理论的辨别框架U，将遥感图像的所有像素作为检测对象，定义辨别框架U为：

U＝{Y，N}

其中Y表示变化类，N表示未变化类，因此2^U的非空子集为：{Y}，{N}，{Y，N}；

纹理特征的基本概率赋值函数(BPAF1)为：

m₁({Y})＝(1.0-S_t)*α₁        (c)

m₁({N})＝S_t*α₁              (d)

m₁({Y，N})＝1.0-α₁          (e)

其中α₁为该证据对判别的信任度；

梯度特征的基本概率赋值函数(BPAF2)为：

m₂({Y})＝(1.0-S_g)*α₂        (f)

m₂({N})＝S_g*α₂              (g)

m₂({Y，N})＝1.0-α₂          (h)

其中α₂为该证据对判别的信任度；

边缘特征的基本概率赋值函数(BPAF3)为：

m₃({Y})＝(1.0-S_e)*α₃        (i)

m₃({N})＝S_e*α₃              (j)

m₃({Y，N})＝1.0-α₃        (k)

其中α₃为该证据对判别的信任度；

E、证据合成并获得变化区域

获得两时相遥感图像某工作窗口内的边缘、纹理、梯度的各自BPAF后，采用下列公式对BPAF进行证据合成，得到合成后的BPAF(BPAF4)：

$m=m_1\⊕m_2\⊕m_3---\left(1\right)$

其中符号

代表D-S证据理论的正交和运算；

最后对未变化类(N)的BPAF4和变化类(Y)的BPAF4分别设置如下阈值，判断某个像素是否发生变化：

支持变化类的BPAF4的值大于0.25或者支持未变化类的BPAF4的值小于0.7；如满足条件，该像素是变化像素，否则为非变化像素；在整幅遥感图像上滑动工作窗口，对整幅图像的所有像素均进行以上处理，则得到所有变化区域，从而实现了变化检测。

Images