简介这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的学习者和科研人员聚焦站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题提出小波多尺度分解与奇异谱分析SSA相结合的建模思路。资源包内含1个docx文件约442KB正文系统梳理了小波多分辨率分析的低频概貌与高频细节分离机制、SSA构造时滞矩阵与奇异值分解的四步流程以及两者优势互补提取周期性与长期性变化信息的原理并给出全球11个测站20年GPS垂向坐标序列的实验验证。读者可据此理解非线性位移振幅达1至2厘米时对ITRF框架精度的影响掌握从噪声中分离周年、半周年等有用分量的具体方法为坐标时间序列建模、残差改正与精度提升提供可借鉴的技术路线。目前已有186人学习。1. 小波多尺度分解遇上 SSAGNSS 站坐标时间序列里那 1~2 cm 的非线性抖动到底怎么拆做 GNSS 高精度数据处理的人迟早会撞上同一个问题ITRF 框架下基准站的历元坐标和速度场名义上已经到毫米级可你把自己站点的垂向坐标时间序列拉出来一看周年、半年、季节性的起伏叠在一起振幅轻松到 1~2 cm线性速度根本描述不了。这部分非线性运动不处理后续速度场估计、参考框架维持、地壳形变解释都会带着系统性偏差。这份文档给出的思路很直接先用小波多尺度分解把原始序列拆成低频概貌加多层高频细节再对每一层单独做奇异谱分析SSA按特征值贡献率截取前几阶重构最后把各层拟合结果叠回去。它面向的是手里已经有 IGS 站坐标时间序列、想把这套组合方法跑通并复现精度提升的从业者尤其是做垂向非线性建模和周期项提取的人。文档用全球 11 个测站 1999—2018 年近 20 年的周采样数据做了验证结论是相比纯 SSARMSE 和 MAE 各降了约 26.5% 和 25.5%。2. 小波多尺度分解dbN 小波怎么选、分解层数为什么卡在 3 层小波多分辨率分析的核心是把信号投影到一串子空间里每一层拆成低频近似和高频细节。文档里给的分解关系是 S A3 D3 D2 D1继续往下拆就把 A3 再分成 A4 和 D4整个过程几乎无损。放到 GNSS 坐标时间序列这个场景低频部分承载趋势项和主要周期项高频部分装的是随机项和短周期细节。这一步做得好不好直接决定后面 SSA 在每层上能不能把周期项干净地拎出来。2.1 小波基函数选型为什么是 dbN 而不是 Haar 或 symN文档表 2 把常用小波基函数的支撑长度、消失矩阶数、对称性和特点列了一遍。Haar 小波支撑长度 1、消失矩 1时域上不连续频率局部性差拿来做坐标序列基本是自找麻烦。symN 小波近似对称、能减少相位失真但它更偏图像处理领域。bior 小波不对称虽然线性相位性在信号重构里常用但对这种以周期提取为目标的序列不是首选。dbN 小波支撑长度 2N、消失矩 N 阶、近似对称光滑性随 N 增大而增强处理坐标时间序列优势明显。文档最终实验用的是 db4兼顾正交性和紧支撑性。选 db4 不是拍脑袋。db 系列的正交性保证分解后各层能量不串紧支撑性保证局部突变不会被抹到整条序列上。N 太小比如 db1 就是 Haar频率局部性差N 太大支撑变长、边界效应加重对只有 20 年周采样的序列不划算。我一般会先在 db2 到 db6 之间扫一遍看哪一档的细节层里周期成分最干净再定下来。2.2 分解层数2 到 6 层的 RMSE/MAE 实测对比分解层数不是越多越好。文档表 3 给了 BJFS 站垂向序列在 2~6 层下的拟合精度分解层数RMSE/mmMAE/mm22.151.6831.881.4941.801.4351.801.4461.801.44从 3 层往上精度基本不再变化4 层和 3 层差异很小。但层数增加会带来计算误差累积和处理时间上升坐标时间序列数据量本来就大效率必须考虑。文档的结论是分解层数一般选 3 层。这个判断我认同3 层已经把趋势和主要周期分到了 a3 和 d3d1、d2 里随机项占主导再往下拆收益递减。2.3 三层分解的落地步骤与参数按文档流程对原始序列 S 做 db4 三层分解得到 a3、d3、d2、d1。用 Python 的 PyWavelets 复现大致是这样import pywt import numpy as np # series: 一维 GNSS 垂向坐标时间序列采样间隔一周 # waveletdb4mode 用 periodization 处理边界避免端点失真 coeffs pywt.wavedec(series, waveletdb4, modeperiodization, level3) # coeffs 顺序: [a3, d3, d2, d1] a3, d3, d2, d1 coeffs # 各层单独重构回原始长度便于后续逐层 SSA def reconstruct_one(coeffs, idx, waveletdb4, modeperiodization): tmp [np.zeros_like(c) for c in coeffs] tmp[idx] coeffs[idx] return pywt.waverec(tmp, waveletwavelet, modemode) a3_rec reconstruct_one(coeffs, 0) d3_rec reconstruct_one(coeffs, 1) d2_rec reconstruct_one(coeffs, 2) d1_rec reconstruct_one(coeffs, 3) # 验证无损性: a3_rec d3_rec d2_rec d1_rec 应约等于 series resid series - (a3_rec d3_rec d2_rec d1_rec) print(重构误差 max:, np.max(np.abs(resid)))逻辑说明wavedec返回的是各层系数直接拿来做 SSA 长度对不上所以用reconstruct_one把每一层单独重构回原始长度。modeperiodization是关键参数坐标序列首尾不连续用默认的对称延拓会在两端引入伪影周期化模式对周期信号更稳。重构误差那行是自检正常应该在 1e-10 量级如果明显偏大说明层数或模式选错了。参数说明level3对应文档结论waveletdb4对应文档实验选择mode可按数据端点质量在periodization和symmetric之间试端点干净优先前者。3. SSA 逐层重构窗口长度 L 取 52、重构阶次 K 怎么定SSA 的四步——构造时滞矩阵、SVD、分组、对角平均化——文档写得很完整。落到实操真正卡人的是两个参数窗口长度 L 和重构阶次 K。文档明确 L 不宜超过数据长度 N 的 1/3若有先验周期则取周期的公倍数。周采样下已知周年和半年周期最小公倍数 52所以 L 取 52。K 靠奇异值贡献率定太小会把信号当噪声剔掉太大又把噪声当信号提出来。3.1 时滞矩阵与 SVD 的实现要点对长度为 N 的序列 x取窗口 LK N - L 1构造 L×K 的 Hankel 矩阵 X副对角线元素相等。然后对 X 做 SVD得到奇异值从大到小排列的奇异谱。文档式 (6) 给了贡献率定义某一组特征值之和除以全部特征值之和。这一步实现上没什么玄学但矩阵规模要注意——N 是 20 年周采样约 1040 个点L52 时 X 是 52×989SVD 秒级完成不用担心。def ssa_decompose(x, L): N len(x) K N - L 1 # 构造轨迹矩阵行 i 列 j 对应 x[ij] X np.column_stack([x[j:jL] for j in range(K)]) # SVD: X U diag(s) Vt U, s, Vt np.linalg.svd(X, full_matricesFalse) # 贡献率 contrib s**2 / np.sum(s**2) return U, s, Vt, contrib def ssa_reconstruct(U, s, Vt, idxs): # 按选定分量重建轨迹矩阵 Xr (U[:, idxs] * s[idxs]) Vt[idxs, :] L, K Xr.shape N L K - 1 # 对角平均化 out np.zeros(N) cnt np.zeros(N) for i in range(L): for j in range(K): out[ij] Xr[i, j] cnt[ij] 1 return out / cnt逻辑说明ssa_decompose里np.column_stack构造的就是文档式 (1) 的 Hankel 矩阵第 j 列是 x[j] 到 x[jL-1]。SVD 返回的 s 平方归一化就是贡献率。ssa_reconstruct里对角平均化对应文档式 (7)把矩阵沿副对角线取平均还原成一维序列。选 idxs 就是选前 K 阶。参数说明L52 是文档给定值idxs 传range(K)即取前 K 阶。注意full_matricesFalse能省大量内存L 和 K 差很大时尤其明显。3.2 贡献率定阶BJFS 站前 14 阶的实测分布文档表 1 给了 BJFS 站垂向序列前 14 阶贡献率这是定 K 的直接依据阶次贡献率/%阶次贡献率/%132.3081.20230.4591.1334.71101.0644.68110.9354.19120.8361.55130.7771.36140.77RRC1 和 RRC2 贡献率接近且都很大是一对近似相等的特征值对应同周期同振幅的分量。RRC3、RRC4、RRC5 次之。从第 6 阶开始贡献率掉到 1.5% 以下。文档对前 6 阶做 FFT 后发现RRC1RRC2 合并是周期 1 a、振幅 4.76 mm 的周年项RRC5 分别与 RRC3、RRC4 对应 0.5 a振幅 0.95 mm和 9 a振幅 0.88 mm周期项RRC6 是 0.3 a 的季节项且振幅很小。所以趋势项和主要周期项集中在前 5 阶K 取 5。这里有个容易翻车的点近似相等的特征值必须成对处理。RRC1 和 RRC2 单独看都是混频的合并才是完整的周年项。判断方法就是看贡献率是否接近、FFT 主频是否一致。3.3 逐层 SSA 与叠加拿到 a3、d3、d2、d1 后对每一层分别跑 SSA各自取贡献率大的前 5 阶重构再把四层的拟合值相加得到最终拟合序列。文档图 8~图 11 显示a3 和 d3 的重构序列与原始层序列大部分时段重合d1、d2 随机项多、周期项少各阶贡献率差异不大拟合效果弱于其他层。最终叠加结果残差振幅从纯 SSA 的约 3 mm 降到约 2 mm。layers [a3_rec, d3_rec, d2_rec, d1_rec] fitted np.zeros_like(series) for layer in layers: U, s, Vt, contrib ssa_decompose(layer, L52) # 取前 5 阶重构K 可按各层贡献率拐点微调 rec ssa_reconstruct(U, s, Vt, idxslist(range(5))) fitted rec residual series - fitted rmse np.sqrt(np.mean(residual**2)) mae np.mean(np.abs(residual)) print(fRMSE{rmse:.2f} mm, MAE{mae:.2f} mm)逻辑说明每层独立定阶是这套方法的关键文档结语第 3 条特别强调不宜对所有测站取相同特征值个数。代码里 K 固定为 5 只是 BJFS 站的取值换站要重新看贡献率拐点。fitted是四层重构之和对应文档式 (12)。参数说明L 保持 52idxs 的阶数按各层贡献率曲线拐点定a3、d3 通常 5 阶够d1、d2 可适当增减。RMSE/MAE 按文档式 (13) 计算。4. 避坑与排查这套组合方法最容易翻车的五个地方4.1 现象重构序列和原始序列对不上误差远超预期原因小波分解用了默认边界延拓模式坐标序列首尾不连续对称延拓在两端造出伪影重构时误差被放大。解决改用modeperiodization并在分解后立刻做一次无损性自检a3d3d2d1与原始序列的差应在 1e-10 量级偏大就换模式重来。4.2 现象SSA 拟合残差里还残留明显周期原因K 取小了短周期项被当噪声剔除。文档图 6 显示纯 SSA 残差频谱里还有半年以下周期。解决先对残差做 FFT看残留主频落在哪回到贡献率表把对应阶次纳入重构或者直接上小波SSA 组合让短周期在细节层里被单独处理。4.3 现象换了测站精度不升反降原因照搬了 BJFS 站的 K5。文档结语明确说各测站前 K 阶贡献率大小存在差异地理位置不同规律不同。解决每个站单独画贡献率曲线找拐点定 K低纬和高纬站的周年、半年振幅差异明显定阶不能一刀切。4.4 现象分解层数加到 5、6 层精度没变但跑得越来越慢原因文档表 3 已经说明 3 层以上精度基本不变多出来的层只增加计算量和误差累积。解决层数卡 3 层把精力花在每层的 K 上收益比加层大得多。4.5 现象近似相等的特征值被拆开单独重构周期项振幅不对原因RRC1 和 RRC2 这类成对特征值对应同一物理分量单独重构会混频。解决定阶时先看贡献率是否成对接近再用 FFT 验证主频是否一致一致就合并重构文档里周年项就是 RRC1RRC2 合并得到的。5. 进阶把 26% 的精度提升稳住以及跨站自适应定阶的一个习惯文档给出的 26.5% RMSE 降幅和 25.5% MAE 降幅是在 11 个测站上算出来的从 9.2°N 的 ADIS 到 67.4°S 的 MAW1纬度跨度覆盖低中高。表 4 里纯 SSA 的 RMSE 从 2.12 mm 到 2.94 mm小波SSA 降到 1.46 mm 到 2.23 mm没有一个站例外。这说明方法本身稳但前提是每站参数对。想把这份提升在自己数据上复现我一般会走这么一套验证流程。先固定小波侧db4、3 层、periodization这三项在文档里已经验证过不用反复调。然后对每一层单独画贡献率曲线找拐点。拐点判断有个土办法把相邻两阶贡献率做比值比值突然掉一个量级的位置就是拐点。BJFS 站是第 5 阶到第 6 阶从 4.19% 掉到 1.55%拐点清晰。如果曲线平滑没有明显拐点说明这层以随机项为主K 可以取小甚至只取 1~2 阶d1、d2 经常是这种情况。再就是对残差做 FFT 闭环验证。拟合完不算完把残差频谱拉出来看还有没有超过噪声底的周期峰。文档图 13 显示组合方法残差里季节、月周期影响明显减小这就是验收标准。如果残差里还有峰回到对应层的 K 上补阶。跨站自适应定阶这件事文档结语把它列为待探讨方向实操里我的习惯是建一张站-层-K 的对照表每处理一个新站先跑一遍贡献率统计把 K 记进去。同一区域、纬度相近的站 K 往往接近可以拿来当初始值但最终仍以本站贡献率拐点为准。这套流程走下来26% 的提升不是碰运气是可复现的。从那以后我每次拿到新站的坐标时间序列都强制先跑一遍小波无损性自检和逐层贡献率统计再动 SSA 的 K省得回头在残差里找周期找到怀疑人生。希望帮到你。本文还有配套的精品资源点击获取 小波多尺度分解结合SSA：GNSS坐标时间序列非线性抖动拆解与精度提升
简介这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的学习者和科研人员聚焦站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题提出小波多尺度分解与奇异谱分析SSA相结合的建模思路。资源包内含1个docx文件约442KB正文系统梳理了小波多分辨率分析的低频概貌与高频细节分离机制、SSA构造时滞矩阵与奇异值分解的四步流程以及两者优势互补提取周期性与长期性变化信息的原理并给出全球11个测站20年GPS垂向坐标序列的实验验证。读者可据此理解非线性位移振幅达1至2厘米时对ITRF框架精度的影响掌握从噪声中分离周年、半周年等有用分量的具体方法为坐标时间序列建模、残差改正与精度提升提供可借鉴的技术路线。目前已有186人学习。1. 小波多尺度分解遇上 SSAGNSS 站坐标时间序列里那 1~2 cm 的非线性抖动到底怎么拆做 GNSS 高精度数据处理的人迟早会撞上同一个问题ITRF 框架下基准站的历元坐标和速度场名义上已经到毫米级可你把自己站点的垂向坐标时间序列拉出来一看周年、半年、季节性的起伏叠在一起振幅轻松到 1~2 cm线性速度根本描述不了。这部分非线性运动不处理后续速度场估计、参考框架维持、地壳形变解释都会带着系统性偏差。这份文档给出的思路很直接先用小波多尺度分解把原始序列拆成低频概貌加多层高频细节再对每一层单独做奇异谱分析SSA按特征值贡献率截取前几阶重构最后把各层拟合结果叠回去。它面向的是手里已经有 IGS 站坐标时间序列、想把这套组合方法跑通并复现精度提升的从业者尤其是做垂向非线性建模和周期项提取的人。文档用全球 11 个测站 1999—2018 年近 20 年的周采样数据做了验证结论是相比纯 SSARMSE 和 MAE 各降了约 26.5% 和 25.5%。2. 小波多尺度分解dbN 小波怎么选、分解层数为什么卡在 3 层小波多分辨率分析的核心是把信号投影到一串子空间里每一层拆成低频近似和高频细节。文档里给的分解关系是 S A3 D3 D2 D1继续往下拆就把 A3 再分成 A4 和 D4整个过程几乎无损。放到 GNSS 坐标时间序列这个场景低频部分承载趋势项和主要周期项高频部分装的是随机项和短周期细节。这一步做得好不好直接决定后面 SSA 在每层上能不能把周期项干净地拎出来。2.1 小波基函数选型为什么是 dbN 而不是 Haar 或 symN文档表 2 把常用小波基函数的支撑长度、消失矩阶数、对称性和特点列了一遍。Haar 小波支撑长度 1、消失矩 1时域上不连续频率局部性差拿来做坐标序列基本是自找麻烦。symN 小波近似对称、能减少相位失真但它更偏图像处理领域。bior 小波不对称虽然线性相位性在信号重构里常用但对这种以周期提取为目标的序列不是首选。dbN 小波支撑长度 2N、消失矩 N 阶、近似对称光滑性随 N 增大而增强处理坐标时间序列优势明显。文档最终实验用的是 db4兼顾正交性和紧支撑性。选 db4 不是拍脑袋。db 系列的正交性保证分解后各层能量不串紧支撑性保证局部突变不会被抹到整条序列上。N 太小比如 db1 就是 Haar频率局部性差N 太大支撑变长、边界效应加重对只有 20 年周采样的序列不划算。我一般会先在 db2 到 db6 之间扫一遍看哪一档的细节层里周期成分最干净再定下来。2.2 分解层数2 到 6 层的 RMSE/MAE 实测对比分解层数不是越多越好。文档表 3 给了 BJFS 站垂向序列在 2~6 层下的拟合精度分解层数RMSE/mmMAE/mm22.151.6831.881.4941.801.4351.801.4461.801.44从 3 层往上精度基本不再变化4 层和 3 层差异很小。但层数增加会带来计算误差累积和处理时间上升坐标时间序列数据量本来就大效率必须考虑。文档的结论是分解层数一般选 3 层。这个判断我认同3 层已经把趋势和主要周期分到了 a3 和 d3d1、d2 里随机项占主导再往下拆收益递减。2.3 三层分解的落地步骤与参数按文档流程对原始序列 S 做 db4 三层分解得到 a3、d3、d2、d1。用 Python 的 PyWavelets 复现大致是这样import pywt import numpy as np # series: 一维 GNSS 垂向坐标时间序列采样间隔一周 # waveletdb4mode 用 periodization 处理边界避免端点失真 coeffs pywt.wavedec(series, waveletdb4, modeperiodization, level3) # coeffs 顺序: [a3, d3, d2, d1] a3, d3, d2, d1 coeffs # 各层单独重构回原始长度便于后续逐层 SSA def reconstruct_one(coeffs, idx, waveletdb4, modeperiodization): tmp [np.zeros_like(c) for c in coeffs] tmp[idx] coeffs[idx] return pywt.waverec(tmp, waveletwavelet, modemode) a3_rec reconstruct_one(coeffs, 0) d3_rec reconstruct_one(coeffs, 1) d2_rec reconstruct_one(coeffs, 2) d1_rec reconstruct_one(coeffs, 3) # 验证无损性: a3_rec d3_rec d2_rec d1_rec 应约等于 series resid series - (a3_rec d3_rec d2_rec d1_rec) print(重构误差 max:, np.max(np.abs(resid)))逻辑说明wavedec返回的是各层系数直接拿来做 SSA 长度对不上所以用reconstruct_one把每一层单独重构回原始长度。modeperiodization是关键参数坐标序列首尾不连续用默认的对称延拓会在两端引入伪影周期化模式对周期信号更稳。重构误差那行是自检正常应该在 1e-10 量级如果明显偏大说明层数或模式选错了。参数说明level3对应文档结论waveletdb4对应文档实验选择mode可按数据端点质量在periodization和symmetric之间试端点干净优先前者。3. SSA 逐层重构窗口长度 L 取 52、重构阶次 K 怎么定SSA 的四步——构造时滞矩阵、SVD、分组、对角平均化——文档写得很完整。落到实操真正卡人的是两个参数窗口长度 L 和重构阶次 K。文档明确 L 不宜超过数据长度 N 的 1/3若有先验周期则取周期的公倍数。周采样下已知周年和半年周期最小公倍数 52所以 L 取 52。K 靠奇异值贡献率定太小会把信号当噪声剔掉太大又把噪声当信号提出来。3.1 时滞矩阵与 SVD 的实现要点对长度为 N 的序列 x取窗口 LK N - L 1构造 L×K 的 Hankel 矩阵 X副对角线元素相等。然后对 X 做 SVD得到奇异值从大到小排列的奇异谱。文档式 (6) 给了贡献率定义某一组特征值之和除以全部特征值之和。这一步实现上没什么玄学但矩阵规模要注意——N 是 20 年周采样约 1040 个点L52 时 X 是 52×989SVD 秒级完成不用担心。def ssa_decompose(x, L): N len(x) K N - L 1 # 构造轨迹矩阵行 i 列 j 对应 x[ij] X np.column_stack([x[j:jL] for j in range(K)]) # SVD: X U diag(s) Vt U, s, Vt np.linalg.svd(X, full_matricesFalse) # 贡献率 contrib s**2 / np.sum(s**2) return U, s, Vt, contrib def ssa_reconstruct(U, s, Vt, idxs): # 按选定分量重建轨迹矩阵 Xr (U[:, idxs] * s[idxs]) Vt[idxs, :] L, K Xr.shape N L K - 1 # 对角平均化 out np.zeros(N) cnt np.zeros(N) for i in range(L): for j in range(K): out[ij] Xr[i, j] cnt[ij] 1 return out / cnt逻辑说明ssa_decompose里np.column_stack构造的就是文档式 (1) 的 Hankel 矩阵第 j 列是 x[j] 到 x[jL-1]。SVD 返回的 s 平方归一化就是贡献率。ssa_reconstruct里对角平均化对应文档式 (7)把矩阵沿副对角线取平均还原成一维序列。选 idxs 就是选前 K 阶。参数说明L52 是文档给定值idxs 传range(K)即取前 K 阶。注意full_matricesFalse能省大量内存L 和 K 差很大时尤其明显。3.2 贡献率定阶BJFS 站前 14 阶的实测分布文档表 1 给了 BJFS 站垂向序列前 14 阶贡献率这是定 K 的直接依据阶次贡献率/%阶次贡献率/%132.3081.20230.4591.1334.71101.0644.68110.9354.19120.8361.55130.7771.36140.77RRC1 和 RRC2 贡献率接近且都很大是一对近似相等的特征值对应同周期同振幅的分量。RRC3、RRC4、RRC5 次之。从第 6 阶开始贡献率掉到 1.5% 以下。文档对前 6 阶做 FFT 后发现RRC1RRC2 合并是周期 1 a、振幅 4.76 mm 的周年项RRC5 分别与 RRC3、RRC4 对应 0.5 a振幅 0.95 mm和 9 a振幅 0.88 mm周期项RRC6 是 0.3 a 的季节项且振幅很小。所以趋势项和主要周期项集中在前 5 阶K 取 5。这里有个容易翻车的点近似相等的特征值必须成对处理。RRC1 和 RRC2 单独看都是混频的合并才是完整的周年项。判断方法就是看贡献率是否接近、FFT 主频是否一致。3.3 逐层 SSA 与叠加拿到 a3、d3、d2、d1 后对每一层分别跑 SSA各自取贡献率大的前 5 阶重构再把四层的拟合值相加得到最终拟合序列。文档图 8~图 11 显示a3 和 d3 的重构序列与原始层序列大部分时段重合d1、d2 随机项多、周期项少各阶贡献率差异不大拟合效果弱于其他层。最终叠加结果残差振幅从纯 SSA 的约 3 mm 降到约 2 mm。layers [a3_rec, d3_rec, d2_rec, d1_rec] fitted np.zeros_like(series) for layer in layers: U, s, Vt, contrib ssa_decompose(layer, L52) # 取前 5 阶重构K 可按各层贡献率拐点微调 rec ssa_reconstruct(U, s, Vt, idxslist(range(5))) fitted rec residual series - fitted rmse np.sqrt(np.mean(residual**2)) mae np.mean(np.abs(residual)) print(fRMSE{rmse:.2f} mm, MAE{mae:.2f} mm)逻辑说明每层独立定阶是这套方法的关键文档结语第 3 条特别强调不宜对所有测站取相同特征值个数。代码里 K 固定为 5 只是 BJFS 站的取值换站要重新看贡献率拐点。fitted是四层重构之和对应文档式 (12)。参数说明L 保持 52idxs 的阶数按各层贡献率曲线拐点定a3、d3 通常 5 阶够d1、d2 可适当增减。RMSE/MAE 按文档式 (13) 计算。4. 避坑与排查这套组合方法最容易翻车的五个地方4.1 现象重构序列和原始序列对不上误差远超预期原因小波分解用了默认边界延拓模式坐标序列首尾不连续对称延拓在两端造出伪影重构时误差被放大。解决改用modeperiodization并在分解后立刻做一次无损性自检a3d3d2d1与原始序列的差应在 1e-10 量级偏大就换模式重来。4.2 现象SSA 拟合残差里还残留明显周期原因K 取小了短周期项被当噪声剔除。文档图 6 显示纯 SSA 残差频谱里还有半年以下周期。解决先对残差做 FFT看残留主频落在哪回到贡献率表把对应阶次纳入重构或者直接上小波SSA 组合让短周期在细节层里被单独处理。4.3 现象换了测站精度不升反降原因照搬了 BJFS 站的 K5。文档结语明确说各测站前 K 阶贡献率大小存在差异地理位置不同规律不同。解决每个站单独画贡献率曲线找拐点定 K低纬和高纬站的周年、半年振幅差异明显定阶不能一刀切。4.4 现象分解层数加到 5、6 层精度没变但跑得越来越慢原因文档表 3 已经说明 3 层以上精度基本不变多出来的层只增加计算量和误差累积。解决层数卡 3 层把精力花在每层的 K 上收益比加层大得多。4.5 现象近似相等的特征值被拆开单独重构周期项振幅不对原因RRC1 和 RRC2 这类成对特征值对应同一物理分量单独重构会混频。解决定阶时先看贡献率是否成对接近再用 FFT 验证主频是否一致一致就合并重构文档里周年项就是 RRC1RRC2 合并得到的。5. 进阶把 26% 的精度提升稳住以及跨站自适应定阶的一个习惯文档给出的 26.5% RMSE 降幅和 25.5% MAE 降幅是在 11 个测站上算出来的从 9.2°N 的 ADIS 到 67.4°S 的 MAW1纬度跨度覆盖低中高。表 4 里纯 SSA 的 RMSE 从 2.12 mm 到 2.94 mm小波SSA 降到 1.46 mm 到 2.23 mm没有一个站例外。这说明方法本身稳但前提是每站参数对。想把这份提升在自己数据上复现我一般会走这么一套验证流程。先固定小波侧db4、3 层、periodization这三项在文档里已经验证过不用反复调。然后对每一层单独画贡献率曲线找拐点。拐点判断有个土办法把相邻两阶贡献率做比值比值突然掉一个量级的位置就是拐点。BJFS 站是第 5 阶到第 6 阶从 4.19% 掉到 1.55%拐点清晰。如果曲线平滑没有明显拐点说明这层以随机项为主K 可以取小甚至只取 1~2 阶d1、d2 经常是这种情况。再就是对残差做 FFT 闭环验证。拟合完不算完把残差频谱拉出来看还有没有超过噪声底的周期峰。文档图 13 显示组合方法残差里季节、月周期影响明显减小这就是验收标准。如果残差里还有峰回到对应层的 K 上补阶。跨站自适应定阶这件事文档结语把它列为待探讨方向实操里我的习惯是建一张站-层-K 的对照表每处理一个新站先跑一遍贡献率统计把 K 记进去。同一区域、纬度相近的站 K 往往接近可以拿来当初始值但最终仍以本站贡献率拐点为准。这套流程走下来26% 的提升不是碰运气是可复现的。从那以后我每次拿到新站的坐标时间序列都强制先跑一遍小波无损性自检和逐层贡献率统计再动 SSA 的 K省得回头在残差里找周期找到怀疑人生。希望帮到你。本文还有配套的精品资源点击获取