📚 课程目录第3章波浪时间序列的频域分析基础

3. 波浪时间序列的频域分析基础

浮标或测波雷达记录到的是波面高度随时间变化的序列 η(t)，它看起来杂乱无章，但海洋工程真正关心的往往是“能量集中在哪些周期上”。频域分析就是把这条时间曲线拆解成不同频率正弦波的叠加，从而得到波浪能量随频率的分布。本节围绕傅里叶变换原理、功率谱密度定义、周期图法、窗函数选择与谱泄漏抑制展开，并给出可直接在 Python 中复现的示例流程。

从时域到频域：离散傅里叶变换

设采样频率为 f_s，采样间隔 Δt = 1/f_s，一段长度为 N 的波面序列为 η[n]，n = 0,1,…,N−1。离散傅里叶变换（DFT）定义为：

X[k] = Σ(n=0→N−1) η[n] · exp(−j·2πkn/N),  k = 0,1,…,N−1

X[k] 是复数，其模值 |X[k]| 表示第 k 个频点上的幅值贡献，相位 arg X[k] 表示该分量的初相。第 k 个频点对应的物理频率为 f_k = k·f_s/N，频率分辨率 Δf = f_s/N = 1/T，其中 T = N·Δt 为记录时长。这带来一个关键工程约束：要分辨 0.01 Hz 的频率间隔，记录长度至少需要 100 s；而最高可分析频率受奈奎斯特频率 f_s/2 限制，因此采样率必须大于波浪有效最高频率的两倍。

实际计算使用快速傅里叶变换（FFT），其结果为双边谱，包含正负频率。对实值信号，负频率部分与正频率共轭对称，物理上只取 k = 0 到 N/2 的单边谱，并把非直流分量的幅值乘以 2 以保持能量一致。

功率谱密度的定义与量纲

功率谱密度 S(f) 描述单位频率带宽内的波面方差贡献，满足：

m0 = ∫(0→∞) S(f) df = Var(η) = σ²

其中 m0 是谱的零阶矩，等于波面方差 σ²。由 m0 可导出有效波高 H_s = 4√m0，这是谱分析与波浪统计之间的桥梁。S(f) 的量纲为 m²/Hz（若 η 以米计），而离散谱估计值 S_k 的量纲为 m²，二者关系为 S(f_k) ≈ S_k/Δf。很多初学者把离散谱值直接当作谱密度使用，导致 H_s 计算偏差，这是必须避免的错误。

周期图法：最直接的谱估计

周期图法（Periodogram）直接用 FFT 结果的平方作为谱估计：

S_k = (Δt / N) · |X[k]|²,  k = 1,…,N/2−1
S_k(单边) = 2 · S_k

周期图法实现简单，但它是不一致估计：随着 N 增大，谱估计的方差并不会减小，谱线依然剧烈起伏。工程上常用两种改进：一是分段平均（Welch 法），把长序列分成若干重叠段，分别求周期图再平均，以降低方差、牺牲部分分辨率；二是平滑窗（如 Hanning 窗加权平均），在频域对相邻谱线做滑动平均。示例中若采样率 f_s = 2 Hz、记录 1024 s（N = 2048），则 Δf ≈ 0.00098 Hz，分段为 8 段、50% 重叠后每段 512 点，Δf 变为约 0.0039 Hz，方差显著下降，适合观察谱峰形状。

窗函数与谱泄漏抑制

DFT 隐含假设序列是周期延拓的。若记录首尾不连续，截断处会产生突变，能量从真实频率“泄漏”到邻近频点，形成谱泄漏，表现为谱峰变宽、旁瓣抬高，甚至淹没弱分量。解决办法是在 FFT 前乘以窗函数 w[n]，使两端平滑衰减到零。

常用窗函数各有取舍：矩形窗主瓣最窄、频率分辨率最高，但旁瓣最高（约 −13 dB），泄漏最严重；Hanning 窗旁瓣约 −31 dB，主瓣稍宽，适合一般海浪谱分析；Hamming 窗旁瓣约 −43 dB，但远端旁瓣衰减慢；Blackman 窗旁瓣可低至 −58 dB，主瓣更宽，适合强动态范围场景。选择原则是：关注谱峰位置与形状时优先 Hanning；需要压制强分量旁瓣对弱分量的干扰时用 Blackman；做幅值精确校准时注意窗函数的相干增益修正（Hanning 窗相干增益约 0.5，需对幅值乘以 2 补偿）。

上图把本节的处理链与谱泄漏的直观差异放在一起：左侧是时域序列经加窗、FFT 到谱密度的主流程，下方对比了矩形窗与 Hanning 窗在同一谱峰上的表现。可以看到，未加窗时谱峰两侧出现明显旁瓣，弱分量容易被掩盖；加窗后旁瓣被压低，谱峰轮廓更接近真实能量分布。

一个可复现的示例流程

假设有一段模拟波面数据（示例参数，非实测）：采样率 2 Hz，时长 1024 s，主周期约 8 s（对应频率约 0.125 Hz），叠加少量 4 s 分量与随机噪声。可按以下步骤处理：

import numpy as np
from scipy.signal import welch

fs = 2.0
t = np.arange(0, 1024, 1/fs)
eta = (1.2*np.sin(2*np.pi*0.125*t)
       + 0.4*np.sin(2*np.pi*0.25*t + 0.7)
       + 0.15*np.random.randn(len(t)))

# Welch 法：Hanning 窗，段长 512，50% 重叠
f, S = welch(eta, fs=fs, window='hann',
             nperseg=512, noverlap=256,
             detrend='constant', scaling='density')

m0 = np.trapz(S, f)
Hs = 4*np.sqrt(m0)
print(f"Δf = {f[1]-f[0]:.4f} Hz, m0 = {m0:.3f} m², Hs = {Hs:.2f} m")

运行后应观察到：谱峰出现在约 0.125 Hz 与 0.25 Hz 处，前者明显更高；由积分得到的 m0 与输入分量的方差之和（约 0.72 + 0.08 + 噪声方差）量级一致，H_s 落在合理范围。若把 window 改为 'boxcar'，会看到谱峰两侧旁瓣抬高、峰形变宽，这正是谱泄漏的直接证据。判断谱估计是否可信，可检查三点：谱峰频率是否与目视周期吻合、m0 是否与波面方差接近、分段数增加后谱形是否趋于稳定。

做谱分析前先去除均值（detrend），否则直流分量会抬高低频端；同时确认采样率满足奈奎斯特条件，避免高频折叠污染谱形。若记录较短，优先保证频率分辨率，不要盲目增加分段数。

掌握上述原理后，后续的波浪谱模型拟合、谱矩计算与设计波参数推导，都建立在这一步得到的可靠 S(f) 之上。频域分析的价值不在于画出好看的曲线，而在于让每一个谱峰、每一段能量都能追溯到明确的物理频率与统计量。