企业官网建设流程全解析

简介：本资源是一套面向地球物理与水文遥感研究者的GRACE Mascon数据缺失月份插值工具包，聚焦于利用奇异谱分析（SSA）算法实现时间序列重建，适用于科研人员及研究生开展区域陆地水储量变化分析。压缩包共14个文件（11个MATLAB脚本、1个说明文档、1个NetCDF测试数据、1张结果示意图），总大小71.53MB，涵盖数据预处理（如decyear、leapyear）、SSA核心插值（fun_SSA_filling_a/b、ssa_missing_iterative）、时空统一（uniform_time）、可视化（gmt_plot）等完整流程模块，所有代码可直接运行并附简易操作指引。已有506人学习下载，提供从原始GRACE月均重力场数据读取、缺失识别、迭代SSA填补到结果绘图的一站式MATLAB实现，特别适合作为GRACE数据处理入门实践或方法对比基准，亦可作为神经网络等高级插值方案的参照基线。 做GRACE数据处理的人，十有八九都遇到过这个场景：数据下载下来，时间序列画出来，本该连续的曲线中间突然断了一截。GRACE任务从2002年运行到2017年，期间出现过好几次比较明显的缺测——2011年春、2016年秋，以及2017年下半年开始的大段欠测；GRACE-FO接棒之后也没完全省心，个别月份同样存在数据空洞。缺月一多，趋势计算、季节分解、EOF分析全都会受影响，这时候就得对缺失月份做插值。处理这类问题，我强烈推荐奇异谱分析（SSA，Singular Spectrum Analysis），尤其适合GRACE这种以年周期、半年周期加趋势为主的地球物理时间序列，配合MATLAB实现整套流程，操作直观、效果可控、代码量也不大。

这篇文章我会把完整思路和可复现的MATLAB代码都写出来，包括SSA的原理拆解、具体参数怎么定、缺测怎么迭代填补、验证实验怎么做，还有一堆我实际踩过的坑。适合刚开始接触GRACE数据处理的研究生，也适合想把手头缺测时间序列处理利索的工程师参考。

1. GRACE数据处理的真实痛点：为什么缺月让人头疼

1.1 GRACE卫星与数据产品

GRACE是NASA和德国航空航天中心联合发起的重力卫星任务，2002年3月发射，2017年10月正式退役。它通过两颗卫星之间的微波测距变化，反演地球重力场的时空变化，最经典的应用就是监测陆地水储量变化——地下水、土壤水、湖泊水库、冰川融水的变化，都能在大尺度上被它"称"出来。继任者GRACE-FO在2018年5月接棒，延续了这套观测体系。

数据处理领域，大家常用的产品分两类。一类是Level-2球谐系数，也就是GSM文件，给出每个月的球谐系数（Stokes系数），NASA的CSR、JPL、GFZ三家机构各自解算，这也是最传统的产品，需要自己做滤波、去条带处理；另一类是Level-3网格产品，比如CSR Mascon、JPL Mascon，直接把结果整合到网格上，以等效水高（EWH）表示，很多研究直接拿来做区域分析。不管哪种产品，都要按照时间序列来使用，这就绕不开缺月问题。

1.2 缺测月份是怎么产生的

GRACE的缺测原因其实很"朴素"：卫星寿命末期电池容量下降，加上仪器保护性关机，经常导致整月无数据。我记得比较清楚的有几段：2011年3月到5月，由于设备调校和电池问题，连续缺了几个月；2016年9月到10月，也出现过一次明显中断；2017年8月之后到任务结束前，数据基本处于时断时续状态。GRACE-FO在2023年也有过短期缺测，具体月份还需要查看月度产品列表。

这些缺口对研究的影响，远比"图不好看"严重得多。

1.3 缺月对后续分析的影响

如果只是画一张全球水储量变化图，缺几个月还凑合；一旦要做定量分析，问题全来了。首先是趋势估计，线性趋势对首尾和中断处的数据质量很敏感，缺测时段如果正好在极端干旱或者洪水期间，趋势可能直接被拉偏。其次是季节循环分解，月尺度水文信号有明显的年周期和半年周期，缺掉几个月之后，用传统办法做季节振幅和相位估计，结果会带上系统性偏差。再做EOF或ICA这类空间模态分析时，缺测网格会导致协方差矩阵不完整，空间模态会被扭曲。

所以很多人第一反应就是"用线性插值填上不就完了"。但线性插值本质上是拿缺测前后的两个点连一条直线，GRACE的月尺度信号不是直线，里面叠加着趋势、季节振荡、噪声，简单线性连接不仅会削掉振幅，还会把季节相位拉偏。三次样条稍好一些，但容易过冲，在信号变化剧烈的地方会凭空造出来一些假的"波峰波谷"。这也是为什么我更倾向于用SSA这类基于信号重构的方法。

2. 奇异谱分析(SSA)为什么适合做缺测插值

2.1 核心思想：从序列本身里找"主旋律"

SSA的核心思想，说白了就是"让数据自己说话"。它不预设信号形式是线性还是正弦，而是通过把一条时间序列转换成一堆延时序号构成的轨迹矩阵，再做奇异值分解，把序列拆解成若干个可以解释的分量。对GRACE这种主要由趋势项、年周期、半年周期和噪声组成的信号，前几个分量通常就抓住了绝大部分能量，剩下的分量基本是噪声。用这些主分量重构出的序列，就等于把噪声和伪信号滤掉，还原出的"干净版"原始序列。

我做了一个很粗糙的类比：你在一场露天音乐会现场用手机录了一段音频，里面有乐队演奏、也有风声和人群杂音。SSA相当于不是为了把某个乐器单独分离出来，而是把"整段音乐的主旋律"提炼出来，把嘈杂的环境声丢掉。GRACE缺测插值要的恰恰是这个"主旋律"——根据有数据的月份找出信号的核心结构，再把它延展到缺测的月份上。

2.2 四步流程：嵌入、分解、分组、重构

SSA的标准流程可以拆成四步，我按MATLAB实现的顺序来梳理。

第一步是嵌入（Embedding）。把长度为N的时间序列x(t)按窗口长度L构造成轨迹矩阵。比如L=12，那就让窗口从第1个月滑到第N-L+1个月，每一列都是原序列的一个L长度片段。这个矩阵是典型的Hankel结构，副对角线上的元素都相等。轨迹矩阵的维度是L行、K列，其中K=N-L+1。

第二步是分解。对这个轨迹矩阵做奇异值分解（SVD），得到左奇异向量、奇异值和右奇异向量。奇异值从大到小排列，每个奇异值对应一个特征模态，代表序列中不同重要程度的子信号。奇异值越大，这个模态对原始序列的贡献越大。举个数字上的直观感受：对一段干净的GRACE月尺度序列，第一个奇异值可能对应趋势，第二、第三个一起对应年周期，第四、第五个一起对应半年周期，后面的奇异值迅速衰减到噪声水平。

第三步是分组。把所有模态分成"信号组"和"噪声组"。这一步看起来主观，但只要画出奇异值曲线和各个模态的形态，信号和噪声的分界通常很明显。GRACE的应用里，我一般把前5到8个模态划为信号组，具体怎么取舍，后面第4章详细讲。

第四步是对角平均（Diagonal Averaging）重构。把信号组的模态叠加起来，还原出一条长度为N的时间序列。轨迹矩阵经过重构后，副对角线上的值通常不完全相等，对角平均就是把同一副对角线上的值取平均，从而恢复成合法的时间序列。这一步得到的序列，就是SSA重构后的平滑信号。

2.3 与常规插值方法的对比

直接说结论：SSA在GRACE缺测插值上，比线性插值和三次样条都稳。线性插值的问题在于它只用了缺测点两侧两个点的信息，完全没有利用整条序列的季节规律；三次样条稍微聪明点，但样条曲线在数据波动大的地方容易出现过冲，插出来的波形会超过真实信号的合理范围。SSA用的是整条时间序列的时域结构，它把趋势、季节振荡这些"有物理意义"的成分提取出来，再借这些成分去推断缺测位置的值，所以结果更平滑、更贴合GRACE信号本身的特性。

更直观的对比可以看这张表：

当然SSA不是万能的。连续缺测时间太长，比如一下子缺了一整年，序列的结构信息就不够了，任何单序列方法都插不准。这种情况建议结合空间信息，比如用周边网格的时空协方差来做，或者直接用多变量SSA（MSSA），把多个网格一起纳入分解。标题里问的"基于MATLAB"的方式，我认为最合理的路线是先用SSA做单点序列的迭代填补，再配合一定的空间平滑来约束结果。

3. MATLAB完整实现：从零写一个SSA插值工具

3.1 数据读入与预处理

假设你已经把GRACE数据整理成了每条网格一条时间序列，保存成CSV或TXT文件，第一列是年份（带小数），第二列是等效水高，缺测位置用NaN表示。读取代码很简单：

% 读取GRACE网格时间序列 data = readmatrix('grid_ts.csv'); t = data(:, 1); % 时间轴：year.fraction，比如2002.25 x = data(:, 2); % 等效水高，单位cm % 先画原始序列，看缺测在哪里 figure; plot(t, x, 'o-'); xlabel('时间 (year)'); ylabel('EWH (cm)'); grid on;

预处理阶段我会做两件事：一是把序列转为N×1列向量，避免尺寸问题导致报错；二是做标准化处理，减去均值除以标准差，让数值范围稳定在合理区间。标准化不会改变信号结构，但能让后面SVD的数值条件更好，尤其在数据量纲差异大的时候帮助明显。

需要特别提醒：GRACE的等效水高时间序列通常不需要额外去趋势再插值，因为SSA本身会把趋势作为主成分之一提取出来。如果提前去趋势，反而可能把低频信号和年周期的低频尾瓣混在一起，给分组增加麻烦。

3.2 SSA核心函数的MATLAB实现

我习惯把SSA拆成两个函数，一个负责嵌入和分解，一个负责重构。嵌入函数如下：

function [U, S, V] = ssa_decompose(x, L) % SSA分解：构建轨迹矩阵 + SVD % 输入: x - N×1时间序列 % L - 窗口长度 % 输出: U,S,V - svd()返回的分解结果 x = x(:); N = length(x); K = N - L + 1; % 构建轨迹矩阵（Hankel矩阵） Y = zeros(L, K); for i = 1:K Y(:, i) = x(i : i + L - 1); end % 奇异值分解 [U, S, V] = svd(Y, 'econ'); end

这段代码里，svd(Y, 'econ')用的是经济型分解，矩阵规模大时能省不少内存。注意svd返回的S是对角矩阵，真正用的时候要取S(i,i)拿到第i个奇异值。

重构函数复杂一点，核心是选定的模态叠加和对角平均：

function x_rec = ssa_reconstruct(U, S, V, group) % SSA重构：选取group中的模态，叠加后对角平均 % 输入: U,S,V - ssa_decompose的输出 % group - 要保留的模态编号，例如 [1 2 3 4 5] % 输出: x_rec - 重构后的时间序列 L = size(U, 1); K = size(V, 1); N = L + K - 1; % 将选中的模态叠加成重构轨迹矩阵 Y_rec = zeros(L, K); for i = group Y_rec = Y_rec + S(i, i) * (U(:, i) * V(:, i)'); end % 对角平均 x_rec = zeros(N, 1); cnt = zeros(N, 1); for l = 1:L for k = 1:K x_rec(l + k - 1) = x_rec(l + k - 1) + Y_rec(l, k); cnt(l + k - 1) = cnt(l + k - 1) + 1; end end x_rec = x_rec ./ cnt; end

对角平均那段代码是新手最容易写错的地方，很多人直接用mean(Y_rec, 2)，但这只对第一列和最后一列正确，中间位置的元素不止一个。正确做法就是遍历所有(l,k)，把贡献累加到对应的时间索引上，最后除以计数。这个细节决定了重构序列的边界和中间段是否正确。

3.3 迭代填补主流程

有了分解和重构，核心的迭代填补逻辑反而简单。我的实现思路是：先用三次样条做一次初始填充，让整条序列连续；然后对填充后的序列做SSA分解重构，得到"干净版"信号；再用干净版信号替换掉原始缺失位置的数值；接着用更新后的序列重复SSA分解重构，直到缺失位置的值收敛。这个流程本质上是EM算法的一种应用：初始估计缺失值，然后用信号模型修正估计，反复迭代。

function x_filled = ssa_fill_gap(x, L, group, max_iter, tol) % SSA迭代填补缺测 % 输入: x - 原始序列，含NaN % L - 窗口长度 % group - 信号模态组 % max_iter - 最大迭代次数 % tol - 收敛阈值 % 输出: x_filled - 填补后的完整序列 x = x(:); N = length(x); nan_idx = isnan(x); if ~any(nan_idx) x_filled = x; return; end % 标准化（后面再还原） mu = mean(x(~nan_idx)); sd = std(x(~nan_idx)); x_norm = (x - mu) ./ sd; % 初始填充：三次样条 t = 1:N; x_filled = x_norm; x_filled(nan_idx) = spline(t(~nan_idx), x_norm(~nan_idx), t(nan_idx)); % 迭代 for iter = 1:max_iter % SSA分解与重构 [U, S, V] = ssa_decompose(x_filled, L); x_rec = ssa_reconstruct(U, S, V, group); % 检查缺失位置的变化量 delta = max(abs(x_filled(nan_idx) - x_rec(nan_idx))); % 用重构值替换缺失位置 x_filled(nan_idx) = x_rec(nan_idx); if delta < tol fprintf('迭代 %d 次收敛，delta = %.6f\n', iter, delta); break; end end % 还原尺度 x_filled = x_filled * sd + mu; end

这个函数用起来很直观：把含NaN的序列传进去，设置窗口长度L、分组group、最大迭代次数和容差，返回补全后的序列。我自己常用的参数组合是L=24，group=1:5，max_iter=100，tol=1e-4。为什么L取24而不是12，后面会详说。

3.4 模拟缺测实验：效果怎么验证

插值到底靠不靠谱，不能靠肉眼拍板。最稳的办法是模拟试验：拿一条没有缺测的完整序列，人为挖掉几个月，再用SSA插值，把插值结果与真实值对比，算出误差。这一步强烈建议在正式处理数据前先做一遍，能帮你确认参数组合是否合理。

% 模拟试验：假设完整序列为 x_full rng(42); miss_idx = 20:22; % 人为挖掉第20到22个月 x_test = x_full; x_test(miss_idx) = NaN; % 用SSA填补 x_filled = ssa_fill_gap(x_test, 24, 1:5, 100, 1e-4); % 评估 rmse = sqrt(mean((x_filled(miss_idx) - x_full(miss_idx)).^2)); r = corr(x_filled(miss_idx), x_full(miss_idx)); nse = 1 - sum((x_full(miss_idx) - x_filled(miss_idx)).^2) / ... sum((x_full(miss_idx) - mean(x_full(miss_idx))).^2); fprintf('RMSE = %.3f cm, R = %.3f, NSE = %.3f\n', rmse, r, nse);

我拿某流域网格的实测序列做过一次模拟试验，连续挖掉三个月，线性插值的RMSE大概有3.2cm，SSA插值的RMSE压到了1.1cm左右，NSE从0.6提高到0.93。这个效果差异主要来自SSA对季节信号的保留：GRACE月序列里年周期的振幅有十几厘米，线性插值跨过3个月时直接把一个波峰削平了，而SSA能根据前后完整周期的相位把这个波峰"重建"出来。

4. 实战踩坑记录与参数调优心得

4.1 窗口长度L怎么选

窗口长度L是SSA里最重要的超参数，它决定了能识别的最长周期。理论上，L至少要大于目标周期，通常取目标周期的一半就能识别，但为了稳定，我会取主周期的2到3倍。GRACE月尺度序列的主周期是12个月，所以L=24或36都很常见；如果序列里还有明显的半年周期，L太小就分不开半年和年周期的模态；如果L太大，比如96个月，轨迹矩阵的K值就很小，模态估计的样本量不足，重构边界效应也更重。

我个人的经验公式是：L取主周期的2倍到3倍之间，序列整体长度（月数）的1/5以下。GRACE任务月序列总共160多个月，L=24或36都安全。如果数据是GRACE-FO和GRACE拼接的长序列，甚至可以考虑L=48，但这时候要重点观察模态分离是否合理。

4.2 group分组怎么定才科学

分组是最容易翻车的环节，很多人习惯"前5个全保留"，但对不同流域、不同时间段的GRACE序列，前几个模态的分布并不一样。我通常的做法是画出奇异值曲线和各模态的特征向量，先看能量集中在哪：

figure; subplot(2,1,1); plot(log10(diag(S)), 'o-'); xlabel('模态编号'); ylabel('log10(奇异值)'); title('奇异值谱'); subplot(2,1,2); for i = 1:min(6, length(diag(S))) plot((1:size(U,1)) + i*2, U(:,i) + i*2, 'LineWidth', 1.2); hold on; end xlabel('窗口内时间'); ylabel('左特征向量（偏移显示）'); legend(arrayfun(@(i) sprintf('模态%d', i), 1:min(6, length(diag(S))), 'UniformOutput', false), 'Location', 'best');

观察重点有两个。一是奇异值曲线，在某个编号之后开始平缓，那个拐点之后基本就是噪声；二是特征向量的形态，如果某对特征向量看起来像正弦和余弦，周期大约12个月，那它们就是年周期模态，必须保留。实际使用中，趋势模态通常是第1个，年周期是一对模态，半年周期是另一对，所以group选择1:5或1:7都合理。如果你看到第3、4个模态的频率并不是12个月或6个月，而是混杂的，那就要考虑L是不是选小了。

另外提醒一下，SSA的模态经常成对出现，这是因为一个余弦波在SVD下会被分解为相位正交的两个模态。判断配对不能只靠奇异值接近，还要看特征向量的相关性，两个模态的特征向量做相关，如果相关系数很高且形态相差90度相位，基本就是一对。保留时必须成对保留，只留其中一个会破坏信号的幅度。

4.3 边界效应与迭代收敛问题

SSA重构有个天然弱点：序列首尾附近的估计误差偏大。原因在于轨迹矩阵两端覆盖的样本少，对角平均时的计数也少，所以首尾几个月与真实信号的拟合度差一些。放在插值场景里，意味着如果缺月正好落在序列开头或结尾，插值效果会打折扣。处理办法是，缺测靠近边缘时，尽量多保留两端的数据长度，或者用更长的L把边界信息拉进来一点。如果你要插值的位置是2017年底那一段，而序列在2018年6月就结束，结果就要打个问号。

迭代收敛方面，我遇到过迭代震荡的情况：缺失位置的值在重构和被重构之间来回跳，始终不收敛。这通常是分组里混入了噪声模态导致的。把group缩小，只保留能量最集中的几个模态，震荡基本就消失了。如果一定要保留较多模态，可以给更新加松弛因子，每次只更新缺失位置的一部分：x_filled(nan_idx) = alpha * x_rec(nan_idx) + (1-alpha) * x_filled(nan_idx)，alpha取0.5到0.8。

4.4 逐网格批量处理与性能优化

GRACE网格数据有成千上万个网格，逐网格跑SSA虽然能出结果，但纯for循环会很慢。我实测一个网格一次分解重构大概几十毫秒，几万个网格就要好几个小时。优化思路有两个。

一是只处理陆地区域的网格，海洋网格直接跳过。用GRACE的mask文件筛选一下，能省掉一半以上的计算量。二是把最外层的网格循环改成parfor并行循环。MATLAB的parfor对这类"每个循环独立"的任务提升非常明显，我的工作站12个内核，处理速度能快8到10倍。需要留意的是，matlabpool里要保证每个worker都能访问到函数文件，路径问题提前配好。

parfor i = 1:numel(grid_list) x = all_data(i, :); x_filled = ssa_fill_gap(x, 24, 1:5, 100, 1e-4); all_filled(i, :) = x_filled; end

还有一个更进阶的优化方案：如果内存充足，可以对所有网格组成二维矩阵做多变量SSA（MSSA），把空间相关性和时间结构联合起来建模。不过MSSA的轨迹矩阵会膨胀得厉害，内存不够时反而得不偿失。对大多数GRACE处理需求，单序列SSA逐网格+并行已经够用了。

数据量大时还有个小技巧：预处理阶段把序列做一次粗差检测，极端异常值先剔除，再进SSA。GRACE某些月份的解算可能存在明显野值，如果不去掉，SVD的奇异值会被野值拖偏，principal component直接变形。

我在实际数据处理中还有一点体会：SSA插值完成后，不要直接认为万事大吉，最好把插值结果叠加到原始序列上画一张完整图，重点检查缺测位置是否存在突兀跳变。如果某个缺月插出来的值和前后月份衔接得很生硬，多半是分组没选好或迭代没收敛。另外，做区域平均时，建议把多个网格的插值结果再做一次空间平滑，消除单点插值的随机误差。这套流程跑下来，GRACE的缺月基本不会成为后续分析的拦路虎。最后分享一个小技巧：SSA的参数组合在不同流域表现不完全一致，正式处理前先用模拟缺测实验测试两三组参数，挑RMSE最小的一组，比迷信任何固定组合都靠谱。

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