水庫水體面積提取方法

李亞琦，余宇峰

(河海大學 計算機與信息學院，江蘇 南京 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算法有很大的改善，能夠在很大程度上反應出水庫的面積變化和邊界特征，并且無需進行人工標注，在提取精度和提取速度上都有很好的表現。