做城市体检、设施可达性或公共资源配置研究时很多人会遇到一个看似基础但实际很棘手的问题从 WorldPop、LandScan 等公开产品下载的 1km 网格人口数据空间分布非常漂亮能看出人口集聚的态势可一旦把网格人口按区县边界汇总数字往往和统计年鉴对不上。少的差百分之几多的时候差百分之二三十。这个差异直接影响后续计算人均医院床位数、人均公园面积、15 分钟生活圈覆盖率只要分母一偏结论就全偏了。本文要讲清楚的就是一件事如何用统计年鉴的行政区划人口数据对网格人口数据进行修正让空间分布保持精细同时总量对齐官方口径。核心判断是网格人口数据真正可信的是“相对空间结构”不是“绝对总量”修正的本质是用年鉴总量约束网格数据的权重分配实现“总量一致、分布保持”。读完本文你可以掌握三条修正路径并直接用 Python 跑通一套完整的修正流程。1. 这篇文章真正要解决的问题先说清楚读者群体。如果你正在做这类工作本文对你非常有用使用网格人口数据做空间分析但需要把结果汇总到行政区划进行汇报研究中发现网格人口总量和省、市、县统计年鉴公布的人口差异很大需要在不同年份、不同数据来源之间做人口数据的一致性校正要做高质量的专题地图或人口空间化产品不只是简单展示。很多人一开始的处理方式是“直接用年鉴数据代替网格数据”。这样做的问题很明显统计年鉴只有行政区划尺度的人口总量无法表达人口在区县内部的空间差异。比如一个县的年鉴人口是 50 万人但人口密集区到底在县城、乡镇还是沿交通线分布年鉴数据完全看不出来。另一种处理方式是“直接用网格数据”。这样空间分布是有了但总量可能偏差很大。偏差来源很多产品训练数据年份可能较早模型对不同区域的人口分布估计不同低密度区的反演误差在加总后被放大。于是就会出现“网格汇总 45 万人年鉴公布 55 万人”的矛盾。第三种处理方式也是本文要展开的以年鉴人口为“总量约束”以网格人口为“空间权重”把年鉴人口重新分配到网格单元上。这样修出来的结果既保留了网格数据的空间细节又符合统计年鉴的官方总量。这套思路在学术上对应 dasymetric mapping分区密度制图和 pycnophylactic interpolation质量守恒插值在实际项目中也是最常用的做法。换句话说修正是为了解决三类矛盾空间尺度矛盾年鉴在行政单元上网格在像元上需要一个桥接机制总量精度矛盾网格数据总量不可控需要引入权威总量校准分析结果可靠性矛盾只有口径一致横向比较和纵向时序分析才有意义。2. 核心概念与原理说明2.1 网格人口数据网格人口数据是一个将人口数量按照规则网格通常是 1km × 1km 或 100m × 100m分布的栅格数据集。每个像元的值表示该区域内的估计人口数。常见的公开产品包括WorldPop基于随机森林等机器学习方法结合人口普查、土地利用、夜间灯光等数据生成LandScan由美国橡树岭国家实验室开发使用多源遥感和社会经济数据建立空间概率模型国内也有若干商业和科研机构发布的中国高精度人口网格数据。这类数据的特点是空间分辨率高能表达人口在行政区划内部的不均匀分布。缺点也很明显它是模型估算结果不同产品之间的总量差异可能很大甚至同一产品不同年份之间的可比性也需要验证。2.2 统计年鉴人口数据统计年鉴人口数据是官方发布的、以行政区划为统计口径的人口数据。常见来源包括国家统计年鉴、省统计年鉴、市统计年鉴和县统计年鉴。这类数据的特点是权威、口径统一、年份连续契合行政管理和社会经济分析需求。但它的空间尺度偏粗无法表达行政单元内部的人口分布细节。这里还要留意一个容易混淆的问题统计年鉴中的人口通常有常住人口、户籍人口、现有人口等不同口径。国家和省级统计年鉴多采用常住人口口径而部分县级资料可能采用户籍人口口径。修正前必须确认年鉴口径和网格产品口径是否一致。网格人口产品多数模拟的是“常住类型人口分布”所以优先匹配常住人口口径。2.3 修正的数学本质假设行政单元 i 的年鉴人口为 (Y_i)网格数据在单元 i 内的汇总值为 (G_i)。最简单的比例修正是[ P_j P_j \times \frac{Y_i}{G_i} ]其中 (P_j) 是网格 j 的原始值。如果网格 j 完全位于行政单元 i 内它就按照单元系数进行缩放。这背后的思想是质量守恒。修正后每个行政单元内的网格汇总人口严格等于年鉴人口[ \sum_{j \in i} P_j Y_i ]同时人口在单元内部的相对空间分布保持原网格的结构。这个“保结构”的特性很重要它让修正结果既好看又可用。3. 常见的修正策略对比与选择在具体实现之前先看整体策略。不同策略的差别主要在于修正尺度的层次、是否引入额外辅助数据、计算复杂度。修正策略核心思路优势局限适用场景全局比例修正全省网格汇总值与全省年鉴人口求一个总系数再统一缩放实现简单速度极快无法校正省内各市县的结构偏差只有一个总量数据时应急使用行政区划分层修正按地市或区县分别计算修正系数再逐层缩放省市级总量完全对齐方法直观适合常规项目需要对应级别的行政边界和年鉴数据手头有完整的市县级人口资料推荐首选多级嵌套修正先在省级修正再在市、县两级对上一级结果进行再修正各级总量同时满足空间连续性更好实现复杂度高需要多级边界和数据需要同时满足省、市、县三级口径的项目辅助变量分区密度制图引入土地利用、POI、夜间灯光等数据重新构造空间权重空间分布更接近真实人口态势精度上限更高需要额外数据源模型参数调试成本高科研或高精度人口空间化项目更稳妥的判断是对于大部分业务数据分析场景“行政区划分层修正”性价比最高。它不需要额外数据逻辑容易向非技术团队解释而且修正后结果在多级行政单元上都有明确含义。如果你在做一个省级人口空间化产品建议采用“省级→市级→县级”逐层嵌套修正。4. 环境准备与前置条件4.1 运行环境本文示例基于 Python 3.9 及以上版本使用 geopandas 处理矢量行政边界rasterio 处理网格人口栅格pandas 处理统计年鉴表格数据。建议用 conda 创建一个独立环境避免与已有项目冲突。# 文件路径environment.yml name: pop-correction channels: - conda-forge dependencies: - python3.9 - geopandas0.14.3 - rasterio1.3.9 - pandas2.1.0 - openpyxl3.1.2 - jupyter创建并激活环境conda env create -f environment.yml conda activate pop-correction以上版本号仅为示例。实际使用请以你本机 Python 环境和库版本为准。核心 API 在 rasterio 1.3 和 geopandas 0.14 上都已验证可用版本相近的旧版本通常也能运行。4.2 数据准备建议把所有数据集中在一个项目目录中pop_correction/ │-- data/ │ ├── census.xlsx # 统计年鉴人口数据 │ ├── boundary.shp # 行政区划边界 │ └── pop_grid.tif # 网格人口数据 │-- scripts/ │ └── correct_pop.py └── output/ └── pop_corrected.tif # 修正后的栅格需要检查三件事行政区划边界的坐标系必须和网格人口栅格的坐标系一致。如果不一致先统一投影。推荐使用 Albers 等面积投影或 UTM 投影因为人口密度计算对面积准确性敏感。栅格数据的 NoData 值必须明确。WorldPop 等产品一般用 0 表示无人区但部分产品用 -9999需要在读取时确认。统计年鉴表格中必须包含可关联的行政区划代码或名称且和边界矢量中的对应字段值完全一致。5. 核心实现从汇总到修正这一部分给出可以直接跑通的代码。示例以“某省所有区县”为处理单元目标是用统计年鉴中的区县人口修正现有网格人口。代码中所有路径、字段名都按常见情况编写你需要替换成实际数据中的字段。5.1 读取统计年鉴与行政区划# 文件路径scripts/correct_pop.py import pandas as pd import geopandas as gpd import rasterio import numpy as np from rasterio.features import rasterize # 1. 读取统计年鉴数据 # 假设 Excel 中包含三列ADM_CODE行政区划代码、NAME名称、YEARBOOK_POP年鉴人口 yearbook pd.read_excel(data/census.xlsx, sheet_name0) print(yearbook.head()) # 2. 读取行政区划边界 boundary gpd.read_file(data/boundary.shp) print(boundary.crs) print(boundary.columns.tolist()) # 3. 读取网格人口数据的属性信息 pop_src rasterio.open(data/pop_grid.tif) pop_array pop_src.read(1).astype(np.float64) print(栅格形状:, pop_array.shape) print(栅格 CRS:, pop_src.crs) print(NoData:, pop_src.nodata)这段代码做了三件基础工作把年鉴读成 DataFrame把行政边界读成 GeoDataFrame把网格人口读成 numpy 数组。真正容易踩坑的地方是字段名对齐。很多统计年鉴下载下来的 Excel 里带有表头、说明行、合并单元格直接读取会导致数据类型混乱。建议先用 Pandas 清理只保留行政区划代码和人口数字两列并且把代码列转为字符串再与边界数据进行关联。5.2 将行政边界栅格化建立像元与区县的对应关系后续修正的关键是要知道每个像元落在哪个区县内。最简单高效的方式不是逐区县裁剪栅格而是把所有区县一次性栅格化生成一张“区县 ID 栅格”。# 4. 构建区县 ID 栅格 # 假设 boundary 中已有 ADM_CODE 字段yearbook 中也有 ADM_CODE # 先按 ADM_CODE 排序保证一一对应 boundary boundary.sort_values(ADM_CODE).reset_index(dropTrue) yearbook yearbook.sort_values(ADM_CODE).reset_index(dropTrue) # 为每个区县分配从 1 开始递增的 ID boundary[ZONE_ID] np.arange(1, len(boundary) 1) # 栅格化把区县矢量转成与人口栅格同一形状的 ID 栅格 zone_id_raster rasterize( shapes[(geom, zone_id) for geom, zone_id in zip(boundary.geometry, boundary.ZONE_ID)], out_shape(pop_src.height, pop_src.width), transformpop_src.transform, fill0, # 0 表示不在任何区县内的像元 dtypeint32, )这一步的价值在于把复杂的“栅格-矢量叠加”问题转化成了批量数组运算。之后对每个像元通过ZONE_ID就能找到它所属的区县再通过区县索引找到对应的年鉴人口和网格汇总人口。5.3 计算每个区县的网格汇总人口和修正系数使用np.bincount按区县 ID 快速汇总人口可以避免逐区县循环裁剪的耗时操作。# 5. 计算每个区县的原始网格人口汇总 flat_zone zone_id_raster.ravel() flat_pop pop_array.ravel() # bincount 返回第 i 个位置的数值表示 ZONE_IDi 的所有像元人口之和 zone_sum np.bincount(flat_zone, weightsflat_pop, minlengthlen(boundary) 1) # 6. 构造每个区县的年鉴人口数组 zone_target np.zeros(len(boundary) 1) for i, row in boundary.iterrows(): zone_target[i 1] yearbook.loc[yearbook[ADM_CODE] row[ADM_CODE], YEARBOOK_POP].values[0] # 7. 计算修正系数防止除零 zone_ratio np.divide( zone_target, zone_sum, outnp.ones_like(zone_target, dtypenp.float64), wherezone_sum 0, )zone_sum和zone_target都是长度为 N1 的一维数组N 是区县数量。下标为 0 的位置留给边界外像元下标从 1 到 N 对应每个区县。np.divide中的where参数确保人口汇总为 0 的区县不会被除零而是保留系数 1。5.4 生成修正后的栅格并输出# 8. 应用修正系数到每个像元 corrected_flat flat_pop * zone_ratio[flat_zone] corrected corrected_flat.reshape(pop_src.shape) # 9. 输出修正后的 GeoTIFF output_path output/pop_corrected.tif with rasterio.open( output_path, w, driverGTiff, heightpop_src.height, widthpop_src.width, count1, dtypefloat32, crspop_src.crs, transformpop_src.transform, nodatapop_src.nodata, ) as dst: dst.write(corrected.astype(np.float32), 1) print(修正结果已保存:, output_path)核心逻辑是flat_pop * zone_ratio[flat_zone]。flat_pop是原始每个像元的人口zone_ratio[flat_zone]是该像元所属区县的修正系数。乘法完成后每个区县内的像元都乘上了同一个系数因此区县内汇总人口会自动等于年鉴值。如果你有地级市和区县两级数据可以把上面的过程封装成函数按“市级修正→区县级再修正”的顺序执行两次。第二次修正时输入栅格是第一次修正后的结果输入边界是区县级边界这样最终结果能同时满足市和县两级总量。5.5 扩展引入辅助变量的分区权重修正如果手头有土地利用分类、夜间灯光或 POI 密度数据可以升级为分区密度制图。思路是在区县内部不直接用原始网格人口作为权重而是用“辅助变量权重”作为空间分配依据。# 假设 land_weight.tif 是一个权重栅格城镇像元权重为 3乡村像元权重为 1 land_src rasterio.open(data/land_weight.tif) land_weight land_src.read(1).astype(np.float64) # 计算每个区县内的权重总量 flat_weight land_weight.ravel() zone_weight_sum np.bincount(flat_zone, weightsflat_weight, minlengthlen(boundary) 1) # 修正系数该区县年鉴人口 / 该区县权重总量 zone_density_ratio np.divide( zone_target, zone_weight_sum, outnp.ones_like(zone_target, dtypenp.float64), wherezone_weight_sum 0, ) # 新人口 权重 * 系数总量对齐后分布依权重变化 dasymetric_pop land_weight * zone_density_ratio[flat_zone].reshape(land_weight.shape)这种方式适合做高精度人口空间化研究场景。它牺牲了原始网格数据的人口分布结构换取了一个更接近“真实人口承载空间”的分布结果。工程项目的上限更高但也需要更谨慎的权重设定并且要对每一类辅助数据做质量评估。6. 运行结果与效果验证6.1 运行脚本在项目根目录执行python scripts/correct_pop.py如果一切正常output/pop_corrected.tif会生成修正后的人口栅格。6.2 验证修正效果修正是否成功不能只看“有没有输出文件”必须做定量验证。验证的核心指标有两个修正后各区县人口汇总与年鉴人口的相对误差修正前后各区县内部的空间人口格局是否保持。# 验证脚本 corrected_src rasterio.open(output/pop_corrected.tif) corrected_array corrected_src.read(1).astype(np.float64) flat_corrected corrected_array.ravel() corrected_zone_sum np.bincount(flat_zone, weightsflat_corrected, minlengthlen(boundary) 1) check_df pd.DataFrame({ ADM_CODE: boundary[ADM_CODE], NAME: boundary[NAME], YEARBOOK_POP: [zone_target[i 1] for i in range(len(boundary))], GRID_RAW: [zone_sum[i 1] for i in range(len(boundary))], GRID_CORRECTED: [corrected_zone_sum[i 1] for i in range(len(boundary))], }) check_df[ERROR_RAW] (check_df[GRID_RAW] - check_df[YEARBOOK_POP]) / check_df[YEARBOOK_POP] * 100 check_df[ERROR_CORRECTED] (check_df[GRID_CORRECTED] - check_df[YEARBOOK_POP]) / check_df[YEARBOOK_POP] * 100 print(check_df.head(10)) print(修正前最大误差: {:.2f}%.format(check_df[ERROR_RAW].abs().max())) print(修正后最大误差: {:.2f}%.format(check_df[ERROR_CORRECTED].abs().max()))预期结果修正后每个区县的GRID_CORRECTED都和YEARBOOK_POP完全相等ERROR_CORRECTED接近 0。如果出现非 0 误差通常是浮点精度问题误差会在 1e-6 量级。还需要检查空间格局。最直接的检查方式是把修正前后的栅格在 GIS 软件中叠加显示观察人口高值区是否出现在原来的高值区低值区是否仍然保持低值。更严格的做法是计算修正前后所有非零像元的 Spearman 秩相关系数一般要求大于 0.9。如果显著下降说明修正过程破坏了原始数据的空间结构。from scipy.stats import spearmanr mask (pop_array 0) (corrected_array 0) (flat_zone.reshape(pop_src.shape) 0) corr, _ spearmanr(pop_array[mask], corrected_array[mask]) print(修正前后像元秩相关系数: {:.4f}.format(corr))如果运行失败第一步先看两个位置zone_sum中是否有大量 0。如果有说明边界和栅格没有套合上优先检查 CRS年鉴数据和边界数据关联后是否存在空值。用yearbook.isnull().sum()检查常见原因是行政区划代码前后年份不一致。7. 常见问题与排查思路网格人口修正看起来简单实际操作中问题往往出现在数据对接和尺度转换阶段。下面整理了几个高频问题。问题现象可能原因排查方式解决方案修正后区县汇总人口与年鉴差异仍然很大行政区划边界与年鉴统计口径不一致栅格像元没有落到边界内打印每个区县的 ZONE_ID、zone_sum、zone_target统一边界数据来源和年份重新检查关联字段栅格化后大量像元 ZONE_ID 为 0边界要素或坐标系有问题在 GIS 中叠加显示边界和栅格检查 boundary.crs 与 pop_src.crs将边界和栅格重投影到同一坐标系建议使用等面积投影某区县原始网格汇总人口为 0区县面积很小但栅格分辨率低或栅格 NoData 值处理错误查看该区县边界内的原始像元值和 NoData 设置将 NoData 转为 0对极小区域采用更高分辨率辅助数据人口高值区出现不合理尖峰原始网格本身在某个像元存在异常高值比例修正直接放大它绘制修正前后分位数分布图先对原始网格做分位数截断或空间平滑再执行修正不同口径人口常住、户籍混用年鉴字段含义不清楚查看统计年鉴表头说明统一采用常住人口口径如无法统一在报告中显式注明修正后结果与邻域省份边界人口突变过大省级边界两侧分别做了独立的比例修正检查省界两侧人口密度变化如研究跨省建议以全国为整体统一全局比例处理后再做省级分层修正8. 最佳实践与工程建议结合实际项目经验这里整理几条容易被忽视但非常重要的建议。8.1 先做差异诊断再决定修正方案拿到数据后第一步不要急着写修正代码。先计算每个区县的“网格汇总人口 / 年鉴人口”系数并绘制地图。如果系数空间分布呈现明显的区域聚集性比如东部系数普遍高于西部说明网格产品存在系统性偏差分层修正的收益会很大。如果系数随机分布且都接近 1说明原始网格总量已经比较准确简单的全局修正即可不需要做复杂处理。8.2 行政区划代码是命脉统计年鉴中的行政区划代码和边界矢量中的代码必须完全一致。但在实际项目中经常遇到县级市升级、代码变更、统计口径调整的情况。建议把所有代码统一处理为字符串并做一次全量比对assert set(yearbook[ADM_CODE]) set(boundary[ADM_CODE]), 行政区划代码不一致如果存在不一致优先以边界文件为准将年鉴数据按名称或新旧代码映射后重新对齐。8.3 保留原始数据备份修正过程可重复修正过程本身就是一个数据加工流程必须保证可复现。建议原始网格人口、原始年鉴、原始边界分别只读不做原地修改修正代码统一放入脚本目录输入输出路径使用参数配置在输出栅格的元数据或配套 README 中记录数据来源、下载日期、修正方法和代码版本。8.4 修正后的栅格仍需要做质量报告最终交付时建议附带一个质量说明文件包括各区县修正前后误差表、修正系数分布直方图、像元秩相关系数。这不仅能帮助他人理解结果也是保护自己成果的重要方式因为网格人口修正后的产品经常会被用于规划建议、论文发表和项目汇报必须有可信的精度描述。8.5 注意生产环境的计算效率如果你的研究范围是全省、全国网格可能达到上亿像元。np.bincount虽然快但仍需要确保内存足够。建议先做小范围测试确认流程无误后再对全量数据执行。必要时可以用rasterio的窗口读入机制分块处理大栅格。8.6 不要盲目追求“完全对齐”比例修正可以做到区县总量完全对齐但如果继续对齐到乡镇、街道级误差和数据需求量都会大幅上升。业务上要明确修正的精度目标。如果手中没有可靠的乡镇级年鉴数据强行把比例修正做到乡镇级只会放大噪声。9. 总结与后续学习方向网格人口数据与统计年鉴人口之间的矛盾本质是“空间精细度”和“统计权威性”之间的矛盾。本文给出的核心解法是以统计年鉴为总量约束以网格数据的空间分布为权重通过质量守恒方式进行再分配。其中最简单实用的方案是行政区划分层比例修正它既能做到区县总量对齐又能保留原始网格的空间结构。建议你先用一个小区域跑通完整流程再扩展到全量数据。下一步可以沿着两个方向继续深入一是引入更多辅助数据比如土地利用、POI、夜间灯光和高分辨率建筑数据把简单的比例修正升级为分区密度制图模型二是把修正流程封装成可复用工具把数据读取、系数计算、栅格输出、质量报告串成一条自动化流水线。如果你的项目中还存在人口数据年份不一致、边界精度不足、乡镇级数据缺失等问题建议在动手修正之前先把数据和口径问题整理成清单逐项确认后再执行算法。毕竟人口数据的每一步处理都直接关系到后续规划和评估结论的可靠性值得充分验证后再交付。