尧图网站建设 
简介这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的科研人员和测绘工程技术人员聚焦GPS站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题提出小波多尺度分解与奇异谱分析SSA相结合的处理思路。资源包内仅含1个docx文件约442KB内容围绕小波多分辨率分解、SSA时滞矩阵构建与奇异值分解、分组重构等原理展开并给出全球11个测站20年垂直坐标序列的实验验证。读者可从中获取两种方法互补提取周期性与长期性变化信息的具体流程理解如何削弱短周期项被当作噪声忽略的影响进而提高坐标时间序列建模精度。目前已有186人学习适合希望深入掌握GNSS坐标非线性建模方法、提升ITRF框架精度认知的读者参考。1. 小波多尺度分解遇上 SSAGNSS 站坐标时间序列里那 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文档统计了几种常用小波基的特性直接抄成表小波基支撑长度消失矩阶数对称性特点Haar11对称时域不连续频率局部性差dbN2NN近似对称光滑性随 N 增加而增加性能较好symN2NN近似对称减少相位失真更适合图像处理Meyer有限长度—对称分析会产生失真bior重构 2Nr1 / 分解 2Nd1Nr-1不对称具线性相位重构中常用结论是 dbN 小波在时间序列分析上特性最好处理坐标时间序列优势明显所以一般选 dbN。文档实验里具体用的是 db4因为它正交且高度紧支撑。这里有个容易忽略的点dbN 的 N 不是越大越好N 大光滑性好但支撑变长、边界效应更明显对只有 20 年、周采样的序列来说 db4 是个稳妥的折中。2.3 分解层数2 到 6 层的精度对比分解层数不是拍脑袋定的。文档以 BJFS 站垂向序列为例跑了 2~6 层统计 RMSE 和 MAE分解层数23456RMSE/mm2.151.881.801.801.80MAE/mm1.681.491.431.441.44超过 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, levellevel) # coeffs 结构: [cA3, cD3, cD2, cD1] cA3, cD3, cD2, cD1 coeffs # 逐层重构得到与原序列等长的低频/高频分量 a3 pywt.upcoef(a, cA3, wavelet, levellevel, takelen(series)) d3 pywt.upcoef(d, cD3, wavelet, level3, takelen(series)) d2 pywt.upcoef(d, cD2, wavelet, level2, takelen(series)) d1 pywt.upcoef(d, cD1, wavelet, level1, takelen(series)) # 验证无损性: a3 d3 d2 d1 应约等于 series recon a3 d3 d2 d1 print(np.max(np.abs(recon - series)))逻辑说明wavedec返回的是系数不是等长序列必须用upcoef逐层重构回原长度才能做后续 SSA 和叠加。takelen(series)保证每层长度对齐否则叠加时长度不匹配直接报错。参数上level3对应文档结论waveletdb4对应选型结论。最后那行验证是必须做的——如果重构和原序列差得离谱八成是层数或take写错了。3. SSA 四步走时滞矩阵、SVD、分组、对角平均化3.1 标准 SSA 的四个步骤对一维序列 x(x1…xn)标准 SSA 分四步第一步构造时滞矩阵。选窗口长度 L1LN/2一般不超过数据长度 N 的 1/3如果已知周期特征L 取周期的公倍数。轨迹矩阵 X 是 L×K 阶KN-L1副对角线元素相等是个 Hankel 矩阵。第二步 SVD。对 X 做奇异值分解 XU·Λ^(1/2)·V^T得到奇异值 √λ1≥√λ2≥…≥√λd≥0这就是奇异谱。X 可写成 XX1X2…Xd。第三步分组。把下标 {1,2…d} 分成 M 个不相交集合每个集合对应的矩阵相加。每个分组的贡献度 η_I Σλi / Σλi∈I。第四步对角平均化。把分组后的矩阵还原成长度 N 的新序列即重建成分 RC所有 RC 之和等于原序列。截前 K 个贡献大的成分近似原序列x̂ z1z2…zK。3.2 窗口长度 L 怎么定52 的来历文档数据是周采样已知周期有周年和半年。周年 52 周、半年 26 周最小公倍数就是 52所以 L 取 52。这不是随便凑的数——L 取周期公倍数能让同一周期的信号在时滞矩阵里对齐SVD 时更容易聚成一对近似相等的特征值。如果你换数据采样率这个数要跟着重算比如日采样周年是 365。3.3 重构阶次 K 怎么定看贡献率拐点K 太小后面信号被当噪声剔掉K 太大噪声被当信号提出来。文档以 BJFS 站为例统计前 14 阶贡献率阶次贡献率/%阶次贡献率/%132.3081.20230.4591.1334.71101.0644.68110.9354.19120.8361.55130.7771.36140.77RRC1 和 RRC2 贡献率最大且近似相等说明是一对同周期同振幅的分量RRC3、RRC4、RRC5 次之第 6 阶开始明显掉下来。文档选前 6 阶做 FFT 提周期发现 RRC1RRC2 合成 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:iL] for i in range(K_traj)]) # 2) SVD U, s, Vt np.linalg.svd(X, full_matricesFalse) # 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[ab] Xi[a, b] cnt[ab] 1 rc / cnt recon rc return recon # L 取周期公倍数 52K 取前 5 阶 fit ssa_reconstruct(series, L52, K5) residual series - fit print(残差振幅(mm):, np.max(np.abs(residual)))逻辑说明时滞矩阵用column_stack按滑窗拼出来天然是 Hankel 结构。SVD 用full_matricesFalse省内存。对角平均化那段双重循环是 SSA 的标准还原公式cnt记录每个位置被累加的次数用于归一化——这一步漏了归一化还原序列会整体偏大。参数 L52、K5 直接对应文档结论。跑完看残差振幅文档里纯 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, waveletdb4, level3, L52, K5): 小波多尺度分解 分层 SSA 重构 N len(series) # 小波分解 coeffs pywt.wavedec(series, wavelet, levellevel) # 逐层重构为等长分量 layers [] layers.append(pywt.upcoef(a, coeffs[0], wavelet, levellevel, takeN)) for i in range(1, level 1): layers.append(pywt.upcoef(d, coeffs[i], wavelet, levellevel - i 1, takeN)) # 对每层单独做 SSA 重构 fit_total np.zeros(N) for layer in layers: fit_total ssa_reconstruct(layer, LL, KK) 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(fRMSE{rmse:.2f} mm, MAE{mae:.2f} mm)逻辑说明layers里第一个是低频 a3后面依次是 d3、d2、d1注意upcoef的level参数要随层号递减写错会导致重构尺度错位。每层独立调ssa_reconstruct再累加这就是文档公式 (12) 的实现。参数 L52、K5 沿用前面的结论。文档里 BJFS 站这套流程残差振幅降到 2 mm 左右短周期项影响明显减小。4.4 两种方法的精度对比文档在全球低、中、高纬度选了 11 个测站做对比直接看表测站纬度SSA RMSE/mmSSA MAE/mm小波SSA RMSE/mm小波SSA MAE/mmADIS9.2°N2.591.981.901.44TUVA8.3°S2.942.172.231.71NAUR0.3°S2.822.092.091.61IISC13.1°N2.732.001.981.47WUHN30.5°N2.902.312.111.66STR135.2°S2.501.811.831.38BJFS39.6°N2.481.931.881.49MAC154.3°S2.121.651.651.27HOFN64.2°N2.131.671.661.30KELY66.6°N2.642.051.821.43MAW167.4°S2.191.681.461.14RMSE 和 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没设导致长度不一致。解决低频层levellevel第 i 层高频levellevel-i1逐层核对所有层takelen(series)强制对齐。跑完先做一次无损性验证再往下走。5.3 窗口长度 L 取错导致周期提不出来现象FFT 提周期时周年、半年对不上。原因L 没取周期公倍数或者超过 N/3。解决先确认采样间隔周年周期换算成采样点数取已知周期的最小公倍数同时保证 LN/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_thresh0.5, floor1.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 双指标的验证光看一个数容易自欺。希望帮到你。本文还有配套的精品资源点击获取 小波多尺度分解结合SSA：GNSS坐标时间序列非线性抖动拆解方法
简介这份文档面向从事GNSS数据处理、地壳形变监测与地球动力学研究的科研人员和测绘工程技术人员聚焦GPS站坐标时间序列中非线性运动趋势难以用线性速度完整描述的问题提出小波多尺度分解与奇异谱分析SSA相结合的处理思路。资源包内仅含1个docx文件约442KB内容围绕小波多分辨率分解、SSA时滞矩阵构建与奇异值分解、分组重构等原理展开并给出全球11个测站20年垂直坐标序列的实验验证。读者可从中获取两种方法互补提取周期性与长期性变化信息的具体流程理解如何削弱短周期项被当作噪声忽略的影响进而提高坐标时间序列建模精度。目前已有186人学习适合希望深入掌握GNSS坐标非线性建模方法、提升ITRF框架精度认知的读者参考。1. 小波多尺度分解遇上 SSAGNSS 站坐标时间序列里那 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文档统计了几种常用小波基的特性直接抄成表小波基支撑长度消失矩阶数对称性特点Haar11对称时域不连续频率局部性差dbN2NN近似对称光滑性随 N 增加而增加性能较好symN2NN近似对称减少相位失真更适合图像处理Meyer有限长度—对称分析会产生失真bior重构 2Nr1 / 分解 2Nd1Nr-1不对称具线性相位重构中常用结论是 dbN 小波在时间序列分析上特性最好处理坐标时间序列优势明显所以一般选 dbN。文档实验里具体用的是 db4因为它正交且高度紧支撑。这里有个容易忽略的点dbN 的 N 不是越大越好N 大光滑性好但支撑变长、边界效应更明显对只有 20 年、周采样的序列来说 db4 是个稳妥的折中。2.3 分解层数2 到 6 层的精度对比分解层数不是拍脑袋定的。文档以 BJFS 站垂向序列为例跑了 2~6 层统计 RMSE 和 MAE分解层数23456RMSE/mm2.151.881.801.801.80MAE/mm1.681.491.431.441.44超过 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, levellevel) # coeffs 结构: [cA3, cD3, cD2, cD1] cA3, cD3, cD2, cD1 coeffs # 逐层重构得到与原序列等长的低频/高频分量 a3 pywt.upcoef(a, cA3, wavelet, levellevel, takelen(series)) d3 pywt.upcoef(d, cD3, wavelet, level3, takelen(series)) d2 pywt.upcoef(d, cD2, wavelet, level2, takelen(series)) d1 pywt.upcoef(d, cD1, wavelet, level1, takelen(series)) # 验证无损性: a3 d3 d2 d1 应约等于 series recon a3 d3 d2 d1 print(np.max(np.abs(recon - series)))逻辑说明wavedec返回的是系数不是等长序列必须用upcoef逐层重构回原长度才能做后续 SSA 和叠加。takelen(series)保证每层长度对齐否则叠加时长度不匹配直接报错。参数上level3对应文档结论waveletdb4对应选型结论。最后那行验证是必须做的——如果重构和原序列差得离谱八成是层数或take写错了。3. SSA 四步走时滞矩阵、SVD、分组、对角平均化3.1 标准 SSA 的四个步骤对一维序列 x(x1…xn)标准 SSA 分四步第一步构造时滞矩阵。选窗口长度 L1LN/2一般不超过数据长度 N 的 1/3如果已知周期特征L 取周期的公倍数。轨迹矩阵 X 是 L×K 阶KN-L1副对角线元素相等是个 Hankel 矩阵。第二步 SVD。对 X 做奇异值分解 XU·Λ^(1/2)·V^T得到奇异值 √λ1≥√λ2≥…≥√λd≥0这就是奇异谱。X 可写成 XX1X2…Xd。第三步分组。把下标 {1,2…d} 分成 M 个不相交集合每个集合对应的矩阵相加。每个分组的贡献度 η_I Σλi / Σλi∈I。第四步对角平均化。把分组后的矩阵还原成长度 N 的新序列即重建成分 RC所有 RC 之和等于原序列。截前 K 个贡献大的成分近似原序列x̂ z1z2…zK。3.2 窗口长度 L 怎么定52 的来历文档数据是周采样已知周期有周年和半年。周年 52 周、半年 26 周最小公倍数就是 52所以 L 取 52。这不是随便凑的数——L 取周期公倍数能让同一周期的信号在时滞矩阵里对齐SVD 时更容易聚成一对近似相等的特征值。如果你换数据采样率这个数要跟着重算比如日采样周年是 365。3.3 重构阶次 K 怎么定看贡献率拐点K 太小后面信号被当噪声剔掉K 太大噪声被当信号提出来。文档以 BJFS 站为例统计前 14 阶贡献率阶次贡献率/%阶次贡献率/%132.3081.20230.4591.1334.71101.0644.68110.9354.19120.8361.55130.7771.36140.77RRC1 和 RRC2 贡献率最大且近似相等说明是一对同周期同振幅的分量RRC3、RRC4、RRC5 次之第 6 阶开始明显掉下来。文档选前 6 阶做 FFT 提周期发现 RRC1RRC2 合成 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:iL] for i in range(K_traj)]) # 2) SVD U, s, Vt np.linalg.svd(X, full_matricesFalse) # 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[ab] Xi[a, b] cnt[ab] 1 rc / cnt recon rc return recon # L 取周期公倍数 52K 取前 5 阶 fit ssa_reconstruct(series, L52, K5) residual series - fit print(残差振幅(mm):, np.max(np.abs(residual)))逻辑说明时滞矩阵用column_stack按滑窗拼出来天然是 Hankel 结构。SVD 用full_matricesFalse省内存。对角平均化那段双重循环是 SSA 的标准还原公式cnt记录每个位置被累加的次数用于归一化——这一步漏了归一化还原序列会整体偏大。参数 L52、K5 直接对应文档结论。跑完看残差振幅文档里纯 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, waveletdb4, level3, L52, K5): 小波多尺度分解 分层 SSA 重构 N len(series) # 小波分解 coeffs pywt.wavedec(series, wavelet, levellevel) # 逐层重构为等长分量 layers [] layers.append(pywt.upcoef(a, coeffs[0], wavelet, levellevel, takeN)) for i in range(1, level 1): layers.append(pywt.upcoef(d, coeffs[i], wavelet, levellevel - i 1, takeN)) # 对每层单独做 SSA 重构 fit_total np.zeros(N) for layer in layers: fit_total ssa_reconstruct(layer, LL, KK) 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(fRMSE{rmse:.2f} mm, MAE{mae:.2f} mm)逻辑说明layers里第一个是低频 a3后面依次是 d3、d2、d1注意upcoef的level参数要随层号递减写错会导致重构尺度错位。每层独立调ssa_reconstruct再累加这就是文档公式 (12) 的实现。参数 L52、K5 沿用前面的结论。文档里 BJFS 站这套流程残差振幅降到 2 mm 左右短周期项影响明显减小。4.4 两种方法的精度对比文档在全球低、中、高纬度选了 11 个测站做对比直接看表测站纬度SSA RMSE/mmSSA MAE/mm小波SSA RMSE/mm小波SSA MAE/mmADIS9.2°N2.591.981.901.44TUVA8.3°S2.942.172.231.71NAUR0.3°S2.822.092.091.61IISC13.1°N2.732.001.981.47WUHN30.5°N2.902.312.111.66STR135.2°S2.501.811.831.38BJFS39.6°N2.481.931.881.49MAC154.3°S2.121.651.651.27HOFN64.2°N2.131.671.661.30KELY66.6°N2.642.051.821.43MAW167.4°S2.191.681.461.14RMSE 和 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没设导致长度不一致。解决低频层levellevel第 i 层高频levellevel-i1逐层核对所有层takelen(series)强制对齐。跑完先做一次无损性验证再往下走。5.3 窗口长度 L 取错导致周期提不出来现象FFT 提周期时周年、半年对不上。原因L 没取周期公倍数或者超过 N/3。解决先确认采样间隔周年周期换算成采样点数取已知周期的最小公倍数同时保证 LN/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_thresh0.5, floor1.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 双指标的验证光看一个数容易自欺。希望帮到你。本文还有配套的精品资源点击获取