基于像元二分模型改進(jìn)的水體指數(shù)提取方法研究

丁鵬飛， 劉漢湖

（成都理工大學(xué)國(guó)土資源部地學(xué)空間信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，成都 610059）

青藏高原地區(qū)分布著世界上海拔最高、數(shù)量最多的內(nèi)陸湖泊，其中大于50 km2的湖泊有138個(gè)，在全球變暖背景下，氣溫升高造成的冰川融化和降水增多都影響湖泊面積變化，因此準(zhǔn)確監(jiān)測(cè)這些湖泊水體面積變化對(duì)研究青藏高原地區(qū)氣候和生態(tài)具有非常重要的意義［1-2］. 近年來遙感技術(shù)的進(jìn)步使得監(jiān)測(cè)高海拔地區(qū)湖泊面積變化的研究逐漸增多，歸一化水體指數(shù)模型（Normalized Difference Water Index，NDWI）利用水體的反射波譜曲線和其他地物的不同提取水體得到廣泛運(yùn)用. 徐涵秋［3］針對(duì)NDWI容易將建筑物提取為水體的困難提出改進(jìn)歸一化水體指數(shù)（The Modified Normalized Difference Water Index，MNDWI）. 德吉央宗等［4］利用TM 影像研究青藏高原地區(qū)湖泊變化，探究氣候變化與湖泊變化存在滯后的相關(guān)性. 張麗華等［5］利用MODIS數(shù)據(jù)對(duì)烏梁素海進(jìn)行水體遙感監(jiān)測(cè)，取得較好效果. Xu等［6］利用MODIS數(shù)據(jù)對(duì)長(zhǎng)江中下游地區(qū)四大湖泊鄱陽(yáng)湖、洞庭湖、太湖和巢湖進(jìn)行面積監(jiān)測(cè)，分析了湖泊面積的驅(qū)動(dòng)力，對(duì)該地區(qū)湖泊保護(hù)具有重要意義. 劉家福等［7］利用TM影像分析2000年至2015年中亞地區(qū)主要湖泊和河渠的水體動(dòng)態(tài)變化. 曹萌萌等［8］利用MODIS 數(shù)據(jù)對(duì)洞庭湖面積萎縮速率及水華分布范圍進(jìn)行分析，發(fā)現(xiàn)水華主要集中發(fā)生在東洞庭湖的西部湖灣區(qū). 但是這些研究都是基于像元的定性研究，由于空間分辨率的限制，混合像元總是存在的，對(duì)于地物的識(shí)別和提取的影響很大，混合像元分解模型的提出有效地解決了這類問題［9-10］. 線性分解模型主要原理是假設(shè)每個(gè)像元的總反射率由各地物反射率經(jīng)過線性組合而成，其主要步驟為端元的選擇，但對(duì)于MODIS 數(shù)據(jù)而言，每個(gè)像元代表了250 m×250 m 的實(shí)際范圍，實(shí)際地物的復(fù)雜分布給端元的選擇和豐度提取帶來很大干擾［11-12］.

像元二分模型也是混合像元分解的一種方法，被廣大學(xué)者運(yùn)用于估算某地區(qū)植被覆蓋度的研究中. 張成才等［13］利用像元二分模型對(duì)伏牛山19年的NDVI數(shù)據(jù)進(jìn)行反演，同時(shí)考慮氣溫和降水對(duì)植被覆蓋度的關(guān)系，發(fā)現(xiàn)氣溫是影響植被覆蓋度變化的主要因子. 王金強(qiáng)等［14］基于像元二分模型發(fā)現(xiàn)新疆石河子地區(qū)20年來植被覆蓋度呈現(xiàn)先增加后減少的特征，分析其主要原因是地下水的過度開采和城市化等人為因素. 除了植被，像元二分模型也被運(yùn)用于提取其他地物類型，章瓊［15］等將像元二分模型運(yùn)用到土壤重金屬反演方面，發(fā)現(xiàn)經(jīng)過像元二分模型處理的土壤反射率與重金屬Cr兩者敏感性顯著相關(guān)，取得較好效果. 丹卓宇等［16］利用改進(jìn)的像元二分模型得到落葉松人工林郁閉度分布圖. 趙秀霞等［17］利用像元二分模型對(duì)艾比湖濕地內(nèi)土地覆蓋類型進(jìn)行估算，并研究了植被覆蓋對(duì)氣候的響應(yīng)程度. 像元二分模型運(yùn)用在水體提取方面還較少，本文采用Landsat8 OLI數(shù)據(jù)，分別利用線性分解模型和改進(jìn)歸一化水體指數(shù)（MNDWI）的基礎(chǔ)上的像元二分模型計(jì)算研究區(qū)水體面積，以高分辨率影像為參考，對(duì)兩種分解模型提取的面積差異及區(qū)域適用性進(jìn)行評(píng)價(jià).

1 研究區(qū)與數(shù)據(jù)處理

1.1 研究區(qū)概況

研究區(qū)全部位于青藏高原地區(qū)，位于30°～32°N，84°～91°E之間. 考慮湖泊面積大小以及形狀等要素，本文選擇了研究區(qū)內(nèi)納木錯(cuò)、色林錯(cuò)、當(dāng)惹雍錯(cuò)和扎日南木錯(cuò)等四個(gè)不同類型的湖泊，研究區(qū)位置如圖1所示.①納木錯(cuò)位于西藏自治區(qū)中部，是西藏第二大湖泊，海拔4718 m，是世界上海拔最高的大型湖泊，形狀近似長(zhǎng)方形，其湖水來源主要是流域內(nèi)冰川的融水以及降水形成的徑流，據(jù)統(tǒng)計(jì)納木錯(cuò)每年的出、入水量基本平衡［18-19］. ②位于那曲地區(qū)的色林錯(cuò)是西藏最大湖泊，也是中國(guó)第二大咸水湖，湖面海拔約4530 m，流域內(nèi)許多河、湖串通形成內(nèi)湖泊群，形狀很不規(guī)則，河水主要靠冰雪融水補(bǔ)給. 在1976—2009年時(shí)間段內(nèi)湖泊面積不斷擴(kuò)大，增長(zhǎng)幅度高達(dá)39.39%，色林錯(cuò)面積增長(zhǎng)而引發(fā)的一系列不良后果包括下游牧場(chǎng)草地不斷被淹，以及影響大北線的正常通行［20］. ③當(dāng)惹雍錯(cuò)是中國(guó)第二深的湖泊，該湖南北走向，三面環(huán)山，形狀如同鞋底，海拔約4600 m，最大深度230 m，相較于全盛時(shí)期，湖泊退化嚴(yán)重. ④位于阿里地區(qū)措勤縣境內(nèi)的扎日南木錯(cuò)，屬于東西向構(gòu)造斷陷湖，湖泊形態(tài)不規(guī)則，南北兩岸較窄，東西兩岸地勢(shì)開闊，湖水主要靠冰雪融水補(bǔ)給［21］.

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

本文使用的數(shù)據(jù)都是來自于地理空間數(shù)據(jù)云網(wǎng)站的Landsat8 OLI 的L1T 級(jí)別數(shù)據(jù)，都經(jīng)過了系統(tǒng)輻射校正和幾何校正. Landsat8衛(wèi)星是美國(guó)NASA在2013年2月11日發(fā)射的第八個(gè)Landsat系列衛(wèi)星，其上搭載了兩個(gè)傳感器，分別為陸地成像儀（OLI）和熱紅外傳感器（TIRS），OLI陸地成像儀包括9個(gè)波段，空間分辨率為30 m，其中全色波段分辨率為15 m. 影像獲取時(shí)間及條代號(hào)分別為：①納木錯(cuò). 2015年9月2日（行列號(hào)138/39）；②色林錯(cuò). 2014 年10 月17 日（行列號(hào)139/38）和2015 年11 月19 日（行列號(hào)140/38）；③當(dāng)惹雍錯(cuò).2016年5月8日（行列號(hào)141/38）和2016年5月8日（行列號(hào)141/39）；④扎日南木錯(cuò). 2016年5月18日（行列號(hào)141/39）. 為了增強(qiáng)水陸差異，更好地提取水體，選用波段543生成標(biāo)準(zhǔn)假彩色圖像.

圖1 研究區(qū)位置圖Fig.1 Location map in the study area

由于受衛(wèi)星本身飛行姿態(tài)的改變以及大氣散射等影響會(huì)使遙感影像產(chǎn)生幾何和輻射畸變，在進(jìn)行分類提取研究之前要對(duì)影像進(jìn)行預(yù)處理，主要包括輻射定標(biāo)、大氣校正和區(qū)域裁剪及融合等. 對(duì)輻射定標(biāo)后的數(shù)據(jù)利用FLAASH 大氣校正模塊進(jìn)行校正，消除大氣影響. 對(duì)于在兩景影像上的研究區(qū)，首先利用Seamless Mosaic 對(duì)影像進(jìn)行鑲嵌，再根據(jù)選擇的研究區(qū)進(jìn)行裁剪. 精度驗(yàn)證數(shù)據(jù)選用Google Earth 數(shù)據(jù)，本文采用91 衛(wèi)圖下載的免費(fèi)無水印的與Landsat8 OLI數(shù)據(jù)同時(shí)相的第15級(jí)別的納木錯(cuò)（4.11 m）、色林錯(cuò)（4.06 m）、當(dāng)惹雍錯(cuò)（4.09 m）和扎日南木錯(cuò)（4.1 m）數(shù)據(jù).

2 水體提取方法

2.1 像元二分模型改進(jìn)的MNDWI

水體的反射波譜主要反映在藍(lán)綠波段，在其他波段主要呈吸收，特別是近紅外和中紅外波段為0，這是水體與其他地物在波譜上的主要差異. 裸地和建筑物等地物在近紅外和中紅外波段反射率很高，而植被在近紅外波段反射率約在40%～50%，在藍(lán)波段和紅波段吸收，綠波段反射，因此歸一化水體指數(shù)模型（NDWI）利用綠波段和近紅外波段進(jìn)行計(jì)算，抑制了植被和裸地背景對(duì)水體提取的影響. 而徐涵秋提出的改進(jìn)歸一化水體指數(shù)（MNDWI）用中紅外波段代替原來的近紅外波段，解決了NDWI 容易將建筑物提取為水體的問題，其原理如下：

其中，Green和SWIR分別代表水體在綠波段和近紅外波段的反射率.

像元二分模型是估算植被覆蓋度的模型，利用綠色植物強(qiáng)吸收可見光波段和高反射近紅外波段特點(diǎn)，增強(qiáng)植被信號(hào)，削弱噪音部分，其原理是假設(shè)每個(gè)像元的反射率由純植被部分反射率和非植被部分反射率組成，那么任一像元的反射率都可以表示為植被部分和非植被部分的線性加權(quán)和. 將其運(yùn)用到提取水體方面，則可以將任意像元的反射率表示為水體部分反射率( RW)和非水體部分反射率( RS)的線性加權(quán)之和：

假設(shè)影像中每個(gè)像元的有水體部分的面積比例為fw，即該像元的水體覆蓋指數(shù)，那么非水體部分為1-fw. 如果該像元全部由水體組成，則該像元反射率可以表示為Rwater；如果該像元無水體組成，則反射率為Rsoil. 因此混合像元的水體部分貢獻(xiàn)的信息( RW)可以表示為純水體反射率( Rwater)與像元中水體面積比例( fw)的乘積，而非水體部分貢獻(xiàn)的信息( RS)可以表示為非水體反射率( Rsoil)與1-fw的乘積.

通過解算公式（3）～（5），可以得到計(jì)算水體覆蓋指數(shù)的公式，如下：

根據(jù)像元二分模型原理我們可以將一個(gè)像元的MNDWI值表示為有水體部分地表和沒有水體部分地表組成的形式，因此計(jì)算水體覆蓋指數(shù)的公式可以表示為：

其中：MNDWIwater表示由全部水體組成的MNDWI 值；MNDWIsoil表示由全部非水體組成的MNDWI 值，利用這兩個(gè)具有實(shí)際意義的值可以削弱大氣和土壤背景的影響. 與其他地物不同，MNDWI直方圖呈現(xiàn)出雙峰的特征，可以幫助我們確定MNDWIwater和MNDWIsoil的值，只要確定這兩個(gè)參數(shù)就可以估算水體覆蓋指數(shù)，從而估算研究區(qū)內(nèi)水體面積. 在這里我們統(tǒng)計(jì)MNDWI的最值、均值和標(biāo)準(zhǔn)方差以及像元個(gè)數(shù)和百分比，分別利用頻率分布直方圖和置信區(qū)間找到水體和非水體兩類地物的拐點(diǎn)，作為MNDWIwater和MNDWIsoil的值.

在分類結(jié)果中fw=0 即為黑色區(qū)域代表陸地面積，fw=1即為白色區(qū)域代表水體面積，對(duì)于介于0和1之間的fw則屬于水體與陸地的混合區(qū)域，在這里我們需要設(shè)定一個(gè)閾值來區(qū)分統(tǒng)計(jì)水體和陸地，根據(jù)研究表明當(dāng)fw=標(biāo)準(zhǔn)差（Stedv）時(shí)，效果較好，根據(jù)每個(gè)研究區(qū)的標(biāo)準(zhǔn)差分為水體和非水體，再統(tǒng)計(jì)水體與非水體的面積（如圖2，表1所示）.

圖2 基于MNDWI頻率直方圖確定MNDWIwater和MNDWIsoil的值Fig.2 Determining the values of MNDWIwater and MNDWIsoil based on the MNDWI frequency histograms

表1 基于像元二分模型水體提取結(jié)果標(biāo)準(zhǔn)差Tab.1 Standard deviations of the water extraction results based on the dimidiate pixel model

2.2 線性分解模型

線性光譜解混是處理高光譜影像分類中混合像元問題的一種常用方法，此方法包括端元提取和像元分解兩個(gè)步驟. 端元提取，目的是提取出“純”地物的光譜，目前提取端元的方法主要有兩大類，一是在室內(nèi)通過光譜儀進(jìn)行光譜測(cè)量，第二個(gè)就是通過分類或者主成分分析的方法由遙感圖像得到，其主要方法有PPI算法、N-FINDER方法以及SMACC算法等；其次是混合像元分解，即利用端元的線性組合來表示混合像元，并計(jì)算各端元豐度. 線性混合模型，假定每個(gè)像元的反射率等于各個(gè)端元反射率的線性組合，模型表達(dá)式見式（7）.

其中：i和j分別表示光譜波段和地物類；R( i )是每個(gè)像元在光譜波段i的觀測(cè)反射率；Re( i,j )為第j個(gè)端元在光譜波段i的反射率；ε為誤差.

為了減少工作數(shù)據(jù)量，提高數(shù)據(jù)處理效率對(duì)預(yù)處理后的影像，要對(duì)影像進(jìn)行最小噪聲分離（MNF）變換，利用純凈像元指數(shù)算法和n維可視化工具提取端元波譜，根據(jù)影像識(shí)別每條波譜曲線所代表地類類型，選擇端元波譜分解影像即可得到水體與其他端元等非水體的面積.

3 結(jié)果與分析

3.1 水體提取結(jié)果

圖3 不同模型提取水體結(jié)果分類圖Fig.3 Classification diagrams of water extraction results from different models

利用線性混合像元分解模型得到水體豐度圖（見圖3，從左到右依次為納木錯(cuò)、色林錯(cuò)、當(dāng)惹雍錯(cuò)和扎日南木錯(cuò)的經(jīng)過不同模型提取的水體面積），根據(jù)水體豐度圖設(shè)定豐度大于0.5的像元為水體端元，得到水體與非水體的分類圖（圖3（a））；將改進(jìn)歸一化水體指數(shù)模型閾值設(shè)置為0提取水體結(jié)果分類圖（圖3（b））；基于改進(jìn)歸一化水體指數(shù)的像元二分模型提取到的水體分類圖（圖3（c）），通過統(tǒng)計(jì)可以得到納木錯(cuò)研究區(qū)內(nèi)共有水體像元2 266 022個(gè)，水體面積為2 039.42 km2，占總面積的54.09%；色林錯(cuò)研究區(qū)共有水體像元2 580 678個(gè)，水體面積為2 322.61 km2，占總面積的39.66%；當(dāng)惹雍錯(cuò)研究區(qū)共有水體像元962 222個(gè)，水體面積為866.67 km2，占總面積的28.6%；扎日南木錯(cuò)研究區(qū)共有水體像元1 159 356 個(gè)，水體面積為1 043.42 km2，占總面積的45.89%（如表2所示）.

表2 像元二分模型提取水體結(jié)果Tab.2 Water extracted results by the dimidiate pixel model

3.2 不同模型提取水體比較

分別將利用線性分解模型，MNDWI模型和像元二分模型提取到的水體面積與高分辨率影像提取的面積進(jìn)行對(duì)比（表3）. 分析表3數(shù)據(jù)可以發(fā)現(xiàn)，把兩類方法提取到的水體面積和較高分辨率影像統(tǒng)計(jì)面積相對(duì)比，基于像元二分模型的MNDWI提取的水體面積信息比線性分解更接近. 在納木錯(cuò)地區(qū)利用線性分解提取到的面積為2 022.05 km2，像元二分模型提取到的面積為2 039.42 km2，高分辨率影像統(tǒng)計(jì)面積為2 037.65 km2，像元二分模型提取的結(jié)果更接近實(shí)際面積. 在色林錯(cuò)地區(qū)利用線性了解提取水體面積為2 309.43 km2，像元二分模型提取的水體面積為2 322.61 km2，高分辨率影像統(tǒng)計(jì)水體面積為2 390.56 km2，說明像元二分模型優(yōu)于線性分解，但是由于色林錯(cuò)地區(qū)很多河流湖泊串通合并造成水陸邊界不明顯，混合像元比例增多，再加上分辨率不同等誤差，使得水體提取面積與高分辨率影像有較大差別. 在當(dāng)惹雍錯(cuò)地區(qū)利用線性分解模型提取水體面積為862.54 km2，像元二分模型提取到的面積為866.67 km2，高分辨率影像統(tǒng)計(jì)面積為868.54 km2，說明在該地區(qū)兩種模型提取水體效果差距不大，且都接近于真實(shí)面積. 在扎日南木錯(cuò)地區(qū)利用線性分解模型提取水體面積為1 038.23 km2，像元二分模型提取到的面積為1 043.42 km2，高分辨率影像統(tǒng)計(jì)面積為1 043.42 km2. 總體而言，改進(jìn)歸一化水體指數(shù)（MNDWI）在考慮混合像元問題上，提取水體效果較差；在MNDWI基礎(chǔ)上建立像元二分模型提取水體面積比線性模型提取水體效果較好.

表3 不同方法提取水體面積比較Tab.3 Comparison of the areas of water bodies extracted by different methods 單位：km2

以高分辨影像得到的分類結(jié)果為參照，選用混淆矩陣對(duì)線性分解模型和像元二分模型提取到4個(gè)研究區(qū)水體進(jìn)行精度分析（表4和表5）.

通過比較表4和表5，可以發(fā)現(xiàn)基于MNDWI的像元二分模型提取的納木錯(cuò)、色林錯(cuò)、當(dāng)惹雍錯(cuò)和扎日南木錯(cuò)的Kappa 系數(shù)分別為0.942 6、0.863 2、0.924 5 和0.963 1，相比較于線性分解模型提取的Kappa 系數(shù)0.930 9、0.854 1、0.908 0和0.943 2，分類精度更高. 比較各研究區(qū)Kappa系數(shù)會(huì)發(fā)現(xiàn)像元二分模型和線性分解模型在色林錯(cuò)地區(qū)的分類精度都不是很高，可能是因?yàn)樯皱e(cuò)地區(qū)許多河、湖串通，組成了一個(gè)內(nèi)陸湖群，造成支流邊界模糊，分類達(dá)不到理想效果.

表4 線性分解模型分類結(jié)果精度評(píng)價(jià)Tab.4 Accuracy evaluations of classification results with linear decomposition model

表5 像元二分模型分類結(jié)果精度評(píng)價(jià)Tab.5 Accuracy evaluations of classification results with the dimidiate pixel model

4 討論

1）將像元二分模型從植被方面運(yùn)用到水體領(lǐng)域，與改進(jìn)歸一化水體指數(shù)相結(jié)合提取水體信息，在研究區(qū)內(nèi)取得較好效果，但還需要大量研究來證明方法的可行性.

2）線性分解雖然是現(xiàn)在常用的像元分解方法，受復(fù)雜地物和主觀因素的影響，精確選擇端元和提取端元豐度仍有難度. 像元二分模型基礎(chǔ)上的MNDWI 模型利用頻率直方圖和置信度區(qū)間來準(zhǔn)確確定參數(shù)，減少水體提取過程中的誤差.

3）本文利用像元二分模型改進(jìn)的水體指數(shù)模型提取形狀規(guī)則且面積較大的湖泊精度較高，對(duì)區(qū)域湖泊面積調(diào)查和氣候變化研究提供一定依據(jù)，但對(duì)于細(xì)小面積湖泊的適用性有待進(jìn)一步研究.