CN103955943A - 基于融合变化检测算子与尺度驱动的无监督变化检测方法 - Google Patents

Abstract

本发明公开了一种基于融合变化检测算子与尺度驱动的无监督变化检测方法，主要解决现有技术通常利用滤波来消除图像噪声导致大量图像细节丢失，增加变化检测难度的问题。其实现步骤是：(1)利用变化检测的图像数据生成对数比图像和均值比图像；(2)将对数比图像和均值比图像进行融合生成初步差异图；(3)判断初步差异图上每个像素点可靠的尺度，并将初步差异图上每个像素点所有可靠的尺度进行融合，生成变化检测差异图；(4)对变化检测差异图进行聚类，生成变化检测结果图。本发明在有效抑制噪声的同时，保留了更多的图像细节，提高了变化检测结果的精确度，可用于对SAR图像的变化检测。

Description

基于融合变化检测算子与尺度驱动的无监督变化检测方法

技术领域

本发明属于图像处理技术领域，特别涉及无监督的变换检测的方法，可用于对SAR图像的变化检测。

背景技术

变化检测属于图像处理领域，是指从不同时期的遥感数据中定量分析和确定地表变化的特征与过程。通常处理的是两时相的变化检测问题，其实质是地表特征随时间变化发生的变化而引起的两个时期影像像元光谱响应的变化。现代遥感技术飞速发展，而微波遥感由于具有全天候、全天时的工作能力，对地物具有一定的穿透能力，而且采用侧视方式成像，覆盖面积大，与可见光和红外遥感相比，具有无可比拟的优点。因此，SAR图像日益成为变化检测的主要数据来源。

通常，两时相的SAR图像变化检测的工作流程为：数据选取、数据预处理、变化信息提取和分离、精度评定。其中，数据选取指的是变化前后的遥感图像的选取，数据预处理通常指的是几何校正和辐射纠正。由于SAR图像固有的乘性噪声的影响，基于SAR图像的变化检测通常需要滤波，以消除斑点噪声的影响，所以，数据预处理中也常常包括滤波的过程。预处理之后的无监督的SAR图像变换检测主要包括两个步骤：1.生成差异图；2.分析差异图。其中，生成差异图中常用的算子有差分算子、比值算子、对数比算子、均值算子等；分析差异图通常使用阈值算法或者聚类算法。

以上为常用的SAR图像无监督变化检测方法，存在如下缺点：1.预处理过程中的滤波过程虽然能够有效地消除乘性噪声的影响，却势必造成图像细节信息的丢失，不利于后续对图像信息的提取，对于某些噪声影响差别巨大的两时相图像处理效果尤其不好。2.由于SAR图像的复杂性，单一算子处理得到的差异图像通常具有一定的局限性。3.直接对差异图算子生成的图像进行分析，最终得到的变化检测效果常常依赖于选定的差异图分析方法，且不同的分析方法造成的结果差异巨大。

发明内容

本发明的目的在于针对上述已有技术的问题，提供一种基于融合变化检测算子与尺度驱动的无监督变化检测方法，以减小图像细节信息的丢失，提高变化检测的精确度。

为实现上述目的，本发明的技术方案包括如下：

(1)获取同一地区经过配准和几何校正后的两幅不同时间的图像M1和M2；

(2)对两幅不同时间的图像M1和M2，应用对数比检测算子得到对数比图像ML，应用均值比检测算子得到均值比图像MR；

(3)根据融合规则，将对数比图像ML和均值比图像MR进行融合，得到变化检测的初步差异图E；

(4)对得到的变化检测的初步差异图E进行多尺度二维离散小波分解，分解尺度为N，在每个尺度上舍去高频信息后进行图像重构，得到该尺度的重构图像，由每个尺度上重构图像得到一组图像E'＝{E₁,E₂,...,E_N}；

(5)对于每个尺度上的重构图像，计算其整幅局域变化系数C，以及该重构图像上每个像素点邻域内的局域变化系数L，并将该整幅局域变化系数C与局域变化系数L进行比较，如果C≥L，则判定该重构图像上的像素点属于匀质区域的内部，否则，不属于匀质区域的内部；

(6)利用步骤(5)判断初步差异图E上每一个像素点可靠的尺度：如果存在整数S，1＜S≤N，使得从尺度1到尺度S的重构图像上，位于同一位置的像素点均属于匀质区域的内部，则判定该像素点可靠的尺度为尺度1到尺度S，否则，判定该像素点可靠的尺度为尺度1，并将每个像素点所有可靠的尺度进行融合，最终得到变化检测差异图F；

(7)使用FCM_S1聚类算法对变化检测差异图F进行聚类，得到变化检测结果图M。

本发明与现有技术相比具有以下优点：

1)本发明由于将对数比图像和均值比图像进行融合得到差异图，用该融合差异图与仅使用对数比图像或均值比图像得到的差异图相比，能够保留更多的图像信息；

2)本发明由于对差异图进行多尺度二维离散小波分解，在每个尺度上舍去高频信息后进行图像重构，使高频信息的舍去能够有效抑制乘性噪声对图像的影响，因此不需要对图像进行滤波预处理，能够更好地保持图像的细节。

附图说明

图1是本发明的实现流程图；

图2是常用的无监督SAR图像变化检测方法的实现流程图；

图3是本发明仿真实验中所使用的一组SAR图像数据；

图4是本发明仿真实验得到的变化检测结果图；

图5是对滤波后的图像使用对数比算子进行仿真实验得到的变化检测结果图；

图6是对滤波后的图像使用均值比算子进行仿真实验得到的变化检测结果图；

图7是对滤波后的图像使用融合检测算子进行仿真实验得到的变化检测结果图；

图8是本发明仿真实验中所使用的一组SAR图像数据的参考图像。

具体实施方式

以下结合附图，对本发明做进一步的详细描述。

参照图1，本发明的实现步骤如下：

步骤1：获取同一地区经过配准和几何校正后的两幅不同时间的图像M1和M2。

步骤2：对两幅不同时间的图像M1和M2，应用对数比检测算子得到对数比图像ML，应用均值比检测算子得到均值比图像MR。

2a)通过对数比检测算子生成对数比图像，对数比检测算子为

ML＝|log(M1/M2)|，

其中M1和M2分别为获取的两幅不同时间的图像；

2b)通过均值比检测算子生成均值比图像，均值比检测算子为

MR(i,j)＝1-min(M1(i,j)/M2(i,j),M2(i,j)/M1(i,j))，

其中M1(i,j)为第一幅图像M1中以(i,j)位置处的像素点为中心的邻域均值,M2(i,j)为第二幅图像M2中以(i,j)位置处的像素点为中心的邻域均值。

步骤3：将对数比图像ML和均值比图像MR按如下公式进行融合，得到变化检测的初步差异图E：E＝(ML×MR)/(1-ML-MR+2×ML×MR)。

步骤4：利用多尺度二维离散小波构造一组关于初步差异图E的重构图像E'＝{E₁,E₂,...,E_N}。

4a)对初步差异图E进行多尺度二维离散小波分解，本实例取分解尺度为7，在每个尺度上先产生四个分解系数矩阵，分别为低频部分分解系数LL，垂直方向分解系数HL，水平方向分解系数LH，对角线方向系数HH；再将对角线方向分解系数设置为相同大小的零矩阵Z；

4b)在每一个尺度上，使用低频部分分解系数LL，垂直方向分解系数HL，水平方向分解系数LH以及零矩阵Z进行二维离散小波重构，得到该尺度的重构图像，由所有尺度上的重构图像构成一组关于初步差异图E的重构图像。

步骤5：对于步骤4得到的一组重构图像，逐图判断图像上的每个像素点是否属于匀质区域的内部。

5a)对于每幅图像，计算其整幅局域变化系数：C＝σ/μ，

其中，σ为整幅图像的标准差，μ为整幅图像的均值；

5b)计算每幅图像上每个像素点邻域内的局域变化系数：L＝σ(i,j)/μ(i,j)，

其中，σ(i,j)为图像中以(i,j)位置处的像素点为中心的邻域的标准差，μ(i,j)为图像中以(i,j)位置处的像素点为中心的邻域的均值；

5c)将整幅局域变化系数C与局域变化系数L进行比较，如果C≥L，则判定该重构图像上的像素点属于匀质区域的内部，否则，不属于匀质区域的内部。

步骤6：对于初步差异图E上的每一个像素点，根据步骤5判断该像素点可靠的尺度，并将其所有可靠的尺度进行融合，得到变化检测差异图F。

6a)如果存在整数S，1＜S≤N，且在从尺度1到尺度S的重构图像上，位于同一位置的像素点均在匀质区域的内部，则判定该位置像素点可靠的尺度为尺度1到尺度S，否则，判定其可靠的尺度仅为尺度1；

6b)按照如下公式对每个像素点所有可靠的尺度进行融合，得到变化检测差异图F：

其中，F(i,j)为变化检测差异图F中位于(i,j)处的像素点，E(i,j)为初步差异图E中位于(i,j)处的像素点，E_h(i,j)为初步差异图E的重构图像E_h中位于(i,j)处的像素点，h为分解的尺度。

步骤7：使用FCM_S1聚类算法对变化检测差异图F进行聚类，得到变化检测结果图M。

7a)设置聚类数目c＝2，模糊指数m＝2，阈值ε＝10-6；

7b)随机选取两个聚类中心V₁ ⁰和V₂ ⁰，设置循环次数t＝1；

7c)更新隶属度矩阵U＝{U_ik,1≤i≤2,1≤k≤n}，其中：

$U_{\mathrm{ik}}=\frac{{({\left|\right|X_k-V_i\left|\right|}^2+\mathrm{\α}{\left|\right|{\stackrel{\‾}{X}}_k-V_i\left|\right|}^2)}^{-\frac{1}{(m-1)}}}{\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{({\left|\right|X_k-V_j\left|\right|}^2+\mathrm{\α}{\left|\right|{\stackrel{\‾}{X}}_k-V_j\left|\right|}^2)}^{-\frac{1}{(m-2)}}},$

其中，n是图像像素点的总个数，X_k是第k个像素点的灰度值,是以第k个像素点为中心的邻域均值，V_i是聚类中心，参数α用来调节原始图像和相关的均值滤波图像的权重；

7d)更新聚类中心V＝{V_i,i＝1,2}，其中：

$V_i=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{U}_{\mathrm{ik}}^m(X_k+\mathrm{\α}{\stackrel{\‾}{X}}_k)}{(1+\mathrm{\α})\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{U}_{\mathrm{ik}}^m};$

7e)计算本次迭代生成的聚类中心与上次迭代生成的聚类中心的距离d＝||V^t-V^t-1||，如果满足d＜ε，则停止迭代，并由隶属度矩阵U得到各像素点的分类结果，否则，令t＝t+1返回步骤7c)。

本发明的效果可以通过以下仿真实验进一步说明。

1.仿真实验条件。

本实验采用软件MATLAB7.10.0作为仿真工具，计算机配置为Intel Core i5/2.27G/2G，仿真实验中采用的数据集取自中国山东省东营市黄河入海口数据集，如图3所示。其中3(a)拍摄于2008年6月18日，3(b)拍摄于2009年6月19日。

2.仿真实验内容。

在上述实验条件下，用本发明方法和三种常用的无监督SAR图像变化检测方法对图3所示的数据集进行仿真实验。三种常用方法具有相同的流程，如图2。该三种方法分别使用不同的变化检测算子：1)对数比检测算子，2)均值比检测算子，3)融合检测算子。本仿真实验中选用的滤波方法均为Lee滤波，聚类方法均为FCM_S1聚类。

仿真实验1，使用本发明对图3所示的一组图像进行变化检测，结果如图4。

仿真实验2，使用对数比检测算子对图3所示的一组图像进行变化检测，结果如图5。

仿真实验3，使用均值比检测算子对图3所示的一组图像进行变化检测，结果如图6。

仿真实验4，使用融合检测算子对图3所示的一组图像进行变化检测，结果如图7。

从图4与图5、6、7的对比可见，本发明的结果图4具有最优的视觉效果，对噪声具有更好的鲁棒性。

为了进一步证明本发明方法的优越性，将仿真实验得到的图4、图5、图6、图7分别与图8进行比较，使用客观评价指标对生成的变化检测结果图进行定量评价，评价结果如表1所示。常用的评价指标有漏检数FN，错检数FP，总错误数OE，正确率PCC以及kappa系数。

表1仿真实验结果图的定量评价结果

仿真实验	FN	FP	OE	PCC	Kappa
						1	191	112	303	0.9976	0.8830
2	263	27484	27747	0.7798	0.0532
						3	12	34135	34147	0.7290	0.0531
4	392	354	746	0.9941	0.7163

从表1可以看出，与其它三个仿真实验2、3、4相比，仿真实验1的结果图具有最少的总错误数，最高的正确率以及最高的Kappa系数，表现最优。

综上所述，本发明方法无论在主观效果还是客观指标上，都表现出了较好的性能，其在有效抑制噪声的同时，保留了更多的图像细节，提高了变化检测结果的精确度。

Claims (7)

1.一种基于融合变化检测算子与尺度驱动的无监督变化检测方法，包括如下步骤：

(1)获取同一地区经过配准和几何校正后的两幅不同时间的图像M1和M2；

(2)对两幅不同时间的图像M1和M2，应用对数比检测算子得到对数比图像ML，应用均值比检测算子得到均值比图像MR；

(3)根据融合规则，将对数比图像ML和均值比图像MR进行融合，得到变化检测的初步差异图E；

(4)对得到的变化检测的初步差异图E进行多尺度二维离散小波分解，分解尺度为N，在每个尺度上舍去高频信息后进行图像重构，得到该尺度的重构图像，由每个尺度上重构图像得到一组图像E'＝{E₁,E₂,...,E_N}；

(5)对于每个尺度上的重构图像，计算其整幅局域变化系数C，以及该重构图像上每个像素点邻域内的局域变化系数L，并将该整幅局域变化系数C与局域变化系数L进行比较，如果C≥L，则判定该重构图像上的像素点属于匀质区域的内部，否则，不属于匀质区域的内部；

(6)利用步骤(5)判断初步差异图E上每一个像素点可靠的尺度：如果存在整数S，1＜S≤N，使得从尺度1到尺度S的重构图像上，位于同一位置的像素点均属于匀质区域的内部，则判定该像素点可靠的尺度为尺度1到尺度S，否则，判定该像素点可靠的尺度为尺度1，并将每个像素点所有可靠的尺度进行融合，最终得到变化检测差异图F；

(7)使用FCM_S1聚类算法对变化检测差异图F进行聚类，得到变化检测结果图M。

2.根据权利要求1所述的方法，其中所述步骤(2)中的对数比检测算子，表示为：

ML＝|log(M1/M2)|，

其中M1和M2分别为获取的同一地区经过配准和几何校正后的两幅不同时间的图像。

3.根据权利要求1所述的方法，其中所述步骤(2)中的均值比检测算子，表示为：

MR(i,j)＝1-min(M1(i,j)/M2(i,j),M2(i,j)/M1(i,j))，

其中M1(i,j)为第一幅图像M1中以(i,j)位置处的像素点为中心的邻域均值，M2(i,j)为第二幅图像M2中以(i,j)位置处的像素点为中心的邻域均值。

4.根据权利要求1所述的方法，其中步骤(3)所述的根据融合规则，将对数比图像ML和均值比图像MR进行融合，按如下公式进行：E＝(ML×MR)/(1-ML-MR+2×ML×MR)，其中，E为变化检测的初步差异图。

5.根据权利要求1所述的方法，其中步骤(5)所述的对于每个尺度上的重构图像，计算其整幅局域变化系数C，通过如下公式计算：

C＝σ/μ，

其中，σ为整幅图像的标准差，μ为整幅图像的均值。

6.根据权利要求1所述的方法，其中所述步骤(5)中计算重构图像上每个像素点邻域内的局域变化系数L，通过如下公式计算：

L＝σ(i,j)/μ(i,j)，

其中，σ(i,j)为重构图像中以(i,j)位置处的像素点为中心的邻域的标准差，μ(i,j)为重构图像中以(i,j)位置处的像素点为中心的邻域的均值。

7.根据权利要求1所述的方法，其中步骤(5)所述的使用FCM_S1聚类算法对变化检测差异图F进行聚类，按如下步骤进行：

7a)设置聚类数目c＝2，模糊指数m＝2，阈值ε＝10^-6；

7b)随机选取两个聚类中心V₁ ⁰和V₂ ⁰，设置循环次数t＝1；

7c)更新隶属度矩阵U＝{U_ik,1≤i≤2,1≤k≤n}，其中：

$U_{\mathrm{ik}}=\frac{{({\left|\right|X_k-V_i\left|\right|}^2+\mathrm{\α}{\left|\right|{\stackrel{\‾}{X}}_k-V_i\left|\right|}^2)}^{-\frac{1}{(m-1)}}}{\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{({\left|\right|X_k-V_j\left|\right|}^2+\mathrm{\α}{\left|\right|{\stackrel{\‾}{X}}_k-V_j\left|\right|}^2)}^{-\frac{1}{(m-2)}}},$

其中，n是图像像素点的总个数，X_k是第k个像素点的灰度值,是以第k个像素点为中心的邻域均值，V_i是聚类中心，参数α用来调节原始图像和相关的均值滤波图像的权重；

7d)更新聚类中心V＝{V_i,i＝1,2}，其中：

$V_i=\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{U}_{\mathrm{ik}}^m(X_k+\mathrm{\α}{\stackrel{\‾}{X}}_k)}{(1+\mathrm{\α})\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{U}_{\mathrm{ik}}^m};$

7e)计算本次迭代生成的聚类中心与上次迭代生成的聚类中心的距离d＝||V^t-V^t-1||，如果满足d＜ε，则停止迭代，并由隶属度矩阵U得到各像素点的分类结果，否则，令t＝t+1并返回步骤7c)。