CN109360184A - 结合阴影补偿与决策融合的遥感影像变化检测方法 - Google Patents
结合阴影补偿与决策融合的遥感影像变化检测方法 Download PDFInfo
- Publication number
- CN109360184A CN109360184A CN201810968569.5A CN201810968569A CN109360184A CN 109360184 A CN109360184 A CN 109360184A CN 201810968569 A CN201810968569 A CN 201810968569A CN 109360184 A CN109360184 A CN 109360184A
- Authority
- CN
- China
- Prior art keywords
- scale
- image
- shadow
- remote sensing
- sensing image
- 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.)
- Pending
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 76
- 230000004927 fusion Effects 0.000 title claims abstract description 27
- 238000000034 method Methods 0.000 claims description 63
- 230000001154 acute effect Effects 0.000 claims description 11
- 230000011218 segmentation Effects 0.000 claims description 10
- 230000000007 visual effect Effects 0.000 claims description 6
- 238000000546 chi-square test Methods 0.000 claims description 5
- 239000000571 coke Substances 0.000 claims description 4
- 238000005259 measurement Methods 0.000 claims description 4
- 238000013139 quantization Methods 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 239000013598 vector Substances 0.000 claims description 3
- 230000015572 biosynthetic process Effects 0.000 claims description 2
- 238000003786 synthesis reaction Methods 0.000 claims description 2
- VMXUWOKSQNHOCA-UKTHLTGXSA-N ranitidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CC=C(CN(C)C)O1 VMXUWOKSQNHOCA-UKTHLTGXSA-N 0.000 abstract description 12
- 238000012360 testing method Methods 0.000 description 19
- 238000002474 experimental method Methods 0.000 description 17
- 238000011156 evaluation Methods 0.000 description 13
- 238000007689 inspection Methods 0.000 description 7
- 238000000605 extraction Methods 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 238000012937 correction Methods 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 238000005315 distribution function Methods 0.000 description 2
- 238000010191 image analysis Methods 0.000 description 2
- 241001270131 Agaricus moelleri Species 0.000 description 1
- 241000208340 Araliaceae Species 0.000 description 1
- 241000406668 Loxodonta cyclotis Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 239000013604 expression vector Substances 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000002844 melting Methods 0.000 description 1
- 230000008018 melting Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 238000000926 separation method Methods 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
-
- G06T5/94—
-
- 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
Abstract
本发明公开了结合阴影补偿与决策融合的遥感影像变化检测方法。在面向对象的变化检测框架下，首先提取遥感影像中的地物阴影，然后对多尺度变化检测进行阴影补偿。其中，通过构建一种尺度间互信息最小化的目标函数实现了尺度参数的自适应提取；在此基础上，结合所提出的阴影补偿因子，设计了一种基于D‑S证据理论的决策级多尺度融合策略，并进一步对变化强度的等级进行了划分。本发明能够较好地解决阴影所导致的错检问题，显著提高变化检测精度。
Description
技术领域
本发明属于图像处理技术领域，特别涉及了一种遥感影像变化检测方法。
背景技术
遥感影像变化检测是从不同时期的遥感数据中定量分析和确定地表变化的过程。近年来，随着多时相高分辨率遥感数据的不断积累，如何从中提取和检测城市场景中的变化信息已成为遥感科学和地理信息科学的重要研究课题。与中、低分辨率遥感影像相比，高分辨率遥感影像的光谱、纹理及空间细节信息等都更加丰富。与此同时，空间分辨率的提高也使变化检测面临着更加突出的“同物异谱”和“同谱异物”问题。为此，学者们试图利用面向对象的图像分析OBIA(Object-Based Image Analysis)来提高变化检测的精度。
与传统像素级检测方法相比，对象级变化检测OBCD(Object-Based ChangeDetection)基于对象固有形状及尺寸进行特征提取，对配准误差、噪声等具有更高的鲁棒性。例如，Chen等人提出的OB-HMAD方法以地理对象作为基本单元，依据最小噪声分离变换MNF(Minimum Noise Fraction Rotation)确定分割阈值来提取变化信息，其检测精度显著优于传统像素级方法；佃袁勇等通过多尺度的图像分割、特征提取及融合来描述变化信息，相较于基于单一尺度的变化检测方法错检率显著降低。尽管如此，现有多尺度OBCD方法中大多直接比较地理对象间光谱、纹理特征的差异，而忽略了地物阴影的影响。特别是在城市场景中，各种人造、自然地物阴影普遍存在且密集分布，地物阴影所导致的“伪变化”已成为产生错检的主要来源之一。因此，有必要在多尺度变化检测中引入阴影补偿策略，而其关键在于准确的阴影提取、尺度参数的合理选择，以及光谱、纹理、阴影及尺度信息的有效融合。
发明内容
为了解决上述背景技术提出的技术问题，本发明旨在提供结合阴影补偿与决策融合的遥感影像变化检测方法，解决阴影所导致的错检问题，显著提高变化检测精度。
为了实现上述技术目的，本发明的技术方案为：
结合阴影补偿与决策融合的遥感影像变化检测方法，括以下步骤：
(1)基于矢量量化方法对遥感影像进行分类，基于分类结果，将阴影像素视为出界点，采用卡方检验进行阴影检测；
(2)采用多尺度分割算法提取地理对象；
(3)基于尺度间互信息最小的优化目标自适应提取尺度参数，获取多尺度J-image影像序列；
(4)获取步骤(1)中检测出的阴影和步骤(2)中提取的地理对象在不同时相、相同尺度J-image影像序列中的投影，并在此基础上，分析任意地理对象在不同时相、相同尺度J-image影像序列间的相似程度；
(5)通过结合阴影补偿的决策融合算法实现遥感影像变化检测及变化强度等级的划分。
进一步地，在步骤(1)中，采用卡方检验进行阴影检测的方法如下：
定义卡方分布：
(X-m)TΣ-1(X-m)～χ2(b)
其中，X为随机变量，m和Σ分别为非阴影区域的均值和协方差矩阵，b为自由度，即多光谱影像的波段数目；
在此基础上，采用下式确定阴影像素：
其中，P表示概率，1-γ为置信度，此时卡方值小于的像素为阴影像素。
进一步地，在步骤(2)中，多尺度分割仅在单一时相影像中进行，再将分割结果根据配准获得的像素间匹配关系直接映射到另一时相影像中，从而确定统一的地理对象集合R＝{R1,R2,...RG}，其中G为分割结果中的对象总数。
进一步地，在步骤(3)中，自适应提取尺度参数的过程如下：
(31)依据目视解译在影像中选择尺寸较大的地理对象，统计其包含的像素数量为E；
(32)设定尺度参数的最大值HMAX为大于的最小整数；
(33)依据阴影检测及配准结果，选择在两时相影像中均为非阴影的像素集合，利用该像素集合分别在单一时相影像中计算相邻尺度间的互信息为NMIt(s-1,s)；其中t∈{1,2}，代表不同时相影像；s为尺度序号，其最大值取HMAX；
(34)获得双时相影像相邻尺度间的互信息值：
(35)通过迭代计算遍历所有的尺度参数组合，从而自适应提取NMIall最小值对应的尺度参数集合HNMImin＝{H1,H2,...HD}，D为HNMImin中的尺度总数。
进一步地，在步骤(4)中，设任意地理对象Ri在时相1、2相同尺度J-image中的投影分别为Ri1和Ri2，则Ri1与Ri2间的相似程度：
其中，μi1、μi2、σi1、σi2、σi1i2分别是Ri1与Ri2的均值、标准差、方差和协方差；C1、C2为常数；SSIM越大，则当前尺度下Ri1与Ri2的相似程度越高，发生变化的可能性越小；遍历所有D个尺度，可获得Ri对应多尺度J-image序列下的相似性度量集合SSIMki∈{SSIM1i,SSIM2i,...,SSIMDi}，其中k为尺度参数集合HNMImin中的尺度序号。
进一步地，在步骤(5)中，首先定义识别框架U＝{SL,MA,UN}，将地理对象划分为剧烈变化类SL、显著变化类MA和未变化类UN，则焦元A包括{SL}，{MA}，{UN}，{SL,MA,UN}；对任意地理对象Ri，建立如下基本概率赋值函数：
mki({SL})＝(1-SSIMki)×0.7×αk×λ
mki({MA})＝(1-SSIMki)×0.3×αk×λ
mki({UN})＝SSIMki×αk×λ
mki({SL,MA,UN})＝1-αk×λ
其中，αk用于描述对k尺度所提供证据的信任程度，k∈(1,2,...D)，λ为阴影补偿因子；
然后采用下式对建立的基本概率赋值函数进行多尺度证据合成：
其中，mi(A)为合成后的基本概率赋值函数，Ak为k尺度下的焦元，Vi为归一化常数：
最后采用以下判别规则对变化强度等级进行划分：
若满足以下条件之一，则地理对象Ri属于剧烈变化类SL：
条件1、mi({SL})＞0.8；
条件2、mi({SL})＞0.6且mi({MA})＞0.3；
若满足以下条件之一，则地理对象Ri属于显著变化类MA：
条件1、mi({MA})＞0.7且mi({SL})＞0.1；
条件2、mi({UN)}＜0.1
若按照上述判别条件地理对象Ri不属于剧烈变化类SL和显著变化类MA，则地理对象Ri属于未变化类UN。
进一步地，阴影补偿因子λ用于描述阴影对于变化检测的干扰程度，定义如下：
λ＝l1β1+l2β2+l3β3
其中，β1为Ri在双时相影像中均为非阴影的像素比例，β2为Ri在双时相影像中均为阴影，β3为Ri仅在一个时相影像中为阴影的像素比例，l1、l2、l3分别为β1、β2、β3的权重，l1＞l2＞l3且l1+l2+l3＝1。
进一步地，设定l1＝0.6，l1＝0.3，l3＝0.1。
进一步地，将D个尺度的αk统一设定为0.5。
采用上述技术方案带来的有益效果：
本发明不仅构建了一种基于尺度间互信息最小化的目标函数以自适应提取尺度参数集合，还提出了一种结合阴影补偿因子的多尺度决策级融合策略，进而实现了对变化强度等级的划分。实验证明，本发明能够有效减少阴影所导致的错检问题，显著提高变化检测精度且具有良好的鲁棒性。
附图说明
图1是本发明的方法流程图；
图2是实验数据集1的遥感影像图；
图3是实验数据集2的遥感影像图；
图4是实验数据集1的阴影检测结果图；
图5是实验数据集2的阴影检测结果图；
图6是实验数据集1、2的本方法变化检测结果图；
图7是实验数据集1、2的OB-HMAD变化检测结果图；
图8是实验数据集1、2的CVA-EM变化检测结果图；
图9是实验数据集3的遥感影像图；
图10是实验数据集3的OB-HMAD变化检测结果图；
图11是实验数据集3的本方法变化检测结果图；
图12是实验数据集1、2、3的变化参考图。
具体实施方式
以下将结合附图，对本发明的技术方案进行详细说明。
首先假设多时相影像已经过几何配准与辐射校正预处理，则本发明所提出方法主要包括四个步骤：首先，分别对多时相影像进行阴影检测；其次，选择单一时相影像进行分割，提取统一的地理对象集合；在此基础上，结合阴影检测结果搜寻目标函数最小值，自适应提取尺度参数集合；最后，综合多种特征进行决策融合，获得最终变化检测结果。方法流程如图1所示。
1、阴影检测
高斯分布背景模型认为影像的自然背景符合高斯分布，则人造地物作为出界点可被检测出来。而对于人造地物密集分布的城市场景，如果采用VQ分类方法将人造地物与自然背景分别划分为一类或者几类，分类结果则依然符合高斯分布，阴影则可以作为出界点被检测出来。
基于以上假设，首先采用Bai S等提出的VQ方法对影像进行分类：定义同组滤波器PGF(Peer Group Filter)，将滤波所获得的局部统计特性作为权重，进行VQ初始量化；在此基础上，采用GLA(Generalized Lloyd Algorithm)对矢量量化结果进行分类。基于分类结果，将阴影像素视为出界点，采用卡方检验进行阴影检测。定义卡方分布为：
(X-m)TΣ-1(X-m)～χ2(b) (1)
其中X为随机变量，m和Σ分别为非阴影区域的均值和协方差矩阵，b为自由度(即多光谱影像的波段数目)。在此基础上，采用公式(2)确定阴影像素：
其中，1-γ为置信度，实验中置信度通过试错法取最优值确定。此时，卡方值小于的像素为阴影像素。
2、地理对象提取
通过分割来提取地理对象是进行对象级变化检测的基础，采用发明人之前提出的多尺度分割方法WJSEG，原因在于：相较于知名商业软件eCognition，WJSEG不仅能够准确定位对象边缘，且保持对象轮廓更加完整，更有利提高变化检测过程的方法透明度。
需要说明的是，分割仅在单一时相影像中进行，再将分割结果根据配准获得的像素间匹配关系直接映射到另一时相影像中，从而确定统一的地理对象集合R＝{R1,R2,...RG}，其中G为分割结果中的对象总数。
3、J-image及尺度参数自适应提取
与传统小波、轮廓波变换等相比，J-image序列不仅能够描述不同尺寸局部区域的光谱、纹理复杂程度，还具有对高频信息的方向性不敏感的特点。因此，本发明采用J-image影像序列作为变化检测的多尺度分析平台。
3.1、多尺度J-image序列
单一尺度的J-image由局部区域同质性指标J-value构成的，J-value定义为：令每一个像素的位置zf(x,y)为像素zf的像素值，f∈{1,2...F}，F为影像中的像素总数。确定尺寸为H×H像素的窗口(H即为尺度参数)，Zf为以像素zf为中心的窗口内所有像素的集合，zf(x,y)∈Zf。则J-value可由公式(3)计算得到，其中SfT为Zf中所有像素的总体方差；SfW为对Zf中所有像素以256个灰度级(0～255)为中心分别计算方差，再对方差求和的结果。
J-value＝(SfT-SfW)/SfW (3)
将J-value替换zf的像素值并遍历所有像素，可获得单一尺度的J-image；进而通过改变尺度参数H，可获得多尺度的J-image影像序列。
3.2、尺度参数自适应提取
尺度参数集合的选择是否合理是影响多尺度变化检测效果的关键问题之一，体现在：尺度参数应当与场景中代表性地物的尺寸相近，从而有利于准确描述此类地物的变化信息；若尺度过多，则必然存在大量的冗余信息；若尺度过少，又容易陷入局部最优。由于互信息能够反映尺度间的关联程度，因此互信息越小时尺度间包含的冗余信息越少；同时，阴影在变化检测中作为一种干扰因素，不应当参与尺度间互信息的计算。由此，本发明提出了一种基于尺度间互信息最小化的目标函数优化策略，以实现尺度参数的自适应选择。
Step1：依据目视解译在影像中选择尺寸较大的地理对象(如湖泊、操场等)，统计其包含的像素数量为E；
Step2：设定尺度参数的最大值HMAX为大于的最小整数；
Step3：依据阴影检测及配准结果，选择在两时相影像中均为非阴影的像素集合。利用该像素集合分别在单一时相影像中计算相邻尺度间的互信息为NMIt(s-1,s)。其中t∈{1,2}，代表不同时相影像；s为尺度序号，最大取HMAX。
Step4：根据公式(4)构建目标函数，从而获得双时相影像相邻尺度间的互信息值。
NMIall＝-[NMI1(1，2)+NMI2(1，2)+NMI1(2，3)+NMI2(2，3)+...] (4)
Step5：通过迭代计算遍历所有可能的尺度参数组合，从而自适应提取NMIall最小值对应的尺度参数集合HNMImin＝{H1,H2,...HD}，并记HNMImin中尺度总数为D。
4、多尺度决策融合及变化检测
依据HNMImin及公式(3)，首先计算每个时相原始影像的多尺度J-image序列；同时，由于J-image的尺寸与原始影像相同，依据WJSEG分割及阴影检测结果可直接获得地理对象集合R及阴影在不同时相、相同尺度J-image中的投影；在此基础上，分析任意对象Ri∈R在不同时相、相同尺度J-image间的相似程度；最后，通过多尺度决策级融合实现变化检测及变化强度等级的划分。
4.1、多尺度相似性度量
本发明采用SSIM来度量Ri在不同时相、相同尺度J-image之间的相似程度，原因在于：与传统的直方图匹配及欧氏距离相比，结构自相似性SSIM能够同时从均值、方差和协方差三个方面综合反映向量间的相似性。计算过程如下：
设Ri在时相1、2相同尺度J-image中的投影分别为Ri1和Ri2，则Ri1与Ri2间的SSIM为：
其中，μi1，μi2，σi1，σi2，σi1i2分别是Ri1与Ri2的均值、标准差、方差和协方差。C1，C2为常数，通常取C1＝0.2，C2＝0.8。SSIM越大，则当前尺度下Ri1与Ri2的相似程度越高，发生变化的可能性越小。遍历所有D个尺度，可获得Ri对应多尺度J-image序列下的相似性度量集合SSIMki∈{SSIM1i,SSIM2i,...,SSIMDi}，其中k为尺度参数集合HNMImin中的尺度序号。
4.2、结合阴影补偿的D-S决策融合
基于阴影检测结果与相似性度量集合，本发明提出了一种综合光谱、纹理特征与阴影补偿因子的D-S决策融合策略，从而实现变化检测及变化强度等级的划分。D-S证据理论具有无需先验概率，能够直接利用异源信息作为证据支持的优点，Dempster合成规则如下：
设A为焦元，U是识别框架，由于尺度总数为D，则对于U上的D个mass函数m1,m2...mD的Dempster合成法则为：
其中，k为尺度序号，V为归一化常数：
所提出的融合策略首先定义识别框架U＝{SL,MA,UN}，将对象划分为剧烈变化类SL、显著变化类MA和未变化类UN，则焦元包括{SL}，{MA}，{UN}，{SL,MA,UN}。对任意对象Ri，建立基本概率赋值函数(BPAF，Basic Probability Assignment Function)如下：
mki({SL})＝(1-SSIMki)×0.7×αk×λ (8)
mki({MA})＝(1-SSIMki)×0.3×αk×λ (9)
mki({UN})＝SSIMki×αk×λ (10)
mki({SL,MA,UN})＝1-αk×λ (11)
其中，αk用于描述对某一尺度所提供证据的信任程度，k∈(1,2,...D)。尽管尺度参数H越小越能够准确反映对象的细节特征，但也对噪声及孤立点更加敏感。因此，统一设定参数αk为0.5。λ为阴影补偿因子，用于描述阴影对于变化检测的干扰程度，定义为：
λ＝l1β1+l2β2+l3β3 (12)
其中，β1为Ri在双时相影像中均为非阴影的像素比例，此类像素的变化检测未受到阴影的影响，因此在补偿因子中的权重l1应最高；β2为Ri在双时相影像中均为阴影的像素比例，由于双时相阴影之间可能存在差异从而导致“伪变化”，因此权重l2应较低；β3为Ri仅在一个时相影像中为阴影的像素比例，由于此类像素受到阴影所导致的“伪变化”影响最为显著，因此权重l3应最低。由此，在l1＞l2＞l3且l1+l2+l3＝1的条件约束下，本发明根据2组实验中阴影补偿对最终变化检测结果的影响，采用试错法确定当l1＝0.6，l1＝0.3，l3＝0.1时能够取得较为理想的效果。
根据公式(12)可知，λ越大则阴影对变化检测的干扰越小。在此基础上，对任意对象Ri，采用公式(6)对所构建的BPAF进行多尺度证据合成。最后，采用如下判别规则对变化强度等级进行划分：
1、属于SL类的Ri应满足：mi({SL})＞0.8，或者mi({SL})＞0.6且mi({MA})＞0.3；
2、属于MA类的Ri应满足：mi({MA})＞0.7且mi({SL})＞0.1，或者mi({UN)}＜0.1；
3、否则Ri属于UN类。
以上判别规则的依据在于，对象Ri属于某类的可能性越大，则该类对应的概率分配函数取值应当更高。因此，以Ri属于剧烈变化对象为例，仅依据剧烈变化类概率分配函数所提供的证据，需满足mi({SL})大于阈值0.8；若依据剧烈及显著变化类概率分配函数所提供的证据，则需满足mi({MA})大于阈值0.3且mi({SL})大于阈值0.6，以此类推。需要指出的是，划分规则中的阈值均采用试错法确定，而在实际应用中可根据需求调整剧烈变化SL类的划分阈值，从而为野外勘测等提供更有价值的靶标信息。
5、实验
实验采用多组不同空间分辨率、不同传感器类型的多时相高分辨率遥感影像，并与未经阴影补偿的对象级变化检测方法OB-HMAD方法及传统像素级变化检测方法CVA-EM[16]进行比较分析。实验平台采用Matlab R2014a，处理器为Inter Core i5 3.2GHz，内存为8GB。
5.1、实验数据
实验采用的多时相影像均已经过几何配准与辐射校正。其中，数据集1为中国重庆地区的多光谱QuickBird影像，采集时间分别为2007年9月和2011年8月，空间分辨率为2.4米，图像尺寸为1024×1024像素，如图2所示。图2中的(a)为遥感影像#1，(b)为遥感影像#2。
数据集2为中国江苏南京地区的航空遥感DOM(Digital Ortho-photo Map)影像，采集时间分别为2009年3月和2012年2月，空间分辨率为0.6米，图像尺寸为1024×1024像素，如图3所示。图3中的(a)为遥感影像#3，(b)为遥感影像#4。
如上图所示，实验影像均为典型的城市场景，主要由植被、道路、阴影、建筑物及其他人造目标构成。由于两个数据集中影像的采集时间均为夏末秋初或冬末春初，植被的物候差异对变化检测的影响较小。场景中的典型变化类型主要为由植被变为建筑物，由荒地变为建筑物或植被，以及建筑物重建等。阴影主要存在于高层建筑物的背向阳光的一侧。特别是在数据集2中，由于光照条件存在明显的差异，阴影导致的“伪变化”也更加突出。
5.2、阴影检测与对象提取
阴影检测中，采用试错法确定数据集1中γ＝0.03，数据集2中γ＝0.06。将提取的阴影像素的灰度值设定为255，其他像素保持不变，检测结果如图4、图5所示。图4中的(a)为遥感影像#1的阴影检测结果，(b)为遥感影像#2的阴影检测结果。图5中的(a)为遥感影像#3的阴影检测结果，(b)为遥感影像#4的阴影检测结果。
通过与原始图像对比可以看出，影像#1、#2及#3较为真实的反映了建筑物等地物阴影。在影像#4中，由于阴影所占的面积较大，因此在一定程度上破坏了高斯分布背景模型，检测结果中存在一些漏检及错检现象。尽管如此，与其他三幅影像相比，影像#4中阴影所在区域仅灰度值有所降低，但基本保持了地表原始的纹理特征。同时，由于变化检测所依赖的J-value主要反映了局部区域纹理特征的复杂程度，因此影像#4中的阴影检测误差对变化检测的影响有限。
采用WJSEG分别对影像#2、#3进行分割，WJSEG能够准确地定位对象边缘，同时有效区分不同种类的相邻地物；对于内部光谱特征较为均匀的大尺寸对象，WJSEG保持了此类对象的完整轮廓，仅在个别局部纹理特征复杂的区域存在欠分割或过分割现象。
5.3、多尺度变化检测与融合
实验中设定HMAX＝20，通过对公式(4)求最小值，获得数据集1对应的尺度参数集合为HNMImin＝{5,7,11,14,16}，数据集2为HNMImin＝{5,9,11,13,17}。在此基础上，利用公式(5)计算每个对象在多尺度J-image中的SSIM，最后结合阴影检测结果进行决策融合。通过将不同变化等级的对象采用不同的灰度进行表示，所获得的变化检测结果如图6所示。图6中的(a)为数据集1的变化检测结果，(b)为数据集2的变化检测结果，图中的深灰色表示未变化类UN，浅灰色表示剧烈变化类SL，白色表示显著变化类MA。
基于本发明的对象提取结果，对原始影像采用OB-HMAD与CVA-EM方法进行变化检测的结果分别如图7、图8所示(白色像素代表变化，黑色像素代表未变化)。图6中的(a)为数据集1的OB-HMAD变化检测结果，(b)为数据集2的OB-HMAD变化检测结果。图8中的(a)为数据集1的CVA-EM变化检测结果，(b)为数据集2的CVA-EM变化检测结果。
通过目视分析可以看出，本方法明显优于其他两种方法，体现在：在2组实验中，对于发生变化且没有受到阴影影响的对象，仅有本发明提出的方法做出了准确的判别。OB-HMAD对于由植被变为建筑物的位置，以及由植被变为篮球场的位置均出现了漏检。CVA-EM方法由于以像素作为变化检测的基元，对于新建人造地物均出现了部分错检和漏检。对于受到阴影干扰而未发生变化的对象，三种方法均未发生错检。而对于受到阴影干扰且真实发生变化的对象，OB-HMAD和CVA-EM方法均发生了漏检，对于仅光谱特征存在较大差异但实际并未发生变化地物，仅CVA-EM发生了错检。
5.4、定量精度评价
5.4.1像素级精度评价
首先基于像素定量评价不同对象、像素级变化检测方法的性能：为保证所选择参考像素的准确性及合理分布，通过实地考察与目视分析，采用人工解译方式提取了5000个变化像素和5000个未变化像素来构建变化参考样本集。评价指标包括总体精度、错检率、漏检率、Kappa系数及运行时间，2组实验的精度评价结果如表1、表2所示：
表1数据集1精度评价结果
| 方法 | 总体精度/％ | 错检率/％ | 漏检率/％ | Kappa | 运行时间/秒 |
| 本文方法 | 86.3 | 12.5 | 17.8 | 0.678 | 28.2 |
| OB-HMAD | 79.5 | 19.3 | 24.4 | 0.616 | 22.5 |
| CVA-EM | 71.3 | 30.1 | 26.1 | 0.534 | 7.8 |
表2数据集2精度评价结果
| 方法 | 总体精度/％ | 错检率/％ | 漏检率/％ | Kappa | 运行时间/秒 |
| 本文方法 | 84.4 | 15.1 | 17.9 | 0.633 | 27.9 |
| OB-HMAD | 72.2 | 29.9 | 26.5 | 0.537 | 23.2 |
| CVA-EM | 62.1 | 42.5 | 28.9 | 0.497 | 7.5 |
根据精度评价结果，本方法在2组城市场景的变化检测实验中总体精度分别能够达到86.3％和84.4％，明显优于其他两种方法，与目视分析结果一致。OB-HMAD方法的局限主要在于仅没有考虑阴影因素的影响。像素级方法在2组实验中检测精度均低于70％，且显著低于两种对象级变化检测方法。相比数据集1实验，三种方法在数据集2实验中各精度指标均有所下降，尤其是错检率指标更加明显，其原因主要是数据集2中受到的阴影干扰更加突出。尽管如此，由于本文方法在决策融合时结合了阴影补偿，因此相较于其他两种方法波动较小。在相同实验环境下的运行时间方面，像素级方法CVA-EM的运行时间最短，但检测精度最低；本方法由于引入了阴影检测及补偿策略，运行时间略高于OB-HMAD方法，但检测精度显著提高。
5.4.2、对象级精度评价
为进一步分析阴影对检测精度的影响，实验补充了一组未受阴影影响的多时相高分影像，记为数据集3。该组数据为中国上海地区的多光谱SPOT 5影像，采集时间分别为2009年3月和2012年2月，空间分辨为2.5米，图像尺寸为512×512像素，如图9所示。图9中的(a)为遥感影像#5，(b)为遥感影像#6。设定HMAX＝20，求取的尺度参数集合为HNMImin＝{5,6,7,9,10}，OB-HMAD及本方法变化检测结果如图10、11所示。
在此基础上，基于分割区域的变化情况进行对象级定量精度评价。数据集1、2、3的变化参考图如图12中的(a)、(b)、(c)所示，其中白色代表变化区域，黑色代表未变化区域，本方法获得的剧烈变化及显著变化区域统一归为变化区域。
评价指标包括总体精度、错检率、漏检率、Kappa系数，3组实验的精度评价结果如表3所示：
表3对象级精度评价结果
对比像素、对象级两种评价方式，数据集1、数据集2实验中本方法和OB-HMAD在两种评价方式下获得的精度指标基本一致，即受到阴影干扰时本方法显著优于OB-HMAD；另外，由于数据集3实验未受阴影影响，OB-HMAD的错检率显著降低，但总体精度仍低于本文方法；最后，在3组实验中，本方法的检测精度并未因为阴影干扰的程度不同而产生明显波动，具有良好的鲁棒性。
5.5、实验结论
针对高分辨率遥感影像变化检测中由于阴影所导致的“伪变化”，提出了一种结合阴影检测与多尺度融合的对象级高分辨率遥感影像变化检测方法。该方法不仅构建了一种基于尺度间互信息最小化的目标函数以自适应提取尺度参数集合，还提出了一种结合阴影补偿因子的多尺度决策级融合策略，进而实现了对变化强度等级的划分。实验证明，该方法能够有效减少阴影所导致的错检问题，显著提高变化检测精度且具有良好的鲁棒性。特别是，尽管数据集2相比数据集1阴影导致的“伪变化”更加显著，本方法总体精度仅下降了不到2％，而CVA-EM的总精度下降接近10％，OB-HMAD的总体精度也下降了超过7％，从而进一步证明了所提出的阴影补偿策略是必要且有效的。
实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。
Claims (9)
1.结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，包括以下步骤：
(1)基于矢量量化方法对遥感影像进行分类，基于分类结果，将阴影像素视为出界点，采用卡方检验进行阴影检测；
(2)采用多尺度分割算法提取地理对象；
(3)基于尺度间互信息最小的优化目标自适应提取尺度参数，获取多尺度J-image影像序列；
(4)获取步骤(1)中检测出的阴影和步骤(2)中提取的地理对象在不同时相、相同尺度J-image影像序列中的投影，并在此基础上，分析任意地理对象在不同时相、相同尺度J-image影像序列间的相似程度；
(5)通过结合阴影补偿的决策融合算法实现遥感影像变化检测及变化强度等级的划分。
2.根据权利要求1所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，在步骤(1)中，采用卡方检验进行阴影检测的方法如下：
定义卡方分布：
(X-m)TΣ-1(X-m)～χ2(b)
其中，X为随机变量，m和Σ分别为非阴影区域的均值和协方差矩阵，b为自由度，即多光谱影像的波段数目；
在此基础上，采用下式确定阴影像素：
其中，P表示概率，1-γ为置信度，此时卡方值小于的像素为阴影像素。
3.根据权利要求1所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，在步骤(2)中，多尺度分割仅在单一时相影像中进行，再将分割结果根据配准获得的像素间匹配关系直接映射到另一时相影像中，从而确定统一的地理对象集合R＝{R1,R2,...RG}，其中G为分割结果中的对象总数。
4.根据权利要求1所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，在步骤(3)中，自适应提取尺度参数的过程如下：
(31)依据目视解译在影像中选择尺寸较大的地理对象，统计其包含的像素数量为E；
(32)设定尺度参数的最大值HMAX为大于的最小整数；
(33)依据阴影检测及配准结果，选择在两时相影像中均为非阴影的像素集合，利用该像素集合分别在单一时相影像中计算相邻尺度间的互信息为NMIt(s-1,s)；其中t∈{1,2}，代表不同时相影像；s为尺度序号，其最大值取HMAX；
(34)获得双时相影像相邻尺度间的互信息值：
(35)通过迭代计算遍历所有的尺度参数组合，从而自适应提取NMIall最小值对应的尺度参数集合HNMImin＝{H1,H2,...HD}，D为HNMImin中的尺度总数。
5.根据权利要求1所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，在步骤(4)中，设任意地理对象Ri在时相1、2相同尺度J-image中的投影分别为Ri1和Ri2，则Ri1与Ri2间的相似程度：
其中，μi1、μi2、σi1、σi2、σi1i2分别是Ri1与Ri2的均值、标准差、方差和协方差；C1、C2为常数；SSIM越大，则当前尺度下Ri1与Ri2的相似程度越高，发生变化的可能性越小；遍历所有D个尺度，可获得Ri对应多尺度J-image序列下的相似性度量集合SSIMki∈{SSIM1i,SSIM2i,...,SSIMDi}，其中k为尺度参数集合HNMImin中的尺度序号。
6.根据权利要求5所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，在步骤(5)中，首先定义识别框架U＝{SL,MA,UN}，将地理对象划分为剧烈变化类SL、显著变化类MA和未变化类UN，则焦元A包括{SL}，{MA}，{UN}，{SL,MA,UN}；对任意地理对象Ri，建立如下基本概率赋值函数：
mki({SL})＝(1-SSIMki)×0.7×αk×λ
mki({MA})＝(1-SSIMki)×0.3×αk×λ
mki({UN})＝SSIMki×αk×λ
mki({SL,MA,UN})＝1-αk×λ
其中，αk用于描述对k尺度所提供证据的信任程度，k∈(1,2,...D)，λ为阴影补偿因子；
然后采用下式对建立的基本概率赋值函数进行多尺度证据合成：
其中，mi(A)为合成后的基本概率赋值函数，Ak为k尺度下的焦元，Vi为归一化常数：
最后采用以下判别规则对变化强度等级进行划分：
若满足以下条件之一，则地理对象Ri属于剧烈变化类SL：
条件1、mi({SL})＞0.8；
条件2、mi({SL})＞0.6且mi({MA})＞0.3；
若满足以下条件之一，则地理对象Ri属于显著变化类MA：
条件1、mi({MA})＞0.7且mi({SL})＞0.1；
条件2、mi({UN)}＜0.1
若按照上述判别条件地理对象Ri不属于剧烈变化类SL和显著变化类MA，则地理对象Ri属于未变化类UN。
7.根据权利要求6所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，阴影补偿因子λ用于描述阴影对于变化检测的干扰程度，定义如下：
λ＝l1β1+l2β2+l3β3
其中，β1为Ri在双时相影像中均为非阴影的像素比例，β2为Ri在双时相影像中均为阴影，β3为Ri仅在一个时相影像中为阴影的像素比例，l1、l2、l3分别为β1、β2、β3的权重，l1＞l2＞l3且l1+l2+l3＝1。
8.根据权利要求7所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，设定l1＝0.6，l1＝0.3，l3＝0.1。
9.根据权利要求6所述结合阴影补偿与决策融合的遥感影像变化检测方法，其特征在于，将D个尺度的αk统一设定为0.5。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810968569.5A CN109360184A (zh) | 2018-08-23 | 2018-08-23 | 结合阴影补偿与决策融合的遥感影像变化检测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810968569.5A CN109360184A (zh) | 2018-08-23 | 2018-08-23 | 结合阴影补偿与决策融合的遥感影像变化检测方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109360184A true CN109360184A (zh) | 2019-02-19 |
Family
ID=65349874
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810968569.5A Pending CN109360184A (zh) | 2018-08-23 | 2018-08-23 | 结合阴影补偿与决策融合的遥感影像变化检测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109360184A (zh) |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109859219A (zh) * | 2019-02-26 | 2019-06-07 | 江西理工大学 | 结合相位与光谱的高分遥感影像分割方法 |
| CN111340761A (zh) * | 2020-02-18 | 2020-06-26 | 南京信息工程大学 | 基于分形属性和决策融合的遥感影像变化检测方法 |
| CN111639618A (zh) * | 2020-06-08 | 2020-09-08 | 中国石油大学(华东) | 一种全极化sar影像变化区域精确提取方法 |
| CN113838078A (zh) * | 2021-09-06 | 2021-12-24 | 中国矿业大学(北京) | 采煤塌陷地裂缝的识别与提取方法、装置及存储介质 |
| CN113963222A (zh) * | 2021-10-28 | 2022-01-21 | 中国电子科技集团公司第五十四研究所 | 一种基于多策略组合的高分辨率遥感影像变化检测方法 |
| CN115410096A (zh) * | 2022-11-03 | 2022-11-29 | 成都国星宇航科技股份有限公司 | 卫星遥感影像多尺度融合变化检测方法、介质及电子装置 |
| CN117152619A (zh) * | 2023-10-27 | 2023-12-01 | 广州蓝图地理信息技术有限公司 | 基于高分辨率建筑物遥感影像数据的优化训练方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040120857A1 (en) * | 2002-09-26 | 2004-06-24 | Honeywell Federal Manufacturing & Technologies, Llc | System and method for identifying, reporting, and evaluating presence of substance |
| CN103632363A (zh) * | 2013-08-27 | 2014-03-12 | 河海大学 | 基于多尺度融合的对象级高分辨率遥感影像变化检测方法 |
-
2018
- 2018-08-23 CN CN201810968569.5A patent/CN109360184A/zh active Pending
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040120857A1 (en) * | 2002-09-26 | 2004-06-24 | Honeywell Federal Manufacturing & Technologies, Llc | System and method for identifying, reporting, and evaluating presence of substance |
| CN103632363A (zh) * | 2013-08-27 | 2014-03-12 | 河海大学 | 基于多尺度融合的对象级高分辨率遥感影像变化检测方法 |
Non-Patent Citations (2)
| Title |
|---|
| CHAO WANG 等: "Object-oriented change detection approach for high-resolution remote sensing images based on multiscale fusion", 《JOURNAL OF APPLIED REMOTE SENSING》 * |
| 王超 等: "一种结合阴影补偿的城市高分遥感影像分割方法", 《电子测量与仪器学报》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109859219A (zh) * | 2019-02-26 | 2019-06-07 | 江西理工大学 | 结合相位与光谱的高分遥感影像分割方法 |
| CN111340761A (zh) * | 2020-02-18 | 2020-06-26 | 南京信息工程大学 | 基于分形属性和决策融合的遥感影像变化检测方法 |
| CN111639618A (zh) * | 2020-06-08 | 2020-09-08 | 中国石油大学(华东) | 一种全极化sar影像变化区域精确提取方法 |
| CN111639618B (zh) * | 2020-06-08 | 2021-04-06 | 中国石油大学(华东) | 一种全极化sar影像变化区域精确提取方法 |
| CN113838078A (zh) * | 2021-09-06 | 2021-12-24 | 中国矿业大学(北京) | 采煤塌陷地裂缝的识别与提取方法、装置及存储介质 |
| CN113838078B (zh) * | 2021-09-06 | 2023-06-30 | 中国矿业大学(北京) | 采煤塌陷地裂缝的识别与提取方法、装置及存储介质 |
| CN113963222A (zh) * | 2021-10-28 | 2022-01-21 | 中国电子科技集团公司第五十四研究所 | 一种基于多策略组合的高分辨率遥感影像变化检测方法 |
| CN113963222B (zh) * | 2021-10-28 | 2022-09-02 | 中国电子科技集团公司第五十四研究所 | 一种基于多策略组合的高分辨率遥感影像变化检测方法 |
| CN115410096A (zh) * | 2022-11-03 | 2022-11-29 | 成都国星宇航科技股份有限公司 | 卫星遥感影像多尺度融合变化检测方法、介质及电子装置 |
| CN117152619A (zh) * | 2023-10-27 | 2023-12-01 | 广州蓝图地理信息技术有限公司 | 基于高分辨率建筑物遥感影像数据的优化训练方法 |
| CN117152619B (zh) * | 2023-10-27 | 2024-02-09 | 广州蓝图地理信息技术有限公司 | 基于高分辨率建筑物遥感影像数据的优化训练方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109360184A (zh) | 结合阴影补偿与决策融合的遥感影像变化检测方法 | |
| CN103632363B (zh) | 基于多尺度融合的对象级高分辨率遥感影像变化检测方法 | |
| CN103971115B (zh) | 一种基于NDVI和PanTex指数的新增建设用地图斑自动提取方法 | |
| CN104751478B (zh) | 一种基于多特征融合的面向对象的建筑物变化检测方法 | |
| CN105956557B (zh) | 一种面向对象的时序遥感影像云覆盖区域自动检测方法 | |
| CN107330875B (zh) | 基于遥感图像正反向异质性的水体周边环境变化检测方法 | |
| CN101840581B (zh) | 一种从卫星遥感影像中提取建筑物轮廓的方法 | |
| CN103578110B (zh) | 基于灰度共生矩阵的多波段高分辨率遥感影像分割方法 | |
| CN105335966B (zh) | 基于局域同质性指标的多尺度遥感影像分割方法 | |
| CN105787937B (zh) | 一种基于osm的高分辨遥感影像道路变化检测方法 | |
| CN104361590B (zh) | 一种控制点自适应分布的高分辨率遥感影像配准方法 | |
| CN108596103A (zh) | 基于最佳光谱指数选择的高分辨率卫星遥感影像建筑物提取方法 | |
| CN104915672B (zh) | 一种基于高分辨率遥感图像的矩形建筑物提取方法及系统 | |
| CN103839267B (zh) | 一种基于形态学建筑物指数的建筑物提取方法 | |
| CN110569751B (zh) | 一种高分遥感影像建筑物提取方法 | |
| CN110309780A (zh) | 基于bfd-iga-svm模型的高分辨率影像房屋信息快速监督识别 | |
| CN107230197B (zh) | 基于卫星云图和rvm的热带气旋客观定强方法 | |
| CN104361589A (zh) | 一种基于尺度间映射的高分辨率遥感影像分割方法 | |
| CN110097101A (zh) | 一种基于改进可靠性因子的遥感图像融合与海岸带分类方法 | |
| Yue et al. | Multiscale roughness measure for color image segmentation | |
| CN106128121A (zh) | 基于局部特征分析的车辆排队长度快速检测算法 | |
| CN109871875A (zh) | 一种基于深度学习的建筑物变化检测方法 | |
| CN103077515A (zh) | 一种多光谱图像建筑物变化检测方法 | |
| CN106340005A (zh) | 基于尺度参数自动优化的高分遥感影像非监督分割方法 | |
| CN109859219A (zh) | 结合相位与光谱的高分遥感影像分割方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB02 | Change of applicant information |
Address after: 210032 No. 219 Ning six road, Jiangbei new district, Nanjing, Jiangsu Applicant after: Nanjing University of Information Science and Technology Address before: 211500 Yuting Square, 59 Wangqiao Road, Liuhe District, Nanjing City, Jiangsu Province Applicant before: Nanjing University of Information Science and Technology |
|
| CB02 | Change of applicant information | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190219 |
|
| RJ01 | Rejection of invention patent application after publication |