基于改进SLIC 算法的Landsat 8 OLI 遥感影像水体提取

方 悦，颜 普，陈 杰

（1.安徽建筑大学 电子与信息工程学院，安徽 合肥 230601；2.安徽省国际古建筑智能感知与高维建模联合研究中心，安徽 合肥 230601）

在所有自然灾害中洪涝灾害发生频率很高［1］，通过提取水体信息、精确绘制水体图来观察水体变化，对灾害救援和评估具有重要意义。与实地勘测相比，基于遥感影像提取水体信息可以节约大量劳动力成本、提高水体信息提取精度［2－3］。近年来，随着卫星技术快速发展，遥感影像数据种类越来越多，其中Landsat 8 OLI 遥感影像数据可免费获得，其具有中等空间分辨率，可满足水体提取要求，因此被广泛用于水体提取［4］。

在遥感影像水体提取方法中最常用的是水体指数法［5－6］，但该方法十分依赖阈值选择，需要重复大量实验才可确定最佳阈值。简单线性迭代聚类（SLIC）算法是一种高效的超像素分割算法［7］，其将图像中的像素点聚类到超像素中，具有优异的分割性能［8］。已有多数研究者利用SLIC 算法在遥感影像中提取水体，湛南渝等［9］采用SLIC 算法在哨兵1 号雷达数据中提取水体，并对SLIC 算法提取信息与Otsu 阈值法提取信息进行比较分析，发现SLIC 算法的提取效率与精度更高。但采用SLIC 算法进行像素相似性度量时只考虑影像颜色特征和空间特征，这对复杂遥感影像分割是不利的。遥感影像中往往有很多颜色接近的地物，如水体一般呈现深色调，其与黑色建筑物表面颜色很接近，因此需要加入纹理特征来精确分割此类地物。Lu等［10］将灰度共生矩阵纹理特征纳入SLIC 算法，针对高分二号全色图像、多光谱图像和融合图像验证改进算法的有效性，结果表明，改进算法的分割完整性和正确性都高于SLIC 算法的。此外，SLIC 算法常出现过分割问题，需要利用区域合并方法来解决此问题，区域邻接图（RAG）［11］、近邻图（NNG）［12］和层次聚类法［13］是3 种主要的区域合并方法。相较于区域邻接图和近邻图，层次聚类法的输入参数较少。Wu 等［14］将LBP纹理特征纳入层次聚类法，验证了该方法具有良好的区域合并性能。因此，笔者也选择将纹理特征纳入层次聚类方法进行区域合并。

为将水体与相似颜色地物分割开，笔者结合线性特征增强方法与局部三值模式（LTP）纹理特征［15］提出一种新的纹理特征（E＿LTP），将E＿LTP 纹理特征纳入SLIC 算法，得到新超像素分割方法（F－SLIC），之后对水体超像素进行区域合并得到完整水体，通过提取同一区域不同时间的水体，分析水体变化来检测洪涝灾害，以期为洪涝灾害的监测和预警提供更加可靠的手段。

1 研究区域选取和数据预处理

1.1 研究区域选取

选择3 个区域进行研究，3 个区域涉及黄河及其3条主要支流（分别为位于河南省巩义市的伊洛河、陕西省渭南市的渭河、甘肃省兰州市的湟水），3 个区域的伪彩色图像见图1。各图像都包含具有线性特征的大水体和小水体，小水体较为狭窄，增大了水体分割的难度。

图1 伪彩色图像（5／6／4 波段组合）

1.2 Landsat 8 OLI 遥感影像和Ground truth 图像

选用的Landsat 8 OLI 遥感影像数据源自地理空间数据云网站（https：∥www.gscloud.cn／），数据已经过几何校正和地形校正，具体信息见表1。Landsat 8 OLI遥感影像数据往往包含云层，为了避免云阴影导致的水体分割误差增大，裁剪出无云阴影图像作为研究区域。

表1 研究区域Landsat 8 OLI 遥感影像数据信息

把高分辨率Google EarthTM图像作为辅助信息，对Landsat 8 OLI 遥感影像数据的“真实”水体边界进行人工数字化标定，以获得Ground truth 图像用来评估分割性能。

1.3 数据预处理

仅仅经过几何校正和地形校正的Landsat 8 OLI遥感影像数据是数字图像，在计算水体指数之前需要先对影像数据进行辐射定标，将影像呈现的原始像元值转换为对应像元的地表反射率。因大气吸收和散射机制等因素会影响地表反射率的精确性，需对辐射定标后的数据再进行大气校正处理，大气校正处理模型是Flaash 模型［16］。

2 研究方法

水体提取流程见图2，主要环节有：计算水体指数WBAWI，利用F－SLIC 超像素分割算法对WBAWI 图像进行分割，基于E＿LTP 纹理特征和颜色特征对超像素进行区域合并得到完整水体。

图2 水体提取流程

2.1 水体指数WBAWI

归一化差异水体指数（NDWI）和修正归一化差异水体指数（MNDWI）等传统水体指数对小水体的增强效果较差，因此通过Kautlr－Thomas（KT）变换、LBV 变换、自动水体提取指数（AWEInsh）和HIS 变换提出一个新水体指数（WBAWI）。KT 变换可产生绿度、亮度、湿度3 个分量。LBV 变换可产生地物总辐射水平（L）、地物可见光－红外光辐射平衡（B）和辐射变化矢量（V）3 个分量。与其他地物相比，水体的红外辐射最弱，采用B 分量可以突出水体与其他地物的差异，显著反映水体特征。AWEInsh可以有效突出深色建筑物与水体的差异，减小误分类概率。

自动水体提取指数计算公式为

式中：A为自动水体提取指数AWEInsh，ρ3、ρ5、ρ6、ρ8分别为Landsat 8 OLI 遥感影像数据的波段3（Green 波段）、波段5（NIR 波段）、波段6（SWIR1 波段）、波段8（SWIR2 波段）的反射率。

将湿度分量、B 分量和AWEInsh分别作为RGB 图像的红色通道、绿色通道和蓝色通道，然后对RGB 图像进行HIS 变换，变换得到的I 分量即为WBAWI 图像。兰州区域WBAWI、NDWI、MNDWI 水体指数图像见图3，从图上红色矩形框中可看出，WBAWI 相比其他2 个传统水体指数，小水体与相邻地物的对比度更高，更有利于后续水体分割。

图3 3 种水体指数图像

2.2 F－SLIC 超像素分割

2.2.1 LTP 纹理特征提取

在SLIC 超像素分割聚类过程中度量周边邻域像素与聚类中心像素的相似性时，只考虑颜色特征和空间特征，这对复杂遥感影像的水体分割很不利，因此可在相似性度量中加入纹理特征距离，得到高度相似的超像素块。

LBP 容易受噪声影响［11］，LTP 是由LBP 扩展的三值模式，为了提高描述符的稳健性，选择LTP 算子来提取WBAWI 图像的纹理特征。

LTP 计算公式如下：

式中：xj、xk分别为3×3 邻域的中心像素、周边像素，t为阈值，s（uk）为LTP 编码8 个位置的值，MLTP为LTP编码，s（u1k）、s（u2k）分别为正模式编码、负模式编码8个位置的值，MLTP正、MLTP负分别为正模式编码、负模式编码。

经过实验发现正模式编码得到的纹理特征更适合用来区分水体和非水体，因此本文选用正模式编码提取LTP 纹理特征。

2.2.2 E＿LTP 纹理特征

为了使WBAWI 图像的纹理特征突出线性水体与非水体的差异，需要对WBAWI 图像进行线性特征增强，再提取E＿LTP 纹理特征。

线性特征增强需先使用0°、45°、90°和135°这4个方向的算子对图像进行滤波。滤波后图像的每个像素在4 个方向的最大值所组成的图像就是线性特征增强图像，见图4。对比兰州区域WBAWI 图像的LTP与E＿LTP 纹理特征，从图上红色矩形框内可明显看出，E＿LTP纹理特征的水体与相邻地物的对比度较高，有利于后续水体分割。

图4 兰州区域WBAWI 图像的两种纹理特征

2.2.3 F－SLIC 算法

复杂遥感影像中有很多与水体颜色相似的地物，为了解决水体与相邻地物难分割的难题，结合E＿LTP 纹理特征和SLIC 算法对遥感影像进行超像素分割。

SLIC 算法有2 个输入参数，即预期的超像素个数（初始聚类中心）k和超像素紧密度m。在对灰度图像进行SLIC 超像素分割时，从初始化步骤开始，k个初始聚类中心分别分布在间隔S个像素的规则网格上，网格间隔S＝N／k，其中N为图像的总像素数。为了避免聚类中心处于图像边缘或者噪声上，将聚类中心移到其3×3 邻域最低梯度位置上，之后在聚类中心周围2S×2S区域中搜索相似像素。

SLIC 算法中原像素相似性度量公式为

式中：dgray为灰度距离；gc为聚类中心c 的灰度值；gi为像素i的灰度值；（xc，yc）是聚类中心c 的坐标；（xi，yi）是像素i的坐标；dxy为空间距离；D′s1为像素i与聚类中心c 之间基于灰度特征和颜色特征的距离；Mgray为类内最大灰度距离，设定为输入参数紧密度m；Mxy为类内最大空间距离，Mxy＝S＝N／k；Ds1为D′s1对应的最终的相似性度量。

在F－SLIC 算法中提出一种新相似性度量，在原来只考虑颜色特征和空间特征的基础上，加入纹理特征距离，纹理特征距离计算公式为

式中：Ec为聚类中心c 的E＿LTP 纹理特征，Ei为像素i的E＿LTP 纹理特征，dE＿LTP为E＿LTP 纹理特征距离。

使用类内最大纹理特征距离将纹理特征与颜色、空间特征合并为一个度量值，新相似性度量Ds′对应的最终相似性度量Ds计算公式为

式中：ME＿LTP为类内最大纹理特征距离为像素i与聚类中心c 之间基于灰度、颜色和纹理特征的距离。

每个像素都与距离最小的聚类中心聚类，一遍聚类之后，将当前超像素的平均矢量更新为新的聚类中心，以此类推，重复迭代10 次。通过后处理把很小的超像素分配给相邻超像素以加强连通性。总体来说，F－SLIC 算法相比于传统SLIC 算法不同的是，F－SLIC算法结合E＿LTP 纹理特征改进了聚类相似性度量。

2.3 基于E＿LTP 和灰度特征的区域合并

F－SLIC 超像素分割算法得到的分割结果并不是完整水体提取结果，此分割结果属于过分割，因此需要通过区域合并将相似度高的超像素迭代合并为一个超像素来解决过分割问题。

在基于层次聚类的区域合并算法中，相似性度量一般为颜色距离。但应用于遥感影像时，水体颜色与相邻地物颜色过于相似，仅仅使用颜色距离会导致误合并，因此选择将E＿LTP 纹理特征纳入相似性度量以优化合并性能。

选用余弦相似度计算超像素与相邻超像素之间的灰度特征和纹理特征距离，即把图像表示为矢量，矢量间的余弦距离即为图像间的相似度。首先将两个超像素缩放至大小一致，获得超像素灰度或纹理特征的直方图分布，再将直方图划分为64 个区，每个区包括连续4 个灰度等级，对每个区的4 个值进行求和得到1个数据，最终得到的64 个数据为超像素灰度或纹理特征对应的矢量，这样就把两个超像素的灰度值或纹理特征转化成矢量A和矢量B，见式（15）。按照式（16）计算2 个矢量的余弦距离，余弦值cos（A，B）越接近1，表明夹角越接近0°，2 个矢量越相似，即2 个超像素越相似。

式中：ai为矢量A中的64 个数据；bi为矢量B中的64个数据；Sgray为超像素与相邻超像素的灰度特征余弦距离；Agray、Bgray分别为2 个超像素灰度特征对应的矢量；AE＿LTP、BE＿LTP分别为2 个超像素E＿LTP 纹理特征对应的矢量；SE＿LTP为超像素与相邻超像素的E＿LTP 纹理特征余弦距离；D为相似性度量，同时考虑灰度特征和E＿LTP 纹理特征距离。

超像素只与相邻的余弦相似度最高且达到阈值的超像素合并，通过不断迭代，直到所有余弦相似度较高的超像素合并为一个超像素。把WBAWI 图像中过分割的水体超像素合并之后就可获得完整的水体提取结果，再通过比较同一区域不同时间的水体变化来监测该区域是否发生洪涝灾害。

2.4 算法评价指标

采用边界召回率［17］、欠分割错误率［18］验证FSLIC 算法的有效性和实用性。边界召回率越大，表示算法分割的边界越准确。边界召回率BR计算公式为

式中：TP为算法分割边界像素出现在Ground truth 边界周围两个像素区域内的像素总数，GB为Ground truth边界像素总数。

欠分割错误率EU是算法分割得到的超像素相比Ground truth 分割结果的“溢出”。EU值越小，表示算法分割结果相比Ground truth 分割结果的“溢出”越少，分割结果越准确。EU计算公式为

式中：si为算法分割结果的超像素；gk为Ground truth 分割结果的超像素；为取与si交集最大的gk。

3 结果与分析

3.1 5 种算法对遥感影像分割结果的比较

为了验证F－SLIC 算法对Landsat 8 OLI 遥感影像的分割性能，采用SLIC 算法、融合LTP 纹理特征的SLIC 算法、F－SLIC 算法分别在巩义区域、渭南区域和兰州区域进行超像素分割处理。这3 个研究区域的最佳分割超像素块数分别为3 000、2 000、2 000，为了分割出呈窄线性的小水体，紧密度m取值都较大，分别为50、30、33。

采用不同算法在最佳分割超像素块数和最佳紧密度下得到的超像素分割图像见图5～图7（分割图像上方图为对应颜色方框的局部放大图，绿框、蓝框和黑框从左至右均依次对应SLIC 图像、SLIC ＋LTP 图像、F－SLIC图像的放大图），3 种算法的输入参数统一。在3 个研究区域内，算法性能差异主要体现在小水体上，从绿框和蓝框局部放大图可看出，SLIC＋LTP 算法和F－SLIC算法的水体边界分割性能都明显比SLIC 算法的好，SLIC 算法存在较多的欠分割情况。从黑框局部放大图可看出，F－SLIC 算法的分割性能比SLIC＋LTP算法的稍好一些，SLIC＋LTP 算法也存在欠分割情况。总的来说，F－SLIC 算法的分割性能最好，产生的水体超像素很规整。

图5 采用不同算法得到的巩义区域超像素分割图像

图6 采用不同算法得到的渭南区域超像素分割图像

SLIC、SLIC ＋LTP、F－SLIC 三种算法的边界召回率、欠分割错误率见图8～图10，F－SLIC 算法与SLIC＋LTP 算法的分割性能较接近，均优于SLIC 算法。F－SLIC算法的边界召回率比SLIC 算法提高了1%～4%，F－SLIC 算法的欠分割错误率比SLIC 算法降低了0.05%～0.30%，超像素块数越少，F－SLIC 算法的两项指标值与其他两种算法的差异越明显。

图9 渭南区域不同算法的边界召回率和欠分割错误率

图10 兰州区域不同算法的边界召回率和欠分割错误率

3.2 水体变化监测

采用F－SLIC 算法分割的各区域水体超像素区域合并后的水体提取结果见图11，水体轮廓图箭头左边图像为对应颜色方框的局部放大图。可以看出，提取水体轮廓较完整，因此可用F－SLIC 算法提取同一区域不同时间的水体来比较水体变化。

提取的不同时间的水体轮廓见图12，水体轮廓图箭头左边图像为对应颜色方框的局部放大图。图12（a）中红色轮廓表示提取巩义区域2021 年7 月12 日的水体轮廓，绿色轮廓表示提取巩义区域2021 年9 月30 日的水体轮廓；图12（b）中红色轮廓表示提取渭南区域2020 年1 月29 日的水体轮廓，绿色轮廓表示提取渭南区域2020 年7 月7 日的水体轮廓；图12（c）中红色轮廓表示提取兰州区域2021 年3 月23 日的水体轮廓，绿色轮廓表示提取兰州区域2021 年8 月14 日的水体轮廓。可以看出，3 个研究区域的水体随时间变化都发生不同程度的变化，其中巩义区域和渭南区域的黄河水体变化较明显，对于这些水体发生明显变化的地区可做标记，在下次遇到连续暴雨天气时进行洪涝灾害防控，以降低对附近居民的伤害。

图12 不同时间水体轮廓对比

4 结论

本文将E＿LTP 纹理特征纳入SLIC 算法的相似性度量，提出一种新超像素分割方法（F－SLIC），对Landsat 8 OLI 遥感影像进行超像素分割。经与SLIC算法和SLIC＋LTP 算法对比，F－SLIC 算法在边缘召回率和欠分割错误率上都得到了优化，能够很好地将水体和背景分割开，并解决了小水体难以分割的问题。此外，对水体超像素进行区域合并，解决了超像素过分割问题，得到完整水体轮廓。通过提取同一区域不同时间的水体轮廓，可有效观察水体变化，提前对水体变化明显的区域做标记，为洪涝防控提供参考。