鄭云柯，王俊霖，張世濤

（昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明 650093）

0 引言

水是人類生存必不可少的資源，是大自然賦予人類最寶貴的財富［1］。隨著全球氣候變化和人類活動的不斷增加，地球上水體的分布和特征發(fā)生了顯著變化，給人類社會和自然生態(tài)帶來了嚴重挑戰(zhàn)，水資源的可持續(xù)性問題引起人們越來越多的關(guān)注。因此，及時、準確、全面地監(jiān)測和分析水體信息對于理解水循環(huán)過程、評估水資源狀況、保障水資源的可持續(xù)開發(fā)利用至關(guān)重要。衛(wèi)星遙感技術(shù)因其能夠高效、準確地獲取全球范圍內(nèi)的水體信息而備受關(guān)注。利用衛(wèi)星遙感技術(shù)，可以獲取大量水體信息，從多源、多時相、多尺度的遙感影像中提取和反演水體信息，包括水體范圍、深度、透明度、色度、溫度等參數(shù)，并結(jié)合地面觀測數(shù)據(jù)和數(shù)值模型進行驗證與分析，可以實現(xiàn)對大范圍、復(fù)雜區(qū)域內(nèi)水體變化情況的動態(tài)監(jiān)測和定量評估，為制定合理、有效的水資源開發(fā)利用和管理措施提供科學(xué)依據(jù)［2］。如今，快速、準確地從衛(wèi)星遙感影像上提取水體信息已成為水資源調(diào)查、水資源宏觀監(jiān)測及濕地保護的重要手段。

1 相關(guān)研究

目前，利用遙感技術(shù)監(jiān)測湖泊主要集中在湖泊的水質(zhì)［3］、水位［4］、面積［5］等方面。常見的水體信息提取方法有單波段閾值法［6］、多波段譜間關(guān)系法［7］和水體指數(shù)法等。其中，水體指數(shù)法是一種簡單、快速、高效的水體信息識別方法，如Mcfeeters［8］在1996 年提出歸一化差異水體指數(shù)（Normalized Difference Water Index，NDWI），利用綠光波段和近紅外波段增強了影像中的水體特征；徐涵秋［9］在NDWI 的基礎(chǔ)上，提出改進的歸一化差異水體指數(shù)（Modified Normalized Difference Water Index，MNDWI），采用中紅外波段代替近紅外波段，取得了比NDWI 更好的提取效果；閆霈等［10］提出增強型水體指數(shù)（Enhanced Water Index，EWI），相較于NDWI 和MNDWI，能更有效地區(qū)分河道與背景噪音。針對不同研究區(qū)的水體特性，以上水體指數(shù)法都取得了較好效果。

湖泊作為具有一定面積的相對封閉水域，為了獲取湖泊的矢量邊界，常用的方法包括手工數(shù)字化方法、區(qū)域分割法［11］與邊緣檢測法［12］。當前國內(nèi)外學(xué)者對水邊界提取的研究主要集中在海岸線與湖岸線，如Alesheik 等［13］利用閾值分割法對烏爾米耶湖岸線進行提?。焕钚忝返龋?4］利用Canny 算子對渤海灣海岸帶進行提取以監(jiān)測其時空變化；Karantzalos 等［15］基于遙感影像，分別用Laplacian 算子和 Canny 算子提取海岸線；申家雙等［16］在分析了邊緣檢測算法用于影像水邊界提取的優(yōu)缺點后，提出將Canny 算子與GAC 模型結(jié)合提取影像水邊界的方法；Abolhassani 等［17］利用邊緣檢測算法提取了美國東海岸線；魏東嵐等［18］在MATLAB 平臺上利用邊緣檢測算法、閾值分割法及小波變換法提取海岸線，結(jié)果表明邊緣檢測法中的Canny 算子能準確提取出圍填區(qū)域的海岸線。

已有研究主要使用Canny 邊緣檢測算子提取海岸線，且Canny 算子的閾值多是通過大量實驗獲取，這是Canny算子自身的不足：高低閾值需人為設(shè)定，容易導(dǎo)致檢測中出現(xiàn)大量虛假邊緣［19］。因此，本文在已有研究基礎(chǔ)上，以湖泊為研究對象，將Canny 邊緣檢測算子與遙感水體指數(shù)法相結(jié)合以提取湖岸邊界。針對Canny 算子的不足，利用類間方差最大閾值分割法來計算其閾值，使Canny 算子具有較理想的閾值，從而實現(xiàn)自動化提取湖泊水邊界。

2 研究區(qū)概況

洱海位于云南省大理市郊區(qū)，地跨大理市和洱源縣。洱海是云貴高原九大湖泊之一，也是云南省第二大淡水湖，海拔1 980 m，經(jīng)緯度為100°05'-100°18'E，25°36'-25°58'N。洱海北起洱源縣，南至大理市下關(guān)鎮(zhèn)，湖泊呈東西窄、南北長的條帶狀，因湖泊形狀像耳朵而得名。洱海水域面積250 km2，平均水深10.5 m，平均水位1 974 m，湖水主要靠河流補給。

3 數(shù)據(jù)來源及預(yù)處理

本文采用的遙感數(shù)據(jù)是從USGS 網(wǎng)站上下載的大理市洱海的Landsat-8 OIL 影像，軌道號/行號為131/42，選取影像的云量少于3%，遙感影像獲取時間為2019 年11 月29號，處于湖泊的枯水期。Landsat-8 OIL數(shù)據(jù)共有9個波段，除全色波段空間分辨率為15 m外，其余波段都為30 m。另外，為驗證水體邊界提取精度，從中國資源衛(wèi)星應(yīng)用中心獲取了一幅時期相近的大理市洱海的2 m 全色/8 m 多光譜GF-6號影像，遙感影像獲取時間為2019年12月7號。

對Landsat-8 OIL 影像進行輻射定標、大氣校正、圖像裁剪等一系列預(yù)處理，對GF-6 號影像進行2 m 全色波段和8 m 多光譜融合，形成空間分辨率為2 m 的融合影像。

4 研究方法

4.1 水體信息增強

國內(nèi)外學(xué)者根據(jù)水體在藍綠波段吸收率較低、近紅外波段吸收率較高這一光譜特征，開展了大量基于遙感影像的表面水體自動提取算法研究。本文將利用不同的影像波段，分別采用歸一化差分水體指數(shù)法、改進的歸一化差分水體指數(shù)法和增強型水體指數(shù)法來增強洱海表面水體信息，如表1所示。

Table 1 Model of the water body index表1 水體指數(shù)模型

4.2 水體邊緣檢測算法

傳統(tǒng)的邊緣檢測算法有Roberts 算子、Prewitt 算子、Sobel 算子等，這些算法簡單，雖容易實現(xiàn)，但處理噪聲的能力較差，裂紋邊緣識別不完整，還容易出現(xiàn)偽邊緣現(xiàn)象。能否選擇合適的邊緣檢測算子，將直接影響到結(jié)果的精度。筆者對前人的研究結(jié)果進行比較，發(fā)現(xiàn)使用Canny 算子處理邊緣和噪聲的效果最佳［20］。因此，本文選用Canny算子對水體增強的遙感影像進行邊緣檢測，主要分為以下4個步驟：

（1）利用高斯函數(shù)對圖像進行平滑處理。高斯函數(shù)表達式為：

平滑圖像后的 結(jié)果I(x，y) 表達式為：I(x，y)=G(x，y)*f(x，y)，其中G(x，y)為高斯函數(shù)，f(x，y)為原圖像。

（2）對平滑后的圖像進行梯度幅值和方向的計算?；?×2 模板，通過X、Y方向像素的一階導(dǎo)數(shù)來確定梯度幅值。設(shè)Fx(x，y)為X方向的偏導(dǎo)數(shù)，F(xiàn)y(x，y)為Y方向的偏導(dǎo)數(shù)，則梯度幅值M(x，y)和方向θ的表達式為：

（3）對梯度幅值的非極大值進行抑制處理，判斷梯度幅值在其八鄰域內(nèi)是否為最大值，若為最大值則為邊緣，否則置零。

（4）對梯度選取兩個門限閾值，即高門限閾值TH和低門限閾值TL，兩者通常的關(guān)系為［21］：

取出經(jīng)過非極大值抑制后圖像中的最大梯度幅值，標記大于高門限閾值的點（邊緣點），將小于低門限閾值的點置0，即可提取出完整邊緣。

4.3 Canny算子閾值設(shè)定

由Canny 算法步驟可知，門限閾值的選取是進行圖像邊緣提取的關(guān)鍵，而Canny 算子需要人為預(yù)先設(shè)定高低閾值，因此需要很多次反復(fù)實驗才能找到合適的閾值。為解決閾值需要人為預(yù)先設(shè)定的問題，本文利用類間方差最大閾值分割算法來確定Canny算子的門限閾值。

類間方差最大閾值分割算法別稱大津（Otsu）法［22］，是由日本學(xué)者大津在1979 年提出的一種自適應(yīng)閾值確定法。其數(shù)學(xué)描述為對于給定的圖像L(x，y)，設(shè)圖像的灰度級為M，灰度值h像素點的像素數(shù)為n，則其出現(xiàn)的概率為［23］：

采用閾值t為界限把圖像區(qū)分為兩個區(qū)域，假設(shè)背景像素區(qū)域為A1，前景像素區(qū)域為A2，則兩個區(qū)域的總概率為：

設(shè)圖像整體的灰度均值為η，則前景像素區(qū)域A1與背景像素區(qū)域A2的灰度均值為：

即當σ2為最大值時，對應(yīng)的t為選擇的最佳閾值。

根據(jù)上述Otsu 算法原理確定一個分割閾值t，將選取的分割閾值t作為Canny 算子的高閾值TH，再用此高閾值乘以一個比例因子0.5 作為其低閾值TL［24-25］。得到閾值的表達式如下：

由此Canny 算子的高低閾值選擇只與分割閾值t有關(guān)，不需要再人為多次反復(fù)實驗設(shè)定，這就為Canny 算子中閾值的確定問題提供了一個較好的解決方法，也增強了Canny算子的分割能力和自適應(yīng)性。

總體技術(shù)路線如圖1所示。

Fig.1 Technical flow圖1 技術(shù)流程

5 結(jié)果與分析

5.1 水體邊界自動提取

5.1.1 水體灰度差增強

如圖2 所示，NDWI、MNDWI 和EWI 水體指數(shù)法都能有效拉伸水體和非水體之間的灰度差異。從目視解譯的效果來看，3 種方法拉伸灰度差異的效果各不相同。相比于MNDWI 法和EWI 法，NDWI 法對洱海南邊的出水口即西洱河的提取有所欠缺。

Fig.2 Water body enhancement effect圖2 水體增強效果

5.1.2 邊緣矢量化

先對3 幅水體增強影像進行歸一化處理，由大津法的數(shù)學(xué)原理可知其特性［26-27］：當目標地物與背景地物的面積比例懸殊時，灰度直方圖可能會無明顯“雙峰”或呈“多峰”形態(tài)，此時使用大津法的分割效果不佳。為得到最優(yōu)閾值，需滿足灰度頻率直方圖呈“雙峰”的情況［28］。通過圖2 可看出目標地物水體與背景地物非水體的面積比例相差不大，研究發(fā)現(xiàn)3 種水體指數(shù)法的灰度頻率直方圖均呈“雙峰”分布，如圖3 所示。其中，橫坐標Data Value 代表歸一化后的灰度像素值。

Fig.3 Frequency histogram of water index method圖3 水體指數(shù)法頻率直方圖

然后，使用Ostu 算法得到水體增強影像的分割閾值（保留兩位有效數(shù)字），最后利用關(guān)系式（7）得到Canny 算子的高閾值TH 和低閾值TL，結(jié)果如表2所示。

Table 2 Thresholds of various water index methods表2 各類水體指數(shù)法閾值

將提取好的柵格影像結(jié)果轉(zhuǎn)換為矢量數(shù)據(jù)，再與原始遙感影像疊加顯示。如圖4 所示，整體上，水體邊界提取的邊緣連續(xù)性及去噪效果都達到最佳。3 種水體指數(shù)法均可提取出洱海的水體邊界，對于洱海的北部、西部、東部的非目標地物類，也均提取了其邊緣。就北部的目視效果而言，MNDWI+Canny 算子提取的非目標地物略少一點。相比其它兩種水體指數(shù)法，NDWI+Canny 算子還提取出了位于洱海西南部的蒼山積雪邊緣及南部的部分建筑邊緣。因此，選取合適的邊緣檢測閾值，可減少目標地物的非必要干擾信息，為后續(xù)的處理節(jié)省時間。

Fig.4 Vector boundary overlay image of water body圖4 水體矢量邊界疊加影像

5.1.3 邊緣細部提取效果

選取洱海湖岸的濕地、河灘及洱海出水口3 種典型地區(qū)進行分析，濕地是水體、植被與裸地等按不同比例混合組成的一種土地形式。洱海湖岸周圍存在著不少濕地區(qū)，比較著名的有海舌濕地公園、洱海月濕地公園等。以海舌濕地為例，分析其濕地邊界提取效果。如圖5 所示，運用NDWI+Canny 算子提取出的海舌濕地邊緣出現(xiàn)明顯的不連續(xù)現(xiàn)象，且產(chǎn)生了細微的偽邊緣現(xiàn)象。其將濕地的近水側(cè)岸線和近陸側(cè)岸線均提取了出來，整體上不能很好地提取出完整岸線。運用MNDWI+Canny 算子提取出的海舌濕地邊緣連續(xù)、完好，可以看出其將湖岸濕地劃分給了水域。運用EWI+Canny 算子提取出的海舌邊緣同樣出現(xiàn)了不連續(xù)現(xiàn)象及少量斷點，位置不夠準確，可以看出其將湖岸濕地劃分給了陸地。

Fig.5 Sea Tongue wetland boundary圖5 海舌濕地邊界

河灘是由于泥沙沉積而形成的天然灘涂土地，選取洱海湖岸的某處河灘進行對比分析。如圖6 所示，3 種方法所提取的邊緣都是連續(xù)的。NDWI+Canny 算子在大河灘（右下方）處和小河灘（右上方）處均有向外膨脹的現(xiàn)象。MNDWI+Canny 算子在提取河灘邊界時有較好效果，可以看出其將河灘也劃分給了水域。EWI+Canny 算子則在大河灘處多提取了一處偽邊緣，在小河灘處出現(xiàn)了向外膨脹的現(xiàn)象。

Fig.6 River beach boundary圖6 河灘邊界

位于洱海西南方向的湖岸出水口稱為西洱河，洱海與西洱河以大理市的興盛大橋為界，該地附近還有著名的洱海月濕地公園，選取該區(qū)域進行分析。如圖7 所示，3 種方法都能有效、連續(xù)地提出洱海出水口邊界，但NDWI+Canny算子和EWI+Canny 算子多提取了大橋邊緣及西洱河（橋梁左側(cè)）的部分邊緣，且兩種方法多提取的部分幾乎是一致的。MNDWI+Canny 算子則恰好提取了出水口邊界。宏觀洱海月濕地的提取效果（位于出水口的左上方），3 種方法并無明顯差異，均能有效、連續(xù)地提取出邊緣；微觀洱海月濕地的提取效果，NDWI+Canny 算子和EWI+Canny 算子都出現(xiàn)了向外膨脹的現(xiàn)象，而MNDWI+Canny 算子提取的邊緣效果較好。

5.1.4 邊界整體提取結(jié)果

由于影像分辨率的問題，混合像元的存在使Canny 算子無法完全克服噪聲，導(dǎo)致提取的邊界存在微小的不連續(xù)情況。疊加原始影像，將存在的微小邊緣斷裂進行連接，并剔除非目標地物的邊緣信息，輸出3 種水體指數(shù)法的洱海水域矢量圖，如圖8所示。

Fig.8 Automatic extraction results of water boundary vector data圖8 水體邊界矢量數(shù)據(jù)自動提取結(jié)果

5.2 水域邊界精度評價

5.2.1 目視判別

將提取的洱海邊界與時期相近的GF-6 號2m 融合影像進行疊加，如圖9 所示。目視檢查發(fā)現(xiàn)NDWI、MNDWI和EWI 都能有效提取出邊緣，MNDWI 更接近精確的邊緣位置，NDWI和EWI都略微有向外膨脹的現(xiàn)象。

Fig.9 Gf-6 2m fusion image of superimposed water boundary vector data圖9 疊加水體邊界矢量數(shù)據(jù)的GF-6號2m融合圖像

5.2.2 空間統(tǒng)計分析

基于空間緩沖區(qū)統(tǒng)計法對提取精度進行量化分析，從統(tǒng)計學(xué)角度提出3 點準則：①向水緩沖區(qū)B 的最大、最小及平均值都要小于向陸緩沖區(qū)A；②B 內(nèi)標準差要小于A；③兩個緩沖區(qū)的標準差Δσ=σA-σB要盡可能大［29］。

因水體對近紅外波段的吸收率較強，故這里選取Landsat-8 OIL 近紅外波段影像作為統(tǒng)計的基礎(chǔ)影像。以水邊界為基準分別在水陸兩側(cè)建立緩沖區(qū)A 和B，緩沖范圍為2 個像素值，即60 m。再統(tǒng)計緩沖區(qū)域的像素灰度特征（DN）值，如表3所示。

Table 3 Characteristic statistics of pixel DN value in buffer表3 緩沖區(qū)內(nèi)像素DN值特征統(tǒng)計

基于表3 的特征值分析，洱海水體邊界提取位置較準確的是MNDWI+Canny 算子。

6 結(jié)語

本文選用2019 年的OIL 影像作為研究數(shù)據(jù)，以洱海為研究對象，分別運用歸一化差分水體指數(shù)、改進的歸一化差分水體指數(shù)和增強型水體指數(shù)對遙感影像進行水體信息灰度增強。分析3 種水體指數(shù)法的灰度頻率直方圖，使用大津法得到水體灰度影像的分割閾值。將分割閾值作為Canny 算子的高門限閾值，結(jié)合Canny 算子高低門限的閾值關(guān)系式，得到低門限閾值以自動提取洱海水體邊緣信息，進而疊加原始影像進行典型地區(qū)分析。將提取的水邊界疊加在GF-6 號2 m 分辨率的融合影像上進行目視判別，再經(jīng)空間緩沖區(qū)分析得到：在NDWI、MNDWI 與EWI 3種水體指數(shù)法中，MNDWI 與Canny 算子結(jié)合的提取效果最好。

另外，本文中存在些許邊緣不連續(xù)的情況，這是由于水體指數(shù)法大多采用具有中紅外波段的遙感影像，混合像元的存在對提取結(jié)果仍具有一定干擾性。Canny 算子無法完全克服噪聲影響，今后可將混合像元分解或?qū)⑾伻核惴ǖ燃夹g(shù)應(yīng)用于研究中。