吴甜甜，陈 煜，段忠东
（哈尔滨工业大学（深圳）土木与环境工程学院，广东深圳 518055）
热带气旋是最具破坏性的天气现象之一，发生在热带海洋上的一种具有暖中心结构的强烈气旋性涡旋，总是伴有狂风暴雨，常给受影响地区造成严重的灾害［1］。我国和东亚地区将这种强热带气旋称为台风，大西洋地区称为飓风，印度洋地区称之为热带风暴［2］。据统计，1979-2008年平均每年有7个台风登陆我国，期间每年由台风引起的直接经济损失最高达到了814.7亿元。尤其在经济发达且人口密集的沿海地区，台风灾害更应引起重视［3］。台风危险性的定量分析是灾害损失估计以及防灾减灾的基础，有助于科学地制定抗风减灾政策［4］。台风危险性分析是指，根据台风的气象特征以及台风的历史记录预测某一特定地点或区域在一定时期内可能遭受到台风的概率［5］。台风危险性分析在个别事件、季节预测和气候适应的背景下，对政府、工业界、金融机构、非政府组织，以及个人家庭都很重要。由于历史记录有限且早期观测质量不佳，台风灾害估计的常用方法包括基于历史观测资料的统计方法和基于数值模型模拟的方法。大多数巨灾模型（模型代表热带气旋危险性以及保险资产的经济损失）使用统计方法生成与历史数据相似的合成热带气旋［6］（例如，AIR WORLDWIDE，2015年）。这种方法通常模拟每个热带气旋的整个演变过程，包括起源、路径、强度和登陆。热带气旋造成的财产破坏和生命损失取决于其强度及登陆的位置，所以对路径和强度模拟尤为重要［6］。
在模拟热带气旋路径的同时，通常需要对热带气旋路径各个位置上的强度进行估算。目前估算强度的方法主要可分为样本统计法和数值模拟法。关于热带气旋路径模拟的误差在过去几十年中大大减少，然而对热带气旋强度模拟的提高则相对较少［7］。此外，由于影响强度变化的物理过程较为复杂，使得统计强度模型仍然具有一定的竞争力，在计算效率上也是台风危险性分析的首选。这也是为什么即使如今的台风天气数值模型已经有了较高的精度，而美国国家飓风中心（National Hurricane Center，NHC）仍然将统计模型作为飓风强度的预测手段之一［8-9］。
一些模型使用根据历史热带气旋数据训练的随机模型来模拟风暴的强度，这些随机模型可以是热带气旋的位置、风暴特征或环境参数的函数。SHIPS模型［10］利用多元线性回归将海气环境因子与强度变化联系起来。Rumpf等［11-12］假设热带气旋在相同区域具有相似的强度特征，统计得到路径各个位置的最大中心风速的概率密度分布，然后随机抽取最大风速。不过这种方法不考虑热带气旋受实际海洋表面温度、海气结构以及垂直风切等的影响。Vickery等［13］采用自回归模型并引入了海洋表面温度，其后进一步将海洋混合温度和风切效应引入自回归模型。这些模型对历史数据的使用或对环境参数的选择可能会限制其在当前气候下的应用。
Emanuel等［14］用二维耦合飓风强度预测系统（CHIPS）估算了北大西洋热带气旋的强度。与样本统计方法相比，虽然这种方法更具大气物理学基础，但是也要消耗更多的计算时间。随后Emanuel等［15］基于热带气旋强度演化的物理机制，提出了一个轴对称动力学模型来模拟热带气旋强度的演化。由于模型更完整地表示热带气旋的强度对环境的依赖性，使得该方法非常适合量化热带气旋强度在人为气候变化影响下的强度演化。Lee等［16］在Emanuel模型的启发下，开发一种新的统计动力学降尺度的多元自回归强度模型，该模型采用的方法是基于观察到的热带气旋强度的演化与大尺度环境场之间的关系，因此该模型较Emanuel模型更多地依赖于统计关系，但较Hall［17］的统计模型更具有物理意义。
由于台风危险性分析需要在短时间内模拟大量的台风强度样本，复杂的数值模型无法满足这个要求，显然统计强度模型在计算效率上是台风危险性分析的首选。传统的统计强度模型是直接将强度以及大量海气环境参数进行线性回归，但是线性回归模型往往需要拟合大量的环境参数才能有效模拟台风的强度演变［18］。
总的来说，热带气旋的强度会受到大气和海洋环境的影响，比如海表温度、垂直风切变、中层湿度等，热带气旋的强度与所处位置的大气与海洋结构有着密切的关系。因此，文中将进一步研究热带气旋强度与外部大气海洋环境之间的关系。为了克服已有线性回归方法的不足，并在热带气旋生命周期内强度的演变中引入增长控制项和其所能达到的最大风速的上限，文中提出了一个基于阻滞增长模型的热带气旋强度模型。该模型可以考虑部分台风强度演化的物理机制，它可以将台风强度的上限控制在最大潜在强度之内，而无需像线性回归模型人工截断台风的强度上限。此外，该强度模型的增长项中考虑台风强度与大气海洋环境环境的演变关系，从而进一步简化方程所需的环境参数。
1.1 数据来源
为了建立西北太平洋热带气旋的强度模型，模拟的热带气旋生成范围为100°～180°E，5°～35°N之内的海域，西北太平洋热带气旋主要发生在该区域，文中采用美国联合预警中心（Joint Typhoon Warning Center，JTWC）的热带气旋最佳路径数据集和欧洲天气预报中心（ERA-interim）的大气海洋环境参数数据作为建模的初始数据。
JTWC统计了所有进入西北太平洋海域（含南海，赤道以北，180°经度以西）的热带气旋的最佳路径数据集，统计数据包括热带气旋每6 h的中心位置、最大持续风速、中心最小压力、最大风速半径以及不同的风速半径等参数。由于1970年之后卫星观测技术的进步以及在气象学的普遍应用使得热带气旋的统计资料更加准确、全面、可靠，因此我们选择于1979-2018年生成的热带气旋作为样本数据建模，共计40年。由于统一将热带气旋初始生成位置定义为其首次达到15 m·s-1的最大持续风速的位置，筛选共计得到40年间西北太平洋区域台风有效记录发生数为1 137条。
模型选用的4个参数600 hPa气压层的相对湿度H（%）、海洋表面温度T（℃）、200 hPa气压层与850 hPa气压层之间的垂直风切变Vshear（m·s-1）、500 hPa的垂直速度ω（m·s-1），均来自ERA-Interim气候再分析数据的1979-2018年的月平均数据，该数据水平分辨率为1°×1°；
其中的大气环境风速和大气环境温度的垂直分辨率为1 000 hPa到50 hPa，共计20个垂直气压层。在模拟过程中将所需的大气和海洋环境参数插值到热带气旋中心所在位置。
1.2 基于阻滞增长模型热带气旋强度模型的建立
阻滞增长模型（Logistic人口模型）最早被应用于生物种群繁衍研究中，它考虑了自然资源、环境条件等因素对种群增长的阻滞作用。它是对马尔萨斯模型的改进，增加了一个附加项［19］。模型的阻滞作用体现在增长率α上，其微分方程的表达如式（1）所示，
式中：α是人口自然增长率；
β为生命系数；
N(t)为t时刻的人口总量；
N(t0)为初始时刻t0时的人口总量。微分方程右式的第一项代表增长项，第二项是考虑到现有资源有限的限制增长项［20］。事实上，式（1）是一个伯努利（Bernoulli）方程的初值问题［21］。在给定初始条件即可得到其解析解为：
当t→∞时，则N(t)→M=α/β，M称为环境承载力，表示环境可以容纳的最大种群数量。
Demaria［8］将阻滞增长模型应用于热带气旋强度的预测，其利用单个热带气旋从生成到本时刻的强度演变来拟合阻滞增长模型的系数，用于预测下一时刻的强度。不同于Demaria对单个热带气旋强度的预测，台风危险性更加关注的是热带气旋强度整体的演变特征，具体建模方法如下：
首先我们将式（1）中的因变量种群数量N(t)替换为热带气旋的强度V(t)，它是一个关于时间的函数。从风致表面热交换的不稳定性来看热带气旋的强化过程，其表面能量取决于当时刻的表面风速。因此，风致表面热交换的过程可由式（1）微分方程右式的第一项增长项表示［9］。由于热带气旋的强化过程不能无限期地继续下去，台风强度的大小存在上界即潜在强度，这个过程可用微分方程的第二项限制增长项来表示。考虑到台风发生发展的4个必要条件［1］：
（1）热力条件：海面水温必须在26℃以上，暖洋面有利于水汽和热量的向上输送，这是扰动形成暖心结构的基础，也是形成台风的主要能源。
（2）初始扰动：低层的初始扰动中，由于摩擦辐合产生的上升运动，可以使得大气的不稳定能量释放出来，转化为发展台风的动能。
（3）一定的地转偏向力：能使辐合气流逐渐形成强大的逆时针旋转的水平涡旋。由涡度方程得知扰动必须位于距赤道一定距离以外的地带才能产生较大的正涡度。
（4）对流层风速的垂直风切变要小：能使对流层上下空气的相对运动很小，可以很快地形成暖中心结构，最后形成台风。
热带气旋强度的演变与海洋表面温度、海气结构以及垂直风切强度等密切相关。因此，假定增长率α是一个与海气环境参数线性相关的函数，选取4个主要海气环境影响参数，即600 hPa气压层的相对湿度H（%）、海洋表面温度T（℃）、200 hPa气压层与850 hPa气压层之间的垂直风切变Vshear（m·s-1）、500 hPa的垂直速度ω（m·s-1）。于是，基于阻滞增长模型建立热带气旋的强度模型表示为：
因此给定某时刻条件下v(t=t0)=v0，可以得到热带气旋强度随时间变化的表达式为：
基于阻滞增长的强度模型估算海域热带气旋的强度的思路如下：将热带气旋当前时刻（t0时刻）的强度作为其在下一个时刻（t时刻）的强度的初值，即v0=v(t0)，然后通过迭代式（5）便可估算出热带气旋在各个时刻的强度。
为了确定式（5）中的参数(a1,a2,a3,a4,β)，将1979-2018年西北太平洋海域的4个海气环境参数的月平均数据、各时刻热带气旋强度的样本数据代入式中，通过多元非线性拟合得到其参数估计值。
热带气旋的强度是指热带气旋中心附近最大平均风力。从1989年开始，我国采用世界气象组织规定的统一标准，对热带气旋强度等级进行划分，如表1所示。
表1 热带气旋强度等级划分［22］Table 1 Intensity classification of tropical cyclone
图1 以台风“Tim”为例强度模型的划分Fig.1 Division of intensity model for typhoon"Tim"
考虑到热带气旋强度的演变是一个循序渐进的过程，呈现先逐步增强再逐渐减弱的大体趋势。因此文中将1979-2018年40年的JTWC最佳路径数据集的热带气旋强度样本，依据热带气旋强度是增强（发生发展阶段，即下一时刻强度大于等于上一时刻强度）还是减弱（衰减阶段，即下一时刻强度小于上一时刻强度）分为两大类，然后对于每个类别的样本数据依据台风强度等级划分为5类，分别为热带风暴强度（15～25 m·s-1）、强热带风暴强度（25～33 m·s-1）、台风强度（33～42 m·s-1）、强台风强度（42～51 m·s-1）、超强台风强度（大于51 m·s-1），建立10个阻滞增长强度模型。图1所示为强度划分的示例图（台风历史样本：编号1 994 060台风Tim），最后通过多元非线性拟合得到各模型的参数。
图2 1979-2018年台风实际强度和模型模拟强度的平均值及各时刻样本个数分布Fig.2 Typhoon actual and simulated average intensity and number of samples from 1979 to 2018
确定各模型的参数后，以15.4 m·s-1作为每条热带气旋初始时刻的强度值，并分别提取该时刻热带气旋所在位置的海气环境月平均数据，代入基于阻滞增长的强度模型中便可估算出下一时刻热带气旋的强度；
然后将估算的强度作为下一时刻的初始强度，提取该时刻热带气旋所处位置的环境变量，根据强度的增强或减弱趋势以及强度大小使用相应的强度模型估算下一时刻的强度，如此迭代便可模拟一场台风生命周期内的强度演变。
如图2所示为通过1979-2018年台风实际强度和模型模拟的结果，对比了热带气旋在其生命周期各个时刻整体强度的平均值（实际强度来源：JTWC最佳路径数据集）以及各时刻样本个数分布，图中的蓝柱表示各时刻的样本个数，可以发现大部分台风的生命周期在300 h以内。从结果直观来看，两条曲线的变化趋势和大小都较为接近，尤当样本数据量较大时（前120 h），模拟结果与再分析资料数据基本一致。在热带气旋发展390 h后，模型模拟的强度变化趋势与实际记录一致，但是大小存在一定误差，可能是数据样本数据过少导致系数拟合效果不佳。
通过计算均方根误差（RMSE）和平均误差百分比（MAPE）对模型进行客观的评价，结果如表2所示。结果表明，文中提出的强度模型对热带气旋在其390 h以前的强度模拟与实际记录较为吻合，390 h后的结果由于数据偏少，存在一定误差。
表2 1979-2018年台风实际强度和模型模拟强度的RMSE和MAPETable 2 RMSE and MAPE of typhoon actual and simulated intensity from 1979 to 2018
接下来，我们对各强度模型进行验证。图3～图7对比了各强度范围的实际和模拟热带气旋在其生命周期各时刻的整体强度均值，包括强度增长阶段和衰减阶段。事实上，从热带气旋整体平均强度演变来看，在300 h之内其强度从生成开始逐渐增强至最大强度，然后逐渐衰减。从图3（a）的热带风暴强度增长阶段模拟结果来看，在48 h之内的实际强度和模拟结果较为吻合，但是超过48 h模拟平均强度相较于实际结果偏大2～4 m·s-1。由于发展超过48 h的热带风暴其强度增加幅度一般不大，而模型中仍将其视作前48 h的热带风暴，类似于指数分布的无记忆性，因而模拟结果偏大。不同于热带风暴，图5（a）～图7（a）中的台风、强台风以及超强台风的生命前60 h模拟结果相较于实际偏低，这是由于台风，强台风以及超强台风生命早期样本数量较少，所以模拟早期快速增强的能力有所欠缺。由于热带气旋在其生命周期内任何时刻强度都有可能衰减，所以在衰减阶段强度的模拟，除了图7（b）超强台风早期衰减样本数量极少有较大误差外，其余衰减强度的模拟均与实际结果较为吻合。进一步，表3计算了各模型模拟结果的均方根误差（RSME）和平均误差百分比（MAPE），可以发现本模型估算的强热带风暴，台风和强台风的强度与实际记录较为吻合，其RSME均不超过2.5 m·s-1。对于热带气旋风暴和超强台风，本模型模拟的强度仍然有所不足，其MAPE在10%左右。
图3 热带风暴阶段模拟结果与实际平均强度的对比Fig.3 Comparison between simulated and actual mean intensity of tropical storm
图4 强热带风暴阶段模拟结果与实际平均强度的对比Fig.4 Comparison between simulated and actual mean intensity of strong tropical storm
图5 台风强度阶段模拟结果与实际平均强度的对比Fig.5 Comparison between simulated and actual mean intensity of typhoon
表3 1979-2018年台风实际强度和模型模拟的强度的RMSE和MAPE Table 3 RMSE and MAPE of typhoon actual and simulated intensity from 1979 to 2018
图6 强台风阶段模拟结果与实际平均强度的对比Fig.6 Comparison between simulated and actual mean intensity of strong typhoon
图7 超强台风阶段模拟结果与台风实际强度的对比Fig.7 Comparison between simulated and actual mean intensity of super typhoon
此外，从表3的对比结果来看，热带气旋强度衰减阶段各模型模拟结果均优于在强度增强阶段的模型模拟结果；
并且发现无论在强度增强阶段还是在强度衰减阶段，强热带风暴、台风以及强台风强度模型模拟结果误差较小，而热带风暴和超强台风强度模型模拟效果欠佳。由图7发现，超强台风强度模型误差主要来自前60 h，因为只有少量的热带气旋能在其生命周期的前60 h发展到超强台风，所以由于数据样本过少，导致拟合效果不佳。
我们通过分析模型的残差来考察模型假设的合理性及模型的可靠性。残差在数理统计中是指实际观察值与估计值（拟合值）之间的差值。图8所示为本模型的残差图。通过观测残差图，发现残差的分布基本符合正态分布，说明选用的模型比较合适，而且带状区域的宽度较窄，说明模型的拟合效果较好。为了进一步检验文中建立的强度模型对西北太平洋台风个例的模拟效果，这里我们利用5条历史台风记录作为检验案例，台风实例的相关信息如表4所示。
第1个模拟的案例是1979年的台风“爱丽斯（Alice）”，它于1月1日12时在太平洋海域（169.6°E，3.9°N）形成并发展，初始观测强度为15.4 m·s-1，随后逐渐向西南方向移动，于1月14日18时消亡。该台风的生成发展消亡均在太平洋海域，其强度经历了2次增强2次衰减。我们将初始强度15.4 m·s-1代入模型中，逐步估算台风生命周期内各个时刻的强度。模拟结果如图9所示，可以发现不仅模拟的强度大小和历史记录非常吻合，而且还模拟出了2次增强2次衰减的强度变化，模拟的平均绝对误差仅1.86 m·s-1。
图8 模型的残差图Fig.8 Model residuals
表4 5个台风案例说明Table 4 Description of five typhoon events
第2个模拟案例是台风“卢碧（Lupit）”，于2009年10月14日12时在太平洋海域（151.0°E，10.6°N）形成并发展，10月27日0时消亡。该台风的生成发展消亡均在太平洋海域。同样，我们将初始强度15.4 m/s代入基于阻滞增长模型的热带气旋强度，便可以逐步的估算台风生命周期内每个时刻的强度，模拟结果如图10所示，可以发现模拟结果和历史记录较为吻合，模拟的平均绝对误差为2.2 m·s-1。
上述检验的2个台风“温比亚”和台风“卢碧”仅在西北太平洋海域生成发展至消亡，并未登陆，接下来将验证几条登陆台风的模拟效果。
图9 台风“爱丽斯（Alice）”生命周期内的强度演变与模型模拟的结果Fig.9 Comparison between simulated and actual intensity of typhoon"Alice"in its life cycle
图10 台风“卢碧（Lupit）”生命周期内的强度演变与模型模拟的结果Fig.10 Comparison between simulated and actual intensity of typhoon"Lupit"in its life cycle
图11 台风“灿鸿（Chan-hom）”生命周期内的强度演变与模型模拟的结果Fig.11 Comparison between simulated and actual intensity of typhoon"Chan-hom"in its life cycle
图12 台风“山竹（Mangkhut）”生命周期内的强度演变与模型模拟的结果Fig.12 Comparison between simulated and actual intensity of typhoon"Mangkhut"in its life cycle
第3个模拟的案例是台风“灿鸿（Chan-hom）”，于2015年6月30日20时在太平洋海域（159.5°E，9.7°N）形成并发展，并于2015年7月11日16时登陆浙江省舟山市。我们将初始强度15.4 m·s-1代入基于阻滞增长模型的热带气旋强度模型中，便可逐步估算出台风生命周期内每个时刻的强度，模拟结果如图11所示，强度模拟的平均绝对误差1.80 m·s-1。
第4个模拟的案例是台风“山竹（Mangkhut）”，于2018年9月7日20时在太平洋海域（165.3°E，12.8°N）形成并发展，并于2018年9月15日登陆菲律宾北部，次日在广东台山海宴镇登陆。我们将初始强度15.4 m·s-1代入基于阻滞增长模型的热带气旋强度模型中，便可逐步估算出台风生命周期内每个时刻的强度，模拟结果如图12所示，强度模拟的平均绝对误差2.98 m·s-1。
通过对比5个历史台风案例发现该模型可以准确地模拟单个热带气旋在其生命周期内的强度变化，包括登陆之前或登陆后的热带气旋强度演变、长生命周期或短生命周期热带气旋的强度演变，均能模拟出强度的增强或减弱，模拟结果比较理想。但是对于有强度快速增强和快速减弱阶段的热带气旋，强度的演变模拟效果不佳。热带气旋强度的快速增强，Brand将其定义为在24 h内热带气旋的强度变化大于25.72 m·s-1［23］。如图13所示为台风“Vera”的实际强度演变与模拟结果的对比结果：在该台风发展的第30～60 h期间，实际强度从36～72 m·s-1，处于热带气旋强度快速增强阶段。而模拟的强度演变是从36～60 m·s-1，尤其在第42 h，误差达到最大为20 m·s-1（实际强度为69 m·s-1，模拟强度为49 m·s-1）。模型对于热带气旋强度的快速增强或减弱模拟能力有所欠缺。
图13 台风“Vera”生命周期内的强度演变与模型模拟的结果Fig.13 Comparison between simulated and actualintensityoftyphoon"Vera"initslifecycle
文中基于阻滞增长模型建立了热带气旋强度模型。类比种群动力学，将种群数量限制在零和一个上界之间，把阻滞增长模型应用于热带气旋强度模型。同时在阻滞增长强度模型中的强度增长率中纳入4个海气环境因素的影响，并依据热带气旋强度发展阶段以及强度等级划分为10类，分别建立阻滞增长强度模型。主要结论为：
（1）相较于数值模型，本模型模拟热带气旋强度演变时具有较高的计算效率；
相较于传统回归统计模型，仅需要较少的参数且无需人为截断热带气旋强度的上限。整体平均模拟结果与历史较为吻合。
（2）从整体来看，本模型模拟的热带气旋强度结果更接近于历史平均结果，当强度变化波动较大时模拟效果不佳，这可能是由于部分阶段的强度样本数据过少导致系数拟合不佳或强度模拟过程中未考虑误差项。
（3）通过对比1979-2018年西北太平洋海域台风实际强度和模型模拟的结果，发现本模型可以较为准确地模拟单个热带气旋在其生命周期内的强度变化，但对于强度快速增强或减弱阶段的热带气旋强度的模拟，效果仍有所欠缺。
未来，我们将从2个方面对强度模型进行改进：一是考虑不同生命阶段热带气旋的强度演变规律；
二是在模型中考虑强度的随机误差项的影响。