计算整层水汽通量散度是先积分后散度还是先散度后积分？

admin管理员组

文章数量:1444586

计算整层水汽通量散度是先积分后散度还是先散度后积分？

先积分后散度还是先散度后积分？

温馨提示

由于可视化代码过长隐藏，可点击运行Fork查看 若没有成功加载可视化图，点击运行可以查看 ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可

前言

在气象学中，整层水汽通量散度（Integrated Moisture Flux Divergence, IMFD）是评估大气柱内水分输送及变化的重要指标

实际上我对于计算整层水汽通量散度是先积分后散度还是先散度后积分有困惑，

究竟两者有什么差异，还是说只是单纯等价，一个快一个慢呢

ps :下面结果基于era5数据，不等于所有数据适用

项目目标

本项目的主要目标包括：

1. 计算整层水汽通量散度：使用ERA5再分析数据，结合metpy、xarray等库，计算从地面至一定高度（例如1000 hPa至300 hPa）之间的整层水汽通量散度。
2. 对比先积分后散度与先散度后积分差异

---

1. 先积分后散度

代码语言：javascript代码运行次数：0运行复制

%%time
import xarray as xr
from metpy.units import units
import numpy as np
import metpy.calc as mpcalc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# 读取 ERA5 数据
ds = xr.open_dataset('/home/mw/input/era58091/ERA5-2023-08_pl.nc')
levels = [1000, 925, 850, 700, 600, 500, 400, 300]  # 定义需要整合的高度层次
times = ds['time']  # 获取时间维度

# 提取所需变量
q = ds['q'].sel(level=levels, time=times[0]) * units('kg/kg')  # 比湿
u = ds['u'].sel(level=levels, time=times[0]) * units('m/s')    # 经向风速
v = ds['v'].sel(level=levels, time=times[0]) * units('m/s')    # 纬向风速

# 计算各层的水汽通量
qx = q * u / 9.81# 经向水汽通量
qy = q * v / 9.81# 纬向水汽通量

# 进行垂直积分获得整层水汽通量
integrated_qx = -qx.integrate(coord='level')
integrated_qy = -qy.integrate(coord='level')

# 计算网格间距
dx, dy = mpcalc.lat_lon_grid_deltas(ds['longitude'], ds['latitude'])

# 计算散度
div_integrated_q = mpcalc.divergence(integrated_qx, integrated_qy, dx=dx, dy=dy)
div_integrated_q = div_integrated_q.rename('integrated_moisture_flux_divergence')*10

def plot_integrated_moisture_flux_divergence(div_q, extent=None):
    fig = plt.figure(figsize=(10, 8), dpi=300)
    ax = plt.axes(projection=ccrs.PlateCarree())
    
    # 如果提供了extent参数，则设置地图范围
    if extent isnotNone:
        ax.set_extent(extent, crs=ccrs.PlateCarree())

    im = ax.pcolormesh(div_q.longitude, div_q.latitude, div_q * 1e5, 
                       cmap='RdBu', vmin=-30, vmax=30, transform=ccrs.PlateCarree())
    
    # 添加色标
    cbar = fig.colorbar(im, ax=ax, orientation='horizontal', pad=0.05)
    cbar.set_label('Divergence 10-5 *(g•(cm²•s)-1)', fontsize=12)
    
    # 添加地图特征
    ax.add_feature(cfeature.COASTLINE)
    ax.set_xlabel('Longitude', fontsize=12)
    ax.set_ylabel('Latitude', fontsize=12)
    ax.set_title('Integrated Moisture Flux Divergence', fontsize=16)

    plt.show()

# 定义需要显示的地图范围 [最小经度, 最大经度, 最小纬度, 最大纬度]
extent = [90, 140, 15, 55]  # 例如：东亚地区

# 调用绘图函数，并传入extent参数
plot_integrated_moisture_flux_divergence(div_integrated_q, extent=extent)

代码语言：javascript代码运行次数：0运行复制

CPU times: user 1.25 s, sys: 149 ms, total: 1.4 s
Wall time: 1.31 s

2. 先散度后积分

代码语言：javascript代码运行次数：0运行复制

%%time
import xarray as xr
from metpy.units import units
import numpy as np
import metpy.calc as mpcalc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

def plot_integrated_moisture_flux_divergence(div_q, extent=None):
    fig = plt.figure(figsize=(10, 8), dpi=300)
    ax = plt.axes(projection=ccrs.PlateCarree())
    
    # 如果提供了extent参数，则设置地图范围
    if extent isnotNone:
        ax.set_extent(extent, crs=ccrs.PlateCarree())

    im = ax.pcolormesh(div_q.longitude, div_q.latitude, div_q * 1e5, 
                       cmap='RdBu', vmin=-30, vmax=30, transform=ccrs.PlateCarree())
    
    # 添加色标
    cbar = fig.colorbar(im, ax=ax, orientation='horizontal', pad=0.05)
    cbar.set_label('Divergence 10-5 *(g•(cm²•s)-1)', fontsize=12)
    
    # 添加地图特征
    ax.add_feature(cfeature.COASTLINE)
    ax.set_xlabel('Longitude', fontsize=12)
    ax.set_ylabel('Latitude', fontsize=12)
    ax.set_title('Integrated Moisture Flux Divergence', fontsize=16)

    plt.show()

# 读取 ERA5 数据
ds = xr.open_dataset('/home/mw/input/era58091/ERA5-2023-08_pl.nc')
levels = [1000, 925, 850, 700, 600, 500, 400, 300]  # 定义需要整合的高度层次
times = ds['time']  # 获取时间维度

# 提取所需变量并确保level顺序正确（低层→高层）
ds_sorted = ds.sortby('level', ascending=False)
q = ds_sorted['q'].sel(level=levels, time=times[0]) * units('kg/kg')  # 比湿
u = ds_sorted['u'].sel(level=levels, time=times[0]) * units('m/s')    # 经向风速
v = ds_sorted['v'].sel(level=levels, time=times[0]) * units('m/s')    # 纬向风速

# 计算网格间距（仅在水平方向计算一次）
dx, dy = mpcalc.lat_lon_grid_deltas(ds_sorted['longitude'], ds_sorted['latitude'])

# 逐层计算水汽通量散度
divergence_per_level = []
for lev in q.level.values:
    q_lev = q.sel(level=lev)
    u_lev = u.sel(level=lev)
    v_lev = v.sel(level=lev)
    
    qx_lev = q_lev * u_lev / 9.81# 经向水汽通量
    qy_lev = q_lev * v_lev / 9.81# 纬向水汽通量
    
    # 计算单层散度
    div_lev = mpcalc.divergence(qx_lev, qy_lev, dx=dx, dy=dy)
    divergence_per_level.append(div_lev)

# 合并各层散度并垂直积分
divergence_da = xr.concat(divergence_per_level, dim='level')
integrated_divergence = -divergence_da.integrate(coord='level') * 10# 符号修正与单位调整
integrated_divergence = integrated_divergence.rename('integrated_moisture_flux_divergence')

# 调用绘图函数
plot_integrated_moisture_flux_divergence(integrated_divergence, extent=[90, 140, 15, 55])

代码语言：javascript代码运行次数：0运行复制

CPU times: user 1.69 s, sys: 155 ms, total: 1.84 s
Wall time: 1.75 s

两种方法的对比分析

特征	先积分通量 → 计算散度	先计算散度 → 垂直积分
数学等价性	理论上等价（积分与微分交换顺序）	理论上等价
数值误差来源	对平滑后的通量场计算散度，可能放大噪声	对每层散度独立计算，噪声被积分平滑
层级顺序敏感性	高度敏感（错误顺序直接导致符号错误）	较不敏感（逐层计算可独立验证符号）
台风区域合理性	易因层级顺序错误导致辐散（正值）异常	更可靠（可逐层检查低层辐合）
计算效率	较快（仅需一次积分和一次散度计算）	较慢（需逐层计算散度）
推荐场景	数据层级顺序已验证正确时	调试或层级顺序不确定时

---

关键差异解释

层级顺序敏感性：

• 先积分通量时，若 level 未按气压递减排序（如默认升序），积分方向会反向，导致散度符号错误。
• 先计算散度时，即使 level 顺序错误，只需在垂直积分时调整符号即可（通过 -divergence_da.integrate()）。

物理合理性验证：

• 在台风低层（如 850 hPa），单层散度应表现为 辐合（负值）。先计算散度的方法允许逐层检查：# 检查 850 hPa 散度符号 div_850 = divergence_da.sel(level=850) plot_integrated_moisture_flux_divergence(div_850, extent=[90, 140, 15, 55])

---

结论

• 优先使用先计算散度 → 垂直积分，尤其在调试阶段或数据层级顺序不确定时。

当你对于数据结构不熟悉时，最好是单独计算一层先测试，批量是往后的事情

从结果而言era5数据无所谓先后，两者皆可，先积分后散度因为没有套for循速度更快罢了

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-03-23，如有侵权请联系 cloudcommunity@tencent 删除变量地图可视化数据调试

本文标签： 计算整层水汽通量散度是先积分后散度还是先散度后积分