一种基于Gerstner模型的波浪折射模拟方法

摘要

本发明公开了一种基于Gerstner模型的波浪折射模拟方法，包括：地形预处理，针对DEM地形数据，提取岸线点并计算地形点到岸线的方向和距离，再利用数据构建地形四叉树LOD层次结构；波浪折射的模拟绘制，在原始Gerstner理论模型的基础上，应用地形深度和地形点到岸线的距离模拟计算波浪波动高度，应用地形点到岸线的距离和方向模拟计算波浪的波动方向，最终得到实时高效的浅水域波浪折射的模拟绘制。本发明计算简单效率高，保证实时绘制的效率，同时保证了绘制真实感，非常适合应用于实时绘制系统。

说明书

技术领域

本发明涉及浅水域波浪折射效果的模拟绘制，尤其涉及一种基于 Gerstner模型的波浪折射模拟方法。

背景技术

今年来，在计算机图形学和虚拟现实领域中，水面的实时真实感绘制 一直是具有挑战性的研究热点之一。水面的模拟绘制在军事、能源、影视、 广告等各领域之中都有着广泛的用途，人们对于水面模拟绘制的研究也在 不断地深入，其中浅水域水面的模拟绘制是非常重要并且非常困难的方向 之一。受到水底地形以及构筑物的影响，浅水域的水面会出现比深水域水 面更加复杂的波动现象，如：波浪浅化（shoaling）、波浪折射（wave  refraction）、波浪反射（wave reflection）和破浪（wave break）等，因此， 浅水域波浪的实时真实感模拟绘制是非常困难的。

浅水域波浪折射的模拟绘制方法大致可以分为以下两类：

（一）基于参数方程的方法

Fournier和Peachey提出了第一个浅滩波浪折射模拟算法，并通过计算 一系列不同地形点对波浪波动相位的影响来模拟波浪折射现象，这种方法 不适用于有岛和海湾的地形。Ts’o利用波向跟踪（wave tracing）算法计算 波浪折射效果，并应用Snell折射定律处理波浪方向的变化来模拟浪折射， 波跟踪类似于光线跟踪，波在不同介质的传播过程中会因为折射改变波的 传播方向，波向跟踪算法假设在深水区域有一系列平行波作为入射波，然 后跟踪这些波进入浅水区域后的传播方向，并将每个点的传播方向和波动 变化信息保存在表格中。赵欣等基于小振幅波理论，提出了一种基于相位 的波向计算方法模拟浅滩浪折射，并且可以通过参数控制来模拟卷浪等效 果。Gonzato等改进了波追踪算法，在射线分开太大的地方新的射线以增 加精度。

该类方法模拟流程基本相同：首先，预处理地形，一般多采用波向跟 踪算法；然后，结合地形数据应用Gerstner理论来模拟浅水域波浪的折射 效果。该类方法在地形预处理阶段计算出了地形对波浪波动的影响，并且 在模拟波浪阶段采用高效的Gerstner理论模拟波浪波动。因此，该类方法 计算简单，高效，广泛应用于实时绘制系统。

（1）波向跟踪算法

波向跟踪算法类似于光线跟踪算法，它假设在深水区域有一系列平行 波作为入射波，这些平行波进入浅水域时速度和方向会发生变化，算法应 用Snell折射定律跟踪计算入射波的传播方向，并将入射波在每个地形点位 置的波动信息进行预处理，保存在表格中。该算法可以得到很好的波浪折 射效果，但是算法在实际应用中有两个不足：首先，算法计算量大，难以 应用于大规模数据处理；其次算法是以特定的波浪作为入射波来预处理地 形的，所以当入射波属性发生变化时，必须重新对地形进行预处理。

（2）Gerstner理论模型

Gernster理论的参数波动模型如下面公式所示，其中λ表示波长，A表 示浪高，表示波浪卷曲度，T表示波浪周期，h表示水深，K表示波数， ω表示角速度，β表示水底地形坡度。

$\left(\begin{array}{c}x=x_0-R_0{e}^{{\mathrm{Kz}}_0}\sin ({\mathrm{Kx}}_0-\mathrm{\omega t})\\ z=z_0+R_0{e}^{{\mathrm{Kz}}_0}\cos ({\mathrm{Kx}}_0-\mathrm{\omega t})\end{array}\right)$

为了简化公式，我们可以假设其他一些参数可以应用物理 学原理计算得到：$\lambda =\frac{2\pi }{K}=\frac{{\mathrm{gT}}^2}{2\pi },A=2R·$

Gerstner理论模型只适用于在水深无限的情况，但是在模拟浅滩水面 波动时水深是必须考虑的因素，根据研究表明，水波进入浅水区域（h<λ/2） 后，周期不变，波长减小。以下所有在无限水深情况下的变量都用∞下标 表示。考虑水深的波长变化公式如下：

$\frac{\lambda }{{\lambda }_{\infty }}=\tanh \left(\mathrm{Kh}\right)$

其中h表示x点的当前水深。

水波进入浅水区域，波浪的浪高也会发生变化。考虑水深的浪高计算 公式如下：

${\left(\frac{A}{A_{\infty }}\right)}^2=\frac{1}{\tanh \left(\mathrm{Kh}\right)(1+\frac{2\mathrm{Kh}}{\sinh \left(2\mathrm{Kh}\right)})}$

（二）基于物理的方法

基于物理的构建方法是通过求解流体力学方程来模拟水的波动效果， 并且通过修改初值以及边界条件来控制水面的流动，该方法能准确的模拟 浅滩水面波动，但是计算量大，很难满足实时绘制系统的需求。Thon应用 简化的2D Navier-Stokes方程模拟水面的水平运动，而水面的高层则用随机 函数扰动。Kass和Miller应用欧拉法求解二维浅水方程模拟水面波动。李 永进根据各种海浪区域的波动和绘制特点，分别应用不同的海浪波动模型 模拟四种不同的海面运动。

相对于第一类方法，该类方法能非常逼真的模拟浅滩水面波动，但是 计算复杂，计算量大，不适合实时绘制系统。

发明内容

本发明着重解决了地形预处理和波浪折射模拟两个问题：在地形预处 理方面，本发明设计了一种更适于波浪实时绘制的地形预处理方法，相对 于波向跟踪算法它有更快的预处理速度，并且当波浪发生改变时不需要重 新预处理；在波浪折射模拟方面，本发明在原有Gerstner理论模型的基础 上，综合考虑了水深和岸线等信息，高效地实现了波浪折射的实时模拟绘 制。

一种基于Gerstner模型的波浪折射模拟方法，分为以下步骤：

第一步：针对DEM地形数据进行地形预处理

1）提取DEM地形数据中的岸线点；

2）针对DEM地形数据中的每个地形点作如下计算：

a）搜索距离所述地形点最近的岸线点；

b）计算所述地形点到最近的岸线点的距离，记为dis；

c）计算所述地形点到最近的岸线点的方向，记为dir；

3）利用地形点的高程h、方向dir和距离dis，构建地形四叉树 LOD层次结构；

第二步：浪折射效果模拟绘制

1）在GPU中预设定水面几何的顶点，根据所述的顶点在四叉树 LOD层次结构中确定分辨率相适应的匹配层次；

2）在所述匹配层次中确定与每个顶点对应的地形高程h、方向 dir和距离dis；

3）根据第二步步骤2）中的高程h和距离dis，以及原始波浪高 度计算出椭圆波动轨迹的短轴和长轴；根据第二步步骤2） 中的高程h、方向dir和波长λ，以及水面几何中波浪的自由 波动方向计算折射波浪的波动方向；最后应用Gerstner理论 模型计算波浪在浅水域的波动，得到浅水域波浪折射效果；

4）应用Gerstner理论模型计算波浪在深水域的自由波动；

5）根据第二步步骤2）中的高程h，波浪折射效果和深水域的 自由波动，最终模拟出浅水域到深水域无缝过渡的波浪折射 效果。

在所述第一步的步骤1）中，在DEM地形数据中选取某一地形点(x, y），所述岸线点的提取方法如下：

a）若地形点(x,y）的高程h为0，则表示该地形点的位置(x,y）为所 述的岸线点；

b）若地形点(x,y）的高程h不为0，则读取地形点(x,y）相邻两点(x-1, y）和(x,y-1)的高程h₁₀和h₀₁，然后判断h与h₁₀的关系，如果h和h₁₀同 号，表示两者之间没有岸线点，异号则根据高程差值求出两点之间的岸线 点；再判断h与h₀₁的关系，如果h和h₀₁同号，表示两者之间没有岸线点， 异号则根据高程差值求出两点之间的岸线点。

在所述第一步的步骤3）中，根据每个地形点的高程h、方向dir和距 离dis，构建新的DEM地形数据，然后计算相邻四个地形点高程h和方向 dir的平均值，并结合四个地形点中距离dis的最小值，生成四叉树LOD 层次结构的上一层地形数据，最后形成层次完全的四叉树LOD层次结构。

在所述第二步的步骤3）中，椭圆轨迹的短轴和长轴计算公式如下：

$S_x=\frac{2\mathrm{\pi h}}{\cos \beta }$

$\tan \beta =\frac{h}{\mathrm{dis}}$

$S_z=A\frac{2\pi \times 4\lambda }{S_x(1+\mathrm{\Delta h}-\tan \mathrm{\beta \Delta h})}$

Δh=λ/2-h

式中：h为当前地形的高程，dis为地形点到最近岸线点的距离，A 为波浪高度，S_x为椭圆短轴，S_z为椭圆长轴。

在所述第二步的步骤3）中，折射波浪的波动方向k′计算公式如下：

k′=k×kBlend+dir×(1-kBlend)

$\mathrm{kBlend}=\sqrt{\frac{h}{\lambda /2}}$

式中：k为波浪的自由波动方向，dir为地形点到最近岸线点的方向， h为当前水深，λ为波长，k′为折射波浪的波动方向。

在所述第二步的步骤3）中，波浪折射效果的计算公式如下：

$\mathrm{ShallowWaterWave}\left(\begin{array}{c}\mathrm{kBlend}=\sqrt{\frac{h}{\lambda /2}}\\ \tan \beta =\frac{h}{\mathrm{dis}}\\ k\text{'}=k\times \mathrm{kBlend}+\mathrm{dir}\times (1-\mathrm{kBlend})\\ K=\frac{2\mathrm{\pi k}\text{'}}{\lambda }\\ S_x=\frac{2\mathrm{\pi h}}{\cos \beta }\\ S_z=A\frac{2\pi \times 4\lambda }{S_x(1+\mathrm{\Delta h}-\tan \mathrm{\beta \Delta h})}\\ \alpha ={\mathrm{Kx}}_0-\mathrm{\omega t}\\ {\mathrm{Sdpos}}_x=S_x\sin \alpha \sin \beta +S_z\cos \alpha \cos \beta \\ {\mathrm{Sdpos}}_z=S_z\cos \alpha \sin \beta -S_x\sin \alpha \cos \beta \\ \mathrm{\Delta h}=\lambda /2-h\end{array}\right)$

式中：A表示波浪的高度，λ为波浪的波长，β为当前地形的斜率，

K为波数，ω表示角速度，k′为折射波浪的波动方向，x₀表示质点的初始 水平位置，t表示时间，Sdpos_x表示应用波浪折射模拟x₀质点在t时刻的水 平位移，Sdpos_z应用波浪折射模拟x₀质点在t时刻的垂直位移。

在所述第二步的步骤4）中，波浪在深水域的自由波动计算公式为：

$\mathrm{DeepWaterWave}\left(\begin{array}{c}K=\frac{2\mathrm{\pi k}}{\lambda }\\ \alpha ={\mathrm{Kx}}_0-\mathrm{\omega t}\\ {\mathrm{Ddpos}}_x=A\sin \alpha \\ {\mathrm{Ddpos}}_z=A\cos \alpha \end{array}\right)$

式中：K为波数，k表示波浪的自由波动方向，λ为波浪的波长，A表 示波浪的高度，ω表示角速度，x₀表示质点的初始水平位置，t表示时间， Ddpos_x为应用深水域自由波动方程模拟x₀质点在t时刻的水平位移，Ddpos_z为应用深水域自由波动方程模拟x₀质点在t时刻的垂直位移。

在所述第二步的步骤5）中，模拟浅水域波浪折射效果的公式为：

$\mathrm{FinalWave}\left(\begin{array}{c}\mathrm{Waveblend}\left(\begin{array}{c}h>\lambda /2\Rightarrow 0\\ \lambda /20\le h\le \lambda /2\Rightarrow \frac{\lambda /2-h}{\lambda /2-\lambda /20}\\ h<\lambda /20\Rightarrow 1\end{array}\right)\\ x=x_0-({\mathrm{Sdpos}}_x(1-\mathrm{Waveblend})+{\mathrm{Ddpos}}_x\times \mathrm{WaveBlend})\\ z=z_0+({\mathrm{Sdpos}}_z(1-\mathrm{WaveBlend})+{\mathrm{Ddpos}}_z\times \mathrm{Waveblend})\end{array}\right)$

式中：h为当前水深，λ为波长，x₀表示质点的初始水平位置，z₀表示 质点的初始垂直位置（即水平面），x为波动模拟后质点的水平位置，z为 波动模拟后质点的垂直高度。

本发明设计的技术方案为：

（一）在地形预处理方面，本发明针对DEM地形数据，预计算了地形 点到岸线的方向和距离，这些信息的计算量小，非常适合大规模地形数据 的预处理；并且这些信息的计算只与地形本身相关，与波的属性无关，所 以当波的波动属性发生变化时，无需重新预处理的地形，这非常有利于不 同水面波浪的波动模拟。

本发明首先分析了波浪折射现象的形成原理，并得出，当（h<λ/2） 时，水底地形会使波长和波速变小，并且水深越小的地方波浪的波长和波 速变得越小，因此从水深为λ/2的等深线开始到水深为0等深线（即为岸线） 为止波浪折射现象会变得越来越明显，最终表现为波浪线与岸线重合，这 时波浪的传播方向与岸线垂直。基于以上分析，本发明在地形预处理阶段 首先提取出了岸线，然后结合岸线信息预计算了地形点深度h、地形点到 最近岸线点的方向dir和距离dis三个参数，其中地形点深度h用于判断水底 地形是否影响波浪传播并且计算波浪在地形影响后的波长和波速；地形点 到最近岸线点的方向dir决定了波浪折射现象模拟时，波浪的最终传播方 向；地形点到最近岸线点的距离dis决定了波浪原始传播方向到dir的趋向 程度，当距离dis为0时波浪传播方向将完全趋向于dir；距离dis还影响波浪 波高的变化。从以上阐述可以得出，相对于波向跟踪算法，本发明的地形 预处理方法计算简单,保证了计算效率，并且地形预处理完全独立于具体波 浪，当波浪属性发生变化时，无需重新预处理，因此该方法非常适用于大 规模地形数据的实时绘制系统。

（二）在波浪折射的模拟绘制方面，本发明在原始Gerstner理论模型 的基础上，应用地形深度和地形点到岸线的距离模拟计算波浪椭圆波动轨 迹参数，应用地形点到岸线的距离和方向模拟计算波浪的波动方向，最终 得到实时高效的浅水域波浪折射的模拟绘制。

附图说明

图1为浅水域波浪折射模拟绘制流程图；

图2a和图2b为地形预处理算法示意图；

图3为波浪折射效果模拟结果图。

具体实施方式

本发明模拟波浪折射的算法流程如图1所示，其实现方法具体分为以 下几步：

（1）岸线的提取

地形的预处理首先需要知道地形的岸线信息，因此，本发明首先针对 DEM地形数据提取岸线点，岸线点的提取流程如下：

如果当前地形点(x,y）的高程h为0，则将该地形点的位置（x,y）加 入到岸线点集合中；如果h不为0，则读取相邻的左上两顶点(x-1,y）和(x, y-1)的高程h₁₀和h₀₁，然后判断h与h₁₀的关系，如果h和h₁₀同号，两者 之间没有岸线点，异号则根据高程差值出高程点(x’,y)加入到岸线点集合 中，对h₀₁也作相同的处理。岸线点提取示意图如图2a所示。

（2）地形数据计算

在已知岸线的基础上，对DEM地形中的每个地形点查找距离地形点 最近的岸线点，并计算该岸线点到地形点的方向dir以及两者之间的距离 dis，然后将地形高程h、方向dir和距离dis保存到文件中，生成扩展的 DEM地形数据，后续的计算均基于扩展的DEM地形数据，地形数据处理 示意图如图2b所示。然后进一步通过计算相邻4个地形点高程h和方向 dir的平均值，以及计算4个地形点中dis的最小值，形成四叉树LOD结 构中的上一层次地形数据。最后形成层次完全的预处理地形的四叉树LOD 层次结构。

（3）浅水域的判断与过渡

根据研究表明，波浪在深水区域是可以自由波动的，但对于的波浪， 当水深h<λ/2时，水底地形开始影响波浪波动，当水深h<λ/20时，波浪 为浅水波动，因此，本发明在这一理论的基础上，判断波浪何时出现浅水 波动，并且认为当λ/20≤h<λ/2时，波浪的自由波动和浅水波动按一定系数 WaveBlend融合，WaveBlend为0是完全为自由波动，越接近1则越趋向 于浅水波动，融合系数WaveBlend的计算公式如下：

公式1

在公式1中，h为水深，λ为波浪的波长。

（4）椭圆轨迹计算

当波浪进入浅水域时，由于水底地形的摩擦作用波浪的波动能量会减 小，波浪的波高、波长以及波速等属性都也会随之发生变化，并且波浪中 质点的运动轨迹也会由圆形变为椭圆形。本发明在预处理地形的基础上模 拟了这一系列浅水域波浪的变化，具体公式如下：

$\tan \beta =\frac{h}{\mathrm{dis}}$公式2

$S_x=\frac{2\mathrm{\pi h}}{\cos \beta }$公式3

$S_z=A\frac{2\pi \times 4\lambda }{S_x(1+\mathrm{\Delta h}-\tan \mathrm{\beta \Delta h})}$公式4

其中，Δh=λ/2-h

在公式2中，β为当前地形的斜率，h为当前水深，dis为地形点到最 近岸线点的距离；

在公式3中，β为当前地形的斜率，h为当前水深，S_x为椭圆轨迹的 短轴；

在公式4中，A表示波浪的高度，λ为波浪的波长，β为当前地形的 斜率，S_z为椭圆轨迹的长轴，且(1+Δh-tanβΔh)这一系数模拟了传播过程 中能量的减少。

（5）折射波浪波动方向k′计算

由波浪折射现象的原理可知，波浪会因为水底地形的影响改变传播方 向，从深水域的自由波动方向逐渐过渡到与岸线垂直的方向，本发明在地 形预处理数据的基础上，结合波浪的自由波动方向，模拟了波浪波动方向 逐渐变化的过程。在现有理论的基础上，通过观察与实验得出：波浪从水 深为λ/2时传播方向开始改变，水深为0时传播方向与岸线垂直，具体计 算公式如下：

$\mathrm{kBlend}=\sqrt{\frac{h}{\lambda /2}}$公式5

k′=k×kBlend+dir×(1-kBlend)公式6

在公式5中：kBlend为波动方向混合系数，h为当前水深，λ为波浪的 波长；

在公式6中：k′为折射波浪的波动方向，k为自由波动方向（设定值）， dir为地形点到最近岸线点的方向。

（6）波浪折射模拟

在步骤（4）和步骤（5）的基础上，应用Gerstner理论，可以计算出 波浪折射效果，计算公式见公式7。然后根据步骤（3）中的计算混合波浪 折射效果和波浪自由波动，波浪的自由波动可由Gerstner模型直接计算得 到，计算公式见公式8，最后进行混合计算，计算公式见公式9。最后的 模拟绘制结果如图3所示。

$\mathrm{ShallowWaterWave}\left(\begin{array}{c}\mathrm{kBlend}=\sqrt{\frac{h}{\lambda /2}}\\ \tan \beta =\frac{h}{\mathrm{dis}}\\ k\text{'}=k\times \mathrm{kBlend}+\mathrm{dir}\times (1-\mathrm{kblend})\\ K=\frac{2\mathrm{\pi k}\text{'}}{\lambda }\\ S_x=\frac{2\mathrm{\pi h}}{\cos \beta }\\ S_z=A\frac{2\pi \times 4\lambda }{S_x(1+\mathrm{\Delta h}-\tan \mathrm{\beta \Delta h})}\\ \alpha ={\mathrm{Kx}}_0-\mathrm{\omega t}\\ {\mathrm{Sdpos}}_x=S_x\sin \alpha \sin \beta +S_z\cos \alpha \cos \beta \\ {\mathrm{Sdpos}}_z=S_z\cos \alpha \sin \beta -S_x\sin \alpha \cos \beta \\ \mathrm{\Delta h}=\lambda /2-h\end{array}\right)$公式7

在公式7中：A表示波浪的高度，λ为波浪的波长，β为当前地形的 斜率，K为波数，ω表示角速度，k′为折射波浪的波动方向，x₀表示质点 的初始水平位置，t表示时间，Sdpos_x表示应用波浪折射模拟x₀质点在t时 刻的水平位移，Sdpos_z应用波浪折射模拟x₀质点在t时刻的垂直位移。

$\mathrm{DeepWaterWave}\left(\begin{array}{c}K=\frac{2\mathrm{\pi k}}{\lambda }\\ \alpha ={\mathrm{Kx}}_0-\mathrm{\omega t}\\ {\mathrm{Ddpos}}_x=A\sin \alpha \\ {\mathrm{Ddpos}}_z=A\cos \alpha \end{array}\right)$公式8

在公式8中，K为波数，k表示波浪的自由波动方向，λ为波浪的波 长，A表示波浪的高度，ω表示角速度，x₀表示质点的初始水平位置，t 表示时间，Ddpos_x应用深水域自由波动方程模拟x₀质点在t时刻的水平位 移，Ddpos_z应用深水域自由波动方程模拟x₀质点在t时刻的垂直位移。

$\mathrm{FinalWave}\left(\begin{array}{c}\mathrm{Waveblend}\left(,\begin{array}{c}h>\lambda /2\Rightarrow 0\\ \lambda /20\le h\le \lambda /2\Rightarrow \frac{\lambda /2-h}{\lambda /2-\lambda /20}\\ h<\lambda /20\Rightarrow 1\end{array}\right)\\ x=x_0-({\mathrm{Sdpos}}_x(1-\mathrm{WaveBlend})+{\mathrm{Ddpos}}_x\times \mathrm{WaveBlend})\\ z=z_0+({\mathrm{Sdpos}}_z(1-\mathrm{WaveBlend})+{\mathrm{Ddpos}}_z\times \mathrm{WaveBlend})\end{array}\right)$公式9

在公式9中：h为当前水深，λ为波长，x₀表示质点的初始水平位置， z₀表示质点的初始垂直位置（即水平面），x为波动模拟后质点的水平位置， z为波动模拟后质点的垂直高度。