简介本资源是一套面向地球物理与水文研究者的GRACE Mascon数据缺失月份插值工具包聚焦于利用奇异谱分析SSA算法实现时间序列重建适用于科研人员及研究生开展区域陆地水储量变化分析。压缩包共14个文件11个MATLAB脚本、1个说明文档、1个NetCDF测试数据、1张结果示意图总大小71.53MB其中核心脚本涵盖时间标准化、年份转换、SSA迭代填充、重力场可视化等全流程模块支持直接运行并复现插值结果。已有506人学习下载配套简易文档与专栏延伸内容便于理解SSA参数敏感性如窗口长度、滞后阶数对重力变化敏感区插值精度的影响。用户可快速掌握GRACE数据预处理—缺失填补—结果验证的技术链路并对比神经网络等多元插值方法的适用边界。 做GRACE水文信号处理的这几年我几乎每次都要面对同一个问题数据缺月份。GRACE卫星的月度重力场解算结果总有几个月份是空白的早期任务时期尤其严重。你要算陆地水储量变化趋势缺了月份的序列很难直接扔给趋势分析你要做季节周期的提取断点更是让人头疼。最开始我的做法很粗暴拿线性插值把坑补上就完事后来发现这样补出来的信号在季节性变化比较剧烈的区域完全没法看。后来逐步接触了奇异谱分析SSA用SSA做缺失数据的插值和信号重建效果好了不止一个档次。这篇文章就把我整套处理流程写出来从GRACE数据读取、缺失月份插值到基于MATLAB的SSA实现全部捋一遍希望能帮到正在踩坑的同行。做这个项目的起因很具体我手上有一份2002到2017年的区域GRACE数据用的是RL06 mascon产品中间缺了大概7个月。目标是提取区域陆地水储量的长期趋势和季节振幅同时评估那几年的干旱事件。这个需求看起来简单但实际操作中牵涉到的细节远比想象中多数据格式怎么读、网格怎么裁剪、缺失月份怎么补、补完之后怎么验证、SSA的窗口参数到底怎么选。这篇文章是完整的复盘记录适合刚接触GRACE数据、或者已经拿到数据但卡在预处理环节的研究生和从业者。1. 项目概述GRACE数据为什么需要插值和SSA1.1 GRACE数据的典型缺失场景GRACEGravity Recovery and Climate Experiment重力恢复与气候实验卫星从2002年运行到2017年与后续的GRACE-FO任务共同构成了近20年连续的地球重力场观测序列。原理上就是两颗卫星在同轨道上前后跟飞通过测量星间距离变化反演地球重力场的时空变化再换算成陆地水储量异常。这个数据在干旱监测、地下水可持续性评估、冰川消融速率估算这些方向上都是核心输入。但数据并不是完整的月度序列。即使在正常运行期也会因为卫星平台的能源限制、姿控系统需要调整、星上仪器校准、轨道机动等原因跳过某些月份。比如2011年上半年就有连续几个月的缺口2016年底到2017年初也有明显断档这还不包括早期2002到2004年那段频繁切换观测模式的时间。GRACE-FO在2018年发射之后同样面临电池管理问题数据缺口依然存在。具体到我做的这份数据原始文件是按月度netCDF格式存储的全球0.25度网格时间覆盖2002年4月到2017年6月共183个月但实际只有176个有效月份缺失了7个月的观测值。这个缺失比例在GRACE任务里算是正常的不算严重但7个连续或者分散的坑足以让后续的时间序列分析出现偏差。1.2 为什么选择SSA而不是直接线性插值先说个我踩过的坑。第一次处理缺失月份我图省事用了MATLAB自带的fillmissing函数指定spline方法。补完之后趋势看起来没问题但把季节循环信号提取出来一对比发现季风区的振幅被明显压低了相位也偏移了一两个月。原因很简单三次样条插值本质是基于局部曲线拟合它对缺失段前后的数据点做光滑连接完全没有利用GRACE时间序列本身的强周期特性。GRACE陆地水储量信号有一个非常明确的季节循环通常一年一个峰谷还有一个缓慢的长期趋势外加一些年际波动。这种信号结构样条插值是不懂的。SSA的思路完全不同。它先把一条时间序列通过嵌入操作变成一个轨迹矩阵然后对这个矩阵做奇异值分解把原序列分解成若干个可以解释的成分长期趋势、季节周期、年际波动、噪声。在缺失数据场景下SSA可以利用整条序列的结构信息来推断缺失位置的合理取值。换句话说线性插值只看缺失点附近两个点SSA看的是全序列的统计结构补出来的值既符合局部走向也符合全局的季节规律。对比一下我做过的几种插值方法的效果线性插值在短缺口1-2个月场景下问题不大但3个月以上就会明显失真三次样条平滑但过冲在信号变化剧烈的区域会补出离谱的极值DCT-PLS离散余弦变换-惩罚最小二乘效果不错但参数敏感性较高基于PCA重建本质上是把多变量信息压缩后重建如果只是单点时间序列就用不上。SSA是单变量时间序列缺失插值里的综合最优解不需要外部先验信息原理透明、可解释性强而且能顺便完成信号分解一石二鸟。1.3 MATLAB在处理这类任务中的角色选择MATLAB首先是因为GRACE数据处理社区有大量现成的MATLAB工具比如常用的read_gracenetcdf、grace_mascon这类脚本网上能找到很多版本省去了从头解析netCDF格式的时间。其次是MATLAB的矩阵操作和内置函数对这个工作流太友好了读netCDF一行代码、区域掩膜是逻辑索引、矩阵转置和切片随手就来画图出图也方便。唯一要提醒的是处理全球网格逐点SSA时MATLAB的循环效率是短板后面我会讲到怎么用矢量化解决。2. GRACE数据获取与预处理要点2.1 数据源选择与格式解析GRACE数据产品按处理机构分成几个主要来源CSR美国德州大学奥斯汀分校空间研究中心、JPL美国喷气推进实验室、GFZ德国地学研究中心以及后来加入的GSFC美国宇航局戈达德太空飞行中心。产品形式又分为球谐系数SH和网格化 mascon 两类。SH系数适合做全球球谐分析、需要自己滤波的场合mascon则是直接解算到空间网格上包含了对信号泄漏的约束处理使用上更方便尤其适合区域水文研究。我使用的是JPL RL06 mascon版本网格分辨率0.25度变量名为lwe_thickness单位是等效水高厘米含义是相对于2004到2009年基准期的异常值。netCDF文件里的维度包括time、latitude、longitude。读取的时候用MATLAB的ncread函数就能搞定要注意的就是坐标系是0到359.75还是-180到179.75不同版本不一样建议统一转换到-180到180。这个细节很基础但容易出错我见过不少人在区域裁剪时把经纬度范围弄反最后画出来的图位置完全不对。如果你拿到的是SH系数产品还需要进行高斯滤波和泄漏校正流程会长很多。所以对没有特殊需求的场景我建议直接用mascon版本预处理工作量小一个量级。2.2 区域裁剪与陆地掩膜处理做区域研究时不需要全球网格。我的研究区是华北平原某区域经纬度范围大约在东经112°到119°、北纬34°到40°。裁剪的时候用逻辑索引直接提取子网格就行。但一个容易被忽略的问题是陆地掩膜GRACE观测的是地球总重力场变化但某个像元如果位于海洋或湖泊上信号含义就完全不同。如果做的是陆地水储量分析一定要先加载一份陆地掩膜数据把非陆地像元排除掉。我用的掩膜来自mascon产品自带的mask变量网格为0或1。处理逻辑很直观mask(lon_range, lat_range) 1才保留对应像元。另外要注意的是GRACE数据在赤道附近和沿海地区的信噪比偏低即使做了掩膜沿海网格也可能混入海洋信号。这种情况只能在后续分析中通过空间平滑或者灵敏度加权来缓解没有完美方案。裁剪完成后我一般会做一次简单的空间平均或多像元平均。区域平均的好处是抵消了部分随机噪声让时间序列更干净。具体做法是求出研究区内所有陆地像元在每个月份的平均等效水高得到一个长度为月份数的一维时间序列。后续的插值和SSA都在这条序列上进行。如果你更关心空间分布比如要做逐像元的趋势图那就要在每个像元上分别跑一遍相同流程计算量会大很多这个我在第5章的实操部分再展开。2.3 趋势项、季节项与年际信号的分离思路拿到完整的时间序列之后正式做SSA之前有一个概念要提前建立GRACE陆地水储量信号的时间结构大致包括四个部分——长期趋势比如长期的超采导致地下水持续下降、季节项降水和蒸散发的年内循环、年际异常ENSO这类气候现象造成的波动、残差噪声观测误差和短时气象扰动。SSA在处理这些成分时效果很自然经过合理选择窗口长度和重构阶数趋势成分会被提取为一阶或者二阶的缓慢变化分量季节项会被分解为一组周期相近的正交分量噪声则集中落在贡献率很低的高阶奇异值上。这个信号分离思路是整个插值方案的理论基石。因为当数据缺失时SSA之所以能“猜”出合理值靠的就是先用完整月份的数据把趋势和周期结构学出来再去预测缺失段。你在脑子里把这个逻辑理顺了后面调参数时就知道自己在干什么而不是拿着工具箱盲目试。3. 缺失月份插值方法对比与SSA原理3.1 常见插值方法的优缺点对比在正式讲SSA之前我把实践过的几种方法放在一起做个对比。这个表格可以当成以后选型的快速参考。方法原理优点缺点适用场景线性插值相邻点直线连接简单、稳定、无参数忽略季节结构长缺口失真1-2个月临时填充三次样条分段三次多项式光滑连接平滑、连续可导可能过冲震荡严重信号平缓区域的临时补值DCT-PLS离散余弦变换惩罚最小二乘能处理长缺口自动化程度高平滑参数需反复尝试长连续缺失季节性分解插值先估计季节项再补残差利用季节性先验季节项估计本身受缺失影响季节规律极为稳定的区域SSA迭代插值轨迹矩阵SVD重构迭代结构自适应无需先验窗口参数需要经验单变量时间序列综合场景SSA胜在不需要你提前假定趋势是线性的还是季节项是标准的正弦波。它直接从数据本身提取可重复的结构这一步在GRACE数据上比任何需要先验假设的方法都更稳妥。3.2 SSA基本数学原理奇异谱分析Singular Spectrum Analysis的核心步骤可以拆成四步嵌入、分解、分组、重构。嵌入操作是把原始时间序列转换成一个轨迹矩阵。假设序列长度为N选择一个窗口长度L1 L N令K N - L 1得到一个L行K列的矩阵X[ X \begin{bmatrix} x_1 x_2 \cdots x_K \ x_2 x_3 \cdots x_{K1} \ \vdots \vdots \ddots \vdots \ x_L x_{L1} \cdots x_N \end{bmatrix} ]轨迹矩阵的每一列都是原序列的一段长度为L的子序列。这个矩阵的一个重要性质是汉克尔结构——副对角线上的元素相等。这是后面重构时用来还原时间序列的关键。分解阶段对这个轨迹矩阵做奇异值分解[ X U \Sigma V^T ]其中U是L×L的正交矩阵V是K×K的正交矩阵Σ是对角矩阵对角线元素是奇异值。奇异值按大小降序排列每个奇异值对应的(U_i)和(V_i)构成一对EOF和主成分。这一步在MATLAB里就是svd函数一行的事但理解背后的几何意义很有帮助SVD把轨迹矩阵的行空间和列空间分解成互相正交的方向这些方向分别对应着原序列在不同时间尺度上的信号模式。分组阶段根据奇异值的大小和对应的频率特征把分解出的成分归类。比如前两个奇异值对应的EOF呈现缓慢变化的那一类就是趋势项接下来频域上呈现年周期的就是季节项后面奇异值很小、没有明显周期的基本就是噪声。重构阶段是把选定的成分组合起来通过对角平均diagonal averaging把矩阵还原成时间序列。对角平均的意义在于由于轨迹矩阵有汉克尔结构还原回去的时候对副对角线上的元素取平均就得到了原序列对应位置的重构值。对每个选定的成分组做一次对角平均就得到了对应的重构成分RC。把所有R C加起来就得到完整的重构序列。3.3 Gap-Filling SSA插值的迭代框架有了上面的基础处理缺失数据的框架就顺理成章了。这个方法在文献里常被称为Iterative SSA gap-filling思路非常直观第一步用简单的线性插值补上缺失位置的初值得到一个完整的临时序列。第二步对整个临时序列进行SSA分解选择前k个主成分重构出信号。第三步把重构序列中缺失位置的值提取出来替换掉临时序列里对应的旧值。第四步重复第二步和第三步直到缺失位置的值变化足够小也就是收敛了。为什么这个迭代过程会收敛因为每做一次重构信号中的结构信息就在被更精确地刻画而被噪声污染的高阶分量被剔除缺失位置的值会逐步靠近信号的真实内在结构。实际操作中迭代次数通常控制在50次以内就能收敛。要注意的是这个方法假设缺失比例不能太高经验上限大概是总长度的30%左右。如果你的数据缺口超过三分之一SSA能用来学结构的完整观测太少补出来的结果可靠性就会明显下降。4. 基于MATLAB的核心实现与参数选择4.1 SSA分解函数完整代码先给一个通用的SSA分解函数。这个函数输入时间序列x和窗口长度L输出全部重构成分RC、奇异值eigenvalues和对应的EOF向量。function [RC, eigenvalues, EOFs] ssa_decompose(x, L) % SSA分解函数 % 输入: % x - 长度为N的一维列向量或行向量 % L - 嵌入窗口长度建议 N/4 ~ N/3 % 输出: % RC - N×L矩阵每一列是一个重构成分 % eigenvalues - 奇异值降序排列 % EOFs - L×L矩阵每一列是对应时间EOF模式 x x(:); % 转为行向量 N length(x); K N - L 1; % 构建轨迹矩阵 X zeros(L, K); for i 1:K X(:, i) x(i:iL-1); end % 奇异值分解 [U, S, V] svd(X, econ); eigenvalues diag(S); EOFs U; % 计算各重构成分 RC zeros(N, L); for k 1:L % 第k个成分对应的轨迹矩阵 Xk U(:, k) * S(k, k) * V(:, k); % 对角平均还原为时间序列 RC(:, k) diag_averaging(Xk, N, L, K); end end function y diag_averaging(Xk, N, L, K) % 对角平均还原时间序列 y zeros(1, N); count zeros(1, N); for i 1:L for j 1:K y(ij-1) y(ij-1) Xk(i, j); count(ij-1) count(ij-1) 1; end end y y ./ count; end这段代码里两点值得注意。一是对svd使用了econ选项它只计算前K列如果K L变量较少时能省不少内存。二是对角平均要除以计数count而不是直接平均因为原序列首尾位置的元素在轨迹矩阵中出现的次数不同中间位置的元素出现的次数多必须加权平均才能正确还原。4.2 缺失数据迭代插值的实现下面是迭代gap-filling的实现。核心逻辑是用线性插值做初始化然后循环调用SSA分解、用前k个成分重构、更新缺失点。function [x_filled, x_recon, rms_change] ssa_gapfilling(x, L, k, maxiter, tol) % 基于SSA的缺失数据迭代插值 % 输入: % x - 原始序列缺失位置用NaN标记 % L - 嵌入窗口长度 % k - 重构阶数即保留前k个成分 % maxiter - 最大迭代次数 % tol - 收敛阈值 % 输出: % x_filled - 填充完成的时间序列 % x_recon - 最终SSA重构信号 % rms_change - 每次迭代缺失位置的变化量 x x(:); N length(x); miss_idx isnan(x); t 1:N; % 初始化线性插值 if all(miss_idx) % 全缺失时直接报错 error(输入序列全部为NaN无法插值); end x_filled x; x_filled(miss_idx) interp1(t(~miss_idx), x(~miss_idx), t(miss_idx), linear); % 迭代 rms_change zeros(maxiter, 1); for iter 1:maxiter % 对当前完整序列做SSA分解 [RC, ~, ~] ssa_decompose(x_filled, L); % 用前k个成分重构 x_recon sum(RC(:, 1:k), 2); x_recon x_recon(:); % 更新缺失位置 x_new x_filled; x_new(miss_idx) x_recon(miss_idx); % 收敛判断 change x_new(miss_idx) - x_filled(miss_idx); rms_change(iter) sqrt(mean(change.^2)); if rms_change(iter) tol x_filled x_new; break; end x_filled x_new; end end实际调用的时候我建议先做一次参数敏感性测试把maxiter设置成100tol设置成1e-4。迭代过程通常在30到60次收敛。如果跑到100次还是没收敛多半是k选得太大导致重构序列里混入了噪声或者是缺失比例太高需要降低k值或者考虑分段处理。4.3 嵌入窗口L与重构阶数k的选择经验这是SSA应用中最容易卡住的地方我把自己的经验整理成几条可操作的规则。嵌入窗口L的选择直接决定了SSA能分辨的频率范围。L太小周期信号无法在轨迹矩阵里完整展开L太大轨迹矩阵尺寸膨胀低阶奇异值拖尾严重计算量也大。经验法则是取序列长度的1/4到1/3。对GRACE这种月分辨率数据L取36个月左右比较合适因为三年窗口能覆盖完整的季节循环和大部分年际波动。如果序列足够长超过200个月还可以考虑L48或60能更好地解析两年以上尺度的信号但代价是边界效应更明显。重构阶数k的选择是SSA应用里最重要的参数。理论上k应该等于信号成分的个数但实际观测中噪声和信号并不会完全分离很难严格确定。我常用的方法有两种。第一种是奇异值谱拐点法画出奇异值大小随序号变化的曲线找明显的拐点拐点之后的奇异值趋于平缓这部分就是噪声。通常GRACE序列的奇异值谱中第一个拐点出现在第1到第3个奇异值处对应趋势第二个拐点在季节频率几个分量之后可能是第5到第10个。第二种是试错法分别用几组k值比如k5、10、15、20做插值然后计算残差平方和和重构序列与原始观测在非缺失位置的相关系数。k太小会欠拟合去掉有效信号k太大会把噪声当信号出现振荡毛刺。我自己的经验是处理GRACE月度水储量数据时k取10到12通常表现不错。这个区间能把趋势项、年度周期项和主要年际项全部纳入同时过滤掉大部分随机噪声。如果研究区水文信号较强比如季风区k可以适当加大如果信号弱比如干旱半干旱区k要小一些防止噪声主导。5. 完整实操一个区域的GRACE数据处理全流程5.1 读取netCDF并生成区域平均时间序列现在把前面的模块串成一个完整流程。第一步读取GRACE mascon数据。% 读取GRACE mascon netCDF文件 filename GRACE_JPL_MASCON_RL06_v02.nc; lwe ncread(filename, lwe_thickness); % 月平均等效水高 lon ncread(filename, lon); lat ncread(filename, lat); time ncread(filename, time); % 通常为月序数或日期 mask ncread(filename, mask); % 海洋为0陆地为1 % 坐标转换有些产品lon范围是0~360 lon(lon 180) lon(lon 180) - 360;第二步区域裁剪并取陆地掩膜。我研究区设为华北平原。lon_idx lon 112 lon 119; lat_idx lat 34 lat 40; lwe_region lwe(lon_idx, lat_idx, :); mask_region mask(lon_idx, lat_idx); % 只保留陆地像元 lwe_region_vals lwe_region(:,:,:); land_mask_3d repmat(mask_region, 1, 1, size(lwe_region, 3)); lwe_region_vals(~land_mask_3d) NaN;第三步求区域平均。由于存在掩膜和可能的缺测不能用简单的mean要逐月对非NaN像元求均值。nt size(lwe_region, 3); region_ts zeros(nt, 1); for t 1:nt tmp lwe_region_vals(:,:,t); region_ts(t) mean(tmp(~isnan(tmp)), omitnan); end % 把缺测月份设置为NaN后续交给插值函数 region_ts(region_ts 0) NaN;注意最后这一步region_ts 0在GRACE mascon里通常是无效值或海洋遮挡要主动转成NaN避免把真实0值当成有效数据。5.2 缺失月份识别与SSA插值执行生成区域平均序列之后先看一下缺失月份分布再调用插值函数。% 识别缺失月份 miss_idx find(isnan(region_ts)); fprintf(缺失月份数量: %d\n, length(miss_idx)); disp(miss_idx); % 设置SSA参数 L 36; % 嵌入窗口3年 k 12; % 保留前12个成分 maxiter 100; tol 1e-4; % 执行迭代插值 [x_filled, x_recon, rms_change] ssa_gapfilling(region_ts, L, k, maxiter, tol);迭代过程中rms_change会逐渐下降并趋于平稳。我处理的结果是迭代24步后rms_change从最初的0.9厘米降到0.01厘米以下28步后低于阈值。填充出来的缺失月份结果和前后月份衔接非常自然春季回升、夏季达到峰值、秋冬回落的季节特征完整保留。这一点用线性插值很难做到尤其在跨度达到两个月的缺失段上。5.3 交叉验证用什么指标说明插值效果插值不是补完就结束了必须验证。我自己最常用的是留一法交叉验证把原本完整的月份人为设为NaN用SSA插值再和真实值比较。实际操作中我在区域时间序列里随机选了12个完整月份做掩膜测试重复了10次。结果显示SSA插值的平均绝对误差在0.5到1.2厘米等效水高之间均方根误差约0.8厘米。相比之下三次样条在同一测试中的均方根误差为1.5厘米。对 GRACE 数据来说这个精度已经相当可接受。交叉验证的代码思路很简单% 在完整序列上随机挖掉12个月测试插值误差 rng(42); n_test 12; test_idx randperm(sum(~isnan(region_ts)), n_test); x_test region_ts; x_test(test_idx) NaN; [x_test_filled, ~, ~] ssa_gapfilling(x_test, L, k, maxiter, tol); err x_test_filled(test_idx) - region_ts(test_idx); rmse sqrt(mean(err.^2)); mae mean(abs(err)); fprintf(测试集RMSE: %.2f cm, MAE: %.2f cm\n, rmse, mae);如果RMSE明显高于该区域的GRACE数据本身噪声水平通常2~3厘米就说明插值引入了额外误差。这时候优先检查k是否过大或者L是否过短而不是怀疑方法本身。5.4 结果分解提取趋势和季节成分插值完成之后可以再用SSA对完整序列做一次最终分解获得趋势项和季节项。% 对填充后的序列重新做SSA分解 [RC, eigenvalues, EOFs] ssa_decompose(x_filled, L); % 趋势项通常取第一个和第二个RC相加 trend RC(:,1) RC(:,2); % 季节项取周期为12个月附近的成分 % 先看前几个EOF的傅里叶谱识别季节成分序号 seasonal_comps [3 4 5 6]; % 根据实际情况调整 seasonal sum(RC(:, seasonal_comps), 2); % 残差项 residual x_filled(:) - trend - seasonal;注意这里RC的列顺序是奇异值降序的但成分顺序并不严格对应时间尺度。比如第3和第4列可能组合成年度周期第5和第6列也可能是年周期的具体调制形式。判别方法是分别画出每个RC的波形和傅里叶频谱看它集中分布在哪个频带。这一步多花一点时间后面做物理解释时就不会出错。6. 常见问题与排查技巧实录6.1 SSA重构序列两端振荡怎么办这个问题几乎每个用SSA的人都会遇到。因为轨迹矩阵对角平均之后序列头部和尾部能用于平均的矩阵元素个数比中间少重构误差自然更大。振荡幅度在信号变化剧烈时会变得非常明显。两个可行的处理策略。第一在SSA分解之前对原始序列做端点扩展比如用镜像法延长序列两端各L个月分析完再把延展部分切掉。这个方法能有效削弱边缘效应但会引入额外的计算量。第二如果只是插值场景缺失位置在序列中间这个问题影响不大如果缺的就是两端就要特别小心最好用交叉验证确认插值结果是否合理。我自己的经验是对于端点缺失不要用SSA得到的重构值直接作为最终结果而是结合区域内的空间相邻网格信号来相互校正。6.2 插值结果出现负值或异常值怎么处理GRACE等效水高是相对于基准期的异常值本身可能为负所以看到负值不要慌。但如果插值结果出现明显超出物理范围的值比如一个干旱区的网格序列突然补出一个20厘米的尖峰基本就是SSA参数选得不对。排查步骤先检查k是否过大导致高频噪声成分混入重构信号再检查L是否过小导致周期结构没有被正确捕获第三看迭代是否收敛如果rms_change在后期来回震荡说明陷入了不稳定循环。最后一步如果还是异常把该网格的原始数据和相邻网格对比看是不是GRACE数据本身有这个尖峰可能是区域真实水文事件比如极端降水。6.3 海量网格点逐点计算太慢的加速方案如果你不是做区域平均而是对全球0.25度网格每个陆地点都跑一遍SSA插值那计算量非常大。全球陆地像元大约几十万个每个像元跑一次SSA迭代在普通PC上用循环可能要跑几天。几个实用的加速思路。第一不要用svd处理完整轨迹矩阵改用svds只求前k个奇异值对应的特征子空间能省一大截时间。第二利用parfor做并行循环MATLAB的并行计算工具箱在这里非常有用我实测8核并行能加速4到5倍。第三如果某个像元附近没有缺失月份直接跳过插值只对存在缺失的像元调用SSA插值函数很多区域其实只有少量网格点有缺失。第四ssa_decompose中对每个成分循环做对角平均时可以整体写成矩阵运算替代双层循环但代码可读性会差一些视你的具体需求取舍。6.4 SSA分解后的趋势项不单调正常吗很多人在做趋势提取时有一个预期趋势项应该是一条单调递增或递减的直线。但SSA提取的“趋势”不是数学意义上的直线趋势而是序列中的低频缓变分量。它可能是多段式的比如早期缓慢下降、中期平稳、后期加速下降。这实际上是水文信号的真实反应恰恰是SSA优于简单线性回归的地方它不强迫趋势是线性的。对于这种低频分量我习惯的做法是先提取RC1和RC2看看形态如果它们确实反映了长期的缓慢变化就合并为趋势项再进行一次线性回归计算趋势速率。这样既得到了非线性趋势形态又可以得到一个便于比较的线性速率指标用来和其他研究发表的速率对比。根据我个人实操的体会GRACE数据处理里最值得投入时间研究的环节不是GRACE数据获取也不是基础的矩阵运算而是缺失月份的插值策略。这个环节决定了后续一切时间序列分析的可信度。SSA迭代插值并不是一个特别复杂的算法但在GRACE这种具有强季节节律的数据上效果显著优于线性方法和样条方法而且顺便完成了信号分解一步解决两个需求。把这个流程跑通之后不管是做区域水储量变化评估还是为水文模型提供校准数据都会顺手很多。另外再分享一个小技巧不管用什么插值方案都要把参数记录在代码注释里并且每次换数据区域先做一轮交叉验证不要一路沿用旧参数。GRACE不同区域、不同时间段的信号结构差异很大参数需要跟着数据走。希望这篇文章能帮你少走一些弯路。本文还有配套的精品资源点击获取