基于多變量灰色模型的臺風強度模擬方法*

孫建鵬，鄭仕豪，馬蕭崗，黃文鋒

(1.西安建筑科技大學 土木工程學院，陜西 西安 710055；2.合肥工業大學 土木與水利工程學院，安徽 合肥 230009)

西北太平洋海域是臺風事件的頻發地區，每年在該地區發生的臺風數量占全球總數36%以上；中國毗鄰西北太平洋，因此成為了世界上遭受臺風侵擾最為嚴重的國家之一。每次臺風事件的發生，都會給途經區域主要建筑物、構筑物的結構安全帶來挑戰，并造成重大經濟損失及人員傷亡[1-3]。因此，亟須對我國沿海地區及近海海域開展臺風危險性分析，對這些地區極值風速進行科學預測，為重點工程的工程抗風設計、防災減災政策制定等提供合理依據，以期把風致災害的損失降到最低。

目前，國內外學者對于臺風危險性分析的研究方法主要有三種[4]。一是以實測臺風強度為依據，通過選取合適的極值分布模型對歷史臺風最大風速序列擬合，得到單個研究點不同重現期內超越風速大小，該方法普遍應用于行業規范當中[4]。二是首先對研究點一定范圍內歷史臺風關鍵參數進行統計分析，通過Monte Carlo數值模擬與風場數值模型相結合的方式得到最大風速序列，最后計算出小區域范圍內極值風速[5-9]。三是第二種方法的進一步發展，該方法首先采用臺風全路徑模擬技術生成大量能夠反映實際分布特征的臺風樣本，然后結合風場數值模型推算出不同地區的極值風速大小[10-15]。與前兩種方法相比，第三種方法可以勝任全海域、大范圍臺風路徑模擬及全球氣候變化背景下的臺風極值風速預測工作，因此近些年來得到了更好的發展。隨機臺風強度的準確模擬是保證該方法計算精度的關鍵。

由于臺風強度變化的物理機制較為復雜，導致傳統的統計強度模型仍然具有一定的競爭力。HALL等[10]首先提取臺風歷史路徑中心最大風速時間序列，然后通過放縮手段，將其設置于隨機模擬路徑之上。RUMPF等[11]則假定同一地區強度特征相似，并對其進行統計，得到中心最大風速概率密度分布，并從中隨機選取最大風速值。以上兩種強度模擬方法均未考慮到環境因素對強度變化的影響。

VICKERY等[12]首先采用海洋表面溫度作為自變量，建立了自回歸強度模型；在此基礎上，VICKERY等[16]進一步在模型中考慮了海洋混合溫度和風切效應對臺風強度變化的影響。EMANUEL[13]則基于熱帶氣旋強度演變物理機制，建立了一種軸對稱動力學模型，該模型可用于量化人為氣候變化影響下的強度演化。LEE等[14]基于EMANUEL的建模思路發展了一種統計-動力學降尺度多元自回歸強度模型，該模型物理機制更加簡潔。陳煜等[15]在原有海洋表面溫度的基礎上基于BISTER和EMANUEL[17]提出的潛在強度理論，引入了潛在強度和相對強度，建立了二階滯后自回歸海洋強度模型。吳甜甜等[18]則基于生物種群研究中的阻滯增長模型，建立了基于阻滯增長模型的臺風強度數值模型。

雖然現階段主流的多元自回歸強度模型及阻滯增長強度模型相較于傳統的統計強度模型在模擬精度上有所提升，且更具物理意義，但該類強度模型模擬精度嚴重依賴于回歸模型的選取及參數擬合效果，這也是該類模型結構形式雖然不斷改變，但模擬精度并無較大提升的原因。本文將結合自回歸模型自身缺陷，采用多變量灰色預測優化模型，分析臺風強度與多個環境因素之間的灰色關聯，并選取與強度變化灰色關聯最強的幾個環境因素為模型自變量，建立起一個新的多變量灰色強度模型；并將該模型應用于西北太平洋地區臺風強度模擬工作，為該地區臺風危險性分析提供理論依據。

1 數據

1.1 數據來源

本文進行臺風強度模型建模所采用的歷史臺風強度記錄來源于中國氣象局(CMA)熱帶氣旋資料中心[19]，該數據集記錄了1949年以來發生在西北太平洋海域的所有熱帶氣旋行進路徑及各位置強度大小，相鄰記錄時間間隔為6 h，數據集詳細說明見表1。本研究所取時間段為1979-2018年，共計1 207條臺風數據。

西北太平洋地區海洋及大氣環境歷史數據呈現出明顯的季節性及年際變化特征。因此，本文所取環境參數歷史數據源自歐洲中期天氣預報中心(ECMWF)的ERA5全球再分析資料月平均數據。本文建模所涉及的環境參數包括：海洋表面溫度、相對濕度、環境風速、垂直速度、海洋表面氣壓、各氣壓層壓力值、各氣壓層溫度及水蒸氣混合比。該數據集所有環境參數的水平尺度(經向及緯向)分辨率均為0.25°×0.25°；垂直尺度方面，除海洋表面溫度及海洋表面氣壓外，其余參數均有37層氣壓等級(從1 hPa至1 000 hPa)。

1.2 數據預處理

為了統一臺風起點強度值，本文將所取時間段內所有臺風的初始強度值取為首次達到熱帶低壓強度等級(即10.8 m/s)時對應的強度值。同時，在臺風強度模型建模工作中發現，CMA所記錄的臺風存在一部分未命名臺風，這些臺風的強度演變趨勢較大一部分與臺風強度演變的一般趨勢不符，為了避免這一部分臺風強度記錄對強度模型的模擬效果造成影響，本文對這一部分未命名臺風記錄進行篩除，經預處理后得到502條臺風數據。

另外，我們取850 hPa和200 hPa氣壓層環境風矢量來表征對流層垂直風切變，具體計算公式如下：

(1)

式中：vshear為垂直風切變(m/s)；u850、u200分別為850 hPa和200 hPa氣壓層環境風速的徑向分量；v850、v200分別為850 hPa和200 hPa氣壓層環境風速的緯向分量。

2 方法

臺風強度演變是一個在多因素影響下的復雜系統，它們之間動力關系模糊且強度動態變化過程存在極強的隨機性，與灰色系統[21]特征相吻合。因此，臺風強度可以采用灰色系統理論相關模型進行模擬。

2.1 多變量灰色預測優化模型理論

多變量灰色預測優化模型(NSGM(1，N)模型)在傳統多變量灰色預測模型GM(1，N)模型的基礎上，利用差分模型代替影子方程，對傳統模型結構進行了優化，提高了模擬精度[22]。

該模型首先需對系統特征序列(因變量序列)X1(0)和與其相關性較高的n-1個相關因素數據序列(自變量序列)Xi(0)(i=2，3，…，n)求一階累加生成序列Xj(1)(j=1，2，…，n)，之后求得X1(1)的緊鄰均值生成序列Z1(1)，含一階差分方程及多個變量的新結構灰色預測模型結構表達式為(k=1，2，…，m)：

(2)

依據最小二乘法對上述模型進行參數估計，易證明模型參數列p=[b2，b2，…，b2，a，h1，h2]的最小二乘估計滿足：

(3)

其中：

(4)

(5)

進一步可以得到，當k=1，2，…，m時，NSGM(1，N)模型的時間響應式為：

(6)

NSGM(1，N)模型的累減生成式為：

(7)

式中：μ1=1/(1+0.5a)；μ2=(1-0.5a)/(1+0.5a)；μ2=h1/(1+0.5a)；μ4=(h2-h1)/(1+0.5a)。

2.2 灰色關聯度分析

過去幾十年的研究表明，臺風強度的發展變化取決于其周圍環境，這些環境因素主要包括海洋表面溫度、大尺度環流、中層大氣濕度、垂直風切變、海浪及飛沫等等，它們對熱帶氣旋強度變化的影響各有不同[23-26]。本研究選取與臺風強度演變關聯性最強的5個環境因素(海洋表面溫度、600 hPa氣壓層相對濕度、垂直風切變、500 hPa氣壓層垂直速度、最大潛在強度)作為待定自變量，并基于灰色關聯度理論[21]，根據式(6)計算它們與臺風強度值之間的灰色相對關聯度大小：

(8)

式中：Ri0(i=1，2，…，5)分別為海洋表面溫度、相對濕度、垂直風切、垂直速度及最大潛在強度與臺風強度值間的灰色相對關聯度；S0為臺風強度值序列經過無量綱化及始點零化像轉化之后，序列逐次映射到平面坐標系中與橫軸形成的閉合圖形面積；Si分別為上述5個環境因素序列經相同轉化后的閉合圖形面積。得到臺風計算結果分布情況(圖1)。

圖1 環境參數與臺風強度相對灰色關聯度分布

可以看出，垂直風切變和垂直速度這兩個參數與臺風強度的灰色相對關聯度要明顯大于海洋表面溫度、相對濕度和最大潛在強度。其中，垂直風切變與臺風強度關聯最為密切，平均灰色相對關聯度為0.771 8；垂直速度次之，平均灰色相對關聯度為0.745 1；其余三個參數與臺風強度關聯度均低于0.6。整體來看，垂直風切變和垂直速度與臺風強度的灰色相對關聯度相較于其余三個參數高出30%以上。造成這一現象的主要原因，可能是由于強度模型中環境參數取值為臺風中心位置處的參數值，這導致模型與垂直梯度參數關聯性相較水平梯度參數更強，而垂直風切變和垂直速度均為垂直梯度參數。

2.3 臺風多變量灰色強度模型建立

根據上述環境參數與臺風強度之間灰色相對關聯度分析結果，選取垂直風切變和垂直速度作為臺風強度模型的自變量，建立起NSGM(1，3)強度模型。

根據NSGM(1，N)模型的求解過程，可以推導出臺風發展至第k步(k=1，2，…，m；m為臺風總步數)模擬強度的差分模型為：

(9)

進一步，得到臺風發展至第k步的模擬強度計算值為：

(10)

3 結果與分析

3.1 經典案例模擬及檢驗

為體現NSGM(1，3)強度模型對某一具體臺風的強度模擬效果，我們選取2018年的第22號臺風“山竹(Mangkhut)”作為經典案例，對其強度變化過程進行模擬。該超強臺風于2018年9月7日20：00左右形成于西北太平洋洋面，9月11日12時左右中心附近最大風速達到最大值65 m/s并維持一段時間，9月16日17：00左右從廣東登陸，登陸時中心附近最大風速達到了14級[27]。該臺風的移動路徑及強度等級變化見圖2。

圖2 超強臺風“山竹”移動路徑及強度等級變化(注：該圖基于國家測繪地理信息局標準地圖服務網GS(2016)1561號標準地圖制作)

依據NSGM(1，3)強度模型的計算過程，通過該臺風活動期間各位置處的環境參數歷史數據，模擬出了該條臺風的強度變化過程(圖3)。

圖3 臺風“山竹”強度變化模擬結果

從強度變化趨勢上看，NSGM(1，3)強度模型對超強臺風“山竹”模擬與歷史保持了一致；臺風強度從生成起逐漸上升，在其生命周期第5 d強度達到最大，維持一段時間后強度逐漸衰弱；且模擬的臺風強度變化曲線較歷史更為平滑。從模擬誤差來看，殘差值在該臺風的整個生命周期內交替出現正負值，最大殘差絕對值出現在第42步時，其值為16.17 m/s，殘差絕對值的平均值為5.13 m/s，維持在較低水平；相對誤差在臺風發展的初始階段及末尾階段相對較大，最大相對誤差出現在第52步時，其他階段維持在較低水平，平均相對誤差為17.28%。總體來講，NSGM(1，3)強度模型較好地還原了臺風“山竹”的強度變化過程。

進一步，采用后驗差檢驗法，對臺風“山竹”的模擬效果進行檢驗分析，其步驟如下：

(11)

(12)

3)計算后驗差比值C及小誤差概率P：

(13)

(14)

4)對照后驗差精度表(表2)，確定模擬精度。

表2 后驗差精度對比表

經計算，采用NSGM(1，3)強度模型對臺風“山竹”強度模擬的后驗差比值為0.349 1，小誤差概率為0.960 8，根據兩者推斷出對臺風“山竹”強度的模擬精度達到了最高等級：一級(好)。

3.2 西北太平洋地區整體模擬及檢驗

利用NSGM(1，3)強度模型對西北太平洋地區所有臺風進行強度模擬，采用最小二乘估計求得該地區強度模型的系數向量，各系數具體估計值見表3。

表3 NSGM(1，3)臺風強度模型系數表

圖4為1979-2018年臺風歷史強度與模擬強度對比圖，對比了臺風在其生命周期內各步時(6 h/步)整體強度的平均值以及對應歷史記錄數量。可以看出，歷史平均強度變化總體上存在兩個拐點(圖4劃線標識處)，將整個變化過程清晰地分為三個階段：第一階段強度快速增強并達到峰值、第二階段強度緩慢下降、第三階段強度曲折波動并最終下降至最小值；兩個拐點分別出現在第24步時(第6 d)和第56步時(第14 d)。這三個階段的臺風歷史記錄數量也在逐漸減少；在第三階段，歷史記錄數量已經低于30條臺風，計算平均強度時第三階段的樣本數量最少，這也是導致平均強度值在該階段波動最大的主要原因。第一階段模擬效果相對最好，模擬強度值與歷史強度接近；第二階段模擬強度明顯高估了臺風強度；第三階段模擬結果與歷史則存在較大差異。總體來看，模型較好地還原了臺風強度變化趨勢。

圖4 1979-2018年西北太平洋地區臺風強度模擬對比

為了進一步量化NSGM(1，3)強度模型模擬效果，分別計算出模擬結果的均方根誤差(RMSE)和平均絕對百分比誤差(MAPE)，結果見表4。

表4 NSGM(1，3)臺風強度模型各階段模擬誤差計算表

可以看出，模型在第一階段的模擬效果最好，該階段RMSE為1.57 m/s，MAPE為5.47%，均保持在整個生命周期最低水平；至第二階段結束，模擬誤差增大，RMSE和MAPE分別增加至3.18 m/s、10.43%；第三階段誤差則遠大于前兩階段誤差。導致這一現象的原因主要來源于后期階段臺風樣本數量的快速下降，這也說明模型對長生命周期臺風強度演變的后期階段模擬效果較差。從模型的整體誤差來看，RMSE為6.70 m/s，MAPE為22.23%，均維持在較低水平，模型整體模擬效果良好。

為了更直觀地反映出西北太平洋海域強度模擬整體效果，我們將西北太平洋區域(0°～50°N、100°～180°E)網格化，分辨率大小為2.5°×2.5°，分別統計各網格內CMA最佳路徑數據集記錄臺風以及模擬臺風強度的平均值和變異系數大小，并進行對比分析。統計得到的各網格內CMA最佳路徑數據集歷史記錄和模擬的臺風強度平均值和變異系數如圖5和圖6所示。從臺風強度平均值的分布特征來看，本文提出的強度模型很好地還原了西北太平洋地區的臺風強度平均值分布特征。強度平均值整體呈現出中心強、四周弱的分布特征；強度平均值高值區域的緯度區間主要集中在北回歸線±10°N，和120°～150°E范圍之內；外圍區域強度平均值相對較小，這一特征沿著中國東部海岸線表現最為明顯。從臺風強度平均值大小來看，強度模型模擬結果相較于歷史值則普遍偏低，尤其在中心的高值區域；外圍低強度區域模擬值偏低情況則在低緯度地區(10°N以南)表現更為明顯。

圖5 西北太平洋海域臺風強度平均值分布圖

圖6 西北太平洋海域臺風強度變異系數分布圖

從強度變異系數的分布特征來看，強度模型也較好地還原了西北太平洋地區的臺風強度變異系數分布特征，該地區強度變異系數整體呈現出北低南高、西低東高的分布特征。強度變異系數的高值區域集中分布在5°～25°N、125°～165°E范圍之內，西北太平洋地區的北部區域及中國東部海岸附近的強度變異系數相對偏低。從強度變異系數大小來看，模擬結果在各區域均存在差異，但差異不大，模擬結果可以反映西北太平洋地區的強度變異系數大小整體情況。

4 結論及討論

本文基于多變量灰色預測優化模型，建立了一種臺風多變量灰色強度模型；利用該模型對臺風“山竹”的強度演變過程進行了還原，并對西北太平洋地區整體臺風強度的變化趨勢及分布特征進行了模擬，同時與歷史數據進行對比分析。研究的主要結論如下：

1)垂直風切變和垂直速度與臺風強度變化的灰色關聯最為密切，兩者與臺風強度灰色相對關聯度分別達到0.771 8、0.745 1。

2)NSGM(1，3)強度模型對臺風“山竹”強度模擬的后驗差比值為0.349 1，小誤差概率為0.960 8，后驗差檢驗模擬精度等級為一級，模型能夠很好地還原單個臺風的強度變化。

3)NSGM(1，3)強度模型的整體模擬效果較好，尤其是在臺風強度變化的第一階段，但該模型對長生命周期臺風后期階段強度變化模擬效果較差。

4)西北太平洋地區臺風強度平均值呈現中心強、四周弱的分布特征，變異系數呈現北低南高、西低東高的分布特征；強度模型模擬結果也體現出該分布特征，且在中國東部沿海及近海海域更為接近。

本文所建模型可為西北太平洋地區臺風危險性分析提供強度參數，但由于本文研究的地域范圍極廣，且僅選取了兩個主要環境因素建立模型，忽略了不同區域在環境因素與強度相互影響上的差異，在對具體某一區域進行計算時就會不可避免地產生誤差，應用時還需依據實際情況具體分析。同時，該模型在模擬長生命周期臺風后期階段的強度變化上存在較大誤差，這是該模型的主要缺陷之一，可以結合灰色系統理論，對模型結構作進一步優化，以實現對長生命周期甚至超長生命周期臺風強度變化的準確模擬。