水库水体面积提取方法

李亚琦，余宇峰

(河海大学 计算机与信息学院，江苏 南京 211100)

0 引 言

随着遥感技术的不断发展，遥感影像的获取变得越来越便捷，而遥感影像水体提取技术也在水文监测、环境保护等方面得到了广泛的应用[1,2]。目前，较为常用的水体提取方法主要包括监督方法和非监督方法两种类型。监督方法主要包括决策树算法[3]、支持向量机(support vector machine,SVM)[4]和深度学习方法[5]，监督方法主要通过对已标注样本的学习来实现图像中水体与非水体的分类，经过对足够多的样本进行学习，监督方法可以通过学习到的水体特征对图像进行水体提取。使用监督方法进行水体提取的效果良好，精度也比较高，但是要使用监督法进行遥感图像的水体提取首先就需要获得足够多的样本数据进行训练，而遥感图像本身的像辐较大，这就为数据样本的训练带来了一些困难，而且在进行数据样本标注时往往需要专业人员对遥感图像进行目视解译，这种人工解译的方法不仅需要大量的人力，也会花费大量的时间。非监督方法主要包括水体指数法[6,7]、波段阈值法[8,9]和谱间关系法[10]。非监督方法计算简单，运算速度快，但是在水体提取时容易受到水体阴影、植被等的影响，存在一定的误提取现象，不能很好满足水体提取的精度需求。

针对上述水体提取方法存在的一些不足，本文提出一种基于遗传算法和改进Otsu算法的水体信息提取方法，并使用最大连通域算法对图像中的水库水体信息进行提取。

1 相关方法研究

1.1 NDWI水体指数法

NDWI水体指数法主要利用遥感图像中水体在不同波段间反射率和吸收率的差异来凸显水体信息，水体在近红外波段的吸收率最强[11]，反射率最弱，而在绿光波段的吸收率最弱，反射率最强，植被则恰恰相反，NDWI水体指数法就利用这种差异，用绿光波段的DN值和近红外波段的DN值进行归一化差值处理，从而区分出图像中的水体与非水体。NDWI水体指数的计算公式为

(1)

其中，ρ(Green)表示绿光波段的反射率，ρ(NIR)表示近红外波段的反射率。NDWI水体指数能够有效抑制非水体光谱对水体提取的影响，使水体信息较明显的区别于其它地物信息，有利于后续对水体的进一步提取。

1.2 Otsu算法

Otsu算法又叫做大津法或最大类间方差法，是日本学者大津(Nobuyuki Otsu)提出的一种求取全局阈值的算法，该算法计算简单且速度快，受到图像亮度和对比度的影响较小。

(2)

假设存在阈值K，这个阈值K将图像分为两部分，即灰度值大于等于K的部分和灰度值小于K的部分，将灰度值在0到K的部分作为A类，灰度值在K+1到L-1部分的作为B类，则这两部分出现的概率分别为

(3)

A类像素和B类像素的灰度均值为

(4)

图像总体灰度均值为

(5)

类A和类B的类间方差为

σ2=PA(μA-μ)2+PB(μB-μ)2

(6)

Otsu算法的原理就是遍历所有的灰度值，找出使得σ2最大的灰度值作为K值，即最佳阈值，根据最佳阈值将图像分为目标像素和背景像素。但是当图像较大时，需要遍历的灰度值数量较多，计算效率就会降低，且Otsu算法对图像噪声较为敏感，会对最终的结果产生影响。

2 水库水体面积提取方法

2.1 改进的Otsu算法

传统的Otsu算法只考虑了灰度直方图呈现单峰分布的图像，当图像的灰度直方图呈现双峰或多峰分布时，阈值分割的效果并不理想。通过对遥感图像光谱特征和NDWI水体指数的研究，发现NDWI水体指数处理后的图像灰度直方图多呈现出双峰分布的特征。为了使得Otsu算法更加适用于遥感图像水体信息的提取，对Otsu算法进行了改进，并使用遗传算法进行优化。

为了使得Otsu算法更加适用于遥感图像水体信息的提取，在计算最佳阈值时，设置最佳阈值的数量为2，将计算出的两个最佳阈值作为图像分割的依据。在进行图像分割时并不直接利用最佳阈值，而是先采用滑动窗口判断窗口内的像素点，若窗口内像素点的灰度值位于两个最佳阈值之间，则该像素点不使用最佳阈值进行分割，而是利用自适应阈值分割算法进行判断。自适应阈值分割算法更加适用于局部阈值分割，可以根据邻域内的像素点进行分割。

遗传算法是模拟达尔文进化的一种优化算法，采用概率化的参数寻优方法，可以自适应进行寻优操作。使用遗传算法训练改进后的Otsu算法，基本思想就是将Otsu算法的类间方差作为遗传算法的适应度函数，类间方差的值越大，适应度越高，所取得的结果也就越好。首先将灰度值进行编码形成种群，再经过选择、交叉和变异操作，用适应度函数来表示个体的优劣程度。使用遗传算法训练改进的Otsu算法的具体步骤为：

(1)将原始遥感图像使用NDWI水体指数方法进行处理，对处理后的图像进行灰度值计算，计算公式为：Grey=0.299*Red+0.587*Green+0.114*Blue， 其中，像素点的颜色由RGB(Red,Green,Blue)组成，Red表示红色，Green表示绿色，Blue表示蓝色。之后进行灰度值编码和种群初始化，对灰度值进行二进制编码，形成R个初始种群；计算适应度，将Otsu算法的最大类间方差公式作为适应度函数，类间方差越大，适应度就越好；

(3)迭代操作，进行多次迭代，直到适应度满足设定的最优解；选择最佳阈值，将满足最优解的两个个体作为最优个体，解码两个最优个体求出对应的灰度值K1和K2(K1

(4)阈值分割，经过NDWI水体指数处理过的图像已经能够强化一部分水体特征，因此，灰度值小于K1或者灰度值大于K2的像素点可以直接进行阈值分割。但是灰度值的K1和K2之间像素点可能会存在阴影和水体混杂在一起的情况，此时直接进行阈值分割无法准确的将水体和阴影区分开，因此，引入滑动窗口和自适应阈值算法进行判断；

(5)设置M×M窗口的大小(本文设置M=3，即判断该像素点周围的8个像素点)，若窗口内的像素点均为水体像素点或非水体像素点，则将该像素点也判断为水体像素点或非水体像素点。其它情况则保持像素点的灰度值不变，采用自适应阈值算法进行局部阈值判断。自适应阈值算法通过计算邻域内的高斯滤波确定阈值，可以根据图像不同区域的特征进行局部的阈值计算，更加适用于水体边界的提取。

遥感图像的水体指数直方图大多呈现双峰分布，因此使用遗传算法和最大类间方差公式计算出的两个阈值能够更好体现出图像中的水体特征，之后使用滑动窗口遍历图像，进一步区分水体和阴影，能够很好改善误提取和漏提取现象。最后使用自适应阈值算法进行局部阈值判断，可以使水体边界更加明显。经过改进后的Otsu算法，考虑到了遥感影像中水体的光谱特征，更加适用于遥感图像的水体信息提取，能够更好区分水体和其它地物。经过改进后的Otsu算法，考虑到了遥感图像中水体的光谱特征，更加适用于遥感图像的水体信息提取，能够更好区分水体和其它地物。

2.2 水库水体面积提取方法

在一幅遥感影像中会包含很多的水体信息，如河流、塘坝等，如果要对水库水体进行面积提取，即使对水库部分进行截取，也会有很多其它的水体存在，为了能够更加方便地研究水库水体信息的变化，并且计算水库水体面积提取的精度，本文采用最大连通域算法进行水库的水体提取。最大连通域算法是对二值图像中的连通域进行标记的算法，将图像中的水体信息作为目标像素，把相邻像素中具有相同像素值的像素点作为一个连通区域，每找到一个连通区域就对此连通区域赋予唯一的标识，最后对比所有连通区域中包含的像素个数，找出最大的连通区域，即为水库的水体信息。最大连通域算法处理图像的具体过程如图1所示，算法具体步骤为：

(1)输入二值图像F={fij},fij表示像素点(i,j)的灰度值，设置初始COMP=1，COMP表示连通域的唯一标识，初始COMP=1表示图像有一个初始连通域即图像本身。假设目标像素的灰度值为1，且初始图像仅由灰度值为0或1的像素点组成；

(2)对图像中的像素点逐一扫描，当fij=1时，令COMP+1并令fij=COMP，对像素点(i,j)的8邻域进行顺时针扫描，若有非0像素，则令该像素点的灰度值等于COMP，再以此像素点为中心进行8邻域扫描，即进行深度优先遍历。当一次深度优先遍历结束时，该连通域的所有像素点的灰度值均为此时的COMP，即第COMP-1个连通域；当fij=0时，则跳过该像素点，继续扫描下一像素点；

(3)当一个连通域遍历结束后，回到初始的像素点(i,j)的下一个像素点即(i,j+1)，继续重复第(2)步的操作，直到所有的像素点都被遍历到，即fij>1或fij=0；

(4)当找到所有的连通域后，计算每个连通域中像素点的个数，即为连通域的大小，找到最大的连通域并进行标记，将非最大连通域内的像素点的灰度值置0。

图1 最大连通域算法处理图像的具体过程

通过遍历图像中所有的像素点，将每一个水体区域都进行连通域的标识，最后找到连通域最大的水体区域，即为需要的目标水库的水体区域。

2.3 水库面积计算方法

图像水库水体面积的计算方法主要是利用图像中水体的像元数量和图像的空间分辨率得到，首先对提取水体后的图像进行二值化处理，采用自适应阈值算法计算图像阈值，若像素点的灰度值大于阈值，则将其设为255，若小于阈值，则将其设为0，即使得图像中像素点的灰度值仅包含0和255，其中0表示非水体像素点，255表示水体像素点。对图像进行遍历，统计水体像素点的个数，记为t，则图像中水体的面积计算公式为

s=t·r2

(7)

其中，s表示水库水体面积，t为图像中水体像元的个数，r表示所用遥感图像的空间分辨率。

3 实验过程与结果分析

3.1 数据获取与数据预处理

实验数据使用Landsat-8OLI遥感影像，Landsat-8是Landsat系列的第八颗卫星，其携带的OLI陆地成像仪含有9个波段，能获取丰富的遥感图像信息。其中Band 5 NIR(近红外波段)的波长范围为0.845 μm～0.885 μm，Band 4 Red(红波段)，波长范围为0.630 μm～0.680 μm，Band 3 Green(绿波段)波长范围为0.525 μm～0.600 μm，本文使用ENVI5.3软件将这3个波段进行融合，得到标准假彩色图像。标准假彩色图像能够更好区分出植被和水体，更加有利于水体信息的提取。之后使用ENVI5.3软件对图像进行辐射定标和大气校正处理，辐射定标用来将遥感影像的DN值转换为反射率，大气校正是为了消除或减小大气分子和气溶胶的散射与吸收对地物反射率的影响。使用NDWI水体指数对图像进行初始水体提取，水体指数处理后的图像水体与地物的区分更加明显，更加有利于后续算法的处理。初始遥感图像处理的基本流程如图2所示。

图2 遥感图像处理过程

本文共进行两组实验，分别选取石梁河水库和小塔山水库2018年8月的遥感影像，石梁河水库总库容5.31亿m3，最大库面积92 km2，是江苏省内最大的人工水库，小塔山水库始建于1958年，库容量3亿m3，兴利水位时的水域面积约为24 km2，是苏北地区第二大人工湖。将遥感影像中包含的水库的部分单独截取出来，进行水体提取，图像中包含水库的全貌和其它的一些地物信息，这样进行水体提取时能够更好获得水库的水体信息。

3.2 实验结果分析

对实验数据分别使用NDWI水体指数法、谱间关系法、直方图法、Otsu算法和基于遗传算法的改进Otsu算法(即本文方法)进行水体提取，第一组实验使用石梁河水库的遥感影像，水体提取结果对比如图3所示，图3(a)为使用ENVI软件处理后的原始图像，图3(b)为使用NDWI水体指数法提取水体的结果，从图中可看出，NDWI水体指数对大范围的水体提取效果较好，但是细小的水体提取不够准确。图3(c)为谱间关系法水体提取结果，图中黑色部分为水体，此方法提取的水体边界信息较为模糊，不能很好突出边界信息，且结果中存在大量的噪声。图3(d)为直方图法的水体提取结果，此方法在处理细小水体和水体阴影时效果不好。图3(e)为Otsu算法水体提取结果，从图中可以明显看出，Otsu算法存在误提取现象，无法将阴影和水体准确的区分开来；图3(f)为本文方法水体提取结果，从图中可以看出，本文方法不仅可以准确提取出水体信息，在水体边缘的提取上要比其它两种方法准确，阴影和水体的区分也更加明确。从图3中可看出，本文方法在水体提取上能够取得良好的效果，在提取细节上要优于其它几种方法的提取效果，水体边界的处理也更加精确，误提取的现象也得到了明显的改善。

图3 石梁河水库水体提取效果对比

第二组实验采用小塔山水库2018年8月的遥感影像，提取结果对比如图4所示，图4(a)为小塔山水库使用ENVI软件处理后的原始遥感影像，相比于石梁河水库的水体，小塔山水库的水体与其它地物之间的光谱差异更小，提取难度也较大。图4(b)为小塔山水库NDWI水体指数的提取结果，从图中可看出，NDWI水体指数法存在着阴影和水体模糊不清的情况，在水体边界的提取上精度不高；图4(c)为谱间关系法水体提取结果，其中黑色部分为水体，图4(d)为直方图法水体提取结果，从图中可以看出，两种方法都存在着误提取现象，在细小水体和水体阴影的区分上不够精确。图4(e)为Otsu算法的提取结果，从图中可看出，Otsu算法在水体边界的提取上比NDWI水体指数法要更准确更清晰，但是却存在着一定的误提取现象，将一些光谱特征近似水体的阴影部分也当作水体提取出来；图4(f)为本文方法提取结果，从图中可看出，在水体边界的提取上比其它几种方法都要好，能够区分出光谱特征近似于水体的其它地物，提取的准确度好于其它几个算法。

通过以上两组实验可以看出本文提出的水体提取方法能够满足不同光谱特征的遥感影像，能够适用于具有不同特征的水库，且提取的效果较好，可以满足水库的面积监测、水体信息管理等方面的需求。

图4 小塔山水库水体提取效果对比

对两组实验中本文方法提取的水体结果图进行水库水体面积提取，使用本文2.2中的水库水体面积提取方法提取水库水体信息，使用最大连通域算法找出图像中最大的连通域，此连通域即为水库的水体面积，提取结果如图5所示，从图中可以看出，使用此方法进行水库水体面积提取的效果良好，能够很好的将非水库水体去除，清晰展现出水库的水体边界和水体面积信息，可以很方便的对水库的水体面积和边界的变化进行比较分析。当需要进一步计算水库的水体面积时，可直接对水库的水体像元数量进行统计。

图5 水库水体提取结果

对第一组实验石梁河水库进行水库水体的面积计算，根据选用的遥感影像的成像时间，查阅石梁河水库的水位数据，得到所用遥感图像成像时石梁河水库的平均水位为23.97 m，再根据石梁河水库的库容曲线和对应水库面积的关系，计算出水库的水体面积约为45.61 km2，将此面积作为水库的实测面积，进行实验对比和误差计算。

图像水库水体面积的计算方法主要是利用图像中水体的像元数量乘以图像分辨率，Landsat-8的遥感影像分辨率为30 m，则一个像素点所表示的面积为900 m2，统计出图像中水体像元的个数，就可以得出图像中的水库水体面积，其结果对比见表1。从表1中可以看出，5种方法算得的面积的绝对误差都在2 km2以内，NDWI水体指数方法所提取的水库面积的相对误差为3.19%，谱间关系法和直方图法所提取的水库水体面积的相对误差分别为1.2%和3.65%，Otsu算法提取的水库水体面积的相对误差为1.4%，本文方法提取的水库水体面积的相对误差仅为0.39%，相较于前几种方法水库水体面积计算的精度有明显的提升。

表1 石梁河水库水体面积计算结果对比

对第二组实验小塔山水库进行水库水体面积计算，小塔山水库的实测面积通过使用原始遥感影像进行目视解译获取，使用ENVI软件获得水库水体的面积约为21.69 km2，将此面积值作为水库的实测面积。5种方法得出的面积值见表2，从表2中可看出，5种方法计算出的面积绝对误差都在2 km2以内，能够满足实际的需要，NDWI水体指数法计算出的面积值得相对误差为5.98%，谱间关系法和直方图法所提取的水库水体面积的相对误差分别为8.23%和7.46%，Otsu算法计算出的面积值相对误差为3.12%，本文方法计算出的面积值相对误差仅为0.36%，因此，本文方法在计算面积的准确度上比其它几种方法有较明显的提升，能够满足实际应用的需求。

表2 小塔山水库水体面积计算结果对比

4 结束语

本文通过对遥感影像的光谱特征进行分析，提出一种基于遗传算法和改进Otsu算法的水体提取方法，并使用最大连通域算法提取水库的水体信息。利用ENVI软件对遥感影像进行波段融合和预处理，使用NDWI水体指数对图像进行初始水体提取，然后用遗传算法优化改进Otsu算法，首先利用遗传算法根据最大类间方差公式计算出两个最佳阈值，之后根据最佳阈值采用滑动窗口遍历图像像素点进行阈值判断，再利用自适应阈值算法进行局部阈值分割，进一步细化水体信息，最后采用最大连通域算法进一步提取水库水体信息，根据计算出的水库面积判断水库水体提取精度。本文方法在水体边界提取和去除阴影方面相比于水体指数法、谱间关系法、直方图法和传统的Otsu算法有很大的改善，能够在很大程度上反应出水库的面积变化和边界特征，并且无需进行人工标注，在提取精度和提取速度上都有很好的表现。