CN104463881A - 一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其实现步骤主要包括：1)输入同一地区的已配准的两时相多光谱图像集，对其进行维纳滤波去噪并归一化处理；2)将处理后的图像集转换为相对地物光谱反射率图像集；3)计算光谱反射率变化模值，得到光谱反射率变化模值图；4)计算光谱反射率角值，得到光谱反射率角制图；5)计算光谱反射率的邻域差异值，得到邻域差异图；6)分别对光谱反射率变化模值图、光谱反射率角制图和邻域差异图进行聚类，得到二值图；7)对所得二值图进行基于邻域概率的融合，得到变化检测结果图。本发明无需人工参与，检测精确度高，可用于城区扩展监测、森林和植被变化监测等领域。

Description

一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法

技术领域

本发明属于光学遥感图像处理技术领域，具体为一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，可用于土地利用、植被覆盖等方面的变化监测。

背景技术

遥感图像的变化检测是通过分析和提取同一地区不同时相的遥感图像之间存在的电磁波谱特征差异或空间结构特征差异，来识别物体的状态变化或现象变化的过程。在国民经济和国防建设的诸多领域已得到广泛应用，如农业调查、森林和植被变化监测等。由于不同地物类型的变化可能会反映在不同的光谱范围内，而多光谱遥感图像数据具有从可见光到红外光波段的多个接收频段，丰富的光谱信息增加了识别多种类型变化的可能性和可信度，因此多光谱遥感图像被广泛用于地物的变化检测。

多光谱遥感图像变化检测最为常见的是基于变化矢量分析(CVA)的变化检测，一般是先通过CVA方法获取差异图，然后对差异图进行变化与未变化的分类。由于同时考虑了多个波段图像的变化信息，CVA相对于单一波段的变化检测增强了检测结果的可靠性。然而，CVA在计算光谱变化矢量的模值时只是简单地将各波段差异信息叠加，由于配准误差等原因，对于没有发生变化的区域，两时相图像对应位置的灰度值也存在一定的偏差，故而简单地对CVA方法获取的差异图像进行分类会导致变化检测结果中存在很多的伪变化信息。在变化检测中引入光谱角制图(SAM)可以降低伪变化信息，SAM方法通过计算两时相光谱矢量之间的角度来定义两者的差异性，然后对两时相差异图进行变化与未变化的分类。SAM方法利用了角度这一参数，减少了光照对不同时相图像间整体亮度差异的影响，从而提高了检测准确率。然而，该方法用于多光谱图像时，由于其波段信息较少，导致不同时相图像间对应地物光谱矢量之间的差异不够大，因此对SAM方法获取的差异图像进行分类所得变化检测结果中，虽然降低了伪变化信息，但同时引入了很多的漏检。

上述的变化检测方法都利用了图像的光谱信息，没有考虑像素的邻域信息，使得变化检测结果不理想。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，综合考虑两时相图像像素的光谱信息和邻域信息，并在此基础上结合上述技术，从而提高检测结果的正确率。

为实现上述目的，本发明提出了一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其技术方案是：

1、一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其特征在于，包括如下步骤：

(1)输入在两个时相获取的同一地区的两个多光谱图像集：I₁＝{A₁ ^b}和I₂＝{A₂ ^b}，其中，I₁为时相一多光谱图像集,I₂为时相二多光谱图像集,A_t ^b为两个多光谱图像集中的每一幅单波段图像，上标b表示波段序号，b＝1,2,…,B，B为总波段数，下标t为时相序号，t＝{1,2}，每一幅单波段图像A_t ^b均由m行n列像素构成；

(2)对两时相多光谱图像集I₁和I₂分别进行维纳滤波去噪，并归一化处理，得到两时相多光谱归一化图像集和

(3)分别对归一化图像集和中每一幅单波段图像，采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到两个时相每一波段的相对地物光谱反射率图像R₁ ^b和R₂ ^b，并将相同时相的每一波段相对地物光谱反射率图像分别构成两个时相的相对地物光谱反射率图像集R₁和R₂；

(4)分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率变化模值mr(x,y)，得到光谱反射率变化模值图Mr＝{mr(x,y)}，其中x为行序号，y为列序号；

(5)分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y)，得到光谱反射率角制图Sr＝{sr(x,y)}；

(6)分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y)，得到邻域差异图Qr＝{qr(x,y)}；

(7)用谱聚类方法将光谱反射率变化模值图Mr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM1；

(8)用谱聚类方法将光谱反射率角制图Sr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM2；

(9)用谱聚类方法将邻域差异图Qr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM3；

(10)对CM1、CM2和CM3采用基于邻域概率的融合方法，得到最终变化检测结果图CM。

上述步骤(5)中计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y)，按如下公式计算：

对两时相各波段的相对地物光谱反射率图像{R₁ ^b}和{R₂ ^b}中的任意空间位置像素点(x,y)，该点的光谱反射率角值sr(x,y)的计算公式如下；

$\mathrm{sr}(x,y)=\mathrm{ar}\cos \left(\frac{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}(r_1^b(x,y)\·r_2^b(x,y))}{\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left(r_1^b(x,y)\right)}^2}\·\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left(r_2^b(x,y)\right)}^2}}\right).$

上述步骤(6)中计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y)，通过如下步骤进行：

(6a)对第t时相第b波段的相对地物光谱反射率图像{R_t ^b}中的任意空间位置像素点(x,y)，计算该点的k×k邻域内的相对地物反射率值的均值u_t ^b(x,y)：

$u_t^b(x,y)=\frac{1}{k\·k}\underset{v=-c}{\overset{c}{\mathrm{\Σ}}}\underset{w=-c}{\overset{c}{\mathrm{\Σ}}}{r}_t^b(x+v,y+w);$

其中c＝(k-1)/2，k的取值为奇数，取k＝3或k＝5，v和w分别表示以中心像素(x,y)为坐标原点时其对应的k×k邻域内像素点的横坐标值与纵坐标值；

(6b)计算两时相第b波段像素点(x,y)的k×k邻域内的邻域变化值g^b(x,y)；

$g^b(x,y)=\frac{1}{k\·k}\underset{v=-c}{\overset{c}{\mathrm{\Σ}}}\underset{w=-c}{\overset{c}{\mathrm{\Σ}}}{(r_1^b(x+v,y+w)-r_2^b(x+v,y+w))}^2;$

(6c)计算像素点(x,y)的邻域差异值qr(x,y)；

$\mathrm{qr}(x,y)=(1-\frac{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\min (r_1^b(x,y)-u_1^b(x,y),r_2^b(x,y)-u_2^b(x,y))}{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\max (r_1^b(x,y)-u_1^b(x,y),r_2^b(x,y)-u_2^b(x,y))})\·\left(\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\frac{g^b(x,y)}{u_1^b(x,y)}\right);$

上述步骤(10)中所述的对CM1、CM2和CM3采用基于邻域概率的融合方法，得到最终变化检测结果图CM，按如下步骤进行：

(10a)计算CM1中像素点(x,y)的k×k邻域内变化类的个数c1(x,y)和非变化类的个数nc1(x,y)，相对应的概率为p_c1(x,y)和p_n1(x,y)；

$p_{c1}(x,y)=\frac{c1(x,y)}{c1(x,y)+\mathrm{nc}1(x,y)};$

$p_{n1}(x,y)=\frac{\mathrm{nc}1(x,y)}{c1(x,y)+\mathrm{nc}1(x,y)};$

(10b)计算CM2中像素点(x,y)的k×k邻域内变化类的个数c2(x,y)和非变化类的个数nc2(x,y)，相应的概率为p_c2(x,y)和p_n2(x,y)；

$p_{c2}(x,y)=\frac{c2(x,y)}{c2(x,y)+\mathrm{nc}2(x,y)};$

$p_{n2}(x,y)=\frac{\mathrm{nc}2(x,y)}{c2(x,y)+\mathrm{nc}2(x,y)};$

(10c)计算CM3中像素点(x,y)的k×k邻域内变化类的个数c3(x,y)和非变化类的个数nc3(x,y)，相应的概率为p_c3(x,y)和p_n3(x,y)；

$p_{c3}(x,y)=\frac{c3(x,y)}{c3(x,y)+\mathrm{nc}3(x,y)};$

$p_{n3}(x,y)=\frac{\mathrm{nc}3(x,y)}{c3(x,y)+\mathrm{nc}3(x,y)};$

(10d)计算像素点(x,y)属于变化类的混合概率cw(x,y)和属于非变化类的混合概率nw(x,y)；

cw(x,y)＝c1(x,y)·p_c1(x,y)+c2(x,y)·p_c2(x,y)+c3(x,y)·p_c3(x,y)；

nw(x,y)＝nc1(x,y)·p_n1(x,y)+nc2(x,y)·p_n2(x,y)+nc3(x,y)·p_n3(x,y)；

(10e)对像素点(x,y)，如果cw(x,y)大于等于nw(x,y)，则该像素点为变化类，否则，该像素点为非变化类；

(10f)对CM1、CM2和CM3中的所有m×n个像素均重复步骤10a)至步骤10e)，得到最终变化检测结果图CM。

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

1、本发明综合利用了两时相图像像素的光谱信息和邻域信息，提出并采用了光谱反射率邻域差异值作为衡量地物变化与否的新指标，提高了检测的能力；另外，本发明将图像灰度信息转化为光谱反射率信息，降低了不同时相图像间存在的整体亮度差异对检测稳定性和准确性的影响。

2、本发明利用了光谱反射率变化模值图、光谱反射率角制图和邻域差异图在变化检测中各自的优点，提出并使用了基于邻域概率的融合方法，提高了变化检测结果的正确率。

附图说明

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

图2是用于实验的第一组多光谱图像集和对应的变化参考图像；

图3是用于实验的第二组多光谱图像集和对应的变化参考图像；

图4是本发明与对比方法仿真第一组多光谱图像集得到的变化检测结果图；

图5是本发明与对比方法仿真第二组多光谱图像集得到的变化检测结果图。

具体实施方式

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

步骤1，输入在两个时相获取的同一地区的两个多光谱图像集I₁和I₂。

输入在两个时相获取的同一地区的两个多光谱图像集：I₁＝{A₁ ^b}和I₂＝{A₂ ^b}，其中，A_t ^b为两个多光谱图像集合中的每一幅单波段图像，上标b表示波段序号，b＝1,2,…,B，B为总波段数，下标t为时相序号，t＝{1,2}，每一幅单波段图像A_t ^b均由m行n列像素构成。

步骤2，对两个时相多光谱图像集I₁和I₂分别进行维纳滤波去噪并归一化处理，得到归一化图像集和

2.1)将两时相多光谱图像集I₁和I₂中每一幅图像A_t ^b的灰度值区间从0～255转换到0～1之间，然后采用窗口大小为3×3像素的维纳滤波进行去噪，得到相应的去噪后图像 ${{\stackrel{..}{A}}_t}^b=\left\{{{\stackrel{..}{A}}_t}^b(x,y)\right|x=\mathrm{1,2},...,m,y=\mathrm{1,2},...,n\},$ 其中，为去噪后图像中的像素点(x,y)的灰度值；

2.2)对去噪后图像进行归一化处理，得到归一化图像 其中，为归一化图像中的像素点(x,y)的灰度值，其计算公式如下：

${\hat{A}}_t^b(x,y)=\frac{{\stackrel{..}{A}}_t^b(x,y)-\min \left(\right\{{\stackrel{..}{A}}_t^b(x,y)\left\}\right)}{\max \left(\right\{{\stackrel{..}{A}}_t^b(x,y)\left\}\right)-\min \left(\right\{{\stackrel{..}{A}}_t^b(x,y)\left\}\right)},$

其中，表示取去噪后图像中所有灰度值中的最大值，表示取去噪后图像中所有灰度值中的最小值，x为行序号，y为列序号；

2.3)用时相一和时相二的归一化图像分别构成归一化图像集和即

步骤3，分别对归一化图像集和中每一幅单波段图像，采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到两个时相每一波段的相对地物光谱反射率图像R₁ ^b和R₂ ^b，构成两个时相的相对地物光谱反射率图像集R₁和R₂。

3.1)对归一化图像集和中第t时相第b波段的归一化图像将该图像中的像素点(x,y)的灰度值转换为相对地物光谱反射率值r_t ^b(x,y)的对数：

$\begin{array}{c}\log r_t^b(x,y)=\log {\hat{A}}_t^b(x,y)-\log \frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\hat{A}}_t^b(x,y)-\log \frac{1}{n^2}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}{\hat{A}}_t^b(x,y)\\ +\log \frac{1}{m\×n\×B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\underset{y=1}{\overset{n}{\mathrm{\Σ}}}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}{\hat{A}}_t^b(x,y);\end{array}$

其中，x和y分别为归一化图像中像素的行序号和列序号，x＝1,2,…,m，y＝1,2,…,n；

3.2)对logr_t ^b(x,y)做指数变换，得到像素点(x,y)的相对地物光谱反射率值r_t ^b(x,y)；

3.3)对第t时相第b波段的归一化图像的所有m×n个像素点均重复步骤3.1)和步骤3.2)，得到第t时相第b波段的相对地物光谱反射率图像R_t ^b＝{r_t ^b(x,y)}，t＝{1,2}；

3.4)对时相一所有B个波段的归一化图像重复步骤3.1)至步骤3.3)，得到时相一的B个波段的相对地物光谱反射率图像{R₁ ^b}，构成时相一的相对地物光谱反射率图像集R₁＝{R₁ ^b|b＝1,2,…,B}；

3.5)对时相二所有B个波段的归一化图像重复步骤3.1)至步骤3.3)，得到时相二的B个波段的相对地物光谱反射率图像{R₂ ^b}，构成时相二的相对地物光谱反射率图像集R₂＝{R₂ ^b|b＝1,2,…,B}。

步骤4，分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率变化模值mr(x,y)，得到光谱反射率变化模值图Mr＝{mr(x,y)}。

4.1)将时相一各波段的相对地物光谱反射率图像{R₁ ^b}与时相二对应波段的相对地物光谱反射率图像{R₂ ^b}中任意一点(x,y)的像素值作差，得到该像素点的光谱反射率变化量Dr(x,y)＝{dr^b(x,y)|b＝1,2,…,B}，其中dr^b(x,y)＝r₁ ^b(x,y)–r₂ ^b(x,y)；

4.2)计算像素点(x,y)的光谱反射率变化模值mr(x,y)；

$\mathrm{mr}(x,y)=\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left({\mathrm{dr}}^b(x,y)\right)}^2};$

4.3)对两时相的相对地物光谱反射率图像集R₁和R₂中的所有m×n个像素均重复步骤4.1)和步骤4.2)，得到光谱反射率变化模值图Mr＝{mr(x,y)│x＝1,2,…,m,y＝1,2,…,n}。

步骤5，分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y)，得到光谱反射率角制图Sr＝{sr(x,y)}。

5.1)对两时相各波段的相对地物光谱反射率图像{R₁ ^b}和{R₂ ^b}中的任意空间位置像素点(x,y)，计算该点的光谱反射率角值sr(x,y)；

$\mathrm{sr}(x,y)=\mathrm{ar}\cos \left(\frac{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}(r_1^b(x,y)\·r_2^b(x,y))}{\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left(r_1^b(x,y)\right)}^2}\·\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left(r_2^b(x,y)\right)}^2}}\right);$

5.2)对两时相的相对地物光谱反射率图像集中的所有m×n个像素均重复步骤5.1)，得到光谱反射率角制图Sr＝{sr(x,y)│x＝1,2,…,m,y＝1,2,…,n}。

步骤6，分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y)，得到邻域差异图Qr＝{qr(x,y)}。

6.1)对第t时相第b波段的相对地物光谱反射率图像{R_t ^b}中的任意空间位置像素点(x,y)，计算该点的k×k邻域内的相对地物反射率值的均值u_t ^b(x,y)；

$u_t^b(x,y)=\frac{1}{k\·k}\underset{v=-c}{\overset{c}{\mathrm{\Σ}}}\underset{w=-c}{\overset{c}{\mathrm{\Σ}}}{r}_t^b(x+v,y+w);$

其中c＝(k-1)/2，k的取值为奇数，k＝3或5，本发明实例中k＝3，v和w分别表示以中心像素(x,y)为坐标原点时其对应的k×k邻域内像素点的横坐标值与纵坐标值；

6.2)计算两时相第b波段像素点(x,y)的k×k邻域内的邻域变化值g^b(x,y)；

$g^b(x,y)=\frac{1}{k\·k}\underset{v=-c}{\overset{c}{\mathrm{\Σ}}}\underset{w=-c}{\overset{c}{\mathrm{\Σ}}}{(r_1^b(x+v,y+w)-r_2^b(x+v,y+w))}^2;$

6.3)计算像素点(x,y)的邻域差异值qr(x,y)；

$\mathrm{qr}(x,y)=(1-\frac{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\min (r_1^b(x,y)-u_1^b(x,y),r_2^b(x,y)-u_2^b(x,y))}{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\max (r_1^b(x,y)-u_1^b(x,y),r_2^b(x,y)-u_2^b(x,y))})\·\left(\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\frac{g^b(x,y)}{u_1^b(x,y)}\right);$

6.4)对两时相的相对地物光谱反射率图像集中的所有m×n个像素均重复步骤6.1)、步骤6.2)和步骤6.3)，得到邻域差异图Qr＝{qr(x,y)│x＝1,2,…,m,y＝1,2,…,n}。

步骤7，用谱聚类方法将光谱反射率变化量的模值图Mr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM1。

7.1)将光谱反射率变化模值图Mr按列变成大小为M×1的矩阵XM＝{mr_i│i＝1,2,…,M}，其中M＝m×n，mr_i为光谱反射率变化模值图Mr中第i个模值；

7.2)计算XM中任意两个值间的相似度wm(i,j)，得到模值相似度矩阵WM＝{wm(i,j)│i＝1,2,…,M,j＝1,2,…,M}，其中wm(i,j)的计算公式如下：

$\mathrm{wm}(i,j)=e^{-{({\mathrm{mr}}_i-{\mathrm{mr}}_j)}^2/2};$

其中mr_i为光谱反射率变化模值图Mr中第i个模值，mr_j为光谱反射率变化模值图Mr中第j个模值；

7.3)将模值相似度矩阵WM中的每一列元素相加，得到的结果放在模值相似度矩阵WM相应的对角线上，并令模值相似度矩阵WM中其它位置的元素为0，得到模值对角矩阵DM，再令LM＝DM-WM；

7.4)求LM的前K个递增的特征值对应的特征向量，把K个特征向量放在一起构造一个M×K的矩阵MM，本发明实例中K＝2；

7.5)将MM的每个行向量看成是K维空间中的一点，采用K-means算法进行聚类；

7.6)若MM某一行被归为变化类，则光谱反射率变化模值图Mr相应位置处的像素点也归为变化类，若MM某一行被归为非变化类，则光谱反射率变化模值图Mr相应位置处的像素点也归为非变化类，从而得到初始变化检测结果图CM1。

步骤8，用谱聚类方法将光谱反射率角制图Sr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM2。

8.1)将光谱反射率角制图Sr按列变成大小为M×1的矩阵XS＝{sr_i│i＝1,2,…,M}，其中M＝m×n，sr_i为光谱反射率角制图Sr中第i个光谱角值；

8.2)计算XS中任意两个值间的相似度ws(i,j)，得到光谱角相似度矩阵WS＝{ws(i,j)│i＝1,2,…,M,j＝1,2,…,M}，其中ws(i,j)的计算公式如下：

$\mathrm{ws}(i,j)=e^{-{({\mathrm{sr}}_i-{\mathrm{sr}}_j)}^2/2};$

其中sr_i为光谱反射率角制图Sr中第i个光谱角值，sr_j为光谱反射率角制图Sr中第j个光谱角值；

8.3)将光谱角相似度矩阵WS中的每一列元素相加，得到的结果放在光谱角相似度矩阵WS相应的对角线上，并令光谱角相似度矩阵WS中其它位置的元素为0，得到光谱角对角矩阵DS，再令LS＝DS-WS；

8.4)求LS的前K个递增的特征值对应的特征向量，把K个特征向量放在一起构造一个M×K的矩阵MS，本发明实例中K＝2；

8.5)将MS的每个行向量看成是K维空间中的一点，采用K-means算法进行聚类；

8.6)若MS某一行被归为变化类，则光谱反射率角制图Sr相应位置处的像素点也归为变化类，若MS某一行被归为非变化类，则光谱反射率角制图Sr相应位置处的像素点也归为非变化类，从而得到初始变化检测结果图CM2。

步骤9，用谱聚类方法将邻域差异图Qr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM3。

9.1)将邻域差异图Qr按列变成M×1的矩阵XQ＝{qr_i│i＝1,2,…,M}，其中M＝m×n，qr_i为邻域差异图Qr中第i个邻域差异值；

9.2)计算XQ中任意两个值间的相似度wq(i,j)，得到邻域差相似度矩阵WQ＝{wq(i,j)│i＝1,2,…,M,j＝1,2,…,M}，其中M＝m×n，wq(i,j)的计算公式如下：

$\mathrm{wq}(i,j)=e^{-{({\mathrm{qr}}_i-{\mathrm{qr}}_j)}^2/2};$

其中qr_i为邻域差异图Qr中第i个邻域差异值，qr_j为邻域差异图Qr中第j个邻域差异值；

9.3)将邻域差相似度矩阵WQ中的每一列元素相加，得到的结果放在邻域差相似度矩阵WQ相应的对角线上，并令邻域差相似度矩阵WQ中其它位置的元素为0，得到邻域差对角矩阵DQ，再令LQ＝DQ-WQ；

9.4)求LQ的前K个递增的特征值对应的特征向量，把K个特征向量放在一起构造一个M×K的矩阵MQ，本发明实例中K＝2；

9.5)将MQ的每个行向量看成是K维空间中的一点，采用K-means算法进行聚类；

9.6)若MQ某一行被归为变化类，则把邻域差异图Qr相应位置处的像素点也归为变化类，若MQ某一行被归为非变化类，则把邻域差异图Qr相应位置处的像素点也归为非变化类，从而得到初始变化检测结果图CM3。

步骤10，对CM1、CM2和CM3采用基于邻域概率的融合方法，得到最终变化检测结果图CM。

10.1)计算CM1中像素点(x,y)的k×k邻域内变化类的个数c1(x,y)和非变化类的个数nc1(x,y)，相对应的概率为p_c1(x,y)和p_n1(x,y)；

$p_{c1}(x,y)=\frac{c1(x,y)}{c1(x,y)+\mathrm{nc}1(x,y)};$

$p_{n1}(x,y)=\frac{\mathrm{nc}1(x,y)}{c1(x,y)+\mathrm{nc}1(x,y)};$

10.2)计算CM2中像素点(x,y)的k×k邻域内变化类的个数c2(x,y)和非变化类的个数nc2(x,y)，相应的概率为p_c2(x,y)和p_n2(x,y)；

$p_{c2}(x,y)=\frac{c2(x,y)}{c2(x,y)+\mathrm{nc}2(x,y)};$

$p_{n2}(x,y)=\frac{\mathrm{nc}2(x,y)}{c2(x,y)+\mathrm{nc}2(x,y)};$

10.3)计算CM3中像素点(x,y)的k×k邻域内变化类的个数c3(x,y)和非变化类的个数nc3(x,y)，相应的概率为p_c3(x,y)和p_n3(x,y)；

$p_{c3}(x,y)=\frac{c3(x,y)}{c3(x,y)+\mathrm{nc}3(x,y)};$

$p_{n3}(x,y)=\frac{\mathrm{nc}3(x,y)}{c3(x,y)+\mathrm{nc}3(x,y)};$

10.4)计算像素点(x,y)属于变化类的混合概率cw(x,y)和属于非变化类的混合概率nw(x,y)；

cw(x,y)＝c1(x,y)·p_c1(x,y)+c2(x,y)·p_c2(x,y)+c3(x,y)·p_c3(x,y)；

nw(x,y)＝nc1(x,y)·p_n1(x,y)+nc2(x,y)·p_n2(x,y)+nc3(x,y)·p_n3(x,y)；

10.5)对像素点(x,y)，如果cw(x,y)大于等于nw(x,y)，则该像素点为变化类，否则，该像素点为非变化类；

10.6)对CM1、CM2和CM3中的所有m×n个像素均重复步骤10.1)至步骤10.5)，得到最终变化检测结果图CM。

本发明的效果可通过如下仿真实验结果进一步说明：

1.实验数据

本发明仿真实验所使用的两组数据分别是从1999年11月20日和2001年3月14日拍摄的一幅哥伦比亚比查达省的增强专题制图仪ETM图像集中选取的大小分别为512×547和499×290像素的两个区域，各时相各波段图像的灰度级均为256级。增强专题制图仪图像集的第6、第8波段与其他波段分辨率不同，因此只选取了除第6、第8波段外的第1至第7共6个波段图像。

第一组数据为上述第一个区域的多光谱图像集的两个时相各个波段图像和对应的变化参考图像。如图2所示，其中：

图2(a1)为1999年11月20日拍摄的多光谱图像集中的第1波段图像；

图2(a2)为1999年11月20日拍摄的多光谱图像集中的第2波段图像；

图2(a3)为1999年11月20日拍摄的多光谱图像集中的第3波段图像；

图2(a4)为1999年11月20日拍摄的多光谱图像集中的第4波段图像；

图2(a5)为1999年11月20日拍摄的多光谱图像集中的第5波段图像；

图2(a6)为1999年11月20日拍摄的多光谱图像集中的第7波段图像；

图2(b1)为2001年3月14日拍摄的多光谱图像集中的第1波段图像；

图2(b2)为2001年3月14日拍摄的多光谱图像集中的第2波段图像；

图2(b3)为2001年3月14日拍摄的多光谱图像集中的第3波段图像；

图2(b4)为2001年3月14日拍摄的多光谱图像集中的第4波段图像；

图2(b5)为2001年3月14日拍摄的多光谱图像集中的第5波段图像；

图2(b6)为2001年3月14日拍摄的多光谱图像集中的第7波段图像；

图2(c)为对应的变化参考图像，该图中白色的像素代表变化类，黑色的像素代表非变化类。

第二组数据为上述第二个区域的多光谱图像集的两个时相各个波段图像和对应的变化参考图像。如图3所示，其中：

图3(a1)为1999年11月20日拍摄的多光谱图像集中的第1波段图像；

图3(a2)为1999年11月20日拍摄的多光谱图像集中的第2波段图像；

图3(a3)为1999年11月20日拍摄的多光谱图像集中的第3波段图像；

图3(a4)为1999年11月20日拍摄的多光谱图像集中的第4波段图像；

图3(a5)为1999年11月20日拍摄的多光谱图像集中的第5波段图像；

图3(a6)为1999年11月20日拍摄的多光谱图像集中的第7波段图像；

图3(b1)为2001年3月14日拍摄的多光谱图像集中的第1波段图像；

图3(b2)为2001年3月14日拍摄的多光谱图像集中的第2波段图像；

图3(b3)为2001年3月14日拍摄的多光谱图像集中的第3波段图像；

图3(b4)为2001年3月14日拍摄的多光谱图像集中的第4波段图像；

图3(b5)为2001年3月14日拍摄的多光谱图像集中的第5波段图像；

图3(b6)为2001年3月14日拍摄的多光谱图像集中的第7波段图像；

图3(c)为对应的变化参考图像，该图中白色的像素代表变化类，黑色的像素代表非变化类。

2.对比试验

为了说明本发明的有效性，本发明与如下三个对比方法进行了对比：

对比方法1采用的是变化矢量分析的方法。本发明与对比方法1对比是为了验证本发明方法对抑制由于配准噪声等引起的伪变化信息的有效性。

对比方法2采用的是光谱角制图的方法。本发明与对比方法2对比是为了验证本发明方法在抑制漏检信息方面的优势，具有更加完整地检测变化信息的能力。

对比方法3采用的是本发明提出的邻域差异图方法，邻域差异图方法相比于变化矢量分析方法和光谱角制图方法，其检测变化信息的能力有所提高，但在抑制漏检信息方面表现不突出。本发明与对比方法3对比是为了验证本发明方法能综合利用变化矢量分析方法、光谱反射率角制图方法和邻域差异图方法的优点，有效降低漏检率，提高变化检测准确率。

其中上述各方法中所采用的聚类方法均为谱聚类。

3.实验内容与分析

由于图像中存在的厚云及其阴影会对变化检测的准确性带来影响，因此在利用本发明方法对图2所示的实验数据进行变化检测之前，先对图2和图3两个地区时相一的各波段图像中存在的厚云及其阴影进行了检测和剔除。

实验1，用本发明与对比方法1、对比方法2和对比方法3对第一组多光谱图像集进行变化检测，结果如图4所示，其中：

图4(a)为对比方法1的变化检测结果图；

图4(b)为对比方法2的变化检测结果图；

图4(c)为对比方法3的变化检测结果图；

图4(d)为本发明方法的变化检测结果图。

实验2，用本发明与对比方法1、对比方法2和对比方法3对第二组多光谱图像集进行变化检测，结果如图5所示，其中：

图5(a)为对比方法1的变化检测结果图；

图5(b)为对比方法2的变化检测结果图；

图5(c)为对比方法3的变化检测结果图；

图5(d)为本发明方法的变化检测结果图。

从图4和图5可以看出，对比方法1的检测结果图中除了变化区域，整幅图还分布着面积很小的杂点，存在较多的伪变化区域；对比方法2与对比方法3的检测结果图中的变化区域的整体形状保持较好，但丢失了部分边缘，产生了较多的漏检；本发明与三种对比方法相比，结果图中变化区域的形状保持好，保留了较完整的边缘信息，且伪变化信息较少。

本发明采用虚警数、漏检数、总错误数和正确率四个指标来定量评价变化检测效果的好坏，其中前三个指标越低，变化检测效果越好，最后一个指标越高，变化检测效果越好。

表1列出了本发明方法与三种对比方法对第一组多光谱图像集和第二组多光谱图像集进行变化检测的定量评价结果。

表1

从表1中可以看出，在以上方法的两组检测结果对比中，本发明所提出的对比方法3(邻域差异图法)相比于对比方法1和对比方法2，其检测正确率有所提高，原因是邻域差异图法综合考虑图像像素的光谱信息和邻域信息；本发明方法相对于其他三种方法，其漏检数和总错误数都最低，正确率最高，原因是本发明结合了对比方法1、对比方法2和对比方法3的优点；而本发明的虚警数高于对比方法2和对比方法3，但低于对比方法1。

从整体上看，本发明能够获得比较准确的变化信息，使得检测结果的漏检率大大降低，能保持较完整的变化区域的边缘信息。因此无论是从视觉效果还是从定量评价指标来看，本发明方法都较优。

本实施方式中没有详细叙述的部分属本行业的公知的常用手段，这里不一一叙述。以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (4)

1.一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其特征在于，包括如下步骤：

(1)输入在两个时相获取的同一地区的两个多光谱图像集：I₁＝{A₁ ^b}和I₂＝{A₂ ^b}，其中，I₁为时相一多光谱图像集,I₂为时相二多光谱图像集,A_t ^b为两个多光谱图像集中的每一幅单波段图像，上标b表示波段序号，b＝1,2,…,B，B为总波段数，下标t为时相序号，t＝{1,2}，每一幅单波段图像A_t ^b均由m行n列像素构成；

(2)对两时相多光谱图像集I₁和I₂分别进行维纳滤波去噪，并归一化处理，得到两时相多光谱归一化图像集和

(3)分别对归一化图像集和中每一幅单波段图像，采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到两个时相每一波段的相对地物光谱反射率图像R₁ ^b和R₂ ^b，并将相同时相的每一波段相对地物光谱反射率图像分别构成两个时相的相对地物光谱反射率图像集R₁和R₂；

(4)分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率变化模值mr(x,y)，得到光谱反射率变化模值图Mr＝{mr(x,y)}，其中x为行序号，y为列序号；

(5)分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y)，得到光谱反射率角制图Sr＝{sr(x,y)}；

(6)分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y)，得到邻域差异图Qr＝{qr(x,y)}；

(7)用谱聚类方法将光谱反射率变化模值图Mr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM1；

(8)用谱聚类方法将光谱反射率角制图Sr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM2；

(9)用谱聚类方法将邻域差异图Qr中的像素聚为两类，分别为变化类和非变化类，得到初始变化检测结果图CM3；

(10)对CM1、CM2和CM3采用基于邻域概率的融合方法，得到最终变化检测结果图CM。

2.根据权利要求1所述的一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其特征在于，步骤(5)中计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y)，按如下公式计算：

对两时相各波段的相对地物光谱反射率图像{R₁ ^b}和{R₂ ^b}中的任意空间位置像素点(x,y)，该点的光谱反射率角值sr(x,y)的计算公式如下；

$\mathrm{sr}(x,y)=\mathrm{ar}\cos \left(\frac{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}(r_1^b(x,y)\·r_2^b(x,y))}{\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left(r_1^b(x,y)\right)}^2}\·\sqrt{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\left(r_2^b(x,y)\right)}^2}}\right).$

3.根据权利要求1所述的一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其特征在于，步骤(6)中计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y)，通过如下步骤进行：

(6a)对第t时相第b波段的相对地物光谱反射率图像{R_t ^b}中的任意空间位置像素点(x,y)，计算该点的k×k邻域内的相对地物反射率值的均值u_t ^b(x,y)：

$u_t^b(x,y)=\frac{1}{k\·k}\underset{v=-c}{\overset{c}{\mathrm{\Σ}}}\underset{w=-c}{\overset{c}{\mathrm{\Σ}}}{r}_t^b(x+v,y+w);$

其中c＝(k-1)/2，k的取值为奇数，取k＝3或k＝5，v和w分别表示以中心像素(x,y)为坐标原点时其对应的k×k邻域内像素点的横坐标值与纵坐标值；

(6b)计算两时相第b波段像素点(x,y)的k×k邻域内的邻域变化值g^b(x,y)；

$g^b(x,y)=\frac{1}{k\·k}\underset{v=-c}{\overset{c}{\mathrm{\Σ}}}\underset{w=-c}{\overset{c}{\mathrm{\Σ}}}{(r_1^b(x+v,y+w)-r_2^b(x+v,y+w))}^2;$

(6c)计算像素点(x,y)的邻域差异值qr(x,y)；

$\mathrm{qr}(x,y)=(1-\frac{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\min (r_1^b(x,y)-u_1^b(x,y),r_2^b(x,y)-u_2^b(x,y))}{\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\max (r_1^b(x,y)-u_1^b(x,y),r_2^b(x,y)-u_2^b(x,y))})\·\left(\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\frac{g^b(x,y)}{u_1^b(x,y)}\right).$

4.根据权利要求1所述的一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法，其特征在于，步骤(10)中所述的对CM1、CM2和CM3采用基于邻域概率的融合方法，得到最终变化检测结果图CM，按如下步骤进行：

(10a)计算CM1中像素点(x,y)的k×k邻域内变化类的个数c1(x,y)和非变化类的个数nc1(x,y)，相对应的概率为p_c1(x,y)和p_n1(x,y)；

$p_{c1}(x,y)=\frac{c1(x,y)}{c1(x,y)+\mathrm{nc}1(x,y)};$

$p_{n1}(x,y)=\frac{\mathrm{nc}1(x,y)}{c1(x,y)+\mathrm{nc}1(x,y)};$

(10b)计算CM2中像素点(x,y)的k×k邻域内变化类的个数c2(x,y)和非变化类的个数nc2(x,y)，相应的概率为p_c2(x,y)和p_n2(x,y)；

$p_{c2}(x,y)=\frac{c2(x,y)}{c2(x,y)+\mathrm{nc}2(x,y)};$

$p_{n2}(x,y)=\frac{\mathrm{nc}2(x,y)}{c2(x,y)+\mathrm{nc}2(x,y)};$

(10c)计算CM3中像素点(x,y)的k×k邻域内变化类的个数c3(x,y)和非变化类的个数nc3(x,y)，相应的概率为p_c3(x,y)和p_n3(x,y)；

$p_{c3}(x,y)=\frac{c3(x,y)}{c3(x,y)+\mathrm{nc}3(x,y)};$

$p_{n3}(x,y)=\frac{\mathrm{nc}3(x,y)}{c3(x,y)+\mathrm{nc}3(x,y)};$

(10d)计算像素点(x,y)属于变化类的混合概率cw(x,y)和属于非变化类的混合概率nw(x,y)；

cw(x,y)＝c1(x,y)·p_c1(x,y)+c2(x,y)·p_c2(x,y)+c3(x,y)·p_c3(x,y)；

nw(x,y)＝nc1(x,y)·p_n1(x,y)+nc2(x,y)·p_n2(x,y)+nc3(x,y)·p_n3(x,y)；

(10e)对像素点(x,y)，如果cw(x,y)大于等于nw(x,y)，则该像素点为变化类，否则，该像素点为非变化类；

(10f)对CM1、CM2和CM3中的所有m×n个像素均重复步骤10a)至步骤10e)，得到最终变化检测结果图CM。