基于GPU粒子系统的可交互水波实时模拟

张艺衡，李晓峰，2，刘小玲，杨智宇

（1.四川大学 计算机学院，四川 成都610064；2.四川大学 视觉合成图形图像技术国防重点学科实验室，四川 成都610064）

0 引 言

对于可交互水波等流体现象的模拟，目前比较典型的方法是数值求解纳维斯托克方程。该方法虽然能够模拟出具有高真实感的水波，但其数值求解过程需要大量的计算资源，不适用于大规模场景中的可交互水波实时模拟。另外，基于波动方程数值求解的水波模拟也是目前应用比较广泛的算法之一，与基于纳维斯托克方程方法类似的是，其数据求解的过程依然会产生庞大的计算开销，算法的绘制效率偏低。这2类方程相较而言，波动方程更加简单，对于流体的表面模拟也能达到比较真实的效果。因此，本文将波动方程作为算法的理论基础，并对该方程进行近似建模，构建水面高度场，并结合粒子系统进行波浪的实时更新，以实现具有高效率和高真实感的可交互水波实时模拟。

1 相关工作

在水体模拟的早期研究中，主要是针对水面形态及运动状态直接进行显示建模，这类方法一般采用三角函数、噪音函数和傅里叶合成公式来进行水表面的模拟建模。Schneider等人利用Perlin噪音函数实现了基于GPU的实时水波模拟。Tessendorf提出采用傅里叶合成法进行海浪运动的模拟，并且该方法成为目前模拟环境波的主流算法之一。但是，该算法不能模拟水面浮游固体与水面之间的交互运动。为此，Tessendorf提出基于该算法的改进算法，实现了浮游固体对水面的单向交互。水面高度场构建法也是水体模拟的主流算法之一，这类算法主要以浅水方程作为模拟的理论基础，出现了许多研究成果。T.R.Hagen等人采用半拉格朗日求解法求解浅水方程构造水面高度图，但该方法在效果逼真度方面的损耗较大。为此，B.M.Kim等人对该算法进行了改进，取得了一定程度的算法逼真度的提高。Daniel等人提出了一种基于浅水方程的自适应网格模拟算法，实现了大规模水面的模拟［1］。Hyokwang等人采用二维SPH实现了水面高度图的创建［2］。Cem Yukel提出WP（wave particles）概念，将水波离散成满足波动方程的波浪粒子，以构建水面高度场［3］。Hilko等人对WP算法进行了扩展，实现了水流效果的模拟［4］。另外，欧拉网格法在水波模拟方面也出现了许多研究成果。Enright等人提出了基于欧拉网格的PLSM算法，Schroeder介绍了对该算法的一些改进，在一定程度上降低了计算复杂度［5］。Andrew等人采用超网格单元 （tall grid cells）思想，将N－S方程与高度场相结合，减少了单元网格的数量，实现了固体与水之间的交互［6］。四面体网格在欧拉法流体模拟中也得到了应用［7－9］。除此之外，拉格朗日粒子法在水模拟方面也得到了非常广泛的应用，如经典SPH流体模拟法。Becker等人基于SPH算法，并对其进行了改进，实现了体积近似恒定的流体模拟［10］。Solenthaler等人则提出了预期修正SPH算法 （PCISPH），保持稳定时间步长和迭代求解，使粒子系统趋向于不可压缩状态［11］。Clavet等人实现了基于拉格朗日粒子的粘性流体的实时交互模拟［12］。拉格朗日粒子法易于表达和实现，相对于欧拉网格法，粒子运动具有自由和随机性，不需要对空间进行有限网格的定义。但拉格朗日法随着粒子数量的增加，其计算量往往会随之增大。

2 算法思想及设计

2.1 高度场离散化

对于水面波浪的模拟，本文采用表面高度场的方法进行。对于一个完整的水面高度场，通常包含了数量不定的多个波浪高度信息，对整体高度场进行直接构建的方法可行性不高。因此，本文采用高度场离散化的思想来构建水面高度场。即将整体高度场以波浪为单位进行划分，更进一步将单个波浪划分成多个波段，形成多个小的局部高度场，从而将问题简化为对局部高度场的独立构建。如图1所示，对于左图中的完整波浪，本文将其离散为右图中方框标记的局部高度场，并对局部高度场进行一定规则的排列控制，以形成左图所示的完整平滑波浪高度场。

图1 高度场离散化

2.2 局部高度场建模

本文采用离散化的构建思想，将完整水波离散为多个离散小型水波波段，并结合粒子系统思想，每个波段对应一个粒子。通过更新粒子属性的方式，达到模拟跟踪水波运动形态的目的。在现实生活中，水波由横波和纵波2个部分组成，分别对应水波的垂直偏移和水平偏移。与此类似，在本文中，每个局部波段对应一个局部高度场，由垂直偏移和水平偏移两部分组成，分别对应一个局部水平偏移函数和一个局部垂直偏移函数。

假设局部垂直偏移函数为D（X，t），则该函数表示在t时刻，水面上X位置处的垂直偏移，即相对于水面水平位置的高度值。假设水面高度场为函数 H（x，t），则 H（x，t）的计算公式如下所示

式中：h0——水面的水平高度，φz（X，t）——t时刻，p（x，y）位置处由水波引起的垂直偏移，Di——第i个粒子表示的垂直偏移。对于D（X，t），本文提出一种径向定义方式，即将每个局部高度场构造为底面为圆形的弧面，每个弧面以波浪圆心为中心，随着波浪的传播进行径向运动，以确保各个弧面之间的相对独立，以充分利用粒子系统进行整体波形及运动状态的模拟。具体构造方式如下所示

式中：ai——局部高度场对应波段的振幅，ri——弧面底部圆面的半径，｜x－xi（t）｜——圆面上某点与圆心之间的距离。通过这种方式构造的局部高度场，以一定的间隔进行弧形排列，即能构造出一个完整逼真的波形高度场。

局部垂直偏移函数反应了波形在垂直方向上的形态，局部水平偏移函数则模拟了波浪运动过程中在水平方向上的变化，如下所示

2.3 波粒子系统

本文算法的关键思想在于离散水面高度场，通过控制各个离散高度场的高度变化及运动方向控制，以达到模拟真实水波的形态变化及运动规律。对此，本文引入粒子系统，进行水波运动变化的实时跟踪。即每个局部高度场对应一个粒子，粒子的属性存储了各个高度场的信息，包括产生时间，产生位置，扩散角，传播方向以及振幅。随着水波的传播扩散，通过更新波粒子的属性来模拟水波的运动规律及形态变化。但是，随着水波的扩散，波粒子进行径向运动，粒子之间的距离也随之增大，当距离达到一定阈值时，水波波面会出现相对分离的现象，即出现失真的情况，如图1右图所示。为了避免这样的失真现象，随着水波的扩散，本文通过分裂波粒子的方式，在粒子之间添加新粒子从而控制相邻粒子之间的间隔距离，如图2所示。为了提高算法效率，本文提出一种基于母粒子的更新方法。

图2 波粒子分裂

2.4 基于母粒子的粒子更新

在本文算法中，一个完整的水波由多个波粒子所代表的圆底弧面，以一定的均匀间隔排列形成。即同一个水波上的波粒子，具有相等的振幅及运动速度。因此，同一水波上的波粒子也具有同时分裂和同时消亡的特征。据此，本文采用二维的方式组织粒子，即每个水波对应一个母粒子，在进行粒子系统更新时，只需要检测母粒子，当母粒子达到分裂或消亡的条件时，则分裂或消亡母粒子下的所有子粒子。除此之外，粒子振幅的更新也只需要针对母粒子进行实时更新计算，所有子粒子与母粒子具有相等的振幅，从而避免了遍历更新所有的波粒子，极大地降低了计算量，实现高效率的粒子系统更新。至于粒子的实时坐标位置，本文通过对粒子属性进行巧妙的设计，通过设置粒子的产生时间，运动速率及扩散角，使其可充分利用可编程着色语言GLSL，在GPU上进行高效快速的并行计算。避免了在CPU上对每个粒子进行坐标位置的实时计算，利用GPU强大的并行计算能力达到高效的计算效率。

2.5 水的光学特性模拟

为了实现具有高真实感的水面效果，本文基于水面的基本光照模拟，引入环境贴图技术。基于基本光照物理模型，水面光照模型如图3所示。水面基本的光照效果主要是折射和反射，且其计算方法基于真实水面的折射和反射物理定律。

在现实中，水面对周围环境具有非常强的反射性，对水面反射性的模拟对水面的真实感具有非常重要的作用。在虚拟仿真领域，对水面的反射效果模拟最为精确且经典的技术是光线跟踪技术，但是该技术的计算开销非常大，不适用于实时绘制系统中。为此，本文采用立方体环境贴图技术进行模拟。该技术不仅可以充分利用GLSL的可编程属性高效地模拟出高质量的光照效果和环境贴图，并且基于该技术可以比较方便地生成远反射效果。

3 算法实现

第2小节介绍了本文的算法思想，本节将介绍波粒子系统与高度场结合模拟水波的算法实现过程。本文的波粒子系统主要负责根据当前水面上的交互信息产生新粒子，并实时更新系统中波粒子的位置，处理粒子传播过程中的分裂和反弹。反弹主要针对波浪遇到水体边界时的处理。对于粒子系统的迭代更新，传统的基于CPU的更新方式实现简单，但效率很低。为了充分利用GPU强大的并行计算能力，本文对粒子属性进行了巧妙的设计，使得粒子的更新处理可以全部在GPU上进行，获得了非常高效的计算效率。

3.1 母波粒子的产生

对于本文基于波粒子系统的水波，水波从产生到扩散，最后到消失的过程实际上对应了波粒子的生成，分裂和消失的过程。若水面上某处发生了交互，将对应一个波粒子的生成，其属性包括：生成时间、生成位置、扩散角、传播方向和振幅。通过这几个属性的设计，使得粒子的更新可以完全通过GPU来进行，具体方式将在3.2小节进行介绍。

3.2 波粒子更新

波粒子的运动基于二维空间且粒子之间相对独立，不存在交互。除此之外，没有其它外力作用在该粒子系统上，因此波粒子具有相对恒定的运动速率，并且除了与水体边界发生的碰撞反弹之外，粒子的运动方向保持不变。基于本文对粒子属性的设计，不需要每帧对波粒子的位置进行更新存储，而是通过GLSL实时计算获得，如下所示

式中：v——波粒子的运动速率， ——波粒子的传播方向。由于粒子数量随着波浪的传播呈指数级增长，甚至导致其更新存储所需的CPU缓存大于系统的CPU缓存，以致算法异常退出。通过利用GPU实时并行计算的方式，避免了这部分内存的读写开销，而消耗小部分式 （4）所需的乘法和加法开销，确保了算法的高效性和稳定性。

在现实中，水波在传播的过程中，由于运动扩散和阻尼导致其振幅的不断降低直至降低到0。为此，本文采用式（5）来进行模拟

式中：δ——阻尼系数，t——当前时刻，t0——波浪粒子的产生时刻。

3.3 波粒子分裂

随着水波的传播，同一水波上相邻波粒子之间的间隔增大，为了确保波形的平滑连续，本文对相邻波粒子之间间距大于波粒子半径1／2的粒子进行分裂处理，分裂出2个新粒子，分布在其左右两边。而粒子一旦产生，其分裂时间可以通过公式进行计算，如下所示

因此，水波传播过程中，不需要计算粒子之间的间隔，避免了粒子之间的关联，从而将复杂问题转化成了简单的时间检测。除此之外，本文基于粒子的二维粒子系统，极大地减少了检测次数，达到了高效的检测效率。

在波粒子分裂为3个波粒子的过程中，新粒子振幅为分裂粒子振幅的1／3，同时，扩散角缩小为原扩散角的1／3，如图4所示。即保证了波浪的波阵面的平滑连续，也模拟了现实中水波振幅随着传播而降低的现象。

3.4 交互水面绘制

交互水面的绘制则结合波粒子系统对水面的形态采用GLSL完成最终的绘制。本文设计了4个GPU自定义绘制管线来完成，即4对GLSL。第一个绘制管线根据CPU传入的粒子信息创建粒子信息图，如图5所示，该图存储粒子的振幅及传播方向信息，并传入第2个绘制管线中；在第二个绘制管线中，结合第1个管线产生的粒子信息图，采用式 （2）对信息图中的振幅在X轴方向上进行过滤；在第3个绘制管线中，则对信息图在Y轴方向上进行过滤；第4个绘制管线则结合两个方向上过滤出来的高度图，结合式 （3）进行圆形平滑过滤，使得高度场平滑连续，满足波形特征。

图5 粒子信息

4 实验结果及分析

本文采用基于母粒子的波浪粒子系统，并且利用GPU进行波浪粒子属性的更新，取得了非常高效的绘制效率，如表1和图6所示。表1是一个帧率表，其中，行表头表示水波数量，列表头表示水波同时分裂的次数，为了便于统计和控制，测试水波都是同时产生的，因此分裂也是同时进行。本文实验采用的屏幕渲染尺寸为1440×900，水面网格尺寸为1024×1024。

表1 不同波浪数量和分裂次数下的运行帧率／（ftp／s）

效果图如图7和图8所示。从图中可以看出，本文算法实现的水面具有较高的真实感，且水波波形也具有非常高的逼真度。

5 结束语

本文设计了一个基于GPU的实时可交互模拟算法，对水面高度场进行离散化处理，基于现实中水波的组成，设计了局部垂直和局部水平偏移函数，进行水波波形的模拟。同时，通过离散化处理后，结合粒子系统，实现了水波的动态实时更新，而基于母粒子的波粒子系统的提出和应用以及波粒子属性的合理设计，极大地降低了遍历开销，并将粒子的更新全部放到GPU上进行，实现了高效地更新计算。为了提高水面的真实感，本文结合水的光学特性，基于水面光照模型，对水面进行了光照计算，并采用立方体贴图技术模拟水面对周围环境的反射效果。实验结果表明，本文算法实现了高效稳定的可交互水波实时模拟。

图7 鼠标滑动触发的水波效果

图8 小船在水面运动导致的水波

［1］Kallin D.Real－time large scale fluids for games ［C］／／SIGRAD，2008：31－35.

［2］Lee H，Han S.Solving the shallow water equations using 2D SPH particles for interactive applications ［J］.The Visual Computer，2010，26 （6－8）：865－872.

［3］Yuksel C，House D H，Keyser J.Wave particles ［J］.ACM Transactions on Graphics（TOG），2007，26 （3）：991－998.

［4］Cords H.Moving with the flow：Wave particles in flowing liquids［C］／／International Conference in Central Europe on Computer Raphics，Visualization and Computer Vision，2008：145－152.

［5］Craig Schroeder.Semi－implicit surface tension formulation with a Lagrangian surface mesh on an Eulerian simulation grid ［J］.Journal of Computational Physics，2012，231 （4）：2092－2095.

［6］Andrew Selle，Ronald Fedkiw，ByungMoon Kim，et al.An unconditionally stable MacCormack method ［J］.Journal of Scientific Computing，2008，2008 （35）：350－371.

［7］Zhang Yizhong，Wang Huamin，Wang Shuai，et al.A deformable surface model for real－time water drop animation ［J］.IEEE Transactions on Visualization and Computer Graphics，2012，18 （8）：1281－1289.

［8］Chentanez N，Feldman B E，Labelle F，et al.Liquid simulation on lattice－based tetrahedral meshes ［C］／／Proceedings of the ACM SIGGRAPH／Eurographics Symposium on Computer Animation.Eurographics Association，2007：219－228.

［9］Nipun Kwatra，Chris Wojtan，Mark Carlson，et al.Fluid simulation with articulated bodies ［J］.IEEE Transactions on Visualization and Computer Graphics，2010，16 （1）：70－80.

［10］Becker M，Teschner M.Weakly compressible SPH for free surface flows ［C］／／Proceedings of the ACM SIGGRAPH／Eurographics Symposium on Computer Animation.Eurographics Association，2007：209－217.

［11］Solenthaler B，Pajarola R.Predictive－corrective incompressible SPH ［C］.ACM Transactions on Graphics.ACM，2009，28 （3）：401－406.

［12］NVIDIA.Physx page ［EB／OL］.［2010－12－15］.http：／／www.nvidia.com／object／physx＿new.html.