0  引言

水庫(kù)和湖泊是地表水體的主要存在形式，對(duì)河流下游地區(qū)調(diào)蓄洪水、維持水沙平衡具有重要作用。湖泊及水庫(kù)蓄水量的研究一直是防洪調(diào)蓄研究的重點(diǎn)，近年來(lái)遙感技術(shù)和地理信息系統(tǒng)技術(shù)越來(lái)越多地被國(guó)內(nèi)外的研究人員應(yīng)用于湖泊、水庫(kù)監(jiān)測(cè)研究中。

隨著遙感和航空數(shù)據(jù)的種類不斷豐富，這些數(shù)據(jù)在水體蓄水量監(jiān)測(cè)中得到了廣泛的應(yīng)用。水體面積、水下及周邊地形是計(jì)算水庫(kù)庫(kù)容的重要因素。不少學(xué)者提出了利用遙感圖像的歸一化指數(shù)和譜間關(guān)系來(lái)提取水體面積。MCF S K等^［1］利用綠波段(Green)和近紅外波段(NIR)構(gòu)建歸一化差異水體指數(shù)NDWI提取水體，消除陰影對(duì)水體的提取影響；為改善NDWI對(duì)建筑用地及居民地的誤分現(xiàn)象，徐涵秋等^［2］提出利用Landsat TM數(shù)據(jù)中的第5波段(MIR)替代NIR，構(gòu)建了改進(jìn)歸一化差異指數(shù)MNDWI，效果好于NDWI指數(shù)法；趙紫薇等^［3］結(jié)合NDWI與MNDWI提取水體信息的特點(diǎn)提出了一種新的水體指數(shù)NMWI；楊存建等^［4］對(duì)水體及背景物在影像上各個(gè)波段的發(fā)射率光譜特征曲線進(jìn)行綜合分析發(fā)現(xiàn)水體的TM數(shù)據(jù)存在(TM2+TM3＞TM4+TM5)的關(guān)系；汪金花等^［5］在上述關(guān)系中加入(TM4/TM2＜0.88)使模型能夠較好地區(qū)分水體和居民地；劉桂林等^［6］將反映地物的濕度信息KT3與TM波段結(jié)合建立譜間關(guān)系提取水體，精度高于指數(shù)法。

目前，結(jié)合遙感數(shù)據(jù)監(jiān)測(cè)水體庫(kù)容主要有以下兩種思路。一是用于庫(kù)容曲線的校核^［7-11］；二是用于三維分析，直接測(cè)算水體庫(kù)容^［12-14］。本文借鑒前述成果，利用指數(shù)法從Landsat 8 OLI影像中提取水體，結(jié)合資源三號(hào)衛(wèi)星獲取的DEM數(shù)據(jù)實(shí)現(xiàn)對(duì)遼寧省撫順市大伙房水庫(kù)動(dòng)態(tài)蓄水量的反演，并與水文資料進(jìn)行對(duì)比。

1  反演原理和方法

1.1  水體提取

格式歸一化差異水體指數(shù)(NDWI)、修正的歸一化差異水體指數(shù)(MNDWI)是常用的水體提取方法，在某些研究區(qū)歸一化差異植被指數(shù)(NDVI)也被用來(lái)提取水體信息。

歸一化差異水體指數(shù)NDWI是利用植被、裸地與水體在可見(jiàn)光波段與近紅外波段的反射差異突出水體信息，計(jì)算公式如式(1)所示：

式中，Green，NIR分別對(duì)應(yīng)Landsat 8 OLI影像的Band3、Band5。

修正的歸一化差異水體指數(shù)MNDWI是為了進(jìn)一步對(duì)水體和建筑物進(jìn)行區(qū)分而提出的，計(jì)算公式如式(2)所示：

式中，MIR為中紅外波段，對(duì)應(yīng)Landsat 8 OLI影像的Band6。

在水體提取中，若采用同一閾值，會(huì)造成較大錯(cuò)分、漏分現(xiàn)象。例如NDWI和MNDWI指數(shù)法，一般采用0作為閾值，但在實(shí)際影像提取中，須采用較大或較小的閾值才能獲得理想的效果。研究中采用MNDWI指數(shù)和紅外波段進(jìn)行水體提取，利用OTSU法確定閾值^［15］。圖1顯示的是本研究提取的水體區(qū)域。

1.2 立體像對(duì)提取

立體像對(duì)提取DEM是基于雙像立體測(cè)圖原理。雙像立體測(cè)圖，是利用一個(gè)立體像片對(duì)，在恢復(fù)它們的內(nèi)、外方位元素后，重建與地面相似的幾何模型，并對(duì)該模型進(jìn)行量測(cè)的一種攝影測(cè)量方法。如圖2所示，像片p1，p2和攝站點(diǎn)的位置S1，S2是真實(shí)的攝影過(guò)程，根據(jù)攝影過(guò)程的可逆性，恢復(fù)它們之間的相互位置關(guān)系，找到同名點(diǎn)，即可構(gòu)成與地面完全相同的幾何模型；在實(shí)際測(cè)量中，將攝站S2移動(dòng)到S2′的位置后的模型與真實(shí)的地面模型相似，模型的比例尺為 1:m=b:B，實(shí)現(xiàn)攝影過(guò)程的反轉(zhuǎn)^［16］。本研究采用ENVI軟件對(duì)資源三號(hào)立體衛(wèi)星數(shù)據(jù)進(jìn)行DEM提取。DEM Extraction是ENVI的DEM自動(dòng)提取擴(kuò)展模塊，它能夠簡(jiǎn)單、快速地從掃描、數(shù)字航空影像，或者沿軌道方向、垂直軌道方向的推掃式衛(wèi)星傳感器創(chuàng)建DEM。

DEM Extraction模塊的提取流程如圖3所示。

利用該模塊進(jìn)行DEM提取時(shí)，需提供研究區(qū)的最大、最小高程進(jìn)行精度控制。本研究利用ASTGTM2數(shù)據(jù)進(jìn)行測(cè)算。圖4為DEM的效果對(duì)比圖。

2   水庫(kù)庫(kù)容反演

2.1   數(shù)據(jù)選擇與處理

2.1.1   Landsat 8數(shù)據(jù)

Landsat 8衛(wèi)星于 2013 年 2月 11日發(fā)射成功，2013年5月30日開(kāi)始免費(fèi)提供數(shù)據(jù)。該衛(wèi)星共攜帶有OLI和TIRS兩個(gè)載荷，其中OLI陸地成像儀設(shè)有9個(gè)波段(Band1~Band9)，TIRS多光譜波段圖像空間分辨率為30 m，全色波段圖像的空間分辨率為15 m，幅寬為185 km，數(shù)據(jù)更新周期為16天。相比于ETM+數(shù)據(jù)，OLI傳感器將Band5的波段范圍調(diào)整為0.845~0.885 μm，排除了0.825 μm處的水汽吸收特征；全色波段Band8范圍變窄，更有利于區(qū)分出植被信息；此外，新增了用于海岸觀測(cè)的藍(lán)波段Band1和用于云檢測(cè)的短波紅外波段Band9。Landsat 8影像詳細(xì)的技術(shù)參數(shù)如表1所示。

本研究數(shù)據(jù)來(lái)源于中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)(http://www.gscloud.cn)，共下載了27幅影像，條帶號(hào)119/31，所下載影像的成像時(shí)間涉及從2013年至2016年，研究區(qū)內(nèi)的云層覆蓋率為5%以下，影像質(zhì)量較好。數(shù)據(jù)列表如表2所示。

2.1.2   資源三號(hào)數(shù)據(jù)

資源三號(hào)衛(wèi)星發(fā)射于2012年 1月9日，是我國(guó)當(dāng)時(shí)第一顆高分辨率立體測(cè)圖衛(wèi)星。該衛(wèi)星搭載了四臺(tái)光學(xué)相機(jī)，包括2臺(tái)分辨率為3.5 m×3.7 m的前、后視TDI CCD相機(jī)，1臺(tái)分辨率為2.1 m的正視全色相機(jī)和1臺(tái)分辨率5.8 m多光譜相機(jī)。前、后視相機(jī)與正視相機(jī)夾角為±22°，對(duì)應(yīng)的基高比為0.85~0.95，幅寬為52 km，滿足1∶50 000比例尺立體測(cè)圖要求和1∶25 000比例尺地形圖更新需求。衛(wèi)星過(guò)境一次即可獲得三視立體影像和多光譜影像，影像可以組成同軌立體，重訪周期為5天。

本研究所采用數(shù)據(jù)購(gòu)買于遙感集市(www.rscloudmart.com)。雖然資源三號(hào)衛(wèi)星數(shù)據(jù)更新周期為5天，但因水位消漲、植被的影響，庫(kù)區(qū)地貌因時(shí)而異，圖像鑲嵌亦會(huì)造成不必要的誤差，因而選擇影像時(shí)以一幅影像能完全覆蓋研究區(qū)為原則。同時(shí)考慮云霧天氣對(duì)圖像質(zhì)量的影響，這就造成可選數(shù)據(jù)較少，從可選數(shù)據(jù)中進(jìn)一步選取水位最低所對(duì)應(yīng)時(shí)相的影像。共篩選出2014年6月24日和2015年6月28日兩天的影像，其中，6月24日的影像為前視、正視與后視，6月28日為前視與后視。具體信息見(jiàn)表3。

2.2   水體面積和完整DEM提取

經(jīng)過(guò)ENVI和Arcgis處理后得到的水體面積如表4所示。

其中日期為2016/04/22、2015/03/19、2014/03/16、2013/04/14的四幅影像為冰面，日期2013/05/16的影像水面為薄冰。因水體狀態(tài)的改變，與液態(tài)水水力條件不同，導(dǎo)致其水體淹沒(méi)范圍大于相同水位時(shí)的液態(tài)水淹沒(méi)范圍。因而，這5幅影像與其他時(shí)相的影像不具有可比性，將不予采用。

對(duì)篩選后的水體面積和相應(yīng)水位數(shù)據(jù)進(jìn)行曲線擬合，結(jié)果如圖5所示。

利用資源三號(hào)提取的DEM水面值有異常，研究采用日期為2014/06/04的影像，其水體面要素文件的淹沒(méi)范圍大于DEM異常值區(qū)域范圍。利用該影像的水面插值結(jié)果對(duì)數(shù)據(jù)異常區(qū)域進(jìn)行替換，即可得到相對(duì)最低水位時(shí)研究區(qū)完整的DEM。處理結(jié)果如圖6所示。

2.3  反演結(jié)果驗(yàn)證

在ArcGIS中將水面高程面要素文件轉(zhuǎn)換為柵格格式。利用對(duì)應(yīng)時(shí)相的面要素文件對(duì)處理后的DEM裁剪，得到水體淹沒(méi)部分的DEM。此時(shí)，水面高程文件和DEM像元是一一對(duì)應(yīng)的，整個(gè)水域由一個(gè)個(gè)緊密排列的四棱柱構(gòu)成。通過(guò)計(jì)算每個(gè)水柱的體積再進(jìn)行累計(jì)，即可得到相對(duì)最低水位以上的水體體積。在ENVI中利用BAND MATH工具將兩個(gè)柵格文件數(shù)據(jù)相減，計(jì)算整幅影像的像元值之和即為蓄水量變化值。

利用從大伙房水庫(kù)管理局收集的水位-庫(kù)容值，對(duì)其按三次關(guān)系進(jìn)行擬合，可以建立水位-庫(kù)容曲線，見(jiàn)圖7。

本研究以水位-庫(kù)容曲線為依據(jù)，計(jì)算各個(gè)時(shí)相庫(kù)容相比于2014/06/04成像時(shí)的庫(kù)容增量ΔV1。將利用庫(kù)容曲線計(jì)算的庫(kù)容變化值和利用本研究方法計(jì)算的蓄水量變化值進(jìn)行比較，如表5所示。通過(guò)計(jì)算二者差值發(fā)現(xiàn)，利用該方法計(jì)算的蓄水量變化值ΔV2與利用庫(kù)容曲線計(jì)算的庫(kù)容變化量ΔV1差值占ΔV1的10%~20%之間。

經(jīng)分析，造成誤差主要有兩方面原因:一是通過(guò)庫(kù)容曲線計(jì)算的水體庫(kù)容假設(shè)水面為水平面，而實(shí)際水面并非水平面；二是本研究方法受遙感影像精度、DEM和多源遙感配準(zhǔn)精度、插值方法等因素的影響，造成與實(shí)際蓄水量變化的誤差。

3 結(jié)論

本研究使用Landsat 8數(shù)據(jù)，采用MNDWI指數(shù)進(jìn)行水體提取，便于利用OTSU法確定閾值。當(dāng)水面結(jié)冰或存在薄冰時(shí)，水體面積及邊界線與其他時(shí)相不具有可比性，故不予討論。根據(jù)提取的水體，可以估算大伙房水庫(kù)的面積；根據(jù)資源三號(hào)衛(wèi)星立體像對(duì)提取的DEM，因其受天氣影響較大，大面積水域中DEM數(shù)值異常，通過(guò)對(duì)其進(jìn)行圖像配準(zhǔn)和異常值處理，結(jié)合水面插值結(jié)果得到完整的DEM用于蓄水量監(jiān)測(cè)。

利用大伙房水庫(kù)管理局的水位和庫(kù)容值建立庫(kù)容曲線，將日期為2014/06/04的一幅影像的水位作為相對(duì)最低水位，得到不同時(shí)相水庫(kù)的庫(kù)容，計(jì)算各個(gè)時(shí)相庫(kù)容相比于最低水位時(shí)的庫(kù)容增量ΔV2。將庫(kù)容曲線計(jì)算得到的庫(kù)容增量ΔV1與ΔV2相比，得到相對(duì)誤差在10%~20%之間。

由于受水質(zhì)、混合像元、成像條件等因素的影響，水體光譜信息復(fù)雜多變，無(wú)法在缺少人工干預(yù)的情況下獲得較高精度的水體信息。下一步應(yīng)對(duì)大面積水域、復(fù)雜水體信息提取方法進(jìn)行研究，提高水體提取的自動(dòng)化程度和分類精度。此外，在不借助實(shí)地勘測(cè)或水文資料的條件下，利用遙感數(shù)據(jù)僅能實(shí)現(xiàn)蓄水量變化值的監(jiān)測(cè)。今后的研究可借助激光雷達(dá)技術(shù)對(duì)水下及周邊地形進(jìn)行勘測(cè)，在不借助實(shí)地勘測(cè)的條件下，實(shí)現(xiàn)對(duì)大范圍水域總蓄水量的監(jiān)測(cè)。

參考文獻(xiàn)

［1］ McFEETERS S K.The use of the normalized difference water index (NDWI) in the delineation of open water features［J］.International Journal of Remote Sensing,1996,17(7):1425-1432.

［2］ 徐涵秋.利用改進(jìn)的歸一化差異水體指數(shù)(MNDWI)提取水體信息的研究［J］.遙感學(xué)報(bào),2005,9(5):589-595.

［3］ 趙紫薇.基于OSTU算法利用新型水體指數(shù)進(jìn)行Landsat數(shù)據(jù)自適應(yīng)閾值水體自動(dòng)提取研究［J］.測(cè)繪與空間地理信息,2016,39(9):57-60.

［4］ 楊存建,周成虎.利用Radarsat SWA SAR和Landsat TM的互補(bǔ)信息確定洪水水體范圍［J］.自然災(zāi)害學(xué)報(bào),2001,10(2):79-83.

［5］ 汪金花,張永彬,孔改紅.譜間關(guān)系法在水體特征提取中的應(yīng)用［J］.礦山測(cè)量,2004,(4):30-32.

［6］ 劉桂林,張落成,劉劍, 等.基于Landsat TM影像的水體信息提?。跩］.中國(guó)科學(xué)院大學(xué)學(xué)報(bào),2013,30(5):644-650.

［7］ 江惠芳,金建樂(lè).衛(wèi)星遙感技術(shù)在緊水灘/石塘水庫(kù)水位庫(kù)容曲線復(fù)核中的應(yīng)用［J］.西北水電,2015,(5):21-24.

［8］ 陳曦,裴毅,姚幫松, 等.水位庫(kù)容曲線的衛(wèi)星影像測(cè)定方法研究［J］.人民長(zhǎng)江,2013,44(20):25-28.

［9］ 丁志雄.DEM與遙感相結(jié)合的水庫(kù)水位面積曲線測(cè)定方法研究［J］.水利水電技術(shù),2010,41(1):83-86.

［10］ 丁志雄,顏廷松,屈吉鴻.多源遙感影像在水庫(kù)水位-庫(kù)容曲線復(fù)核中的應(yīng)用［J］.華北水利水電大學(xué)學(xué)報(bào)(自然科學(xué)版),2012,33(4):32-35.

［11］ 陳曉玲,陸建忠,蔡曉斌,等.基于空間信息技術(shù)的堰塞湖庫(kù)容分析方法研究［J］.遙感學(xué)報(bào),2008,12(6):885-892.

［12］ 朱長(zhǎng)明,張新,路明,等.湖盆數(shù)據(jù)未知的湖泊動(dòng)態(tài)庫(kù)容遙感監(jiān)測(cè)方法［J］.測(cè)繪學(xué)報(bào),2015,44(3):309-315.

［13］ 劉東,李艷.基于遙感技術(shù)的鄱陽(yáng)湖面積庫(kù)容估算［J］.遙感信息,2012(2):57-61.

［14］ 李紀(jì)人,黃詩(shī)峰.“3S”技術(shù)水利應(yīng)用指南［M］.北京：中國(guó)水利水電出版社,2002.

［15］ 武國(guó)瑛,鄧正棟,陳一村.改進(jìn)的OTSU法在遙感水體信息提取中的應(yīng)用［J］.微型機(jī)與應(yīng)用,2017,36(20):17-18,22.

［16］ 張祖勛.數(shù)字?jǐn)z影測(cè)量學(xué)［M］.武漢：武漢大學(xué)出版社,2012.

（收稿日期：2018-10-20）

作者簡(jiǎn)介：

姜曉晨(1994-)，男，碩士研究生，主要研究方向：水深水質(zhì)遙感。

鄧正棟(1960-)，通信作者，男，博士，教授，主要研究方向：野戰(zhàn)給水保障理論。E-mail:dengzdong@sina.com。

武國(guó)瑛(1992-)，男，碩士研究生，主要研究方向：水深水質(zhì)遙感。