基于地基雷达图像的无监督变化检测

黄平平 任慧芳 谭维贤 段盈宏 徐伟 刘方

引用本文:
Citation:

基于地基雷达图像的无监督变化检测

    作者简介:
    黄平平(1978–),男,山东海阳人,博士,教授,2010年获中国科学院电子学研究所博士学位,现任内蒙古工业大学信息工程学院副院长,自治区雷达技术与应用重点实验室主任,“草原英才”创新团队负责人,全国“工人先锋号”负责人。兼任中央军委装备发展部某专家组成员、中国电子学会信号处理分会委员、中国电子教育学会研究生教育分会理事。获自治区优秀共产党员、自治区五一劳动奖章、自治区青年科技奖、自治区自然科学一等奖、自治区科技进步一等奖等荣誉与奖励。入选国家“百千万人才工程”、自治区突出贡献专家、自治区“草原英才”、自治区自然科学杰出青年、高等学校青年科技英才等人才工程。近几年,共主持国家级项目8项,其他各类科研项目20余项,获得授权国家发明专利30余项,发表学术论文130余篇。主要从事结合国家重大需求及内蒙古自治区经济发展需求,新体制雷达系统设计、雷达信号处理和微波遥感应用等方面的研究。E-mail: hpp@imut.edu.cn;
    任慧芳(1992–),女,内蒙古呼和浩特人,现于内蒙古工业大学信息工程学院雷达技术与应用重点实验室攻读硕士学位,主要研究方向为微变监测雷达和遥感图像变化检测。E-mail: Rhf_412922@163.com;
    谭维贤(1981–),男,湖北恩施人,博士,教授,硕士生导师,2009年获中国科学院电子学研究所工学博士学位,2015年至今,内蒙古工业大学雷达技术研究所、内蒙古自治区雷达技术与应用重点实验室任教。入选内蒙古自治区“草原英才”工程,主要从事结合国家重大需求及内蒙古自治区经济发展需求,雷达系统技术、微变监测雷达、雷达信号处理和微波遥感等方面的研究。E-mail: wxtan@imut.edu.cn


    .
    通讯作者: 黄平平 hpp@imut.edu.cn
  • 基金项目:

    国家自然科学基金重点项目(61631011),内蒙古自治区财政厅创新引导项目(KCBJ2017, KCBJ2018014),内蒙古自治区科技重大专项和科技计划项目,装备预研领域基金一般领域基金(JZX7Y20190253040901, JZX7Y20190253041401)

  • 中图分类号: TP753

Unsupervised Change Detection Using Ground-based Radar Image

    Corresponding author: HUANG Pingping, hpp@imut.edu.cn
  • Fund Project: The Key Projects of National Natural Science Foundation of China (61631011), The Innovation Guidance Project of Finance Department of Inner Mongolia Autonomous Region (KCBJ2017, KCBJ2018014), The Major Science and Technology Projects and Science and Technology Programs of Inner Mongolia Autonomous Region, Equipment Pre-research Field Fund General Field Fund (JZX7Y20190253040901, JZX7Y20190253041401)

    CLC number: TP753

  • 摘要: 地基雷达是近20几年逐渐发展成熟的微波遥感成像技术,目前已广泛应用于滑坡、崩塌等地质灾害的监测中。地基雷达通过干涉测量原理可以监测到目标区域发生的微小形变,然而受人为因素、地质因素、气象因素等影响,导致雷达图像失相干严重,给长期定量化监测带来较大的难度。因此,迫切需要在定量监测的基础上,进一步开展变化检测方面的应用,为长期全面了解监测区域的动态变化提供有效信息。针对上述问题,该文提出了一种基于改进的模糊C均值聚类(FCM)算法对地基雷达图像进行无监督变化检测,该方法首次利用相干系数图和均值对数比值图进行非下采样轮廓波变换(NSCT)和局部能量法得到合成差异图,然后利用主成分分析(PCA)提取合成差异图中每个像素的特征向量,根据地基雷达图像特点对FCM进行改进,通过改进的FCM对每个像素的特征向量进行聚类得到最终的变化检测结果。利用地基雷达LSA对中国西南某省出现的堰塞体的治理过程进行监测,获取监测区域的地基雷达图像,监测过程中受降水等影响监测体出现滑坡,使用该文方法对其进行变化检测,结果表明该文方法更容易进行聚类分割,变化检测结果在保留变化区域的同时噪声点明显减少。
  • 图 1  地基雷达图像无监督变化检测流程图

    Figure 1.  Flow chart for unsupervised change detection based on ground-based radar image

    图 2  监测现场图

    Figure 2.  Monitoring site images

    图 3  实验数据

    Figure 3.  Experimental data

    图 4  初始差异图

    Figure 4.  Initial difference images

    图 5  差异图对比

    Figure 5.  Difference images comparison

    图 6  聚类结果对比

    Figure 6.  Clustering result comparison

    图 7  监测现场光学图

    Figure 7.  Optical image of monitoring site

    图 8  变化检测结果与地形3维数据融合图

    Figure 8.  Effect of fusion of change detection results and UAV tilt photography

    表 1  改进的FCM算法计算步骤

    Table 1.  Improved FCM algorithm calculation steps

     算法:改进的FCM算法计算步骤
     输入:融合差异图每个像素的特征向量${{ Y}_\eta }$
     (1) 初始化隶属度$u$,满足式(10);
     (2) 通过式(12)计算聚类中${C_i}$
     (3) 根据式(17)计算目标函数$J$
     (4) 迭代计算:当迭代次数$ \le q$时,
        根据式(18)更新隶属度矩阵,返回步骤(2)继续执行
        当迭代次数$ > q$时,
        停止迭代;
     输出:隶属度矩阵${{ u}_{\eta 1} }$${{ u}_{\eta 2} }$
    下载: 导出CSV

    表 2  地基雷达LSA系统参数

    Table 2.  Ground-based radar LSA system parameters

    参数数值参数数值
    频率范围16.5~17.5 GHz监测范围(方位向)60°× 30°
    形变监测精度0.1 mm脉冲重复频率优于500 Hz
    作用距离5 km位移观测周期240~600 s
    图像分辨率0.3 m × 0.0054 rad供电电源AC 220 V/DC 24 V
    设备总重$ \le 150\;{\rm{kg} }$工作温度–30~ 50 ℃
    下载: 导出CSV
  • [1] WANG Yan, DU Lan, and DAI Hui. Unsupervised SAR image change detection based on SIFT keypoints and region information[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(7): 931–935. doi: 10.1109/LGRS.2016.2554606
    [2] GONG Maoguo, SU Linzhi, JIA Meng, et al. Fuzzy clustering with a modified MRF energy function for change detection in synthetic aperture radar images[J]. IEEE Transactions on Fuzzy Systems, 2014, 22(1): 98–109. doi: 10.1109/TFUZZ.2013.2249072
    [3] GAO Feng, DONG Junyu, LI Bo, et al. Automatic change detection in synthetic aperture radar images based on PCANet[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(12): 1792–1796. doi: 10.1109/LGRS.2016.2611001
    [4] GONG Maoguo, CAO Yu, and WU Qiaodi. A neighborhood-based ratio approach for change detection in SAR images[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(2): 307–311. doi: 10.1109/LGRS.2011.2167211
    [5] GONG Maoguo, ZHOU Zhiqiang, and MA Jingjing. Change detection in synthetic aperture radar images based on image fusion and fuzzy clustering[J]. IEEE Transactions on Image Processing, 2012, 21(4): 2141–2151. doi: 10.1109/TIP.2011.2170702
    [6] SUMAIYA M N and KUMARI R S S. Logarithmic mean-based thresholding for SAR image change detection[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(11): 1726–1728. doi: 10.1109/LGRS.2016.2606119
    [7] 杨祥立, 徐德伟, 黄平平, 等. 融合相干/非相干信息的高分辨率SAR图像变化检测[J]. 雷达学报, 2015, 4(5): 582–590. doi: 10.12000/JR15073YANG Xiangli, XU Dewei, HUANG Pingping, et al. Change detection of high resolution SAR images by the fusion of coherent/incoherent information[J]. Journal of Radars, 2015, 4(5): 582–590. doi: 10.12000/JR15073
    [8] ELGUEBALY T and BOUGUILA N. A Bayesian approach for SAR images segmentation and changes detection[C]. The 25th Biennial Symposium on Communications, Kingston, Canada, 2010: 24–27. doi: 10.1109/BSC.2010.5473011.
    [9] 邵宁远, 邹焕新, 陈诚, 等. 面向变化检测的SAR图像超像素协同分割算法[J]. 系统工程与电子技术, 2019, 41(7): 1496–1503. doi: 10.3969/j.issn.1001-506X.2019.07.09SHAO Ningyuan, ZOU Huanxin, CHEN Cheng, et al. Change detection oriented superpixel cosegmentation algorithm for SAR images[J]. Systems Engineering and Electronics, 2019, 41(7): 1496–1503. doi: 10.3969/j.issn.1001-506X.2019.07.09
    [10] CARINCOTTE C, DERRODE S, and BOURENNANE S. Unsupervised change detection on SAR images using fuzzy hidden Markov chains[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(2): 432–441. doi: 10.1109/TGRS.2005.861007
    [11] BEZDEK J C, EHRLICH R, and FULL W. FCM: The fuzzy c-means clustering algorithm[J]. Computers & Geosciences, 1984, 10(2/3): 191–203. doi: 10.1016/0098-3004(84)90020-7
    [12] PHAM D L. Spatial models for fuzzy clustering[J]. Computer Vision and Image Understanding, 2001, 84(2): 285–297. doi: 10.1006/cviu.2001.0951
    [13] TRAN T N, WEHRENS R, and BUYDENS L M C. Clustering multispectral images: A tutorial[J]. Chemometrics and Intelligent Laboratory Systems, 2005, 77(1/2): 3–17. doi: 10.1016/j.chemolab.2004.07.011
    [14] LI Zhenfang, BAO Zheng, LI Hai, et al. Image autocoregistration and InSAR interferogram estimation using joint subspace projection[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(2): 288–297. doi: 10.1109/TGRS.2005.860984
    [15] HYDE P, DUBAYAH R, WALKER W, et al. Mapping forest structure for wildlife habitat analysis using multi-sensor (LiDAR, SAR/InSAR, ETM+, Quickbird) synergy[J]. Remote Sensing of Environment, 2006, 102(1/2): 63–73. doi: 10.1016/j.rse.2006.01.021
    [16] ADU J, GAN Jianhong, WANG Yan, et al. Image fusion based on nonsubsampled contourlet transform for infrared and visible light image[J]. Infrared Physics & Technology, 2013, 61: 94–100. doi: 10.1016/j.infrared.2013.07.010
    [17] WU Chao and WU Yiquan. Multitemporal images change detection using nonsubsampled contourlet transform and kernel fuzzy C-means clustering[C]. The 2nd International Symposium on Intelligence Information Processing and Trusted Computing, Wuhan, China, 2011: 22–23. doi: 10.1109/IPTC.2011.31.
    [18] 李青松, 覃锡忠, 贾振红, 等. 非下采样Contourlet域融合和参数化内核图割的SAR图像无监督水灾变化检测[J]. 中国图象图形学报, 2014, 19(6): 971–978. doi: 10.11834/jig.20140619LI Qingsong, QIN Xizhong, JIA Zhenhong, et al. Unsupervised detection of flood changes with SAR images combining nonsubsampled Contourlet domain fusion and parametric kernel graph cuts[J]. Journal of Image and Graphics, 2014, 19(6): 971–978. doi: 10.11834/jig.20140619
    [19] 刘陆洋, 贾振红, 杨杰, 等. 利用双差异图和PCA的SAR图像变化检测[J]. 计算机工程与设计, 2019, 40(7): 2002–2006.LIU Luyang, JIA Zhenhong, YANG Jie, et al. SAR image chang detective using double difference images and PCA algorithm[J]. Computer Engineering and Design, 2019, 40(7): 2002–2006.
    [20] 慕彩红, 霍利利, 刘逸, 等. 基于小波融合和PCA-核模糊聚类的遥感图像变化检测[J]. 电子学报, 2015, 43(7): 1375–1381. doi: 10.3969/j.issn.0372-2112.2015.07.019MU Caihong, HUO Lili, LIU Yi, et al. Change detection for remote sensing images based on wavelet fusion and PCA-Kernel fuzzy clustering[J]. Acta Electronica Sinica, 2015, 43(7): 1375–1381. doi: 10.3969/j.issn.0372-2112.2015.07.019
    [21] EL-DAWLATLY S, OSMAN H, and SHAHEIN H I. New spatial FCM approach with application to SAR target clustering[C]. The 8th International Conference on Signal Processing, Beijing, China, 2006. doi: 10.1109/ICOSP.2006.345896.
  • [1] 杨祥立徐德伟黄平平杨文 . 融合相干/非相干信息的高分辨率SAR图像变化检测. 雷达学报, doi: 10.12000/JR15073
    [2] 郑瑾尤红建 . 基于Radon 变换和Jeffrey 散度的SAR 图像变化检测方法. 雷达学报, doi: 10.3724/SP.J.1300.2012.10068
    [3] 浮瑶瑶柳彬张增辉郁文贤 . 基于词包模型的高分辨率SAR 图像变化检测与分析. 雷达学报, doi: 10.3724/SP.J.1300.2014.13134
    [4] 徐真王宇李宁张衡张磊 . 一种基于CNN的SAR图像变化检测方法. 雷达学报, doi: 10.12000/JR17075
    [5] 冷英李宁 . 一种改进的变化检测方法及其在洪水监测中的应用. 雷达学报, doi: 10.12000/JR16139
    [6] 冀广宇董勇伟卜运成李焱磊周良将梁兴东 . 基于目标相干性表征差异的多波段SAR相干变化检测方法. 雷达学报, doi: 10.12000/JR18020
    [7] 赵军香梁兴东李焱磊 . 一种基于似然比统计量的SAR相干变化检测. 雷达学报, doi: 10.12000/JR16065
    [8] 赵耀东吕晓德李纪传向茂生 . 无源雷达多普勒谱分析实现动目标检测的方法. 雷达学报, doi: 10.3724/SP.J.1300.2012.20081
    [9] 赵永科吕晓德 . 一种联合优化的无源雷达实时目标检测算法. 雷达学报, doi: 10.12000/JR14005
    [10] 刘振魏玺章黎湘 . 一种新的随机PRI脉冲多普勒雷达无模糊MTD算法. 雷达学报, doi: 10.3724/SP.J.1300.2013.10063
    [11] 赵勇胜赵拥军赵闯 . 联合角度和时差的单站无源相干定位加权最小二乘算法. 雷达学报, doi: 10.12000/JR15133
    [12] 滑文强王爽侯彪 . 基于半监督学习的SVM-Wishart极化SAR图像分类方法. 雷达学报, doi: 10.12000/JR14138
    [13] 邢艳肖张毅李宁王宇胡桂香 . 一种联合特征值信息的全极化SAR图像监督分类方法. 雷达学报, doi: 10.12000/JR16019
    [14] 钟能杨文杨祥立郭威 . 基于混合Wishart模型的极化SAR图像非监督分类. 雷达学报, doi: 10.12000/JR16133
    [15] 邹焕新罗天成张月周石琳 . 基于组合条件随机场的极化SAR图像监督地物分类. 雷达学报, doi: 10.12000/JR16109
    [16] 滑文强王爽郭岩河谢雯 . 基于邻域最小生成树的半监督极化SAR图像分类方法. 雷达学报, doi: 10.12000/JR18104
    [17] 尹嫱李洋黄平平林赟洪文 . 土壤湿度变化引起的干涉SAR相干性损失分析(英文). 雷达学报, doi: 10.12000/JR15075
    [18] 张向荣于心源唐旭侯彪焦李成 . 基于马尔科夫判别谱聚类的极化SAR图像分类方法. 雷达学报, doi: 10.12000/JR19059
    [19] 胡程邓云开田卫明曾涛 . 地基干涉合成孔径雷达图像非线性大气相位补偿方法. 雷达学报, doi: 10.12000/JR19073
    [20] 曾丽娜周德云李枭扬张堃 . 基于无训练单样本有效特征的SAR目标检测. 雷达学报, doi: 10.12000/JR16114
  • 加载中
图(8)表(2)
计量
  • 文章访问数:  1068
  • HTML浏览量:  208
  • PDF下载量:  40
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-01-11
  • 录用日期:  2020-03-22
  • 网络出版日期:  2020-04-08

基于地基雷达图像的无监督变化检测

    通讯作者: 黄平平 hpp@imut.edu.cn
    作者简介:
    黄平平(1978–),男,山东海阳人,博士,教授,2010年获中国科学院电子学研究所博士学位,现任内蒙古工业大学信息工程学院副院长,自治区雷达技术与应用重点实验室主任,“草原英才”创新团队负责人,全国“工人先锋号”负责人。兼任中央军委装备发展部某专家组成员、中国电子学会信号处理分会委员、中国电子教育学会研究生教育分会理事。获自治区优秀共产党员、自治区五一劳动奖章、自治区青年科技奖、自治区自然科学一等奖、自治区科技进步一等奖等荣誉与奖励。入选国家“百千万人才工程”、自治区突出贡献专家、自治区“草原英才”、自治区自然科学杰出青年、高等学校青年科技英才等人才工程。近几年,共主持国家级项目8项,其他各类科研项目20余项,获得授权国家发明专利30余项,发表学术论文130余篇。主要从事结合国家重大需求及内蒙古自治区经济发展需求,新体制雷达系统设计、雷达信号处理和微波遥感应用等方面的研究。E-mail: hpp@imut.edu.cn;
    任慧芳(1992–),女,内蒙古呼和浩特人,现于内蒙古工业大学信息工程学院雷达技术与应用重点实验室攻读硕士学位,主要研究方向为微变监测雷达和遥感图像变化检测。E-mail: Rhf_412922@163.com;
    谭维贤(1981–),男,湖北恩施人,博士,教授,硕士生导师,2009年获中国科学院电子学研究所工学博士学位,2015年至今,内蒙古工业大学雷达技术研究所、内蒙古自治区雷达技术与应用重点实验室任教。入选内蒙古自治区“草原英才”工程,主要从事结合国家重大需求及内蒙古自治区经济发展需求,雷达系统技术、微变监测雷达、雷达信号处理和微波遥感等方面的研究。E-mail: wxtan@imut.edu.cn


  • 内蒙古工业大学信息工程学院 呼和浩特 010051
  • 内蒙古自治区雷达技术与应用重点实验室 呼和浩特 010051
基金项目:  国家自然科学基金重点项目(61631011),内蒙古自治区财政厅创新引导项目(KCBJ2017, KCBJ2018014),内蒙古自治区科技重大专项和科技计划项目,装备预研领域基金一般领域基金(JZX7Y20190253040901, JZX7Y20190253041401)

摘要: 地基雷达是近20几年逐渐发展成熟的微波遥感成像技术,目前已广泛应用于滑坡、崩塌等地质灾害的监测中。地基雷达通过干涉测量原理可以监测到目标区域发生的微小形变,然而受人为因素、地质因素、气象因素等影响,导致雷达图像失相干严重,给长期定量化监测带来较大的难度。因此,迫切需要在定量监测的基础上,进一步开展变化检测方面的应用,为长期全面了解监测区域的动态变化提供有效信息。针对上述问题,该文提出了一种基于改进的模糊C均值聚类(FCM)算法对地基雷达图像进行无监督变化检测,该方法首次利用相干系数图和均值对数比值图进行非下采样轮廓波变换(NSCT)和局部能量法得到合成差异图,然后利用主成分分析(PCA)提取合成差异图中每个像素的特征向量,根据地基雷达图像特点对FCM进行改进,通过改进的FCM对每个像素的特征向量进行聚类得到最终的变化检测结果。利用地基雷达LSA对中国西南某省出现的堰塞体的治理过程进行监测,获取监测区域的地基雷达图像,监测过程中受降水等影响监测体出现滑坡,使用该文方法对其进行变化检测,结果表明该文方法更容易进行聚类分割,变化检测结果在保留变化区域的同时噪声点明显减少。

English Abstract

    • 我国地质构造复杂、地形地貌起伏变化大,极易发生滑坡、崩塌、泥石流等地质灾害。此外,人类的生产、生活活动也很容易诱发滑坡灾害的发生,大量岩土边坡的失稳或滑动屡见不鲜,给工程建设和人民生命财产带来了严重的危害和损失。灾害发生区域的变化检测分析对安全生产有着重要意义。从遥感数据中获取灾害信息是一种十分重要的灾害研究手段,而快速有效的变化检测方法对于掌握灾害情况尤为重要。地基雷达技术具有全天时、全天候、大范围、远距离、非接触等特点,在地质灾害监测、预警方面具有巨大优势。目前,地基雷达已广泛应用于滑坡、崩塌等地质灾害的监测,利用干涉测量原理地基雷达可以监测到目标区域发生的微小形变量,根据累计形变量、形变速度、形变加速度等对地质灾害进行实时监测。然而,在长期施工建设监测中,部分区域的变化并不会引起滑坡等灾害,但受人为因素、地质因素、气象因素等影响,导致监测区域的雷达图像失相干较为严重,给长期定量化监测带来了较大的难度。因此,迫切需要在定量监测的基础上,进一步开展变化检测方面的应用,为长期全面了解监测区域的动态变化提供有效信息。针对上述问题,本文将研究一种基于地基雷达图像的变化检测方法。

      一般来说,图像变化检测技术可以分为两大类:监督式和无监督式。监督变化检测方法,需要先验信息;而无监督变化检测方法直接自动比较多时相的遥感图像,不需要任何额外信息,降低人为误差的影响,符合实际应用中先验变化信息缺失的现实情况,因此得到了广泛地研究[1]。雷达图像中的无监督变化检测一般包括3个步骤:(1)图像预处理;(2)获取差分图像;(3)分析差异图像及后处理[2,3]

      近年来雷达图像变化检测的研究主要集中在机载和星载合成孔径雷达(Synthetic Aperture Radar, SAR)图像上,基于星载SAR图像中的无监督变化检测,龚茂国等人[4]提出了一种基于邻域的比值法来产生差异图,虽然该算法优于产生差异图的传统方法,如差值法和比值法,但是其检测结果受噪声干扰严重;为了提高变化检测结果的准确率,使用均值比图像和对数比图像的互补信息来生成差异图,取得了较好的检测结果[5]。Sumaiya等人[6]采用了一种经济且简单的对数均值阈值进行星载SAR图像变化检测。杨祥立等人[7]运用D-S证据理论(Dempster-Shafer evidence theory, D-S)融合高分辨率SAR影像的相干/非相干差异特征进行变化检测,采用美国旧金山地区港口附近的TerraSAR-X数据切片进行实验验证。此外,在图像变化检测中一个重要的步骤是图像的分割聚类,Elguebaly等人[8]提出了一种基于广义高斯混合模型(Generalized Gaussian Mixture models, GGM)的高效无监督图像分割和变化检测算法,论证了贝叶斯GGM在图像分割中的应用潜力。邵宁远等人[9]提出了一种面向变化检测的SAR图像超像素协同分割算法,解决了变化检测中存在的多时相图像边界和空间对应关系不一致的问题。Carincotte等人[10]提出了一种新的模糊隐马尔可夫链(Hidden Markov ChainS, HMCS),从而解决了模糊变化检测的统计方法。而最著名的模糊聚类算法是模糊C均值聚类(Fuzzy C-Means, FCM)分割算法[11],它主要的缺点是在成像应用中缺乏对空间环境信息的整合,导致其对噪声和其它人为因素较为敏感。而FCM的变形,如鲁棒FCM(Robust FCM, RFCM)[12]和空间条件FCM(the Spatial Conditional FCM, SCFCM)[13]在星载SAR图像聚类中都得到了较好的效果。

      星载SAR适用于大范围对地成像,其成像范围可达数千平方公里,基于星载SAR图像变化检测适用于监测区域内发生的大面积、特征明显的变化。地基雷达通常针对特定的目标区域成像,多数情况下发生变化的特征不够明显,将星载SAR图像变化检测技术应用到地基雷达图像变化检测中,会造成差异图群分性比较差和聚类分割结果存在较多噪声点等缺陷。本文提出了一种无监督地基雷达图像变化检测方法,本方法通过对相干系数图和均值对数比值图进行非下采样轮廓波变换(NonsubSampled Contourlet Transform, NSCT)和局部能量法得到合成差异图,得到的合成差异图具有较好的群分性。然后采用主成分分析法(Principal Component Analysis, PCA)提取特征向量,增强抗噪性。最后使用改进的模糊C均值聚类(FCM)对特征向量进行聚类得到最终的变化检测结果。文章安排如下:本文第2部分提出了一种地基雷达图像变化检测方法;第3部分利用地基雷达实测数据进行结果分析并验证该方法的有效性;第4部分对全文做出总结。

    • 本文所提基于地基雷达图像的无监督变化检测方法主要分为4个步骤:

      (1) 初始差异图获取。先取两幅已配准的地基雷达图像${X_1}$${X_2}$,通过相干系数算子以及均值对数比算子生成相干系数差异图和均值对数比差异图;

      (2) NSCT变换。对得到的两幅初始差异图,进行NSCT变换并融合低频系数,然后使用融合的低频系数分别与初始差异图的高频系数进行NSCT反变换得到两幅融合结果图,对这两幅融合结果图进行像素能量合成得到最终的合成差异图;

      (3) PCA提取特征向量。利用PCA算法提取出合成差异图的像素特征;

      (4) 改进的FCM聚类。最后通过一种加相似性惩法项的模糊C-均值聚类算法(FCM)聚类,获取变化检测结果。

      本文的算法流程如图1所示。

      图  1  地基雷达图像无监督变化检测流程图

      Figure 1.  Flow chart for unsupervised change detection based on ground-based radar image

    • 由于地基雷达图像自身存在乘性斑点噪声,对数比值法能够将乘性斑点噪声转化为加性噪声,便于有效去除斑点噪声。此外,对数运算作用于两幅地基雷达图中能够对像素的灰度值的比值做非线性拉伸,可以在一定程度上增强变化和非变化区域的对比。为了充分考虑像素点周围的特征,采用均值对数比获取的差异图为${X_{\rm{d}}}$,计算公式为

      $ {X_{\rm{d}}}(i,j) = {\rm{log}}\left(\frac{{{m_{\rm{1}}}(i,j) + 1}}{{{m_{\rm{2}}}(i,j) + 1}}\right) $

      其中,${m_1}(i,j)$${m_2}(i,j)$分别为图像${X_1}$${X_2}$在像素点$(i,j)$处的$3 \times 3$邻域均值。

      雷达图像之间的相干系数用来衡量影像之间的相似程度,它提供了检测区域重要的变化信息。采用相干系数图,可以得到变化区域非常明显和背景信息较为平滑的差异图[14,15]。图像${X_1}$${X_2}$的相干系数定义如式(2)所示

      $ \begin{split} & {X_{\rm{m}}}(i,j) = \\ & \frac{{\displaystyle\sum {\left| {{X_1}(i,j) \times X_2^*(i,j)} \right|} }}{{\sqrt {\displaystyle\sum {\left| {{X_1}(i,j) \times X_1^*(i,j)} \right|} \times \displaystyle\sum {\left| {{X_2}(i,j) \times X_2^*(i,j)} \right|} } }} \end{split} $

      式中,$*$表示共轭复数,${X_1}$${X_2}$为经过配准滤波的地基雷达图像,相干系数${X_{\rm{m}}}(i,j)$的值在$0{\text{~}} 1$之间,0表示失相干,1表示完全相干。

      采用对数比值图和相干系数图相结合的方法,综合利用各自的优势,能够得到群分性较好的融合差异图,为后续的图像聚类分割提供重要依据。

    • NSCT是基于非下采样金字塔滤波器组(Non Subsamlped Pyramid Filter Banks, NSPFB)和非下采样定向滤波器组(Non Subsampled Directional Filter Banks, NSDFB)进行的,每个子带图像与原始图像在形状和大小上都是完全一样的,因此,NSCT具有多分辨率、局部化、方向性、各向异性和位移不变性等特点[16]。该方法能较好地捕捉图像的几何细节,保持目标信息和边缘信息[17]

      基于NSCT的图像融合算法通常将图像进行NSCT分解,针对不同频率域的特点选择不同的融合规则。对数比值图${X_{\rm{d}}}$和相干系数图${X_{\rm{m}}}$的融合结果图为${F_{\rm{d}}}$,具体实现步骤如下:

      (1) NSCT变换。轮廓波变换使用“9-7”型拉普拉斯金字塔滤波器,使用“pkva”作为方向组滤波器,每个金字塔级方向滤波器组分解层数为4。经NSCT变换得到差异图${X_{\rm{d}}}$${X_{\rm{m}}}$对应的低频系数${L_{\rm{d}}}$${L_{\rm{m}}}$,高频系数${H_{\rm{d}}}$${H_{\rm{m}}}$,高、低频系数反映了高、低频图像信息。

      (2) 加权平均法。由于图像的低频信息通常反映图像的结构信息[18],因此在融合过程中,为了凸显变化区域,针对低频系数采取加权平均的方式,对差异图${X_{\rm{d}}}$${X_{\rm{m}}}$经NSCT变换得到的低频系数${L_{\rm{d}}}$${L_{\rm{m}}}$进行融合,得到融合后的低频系数${L_{\rm{f}}}$,计算过程为[19]

      ${L_{\rm{f}}} = \alpha {L_{\rm{m}}} + (1 - \alpha ){L_{\rm{d}}}$

      其中,$\alpha (\alpha \ge 0)$为加权平均的融合系数,本文实验中$\alpha = 0.3$

      (3) NSCT反变换。用融合的低频系数${L_{\rm{f}}}$分别与差异图${X_{\rm{d}}}$${X_{\rm{m}}}$的高频系数进行NSCT反变换,得到融合结果图${F_{\rm{d}}}$${F_{\rm{m}}}$[19]

      因为地基雷达图像的相干系数图的背景相对平坦,可以采用局部能量法抑制背景信息[5],得到最终的合成差异图${F_{\rm{h}}}$

      $ {F_{\rm{h}}}(i,j) = \left\{ \begin{aligned} & {F_{\rm{m}}}(i,j),\ \ {{E_{\rm{m}}}(i,j) \le } {E_{\rm{d}}}(i,j)\\ &{F_{\rm{d}}}(i,j),\ \ {\text{其他}} \end{aligned} \right. $

      $ {E_{({\rm{m,d)}}}}(i,j) = \sum\limits_{q \in {L_{i,j}}} {{{\left[ {{F_{{\rm{(m,d)}}}}(q)} \right]}^2}} \hspace{40pt} $

      其中,${E_{\rm{d}}}(i,j)$${E_{\rm{m}}}(i,j)$分别表示差异图${F_{\rm{d}}}$${F_{\rm{m}}}$的局部能量,${L_{i,j}}$表示以像素$(i,j)$为中心,大小为$3 \times 3$的邻域,${F_{\rm{d}}}(q)$${F_{\rm{m}}}(q)$表示邻域内第$q$个像素值[10]

    • 主成分分析(PCA),也称Karhunen-Leve变换(简称K-L变换),是一种线性变换。它是一种均方误差最小的最佳正交变换,是在统计特征基础上的多维线性变换。本文采用PCA方法对融合差异图进行特征提取,由于采用分块的思想,使得该方法具有一定的抗噪能力[20]

      将合成差异图分为$n$个大小为$3 \times 3$的相邻的小块,小块之间互不重叠。每个小块可以看作一个矩阵,然后将每个小块矩阵转换为列向量${{{P}}^t}$,其中$t = 1,2, ··· ,n$,并计算全部列向量的均值向量

      ${{m}} = \frac{1}{n}\sum\limits_{t = 1}^n {{{{P}}^t}} $

      将所有列向量重构为${3^2} \times n$的矩阵${{S}}$,计算${{S}}$的协方差矩阵${{C}}$

      ${{C}} = \frac{1}{n}\sum\limits_{t = 1}^n {({{{P}}^t} - {{m}}){{({{{P}}^t} - {{m}})}^{\rm{T}}}} $

      其中,${\rm{T}}$表示矩阵转置。协方差矩阵${{C}}$是一个大小为${3^2} \times {3^2}$的矩阵。对协方差矩阵进行特征值分解,求出特征值和特征向量,按特征值从大到小的顺序排列,选出对应的特征向量,这些特征向量形成一组正交基,取特征向量的前$r$个列向量构成矩阵${{A}}$, ${{A}}$是一个大小为${3^2} \times r$的矩阵。取合成差异图${F_{\rm{h}}}$中每个像素所在的$3 \times 3$邻域小块,将每个小块转化为列向量${{{Y}}_\eta }$, $\eta = 1,2, ··· ,k$, $k$为合成差异图${F_{\rm{h}}}$中像素的总个数。将所有的列向量映射到矩阵${{A}}$上,即

      ${{{Y}}_\eta } = {{{A}}^{\rm{T}}}({{{V}}_\eta } - {{m}})$

      其中,${{{Y}}_\eta }$为每个像素代表的特征向量,${{{A}}^{\rm T}}$${{A}}$的转置矩阵,这样每个像素都用一个$R$维的特征向量来表示,$\left[ {{Y_1}} \ \ {{Y_2}} \ \ ··· \ \ {{Y_k}} \right]$构成特征向量空间${{Q}}$${{Q}}$即为一个$r \times k$矩阵。

    • FCM是一种无监督的模糊聚类方法,传统FCM[11]目标函数的形式为

      $ {J_{{\rm{FCM}}}} = \sum\limits_{\eta = 1}^k {\sum\limits_{i = 1}^c {{{({{ u}_{\eta i}})}^m}{{({{{Y}}_\eta } - {C_i})}^2}} } $

      其中,${{ u}_{\eta i}}$是像素点${{{Y}}_\eta }$$i$类的隶属度矩阵,$c$为分类数,$m$为模糊权重(一般为2), $\eta $为图像的像素点个数,${C_i}$为聚类中心。

      ${{ u}_{\eta i}}$表示${{{Y}}_\eta }$隶属于第$i$类的概率,介于0和1之间。因此,约束必须满足

      $ \sum\limits_{i = 1}^c {{{ u}_{\eta i}} = 1} $

      考虑到式(10)中的约束,通过拉格朗日数乘算子计算当式(9)中的目标函数达到最小值时,${{ u}_{\eta i}}$

      $ {{ u}_{\eta i}} = \frac{1}{{\displaystyle\sum\limits_{s = 1}^c {{{\left( {\dfrac{{{{\left\| {{{{Y}}_\eta } - {C_i}} \right\|}^2}}}{{{{\left\| {{{{Y}}_\eta } - {C_s}} \right\|}^2}}}} \right)}^{1/{(m-1)}}}} }} $

      聚类中心${C_i}$

      $ {C_i} = \frac{{\displaystyle\sum\limits_{\eta = 1}^k {{ u}_{\eta i}^m{{{Y}}_\eta }} }}{{\displaystyle\sum\limits_{\eta = 1}^k {{ u}_{\eta i}^m} }} $

      当用于图像聚类时,FCM没有考虑到相邻图像像素之间的空间关系[21]。采用变形的FCM为鲁棒模糊C-均值(RFCM)[12]的目标函数采取了如式(13)所示形式

      $ \begin{split} {J_{{\rm{RFCM}}}} =\,& \sum\limits_{\eta = 1}^k \sum\limits_{i = 1}^c { u}_{\eta k}^m{{\left\| {{{{Y}}_\eta } - {C_i}} \right\|}^2} \\ & + \frac{\beta }{2}\sum\limits_{\eta = 1}^k {\sum\limits_{i = 1}^c {{ u}_{\eta k}^m\sum\limits_{l \in {{{N}}_\eta }} {\sum\limits_{p \in {{{P}}_i}} {{ u}_{lp}^m} } } } \end{split} $

      其中,${{{N}}_\eta }$是像素点${{{Y}}_\eta }$的空间邻域集合,${{{P}}_i}$是除$i$类外的所有类的集合,$\beta $是控制附加空间惩罚项效果的正偏差。

      增加的空间惩罚项鼓励像素点在同一类中拥有有较高的聚类资格。最小化式(13)得到隶属度函数为

      $ {{ u}_{\eta i}} \!=\! \frac{{{{\left( {{{\left\| {{{{Y}}_\eta } - {C_i}} \right\|}^2} + \beta \displaystyle\sum\limits_{l \in {{{N}}_\eta }} {\displaystyle\sum\limits_{p \in {{{P}}_i}} {{ u}_{lp}^m} } } \right)}^{1/{(m-1)}}}}}{{\displaystyle\sum\limits_{s = 1}^c {{{\left(\!\! {{{\left\| {{{{Y}}_\eta } \!-\! {C_s}} \right\|}^2} \!+\! \beta\! \displaystyle\sum\limits_{l \in {{{N}}_\eta }} {\displaystyle\sum\limits_{p \in {{{P}}_i}}\! {{ u}_{lp}^m} } } \!\right)}^{1/{(m-1)}}}} }} $

      ${{{C}}_i}$的计算过程为式(12),因为惩罚项不依赖于它。

      FCM的另一个变型也在其目标函数中包含了空间关系,它是在目标函数中引入的空间条件FCM(也称SCFCM)[13]。具体来说,它的目标函数采用如式(15)的形式

      $ {J_{{\rm{SCFCM}}}} = \sum\limits_{\eta = 1}^k {\displaystyle\sum\limits_{i = 1}^c {{ u}_{\eta k}^m{{\left\| {{{{Y}}_\eta } - {C_i}} \right\|}^2}\frac{{\displaystyle\sum\limits_{l \in {{{N}}_\eta }} {{{\rm e}^{{\rm{ - }}\beta {{{u}}_{li}}}}} }}{{\displaystyle\sum\limits_{l \in {{{N}}_\eta }} {\displaystyle\sum\limits_{r = 1}^c {{{\rm e}^{{\rm{ - }}\beta {{{u}}_{lr}}}}} } }}} } $

      最小化式(15)得到隶属度函数为

      $ {{ u}_{\eta i}} = \frac{{{{\left(\! {{{\left\| {{{{Y}}_\eta } - {C_i}} \right\|}^2}\frac{{\displaystyle\sum\limits_{l \in {{{N}}_\eta }} {{{\rm e}^{{\rm{ - }}\beta {{{u}}_{li}}}}} }}{{\displaystyle\sum\limits_{l \in {{{N}}_\eta }} {\displaystyle\sum\limits_{r = 1}^c {{{\rm e}^{{\rm{ - }}\beta {{{u}}_{lr}}}}} } }}} \!\right)}^{1/{(m-1)}}}}}{{\displaystyle\sum\limits_{s = 1}^c {{{\left(\!\! {{{\left\| {{{{Y}}_\eta } \!-\! {C_s}} \right\|}^2}\frac{{\displaystyle\sum\limits_{l \in {{{N}}_\eta }} {{{\rm e}^{{\rm{ - }}\beta {{{u}}_{ls}}}}} }}{{\displaystyle\sum\limits_{l \in {{{N}}_\eta }} {\displaystyle\sum\limits_{r = 1}^c {{{\rm e}^{{\rm{ - }}\beta {{{u}}_{lr}}}}} } }}} \!\!\right)}^{1/{(m-1)}}}} }} $

      ${{{C}}_i}$的计算过程为式(12)。

    • 为了弥补FCM和其变形RFCM和SCFCM分割率低以及对于地基雷达图像噪声敏感等问题。本文对FCM的目标函数进行改进,在原始的FCM的目标函数中加入相似性惩罚项。目的是使图像中相类似的像素点以更大概率聚类到相同的类,提高聚类效果并抑制噪声点。改进后的聚类算法的目标函数的表达式为

      $ \begin{split} J =\,& \displaystyle\sum\limits_{\eta = 1}^k \sum\limits_{i = 1}^c {{({{ u}_{\eta i}})}^m}{{({{{Y}}_\eta } - {C_i})}^2} - \beta {{({{ u}_{\eta i}})}^m}\\ & \cdot\ln \left[\sum\limits_{r = 1}^c \sum\limits_{l \in {{{N}}_\eta }} 1 - {{({{ u}_{lr}})}^m}\right] \end{split} $

      其中,${{{Y}}_\eta }$为输入每个像素的特征向量,分类数$c = 2$(变化类和非变化类),$m$模糊权重(一般为2), $\beta $控制惩罚项的积极因子,本文中$\beta = 0.2615$${{{N}}_\eta }$${{{Y}}_\eta }$及其周围邻域的一个集合,$l$为此集合中的一个元素。

      改进的目标函数在保留传统FCM目标函数的基础上,添加了相似性惩罚项。需要强调的是,在式(17)的邻域中,邻域是以相似性而不是以空间定义的,这样可以找到每个输入图像的最相似邻域,以鼓励相似的像素以更大的概率聚类到同一类中。通过拉格朗日数乘算子计算当式(17)中的目标函数达到最小值时,${{ u}_{\eta i}}$

      $ \begin{split} & {{ u}_{\eta i}} = \\ & \frac{1}{{\displaystyle\sum\limits_{s = 1}^c {\frac{{{{\left( {{{\left\| {{{{Y}}_\eta } - {C_i}} \right\|}^2} - \beta \ln \left[\displaystyle\sum\limits_{r = 1}^c \displaystyle\sum\limits_{l \in {{{N}}_\eta }} 1 - {{({{ u}_{lr}})}^m}\right] } \right)}^{\frac{1}{{m - 1}}}}}}{{{{\left( {{{\left\| {{{{Y}}_\eta } - {C_s}} \right\|}^2} - \beta \ln \left[\displaystyle\sum\limits_{r = 1}^c \displaystyle\sum\limits_{l \in {{{N}}_\eta }} 1 - {{({{ u}_{lr}})}^m}\right] } \right)}^{\frac{1}{{m - 1}}}}}}} }} \end{split} $

      由于相似度惩罚项不包含${C_i}$,聚类中心采用式(12)中给出的FCM算法的相同形式。改进的FCM聚类算法计算过程如表1所示。

       算法:改进的FCM算法计算步骤
       输入:融合差异图每个像素的特征向量${{ Y}_\eta }$
       (1) 初始化隶属度$u$,满足式(10);
       (2) 通过式(12)计算聚类中${C_i}$;
       (3) 根据式(17)计算目标函数$J$;
       (4) 迭代计算:当迭代次数$ \le q$时,
          根据式(18)更新隶属度矩阵,返回步骤(2)继续执行
          当迭代次数$ > q$时,
          停止迭代;
       输出:隶属度矩阵${{ u}_{\eta 1} }$和${{ u}_{\eta 2} }$

      表 1  改进的FCM算法计算步骤

      Table 1.  Improved FCM algorithm calculation steps

      通过比较每个像素的隶属度${{ u}_{\eta i}}$的大小,${{ u}_{\eta 1}}$${{ u}_{\eta 2}}$分别代表第$\eta $个像素对应两类隶属度,采用如式(19)所示的判决条件对${{ u}_{\eta 1}}$${{ u}_{\eta 2}}$进行二值化处理

      $ R(i,j) = \left\{ \begin{aligned} & 1,\ \ {{{ u}_{\eta 1}} \le {{ u}_{\eta 2}}}\\ & 0,\ \ {\text{其他}} \end{aligned} \right. $

      $R$为最终变化检测结果图。

    • 我国西南某省发生地震之后山体严重垮塌,形成堰塞体,为确保堰塞体整治工程安全、顺利进行,在现场部署地基雷达LSA对堰塞体治理施工过程进行监测。经现场调研,地基雷达LSA布置在距堰塞体300 m正对面,雷达的监测范围完全可以有效覆盖堰塞体。表2为地基雷达LSA系统参数,图2所示为监测现场光学图。监测期间,施工地出现持续性强降水,导致边坡发生滑坡。如图3(a)为滑坡前雷达图,数据获取时间为2019年7月15日,图3(b)为滑坡后雷达图,数据获取时间为2019年7月19日。截取滑坡前后数据切片对进行实验,如图3(c)图3(d)所示。

      参数数值参数数值
      频率范围16.5~17.5 GHz监测范围(方位向)60°× 30°
      形变监测精度0.1 mm脉冲重复频率优于500 Hz
      作用距离5 km位移观测周期240~600 s
      图像分辨率0.3 m × 0.0054 rad供电电源AC 220 V/DC 24 V
      设备总重$ \le 150\;{\rm{kg} }$工作温度–30~ 50 ℃

      表 2  地基雷达LSA系统参数

      Table 2.  Ground-based radar LSA system parameters

      图  2  监测现场图

      Figure 2.  Monitoring site images

      图  3  实验数据

      Figure 3.  Experimental data

      采用均值对数比值算子和相干系数算子对图3(c)图3(d)的实验数据进行处理,得到均值对数比值图和相干系数图,如图4所示。将相干系数图和对数比值图经过NSCT处理融合(方法1)产生的合成差异图${F_h}$与通过邻域均值比和对数比值图经过NSCT处理融合(方法2)产生的差异图进行对比,如图5(a)图5(b)所示,其中colorbar代表图像像素值域,方法2为星载SAR图像变化检测生成差异图最常采用的方法[19]。将两种方法产生的差异图显示在3维图中,3维图如图5(c)图5(d)所示,X, Y轴代表图像像素的位置信息,Z轴代表差异度值越大发生变化的可能性越大,可以明显发现通过相干系数图和对数比值图经过NSCT处理融合的差异图群分性最好,这就代表了本文产生差异图的方法更容易进行聚类分割,进而区分变化区域与非变化区域。

      图  4  初始差异图

      Figure 4.  Initial difference images

      图  5  差异图对比

      Figure 5.  Difference images comparison

      对两时相的地基雷达图像进行变化检测,本文方法产生的合成差异图采用改进的FCM算法与传统的FCM算法的变形(包括RFCM和SCFCM)进行聚类结果对比,由于采用FCM算法进行聚类,无法得出聚类结果,在本文中不作分析。其结果如图6所示,图6(a)为RFCM聚类结果,图6(b)为SCFCM聚类结果,图6(c)为采用本文改进的FCM聚类算法得到的变化检测结果,即式(19)的变化检测结果$R$,其中蓝色的像素值为0,表示变化区域;黄色的像素值为1,表示未变化区域。

      图  6  聚类结果对比

      Figure 6.  Clustering result comparison

      图6(a)中可以明显地看到,RFCM算法聚类结果含有较多噪声点,无法聚类出结果。图6(b)中SCFCM算法去除较多噪声点,但依然有少量噪声存在,如图6(b)中红色矩形区域。而图6(c)中本文算法很好地滤除了周围的噪声点,得到的变化检测结果更加准确。通过实地考察,证实实验数据对应区域发生明显的滑坡,图7所示为监测现场光学图,图7(b)中黄色勾选区域发生了明显的滑坡。此外,图6(c)图6(b)相比,椭圆红框中检测到的变化区域存在较大差异,本文方法检测到该部分并没有发生明显的变化区域,在滑坡现场实地考察过程中发现滑坡前后此处为地质条件较为稳定的岩石,并没有随着滑坡的发生有所变化,只有两块岩石中间的部分小区域发生浮土下滑,如图7(b)中最小的黄色框区域,由此可确定检测到的变化结果与实地考察结果相符。地基雷达图像变化检测结果与地形3维数据融合结果如图8所示,根据融合图,可以准确地定位变化区域,如图8中红色框区域,通过融合图,进一步验证了本文方法的有效性,在实际应用中,可以为安全生产提供有效监测数据。

      图  7  监测现场光学图

      Figure 7.  Optical image of monitoring site

      图  8  变化检测结果与地形3维数据融合图

      Figure 8.  Effect of fusion of change detection results and UAV tilt photography

    • 对于地基雷达图像采用均值比以及对数比都无法得到群分性较好的差异图,本文提出使用相干系数与均值对数比得到群分性较好的差异图,使得地基雷达图像在变化检测过程中更容易进行聚类分割。此外,本文对传统FCM聚类算法进行改进,使图像中相类似的像素点以更大概率聚类到相同的类,以提高聚类效果并抑制噪声点。通过将所提出的方法应用于地基雷达实测数据,并将改进的FCM与FCM的另外两种变形RFCM和SCFCM聚类算法进行聚类结果对比,结果表明,本文所提方法具有更好的聚类效果,并且变化检测结果在保留变化区域的同时噪声点明显减少。未来的工作中,拟进一步探索地基雷达图像的无监督变化检测方法,以得到更为精确的变化检测结果,为安全生产、排险、救援提供数据参考。

参考文献 (21)

目录

    /

    返回文章
    返回