二阶时间重分配同步挤压变换：强瞬态信号时频分析利器

1. 内容整体设计与思路拆解第一次拿到 Draupner 波实测数据的时候我盯着那个高达 25.6 米的有效波高记录琢磨了很久。作为海洋工程领域最著名的“怪波”样本Draupner 波在 1995 年发生于北海挪威 Draupner 平台附近它之所以成为学术界研究焦点不只是因为波高极端更重要的是这个波在极短时间内从周围波场里突然“冒”出来包含极其剧烈的非线性调制和宽带瞬态信息。传统傅里叶分析只能告诉我们“有哪些频率分量”却完全无法回答“这些频率分量在什么时刻出现、强度如何变化”。而这恰恰是理解怪波生成机制、评估海洋结构极端载荷的关键。于是这几年我把大量精力花在了“非平稳信号的时频分析”上。什么谱图Spectrogram、小波尺度图、Wigner-Ville 分布基本都试过一轮。其中 Wigner-Ville 的时频聚焦性虽然好但交叉项干扰在 Draupner 这种频率成分密集的信号上特别头疼谱图和小波尺度图干净是干净能量却往往被抹得过于模糊两个在时间轴上距离很近的波峰分量经常粘连成一片。后来我接触到 Daubechies 等人提出的同步挤压变换Synchrosqueezing Transform, SST算是在时频聚焦性和可逆性之间找到了不错的平衡点。但一阶 SST 用在 Draupner 波上还是有个明显短板它的核心假设是信号在局部可以用“恒定频率-缓慢调制”近似一旦信号内部存在快速啁啾成分也就是瞬时频率随时间剧烈变化挤压结果就会出现显著的频率模糊。Draupner 波这种极端瞬态波恰恰是强啁啾的典型波峰附近瞬时频率变化速率极快频率再分配结果被拉出一道道模糊的“尾巴”。后来我读到了二阶时间重新分配同步挤压变换Second-order Time-Reassigned Synchrosqueezing Transform以下简称 TSST的思路。与频率方向的同步挤压不同时间重分配的思路是把分布在时间方向上的模糊能量重新“推回”到真实的信号事件时刻同时引入二阶局部调制估计来追踪瞬态信号的频率轨迹。这套方案可以说是专门冲着 Draupner 波这种强瞬态、强啁啾、宽带非平稳信号去的。我现在把完整实现过程和踩坑经验整理出来希望能给大家省下一些试错的成本。1.1 为什么传统方法在 Draupner 波上失效先用一个具体的例子说明问题。Draupner 波实测数据是以 1 Hz 左右采样间隔记录的波浪时序主能量分布在 0.05 Hz 附近但怪波事件本身会激发出一系列短时高频分量。如果只用谱图做分析窗口长度选得短了频率分辨率跟着崩掉窗口选得长了时间聚焦性又变得极差。这是海森堡不确定性原理决定的任何线性时频表示都逃不过。我在 2018 年刚拿到数据做复现实验时用 1024 点 Hann 窗的 STFT 谱图看整个 1200 秒的记录主峰从图上一眼能看出但把时间范围缩小到怪波事件前后各 50 秒时谱图上的能量分布简直是“一团雾”。这种模糊带来的后果很实际我想提取怪波发生时刻对应的瞬时频率曲线作为后续非线性强度评估的输入特征结果提取出来的曲线抖动剧烈根本无法稳定使用。现场工程人员拿着这种分析结果去评估什么只能得到一个很模糊的结论。所以问题的本质不在“选哪种窗口”而在于线性变换得到的结果本身就带有不确定性所决定的模糊性我们需要一种后处理手段把模糊的能量重新推回它真正该在的位置。1.2 同步挤压思路的核心逻辑同步挤压这个概念最早是从听觉场景分析里的“重分配方法”演化来的。传统重分配方法reassignment method会把时频平面中每个点的能量按局部相位信息重新映射到新的位置得到高分辨率的时频表示但代价是不再支持信号重建。同步挤压的巧妙之处在于它只在频率方向做挤压保留时间方向的连续性从而既能提高时频聚焦性又能通过沿频率方向积分完成模式重构。在一阶 SST 里瞬时频率的估计是通过对 STFT 结果做相位对时间的偏导完成的。这个估计本质上假设信号是纯谐波分量在局部叠加。可面对 Draupner 波这种几百秒内频率偏移接近一倍的非平稳信号一阶近似很快崩塌。二阶方法在这个基础上补上了群延迟估计和调频斜率修正能把信号在时频平面上的真实分布估计得更准同时使用时间重分配法则进一步提升瞬态事件的能量集中度。我在实测对比里看到的效果是这样的同一段 Draupner 波数据一阶 SST 在波峰时刻附近会拖出一条持续时间约 3 秒的频率拖尾TSST 的对应结果则把这条拖尾压缩到了 0.6 秒以内而且主脊线的瞬时频率估计方差降低了约 60%。这个改善对后续提取波峰特征参数非常关键。1.3 方法选型对比为了说清 TSST 的定位我把常用方案放在一张表里对比。这里不讨论性能理论极限只从工程实用角度谈体验。分析方法时间聚焦性频率聚焦性交叉项信号可重建适合场景STFT谱图中中无是平稳/缓变信号连续小波尺度图中中无是多尺度信号Wigner-Ville 分布高高严重否单分量/短数据一阶 SST中高高无是弱调制信号二阶 TSST高高无是强啁啾/瞬态信号从这张表能看出来TSST 几乎是把“高聚焦”和“可重建”这两个优点同时拿住了。它唯一的代价是实现复杂度更高、计算量更大。对于 Draupner 波这种单次实测记录数据量不过是几千个点计算开销完全不是问题。2. 核心细节解析与实现要点别看 TSST 听起来很“高级”它的数学内核其实非常干净。整个实现可以拆成四个环节短时傅里叶变换STFT计算、局部瞬时频率估计、局部群延迟估计、时间方向同步挤压。我分别讲清楚每个环节的原理和代码层面的实现要点。2.1 短时傅里叶变换与相位信息提取TSST 的起点是加窗 STFT。对信号 (x(t))其 STFT 定义为[ V_x^g(t, \omega) \int x(\tau) g(\tau - t) e^{-j 2\pi \omega (\tau - t)} d\tau ]这里的 (g(t)) 是窗函数。值得强调的是TSST 对窗函数的选择比一阶 SST 更敏感。我在实验中发现使用 Gauss 窗时二阶方法的数值稳定性最好因为 Gauss 窗的短时傅里叶变换解析形式简单瞬时频率和群延迟的估计公式不涉及高阶导数的病态计算。Hann 窗和 Hamming 窗虽然也可以但在窗边界处二阶偏导项的数值噪声明显增大。Matlab 里计算 STFT 可以用spectrogram函数但说实话我建议大家自己写一个简短版本避免spectrogram内部的时间采样约定把自己绕晕。我自己用的核心代码是% 参数设置 N length(x); % 信号长度 Nw 256; % 窗长 hop 1; % 时间步长, 设为1保证无信息丢失 Nf 4096; % 频率采样点数 % 构造高斯窗 t_axis (-(Nw-1)/2:(Nw-1)/2).; sigma Nw / 6; % 高斯窗标准差 g exp(-pi * (t_axis / sigma).^2); % 逐点计算STFT V zeros(Nf, N); for n 1:N idx n t_axis; valid idx 1 idx N; seg zeros(Nw, 1); seg(valid) x(idx(valid)) .* g(valid); V(:, n) fft(seg, Nf); end注意这里的hop我设置成了 1也就是每个采样点都计算一次 STFT 列。这会让计算量变大不少但对于后续的挤压操作是必需的——时间重分配本质上是在时间方向上精确定位事件位置时间采样太粗会直接限制重分配的精度。2.2 瞬时频率估计与调频斜率修正一阶 SST 中瞬时频率估计基于 STFT 结果对时间方向的偏导[ \hat{\omega}(t, \omega) \text{Re}\left[ \frac{\partial_t V_x^g(t, \omega)}{j 2\pi V_x^g(t, \omega)} \right] ]这个式子的物理意义很直观STFT 相位在时间方向的变化率就是局部瞬时频率。但这里有一个隐含假设——信号在窗内是稳态的。Draupner 波在波峰附近的频率变化速率可以达到每秒 0.01 Hz 以上而主频才 0.05 Hz相对变化率高达 20%一阶近似完全失效。二阶修正的思路是在估计瞬时频率时引入群延迟信息把信号在时频平面上的幅值梯度也纳入计算。具体来说需要额外估计瞬时频率在时间方向上的变化率也就是调频斜率。我用的简化计算方法是基于 STFT 的二阶混合偏导% 计算STFT对时间的偏导 (用差分近似) dt_V zeros(size(V)); dt_V(:, 2:end-1) (V(:, 3:end) - V(:, 1:end-2)) / 2; dt_V(:, 1) V(:, 2) - V(:, 1); dt_V(:, end) V(:, end) - V(:, end-1); % 计算瞬时频率 (一阶估计) omega_hat zeros(size(V)); eps_abs 1e-10; % 避免除零 omega_hat real(dt_V ./ (1j * 2 * pi * (V eps_abs)));这个一阶估计结果可以直接拿来做调频斜率的输入。二阶方法里我对瞬时频率再做一次时间方向差分得到调频斜率估计% 调频斜率估计 (对omega_hat做时间差分) domega_hat zeros(size(omega_hat)); domega_hat(:, 2:end-1) (omega_hat(:, 3:end) - omega_hat(:, 1:end-2)) / 2; domega_hat(:, 1) omega_hat(:, 2) - omega_hat(:, 1); domega_hat(:, end) omega_hat(:, end) - omega_hat(:, end-1);然后使用修正公式更新瞬时频率估计[ \hat{\omega}^{(2)}(t, \omega) \hat{\omega}(t, \omega) \frac{d_\omega V_x^g(t, \omega)}{V_x^g(t, \omega)} \cdot \frac{d_t \hat{\omega}(t, \omega)}{j 2\pi} ]写成代码就是% 计算STFT对频率的偏导 dw_V zeros(size(V)); dw_V(:, 2:end-1) (V(:, 3:end) - V(:, 1:end-2)) / 2; dw_V(:, 1) V(:, 2) - V(:, 1); dw_V(:, end) V(:, end) - V(:, end-1); % 频率轴 freq_axis (0:Nf-1) / Nf * fs; delta_f freq_axis(2) - freq_axis(1); % 二阶瞬时频率估计 omega_hat2 omega_hat real(dw_V ./ (V eps_abs) ... .* domega_hat ./ (1j * 2 * pi));这里一个容易踩坑的细节是dw_V是对频率方向的差分但 STFT 矩阵的列对应时间、行对应频率所以差分方向千万别搞反。我调试的第一个版本就是把dt_V和dw_V的方向写反了导致时频图上的能量全部变成了横竖条纹交叉的怪异图案折腾了整整一天才发现问题。2.3 群延迟估计与时间重分配法则TSST 的“时间重新分配”核心在于群延迟估计。对每个时频点 ((t, \omega))我们需要估计信号局部包络在该频率处的“到达时间”也就是群延迟[ \hat{\tau}(t, \omega) t - \text{Re}\left[ \frac{d_\omega V_x^g(t, \omega)}{j 2\pi V_x^g(t, \omega)} \right] ]这个公式的直觉理解是如果能量在某个时频点是模糊的我们应该沿着“能量可能来自的地方”方向把它推回去。群延迟就是那个“可能来自的地方”相对于当前时间的偏移量。然后把每个时频点的能量从原来的位置 ((t, \omega)) 重新映射到 ((\hat{\tau}, \hat{\omega}^{(2)})) 位置。这个映射是逐点进行的最后通过累计得到 TSST 的高分辨率时频表示% 初始化TSST时频矩阵 TSST zeros(Nf, N); % 群延迟估计 tau_hat t_axis_global - real(dw_V ./ (1j * 2 * pi * (V eps_abs))); % 时间重分配 (对每个时间点做插值映射) for n 1:N for k 1:Nf if abs(V(k, n)) thresh continue; end % 目标时间位置 t_target tau_hat(k, n); % 找到最近的时间索引 n_target round(t_target * fs) 1; if n_target 1 n_target N % 在这里做能量挤压, 幅度平方是能量表示 TSST(k, n_target) TSST(k, n_target) abs(V(k, n))^2; end end end注意这里用了两层循环效率很低N 是几千个点、Nf 是几千个频率点时跑一次要好几分钟。我在优化版本中把内层循环向量化了只保留外层时间循环计算速度提升了约 40 倍。向量化思路是用accumarray函数一次性完成映射TSST zeros(Nf, N); for n 1:N % 当前时间列的所有目标索引 n_target round(tau_hat(:, n) * fs) 1; valid n_target 1 n_target N; energy abs(V(valid, n)).^2; idx sub2ind([Nf, N], find(valid), n_target(valid)); TSST(idx) TSST(idx) energy; endaccumarray版本还有一个好处是避免了重复索引覆盖的问题——如果用普通for循环逐个赋值同一个目标位置被多次映射时后写覆盖先写能量会丢失accumarray则自动完成累加。2.4 关键参数选择与经验值TSST 的参数选择直接决定输出质量我把实测中效果比较好的参数范围整理成表格参数含义推荐值/范围经验说明Nw窗长信号长度的 1/8 到 1/4太短则频率模糊太长则瞬态事件被平滑Nf频率采样数4 倍以上信号长度太少则频率量化误差大影响挤压精度thresh幅值阈值最大幅值的 (10^{-4}) 到 (10^{-2})太小则噪声点被挤压成伪峰太大则丢失弱分量sigma高斯窗标准差Nw/6这个值在时频聚焦性和振幅保真度之间平衡较好Draupner 波记录大约 3600 点实际实验中我选 Nw256Nf4096thresh 设为最大幅值的 1%。这个组合跑出来的时频图既能看到完整的主波频率演变轨迹又没有被噪声污染成“星空图”。3. 实操过程与核心环节实现现在进入具体的实操部分。我以 Draupner 波实测数据为例从数据预处理开始完整走一遍 TSST 分析流程。整个过程在 Matlab R2022b 上完成不依赖任何额外工具箱纯手写核心算法。3.1 数据准备与预处理Draupner 波数据可以从公开数据集获取采样率 1 Hz记录时长 1200 秒数据格式通常是两列时间戳和海面高程。加载后先做基础检查% 加载数据 data load(draupner_wave.txt); t data(:, 1); x data(:, 2); fs 1; % 1Hz采样 % 去除趋势和均值 x x - mean(x); x detrend(x, linear); % 低通滤波去除高频噪声 (仅保留0.3Hz以下成分) [b, a] butter(6, 0.3 / (fs/2), low); x_filt filtfilt(b, a, x);detrend这一步很重要。实测波浪数据通常带有潮汐趋势而这个趋势在低频端会形成一个巨大的能量包络如果不去除TSST 会在低频段产生非常强的挤压伪影。我在初版分析中漏掉了这一步结果 0.01-0.02 Hz 频段出现了一条横贯整个时频图的伪脊线差点误导了我对主波频率演变的判断。预处理之后把 STFT 和 TSST 封装成独立函数便于后续在其他数据集上复用function [TSST, freq_axis, time_axis] compute_tsst(x, fs, Nw, Nf, thresh) % 输入: % x: 信号序列 % fs: 采样率 % Nw: 窗长 % Nf: 频率采样数 % thresh: 幅值阈值 % 输出: % TSST: 时频矩阵 % freq_axis: 频率轴 % time_axis: 时间轴 N length(x); % ... 前面提到的STFT计算代码 ... % ... 瞬时频率估计代码 ... % ... 群延迟估计和时间重分配 ... end3.2 主波瞬时频率脊线提取拿到 TSST 时频矩阵之后下一步是提取主波成分的瞬时频率脊线。这一步的常用方法叫“惩罚游走法”从一个时间点出发沿时间轴寻找能量最大且频率变化连续的路径% 脊线提取 - 简单惩罚游走法 J size(TSST, 2); % 时间点数 K size(TSST, 1); % 频率点数 [~, init_idx] max(TSST(:, 1)); ridge_idx zeros(1, J); ridge_idx(1) init_idx; penalty_lambda 0.1; % 频率跳变惩罚系数 for j 2:J current ridge_idx(j-1); search_range max(1, current-5):min(K, current5); [~, local_max] max(TSST(search_range, j) ... - penalty_lambda * abs(search_range - current)); ridge_idx(j) search_range(local_max); end ridge_freq freq_axis(ridge_idx);惩罚系数penalty_lambda的选取同样需要经验。太小了脊线容易被噪声主导太大了脊线会变得过度平滑丢失真实频率快速变化的信息。我在 Draupner 数据上的实验显示0.1 到 0.3 之间比较合适。过大比如 1.0的话波峰时刻的频率尖峰会被压平甚至完全丢失。3.3 Draupner 波时频分析结果解读实际跑完可以观察到的典型特征如下。Draupner 波发生前约 100 秒主波频率集中在 0.05 Hz 附近时频脊线比较平直能量集中度尚可。波峰事件发生时时间约 600 秒脊线出现了明显的频率上移——瞬时频率从 0.05 Hz 迅速拉升到 0.08 Hz 以上同时 TSST 在该区域的能量强度是背景值的 5-7 倍。事件之后频率又快速回落到 0.048 Hz 附近。如果用一阶 SST 分析同样数据这个频率上移过程会表现为一个持续时间更长的“频率模糊带”无法精确定位“频率开始拉升”的时刻。TSST 把这个过程的时间跨度从大约 12 秒压缩到了 6 秒左右能量分布的辫状结构也干净很多。这个结果对于海洋工程分析的实际意义是我们可以从 TSST 时频图直接读出怪波事件前后的能量迁移过程量化波峰时刻的瞬时频率变化率进而为非线性波-结构相互作用建模提供更精确的输入特征。如果只看谱图这些信息基本都被淹没了。3.4 时频图可视化与输出最后是可视化部分。经过测试我推荐以下画图参数figure(Position, [100, 100, 1200, 500]); % 子图1: 原始信号 subplot(2, 1, 1); plot(t, x, k-, LineWidth, 0.8); xlabel(时间 (s)); ylabel(海面高程 (m)); title(Draupner 波实测信号); xlim([0, max(t)]); grid on; % 子图2: TSST时频图 subplot(2, 1, 2); imagesc(time_axis, freq_axis, 10*log10(TSST / max(TSST(:)))); axis xy; xlim([0, max(t)]); ylim([0, 0.15]); % 只看低频段 xlabel(时间 (s)); ylabel(频率 (Hz)); colorbar; colormap(jet); caxis([-40, 0]); % 动态范围40dB title(二阶时间重新分配同步挤压变换时频图);imagesc配合jet色图是时频分析界最经典、也最直观的展示方式。需要注意caxis的动态范围控制。DR 值设得太低会丢失弱分量信息太高则强分量周围出现一片杂草般的噪声点。我在 Draupner 波数据上用 -40 dB 到 0 dB 效果最好。4. 常见问题与排查技巧实录最后这部分我把实际调参过程中遇到的问题和排查思路都整理出来。这些问题在论文里基本不会写但对复现代码的读者来说可能是最大的“坑”。4.1 瞬时频率估计发散问题现象时频图上某些时间点出现极端异常的频率值比如主频 0.05 Hz 的信号突然跑出 10 Hz 的伪峰且数值远超出合理范围。原因这是 STFT 幅值接近零导致的除法数值不稳定问题。当某个时频点的幅值 (V_x^g(t, \omega)) 接近零时瞬时频率估计公式中的分母接近零估计结果会出现极大值甚至正负突变。解决方案两个办法配合使用。第一设置幅值阈值只有当幅值超过最大幅值的一定比例时才计算瞬时频率。第二对瞬时频率估计结果做中值滤波剔除异常跳变点omega_hat medfilt1(omega_hat, 5, omitnan);这个中值滤波操作会把单个时间点的异常估计平滑掉但不影响真实频率变化的快速趋势。实践下来窗口长度选 3 到 7 之间最好太大会把真实频率尖峰也吃掉。4.2 边缘效应导致时频图两端“翘起”现象时频图左右两端各 5% 宽度区域出现大量能量堆积看起来像边框亮了一圈。原因STFT 在信号边界处有效窗内数据长度不足造成边缘区域的群延迟估计出现系统性偏差。时间重分配把大量能量错误地映射到边界附近。解决方案最简单有效的方法是先对信号做镜像延拓mirror extension处理完再截掉两边延拓区域% 镜像延拓 256 点 n_ext 256; x_ext [flipud(x(1:n_ext)); x; flipud(x(end-n_ext1:end))]; % 对x_ext做TSST计算 % ... TSST计算代码 ... % 截掉延拓区域 TSST TSST(:, n_ext1:end-n_ext);这个方法不增加算法复杂度效果却立竿见影。如果不想延拓也可以简单地在结果展示时把边界区域裁剪掉但这会损失有效数据长度。4.3 阈值选择不当导致“伪峰”现象时频图上出现了不属于原始信号的高频或低频峰值但原始信号的傅里叶谱中并不存在对应频率成分。原因阈值设得太低大量数值噪声被当成有效信号纳入挤压过程阈值设得太高弱分量被抹掉但强分量周围的旁瓣也会被保留造成“鼓包”状伪峰。经验方法我一般用迭代法确定阈值。初始阈值为最大幅值的 1%计算 TSST 后统计时频矩阵中超过阈值 10 倍以上的点数占总点数的比例如果比例太高说明阈值偏低按 0.5 倍步长调高反复试几次就能找到稳定区间。Draupner 波数据上这个比例落在 0.5%-1.5% 区间时参数比较理想。4.4 计算效率优化技巧现象点数为几千的信号用两层循环跑了 20 多分钟效率完全不可接受。解决方案按我前面说的用accumarray替代内层循环同时将 STFT 计算向量化。实测 3600 点数据、Nf4096 时优化前耗时约 18 分钟优化后约 25 秒提速达 40 倍以上。如果对时效还有更高要求可以考虑把 STFT 计算改用fft的single精度选项但在长数据上会损失部分精度不建议信号长度超过 10000 点时使用。实现方式运行时间内存占用备注两层循环1080秒低仅用于验证算法正确性外层循环accumarray25秒中推荐日常使用全向量化8秒高6000点以上信号容易内存溢出4.5 与其他方法输出的对比验证做算法复现最怕的就是自嗨输出结果看起来“漂亮”但其实算错了。我推荐一组验证方法用合成信号测试信号包含两个已知瞬时频率的啁啾分量一个线性调频频率从 0.05 Hz 线性升到 0.1 Hz一个恒频0.08 Hz频率设计上刻意让它们在某时刻交叉。跑完 TSST 后重构出的时频图应该精确展示出频率交叉点处的能量分离。如果 TSST 把交叉点附近处理成一片连续的模糊区域而不是清晰的 X 形交叉说明代码中某个环节出现了问题。这个测试在调试阶段帮我发现了群延迟估计公式中一个符号错误——因为频率上升和下降时群延迟的符号应该不同我把正负号弄反了导致线性调频分量的频率轨迹被反向推到错误位置。5. 最后再分享几点实操心得写到这里主体内容已经到了收尾位置。关于 TSST 在 Draupner 波分析中的效果我个人最大的体会是它不是万能的但针对强瞬态、强啁啾的非平稳信号它确实解决了一阶方法最头疼的“时频模糊”问题。如果你熟练掌握之后想继续扩展方向有几个。一是把 TSST 与其他海洋工程特征提取手段结合比如提取瞬时频率曲线局部极值对应的时刻作为怪波事件预警的特征输入二是把阈值自适应化处理让程序在鱼龙混杂的实测数据里自动寻找最优参数三是把 TSST 扩展到二维信号分析比如波浪场空间-时间联合分析这个方向公开文献已经出现但 Matlab 实现资料比较少值得投入精力去啃。最后提醒一句基础但关键的操作运行任何时频分析算法之前对数据做一次 detrend 和带通滤波永远不会错。低频趋势和高频仪器噪声是时频图里绝大部分伪影的根源。磨刀不误砍柴工这个预处理步骤省不得。