近年来, 由于气候变化, 极端天气增多, 洪灾仍是给人类社会带来损失最大的自然灾害^[1]. 水库作为一种重要的水利枢纽, 对拦截洪水、保障供水、调配水资源等工作具有积极作用, 也为社会经济发展和生态文明建设提供着不可忽略的助力. 水库蓄水量是水库健康运行的重要指标之一, 与区域的气候、湿度、降水等密切相关. 因此, 研究水库蓄水量对合理保护和利用水资源以及制定相关政策十分必要^[2].

传统的水库蓄水量测量方法主要依靠人力到现场布设各种监测仪器来获取数据^[3], 通过绘制库区地图并将其分割成等高线、三角网等, 进行水域面积的提取, 然后结合等高距, 计算水库蓄水量. 传统方法不仅耗费大量的人力、物力和财力, 而且水域面积的测量精度受制图精度的影响较大, 所利用的仪器也常常会因为天气恶劣和日常消耗而损坏, 从而导致数据的缺失. 相对于传统方法, 遥感方法具有监测范围广、收集信息快和工作平台多样化等优势. 单一遥感方法主要利用遥感技术对水库的一种参数进行提取^[4], 其他参数则采用现场调查的方式获取, 这种方法研究手段单一, 在同一时期获取的可用遥感数据少, 精度相对较低. 而对于我国国土面积辽阔、水库分布零散的现状而言, 并不是所有的水库都适合到现场开展调查. 例如青藏高原, 由于海拔高、气候恶劣、地形地貌复杂等原因, 会在一定时期造成人类难以到达的无人区, 使之不宜开展现场调查. 采用多源遥感能够丰富数据源, 使不同遥感技术优势互补, 获取比单一遥感更好的监测效果. 因此, 基于多源遥感信息进行水库蓄水量估算研究具有十分重要的现实意义.

目前已有许多基于多源遥感信息的水资源监测研究. 2016年, 文献[5]基于Envisat、Jason-2和Terra/Aqua MODIS等多源卫星遥感数据, 获取了河流水位和河宽数据, 利用改进的曼宁公式估算了亚马逊和尼罗河的径流量. 2017年, 文献[6]利用合成孔径雷达(Synthetic Aperture Radar, SAR)影像、光学遥感影像和机载激光雷达(Light Detection And Ranging, LiDAR)数据, 建立了3个数据集的融合分类模型, 采用多级决策树, 通过最小化单个数据集的模型之间的差异实现水体的高精度提取. 2019年, 文献[7]综合运用卫星测高技术和光学遥感岸线监测技术, 研制了高时间分辨率的2000年−2017年青藏高原湖泊水位、水量变化数据集. 多源遥感卫星技术已在地表水资源监测研究中得到广泛应用, 为地表水资源监测提供了可靠的理论和技术支撑, 但目前基于多源遥感信息的水库蓄水量估算研究非常少见, 基于多源遥感信息的水库蓄水量估算系统更是鲜有报道.

近年来, 已有研究人员开发了水库防洪调度^[8,9]、库容计算^[10]、水库地质灾害视频监测预警^[11,12]等多种水库信息管理系统, 但是这些系统所利用的数据多是地面实测数据或是第三方提供的数据, 很少利用多源遥感技术直接在软件系统中提取水库参数. 2018年, 文献[13]设计和实现了一款基于SAR遥感监测的水库参数计算系统, 利用从SAR影像提取的水域面积、从国家水利部“全国水雨情信息网”获取的实测水位和标准蓄水量等数据, 通过数据同化, 对丹江口水库的面积和蓄水量等参数进行了估算研究, 但该文的蓄水量估算模型仅利用SAR遥感影像自动提取水域面积参数, 模型所需的实测水位等关键参数仍需从第三方网站获取, 对实测数据的依赖较大, 而且受限于网络传输速率等方面的限制, 一旦缺乏数据支持(如重大灾害导致数据网络中断等), 该模型和系统将无法独立工作. 因此, 开发出一款具有独立工作能力的多源遥感水库蓄水量估算系统具有十分重要的研究意义和应用价值.

针对上述问题, 本文利用星载雷达高度计(Rader Altimetry, RA)^[14]提取水库水位, 利用优势互补的星载多光谱成像仪(Multi-Spectral Instrument, MSI)影像和星载SAR影像提取水库的水域面积, 采用数据同化法建立了多源遥感水域面积差异的定量函数关系, 构建了多源遥感水库蓄水量估算模型, 并基于该模型, 设计和实现了一款多源遥感水库蓄水量估算系统.

1 基础理论与相关技术

本系统在水库参数提取中, 基于多源遥感技术获取水位和水域面积等参数, 利用数据同化法建立了多源遥感水库蓄水量估算模型, 其技术路线如图1所示. 在系统实现上, 采用Matlab和Java混合编程技术来完成数据处理功能; 利用虚拟水文站点技术使水位以点标形式呈现在百度地图API上; 使用ECharts技术将系统中数字化信息转化成可视化图表信息.

1.1 多源遥感提取水库参数 1.1.1 提取遥测水位

系统利用RA数据提取水库的遥测水位. 提取过程中, 要进行一系列的误差校正, 包括干对流层Range_dt、湿对流层Range_wt、电离层Range_i、固体潮Range_set和极潮Range_pt等测距误差校正, 之后需要进行水面椭球高Geo_EGM96转化, 使遥测水位与实测水位一样, 都以中国黄海水准面为基准, 如式(1)、式(2)所示, 其中H为遥测水位, Altitude为卫星到地球椭球面距离, Range为卫星到水面距离, Range_Corrected为误差校正距离.

$ H = Altitude - Range - Rang{e_{\rm Corrected}} $

(1)

$ \begin{split} Rang{e_{\rm Corrected}} =& Rang{e_{dt}} + Rang{e_{wt}} + Rang{e_i}\\ & + Rang{e_{\rm set}} + Rang{e_{pt}} + Ge{o_{\rm EGM96}} \end{split} $

(2)

1.1.2 提取水域面积

系统使用MSI影像和SAR影像作为获取水域面积的数据源, 对多源遥感数据进行预处理, 利用水陆分割算法提取水库的水体轮廓, 根据水体轮廓计算面积^[15]. 水域面积提取原理为: 先统计整幅影像中总像元个数M, 再统计水库水域面积的像元个数N, 设整幅影像对应到现实地面中的实际大小为S, 则实际水域面积P的计算公式如式(3).

$ P = \frac{{NS}}{M} $

(3)

1.1.3 估算蓄水量

多源遥感水域面积提取后, 对其精度进行定量分析, 基于数据同化技术, 构建两者水域面积差异的定量函数关系, 以SAR影像为基准, 修正水域面积差异, 构建“遥测水位-多源遥感水域面积”模型和“多源遥感水域面积-标准蓄水量”模型, 然后将上述两种模型进行数学转化, 建立本系统所使用的多源遥感水库蓄水量估算模型, 如式(4)所示, 其中Impoundage为估算蓄水量, H为遥测水位.

$ \begin{split} Impoundage =& 0.000\;024\;052\;74{H^4} - 0.019\;984\;11{H^3} \\ & + 6.030\;527\;67{H^2} - 780.828\;29H \\ & + 36\;821.0163 \\ \end{split}\!\!\!\!\!\!\!\! $

(4)

在上述过程中, 标准蓄水量用于辅助建模, 模型建成后, 在模型使用过程中无需再使用标准蓄水量, 克服了当前蓄水量估算模型建立之后仍需输入实测数据的不足, 在一定程度上完善了实测数据缺失情况下的研究工作, 具有一定现实意义.

1.2 混合编程

系统采用C/S (Client/Server, 客户端-服务器)架构, 使用Matlab和Java混合编程技术来完成数据处理功能. 利用Matlab强大的影像处理功能, 进行遥感影像水陆分割算法以及RA提取遥测水位算法的编程实现. 利用Java统计水陆分割之后水库水体轮廓的像素信息, 进而获得具体水域面积的值. 在Java与Matlab混合编程过程中, 先将各种处理算法的Matlab程序打包, 通过设置调用方法的Class类名导入Java工程中, 然后配置MCR (Matlab Compiler Runtime, Matlab编译运行)环境, 并安装MCR插件, 即可完成相关的调用. 系统调用框架图如图2所示.

1.3 虚拟水文站点

处理RA原始数据时, 仅根据处理之后遥测数据的经纬度数值, 难以判断所提取位置是否在水库水体上. 系统利用百度地图API抽取RA数据中每一个数值数据的经纬度, 以点标形式将其显示在百度地图上, 从而达到虚拟水文站的效果. 百度地图API使用独立的坐标体系, 因此在地图上使用经纬度坐标前还需要调用百度地图API接口进行坐标转换, 之后就可以以点标形式显示在地图上, 如图3所示, 用户能够直观看到RA数值数据的位置信息, 并可以通过点击水库水体上的各个点标, 来查看相应的水位值.

1.4 ECharts数据可视化

ECharts是一个使用JavaScript实现的开源可视化库, 不仅提供了折线图、雷达图、和弦图等各种常用的图表类型来满足不同用户的需求, 还提供了一些相对专业的图表来满足专业性较高的要求. 本系统中ECharts可视化处理流程如图4所示.

2 系统设计

本文从系统框架、流程以及数据库等3个方面对系统的总体设计进行分析与介绍.

2.1 系统框架设计

系统分为数据管理、数据处理和结果展示等3个功能模块, 模块结构如图5所示.

(1)数据管理模块

本模块主要包括水文数据管理、遥感数据管理和模型管理等功能. 水文数据管理主要是记录水文站的实测数据; 遥感数据管理主要进行RA数据、MSI影像数据和SAR影像数据的上传、检索以及处理后的下载与删除等操作; 模型管理主要实现系统内各算法和模型的增加、删除和修改等功能.

(2)数据处理模块

本模块主要实现各种数据处理功能, 包括利用RA提取算法提取遥测水位; 利用水陆分割算法对预处理后的SAR影像和MSI影像进行水陆分割, 从而得到水体轮廓; 调用面积统计算法计算具体的水域面积值; 结合遥测水位、水域面积、标准蓄水量, 建立蓄水量估算模型, 并估算水库蓄水量.

(3)结果展示模块

本模块主要实现水库水位、水域面积和蓄水量等结果的展示功能, 利用建好的相关数据库, 实现水位、水域面积和蓄水量之间的相互匹配, 并以曲线图和表格等方式予以展示.

2.2 系统流程设计

系统业务流程如图6所示, 通过数据输入、数据处理和结果展示等3大步骤, 实现基于多源遥感数据的水库蓄水量估算功能.

(1)数据输入

该步骤用于输入模型处理所需的各种数据, 主要包括实测数据、RA数据和遥感影像数据等3种类型.

实测数据: 用户上传每日记录的Excel文件, 或者直接在系统里利用爬虫技术爬取互联网上的相关数据, 包括实测水位和蓄水量.

RA数据: 从网上下载原始的RA数据, 经过解压之后存储到数据库里.

遥感影像数据: 包括预处理后的MSI影像和SAR影像, 上传时需要填写时间、卫星名称、所属水库等相关信息.

(2)数据处理

该步骤通过调用各种算法来计算水位、水域面积、蓄水量等水库参数, 是水库蓄水量估算模型的主要实现步骤.

水位计算: 遥测水位是调用RA提取算法处理RA原始数据而得到的水位, 在百度地图中显示了遥测水位对应在RA卫星星下点的位置, 方便用户选取合适的虚拟水文站点.

水域面积计算: 输入系统的SAR影像已进行斑点滤波和辐射定标等预处理操作, 在系统中只需调用水陆分割算法区分水体和非水体就能完成水库轮廓提取; 输入系统的MSI影像已进行水体信息增强、水体与非水体区分等预处理步骤, 在系统中仅需进行连通域识别, 以去除非水体小连通域. 通过以上步骤获得MSI影像和SAR影像中的水体轮廓后, 再调用面积统计算法, 计算水库水域面积值.

蓄水量估算: 基于数据同化技术, 将遥测水位、水域面积和标准蓄水量输入蓄水量估算模型, 估算水库蓄水量, 并绘制多种过程曲线图和表格来展示结果.

(3)结果展示

该步骤主要以曲线图和表格等方式展示各种处理结果. 系统里的过程曲线主要有“水位-日期”、“水位-水域面积”和“水位-蓄水量”等过程曲线, 其中, 当用户点击“水位-日期”过程曲线图时, 在“水位-水域面积”和“水位-蓄水量”等过程曲线图里相应的标记会随之改变, 在结果表格中会显示具体的数值, 包括日期、实测水位、遥测水位、水域面积和蓄水量等.

2.3 系统数据库设计

系统数据库用于存放系统处理过程中所涉及到的各种数据. 该系统的实体-联系图(Entity Relationship diagram, E-R图)如图7所示. 根据系统E-R图, 共设计了水库信息表、实测水位结果表、RA原始数据表、遥测水位结果表、MSI影像信息表、SAR影像信息表、影像处理算法表、水域面积表、蓄水量模型表、水库蓄水量表等共计10个数据表.

3 系统实现

系统采用C/S架构, 客户端提供数据管理、功能选择、结果显示等用户交互界面, 服务器端提供数据处理、参数计算等后台响应功能.

3.1 数据管理模块实现

数据管理模块主要对系统中各种数据进行管理, 实现实测数据、遥测水位、水域面积、影像管理、模型管理、水库管理等数据的上传、查询、修改、下载和删除等功能. 在图8所示的系统首页中, 点击左侧的“数据管理”按钮, 即可进入图9所示的数据管理页面, 对系统数据进行各种管理操作. 其中, SAR影像的上传功能如图10所示, 用户填好相关信息之后, 选择SAR影像路径或直接把SAR影像文件拖拽到上传控件, 点击“下一步”按钮, 客户端即可调用文件传输接口, 以流式文件方式向服务器端上传SAR影像文件. 影像下载功能实现方式类似, 在图9数据管理页面中点击右侧“下载”按钮后, 客户端会向服务器端发送SAR影像记录ID, 服务器端收到请求后会将SAR影像流式文件传输至客户端, 并弹出如图11所示的下载任务确认窗口, 用户可以选择存放SAR影像的路径, 点击“下载”按钮后, 完成SAR影像的下载.

3.2 数据处理模块实现

数据处理模块主要包含水位提取、水域面积提取、蓄水量估算等3个功能. 系统通过调用RA提取算法的Matlab包来处理RA的原始数据进而提取水库的遥测水位; 调用水陆分割算法的Matlab包实现对遥感影像中水体轮廓的提取, 然后调用面积统计算法的Java程序对水体轮廓所代表的水域面积进行计算; 最后结合数据同化算法估算出水库蓄水量.

(1)水位提取

系统中使用的水位数据主要有两种, 一种是水文站的实测水位, 另一种是从RA数据中提取的遥测水位, 其中实测水位的作用是对遥测水位进行验证. 对于水文站实测水位, 系统提供了两种获取方式, 一种是人工导入方式, 由用户每天登录国家水利部“全国水雨情信息网”进行记录, 并以Excel文件形式人工导入系统; 另一种是自动爬取方式, 系统利用爬虫技术自动从“全国水雨情信息网”爬取水文站实测水位及官方公布的标准蓄水量, 并自动存入系统数据库. 图12是利用爬虫技术爬取实测水位的开始页面, 用户点击图中右上角“获取今日实测数据”按钮后, 系统将自动爬取相关数据并直接存入系统数据库.

(2)水域面积提取

SAR影像和MSI影像的水域面积提取过程类似, 这里以SAR影像提取水域面积过程为例予以说明. 首先, 设置调用方法的class类名与Program()方法名, 将.m程序打包成.jar文件, 与JavaBuilder.jar文件一起导入Java工程. 然后, 调用接口访问上传至系统的SAR影像, 若影像为空, 则用户需再次提交SAR影像; 非空时, 先随机初始化模糊矩阵并根据SAR影像建立目标函数, 然后根据目标函数建立聚类中心, 之后更新模糊矩阵并分割出水体轮廓, 再统计水体轮廓和整个SAR影像像素点个数. 最后, 利用经纬度信息计算出SAR影像水域面积. 水陆分割后的SAR影像数据示例如图13所示, SAR影像水域面积提取结果示例如图14所示.

(3)蓄水量估算

基于不同数据源数据的不同组合方式, 系统提供了多种面积提取模型和蓄水量估算模型, 方便用户针对目标水库实际情况对比分析和择优选出精度更高的蓄水量估算模型. 对于系统提供的各种模型, 用户可以进行添加、修改和删除等编辑操作, 模型编辑界面如图15所示. 用户对编辑界面中的各模型参数进行调整并保存后, 系统会使用保存的各模型参数进行后续计算. 各模型计算结果对比示例如图16所示.

3.3 结果展示模块实现

结果展示模块调用数据库中的数据并以图表形式将其绘制在系统页面上, 实现“水位-日期”、“水位-水域面积”和“水位-蓄水量”等数据的可视化展示. 系统读取数据库里的数据后, 利用ECharts可视化技术将数字化信息转化成可视化的图表信息, 最后把转化后的图表信息呈现在系统的前端页面上, 如图17所示.

4 系统测试

本文以丹江口水库为例进行系统实例测试, 以从星载RA数据提取的遥测水位、从星载MSI和SAR数据提取的水域面积, 作为多源遥感水库蓄水量估算模型的输入参数, 将输出的蓄水量与“全国水雨情信息网”上公布的标准蓄水量进行对比, 结果如图18所示, 日期范围是2019年1月1日至2019年12月31日. 由图18可知, 本系统所使用的多源遥感水库蓄水量估算模型得到的蓄水量和标准蓄水量的波动趋势具有较高吻合度. 在7~10月的汛期, 丹江口水库地区降雨较多, 蓄水量增长迅速, 在11月17日左右达到最高峰, 然后蓄水量开始慢慢减少, 逐渐进入枯水期, 在5月初达到最低值. 整体上, 随着旱季和汛季的转换, 丹江口水库的蓄水量呈现规律性变化.

采用拟合优度(R-square, R²)和均方根误差(Root Mean Squared Error, RMSE)对多源遥感水库蓄水量估算模型得到的蓄水量和标准蓄水量进行定量分析. 分析结果显示R²=0.99, 说明该模型所得蓄水量与标准蓄水量具有很高的吻合度, 能够真实地反映水库蓄水量的变化情况; 此外RMSE=1.055, 表明该模型估算的水库蓄水量标准差在1.055亿m³以内, 相对误差控制在1%以内.

5 结论与展望

本文开发了一款基于多源遥感技术的水库蓄水量估算系统, 在数据同化法的基础上, 通过建立多源遥感水库蓄水量估算模型, 实现水库蓄水量估算功能. 以丹江口水库为例进行实验, 结果表明本文建立的估算模型所估算的蓄水量与官方公布的标准蓄水量具有很高的吻合度, 能够真实地反映水库蓄水量的变化情况. 该系统能够以多源遥感方式获取水位、水域面积、蓄水量等水库参数, 可以很好地弥补目前单一遥感水库蓄水量估算系统的不足, 有望应用于防汛抗旱、水库调度等领域, 为管理部门科学化、高效化管理水库提供技术支撑.

本文所提出的水库蓄水量估算模型在建模阶段仅使用了数据同化法, 为了使蓄水量模型的估算结果更加精确, 在未来的研究工作中, 还可以尝试使用水文公式法、积分法、均值法等其他估算方法, 通过对比分析与实验验证, 构建和选用更合适的蓄水量估算模型, 以进一步提高蓄水量估算精度. 同时, 在系统的功能实现方面, 还可以添加出/入库径流量监测、水库水质监测等功能, 为多角度、全方位、立体化的水库监测研究提供更多的技术支撑.