本文介绍了一种用于未勘测水库合成 bathymetry（水下地形）的 hydrologically conditioned cubic spline interpolation（水文约束的三次样条插值）方法。该方法由 Camilo J. Bastidas Pacheco 和 Denny Rivas 提出，发表于 Journal of Hydrologic Engineering（2023 年 12 月，第 28 卷第 12 期，文章编号 04023038，DOI: 10.1061/JHyeff.HEENG-5885）。作者单位分别为 Utah State University 土木与环境工程系及 Utah Water Research Laboratory，以及独立水资源与灾害风险建模顾问。本研究针对全球大量水库缺乏 bathymetric data（水下地形数据）的现实问题，发展了一种仅依赖水库周围数字地形模型（DTM）即可生成合成水下地形的低成本和高效方法。

全球范围内，水坝在水文、生态、社会和经济系统中扮演双重角色：一方面提供供水、发电、防洪等功能，另一方面又会阻断生态廊道、改变栖息地，甚至引发溃坝等风险。然而，大量水库，尤其是发展中国家、偏远地区或冲突地区的水库，缺乏基本的水下地形数据和库容信息。已有的全球数据库如 GOODD、Global DamWatch、GeoDAR 多只记录了水库位置、坝高或总库容等属性，极少包含详细 bathymetry。传统的水下地形获取方式包括船载声呐测量、SAR 影像、LiDAR 和线性空间插值等，但船载测量成本高且受地理与安全条件限制，SAR 对水位变化不充分的水体存在盲区，LiDAR 对水深和浊度敏感，线性插值则易导致最大深度的高估。针对这些不足，作者开发了基于三次样条空间插值的水下地形估算方法，其核心是使用 SAGA-GIS 中的 Close Gaps with Splines 算法，并加入水文约束以控制深度和纵向连续性。方法所需输入仅为水库周边 DTM，适用于数据严重匮乏地区。

研究区包括分布在美国西部七个州的 12 座水库，面积从 5.78 km²（Gibson）到 64.12 km²（Theodore Roosevelt Lake），最大深度从 23.21 m（Beulah）到 100.51 m（Anderson Ranch）。地质背景包括沉积岩、火成岩以及二者的混合。研究所用 DTM 来自 USGS 3D Elevation Program（3DEP）的 ¹⁄₃ 弧秒（约 10 m）分辨率产品，验证数据则来自美国垦务局（USBR）的船载 bathymetric surveys。作者将所有数据统一到 UTM 投影、NAD83 水平基准和 NAVD88 垂直基准下，DTM 通过双线性重采样转换为 10 m 分辨率。

方法流程包括四个主要步骤。第一步是创建 gap（空值区）以去除 DTM 中的水面高程。静态水体在 DTM 中表现为一组高程一致的平坦像元，这些像元不包含水下地形信息，需要以 mask polygon（掩膜多边形）覆盖并转换为空值。掩膜需满足四条准则：完全覆盖水体平坦像元；延伸覆盖坝体、泄水建筑物、桥梁等人工结构；向下游延伸至未受扰动的河道，以便插值结果与下游河道拓扑衔接；在水体边界外再外扩 3–4 个像元，以消除岸线过渡像元导致的坡度失真。第二步是无水文约束的初步插值。使用 SAGA-GIS 的 Close Gaps with Splines 算法对含空值的 DTM 进行填充。该算法采用三次样条，在二维空间中使用 Moore 八邻域方案（包含正交和对角方向），通过多次迭代求解 16 个多项式系数，保证插值面在已知点处精确且一阶导数连续。唯一标定参数是驱动插值的节点数 n，节点数影响插值面的刚性与柔性，从而影响生成的库容。第三步是构建库区栅格河网。首先对初步插值结果进行水文校正，采用 SAGA-GIS 中的 Sink Removal 算法进行深切而非填洼处理，以避免恢复水体表面高程。之后基于 D8 流向算法计算流向和汇流累积量，提取最低高程剖面（thalweg）作为主河道或支流。主河道沿水流方向从入库口到出库口按线性插值确定高程；支流高程以与主河道交汇处为下游端进行线性插值，采样间距为 30 m。将河道栅格点与空值 DTM 合并，形成 stream-DTM。第四步是水文约束的最终插值。将 stream-DTM 再次输入 Close Gaps with Splines 算法，由于河道高程已知，插值结果在保持样条曲率的同时受河道深度控制，生成水文一致的合成水下地形。最后应用 5×5 窗口的样条重采样滤波，先上采样后下采样，以平滑表面并保留主要形态特征。

评估采用三种方式：总库容误差、深度–库容曲线比较以及逐像元深度统计。性能指标包括 bias（偏差）、Kling–Gupta efficiency（KGE）和 RMSE。KGE 综合考虑偏差、相关性和方差，相较于 Nash–Sutcliffe efficiency 更适合评价模型整体性能。作者采用 KGE ≥ 0.5 和 bias 绝对值 ≤ 25% 作为满意阈值。敏感性分析通过改变节点数 n 从 2 到 20，分析不同地质背景下性能指标的变化。结果显示，优化后 12 座水库的合成 bathymetry 全部达到满意性能：bias 介于 −8.8% 到 +1.9%，KGE 介于 0.53 到 0.90 之间，RMSE 介于 3.53 m 到 12.73 m。在逐像元比较中，性能最佳的水库包括 Gibson（KGE=0.90）和 Prineville（KGE=0.89），最低的为 Elephant Butte 和 Warm Springs（KGE 均为 0.53）。从总库容误差角度看，除 Millerton Lake 和 Warm Springs 等少数水库外，大部分误差在 ±15% 以内。深度–库容曲线比较显示 KGE 值从 0.71 到 0.94，R² 接近 1，表明合成曲线与实测曲线在库容随水深变化趋势上高度一致。值得注意的是，深度–库容曲线的 R² 普遍较高，但作者指出该指标在单调曲线比较中阈值较低，KGE 更具参考意义。

敏感性分析揭示了一个重要规律：KGE 和 bias 曲线根据下伏地质材料分为两类，而与水库尺寸无显著关系。对于火成岩或冰川沉积物（如 Anderson Ranch、Millerton Lake、Prineville、Tiber、Gibson），KGE 随节点数增加而单调上升，在 14 个节点时达到最优；而沉积岩或沉积岩与火成岩混合区（如 Elephant Butte、El Vado、Lake Pleasant、Warm Springs 等）中，KGE 在约 4 个节点时达到峰值，随后随节点数增加而下降。这种现象可解释为：火成岩地形多具不规则起伏和更复杂的局部地貌，需要更多节点以增强插值面灵活性；沉积岩区则多形成宽阔平坦的谷底，过多节点易导致过拟合和表面畸变。作者据此提出全局性参数选择准则：地质为火成岩或冰川沉积物时选择 n=14，沉积岩或混合岩性时选择 n=4。该准则仅需地质图即可确定最佳插值参数，为无实测数据地区的应用提供了可操作依据。

研究还讨论了影响插值精度的其他因素，包括泥沙淤积速率、船测数据距离水库建设的时间差、土地利用变化以及超出岩性之外的局地地貌过程。由于方法假设水库底部形态与周围地形保持连续，未考虑建库后的淤积影响，因此合成 bathymetry 更接近建库前的地形状态。作者认为，通过调整节点数或移动 thalweg 高程，可以生成地形上界和下界，从而间接反映淤积带来的不确定性，但这些方向仍需进一步研究。研究也未探讨不同 DTM 分辨率对插值结果的影响，建议在不同数据源情况下将 DTM 重采样至 10 m 以保持一致性。

该方法的应用价值体现在：可服务于早期水文、水环境、水库溃坝与拆除等风险评估场景，在预可行性或概念设计阶段能够提供合理的库容和地形信息；也可与 SAR、LiDAR 等遥感数据互补，当遥感技术仅能获取浅水区深度时，本方法可提供深水区估计；此外，还可先以合成 bathymetry 识别高复杂度或高精度需求区域，再配置有限的船载测量资源，形成混合测量模式，从而减少测量活动的时间和成本。对位于偏远、危险或有社会冲突地区的水库，该方法具有突出的实用意义。

本研究的亮点在于：第一，提出了一种以三次样条插值为核心的水下地形合成方法，相较于常用的线性插值具有更好的曲率适应性和深度控制能力；第二，在插值流程中加入水文约束，通过自提取的 thalweg 控制纵向坡度与连续性，避免了人工洼地和峰值，使结果更符合自然地形；第三，通过敏感性分析揭示了最优节点数与地质岩性之间的关系，建立了基于地质背景的参数选择规则，将大尺度地质过程与小型水库尺度插值行为联系起来；第四，方法仅需公开 DTM 数据，计算过程无需专门设备，时间与经费成本低，适合数据稀缺地区的大规模应用。总体而言，该研究为解决全球范围水库 bathymetry 数据缺失问题提供了一种可操作、可推广的技术路径，其结论对于水文建模、灾害风险评估和水资源规划具有实际的科学与工程价值。