这个系列写到第六篇了。前面几篇我们聊过NetCDF读取、时间维处理、重插值、气候态计算这些基础操作今天这篇我想认真聊聊一个很多人第一眼觉得“不就是求平均嘛”但实际上一不留神就会出问题的环节——数据的面积加权。做气象的人都清楚我们每天经手的再分析资料、模式输出、卫星反演产品绝大多数存的是经纬度格点。只要你是用经纬度网格做区域平均、全球平均、或者把网格数据聚合到流域尺度就必须面对“每个格点代表的地面面积不一样”这件事。这篇文章我会把面积加权从原理、Matlab实现到验证方法、进阶场景和常见坑一次讲明白希望能帮你少走点弯路。1. 为什么区域平均必须算上“每个格点面积不同”这笔账1.1 经纬度网格的本质格点面积随纬度变化先看一个最基础的问题一个按 1°×1° 划分的经纬度网格赤道上的一个格子和北纬 60° 的一个格子覆盖的地表面积一样吗答案是不一样而且差别非常大。经度方向的两条经线从赤道向两极靠拢所以纬度越高单个格点在纬向的实际距离就越短。计算一个格点面积的标准公式是[ \Delta A R^2 \cos\phi \cdot \Delta\lambda \cdot \Delta\phi ]其中 (R) 是地球半径(\phi) 是纬度(\Delta\lambda) 和 (\Delta\phi) 分别是经向和纬向的网格步长。对于等经纬度网格(\Delta\lambda) 和 (\Delta\phi) 固定因此面积权重只和纬度的余弦有关。我做个表方便你直观感受。假设地球半径按 6371 km 算1°×1° 网格在不同纬度处的近似面积如下纬度°N(\cos(\text{lat}))1°×1° 网格近似面积(km^2)与赤道面积之比01.000123651.00300.866107100.87450.70787450.71600.50061830.50800.17421470.17北纬 60° 的格子面积只有赤道的一半到了极区附近已经不足赤道的五分之一。如果你拿着一份全球数据漫不经心地做算术平均就等于把极区的格点当作和赤道格点一样重要这显然不合理。1.2 算术平均和面积加权的结果能差多少很多人会问这个误差真的重要吗我用实际数据试过答案是“看情况但通常不能忽略”。先说温度。曾经有一段时间我处理 ERA5 的 2m 气温做全球平均直接用mean(temp_global(:))算出来比用面积加权平均的高出大约 0.30.5℃。原因不复杂高纬度地区格点面积小但数量却不比低纬度少如果某个季节北极异常偏暖算术平均会把这个偏暖的贡献放大加权平均则会因为高纬面积权重的降低而更接近“真实的地表平均温度”。降水更明显。降水主要集中在热带和副热带热带格点的面积权重本来就大如果你用算术平均去算一个横跨南北半球的区域平均降水结果往往会偏向中高纬度那些面积小但数值不一定小的噪声格点。我做过一次 30°S30°N 平均降水的测试算术平均和面积加权平均能差出 0.150.3 mm/day。对于研究气候变化趋势来说这种量级的偏差完全可能掩盖真实的信号。1.3 什么时候可以偷懒不用面积加权我写这么多不是为了吓唬你实际上有些场景确实可以近似处理。你研究的是一个小区域比如单个城市、一个省的局部范围纬度跨度不超过 23°面积差异不明显此时算术平均和面积加权平均差别很小。数据本身已经投到等面积投影网格上比如很多区域模式输出用了兰伯特投影或极地立体投影这时网格面积已经均匀继续做面积加权反而是错。你用的产品已经在生产阶段做过面积平均比如某些官方发布的指数序列文档里通常写清楚了“area-weighted average”这时候不需要再重复加权。但我的习惯是除非你非常有把握否则一律加权。面积加权在 Matlab 里实现起来成本很低没必要拿结果准确性去赌。2. Matlab 里构造面积权重矩阵的几种做法2.1 最常用的中心点余弦权重公式最常见的做法就是用格点中心纬度计算权重% lat 是网格纬度向量单位度 weight cos(deg2rad(lat(:)));这份权重本身是 1×nlat 的向量最后要扩展成二维矩阵参与运算% 假设 data 维度是 lon × lat× time且 lat 对应第二维 nlon numel(lon); nlat numel(lat); weight2d repmat(weight, [1, nlon]); % 维度变成 nlon × nlat % 或者直接用隐式扩展MATLAB R2016b 以后都支持 weight2d repmat(cos(deg2rad(lat)), nlon, 1);然后计算加权平均就很简单area_mean sum(data .* weight2d, all) ./ sum(weight2d, all);这里有一个非常容易被忽略的细节cos(deg2rad(lat))得到的是“相对面积权重”不是绝对面积。因为我们计算的是平均值的分子分母权重比最终归一化的时候绝对面积常数会被约掉所以直接用余弦权重是没问题的。2.2 如果网格提供的是“中心坐标”而非“边界坐标”如何更精确很多 NetCDF 数据里的lat变量是格点中心纬度不是格点边界。此时上面的余弦权重已经足够精确因为网格步长固定时面积公式中的余弦项取中心纬度是二阶精度。但如果你追求更高精度或者网格不均一比如高斯网格、变分辨率网格更好的做法是显式计算每个格点的面积权重R 6371000; % 地球平均半径单位米 dlon abs(lon(2) - lon(1)) * pi / 180; % 由格点中心纬度构造格点边界纬度 dlat abs(lat(2) - lat(1)); lat_edges [lat(1) - dlat/2, (lat(1:end-1) lat(2:end))/2, lat(end) dlat/2]; lat_edges lat_edges * pi / 180; % 每个格点的面积单位m^2 lon_width dlon; % 经向宽度弧度度转换为弧度已在前面对 dlon 处理 area R^2 * lon_width .* (sin(lat_edges(2:end)) - sin(lat_edges(1:end-1))); area repmat(area(:), nlon, 1); % 扩展成 nlon × nlat这个公式来自球面面积微元的积分在纬度 (\phi_1) 到 (\phi_2) 之间、经度 (\lambda_1) 到 (\lambda_2) 之间的面积是[ A R^2(\lambda_2 - \lambda_1)(\sin\phi_2 - \sin\phi_1) ]用这个方式算出来的权重天然包含了 (1^\circ) 网格在边界上的“棱角效应”比简单乘 (\cos(\text{中心纬度})) 更接近真实面积。比如在 60°N 附近精确面积法和中心余弦法的差异虽然不大但当你统计长时间序列时这种系统性误差会被累积放大。2.3 权重矩阵的存储与复用权重只依赖网格经纬度不依赖时间。如果你在循环里反复读数据、反复计算权重那就是浪费。我一般会在第一次处理某个数据集时把权重矩阵存成一个.mat文件save(area_weights_era5_0.25x0.25.mat, weight2d, lat, lon, area, -v7.3);之后只要数据网格不变就直接加载load(area_weights_era5_0.25x0.25.mat, weight2d);这里提醒一句如果是全球 0.25°×0.25° 的数据weight2d是 1440×720 的 double 数组大约 8.3 MB保存和加载都没什么压力。但如果你同时存了多个变量建议用-v7.3格式后续扩展大数组更稳妥。3. 实际动手用面积加权算一个区域的逐月平均降水3.1 数据准备与目标区域截取以一份常见的全球逐月降水数据为例。假设文件是precip.mon.mean.nc变量名precip维度顺序是(lon, lat, time)其中lon从 0 到 359.5lat从 90 到 -90有些数据纬度是递减的这个后面会专门说。我先用 ncread 读取需要的部分lon ncread(precip.mon.mean.nc, lon); lat ncread(precip.mon.mean.nc, lat); precip ncread(precip.mon.mean.nc, precip); % 目标区域例如中国大陆周边的 20°N-45°N, 75°E-135°E lonlim [75, 135]; latlim [20, 45]; ilon find(lon lonlim(1) lon lonlim(2)); ilat find(lat latlim(1) lat latlim(2)); precip_sub precip(ilon, ilat, :); lon_sub lon(ilon); lat_sub lat(ilat);这里要特别注意经度范围。中国区域在 75°E-135°E 这个区间一般不会跨 0° 或 360° 分界所以直接find就行。如果目标区域跨国际日期变更线比如 120°E-120°W你大概率需要把经度统一到 0360 再选后面我会专门写这个坑。3.2 面积加权平均的完整实现步骤有了precip_sub和lat_sub接下来就进入正题。假设precip_sub维度是nlon × nlat × ntime纬度方向是递增还是递减要先确认好。我这里写成通用代码核心逻辑是把权重广播到三维并正确处理缺测nlon numel(lon_sub); nlat numel(lat_sub); ntime size(precip_sub, 3); % 面积权重先算 1 × nlat 的向量 weight_1d cos(deg2rad(lat_sub(:))); % 扩展为 nlon × nlat × 1方便和三维数据直接相乘 weight_3d repmat(weight_1d, [nlon, 1, 1]); weight_3d repmat(weight_3d, [1, 1, ntime]);repmat扩展三维数组在 ntime 很大时容易爆内存。1440×720×120 个 double 已经接近 1 GB所以更推荐用隐式扩展配合向量化% 设定有效数据掩膜 valid ~isnan(precip_sub); precip_sub(~valid) 0; % 分子加权后的累积值分母有效权重和 weight_3d repmat(weight_1d, [nlon, 1]); % 先只做二维扩展 weight_3d reshape(weight_3d, [nlon, nlat, 1]); data_score sum(precip_sub .* weight_3d .* valid, [1, 2]); weight_score sum(weight_3d .* valid, [1, 2]); area_avg squeeze(data_score ./ weight_score);这段代码的逻辑是把缺测值先换成 0同时用valid记录哪些位置原本有数据。分子只累计“有数据的位置”的加权值分母也只累计“有数据的位置”的权重。如果整个区域在某个时次全部缺测weight_score会是 0结果自然变成 NaN不会出现除零问题。如果你觉得这样绕也可以用循环逐时次算慢一点但更好理解area_avg nan(ntime, 1); for t 1:ntime x precip_sub(:, :, t); w repmat(weight_1d, nlon, 1); mask ~isnan(x); x(~mask) 0; area_avg(t) sum(x .* w, all) / sum(w .* mask, all); end这种逐时间循环在 MATLAB 里确实不快但胜在不会因为三维repmat把内存打爆。如果你嫌慢可以每 10 年一个块处理平衡内存和速度。3.3 结果验证拿算术平均对比再用“全 1 测试”兜底拿到area_avg之后不要急着用。我最常用的验证方法是“全 1 测试”把precip_sub全部替换成 1正确的结果应该仍然是 1。如果算出来不是 1说明权重或维度匹配有问题。test_data ones(nlon, nlat, ntime); test_avg squeeze(sum(test_data .* weight_3d .* valid, [1, 2]) ./ sum(weight_3d .* valid, [1, 2])); disp(unique(round(test_avg, 10)));另一个更实际的验证是和算术平均对比。下面是同一组目标区域数据我用两种算法得到的不同月份均值差异示意月份算术平均mm/day面积加权平均mm/day差异mm/day1月2.872.74-0.134月3.213.300.097月5.425.18-0.2410月3.083.120.04差异看起来不大但如果你是做长时间趋势分析这种系统性偏差会影响统计显著性。尤其是降水这类空间分布极不均匀的要素高分辨率数据上的差异会更明显。如果你有机会用 CDOClimate Data Operators可以交叉验证一下。CDO 的fldmean默认是面积加权平均使用的正是网格面积权重和我们的算法结果应该非常接近。我用 0.25° 数据对比过两种方法算出来的区域平均序列相关系数无限接近 1数值差通常小于 1e-5。这种独立工具交叉验证是最让人放心的。4. 面积加权不止用于区域平均几个高价值场景4.1 全球平均温度与温度异常序列全球平均温度序列是气候监测的核心指标几乎所有机构发布的全球温度产品都会明确说明自己用了面积加权。原因前面已经说过如果直接用算术平均高纬度格点权重过大而高纬度又是升温最快的区域这会显著高估全球增暖幅度。用 Matlab 计算全球平均温度的代码和区域平均几乎一样只不过区域范围变成全球lat ncread(tas_aiwg_1x1.nc, lat); lon ncread(tas_aiwg_1x1.nc, lon); tas ncread(tas_aiwg_1x1.nc, tas); weight cos(deg2rad(lat)); weight2d repmat(weight(:), nlon, 1); % 逐时次加权平均 ntime size(tas, 3); tas_global nan(ntime, 1); for t 1:ntime x tas(:, :, t); mask ~isnan(x); x(~mask) 0; tas_global(t) sum(x .* weight2d, all) / sum(weight2d .* mask, all); end这里再强调一点如果你用的是 NCEP/NCAR 再分析数据它的网格可能是高斯格点。高斯格点虽然不需要余弦权重但官方会给一个一维权重数组gweight你直接用它当作weight去扩展二维矩阵即可不要再乘余弦。否则等于对高纬度权重进行了二次压缩结果会有偏。4.2 三维场逐层面积加权再垂直积分很多时候我们不只是算地面变量还会处理 3D 场比如温度、位势高度、比湿、风速等。面积加权在这类数据里只作用于水平维度lon-lat垂直维度和时间维度都不参与。假设数据data维度是(lon, lat, lev, time)你想得到每个高度层上的区域平均nlev size(data, 3); ntime size(data, 4); plev_avg nan(nlev, ntime); weight2d repmat(cos(deg2rad(lat)), nlon, 1); for k 1:nlev x squeeze(data(:, :, k, :)); % nlon × nlat × ntime % 这里可以调用前面写好的加权函数 for t 1:ntime x_t x(:, :, t); mask ~isnan(x_t); x_t(~mask) 0; plev_avg(k, t) sum(x_t .* weight2d, all) / sum(weight2d .* mask, all); end end得到各层的区域平均后再做垂直积分比如整层水汽通量、位势厚度注意面积加权和垂直积分要分层完成不能把三维数组直接全部乘一个标量权重因为垂直层次之间没有面积关系。我见过有人把权重矩阵扩展成三维后对整个三维数组求和这其实是把垂直层也当成水平维度平均掉了结果会完全偏离定义。4.3 流域面雨量估算面积加权思想的变体水文上的面雨量和气象上的区域平均本质是同一个问题只是权重来源不同。气象上通常用经纬度格点面积作为权重水文上常用泰森多边形法或多边形裁剪法把流域内每个站点/格点的控制面积求出来再做加权平均。在 Matlab 里如果你有流域边界和一份格点降水数据最简单的做法是先用inpolygon生成流域掩膜in inpolygon(lon2d, lat2d, basin_lon, basin_lat); mask double(in); % 权重 格点面积 × 流域掩膜 weight_basin weight2d .* mask; precip_basin sum(precip .* weight_basin, all) / sum(weight_basin .* ~isnan(precip), all);mask的存在相当于把流域之外的格点权重清零。如果你用的是高分辨率格点数据这种方法比先插值到站点、再做泰森多边形要简单而且误差更可控。我个人的经验是当流域形状复杂、跨越多个纬度带时这种基于掩膜的面积加权面雨量比简单算术平均稳定很多尤其适合做长时间序列的水文变化分析。5. 绕过这些坑面积加权才算真正落地5.1 NaN 数据的处理顺序面积加权最怕的就是缺测值。很多初学者会这么写bad mean(data .* weight2d, all) / mean(weight2d, all);如果数据里存在 NaNdata .* weight2d的对应位置也变成 NaN整个mean直接就变成 NaN。即便是用了omitnan也仅仅是忽略了该位置的分子分母的权重没被同时忽略结果照样有偏。正确的方式是先构造掩膜再用掩膜同时更新分子和分母mask ~isnan(data); data_safe data; data_safe(~mask) 0; numer sum(data_safe .* weight2d, all); denom sum(weight2d .* mask, all); result numer / denom;这也是我前面在代码里反复强调weight2d .* mask的原因。这一步看似简单但却是我帮不少人调代码时发现的最常见问题。5.2 纬度方向、经度跨越特殊边界的处理很多数据集里的纬度是从 90 到 -90 递减的例如 NCEP/NCAR。如果你直接用lat -30 lat 30去选在递减向量上并不会出错但如果你假定纬度递增权重计算时cos(deg2rad(lat))的方向就会反导致后续repmat出来的权重矩阵和数据的纬度维不匹配。我的做法是无论输入怎样先统一成递增方向if lat(1) lat(end) lat flipud(lat); data flip(data, 2); % 假设纬度是第二维 end经度边界更隐蔽。比如你想选太平洋区域 120°E-120°W在 -180 到 180 的表达里不能简单用一个lon 120 lon -120因为没有任何数能同时满足。这时最好把经度统一到 0360 的范围内lon360 mod(lon, 360); data360 cat(1, data(lon360 120, :, :), data(lon360 240, :, :)); % 因为 120°W 240°E或者反过来把数据重复拼接一次再做掩膜。总之先搞清楚数据的经度区间再决定截取方式不要在find这一步图省事。5.3 内存与性能优化心得处理全球高分辨率数据时内存往往是比速度更先遇到的瓶颈。我建议几条实战经验不要动不动就repmat三维大数组。能用二维权重就不扩三维能隐式扩展就用隐式扩展或者干脆循环。权重矩阵用单精度存。面积权重的精度不需要 double单精度足够内存直接减半。只是注意single与double混合运算时会自动转成 double别让白省的内存在运算时又涨回去。时间循环如果太慢可以每 1020 个时次一组分批向量化。比如repmat权重到(nlon, nlat, 10)然后一次性处理 10 个时次再进入下一组。用squeeze和permute之前先想清楚维度顺序。我最怕看到一堆squeeze把原本清晰的维度搞乱最后排错比写代码还累。5.4 快速 debug 技巧最后分享几个我调试时必用的自检手段症状可能原因检查方法加权平均结果和算术平均几乎一致权重没有参与运算或权重全为常数打印min(weight(:))和max(weight(:))看是否随纬度变化高纬度区域权重异常偏大cos和deg2rad顺序错了或把经度当成了维度画imagesc(lon, lat, weight2d)看权重沿纬度是否平滑递减结果全是 NaN分子分母的掩膜没设置好或有 NaN 污染检查sum(valid(:))是否为 0纬度方向和想象相反数据 lat 递减但代码按递增处理统一 flip 之后再算一次对比结果我习惯在正式跑任何区域平均之前先画一张权重图看一眼高纬度颜色是否明显比低纬度浅。如果权重图不对后面根本不用往下算。面积加权这件事说难不难说简单又总有人在不该出问题的地方栽跟头。我个人最深的体会是权重矩阵一定要复用别每次循环都去重新生成缺测掩膜一定要记得参与分母最终结果一定要用全 1 测试兜底。把这三点做到位你在 Matlab 里算出来的面积加权平均基本就不会有方向性错误了。如果后续遇到更奇葩的网格类型也欢迎回来一起讨论。