CN102750705A - 基于图像融合的光学遥感图像变化检测 - Google Patents

Abstract

本发明公开一种基于图像融合的光学遥感图像变化检测方法，主要解决现有变化检测技术检测结果精度不高的问题。其实现过程为：对同一地域不同时间获取的两幅光学遥感图像预处理后，分别构造差值法差异图和比值法差异图；对这两幅差异图分别进行N层小波分解后获取各分解层高频带和低频带的小波系数；对高频带和低频带的小波系数分别用不同的融合算子融合处理，得到高频带和低频带的融合小波系数；对高频带和低频带的融合小波系数进行逆变换得到融合后差异图；运用模糊局部C均值聚类法分割融合后的差异图，得到变化检测结果。本发明通过提高光学遥感图像的差异图性能，提高了检测结果的精度，可用于自然灾害评估和环境检测。

Description

基于图像融合的光学遥感图像变化检测

技术领域

本发明属于图像处理领域，涉及光学遥感图像变化检测，可用于环境变化中的湖泊水位动态检测，农作物生长状态的动态检测，军事侦察，自然灾害评估等诸多领域中，提高光学遥感图像变化检测结果的精确度。

背景技术

遥感图像的变化检测技术是遥感图像研究的重要组成部分，它是通过对比分析同一地域不同时刻获得的多时相遥感影像，根据图像之间的差异来得到人们所需要的地物或目标随时间发生的变化信息。变化检测技术可以检测出不同时期图像灰度值或局部纹理之间的变化，在此基础上获得感兴趣目标在形状、位置、数量及其它属性的变化情况。这些变化可能是由图像场景的真实变化引起的，也可能是由入射角、大气条件、传感器精度、地面湿度等变化引起的。光学遥感图像的变化检测是检测同一场景的两幅或多幅光学遥感图像之间随时间变化发生的变化信息。变化检测过程主要包括三部分：（1）图像预处理（2）构造差异图（3）对差异图进行分析。其中图像的预处理主要包括几何配准和辐射校正，消除图像的几何误差，以此达到同一区域不同图像的地理坐标的匹配，消除传感器自身引起的，大气辐射引起的辐射噪声。构造差异图通常采用图像差值法和比值法构造，这两种方法构造的差异图进行变化检测操作简单，对差异图进行分析是利用阈值或聚类的分割方法生成最终的变化检测结果图。

差值法是对多时相图像中对应的像素点灰度值进行相减，结果代表两个图像随时间发生的变化情况。如果差值为0，则判定该像素点的信息未发生变化，否则，判定像素点的信息发生了变化。比值法是计算已配准的多时相图像对应像素点灰度值的比值，结果代表两个图像随时间发生的变化情况。如果一个像素点没有发生变化，则比值应该接近于1，反之，则远大于1或远小于1。

尽管用差值法构造差异图和比值法构造差异图这两种变化检测方法都可以有效的检测出光学遥感图像的变化信息，但这两种方法仍然存在各自的不足：差值法容易受图像成像质量、波谱特征等客观条件影响，易产生“伪变化”信息；不能完全反映地面地物的辐射能量变化，如灰度值从50到30与250到220的变化，地面相应处的能量变化并不等量，不能等同看待；另外利用差值法构造的差异图对光学遥感图像进行变化检测时，检测结果漏检率较高。图像比值法在增强变化信息、抑制背景信息的同时可以帮助减少大气条件造成的影响，但有时会过于夸大部分的变化，如对灰度值从200到20与从20到10的变化，比值法不能区分出来，而用差值法得到的差异却非常大；此外，用比值法构造的差异影像进行光学遥感图像变化检测时，检测结果具有较高的误检率。

发明内容

本发明的目的在于针对上述差值法构造差异图和比值法构造差异图这两种变化检测方法各自的不足，提出了一种基于图像融合的光学遥感图像变化检测方法，以全面反映地面地物的辐射能量变化，减少“伪变化”信息，增强变化信息、抑制背景信息，提高检测正确率和检测精度。

实现本发明目的的技术思路是：采用基于离散小波变换的图像融合方法，对用差值法构造的差异图和用比值法构造的差异图进行融合，实现对光学遥感图像的变化检测。具体步骤包括如下：

(1)对同一地域不同时间获取的两幅光学遥感图像进行滤波去噪，辐射校正和几何配准的预处理，得到预处理后的两幅图像X₁,X₂；

(2)利用预处理后的两幅图像X₁,X₂，分别构造差值法差异图X_D和比值法差异图X_R：

X_D＝|X₁-X₂|，

X_R＝X₁/X₂；

(3)对构造出的差值法差异图X_D和比值法差异图X_R分别进行N层小波分解，得到每幅差异图在N个分解层上高频带和低频带上的多尺度分解的小波系数，N＝3；

(4)对高频带和低频带的小波系数用不同的融合算子进行融合处理，即对低频带小波系数采用取平均值的方法进行融合，得到低频带融合小波系数，对高频带小波系数采用选能量最小的方法进行融合处理，得到高频带融合小波系数；

(5)对低频带融合小波系数和高频带融合小波系数进行小波逆变换，所得的重构图像即为融合后差异图X_F；

(6)运用模糊局部C均值聚类方法对融合后差异图X_F进行图像分割，生成变化检测结果图，完成对两幅光学遥感图像变化信息的最终检测。

本发明与现有的技术相比具有如下优点:

1、本发明由于采用了基于小波分解的图像融合技术，在时间和空间上分解频率，可以更容易的提取细节信息，可以很好的保存图像的细节信息，从而得到正确率更高的变化检测结果。

2、本发明由于利用小波分解的方向性，针对人眼对不同方向的高频分量具有不同的分辨率这一视觉特性，可以获得视觉效果更佳的融合图像。

3、本发明由于选择不同的合适的融合规则融合图像，可以更大程度上的抑制变化区域的背景信息，增强变化区域的变化信息，使差值法的漏检率和比值法的误检率互补，获得了较高的检测精度。

附图说明

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

图2是本发明中基于小波变换的图像融合子流程图；

图3是本发明仿真使用的Mexico地区两幅光学遥感图像和一幅标准参考图；

图4是本发明分别使用差值运算、比值运算和图像融合法构造的Mexico地区三幅光学遥感差异图；

图5是用本发明和现有的差值法、比值法对Mexico地区光学遥感图像进行变化检测的仿真结果图；

图6是本发明仿真使用的意大利Sardinia地区两幅光学遥感图像和一幅标准参考图；

图7是本发明分别使用差值运算、比值运算和图像融合法构造的意大利Sardinia地区三幅光学遥感图像差异图；

图8是用本发明和现有的差值法、比值法对意大利Sardinia地区光学遥感图像进行变化检测的仿真结果图。

具体实施方式

以下参照附图，对本发明的实现方案及优点进行详细描述。

实施例1,

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

步骤1，在Mexico同一地域的不同时间获取两幅光学遥感图像，将这两幅光学遥感图像进行滤波去噪，辐射校正与几何配准的预处理，得到处理后的两幅图像X₁,X₂，其中预处理后得到的图像X₁如图3(a)所示，预处理后得到的图像X₂如图3(b)所示。

通过预处理可以消除图像的几何误差，以达到对同一区域不同图像的地理坐标的匹配，消除传感器自身引起的噪声和大气辐射引起的辐射噪声。

步骤2，用预处理后的两幅图像X₁,X₂，分别构造差值法差异图X_D和比值法差异图X_R。

(2.1)将如图3(a)所示的图像X₁中位于i行j列的像素点(i,j)的灰度值X₁(i,j)和相对应的图3(b)所示的图像X₂中位于i行j列的像素点的灰度值X₂(i,j)，通过差值运算X_D(i,j)＝|X₁(i,j)-X₂(i,j)|得到差值法差异图X_D中位于i行j列的像素点(i,j)的灰度值X_D(i,j)，若X_D(i,j)为0，则表示图像X₁中该像素点没有随时间的变化而发生变化，否则，认为图像X₁中的该像素点发生了变化；对图像X₁和图像X₂中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行差值运算，得到如图4(a)所示的差值法差异图X_D；

(2.2)将如图3(a)所示的图像X₁中位于i行j列的像素点(i,j)的灰度值X₁(i,j)和相对应的图3(b)所示的图像X₂中位于i行j列的像素点的灰度值X₂(i,j)，通过比值运算X_R(i,j)＝X₁(i,j)/X₂(i,j)，得到比值法差异图X_R中位于i行j列的像素点(i,j)的灰度值X_R(i,j)，如果比值接近于1,表示图像X₁中位于i行j列的像素点没有随时间的变化而发生变化，反之，比值若远大于或远小于1，则认为该像素点发生了变化；对图像X₁和图像X₂中每个位于i行j列的像素点的灰度值从左到右，从上到下都进行比值运算，构造如图4(b)所示的比值法差异图X_R。

步骤3，分别对差值法差异图X_D和比值法差异图X_R进行N层小波分解，得到每幅差异图在N个分解层上高频带和低频带上的多尺度分解的小波系数，N=3。

(3.1)在第一层，即N＝1，将如图4(a)所示的差值法差异图X_D分解为一个低频带子图像和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数

将如图4(b)所示的比值法差异图X_R分解为一个低频带子图像

和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数

其中LL表示低频带，ε＝LH,HL,HH，LH表示垂直方向的高频带，HL表示水平方向的高频带，HH表示对角线方向的高频带；

(3.2)在第二层，即N=2，将如图4(a)所示的差值法差异图X_D经步骤(3a)分解所得的低频带子图像

继续分解为一个低频带子图像

和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数

将如图4(b)所示的比值法差异图X_R经步骤(3a)分解所得的低频带子图像

继续分解为一个低频带子图像和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数

(3.3)在第三层，即N=3，将如图4(a)所示的差值法差异图X_D经步骤(3b)分解所得的低频带子图像

继续分解为一个低频带子图像

和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数

将如图4(b)所示的比值法差异图X_R经步骤(3b)分解所得的低频带子图像

继续分解为一个低频带子图像

和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数

步骤4，对高频带和低频带的小波系数用不同的融合算子进行融合处理，得到高频带、低频带融合小波系数。

(4.1)对低频带小波系数采用取平均值的方法进行融合，得到低频带融合小波系数通过如下公式进行：

$W_{\mathrm{LL}}^F(i,j)=\frac{W_{\mathrm{LL},3}^D(i,j)+W_{\mathrm{LL},3}^R(i,j)}{2}$

其中，D指代差值法差异图的子图像，R指代比值法差异图的子图像，F指代融合后差异图的子图像，(i,j)表示图像中像素点的位置为i行j列，LL表示低频带，

指融合后差异图像的低频带小波系数

在像素点(i,j)的值，

指差值法差异图在第三分解层的低频带小波系数

在像素点(i,j)的值，

指比值法差异图在第三分解层的低频带小波系数

在像素点(i,j)的值。

由于分解后的低频带主要代表了源图像的轮廓特征，可以明确的反应两个不同时刻的源图像之间的变化信息，所以对低频带小波系数采用取平均值的方法进行融合，可以最大程度的增强变化区域的边缘特征。

(4.2)对高频带小波系数采用选能量最小的方法进行融合处理，得到高频带融合小波系数：

$W_{\mathrm{\&epsiv;},N}^F(i,j)=\left\{\begin{array}{cc}{W}_{\mathrm{\&epsiv;},N}^D(i,j),& E_{\mathrm{\&epsiv;},N}^D(i,j)<E_{\mathrm{\&epsiv;},N}^R(i,j)\\ W_{\mathrm{\&epsiv;},N}^R(i,j),& E_{\mathrm{\&epsiv;},N}^D(i,j)\&GreaterEqual;E_{\mathrm{\&epsiv;},N}^R(i,j)\end{array}\right.,$

其中，N＝1,2,3，D指差值法差异图的子图像，R指比值法差异图的子图像，F指融合后差异图的子图像，(i,j)表示图像中像素点的位置为i行j列，为融合后差异图像中第N分解层的高频带小波系数在像素点(i,j)的值，

为融合后差异图像中第N分解层的高频带小波系数

在像素点(i,j)的值，

为融合后差异图像中第N分解层的高频带小波系数

在像素点(i,j)的值，ε＝LH,HL,HH，LH表示垂直方向的高频带，HL表示水平方向的高频带，HH表示对角线方向的高频带；

表示差值法差异图X_D第N分解层的局部能量系数，其中

表示以(i,j)为中心的局部窗口

中第k个小波系数值；

表示比值法差异图X_R第N分解层的局部能量系数，其中

表示以(i,j)为中心的局部窗口

中第k个小波系数值。

高频带主要反应源图像的细节信息，比如具有跳跃性的边缘和线条特征，所以，本发明选择局部区域能量最小的小波系数作为融合后的小波系数，以达到最大程度的抑制背景信息的目的。这种融合方法可以将如图4(a)所示的差值法差异图和如图4(b)所示的比值法差异图的高频带中的均匀区域能很好的融合在一起，减小比值运算构造差异图的变化检测方法对变化区域造成的夸大效应。

步骤5，对低频带融合小波系数和高频带融合小波系数进行小波逆变换，所得的重构图像即为如图4(c)所示的融合后差异图X_F。

上述步骤3-步骤5的内容，如图2所示。

步骤6，运用模糊局部C均值聚类即FLICM方法对融合后差异图X_F进行分割，获得如图5(c)所示的变化检测结果图，完成对两幅光学遥感图像变化信息的最终检测。

实施例2

参照图1和图2，本实例的实现步骤如下：

步骤一，对同一地域即意大利Sardinia地区不同时间获取的两幅光学遥感图像进行滤波去噪，辐射校正与几何配准的预处理，得到处理后的两幅图像X′₁,X′₂，其中预处理后得到的图像X′₁如图6(a)所示，预处理后得到的图像X′₂如图6(b)所示。

步骤二，利用预处理后的如图6(a)所示的图像X′₁和图6(b)所示的图像X′₂，用差值运算构造如图7(a)所示的差值法差异图X′_D，用比值运算构造如图7(b)所示的比值法差异图X′_R，具体实施步骤如实施例1中的步骤2。

步骤三，对如图7(a)所示的差值法差异图X′_D和如图7(b)所示的比值法差异图X′_R分别进行N层小波分解，得到每幅差异图在N个分解层上高频带和低频带上的多尺度分解的小波系数，N＝3，具体实施步骤如实施例1中的步骤3。

步骤四，对高频带和低频带的小波系数用不同的融合算子进行融合处理，即对低频带小波系数采用取平均值的方法进行融合，得到低频带融合小波系数，对高频带小波系数采用选能量最小的方法进行融合处理，得到高频带融合小波系数，具体实施步骤如实施例1中的步骤4。

步骤五，对低频带融合小波系数和高频带融合小波系数进行小波逆变换，所得的重构图像即为如图7(c)所示的融合后差异图X′_F。

上述步骤三-步骤五的内容，如图2所示。

步骤六，运用模糊局部C均值聚类即FLICM的方法对融合后的差异图X′_F进行分割，获得如图8(c)所示的变化检测结果图，完成对两幅光学遥感图像变化信息的最终检测。

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

1、仿真条件：

本发明的仿真是在主频1.87GHZ的Intel Pentium CPU P6000、内存2.00GB的硬件环境和MATLAB R2009a的软件环境下进行的。

2、仿真参数

对于具有参考图的实验仿真图，可进行定量的变化检测结果分析：

①漏检数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检数；

②误检数：统计实验结果图中未发生变化区域的像素个数，与参考图中未发生变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为误检个数；

③总检测错误数：漏检数和误检数的和；

④

3、仿真实验内容与结果分析

将本发明中图像融合构造的差异图，与差值运算、比值运算构造的差异图，用模糊局部C均值聚类即FLICM方法进行分割得到变化检测结果图对比和定量分析，通过两组真实遥感数据集的仿真实验完成。

实验1.Mexico地区遥感图像数据

这里将本发明中涉及的变化检测方法应用在如图3所示的Mexico地区光学遥感图像上，并与现有的差值法变化检测方法和比值法变化检测方法进行对比分析。

如图3所示的Mexico地区三幅图像，大小均为512×512像素，其中图3(a)表示2000年4月Mexico地区的地貌信息；图3(b)表示2002年5月Mexico地区的地貌信息。图3(c)为变化检测的标准参考图。

用图3(a)和图3(b)所示的两幅光学遥感图像分别使用差值运算、比值运算和图像融合法构造三幅差异图，如图4所示，其中图4(a)为差值法差异图；图4(b)为比值法差异图；图4(c)为融合后差异图。

将图4所示的三幅差异图分别用FLICM聚类方法分割得到变化检测结果图如图5所示，其中图5(a)为差值法的变化检测结果图；图5(b)为比值法的变化检测结果图；图5(c)为本发明的变化检测结果图。

从图5可以看出，差值法的变化检测结果图存在部分像素点的漏检问题，比值法的变化检测结果图则有较多的斑点，同时也存在误检问题，而本发明的变化检测结果图则轮廓清晰，比较接近标准参考图中的实际变化区域轮廓。

将图5所示的三种方法的性能指标进行对比，如表1所示：

表1    Mexico地区三种方法变化检测结果性能比较

从表1中可以看出，本发明方法的变化检测结果误检数和漏检数比较均衡，处于中间水平，与差值法和比值法的变化检测结果相比，总错误数最低，正确率也有了明显提高，显示了本发明方法的优越性。

实验2.意大利Sardinia地区的真实光学遥感图像

将本发明中涉及的变化检测方法应用在如图6所示的意大利Sardinia地区光学遥感图像上，并与现有的差值法变化检测方法和比值法变化检测方法进行对比分析。

图6中的三幅图像，大小均为412×300像素，其中图6(a)表示1995年9月意大利Sardinia地区的地貌信息；图6(b)表示1996年7月意大利Sardinia地区的地貌信息；图6(c)为变化检测的标准参考图。

用图6(a)和图6(b)所示的两幅光学遥感图像分别使用差值运算、比值运算和图像融合法构造三幅差异图，如图7所示，其中图7(a)为差值法差异图；图7(b)为比值法差异图；图7(c)为融合后差异图。

将图7所示的三幅差异图分别用FLICM聚类方法分割得到变化检测结果图，如图8所示。其中图8(a)为差值法的变化检测结果图；图8(b)为比值法的变化检测结果图；图8(c)为本发明方法变化检测结果图。

从参考图8可以看出，差值法变化检测结果图中许多小斑点，而本发明方法的变化检测结果图轮廓很清晰，检测出的变化区域轮廓也最接近标准参考图中的实际变化区域轮廓。

将图8所示的三种方法的性能指标进行对比，如表2所示：

表2    意大利Sardinia地区变化检测结果性能比较

从表2可以看出，差值法的变化检测结果图的漏检数高达1038，比值法的变化检测结果图的误检数高达3102，而本发明方法的变化检测结果图的漏检数降为586，误检数为1370，从而降低了变化检测结果的总错误数，使正确率提升到98.42%。

Claims (4)

1.一种基于图像融合的光学遥感图像变化检测方法，包括如下步骤：

(1)对同一地域不同时间获取的两幅光学遥感图像进行滤波去噪，辐射校正和几何配准的预处理，得到预处理后的两幅图像X₁,X₂；

(2)利用预处理后的两幅图像X₁,X₂，分别构造差值法差异图X_D和比值法差异图X_R：

X_D＝|X₁-X₂|，

X_R＝X₁/X₂；

(3)对构造出的差值法差异图X_D和比值法差异图X_R分别进行N层小波分解，得到每幅差异图在N个分解层上高频带和低频带上的多尺度分解的小波系数，N＝3；

(4)对高频带和低频带的小波系数用不同的融合算子进行融合处理，即对低频带小波系数采用取平均值的方法进行融合，得到低频带融合小波系数，对高频带小波系数采用选能量最小的方法进行融合，得到高频带融合小波系数；

(5)对低频带融合小波系数和高频带融合小波系数进行小波逆变换，所得的重构图像即为融合后差异图X_F；

(6)运用模糊局部C均值聚类方法对融合后差异图X_F进行图像分割，生成变化检测结果图，完成对两幅光学遥感图像变化信息的最终检测。

2.根据权利要求1所述的基于图像融合的光学遥感图像变化检测方法，其中步骤（3）中所述的对构造出的差值法差异图X_D和比值法差异图X_R分别进行N层小波分解，按如下步骤进行：

(3a)在第一层，即N＝1，将差值法差异图X_D分解为一个低频带子图像和三个高频带子图像

得到一个低频带小波系数和三个高频带小波系数将比值法差异图X_R分解为一个低频带子图像

和三个高频带子图像得到一个低频带小波系数

和三个高频带小波系数其中LL表示低频带，ε＝LH,HL,HH，LH表示垂直方向的高频带，HL表示水平方向的高频带，HH表示对角线方向的高频带；

(3b)在第二层，即N＝2，将差值法差异图X_D经步骤(3a)分解所得的低频带子图像

继续分解为一个低频带子图像

和三个高频带子图像

得到一个低频带小波系数

和三个高频带小波系数将比值法差异图X_R经步骤(3a)分解所得的低频带子图像

继续分解为一个低频带子图像

和三个高频带子图像得到一个低频带小波系数

和三个高频带小波系数

(3c)在第三层，即N＝3，将差值法差异图X_D经步骤(3b)分解所得的低频带子图像

继续分解为一个低频带子图像

和三个高频带子图像得到一个低频带小波系数

和三个高频带小波系数

将比值法差异图X_R经步骤(3b)分解所得的低频带子图像继续分解为一个低频带子图像和三个高频带子图像

得到一个低频带小波系数和三个高频带小波系数

3.根据权利要求1所述的基于图像融合的光学遥感图像变化检测方法，其中步骤（4）中所述的对低频带小波系数采用取平均值的方法进行融合，通过如下公式进行：

$W_{\mathrm{LL}}^F(i,j)=\frac{W_{\mathrm{LL},3}^D(i,j)+W_{\mathrm{LL},3}^R(i,j)}{2}$

其中，D指代差值法差异图的子图像，R指代比值法差异图的子图像，F指代融合后差异图的子图像，(i,j)表示图像中像素点的位置为i行j列，LL表示低频带，

指融合后差异图像的低频带小波系数

在像素点(i,j)的值，

指差值法差异图在第三分解层的低频带小波系数在像素点(i,j)的值，

指比值法差异图在第三分解层的低频带小波系数

在像素点(i,j)的值。

4.根据权利要求1所述的基于图像融合的光学遥感图像变化检测方法，其中步骤（4）中所述的对高频带小波系数采用选能量最小的方法进行融合，是选择能量值较小者的小波系数作为融合后的小波系数，通过如下公式进行：

$W_{\mathrm{\&epsiv;},N}^F(i,j)=\left\{\begin{array}{cc}{W}_{\mathrm{\&epsiv;},N}^D(i,j),& E_{\mathrm{\&epsiv;},N}^D(i,j)<E_{\mathrm{\&epsiv;},N}^R(i,j)\\ W_{\mathrm{\&epsiv;},N}^R(i,j),& E_{\mathrm{\&epsiv;},N}^D(i,j)\&GreaterEqual;E_{\mathrm{\&epsiv;},N}^R(i,j)\end{array}\right.$

其中，N＝1,2,3，D指代差值法差异图的子图像，R指代比值法差异图的子图像，F指代融合后差异图的子图像，(i,j)表示图像中像素点的位置为i行j列，

表示融合后差异图像中第N分解层的高频带小波系数

在像素点(i,j)的值，

表示融合后差异图像中第N分解层的高频带小波系数

在像素点(i,j)的值，

表示融合后差异图像中第N分解层的高频带小波系数

在像素点(i,j)的值，ε＝LH,HL,HH，LH表示垂直方向的高频带，HL表示水平方向的高频带，HH表示对角线方向的高频带；

表示差值法差异图X_D第N分解层的局部能量系数，其中

表示以(i,j)为中心的局部窗口

中第k个小波系数值；

表示比值法差异图X_R第N分解层的局部能量系数，其中表示以(i,j)为中心的局部窗口

中第k个小波系数值。

Images