做水文分析或者区域资源评估的朋友肯定绕不开华北地区的水资源问题。这一片地方降水在时间和空间上分布极不均匀地下水又是长期超采的状态很多研究工作到最后都会卡在数据这关。我整理了一套2005–2016年间、围绕水资源承载力主题的华北地区降水与地下水要素数据这段时间恰好覆盖了华北平原地下水变化最剧烈的时期也是水资源承载力相关研究的关键阶段。这篇内容把我整理数据时的思路、数据结构的设计、实操处理流程以及踩过的坑全部摊开来说给做相关方向的朋友一份可以直接参考的实践笔记。1. 这套华北降水与地下水数据为什么值得用1.1 数据背后的现实需求常住华北地区的人对水资源的紧张程度都有体感但做分析的人需要更精确的判断。水资源承载力这个概念通俗讲就是一块区域在不破坏生态系统、不造成地下水持续下降的前提下能支撑多少人口和多少经济活动。而降水量和地下水要素恰好是这个计算体系中最硬的两个输入参数。降水是地表水和地下水的补给来源地下水又是华北很多城市的主要供水水源两者必须放在一起看。长期降水偏少导致补给不足地下水开采量又持续大于补给量水位自然就逐年下降。分析这一段时期的降水与地下水数据不只是算一个平均值更重要的是看清变化趋势、响应关系和空间差异这是承载力评估的前提。这套数据把2005到2016一共12年的降水序列和地下水埋深、储量变化等要素做了统一整理可以支撑从单点时间序列分析到区域空间制图的各种用途。1.2 2005—2016年这个时段为什么特殊选择这个时间窗口不是随意的。2005年前后华北地区地下水监测网络已经相对完善实测井数据质量比较稳定。而在2016年以后南水北调工程通水的影响开始逐渐显现地下水开采格局发生了明显变化如果从承载力研究的角度讲2005到2016正好是一段“超采压力持续积累、水位持续下降”的代表性时期用来研究自然降水和地下水的响应关系非常合适。在这个12年里年降水量波动明显有些年份偏涝、有些年份偏旱。地下水面对降水的响应往往存在滞后所以数据整理的时间跨度至少要覆盖一个完整的干湿波动周期才能捕捉到埋深对降水的滞后响应。这也是我选择从2005年而非更早年份开始的原因之一再往前早几年的站点数据完整性会差很多处理起来需要大量插补对精度影响较大。1.3 数据适合谁来用、能解决什么问题这套数据对以下几类人比较有用做区域水资源承载力评价的科研人员、做农业灌溉需水分析的农业气象工程师、做国土空间规划中水资源约束条件分析的设计人员以及做相关课程设计或毕业论文的学生。比起自己去各个平台找散装数据再手工对齐这套数据最大的价值在于统一了时间口径、空间范围和单位体系。拿到手之后可以直接进行趋势分析、空间对比、或者作为模型输入的驱动数据省掉了大部分数据清洗的时间。2. 数据要素解析降水与地下水要素到底包含什么2.1 降水要素的层级结构降水数据不是一个单一年值那么简单我整理了三个层级。第一个层级是月尺度降水序列这是最基础的数据。华北地区降水主要集中在7到9月月尺度数据能清楚反映季节性差异也为后续计算年内分配不均匀系数提供了基础。第二个层级是年尺度降水总量和生长季降水总量年降水量是判断当年干湿状况的核心指标。第三个层级是降水距平序列以2005到2016年的平均值为基准计算逐年的偏离程度这样能直观看出哪几年属于少水年、哪几年属于多水年。降水数据在原始状态下通常是站点观测值分布在气象站和水文站上而我整理的数据集中提供的是经过空间插值后的栅格数据分辨率统一为0.05度大约5公里左右这个精度在区域尺度分析中够用又不会因为分辨率过细而产生大量假精度问题。2.2 地下水要素的主要构成地下水这部分内容比降水复杂得多。首先是地下水埋深数据也就是地下水面到地面的距离单位是米。这是衡量地下水超采程度最直接的指标埋深越大说明水位越深开采漏斗越明显。华北平原很多区域在2010年以后埋深已经超过30米甚至一些集中开采区超过50米。其次是地下水储量变化数据我综合了实测井数据和重力卫星反演的结果。GRACE重力卫星可以监测大尺度地下水储量的时空变化趋势以等效水高表示单位是厘米。实测井数据反映的是点位水位变化卫星数据反映的是区域储量变化两者结合能弥补空间代表性不足的问题。另外还有一组衍生数据降水入渗补给估算量和地下水可开采系数。前者是基于降水、土壤质地、包气带岩性估算的后者是实际开采量与可开采量的比值当系数大于1时说明处于超采状态。这两个指标把降水和地下水的关联性串了起来是水资源承载力计算中最常用的中间参数。2.3 数据文件组织和字段命名规范数据文件的组织看起来是个小事情但实际使用时差别很大。我把数据整理成三类文件逐年栅格文件GeoTIFF格式命名格式为pre_year.tif、gwd_year.tif统计表格CSV格式包含兴趣点或行政区的逐年降水、埋深、储量变化数据元数据说明文件PDF格式说明数据来源、精度、投影信息和单位字段命名的核心原则是自解释。比如降水量字段命名PRCP_MM地下水埋深字段命名为GWD_MGRACE储量异常命名为GRACE_CM_WATER。看到字段名就能直接知道它的含义和单位不需要反复查说明文档。字段名含义单位时间粒度PRCP_MM降水量毫米月/年PRCP_ANOM降水距平毫米年GWD_M地下水埋深米月/年GRACE_CM储量异常厘米月RECHARGE_MM降水入渗补给量毫米年OVERDRAFT_RATIO地下水开采系数无量纲年3. 数据获取与预处理实操从原始资料到标准数据集3.1 基础资料来源与选取逻辑降水数据的原始来源主要是国家级地面气象观测站日值数据这部分比较容易获取常规的数据共享平台都提供。但这里有一个关键问题华北地区面积大站点数量虽然不少但分布并不均匀西部山区站点稀、东部平原站点密直接使用站点数据进行区域统计会产生偏差。地下水数据的原始来源相对复杂一些不同省份的监测井编号体系不统一测量频率也有差异。有的井是逐月测量有的井是一年测两次整理时需要进行时间统一。另外开采量的数据往往统计到县级行政单元要和空间数据关联还需要做行政区编码匹配。所以我在整理时采用了“站点数据优先、遥感数据补充、统计年鉴校验”的三源思路。水文年鉴和地下水通报用于校验总量级GRACE数据用于校验趋势级尽量确保数据不仅数值对得上变化趋势也有依据。3.2 空间插值处理流程从站点观测值到栅格空间分布这一步是整个流程中需要谨慎操作的地方。我采用的是协同克里金方法将海拔作为协变量。华北地区降水受地形影响较明显太行山迎风坡降水明显高于山前平原如果只用普通克里金插值而不考虑海拔燕山和太行山一带的降水会被严重低估。处理流程大致如下第一步将所有站点数据的坐标系统一转换为CGCS2000地理坐标系第二步提取每个站点的海拔信息第三步在Python中调用相关空间插值库设置变异函数模型为指数模型搜索半径设置为基础变程的1.5倍第四步输出0.05度分辨率的栅格再统一裁剪到华北地区研究边界第五步和独立验证站点数据进行交叉验证计算平均误差和均方根误差。注意交叉验证时不能用参与插值的站点直接验否则误差会系统偏低看起来精度很好实际使用时会出问题。需要预留大约15%的站点作为验证集等插值完成后再用验证集数据评估。地下水埋深的插值逻辑和降水不太一样。地下水埋深的空间变化受开采格局和含水层结构的影响更大地形因素相关性弱。我直接采用普通克里金插值同时参考已有的水文地质分区图作为约束在插值结果中对不同水文地质子区设置了允许范围避免出现明显的异常突变。3.3 时间序列对齐与异常值处理时间对齐是数据整理中最容易出问题的地方。降水数据涉及闰年2月29号的处理地下水数据涉及不同月份的测量缺测如果处理不当时间序列会出现错位。我在处理降水数据时将所有的日期统一为公历日并对2月29日的数据进行了独立处理在计算多年平均时剔除了这部分避免对月均值产生干扰。地下水埋深数据方面部分监测井存在连续数月的缺测。我没有直接采用简单的线性插值而是结合附近同类监测井的变化速率进行修复。具体做法是找到待修复井周边同一含水层组的参考井计算出它们在缺测期间的井水位变化速率再应用到待修复井上这样可以更好地保持井点本身的动态特征。3.4 基于Python的批量处理示例这里提供一个简化版的处理思路方便实际动手操作。import rasterio import numpy as np import pandas as pd # 读取逐年降水栅格计算区域平均值 def zonal_mean_precip(raster_path, region_mask): with rasterio.open(raster_path) as src: data src.read(1) profile src.profile # 将无数据区域设为NaN data[data profile.get(nodata, -9999)] np.nan # 应用研究区域掩膜 data_masked np.where(region_mask 1, data, np.nan) return np.nanmean(data_masked) region_mask np.load(north_china_mask.npy) years range(2005, 2017) monthly_data [] for year in years: for month in range(1, 13): raster_file fpre_{year}_{month:02d}.tif mean_val zonal_mean_precip(raster_file, region_mask) monthly_data.append({year: year, month: month, pre_mean: mean_val}) df pd.DataFrame(monthly_data) # 生成年汇总 annual_precip df.groupby(year)[pre_mean].sum() * 0.05 * 0.05 * 111 * 111 / 100这段代码的思路非常简单就是批量读取栅格使用研究区掩膜过滤然后做区域平均。需要注意的是最后乘以栅格单元面积时要考虑投影方式。0.05度在赤道附近对应的地面距离和在高纬度地区是不同的如果直接在华北区域用同一个固定系数会带来一定误差。严谨的做法是先计算每个栅格单元的实际面积。4. 基于这套数据的水资源承载力评估思路4.1 水量平衡与承载力计算框架有了降水和地下水数据之后承载力评估的计算框架就顺理成章了。基本思路是从水量平衡出发把区域可利用水资源量和人类用水需求放在一个框架里比较。公式可以简化为区域可利用水资源量 当地地表产水量 × 可利用系数 地下水可开采量 外调水可用量其中地表产水量主要受降水控制通过降水减去实际蒸散发加上地表径流输出得到地下水可开采量则通过可开采系数乘以地下水资源总量获得。承载人口规模的计算公式为承载人口规模 区域可利用水资源量 / 人均综合用水定额华北地区的人均综合用水定额一般取200到300立方米每人每年包含生活、生产和生态用水。当然这个定额要根据具体区域的产业结构调整农业占比高的区域定额会明显偏大。我在实际测算中发现用这套框架算出来的承载力和实际常住人口数量之间往往存在不小差距这正是“超载”的直接体现。当测算承载力低于实际人口规模时说明该区域在不增加外调水的情况下水资源系统已经承受了额外压力。4.2 降水入渗补给估算降水入渗补给是连接降水和地下水的桥梁。降水落到地面后一部分形成地表径流流走一部分被土壤吸收蒸发还有一部分入渗补给含水层。降水入渗补给量一般通过入渗系数法估算Q入渗 P × α × A其中P是年降水量A是计算区域面积α是降水入渗补给系数。α的取值不是固定的它和包气带岩性、地下水埋深、降水量大小都有关系。华北平原的浅层粉质黏土区域α一般在0.10到0.18之间砂质土壤区域稍大一些可以达到0.2以上。需要注意的是当埋深很浅时降水入渗在短期内可能表现为水位上升显著但当埋深很深时包气带厚度大入渗的补给路径也长响应会有显著的滞后。在实际计算时我会先用埋深数据对整个区域进行分区不同埋深区间采用不同的α值这样估算结果会比用一个固定系数更合理。4.3 地下水超采程度判断判断一个区域有没有超采最直接的指标是地下水位动态。如果地下水埋深在较长时期内持续增加比如连续5年以上的整体下降那基本可以判定开采量大于补给量。但我更建议用开采系数来判断超采状态而不是只看水位。因为水位下降除了受开采影响还和降水偏少、补给减少有关。开采系数是实际开采量和可开采量之间的比值开采系数 实际开采量 / 可开采量当开采系数大于1时处于超采状态。华北地区很多县市在这个时段内的开采系数常年维持在1.2以上个别集中开采区甚至达到2左右。把开采系数和降水距平放在一起看能发现很有意思的现象连续两个枯水年之后超采程度往往明显加剧这和开采结构调整有直接关系。5. 常见问题与排查技巧实录5.1 数据时间序列对不齐怎么办实际使用中最常遇到的问题是降水数据有12个月记录但地下水数据因为某些月份缺测导致时间序列不完整。直接删掉缺测月份会让年统计值失去可比性。我的处理方式是如果某年缺测不超过2个月用该井的多年均值插补超过2个月则对这一年标注为数据不足计算时不参与趋势检验。不建议用回归插补强行补齐太多月份这会让数据本身的信息失真。5.2 空间分辨率不匹配怎么统一降水栅格和地下水栅格有时候分辨率不一样比如降水数据是0.05度地下水数据原始可能是1公里分辨率。直接重采样会导致边界模糊或者信息丢失。统一做法是以项目中用途更精细的那个数据为基准。如果是做承载力分区计算建议保持行政区尺度先将各栅格做区域统计得到区县平均值再参与计算尽量避免在栅格层面做重采样。因为重采样在区域均值上变化不大但在极值区域比如超采漏斗中心会明显拉低峰值。5.3 坐标系统不一致导致的偏移现象华北地区的数据来源多样有的用WGS84有的用CGCS2000如果混用而不做转换在区域分析时会出现几百米的偏移。几百米在高精度站点对比时可能造成站点落在隔壁栅格里影响提取值。排查方法很简单将所有数据加载到GIS软件里叠加行政边界查看边界处栅格像元是否对齐。如果不齐统一转换为CGCS2000。转换后再检查一次。5.4 常见问题的速查速解表问题表现可能原因处理方法降水栅格一片空白投影信息缺失、裁剪范围错误检查投影重新裁剪地下水埋深值异常偏大单位混淆厘米写了米核对原始资料单位区域年降水序列和年鉴对不上栅格单元面积计算错误按每栅格实际面积加权GRACE数据趋势和实测井相反使用了不同版本解算结果统一使用同一机构版本插值结果出现中间凹陷克里金搜索半径过小增大搜索半径并检查变程5.5 避坑心得别迷信高分率在处理数据的过程中我最大的体会是不要盲目追求高分辨率。很多时候数据源本身的精度就有限比如GRACE的原始分辨率其实很粗强行重采样到0.05度并不会提高真实精度只会让影像看起来更细腻而已。区域尺度的水资源承载力分析核心在于趋势判断和相对空间差异而不是单个栅格的绝对数值。一套分辨率稍低但时间连续、空间覆盖完整的数据比一套分辨率高但到处填补的数据实用得多。另外做水资源承载力分析时建议始终保留一份原始数据备份。所有预处理步骤都基于备份操作不要直接在原始数据上改。预处理脚本要保留这样才能保证分析过程可以被回溯数据更新时也能快速重新生成新版本。这套数据用下来尤其是在分析降水对地下水的补给响应时我确实感受到区域水循环系统的复杂性。两年的降水变化看不太出明显趋势放到10年以上的尺度才能看清规律。数据整理的功夫很多时候花在看不见的地方但最后分析结果的可靠性恰恰就来自这些基础工作。