CN103413341A - 一种基于Gerstner模型的波浪折射模拟方法 - Google Patents

Abstract

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

Description

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

技术领域

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

背景技术

今年来，在计算机图形学和虚拟现实领域中，水面的实时真实感绘制一直是具有挑战性的研究热点之一。水面的模拟绘制在军事、能源、影视、广告等各领域之中都有着广泛的用途，人们对于水面模拟绘制的研究也在不断地深入，其中浅水域水面的模拟绘制是非常重要并且非常困难的方向之一。受到水底地形以及构筑物的影响，浅水域的水面会出现比深水域水面更加复杂的波动现象，如：波浪浅化（shoaling）、波浪折射（waverefraction）、波浪反射（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{\ωt})\\ z=z_0+R_0{e}^{{\mathrm{Kz}}_0}\cos ({\mathrm{Kx}}_0-\mathrm{\ωt})\end{array}\right.$

为了简化公式，我们可以假设

其他一些参数可以应用物理学原理计算得到： $\mathrm{\λ}=\frac{2\mathrm{\π}}{K}=\frac{{\mathrm{gT}}^2}{2\mathrm{\π}},A=2R\·$

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

$\frac{\mathrm{\&lambda;}}{{\mathrm{\&lambda;}}_{\&infin;}}=\tanh \left(\mathrm{Kh}\right)$

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

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

${\left(\frac{A}{A_{\&infin;}}\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 \mathrm{\&beta;}}$

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

$S_z=A\frac{2\mathrm{\&pi;}\&times;4\mathrm{\&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}{\mathrm{\&lambda;}/2}}$

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

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

$\mathrm{ShallowWaterWave}\left\{\begin{array}{c}\mathrm{kBlend}=\sqrt{\frac{h}{\mathrm{\&lambda;}/2}}\\ \tan \mathrm{\&beta;}=\frac{h}{\mathrm{dis}}\\ k\text{'}=k\&times;\mathrm{kBlend}+\mathrm{dir}\&times;(1-\mathrm{kBlend})\\ K=\frac{2\mathrm{\&pi;k}\text{'}}{\mathrm{\&lambda;}}\\ S_x=\frac{2\mathrm{\&pi;h}}{\cos \mathrm{\&beta;}}\\ S_z=A\frac{2\mathrm{\&pi;}\&times;4\mathrm{\&lambda;}}{S_x(1+\mathrm{\&Delta;h}-\tan \mathrm{\&beta;\&Delta;h})}\\ \mathrm{\&alpha;}={\mathrm{Kx}}_0-\mathrm{\&omega;t}\\ {\mathrm{Sdpos}}_x=S_x\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}+S_z\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ {\mathrm{Sdpos}}_z=S_z\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-S_x\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ \mathrm{\&Delta;h}=\mathrm{\&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}}{\mathrm{\&lambda;}}\\ \mathrm{\&alpha;}={\mathrm{Kx}}_0-\mathrm{\&omega;t}\\ {\mathrm{Ddpos}}_x=A\sin \mathrm{\&alpha;}\\ {\mathrm{Ddpos}}_z=A\cos \mathrm{\&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>\mathrm{\&lambda;}/2\&DoubleRightArrow;0\\ \mathrm{\&lambda;}/20\&le;h\&le;\mathrm{\&lambda;}/2\&DoubleRightArrow;\frac{\mathrm{\&lambda;}/2-h}{\mathrm{\&lambda;}/2-\mathrm{\&lambda;}/20}\\ h<\mathrm{\&lambda;}/20\&DoubleRightArrow;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 \mathrm{\&beta;}=\frac{h}{\mathrm{dis}}$ 公式2

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

$S_z=A\frac{2\mathrm{\&pi;}\&times;4\mathrm{\&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}{\mathrm{\&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}{\mathrm{\&lambda;}/2}}\\ \tan \mathrm{\&beta;}=\frac{h}{\mathrm{dis}}\\ k\text{'}=k\&times;\mathrm{kBlend}+\mathrm{dir}\&times;(1-\mathrm{kblend})\\ K=\frac{2\mathrm{\&pi;k}\text{'}}{\mathrm{\&lambda;}}\\ S_x=\frac{2\mathrm{\&pi;h}}{\cos \mathrm{\&beta;}}\\ S_z=A\frac{2\mathrm{\&pi;}\&times;4\mathrm{\&lambda;}}{S_x(1+\mathrm{\&Delta;h}-\tan \mathrm{\&beta;\&Delta;h})}\\ \mathrm{\&alpha;}={\mathrm{Kx}}_0-\mathrm{\&omega;t}\\ {\mathrm{Sdpos}}_x=S_x\sin \mathrm{\&alpha;}\sin \mathrm{\&beta;}+S_z\cos \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ {\mathrm{Sdpos}}_z=S_z\cos \mathrm{\&alpha;}\sin \mathrm{\&beta;}-S_x\sin \mathrm{\&alpha;}\cos \mathrm{\&beta;}\\ \mathrm{\&Delta;h}=\mathrm{\&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}}{\mathrm{\&lambda;}}\\ \mathrm{\&alpha;}={\mathrm{Kx}}_0-\mathrm{\&omega;t}\\ {\mathrm{Ddpos}}_x=A\sin \mathrm{\&alpha;}\\ {\mathrm{Ddpos}}_z=A\cos \mathrm{\&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>\mathrm{\&lambda;}/2\&DoubleRightArrow;0\\ \mathrm{\&lambda;}/20\&le;h\&le;\mathrm{\&lambda;}/2\&DoubleRightArrow;\frac{\mathrm{\&lambda;}/2-h}{\mathrm{\&lambda;}/2-\mathrm{\&lambda;}/20}\\ h<\mathrm{\&lambda;}/20\&DoubleRightArrow;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为波动模拟后质点的垂直高度。

Claims (8)

1.一种基于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，波浪折射效果和深水域的自由波动，最终模拟出浅水域到深水域无缝过渡的波浪折射效果。

2.如权利要求1所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第一步的步骤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.如权利要求1所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第一步的步骤3)中，根据每个地形点的高程h、方向dir和距离dis，构建新的DEM地形数据，然后计算相邻四个地形点高程h和方向dir的平均值，并结合四个地形点中距离dis的最小值，生成四叉树LOD层次结构的上一层地形数据，最后形成层次完全的四叉树LOD层次结构。

4.如权利要求1所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第二步的步骤3)中，椭圆波动轨迹的短轴和长轴计算公式如下:

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

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

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

Δh=λ/2-h

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

5.如权利要求4所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第二步的步骤3）中，折射波浪的波动方向k′计算公式如下：

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

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

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

6.如权利要求5所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第二步的步骤3）中，波浪折射效果的计算公式如下：

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

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

7.如权利要求6所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第二步的步骤4）中，波浪在深水域的自由波动计算公式为：

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

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

8.如权利要求7所述的基于Gerstner模型的波浪折射模拟方法，其特征在于，在所述第二步的步骤5）中，模拟浅水域波浪折射效果的公式为：

$\mathrm{FinalWave}\left\{\begin{array}{c}\mathrm{WaveBlend}\begin{array}{c}\left\{\begin{array}{c}h>\mathrm{\&lambda;}/2\&DoubleRightArrow;0\\ \mathrm{\&lambda;}/20\&le;h\&le;\mathrm{\&lambda;}/2\&DoubleRightArrow;\frac{\mathrm{\&lambda;}/2-h}{\mathrm{\&lambda;}/2-\mathrm{\&lambda;}/20}\\ h<\mathrm{\&lambda;}/20\&DoubleRightArrow;1\end{array}\right.\end{array}\\ 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为波动模拟后质点的垂直高度。

Images