CN106855947A - 基于核互模态因素分析核融合的多光谱图像变化检测方法 - Google Patents

Abstract

本发明公开了基于核互模态因素分析核融合的多光谱图像变化检测方法，主要解决现有差值核函数不能充分挖掘及利用图像多种特征之间的相关性而导致检测精度低下的问题。其实现过程是：1)将两时刻三通道光谱图像变换至HIS空间，得到两时刻图像的色调、饱和度、强度信息；2)根据1)的结果提取两时刻图像颜色特征和强度纹理特征；3)利用所提特征求取颜色差值核函数矩阵和强度纹理差值核函数矩阵；4)对3)获得的两个矩阵加权融合，构造合成核函数矩阵，并自适应选取最优的加权系数；5)将合成核函数矩阵输入支持向量机SVM进行检测，得到变化检测结果。本发明检测精度高，结果稳定，运算量低，可用于多光谱图像变化检测。

Description

基于核互模态因素分析核融合的多光谱图像变化检测方法

技术领域

本发明属于图像处理技术领域，特别涉及一种图像变化检测方法，可用于对多光谱图像进行地物状态变化监测与评估。

背景技术

多光谱图像是遥感卫星获取的在多个波段上的对地观测数字图像。多光谱图像变化检测通过分析同一地区不同时刻观测到的多幅多光谱图像来确定地物变化的特征及过程。多光谱图像具有丰富的光谱信息，能够很好地反映地物特性。多光谱图像变化检测在军事和民事领域中具有突出的作用。

最常用的多光谱图像变化检测方法是变化矢量分析法CVA，该方法将多通道的原始两时刻图像定义为向量，将两时刻的向量作差得到变化向量，再对变化向量进行分析从而检测出两时刻图像是否发生变化。CVA方法在多光谱图像变化检测中得到了成功的应用，近年各国学者对其进行了多方面的改进。

J Chen定义了后验概率空间中的CVA算法，该算法降低了传统CVA对几何配准和辐射校正的严格要求，不同变化所得到的结果都处于同一尺度空间，对单门限检测法的鲁棒性更好。

F Bovolo提出了极坐标空间中的CVA算法，该算法从理论上分析了极坐标空间中变化与非变化像素点的分布及其统计特性，给出了极坐标空间中变化矢量的严格定义及分析。此后，F Bovolo又提出了2维CVA和多级CVA算法。但是上述的CVA方法及其改进方法，不能够实现对多光谱图像多种特征的自适应融合，不能充分利用多种特征中的信息，导致运算量比较大，计算效率比较低。

除了CVA方法，其它方法也在多光谱图像变化检测中得到了成功的应用。例如，Hichri等在2013年提出了一种实现交式互分割的多光谱图像变化检测方法，该方法使用支持向量机SVM对已知标记的像素进行训练，使用训练结果指导图像分割得到初始变化检测图像，最后使用马尔科夫随机场MRFs调整初始变化检测结果图，得到最终的变化检测结果图。但是该方法实现时需要两个步骤，即先对已知标签样本进行训练得到初分结果，再进行优化，使得该算法流程比较复杂。

近年来新发展起来的是基于核方法的多光谱图像变化检测算法。GustavoCamps-Valls在2008年首先提出将核方法应用于遥感图像变化检测，通过提取图像的强度信息和纹理信息并构造强度纹理差值核函数实现了遥感图像变化检测。但是该方法处理多光谱图像时，没有考虑多种特性之间的相关性，不能充分挖掘和融合多种图像特征，因而其检测精度有待进一步提高。

发明内容

本发明的目的是在针对上述已有算法的不足，提出一种基于核互模态因素分析核融合的多光谱图像变化检测方法，通过分析图像多种特征之间的相关性，充分挖掘和融合图像特征，从而降低运算量，提高运算效率，提高检测精度。

为实现上述目的，本发明的实现步骤包括如下：

(1)给定原始的预处理后的三通道的相同尺寸两时刻光谱图像，记为{R_t,G_t,B_t}，t＝t₀,t₁，其中，R_t,G_t,B_t分别表示原RGB空间中的三通道光谱图像，使用RGB-HIS变换将该三通道图像变换至HIS空间，得到HIS通道的两时刻光谱图像{H_t,I_t,S_t}，其中H_t为t时刻光谱图像的色调分量，S_t为t时刻光谱图像的饱和度分量，I_t为t时刻光谱图像的强度分量；

(2)利用t₀,t₁两时刻光谱图像的色调饱和度和强度分量，提取光谱图像特征；

2a)分别提取分量光谱图像的灰度值和分量光谱图像的灰度值，将其堆砌成两时刻光谱图像的颜色特征和

2b)提取强度分量光谱图像的灰度特征和纹理特征将其堆砌成强度纹理特征

(3)利用两时刻光谱图像的颜色特征求取颜色差值核函数矩阵K_Col；

(4)利用两时刻光谱图像的强度纹理特征求取强度纹理差值核函数矩阵K_ST；

(5)利用基于核互模态因素分析算法KCFA分别求取颜色差值核函数矩阵K_Col的权重系数w_c和强度纹理差值核函数矩阵K_ST的权重系数w_s；

(6)利用步骤5)、步骤4)和步骤3)的结果，构造合成核函数矩阵K_F：K_F＝w_c·K_Col+w_s·K_ST；

(7)将步骤6)得到的合成核函数矩阵K_F输入到支持向量机SVM中进行检测，得到二值变化检测结果，即完成对多光谱图像的变化检测。

本发明具有如下优点：

1)本发明由于利用了核方法，可以将实际多光谱图像处理中存在的非线性不可分问题转化为高维空间中的线性可分问题。

2)本发明通过将原RGB空间的伪彩色图像变换至HIS空间，充分提取图像的颜色特征和强度纹理特征，有效突出了变化区域及其边界；并将颜色信息和强度纹理信息进行融合，最大限度挖掘两种特征之间的相关性，有效实现了核空间中的特征融合。

3)本发明能全面地提取图像特征，可以显著提高变化检测结果的精度。

仿真结果表明，本发明与现有的差值核方法相比，运算量低，检测精度高，检测性能更为优越。

附图说明

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

图2是用本发明对日本地震灾害实测多光谱图像的变化检测结果图；

图3是用本发明对格洛斯特洪水灾害实测多光谱图像的变化检测结果图；

图4是用本发明对哈萨克斯坦海岸受流沙冲击前后实测多光谱图像的变化检测结果图；

图5是用本发明对塔霍湖火灾实测多光谱图像的变化检测结果图。

具体实施方式

参照图1，本发明的具体实施过程如下：

步骤1.将预处理后的光谱图像变换至HIS空间。

给定原始的预处理后的三通道的相同尺寸两时刻光谱图像，记为{R_t,G_t,B_t}，t＝t₀,t₁，其中，R_t,G_t,B_t分别表示原RGB空间中的三通道光谱图像；

使用RGB-HIS变换将该三通道光谱图像变换至HIS空间，得到HIS通道的两时刻光谱图像{H_t,I_t,S_t}，其中H_t为t时刻光谱图像的色调分量，S_t为t时刻光谱图像的饱和度分量，I_t为t时刻光谱图像的强度分量。

其中，RGB空间中的R代表光谱图像的红色分量，G代表光谱图像的绿色分量，B代表光谱图像的蓝色分量；HIS空间中的H代表光谱图像的色调分量，S代表光谱图像的饱和度分量，I光谱图像的代表强度分量。

步骤2.提取光谱图像的颜色特征和强度纹理特征。

利用t₀,t₁两时刻光谱图像的色调饱和度和强度分量，提取光谱图像特征，其步骤如下：

2a)提取t₀时刻光谱图像分量的灰度值和分量的灰度值用这两个灰度值堆砌构成t₀时刻光谱图像的颜色特征提取t₁两时刻光谱图像分量的灰度值和分量的灰度值用这两个灰度值堆砌构成t₁时刻光谱图像的颜色特征表示如下：

2b)提取t₀时刻光谱图像的强度分量的灰度值和纹理特征用该灰度值和纹理特征堆砌构成t₀时刻光谱图像强度纹理特征提取t₁时刻光谱图像分量灰度值和纹理特征用该灰度值和纹理特征堆砌构成t₁时刻光谱图像的强度纹理特征表示如下：

步骤3.构造颜色差值核函数矩阵K_Col。

利用t₀时刻光谱图像的颜色特征和t₁时刻光谱图像的颜色特征构造颜色差值核函数矩阵K_Col，其步骤如下：

3a)从t₀时刻光谱图像中任意选取两个像素点和从t₁时刻光谱图像中选取与和位置相同的两个像素点，记为和

3b)从中找出像素点对应的颜色特征和像素点对应的颜色特征从中找出像素点对应的颜色特征和像素点对应的颜色特征

3c)利用3b)获得的颜色特征求取颜色插值核函数的值：

其中，K(·)表示高斯径向基RBF核函数，ζ代表核宽参数，exp为指数函数；

3d)在t₀时刻光谱图像中找出没有被选取过的任意两个像素以及中相同位置的像素点按照与3b)相同的方法找到与对应的颜色特征再按照与3c)相同的方法计算颜色差值核函数的值，

其中p是像素点的标号，1≤p≤N-1，N是光谱图像像素点个数；

3e)重复步骤3d)，直到所有像素点被选取过，得到颜色差值核函数的全部值，构成颜色差值核函数矩阵K_Col。

步骤4.构造强度纹理差值核函数矩阵K_ST。

利用t₀时刻光谱图像强度纹理特征和t₁时刻光谱图像强度纹理特征构造强度纹理差值核函数矩阵K_ST，其步骤如下：

4a)从t₀时刻光谱图像中任意选取两个像素点和从t₁时刻光谱图像中选取与和位置相同的两个像素点，记为和

4b)从中找出像素点对应的强度纹理特征和像素点对应的强度纹理特征从中找出像素点对应的强度纹理特征和像素点对应的强度纹理特征

4c)利用从4b)获得的求取强度纹理差值核函数的值：

其中，K(·)表示高斯径向基RBF核函数，ζ代表核宽参数，exp为指数函数；

4d)在t₀时刻光谱图像中找出没有被选取过的任意两个像素以及中相同位置的像素点按照与4b)相同的方法找到与对应的强度纹理特征再按照与4c)相同的方法计算强度纹理差值核函数的值，

其中q是像素点的标号，1≤q≤N-1，N是光谱图像像素点个数；

4e)重复步骤4d)，直到所有像素点被选取过，得到强度纹理差值核函数的全部值，构成强度纹理差值核函数矩阵K_ST。

步骤5.求取权重系数w_c和w_s。

利用步骤3得到的颜色差值核函数矩阵K_Col和步骤4得到的强度纹理差值核函数矩阵K_ST，利用基于核互模态因素分析算法KCFA分别求取颜色差值核函数矩阵K_Col的权重系数w_c和强度纹理差值核函数矩阵K_ST的权重系数w_s，其具体步骤如下：

5a)通过特征值分解法得到由K_ST*K_Col特征向量构成的第一投影矩阵V和特征值构成的第一对角矩阵λ，这三者的关系如下：

(K_ST*K_Col)*V＝λ*V；

5b)提取第一投影矩阵V的第一列，记为α，定义颜色差值核函数矩阵K_Col融合系数w_c：

其中，T表示矩阵转置，K_Col是颜色差值核函数矩阵；

5c)通过特征值分解法得到由K_Col*K_ST特征向量构成的第二投影矩阵U和特征值构成的第二对角矩阵λ′，这三者的关系如下：

(K_Col*K_ST)*U＝λ′*U；

5d)提取第二投影矩阵U的第一列，记为β，定义强度纹理插值核函数矩阵K_ST的融合系数w_s：

其中，K_ST是强度纹理差值核函数矩阵。

步骤6.构造合成核函数矩阵K_F。

利用步骤5得到的权重系数w_c和w_s，以及步骤3得到的颜色差值核函数矩阵K_Col和步骤4得到的强度纹理差值核函数矩阵K_ST构造合成核函数矩阵：

K_F＝w_c·K_Col+w_s·K_ST。

步骤7.完成光谱图像的变化检测。

将步骤6得到的合成核函数矩阵K_F输入到支持向量机SVM中进行检测，得到二值变化检测结果，即完成对光谱图像的变化检测。

本发明效果可以通过以下实验进一步证实：

1.实验条件

实验仿真环境为：MATLAB R2009b，Intel(R)Core(TM)I5-4590 CPU@3.3GHz，Window 7旗舰版。

2.实验内容

实验1，应用差值核方法和本发明对日本地震灾害真实多光谱图像进行变化检测，结果如图2，其中，图2(a)是实测日本地震灾害的第一时刻Landsat5TM图像，图2(b)是实测日本地震灾害的第二时刻Landsat5TM图像，图2(c)是变化检测结果参考图，图2(d)是差值核方法的变化检测结果，图2(e)是本发明的变化检测结果；

实验2，应用差值核方法和本发明对格洛斯特洪水灾害实测多光谱图像进行变化检测，结果如图3，其中，图3(a)是实测格洛斯特洪水灾害的第一时刻SPOT-XS图像，图3(b)是实测格洛斯特洪水灾害的第二时刻SPOT-XS图像，图3(c)是变化检测结果参考图，图3(d)是差值核方法的变化检测结果，图3(e)是本发明的变化检测结果；

实验3，应用差值核方法和本发明对哈萨克斯坦海岸受流沙冲击前后实测多光谱像进行变化检测，结果如图4，其中，图4(a)是实测哈萨克斯坦海岸受流沙冲击前的第一时刻Landsat7ETM+图像，图4(b)是实测城市受洪水哈萨克斯坦海岸受流沙冲击后的第二时刻Landsat7ETM+图像，图4(c)是变化检测结果参考图，图4(d)是差值核方法的变化检测结果，图4(e)是本发明的变化检测结果；

实验4，应用差值核方法和本发明对塔霍湖火灾实测多光谱图像进行变化检测，结果如图5，其中，图5(a)是实测塔霍湖火灾的第一时刻Landsat5TM图像，图5(b)是实测塔霍湖火灾的第二时刻Landsat5TM图像，图5(c)是变化检测结果参考图，图5(d)是差值核方法的变化检测结果，图5(e)是本发明的变化检测结果。

将上述四个实验的结果进行比较，如表1所示，选取检测精度、Kappa系数作为指标参数，评价本发明检测精度，验证本发明检性能。

表1 本发明方法与差值核方法检测结果精度比较

将四个实验中本发明方法与差值核方法的运行时间进行比较，如表2所示，用以定量的说明本发明与差值核方法变化检测相比的速度优势。

表2 本发明方法与差值核方法运行时间对比

3.实验结果分析

从图2、图3、图4，图5和表1、表2可以看到，与基于差值核的变化检测方法相比，本发明的检测性能更为优越，在变化区域具有更好的一致性，且能突出变化边界，检测精度高，运算量低，检测速度快。这是因为本发明充分提取了图像的颜色特征和强度纹理特征，有效突出了变化区域及其边界，并最大限度的挖掘颜色特征和强度纹理两种特征之间的相关性，降低了运算量，有效实现了核空间中的特征融合。

Claims (7)

1.一种基于核互模态因素分析核融合的多光谱图像变化检测方法，包括如下步骤：

1)给定原始的预处理后的三通道的相同尺寸两时刻光谱图像，记为{R_t,G_t,B_t}，t＝t₀,t₁，其中，R_t,G_t,B_t分别表示原RGB空间中的三通道光谱图像，使用RGB-HIS变换将该三通道图像变换至HIS空间，得到HIS通道的两时刻光谱图像{H_t,I_t,S_t}，其中H_t为t时刻光谱图像的色调分量，S_t为t时刻光谱图像的饱和度分量，I_t为t时刻光谱图像的强度分量；

2)利用t₀,t₁两时刻光谱图像的色调饱和度和强度分量，提取光谱图像特征；

2a)分别提取分量光谱图像的灰度值和分量光谱图像的灰度值，将其堆砌成两时刻光谱图像的颜色特征和

2b)提取强度分量光谱图像的灰度特征和纹理特征将其堆砌成强度纹理特征

3)利用两时刻光谱图像的颜色特征求取颜色差值核函数矩阵K_Col；

4)利用两时刻光谱图像的强度纹理特征求取强度纹理差值核函数矩阵K_ST；

5)利用基于核互模态因素分析算法KCFA分别求取颜色差值核函数矩阵K_Col的权重系数w_c和强度纹理差值核函数矩阵K_ST的权重系数w_s；

6)利用步骤5)、步骤4)和步骤3)的结果，构造合成核函数矩阵K_F：K_F＝w_c·K_Col+w_s·K_ST；

7)将步骤6)得到的合成核函数矩阵K_F输入到支持向量机SVM中进行检测，得到二值变化检测结果，即完成对多光谱图像的变化检测。

2.根据权利要求1所述的方法，其中步骤2a)中堆砌t₀时刻颜色特征和t₁时刻颜色特征表示如下：

${\mathrm{Col}}_{t_0}=\{{\mathrm{hH}}_{t_0},{\mathrm{hS}}_{t_0}\}$

${\mathrm{Col}}_{t_1}=\{{\mathrm{hH}}_{t_1},{\mathrm{hS}}_{t_1}\}$

其中，和分别为提取出的t₀时刻图像色调分量的灰度值和光谱图像饱和度分量的灰度值；和分别为提取出的t₁时刻光谱图像色调分量的灰度值和光谱图像饱和度分量的灰度值。

3.根据权利要求1所述的方法，其中步骤2b)中堆砌t₀时刻强度纹理特征和t₁时刻强度纹理特征表示如下：

${\mathrm{ST}}_{t_0}=\{{\mathrm{hI}}_{t_0},{\mathrm{\ωI}}_{t_0}\}$

${\mathrm{ST}}_{t_1}=\{{\mathrm{hI}}_{t_1},{\mathrm{\ωI}}_{t_1}\}$

其中和分别表示t₀时刻光谱图像强度分量的灰度值和t₁时刻光谱图像强度分量的灰度值；和分别表示t₀时刻光谱图像的强度分量的纹理特征和t₁时刻光谱图像的强度分量的纹理特征。

4.根据权利要求1所述方法，其中步骤3)中利用两时刻光谱图像的颜色特征构造颜色差值核函数矩阵K_Col，按如下步骤进行：

3a)从t₀时刻光谱图像中任意选取两个像素点和从t₁时刻光谱图像中选取与和位置相同的两个像素点，记为和

3b)从中找出像素点对应的颜色特征和像素点对应的颜色特征从中找出像素点对应的颜色特征和像素点对应的颜色特征

3c)利用从3b)中获得的颜色特征求取颜色插值核函数的值：

其中，K(·)表示高斯径向基RBF核函数，ζ代表核宽参数，exp为指数函数；

3d)在t₀时刻光谱图像中找出没有被选取过的任意两个像素以及中相同位置的像素点按照与3b)相同的方法找到与对应的颜色特征再按照与3c)相同的方法计算颜色差值核函数的值，其中p是像素点的标号，1≤p≤N-1，N是光谱图像像素点个数；

3e)重复步骤3d)，直到所有像素点被选取过，得到颜色差值核函数的全部值，构成颜色差值核函数矩阵K_Col。

5.根据权利要求1所述方法，其中步骤4)中利用两时刻光谱图像强度纹理特征构造强度纹理差值核函数矩阵K_ST，按如下步骤进行：

4a)从t₀时刻光谱图像中任意选取两个像素点和从t₁时刻光谱图像中选取与和位置相同的两个像素点，记为和

4b)从中找出像素点对应的强度纹理特征和像素点对应的强度纹理特征从中找出像素点对应的强度纹理特征和像素点对应的强度纹理特征

4c)利用从4b)中获得的求取强度纹理差值核函数的值：

其中，K(·)表示高斯径向基RBF核函数，ζ代表核宽参数，exp为指数函数；

4d)在t₀时刻光谱图像中找出没有被选取过的任意两个像素以及中相同位置的像素点按照与4b)相同的方法找到与对应的强度纹理特征再按照与4c)相同的方法计算强度纹理差值核函数的值，其中q是像素点的标号，1≤q≤N-1，N是光谱图像像素点个数；

4e)重复步骤4d)，直到所有像素点被选取过，得到强度纹理差值核函数的全部值，构成强度纹理差值核函数矩阵K_ST。

6.根据权利要求1所述方法，其中步骤5)中求取颜色差值核函数矩阵K_Col的权重系数w_c和强度纹理差值核函数矩阵K_ST的权重系数w_s，按照如下步骤进行：

5a)通过特征值分解法得到由K_ST*K_Col特征向量构成的第一投影矩阵V和特征值构成的第一对角矩阵λ，K_ST*K_Col、V、λ这三者的关系如下：

(K_ST*K_Col)*V＝λ*V；

5b)提取第一投影矩阵V的第一列，记为α，定义颜色差值核函数矩阵K_Col融合系数w_c：

$w_c=\frac{1}{\sqrt{{\mathrm{\α}}^T{K}_{Col}\mathrm{\α}}}{V}^T;$

其中，T表示矩阵转置，K_Col是颜色差值核函数矩阵；

5c)通过特征值分解法得到矩阵K_Col*K_ST特征向量构成的第二投影矩阵U和特征值构成的第二对角矩阵λ′，K_Col*K_ST、U、λ′这三者的关系如下：

(K_Col*K_ST)*U＝λ′*U；

5d)提取第二投影矩阵U的第一列，记为β；定义强度纹理插值核函数矩阵融合系数w_c：

$w_s=\frac{1}{\sqrt{{\mathrm{\β}}^T{K}_{ST}\mathrm{\β}}}{U}^T;$

其中，K_ST是强度纹理差值核函数矩阵。

7.根据权利要求1所述方法，其中步骤6)中利用步骤5)、步骤4)以及步骤3)的结果构造合成核函数矩阵K_F，按照如下公式进行：

K_F＝w_c·K_Col+w_s·K_ST。