《自然》：陆地蒸散量研究的统计失误被撤回-空间加权平均的计算方法

文章目录

前言
一、空间加权平均的计算方法
二、代码
- 1.Python 实现
- 2.MATLAB代码

前言

In this article, we calculated global land evapotranspiration for 2003 to 2019 using a mass-balance approach. To do this, we calculated evapotranspiration as the residual of the water balance, using an ensemble of datasets for precipitation, discharge and total water storage change. We made an error in calculating the global mean precipitation: we used arithmetic averaging to calculate the mean, instead of calculating a spatially weighted mean to account for the changing grid box size with latitude. As a result, the magnitudes of the global mean precipitation time series were underestimated. This impacted the subsequent calculation of global mean evapotranspiration, resulting in the mean evapotranspiration values being underestimated and altering some results. We are therefore retracting this article. We thank Ning Ma and others for bringing this error to our attention.
《自然》期刊于2021年5月26日发表了一篇题为“A 10 per cent increase in global land evapotranspiration from 2003 to 2019”的论文，研究显示2003年至2019年间全球陆地蒸散量增加了10% ± 2%。该研究曾为全球陆地蒸散量研究提供了重要数据和方向，但因统计错误于2022年2月24日被撤稿。问题在于作者在计算全球平均降水量时，错误地使用了算术平均值，而非考虑纬度变化的空间加权平均值。这导致全球平均降水时间序列被低估，进而影响了全球平均蒸散量的计算，使其被低估并改变了部分研究结果。这一错误削弱了研究结论的可靠性，误导了对全球陆地蒸散量变化趋势的判断，并影响了后续研究的方向。

一、空间加权平均的计算方法

在这里插入图片描述

二、代码

1.Python 实现

假设你有一个全球降水数据集 precip，其中：
纬度数组 lat（单位：度）。
降水数据 precip（形状为 [lat, lon]）。

import numpy as np# 生成示例纬度数组（假设等间距）
lat = np.linspace(90, -90, 180)  # 180 纬度点，示例为 1° 间隔
lon = np.linspace(-180, 180, 360)  # 360 经度点# 生成降水数据（示例数据）
precip = np.random.rand(len(lat), len(lon))  # 生成随机降水数据# 计算纬度权重
lat_radians = np.deg2rad(lat)  # 转换为弧度
weights = np.cos(lat_radians)  # 计算权重（cos(latitude)）# 进行加权平均（按纬度加权）
weighted_precip = np.average(precip, axis=0, weights=weights[:, np.newaxis])# 计算全球平均降水
global_mean_precip = np.mean(weighted_precip)print("全球加权平均降水量:", global_mean_precip)

2.MATLAB代码

% 生成示例纬度（1° 间隔）
lat = linspace(90, -90, 180);  
lon = linspace(-180, 180, 360);% 生成随机降水数据
precip = rand(length(lat), length(lon));% 计算纬度权重
weights = cosd(lat); % cos(latitude)，纬度单位为度
weights = weights / sum(weights); % 归一化权重% 计算加权平均
global_mean_precip = sum(sum(precip .* weights', 1)) / length(lon);disp(['全球加权平均降水量: ', num2str(global_mean_precip)]);

总结

为什么要加权？
由于球体表面的纬度网格面积不同，高纬度网格的面积小于低纬度网格，直接算术平均会导致高纬度地区被过度代表。
如何加权？每个网格的权重 wi=cos⁡(ϕi)wi=cos(ϕi)（纬度角的余弦）。
在计算全球平均时，用加权平均，而非算术平均。
这样可以得到更准确的全球平均降水值，避免文章中的计算错误！
在这里插入图片描述 # 生成示例纬度数组
import numpy as np
lat = np.linspace(90, -90, 180) # 180 纬度点
lon = np.linspace(-180, 180, 360) # 360 经度点
precip = np.random.rand(len(lat), len(lon)) # 随机降水数据
错误的算术平均
arithmetic_mean = np.mean(precip)
正确的空间加权平均
lat_radians = np.deg2rad(lat) # 转换为弧度
weights = np.cos(lat_radians) # 计算权重（纬度的 cos 值）
weights = weights / np.sum(weights) # 归一化
按纬度加权
weighted_mean = np.sum(precip * weights[:, np.newaxis]) / np.sum(weights)
print(f"错误的算术平均: {arithmetic_mean:.4f}“)
print(f"正确的加权平均: {weighted_mean:.4f}”)