简介面向水文气象与GIS分析人员基于MATLAB实现的泰森多边形法面平均雨量计算工具可依据离散气象站点坐标与雨量观测值通过构建Voronoi图并按多边形面积加权估算区域面平均雨量。包体共3个文件均为.m脚本压缩包仅3KB涵盖主流程、多边形面积求交与泰森多边形构建等核心函数结构精简、便于二次开发。读者可参照其计算思路将站点分布、雨量数据替换为自己的研究对象快速完成面雨量估算与空间权重分析代码中面积求交与泰森多边形生成模块可直接复用适合水文、气象、环境类专业学生及科研人员用于课程设计、论文实验或工程实践。已有921人学习下载内含汉江流域站点数据文件上手门槛低能帮助理解泰森多边形法从建网到加权平均的完整流程。1. 泰森多边形法求面平均雨量从站点观测到流域降雨量的第一课拿到一个流域内十几个雨量站的日雨量要回答“这场雨流域平均下了多少毫米”看起来是个小学算术题但站点在山前和山顶的分布并不均匀直接平均会低估或高估。泰森多边形法Thiessen Polygon按每个站点控场范围切分平面再把面积作权重对站点雨量做加权平均是水文计算里最常用的面平均雨量估算方法之一。它不依赖地形、不依赖显著性检验只要有站点坐标和边界文件就能自动化复现适合做融雪预报、洪水预警、农业墒情里“全流域到底下了多少雨”这类问题。哪怕你已经做过五年以上水文数据也会在看权重、剪边界、处理坐标投影时踩到新坑。2. 泰森多边形法原理垂直平分线与面积权重如何从离散点长出一块控场范围2.1 构造规则平面上的点都归最近站点管泰森多边形的数学名字是Voronoi图。给定一组站点平面被划分成若干多边形每个多边形内任意一点到对应站点的距离都小于它到其他站点的距离。连接相邻站点作连线的垂直平分线这些分界线围成的封闭区域就是一个泰森多边形。所有多边形铺满整个平面。对站点 ( P_i )它的控场区域定义为[ R_i { X \in \mathbb{R}^2 \mid d(X, P_i) \le d(X, P_j), \forall j \ne i } ]这个定义非常反直觉的一点是它完全不考虑地形高低也不管站点海拔差异只依据平面距离。所以流域内有一道山梁雨量站落在山前和山后可能差很多但泰森多边形仍然按最近距离切分山梁处的“真实降雨渐变”并不会被体现。这个特性决定了泰森多边形法适合站网比较密、地形缓变的区域也决定了它必须借助“面积权重”来修正站网稀疏带来的偏差。2.2 面积权重的数学含义不是给站点打分的权重而是给平面分区的面积面平均雨量的计算公式很直接[ \bar{P} \frac{\sum_{i1}^{n} P_i \cdot A_i}{\sum_{i1}^{n} A_i} ]其中 ( P_i ) 是站点雨量( A_i ) 是站点 i 的多边形控制面积。如果所有站点控制面积相同这个式子退化成算术平均如果某站位于流域边缘控制面积可能很大它的雨量就会以更高的权重进入最终结果。注意这里的 ( A_i ) 必须是“流域边界裁剪后的有效面积”不是生成泰森多边形时那片无限延伸的原始多边形面积。假设流域边界是个不规则的狭长山谷流域外一个雨量站控制的区域可能很大但裁剪进流域内的只有一小角权重自然小。很多人在这一步直接用ArcGIS里生成的多边形属性面积忽略了裁剪导致边界站权重虚高最终面雨量偏大。我一般会在计算面积之前先做一次空间求交再统计面积避免这类问题。2.3 泰森多边形法 vs 算术平均 vs 等值线法一张表看出选谁方法核心假设优点缺点适用场景算术平均站点均匀代表全流域计算极简单站点分布不均时偏差大站网均匀、地势平缓等值线法雨量在空间上连续渐变能反映地形和雨量梯度手工勾绘主观、难以自动化有经验的预报员做人工分析泰森多边形法任意点降雨等于最近站点的降雨客观、可编程、面积权重平均忽略地形和雨量渐变站点稀疏误差大自动化报汛、多时段批量计算从表里能看出泰森多边形法本质上是“最近邻插值”的空间分区版本。它比算术平均更尊重站点疏密比等值线法更可复现。对于不需要精细内部结构、只要一个流域总雨量的业务来说泰森多边形是性价比最高的起点。为了让“垂直平分线”这个构造规则直接可操作可以看下面这段用scipy.spatial.Voronoi求站点凸包内分区的最小代码import numpy as np from scipy.spatial import Voronoi # 站点平面坐标经度、纬度已转投影 stations np.array([ [500340.2, 3274100.1], [501002.8, 3274800.6], [499870.5, 3275200.3], ]) vor Voronoi(stations) print(分区原点索引:, vor.point_region)scipy.spatial.Voronoi默认在无限平面构造完整Voronoi图point_region返回每个站点对应的region索引通过vor.regions可以拿到构成该区域的所有顶点编号。但要注意这个图没有边界边缘区域会延伸到无穷远。所以实际做流域泰森多边形还必须用shapely把Voronoi结果裁剪到流域边界而不是直接拿scipy输出当面积用。下一章就解决这个从原始图到流域有效多边形的衔接问题。3. 站点数据准备与泰森多边形生成从CSV坐标到流域面层3.1 输入数据至少要有什么做泰森多边形输入分为空间层和属性层。空间层是站点坐标和流域边界属性层是每个站点对应时段雨量。站点坐标的常见格式是CSVsite_id, x, y S01, 500340.2, 3274100.1 S02, 501002.8, 3274800.6 S03, 499870.5, 3275200.3这里的x, y最好已经是投影坐标如UTM、Albers单位是米。如果源数据只有经纬度要先用pyproj或GDAL转换成投影坐标否则“距离”和“面积”以经纬度为单位没有物理意义泰森多边形的构造也会因经线收敛而变形。流域边界推荐用GeoJSON或Shapefile至少是一个闭合面要素。我会习惯先把边界单独存成一个面层后面所有裁剪都复用它。3.2 不开代码用QGIS点几下生成并裁剪如果你只是临时算一次面雨量不追求脚本化QGIS是最高效的路径。操作步骤导入站点点和流域面两个图层确保二者坐标系一致建议先重投影到同一投影坐标系。在菜单矢量→几何工具→Voronoi多边形打开Voronoi对话框。输入点图层缓冲区设置选“覆盖范围”缓冲区域百分比填0或一个小值比如1。勾选“添加几何属性”可以输出每个多边形的面积字段。运行后得到一个覆盖所有点凸包的泰森多边形图层。用裁剪工具以流域面为裁剪范围对泰森多边形图层做裁剪。注意QGIS的Voronoi工具基于凸包范围生成不会自动把流域外更远处的站点纳入计算。如果流域面积大、周边站稀疏生成的边界多边形在边缘处可能与流域边界相交不准确。此时我会先手动扩大站点输入范围把流域周边50 km以内的站点都放进去再裁剪。3.3 用Python可控生成Voronoi shapely裁剪下面这段代码是“从站点CSV和边界GeoJSON生成流域泰森多边形”的标准流程支持任意投影坐标import geopandas as gpd import pandas as pd from scipy.spatial import Voronoi from shapely.geometry import Polygon from shapely.ops import unary_union # 1. 读站点 st_df pd.read_csv(stations.csv) points [ (row[x], row[y]) for _, row in st_df.iterrows() ] # 2. 读流域边界 basin gpd.read_file(basin.geojson) basin_poly basin.geometry.unary_union # 3. 构造Voronoi vor Voronoi(points) # 4. 提取有限区域多边形并裁剪 polys [] valid_ids [] for idx, region_idx in enumerate(vor.point_region): region vor.regions[region_idx] # -1表示开放区域需要跳过或改写 if -1 in region or len(region) 0: continue coords [ vor.vertices[v] for v in region ] poly Polygon(coords).intersection(basin_poly) if poly.is_empty or poly.area 1e-6: continue polys.append(poly) valid_ids.append(st_df.iloc[idx][site_id]) # 5. 输出面层 result gpd.GeoDataFrame({site_id: valid_ids, geometry: polys}, crsEPSG:3857) result[area_m2] result.geometry.area result.to_file(thiessen_clipped.geojson, driverGeoJSON)逻辑说明vor.point_region中每个站点对应一个区域索引region里的顶点坐标从vor.vertices取。-1代表该区域延伸到无穷远这在流域计算中不合法直接跳过。intersection把无限多边形裁剪到流域内并去掉没交集的区域。最后用GeoDataFrame存储结果crs要和输入站点坐标的坐标系保持一致。参数说明EPSG:3857是Web Mercator适合做显示不适合做面积统计。如果要得到以米为单位的真实面积应该用站点的原始投影坐标系例如中国区常用EPSG:4528CGCS2000 / Gauss-Kruger zone 18或美国地区用EPSG:26910等。代码里的crs可以按需替换但必须与stations.csv里x/y的实际坐标系一致否则后续面积计算会失真。3.4 生成环节的三个关键参数参数作用推荐值/做法坐标系决定距离和面积的真实性统一到站点所在区域的投影坐标系不使用经纬度Voronoi区域是否包含开放区开放区会导致边界外多边形无界跳过region中-1或设置有限包围盒裁剪边界把无限平面收进流域范围使用流域边界的unary_union避免多个子面不合并在做任何面积计算之前先打印一下result[area_m2].sum()和流域边界面积作对比比值在1 ± 0.001内才说明裁剪合理。4. 用Python代码算面平均雨量权重、循环与批量输出4.1 核心逻辑面积权重只算一次雨量序列重复加权泰森多边形法的计算可以分为两步空间计算和时序计算。空间计算把站点坐标和流域边界转化为每个站点对应的有效面积权重这部分只依赖站点布局不依赖降雨数值。时序计算则是把每天的雨量数组按权重做加权平均适合批量处理多年日雨量或逐小时雨量。空间计算在上一章已经完成现在只需要把area_m2转成权重weights result[area_m2] / result[area_m2].sum()然后每个时段的雨量序列如果对应站点顺序一致面平均雨量就是arr rain_values # 长度等于站点数单位mm mean_rain (arr * weights).sum()4.2 完整代码从裁剪好的泰森面层到多时段面雨量表假设你已经有裁剪后的thiessen_clipped.geojson以及一张雨量长表rainfall.csv结构是date, site_id, rain_mm。下面代码计算每天的面平均雨量import geopandas as gpd import pandas as pd import numpy as np # 读泰森面层 thiessen gpd.read_file(thiessen_clipped.geojson) thiessen thiessen.to_crs(EPSG:4528) # 确保投影面积单位 thiessen[weight] thiessen.geometry.area / thiessen.geometry.area.sum() # 读雨量长表 rain pd.read_csv(rainfall.csv, parse_dates[date]) # 把权重映射到站点ID w_map dict(zip(thiessen[site_id], thiessen[weight])) # 按日期分组计算加权平均 records [] for date, grp in rain.groupby(date): grp grp.set_index(site_id) valid grp[rain_mm].dropna() # 只使用有权重值且雨量非空的站点 weights_series pd.Series({sid: w_map[sid] for sid in valid.index}) weights_series weights_series / weights_series.sum() mean_rain (valid * weights_series).sum() records.append({date: date, area_mean_rain: round(mean_rain, 2)}) daily_rain pd.DataFrame(records) daily_rain.to_csv(daily_area_mean.csv, indexFalse)逻辑说明先用投影后的几何面积算出每个站点的权重weight再按日期分组。分组内可能存在某个站当天缺测此时不能机械地使用全局权重而是要把缺失站点的权重剔除把剩余站点权重重新归一化。代码里weights_series weights_series / weights_series.sum()做的正是这件事。参数说明to_crs(EPSG:4528)因人而异关键是确保几何面积单位为平方米。groupby(date)对大时间序列内存友好几千个日期的数据也能跑。输出顺序按日期排序方便直接进入预报流程。4.3 权重重归一化的坑缺测站点多怎么办上面代码在缺测站点较多时会暴露一个假设缺测站的降雨量被“近似看作与该站坐标无关”用剩余站的权重重新分配。这个处理方式在缺测站占比小于30%时是可接受的但如果某天有三分之一站点缺测泰森多边形的“控场范围”实际上已经变了理想做法是为每个时段重新计算只含可用站点的泰森多边形。不过那样计算量大且不稳定业务上更常见的是按固定权重再用缺测站周边站插补。实际项目中我一般会额外输出一列“有效站点数”和“有效权重占比”records.append({ date: date, area_mean_rain: round(mean_rain, 2), n_used: len(valid), eff_weight: round(weights_series.sum(), 4), })eff_weight表示实际参与加权计算的面积权重之和如果是1.0则代表所有站都用上了。当这个值低于0.7时结果的可信度要打折扣需要在报表里标记提醒。4.4 多站点多时段计算的速度优化如果站点500个、时间序列3650天逐日循环调用Pandas的groupby会有一定开销但一般几十秒内能完成。如果再快点可以先把站点雨量变成宽表然后用矩阵乘法wide rain.pivot(indexdate, columnssite_id, valuesrain_mm) wide wide.fillna(0) weight_row pd.Series({sid: w_map[sid] for sid in wide.columns}) weight_row weight_row / weight_row.sum() daily_mean wide.dot(weight_row)注意矩阵乘法里fillna(0)会把缺测当作0这是错的。更安全的做法是把缺测站对应权重置0并重新归一化def safe_dot(row): w weight_row.dropna() w w * row.notna().astype(int) w w / w.sum() return (row * w).sum() daily_mean wide.apply(safe_dot, axis1)row.notna()用于把缺测列权重清零再归一化避免把缺测当0。这种方式虽然用了apply但比纯Python循环快得多时间成本也可接受。5. 泰森多边形法在应用中的边界效应、站点稀疏与验证方法5.1 为什么流域外的站点也要参与生成泰森多边形流域边界不是一个规则的圆当边界上存在“内凹”地形时如果只用流域内站点构造Voronoi流域边界附近可能会出现大片区域不属于任何站点因为该区域到所有流域内站点的距离可能都比到流域外某站远。解决办法是先取流域周边一圈站点一般用缓冲区查询一起生成Voronoi再裁剪到流域内。具体操作buffer basin.geometry.buffer(50000) # 向外50km buffered gpd.GeoDataFrame(geometrybuffer, crsbasin.crs) near_stations gpd.sjoin(st_df, buffered, predicatewithin)buffer(50000)以米为单位需要投影坐标sjoin把缓冲区内的站点全部选出来。之后用这批站点构造Voronoi再裁剪。这样生成的多边形在边界处都会归属于真实最近站不会留下“无主地”。5.2 新增站点权重必须重新算泰森多边形的边界由相邻站点决定。新增一个站点后最靠近它的几个站点控制区域会被重新瓜分面积权重发生变化甚至原本互不相邻的站点也可能因为新站插入而隔开。常见错误是只在站点列表里追加一行雨量却沿用旧的weight字段。正确做法是把新增站点坐标加入站点表重新走第3章的Voronoi生成和裁剪流程再计算新权重。如果站点移动超过一定距离比如500米也应该重建权重。5.3 验证方法交叉检验和与算术平均的对比验证泰森多边形法的结果是否合理最简单的办法是做一个“留一法”交叉检验。假设有N个站每次拿走一个站用剩下N-1个站构造泰森多边形对被拿走的站点位置进行插值得到该位置的面平均雨量估计值与真实值对比。插值方式就是把“被插值点”当作一个临时站点加入Voronoi生成找到它所在的多边形归属的站点雨量然后用那个站点雨量作为估计值泰森法本质是最近邻。或者更简化直接用另一个空间插值方法如反距离加权IDW做对比。在实际业务中我常做的一个检验是把泰森法与算术平均的年总量做对比daily_rain[pet_mean] wide.mean(axis1) daily_rain[diff_ratio] (daily_rain[area_mean_rain] - daily_rain[pet_mean]) / daily_rain[pet_mean]diff_ratio的绝对值如果常年小于0.1说明站点分布均匀两种方法差异不大如果常年大于0.2说明站点空间分布极不均匀泰森法自主调整权重的意义就体现出来了。还可以按季度或汛期统计平均差异检查泰森结果是否在某类环流形势下系统性偏高或偏低。5.4 常见错误排查面积和、空多边形、CRS漂移症状可能原因排查方法裁剪后面积之和比流域面积小0.5%以上部分Voronoi区域被跳过或裁剪边界有裂缝检查vor.point_region中是否有-1区域被误跳过某个站点面积权重为0该站点离流域边界太远裁剪后无交集打印该站坐标到流域边界的最小距离面积单位异常巨大或极小使用了经纬度坐标未投影打印坐标范围确认是不是经纬度值如100~200把thiessen_clipped.geojson在GIS里叠加流域边界看一遍能发现80%的问题。剩下的是逻辑问题比如站点ID在雨量表里和面层里对不上需要先做set(rain.site_id) - set(thiessen.site_id)找差异。6. 进阶技巧把权重矩阵预计算好批量雨量秒级出结果日常做预报会面临“站点没变、每天来新雨量”的场景这时候反复读取GeoJSON和计算面积很浪费。正确做法是把站点权重一次性保存成CSV或JSON之后每次只需要读权重向量与当日雨量向量做点积。# 第一次生成并保存权重 weight_out thiessen[[site_id, weight]] weight_out.to_csv(area_weights.csv, indexFalse) # 之后每次计算 weights pd.read_csv(area_weights.csv).set_index(site_id)[weight] rain_today pd.Series({sid: 12.3, S02: 8.6, S03: 15.2}) # 举例 valid_id rain_today[rain_today.notna()].index w weights[valid_id].copy() w w / w.sum() area_mean (rain_today[valid_id] * w).sum() print(f面平均雨量: {area_mean:.2f} mm)这个思路的关键是权重向量与雨量向量必须按站点ID对齐。如果站点顺序经常变化不要依赖位置索引始终用pandas.Series的对齐机制。另外当出现缺测时仍然要做权重重归一化这一点在批量生成中容易被遗漏。更进一步的实战技巧把面平均雨量连接到GeoDataFrame用不同颜色渲染每个泰森多边形内的降雨量生成面雨量分布图比只看一个平均数字更能发现“哪个站拉高了全流域平均”。输出PNG或GeoJSON之后再做同比或环比分析这算是泰森多边形法在业务里最常见的下钻场景。本文还有配套的精品资源点击获取