CN104463881B - 一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 - Google Patents
一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 Download PDFInfo
- Publication number
- CN104463881B CN104463881B CN201410770398.7A CN201410770398A CN104463881B CN 104463881 B CN104463881 B CN 104463881B CN 201410770398 A CN201410770398 A CN 201410770398A CN 104463881 B CN104463881 B CN 104463881B
- Authority
- CN
- China
- Prior art keywords
- mrow
- neighborhood
- msubsup
- variation
- spectral
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 230000003595 spectral effect Effects 0.000 title claims abstract description 139
- 230000008859 change Effects 0.000 title claims abstract description 116
- 238000000034 method Methods 0.000 title claims abstract description 73
- 238000002310 reflectometry Methods 0.000 title claims abstract description 68
- 238000001228 spectrum Methods 0.000 title claims abstract description 18
- 230000004927 fusion Effects 0.000 title claims abstract description 9
- 238000001514 detection method Methods 0.000 claims abstract description 75
- 238000012545 processing Methods 0.000 claims abstract description 6
- 238000001914 filtration Methods 0.000 claims abstract description 5
- 239000011159 matrix material Substances 0.000 claims description 42
- 239000013598 vector Substances 0.000 claims description 14
- 238000004422 calculation algorithm Methods 0.000 claims description 6
- 238000007500 overflow downdraw method Methods 0.000 claims description 5
- 238000010586 diagram Methods 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 101000986989 Naja kaouthia Acidic phospholipase A2 CM-II Proteins 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 3
- 239000008186 active pharmaceutical agent Substances 0.000 claims description 2
- 230000002123 temporal effect Effects 0.000 claims description 2
- 238000012544 monitoring process Methods 0.000 abstract description 4
- 238000012360 testing method Methods 0.000 abstract description 2
- 230000008878 coupling Effects 0.000 abstract 1
- 238000010168 coupling process Methods 0.000 abstract 1
- 238000005859 coupling reaction Methods 0.000 abstract 1
- 238000012733 comparative method Methods 0.000 description 11
- 238000004458 analytical method Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 238000013507 mapping Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000011158 quantitative evaluation Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10036—Multispectral image; Hyperspectral image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30168—Image quality inspection
Landscapes
- Engineering & Computer Science (AREA)
- Quality & Reliability (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法,其实现步骤主要包括:1)输入同一地区的已配准的两时相多光谱图像集,对其进行维纳滤波去噪并归一化处理;2)将处理后的图像集转换为相对地物光谱反射率图像集;3)计算光谱反射率变化模值,得到光谱反射率变化模值图;4)计算光谱反射率角值,得到光谱反射率角制图;5)计算光谱反射率的邻域差异值,得到邻域差异图;6)分别对光谱反射率变化模值图、光谱反射率角制图和邻域差异图进行聚类,得到二值图;7)对所得二值图进行基于邻域概率的融合,得到变化检测结果图。本发明无需人工参与,检测精确度高,可用于城区扩展监测、森林和植被变化监测等领域。
Description
技术领域
本发明属于光学遥感图像处理技术领域,具体为一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法,可用于土地利用、植被覆盖等方面的变化监测。
背景技术
遥感图像的变化检测是通过分析和提取同一地区不同时相的遥感图像之间存在的电磁波谱特征差异或空间结构特征差异,来识别物体的状态变化或现象变化的过程。在国民经济和国防建设的诸多领域已得到广泛应用,如农业调查、森林和植被变化监测等。由于不同地物类型的变化可能会反映在不同的光谱范围内,而多光谱遥感图像数据具有从可见光到红外光波段的多个接收频段,丰富的光谱信息增加了识别多种类型变化的可能性和可信度,因此多光谱遥感图像被广泛用于地物的变化检测。
多光谱遥感图像变化检测最为常见的是基于变化矢量分析(CVA)的变化检测,一般是先通过CVA方法获取差异图,然后对差异图进行变化与未变化的分类。由于同时考虑了多个波段图像的变化信息,CVA相对于单一波段的变化检测增强了检测结果的可靠性。然而,CVA在计算光谱变化矢量的模值时只是简单地将各波段差异信息叠加,由于配准误差等原因,对于没有发生变化的区域,两时相图像对应位置的灰度值也存在一定的偏差,故而简单地对CVA方法获取的差异图像进行分类会导致变化检测结果中存在很多的伪变化信息。在变化检测中引入光谱角制图(SAM)可以降低伪变化信息,SAM方法通过计算两时相光谱矢量之间的角度来定义两者的差异性,然后对两时相差异图进行变化与未变化的分类。SAM方法利用了角度这一参数,减少了光照对不同时相图像间整体亮度差异的影响,从而提高了检测准确率。然而,该方法用于多光谱图像时,由于其波段信息较少,导致不同时相图像间对应地物光谱矢量之间的差异不够大,因此对SAM方法获取的差异图像进行分类所得变化检测结果中,虽然降低了伪变化信息,但同时引入了很多的漏检。
上述的变化检测方法都利用了图像的光谱信息,没有考虑像素的邻域信息,使得变化检测结果不理想。
发明内容
本发明的目的在于针对上述已有技术的不足,提出一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法,综合考虑两时相图像像素的光谱信息和邻域信息,并在此基础上结合上述技术,从而提高检测结果的正确率。
为实现上述目的,本发明提出了一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法,其技术方案是:
1、一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法,其特征在于,包括如下步骤:
(1)输入在两个时相获取的同一地区的两个多光谱图像集:I1={A1 b}和I2={A2 b},其中,I1为时相一多光谱图像集,I2为时相二多光谱图像集,At b为两个多光谱图像集中的每一幅单波段图像,上标b表示波段序号,b=1,2,…,B,B为总波段数,下标t为时相序号,t={1,2},每一幅单波段图像At b均由m行n列像素构成;
(2)对两时相多光谱图像集I1和I2分别进行维纳滤波去噪,并归一化处理,得到两时相多光谱归一化图像集和
(3)分别对归一化图像集和中每一幅单波段图像,采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值,得到两个时相每一波段的相对地物光谱反射率图像R1 b和R2 b,并将相同时相的每一波段相对地物光谱反射率图像分别构成两个时相的相对地物光谱反射率图像集R1和R2;
(4)分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率变化模值mr(x,y),得到光谱反射率变化模值图Mr={mr(x,y)},其中x为行序号,y为列序号;
(5)分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y),得到光谱反射率角制图Sr={sr(x,y)};
(6)分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(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),按如下公式计算:
对两时相各波段的相对地物光谱反射率图像{R1 b}和{R2 b}中的任意空间位置像素点(x,y),该点的光谱反射率角值sr(x,y)的计算公式如下;
上述步骤(6)中计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y),通过如下步骤进行:
(6a)对第t时相第b波段的相对地物光谱反射率图像{Rt b}中的任意空间位置像素点(x,y),计算该点的k×k邻域内的相对地物反射率值的均值ut b(x,y):
其中c=(k-1)/2,k的取值为奇数,取k=3或k=5,v和w分别表示以中心像素(x,y)为坐标原点时其对应的k×k邻域内像素点的横坐标值与纵坐标值;
(6b)计算两时相第b波段像素点(x,y)的k×k邻域内的邻域变化值gb(x,y);
(6c)计算像素点(x,y)的邻域差异值qr(x,y);
上述步骤(10)中所述的对CM1、CM2和CM3采用基于邻域概率的融合方法,得到最终变化检测结果图CM,按如下步骤进行:
(10a)计算CM1中像素点(x,y)的k×k邻域内变化类的个数c1(x,y)和非变化类的个数nc1(x,y),相对应的概率为pc1(x,y)和pn1(x,y);
(10b)计算CM2中像素点(x,y)的k×k邻域内变化类的个数c2(x,y)和非变化类的个数nc2(x,y),相应的概率为pc2(x,y)和pn2(x,y);
(10c)计算CM3中像素点(x,y)的k×k邻域内变化类的个数c3(x,y)和非变化类的个数nc3(x,y),相应的概率为pc3(x,y)和pn3(x,y);
(10d)计算像素点(x,y)属于变化类的混合概率cw(x,y)和属于非变化类的混合概率nw(x,y);
cw(x,y)=c1(x,y)·pc1(x,y)+c2(x,y)·pc2(x,y)+c3(x,y)·pc3(x,y);
nw(x,y)=nc1(x,y)·pn1(x,y)+nc2(x,y)·pn2(x,y)+nc3(x,y)·pn3(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,输入在两个时相获取的同一地区的两个多光谱图像集I1和I2。
输入在两个时相获取的同一地区的两个多光谱图像集:I1={A1 b}和I2={A2 b},其中,At b为两个多光谱图像集合中的每一幅单波段图像,上标b表示波段序号,b=1,2,…,B,B为总波段数,下标t为时相序号,t={1,2},每一幅单波段图像At b均由m行n列像素构成。
步骤2,对两个时相多光谱图像集I1和I2分别进行维纳滤波去噪并归一化处理,得到归一化图像集和
2.1)将两时相多光谱图像集I1和I2中每一幅图像At b的灰度值区间从0~255转换到0~1之间,然后采用窗口大小为3×3像素的维纳滤波进行去噪,得到相应的去噪后图像 其中,为去噪后图像中的像素点(x,y)的灰度值;
2.2)对去噪后图像进行归一化处理,得到归一化图像 其中,为归一化图像中的像素点(x,y)的灰度值,其计算公式如下:
其中,表示取去噪后图像中所有灰度值中的最大值,表示取去噪后图像中所有灰度值中的最小值,x为行序号,y为列序号;
2.3)用时相一和时相二的归一化图像分别构成归一化图像集和即
步骤3,分别对归一化图像集和中每一幅单波段图像,采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值,得到两个时相每一波段的相对地物光谱反射率图像R1 b和R2 b,构成两个时相的相对地物光谱反射率图像集R1和R2。
3.1)对归一化图像集和中第t时相第b波段的归一化图像将该图像中的像素点(x,y)的灰度值转换为相对地物光谱反射率值rt b(x,y)的对数:
其中,x和y分别为归一化图像中像素的行序号和列序号,x=1,2,…,m,y=1,2,…,n;
3.2)对logrt b(x,y)做指数变换,得到像素点(x,y)的相对地物光谱反射率值rt b(x,y);
3.3)对第t时相第b波段的归一化图像的所有m×n个像素点均重复步骤3.1)和步骤3.2),得到第t时相第b波段的相对地物光谱反射率图像Rt b={rt b(x,y)},t={1,2};
3.4)对时相一所有B个波段的归一化图像重复步骤3.1)至步骤3.3),得到时相一的B个波段的相对地物光谱反射率图像{R1 b},构成时相一的相对地物光谱反射率图像集R1={R1 b|b=1,2,…,B};
3.5)对时相二所有B个波段的归一化图像重复步骤3.1)至步骤3.3),得到时相二的B个波段的相对地物光谱反射率图像{R2 b},构成时相二的相对地物光谱反射率图像集R2={R2 b|b=1,2,…,B}。
步骤4,分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率变化模值mr(x,y),得到光谱反射率变化模值图Mr={mr(x,y)}。
4.1)将时相一各波段的相对地物光谱反射率图像{R1 b}与时相二对应波段的相对地物光谱反射率图像{R2 b}中任意一点(x,y)的像素值作差,得到该像素点的光谱反射率变化量Dr(x,y)={drb(x,y)|b=1,2,…,B},其中drb(x,y)=r1 b(x,y)–r2 b(x,y);
4.2)计算像素点(x,y)的光谱反射率变化模值mr(x,y);
4.3)对两时相的相对地物光谱反射率图像集R1和R2中的所有m×n个像素均重复步骤4.1)和步骤4.2),得到光谱反射率变化模值图Mr={mr(x,y)│x=1,2,…,m,y=1,2,…,n}。
步骤5,分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y),得到光谱反射率角制图Sr={sr(x,y)}。
5.1)对两时相各波段的相对地物光谱反射率图像{R1 b}和{R2 b}中的任意空间位置像素点(x,y),计算该点的光谱反射率角值sr(x,y);
5.2)对两时相的相对地物光谱反射率图像集中的所有m×n个像素均重复步骤5.1),得到光谱反射率角制图Sr={sr(x,y)│x=1,2,…,m,y=1,2,…,n}。
步骤6,分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y),得到邻域差异图Qr={qr(x,y)}。
6.1)对第t时相第b波段的相对地物光谱反射率图像{Rt b}中的任意空间位置像素点(x,y),计算该点的k×k邻域内的相对地物反射率值的均值ut b(x,y);
其中c=(k-1)/2,k的取值为奇数,k=3或5,本发明实例中k=3,v和w分别表示以中心像素(x,y)为坐标原点时其对应的k×k邻域内像素点的横坐标值与纵坐标值;
6.2)计算两时相第b波段像素点(x,y)的k×k邻域内的邻域变化值gb(x,y);
6.3)计算像素点(x,y)的邻域差异值qr(x,y);
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={mri│i=1,2,…,M},其中M=m×n,mri为光谱反射率变化模值图Mr中第i个模值;
7.2)计算XM中任意两个值间的相似度wm(i,j),得到模值相似度矩阵WM={wm(i,j)│i=1,2,…,M,j=1,2,…,M},其中wm(i,j)的计算公式如下:
其中mri为光谱反射率变化模值图Mr中第i个模值,mrj为光谱反射率变化模值图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={sri│i=1,2,…,M},其中M=m×n,sri为光谱反射率角制图Sr中第i个光谱角值;
8.2)计算XS中任意两个值间的相似度ws(i,j),得到光谱角相似度矩阵WS={ws(i,j)│i=1,2,…,M,j=1,2,…,M},其中ws(i,j)的计算公式如下:
其中sri为光谱反射率角制图Sr中第i个光谱角值,srj为光谱反射率角制图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={qri│i=1,2,…,M},其中M=m×n,qri为邻域差异图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)的计算公式如下:
其中qri为邻域差异图Qr中第i个邻域差异值,qrj为邻域差异图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),相对应的概率为pc1(x,y)和pn1(x,y);
10.2)计算CM2中像素点(x,y)的k×k邻域内变化类的个数c2(x,y)和非变化类的个数nc2(x,y),相应的概率为pc2(x,y)和pn2(x,y);
10.3)计算CM3中像素点(x,y)的k×k邻域内变化类的个数c3(x,y)和非变化类的个数nc3(x,y),相应的概率为pc3(x,y)和pn3(x,y);
10.4)计算像素点(x,y)属于变化类的混合概率cw(x,y)和属于非变化类的混合概率nw(x,y);
cw(x,y)=c1(x,y)·pc1(x,y)+c2(x,y)·pc2(x,y)+c3(x,y)·pc3(x,y);
nw(x,y)=nc1(x,y)·pn1(x,y)+nc2(x,y)·pn2(x,y)+nc3(x,y)·pn3(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 (1)
1.一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法,其特征在于,包括如下步骤:
(1)输入在两个时相获取的同一地区的两个多光谱图像集:I1={A1 b}和I2={A2 b},其中,I1为时相一多光谱图像集,I2为时相二多光谱图像集,At b为两个多光谱图像集中的每一幅单波段图像,上标b表示波段序号,b=1,2,…,B,B为总波段数,下标t为时相序号,t={1,2},每一幅单波段图像At b均由m行n列像素构成;
(2)对两时相多光谱图像集I1和I2分别进行维纳滤波去噪,并归一化处理,得到两时相多光谱归一化图像集和
(3)分别对归一化图像集和中每一幅单波段图像,采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值,得到两个时相每一波段的相对地物光谱反射率图像R1 b和R2 b,并将相同时相的每一波段相对地物光谱反射率图像分别构成两个时相的相对地物光谱反射率图像集R1和R2;
(4)分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率变化模值mr(x,y),得到光谱反射率变化模值图Mr={mr(x,y)},其中x为行序号,y为列序号;
(5)分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的光谱反射率角值sr(x,y),得到光谱反射率角制图Sr={sr(x,y)};
该步骤通过如下步骤进行:
5.1)对两时相各波段的相对地物光谱反射率图像{R1 b}和{R2 b}中的任意空间位置像素点(x,y),计算该点的光谱反射率角值sr(x,y);
<mrow>
<mi>s</mi>
<mi>r</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>arccos</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>B</mi>
</munderover>
<mrow>
<mo>(</mo>
<msubsup>
<mi>r</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>&CenterDot;</mo>
<msubsup>
<mi>r</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msqrt>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>B</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>r</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
<msqrt>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>B</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>r</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
5.2)对两时相的相对地物光谱反射率图像集中的所有m×n个像素均重复步骤5.1),得到光谱反射率角制图Sr={sr(x,y)│x=1,2,…,m,y=1,2,…,n};
(6)分别计算两个时相的相对地物光谱反射率图像集R1和R2中各波段对应空间位置像素(x,y)的邻域差异值qr(x,y),得到邻域差异图Qr={qr(x,y)};
该步骤通过如下步骤进行:
6.1)对第t时相第b波段的相对地物光谱反射率图像{Rt b}中的任意空间位置像素点(x,y),计算该点的k×k邻域内的相对地物反射率值的均值ut b(x,y);
<mrow>
<msubsup>
<mi>u</mi>
<mi>t</mi>
<mi>b</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mi>k</mi>
<mo>&CenterDot;</mo>
<mi>k</mi>
</mrow>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>v</mi>
<mo>=</mo>
<mo>-</mo>
<mi>c</mi>
</mrow>
<mi>c</mi>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>w</mi>
<mo>=</mo>
<mo>-</mo>
<mi>c</mi>
</mrow>
<mi>c</mi>
</munderover>
<msubsup>
<mi>r</mi>
<mi>t</mi>
<mi>b</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>+</mo>
<mi>v</mi>
<mo>,</mo>
<mi>y</mi>
<mo>+</mo>
<mi>w</mi>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
其中c=(k-1)/2,k的取值为奇数,k=3或5,v和w分别表示以中心像素(x,y)为坐标原点时其对应的k×k邻域内像素点的横坐标值与纵坐标值;
6.2)计算两时相第b波段像素点(x,y)的k×k邻域内的邻域变化值gb(x,y);
<mrow>
<msup>
<mi>g</mi>
<mi>b</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mi>k</mi>
<mo>&CenterDot;</mo>
<mi>k</mi>
</mrow>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>v</mi>
<mo>=</mo>
<mo>-</mo>
<mi>c</mi>
</mrow>
<mi>c</mi>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>w</mi>
<mo>=</mo>
<mo>-</mo>
<mi>c</mi>
</mrow>
<mi>c</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>r</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>+</mo>
<mi>v</mi>
<mo>,</mo>
<mi>y</mi>
<mo>+</mo>
<mi>w</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>r</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>+</mo>
<mi>v</mi>
<mo>,</mo>
<mi>y</mi>
<mo>+</mo>
<mi>w</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>;</mo>
</mrow>
6.3)计算像素点(x,y)的邻域差异值qr(x,y);
<mrow>
<mi>q</mi>
<mi>r</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>B</mi>
</munderover>
<mi>m</mi>
<mi>i</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>r</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>u</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>,</mo>
<msubsup>
<mi>r</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>u</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>B</mi>
</munderover>
<mi>max</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>r</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>u</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>,</mo>
<msubsup>
<mi>r</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<msubsup>
<mi>u</mi>
<mn>2</mn>
<mi>b</mi>
</msubsup>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<mrow>
<mo>(</mo>
<mfrac>
<mn>1</mn>
<mi>B</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>b</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>B</mi>
</munderover>
<mfrac>
<mrow>
<msup>
<mi>g</mi>
<mi>b</mi>
</msup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msubsup>
<mi>u</mi>
<mn>1</mn>
<mi>b</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>;</mo>
</mrow>
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={mri│i=1,2,…,M},其中M=m×n,mri为光谱反射率变化模值图Mr中第i个模值;
7.2)计算XM中任意两个值间的相似度wm(i,j),得到模值相似度矩阵WM={wm(i,j)│i=1,2,…,M,j=1,2,…,M},其中wm(i,j)的计算公式如下:
<mrow>
<mi>w</mi>
<mi>m</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>mr</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>mr</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msup>
<mo>;</mo>
</mrow>
其中mri为光谱反射率变化模值图Mr中第i个模值,mrj为光谱反射率变化模值图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={sri│i=1,2,…,M},其中M=m×n,sri为光谱反射率角制图Sr中第i个光谱角值;
8.2)计算XS中任意两个值间的相似度ws(i,j),得到光谱角相似度矩阵WS={ws(i,j)│i=1,2,…,M,j=1,2,…,M},其中ws(i,j)的计算公式如下:
<mrow>
<mi>w</mi>
<mi>s</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>sr</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>sr</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msup>
<mo>;</mo>
</mrow>
其中sri为光谱反射率角制图Sr中第i个光谱角值,srj为光谱反射率角制图Sr中第j个光谱角值;
8.3)将光谱角相似度矩阵WS中的每一列元素相加,得到的结果放在光谱角相似度矩阵WS相应的对角线上,并令光谱角相似度矩阵WS中其它位置的元素为0,得到光谱角对角矩阵DS,再令LS=DS-WS;
8.4)求LS的前K个递增的特征值对应的特征向量,把K个特征向量放在一起构造一个M×K的矩阵MS;
8.5)将MS的每个行向量看成是K维空间中的一点,采用K-means算法进行聚类;
8.6)若MS某一行被归为变化类,则光谱反射率角制图Sr相应位置处的像素点也归为变化类,若MS某一行被归为非变化类,则光谱反射率角制图Sr相应位置处的像素点也归为非变化类,从而得到初始变化检测结果图CM2;
(9)用谱聚类方法将邻域差异图Qr中的像素聚为两类,分别为变化类和非变化类,得到初始变化检测结果图CM3;
该步骤通过如下步骤进行:
9.1)将邻域差异图Qr按列变成M×1的矩阵XQ={qri│i=1,2,…,M},其中M=m×n,qri为邻域差异图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)的计算公式如下:
<mrow>
<mi>w</mi>
<mi>q</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>qr</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>qr</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msup>
<mo>;</mo>
</mrow>
其中qri为邻域差异图Qr中第i个邻域差异值,qrj为邻域差异图Qr中第j个邻域差异值;
9.3)将邻域差相似度矩阵WQ中的每一列元素相加,得到的结果放在邻域差相似度矩阵WQ相应的对角线上,并令邻域差相似度矩阵WQ中其它位置的元素为0,得到邻域差对角矩阵DQ,再令LQ=DQ-WQ;
9.4)求LQ的前K个递增的特征值对应的特征向量,把K个特征向量放在一起构造一个M×K的矩阵MQ;
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),相对应的概率为pc1(x,y)和pn1(x,y);
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>c</mi>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mi>c</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>n</mi>
<mi>c</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>n</mi>
<mn>1</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mi>n</mi>
<mi>c</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>n</mi>
<mi>c</mi>
<mn>1</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
10.2)计算CM2中像素点(x,y)的k×k邻域内变化类的个数c2(x,y)和非变化类的个数nc2(x,y),相应的概率为pc2(x,y)和pn2(x,y);
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>c</mi>
<mn>2</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mi>c</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>n</mi>
<mi>c</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>n</mi>
<mn>2</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mi>n</mi>
<mi>c</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>n</mi>
<mi>c</mi>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
10.3)计算CM3中像素点(x,y)的k×k邻域内变化类的个数c3(x,y)和非变化类的个数nc3(x,y),相应的概率为pc3(x,y)和pn3(x,y);
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>c</mi>
<mn>3</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mi>c</mi>
<mn>3</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mn>3</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>n</mi>
<mi>c</mi>
<mn>3</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
<mrow>
<msub>
<mi>p</mi>
<mrow>
<mi>n</mi>
<mn>3</mn>
</mrow>
</msub>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mi>n</mi>
<mi>c</mi>
<mn>3</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mn>3</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>n</mi>
<mi>c</mi>
<mn>3</mn>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>;</mo>
</mrow>
10.4)计算像素点(x,y)属于变化类的混合概率cw(x,y)和属于非变化类的混合概率nw(x,y);
cw(x,y)=c1(x,y)·pc1(x,y)+c2(x,y)·pc2(x,y)+c3(x,y)·pc3(x,y);
nw(x,y)=nc1(x,y)·pn1(x,y)+nc2(x,y)·pn2(x,y)+nc3(x,y)·pn3(x,y);
10.5)对像素点(x,y),如果cw(x,y)大于等于nw(x,y),则该像素点为变化类,否则,该像素点为非变化类;
10.6)对CM1、CM2和CM3中的所有m×n个像素均重复步骤10.1)至步骤10.5),得到最终变化检测结果图CM。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410770398.7A CN104463881B (zh) | 2014-12-12 | 2014-12-12 | 一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410770398.7A CN104463881B (zh) | 2014-12-12 | 2014-12-12 | 一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104463881A CN104463881A (zh) | 2015-03-25 |
CN104463881B true CN104463881B (zh) | 2018-03-16 |
Family
ID=52909863
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410770398.7A Expired - Fee Related CN104463881B (zh) | 2014-12-12 | 2014-12-12 | 一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104463881B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106841055B (zh) * | 2017-03-22 | 2019-06-14 | 浙江大学 | 一种重构艺术绘画光谱图像的训练样本选择方法 |
CN107992891B (zh) * | 2017-12-01 | 2022-01-25 | 西安电子科技大学 | 基于光谱矢量分析多光谱遥感图像变化检测方法 |
CN108021887B (zh) * | 2017-12-05 | 2019-10-01 | 中国科学院遥感与数字地球研究所 | 基于空间光谱差比参量的遥感影像分析方法及应用 |
CN110706188B (zh) * | 2019-09-23 | 2020-10-23 | 北京航天宏图信息技术股份有限公司 | 一种影像融合方法、装置、电子设备及存储介质 |
CN110779876B (zh) * | 2019-11-07 | 2022-12-02 | 长光禹辰信息技术与装备(青岛)有限公司 | 一种疫木识别方法、装置、设备及计算机可读存储介质 |
CN113627357B (zh) * | 2021-08-13 | 2022-07-08 | 哈尔滨工业大学 | 一种高空间-高光谱分辨率遥感图像本征分解方法及*** |
CN113902717B (zh) * | 2021-10-13 | 2022-04-26 | 自然资源部国土卫星遥感应用中心 | 一种基于光谱库的星载高光谱农田裸土目标识别方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102663752A (zh) * | 2012-04-11 | 2012-09-12 | 南京理工大学 | 一种sam加权kest高光谱异常检测算法 |
CN103198483A (zh) * | 2013-04-07 | 2013-07-10 | 西安电子科技大学 | 基于边缘和光谱反射率曲线的多时相遥感图像配准方法 |
CN103226832A (zh) * | 2013-05-07 | 2013-07-31 | 西安电子科技大学 | 基于光谱反射率变化分析的多光谱遥感影像变化检测方法 |
-
2014
- 2014-12-12 CN CN201410770398.7A patent/CN104463881B/zh not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102663752A (zh) * | 2012-04-11 | 2012-09-12 | 南京理工大学 | 一种sam加权kest高光谱异常检测算法 |
CN103198483A (zh) * | 2013-04-07 | 2013-07-10 | 西安电子科技大学 | 基于边缘和光谱反射率曲线的多时相遥感图像配准方法 |
CN103226832A (zh) * | 2013-05-07 | 2013-07-31 | 西安电子科技大学 | 基于光谱反射率变化分析的多光谱遥感影像变化检测方法 |
Non-Patent Citations (5)
Title |
---|
Fusion of Difference Images for Change Detection in Urban Areas;Sicong Liu 等;《ResearchGate》;20110430;第1-4页 * |
基于RS的现代黄河三角洲土地利用/覆被变化检测方法研究;李翔宇;《万方学位论文数据库》;20090115;正文第39页 * |
基于treelet的遥感图像变化检测方法研究;曹娟;《万方学位论文数据库》;20120731;正文第51-52页 * |
基于谱聚类的SAR图像变化检测;李泽敏;《万方学位论文数据库》;20120731;正文第18-21页 * |
多波段信息融合的遥感影像变化检测;魏立飞 等;《武汉大学学报 信息科学版》;20140131;第39卷(第1期);第8-11页 * |
Also Published As
Publication number | Publication date |
---|---|
CN104463881A (zh) | 2015-03-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104463881B (zh) | 一种基于光谱反射率邻域差异图和邻域概率融合的多光谱遥感影像变化检测方法 | |
CN108389188B (zh) | 一种稀疏高光谱异常目标检测方法 | |
CN103226832B (zh) | 基于光谱反射率变化分析的多光谱遥感影像变化检测方法 | |
CN103440505B (zh) | 空间邻域信息加权的高光谱遥感图像分类方法 | |
CN109902715B (zh) | 一种基于上下文聚合网络的红外弱小目标检测方法 | |
CN102819740B (zh) | 一种单帧红外图像弱小目标检测和定位方法 | |
CN102254319A (zh) | 一种多层次分割的遥感影像变化检测方法 | |
CN106339674A (zh) | 基于边缘保持与图割模型的高光谱影像分类方法 | |
CN103226826B (zh) | 基于局部熵视觉注意模型的遥感图像变化检测方法 | |
CN109712149B (zh) | 一种基于小波能量和模糊c均值的图像分割方法 | |
CN103093478B (zh) | 基于快速核空间模糊聚类的异源图像粗大边缘检测方法 | |
Tang et al. | A multiple-point spatially weighted k-NN method for object-based classification | |
CN111046800A (zh) | 一种基于低秩与稀疏分解的高光谱图像异常目标检测方法 | |
CN109544494B (zh) | 一种人体安检中被动毫米波图像与可见光图像的融合方法 | |
CN115205590A (zh) | 一种基于互补集成Transformer网络的高光谱图像分类方法 | |
Liang et al. | Maximum likelihood classification of soil remote sensing image based on deep learning | |
Zhang et al. | Multi-temporal cloud detection based on robust PCA for optical remote sensing imagery | |
CN116664892A (zh) | 基于交叉注意与可形变卷积的多时相遥感图像配准方法 | |
CN103218823B (zh) | 基于核传播的遥感图像变化检测方法 | |
Chen et al. | Scene segmentation of remotely sensed images with data augmentation using U-net++ | |
CN103093472B (zh) | 基于双字典交叉稀疏表示的光学遥感图像变化检测方法 | |
CN110276782B (zh) | 一种结合空谱特征和相关滤波的高光谱目标跟踪方法 | |
Su et al. | Enhancing concealed object detection in Active Millimeter Wave Images using wavelet transform | |
CN112784777B (zh) | 基于对抗学习的无监督高光谱图像变化检测方法 | |
CN113421198A (zh) | 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180316 Termination date: 20181212 |