基于多变量灰色模型的台风强度模拟方法*

孙建鹏，郑仕豪，马萧岗，黄文锋

(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)西北太平洋地区台风强度平均值呈现中心强、四周弱的分布特征，变异系数呈现北低南高、西低东高的分布特征；强度模型模拟结果也体现出该分布特征，且在中国东部沿海及近海海域更为接近。

本文所建模型可为西北太平洋地区台风危险性分析提供强度参数，但由于本文研究的地域范围极广，且仅选取了两个主要环境因素建立模型，忽略了不同区域在环境因素与强度相互影响上的差异，在对具体某一区域进行计算时就会不可避免地产生误差，应用时还需依据实际情况具体分析。同时，该模型在模拟长生命周期台风后期阶段的强度变化上存在较大误差，这是该模型的主要缺陷之一，可以结合灰色系统理论，对模型结构作进一步优化，以实现对长生命周期甚至超长生命周期台风强度变化的准确模拟。