在海洋工程和船舶设计里，随机波浪谱分析是一门绕不开的基础功课。拿到一条波面升高时间序列，如果只做时域统计，等于把五颜六色的频谱强行压成一团灰度。同样的有效波高和方差，能量集中在0.1Hz附近的涌浪与散布在0.3~0.5Hz的风浪，作用到船体和平台上是完全不同的两种结果。所以把时间序列转成频率上的能量分布，是耐波性计算、疲劳分析、波浪能装置设计的第一步。这篇文章要讲的方法叫自相关函数法，也就是常说的Blackman-Tukey法思路：先算波面序列的自相关，再做傅里叶变换得到功率谱密度。我会给出一套能在MATLAB里直接跑通的核心代码，覆盖数据预处理、谱估计、理论谱验证和参数提取，再把实际项目中容易踩的坑一并交代清楚。

1. 随机波浪为什么非做谱分析不可：从时域到频域的思考

1.1 一条波高记录到底能看出什么

假设某个海上测站记录了60分钟波面，采样率4Hz，得到一万四千多个数据点。把它画出来，大部分人会得出三个结论：最大值大概几米、均值接近零点、序列像乱码一样上下跳。但这些结论对工程判断几乎没用。真正决定结构响应的，是波浪能量在频域上的分布。

举个例子：有效波高都是2.5m的两段波浪，第一段是远距离涌浪，能量集中在0.08~0.12Hz；第二段是当地风浪，能量分布在0.2~0.5Hz。对一条浮式结构来说，前者可能激发明显的低频共振，后者则主要影响局部构件疲劳。两段波的时域图放在一起，肉眼几乎分辨不出来，都是"上上下下的曲线"，但频域谱画出来之后差异极其明显。

这也是为什么国际通用的设计波浪谱都会给出谱型参数，而不只是给一个有效波高。海上平台桩基疲劳、船舶耐波性评估、波浪能装置年发电量估算，都要靠频谱形状去推算。如果只拿时域统计算完Hs就收工，等于把一个随机过程的全部信息压成一个数字，后续共振频率和疲劳热点全都无从谈起。

1.2 频域谱在工程里具体贡献了什么

功率谱密度在工程计算里几乎是标准输入。最直接的用途是耐波性频域分析：把波浪谱乘上响应幅值算子（RAO），就能得到运动或载荷的响应谱，再统计预报极值运动幅值和加速度。疲劳分析里，谱矩方法——也就是m0、m1、m2这些频率矩——可以直接估算疲劳损伤率，省去逐波浪计数的大量时域模拟。波浪能装置设计需要知道有效波高和能谱周期的联合分布，才能估算年度产能。系泊与定位分析也要靠谱密度推算低频慢漂力。

可能有人会问：现在MATLAB里不是有pwelch这类现成函数吗，为什么还要自己写自相关函数法？我的体会是：pwelch的默认窗、默认分段和默认归一化都会影响结果，单位问题也常让人看懵。自相关函数法虽然代码多几行，但每一步都透明可控，截断比例、窗函数、归一化系数都能清楚知道作用在哪。而且自己实现一遍之后，再回头看内置工具的输出就不会被黑盒带到沟里去。

2. 自相关函数法的核心原理：维纳-辛钦定理与有偏估计的取舍

2.1 自相关函数为什么能“看见”周期成分

自相关函数的定义是 R(τ)=E[η(t)η(t+τ)]。对随机波浪来说，就是把一条记录在时间轴上错开一个滞后τ后，与原始序列相乘再求平均。这个值高，说明相距τ的一对采样点在统计上步调一致；值低，说明它们之间基本没有规律性关联。

随机海浪不是白噪声，它内部存在与波浪周期相关的记忆性，因此自相关函数会在某些滞后位置呈现峰值。更直观地理解：任何一个窄带随机过程，看起来都像是按某个主周期振荡、幅度缓慢变化的波形，自相关函数保留了这个主周期信息。算出自相关之后再做傅里叶变换，就能把埋藏在时间统计里的周期成分映射成频域上的能量集中区，这正是维纳-辛钦定理做的事。

实际估计自相关时，我们用的是离散形式：R[m] = (1/N)·Σ η[n]·η[n+|m|]，求和范围是n=0到N-|m|-1。注意参与平均的项数只有N-|m|个，滞后越大，尾部噪声越明显。这个细节在后面选有偏还是无偏时会起决定性作用。

2.2 维纳-辛钦定理和谱的物理单位

功率谱密度可以写成 S(f)=∫R(τ)e^{-i2πfτ}dτ，逆变换为 R(τ)=∫S(f)e^{+i2πfτ}df。对实数平稳过程，S(f)是偶函数，一般只看单边正频率并把幅度加倍。

单位是最容易出错的地方。自相关R(τ)的单位是m²，τ单位是秒，积分之后S(f)的量纲就是m²·s，等价于m²/Hz。频率轴上如果取等差Δf，那么每根谱线面积近似为S(f)Δf，把所有谱线面积加起来，应该等于波形平方的平均值，去均值后就是方差。这条面积守恒性质是验证程序对错最有力的工具。

谱值的物理含义也别搞混：S(f)df是频率落在f到f+df之间的波能，谱值越大，该频段波浪能量越大。谱值单位是m²/Hz，比如0.1Hz处有1 m²/Hz的谱值时，该频带对总方差的贡献要乘带宽，而不是把1直接当成波高平方。初学者拿谱峰值去和Hs比大小，常被数值吓一跳，本质就是单位没理清。

2.3 有偏估计和无偏估计到底该选谁

用MATLAB的xcorr时，可以指定'biased'或'unbiased'。有偏估计统一除以N，无偏估计在每个滞后分别除以有效样本数N-|m|。听起来无偏更"准确"，为什么工程上默认用有偏？

因为滞后越大，N-|m|越小，单点估计的方差会急剧放大。举一个具体场景：N=10000，滞后到8000时，只剩下2000个乘积项参与平均，噪声导致的波动会被成倍放大。把这些高方差自相关值拿去做傅里叶变换，谱估计会被尾部随机抖动污染，曲线出现大量假峰。有偏估计虽然在大滞后处引入了偏差，但整体方差小得多，总误差通常更可控。

小样本且截断长度很短时，两者差别不大。但由于本文后面要把自相关截断到N/4~N/2的量级，尾部样本数量仍然不小，我建议统一使用有偏估计，这也是自相关函数法里的常见约定。

3. MATLAB实现：从波面序列到功率谱密度

3.1 数据预处理：均值、趋势和野值别急着跳过

实测数据拿到手，不能直接扔进xcorr。第一件事是剔除野值、去均值和按需去趋势。

野值通常来自测波仪器偶发噪声，表现为单点或多个点突然跳到异常大或异常小的值。简单有效的做法是用局部中位数做基准，偏离超过3倍标准差就判为野值，再用线性插值填补。第二件事是去均值，这一步必须做。如果不减均值，R(0)里包含均值平方，会使谱在零频附近出现一个巨大的尖峰，并泄漏到相邻低频段。

eta_raw = load('buoy_eta.mat').eta; % 实测波面序列，单位m fs = 4; % 采样率 dt = 1 / fs; t = (0:length(eta_raw)-1) * dt; % 1) 野值剔除：与局部中位数偏离过大则插值填补 med0 = movmedian(eta_raw, 201); % 201点左右局部中位数 sd0 = std(eta_raw - med0); bad = abs(eta_raw - med0) > 3 * sd0; eta_raw(bad) = nan; eta = fillmissing(eta_raw, 'linear'); % 2) 去均值 eta = eta - mean(eta);

趋势处理要更谨慎。仪器漂移可能带来线性趋势，但真实的缓变水位也表现为低频信号。我的习惯是先做均值去除，再看时域图判断是否存在明显漂移；确认是仪器漂移才去趋势，否则保留低频真实信息。去趋势可以用polyfit拟合一次多项式后减去，但对天然含涌浪成分的波浪数据，过度去趋势会把真实涌浪能量削减掉。

3.2 核心函数：自相关、截断加窗、FFT

下面是自相关函数法的主函数。逻辑并不复杂：先算有偏自相关，再截断加窗，最后对称延拓做FFT。

function [f, S] = waveSpectrumAutocorr(eta, fs, Lratio) % 基于自相关函数法的随机波浪谱分析 % eta : 波面升高序列，单位m % fs : 采样率，单位Hz % Lratio: 自相关截断比例，默认0.5 eta = eta(:); eta = eta - mean(eta); N = length(eta); dt = 1 / fs; if nargin < 3 Lratio = 0.5; end % 计算有偏自相关，输出长度为2N-1 R = xcorr(eta, 'biased'); R_pos = R(N:end); % 取lag=0到N-1的单边部分 % 截断并加窗 L = max(10, round(Lratio * N)); L = min(L, N - 1); R_trunc = R_pos(1:L+1); w = hamming(L+1); w = w / w(1); % 保证lag=0处权重为1 R_trunc = R_trunc .* w; % 对称延拓后FFT R_sym = [flip(R_trunc(2:end)), R_trunc]; nfft = 2^nextpow2(length(R_sym)); S_fft = abs(fft(R_sym, nfft)) * dt; % 转单边谱 nhalf = nfft / 2; S = S_fft(1:nhalf+1); S(2:end-1) = 2 * S(2:end-1); f = (0:nhalf) * fs / nfft; end

有几个点值得展开。hamming窗在lag=0处的值不是1，而约等于0.08，直接乘会把R(0)压小，算出来的总方差会严重偏低。所以必须执行w = w / w(1)，把窗在零滞后处的权重归一为1，这样才能保证谱面积守恒。

R_sym的构造也很关键：把单边自相关R0,R1,...,RL翻折成RL,...,R1,R0,R1,...,RL的偶序列，再做FFT。由于取绝对值，FFT的相位偏移不会影响幅度谱，频率轴的起点也能对应到0频。零填充到nfft只是让频率网格更细，并不会提高真实分辨率，真实分辨率由截断长度L决定，近似为fs/L。

3.3 单边谱归一化：2倍因子和频率轴

函数返回的S已经是单边功率谱密度，单位m²/Hz。实数随机信号的谱关于零频对称，负频率和正频率各占一半能量，因此只看正频率时要把幅度加倍。但直流分量f=0不翻倍，Nyquist频率处也要按约定处理，代码里S(2:end-1)=2*S(2:end-1)正好避开了这两个端点。

调用之后一定要做面积校验：

[f, S] = waveSpectrumAutocorr(eta, fs, 0.25); df = f(2) - f(1); m0 = sum(S) * df; fprintf('m0 = %.6f, var(eta) = %.6f\n', m0, var(eta));

如果两个值接近，说明实现正确。如果m0明显小于var(eta)，第一嫌疑是窗没做零滞后归一化；第二嫌疑是单边谱的2倍因子漏了。这个方法我几乎每个项目都会跑一遍，比看曲线是否符合理论谱更直接。

4. 用理论谱检验实现：JONSWAP谱对比与谱参数提取

4.1 先用合成波浪验证算法

不要一上来就分析实测数据。先用已知谱形的目标谱生成人工波面，验证算法能还原出同样的谱形。生成办法是线性叠加法：把目标谱分成若干频段，每个频段生成一个余弦波，振幅取sqrt(2·S(f)·Δf)，相位随机。

Hs = 2.5; Tp = 8.0; fs = 4; dt = 1/fs; Tdur = 3600; t = (0:dt:Tdur-dt).'; Nt = length(t); fmin = 0.01; fmax = 0.6; nf = 500; freqs = linspace(fmin, fmax, nf).'; df = freqs(2) - freqs(1); % 简化JONSWAP谱（验证用，非严格标定版本） fp = 1 / Tp; sigma = 0.07 * (freqs <= fp) + 0.09 * (freqs > fp); S0 = 0.3125 * Hs^2 * Tp .* (freqs / fp).^(-5) .* exp(-1.25 * (freqs / fp).^(-4)); gamma = 3.3; S_target = S0 .* gamma .^ exp(-0.5 * ((freqs - fp) ./ (sigma * fp)).^2); S_target(freqs <= 0) = 0; amp = sqrt(2 * S_target * df); phase = 2 * pi * rand(nf, 1); eta = zeros(Nt, 1); for k = 1:nf eta = eta + amp(k) * cos(2*pi*freqs(k)*t + phase(k)); end [f_est, S_est] = waveSpectrumAutocorr(eta, fs, 0.25); figure; loglog(freqs(freqs > 0.01), S_target(freqs > 0.01), 'r-', ... f_est(f_est > 0.01), S_est(f_est > 0.01), 'b-'); xlabel('f (Hz)'); ylabel('S(f) (m^2/Hz)'); legend('理论JONSWAP', '自相关法估计');

注意：这个JONSWAP写法是简化验证版本，能量尺度因子没有做严格归一化，生成的Hs会有一定偏差，但不影响检验谱形和谱峰位置。正式工程中要换成自己项目标定的谱型函数。

曲线大体吻合，尤其谱峰位置和总趋势一致，算法就算通过了。更严格的做法是比较m0与var(eta)，再检查Hs估计值是否接近目标值。

4.2 从谱曲线提取工程参数

随机波浪谱分析最终要落到几个工程参数上。以频率f（单位Hz）为自变量定义谱矩：m_k = Σ f^k·S(f)·Δf。

常用参数有这几个：有效波高Hs=4√m0；谱峰周期Tp=1/f_peak；平均能量周期Tm01=m0/m1；平均跨零周期Tz=√(m0/m2)。

df = f_est(2) - f_est(1); m0 = sum(S_est) * df; m1 = sum(f_est .* S_est) * df; m2 = sum(f_est.^2 .* S_est) * df; Hs_est = 4 * sqrt(m0); [~, imax] = max(S_est); Tp_est = 1 / f_est(imax); Tm01_est = m0 / m1; Tz_est = sqrt(m0 / m2);

验证时，Hs_est会接近时域统计的4倍标准差。如果偏差明显，先查谱面积守恒，再查频率网格是否覆盖了主要能量频带。m2对高频段非常敏感，谱在Nyquist附近的数值噪声会被平方放大。如果Tz明显偏小，说明高频端有噪声尾巴，应该在积分前截到合理频带，比如只积到0.5Hz。

双峰谱的情况要额外小心。涌浪峰和风浪峰分开时，总谱矩参数仍然可算，但峰周期只取全局最大峰容易忽略次峰。建议分别提取两个峰的位置和频带范围，再结合工程目的决定用哪个峰周期。

4.3 谱估计的置信度：单点尖峰别当真

谱估计本质上是随机量。直接FFT周期图的相对标准差非常大，可以接近100%；自相关函数法通过截断与加窗压低了方差，但代价是频率分辨率下降。截断比例Lratio取0.25时，等效自由度大致为2N/L，长记录下相对标准差仍然不小。

所以谱曲线上出现细尖峰时，不要立刻当成确定的物理信号。我处理数据时多次遇到这样的情况：某频率处出现一个明显谱峰，看着像是涌浪能量集中，但换一段记录重算之后峰完全消失。判断谱峰是否可信，至少要看它能否在多个分段子样本中稳定出现，或者跟理论谱形态做对比。更严格的做法是用χ²分布计算置信区间，工程上先把握这个判断原则比背公式更重要。

5. 实测数据分析的常见坑：采样、截断、平滑与单位换算

5.1 采样率与记录时长怎么匹配

波浪谱的工程关注范围通常在0.02~0.5Hz，Nyquist频率留2~3倍裕量，采样率2~4Hz已经足够。常见误区是采样率拉得很高、记录却很短。比如fs=10Hz只记录5分钟，低频端只能分辨到1/300Hz，大约0.0033Hz，这并不算问题；真正的问题是0.05Hz的涌浪峰需要至少20秒以上的时间窗口才能分开，5分钟记录对更低频成分几乎无法有效分辨。

经验值是记录时长至少包含200个跨零周期。风浪典型周期6~10s，对应20~30分钟以上；涌浪周期更长，可能需要1小时记录。如果只有短记录，就不要试图解析过低频段的细节，否则算出来的谱峰周期误差大到没法用。

采样率与记录时长本质上是互斥的：要低频分辨就拉长记录，要高频覆盖就提高采样率。接到数据先问清楚用户关心的频段是什么，再定处理策略，不要盲目追求高采样率。

5.2 窗函数和平滑的平衡：分辨率与方差的取舍

Lratio是自相关函数法里最关键的旋钮。L取得大，相关滞后保留多，频率分辨率高，但尾部方差也大；L取得小，曲线更平滑，但峰会被抹平，相邻频段的能量互相泄漏。

我的一般取值在0.25~0.5。谱看起来毛刺太多就降到0.2；峰被压得很平就升到0.6。窗函数方面，Hamming和Hann比矩形窗旁瓣低，峰更干净，但主瓣稍宽。做疲劳和极值预报时，主峰被抹平往往比小旁瓣泄漏更有害，所以默认Hamming就够了，不建议用Blackman这类主瓣太宽的窗。

平滑思路还有另一条路：分段Welch平均。把长记录切成多段分别算周期图再平均，可以有效降方差。工程中也可以先做分段Welch，再对结果做自相关法意义的平滑，曲线会很稳，只是计算链条长一些。

5.3 零频、趋势与真实低频成分的取舍

实测波面数据里往往混入水位变化，比如潮汐和风增水。直接对原始序列做谱分析，零频附近会出现巨大能量，并泄漏到0.05Hz以下，把真实涌浪峰淹没。

随机波浪分析关心的是扰动波，频率下限通常取0.02Hz以上。潮汐周期以小时计，对应频率约0.0001Hz量级，属于趋势项，应滤波掉。我的做法是先做高通滤波，截止频率设在0.01~0.02Hz，然后再做谱分析。如果MATLAB有Signal Processing Toolbox，可以用highpass，没有的话自己搭一个二阶IIR高通。要注意的是，滤波会带来瞬态效应，处理后的首尾各弃掉约10倍截止周期长度的数据，否则边缘振荡会污染谱估计。如果记录足够长且只关心涌浪以上频段，简单去均值也能接受。

5.4 周期图法、Welch法和自相关函数法的关系

最后把几条谱估计路线理顺。直接FFT周期图是最原生的无平滑估计，方差极大；Welch法通过分段平均降方差；自相关函数法通过截断自相关来实现平滑。两者在长数据极限下估计的是同一个功率谱密度，差别主要在平滑策略。

工程上，用pwelch算得更快，但自相关函数法每一步都可控可解释，适合排查归一化和泄漏问题。我建议的实践是：同一段数据分别用自相关函数法和pwelch算一遍，叠图对比。如果趋势一致但数值相差常数倍，几乎可以肯定是单边谱归一化或窗能量校正出了问题。如果曲线形状都对不上，再回头检查数据预处理和截断参数。这个方法看起来多一道工序，但能省下大量基于错误谱做后续计算的返工时间。

我自己把这个算法真正跑通之后，习惯在每个项目里保留一张校验图：把估计谱和对应的理论JONSWAP谱叠在log-log坐标里，低频异常泄漏和高频噪声尾巴一眼就能看出来。凡是Hs偏差超过5%，或者谱峰周期偏差超过0.5s，我都会先回头查预处理和窗归一化，而不是急着调平滑参数。这套流程说起来简单，做起来需要耐心，但一旦把"自相关→谱→参数"的链路摸清，后面换数据集、换工程对象，都只是改改输入和单位的事。