波浪是影响船舶与海洋工程结构安全的重要环境条件^[1]. 在结构设计、评估过程中，长期风浪可离散为一系列满足各态历经性假设的短期平稳风浪. 短期风浪的特征参数可以通过波高周期的联合分布得到精确描述^[1-2]. 因此，针对短期风浪波高周期的联合分布开展研究，具有重要的科学意义和实用价值.

1975年Longuet-Higgins提出周期呈对称状态的短期风浪波高周期联合分布模型^[3]. 这种对称状态与实际海况不太相符^[4]. 该模型在1983年被提出者修正^[5] (以下简称L-H 1983模型)，修正后的联合概率密度函数不再关于周期对称，但由此推导出的波高分布不再是瑞利(Rayleigh)分布. 孙孚提出新的联合分布模型^[6] (以下简称孙孚模型)，由此式推导出的波高分布仍然是瑞利分布. Zheng等^[7]对L-H 1983模型和孙孚模型进行修正，修正后的2种模型 (以下分别简称L-H 1983-Zheng模型、孙孚-Zheng模型)在周期均值的计算上和实际海浪更相符. 此外，Cavanié等^[8-9]给出的CNEXO模型需要使用波浪谱的四阶矩.

上述模型是在满足各态历经性假设的平稳高斯随机过程理论的基础上，结合一定的理论近似推导出来的. 潘锦嫦等^[10-11]的研究表明，这些模型不能满足精确预报的需求. 因此，有必要发展新的、更加精确的波高周期联合分布模型. 基于相同的理论基础，本研究在广泛数值分析的基础上，采用条件概率方法提出新的短期风浪波高周期联合分布模型，以实测波浪谱、常用的2种理论波浪谱〔Pierson-Moscowitz (P-M)谱和布氏(Bretschneider)谱^[12]〕为靶谱，通过模拟方法得到波高周期的经验联合分布，验证所提模型的正确性以及有效性.

CNEXO模型^[8-9]给出以波浪谱0、2、4阶谱矩：m₀、m₂、m₄为参数的波高H周期T的联合概率密度函数 $ {f_{\rm{C}}}(H,T) $，表达式为

(1) $ \begin{split} {f_{\text{C}}}(H,T) &= \dfrac{{{H^2}}}{{8{m_4}\sqrt {2\text{π} {m_0}\left( {1 - {\alpha ^2}} \right)} }}\dfrac{{32{{\text{π}}^4}}}{{{T^5}}} \times \\ &\exp \left[ { - \dfrac{{\left( {{m_4}{T^4} + 16{{\text{π}}^4}{m_0}{T^4} - 8{{\text{π}}^2}{T^2}{m_2}} \right) {H^2}}}{{8\left( {{m_0}{m_4} - m_2^2} \right){T^4}}}} \right]. \\[-8pt] \end{split}$

式中： $\alpha $为带宽参数， $\alpha = {{{m_2}} \mathord{\left/ {\vphantom {{{m_2}} {\sqrt {{m_0}{m_4}} }}} \right. } {\sqrt {{m_0}{m_4}} }}$。

L-H 1983模型^[5]给出以波浪谱0、1、2阶谱矩：m₀，m₁，m₂为参数的波高周期的联合概率密度函数 $ {f_{{\text{L - H}}}}(H,T) $，表达式为

(2) $ \begin{split} {f_{{\text{L - H}}}}(H,T) =& L\left( \nu \right)\frac{{\sqrt {\text{π }} }}{{4\nu \sqrt {2{m_0}} {m_1}}}\frac{{{H^2}}}{{{T^2}}} \times \\ &\exp \left\{ { - \frac{{{H^2}}}{{8{m_0}}}\left[ {1+\frac{1}{{{\nu ^2}}}{{\left( {\frac{{2{\text{π}}{m_0}}}{{T{m_1}}} - 1} \right)}^2}} \right]} \right\}, \\ & L\left( \nu \right) = 2{\left[ {1+{{\left( {1+{\nu ^2}} \right)}^{ - \tfrac{1}{2}}}} \right]^{ - 1}}. \\ \end{split} $

式中：v为Longuet-Higgins定义的谱宽参数^[3]， $\nu = \sqrt {{{{m_0}{m_2}} \mathord{\left/ {\vphantom {{{m_0}{m_2}} {m_1^2}}} \right. } {m_1^2}} - 1}$。

孙孚模型^[6]给出波高周期的联合概率密度函数 $ {f_{\text{S}}}(H,T) $，表达式为

(3) $ \begin{split} {f_{\text{S}}}(H,T) = & \frac{{\sqrt {\text{π }} }}{{4\nu \sqrt {2{m_0}} {m_1}}}\frac{{{H^2}}}{{{T^2}}}\left[ {1+\exp \left( { - \frac{{{\text{π}}{H^2}}}{{{m_1}{\nu ^2}T}}} \right)} \right] \times \\ & \exp \left\{ { - \frac{{{H^2}}}{{8{m_0}}}\left[ {1+\frac{1}{{{\nu ^2}}}{{\left( {1 - \frac{{2{\text{π }}{m_0}}}{{T{m_1}}}} \right)}^2}} \right]} \right\} . \\ \end{split} $

L-H 1983-Zheng模型^[7]给出波高周期的联合概率密度函数 $ {f_{{\text{LZ}}}}(H,T) $，表达式为

(4) $ \begin{split} {f_{{\text{LZ}}}}\left( {H,T} \right) =& \frac{{L\left( \nu \right)}}{{\sqrt {1+{\nu ^2}} }}\frac{{{\text{π }}\sqrt {2{\text{π}}{m_0}} }}{{4\nu m_1^2}}\frac{{{H^2}}}{{{T^3}}} \times \\ & \exp \left\{ { - \frac{{{H^2}}}{{8{m_0}}}\left[ {1+\frac{1}{{{\nu ^2}}}{{\left( {\frac{{2{\text{π}}{m_0}}}{{T{m_1}}} - 1} \right)}^2}} \right]} \right\} .\\ \end{split} $

孙孚-Zheng模型^[7]给出波高周期的联合概率密度函数 $ {f_{{\text{SZ}}}}(H,T) $，表达式为

(5) $ \begin{split} {f_{{\text{SZ}}}}\left( {H,T} \right) = & \frac{1}{{\sqrt {1+{\nu ^2}} }}\frac{{{\text{π}}\sqrt {2\text{π} {m_0}} }}{{4\nu m_1^2}}\frac{{{H^2}}}{{{T^3}}}\left[ {1+\exp \left( { - \frac{{{\text{π}}{H^2}}}{{{\nu ^2}{m_1}T}}} \right)} \right] \times \\ & \exp \left\{ { - \frac{{{H^2}}}{{8{m_0}}}\left[ {1+\frac{1}{{{\nu ^2}}}{{\left( {1 - \frac{{2{\text{π}}{m_0}}}{{T{m_1}}}} \right)}^2}} \right]} \right\} . \\[-7pt] \end{split} $

波高周期的联合概率密度函数可以表示为

(6) $ f(H,T) = f(H) \times f\left( {T\left| H \right.} \right) . $

式中： $ f(H) $为波高的概率密度函数， $f(T|H) $为已知波高时周期的条件概率密度函数. 波高H服从瑞利分布^[1,3,6,13].

(7) $ f(H) = \frac{H}{{4{m_0}}}\exp \left( { - \frac{{{H^2}}}{{8{m_0}}}} \right) . $

由式(7)可以推导得到波浪谱的零阶矩 $ {m_0} $与有义波高 $ {H_{{\rm{s}}}} $存在如下的关系^[1,4]：

(8) $ {m_0} = {{H_{{\rm{s}}}^{2}} \mathord{\left/ {\vphantom {{H_{s}^{2}} {16}}} \right. } {16}} . $

对于条件概率密度函数 $ f\left( {T\left| H \right.} \right) $，文献[13]建议使用高斯分布. 本研究通过大量文献、数据分析后发现， $ f\left( {T\left| H \right.} \right) $为具有一定偏态的单峰分布. 因此，可以通过构造新的 $ f\left( {T\left| H \right.} \right) $函数，得到新的波高周期联合分布模型. 基于大量数据分析，将 $ f\left( {T\left| H \right.} \right) $表示为以下一般形式：

(9) $ f\left( {T\left| H \right.} \right) = \frac{c}{T}{\left( {\ln T - b} \right)^{e - 1}}\exp \left[ { - {{\left( {\frac{{\ln T - b}}{a}} \right)}^d}} \right] . $

式中： $ a、b、c、d和e $均为与波高有关的参数，其中 $ d、e $可以分别固定为数值2和1. 根据概率密度函数积分为1的条件，进一步确定式(9)只有2个独立参数. 以式(6)~(9)为基础，得到新的波高周期联合概率密度函数 $ {f_{{\text{N}}}}(H,T) $，表达式为

(10) $ {f_{{\text{N}}}}(H,T) = \frac{{4H}}{{\sqrt {2{\text{π} }} H_{{\rm{s}}}^2 \cdot T \cdot {B_{{\rm{h}}}}}} \exp \left\{ { - \frac{{{{\left( {\ln T - {A_{{\rm{h}}}}} \right)}^2}}}{{2B_{{\rm{h}}}^{2}}} - \frac{{2{H^2}}}{{H_{{\rm{s}}}^{2}}}} \right\} . $

式中： ${A_{{\rm{h}}}}$、 ${B_{{\rm{h}}}}$均为波高的函数，由回归分析确定：

(11) $ {A_{{\rm{h}}}} = - 2\exp \left( { - 3.7\frac{H}{{{H_{{\rm{s}}}}}}} \right)+0.97\ln {T_{{\rm{z}}}}+0.331\;7 , $

(12) $ {B_{{\rm{h}}}} = - 0.47{\left( {\frac{H}{{{H_{{\rm{s}}}}}}} \right)^4}\exp \left( { - 1.7\frac{H}{{{H_{{\rm{s}}}}}}} \right)+0.28 . $

式(10)~(12)的输入参数为有义波高 ${H_{{\rm{s}}}}$和平均跨零周期 ${T_{{\rm{z}}}}$，可以根据波浪谱推导得到，也可以从波面高程时域数据中得到. 同时，这2个参数也是长期风浪统计资料(波浪散布图，wave scatter diagram)提供的数据^[2,14]. 因此，以 ${H_{{\rm{s}}}}$、 ${T_{{\rm{z}}}}$作为模型输入参数，便于在实际工程中使用.

在已知波浪谱 $ S\left( \omega \right) $的情况下，波面高程 $ \eta \left( t \right) $可以采用余弦波叠加法进行模拟^[15]. 将波浪谱 $ S\left( \omega \right) $在其频率上界 ${\omega _{{\rm{u}}}}$和频率下界 ${\omega _{{\rm{l}}}}$之间分成 $ N $等份. 令 $\Delta \omega = \left( {{\omega _{{\rm{u}}}} - {\omega _{{\rm{l}}}}} \right)/N$、 ${\omega _i} = {\omega _{\rm{l}}}+\left( {i - 1/2} \right)\Delta \omega$，则

(13) $ \eta \left( t \right) = \sum\limits_{i = 1}^N {\sqrt {2S\left( {{\omega _i}} \right)\Delta \omega } } \cos \; \left( {{\omega _i}t+{\theta _i}} \right) . $

式中： $ S\left( {{\omega _i}} \right) $为频率 $ {\omega _i} $对应的波浪谱值， $ {\theta _i} $为在 $ \left[ {0, \; 2{\text{π} }} \right] $均匀分布的 $ N $个随机数. 对式(13)进一步推导，采用离散傅里叶逆变换技术对波面高程时历进行数值模拟^[16]，流程图如图1所示. 模拟得到的波面高程 $ \eta \left( t \right) $满足各态历经性假设的平稳高斯随机过程^[1]，这与本研究提到的各种模型具有相同的理论基础. $ \eta \left( t \right) $的所有样本服从高斯分布，峰值、谷值样本服从Rice分布^[17]. 对 $ \eta \left( t \right) $的模拟数据进行高斯分布、Rice分布检验，确保仿真数据的正确性.

在得到满足要求的波面高程 $ \eta \left( t \right) $数据后，采用上跨零点法^[15]得到一系列单个波的波高和周期数据，再根据概率统计方法得到波高周期的经验联合概率密度. 对于给定的波浪谱，仿真得到2.0×10⁷个通过检验的 $ \eta \left( t \right) $样本. 经过上跨零点法处理后得到4.2×10⁵个单个波的波高和周期数据. 根据波高和周期仿真数据得到经验联合概率密度，并通过对同一波浪谱进行多次模拟的方式，确保经验联合概率密度的准确性.

1）P-M谱为单参数波浪谱，适用于充分发展的波浪^[1,4,12,15]，表达式为

(14) $ S(\omega ) = 0.008\;1\frac{{{g^2}}}{{{\omega ^5}}}\exp \left[ { - 0.032\;4{{\left( {\frac{g}{{{\omega ^2}{H_{{\rm{s}}}}}}} \right)}^2}} \right] . $

式中： $ g $为重力加速度. 本研究提出的波高周期联合概率密度函数需要的输入参数为 ${H_{{\rm{s}}}}$、 ${T_{{\rm{z}}}}$. 为了便于应用，根据平均跨零周期 ${T_{{\rm{z}}}}$的定义：

(15) $ {T_{{\rm{z}}}} = 2{\text{π} }\sqrt {{{{m_0}} \mathord{\left/ {\vphantom {{{m_0}} {{m_2}}}} \right. } {{m_2}}}} . $

可以推导出P-M谱的平均跨零周期表达式为

(16) $ {T_{\text{z}}} = \sqrt {{{200} \mathord{\left/ {\vphantom {{200} 9}} \right. } 9}} {{\text{π} }^{3/4}}\sqrt {{{{H_{{\rm{s}}}}} \mathord{\left/ {\vphantom {{{H_{s}}} g}} \right. } g}} . $

在已知有义波高的情况下，利用式(16)、(10)可以给出波高周期的联合概率密度函数.

2）布氏谱是二参数谱，不仅适用于充分发展的波浪，也适用于成长中的波浪^[1,4,12,15]，表达式为

(17) $ S\left( \omega \right) = \frac{5}{{16}}\frac{{\omega _{{\rm{p}}}^4}}{{{\omega ^5}}}H_{{\rm{s}}}^{2}\exp \left[ { - 1.25{{\left( {\frac{{{\omega _{{\rm{p}}}}}}{\omega }} \right)}^4}} \right] . $

根据式(15)可以推导出布氏谱的谱峰频率 ${\omega _{{\rm{p}}}}$与平均跨零周期 ${T_{{\rm{z}}}}$的关系，表达式为

(18) $ {\omega _{{\rm{p}}}} = \frac{1}{{{T_{{\rm{z}}}}}}{\left( {\frac{{64}}{5}{{\text{π} }^3}} \right)^{0.25}} . $

在已知有义波高和平均跨零周期的情况下，根据式(18)、(10)可以给出波高周期的联合概率密度函数.

为了检验各模型与经验联合概率密度的匹配程度，基于均方根误差提出模型有效性判别方法，表达式为

(19) $ {{\rm{RMSE}}} = \frac{1}{{n}}\sqrt {\sum\limits_{i = 1}^{n} {\sum\limits_{j = 1}^{n} {{{\left[ {{f_{\text{E}}}\left( {{H_i},{T_j}} \right) - {f_{\text{M}}}\left( {{H_i},{T_j}} \right)} \right]}^2}} } } , $

(20) $ {H_i} = {H_{\min }}+\left( {i - \frac{1}{2}} \right)\frac{{{H_{\max }} - {H_{\min }}}}{{n}} , $

(21) $ {T_i} = {T_{\min }}+\left( {i - \frac{1}{2}} \right)\frac{{{T_{\max }} - {T_{\min }}}}{{n}} . $

式中： $ {n} $为任意整数，取 $ {n} $=80； $ {f_{\text{E}}}\left( {{H_i},{T_j}} \right) $为经验联合概率密度函数值； $ {f_{\text{M}}}\left( {{H_i},{T_j}} \right) $为各模型的联合概率密度函数值； $ {H_{\max }} $、 $ {H_{\min }} $分别为单个波的最大波高和最小波高； $ {T_{\max }} $、 $ {T_{\min }} $分别为单个波的最大周期和最小周期. 由式(19)可知，RMSE数值越小，说明理论模型的分布与经验联合分布越接近.

根据全球波浪数据GWS中的4个海区（8、9、15和16）长期统计资料进行仿真^[14]，有义波高取值范围为[0.5, 13.5] m. 根据式(19)计算各模型的RMSE，结果分别如图2和表1所示. RMSE数值随有义波高的增大而不断减小. 随着有义波高 $ {H_{s}} $的增大， $H_{\rm{max}}-H_{\rm{min}}$和 $T_{\rm{max}}-T_{\rm{min}}$的数值不断增大，导致概率密度数值减小，RMSE数值也减小. 不同有义波高仿真的算例之间不具有比较意义. 对于不同的有义波高，各模型的RMSE具有相同的变化规律，按照RMSE从大到小排序为CNEXO模型^[8-9]、孙孚-Zheng模型、L-H 1983-Zheng模型^[7]、L-H 1983模型^[5]和孙孚模型^[6]、本研究模型，说明在不同的有义波高下，本研究模型比其他模型更接近仿真结果.

RMSE不能完全反映各模型在细节上的差异. 对比分析仿真结果与各模型概率密度函数发现，同一模型在不同海况中得到的结果相似. 选取 $ {H_{{\rm{s}}}} $=3 m时的结果进行展示, 如图3所示. 由图可知，本研究模型与仿真数据的匹配度最高，CNEXO模型^[8-9]与仿真数据的匹配度最低，L-H 1983模型^[5]和孙孚模型^[6]的分布形态较接近，Zheng修正的2种模型^[7]分布形态较接近. 如果不同模型的分布形态接近，那么RMSE (图2和表1)必然会相差不大. 分析图3还可以发现，各模型与仿真数据的匹配程度主要受周期T的影响. 对各模型的联合概率密度函数进行积分，得到如图4所示的周期的概率密度函数. 结合图3、4可以发现，L-H 1983模型^[5]、孙孚模型^[6]以及Zheng修正的2种模型^[7]存在大周期的情况，即模型预报存在周期较大的波，但在仿真数据中这种情况没有出现，Zheng等^[7]提出的模型修正在一定程度上改善了这种情况. CNEXO模型^[8-9]预报的周期整体偏小. 本研究模型无论是波高周期联合概率密度函数，还是周期概率密度函数都与仿真结果匹配较好. 分析结果表明，使用RMSE度量各模型与仿真结果的匹配度完全合理.

根据全球波浪数据(GWS)的 4个海区（8、9、15和16）长期统计资料进行仿真^[14]. 有义波高取值范围为[0.5, 13.5] m，平均跨零周期为[0.5, 15.5] s，进行280种海况的仿真. 在每种海况下，根据式(19)计算各模型的RMSE. $ {H_{\rm{s}}} $=0.5、6、10 m时，各模型的RMSE分别如图5和表2~4所示. 其中T_p为谱峰周期，T_p=2π/ω_p. 基于布氏谱仿真数据对比所得出的结果和P-M谱仿真数据所得出的结果一致，说明本研究模型对于这2种理论谱仿真的波高周期联合分布都可以给出较好的预报.

实测波浪谱与理论波浪谱存在差异. 为了进一步检验本研究模型的合理性，选取美国国家数据浮标中心2个不同浮标站具有代表性的4个实测波浪谱数据进行仿真分析. 实测波浪谱数据的基本情况如下. 1) 44066站(水深为78 m)：2021年11月26日18:00数据，有义波高为1.1 m，谱峰周期为6 s，无涌浪情况；2021年11月26日22:00数据，有义波高为1.3 m，谱峰周期为7 s，无涌浪情况. 2) 42098站(水深为12 m)：2021年11月27日7:00数据，有义波高为1 m，谱峰周期为5 s，无涌浪情况；2021年11月27日8:00数据，有义波高为1.1 m，谱峰周期为5 s，无涌浪情况. 将实测波浪谱数据绘图表示，如图6所示. 可以看出，实测波浪谱即使在仅风浪的情况下也会呈现出多峰形态，与经常使用的适用于风浪的单峰理论波浪谱有很大区别.

按照图6中的4个实测波浪谱进行仿真，并根据式(19)计算各模型的RMSE，结果如表5所示. 可以看出，对于每个实测波浪谱， CNEXO模型^[8-9]的RMSE最大，本研究模型RMSE最小；L-H 1983模型和孙孚模型的RMSE接近，Zheng修正的2种模型的RMSE接近. 各模型的RMSE接近这一结论与2种理论谱仿真数据得出的结论一致. 由表5可以发现新的现象：对于42098浮标站的实测波浪谱，Zheng修正的2种模型^[7]与仿真数据的匹配程度低于本研究模型、L-H 1983模型、孙孚模型. 对于44066浮标站的实测波浪谱，Zheng修正的2种模型与仿真数据的匹配程度仅次于本研究模型. 其中以44066浮标站实测波浪谱仿真数据所得出的结论与2种理论波浪谱得出的结论有差异.

为了进一步确认RMSE对各联合模型与仿真数据匹配度的描述，同时确认周期T是影响各模型匹配度的主要因素，取44066站2021年11月26日18:00波浪谱仿真数据与各模型波高周期联合概率密度函数进行对比分析，结果如图7所示，与各模型周期T概率密度函数的比较如图8所示. 由图7、8可以发现，L-H 1983模型^[5]、孙孚模型^[6]存在大周期的情况，Zheng修正的2种模型^[7]有效改善了这种情况，并在一定程度上提高了与仿真数据的匹配度. 结合图4、8可以发现，周期T概率密度函数的尾部(大周期)越接近仿真数据，模型的匹配度就越高. 由图7可以看出，L-H 1983模型^[5]和孙孚模型^[6]比较接近，Zheng修正的2种模型^[7]比较接近. 这与2种理论波浪谱给出的结论一致. 分析结果表明，本研究模型无论是波高周期联合概率密度函数还是周期概率密度函数，都与仿真结果匹配较好.

在以上理论谱和实测谱的分析过程中发现，L-H 1983模型^[5]和孙孚模型^[6]分布形态比较接近. 这主要是由于这2种模型采用的方法基本相同，不同点仅在于对负角速度情况的处理上. 虽然2种模型的公式表达形式存在一定差异，但在数值上比较接近. 潘锦嫦等^[10-11]对这2种模型进行分析后也得出类似结论. L-H 1983-Zheng模型、孙孚-Zheng模型因为提出者对模型进行了相同修正，所以在数值上比较接近. 对于深水情况，波浪谱的谱宽参数 $ \nu $一般情况下为0.3~0.5，均值为0.42^[11]. 本研究使用的各种波浪谱的谱宽为0.33~0.46，基本和实际海况相当. 说明本研究模型应用于实际海况具有一定的准确性. 对于浅水情况，波面高程 $ \eta \left( t \right) $不再服从高斯分布^[18]，可以采用如Hermite变换、Rychlic变换的方法^[19-20]，将本研究模型推广到非高斯情况. 由图4、8分析结果可以发现，本研究模型得出的周期分布预报精度高于现有的其他周期模型，在船舶耐波性设计、海洋工程结构振动特性分析方面有重要应用价值.

本研究提出新的短期风浪波高周期的联合分布模型，通过数值仿真验证了所提模型的准确性. 将所提模型与文献中常用的5种模型进行对比分析，得出以下结论：1) 本研究模型的精度较高，使用方便，但基于实测谱的计算精度略低于基于理论谱的计算精度；2) 由本研究模型推导出周期的概率密度函数具有较高的计算精度；3) L-H 1983模型的计算结果和孙孚模型的计算结果较接近，L-H 1983-Zheng模型的计算结果与孙孚-Zheng模型的计算结果较接近；4) L-H 1983-Zheng模型、孙孚-Zheng模型在短期风浪海况下，对波高周期联合分布的预报精度提高不明显. 此外，本研究模型可以推广到非高斯海况的应用场景中. 在后续的研究中，计划分析条件密度概率函数的一般表达公式，进一步提高模型的计算精度.