遙感圖像污染監(jiān)測分析的Grab-Cut算法實現(xiàn)

宋占杰 ，湯 楊，王兵振 ，張東亮，田 震

（1．天津大學(xué) 數(shù)學(xué)學(xué)院，天津 300354；2．天津大學(xué) 電視與圖像信息研究所，天津 300072；3．國家海洋技術(shù)中心，天津 300112）

渤海是中國半封閉內(nèi)海，自凈能力差，周邊又都是重工業(yè)發(fā)達城市，目前海水污染非常嚴重。據(jù)不完全統(tǒng)計，每年排入渤海的污染物達70多萬噸，占全國入海污染物總量的47．7%。渤海沿岸直排排污口攜帶大量懸浮顆粒物、營養(yǎng)鹽、重金屬入海，渤海的海洋生態(tài)環(huán)境毀壞嚴重、壓力巨大。

天津緊鄰渤海，并且位于渤海的核心區(qū)，與渤海其它海域相比，天津周邊海域海洋自凈能力更加脆弱。同時，天津又是重工業(yè)較為發(fā)達的城市，每年向渤海排放大量污水，使得天津周邊的海洋污染日趨嚴重。圖1是2014年渤海監(jiān)測的陸源入海排污口狀況。圖中監(jiān)測點的顏色越紅，污染越嚴重。圖中顯示，天津市周邊污染十分嚴重的排污口主要有3個，對渤海灣造成了嚴重的污染。

圖1 2014年渤海監(jiān)測的陸源入海排污口狀況

傳統(tǒng)的入海排污口監(jiān)測考察主要依靠船只人工采樣和樣本分析，這種方法耗時費力，受環(huán)境的制約較大。隨著遙感技術(shù)的快速發(fā)展，以及大量高精度遙感衛(wèi)星的升空，高分辨率、多頻段、價格低廉的遙感數(shù)據(jù)的獲取成為現(xiàn)實。這些遙感數(shù)據(jù)的使用為我們更好地對天津周邊的海洋污染進行追蹤提供了便利，且一些圖像分割算法能夠幫助我們對污染區(qū)域進行監(jiān)測和量化分析。

圖像分割是把圖像分割成一些連續(xù)的、非重疊的區(qū)域。這種迭代的圖像分割算法目前在計算機視覺方面應(yīng)用越來越廣泛[1-3]。Boykov[4]第一個在圖像分割領(lǐng)域中提出馬爾科夫隨機場模型，基于這個模型，提出了Graph-Cut圖像分割算法。只需在背景和前景中標記簡單的幾筆，此算法能夠在短時間內(nèi)分割出想要的圖形。在Boykov和Jolly研究的基礎(chǔ)上，又出現(xiàn)了大量基于Graph-Cut和用戶筆畫的算法模型[4-5]。Slabaugh[6]提出了一種橢圓形狀先驗?zāi)Ｐ?，此算法主要?yīng)用于對醫(yī)學(xué)中共振圖像的血管和淋巴結(jié)的提取。Rother[7]提出了Grab-Cut算法，此算法的迭代速度更快，而且所需的用戶筆畫更少，適用范圍廣泛，能夠應(yīng)用在很多領(lǐng)域。本文主要進一步利用改進的Grab-Cut分割算法對天津沿岸的排污口進行監(jiān)測。

本研究利用美國USGS Landsat-8衛(wèi)星的遙感圖像，結(jié)合改進的Grab-Cut圖像分割算法，對天津周邊渤海灣的排污口進行了監(jiān)測，并對污染區(qū)域進行了時間和空間上的分析，以期能對相關(guān)部門關(guān)于河流污染的治理提供參考。

1 Grab-Cut算法

1.1 Graph Cut算法

對于每一幅遙感圖像z=(z1,…,zn,…zN)，zn表示圖像中的每一個像素點的灰度值。定義α=(α1,…,αn,…αN)，其中 αn∈{0,1}，0 代表像素點屬于背景，1 代表像素點屬于要分離出的前景。本文中，在算法進行前會要求使用者在圖形中初步選定一個矩形框，矩形框內(nèi)αn=1，矩形框外αn=0。

針對αn=1的區(qū)域和αn=0的區(qū)域，即前景與背景，分別建立像素灰度值的頻率分布直方圖：

前景和背景的頻率分布直方圖中，總會出現(xiàn)一個或兩個所占比例較大的灰度值，這個灰度值就是前景和背景中的主灰度值，前景和背景的區(qū)分依賴于這個主灰度值。

利用每個圖像像素點的灰度值向量z=(z1,…,zn,…zN)，和區(qū)分前景和背景的向量 α=(α1,…,αn,…αN)，以及前景和背景的灰度頻率分布直方圖θ，可以建立最小能量函數(shù)：

最后優(yōu)化出的向量α?就能夠較好地分離出前景和背景。

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

上式U用來衡量α的正確度，表達式為：

式中：[…]代表一種判斷，滿足括號內(nèi)的條件為1，不滿足括號內(nèi)的條件為0。dis(m,n)表示兩點之間的歐幾里德距離。依據(jù) Boykov[4]，β=(2＜(zm-zn)2＞)-1，＜…＞表示均值。(m,n)∈C中的C表示像素點與鄰像素點配對組成的集合，對于某個像素點而言，其鄰像素點為其前后左右4個以及斜對角的4個共8個像素點，對于一個具體的像素點而言，C中的元素共有8個。基于大量的訓(xùn)練集測試，γ的值選為50較為合適(見文獻[9])。

1.2 Grab-Cut算法模型的建立

傳統(tǒng)的Graph Cut算法不需要使用像素點的RGB信息，它基于像素點的灰度值直接進行計算。為了更好地使用這些顏色信息，建立了2個高斯混合模型（GMMs）[8-9]，一個對應(yīng)于背景，一個對應(yīng)于要分割出的前景。為了更好地應(yīng)用高斯混合模型，對每個像素 zn賦值 kn，kn∈{1,…,K}，這樣，對于一整幅圖像，有 k={k1,…,kn,…,kN}。

借鑒上述Graph Cut算法，建立如下新的Grab-Cut算法最小能量函數(shù)：

V用來衡量整個分割的穩(wěn)定性，表達式為：

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

其中 α=(α1,…,αn,…αN)，αn∈{0,1}，0 代表背景，1代表想要分離出的前景。類似于Graph Cut算法，式中U衡量α的正確度，V衡量穩(wěn)定性。分別定義為如下形式：

其中 p(…)代表高斯分布，π(…)為 αn與 kn的混合加權(quán)系數(shù)。

參數(shù) γ=50，β=(2＜‖zm-zn‖2＞)-1，其中的＜…＞代表均值；C代表某一像素點與其鄰像素點組合形成的集合[7]。

Grab-Cut算法會迭代計算去優(yōu)化能量函數(shù)，以求得最優(yōu)解，算法會對公式(6)中的向量α自動逼近，使分割的邊緣部分清晰準確。

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

Grab-Cut算法在以前經(jīng)常用來處理一些小尺寸照片，很少用在尺寸很大的高分辨率遙感圖像上。

本文利用Grab-Cut算法來處理大尺寸的遙感圖像，并用于海洋污染的監(jiān)測。因為對整個遙感圖像進行圖像分割的效率低且耗時長，我們基于對污染區(qū)域的認識，劃定感興趣的區(qū)域，然后再對這些區(qū)域進行圖像分割，從而產(chǎn)生更好的分割效果。

2 遙感圖像的獲取

2013年2月11日，NASA成功發(fā)射了Landsat-8衛(wèi)星，衛(wèi)星上新增兩個主要載荷：運營性陸地成像儀OLI和熱紅外傳感器TIRS，2013年5月30日開始向全球提供免費數(shù)據(jù)下載。兩個主要載荷所收集的數(shù)據(jù)包含11個波段的信息，不同波段的組合可以得到不同的圖像，以用于不同的研究。本文的數(shù)據(jù)只用到其中的紅、綠、藍三個波段。

3 天津周邊渤海灣入海排污口的監(jiān)測

3.1 天津入海排污口的監(jiān)測

天津市河流流域眾多，包括薊運河、海河、獨流減河、北排河等等。天津又是重工業(yè)較為發(fā)達的城市，大量的重工業(yè)企業(yè)聚集在河流周邊，這些企業(yè)排放的污水大部分流入到這些河流當中。每年，這些河流攜帶大量的工業(yè)生活污染物排入渤海，給渤海的生態(tài)環(huán)境造成了巨大的破壞。

圖2是一幅由美國USGS landsat-8衛(wèi)星于2014年3月21日拍攝的天津市周邊渤海灣的遙感圖像。從圖中可以看出，天津周邊渤海灣的污染非常嚴重，黑色的區(qū)域均為河流排向渤海的污染物所形成的異常圖。天津周邊淺海區(qū)域的海藻繁殖較為嚴重，整體呈現(xiàn)藍綠色，所以一旦有污水排入會與這些海藻形成鮮明的對比。

天津沿海地區(qū)主要分為3個流域，自北向南分別為薊運河、潮白新河、永定新河流域；海河流域；獨流減河、子牙新河、北排河流域。本文對這3個主要流域，利用改進的Grab-Cut算法來分割出各個流域排污口的污染區(qū)域，如圖3。圖中綠色方框表示“感興趣區(qū)域”，紅色線條表示“標記的前景”，藍色線條表示“標記的背景”。

圖2 天津周邊渤海灣污染情況

圖3 改進的Grab-Cut算法分割排污口污染區(qū)域

圖中的結(jié)果顯示，算法的分割效果非常好，分割出了污染區(qū)域。能夠看出各個流域入海口的污染比較嚴重，污染面積比較大。利用比例尺，估算出了各個污染區(qū)域面積的大小。其中，薊運河、潮白新河、永定新河流域和海河流域的污染最嚴重，污染面積分別約為259 km2和217 km2。薊運河、潮白新河、永定新河流域分布區(qū)域廣闊，流域內(nèi)的河流源起河北，流經(jīng)天津市重工業(yè)發(fā)達的濱海新區(qū)，最后流入渤海；海河流經(jīng)人口較多的南開區(qū)、和平區(qū)以及東麗區(qū)，最后經(jīng)過重工業(yè)發(fā)達的濱海新區(qū)注入渤海，河流中攜帶大量的生活污染物和工業(yè)污染物。這些原因造成了這兩個流域入海口的污染非常嚴重。獨流減河、子牙新河、北排河流域流經(jīng)近幾年發(fā)展迅速的津南區(qū)，區(qū)域中的重工業(yè)發(fā)展迅速，河流也有較大的污染，但相比于上面兩個流域污染相對較輕，總污染面積約為116 km2。

3.2 天津入海排污口污染區(qū)域隨時間的變化

針對已有的改進算法作為典型實例，本文分析了2014年到2016年這3年，上述排污口污染區(qū)域的變化情況。所選的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顯示，各個排污口的污染區(qū)域都有變化。利用比例尺，估算出每個區(qū)域的污染面積。各個區(qū)域的污染面積隨時間的變化情況見表1，折線圖見圖6。

表1 各排污口面積隨時間的變化情況

圖6 各排污口面積隨時間的變化情況

對比污染面積折線圖(圖6)以及各個排污口污染區(qū)域變化情況圖(圖5)，薊運河、潮白新河、永定新河流域（圖5左邊一列）和獨流減河、子牙新河、北排河流域(圖5右邊一列)排污狀況有較大改善。但海河流域（圖5中間一列）的排污狀況雖然有所減輕，但相比以上兩個區(qū)域，一直非常嚴重，這間接說明海河流經(jīng)的濱海新區(qū)重工業(yè)比較發(fā)達，大量重工業(yè)企業(yè)污染物和生活污水排向海河流域，導(dǎo)致海河的污染十分嚴重。另一個原因是天津港位于海河流域下游，港內(nèi)海水與外界交換很少，自凈能力非常差。因此，對海河周邊重工業(yè)企業(yè)的排污治理刻不容緩。

總體而言，天津市對于工業(yè)排污的治理有一定成效，近3年排污口對渤海的污染明顯減輕，但仍需進一步嚴加管理。

3.3 其它圖像分割算法與“Grab-Cut”算法的比較

本文選擇了一些圖像分割算法與Grab-Cut算法進行比較，從而分析不同的圖像分割算法分割效果的優(yōu)劣。這里使用“薊運河、潮白新河、永定新河流域入海排污口”圖像，分別使用PhotoShop中的“Quick Search”工具和“Magnetic Lasso”工具，以及“Graph Cut算法”這3種圖像分割算法與“改進的Grab-Cut算法”進行分割效果比對（圖7）。

圖7 幾種圖像分割算法的比較

圖7所列的方法，“Quick Search”和“Magnetic Lasso”是PhotoShop圖像處理軟件中較為常用的圖像分割算法，“Quick Search”在用戶指定前景后，會自動識別與前景相似的區(qū)域，但分割效果較差，可能誤將大量無關(guān)區(qū)域識別為要分割的前景；“Magnetic Lasso”可以較好地分割出前景，但其缺點是需要用戶比較詳細地標記出邊界，這種方法繁瑣、費時且效率低下。有效的圖像分割算法應(yīng)該是只需要簡單的幾筆標記就能夠完成圖像分割；已有的“Graph Cut”算法雖然標記較少，但是在圖像的邊緣處光滑性差，準確率低，難得到一塊完整的分割結(jié)果；改進的“Grab-Cut”算法能夠較好地區(qū)別要分割的背景和前景，在較少的用戶標記下，能夠很快地分割出所需區(qū)域。綜合比較，改進的“Grab-Cut”算法的圖像分割效果在已知算法中效果最優(yōu)。

4 結(jié)論

對比于傳統(tǒng)的“Quick Search”和“Magnetic Lasso”圖像分割算法，改進的“Grab-Cut”算法的用戶標記要更少；而已有的“Graph Cut”算法雖然標記較少，但圖像邊緣處光滑性差，準確率低，很難得到一塊完整的分割結(jié)果。而改進的“Grab-Cut”算法在很少的用戶標記下，能夠很快地分割出所需區(qū)域。綜合來看，改進的“Grab-Cut”算法的圖像分割算法相較于前幾種而言效果最好。我們將改進的“Grab-Cut”算法應(yīng)用于渤海灣的污染監(jiān)測。發(fā)現(xiàn)近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.