后向轨迹可以帮助我们判断气团大致从哪里来，但当问题进一步变成“上游某个区域可能贡献了多少浓度”，几条轨迹线就不够了。这个时候需要用拉格朗日粒子扩散模型计算后向足迹，再把足迹与排放清单结合起来。
这一步看起来只是两个数组相乘，真正容易出问题的地方却集中在单位、网格面积、时间窗口和足迹归一化上。任何一处差一个数量级，最后的 PSC（Potential Source Contribution） 都可能看起来合理，却无法解释。
Summary
潜在源贡献计算可以写成“后向足迹 × 网格排放速率”。计算前要先确认足迹已经按释放质量归一化，并统一到
h/m³；排放则统一到μg/h/grid。两者相乘得到各网格对受体的潜在浓度贡献，再对空间和时间求和得到总贡献。
LPDM 是什么
LPDM 是 Lagrangian Particle Dispersion Model 的缩写，也就是拉格朗日粒子扩散模型。
这类模型会在气象场中释放大量计算粒子，再根据平均风场、湍流和垂直运动等过程更新粒子的位置。单条后向轨迹给出的是一条理想化空气质点路径；LPDM 同时跟踪大量粒子，可以描述气团在输送过程中的扩散，并把结果统计到规则网格上。最后得到的通常是一片连续分布，可以进一步与网格化排放清单对应。
HYSPLIT 是常用的大气输送与扩散模型之一。NOAA Air Resources Laboratory 对它的介绍中提到，HYSPLIT 的计算结合了随空气质点移动的拉格朗日方法和固定浓度网格的欧拉方法，可以进行轨迹、输送、扩散、化学转化与沉降模拟。我们使用的 LPDM 程序基于 HYSPLIT 发展而来，所以下面的 PSC 计算会以这套后向足迹输出为例，但同样的单位和匹配逻辑也适用于其他能够输出源—受体敏感度的拉格朗日模型。
后向足迹是怎样得到的
PSC 使用的是后向扩散结果。普通后向轨迹只有路径端点，不能直接与排放清单逐格相乘。获得后向足迹时，大致需要经过下面几个环节：
- 把观测站点或关注位置设为受体，同时确定受体时间、高度和向后追踪时长；
- 从受体位置释放一组携带已知示踪质量的计算粒子，在气象场驱动下反向积分；
- 统计粒子在上游各网格和高度层中的停留或浓度敏感度，形成随经纬度分布的后向足迹；
- 将模型原始输出转换成规则网格数据，并按释放质量归一化，得到后续可以与排放速率相乘的足迹。
实际计算时还要选择输出网格、垂直层和累计时间窗。例如同一次后向模拟可能分别输出 6 小时、12 小时、1 天或更长时间的累计足迹。这里先得到的是受体对上游排放的敏感度，真实排放强度要到 PSC 计算中才会加入。
Important
PSC 的输入应当是后向扩散产生的网格化足迹，并且已经按释放质量归一化。普通轨迹端点图、轨迹频率图和未归一化的示踪质量不能直接代替这个结果。
先把 PSC 想清楚
这里所说的 PSC 是 Potential Source Contribution，也就是潜在源贡献。它与轨迹统计里常见的 PSCF（Potential Source Contribution Function）不是同一个量。PSCF 通常根据轨迹端点和高浓度样本出现的频率判断潜在源区；这里的 PSC 则把 LPDM 后向足迹与排放清单结合起来，结果可以落到浓度单位。
后向足迹可以理解为源—受体敏感度：某个上游网格的值越高，说明该网格排放对受体浓度的影响潜力越大。NOAA Air Resources Laboratory 对 HYSPLIT footprints 的说明也指出，这类场用于估计观测浓度对上风向地表通量的敏感度。HYSPLIT footprints 官方说明
如果把网格记为 i,j，时间记为 t，最基本的计算关系是：
PSC_{i,j}=\sum_t F_{t,i,j}\,E_{t,i,j}C_{receptor}=\sum_i\sum_j PSC_{i,j}其中：
F是归一化后的 LPDM 足迹，常用单位为h\,m^{-3}；E是每个网格的排放速率，单位为\mu g\,h^{-1}\,grid^{-1}；PSC_{i,j}是单个网格对受体的潜在浓度贡献，单位为\mu g\,m^{-3}；C_{receptor}是对全部源区求和后的受体总贡献。
Important
有些后处理文件把足迹单位写成
mass·h/m³。只有当模拟释放的是单位质量，或者结果已经除以释放质量时，才能把它作为h/m³与排放速率相乘。如果文件里仍保留实际释放的示踪质量，需要先除以释放质量。不能只看数值大小判断。
HYSPLIT 的前向浓度计算允许使用任意一致的质量单位，官方示例中输入 kg/h 会得到 kg/m³ 输出。排放与输出单位说明 后向计算虽然换成了源—受体敏感度，单位闭合仍然要遵守同一逻辑。
把两类排放清单统一到网格排放速率
PSC 计算最稳妥的中间单位是 μg/h/grid。无论原始清单给的是月总量还是面源通量，都先转换到“每个网格、每小时排放多少微克”，再与足迹相乘。
月总排放量：Mg/month
如果清单已经给出每个网格的月总排放量，单位为 Mg/month/grid，换算关系为：
E_{\mu g/h}=E_{Mg/month}\times\frac{10^{12}}{D\times24}D 是当月实际天数。以 30 天为例：
E_{\mu g/h}=E_{Mg/month}\times\frac{10^{12}}{720}这里的数量级很容易写错。Mg 是 megagram，也就是 10^6 g，因此 1\,Mg=10^{12}\,\mu g。如果清单实际单位是 kg/month，换算因子是 10^9。计算前最好直接查看 NetCDF 变量的 units 属性和清单说明。
月排放量已经是单格总量时，不需要再乘网格面积。还要使用对应月份的真实天数，2 月、大小月和闰年都不应固定除以 720 小时。
面源排放速率：kg/m²/s
如果清单给的是单位面积、单位时间的排放通量 q，需要先乘网格面积：
E_{\mu g/h/grid}=q_{kg/m^2/s}\times A_{grid}\times10^9\times3600也可以写成：
E_{\mu g/h/grid}=q_{kg/m^2/s}\times A_{grid}\times3.6\times10^{12}对于规则经纬度网格，不能在整个区域都使用同一个面积。经度方向的实际距离会随纬度缩短。设地球半径为 R，网格经纬度跨度分别为 \Delta\lambda 和 \Delta\phi，中心纬度为 \phi，更合适的球面面积公式是：
A=R^2\Delta\lambda\left[\sin\left(\phi+\frac{\Delta\phi}{2}\right)-\sin\left(\phi-\frac{\Delta\phi}{2}\right)\right]公式中的角度都要先转成弧度。以 0.1° × 0.1° 网格为例，赤道附近面积约为 1.24\times10^8\,m^2，30°N 附近约为 1.07\times10^8\,m^2，45°N 附近只剩约 8.74\times10^7\,m^2。研究区跨度较大时，固定使用 1.23\times10^8\,m^2 会形成明显的纬向偏差。
Note
“总量”和“通量”必须分开处理。
Mg/month/grid属于单格总量，只需要做质量与时间换算；kg/m²/s属于面源通量，必须乘网格面积。两者混用会直接多乘或少乘一次面积。
从 LPDM 输出到 PSC
先还原足迹数值
在这套基于 HYSPLIT 的 LPDM 后处理文件中，足迹常以 10 的幂指数保存。例如 -13.6687 表示：
F=10^{-13.6687}这里要单独处理文件里的零值。如果 0.0000 代表“该网格没有足迹”，它必须继续保留为 0，不能执行 10^0 后变成 1。还应确认后处理输出的缺测值、填充值和真实指数 0 是否有明确区分。
同一个文件可能同时给出 6 小时、12 小时、1 天和 7 天等累计足迹。这些列表示不同的回溯时间窗，彼此包含，不能把它们相加。研究使用 72 小时后向模拟，就应选择与 72 小时定义一致的足迹结果，而不是把多个累计列拼在一起。
再做时空对齐
两个数据都是 0.1° 分辨率，并不代表格点已经对应。至少要检查：
- 经纬度表示的是格点中心还是边界；
- 纬度是从南到北还是从北到南；
- 经度采用
0–360°还是-180–180°； - 两个网格的起点、终点和格点数量是否完全一致；
- LPDM 时间使用 UTC，排放清单使用 UTC、本地时还是月平均；
- 重网格前的排放变量是总量还是通量。
如果网格不同，排放总量适合使用守恒重网格，保证区域总排放在转换前后基本一致。直接用普通插值处理单格总量，可能在改变分辨率时顺带改变总排放。
时空一致后，逐格计算：
PSC_{i,j}=F_{i,j}\times E_{i,j}得到的二维场可以直接画潜在贡献分布。对研究区域求和后，才是该排放清单在给定假设下对受体的总浓度贡献。
质量浓度转为体积混合比
气态污染物有时需要把 \mu g/m^3 转成 ppbv。在给定温度和气压下：
C_{ppbv}=C_{\mu g/m^3}\times\frac{10^3RT}{M p}其中 M 是物种摩尔质量，单位为 g/mol；T 为 K；p 为 Pa。若使用摩尔体积 V_m，公式可写成：
C_{ppbv}=C_{\mu g/m^3}\times\frac{V_m\,(L/mol)}{M\,(g/mol)}在 25 ℃、1 atm 下，V_m 约为 24.45 L/mol；在 0 ℃、1 atm 下约为 22.41 L/mol。摩尔体积应放在分子，摩尔质量放在分母。若要与观测严格比较，最好使用观测所采用的参考温压，或者直接带入实际 T、p，并在方法中写清楚。
这个换算只适合有明确摩尔质量的气体。颗粒物质量浓度通常保留为 \mu g/m^3；NOx、NMVOC 等汇总物种还要先确认清单采用的质量基准，不能随意指定一个摩尔质量。
Python 计算与结果检查
一段可复用的计算框架
下面的代码只保留计算骨架。变量名、维度名和缺测值需要按自己的文件调整。
import calendar
import numpy as np
import xarray as xr
R_EARTH = 6_371_000.0 # m
def regular_latlon_grid_area(lat, dlat=0.1, dlon=0.1):
"""返回规则经纬网格中每个纬度带的单格面积，单位 m2。"""
lat = np.asarray(lat)
phi = np.deg2rad(lat)
dphi = np.deg2rad(dlat)
dlambda = np.deg2rad(dlon)
return (
R_EARTH**2
* dlambda
* (np.sin(phi + dphi / 2) - np.sin(phi - dphi / 2))
)
def restore_footprint(exponent):
"""按当前后处理约定，把 0 保留为无足迹，其余值还原为 10 的指数。"""
return xr.where(exponent == 0.0, 0.0, 10.0**exponent)
def monthly_total_to_rate(emission_Mg_month, year, month):
"""Mg/month/grid 转为 μg/h/grid。"""
days = calendar.monthrange(year, month)[1]
return emission_Mg_month * 1.0e12 / (days * 24)
def surface_flux_to_rate(emission_flux, dlat=0.1, dlon=0.1):
"""kg/m2/s 转为 μg/h/grid。"""
area_1d = xr.DataArray(
regular_latlon_grid_area(emission_flux.lat, dlat, dlon),
coords={"lat": emission_flux.lat},
dims=("lat",),
)
return emission_flux * area_1d * 1.0e9 * 3600.0
# 选择与研究回溯时长一致的一列，例如 72 h 对应的累计足迹
footprint = restore_footprint(lpdm_exponent)
# 根据清单单位二选一，另一行保持注释
emission_ug_h = monthly_total_to_rate(emission_Mg_month, year, month)
# emission_ug_h = surface_flux_to_rate(emission_flux, dlat=0.1, dlon=0.1)
# 网格若不能严格对齐，这里会直接报错，避免静默错位
footprint, emission_ug_h = xr.align(
footprint,
emission_ug_h,
join="exact",
)
psc_grid = footprint * emission_ug_h
psc_total = psc_grid.sum(dim=("lat", "lon"), skipna=True)xr.align(…, join="exact") 适合在网格本应完全一致时做最后一道检查。如果两套数据原本就不是同一网格，应先明确重网格方法，再进入这一步。
结果检查与适用边界
一个示例结果：
PSC 图画出来以后，先不要急着解释高值区。建议依次检查下面几项。
- 足迹是否已经除以释放质量，最终可按
h/m³使用 - 指数型输出是否已还原，零值有没有被错误地变成 1
- 月排放量中的
Mg是否按10^{12}\,\mu g换算 - 面源通量是否乘了逐纬度变化的网格面积
- 足迹和排放的经纬度、顺序、时间与物种是否一致
- 选择的是单个累计时间窗，没有把 6 h、12 h、1 d、7 d 多列相加
- 逐格 PSC 求和后，结果单位是否仍能闭合为
\mu g/m^3 - 转换为
ppbv时，温压条件和摩尔质量是否写清楚
还要注意，PSC 中的“潜在”两个字不能省略。这个计算依赖几个重要假设：输送响应近似线性，排放清单能够代表回溯时段，足迹与排放在高度和时间上匹配，目标物种在输送过程中没有超出模型设置的化学损失或生成。
因此，惰性或近似惰性的示踪物通常更容易解释。若直接把 NOx 排放与受体 O3 浓度相乘，或者忽略 SO₂ 氧化、颗粒物二次生成与沉降，得到的数值不能当成真实浓度贡献。高架点源也不能简单塞进地面排放层，需要做垂直分配，并使用与释放高度对应的足迹。
Warning
PSC 适合回答“在当前输送场和排放清单下，哪些区域具有较高贡献潜力”。它不能自动修正排放清单偏差，也不会替代化学输送模式。与观测对比时，还要考虑背景浓度、清单遗漏、气象误差、粒子统计噪声以及化学和沉降过程。
总结
从 LPDM 足迹到 PSC，核心流程只有四步：确认足迹归一化，统一排放为单格每小时排放量，完成时空对齐，逐格相乘后再按研究问题求和。
公式不复杂，真正决定结果是否可信的是每一步的单位与边界。尤其是 Mg 的数量级、经纬网格面积、指数输出中的零值，以及 ppbv 换算方向，最好在脚本里显式写出并保留中间结果。这样后面更换月份、物种或排放清单时，整条计算链仍然能复查。