留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于混合像元分解与EM算法的中低分辨率遥感影像变化检测

吴柯 何坦 杨叶涛

吴柯, 何坦, 杨叶涛. 基于混合像元分解与EM算法的中低分辨率遥感影像变化检测[J]. 武汉大学学报 ● 信息科学版, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
引用本文: 吴柯, 何坦, 杨叶涛. 基于混合像元分解与EM算法的中低分辨率遥感影像变化检测[J]. 武汉大学学报 ● 信息科学版, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
WU Ke, HE Tan, YANG Yetao. Change Detection Method Based on Pixel Unmixing and EM Algorithm for Low and Medium Resolution Remote Sensing Imagery[J]. Geomatics and Information Science of Wuhan University, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
Citation: WU Ke, HE Tan, YANG Yetao. Change Detection Method Based on Pixel Unmixing and EM Algorithm for Low and Medium Resolution Remote Sensing Imagery[J]. Geomatics and Information Science of Wuhan University, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546

基于混合像元分解与EM算法的中低分辨率遥感影像变化检测

doi: 10.13203/j.whugis20160546
基金项目: 

国家自然科学基金 61372153

中央高校基本科研业务费专项 CUGL140410

中央高校基本科研业务费专项 26420160125

详细信息
    作者简介:

    吴柯, 博士, 副教授, 主要从事遥感影像信息处理、遥感地质、变化检测方面的研究。tingke2000@126.com

  • 中图分类号: P237

Change Detection Method Based on Pixel Unmixing and EM Algorithm for Low and Medium Resolution Remote Sensing Imagery

Funds: 

The National Natural Science Foundation of China 61372153

the Fundamental Research Funds for the Central Universities CUGL140410

the Fundamental Research Funds for the Central Universities 26420160125

More Information
    Author Bio:

    WU Ke, PhD, associate professor, specializes in remote sensing image processing, remote sensing geology, and change detection. E-mail: tingke2000@126.com

  • 摘要: 中低分辨率遥感影像中广泛存在的混合像元极大地限制了变化检测结果的精度。基于混合像元分解技术,能够深入到像元内部,比较不同端元的组成分差异影像,然后获取亚像元级别的变化信息。如何从差异影像中确定合适的变化阈值,从而准确地判断变化是否发生,是一个难点问题。在高斯模型分布假设的情况下,采用最大期望法(expectation maximization,EM)自动提取最佳阈值,完成自适应的变化检测过程。选择了两种典型的阈值选择方法与该方法进行比较,结果证明基于EM的自适应变化检测方法可以更准确地提取变化信息,具有较好的稳健性。
  • 图  1  基于EM的阈值选择算法

    Figure  1.  The Threshold Selection Scheme Based on EM Algorithm

    图  2  原始影像TM1和TM2以及模拟影像TM3

    Figure  2.  Original Images TM1 and TM2 and Simulated Image TM3

    图  3  水体、植被和城区的组成分及对应的差值图像

    Figure  3.  The Fractional Images of Water, Vegetation and Urban City and the Corresponding Difference Images

    图  4  3种地物差异影像的变化曲线

    Figure  4.  The Change Curves of Three Fractional Difference Images

    图  5  3种方法真实数据变化检测的结果比较

    Figure  5.  Comparison of Three Methods for Real Data

    图  6  3种方法在TM影像的两个局部区域比较

    Figure  6.  Comparison of Two Parts of the TM Images Based on Three Different Methods

    图  7  3种方法模拟数据变化检测的结果比较

    Figure  7.  Comparison of Three Methods for Simulated Data

    表  1  3种方法的变化率比较

    Table  1.   Comparison of Rate of Change Based on Three Methods

    类别 方法1 方法2 本文方法
    阈值 变化率/% 阈值 变化率/% 阈值 变化率/%
    水体 0.21 1.30 0.25 0.90 0.28 0.50
    植被 0.09 3.90 0.11 3.30 0.13 3.10
    城区 0.21 6.12 0.29 3.91 0.37 2.41
    总体变化率/% 11.32 8.11 6.01
    下载: 导出CSV

    表  2  3种方法模拟数据的定量比较

    Table  2.   The Quantitative Comparison of Three Methods for Simulated Data

    方法 类别 TR/% FAR/%
    水体 94.3 5.7
    方法1 植被 84.0 6.0
    建筑物 52.1 47.9
    总体 71.928.1
    水体 99.9 0.1
    方法2 植被 86.9 13.1
    建筑物 70.9 29.1
    总体 84.0 16.0
    水体 100 0
    本文方法 植被 95.2 4.8
    建筑物 76.9 23.1
    总体 89.5 10.5
    下载: 导出CSV
  • [1] 杜培军, 柳思聪.融合多特征的遥感影像变化检测[J].遥感学报, 2012, 16 (4): 663-677 http://d.old.wanfangdata.com.cn/Periodical/ygxb201204001

    Du Peijun, Liu Sicong. Change Detection from Multi-temporal Remote Sensing Images by Integrating Multiple Features[J]. Journal of Remote Sensing, 2012, 16(4): 663-677 http://d.old.wanfangdata.com.cn/Periodical/ygxb201204001
    [2] Howarth P J, Wickware G M. Procedures for Change Detection Using Landsat Digital Data[J]. International Journal of Remote Sensing, 1981, 2(3): 191-277 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1080/01431168108948362
    [3] Mas J F. Monitoring Land-Cover Changes: A Comparison of Change Detection Techniques[J]. International Journal of Remote Sensing, 1999, 20(1):139-152 doi:  10.1080/014311699213659
    [4] 陈雪, 马建文, 戴芹.基于贝叶斯网络分类的遥感影像变化检测[J].遥感学报, 2005, 9(6): 667-672 http://d.old.wanfangdata.com.cn/Periodical/ygxb200506007

    Chen Xue, Ma Jianwen, Dai Qin. Remote Sensing Change Detection Based on Bayesian Networks Classifications[J]. Journal of Remote Sensing, 2005, 9(6): 667-672 http://d.old.wanfangdata.com.cn/Periodical/ygxb200506007
    [5] 武辰, 杜博, 张良培.基于独立成分分析的高光谱变化检测[J].遥感学报, 2012, 16(3): 545-561 http://d.old.wanfangdata.com.cn/Periodical/ygxb201203009

    Wu Chen, Du Bo, Zhang Liangpei. Hyperspectral Change Detection Based on Independent Component Analysis[J]. Journal of Remote Sensing, 2012, 16(3): 545-561 http://d.old.wanfangdata.com.cn/Periodical/ygxb201203009
    [6] 李松, 安裕伦, 华厚强.基于遥感变化检测的石漠化信息自动提取[J].遥感技术与应用, 2012, 27(1): 149-153 http://d.old.wanfangdata.com.cn/Periodical/ygjsyyy201201021

    Li Song, An Yulun, Hua Houqiang. Automated Method Based on Change Detection for Extracting Karst Rock Desertification Information Using Remote Sensing[J]. Remote Sensing Technology and Application, 2012, 27(1): 149-153 http://d.old.wanfangdata.com.cn/Periodical/ygjsyyy201201021
    [7] 魏立飞, 钟燕飞, 张良培, 等.多波段信息融合的遥感影像变化检测[J].武汉大学学报·信息科学版, 2014, 39(1): 8-11 http://ch.whu.edu.cn/CN/abstract/abstract2859.shtml

    Wei Lifei, Zhong Yanfei, Zhang Liangpei, et al. Remote Sensing Image Change Detection Based on Multi-band Information Fusion[J]. Geomatics and Information Science of Wuhan University, 2014, 39(1): 8-11 http://ch.whu.edu.cn/CN/abstract/abstract2859.shtml
    [8] Tong S S, Pham T L, Pham V C. Land Cover Change Analysis Using Change Vector Analysis Method in Duy Tien District, Ha Nam Province in Vietnam[C]. 7th FIG Regional Conference TS 1G-Remote Sensing for Sustainable Development, Hanoi, Vietnam, 2009
    [9] Turner B L, Meyer W B, Skole D L. Global Land-Use/Land-Cover Change: Towards an Integrated Study[J]. Ambio, 1994, 23(1): 91-95
    [10] Haertel V, Shimabukuro Y E, Almeida-Filho R. Fraction Images in Multitemporal Change Detection[J]. International Journal of Remote Sensing, 2004, 25(23):5473-5489 doi:  10.1080/01431160412331269751
    [11] Anderson L O, Shimabukuro Y E, Arai E. Cover: Multitemporal Fraction Images Derived from Terra MODIS Data for Analysing Land Cover Change over the Amazon Region[J]. International Journal of Remote Sensing, 2005, 26(11):2251-2257 doi:  10.1080/01431160310001620795
    [12] Lu D, Batistella M, Moran E. Multitemporal Spectral Mixture Analysis for Amazonian Land-Cover Change Detection[J]. Canadian Journal of Remote Sensing, 2004, 30(1):87-100 doi:  10.5589/m03-055
    [13] 赵春晖, 齐滨, 王玉磊.一种改进的N-FINDR高光谱端元提取算法[J].电子与信息学报, 2012, 34(2):500-503 http://d.old.wanfangdata.com.cn/Periodical/dzkxxk201202041

    Zhao Chunhui, Qi Bin, Wang Yulei. An Improved N-FINDR Hyperspectral Endmember Extraction Algorithm[J]. Journal of Electronics and Information Technology, 2012, 34(2):500-503 http://d.old.wanfangdata.com.cn/Periodical/dzkxxk201202041
    [14] Heinz D C, Chang C. Fully Constrained Least Squares Linear Spectral Mixture Analysis Method for Material Quantification in Hyperspectral Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 39(3):529-545
    [15] Yang H L, Peng J H, Xia B R, et al. An Improved EM Algorithm for Remote Sensing Classification[J]. Chinese Science Bulletin, 2013, 58(9): 1060-1071 doi:  10.1007/s11434-012-5485-4
    [16] 吴柯, 牛瑞卿, 王毅.基于PCA与EM算法的多光谱遥感影像变化检测研究[J].计算机科学, 2010, 37(3): 275-282 http://d.old.wanfangdata.com.cn/Periodical/jsjkx201003071

    Wu Ke, Niu Ruiqing, Wang Yi. Change Detection of Multi-spectral Remote Sensed Images Based on PCA and EM Algorithm[J].Computer Science, 2010, 37(3): 275-282 http://d.old.wanfangdata.com.cn/Periodical/jsjkx201003071
    [17] 佃袁勇.基于遥感影像的变化检测研究[D].武汉: 武汉大学, 2005

    Dian Yuanyong. Analysis of Change Detection Based on Remote Sensed Images[D]. Wuhan: Wuhan University, 2005
    [18] 李雪, 舒宁, 王琰, 等.利用土地利用状态转移分析的变化检测[J].武汉大学学报·信息科学版, 2011, 36(8): 951-952 http://ch.whu.edu.cn/CN/abstract/abstract618.shtml

    Li Xue, Shu Ning, Wang Yan, et al. Change Detection Based on Land-Use Status Transition Analysis[J]. Geomatics and Information Science of Wuhan University, 2011, 36(8): 951-952 http://ch.whu.edu.cn/CN/abstract/abstract618.shtml
    [19] 吴炜, 沈占锋, 吴田军, 等.联合概率密度空间的遥感自适应变化检测方法[J].测绘学报, 2016, 45(1):74-79 http://d.old.wanfangdata.com.cn/Periodical/chxb201601011

    Wu Wei, Shen Zhanfeng, Wu Tianjun, et al. Joint Probability Space Based Self-Adaptive Remote Sen- sing Change Detection Method[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(1):74-79 http://d.old.wanfangdata.com.cn/Periodical/chxb201601011
  • [1] 季顺平, 田思琦, 张驰.  利用全空洞卷积神经元网络进行城市土地覆盖分类与变化检测 . 武汉大学学报 ● 信息科学版, 2020, 45(2): 233-241. doi: 10.13203/j.whugis20180481
    [2] 眭海刚, 冯文卿, 李文卓, 孙开敏, 徐川.  多时相遥感影像变化检测方法综述 . 武汉大学学报 ● 信息科学版, 2018, 43(12): 1885-1898. doi: 10.13203/j.whugis20180251
    [3] 魏立飞, 钟燕飞, 张良培, 李平湘.  多波段信息融合的遥感影像变化检测 . 武汉大学学报 ● 信息科学版, 2014, 39(1): 8-11.
    [4] 王楠, 张良培, 杜博.  最小光谱相关约束NMF的高光谱遥感图像混合像元分解 . 武汉大学学报 ● 信息科学版, 2014, 39(1): 22-26.
    [5] 佃袁勇, 方圣辉, 姚崇怀.  一种面向地理对象的遥感影像变化检测方法 . 武汉大学学报 ● 信息科学版, 2014, 39(8): 906-912. doi: 10.13203/j.whugis20130053
    [6] 申邵洪, 万幼川, 龚浩, 赖祖龙.  遥感影像变化检测自适应阈值分割的Kriging方法 . 武汉大学学报 ● 信息科学版, 2009, 34(8): 902-905.
    [7] 董明, 张海涛, 祝晓坤, 陈春希.  基于遥感影像的地图道路网数据变化检测研究 . 武汉大学学报 ● 信息科学版, 2009, 34(2): 178-182.
    [8] 万幼川, 申邵洪, 张景雄.  基于概率统计模型的遥感影像变化检测 . 武汉大学学报 ● 信息科学版, 2008, 33(7): 669-672.
    [9] 尤红建, 傅琨.  基于分布模型差异的SAR变化检测 . 武汉大学学报 ● 信息科学版, 2008, 33(5): 454-456.
    [10] 袁修孝, 宋妍.  基于边缘特征匹配的遥感影像变化检测预处理方法 . 武汉大学学报 ● 信息科学版, 2007, 32(5): 381-384.
    [11] 袁修孝, 季顺平.  POS数据辅助的航空影像变化检测方法研究 . 武汉大学学报 ● 信息科学版, 2007, 32(4): 283-286.
    [12] 马国锐, 眭海刚, 李平湘, 秦前清.  一种遥感影像核变化检测方法 . 武汉大学学报 ● 信息科学版, 2007, 32(7): 597-600.
    [13] 张晓东, 李德仁, 龚健雅, 秦前清.  遥感影像与GIS分析相结合的变化检测方法 . 武汉大学学报 ● 信息科学版, 2006, 31(3): 266-269.
    [14] 江利明, 廖明生, 张路, 林珲.  顾及空间邻域关系的多时相SAR影像变化检测 . 武汉大学学报 ● 信息科学版, 2006, 31(4): 312-315.
    [15] 徐丽华, 江万寿.  顾及配准误差的遥感影像变化检测 . 武汉大学学报 ● 信息科学版, 2006, 31(8): 687-690.
    [16] 苏国中, 张剑清, 陈炳贵.  GIS知识引导的正射影像变化检测及更新 . 武汉大学学报 ● 信息科学版, 2005, 30(8): 664-667.
    [17] 李全, 李霖, 赵曦.  基于Landsat TM影像的城市变化检测研究 . 武汉大学学报 ● 信息科学版, 2005, 30(4): 351-354.
    [18] 张路, 廖明生, 盛辉.  基于正交变换的多通道遥感影像变化检测 . 武汉大学学报 ● 信息科学版, 2004, 29(5): 456-460,469. doi: 10.13203/j.whugis2004.05.018
    [19] 李德仁.  利用遥感影像进行变化检测 . 武汉大学学报 ● 信息科学版, 2003, 28(S1): 7-12.
    [20] 方针, 张剑清, 张祖勋.  基于城区航空影像的变化检测 . 武汉大学学报 ● 信息科学版, 1997, 22(3): 240-244.
  • 加载中
图(7) / 表(2)
计量
  • 文章访问数:  334
  • HTML全文浏览量:  15
  • PDF下载量:  168
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-09-18
  • 刊出日期:  2019-04-05

基于混合像元分解与EM算法的中低分辨率遥感影像变化检测

doi: 10.13203/j.whugis20160546
    基金项目:

    国家自然科学基金 61372153

    中央高校基本科研业务费专项 CUGL140410

    中央高校基本科研业务费专项 26420160125

    作者简介:

    吴柯, 博士, 副教授, 主要从事遥感影像信息处理、遥感地质、变化检测方面的研究。tingke2000@126.com

  • 中图分类号: P237

摘要: 中低分辨率遥感影像中广泛存在的混合像元极大地限制了变化检测结果的精度。基于混合像元分解技术,能够深入到像元内部,比较不同端元的组成分差异影像,然后获取亚像元级别的变化信息。如何从差异影像中确定合适的变化阈值,从而准确地判断变化是否发生,是一个难点问题。在高斯模型分布假设的情况下,采用最大期望法(expectation maximization,EM)自动提取最佳阈值,完成自适应的变化检测过程。选择了两种典型的阈值选择方法与该方法进行比较,结果证明基于EM的自适应变化检测方法可以更准确地提取变化信息,具有较好的稳健性。

English Abstract

吴柯, 何坦, 杨叶涛. 基于混合像元分解与EM算法的中低分辨率遥感影像变化检测[J]. 武汉大学学报 ● 信息科学版, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
引用本文: 吴柯, 何坦, 杨叶涛. 基于混合像元分解与EM算法的中低分辨率遥感影像变化检测[J]. 武汉大学学报 ● 信息科学版, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
WU Ke, HE Tan, YANG Yetao. Change Detection Method Based on Pixel Unmixing and EM Algorithm for Low and Medium Resolution Remote Sensing Imagery[J]. Geomatics and Information Science of Wuhan University, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
Citation: WU Ke, HE Tan, YANG Yetao. Change Detection Method Based on Pixel Unmixing and EM Algorithm for Low and Medium Resolution Remote Sensing Imagery[J]. Geomatics and Information Science of Wuhan University, 2019, 44(4): 555-562. doi: 10.13203/j.whugis20160546
  • 遥感影像的变化检测一直是监测和管理自然资源和城市发展的重要手段之一[1],主要有3种模式。一是遥感影像直接比较法,例如图像比值法[2]、归一化植被指数(normalized difference vegetation index,NDVI)[3]等;二是分类后比较变化检测法,首先对两个时相的影像运用统计表的分类体系各自做分类,再对分类后的结果比较,进而得出变化的区域范围,如贝叶斯网络[4]、支持向量机(support vector machine,SVM)[5]、面向对象分类[6]等;三是组合分析法,比较流行的是光谱向量分析(change vector analysis,CVA)和马尔可夫随机场模型(Markov random field,MRF)[7-8]两种代表性方法。前者是在相对辐射归一化校正的基础上进行的,运用一个向量来描述从第一时相到第二时相的光谱变化的方向和量级;后者利用融合所有波段的变化信息,通过迭代计算得到最优的MRF模型,获得最完整准确的检测结果。

    由于空间分辨率的不足,影像中普遍存在的混合像元(即单个像元中包含多种地物类别)成为变化检测技术发展的一大障碍。传统的检测方法基本上是以单个像元作为对象,无法深入到像元内部去完成亚像元级的检测。高分辨率的遥感数据虽然能够更精细地判断混合像元内部的分布,但是考虑到运行成本、快速反应等方面的因素,实际应用中的变化检测仍然是以中、低分辨率数据为主[9]。因此,研究者提出了亚像元变化检测模型,将像元级的变化检测转变为亚像元级,深入到像元内部完成变化检测,从而获取更多的变化信息[10]。这类方法主要是采用变化检测的第二种模式,即分类后变化检测法。在这个过程中,如何确定变化阈值一直是一个难点问题,目前的研究工作未有专门的涉及[11]。传统的阈值选择方法一般是根据人工经验设定,即不断地调整参数和评估结果,直到达到满意的检测效果为止。这种方法需要耗费大量的时间,同时没有明确的物理意义;另一种方法是根据差异图像的统计学规律来设定。假定不变的像素集中在差异图像的直方图中心,而变化的像元分布在直方图的两侧,通过计算图像的标准差σ,以2σ作为变化的阈值[12]。这种方法主观性较强,没有自适应的调整过程且容易发生错判,可能导致不稳定的变化检测结果。

    鉴于以上背景,本文引入一种最大期望算法(expectation maximization,EM)来确定变化阈值。假设差异图像满足高斯统计分布,可以基于最小错误概率的贝叶斯决策方法来区别未变化与变化的类别。高斯分布模型的参数通过EM迭代算法自动获取,满足变化和未变化类别的数学期望均达到最大时,定义这两类的条件概率相等,从而求解阈值。本文通过与两种传统阈值挑选方法进行比较,证明EM自动提取阈值的算法能够获得最佳的检测效果,有效提高亚像元变化检测的精度。

    • 基于混合像元分解的变化检测主要可以分为以下3个步骤。

      1) 提取影像的目标端元类别。由于N-FINDR算法具有简单的结构及良好的提取效果,本文选用其来提取影像的目标端元[13]。该算法认为遥感图像中所有的像元光谱构成一个凸面点集,端元位于凸面体的顶点,混合像元则位于凸面体表面或凸面体内部。在该凸面体的点集中,由凸面体顶点(端元)组成的凸面体被认为是体积最大的。因此,寻找目标端元实际上就是寻找凸面体体积最大的顶点。在算法的实际操作中,随机选择初始端元的方式大大增加了运算时间,且容易造成算法无法收敛。因此,本文利用人工提取初始端元的方式来代替随机选择。假设人工提取的端元集合为p,需要从中挑选k(k < p)个端元构成初始端元集合E,可表示为:

      $$ \mathit{\boldsymbol{E}} = \left[ {\begin{array}{*{20}{c}} 1&1& \cdots &1\\ {{e_1}}&{{e_2}}& \cdots &{{e_k}} \end{array}} \right] $$ (1)

      式中,ei是第i个端元向量,i=1~k。凸面体体积V的计算公式为:

      $$ V\left( \mathit{\boldsymbol{E}} \right) = \frac{1}{{\left( {k - 1} \right)!}}{\rm{abs}}\left( {\left| \mathit{\boldsymbol{E}} \right|} \right) $$ (2)

      式中,abs( )表示取绝对值。

      为了获取最纯净的端元,每一次迭代都从集合p中选择新的k个端元到式(2)中计算,体积大的端元被选中并保留下来。整个过程不需要任何参数输入,直到找到最大的体积为止。这时,凸面体的顶点被认为是所提取的端元。

      2) 计算目标端元所对应的组份百分比[14]。依据线性光谱混合理论,影像中任意一个像元的光谱是所包含的各端元光谱的线性组合。设m表示波段数,k表示选定的端元组成分的个数,r是大小为m×1的像元光谱矢量,S=[s1 s2sk]是大小为m×k的端元光谱矩阵,α=[α1 α2αk]Tk维矢量,其各分量元素为对应像元组分,nm维随机噪声,且满足n~N(0, 1)的独立分布。线性混合像元模型为:

      $$ \mathit{\boldsymbol{r}} = \mathit{\boldsymbol{S\alpha }} + \mathit{\boldsymbol{n}} $$ (3)

      在端元光谱矩阵S已知的情况下,按照最小二乘算法得出组分α的最优估计:

      $$ \mathit{\boldsymbol{\hat \alpha = }}{\left( {{\mathit{\boldsymbol{S}}^{\rm{T}}}\mathit{\boldsymbol{S}}} \right)^{ - 1}}{\mathit{\boldsymbol{S}}^{\rm{T}}}\mathit{\boldsymbol{r}} $$ (4)

      实际应用中,各组份值必须加两个条件的限制,即:

      $$ \sum\limits_{j = 1}^k {{\alpha _j}} = 1\;及\;{\alpha _j} \ge 0 $$ (5)

      2) 构造差异影像。在获取不同端元组份的基础上,对于变化前后(假设两期时间表示为T1T2)的组成分影像进行差值运算。不同的端元可以构造出不同的差异影像,例如第j类端元的组份的差异影像为:

      $$ {x_k} = {\alpha _{j,T1}} - {\alpha _{j,T2}} $$ (6)

      式中,xk>0表示第k类端元的组成分比例增加,xk < 0表示组成分比例减少。由于本文的研究只关注获取变化的区域,并不需要了解具体变化类别的趋势,因此,可将式(6)修改为:

      $$ {x_k} = \left| {{\alpha _{k,T1}} - {\alpha _{k,T2}}} \right| $$ (7)
    • 首先,把获取的差异影像分成未变化和变化两类,分别用wnwc表示。贝叶斯公式实质是通过观测差异影像中的像元值x把状态的先验概率p(wi)转化为状态的条件概率p(wi|x):

      $$ p\left( {{w_i}\left| x \right.} \right) = \frac{{p\left( {{w_i}} \right)p\left( {x\left| {{w_i}} \right.} \right)}}{{p\left( x \right)}} $$ (8)

      式中,in, cp(wi)是wi的先验概率分布;p(x|wi)是差异影像中的像元xwi中出现的似然概率;p(x)是全概率函数,由未变化和变化两类概率密度函数混合组成:

      $$ p\left( x \right) = p\left( {x\left| {{w_n}} \right.} \right)p\left( {{w_n}} \right) + p\left( {x\left| {{w_c}} \right.} \right)p\left( {{w_c}} \right) $$ (9)

      对于未变化的区域,差异影像是不同时相的组成分比例的微小差异和一些随机噪声。这种差异近似等于0,噪声近似服从高斯分布;而对于发生变化的区域,由于两时相的丰度影像近似服从高斯分布,根据中心极限定理,差异影像中的对应像元也被认为近似服从高斯分布。因此,似然概率密度函数可以表示为:

      $$ p\left( {x\left| {{w_i}} \right.} \right) = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}\sigma _i^2} }}\exp \left\{ { - \frac{{{{\left( {x - {m_i}} \right)}^2}}}{{2\sigma _i^2}}} \right\} $$ (10)

      式中,mi为均值;σ2为方差。未变化的像元类wn和变化的像元类wc的均值和方差分别用mnσn2mcσc2来表示。

      然后,基于最小错误率的贝叶斯决策规则,比较条件概率p(wi|x)的大小。若p(wn|x)>p(wc|x),则把x归类于未变化类wn;若p(wn|x) < p(wc|x),则把x归类于变化类wc;若p(wn|x)=p(wc|x),则认为x是未变化类与变化类两个集合的交点,记为差值影像的最佳阈值T,因此,利用式(8)将该等式变化为:

      $$ p\left( {{w_n}} \right)p\left( {T\left| {{w_n}} \right.} \right) = p\left( {{w_c}} \right)p\left( {T\left| {{w_c}} \right.} \right) $$ (11)

      该方程中的参数需要从差异影像中估计得到。

    • EM算法是一个迭代的过程,无需任何外来数据和先验知识,仅从观测数据本身得到参数的估计值。每次迭代由求期望值和期望最大化两个步骤组成,前者根据待估计参数的当前值,从观测数据中直接估计概率密度的期望值,后者通过最大化这一期望来更新参数的估计量,这两步在整个迭代过程中依次交替进行,直至迭代过程收敛[15]。以未变化像元wn为例,采用EM算法来估计wn的统计分布参数[16]

      $$ {p^{t + 1}}\left( {{w_n}} \right) = \frac{{\sum\limits_x {\frac{{{p^t}\left( {{w_n}} \right){p^t}\left( {x\left| {{w_n}} \right.} \right)}}{{{p^t}\left( x \right)}}} }}{{I \times J}} $$ (12)
      $$ m_n^{t + 1}\left( {{w_n}} \right) = \frac{{\sum\limits_x {\frac{{{p^t}\left( {{w_n}} \right){p^t}\left( {x\left| {{w_n}} \right.} \right)}}{{{p^t}\left( x \right)}}x} }}{{\sum\limits_x {\frac{{{p^t}\left( {{w_n}} \right){p^t}\left( {x\left| {{w_n}} \right.} \right)}}{{{p^t}\left( x \right)}}} }} $$ (13)
      $$ \begin{array}{*{20}{c}} {{{\left( {{\sigma ^2}} \right)}^{t + 1}}\left( {{w_n}} \right) = }\\ {\frac{{\sum\limits_x {\frac{{{p^t}\left( {{w_n}} \right){p^t}\left( {x\left| {{w_n}} \right.} \right)}}{{{p^t}\left( x \right)}}{{\left( {x - m_n^t} \right)}^2}} }}{{\sum\limits_x {\frac{{{p^t}\left( {{w_n}} \right){p^t}\left( {x\left| {{w_n}} \right.} \right)}}{{{p^t}\left( x \right)}}} }}} \end{array} $$ (14)

      式(12)~(14)估计的分别是先验概率pt+1(wn)、均值mnt+1(wn)

      和方差(σ2)t+1(wn)。式中,tt+1分别代表了当前和下一次的迭代次数;IJ分别代表了差异影像的行数和列数;x表示差异影像中的像元值;p(x|wn)是类条件概率;p(x)是全概率。同理,变化类wc的相关参数先验概率p(wc)、均值mc和方差σc2的估计方法同上。

      根据差异影像的自身特点,本文采用了一种选取合理阈值以构造高置信度像元子集,从中估计待定参数初始值的策略[17]。在差异影像的直方图中,未变化的像元值较小,主要集中在直方图左边界附近,而变化类别的像元值较大,分布在靠近直方图右边界的范围上。基于这一特点,可以在直方图上选取左右各一个阈值TnTc。选取方法如式(15)至式(17),其中,λ∈(0, 1)为权重调节因子,md为中间值。

      $$ {T_n} = \left( {1 - \lambda } \right){m_d} $$ (15)
      $$ {T_c} = \left( {1 + \lambda } \right){m_d} $$ (16)
      $$ {m_d} = \frac{1}{2}\left( {\max \left( x \right) - \min \left( x \right)} \right) $$ (17)

      以两个像元子集Sn={x|x < Tn}和Sc={x|x>Tc}分别作为未变化和变化两类像元的初始典型样本集,而没有被这两个子集所包含的像元则被记为无标识样本,可以将该样本记为Su={x|TnxTc}。然后分别在这两个子集上按式(18)至式(20)计算先验概率、均值以及方差的初始值:

      $$ p\left( {{w_n}} \right) = \frac{{\left\| {{S_n}} \right\|}}{{IJ}} $$ (18)
      $$ {m_n} = \frac{{\sum\limits_{x \in {S_n}} x }}{{\left\| {{S_n}} \right\|}} $$ (19)
      $$ \sigma _n^2 = \frac{{\sum\limits_{x \in {S_n}} {{{\left( {x - {m_n}} \right)}^2}} }}{{\left\| {{S_n}} \right\|}} $$ (20)

      变化像元的先验概率、均值和方差初始值的计算方法同上。得到分布参数的初始估计后,将初始值代入式(12)至式(14)迭代计算,调整参数。整个过程中,无标识类样本子集Su逐步减少,未变化类与变化类的样本子集SnSc则逐步增加。当迭代终止时,这些参数的值达到最大的数学期望。无标识类样本子集变为0,而未变化类与变化类的样本子集拥有一个交点,满足两类条件概率相等,记为最佳阈值T。根据式(11),可以得到关于T的方程:

      $$ \begin{array}{*{20}{c}} {\left( {\sigma _n^2 - \sigma _c^2} \right){T^2} + 2\left( {{m_n}\sigma _c^2 - {m_c}\sigma _n^2} \right)T + m_c^2\sigma _n^2 - }\\ {m_n^2\sigma _c^2 - 2\sigma _c^2\sigma _n^2\ln \left[ {\frac{{{\sigma _n}p\left( {{w_n}} \right)}}{{{\sigma _c}p\left( {{w_c}} \right)}}} \right] = 0} \end{array} $$ (21)

      求解阈值T后,根据以下判断准则得到变化检测结果图:若x>T,则xwc,归属到变化像类,赋值为1,否则,xwn,归属到未变化像元类,赋值为0。整个EM阈值选择法的流程可用图 1表示。

      图  1  基于EM的阈值选择算法

      Figure 1.  The Threshold Selection Scheme Based on EM Algorithm

    • 实验选择广东省深圳市南部地区两个不同时期6个波段的TM影像。两期影像的获取时间分别是2000年9月5日和2002年11月12日,大小是400×400像素。经过严格的几何配准后,将它们分别定义为TM1和TM2,作为实验的真实数据,分别如图 2(a)图 2(b)所示。为了定量比较,本文还构造一个模拟数据,即在原始TM1影像中截取3个10像素×10像素的纯净地块,分别代表水体、植被、建筑3种不同的类别,把它们的分布地块进行交换覆盖,定义为TM3,如图 2(c)所示。模拟数据是通过变换纯净地块得到的,因此,在TM1和TM3两幅影像中只有被交换的覆盖区是变化的,交换区以外的部分完全一致。

      图  2  原始影像TM1和TM2以及模拟影像TM3

      Figure 2.  Original Images TM1 and TM2 and Simulated Image TM3

    • 由于TM1和TM2两幅影像的面积较小,而且间隔的时间不长,因此,可以认为两期影像的纯净端元的数目和类别均保持一致。首先,将数据进行MNF(minmum noise fraction)变换,结合N-维光谱可视化分析,人工挑选出多个纯净像元作为N-FINDR的初始端元集,代入式(1)、式(2)中进行迭代,生成植被、水体及建筑物的3种典型地物端元光谱。然后,对两幅影像分别进行限定性的混合像元分解,并利用式(7)作差值,最终获得水体、植被和城区的3个类别的差异影像。如图 3(a)~3(c)分别代表水体在不同时间的分解结果以及对应的差异影像。同样,图 3(d)~3(f)图 3(g)~3(i)分别代表植被与城区的分解结果及对应的差异影像。组成分图像上,越亮的部分代表该地物所占的组分比例高,反之,代表所占的组分比例低。对比图 3(c)图 3(f)图 3(i)这3个类别的差异图像,可以发现,每一种类别都包含许多变化的组分信息,而且拥有相似的结构特点。除此之外,差异图像中还包含许多伪变化的信息。

      图  3  水体、植被和城区的组成分及对应的差值图像

      Figure 3.  The Fractional Images of Water, Vegetation and Urban City and the Corresponding Difference Images

      采用两种自动提取阈值的方法与本文方法进行对比:以文献[12]的方法计算差异影像中不变像元的标准差,取值2倍作为阈值,记为方法1;以文献[18]中提到的联合概率密度方法,将变化定义在联合概率密度空间,以此作为相似性度量,并采用局部混合信息指标自动确定阈值,记为方法2。在本文方法中,首先求取每个类别变化的组分比例差异图像相应的变化曲线,如图 4(a)~4(c)所示。根据变化曲线的结果分别计算初始阈值TnTc。取λ为0.5,得出类别1~3的初始阈值分别为0.673和2.349、0.385和0.756以及0.623和2.480。然后,根据TnTc选择像元的子集,求取EM算法的4个初始参数。最后,代入式(12)至式(14)进行迭代计算,更新参数直到收敛,计算出最佳阈值。根据阈值对差异影像进行总体分割后得到3个类别的变化区域。

      图  4  3种地物差异影像的变化曲线

      Figure 4.  The Change Curves of Three Fractional Difference Images

      为了方便比较,将3个类别的变化合并在一个图里面显示。图 5(a)~5(c)分别代表利用3种方法得到的变化检测结果图。其中,白色的区域表示变化部分,黑色表示没有变化的部分。目视对比可以发现,3种结果的差异比较明显:方法1检测的变化范围最大,基本覆盖了整个图像,而方法2检测的变化要小一些,本文方法检测出来的变化范围最小,主要覆盖了混合地物比较集中的区域。由于缺乏真实验证数据,通过放大两个局部区域,直观地说明检测效果(图 6)。第一个局部放大区选择原始影像的右上部分,即植被覆盖茂盛的地区,对图 6(a)图 6(b)两个时相影像进行对比观察,可以发现植被减少,变成了裸土与道路。方法1可以检测出变化,但是有许多几何形态各异的离散点,且伪变化较多(图 6(c))。方法2并未检测出变化(图 6(d))。而运用本文提出的方法(图 6(e))可以较好地检测出明显的、完整的变化几何形状。第二个局部放大区选择原图像的左下部分,即水体、城区与植被混杂的地区,从图 6(f)图 6(g)两个时相的影像上可以很清晰地看到沿岸水体所发生的变化。由于河道沿岸的冲刷作用,导致该部分的水体中泥沙含量有明显的增加,造成混合区域明显增大。方法1检测的误差比较大(图 6(h)),除了检测出沿岸水体的变化之外,还夹杂着一些城区和植被的伪变化,方法2能够缓解这种伪变化的发生(图 6(h)),本文方法则较好地反映了真实的变化情况,主要反映出沿岸水体的泥沙增加(图 6(j))。

      图  5  3种方法真实数据变化检测的结果比较

      Figure 5.  Comparison of Three Methods for Real Data

      图  6  3种方法在TM影像的两个局部区域比较

      Figure 6.  Comparison of Two Parts of the TM Images Based on Three Different Methods

      利用状态转移矩阵对两种方法检测的结果统计变化量[18],见表 1。本文方法计算得出的水体、植被及建筑物3类地物的变化阈值分别为0.28、0.13和0.37,均高于前两种检测方法。从表 1的数据可以看出,每类地物的检测变化率中,方法1最大,而本文方法的检测变化率最小,这与目视效果是相同的。综合计算,3种方法的检测总体变化量分别为11.32%、8.11%以及6.01%,本文方法与最大的检测变化结果之间相差近两倍。

      表 1  3种方法的变化率比较

      Table 1.  Comparison of Rate of Change Based on Three Methods

      类别 方法1 方法2 本文方法
      阈值 变化率/% 阈值 变化率/% 阈值 变化率/%
      水体 0.21 1.30 0.25 0.90 0.28 0.50
      植被 0.09 3.90 0.11 3.30 0.13 3.10
      城区 0.21 6.12 0.29 3.91 0.37 2.41
      总体变化率/% 11.32 8.11 6.01
    • 模拟数据选择图 1的TM1和TM3作为两期不同的变化影像,主要观察交换的纯净地块的检测情况。本实验采取与真实数据同样的步骤,依次检测出水体、植被及建筑的变化结果后,合并在同一幅图像里显示,如图 7(a)~7(c)。3种方法均可以完整地将模拟交换的所有纯净地块(方形区域)检测出来,而且没有出现漏检的情况。但是每种方法的检测图像中明显存在一些误检的像元。其中,方法1的误检像元最多,而本文方法的误检像元最少。

      图  7  3种方法模拟数据变化检测的结果比较

      Figure 7.  Comparison of Three Methods for Simulated Data

      定量统计两种方法的检测精度[19]。假设P代表目标物的数目,N代表非目标物的数目,TP表示算法检测为目标物且经验证为目标物的数目,FP是算法检测为目标物而经验证为非目标物的数目(误检)。在此基础上定义正确率TR、误检率FAR两个指标,即:

      $$ \left\{ \begin{array}{l} {\rm{TR}} = \frac{{{\rm{TP}}}}{{{\rm{TP}} + {\rm{FP}}}} \times 100\% \\ {\rm{FAR = }}\frac{{{\rm{FP}}}}{{{\rm{TP + FP}}}} \times 100\% \end{array} \right. $$

      采用上述3种方法实验,以像素为单位,统计这两个指标评价检测精度。表 2中3种方法的检测率分别为71.9%、84.0%及89.5%,与目视结果相符。方法1有最高的误检率,达到28.1%,尤其是对建筑物的误检达到了47.9%;方法2的检测率比前者有所提高;而本文方法对所有变化区域的误检率最小,仅为10.5%。水体、植被和城区3种地物类别的误检率分别为0、4.8%以及23.1%。这说明本文方法能够有效区分变化较小的地物,具有更好的灵敏性和稳健性。

      表 2  3种方法模拟数据的定量比较

      Table 2.  The Quantitative Comparison of Three Methods for Simulated Data

      方法 类别 TR/% FAR/%
      水体 94.3 5.7
      方法1 植被 84.0 6.0
      建筑物 52.1 47.9
      总体 71.928.1
      水体 99.9 0.1
      方法2 植被 86.9 13.1
      建筑物 70.9 29.1
      总体 84.0 16.0
      水体 100 0
      本文方法 植被 95.2 4.8
      建筑物 76.9 23.1
      总体 89.5 10.5
    • 中、低空间分辨率的遥感影像存在大量的混合像元,传统的变化检测算法往往无法获得良好的检测效果。本文提出了一种结合混合像元分解与EM算法的遥感影像亚像元变化检测方法。首先利用混合像元分解技术,获取像元内部不同组分的百分比后,再进行差值计算,通过EM算法确定变化阈值,完成亚像元级的变化检测。EM阈值选择方法将差异影像的直方图定义为符合高斯分布的概率密度空间,充分考虑了不同端元组成分比例的变化特点,通过迭代优化后确定最优阈值进行识别。实验证明,该算法在正确率和误检率方面具有一定的优势,可以明显地提高检测精度并维持变化地物较好的几何特性,尤其对于一些混合情况比较复杂的区域,检测的效果优于传统变化检测方法。下一步将选择不同分辨率与不同尺度的影像进行实验,同时探索监督型的检测模式,进一步探索分析该方法的适用性。

参考文献 (19)

目录

    /

    返回文章
    返回