基于光学和SAR遥感图像融合的洪灾区域检测方法

王志豪 李刚 蒋骁

引用本文:
Citation:

基于光学和SAR遥感图像融合的洪灾区域检测方法

    作者简介:
    王志豪(1994–),男,清华大学电子工程系在读硕士生,目前研究方向雷达遥感图像处理;
    李 刚(1979–),男,2002年和2007年于清华大学电子系分别获得学士、博士学位,现为清华大学电子系教授、博士生导师,研究方向包括雷达成像、目标识别、稀疏信号处理、分布式信号处理等、信息融合等;
    蒋 骁(1989–),男,2010年和2013年分别于北京邮电大学通信工程系和清华大学电子工程系获得学士、硕士学位。现为清华大学电子工程系在读博士研究生。研究方向包括遥感图像融合与变化检测.
    通讯作者: 李刚 gangli@tsinghua.edu.cn
  • 基金项目:

    国家自然科学基金(61790551, 61925106),民用航天十三五预研项目(D010305)

  • 中图分类号: TN959.3

Flooded Area Detection Method Based on Fusion of Optical and SAR Remote Sensing Images

    Corresponding author: LI Gang, gangli@tsinghua.edu.cn
  • Fund Project: The National Natural Science Foundation of China (61790551 and 61925106), Civil Space Advance Research Program of China (D010305)

    CLC number: TN959.3

  • 摘要: 基于光学和合成孔径雷达(SAR)图像融合的洪灾区域检测方法可以全天候、高时效地检测洪灾区域。由于SAR图像中存在大量随机分布的相干斑噪声,传统洪灾区域检测方法的检测结果存在较高的虚警率。该文在模糊C均值聚类方法(FCM)的基础上提出了分级聚类算法(H-FCM),该方法将洪灾后的SAR图像与洪灾前的光学图像融合。基于融合图像,该方法利用提出的分级聚类模型获得洪灾区域的初步检测结果。此外,该算法在利用所提出的区域生长算法获得洪灾前河流位置后,将其作为初步检测结果的空间约束,进一步筛除疑似洪灾区域,并显著地提升了检测性能。该文的实验数据包括1999年英国格洛斯特洪灾前后的遥感图像和2019年中国南昌洪灾前后的遥感图像。通过对比实验,H-FCM算法的有效性得到验证。
  • 图 1  利奇马台风的GF4图像(资源卫星中心)

    Figure 1.  The GF4 image of Typhoon Lekima (Resource Satellite Centre)

    图 2  利奇马台风的GF3图像(资源卫星中心)

    Figure 2.  The GF3 image of Typhoon Lekima (Resource Satellite Centre)

    图 3  典型地物的光谱特征曲线[23]

    Figure 3.  Spectral characteristic curves of typical features[23]

    图 4  洪灾后的SAR图像Ea[10]

    Figure 4.  SAR image Ea after the flood[10]

    图 5  洪灾前的近红外波段影像Ei[10]

    Figure 5.  Near-infrared image Ei before the flood[10]

    图 6  方向矩阵Jd1(向下)

    Figure 6.  Direction matrix Jd1(down)

    图 7  方向矩阵Jd2(向上)

    Figure 7.  Direction matrix Jd2(up)

    图 8  方向矩阵Jd3(向左)

    Figure 8.  Direction matrix Jd3(left)

    图 9  方向矩阵Jd4(向右)

    Figure 9.  Direction matrix Jd4(right)

    图 10  第2类阈值图Pt

    Figure 10.  Threshold graph Pt

    图 11  4连通区域的定义

    Figure 11.  Definition of four-connected regions

    图 12  初始生长点位置

    Figure 12.  Initial growth point location

    图 13  初步生长结果

    Figure 13.  Preliminary growth results

    图 14  道路的识别结果

    Figure 14.  Road recognition results

    图 15  洪灾前的河流位置Eb

    Figure 15.  River location Eb before the flood

    图 16  融合后的灾后SAR图像Fu

    Figure 16.  Fusion post-disaster SAR image Fu

    图 17  聚类模型Gm7

    Figure 17.  Clustering model Gm7

    图 18  聚类模型的稀疏度分析

    Figure 18.  Analysis of sparsity of clustering model

    图 19  初步聚类结果Gr7

    Figure 19.  Preliminary clustering results Gr7

    图 20  洪灾区域空间约束曲线

    Figure 20.  Spatial constraint curve of flooded area

    图 21  最终检测结果FF

    Figure 21.  Final experimental result FF

    图 22  洪灾后的NDVI图像[10]

    Figure 22.  NDVI images after the flood[10]

    图 23  人工标注的真值图

    Figure 23.  Manually labeled truth map

    图 24  像素值分布3维图

    Figure 24.  3-D map of the distribution of pixel values

    图 25  Ea的边缘轮廓

    Figure 25.  Edge profile of Ea

    图 26  4种算法的洪灾区域检测结果

    Figure 26.  Detection results of flooded area based on four algorithms

    图 27  人工标注的洪灾区域

    Figure 27.  Manually labeled flooded areas

    图 28  洪灾前后哨兵1,哨兵2数据

    Figure 28.  Sentinel 1, 2 data before and after flood

    图 29  4种算法的洪灾区域检测结果

    Figure 29.  Detection results of flooded area based on four algorithms

    图 30  RL滤波后4种算法的洪灾区域检测结果

    Figure 30.  Detection results of flooded area based on four algorithms after RL filtering

    表 1  FCM算法

    Table 1.  FCM algorithm

     (1) 给定聚类中心数目c,隶属度因子m,迭代次数t
     (2) 初始化和为1的隶属度矩阵U
     (3) 由式(4)确定c个聚类中心Ci
     (4) 计算目标函数J的值;
     (5) 将Ci代入式(3),更新隶属度矩阵U
     (6)重复步骤(3)—步骤(5),迭代次数为t
    下载: 导出CSV

    表 2  区分河流和道路的区域生长算法

    Table 2.  Region growing algorithm to distinguish rivers and roads

     (1) 应用FCM算法获得图像的八分类聚类中心C
     (2) 获得阈值矩阵T,并获得小于T2的二值图Pt
     (3) 通过卷积消除Pt中的区块噪声,并确定初始生长点;
     (4) 设定初始生长方向矩阵$ {\bf{J}}{{\bf{d}}_i},i = 1,2,3,4$,如图6图9
     (5) 每次迭代都以初始生长点为中心,取5×5区域中的最小值点
       作为下一个生长点,按同一方向迭代直至生长到全图边界;
     (6) 将每个方向的生长结果累加,得到初步生长结果;
     (7) 利用霍夫变换[26]提取初步生长结果中的道路;
     (8) 将初步生长结果去除道路部分,获得河流位置。
    下载: 导出CSV

    表 3  融合后SAR图像洪灾区域的检测结果

    Table 3.  Detection results of flooded area in fusioned SAR image

    方法H-FCMWASnakeH-Kmeans
    Righta0.89460.43460.67570.8837
    Missa0.15230.06800.16070.1519
    Wra0.25130.70010.50700.2634
    Kappa0.60920.18400.29640.5837
    Time(s)51.982.722.5373.27
    Iter(次)251125
    下载: 导出CSV

    表 4  融合后SAR图像洪灾区域的检测结果

    Table 4.  Detection results of flooded area in the fusioned SAR image

    方法H-FCMWASnakeH-Kmeans
    Righta0.94140.62300.64260.9391
    Missa0.25200.05550.06080.2434
    Wra0.28860.46900.45180.2836
    Kappa0.69110.26290.27490.6848
    Time(s)30.592.823.1479.13
    Iter(次)251125
    下载: 导出CSV

    表 5  RL滤波后SAR图像的洪灾区域检测结果

    Table 5.  Detection results of flooded area in the SAR image after RL filtering

    方法RL_H-FCMRL_WARL_SnakeRL_H-Kmeans
    Righta0.89960.43910.76150.8899
    Missa0.17500.06600.25980.1351
    Wra0.26530.69330.49540.2418
    Kappa0.61470.18630.34600.6032
    Time(s)646.04593.72593.59745.62
    Iter(次)251125
    下载: 导出CSV
  • [1] SABINS F F. Remote sensing for mineral exploration[J]. Ore Geology Reviews, 1999, 14(3/4): 157–183. doi: 10.1016/S0169-1368(99)00007-4
    [2] 解金卫, 李真芳, 王帆, 等. 基于幅相不一致准则的建筑物SAR层析成像[J]. 雷达学报, 待出版. doi: 12000/JR19062.XIE Jinwei, LI Zhenfang, WANG Fan, et al. SAR tomography imaging for buildings using an inconsistency criterion for amplitude and phase[J]. Journal of Radars, in press, doi: 12000/JR19062.
    [3] BERNI J A J, ZARCO-TEJADA P J, SUÁREZ L, et al. Thermal and narrowband multispectral remote sensing for vegetation monitoring from an unmanned aerial vehicle[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(3): 722–738. doi: 10.1109/TGRS.2008.2010457
    [4] SAWAYA K E, OLMANSON L G, HEINERT N J, et al. Extending satellite remote sensing to local scales: Land and water resource monitoring using high-resolution imagery[J]. Remote Sensing of Environment, 2003, 88(1/2): 144–156. doi: 10.1016/j.rse.2003.04.006
    [5] DONG Laigen and SHAN Jie. A comprehensive review of earthquake-induced building damage detection with remote sensing techniques[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 84: 85–99. doi: 10.1016/j.isprsjprs.2013.06.011
    [6] SANYAL J and LU X X. Application of remote sensing in flood management with special reference to monsoon Asia: A review[J]. Natural Hazards, 2004, 33(2): 283–301. doi: 10.1023/B:NHAZ.0000037035.65105.95
    [7] 吴一全, 王志来. 基于联合稀疏表示的复Contourlet域SAR图像与红外图像融合(英文)[J]. 雷达学报, 2017, 6(4): 349–358. doi: 10.12000/JR17019WU Yiquan and WANG Zhilai. SAR and infrared image fusion in complex contourlet domain based on joint sparse representation[J]. Journal of Radars, 2017, 6(4): 349–358. doi: 10.12000/JR17019
    [8] BELL J R, SCHULTZ L A, JONES M, et al. Using optical remote sensing and synthetic aperture radar for near-real-time response to the central U.S. flooding in April-May 2017[C]. The 98th 2018 American Meteorological Society Meeting, Austin, Texas, 2018.
    [9] KUSSUL N, SHELESTOV A, and SKAKUN S. Flood Monitoring from SAR Data[M]. KOGAN F, POWELL A, and FEDOROV O. Use of Satellite and In-Situ Data to Improve Sustainability. Dordrecht: Springer, 2011: 19-29. doi: 10.1007/978-90-481-9618-0_3.
    [10] LIU Zhunga, ZHANG Li, LI Gang, et al. Change detection in heterogeneous remote sensing images based on the fusion of pixel transformation[C]. The 2017 20th International Conference on Information Fusion, Xi’an, China, 2017: 1–6. doi: 10.23919/ICIF.2017.8009656.
    [11] LIU Zhunga, LI Gang, MERCIER G, et al. Change detection in heterogenous remote sensing images via homogeneous pixel transformation[J]. IEEE Transactions on Image Processing, 2018, 27(4): 1822–1834. doi: 10.1109/TIP.2017.2784560
    [12] 李丹, 吴保生, 陈博伟, 等. 基于卫星遥感的水体信息提取研究进展与展望[J/OL]. 清华大学学报: 自然科学版. http://kns.cnki.net/kcms/detail/detail.aspx?doi=10.16511/j.cnki.qhdxxb.2019.22.038, 2020.LI Dan, WU Baosheng, CHEN Bowei, et al. Review of water body information extraction based on satellite remote sensing[J/OL]. Journal of Tsinghua University: Science and Technology. http://kns.cnki.net/kcms/detail/detail.aspx?doi=10.16511/j.cnki.qhdxxb.2019.22.038, 2020.
    [13] OTSU N. A threshold selection method from gray-level histograms[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1979, 9(1): 62–66. doi: 10.1109/TSMC.1979.4310076
    [14] 窦建方, 陈鹰, 翁玉坤. 基于序列非线性滤波SAR影像水体自动提取[J]. 测绘通报, 2008(9): 37–39, 45.DOU Jianfang, CHEN Ying, and WENG Yukun. Automatic water body extraction from SAR images based on sequence non-linear filter[J]. Bulletin of Surveying and Mapping, 2008(9): 37–39, 45.
    [15] LEE J S and JURKEVICH I. Coastline detection and tracing in SAR images[J]. IEEE Transactions on Geoscience and Remote Sensing, 1990, 28(4): 662–668. doi: 10.1109/TGRS.1990.572976
    [16] DESCOMBES X, MOCTEZUMA M, MAÎTRE H, et al. Coastline detection by a Markovian segmentation on SAR images[J]. Signal Processing, 1996, 55(1): 123–132. doi: 10.1016/S0165-1684(96)00125-9
    [17] 滑文强, 王爽, 郭岩河, 等. 基于邻域最小生成树的半监督极化SAR图像分类方法[J]. 雷达学报, 2019, 8(4): 458–470. doi: 10.12000/JR18104HUA Wenqiang, WANG Shuang, GUO Yanhe, et al. Semi-supervised PolSAR image classification based on the neighborhood minimum spanning tree[J]. Journal of Radars, 2019, 8(4): 458–470. doi: 10.12000/JR18104
    [18] 赵娟萍, 郭炜炜, 柳彬, 等. 基于概率转移卷积神经网络的含噪标记SAR图像分类[J]. 雷达学报, 2017, 6(5): 514–523. doi: 10.12000/JR16140ZHAO Juanping, GUO Weiwei, LIU Bin, et al. Convolutional neural network-based SAR image classification with noisy labels[J]. Journal of Radars, 2017, 6(5): 514–523. doi: 10.12000/JR16140
    [19] HARTIGAN J A and WONG M A. Algorithm AS 136: A k-means clustering algorithm[J]. Journal of the Royal Statistical Society. Series C (Applied Statistics) , 1979, 28(1): 100–108.
    [20] 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
    [21] VANNOTE R L, MINSHALL G W, CUMMINS K W, et al. The river continuum concept[J]. Canadian Journal of Fisheries and Aquatic Sciences, 1980, 37(1): 130–137. doi: 10.1139/f80-017
    [22] 方佳佳, 王烜, 孙涛, 等. 河流连通性及其对生态水文过程影响研究进展[J]. 水资源与水工程学报, 2018, 29(2): 19–26. doi: 10.11705/j.issn.1672-643X.2018.02.04FANG Jiajia, WANG Xuan, SUN Tao, et al. Review of research on river connectivity and its impact on eco -hydrological process[J]. Journal of Water Resources and Water Engineering, 2018, 29(2): 19–26. doi: 10.11705/j.issn.1672-643X.2018.02.04
    [23] 童庆禧, 田国良. 中国典型地物波谱及其特征分析[M]. 北京: 科学出版社, 1990.TONG Qingxi and TIAN Guoliang. Spectra and Analysis of Typical Earth Objects of China[M]. Beijing: Science Press, 1990.
    [24] CURCIO J A and PETTY C C. The near infrared absorption spectrum of liquid water[J]. Journal of the Optical Society of America, 1951, 41(5): 302–304. doi: 10.1364/JOSA.41.000302
    [25] 章毓晋. 图像工程[M]. 4版. 北京: 清华大学出版社, 2018: 51–53.ZHANG Yujin. Image Engineering[M]. 4th ed. Beijing: Tsinghua University Press, 2018: 51–53.
    [26] DUDA R O and HART P E. A Generalized Hough Transformation for Detecting Lines in Pictures[M]. Artificial Intelligence Group, SRI International, 1970.
    [27] XU Hanqiu. Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery[J]. International Journal of Remote Sensing, 2006, 27(14): 3025–3033. doi: 10.1080/01431160600589179
    [28] LI Gang and BURKHOLDER R J. Hybrid matching pursuit for distributed through-wall radar imaging[J]. IEEE Transactions on Antennas and Propagation, 2015, 63(4): 1701–1711. doi: 10.1109/TAP.2015.2398115
    [29] WANG Xueqian, LI Gang, QUAN Chen, et al. Distributed detection of sparse stochastic signals with quantized measurements: The generalized Gaussian case[J]. IEEE Transactions on Signal Processing, 2019, 67(18): 4886–4898. doi: 10.1109/TSP.2019.2932884
    [30] BATES P D, HORRITT M S, ARONICA G, et al. Bayesian updating of flood inundation likelihoods conditioned on flood extent data[J]. Hydrological Processes, 2004, 18(17): 3347–3370. doi: 10.1002/hyp.1499
    [31] MCEWEN L J, KRAUSE F, JONES O, et al. Sustainable flood memories, informal knowledge and the development of community resilience to future flood risk[J]. WIT Transactions on Ecology and the Environment, 2012, 159(12): 253–264. doi: 10.2495/FRIAR120211
    [32] RAMBABU C, CHAKRABARTI I, and MAHANTA A. Flooding-based watershed algorithm and its prototype hardware architecture[J]. IEE Proceedings-Vision, Image and Signal Processing, 2004, 151(3): 224–234. doi: 10.1049/ip-vis:20040397
    [33] DE ROO A, VAN DER KNIJFF J, HORRITT M, et al. Assessing flood damages of the 1997 Oder flood and the 1995 Meuse flood[C]. The 2nd International Symposium on Operationalization of Remote Sensing, Enschede, The Netherlands, 1999: 16–20.
    [34] LEE J S. Refined filtering of image noise using local statistics[J]. Computer Graphics and Image Processing, 1981, 15(4): 380–389. doi: 10.1016/S0146-664X(81)80018-4
  • [1] 赵逸超朱宇涛粟毅杨猛 . 用于线阵三维SAR成像的二维快速ESPRIT算法. 雷达学报, doi: 10.12000/JR15065
    [2] 邢孟道林浩陈溅来孙光才严棒棒 . 多平台合成孔径雷达成像算法综述. 雷达学报, doi: 10.12000/JR19102
    [3] 李晓峰张彪杨晓峰 . 星载合成孔径雷达遥感海洋风场波浪场. 雷达学报, doi: 10.12000/JR20079
    [4] 王松张福博陈龙永梁兴东 . 基于空间聚类种子生长算法的阵列干涉SAR点云滤波. 雷达学报, doi: 10.12000/JR18006
    [5] 周春晖李飞李宁郑慧芳王翔宇 . 改进的基于特征子空间的SAR图像射频干扰抑制算法. 雷达学报, doi: 10.12000/JR17025
    [6] 黄岩赵博陶明亮陈展野洪伟 . 合成孔径雷达抗干扰技术综述. 雷达学报, doi: 10.12000/JR19113
    [7] 卫扬铠曾涛陈新亮丁泽刚范宇杰温育涵 . 典型线面目标合成孔径雷达参数化成像. 雷达学报, doi: 10.12000/JR19077
    [8] 毛永飞汪小洁向茂生 . 机载干涉SAR 区域网三维定位算法. 雷达学报, doi: 10.3724/SP.J.1300.2012.20107
    [9] 杜康宁邓云凯王宇李宁 . 基于多层神经网络的中分辨SAR图像时间序列建筑区域提取. 雷达学报, doi: 10.12000/JR16060
    [10] 吴一戎洪文张冰尘蒋成龙张柘赵曜 . 稀疏微波成像研究进展(科普类). 雷达学报, doi: 10.3724/SP.J.1300.2014.14105
    [11] 金添 . 叶簇穿透合成孔径雷达增强成像方法. 雷达学报, doi: 10.12000/JR15114
    [12] 吴玉峰叶少华冯大政 . 基于方位相位编码的脉内聚束SAR成像方法. 雷达学报, doi: 10.12000/JR17114
    [13] 唐江文邓云凯王宇赵硕李宁 . 高分辨率滑动聚束SAR BP成像及其异构并行实现. 雷达学报, doi: 10.12000/JR16053
    [14] 赵团邓云凯王宇李宁王翔宇 . 基于扇贝效应校正的改进滑动Mosaic全孔径成像算法. 雷达学报, doi: 10.12000/JR16014
    [15] 任笑真杨汝良 . 一种基于幅度和相位迭代重建的四维合成孔径雷达成像方法. 雷达学报, doi: 10.12000/JR15135
    [16] 张金松邢孟道孙光才 . 一种基于密集深度分离卷积的SAR图像水域分割算法. 雷达学报, doi: 10.12000/JR19008
    [17] 张文彬邓云凯王宇 . 星地双基合成孔径雷达聚束模式快速BP算法. 雷达学报, doi: 10.3724/SP.J.1300.2013.13031
    [18] 张问一胡东辉丁赤飚 . 基于FABEMD 和Goldstein 滤波器的SAR 舰船尾迹图像增强方法. 雷达学报, doi: 10.3724/SP.J.1300.2012.20059
    [19] 张新征谭志颖王亦坚 . 基于多特征-多表示融合的SAR图像目标识别. 雷达学报, doi: 10.12000/JR17078
    [20] 刘永才王伟潘小义代大海 . 基于延迟-移频的SAR 有源欺骗干扰有效区域研究. 雷达学报, doi: 10.3724/SP.J.1300.2012.13001
  • 加载中
图(30)表(5)
计量
  • 文章访问数:  753
  • HTML浏览量:  424
  • PDF下载量:  140
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-11-05
  • 录用日期:  2020-02-18
  • 网络出版日期:  2020-03-07

基于光学和SAR遥感图像融合的洪灾区域检测方法

    通讯作者: 李刚 gangli@tsinghua.edu.cn
    作者简介:
    王志豪(1994–),男,清华大学电子工程系在读硕士生,目前研究方向雷达遥感图像处理;
    李 刚(1979–),男,2002年和2007年于清华大学电子系分别获得学士、博士学位,现为清华大学电子系教授、博士生导师,研究方向包括雷达成像、目标识别、稀疏信号处理、分布式信号处理等、信息融合等;
    蒋 骁(1989–),男,2010年和2013年分别于北京邮电大学通信工程系和清华大学电子工程系获得学士、硕士学位。现为清华大学电子工程系在读博士研究生。研究方向包括遥感图像融合与变化检测
  • 清华大学电子工程系 北京 100084
基金项目:  国家自然科学基金(61790551, 61925106),民用航天十三五预研项目(D010305)

摘要: 基于光学和合成孔径雷达(SAR)图像融合的洪灾区域检测方法可以全天候、高时效地检测洪灾区域。由于SAR图像中存在大量随机分布的相干斑噪声,传统洪灾区域检测方法的检测结果存在较高的虚警率。该文在模糊C均值聚类方法(FCM)的基础上提出了分级聚类算法(H-FCM),该方法将洪灾后的SAR图像与洪灾前的光学图像融合。基于融合图像,该方法利用提出的分级聚类模型获得洪灾区域的初步检测结果。此外,该算法在利用所提出的区域生长算法获得洪灾前河流位置后,将其作为初步检测结果的空间约束,进一步筛除疑似洪灾区域,并显著地提升了检测性能。该文的实验数据包括1999年英国格洛斯特洪灾前后的遥感图像和2019年中国南昌洪灾前后的遥感图像。通过对比实验,H-FCM算法的有效性得到验证。

English Abstract

    • 遥感技术为地理勘探[1]、建筑测绘[2]、植被监测[3]、水域监测[4]、灾区探测[5]等提供了重要的技术支持。在洪灾区域检测中,光学遥感能够提供丰富的光谱信息和较高的图像分辨率,而雷达遥感则具备全天候的洪灾区域检测能力[6]。异质图像融合已经在复Contourlet域得到了实现[7],因此结合光学遥感数据和雷达遥感数据能实现洪灾区域的全天候、高时效检测[8]

      当前,学者们提出许多基于遥感图像的洪灾区域检测方法。其中,应用人工神经网络和自组织映射网络[9]方法处理洪灾前后的同质SAR图像时能大致区分洪灾区域,但是受到大量相干斑噪声的影响,检测结果仍然存在很高的虚警率。为了结合光学遥感图像和SAR图像的成像优势,一些学者提出了基于异质图像的洪灾区域检测方法,例如像素变换法[10]和特征空间变换法[11]。上述方法对洪灾前的SAR图像和洪灾后的光学图像有良好的检测效果,但是对洪灾前的光学图像和洪灾后的SAR图像的检测效果较差。

      针对SAR图像中洪灾区域的检测难点,本文提出了H-FCM算法。该算法的主要贡献如下:

      (1) 提出SAR图像的分级聚类模型;

      (2) 将洪灾前的近红外波段光学遥感图像中提取出的河流位置融合到洪灾后的SAR图像中,并作为洪灾区域检测的空间约束,进一步提升了检测结果的准确率。

      文章结构如下:第2节分析了近年来洪灾区域检测方法的发展和存在的问题。第3节介绍了分级C均值聚类算法Hierarchical Fuzzy C-Means (H-FCM)算法的原理、步骤以及参数设定对算法性能的影响。第4节应用包括H-FCM算法的4种算法对不同类型遥感数据进行洪灾区域检测,并分析了图像降噪处理对4种算法检测性能的影响。第5节系统总结了全文的主要工作和存在的不足。

    • 目前,学者们提出的基于雷达遥感图像的水体检测方法包括灰度阈值分割法、滤波法、机器学习法和结合辅助信息的提取方法等[12];基于光学遥感图像的提取水体方法包括基于像元分类的阈值法和基于目标分类的分类法等[12]。洪灾是由于大雨、暴雨或者持续降雨造成的,此时无法通过实时获取的光学遥感数据检测地面的洪灾区域。

      以2019年8月9日对利奇马台风实时追踪的GF4数据(见图1)为例,光学图像无法显示地表的洪灾区域。但是通过适当的卫星调度方案可以短时间内获取雷达遥感数据。以2019年8月8日对利奇马台风实时追踪的GF3数据(见图2)为例,SAR图像可以清晰地显示地面受灾情况。

      图  1  利奇马台风的GF4图像(资源卫星中心)

      Figure 1.  The GF4 image of Typhoon Lekima (Resource Satellite Centre)

      图  2  利奇马台风的GF3图像(资源卫星中心)

      Figure 2.  The GF3 image of Typhoon Lekima (Resource Satellite Centre)

      应用阈值分割法中的最大类间方差算法[13]对灾后的SAR图像进行洪灾区域检测时,由于受到相干斑噪声的影响,检测精度较低。为克服相干斑噪声的影响,学者们提出了滤波法[14]、边界追踪法[15]、Markov分割法[16]、基于邻域最小生成树的半监督分类法[17]、基于概率转移卷积神经网络的分类法[18]等。但是洪灾区域检测问题的核心在于保证高检测准确率的同时,尽可能缩短检测时间,而上述方法的计算复杂度普遍较高,并不能满足洪灾区域检测的高时效性需求。

      直接利用K-means算法[19]、C均值聚类方法(FCM)算法[20]等无监督聚类算法对洪灾后的SAR图像进行洪灾区域检测可以缩短检测时间,但是检测结果的分类准确率很低。具体原因如下:根据灰度阈值分割法的原理,在SAR图像中,水体的散射值较低,整体呈现为暗区。在没有发生暴雨或者持续降雨的情况下,可以通过分析SAR图像的整体像素值分布,应用FCM算法等聚类算法确定水体的分割阈值。根据分割阈值,将大于阈值的标记为非水体,小于阈值的标记为水体。但是洪水灾害发生时伴随着暴雨或者持续降雨,导致很多未受到洪灾影响的区域在SAR图像中也呈现为暗区。因此,直接通过阈值法获得的检测结果存在较高的虚警率。

    • FCM算法的步骤总结为表1。以洪灾后的大小为n1×n2像素的SAR图像为例,将其所有像素点的像素值记作2维矩阵X。将X分成c类,则可以得到c个聚类中心。将每一个像素点作为单个样本xj,则xj对于第i类聚类中心的隶属度可记作uij, FCM算法的目标函数、约束条件的定义式分别为式(1)和式(2)

       (1) 给定聚类中心数目c,隶属度因子m,迭代次数t
       (2) 初始化和为1的隶属度矩阵U
       (3) 由式(4)确定c个聚类中心Ci
       (4) 计算目标函数J的值;
       (5) 将Ci代入式(3),更新隶属度矩阵U
       (6)重复步骤(3)—步骤(5),迭代次数为t

      表 1  FCM算法

      Table 1.  FCM algorithm

      $ J = {\rm{}}\sum\limits_{i = 1}^c \sum\limits_{j = 1}^n u_{ij}^m{\left\| {{x_j} - {C_i}} \right\|^2} $

      $ \sum\limits_{i = 1}^c {u_{ij}} = 1 \hspace{70pt} $

      其中,Ci表示第i类聚类中心;m表示隶属度因子;n=n1×n2; $ i = 1,2, ··· ,c$; $ j = 1,2, ··· ,n$

      结合式(2)的约束条件,利用拉格朗日乘数法,将J分别对uijCi求导,并令其结果等于0,可得uijCi的迭代式分别为式(3)和式(4)

      $ {u_{ij}} = \frac{1}{{\displaystyle\sum\limits_{k = 1}^c {{{\left( {\frac{{\left\| {{x_j} - {C_i}} \right\|}}{{\left\| {{x_j} - {C_k}} \right\|}}} \right)}^{\frac{2}{{m - 1}}}}} }} $

      $ {C_i} = \frac{{\displaystyle\sum\limits_{j = 1}^n \left( {{x_j}u_{ij}^m} \right)}}{{\displaystyle\sum\limits_{j = 1}^n u_{ij}^m}} \hspace{45pt} $

    • 本文提出的H-FCM算法致力于缩短SAR图像洪灾区域的检测时间、降低洪灾区域的错误检测概率以及抑制SAR图像中存在的大量随机分布的相干斑噪声对检测结果的影响。根据河流连续体(River Continuum Concept, RCC)[21]的概念,在已知的河流系统内,自上游到下游的物理变量的梯度变化是连续的[22]。因此发生洪水灾害之后,若地理环境没有发生剧烈变化,洪水的形态仍然符合RCC。并且,根据典型地物的光谱特征曲线图[23](见图3),水体对近红外波段的电磁波吸收能力比大部分地物强。

      图  3  典型地物的光谱特征曲线[23]

      Figure 3.  Spectral characteristic curves of typical features[23]

      H-FCM算法的应用前提为可获得洪灾前的光学图像以及洪灾后的SAR图像。该算法通过对遥感图像的像素强度划分进行分级聚类。以uint8格式的图像数据为例,应用FCM算法获得待测数据像素值的8个聚类中心,并且由低至高排列。根据图3, SAR图像中代表水体的像素点的像素值主要分布在0~64之间,即处于低像素值区间。但是由于发生洪灾时,SAR图像中存在大量随机分布的相干斑噪声以及低洼地区存在的积水,直接用阈值法进行洪灾区域检测的准确率较低。

      H-FCM算法采用反向聚类的方法,考虑到在SAR图像中,代表水体的像素点几乎不可能处于高像素值区间。H-FCM算法在聚类过程中将处于最高级的像素值区间的像素点记为1,其余像素点记为0,对所有零像素点进行区域聚类,获得初步的洪灾区域检测结果。又由RCC原理,洪灾发生后,主要受灾区域与洪灾前的河流相连通,利用灾害前的河流位置作为空间约束可以得到近似真实洪灾区域的检测结果。因此,应用H-FCM算法检测洪灾发生后的SAR图像的洪灾区域,可以降低相干斑噪声对检测结果的影响。

      H-FCM算法主要分为4步。首先,应用FCM算法获得洪灾后SAR图像的聚类中心,进而设置分类阈值,并确定分级聚类模型;其次,提取灾害前近红外波段光学遥感数据中的河流位置,并将其融合到洪灾后的SAR图像中;再次,以聚类模型为基础进行分级聚类,获得洪灾区域的初步聚类结果;最后,通过空间约束矩阵对初步聚类结果进行约束,进一步提升洪灾区域检测的准确率。

    • 利用FCM算法获得洪灾后SAR图像的八分类的聚类中心矩阵C和阈值矩阵T,并由此获得7个聚类模型

      $ {{C}} = \left[ {{C_1}{C_2} ··· {C_8}} \right] $

      其中,1维矩阵C表示聚类中心矩阵,包含8个按照升序排列的聚类中心

      $ {{T}} = \left[ {{T_1}{T_2} ··· {T_7}} \right] \hspace{11pt} $

      $ {T_q} = \frac{{({C_q} + {C_{q + 1}})}}{2} $

      其中,$ {T_0} = 0,{\rm{}}q = 1, 2,··· ,7$

      以1999年英国格洛斯特洪灾后获得的SAR图像Ea[10]为例,其大小为n1×n2(见图4)。遍历Ea中的所有像素点,若其像素值大于Tq–1且小于Tq则在对应位置标记为1,否则标记为0,将所得到的矩阵记作聚类模型Gmq,如式(8)

      图  4  洪灾后的SAR图像Ea[10]

      Figure 4.  SAR image Ea after the flood[10]

      $ {\left( {{\bf{G}}{{\bf{m}}_q}} \right)_{mn}} = \left\{ \begin{aligned} & {1,\; {T_{q - 1}} < {{\left( {{{{E}}_{\rm{a}}}} \right)}_{mn}} < {T_q}}\\ & {0,\; {\rm{otherwise}}} \end{aligned} \right. $

      其中,“(Gmq)mnGmq”在本文中均代表(Gmq)mnGmq矩阵中第m行,第n列的元素,下文符号“∈”的含义与此处一致;(Ea)mnEa; $ q = 1, 2,··· ,7$

    • 由于SAR图像的成像模式一般为宽幅扫描模式,同时受到SAR图像中相干斑噪声的影响,在没有灾害前河流信息的情况下直接采用无监督聚类的方法来检测Ea中的洪灾区域存在较高的虚警率,导致检测结果的准确率较低。H-FCM算法将洪灾前河流的空间信息融合到洪灾后的SAR图像中,通过局部空间约束,可以提高洪灾区域检测的准确率。但是在图4中,Ea的河流位置有一部分被背景地物或者噪声覆盖,另一部分和洪灾区域重合,导致大部分的河流位置已经无法识别。为了获得洪灾前的河流位置信息,H-FCM算法进一步检测洪灾前的光学图像。

      以1999年英国格洛斯特洪灾前的SPOT图像中提取的近红外波段图像Ei[10]为例(见图5)。由于不同河段的河流泥沙含量不同,造成不同河段对近红外波段电磁波的吸收能力不同[24],通过阈值法无法直接获得所有河流区域的空间位置。考虑到洪灾前的水体表面近似为平滑面,可以应用区域生长算法[25]提取洪灾前近红外波段遥感图像的河流位置。本文提出的区域生长算法为表2,根据图5中近红外波段的光学遥感图像可以利用该算法获得洪灾前的河流位置。

       (1) 应用FCM算法获得图像的八分类聚类中心C
       (2) 获得阈值矩阵T,并获得小于T2的二值图Pt
       (3) 通过卷积消除Pt中的区块噪声,并确定初始生长点;
       (4) 设定初始生长方向矩阵$ {\bf{J}}{{\bf{d}}_i},i = 1,2,3,4$,如图6图9
       (5) 每次迭代都以初始生长点为中心,取5×5区域中的最小值点
         作为下一个生长点,按同一方向迭代直至生长到全图边界;
       (6) 将每个方向的生长结果累加,得到初步生长结果;
       (7) 利用霍夫变换[26]提取初步生长结果中的道路;
       (8) 将初步生长结果去除道路部分,获得河流位置。

      表 2  区分河流和道路的区域生长算法

      Table 2.  Region growing algorithm to distinguish rivers and roads

      图  5  洪灾前的近红外波段影像Ei[10]

      Figure 5.  Near-infrared image Ei before the flood[10]

    • 图  6  方向矩阵Jd1(向下)

      Figure 6.  Direction matrix Jd1(down)

      图  7  方向矩阵Jd2(向上)

      Figure 7.  Direction matrix Jd2(up)

      图  8  方向矩阵Jd3(向左)

      Figure 8.  Direction matrix Jd3(left)

      图  9  方向矩阵Jd4(向右)

      Figure 9.  Direction matrix Jd4(right)

      由于步骤(2)获得的第2类阈值图Pt(见图10)中存在大量的区块噪声的干扰,无法从中直接获得灾害前的具体河流的空间位置。因此,需要筛除大面积区块性的非河流区域。像素点的面积定义为:将像素值为1的像素点的面积记作1,像素值为0的像素点面积记作0。连通区域的面积定义为:将连通区域内所有非零像素点的面积累加得到该区域的面积。4连通区域的定义为目标像素点(红色方块)的上、下、左、右4个像素点位置(见图11)。

      图  10  第2类阈值图Pt

      Figure 10.  Threshold graph Pt

      图  11  4连通区域的定义

      Figure 11.  Definition of four-connected regions

      为了充分消除区块噪声的影响,并保留河流区域的完整性,步骤(3)的卷积操作分为两步。首先,确定大小为α×α的全1卷积核Wα,将PtWα卷积得到初步消除区块噪声的结果P1,如式(9)和式(10);其次,将P1Wβ卷积得到基本消除区块噪声的结果P2,如式(11)和式(12)。为保证卷积之后的输出结果的尺寸与原图像大小保持一致,本文所有的卷积模式均采取matlab中conv2函数的‘same’模式进行卷积

      $ {{{P}}_{{\rm{t}}1}} = {{{P}}_{\rm{t}}}*{{{W}}_{\rm{\alpha }}} $

      其中,“*”表示以‘same’模式进行卷积

      $ {\left( {{{{P}}_1}} \right)_{ij}} = \left\{ \begin{aligned} & {1,\; {{\left( {{{{P}}_{{\rm{t}}1}}} \right)}_{ij}} < 0.6 \cdot {{\rm{\alpha }}^2}}\\ & { 0,\; {\rm{otherwise}}} \end{aligned} \right. $

      其中,(Pt1)ijPt1; (P1)ijP1

      $ {{{P}}_{{\rm{t}}2}} = {{{P}}_1}*{{{W}}_{\rm{\beta }}} $

      其中,$ {\rm{\beta }} ={\alpha }/{2}$, Wβ为大小为β×β的全1卷积核

      $ {\left( {{{{P}}_2}} \right)_{ij}} = \left\{ \begin{aligned} & {1,\; {{\left( {{{{P}}_{{\rm{t}}2}}} \right)}_{ij}} < 0.6 \cdot {{\rm{\beta }}^2}}\\ & {0,\; {\rm{otherwise}}} \end{aligned} \right. $

      其中,(P2)ijP2, (Pt2)ijPt2

    • 参数α一般为固定的大小,当选择的α较大或者较小时,在筛除区块噪声的同时,也筛除了河流。本文固定参数α=18的卷积核Wα初步筛除大面积聚点,获得P1;其次,固定参数β=9的卷积核Wβ进一步筛除小面积聚点,获得P2。利用式(13)至式(15)确定河流区域种子生长点的位置Po,在Po的每一个4连通区域中等间距选择10个像素点作为初始生长点(见图12)。

      图  12  初始生长点位置

      Figure 12.  Initial growth point location

      $ {{{P}}_{\rm{o}}} = \sum\limits_{j = 1}^{{n_{\rm{a}}}} \left( {{\rm{de}}{{\rm{t}}_j} \cdot {\bf{Ob}}{{\bf{j}}_j}} \right) $

      其中,Objj表示P2中的第j个4连通区域;P2中4连通区域总数为na;每个连通区域的背景大小均为n1×n2; $ j = 1,2,··· ,{n_{\rm{a}}}$

      $ {\rm{de}}{{\rm{t}}_j} = \left\{ {\begin{array}{*{20}{c}} {1,{\rm{Are}}{{\rm{a}}_j} > 500}\\ {0,{\rm{Are}}{{\rm{a}}_j} \le 500} \end{array}} \right. $

      $ {\rm{Are}}{{\rm{a}}_j} = \sum\limits_{m = 1}^{{n_1}} \sum\limits_{q = 1}^{{n_2}} {({\bf{Ob}}{{\bf{j}}_j})_{mq}}\hspace{7pt} $

      其中,Areaj表示P2中第j个4连通区域Objj的面积。

    • 根据图12中获得初始生长点,按照步骤(5)进行4个方向的区域生长,最后将各个方向的生长结果累加得到初步生长结果(见图13)。但是由于本例中道路和河流交错,当道路和河流宽度相近且没有先验测绘信息时,直接通过机器识别提取并区分河流和道路的空间位置是非常困难的。所以在步骤(7)中利用霍夫变换[26]检测初步生长结果中属于公路的长距离直线段区域,进而识别道路所在区域(见图14)。最后,将初步生长结果去除道路区域后,得到河流位置Eb(见图15)。

      图  13  初步生长结果

      Figure 13.  Preliminary growth results

      图  14  道路的识别结果

      Figure 14.  Road recognition results

      图  15  洪灾前的河流位置Eb

      Figure 15.  River location Eb before the flood

    • 将洪灾前的河流位置Eb和洪灾后的SAR图像Ea融合。融合规则为将Ea中与Eb中灾前河流位置对应的所有像素点的像素值记作T1,其余像素点的像素值保持不变,如式(16),将融合后的灾后SAR图像记作Fu(见图16)

      图  16  融合后的灾后SAR图像Fu

      Figure 16.  Fusion post-disaster SAR image Fu

      $ {\left( {{\bf{Fu}}} \right)_{ij}} = \left\{ {\begin{array}{*{20}{c}} {{T_1},\;\;\;\;\;\;\;\;{{({{{E}}_{\rm{b}}})}_{ij}} = 1}\\ {{{\left( {{{{E}}_{\rm{a}}}} \right)}_{ij}},\;\;\;\;{\rm{otherwise}}} \end{array}} \right. $

      其中,(Fu)ijFu, (Eb)ijEb, (Ea)ijEa

    • 光学图像的水体和非水体的区分可以通过计算归一化水指数(Normalised Difference Water Index, NDWI)[27]检测,但是SAR图像中不含多通道图像信息,不能直接计算其NDWI来区分水体和非水体。通过观察融合后SAR图像的高像素值区间Gm7(见图17),可以发现真实洪灾区域的像素点在其中几乎全部表现为零像素值。同时,根据RCC原理,真实洪灾区域在整个图像空间呈现较高的区块性的稀疏度。即代表洪灾区域的非零像素值区域是一个连通区域,其余零像素值区域均为非洪灾区域,所以作为洪灾区域检测的聚类模型应当具备较高的稀疏度。Gm7是一个反向模型,洪灾区域基本呈现为零像素区域,非洪灾区域基本呈现为非零像素的散点分布。

      图  17  聚类模型Gm7

      Figure 17.  Clustering model Gm7

      在目标检测中,文献[28]定义图像稀疏度为非零像素点个数K1,考虑到洪灾区域检测的目标地物覆盖范围较广,本文进一步引入文献[29]中图像稀疏占比的定义。即模型的稀疏占比为待测图像中非零像素点个数和所有像素点个数的比例,稀疏占比K2的定义式为公式(17)。聚类模型的稀疏度分析,以Gm7聚类模型为例,检测结果如图18

      图  18  聚类模型的稀疏度分析

      Figure 18.  Analysis of sparsity of clustering model

      $ {K_1} = \sum\limits_{i = 1}^{{n_1}} \sum\limits_{j = 1}^{{n_2}} {\left( {{\bf{G}}{{\bf{m}}_7}} \right)_{ij}},{K_2} = \frac{{{K_1}}}{{{n_1} \times {n_2}}} $

      通过分析聚类模型Gm7的稀疏度,得到其稀疏占比为0.10905,可以确认该模型本身具有较高的稀疏度。并且该模型的非零像素点代表最高级像素值的分布,实际洪灾区域在该模型中基本处于零像素值区域。对聚类模型Gm7进行卷积运算后,根据卷积结果${\bf {Gt}}_7^k $得到每一个像素点周围非零像素点的总数,并由此判断该像素点是否为洪灾区域。

      $ {\bf{Gt}}_7^k = {\bf{G}}{{\bf{m}}_7}*{{{W}}_k} $

      其中,Wk表示大小为k×k的全1卷积核。

      H-FCM算法对聚类系数$\varphi $和卷积系数k的选择较为敏感。当选择的$\varphi $过大时,聚类结果包含较多的非洪灾区域;当选择的$\varphi $较小时,聚类结果舍弃了较多真实洪灾区域。相应地,当选择的k过大或者过小时,聚类结果都将包含较多的非洪灾区域。为使这两类参数选择达到最优检测性能,应当在选择最佳聚类模型的同时,依据模型的稀疏度来确定相应的聚类系数和卷积系数。考虑到SAR图中大量随机分布的相干斑噪声,为减弱噪声对检测性能的影响,$ \varphi ,k$的选择要比实际的稀疏占比更小,$ \varphi ,k$K2的关系定义为式(19)。

      $ \varphi = \frac{{{K_2}}}{2},k = \varphi \cdot \sqrt {{n_1} \times {n_2}} \hspace{35pt} $

      $ {\left( {{\bf{G}}{{\bf{r}}_7}} \right)_{ij}} = \left\{ {\begin{array}{*{20}{c}} {1,{{\left( {{\bf{Gt}}_7^k} \right)}_{ij}} < \varphi \cdot {{\left( k \right)}^2}}\\ {0,{{\left( {{\bf{Gt}}_7^k} \right)}_{ij}} \ge \varphi \cdot {{\left( k \right)}^2}} \end{array}} \right. $

      其中,(Gr7)ijGr7; (${\bf {Gt}}_7^k $)ij${\bf {Gt}}_7^k $

      根据卷积结果Gt7k,用阈值法进行二分类,得到初步聚类结果Gr7,即洪灾区域的初步检测结果(见图19)。观察图19的初步聚类结果可以发现,反向聚类的方法虽然减少了洪灾区域相干斑噪声的影响,但是由于洪灾发生时地表坑洼处存在积水,初始聚类的结果还存在大量区块的虚警。在3.4节中,根据融合后SAR图像中的灾害前河流位置信息,对洪灾区域进行空间约束,进一步提升检测性能。

      图  19  初步聚类结果Gr7

      Figure 19.  Preliminary clustering results Gr7

    • H-FCM算法在对洪灾区域的空间约束上提出了一种参数自适应的空间约束方法,有效地筛除疑似洪灾区域,并提升检测性能。由于洪灾区域的水体满足RCC原理,定义洪灾区域主要范围为RtRt表示初步聚类结果Gr7的4连通区域面积最大的连通区域。

      由于真实洪灾区域在整个图像空间存在较高的稀疏度,即非洪灾区域处于零像素值区域,洪灾区域呈现为非零像素值的连通区域。此时,等间距增大空间约束范围会增加计算复杂度。引入约束矩阵${\bf{C}}{{\bf{S}}^{{\beta _p}}} $作为初步聚类结果的空间约束,当约束系数βp不断减少时,整体的空间约束范围按照e指数级不断增大。由于真实洪灾区域在全局上具备较高的稀疏度,约束系数βp按照等间隔减小的过程,对应约束范围呈e指数级的扩大。因此,式(21)的搜索规则符合洪灾区域的区块稀疏分布特点,既能保证获得适应待检测图像的约束系数,又不增加额外的运算量。

      约束系数βp的自适应选择方法,主要是依据两个因素:

      (1) 根据RCC原理,洪灾发生时,洪灾区域在遥感图像上满足4连通区域的连通性;

      (2) 远离洪灾前河流位置的低洼地区在洪灾发生时与洪灾区域相邻接的概率较低。

      $ {\beta _p} = {\beta _{p - 1}} - \frac{{{\rm{e}} - 1}}{{100{\rm{e}}}},{\beta _0} = 1 $

      其中,e为自然常数,$ p = 1,2,··· ,100$

      $ {\left( {{\bf{C}}{{\bf{S}}^{{\beta _p}}}} \right)_{ij}} = \left\{ \begin{aligned} & {1,\exp \left( {\frac{{ - {d_{ij}}}}{{{d_{{\rm{max}}}}}}} \right) > {\beta _p}}\\ & {0,\exp \left( {\frac{{ - {d_{ij}}}}{{{d_{{\rm{max}}}}}}} \right) \le {\beta _p}} \end{aligned} \right. $

      其中,$ {\left( {{\bf{C}}{{\bf{S}}^{{\beta _p}}}} \right)_{ij}} \in {\bf{C}}{{\bf{S}}^{{\beta _p}}}$; dij表示Ea中第i行,第j列像素点与灾害前河流位置(即标记为T1的像素点)的最小曼哈顿距离,dmax表示所有dij元素中的最大值。

      本文采用曼哈顿距离来表示像素点之间的空间距离,可以在保证检测性能的基础上,降低算法的运算量。洪灾区域的空间约束式$ \exp \left( {{{ - {d_{ij}}}}/{{{d_{{\rm{max}}}}}}} \right)$的取值范围为(0.367, 1](见图20)。并且当dij越小时,$ \exp \left( {{{ - {d_{ij}}}}/{{{d_{{\rm{max}}}}}}} \right)$的值越接近1,当dij越大时,$ \exp \left( {{{ - {d_{ij}}}}/{{{d_{{\rm{max}}}}}}} \right)$的值越接近0.367,确保距离灾害前河流位置越近的像素点被检测为洪灾区域的概率越大。

      图  20  洪灾区域空间约束曲线

      Figure 20.  Spatial constraint curve of flooded area

      将洪灾区域主要范围Rt加上空间约束${\bf{C}}{{\bf{S}}^{{\beta _p}}} $,可以进一步提升洪灾区域检测结果的准确率。由式(23)可以获得约束系数为βp时的洪灾区域检测结果$ {{\bf{Out}}^{{\beta _p}}}$

      $ {\bf{Ou}}{{\bf{t}}^{{\beta _p}}} = {\bf{Rt}} \cdot {\bf{C}}{{\bf{S}}^{{\beta _p}}} $

      其中,“A·B”表示矩阵A与矩阵B的点乘运算,即矩阵A和矩阵B对应点相乘。

      约束系数βp自适应选择的目的是判定疑似灾区S是否为真正灾区,以获得洪灾区域的最终检测结果。所以在βp的搜索过程中,若对应的疑似灾区的面积增长比率远小于聚类模型本身的稀疏度时,则判定疑似灾区为未受灾区域

      $ {{S}} = {\bf{Ou}}{{\bf{t}}^{{\beta _p}}} - {\bf{Ou}}{{\bf{t}}^{{\beta _{p - 1}}}} $

      其中,“–”表示矩阵对应点像素值作差

      $ {\rm{rati}}{{\rm{o}}_p} = \frac{{\displaystyle\sum\limits_{m = 1}^{{n_1}} \displaystyle\sum\limits_{q = 1}^{{n_2}} {{\left( { S} \right)}_{mq}}}}{{\displaystyle\sum\limits_{m = 1}^{{n_1}} \sum\limits_{q = 1}^{{n_2}} {{\left( {{\bf{Ou}}{{\bf{t}}^{{\beta _{p - 1}}}}} \right)}_{mq}}}} \hspace{5pt} $

      $ {\beta _p} = \left\{ \begin{aligned} & {{\beta _{\rm{a}}},\; {\rm{rati}}{{\rm{o}}_p} \approx {{{K}}_2}}\\ & {{\beta _{\rm{b}}},\; {\rm{rati}}{{\rm{o}}_p} < 0.01 \times {{{K}}_2}} \end{aligned} \right. $

      其中,“≈”近似的误差为±0.001; $ p = 1,2, ··· ,100$

      $ {\bf{FF}} = \sum\limits_{j = 1}^l \left[ {{{{S}}_j} \cdot {\rm{d}}{{\rm{t}}_j}} \right] + {\bf{Ou}}{{\bf{t}}^{{\beta _{\rm{a}}}}} $

      其中,Sj表示S中的第j个4连通区域,其背景大小为n1×n2

      $ {\rm{d}}{{\rm{t}}_j} = \left\{ \begin{aligned} & {1,\; {{\rm{\alpha }}_1} < \dfrac{{{A_j}}}{{{{{l}}_j}}} < {{\rm{\alpha }}_2}}\\ & {0,\; {\rm{otherwise}}} \end{aligned} \right. $

      其中,AjSj的面积,$ {{\rm{\alpha }}_1},{{\rm{\alpha }}_2}$为判定系数

      $ {{{l}}_j} = \sum\limits_{m = 1}^{{n_1}} \sum\limits_{q = 1}^{{n_2}} {\left( {{\bf{C}}{{\bf{H}}_j}} \right)_{mq}} $

      其中,lj表示重合边界长度;(CHj)mqCHj; CHjSj${{\bf{Out}}^{{\beta _{\rm a}}}} $的重合边界,其宽度为2个像素单元。

      将得到的约束系数βa, βb带入到式(24),获得待定的疑似灾区S。根据RCC原理,若疑似的4连通区域为真实洪灾区域,则其像素面积AjSj${{\bf{Out}}^{{\beta _{\rm a}}}} $邻接的边长lj的比值应处于$ \left( {{{\rm{\alpha }}_1},{{\rm{\alpha }}_2}} \right)$区间内。当面积边长比过大时,说明疑似灾区与真实洪灾区域相邻接的边界较短,可判定为非灾区;当面积边长比过小时,说明疑似灾区面积较小,属于区块噪声的可能性比较大。最后,根据每一个疑似区域的判定结果dtj确定该区域是否为洪灾区域,并将判定结果与${{\bf{Out}}^{{\beta _{\rm a}}}} $累加得到最终的洪灾区域检测结果FF

    • (1) 正确检测概率。Righta表示所有像素点被正确检测为变化的概率与被正确检测为未变化的概率之和

      $ {\rm{Righ}}{{\rm{t}}_{\rm{a}}} = {\rm{}}\frac{{{N_{{\rm{r}}\_{\rm{c}}}} + {N_{{\rm{r}}\_{\rm{uc}}}}}}{{{N_{\rm{t}}}}} $

      其中,Nr_c表示正确检测为变化的像素点总数,Nr_uc表示正确检测为未变化的像素点总数,Nt表示像素点总数。

      (2) 遗漏检测概率。Missa表示所有实际发生变化的像素点没有被检测为发生变化的概率

      $ {\rm{Mis}}{{\rm{s}}_{\rm{a}}} = {\rm{}}\frac{{{N_{\rm{m}}}}}{{{N_{\rm{c}}}}} $

      其中,Nc表示实际发生变化的像素点总数;Nm表示实际变化的像素点未被检测为变化的像素点总数。

      (3) 错误检测概率。Wra表示所有像素点被错误检测的概率,即实际未发生变化的像素点被错误检测为发生变化的概率与实际发生变化的像素点被检测为未变化的概率之和

      $ {\rm{W}}{{\rm{r}}_{\rm{a}}} = {\rm{}}\frac{{{N_{{\rm{f}}\_{\rm{u}}}}}}{{{N_{{\rm{un}}\_{\rm{c}}}}}} + \frac{{{N_{{\rm{f}}\_{\rm{c}}}}}}{{{N_{\rm{c}}}}} $

      其中,Nf_u表示实际发生变化的像素点被检测为未变化的总数;Nun_c表示实际未发生变化的像素点总数,Nf_c表示实际未发生变化像素点被检测为变化的总数。

      (4) 卡帕系数

      $ {\rm{Kappa}} = {\rm{}}\frac{{{p_{\rm{o}}} - {p_{\rm{e}}}}}{{1 - {p_{\rm{e}}}}} $

      其中,po表示总体分类精度。

      $ {p_{\rm{o}}} = {\rm{}}\frac{{{N_{{\rm{r}}\_{\rm{c}}}} + {N_{{\rm{r}}\_{\rm{uc}}}}}}{{{N_{\rm{t}}}}} \hspace{50pt} $

      $ {p_{\rm{e}}} = {\rm{}}\frac{{{N_{{\rm{d}}\_{\rm{c}}}} \cdot {N_{\rm{c}}} + {N_{{\rm{d}}\_{\rm{uc}}}} \cdot {N_{{\rm{r}}\_{\rm{uc}}}}}}{{{N_{\rm{t}}} \cdot {N_{\rm{t}}}}} $

      Nd_c表示检测为发生变化的像素点总数,Nd_uc表示检测为未变化的像素点总数。

      (5) 计算时间。以计算时间衡量算法的时间复杂度,计算时间定义为Time,其单位为“s”(秒)。

      (6) 迭代次数。以运算的迭代次数体现算法的运算复杂度。迭代次数定义为Iter,其单位为“次”,表示对应算法获得聚类中心的迭代次数。

    • 本文的实验数据是在确定目标灾害区域地理坐标的基础上,获得洪灾前遥感图像和洪灾后遥感图像的历史数据。实验数据的预处理过程主要分为4步。首先,利用ENVI5.3软件的Orthorectification模块对洪灾前后的遥感数据分别进行正射矫正;其次,利用ENVI5.3软件Registration模块中的Image to Image方式,以洪灾前的光学图像为基准,选择关键角点为控制点对同地区的其它遥感数据进行配准;再次,裁剪相同大小配准后遥感数据作为实验数据;最后,利用ENVI软件获得洪灾后光学图像的NDVI图像,通过人工标注的方法确定真实洪灾区域,作为各类算法检测结果性能指标评定的参照对象。

      H-FCM算法的参数选择,只需要确定判定系数α1, α2。卷积系数k、聚类系数$\varphi $可以通过分析聚类模型的稀疏度获得,约束系数βa, βb可以通过自适应选择获得。判定系数的范围是通过统计英国格洛斯特的历史洪灾数据的贝叶斯更新结果[30]获得的最优区间,在H-FCM算法中被设置为固定的区间。并且该区间定义在4.3节中实验1和实验2的结果中均获得了较好的检测性能。

      为了验证H-FCM算法在洪灾检测领域的普适性,本文通过收集案例中的河流历史洪灾信息,利用历史实测数据来确定不同河段发生洪灾的空间约束系数。基于该河流的历史灾害数据进行贝叶斯更新[30],可以发现某一河段附近地区受洪灾影响的范围与该河段的河流宽度、河床高度、河道转角系数密切相关。参考英国格洛斯特的历史灾害数据[31],可以得到洪灾区域的约束系数βa, βb与河流平均宽度lm的近似对应关系

      $ {\beta _{\rm{a}}} = {\rm{}}\exp \left( {\frac{{ - 100\sqrt 2 {l_{\rm{m}}}}}{{{d_{{\rm{max}}}}}}} \right){\rm{}} $

      $ {\beta _{\rm{b}}} = {\rm{}}\exp \left( {\frac{{ - 200\sqrt 3 {l_{\rm{m}}}}}{{{d_{{\rm{max}}}}}}} \right){\rm{}} $

      其中,lm通过计算图15的平均河宽得到。

      计算实际数据得到dmax=1338, lm=4.5,其单位均为像素单元。依据真实灾害数据统计的近似对应关系得到的βa=0.8997, βb=0.7719。对比自适应选择的约束系数,βa的选择结果与洪灾的统计历史数据较为吻合,但是βb的选择结果与洪灾的统计历史数据有一部分偏差。这是由于该流域的历史洪灾规模不同,本例所体现的洪灾区域要大于历史平均水平,所以本例的βb的自适应选择结果要小于根据历史统计得出的结果,对应实际更大的洪灾区域。

    • 为充分论证H-FCM算法对洪灾前光学图像和洪灾后SAR图像的变化检测性能,在对比实验中加入了图像分割算法—分水岭算法(Watershed Algorithm, WA)[32],边缘检测算法(Snake)[33], H-Kmeans算法等3种方法作为对比实验方法。此外,将1999年英国格洛斯特洪灾数据、2019年中国南昌附近洪灾的洪灾数据等两个场景作为实验数据。

    • 以1999年英国格洛斯特洪灾前的近红外波段图像和洪灾后的SAR图像[10]为例,检验H-FCM算法的有效性。将融合后的SAR图像(见图16)记作Fu,其大小为2359×1318像素,采用的聚类模型为Gm7。设定判定系数α1=120, α2=160,计算获得其卷积系数k=85,聚类系数$\varphi $=0.054。通过自适应选择得到的约束系数βa=0.9199, βb=0.7199,最终的检测结果为图21

      图  21  最终检测结果FF

      Figure 21.  Final experimental result FF

      首先,将实验1中洪灾后的NDVI图[10] (见图22)进行人工标注获得洪灾区域的真值图(见图23)。WA算法的洪灾区域检测方法是根据整体像素值分布(见图24),将所有像素点进行二值化,然后求出零值与最近非零值的距离,画出分水岭,进而对图像进行二分类。Snake算法获得Ea的边缘轮廓(见图25),按照式(18)和式(20)对边缘轮廓检测结果进行卷积运算得到最终洪灾区域检测结果。H-Kmeans算法是在Kmeans聚类(八分类)的基础上,进行分级聚类,从而得到洪灾区域的检测结果。

      图  22  洪灾后的NDVI图像[10]

      Figure 22.  NDVI images after the flood[10]

      图  23  人工标注的真值图

      Figure 23.  Manually labeled truth map

      图  24  像素值分布3维图

      Figure 24.  3-D map of the distribution of pixel values

      图  25  Ea的边缘轮廓

      Figure 25.  Edge profile of Ea

      其次,将H-FCM算法、WA算法、Snake算法以及H-Kmeans算法获得的洪灾区域检测结果与真值图(图23)比较,得到4种算法的检测结果与真实灾区的比较结果(见图26)。在图26中,白色部分表示实际洪灾区域被检测为洪灾区域;绿色部分表示实际洪灾区域被检测为非洪灾区域;红色区域表示实际非洪灾区域被检测为洪灾区域;黑色部分表示实际非洪灾区域被检测为非洪灾区域。本文后续实验的检测比较结果与图26的表示方法一致。

      图  26  4种算法的洪灾区域检测结果

      Figure 26.  Detection results of flooded area based on four algorithms

      根据4.1节中对6项评价指标的规定,按照4种算法分别用matlab进行计算分析,得到每一种算法检测结果对应的性能指标(见表3)。观察表3的各项评价指标可以发现H-FCM算法有最高的正确检测概率、最低的错误检测概率以及最高的Kappa系数。H-Kmean算法虽然有着与H-FCM算法接近一致的检测性能,但是其计算成本较高。综合考虑洪灾区域检测的高准确检测率和高时效性的要求,H-FCM算法具有相对最优的检测性能。

      方法H-FCMWASnakeH-Kmeans
      Righta0.89460.43460.67570.8837
      Missa0.15230.06800.16070.1519
      Wra0.25130.70010.50700.2634
      Kappa0.60920.18400.29640.5837
      Time(s)51.982.722.5373.27
      Iter(次)251125

      表 3  融合后SAR图像洪灾区域的检测结果

      Table 3.  Detection results of flooded area in fusioned SAR image

    • 实验2的数据来源于欧洲航天局的开源数据,分别为哨兵1号和哨兵2号卫星在2019年6月对中国南昌附近的观测数据(见图28)。按照H-FCM算法,首先,将图28(a)中提取的的洪灾前河流空间位置与图28(d)中洪灾后的SAR图像融合;其次,对融合图像进行分级聚类;最后,对聚类结果加以空间约束并获得检测结果。根据图28(c)的洪灾后的光学数据,人工标注真实洪灾区域范围(见图27)。比较4种检测方法对融合后SAR图像的洪灾区域检测结果(见图29)和性能指标(见表4),可以得出H-FCM算法具有相对最优的检测性能。

      图  27  人工标注的洪灾区域

      Figure 27.  Manually labeled flooded areas

      方法H-FCMWASnakeH-Kmeans
      Righta0.94140.62300.64260.9391
      Missa0.25200.05550.06080.2434
      Wra0.28860.46900.45180.2836
      Kappa0.69110.26290.27490.6848
      Time(s)30.592.823.1479.13
      Iter(次)251125

      表 4  融合后SAR图像洪灾区域的检测结果

      Table 4.  Detection results of flooded area in the fusioned SAR image

      图  28  洪灾前后哨兵1,哨兵2数据

      Figure 28.  Sentinel 1, 2 data before and after flood

      图  29  4种算法的洪灾区域检测结果

      Figure 29.  Detection results of flooded area based on four algorithms

      实验2的参数设定如下,将融合后的SAR图像记作NCa,其大小为1025×1025像素,采用的聚类模型为Gm7;设定判定系数α1=120, α2=160;卷积系数k,聚类系数$\varphi $和约束系数βa, βb均基于NCa由式(19)和式(26)获得。

    • Refined Lee filter(RL)[34]在滤除相干斑噪声的同时能较好地保持SAR图像的边缘特征。为分析图像降噪处理对上述算法检测性能的影响,首先对实验1中洪灾后的SAR图像进行RL滤波降噪,再利用上述4种方法对降噪处理后的图片进行洪灾区域检测。定义新的检测方法分别为RL_H-FCM算法、RL_WA算法、RL_Snake算法、RL_H-Kmeans算法,其检测结果如图30所示,检测性能评价指标如表5所示。

      方法RL_H-FCMRL_WARL_SnakeRL_H-Kmeans
      Righta0.89960.43910.76150.8899
      Missa0.17500.06600.25980.1351
      Wra0.26530.69330.49540.2418
      Kappa0.61470.18630.34600.6032
      Time(s)646.04593.72593.59745.62
      Iter(次)251125

      表 5  RL滤波后SAR图像的洪灾区域检测结果

      Table 5.  Detection results of flooded area in the SAR image after RL filtering

      图  30  RL滤波后4种算法的洪灾区域检测结果

      Figure 30.  Detection results of flooded area based on four algorithms after RL filtering

      比较表3表5中上述4种算法对应的性能指标,可以发现RL_Snake算法的检测性能提升最为明显。这是由于Snake算法是先获得图像的边缘特征,再对边缘轮廓的内外进行区分。通过RL滤波可以减少相干斑噪声的影响,使得图像整体像素值的均值和方差更稳定,有利于提取图像的边缘特征。但是通过RL滤波对H-FCM算法、WA算法和H-Kmeans算法的性能提升有限,反而极大地增加了算法的计算复杂度。总体而言,利用H-FCM算法直接对洪灾前的光学图像和洪灾后的SAR图像进行洪灾区域检测时,其检测性能与降噪处理后的图像检测性能几乎一致,而直接检测的运算时间大幅度减少。所以应用H-FCM算法直接处理配准后的光学和SAR图像,既可以缩短洪灾区域的检测时间又可以保证相对较高的检测准确率。

    • 本文提出的H-FCM算法综合考虑洪灾区域检测的高准确率和高时效性的需求,既提升了SAR图像的洪灾区域检测精度又缩短了检测时间。H-FCM算法区别于传统的无监督聚类算法,没有直接对待处理的图像进行无监督聚类,而是提出了一种反向聚类的方法。首先,H-FCM算法将洪灾前光学图像的河流位置检测结果与洪灾后的SAR图像融合;其次,在融合图像的最高级像素值区间对零像素值点进行分级聚类;最后,H-FCM算法在对融合图像的洪灾区域进行分级聚类的基础上,将洪灾前的河流位置作为空间约束,进一步筛除疑似洪灾区域以提升检测性能。

      H-FCM算法将洪灾前的光学图像和洪灾后的SAR图像融合后进行洪灾区域检测,实现了全天侯、高准确率、高时效的洪灾区域检测要求。虽然H-FCM算法在一定程度上解决了SAR图像相干斑噪声给洪灾区域检测造成的问题,但是还有部分灾区的细节没有精确检测到,这有待在未来的工作中解决。

参考文献 (34)

目录

    /

    返回文章
    返回