基于多源数据融合的水库DEM构建方法研究
杨瑾 1崔海涛 1肖洲 1杨建英 2庄妍 3
1. 北京市水利规划设计研究院2. 北京林业大学3. 北京市密云水库管理处
摘要：针对水库水下地形获取难度大、成本高的问题，提出了一种利用长时间序列影像和水位数据构建水库数字高程模型的方法，以密云水库为研究区，基于Google earth engine云平台收集了2014—2023年370景Sentinel-1和297景Sentinel-2影像，结合SNIC分割和OSTU阈值算法进行面向对象的水体范围提取。根据影像采集当日的水位，采用K-Means方法将影像集聚类为45组，对每个分组内的水体提取结果进行合并，得到45个水位下的库容曲线。通过库容曲线提取了467 601个高程点，采用创建TIN的空间插值方法构建5 m分辨率的DEM,基于DEM计算的库容与密云水库管理处观测库容基本吻合（R2=0.999 8）,在此基础上，探究了密云水库“水位-面积-库容”关系模型和水深分布特征。结果表明，密云水库呈现“西南深、东北浅”“白河入口深、潮河入口浅”的特征，水位与面积呈强线性相关关系（R2=0.994 7）、水位与库容呈明显的二次多项式关系（R2=1）,根据建立的关系式计算面积相对误差约2%,计算库容的相对误差小于1%,可为水库的监测和管理提供有效数据。
引用：[1] 杨瑾,崔海涛,肖洲,等. 基于多源数据融合的水库DEM构建方法研究[J]. 测绘科学, 2024, 49 (07): 184-194.
DOI:10.16251/j.cnki.1009-2307.2024.07.019.
0 引言
新阶段水利行业亟需进行数字化和智慧化升级，构建水库数字孪生流域，从而实现水库安全管理、水资源调配等业务智能化 [1] 。其中，表征水库“三维空间”的DEM是“数字孪生”的基础底板数据 [2] ,主要采用声呐测深技术获取高精度的水下地形，并实现了单波束到多波束、有人船到无人船技术转变，使得水下地形数据获取实现降本增效 [3-4] ,但对于大范围的水库，通过无人船观测水下地形，仍需投入较多人力、物力。
随着遥感技术的不断发展以及时间的推移，完成了水库不同水位监测影像的积累，为水库三维模型的构建提供了新的思路。文献 [5]利用6 424个内部水深点位和Landsat影像提取不同水位的湖岸线，通过创建TIN的方法构建摩苏尔(Mosul)湖泊DEM。文献 [6]从时间序列Sentinel-2影像中提取124.35～144.62 m水位范围内的库容曲线，并对比ANUDEM和TIN两种方法基于不同水位库容曲线构建水库DEM的精度，结果表明TIN方法优于ANUDEM。文献 [7]利用时序Sentinel-2和ICESat-2数据，采用Kriging插值构建台特玛湖盆DEM,高程中误差约为0.103 m, DEM提取的等高线与实际水面线高度吻合。利用水面等高特性，通过长时间序列的历史影像提取湖泊或水库的等高线，采用空间内插方法构建水库DEM,进而将问题转为提取不同“水位”的“水体范围”。
目前，通过遥感手段提取水位主要采用微波雷达测高和激光测高数据，例如，激光测高卫星ICESat-1/2 [8-9] 、星载雷达干涉测高卫星SWOT [10] ,但对于较大型水库，一般都设立水文观测站点，可实现水位实时动态变化监测。水体范围提取主要采用被动光学遥感和主动合成孔径雷达两种遥感数据，是应用较早和发展较成熟的技术，其中，遥感影像空间分辨率是影响提取范围精度的关键参数，适合采用较高分辨率的影像 [11] 。被动光学遥感技术发展较早，具有数据源丰富和处理技术成熟的优势，但易受云雨天气的影响。主动合成孔径雷达具有较强穿透性，可以有效避免云雨天气下数据失效的问题，同时，能够获取较高空间分辨率影像，对于水体边界提取具有优势，但其技术发展较晚，存在历史数据不足，处理技术尚未成熟，影像存在大量噪点等问题 [12] 。利用上述两种传感器数据进行水体提取时，均会受地形的影响 [13] ,光学影像中山体阴影与水体产生相似的光谱特征，在水体提取中会存在将山体阴影错提为水体的现象。利用合成孔径雷达数据提取水体时，依据水体在微波波段的后向散射强度，其本身会包含较多的地形信息，地形平坦的水体与其他平坦的水田、道路等地物会产生相似的后向散射强度，会出现将平坦的地物错提为水体的现象。
单一传感器技术获取的影像无法冲破其本身的技术瓶颈，通过多模态数据融合技术，可以突破常规手段，助力水下地形的提取 [14] 。因此，本文参考上述分析，以密云水库为研究对象，基于GEE云平台，收集2014—2023年长时间序列Sentinel-1和Sentinel-2影像数据，结合云量掩膜、SNIC图像分割以及OSTU阈值分割算法，提取不同日期对应的水体范围。根据水文站实测水位数据，采用K-means聚类法，依据水位相近原则，将不同日期水体范围提取结果进行分组，并对每个分组内的结果进行融合，以提高水体提取结果的鲁棒性，从而，提取出不同水位的库容曲线，最后，采用空间内插创建TIN的方法构建水库DEM [15] 。
密云水库(40°26′～40°35′N,116°48′～117°06′E)位于华北平原燕山丘陵之中，是华北最大的调蓄水库，有“燕山明珠”之称。水库于1960年9月建成，面积约180 km2,库容4.375×109 m3,分为东、西两个库区。密云水库有潮河和白河两大支流，其中，白河起源于河北省沽源县，经赤城县、延庆区、怀柔区入库，潮河起源于河北省丰宁县，经滦平县，自古北口入库。潮白河流域具有明显的暖温带季风性半湿润半干旱气候特征 [16-17] ,干湿变化明显，夏季降雨多，易发生洪灾，冬季则雨雪少，易发生旱灾。
1.2 数据源及预处理
本文通过GEE云平台收集密云水库2014—2023年的Sentinel-2 MSI和Sentinel-1 SAR影像 [13,18] 。对影像产品数据集进行筛选，筛选原则：①避免使用冰冻和消融期的影像；②避免使用山体阴影较多的影像；③确保研究区内Sentinel-2影像不受云雨影响；④确保研究区内影像的完整性。根据影像筛选原则，制定具体筛选方法，通过对比不同时期的影像，发现5月水库水体已完全消融，10月之后影像山体阴影较多，确定利用GEE中ee.filterDate和ee.filterBounds方法筛选研究区范围内5—10月的影像。结合“QA60”波段对每景Sentinel-2影像逐像元进行云和卷云的掩膜处理 [19] 。最后，筛选出研究区范围内有效像元大于98%影像，共收集了370景Sentinel-1影像和297景Sentinel-2影像(图1)。
2 研究方法
本文利用密云水库不同水位的Sentinel-1和Sentinel-2影像构建DEM,实现水库三维空间特征分析，主要包括以下内容。
1)水体范围提取。采用SNIC影像分割和OSTU阈值法进行面向对象的水体提取。
2)不同水位库容曲线提取。在每景影像水体范围提取基础上，结合影像采集当日水文站实测水位，采用K-means聚类法对影像进行分组，并将同一分组内影像的水体提取结果进行合并，提取出不同水位库容曲线。
3)DEM构建。根据不同水位库容曲线，利用创建TIN空间内插方法生成DEM,并从水体面积和库容两方面进行精度评价。
4)密云水库三维特征分析。基于生成的DEM,分析密云水库水位、面积、库容的关系，以及水深的分布特征。
2.1 水体提取算法
为实现长时间序列数据的高效收集和处理，本文基于GEE云平台进行影像预处理和水体提取 [20] ,具体流程详见图3。因Sentinel-1/2影像存在噪点和水陆交界处存在混合像元，提取的水体存在空洞和碎屑多边形 [21] 。针对此问题，本文图像分割采用简单非迭代聚类SNIC方法，该方法由文献 [22]在简单线性迭代聚类算法(simple linear iterative clustering, SLIC)基础上进行优化，采用一个优先级队列替换SLIC算法的k-means迭代聚类过程，从而减少内存用量，更快速地生成超像素，并集成于GEE平台，可快速实现影像面向对象分析。其原理是根据影像中不同波段的权重、紧致度和平滑度进行分割，分割为相似的均质斑块对象，关键参数是seeds和compactness, seeds作为初始种子点间隔聚类，值越大，代表分割对象聚类中心越远，分割的对象也越大，compactness决定分割对象的紧致程度，值越大，分割对象越规则。针对本研究，经试验，设置seeds参数为10、compactness参数为0.5时，分割效果最佳。在此基础上，分别求取分割对象内Sentinel-1影像VV极化和Sentinel-2近红外波段(B8)的平均值，采用OSTU法(最大类间方差法)确定分割阈值 [23] ,进而提取出每景影像的水体范围。其中，OSTU法是基于最小二乘法原理的一种自适应图像分割阈值确定算法，设定某一阈值把影像像元分配为水体和非水体，计算类间方差，当方差达到最大值时，此阈值则为最优水体提取分割阈值。
结合水文站观测水位，将影像获取当日水位赋值给每景影像，采用K-Means聚类法根据水位属性进行影像分组 [24] ,然后，对每个分组内多景影像的水体提取结果进行合并，合并算法是：以像元为单元进行合成，若对应像元内提取为水体的概率≥60%,则认定为水体，反之为非水体(图3),计算见式(1)。
2.2 DEM精度评价方法
本文参考文献 [6,25]的研究，从水体面积和库容两个方面对构建的DEM进行精度评价，选取平均误差(ME)、平均绝对误差(MAE)和中误差(RMSE)分析模拟值与真实值误差的大小，选取绝对系数(R2)描述模拟值与真实值之间的拟合关系，相关计算见式(2)～式(5)。
通过目视判读，发现Sentinel-1和Sentinel-2影像提取水体结果中存在颗粒状的“杂质像元”,为改善此现象，本文采用面向对象分析方法进行水库水面提取，根据像元之间光谱属性的相似性进行空间聚类，将像元聚类为相似的斑块。通过对比Sentinel-1和Sentinel-2面向像元和对象的水体结果(图4),发现面向对象分析方法提取的水体更完整，可显著减少面向像元提取结果中存在“细碎水体”的现象，有利于提高后续水体范围的提取精度。
3.1.2 影像分组情况分析
通过影像筛选，共收集了密云水库370景Sentinel-1和297景Sentinel-2影像。根据影像获取当日水位，采用K-Means聚类法，将影像划分为45组，计算每组影像对应平均水位。由图5可知，Sentinel-1在45个水位分组内均有影像，每个分组内影像个数在1～29景不等，在水位149.13 m左右收集Sentinel-1影像数量最多；Sentinel-2影像在39个水位分组内有影像分布，每个分组内影像个数在1～53景不等，在水位151.52 m左右收集Sentinel-2影像最多。相比之下，Sentinel-1影像采集的数量较多，且覆盖的水位范围更广，更适用于提取不同水位库容曲线。
3.1.3 Sentinel-1和Sentinel-2水体提取结果对比
将Sentinel-1和Sentinel-2(图6)数据提取的水体结果进行对比，发现通过Sentinel-2提取时，会将水库周边的山体阴影误提为水体，而Sentinel-1提取水体结果可有效避免此现象，提取水体边界与实际更吻合。
结合高分辨率影像对Sentinel-1和Sentinel-2的提取结果进行分析，如图7所示，发现Sentinel-1提取的水体范围偏小，造成此现象的原因是Sentinel-1具有一定的穿透性，对于浅水面，Sentinel-1会穿透水面，得到非水体表面的反射特征。Sentinel-2影像穿透性差，反映的是表层水体的反射特征，能够识别出浅水面的水体。在地形比较平坦的地方，Sentinel-2提取的水体边界较准确，Sentinel-1提取的水体结果在一些细小的分支处无法提取到边，此时，Sentinel-2的水体的提取结果更符合实际水体的边界。然而，由图7发现，在水库边界地形复杂或者植被茂密的地区会存在将山体阴影或植被茂密覆盖区域错提取为水体，此时，Sentinel-1提取结果更符合实际。
针对每个分组内的影像数据集，以像元为单元对Sentinel-1与Sentinel-2提取结果进行融合，分别计算每个像元通过Sentinel-1和Sentinel-2影像识别为水体的概率，然后求得两种数据源提取为水体概率的平均值，将概率≥60%的识别为水体。从图7(绿色)中可以看出，Sentinel-1与Sentinel-2融合实现两种数据优势互补，有效避免Sentinel-1无法识别出浅层水面、Sentinel-2将山体阴影误提为水体的现象，提高了水体范围提取精度。因此，本文采用Sentinel-1与Sentinel-2融合的算法提取45个水位下的库容曲线，如图8所示。
3.2 密云水库DEM建立及精度分析
首先，从不同水位的库容曲线中提取高程点，利用Arcmap中要素折点转点工具将库容曲线折点转换为点要素，结合高分辨率影像进行目视检查，剔除异常点，最终得到467 601个高程点，每个水位下高程点个数在6 866～14 469不等(表1),且随水位增加个数呈现增多的趋势。然后，通过空间内插创建TIN的方法得到密云水库的不规则三角网(TIN),最后，采用TIN转栅格工具生成5 m分辨率的DEM(图9),其中，生成的DEM与密云水库管理处水文站观测水位采用的高程基准一致。
为了进一步验证DEM的精度，基于本文构建的密云水库DEM,计算不同水位下的水体面积和库容，并与Sentinel-1与Sentinel-2影像融合提取的水体面积和密云水库管理处观测库容(参考库容)进行对比分析。由于本文构建的DEM缺少132.98 m水位下的数据，基于DEM计算库容是以132.98 m为基准的库容，而密云水库管理处观测数据反映的以真实库底为基准的库容。首先，将以132.98 m为基准的库容转换为以133 m为基准的库容，然后，根据水文站观测的库底到133 m的库容量，将以133 m为基准的库容转换为以真实库底为基准的库容，具体转换关系见式(6)。
通过影像提取面积对DEM建模精度进行验证(图10),结果表明，DEM计算面积与影像提取水体范围面积呈明显线性关系(R2=0.989 9),验证点分布在y=x两侧，中误差(RMSE)约3.01 km2,说明本文DEM建模方法的拟合精度较好，验证了空间内插创建TIN方法的可行性。
通过影像提取面积对DEM建模精度进行验证(图10),结果表明，DEM计算面积与影像提取水体范围面积呈明显线性关系(R2=0.989 9),验证点分布在y=x两侧，中误差(RMSE)约3.01 km2,说明本文DEM建模方法的拟合精度较好，验证了空间内插创建TIN方法的可行性。
基于密云水库管理处观测库容(参考库容)进行DEM精度验证(图11),结果表明，两者变化趋势基本一致，呈现强线性相关关系(R2=0.999 8)。DEM计算的库容均小于密云水库管理处观测数据，平均误差(ME)与平均绝对误差(MAE)大小均约9.8×107 m3(表2),中误差(RMSE)约1.14×108 m3,平均绝对误差(MAE)与库容成明显的正相关性，库容增大，MAE也随之增加。图12显示了MAE与DEM计算库容之间的关系，结果表明，两者呈明显的二次多项式关系(R2=0.992 5)。
根据DEM计算库容的误差关系模型建立修正公式(式7),经修正之后库容与密云水库管理处观测库容(参考库容)高度吻合，绝对系数R2为1,两者之间基本无偏差(图13)。
3.3 密云水库水位与水体面积、库容关系
利用DEM计算库容、水体面积并探究两者关系，从图14明显发现两者呈二次多项式关系(R2=0.998 5),根据此关系，结合计算的修正库容对水体面积进行修正，见式(8)。
密云水库水位与水体面积存在明显的线性关系，R2为0.994 7,截距为-567.1,斜率约为4.753,见图15。水位每抬升1 m, 水体面积约增加4.753 km2,见式(9)。
根据密云水库管理处观测数据得知密云水库水位为146.31、151.83 m对应的水面面积130.95、157.83 km2,水库由146.31 m上升到151.83 m时面积增加26.88 km2。通过式(9)计算的两个水位水面面积分别是128.31、154.55 km2,面积增量26.24 km2,均比实际少2.64、3.28、0.64 km2,计算面积及其变化量的平均相对误差约2%。
水位与库容呈二次多项式关系，见图16。
3.4 密云水库水深空间特征分析
水深作为表征水库三维特征的重要指标 [26] ,是水库管理的重要参考数据，本文利用构建的DEM对密云水库的水深分布情况进行分析，以155.16 m为水平面计算水深(图17),通过重分类方法对水深进行分级，分为0～5 m(12%)、5～10 m(14%)、10～15 m(16%)、15～20 m(12%)以及>20 m(46%)五级，发现密云水库近一半区域的水深大于20 m, 集中分布在水库的西南区域，而东北区域较浅。白河支流入口区域的水深较深，潮河支流入口区域水深较浅。
4 结束语
本文基于GEE云平台和多源数据，以密云水库为研究区域，探究利用长时间序列遥感影像和水位数据构建水库“三维”模型(DEM)的可行性。水体范围作为DEM构建的关键数据，该文对其提取方法进行探究，最终提出一套不同水位库容曲线提取流程。利用影像提取的水体范围和密云水库管理处的观测数据，从面积和库容两个方面对DEM精度进行分析。此外，通过DEM探究密云水库水位、水体面积、库容的三者的关系以及水深分布特征，得到4点结论。
1)在水体范围提取时，面向对象方法可避免影像中“噪声”干扰，提高库容曲线的提取精度，在地形平坦区域，光学遥感(Sentinel-2)水体提取结果优于合成孔径雷达(Sentinel-1),而在地形陡峭或植被茂密的区域，Sentinel-1水体提取结果优于Sentinel-2。通过融合两者提取的水体结果可以提高水体范围提取精度，有效避免Sentinel-2将山体阴影误提为水体、Sentinel-1无法识别出浅层水面的现象。
2)利用密云水库管理处观测数据验证了本文构建的DEM能够准确地表达水库的三维特征，但DEM计算的库容相对较小，且平均绝对误差与库容呈二次多项式的正相关关系，随着水位抬升，库容变大，对应DEM计算的库容误差也随之增大，分析其原因发现，通过影像提取水体范围偏小，导致提取不同水位库容曲线向内收缩，通过其构建的DEM计算的水体面积、库容会普遍较小。
3)密云水库水位、水体面积、库容三者之间存在明显的相关关系，水体面积与库容(R2=0.998 5)、水位与库容(R2=1)呈强二次多项式关系，水位与水体面积(R2=0.994 7)呈强线性相关关系，水位每抬升1 m, 水体面积增加4.753 km2。通过建立的水位与水体面积、水位与库容关系式，计算某一水位下的水体面积相对误差约2%,计算库容相对误差<1%。
4)根据DEM分析密云水库水深分布特征，发现呈现“西南深、东北浅”的分布特征，白河支流入口处较深，潮河入口处较浅。
本文提出一种利用多时相“二维”影像和水位构建水库“三维”数据的方法，采用此方法提取了密云水库132.98～155.16 m水位范围内的DEM,并验证该方法的可行性，一定程度上提高水下地形获取效率，同时也为其他水库DEM构建提供思路。本文提出的构建水库DEM存在两点局限性：①仅能构建可获取的库容曲线水位范围内的DEM;②DEM构建精度取决于水体范围的提取精度，受影像数据分辨率(10 m)限制，本文中无法识别面积小于100 m2水体，影响DEM构建精度。在今后研究中，应考虑使用更高分辨率的影像数据，提高水体范围的提取精度，例如，国产高分系列、资源系列卫星 [17] 。同时，针对某一水位范围内因无法获取库容曲线，导致构建DEM缺失的问题，可结合无人船水下地形测量技术进行补测，从而提取水库完整的“三维”数据。
基于多源数据融合的雅丹体三维建模
西宁市产业园区时空分布研究，基于多源数据
一个城市路网生成模型，基于多源大数据
基于多源遥感数据的森林火灾应急监测
