基于多源數(shù)據(jù)的臺風風暴潮概率預報研究：數(shù)值預報系統(tǒng)

郭文云，安佰超，裘誠，李鋮，李丕學，葛建忠，丁平興

（1.上海海事大學海洋科學與工程學院，上海201300；2.上海市海洋監(jiān)測預報中心，上海200062；3.華東師范大學河口海岸學國家重點實驗室，上海200062）

1 引言

臺風引起的風暴潮災害每年都給我國造成大量經(jīng)濟損失[1-2]。開展風暴潮預報是防止（減少）風暴潮災害并且保障人民生命財產(chǎn)安全的有效手段。隨著SPLASH（Special Program to List Amplitudes of Surges from Hurricanes）、SLOSH（Sea，Lake，and Overland Surges from Hurricanes model）和ADCIRC（The ADvanced CIRCulation model）等模型的應用，數(shù)值模型已成為風暴潮預報的重要方法[3-9]。風暴增水強度對臺風路徑和強度等參數(shù)很敏感，但臺風路徑和強度預報存在很大的不確定性。據(jù)統(tǒng)計，2016年中央氣象臺臺風24 h、48 h和72 h預報的平均路徑誤差分別為76.2 km、147.3 km和244.7 km，平均風速誤差分別為5.4 m/s、7.3 m/s和7.6 m/s[10]。這些臺風預報誤差給風暴潮數(shù)值預報工作帶來挑戰(zhàn)。充分考慮臺風預報中的誤差，進行多個臺風可能路徑及可能強度的風暴潮集合數(shù)值預報，從而提供更多更全面的預報信息，可以有效彌補這一不足[11-15]。對集合預報中樣本的發(fā)生概率進行仔細估計，稱為概率預報。

目前，多個國家的風暴潮業(yè)務化預報系統(tǒng)已經(jīng)開始采用集合（概率）預報技術。美國國家颶風中心（National Hurricane Center）使用的P_surge風暴潮概率預報模型考慮了臺風路徑、臺風中心附近最大風速和臺風大小的不確定性，根據(jù)歷史預報誤差統(tǒng)計生成27個臺風樣本來進行風暴潮概率預報[16]。澳大利亞氣象局（Bureau of Meteorology Australia）則基于特定方法生成1 000條可能的臺風路徑，然后從中隨機選取200條來進行風暴潮集合預報[17]。

我國在風暴潮集合預報研究上也取得不少進展。國家海洋環(huán)境預報中心已經(jīng)開發(fā)了一套風暴潮概率預報模式，并實現(xiàn)了業(yè)務化運行。該模式主要考慮了臺風路徑的不確定性，計算最大可能情況下的風暴增減水極值。利用該集合預報模式已對0814號強臺風“黑格比”進行數(shù)值預報實驗[18]。王培濤等[12]基于預報路徑衍生偏快、偏慢、偏左和偏右4條路徑對福建沿海進行了風暴潮集合預報。陳永平等[19]在5條可能路徑（包括預報路徑）的基礎上進一步考慮了3個可能的最大風速，使用15個臺風樣本進行風暴潮集合計算。以上結果都顯示，采用集合預報技術可以有效降低風暴潮預報誤差。目前我國的研究大多還停留在集合預報上，只是通過構造不同臺風樣本來進行多個情景的預報，并通過分析預報增水平均值來給出預報增水幅度，在統(tǒng)計學上缺乏實際概率意義，并未達到真正的概率預報。事實上，不同臺風樣本的發(fā)生概率應該有所不同。

郭文云等[1]建立了一個臺風集合構造方案。該方案可構造27個臺風樣本（9條可能路徑和3個可能最大風速），并基于歷史預報誤差的統(tǒng)計分布來確定不同臺風樣本的發(fā)生概率。該集合構建方案可以較為全面地覆蓋可能的情景。本文基于這一臺風集合構造方案，采用中尺度天氣預報（Weather Research and Forecasting，WRF）模式和有限體積海岸海洋模式（Finite Volume Coastal Ocean Model，F(xiàn)VCOM）建立長江口臺風風暴潮集合預報系統(tǒng)，并以1909號臺風“利奇馬”為例，詳細說明了概率預報提供的豐富的風暴潮概率預報信息。這些信息可以為風暴潮防災減災工作提供詳細的數(shù)據(jù)參考。

2 長江口風暴潮數(shù)值模型系統(tǒng)

臺風是一種局地的強非線性氣象系統(tǒng)，其最大風速半徑一般為幾十公里，對風暴潮進行數(shù)值模擬需要精細地刻畫臺風的風場。臺風中心具有快速移動性，因此需要較大的模型區(qū)域來準確模擬風暴增水過程。本文采用高分辨率臺風風場疊加大尺度背景風場的方法來提供風暴潮模型所需的表面風應力輸入。其中臺風風場使用經(jīng)典的藤田公式來描述，大尺度背景風場由三重嵌套的WRF中尺度氣象模型提供，而風暴潮模型則選用國際先進的FVCOM三維模型進行兩重嵌套。該模型采用無結構網(wǎng)格，可以精確地刻畫長江口復雜的岸線和地形條件。

2.1 WRF氣象模型

WRF模式具有計算高效、可移植性強、方便且易于維護等優(yōu)點，適合于大氣動力學研究和業(yè)務化預報等工作。為提高WRF的計算效率，同時保證計算范圍足夠大，以涵蓋臺風風暴潮的生成、發(fā)展和消亡的全過程，因此采用三重嵌套和雙向數(shù)據(jù)交換的方法。模式的計算區(qū)域如圖1所示。WRF大區(qū)覆蓋從15°～51°N，109°～149°E的區(qū)域；中區(qū)覆蓋從23°～41°N，117°～129°E的區(qū)域；而小區(qū)覆蓋長江口和杭州灣沿海地區(qū)及江浙滬部分地區(qū)。各區(qū)域的分辨率分別為81 km、27 km和9 km。

圖1 長江口風暴潮預報系統(tǒng)中的三重嵌套WRF氣象模式區(qū)域

WRF模型在進行后報時采用美國國家環(huán)境預報中心（National Centers for Environmental Prediction，NCEP）提供的經(jīng)過同化后的FNL（Final Operational Global Analysis）數(shù)據(jù)，空間分辨率為1°×1°，時間間隔為6 h；預報時則采用NCEP提供的全球預報產(chǎn)品GFS（Global Forecast System），空間分辨率為0.5°×0.5°，時間間隔為3 h。

2.2 臺風模型

本文采用對稱臺風氣壓場模型。模型需要用到的臺風參數(shù)包括臺風位置、臺風中心氣壓P0、臺風中心附近最大風速Vmax和最大風速半徑Rmax。其中最大風速半徑無法在預報臺風中給出，因此使用經(jīng)驗公式[20]進行估計：

臺風氣壓場使用藤田公式[21]來估計：

式中：P0為臺風中心氣壓；P∞為臺風外圍氣壓；r是計算點至氣壓中心的距離；R0是表征臺風系統(tǒng)特征的參數(shù)，可根據(jù)最大風速半徑調整。

氣旋的風速分布根據(jù)梯度風關系[22-23]參照氣壓分布給出：

式中：f=2Ωsinφ為柯氏參數(shù)；Ω為地球自轉角速度；φ為緯度；ρa為空氣密度。臺風移行產(chǎn)生的風場采用宮崎正衛(wèi)公式[24-25]：

式中：Vx和Vy分別是臺風中心移動速度的正東分量和正北分量。

2.3 風場融合

合成后的臺風風場可以表示為：

W→=c1W1[-sin(?+β),cos(?+β)]+c2W→2（5）

式中：?是計算點與臺風中心連線與正東方向的夾角；β是梯度風與海面風的夾角；c1和c2是訂正系數(shù)。

由對稱模型風場計算所得的風場與背景風場之間的合成，根據(jù)計算點與臺風中心的距離r與判別距離R1和R2的關系進行插值合成：

式中：a=(r-R1)/(R2-R1)。

表面風應力的計算采用二次率參數(shù)化形式：

式中：ρa為空氣密度，取1.23 kg/m3；W→為風矢量；CD為拖曳系數(shù)，隨風速大小而變化：

2.4 風暴潮模型

長江口風暴潮模型采用葛建忠等[26-27]建立的東中國海-長江口FVCOM嵌套模型（見圖2）。東中國海大區(qū)模型覆蓋我國渤海、黃海、東海及其鄰近海域。長江口精細化模型覆蓋了長江各個入海口以及杭州灣和舟山群島。長江口上界取至潮流界江陰以上200 km左右的大通；錢塘江上界取至老鹽倉附近；外海開邊界東至124.5°E，北至34.5°N，南至28.5°N。網(wǎng)格分辨率在外海陸架海域為0.5～5 km，在長江口則精細到200 m左右。模型開邊界基于葛建忠等[26]的數(shù)據(jù)調整得到，在8個主要分潮（M2、S2、K1、O1、N2、K2、P1和Q1）的基礎上，加入兩個長周期分潮（Mf和Mm）和3個淺水分潮（M4、MS4和MN4）。所加入的分潮都是在OSU TOPEX/POSEIDON全球海洋潮汐模式[28]的基礎上進行適當調整后的結果。

2.5 模型驗證

2018年7月18日20時（北京時，下同），臺風“安比”（1810）在西北太平洋洋面上生成（見圖3a）。2018年7月22日12時30分前后，臺風“安比”在上海市崇明島沿海登陸，登陸時中心附近最大風力達10級（28 m/s），中心最低氣壓為982 hPa，成為1949年以來第三個直接登陸上海的臺風。2018年7月24日夜間，臺風“安比”（熱帶風暴級）減弱變性為溫帶氣旋，之后強度持續(xù)減弱，中央氣象臺于25日02時對其停止編號。

圖2 長江口風暴潮模型的東中國海-長江口嵌套網(wǎng)格

2019年8月4日，臺風“利奇馬”（1909）生成（見圖3b），并于2019年8月10日01時45分左右在浙江省溫嶺市城南鎮(zhèn)沿海登陸，登陸時中心附近最大風力16級（52 m/s），屬超強臺風，是1949年以來登陸我國的第五個超強臺風。隨后臺風“利奇馬”縱穿浙江和江蘇兩省并移入黃海海面，又于8月11日20時50分許在山東省青島市黃島區(qū)沿海再次登陸，此后其移入渤海海面并不斷減弱，最終于8月13日14時被中央氣象臺停止編號。截至2019年8月14日10時，臺風“利奇馬”共造成我國57人死亡，直接經(jīng)濟損失537.2億元。金山嘴潮位站的水位觀測數(shù)據(jù)顯示，臺風“利奇馬”造成該站最大增水超過1.0 m。

采用本文構建的風暴潮模型對臺風“安比”和“利奇馬”進行后報檢驗，所得南門、堡鎮(zhèn)、蘆潮港和金山嘴4站的增水過程對比如圖4所示。本文構建的風暴潮模型不僅可以準確地刻畫臺風“安比”和“利奇馬”引起的增水過程，對臺風過后的減水過程也能準確模擬。南門、堡鎮(zhèn)和蘆潮港的平均增水誤差都小于15 cm。根據(jù)觀測數(shù)據(jù)，金山嘴站在兩個臺風過后都存在較強余振（振幅超過40 cm），本文的模擬結果雖然也能體現(xiàn)余振效果，但對余振的刻畫還不夠精細，使得該站的平均誤差超過20 cm。這可能是由于臺風登陸后風場發(fā)生變形，而本文所采用的對稱臺風模型難以刻畫這種變形。模型模擬的最大增水時間與觀測基本一致，最大增水幅度比觀測稍強，從預報安全性角度上是可以接受的。此外，臺風“利奇馬”引起的長江口增水都帶有雙峰特征，本文的模型也能反映出類似的特征。

3 預報實驗

本文的主要目的是建立長江口概率預報系統(tǒng)。相對于單一預報，概率預報系統(tǒng)能提供更多的風暴增水預報結果。如何從這些結果中挖掘出所需的風暴增水信息，也是風暴潮概率預報的一項重要工作。

仍然以臺風“利奇馬”為例，選取蘆潮港站最大增水發(fā)生前24 h左右（2019年8月9日14時）為預報時刻進行虛擬預報實驗。采用郭文云等[1]的臺風集合構造方案，先根據(jù)中國中央臺、中國香港臺、中國臺灣臺、美國臺、日本臺和韓國臺的預報路徑生成一條誤差更小的分析路徑（見圖3中紅色線），然后基于分析路徑及其統(tǒng)計誤差生成兩個概率圓上的8條衍生路徑（見圖3），路徑集合覆蓋了實測路徑的運移范圍。類似方法可以得到預報臺風最大風速集合（見表1）。其中48 h和72 h預報時效最大風速都覆蓋到了實測最大風速，但24 h預報最大風速沒有覆蓋到實測最大風速，實測風速弱于概率偏小的預報最大風速。

圖3 兩個臺風的實測路徑和臺風路徑樣本構建（注：臺風強度標注只用于實測路徑）

以上路徑集合和最大風速集合相互交叉形成27個臺風樣本，每個臺風樣本的發(fā)生概率都根據(jù)分析臺風的誤差分布確定。對以上27個臺風樣本分別進行風暴潮預報，可得到27個風暴潮預報結果。

圖4 兩個臺風過程中4站的實測和后報增水過程對比

3.1 最大和最小可能增水分布

將以上27個預報中的水位數(shù)據(jù)輸出，那么對于任意時刻任意空間點，都存在27個可能的風暴水位。同時進行只考慮天文潮汐而不考慮臺風過程的模擬，相減即可得到27個不同的增水場。其中最大增水值屬最壞情況，最小增水值屬最好情況。從預報24 h后即2019年8月10日14時的最大和最小可能增水分布圖上可以看到（見圖5），此時長江口4條分支和杭州灣內的最大可能增水都超過2.0 m，但最小可能增水還不到0.5 m；長江口門附近甚至還有可能出現(xiàn)減水情況。最大和最小可能增水之間存在巨大差異，充分說明臺風預報中的路徑和風速等參數(shù)的誤差會使風暴潮預報增水產(chǎn)生巨大偏差，也體現(xiàn)了建立風暴潮概率預報系統(tǒng)的重要性。若采用確定性的風暴潮單一預報系統(tǒng)，以上信息都無法提供。

3.2 平均增水及其標準偏差分布

結合以上27個臺風樣本的風暴潮預報結果及其發(fā)生概率，可以計算不同時刻不同空間點的平均增水值及其標準差：

式中：ηi為第i個數(shù)值預報實驗的增水值；pi為該數(shù)值實驗所用臺風樣本的發(fā)生概率。

圖6為24 h預報時效（即2019年8月10日14時）的平均增水分布及其標準差。圖中可以看出，此時東中國海的增水并不明顯，一般在0.5 m以下；沿長江口向陸以及沿杭州灣向內，平均增水值不斷增大，長江口4條分支和杭州灣內的平均增水都超過1.5 m，但實測增水為1.0 m左右，這與臺風預報最大風速較實際最大風速顯著偏大有關。即使是概率偏弱（20%）的臺風最大風速都比實測最大風速更大（見表1）。該時刻的預報增水標準偏差也呈現(xiàn)外海較小近岸較大的空間分布特征，在長江口內和杭州灣內，大部分區(qū)域的標準偏差大于0.5 m，表明該時刻預報增水的不確定性較大。這些信息對于風暴潮的防災減災工作都具有重要的參考價值。

圖5 臺風“利奇馬”24 h預報時效的可能增水（預報時間為2019年8月9日14時；單位：m）

圖6 臺風“利奇馬”24 h預報時效結果（預報時間為2019年8月9日14時，單位：m）

為方便對比，本文同時也將分析路徑的預報增水結果（代表確定性單一預報結果）呈現(xiàn)于圖6c。其總體增水場分布與預報平均增水分布相似，但在長江口內和杭州灣內的增水值略微更大。只根據(jù)這單一預報結果，我們無從得知該時刻增水的不確定性。由于預報臺風路徑存在明顯誤差，勢必導致預報風暴潮增水的不確定性。相對于單一預報系統(tǒng)，概率預報系統(tǒng)則可以提供這些信息，充分體現(xiàn)了建立風暴潮概率預報系統(tǒng)的必要性。

3.3 概率增水值和增水概率值

風暴潮概率預報系統(tǒng)可以得到某個站點在未來某個時刻的多個增水值及其發(fā)生概率。將這些增水值從大到小排列并計算其累積可得到該站點在某時刻的增水概率圖。圖7為24 h預報時效的蘆潮港站臺風增水累積概率圖，在所有可能情況中，該時刻蘆潮港站都處于增水狀態(tài)。對于某個固定的概率，其增水值可以達到多大，此為概率增水值。蘆潮港站在該時刻的20%概率增水值為1.65 m（見圖7a）。

風暴潮可以根據(jù)其最大增水值來劃分等級，一般以1.0 m、1.5 m、2.0 m和2.5 m為界將風暴潮分為一般、中等、較大、大和特大5個等級[29]。我們同樣希望得到對于某個固定的增水值，其發(fā)生概率如何，此為增水概率值。根據(jù)圖7b可以確定蘆潮港站（高于）1.0 m增水概率為84.9%。

依據(jù)以上對概率增水值和增水概率值的定義，我們可以根據(jù)概率預報結果得到未來某個時刻的概率增水分布及增水概率分布。圖8a為2019年8月10日14時預報的20%概率增水分布。圖中可以看到外海的增水值與平均增水沒有太大差別，表明臺風預報誤差對外海的風暴增水影響較小，而在長江口內及杭州灣內，20%概率增水超過2.0 m。

圖8b為1.0 m增水概率分布，即預報該時刻增水高于1.0 m的概率。在外海，這一概率都為0，表明外海增水不可能超過1.0 m；一旦進入河道或杭州灣內，由于地形縮窄，潮能匯聚，1.0 m增水概率迅速增加，直至達到100%。

3.4 站點時間序列

將以上27個模擬中的重要站點數(shù)據(jù)輸出，可以得到27條不同的風暴過程水位線（見圖9a）。結合不同臺風樣本的發(fā)生概率，可以計算站點在指定發(fā)生概率范圍內的增水范圍。圖9b呈現(xiàn)了蘆潮港站預報的100%和70%概率范圍的增水，其中垂向虛線指示預報時刻。圖中可以看到隨著預報時效的延伸，預報的增水100%范圍先隨著臺風預報誤差的增長逐漸加大；然后臺風逐漸遠離預報區(qū)域，即使臺風預報誤差更大，但由于臺風對預報區(qū)域的水位影響逐漸減小，增水100%范圍也逐漸減小。2019年8月11日00時，由于臺風運動到離長江口較近的位置，臺風預報誤差也較為顯著，此時的增水100%范圍達到最大。蘆潮港站最大可能增水幾乎可以達到2 m，最小可能增水卻只有-1.80 m，顯示出顯著的風暴潮預報的不確定性。預報增水70%范圍同樣也呈現(xiàn)先增大后減小的特征，但其最大范圍比100%范圍顯著減小，其發(fā)生的時間也與100%范圍的略有差異。

圖7 2019年8月9日14時蘆潮港站的臺風增水累積概率圖

圖8 臺風“利奇馬”24 h預報時效結果（預報時間為2019年8月9日14時)

圖9 臺風“利奇馬”影響下不同站點的預報結果（預報時間為圖中虛線所示的2019年8月9日14時）

圖9顯示本次預報的最大增水值比實測的最大增水顯著偏大，這是由于臺風預報的最大風速顯著偏大所致；特別是在24 h預報時效，即使是概率偏小的最大風速，依然大于實測的最大風速。

表2 4個潮位站概率預報結果統(tǒng)計檢驗（概率預報為平均值加1個標準差）

3.5 集合預報精度檢驗

采用圖4中臺風“利奇馬”期間南門、堡鎮(zhèn)、蘆潮港和金山嘴4個潮位站的逐時觀測數(shù)據(jù)對本文建立的概率預報系統(tǒng)進行檢驗。最大增水及其發(fā)生時刻是風暴潮預報的重要參數(shù)。本文對其誤差進行統(tǒng)計，同時將分析路徑作為確定性單一預報的代表，所得結果見表2。

臺風“利奇馬”活動期間，南門、堡鎮(zhèn)、蘆潮港和金山嘴4站的實測最大增水分別為0.84 m、0.66 m、0.84 m和1.28 m。9日14時對臺風“利奇馬”近中心最大風速的24 h預報風速顯著偏大，對于24 h預報，即使是概率偏弱的風速（30.4 m/s）也比實測風速（28 m/s）更大（見表1），因此，預報實驗所得最大增水的估計也較實測最大增水顯著偏大。27個集合樣本所得到的4站最大增水分別為0.80～3.33 m、0.73～2.88 m、1.24～2.47 m和1.73～3.11 m。集合平均增水過程的最大值（見圖9紅線）分別為1.87 m、1.50 m、1.90 m和2.56 m，都較實測最大增水偏大約1.0 m。相對于確定性單一預報，本文建立的概率預報系統(tǒng)對最大增水的估計都更接近實測最大增水，表明本文建立的概率預報系統(tǒng)可以通過集合消減來提高風暴潮的預報精度。不僅如此，概率預報系統(tǒng)還能夠提供更多的可用信息，如表2中所示的各站最大增水的標準誤差。誤差越大，表明臺風路徑和強度等參數(shù)中的預報誤差對該站的最大增水影響越顯著。對于這些海域，實際業(yè)務中需要特別注意臺風實際路徑變化的影響。

臺風“利奇馬”期間，南門和堡鎮(zhèn)站的實測增水在8月10日17時達到最大，蘆潮港和金山嘴站則在14時達到最大。概率預報和確定性預報對最大增水時刻（約24 h預報時效）的估計誤差為2 h左右（見表2）。概率預報和確定性預報的最大增水時刻誤差沒有顯著差異，但采用概率預報技術可以得到更多的信息。例如：雖然預報南門和堡鎮(zhèn)站的最大增水時刻與實測只相差約1.5 h，但其標準差分別為5.14 h和6.38 h，比蘆潮港和金山嘴的標準誤差明顯更大。因此，如果臺風的移動路徑、移動速度和臺風強度等參數(shù)出現(xiàn)偏差，這兩站的最大增水時刻有可能變化更大。

4 結論與展望

本文采用郭文云等[1]的臺風集合構造方案，基于WRF氣象模型和FVCOM三維海洋模型，建立了一套適用于長江口的風暴潮概率預報系統(tǒng)。為兼顧臺風風速空間變化強烈和移動范圍大的特點，WRF模型采用三重嵌套模式，F(xiàn)VCOM海洋模型采用兩重嵌套，在保證精度的前提下提高計算效率。系統(tǒng)采用9條可能路徑和3個可能最大風速進行27個風暴過程模擬。每個模擬風暴過程的發(fā)生概率都根據(jù)歷史預報誤差分布進行確定。

風暴潮概率預報較單一預報的優(yōu)點是提供了多個可能的預報增水值及各種豐富的預報信息。本文建立了一套詳細分析風暴潮概率預報結果的產(chǎn)品，如最大和最小可能增水值、平均增水值及其標準偏差、概率增水值及增水概率值等。基于以上分析，我們可以從風暴潮概率預報結果中得到更全面的信息，為防災減災工作提供更多數(shù)據(jù)參考。當然，本文所建立的長江口概率預報方案還存在一些可改進之處。其一，預報方案只考慮了預報臺風的路徑及最大風速誤差，沒有考慮預報臺風中心氣壓的誤差，而臺風中心氣壓直接影響到臺風的最大風速半徑（見式（1）），是決定臺風大小的重要參數(shù)，因此，臺風中心氣壓的預報誤差也應該在風暴潮概率預報系統(tǒng)中加以考慮；其二，本文所建立的概率預報系統(tǒng)采用了9條可能路徑，使得預報系統(tǒng)運行過程中需要用到較多的計算資源，需要進一步做不同臺風集合方案的對比，在計算資源與預報精度上取得合理的平衡；其三，本文所采用的概率預報方案是基于歷史統(tǒng)計得來的，實際預報過程中，不同臺風的預報誤差也有所不同（如對于轉向臺風及走向異常的臺風，其預報誤差往往更大）。所以，加強與臺風預報部門的合作，爭取采用實時的臺風預報誤差數(shù)據(jù)是改進風暴潮概率預報系統(tǒng)的根本途徑。