基于資源三號衛星與Landsat 8 OLI的水庫庫容估算*

姜曉晨，鄧正棟，武國瑛

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

0 引言

水庫和湖泊是地表水體的主要存在形式，對河流下游地區調蓄洪水、維持水沙平衡具有重要作用。湖泊及水庫蓄水量的研究一直是防洪調蓄研究的重點，近年來遙感技術和地理信息系統技術越來越多地被國內外的研究人員應用于湖泊、水庫監測研究中。

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

目前，結合遙感數據監測水體庫容主要有以下兩種思路。一是用于庫容曲線的校核[7-11]；二是用于三維分析，直接測算水體庫容[12-14]。本文借鑒前述成果，利用指數法從Landsat 8 OLI影像中提取水體，結合資源三號衛星獲取的DEM數據實現對遼寧省撫順市大伙房水庫動態蓄水量的反演，并與水文資料進行對比。

1 反演原理和方法

1.1 水體提取

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

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

(1)

式中，Green，NIR分別對應Landsat 8 OLI影像的Band3、Band5。

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

(2)

式中，MIR為中紅外波段，對應Landsat 8 OLI影像的Band6。

在水體提取中，若采用同一閾值，會造成較大錯分、漏分現象。例如NDWI和MNDWI指數法，一般采用0作為閾值，但在實際影像提取中，須采用較大或較小的閾值才能獲得理想的效果。研究中采用MNDWI指數和紅外波段進行水體提取，利用OTSU法確定閾值[15]。圖1顯示的是本研究提取的水體區域。

圖1 水體提取結果

1.2 立體像對提取

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

圖2立體像對提取DEM原理圖

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

圖3 DEM提取流程

利用該模塊進行DEM提取時，需提供研究區的最大、最小高程進行精度控制。本研究利用ASTGTM2數據進行測算。圖4為DEM的效果對比圖。

圖4 DEM效果對比

2 水庫庫容反演

2.1 數據選擇與處理

2.1.1Landsat8數據

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

表1 Landsat 8 OLI傳感器波段參數信息

本研究數據來源于中國科學院計算機網絡信息中心地理空間數據云平臺(http://www.gscloud.cn)，共下載了27幅影像，條帶號119/31，所下載影像的成像時間涉及從2013年至2016年，研究區內的云層覆蓋率為5%以下，影像質量較好。數據列表如表2所示。

2.1.2資源三號數據

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

表2 Landsat 8 OLI數據列表

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

表3 資源三號衛星數據列表

2.2 水體面積和完整DEM提取

經過ENVI和Arcgis處理后得到的水體面積如表4所示。

表4 水庫水位-面積表

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

對篩選后的水體面積和相應水位數據進行曲線擬合，結果如圖5所示。

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

圖5 水位-面積曲線

圖6 研究區完整DEM

2.3 反演結果驗證

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

利用從大伙房水庫管理局收集的水位-庫容值，對其按三次關系進行擬合，可以建立水位-庫容曲線，見圖7。

圖7 水位-庫容曲線

本研究以水位-庫容曲線為依據，計算各個時相庫容相比于2014/06/04成像時的庫容增量ΔV1。將利用庫容曲線計算的庫容變化值和利用本研究方法計算的蓄水量變化值進行比較，如表5所示。通過計算二者差值發現，利用該方法計算的蓄水量變化值ΔV2與利用庫容曲線計算的庫容變化量ΔV1差值占ΔV1的10%～20%之間。

表5 蓄水量變化結果對比

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

3 結論

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

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

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