基于波环粒子的实时水波仿真方法

顾浩杰，张 军，2*

（1.江南大学 人工智能与计算机学院，江苏 无锡 214122；2.江苏省媒体设计与软件技术重点实验室，江苏 无锡 214122）

0 引言

流体模拟一直是计算机图形学中极具挑战的研究问题，其计算复杂度远超现有通用计算机的实时计算能力。与此同时，由于自然场景的虚拟仿真中不可避免涉及流体对象，因此流体模拟吸引着包括互动游戏，媒体艺术以及增强现实（Augmented Reality，AR）与虚拟现实（Virtual Reality，VR）等领域的大量开发人员。

水波模拟是流体模拟研究中一个较为特殊的研究点，直接采用流体模拟通用方法求解纳维-斯托克斯方程（Navier-Stokes equations，NS），很难实现高逼真度的实时计算性能。伴随着游戏领域出现大量户外开放环境仿真的需求，水波的实时仿真和渲染技术越来越受到行业人员的重视。

自1980 年以来，计算机图形学领域的学者已经提出过多种模拟水波方法。这些方法的主要策略是采用微幅波假设［1］对水波问题进行线性化处理，以近似求解满足不可压缩条件的NS。根据具体模拟方法的不同，可以将这些方法细分为欧拉网格法［2-3］和拉格朗日粒子法［4-5］。

欧拉网格法计算开销相对较少，但会受CFL（Courant-Friedrichs-Lewy）条件限制。该方法将空间分割为离散的元胞，然后追踪流过元胞的被模拟材料，并将振幅等模拟信息记录在每个网格单元里。2015 年，Jeschke 等［6］提出一种高效的波前追踪算法模拟水波与复杂环境的交互效果，并利用多值函数插值的方法避开了奈奎斯特限制。这种模拟水波的方法需要进行大量的预计算，因此不适用于实时交互的可变场景。为了模拟出真实的频散效果，Canabal 等［7］在求解偏微分方程数值解的基础上使用频散核，表现了丰富的水波动态细节，并通过金字塔过滤器和阴影卷积运算有效提高了计算效率以及与障碍物碰撞现象的几何计算效率。Jeschke等［8］基于小波的水波离散化，使水面波的模拟既有丰富的水波细节，又可以高效处理碰撞关系。

相较于欧拉网格法，拉格朗日粒子法避免了网格分辨率的限制，可以在任意不规则区域实现水波仿真，并避免了无水波区域的“空计算”问题。拉格朗日方法将水波分割为离散、独立的元素，然后分别追踪其中每个元素的运动过程。每个离散的元素在内存空间中有一个独立的数据结构体，负责携带相关水波的模拟数据。2007 年，Yuksel 等［9-10］通过拉格朗日粒子法模拟水波的反射效果，并根据体积守恒的规则推算出物体入水与出水过程中水波在物体周围的位置与振幅。该方法渲染出的水波具有一定视觉逼真性，计算性能也达到了实时交互水平，是基于粒子系统实时模拟水波交互场景的经典算法。

Yuksel 等［9-10］的方法在水波物理特性的细节（比如频散现象）方面仿真精度非常有限，经不起仔细观察。对此，Jeschke 等［11］引入“波包”的概念，不仅能有效模拟出水波的高频波细节，而且使水波具有物体上更加精确的波速。最近，Skrivan 等［12］基于波包的概念，提出水波曲线方法，通过插值函数与多个控制点代表曲线的方式，成功地在三维动态流体表面上模拟出丰富的水波细节。虽然这种方法模拟出的水波更为细致，但其计算复杂度较大，无法达到实时绘制的要求。

拉格朗日粒子方法面临的主要问题是粒子数量在水波发展初期和后期的差异非常大，导致计算量不稳定。初始水波较为简单，少量粒子即可完成仿真需求；而后期的水波环半径急剧增大，且发生复杂碰撞后会产生更多新水波，必须使用大量粒子才能叠加产生视觉上无裂痕的仿真效果。随着大量新粒子被添加到仿真系统，计算量会显著提升，渲染帧率出现大幅下降，严重影响开发人员对帧率的控制能力。

本文将有效减少水波模拟中粒子数量作为解决上述问题的途径。而减少粒子数量的一个直接方法就是增大粒子的覆盖面积，以避免海量微小粒子堆积在一起（如图1 所示）。为此，本文提出以波环作为基础粒子的方法（如图2（b）所示），动态提高粒子覆盖面积，从而显著减少水波模拟过程中所需的粒子数量。为解决由此带来的水波反射时碰撞分析的复杂性问题，本文采用镜像波源的方法。该方法有效解决了波环粒子作为圆环无法模拟碰撞反射的问题，并同时解决了传统小粒子处理反射可能出现的水波断裂问题。实践表明，本文方法在保持与波包方法近似的逼真度的同时，大幅减少了粒子使用数量，并有效控制了粒子增长规模，能在低端PC 机、移动平台等低性能硬件上实现交互级水波仿真应用。

图1 微小粒子堆积现象Fig.1 Accumulation phenomenon of small particles

1 算法理论

由于不同频率、不同方向传播的水波彼此不相干，所以水面可认为是大小不等、方向与频率各异的水波的叠加组合——波包（如式（1））。相较于传统的采用邻域传播思想对水波进行动态造型的方法［13］，使用拉格朗日粒子法并引入波包的概念，往往可以使模拟的水波更接近真实的水波。波包是局限在空间的某有限范围区域内的波动，波包整体按群速度cg在空间移动，内部按相速度cp移动。

在小粒子上引用波包，需要在二维条件下进行水面高度的计算（如图2（a）的矩形块）。本文利用波环粒子表示波包的圆环结构（如图2（b）），有效地将水面高度的计算问题维持在一维条件下，能大幅度减少模拟时所需的粒子数量，降低计算开销。

如图2（a）所示，代表相同大小的水波时，所需矩形粒子的数量远多于所需波环粒子的数量。随着水波继续发展，其波环半径持续增大，周长也相应增加。为保持波包与波包相互重叠（按照香农采样定理，必须重叠一半才能保持视觉连续性），新的波包必须持续加入到场景中，大幅延长了计算时间；而本文的波环只需修正其内部参数，波环数量始终保持不变，计算量保持稳定。

图2 矩形粒子与波环粒子模拟水波对比Fig.2 Comparison of rectangular particle and wave annulus particle in simulating water waves

依据波环粒子作为基础粒子的方式，本文算法结构如图3 所示。

图3 本文算法结构Fig.3 Framework of the proposed algorithm

首先，水面网格点的高度计算可以简化成网格点上波环粒子高度的叠加，因此高度函数可以近似表达为：

其中：N是网格点上对应的波环粒子数量，aj是波环粒子的振幅，kj是波环粒子的代表波数，cg是波环粒子的半径扩张速度，cp是波环粒子的内部相速度。ϕ(x)代表核函数，因为波包需要一个核函数来描述其在空间中起局部作用，并且其傅里叶变换Φ(k)在波数坐标系中起局部作用。Φ(k)越宽，波包则越狭窄，一般使用高斯函数作为波包的核函数。此外，波环粒子的宽度l选择为3 倍波环粒子的代表波长λj=2π/kj。

其次，由于波环继承了波包的概念，水波的频散变得易于模拟。通常使用ω(k)来描述频散关系，相速度由此表示为ω/k，群速度则表示为∂ω/∂k。角频率ω与波数k的关系一般表示为：

其中：g是重力，σ是表面张力，ρ是水的密度，k是波数，h是水深。由于波环粒子无法直接反映各处的水深带来的影响，因此需利用波环上多点的水深取平均的方式来近似水深对波环的影响。

波的相速度cp是波的相位在空间中传递的速度，描述波包中某一单频波的相位移动速度，其计算公式为：

式（3）正确描述了当水深h很大（深水波）或者波数k很大（表面张力波）时，相速度与水深无关；其他情况，相速度受水深影响。

群速度cg是能量或信息顺着波动传播的速度，其计算公式为：

根据式（4），当不同波长的水波表现出不同的传播速度时，造成混合各种波长的波环渐渐地在空间中分离开来，从而实现频散的效果。

式（5）即为最终水面高度的计算公式。此外，为了便于计算，本文将局部坐标范围控制在[-1，1]。通过对不同频段的水波进行叠加，可以得到较为真实的水波，如图4、5 分别展示了多频段一维波函数叠加的效果与三维下多频段波环叠加的效果。

图4 多频段一维波函数叠加Fig.4 Multiband one-dimensional wave function stacking

图5 多频段波环叠加Fig.5 Multiband band wave annulus particles stacking

最后，在传统的小粒子模拟水波的方法中，面对水波的碰撞反射，往往让入射水波和反射水波遵守光学中的镜面反射定律［14］，并由此更新碰撞后的粒子行进方向与位置。相较于小粒子，本文的波环粒子是代表一个圆环上的水波，而非水波中的局部，所以无法直接通过更新传播方向进行处理。综合考虑性能和精确性后，本文采取镜像虚拟波环的形式处理水波反射现象。当发生波环与容器壁发生碰撞后，以碰撞边界为镜面，在镜像波源处产生相同大小的新波环粒子，来模拟碰撞后的反射效果。此外，由于频散模拟会产生大量的同波源的波环，因此波环的数量远远大于波源的数量。使用波环检测是否与障碍物边界发生碰撞，需要大量成本来计算距离dki(k=1，2，…，p，i=1，2，…，n)，其中k为波环数量，n为边界数量。为了减少计算成本，本文采用障碍物边界检测碰撞的方法。如图6 所示，当且仅当出现新波源时，需更新各边界的动态表；边界依据波环粒子的波源ID 查找应该发生碰撞的距离，并与波环粒子当前时刻的半径进行对比以检测碰撞，从而达到避免重复计算的目的。

图6 边界动态表的更新Fig.6 Update of boundary dynamic table

2 水波环粒子的具体实现

本文算法的计算流程如图7 所示，在CPU 上处理水波的行为状态，在GPU 上处理水面网格的高度计算以及最终的水面绘制。

图7 本文算法具体实现流程Fig.7 Implementations flow of the proposed algorithm

水波行为状态处理需要每一“时步”调整波环粒子的属性，以保证水波物理特性以及位置的正确。水面高度的计算则需要水面网格点根据不同波环粒子的波源、外环的半径以及波环宽度来判断自身在哪些波环内，并叠加高度得到最终改点高度，如图8 所示，其中，Oi是波源，ri是波环粒子外环半径，li是波环宽度。

图8 水面网格的高度计算Fig.8 Height calculation of water surface grid

2.1 水波频散效果模拟

因为水波物理属性的复杂性以及可视性，水波频散效果对于水波的模拟十分重要。水波波长的不同，会使其传播速度也不相同，因此同一波源所产生的各种不同频率波在向外传播时会发生分离的现象。为了使水波具有较为真实的频散效果，本文采用频散细分方法，为每个波环粒子分配一个具体的波数范围[kmin，kmax]，使波环粒子的代表波数kj=(kmax-kmin)/2。根据式（4）可知，群速度是关于波数k的函数，群中的一些波比其他的波传播得更快，会将数据包的一部分推到平均群速度之前；而较慢的波则会将数据包的另一部分拖到后面。因此当分别以cg(kmin)与cg(kmax)行进的虚拟距离相差大于0.3 倍的波环宽度时，生成一个同波源的新波环。两个波环的波数范围分别调整为[kmin，kj]与[kj，kmax]，调整后两个波环的代表波数不同，传播速度也随之发生改变，从而模拟出频散效果。

此外，当发生频散后波环粒子细分为两个新波环粒子的情况，遵循能量守恒定律，需要对波环粒子的振幅与能量进行重新计算。振幅与能量满足：

其中：E为波环粒子能量，a为波环粒子振幅，S为波环面积。

假设发生频散后，两个新波环陡度是一样，则可以得到两个新波环粒子的能量之比：

其中：Ei是对应新波环粒子的能量；li是对应新波环粒子的波环宽度；ki是对应新波环粒子的代表波数；r是波环粒子的外环半径，因为细分是对波数的划分，因此两个波环粒子的外环半径相同。

图9、10 分别展示了有无频散对水波模拟的影响。通过对比可以明显看出，图10 中不同波长的水波是以不同的传播速度行进的，这使模拟的水波更为真实。

图9 没有频散现象的单波源模拟Fig.9 Simulation of single wave source without dispersion

图10 具有频散现象的单波源模拟Fig.10 Simulation of single wave source with dispersion

2.2 水波环反射计算

在复杂环境中，水波与障碍物发生碰撞会发生反射。当发生碰撞时，为了体现障碍物的表面粗糙性，可以将碰撞部分的水波的频谱拉伸到高频。针对不同计算平台的性能差异，本文通过设定全局碰撞变量γ控制波环的反射次数。γ值越大，反射次数越多，但计算量也会随之增加。由于本文为减少重复计算，记录各障碍物边界到各波源距离，所以边界检测波环粒子时，利用波环上一“时步”的外环半径与当前“时步”的外环半径进行处理。碰撞处理的算法如下所示。

为验证频散导致波环粒子增多而不会增加波源数目，表1 对不同时间段两处水波产生至不断扩散与频散的过程进行部分记录，并由表2 给出所使用的边界动态表的并集。

表1 波环粒子的半径变化Tab.1 Radius variation of wave annulus particles

表2 边界动态表的并集Tab.2 Union of boundary dynamic tables

图11 展示了本文方法所模拟的碰撞效果，此外，也可以清楚看到不同的γ造成的帧数与效果上的差别，相较于γ=1，当γ=2 时，有明显的二次反射，但粒子数量也明显增多，导致帧数下降。

图11 本文方法的碰撞效果模拟Fig.11 Collision simulation of the proposed method

2.3 水波扩散效果模拟

水波向外扩散的过程中，波环粒子的面积会发生改变，根据能量守恒定律，波环的振幅需要随之改变；所以，每个“时步”需要使用式（6）重新计算每个波环粒子的振幅。

此外，为了增加水波模拟的真实性，本文在对水波扩散的模拟中，对波环能量引入阻尼项。虽然现实中存在各种因素会使水波损失能量，但很多因素难以在较小的计算成本下模拟出；因此，本文主要考虑粘性势流、表面污染以及水波破碎所引起的能量耗散。粘性势流与表面污染可以通过对能量的衰减进行表示［15-16］，公式为：

其中：ν是水的粘度，一般约为10-6m2/s；a是波环的振幅；ξ≥0 为阻尼系数，用于调整数值耗散的程度。

水波的稳定性的一个重要参考指标是波陡，当波陡达一定限度，波峰就会破碎。由于无法直接模拟出水波的破碎的效果，因此通过对波陡大于0.142［17］的水波移除其一部分能量的方式进行近似的模拟。

为验证水波扩散过程中阻尼带来的影响，通过追踪单个波环粒子的振幅变化，从表3 可以看出振幅衰减速度的不同。为体现在视觉效果上的差别，图12 和图13 分别展示无阻尼和使用本文的阻尼项所模拟的水波扩散效果。通过对比可以看到，增加阻尼后，有效让振幅得到了渐进的指数衰减，并且让短波相较于长波得到较快衰减，使水波在扩散的过程中更为真实。

表3 不同阻尼下单粒子的振幅对比Tab.3 Comparison of amplitudes of single particles with different damping

图12 无阻尼的水波扩散效果Fig.12 Wave diffusion without damping

图13 本文方法的水波扩散效果Fig.13 Wave diffusion of proposed method

3 实验与结果分析

本文实验平台为普通笔记本计算机，主要硬件配置为Intel Core i7-6700HQ，显卡为GeForce GTX 1060M。

为体现水波模拟的效果与实时性，本文将实验结果与主流的文献［11］方法进行横向对比，其中本文的全局碰撞变量γ设定为3，阻尼系数ξ设定为1。

通过对粒子边界的绘制（如图14）可以看出，在模拟相同大小的水波时，本文所需的波环粒子数量远远小于文献［11］所需的矩形粒子数量。

图14 水波模拟所需粒子数对比Fig.14 Comparison of number of particles required for water wave simulation

图15 展示了本文方法和文献［11］方法水波碰撞反射后的扩散效果。可以从图中长方形框中看到，两种方法都能产生较为清晰的二次反射。此外，由于文献［11］中采用矩形粒子，所以在发生碰撞反射后，改变方向后已细分的矩形粒子无法快速再次左右细分，会产生如图15（b）椭圆形框中水波断裂的现象；而本文是以镜像波源产生新波环粒子的方式处理碰撞反射，因此避免了水波的断裂。

图15 水波碰撞反射效果对比Fig.15 Comparison of water wave collision reflection effect

图16 展示了文献［18］、本文以及文献［11］对雨滴涟漪的模拟效果，可以观察到，本文方法与文献［11］方法相较于文献［18］方法具有更为真实的频散特性。

图16 雨滴涟漪模拟效果对比Fig.16 Comparison of raindrop ripple simulation effect

为了对比本文方法和文献［11］方法在更新粒子数据上的成本，本文追踪在复杂环境下单波源水波产生到消逝的过程，记录其在CPU 端处理水波行为以及各种状态更新所需的计算成本。由于模拟计算时波包和波环粒子数量是动态变化的，造成CPU 端消耗时间波动较大，对消耗时间采用滑动平均，再取对数的方式进行对比分析。如图17 所示，可以明显看出，相较于文献［11］方法，本文方法在这方面的计算成本一直维持在2 ms 以下，对粒子数据更新的开销低且稳定。

图17 单波源水波CPU 端计算时间对比Fig.17 Comparison of calculation time of single wave source on CPU

在上述追踪在复杂环境下单波源水波产生到消逝的过程中，还同时记录了模拟时的帧率差别，图18 展示了本文方法和文献［11］方法模拟时的帧率。可以发现，本文方法帧数始终高于文献［11］方法，且维持在60 帧以上。

图18 复杂环境下单波模拟帧率对比Fig.18 Comparison of single wave simulation FPS in complex environment

为了进一步定量分析本文方法的实时性，表4 给出在复杂环境下同时产生多波源水波的情况进行追踪的时间统计结果。通过对比可以发现，本文方法在同时有多个波源水波产生至消逝的过程中，能保持流畅稳定，而文献［11］方法则会在中间时段出现明显的卡顿现象。

表4 复杂环境下多波源模拟帧率对比Tab.4 Comparison of multi-wave source simulation FPS in complex environment

最后，图19 给出本文方法实际水池场景中的具体应用效果。

图19 水池场景中多波源水波模拟Fig.19 Simulation of multi-wave source water wave in pool scene

4 结语

本文利用水波动态过程可以近似为中心平静且呈环形扩展的特点，将波环作为基础粒子用于水波实时仿真模拟。该方法继承了波包算法的物体特性逼真优点，并显著降低了水波模拟过程中对粒子数据更新的计算成本。通过使用镜像波源的方法，新波环算法可以正确模拟出水波碰撞后的反射波，有效避免复杂碰撞检测步骤。通过大量对比实验表明，本文方法的渲染结果在保持足够水波动态细节的情况下，使用更少的粒子数量，且不会出现波纹断裂现象。最终，本文提出的方法在普通PC 上可以达到更高的渲染帧数，从而更加适合移动、廉价硬件平台部署虚拟仿真应用。

由于波环在碰撞边界后会继续向外传播，故本文方法不能对水面漂浮的小型障碍物实现水波的阻挡效果，这是算法今后的改进方向。