遥感图像污染监测分析的Grab-Cut算法实现

宋占杰 ，汤 杨，王兵振 ，张东亮，田 震

（1．天津大学 数学学院，天津 300354；2．天津大学 电视与图像信息研究所，天津 300072；3．国家海洋技术中心，天津 300112）

渤海是中国半封闭内海，自净能力差，周边又都是重工业发达城市，目前海水污染非常严重。据不完全统计，每年排入渤海的污染物达70多万吨，占全国入海污染物总量的47．7%。渤海沿岸直排排污口携带大量悬浮颗粒物、营养盐、重金属入海，渤海的海洋生态环境毁坏严重、压力巨大。

天津紧邻渤海，并且位于渤海的核心区，与渤海其它海域相比，天津周边海域海洋自净能力更加脆弱。同时，天津又是重工业较为发达的城市，每年向渤海排放大量污水，使得天津周边的海洋污染日趋严重。图1是2014年渤海监测的陆源入海排污口状况。图中监测点的颜色越红，污染越严重。图中显示，天津市周边污染十分严重的排污口主要有3个，对渤海湾造成了严重的污染。

图1 2014年渤海监测的陆源入海排污口状况

传统的入海排污口监测考察主要依靠船只人工采样和样本分析，这种方法耗时费力，受环境的制约较大。随着遥感技术的快速发展，以及大量高精度遥感卫星的升空，高分辨率、多频段、价格低廉的遥感数据的获取成为现实。这些遥感数据的使用为我们更好地对天津周边的海洋污染进行追踪提供了便利，且一些图像分割算法能够帮助我们对污染区域进行监测和量化分析。

图像分割是把图像分割成一些连续的、非重叠的区域。这种迭代的图像分割算法目前在计算机视觉方面应用越来越广泛[1-3]。Boykov[4]第一个在图像分割领域中提出马尔科夫随机场模型，基于这个模型，提出了Graph-Cut图像分割算法。只需在背景和前景中标记简单的几笔，此算法能够在短时间内分割出想要的图形。在Boykov和Jolly研究的基础上，又出现了大量基于Graph-Cut和用户笔画的算法模型[4-5]。Slabaugh[6]提出了一种椭圆形状先验模型，此算法主要应用于对医学中共振图像的血管和淋巴结的提取。Rother[7]提出了Grab-Cut算法，此算法的迭代速度更快，而且所需的用户笔画更少，适用范围广泛，能够应用在很多领域。本文主要进一步利用改进的Grab-Cut分割算法对天津沿岸的排污口进行监测。

本研究利用美国USGS Landsat-8卫星的遥感图像，结合改进的Grab-Cut图像分割算法，对天津周边渤海湾的排污口进行了监测，并对污染区域进行了时间和空间上的分析，以期能对相关部门关于河流污染的治理提供参考。

1 Grab-Cut算法

1.1 Graph Cut算法

对于每一幅遥感图像z=(z1,…,zn,…zN)，zn表示图像中的每一个像素点的灰度值。定义α=(α1,…,αn,…αN)，其中 αn∈{0,1}，0 代表像素点属于背景，1 代表像素点属于要分离出的前景。本文中，在算法进行前会要求使用者在图形中初步选定一个矩形框，矩形框内αn=1，矩形框外αn=0。

针对αn=1的区域和αn=0的区域，即前景与背景，分别建立像素灰度值的频率分布直方图：

前景和背景的频率分布直方图中，总会出现一个或两个所占比例较大的灰度值，这个灰度值就是前景和背景中的主灰度值，前景和背景的区分依赖于这个主灰度值。

利用每个图像像素点的灰度值向量z=(z1,…,zn,…zN)，和区分前景和背景的向量 α=(α1,…,αn,…αN)，以及前景和背景的灰度频率分布直方图θ，可以建立最小能量函数：

最后优化出的向量αˆ就能够较好地分离出前景和背景。

其中，E(α,θ,z)被定义为 Gibbs能量：

上式U用来衡量α的正确度，表达式为：

式中：[…]代表一种判断，满足括号内的条件为1，不满足括号内的条件为0。dis(m,n)表示两点之间的欧几里德距离。依据 Boykov[4]，β=(2＜(zm-zn)2＞)-1，＜…＞表示均值。(m,n)∈C中的C表示像素点与邻像素点配对组成的集合，对于某个像素点而言，其邻像素点为其前后左右4个以及斜对角的4个共8个像素点，对于一个具体的像素点而言，C中的元素共有8个。基于大量的训练集测试，γ的值选为50较为合适(见文献[9])。

1.2 Grab-Cut算法模型的建立

传统的Graph Cut算法不需要使用像素点的RGB信息，它基于像素点的灰度值直接进行计算。为了更好地使用这些颜色信息，建立了2个高斯混合模型（GMMs）[8-9]，一个对应于背景，一个对应于要分割出的前景。为了更好地应用高斯混合模型，对每个像素 zn赋值 kn，kn∈{1,…,K}，这样，对于一整幅图像，有 k={k1,…,kn,…,kN}。

借鉴上述Graph Cut算法，建立如下新的Grab-Cut算法最小能量函数：

V用来衡量整个分割的稳定性，表达式为：

E（α,k,θ,z）定义为如下的 Gibbs能量：

其中 α=(α1,…,αn,…αN)，αn∈{0,1}，0 代表背景，1代表想要分离出的前景。类似于Graph Cut算法，式中U衡量α的正确度，V衡量稳定性。分别定义为如下形式：

其中 p(…)代表高斯分布，π(…)为 αn与 kn的混合加权系数。

参数 γ=50，β=(2＜‖zm-zn‖2＞)-1，其中的＜…＞代表均值；C代表某一像素点与其邻像素点组合形成的集合[7]。

Grab-Cut算法会迭代计算去优化能量函数，以求得最优解，算法会对公式(6)中的向量α自动逼近，使分割的边缘部分清晰准确。

1.3 进一步改进的Grab-Cut算法

Grab-Cut算法在以前经常用来处理一些小尺寸照片，很少用在尺寸很大的高分辨率遥感图像上。

本文利用Grab-Cut算法来处理大尺寸的遥感图像，并用于海洋污染的监测。因为对整个遥感图像进行图像分割的效率低且耗时长，我们基于对污染区域的认识，划定感兴趣的区域，然后再对这些区域进行图像分割，从而产生更好的分割效果。

2 遥感图像的获取

2013年2月11日，NASA成功发射了Landsat-8卫星，卫星上新增两个主要载荷：运营性陆地成像仪OLI和热红外传感器TIRS，2013年5月30日开始向全球提供免费数据下载。两个主要载荷所收集的数据包含11个波段的信息，不同波段的组合可以得到不同的图像，以用于不同的研究。本文的数据只用到其中的红、绿、蓝三个波段。

3 天津周边渤海湾入海排污口的监测

3.1 天津入海排污口的监测

天津市河流流域众多，包括蓟运河、海河、独流减河、北排河等等。天津又是重工业较为发达的城市，大量的重工业企业聚集在河流周边，这些企业排放的污水大部分流入到这些河流当中。每年，这些河流携带大量的工业生活污染物排入渤海，给渤海的生态环境造成了巨大的破坏。

图2是一幅由美国USGS landsat-8卫星于2014年3月21日拍摄的天津市周边渤海湾的遥感图像。从图中可以看出，天津周边渤海湾的污染非常严重，黑色的区域均为河流排向渤海的污染物所形成的异常图。天津周边浅海区域的海藻繁殖较为严重，整体呈现蓝绿色，所以一旦有污水排入会与这些海藻形成鲜明的对比。

天津沿海地区主要分为3个流域，自北向南分别为蓟运河、潮白新河、永定新河流域；海河流域；独流减河、子牙新河、北排河流域。本文对这3个主要流域，利用改进的Grab-Cut算法来分割出各个流域排污口的污染区域，如图3。图中绿色方框表示“感兴趣区域”，红色线条表示“标记的前景”，蓝色线条表示“标记的背景”。

图2 天津周边渤海湾污染情况

图3 改进的Grab-Cut算法分割排污口污染区域

图中的结果显示，算法的分割效果非常好，分割出了污染区域。能够看出各个流域入海口的污染比较严重，污染面积比较大。利用比例尺，估算出了各个污染区域面积的大小。其中，蓟运河、潮白新河、永定新河流域和海河流域的污染最严重，污染面积分别约为259 km2和217 km2。蓟运河、潮白新河、永定新河流域分布区域广阔，流域内的河流源起河北，流经天津市重工业发达的滨海新区，最后流入渤海；海河流经人口较多的南开区、和平区以及东丽区，最后经过重工业发达的滨海新区注入渤海，河流中携带大量的生活污染物和工业污染物。这些原因造成了这两个流域入海口的污染非常严重。独流减河、子牙新河、北排河流域流经近几年发展迅速的津南区，区域中的重工业发展迅速，河流也有较大的污染，但相比于上面两个流域污染相对较轻，总污染面积约为116 km2。

3.2 天津入海排污口污染区域随时间的变化

针对已有的改进算法作为典型实例，本文分析了2014年到2016年这3年，上述排污口污染区域的变化情况。所选的4幅遥感图像的拍摄时间分别为2014年2月21日、2015年3月24日、2016年3月26日以及2016年12月23日。

图4展示了2014-2016年天津渤海湾整体的污染变化情况。图5展示了各个排污口在2014-2016年局部的污染变化情况。

图4 2014-2016年天津渤海湾污染整体变化情况

图5 2014-2016年天津渤海湾各个排污口污染局部变化情况

图4～图5显示，各个排污口的污染区域都有变化。利用比例尺，估算出每个区域的污染面积。各个区域的污染面积随时间的变化情况见表1，折线图见图6。

表1 各排污口面积随时间的变化情况

图6 各排污口面积随时间的变化情况

对比污染面积折线图(图6)以及各个排污口污染区域变化情况图(图5)，蓟运河、潮白新河、永定新河流域（图5左边一列）和独流减河、子牙新河、北排河流域(图5右边一列)排污状况有较大改善。但海河流域（图5中间一列）的排污状况虽然有所减轻，但相比以上两个区域，一直非常严重，这间接说明海河流经的滨海新区重工业比较发达，大量重工业企业污染物和生活污水排向海河流域，导致海河的污染十分严重。另一个原因是天津港位于海河流域下游，港内海水与外界交换很少，自净能力非常差。因此，对海河周边重工业企业的排污治理刻不容缓。

总体而言，天津市对于工业排污的治理有一定成效，近3年排污口对渤海的污染明显减轻，但仍需进一步严加管理。

3.3 其它图像分割算法与“Grab-Cut”算法的比较

本文选择了一些图像分割算法与Grab-Cut算法进行比较，从而分析不同的图像分割算法分割效果的优劣。这里使用“蓟运河、潮白新河、永定新河流域入海排污口”图像，分别使用PhotoShop中的“Quick Search”工具和“Magnetic Lasso”工具，以及“Graph Cut算法”这3种图像分割算法与“改进的Grab-Cut算法”进行分割效果比对（图7）。

图7 几种图像分割算法的比较

图7所列的方法，“Quick Search”和“Magnetic Lasso”是PhotoShop图像处理软件中较为常用的图像分割算法，“Quick Search”在用户指定前景后，会自动识别与前景相似的区域，但分割效果较差，可能误将大量无关区域识别为要分割的前景；“Magnetic Lasso”可以较好地分割出前景，但其缺点是需要用户比较详细地标记出边界，这种方法繁琐、费时且效率低下。有效的图像分割算法应该是只需要简单的几笔标记就能够完成图像分割；已有的“Graph Cut”算法虽然标记较少，但是在图像的边缘处光滑性差，准确率低，难得到一块完整的分割结果；改进的“Grab-Cut”算法能够较好地区别要分割的背景和前景，在较少的用户标记下，能够很快地分割出所需区域。综合比较，改进的“Grab-Cut”算法的图像分割效果在已知算法中效果最优。

4 结论

对比于传统的“Quick Search”和“Magnetic Lasso”图像分割算法，改进的“Grab-Cut”算法的用户标记要更少；而已有的“Graph Cut”算法虽然标记较少，但图像边缘处光滑性差，准确率低，很难得到一块完整的分割结果。而改进的“Grab-Cut”算法在很少的用户标记下，能够很快地分割出所需区域。综合来看，改进的“Grab-Cut”算法的图像分割算法相较于前几种而言效果最好。我们将改进的“Grab-Cut”算法应用于渤海湾的污染监测。发现近3年天津周围渤海湾的污染非常严重，但近1年来污染情况有所改善。

[1]Kolmogorov V,Zabih R.What energy functions can be minimized via graph cuts?[J].Pattern Analysis&Machine Intelligence IEEE Transactions on,2004,26(2):147-159.

[2]Boykov Y,Veksler O,Zabih R.Fast approximate energy minimization via graph cuts[J].Pattern Analysis&Machine Intelligence IEEE Transactions on,2001,23(11):1222-1239.

[3]Boykov Y,Funka-Lea G.Graph cuts and efficient N-D image segmentation[J].International Journal of Computer Vision,2006,70(2):109-131.

[4]Boykov Y Y,Jolly M P.Interactive graph cuts for optimal boundary®ion segmentation of objects in N-D images[J].2001,1:105-112.

[5]Grady L.Random walks for image segmentation.[J].Pattern Analysis&Machine Intelligence IEEE Transactions on,2006,28(11):1768-1783.

[6]Slabaugh G,Unal G.Graph cuts segmentation using an elliptical shape prior[C]//IEEE International Conference on Image Processing.IEEE,2005:II-1222-5.

[7]Rother B C.Kolmogorov V and Blake A(2004)"GrabCut":interactive foreground extraction using iterated graph cuts[C]//ACM Trans Graph.2010.

[8]Ruzon M A,Tomasi C.Alpha estimation in natural images[J].2000,1:18-25.

[9]Blake A,Rother C,Brown M,et al.Interactive Image Segmentation Using an Adaptive GMMRF Model[M]//Computer Vision-ECCV 2004.Springer Berlin Heidelberg,2004:428-441.