基于资源三号卫星与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%之间。

由于受水质、混合像元、成像条件等因素的影响，水体光谱信息复杂多变，无法在缺少人工干预的情况下获得较高精度的水体信息。下一步应对大面积水域、复杂水体信息提取方法进行研究，提高水体提取的自动化程度和分类精度。此外，在不借助实地勘测或水文资料的条件下，利用遥感数据仅能实现蓄水量变化值的监测。今后的研究可借助激光雷达技术对水下及周边地形进行勘测，在不借助实地勘测的条件下，实现对大范围水域总蓄水量的监测。