企业官网建设流程全解析

简介：面向Matlab使用者的经验正交函数分析代码包，专门处理海表温度（SST）数据，并探索其与大气环流指数之间的滑动相关关系。压缩包内共两个文件，均为m脚本，整体大小仅3KB，包含数据读取与预处理、基于主成分或奇异值分解的经验正交函数计算、空间模态与时间系数输出、以及滑动相关分析等模块，结构紧凑，便于借鉴和按需二次开发。已有1048人学习下载，适合气候学、海洋学方向的学生及科研人员作为快速上手的参考范例，尤其适用于厄尔尼诺、北极涛动等典型气候模式分析。通过运行与研读代码，读者可以理解海表温度主要空间模态的提取流程，掌握经验正交函数与环流指数在连续时间窗口内的相关分析思路，为气候诊断与预测研究提供可用工具。

1. 为什么海温分析总从 EOF 开始：把 SST EOF 分析落到一包能跑的代码上

面对一份几十年的月平均海温格点场，最先想回答的两个问题几乎总是：哪些海域在同步地增暖或者变冷？哪个时间段里异常信号最强？EOFMatlab 这套代码给出的正是 SST EOF 分析的标准答案——把海温场做经验正交函数分解，得到一组按方差贡献从大到小排列的空间模态，以及每个模态对应的时间系数。做海洋再分析资料诊断、模式评估、海温指数定义的人，几乎都会把它当成第一工具。这套流程不挑资料：气候态海温场、业务化预报产品、同化系统输出都可以作为输入。EOF 原理本身不难，真正容易翻车的是预处理顺序、维度权重和结果判读，这些在后面的章节逐个展开。

2. 经验正交函数是怎么回事：协方差矩阵里的时空分解

2.1 从资料场到距平场：先减去平均，再做分解

EOF 全称经验正交函数，本质是把一个「空间点 × 时间点」的异常场，分解成少数几个正交的空间型态与时间序列。很多人第一次跑 EOF 时直接用原始 SST，结果第一模态几乎就是气候平均态：海洋温度高的地方载荷恒为正，时间序列几乎是一条直线。这不算做错，但说明没理解 EOF 是分析「异常」而不是「绝对状态」的工具。标准做法是先让每个格点的时间序列减掉自己的时间平均，再做分解。

下面这段代码演示的是数据排布方式，不是真实资料，只是让你先看清矩阵形状：

% 构造演示数据：16 个格点，48 个时次 rng(42); X = randn(16, 48); clim = mean(X, 2); % 每个格点的气候平均，16x1 Xa = X - clim; % 距平场，16x48 % 检查每行均值是否接近机器零 disp(mean(Xa, 2));

这里的约定是行代表空间格点（经度×纬度展平后的编号），列代表时间。之所以把空间展平，是因为协方差统计关心的是「格点之间的同步关系」。mean 沿第 2 维度做，得到每个格点自己的平均值，减法之后每一行就不再携带该点的平均态信息。如果资料是逐月多年覆盖，这里减的是总平均；到第 4 章会把更精细的季节循环也扣掉。

2.2 特征值、空间模态与时间系数：三个结果对应什么

距平场 Xa 的协方差矩阵这样定义：C = Xa * Xa' / (nt-1)。C 是一个「空间格点数 × 空间格点数」的方阵，第 i 行第 j 列描述了格点 i 与格点 j 的异常在海温起伏上是否同步。对 C 做特征分解，特征向量就是我们要的 EOF 空间模态；把 Xa 投影到这些特征向量上，得到的系数序列就是时间系数（PC）。特征值则告诉你这个模态在整个场里解释了多少异常方差。

三者对应关系和常用说法，列成一张表：

一个海温场里，第一模态往往与大尺度海温异常型对应，因为它覆盖区域广、时间持续性又强，对总方差的贡献通常最大；第二模态则可能是年代际的南北不对称型或中纬度带状模态，要看具体资料时段。PC 序列里的高值年份，就是该空间型表现最强的年份。判断一个模态是否可信，先看它的方差占比，再看 PC 与已知物理信号的对应关系。

2.3 为什么选择 SVD 而不是直接解协方差特征向量

直接从 C 出发做 eig 在数学上等价，但在实际数据上我不推荐。第一次算 EOF 时我见过某开发者直接写[V, D] = eig(Xa * Xa')，小区域资料没问题，换成全球 1 度格点海温后，C 的维度是 64800×64800，矩阵乘法就要数 GB 内存。改用 SVD，Xa = U * S * V'，左奇异向量 U 就是空间模态，奇异值平方再除以自由度就是对应方差，不需要显式构造协方差矩阵：

[U, S, V] = svd(Xa, 'econ'); % Xa 是空间x时间 EOF = U; % 左奇异向量，列是空间模态 PC = V * S; % 时间系数，行对应时间 var_explained = diag(S).^2 / (size(Xa,2)-1);

'econ' 选项的意思是只计算真正需要的列或行，避免把巨大矩阵全部算出来。注意 PC 的写法：左乘 S 是为了让时间系数保持与原场相同的量纲，否则奇异值全压在 EOF 上，画图时数值会显得很别扭。后面第 3 章会给出一个更省内存的版本，直接用 Xa 与 V 反推 EOF，避免把 U 完整存下来，这也是处理全球海温场时的关键优化。

提示：直接解协方差矩阵在数学上没问题，但在全球格点资料上内存和时间成本都很高。SVD 是经验正交函数的标准实现路径。

3. 把 EOFMatlab 核心函数拆开看：输入输出与一个可跑的接口

3.1 函数接口设计：一个函数搞定权重、中心化与截断

这套代码的核心就一个函数：eof_sst.m。它把权重、中心化、SVD、符号对齐、重建集成在一起，调用者只需要把「空间×时间」的距平矩阵和纬度权重向量传进去。先看完整实现，再讲参数怎么改：

function [EOF, PC, variance, recon] = eof_sst(X, w, K) % EOF_SST 经验正交函数分解 % 输入: % X : 空间格点 x 时间 的距平场（可选，内部会再中心化） % w : 空间格点 x 1 的权重向量，如 sqrt(cosd(lat)) % K : 保留模态数，默认 10 % 输出: % EOF : 空间 x K 空间模态 % PC : 时间 x K 时间系数 % variance : K x 1 每个模态解释方差百分比 % recon : 空间 x 时间 前 K 模态重建场 if nargin < 2 || isempty(w), w = ones(size(X,1), 1); end if nargin < 3 || isempty(K), K = 10; end % 1. 中心化：每行减时间平均 Xa = X - mean(X, 2); % 2. 纬度加权：权重广播到每一列 Xw = Xa .* w; % 3. 奇异值分解：只用 V 和 S，省掉 U 的存储 [~, S, V] = svd(Xw, 'econ'); % 取前 K 个奇异值 S_k = S(1:K, 1:K); V_k = V(:, 1:K); % 4. 空间模态 U_k = Xw * V_k * inv(S_k) EOF = Xw * V_k / S_k; % 5. 时间系数 PC = V_k * S_k; % 6. 方差贡献百分比 var_total = sum(diag(S).^2) / (size(X,2) - 1); variance = diag(S_k).^2 / (size(X,2) - 1) ./ var_total * 100; % 7. 符号对齐：让每个模态的最大绝对载荷为正 for k = 1:K [~, idx] = max(abs(EOF(:,k))); if EOF(idx, k) < 0 EOF(:,k) = -EOF(:,k); PC(:,k) = -PC(:,k); end end % 8. 用前 K 模态重建 recon = EOF * PC'; % 还原权重缩放 recon = recon ./ w; end

这段代码把第 2 章的两个数学结论直接焊进了接口里。第 3 步svd返回~代表丢弃左奇异矩阵 U；没有 U，后续用公式EOF = Xw * V_k / S_k反推左奇异向量。这里用的是矩阵右除，数学含义是乘上 S_k 的逆。为什么要多此一举？全球 1 度海温资料里，空间格点约六万，时间序列几百，U 如果是 64800×几百，每个元素 8 字节，就是上百 MB 的显式存储。换掉之后内存占用直线下降，这是处理大数据场时关键的一笔。

第 7 步符号对齐值得单独说一句。特征向量的符号本身没有物理意义，同一个模态可以整体乘 -1，画出来完全是翻转的冷暖分布。如果不做约束，同一台机器跑十次可能得到冷暖区互换的图。代码里的策略是找到每个模态中绝对值最大的格点，强制这个格点的载荷为正，这样绝大多数情况下空间型正负能够视觉上保持一致。

3.2 参数说明与调用示例：K 取多少、权重怎么传

函数三入参 X、w、K，行为各有边界：

调用层面，建议先打印方差贡献，再决定到底保留几个模态：

load sst_anom.mat; % 假设里面是 Xa 和 lat_vec w = sqrt(cosd(lat_vec)); % 一维纬度权重 % 注意 Xa 的格点行顺序里每一行对应一个纬度，w 要广播到每行 if size(Xa, 1) == length(lat_vec) w_all = w; end [EOF, PC, variance] = eof_sst(Xa, w_all, 8); for k = 1:8 fprintf('EOF%02d : %.2f%%\n', k, variance(k)); end

上面这段假定距平场已经整理成行，但实际资料往往是经纬网格，所以第 4 章会给出完整的从网格到矩阵的 reshape 流程。K 取值没有绝对标准，我的习惯是先放 8 到 12，然后看累计方差贡献，前几个模态能到 70% 以上就可以收了；K 越大并不能说明进步，超过资料时间长度后还会出现伪模态。

3.3 结果后处理：模态顺序、方差贡献与符号约定

SVD 返回的顺序天然按奇异值从大到小排列，所以 EOF 第一列是第一模态，第二列是第二模态，PC 的列序与 EOF 一致。这里最容易把人绕晕的是 PC 的行列方向。函数里 PC 是「时间×K」，也就是说PC(:,1)画出来就是第一模态的时间曲线。如果你换一种实现方式把 PC 写成「K×时间」，出图前要先转置，这样的事故我见过不止一次。

另一个对新手极有帮助的验证习惯：用第一模态的 PC，与目标海域平均海温异常序列做相关。比如你关心热带诊断海域的信号，把该区域的面积加权平均距平序列取出来，corrcoef(PC(:,1), 区域序列) 常常能到 0.8 以上。如果相关系数低得离谱，大概率不是模态错，而是前面预处理里哪一步把符号或者相位弄反了。

4. SST EOF 分析全流程实战：从月平均资料到前三个模态

4.1 读入海温资料与提取海洋格点

我们以一份常见的月平均海温 NetCDF 文件为例，变量一般是经度 lon、纬度 lat、时间 time、海温 sst，陆地点为 NaN。第一步不是急着算 EOF，而是把海洋格点挑出来，陆地格点在后续所有步骤里都属于噪音源：

% 读取 NetCDF 变量 lon = ncread('sst_mon.nc', 'lon'); lat = ncread('sst_mon.nc', 'lat'); sst = ncread('sst_mon.nc', 'sst'); % lon x lat x time time = ncread('sst_mon.nc', 'time'); % 时间维长度 nt = size(sst, 3); nx = length(lon); ny = length(lat); % 构造经纬网格，方便之后 reshape [LON, LAT] = meshgrid(lon, lat); % 用多年平均识别稳定海洋格点 clim_mean = mean(sst, 3, 'omitnan'); ocean = ~isnan(clim_mean); % 提取海洋格点下标 idx_ocean = find(ocean); nspace = length(idx_ocean);

我特意用气候平均而非第一个时次来判断海洋格点，原因在于海冰边缘：单一时次里极区会出现大片 NaN，但用多年平均会把偶发的冰区缺测也纳入筛选。取到 ocean 逻辑矩阵后，后续所有距离计算和绘图都以它为底图。注意经度排序：如果资料里经度是 0 到 360，后续画图时要提前想好区域范围。

4.2 扣除季节循环、处理趋势与 NaN

逐月资料必须先做季节性扣减。海温的季节循环强度远大于年际信号，不扣掉它，第一模态会被季风性增暖这类纯气候态信号占据。做法是逐月算气候态，再按月份原位相减：

% 气候态：每个月份 30 年平均 clim = zeros(nx, ny, 12); for m = 1:12 clim(:,:,m) = mean(sst(:,:,m:12:nt), 3, 'omitnan'); end % 月距平 sst_anom = sst; for m = 1:12 sst_anom(:,:,m:12:nt) = sst(:,:,m:12:nt) - clim(:,:,m); end

这一步之后，每个格点的时间序列已经扣掉了多年平均和季节循环。要不要再扣线性趋势，看你的研究目的：诊断海温模态与气候态变化时，保留趋势能揭示长期增暖信号；检测年际变率时，趋势往往会污染前两个模态，建议用 polyfit 把每个格点的线性趋势去掉。我的做法是默认保留，但在计算年际异常指数时另存一份去趋势版本。

上面代码段里sst_anom(:,:,m:12:nt)这种索引方式假设第 1 时次是 1 月，如果资料从 2 月开始，先对齐月份再执行。为了省事，也可以用 reshape 把时间维重排成 ntime/12×12，再按列直接减对应月份气候态，两种办法结果一致。

NaN 的处理原则是宁缺勿填。陆地点、冰区点在这一步已经通过 ocean mask 排除；海温资料偶尔在个别时次有缺测，常见做法是只保留所有时次都完整的格点。插值填补极区缺测会让 EOF 出现沿插值边界的条带状伪模态，我不建议对海温场做全局插值。

4.3 构造加权数据矩阵并运行 EOF

全球 1 度格点下，高纬 1 格点代表的面积远小于低纬。海冰区变异又大，如果不做纬度加权，EOF 前几个模态会被极区抢走，热带信号被严重压制。标准做法是乘 sqrt(cosd(lat)) 面积权重，然后恢复成统一的矩阵：

% 面积权重：纬度余弦开方 w_lat = sqrt(cosd(lat)); % ny x 1 W = repmat(w_lat, 1, nx); % 转成 ny x nx，对齐 LAT % 给所有海洋格点取权重 w_vec = W(idx_ocean); % nspace x 1 % 把海洋格点距平场重排成 空间 x 时间 Xa = reshape(sst_anom, nx*ny, nt); Xa = Xa(idx_ocean, :); % 只留海洋格点 % 检查每行是否干净 assert(all(isfinite(Xa(:))), 'Xa 里还有 NaN，检查 ocean mask'); % 跑 EOF [EOF, PC, variance] = eof_sst(Xa, w_vec, 6); % 打印累计方差 cumvar = cumsum(variance); fprintf('前4个模态累计方差: %.1f%%\n', cumvar(4));

W 的构造这里有个不仔细看就会搞错的细节：LON/LAT 来自 meshgrid，LAT 的每一行都是相同的纬度；直接取W(idx_ocean)前要确认 idx_ocean 是按哪个维度顺序 find 出来的。我这里统一用LAT(:)的维度顺序，所以权重矩阵也要先repmat成 ny×nx。另一种更稳妥的方式是把经度纬度和权重一起整理到向量后，再做一次 find，就不会错位。如果输出显示累计方差过低，先检查这一行。

eof_sst内部已经做过中心化和加权，所以传入的 Xa 不必预先再减总平均。但距离标准做法还有一步：第 4.2 节的季节循环处理已经把原始场变成了距平场，时间平均减掉后剩下的是季节内和年际异常，层次是分开的。

4.4 出图与物理解读：第一模态、海温异常型与时间序列

得到 EOF 后，先把它还原成经纬网格再画图。这里针对海洋格点的矩阵，陆地部分填 NaN：

% 把第一模态放回网格 eof_map = nan(nx, ny); eof_map(idx_ocean) = EOF(:,1); eof_map = eof_map'; % 画图：经度为 X 轴，纬度为 Y 轴 figure; contourf(LON, LAT, eof_map, 21, 'LineStyle', 'none'); colorbar; title('EOF1 空间模态'); % 用你常用的海岸线绘图工具箱叠加边界 % 时间系数 figure; plot(time, PC(:,1), 'LineWidth', 1.2); ylabel('PC1');

解读时我的顺序是：先看第一模态的方差占比，大于 15% 说明场里有清晰的准周期信号；再看高载荷区落在哪个海域，结合 PC 的峰值年份判断它对应什么海温异常型。如果第一模态空间型分裂成两个不相关的极值区、PC 又出现明显的半年周期，多半是季节循环没扣干净，回到 4.2 重新做。

一个实操技巧：把 PC 标准化之后与区域平均海温异常序列画在同一张图里，尺度对齐后可以直观看出位相是否一致。若两者峰值差几个月，通常是资料时间轴定义问题，检查 NetCDF 的 time 变量的参考日期，不要默认它是公历 1 月。

注意：EOF 空间载荷数值是相对大小，不是温度异常本身。某个格点载荷 0.8 不代表该处有 0.8℃ 异常，只代表它和该模态的同步程度高。

5. SST EOF 避坑笔记：五个很容易翻车的地方

EOF 在统计上是成熟方法，但实跑时坑几乎都埋在数据侧。以下五条都是实际跑 SST 时反复出现的血泪经验，按「现象 → 原因 → 解决」格式整理。

5.1 模态符号随机翻转，同一份资料两次结果冷暖区互换

现象：同一个脚本跑两次，第一模态的暖中心变成冷中心，PC 序列也整体反转。 原因：特征向量方向不唯一，整体乘负号仍是合法解。SVD 得到的 U 列没有固定正负规则。 解决：按第 3 章函数里的方式做符号对齐。更稳妥的进阶做法是选定一个物理含义明确的参考格点或区域，计算 PC 与该区域平均序列的相关系数，相关系数为负时把整个模态与 PC 同时取负。后者能同时对齐符号与物理意义。

5.2 相邻模态排序不稳定：特征值太接近，第二第三模态互换

现象：资料长度微调几个月，第二模态和第三模态的出现顺序就换位，空间型也随之大变。 原因：EOF 的模态顺序按特征值大小排列，当两个特征值落入彼此的置信区间时，模态方向在数学上可以旋转混合，排序和形状都缺乏稳健性。 解决：引用经典的特征值不确定度判据，当相邻特征值间隔超过置信区间时才认为模态分离可靠。具体做法是计算λ ± λ * sqrt(2/nt)的区间，若两区间重叠，不强行解释这两个模态的独立物理意义，改用旋转 EOF 或只取第一模态。

5.3 高纬海冰区抢戏：没做纬向面积加权

现象：第一模态空间型主载荷出现在南北极海冰边缘，热带信号微弱，PC 高频抖动明显。 原因：高纬一个格点的物理面积远小于低纬，海冰区的年际变率又大，未加权的协方差统计被高纬主导。 解决：在构造数据矩阵时乘sqrt(cosd(lat))。函数内置权重入口，只要外部传入正确的w_vec就不会踩。检查方法是打印每个模态的方差占比，如果第一模态超过 30% 且高纬载荷突出，多半是权重没传或者传成了全 1。

5.4 陆地点被填成 0，空间型出现方形假象

现象：EOF1 的空间型里大陆形状清晰可见，海岸线附近出现方框状高载荷带。 原因：某些代码在处理 NaN 时习惯性sst(isnan(sst)) = 0。陆地被当成 0℃ 的稳定值参与协方差统计，与海洋距平完全不同步，于是作为「额外模态」被分解出来。 解决：坚持用海洋 mask 提取格点，保证进入 eof_sst 的矩阵里没有 NaN，也不需要填 0。顺手用第 4 节里的assert(all(isfinite(Xa(:))))做入口校验，一劳永逸。

5.5 第一模态竟然是完整的季节循环

现象：PC1 周期为一年，EOF1 看起来就是全球冬季海温距平分布，方差贡献甚至超过 50%。 原因：季节循环没扣干净。只减了年总平均没按月扣气候态，或者扣气候态时月份索引错位。 解决：回到 4.2 节，确认每个时次的月份索引与气候态一一对应，并用谱分析检查 PC 曲线在 12 个月周期上是否还有显著峰值。常见做法是用fft快速检查 PC 主周期，如果 12 个月周期振幅仍然明显，说明季节扣减有残留。

6. 再进一步：模态重建验证与滑动窗口 EOF

6.1 用前 K 个模态重建原场，验证分解完整性

eof_sst 返回的 recon 本身就是一个现成的验证工具：把前 K 个模态乘回时间系数，再还原权重，就能得到重建场。对每个格点计算重建场与原始距平场的相关系数，接近 1 表示前 K 个模态已涵盖绝大部分信号；明显偏低，则说明场里还存在较强的局地变率。

% 重建前 4 个模态 [EOF, PC, variance, recon] = eof_sst(Xa, w_vec, 4); % 逐个格点算相关系数 rho = zeros(size(Xa,1), 1); for i = 1:size(Xa,1) tmp = corrcoef(Xa(i,:), recon(i,:)); rho(i) = tmp(1,2); end % 看低相关格点分布在哪 bad = find(rho < 0.8); fprintf('低相关格点数: %d\n', length(bad));

我习惯把 rho 低值画回经纬网格。如果它们集中在中纬度西边界流区，属于正常现象——强流区局地变率太强，模态分解本来就不容易完全覆盖；但如果大面积出现低相关，就要回头查资料质量。

6.2 滑动窗口 EOF：看模态解释方差如何随时间演变

固定全时段 EOF 假设模态空间型不随时间改变。现实里海温主模态的强度有年代际变化，更进阶的做法是滑动窗口 EOF，每次只取一个窗口长度的数据做分解，记录窗口中心年份的方差贡献与空间型变化：

win = 360; % 窗口长度，360 个月约 30 年 step = 12; % 每年滑动 12 个月 nwin = floor((nt - win)/step) + 1; var1 = zeros(nwin, 1); yr_center = zeros(nwin, 1); for i = 1:nwin idx = (i-1)*step + (1:win); Xw = Xa(:, idx); [~, S, ~] = svd(Xw, 'econ'); var1(i) = S(1,1)^2 / sum(diag(S).^2) * 100; yr_center(i) = time(floor(idx(1) + win/2)); end plot(yr_center, var1); xlabel('窗口中心年份'); ylabel('EOF1 方差贡献(%)');

窗口长度的选择有讲究：太短，样本数不足会让 EOF 不稳定；太长，又回到全时段静态结果。30 年窗口用于年际到年代际信号是常见的折中，步长取 12 个月能保证窗口间有较大重叠，曲线平滑。窗口内如果出现强海温异常事件，第一模态方差贡献会显著抬升，这正是滑动 EOF 想捕捉的信息。

从那以后，我每次跑 SST EOF 都强制走一遍「海洋 mask → 季节循环 → 纬度权重 → 符号对齐 → 重建相关性核对」五步，前两步省事、后三步防坑。这套流程里每一步都在前面对应了一个具体代码块，直接用即可。希望帮到你。

本文还有配套的精品资源，点击获取