GRACE卫星重力数据缺失月份插值：基于奇异谱分析(SSA)的MATLAB实现

简介：本资源是一套面向地球物理与水文研究者的GRACE卫星Mascon数据缺失月份插值工具包，聚焦于利用奇异谱分析（SSA）算法实现时间序列重建，适用于缺乏多源辅助数据、难以构建神经网络模型的中小尺度研究场景。压缩包共14个文件，含11个MATLAB核心脚本（如ssa_missing_iterative.m、fun_SSA_filling_b.m等，覆盖数据预处理、SSA分解重构、时序对齐与可视化全流程）、1个NetCDF格式实测GRACE Mascon数据、1个说明文档及1张结果示意图，整体大小71.53MB。已有506人学习下载，资源提供完整可运行代码、内置测试数据及简易操作指引，开箱即用；程序模块划分清晰，包含闰年处理（leapyear.m）、十进制年转换（decyear.m）、地理坐标提取（get_coord_n.m）等实用工具函数，显著降低GRACE数据处理门槛，尤其适合初涉重力场反演与长期水储量变化分析的科研人员快速上手。 做GRACE卫星重力数据处理，尤其是用等效水高（EWH）时间序列研究区域水储量变化时，我最常被问的问题不是怎么滤波，而是：数据缺月份了怎么办。GRACE数据从2002年持续到2017年，后面GRACE-FO又接上，但期间因为卫星轨道维持、电池老化、日食季供电不足这些原因，月度产品断档非常普遍。你去下载任何一家机构的月度网格，拿到的序列一定不是等间隔的。时间序列分析、季节性提取、趋势显著性检验，全都要求数据连续，所以缺失月份插值就成了绕不开的一步。而插值方法里，奇异谱分析（SSA）又是一个被低估的好工具——它不仅能补空缺，还能把趋势、季节信号和噪声拆开，非常契合GRACE这类以低频信号为主的序列。这篇博文我会把整个流程讲透，包括为什么用SSA、窗口长度怎么选、MATLAB代码怎么写、插完怎么验证，把我实际踩过的坑也一并交代。

1. GRACE月度序列的常见缺口：先从数据源头说起

1.1 GRACE数据长什么样，为什么偏偏缺月份

GRACE卫星通过精确测量两颗卫星之间的微波测距变化，反演全球重力场的时间变化。对做水文、大地测量的人来说，最常见的产品就是月度全球重力场网格，每个格点的数值代表该月相对于多年平均的等效水高异常（Equivalent Water Height, EWH），单位通常给到cm，有些产品给mm。空间分辨率一般是1°×1°或者0.5°×0.5°，时间分辨率就是一个月一个值。所以如果你研究某个流域、某个含水层，提取出来的原始时间序列，本质上就是200多个月的高斯滤波后水储量异常。

但真实数据不会是完美连续的。我整理过几套数据，缺失月份多的时候能占到总长度的20%以上。GRACE出现缺失的原因主要集中在几个方面：

• 日食季（eclipse season）期间卫星长期处于地球阴影区，太阳能供电不足以支持全部载荷连续工作，测量任务被迫降级或中断。
• 电池性能随任务年限加深不断衰减，2016年以后这个问题尤其严重，导致大量月份没有正式发布产品。
• 卫星姿控、星间测距系统在部分时段进入校准或安全模式，数据段长度不足以生成可靠的月度解。
• 不同数据处理中心对质量较差的月份处理策略不一致，有些机构直接标记为缺失，有些给出但建议不要用。

我在实际处理中还会遇到一个更隐蔽的事情：数据文件明明存在，但某些格点内的数值被设为特殊填充值（比如-9999或99999），这本质上也是缺失。所以拿到数据后的第一步永远是统一质量控制，而不是直接插值。把质量标记和填充值全部转成NaN，再做后续处理，否则插值器会把-9999这种值当成真实信号，结果全乱。

1.2 缺失数据带来的连锁麻烦

缺失月份不处理，最直接的问题就是用不了现成的时间序列工具。MATLAB里很多函数比如fft、filter、movmean、arima，尽管有的能容忍NaN，但输出会变得不可靠，自回归类的模型基本不工作。

更麻烦的是偏差问题。GRACE序列本身有很强的季节循环，如果缺口不随机，而是系统性集中在某些时段，比如2016年之后连续缺了十几个月份，那直接求年平均会导致这年冬季权重偏低、夏季权重偏高，算出来的年际变化就是歪的。做趋势估计的时候，缺失位置不同，最小二乘拟合出的斜率可以差出30%甚至更多，这个我在实验里反复验证过。所以插值不是简单的"不美观"问题，而是直接关系到后续所有定量结论的正确性。

对GRACE这种长序列低频信号，理想插值需要做到两件事：一是恢复被缺失时段覆盖掉的趋势和季节信号，二是不要引入过多人为的高频振荡。普通的局部插值方法很难同时满足这两点，这就是SSA这类全局重构方法的价值所在。

2. 缺失月份插值的三条路线与选型逻辑

2.1 最直接的路线：通用插值函数

MATLAB里有interp1、fillmissing这些现成函数，用起来几行代码就搞定，也是大多数人的第一选择。线性插值假设相邻值之间是线性过渡，适合缺口很短、序列本身波动不剧烈的情况。但我实测下来，GRACE月值序列在一年内的变化幅度通常有10-20cm等效水高，如果连续缺3个月以上，线性插值会系统性地削平峰谷，把季节振幅明显低估。

PCHIP（分段三次Hermite插值）比线性好一些，它不会在极值附近产生过冲，适合有周期性变化的序列，能够比较好地保留缺口的凹槽或峰顶形态。样条插值（spline）看起来最平滑，但很危险，因为它要求二阶导数连续，遇到较长缺口会在端点附近出现明显过冲（overshoot），插出一些物理上不可能的负值或巨大峰值。对GRACE水储量这种物理量，样条插出的离谱值一旦进入后续分析，往往要花很长时间才能发现。

所以如果只是零星缺一两个月，我推荐用fillmissing(ts, 'pchip')，速度快、无参数、结果稳健。但如果你连续缺了3个月以上，或者整个任务末期缺了大半年，通用插值基本不够用，需要用全局方法。

2.2 更聪明的路线：基于SSA的迭代重建

SSA插值的思路和局部插值完全不同。它先假设真实信号可以表示为少数几个主要成分的叠加（趋势+逐年周期+半年周期+少量噪声），然后通过矩阵分解把时间序列拆成这些成分，再用主要成分去估计缺失位置的值。

具体到缺失填充，比较成熟的做法是迭代式填充，流程很简单：

1. 先用线性或PCHIP给缺失点补一个初始值，得到完整序列。
2. 对完整序列做SSA分解，得到前d个主要成分的重构序列。
3. 用重构序列去更新缺失点的旧值。
4. 重复进行SSA分解和更新，直到缺失点数值变化小于阈值。

这个流程执行下来，缺失值会逐步收敛到一个与整体信号结构自洽的状态。因为SSA重建过程中，非缺失点的信息也参与了每个主成分的估计，相当于用全序列的数据来约束缺失段，比只用缺口边界局部的信息去猜要稳健得多。我后面会给完整可跑的MATLAB代码。

2.3 进阶路线：时空联合插值（DINEOF等）

如果研究区域不是一个格点，而是整个网格场，还可以用DINEOF（Data Interpolating Empirical Orthogonal Functions）这类方法。它的原理是利用时空场的经验正交函数结构来填充数据，本质上可以看成SSA在高维空间里的扩展。DINEOF在海洋遥感数据（SST、叶绿素）里用得很多，插值同时能保持空间连续性，对大面积云覆盖导致的缺测特别好用。

但对GRACE月度网格来说，有一个现实问题：如果是某个月份整个全球数据都缺失，DINEOF也一样无从下手，因为它需要利用已有的时间协方差结构。而且DINEOF的调参过程（判断截断模态数、迭代收敛标准）比SSA更复杂，对刚接触数据处理的同学没那么友好。我通常的建议是：先用SSA把单点时间序列补好，如果后续要做空间场分析，再考虑DINEOF这种更重的工具。

2.4 到底怎么选：一张选型表

我把各种情况下的推荐方案总结如下。

表格之外还有一点需要强调：无论选哪种插值，最终都要做交叉验证，把已知的完整月份故意挖去，用剩余数据重建，再与真值比较。不做验证的插值，参数随便一调都会看起来很漂亮，但真实误差你心里没底。

3. 奇异谱分析(SSA)的原理与参数选法

3.1 SSA的一次完整计算过程

SSA的核心思想，是找一个合适的窗口把一维时间序列转成一个"轨迹矩阵"，然后对这个矩阵做奇异值分解，再用特征值较大的那几个分量重建原序列。

我用一个日常类比来帮助理解。假设你有一盘混合音频录音，里面有人的说话声、背景音乐和电流底噪。SSA做的事情就是把这盘录音切成许多等长的片段窗口，然后分析这些片段中重复出现的模式——说话声和音乐是有规律的，底噪是随机的。通过矩阵分解，能把这三种成分在数学上分离开来。GRACE序列也是类似的混合体：长期趋势（如地下水持续减少）、年周期（降水补给和蒸散发的季节变化）、半年周期、随机噪声，它们各自的"规律性"不同，SSA就按这种规律性强弱把它们排序。

数学过程可以拆成四步：

第一步，嵌入。给定序列x(t)，长度N，选择窗口长度L（1 < L < N），构造轨迹矩阵X。矩阵的每一列是长度为L的滑动窗口片段，总共有K = N - L + 1列，所以X是L×K的矩阵：

X = [x(1) x(2) ... x(K) x(2) x(3) ... x(K+1) ... x(L) x(L+1) ... x(N)]

第二步，奇异值分解。对X做SVD分解：X = UΣV^T。U是L×L正交矩阵，V是K×K正交矩阵，Σ对角线上的元素就是奇异值σ，按从大到小排列。奇异值平方代表该成分对总方差的贡献。

第三步，分组。把奇异值/奇异向量按贡献大小分为若干组。通常第一大奇异值对应趋势项，之后每两个奇异值一组对应一个周期的正弦/余弦对（比如年周期、半年周期），剩下的都是噪声项。只保留前d个主成分，其余置零。

第四步，重构。把保留下来的分量合回去，得到去噪后的矩阵，再沿反对角线做平均（对角平均），重新变回一条与原始序列等长的时间序列。对角平均是必需的，不然无法把矩阵还原成一维序列。

整条链路里，没有假设序列是线性的，也没有假设窗口内的局部关系，所以它对非平稳、周期成分复杂的序列比简单插值更合适。GRACE序列恰恰满足这种可分解性，这是SSA能在这类问题上发挥作用的前提。

3.2 窗口长度L怎么定

窗口长度L是SSA最重要的参数，直接决定分解结果，因为L本质上是你要捕捉的波动的最长时间尺度。如果L太小，比如小于12，就无法在窗口内完整容纳一个年周期的形态，年周期会被撕裂到多个分量里；如果L太大，轨迹矩阵行数太多，分解出的分量数量膨胀，模态混叠和边缘效应都会变严重。而且L取太大后，每个分量对应的窗口片段太少，统计意义下降，重构也会退化。

我的经验是分两种场景处理。对于GRACE月度序列，一个自然年周期是12个月，如果要捕捉这个周期，L至少要大于12，最好取24或36，这样窗口内能包含两到三个完整年周期，周期信号在奇异值分解中可以更稳定地配对。对于长度只有150-200个月的GRACE序列，L=36已经占序列长度的五分之一到四分之一，用它做滞后矩阵的列数还有一百多，足够完成SVD。如果序列更长（比如2002到2024年积累到260个月），L=48也是可以考虑的，但收益不明显，反而让重构端点误差更宽。

还有个实用技巧：如果只想捕捉年周期而忽略半年周期，L取24就够；如果还想把半年周期也稳定分离，L≥36更合适。因为窗口长度必须大于最大目标周期的两倍，才能保证基频和谐波不在SVD里发生严重混叠。这是我在对比很多组实验后得出的结论。

3.3 如何识别信号分量：特征值谱与累计方差贡献率

选定L后，怎么决定保留几个主成分（也就是d的取值）？最直观的工具是特征值谱图。把奇异值平方（等价于特征值）按降序画成折线图，然后看拐点在哪里：特征值从某个位置开始变得很小且下降平缓，这个位置就是信号和噪声的分界。GRACE序列的特征值谱通常长这样：第1个特征值特别大，对应长期趋势；第2、3个特征值一组，代表年周期；第4、5个一组，代表半年周期；从第6或第7个开始特征值大幅缩小，进入缓慢衰减的"长尾"，这部分基本就是噪声。

另一个判断标准是累计方差贡献率。保留前d个主成分后，它们对轨迹矩阵总方差的贡献比例，用公式表示就是前d个奇异值平方之和除以全部奇异值平方之和。对GRACE典型序列，前5-6个主成分通常能解释85%-95%的方差，剩余的都是噪声和局地异常。你要是保留太多成分，比如d=20，那噪声也被当成信号重建回去了，插值效果会恶化；保留太少，又可能把半年周期甚至部分季节细节丢掉，插值结果会过于平滑。

我看过一些新手的做法是把d固定成5，然后所有格点、所有月份都套用，这是不推荐的。因为不同格点的时间序列特征差异很大：干旱区序列噪声小、趋势和年周期干净；而某些季节积水区、冰盖边缘区，信号复杂度高，d=5反而不够。稳妥做法是每根时间序列先做个快速SSA分解，看特征值谱，再决定保留个数。批量格点处理时，可以按贡献率阈值（比如解释85%方差）自动选d，这样既省事又比固定d更稳妥。

3.4 为什么SSA对缺失填充天然有效

把缺失点补上，本质上是在问：这些缺失月份的值，最可能是什么？局部插值只看缺失点附近的几个点，而SSA问的是全序列的标准答案。如果序列确实由趋势和少数周期成分组成，那么任何位置的值都理应满足这些成分在时间和振幅上的约束。缺失点的合理值，就是让整条序列在低维空间里最"自洽"的那个值。

从这个角度看，SSA插值对长缺口的优势就很清楚了。假设2015年7月到2016年3月连续缺失，线性插值只能连一条从数据起点到终点的直线或曲线，把真实季节循环压制掉。SSA迭代插值则会借助其他年份同期的信息——年周期被成功分解出来后，重建信号会自然地把这个窗口内的季节峰谷也恢复出来。当然，如果这个窗口内有特殊气象事件（比如极端干旱），SSA无法从历史信息中"发明"出这种异常，它会给出一个平滑的推定值。对这个缺点，我建议做完SSA插值后，再用独立的地面水文观测或再分析数据对插入值做合理性检查。

4. 基于MATLAB的完整实操流程

4.1 读取真实GRACE数据并整理时间轴

GRACE月度网格产品通常以NetCDF格式发布，常见变量名有lwe_thickness、weird等，不同中心产品命名不一样，读取前先用ncinfo查一下文件结构最稳妥。下面是读取代码的骨架，以1°×1°网格为例：

filename = 'GRACE_200204_201706_lwe_thickness.nc'; lon = ncread(filename, 'lon'); lat = ncread(filename, 'lat'); time = ncread(filename, 'time'); % 单位通常是 days since 2002-01-01 lwe = ncread(filename, 'lwe_thickness'); % 维度一般是 lon x lat x time % 时间轴转换 t0 = datetime('2002-01-01'); time_axis = t0 + days(time(1:size(lwe,3))); % 质量控制和填充值处理：把无效值全部转成NaN lwe(lwe <= -9999) = NaN; lwe(abs(lwe) > 10000) = NaN;

这段代码有一些细节值得说。第一，不要对time数组假设默认坐标，因为有些产品的时间基准不是2002-01-01，而是其他日期，建议先看ncreadatt(filename, 'time', 'units')确认。第二，原始网格里通常把陆地上的非反演区域设置为填充值，我的经验是直接把这些位置记为NaN，后续所有统计工具对NaN的处理会更可控。第三，对于squeeze过来的二维网格，务必确认数据维度顺序是lon × lat × time，不同产品可能有time × lon × lat，读出来就索性地用squeeze(lwe(i,j,:))这种按实际索引取数，避免混乱。

提取单个格点序列也很直接，比如我常研究华北平原某个格点，大致经纬度为(115°E, 35°N)，可以用以下方式定位最近格点：

idx_lon = find(abs(lon - 115) == min(abs(lon - 115))); idx_lat = find(abs(lat - 35) == min(abs(lat - 35))); ts = squeeze(lwe(idx_lon, idx_lat, :));

记得把ts也过一遍质量控制：如果原始变量和掩膜已经读了，这里主要将无效填充值转成NaN。

4.2 编写SSA分解函数

我在项目里用的SSA分解函数，核心就是前面讲到的嵌入、SVD和重构三步。这里的diag_averaging是对角平均操作，负责把轨迹矩阵还原为时间序列。代码可以直接复制使用。

function [SIG, recon, relvar] = ssa_decompose(x, L) % SSA分解：对输入时间序列x进行奇异谱分析 % 输入： % x - 列向量，长度N，不含NaN（调用前先做好填充） % L - 窗口长度 % 输出： % SIG - N x L 矩阵，列向量为每个主成分对应的时间序列分量 % recon- 所有成分叠加后的重建序列 % relvar-每个分量的方差贡献比例 x = x(:); N = length(x); K = N - L + 1; % 嵌入：构造轨迹矩阵 X = zeros(L, K); for i = 1:K X(:, i) = x(i:i+L-1); end % SVD分解 [U, S, V] = svd(X, 'econ'); sigma = diag(S); % 逐个成分重构 SIG = zeros(N, L); for i = 1:L Xi = sigma(i) * U(:, i) * V(:, i)'; SIG(:, i) = diag_averaging(Xi, N); end relvar = sigma.^2 / sum(sigma.^2); recon = sum(SIG, 2); end function s = diag_averaging(X, N) % 对角平均：将轨迹矩阵转换为等长时间序列 [L, K] = size(X); idx = (1:L)' + (1:K) - 1; % 每个元素对应的对角序号 s = accumarray(idx(:), X(:)) ./ accumarray(idx(:), ones(numel(X), 1)); s = s(1:N); end

这里用accumarray做对角平均是整个代码的精华点，速度比循环快得多，在批量处理几百个格点时能明显感受到差距。写到这里你可能会问，为什么轨迹矩阵是L×K而不是K×L？两种构造方式在数学上等价，只要SVD后重构时对角平均对应上就行。我习惯用L行、K列，这样奇异向量U的维度是L×L，和窗口长度对应，更容易画图解释每个主成分的形态。

4.3 迭代SSA缺失填充函数

核心的缺失填充函数如下。它的循环逻辑是：先对NaN点做PCHIP初始填充，然后不断调用ssa_decompose，用前d个主成分的重建结果更新缺失点，直到两次迭代差异足够小。

function [x_filled, info] = ssa_fill_missing(x, L, d, maxIter, tol) % 基于迭代SSA的缺失时间序列插值 % 输入： % x - 输入时间序列，包含NaN表示缺失 % L - SSA窗口长度 % d - 保留的主成分个数 % maxIter - 最大迭代次数 % tol - 收敛阈值 % 输出： % x_filled - 插值后的完整序列 % info - 结构体，包含迭代过程信息 x = x(:); N = length(x); miss = isnan(x); t = (1:N)'; % 初始填充：PCHIP插值 x_init = x; if any(miss) x_init(miss) = interp1(t(~miss), x(~miss), t(miss), 'pchip'); end x_old = x_init; history = zeros(maxIter, 1); for iter = 1:maxIter % SSA分解 SIG = ssa_decompose(x_old, L); recon_main = sum(SIG(:, 1:min(d, L)), 2); % 前d个主成分重建 % 更新缺失点 x_new = x_old; x_new(miss) = recon_main(miss); % 计算缺失点更新量 delta = norm(x_new(miss) - x_old(miss)) / (norm(x_old(miss)) + eps); history(iter) = delta; x_old = x_new; if delta < tol break; end end x_filled = x_old; info = struct('iter', iter, 'delta_history', history(1:iter), ... 'final_recon', recon_main); end

这段代码有几个关键点需要解释。

一方面，我迭代时更新的只是缺失点，非缺失点始终保留原始观测值。这是有意为之，因为SSA重建结果本身会有轻微偏差，如果拿重建值覆盖全部观测值，等于把去噪无条件的强加了，会损失真实观测中的有效信息。只更新缺失点，相当于把SSA当成一个"向导"，而非"替代者"。

另一方面，收敛判据只看缺失点值的相对变化，不看整条序列的变化。因为非缺失点本来就不变，整条序列变化反而不敏感。迭代次数一般不需要太多，我实测10-30次内就能收敛，tol设到1e-5已经足够。设太小不必要，徒增计算时间。

关于初始填充的选择，PCHIP比线性好，主要因为它更接近真实序列的峰谷形态。如果你愿意，也可以用别的方法初始化，比如气候态平均值填充（该月份的多年平均值），然后迭代SSA也能收敛。但我建议尽量用PCHIP，因为它在短期波动上更合理，初始值越好，迭代收敛越快。

4.4 效果评估：挖洞测试（交叉验证）

插值结果是用来做后续分析的，不验证就交出去，等于把一个未标定的传感器装进系统里。对GRACE这种有大量真实观测月份的序列，最实用的验证方法是挖洞测试：

1. 取出原始时间序列中所有有效月份。
2. 随机挖走其中10%的有效值，标记为NaN。
3. 用同一套SSA参数做插值。
4. 把插值结果和挖走前的真实值比较，计算RMSE、相关系数、偏差等指标。
5. 重复多次（比如50次），看统计稳定性。

代码写起来也不复杂：

rng(2024); ts_valid = ts(~isnan(ts)); valid_idx = find(~isnan(ts)); test_num = max(3, round(0.1 * length(valid_idx))); rmse_list = zeros(50, 1); corr_list = zeros(50, 1); for k = 1:50 test_idx = randsample(valid_idx, test_num); ts_test = ts; ts_test(test_idx) = NaN; ts_fill = ssa_fill_missing(ts_test, L, d, 50, 1e-5); truth = ts(test_idx); pred = ts_fill(test_idx); rmse_list(k) = sqrt(mean((pred - truth).^2)); corr_list(k) = corr(pred, truth); end fprintf('RMSE: %.2f ± %.2f cm\n', mean(rmse_list), std(rmse_list)); fprintf('Corr: %.3f ± %.3f\n', mean(corr_list), std(corr_list));

我做完挖洞测试后还会做一件事：画"预估误差随缺失长度变化"的关系图。做法是故意挖掉长度为1、3、6、12个月的连续段，统计插值RMSE，看误差增长曲线。这个图能帮你判断哪些月份的缺失补出来可信，哪些只能当参考。如果连续缺12个月补出来的RMSE已经超过信号振幅的一半，那这段插值结果在定量分析里就必须降权或者剔除。

4.5 批量处理多个格点的效率优化

GRACE网格全球有6万多个格点，如果每个格点都调用ssa_fill_missing，在MATLAB里循环跑会非常慢。我实际处理时做了两件事。

第一，只对有效格点计算。陆地冰盖、海洋、沙漠里很多格点要么长期为NaN，要么是纯噪声，先做一个简单的有效月份数量统计，比如有效月份占比小于60%的格点直接跳过，只插值有效月份占比高的格点，计算量能省一半以上。

第二，把ssa_fill_missing的SVD部分向量化，或者用parfor并行。最简单的改动就是把外层循环改成parfor，但要注意parfor要求循环体内的代码不能依赖前次迭代结果、不能修改共享变量，我的函数是纯计算，完全满足。下面是批量处理的示意代码：

% 假设 lwe 是 lon x lat x time 的三维数组 [LonN, LatN, T] = size(lwe); L = 36; d = 6; lwe_filled = lwe; valid_ratio = squeeze(sum(~isnan(lwe), 3) / T); parfor i = 1:LonN for j = 1:LatN if valid_ratio(i, j) < 0.6 continue; end ts = squeeze(lwe(i, j, :)); if sum(~isnan(ts)) < L continue; % 有效点数太少，SSA分解无意义 end ts_fill = ssa_fill_missing(ts, L, d, 50, 1e-5); lwe_filled(i, j, :) = ts_fill; end end

这里有个注意事项：parfor循环里每次调用ssa_fill_missing都会做SVD，计算量仍然比较大。如果机器内存足够，可以先提取出所有需要插值的格点序列到一个大的二维数组，再用parfor并行处理，这样能减少一部分重复读取三维数组的开销。我试过同一批数据，从串行改成parfor，在8核机器上大概能快三到四倍。

5. 模拟和真实序列的效果对比与细节提醒

5.1 用模拟数据测试不同插值方法的误差

我先构造一条接近GRACE行为的模拟时间序列，用来系统比较各方法的精度。模拟序列长度设为180个月，包含一个线性趋势、一个年周期、一个半年周期，再叠加高斯白噪声：

t = (1:180)'; trend = -0.03 * t; % 长期趋势 season = 3 * sin(2*pi*t/12) + 1.5 * cos(2*pi*t/12); semi = 0.8 * sin(4*pi*t/12 + 0.5); noise = 0.7 * randn(size(t)); ts_true = trend + season + semi + noise;

接下来我按三种情景挖洞：

• 情景A：随机挖15个孤立月份；
• 情景B：连续挖掉6个月（比如第100-105月）；
• 情景C：挖掉两个长段（第40-48月、第130-140月）。

然后分别用PCHIP插值和迭代SSA插值（L=36，d=6）重建，计算与真值的RMSE。结果如下表：

从这个模拟结果能明显看出，随机零散缺失时，PCHIP其实够用，SSA优势不显著；但一旦出现连续缺口，PCHIP会把缺口内季节循环削平，RMSE显著恶化，而SSA由于利用了全序列的周期信息，误差增长慢得多。这个结果也解释了为什么我不建议所有情况都无脑上SSA——计算量高、参数要调，但对零散缺失并没有比PCHIP好多少。合理的策略是：先统计缺失结构，再决定方法。

5.2 真实GRACE格点序列的插值效果

拿真实的GRACE格点序列来看，比如某华北地区格点，原始序列从2002年4月到2017年6月，中间缺了大约24个月，其中有两个接近连续的缺口，分别出现在2012年（约3个月）和2016年末到2017年初（约8个月）。我用PCHIP和SSA分别插值后发现，SSA补出的2016年末-2017年初那一段呈现出一个较平滑的下降-反弹过程，与GRACE-FO后续观测的趋势衔接自然；而PCHIP补出来的是一段近似直线，明显削弱了季节性。

更关键的差异在后验验证里。我把2014-2015年这段已有观测的月份挖掉6个点做测试，SSA补出来的RMSE是1.2cm，PCHIP是1.8cm。由于该时段本身包含一次强降水恢复信号，PCHIP把恢复过程拉平了，SSA则较好地保留了恢复幅度。这正是真实序列中值得关注的差异——那段时间恰好对研究干旱恢复有重要意义，插值方法选错了，结论会有实质差异。

当然，真实数据插值永远不等于真实观测。对SSA补出的重要信号（比如某年汛期水位是否恢复），我都会去查阅该区域的地面水文站数据或再分析降水数据，做一次"旁证"再写进论文。

5.3 SSA插值容易翻车的四个细节

第一个是边缘效应。SSA重构在序列首尾两端误差较大，因为末尾段的窗口片段数量少，SVD分解对这段的约束不足。如果缺失发生在序列末端，插值结果不确定性很大。我的处理办法是把插值结果的两端（各L个月）单独标记，在做趋势分析时不对这两年过度解读。如果非要用末端插值结果，建议用较短窗口（L=24）重新做一个敏感性测试。

第二个是主成分数d的敏感性。同样一组数据，d从4调到8，插值结果在长缺口的差异可以达到1-2cm。不要觉得d越大越精细，噪声多的高频成分会被当成有效信号。我一般做法是：先看特征值谱，d取在拐点前一位；再做一个d从4到10的敏感性扫描，看趋势和季节振幅变化大不大。如果某格点对不同d的结果差异特别大，说明这个格点信号本身就不稳定，插值结果要谨慎处理。

第三个是窗口长度L与序列长度的关系。如果序列太短（有效月份只有80多个月），L取36会造成K很小（60左右），每个分量的估计方差增大。这时把L降到24会更稳。一个粗略的下限是：L最好不超过N/3，同时L要大于目标周期的两倍。两个条件冲突时，优先保证目标周期，再放宽到L≈N/3。

第四个是数据里有强异常值没有先行剔除。GRACE网格产品在某些月份会受海洋混叠误差、冰盖信号泄漏影响，让孤立月份出现特别大或特别小的异常值。这些值如果不去掉就参与SSA分解，会污染整个主成分结构，插值结果也会被连带扭曲。我的建议是：插值前先做一遍绝对离群值检测（比如超过五年滑动标准差4倍以上的点标记为缺失），再走SSA流程。当然，对于真实地球物理信号，异常值不一定就是噪声，需要结合具体事件判断，不能机械剔除。

这套SSA插值流程我从最早接触GRACE数据处理时就开始用，后来做过水文干旱研究、地下水储量趋势分析，凡是涉及月度序列补全的，都用这套方法统一处理。它的优点在于把"插值"和"去噪"合在了一步，省事，也避免了两阶段处理带来的信号交叉干扰。实际动手算的时候，我建议你先拿一个格点跑通流程，把PCHIP和SSA结果都画出来对比，再决定要不要对整片区域批量处理。代码、参数都是可以复现的，但数据里那些说不清道不明的异常，终究需要你亲眼盯着图看一遍。

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