長江口海域風(fēng)浪場的計算與特征分析

耿浩博，張洪生，洪楊彬,，胡國棟

(1.上海海事大學(xué) 海洋科學(xué)與工程學(xué)院，上海 201306；2.浙江省嵊州市水利水電局，浙江 嵊州 312400；3.長江水利委員會長江口水文水資源勘測局，上海 200136)

長江口獨(dú)有的三級分汊、四口入海的地貌特征，使其具有復(fù)雜的風(fēng)浪特征。數(shù)值模擬是研究風(fēng)浪場的常用方法，但以往在利用數(shù)值模擬方法研究長江口水域的風(fēng)浪場時，要么只是數(shù)值模擬臺風(fēng)風(fēng)浪場[1-4]，要么是對包含該區(qū)域的西北太平洋相關(guān)海域的風(fēng)浪場進(jìn)行整體的數(shù)值模擬和分析[5-8]。查閱公開發(fā)表的文獻(xiàn)，未發(fā)現(xiàn)利用數(shù)值模擬通過計算連續(xù)多年的風(fēng)浪場，進(jìn)而統(tǒng)計分析長江口水域風(fēng)浪特征的相關(guān)研究。

數(shù)值計算方法在20 世紀(jì)中葉被引入到海洋學(xué)研究。早期的海浪模型因沒有考慮波浪內(nèi)部的復(fù)雜性，因此模擬精度不高。隨著波浪理論和計算機(jī)技術(shù)的發(fā)展，陸續(xù)提出了三代海浪數(shù)值模型[9]，其中WAVEWATCH III[10]（以下簡稱WW III）和SWAN[10-11]（Simulating Waves Nearshore）都屬于第三代海浪模型。其特點(diǎn)是所有源項進(jìn)行參數(shù)化之前，不對波譜加以限定，因此第三代海浪模型具有較強(qiáng)的適應(yīng)性。

影響長江口海域波浪的因素復(fù)雜，除了當(dāng)?shù)氐娘L(fēng)場之外，外海傳播過來的涌浪也是重要影響因素。為了提高風(fēng)浪場模擬的精度，需要選擇合適的計算范圍。本文通過敏感性分析來分別確定WW III 模型和SWAN 模型的計算范圍，其中確定SWAN 模型的計算范圍時以WW III 模型的輸出結(jié)果作為邊界條件。

本文以ECMWF 數(shù)據(jù)資料作為輸入風(fēng)場，以ETOPO1 地形數(shù)據(jù)和實測水深資料相結(jié)合作為計算水深，通過將WW III 和SWAN 模型進(jìn)行嵌套，計算長江口水域20 年（1996 年1 月至2015 年12 月）的風(fēng)浪場；并選取1 個代表性站點(diǎn)，根據(jù)該站點(diǎn)已有的觀測資料對計算結(jié)果加以修正；統(tǒng)計分析修正后的計算結(jié)果并利用風(fēng)場數(shù)據(jù)資料，進(jìn)而得出所選站點(diǎn)的風(fēng)場和風(fēng)浪場特征。

1 第三代海浪模型簡介

WW III 是目前比較成熟的海浪數(shù)值模型之一，具有穩(wěn)定性好、計算精度高等優(yōu)點(diǎn)。SWAN 模型來源于WW III 等第三代海浪數(shù)值模型，采用隱式格式離散控制方程。與WW III 相比，SWAN 模型更適合地形多變的河口海岸水域風(fēng)浪場的數(shù)值模擬，這是因為該模型能夠有效地模擬波浪在淺水域傳播變形過程中的特征，因此很多學(xué)者將兩個模型進(jìn)行嵌套計算[5,12-14]，其中文獻(xiàn)[12]對兩個模型進(jìn)行嵌套計算的精度進(jìn)行了詳細(xì)論證。因此，使用這兩種海浪模型進(jìn)行嵌套計算是合理有效的。

在球坐標(biāo)系下WW III 和SWAN 模型的控制方程均可表示為：

式中：t 為時間；σ 為相對頻率；θ 為波向；λ 和φ 分別代表經(jīng)度和緯度；cλ，cφ，cσ和cθ分別代表波作用量N 在λ，φ，σ 和θ 方向上的相應(yīng)傳播速度；為各物理過程中產(chǎn)生的源項，包括風(fēng)能輸入、波-波非線性相互作用和耗散項（白浪、水深變淺引起的破碎和底摩擦等），在WW III 和SWAN 模型中的形式略有不同，詳見文獻(xiàn)[10-11]。

2 模型建立與驗證

風(fēng)場數(shù)據(jù)的選擇對風(fēng)浪場模擬的精度有著重要影響。即使風(fēng)場數(shù)據(jù)為實測數(shù)據(jù)，觀測站點(diǎn)是否具有代表性也影響著最終結(jié)果的精度[12,15]。目前對海浪進(jìn)行數(shù)值模擬時，主要使用ECMWF（European Centre for Medium-Range Weather Forecasts，歐洲中期天氣預(yù)報中心）提供的風(fēng)場。ECMWF 風(fēng)場數(shù)據(jù)是當(dāng)前高精度氣象預(yù)報產(chǎn)品的代表[16]。文中采用分辨率為0.125°×0.125°的ECMWF 風(fēng)場作為外強(qiáng)迫。

地形資料采用ETOPO1（1 Minute Gridded Global Relief Data Collection）海底地形數(shù)據(jù)和實測數(shù)據(jù)。ETOPO1 數(shù)據(jù)因其分辨率高達(dá)1′×1′，且包含冰面和基巖的情況，可以較好反映海底地形，所以是目前進(jìn)行風(fēng)浪場數(shù)值模擬時常用的地形數(shù)據(jù)[17]?？紤]到長江口內(nèi)水深地形復(fù)雜，因此在SWAN 模型的計算范圍內(nèi)，部分采用實測水深數(shù)據(jù)以體現(xiàn)長江口深水航道等水工結(jié)構(gòu)物的影響，從而提高計算的精度。

長江口深水航道位于長江口南港北槽。北槽附近的牛皮礁（坐標(biāo)為122°15′06″E，31°08′12″N，記為點(diǎn)P）具有較強(qiáng)的代表性（見圖1），作為計算對象可以為今后航道整治和通航管理等提供參考。

圖1 研究范圍和代表性位置示意Fig.1 Sketch of study scope and representative location

在進(jìn)行風(fēng)浪場的數(shù)值模擬時，需要先確定模型的計算域和網(wǎng)格分辨率，因為這兩者會嚴(yán)重影響計算精度和計算效率。如果計算域太小，就不能充分考慮在其他海域生成和傳播而來的風(fēng)浪，從而使計算結(jié)果產(chǎn)生較大的誤差；如果計算域過大則又會增加不必要的計算時間。在確定計算網(wǎng)格分辨率時，同樣需要遵循在滿足計算精度的前提下，盡可能節(jié)省計算時間的原則。因為篇幅原因，文中只概述確定WW III 的計算范圍和嵌套計算兩個模型時網(wǎng)格類型的過程。確定SWAN 模型計算范圍的方法與WW III 類似，故不再贅述。模型參數(shù)設(shè)定為手冊的建議默認(rèn)值。

2.1 計算范圍的確定

首先基于往年的臺風(fēng)路徑和寒潮影響范圍，在長江口海域劃定一個大范圍區(qū)域（118°E～136°E，23°N～38°N）來確定WW III 的計算范圍，再通過算例來討論計算范圍邊界的具體位置。經(jīng)過敏感性分析最終確定WW III 的計算范圍為119°E～135°E，22°N～39°N；SWAN 的計算范圍為121°E～124°E，30°N～33°N。圖1（a）為計算范圍示意圖。

鑒于影響東邊界和南邊界所在區(qū)域海浪大小的主要因素為臺風(fēng)，文中以2013 年12 號臺風(fēng)“潭美”影響期間點(diǎn)P 的顯著波高為計算指標(biāo)，以1°為間距來確定東邊界和南邊界。而影響北邊界和西邊界所在區(qū)域海浪大小的主要因素為寒潮，所以以2010 年1 月發(fā)生寒潮期間點(diǎn)P 的顯著波高為計算指標(biāo)，以1°為間距來確定北邊界和西邊界。

在137°E 到132°E 之間、以1°為間距分別假定為東邊界，并保持SWAN 的計算范圍不變。WW III 的計算網(wǎng)格為6′×6′，SWAN 的計算網(wǎng)格為2′×2′。通過比較東邊界不同經(jīng)度時點(diǎn)P 的顯著波高，以確定東邊界的具體位置，計算結(jié)果如圖2 所示。由圖2 可見，東邊界為138°E～133°E 時波高的計算結(jié)果基本重合，而取132°E 時波高發(fā)生了較明顯變化。這說明在東邊界取為132°E 時，計算區(qū)域外的風(fēng)浪作用沒有被充分考慮而導(dǎo)致計算結(jié)果不準(zhǔn)確；而在東邊界取為137°E～133°E 時，對于外海影響的考慮較為充分。為了慎重起見，在分析了5 年（2011—2015 年）影響中國東南沿海的西北太平洋臺風(fēng)運(yùn)動軌跡后，取東邊界為135°E。

圖2 東邊界取不同經(jīng)度時P 點(diǎn)的顯著波高對比Fig.2 Comparisons of the calculated significant wave heights of point P with different east boundaries

用同樣的方法確定其余邊界，最終確定WW III 的計算范圍為119°E～135°E，22°N～39°N；SWAN 的計算范圍為121°E～124°E，30°N～33°N。

2.2 計算網(wǎng)格的確定

長江口灘槽相間的地形對海浪的傳播有很大的影響，因此計算網(wǎng)格對計算精度會有很大的影響。WW III 的計算網(wǎng)格分辨率設(shè)置為6′×6′；SWAN 的計算網(wǎng)格分別設(shè)置為結(jié)構(gòu)化網(wǎng)格（分辨率為2′×2′）和非結(jié)構(gòu)化網(wǎng)格（分辨率略低于0.5′×0.5′）。以2014 年8 號臺風(fēng)“浣熊”作用期間點(diǎn)P 的顯著波高為計算指標(biāo)，將WW III 和SWAN 的嵌套計算網(wǎng)格對應(yīng)組合，分別進(jìn)行計算，并與觀測資料加以對比（見圖3）。用于對比的觀測資料為牛皮礁水文站的現(xiàn)場觀測數(shù)據(jù)。

圖3 不同網(wǎng)格組合下P 點(diǎn)的顯著波高對比Fig.3 Comparisons of significant wave heights of point P for different combinations of grid mesh

由圖3 可見，采用非結(jié)構(gòu)網(wǎng)格計算的顯著波高與觀測值較為接近，86.7%的差值在0.5 m 以內(nèi)，采用結(jié)構(gòu)網(wǎng)格計算的結(jié)果則與觀測值有較大差距，最大差值達(dá)到了1.6 m 以上。因為計算的時間段內(nèi)包括臺風(fēng)過程，所以該算例具有一定參考價值，可以較好地反映計算結(jié)果的精度。

基于以上分析，為了在保證精度的前提下提高計算效率，采用分辨率為6′×6′的WW III 模型網(wǎng)格和SWAN 模型非結(jié)構(gòu)網(wǎng)格來嵌套計算長江口海域的風(fēng)浪場。SWAN 模型計算范圍內(nèi)共有51 618 個三角形網(wǎng)格，26 103 個計算節(jié)點(diǎn)，與文獻(xiàn)[18-20]所采用的網(wǎng)格分辨率近似。

2.3 計算結(jié)果的驗證及修正

由于篇幅所限，圖4 僅給出了2014 年7 月1 至28 日期間顯著波高的觀測值和計算值對比結(jié)果。由圖4可見，計算結(jié)果較好地模擬了波高的變化趨勢，兩者吻合程度較好。圖5 為2012 年到2015 年，共計4 年的顯著波高的散點(diǎn)圖。根據(jù)文獻(xiàn)[12]中提出的驗證方法對計算結(jié)果進(jìn)行修正。綜合4 年的結(jié)果來看，顯著波高絕對誤差為0.16 m，相對誤差為18.63%，相關(guān)系數(shù)為0.82。依據(jù)算出的皮爾遜相關(guān)系數(shù)對兩者進(jìn)行線性相關(guān)分析，當(dāng)相關(guān)系數(shù)為0.8～1.0 為極強(qiáng)相關(guān)。所以，顯著波高的計算值與觀測值兩者為極強(qiáng)相關(guān)，在以往的研究中也得到了類似的結(jié)論[12]。因此，可以通過線性擬合對計算結(jié)果進(jìn)行修正。

圖5 顯著波高的計算值和觀測值散點(diǎn)Fig.5 Plot of the observed and calculated values of significant wave heights

根據(jù)散點(diǎn)圖的擬合結(jié)果對顯著波高做了如下修正：

式中：Hm為顯著波高的修正值；Hc為顯著波高的原始計算值。

修正后的顯著波高見圖4。由圖4 可以看出，計算值修正后的結(jié)果要優(yōu)于修正前。其中顯著波高絕對誤差為0.05 m，相對誤差為8.24%。總體來看修正后波高的計算結(jié)果與觀測值的吻合程度較為理想，計算結(jié)果可以有效地反映出臺風(fēng)過程所引起的波高變化。

圖6 為波周期的觀測值和計算值對比，圖中波周期的計算值為平均波周期（TM01），為由波譜的零階矩和一階矩所求得的平均周期；波周期的觀測值為譜峰所對應(yīng)的周期。因為兩者并沒有直接對應(yīng)關(guān)系，所以本文沒有對周期進(jìn)行修正。由圖6 可以看到，兩者的變化趨勢一致，但是計算值明顯小于觀測值。其原因一方面與海浪模型本身有關(guān)，在以往的研究中也得到了波周期偏小的結(jié)論[12]；另一方面也符合在文獻(xiàn)[21]中對兩種不同類型波周期大小的分析，即通常譜峰周期會大于平均波周期。

圖6 波周期的觀測值和計算值對比Fig.6 Comparisons of the observed and calculated values of wave periods

3 特征點(diǎn)的風(fēng)浪特性統(tǒng)計與分析

3.1 風(fēng)場特性統(tǒng)計與分析

利用ECMWF 風(fēng)場數(shù)據(jù)，對點(diǎn)P 總計20 年的風(fēng)場數(shù)據(jù)按16 個方向進(jìn)行統(tǒng)計分析。表1 所示為點(diǎn)P 的多年各月平均風(fēng)速；圖7 為點(diǎn)P 的多年風(fēng)玫瑰圖。點(diǎn)P 多年平均風(fēng)速為6.2 m/s。由表1 和圖7 可見，全年常風(fēng)向及大風(fēng)向均為NW～NE，秋冬兩季6 級以上大風(fēng)的發(fā)生頻率較高，且平均風(fēng)速也大于春夏兩季。

表1 點(diǎn)P 多年平均風(fēng)速統(tǒng)計Tab.1 Average wind speeds for point P

圖7 點(diǎn)P 多年風(fēng)玫瑰圖Fig.7 Wind rose plot of point P

由圖7 可見，點(diǎn)P 常風(fēng)向為N 向，WSW 向和W 向風(fēng)頻率最低。風(fēng)向季節(jié)特征明顯，冬季處于東北向季風(fēng)時期，偏北風(fēng)盛行，強(qiáng)風(fēng)向多為NW～N；夏季受副熱帶高壓影響，偏南風(fēng)較多，風(fēng)向以SE～S 為主；春秋兩季為季風(fēng)轉(zhuǎn)換時期，春季風(fēng)向以SE～S 向為主，N～NNE 向為次；秋季時冬季季風(fēng)已開始增強(qiáng)，因此秋季風(fēng)向以N～ENE 向為主。

3.2 波高和波向的聯(lián)合分布

分別統(tǒng)計全年、各季度和各月修正后的點(diǎn)P 多年平均H1/10波高和波向特征。表2 為點(diǎn)P 的波高與波向聯(lián)合分布表（全年）。由于篇幅原因，其他聯(lián)合分布表未列出。圖8 為點(diǎn)P 的多年各級波高玫瑰圖。

表2 點(diǎn)P 多年顯著波高與波向聯(lián)合分布(全年)Tab.2 Joint probability distribution of significant wave heights and directions for point P 單位：%

圖8 點(diǎn)P 多年各級波高玫瑰圖Fig.8 Wave rose plot of point P

由表2 和圖8 可見點(diǎn)P 的H1/10波高主要分布在1.5 m 以下的波級，且集中分布于0.7～1.2 m，常浪向為NE 向，頻率為16.3%。波高季節(jié)性變化明顯，秋冬兩季大浪出現(xiàn)的頻率較高，這是由于秋冬兩季處于冬季季風(fēng)時期，此時的風(fēng)速為全年最大。全年浪向主要分布于NE～S 連線的右側(cè)。春季浪向以NNE～SSE 為主，夏季以ENE～S 為主，秋季以NNE～ENE 為主，冬季集中于NNW～NE。與文獻(xiàn)[6-7]中所得到的波高與波向的季節(jié)特征相一致，主要原因均是受風(fēng)場的影響。對比圖7 和8 可見，浪向的季節(jié)特征與風(fēng)向并不完全相同。例如，雖然夏冬兩季都有相當(dāng)比例的偏西風(fēng)，但同向來浪頻率卻明顯較低；雖然春季的偏東向風(fēng)不多，但卻有較多的同向風(fēng)浪。這是因為風(fēng)浪除了受風(fēng)場影響外，外海涌浪和地形也是重要的影響因素。

3.3 波高和周期的聯(lián)合分布

表3 點(diǎn)P 顯著波高與波周期平均值的統(tǒng)計Tab.3 The table of averaged significant wave heights and mean wave periods for point P

表4 點(diǎn)P 多年顯著波高與平均波周期聯(lián)合分布(全年)Tab.4 Joint probability distribution of significant wave heights and mean wave periods for point P 單位：%

3.4 不同重現(xiàn)期的風(fēng)浪計算

根據(jù)計算的20 年的波高，假設(shè)點(diǎn)P 風(fēng)浪極值滿足耿貝爾極值I 型分布，分別計算了點(diǎn)P 不同重現(xiàn)期的H1/10極值波高，并對風(fēng)浪極值進(jìn)行了K-S 檢驗。檢驗結(jié)果表明，當(dāng)置信度為0.95 時，接受點(diǎn)P 風(fēng)浪極值滿足耿貝爾極值I 型分布的假設(shè)。圖9 為點(diǎn)P 不同重現(xiàn)期的H1/10極值波高圖，其中點(diǎn)P 的50 年一遇極值波高為7.1 m。

因為極值波高的最大來浪向為E 向，所以計算了鄰近海域E 向的極值波高分布，圖10 為點(diǎn)P 鄰近海域50 年一遇E 向極值波高分布圖。從圖10 可以看出，極值波高由東向西遞減，這是由于東邊是較開闊的外海，其水深和風(fēng)速都較大，而無論是水深還是臺風(fēng)吹程都會影響波高的分布規(guī)律。

圖9 點(diǎn)P 不同重現(xiàn)期的H1/10 極值波高Fig.9 Plot of extreme significant wave heights of point P for different return periods

圖10 點(diǎn)P 鄰近海域50 年一遇E 向極值波高分布（單位：m）Fig.10 Extreme significant wave heights of 50 years return distribution (unit: m)

4 結(jié) 語

以ECMWF 數(shù)據(jù)資料作為輸入風(fēng)場，以ETOPO1地形數(shù)據(jù)和實測水深資料相結(jié)合作為計算水深，通過將WW III 和SWAN 模型進(jìn)行嵌套，模擬計算了長江口海域20 年的風(fēng)浪，并根據(jù)已有的觀測數(shù)據(jù)對計算結(jié)果加以修正。根據(jù)風(fēng)場數(shù)據(jù)資料和修正后的波高，統(tǒng)計分析了長江口水域具有代表性的點(diǎn)位—牛皮礁的風(fēng)場和風(fēng)浪場特性。

（1）牛皮礁處常風(fēng)向為N 向，WSW～W 向來風(fēng)最少。風(fēng)向和風(fēng)速隨季節(jié)變化明顯，冬季以北向季風(fēng)為主；夏季以偏南風(fēng)為主；春秋兩季為季風(fēng)轉(zhuǎn)換時期，春季東南風(fēng)稍多；秋季東北風(fēng)稍多。秋冬兩季受冬季季風(fēng)影響風(fēng)速大于春夏兩季。

（2）受風(fēng)場影響，牛皮礁處波高與波周期也有秋冬大，春夏小的特征。牛皮礁處常浪向為NE 向，出現(xiàn)頻率為16.3%，W 和WSW 向浪出現(xiàn)的頻率最小，為0.4%；季節(jié)特征明顯，不過在地形和外海涌浪的影響下浪向與風(fēng)向略有不同，浪向春季主要分布在NNE～SSE，夏季主要分布在ENE～S，秋季主要分布在NNE～ENE，冬季集中分布在NNW～NE。

（3）計算了牛皮礁處在不同重現(xiàn)期的H1/10極值波高和鄰近海域50 年一遇E 向極值波高分布，結(jié)果顯示50 年一遇的極值波高為7.1 m，附近海域極值波高呈現(xiàn)自西向東逐漸增大的特征。

本文的數(shù)值計算結(jié)果較好地體現(xiàn)了風(fēng)場、地形和外海涌浪的作用，而且依據(jù)實測資料進(jìn)行了修正。因此計算結(jié)果是合理的，統(tǒng)計分析和討論所得出的結(jié)論是可信的。但實際的風(fēng)浪生成和傳播會受到多種因素的共同影響，因此風(fēng)浪的數(shù)值模擬還可以進(jìn)一步發(fā)展和完善。比如第三代海浪數(shù)值模型在周期的計算方面還存在一些不足，模型本身需要改進(jìn)；此外，水位、潮流和風(fēng)暴潮等都是影響風(fēng)浪的重要因素，且長江口的地形因為工程建設(shè)原因在不斷發(fā)生變化，在日后的研究中需要盡可能地將這些因素一并考慮，并采用分辨率更高的風(fēng)場和地形資料，以得到更加符合實際的結(jié)果。