虚幻引擎 Gerstner Waves -GPU Gems 从物理模型中实现有效的水体模拟

1.1 目标与范围

我们从简单的正弦函数开始，然后逐步过渡到更复杂的函数，以适应需要。
本章主要解释系统参数的物理意义，表明将水表面近似为正弦波的总和并不像人们通常认为的那样是随意的。我们特别关注将基本模型转换为实际实现所需的数学方法，因为这些数学方法是扩展实现的关键。
该系统适用于从小池塘到海洋的各种水体，从海湾或岛屿的角度进行观察。尽管这不是一个严格的物理模拟，但它确实能够提供令人信服、灵活且动态的水体渲染效果。由于该模拟完全在GPU上运行，因此不会与AI或物理争夺CPU资源。
由于系统参数具有一定的物理基础，因此它们比通过反复试验得出的参数更容易进行脚本编写。
将整个系统设置为动态的——除了其组件波浪之外——为系统增添了额外的生命力。

1.2 "正弦和"近似法

我们运行了两个表面模拟：一个是表面网格的几何起伏模拟，另一个是该网格上法线贴图中的涟漪模拟。这两个模拟本质上是相同的。水面的高度由简单周期波的总和来表示。我们从叠加正弦函数开始，随着进程发展，逐渐过渡到更有趣的波形。
正弦之和提供了一个连续函数，用于描述所有点上水面的高度和表面方向。在处理顶点时，我们根据每个顶点的水平位置对该函数进行采样，使网格符合其细分到连续水面的限制。在几何分辨率以下，我们将这种技术继续应用到纹理空间中。我们通过在渲染到渲染目标纹理的过程中，通过简单的像素着色器操作对正弦之和的近似法线进行采样，生成表面的法线贴图。每一帧渲染我们的法线贴图，让我们有限的正弦波集可以独立移动，极大地增强了渲染的真实感。
实际上，我们水体纹理中的细微波纹在很大程度上决定了我们模拟的真实感。我们波浪表面的几何波动为我们提供了一个更微妙的框架，用于呈现这种纹理。因此，我们对选择几何波浪还是纹理波浪有不同的标准。

1.2.1选择波形

我们需要一组参数来定义每个波形。如图1-2所示，这些参数包括：

单个波函数的参数

图1-2 单个波函数的参数

Wavelength _波长（ $L$ ）

世界空间中波峰到波峰之间的距离。
$L$ 与频率 $w$ 之间的关系为 $\frac{2}{L}$ 。

Amplitude _振幅（ $A$ ）

从水平面到波峰的高度。

Speed _速度（S）

波峰每秒前进的距离。
将速度表示为相位常数 $\varphi$ 更为方便。
$\varphi = S \times \frac{2}{L}$

Direction _方向（D）

波前垂直方向上的水平向量，该向量指向波峰的运动方向。

然后，每个波的状态可以定义为水平位置 (x, y) 和时间 (t) 的函数：

Equation 1 _公式1

$W_{i}\left(x,y,t\right)= A_{i}\times\sin\left(\mathbf{D}_{i} \cdot\left(x,y\right)\times w_{i} + t\times\varphi_{i}\right)$
而整个表面是：

Equation 2 _公式2

$H\left(x,y,t\right)=\sum\left(A_{i}\times\sin\left(\mathbf{D}_{i} \cdot\left(x,y\right)\times w_{i} + t\times\varphi_{i}\right)\right)$

总和所有波的 $i$ 。

为了在场景的动态中提供变化，我们将在一定的约束条件下随机生成这些波浪参数。随着时间的推移，我们将持续地将一个波浪渐隐，然后再以一组不同的参数将其渐显。事实证明，这些参数是相互依赖的。必须小心地为每个波浪生成一整套参数，这些参数需要以一种令人信服的方式组合在一起。

1.2.2 法线和切线

因为我们有一个关于表面的显式函数，我们可以直接计算任何给定点的表面方向，而不需要依赖有限差分技术。我们的副法线向量 B 和切线向量 T 是分别沿 x 和 y 方向的偏导数。对于二维水平面中的任何点 (x, y)，表面上的三维位置 P 是：

Equation 3 _公式3

$\mathbf{P}\begin{pmatrix}x,y,t\end{pmatrix}=\begin{pmatrix}x,y,H\begin{pmatrix}x,y,t\end{pmatrix}\end{pmatrix}$
那么，x方向的偏导数是：

Equation 4a _公式4a

$\mathbf{B}\left(x,y\right)=\left({\frac{\partial x}{\partial x}},{\frac{\partial y}{\partial x}},{\frac{\partial}{\partial x}}\left(H\left(x,y,t\right)\right)\right)$

Equation 4b _公式4b

$\mathbf{B}\left(x,y\right)=\bigg(1,0,\frac{\partial}{\partial x}\big(H\left(x,y,t\right)\big)\bigg)$

同样地，切线向量为：

Equation 5a _公式5a

$\mathbf{T}\left(x,y\right)=\left({\frac{\partial x}{\partial y}},{\frac{\partial y}{\partial y}},{\frac{\partial}{\partial y}}\left(H\left(x,y,t\right)\right)\right)$

Equation 5b _公式5b

$\mathbf{T}\left(x,y\right)=\left(0,1,{\frac{\partial}{\partial y}}\left(H\left(x,y,t\right)\right)\right)$

法线由副法线和切线的叉积给出，表示为：

Equation 6a _公式6a

$\mathbf{N}(x,y)=\mathbf{B}(x,y)\times\mathbf{T}(x,y)$

Equation 6b _公式6b

$\mathbf{N}\left(x,y\right)=\left(-\frac{\partial}{\partial x}\Big(H\left(x,y,t\right)\Big),-\frac{\partial}{\partial y}\Big(H\left(x,y,t\right)\Big),1\right)$

在我们将函数 $H$ 的偏导数代入之前，请注意方程 3-6 中的公式是多么方便。两个偏导数的计算为我们提供了切空间基的九个分量。这是我们使用高度场来近似我们表面的直接后果。
即 $P (x, y) . x = x$ 和 $P (x, y) . y = y$ ，在偏导数中变成了0和1。这仅适用于这样的高度场，但对我们选择的任何函数 $H (x, y, t)$ 都是通用的。

对于 选择波形 章节中描述的高度函数，其偏导数特别容易计算。因为函数的导数是各部分导数之和：

Equation 7 _公式7

$\begin{aligned} {\frac{\partial}{\partial x}}\big(H(x,y,t)\big)& =\sum\biggl(\frac{\partial}{\partial x}\bigl(W_{i}\bigl(x,y,t\bigr)\bigr)\biggr) \\ &=\sum\left(w_{i}\times\mathbf{D}_{i}.x\times A_{i}\times\cos\left(\mathbf{D}_{i}\cdot(x,y)\times w_{i}+t\times\varphi_{i}\right)\right) \end{aligned}$

覆盖所有波浪 $i$ 。

关于通过直接叠加正弦波生成的波浪的一个常见抱怨是它们有太多的“滚动”，而真实的波浪拥有更尖锐的峰和更宽的谷。事实证明，有一个简单的正弦函数变体可以非常可控地产生这种效果。我们将我们的正弦函数偏移使其非负，并将其提高到指数 k。该函数及其相对于 x 的偏导数是：

Equation 8a _公式8a

$W_{i}\left(x,y,t\right)=2A_{i}\times\left(\frac{\sin\left(\mathbf{D}_{i}\cdot\left(x,y\right)\times w_{i}+t\times\varphi_{i}\right)+1}{2}\right)^{k}$

Equation 8b _公式8b

$\begin{aligned}\frac{\partial}{\partial x}\big(W_{i}\left(x,y,t\right)\big)&=k\times\mathbf{D}_{i}.x\times w_{i}\times A_{i}\times\left(\frac{\sin\left(\mathbf{D}_{i}\cdot\left(x,y\right)\times w_{i}+t\times\varphi_{i}\right)+1}{2}\right)^{k-1}\\&\quad\times\cos\left(\mathbf{D}_{i}\cdot\left(x,y\right)\times w_{i}+t\times\varphi_{i}\right).\end{aligned}$

图 1-3 展示了作为功率常数 k 的函数生成的波形。这是我们实际用于纹理波的函数，但为了简化，我们继续用我们简单的正弦和来表达波浪，并且我们注意到我们必须在何处考虑我们在基础波形中的变化。

图1-3

图1-3 不同波形的示例

1.2.3 几何波形^{GeometricWaves}

我们将自己限制在四种几何波形中。增加更多的波形并不涉及新概念，只是增加了更多相同的顶点着色器指令和常数。

方向性或圆形^{DirectionalOrCircular}

我们可以选择圆形或方向性波形，如图1-4所示。方向性波形需要的顶点着色器指令稍微少一些，但具体选择取决于正在模拟的场景。
在这里插入图片描述

图1-4 方向波与环形波

对于方向性波浪，方程1中的每一个 $D_i$ 在波浪的整个生命周期内都是常数。
对于圆形波浪，方向必须在每个顶点处计算，简单来说就是从波浪中心 $C_i$ 到顶点的规范化向量：

$\mathbf{D}_{i}\left(x,y\right)=\left({\frac{\left(x,y\right)-\mathbf{C}_{i}}{\left|\left(x,y\right)-\mathbf{C}_{i}\right|}}\right)$

对于大面积的水体，方向性波浪通常更受青睐，因为它们更好地模拟了风驱动的波浪。对于波浪来源不是风的小型水池（如瀑布底部），圆形波浪更为合适。圆形波浪还有一个好处是它们的干涉图案永不重复。

这两种类型的波浪的实现非常相似。
对于方向性波浪，波浪方向是从关于风向的某个范围内随机抽取的。
对于圆形波浪，波浪中心是从某个有限范围内随机抽取的（例如瀑布击中水面的线）。

接下来的讨论将重点放在方向性波浪上。

格斯特纳波^{GerstnerWaves}

为了有效的模拟，我们需要控制波浪的陡峭程度。如前所述，正弦波呈现出圆润的外观——这可能正是我们想要的平静、田园诗般的池塘效果。但对于粗糙的海面，我们需要形成更尖锐的波峰和更宽阔的波谷。我们可以使用公式8a和8b来实现所需的形状，但我们选择了相关的格斯特纳波。格斯特纳波函数最初是在计算机图形学出现很久之前就被开发出来，用于在物理基础上模拟海水。因此，格斯特纳波提供了一些表面运动的微妙变化，这些变化非常令人信服而不显突兀。（详细描述参见Tessendorf 2001）。
我们选择格斯特纳波，因为它们有一个常被忽视的特性：通过将顶点向每个波峰移动，它们形成了更尖锐的波峰。由于波峰是表面上最尖锐（即最高频）的特征，这正是我们希望顶点集中的地方，如图1-5所示。

GerstnerWaves

图1-5 GerstnerWaves

这是GerstnerWaves函数：

Equation 9 _公式9

$\mathbf{P}\left(x,y,t\right)=\left(\begin{array}{l}{{x+\sum\left(Q_{i}A_{i}\times\mathbf{D}_{i}.x\times\cos\left(w_{i}\mathbf{D}_{i} \cdot\left(x,y\right)+\varphi_{i}t\right)\right),}}\\{{y+\sum\left(Q_{i}A_{i}\times\mathbf{D}_{i}.y\times\cos\left(w_{i}\mathbf{D}_{i} \cdot\left(x,y\right)+\varphi_{i}t\right)\right),}}\\{{\sum\left(A_{i}\sin\left(w_{i}\mathbf{D}_{i} \cdot\left(x,y\right)+\varphi_{i}t\right)\right),}}\end{array}\right)$

这里 $Q_i$ 是一个控制波浪陡峭程度的参数。
对于单个波浪 $i$ ， $Q_i=0$ 产生常见的滚动正弦波，而 $Q_i = \frac{1}{w_i A_i}$ 产生尖锐的波峰。应避免使用较大的 $Q_i$ 值，因为它们会在波峰上方形成环。
实际上，我们可以将 $Q$ 作为“陡峭程度”参数留给制作艺术家来指定，允许范围是 $0$ 到 $1$ ，并使用 $Q_i = \frac{Q}{w_i A_i \times \text{numWaves}}$ 来变化，从完全平滑的波浪到我们能产生的最尖锐的波浪。

值得注意的是，公式3和公式9之间唯一的区别是顶点的侧向移动。高度是相同的。这意味着我们不再有一个严格的高度函数。即， $\mathbf{P}(x,y,t).x \neq x$ 然而，该函数仍然容易微分，并且有一些方便的项消去。为了省去推导过程，作为读者的练习，我们看到切向空间基向量是：

Equation 10 _公式10

$\mathbf{B}= \begin{pmatrix} 1-\sum\bigl(Q_{i}\times\mathbf{D}_{i}.x^{2}\times WA\times S()\bigr),\\ -\sum\bigl(Q_{i}\times\mathbf{D}_{i}.x\times\mathbf{D}_{i}.y\times WA\times S()\bigr),\\ \sum\bigl(\mathbf{D}_{i}.x\times WA\times C()\bigr) \end{pmatrix}$

Equation 11 _公式11

$\mathbf{T}= \begin{pmatrix} -\sum\bigl(Q_{i}\times\mathbf{D}_{i}.x\times\mathbf{D}_{i}.y\times WA\times S()\bigr),\\ 1-\sum\bigl(Q_{i}\times\mathbf{D}_{i}.y^{2}\times WA\times S()\bigr),\\ \sum\bigl(\mathbf{D}_{i}.y\times WA\times C()\bigr) \end{pmatrix}$

Equation 12 _公式12

$\mathbf{N}= \begin{pmatrix} -\sum\bigl(\mathbf{D}_{i}.x\times WA\times C()\bigr),\\ -\sum\bigl(\mathbf{D}_{i}.y\times WA\times C()\bigr),\\ \sum\bigl(Q_{i}\times WA\times S()\bigr) \end{pmatrix}$

说明:
$WA=w_{i}\times A_{i},$
$S()=\sin\left(w_{i}\times\mathbf{D}_{i}\cdot\mathbf{P}+\varphi_{i}t\right)$
$C()=\cos\bigl(w_{i}\times\mathbf{D}_{i}\cdot\mathbf{P}+\varphi_{i}t\bigr)$

这些公式不像4b、5b和6b方程那样简洁明了，但它们计算起来非常高效。

在形成波峰环的背景下，仔细观察法线的 $z$ 分量证明了这一点非常有趣。虽然 Tessendorf (2001) 从流体动力学的纳维-斯托克斯¹描述和“李变换技术²”中推导出他的“切碎效应³”，最终结果是在频率域中表达的格斯特纳波的一个变体。
在频率域中，可以避免并检测到波顶的环形，但在空间域中，我们可以清楚地看到正在发生的事情。
当 $Q_i \times w_i \times A_i$ 的和大于 $1$ 时，我们法线的 $z$ 分量在峰值处可能变为负值，因为我们的波浪会在自身上方形成环。
只要我们选择的 $Q_i$ 使得这个和始终小于或等于 $1$ ，我们将形成尖锐的峰值但永远不会形成环。

1.2.4参数解释

波长和速度 ^Wavelength ^Speed

我们首先选择合适的波长。与其追求现实世界的分布，不如最大化我们能承担的少数波浪的效果。
相似长度的波浪的叠加突出了水面的动态性。
因此，我们选择一个中值波长，并在该长度的一半到两倍之间生成随机波长。中值波长在创作过程中被编写，它可以随时间变化。
例如，在暴风雨期间，波浪可能比晴朗平静时显著更大。
注意，我们不能改变活跃波浪的波长。即便是逐渐改变，波浪的波峰也会向原点扩展或收缩，这看起来非常不自然。因此，我们改变当前的平均波长，随着时间的推移，当波浪逐渐消退，它们将基于新的长度重新生成。方向也是如此。

根据波长，我们可以轻松计算它在表面上的传播速度。水的色散关系（参见 Tessendorf 2001），忽略高阶项，给出：