海水仿真在虚拟训练系统、3D游戏中有着广泛的应用。与地形绘制算法相比，动态的海水仿真算法要更为复杂，同时兼顾到海水仿真效果的真实性和实时性是一个很具有挑战性的任务，但是我们勇敢的接受了这个挑战。在这一章节里面，首先简单的介绍了正弦波动的情况，然后引入Gestner-wave以及Gestner-wave海水仿真算法，最后我们再通过一个巧妙而简单的推导出基于二维快速傅利叶变换的海水仿真算法。为了不至于曲高和寡，我这里没有涉及复杂的流体力学内容，而对于没有学过信号处理，不了解快速傅立叶变换的朋友，这也没有关系，我这里对算法的推导处理得很巧妙，你不需要任何信号处理的知识，只需要简单的初等数学知识就可以看明白。
1.言归正传，我们首先来介绍正弦波动的基本概念，以下这些参数是确定正弦波动的主要参数：
波动振幅（A）：波峰到波谷的高度差的二分之一
波长（L）：两个相邻波峰之间的距离
空间角频率向量（w）：空间角频率向量w的方向和波传播的方向相同，空间角频率向量的大小|w|和波长L的关系为：|w|=2*PI/L。
波速（S）：每秒钟波峰移动的距离。
时间角频率（wt）：时间角频率wt=S*2*PI/L。
波动传播方向（D）：波峰移动的方向。
初始相位（FI）：初始相位。
根据上述参数可以写出正弦波动的波动方程：
Y（x，y，t）=A*cos（w *（x，y）+ wt*t + FI）
根据正弦波动的等相位面，可以将正弦波动分为方向波和圆形波，等相位面为平行的平面的正弦波动，为方向波。而等相位面为一系列圆柱面的波动为圆形波。一般来说在风驱动的广阔的水体适合采用方向波，而对于一些较小的池塘，输入激励近似为点激励的情况下（例如扔一个小石子到河里面去，激起一圈一圈的波纹）则适合采用圆形波。在本章里面，主要涉及的是方向波。
2.Gestner-wave是最早的用于计算机图形学海水仿真的方法（1986年Fourier和Reeves）。P为海水表面上任意一点，没有经过波动时P点的位置为 （x0，0，z0），下列表达式为Gestner-wave波动计算公式：
（x，z）=（x0，z0）+ Q*A*sin（w *（x0，y0）+ wt*t + FI） （1）
y = A*cos（w *（x0，y0）+ wt*t + FI） （2）
这里波浪传播的时间角频率和空间角频率存在一个直接关系（在深海情况下，wt=sqrt（g*w），g为重力加速度），与前面的正弦波动稍有不同的是Gestner-wave在水平方向也有波动，而这种波动产生了波浪挤压的效果，而Q为控制波浪的陡峭程度，但是最好不要使Q超过1，否则会在表面出现环的效果。有时，为了模拟出更为真实的海水效果，采用多个正弦波叠加的方法来模拟复杂的波动效果。
3.单纯的基于Gestner-Wave的海洋仿真算法无法模拟出真实的海水效果，一个能够真实模拟出海水效果的经典算法是基于快速傅利叶变换的海水仿真算法（Jerry Tessendorf，SIGGRAPH Course Notes）。
下面我们来通过推导这个算法。首先分析一维的情况，考查一个长度为Ls的一维区域AB的海水波动情况，如下图所示。
这里有两个前提条件：
a：整个区域的波动情况具有周期为L的空间周期性。
b：我们希望使用多个Gestner-Wave的叠加来合成出整个波动效果。由于整个区域的波动效果具有周期为Ls的空间周期性，因此空间周期Ls必须为每个叠加的Gestner-wave的波长L的整数倍，如下图所示。
波长等于周期Ls
波长为空间周期的二分之一
波长为空间周期的三分之一
波长为空间周期的四分之一
。。。
波长为空间周期的n分之一
对于Gestner-wave，在已知波长L的情况下，可以确定出空间角频率w和时间角频率wt可以通过以下两个式子来确定。
w = 2*PI/L. （3）
wt = sqrt（g*w）. （4）
其中L=Ls/k；k为正整数。
在确定这两个参数之后，对于确定整个波动方程，还缺振幅和初始相位这两个参数。而对于振幅和初始相位这两个参数我们无法通过解析的方法来求出，对于怎样来获取这两个参数，我们将会在后面的方法里面介绍。
将这n个正弦波进行叠加之后得到一个复合的波动效果，如下式：
（6）
通过欧拉公式我们可以将上式改写为傅立叶变换的形式：
（7）
即
也就是
（9）
如果得到了初相位以及振幅的信息，我们就可以通过快速傅立叶变化的形式来进行整个海水波动计算。
在考察完一维的情况后，来考察二维的情况。选择一个长为Lx，宽为Ly的矩形区域，如下图所示。
a：整个区域的波动情况具有一定的空间周期性，在x方向上具有周期为Lx的空间周期性，在z方向上具有周期为Ly的空间周期性。
b：使用多个Gestner-Wave的叠加来合成出整个波动效果。
考虑一个沿AB方向传播的正弦波，要使得这个波动在x方向上具有Lx的空间周期性，则相当于点A和点C的相位差是2PI的整数倍。CE垂直于波传播的方向AB，因此C，E处在同一波阵面上，具有相同的相位。于是要求A和E的相差为2PI的整数倍，也就是AE的长度为正弦波的波长L的整数倍。
AC的长度为Lx，则AE的长度为Lx*cos（a）。
AE的长度为波长L的整数倍，即
Lx*cos（a）= k1*L， k1为整数。 （10）
同时波动在y方向也具有周期为Ly的空间周期性。
同样可以得到
Ly*sin（a）= k2*L，k2为整数。 （11）
其中k1和k2都为整数（可以是负整数）。
反过来，而对于任意一个给定的整数对（k1，k2），可以得到一个沿方向
（cos（a），sin（a））传播，波长为L的Gestner-Wave波满足我们空间周期性的要求。
在给定k1，k2的情况下，通过求解式（10）（11）方程组。
可以得到
Tan（a）= k2*Lx/（k1*Ly）。 （12）
a = atan（ k2*Lx/（k1*Ly））或者 a = PI + atan（ k2*Lx/（k1*Ly）） ， 这里把a的值范围限定在-PI<a<=PI.
这里a有两个值，而且相差180度，这个含义是，在一条直线上，有两个Gestner-Wave的传播方向在这条直线上，而且这两个波的传播方向相反。
根据a值可以解算出波长L。
L=1/sqrt（（k1/Lx）*（k1/Lx）+（k2/Lx）*（k2/Lx）） （13）
进一步我们可以求出一对空间角频率的向量。
w =（2*PI*k1/Lx，2*PI*k2/Ly）以及 -w =（-2*PI*k1/Lx，-2*PI*k2/Ly）
而根据传播关系，我们可以计算出时间角频率。
wt = sqrt（g*|w|）。
根据上述参数就可以写出（k1，k2）确定的一对Gestner-Wave
Y1= A0 * cos（ w *（x，y） + wt*t + FI0）
Y2= A1 * cos（ w *（x，y） - wt*t + FI1） （14）
最后我们选取多个这样的Gestner-Wave波进行叠加，选取整数对的范围是-N/2<k1<N/2, -M/2<k2<M/2. 这样的取值范围使得我们尽可能多的得到覆盖到任意方向上的波动。复合波动的计算公式为
（15）
再次使用欧拉公式将这个式子转化为二维福利叶变换的形式
（16）
这个公式就是使用二维快速傅立叶变换来计算海水波动的公式。 在后面的章节里，将要介绍使用海洋统计学频谱结合上述使用傅立叶变换描述的海水波动计算公式来进行海水模拟。
下面是本文的pdf格式下载
海水仿真算法论文（一）