基于GF-1影像的山區(qū)河道豐、枯水期水體提取方法改進(jìn)

趙程銘，董曉華，薄會(huì)娟，章程焱，張慶玉

（1.三峽大學(xué)水利與環(huán)境學(xué)院，湖北宜昌 443002；2.三峽庫(kù)區(qū)生態(tài)環(huán)境教育部工程研究中心，湖北宜昌 443002；3.水資源安全保障湖北省協(xié)同創(chuàng)新中心，武漢 430072）

0 引 言

從遙感影像中準(zhǔn)確提取水體信息對(duì)于水資源的調(diào)查、規(guī)劃和保護(hù)具有重要意義［1］。黃柏河?xùn)|支流域位于長(zhǎng)江北岸，是湖北省宜昌市的重要水源地，但該流域同時(shí)擁有全國(guó)儲(chǔ)量第二的磷礦資源，磷礦開(kāi)采導(dǎo)致位于黃柏河?xùn)|支上游的玄廟觀水庫(kù)與天福廟水庫(kù)出現(xiàn)了不同程度的水華，且兩座水庫(kù)都處于中等營(yíng)養(yǎng)水平［2］，為遏制污染物向下游水庫(kù)的進(jìn)一步惡化蔓延，需要對(duì)該流域進(jìn)行嚴(yán)格監(jiān)測(cè)；而天福廟水庫(kù)與下游西北口水庫(kù)之間河道水體地處山區(qū)，地形蜿蜒曲折，河道細(xì)小狹長(zhǎng)，水體實(shí)地采樣難度增大。當(dāng)前，遙感是一項(xiàng)迅速發(fā)展的技術(shù)，可以實(shí)現(xiàn)在人跡罕至地區(qū)對(duì)大范圍的實(shí)時(shí)數(shù)據(jù)進(jìn)行長(zhǎng)時(shí)間跨度連續(xù)收集，為地方、區(qū)域以及全球尺度的水環(huán)境變化提供低成本且可靠的信息［3］。利用遙感影像對(duì)流域水體進(jìn)行精確提取，建立水質(zhì)反演模型，對(duì)監(jiān)測(cè)流域水體富營(yíng)養(yǎng)化情況具有明顯優(yōu)勢(shì)［4］；因此基于遙感影像消除噪聲干擾，準(zhǔn)確提取山區(qū)河道水體信息，是合理、正確地應(yīng)用遙感技術(shù)分析水環(huán)境區(qū)域變化的重要基礎(chǔ)。

目前，利用多光譜遙感影像對(duì)大面積水體進(jìn)行自動(dòng)、快速提取已被許多國(guó)內(nèi)外學(xué)者廣泛研究。例如：Work 和Gilmer 等［5］利用單波段閾值法從遙感影像中提取水體；McFeeters［6］提出歸一化差分水體指數(shù)（英文全稱(chēng)，NDWI）算法，利用綠波段和近紅外波段來(lái)描述開(kāi)闊水體；徐涵秋等［7］采用改進(jìn)的NDWI（英文全稱(chēng)，MNDWI）算法，用短波紅外波段（英文全稱(chēng)，SWIR）代替近紅外波段增強(qiáng)水體反射率特征；Qiao 等［8］提出了一種自適應(yīng)的水體提取方法，該方法將水體特征進(jìn)行分層提取；Feyisa 等［9］提出自動(dòng)水體提取指數(shù)（英文全稱(chēng)，AWEI），旨在通過(guò)波段相加、差分運(yùn)算，提高水體像元與非水體像元之間的可分性；陳文倩等［10］提出陰影水體指數(shù)（英文全稱(chēng)，SWI），使山體陰影與水體之間的區(qū)分度顯著提高。上述水體提取方法被廣泛應(yīng)用于提取遙感影像中如水庫(kù)、湖泊等面積較大水域，然而對(duì)于狹長(zhǎng)且處于復(fù)雜地形山區(qū)環(huán)境下的水體，僅利用光譜特性構(gòu)建的水體提取方法并不能全面地避免環(huán)境噪聲，具有一定的局限性［11］。針對(duì)山區(qū)狹長(zhǎng)水體提取，Goumehei 等［12］利用研究區(qū)DEM 高程數(shù)據(jù)消除高海拔山體陰影和暗面影響，相較于僅使用MNDWI 水體指數(shù)精度提高4.88%；朱長(zhǎng)明等［13］在DEM 生成水系數(shù)據(jù)輔助下，采用Landsat 遙感數(shù)據(jù)，提出一種基于自適應(yīng)閾值水體指數(shù)空間優(yōu)化迭代算法，水體提取正確率達(dá)到90%左右；李艷華等［14］基于國(guó)產(chǎn)GF-1 遙感數(shù)據(jù)應(yīng)用基于規(guī)則的面向?qū)ο蠓椒?，輔以DEM 數(shù)據(jù)使山區(qū)細(xì)小水體提取總體精度達(dá)到93.5%，Kappa系數(shù)達(dá)到0.87；薛源等［15］利用GF-1 遙感影像與DEM 高程數(shù)據(jù)相結(jié)合，以總體判別精度為89.5%對(duì)山區(qū)細(xì)小河流邊界進(jìn)行識(shí)別。通過(guò)DEM 數(shù)據(jù)進(jìn)行水文分析對(duì)山區(qū)細(xì)小水體提取邊界進(jìn)行修正的方法雖以較為成熟，但并未確定山體陰影分布，無(wú)法消除與水體相鄰近的陰影噪聲；同時(shí)，喬丹玉等人研究表明，當(dāng)水體處于不同地物背景條件下，提取背景的差異決定了水體提取方法的效率［16］；地物特征的變化是造成提取背景差異的主要原因，在不同水文期典型月，豐、枯水期地物特征呈明顯季節(jié)性變化，如豐水期典型月多處于夏、秋雨水充沛季節(jié)，植被覆蓋率高，地物多以茂盛植被為主；枯水期典型月則多處于春、冬少雨季節(jié)，森林植被稀疏，地物多由裸露山體組成。

因此，本文首先對(duì)比傳統(tǒng)高效的單波段閾值法、NDWI水體指數(shù)和SWI 陰影水體指數(shù)在豐、枯水期不同地物特征背景下水體提取的應(yīng)用效果，確定豐、枯水期典型月水體提取的最適方法；其次對(duì)豐、枯水期水體提取結(jié)果進(jìn)行綜合改進(jìn)處理，針對(duì)山體陰影、水系周?chē)闵⑺w小斑塊以及水體提取結(jié)果斷線現(xiàn)象等問(wèn)題，本文引入遙感成像時(shí)太陽(yáng)高度和太陽(yáng)輻射角度，對(duì)實(shí)驗(yàn)區(qū)中逐個(gè)像元進(jìn)行亮度識(shí)別，得到實(shí)驗(yàn)區(qū)山體陰影圖層，確定實(shí)驗(yàn)區(qū)山區(qū)陰影準(zhǔn)確位置信息，并使用DEM數(shù)據(jù)計(jì)算流域坡度和坡向生成流域河網(wǎng)水系對(duì)水體邊界進(jìn)行校正，最后利用膨脹濾波和Pavlidis 異步細(xì)化填充水體斷線現(xiàn)象處，獲到豐、枯水期山區(qū)河道水體準(zhǔn)確分布信息，以期為后續(xù)流域水資源監(jiān)控與管理提供可靠支撐。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

黃柏河流域位于宜昌市中部，地處東經(jīng)111°04′～111°30′，北緯30°43′～31°29′之間典型峽谷型河流。黃柏河由東、西兩支組成，其中東支發(fā)祥于宜昌市夷陵區(qū)黑良山山脈，是黃柏河的主流，其長(zhǎng)度為130 km；西支發(fā)源于夷陵區(qū)的五郎寨，長(zhǎng)度為78 km［17］。兩支在夷陵區(qū)兩河口匯合后在葛洲壩壩址上游流入長(zhǎng)江。黃柏河流域總面積1 902 km2，年平均氣溫16.9 ℃，屬于亞熱帶季風(fēng)氣候，降水分布具有明顯的季節(jié)性，暴雨、大暴雨集中在6-8月份出現(xiàn)，其3個(gè)月降水量幾乎占全年降水的一半［18］，根據(jù)黃柏河?xùn)|支流域2016年的月降雨量及月流量情況（圖2），選取7月為豐水期典型月，2月為枯水期典型月。黃柏河?xùn)|支流域內(nèi)水體主要以水庫(kù)與河道水體為主，依次建有玄廟觀水庫(kù)、天福廟水庫(kù)、西北口水庫(kù)、尚家河水庫(kù)四座梯級(jí)水庫(kù)。本文選取天福廟水庫(kù)與西北口水庫(kù)之間山區(qū)河道水體為實(shí)驗(yàn)對(duì)象，研究區(qū)地理位置及實(shí)驗(yàn)區(qū)遙感影像如圖1所示。

圖1 研究區(qū)概況圖Fig.1 Overview of the study area

圖2 黃柏河?xùn)|支流域2016年降雨-流量變化圖Fig.2 Rain-flow variation in the East Branch of Huangbai River basin in 2016

1.2 數(shù)據(jù)來(lái)源

國(guó)產(chǎn)GF-1 衛(wèi)星于2013年發(fā)射，是中國(guó)高分辨率對(duì)地觀測(cè)系統(tǒng)中的第一顆衛(wèi)星；GF-1 號(hào)配有兩臺(tái)2 m 分辨率全色/8 m 分辨率多光譜復(fù)合PMS（Panchromatic and Multispectral）傳感器，四臺(tái)16 m 分辨率多光譜WFV（Wide Field of View）傳感器，重返周期4 d。本研究采用GF-1（WFV）遙感影像為數(shù)據(jù)源，數(shù)據(jù)來(lái)源于中國(guó)資源衛(wèi)星應(yīng)用中心（http：//www.cresda.com/CN/），WFV傳感器光譜波段信息如表1 所示；選擇成像時(shí)間為2016年7月14日與2016年2月18日的兩景遙感影像分別作為黃柏河?xùn)|支流域豐水期和枯水期進(jìn)行水體提取實(shí)驗(yàn)。同時(shí)使用到分辨率為30 m×30 m 的數(shù)字高程數(shù)據(jù)（ASTER GDEM），數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云（http：//www.gscloud.cn/search）。

表1 GF-1衛(wèi)星WFV傳感器光譜波段信息Tab.1 Spectral band information of GF-1 satellite WFV sensor

1.3 數(shù)據(jù)預(yù)處理

GF-1 遙感影像預(yù)處理主要包括正射校正、輻射定標(biāo)、大氣校正和影像裁剪。正射校正選擇與待校正高分影像時(shí)相相近的Landsat8-OLI 全色影像為參照，采用ENVI 軟件中RPC（有理多項(xiàng)系數(shù)）正交校正模塊進(jìn)行；繼而將完成正射校正的遙感影像進(jìn)行輻射定標(biāo)，將像元亮度DN 值（Digital Number）轉(zhuǎn)換為具有實(shí)際物理意義的表觀反射率；最后使用ENVI 中FLAASH（Fast Line-of-Sight Atmospheric Analysis of Spectural Hypercubes）模塊進(jìn)行大氣校正。同時(shí)裁剪出研究區(qū)域相應(yīng)的ASTER GDEM 30m數(shù)據(jù)。

1.4 水體提取方法

水體提取方法通?？梢苑譃閱尾ǘ伍撝捣ê投嗖ǘ嗡惴▋煞N方法。單波段閾值法通常是依據(jù)多光譜圖像選擇出水體的特征波段，然后確定該波段的合適閾值，以區(qū)分其他典型地物；多波段組合算法則是利用不同波段之間的反射差異，通過(guò)不同的波段組合來(lái)增強(qiáng)水體與其他典型地物之間的差異。目前，常用的波段組合水體指數(shù)計(jì)算方法有NDWI、改進(jìn)NDWI（MNDWI）、EWI、SWI 和AWEI 等，由于高分一號(hào)衛(wèi)星數(shù)據(jù)的波段數(shù)量有限，因此本文選擇性能較好的單波段閾值法、NDWI水體指數(shù)和SWI陰影水體指數(shù)。水體提取流程如圖3所示。

圖3 水體提取流程圖Fig.3 Flow chart of water extraction

（1）單波段閾值法。單波段閾值法主要是依據(jù)不同地物反射率差異設(shè)定合適的劃分閾值對(duì)水體進(jìn)行分割。通過(guò)對(duì)實(shí)驗(yàn)區(qū)中水體、山體陰影、裸露山體和植被等典型地物進(jìn)行光譜特征分析（圖4），可以發(fā)現(xiàn)裸露山體、植被和山體陰影在Band1、Band2、Band3、Band4 具有相似的變化趨勢(shì)，同時(shí)水體在Band2、Band3、Band4 與其他典型地物無(wú)交集，并且其表觀反射率值普遍低于其他地物，呈逐漸下降趨勢(shì)，在Band4（NIR）處達(dá)到反射率最低值，這一光譜特點(diǎn)是區(qū)分水體與其他典型地物的關(guān)鍵。因此，本文在兩幅GF-1遙感影像數(shù)據(jù)的近紅外波段（NIR）建立一個(gè)閾值進(jìn)行水體提取。

圖4 水體及典型地物光譜曲線圖Fig.4 Spectral curves of water bodies and typical ground objects

（2）NDWI（Normalized Difference Water Index）水體指數(shù)閾值法。

式中：ρGreen和ρNIR分別為綠波段和近紅外波段的反射率數(shù)值。

（3）SWI（Shadow Water Index）陰影水體指數(shù)閾值法。

式中：ρGreen為藍(lán)波段反射率數(shù)值。

1.5 水體提取結(jié)果綜合改進(jìn)處理

（1）水體邊界修正。通過(guò)上述3 種水體提取方法劃定合適的閾值，確保水體得到完整提取，獲得的最佳水體提取結(jié)果不可避免會(huì)存在少量陰影和零散水體小斑塊等噪聲干擾；因此本文引入遙感成像時(shí)太陽(yáng)輻射角度和太陽(yáng)高度（太陽(yáng)高出地平線的角度或坡度），結(jié)合DEM 高程數(shù)據(jù)（ASTER GDEM 30 m）計(jì)算得到的實(shí)驗(yàn)區(qū)坡度和坡向，對(duì)實(shí)驗(yàn)區(qū)中逐個(gè)像元進(jìn)行亮度識(shí)別，得到包含山坡陰影和斜坡暗影等信息的山體陰影圖層（圖5），確定實(shí)驗(yàn)區(qū)山體陰影準(zhǔn)確分布信息，使山體陰影得到全面可視化，并將其與上述提取的水體相疊加，去除陰影噪聲干擾；并依據(jù)DEM高程數(shù)據(jù)獲得的實(shí)驗(yàn)區(qū)河網(wǎng)水系，確定河道水體整體分布，通過(guò)高程落差限定，去除實(shí)驗(yàn)區(qū)河道干流水體周邊零散水體小斑塊，對(duì)經(jīng)過(guò)最佳閾值劃分的河道干流水體提取邊界進(jìn)行修正。

圖5 實(shí)驗(yàn)區(qū)山體陰影與河網(wǎng)水系疊加效果圖Fig.5 Overlay effect of mountain shadow and river network in the experimental area

（2）水體斷線處填充。由于實(shí)驗(yàn)區(qū)水體形態(tài)分布呈狹小細(xì)長(zhǎng)狀，同時(shí)受豐、枯不同水文期水量影響，在提取過(guò)程中部分水體存在明顯的斷線現(xiàn)象。膨脹濾波算法被廣泛應(yīng)用于填充比自身結(jié)構(gòu)元素小的像元結(jié)構(gòu)中［19］，通過(guò)統(tǒng)計(jì)出斷線處像元數(shù)量級(jí)，設(shè)定相同數(shù)量級(jí)的結(jié)構(gòu)元素，即可通過(guò)公式（3）將斷線處加以膨脹連接。Pavlidis 異步細(xì)化算法可以在保持原影像的像元結(jié)構(gòu)基礎(chǔ)上，識(shí)別并標(biāo)記像元骨架，將影像上多像元寬度的直線、曲線沿中心軸線細(xì)化成適當(dāng)寬度［20］。使經(jīng)過(guò)填充后的膨脹水體細(xì)化為實(shí)際寬度的水體。

式中：f(s，t)為圖像輸入函數(shù)；b(x，y)為結(jié)構(gòu)元素；Df、Db分別是f和b的定義域。

1.6 精度驗(yàn)證

總體精度［21］（Overall accuracy，OA）表示所涉及到的所有像元分類(lèi)的正確性；Kappa系數(shù)［22］可以從整體上對(duì)兩幅圖中像元點(diǎn)的一致性進(jìn)行評(píng)價(jià)，針對(duì)實(shí)驗(yàn)區(qū)水體信息提取精度進(jìn)行全面定量評(píng)估，本文通過(guò)從豐枯水期水體提取結(jié)果影像中隨機(jī)抽取100 個(gè)像元點(diǎn)與參照樣本之間建立誤差矩陣，計(jì)算得到總體精度（OA）和Kappa系數(shù)；參照樣本數(shù)據(jù)來(lái)源于同時(shí)期Google Earth 提供的研究區(qū)高空間分辨率影像。Kappa系數(shù)數(shù)值與質(zhì)量分級(jí)關(guān)系見(jiàn)表2［23］。

表2 Kappa系數(shù)數(shù)值與質(zhì)量分級(jí)關(guān)系表Tab.2 Relationship between Kappa coefficient value and quality grade

式中：Oreal為真實(shí)水體信息提取區(qū)域像元數(shù)量；Ototal為標(biāo)準(zhǔn)參照?qǐng)D像中的真實(shí)水體信息像元數(shù)量；K為Kappa系數(shù)；N為總像元數(shù)量；xi+、x+j分別為第i行和第j列總像元數(shù)量；xij為第i行第j列單元上的像元數(shù)量。

2 結(jié)果與分析

單波段閾值法、NDWI 水體指數(shù)和SWI 陰影水體指數(shù)3 種水體提取方法在實(shí)驗(yàn)區(qū)豐、枯水期的水體提取結(jié)果如圖6所示。對(duì)比3 種水體提取方法，在豐水期典型月，水體提取結(jié)果［圖6（a）、6（b）、6（c）］受山體陰影影響較小，歸因于在豐水期實(shí)驗(yàn)區(qū)處于夏季，水量充沛，地物生長(zhǎng)茂盛，僅有少量植被暗影會(huì)被遙感衛(wèi)星所捕捉，通過(guò)單波段閾值法即可獲得良好的提取效果。在枯水期典型月，NDWI 水體指數(shù)［圖6（e）］和SWI 陰影水體指數(shù)［圖6（f）］相比單波段閾值法能較好地提取河道水體，提取結(jié)果受少量山體陰影和零散水體小斑塊影響，而單波段閾值法［圖6（d）］則會(huì)將大量的山體陰影誤認(rèn)為水體被提取；主要因?yàn)樵诳菟趯?shí)驗(yàn)區(qū)正直冬季，森林植被較稀疏，地形連綿起伏，在遙感影像中從低山河谷地帶到高山區(qū)極易存在大量的山體陰影。

圖6 豐（上）、枯（下）水期3種方法水體提取結(jié)果圖Fig.6 Water extraction results of the three methods in the full（top）and dry（bottom）water periods

通過(guò)計(jì)算OA和Kappa系數(shù)（表3）對(duì)3 種水體提取方法在豐、枯水期提取結(jié)果進(jìn)行精度評(píng)價(jià)，表明在豐水期對(duì)水體提取的最適方法為單波段閾值法，Kappa系數(shù)達(dá)到“極好”水平（Kappa系數(shù)=0.80）；枯水期的最適方法為SWI 陰影水體指數(shù)，Kappa系數(shù)同樣達(dá)到“極好”水平（Kappa系數(shù)=0.86）。

表3 豐、枯水期水體提取及綜合改進(jìn)處理結(jié)果精度評(píng)價(jià)結(jié)果Tab.3 Precision evaluation results of water extraction and comprehensive improved treatment in wet and dry seasons

為進(jìn)一步去除噪聲干擾，提高山區(qū)河道水體提取精度，對(duì)豐、枯水期最佳水體提取結(jié)果進(jìn)行綜合改進(jìn)處理，依據(jù)山體陰影圖層和河網(wǎng)水系對(duì)豐枯水期水體提取邊界做修正處理，結(jié)果分別如圖7（a）、7（f），陰影噪聲明顯改善。然而噪聲去除后的山區(qū)河道水體，受混合像元效應(yīng)影響，在河道狹窄、轉(zhuǎn)彎處部分水體斷線現(xiàn)象較為嚴(yán)重（紅色框線），導(dǎo)致局部水體信息丟失，其中圖7（b）、7（g）分別為豐、枯水期水體斷線處放大效果圖；針對(duì)水體斷線處本文使用3×3 的變換單位矩陣進(jìn)行膨脹濾波處理，經(jīng)過(guò)3 次迭代［圖7（c）、7（h）］，水體斷線處得到連接，在膨脹濾波處理基礎(chǔ)上進(jìn)行Pavlidis 異步細(xì)化處理，得到豐枯水期水體提取最終結(jié)果［圖7（d）、7（i）］，從斷線填充效果局部放大圖7（e）、7（j）中可以看到斷線處連接良好。

圖7 豐（上）、枯（下）水期水體提取綜合改進(jìn)處理效果圖Fig.7 Effect diagram of comprehensive improved treatment of water extraction during abundant（top）and dry（bottom）periods

經(jīng)綜合改進(jìn)處理后，豐水期水體提取結(jié)果總體精度從89.54%提高到99.52%，Kappa系數(shù)從0.80 提高到0.98；枯水期總體精度從91.65%提高到99.27%，Kappa系數(shù)從0.86 提高到0.97；水體提取精度達(dá)到“極好”最優(yōu)精度等級(jí)。為進(jìn)一步驗(yàn)證水體提取效果，將綜合改進(jìn)后的面狀水體與真彩色原始遙感影像進(jìn)行疊加，圖8兩幅影像清楚地顯示了豐、枯水期水體提取結(jié)果與實(shí)際水體邊界基本吻合。通過(guò)統(tǒng)計(jì)豐、枯水期水體提取結(jié)果像元個(gè)數(shù)，得知水面面積呈下降趨勢(shì)，在豐水期水面面積達(dá)到350 km2，當(dāng)枯水期到來(lái)時(shí)水面面積減小到194 km2。

圖8 豐（左）、枯（右）水期水體提取結(jié)果疊加圖Fig.8 Superposition diagram of water extraction results in abundant（left）and dry（right）periods

3 結(jié) 論

本文基于GF-1遙感數(shù)據(jù)提取了2016年豐枯水期典型月黃柏河?xùn)|支流域上游天福廟水庫(kù)與下游西北口水庫(kù)之間山區(qū)河道水體，針對(duì)不同水文期確定了最佳水體提取方法，并對(duì)豐、枯水期最佳提取結(jié)果進(jìn)行綜合改進(jìn)處理，得到以下結(jié)論。

（1）通過(guò)對(duì)比單波段閾值法、NDWI 水體指數(shù)、SWI 陰影水體指數(shù)三種水體提取方法在不同水文期的提取效果，發(fā)現(xiàn)在豐水期，水體提取干擾噪聲主要以植被暗影為主，單波段閾值法為最佳提取方法，水體提取結(jié)果總體精度達(dá)到89.54%，Kappa系數(shù)為0.80；在枯水期，實(shí)驗(yàn)區(qū)遙感影像中存在大量山體陰影，SWI 陰影水體指數(shù)為最佳提取方法，水體提取結(jié)果總體精度達(dá)到91.65%，Kappa系數(shù)為0.86。

（2）基于山體陰影圖層、河網(wǎng)水系、膨脹濾波和Pavlidis 異步細(xì)化算法對(duì)豐、枯水期水體最佳提取結(jié)果做綜合改進(jìn)處理后，豐水期水體提取總體精度提高到99.52%，Kappa系數(shù)提高到0.98；枯水期水體提取總體精度提高到99.27%，Kappa系數(shù)提高到0.97；水體提取精度均達(dá)到“極好”最優(yōu)精度等級(jí)。