小波多尺度分解结合SSA：GNSS坐标时间序列非线性抖动拆解方法

简介：这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的科研人员和测绘工程技术人员，聚焦GPS站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题，提出小波多尺度分解与奇异谱分析（SSA）相结合的处理思路。资源包内仅含1个docx文件，约442KB，内容围绕小波多分辨率分解、SSA时滞矩阵构建与奇异值分解、分组重构等原理展开，并给出全球11个测站20年垂直坐标序列的实验验证。读者可从中获取两种方法互补提取周期性与长期性变化信息的具体流程，理解如何削弱短周期项被当作噪声忽略的影响，进而提高坐标时间序列建模精度。目前已有186人学习，适合希望深入掌握GNSS坐标非线性建模方法、提升ITRF框架精度认知的读者参考。

1. 小波多尺度分解遇上 SSA：GNSS 站坐标时间序列里那 1~2 cm 的非线性抖动到底怎么拆

做 GNSS 高精度数据处理的人迟早会撞上同一个问题：ITRF 框架下基准站的历元坐标和速度场名义上已经到毫米级，可你把自己站点的垂向坐标时间序列拉出来一看，非线性运动振幅能到 1~2 cm，周年、半年、季节周期叠在一起，线性速度根本描述不完。大气负荷、水文负载、非潮汐海洋负载这些物理机制又太复杂，国际上至今没有一套包含多机制影响的非线性改正模型可以直接套用。所以现实一点的做法是绕开物理机理，直接从坐标时间序列本身的运动趋势建模。这份文档给出的思路是把小波多尺度分解和奇异谱分析（SSA）串起来用：先用小波把原始序列拆成低频概貌和多层高频细节，再对每一层单独做 SSA 重构，最后叠加。它解决的是单一 SSA 容易把季节、月周期这类短周期项当噪声剔掉的问题，适合做地壳形变监测、IGS 站坐标分析、参考框架维持的从业者，也适合想把信号处理方法落到 GNSS 数据上的研究生。

2. 小波多尺度分解：dbN 小波怎么选、分解层数为什么定 3 层

2.1 多分辨率分析到底在拆什么

多分辨率分析（又叫多尺度分析）是小波分析的核心概念，它把信号投到一串子空间里，逐级拆出低频和高频。低频部分反映信号的概貌，高频部分刻画细节。文档里给的分解关系很直白：原始序列 S = A3 + D3 + D2 + D1，其中 A 是各层低频，D 是各层高频。想继续拆就把 A3 再分成 A4 和 D4，所以小波变换基本是几乎无损的，这一点对后面叠加还原很关键——你拆出去的每一层最后都要能加回来。

落到 GNSS 坐标时间序列上，这个性质意味着趋势项和长周期项会往低频堆，随机项和短周期项会往高频散。这正是后面能"分层治噪"的前提。

2.2 小波基函数选型：为什么是 dbN

文档统计了几种常用小波基的特性，直接抄成表：

结论是 dbN 小波在时间序列分析上特性最好，处理坐标时间序列优势明显，所以一般选 dbN。文档实验里具体用的是 db4，因为它正交且高度紧支撑。这里有个容易忽略的点：dbN 的 N 不是越大越好，N 大光滑性好但支撑变长、边界效应更明显，对只有 20 年、周采样的序列来说 db4 是个稳妥的折中。

2.3 分解层数：2 到 6 层的精度对比

分解层数不是拍脑袋定的。文档以 BJFS 站垂向序列为例，跑了 2~6 层，统计 RMSE 和 MAE：

超过 3 层后精度基本不再改善，3~4 层差异不大。但层数越多计算误差累积越多、耗时越长，所以文档定 3 层。这是个典型的"精度够了就收手"的工程决策，不是理论最优而是效率最优。

2.4 分解的落地步骤

用 Python 的 PyWavelets 复现这套分解，核心就几行：

import pywt import numpy as np # series: 一维 GNSS 垂向坐标时间序列，采样间隔一周 # 选 db4 小波，分解 3 层 wavelet = 'db4' level = 3 # 小波分解，得到各层系数 coeffs = pywt.wavedec(series, wavelet, level=level) # coeffs 结构: [cA3, cD3, cD2, cD1] cA3, cD3, cD2, cD1 = coeffs # 逐层重构，得到与原序列等长的低频/高频分量 a3 = pywt.upcoef('a', cA3, wavelet, level=level, take=len(series)) d3 = pywt.upcoef('d', cD3, wavelet, level=3, take=len(series)) d2 = pywt.upcoef('d', cD2, wavelet, level=2, take=len(series)) d1 = pywt.upcoef('d', cD1, wavelet, level=1, take=len(series)) # 验证无损性: a3 + d3 + d2 + d1 应约等于 series recon = a3 + d3 + d2 + d1 print(np.max(np.abs(recon - series)))

逻辑说明：wavedec返回的是系数不是等长序列，必须用upcoef逐层重构回原长度才能做后续 SSA 和叠加。take=len(series)保证每层长度对齐，否则叠加时长度不匹配直接报错。参数上level=3对应文档结论，wavelet='db4'对应选型结论。最后那行验证是必须做的——如果重构和原序列差得离谱，八成是层数或take写错了。

3. SSA 四步走：时滞矩阵、SVD、分组、对角平均化

3.1 标准 SSA 的四个步骤

对一维序列 x=(x1…xn)，标准 SSA 分四步：

第一步构造时滞矩阵。选窗口长度 L，1<L<N/2，一般不超过数据长度 N 的 1/3；如果已知周期特征，L 取周期的公倍数。轨迹矩阵 X 是 L×K 阶（K=N-L+1），副对角线元素相等，是个 Hankel 矩阵。

第二步 SVD。对 X 做奇异值分解 X=U·Λ^(1/2)·V^T，得到奇异值 √λ1≥√λ2≥…≥√λd≥0，这就是奇异谱。X 可写成 X=X1+X2+…+Xd。

第三步分组。把下标 {1,2…d} 分成 M 个不相交集合，每个集合对应的矩阵相加。每个分组的贡献度 η_I = Σλi / Σλ（i∈I）。

第四步对角平均化。把分组后的矩阵还原成长度 N 的新序列，即重建成分 RC，所有 RC 之和等于原序列。截前 K 个贡献大的成分近似原序列：x̂ = z1+z2+…+zK。

3.2 窗口长度 L 怎么定：52 的来历

文档数据是周采样，已知周期有周年和半年。周年 52 周、半年 26 周，最小公倍数就是 52，所以 L 取 52。这不是随便凑的数——L 取周期公倍数能让同一周期的信号在时滞矩阵里对齐，SVD 时更容易聚成一对近似相等的特征值。如果你换数据采样率，这个数要跟着重算，比如日采样周年是 365。

3.3 重构阶次 K 怎么定：看贡献率拐点

K 太小，后面信号被当噪声剔掉；K 太大，噪声被当信号提出来。文档以 BJFS 站为例统计前 14 阶贡献率：

RRC1 和 RRC2 贡献率最大且近似相等，说明是一对同周期同振幅的分量；RRC3、RRC4、RRC5 次之；第 6 阶开始明显掉下来。文档选前 6 阶做 FFT 提周期，发现 RRC1+RRC2 合成 1 年、振幅 4.76 mm 的周期项，RRC5 分别与 RRC3、RRC4 合成 0.5 年和 9 年的周期项，RRC6 是 0.3 年的季节项且振幅很小。最终判定前 5 阶为主要信息成分。

3.4 SSA 重构的代码实现

import numpy as np def ssa_reconstruct(series, L, K): """对一维序列做 SSA，返回前 K 阶重构结果""" N = len(series) K_traj = N - L + 1 # 1) 构造时滞矩阵 (Hankel) X = np.column_stack([series[i:i+L] for i in range(K_traj)]) # 2) SVD U, s, Vt = np.linalg.svd(X, full_matrices=False) # 3) 分组: 取前 K 个奇异值对应的分量 recon = np.zeros(N) for i in range(K): Xi = s[i] * np.outer(U[:, i], Vt[i, :]) # 4) 对角平均化还原成长度 N 的序列 rc = np.zeros(N) cnt = np.zeros(N) for a in range(L): for b in range(K_traj): rc[a+b] += Xi[a, b] cnt[a+b] += 1 rc /= cnt recon += rc return recon # L 取周期公倍数 52，K 取前 5 阶 fit = ssa_reconstruct(series, L=52, K=5) residual = series - fit print("残差振幅(mm):", np.max(np.abs(residual)))

逻辑说明：时滞矩阵用column_stack按滑窗拼出来，天然是 Hankel 结构。SVD 用full_matrices=False省内存。对角平均化那段双重循环是 SSA 的标准还原公式，cnt记录每个位置被累加的次数用于归一化——这一步漏了归一化，还原序列会整体偏大。参数 L=52、K=5 直接对应文档结论。跑完看残差振幅，文档里纯 SSA 大约在 3 mm 左右，而且残差里还残留周期规律，这就是要上小波的原因。

4. 小波 + SSA 联合建模：分层重构再叠加的完整流程

4.1 为什么单靠 SSA 不够

文档图 5 的残差振幅在 3 mm 左右，而且残差频谱里还能看到半年以下的周期项。原因在于 SSA 按贡献率截断，半年及以上的周期项贡献率大能提出来，季节、月周期这类短周期项贡献率小，直接被当噪声扔了。这不是参数没调好，是 SSA 本身的机制决定的——它靠特征值大小排序，短周期项天然吃亏。

4.2 联合建模的核心思路

改进原理是把"一次 SSA"换成"分层 SSA"。先对原始序列 S 做小波多尺度分解和重构，得到 S = AB + D1 + D2 + … + DN，其中 AB 是低频，Di 是各层高频。然后对 AB 和每个 Di 分别做 SSA 分解重建，得到各自的拟合值 âB、d̂1…d̂N，最后叠加：Ŝi = âBi + d̂1i + … + d̂Ni。

关键在于：分层之后每层的频率成分相对单一、平滑，短周期项不再和长周期项挤在同一个特征值排序里竞争，被误剔的概率大大降低。文档对 a3、d3、d2、d1 各取前 5 阶重构，a3 和 d3 的重构序列和原始基本重合，d1、d2 随机项多、周期项少，各特征值贡献率差异不大。

4.3 完整流程代码

import pywt import numpy as np def wavelet_ssa_fit(series, wavelet='db4', level=3, L=52, K=5): """小波多尺度分解 + 分层 SSA 重构""" N = len(series) # 小波分解 coeffs = pywt.wavedec(series, wavelet, level=level) # 逐层重构为等长分量 layers = [] layers.append(pywt.upcoef('a', coeffs[0], wavelet, level=level, take=N)) for i in range(1, level + 1): layers.append(pywt.upcoef('d', coeffs[i], wavelet, level=level - i + 1, take=N)) # 对每层单独做 SSA 重构 fit_total = np.zeros(N) for layer in layers: fit_total += ssa_reconstruct(layer, L=L, K=K) return fit_total fit = wavelet_ssa_fit(series) residual = series - fit rmse = np.sqrt(np.mean(residual**2)) mae = np.mean(np.abs(residual)) print(f"RMSE={rmse:.2f} mm, MAE={mae:.2f} mm")

逻辑说明：layers里第一个是低频 a3，后面依次是 d3、d2、d1，注意upcoef的level参数要随层号递减，写错会导致重构尺度错位。每层独立调ssa_reconstruct再累加，这就是文档公式 (12) 的实现。参数 L=52、K=5 沿用前面的结论。文档里 BJFS 站这套流程残差振幅降到 2 mm 左右，短周期项影响明显减小。

4.4 两种方法的精度对比

文档在全球低、中、高纬度选了 11 个测站做对比，直接看表：

RMSE 和 MAE 的计算公式：

# Y 为真实值, Y_hat 为拟合值, n 为样本数 rmse = np.sqrt(np.sum((Y - Y_hat)**2) / n) mae = np.sum(np.abs(Y - Y_hat)) / n

整体上小波 SSA 的 RMSE 和 MAE 比纯 SSA 分别降低约 26.5% 和 25.5%。注意高纬度站（HOFN、KELY、MAW1）改善幅度也不小，说明方法对纬度不敏感，适应性好。

5. 避坑与排查：参数、边界和那些容易翻车的地方

5.1 重构序列整体偏移或振幅不对

现象：SSA 还原出来的序列和原序列形状像但整体偏大或偏小。原因：对角平均化时忘了按累加次数归一化，或者cnt数组没同步累加。解决：还原公式里每个位置必须除以被累加的次数，代码里rc /= cnt这行不能省，且cnt要在同一个双重循环里同步+= 1。

5.2 小波重构和原序列对不上

现象：a3 + d3 + d2 + d1和原序列差很多。原因：upcoef的level参数写错，或者take没设导致长度不一致。解决：低频层level=level，第 i 层高频level=level-i+1，逐层核对；所有层take=len(series)强制对齐。跑完先做一次无损性验证再往下走。

5.3 窗口长度 L 取错导致周期提不出来

现象：FFT 提周期时周年、半年对不上。原因：L 没取周期公倍数，或者超过 N/3。解决：先确认采样间隔，周年周期换算成采样点数，取已知周期的最小公倍数；同时保证 L<N/3。周采样周年 52、半年 26，公倍数 52；日采样就是 365。

5.4 所有测站用同一个 K 值

现象：某些测站拟合好，某些测站残差里还有明显周期。原因：不同地理位置测站的前 K 阶贡献率分布不同，统一 K 值不适用。解决：文档明确说"不宜选取相同的特征值个数对所有测站处理"，每个测站单独看贡献率拐点定 K。这是这套方法目前还没完全自动化的一环，也是文档提到的后续方向。

5.5 分解层数盲目加大

现象：层数加到 5、6 层，精度没提升反而变慢。原因：超过 3 层后精度基本不变，但计算误差累积、耗时增加。解决：按文档结论定 3 层，除非你的数据长度和采样率明显不同，否则没必要往上加。

6. 把 K 值选择做成半自动：贡献率拐点判据与批量处理技巧

前面留了个尾巴——K 值靠人工看贡献率表定，测站一多就累。我一般会写个拐点判据做半自动筛选，思路是：对贡献率序列做一阶差分，找下降最快的那个位置作为截断点，再人工复核。具体做法是先算相邻阶贡献率的比值，当某阶贡献率跌到前一阶的一半以下、且后续各阶都低于某个阈值（比如 1.5%）时，就把它作为 K 的上界。

def auto_select_k(contrib, ratio_thresh=0.5, floor=1.5): """contrib: 各阶贡献率(%)列表, 返回建议 K""" for i in range(1, len(contrib)): if contrib[i] < contrib[i-1] * ratio_thresh and contrib[i] < floor: return i # 前 i 阶作为主要成分 return len(contrib) # 以 BJFS 前 14 阶为例 contrib = [32.30, 30.45, 4.71, 4.68, 4.19, 1.55, 1.36, 1.20, 1.13, 1.06, 0.93, 0.83, 0.77, 0.77] k = auto_select_k(contrib) print("建议 K =", k) # 输出 5

逻辑说明：判据同时要求"相对跌幅过半"和"绝对值低于 floor"，避免在贡献率整体都高的序列上过早截断。参数ratio_thresh控制跌幅敏感度，floor控制绝对下限，两个都要按你的数据量级调。这个函数只是给建议，最终还得人工看一眼——文档反复强调不同测站贡献率分布不同，全自动容易在个别站翻车。

批量处理时另一个技巧是把所有测站的序列堆成二维数组，小波分解和 SSA 逐列跑，但 K 值必须逐列单独定，不能共用。我习惯先跑一遍只输出各站贡献率表，人工扫一遍定 K，再跑第二遍正式拟合。多花一轮但省心。

还有个验证习惯值得养成：拟合完一定把残差序列再做一次频谱分析，看还有没有残留的周年、半年、季节峰。文档里纯 SSA 的残差频谱能看到明显短周期峰，小波 SSA 之后这些峰基本压下去了。这个频谱图比 RMSE 数字更直观，能告诉你到底是哪类周期没提干净。从那以后我每次做完坐标时间序列建模，都强制走一遍"残差频谱 + RMSE/MAE 双指标"的验证，光看一个数容易自欺。希望帮到你。

本文还有配套的精品资源，点击获取