CN103810710A - 基于半监督降维和显著图的多光谱图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于半监督降维和显著图的多光谱图像变化检测方法，解决现有技术不能准确反映多光谱图像光谱特性关系的缺点，其实现步骤为：(1)输入两时相多光谱图像；(2)预处理；(3)生成对应波段差异图；(4)半监督降维；(5)生成显著图；(6)K均值聚类；(7)生成变化检测结果图；(8)输出变化检测结果图。本方法既可以较好的保持变化区域的边缘信息，又可以较好的兼顾变化检测结果中的漏检信息和虚警信息，具有较好的实时性和较高的变化检测结果精度。本发明可以应用于城区扩展监测、森林和植被灾害监测、农作物生长状态的动态监测、军事侦察等领域。

Description

基于半监督降维和显著图的多光谱图像变化检测方法

技术领域

本发明属于图像处理技术领域，更进一步涉及遥感图像处理技术领域中的一种基于半监督降维和显著图的多光谱图像变化检测方法。本发明可应用于土地利用、植被覆盖、水资源以及矿产资源的勘探等方面，能够快速的检测出两时相多光谱遥感图像的变化信息。

背景技术

变化检测通过分析同一地区不同时期拍摄的多幅遥感图像，检测出该地区的地物随时间发生变化的信息，其在国民经济和国防建设的很多领域得到了应用，是遥感图像信息处理中的热点问题。不同地物类型的变化可能反映在不同的光谱范围内，而多光谱遥感图像数据具有从可见光到红外光波段的多个接收频段，丰富的光谱信息增加了识别多种类型变化的可能性和可信度，因此多光谱遥感图像在变化检测中的应用越来越广泛。

北方工业大学在其申请的专利“一种多光谱图像建筑物变化检测方法”（专利申请号：201210587576.3，公开号：CN103077515A）中提出了一种应用多光谱图像检测建筑物变化的变化检测方法。该方法实现的具体步骤是，首先，对两时相多光谱图像进行基于像素比值法的变化检测，得到建筑物变化候选区域；其次，从第一时相和第二时相中将建筑物候选区域分割出来，对去除建筑物候选区域的第一时相图像和第二时相图像逐像素比较得到建筑物变化区域；最后，通过得到的建筑物变化区域指导原始两时相多光谱图像得到最终的变化检测结果。该方法虽然能够较好地利用多个光谱波段的信息，但仍存在的不足是，该方法对于单次变化区域确定的精度要求较高，容易造成误检过多的情况。

河海大学在其专利申请“一种多光谱遥感影像变化检测方法”（专利申请号：201310097259，公开号：CN103218807A）中提出了一种基于支持矢量机（Support vectormachine,SVM）混合核函数的多光谱遥感影像变化检测方法。该方法实现的具体步骤是，首先，采用主成分分析（Principal component analysis,PCA）变换和相关系数融合法相结合的方式构造差异影像；其次，通过支持矢量机对从差异影像中提取的特征进行训练，得到支持矢量机混合核分类器；最后，采用该分类器对差异影像进行分类得到最终的变化检测结果。该方法虽然减少了波段间的冗余信息，增强了差异信息，但仍存在的不足是，该方法中采用的主成分分析变换由于提取的仅为主要成分的信息，容易造成变化信息的丢失。此外，支持矢量机方法对参与训练的样本的要求比较高，变化检测结果的随机性较大。

Marchesi和Bruzzone发表的“ICA and kernel ICA for change detection inmultispectral remote sensing”（IEEE International Geoscience and Remote SensingSymposium,2009,980-983）论文中提出一种基于独立成分分析（Independent componentanalysis,ICA）的多光谱遥感图像变化检测方法。该方法实现的具体步骤是，首先，对各波段差异图像作独立成分分析，使不同类型的变化区域分布在若干相互独立的分量图像上；其次，人工选择几个分量图像用于进行后续分类处理；最后，采用分类分割方法得到最终的变化检测结果。该方法虽然能够从多时相图像中分离出具有非高斯分布的变化区域，但仍存在的不足是，该方法需要人工选取降维后所需要的分量图像，不能实现自动化。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于显著图和半监督降维的多光谱遥感图像变化检测方法。以避免不同时相图像间整体亮度差异的影响，以及变化信息选择误差造成变化信息检测的错误，提高变化检测的精度。

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

(1)输入两时相多光谱图像：

同时输入同一地区、不同时刻的两组多光谱图像，任选其中一组多光谱图像作为第1时相图像，另一组多光谱图像作为第2时相图像；

(2)预处理：

(2a)任选第1时相图像中的一个光谱波段图像作为参考图像，采用加权辐射校正方法，对该参考图像在第2时相图像中对应的光谱波段图像进行辐射校正处理；

(2b)重复步骤(2a)，直至处理完第2时相中的每一幅图像；

(2c)对第1时相图像和加权辐射校正处理后的第2时相图像中的每一幅图像，均采用窗口大小为3×3的维纳滤波进行去噪，得到预处理后的第1时相图像和预处理后的第2时相图像；

(3)生成对应波段差异图：

3a)用预处理后的第1时相图像灰度值减去与其对应的第2时相图像对应波段、对应位置的灰度值，并对差值取绝对值，得到多幅对应波段差异图；

(4)半监督降维：

(4a)求多幅对应波段差异图中每幅对应波段差异图灰度值的平方，得到多幅差异平方图，将多幅差异平方图对应位置灰度值相加，得到一幅差异组合图；

(4b)构造正约束集和负约束集；

(4c)将步骤(3)中的多幅对应波段差异图表示为多幅对应波段差异图矩阵，将多幅对应波段差异图矩阵全部按列表示一维矢量，得到多个一维矢量，将多个一维矢量按行组合成一个多维矢量；

(4d)采用半监督降维方法对多维矢量进行降维，得到一个一维矢量；

(4e)将一维矢量按列表示为尺寸大小与差异组合图尺寸大小相同的一个矩阵，用该矩阵表示半监督降维后差异图；

(5)生成显著图：

(5a)用窗口大小为3×3的高斯滤波，对半监督降维后差异图进行滤波处理，得到滤波后图像；

(5b)计算滤波后图像的灰度均值；

(5c)对半监督降维后差异图的每个像素值均减去滤波后图像的灰度均值，得到一幅显著图；

(6)K均值聚类：

(6a)采用K均值聚类方法，将半监督降维后差异图的像素聚为两类，得到标记图1；

(6b)采用K均值聚类方法，将显著图的像素聚为两类，得到标记图2；

(7)生成变化检测结果图：

(7a)用结构元素为20的膨胀运算，对标记图2进行膨胀，得到膨胀图像；

(7b)将膨胀图像中表示膨胀位置的值替换为标记图2中对应位置的值，将膨胀图像中表示未膨胀位置的值替换为标记图1中对应位置的值，得到替换后的膨胀图像；

(7c)将替换后的膨胀图像表示为变化检测结果图；

(8)输出变化检测结果图。

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

第一、由于本发明使用了加权辐射校正的方法进行预处理，克服了现有技术中由于预处理精度不高造成变化检测结果精度较低的缺点，使得本发明中变化检测结果的虚警数和误检数减少，提高了变化检测结果的精度。

第二、由于本发明使用了半监督降维的方法获取差异图，克服了现有技术中多光谱图像差异图可分性不高的缺点，使得本发明中差异图的可分性更高，提高了变化检测结果的精度。

第三、由于本发明使用计算显著图的方法对变化检测结果进行调整，克服了现有技术中变化检测结果误检较多的缺点，提高了变化检测结果的精度。

附图说明

图1本发明的流程图；

图2本发明实施例中采用的第一组多光谱图像集；

图3本发明实施例中采用的第二组多光谱图像集；

图4本发明的仿真图。

具体实施方式

下面结合附图对本发明做进一步的描述。

结合附图1，对实现本发明的具体步骤描述如下：

步骤1，输入两时相多光谱图像。

同时输入同一地区、不同时刻的两组多光谱图像，任选其中一组多光谱图像作为第1时相多光谱图像，另一组多光谱图像作为第2时相多光谱图像。本发明的实施例中，所选取图像如图2和图3所示。图2中所有的图像和图3中所有的图像均是Landsat7卫星在哥伦比亚比查达省的某一河流区域拍摄的图像。图2(a)是1999年11月20日Landsat7卫星中的蓝色波段拍摄的图像，在本发明中用于表示第1时相第1个光谱波段的图像。图2(b)是1999年11月20日Landsat7卫星中的绿色波段拍摄的图像，在本发明中用于表示第1时相第2个光谱波段的图像。图2(c)是1999年11月20日Landsat7卫星中的红色波段拍摄的图像，在本发明中用于表示第1时相第3个光谱波段的图像。图2(d)是1999年11月20日Landsat7卫星中的近红外波段拍摄的图像，在本发明中用于表示第1时相第4个光谱波段的图像。图2(e)是1999年11月20日Landsat7卫星中的中红外波段拍摄的图像，在本发明中用于表示第1时相第5个光谱波段的图像。图2(f)是1999年11月20日Landsat7卫星中的短波波段拍摄的图像，在本发明中用于表示第1时相第6个光谱波段的图像。图3(a)是2001年3月14日Landsat7卫星中的蓝色波段拍摄的图像，在本发明中用于表示第2时相第1个光谱波段的图像。图3(b)是2001年3月14日Landsat7卫星中的绿色波段拍摄的图像，在本发明中用于表示第2时相第2个光谱波段的图像。图3(c)是2001年3月14日Landsat7卫星中的红色波段拍摄的图像，在本发明中用于表示第2时相第3个光谱波段的图像。图3(d)是2001年3月14日Landsat7卫星中的近红外波段拍摄的图像，在本发明中用于表示第2时相第4个光谱波段的图像。图3(e)是2001年3月14日Landsat7卫星中的中红外波段拍摄的图像，在本发明中用于表示第2时相第5个光谱波段的图像。图3(f)是2001年3月14日Landsat7卫星中的短波波段拍摄的图像，在本发明中用于表示第2时相第6个光谱波段的图像。两时相多光谱图像的大小均为499×290像素。

步骤2，预处理。

任选第1时相图像中的一个光谱波段图像作为参考图像，采用加权辐射校正方法，对该参考图像在第2时相图像中对应的光谱波段图像进行辐射校正处理。重复上述步骤，直至处理完第2时相中的每一幅图像。对第1时相图像和加权辐射校正处理后的第2时相图像中的每一幅图像均采用窗口大小为3×3的维纳滤波进行去噪，得到预处理后的第1时相图像和预处理后的第2时相图像。

加权辐射校正方法的具体步骤如下：

第一步，按照下式，计算待校正图像的加权灰度均值：

$\mathrm{\μ}=\frac{\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y)Z(x,y)}{\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y)}$

其中，μ表示待校正图像的加权灰度均值，m表示待校正图像的行尺寸大小，n表示待校正图像的列尺寸大小，x表示待校正图像的行序号，y表示待校正图像的列序号，w(x,y)表示待校正图像在(x,y)位置处的权重值，Z(x,y)表示待校正图像的在(x,y)位置处的灰度值；

第二步，按照下式，计算待校正图像的加权灰度方差：

$\mathrm{\σ}=\frac{1}{\frac{m\×n-1}{m\×n}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y)}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y){(Z(x,y)-\mathrm{\μ})}^2$

其中，σ表示待校正图像的加权灰度方差，m表示待校正图像的行尺寸大小，n表示待校正图像的列尺寸大小，w(x,y)表示待校正图像在(x,y)位置处的权重值，x表示待校正图像的行序号，y表示待校正图像的列序号，Z(x,y)表示待校正图像在(x,y)位置处的灰度值，μ表示待校正图像的加权灰度均值；

第三步，按照下式，计算待校正图像的辐射校正灰度值：

$I(x,y)=\frac{Z(x,y)-\mathrm{\μ}}{\mathrm{\σ}}\×\mathrm{\α}+\mathrm{\δ}$

其中，I(x,y)表示待校正图像在(x,y)位置处经过辐射校正后得到的灰度值，Z(x,y)表示待校正图像在(x,y)位置处的灰度值，x表示待校正图像的行序号，y表示待校正图像的列序号，μ表示待校正图像的加权灰度均值，σ表示待校正图像的加权灰度方差，α表示参考图像的灰度均值，δ表示参考图像的灰度方差；

第四步，重复第一步至第三步，直至处理完待校正图像中的所有像素点，得到加权辐射校正处理后的图像。

步骤3，生成对应波段差异图组。

用预处理后的第1时相图像灰度值减去与其对应的第2时相图像对应波段、对应位置的灰度值，并对差值取绝对值，得到多幅对应波段差异图。

步骤4，半监督降维。

求多幅对应波段差异图中每幅对应波段差异图灰度值的平方，得到多幅差异平方图，将多幅差异平方图对应位置灰度值相加，得到一幅差异组合图。根据差异图组合图中的像素构造正约束集和负约束集。将步骤3中的多幅对应波段差异图表示为多幅对应波段差异图矩阵，将多幅对应波段差异图矩阵全部按列表示一维矢量，得到多个一维矢量，将多个一维矢量按行组合成一个多维矢量。采用半监督降维方法对多维矢量进行降维，得到一个一维矢量。将一维矢量按列表示为尺寸大小与差异组合图尺寸大小相同的一个矩阵，用该矩阵表示半监督降维后差异图。

构造正约束集和负约束集的具体步骤如下：

第一步，采用K均值聚类算法，将差异组合图中的像素分为三类，得到初始类别标记图，将第一类作为变化类，第二类作为不确定类，第三类作为未变化类。

第二步，将初始类别标记图按列表示为一维矢量，得到初始类别标记矢量。

第三步，随机抽取初始类别标记矢量中变化类和未变化类的列序号，用属于同一类别的两个列序号构造正约束，用属于不同类别的两个列学号构造负约束。

第四步，重复第三步，将多个正约束组合为一个二维的正约束集，将多个负约束组合一个二维的负约束集。

半监督降维方法的具体步骤如下：

第一步，按照下式，计算多维矢量中任意两个元素的权重值：

$S_{\mathrm{xy}}=\left\{\begin{array}{cc}\frac{1}{d^2}+\frac{1}{r}& \mathrm{if}(x,y)\∈C\\ \frac{1}{d^2}-\frac{\mathrm{\β}}{t}& \mathrm{if}(x,y)\∈M\\ \frac{1}{d^2}& \mathrm{otherwise}\end{array}\right.$

其中，S_xy表示多维矢量中第x个元素和第y个元素的权重值，x表示多维矢量的列序号，y表示多维矢量的列序号，d表示步骤(4b)中正约束个数和负约束个数的总和，r表示步骤(4b)中负约束的个数，t表示步骤(4b)中正约束的个数，β表示尺度参数，C表示负约束集，M表示正约束集。

第二步，重复第一步，直至处理完多维矢量中的所有元素，用第x个元素和第y个元素的权重值表示矩阵的第x行第y列的值，得到一个相似性矩阵。

第三步，按照下式，计算相似性矩阵中节点的度：

$d_i=\underset{j}{\mathrm{\Σ}}{S}_{\mathrm{ij}}$

其中，d_i表示相似性矩阵中第i个节点的度，i表示相似性矩阵的行序号，S_ij表示相似性中第i行第j列的值，j表示相似性矩阵的列序号。

用对角矩阵表示相似性矩阵的度矩阵如下：

其中，D表示相似性矩阵的度矩阵，d₁表示相似性矩阵中第1个节点的度，d₂表示相似性矩阵中第2个节点的度，d_Z表示相似性矩阵中第Z个节点的度。

第四步，按照下式，计算图拉普拉斯矩阵：

L=D-F

其中，L表示图拉普拉斯矩阵，D表示相似性矩阵的度矩阵，F表示相似性矩阵。

第五步，计算HLH^T的特征值及特征值所对应的特征矢量，其中H表示多维矢量，L表示图拉普拉斯矩阵，H^T表示多维矢量H的转置，T表示转置符号。

第六步，按照下式，计算多维矢量经过降维处理后的一维矢量：

P=H×W

其中，P表示多维矢量经过降维处理后的一维矢量，H表示多维矢量，W表示特征值中最大特征值所对应的特征矢量。

步骤5，生成显著图。

用窗口大小为3×3的高斯滤波，对半监督降维后差异图进行滤波处理，得到滤波后图像。计算滤波后图像的灰度均值。对半监督降维后差异图的每个像素值均减去滤波后图像的灰度均值，得到一幅显著图。

步骤6，K均值聚类。

采用K均值聚类算法，将半监督降维后差异图的像素聚为两类，得到标记图1。采用K均值聚类算法，将显著图的像素聚为两类，得到标记图2。

步骤7，生成变化检测结果图。

用结构元素为20的膨胀运算，对标记图2进行膨胀，得到膨胀图像。将膨胀图像中表示膨胀位置的值替换为标记图2中对应位置的值，将膨胀图像中表示未膨胀位置的值替换为标记图1中对应位置的值，得到替换后的膨胀图像。将替换后的膨胀图像表示为变化检测结果图。

步骤8，输出变化检测结果图。

下面结合仿真图对本发明的效果做进一步说明。

1.仿真实验条件：

仿真实验在CPU为Intel(R)Pentium(R)Dual、主频2.16GHz，内存为2G的WINDOWS XP系统上用MATLAB7.14.0.0软件进行。

本发明仿真实验所用数据为一组真实多光谱遥感数据，其样例图像如图2和图3所示。图2中所有的图像和图3中所有的图像均是Landsat7卫星在哥伦比亚比查达省的某一河流区域拍摄的图像。图2中所有的图像和图3中所有的图像均是Landsat7卫星在哥伦比亚比查达省的某一河流区域拍摄的图像。图2(a)是1999年11月20日Landsat7卫星中的蓝色波段拍摄的图像，在本发明中用于表示第1时相第1个光谱波段的图像。图2(b)是1999年11月20日Landsat7卫星中的绿色波段拍摄的图像，在本发明中用于表示第1时相第2个光谱波段的图像。图2(c)是1999年11月20日Landsat7卫星中的红色波段拍摄的图像，在本发明中用于表示第1时相第3个光谱波段的图像。图2(d)是1999年11月20日Landsat7卫星中的近红外波段拍摄的图像，在本发明中用于表示第1时相第4个光谱波段的图像。图2(e)是1999年11月20日Landsat7卫星中的中红外波段拍摄的图像，在本发明中用于表示第1时相第5个光谱波段的图像。图2(f)是1999年11月20日Landsat7卫星中的短波波段拍摄的图像，在本发明中用于表示第1时相第6个光谱波段的图像。图3(a)是2001年3月14日Landsat7卫星中的蓝色波段拍摄的图像，在本发明中用于表示第2时相第1个光谱波段的图像。图3(b)是2001年3月14日Landsat7卫星中的绿色波段拍摄的图像，在本发明中用于表示第2时相第2个光谱波段的图像。图3(c)是2001年3月14日Landsat7卫星中的红色波段拍摄的图像，在本发明中用于表示第2时相第3个光谱波段的图像。图3(d)是2001年3月14日Landsat7卫星中的近红外波段拍摄的图像，在本发明中用于表示第2时相第4个光谱波段的图像。图3(e)是2001年3月14日Landsat7卫星中的中红外波段拍摄的图像，在本发明中用于表示第2时相第5个光谱波段的图像。图3(f)是2001年3月14日Landsat7卫星中的短波波段拍摄的图像，在本发明中用于表示第2时相第6个光谱波段的图像。两时相多光谱图像的大小均为499×290像素。两时相多光谱图像之间发生的变化是由水流冲击河床造成的河流区域变化所致，包括10015个变化像素和134695个未变化的像素。

2.仿真内容及分析：

本发明与基于变化矢量分析的多光谱图像变化检测方法进行对比，用于验证本发明获取的差异图在抑制噪声冗余和减少整体亮度差异干扰方面的优越性，将该方法简记为CVA法。本发明与基于主成分分析的多光谱图像变化检测方法进行对比，用于验证本发明方法中采用的半监督降维方法在获取高可分性差异图方面的优越性，将该方法简记为PCA方法。本发明与基于半监督降维的多光谱图像变化检测方法进行对比，用于验证本发明方法中采用的显著图调整方法在提高变化检测结果精度中的优越性，将该方法简记为半监督降维方法。

图4为本发明的仿真结果图。其中，图4(a)为现有技术采用CVA方法的效果图，图4(b)为现有技术采用PCA方法的效果图，图4(c)为现有技术采用半监督降维方法的效果图，图4(d)为本发明的效果图。

从图4可以看出，采用现有技术采用CVA方法、PCA方法或半监督降维方法，对多光谱图像进行变化检测处理得到的结果图中均存在着大面积的伪变化区域，而且在边缘保持方面的效果也不佳。通过本发明得到的变化检测结果图存在较少的误检信息，较好地保持了变化区域的边缘信息。

本发明通过下表中的误检数、漏检数、总错误数和正确率这四个指标评价变化检测效果的好坏，其中误检数、漏检数和总错误数的值越低，变化检测效果越好，而正确率的值越高，变化检测效果越好。

方法	误检数	漏检数	总错误数	正确率
					CVA方法	7867	754	8621	94.04%
PCA方法	8141	1097	9238	93.62%
					半监督降维方法	2398	1061	3459	97.61%
本发明	921	1379	2300	98.41%

从上表可以看出，相对于其他三种方法，本发明的误检数和总错误数是最少的，正确率是最高的，即本发明方法能够有效地抑制不同时相图像间整体亮度差异及噪声对多光谱图像的影响，提高了变化检测结果的精度。

Claims (4)

1.基于半监督降维和显著图的多光谱图像变化检测方法，包括如下步骤：

(1)输入两时相多光谱图像：

同时输入同一地区、不同时刻的两组多光谱图像，任选其中一组多光谱图像作为第1时相图像，另一组多光谱图像作为第2时相图像；

(2)预处理：

(2a)任选第1时相图像中的一个光谱波段图像作为参考图像，采用加权辐射校正方法，对该参考图像在第2时相图像中对应的光谱波段图像进行辐射校正处理；

(2b)重复步骤(2a)，直至处理完第2时相中的每一幅图像；

(2c)对第1时相图像和加权辐射校正处理后的第2时相图像中的每一幅图像，均采用窗口大小为3×3的维纳滤波进行去噪，得到预处理后的第1时相图像和预处理后的第2时相图像；

(3)生成对应波段差异图：

3a)用预处理后的第1时相图像灰度值减去与其对应的第2时相图像对应波段、对应位置的灰度值，并对差值取绝对值，得到多幅对应波段差异图；

(4)半监督降维：

(4a)求多幅对应波段差异图中每幅对应波段差异图灰度值的平方，得到多幅差异平方图，将多幅差异平方图对应位置灰度值相加，得到一幅差异组合图；

(4b)构造正约束集和负约束集；

(4c)将步骤(3)中的多幅对应波段差异图表示为多幅对应波段差异图矩阵，将多幅对应波段差异图矩阵全部按列表示一维矢量，得到多个一维矢量，将多个一维矢量按行组合成一个多维矢量；

(4d)采用半监督降维方法对多维矢量进行降维，得到一个一维矢量；

(4e)将一维矢量按列表示为尺寸大小与差异组合图尺寸大小相同的一个矩阵，用该矩阵表示半监督降维后差异图；

(5)生成显著图：

(5a)用窗口大小为3×3的高斯滤波，对半监督降维后差异图进行滤波处理，得到滤波后图像；

(5b)计算滤波后图像的灰度均值；

(5c)对半监督降维后差异图的每个像素值均减去滤波后图像的灰度均值，得到一幅显著图；

(6)K均值聚类：

(6a)采用K均值聚类方法，将半监督降维后差异图的像素聚为两类，得到标记图1；

(6b)采用K均值聚类方法，将显著图的像素聚为两类，得到标记图2；

(7)生成变化检测结果图：

(7a)用结构元素为20的膨胀运算，对标记图2进行膨胀，得到膨胀图像；

(7b)将膨胀图像中表示膨胀位置的值替换为标记图2中对应位置的值，将膨胀图像中表示未膨胀位置的值替换为标记图1中对应位置的值，得到替换后的膨胀图像；

(7c)将替换后的膨胀图像表示为变化检测结果图；

(8)输出变化检测结果图。

2.根据权利要求1所述的基于半监督降维和显著图的多光谱图像变化检测方法，其特征在于，步骤(2a)中所述的加权辐射校正方法，按照如下步骤进行：

第一步，按照下式，计算待校正图像的加权灰度均值：

$\mathrm{\μ}=\frac{\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y)Z(x,y)}{\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y)}$

其中，μ表示待校正图像的加权灰度均值，m表示待校正图像的行尺寸大小，n表示待校正图像的列尺寸大小，x表示待校正图像的行序号，y表示待校正图像的列序号，w(x,y)表示待校正图像在(x,y)位置处的权重值，Z(x,y)表示待校正图像的在(x,y)位置处的灰度值；

第二步，按照下式，计算待校正图像的加权灰度方差：

$\mathrm{\σ}=\frac{1}{\frac{m\×n-1}{m\×n}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y)}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}w(x,y){(Z(x,y)-\mathrm{\μ})}^2$

其中，σ表示待校正图像的加权灰度方差，m表示待校正图像的行尺寸大小，n表示待校正图像的列尺寸大小，w(x,y)表示待校正图像在(x,y)位置处的权重值，x表示待校正图像的行序号，y表示待校正图像的列序号，Z(x,y)表示待校正图像在(x,y)位置处的灰度值，μ表示待校正图像的加权灰度均值；

第三步，按照下式，计算待校正图像的辐射校正灰度值：

$I(x,y)=\frac{Z(x,y)-\mathrm{\μ}}{\mathrm{\σ}}\×\mathrm{\α}+\mathrm{\δ}$

其中，I(x,y)表示待校正图像在(x,y)位置处经过辐射校正后得到的灰度值，Z(x,y)表示待校正图像在(x,y)位置处的灰度值，x表示待校正图像的行序号，y表示待校正图像的列序号，μ表示待校正图像的加权灰度均值，σ表示待校正图像的加权灰度方差，α表示参考图像的灰度均值，δ表示参考图像的灰度方差；

第四步，重复第一步至第三步，直至处理完待校正图像中的所有像素点，得到加权辐射校正处理后的图像。

3.根据权利要求1所述的基于半监督降维和显著图的多光谱图像变化检测方法，其特征在于，步骤(4b)所述的构造正约束集和负约束集，按照如下步骤进行：

第一步，采用K均值聚类方法，将差异组合图中的像素分为三类，得到初始类别标记图，将第一类作为变化类，第二类作为不确定类，第三类作为未变化类；

第二步，将初始类别标记图按列表示为一维矢量，得到初始类别标记矢量；

第三步，随机抽取初始类别标记矢量中变化类和未变化类的列序号，用属于同一类别的两个列序号构造正约束，用属于不同类别的两个列学号构造负约束；

第四步，重复第三步，得到多个正约束和多个负约束；

第五步，将多个正约束组合为一个二维的正约束集，将多个负约束组合一个二维的负约束集。

4.根据权利要求1所述的基于半监督降维和显著图的多光谱图像变化检测方法，其特征在于，步骤(4d)所述的半监督降维方法，按照如下步骤进行：

第一步，按照下式，计算多维矢量中任意两个元素的权重值：

$S_{\mathrm{xy}}=\left\{\begin{array}{cc}\frac{1}{d^2}+\frac{1}{r}& \mathrm{if}(x,y)\∈C\\ \frac{1}{d^2}-\frac{\mathrm{\β}}{t}& \mathrm{if}(x,y)\∈M\\ \frac{1}{d^2}& \mathrm{otherwise}\end{array}\right.$

其中，S_xy表示多维矢量中第x个元素和第y个元素的权重值，x表示多维矢量的列序号，y表示多维矢量的列序号，d表示步骤(4b)中正约束个数和负约束个数的总和，r表示步骤(4b)中负约束的个数，t表示步骤(4b)中正约束的个数，β表示尺度参数，C表示负约束集，M表示正约束集；

第二步，重复第一步，直至处理完多维矢量中的所有元素，用第x个元素和第y个元素的权重值表示矩阵的第x行第y列的值，得到一个相似性矩阵；

第三步，按照下式，计算相似性矩阵中节点的度：

$d_i=\underset{j}{\mathrm{\Σ}}{S}_{\mathrm{ij}}$

其中，d_i表示相似性矩阵中第i个节点的度，i表示相似性矩阵的行序号，S_ij表示相似性中第i行第j列的值，j表示相似性矩阵的列序号；

用对角矩阵表示相似性矩阵的度矩阵如下：

其中，D表示相似性矩阵的度矩阵，d₁表示相似性矩阵中第1个节点的度，d₂表示相似性矩阵中第2个节点的度，d_Z表示相似性矩阵中第Z个节点的度；

第四步，按照下式，计算图拉普拉斯矩阵：

L=D-F

其中，L表示图拉普拉斯矩阵，D表示相似性矩阵的度矩阵，F表示相似性矩阵；

第五步，计算HLH^T的特征值及特征值所对应的特征矢量，其中H表示多维矢量，L表示图拉普拉斯矩阵，H^T表示多维矢量H的转置，T表示转置符号；

第六步，按照下式，计算多维矢量经过降维处理后的一维矢量：

P=H×W

其中，P表示多维矢量经过降维处理后的一维矢量，H表示多维矢量，W表示特征值中最大特征值所对应的特征矢量。

Images