随机波数值模拟方法1 概述研究海浪及其对工程的作用有三种途径：一是现场观测研究；二是在实验室内进行模拟研究；三是理论分析研究。

由于海浪的复杂多变性，加上现场环境恶劣，进行现场观测需花费大量的人力物力；理论研究目前也有较大的局限性，特别是对于不规则波浪，很多问题有赖于室内的模拟研究。

模拟研究的方法可分为两大类。

开始是在水槽或水池内利用风或造波机进行物理模拟，亦即进行波浪模型试验。

在人们的精心设计下，可以把负责的现象分解为多个简单的模型，然后再把成果综合起来。

过去已取得了大量的研究成果，目前仍是主要的研究方法之一。

随着电子计算机的发展和普及，海浪的数值模拟得到迅速的发展，它具有经济方便等优点，日益受到人们的重视和广泛的应用。

天然海浪是很复杂的，人们对它的认识和研究过程是由简到繁，由浅入深，及即由单向规则波—斜向规则波—单向不规则波—多向不规则波。

2 不规则波浪的数值模拟—模拟频谱单向不规则波浪的数值模拟方法，大多建立在线性波浪理论的基础上。

本文主要介绍利用线性叠加法和线性过滤法进行二维不规则波的模拟。

2.1 线性波浪叠加法在工程中，如果已经得到了特征波的波参数如有效波高H s、周期T 等参数，如何得到一列不规则波面时间历程呢？一般通过模拟靶谱法来完成。

将有效波高H s、周期T 等参数代入某波浪频谱形式中，得到的海浪谱即为靶谱。

现在要模拟某波面不规则波面时间历程，使得模拟的波谱同靶谱一致。

平稳海况下的海浪可视为平稳的具有各态历经性的随机过程，波动可以看作无限多个振幅不等、频率不等、初相位随机的简谐余弦波叠加而成，即Mt a i cos k i x i t i (1i 1式中，t 为波动水面相对于静水面的瞬时高度；a i 为第i 个组成波的振幅；k i, i为第i个组成波的波数和圆频率；k i 2 L i , i 2 T iL ，T 分别为波长、周期；x，t 分别表示位置和时间，通常固定位置，可取x=0；i 为第i 个组成波的初位相，此处取在（0,2 π）范围内均布的随机数。

通过频谱来模拟海浪，设欲模拟的对象谱（靶谱）S 的能量绝大部分分布在L ~ H 范围内其余部分可忽略不计。

把频率范围划分为M 个区间，其间距为i i i1，取?i i1 i 2，则第i个组成波的振幅为(2)则将代表M 个区间内波能的M 个余弦波动叠加起来，即得海浪的波面：Mt 2S ?i i cos %i t ii1式中，~i 为第i 个组成波的代表频率。

用波浪叠加法模拟海浪时应注意以下几点：2.1.1 频谱范围L ~ H 的选取频谱范围L ~ H 的选取，取决于所要求的精度。

设在频谱高低侧各允许略去总能量的部分（例如取0.2%），对于可积分的谱，易于确定L 和H(3)若采用公式SL 0.785 exp3.114H s2表示的P-M 谱，可以得到3.11H s2 ln1/4，H1/43.11H s2 ln 1 (4)a i对不可积分的谱，可以采用数值计算的方法来确定H 。

首先采用数值积分的方法计算波浪频谱的总能量E，然后计算对应每个频率i 的累积能量E i ，则E i/E 对应的频率即为下限L，E i/ E 1 对应的频率即为上限H 应该看到，在M 一定的情况下，不恰当地增大谱频范围，反而会使精度下降。

一般取谱峰频率的3~4倍作为H 已足够。

图1 划分波谱的频率区间示意图2.1.2 频率区间的划分划分频率区间的方法，有等分频率和等分能量法。

2.1.2.1 等分频率法下面简要介绍下等分频率法。

取H L M （一般取M=50~100）。

但若采用式?i i1 i 2中的?i作为i 区间的代表频率，则由式（3）模拟所得的波浪将以周期 2 重复出现，除非值足够小；否则与实际的海浪情况不符。

应在各区间内部随机选取频率作为该区间的代表频率~i。

~i的选取方法对模拟结果有相当的影响。

由于波能集中在谱峰部，如M 值较小；只有少数位于谱峰处的组成波起主要作用，可能产生较大的误差。

2.1.2.2 等分能量法定义累积谱为ES 0d(5)如果按照等分能量法分成N 份，则分界频率i可以用下式来确定。

E iiE im 0 (6)NN对 P-M 等分可积分的谱， 则B1/4(7)iln N / i各组成波的振幅 a 相等2m 0 Na i 2S ?ii(8)此时式 ta n cosi1ntn变为t2m 0 NcosN i 1 ?i ti(9)?i 1 i2(10)2.1.3 随机相位的选取随机初位相 i 应在0~ 2 区间内均布。

如组成波数 M 不很大，则由计算机 产生的随机数往往不够均布，影响模拟结果。

我们采用人造的比较均匀的随机数，模拟结果较好。

合田采用M=200，由计算机产生随机数（每次不同）进行多次重复计算，对结果进行统计分析，取其 特征值。

2.2 线性过滤法应用线性滤波法模拟海浪的基本思路是： 以白噪声为一线性系统的输入， 通 过选择适当的系统函数使该系统输出的谱恰恰等于靶谱海浪等随机过程由多种不同频率的成分组成， 他们可以通过不同的滤波器分 离开来。

如图 2 所示，只有高频信号能通过高通滤波器， 通过低通滤波器的是低(15)频信号，允许一定频率范围内的信号通过的滤波器称做带通滤波器具有如图 3 中所示传递函数的滤波器称做成型滤波器。

这些滤波器可以是数 字式的，也可以由硬件组成。

线性系统的输入谱S * f 和输出谱 S *yy f 之间存在 下列关系：2S *yy fT f 2S *xxf (11)白噪声的谱密度为常数，且可等于 1，如将它作为输入，通过按靶谱设计的 成型滤波器后， 即可得到谱形符合靶谱的随机波浪。

因此线性过滤法的关键在于 靶谱设计过滤器。

过滤器的选择。

输入白噪声的谱 S *xx f 1 ，要模拟的波浪靶谱为 S * f (双 xx侧谱)，由上式得过滤器的传递函数为T fS * f S f 2 (12)在时域，线性系统的输入。

输出函数间有关系即t x t h d(13)h 是脉冲响应函数，也是过滤器的权函数，其傅里叶变换即为传递函数T(f)，即T f e i2 f df写成离散形式：La j x t j t t 0, t,2 t ,..., N t jL(14)图2 滤波原理示意图模拟不规则波浪。

将上式代入可得到所要的波面。

为便于计算，把它改写成Lt A0 x L t A j x L i j x L i j (1图 3 用过滤法模拟波浪示意图式中，x L t相当于x t L t ，可取L=20~30白噪声x(t) 可用一系列独立的正态分布的变量x1,x2,... 来接近，这些变量的均值为零，方差为1。

可按下式得到：nx k 2RAN i 1 3 n ，k 1, 2,3,... N 2Li1RAN i 为在( 0,1)区间内均布的伪随机数，一般计算机可直接产生。

可取n=30~50。

3 程序实现3.1 程序一：线性波浪叠加法模拟频谱%% 不规则波浪的数值模拟—模拟频谱%% 线性波浪叠加法t=0:0.01:1000; % 时间间隔x=0; % 初始尾椎a=abs(randn(1,3));% 幅值T=abs(randn(1,3)); % 周期L=abs(randn(1,3)); % 波长c=abs(rand(1,3))*2*pi; % 初相位A=0;for i=1:3A=A+a(i)*cos(2*pi*x/L(i) -2*pi*t/T(i)+c(i));end plot(A);set(gca,'xlim' ,[0 100000]);set(gca,'ylim' ,[ -5 5]); xlabel('\itTime');ylabel('\itAltitude' ); title( '\bf 线性叠加法模拟频谱');grid on %% end图4 线性波浪叠加法模拟频谱图（N=3）图7 线性波浪叠加法模拟频谱图（N=200 ）图 5 线性波浪叠加法模拟频谱图（N=10）图 6 线性波浪叠加法模拟频谱图（N=50）3.2 程序二：线性过滤法%% 随机脉冲函响应数num=[0,0,25]; den=[1,4,25];[y,~]=impulse(num,den); plot(randn(1)*y,'r:','linewidth' ,1.5); set(gca',xlim' ,[0 127]);ylabel('\itAmplitude' );title( '\bfRandom Impulse Responcse)' grid on图8 线性波浪叠加法模拟频谱图（N=1000 ）T=fft(y,127); plot(real(T),':','linewidth',1.5);set(gca',xlim' ,[0 127]); ylabel('\itAmplitude(Real/Imag)' ); hold onplot(imag(T),'r:','linewidth',1.5); title( '\bfTransfer Function'); grid on图 9 生成随机脉冲响应函数图 %% 生成传递函数M=127;dt=1/(127*2);t=0:dt:1000;F=1/(2*dt);df=F/M;L=20;A=0;for j=1:Lfor i=1:MA=A+1/M*real(T(i*df))*cos(1*pi*i/M); B(j)=A;A=0endend stem(B,'filled' ); ylabel('\itAmplitude' ); xlabel('\itj' );title( '\bfA(j)' );grid on图 10 生成传递函数实部、虚部图 %% A(j)%% 白噪声的模拟x=normrnd(0,1,1,167);plot(x, 'linewidth' ,1.5);ylabel('\itAmplitude' );title( '\bfWhite Gaussian Noise)' ;set(gca,'xlim' ,[0 167]); grid on图11 生成A(j) 数值图图12 模拟高斯白噪声图%% 不规则波浪的模拟A0=randn(1);N1=0;N2=0;for i=1:Mfor j=1:L if (i-j<=0) N1=N1+A(1)*(x(i+j)); elseN2=N2+A(1)*(x((i+j))+x((i -j)));endendendN=N1+N2;t1=t+20*0.039;z=round(t1);H=A0*x(z)+N; for i=1:25400;if (H(i+1)==H(i))H(i+1)=H(i); H(i)=0;end end H=H(H~=0);plot(H, 'linewidth',1.5); ylabel('\itAmplitude' );xlabel( '\itTime' ); title( '\bf 线性过滤法模拟频谱');set(gca,'xlim',[0 100]); grid on %% end图13 线性滤波法模拟频谱图3 线性海浪数值模拟的其他方法除上述的海浪模拟的主要方法外，还有一些其他方法。