CN103226832B - 基于光谱反射率变化分析的多光谱遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于光谱反射率变化分析的多光谱遥感影像变化检测方法，主要解决现有技术对不同时相图像间的整体亮度差异敏感的问题。其实现步骤为：（1）输入同一地区的已配准的两时相多光谱图像集，对其进行维纳滤波去噪，并归一化处理；（2）将处理后的图像集转换为相对地物光谱反射率图像集；（3）计算光谱反射率变化量的方差和模，得到光谱反射率变化量的方差图和模值图，并分别进行增强，得到增强方差图和增强模值图；（4）分别对增强方差图和增强模值图进行分割，并融合分割所得的二值图，得到变化检测结果图。本发明无需人工参与，检测精确度高，可用于土地利用与覆盖监测、植被覆盖监测、水资源及矿产资源监测。

Description

基于光谱反射率变化分析的多光谱遥感影像变化检测方法

技术领域

本发明属于光学遥感图像处理技术领域，涉及两时相多光谱遥感图像间存在较大的整体亮度差异情况下的无监督变化检测，可用于土地利用、植被覆盖、水资源以及矿产资源等方面的变化监测。

背景技术

随着遥感图像数据的公开化程度和处理技术的不断提高，利用多时相遥感图像数据进行土地利用、植被覆盖、水资源以及矿产资源等方面的变化监测已经越来越普及，这在环境、农林、水利和矿产等诸多方面具有重要意义。由于不同地物类型的变化可能会反映在不同的光谱范围内，而多光谱遥感图像数据具有从可见光到红外光波段的多个接收频段，丰富的光谱信息增加了识别多种类型变化的可能性和可信度，因此多光谱遥感图像被广泛用于地物的变化检测。

一幅完整的多光谱图像具有多个不同波段的图像。目前人们在对其进行变化检测时往往利用其中某一个波段或者多个波段的组合。其中，利用单一波段的变化检测算法通常是在某一波段上构造差异图像并进行分类，故对该类方法的研究通常是对差异图像的增强方法、自适应阈值或其他分割方法、聚类方法等方面的改进。该类方法原理简单，就单一波段的变化检测精确度而言能够达到很高。但由于不同地物类型的变化可能会反映在不同的光谱范围内，因此这类方法在进行变化检测前需要先根据待检测变化的类型和以往经验选择适合的波段，而对于复杂地物这样的选择并不总是准确、有效；如果待检测的变化信息分布在不同波段上，则无论选择哪一个波段都无法完全检测到所有变化信息。因此，有效利用多个波段图像的变化信息是解决多光谱变化检测问题的唯一有效途径。

目前，利用多个波段光谱图像进行变化检测的方法主要有变化矢量分析CVA、主成分分析PCA和独立成分分析ICA。标准变化矢量分析先将两幅不同时相的多光谱图像相同波段对应位置处的像素作差得到光谱变化矢量SCVs，在计算其模值后通过阈值分割得到变化类与非变化类的分类结果。由于同时考虑了多个波段图像的变化信息，变化矢量分析相对于单一波段的变化检测增强了检测结果的可靠性和检测多类型变化的能力。然而，变化矢量分析在计算光谱变化矢量的模值时只是简单地将各波段差异信息叠加，忽略了波段之间的差异性和相关性；同时由于易受到不同时相图像间整体亮度差异的影响而使得光谱变化矢量模值图的背景亮度与变化类目标的亮度非常接近，大大降低了阈值分割的准确性和稳定性。尽管有很多学者对变化矢量分析进行了进一步研究，但由于受到其整体框架及思想的局限，都只是对光谱变化矢量方向的描述、分割阈值的确定以及表示空间等方面的改进，上述缺陷并没有得到改善。相对地，主成分分析用于多光谱图像变化检测时可以减小波段间的冗余信息和增强差异信息，通常采用两种方式进行——直接对两时相图像做主成分分析和对各波段差异图像做主成分分析。然而，少数波段的变化信息经过主成分分析变换后多分布在次主成分图像上，容易造成变化信息的丢失。另外，主成分分析与变化矢量分析相似，容易受到不同时相图像间整体亮度差异的影响而使得检测准确率大幅下降，即使先采用辐射校正的方法进行改善，仍不能从根本上解决这一问题。与主成分分析非常相似，独立成分分析用于变化检测时也通常采用直接对两时相图像作独立成分分析或对各波段差异图像作独立成分分析两种策略。独立成分分析与主成分分析相比虽然对配准误差和亮度差异具有一定鲁棒性，但独立成分分析变换使不同类型的变化区域分布在若干相互独立的分量图像上，通常需要依靠经验和先验知识人工选择某几个分量图像进行后续分类，造成完全检测变化信息的困难。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于光谱反射率变化分析的多光谱遥感图像无监督变化检测方法，以避免不同时相图像间整体亮度差异的影响及人工选择光谱波段造成变化信息检测的不完整，实现变化检测的完全自动化。

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

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

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

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

（4）分别计算两个时相的相对地物光谱反射率图像集R₁和R₂中各波段对应空间位置像素(x,y)的光谱反射率变化量Dr(x,y)、光谱反射率变化量的方差vr(x,y)和光谱反射率变化量的模值mr(x,y)，得到光谱反射率变化量的方差图Vr={vr(x,y)|x=1,2,…,m,y=1,2,…,n}和光谱反射率变化量的模值图Mr={mr(x,y)|x=1,2,…,m,y=1,2,…,n}，其中x为列序号，y为行序号；

（5）对光谱反射率变化量的模值图Mr做增强处理，得到增强模值图AMr：

（5a）将时相1归一化图像集的各波段归一化图像与时相2归一化图像集对应波段的归一化图像中任意一点(x,y)的像素值作差，得到该像素点的光谱变化量S(x,y)={dA^b(x,y)|b=1,2,…,B}，其中

（5b）计算像素点(x,y)的光谱变化矢量S(x,y)的模M(x,y)：

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

（5c）对两时相归一化图像集和中的所有m×n个像素均重复步骤（5a）和（5b），得到光谱变化矢量的模值图Md={M(x,y)|x=1,2,…,m,y=1,2,…,n}；

（5d）将光谱变化矢量的模值图Md与光谱反射率变化量的模值图Mr对应空间位置的像素值相乘，得到增强模值图AMr={AMr(x,y)|x=1,2,…,m,y=1,2,…,n}，其中AMr(x,y)为增强模值图中(x,y)点的像素值，AMr(x,y)=M(x,y)×mr(x,y)；

（6）对光谱反射率变化量的方差图Vr做增强处理，得到增强方差图AVr：

（6a）将时相1归一化图像集的第b波段归一化图像与时相2归一化图像集的第b波段归一化图像对应空间位置的像素值作差并取绝对值，得到该波段的差值差异图DI^b={|dA^b(x,y)|│x=1,2,…,m,y=1,2,…,n}，其中|·|表示取绝对值操作；

（6b）对所有B个波段均重复步骤（6a），得到所有B个波段的差值差异图DI¹、DI²、…、DI^b、…、DI^B；

（6c）将所有B个波段的差值差异图DI¹、DI²、…、DI^b、…、DI^B以及增强模值图AMr均采用最大类间方差法进行阈值分割，分别得到B个波段的差值差异图的二值分割图DS¹、DS²、…、DS^b、…和DS^B和增强模值图的二值分割图AMs；

（6d）计算各波段差值差异图的二值分割图DS^b与增强模值图的二值分割图AMs的差异度Sim^b：

${\mathrm{Sim}}^b=\underset{y=1}{\overset{n}{\mathrm{\Σ}}}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}|{\mathrm{DS}}^b(x,y)-\mathrm{AMs}(x,y)|,$

其中，DS^b(x,y)为差值差异图的二值分割图DS^b中(x,y)点的像素值，AMs(x,y)为增强模值图的二值分割图AMs中(x,y)点的像素值；

（6e）将所有波段的差异度{Sim^b}中最小值所对应的波段序号记为a，计算第a个波段的差值差异图DI^a与光谱反射率变化量的方差图Vr在(x,y)点的像素值的几何平均值，得到一幅增强方差图AVr={AVr(x,y)|x=1,2,…,m,y=1,2,…,n}，其中，AVr(x,y)为增强方差图在点(x,y)的像素值， $\mathrm{AVr}(x,y)=\sqrt{\mathrm{vr}(x,y)\×{\mathrm{DI}}^a(x,y)};$

（7）对增强方差图AVr采用最大类间方差法进行阈值分割，得到增强方差图的二值分割图AVs；

（8）将增强方差图的二值分割图AVs和增强模值图的二值分割图AMs对应空间位置像素点的值逐一做与运算，得到差模交集图VM；

（9）对差模交集图VM中的变化类像素进行补充检测，得到光谱反射率变化粗检图Cm；

（10）对增强方差图AVr与增强模值图AMr对应空间位置像素点的值计算几何平均值，得到差模平均图VMA；

（11）将光谱反射率变化粗检图Cm中所有像素值为1的像素点作为变化类种子点，在差模平均图VMA中利用区域生长法对变化类种子点进行区域生长，得到变化检测结果图CM。

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

a、本发明与现有的多光谱遥感图像变化检测方法相比，提出并采用了光谱反射率变化量的方差和模作为衡量地物变化与否的两个新指标，提高了检测多种类型变化的能力；另外，本发明采用图像灰度信息对光谱反射率变化量的方差图和模值图进行增强处理，消除了不同时相图像间存在的整体亮度差异对检测稳定性和准确性的影响。

b、本发明与现有基于独立成分分析ICA或主成分分析PCA的多光谱遥感图像变化检测方法相比，由于不产生多个成分图像，所有变化信息都存在于光谱反射率变化量的方差图和光谱反射率变化量的模值图中，因此不存在需要人工选取某成分图像的问题，实现了完全自动化。

附图说明

图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均由n行m列像素构成。

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

2.1）将两时相多光谱图像集I₁和I₂中每一幅图像A_t ^b的灰度值区间从0～255转换到0～1之间，然后采用窗口大小为3×3像素的维纳滤波进行去噪，得到相应的去噪后图像

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）用时相1和时相2的归一化图像分别构成归一化图像集和即

${\hat{I}}_1=\{{{\hat{A}}_1}^1,{{\hat{A}}_1}^2,\·\·\·,{{\hat{A}}_1}^b,\·\·\·{{\hat{A}}_1}^B\},$

${\hat{I}}_2=\{{{\hat{A}}_2}^1,{{\hat{A}}_2}^2,\·\·\·,{{\hat{A}}_2}^b,\·\·\·{{\hat{A}}_2}^B\}.$

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

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

$\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)$

其中，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）对时相1所有B个波段的归一化图像重复步骤3.1）至步骤3.3），得到时相1的B个波段的相对地物光谱反射率图像{R₁ ^b}，构成时相1的相对地物光谱反射率图像集R₁={R₁ ^b|b=1,2,…,B}；

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

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

4.1）将时相1各波段的相对地物光谱反射率图像{R₁ ^b}与时相2对应波段的相对地物光谱反射率图像{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)的光谱反射率变化量Dr(x,y)的方差vr(x,y)：

$\mathrm{vr}(x,y)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{[{\mathrm{dr}}^b(x,y)-\mathrm{ar}(x,y)]}^2;$

其中，ar(x,y)为像素点(x,y)的所有B个波段光谱反射率变化量{dr^b(x,y)}的均值：

$\mathrm{ar}(x,y)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{dr}}^b(x,y);$

4.3）计算像素点(x,y)的光谱反射率变化量Dr(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.4）对两时相的相对地物光谱反射率图像集R₁和R₂中的所有m×n个像素均重复步骤4.1）和步骤4.2），得到光谱反射率变化量的方差图Vr={vr(x,y)|x=1,2,…,m,y=1,2,…,n}；

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

步骤5，对光谱反射率变化量的模值图Mr做增强处理，得到增强模值图AMr。

5.1）将时相1归一化图像集的各波段归一化图像与时相2归一化图像集对应波段的归一化图像中任意一点(x,y)的像素值作差，得到该像素点的光谱变化量S(x,y)={dA^b(x,y)|b=1,2,…,B}，其中

5.2）计算像素点(x,y)的光谱变化矢量S(x,y)的模M(x,y)：

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

5.3）对两时相所有B个波段的归一化图像集和中的所有m×n个像素均重复步骤5.1）和5.2），得到光谱变化矢量的模值图Md={M(x,y)|x=1,2,…,m,y=1,2,…,n}；

5.4）将光谱变化矢量的模值图Md与光谱反射率变化量的模值图Mr对应空间位置的像素值相乘，得到增强模值图AMr={AMr(x,y)|x=1,2,…,m,y=1,2,…,n}，其中AMr(x,y)为增强模值图中(x,y)点的像素值，AMr(x,y)=M(x,y)×mr(x,y)。

步骤6，对光谱反射率变化量的方差图Vr做增强处理，得到增强方差图AVr。

6.1）将时相1归一化图像集第b波段归一化图像与时相2归一化图像集第b波段的归一化图像对应空间位置的像素值作差并取绝对值，得到该波段的差值差异图DI^b={|dA^b(x,y)|│x=1,2,…,m,y=1,2,…,n}，其中|·|表示取绝对值操作；

6.2）对所有B个波段均重复步骤6.1），得到所有B个波段的差值差异图DI¹、DI²、…、DI^b、…、DI^B；

6.3）将所有B个波段的差值差异图DI¹、DI²、…、DI^b、…、DI^B以及增强模值图AMr均采用最大类间方差法进行阈值分割，分别得到B个波段的差值差异图的二值分割图DS¹、DS²、…、DS^b、…和DS^B和增强模值图的二值分割图AMs；

6.4）计算各个波段的差值差异图的二值分割图DS^b与增强模值图的二值分割图AMs的差异度Sim^b；

${\mathrm{Sim}}^b=\underset{y=1}{\overset{n}{\mathrm{\Σ}}}\underset{x=1}{\overset{m}{\mathrm{\Σ}}}|{\mathrm{DS}}^b(x,y)-\mathrm{AMs}(x,y)|,$

其中，DS^b(x,y)为差值差异图的二值分割图DS^b中(x,y)点的像素值，AMs(x,y)为增强模值图的二值分割图AMs中(x,y)点的像素值；

6.5）将所有波段的差异度{Sim^b}中最小值所对应的波段序号记为a，计算第a个波段的差值差异图DI^a与光谱反射率变化量的方差图Vr在(x,y)点的像素值的几何平均值，得到一幅增强方差图AVr={AVr(x,y)|x=1,2,…,m,y=1,2,…,n}，其中，AVr(x,y)为增强方差图在点(x,y)的像素值：

$\mathrm{AVr}(x,y)=\sqrt{\mathrm{Vr}(x,y)\×{\mathrm{DI}}^a(x,y)}.$

步骤7，对增强方差图AVr采用最大类间方差法进行阈值分割，得到增强方差图的二值分割图AVs。

步骤8，将增强方差图的二值分割图AVs和增强模值图的二值分割图AMs对应空间位置像素点的值逐一进行与运算，得到差模交集图VM。

步骤9，对差模交集图VM中的变化类像素进行补充检测，得到光谱反射率变化粗检图Cm。

9.1）将差模交集图VM中第1行第1列的像素点(1,1)作为当前像素，记为(x_z,y_z)；

9.2）判断当前像素(x_z,y_z)的像素值是否为1，如果为1，则执行步骤9.4）；否则执行步骤9.3）；

9.3）对当前像素(x_z,y_z)的行、列序号分别与图像总行数和总列数进行比较：

若当前像素(x_z,y_z)的列序号x_z<m时，则行序号y_z保持不变，将列序号x_z加1，返回步骤9.2），其中m为图像总列数；

若当前像素的列序号x_z=m且行序号y_z<n时，则对行序号y_z加1，且将列序号x_z置为1，返回步骤9.2），其中n为图像总行数；

若当前像素的列序号x_z=m且行序号y_z=n时，得到光谱反射率变化粗检图Cm；

9.4）在增强方差图的二值分割图AVs和增强模值图的二值分割图AMs中，分别选取列序号范围为 行序号范围为 的两个K×K大小的图像块，将这两个图像块中所有像素值为1的像素点在差模交集图VM中对应位置像素的值赋为1，返回步骤（9c），其中图像块尺寸K的取值范围为7或9，本发明实例中K=9。

步骤10，计算增强方差图AVr在(x,y)点的像素值AVr(x,y)与增强模值图AMr对应空间位置的像素值AMr(x,y)的几何平均值得到一幅差模平均图VMA={VMA(x,y)│x=1,2,…,m,y=1,2,…,n}。

步骤11，将光谱反射率变化粗检图Cm中所有像素值为1的像素点作为变化类种子点，在差模平均图VMA中利用区域生长法对变化类种子点进行区域生长，得到变化检测结果图CM。

11.1）在光谱反射率变化粗检图Cm中所有像素值为1的像素点中任意选一个像素点作为变化类种子点；

11.2）对已选中的变化类种子点在差模平均图VMA中空间位置对应的像素点，计算以该像素点为中心的8-连通的邻域像素点与该中心像素点的像素值之差的绝对值，对于绝对值小于邻域判断阈值Th的邻域像素点，如果该邻域像素点在光谱反射率变化粗检图Cm中对应的像素点不是变化类种子点，则将该邻域像素点作为一个新的变化类种子点加入光谱反射率变化粗检图中；否则返回步骤11.1）；其中，区域生长的邻域判断阈值Th的取值范围为[0.08,0.1]，本实例中取Th=0.09；

11.3）对光谱反射率变化粗检图中的每一个变化类种子点和新加入的变化类种子点重复步骤11.1）-步骤11.2），直到不再有新的变化类种子点加入光谱反射率变化粗检图为止，得到变化检测结果图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采用的是变化矢量分析CVA的方法。本发明与对比方法1对比是为了验证本发明方法对抑制不同时相图像间整体亮度差异的干扰的有效性。

对比方法2采用的是主成分分析PCA的方法，利用主成分分析PCA得到的第一主成分图像并利用最大类间方差法OTSU进行阈值分割得到变化检测结果图。本发明与对比方法2对比是为了验证本发明方法对不同时相图像间整体亮度差异的干扰具有更好的抑制效果，以及具有更加完整地检测变化信息的能力。

对比方法3采用的是独立成分分析ICA的方法，对6个波段的差异图像先利用独立成分分析ICA得到6个成分图像，再由人工选取其中某一个成分图像利用最大类间方差法OTSU进行阈值分割得到变化检测结果图。本发明与对比方法3对比是为了验证本发明方法具有更加完整地检测变化信息的能力，并且是完全自动化的。

其中上述各方法中所采用的自适应阈值分割方法均采用最大类间方差法，旨在排除因采用不同阈值分割方法对检测结果产生的影响。

3.实验内容与分析

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

实验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中可以看出，在以上四种方法的两组检测结果对比中，本发明方法相对于其他三种方法，其虚警数和虚警率都最低，正确率最高；而本发明的漏检数和漏检率略高于变化矢量分析CVA和主成分分析PCA，但低于独立成分分析ICA。从整体上看，本发明方法能够有效地抑制不同时相图像间整体亮度差异的干扰，使得检测结果的虚警率大大降低，同时又能基本检测到所有的变化信息，变化区域的形状和边缘保持较完整。因此无论是从视觉效果还是从定量评价指标来看，本发明方法都相比这三种方法具有一定的提升。

Claims (2)

1.一种基于光谱反射率变化分析的多光谱遥感影像变化检测方法，包括步骤如下：

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

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

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

(4)计算光谱反射率变化量的方差图Vr和光谱反射率变化量的模值图Mr：

(4a)将时相1第b波段的相对地物光谱反射率图像R₁ ^b与时相2第b波段的相对地物光谱反射率图像R₂ ^b中任意一点(x,y)的像素值作差，得到第b波段的反射率差值dr^b(x,y)：

dr^b(x,y)＝r₁ ^b(x,y)–r₂ ^b(x,y)；

(4b)对所有B个波段均重复步骤(4a)，得到B个波段的反射率差值，并将它们共同构成该像素点的光谱反射率变化量Dr(x,y)＝{dr^b(x,y)|b＝1,2,…,B}；

(4c)通过下式计算光谱反射率变化量的方差：

$\mathrm{vr}(x,y)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{[{\mathrm{dr}}^b(x,y)-\mathrm{ar}(x,y)]}^2,$

其中，ar(x,y)为像素点(x,y)的所有B个波段光谱反射率变化量{dr^b(x,y)}的均值： $\mathrm{ar}(x,y)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\mathrm{dr}}^b(x,y);$

(4d)通过如下公式计算光谱反射率变化量的模值；

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

(4e)根据光谱反射率变化量的方差vr(x,y)和光谱反射率变化量的模值mr(x,y)，得到光谱反射率变化量的方差图Vr＝{vr(x,y)|x＝1,2,…,m,y＝1,2,…,n}和光谱反射率变化量的模值图Mr＝{mr(x,y)|x＝1,2,…,m,y＝1,2,…,n}，其中x为列序号，y为行序号；

(5)对光谱反射率变化量的模值图Mr做增强处理，得到增强模值图AMr：

(5a)将时相1归一化图像集的各波段归一化图像与时相2归一化图像集对应波段的归一化图像中任意一点(x,y)的像素值作差，得到该像素点的光谱变化量S(x,y)＝{dA^b(x,y)|b＝1,2,…,B}，其中

(5b)计算像素点(x,y)的光谱变化矢量S(x,y)的模M(x,y)：

$M(x,y)=\sqrt{\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\left({\mathrm{dA}}^b(x,y)\right)}^2};$

(5c)对两时相归一化图像集和中的所有m×n个像素均重复步骤(5a)和(5b)，得到光谱变化矢量的模值图Md＝{M(x,y)|x＝1,2,…,m,y＝1,2,…,n}；

(5d)将光谱变化矢量的模值图Md与光谱反射率变化量的模值图Mr对应空间位置的像素值相乘，得到增强模值图AMr＝{AMr(x,y)|x＝1,2,…,m,y＝1,2,…,n}，其中AMr(x,y)为增强模值图中(x,y)点的像素值，AMr(x,y)＝M(x,y)×mr(x,y)；

(6)对光谱反射率变化量的方差图Vr做增强处理，得到增强方差图AVr：

(6a)将时相1归一化图像集的第b波段归一化图像与时相2归一化图像集的第b波段归一化图像对应空间位置的像素值作差并取绝对值，得到该波段的差值差异图DI^b＝{|dA^b(x,y)|│x＝1,2,…,m,y＝1,2,…,n}，其中|·|表示取绝对值操作；

(6b)对所有B个波段均重复步骤(6a)，得到所有B个波段的差值差异图DI¹、DI²、…、DI^b、…、DI^B；

(6c)将所有B个波段的差值差异图DI¹、DI²、…、DI^b、…、DI^B以及增强模值图AMr均采用最大类间方差法进行阈值分割，分别得到B个波段的差值差异图的二值分割图DS¹、DS²、…、DS^b、…和DS^B和增强模值图的二值分割图AMs；

(6d)计算各波段差值差异图的二值分割图DS^b与增强模值图的二值分割图AMs的差异度Sim^b：

${\mathrm{Sim}}^b=\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{m}{\mathrm{\&Sigma;}}}|{\mathrm{DS}}^b(x,y)-\mathrm{AMs}(x,y)|,$

其中，DS^b(x,y)为差值差异图的二值分割图DS^b中(x,y)点的像素值，AMs(x,y)为增强模值图的二值分割图AMs中(x,y)点的像素值；

(6e)将所有波段的差异度{Sim^b}中最小值所对应的波段序号记为a，计算第a个波段的差值差异图DI^a与光谱反射率变化量的方差图Vr在(x,y)点的像素值的几何平均值，得到一幅增强方差图AVr＝{AVr(x,y)|x＝1,2,…,m,y＝1,2,…,n}，其中，AVr(x,y)为增强方差图在点(x,y)的像素值， $\mathrm{AVr}(x,y)=\sqrt{\mathrm{vr}(x,y)\&times;{\mathrm{DI}}^a(x,y)};$

(7)对增强方差图AVr采用最大类间方差法进行阈值分割，得到增强方差图的二值分割图AVs；

(8)将增强方差图的二值分割图AVs和增强模值图的二值分割图AMs对应空间位置像素点的值逐一做与运算，得到差模交集图VM；

(9)对差模交集图VM中的变化类像素进行补充检测，得到光谱反射率变化粗检图Cm:

(9a)将差模交集图VM中第1行第1列的像素点(1,1)作为当前像素，记为(x_z,y_z)；

(9b)判断当前像素(x_z,y_z)的像素值是否为1，如果为1，则执行步骤(9d)；否则执行步骤(9c)；

(9c)将当前像素(x_z,y_z)的行、列序号分别与图像总行数和总列数进行比较：

若当前像素(x_z,y_z)的列序号x_z<m时，则行序号y_z保持不变，将列序号x_z加1，返回步骤(9b)，其中m为图像总列数；

若当前像素的列序号x_z＝m且行序号y_z<n时，则对行序号y_z加1，且将列序号x_z置为1，返回步骤(9b)，其中n为图像总行数；

若当前像素的列序号x_z＝m且行序号y_z＝n时，得到光谱反射率变化粗检图Cm；

(9d)在增强方差图的二值分割图AVs和增强模值图的二值分割图AMs中，分别选取列序号范围为行序号范围为 的两个K×K大小的图像块，将这两个图像块中所有像素值为1的像素点在差模交集图VM中对应位置像素的值赋为1，返回步骤(9c)，其中图像块尺寸K的取值范围为7或9；

(10)对增强方差图AVr与增强模值图AMr对应空间位置像素点的值计算几何平均值，得到差模平均图VMA；

(11)将光谱反射率变化粗检图Cm中所有像素值为1的像素点作为变化类种子点，在差模平均图VMA中利用区域生长法对变化类种子点进行区域生长，得到变化检测结果图CM。

2.根据权利要求1所述的基于光谱反射率变化分析的多光谱遥感影像变化检测方法，其中所述步骤(11)中将光谱反射率变化粗检图Cm中所有像素值为1的像素点作为变化类种子点，在差模平均图VMA中利用区域生长法对变化类种子点进行区域生长，按如下步骤进行：

(11a)在光谱反射率变化粗检图Cm中所有像素值为1的像素点中任意选一个像素点作为变化类种子点；

(11b)对已选中的变化类种子点在差模平均图VMA中空间位置对应的像素点，计算以该像素点为中心的8-连通的邻域像素点与该中心像素点的像素值之差的绝对值，对于绝对值小于邻域判断阈值Th的邻域像素点，Th的取值范围为[0.08,0.1]，如果该邻域像素点在光谱反射率变化粗检图Cm中对应的像素点不是变化类种子点，则将该邻域像素点作为一个新的变化类种子点加入光谱反射率变化粗检图中；

(11c)对光谱反射率变化粗检图中的每一个变化类种子点和新加入的变化类种子点重复步骤(11a)-步骤(11b)，直到不再有新的变化类种子点加入光谱反射率变化粗检图为止，得到变化检测结果图CM。