MATLAB实现EOF与REOF：从SVD分解到Varimax旋转的完整流程

简介一套面向气象、海洋等地球科学数据分析的MATLAB工具包聚焦EOF、REOF、SVD与CCA四种常用统计方法适合需要做多变量降维、时空模式识别或区域异质性研究的学生和科研人员使用。压缩包内含4个.m文件分别对应经验正交函数、区域经验正交函数、奇异值分解和典型相关分析的实现整体仅4KB代码精炼便于阅读和二次修改。已有323人学习下载。通过这套代码可以快速掌握EOF/REOF的完整处理流程从数据预处理、SVD分解、EOF与负荷计算到按空间划分子区域、分别分析并组合结果同时理解CCA在两变量集相关性分析中的应用。脚本既可独立运行也可串联使用方便根据实际数据和研究目标灵活调整保留模态数或区域划分方式是地学统计入门与日常分析的实用工具。1. 这串文件名背后的 EOF/REOF 分析流程看到eof,reof等.rar这类资源很多人的第一反应是解压、改名、直接跑REOF.m。真正决定分析质量的不是 RAR 里那几个.m文件而是你对 EOF、REOF 和 SVD 三者关系的理解。经验正交函数Empirical Orthogonal Function是一套把“空间点 × 时间序列”矩阵分解为空间型与时间系数的降维工具实际计算大多落在奇异值分解SVD上REOF 则是基于 EOF 前 k 个模态再做旋转让空间模态更集中、更便于物理解释。这套流程常见于气象海洋领域的网格数据、遥感影像序列、工业传感阵列与任何“批量位置 连续时间”的测量数据。适合的人群是拿到 MATLAB 脚本但不知道参数怎么设的数据工程师以及被要求把 EOF/REOF 结果写进报告却说不清 SVD 与前几个奇异值关系的分析师。下面直接从数值实现开始把每一步用最小可运行代码写出来。2. EOF 的数值骨架为什么用 SVD 而不是直接算协方差特征分解针对网格化数据做 EOF 时常见做法是先把观测场排成二维矩阵 X然后用svd(X,econ)一次拿到左奇异向量、奇异值和右奇异向量。之所以不走协方差矩阵eig不是因为两者结果不同而是因为数值可靠性。协方差矩阵 C Xc*Xc/(n-1) 是 m×m 实对称矩阵理论上特征值非负但网格数据量大、量纲差异明显时C 的条件数会放大舍入误差eig偶尔给出微小的负特征值后续求方差贡献时开根号会直接报错或产生 NaN。SVD 直接作用在 Xc 上奇异值永远不会为负“经济型分解”同时控制了内存占用是写 EOF 代码时的默认选项。2.1 数据矩阵的排列约定与 SVD 三因子的物理解读在气候分析里EOF 的数据矩阵一般写成X(m,n)m 是空间格点数n 是时间观测数时间维固定在第二维。这个约定与很多从 Fortran 转过来的脚本不同旧脚本常把三层do循环的 (经度,纬度,时间) 直接写进文件进 MATLAB 后按列优先展开空间维被拆散十个粗心的人里有八个会算完才发现第一模态对应的是“半个地球拼在一起”。保持时间维在后的好处从 SVD 三因子里直接看得出来| SVD 返回项 | 维度 | EOF 语义 | | --- | --- | --- | | U | m×min(m,n) | 空间载荷向量EOF 的空间型 | | S | min(m,n)×min(m,n) | 奇异值平方后对应特征值 | | V | n×min(m,n) | 归一化时间系数主成分方向 |左奇异向量 U 用于画空间型右奇异向量 V 用于画时间系数S 对角元用于计算各模态方差贡献。理解了这三者的角色就不会再把REOF.m里的旋转对象搞错。需要强调一点svd(Xc,econ)返回的列数取决于 m 与 n 的较小值当空间点数远大于时间样本数时返回的列数是 n此时第 n 个奇异值之后的空间变化根本进不了分解时间样本长度直接决定可解的模态数量。2.2 用 MATLAB 的 svd() 跑通一个最小 EOF 代码手头没有真实观测数据时我通常先合成一组含两个空间模态的信号这样能验证代码对奇异值分解的每个中间量都认识。下面的脚本不需要任何附加工具箱生成 60 个空间点、200 个时间样本的场噪声方差约为信号能量的四分之一。% 合成数据两个空间模态各自有独立时间系数 rng(2024); m 60; % 空间格点数经向或站点数 n 200; % 时间样本数 lon (1:m); % 模式1高斯型空间分布时间系数低频振荡 pattern1 exp(-((lon-20).^2)/50); pc1 sin(2*pi*(1:n)/40) 0.2*randn(1,n); % 模式2带正负位相的两个中心时间系数高频 pattern2 exp(-((lon-45).^2)/30) - exp(-((lon-12).^2)/20); pc2 cos(2*pi*(1:n)/15) 0.3*randn(1,n); X pattern1 * pc1 pattern2 * pc2; X X 0.5 * randn(m, n); % 附加噪声 % EOF 核心分解 Xc X - mean(X, 2); % 去掉各格点的时间均值 [U, S, V] svd(Xc, econ); % 方差贡献 lambda diag(S).^2 / (n - 1); frac lambda / sum(lambda); % 取前两个模态做重构检验 k 2; X_rec U(:, 1:k) * S(1:k, 1:k) * V(:, 1:k);这段代码的参数和逻辑说明如下mean(X, 2)沿第二维求均值得到每个格点的时间平均值Xc保存的是距平场。如果对整个矩阵只减去一个全局标量格点间的气候均值差异仍然保留第一模态会被平均水平空间分布占据。svd(Xc,econ)中的econ是经济型分解。m60、n200 时完整 SVD 会返回 60×60 的 U、60×200 的 S、200×200 的 Vecon把 S 和 V 的零空间部分截掉返回 60×60 的 S 和 200×60 的 V。网格数据跑到几万格点时这个参数省下的内存非常可观。diag(S).^2/(n-1)把奇异值还原成协方差矩阵的特征值因为 Xc USV而 XcXc/(n-1) US*S*U/(n-1)。漏掉/(n-1)的REOF.m脚本很常见表现是方差贡献偏大一个量级。重构校验X_rec是判断 EOF 分解是否正确的通用手段。用前两个模态重建的场均与Xc的空间趋势明显不符时先检查数据生成逻辑而不是怀疑 SVD 算法。2.3 EOF 的预处理参数去均值、去趋势与空间权重把代码跑通后真正需要跟业务交互的是进入svd之前的三个参数选择。第一是去均值。前面用的是格点时间均值更严格的气象分析会强调“气候态”即每个格点在相同日期上的多年平均值。用总平均代替逐日或逐月气候态会留下明显的年循环使第一模态总是一个单调变化的信号。第二是标准化。单独分析一种变量比如只做海温通常不必对行做标准化但要同时分析温度、风速、气压这类量纲不同的变量时必须对每一行做 z-score否则数值大的变量会主导整个 EOF 分解而这不代表它的物理重要性。第三是空间权重。对等经纬度网格高纬度格点代表的实际面积小一般先乘以sqrt(cos(lat))权重再进 SVD。很多人用纬度平均面积去修正结果但修正动作若不放在分解前而放在分解后会直接改变模态的形状。这三个参数没有绝对对错但要在报告里写清楚因为 EOF 对预处理相当敏感同一套数据可能因为这三种选择的不同而得到差异明显的第一模态。3. REOF 步骤的实质Varimax 旋转、载荷符号与模态可解释性REOF 是 EOF 的旋转延伸不是新算法。REOF.m拿 EOF 前 k 个载荷向量做旋转得到空间上更集中的新载荷。最常见的准则是 Varimax它最大化载荷平方的方差。这么做的原因是 EOF 的正交性约束很强第二个模态必须与第一个模态在空间上正交这会导致模态变成“两个区域反相”的全球型图案区域化特征反而被抹掉。旋转放松正交约束后载荷的平方分布更趋于集中在少数格点上更接近人们对特定气候区的直观判断。注意旋转不改变原始场的信息量只改变前 k 个模态的解释方式。3.1 为什么旋转的是载荷而不是时间系数搞清旋转对象是所有REOF.m能否正确工作的分水岭。EOF 计算得到空间载荷 U 与时间系数 V 的乘积旋转通常只作用于载荷 A U(:,1:k)。旋转矩阵 R 是 k×k 矩阵新载荷 A_rot A*R。时间系数不能用旧 V 直接代替需要用 A_rot 对原始距平场做最小二乘投影即 PC_rot A_rot \ Xc。如果某个脚本旋转后只画了新载荷时间曲线还是老 V 的列它的时空配对就错位了。判断脚本是否可靠可以看它有没有求A_rot \ Xc这一步。3.2 MATLAB 里用 rotatefactors 做 Varimax 旋转与参数设置MATLAB 自带rotatefactors函数可以直接对前 k 个载荷做正交旋转。一个最小用法如下% 用前 k 个 EOF 载荷做 Varimax 旋转 k 4; % 旋转模态数需要根据谱确定 A_raw U(:, 1:k); % 原始空间载荷矩阵 (m*k) % Normalize 开启 Kaiser 归一化对大尺度场更稳健 A_rot rotatefactors(A_raw, Method, varimax, Normalize, true); % 求旋转后的时间系数最小二乘投影 Xc X - mean(X, 2); PC_rot A_rot \ Xc; % 左除等价于求最小二乘解 % 旋转后各模态实际方差贡献 var_rot sum(PC_rot.^2, 2) / (size(X, 2) - 1); frac_rot var_rot / sum(var(Xc, 0, 2));rotatefactors默认把每列当作因子、每行当作观测所以A_raw的行是空间格点、列是模态。Normalize, true对应 Kaiser 归一化旋转前先消除各格点载荷向量模长差异带来的尺度干扰旋转完再恢复对前几个模态方差贡献相差较大的场更合适。A_rot \ Xc返回 k×n 矩阵每一行是一个旋转模态的时间系数这一步保证了新时间曲线与新空间载荷匹配。frac_rot用旋转后时间系数的方差除以原始场总方差是 RE 模态的真实贡献不要直接沿用frac(1:k)。如果不想用工具箱也可以手工实现 Varimax对 A 迭代执行“计算载荷平方方差梯度、按 Givens 旋转更新、直到目标函数变化小于阈值”但手工实现需要处理迭代步长和正交约束业务中基本没有理由重复造轮子。3.3 REOF 的模态数、方差贡献与符号约定旋转模态数是 REOF 里争议最多的参数。取少了区域特征混合不充分取多了会把噪声也旋转出“看起来像模态”的图案。经验上先用 North 准则圈定显著 EOF 模态的个数特征值 λ 的采样误差约为 λ·sqrt(2/N)N 是时间样本数当相邻两个特征值的差值小于这个误差区间时两个模态不可区分旋转时要慎重。之后在前后两个整数上各做一次旋转比较空间载荷的相关性载荷对 k 的变化不敏感的区间就是可靠区间。不要为了追求“前两个模态解释 80% 方差”而强行把 k 设成 2因为旋转模态上的方差会重新分配这个数字并不代表区域模态的独立性。对比项EOFREOF模态正交性严格正交旋转后无正交约束空间分布可能全球型、载荷分散更集中、区域化明显方差贡献随奇异值递减需根据旋转后时间系数重新计算时间系数V 与 S 的乘积对旋转载荷做最小二乘投影主要用途查看总体变率识别可物理解释的区域模态符号约定也常被忽视EOF/REOF 的载荷符号是任意的同一分析用不同版本 MATLAB 计算得到的第二模态可能整体反号。写报告前要做一个符号固定指定某个区域中心点或某个载荷绝对值最大的格点的符号为正然后对整个模态和时间系数同时乘以 -1。缺少这一步两台机器上跑同一份数据画出来的图颜色会相反分析结论却完全一致容易引起不必要的返工。4. 把 RAR 里的 EOF/REOF 脚本落地成可复现流程拿到eof,reof等.rar这类资源后先别急着运行。RAR 包里的文件命名往往带着不同历史阶段的分析痕迹REOF.m、eof.m、plot_EOF.m、data.mat可能来自不同时期函数接口未必能直接串起来。我一般先把 RAR 解压到一个独立工作目录在 MATLAB 中用addpath指向该目录然后逐文件查看函数签名确认输入输出变量名后再组装流程。可以关注reof.m是否包含rotatefactors调用以及旋转后的时间系数是否重新投影过这两点最容易暴露“伪 REOF”。4.1 从三维网格到二维矩阵RAR 数据文件的导入与重构先看包里文件类型如果是.mat文件直接whos -file data.mat查看变量维度如果是.txt或.dat用readmatrix读取但要注意文件编码。很多早期脚本用 GBK 编码保存说明文字直接readmatrix不影响数值列但用textscan或fopen读取时若报出unexpected end of file或premature EOF优先检查文本文件是否完整解压其次看文件行数是否与脚本里的n一致。这类错误不是算法问题是数据读取层的问题。三维网格数据一般长这样% 假设数据为 (lon, lat, time)维度 [nlon, nlat, nt] [nlon, nlat, nt] size(X3d); % 将空间维拉平时间维保持最后一维 X2d reshape(X3d, nlon*nlat, nt);如果读进来发现维度是 (time, lon, lat)要先permute(X3d, [2 3 1])再reshape。这个 permute 步骤看着不起眼却决定后续空间图的绘制能否对上坐标。拉平后建议保留lon2d repmat(lon, nlat, 1)和lat2d repmat(lat, 1, nlon)之类的网格坐标数组防止画图时重新造坐标出现转置错误。4.2 一个通用的 EOFREOF 流水线脚本骨架在真实数据分析中我一般会写一个函数把 EOF、REOF 和方差贡献封装起来参数只留一个旋转模态数。这样无论换数据还是调参都不必改动分解核心代码。function [eof_load, eof_pc, eof_frac, ... reof_load, reof_pc, reof_frac] eof_reof_pipeline(X, kRot) % X: 空间格点×时间样本的距平场建议先在外面去掉气候态 % kRot: 参与旋转的模态数 Xc X - mean(X, 2); [U, S, V] svd(Xc, econ); lam diag(S).^2 / (size(X, 2) - 1); eof_frac lam / sum(lam); eof_load U(:, 1:kRot); eof_pc S(1:kRot, 1:kRot) * V(:, 1:kRot); % 保留振幅的时间系数 reof_load rotatefactors(eof_load, Method, varimax, Normalize, true); reof_pc reof_load \ Xc; reof_frac sum(reof_pc.^2, 2) / (size(X, 2) - 1) / sum(eof_frac); end函数里有两个细节。第一eof_pc没有直接用V(:,1:kRot)而是乘了S(1:kRot,1:kRot)因为V的列是归一化的乘上奇异值后时间序列才保留真实振幅不乘奇异值的后果是画时间系数图时数值普遍小于 0.1看不出发散趋势差异。第二reof_frac的分母用sum(eof_frac)与var(Xc(:))等价但前者更不容易因为单位换算出错。旋转模态的方差贡献不是lam(1:kRot)原样搬运这一点务必在输出阶段重新计算。4.3 三个必调参数与两个排错点必调参数至少有三个整理成表参数建议做法作用kRot依据 North 准则或特征值谱转折点选取控制旋转模态个数NormalizetrueKaiser 归一化提高旋转结果的可靠性左除投影reof_pc reof_load \ Xc保证旋转后的时空配对正确排错点一输入矩阵包含 NaN。SVD 对 NaN 零容忍任何一格缺测都会让对应奇异值变 NaN。不要用 0 填充缺测这会人为压低估方差使前几个模态虚高更常见做法是用该格点序列的多年平均插补或直接删除缺测行。排错点二奇异值非零但特征值出现负值。若有人坚持使用eig(Xc*Xc)并开根号负特征值会直接报错改用svd后此问题消失。若svd依然报出收敛警告检查数据是否含 Inf以及是否在分解前对行列做了不合理的标准化。5. 三个必须掌握的验证手法重构残差、符号固定与分半鲁棒性检验最后一章落在“怎么判断算对”上。EOF 和 REOF 没有真实标签可以对照输出图再漂亮也可能在符号、模态数或空间权重上出偏差所以我通常固定做三个验证。5.1 重构残差检验用前 k 个模态重建距平场残差里若还残留明显的空间主模态说明 k 取小了。判断残差是否随机可以对残差再做一次 SVD看第一模态方差占比是否显著高于后续模态如果还是“悬崖式”下降继续加模态。% 残差 EOF检查残余空间结构 Xc X - mean(X, 2); X_rec U(:, 1:k) * S(1:k, 1:k) * V(:, 1:k); resid Xc - X_rec; [Ur, Sr] svd(resid, econ); frac_resid diag(Sr).^2 / (size(X, 2) - 1) / sum(diag(Sr).^2); bar(frac_resid(1:min(10, numel(frac_resid))));5.2 符号固定EOF 的载荷符号是任意的换机器、换 SVD 库都可能整体反号。实务中我一般先指定一个物理基准点比如某个已知异常中心的格点索引检查该点在每个模态的载荷符号如果为负就把该模态的载荷和时间系数同时乘以 -1for j 1:k if eof_load(anchor_idx, j) 0 eof_load(:, j) -eof_load(:, j); eof_pc(j, :) -eof_pc(j, :); end end注意旋转前后的模态必须使用同一套符号规则否则 REOF 的空间型和时间系数会错配。5.3 分半鲁棒性检验把时间序列切成前后两半各跑一次 EOF/REOF比较对应模态的空间相关是判断旋转模态是否可信的最快方法。对应模态相关系数高于 0.6 才认为可靠如果某个模态在前后半段里变形明显它很可能由几个特征值极为接近的 EOF 混合而成不宜在结论中重点讨论。n size(X, 2); half floor(n/2); [eof1_l, ~, ~, reof1_l] eof_reof_pipeline(X(:, 1:half), kRot); [eof2_l, ~, ~, reof2_l] eof_reof_pipeline(X(:, half1:end), kRot); corr_mat corrcoef([reof1_l, reof2_l]); % 检查 corr_mat(1:kRot, kRot1:2*kRot) 中对应模态的匹配程度这三个验证都不需要额外工具箱却能在输出报告前拦截绝大多数错误。实际项目里三分之二的 EOF/REOF 问题出在 SVD 之前的预处理而不是算法本身把重构残差、符号和分半一致性做成脚本里的固定步骤比换更复杂的旋转算法更值得先做。本文还有配套的精品资源点击获取