简介这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的科研人员和测绘专业学生聚焦GPS站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题提出小波多尺度分解与奇异谱分析SSA相结合的建模思路。文档系统梳理了小波多分辨率分析的低频概貌与高频细节分离原理、SSA构造时滞矩阵与奇异值分解的四步流程并说明两者结合可优势互补从噪声中提取周期性、长期性变化信息。文中以全球11个测站20年1999—2018年GPS垂向坐标时间序列为实例验证该方法对提高坐标精度的有效性同时讨论了大气负荷、水文负载等非线性因素及ITRF框架精度问题。资源包内含1个docx文档约442KB结构完整、公式与图示清晰适合作为相关方向建模与实验的参考材料。目前已有186人学习。1. 小波多尺度分解与奇异谱分析GNSS 坐标时间序列里到底谁在“抖”GNSS 站坐标时间序列看起来就是三条曲线——东向、北向、垂向随时间变化但真正做过数据处理的人都知道这三条曲线里混着太多东西构造运动趋势、周年和半周年周期、潮汐混叠、多路径残余、天线相位中心变化、甚至换接收机那天的一个阶跃。直接拿原始序列去做速度场估计噪声一大线性项就被淹了。小波多尺度分解和奇异谱分析SSA就是用来把这条曲线拆开的两种手段小波按频率把信号分层SSA 按时间上的相关结构把振荡和趋势分离。两者结合能让你在 GNSS 坐标时间序列里把周期项、趋势项和噪声项分得比单纯滤波干净得多。这套方法适合做地壳形变分析、站速度场解算、参考站稳定性监测的从业者尤其是那些已经拿到 RINEX 观测文件、跑完基线解算、手里有一堆 .pos 或 .neu 坐标序列的人。下面按“先立住原理、再动手复现、最后说坑”的顺序讲清楚。2. 小波多尺度分解在 GNSS 坐标序列上的选型与实现2.1 为什么 GNSS 坐标序列适合小波而不是简单滑动平均GNSS 坐标时间序列的采样率常见有 30 秒、5 分钟、天解几种。日解序列里周年和半周年周期是主要周期项周期分别约 365 天和 182 天高频部分则混着多路径、观测噪声和短周期潮汐残余。滑动平均只有一个尺度窗口一长就把周年项抹平窗口一短又压不住高频噪声。小波分解的本质是把序列投影到一组由母小波伸缩平移得到的基函数上得到不同尺度上的系数每个尺度对应一个频带。对 GNSS 日解序列常用 db4 或 sym4 小波做 5 到 8 层分解层数由序列长度和你要保留的最低频决定。选小波时看两个指标消失矩阶数和支撑长度。消失矩越高对多项式趋势的压制越好支撑越短边界效应越小。GNSS 序列里趋势项接近线性db4 的 4 阶消失矩足够sym4 的对称性更好重构时相位失真小。我一般先用 sym4 试边界用对称延拓避免端点附近出现虚假振荡。2.2 用 PyWavelets 做多尺度分解的最小可跑代码下面这段代码读取一个三列文本坐标序列时间、东向、北向、垂向对垂向做 6 层 sym4 分解输出各层系数并重构趋势项和周期项。数据格式假设是空白分隔第一列是简化儒略日或年积日后面三列是坐标值单位米。import numpy as np import pywt import matplotlib.pyplot as plt # 读取坐标序列假设文件为 neu.txt四列t, east, north, up data np.loadtxt(neu.txt) t data[:, 0] up data[:, 3] - np.mean(data[:, 3]) # 去均值避免重构时直流分量干扰 # 小波分解参数 wavelet sym4 level 6 mode symmetric # 边界延拓方式减少端点效应 # 多尺度分解 coeffs pywt.wavedec(up, wavelet, modemode, levellevel) # coeffs[0] 是近似系数最低频coeffs[1:] 是各层细节系数从高频到低频 # 重构趋势项只保留近似系数细节置零 coeffs_trend [coeffs[0]] [np.zeros_like(c) for c in coeffs[1:]] trend pywt.waverec(coeffs_trend, wavelet, modemode) trend trend[:len(up)] # 重构长度可能多一两个点截断 # 重构高频噪声只保留第一层细节 coeffs_noise [np.zeros_like(coeffs[0])] [coeffs[1]] \ [np.zeros_like(c) for c in coeffs[2:]] noise pywt.waverec(coeffs_noise, wavelet, modemode) noise noise[:len(up)] # 周期项 原始 - 趋势 - 高频噪声 periodic up - trend - noise # 输出各层能量占比帮助判断分解层数是否合理 energy [np.sum(c**2) for c in coeffs] energy_ratio [e / sum(energy) for e in energy] print(各层能量占比近似层在前:, np.round(energy_ratio, 4)) # 简单绘图检查 plt.figure(figsize(10, 6)) plt.plot(t, up, label原始, alpha0.5) plt.plot(t, trend, label趋势, linewidth2) plt.plot(t, periodic, label周期项, linewidth1) plt.legend() plt.savefig(decomp.png, dpi150)逻辑说明pywt.wavedec返回的系数列表长度是 level1第一个是近似系数后面依次是第 level 层到第 1 层的细节系数。重构时把不需要的层置零再waverec就得到对应频带的成分。参数方面level的选择有个经验公式level floor(log2(N/(L-1)))N 是序列长度L 是小波滤波器长度。sym4 的 L 是 8如果序列有 2000 个历元log2(2000/7) 约 8.1取 6 到 7 层都合理。mode选 symmetric 是因为 GNSS 序列两端没有周期性零填充会在端点产生突变对称延拓更平滑。能量占比表用来判断如果近似层能量占比超过 95%说明分解层数太少周期项还没被分出来如果第一层细节能量占比超过 30%说明高频噪声很强可能需要先做粗差剔除。2.3 分解层数和边界效应的参数怎么定层数不是越多越好。层数过多低频近似系数点数太少重构趋势时端点会翘。我一般按“最低频成分至少保留两个完整周期”来定如果你关心周年项最低频层的中心周期要大于 2 年对应尺度约 2*365 天。对日解序列6 层 sym4 分解后最低频层对应的周期大约在 1 到 2 年量级刚好把周年和半周年分开。边界效应方面symmetric模式在端点处镜像延拓重构后前 L/2 个点和后 L/2 个点仍可能有偏差处理时可以把两端各 15 个历元标记为不可用或者用pywt.wavedec的modeperiodization做周期延拓但后者要求序列本身近似周期GNSS 趋势项会让它产生阶跃。实际项目中我倾向 symmetric 加两端剔除简单可靠。3. 奇异谱分析把周期项从噪声里“抠”出来的嵌入与重构3.1 SSA 的轨迹矩阵、嵌入维数和重构阶数SSA 的核心是把一维序列嵌入到高维空间构造轨迹矩阵再做奇异值分解最后按奇异值大小分组重构。对 GNSS 坐标序列SSA 擅长分离出振幅和相位都随时间变化的振荡分量比如周年项的年际调制。关键参数有两个窗口长度 L也叫嵌入维数和重构阶数 r。L 一般取序列长度的三分之一到二分之一但不要超过 N/2。对 2000 个历元的日解序列L 取 365 到 600 都常见。L 越大频率分辨率越高但计算量按 L^2 增长。重构阶数 r 决定保留多少个奇异值对应的成分。如果只关心趋势和周年项r 取 2 到 4 就够第一对奇异值通常对应趋势第二对对应周年第三对对应半周年。判断依据是奇异值谱的“拐点”——奇异值从大到小排列前几个下降很快后面趋于平缓拐点前的就是信号主导成分。3.2 用 Python 实现 SSA 并重构趋势与周年项下面代码对垂向序列做 SSA窗口 L365重构前 4 个成分输出趋势和周年项。注意 SSA 重构需要对角平均把轨迹矩阵还原成一维序列。import numpy as np def ssa_decompose(series, L, r): series: 一维序列 L: 窗口长度 r: 重构阶数保留前 r 个奇异值 返回: 重构后的序列趋势周期以及残差 N len(series) K N - L 1 # 构造轨迹矩阵 X np.zeros((L, K)) for i in range(K): X[:, i] series[i:iL] # SVD U, s, Vt np.linalg.svd(X, full_matricesFalse) # 用前 r 个成分重构轨迹矩阵 X_rec np.zeros_like(X) for i in range(r): X_rec s[i] * np.outer(U[:, i], Vt[i, :]) # 对角平均还原一维序列 rec np.zeros(N) count np.zeros(N) for i in range(L): for j in range(K): rec[ij] X_rec[i, j] count[ij] 1 rec rec / count residual series - rec return rec, residual, s # 使用示例 data np.loadtxt(neu.txt) up data[:, 3] - np.mean(data[:, 3]) L 365 r 4 rec, residual, s ssa_decompose(up, L, r) print(前 10 个奇异值:, np.round(s[:10], 4)) print(前 4 个成分方差贡献:, np.round(np.sum(s[:4]**2)/np.sum(s**2), 4)) # 保存结果 np.savetxt(ssa_rec.txt, np.column_stack([data[:, 0], rec, residual]))逻辑说明轨迹矩阵构造时L 决定嵌入维数KN-L1 是列数。SVD 之后奇异值从大到小排列前 r 个对应主要成分。对角平均是 SSA 重构的关键步骤把轨迹矩阵中同一时间索引的元素取平均还原出一维序列。参数 L 的选择影响频率分辨率和计算时间L365 对日解序列意味着窗口覆盖一年能分辨周年和半周年L 太小则低频成分分不开L 太大则 SVD 耗时明显增加。r 的选择看奇异值谱如果前 4 个奇异值之和占总方差 80% 以上取 r4 合理如果前 2 个就占 90%说明趋势和周年主导取 r2 也行。残差项里通常还混着未建模的高频噪声和短周期项可以再对小波分解的高频层做一次 SSA或者直接当噪声处理。3.3 小波和 SSA 串起来用的顺序问题先小波还是先 SSA结果不一样。我一般先小波分解把最高频的两层细节去掉相当于预滤波再做 SSA。原因是 SSA 对高频噪声敏感噪声太强时奇异值谱的拐点不明显重构阶数不好定。小波预滤波后序列信噪比提高SSA 的奇异值谱更干净。反过来先 SSA 再小波也可以但 SSA 会把噪声分散到多个成分里小波再分解时层数不好控制。串行顺序没有绝对对错但先小波后 SSA 在 GNSS 日解序列上更稳。如果序列里有明显的阶跃换天线、换接收机无论哪种顺序都要先做阶跃探测和改正否则小波和 SSA 都会把阶跃当成低频成分污染趋势项。4. 避坑与排查GNSS 坐标序列做小波 SSA 时最容易翻车的五件事4.1 现象重构趋势项两端翘起中间正常原因小波分解的边界延拓方式与序列端点不匹配。GNSS 序列两端没有周期性零填充或周期延拓都会在端点产生突变重构时低频系数受污染。解决改用symmetric模式并把重构结果两端各 15 到 30 个历元剔除如果序列长度允许先做端点镜像延拓再分解分解完再截断。4.2 现象SSA 奇异值谱没有明显拐点重构阶数定不下来原因序列里高频噪声太强或者存在未改正的阶跃和粗差。噪声和阶跃会让奇异值缓慢下降拐点模糊。解决先做粗差探测3σ 或 IQR 准则再对小波高频层做软阈值去噪然后重新做 SSA。如果阶跃明显用ruptures库或手动分段先改正阶跃再分析。4.3 现象周年项振幅被低估周期项重构后比原始小一圈原因小波分解层数不够周年项被分到了近似层和细节层的交界处重构时被部分置零。或者 SSA 的窗口 L 太小频率分辨率不足以把周年和半周年分开。解决增加小波分解层数让周年项落在中间层SSA 的 L 至少取 365最好取 500 以上。检查方法对重构的周期项做 FFT看周年和半周年峰值是否与原始序列一致。4.4 现象换接收机那天之后趋势项整体平移原因坐标序列存在阶跃小波和 SSA 都把阶跃当成低频趋势的一部分重构后趋势项出现台阶。解决在做任何分解之前先做阶跃探测。常用方法有基于中位数的滑动窗口检验、贝叶斯阶跃模型。检测到阶跃后在序列中减去阶跃量再分解。阶跃改正不彻底后面所有分析都是白做。4.5 现象分解后残差里还有明显周期但振幅很小原因GNSS 坐标序列里除了周年和半周年还有周日、半日潮汐混叠产生的短周期项以及多路径的重复周期。小波层数少时这些短周期项被归入高频细节直接当噪声去掉了。解决如果分析目标包含短周期增加小波分解层数把对应频带单独重构或者对残差再做一次 SSA专门提取短周期成分。不要默认残差就是白噪声先看它的功率谱。5. 进阶技巧用交叉验证定参数别靠肉眼调小波层数和 SSA 的 L、r 这三个参数靠肉眼调容易过拟合。我现在的习惯是留出一段独立验证序列把前 80% 的数据用来分解和重构后 20% 不参与参数选择只用来算重构误差。具体做法是对每一组候选参数小波层数 4 到 8SSA 的 L 取 200、365、500r 取 2 到 6在前 80% 上做分解用重构的趋势和周期项外推后 20%算 RMSE。RMSE 最小的那组参数就是当前序列的最优参数。这个方法比看奇异值拐点更客观尤其适合多站批量处理。另一个技巧是分方向处理。东向、北向、垂向三个分量的噪声特性不一样垂向通常噪声最大周期项振幅也最大水平方向噪声小但构造信号可能更强。不要三个方向用同一套参数。我一般对垂向用更大的 L 和更多的分解层数水平方向适当减少。下面这张表是我在某区域连续运行站上总结的经验参数范围供参考。分量小波层数sym4SSA 窗口 L重构阶数 r备注东向5 到 6300 到 4002 到 4水平噪声小层数不宜过多北向5 到 6300 到 4002 到 4同东向注意构造信号垂向6 到 7400 到 6003 到 5噪声大周期振幅大需更多层最后说一个验证重构是否合理的硬指标把重构的趋势项做线性拟合得到的斜率应该和用原始序列做抗差最小二乘估计的速度在 1 mm/yr 以内一致。如果差太多说明分解把趋势项也当成周期项去掉了或者阶跃没改干净。这个检查我每次批量处理完都会跑一遍比看图的后悔药管用。希望帮到你。本文还有配套的精品资源点击获取