水体方案综述 - 波形篇

作者: 离原春草 | 来源:发表于2024-07-26 09:12 被阅读0次

波形模拟解决的是顶点在每一帧的tick中如何更新自己的transform的问题，按照[1]的总结，目前有多种模拟方案（相信这里也不是全部的方案）：
1.线性波形叠加方案

a.正弦波（Sinusoids Wave）[Max 1981]
b.Gerstner波（Gerstner Wave）[Fournier 1986]

2.统计模型方案

a.快速傅立叶变换（Fast Fourier Transform，简称FFT）[Mastin 1987]
b.空间-频谱混合方法（Spatial -Spectral Approaches） [Thon 2000]

3.波动粒子方案

a.波动粒子（Wave Particles）[Yuksel 2007]
b.水波小包方法（Water Wave Packets）[Jeschke 2017]
c.水面小波方法（（Water Surface Wavelets）[Jeschke 2018]

4.物理方程方案

a.欧拉方法（Eulerian approaches）[Kass 1990]
b.拉格朗日方法（Lagrangian approaches）[Stam 1995]
c.欧拉-拉格朗日混合方法（Hybrid approaches）[Brien 1995]

波动方程（Water equation）方案

a. Water equation
b. Shallow water equation

1. 线性波形叠加算法

这里所谓的线性波形叠加指的不是线性波形的叠加，而是波形的线性叠加，即通过类似如下的公式实现的波形模拟：

$\sum_{i=1}^n A_i * Wave_i$

水波模拟常用的波形函数有两种，分别是正弦波(Sinusoids Wave）和Gerstner波(Gerstner Wave)。

1.1 正弦波

Max在1981年提出了将不同振幅的正弦函数累加起来以模拟水波的想法，对应到图形学，我们其实需要的是mesh上的每个顶点的高度y与其平面坐标(x, z)以及时间t之间的关系，用公式表示的话，大概是如下的形式：
$y=h(x,z,t)=y_0+\sum_{i=1}^n A_i sin(k_{i_x}x+k_{i_z}z+\omega_i t)$

其中：
● $y_0$ 可以看成是水波的原始高度（海平面）
●n是叠加波形的数目
● $A_i$ 是每道正弦波的振幅，即最大的起伏落差
● $k_{i_x}, k_{i_z}$ 则是正弦波分别在x、z两个方向上的（空间）频率（即波长的倒数乘以 $2*\pi$ ）
● $\omega_i$ 则是正弦波的时间频率

关于波形公式的介绍，可以参考Rendering Water With Sine Waves视频中的讲解，里面做了从零到一的介绍。

按照这个公式，我们可以看到，顶点并不会跟随波形而产生水平方向上的移动，只是会上上下下做简谐振动（现实中，水体粒子其实并不只是上上下下，而是一半时间在跟随波形往前，另一半时间在逆着波形往后，这一点在后续的Gerstner波模拟得比较好）。

正弦波的特点是平滑，圆润，适合表达如池塘一样平静的水面，下图是Unity下实现的基于正弦波的水体效果

1.2 Gerstner Wave

Gerstner波(Gerstner Wave，以发现者名字命名)也称为Trochoidal Wave（余摆波，以形状命名），是流体动力学中周期性表面重力波（periodic surface gravity waves）欧拉方程的精确解，由Gerstner在1802 年初次发现，并在1986年由Fournier等人引入水体渲染领域。

Gerstner波的特点是波峰尖锐，波谷平缓，在波峰两侧有收紧的趋势，与真实海洋表面比较接近，且计算量可控，性价比高。前面提到的正弦波虽然计算简单，但是其形状并不是很能匹配因为强风吹拂而形成的真实波形效果，模拟接近度比较高的时Stokes波形函数，但是这种函数计算过于复杂，在图形学中，常常用Gerstner波作为替代方案。

正弦波 vs Gerstner波，前者上上下下，后者在给定半径的圆上匀速运行

为了方便理解，我们先看2D平面下的顶点运动计算公式：
$x = x_0 + A * cos(k*x_0 + \omega * t) \\ y = A * sin(k*x_0 + \omega * t)$

其中:
● $x_0$ 是顶点在X轴上的原始位置
●A是波形振幅
●k是波形的（空间）频率（即波长的倒数乘以 $2*\pi$ ）
● $\omega_i$ 则是正弦波的时间频率

经由上述公式，我们得到的波形效果如下图所示：

用顶点来模拟，可以看到，每个顶点都是在做圆周运动（每个顶点的圆周周长不同）：

上面的公式给的是单个正弦波的计算逻辑，实际上，我们会叠加多个不同波长、不同振幅的波来实现更为真实的模拟，此外，在实际的水波模拟中，我们也需要将公式从2D扩展到3D：
$\vec x = \vec x_0 + \sum_{i=1}^n \frac{\vec k_i}{|\vec k_i|} A_i * cos(\vec k_i*\vec x_0 + \omega_i * t + \phi_i) \\ y = \sum_{i=1}^n \frac{\vec k_i}{|\vec k_i|} A_i * sin(\vec k_i*\vec x_0 + \omega_i * t + \phi_i)$

其中，上方带有箭头的都是2D向量，里面的两个分量分别是x & z方向上的坐标。

将这个公式换算成更有物理意义的格式：
$P(x,y,t)= \begin{cases} x+\sum_i Q_i * A_i * \vec{D_i}.x * cos(w_i \vec{D_i}+\phi_i t) \\ y+\sum_i Q_i * A_i * \vec{D_i}.y * cos(w_i\vec{D_i}+\phi_i t) \\ \sum_i A_i * sin(w_iD_i \cdot (x, y)+\phi_i t) \end{cases}$

这个公式代表的是水平面上(x, y)点在t时刻的三维坐标，这个公式相对于上面的公式增加了一些变量，这里借用[2]中的解释来做一下阐述：

波长( $L$ ):波与波峰之间的波峰距离,
频率( $w$ ): $w$ 与波长 $L$ 有关， $w=\frac{2}{L}$ 。
振幅 ( $A$ ):从水平面到波峰的高度。
速度( $S$ ):波峰每秒向前移动的距离。将速度表示为相位常数很方便，其中 $\phi=\frac{2S}{L}$ 。
方向 ( $\vec D$ ):垂直于波峰沿波峰行进的波前的水平矢量。
波浪陡度( $Q$ ): $Q$ 如果为0，则产生的是通常的正弦波；而如果 $Q =\frac{1}{wA}$ ，则生成的是尖锐的波峰；如果继续提高 $Q$ ，则会导致波峰形成环形。
时间变量(t):不断变化的时间变量，用以产生波移动的效果。

Gerstner Wave的优点有：

不需要进行贴图的采样，即可计算出各个顶点在任意时刻的偏移
计算简单
可以很容易就能实现效果在多端的同步

Gerstner Wave的不足有：

给定X-Z坐标，想要查询到该点的水面高度，需要通过数值迭代的方式来实现
不支持波形与场景人物的交互效果

目前应用了该波形方案的项目有：

《Shoestring》的Wii平台
Unity的Crest Ocean方案

2. 统计模型方法

2.1 快速傅立叶变换（Fast Fourier Transform , FFT）

FFT是目前电影业广泛采用的海洋表面渲染解决方案，这个方法的特点是真实感出色，全局可控，细节丰富，但计算量也相对较大。

傅里叶变换是一种用于信号在时域和频域之间的线性积分变换（时域转频域称为傅里叶变换，频域转时域则是傅里叶逆变换），而快速傅里叶变换则是一种可以在O（nlogn）时间内完成离散傅里叶变换（Discrete Fourier transform，DFT）的高效、快速计算方法集的统称。

水体的波动是一种时域的效果展示（即随着时间而波动），假设我们可以获取到水波的频域信息（频谱），并将之分解为若干正余弦（实际上是复指数，不过两者都是同一种事物的不同表达形式而已，不必纠结）信号，就能够反过来还原出（逼近）时域信号，这就是FFT方案的基本原理：
1.通过理论推导或实际测量，可以得到一张用于描述海洋表面的海浪频谱
a.最常见的频谱为[Tessendorf 2001]使用的Phillips频谱
2.该频谱可以转化为若干离散的正余弦信号（复指数信号）
3.对上述信号执行逆FFT，就可以在计算机中将结果从频域转换到空间域
4.经过运算生成位移贴图（displacement map），即海洋的高度场
5.基于位移贴图，我们可以导出表面法线贴图，以及如代表浮沫区域的Folding Map的其他相关数据

更多的实现细节，可以参考此前《【GDC 2015】An Introduction to Realistic Ocean Rendering through FFT》中的分析

疑问：如何实现水体的动态效果呢？需要多张频域贴图组成的序列帧吗？那如何实现水体与场景交互后的扰动呢？

采用FFT的水体渲染方法从90年代开始广泛用于电影制作（离线渲染）从2000年代开始广泛用于游戏（实时渲染）。
离线渲染和实时渲染的选择，主要在于当时硬件计算能力可以承受多少分辨率的高度图的实时运算。
早期的游戏，如Crysis，由于硬件计算量的限制，采用了64 x 64的高度图分辨率。而由于硬件的发展，目前512 x 512的分辨率的计算量已经在实时渲染中较为普遍。
电影工业中采用FFT生成的高度图，由于可以采用离线渲染，以及品质的要求，分辨率一般较大.
早在1997年的《泰坦尼克号》的海洋渲染的渲染，就已经采用了2048 x 2048分辨率的高度图。

FFT方案的优点有：

波形模拟的效果非常真实（电影界使用的方案）
可以很容易就能实现效果在多端的同步

FFT方案的不足有：

计算逻辑相对复杂
不支持波形与场景人物的交互效果
缺少边界效果

目前应用了该波形方案的项目有：

《Shoestring》的XBox等主机平台
《刺客信条III》
《盗贼之海》
《War Thunder》
电影《泰坦尼克号》 & 《海洋奇缘》

3. 波动粒子方法

3.1 波动粒子（Wave Particle）

波动粒子方法指的是将粒子作为水体的最小组成单元，通过对粒子的驱动来实现水体的波动模拟的方法，这种方法将动态模拟三维水波的复杂度降维到模拟在平面上运动的粒子系统的级别。

这个方法最初由Yuksel于2007年[Yuksel 2007]提出，其特点兼具了波形叠加方法的灵活性和基于FFT方法的稳定性和视觉细节，具有较好的可控制性、性能，以及出色的交互表现

参考[3]，对Wave Particle的基本原理做一下描述。

首先，假设我们将水体看成是一个连续的高度场，水平面上某一点的坐标为 $\vec p = (x, y)$ ，该点的高度为 $z(x, y)$ ，那么当外力（比如某个与水面发生交互的物件对水面施加的作用力）与水面发生交互，把水体看成是多个细分的单元组成的，那么就能很自然的想到，这个交互会在水面上的各个（单元的）位置处形成对应的高度偏差（deviations），这些偏差最后转化以v（没有方向？）为扩散速度的水波，并且这个水波将满足如下的二维也是二阶的波动方程：
$\frac{\partial^2 z}{\partial x^2} + \frac{\partial^2 z}{\partial y^2} = \frac{1}{v^2} \frac{\partial^2 z}{\partial t^2}$

为了得到上述方程的解析解，我们这里先做如下的格式转换：
$z(\vec p, t) = z_0 + \eta_z(\vec p, t)$

即我们将每个单元的高度表达为原始高度 $z_0$ 跟该点的高度偏差 $\eta_z(\vec p, t)$ 之和。

而该点的高度偏差又可以表示为从各个方向扩散过来的水波（局部）高度偏差之和：
$\eta_z(\vec p, t) = \sum_i D_i(\vec p, t)$

而通过将局部偏差与一个以速度 $v$ 进行扩散的wave particle（水波粒子）关联起来，每一道水波的高度偏差又可以表达成如下的局部偏差公式：
$D_i(\vec p, t) = A_i w_i(\vec p - \vec p_i(t))$

需要注意，这里的粒子，不是一个细粒度的粒子，可以理解为一个基本的水波形状。这个公式中 $A_i$ 是该水波的振幅， $w_i$ 是该水波的波形函数，输入的自变量是点 $\vec p$ 跟扩散的水波粒子 $\vec p_i$ 在时刻t的有向距离：

对于每个水波粒子来说，其偏差公式是不会发生变化的
水波粒子的位置、扩散方向以及其他的一些属性，共同决定了水波的形状

通过这种方法，我们就可以将复杂的流体模拟计算转化为了平面上粒子的移动计算，从而大大实现了计算的简化。

接下来我们需要的就是找到一个符合需要的局部偏差公式，或者叫波形函数，这个函数需要满足几点要求：

能够在水面上构造波阵面（或者说通俗一点，是实现水面高度场的起伏，生成水波波形）
构造的波形要足够的逼真，且能在水面上做整体的移动

在2维空间中，很容易就让人想到正弦波，基于正弦波，我们可以构造一个波形函数：
$w_i(u) = \frac{1}{2}(cos(\frac{2\pi u}{l_i})+1)\prod \frac{u}{l_i}$
其中 $l_i$ 是波长，而 $\prod$ 这里是矩形函数（rectangle function）：
$\prod(x)= \begin{cases} 1& {|x|<\frac{1}{2}}\\ \frac{1}{2}& {|x|=\frac{1}{2}}\\ 0&{otherwise} \end{cases}$

这个函数除了波形漂亮之外，还有如下的几个优点：

在一个有限区域内的数据是非0的（波形有起伏），且在边界区域有0值的一阶导数，从而可以平滑跟安静的水面衔接起来
基于上面的特性，我们可以通过调整振幅的正负号，同时设置合理的间距来得到连续的波形效果

单个水波粒子在三维空间中会构造出连续的波阵面，与其通过一个波阵面来模拟各种复杂的水波效果，这里的做法是通过多个波阵面的组合、混合来实现不同的效果，且这种做法可以保障前面列举的两个优点，这个混合可以通过如下的公式来表达：
$D_i(\vec p, t) = A_i w_i(u) B_i(v)$
其中
$u = \hat u_i \cdot (\vec p - \vec p_i(t))$
$v = \hat u_i^{\perp} \cdot (\vec p - \vec p_i(t))$

这里 $\hat u_i$ 跟 $\hat u_i^{\perp}$ 都是水波粒子局部坐标系中的2D向量， $\hat u_i$ 是水波粒子的传播方向向量，而 $\hat u_i^{\perp}$ 则是与传播方向垂直的方向向量，这两个向量与位置差值向量的点乘结果就得到了 $u, v$

$B_i$ 混合或者说权重函数，这个函数的选择没有明确的约束，相对主观，只要能够满足如下要求即可:

取值是有限的，非无穷大的
（针对不同v值，对应的函数值的）累加值等于1

可以看到波形函数的自变量是u，也就是说这个函数对应的是传播方向上的波形起伏效果，下图给出了多个水波粒子

(a) 单个水波粒子组合效果 (b) 多个水波粒子合并效果

每个水波粒子中存储的属性，除了位置之外，还有振幅跟半径等信息。当相邻两个粒子的距离大于半径的一半时，为了保障效果的平滑，需要对其中的一个粒子进行拆分（subdivision），从1个变成3个，新增的3个粒子会基于能量守恒的考虑，其属性会与此前的粒子的属性存在较强关联：
1.新粒子的振幅和扩散角度变为父粒子的三分之一
2.新粒子的半径与父波动粒子保持相同

Yuksel提出的原版Wave Particles的波动粒子源是基于输入的point source，而《神秘海域3》中则是在环形区域中放置随机分布的粒子源（速度跟位置都随机），以实现对开放水域的混沌运动的近似模拟：

《神秘海域4》中则采用了多分辨率Wave Particles方案，高低频信号叠加，进一步改进了表现与性能

Wave Particle方案的优点有：

支持与场景的交互
具有很好的边界效果
模拟效果是稳定的（相对于Water equation以及shallow water equation来说）

Wave Particle方案的不足有：

不支持大尺寸的流动效果
网络同步可能会存在风险？

目前应用了该波形方案的项目有：

《神秘海域3》 & 《神秘海域4》
《Shoestring方案》

3.2 水波小包方法（Water Wave Packets）[SIGGRAPH 2017]

沿着波动粒子模拟的思路， Jeschke和Wojtan[2017]在SIGGRAPH 2017介绍了以theoretical group speed（理论群速度）传播的Water wave packets（水波小包）方法。

此方案在保留了基于频谱方案数值稳定和可得到理论准确的速度的特点，同时通过将全局余弦波分解成一系列更短的波分量以降低基于频谱的方法的复杂性

3.3 水面小波方法（Water Surface Wavelets）

在Siggraph 2018上，Jeschke接着提出了针对Water Wave Packets的改进方法：Water Surface Wavelets。

Water Surface Wavelets基于欧拉方法实现，自由度与空间区域有关，与波动本身无关。因此，该方法可以和GPU更好的结合，不过还没听说有项目尝试过

4. 基于物理的方法

基于物理的方法指的是采用类似物理模拟（流体模拟）的方式来实现水体表面的运动驱动的方法，这类方法一般消耗较高，多用于离线渲染等场景。

基于物理模型的水体模拟方法的核心是对Navier-Stokes方程（Navier-Stokes Equations,NS方程）进行求解。

Navier-Stokes方程是一组描述液体和空气之类的流体物质的方程，描述作用于液体任意给定区域的力的动态平衡。除了水体模拟之外，其还可以用于模拟天气，水流，气流，恒星的运动。

Navier-Stokes方程如下：

$\frac{\partial u}{\partial t} + u \cdot \nabla u = -\frac{\nabla p}{\rho}+g+v \nabla^2u \\ \nabla \cdot u = 0$

其中：

● u为流体速度矢量（fluid velocity vector）

● p为流体压力（fluid pressure）

● $\rho$ 为流体密度（fluid density）

● v为运动粘度系数（kinematic viscosity coefficient）

● $\nabla$ 为梯度微分（gradient differential）

● $\nabla^2$ 拉普拉斯算子（Laplacian operators）

求解Navier-Stokes方程常用的方法有欧拉方法（Eularian Method）和拉格朗日方法（Lagrangian Method）两种：

1.欧拉方法（Eularian Method）是一种基于网格的方法。
a.它从研究流体所占据的空间中各个固定点处的运动着手，分析被运动流体所充满的空间中每一个固定点上流体的速度、压强、密度等参数随时间的变化，以及由某一空间点转到另一空间点时这些参数的变化。

2.拉格朗日方法（Lagrangian Method）是一种基于粒子的方法。
a.它从分析流体各个微粒的运动着手，即研究流体中某一指定微粒的速度、压强、密度等参数随时间的变化，以及研究由一个流体微粒转到其他流体微粒时参数的变化，以此来研究整个流体的运动。
b.最常用的拉格朗日方法是光滑粒子流体力学(Smoothed Particle Hydrodynamics，SPH)方法，其核心渲染思想为流体模拟产生粒子，然后多边形化粒子以产生波。

基于SPH方法的水体渲染表现

除了独立的两种方法之外，还有结合两者的欧拉-拉格朗日混合方法（Eularian-Lagrangian Hybrid approaches），其主要思想是使用欧拉方法来模拟流体的主体，并使用拉格朗日方法来模拟诸如泡沫，喷雾或气泡之类的细小细节

另外，也可以采用bake to flipbook方法，将离线的流体模拟，烘焙成flipbook帧动画，用于实时渲染。

Houdini下的流体模拟

5. 波动方程方案

5.1 Water euqation

基于波动方程的方法是基于这样一个假设提出的：水体表面是一个高度场，水体的波动永远维持在垂直方向，水体较浅，只支持小尺寸的波动，对水体的交互压力也局限于垂直方向，同时放弃非线性的部分。

这种方案也是通过2D数组来存储高度数据的

这个方案的基本思路是：

用2D数组存储水体高度数据
相邻数据节点（顶点？）通过pipe（管道）来连接
pipe内部的流速则是基于高度来更新
之后又基于流速来完成对顶点高度的更新
这里应该还要补一点，交互产生的外力，会用于更新流速或者顶点的高度，以实现上述迭代的启动

这个方法的优点是：

计算很快
支持与场景的交互
具有边界效果

方案的不足是：

不支持旋涡效果
不支持大尺寸的流动效果
并不是永远稳定的

5.2 Shallow Water Equation

水方程方案是波动方程的子集，与波动方程的假设的差异是：

需要考虑非线性部分
不存储高度，而是存储速度

这个方案相对于波动方程方案来说，差异是支持了旋涡效果，依然不是永远稳定。

参考

[1]. 真实感水体渲染技术总结
 [2]. UE5入门学习笔记（一）——Gerstner波
 [3]. 【Siggraph 2017】Wave Particles（波粒子）
[4]. Wave Particles - Siggraph 2017 pdf

水体方案综述 - 波形篇

1. 线性波形叠加算法

1.1 正弦波

1.2 Gerstner Wave

2. 统计模型方法

2.1 快速傅立叶变换（Fast Fourier Transform , FFT）

3. 波动粒子方法

3.1 波动粒子（Wave Particle）

3.2 水波小包方法（Water Wave Packets）[SIGGRAPH 2017]

3.3 水面小波方法（Water Surface Wavelets）

4. 基于物理的方法

5. 波动方程方案

5.1 Water euqation

5.2 Shallow Water Equation

参考

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

水体渲染