4. 功率谱估计方法对比与实现
在海洋工程中，我们拿到的往往不是无限长的连续波面记录，而是一段有限时长、离散采样的波浪数据。从这段有限数据出发去估计真实谱，本质上是一个“统计推断”问题：数据越短、谱越宽，估计的方差就越大。经典周期图法虽然计算简单，但其估计方差并不随数据长度增加而减小，谱线剧烈起伏，直接用于设计谱拟合会带来很大偏差。本节围绕三种在工程与科研中常用的改进方法——Welch法、多锥度法（Multitaper）和最大熵法（MEM）——讨论它们的原理差异、参数选择对结果的影响，并给出可复现的实现思路与对比判读方法。
三种方法的核心思想
Welch法属于“分段平均”路线。它把长度为 N 的序列分成若干段（可重叠），对每段加窗后计算周期图，再对所有段的周期图取平均。平均使方差下降，代价是频率分辨率由段长决定。设段长为 L、重叠率为 50%，则分段数约为 K≈2N/L，方差大致降为单段周期图的 1/K。工程上常用 Hann 窗配合 50% 重叠，因为 Hann 窗旁瓣衰减快，重叠一半可补偿窗函数两端的信息损失。
多锥度法走的是“正交锥度加权”路线。它不切分数据，而是对整段数据同时施加 K 个相互正交的 Slepian 锥度序列（离散扁球序列），得到 K 个近似独立的谱估计再平均。其关键优势是：在给定分辨率带宽 NW 下，能在“低泄漏”和“低方差”之间取得很好的折中，特别适合短记录、宽峰或存在强低频能量的海浪数据。参数上，时间带宽积 NW 通常取 2.5～4，锥度数 K 取 2NW−1 左右，K 越大方差越小但偏差略增。
最大熵法则是“参数模型”路线。它假设序列可由一个 p 阶自回归（AR）模型生成，通过 Burg 法或 Yule-Walker 方程求解 AR 系数，再由系数直接写出谱密度。MEM 的突出特点是对短数据、窄带峰的分辨能力很强，谱峰尖锐、无窗函数泄漏；但阶数 p 选择敏感：p 过小谱过于平滑、峰被抹平，p 过大则出现虚假分裂峰。对海浪这类宽峰谱，p 一般取 N/3～N/2 量级需谨慎，实际更常用 20～60 的适中阶数并配合稳定性检验。
参数选择如何改变谱形
三种方法的参数敏感性各不相同，理解这一点比记住“用哪个函数”更重要。
对 Welch 法，段长 L 是分辨率与方差的“总开关”。若采样频率 fs=2 Hz、目标分辨 0.01 Hz，则需 L≈100 s，即约 200 个采样点；若记录总长只有 600 s，则只能分 6 段左右，方差改善有限。此时若强行增大 L 换取分辨率，段数骤减，谱线又会变得毛糙。实际做法是先按目标谱峰宽度确定 L，再检查剩余段数是否足以把方差压到可接受范围；不足时改用多锥度法往往更划算。
对多锥度法，NW 决定分辨率带宽，K 决定方差。NW 增大使分辨率带宽变宽、谱更平滑；K 增大进一步降方差，但当 K 接近 2NW 时各锥度间相关性上升，收益递减。海浪谱主峰通常较宽，取 NW=3、K=5 是稳健起点；若关注的是涌浪窄峰，可适当减小 NW 以保留峰形。
对最大熵法，阶数 p 是最需要试验的参数。建议做法是：对 p 从 10 到 80 扫描，观察谱峰位置是否稳定、是否出现无物理依据的分裂峰，并结合 AR 模型反射系数的收敛情况判断。若 p 增大后主峰位置漂移或冒出多个尖峰，说明已过拟合，应回退。
可复现的实现示例
下面以 Python 为例给出三种方法的核心调用。示例参数为教学假设：fs=2 Hz，记录 1024 点（约 512 s），目标分辨约 0.02 Hz。
import numpy as np
from scipy import signal
fs = 2.0
x = ... # 去均值、去线性趋势后的波面序列，长度 1024
# 1) Welch：段长 128 点 → 分辨约 0.0156 Hz，50% 重叠
f_w, P_w = signal.welch(x, fs=fs, window='hann',
nperseg=128, noverlap=64,
detrend='linear', scaling='density')
# 2) 多锥度：NW=3, K=5
from scipy.signal.windows import dpss
NW, K = 3.0, 5
tapers = dpss(len(x), NW, Kmax=K)
P_mt = np.zeros(len(x)//2 + 1)
for w in tapers:
f_m, P_seg = signal.periodogram(x * w, fs=fs, scaling='density')
P_mt += P_seg
P_mt /= K
# 3) 最大熵：Burg 法求 AR 系数，再算谱
from spectrum import pburg, arma2psd
p = 40
ar_model = pburg(x, p, NFFT=1024, sampling=fs)
f_mem, P_mem = ar_model.frequencies(), ar_model.psd实现时有几个易错点值得强调：Welch 与多锥度都应先去除均值与线性趋势，否则低频段会被趋势能量严重污染；多锥度法必须对每个锥度分别做周期图再平均，不能先平均锥度再加窗；MEM 得到的谱是连续解析式，需在关心的频段上离散采样，且其谱密度单位要与 Welch 的 density 缩放保持一致，否则三者无法直接叠图比较。
实战对比与结果判读
把三条谱曲线画在同一张半对数图上，按以下顺序判读：先看主峰频率是否一致，若某方法峰位明显偏移，多半是分辨率或阶数设置不当；再看峰高与峰宽，Welch 通常峰最平缓，MEM 峰最尖锐，多锥度居中；最后看高频尾部与低频端，若某方法在低频异常抬高，应怀疑泄漏或趋势未去除。对海浪数据，还建议把估计谱与 JONSWAP 或 Pierson-Moskowitz 目标谱叠加，观察偏差集中在哪个频段——这比单纯比较曲线形状更有工程意义。