一种用于从卫星图像中识别加利福尼亚州圣华金河谷含水层 recharge 位置的分类学习流程

《Journal of Hydrology: Regional Studies》：A machine learning pipeline for identifying infiltration managed aquifer recharge locations from satellite imagery in the San Joaquin Valley, California

　　双宇|Ate Visser|Indrasis Chakraborty|Gemma J. Anderson|杨柳|Jory Chapin Lerback
美国劳伦斯利弗莫尔国家实验室

**摘要**
**研究区域**
本研究聚焦于美国加利福尼亚州中央谷地的一个农业区域，该地区广泛实施了人工含水层补给（Managed Aquifer Recharge, MAR）措施，以缓解日益增长的用水需求和气候变率带来的地下水枯竭问题。

**研究重点**
本研究开发了一种基于深度学习和机器学习的框架，利用卫星图像和环境数据来识别人工含水层补给的位置。该框架整合了来自Sentinel-2卫星图像的地表水检测、水体的地理空间划定、水体动态的时空追踪，以及基于气象、环境和地形变量的监督分类方法。该框架应用于弗雷斯诺市西南方向2379平方公里的研究区域，共检测到765个水体，其中139个是根据公开数据和专家知识确认的人工含水层补给站点。分类模型的准确率达到0.94，F1分数为0.85。

**对该区域的新水文见解**
特征重要性分析表明，农田、归一化植被指数（NDVI）和蒸发量是影响人工含水层补给效果的主要因素。值得注意的是，该框架指出，在人工含水层补给系统中，人工水资源管理可能会破坏甚至逆转地表水面积与降雨量之间的预期正相关关系。这些发现为现有人工含水层补给设施的特性提供了可物理解释的见解，并展示了所提出框架在日益复杂的气候和水文不确定性条件下，作为可复制、可解释且具有潜在通用性的工具，用于数据驱动的人工含水层补给位置识别和清单编制的潜力。

**1. 引言**
由于气候变化导致的长期干旱、城市化和人口增长，全球范围内对水资源的的需求不断增加，给地下水资源带来了巨大压力（Bouwer, 2002; Vanderzalm et al., 2010）。人工含水层补给（MAR）已成为一种越来越常见的方法，用于增加地下水资源，以补充自然或非人为的自然含水层补给（Dillon et al., 2009; Maliva and Missimer, 2012; O'Leary et al., 2012）。MAR通过增加地下水资源并在地表水供应不足时起到缓冲作用，提高了干旱抵御能力（Yuan et al., 2023）。它还有助于防止过度抽取地下水的地区发生地面沉降，并通过利用地表水或雨水补充含水层来改善水质。全球范围内已实施多种MAR技术，主要包括渗透型MAR（infiltration-MAR）和注入型MAR（injection-MAR）。渗透型MAR涉及对土地表面的改造以增强补给效果，例如铺水方法、河道改造以及雨水和径流的收集；而注入型MAR则包括通过井、涵洞和钻孔进行补给（Huber and Scheibler, 2013; Hannappel et al., 2014; Ringleb et al., 2016）。基于渗透原理的MAR系统依赖于适宜的地下水文条件，这些条件会影响补给能力、地下水储存量以及补给操作的总体效果（Sinshaw et al., 2026）。本研究的具体目标是识别当前正在运行的人工含水层补给地点。

**2. 数据**
**2.1. 研究区域**
加州（美国）的严重干旱和地下水过度抽取导致了地下水位下降和水资源短缺。该州已启动多项MAR项目，以增加地下含水层的储水量（CNRA, 2020; Sheik et al., 2024）。在加利福尼亚州，中央谷地受不可持续地下水抽取的影响最为严重。中央谷地是一条长700公里、宽80公里的河谷，西邻海岸山脉，北接喀斯喀特山脉，东靠内华达山脉。该地区具有地中海气候特征，年平均降水量从南部的约150毫米到北部的800毫米不等（PRISM Climate Group, http://prism.oregonstate.edu）。中央谷地90%的区域的年潜在蒸发量超过年降水量（Azhar et al., 2025），范围在706毫米至957毫米之间，其中降水量最少的地区蒸发量最高（Wolock, 2003）。中央谷地地下水系统分为萨克拉门托含水层盆地（旧金山湾三角洲以北）和圣华金含水层盆地（三角洲以南）（CNRA, 2020）。地表水贡献了圣华金含水层盆地新补给水量的一半（Visser et al., 2018）。本研究关注的区域位于弗雷斯诺市西南方向，地理范围为北纬36.5°至36.9°，西经120.1°至119.5°（图1）。该区域包括由弗雷斯诺灌溉区（FID）、弗雷斯诺都市防洪区（FMFCD）和联合灌溉区（CID）实施的人工含水层补给操作。研究区域内还有三条主要水道：北部的圣华金河、东南部的金斯河以及西南部的弗雷斯诺沼泽。其他地表水特征包括大型污水处理设施、灌溉渠道和防洪 Basin。该地区的地表层由未固结的冲积沉积物组成，这种地形有利于基于渗透原理的人工含水层补给（Ferriz, 2001）。

**2.2. 遥感数据**
我们收集了2022年和2023年4月1日至10月3日期间Sentinel-2卫星（Drusch et al., 2012）的遥感图像。数据采集主要集中在夏季，因为人工含水层补给设施在春季径流期间可能会被淹没，而夏季的云层覆盖通常较低。晚春和夏季时，周围的自然水体和农业区往往会经历季节性干燥，而人工含水层补给区由于持续补水而保持淹没状态，从而形成明显的地表水动态对比，有助于提高补给位置的识别效率。2022年是一个干旱年份，弗雷斯诺的降水量仅为170毫米，而2023年的降水量为373毫米（PRISM Group, 2025）。数据采集间隔为5至10天，共计55个观测日期。每幅图像包含六个光谱波段：02波段（蓝色）、03波段（绿色）、04波段（红色）、08波段（近红外，NIR）、11波段（短波红外SWIR1）和12波段（短波红外SWIR2）。这些遥感数据被用于WatNet模型（第3.1节）以辅助水体检测。此外，同一数据集还被用于计算每个网格单元的归一化植被指数（NDVI）。NDVI的计算公式为：NDVI = (NIR - RED) / (NIR + RED)，其中NIR表示第8波段，RED表示第4波段。本研究使用了Sentinel-2 Level-2A地表反射率产品。通过Sentinel-2数据提供的质量评估波段移除了云层和云影像素，只保留了无云的观测数据进行分析。所有光谱波段在使用之前都被重新采样到统一的空间分辨率。处理后的影像被用于地表水体检测和NDVI计算。

2.3. 地形数据
海拔数据来自美国地质调查局（USGS）的3D海拔计划（https://www.usgs.gov/3d-elevation-program）。空间分辨率为1/3弧秒，相当于研究区域大约10米。

2.4. 气象数据
2022-2023年的每小时网格化气象数据来自ERA5-Land再分析数据集（https://cds.climate.copernicus.eu/datasets/reanalysis-era5-land?tab=download）。空间分辨率为0.25° × 0.25°。选取的变量包括：降水量、2米高度的空气温度、开阔水面的蒸发量（不包括海洋）、植被蒸腾作用引起的蒸发量、潜在蒸发量和总蒸发量。

2.5. 环境数据
我们整合了Wolock等人（2003年）、Reitz等人（2017年）、Falcone（2015年）和Faunt（2009年）的研究中的每日环境变量，包括有关农田、径流系数、蒸散量、低强度城市区域、饱和导率、地表水体、半开发土地、牧场、放牧区以及开发/城市化土地的信息。这些变量的空间分辨率根据数据来源的不同，范围在大约0.7平方公里到3.6平方公里之间。

由于这些环境和气象数据集的空间分辨率不同，所有变量都基于包含多边形质心的网格单元与每个水体多边形进行了空间匹配。因此，每个水体都分配了来自每个数据集的相应网格值。这种方法使得具有不同空间分辨率的多源数据集能够整合到一个统一的建模框架中。然而，这可能会引入一些不确定性，特别是对于像ERA5-Land这样的低分辨率数据集，因为多个水体可能共享同一个网格单元，从而具有相同的气象数据。在这种情况下，气象变量主要代表背景水文气候条件，并捕捉与人工湖（MAR）出现相关的更广泛的水文气候变化和环境条件。

3. 方法
我们开发了一个基于深度学习和机器学习的框架，通过整合遥感、气象、地形和环境数据来识别当前正在运行的人工湖位置。工作流程包括以下主要组成部分（图2）：
（1）使用遥感影像检测水体；
（2）使用Segment Geospatial模型（Wu, Osco, 2023）定义水体边界；
（3）为每个水体分配唯一的ID号码，并跟踪其随时间的变化；
（4）根据预定义的标准手动将每个水体标记为人工湖或非人工湖；
（5）训练分类模型以预测人工湖的位置。

### 3.1 使用遥感影像检测水体
首先使用预训练的深度学习模型（WatNet）处理遥感影像以检测地表水体。WatNet是一种专为遥感影像设计的高精度水体提取方法（Luo等人，2021年）。它包括三个主要组件：全局多注意力融合模块（Global Multi-Attention Fusion）、水体前向网络模块（Water Forward Network）和边缘聚焦注意力模块（Edge Focus Attention）。通过这些组件的整合，WatNet能够有效地捕捉细小的水体，并在复杂背景（如密集植被、地形变化或云层干扰）中增强水体的边界划分。在这项工作中，我们使用了Luo（2021年）开发的预训练WatNet模型（可从https://github.com/xinluo2018/WatNet获取）。输入影像覆盖2022-2023年的55天数据，每个场景包含2078行纬度和2500列经度，以及六个光谱波段。输出是一个二值分类图（以下简称W，其中W(i, j) = 1表示像素位置(i, j)存在地表水，W(i, j) = 0表示没有地表水）。

### 3.2 通过基于深度学习的分割来确定水体边界
第3.1节生成的二值W仅指示了地表水的存在与否，但没有提供关于各个水体边界或中心位置的信息。为了提取边界和位置信息，我们使用了segment-geospatial框架（Wu, Osco, 2023）。Segment-geospatial是一个开源的、预训练的Python框架，基于Meta AI的Segment Anything Model（SAM，Kirillov等人，2023年）。它提供了一种高效且可扩展的方法，用于从高分辨率遥感影像和地理空间栅格数据集中提取对象级别的特征。

在我们的工作流程中，对于每个观测日，二值水体分类图W首先被划分为16×16的网格瓦片，共得到256个不重叠的子图像（图3a）。这种分块策略对于分割精度至关重要：直接在整个图像上应用分割方法会由于对象数量多和空间复杂性高而降低精度。通过将W划分为较小的瓦片，我们可以减少每个瓦片中的水体数量，使模型能够更好地捕捉细小和破碎的水体。

### 3.3 为每个水体分配唯一标识符并跟踪其随时间的变化
水体的表面积和中心位置会随着时间的推移因降水或人工湖设施上地表水的应用而显著变化。准确识别和跟踪多个观测日期中的同一水体对于构建一致的输入特征和支持人工湖位置评估至关重要。我们选择2023年4月1日作为“参考”日期，因为该日期经历了强降水事件，并且在所有55个观测日期中显示出最多的水体数量（765个多边形）。选择检测到最多水体的日期有助于最小化遗漏，并为后续跟踪提供最完整的基础。然而，这种方法可能会错过一些在研究期间仅出现一次的小型或短暂存在的水体。这样的短暂水体无法为时间序列分析提供一致的时间信息，因此预计对分类结果的影响有限。在这个日期检测到的每个多边形都被分配了一个唯一的标识符，从reference_001到reference_765。下一步是跟踪每个参考水体在所有剩余观测日期中的空间位置和表面积变化。

对于每个剩余的日期，将检测到的所有水体多边形与“参考”日期的765个参考多边形进行比较，以建立空间对应关系。例如，在2022年4月5日，检测到324个水体多边形，并暂时标记为001-324。然后检查每个多边形是否与765个参考多边形有重叠。如果有重叠——例如多边形001与reference_034之间有重叠——则将多边形001重新标记为reference_034，并将其面积记录为2022年4月5日reference_034的面积。没有重叠的多边形被排除，因为它们代表参考日期不存在的水体。

然后，每个瓦片通过segment-geospatial模型处理，生成对应于各个水体的实例级别掩膜。生成的256个瓦片级别的掩膜集被拼接起来，重建与W的空间范围一致的全尺寸分割图。这些实例掩膜随后被转换为地理参考的矢量格式（例如GeoJSON）以进行空间分析。由于分块过程，一些跨瓦片边界的水体可能会被分割成多个相邻的掩膜（图3d）。为了纠正这种边缘效应，将对应于同一水体的相邻掩膜在瓦片边界处进行空间合并（图3e），确保在整个研究区域内物体边界的一致性和准确性。一部分划分出的水体多边形与相应的WatNet地表水检测输出和高分辨率卫星影像进行了视觉比较，以评估多边形划分的准确性。大多数主要水体被正确划分，尽管对于小型或破碎的水体以及受到瓦片边界合并影响的多边形仍存在一些不确定性。这些不确定性可能会传递到衍生特征（如水体面积）上。然而，由于受影响的水体通常较小，且边界不确定性一般在几个像素范围内，因此对总体分类结果的影响预计是有限的。

该过程应用于所有55个观测日期，生成了55个文件，每个文件包含特定日期识别出的所有水体的多边形表示。

### 3.4 根据预定义的标准将每个水体标记为人工湖或非人工湖
为了开发一个监督分类模型，需要足够数量的带标签样本。我们将“参考”日期（reference_001到reference_765）的水体多边形导入ArcGIS。首先将这些水体多边形与该地区两个灌溉区（Fresno灌溉区（FID）和Consolidated灌溉区（CID）运营的补给设施地图进行交叉参考，这些设施被分类为人工湖。Fresno大都会洪水控制区（FMFCD）绘制的城市滞洪区被排除在已知的人工湖设施之外。其余的水体多边形通过在Google Earth上进行视觉检查，手动分类为人工湖或非人工湖。为了提高透明度和可重复性，我们基于多个数据源和视觉解释标准制定了一套一致的标记协议（表1）。人工湖位置通常具有矩形形状（1），周围有明显的堤坝（2）。堤坝和盆地区域（干燥时）清除植被（3），以增强补给作用并限制蒸散量。人工湖位置通常由多个沿运河分布的盆地组成（4）。非人工湖类别包括娱乐用水体（如公园或高尔夫球场中的水体）、运河、河流、洪泛区以及天然淹没的土地、滞洪区（类似于人工湖位置，但通常植被更茂盛）、污水处理厂池塘、粪便池（靠近集中动物饲养场）和地下水再充氧池（与地下水井和沙滤池相关）。

### 3.5 构建分类模型以预测人工湖位置
为了区分人工湖和非人工湖水体，我们通过整合水体动态、地形和环境特征来构建一个监督分类模型。工作流程包括四个主要步骤：（1）特征提取，（2）训练-测试划分，（3）模型构建，以及（4）模型评估和优化。以下小节将详细描述每个步骤。

3.5.1 特征提取
我们整理了29个输入特征，并将其分为三类。

第1类：地形特征。对于每个水体多边形，我们提取了三个地形特征：海拔高度，定义为多边形中心的海拔值；平均面积，表示所有55个观测日期中水体的平均表面积；以及海拔面积指数，计算方法为海拔高度与平均面积的乘积。

第2类：环境特征。对于每个水体多边形，我们从美国的水文景观区域（Wolock, 2003）中提取了12个环境变量，以及通过经验回归方程（Reitz, 2017）得出的有效补给量、径流系数和饱和导水性，还有Falcone（2015）提供的人为土地利用信息。计算了多边形中心周围500米半径内的土地利用类型比例，包括耕地、低强度城市、地表水、半开发土地、牧场和已开发/城市化区域。对于每个多边形，河流距离是通过中心网格单元与最近主要河流之间的欧几里得距离来计算的（Faunt, 2009）。

第3类：气象特征。气象特征是从ERA5-Land再分析数据集中提取的。对于以下每个气象变量：降水量、地面以上2米处的空气温度、开放水面蒸发量、植被蒸腾作用产生的蒸发量、潜在蒸发量和NDVI（在第2.2节中定义），我们得出了两种类型的指标：
（1）平均值。对于每一天，计算水体边界内的每个变量的空间平均值，然后在所有55个观测日期上取平均值。对于降水量，为了考虑水文响应的延迟，我们计算了包括观测日期及其前29天在内的30天窗口内的平均降水量，并在每个水体多边形内汇总。使用更短的窗口（例如14天）进行的敏感性测试得到了类似的结果，因此选择30天窗口来表示月度时间尺度上的前驱条件，并提供更稳定的累积湿度指标。
（2）与每个多边形内淹没面积的相关性。对于每个变量，我们计算了55个观测日期内多边形内淹没面积与气象变量之间的皮尔逊相关系数（r），从而得到每个多边形的一个相关值（见公式1）。对于降水量，为了考虑水文响应的延迟，我们计算了淹没面积与30天前驱降水量之间的相关性。

表2总结了所有特征的名称、缩写和类别分配。所有输入特征都使用最小-最大缩放进行了标准化，以确保不同单位和量级的变量之间的可比性。

3.5.2 训练-测试划分
对于每个类别，70%的样本被随机分配到训练集，30%分配到测试集。测试了多种随机种子进行数据划分，模型性能在不同划分下没有显著差异。

3.5.3 分类模型
我们评估了scikit-learn库中八种广泛使用的机器学习分类模型，包括K-最近邻（KNN；Cover, Hart, 1967）、线性支持向量机（Linear SVM；Cortes, Vapnik, 1995）、带RBF核的支持向量机（RBF SVM；Cortes, Vapnik, 1995）、决策树分类器（Decision Tree；Quinlan, 1986）、随机森林分类器（Random Forest；Breiman, 2001）、自适应提升分类器（AdaBoost；Freund, Schapire, 1997）、朴素贝叶斯分类器（Naive Bayes；Langley et al., 1992）和二次判别分析（QDA；McLachlan, 2004）。模型超参数在表3中总结。除非另有说明，所有模型都使用了默认超参数。

表3. 分类模型的参数设置和评估指标。我们报告了以下指标：真阴性（TN）、假阳性（FP）、假阴性（FN）、真阳性（TP）、准确性、精确度、召回率和F1分数。

3.5.4 模型评估和优化
模型性能使用两个标准分类指标进行评估：准确性和F1分数。这些指标用于评估预测标签与真实标签之间的质量。
（1）准确性衡量正确分类实例的总体比例。其计算公式为：
（2）准确性 = （TP + TN）/ (TP + TN + FP + FN)
其中，TP是真阳性（正确分类为MAR的站点）；TN是真阴性；FP是假阳性（错误分类为MAR的站点）；FN是假阴性。更高的准确性表示更好的整体预测性能。
（2）F1分数是精确度和召回率的调和平均值，特别适用于不平衡的数据集。其定义为：
（3）F1分数 = 2 × （精确度 × 召回率）/ (精确度 + 召回率)
其中，精确度 = TP / (TP + FP)；召回率 = TP / (TP + FN)；精确度是正确预测为正的实例（真阳性）在所有预测为正的实例（真阳性和假阳性）中的比例；召回率是正确预测为正的实例在数据集中所有实际正实例（真阳性和假阴性）中的比例。F1分数平衡了这两个方面，在假阳性和假阴性的影响不同时尤为重要。

通过交叉验证和网格/随机搜索对分类算法进行了超参数调整，以最大化F1分数。

4. 结果
4.1 水体分割和标注
为了识别水体，使用六波段遥感图像作为预训练的WatNet模型的输入。图3b展示了WatNet检测水体的一个例子。WatNet有效地捕捉了大规模河流系统，如西北部的圣华金河和西南部的弗雷斯诺沼泽，以及研究区域内分散的小型水体。许多这些较小的水体位于植被茂密或地形复杂的区域，但WatNet仍能高精度地检测到它们，并且回界清晰明确。
为了定义单个水体的边界，使用segment-geospatial框架进一步处理WatNet的输出，生成代表单一水体及其相应边界的多边形（图3a）。该框架有效捕捉了大小水体的形态。例如，图3c中的水体（ID = 79）形状特征保存完好，而附近较小的水体（ID = 82）也被正确识别。
在少数情况下，segment-geospatial将一个水体分割成多个多边形。为了纠正这种情况，增加每个多边形的大小（例如扩大2米），然后识别并合并重叠的多边形，假设它们代表同一个水体（图2d和2e）。这一程序确保了在整个研究区域内每个水体都能唯一且一致地被识别。
最终，在研究区域内共识别出765个水体，包括河流、湖泊和水库。结合ArcGIS、土地利用数据和专家领域知识，我们手动标注了每个水体。其中，139个被标记为MAR站点，626个被标记为非MAR站点。
这些标签用作训练和测试分类模型的真实标签。为了构建训练和测试数据集，MAR和非MAR样本分别随机分为70%的训练集和30%的测试集。因此，训练集包含97个MAR和438个非MAR多边形，而测试集包含42个MAR和188个非MAR多边形。
在图3a中，灰色多边形代表用于分类器模型训练的MAR样本，红色多边形代表用于测试的MAR样本。样本的空间分布是随机化的，并在研究区域内均匀分布，以提高分类模型的泛化能力。

4.2 基线机器学习分类
八种机器学习分类器的分类结果总结在表3中。
在所有模型中，AdaBoost取得了最高的整体性能，准确率为0.90，F1分数为0.70。这种强劲的性能归功于AdaBoost的迭代重加权机制，该机制结合了多个弱学习器并强调错误分类的实例。这种方法在处理类别不平衡和捕捉异构特征空间中的细微模式时特别有效。高的F1分数进一步突显了AdaBoost在平衡精确度和召回率方面的稳健性——这是MAR分类的重要要求，因为假阳性和假阴性都可能导致不理想的水资源规划决策。KNN分类器和RBF SVM分类器也表现良好，准确率分别为0.89和0.90。
相比之下，线性SVM模型的性能最差。该模型没有产生任何真阳性或假阳性，表明完全无法识别任何MAR位置，并倾向于将所有样本分类为非MAR。这表明在该特征空间中，当MAR发生与环境和水文气象特征之间的关系高度非线性时，该模型缺乏足够的区分能力。

4.3 AdaBoost
由于默认的AdaBoost在评估模型中取得了最佳的分类性能，我们选择了它进行进一步的超参数调整和特征选择。超参数调整的目标是开发一个更准确、更稳健的分类器，能够可靠地区分MAR和非MAR水体。

4.3.1 特征选择
我们对一系列超参数值进行了网格搜索。Boost迭代次数（即估计器）分别测试了50、80、100、150、200、250和300次。学习率在0.01、0.05、0.07、0.1、0.2、0.5、0.8和1之间变化。基础估计器的最大深度分别尝试了1、2、3、4、5、8和10个分支。每个配置都使用五折交叉验证进行训练，选择F1分数最高的配置作为最优配置。最佳配置包括200个估计器、学习率1和最大树深度2。在这种设置下，模型在测试集上的准确率为0.93，F1分数为0.81。然而，这导致了相对较高的假阴性率26.2%，表明相当一部分MAR站点被错误分类为非MAR。
上述模型使用了所有29个可用输入特征进行训练。为了进一步提高分类性能并降低假阴性率，我们进行了特征选择，去除了信息量较小的变量，仅保留了对MAR站点预测影响较大的特征。特征重要性得分来自优化的AdaBoost模型，并据此对所有特征进行了排名。结果如图4和图5b所示。在气象变量中，降水量和蒸发量最具影响力。对于环境特征，耕地面积和径流量排名最高。

下载：下载高分辨率图像（163KB）
下载：下载全尺寸图像
图4. (a) 经过超参数调整（包括29个特征）后的AdaBoost分类模型的混淆矩阵，总体准确率为0.93，F1分数为0.81。(b) AdaBoost模型的特征重要性排名。

下载：下载高分辨率图像（264KB）
下载：下载全尺寸图像
图5. (a) 在测试集上评估的包含14个特征的最优AdaBoost分类模型的混淆矩阵。(b) 使用与(a)相同的特征集和超参数配置，模型用不同的随机种子重新训练了300次，以评估准确性和F1分数的变化性。(c) 前N个特征的准确率和F1分数。橙色条形代表准确性，蓝色条形代表F1分数。
为了确定模型构建的最佳输入特征数量，我们测试了一系列包含前n个特征的AdaBoost模型，其中n=10到29。对于每个模型，应用了以下步骤：(1) 根据图4b中的重要性排名选择前n个特征；(2) 在之前描述的相同超参数范围内进行网格搜索；(3) 确定在测试集上具有最高F1分数的配置；(4) 分析最优模型的性能和特征重要性。最终模型使用了14个特征进行训练，这些特征代表了地形、气象和环境条件，如表2所列。该模型获得了最高的测试F1分数0.85和0.94的准确率。该模型包含300个估计器，学习率为0.5，最大树深度为2。假阳性率为3.7%，假阴性率为14.3%。这两种类型的错误发生幅度相当，并且在整个研究区域内空间分布均匀，这突显了模型对于MAR和非MAR样本之间类别不平衡的泛化能力和鲁棒性。

4.3.2 分类评估
为了评估模型的不确定性和鲁棒性，我们使用不同的随机种子和初始权重配置重复了300次训练过程。结果显示运行之间的方差较低，准确率在0.90到0.94之间（平均值=0.92，标准差=0.01），F1分数在0.73到0.83之间（平均值=0.78，标准差=0.04）。这些结果证实了模型的稳定性和可重复性，以及其对随机初始化和数据划分的鲁棒性。然而，模型对非MAR地点的分类准确率略高于MAR地点，这可能是由于标记的MAR样本数量有限。预计未来通过纳入更多的MAR训练数据可以进一步提高性能。尽管如此，当前的结果仍然稳健，并且客观地反映了模型在现有数据限制下的有效性。

4.4 物理评估
机器学习方法在水文研究中的一个常见局限是其“黑盒”特性，这使得物理解释和机制分析变得具有挑战性。这种缺乏可解释性可能会降低科学的透明度，并限制模型输出的可信度。然而，AdaBoost算法提供了特征重要性的内在度量，为解释性提供了一条途径。AdaBoost中的特征重要性来自于每个特征对减少所有弱学习器分类错误的加权贡献。它反映了每个特征对模型准确性的累积改进（Friedman，2001年；Halder等人，2024年）。这有助于识别主要的物理驱动因素，并支持进一步分析影响MAR地点分类的潜在机制。
基于模型中所有特征的重要性排名，我们确定了14个最有影响力的变量（图4b；表2）。这些特征在区分MAR和非MAR区域方面被证明是重要的，并且可以根据它们的时间特性分为两类。
第一类包括时间相关特征，这些特征反映了每个多边形内淹没区域与气象特征之间的时间关系。第二类包括静态或时间平均变量，包括地形、土地利用和气候特征。

4.4.1 第一类：时间相关特征
MAR地点显示淹没区域与气象变量（如降水、NDVI和植被蒸腾蒸发）之间的相关性较弱或不显著（例如，图6a-c）。这种解耦可能与工程的补给基础设施（即渗透池）的存在有关——这些基础设施独立于短期气象变化来稳定水位。在这些系统中，即使在降雨事件期间或之后，由于活跃的渗透、存储操作或重新分配，地表水范围也可能保持不变甚至减少，从而减少了地表水文与气象驱动之间的明显相关性。
下载：下载高分辨率图像（419KB）
下载：下载全尺寸图像
图6. 例如MAR和非MAR多边形的地表水面积和气象变量的时间变化。面板（a）-（c）显示了2022-2023年间55天内一个代表性MAR地点的地表水面积（蓝线，左y轴）和降水、植被蒸腾蒸发以及NDVI（红线，右y轴）的日变化。面板（d）-（f）显示了一个代表性的非MAR地点的相同变量。面板（g）和（h）分别显示了2023年4月6日至4月21日期间MAR和非MAR多边形的地表水范围变化。示例MAR多边形的中心位于36.863°N，119.784°W，示例非MAR多边形的中心位于36.828°N，119.713°W。
相比之下，我们观察到非MAR区域（图6d-f）的淹没区域与气象变量显示出强烈的相关性。这些模式与水文预期一致：降水量增加通常会导致地表水范围扩大，而周围区域NDVI升高和植被蒸腾增强往往与由于蒸发增加和植物水分吸收导致地表水覆盖减少有关。这些相反的效果导致了统计上显著的负相关。
图6g-h展示了这一时期地表水动态的一个代表性例子，说明了2023年4月6日至4月21日之间水体范围的变化。在此期间，如图6a所示，降水量显著下降。在MAR地点（图6g），地表水面积相对稳定，反映了水文缓冲或管理的补给作用。相比之下，非MAR地点（图6h）的地表水范围大幅减少，缩小到大约原来的75%，突出了未管理系统对气象变化的更大敏感性。这些模式与工程水管理系统的行为一致。

4.4.2 第二类：静态或时间平均特征
在第二类特征中，我们研究了两个代表性气象变量的空间模式：平均降水量和平均温度（图7）。选择这些变量是因为它们在分类模型中的重要性评分较高，并且与水文过程直接相关。
下载：下载高分辨率图像（148KB）
下载：下载全尺寸图像
图7. 平均降水（a）和平均温度（b）的空间分布。叠加的多边形表示水体位置。红色多边形代表MAR地点，灰色多边形代表非MAR地点。
如图7a所示，平均降水量从研究区域的西南向东北方向增加。这种梯度受到研究区域东北部山脉的影响。图7b显示平均温度从西南向东北呈下降趋势，这与纬度梯度以及地形效应一致。这些空间模式反映了潜在的气候异质性，并有助于解释分类模型捕捉到的区域补给潜力差异。由于气象数据的空间分辨率相对较低，每个网格单元可能包含多个MAR和非MAR地点。因此，结合额外的环境和地形特征对于补充粗粒度的气象信号和提高MAR地点分类的准确性非常重要。
接下来，我们比较了MAR和非MAR地点之间关键预测变量的分布（图8）。一些变量显示出统计上的显著差异，进一步揭示了与MAR地点相关的物理和地理特征。
下载：下载高分辨率图像（297KB）
下载：下载全尺寸图像
图8. 比较MAR和非MAR地点关键预测变量的箱形图。变量包括：（a）耕地，（b）平均NDVI，（c）平均降水量，（d）平均潜在蒸发量，（e）植被蒸腾平均蒸发量，（f）平均2米空气温度，（g）高程-面积积，（h）径流系数，（i）多边形内的平均淹没面积，（j）到最近主要河流的距离，以及（k）开发/城市化土地。每个箱形图中的五条水平线分别代表最大值、上四分位数、中位数、下四分位数和最小值。每个面板中显示的p值反映了MAR和非MAR地点的分布是否显著不同；p值低于0.05表示统计学上的显著差异。p值是使用Mann-Whitney U检验计算的，该检验不假设各组之间的方差相等。
耕地（图8a）在MAR地点周围500米半径内的值显著较高，这支持了MAR基础设施通常位于灌溉农业区，这些地区对水资源的需求较高且土地管理更为密集的观点。同样，MAR多边形的平均淹没面积通常大于非MAR多边形（图8i）。这可能是由于工程水管理基础设施的稳定作用，使得MAR多边形保持淹没状态，而其他自然非MAR多边形则受到短期气象变化的影响。这种差异也区分了MAR多边形与其他特征，如用于地下水抽取的蓄水池或娱乐池塘。
从水文气象因素来看，MAR地点与较高的平均降水量（图8c）和较高的平均温度（图8f）相关。结合空间地图（图7）的信息，在该区域的东南部发现了较高的平均降水量和较高温度。这些差异可以归因于 Consolidated Irrigation District（CID）运营的设施，这可能反映了该研究区域的特定模式。
从地形角度来看，我们没有发现MAR地点比非MAR地点更靠近主要河流（图8j）。其他两个包含在最佳性能AdaBoost模型中的变量（径流效率或开发面积比例）在MAR和非MAR多边形的并排比较中没有显示出显著差异。
总体而言，该模型不仅实现了强大的分类性能，还成功捕捉到了MAR和非MAR地点之间的有意义的物理区别。这些发现强调了整合地理空间、水文和土地利用特征对于理解现有MAR系统的特征和改善MAR库存开发的重要性。

5. 讨论与结论
我们开发了一个基于深度学习和机器学习的流程来识别MAR地点，整合了多源数据集，包括卫星遥感、气象再分析、环境指标和地形信息。在应用深度学习模型检测和划定55个观测日期卫星图像中的水体后，我们基于14个特征构建了一个监督分类模型。在八个基准机器学习分类器中，AdaBoost取得了最佳性能，在测试集上的准确率为0.94，F1分数为0.85。除了强大的分类性能外，该分类模型在有限标记数据的情况下也展示了稳定的性能。
特征重要性分析揭示了在区分MAR和非MAR地点方面起关键作用的变量。MAR地点通常与耕地和较大的淹没区域相关，这与需要大面积表面区域来支持基于渗透的补给操作一致。与耕地的关联表明，MAR位于地下水支持灌溉农业操作的地区。这些主要预测因子与San Joaquin Valley中密集灌溉农业的支持性水管理基础设施一致，其中工程渠道、水库和补给设施的实施受到灌溉需求和补充地下水含水层的需求的影响。相比之下，由于特征重要性较低，最终模型中没有包括有效补给的估计（Reitz等人，2017年），这表明其他因素在MAR地点的定位中起着更重要的作用。在Central Valley，补给通过降水、河流渗透和灌溉回流进行，管理补给用于补充自然补给并稳定地下水资源。该研究区域的MAR系统实施在渗透性较高的地区，补给效率得到提高以满足地下水需求。分析表明，由于工程水管理系统的作用，MAR地点通常与短期气候变化（如降水、NDVI）脱钩；例如，即使在降水期间，地表水范围也可能保持稳定。相比之下，非MAR水体主要受自然水文过程的控制，因此对气象驱动的反应更强烈。MAR地点并不显著地靠近主要河流。这可以通过广泛的灌溉基础设施来解释，这些基础设施可以将地表水输送到合适的补给地点。
通过评估输入特征的重要性，我们的框架提高了可解释性，实现了模型输出与物理过程之间的联系。我们的结果表明，工程补给系统可能改变了降水量和地表水范围之间的典型关系，这有助于解释为什么这些变量在分类模型中是有用的预测因子。这些关系应被视为对模型行为和特征重要性的解释，而不是独立水文因果关系的证据。此前在中谷地地区的研究主要集中在地下水补给过程（Kocis等人，2017年）、地下水储量变化（Famiglietti等人，2011年）以及人工补给（MAR）适宜性评估（Dahlke等人，2018年）上。然而，由于缺乏全面且空间明确的数据集，区域尺度的MAR设施测绘和分析受到了限制。相比之下，本研究通过多源遥感和机器学习技术来识别当前正在运行的MAR设施，提供了一个数据驱动的框架来绘制现有的MAR操作图谱，并支持数据稀缺地区MAR目录的开发。虽然之前的研究集中在卫星图像像素尺度上识别被淹没的区域（Farooq, Manocha, 2025年），但我们的方法结合了环境特征，能够区分不同类型的水体，并可以扩展到识别卫星图像中的其他人工和自然水体特征，如灌溉水库和自然季节性湿地。该框架经过对现代基础设施的训练后，还可以根据新的卫星图像帮助更新水体特征的多边形地图，例如美国地质调查局的国家水文数据集（NHD）。

该模型是根据加利福尼亚州中谷地地区的数据开发和验证的，其适用于其他水文气候、地质和社会经济背景的情况仍有待评估。由于研究区域内存在空间自相关性和共同的水资源管理特征，所报告的性能应被视为区域内分类性能。在应用于其他区域时，需要考虑环境和水资源管理条件的差异，并可能在应用前进行区域重新校准。此外，由于缺乏官方的MAR设施目录，MAR和非MAR设施的分类依赖于专家解释和辅助数据集；因此，在缺乏独立真实验证数据的地区，可能存在一定的不确定性。尽管如此，这种方法有助于创建连贯的大规模MAR位置数据集，进而可用于研究当前（和潜在的）MAR设施的景观特征。此外，由于数据有限，当前模型中排除了几个地下水文因素，如地下水位深度、含水层储水量和源水质量，尽管这些因素对MAR识别具有重要意义。模型中使用的多个变量，包括饱和水力传导率、径流系数、蒸散量和距离河流的距离，在物理上与渗透能力、径流生成和潜在补给路径有关，因此提供了关于区域尺度上补给相关条件的间接信息（O'Geen等人，2015年）。因此，模型捕捉的地表表达和土地利用特征可以解释为地下条件和人为管理补给操作的综合信号。所提出的框架主要应被视为识别现有运行中的MAR设施的方法，而不是直接评估水文地质适宜性或地下水补给潜力。未来的研究应旨在结合这些水文地质和制度因素以及利益相关者的意见，以便在日益增长的气候和水资源压力下制定更全面和可行的MAR规划。

通过将所提出的流程应用于其他区域，并在必要时进行区域重新校准和重新训练，可以开发出更全面的现有运行中MAR位置的数据集。这样的数据集可以用于分析现有MAR位置的相关地理和地下特征。如果可行，还可以包括机载电磁（AEM）数据作为额外的地下信息来源。结合详细的水文地质数据，得到的MAR目录可以作为未来MAR适宜性和补给潜力评估的基础。专家知识和详细的地下数据（例如AEM）在当地尺度上的MAR相关分析中可能最有价值；在区域尺度上，地形、土地利用和水文地质更为重要；而在更广泛的尺度上，气候（例如年际和年内降水变化）、土地利用和靠近地表水源（例如由山区流域补给的河流）可能是MAR潜力的关键驱动因素。

作者贡献声明：
Indrasis Chakraborty：撰写 – 审稿与编辑、验证、软件、方法论、调查
Gemma J. Anderson：撰写 – 审稿与编辑、项目管理和组织、方法论、调查
Ate Visser：撰写 – 审稿与编辑、监督、资源管理、项目管理和组织、正式分析、概念化
Shuang Yu：撰写 – 审稿与编辑、初稿撰写、可视化、验证、软件、方法论、调查、正式分析
Yang Liu：软件、方法论
Jory Chapin Lerback：撰写 – 审稿与编辑、方法论、概念化