分形算法调和的海浪模拟方法

印桂生，王海玲，张菁，陈怀友

(哈尔滨工程大学 计算机科学与技术学院，黑龙江 哈尔滨 1 50001)

随着虚拟现实和计算机图形学技术的发展，动态海面作为三维游戏系统、虚拟现实系统、军事仿真系统的重要组成部分得到了多方的重视．众所周知，海面是一个形态复杂，变化不规则的，这给高真实感的波浪建模带来了一定难度，因此，如何在满足实时性要求的基础上进行高真实感的海浪建模成为了研究重点．目前，水面建模方法主要分为基于物理［1］、基于几何构造［2］和面向视觉效果［3］3 类方法．基于物理的方法遵行水体的运动原理，通过求解流体力学方程——Navier Stokes方程(NSF)来得到水面的速度场、压力场景，从而生成波浪运动形态，其方法公式求解复杂，计算量巨大，实时性较差，多用于数值分析和工程应用．为了避开大量的物理方程计算，学者们开始从波浪的几何形状入手，利用数学方法构造波浪形态，这种方法的构造函数相对简单，计算速度快，得到了游戏、影视特效等领域的青睐，但仍存在着波浪模型单一，波动规则的弊端．面向视觉效果的方法是采用纹理贴图或者图像叠加的方法来合成水面效果，这类方法是迎合实时性要求很高的产物，它形成的波形局限性较大，近视点观察时，视差效果消失，真实感降低．

虚拟现实系统对水面绘制的实时性和真实性要求较高，为了解决这一问题，提出一种分形算法调和的海面模拟方法．该算法将海浪模型分解为大尺度基波、小尺度细节和视觉效果三部分，通过余弦波叠加对基波建模，有效地降低了计算量;算法的核心是引入MPRD算法，充分利用分形算法自相似性和fBm特性来拟合海面的分形细节．文中不是简单地应用分形算法，而是通过对海浪的分析，提取具有海面特征的分形随机变量，用于优化分形参数，控制细节的变化;同时在多分辨率网格模型基础上，结合视点距离、人眼特性及视点变化速度等因素生成不同递归深度的细节层次，并根据设定的标准将基波与分形细节进行调和，合成的波形不仅具有各向同性的波浪特征，还具有各向异性的细节特性，实现了视相关的海面建模．算法力求完成实时的高真实感海面建模，在建模的过程中，利用GPU硬件特性进行海面视觉效果模拟，有效地保证了算法的绘制速度．

1 大尺度基波建模

根据线性波浪理论［4-5］，采用如下大尺度基波模型:

式中:H(x，z，t)表示多个余弦波组成的多向随机波;Ai是波幅;ki是波数，Ki=2π/Li，Li是波长;θi为波浪传播的方向角，θi∈［-π，π］;φi随机初始相位，φi∈［0，2π］;ωi为第 i项角频率，ω2=gKtanh(Kh)，对于深水波，=gK，h是水深，g 重力常数．Ls是平均波长，决定了高频波和低频波的比例;ħ为波形测度系数，决定海面形态．

式中:φ为风向角度，λ是波浪扩散程度的调整系统．为了提高波浪的动力学特性，在求解各种参数后，将在约束中随机产生某些参数，随着时间参数t的变化，生成大尺度波浪形态．

2 波浪的分形细节建模

分形理论(FT)是由Mandelbrot首先引入自然科学领域［6］．FT的自相似性具有膨胀/收缩对称性，适合描述各向同性的现象;而其自仿射性的膨胀/收缩的对称性不同，适合描述各向异性的现象，这些性质与海面波动特征相同，因此引入分形算法来添加波形细节具有很好的理论支撑．

2．1 MPRD算法的引入

MPRD是一种快速细分方法，假如海面网格的大小为(2n+1)×(2n+1)，网格间距为l00，对任一网格元，假设 Hij、Hij+1、Hi+1j、Hi+1j+1是其 4 个顶点的高度值，网格的MPRD分形过程为:首先进行网格元中心点的插值，网格元的中心点高度值Hic为

式中:Hit、Hir、Hib、Hil是与 Hic相邻的网格元中心点高度值，Ric、Rit、Rir、Rib、Ril是自定义的随机值位移，此时，网格间距为

式(5)是式(2)的特殊形式，说明MPRD可用于海面形态模拟过程．MPRD完成细分后，网格数量变为原来的22i倍，i=1，2，…，n 为递归深度，随着递归深度增加，分形细节越来越粗糙．从上面的公式可以看出，参与细分的数据有已知的顶点高度值和随机位移量，并者体现了点与点的相关性，而后者体现了点与点的差异性．

2．2 特征约束的随机变量

随机变量决定分形波浪的形态特征和多样性，以往的随机变量形式难以保证分形的局部稳态，因此，设计一个好的随机变量至关重要．由 FT知，MPRD满足fBm的幂指数规律:

那么，细分的点与任一相邻网格点满足上式关系，式(4)可表述为如下方式:

将式(7)展开，得

式中:H，σ是分形特征参数．由式(8)可求出随机变量的形式为

构建MPRD的最大难点在于消除分形曲面的折痕，保证变量取值的合理性十分关键，由式(9)及分形变量特征，定义MPRD的随机变量如下:

式中:l是调节函数，控制变量的随机性，Gauss为服从N(0，1)的高斯函数．

余弦叠加的波浪具有各向同性和各向异性，若估计出频域和方向域的变化程度，则可修正分形参数，控制分形结果．文中引入水面波向、视向及水面熵特征参数来改进分形特征，对于分形特征取值窗口W为M×N的海面网格，特征提取的公式如下．

1)波向特征值Hh．

式中:Δh 为采样步长，P(x，y)为海面坐标(x，y)处的高度值，mh(Δh)为分形窗口的相对高度值，保证了分形参数与领域顶点信息的相关性．不同的采样步长，可以得到对应的mh(Δh)值，根据随机过程属性有:

将式(12)取对数:

通过最小二乘法的逼近原理，对不同Δh和mh(Δh)拟合(lg‖Δh‖，lgmh)，使用拟合直线的斜率求出参数Hh的估计．

2)视向特征值Hv．

如上述方法，可求出参数Hv的估计．

3)水面熵特征值Hs．

为了体现分形区域内所有顶点的相关性，确保分形稳态性，对网格窗口 作整体度量，定义水面熵Hw:

hij是相对高度分布，∂是调节值，Pl=min{Pi;i，j=1，2，…，m}，同理可求出 Hs，σs．

FT指出，参数H可以反映曲面变形程度，其值越小，曲面越粗糙．为了计算简单，改进的分形特征表示如下:

2．3 多分辨分形

人眼在观察场景时存在视觉敏感区，且由于视点的运动性，海面细节无需过高．因此，方法在文献［7］的基础上改进，为不同区域的海面添加不同分形细节．

图1 视域划分Fig．1 View field division

首先根据视截体判定视场，如图1所示，V是视点位置，V'是视点垂足，V0是其在场景中的投影，he是视点高度，Pn是最近可视点，dn是最近可视点距离，Pf是最远可视点，df为可见的最远点距离，α视向角，视场角为fO．计算出可视区域后，假设Vs将进一步划分层次级别，而对Vs以外的区域直接剪裁．在可视区域内，顶点与视点的距离不同，视觉敏感度也不同，距离越近，敏感性越强，海面细节也应相对增加．则视域的敏感分层可表达如下:

式中:Δθ是不同敏感度间的递增角度，η为敏感层次数，Δθ=f0/2η，V0=Pn+(Pf-Pn)cos α，称为敏感点，rO是最高级的敏感区域半径，ri是不同LOD层的半径．

由Nyquist定理，设定网格与波长关系是1/4min(L/Li，i=1，2，…，n)≤L00＜ 1/2min(L|Li，i=1，2，…，n)，此外，视点的运行速度也是影响分形的重要因素，视点运动很快，对海面区域认知性较少，不需要对海面增加分形细节．因此，分形递归的影响函数可定义为

式中:ν是视点单位时间内的变化速度，θh是视点左右的变化角度，θv是前后移动时引起的视角变化．ν值越小，说明视点变化比较快，视觉效果分辨率较低，不用进行分形变化．μ(p，ε)分形的测度函数，A为波浪振幅，ε=κ×Δτ尺度的测度，Δτ采样间隔，f(ri，di)是与视点相关的调整因子，νf、vn分别是视域的调整系数，di是视线与场景的距离．

海面细分时，通过设置阈值来控制分形变化．正弦波和与分形算法细节加权求和的海面形态方程可写为 ζ=H+ γHf，ζ是组成波在时间 t时点(x，z)处的波形;H代表海浪谱波形;Hf表示分形细节结果;γ =cos(α +idθ)(i=0，1，2，…，η -1)是波浪调和系数，取值为［0，1］间的任意实数．可以证明，函数ζ是周期为2π的周期函数，当γ由0连续变化到1时，海浪也由细节丰富的波浪逐渐过渡到基于谱的波形．

3 面向GPU的海面视觉效果

光照是提高海面的视觉效果的重要因素，基于物理方法的实时计算海面光照属性计算量大．为了进一步提高计算效率，算法利用GPU技术特性实现海面的纹理贴图和光照效果．文中海面的视觉效果主要由3部分组成，近似表示为

式中:Ws海面基色，Wrfl是水面的反射效果，Wrfr是水面的折射效果，ψ是菲涅尔系数，表示了光在表面发生反射和折射现象的概率，决定两者对海面作用的相对比重，ψ∈［0．1，1．0］，其中，ζ= ηa/ηw-sin2α，ζ'=(cos α(ζ+cos α)-1)2，ζ″=(cos α(ζ-cos α)+1)2，ηa和ηw分别是空气和水的折射率，α是入射光与海面法向的夹角．

3．1 菲涅尔系数

海面是连续变化的，若实时计算菲涅尔系数会比较费时，如果将ψ值的求解在预处理阶段完成，则节省了大量实时运算时间．通过顶点shader获得海面顶点信息，将其转换到切空间坐标中，记录不同视向时的所有入射角，然后将这些信息保存在纹理中，然后制作成菲涅尔系数的纹理图，图像中的信息作为像素shader的输入，像素shader程度过程中针对海面顶点进行视觉效果融合．

3．2 海面基色

基色是海面的自身颜色，没有任何外界环境的影响，文中采用Wang tile贴图方法［8］进行映射海面的基色纹理．根据顶点shader记录的波浪高度值及视点向量，对纹理进行寻址，然后在像素shader中计算对应视点的海面顶点颜色．

3．3 反射效果

海面的反射效果采用平面反射贴图实现，以海面为裁减平面，将海面以上的场景渲染为反射贴图．然后由法线图得到法向量，根据视线方向计算入射光，并以此为依据计算投影纹理坐标，采样颜色信息，动态调整纹理．

3．4 折射现象

文中在反射贴图的基础上，计算折射网格，把海平面作为折射掩码记录到帧缓存的Alpha通道中，用以界定折射区域，使用法线图，施加一个很小的扰动值对纹理坐标进行位移，最后使用扰动后的纹理坐标查找对应的反射贴图中的值来实现折射效果．

4 实验结果分析

本文方法的硬件实验平台是奔腾4、3．00 GHz CPU、1GB内存和nVIDIA GeForce 6800显卡;海面模拟程序运行于Window XP系统下，软件环境是Visual studio2005、DirectX图形库和CG语言．

图2为引入分形算法的海浪模拟实验，图2(a)与图2(b)是固定参数σ，τ后，分形特征参数H对海面形态的影响，可以看出，H控制着海面波浪的起伏程度，H越小，波浪起伏幅度越小．图2(c)和图2(d)是固定参数H，τ后的海面形态的变化情况，可以看出，σ控制着分形值的分布，σ越小，波浪形态分布越均匀，说明只要合理的控制H，σ，就可生成更真实的波浪模型．图2(e)～(h)是根据波浪特性对分形参数进行特征约束后的实验效果，实验表明，窗口取值和分形影响函数的阈值都直接影响着分形效果，阈值越大，分形递归深度越大，波浪形态越多样．

图3是在小波浪级别情况下加入光照视觉效果的海面场景，图3(a)远视点的海面光照效果，图3(b)是近视点海面光照效果，在远点时，视角较大，反射光较多，折射光比例较少，海面颜色呈现偏蓝效果;而近视点时，视角较小，此时折射光相比反射光比较要大，海面颜色出现青色效果．

图3 海面的光照效果Fig．3 Wave surface with optic

在海面绘制效率方面，图4给出了3种海面模拟方法的绘制速度对比，实验数据是对100 s内的速率进行采样而得，可以看出，海面绘制的初始阶段，由于加载数据，初始化等，绘制帧率较低，随后绘制速度逐渐平稳．其中，Tessendorf方法和文中方法达到了较好的绘制速度，但文中方法比Tessendorf方法的绘制效率更高，在普通运行环境中，平均帧率可达60 f/s以上．

图4 波浪模拟方法的帧率对比Fig．4 Frame ratio of algorithms

在海面模拟的性能方法，图5给出了不同规模海面模拟，GPU加速与否的性能对比．图5(a)海面网格数据为64×64，大尺度波的叠加数为4，分形影响函数的阈值为0．3，实线表示利用GPU加速的模拟效率，虚线表示仅采用CPU的效率，由图可知，GPU加速的绘制速率更高，速率高于后者将近1．15倍．图5(b)是不同网格数量下的绘制速率，当波面网格数量小于256×256时，文中方法仍然能达到较好实时性，说明本文方法可以满足中等区域的海浪模拟．当网格规模大于256×256，CPU的绘制速率下降较快，这说明模拟规模对绘制效率影响较大，因此，采用GPU技术可有效提高绘制速度．

图5 算法的性能对比Fig．5 Comparison of algorithm efficiency

表1列出了不同LOD层数在不同分形阈值控制下的模拟效率．增加波数或增加LOD层数都会增加系统计算量，且分形函数也影响绘制速度．由表中数据看出，在网格规模小于256×256，波数为4时，方法可以得到高真实感的海面效果，且具有很好的绘制速度．

海面网格数LOD×4 LOD×5 0．3 0．5 0．3 0．5 64 ×64 67．42 60．27 63．19 60．23 128 ×128 60．26 57．31 56．24 51．27 256 ×256 52．64 47．39 49．82 46．61 512 ×512 40．83 37．68 37．93 34．37

5结论

提出一种高效的海浪模拟方法，将波浪分为大尺度基波和小尺度细节，通过分形算法增加海面细节，同时利用GPU实现海面的视觉效果．该算法有以下优点:

1)将海浪谱建模方法与分形算法相结合，实时了真实的波浪变形，并可实时控制变形程度．

2)对分形算法进行优化，利用海浪谱的频率特性、方向特性和海面曲度控制分形参数，并结合视相关、人眼敏感特性及视点变化等因素控制细节级别．

3)面向GPU实现光照效果，充分利用GPU的加速特性提高算法效率．文中将菲涅尔系数的计算结果存储为纹理图，减少了系统计算量．实验表明，文中算法可以实现实时绘制的高真实感波浪场景．

［1］CHEN J X，LOBO N V．Toward interactive-rate simulation of fluids with moving obstacles using Navier-Stokes equations［J］．Graphics Models Image Processing，1995，57(2):107-116．

［2］TESSENDORF J．Simulating ocean water［J］．Radiosity of the Ocean Environment，2004:2-25．

［3］PHARR M，FERNANDO R．GPU Gems 2:programming techniques for high-performance graphics and general-purpose computation［M］．［s．l．］:Addison-Wesley Professional，2005:283-294．

［4］俞聿修．随机波浪及其工程应用［M］．3版．大连:大连理工大学出版社，2003:131-137．YU Yuxiu．Random wave and its applications to engineering［M］．3rd．ed．Dalian:Dalian University of Technology Press，2003:131-137．

［5］HENRY D．On Gerstner's water wave［J］．Journal of Nonlinear Mathematical Physics，2008，15(2):87-95．

［6］MARTYN T．Realistic rendering 3D IFS fractals in realtime with graphics accelerators［J］．Computers＆ Graphics，2010，34:167-175．

［7］JOHANSON C．Real-time water rendering［D］．Lund:Lund University，2004:14-25．

［8］韩建伟，王青，周昆，等．基于Wang Tiles的几何纹理合成［J］．软件学报，2009，20(12):3254-3264．HAN Jianwei，WANG Qing，ZHOU Kun，et al．Wang Tiles based geometric texture synthesis［J］．Journal of Software，2009，20(12):3254-3264．