基于多源遥感技术的水库特征参数反演研究

李闯　李兵　王宁

摘要：以丹江口水库为例，采用Sentinel-1A合成孔径雷达（SAR）和Sentinel-2A/B多光谱成像仪（MSI）联合提取水域面积，通过建立多源传感器（SAR和MSI）所提取水域面积差异的定量函数关系，修正二者水域面积之间的差异，并将修正后的水域面积与Sentinel-3A雷达高度计（RA）所提取的水位相结合，开展水库蓄水量反演研究，最后以“全国水雨情网”公布的实测数据等为参考，评估水库特征参数的反演精度。结果表明：由RA数据提取的遥测水位与实测水位相比，R2达到0.99以上，RMSE达到0.093 3，遥测水位与实测水位的平均误差仅有9.33 cm;相对于单一MSI影像或SAR影像所提取的水域面积，二者水域面积联合参与建立的“遥测水位-多源遥感水域面积”模型的R2和RMSE均有所优化，其中，R2分別提高了0.26%和0.60%，RMSE分别降低了2.63和1.29;由遥测水位和多源遥感水域面积结合水库蓄水量模型反演的蓄水量与标准蓄水量相比，R2达到0.99，RMSE达到1.177，相对误差平均为1.8%。研究结果表明：利用多源遥感技术开展水库特征参数的反演研究，能够更加准确地反映水库的运行状态。

关键词：水库特征参数; 多源遥感; 雷达高度计; 多光谱成像仪; 合成孔径雷达

中图法分类号： TN958

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2021.11.015

0引 言

水库是陆地水资源的聚集地，水库特征参数主要包括水库的水位、水域面积和蓄水量等。这些参数直接反映着水库健康运行的状态，对区域的气候、湿度、降水等有着重要的影响。研究水库特征参数有助于掌握水库的运行状态，对拦截洪水、保障供水、调配水资源等具有积极作用，因此，开展水库特征参数的反演研究具有十分重要的现实意义[1]。

传统方法主要依靠人力到现场布设各种监测仪器并绘制库区地图等来反演水库特征参数;单一遥感方法主要利用一种遥感技术对水库的一种参数进行提取，其他参数仍采用现场调查的方式。李书杰[2]利用人力进行外业测量获得水库的水位并绘制了库区地图，通过对内业资料的整理，利用等高线、方格网等求得水库的水域面积，受制图过程中人为操作产生的精度损失影响，所反演的蓄水量只是粗略的结果。苏业助等[3]为了科学化管理水库，减少因外业测量和内业资料整理而造成的数据偏差，根据蓄水量平衡原理修正了蓄水量反演曲线，提高了蓄水量反演精度。孙建芸等[4]基于我国GF-1卫星数据和实测水文数据，利用归一化水体指数（NDWI）提取了丹江口水库的水域面积，并根据水位、水域面积和蓄水量三者之间的相关性，分析了丹江口水库的供水情况，是遥感数据与实测数据相结合的应用。王庆健[5]利用Sentinel-1A提取的水域面积、“全国水雨情网”公布的实测水位和标准蓄水量反演了水库参数，并通过蓄水量平衡公式计算水库的渗透量等，最后，对影响蓄水量的各种因素进行了分析。这些研究采集信息的时间长、工作量大、耗资多，所布设的仪器也常常会因为恶劣天气和日常消耗而损坏，从而导致数据的缺失;一段时期内，获取的可用遥感数据少，并且在相关模型建立之后仍需输入实测数据才能获得参数结果，对实测数据的依赖较大，精度相对较低。

相对于传统方法和单一遥感方法，多源遥感方法能够使各遥感技术优势互补，数据源丰富，水库特征参数反演精度高。本文利用星载的雷达高度计（Rader Altimetry，RA）提取水库的水位，以及星载的多光谱成像仪（Multi-Spectral Instrument，MSI）和合成孔径雷达（Synthetic Aperture Radar，SAR）联合提取水库的水域面积。提取水域面积利用的是卫星有效载荷生产的多源遥感数据，会存在一定的误差。误差来源主要有3个方面：① 不同卫星所搭载的传感器本身会有一定误差;② 不同卫星过境同一水库的时间不同，会增加以小时、天等为时间单位的计算误差;③ 不同卫星传感器的成像机理不同，所生产的数据需要采用特定的处理方法，这些处理方法精度不一，会造成一定的误差[6-7]。虽然可以采用轨迹间和卫星间等偏差的先验调整，但不同卫星的修正方法使得水库特征参数的计算变得繁琐与复杂，因此，综合考虑多源遥感技术在水库特征参数上的研究并不多[8]。本文不去深究修正上述3种误差的具体方法，而是在验证水域面积精度的基础上，通过建立差异函数直接对提取的水域面积结果进行修正，减少了中间过程中的复杂计算，节约了大量资源，为利用多源遥感技术开展水库特征参数的反演研究提供了有益参考。

1研究区域与数据资料

1.1研究区域

丹江口水库（32°36′～33°48′N、110°59′～111°49′E）位于汉江中上游，如图1所示，是国家一级水源保护区和国家级生态文明示范区，作为我国南水北调中线工程的水源地，保障了河南、河北、北京、天津四省（市）的生产和生活用水[9]。

丹江口水库从1958年9月开工，到1973年结束初期工程的建设，自建库以来发挥了巨大的经济、社会和生态效益。2013年8月完成对二期工程的验收，设计坝顶高程为176.6 m，正常水位为170 m，水域面积可达1 050 km2，蓄水量可达290.5亿m3，2014年12月正式向南水北调中线一期工程输水[10-11]。

1.2数据资料

采用欧洲航天局的Sentinel系列卫星作为遥感数据源，获取地址为：https：∥scihub.copernicus.eu/dhus/，采用“全国水雨情网”上公布的实测数据作为地面验证的数据源，获取地址为：http：∥xxfb.mwr.cn/，所有数据的获取日期为2019年1月1日至2019年12月31日，其中，Sentinel-1A SAR影像29景、Sentinel-2A/B MSI影像14景和Sentinel-3A RA数据13份等。

重访周期为12 d的Sentinel-1A卫星和重访周期为27 d的Sentinel-3A卫星是发射电磁波的主动工作模式，所搭载的合成孔径雷达，不受云雾和天气的干扰，具有全天时、全天候的特点。重访周期为5 d的Sentinel-2A/B卫星是反射光波的被动工作模式，所搭载的多光谱成像仪，可提供媲美人眼的光学数据，相对于灰度值的成像模式，其对地物目标的光谱差异性具有更加细微的捕捉能力[12-13]。

2研究方法

2.1水位反演

利用Sentinel-3A雷达高度计所生产的RA数据反演水库水位，理想情况下，RA测量结果为卫星到水面的瞬时距离[14]，如式（1）所示。

h=cτ2（1）

式中：h为RA到水面的高度;c为光速;τ为从RA发射信号到接收信号的往返时间。

RA在测量过程中会受许多因素的影响，例如大气、电离层、波浪等，实际计算时必须消除这些影响，主要有干对流层Rangedt、湿对流层Rangewt、电离层Rangei、固体潮Rangeset、极潮Rangept以及水面椭球高GeoEGM96的修正，如式（2）所示。受误差影响和参考水准面的选取，该文通过Ice-1波形重跟踪算法处理Sentinel-3A RA数据获得以中国黄海水准面为基准的遥测水位，如式（3）所示，Range Cor为总误差修正量，Height为水位，Altitude为卫星到地球椭球面距离，Range为卫星到水面距离。

2.2水域面积反演

2.2.1MSI影像反演水域面积

由MSI影像反演水域面积的处理流程如图2所示，包括影像预处理、水体信息增强、水陆分割、后处理，以及水域面积计算等步骤。

其中，水体信息增强是为了加强水体与非水体之间的对比度，以便提高后续处理的精度;水陆分割用于水体提取，将影像特征相似的非水体区分开来;后处理主要用于对分割結果中可能存在的干扰信息进行去除，例如桥梁、浮岛、池塘等地物的影响，并获得完整水体轮廓。

图3展示了Sentinel-2A/B卫星所生产的MSI影像反演水域面积过程中关键环节的处理结果。采用NDWI指数对图3（a）所示的MSI影像图以地物绿波段和近红外波段反射率的差异性增强水体信息[15]，处理结果如图3（b）所示。水体信息增强后的二值图中，水体呈现较低灰度值的黑色，非水体呈现较高灰度值的灰色，具有良好的区分度。水陆分割利用决策树将具有相似性的水体与地形阴影进行区分，结果如图3（c）所示。图3（c）中可以明显看出存在非水体的干扰信息，可以通过设置连通域阈值方式对小连通域干扰物进行去除，从而获得水体轮廓的提取结果，如图3（d）所示。

基于提取到的水体，采用如式（4）所示的像元统计方法进行计算，可获得水域面积Smw。

式中：M为影像中像元总数目;N为提取到的水体像元数目;Sim为影像所对应的地表区域面积，km2;Smw为最终反演的水域面积，km2。

2.2.2SAR影像反演水域面积

SAR影像反演水域面积的处理流程如图4所示。首先要进行辐射定标、多视处理、斑点滤波、几何校正等处理，以提高影像质量，为获取更加精准水体信息奠定基础[16-17]。

预处理后的Sentinel-1A SAR影像如图5（a）所示。由细节放大图可知，该影像具有较高分辨率。水陆分割是一个区分水体与非水体的二分类过程，模糊C均值算法（Fuzzy C-Means，FCM）是一种常用分类方法，该方法以随机聚类中心开始，经过不断的迭代优化完成地物类型归属划分[18]。经水陆分割和小连通域去除的水体轮廓提取结果如图5（b）所示。

采用FCM算法进行水陆分割的关键是确定影像中水体与非水体的临界灰度值信息，从而获得更加准确的聚类中心。因此，该文对SAR影像像元的灰度值进行了统计，其直方图如图5（c）所示。

可以明显地看到直方图中出现了2个波峰，第一个波峰主要为SAR影像中的水体，其像元灰度值分布在0～8内，第二个波峰以及剩余灰度值主要代表SAR影像中的非水体。灰度值8为水体与非水体的临界值，为了避免漏提浅水水体，结合灰度统计直方图和实践研究，应将FCM算法的一个关键聚类中心设置为8。

2.3多源遥感水域面积修正研究

MSI影像的观测周期较短，影像信息丰富，但受天气和光照的影响较大，故MSI影像实际可用的数据较少。SAR影像可全天时、全天候工作，但其重返周期较长，加上受限于军事保密、技术瓶颈和成本控制，获取到合适的SAR影像数据也相对较少。针对上述情况，将MSI影像和SAR影像联合应用于研究，可以丰富数据源、提高实验精度。

MSI影像和SAR影像属于异源传感器影像，前者是被动成像模式，利用光波的反射成像，后者是主动成像模式，利用电磁波的散射成像，两者具有不同的成像机理，出现了“同物异谱”现象。为了获得后续更精确的蓄水量反演结果，必须构建两者水域面积差异的定量函数关系，修正上述误差。

在构建定量函数关系之前，需要先验证水域面积的提取精度。水体轮廓对水域面积的提取精度有着至关重要的影响，只要验证水体轮廓的提取精度即可说明水域面积的提取精度。利用遥感数据提取水体轮廓，手动提取法精度最高，但其效率低下。以手动提取的水体轮廓为验证数据，通过漏警率、虚警率、准确率和Kappa系数等指标，对提取的水体轮廓进行精度验证，其定量分析结果如表1所列。可知，2种影像水体轮廓的漏警率和虚警率都较低，准确率优于0.97，Kappa系数高于0.91，表明提取的水体轮廓足够用于后续研究。

对SAR水域面积和MSI水域面积进行指数、幂次、多项式等不同函数的拟合，并选取最优函数来建立多源遥感水域面积差异的定量函数关系，所得定量分析的结果如表2所列。

根据R2越大越好、RMSE越小越好和多项式次数越低计算越简单的定量分析原则，四次多项式函数具有最小的RMSE、较大的R2和较低的多项式次数。因此，采用四次多项式函数表示MSI影像和SAR影像所提取水域面积之间的差异。

四次多项式函数的表达式如式（5）所示。

F（x）=ax4+bx3+cx2+dx+e （5）

式中：x为水位，m;F（x）为同一水位下MSI影像和SAR影像提取水域面积的差异值，m2;a、b、c、d、e为多项式系数，拟合后的值分别为-0.000 476，0.302 6，-72.1，7 632，-0.000 030 28。

丹江口水库位于北半球和南半球冷暖气流的交汇处，易受副热带高压北进摆动和稳定南退的影响，天气变化剧烈，水库上空云雾天气的概率较大。对于选取的部分MSI影像，不可避免地存在少量的云雾，其遮挡效应在一定程度上影响了水域面积提取的精度。相比之下，SAR影像具有全天时、全天候的特点，不存在上述误差。另外，从水体轮廓提取结果来说，SAR影像所提取水体轮廓精度要优于MSI影像。因此，考虑到上述情况，该文以SAR影像为基准，修正水域面积差异，即将MSI影像提取的水域面积值P（x），加上面积误差修正值F（x），得到修正后的MSI水域面积S（x），如式（6）所示。修正后的MSI水域面积与SAR水域面积联合在一起称为多源遥感水域面积。

S（x）=P（x）+F（x）（6）

2.4水库蓄水量反演

水库蓄水量是无法直接观测的量，一般通过水位和水域面积进行间接反演。该文把水库理想化为棱台，获得邻水位或者相邻横断面之间的一小部分蓄水量，并把这些小蓄水量累加起来，再加上死水位对应的蓄水量就能求得整个水库的总蓄水量[19]，如式（7）～（9）所示。

3结果与分析

3.1水位结果及分析

根据Sentinel-3A RA数据提取遥测水位的获取日期查找到对应的实测水位，然后将遥测水位和实测水位进行对比，结果如图6所示，日期范围是2019年1月1日至12月31日。

由图6可知，2019年丹江口水库的水位变化范围为150～167m。随着旱季和汛季的转换，水位呈现规律性波动，在7～10月的汛期，水位上涨迅速，11月之后将逐渐进入枯水期，水位开始慢慢下降，并且由Sentinel-3A RA提取的遥测水位与实测水位的波动趋势相吻合。

采用R2和RMSE对遥测水位的精度进行评价，得到R2=0.99，RMSE=0.0933，表明提取的遥测水位具有很高的精度，遥测水位与实测水位的平均误差仅有9.33 cm。

3.2水域面积结果及分析

利用星载MSI影像和星载SAR影像分别提取的MSI水域面积和SAR水域面积，结果如图7～8所示。

由图7～8可知：2019年丹江口水库的水域面积在5月初达到最小值，7月开始明显增大，在11月中旬达到最大值，然后再次回落，最大水域面积与最小水域面积相差约250 km2。随着时间推移，整体上MSI水域面积和SAR水域面积变化趋势较为吻合，MSI水域面积和SAR水域面积的准确率分别达到了97.95%和98.95%。

3.3建立“水位-水域面积”模型

由2.4小节可知：只要获得水位和水域面积就能在不依靠地面实测数据的情况下获得蓄水量。对于水位和水域面积，已知其中一个参数，想要快速获得另外一个参数，需要建立“水位-水域面积”关系曲线。该关系曲线建立之后，再通过多源遥感技术获得水位或水域面积任何一个参数，就能获得包括蓄水量在内的另外两个参数。

以多源遥感技术获取的遥测水位、MSI水域面积、SAR水域面积以及修正后的多源遥感水域面积分别建立“遥测水位-MSI水域面积”关系曲线、“遥测水位-SAR水域面积”

关系曲线和“遥测水位-多源遥感水域面积”关系曲线。以“遥测水位-多源遥感水域面积”关系曲线为例，建立过程如下：

首先构建“遥测水位-多源遥感水域面积”关系曲线的指数函数、幂函数和多项式函数的数学模型，由于四次及四次以上多项式函数出现过拟合现象，选择指数函数、冪函数和四次以下多项式函数对拟合效果作进一步验证，如表3所列。

由表3可知：指数函数的R2最小，RMSE最大;幂函数与指数函数的R2和RMSE最接近;三次多项式函数的R2最大，但RMSE不是最小，其R2只比二次多项式函数的R2多0.0003，二次多项式函数的RMSE最小。根据R2越大越好，RMSE越小越好的定量分析原则，二次多项式函数最适合描述遥测水位与多源遥感水域面积之间的关系。

得到“遥测水位-多源遥感水域面积”模型如式（10）所示，关系曲线如图9（a）所示;同理可得“遥测水位-MSI水域面积”模型如式（11）所示，关系曲线如图9（b）所示;“遥测水位-SAR水域面积”模型如式（12）所示，关系曲线如图9（c）所示。

“水位-水域面积”关系曲线建立之后，各模型最优函数的R2和RMSE如表4所列。

由表4可知：各模型本身也已经具有较高的R2和较低的RMSE，相比单独使用MSI影像或SAR影像，联合使用两种遥感影像后，“遥测水位-多源遥感水域面积”模型的R2分别提高了0.26%和0.60%，RMSE分别降低了2.63和1.29，即R2和RMSE均有所优化。

3.4蓄水量结果及分析

利用多源遥感技术反演得到水库的遥测水位和多源遥感水域面积，进而得到水库的蓄水量，将反演的蓄水量与“全国水雨情网”上公布的标准蓄水量进行对比，对比结果如图10所示，日期范围是2019年1月1日至12月31日。

由图10可知，反演蓄水量和标准蓄水量的波动趋势具有较高吻合度。在7～10月的汛期，相关地区降雨较多，丹江口水库的蓄水量增长迅速，在11月17日左右达到最高峰，然后开始慢慢减少，逐渐进入枯水期，在5月初达到最低值。整体上，随着旱季和汛季的转换，其蓄水量也呈现规律性的变化。

采用R2和RMSE进行定量分析。分析结果表明，R2=0.99，说明反演蓄水量与标准蓄水量具有极高的吻合度，能够真实地反映水库蓄水量的变化情况;另外，RMSE=1.177，表明反演的水库蓄水量标准差在1.177亿m3以内，相对误差平均为1.8%。

4结 语

本文以丹江口水库为例，利用多源遥感技术反演了水库参数，其中，利用RA数据提取了水库的水位，利用多源传感器生产的MSI影像和SAR影像提取水库的水域面积。通过构建多源传感器水域面积之间的定量函数关系，修正了二者水域面积之间的误差，最后，基于修正后多源遥感水域面积，结合蓄水量模型反演水库的蓄水量。

实验结果表明，利用多源遥感技术反演的水库特征参数能够真实地反映水库健康运行的状态。在水库特征参数反演过程中，除了利用地面实测数据对反演结果进行验证之外，不再依赖实测数据的支持就能达到对水库特征参数的反演，降低了成本，节约了资源，有效补充了实测数据缺失情况下的研究工作。同时，利用多源遥感技术反演水库特征参数能够达到对水库异常预警和日常管理的最大效益化，特别是对堰塞湖或出现大面积违法的围湖造田等的预警。该文关于水库特征参数的反演方法不仅能够克服当前传统方法和单一遥感方法的不足，也可以进一步推广到其它大型水库特征参数的反演工作当中。未来还可以考虑采用更多类型的遥感数据和更加优良的水库蓄水量模型等参与反演。

参考文献：

[1]周建中，顿晓晗，张勇传.基于库容风险频率曲线的水库群联合防洪调度研究[J].水利学报，2019，50（11）：1318-1325.

[2]李书杰.水库库容的量算精度[J].东北水利水电，1997（7）：34-36.

[3]苏业助，洪为善.根据水量平衡原理修正水库库容曲线方法[J].人民长江，1999，30（12）：39-41.

[4]孙建芸，袁琳，王新生，等.基于GF_1卫星的丹江口水库水面面积-蓄水量-水位相关性研究[J].南水北调与水利科技，2017，15（5）：89-96.

[5]王庆健.基于SAR监测的水库主要参数计算系统的设计与实现[D].开封：河南大学，2018.

[6]BONNEFOND P，EXERTIER P，LAURAIN O，et al.Absolute calibration of Jason-1 and Jason-2 altimeters in Corsica during the for- mation flight phase[J].Marine Geodesy，2010，33（S1）：80-90.

[7]MERTIKAS S P，IOANNIDES R T，TZIAVOS I N，et al.Statistical models and latest results in the determination of the absolute bias for the radar altimeters of Jason satellites using the gavdos facility[J].Marine Geodesy，2010，33（S1） ：114-149.

[8]王文种，黄对，刘九夫，等.基于Landsat与Sentinel-3A卫星数据的当惹雍错1988～2018年湖泊水位-水量变化及归因[J].湖泊科学，2020，32（5）：1552-1563.

[9]王元超，王旭，雷晓辉，等.丹江口水库入库径流特征及其演变规律[J].南水北调与水利科技，2015，13（1）：15-19.

[10]吴川，张玉龙，许秀贞，等.基于Landset TM/ETM和HJ-1A/B影像的丹江口水库水域变化监测研究[J].长江流域资源与环境，2013，22（9）：108-114.

[11]尹志杰，王容，赵兰兰，等.丹江口水库以上流域夏秋汛期来水及相关性分析[J].水文，2019，39（4）：74-77.

[12]杨魁，杨建兵，江冰茹.Sentinel-1卫星综述[J].城市勘测，2015（2）：26-29.

[13]田颖，陈卓奇，惠凤鸣，等.欧空局哨兵卫星Sentinel-2A/B数据特征及应用前景分析[J].北京师范大学学报（自然科学版），2019，55（1）：57-65.

[14]郭华东.雷达对地观测理论与应用[M].北京：科学出版社，2000.

[15]GAO B C.NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space[J].Remote Sensing of Environment，1995，58（3）：257-266.

[16]潘習哲.星载SAR图像处理[M].北京：科学出版社，1996.

[17]王超.SAR图像处理新进展[J].中国图象图形学报，2009，14（1）：1-2.

[18]BEZDEK J C，EHRLICH R，FULL W.FCM：The fuzzy c-means clustering algorithm[J].Computers & Geoences，1984，10（2-3）：191-203.

[19]张国庆.青藏高原湖泊变化遥感监测及其对气候变化的响应研究进展[J].地理科学进展，2018，（2）：214-223.

（编辑：黄文晋）

Abstract：Take Danjiangkou Reservoir as an example，the water area was extracted by using the Sentinel-1A Synthetic Aperture Radar （SAR） and Sentinel-2A/B Multi-Spectral Imager （MSI），the difference of water areas by SAR and by MSI was corrected by establishing a quantitative function relation of these areas，and the revised water area was combined with the inverted water level of Sentinel-3A Radar Altimeter （RA） to calculate the reservoir volume.Finally，the retrieval precision of the reservoir parameters was evaluated according to the measured data published by the National Water and rain regime network.The results showed that R2 was above 0.99 and RMSE was 0.0933，which meant that the average error between inverted water level and measured water level was only 9.33 cm.Compared with the water area extracted from single MSI image or SAR image，the R2 and RMSE of the "remote sensing water level-Multi-Source Remote Sensing water area" model were optimized，especially R2 increased by 0.26% and 0.60% ，while RMSE decreased by 2.63 and 1.29，respectively.Compared with the standard water storage，for the novel model，R2 reached 0.99 and RMSE reached 1.177，respectively，and the relative error was 1.8% on average.The result shows that the reservoir characteristic parameters inversion using multi-source remote sensing technology can more accurately reflect the operation of the reservoir.

Key words：reservoir characteristic parameters;multi-source remote sensing;Rader Altimetry;Multi-Spectral Instrument;Synthetic Aperture Radar