的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插值的思路和局部插值完全不同。它先假设真实信号可以表示为少数几个主要成分的叠加趋势逐年周期半年周期少量噪声然后通过矩阵分解把时间序列拆成这些成分再用主要成分去估计缺失位置的值。具体到缺失填充比较成熟的做法是迭代式填充流程很简单先用线性或PCHIP给缺失点补一个初始值得到完整序列。对完整序列做SSA分解得到前d个主要成分的重构序列。用重构序列去更新缺失点的旧值。重复进行SSA分解和更新直到缺失点数值变化小于阈值。这个流程执行下来缺失值会逐步收敛到一个与整体信号结构自洽的状态。因为SSA重建过程中非缺失点的信息也参与了每个主成分的估计相当于用全序列的数据来约束缺失段比只用缺口边界局部的信息去猜要稳健得多。我后面会给完整可跑的MATLAB代码。2.3 进阶路线时空联合插值DINEOF等如果研究区域不是一个格点而是整个网格场还可以用DINEOFData Interpolating Empirical Orthogonal Functions这类方法。它的原理是利用时空场的经验正交函数结构来填充数据本质上可以看成SSA在高维空间里的扩展。DINEOF在海洋遥感数据SST、叶绿素里用得很多插值同时能保持空间连续性对大面积云覆盖导致的缺测特别好用。但对GRACE月度网格来说有一个现实问题如果是某个月份整个全球数据都缺失DINEOF也一样无从下手因为它需要利用已有的时间协方差结构。而且DINEOF的调参过程判断截断模态数、迭代收敛标准比SSA更复杂对刚接触数据处理的同学没那么友好。我通常的建议是先用SSA把单点时间序列补好如果后续要做空间场分析再考虑DINEOF这种更重的工具。2.4 到底怎么选一张选型表我把各种情况下的推荐方案总结如下。场景推荐方案理由缺1-2个孤立月份PCHIP插值简单、稳健、不会明显改变谱结构连续缺3-6个月迭代SSA插值L24或36能利用全序列趋势和周期信息降低振幅低估连续缺半年以上迭代SSA插值NCEP/ERA5辅助数据交叉验证长缺口本身不确定性大必须做敏感性分析整场网格有空间孔洞DINEOF或MSSA多变量SSA利用时空协方差保持空间相关性序列末端缺失SSA插值要格外小心重构在两端有较大边缘误差需结合外推或降低权重表格之外还有一点需要强调无论选哪种插值最终都要做交叉验证把已知的完整月份故意挖去用剩余数据重建再与真值比较。不做验证的插值参数随便一调都会看起来很漂亮但真实误差你心里没底。3. 奇异谱分析(SSA)的原理与参数选法3.1 SSA的一次完整计算过程SSA的核心思想是找一个合适的窗口把一维时间序列转成一个轨迹矩阵然后对这个矩阵做奇异值分解再用特征值较大的那几个分量重建原序列。我用一个日常类比来帮助理解。假设你有一盘混合音频录音里面有人的说话声、背景音乐和电流底噪。SSA做的事情就是把这盘录音切成许多等长的片段窗口然后分析这些片段中重复出现的模式——说话声和音乐是有规律的底噪是随机的。通过矩阵分解能把这三种成分在数学上分离开来。GRACE序列也是类似的混合体长期趋势如地下水持续减少、年周期降水补给和蒸散发的季节变化、半年周期、随机噪声它们各自的规律性不同SSA就按这种规律性强弱把它们排序。数学过程可以拆成四步第一步嵌入。给定序列x(t)长度N选择窗口长度L1 L N构造轨迹矩阵X。矩阵的每一列是长度为L的滑动窗口片段总共有K N - L 1列所以X是L×K的矩阵X [x(1) x(2) ... x(K) x(2) x(3) ... x(K1) ... x(L) x(L1) ... 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序列L36已经占序列长度的五分之一到四分之一用它做滞后矩阵的列数还有一百多足够完成SVD。如果序列更长比如2002到2024年积累到260个月L48也是可以考虑的但收益不明显反而让重构端点误差更宽。还有个实用技巧如果只想捕捉年周期而忽略半年周期L取24就够如果还想把半年周期也稳定分离L≥36更合适。因为窗口长度必须大于最大目标周期的两倍才能保证基频和谐波不在SVD里发生严重混叠。这是我在对比很多组实验后得出的结论。3.3 如何识别信号分量特征值谱与累计方差贡献率选定L后怎么决定保留几个主成分也就是d的取值最直观的工具是特征值谱图。把奇异值平方等价于特征值按降序画成折线图然后看拐点在哪里特征值从某个位置开始变得很小且下降平缓这个位置就是信号和噪声的分界。GRACE序列的特征值谱通常长这样第1个特征值特别大对应长期趋势第2、3个特征值一组代表年周期第4、5个一组代表半年周期从第6或第7个开始特征值大幅缩小进入缓慢衰减的长尾这部分基本就是噪声。另一个判断标准是累计方差贡献率。保留前d个主成分后它们对轨迹矩阵总方差的贡献比例用公式表示就是前d个奇异值平方之和除以全部奇异值平方之和。对GRACE典型序列前5-6个主成分通常能解释85%-95%的方差剩余的都是噪声和局地异常。你要是保留太多成分比如d20那噪声也被当成信号重建回去了插值效果会恶化保留太少又可能把半年周期甚至部分季节细节丢掉插值结果会过于平滑。我看过一些新手的做法是把d固定成5然后所有格点、所有月份都套用这是不推荐的。因为不同格点的时间序列特征差异很大干旱区序列噪声小、趋势和年周期干净而某些季节积水区、冰盖边缘区信号复杂度高d5反而不够。稳妥做法是每根时间序列先做个快速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:iL-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这种有大量真实观测月份的序列最实用的验证方法是挖洞测试取出原始时间序列中所有有效月份。随机挖走其中10%的有效值标记为NaN。用同一套SSA参数做插值。把插值结果和挖走前的真实值比较计算RMSE、相关系数、偏差等指标。重复多次比如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插值L36d6重建计算与真值的RMSE。结果如下表情景PCHIP RMSE (cm)SSA插值 RMSE (cm)PCHIP相关系数SSA相关系数A随机缺失15个点0.420.380.9720.981B连续缺失6个月1.871.090.8360.933C两个长段缺失2.611.370.7540.902从这个模拟结果能明显看出随机零散缺失时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.2cmPCHIP是1.8cm。由于该时段本身包含一次强降水恢复信号PCHIP把恢复过程拉平了SSA则较好地保留了恢复幅度。这正是真实序列中值得关注的差异——那段时间恰好对研究干旱恢复有重要意义插值方法选错了结论会有实质差异。当然真实数据插值永远不等于真实观测。对SSA补出的重要信号比如某年汛期水位是否恢复我都会去查阅该区域的地面水文站数据或再分析降水数据做一次旁证再写进论文。5.3 SSA插值容易翻车的四个细节第一个是边缘效应。SSA重构在序列首尾两端误差较大因为末尾段的窗口片段数量少SVD分解对这段的约束不足。如果缺失发生在序列末端插值结果不确定性很大。我的处理办法是把插值结果的两端各L个月单独标记在做趋势分析时不对这两年过度解读。如果非要用末端插值结果建议用较短窗口L24重新做一个敏感性测试。第二个是主成分数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结果都画出来对比再决定要不要对整片区域批量处理。代码、参数都是可以复现的但数据里那些说不清道不明的异常终究需要你亲眼盯着图看一遍。本文还有配套的精品资源点击获取
简介本资源是一套面向地球物理与水文研究者的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插值的思路和局部插值完全不同。它先假设真实信号可以表示为少数几个主要成分的叠加趋势逐年周期半年周期少量噪声然后通过矩阵分解把时间序列拆成这些成分再用主要成分去估计缺失位置的值。具体到缺失填充比较成熟的做法是迭代式填充流程很简单先用线性或PCHIP给缺失点补一个初始值得到完整序列。对完整序列做SSA分解得到前d个主要成分的重构序列。用重构序列去更新缺失点的旧值。重复进行SSA分解和更新直到缺失点数值变化小于阈值。这个流程执行下来缺失值会逐步收敛到一个与整体信号结构自洽的状态。因为SSA重建过程中非缺失点的信息也参与了每个主成分的估计相当于用全序列的数据来约束缺失段比只用缺口边界局部的信息去猜要稳健得多。我后面会给完整可跑的MATLAB代码。2.3 进阶路线时空联合插值DINEOF等如果研究区域不是一个格点而是整个网格场还可以用DINEOFData Interpolating Empirical Orthogonal Functions这类方法。它的原理是利用时空场的经验正交函数结构来填充数据本质上可以看成SSA在高维空间里的扩展。DINEOF在海洋遥感数据SST、叶绿素里用得很多插值同时能保持空间连续性对大面积云覆盖导致的缺测特别好用。但对GRACE月度网格来说有一个现实问题如果是某个月份整个全球数据都缺失DINEOF也一样无从下手因为它需要利用已有的时间协方差结构。而且DINEOF的调参过程判断截断模态数、迭代收敛标准比SSA更复杂对刚接触数据处理的同学没那么友好。我通常的建议是先用SSA把单点时间序列补好如果后续要做空间场分析再考虑DINEOF这种更重的工具。2.4 到底怎么选一张选型表我把各种情况下的推荐方案总结如下。场景推荐方案理由缺1-2个孤立月份PCHIP插值简单、稳健、不会明显改变谱结构连续缺3-6个月迭代SSA插值L24或36能利用全序列趋势和周期信息降低振幅低估连续缺半年以上迭代SSA插值NCEP/ERA5辅助数据交叉验证长缺口本身不确定性大必须做敏感性分析整场网格有空间孔洞DINEOF或MSSA多变量SSA利用时空协方差保持空间相关性序列末端缺失SSA插值要格外小心重构在两端有较大边缘误差需结合外推或降低权重表格之外还有一点需要强调无论选哪种插值最终都要做交叉验证把已知的完整月份故意挖去用剩余数据重建再与真值比较。不做验证的插值参数随便一调都会看起来很漂亮但真实误差你心里没底。3. 奇异谱分析(SSA)的原理与参数选法3.1 SSA的一次完整计算过程SSA的核心思想是找一个合适的窗口把一维时间序列转成一个轨迹矩阵然后对这个矩阵做奇异值分解再用特征值较大的那几个分量重建原序列。我用一个日常类比来帮助理解。假设你有一盘混合音频录音里面有人的说话声、背景音乐和电流底噪。SSA做的事情就是把这盘录音切成许多等长的片段窗口然后分析这些片段中重复出现的模式——说话声和音乐是有规律的底噪是随机的。通过矩阵分解能把这三种成分在数学上分离开来。GRACE序列也是类似的混合体长期趋势如地下水持续减少、年周期降水补给和蒸散发的季节变化、半年周期、随机噪声它们各自的规律性不同SSA就按这种规律性强弱把它们排序。数学过程可以拆成四步第一步嵌入。给定序列x(t)长度N选择窗口长度L1 L N构造轨迹矩阵X。矩阵的每一列是长度为L的滑动窗口片段总共有K N - L 1列所以X是L×K的矩阵X [x(1) x(2) ... x(K) x(2) x(3) ... x(K1) ... x(L) x(L1) ... 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序列L36已经占序列长度的五分之一到四分之一用它做滞后矩阵的列数还有一百多足够完成SVD。如果序列更长比如2002到2024年积累到260个月L48也是可以考虑的但收益不明显反而让重构端点误差更宽。还有个实用技巧如果只想捕捉年周期而忽略半年周期L取24就够如果还想把半年周期也稳定分离L≥36更合适。因为窗口长度必须大于最大目标周期的两倍才能保证基频和谐波不在SVD里发生严重混叠。这是我在对比很多组实验后得出的结论。3.3 如何识别信号分量特征值谱与累计方差贡献率选定L后怎么决定保留几个主成分也就是d的取值最直观的工具是特征值谱图。把奇异值平方等价于特征值按降序画成折线图然后看拐点在哪里特征值从某个位置开始变得很小且下降平缓这个位置就是信号和噪声的分界。GRACE序列的特征值谱通常长这样第1个特征值特别大对应长期趋势第2、3个特征值一组代表年周期第4、5个一组代表半年周期从第6或第7个开始特征值大幅缩小进入缓慢衰减的长尾这部分基本就是噪声。另一个判断标准是累计方差贡献率。保留前d个主成分后它们对轨迹矩阵总方差的贡献比例用公式表示就是前d个奇异值平方之和除以全部奇异值平方之和。对GRACE典型序列前5-6个主成分通常能解释85%-95%的方差剩余的都是噪声和局地异常。你要是保留太多成分比如d20那噪声也被当成信号重建回去了插值效果会恶化保留太少又可能把半年周期甚至部分季节细节丢掉插值结果会过于平滑。我看过一些新手的做法是把d固定成5然后所有格点、所有月份都套用这是不推荐的。因为不同格点的时间序列特征差异很大干旱区序列噪声小、趋势和年周期干净而某些季节积水区、冰盖边缘区信号复杂度高d5反而不够。稳妥做法是每根时间序列先做个快速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:iL-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这种有大量真实观测月份的序列最实用的验证方法是挖洞测试取出原始时间序列中所有有效月份。随机挖走其中10%的有效值标记为NaN。用同一套SSA参数做插值。把插值结果和挖走前的真实值比较计算RMSE、相关系数、偏差等指标。重复多次比如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插值L36d6重建计算与真值的RMSE。结果如下表情景PCHIP RMSE (cm)SSA插值 RMSE (cm)PCHIP相关系数SSA相关系数A随机缺失15个点0.420.380.9720.981B连续缺失6个月1.871.090.8360.933C两个长段缺失2.611.370.7540.902从这个模拟结果能明显看出随机零散缺失时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.2cmPCHIP是1.8cm。由于该时段本身包含一次强降水恢复信号PCHIP把恢复过程拉平了SSA则较好地保留了恢复幅度。这正是真实序列中值得关注的差异——那段时间恰好对研究干旱恢复有重要意义插值方法选错了结论会有实质差异。当然真实数据插值永远不等于真实观测。对SSA补出的重要信号比如某年汛期水位是否恢复我都会去查阅该区域的地面水文站数据或再分析降水数据做一次旁证再写进论文。5.3 SSA插值容易翻车的四个细节第一个是边缘效应。SSA重构在序列首尾两端误差较大因为末尾段的窗口片段数量少SVD分解对这段的约束不足。如果缺失发生在序列末端插值结果不确定性很大。我的处理办法是把插值结果的两端各L个月单独标记在做趋势分析时不对这两年过度解读。如果非要用末端插值结果建议用较短窗口L24重新做一个敏感性测试。第二个是主成分数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结果都画出来对比再决定要不要对整片区域批量处理。代码、参数都是可以复现的但数据里那些说不清道不明的异常终究需要你亲眼盯着图看一遍。本文还有配套的精品资源点击获取
