臺風浪集合預報方法研究

劉安琪，李鋮，郭文云，葛建忠*

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

0 引言

臺風期間，受臺風影響的海域的環流、水位、波浪等要素通常會發生劇烈變化，其中水位的異常變化會導致風暴潮現象的發生，而劇烈的海表風場會激發產生巨大的波浪，即臺風浪現象[1]。與風暴潮相比，臺風浪作用周期更短，但產生的瞬時破壞作用更大，會破壞近岸工程，并威脅人們的生命以及財產安全。2019年，我國近海共發生有效波高4.0 m（含）以上的災害性海浪過程39次，因災直接經濟損失達0.34 億元，死亡（含失蹤）22 人。臺風浪在向近岸傳播的過程中，會與復雜海底地形以及不規則岸線發生非線性相互作用[2]，使得波浪預測難度增大，另外近岸波浪觀測資料的缺乏也限制了區域波浪預報產品預報精度的提升。

對臺風浪的預報，國內外很多學者從不同的方向做了很多努力，取得了豐碩的成果。CHAWLA等[3]利用了美國國家環境預報中心（National Centers for Environmental Prediction，NCEP）的CFSR（Climate Forecast System Reanalysis）再分析數據，使用第三代風浪模型（WAVEWATCH III）進行了31 a 的后向預報，分析發現，在南半球，CFSR對風速的過度預測對波浪預測的準確性有重要影響。ALVES 等[4]基于NCEP 和美國海軍數字氣象和海洋學研究中心（Fleet Numerical Meteorology and Oceanography Center，FNMOC）的數值模式對有效波高等要素進行綜合概率預報，結果表明多模式聯合預報方式的準確性比單一數值模式有了很大的提升。FAN 等[5]使用多個全球大氣模式數據作為波浪模型輸入得到的有效波高，與空間上不同近岸浮標的觀測值的吻合程度各不相同。PAN 等[6]提出了一種優化模型-集合波高的算法，計算了不同風模式輸入波浪模型產生的波高，該算法在空間上給不同風模式賦不同的權重，調整后的波高在近岸范圍內與浮標觀測值的吻合較好。ZIEGER 等[7]使用歐洲中期天氣預報中心集合預報模式（The European Centre for Medium-Range Weather Forecasts-Ensemble Prediction System，ECMWF-EPS）中經過矯正的大氣風場生成的集合風場進行波浪模擬，并對比了近岸區域的浮標觀測數據，部分風場的結果和觀測的有效波高峰值吻合較好。

綜上可見，近岸波浪預報結果的準確性在很大程度上依賴于風場的準確性。近岸波浪由于各物理過程的非線性作用導致不確定性很強，很難找到一種傳統模式或者參數化方案在整個預報時段內都表現較好。在實際預報中，既要關注預報結果的準確性，也要關注預報結果的概率區間預測。本文旨在提出一種新的臺風路徑集合方法開展波浪預報實驗，并給出臺風浪的極值范圍作為預報參考。

集合風場目前主要有以下幾種生成方法：一種是使用單一業務預報中心的數據作為數據源（如應用較多的增長繁殖法[8]）以及特定業務中心集合預報系統中的集合成員[7]，這種方法的優勢在于數據一致性較好，不足之處是會引入單一模式的固有誤差；另一種是使用多業務預報中心數據成員，基于“共識”的方式挑選出其中較好的成員[9]；第三種是基于不同業務中心歷史預報誤差分配權重，使用“概率圓”的方法進行集合概率預報[10]，這種方法能避免單一模式的偶然誤差，不足在于各預報源的“好壞”一般由其歷史預報誤差確定，無法在預報過程中實時變動；還有一種是使用研究區域內歷史臺風運動數據作為數據源，并基于臺風的歷史最佳路徑集數據確定所有狀態的概率轉移情況，從而在實際預報過程中對臺風移動的路徑進行模擬[11]，這種方法的不足之處在于對實時預報數據的利用較少，且狀態轉移概率情況受限于研究區域。這些方法都只聚焦于歷史數據或實時預報數據，沒有同時兼顧到兩者。

本文使用實況臺風數據融合背景分析風場，驅動近岸波浪模式以及海洋模式進行浪-流耦合模擬，使用后報方式驗證數值模式；在數值模式精度良好的基礎上，提出一種新的臺風預報路徑集合生成方法，考慮研究區域內歷史臺風數據和業務中心實時預報數據，構建集合路徑并形成集合風場，開展近岸區域波浪要素計算的數值實驗，并考慮集合路徑情形下有效波高計算結果的覆蓋率，以此提供一種新的臺風浪概率數值預報方法。

1 研究區域和研究方法

1.1 研究區域與數據概況

長江口區域處于亞熱帶季風氣候區域，最大風速大多發生在夏秋季的臺風時期。該區域內的波浪以風浪為主[12]，由于入海口分叉較多，地形和岸線較復雜，區域內不同位置的波高差異性很大，最大波高大多出現在臺風期間。據統計，平均每年有2～3 個臺風對該區域產生直接影響[13]。長江口區域存在諸多重要的近岸站點（見圖1）。本文采用H1與H2兩浮標平臺的實測波高數據，采樣頻率為1 h，其中H1位于杭州灣北岸近岸區域（121°33′12.05″E，30°44′52.84″N），H2 位于長江口外30 m 水深處（122°48′29.29″E，30°58′12.19″N），兩個站位兼顧了外海與近岸的波浪差異及變化。

圖1 研究區域及主要岸基及浮標站點Fig.1 Study area and coastal and buoy observation sites

1.2 數值模式

本文選擇數值模式作為研究方法，波浪模式選用第三代海浪數值模式（Simulating Waves Nearshore，SWAN），考慮風能輸入、波浪折射、三波和四波相互作用、深度誘導波破碎、底摩擦等因素；風場輸入采用JANSSEN 風增長方式，白冠耗散使用KOMEN 模式，底摩擦采用JONSWAP 模式進行參數設置?？紤]到潮汐作用，除為波浪模式輸入風場外，還必須為其輸入潮流和潮位作為驅動要素，其由海洋模式進行計算，從而進行高精度的浪流耦合計算（見圖3）。

海洋模式選用有限體積海洋數值模式（Finite-Volume Community Ocean Model，FVCOM）。它是一種無結構三角形網格架構、有限體積自由表面的三維海洋數值模型，其原始方程包括動量、質量連續方程，以及溫度、鹽度和密度方程。和結構網格海洋模型相比，三角形網格在岸線相接處有局部加密效果，且過渡更加平滑，在長江口杭州灣區域已廣泛被使用[14-15]。

為提高預報時效并保證計算精度，SWAN 模式采用大小區域嵌套計算的方式。大區域覆蓋了研究區域內臺風從產生到對近岸造成影響的整個過程，能給小區域提供合理的計算邊界（見圖2a）。為兼顧計算效率，采用正交曲線網格，覆蓋范圍為23°～42°N，115°～130°E，包括了東海、黃海、渤海和部分西北太平洋洋面，開邊界處分辨率為3 km 左右。近岸區域，由于岸線和地形的復雜性很高，對波浪要素的影響十分顯著。為更加精細地刻畫各種局地變化，模擬近岸波浪，保證計算結果的準確性，小區域（見圖2b）采用高分辨率無結構三角形網格，主要包括長江口以及杭州灣，河流上邊界延伸至大通站，在河道以及沿岸處進行加密，分辨率為100 m 左右。在預報計算時，大區域計算較快，將其波浪譜結果作為小區域的計算邊界，并在小區域內進行精細化計算，從而在保證計算結果準確性的同時兼顧計算效率。浪-流耦合計算框架見圖3。

圖2 研究區域內SWAN模式計算網格Fig.2 SWAN model grid in the study area

圖3 浪-流耦合計算框架Fig.3 Wave-flow coupled computing framework

2 集合風場的構建

2.1 集合路徑的形成

本文中以24 h、48 h、72 h 3個預報時刻為代表，構建集合路徑并形成集合風場。集合路徑的構建采用以下步驟：首先形成某一時刻的預報位置集合，并為它們計算分配相應的權重，再將3個代表性預報時刻各自位置集合進行組合，形成所有的集合路徑，并計算出每條路徑的權重。整個流程的重點是對某一時刻各預報位置的確定和權重分配，本文采用如下方法：統計業務預報中心（中央氣象臺）歷史臺風位置預報誤差，并形成當前預報時刻位置偏差的概率分布情況，選取兩個特征概率從而確定兩個特征偏差距離；同時，使用臺風歷史最佳路徑集數據，基于面積相似法[16]優選出若干相似路徑并綜合距離相似及時間相似搜尋出預報時刻相似路徑上的對應位置，統計相似位置相對當前時刻預報位置的方向，這些方向代表著臺風移動方向以及移動速度的不同趨勢，使用聚類分析的方法確定若干個主要聚類方向以及各自的權重；將業務預報位置作為中心，偏差距離作為特征半徑，聚類方向結果作為特征方向，形成某一預報時刻的特征預報位置集合。

2.1.1 代表性位置偏差計算

本文搜集了中央氣象臺2008—2020 年東海區域共318場臺風的歷史預報數據以及對應的實際臺風過程數據（網址：http://typhoon.zjwater.gov.cn）；以實時起始時刻和預報時刻的位置形成的區域篩選出此區域內相同預報時效的歷史預報數據，統計各預報時效下歷史位置預報偏差的累計概率情況。預報偏差EDis為根據兩點經緯度計算出的距離誤差，公式為：

式中：R為地球半徑;a為兩點緯度之差；b為兩點經度之差；φ1和φ2為兩點緯度。

各預報時效位置偏差的累計概率曲線見圖4。從圖中可以看出，隨著累計概率的增大，各預報時刻位置偏差增大的速率都越來越大。本文選取70%和95%概率為典型代表作為特征概率，確定各預報時刻兩個特征偏差距離作為各預報時刻的兩個特征半徑，結果見表1。

表1 特征偏差距離（單位：km）Tab.1 Feature deviation distance(unit：km)

圖4 各預報時刻位置偏差的累計概率Fig.4 Cumulative probability of position deviation at each forecast moment

2.1.2 主要聚類方向計算

搜集了1949—2020 年西太平洋區域臺風最佳路徑集資料（網址：http：//tcdata.typhoon.org.cn）[17-18]，對實際臺風過程中的各預報時刻，基于面積相似法[16]，以歷史臺風軌跡和實時臺風軌跡上所有路徑點為分割點將整個面積分成L段，采用積分方式計算歷史臺風路徑和當前路徑圍成的面積，并作為相似度依據。計算公式為：

式中：S為總面積；Si為面積微元。

采用這種方式，當預報臺風、預報起始時刻以及預報時效發生變化時都會導致相似路徑成員發生變化。獲取所有相似路徑后，選擇時間和距離最近的對應預報時刻的相似位置，將業務預報位置和所有的相似預報位置投影為二維平面上的點，計算出所有相似位置相對于業務預報位置的方向角（0°～359°）。將這些方向角作為數據，使用一維均值聚類方法形成特征聚類方向。計算方法為：

式中：樣本共劃分k簇；Ci為第i簇；ui為簇Ci的均值向量，也稱為質心；x為樣本點；E為平方誤差。

式（3）中由于k值不同，計算出的E也不同。根據經驗，一般將k從2～10 賦值試算E，隨著k增大，E呈下降趨勢，以E在下降過程中梯度最大處的k值作為最優的聚類個數（簇個數）。每個簇質心的方向角代表一個特征方向，該簇包含的方向角數量在所有方向角數量中的占比為該特征方向的權重。

2.1.3 形成集合路徑

獲得各預報時刻的特征半徑和特征方向后，對某一特定預報時刻而言，以業務預報位置為中心可以形成特定時刻所有的特征成員位置；從每個預報時刻中各取一個成員位置，再連接所有選出的成員位置作為一條路徑，集合路徑的數量為各預報時刻成員位置進行組合的所有情況，每條路徑的概率也隨之確定（見圖5）。具體流程如下：

圖5 臺風集合路徑（0～48 h）Fig.5 Typhoon track ensembles(0～48 h)

①對特定預報時刻形成特征半徑及特征方向，獲得各特征方向的權重。以24 h 為例，基于歷史業務預報數據確定指定概率下特征半徑為Rt1；基于相似路徑數據，使用聚類分析方法收斂到最佳聚類個數Kt1（Kt1為式（3）中最優k取值，不同數據集最優值不同），每個聚類的中心作為一個特征方向，每個聚類樣本量占總樣本的比例作為特征方向的權重（Wt11,Wt12…,Wt1Kt1），對應的特征點空間位置為Pt11,Pt12…,Pt1Kt1。

②對24 h、48 h和72 h預報時刻，從每一個時刻選擇一個特征位置形成一條成員路徑，例如Pt1Kt1→Pt2Kt2→Pt3Kt3，相應地，這條路徑的權重為Wt1Kt1·Wt2Kt2·Wt3Kt3，3 個預報時刻一共可以形成Kt1·Kt2·Kt3條路徑，計算每條路徑的權重并進行集中構成集合路徑。

2.2 集合風場的形成

在具有確定性的臺風實況路徑的后報驗證中，將臺風風場模型融合背景風場作為風場輸入，而預報實驗中僅考慮臺風模型風場。

臺風風眼處風速較小且向外迅速增加，超過眼壁后風速又急劇減小直至消失。本文采用對稱風場模型刻畫臺風風場，其氣壓場由藤田公式[19]給出：

式中：r為距離臺風中心的距離；P為空氣質點的氣壓，Pe為環境氣壓，Pc為臺風中心氣壓；R為最大風速半徑。通過藤田公式可以計算出整個臺風風域內的氣壓分布情況。

臺風風速根據梯度風關系結合氣壓分布給出[20]，計算公式如下：

式中：f= 2Ωsinφ為科氏參數，Ω 為地球自轉角速度，φ為緯度；ρa為空氣密度。

在后報驗證中，本文使用CFSR 作為背景風場。CFSR 是第三代再分析產品，它是一個全球性的、高分辨率的、大氣-海洋-陸地-表面-海冰耦合系統，旨在提供這一時期這些耦合域狀態的最佳估計。該數據經緯度分辨率可以達到0.5°，時間分辨率可以達到1 h。兩種風場融合公式為[21]：

式中：WI為臺風模型風場計算風速；WQ為CFSR 風場風速；r為網格點與臺風中心的距離；e、c為系數；R為最大風速半徑，計算公式為[21]：

式中：Pc為臺風中心氣壓。

3 應用示例

本文以2019 年9 號臺風“利奇馬”為例，使用實況臺風路徑對數值模式進行后報驗證；選擇2019年8 月9 日00 時（北京時，下同）作為起始預報時刻進行72 h的虛擬集合預報實驗。

臺風“利奇馬”于8月7日05時被中央氣象臺升格為臺風，23 時升級為超強臺風，此后臺風一直向西北方向移動，靠近浙江沿岸地區，8 月10 日01 時45分左右在浙江省溫嶺市城南鎮沿海登陸，此時中心附近最大風力達到16級，中心最低氣壓為930 hPa，隨后臺風繼續移動至黃海海域，8月11日20時50分左右在青島市沿海再次登陸，此時中心最大風力仍有9 級，中心最低氣壓為980 hPa，此后臺風移動至渤海海面并不斷減弱直至消亡。

自起始預報時刻，構建24 h、48 h、72 h 的集合路徑。選取從臺風自洋面生成位置—各時刻業務預報位置區段作為實時軌跡，并與歷史軌跡做相似度計算，計算出24 h、48 h、72 h 預報時刻情形下相似度最高的各50 條路徑（見圖6）。從圖中可以看出，不同預報時刻的實時軌跡不同，相似路徑的選擇結果也不同。獲取相似路徑后，采用聚類分析方法確定各預報時刻的特征方向以及權重（見表2）。

表2 各預報時刻特征方向及權重Tab.2 Feature directions and weights for each forecast moment

圖6 各預報時刻的相似路徑及相似位置Fig.6 Similar tracks and similar positions for each forecast time

形成各預報時刻預報位置以及集合成員位置后（見圖7），采用2.1節中的方法形成48條（Kt1·Kt2·Kt3=4×3×4）集合路徑，并計算出每一條成員路徑的權重。

圖7 各預報時刻特征位置及實況路徑（0～72 h）Fig.7 Location of features at each forecast time and the real track(0～72 h)

考慮到模式冷啟動計算的不穩定性，FVCOM模式和SWAN 模式都提前7 d 開始運行。計算完成后，驗證部分站點有效波高的分布情況，給出各預報時刻下部分重點近岸站點以及整個研究區域有效波高的概率分布預測作為參考。

3.1 后報檢驗與數值預報實驗

從與H1、H2 站點實測波高數據的對比結果來看（見圖8），后報檢驗（實況路徑）的效果較好，自起始預報時刻～72 h 預報時段，H1 和H2 站點波高的均方根誤差平均值分別為0.31和0.60。使用集合路徑進行數值預報實驗的結果符合實際波高的基本發展趨勢，在臺風離岸較近時間段內（24～72 h），集合結果基本能夠覆蓋實際波高的波動范圍，且在一定程度上能夠反映出波高的極值。值得注意的是，H1 站點位于杭州灣靠內區域，在24～48 h 時段內，集合預報實驗結果都小于實際值（見圖8a），且兩種概率路徑（70%和95%）區間值的差異性小于H2 站點，而在48～72 h 時段內，波高集合的結果基本覆蓋了實際值；H2 站點位于長江口外，集合預報實驗在24 h 附近的波高預測偏大（見圖8b），在48 h 附近的波高集合結果基本覆蓋了實際值，在此時段內，95%概率路徑結果的覆蓋效果更好。

圖8 H1（a）與H2（b）站點有效波高的驗證Fig.8 Validation of effective wave height at H1(a)and H2(b)sites

為量化集合預報在兩特征概率情形下各預報時段內對實際有效波高的預測情況，我們對每個實測站位計算其覆蓋指數，將指定時段內處于集合預報預測區間范圍內的實測波高點數量占該時段內總實測波高點數量的比例作為指定時段的覆蓋指數。計算公式為：

式中：CR 表示覆蓋指數；Oin表示預報時段內觀測點含于集合預報實驗結果上下限范圍內；Oall為預報時段內所有觀測點的數量。作為集合預報結果覆蓋范圍的參考，CR 的波動范圍為0～100%，其值越大越好。CR 的計算結果見表3。從表中可以看出，在0～24 h，臺風離岸較遠，集合成員離業務預報位置較近，集合預報實驗的準確性非常依賴業務預報的準確性，集合結果的發散程度小，因此，該時段由于業務預報效果較差，導致兩站的概率路徑的CR都為0，其中H1站位集合結果區間范圍??；在24～48 h預報時段的CR 較高，兩概率路徑下的CR 分別為39%和52%，其中H2 站位的波高峰值結果偏大，兩概率路徑的CR分別為22%和54%；在48～72 h預報時段的兩概率路徑的CR有所上升，分別為74%和94%。

表3 各站位集合預報的覆蓋指數（單位：%）Tab.3 Cover rate of ensemble forecasts at each station(unit：%)

3.2 重要沿岸站點各預報時刻有效波高預測的概率分布參考

在實測站點數據驗證良好的基礎上，分析重要沿岸以及近岸浮標站點各預報時刻實況路徑有效波高的分布情況，并給出集合預報路徑相對實況路徑有效波高預測的概率分布參考，結果見圖9，其中站點自江蘇沿岸開始，經沿北支、南支到杭州灣北岸區域，外加部分口外站點（見圖1）。近岸站點波高預測概率分布在70%和95%概率路徑下表現出的特征有所不同。95%路徑在24 h 預報時刻顯示出較大的概率分布（見圖9b），這一現象在沿岸及近岸站點如H5、堡鎮、趙家溝、蘆潮港有所體現，24 h臺風距離站點較遠，95%路徑的部分集合成員距離上述站點更近；而70%路徑在48 h 預報時刻顯示出較大的概率分布（見圖9a），該現象在H2、綠華山等部分口外站點有所體現，48 h 臺風位于站點附近，70%路徑的部分集合成員距離上述站點更近；兩種概率路徑在72 h預報時刻給出的概率分布情況基本一致。從空間分布來看，H2、嵊山等離岸稍遠的站點在整個時間序列上表現出相對較大的波高預測值，而徐六徑、崇西、堡鎮等站點波高預測相對較小。

圖9 重要沿岸站點有效波高預測概率分布Fig.9 Predicted probability distribution of significant wave height at important coastal stations

3.3 有效波高空間預測的概率分布參考

為分析研究區域空間上有效波高的變化情況，圖10 和圖11 分別給出70%與95%概率路徑集合風場情形下各預報時刻有效波高的空間分布情況。圖中為分別計算所有集合成員有效波高并減去實況后報的結果，并根據權重加權分別累計差值大于0和小于0的增減量分布。

圖10 70%集合風場情形各預報時刻有效波高空間分布Fig.10 The spatial distribution of significant wave height at each forecast moment for the 70%ensemble wind field case

圖11 95%集合風場情形各預報時刻有效波高空間分布Fig.11 The spatial distribution of significant wave height at each forecast moment for the 95%ensemble wind field case

在臺風中心移向研究區域的時段，與實況后報相比，集合預報的波高結果在臺風中心附近的遠岸區域表現為正增量，而在近岸處表現為負增量（見圖10b、10c、11b、11c），這表明集合預報對波高空間場有重塑作用，導致空間上波高局地分布的差異。兩概率路徑情形的正增量最高都超過了4 m（見圖10b、11b），兩者負增量在近岸區域最大為-1 m 左右（見圖10c、11c）。當臺風中心漸遠時，較大概率路徑情形的正增量區域（見圖11e）相對較小概率路徑情形（見圖10e）偏移近岸研究區域更遠，最大值均為1 m 左右，而較小概率路徑情形的負增量區域總體相對偏移近岸研究區域更遠（見圖10f、11f）。當臺風中心離研究區域更遠時，兩者都不再有顯著的正增量貢獻（見圖10h、11h），且負增量空間分布的差異不大（見圖10i、11i）。從整個時序來看，空間中波高正增量區域基本和臺風中心位置保持一致，負增量區域滯后于臺風中心位置，該特性使得集合預報對部分口外站位波高峰值后的下降時段刻畫更好（見圖8b，24～72 h）。對比兩種路徑概率情形可以看出，范圍偏小的概率路徑（70%）集合預報波高加權結果的空間分布局地差異更大，波高隨空間分布的梯度更大。

4 結論

針對近岸區域波高受到很多因素的影響導致的預測難度較大的問題，本文提出了一種新型的臺風集合路徑的方法?；跉v史臺風數據和業務中心實時預報數據，采用統計學的方法給出任意實時情況下臺風各預測移動方向的概率，使用動態權重生成臺風預報路徑集合并進行數值預報實驗，得出實驗結果的概率分布。以長江口及杭州灣附近為研究區域，結合數值模式對臺風過境期間近岸區域的有效波高進行后報驗證，在驗證基礎上應用該新型集合路徑生成方法進行預報試驗，其中在48～72 h時段內，近岸區域H1 浮標集合預報結果的覆蓋指數能達到52%左右，H2浮標的覆蓋指數能達到94%左右。該方法可適用于臺風期間風暴潮、臺風浪等要素的預報，可為近岸海域的復雜動力要素預報工作提供新的集合預報技術。