环境健康数据分析实战：从PM2.5暴露到归因死亡数的Python计算流程

在环境健康与公共卫生领域，细颗粒物污染对全球疾病负担的影响一直是研究热点。近期一项覆盖全球范围的研究指出，超细颗粒物每年可能导致近200万例过早死亡，这一结论再次将公众视线聚焦于空气污染的微观危害。对于从事环境数据分析、公共卫生政策研究或相关领域开发的工程师和研究者而言，理解这一结论背后的数据来源、分析方法和潜在的技术实现路径，具有重要的现实意义。本文将从一个技术实践者的视角，拆解此类全球健康影响评估研究可能涉及的数据处理、模型构建与结果分析流程，并提供一套可复现的数据分析框架示例。

1. 研究背景与核心概念解析

1.1 什么是超细颗粒物？

超细颗粒物通常指空气动力学直径小于或等于0.1微米的颗粒物。与更为人熟知的PM2.5相比，其粒径更小，数量浓度更高，表面积更大。由于尺寸极小，它们能够穿透人体肺泡屏障，直接进入血液循环系统，并可能抵达其他器官，因此其健康风险备受关注。在环境监测与研究中，超细颗粒物浓度常通过特殊仪器测量，数据获取和处理比常规PM2.5更为复杂。

1.2 全球疾病负担研究的方法论

“过早死亡”或“疾病负担”的归因分析，是环境流行病学的核心。其基本逻辑是：通过构建暴露-反应关系模型，估算在特定污染水平下，相较于一个理论上的最低风险水平，所额外导致的健康结局（如死亡、发病）数量。这类研究通常依赖于几类关键数据：

1. 全球暴露数据：来自卫星遥感反演、地面监测站网络和大气化学传输模型的融合数据产品。
2. 基线健康数据：全球或各国的人口、死亡率、疾病发病率数据，例如来自世界卫生组织或全球疾病负担研究。
3. 暴露-反应关系系数：来自长期队列研究或Meta分析的统计学参数，表示污染浓度每增加一个单位，特定健康风险增加的百分比。

1.3 技术挑战与价值

从技术实现角度看，完成这样一项全球评估面临多重挑战：多源异构数据的对齐与融合、高分辨率时空数据的处理、复杂统计模型的计算、结果的不确定性量化等。掌握相关的数据处理与分析技能，不仅是理解这类报告的前提，更是参与相关研究或开发环境健康预警系统的基础。

2. 环境准备与数据分析栈

为了模拟此类研究的核心分析步骤，我们需要搭建一个轻量化的数据分析环境。以下配置以Python生态为核心，适合进行数据探索、统计建模和可视化。

• 操作系统：Windows 10/11, macOS, 或 Linux (Ubuntu 20.04+) 均可。
• 编程语言：Python 3.8 或以上版本。
• 核心工具包：
  • pandas&numpy: 用于数据清洗、整理和数值计算。
  • geopandas&rasterio: 用于处理地理空间数据（如栅格格式的污染浓度图）。
  • xarray: 非常适合处理具有经纬度、时间维度的网格化科学数据。
  • statsmodels&scipy: 用于构建统计模型和进行假设检验。
  • matplotlib&seaborn: 用于数据可视化。
  • jupyter lab: 提供交互式分析环境，便于分步探索。
• 数据：我们将使用公开的模拟数据集进行演示，避免处理真实的巨量全球数据。

环境搭建命令： 建议使用conda创建独立环境以管理依赖。

# 创建并激活名为‘env_health’的conda环境 conda create -n env_health python=3.9 conda activate env_health # 安装核心数据分析库 conda install -c conda-forge pandas numpy matplotlib seaborn jupyterlab conda install -c conda-forge geopandas rasterio xarray pip install statsmodels

3. 核心分析流程拆解

一项完整的归因分析，在技术上可以简化为几个关键步骤。理解每一步的技术实现，比记住最终数字更重要。

3.1 数据获取与预处理

全球暴露数据通常是NetCDF或GeoTIFF格式的栅格数据，包含经纬度网格和每个格点的浓度值。健康基线数据则多为表格数据，需要与空间数据进行关联。

关键技术点：

• 空间对齐：将不同分辨率、不同投影的栅格数据重采样到统一网格。
• 人口加权：健康影响与受影响人口数量直接相关。需要将高分辨率人口分布数据与污染浓度数据叠加，计算人口加权平均暴露水平。
• 缺失值处理：对于监测数据缺失的区域，需要使用空间插值或模型数据填补。

3.2 暴露-反应关系模型的应用

这是归因计算的核心。通常采用对数线性关系模型（如Cox比例风险模型的近似）。归因分数（AF）的计算公式可简化为：

AF = (RR - 1) / RR其中，RR（相对风险）=exp(β * (C - C0))

• β: 暴露-反应关系系数（来自文献）。
• C: 实际暴露浓度。
• C0: 理论最低风险暴露水平。

为什么用这个模型？因为它能量化在特定暴露水平下，疾病风险相较于理想水平的超额部分，且在许多环境流行病学研究中被验证。

3.3 归因死亡数计算

将归因分数与基线死亡数结合：归因死亡数 = 基线死亡数 * AF

这一步需要在每个空间单元（如国家、网格）上分别计算，然后汇总到全球。

3.4 不确定性分析

任何模型结果都有不确定性。通常采用蒙特卡洛模拟方法，对关键参数（如β系数、基线死亡率）在其概率分布内进行多次随机抽样，重复整个计算过程，最终得到归因死亡数的置信区间。

4. 完整实战案例：模拟城市群PM2.5归因分析

我们以一个简化的模拟案例，演示如何为一个虚构的城市群计算PM2.5暴露导致的归因死亡数。本例聚焦于技术流程，数据均为模拟生成。

4.1 创建项目结构与模拟数据

首先，创建项目目录并初始化Jupyter Notebook或Python脚本。

# 文件：simulation_data.py import numpy as np import pandas as pd # 模拟生成5个城市的数据 np.random.seed(42) # 确保结果可复现 cities = ['City_A', 'City_B', 'City_C', 'City_D', 'City_E'] # 模拟数据：年均PM2.5浓度 (μg/m³), 人口(百万), 基线呼吸系统疾病死亡率(每10万人) sim_data = pd.DataFrame({ 'city': cities, 'pm25': np.random.uniform(20, 80, 5), # 浓度在20-80之间 'population': np.random.uniform(1, 10, 5), # 人口1-10百万 'baseline_mortality': np.random.uniform(50, 150, 5) # 基线死亡率 }) print("模拟城市数据：") print(sim_data)

4.2 定义核心计算函数

我们将归因计算的关键步骤封装成函数。

# 文件：attribution_calculation.py import numpy as np def calculate_attribution(pm25_concentration, baseline_deaths, beta, counterfactual=5.0): """ 计算单个区域的归因死亡数。 参数: pm25_concentration (float): PM2.5年均浓度 (μg/m³). baseline_deaths (float): 该疾病的基线死亡人数. beta (float): 暴露-反应关系系数，表示浓度每增加10μg/m³，相对风险的对数增加值. counterfactual (float): 理论最低风险浓度水平 (μg/m³). 常用5.0或2.4. 返回: tuple: (归因分数, 归因死亡数) """ # 计算相对风险 rr = np.exp(beta * (pm25_concentration - counterfactual) / 10.0) # 计算归因分数 af = (rr - 1) / rr if rr > 1 else 0.0 # 计算归因死亡数 attributable_deaths = baseline_deaths * af return af, attributable_deaths # 示例：使用一个来自文献的β系数（例如，针对心肺疾病死亡） # 假设β=0.156，表示PM2.5每增加10μg/m³，死亡风险增加约16.9% (exp(0.156)-1) BETA = 0.156 COUNTERFACTUAL = 5.0

4.3 应用计算并汇总结果

将计算函数应用到每个城市的数据上。

# 文件：main_analysis.py import pandas as pd from attribution_calculation import calculate_attribution, BETA, COUNTERFACTUAL from simulation_data import sim_data # 计算每个城市的基线死亡人数（基线死亡率 * 人口） sim_data['baseline_deaths'] = (sim_data['baseline_mortality'] * sim_data['population'] * 10) # 注意单位转换：每10万人 -> 实际人数 # 应用归因计算 results = [] for idx, row in sim_data.iterrows(): af, ad = calculate_attribution(row['pm25'], row['baseline_deaths'], BETA, COUNTERFACTUAL) results.append({ 'city': row['city'], 'pm25': row['pm25'], 'population_millions': row['population'], 'baseline_deaths': round(row['baseline_deaths'], 1), 'attributable_fraction': round(af, 4), 'attributable_deaths': round(ad, 1) }) results_df = pd.DataFrame(results) print("\n归因分析结果：") print(results_df.to_string(index=False)) # 汇总总归因死亡数 total_attributable_deaths = results_df['attributable_deaths'].sum() print(f"\n在该模拟场景下，这5个城市由PM2.5暴露导致的归因死亡数估算为：{total_attributable_deaths:.1f} 例")

4.4 结果可视化

使用matplotlib生成直观的图表。

# 文件：visualization.py import matplotlib.pyplot as plt import seaborn as sns sns.set_style("whitegrid") fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # 子图1：各城市PM2.5浓度与归因死亡数散点图 ax1 = axes[0] scatter = ax1.scatter(results_df['pm25'], results_df['attributable_deaths'], s=results_df['population_millions']*100, alpha=0.6, # 点大小代表人口 c=results_df['attributable_fraction'], cmap='Reds') ax1.set_xlabel('PM2.5 Concentration (μg/m³)') ax1.set_ylabel('Attributable Deaths') ax1.set_title('PM2.5 vs. Attributable Deaths (Bubble size=Population)') plt.colorbar(scatter, ax=ax1, label='Attributable Fraction') # 在点上标注城市名 for i, row in results_df.iterrows(): ax1.annotate(row['city'], (row['pm25'], row['attributable_deaths']), textcoords="offset points", xytext=(0,5), ha='center', fontsize=9) # 子图2：归因死亡数城市分布条形图 ax2 = axes[1] bars = ax2.bar(results_df['city'], results_df['attributable_deaths'], color='steelblue') ax2.set_xlabel('City') ax2.set_ylabel('Attributable Deaths') ax2.set_title('Distribution of Attributable Deaths by City') # 在柱子上添加数值标签 for bar in bars: height = bar.get_height() ax2.text(bar.get_x() + bar.get_width()/2., height + 0.5, f'{height:.1f}', ha='center', va='bottom', fontsize=10) plt.tight_layout() plt.savefig('attribution_analysis_results.png', dpi=300) plt.show()

4.5 运行与解读

运行上述脚本后，你会得到数据表格和两张图表。图表1显示了污染浓度、人口规模与归因死亡数的关系，通常可见浓度越高、人口越多，归因死亡数越高。图表2直观对比了各城市的归因负担。这个简化流程清晰地展示了从原始数据到健康影响评估结果的技术路径。

5. 常见问题与排查思路

在实际进行类似数据分析时，你可能会遇到以下问题：

6. 最佳实践与工程建议

将学术研究方法转化为稳健、可复现的分析流程，需要遵循以下工程实践：

1. 数据版本控制：使用DVC或Git LFS管理大型的原始栅格数据和中间处理结果。确保每次分析都能追溯到特定的数据版本。
2. 配置化参数管理：将所有关键参数（如β系数、理论最低风险浓度、疾病编码）放在独立的配置文件（如config.yaml或params.json）中，避免硬编码在脚本里。

   # config.yaml 示例 exposure_response: pm25_mortality: beta: 0.156 beta_se: 0.023 distribution: lognormal counterfactual: 5.0 diseases: - code: RES name: Respiratory Diseases baseline_file: data/baseline_respiratory.csv

3. 模块化代码设计：如示例所示，将数据读取、核心计算、可视化分离成不同模块或函数。这提高了代码的可读性、可测试性和复用性。
4. 不确定性量化是必须环节：任何点估计结果都必须附上不确定性范围（如95%置信区间）。使用概率分布描述参数不确定性，并采用蒙特卡洛模拟进行传播。
5. 敏感性分析：报告结果对关键假设的敏感性。例如，改变理论最低风险浓度从5.0到2.4 μg/m³，观察归因死亡数如何变化。这能增强结论的可靠性。
6. 文档与注释：在代码中详细注释数据来源、公式出处、单位换算过程。撰写README说明整个项目的运行环境、步骤和输出文件含义。
7. 可视化规范：地图可视化时，使用科学、客观的色带。避免使用可能误导读者的色带。在图中明确标注数据来源、处理方法和不确定性信息。

通过这个完整的从概念到代码的梳理，我们不仅理解了“超细颗粒物导致过早死亡”这一结论是如何从数据中产生的，更掌握了一套可以应用于类似环境健康影响评估项目的技术框架。这套方法的核心在于严谨的数据处理、清晰的模型实现和全面的不确定性考量。