大规模海浪交互模拟＊

刘 磊 丁剑飞 李 飞 杨 波 张 学

（石家庄机械化步兵学院 石家庄 050083）

1 引言

海浪的建模与绘制在电影、游戏、数字海洋、虚拟现实和虚拟战场环境等领域得到了广泛的应用。海浪本身具有实时性、高动态性、多变性、多维性和模糊边界环境性等方面因素，这都使得基于真实海洋环境仿真的研究还比较不成熟。随着图形处理硬件（GPU）的迅猛发展，增加海浪的真实感和实现海浪的快速绘制是一个急需并且可以解决的问题。同时，小范围的海浪模拟研究相对来说比较成熟，但考虑海面与船舶等水面运动物体双向交互作用的大规模海洋场景的仿真还有很广阔的研究价值。

关于海洋场景的仿真，有基于物理的方法，如比较经典的SPH 模型［1～3］，它能较好地模拟物体和水体的交互作用，但只适合范围较小、非实时的海洋场景仿真；基于高度场的仿真算法［4～5］，能较好地仿真海洋的表面，但是很少有作者对物体与水体的交互作用进行考虑；有基于粒子系统的仿真算法，能取得一定的表面效果，但是生成的海洋场景具有人工因素，缺乏相应的物理基础，从而缺乏真实性。

本文提出了一种海浪交互模拟的有效方法。首先采用复合动态网格来仿真海面网格；其次，采用基于频谱统计的FFT算法对非交互波浪进行仿真，求得非交互波浪的高度场；接着，从刚体力学角度分析了交互波产生的位置和时机，通过引入仿真引擎和求解二维波动方程来实现对交互波的仿真，求得交互波的高度场。最后对两种波的高度场进行叠加，从而实现对海洋场景的整体仿真。

2 海面波浪仿真

2.1 波浪分类

我们把海洋上的波浪分为两类，一类是交互波浪，就是因为和物体交互产生的波浪；另一类是非交互波浪，就是不存在和物体的交互。

2.2 波浪模拟

2.2.1 非交互波浪模拟

便可以计算梯度向量，微分得到：

用这种FFT方法生成的水面的高度图如图1所示。

图1 高度场

图2 RGB相关梯度图

图3 实际海面效果图

图4 加入光照效果图

为了在GPU中实现更高效的计算，我们在上述结果的基础上，直接由高度图生成与其相对应的R、G、B三个通道梯度图。这样不仅可以减少法向的计算量，还可以减少纹理所占用的内存和显存空间。其效果如图2所示。最终渲染结果如图3所示。对于高光部分的海面效果，我们利用了色调映射的方法增强了光照的效果，使得能在非常亮的地方看到高对比度的图像如图4。

2.2.2 交互波浪模拟

交互波浪是由海面上的物体和海面交互而产生的波浪。它是由运动对象压迫水体产生的。对于这种波浪的模拟，我们一般采用求解基于流体的Navier-Stokes方程组来进行模拟。由于Navier-Stokes方程组是包含连续性方程、动量守恒方程和能量守恒方程的方程组，求解起来比较复杂。我们对其进行简化，采用求解二维的Navier-Stokes方程来模拟交互波浪。该二维波动方程只考虑水平方向速度分布，而忽略垂直方向速度来求解水面高度场，这是应用比较广泛的一种形式。

恒定张力作用下的二维表面质点的运动一般用二维波动偏微分方程表示，假设三维空间坐标中XZ平面为水平面，水面高度场为Y轴正向方向，则含粘性阻尼力的二维波动方程如下：

上式中，高度场y（x，z，t）为水平面坐标（x，z）和时间t的偏微分方程，c、μ分别是水波速度和阻尼系数。求解二维波动方程的方法有很多，我们采用有限差分方法（FMD）来对其进行求解。因为有限差分方法过程直观且容易采取GPU进行求解。我们把二维水平面进行离散，划分为规格均匀网格，每一个网格的坐标为（i，j），对于第k帧其中包含的值为y＝（xi，zj，tk）。应用有限差分方法，我们可以近似求得上式中的各个组成部分：

3 船行波的模拟

3.1 刚体力学

我们知道，视觉的真实性不单单依赖于液体仿真，而且更加依赖于液体和物体交互的可信度。我们把交互分为：1）物体对液体的交互；2）液体对物体的交互。前者描述了漂浮物体的受力和它对液体运动的影响。后者描述了液体对运动的船只产生的作用或对物体的影响。由于液体对船的影响相对来说较小，并且研究起来涉及比较复杂的力学分析，所以我们不予分析。

其中h（x，y）：ℝ2→ℝ描述了液体在位置（x，z）的高度，并且x＝（xx，xy，xz）T。Mtrot根据角加速度随时间步长Δt做如下变化：

其中ΔMrot是Δt·‖M‖关于矢量M／‖M‖的转动矩阵（M是角动量）。一个粒子集合的角动量等于每个粒子角动量的和。

因此，ΔMrot可以由每个粒子角动量的总和所决定。因此，我们可以仿真具有不同物理参数的任意形状的物体。

下面我们来研究一下其稳定性。由于为求得高性能，我们对刚体物体只使用了一个相对较小的离散化，所以当一个粒子通过水面时，在水面，浮力的不连续会导致力的瞬间变化（突变）。为了避免这样的边缘效果，浮力随水下距离dt线性变化，这样浮力就会在水面变化为0。这样刚体的运动就会变得平滑，即使对于稀疏的样品物体，亦是如此。

3.2 船的仿真

因为我们的方法使用粒子来模拟刚体物体，所以在仿真中我们可以控制物体为任意的形状。明显的，我们的性能会随着粒子数目的增多而降低。为了减少粒子的数目，我们仅仅模拟船身。因为物理计算建立在每一个粒子的基础上，所以这种方法会在非交互波的的基础上产生真实的运动。

3.3 交互波产生的时机和位置

因为我们不仿真液体流，所以来自刚体的波浪产生仅仅当物体下沉或运动时才会发生。对于下沉物体的情形，当物体正在横穿水面时，一种拥有物体和水面交互形状的波会产生。为了决定波的开始时间，每一个物体需要知道它的状态s：（s＝0）在水面之上，（s＝1）在水面之下。这个状态在每一时间步刷新，并且当物体开始横穿表面时，波会产生。

如果对于一个给定的刚体物体，式（12）等于1，那么物体已经横穿了水表面，并且产生了波。因此，物体与表面平面的交互投影到仿真网格上，根据影响的强度来确定幅度。

对于由运动的物体引起而产生的波，我们充分利用下面的物理特性。对于一个运动船体，船头的压力比船尾大。这样，在船的前面会产生正向的波，在船的后面会产生负向的波。产生波的振幅依赖于船的速度。为了确定船在一个时间步长内通过的区域，我们把当前和先前时间步长内的位置和方向结合起来。对每一时间步长我们产生一个包括船体和水面交互的图。因为计算精确的交互区域花销太大，因此我们做两个假设：1）我们近似表面为平面。因为非交互波和海洋表面的面积比起来通常很小，因此这是一个合理的假设。然而在具有大量非交互波的场景中，这个假设可能会成为一个瓶颈并且导致加工的痕迹。在最终的绘制中，这个差异几乎是不可见的，因为船在水中，并且随船体产生的波浪叠加到非交互波上面。2）我们假设每一个船的切面在y＝常数内。y＜0是y＝0切面的子集。切面定义为多边形船和平面y＝常数的横切面。对于大多数水面的小船，甚至双体船来说，这个假设都是有效的。对于被拉出水面的比较复杂的漂浮物体来说，需要采用一种更加高级的方法。例如，通过使用两个平面y＝±ε进行剔除，横切面可以被精确地提取出来。因此，即使对拥有复杂形状的正在剧烈运动的物体来说，它们也将在边界层横断面产生波。在横断面st的纹理产生以后，把它和来自先前时间步的对应的纹理相结合。考虑纹理作为一个集合，三个不同的波产生状态可以如下建立：

图5 船行波的叠加

1）头波：st＼st－1

2）船体：Bt＝st∩st－1

3）尾波：st－1＼st

为了产生波，波浪强度可能从实际仿真网格值中添加或剔除。因此，当仿真的船和存在的波相互影响时，产生的波和存在的波便会相互干涉，考虑这是必要的。如果这种效果不需要，那么仿真网格中的值可以很简单地被替换。如果被替换，那么它们将轻微地变化，并且和前面的变化相一致。这样，在船外形上的人为加工感会降低并且会产生头波。如果船速比波的传播速度快，那么头波会始终保持在前面。为了更加逼真，这个区域可以通过如今电脑游戏熟知的飞溅效果加以改善。

3.4 交互波仿真

对于头波和尾波，我们使用干扰引擎来对其进行仿真。为了能大体模拟上述船行波效果，我们设计一个较大的船头波干扰引擎和m个较小的船尾波引擎，并且设定第i个船尾波干扰引擎的位置ez为：

其中，w为船的宽度。在具体的实验中，干扰引擎的幅度和广度是由船舶的重量大小决定的，因此需要根据船体的形状来构造船头和船尾波干扰引擎，在船行驶过程中对水面高度场进行干扰。下面我们分别对船头波和船尾波设计不同的平滑扰动模板来实现船行波的平滑。现实中，船在行驶过程中，从船头向后左右两侧形成两股大波，类似于人字形波。所以我们应用一个幅度和范围较大的扰动模板算子D（x，z）。并且通过在网格的不同位置设置不同的扰动系数来实现人形波的平滑。我们如下设定D（x，z）。

相应的，在船尾处，我们设计m个范围和幅度相对较小的扰动引擎。对于其扰动平滑模板，我们把网格中心点扰动系数设计为0.3，与中心点相邻的四个点系数设计为0.15。则船尾波引擎模板矩阵算子F（x，z）为：

图6 船只左转弯时海洋场景仿真

图7 船只右转弯时的海洋场景仿真

在完成了上述的仿真后，我们采用标准的光照模型和立方体映射以及文献［6～7］中的反射和折射绘制技术来实现对水体中物体的折射和反射现象的仿真。最终的仿真结果如图6和图7。

4 实验

在基于上述模型的理论基础上，我们进行了仿真实验工作。本文所采用的实验平台为Intel Xeon CPU 16核2.27GHz，主存为12G，显卡采用的是GeForce GTX295，显卡内存为1800M，操作系统为 Windows XP，整个实现基于 OpenGL／Cg。在绘制时，我们还加入了色调映射（Tone Map）［8］和高动态范围（HDR）［9］图像处理技术等特效来加强海面场景绘制时的光照效果的真实感。我们得到了有关的实验数据如表1和表2。

表1 不同网格步长下的帧速率

表2 相同网格尺寸下各部分仿真时间消耗所占百分比

对于表1来说，我们可以看出帧速率随网格步长的增大而增大。也就是说帧速率随着分辨率的增加而降低，但总体来看，帧速率相对来说还是比较高的，能满足实时性要求。对于表2来说，表面提取和渲染所花费的时间较多，尽管它能在GPU上有效地执行。我们的方法需要不断地从GPU中读回。因为刚体仿真是在CPU上执行的，所以包围刚体的局部区域需要不断地从GPU读回来实现刚体物体的仿真。

5 结论与展望

本文给出了一种仿真拥有交互物体的大规模海洋场景的方法，采用复合动态网格来仿真海面网格；采用基于频谱统计的FFT算法求得非交互波浪的高度场；从刚体力学角度分析了交互波产生的位置和时机，通过引入仿真引擎和求解二维波动方程来求得交互波的高度场。最后对两种波的高度场进行叠加，从而实现对海洋场景的整体仿真。将来可能改进的工作：在需要的区域，利用基于3D Navier Stoke方程［10］来实现对流的仿真；仿真水对船的交互；添加3D液体效果，如液体飞溅等；进一步利用GPU提高仿真的效率。

［1］He Yan，Zhang ye Wang，Jian He，et al.Real-time fluid simulation with adaptive SPH［J］.Computer Animation and Virtual Worlds，2009，20（5）：417～426

［2］W.T.Reeves.Particle systems：a technique for modeling a class of fuzzy objects［C］／／Proc.of ACM SIGGRAPH，1983：359～376

［3］M.Muller，B.Solenthaler，R.Keiser，et al.Particlebased fluid interaction［C］／／Proc.of the ACM SIGGRAPH／Eurographics Symposium on Computer Animation，2005：237～244

［4］M.Kass，G.Miller，Rapid.Stable Fluid Dynamics for Computer Graphics.Computer Graphics（Proc.of SIGGRAPH'90），1990：49～57

［5］Gary A.Mastin，Peter A.Watterger，John F.Mareda.Fourier Synthesis of Ocean Scenes.IEEE Computer Graphics ＆ Application，1987：16～23

［6］SOUSA T.Generic refraction simulation.In GPU Gems 2，2005：295～305

［7］JOHANSON C.Real-time water rendering-introducing the projected grid concept［C］／／Master of Science Thesis（Lund University），2004：96～104

［8］Frank Losasso，Hugues Hoppe，Geometry Clipmaps.Terrain Rendering Using Nested Regular Grids［C］／／Proceedings of ACM SIGGRAPH，2004，23（3）：769～776

［9］Jeff McGee.Real-Time Water Rendering［C］／／proceedings of CSCE5933Topics in Game Development，April 3，2007

［10］Kass M，Miller G.Rapid，stable fluid dynamics for computer graphics［J］.Computer Graphics，1990，24（4）：49～57