简介这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的学习者与科研人员聚焦站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题提出小波多尺度分解与奇异谱分析SSA相结合的分析思路。资源包内含1个docx文件约442KB正文系统梳理了小波多分辨率分解的低频概貌与高频细节分离机制、SSA构造时滞矩阵与奇异值分解的完整流程以及两者优势互补的建模逻辑并给出全球11个测站20年GPS垂向坐标序列的实验验证。读者可从中获取非线性变化建模的完整方法框架、算法公式推导与实验分析结论理解如何从噪声中提取周年、半周年等周期信号进而提升坐标时间序列精度。目前已有186人学习适合希望深入掌握GNSS坐标时序分析技术的中高级读者参考。1. 从一条“抖得不像话”的坐标序列说起GNSS 站坐标时间序列分析里最让人头疼的不是缺数据而是数据看起来“什么都有”趋势、周年、半周年、构造形变、天线相位中心变化、多路径、未模型化的轨道误差全糊在一条曲线上。你直接拿它去做线性拟合残差里全是结构你直接拿它去做频谱主频又被非平稳性抹开。小波多尺度分解和奇异谱分析SSA之所以常被放在一起用是因为它们解决的是同一件事的两个侧面小波负责把非平稳信号按尺度拆开SSA 负责在每一层里把“可重建的振荡分量”从噪声里拎出来。GNSS 坐标时间序列的常见做法是先用小波把高频多路径和低频趋势分离再对中间频段做 SSA 重构最后把各分量叠回去做速度场或形变解释。这套流程适合两类人一类是做地壳形变、沉降监测、站速度估计的从业者另一类是被“坐标序列预测”或“gnss时间序列预测”需求推着走、但发现 ARIMA 一上就翻车的工程师。下面按“先立住概念、再跑通最小流程、最后讲坑”的顺序展开。2. 小波多尺度分解在 GNSS 坐标序列里到底拆什么2.1 为什么不能直接对原始序列做 SSASSA 的核心是把一维序列嵌入成轨迹矩阵再做奇异值分解按奇异值大小分组重构。它对“近似平稳、由若干振荡加噪声组成”的序列效果很好。但 GNSS 站坐标序列有三个不平稳来源长期线性趋势、阶跃天线更换、地震同震、以及方差随尺度变化的噪声。如果你把原始序列直接丢给 SSA趋势会占据第一个奇异值对把真正的周年项挤到后面阶跃会被当成一个超强“振荡”重构出来就是一条假周期。常见做法是先做一阶差分或去线性趋势但差分会放大高频噪声去趋势又可能把低频构造信号一起削掉。小波多尺度分解的价值就在这里它用一组不同尺度的基函数把序列展开你可以选择在哪些尺度上做 SSA而不是一刀切。2.2 小波基、分解层数与边界延拓的选型工程上最常用的是离散小波变换DWT基函数选 sym8 或 db4。sym8 的对称性较好重构时相位失真小适合坐标序列这种需要保留波形形状的场景db4 更短计算快但边界效应更明显。分解层数按采样率和目标频段定GNSS 日解序列采样间隔 1 天Nyquist 频率 0.5 cycle/day周年项在 1/365 cycle/day。如果你关心周年和半周年分解到第 6 层左右细节层 d4~d6 大致覆盖 16~128 天周期逼近层 a6 就是更长周期的趋势。边界延拓必须做否则两端会出现大幅摆动常见用对称延拓symmetric或周期延拓periodic。坐标序列不是严格周期的对称延拓更稳。2.3 用 PyWavelets 跑通一次多尺度分解import numpy as np import pywt import matplotlib.pyplot as plt # 假设 data 是长度 N 的 GNSS 日解坐标序列已去掉明显粗差 # 采样间隔 1 天单位毫米 N len(data) wavelet sym8 level 6 # 对称延拓减少边界效应 coeffs pywt.wavedec(data, wavelet, modesymmetric, levellevel) # coeffs[0] 是逼近系数 a6coeffs[1] 是 d6... coeffs[6] 是 d1 # 重构各层细节便于逐层分析 details [] for i in range(1, level 1): c [np.zeros_like(coeffs[j]) for j in range(len(coeffs))] c[i] coeffs[i] details.append(pywt.waverec(c, wavelet, modesymmetric)[:N]) approx pywt.waverec([coeffs[0]] [np.zeros_like(coeffs[j]) for j in range(1, len(coeffs))], wavelet, modesymmetric)[:N]这段代码的逻辑是wavedec把序列拆成一层逼近系数加多层细节系数重构时只保留某一层系数、其余置零就能得到该尺度上的分量。参数上modesymmetric控制边界延拓方式level6决定最粗尺度。跑完后先看approx是否还残留周年波动——如果周年项跑到逼近层里说明层数不够加到 7 或 8。再看details[-1]d1是不是几乎全是白噪声如果是说明高频层可以整层丢弃不必送进 SSA。2.4 分解后怎么判断哪一层该进 SSA不是所有层都值得做 SSA。d1、d2 通常以多路径和观测噪声为主SSA 重构出来的“周期”多半是噪声的随机起伏强行解释就是玄学。a6 或 a7 是趋势和长周期项SSA 对趋势不敏感做了也白做。真正值得做 SSA 的是中间层比如 d4~d6这里往往混着周年、半周年、以及一些区域性的非构造信号。判断方法很简单对每一层细节做自相关如果自相关在滞后 365 天附近有明显峰值说明周年项主要落在这一层如果自相关快速衰减到零这层就是噪声主导跳过。这一步不做后面 SSA 的窗口长度和分组数就没有依据。3. 奇异谱分析窗口长度、分组与重构的工程参数3.1 轨迹矩阵的构造与窗口长度 L 的取值SSA 第一步是嵌入给定窗口长度 L把长度 N 的序列构造成 L×(K) 的轨迹矩阵KN-L1。L 的选择直接决定你能分离出的周期下限和上限。经验规则是 L 取目标周期长度的 1~2 倍。如果你要分离周年项L 至少取 365最好 500~730。L 太小周年和半周年在奇异值谱上会挤在一起分不开L 太大计算量上去而且趋势会污染更多分量。对 GNSS 日解序列我一般先取 L365 跑一遍看奇异值谱如果前两对奇异值对应的重构分量明显是周年就固定如果分不开加到 547 或 730。注意 K 必须大于 L否则轨迹矩阵秩不够所以序列长度至少要有 2L 以上做周年分析时序列最好 3 年以上。3.2 奇异值分组怎么把“成对”的分量认出来SSA 重构时一个实振荡分量通常对应两个相邻的奇异值因为正弦和余弦成对出现它们的重构分量频率相同、相位差 90 度。分组就是把这些成对的奇异值归到一组再重构。工程上不要只看奇异值大小要看重构分量的波形和频谱。具体做法对每个奇异值单独重构一个分量画出来做 FFT。如果第 i 和第 i1 个分量的主频一致、振幅接近就归为一组。周年项一般落在第 2、3 个奇异值第 1 个常是趋势半周年在第 4、5 个附近。如果第 1 个奇异值重构出来是趋势而不是振荡说明去趋势没做干净回到小波那一步把 a6 去掉再进 SSA。3.3 用 Python 实现 SSA 并重构周年分量import numpy as np def ssa_decompose(x, L): N len(x) K N - L 1 # 构造轨迹矩阵 X np.column_stack([x[i:iL] for i in range(K)]) # SVD U, s, Vt np.linalg.svd(X, full_matricesFalse) return U, s, Vt, L, K def ssa_reconstruct(U, s, Vt, L, K, groups): # groups: list of lists, 每个子列表是一组奇异值索引 N L K - 1 rec np.zeros(N) for g in groups: Xg np.zeros((L, K)) for i in g: Xg s[i] * np.outer(U[:, i], Vt[i, :]) # 反对角平均还原一维序列 for k in range(N): idx [j for j in range(K) if 0 k - j L] rec[k] np.mean([Xg[k - j, j] for j in idx]) return rec # 对某一层细节 d 做 SSA U, s, Vt, L, K ssa_decompose(d, L365) # 先看前 10 个奇异值 print(s[:10]) # 假设周年在第 2、3 个索引 1、2半周年在第 4、5 个索引 3、4 annual ssa_reconstruct(U, s, Vt, L, K, [[1, 2]]) semiannual ssa_reconstruct(U, s, Vt, L, K, [[3, 4]])ssa_decompose里L365是窗口长度K自动算出。ssa_reconstruct的groups参数是分组索引注意 Python 从 0 开始所以第 2、3 个奇异值对应索引 1、2。反对角平均那一步是 SSA 还原一维序列的标准操作不能省否则重构序列会有相位偏移。跑完后把annual和semiannual叠回小波分解的其他层就得到去噪后的坐标序列。参数上如果s[:10]里前两个奇异值远大于后面说明趋势没去干净如果第 2、3 个和第 4、5 个大小接近说明周年和半周年能量相当分组时要小心别把半周年并进周年。3.4 重构分量叠回去之前要检查什么叠回去之前至少做三件事。第一检查重构的周年分量振幅和相位是否逐年稳定。GNSS 坐标序列的周年项振幅通常几毫米如果某一年突然跳到十几毫米多半是那段时间有未模型化的阶跃或数据中断需要单独处理。第二检查残差序列的 RMS 是否比原始序列明显下降但不要追求降得越低越好——降得太狠说明你把噪声也当信号重构了残差里会留下明显的白噪声特征反而说明过拟合。第三把重构分量和原始序列叠画看周年峰值位置是否对齐。如果错开半个月以上检查小波重构时的边界延拓和 SSA 的反对角平均是否一致。4. 避坑与排查这套流程最容易翻车的五个地方4.1 现象重构后周年项振幅逐年漂移原因小波分解时用了periodic延拓而坐标序列两端并不周期边界系数被强行扭曲重构后误差向内部传播。解决改用symmetric或zero延拓并在分解前把序列两端各截掉 30 天重构后再补回避免边界污染进入 SSA 窗口。4.2 现象SSA 奇异值谱里前两个分量频率相同但相位差不是 90 度原因序列里有未去除的阶跃阶跃在 SSA 里表现为一个低频强分量和趋势混在一起破坏了振荡分量的成对性。解决进 SSA 之前先做阶跃检测常用方法是对一阶差分做中位数绝对偏差检验超过阈值的点标记为阶跃分段去均值后再拼接。4.3 现象窗口长度 L 取 365 时周年和半周年在奇异值谱上分不开原因L 刚好等于周年周期轨迹矩阵对周年项的响应最强但半周年周期是 182.5 天L365 时它的嵌入维度不够能量泄漏到相邻奇异值。解决把 L 提到 547 或 730让半周年也有足够的嵌入维度或者先对序列做 2 天采样平均把 Nyquist 降下来再取 L365。4.4 现象重构残差里出现明显 2~3 天周期原因这是 GNSS 日解序列里常见的轨道重复周期残留小波分解时落在 d1 或 d2如果这两层没丢弃而是送进 SSASSA 会把它当成一个“振荡”重构出来叠回去后污染坐标序列。解决d1、d2 直接置零或者在做 SSA 之前对这两层做低通滤波截止频率设在 0.1 cycle/day 以下。4.5 现象同一站点不同坐标分量N、E、U用同一套参数U 方向效果明显差原因U 方向噪声水平通常是水平分量的 2~3 倍且多路径影响更重同样的分解层数和 SSA 窗口长度在 U 方向会把噪声当成信号。解决U 方向单独调参小波分解层数加 1SSA 窗口长度取水平方向的 1.5 倍分组时只保留奇异值明显成对且振幅超过噪声 RMS 3 倍以上的分量。5. 进阶把重构分量拿去做速度场和形变解释时的验证技巧走到这一步你手里已经有一条去噪后的坐标序列和几个分离出来的分量。但“去噪好看”不等于“速度估计更准”。我一般会做两个验证。第一个是分年速度一致性把去噪前后的序列分别按年做线性拟合看去噪后各年速度的离散度是否下降。如果去噪后离散度反而变大说明重构时把真实的构造信号当噪声削掉了需要回退分组把更多奇异值纳入重构。第二个是残差白噪声检验对去噪后的残差做 Ljung-Box 检验如果残差在滞后 30 天以内还有显著自相关说明 SSA 分组漏掉了某个周期分量通常是半年或季节项需要回到奇异值谱重新看成对结构。一个具体技巧是不要一次性把所有中间层都做 SSA。先对 d4 做看重构后的周年振幅是否和已知的区域水文负载或热膨胀模型对得上对得上再把 d5、d6 加进来。对不上说明这一层的周年不是物理信号可能是多路径的年变化强行重构只会让速度场有偏。我自己的习惯是保留一份“未做 SSA”的中间层分量和 SSA 重构后的分量并排画如果两者在非周年频段上的差异超过 1 毫米就检查是不是窗口长度 L 把某个低频构造信号也当成周期分离了。这套流程没有后悔药参数调错一步后面速度场解释就全歪所以每一步都留一份中间结果比事后反推省事得多。希望帮到你。本文还有配套的精品资源点击获取