毛科峰, 陳 希, 李 妍, 蕭中樂

（1. 解放軍理工大學(xué)氣象學(xué)院, 江蘇 南京 211101; 2. 解放軍96631部隊, 北京 102208）

琉球群島海域海浪數(shù)值計算地形處理效應(yīng)及試驗分析

毛科峰1, 陳 希1, 李 妍1, 蕭中樂2

（1. 解放軍理工大學(xué)氣象學(xué)院, 江蘇 南京 211101; 2. 解放軍96631部隊, 北京 102208）

針對琉球群島海域內(nèi)多島嶼復(fù)雜地形對海浪數(shù)值計算的特殊影響, 發(fā)展一種綜合利用水深數(shù)據(jù)和高分辨率海岸線數(shù)據(jù)優(yōu)化計算網(wǎng)格且引入次網(wǎng)格地形效應(yīng)的方法, 并利用 WAVEWATCH-Ⅲ模式進行連續(xù)1個月的數(shù)值模擬試驗。結(jié)果表明: 采用該方法后, 在計算網(wǎng)格分辨率上, 充分考慮了海岸和多島嶼地形對海浪傳播的作用和多島嶼的次網(wǎng)格效應(yīng), 數(shù)值計算結(jié)果有明顯改善。

海浪; 琉球群島; 次網(wǎng)格; 數(shù)值模擬

琉球群島海域, 是大連、青島、廈門、釜山、仁川等港口東出太平洋的必經(jīng)之路, 也是太平洋東岸及大洋洲各國出入東海和黃海的重要水運航道, 在國民經(jīng)濟和軍事方面具有重要戰(zhàn)略意義; 從東南陸架淺海經(jīng)沖繩海槽過渡到東部西北太平洋, 日本九州至琉球、臺灣一線長約600 n mile的水域內(nèi), 有眾多的海峽、水道與太平洋溝通, 這里島嶼多、島礁多、海況各異、地形復(fù)雜。深淺不一的海區(qū)和水道里, 影響海浪的因素很多, 除了風(fēng)這一產(chǎn)生海浪的主要動力機制外, 波高的大小還取決于上風(fēng)向水域的長度和寬度、海底的地形、水位的深淺等諸多因素的影響。因此該海域的海浪數(shù)值計算有一定的特殊性,同時對海浪預(yù)報、海洋工程等方面有重要的現(xiàn)實意義。

目前基于組成波譜能量平衡方程（如式（1））的海浪計算模型發(fā)展到了以WAM[1], SWAN[2],

WAVEWATCH-Ⅲ[3]為代表的第三代海浪數(shù)值模式。在這個能量平衡方程中N即是波作用密度譜, 它是頻率f、傳播方向θ、時間t和地理空間位置的函數(shù),cg為群速, S是能量源匯, 包括風(fēng)攝入波動能量、波浪破碎和白帽破碎時能量的耗散、波與波之間的非線性相互作用引起的渦動黏滯、當水深較淺時底摩擦作用等能量的交換過程。基于（1）式進行數(shù)值計算時, 對譜空間（f,θ）或者（k,θ）和地理空間的計算分辨率選擇是十分重要的問題, 它對海浪計算的精度和效率有較大影響[4]。為了在琉球群島海域內(nèi)更準確地計算海浪, 必須提高海浪模式的空間網(wǎng)格分辨率來反映島嶼地形分布、海底地形及水深變化, 但是分辨率的提高又受到水深地形資料、風(fēng)場資料的分辨率以及計算效率的限制。Tolman和Chawla[4]曾指出某一空間分辨率計算網(wǎng)格未體現(xiàn)出的島嶼群和障礙物是海浪計算和預(yù)報誤差的主要來源之一。為了解決這一問題, Holthuijsen[2]提出了在SWAN模式中考慮次網(wǎng)格障礙物對海浪能量傳播的抑制作用的思想,Tolman[5]將這種方法用在 WAVEWATCH-Ⅲ模式中,考慮島嶼群和局部海冰的次網(wǎng)格效應(yīng)對海浪傳播的影響。但是如何科學(xué)方便地在海浪數(shù)值計算中處理多島嶼的分布和復(fù)雜海岸線等特征, 如何準確地設(shè)計次網(wǎng)格的計算方案還值得深入研究。鑒于此, 本文發(fā)展一種在群島海域內(nèi)綜合利用水深數(shù)據(jù)和高分辨率海岸線數(shù)據(jù)優(yōu)化海浪數(shù)值計算網(wǎng)格且引入島嶼次網(wǎng)格地形效應(yīng)的方法, 利用該方法和WAVEWATCH-Ⅲ模式進行了數(shù)值試驗并利用實測資料分析了數(shù)值試驗結(jié)果。

1 海浪數(shù)值模式和數(shù)值計算方案

琉球群島海域內(nèi)的琉球海脊將東海大陸架東側(cè)的深海區(qū)分割出來, 在東海大陸架和琉球海脊之間,形成狹長的深海區(qū)域——沖繩海槽。琉球海脊在600 m以淺主要以西表—石垣—宮古、沖繩島、奄美大島的島鏈形式存在。在600 m以深, 幾個島鏈相連,除了沖繩島以南的慶良間水道以外, 形成了完整的東北—西南向的海脊, 如圖1所示, 圖中等值線為水深值分布, 單位為 m。根據(jù)琉球海域的特點, 利用WAVEWATCH-Ⅲ模式進行海浪模擬是合適的。WAVEWATCH-Ⅲ模式合理地考慮了風(fēng)浪相互作用、非線性相互作用、耗散及底摩擦等作用, 能比較準確地模擬復(fù)雜的潮流、地形、風(fēng)場環(huán)境下的波浪場, 該模式是以上文的（1）式組成波譜能量平衡方程為基礎(chǔ),球坐標系下該方程可表式為:

式（2）中k為波數(shù), θ為波向, U為平均流速（對水深、時間平均）,為群速, σ為相對頻率,區(qū)別于絕對頻率是經(jīng)度, 是緯度, R為地球半徑, d為水深,是沿緯向、經(jīng)向的流速, S是源匯項。源匯項S包含風(fēng)攝入波動能量項Sin、非線性的波-波相互作用項Snl, 耗散項（白浪效應(yīng)）Sds和底摩擦項Sbot; 即通常情況下, 源匯項S可表示為:

對（2）式進行數(shù)值求解時采用分裂算法, 首先考慮水深在時間上的變化以及對應(yīng)波數(shù)網(wǎng)格上的變化; 這樣不考慮（暫時的）水面變化的影響后, 波數(shù)網(wǎng)格就是不變的, 水深也是準穩(wěn)定的; 然后分步計算（2）式左邊的海浪譜在物理空間上的傳播、波數(shù)空間上的傳播項和（2）式右邊的源函數(shù)項。在球坐標下, 將海浪在地理空間傳播表示為:

本文在 WAVEWATCH-Ⅲ中采用 ULTIMATE QUICKEST數(shù)值格式[6], 這個格式把緯向傳播和徑向傳播分開進行計算, 計算順序是任意, 在 方向上,第i和i－1這兩個格點之間的通量為Fi,?。

式中j, l和m分別是緯度λ, 譜空間θ和k方向上的離散網(wǎng)格計數(shù), n表示時間層, Cu是作用量密度的曲率, C是帶有符號的CFL數(shù),bφ˙代表兩個格點之間格點單元邊界上的傳播速度,iφ˙代表格點上的傳播速度, Nb代表格點單元邊界上的作用量, Ni代表第 i個格點上的作用量。在時, 該格式可得到穩(wěn)定解。故而該格式在 方向上也可表示為:

更改下標和相應(yīng)的增量, 同樣可得 方向上的傳播格式。

圖1 東中國海及琉球群島海域地形Fig. 1 Bathymetry distribution in the China offshore area and Ryukyu area

2 群島海域內(nèi)復(fù)雜地形的處理方案

2.1 優(yōu)化海浪數(shù)值計算網(wǎng)格的方法

WAVEWATCH-Ⅲ的數(shù)值計算網(wǎng)格利用“干濕”屬性來區(qū)分網(wǎng)格點是陸地或是海洋, 海洋上的“濕”網(wǎng)格是有效計算網(wǎng)格, 陸地“干”網(wǎng)格是計算的邊界。如果計算網(wǎng)格中, 有部分網(wǎng)格的“干濕”屬性是虛假的, 那么海浪數(shù)值模擬時必然會造成誤差。本文在設(shè)計琉球群島海域海浪計算網(wǎng)格時, 利用海圖中存儲的一些海岸線、島嶼、數(shù)據(jù), 對網(wǎng)格點進行海陸屬性判別, 優(yōu)化計算網(wǎng)格, 具體方法是:第一步, 將大陸或者島嶼當作閉合的任意多邊形, 得到閉合多邊形數(shù)據(jù); 第二步, 通過判斷計算網(wǎng)格點與這些多邊形的位置關(guān)系, 來區(qū)分每個網(wǎng)格點的海陸屬性, 在多邊形內(nèi)部的格點判別為陸地, 多邊形外的格點判別為海洋。本文引入了一種高效的判斷點與多邊形位置關(guān)系的算法[7]即夾角和法: 將整個大陸和任意島嶼看成一個邊數(shù)為 n的多邊形,其頂點序列為,海浪數(shù)值計算網(wǎng)格作為待判別的點為連接 P和多邊形的各個頂點, 計算其夾角和, 其中順時針向為正, 逆時針方向為負, 如圖2所示。

圖2 點與多邊形的關(guān)系Fig. 2 Ubiety between points and polygons

經(jīng)過這兩個步驟優(yōu)化后的計算網(wǎng)格融合了島嶼和海岸線信息, 判別了計算格點的“干濕”屬性即陸地和海洋屬性, 如圖2所示, 在島嶼閉合多邊形內(nèi)的計算點被判別為“干”屬性點, 之外的為“濕”屬性點。然后在海上的“濕”網(wǎng)格上進行水深插值, 獲得計算網(wǎng)格點上的水深值。

2.2 次網(wǎng)格地形效應(yīng)的計算方案

經(jīng)過優(yōu)化的海浪數(shù)值計算網(wǎng)格, 能夠在計算網(wǎng)格分辨率的水平上充分表現(xiàn)多島嶼的分布和海岸的地形。如果僅采用這種方法, 隨著計算網(wǎng)格分辨率的提高, 能夠描述的島嶼地形分布就越細致, 但是計算網(wǎng)格的分辨率不可能無限提高。在一定計算分辨率水平下, 仍然有比計算網(wǎng)格尺度更小的島嶼等地形特征不能被描述, 如圖3所示有若干個類似圖中A點所指的小島嶼, 在這種計算分辨率水平下, 它周圍的計算格點都是“濕”屬性, 在計算時, 便把它“忽略”了。因此, 有必要在海浪數(shù)值計算中引入次網(wǎng)格地形效應(yīng), 充分考慮計算網(wǎng)格分辨所不能分辨的小島嶼的地形作用。

圖3 計算網(wǎng)格和島嶼地形分布Fig. 3 Island distribution and computing grid

在（8）式中若在第i個計算格點所在的單元格內(nèi),有一個未被網(wǎng)格分辨的小島嶼或障礙物, 如圖 4, 為了在這個計算單元格內(nèi)考慮它對海浪傳播的阻礙作用, 在計算單元邊界的兩個方向上定義能通量穿透度系數(shù)系數(shù)的值域為: 0～1, 取0時表示次網(wǎng)格島嶼的影響作用為封閉邊界, 在這個計算單元的海浪能通量被完全阻礙, 取 1時表示沒有次網(wǎng)格島嶼。于是（8）式可以寫為:

圖4 計算網(wǎng)格和次網(wǎng)格尺度的島嶼Fig. 4 Computation grid and sub-grid island

通過次網(wǎng)格能通量穿透度系數(shù)的引入來抑制計算格點單元之間的能量通量, 體現(xiàn)了次網(wǎng)格尺度島嶼對海浪的阻礙作用。在優(yōu)化過的計算網(wǎng)格上, 利用上述判斷點與多邊形位置關(guān)系的夾角和法, 判別計算網(wǎng)格的“濕”屬性海洋計算點和未被計算網(wǎng)格分辨的小島嶼的位置關(guān)系, 搜索出需要計算能通量穿透度系數(shù)的計算格點單元, 然后求計算格點單元的穿透度系數(shù)如圖 4, 在 方向上在?L為 方向島嶼的寬度,? 為 方向格點單元的寬度, ?Lφ為 方向島嶼的寬度, ? 為 方向格點單元的寬度。

3 數(shù)值試驗與分析

3.1 試驗方案

對琉球群島海域 2004年 10月的海浪過程進行模擬, 并與實測資料對比分析: 計算范圍為 20°～35°N, 116°～132.75°E, 格距為 15′; 模式包括風(fēng)浪相互作用、波波相互作用、白帽耗散、底摩擦作用等物理過程。模擬時段從2004年10月1日02時以120 s的時間步長積分到的10月31日23時。由于該時段內(nèi)有 3次臺風(fēng)過程影響該海域, 故試驗所用的風(fēng)場根據(jù)Holland[8]臺風(fēng)風(fēng)場模型和QuickSCAT/NCEP混合風(fēng)場[9]合成。試驗1采用本文提出的綜合利用水深數(shù)據(jù)和高分辨率海岸線數(shù)據(jù)優(yōu)化計算網(wǎng)格且引入次網(wǎng)格地形效應(yīng)的方法進行模擬, 試驗 2直接用ETOPO2水深資料并且不考慮次網(wǎng)格地形效應(yīng)進行模擬。圖5中給出了琉球群島海域附近的島嶼分布,圓點表示計算網(wǎng)格的陸地點, “×”點表示次網(wǎng)格點,正方形標出的位置為海浪測站的3個站點: NAHA站、NAKA站和NASE站。

圖5 琉球群島海域陸地點和次網(wǎng)格點Fig. 5 Land points and sub-grid points in the Ryukyu offshore area

3.2 結(jié)果分析

這 3個測站的海浪有效波高每兩小時有一個觀測值, 分別將試驗1和試驗2三個測站2004年10月的有效波高模擬結(jié)果的誤差進行統(tǒng)計, 結(jié)果如表 1所示, 試驗1中有效波高的偏度、均方根誤差、平均絕對誤差、均方相對誤差、平均相對誤差均明顯小于試驗2, 其中試驗1較試驗2的偏度降低了0.31 m,平均相對誤差降低約 12%。這說明在綜合利用水深數(shù)據(jù)和高分辨率海岸線數(shù)據(jù)優(yōu)化計算網(wǎng)格和在海浪數(shù)值計算中引入次網(wǎng)格地形效應(yīng)后海浪有效波高的模擬的精度有了較大提高了。這在圖6中 3個測站10月的有效波高模擬值與實測值的散點圖也能得到證明。圖7是兩個試驗方案的模擬結(jié)果和3個浮標資料的時間序列比較圖, 圖 7表明在 NAHA站和NASE站具有同樣的特點: 試驗2的有效波高值較浮標觀測值明顯偏大, 試驗1的較試驗2的值偏小且更接近浮標觀測值, 可見在這兩個浮標站, 試驗1結(jié)果明顯優(yōu)于試驗2; 在NAKA站, 試驗1和2模擬的有效波高值幾乎完全一致, 兩個試驗方案模擬結(jié)果沒有差別, 而且與浮標觀測值十分接近, 有效波高模擬誤差很小。究其原因: 如圖6試驗2所示, NAHA站和 NASE站附近均有若干個小島嶼, 并未被計算

網(wǎng)格分辨為陸地點, 故而在試驗 2中這是導(dǎo)致NAHA站和NASE站有效波高模擬誤差的主要原因之一。在試驗1中引入了島嶼次網(wǎng)格地形效應(yīng), 考慮了這些島嶼存在對海浪的影響, 克服了有效波高模擬偏大的誤差, 這說明本文的次網(wǎng)格作用計算方案能夠體現(xiàn)出比計算網(wǎng)格尺度小的島嶼對海浪傳播明顯的抑制作用。兩種試驗方案結(jié)果的差別在NAKA站卻不存在, 是因為該站附近沒有其他的小島嶼, 因而該點附近也沒有次網(wǎng)格效應(yīng)的影響。如圖 8是兩種試驗方案在計算區(qū)域內(nèi)10月份有效波高的月平均值之差的分布圖, 試驗1與試驗2相比, 具有明顯的負系統(tǒng)偏差, 月平均有效波高之差的最大值達0.7 m,同時與圖 5對比可以發(fā)現(xiàn)月平均有效波高有明顯差異的地方就集中在考慮次網(wǎng)格效應(yīng)的計算網(wǎng)格點附近。因此可以說, 島嶼次網(wǎng)格效應(yīng)對有效波高的模擬誤差有明顯的影響, 但是這種影響也是有局地性的,這可能是島嶼次網(wǎng)格效應(yīng)對涌浪傳播的阻礙和抑制作用導(dǎo)致的。

表1 兩次試驗的誤差Tab. 1 Errors between observed data and computation results for two experiments

圖6 3個測站波高模擬值與實測值的散點圖Fig. 6 Scatter plots of wave heights of observed data at three sites and computation results computed for two experiments

圖7 3個測站10月份的有效波高模擬值與實測值的比較Fig. 7 Comparison of wave heights between model results and observations in October at three sites

圖8 兩種模擬方案10月份平均有效波高的差值分布Fig. 8 Mean difference maps of the wave heights for the month of October between two simulation schemes

4 結(jié)論

1）綜合利用水深數(shù)據(jù)和高分辨率海岸線數(shù)據(jù)來判別海浪數(shù)值計算網(wǎng)格點的海陸類型, 使計算網(wǎng)格設(shè)計更加優(yōu)化合理, 該方法能夠在計算網(wǎng)格的分辨率上充分體現(xiàn)多島嶼的分布特征。

2）在琉球群島海域內(nèi), 對多島嶼復(fù)雜地形是海浪模擬的主要誤差源; 引入島嶼次網(wǎng)格作用計算方案后, 能夠體現(xiàn)出比計算網(wǎng)格尺度小的島嶼對海浪傳播有明顯的抑制作用, 可以提高海浪模擬的精度,通過與浮標實測資料的對比表明, 有效波高模擬的平均相對誤差降低約12%。

3）島嶼次網(wǎng)格效應(yīng)對有效波高的模擬誤差有明顯的影響, 但影響是局地性的, 這可能是島嶼次網(wǎng)格效應(yīng)對涌浪傳播的阻礙和抑制作用導(dǎo)致的。

4）通過琉球群島海域海浪模擬數(shù)值試驗表明本文發(fā)展的優(yōu)化海浪數(shù)值計算網(wǎng)格且引入次網(wǎng)格地形效應(yīng)的方法具有實用性和有效性, 這種方法, 為多島嶼復(fù)雜地形條件下海浪數(shù)值研究打下基礎(chǔ)。

[1] WAMDI Group. The WAM model-A third generation ocean wave prediction model[J]. J Phys Oceanography, 1988, 18（8）: 1 775-1 810.

[2] Holthuijsen L H. SWAN III version 40.45 USER MANUAL[R]. Delft of Netherlands: Delft University of Technology, 2006.

[3] Tolman H L. User manual and system documentation of WAVEWATCH-III version 2.22[EB/OL]. http://NOAA /NWS / NCEP / MMAB Technical Note, 2002-12-18.

[4] Tolman H L. Alleviating the garden sprinkler effect in wind wave models[J]. Ocean Modelling, 2002, 4（3）:269-289.

[5] Hendrik L T. Treatment of unresolved islands and ice in wind wave models[J]. Ocean Modelling, 2003, 5（3）:219-231.

[6] Leonard B P. A stable and accurate convective modeling procedure based on quadratic upstream interpolation[J]. Comput Methods Appl Mech Engng, 1979,19（1）: 59-98.

[7] Feito F R, Tomes J C, Urena L A. Orientation simplicity and inclusion test for planar polygons[J].Com-puter&Graphics, 1995, 19（4）: 595-600.

[8] Holland G J. An Analytic model of the wind and pressure profiles in Hurricanes[J]. Mon Wea Rev, 1980,108（8）: 1 212-1 215.

[9] 李明悝, 侯一筠. 利用 QuikSCAT/NCEP混合風(fēng)場及WAVEWATCH模擬東中國海風(fēng)浪場[J].海洋科學(xué),2005, 29（6）: 9-12.

Received: Mar., 19, 2009

Key words:ocean wave; Ryukyu area; sub-grid; numerical modeling

Abstract:In order to describe complex multi-island terrain and coastlines adequately for modeling ocean wave in sea areas by Ryukyu, an algorithm based on high-resolution bathymetry and coastline data was developed to optimize the design of wave computation grid and introduce sub-grid terrain effect into wave computation.WAVEWATCH-Ⅲ was used to carry on continuous numerical simulation for a month. It was validated that the impact of complex coastline and multi-island on the wave propagation, computation grid resolution, and sub-grid effect of the multi-island effect was adequately address. Application of this algorithm led to much improved computational results.

（本文編輯: 劉珊珊）

Numerical simulation for ocean wave in sea areas by Ryukyu

MAO Ke-feng1, CHEN Xi1, LI Yan1, XIAO Zhong-le2
（1. Institute of Meteorology. PLA University of Science and Technology, Nanjing 211101,China;2. No.96631 army of PLA,Beijing 102208,China）

P731

A

1000-3096（2010）11-0091-06

2009-03-19;

2010-07-10

軍隊重點基金; 衛(wèi)星海洋環(huán)境動力學(xué)國家重點實驗室開放研究基金（SOED1009）

毛科峰（1981-）, 男, 湖南常德人, 解放軍理工大學(xué)氣象學(xué)院博士研究生, 研究方向: 物理海洋, 電話: 025-80831609, E-mail: maomaopla@163.com