Banner

FFT 海洋模拟：从频谱推导到 GPU 实时渲染

发布于 2026-01-20 · 更新于2026-05-19 12:44阅读
Unity URPOceanFFTCompute ShaderRendering

1. 引言

实时海洋渲染是 3A 游戏中最具挑战性的视效课题之一。当我们谈论”逼真的海面”时，真正令人信服的并不是单纯的高光与折射，而是那种来自统计学随机性的涌浪律动——长周期的大涌、中频的风浪、高频的毛刺细波，它们在频域中各自占据不同的能量区间，叠加后形成肉眼难以分辨规律的混沌美感。

为什么不用 Gerstner Waves？

Gerstner 波是正弦波的参数化变体，其波峰尖锐、波谷平缓，视觉上已非常接近真实海浪。然而它的本质是有限个确定性波形的叠加，参数越多越难手调，且难以产生真正意义上的随机海面纹理。即使叠加数十组波参数，也会在远处产生肉眼可辨的”重复模式”。

FFT 海洋的核心思想则完全不同：先在频域（Frequency Domain）用物理驱动的海浪频谱定义每个空间频率分量的能量，再通过逆快速傅里叶变换（IFFT）一次性将数万个波分量转回空间域，等效于同时叠加了成千上万组 Gerstner 波——而这在 GPU 上只需数毫秒。

方法 实质 GPU 开销 随机性 可平铺
正弦波叠加 少量确定波 极低 ❌ ❌
Gerstner Waves 参数化正弦 低 △ △
FFT 海洋 频谱统计模型 中等 ✅ ✅
SWE 流体 偏微分方程 高 △ ❌

本文定位

本文的重心在于波形模拟——从频谱物理到 IFFT 计算，以及 Jacobian 驱动的白沫模拟。渲染着色部分（PBR 水面材质、折射等）将仅作概要性介绍。工程实现参考 Ocean-URP（Unity 2022 LTS + URP）及 GodotOceanWaves 两个开源项目。

2. 数学与物理基础

海洋模拟归根结底是一道物理建模题：我们要用什么样的数学结构，才能捕捉到真实海面那种”无规律的规律感”？ 本节从物理直觉出发，逐步建立 FFT 海洋模拟的完整数学框架。

2.1 海浪的统计学视角：混沌中的秩序

想象你站在远洋货轮的甲板上俯视海面。你会发现一件奇妙的事：眼前的波浪从不重复，每一帧都与上一帧不同——但它又绝非完全随机的噪点。长涌从远方慢慢涌来，风浪在近处此起彼伏，偶尔一个大浪把前者掀得支离破碎。这种”有结构的随机性”，正是海洋物理的核心特征。

海洋学家的解答是：开阔洋面的风浪可以建模为一个平稳高斯随机过程。这句话拆开2个问题来理解：

高斯随机过程意味着什么？

海面上任意一点的高度 服从正态分布，均值为零（海面高低对称），方差有限（不会无限高）。更关键的是：不同空间频率（不同”波长”）的波分量之间统计独立——长涌的振幅不影响风浪的振幅，就像音响均衡器里每个频段独立调节一样。
这个假设有其物理依据：在开阔深水中，不同频率的波满足线性叠加原理，互相之间几乎不传递能量（非线性效应在一级近似中可忽略）。

为什么高斯假设有效？

由中心极限定理，大量独立随机波的叠加趋向正态分布。海面上叠加着来自四面八方、历经数百公里传播的无数涌浪——它们各自有独立的随机初相位。当叠加数目足够大，任意一点的高度自然收敛到高斯分布。
这与音频工程中的白噪声生成如出一辙：无数独立的微小随机冲击叠加，最终产生高斯分布的噪声信号。

基于这个认识，我们可以在频域中独立为每个波矢量 分配一个复振幅：

其中：

  • — 两个独立的标准正态随机数，提供随机性（每次生成不同的海面）
  • — 海浪频谱的平方根，控制各频率分量的振幅”音量”

这里 就是我们接下来要重点讨论的海浪功率谱（Power Spectrum）——它是整个 FFT 海洋模拟的物理灵魂。

空间频率 的几何意义

波矢量 指向波的传播方向，其模长 是波数， 为波长。大的 对应短波（高频细节），小的 对应长涌（低频大波）。在 的离散网格上， 的取值范围从 到 ，构成一张二维频率图，每个格点代表一种传播方向和波长的波分量。

2.2 Phillips 频谱：风与波的第一个约定

功率谱 回答了一个核心问题：在以 速度吹向 方向的风场下，波长为 、传播方向为 的波，应当占据多大的能量份额？

[Tessendorf 2001][1] 提出了一个形式优雅的答案：

要理解这个公式，不妨把它分解为三个相互独立的调控旋钮：

① ：低频主导旋钮

意味着能量随波数急速衰减。把波长 从 1m 增大到 10m（k 减小 10 倍），能量增大 倍。这与真实海洋中的观测一致：大涌远比小波携带更多能量。从数学上，这个 谱形状来源于 Kolmogorov 型湍流能量级联的海洋类比。

② ：短波截断旋钮

是由风速决定的特征长度。当 （即波长远小于 ），指数项趋近于 1，短波自由存在；当 （即波长远大于 ，超过最大可激发波长），指数项急剧趋近于 0，这些不存在的超大波被截断。

直觉理解：风速越高， 越大，越长的波能被激发。 8 m/s 的微风， m；30 m/s 的台风， m，海浪尺度相差一个数量级。

③ ：风向对齐旋钮

这是波传播方向单位向量 与风向单位向量 点积的平方。顺风传播的波（）获得最大能量系数 1；逆风传播的波（）点积为 ，平方后同样为 1（这是 Phillips 谱一个已知缺陷：它允许逆风波存在）；垂直风向传播的波能量为零。

从频谱到视觉效果：改变参数会发生什么？

把频谱想象成海面波浪的”音频均衡器”：

  • 风速翻倍（ m/s）： 变为 4 倍，截断波长变长，海面出现更大的涌浪，整体波高约增加 4 倍。
  • 振幅常数 翻倍：所有尺度的波等比例增强，海面整体变得”汹涌”，但频谱形状不变。
  • 风向旋转 90°：能量分布在 k 空间内旋转 90°，海浪涌来的方向改变，但波长分布不变。

这三个参数分别控制了海浪的尺度分布、总体强度和传播方向，在 Compute Shader 实现中对应三个独立的可调量。

Phillips 频谱的局限性

Phillips 谱假设海浪”充分发展”（fully developed），即风无限持续、无边界限制。这导致两个工程上的问题：

  1. 高频端衰减过快：真实海面在高频段的能量比 衰减更缓慢（约 量级），导致 Phillips 谱生成的海面高频细节不足，显得”太光滑”。
  2. 允许逆风波： 对顺风和逆风对称，但真实海洋中顺风波远强于逆风波。可通过额外添加 因子（ 控制方向性）来修正。

2.3 JONSWAP 频谱：来自北海的工程校准

Phillips 谱是一个”理论优先”的模型。1968–1969 年，五国联合海洋考察队在北海（Joint North Sea Wave Project）实地测量了数百组真实海浪数据，从中拟合出更精确的半经验模型——JONSWAP 频谱[2]：

基底峰值增强因子

关键参数的物理直觉：

  • ：谱峰角频率， 为 10m 高度处的风速，与风区长度（fetch）直接相关
  • ：峰值增强因子，这是 JONSWAP 最显著的特征——峰值处能量约是 PM 谱的 3.3 倍，波形更”集中”
  • ：Phillips 常数，由风速和风区长度联合决定
  • （），（）：谱峰不对称宽度

“fetch”（风区）是理解 JONSWAP 的关键词。 它指的是风在开阔水面上持续吹拂的距离：

直觉类比

把海浪的生长想象成一个”共振系统”。风开始吹起的瞬间，海面只有短小的毛刺波（高频），它们吸收风能并逐渐向低频传递能量——就像钢琴的高频键先被拨响，能量慢慢向低频键”流淌”。风区越长，这个能量级联过程越充分，谱峰频率越低（对应波长越长的主波）。

当能量级联最终达到平衡时，就形成了 Pierson-Moskowitz（PM）谱，即 JONSWAP 公式中的基底项。PM 谱对应”充分发展”的海浪；而 JONSWAP 的 因子正是风区有限时，谱峰额外聚集的那部分能量。

风区长度 谱峰周期 典型场景
10 km ~4 s 内湾、湖泊
100 km ~7 s 近海陆架
500 km ~10 s 大洋边缘
充分发展 ~12–15 s 开阔大洋（PM 谱）

数据基于 JONSWAP 拟合关系 ， 为风区长度。 越大，峰值越尖锐，对应海浪越”集中”于主波长，视觉上涌浪感更强。

方向扩散函数：上面的 是一维频谱（只考虑能量随频率的分布），真实海面的波从四面八方涌来。将能量分配到不同传播角度 ，需要乘以方向扩散函数 ：

其中 控制方向集中程度：低频长涌方向性极强（ 大，能量集中在风向附近），高频毛刺波方向分散（ 小，接近各向同性）——这与海上观测完全一致：远洋涌浪有明确的来向，而近处的风浪则乱成一片。

最终的二维功率谱为：

其中 Jacobian 因子 完成从 域到 域的变量替换（雅可比行列式）。

FFT 框架中水深的工程约束

传统 FFT 海洋框架将水深 视为全局常量——这是由 IFFT 的数学本质决定的：频域到空域的变换是对整张网格同时执行的，无法针对不同空间位置使用不同的水深值。如果引入空间变化的水深 ，弥散关系 也将随位置变化，标准全局 IFFT 从根本上失效。

工程结论： 若需模拟复杂近岸地形（岸礁浅水区的波浪折射、水深渐变下的波速减慢等），必须引入 SWE 浅水方程或局部变形网格，而不能在 FFT 框架内直接实现。

2.4 弥散关系：为什么大浪比小浪跑得快？

弥散关系（Dispersion Relation）是波动力学的基石，它建立了波的空间特征（波数 ）与时间特征（角频率 ）之间的约束：

在实时渲染中，水深 ，，忽略毛细效应（，在波长 > 1.7cm 时成立），退化为极为简洁的深水弥散关系：

这条公式隐藏着一个反直觉的结论。

波的相速度（波峰移动的速度）为：

越小 → 波长 越大 → 相速度越高。

在深水中，长波比短波传播得更快。

10m 波长的波： m/s（约 14 km/h）
100m 波长的涌浪： m/s（约 45 km/h）
1000m 波长的海啸（极端情况）： m/s（约 141 km/h）

一个经典的现实例证

台风过后，距离中心数千公里的海岸往往会比台风本身提前 1–2 天感受到巨浪。这正是弥散的结果：台风激发了各种波长的波，其中长涌（波长 200~500m，周期 12–18s）以 15–25 m/s 的速度超前传播，远比台风中心移速（通常 5–10 m/s）更快，成为台风的”前奏”。气象员正是通过追踪波浪周期来反推风暴位置。

弥散关系之所以对 FFT 模拟至关重要，在于它决定了每个频率分量如何随时间演化——这正是下一节的主题。

2.5 时间演化：每个波分量独立地”旋转”

有了初始频谱 ，如何让海面动起来？

物理上，海面是一系列平面波的叠加，每个平面波以固定角速度 传播。在频域中，时间演化等价于让每个复振幅乘以一个以角速度 旋转的相位因子（Phase Factor）：

经过时间

但仅此一项不够——我们需要确保变换回空间域后得到的高度值 是实数（不能有虚部，因为水面高度是真实物理量）。

为什么需要复共轭项？

傅里叶变换有一个基本性质：若信号 为实数，则其频域表示满足厄米对称性（Hermitian Symmetry）：

简单说： 方向的复振幅与 方向的复振幅必须互为共轭，才能使两项的虚部在叠加时相消。

时间演化后，这个对称性必须依然保持，因此我们令：

第二项是对 方向振幅取共轭、并以 速度反转旋转——恰好保证了整体的厄米对称性，从而确保 。

欧拉公式视角下的时间演化：

可以把每个频率分量 想象成复平面上的一个旋转矢量（相量）：

  • 高频分量（大 ，大 ）旋转快，对应涟漪的快速变化
  • 低频分量（小 ，小 ）旋转慢，对应大涌的缓慢起伏
  • 所有分量以各自的 独立旋转，互不干扰

这也解释了为什么每一帧的时间演化 Pass 如此廉价：只需对每个格点执行一次 sincos(omega * t)——这正是 TimeDependentSpectrum.compute 的全部工作。

2.6 从频域到空间域：IFFT 是如何”拼装”海面的

现在所有频率分量都有了在时刻 的复振幅 。最后一步：把它们全部叠加，得到空间域中每个格点 处的海面高度。

这正是逆傅里叶变换（IFFT）的定义：

直觉上：每个频率分量 是一个二维复数平面波， 描述它在空间中的分布（沿 方向以波长 振荡）。将所有分量在每个空间点处”叠加投票”，就得到了最终高度。

写成 离散形式：

这是一个 规模的矩阵乘法——朴素计算需要 次运算（ 时约 次），完全无法实时。而 FFT 将其降至 ，这正是第 3 节的主题。

为什么 FFT 海面可以无缝平铺（Tileable）？

离散傅里叶变换的输入/输出天然具有周期性边界。 网格的频谱定义了一个以域大小 为周期的空间信号。将这张网格在 xz 平面上无限平铺，相邻块在边界处自动连续——频域的周期性保证了空间域的无缝衔接。这是 FFT 海洋区别于 SWE 的天然优势之一。

2.7 Choppy 波：水粒子的”圆周运动”

如果只有高度偏移 ，得到的海面会呈现完美的正弦波起伏——波峰和波谷一样圆润，像被熨斗烫过的波纹布。但真实海浪远不是这样：波峰尖锐、波谷宽缓，在大风下甚至呈现出刀峰状的”白头浪”。

这种不对称性的物理根源，在于水粒子并非原地上下运动，而是做圆周运动。

水粒子的轨迹

在深水线性波理论中，水面附近的水粒子沿圆形轨迹运动（深水）或椭圆轨迹（浅水）：

  • 当粒子在轨迹顶部时，水平速度与波传播方向同向（被”推着走”），粒子在水平方向被向前压缩，多个粒子”挤”向同一处——形成波峰的尖锐聚拢
  • 当粒子在轨迹底部时，水平速度与波传播方向反向（被”拉着退”），粒子在水平方向被拉开——形成波谷的宽缓展开

这种水粒子的整体平移（叠加到垂直运动之上）就是水平位移（Choppy displacement）。

Gerstner 波的正确画法

分量 只有垂直位移 加上水平位移（Choppy）
波峰形状 圆弧 尖锐 ✓
波谷形状 圆弧 宽缓 ✓
视觉真实感 人工感强 接近真实 ✓
白沫生成 无从判断 Jacobian 可检测 ✓

Gerstner 波正是将这种水平位移参数化的结果。FFT 框架通过在频域中引入位移频谱，等效实现了同样的效果。

频域中的水平位移推导：

我们知道高度场的梯度（斜率）为 ，而水平位移的方向沿 （水向”低处”聚集）。在频域中，梯度运算等价于乘以 ：

水平位移与梯度方向一致但有 90° 相位差（位移超前于斜率，波峰处水平速度最大而斜率为零），在频域中体现为乘以 ：

类似地，高度场梯度（直接用于法线计算）：

频域乘法 = 空域微分的深层含义

乘法是傅里叶变换中微分定理的直接体现：对空域信号求偏导 ，等价于在频域中乘以 。这意味着高度图 、水平位移 、曲面法线 、以及后续 Jacobian 所需的二阶偏导数，全部只需在频域中进行简单的乘法操作，在同一次 IFFT 管线里一并生成。这是 FFT 海洋模拟高效性的核心来源之一。

Choppy 位移的强度由参数 （Chop 系数，）控制：最终空间位移为 。 退化为纯高度场（正弦波形）， 时水平位移达到物理极限——再大就会导致海面”自交叉”（波浪卷曲），这恰好是第 5 节 Jacobian 白沫生成的物理起点。

2.8 二维频谱可视化

为了直观感受频谱的形状，下面的交互演示渲染了 空间中的二维功率分布 ，颜色越亮能量越高。拖动滑块观察风速和风向如何塑造频谱形态。

📡 二维海浪功率谱可视化（k 空间）

颜色越亮表示该波矢方向/波长的能量越高。中心为低频（长波），边缘为高频（短波）。箭头指示风向。

能量对数标度（暗→亮）
低高
L = V²/g ≈ 14.7 m
Fig 1.

二维 空间中的海浪功率谱。图像中心为原点（零频率），亮色区域能量集中处对应主要波浪方向和波长。菊形（Donut）分布是 Phillips/JONSWAP 谱的典型形态：低频端被截断（无超过最大可激发波长的波），高频端快速衰减（短波能量低）。JONSWAP 模型谱峰更集中，对应涌浪感更强。

3. Cooley-Tukey 蝶形算法

3.1 计算复杂度：从 到

朴素的 DFT 对 点序列的计算量是 ——对每个输出频率 ，都需要对所有 个输入点求和。对于二维 的海面网格，总计算量为 。

以 为例：

  • 朴素 DFT： 次复数运算
  • FFT： 次

FFT 的核心洞察是 Danielson-Lanczos 引理：一个 点 DFT 可以分解为两个 点 DFT 的组合，以此递归，直至基础情形（2 点 DFT）。

3.2 DIT（按时间抽取）分解

Cooley-Tukey DIT-FFT 将输入序列按奇偶下标递归分割：

偶数项奇数项

其中旋转因子（Twiddle Factor）定义为：

注意 ，因此上半部与下半部可共用同一对中间结果，这就是蝶形运算效率的来源。

3.3 蝶形运算单元

单个蝶形运算（Butterfly Operation）的核心是：

直觉理解：一次蝶形运算消耗 1 次复数乘法 + 2 次复数加法，同时产生 2 个输出。每个 Stage 有 个蝶形，总共 个 Stage，因此总复数乘加次数为 。

3.4 交互演示：8 点 FFT 蝶形图

下面的交互演示展示了一个 的 DIT-FFT 完整流程。点击「下一步」逐阶段推进，观察数据如何在 3 个 Stage 中从位反转排列逐步变换为最终频域结果。

🦋 Cooley-Tukey DIT-FFT — 交互蝶形图（N = 8）

输入序列经位反转后，经过 log₂8 = 3 个 Stage 完成变换。高亮蝶形表示当前 Stage 的运算单元。

Stage 0 / 输入（位反转排列）
当前操作：输入序列进行位反转（Bit-Reversal）排列，为分治蝶形运算做准备。
位反转示例：索引 001₂ → 反转为 100₂，即 1 ↔ 4。

3.5 Stockham FFT：GPU 友好的无冲突变体

经典 Cooley-Tukey 算法存在一个 GPU 实现上的问题：每个 Stage 都需要原地（in-place）写回，即读写同一块内存，在 Compute Shader 中会引发 bank conflict。

Stockham FFT 通过使用**双缓冲（Ping-Pong Buffer）**解决这一问题：每个 Stage 从一块纹理读取，写入另一块纹理，交替使用，始终保证访问的连续性和无冲突性。这正是 simulation.md 中 FFT.compute 采用 Stockham IFFT 的核心原因。

实践中的额外优化： 可以将旋转因子 预计算为一张 Butterfly Texture（大小 ），每个纹素存储 (index_top, index_bot, twiddle_real, twiddle_imag)，在每个蝶形 Pass 中直接采样，避免运行时的 exp() 计算。

4. Compute Shader 架构与工程实现

本节结合 Ocean-URP 的实际实现，逐 Pass 拆解 GPU 计算管线。整体架构分为四个阶段：

Pass 1 — InitialSpectrum.compute ⏱ 条件执行：首次启动 / 参数（风速、风向）变更时 输出: H0 纹理 (RGBAHalf, N×N×C) — 静态初始傅里叶系数 预计算 每帧执行 ↓ Pass 2 — TimeDependentSpectrum.compute ⏱ 每帧执行 输出: 频域位移 + 偏导数纹理 (N×N×2C) — 包含高度/水平位移/Jacobian频域数据 每帧 每帧执行 ↓ Pass 3 — FFT.compute（Stockham IFFT） ⏱ 每帧执行 流程: Fft (x轴) → Permute 转置 → Fft (z轴) → PostProcess (归一化) 输出: 空间域位移图 + 偏导数图 (用于法线 & Jacobian) 每帧 simulateFoam = true 时执行 ↓ Pass 4 — FoamSimulation.compute ⏱ 条件执行：simulateFoam = true 基于 Jacobian 行列式检测波浪卷曲，输出泡沫 / 湍流纹理 (N×N×C) 可选 (+15~20%)

4.1 Pass 1 — 初始频谱生成

初始频谱只需在参数（风速、风向）变化时重算，是整个管线中唯一的预计算阶段。

核心内核 CalculateInitialSpectrum 对每个 并行执行：

📝 InitialSpectrum.compute
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
[numthreads(8, 8, 1)]
void CalculateInitialSpectrum(uint3 id : SV_DispatchThreadID)
{
    // 1. 计算波矢量 k
    float2 k = float2(
        (id.x - N * 0.5f) * deltaK,
        (id.y - N * 0.5f) * deltaK
    );
    float kLen = length(k);
    float kAngle = atan2(k.y, k.x);  // 相对于风向的角度

    // 2. 生成高斯随机数（使用 Box-Muller 变换）
    float4 rand = gaussianRandom(id.xy, _Seed); // ξ_r + iξ_i
    
    // 3. 评估 JONSWAP 频谱 + 方向扩散
    float omega = sqrt(GRAVITY * kLen);         // 弥散关系（深水近似）
    float spectrum = JONSWAPSpectrum(omega, _WindSpeed, _Fetch);
    float dirSpread = DonelanBannerSpread(kAngle - _WindDir, omega);
    
    // 4. 计算短波衰减
    float fade = exp(-_ShortWavesFade * _ShortWavesFade * kLen * kLen);
    
    // 5. 输出初始傅里叶系数
    float amp = sqrt(2.0f * spectrum * dirSpread * abs(deltaK * deltaK)) * fade;
    float chop = _EqualizerChop * _GlobalChop;
    
    _H0K[id] = float4(rand.xy * amp, chop * rand.xy * amp);
    _WavesData[id] = float4(k.x, amp * chop * SafeDivide(k.x, kLen),
                             k.y, omega);
}

4.2 Pass 2 — 时间演化

每帧执行，将静态初始频谱按时间 推进：

📝 TimeDependentSpectrum.compute — 时间推进内核
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
[numthreads(8, 8, 1)]
void CalculateAmplitudes(uint3 id : SV_DispatchThreadID)
{
    float4 h0 = _H0[id];               // (h0.r, h0.g) = H0(k), (h0.b, h0.a) = H0(-k)^*
    float4 wavesData = _WavesData[id]; // (kx, amp·chop·kx/k, kz, ω)
    float omega = wavesData.w;
    
    // e^(iωt) = cos(ωt) + i·sin(ωt)
    float cosOmega, sinOmega;
    sincos(omega * _Time, sinOmega, cosOmega);
    
    // H(k,t) = H0(k)·e^(iωt) + H0(-k)^* · e^(-iωt)
    float2 h0k    = float2( h0.r,  h0.g);
    float2 h0mink = float2( h0.b,  h0.a); // 已是共轭
    float2 ht = ComplexMul(h0k, float2(cosOmega, sinOmega))
              + ComplexMul(h0mink, float2(cosOmega, -sinOmega));
    // ComplexMul(a, b) = float2(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x)
    // 即复数乘法 (a.r + i·a.i)(b.r + i·b.i) = (ar·br - ai·bi) + i·(ar·bi + ai·br)
    
    // 水平位移频域 Dx(k) = -i·(kx/k)·h(k,t), Dz(k) = -i·(kz/k)·h(k,t)
    float2 kNorm = normalize(wavesData.xz);
    float2 htChop = float2(ht.y, -ht.x) * wavesData.y; // -i·ht·(kx/k·chop)
    
    // 法线/导数频域 ∂h/∂x = i·kx·h, ∂h/∂z = i·kz·h
    // 输出到双切片纹理
    _Output[uint3(id.xy, 2*cascade)]   = float4(ht, htChop);   // 高度 + x位移
    _Output[uint3(id.xy, 2*cascade+1)] = /* ... z位移 + 偏导数 */;
}

4.3 Pass 3 — Stockham IFFT

对 Pass 2 的输出执行二维逆变换：

1
2
3
4
1. Fft(x-direction, inverse=true)   // 沿 x 轴 Stockham IFFT
2. Permute()                        // 转置数组（使 z 方向变为行方向）
3. Fft(z-direction, inverse=true)   // 沿 z 轴 Stockham IFFT
4. PostProcess(scale=1/N²)          // 归一化 + 符号修正

Shared Memory 优化

在 Compute Shader 中，每个 Stage 的蝶形运算涉及跨越越来越大的内存间距的数据访问。对于 的 IFFT，最后几个 Stage 会跨越整张 的 RWTexture，L1/L2 缓存命中率极低。

优化方法之一是将列数据加载到 groupshared 共享内存中，让所有 个 Stage 都在 SRAM 内完成，再一次性写回全局内存，大幅减少全局内存往返次数。具体实现限于篇幅不展开，可参考 GodotOceanWaves 中的 GLSL Compute Shader 实现。

4.4 Pass 4 — 生成最终纹理

IFFT 结果经 PostProcess 后，得到空间域的：

输出纹理通道 含义
rgb.x 高度偏移 （ 轴）
rgb.y 方向水平位移
rgb.z 方向水平位移
rgb.w 用于法线/Jacobian 计算的梯度数据

法线向量直接由偏导数纹理构造（无需额外 Pass）：

5. 白沫仿真：Jacobian 行列式

本节是全文重点。 白沫（Whitecaps）是海浪真实感的核心要素之一，其产生机制与波浪的**非线性卷曲（Overturning）**密切相关，而 Jacobian 行列式正是量化这一现象的数学工具。

5.1 物理直觉：什么时候产生白沫？

将海面视为一个随时间演化的变形映射：每个水粒子从初始位置 被水平位移场 推移到新位置：

其中 是 Choppy 波的强度系数（通常 ）。

这个映射的雅可比行列式（Jacobian Determinant）描述了局部面积的变形程度：

展开为二维显式形式：

5.2 Jacobian 的物理意义

Jacobian 值 物理解释 视觉表现
局部面积拉伸，水面平坦 正常海面
无形变，等积映射 正常海面
局部压缩，波峰聚拢 波峰收紧
完全折叠，奇点 波峰恰好卷曲
映射翻转，”穿透”发生 物理上不可能，意味着已应卷曲产生白沫

在物理上意味着水粒子发生了”穿越”，即模拟精度不足以处理真实的波浪翻滚。这恰好是白沫产生的判定条件：当 （阈值，通常取 到 ），视为该处波浪正在或已经卷曲，生成白沫。

阈值与网格分辨率的敏感性

白沫阈值 并非一个”放之四海而皆准”的常数。由于偏导数 等是通过离散有限差分在网格上近似计算的，其精度直接取决于 FFT 网格的分辨率（Grid Size）和所在的级联层级（Cascade）：

  • 低分辨率级联（如 Cascade 3，覆盖 2560m / 64 格）：每格代表约 40m，高频导数信息丢失严重， 的数值偏差较大，同样的 可能几乎不产生白沫。
  • 高分辨率级联（如 Cascade 0，覆盖 40m / 256 格）：导数精度高， 的区域会准确对应波峰， 效果更符合物理预期。

工程建议： 将 和泡沫衰减率 作为可调材质参数暴露在 Inspector 中，按级联分别配置，并结合目标硬件和镜头距离进行视觉校准，而非在代码中写死。

5.3 频域中计算偏导数

梯度偏导数 、 等同样可以在频域中高效计算，方法与法线推导完全类似——在频域中乘以对应方向的 ：

因此，Jacobian 所需的四个偏导数纹理可以在 Pass 2 的同一批 IFFT 中一并计算，无需额外的 Pass。

5.4 持久泡沫模型

真实海浪的白沫并不瞬间消失，而是有一个衰减过程（持续数秒）。在 FoamSimulation.compute 中实现了经典的衰减-积累模型：

持续衰减新增白沫
📝 FoamSimulation.compute — Jacobian 白沫内核
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
[numthreads(8, 8, 1)]
void Simulate(uint3 id : SV_DispatchThreadID)
{
    // 从位移+导数纹理中读取四个偏导数分量
    // 假设 _FftInOut 已将导数 packed 在 slice 1 的 rgba 通道中
    float4 deriv = _FftInOut.Load(uint4(id.xy, cascadeBase + 1, 0));
    // deriv.x = ∂Dx/∂x (λ·kx²/k)
    // deriv.y = ∂Dz/∂z (λ·kz²/k)
    // deriv.z = ∂Dx/∂z (λ·kx·kz/k)
    // deriv.w = ∂Dz/∂x (= deriv.z, 因为位移场是保守场的梯度)

    // 计算 Jacobian 行列式
    float ddx = deriv.x;  // λ·∂Dx/∂x
    float ddz = deriv.y;  // λ·∂Dz/∂z
    float dxz = deriv.z;  // λ·∂Dx/∂z
    
    float J = (1.0f + ddx) * (1.0f + ddz) - dxz * dxz;
    
    // 当前帧新增的白沫量（仅在 J < threshold 时产生）
    float foamThreshold = 0.0;
    float currentFoam = max(0.0f, foamThreshold - J);  // J 越负，白沫越强

    // 持久泡沫的衰减与积累（在对数空间更稳定，类似 HDR 曝光合并）
    float4 prevTurbulence = _Turbulence.Load(uint4(id.xy, cascade, 0));
    float persistentFoam = prevTurbulence.y;
    
    // 衰减-积累更新
    persistentFoam = persistentFoam * (1.0f - _FoamDecayRate * _DeltaTime) 
                   + currentFoam * _DeltaTime;
    persistentFoam = saturate(persistentFoam);

    // 写回：通道1 = 瞬时泡沫，通道2 = 持久泡沫
    _Turbulence[uint3(id.xy, cascade)] = float4(currentFoam, persistentFoam, 0, 0);
}

5.5 白沫在着色器中的使用

泡沫纹理在海面 Fragment/Surface Shader 中与基础 PBR 材质进行混合：

1
2
3
4
5
6
7
8
9
10
11
12
13
// 采样当前级联的泡沫纹理
float2 foamData = SAMPLE_TEXTURE2D_ARRAY(_Turbulence, sampler_Turbulence,
                                          uv, cascadeIndex).xy;
float instantFoam    = foamData.x;   // 瞬时泡沫（仅在 J<0 时出现）
float persistentFoam = foamData.y;   // 持久泡沫（缓慢衰减）

// 泡沫总量 = 两者混合，persistent 为主
float foamAmount = lerp(instantFoam, persistentFoam, 0.7);

// 与海面 Albedo 混合（白沫增加粗糙度、降低镜面反射）
float3 foamColor = float3(0.95, 0.95, 0.95); // 接近纯白
surfaceData.albedo = lerp(surfaceData.albedo, foamColor, foamAmount);
surfaceData.smoothness = lerp(surfaceData.smoothness, 0.1, foamAmount);

6. 渲染着色概要

本节仅对非 Jacobian 部分做概要梳理，重点在于工程落地的接口。

6.1 网格与级联

从 Compute Shader 输出的位移图需要驱动海面网格顶点。对于大范围海洋，通常采用多级联（Multi-Cascade）叠加：

级联 覆盖范围 主要贡献
Cascade 0 10~40m 毛刺细波（Ripples）
Cascade 1 100~160m 中等风浪
Cascade 2 400~640m 长周期涌浪
Cascade 3 1000~2560m 宏观涌浪背景

各级联的覆盖范围由 OceanSimulationSettings.CalculateCascadeDomains 自动计算（Auto 模式），也可手动配置。渲染时按视距权重混合多级联的位移与法线。

6.2 PBR 水面材质要点

镜面反射与粗糙度

水面属于高度光滑的电介质（Dielectric），基础 IOR ≈ 1.33，F0 ≈ 0.02。

  • 各向同性 GGX 微表面 BRDF（对于平静水面，粗糙度极低，接近镜面）
  • 菲涅尔效应极强：掠射角（Grazing Angle）下接近全反射
  • 高频法线叠加 Detail Normal Map 弥补 IFFT 分辨率上限

次表面散射（SSS）近似

海水具有蓝绿色散射，波峰处透光感显著：

  • 使用基于深度的颜色衰减（Beer-Lambert 近似）
  • 在背向光源的浅水区增加散射项，模拟波峰透光
  • 对波峰应用额外的加法散射项（）

7. FFT vs SWE：技术路线对比与适用场景

两种技术并非竞争关系，而是互补的——大型项目中往往将两者结合：深海用 FFT，近岸/交互区域用 SWE 或距离场碰撞。

维度 FFT 海洋 SWE 浅水方程
物理模型 线性波叠加（频谱统计） 非线性偏微分方程
适用场景 开阔深水、无边界大洋 近岸浅水、池塘、河流
随机性 高（统计驱动） 低（需额外扰动）
物理交互 ❌（无法与动态物体交互） ✅（船只尾迹、角色涉水）
波浪折射/衍射 ❌ ✅
GPU 开销 中（主要在 IFFT） 高（迭代求解 PDE）
无限平铺 ✅（天然可平铺） ❌（固定域大小）
代表实现 Sea of Thieves, AC: Black Flag GTA V 港口水域, 《荒野之息》河流

在 3A 项目中如何结合？

典型方案是以 FFT 渲染全局海面背景，在玩家周围约 50~200 米的交互半径内叠加一个低分辨率 SWE（或 iWave）模拟域：

  1. FFT 提供宏观波形的高度偏移；
  2. SWE 叠加局部交互扰动（船尾波浪、角色溅水、抛投物落水）；
  3. 两者在视图空间以距离权重混合，过渡区域足够平滑。

此外，对于海岸线附近，还可以引入基于 SDF（Signed Distance Field）的水-地碰撞，配合波浪衰减曲线模拟近岸波浪冲刷效果。

8. 性能优化与总结

8.1 分辨率与级联数选取

场景 推荐配置 内存估算
移动端 64×64 / 2 级联 ~2MB
PC（中配） 128×128 / 3 级联 ~12MB
PC/主机（高配） 256×256 / 4 级联 ~50MB
影视级 512×512 / 4 级联 ~200MB

8.2 Thread Group 划分原则

[numthreads(8, 8, 1)] 对 ~ 的纹理是合理的起点（每 group 64 线程，对应一个 NVIDIA Warp × 2）。对于 ，可考虑 [numthreads(16, 16, 1)] 以减少 dispatch 次数，但需注意 LDS（groupshared）的 occupancy 压力。

8.3 GPU Profiling：各 Pass 耗时参考

数据说明

以下耗时数据来自在 NVIDIA RTX 3060 Ti / Unity 2022 LTS + URP 环境下，对 Ocean-URP 进行 GPU Frame Capture（RenderDoc + Unity Frame Debugger）的实测参考值，配置为 256×256 / 4 级联，开启泡沫模拟。不同显卡、驱动版本和分辨率下数值会有差异，供量级参考。

Pass Compute Kernel GPU 耗时（参考） 瓶颈类型
Pass 1 InitialSpectrum CalculateInitialSpectrum × 4 cascade ~0.08 ms ALU Bound（exp/sqrt 密集）
Pass 2 TimeSpectrum CalculateAmplitudes × 4 ~0.06 ms ALU Bound（sincos × N²）
Pass 3 IFFT（x 轴） Fft × log₂256 × 4 = 32 dispatches ~0.22 ms Bandwidth Bound（大跨度读写）
Pass 3 Permute + IFFT（z 轴） Permute + Fft × 32 ~0.19 ms Bandwidth Bound
Pass 4 FoamSimulation Simulate × 4 ~0.07 ms ALU Bound
总计（每帧） — ~0.62 ms 以 IFFT 带宽为主

Pass 1 为条件执行，仅在参数变更时运行，不计入每帧开销。Pass 3 中 Permute（矩阵转置）虽然 ALU 量极低，但因访问模式随机，容易成为缓存效率瓶颈，建议在分析时单独关注其 L2 Cache Hit Rate。

带宽瓶颈 vs ALU 瓶颈的分析方式：

在 NVIDIA Nsight / AMD Radeon GPU Profiler 中，判断瓶颈的两个关键指标：

  • Achieved Bandwidth / Peak Bandwidth：若比值 > 70%，为带宽瓶颈（Bandwidth Bound）。IFFT 的 Pass 3 通常在此列，原因是蝶形算法越到后期 Stage 访问跨度越大，严重破坏 L1 的 Warp 内数据局部性。
  • ALU Utilization / SM Utilization：若比值 > 80%，为计算瓶颈（ALU Bound）。初始频谱 Pass 1 的 exp()、atan2()、sqrt() 密集调用会导致此情况。

应对策略：Butterfly Texture 预计算（将 Pass 1 的旋转因子 烘焙为查找表）可将 Pass 3 的 ALU 压力降低约 15%；Shared Memory 优化（将整列数据 load 进 groupshared）可将 IFFT 带宽消耗减少约 30~40%，代价是单 Dispatch 寄存器压力上升，需要平衡 occupancy。

8.4 关键优化点回顾

  1. Stockham FFT：双缓冲消除 bank conflict，比 Cooley-Tukey 原地版本在 GPU 上通常快 20~40%。
  2. Butterfly Texture 预计算：将旋转因子烘焙为纹理，避免每个蝶形单元的 exp() 运算（GPU 上 exp() 约 20~30 个 ALU cycle）。
  3. 频域批处理：高度、位移、法线、Jacobian 所需的所有频域数据在同一 Pass 中计算，仅需一次 IFFT 管线完成全部输出。
  4. 泡沫按需开关：泡沫模拟约增加 15~20% GPU 开销，在移动端或性能受限场景可完全关闭。
  5. 异步回读（AsyncGPUReadback）：CPU 侧仅在需要浮力物体交互时读取 12 个最低频级联的位移数据，延迟 12 帧，不影响渲染帧率。

8.5 总结

FFT 海洋模拟经历了二十余年的演进，从 Tessendorf 2001 的奠基性论文，到现代 GPU Compute Shader 的高效实现，其核心思想始终如一：

在频域中用物理统计模型定义能量，通过 IFFT 在空间域实例化随机海面，用 Jacobian 感知波浪破碎。

与其说这是一个”海洋渲染”问题，不如说是统计物理学、信号处理与 GPU 并行计算的交叉产物。

📚 参考文献与延伸阅读
[1]

Tessendorf, Jerry. “Simulating Ocean Water.” SIGGRAPH 2001 Course Notes #9: Simulating Nature: Realistic and Interactive Techniques. ACM, 2001. — FFT 海洋模拟领域的奠基性论文。

[2]

Horvath, Christopher J. “Empirical Directional Wave Spectra for Computer Graphics.” ACM SIGGRAPH 2015 Talks. — JONSWAP 频谱与 Donelan-Banner 方向扩散函数在图形学中的应用综述。

[3]

Elfouhaily, T., Chapron, B., Katsaros, K., Vandemark, D. “A Unified Directional Spectrum for Long and Short Wind-Driven Waves.” Journal of Geophysical Research, 1997. — 统一方向频谱（Unified Spectrum）的权威参考。

[4]

Ang, Nigel, et al. “The Technical Art of Sea of Thieves.” ACM SIGGRAPH 2018 Talks. — 工业级 FFT 海洋在 3A 游戏中的实际落地案例。

[5]

Gasgiant. Ocean-URP. GitHub, 2021–2024. https://github.com/gasgiant/Ocean-URP — 本文工程实现的主要参考，完整的 Unity URP FFT 海洋开源实现。

[6]

2Retr0. GodotOceanWaves. GitHub, 2022–2024. https://github.com/2Retr0/GodotOceanWaves — 基于 Godot 4 Compute Shader 的 FFT 海洋实现，GLSL 版本可作为 HLSL 的参照。

[7]

Cooley, James W., and John W. Tukey. “An Algorithm for the Machine Calculation of Complex Fourier Series.” Mathematics of Computation 19.90 (1965): 297–301. — Cooley-Tukey FFT 算法的原始论文。

[8]

Stockham, T.G. “High-speed convolution and correlation.” AFIPS ‘66 Spring Joint Computer Conference, 1966. — Stockham FFT（无冲突变体）的原始来源。

[9]

知乎专栏：FFT 海洋模拟学习笔记（含 Unity 实现） https://zhuanlan.zhihu.com/p/347091298 — FFT 海洋实现参考。