CN116704369A - 一种面向对象的光学与sar遥感影像融合洪水提取方法和系统 - Google Patents
一种面向对象的光学与sar遥感影像融合洪水提取方法和系统 Download PDFInfo
- Publication number
- CN116704369A CN116704369A CN202310671107.8A CN202310671107A CN116704369A CN 116704369 A CN116704369 A CN 116704369A CN 202310671107 A CN202310671107 A CN 202310671107A CN 116704369 A CN116704369 A CN 116704369A
- Authority
- CN
- China
- Prior art keywords
- image
- water body
- sar
- optical
- remote sensing
- 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.)
- Withdrawn
Links
- 230000003287 optical effect Effects 0.000 title claims abstract description 111
- 238000000605 extraction Methods 0.000 title claims abstract description 63
- 230000004927 fusion Effects 0.000 title claims abstract description 21
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 265
- 238000010586 diagram Methods 0.000 claims abstract description 90
- 230000011218 segmentation Effects 0.000 claims abstract description 48
- 238000001914 filtration Methods 0.000 claims abstract description 26
- 238000003066 decision tree Methods 0.000 claims abstract description 18
- 238000000034 method Methods 0.000 claims description 36
- 238000012545 processing Methods 0.000 claims description 22
- 238000007781 pre-processing Methods 0.000 claims description 16
- 238000012937 correction Methods 0.000 claims description 12
- 238000003709 image segmentation Methods 0.000 claims description 11
- 238000002310 reflectometry Methods 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 6
- 239000000284 extract Substances 0.000 claims description 6
- 230000010287 polarization Effects 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 5
- 230000000877 morphologic effect Effects 0.000 claims description 5
- 230000005855 radiation Effects 0.000 claims description 5
- 238000012544 monitoring process Methods 0.000 description 7
- 238000011161 development Methods 0.000 description 4
- 230000018109 developmental process Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000003809 water extraction Methods 0.000 description 4
- 230000002902 bimodal effect Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 241000132092 Aster Species 0.000 description 1
- 238000012300 Sequence Analysis Methods 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012634 optical imaging Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
- G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
- G06V10/7715—Feature extraction, e.g. by transforming the feature space, e.g. multi-dimensional scaling [MDS]; Mappings, e.g. subspace methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
- G06V10/80—Fusion, i.e. combining data from various sources at the sensor level, preprocessing level, feature extraction level or classification level
- G06V10/806—Fusion, i.e. combining data from various sources at the sensor level, preprocessing level, feature extraction level or classification level of extracted features
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Multimedia (AREA)
- Evolutionary Computation (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Remote Sensing (AREA)
- Astronomy & Astrophysics (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种面向对象的光学与SAR遥感影像融合洪水提取方法和系统，提出了利用规则决策树提取光学影像水体图、利用光学水体图辅助OTSU提取SAR影像水体图、以及利用时空滤波对水体图进行去噪的光学与SAR影像融合洪水提取方法，重点解决当前方法中可用遥感影像数量少、洪水图时间分辨率低、SAR影像水体提取精度低、洪水图受噪声影响严重的问题。采用mNDWI和NDVI构建规则决策树实现光学图像水体的高精度提取；采用光学图像辅助的OTSU算法提取SAR图像的水体，解决了传统OTSU算法无法准确定位最优分割阈值的难题，明显提高了SAR图像水体提取的精度；考虑洪水演变特点，构建时空滤波器去除时序水体图中的噪声和奇异值，提高了水体图的精度与可靠性。
Description
技术领域
本发明属于遥感影像处理领域，特别是涉及一种面向对象的光学与SAR遥感影像融合洪水提取方法和系统。
背景技术
利用遥感卫星观测数据对洪水事件进行动态监测是保护地区气候、安全和环境的重要手段。多光谱遥感数据(如中分辨率成像光谱仪MODIS、陆地卫星Landsat数据、Sentinel-2数据)被广泛用于洪水探测。目前已经开发了许多利用多光谱遥感影像进行洪水检测的方法，包括基于阈值的方法、基于指数的方法、光谱角距离法、基于分类的方法等。水体指数和植被指数，如归一化差水指数(NDWI)、修正NDWI(mNDWI)、归一化差植被指数(NDVI)和增强型植被指数(EVI)被广泛用于多光谱光学图像的水体提取。洪水事件往往伴随着多云、多风和多雨等天气，不利的天气条件阻碍了电磁波在光学波谱范围内的传播，从而阻碍了光学传感器的数据采集。考虑到光学成像在洪水监测中的局限性，可以全天候观测的SAR传感器是一个有力的补充甚至替代手段。基于SAR影像的洪水提取大多是通过设置HH、HV、VV或者VH极化波段的后向散射系数阈值，利用阈值对SAR图像进行分割从而提取出水体。
现有技术的缺陷和不足：
1.现有洪水提取方法主要依赖于单传感器影像数据。由于受到卫星重访周期、影像空间分辨率、云污染等因素的影响，单传感器数据对洪水事件的重访频率低、可利用特征单一，无法提供连续准确的洪水发展信息，从而造成可以遥感影像数量不足、专题图时间分辨率低的问题。
2.现有洪水提取方法大多采用基于像素的分析方法。基于像素的方法以像素为单位进行特征提取与分析，仅仅利用了像素尺度的特征(如光谱特征)，较少考虑像素之间的空间关系，或者对空间关系的利用有限；容易受到噪声影响，特别不利于SAR图像的处理，造成水体图精度较低、地物边缘不规则、与实际边缘贴合不准确且存在大量斑点噪声。
3.现有洪水提取方法算法复杂度较高、鲁棒性不足、可移植性差。
发明内容
针对现有算法可用影像数量少、专题图时间分辨率低、SAR水体图精度低等问题，本发明提供了基于规则决策树光学影像水体提取算法、光学辅助SAR影像OTSU水体提取算法和水体图时空滤波去噪的光学与SAR影像融合洪水提取方法，成功提取了高精度、高时间分辨率的洪水图，有效揭示了洪水演变过程。
为了达到上述目的，本发明提供的技术方案是：一种面向对象的光学与SAR遥感影像融合洪水提取方法，包括如下步骤：
步骤1，对SAR遥感图像进行预处理、图像分割与特征提取，得到对象级的遥感图像和用于水体提取的遥感极化后向散射系数特征；对光学图像进行图像分割与特征提取，得到用于水体提取的光学遥感指数特征；
步骤2，利用遥感指数特征构建规则决策树，从光学图像中提取高精度的水体图；
步骤3，结合VV极化波段的后向散射系数特征，利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图；
步骤4，对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图；
步骤5，利用时序水体图进行洪水演变时间动态分析。
进一步的，步骤1中，所述预处理包括对遥感图像进行辐射校正、大气校正、正射校正，以及交叉校准；图像分割包括：采用SNIC超像素分割算法，对光学与SAR图像分别进行超像素分割；特征提取包括：从光学影像上分别提取修正归一化差值水体指数mNDWI、归一化差植被指数NDVI作为特征，用于构建规则决策树提取光学影像上的水体，从SAR图像上选取VV极化波段的后向散射系数作为特征，用于OTSU阈值法提取SAR图像上的水体。
进一步的，步骤2的具体实现方式如下；
(21)采用式(1)和(2)分别计算光学图像上每个图像单元的mNDWI和NDVI指数，得到两种指数图像；
其中，ρGreen、ρSWIR、ρNIR、ρRed分别代表光学影像上绿波段、短驳红外波段、近红外波段和红波段的地表反射率；
(22)采用决策准则mNDWI>NDVI提取水体单元；
(23)采用决策准则mNDWI>0去除上一步结果中的非水体单元；
(24)采用式(4)进行图像二值化处理得到水体地图；
其中，G(i)代表遥感影像上图像单元i处的元素值。
进一步的，步骤3的具体实现方式如下；
(31)计算初始分割阈值；计算SAR图像VV极化波段的灰度直方图，采用OTSU算法自动确定最优分割阈值，把这一阈值定义为初始分割阈值，用变量T0表示；
(32)设置阈值搜索区间；在灰度直方图上选取初始分割阈值T0左、右长度为ε的BC值域范围[T0-ε,T0+ε]，将其设置为搜索区间；
(33)依次取阈值提取SAR图像的水体；以设定的步长在搜索区间内依次取值T(j)，以T(j)为分割阈值按照式(5)提取SAR图像的水体，采用式(4)二值化得到水体图像，并统计水体图像中水体单元的个数；
其中i表示SAR图像的基本处理单元，G(i)表示遥感影像上图像单元的元素值；公式的含义是将后向散射系数值小于阈值的图像单元划分为水体，否则划分为非水体；
(34)循环步骤(33)，直到遍历完整个搜索区间；
(35)光学图像辅助确定最优分割阈值；计算SAR水体结果图中水体单元个数与同期光学遥感图像提取的水体图像中水体单元个数的差值，取差值最小的T(j)为SAR图像水体提取的最优阈值T；
(36)生成最终水体图；利用最优阈值T提取SAR图像上的水体，通过形态滤波去除散斑噪声因素引起的伪水像元，并按式(4)进行二值化处理，得到最终水体图。
进一步的，步骤4的具体实现方式如下；
(41)利用生成的时序水体图像构建三维数据立方体；时序水体图像构成的数据立方体用一个三维矩阵来表示，其中，任一点i可以用一组三维坐标i(x,y,t)来表示，x、y代表像素的二维空间坐标，t代表最小图像单元的时间维坐标，即其所在的水体图像的时间序号；i处的元素值G(i)为1表示水体单元，为0表示非水体单元；
(42)设置空间邻域与时间邻域；以图像元素i为中心点，取周围3×3大小的窗口作为空间邻域，前后长度为3的窗口作为时间邻域，空间邻域和时间邻域共同构成中心元素点i的单元立方体，作为时空滤波的最小处理单元，单元立方体内总共包含3×3×3＝27个元素，中心元素i(x,y,t)的邻域表示为：
空间邻域：时间邻域：[t-1,t,t+1]
(43)时空滤波处理；在以i为中心元素的单元立方体内，按照决策规则对图像值G(i)进行修正，得到修正后的G(i)值；
(44)重复(42)和(43)步骤，直到处理完所有元素点i(x,y,t)，输出最终时序水体图。
进一步的，步骤(43)中的决策规则如下；
当单元立方体内中心元素i在当前时相的元素值G(i)为1时，如果该单元立方体内值为1的元素个数之和大于18，或者中心元素在前后两个时相的G(i)均为1时，则判定当前中心元素i为水体单元，否则判定为非水体单元；当单元立方体内中心元素i在当前时相的元素值G(i)为0时，如果单元立方体内值为1的元素个数之和小于3，则判定中心元素为非水体单元，否则判定为水体单元。
进一步的，步骤5的具体实现方式如下；
从每幅水体图上计算洪水覆盖面积，得到洪水淹没面积时间变化曲线；融合多源遥感影洪水提取结果，得到融合洪水地图，对于每个像素，计算时序二进制洪水图中像素值1出现的数量，并记录出现的日期，最终将它们融合成一个洪水地图，利用该洪地图分析洪水淹没地点、持续时间，评估灾害破坏程度。
本发明还提供一种面向对象的光学与SAR遥感影像融合洪水提取系统，包括如下模块：
预处理模块，对SAR遥感图像进行预处理、图像分割与特征提取，得到对象级的遥感图像和用于水体提取的遥感极化后向散射系数特征；对光学图像进行图像分割与特征提取，得到用于水体提取的光学遥感指数特征；
光学水体图获取模块，用于利用遥感指数特征构建规则决策树，从光学图像中提取高精度的水体图；
SAR水体图获取模块，用于结合VV极化波段的后向散射系数特征，利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图；
时序水体图获取模块，用于对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图；
分析模块，用于利用时序水体图进行洪水演变时间动态分析。
进一步的，光学水体图获取模块的具体实现方式如下；
(21)采用式(1)和(2)分别计算光学图像上每个图像单元的mNDWI和NDVI指数，得到两种指数图像；
其中，ρGreen、ρSWIR、ρNIR、ρRed分别代表光学影像上绿波段、短驳红外波段、近红外波段和红波段的地表反射率；
(22)采用决策准则mNDWI>NDVI提取水体单元；
(23)采用决策准则mNDWI>0去除上一步结果中的非水体单元；
(24)采用式(4)进行图像二值化处理得到水体地图；
其中，G(i)代表遥感影像上图像单元i处的元素值。
进一步的，SAR水体图获取模块的具体实现方式如下；
(31)计算初始分割阈值；计算SAR图像VV极化波段的灰度直方图，采用OTSU算法自动确定最优分割阈值，把这一阈值定义为初始分割阈值，用变量T0表示；
(32)设置阈值搜索区间；在灰度直方图上选取初始分割阈值T0左、右长度为ε的BC值域范围[T0-ε,T0+ε]，将其设置为搜索区间；
(33)依次取阈值提取SAR图像的水体；以设定的步长在搜索区间内依次取值T(j)，以T(j)为分割阈值按照式(5)提取SAR图像的水体，采用式(4)二值化得到水体图像，并统计水体图像中水体单元的个数；
其中i表示SAR图像的基本处理单元，G(i)表示遥感影像上图像单元的元素值；公式的含义是将后向散射系数值小于阈值的图像单元划分为水体，否则划分为非水体；
(34)循环步骤(33)，直到遍历完整个搜索区间；
(35)光学图像辅助确定最优分割阈值；计算SAR水体结果图中水体单元个数与同期光学遥感图像提取的水体图像中水体单元个数的差值，取差值最小的T(j)为SAR图像水体提取的最优阈值T；
(36)生成最终水体图；利用最优阈值T提取SAR图像上的水体，通过形态滤波去除散斑噪声因素引起的伪水像元，并按式(4)进行二值化处理，得到最终水体图。
本发明提出了利用规则决策树提取光学影像水体图、利用光学水体图辅助OTSU提取SAR影像水体图、以及利用时空滤波对水体图进行去噪的光学与SAR影像融合洪水提取方法，重点解决当前方法中可用遥感影像数量少、洪水图时间分辨率低、SAR影像水体提取精度低、洪水图受噪声影响严重的问题。采用mNDWI和NDVI构建规则决策树实现光学图像水体的高精度提取；采用光学图像辅助的OTSU算法提取SAR图像的水体，解决了传统OTSU算法无法准确定位最优分割阈值的难题，明显提高了SAR图像水体提取的精度；考虑洪水演变特点，构建时空滤波器去除时序水体图中的噪声和奇异值，提高了水体图的精度与可靠性。
采用以上技术方案与现有技术相比，具有以下有益效果：
(1)通过融合光学与SAR遥感影像，解决洪水监测过程中由于传感器重访周期长、光学图像获取困难等原因导致的图像数量少、图像质量低、专题图时间分辨率低的问题。
(2)采用光学图像辅助的OTSU算法提取SAR图像的水体，解决了传统OTSU算法无法准确定位最优分割阈值的难题，明显提高了SAR图像水体提取的精度。
(3)针对洪水演变特点，设计了一种时空滤波算法对时序水体图进行后处理，去除了洪水序列图中的噪声和奇异值，提高了洪水图的精度与可靠性。
(4)该方法采用经典且简单的算法及其组合，在获取高精度的洪水图、有效监测洪水动态的同时，具有算法复杂度低、鲁棒性好、可移植性高等优势。
附图说明
图1为本发明的技术流程图。
图2为光学与SAR图像超像素分割结果。
图3为基于规则决策树的光学图像水体提取流程。
图4为洪水期间获取的SAR图像的灰度直方图。
图5为光学辅助的SAR图像OTSU水体提取原理。
图6为时空滤波算法原理与流程。
图7为时空滤波前后的SAR水体图对比。
图8为洪水监测多源遥感时间序列数据集(彩色合成方式为真彩色合成)。
图9为本文方法提取的时序洪水图。
图10为洪水淹没面积时间变化曲线。
图11为多源遥感图像提取的融合洪水地图。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步说明。
本发明提供一种洪水提取方法，用于获取准确连续的洪水发展信息，如图1所示，包括如下步骤1至步骤5。步骤1，遥感图像预处理、分割与特征提取，得到对象级的遥感图像和用于水体提取的特征。步骤2，利用遥感指数特征构建规则决策树，从光学图像中提取高精度的水体图。步骤3，利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图。步骤4，对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图。步骤5，利用时序水体图进行洪水演变时间动态分析。下面对各步骤的具体实现进行详细说明：
1.遥感图像预处理、分割与特征提取。获取覆盖研究区域的Sentinel-1SAR、Landsat-8、Sentinel-2多光谱遥感影像，并对所有影像数据进行预处理、图像分割与特征提取。
预处理：本次使用的是Landsat-8地表反射率图像，在图像预处理过程中，进行了辐射校正、大气校正、正射校正，以及Landsat不同系列传感器之间的交叉校准，最终重采样到20米分辨率。Sentinel-1地距图像(GRD)是后向散射系数(BC)图像集合，使用Sentinel-1工具箱进行了热噪声去除、散斑噪声去除和辐射校准等图像预处理；并使用SRTM DEM进行了地形校正(对纬度大于60°且SRTM不可用的地区，采用ASTER DEM校正)；接着通过对数缩放(即10*log10(x))公式将像素值转换为以分贝单位的BC值；最后使用距离-多普勒算法对BC图像进行精确的地理编码，并重新采样到20m×20m分辨率。本次使用的Sentinel-2数据是正射地图投影的2A级地表反射率数据(Bottom of Atmosphere Reflectance,BOA)，是采用sen2cor工具箱由1C级数据经过辐射校准和大气校正生成的，最终重采样到20米分辨率。将所有遥感影像按照时间顺序排列，得到研究区时间序列遥感影像数据集，如图8所示。本研究使用的所有光学图像均未经过去云和插值处理，因为受到云等污染无法使用的图像都被舍弃，只保留不受污染影响的高质量图像，这么做可以有效降低去云算法和插值处理等操作引入的数据误差，保证结果精度和时间序列分析的可靠性。
图像分割：采用SNIC超像素分割算法，对光学与SAR图像分别进行超像素分割。其中，Sentinel-2(S1)和Landsat-8(L8)多光谱图像分割采用的特征是波段+指数，其中波段选择适用于各种地类识别的蓝、绿、红、近红外、短波红外1、短波红外2，指数选用对水体敏感的mNDWI以及对地表覆盖敏感的NDVI。Sentinel-1SAR(S1)图像分割使用VH、VV两个极化波段。超像素分割效果如图2所示，结果显示，超像素分割后，同一地物内部均质性提高，不同地物之间异质性增强，SAR图像的斑点噪声得到有效抑制，有利于后续图像处理与解译过程。
特征提取：分别采用相应的指数计算公式，从Landsat-8和Sentinel-2光学影像上分别提取修正归一化差值水体指数(mNDWI)、归一化差植被指数(NDVI)作为特征，用于构建规则决策树提取光学影像上的水体。从Sentinel-1SAR图像上选取VV极化波段的后向散射系数作为特征，用于阈值法提取SAR图像上的水体。
2.利用遥感指数特征mNDWI和NDVI构建规则决策树，从光学图像中提取高精度的水体图。提出了规则决策树法光学图像水体提取算法，结合对水体信号敏感的mNDWI指数和对植被信号敏感的NDVI指数构建规则决策树，提取光学遥感图像上的水体单元。算法流程如图3所示，具体实现过程如下：
(1)采用式(1)和(2)分别计算Landsat-8和Sentinel-2地表反射率图像上每个图像单元的mNDWI和NDVI指数，得到两种指数图像；
其中，ρGreen、ρSWIR、ρNIR、ρRed分别代表光学影像上绿波段、短驳红外波段、近红外波段和红波段的地表反射率。
(2)采用决策准则mNDWI>NDVI提取水体单元；
(3)采用决策准则mNDWI>0去除上一步结果中的非水体单元；
(4)采用式(4)进行图像二值化处理得到水体地图。
其中，G(i)代表遥感影像上图像单元i处的元素值。
3.利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图。
图4显示了洪水期间获取的SAR图像的后向散射系数(BackscatteringCoefficient,BC)分布直方图。对于频率分布直方图中有两个峰的SAR图像(即双峰图像)，OTSU算法得到的最佳阈值近似等于两个峰之间的谷值。然而，当SAR图像双峰特征不明显时，算法精度受到限制。提出一种基于OTSU的光学辅助SAR图像水体提取算法，在光学图像水体提取结果的辅助下，借助OTSU算法自动确定分割阈值实现SAR图像的水体提取，以提高SAR图像水体提取的精度。基于OTSU的光学辅助SAR图像水体提取算法的具体实现过程如下：
(1)计算初始分割阈值。计算Sentinel-1SAR图像VV极化波段的灰度直方图，采用OTSU算法自动确定最优分割阈值，把这一阈值定义为初始分割阈值，用变量T0表示；
(2)设置阈值搜索区间。在灰度直方图上选取初始分割阈值T0左、右长度为ε的BC值域范围[T0-ε,T0+ε]，将其设置为搜索区间；
(3)依次取阈值提取SAR图像的水体。以设定的步长(如0.1)在搜索区间内依次取值T(j)，以T(j)为分割阈值按照式(5)提取SAR图像的水体单元，采用式(4)二值化得到水体图像，并统计水体图像中水体单元的个数；
其中i表示SAR影像的图像单元，G(i)表示遥感影像上图像单元的元素值，公式的含义是将SAR图像上后向散射系数值小于阈值T(j)的图像单元划分为水体，否则划分为非水体。
(4)循环步骤(3)，直到遍历完整个搜索区间；
(5)光学图像辅助确定最优分割阈值。计算SAR水体结果图中水体单元个数与同期光学遥感图像提取的水体图像中水体单元个数的差值，取差值最小的T(j)为SAR图像水体提取的最优阈值T；
(6)生成最终水体图。利用最优阈值T提取SAR图像上的水体，通过形态滤波去除散斑噪声等因素引起的伪水像元，并按式(4)进行二值化处理，得到最终水体图。
图5描述了采用光学辅助OTSU进行SAR图像水体提取的两个实例，其中，(a)、(e)为原始Sentinel-1SAR图像；(b)、(f)为同期的Sentinel-2光学图像水体提取结果；(d)、(h)为光学辅助OTSU确定的SAR图像最优水体分割阈值，分别为-16.1和-16.5，表示BC值，单位为db；(c)、(g)为应用最优分割阈值从SAR图像提取的水体图。
4.对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图。设计了一种时空滤波算法，综合利用中心像素的空间邻域和时间邻域信息，对像素时间轨迹中的奇异值进行修正，提高时序洪水图的精度和可靠性。算法的原理与流程如图6所示，具体实现过程如下：
(1)利用生成的时序水体图像构建三维数据立方体。时序水体图像构成的数据立方体可以用一个三维矩阵来表示，其中，任一点i可以用一组三维坐标i(x,y,t)来表示，x、y代表像素的二维空间坐标，t代表最小图像单元的时间维坐标，即其所在的水体图像的时间序号。由于水体图像经过了二值化，因此根据式(4)，i处的元素值G(i)为1(水体单元)或者0(非水体单元)。构建三维矩阵后，就可以采用矩阵操作对数据进行处理。
(2)设置空间邻域与时间邻域。以图像元素i为中心点，取周围3×3大小的窗口作为空间邻域，前后长度为3的窗口作为时间邻域，空间邻域和时间邻域共同构成中心元素点i的单元立方体，作为时空滤波的最小处理单元，单元立方体内总共包含3×3×3＝27个元素。
中心元素i(x,y,t)的邻域可以表示为：
空间邻域：时间邻域：[t-1,t,t+1]
(3)时空滤波处理。在以i为中心元素的单元立方体内，按照图6所示的决策规则对图像值G(i)进行修正，得到修正后的G(i)值。具体来说：当单元立方体内中心元素i在当前时相的元素值G(i)为1时，如果该单元立方体内值为1的元素个数之和大于18，或者中心元素在前后两个时相的G(i)均为1时，则判定当前中心元素i为水体单元，否则判定为非水体单元；当单元立方体内中心元素i在当前时相的元素值G(i)为0时，如果单元立方体内值为1的元素个数之和小于3，则判定中心元素为非水体单元，否则判定为水体单元。
(4)重复(2)和(3)步骤，直到处理完所有元素点i(x,y,t)，输出最终时序水体图。图6展示了时空滤波算法对三景水体图的处理效果，可以看到，滤波后图像的噪声和伪水体单元明显减少，得到如图9所示的最终时序水体图。
5.利用时序水体图进行洪水演变时间动态分析。从每幅水体图上计算洪水覆盖面积，得到洪水淹没面积时间变化曲线，如图10所示，可以清晰观察到洪水淹没面积的变化趋势。融合多源遥感影洪水提取结果，可以得到融合洪水地图，对于每个像素，我们计算时序二进制洪水图中像素值1出现的数量，并记录出现的日期，最终将它们融合成一个洪水地图，如图11所示。利用该地图可以分析洪水淹没地点、持续时间，评估灾害破坏程度。多传感器图像组成的时间序列数据为高时间分辨率的洪水事件跟踪提供了有效工具。通过绘制SAR和光学图像融合的时间序列洪水图，可以清晰分析洪水的发展情况。地图上的不同颜色表示洪水淹没的持续时间，通过这个信息可以了解洪水的破坏程度。
上述技术方案所设计的洪水提取与监测方法，通过SAR和光学数据的结合，相当于增加了遥感传感器对研究区的重访频率，从而增加了洪水期间可用遥感观测图像的数量，受到云污染的光学图像以及低质量的SAR图像可以被弃用，间接提高了研究数据的质量。高时间分辨率的时序洪水图提供了洪水发展和淹没程度的详细信息，多源遥感图像提取的融合洪水地图显示了洪水淹没的持续时间以及淹没范围，显示了多传感器图像综合利用在洪水监测中的优势。
本发明还提供一种面向对象的光学与SAR遥感影像融合洪水提取系统，包括如下模块：
预处理模块，对SAR遥感图像进行预处理、图像分割与特征提取，得到对象级的遥感图像和用于水体提取的遥感极化后向散射系数特征；对光学图像进行图像分割与特征提取，得到用于水体提取的光学遥感指数特征；
光学水体图获取模块，用于利用遥感指数特征构建规则决策树，从光学图像中提取高精度的水体图；
SAR水体图获取模块，用于结合VV极化后向散射系数特征，利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图；
时序水体图获取模块，用于对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图；
分析模块，用于利用时序水体图进行洪水演变时间动态分析。
各模块的具体实现方法和各步骤相同，本发明不予撰述。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
Claims (10)
1.一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于，包括如下步骤：
步骤1，对SAR遥感图像进行预处理、图像分割与特征提取，得到对象级的遥感图像和用于水体提取的遥感极化后向散射系数特征；对光学图像进行预处理、图像分割与特征提取，得到用于水体提取的遥感指数特征；
步骤2，利用遥感指数特征构建规则决策树，从光学图像中提取高精度的水体图；
步骤3，结合VV极化波段的后向散射系数特征，利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图；
步骤4，对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图；
步骤5，利用时序水体图进行洪水演变时间动态分析。
2.如权利要求1所述的一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于：步骤1中，所述预处理包括对遥感图像进行辐射校正、大气校正、正射校正，以及交叉校准；图像分割包括：采用SNIC超像素分割算法，对光学与SAR图像分别进行超像素分割；特征提取包括：从光学影像上分别提取修正归一化差值水体指数mNDWI、归一化差植被指数NDVI作为特征，用于构建规则决策树提取光学影像上的水体，从SAR图像上选取VV极化波段的后向散射系数作为特征，用于OTSU阈值法提取SAR图像上的水体。
3.如权利要求1所述的一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于：步骤2的具体实现方式如下；
(21)采用式(1)和(2)分别计算光学图像上每个图像单元的mNDWI和NDVI指数，得到两种指数图像；
其中，ρGreen、ρSWIR、ρNIR、ρRed分别代表光学影像上绿波段、短驳红外波段、近红外波段和红波段的地表反射率；
(22)采用决策准则mNDWI>NDVI提取水体单元；
(23)采用决策准则mNDWI>0去除上一步结果中的非水体单元；
(24)采用式(4)进行图像二值化处理得到水体地图；
其中，G(i)代表遥感影像上图像单元i处的元素值。
4.如权利要求3所述的一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于：步骤3的具体实现方式如下；
(31)计算初始分割阈值；计算SAR图像VV极化波段的灰度直方图，采用OTSU算法自动确定最优分割阈值，把这一阈值定义为初始分割阈值，用变量T0表示；
(32)设置阈值搜索区间；在灰度直方图上选取初始分割阈值T0左、右长度为ε的BC值域范围[T0-ε,T0+ε]，将其设置为搜索区间；
(33)依次取阈值提取SAR图像的水体；以设定的步长在搜索区间内依次取值T(j)，以T(j)为分割阈值按照式(5)提取SAR图像的水体，采用式(4)二值化得到水体图像，并统计水体图像中水体单元的个数；
其中i表示SAR图像的基本处理单元，G(i)表示遥感影像上图像单元的元素值；公式的含义是将后向散射系数值小于阈值的图像单元划分为水体，否则划分为非水体；
(34)循环步骤(33)，直到遍历完整个搜索区间；
(35)光学图像辅助确定最优分割阈值；计算SAR水体结果图中水体单元个数与同期光学遥感图像提取的水体图像中水体单元个数的差值，取差值最小的T(j)为SAR图像水体提取的最优阈值T；
(36)生成最终水体图；利用最优阈值T提取SAR图像上的水体，通过形态滤波去除散斑噪声因素引起的伪水像元，并按式(4)进行二值化处理，得到最终水体图。
5.如权利要求1所述的一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于：步骤4的具体实现方式如下；
(41)利用生成的时序水体图像构建三维数据立方体；时序水体图像构成的数据立方体用一个三维矩阵来表示，其中，任一点i可以用一组三维坐标i(x,y,t)来表示，x、y代表像素的二维空间坐标，t代表最小图像单元的时间维坐标，即其所在的水体图像的时间序号；i处的元素值G(i)为1表示水体单元，为0表示非水体单元；
(42)设置空间邻域与时间邻域；以图像元素i为中心点，取周围3×3大小的窗口作为空间邻域，前后长度为3的窗口作为时间邻域，空间邻域和时间邻域共同构成中心元素点i的单元立方体，作为时空滤波的最小处理单元，单元立方体内总共包含3×3×3＝27个元素，中心元素i(x,y,t)的邻域表示为：
空间邻域：时间邻域：[t-1,t,t+1]
(43)时空滤波处理；在以i为中心元素的单元立方体内，按照决策规则对图像值G(i)进行修正，得到修正后的G(i)值；
(44)重复(42)和(43)步骤，直到处理完所有元素点i(x,y,t)，输出最终时序水体图。
6.如权利要求5所述的一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于：步骤(43)中的决策规则如下；
当单元立方体内中心元素i在当前时相的元素值G(i)为1时，如果该单元立方体内值为1的元素个数之和大于18，或者中心元素在前后两个时相的G(i)均为1时，则判定当前中心元素i为水体单元，否则判定为非水体单元；当单元立方体内中心元素i在当前时相的元素值G(i)为0时，如果单元立方体内值为1的元素个数之和小于3，则判定中心元素为非水体单元，否则判定为水体单元。
7.如权利要求1所述的一种面向对象的光学与SAR遥感影像融合洪水提取方法，其特征在于：步骤5的具体实现方式如下；
从每幅水体图上计算洪水覆盖面积，得到洪水淹没面积时间变化曲线；融合多源遥感影洪水提取结果，得到融合洪水地图，对于每个像素，计算时序二进制洪水图中像素值1出现的数量，并记录出现的日期，最终将它们融合成一个洪水地图，利用该洪地图分析洪水淹没地点、持续时间，评估灾害破坏程度。
8.一种面向对象的光学与SAR遥感影像融合洪水提取系统，其特征在于，包括如下模块：
预处理模块，对SAR遥感图像进行预处理、图像分割与特征提取，得到对象级的遥感图像和用于水体提取的遥感极化散射系数特征；对光学图像进行图像分割与特征提取，得到用于水体提取的光学遥感指数特征；
光学水体图获取模块，用于利用遥感指数特征构建规则决策树，从光学图像中提取高精度的水体图；
SAR水体图获取模块，用于结合VV极化波段的后向散射系数特征，利用从光学图像提取的水体图辅助SAR图像进行OTSU水体自动提取，获取高精度的水体图；
时序水体图获取模块，用于对光学与SAR影像提取的所有水体图构成的时间序列进行时空滤波，去除水体图中的噪声和奇异值，得到最终的时序水体图；
分析模块，用于利用时序水体图进行洪水演变时间动态分析。
9.如权利要求8所述的一种面向对象的光学与SAR遥感影像融合洪水提取系统，其特征在于：光学水体图获取模块的具体实现方式如下；
(21)采用式(1)和(2)分别计算光学图像上每个图像单元的mNDWI和NDVI指数，得到两种指数图像；
其中，ρGreen、ρSWIR、ρNIR、ρRed分别代表光学影像上绿波段、短驳红外波段、近红外波段和红波段的地表反射率；
(22)采用决策准则mNDWI>NDVI提取水体单元；
(23)采用决策准则mNDWI>0去除上一步结果中的非水体单元；
(24)采用式(4)进行图像二值化处理得到水体地图；
其中，G(i)代表遥感影像上图像单元i处的元素值。
10.如权利要求9所述的一种面向对象的光学与SAR遥感影像融合洪水提取系统，其特征在于：SAR水体图获取模块的具体实现方式如下；
(31)计算初始分割阈值；计算SAR图像VV极化波段的灰度直方图，采用OTSU算法自动确定最优分割阈值，把这一阈值定义为初始分割阈值，用变量T0表示；
(32)设置阈值搜索区间；在灰度直方图上选取初始分割阈值T0左、右长度为ε的BC值域范围[T0-ε,T0+ε]，将其设置为搜索区间；
(33)依次取阈值提取SAR图像的水体；以设定的步长在搜索区间内依次取值T(j)，以T(j)为分割阈值按照式(5)提取SAR图像的水体，采用式(4)二值化得到水体图像，并统计水体图像中水体单元的个数；
其中i表示SAR图像的基本处理单元，G(i)表示遥感影像上图像单元的元素值；公式的含义是将后向散射系数值小于阈值的图像单元划分为水体，否则划分为非水体；
(34)循环步骤(33)，直到遍历完整个搜索区间；
(35)光学图像辅助确定最优分割阈值；计算SAR水体结果图中水体单元个数与同期光学遥感图像提取的水体图像中水体单元个数的差值，取差值最小的T(j)为SAR图像水体提取的最优阈值T；
(36)生成最终水体图；利用最优阈值T提取SAR图像上的水体，通过形态滤波去除散斑噪声因素引起的伪水像元，并按式(4)进行二值化处理，得到最终水体图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310671107.8A CN116704369A (zh) | 2023-06-07 | 2023-06-07 | 一种面向对象的光学与sar遥感影像融合洪水提取方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310671107.8A CN116704369A (zh) | 2023-06-07 | 2023-06-07 | 一种面向对象的光学与sar遥感影像融合洪水提取方法和系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116704369A true CN116704369A (zh) | 2023-09-05 |
Family
ID=87830683
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310671107.8A Withdrawn CN116704369A (zh) | 2023-06-07 | 2023-06-07 | 一种面向对象的光学与sar遥感影像融合洪水提取方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116704369A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117392433A (zh) * | 2023-09-15 | 2024-01-12 | 宁波大学 | 联合sar和光学影像的不同淡水资源类型精细识别方法 |
| CN118097449A (zh) * | 2024-02-23 | 2024-05-28 | 中科星睿科技(北京)有限公司 | 基于多源遥感的洪涝灾害监测方法、装置和电子设备 |
| CN118506207A (zh) * | 2024-07-22 | 2024-08-16 | 广东省科学院广州地理研究所 | 基于sar影像和可见光影像的城市不透水面提取方法 |
-
2023
- 2023-06-07 CN CN202310671107.8A patent/CN116704369A/zh not_active Withdrawn
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117392433A (zh) * | 2023-09-15 | 2024-01-12 | 宁波大学 | 联合sar和光学影像的不同淡水资源类型精细识别方法 |
| CN118097449A (zh) * | 2024-02-23 | 2024-05-28 | 中科星睿科技(北京)有限公司 | 基于多源遥感的洪涝灾害监测方法、装置和电子设备 |
| CN118506207A (zh) * | 2024-07-22 | 2024-08-16 | 广东省科学院广州地理研究所 | 基于sar影像和可见光影像的城市不透水面提取方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Khandelwal et al. | An approach for global monitoring of surface water extent variations in reservoirs using MODIS data | |
| CN116704369A (zh) | 一种面向对象的光学与sar遥感影像融合洪水提取方法和系统 | |
| Shahtahmassebi et al. | Review of shadow detection and de-shadowing methods in remote sensing | |
| Makarau et al. | Haze detection and removal in remotely sensed multispectral imagery | |
| CN111932567B (zh) | 基于卫星影像的冰湖轮廓自动提取方法 | |
| US20140029844A1 (en) | Advanced cloud cover assessment for panchromatic images | |
| CN102750701B (zh) | 针对Landsat TM和ETM图像的厚云及其阴影检测方法 | |
| CN110849814A (zh) | 一种基于多源遥感卫星的遥感影像处理方法 | |
| CN107145891B (zh) | 一种基于遥感影像的水体提取方法及系统 | |
| CN111310771B (zh) | 遥感影像的道路图像提取方法、装置、设备及存储介质 | |
| CN105139396B (zh) | 一种全自动遥感影像云雾检测方法 | |
| CN112418075B (zh) | 一种基于冠层高度模型的玉米倒伏区域检测方法及系统 | |
| CN113744191A (zh) | 一种卫星遥感影像自动云检测方法 | |
| CN114202535A (zh) | 作物种植面积提取方法及装置 | |
| CN110796113B (zh) | 一种基于WorldView-2影像的城市蓝色地物检测方法 | |
| CN116758049A (zh) | 一种基于主被动卫星遥感的城市洪涝立体监测方法 | |
| CN113628098B (zh) | 一种多时相无云卫星影像重构方法 | |
| CN114241333A (zh) | 一种基于多源时序遥感影像精准识别新生滑坡区域的方法 | |
| Wang et al. | Automatic cloud detection in remote sensing imagery using saliency-based mixed features | |
| CN115661675A (zh) | 多云地区地震滑坡遥感识别方法、系统、设备及存储介质 | |
| CN115147728A (zh) | 一种基于雷达数据和光学数据协同的塘坝快速识别方法 | |
| CN113128523B (zh) | 一种基于时间序列遥感影像自动提取珊瑚礁的方法 | |
| CN111696054B (zh) | 基于全极化sar图像的橡胶坝体检测方法 | |
| Wang et al. | Enhancing SDGSAT-1 night light images using a panchromatic guidance denoising algorithm | |
| CN113888421A (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 | ||
| WW01 | Invention patent application withdrawn after publication | ||
| WW01 | Invention patent application withdrawn after publication |
Application publication date: 20230905 |