琉球群島海域海浪數值計算地形處理效應及試驗分析

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

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

琉球群島海域海浪數值計算地形處理效應及試驗分析

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

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

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

海浪; 琉球群島; 次網格; 數值模擬

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

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

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

1 海浪數值模式和數值計算方案

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

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

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

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

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

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

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

2 群島海域內復雜地形的處理方案

2.1 優化海浪數值計算網格的方法

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

圖2 點與多邊形的關系Fig. 2 Ubiety between points and polygons

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

2.2 次網格地形效應的計算方案

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

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

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

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

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

3 數值試驗與分析

3.1 試驗方案

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

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

3.2 結果分析

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

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

表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 結論

1）綜合利用水深數據和高分辨率海岸線數據來判別海浪數值計算網格點的海陸類型, 使計算網格設計更加優化合理, 該方法能夠在計算網格的分辨率上充分體現多島嶼的分布特征。

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

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

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

[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混合風場及WAVEWATCH模擬東中國海風浪場[J].海洋科學,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

軍隊重點基金; 衛星海洋環境動力學國家重點實驗室開放研究基金（SOED1009）

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