简介一套基于平稳随机过程理论开展船舶海上横摇、纵摇运动仿真的MATLAB资料包面向船舶海洋工程、航海技术等专业本科生与硕士生适用于课程设计、毕业设计及教研验证场景。包内共8个文件5个Simulink模型.mdl用于构建船舶运动仿真子系统与网络化仿真模块2个MATLAB脚本.m承担数据处理与仿真流程驱动另附1个txt说明文档便于快速上手整体压缩包仅27KBMATLAB 2019a环境即可运行。目前已有529人学习对初学者而言可对照随机海浪作用下的横摇、纵摇响应分析练习理解平稳随机过程理论如何落实到仿真实现。通过梳理各模型与脚本的调用关系还能掌握运动仿真建模思路、参数设置与结果分析流程为后续船舶耐波性与波浪载荷研究奠定基础。1. 平稳随机过程算船舶横摇纵摇到底在算什么拿到一个耐波性任务第一步通常不是仿真而是先回答一个很实际的问题这条船在 5 级海况下横摇有义单幅是多少度、纵摇多少度、多长时间摇一个来回。海上实际波浪是不规则波一条条规则波去叠既不现实也不可靠工程上把这问题归到平稳随机过程框架里把海浪看作平稳随机过程用谱分析一次算出横摇、纵摇的统计特征值。这个方案能给出有义单幅、平均跨零周期等设计输入但给不了“下一分钟具体摇到几度”这种瞬时值做不了也要知道边界在哪。适合船型开发、舾装设计、海事软件配套的工程师也适合刚入门想做耐波性基线的学生。2. 把海浪当成平稳随机过程从波谱、RAO到响应谱的三段式2.1 为什么“平稳随机”够用波面位移与各态历经性海上波浪在严格意义上是非平稳的风在变、潮在变、船在动。但工程处理有一个“短期海况”假设在 20 分钟到 3 小时这个时间窗口内认为波高统计特征不随时间平移改变均值、方差、相关函数都稳定这就是平稳性再把大量不同频率、随机相位的余弦波叠加成波面得到的就是平稳随机过程的一个现实样本。真正让谱分析能在 MATLAB 里落地的是各态历经性。意思是我不用对同一海况做几百次海洋试验取集合平均只要取一条足够长的时历做时间平均就能近似得到集合平均的结果。这就是用一段 20 分钟模拟时历验证谱矩能成立的根据。窄带随机过程概念在这里很关键横摇、纵摇响应谱往往集中在一个窄频带里时历看起来像“载波频率不变、幅值缓慢变化”的调幅波峰值包络服从 Rayleigh 分布。有义单幅值、最大预期值都是从这条假设延伸出去的。这套框架成立有两个前提一是海况参数在这个窗口里稳定二是船的运动响应是线性的。大角度横摇超过 15~20 度时非线性阻尼会改变谱形纯谱分析就不够准这点后面避坑章节还会提到。2.2 运动响应谱的三段式波谱乘RAO模方船舶对不规则波的响应在平稳随机过程理论里不是直接“波高乘一个系数”而是三个环节串起来海浪谱 S_wave(ω)描述波能按圆频率 ω 的分布单位是 m²·s/rad。RAO幅值响应算子H(iω)船在单位波幅规则波下的运动幅值响应取决于船型、航速、浪向。响应谱 S_response(ω)横摇、纵摇运动的能量谱等于波谱乘 RAO 模的平方。写成式子就是最常用的三段式S_response(ω) S_wave(ω) · |H(iω)|²这里有一个新手常忽略的细节|H(iω)|² 是幅值平方不是 H(iω) 直接乘。RAO 通常由势流软件如 WAMIT、SESAM、AQWA算出来或者由模型试验得到横摇、纵摇都有自己的实部和虚部或幅值、相位。如果你拿到的是复数形式的 RAO做谱分析时只需要幅值部分相位角在谱分析里被丢掉了这不是失误因为功率谱本身只保留幅值信息只有合成时历做短时预报时才需要相位。响应谱算出来后没有直接“读”出运动幅值。要从谱矩入手m_n ∫₀^∞ ωⁿ · S_response(ω) dω其中常用的是零阶矩 m₀ 和二阶矩 m₂。有义单幅值横摇角、纵摇角习惯用单幅就等于 2√m₀平均跨零周期 T_z 2π√(m₀/m₂)。注意波浪的有义波高是 4√m₀因为波高是峰谷差、是双幅运动响应大家习惯报单幅所以用 2√m₀。这个单位约定如果混了结果会差整整一倍。2.3 波谱选型P-M谱、ITTC双参数谱怎么给参数海浪谱的选取直接决定响应谱的绝对量级是整条链路的起点。工程上最常用的两类Pierson-Moskowitz 谱P-M 谱用于充分成长的风浪输入是风速。它以风速 U₁₉.₅ 为参数谱形简单适合开阔海域。但它只有一个输入没法同时控制波高和周期遇到海况预报给了“有义波高 周期”两个量时直接用 P-M 谱会有偏差。ITTC 双参数谱是 P-M 谱的工程改写输入换成更容易得到的两个量有义波高 Hs 和特征周期 T₁。谱形如下S_wave(ω) (Hs² / 4π) · (2π/T₁)⁴ · ω⁻⁵ · exp( −(1/π) · (2π/(T₁·ω))⁴ )这个谱的好处是 Hs 和 T₁ 解耦海况预报、统计资料基本都按这两个量给我实际工程里用这个谱最多。JONSWAP 谱适合风区受限的海域比如近海、港湾它在 P-M 基础上加了谱峰增强因子 γ谱形更尖算出来的横摇响应往往更大更保守。初次跑通程序建议就用 ITTC 双参数谱参数直观、不依赖风向等把链路跑顺了再换 JONSWAP 或方向谱。3. 用MATLAB算横摇纵摇有义值核心代码与参数拆解3.1 主程序框架波谱离散化、RAO、谱矩一脉算下来下面这段是完整的可运行主程序读入海况参数算横摇、纵摇响应谱和统计量最后出图。代码按“参数 → 谱 → RAO → 响应谱 → 统计量 → 绘图”组织复制保存为 ship_roll_pitch_spectrum.m 直接运行。%% ship_roll_pitch_spectrum.m % 基于平稳随机过程理论由海浪谱和船舶RAO计算横摇、纵摇有义单幅值 clear; clc; close all; %% 1. 海况参数ITTC双参数谱 Hs 3.0; % 有义波高m T1 8.5; % 特征周期sITTC谱参数不是跨零周期 w_min 0.05; % 频率下限 rad/s必须覆盖长波低频段 w_max 3.0; % 频率上限 rad/s覆盖高频段即可 Nw 600; % 频率离散点数 %% 2. 频率网格列向量 w linspace(w_min, w_max, Nw); %% 3. 波谱密度 S_wave S_wave Hs^2/(4*pi) * (2*pi/T1).^4 ./ w.^5 ... .* exp( -1/pi * (2*pi./(T1*w)).^4 ); %% 4. 横摇 RAO二阶振荡系统模型 wn_r 0.55; % 横摇固有频率 rad/s zeta_r 0.10; % 横摇阻尼比 k_r 0.60; % 低频增益按船型/装载修正 H_roll k_r ./ sqrt( (1 - (w/wn_r).^2).^2 ... (2*zeta_r*w/wn_r).^2 ); %% 5. 纵摇 RAO同一二阶模型参数按船型换 wn_p 0.85; % 纵摇固有频率 rad/s zeta_p 0.30; % 纵摇阻尼比通常明显大于横摇 k_p 1.00; % 纵摇低频增益接近1 H_pitch k_p ./ sqrt( (1 - (w/wn_p).^2).^2 ... (2*zeta_p*w/wn_p).^2 ); %% 6. 响应谱波谱乘RAO模平方 S_roll S_wave .* H_roll.^2; S_pitch S_wave .* H_pitch.^2; %% 7. 谱矩与有义单幅值 m0_roll trapz(w, S_roll); m2_roll trapz(w, S_roll .* w.^2); m0_pitch trapz(w, S_pitch); m2_pitch trapz(w, S_pitch .* w.^2); % 运动量习惯报有义单幅等于 2*sqrt(m0) A1_3_roll 2 * sqrt(m0_roll); A1_3_pitch 2 * sqrt(m0_pitch); % 平均跨零周期 Tz_roll 2*pi * sqrt(m0_roll/m2_roll); Tz_pitch 2*pi * sqrt(m0_pitch/m2_pitch); fprintf(横摇有义单幅 %.2f\n, A1_3_roll); fprintf(纵摇有义单幅 %.2f\n, A1_3_pitch); fprintf(横摇平均跨零周期 %.2f s\n, Tz_roll); fprintf(纵摇平均跨零周期 %.2f s\n, Tz_pitch); %% 8. 绘图左波浪谱右响应谱 figure(Color,w); subplot(2,1,1); plot(w, S_wave, LineWidth, 1.5); xlabel(\omega (rad/s)); ylabel(S_wave (m^2.s/rad)); title(ITTC双参数海浪谱); grid on; subplot(2,1,2); plot(w, S_roll, LineWidth, 1.5); hold on; plot(w, S_pitch, LineWidth, 1.5); xlabel(\omega (rad/s)); ylabel(S_response); legend(横摇响应谱,纵摇响应谱, Location,best); grid on;逻辑说明程序先建立频率网格再算波谱和两个 RAO然后用点乘计算响应谱。这里全部用列向量和点运算避免 for 循环和行列维数错位。谱矩用 trapz 做梯形积分m₀ 是曲线下面积m₂ 是对 ω² 加权后的面积两个都来自同一根 S_response 曲线。有义单幅和跨零周期由谱矩一次算出不需要合成时历。参数说明w_min 取 0.05 rad/s对应周期约 125 s覆盖涌浪低频段如果只关注风浪取 0.1 可以但结果会丢掉一部分低频面积。Nw 取 600 是最低配置我一般用 10001500尤其波谱峰和 RAO 峰接近时要加密。Hs 和 T1 来自海况统计资料T1 不是普遍意义上有义波高对应周期的 T_z而是 ITTC 谱自身的输入参数二者数值相近但不要混。3.2 三个必调参数频率范围、离散点数、谱参数第一个必调参数是频率范围。横摇固有频率典型在 0.30.7 rad/s波谱峰在 0.41.0 rad/s纵摇固有频率更高一些覆盖 0.053.0 rad/s 基本够用。最高频率要不要再往上提取决于 RAO 数据来源如果 RAO 是势流软件算的通常会给到 2.53.0 rad/s如果只有模型试验数据高频段数据噪声大建议直接截断到试验覆盖范围不要外推。第二个必调参数是阻尼比 zeta 和固有频率 wn。它们决定 RAO 峰的位置和高度是整个计算里最“玄学”的部分。横摇阻尼比小峰高而窄纵摇阻尼比大峰低而宽。实船设计阶段没有衰减试验数据时横摇取 0.060.12纵摇取 0.20.35算出来的结果属于“设计谱平均值”还要靠后续试验修正。这里不值得花时间把二阶模型调成“完美贴合某条试验曲线”第一次跑通更重要。第三个必调参数是 Hs 和 T1。海况资料里如果有实测谱可以直接替代经验谱如果没有用 ITTC 双参数谱把 Hs 和 T1 从海况预报里取出来就行。需要注意T1 是用谱零阶矩和一阶矩定义的与目测周期、跨零周期都有系统差工程预报单里有时给的是峰值周期要先换算回去再代入。换算关系我没法给一个全网通用的常数最可靠的做法是拿到海况统计的原数据来源确认它定义的周期是哪一种。3.3 把纵摇接进来换RAO曲线代码怎么组织纵摇和横摇在程序上区别很小RAO 形状不同而已纵摇固有频率一般比横摇高阻尼比大低频增益接近 1响应谱更宽。代码里我用同一个二阶系统表达式换了 wn、zeta、k 三个参数。实际项目中纵摇 RAO 往往来自水动力软件是一组频率和幅值对应的离散表这时不要把这个表拟合进二阶系统直接插值到统一网格% 水动力软件给出的 RAO 数据表示例格式 % w_rao_pitch 为频率列rao_pitch 为幅值列 wq w; % 全局统一频率网格 H_pitch_data interp1(w_rao_pitch, rao_pitch, wq, pchip, 0);使用 pchip 而不是 linear 作为插值方法是为了 RAO 峰附近不至于被线性折线削平。插值这一步很关键后面波谱网格和 RAO 网格必须完全一致才能点乘不一致时要么先插值统一要么会用错位数据算出毫无意义的谱。我见过有人直接把两个不同频率分辨率的数组点乘MATLAB 报矩阵维度错误后用循环硬算结果响应谱比波谱还大一个量级最后发现是频率网格错位。这个坑后面避坑章节还会说。4. 从谱逆推时历与短时预报FFT点数、窗函数和采样率怎么选4.1 为何要合成时历谱只给统计给不了最大单幅谱分析给出的有义单幅、跨零周期是统计期望值实际问题往往问的是“十分钟内最大横摇多少度”“这个海况能不能安全工作”。回答这类问题需要时历。在平稳随机过程框架里时历可以按“随机相位叠加法”合成把连续的功率谱离散成有限个频带每个频带用幅值 √(2S(ω)Δω) 的余弦波相位取随机均匀分布叠加后就是一个功率谱与目标谱一致的平稳高斯过程。合成时历的价值有三层一是把谱分析结果变成直观的横摇角、纵摇角曲线便于和试航数据对照二是做峰值统计回答“最大多少度”三是后续接时域仿真、控制系统设计时这个时历作为输入。合成时历本质是和谱分析同一套物理模型不要指望它给出超越谱统计的新信息它只是把同一个模型从频域换到时域呈现。4.2 合成时历的参数持续时间、采样率、随机相位下面的函数接收谱密度和频率网格返回一条合成的横摇角时历。也可以直接喂波谱生成波浪时历函数的输入 S 换成其它谱密度即可。function [t, y] synth_rand_process(w, S, duration, dt, seed) % 按给定功率谱合成平稳随机时历 % 输入 % w : 频率网格列向量rad/s % S : 与 w 对应的谱密度 % duration : 时历长度s % dt : 采样间隔s % seed : 随机种子可复现 % 输出 % t : 时间列向量s % y : 合成时历单位与 S 的开方一致 if nargin 5 seed 1; end rng(seed); w w(:); S S(:); dw mean(diff(w)); amp sqrt(2 * S * dw); % 每频带的余弦分量幅值 phase 2 * pi * rand(size(S)); % 均匀随机相位 [0, 2pi) n round(duration / dt); t (0:n-1) * dt; y cos(w * t phase) * amp; % NxM 加权累加得到 N 维时历 y y(:); end逻辑说明余弦分量幅值来自能量谱密度乘以频带宽度再开根号这是严格满足“谱密度守恒”的做法。相位用 rand 均匀分布在 0 到 2π每个频率独立叠加结果自动逼近高斯分布频率点数越多越明显600 个频率足以让时历统计接近理论高斯。返回值 y 的物理单位取决于 SS 是横摇响应谱y 就是横摇角S 是波谱y 就是波面高度。参数说明duration 至少要覆盖最低频率一个完整周期w_min0.05 rad/s 对应 126 s 周期所以最短也要 300 s工程上建议 6001200 s也就是 1020 分钟顺着短期海况的平稳窗口走。dt 由最高频率决定w_max3.0 rad/s 对应 0.477 Hz根据 Nyquist 要求采样频率至少两倍、实际取 510 倍dt 取 0.10.2 s 完全足够。不要为了省内存把 dt 拉到 1 s那样最高频分量会混叠成低频虚假能量。4.3 FFT谱估计与窗函数为什么合成谱和理论谱差一截把合成时历拿去做 FFT再用周期图法估计功率谱会发现估计谱在峰附近有起伏、在低谷有泄漏尤其低频段比理论谱高出一块。这不是程序 bug是有限长时历的谱估计方差和窗函数泄漏。遇到这类问题我一般不用手写 FFT 估计直接用 Signal Processing Toolbox 的 pwelch 加 Hann 窗窗长取时历的一半重叠 50%Fs 1 / dt; % 采样频率 Hz [pxx, f] pwelch(y, hann(round(length(y)/2)), [], [], Fs); omega_est 2 * pi * f; S_est pxx * pi / (2*pi); % 单边谱转换为圆频率谱密度这段里 pwelch 的结果是线频率单边功率谱要除以 2π 转成圆频率谱密度才能和理论谱比。如果没装工具箱用 matlab 自带 periodogram 函数也可以Hann 窗照用。对比时重点看谱峰频率和峰下面积不要看单点值单点波动是正常的。如果不加窗直接对整条时历做 FFT泄漏会让低频段出现明显的抬升这是对窄带信号做短样本谱估计的典型现象。判断估计是否合格不要求曲线重合要求 m₀ 在 5% 以内、峰值频率不偏。这个指标后面第 6 章的验证试验还会再提。5. 船舶横摇纵摇谱分析避坑指南五处容易翻车的地方5.1 中文注释乱码GBK与UTF-8的编码冲突现象从网上下载或同事发来的 m 文件中文注释在编辑器里显示成乱码常见的是“涓€浜涘睆骞”之类的错乱字符程序本身能跑但没法读。原因文件保存编码是 UTF-8而 matlab 旧版本或某些系统区域设置下默认按 GBK 解码或者反过来文件是 GBK 保存而 matlab 按 UTF-8 读取。这个问题在 R2023b 之后仍然存在只是报错形式不同而已不是新装的软件坏了。解决在 matlab 2023b 及以后版本打开文件编辑器里 CtrlA 全选后刷新没用正确做法是“Preferences → Editor/Debugger → Language”把文本编码设为 UTF-8重启编辑器再打开如果文件本来就是 GBK要把编码改为 GBK 或 GB18030。团队协作时建议统一全部文件保存为 UTF-8并在每个文件头部加英文注释说明海况参数单位避免编码问题干扰数值理解。5.2 RAO峰与波谱峰重叠有义值虚高阻尼比是关键现象算出来的横摇有义单幅 8 度经验公式估算只有 4 度差距大得离谱。检查 RAO 固有频率 0.55 rad/s波谱谱峰也恰好落在 0.55 rad/s 附近两个峰重叠响应谱出现一个又高又窄的尖峰。原因平稳随机过程线性理论里响应谱是波谱乘 RAO 模方两个峰重叠时能量被双重放大。这本身没有错但如果 RAO 阻尼比取太小模型谐振峰过于尖锐实际船的非线性阻尼在大幅横摇下会压住峰值线性模型就虚高了。另一个常见原因是横摇 RAO 峰附近频率步长太粗梯形积分把这个窄峰“跳过”了一部分或者恰好跨过峰值。解决先把 zeta 提高到试验或规范建议范围再用 1200 个频率点重新算响应谱峰附近可以单独加密比如在 0.40.7 rad/s 范围内局部加网格再重新 trapezoid 积分。如果虚高仍然存在要怀疑 RAO 本身势流软件无黏性计算对横摇阻尼估计偏低横摇 RAO 峰值甚至会达到 58这时必须用试验或经验阻尼修正。5.3 频率下限设太高低频长波能量被白白丢掉现象同一艘船、同一个海况把 w_min 从 0.05 改成 0.15有义单幅值掉了 8% 以上纵摇变化更明显。原因波谱低频段虽然幅值小但 ω⁻⁵ 使谱随频率降低快速增长0.10.2 rad/s 区间并不小加上横摇 RAO 低频段增益不为零长周期能量会被横摇响应接收。频率下限人为截断造成的误差不是 1%、2%是系统性的。解决频率下限尽量取 0.030.05 rad/s。浅水或受限水域长波传播受地形影响波谱本身在低频会有衰减这时可以合理截断 0.1但要先确认本海区的实测谱低频段趋势不能看到跨零周期对上了就认为频率范围没问题。判断截断误差的方法很简单把 w_min 减半如果结果变化小于 1%说明下限足够否则继续降低。5.4 横摇纵摇RAO频率网格不一致插值对齐是必修课现象横摇 RAO 来自某个软件的 0.13.0 rad/s 共 120 点纵摇 RAO 来自另一个报告的 0.23.5 rad/s 共 95 点直接把两个 RAO 与波谱点乘纵摇那条报“矩阵维度不一致”或者算出来纵摇响应谱形状怪异高频出现锯齿。原因频率网格不一致点乘无法进行用循环硬算时频率错位导致某个频率的波谱能量乘了错误频率的 RAO 值谱峰位置偏移统计量自然不可信。解决统一用目标网格对每一条 RAO 数据用 interp1 插值到全局网格方法选 pchip比例 0 作为超出范围的补值。代码在第 3.3 节已给出。插值后画一条重叠图确认所有曲线的频率轴一致再进入响应谱计算。这个“先对齐频率轴再点乘”的习惯能省掉一大类排查时间。5.5 只报有义值不够跨零周期、谱宽、极值都要给现象报告里只有“横摇有义单幅 3.2 度”这一句话评审问“最大横摇多少”“平均周期多少”“谱峰周期多少”一个都答不上。原因有义值只是谱零阶矩的表现它丢失了谱形信息。两条船一条谱窄一条谱宽有义值相同可最大横摇预期和驾驶体验完全不同跨零周期决定船舶在某一海况下的“摇动节奏”直接影响晕船评价、设备疲劳这些信息都在谱矩里不取出来就白算了。解决至少补三样东西。跨零周期 T_z 2π√(m₀/m₂)谱峰周期 T_p 直接在响应谱里找峰值频率反算最大预期单幅用 Rayleigh 分布估算对 N 个波峰最大单幅约等于有义单幅乘以 √(0.5·ln N) 的修正更稳的做法是合成 600 s 时历直接数极值那个结果更直观也更好向非专业的人解释。6. 用模拟信号验证谱分析程序一小时校核三件事6.1 用高密度积分做基准检验粗网格误差程序写完先别急着换参数先用同一组海况把“参考解”算出来频率网格取 0.024.0 rad/s、2400 点trapz 积分作为基准。然后把主程序频率网格 600 点、范围 0.053.0 的结果跟基准比横摇有义单幅误差应小于 2%。这一步目的是确认离散误差和截断误差总体可接受也顺手验证插值后的 RAO 完整。% 参考解高密度网格 w_ref linspace(0.02, 4.0, 2400); S_ref Hs^2/(4*pi) * (2*pi/T1).^4 ./ w_ref.^5 ... .* exp(-1/pi * (2*pi./(T1*w_ref)).^4); H_ref 0.6 ./ sqrt((1-(w_ref/0.55).^2).^2 (2*0.1*w_ref/0.55).^2); m0_ref trapz(w_ref, S_ref .* H_ref.^2); A_ref 2*sqrt(m0_ref); % 与主程序 A1_3_roll 对比误差超过2%则加密网格6.2 时域统计对谱统计十分钟时历的方差有多大合成一条横摇时历时长 600 sdt 取 0.2 s分别统计时历标准差和有义单幅再和谱矩值对比。这里注意时历有义值按“峰谷差的一半取前三分之一大值平均”算和谱矩 2√m₀ 是两种不同估计允许差 5% 左右。多换几个随机种子跑能看出有限时历的统计波动范围这个波动是真实的不代表程序有问题。我用过 matlab Online 跑这段代码纯基础函数不依赖额外箱子网页版也能出结果。6.3 形成自己的报数习惯希望帮到你最后养成一个习惯每次运行记录三行字——海况参数 Hs、T1频率范围与点数RAO来源与阻尼比。这个习惯救过我很多次隔一个月回看结果时不会把“用 P-M 谱算的”和“用 ITTC 谱算的”混在一起。凡是谱分析结果都主动标注“线性理论、短期海况假设、统计值非瞬时值”让下游用数据的人知道边界。项目里栽过的跟头都是这几类编码乱码、RAO频率没对齐、低频截断太狠、只报有义值不报周期。把这几个点记牢平稳随机过程这套方法在横摇、纵摇预报上能稳定输出可用的设计输入。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/1 1:04:00 · 阅读 1.2k 船舶横摇纵摇谱分析：平稳随机过程到MATLAB实现
简介一套基于平稳随机过程理论开展船舶海上横摇、纵摇运动仿真的MATLAB资料包面向船舶海洋工程、航海技术等专业本科生与硕士生适用于课程设计、毕业设计及教研验证场景。包内共8个文件5个Simulink模型.mdl用于构建船舶运动仿真子系统与网络化仿真模块2个MATLAB脚本.m承担数据处理与仿真流程驱动另附1个txt说明文档便于快速上手整体压缩包仅27KBMATLAB 2019a环境即可运行。目前已有529人学习对初学者而言可对照随机海浪作用下的横摇、纵摇响应分析练习理解平稳随机过程理论如何落实到仿真实现。通过梳理各模型与脚本的调用关系还能掌握运动仿真建模思路、参数设置与结果分析流程为后续船舶耐波性与波浪载荷研究奠定基础。1. 平稳随机过程算船舶横摇纵摇到底在算什么拿到一个耐波性任务第一步通常不是仿真而是先回答一个很实际的问题这条船在 5 级海况下横摇有义单幅是多少度、纵摇多少度、多长时间摇一个来回。海上实际波浪是不规则波一条条规则波去叠既不现实也不可靠工程上把这问题归到平稳随机过程框架里把海浪看作平稳随机过程用谱分析一次算出横摇、纵摇的统计特征值。这个方案能给出有义单幅、平均跨零周期等设计输入但给不了“下一分钟具体摇到几度”这种瞬时值做不了也要知道边界在哪。适合船型开发、舾装设计、海事软件配套的工程师也适合刚入门想做耐波性基线的学生。2. 把海浪当成平稳随机过程从波谱、RAO到响应谱的三段式2.1 为什么“平稳随机”够用波面位移与各态历经性海上波浪在严格意义上是非平稳的风在变、潮在变、船在动。但工程处理有一个“短期海况”假设在 20 分钟到 3 小时这个时间窗口内认为波高统计特征不随时间平移改变均值、方差、相关函数都稳定这就是平稳性再把大量不同频率、随机相位的余弦波叠加成波面得到的就是平稳随机过程的一个现实样本。真正让谱分析能在 MATLAB 里落地的是各态历经性。意思是我不用对同一海况做几百次海洋试验取集合平均只要取一条足够长的时历做时间平均就能近似得到集合平均的结果。这就是用一段 20 分钟模拟时历验证谱矩能成立的根据。窄带随机过程概念在这里很关键横摇、纵摇响应谱往往集中在一个窄频带里时历看起来像“载波频率不变、幅值缓慢变化”的调幅波峰值包络服从 Rayleigh 分布。有义单幅值、最大预期值都是从这条假设延伸出去的。这套框架成立有两个前提一是海况参数在这个窗口里稳定二是船的运动响应是线性的。大角度横摇超过 15~20 度时非线性阻尼会改变谱形纯谱分析就不够准这点后面避坑章节还会提到。2.2 运动响应谱的三段式波谱乘RAO模方船舶对不规则波的响应在平稳随机过程理论里不是直接“波高乘一个系数”而是三个环节串起来海浪谱 S_wave(ω)描述波能按圆频率 ω 的分布单位是 m²·s/rad。RAO幅值响应算子H(iω)船在单位波幅规则波下的运动幅值响应取决于船型、航速、浪向。响应谱 S_response(ω)横摇、纵摇运动的能量谱等于波谱乘 RAO 模的平方。写成式子就是最常用的三段式S_response(ω) S_wave(ω) · |H(iω)|²这里有一个新手常忽略的细节|H(iω)|² 是幅值平方不是 H(iω) 直接乘。RAO 通常由势流软件如 WAMIT、SESAM、AQWA算出来或者由模型试验得到横摇、纵摇都有自己的实部和虚部或幅值、相位。如果你拿到的是复数形式的 RAO做谱分析时只需要幅值部分相位角在谱分析里被丢掉了这不是失误因为功率谱本身只保留幅值信息只有合成时历做短时预报时才需要相位。响应谱算出来后没有直接“读”出运动幅值。要从谱矩入手m_n ∫₀^∞ ωⁿ · S_response(ω) dω其中常用的是零阶矩 m₀ 和二阶矩 m₂。有义单幅值横摇角、纵摇角习惯用单幅就等于 2√m₀平均跨零周期 T_z 2π√(m₀/m₂)。注意波浪的有义波高是 4√m₀因为波高是峰谷差、是双幅运动响应大家习惯报单幅所以用 2√m₀。这个单位约定如果混了结果会差整整一倍。2.3 波谱选型P-M谱、ITTC双参数谱怎么给参数海浪谱的选取直接决定响应谱的绝对量级是整条链路的起点。工程上最常用的两类Pierson-Moskowitz 谱P-M 谱用于充分成长的风浪输入是风速。它以风速 U₁₉.₅ 为参数谱形简单适合开阔海域。但它只有一个输入没法同时控制波高和周期遇到海况预报给了“有义波高 周期”两个量时直接用 P-M 谱会有偏差。ITTC 双参数谱是 P-M 谱的工程改写输入换成更容易得到的两个量有义波高 Hs 和特征周期 T₁。谱形如下S_wave(ω) (Hs² / 4π) · (2π/T₁)⁴ · ω⁻⁵ · exp( −(1/π) · (2π/(T₁·ω))⁴ )这个谱的好处是 Hs 和 T₁ 解耦海况预报、统计资料基本都按这两个量给我实际工程里用这个谱最多。JONSWAP 谱适合风区受限的海域比如近海、港湾它在 P-M 基础上加了谱峰增强因子 γ谱形更尖算出来的横摇响应往往更大更保守。初次跑通程序建议就用 ITTC 双参数谱参数直观、不依赖风向等把链路跑顺了再换 JONSWAP 或方向谱。3. 用MATLAB算横摇纵摇有义值核心代码与参数拆解3.1 主程序框架波谱离散化、RAO、谱矩一脉算下来下面这段是完整的可运行主程序读入海况参数算横摇、纵摇响应谱和统计量最后出图。代码按“参数 → 谱 → RAO → 响应谱 → 统计量 → 绘图”组织复制保存为 ship_roll_pitch_spectrum.m 直接运行。%% ship_roll_pitch_spectrum.m % 基于平稳随机过程理论由海浪谱和船舶RAO计算横摇、纵摇有义单幅值 clear; clc; close all; %% 1. 海况参数ITTC双参数谱 Hs 3.0; % 有义波高m T1 8.5; % 特征周期sITTC谱参数不是跨零周期 w_min 0.05; % 频率下限 rad/s必须覆盖长波低频段 w_max 3.0; % 频率上限 rad/s覆盖高频段即可 Nw 600; % 频率离散点数 %% 2. 频率网格列向量 w linspace(w_min, w_max, Nw); %% 3. 波谱密度 S_wave S_wave Hs^2/(4*pi) * (2*pi/T1).^4 ./ w.^5 ... .* exp( -1/pi * (2*pi./(T1*w)).^4 ); %% 4. 横摇 RAO二阶振荡系统模型 wn_r 0.55; % 横摇固有频率 rad/s zeta_r 0.10; % 横摇阻尼比 k_r 0.60; % 低频增益按船型/装载修正 H_roll k_r ./ sqrt( (1 - (w/wn_r).^2).^2 ... (2*zeta_r*w/wn_r).^2 ); %% 5. 纵摇 RAO同一二阶模型参数按船型换 wn_p 0.85; % 纵摇固有频率 rad/s zeta_p 0.30; % 纵摇阻尼比通常明显大于横摇 k_p 1.00; % 纵摇低频增益接近1 H_pitch k_p ./ sqrt( (1 - (w/wn_p).^2).^2 ... (2*zeta_p*w/wn_p).^2 ); %% 6. 响应谱波谱乘RAO模平方 S_roll S_wave .* H_roll.^2; S_pitch S_wave .* H_pitch.^2; %% 7. 谱矩与有义单幅值 m0_roll trapz(w, S_roll); m2_roll trapz(w, S_roll .* w.^2); m0_pitch trapz(w, S_pitch); m2_pitch trapz(w, S_pitch .* w.^2); % 运动量习惯报有义单幅等于 2*sqrt(m0) A1_3_roll 2 * sqrt(m0_roll); A1_3_pitch 2 * sqrt(m0_pitch); % 平均跨零周期 Tz_roll 2*pi * sqrt(m0_roll/m2_roll); Tz_pitch 2*pi * sqrt(m0_pitch/m2_pitch); fprintf(横摇有义单幅 %.2f\n, A1_3_roll); fprintf(纵摇有义单幅 %.2f\n, A1_3_pitch); fprintf(横摇平均跨零周期 %.2f s\n, Tz_roll); fprintf(纵摇平均跨零周期 %.2f s\n, Tz_pitch); %% 8. 绘图左波浪谱右响应谱 figure(Color,w); subplot(2,1,1); plot(w, S_wave, LineWidth, 1.5); xlabel(\omega (rad/s)); ylabel(S_wave (m^2.s/rad)); title(ITTC双参数海浪谱); grid on; subplot(2,1,2); plot(w, S_roll, LineWidth, 1.5); hold on; plot(w, S_pitch, LineWidth, 1.5); xlabel(\omega (rad/s)); ylabel(S_response); legend(横摇响应谱,纵摇响应谱, Location,best); grid on;逻辑说明程序先建立频率网格再算波谱和两个 RAO然后用点乘计算响应谱。这里全部用列向量和点运算避免 for 循环和行列维数错位。谱矩用 trapz 做梯形积分m₀ 是曲线下面积m₂ 是对 ω² 加权后的面积两个都来自同一根 S_response 曲线。有义单幅和跨零周期由谱矩一次算出不需要合成时历。参数说明w_min 取 0.05 rad/s对应周期约 125 s覆盖涌浪低频段如果只关注风浪取 0.1 可以但结果会丢掉一部分低频面积。Nw 取 600 是最低配置我一般用 10001500尤其波谱峰和 RAO 峰接近时要加密。Hs 和 T1 来自海况统计资料T1 不是普遍意义上有义波高对应周期的 T_z而是 ITTC 谱自身的输入参数二者数值相近但不要混。3.2 三个必调参数频率范围、离散点数、谱参数第一个必调参数是频率范围。横摇固有频率典型在 0.30.7 rad/s波谱峰在 0.41.0 rad/s纵摇固有频率更高一些覆盖 0.053.0 rad/s 基本够用。最高频率要不要再往上提取决于 RAO 数据来源如果 RAO 是势流软件算的通常会给到 2.53.0 rad/s如果只有模型试验数据高频段数据噪声大建议直接截断到试验覆盖范围不要外推。第二个必调参数是阻尼比 zeta 和固有频率 wn。它们决定 RAO 峰的位置和高度是整个计算里最“玄学”的部分。横摇阻尼比小峰高而窄纵摇阻尼比大峰低而宽。实船设计阶段没有衰减试验数据时横摇取 0.060.12纵摇取 0.20.35算出来的结果属于“设计谱平均值”还要靠后续试验修正。这里不值得花时间把二阶模型调成“完美贴合某条试验曲线”第一次跑通更重要。第三个必调参数是 Hs 和 T1。海况资料里如果有实测谱可以直接替代经验谱如果没有用 ITTC 双参数谱把 Hs 和 T1 从海况预报里取出来就行。需要注意T1 是用谱零阶矩和一阶矩定义的与目测周期、跨零周期都有系统差工程预报单里有时给的是峰值周期要先换算回去再代入。换算关系我没法给一个全网通用的常数最可靠的做法是拿到海况统计的原数据来源确认它定义的周期是哪一种。3.3 把纵摇接进来换RAO曲线代码怎么组织纵摇和横摇在程序上区别很小RAO 形状不同而已纵摇固有频率一般比横摇高阻尼比大低频增益接近 1响应谱更宽。代码里我用同一个二阶系统表达式换了 wn、zeta、k 三个参数。实际项目中纵摇 RAO 往往来自水动力软件是一组频率和幅值对应的离散表这时不要把这个表拟合进二阶系统直接插值到统一网格% 水动力软件给出的 RAO 数据表示例格式 % w_rao_pitch 为频率列rao_pitch 为幅值列 wq w; % 全局统一频率网格 H_pitch_data interp1(w_rao_pitch, rao_pitch, wq, pchip, 0);使用 pchip 而不是 linear 作为插值方法是为了 RAO 峰附近不至于被线性折线削平。插值这一步很关键后面波谱网格和 RAO 网格必须完全一致才能点乘不一致时要么先插值统一要么会用错位数据算出毫无意义的谱。我见过有人直接把两个不同频率分辨率的数组点乘MATLAB 报矩阵维度错误后用循环硬算结果响应谱比波谱还大一个量级最后发现是频率网格错位。这个坑后面避坑章节还会说。4. 从谱逆推时历与短时预报FFT点数、窗函数和采样率怎么选4.1 为何要合成时历谱只给统计给不了最大单幅谱分析给出的有义单幅、跨零周期是统计期望值实际问题往往问的是“十分钟内最大横摇多少度”“这个海况能不能安全工作”。回答这类问题需要时历。在平稳随机过程框架里时历可以按“随机相位叠加法”合成把连续的功率谱离散成有限个频带每个频带用幅值 √(2S(ω)Δω) 的余弦波相位取随机均匀分布叠加后就是一个功率谱与目标谱一致的平稳高斯过程。合成时历的价值有三层一是把谱分析结果变成直观的横摇角、纵摇角曲线便于和试航数据对照二是做峰值统计回答“最大多少度”三是后续接时域仿真、控制系统设计时这个时历作为输入。合成时历本质是和谱分析同一套物理模型不要指望它给出超越谱统计的新信息它只是把同一个模型从频域换到时域呈现。4.2 合成时历的参数持续时间、采样率、随机相位下面的函数接收谱密度和频率网格返回一条合成的横摇角时历。也可以直接喂波谱生成波浪时历函数的输入 S 换成其它谱密度即可。function [t, y] synth_rand_process(w, S, duration, dt, seed) % 按给定功率谱合成平稳随机时历 % 输入 % w : 频率网格列向量rad/s % S : 与 w 对应的谱密度 % duration : 时历长度s % dt : 采样间隔s % seed : 随机种子可复现 % 输出 % t : 时间列向量s % y : 合成时历单位与 S 的开方一致 if nargin 5 seed 1; end rng(seed); w w(:); S S(:); dw mean(diff(w)); amp sqrt(2 * S * dw); % 每频带的余弦分量幅值 phase 2 * pi * rand(size(S)); % 均匀随机相位 [0, 2pi) n round(duration / dt); t (0:n-1) * dt; y cos(w * t phase) * amp; % NxM 加权累加得到 N 维时历 y y(:); end逻辑说明余弦分量幅值来自能量谱密度乘以频带宽度再开根号这是严格满足“谱密度守恒”的做法。相位用 rand 均匀分布在 0 到 2π每个频率独立叠加结果自动逼近高斯分布频率点数越多越明显600 个频率足以让时历统计接近理论高斯。返回值 y 的物理单位取决于 SS 是横摇响应谱y 就是横摇角S 是波谱y 就是波面高度。参数说明duration 至少要覆盖最低频率一个完整周期w_min0.05 rad/s 对应 126 s 周期所以最短也要 300 s工程上建议 6001200 s也就是 1020 分钟顺着短期海况的平稳窗口走。dt 由最高频率决定w_max3.0 rad/s 对应 0.477 Hz根据 Nyquist 要求采样频率至少两倍、实际取 510 倍dt 取 0.10.2 s 完全足够。不要为了省内存把 dt 拉到 1 s那样最高频分量会混叠成低频虚假能量。4.3 FFT谱估计与窗函数为什么合成谱和理论谱差一截把合成时历拿去做 FFT再用周期图法估计功率谱会发现估计谱在峰附近有起伏、在低谷有泄漏尤其低频段比理论谱高出一块。这不是程序 bug是有限长时历的谱估计方差和窗函数泄漏。遇到这类问题我一般不用手写 FFT 估计直接用 Signal Processing Toolbox 的 pwelch 加 Hann 窗窗长取时历的一半重叠 50%Fs 1 / dt; % 采样频率 Hz [pxx, f] pwelch(y, hann(round(length(y)/2)), [], [], Fs); omega_est 2 * pi * f; S_est pxx * pi / (2*pi); % 单边谱转换为圆频率谱密度这段里 pwelch 的结果是线频率单边功率谱要除以 2π 转成圆频率谱密度才能和理论谱比。如果没装工具箱用 matlab 自带 periodogram 函数也可以Hann 窗照用。对比时重点看谱峰频率和峰下面积不要看单点值单点波动是正常的。如果不加窗直接对整条时历做 FFT泄漏会让低频段出现明显的抬升这是对窄带信号做短样本谱估计的典型现象。判断估计是否合格不要求曲线重合要求 m₀ 在 5% 以内、峰值频率不偏。这个指标后面第 6 章的验证试验还会再提。5. 船舶横摇纵摇谱分析避坑指南五处容易翻车的地方5.1 中文注释乱码GBK与UTF-8的编码冲突现象从网上下载或同事发来的 m 文件中文注释在编辑器里显示成乱码常见的是“涓€浜涘睆骞”之类的错乱字符程序本身能跑但没法读。原因文件保存编码是 UTF-8而 matlab 旧版本或某些系统区域设置下默认按 GBK 解码或者反过来文件是 GBK 保存而 matlab 按 UTF-8 读取。这个问题在 R2023b 之后仍然存在只是报错形式不同而已不是新装的软件坏了。解决在 matlab 2023b 及以后版本打开文件编辑器里 CtrlA 全选后刷新没用正确做法是“Preferences → Editor/Debugger → Language”把文本编码设为 UTF-8重启编辑器再打开如果文件本来就是 GBK要把编码改为 GBK 或 GB18030。团队协作时建议统一全部文件保存为 UTF-8并在每个文件头部加英文注释说明海况参数单位避免编码问题干扰数值理解。5.2 RAO峰与波谱峰重叠有义值虚高阻尼比是关键现象算出来的横摇有义单幅 8 度经验公式估算只有 4 度差距大得离谱。检查 RAO 固有频率 0.55 rad/s波谱谱峰也恰好落在 0.55 rad/s 附近两个峰重叠响应谱出现一个又高又窄的尖峰。原因平稳随机过程线性理论里响应谱是波谱乘 RAO 模方两个峰重叠时能量被双重放大。这本身没有错但如果 RAO 阻尼比取太小模型谐振峰过于尖锐实际船的非线性阻尼在大幅横摇下会压住峰值线性模型就虚高了。另一个常见原因是横摇 RAO 峰附近频率步长太粗梯形积分把这个窄峰“跳过”了一部分或者恰好跨过峰值。解决先把 zeta 提高到试验或规范建议范围再用 1200 个频率点重新算响应谱峰附近可以单独加密比如在 0.40.7 rad/s 范围内局部加网格再重新 trapezoid 积分。如果虚高仍然存在要怀疑 RAO 本身势流软件无黏性计算对横摇阻尼估计偏低横摇 RAO 峰值甚至会达到 58这时必须用试验或经验阻尼修正。5.3 频率下限设太高低频长波能量被白白丢掉现象同一艘船、同一个海况把 w_min 从 0.05 改成 0.15有义单幅值掉了 8% 以上纵摇变化更明显。原因波谱低频段虽然幅值小但 ω⁻⁵ 使谱随频率降低快速增长0.10.2 rad/s 区间并不小加上横摇 RAO 低频段增益不为零长周期能量会被横摇响应接收。频率下限人为截断造成的误差不是 1%、2%是系统性的。解决频率下限尽量取 0.030.05 rad/s。浅水或受限水域长波传播受地形影响波谱本身在低频会有衰减这时可以合理截断 0.1但要先确认本海区的实测谱低频段趋势不能看到跨零周期对上了就认为频率范围没问题。判断截断误差的方法很简单把 w_min 减半如果结果变化小于 1%说明下限足够否则继续降低。5.4 横摇纵摇RAO频率网格不一致插值对齐是必修课现象横摇 RAO 来自某个软件的 0.13.0 rad/s 共 120 点纵摇 RAO 来自另一个报告的 0.23.5 rad/s 共 95 点直接把两个 RAO 与波谱点乘纵摇那条报“矩阵维度不一致”或者算出来纵摇响应谱形状怪异高频出现锯齿。原因频率网格不一致点乘无法进行用循环硬算时频率错位导致某个频率的波谱能量乘了错误频率的 RAO 值谱峰位置偏移统计量自然不可信。解决统一用目标网格对每一条 RAO 数据用 interp1 插值到全局网格方法选 pchip比例 0 作为超出范围的补值。代码在第 3.3 节已给出。插值后画一条重叠图确认所有曲线的频率轴一致再进入响应谱计算。这个“先对齐频率轴再点乘”的习惯能省掉一大类排查时间。5.5 只报有义值不够跨零周期、谱宽、极值都要给现象报告里只有“横摇有义单幅 3.2 度”这一句话评审问“最大横摇多少”“平均周期多少”“谱峰周期多少”一个都答不上。原因有义值只是谱零阶矩的表现它丢失了谱形信息。两条船一条谱窄一条谱宽有义值相同可最大横摇预期和驾驶体验完全不同跨零周期决定船舶在某一海况下的“摇动节奏”直接影响晕船评价、设备疲劳这些信息都在谱矩里不取出来就白算了。解决至少补三样东西。跨零周期 T_z 2π√(m₀/m₂)谱峰周期 T_p 直接在响应谱里找峰值频率反算最大预期单幅用 Rayleigh 分布估算对 N 个波峰最大单幅约等于有义单幅乘以 √(0.5·ln N) 的修正更稳的做法是合成 600 s 时历直接数极值那个结果更直观也更好向非专业的人解释。6. 用模拟信号验证谱分析程序一小时校核三件事6.1 用高密度积分做基准检验粗网格误差程序写完先别急着换参数先用同一组海况把“参考解”算出来频率网格取 0.024.0 rad/s、2400 点trapz 积分作为基准。然后把主程序频率网格 600 点、范围 0.053.0 的结果跟基准比横摇有义单幅误差应小于 2%。这一步目的是确认离散误差和截断误差总体可接受也顺手验证插值后的 RAO 完整。% 参考解高密度网格 w_ref linspace(0.02, 4.0, 2400); S_ref Hs^2/(4*pi) * (2*pi/T1).^4 ./ w_ref.^5 ... .* exp(-1/pi * (2*pi./(T1*w_ref)).^4); H_ref 0.6 ./ sqrt((1-(w_ref/0.55).^2).^2 (2*0.1*w_ref/0.55).^2); m0_ref trapz(w_ref, S_ref .* H_ref.^2); A_ref 2*sqrt(m0_ref); % 与主程序 A1_3_roll 对比误差超过2%则加密网格6.2 时域统计对谱统计十分钟时历的方差有多大合成一条横摇时历时长 600 sdt 取 0.2 s分别统计时历标准差和有义单幅再和谱矩值对比。这里注意时历有义值按“峰谷差的一半取前三分之一大值平均”算和谱矩 2√m₀ 是两种不同估计允许差 5% 左右。多换几个随机种子跑能看出有限时历的统计波动范围这个波动是真实的不代表程序有问题。我用过 matlab Online 跑这段代码纯基础函数不依赖额外箱子网页版也能出结果。6.3 形成自己的报数习惯希望帮到你最后养成一个习惯每次运行记录三行字——海况参数 Hs、T1频率范围与点数RAO来源与阻尼比。这个习惯救过我很多次隔一个月回看结果时不会把“用 P-M 谱算的”和“用 ITTC 谱算的”混在一起。凡是谱分析结果都主动标注“线性理论、短期海况假设、统计值非瞬时值”让下游用数据的人知道边界。项目里栽过的跟头都是这几类编码乱码、RAO频率没对齐、低频截断太狠、只报有义值不报周期。把这几个点记牢平稳随机过程这套方法在横摇、纵摇预报上能稳定输出可用的设计输入。希望帮到你。本文还有配套的精品资源点击获取