海浪預報三維動畫計算原理與制作方法

高志一，于福江，許富祥

（國家海洋環(huán)境預報中心 北京 100081）

海浪預報三維動畫計算原理與制作方法

高志一，于福江，許富祥

（國家海洋環(huán)境預報中心 北京 100081）

研究了基于線性海浪模型制作三維動畫海浪預報產(chǎn)品的制作方法，并制作了首個三維動畫海浪警報產(chǎn)品?；静襟E為：由線性海浪模型生成一系列不同波高和波長的二維海面高度場，以及這些波面高度場對應的反光系數(shù)場；根據(jù)海浪預報圖中波高的空間分布情況將不同波高的波面高度場拼接成復合的波面高度場，同時拼接出與之對應的復合反光系數(shù)場；由復合的波面高度場和反光系數(shù)場渲染出三維波面，并將渲染后的波面時間序列串接成動畫形成最終產(chǎn)品。三維動畫形式的海浪預報產(chǎn)品能有效地克服海浪有效波高等值線圖信息不直觀的缺點，預報信息易于被廣大用戶和公眾理解。

海浪預報；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

海浪預報服務最初是為滿足海上軍事活動的需要而出現(xiàn)的，迄今已有數(shù)十年的歷史。由于海浪的隨機特性，對它只能進行統(tǒng)計預報。海浪預報產(chǎn)品需要告訴用戶洋面或海區(qū)平均波高的情況。因此海浪預報產(chǎn)品一般采用有效波高等值線圖，該形式的產(chǎn)品實際上是一種傳統(tǒng)涉海行業(yè)中廣泛采用的專業(yè)用圖。然而隨著經(jīng)濟進步和人類海上活動日益頻繁，需要海浪預報服務的用戶群體早已突破傳統(tǒng)涉海行業(yè)。由于新興用戶可能具有不同的專業(yè)背景，對他們來說波高等值線圖不直觀甚至可能產(chǎn)生誤導。如圖1a為一個等高線示意圖，該圖的含義是顏色越淺的等值線圈出區(qū)域的有效波高越大，但一般用戶不容易想象出圖1a對應的浪場波動狀態(tài)，而是容易將它想象成圖1b所示的三維圖景，即“一座水山”。由此可見等值線形式的海浪預報圖已經(jīng)不能滿足面向一般用戶的預報服務需要，新的預報產(chǎn)品應該具有更直觀生動的形式。這就要求新的預報產(chǎn)品不僅具有傳統(tǒng)預報產(chǎn)品的功能，更要便于一般用戶了解海浪預報中的深層次信息，如等高線大值區(qū)的隨機波動具有較大的平均波高和波長。因此三維動畫形式是新型海浪預報產(chǎn)品首選形式，然而目前尚無此類產(chǎn)品面世。

圖1 a: 圖等高線圖。b: 圖為等高線圖對應的高度場Fig.1 a: Contour map.b: The imagined sea surface corresponding to the contour map

計算機運算能力的提高和計算機圖形技術(shù)的進步為制作三維浪預報提供了有利條件。迄今國內(nèi)外關(guān)于水波虛擬仿真已有大量研究，可以制作出逼真的海浪動畫。海浪動畫的制作過程一般包括兩個主要步驟，首先建立波面高度場，然后對波面高度場進行渲染并生成三維海浪動畫。前人根據(jù)不同應用目的開發(fā)了多種生成波面高度場的方法。其中，數(shù)值求解Navier-Stokes方程[1,2]的方法得到的波面極為逼真，但計算成本過高難以實現(xiàn)實時仿真；實際應用中多采用基于海浪模型的波面仿真[3]，這種方法計算量小同時畫面的逼真性較高；基于快速傅立葉變換對給定靶譜的波面仿真是目前的主流方法[4-9]；另外還有基于Gerstner波動模型[3,10]或分形模型[11,12]的仿真方法。得到波面高度場之后需要對其渲染，該步驟非常復雜，結(jié)合了材質(zhì)、用光、運動渲染、圖像顯示和動畫演示等操作，前人針對不同的應用需求也有大量研究[3]在此不作贅述。目前一些先進的商業(yè)動畫軟件如Autodesk公司出品的MAYA和3Ds Max提供了強大的動畫制作功能在上述動畫開發(fā)環(huán)境中可以完成波面建模、渲染、動畫制作一系列工作并最終形成動畫。例如3Ds Max的工具包Dreamscape專門提供了制作海洋和山脈等自然地貌動畫的功能。在要求不高的場合也可以在3Ds Max中定義“平面”基本體并對其添加“噪波”修改器，然后指定噪聲強度來控制波面的振幅，設定波動傳播速度及方向，然后對波面模型作光照調(diào)整、材質(zhì)貼圖等渲染得到波面的動畫[13]。而三維動畫海浪預報產(chǎn)品需要按波高分布情況構(gòu)造包含多種波浪狀態(tài)（多種平均波高和波長）的波面，即需要將不同平均波高和波長的波面拼合成一個波面。前人研究中波面模型大多是一個波面對應著一種波動狀態(tài)（一種平均波高和波長），鄢來斌等[14]和吳晶等[15]為高效地模擬大面積海浪或海浪的變化過程（例如由于人眼觀看越遠處物體越模糊，將由近及遠的波面以分辨率逐漸降低的波面表示并將不同分辨率的波面拼接到一起以節(jié)省計算量）考慮將不同的波面拼接到一起。前人研究未涉及如何按波高等值線劃定區(qū)域?qū)⒉煌ǜ吆筒ㄩL的波面拼接成一個波面，為此需要專門設計一種波面模型的算法。渲染波面模型時若能采用商業(yè)軟件處理將取得更好的效果，但是商業(yè)軟件價格昂貴并且計算量較大，不能滿足業(yè)務化預報對時效性的要求，因此需要自行開發(fā)一個簡單實用的渲染器。

為適應海浪業(yè)務化預報工作的需求，本文基于線性海浪模型研究了三維動畫海浪預報產(chǎn)品的制作方法。首先由方向譜生成一系列隨時間變化的線性波面高度場，并計算波面高度場對應的反射系數(shù)（即制作一個簡單渲染器）。然后按波高等值線圖將平均波高和波長不同的波面高度場拼接成復合波面高度場，同樣拼接出復合的反射系數(shù)場。最后用一般的繪圖軟件生成波面高度場的圖像（即本文制作的渲染器的渲染結(jié)果），并將渲染好的波面時間序列串接成動畫。2009年“芭瑪”臺風預警期間作者根據(jù)上述方法制作了首個三維動畫海浪警報產(chǎn)品，該警報產(chǎn)品能以直觀生動的形式將海浪預警信息傳達給廣大用戶和公眾，這對海洋預報工作具有重要意義。

1 波面高度場的理論模型

本節(jié)基于線性海浪模型采用線性濾波法[16]計算二維波面高度場。Longuet-Higgins將海浪波面高度ζ看作大量余弦波疊加構(gòu)成的[17]

其中f0為譜峰頻率，g為重力加速度。α、σ和γ分別為JONSWAP譜的強度因子、寬度因子和升高因子（取值范圍參考[1]）。方向函數(shù)采用cos2θ，其中θ為方向角。得到方向頻率譜

式中：F-1表示傅立葉逆變換，此處調(diào)用二維快速傅立葉逆變換（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計算各個分量的相位作為初始相位

相當于給調(diào)制的噪聲譜S.SW乘上exp(-i2π fndt)。按上述方法即可得到各種平均波高的波面高度場時間序列。

2 制作復合波面高度場

為反映出大浪的空間分布情況，需要將不同波高和波長的波浪場按波高等值線圖拼接起來。這里需要解決兩個問題：第一，將海浪預報圖轉(zhuǎn)化成高程數(shù)據(jù)，根據(jù)高程來確定各預報區(qū)域所使用的平均波高的波面場；第二，不同波面高度場之間存在的相位差會導致合成波面不連續(xù)，因此需要設計一組平滑函數(shù)使不同波面高度場平滑連接。對第一個問題，需要根據(jù)海浪預報圖的存儲格式來確定相應算法，本文研究中海浪預報圖為位圖格式（Windows 256色bmp文件），只需將色標數(shù)據(jù)轉(zhuǎn)化成數(shù)字高程即可。對第二個問題，平滑函數(shù)應采用鐘形曲面，該曲面的水平截面形狀應和等值線的形狀一致。由于等高范圍可能為多連通區(qū)域，因此難以給出解析形式的平滑函數(shù)[19]，但可以從浪區(qū)邊界開始逐層向內(nèi)（外）搜索，并為各層數(shù)據(jù)點遞增（減）賦值，這樣賦值相當于構(gòu)造窗函數(shù)使不同波高波浪場平滑連接在一起，本文研究發(fā)現(xiàn)采用正弦或指數(shù)形式的函數(shù)平滑效果較好，采用線性函數(shù)時波高等值線附近區(qū)域波面高度場和反光系數(shù)場的漸變性較差。圖3為邊界點算法的示意圖，其中灰色和白色表示波高不同的兩個區(qū)域（如灰色和白色分別表示3 m和4 m的浪區(qū)），每個小方格(i, j)表示預報圖中一個像素。由(i, j)點周圍8個點的色階取值（本文所用256色BMP文件包含256個色階，每個色階分別以整數(shù)0 ~ 255表示）判斷該點是否為邊界點（若某點色階值與其周圍8個點的色階值相等則說明該點為區(qū)域內(nèi)點；反之該點為區(qū)域的邊界點），容易判斷出下標為a的數(shù)據(jù)點為白色區(qū)域的邊界點，將a層顏色更新為灰色同樣算法可以得到a層內(nèi)一層（下標為b）。然后對a，b，c，…層分別賦(sin(Pα)+1)/2，其中α=π/2N，P=0,1,2,3,…,N；第N層以內(nèi)的點都賦1；對a層外部的a’，b’，c’，…層分別賦(sin(-Pα)+1)/2，其中α=π/2N，P=1,2,3,…,N，從而構(gòu)造出圖3中白色 區(qū) 域 的 平 滑 函 數(shù) ， 同 樣 道 理可以構(gòu)造灰色區(qū)域的平滑函數(shù)。按上述方法得到圖3中黑色粗線所標斷面的平滑函數(shù)如圖4

圖3 邊界點搜索算法示意圖。白色和灰色區(qū)域分別表示波高不同的兩個海域，在預報圖中為不同顏色的兩個海域（色階取值不同），圖中每個小方格(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中黑色粗線標示斷面對應的平滑函數(shù)。Fig.4 Smoothing functions for the section labeled by bold line in Fig.3.

下面以國家海洋預報中心2009年10月11日針對臺風“芭瑪”發(fā)布的海浪黃色警報圖（圖5）為例構(gòu)造平滑函數(shù)（圖6）。用圖6所示平滑函數(shù)乘以圖2所示波面得到平滑后的 4m和 3m浪區(qū)（圖7）。同樣辦法處理波高小于3米的區(qū)域。將平滑后的浪區(qū)疊加就得到復合的波面高度場。下一節(jié)介紹復合波浪場的渲染方法。

圖5 2009年10月11日08時發(fā)布“芭瑪”臺風海浪黃色警報圖Fig.5 Yellow warning of typhoon ‘Parma’.UTC00 Oct.11 2009

圖6 4m和3m浪區(qū)對應的平滑函數(shù)Fig.6 Smoothing function for sea areas of 4 m and 3 m significant wave heights

圖7 平滑后4m和3m波高的浪區(qū)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

波面高度場反射系數(shù)的變化是產(chǎn)生水體波動感覺的重要原因，而波面反射系數(shù)直接受波面坡度影響[3]。本節(jié)主要研究波面坡度對波面反射系數(shù)的影響。事實上要渲染出逼真的水波效果面臨很多非常復雜的問題，如水波的顏色由入射光的顏色（如天空的顏色）、水體本身的顏色、水底反射光的顏色等多種因素決定，又如光線在波面產(chǎn)生的多次反射折射的效果等[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為光線在隨機波面發(fā)生的反射和折射過程。光線L1入射到波面ζ上發(fā)生反射，反射光線L2再次射入波面ζ產(chǎn)生二次反射光線L5，兩次折射過程中的折射光線分別為L3和L4。第2節(jié)生成波面高度場坡度變化非常復雜，考慮到計算成本，渲染波面時僅考慮第一次反射和折射的影響。下面簡要介紹波面反射系數(shù)場的計算方法。

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

根據(jù)該式對波面進行渲染。

4 一個實例

國家海洋預報中心2009年10月11日針對臺風“芭瑪”發(fā)布了海浪黃色警報，此次預警過程中作者制作了首個三維動畫海浪警報產(chǎn)品。圖9為三維海浪動畫預報中的一幀畫面。圖中黃色曲線表示有效波高4米的區(qū)域，藍色曲線表示有效波高3 m的區(qū)域，其余為有效波高小于3 m的區(qū)域，與圖5中等值線完全一致。容易看出預報范圍內(nèi)波高 3 m以上的浪區(qū)中波高和波長明顯大于周邊區(qū)域。從動畫預報中可以看到不同波高區(qū)域中波動的傳播過程（動畫產(chǎn)品為GIF文件,可由作者提供）。由此可見三維動畫海浪預報不僅具有波高等值線預報圖的功能，而且具有直觀生動的優(yōu)點，便于一般用戶了解海浪預報的內(nèi)容。另外，在PC機上制作60幀畫面的動畫耗時約10 min至15 min，能夠滿足業(yè)務化預報對快捷性的要求。

圖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 結(jié) 論

本文基于線性海浪模型研究了三維動畫海浪預報產(chǎn)品的制作方法，基本步驟為：以JONSWAP譜為靶譜生成一系列不同波高和波長的二維海面高度場，以及這些波面高度場對應的反光系數(shù)場；根據(jù)海浪有效波高等值線預報圖設計出平滑函數(shù)，將不同波高的波面高度場平滑連接成復合波面高度場，同理連接出與之對應的復合反光系數(shù)場；由復合的波面高度場和反光系數(shù)場渲染三維波面，并將渲染后的波面時間序列串接成動畫即可形成最終產(chǎn)品。三維動畫形式的海浪預報產(chǎn)品能有效地克服傳統(tǒng)海浪預報產(chǎn)品信息不直觀的缺點，預報信息容易被廣大用戶和公眾理解。

致謝：感謝審稿專家和中國海洋大學苗春葆博士以及華東師范大學河口海岸學國家重點實驗室宗海波老師提出的寶貴意見。

[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合成技術(shù)的水波動畫 [J].計算機工程與應用, 2002, 6: 10-11.

[7]皮學賢, 楊旭東, 李思昆.近岸水域的波浪與水面仿真 [J].計算機學報, 2007, 30(2): 324-329.

[8]陳文輝, 談曉軍.大規(guī)模流域內(nèi)水體三維仿真研究 [J].系統(tǒng)仿真學報, 2004, 16(11): 2 409-2 412.

[9]徐利明, 姜昱明.基于譜分析的實時波浪模擬 [J].系統(tǒng)仿真學報, 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]李廣鑫, 丁振國, 詹海生, 等.一種面向虛擬環(huán)境的真實感水波面建模算法 [J].計算機研究與發(fā)展, 2004, 41(9): 1 580-1 585.

[13]王琦.Autodesk 3ds Max 2010標準培訓教材II [M].北京: 人民郵電出版社, 2010.

[14]鄢來斌, 李思昆, 張秀山.虛擬海戰(zhàn)場景中的海浪實時建模與繪制技術(shù)研究 [J].計算機研究與發(fā)展, 2001, 38(5): 568-573.

[15]吳晶, 徐曉剛.大規(guī)模海面的快速漫游 [J].系統(tǒng)仿真學報, 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

十一五國家科技支撐計劃：重大海洋災害預警及應急技術(shù)研究（2006BAC0300）

高志一（1980-），男，博士，研究實習員，主要從事小尺度海氣相互作用與海浪理論研究。電子郵箱：zhiyi_gao@yahoo.com.cn。