做气象数据处理的人应该都有过这种经历拿到一套全球格点数据比如 CMIP6 的月平均降水或者 ERA5 的逐时温度第一反应就是mean(mean(data, 1), 2)然后画一条区域平均的时间序列。但画出来之后发现曲线的数值和论文里总是对不上偏大偏小都有有时候趋势都差出一截。我之前也困惑过很久直到把面积加权这件事彻底搞明白才意识到问题多半出在这里。Matlab 处理气象数据系列写到第六篇这篇专门把“面积加权”聊透。它不是锦上添花的技巧而是处理格点数据时绕不开的一步尤其是涉及全球平均、半球平均、大范围区域平均甚至 EOF 分析、模式评估这类统计量时不做面积加权结果就是错的。这篇会从原理讲到代码再把我实际踩过的坑一并列出来希望能帮你省下几天的折腾时间。1. 为什么气候数据的区域平均不能直接取平均1.1 网格面积不等的“陷阱”到底藏在哪很多人以为全球等经纬度网格上每个格点代表的面积是一样的这是个非常普遍的误解。你可以想象把一个地球仪表面用纬线和经线切成小格子纬线越往两极纬线圈越短所以相同经度间隔对应的实际地面距离就越小。也就是说在赤道附近一个 1° × 1° 的格网代表的地面面积可能有 12300 多平方公里而到了 60°N同样 1° × 1° 的格网面积只剩下一半左右。如果直接对所有格点求算术平均就等于把高纬度那些面积偏小的格点和赤道附近面积很大的格点一视同仁高纬地区的数值会被不成比例地放大。举个例子全球 2.5° × 2.5° 的网格上60°N 以北的格点数量占全球格点数的约六分之一但如果按实际面积算这片区域只占全球总面积的 6.7% 左右。算术平均相当于把高纬区域的实际贡献放大了差不多两倍半。这个偏差可不是误差范围内的小事在讨论全球平均温度、降水总量这类问题时足以改变结论。1.2 面积加权到底在加什么面积加权的思想很简单每个格点的数值乘上它所代表的实际面积然后求和再除以总面积。用公式表达就是[ \bar{X} \frac{\sum_{i1}^{n} X_i \cdot A_i}{\sum_{i1}^{n} A_i} ]其中 (X_i) 是第 i 个格点的数值(A_i) 是这个格点代表的实际面积。这个公式灵魂在于 (A_i)在等经纬度网格上(A_i) 近似正比于该格点纬度的余弦值。推导起来也不难球面上一个格元的面积是 (dA R^2 \cos\phi , d\phi , d\lambda)其中 (\phi) 是纬度(\lambda) 是经度。当经纬度间隔固定时后面那一串 (R^2 d\phi d\lambda) 是常数所以唯一的变量就是 (\cos\phi)。所以在实际计算里最常用的权重就是cosd(lat)。注意这里要强调一下纬度单位如果是度一定要用cosd而不是coscos默认接收弧度你直接用会把整个权重矩阵算错而且这种错很难肉眼发现。1.3 什么时候可以偷懒不做面积加权不是所有情况都需要面积加权过度设计也没必要。我自己判断的标准是看研究区域的纬度跨度和网格分辨率。如果区域很小比如单个省、某个流域纬度跨度只有两三度网格间面积差异不大算术平均的误差很小可以不做。如果是高分辨率网格且区域不大比如 0.25° 数据覆盖一个城市群面积差异同样可以忽略。但如果处理的是全球、半球、或者像 30°N-60°N 这种纬度跨度大的带状区域面积加权就是必须的了。另外特别提一句有些再分析数据或模式输出用的是高斯格点Gaussian grid纬度方向不是等间距的这时候cosd(lat)只能算近似更严谨的做法是按实际格点面积权重来算后面会讲方法。2. 面积加权的三种 Matlab 实现方案2.1 方案一cos(纬度) 快速权重覆盖九成场景对于绝大多数等经纬度网格数据用cosd(lat)构造权重矩阵就足够了。核心代码很短几行就能搞定lat ncread(file, lat); % 纬度向量单位度 lon ncread(file, lon); % 经度向量单位度 % 构造权重矩阵每一行的权重相同因为同一纬度带的面积一致 wgt repmat(cosd(lat(:)), 1, length(lon)); wgt wgt ./ sum(wgt(:)); % 归一化所有权重加起来等于1之后做全球平均就是一句矩阵乘法data_2d reshape(data, size(data,1)*size(data,2), size(data,3)); % 假如 data 是 [lat x lon x time] global_avg wgt(:) * data_2d; % 结果是 1 x time这个方法的核心就是那个repmat它把纬度权重复制到每个经度上形成一个和网格一样大的权重矩阵。归一化不能省否则加权求和的结果带上了面积量纲不是一个“平均”的概念了。2.2 方案二逐格点实际面积权重适合高斯网格或不规则网格cosd(lat)的近似默认了经度间隔和纬度间隔都是恒定的并且网格近似矩形。但高斯网格、或者某些区域模式输出的旋转网格纬向间隔并不规则这时候更可靠的方案是直接计算每个格点的实际面积。计算公式是[ A_{i,j} R^2 \cdot \cos\phi_j \cdot \Delta\phi_j \cdot \Delta\lambda_i ]在 Matlab 里可以写成R 6371000; % 地球平均半径单位米 dphi abs(gradient(lat)) * pi/180; % 纬度方向间隔转弧度 dlambda abs(gradient(lon)) * pi/180; % 经度方向间隔转弧度 [LON, LAT] meshgrid(lon, lat); A R^2 * cosd(LAT) .* dphi .* dlambda; % 每个格点的实际面积gradient是这里比较关键的函数它会根据相邻格点自动计算每个格点对应的间隔内部格点用中心差分边界用单边差分。对规则网格来说结果和直接用常数值一样对高斯网格来说它能体现纬向间隔的变化。实际业务里我通常用这个方案做验证确保方案一的结果没跑偏。2.3 方案三掩膜配合内置函数处理区域统计更省心如果目标是计算某个不规则区域的面积加权平均比如亚洲大陆、某个气候区就要在权重之外叠加一个掩膜矩阵 maskmask 大小和网格一致区域内为 1区域外为 0。% 假设 mask 是逻辑矩阵 [lat x lon] wgt repmat(cosd(lat(:)), 1, length(lon)); wgt(~mask) 0; % 只保留目标区域内的权重 wgt wgt ./ sum(wgt(:)); % 注意只在区域内归一化这里有个细节归一化一定要在置零 mask 之后再做。如果你先全局归一化再乘 mask区域内的权重之和就不等于 1 了算出来的平均值会整体偏小。当初我第一次写就是顺序搞反结果区域平均结果比该有的值低了不少排查了很久才意识到这个排序问题。三种方案的适用场景和复杂度我列成表格方便对照方案适用数据精度代码量备注cosd 权重规则等经纬度网格高极少最常用推荐优先尝试逐格点面积高斯网格、旋转网格最高中等用 gradient 算间隔掩膜加权任意网格 区域统计取决于权重方案中等与方案一或方案二结合3. 从 NetCDF 数据到加权区域平均的完整实操3.1 读取数据并确认维度顺序不管数据来源是 CMIP6、ERA5 还是站点插值产品第一步永远是确认维度顺序。NetCDF 文件里变量最常见的存储顺序是[lon x lat x time]而 Matlab 的矩阵下标习惯是[行 x 列]也就是第一维是纬度、第二维是经度。这一反一正特别容易踩坑。我用的常规流程是这样file pr_Amon_GFDL-ESM4_historical_r1i1p1f1_gn_195001-201412.nc; lat ncread(file, lat); lon ncread(file, lon); pr ncread(file, pr); % 读出来通常是 [lon x lat x time] % 用 permute 把维度调成 [lat x lon x time] pr permute(pr, [2 1 3]);permute(pr, [2 1 3])的意思是把前两维交换第三维时间不变。调完之后size(pr,1)对应 latsize(pr,2)对应 lon这样后续用repmat构造权重矩阵时就不容易对错方向。建议读完之后先做一次 sanity checkdisp([length(lat), length(lon), size(pr,3)])确认维度和预期一致再往下走。这一步多花十秒钟能省掉后面一大半的 debug 时间。3.2 生成区域掩膜经纬度范围、陆地和自定义形状掩膜 mask 的生成方式取决于你想要什么区域。最简单的是经纬度范围式比如我要算热带印度洋-西太平洋暖池区域20°S-20°N100°E-150°E的海表温度平均[LON, LAT] meshgrid(lon, lat); mask (LAT -20 LAT 20) (LON 100 LON 150);更复杂一点要陆地平均就得加载陆地掩膜数据。Matlab 的 Mapping Toolbox 里有landmask函数不过它基于 GSHHS 海岸线数据分辨率有限。实测下来处理全球尺度问题够用处理区域尺度容易在沿海格点出错。更稳妥的做法是读取陆地和海洋分率的文件比如 ERA5 自带的lsm变量用lsm 0.5作为陆地判定。如果你要严格按国界或流域边界那就得用 shapfile 了。基本思路是用shaperead读入边界多边形然后用inpolygon判断每个格点中心是否在多边形内。这个过程在网格很大的时候跑得比较慢但逻辑清晰不容易出错。下面给个示例骨架shape shaperead(some_basin.shp); [LON, LAT] meshgrid(lon, lat); mask inpolygon(LON, LAT, shape.X, shape.Y);3.3 核心函数加权区域平均的完整代码我把常用的逻辑封装成了一个自认为比较健壮的函数放在我的气象数据处理工具集里。它支持时间序列、NaN 自动处理、掩膜可选function avg area_weighted_mean(data, lat, lon, mask) % area_weighted_mean 计算格点数据的面积加权区域平均 % 输入 % data : [lat x lon] 或 [lat x lon x time]单位为任意 % lat : 纬度向量单位度升序或降序均可 % lon : 经度向量单位度 % mask : 可选逻辑矩阵 [lat x lon]目标区域为 true % 输出 % avg : 标量或 1 x time 的区域加权平均序列 % 统一维度到三维 if ndims(data) 2 data reshape(data, size(data,1), size(data,2), 1); end [nlat, nlon, nt] size(data); % 默认 mask 为全区域 if nargin 4 || isempty(mask) mask true(nlat, nlon); end % 构造权重矩阵仅在 mask 内保留权重 wgt repmat(cosd(lat(:)), 1, nlon); wgt(~mask) 0; wgt wgt ./ sum(wgt(:)); % 展平数据逐时间步做加权平均 data_2d reshape(data, nlat*nlon, nt); wgt_vec wgt(:); avg zeros(1, nt); for t 1:nt col data_2d(:, t); valid isfinite(col) (wgt_vec 0); if sum(valid) 0 avg(t) NaN; else w wgt_vec; w(~valid) 0; % 数据为 NaN 的格点权重也清零 avg(t) sum(col(valid) .* w(valid)) / sum(w); end end if nt 1 avg avg(1); end end这里valid isfinite(col) (wgt_vec 0)是我后来加上去的。最初版本只判断isfinite遇到 mask 外权重本来就是 0 的格点也没事但后来发现如果数据里有 NaN而权重不为 0那么sum(col .* w)会因为 NaN 而污染整个结果。所以必须在每个时间步做一次有效点筛选并且把无效格点的权重同步清零这就保证了 Nan 不会拖垮整条序列。3.4 结果验证退化测试是必须做的一步写完加权平均功能先别急着拿到业务数据上跑做一次退化检测能帮你确认代码里有没有低级错误。具体做法是构造一个全为 1 的场面积加权平均结果必须等于 1构造一个纬度值本身作为场全球面积加权平均应该在 0 附近因为纬度分布关于赤道对称同时算术平均会略偏正或偏负具体偏哪边取决于纬度网格轴的取点。% 退化测试 1常数场 test1 ones(length(lat), length(lon), 1); result1 area_weighted_mean(test1, lat, lon); assert(abs(result1 - 1) 1e-12, 退化测试失败常数场未通过); % 退化测试 2线性纬度场 [LAT, ~] meshgrid(lat, lon); test2 permute(LAT, [2 1 3]); % 调整成 [lat x lon x 1] result2 area_weighted_mean(test2, lat, lon); fprintf(纬度场面积加权平均 %.4f°\n, result2);如果第二项结果明显偏离 0比如跑到 10 以上那基本可以断定权重矩阵的方向和网格维度没对齐。以前我有个同事用repmat构造权重的时候把纬度向量按行复制结果权重矩阵是[lon x lat]方向和数据的[lat x lon]完全反了算出来的全球平均温度长期比 ERA5 报告值高 8°C就是这种细小的方向错误导致的。4. 实战中遇到的坑与排查方法4.1 常见问题速查表我把这几年处理气象数据时遇到的高频问题整理成一张速查表覆盖了面积加权的大部分雷区现象可能原因解决办法加权平均结果比文献值差很多权重公式用错cos和cosd混用检查纬度单位统一使用cosd区域平均结果整体偏小先归一化后乘 mask先置零 mask 外权重再归一化时间序列在某个时段突然变 NaN该时段存在 NaN 数据在计算中动态屏蔽 NaN同步权重置零加权结果和算术平均几乎一样可能纬向稀疏或区域太小确认是否真的需要加权若区域跨越 20° 以上则必须加权与论文数值对不上网格方向不同或变量层级、单位不同核对维度和单位必要时 permute计算长时间序列太慢循环体内反复构造权重把权重构造移出循环或改用矩阵乘法4.2 纬度单位搞混cosd 还是 cos这个问题值得单独强调一下因为太容易犯了。cos和cosd的差别就是弧度与角度的差别。很多模式输出的纬度是度像 30、45、60 这种。如果你写了cos(lat)Matlab 会把这些数字当成弧度来算结果荒谬到你根本看不出问题。比如cos(60)算出来是 -0.95而正确的cosd(60)是 0.5。负的权重一旦出现加权平均有时候比算术平均还离谱。建议在构造权重之后直接画个图看一眼plot(lat, wgt(:,1))如果曲线在 -90 到 90 之间都是非负、且在两极接近 0、在赤道最大那基本没问题。如果出现负数或不合理的波动马上检查单位。4.3 NaN 处理权重必须同步清零数据里有 NaN 是气象数据常态地面站插值产品会有缺测模式输出偶尔也会因为地形原因出现无效值。我记得第一次写面积加权时没有处理 NaN直接sum(data .* wgt, all) / sum(wgt)结果整条全球平均序列的前十年全是 NaN排查后发现是南极某个格点在特定时段缺测。后来的处理原则就一条凡是数据为 NaN 的格点在加权平均时视为不存在既要把数据剔除也要把对应权重清零再用剩余权重重新归一化。反映在代码里就是上面函数中的valid和w(~valid) 0那两行。4.4 性能优化思路如果你手上是几百年长度的日数据逐层循环算面积加权确实比较慢。最有效的优化是如果数据里没有 NaN或者你提前做了 NaN 填充那么可以直接用矩阵乘法一步搞定。% 前提data2d 中已将 NaN 填为 0且 wgt_vec 已归一化 data2d(~isfinite(data2d)) 0; avg wgt_vec * data2d; % 1 x time这条语句的内部实现是优化的线代运算比 for 循环快一到两个数量级。如果是数据量太大、内存装不下那就用matfile对象按时间分块读取每块调用加权平均函数最后拼接结果。还有一个思路是降采样重采样到等面积网格再做算术平均但那样会引入插值误差远不如直接加权来得干净。5. 我的习惯做法与扩展建议5.1 把加权平均封装成自己的工具箱函数平时我处理气象数据不会每次都手写一遍面积加权代码。我更建议你把这个能力沉淀成一个独立函数就像上面那个area_weighted_mean一样放到自己的 Matlab 路径下之后任何项目都能直接调用。你可以在函数注释里记录清楚输入输出、单位要求和使用案例这样半年后再看也还能秒懂。如果你处理的数据源比较多比如同时用 CMIP6、ERA5、MERRA-2可以再封装一层接口把不同数据源自动 permute 到统一的[lat x lon x time]维度再统一调用加权函数。这套工作看起来很小但能帮你省掉大量重复劳动也让自己的分析流程更规范、更好复现。5.2 面积加权还能用在哪些分析里除了常规的区域平均面积加权在很多衍生分析里都适用。做 EOF 分析时如果直接用原始格点场计算协方差矩阵高纬网格会被过度代表导致空间模态偏向高纬特征。所以做 EOF 之前应当先给每个格点乘上面积权重的平方根做完 EOF 再把空间模态还原回去这一步在气候统计里面非常关键。计算区域总降水量或总热通量的时候用的是另一个变体总面积求和而不是面积加权平均。也就是说[ \text{Total} \sum_i X_i \cdot A_i ]比如算整个流域的年总降水量就需要把每个格点降水率乘以格点面积再对时间积分。这种情况下cosd(lat)同样可以用来构造面积权重只是最后是按权重乘以数据后再求和而不是除以总权重。还有模式评估和偏差分析。把模式插值到观测网格后如果要计算全球或区域的平均偏差同样需要面积加权。不然北极增暖信号强的模式在算术平均下会显得偏差特别大而实际上它影响的面积可能很有限。5.3 后续可以进一步做的事如果不想每次自己算面积也可以考虑引入成熟的气候数据工具箱。MATLAB 的 Mapping Toolbox 有一些地理统计相关的函数比如areaquad可以计算经纬度四边形面积不过要自己写循环。第三方工具包里Climate Data Toolbox for MATLAB 里提供了areaweight之类的函数处理网格面积时比较方便推荐对工具箱不反感的读者尝试。另一个可以拓展的方向是结合插值把不同分辨率的模式数据先统一插值到同一套网格再做面积加权平均。这个过程在 CMIP6 多模式集合平均里几乎是必经之路。插值本身会引入误差所以通常建议先在原始网格上算好区域平均需要做模式间比较时再插值到公共网格这样能尽量减少信息损失。最后分享一个小技巧我做完任何一套面积加权平均之后一定会顺手把全球平均温度或者降水的时间序列和主流数据集报告值做一次散点对比。比如全球平均温度应该在 14°C 到 15°C 左右如果算出来是 10°C 或者 20°C那一定不是加权的问题就是数据掩膜或单位的问题。这种基准检查虽然土但在实际工作中往往是最快能发现问题的方法。希望这篇关于面积加权的整理能对你有用咱们下一篇继续聊气象数据处理的其他细节。