海浪预报三维动画计算原理与制作方法

高志一，于福江，许富祥

（国家海洋环境预报中心 北京 100081）

海浪预报三维动画计算原理与制作方法

高志一，于福江，许富祥

（国家海洋环境预报中心 北京 100081）

研究了基于线性海浪模型制作三维动画海浪预报产品的制作方法，并制作了首个三维动画海浪警报产品。基本步骤为：由线性海浪模型生成一系列不同波高和波长的二维海面高度场，以及这些波面高度场对应的反光系数场；根据海浪预报图中波高的空间分布情况将不同波高的波面高度场拼接成复合的波面高度场，同时拼接出与之对应的复合反光系数场；由复合的波面高度场和反光系数场渲染出三维波面，并将渲染后的波面时间序列串接成动画形成最终产品。三维动画形式的海浪预报产品能有效地克服海浪有效波高等值线图信息不直观的缺点，预报信息易于被广大用户和公众理解。

海浪预报；FFT；海浪渲染

Abstract: A method of preparing an animation of 3-D wave forecasting is explored.A series of wave surface elevation fields with different mean wave height and wave length are generated on the basis of linear wave model, and also the corresponding reflective fields are generated.A composite wave surface elevation field and the corresponding reflective coefficient field are generated according to the significant wave height contour map.The frames of 3-D wave field are rendered based on these composite wave surface elevation fields.A series of rendered frames are animated at last.The 3-D animation of wave forecasting is much easier to understand for the public.

Keywords: Sea Wave Forecasting; FFT; Sea Wave Render

海浪预报服务最初是为满足海上军事活动的需要而出现的，迄今已有数十年的历史。由于海浪的随机特性，对它只能进行统计预报。海浪预报产品需要告诉用户洋面或海区平均波高的情况。因此海浪预报产品一般采用有效波高等值线图，该形式的产品实际上是一种传统涉海行业中广泛采用的专业用图。然而随着经济进步和人类海上活动日益频繁，需要海浪预报服务的用户群体早已突破传统涉海行业。由于新兴用户可能具有不同的专业背景，对他们来说波高等值线图不直观甚至可能产生误导。如图1a为一个等高线示意图，该图的含义是颜色越浅的等值线圈出区域的有效波高越大，但一般用户不容易想象出图1a对应的浪场波动状态，而是容易将它想象成图1b所示的三维图景，即“一座水山”。由此可见等值线形式的海浪预报图已经不能满足面向一般用户的预报服务需要，新的预报产品应该具有更直观生动的形式。这就要求新的预报产品不仅具有传统预报产品的功能，更要便于一般用户了解海浪预报中的深层次信息，如等高线大值区的随机波动具有较大的平均波高和波长。因此三维动画形式是新型海浪预报产品首选形式，然而目前尚无此类产品面世。

图1 a: 图等高线图。b: 图为等高线图对应的高度场Fig.1 a: Contour map.b: The imagined sea surface corresponding to the contour map

计算机运算能力的提高和计算机图形技术的进步为制作三维浪预报提供了有利条件。迄今国内外关于水波虚拟仿真已有大量研究，可以制作出逼真的海浪动画。海浪动画的制作过程一般包括两个主要步骤，首先建立波面高度场，然后对波面高度场进行渲染并生成三维海浪动画。前人根据不同应用目的开发了多种生成波面高度场的方法。其中，数值求解Navier-Stokes方程[1,2]的方法得到的波面极为逼真，但计算成本过高难以实现实时仿真；实际应用中多采用基于海浪模型的波面仿真[3]，这种方法计算量小同时画面的逼真性较高；基于快速傅立叶变换对给定靶谱的波面仿真是目前的主流方法[4-9]；另外还有基于Gerstner波动模型[3,10]或分形模型[11,12]的仿真方法。得到波面高度场之后需要对其渲染，该步骤非常复杂，结合了材质、用光、运动渲染、图像显示和动画演示等操作，前人针对不同的应用需求也有大量研究[3]在此不作赘述。目前一些先进的商业动画软件如Autodesk公司出品的MAYA和3Ds Max提供了强大的动画制作功能在上述动画开发环境中可以完成波面建模、渲染、动画制作一系列工作并最终形成动画。例如3Ds Max的工具包Dreamscape专门提供了制作海洋和山脉等自然地貌动画的功能。在要求不高的场合也可以在3Ds Max中定义“平面”基本体并对其添加“噪波”修改器，然后指定噪声强度来控制波面的振幅，设定波动传播速度及方向，然后对波面模型作光照调整、材质贴图等渲染得到波面的动画[13]。而三维动画海浪预报产品需要按波高分布情况构造包含多种波浪状态（多种平均波高和波长）的波面，即需要将不同平均波高和波长的波面拼合成一个波面。前人研究中波面模型大多是一个波面对应着一种波动状态（一种平均波高和波长），鄢来斌等[14]和吴晶等[15]为高效地模拟大面积海浪或海浪的变化过程（例如由于人眼观看越远处物体越模糊，将由近及远的波面以分辨率逐渐降低的波面表示并将不同分辨率的波面拼接到一起以节省计算量）考虑将不同的波面拼接到一起。前人研究未涉及如何按波高等值线划定区域将不同波高和波长的波面拼接成一个波面，为此需要专门设计一种波面模型的算法。渲染波面模型时若能采用商业软件处理将取得更好的效果，但是商业软件价格昂贵并且计算量较大，不能满足业务化预报对时效性的要求，因此需要自行开发一个简单实用的渲染器。

为适应海浪业务化预报工作的需求，本文基于线性海浪模型研究了三维动画海浪预报产品的制作方法。首先由方向谱生成一系列随时间变化的线性波面高度场，并计算波面高度场对应的反射系数（即制作一个简单渲染器）。然后按波高等值线图将平均波高和波长不同的波面高度场拼接成复合波面高度场，同样拼接出复合的反射系数场。最后用一般的绘图软件生成波面高度场的图像（即本文制作的渲染器的渲染结果），并将渲染好的波面时间序列串接成动画。2009年“芭玛”台风预警期间作者根据上述方法制作了首个三维动画海浪警报产品，该警报产品能以直观生动的形式将海浪预警信息传达给广大用户和公众，这对海洋预报工作具有重要意义。

1 波面高度场的理论模型

本节基于线性海浪模型采用线性滤波法[16]计算二维波面高度场。Longuet-Higgins将海浪波面高度ζ看作大量余弦波叠加构成的[17]

其中f0为谱峰频率，g为重力加速度。α、σ和γ分别为JONSWAP谱的强度因子、宽度因子和升高因子（取值范围参考[1]）。方向函数采用cos2θ，其中θ为方向角。得到方向频率谱

式中：F-1表示傅立叶逆变换，此处调用二维快速傅立叶逆变换（IFFT）。得到某一时刻的波面高度场。

图2 以JONSWAP谱为靶谱生成的波高分别为4 m和3 m的波面高度场Fig.2 Wave surface elevation field with mean wave height of 4m and 3m generated from JONSWAP spectrum

图2为某时刻t的波面高度场。为制作动画需要随时间演化的波面高度场，即波面高度场的时间序列ζ(tn)（tn= t+ndt，n = 1, 2, 3…N）。仅需要计算(8)式时先求出各个组成波相位随时间的变化再进行傅立叶逆变换即可得到波面高度场随时间的演化。由SW计算各个分量的相位作为初始相位

相当于给调制的噪声谱S.SW乘上exp(-i2π fndt)。按上述方法即可得到各种平均波高的波面高度场时间序列。

2 制作复合波面高度场

为反映出大浪的空间分布情况，需要将不同波高和波长的波浪场按波高等值线图拼接起来。这里需要解决两个问题：第一，将海浪预报图转化成高程数据，根据高程来确定各预报区域所使用的平均波高的波面场；第二，不同波面高度场之间存在的相位差会导致合成波面不连续，因此需要设计一组平滑函数使不同波面高度场平滑连接。对第一个问题，需要根据海浪预报图的存储格式来确定相应算法，本文研究中海浪预报图为位图格式（Windows 256色bmp文件），只需将色标数据转化成数字高程即可。对第二个问题，平滑函数应采用钟形曲面，该曲面的水平截面形状应和等值线的形状一致。由于等高范围可能为多连通区域，因此难以给出解析形式的平滑函数[19]，但可以从浪区边界开始逐层向内（外）搜索，并为各层数据点递增（减）赋值，这样赋值相当于构造窗函数使不同波高波浪场平滑连接在一起，本文研究发现采用正弦或指数形式的函数平滑效果较好，采用线性函数时波高等值线附近区域波面高度场和反光系数场的渐变性较差。图3为边界点算法的示意图，其中灰色和白色表示波高不同的两个区域（如灰色和白色分别表示3 m和4 m的浪区），每个小方格(i, j)表示预报图中一个像素。由(i, j)点周围8个点的色阶取值（本文所用256色BMP文件包含256个色阶，每个色阶分别以整数0 ~ 255表示）判断该点是否为边界点（若某点色阶值与其周围8个点的色阶值相等则说明该点为区域内点；反之该点为区域的边界点），容易判断出下标为a的数据点为白色区域的边界点，将a层颜色更新为灰色同样算法可以得到a层内一层（下标为b）。然后对a，b，c，…层分别赋(sin(Pα)+1)/2，其中α=π/2N，P=0,1,2,3,…,N；第N层以内的点都赋1；对a层外部的a’，b’，c’，…层分别赋(sin(-Pα)+1)/2，其中α=π/2N，P=1,2,3,…,N，从而构造出图3中白色 区 域 的 平 滑 函 数 ， 同 样 道 理可以构造灰色区域的平滑函数。按上述方法得到图3中黑色粗线所标断面的平滑函数如图4

图3 边界点搜索算法示意图。白色和灰色区域分别表示波高不同的两个海域，在预报图中为不同颜色的两个海域（色阶取值不同），图中每个小方格(i,j)表示预报图的一个像素。其中黑色粗线表示穿过海域的一条断面Fig.3 Boundary points algorithm.The white and gray areas represent sea areas with different wave heights in wave forecasting map.Square (i,j) represents a pixel in the forecasting map.The bold line represents a section across the sea area

图4 图3中黑色粗线标示断面对应的平滑函数。Fig.4 Smoothing functions for the section labeled by bold line in Fig.3.

下面以国家海洋预报中心2009年10月11日针对台风“芭玛”发布的海浪黄色警报图（图5）为例构造平滑函数（图6）。用图6所示平滑函数乘以图2所示波面得到平滑后的 4m和 3m浪区（图7）。同样办法处理波高小于3米的区域。将平滑后的浪区叠加就得到复合的波面高度场。下一节介绍复合波浪场的渲染方法。

图5 2009年10月11日08时发布“芭玛”台风海浪黄色警报图Fig.5 Yellow warning of typhoon ‘Parma’.UTC00 Oct.11 2009

图6 4m和3m浪区对应的平滑函数Fig.6 Smoothing function for sea areas of 4 m and 3 m significant wave heights

图7 平滑后4m和3m波高的浪区Fig.7 Smoothed wave fields with 4m and 3m significant wave heights.

3 渲染波面光学效果

图8 光线在波面的反射与折射。其中ζ表示波面，L1、L2和L3分别表示第一处光线入射点处的入射、反射和折射光线，θ i1、θ r2、和θt3分别表示该点入射、反射和折射角， 和 分别表示该点处波面的法向和切向。L2作为入射光线第二次反射和折射的物理量以下标2表示。Fig.8 Reflections and refractions at the wave surface.ζ is the wave surface.L1, L2and L3represent the first incident light, the first reflection and

波面高度场反射系数的变化是产生水体波动感觉的重要原因，而波面反射系数直接受波面坡度影响[3]。本节主要研究波面坡度对波面反射系数的影响。事实上要渲染出逼真的水波效果面临很多非常复杂的问题，如水波的颜色由入射光的颜色（如天空的颜色）、水体本身的颜色、水底反射光的颜色等多种因素决定，又如光线在波面产生的多次反射折射的效果等[5]。 the first refraction respectively.θi1, θr2and θt3represent the first incident angle, the first reflective angle and the first refractive angle.andrepresent the normal and tangential vectors of the first incident light respectively.Symbols with subscript 2 are variables corresponding to the second process of reflection and refraction, in which case L2is the incident light

图8为光线在随机波面发生的反射和折射过程。光线L1入射到波面ζ上发生反射，反射光线L2再次射入波面ζ产生二次反射光线L5，两次折射过程中的折射光线分别为L3和L4。第2节生成波面高度场坡度变化非常复杂，考虑到计算成本，渲染波面时仅考虑第一次反射和折射的影响。下面简要介绍波面反射系数场的计算方法。

三维波面高度场中水平位置x，高度ζ (x,t)的位置矢量可以表示为

根据该式对波面进行渲染。

4 一个实例

国家海洋预报中心2009年10月11日针对台风“芭玛”发布了海浪黄色警报，此次预警过程中作者制作了首个三维动画海浪警报产品。图9为三维海浪动画预报中的一帧画面。图中黄色曲线表示有效波高4米的区域，蓝色曲线表示有效波高3 m的区域，其余为有效波高小于3 m的区域，与图5中等值线完全一致。容易看出预报范围内波高 3 m以上的浪区中波高和波长明显大于周边区域。从动画预报中可以看到不同波高区域中波动的传播过程（动画产品为GIF文件,可由作者提供）。由此可见三维动画海浪预报不仅具有波高等值线预报图的功能，而且具有直观生动的优点，便于一般用户了解海浪预报的内容。另外，在PC机上制作60帧画面的动画耗时约10 min至15 min，能够满足业务化预报对快捷性的要求。

图9 由警报图5制作的三维海浪动画的一帧画面。黄线和蓝线分别表示4m和3m有效波高等值线Fig.9 A frame of the 3-D animation corresponding to the yellow warning in Fig.5.The yellow and blue curves represent the contours of 4m and 3m respectively

5 结 论

本文基于线性海浪模型研究了三维动画海浪预报产品的制作方法，基本步骤为：以JONSWAP谱为靶谱生成一系列不同波高和波长的二维海面高度场，以及这些波面高度场对应的反光系数场；根据海浪有效波高等值线预报图设计出平滑函数，将不同波高的波面高度场平滑连接成复合波面高度场，同理连接出与之对应的复合反光系数场；由复合的波面高度场和反光系数场渲染三维波面，并将渲染后的波面时间序列串接成动画即可形成最终产品。三维动画形式的海浪预报产品能有效地克服传统海浪预报产品信息不直观的缺点，预报信息容易被广大用户和公众理解。

致谢：感谢审稿专家和中国海洋大学苗春葆博士以及华东师范大学河口海岸学国家重点实验室宗海波老师提出的宝贵意见。

[1]Foster N, D Metaxas.Realistic animation of liquids [J].Graphical models and image processing, 1996, 58(5): 471-483.

[2]Kass R , M G.Rapid.Stable fluid dynamics for computer graphics [J].Computer Graphics, 1990, 24(4): 49-57.

[3]Tessendorf J.Simulating Ocean water.http: //graphics.ucsd.Edu / courses/rendering/2005/jdewall/tessendorf.pdf [OL].2004.

[4]Joe S.A simple fluid solver based on the FFT [J].Journal of Graphics Tools, 2001, 6(2): 43-52.

[5]Mastin A G, A P Watterberg, F J Mareda.Fourier synthesis of ocean scenes [J].IEEE Computer Graphics and Applications, 1987: 16-23.

[6]宋利, 周源华, 周军.基于 Fourier合成技术的水波动画 [J].计算机工程与应用, 2002, 6: 10-11.

[7]皮学贤, 杨旭东, 李思昆.近岸水域的波浪与水面仿真 [J].计算机学报, 2007, 30(2): 324-329.

[8]陈文辉, 谈晓军.大规模流域内水体三维仿真研究 [J].系统仿真学报, 2004, 16(11): 2 409-2 412.

[9]徐利明, 姜昱明.基于谱分析的实时波浪模拟 [J].系统仿真学报, 2005, 17(9): 2 092-2 095.

[10]Damien H, F Neyret, M Cani.Interactive animation of ocean waves [C].in Proceedings of the 2002 ACM SIGGRAPH.2002.

[11]Perlin K, E M.Hoffert.Hyper texture [C].in Proceedings of the 1989 ACM SIGGRAPH Conference.1989.New York: ACM.

[12]李广鑫, 丁振国, 詹海生, 等.一种面向虚拟环境的真实感水波面建模算法 [J].计算机研究与发展, 2004, 41(9): 1 580-1 585.

[13]王琦.Autodesk 3ds Max 2010标准培训教材II [M].北京: 人民邮电出版社, 2010.

[14]鄢来斌, 李思昆, 张秀山.虚拟海战场景中的海浪实时建模与绘制技术研究 [J].计算机研究与发展, 2001, 38(5): 568-573.

[15]吴晶, 徐晓刚.大规模海面的快速漫游 [J].系统仿真学报, 2006, 18(9): 2 505-2 507.

[16]文圣常, 余宙文.海浪理论与计算原理 [M].北京: 科学出版社, 1984.

[17]Longuet-Higgins M S .The statistical distribution of the heights of sea waves [J].J.mar.Res, 1952, 11: 245-266.

[18]Hasselmann K.Measurements of wind wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP) [R].1973, Erganzungsbeft Zeit.Deut.Hydr.Zeit.No.12.99.

[19]张洪刚, 陈光, 郭军.图像处理与识别 [M].北京: 北京邮电大学出版社, 2006.

Preparing 3-D animation of sea wave forecasting

Gao Zhi-yi, Yu Fu-jiang, Xu Fu-xiang

(National Marine Environmental Forecasting Center, Beijing 100081, China)

P731.33

A

1001-6932(2011)02-0173-07

2010-4-22；收修改稿日期：2010-6-17

十一五国家科技支撑计划：重大海洋灾害预警及应急技术研究（2006BAC0300）

高志一（1980-），男，博士，研究实习员，主要从事小尺度海气相互作用与海浪理论研究。电子邮箱：zhiyi_gao@yahoo.com.cn。