简介：这份资源面向地球科学、气象与海洋领域的学习者和科研人员，聚焦EOF、REOF、SVD与CCA等常用统计分析方法在MATLAB中的实现，帮助解决多变量数据降维、空间模式识别与区域气候特征提取等问题，适合具备一定MATLAB基础、需要复现或理解相关算法流程的中高级用户。压缩包共4个文件，均为m脚本，整体约4KB，分别对应EOF、REOF、SVD与CCA四类算法的代码实现，便于按需调用与二次修改。目前已有325人学习下载，说明其在相关方向具有一定参考价值。读者可从中获得从数据预处理、SVD分解、EOF与负荷提取，到区域划分、子区域分析及结果组合的完整思路，并理解奇异值分解在EOF中的关键作用，为区域气候变化或环境问题的模式分析提供可运行的脚本基础与排错参考。
1. 从一份 eof,reof 等.rar 说起：REOF 到底在解决 EOF 的什么痛点
如果你手头正好有一份名为eof,reof等.rar的压缩包，里面躺着REOF.m和一堆EOF相关脚本，那你大概率已经踩进了气象、海洋或者遥感领域最经典的一组时空分解工具里。EOF（经验正交函数）大家都不陌生，本质就是对时空场做 SVD 分解，把原始场拆成空间模态和时间系数，用前几个模态解释大部分方差。但真正做过实际分析的人都知道，EOF 有个让人头疼的毛病：前几个模态经常是"混"的，一个模态里同时装着两个物理过程，空间型还随区域漂移，解释起来全靠玄学。
REOF（旋转经验正交函数）就是冲着这个痛点来的。它不改变 EOF 分解的总方差，而是在前若干个模态张成的子空间里做一次正交旋转，让每个模态的空间载荷尽量向少数区域集中，物理意义更干净。这份REOF.m配合SVD和EOF脚本，基本就是一套从原始场到旋转模态的完整链路。这篇文章不假设你手里有源码包，只按标题里这几个关键词，把 REOF 的选型理由、REOF.m的实现步骤、SVD 在其中的角色、以及参数怎么设、坑在哪，一条条讲清楚。适合已经跑过 EOF、想进一步做区域分型的人，也适合刚拿到这份 rar 不知道从哪下手的新手。
2. REOF 的数学底子与 SVD 在其中的真实角色
2.1 EOF 为什么需要旋转：方差集中不等于物理解释清晰
EOF 分解的数学目标非常明确：找到一组正交的空间基，使得投影后的时间系数方差依次最大。用 SVD 写出来就是，给定一个去均值后的时空矩阵 (X)（行是时间，列是空间格点），做
[ X = U \Sigma V^T ]
其中 (V) 的列就是空间模态（EOF），(U\Sigma) 对应时间系数（PC）。前几个模态方差贡献大，这是"最优"的，但最优是针对方差而言，不是针对物理解释。问题出在正交约束上：正交是数学上的方便，不是自然界的规律。真实的大气或海洋过程在空间上往往不是正交的，两个相邻区域的异常可以同时出现，EOF 为了保持正交，只能把它们揉进同一个模态，或者拆成两个空间型相似、时间系数却正交的模态，这就是所谓的模态混合。
我一般会用一个很直观的判断：如果第二和第三模态的空间型长得像一对双胞胎，时间系数还差不多，那基本可以判定这两个模态被 EOF 强行拆开了，这时候就该上 REOF。旋转的目的不是重新分配方差，而是在前 (k) 个模态张成的子空间里换一组基，让每个基向量的载荷更"局部化"。总方差不变，解释方差的前几个模态会重新洗牌，通常第一个旋转模态的解释方差会下降，但空间型会干净很多。
2.2 旋转准则怎么选：Varimax 是默认答案但不是唯一答案
REOF 的核心是旋转准则。最常见的是 Varimax（方差极大旋转），它的目标函数是让每个模态的载荷平方的方差最大化，翻译成人话就是：让大的载荷更大、小的载荷更小，空间型要么强要么弱，不要温吞水。REOF.m里如果只实现了一种旋转，八成就是 Varimax。
Varimax 又分两种：正交旋转和非正交旋转。正交旋转保持模态之间仍然正交，实现简单，结果稳定，是绝大多数论文的默认选择。非正交旋转（比如 Promax）允许模态之间相关，解释上更灵活，但代价是模态不再正交，后续做回归或者合成分析时要小心。我的经验是，第一次做 REOF 就用正交 Varimax，跑通了、结果合理了，再考虑要不要换 Promax。别一上来就追求"更物理"，先把基线跑出来。
还有一个参数是旋转的模态数 (k)。这个没有标准答案，但有经验区间。一般取 EOF 前 10 到 25 个模态做旋转，具体看你的场有多少个有效自由度。取太少，旋转空间不够，模态还是混；取太多，会把噪声模态也拉进来旋转，结果反而不稳定。我通常的做法是：先看 EOF 的方差贡献曲线，找到拐点，拐点之前的模态数再加 5 到 10 个，作为旋转的 (k)。这个值在REOF.m里通常是一个输入参数，改起来很方便。
2.3 SVD 与 EOF 的关系：别把两个东西混为一谈
标题里同时出现了SVD和EOF，很多人会问：到底用 SVD 还是用 EOF？答案是，EOF 的数值实现通常就是靠 SVD。对时空矩阵做 SVD，左奇异向量对应时间系数，右奇异向量对应空间模态，奇异值的平方对应方差。所以SVD是手段，EOF是目的，两者不是并列关系。
但要注意一个坑：SVD 分解出来的模态顺序是按奇异值从大到小排的，这个顺序在旋转之后会被打乱。旋转后的模态不再按方差贡献排序，你需要自己重新算每个旋转模态的解释方差，然后手动排序。REOF.m里如果没做这一步，输出的模态顺序可能是乱的，直接拿去画图会闹笑话。我见过有人把旋转后的第一个模态当成方差最大的模态，结果解释方差只有 8%，还硬说这是主导模态，这就是没重新排序的后果。
另外，做 SVD 之前一定要去均值。EOF 分析的是异常场，不是原始场。如果原始场有很强的气候态，不去均值直接 SVD，第一个模态会是气候态本身，方差贡献可能高达 90% 以上，后面的模态全被压死。这个坑太常见了，但每年都有人踩。
3. 用 REOF.m 跑通一次旋转分解：从数据准备到模态输出
3.1 数据准备：时空矩阵的维度约定与去均值
在 MATLAB 里跑REOF.m之前，你得先把数据整理成它认识的格式。最常见的约定是：行是时间，列是空间。假设你有一个三维数组data，维度是[nt, nlat, nlon]，需要先 reshape 成二维：
% 假设 data 维度为 [nt, nlat, nlon] [nt, nlat, nlon] = size(data); % 重塑为 [nt, nlat*nlon]，行是时间，列是空间 X = reshape(data, nt, nlat*nlon); % 去均值：对每一列减去时间均值 X = X - repmat(mean(X, 1), nt, 1); % 处理缺测：如果某些格点全是 NaN，直接置零或剔除 X(isnan(X)) = 0;这段代码做了三件事：重塑、去均值、缺测处理。去均值那一步用repmat是为了兼容老版本 MATLAB，新版本可以直接用隐式扩展X - mean(X,1)。缺测处理要小心，如果某个格点缺测太多，直接置零会引入虚假信号，更好的做法是把这个格点整列剔除，或者用插值补全。我一般会先统计每个格点的缺测率，超过 30% 的直接不要。
提示：去均值之后，建议再检查一下每一列的均值是否接近零。如果某个格点均值还是很大，说明去均值没做对，后面 SVD 出来的第一个模态会很脏。
3.2 调用 REOF.m 的标准流程与参数含义
假设REOF.m的接口是[L, V, pc, var] = REOF(X, k, rotmode)，其中X是去均值后的时空矩阵，k是旋转模态数，rotmode是旋转准则（比如'varimax'）。一个典型的调用如下：
% 设置旋转模态数，一般取 10 到 25 k = 15; % 调用 REOF，rotmode 用 varimax 正交旋转 [L, V, pc, var] = REOF(X, k, 'varimax'); % L 是旋转后的空间载荷，V 是旋转矩阵，pc 是时间系数，var 是解释方差 % 重新按解释方差排序 [var_sorted, idx] = sort(var, 'descend'); L = L(:, idx); pc = pc(:, idx);这里k=15不是随便写的。如果你的场有 30 年逐月数据，时间样本 360 个，空间格点可能几千个，有效自由度大概在几十的量级，取 15 个模态旋转是合理的。如果k取到 30，会把太多噪声拉进来，旋转结果会碎掉。rotmode选'varimax'是默认，如果REOF.m支持'promax'，可以后面再试。
排序那一步非常关键。旋转之后，var里的解释方差不再是从大到小排的，必须手动sort。我见过太多人忘了排序，直接把第一个模态当主导模态，结果图一画出来，空间型乱七八糟，还以为是方法有问题。
3.3 解释方差重算与模态显著性检验
旋转后的解释方差不能直接用 SVD 的奇异值算，因为旋转改变了基向量，方差分配变了。正确的做法是：对每个旋转模态，用它的空间载荷去投影原始场，得到新的时间系数，再算这个时间系数的方差占总方差的比例。REOF.m如果返回了var，一般已经帮你算好了，但你要确认它是怎么算的。
显著性检验常用的是 North 准则：
[ \lambda_j^{error} = \lambda_j \sqrt{\frac{2}{N^*}} ]
其中 (N^*) 是有效自由度，通常用时间样本数除以某个因子估计。如果相邻两个模态的方差差异小于这个误差，说明这两个模态没有显著分离，解释时要谨慎。这个检验在REOF.m里不一定有，需要自己加。我一般会在旋转之后手动算一遍，把不显著的模态标出来，避免过度解释。
% 假设 lambda 是旋转后的解释方差（已排序），N 是时间样本数 N = size(X, 1); % 有效自由度估计，简单取 N/2，实际可用自相关算 Nstar = N / 2; % North 误差 lambda_err = lambda * sqrt(2 / Nstar); % 判断相邻模态是否显著分离 for i = 1:length(lambda)-1 if (lambda(i) - lambda(i+1)) < lambda_err(i) fprintf('模态 %d 和 %d 未显著分离\n', i, i+1); end end这段代码不是万能的，Nstar的估计有很多方法，但至少能给你一个粗略的判断。如果两个模态的解释方差差得还没误差大，那它们的顺序就没有意义，别硬说第一个比第二个重要。
4. 避坑与排查：REOF 实操中最容易翻车的五个地方
4.1 现象：旋转后模态空间型还是混在一起
原因：旋转模态数 (k) 取太小，或者原始 EOF 前几个模态本身就没分离好。旋转是在子空间里换基，如果子空间本身不够大，旋转也救不了。
解决：先把 (k) 调大，比如从 10 调到 20，看空间型是否变干净。如果还是混，回去检查 EOF 阶段，看看前几个模态的时间系数是不是高度相关，如果是，说明原始场里确实有耦合过程，REOF 也拆不开，这时候要考虑分区域做或者分季节做。
4.2 现象：解释方差排序后第一个模态只有 10%
原因：旋转把方差重新分配了，原本 EOF 第一个模态可能占 30%，旋转后分散到多个模态，每个都不大。这是正常的，不是 bug。
解决：不要拿旋转后的解释方差和旋转前的比。旋转后的解释方差之和等于前 (k) 个 EOF 模态的方差之和，但单个模态的贡献会变小。解释时看空间型的物理意义，不要只看方差数字。
4.3 现象：时间系数和原始场对不上
原因：旋转后的时间系数不是直接来自 SVD，而是用旋转后的空间载荷重新投影得到的。如果REOF.m返回的pc没有做这一步，或者投影时用了错误的载荷，时间系数就会错。
解决：手动验证。取一个旋转模态的空间载荷，和原始场做投影，得到时间系数，和REOF.m返回的pc对比。如果对不上，说明REOF.m的实现有问题，需要自己重写投影那一步。
% 手动投影验证 L1 = L(:, 1); % 第一个旋转模态的空间载荷 pc1_manual = X * L1; % 投影得到时间系数 % 和 REOF.m 返回的 pc(:,1) 对比 corr_coef = corr(pc1_manual, pc(:, 1)); fprintf('手动投影与返回时间系数的相关系数: %.4f\n', corr_coef);如果相关系数接近 1，说明一致；如果差很远，就要查REOF.m的投影逻辑。
4.4 现象：MATLAB 报错 "eof when reading a line" 或编码乱码
原因：这通常不是 REOF 本身的问题，而是数据文件读取时遇到了 EOF（文件结束）或者编码不匹配。MATLAB 2023 之后默认编码变了，老脚本里的中文注释可能乱码，导致解析出错。
解决：检查数据文件是否完整，用fopen时指定编码，比如fopen(filename, 'r', 'n', 'UTF-8')。如果是脚本注释乱码，把文件另存为 UTF-8 无 BOM 格式。这个坑和 REOF 无关，但会卡住整个流程。
4.5 现象：旋转结果每次跑都不一样
原因：Varimax 旋转是一个迭代优化过程，如果初始值随机，或者迭代次数不够，结果可能不稳定。另外，如果k很大，优化空间维度高，也容易陷入局部最优。
解决：固定随机种子，增加迭代次数。在REOF.m里找到旋转迭代的部分，把最大迭代次数从默认的 100 调到 500 或 1000。如果还是不稳定，说明 (k) 太大了，降下来。
5. 进阶：用 CCA 串起 REOF 模态与外部强迫，以及一个我常用的验证习惯
REOF 跑完之后，很多人会停在"画出空间型、写一段物理解释"这一步。但真正让分析站得住脚的，是把旋转模态和外部因子联系起来。标题热词里出现了 CCA（典型相关分析），这正好是 REOF 的下游工具。你可以把 REOF 得到的前几个时间系数作为一个场，把海温、风场或者某个指数作为另一个场，做 CCA，看哪个旋转模态和哪个外部因子耦合最紧。
具体做法是：取 REOF 前 5 到 8 个时间系数，组成矩阵PC_reof，取外部场做同样的去均值和标准化，组成Y，然后调用 MATLAB 的canoncorr：
% PC_reof: [nt, n_modes]，Y: [nt, n_vars] [ A, B, r, U, V, stats ] = canoncorr(PC_reof, Y); % A 和 B 是典型相关系数对应的权重，r 是典型相关系数 % U 和 V 是典型变量 % 看哪个模态在 A 里权重最大，就知道它和外部场关系最紧canoncorr是 MATLAB 自带的，不需要额外工具箱。A的每一列对应一个典型相关模态，看绝对值最大的那几个系数，对应到 REOF 的哪个模态，就能建立联系。r是典型相关系数，一般看前两三个，后面的通常不显著。
但这里有个坑：CCA 对样本量很敏感，如果时间样本少于 30 个，结果基本不可信。另外，CCA 之前一定要做标准化，否则量纲大的变量会主导结果。我一般会先把PC_reof和Y都做 z-score 标准化，再做 CCA。
验证 REOF 结果是否靠谱，我有一个习惯：把时间序列分成两段，前一半做 REOF，后一半做 REOF，看空间型是否相似。如果两段的空间型相关系数在 0.8 以上，说明结果稳定；如果差很多，说明这个模态不稳健，可能是噪声。这个做法比任何显著性检验都直观，虽然土，但管用。
还有一个技巧是，旋转后的空间载荷可以画成填色图，但要注意载荷的符号是任意的。同一个模态，整体乘 -1 还是同一个模态。所以比较不同实验的结果时，要先统一符号，否则会以为两个模态反了。我一般会选一个参考区域，如果这个区域的载荷是负的，就整体乘 -1，保证符号一致。
最后说一个我踩过的坑：REOF.m里如果用了svds而不是svd，对于大规模矩阵会快很多，但svds只算前几个奇异值，旋转需要的模态数如果超过了svds算出来的数量，就会出错。所以如果你要旋转 15 个模态，svds至少要算 15 个以上，最好多算 5 个作为缓冲。这个细节在脚本里往往不显眼，但一旦出错，报错信息很难懂。
希望这些步骤和坑能帮你把这份eof,reof等.rar里的东西真正跑起来，而不是停在解压完看一眼就搁置的状态。
本文还有配套的精品资源，点击获取
