姜曉晨，鄧正棟，武國(guó)瑛

(陸軍工程大學(xué) 國(guó)防工程學(xué)院，江蘇 南京 210007)

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ù)的種類(lèi)不斷豐富，這些數(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)所示：

(1)

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

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

(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 水體提取結(jié)果

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。

圖2立體像對(duì)提取DEM原理圖

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

圖3 DEM提取流程

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

圖4 DEM效果對(duì)比

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

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

2.1.1Landsat8數(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所示。

表1 Landsat 8 OLI傳感器波段參數(shù)信息

本研究數(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，滿(mǎn)足1∶50 000比例尺立體測(cè)圖要求和1∶25 000比例尺地形圖更新需求。衛(wèi)星過(guò)境一次即可獲得三視立體影像和多光譜影像，影像可以組成同軌立體，重訪周期為5天。

表2 Landsat 8 OLI數(shù)據(jù)列表

本研究所采用數(shù)據(jù)購(gòu)買(mǎi)于遙感集市(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。

表3 資源三號(hào)衛(wèi)星數(shù)據(jù)列表

2.2 水體面積和完整DEM提取

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

表4 水庫(kù)水位-面積表

其中日期為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所示。

圖5 水位-面積曲線

圖6 研究區(qū)完整DEM

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。

圖7 水位-庫(kù)容曲線

本研究以水位-庫(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%之間。

表5 蓄水量變化結(jié)果對(duì)比

經(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)化程度和分類(lèi)精度。此外，在不借助實(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è)。