基于迭代阈值分割的星载SAR洪水区域快速提取

苗 添, 曾虹程, 王 贺, 陈 杰

(北京航空航天大学电子信息工程学院, 北京 100191)

0 引 言

洪涝灾害以其发生的频率高、影响的范围广、造成的破坏强等特点一直受到科研人员的高度关注。由于洪灾成因复杂,对其进行准确预测较为困难。因此,在洪灾发生后,及时地进行洪灾区域提取对救灾工作有着重要的意义。

目前,遥感卫星技术被广泛应用于洪灾监测中。其中,光学卫星的图像直观、分辨率高,但受观测区域的气象条件影响较大,在复杂的气象条件下难以胜任检测工作,且仅能在白天工作。而星载合成孔径雷达(synthetic aperture radar, SAR)凭借其全天时、全天候工作,分辨率高,不受昼夜、天气影响,检测能力强、范围广等独特的优势,在近些年来被广泛应用到对地环境监测中。基于上述的优点,本文开展了基于星载SAR图像的洪水灾害检测处理研究。

利用星载SAR进行洪水区域提取主要是通过水体提取和图像变化检测来实现的。当前,较为常见的洪水区域提取方法是:利用Otsu阈值分割算法、K-Means聚类算法等图像分割或图像分类方式提取灾前与灾后水体,再通过相减变化检测来提取洪灾区域。这类方法效果稳定、可靠,提取的准确率较高,但也存在局限:由于需要进行遍历或多次迭代,其耗时都相对较长,检测效率较低。针对这一问题,本文提出一种基于改进迭代阈值分割的SAR图像洪水区域快速提取方法,可在保持较高准确率的同时显著缩短算法的处理时间,提升洪灾区域提取的效率。

1 基于SAR图像的洪灾区域传统提取方法

目前,有许多学者对遥感技术在洪灾检测中的应用进行研究,并取得了一些成果。由于光学卫星成像易受天气及水汽云层的影响,因此多采用SAR图像或光学图像与SAR图像融合的方式进行水体提取与洪灾区域检测。

对SAR来说,由于不同地表的粗糙度相差较大,雷达信号在其上的散射方式不同,导致不同地表对应的后向散射系数不同。雷达信号在水体上以镜面反射为主,而在陆地上以散射为主。根据这一点,可利用图像分割或图像分类进行水体提取,再通过图像变化检测提取洪灾区域。传统的水体提取方法有Otsu阈值分割算法、最大熵阈值分割算法、K-Means算法、期望最大化(expectation-maximum, EM)算法等方法。近年来,随着深度学习等领域的发展, U型网络(U-Net)、全卷积网络(fully convolutional networks, FCN)等深度学习模型也在SAR的水体提取中得到了应用。

Otsu阈值分割算法在水体提取中的应用较为普遍。其原理为寻找阈值,使该阈值分割所得的两部分图像的类间方差达到最大值:

(1)

最大熵阈值分割算法与Otsu阈值分割算法原理类似,区别在于其阈值取值为使两部分图像的总熵值达到最大值时的值:

(2)

式中:为图像的总熵值;、分别代表阈值分割后两部分图像的熵值;、分别为两部分图像中各个像素值所占比例;为阈值对应的像素值等级;为像素值等级个数。利用最大熵阈值分割算法提取水体,同样存在时间复杂度高、易受山体阴影影响等缺点。

在水体提取中,常用的图像分类方法有K-Means算法和EM算法。K-Means算法通过多次迭代,修正聚类中心,将原始数据分为若干不同类别,直到准则函数达到最优值,从而实现图像的分类。EM算法的原理是在概率模型中找到参数最大似然估计或最大后验估计的算法。通过多次迭代,实现对目标的分类。

(3)

式中:为待估测值;为观测变量值;为未知变量值;()为似然函数;()为对似然函数取对数的结果。其原理如下:通过多次迭代,改变,使似然函数最大化,优化分类结果。随着机器学习的发展,这两种非监督性学习在SAR图像水体提取中也得到了较多应用,但是由于需要通过多次迭代来进行提取,因而所需时间往往比阈值分割法更长。

近年来,随着相关理论与技术的完善,深度学习在SAR图像处理中得到了广泛应用。在水体提取中,研究者们根据SAR图像自身的特点对经典的深度学习分割模型进行改进,提出了针对SAR图像的基于U-Net、FCN-8S、双边分割网络(bilateral segmentation network, BiSeNet)等分割网络的水体提取算法,取得了较好的结果。戴牧宸等根据SAR图像特点,对传统的双边网络进行改进,减少了原网络中空间路径的卷积层数,采用深度残差网络(deep residual network, ResNet)作为上下文路径骨干网络对SAR图像进行海陆分割,实验表明,其分割效果较好,且泛化能力较强。Shamsolmoali等基于常规的U-Net网络,针对SAR图像的特点进行改进,将残差网络加入传统的U-Net之中。在上下采样路径中,加入紧密连接的残差网络块,对不同尺度的上下文信息进行聚合。实验表明,该方法较传统的U-Net网络而言,提取准确率得到了提高。

此外,相关学者还基于SAR图像与光学图像融合结果开展洪灾区域提取研究,取得了较好的效果,但超出本文研究范围,不再详细说明。

由于洪水检测的特殊性,对算法的时效性和准确性均有较高的要求,上述算法均需较长时间。针对这一问题,本文将开展基于星载SAR图像的洪灾区域快速提取方法研究。

2 基于改进迭代阈值分割的SAR图像洪水区域快速提取方法

迭代阈值分割算法在光学图像处理中得到了较多的应用,但在SAR图像处理中应用较少。由于SAR在成像过程中会受到相干斑噪声及各类畸变的影响,因此部分像素点会受到影响,像素值因此被改变,进而影响到阈值的选取。为了减少异常点对提取效果的影响,需结合SAR图像特点,对噪声进行抑制。基于迭代阈值分割算法的原理及SAR图像处理的相关理论,对迭代阈值分割算法进行改进,以实现对SAR图像洪水区域的快速提取,如图1所示。

图1 基于迭代阈值分割算法的SAR图像洪水区域提取流程图Fig.1 Flow chart of flood area extraction of SAR image based on iterative threshold segmentation algorithm

首先,对灾前和灾后的SAR图像进行预处理,具体包括裁剪、多视处理、辐射校正、几何校正,保证SAR图像的相干斑噪声、辐射畸变、几何畸变等得到抑制;在此基础上,对预处理后的数据进行分级、高斯拟合、抽样,对异常点带来的影响进行进一步抑制;其次,利用迭代阈值分割算法计算所需阈值,实现对SAR图像进行分割处理;然后,开展形态学滤波处理,抑制图像噪声对水体识别带来的影响;最后,基于相减变化检测,实现洪灾区域的有效提取。

2.1 预处理

阈值分割是完全根据图像像素值进行分割的一项技术,因此图像中存在的各类畸变及异常点都会对分割结果带来影响。SAR在成像过程中会受到各类畸变的影响而产生部分异常点,进而影响后续水体的提取及洪灾区域的检测。因此,在水体提取和图像变化检测之前,需要对原始图像进行预处理,以达到抑制各类噪声与畸变的目的。预处理的流程如图2所示。

图2 预处理流程图Fig.2 Preprocessing flow chart

对灾前与灾后图像进行裁剪,使洪灾区域位于图像中央,方便观测的同时减小图像大小,加快处理速率,如图3所示。再依次对图像进行多视处理、辐射校正、几何校正,分别对图像的相干斑噪声、辐射畸变、几何畸变进行抑制,得到预处理后的结果。

图3 原始图像与裁剪后的图像的范围Fig.3 Range of the original image and the cropped image

2.2 改进的迭代阈值分割算法

预处理结束之后,基于迭代阈值分割算法原理,分别对灾前和灾后研究区域的SAR图像进行水体提取。迭代阈值分割算法是基于无限逼近思想提出的一种阈值分割方法。最初,由Ridler和Calvard提出一种阈值计算框架,通过多次选择不断优化阈值,以达到最佳分割效果,如图4所示。

图4 用于迭代阈值选择的示意图图像处理器Fig.4 Schematic image processor for iterative threshold selection

图4中,比较器对输入图像像素与当前阈值进行比较。选择器根据比较器的结果向转换器发送信号,转换器在接受信号后,将对应的输入图像像素点转换为背景图像像素点或目标图像像素点;阈值平均模块负责分别计算两个区域的阈值并进行平均得到新的阈值。当图像中的每个像素点都经过一次处理后,即完成一次迭代。通过多次迭代,当阈值不再变化或两次阈值的差值小于规定值时,即完成迭代阈值分割。

Trussel对该计算框架进行数学归纳,将阈值计算公式归纳为

(4)

式中:表示第次迭代的阈值;表示灰度等级;()表示对应灰度等级的像素点个数。Magid等又对Trussel的归纳进行了证明。

假设阈值满足:

(5)

式中:、分别为背景部分和目标部分的某一特殊数学表示。定义误差函数:

(6)

式中:为各个像素点的灰度值;()为每个像素值等级所占比例。为使误差函数最小,对其求导,得

(7)

计算可得

(8)

式中:表示背景部分图像的平均灰度等级。同理,可得

=

(9)

式中:表示目标部分图像的平均灰度等级。至此,迭代阈值分割算法的证明完成。

迭代阈值分割算法通过多次迭代不断优化阈值,当相邻两次阈值相同或其差值小到预定范围时,确定最后一次迭代的结果为最终阈值。迭代阈值分割算法流程图如图5所示。

图5 迭代阈值分割算法流程图Fig.5 Flow chart of iterative threshold segmentation algorithm

如图5所示,其计算流程如下:将高斯拟合并采样后的数据输入,将初始阈值设置为所有像素点的均值,并进行初次分割。在分割后,分别计算背景部分与目标部分的平均像素值,并将二者的平均值作为新的阈值, 比较与,当二者不相等时,继续迭代,当二者相等时,完成阈值的计算。

这一算法对于像素值呈双峰状分布的图像有良好的分割效果。常见的Otsu阈值分割算法和最大熵阈值分割算法由于需要对所有像素等级进行遍历,因此其时间复杂度较高,采用常规方法实现上述两种算法,复杂度为()。而迭代阈值分割算法不需要对所有像素值进行遍历,仅通过数次迭代即可达到较好的分割效果,计算速度相对更快,效率更高。

在进行SAR图像水体提取时,由外界影响导致的异常点往往会影响提取效果,针对这一问题,提出一种改进的迭代阈值分割算法。由于水体的后向散射系数较低,陆地的后向散射系数较高,因此在理想状态下,包含较大区域水体的SAR图像后向散射系数分布为双峰状,当山体阴影影响较大时,后向散射系数的分布会产生畸变,使个别等级的像素点个数异常。针对SAR图像自身易受山体阴影影响、存在相干斑噪声等干扰的问题,结合SAR图像中像素值的分布特征,可先对SAR图像的像素值进行分级量化,对量化后的数据进行拟合再抽样,对该类畸变进行抑制。之后,利用抽样所得数据进行迭代阈值分割,最后再进行形态学滤波,对噪点进行抑制。其流程如图6所示。

图6 改进的迭代阈值分割算法流程图Fig.6 Flow chart of improved iterative threshold segmentation algorithm

具体步骤如下:

首先,将预处理后的SAR图像中每个像素点的像素值转化为分贝形式：

=10lg

(10)

式中:为原始像素点的像素值;为转化后dB图像对应像素点的值。对转化后图像的像素点值进行分级、量化,分为个等级,并统计每个等级的像素点个数。之后,对其进行拟合。目前,在变化检测技术发展过程中,基于参数模型的方法应用较广。常用的模型有高斯模型、广义高斯模型、瑞利-高斯模型等。考虑到高斯分布的灵活性和稳定性,且多种涉及被动传感器的数据分布模型均为高斯模型或广义高斯模型,本文中的SAR图像像素值分布图采用高斯函数进行拟合,本文选用高斯函数对数据进行拟合,拟合函数如下:

(11)

式中:()为高斯拟合函数;为原始的离散数据;、、为高斯函数中的常量;为拟合阶数。依据拟合情况调整拟合阶数,以达到最优拟合效果。完成拟合后,进行抽样,将抽样结果存入×1的一维矩阵之中。然后,再对抽样结果进行迭代阈值分割。

迭代阈值分割算法的总体思路为:首先设置初始阈值,再不断对其进行调整、迭代,每一次调整的输入参数为上一次迭代中得到的结果,经过数次迭代之后,最后产生一个理想的阈值。

其算法实现步骤如下:

首先,将图片中的像素点依据其幅值进行排序,将初始阈值设置为所有像素点像素值的平均值:

(12)

然后,按照初始阈值将原始图像分成1、2两部分,分别计算两部分的后向散射系数差值的平均值:

(13)

(14)

式中:代表像素点的值;为未转化前像素点的值;表示其对应的概率;表示阈值对应的像素等级。

接着将二者的平均值作为新的阈值:

(15)

对新的阈值与之前计算的阈值进行对比。若与的值相同,则为所求阈值;若不相同,则继续计算,不断迭代,直至相等或两次阈值的差值小于设定值时,最后一次迭代的阈值即为所求阈值。根据阈值,将原图像转化为二值图像。所得二值图像即为洪水区域的提取结果:

(16)

式中:为原始像素值;为阈值;()为处理后的像素值。

2.3 形态学滤波与洪水区域提取

由于阈值分割法在处理时会忽视图像中蕴含的空间信息,因此部分没有得到完全抑制的异常点会给算法的判断带来影响,从而产生离散的噪点。针对这一问题,在完成水体初步提取后,采用形态学滤波对噪点进行抑制。由于噪点分布较为分散,为了减少形态学滤波对水体提取效果的影响,尽量选择尺寸较小的处理元素进行处理。本文中,选取尺寸为3×3的矩形进行形态学滤波。先进行腐蚀运算,再进行膨胀运算。腐蚀运算以处理单元的中心点为锚点,在原图像中移动处理单元,遍历图像的每一个像素。取结构元素覆盖下的原图对应区域内的所有像素的最小值,用这个最小值替换当前像素值,得到结果;膨胀运算与腐蚀运算相反,取结构元素覆盖下的原图对应区域内的所有像素的最大值,用最大值替换当前像素值,得到结果。经过形态学滤波,图像的噪声得到抑制,灾前与灾后的水体得到提取。

洪水区域提取的最后一步,对灾前、灾后图像进行变化检测。在之前的处理中,实现了对灾前与灾后水体的分别提取。在此基础上,利用相减变化检测可得到灾后相比灾前扩张的水体部分,此部分即为洪灾区域。通过先提取水体再进行变化检测的方式,可抑制零星噪点对变化检测结果的影响。具体流程为:利用灾后水体提取图像减去灾前水体提取图像,对应相随点值相减。对于未变化地区,相减后的结果为0,对于遭受洪涝灾害的地区,相减后结果为255。所得二值图像即为遭受洪灾区域的提取结果。

3 实验结果与分析

3.1 实验数据

2020年7月,鄱阳湖流域发生特大洪水,给当地群众的生命财产安全带来了严重威胁。本文基于星载SAR图像对2020年7月鄱阳湖流域的洪水进行研究。本文选用欧空局哨兵-1卫星的数据进行相关研究。选用的SAR图像为Level-1级别,干涉宽幅(interferometric wide swath, IW)模式、地距多视影像(ground range detected, GRD)类型的单极化成像。本次使用的数据信息如表1所示。

表1 原始图像数据信息Table 1 Original image data information

3.2 实验结果

经过预处理的图像如图7所示。

图7 预处理后的SAR图像Fig.7 Preprocessed SAR image

提取预处理后SAR图像像素点数据,按照其像素值大小,均匀分为500个等级,并统计每个等级像素点个数。之后利用高斯分布对数据进行拟合,对拟合效果进行评估。分别计算不同阶数拟合函数的拟合系数与均方根误差(root mean square error, RMSE),如图8所示。

图8 不同阶数高斯拟合的效果Fig.8 Effect of different orders of Gaussian fitting

由图8可知,数据1与数据3应采用六阶高斯函数,数据2采用五阶高斯函数。高斯拟合后的结果如图9所示。常量取值如表2所示。

图9 原始数据与拟合结果Fig.9 Original data and fitting results

表2 高斯拟合函数参数Table 2 Parameters of Gaussian fitting function

利用迭代阈值分割算法分别计算3幅图像的阈值,计算结果如下:2019年7月20日的SAR图像阈值等级为119,图像中对应像素点的值为-16.18 dB;2020年3月16日的SAR图像阈值等级为109,图像中对应像素点的值为-16.07 dB;2020年7月14日的SAR图像阈值等级为122,图像中对应像素点的值为-15.36 dB。依据阈值对图像进行分割,得到初步提取的结果,再对其进行形态学滤波,滤波前后结果如图10所示。

图10 形态学滤波前后水体提取结果Fig.10 Water area extraction before and after morphological filtering

对提取水体后的结果进行统计:数据1的SAR图像中,提取水体所占比例为23.31%;数据2的SAR图像中,提取水体所占比例为18.45%;数据3的SAR图像中,提取水体所占比例为27.81%。与数据2相比,鄱阳湖流域的水体面积在7月增加了50.73%;即使与数据3,即2019年同期相比,鄱阳湖流域的水体面积也增加了19.30%。

分别对数据3与数据2、数据3与数据1的提取结果进行相减变化检测,并将结果与灾前、灾后叠加形成的RGB (red green blue)图像进行对比,其结果如图11所示。

图11 洪灾区域提取结果对比Fig.11 Comparison of extraction of flood area

在RGB图像中,红色区域即为遭受洪涝灾害地区,由图11可知,通过水体提取和相减变化检测,洪灾区域基本实现了提取。根据相减变化检测结果可知,与2020年3月相比,在2020年7月,鄱阳湖流域的水体面积有明显的增大,鄱阳湖流域的西南部分扩张尤为明显,鄱阳湖的北部区域扩张相对较少;与2019年同期相比,鄱阳湖也扩张了一部分,扩张部分集中在鄱阳湖西部区域和东南区域。

为评估水体提取结果,利用Landsat-8光学卫星开展验证实验。由于2020年7月鄱阳湖流域雨水较多、云层较厚,光学卫星难以观测,因此选用2019年7月Landsat-8的图像与提取结果进行对比,如图12所示。

图12 Landsat-8光学卫星2019年7月鄱阳湖区域图像Fig.12 Landsat-8 optical satellite image of Poyang Lake area in July 2019

除了与迭代阈值分割算法进行对比外,本文还选用Otsu阈值分割算法、K-Means非监督性学习、EM算法的提取结果进行对比,如图13所示。

图13 不同算法的提取结果Fig.13 Extraction results of different algorithms

为验证本文所提方法的有效性,通过采样评估对上述3种提取方式的提取效果进行评估。具体方法为:在光学卫星图像中分别随机选取300个点,分别统计3种提取方式的准确率、虚警率、漏警率、计算时间等。其结果如表3所示。

表3 不同提取方式评估结果Table 3 Evaluation of different extraction methods

根据随机采样后的对比结果可知,对于2019年7月鄱阳湖流域的SAR图像,迭代阈值分割法与Otsu阈值分割法的阈值较为接近,两者的提取效果差距不大;K-Means非监督性学习提取的结果与前二者相比,准确率较低、虚警率较高、漏警率较为接近。EM算法的准确率最高,虚警率和漏警率都较低;在计算时间方面,本文使用处理器为Intel Core i7-9850H CPU的计算机完成实验, 其中迭代阈值分割算法完成对2019年7月鄱阳湖流域SAR图像数据的处理需要1.045 s,而Otsu阈值分割算法需要10.016 s,K-Means非监督性学习需要78.153 s。Otsu阈值分割算法由于需要对所有后向散射系数等级进行一次遍历,因此算法的时间复杂度较高,导致所需时间较长。而对于K-Means算法和EM算法,均需要进行多次迭代,所需时间远大于前两种方法。对于迭代阈值分割算法,在一般情况下,仅需迭代数次即可得出阈值,所需时间大大减少,提取效率也相应得到了提高。

4 结 论

本文提出了一种基于改进迭代阈值分割的SAR图像洪水区域快速提取方法。该方法首先通过预处理与高斯拟合对噪声和畸变进行抑制,之后通过迭代阈值分割算法提取水体,再利用形态学滤波抑制噪声,最后通过变化检测提取洪水区域。与传统的洪水区域提取算法相比,该算法在保证较高提取准确率的同时提高了处理速度,使提取效率得到提升。将该方法应用于2020年7月鄱阳湖流域的洪水灾害研究中,通过预处理、水体提取、变化检测等流程,得到了鄱阳湖流域洪水区域的提取结果。通过将结果与Landsat-8光学卫星图像的对比,该方法的可靠性与准确性得到了验证;通过与传统提取方式的对比实验,表明了迭代阈值分割算法在洪水区域提取中的高效性,为利用SAR图像提取洪水区域提供了新的研究思路。

后续在这方面的主要工作内容为:① 对算法的适用性进行进一步检验,考虑复杂情况下的洪涝区域提取,如城市内涝等;② 考虑结合光学图像信息,改进算法结构,提升准确率。