改进Otsu与形态学相结合的水体信息提取

陈景珏，刘瑞，杨鑫，杨梅，杨远陶

(1.成都理工大学 地球科学学院，成都 610059；2.成都理工大学 地球物理学院，成都 610059)

0 引言

星载雷达卫星具有全天时、全天候、穿透云层的能力，可以利用星载雷达的优势获取多云暴雨情况下的合成孔径雷达(synthetic aperture radar，SAR)影像用于快速洪涝灾害评估。随着SAR数据空间分辨率的提高，出现了种类繁多的超小尺度图像目标，这给图像处理和分析带来了更大的挑战。一方面，由于相干斑点的存在，影响了SAR影像质量；另一方面，建筑物阴影和道路等混乱的目标可能会影响分割结果，因此SAR影像分割是一个公认的难题[1]。

SAR影像分割广泛使用的方法有：最小误差阈值、熵阈值和Otsu等。Otsu也被称为大津法，于1979年由Otsu[2]提出。Otsu因其原理易理解且稳定有效，被广泛用于图像分割[3]。SAR图像中水体和背景内的像元光谱测度各具有同质性，可视为图像中只有水体和背景两类地物[4]，而Otsu可以快速地获取两类目标的分割阈值，因此Otsu可以快速分割水陆边界。基于此特点，Otsu被广泛地应用于SAR图像分割领域[5-8]。虽然Otsu稳定高效，但是在SAR影像中目标与背景的分布通常是不平衡的，这造成Otsu法得到的阈值会偏向方差较大的一类，偏离正确的分割阈值。

SAR卫星是侧式成像，使得图像中包含植物、建筑物和山体所产生的阴影，并且呈现与水体相似的暗色调，造成分割结果中存在一些伪水体[9]。SAR图像相干斑的存在，也会使得分割结果存在杂散点[10]。采用阈值分割法得到的水体分布图，阴影和杂散点是主要的影响对象，剔除它们的影响显得尤为重要，因此对阈值法得到的水体分布图进行优化是准确提取水体的重要技术环节。针对水体在分割结果中大面积连通的特点，利用形态学中的连通区域标定算法除去小面积的杂散点和伪水体，可以有效保留真实水体。但是在海拔较高区域，产生的山体阴影可能比部分真实水体面积要大。在该条件下，使用形态学中连通区标定算法删除伪水体的同时，面积较小的真实水体也会被删除。

综上所述，为解决传统Otsu阈值分割法会偏向方差较大一类的问题，本文提出目标权重改进的Otsu算法。加入数字高程模型(digital elevation model，DEM)进行形态学连通区标定，提出结合高程连通标定算法的形态学优化方法。利用高程连通标定算法删除面积较大的山体阴影，并继续对小面积伪水体进行删除，从而得到精细的水体分布图。

1 研究区域与数据处理

1.1 研究区概况

2020年6—8月长江干流先后发生五次编号洪水，我国27省3 789万人受灾，195.8万人紧急转移安置，直接造成的经济损失高达700亿元，隶属于江西省九江市的彭泽县地区内涝严重，多地村民受灾。连续暴雨、大暴雨引起的山洪暴发或江水倒灌对该地区会造成严重的威胁[11]。

1.2 数据来源

Sentinel-1A卫星于2014年成功发射升空，重访周期为12 d。本文获取了Sentinel-1A两景单视复数(single look complex，SLC)影像，详细参数如表1所示。使用航天飞机雷达地形测绘使命(shuttle radar topography mission，SRTM)90 m分辨率DEM(https：//srtm.csi.cgiar.org/)。

表1 研究数据

1.3 数据预处理

Sentinel-1A影像预处理包括：条带拼接生成SLC数据、多视处理、滤波、地理编码和辐射校正。多视处理牺牲了空间分辨率以提高影像质量[12]。本文使用Frost滤波方法，可以减少斑点噪声的影响，同时保留了边缘信息[13]。地理编码将SAR影像雷达坐标转为地理坐标，辐射校正根据入射角信息计算得出雷达回波后向散射系数。

1.4 哨兵1双极化水体指数计算

贾诗超等[14]提出了哨兵1双极化水体指数(sentinel-1 dual-polarize water index，SDWI)。SDWI将Sentinel-1卫星两种极化影像后向散射系数相乘并扩大十倍，再以自然对数作为函数式。将研究区Sentinel-1A 进行SDWI运算后得到的图像(图1)，经SDWI运算后水体从暗色调转换成了亮色调。本文把Sentinel-1A影像使用SDWI计算得到的图像进行水体提取研究，将ND值量化到256个灰度级。如图2所示，VV和VH极化影像直方图双峰形状基本一致，波谷与第一个波峰起伏较小；SDWI图像的直方图比VV和VH极化影像的直方图双峰现象更为明显。 其中，T1是Otsu计算出地分割阈值；T2是改进Otsu计算出地分割阈值；T3是目视解译地分割阈值。

图1 SDWI图像与湖泊目视解译结果

图2 直方图统计与分割阈值结果

2 研究方法

2.1 改进Otsu算法

周迪等[15]提出了改进Otsu阈值分割方法并命名为KOtsu阈值法。KOtsu考虑了目标类在图像中约占的权重，并结合了类内方差的影响，可以有效地解决Otsu算法会偏向方差大的问题。KOtsu阈值分割方法已成功实践在Lena图像上，本文将该方法应用到SDWI图像，进行阈值分割。经实验发现，类内方差的加入将更依赖权重K的取值，甚至要设置到极低权重值，才能平衡Otsu算法阈值偏向问题。因此，本文对KOtsu进行了改造，如式(1)所示，改造算法没有考虑类内方差，引入了权重K的思想。当使用原始SAR影像进行阈值分割时，水体类在阈值的左边，非水体类在阈值的右边，而经过SDWI运算后，水体类和非水体类则反之。本文提出的改进Otsu与传统的Otsu相比，其主要的变化在于阈值判断公式上的变化，表达如式(1)所示。

T=argt∈{0，1，…，L-1}max(Kω0(t)(μ0(t)-μ(t))2+ω1(t)(μ1(t)-μ(t))2)

(1)

式中：K为目标约占图像的权重，本文权重K取0.2；t为所遍历的分割阈值；L为灰度级， 取值范围为[0,256]；μ(t)为图像总灰度均值；μ0(t)为背景类灰度均值；μ1(t)为目标类灰度均值；ω0(t)为背景类的比例；ω1(t)为目标类的比例；T为最佳分割阈值。

2.2 形态学优化算法

连通区标定算法是形态学中常用的方法之一，该方法在四连通和八连通区域对同质区域进行标定，同时记录标定的坐标(x，y)和标记号i(i=1，2，…，n)，把不符合连通阈值的标记号删除。相干斑造成的杂散点通常由若干像素组成，而水体则是大面积连通区，因此通过设定较小的连通阈值可以把杂散点有效去除。

在Otsu和改进的Otsu的提取结果中，把真实水体、山体阴影和杂散点都分割成了水体，并且山体阴影形成的伪水体比部分真实水体连通面积大，使用连通标定算法把山体阴影去除，部分真实水体也会被随之删除，为解决该问题，本文提出了高程连通区标定算法。

图3(a)为四连通像元表示，图3(b)为八连通像元表示。设改进Otsu算法分割结果用F表示，那么像元值可表示为f(i，j)。F是一幅二值图，用0和1表示，其中0为非水体，1为水体。当f(i，j)=1时，使用四连通或者八连通方法对水体区域进行标定，同时增加了高程h进行判断，即当连通区域的像元满足式(2)时，记录其坐标(i，j)，将F中对应坐标的值进行删除，赋于0值，表达如式(2)所示。

f(i，j)≥h

(2)

式中：i为F的行；j为F的列；h为高程阈值。

图3 连通方法

(3)

当满足式(2)时，连通像元被认为是山体阴影造成误分的伪水体，利用式(3)更新二值图中对应坐标的值，即将值赋给0值，伪水体变更为非水体。

3 结果与分析

基于阈值分割法的水体提取，通常分为两步：第一步，水体粗提取；第二步，水体精提取。文中使用改进Otsu阈值分割法粗提取水体，获得初始水体分布图，然后利用形态学方法优化初始水体分布图，获得最终的水体分布图。

3.1 基于改进Otsu的水体粗提取

图4和图5分别为Otsu提取的水体结果和改进Otsu提取的水体结果。

图4 Otsu水体粗提取

图5 改进Otsu法水体粗提取

由图4和图5可知，两种方法对大面积水体流域都有效地进行了提取，但是在东南和东北山区，Otsu提取的伪水体显然多于改进Otsu。由于水体和非水体在整个研究区域面积分布并不均匀，造成分割的阈值将会偏向非水体，所以Otsu法在东南和东北山区提取的伪水体不仅仅包括了山体阴影还包括了误分的非水体值。直接观察图4和图5可知，改进Otsu提取水体的结果要优于Otsu提取结果。统计两种方法提取水体在研究区的占比，Otsu法提取的水体占整个研究区的59.53%，改进Otsu占19.69%。通过图2中SDWI图像的直方图可知，背景类和目标类是不平衡的，背景类的方差要远大于水体类的方差，此时Otsu算法的计算阈值会偏向背景类，见图2中T1，根据T1阈值进行分割，会出现大量的伪水体。如图1所示，观察SDWI图像可判断出水体在研究区的占比不会超过50%，显然改进Otsu比Otsu提取水体更为合理。

谷鑫志等[16]把目视解译阈值作为水体与非水体的真实阈值，将阈值法计算的阈值与目视解译阈值相比较，差值越小阈值分割精度越高。

目视解译阈值利用ENVI 5.3软件中的 Quick stats 工具和New Raster Color Slice 工具，人机交互式选取了真实阈值，即目视解译的灰度值，目视解译的灰度值见表2以及图2中T3。表2中目视解译阈值是由量化到256灰度级公式通过目视解译灰度值反解计算得出。Otsu和改进Otsu的分割阈值通过遍历图2的直方图，计算出相关的参数后，利用Otsu阈值判断公式和式(1)计算得出。

表2 阈值比较

由表2可知，Otsu法计算的阈值与目视解译阈值相差较大，改进Otsu法计算的阈值比Otsu更接近目视解译阈值，显示了本文改进Otsu阈值分割方法的有效性。综上所述，把改进Otsu水体粗提取图作为初始水体分布图，并使用形态学方法进行优化。

3.2 基于形态学优化的水体精提取

图6显示了采用结合高程连通标定的形态学优化方法对水体的精细提取流程。如图6所示，把改进Otsu法提取的初始水体分布图和DEM作为输入，利用高程连通标定算法，记录大于设置高程阈值的坐标，在粗提取结果中删除被标定坐标位置的伪水体。精提取流程图中，虚线矩形内是文献[4]优化方法，两次连通标定与两次取反的意义：第一次连通标定的作用是删除杂散点和残留的山体阴影；第一次取反把水体转为了背景，非水体转为了目标；第二次连通标定的目的是把水面上分布的船只和杂散点进行删除；第二次取反把真实水体和非水体恢复。

图6 水体精提取流程图

1)高程连通区标定与删除。高程连通标定删除的核心在于高程阈值h的设置。对图5矩形区域，设置了六个高程阈值探索山体阴影、真实水体与设置高程阈值h之间的规则，这六个高程阈值分别为80、90、100、110、120和130 m。处理结果如图7所示，其中图7(a)至图7(f)所示分别为高程阈值80、90、100、110、120和130 m的截取结果，图7(g)为改进Otsu法在矩形区水体粗提取结果。图7(a)至图7(g)椭圆区通过目视解译有三个湖泊，矩形区有两个湖泊。图7(a)椭圆区内，山体阴影和部分杂散点被大量剔除，随着高程阈值的增加,椭圆区域的伪水体数量也逐渐地增加，如图7(b)至图7(f)所示。图7(b)至图7(g)矩形区两个湖泊分别位于左上和右下方，高程阈值为80 m时，两个湖泊被剔除，随着高程阈值的升高，湖泊轮廓也越加明显，提取效果较好的是图7(e)和图7(f)。图7(e)与图7(g)对比可见，矩形右下方湖泊部分被剔除，湖泊轮廓缺失。图7(f)与图7(g)对比可见，矩形区的两个湖泊被完整地保留，并且区域内大部分伪水体被剔除。

图7 规则探索实验

综上所述，设置的阈值越小山体阴影的剔除能力越强，但是不可避免真实水体会被删除或者轮廓不完整，随着高程阈值的增加，山体阴影剔除能力减弱，水体保留效果增强。因此，所设置的高程阈值需要满足两个规则：真实水体能够有效地保留、对山体阴影能有效剔除。

为了获取适合研究区的高程阈值，基于SDWI图像、原始SAR影像，使用ArcGIS软件矢量化图1中目视解译的小中型湖泊并提取了山体阴影。矢量化的湖泊和山体阴影部分区域，如图8所示。

图8 矢量化湖泊与提取的山体阴影

研究区高程范围在-4.546～837.152 m之间。通过遍历研究区高程范围选取高程阈值，对矢量化水体和提取的山体阴影执行高程连通标定删除。为了简化算法的执行时间，高程阈值范围设置为[-5，837]之间，以1为步长，对高程连通标定删除算法进行全局搜索。基于有效保留水体和有效剔除山体阴影两个规则，统计了每一个高程阈值保留水体像素数量和存在山体阴影像素数，如图9所示。图9中，横坐标为高程阈值，纵坐标为像素数，即水体和山体阴影在某个高程阈值被保留的像素个数。从图9可知，真实水体和山体阴影随着高程阈值的增加，水体和山体阴影保留像素数也随之增加。经实验，高程阈值为219 m时，水体被完整保留，见图9中红虚线h1。

图9 高程阈值全局搜索

借鉴机器学习中错误率公式，可以得到水体的损失率以及山体阴影的剔除率。损失率是指真实水体被删除的像素数与总水体像素数之比。剔除率是指山体阴影被删除的像素数与总山体阴影像素数之比。损失率和剔除率可以共用一个公式表示，但是意义却不相同。损失率越大真实水体被保留得越少，剔除率越大山体阴影被删除的能力越强。式(4)即损失率和剔除率的公式。

(4)

式中：S为真实水体或者山体阴影总的像元数；W为处于某个高程阈值时，水体被保留的像元数或者山体阴影存留像元数。

每个高程阈值对应的损失率和剔除率，如图10所示。图10中，p1、h1、h2的值分别为0.1、219、109，当高程阈值为h1时，即高程阈值为219 m时，损失率为0.00%，剔除率为30.33%，可见当取得最小的损失率时，剔除率并不高。当高程阈值取h2时，剔除了率达到了73%，其损失率接近于10%，由此可见，在放弃一定损失率时，能够提高剔除率。表3列出了具有代表性的12个高程阈值与损失率、剔除率的关系。从表3可以看到，当高程阈值为64 m时，剔除率超过了90%，损失率为40.41%;当高程阈值为162 m时，山体阴影的剔除了超过了50%;当高程大于等于832 m时，高程连通标定删除算法失效。

图10 损失率与剔除率

表3 高程阈值与损失率、剔除率

可见，当高程小于832 m时，高程连通标定删除算法对真实水体有一定的保护作用，也兼顾一定的山体阴影剔除能力，但最为有效的高程阈值并不能从表3中得到。利用从表3中的12个高程阈值执行两次连通标定与两次取反的结果计算与矢量化的水体交并比(intersection over union，IoU)。

当IoU取最大值时，即为最佳的高程阈值，表达如式(5)至式(6)所示。

(5)

H=argh∈[h1，h2，…，hn]max(IoU)

(6)

式中：A为矢量化小中型湖泊的矢量化像元；M(·)为水体形态学优化结果；h为高程阈值；l、d为损失率和剔除率；H为最佳高程阈值。

2)两次连通标定和两次取反。如图6所示，两次连通标定与两次取反是水体精提取的处理环节之一。从表3可知，当高程阈值为832 m时，高程连通标定删除算法失效，此时进行两次连通标定和两次取反处理等同于文献[4]优化方法。对表3中12个高程阈值进行两次连通标定与删除，并计算IoU，计算结果如表4所示。表4给出了两次连通标定与两次取反的参数设置。第一次连通值是连通标定算法的连通像元数，当被标记的连通像元小于等于该连通像元时，删除该连通区域。第二次连通值的作用与第一次连通值的作用相同，但是第二次的目的是为了删除水面上的复杂目标，并且水面上复杂目标物，不管处于哪个高程阈值下都不会发生变化，因此该连通值在不同的高程阈值中是不变值。对IoU而言，如图11所示，当高程阈值约为162 m时，得到了最高的IoU。从表3可知，高程阈值为162 m时，损失率为0.72%，剔除率为 50.19%，因此，牺牲较小的损失率，可以极大提高山体阴影的剔除率；剔除率的提升可以改变伪水体的连通区结构，使得可以使用较小的参数删除伪水体，可以有效地保留真实水体。

表4 各高程阈值的IoU与参数设置

图11 交并比曲线

3.3 结果精度评价

本文提出高程连通标定和连通标定算法相结合的形态学优化方法本质上不改变改进Otsu法粗提取水体的边界和大流域连通水体的面积，其作用是把山体阴影、杂散点和水面上船只去除，得到精细的水体分布图。综合Sentinel-1A影像与SDWI图像目视解译了研究区域139个小中型潜在湖泊，并对其进行了矢量化，湖泊位置如图1所示。

从图11可知，随着高程阈值变大，IoU呈先增加后减少的趋势，因此IoU曲线会存在一个极大值点，当高程阈值约为162 m时，IoU达到了最大值，并且第一次连通标定参数设置与剔除率有一定的相关性，如表3、表4所示，随着剔除率的增加，伪水体连通性变差，需要连通删除的阈值也随之变小。通过表4可知，高程阈值选取在108～219 m之间的IoU可以达到0.8以上。可视化高程阈值为162 m时的最终结果和文献[4]结果，如图12所示。

图12 水体提取最终结果

图12(a)为本文选取162 m的高程阈值最终优化结果，图12(b)为文献[4]优化结果。从图12可知，本文方法相比文献[4]优化方法，在大范围水体提取时能够保留更多的细节，对小中型潜在的湖泊有较高提取率。并且从表3可知，高程阈值在74～219 m时，优化结果都高于文献[4]优化方法，这也表明高程连通标定删除是一种有效的优化方法。

4 结束语

本文提出了改进Otsu法和高程连通标定的优化算法，结果表明，改进Otsu法阈值分割精度高于传统Otsu法阈值分割精度，在进行优化时，使用高程连通标定可以很好地保留小中型湖泊，得出结论如下。

使用改进Otsu算法和传统Otsu法对SDWI图像的分割阈值进行了计算。由表2可知，目视解译的分割阈值为226，传统Otsu阈值分割法分割阈值为210，改进Otsu阈值分割方法分割阈值为225。改进Otsu阈值分割结果与目视解译阈值相差1，而传统Otsu阈值分割法相差16，显然改进Otsu法对目标和背景不平衡的情况下，阈值分割效果更显著。

设置了六个高程连通阈值，分别为80、90、100、110、120和130 m，讨论高程连通阈值应满足两个规则。通过矢量化的真实水体和提取的山体阴影，设置步长为1，对高程阈值进行全局搜索，得到当高程阈值大于等于219 m时，可以完整地保留真实水体，但山体阴影剔除能力变弱，当高程大于等于832 m时，高程连通标定算法失效。

所提出的结合高程连通标定算法的形态学优化方法，区别于文献[4]优化方法，是一种新的组合优化方法，通过计算最终结果的IoU，得到了IoU曲线。IoU曲线表明，研究区最佳的高程阈值约为162 m，108～219 m之间的高程阈值IoU在0.8以上。当高程阈值为162 m时，IoU为0.882，文献[4]优化后的IoU为0.608。

本文结合了改进Otsu法和形态学优化方法，提取了彭泽县区域于2020年7月20日的精细水体分布图。本文方法在提取大范围水体时，可以有效解决山体阴影与小型湖泊被一同删除的问题，但仍然有少量小面积湖泊被删除。