CN112883635B - 一种基于随机森林算法的热带气旋全路径模拟方法 - Google Patents

Abstract

本发明提供了一种新型的基于随机森林算法的热带气旋全路径模拟方法。采用基于机器学习的随机森林算法挖掘历史台风关键信息的非线性映射关系，能有效避免传统全路径模拟方法的线性回归分析拟合能力的不足。该方法对基于聚类分析生成的各类路径的整体信息建立模型，可以针对每个模型进行调试以保证学习效果，无需将计算海域划分成网格以进行大量的回归分析。随机森林算法技术手段的引进可一定程度上改进传统方法参数标定没有标准流程的缺陷。与传统模拟方法一致的是，该方法亦能对整个西北太平洋区域的热带气旋生成、行进路径、行进方向、洋面强度发展以及登陆强度衰减进行随机模拟，生成大量符合历史样本特征的热带气旋路径及强度随机事件样本集，据此进行中国沿海地区台风灾害风险的定量精细化评估，从而有效地解决局部区域历史样本不足的问题，为灾害风险的定量精细化评估奠定了基础。该方法能保证较好的准确性和普遍性，可推广到其他区域进行飓风灾害风险评估分析。

Description

一种基于随机森林算法的热带气旋全路径模拟方法

技术领域

本发明属于工程结构抗风性能设计领域，尤其涉及一种基于随机森林算法的热带气旋全路径模拟方法。

背景技术

台风过境常常带来狂风暴雨天气，引起海面巨浪并形成风暴潮危害，严重威胁沿海工程结构安全。据统计，每年平均登陆我国东南沿海的热带气旋包括台风及超强台风的数目约为 9个，虽然西北太平洋热带气旋总数近年有减少的趋势，但气旋强度呈现增大趋势。随着沿海地区经济的迅速发展，人口和建筑密度的迅速增大，随着愈来愈多的风敏感结构如超高层建筑、大跨体育场、大跨桥梁等的兴建，台风引起的工程事故屡见不鲜。因此，鉴于台风灾害可能造成的巨大损失，工程结构抗风性能设计有必要考虑台风作用。但是，目前国内结构设计规范关于风荷载的平均风和脉动风特性部分都是基于空旷场地良态风的统计成果，并未对台风风场作出明确规定。台风风场相对于良态风存在尺度小、频率低、数据采集难度大等特征，导致目前近地台风实测数据尚不充足，无法对台风气候模式作出准确标定。因此，针对沿海重大工程结构设计进行台风灾害风险评估以及重现期台风风荷载的合理估算须基于有限的气象资料，借助合理的数值模拟手段实现。

在工程设计领域，面向灾害风险评估的热带气旋数值模拟方法主要包括模拟圆法和全路径模拟法。模拟圆法基于单个站点或是局部区域范围的过境热带气旋实测数据，应用台风风场模型通过数值计算模拟台风风速。基于局部区域范围的模拟圆法操作步骤少，计算效率高，同时也存在着局限性，即需要统计得到模拟点处各个台风关键参数的概率模型，当高纬度地区台风数据较少时该方法建立的概率模型存在较大不确定性。为应对模拟圆法的不足， Vickery提出了一种全路径模拟方法，即对整个海域进行分区，针对每小块区域内的台风信息作线性回归分析，能够模拟整个洋面的热带气旋发展全过程。该方法在北大西洋和西北太平洋的热带气旋灾害评估中得到一定的应用，也被用于参与制定美国ASCE7设计规范的基本风速。但是Vickery的线性回归模型需要经过大量的模型系数调试以使模拟气旋的关键参数统计值与历史气旋的关键参数统计值相一致。参数调试的方法较为主观，无法形成统一的标准，难以广泛推行。

发明内容

本发明的目的在于，针对现有技术中存在的不足，提供一种能实现精细化评估、结果精确可靠、标准化程度高、适用范围广的面向灾害风险评估的基于随机森林算法的热带气旋全路径模拟方法。

为此，本发明的上述目的通过以下技术方案实现：

第一步，建立起始点模型，包括对历史热带气旋年发生率进行负二项分布拟合并随机生成模拟热带气旋的年发生次数，以及统计热带气旋历史数据样本信息并随机抽样生成起始点信息；

第二步，建立行进模型，基于历史热带气旋的行进特征将其分为四类，基于随机森林回归算法对各类历史气旋的经纬度信息进行学习训练，生成随机森林行进模型，模拟预测热带气旋的经度和纬度以形成热带气旋路径；

第三步，建立强度模型，基于随机森林回归算法针对历史热带气旋的强度影响因子和相对强度数据进行学习训练，生成随机森林洋面强度发展模型。陆地强度衰减模型采用指数衰减函数，并根据热带气旋路径中下一点的是否位于洋面选择相应的强度模型；

第四步，终止判别，基于随机森林分类算法对不同路径类别的历史气旋的终止条件进行学习训练，并判定模拟路径的模拟点是否终止；

第五步，重复前四步，直至生成给定模拟年的所有台风路径，并对路径模拟结果进行检验，对随机生成的热带气旋路径及强度大样本的关键参数进行统计结果检验以及统计规律分析；

第六步，统计进入模拟点250km范围内的台风路径和强度信息，建立台风风场工程模型和边界层模型，包括风场模型关键参数的标定，计算模拟点台风风速，以进行极值分布分析。

在采用上述技术方案的同时，本发明还可以采用或者组合采用如下技术方案：

作为本发明的优选技术方案：第一步中采用负二项分布对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。根据西北太平洋区域的历史热带气旋年发生次数数据特点，对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。起源于洋面上的台风年发生次数宜用负二项分布模拟。

作为本发明的优选技术方案：第一步中基于历史的热带气旋起始点信息计算全路径模拟所需的信息：所在位置的洋面温度、相对强度，进而通过随机抽样确定每一条模拟台风的起始点信息。

作为本发明的优选技术方案：首先对西北太平洋热带气旋的行进特征进行聚类分析，将所有路径样本分成四类，分别计算并统计各类气旋的基本信息，包括各数据点的经纬度、中心最低气压、洋面温度、相对强度等。路径类型包括，类型1：从生成至终止均在洋面且主要在东西方向运行分布；类型2：主要行进路线为西进方向，可能途径菲律宾、中国南海、中国东南沿海诸省和东南亚部分区域；类型3：初始行进路线为西北方向，之后转向东北方向行进，可能途径中国东部沿海诸省、日本和朝鲜半岛；类型4：从生成至终止均在洋面且与类型1相比，热带气旋更多地向北发展。

作为本发明的优选技术方案：第二步中分别对各类数据的经纬度信息进行回归分析。采用基于机器学习的随机森林回归算法，对下一步经纬度做出预测：

式中：RFr_lat和RFr_lon分别为用于纬度预测和经度预测的随机森林回归模型；j＝1,2,…,n，表示采用的路径类别，n表示路径类别总数，本发明中n＝4；lat为纬度，lon为经度；i+1、i、 i-1、i-2、i-3分别表示下一、当前、前一、前二、前三时间步；ε_lat和ε_lon为符合正态分布的零均值随机扰动项。

作为本发明的优选技术方案：第三步中对热带气旋路径中下一点的经纬度进行判断，如果下一点是位于陆地，则选择陆地强度衰减模型，如果下一点是位于洋面，则选择洋面强度发展模型。

作为本发明的优选技术方案：第三步中洋面强度发展模型是根据随机森林回归算法建立洋面强度回归模型，预测洋面上下一点的相对强度，从而确定中心气压：

I_i+1＝RFr_I(lat_i,lon_i,I_i,I_i-1,I_i-2,SST_i+1,SST_i,SST_i-1,SST_i-2,MSLP_i,MSLP_i-1,MSLP_i-2)+ε_I (2)

式中：RFr_I为用于热带气旋强度预测的随机森林回归模型；I为相对强度；SST为洋面温度 (K)；ε_I为符合正态分布的零均值随机扰动项；MSLP为近中心最小海平面气压(hpa)，综合考虑洋面温度、平流层温度、相对湿度等环境因素影响，将中心气压转换为相对强度进行预测。

作为本发明的优选技术方案：第三步中陆地强度衰减模型利用登陆强度指数衰减模型，模拟热带气旋登陆后的强度衰减：

Δp_t＝Δp₀exp(-mt) (3)

式中：Δp₀为登陆时的中心气压差；Δp_t为登陆t小时之后的中心气压差； m＝m₀+m₁Δp₀+ε_m，系数m₀、m₁和随机扰动项ε_m可通过历史登陆热带气旋的气压差变化回归分析得到，综合考虑登陆持续时间、热带气旋强度大小以及地表粗糙度等因素影响，建立登陆强度回归模型，模拟热带气旋登陆后的强度衰减。

作为本发明的优选技术方案：第四步中采用随机森林分类算法对历史数据进行学习训练，对四种路径类别的热带气旋进行终止判别，判断下一点是否为终止点：

end_i+1＝RFc_end,j(lat_i,lon_i,SST_i,MSLP_i,N_i) (4)

式中：end表示热带气旋是否终止的判决结果，属于逻辑标量，“1”为终止，“0”，为不终止；RFc_end表示用于终止判定的随机森林分类模型；N表示热带气旋已行进的时间步。

作为本发明的优选技术方案：第五步中模拟结果检验是对比随机模拟的路径及强度与历史热带气旋数据之间的差异，定量刻画模拟结果的精度。

作为本发明的优选技术方案：第五步中的统计规律分析是将中国东南沿海海岸线划分为若干海岸站点，再统计每一个海岸点250km范围内的模拟和历史热带气旋年均频次、移动速度、朝向以及中心气压差四个关键参数的结果差异，从一定程度上可以反映热带气旋全路径模拟的空间分布效果。采用Kolmogorov–Smirnov test(K-S检验)对历史样本和模拟样本移动速度、朝向以及中心气压差三个关键参数的经验分布进行比较。

作为本发明的优选技术方案：第六步中对于进入模拟点250km范围内的台风路径以及强度样本，采用Yan Meng风场模型计算模拟点的台风风速，具体步骤如下：

(1)台风风场模型关键参数的标定

采用Vickery经验模型估计最大风速半径R_max和Holland气压剖面参数B：

式中：Δp为中心气压差；ψ为纬度；f_c＝2×7.273×10^-5sinψ；

和ε_B为符合正态分布的零均值随机扰动项；

(2)建立气压模型

采用Holland气压模型：

p＝MSLP+Δp exp[-(R_max/r)^B] (7)

式中：p为距台风中心径向距离为r处的海平面压强；R_max为最大风速半径；B为气压剖面参数；

根据式(5)，气压梯度表达式为：

(3)建立平衡方程

中性大气条件下的Navier-Stokes方程为：

式中：V为空气微团的运动速度，Yan Meng风场模型把V视为梯度风速V_g和地表摩擦风速V′两部分矢量叠加，即V＝V_g+V′；f为Coriolis参数；F为边界层摩擦力；

由于径向风压梯度在台风边界层内随高度变化很小，再忽略边界层以上的摩擦力，式(9) 可以分解为梯度层和边界层两部分，分别如式(10)和式(11)：

结合

以及

式(10)和式(11)可进一步化简为：

V′·▽V′+V′·▽V_g+V·▽V′＝-fk×V′+F (13)

在随台风中心移动的二维极坐标中，沿径向和切向分解式(10)得到：

式中：V_θg和V_rg分别是切向和径向的梯度风速；c_r＝-V_T cos(θ-β′)，c_θ＝-V_T sin(θ-β′)；ρ是空气密度，取1.2kg/m³；θ是正东方向与模拟点和台风中心连线夹角，以逆时针旋转为正；β′是正东方向与台风移动方向的夹角，以逆时针旋转为正；

在上层大气环境中，考虑到V_rg远远小于V_θg，忽略式(14)中的前两项，可以求出V_θg；由于V_rg值很小，近似取0仍可满足风速模拟的准确性，因而，演化出梯度风速的解析式，如下式所示：

在近地边界层中，地表摩擦风速的切向和径向分量V′和V_r′小于相应的梯度风速分量，因此，这些分量对θ的一阶导数可假定同样小于相应的梯度风速分量的一阶导数；于是式(11) 可作线性化分解，如下式所示：

上层大气和地表附近的边界条件分别为式(18)和式(19)：

V′|_z′→∞＝0 (18)

式中，V_s表示近地面风速，计算起始点取为z＝h+10即z′＝0处，考虑边界条件后，式(17) 可以求解为：

V″＝Dexp[-(1+i)λz′] (20)

式中：边界层地表复常量D＝D₁+iD₂，引入式(21)所示的无量纲参数ξ和带量纲参数λ，最终可推导出地表摩擦风速的解析式，如式(22)所示：

D₁和D₂可按下式计算：

式中：C_d为阻力系数；k_m为运动黏度，取100m²/s；k为Karman常数，取0.4；平均粗糙单元高度

实测得A＝11.4；零平面位移d＝0.75h；z₁₀设定在平均粗糙单元h高 10m处；z₀为考虑地形和地表粗糙度影响而引入的等效粗糙长度；

(4)计算流程

Yan Meng台风风场模型计算地表风速的流程为：首先在二维极坐标内由式(16)计算梯度风速，然后取梯度风速作为地表风速的初始值，依次代入解析公式(23)、(21)和(22)求得地表摩擦风速，由梯度风速和地表摩擦风速叠加得到新的地表风速，经过多次迭代，直到收敛；

(5)台风极值风速预测及灾害风险评估

通过数值模拟得到台风极值风速序列后，需要用极值概率分布模型对该序列进行拟合，进而预测重现期极值风速；常用的极值概率分布有极值I型分布、极值II型分布和极值III 型分布；当模拟的台风样本足够多时，应当首选经验分布。

本发明提供了一种新型的基于随机森林算法的热带气旋全路径模拟方法。采用基于机器学习的随机森林算法挖掘历史台风关键信息的非线性映射关系，能有效避免传统全路径模拟方法的线性回归分析拟合能力的不足。该方法对基于聚类分析生成的各类路径的整体信息建立模型，可以针对每个模型进行调试以保证学习效果，无需将计算海域划分成网格以进行大量的回归分析。随机森林算法技术手段的引进可一定程度上改进传统方法参数标定没有标准流程的缺陷。与传统模拟方法一致的是，该方法亦能对整个西北太平洋区域的热带气旋生成、行进路径、行进方向、洋面强度发展以及登陆强度衰减进行随机模拟，生成大量符合历史样本特征的热带气旋路径及强度随机事件样本集，据此进行中国沿海地区台风灾害风险的定量精细化评估，从而有效地解决局部区域历史样本不足的问题，为灾害风险的定量精细化评估奠定了基础。该方法能保证较好的准确性和普遍性，可推广到其他区域进行飓风灾害风险评估分析。

附图说明

图1为本发明所提供的方法的流程图。

图2为历史热带气旋年发生次数年际变化图。

图3为历史热带气旋年发生次数概率分布图。

图4a为历史热带气旋密度空间热度图。

图4b为模拟热带气旋密度空间热度图。

图5a为CMA历史路径与模拟路径在各个站点的年发生率对比。

图5b为CMA历史路径与模拟路径在各个站点的朝向对比。

图5c为CMA历史路径与模拟路径在各个站点的移动速度对比。

图5d为CMA历史路径与模拟路径在各个站点的中心气压差对比。

图6为历史路径与模拟路径各关键参数经验分布的K-S检验结果图。

图7为深圳地区台风模拟所得10000年极值风速序列。

图8为深圳地区台风年极值风速经验分布图。

具体实施方式

下面结合具体实施例来对本发明进行进一步说明，但并不将本发明局限于这些具体实施方式。本领域技术人员应该认识到，本发明涵盖了权利要求书范围内所可能包括的所有备选方案、改进方案和等效方案。

参见图1，一种面向灾害风险评估的基于随机森林算法的热带气旋全路径模拟方法，其具体步骤包括：

第一步，建立起始点模型，包括抽样生成模拟热带气旋的年发生次数和生成起始点信息；具体的采用负二项分布对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。根据西北太平洋区域的历史热带气旋年发生次数数据特点，对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。

基于历史的热带气旋起始点信息计算全路径模拟所需的信息：移动速度、朝向、所在位置的洋面温度、相对强度，进而通过随机抽样确定每一条模拟台风的起始点信息。对于行进模型和终止判别，需要将路径分类后在针对各类路径进行回归分析并预测，分类依据为热带气旋的分布与行进特征。路径类型包括，类型1：从生成至终止均在洋面且主要在东西方向运行分布，分布区域在西北太平洋南部海域；类型2：主要行进路线为西进方向，可能途径菲律宾、中国 南海、中国东南沿海诸省和东南亚部分区域；类型3：初始行进路线为西北方向，之后转向东北方向行进，可能途径中国东部沿海诸省、日本和朝鲜半岛；类型4：从生成至终止均在洋面且与类型1相比，热带气旋更多地向北发展，分布区域在西北太平洋全部海域。

第二步，基于随机森林回归算法分别对各类路径数据的经纬度信息特点进行学习训练，生成路径预测回归模型。基于前三步的经纬度信息对下一步经纬度做出预测：

式中：RFr_lat和RFr_lon分别为用于纬度预测和经度预测的随机森林回归模型；j＝1,2,…,n，表示采用的路径类别，n表示路径类别总数，在本实施例中，n＝4；lat为纬度，lon为经度； i+1、i、i-1、i-2、i-3分别表示下一、当前、前一、前二、前三时间步；ε_lat和ε_lon为符合正态分布的零均值随机扰动项。

第三步，建立强度模型，包括洋面强度发展模型和陆地强度衰减模型，并根据热带气旋路径中下一点的经纬度选择相应的强度模型，具体对热带气旋路径中下一点的经纬度进行判断，如果下一点是位于陆地，则选择陆地强度衰减模型，如果下一点是位于洋面，则选择洋面强度发展模型。

对于洋面强度发展模型，综合考虑洋面温度、平流层温度、相对湿度等环境因素影响，将中心气压转换为相对强度。因为气旋强度特征与运行规律关系并不明显，因此采用目标区域内所有路径点样本信息进行训练学习，根据随机森林算法建立洋面强度回归模型，基于前几步的相对强度、洋面温度、中心气压及经纬度信息预测洋面上下一点的相对强度，从而确定中心气压：

I_i+1＝RFr_I(lat_i,lon_i,I_i,I_i-1,I_i-2,SST_i+1,SST_i,SST_i-1,SST_i-2,MSLP_i,MSLP_i-1,MSLP_i-2)+ε_I(2)

式中：RFr_I为用于热带气旋强度预测的随机森林回归模型；I为相对强度；SST为洋面温度(K)；ε_I为符合正态分布的零均值随机扰动项；MSLP为近中心最小海平面气压(hpa)。

对于陆地强度衰减模型，综合考虑登陆持续时间、热带气旋强度大小以及地表粗糙度等因素影响，建立登陆强度回归模型，模拟热带气旋登陆后的强度衰减，一般认为主要跟登陆时间、登陆纬度及登陆中心气压差有关：

Δp_t＝Δp₀exp(-mt) (3)

式中：Δp₀为登陆时的中心气压差；Δp_t为登陆t小时之后的中心气压差； m＝m₀+m₁Δp₀+ε_m，系数m₀、m₁和随机扰动项ε_m可通过历史登陆热带气旋的气压差变化回归分析得到。

第四步，针对各类别路径分别采用随机森林分类算法对历史数据的终止信息进行学习训练，判定各路径类别的热带气旋下一模拟点是否为终止点：

end_i+1＝RFc_end,j(lat_i,lon_i,SST_i,MSLP_i,N_i) (4)

式中：end表示热带气旋是否终止的判决结果，属于逻辑标量，“1”为终止，“0”，为不终止；RFc_end表示用于终止判定的随机森林分类模型；N表示热带气旋已行进的时间步。

第五步，重复上述步骤直至生成给定模拟年的所有路径，并对路径模拟结果进行检验。对随机生成的热带气旋路径及强度大样本关键参数进行统计结果检验以及统计规律分析；具体的检验是对比随机模拟的路径及强度与历史热带气旋数据之间的差异，定量刻画模拟结果的精度。

统计规律分析是将中国东南沿海海岸线划分为若干海岸站点，再统计每一个海岸点250 km范围内模拟的和历史的在热带气旋年均频次、移动速度、朝向以及中心气压差的结果差异，从一定程度上可以反映热带气旋全路径模拟的空间分布效果。采用K-S检验对历史样本和模拟样本移动速度、朝向以及中心气压差三个关键参数的经验分布进行比较。

第六步，建立台风风场工程模型和边界层模型，包括风场模型关键参数的标定，计算模拟点台风风速，进行台风风险灾害评估。

对于进入模拟点250km范围内的台风路径以及强度样本，采用Yan Meng风场模型计算模拟点的台风风速，具体步骤如下：

(1)台风风场模型关键参数的标定

采用Vickery经验模型估计最大风速半径R_max和Holland气压剖面参数B：

式中：Δp为中心气压差；ψ为纬度；f_c＝2×7.273×10^-5sinψ；

和ε_B为符合正态分布的零均值随机扰动项；

(2)建立气压模型

采用Holland气压模型：

p＝MSLP+Δp exp[-(R_max/r)^B] (7)

式中：p为距台风中心径向距离为r处的海平面压强；R_max为最大风速半径；B为气压剖面参数；

根据式(5)，气压梯度表达式为：

(3)建立平衡方程

中性大气条件下的Navier-Stokes方程为：

式中：V为空气微团的运动速度，Yan Meng风场模型把V视为梯度风速V_g和地表摩擦风速V′两部分矢量叠加，即V＝V_g+V′；f为Coriolis参数；F为边界层摩擦力；

由于径向风压梯度在台风边界层内随高度变化很小，再忽略边界层以上的摩擦力，式(9) 可以分解为梯度层和边界层两部分，分别如式(10)和式(11)：

结合

以及

式(10)和式(11)可进一步化简为：

V′·▽V′+V′·▽V_g+V·▽V′＝-fk×V′+F (13)

在随台风中心移动的二维极坐标中，沿径向和切向分解式(10)得到：

式中：V_θg和V_rg分别是切向和径向的梯度风速；c_r＝-V_T cos(θ-β′)，c_θ＝-V_T sin(θ-β′)；ρ是空气密度，取1.2kg/m³；θ是正东方向与模拟点和台风中心连线夹角，以逆时针旋转为正；β′是正东方向与台风移动方向的夹角，以逆时针旋转为正；

在上层大气环境中，考虑到V_rg远远小于V_θg，忽略式(14)中的前两项，可以求出V_θg；由于V_rg值很小，近似取0仍可满足风速模拟的准确性，因而，演化出梯度风速的解析式，如下式所示：

在近地边界层中，地表摩擦风速的切向和径向分量V′和V_r′小于相应的梯度风速分量，因此，这些分量对θ的一阶导数可假定同样小于相应的梯度风速分量的一阶导数；于是式(11) 可作线性化分解，如下式所示：

上层大气和地表附近的边界条件分别为式(18)和式(19)：

V′|_z′→∞＝0 (18)

式中，V_s表示近地面风速，计算起始点取为z＝h+10即z′＝0处，考虑边界条件后，式(17) 可以求解为：

V″＝Dexp[-(1+i)λz′] (20)

式中：边界层地表复常量D＝D₁+iD₂，引入式(21)所示的无量纲参数ξ和带量纲参数λ，最终可推导出地表摩擦风速的解析式，如式(22)所示：

D₁和D₂可按下式计算：

式中：C_d为阻力系数；k_m为运动黏度，取100m²/s；k为Karman常数，取0.4；平均粗糙单元高度

实测得A＝11.4；零平面位移d＝0.75h；z₁₀设定在平均粗糙单元h高 10m处；z₀为考虑地形和地表粗糙度影响而引入的等效粗糙长度；

(4)计算流程

Yan Meng台风风场模型计算地表风速的流程为：首先在二维极坐标内由式(16)计算梯度风速，然后取梯度风速作为地表风速的初始值，依次代入解析公式(23)、(21)和(22)求得地表摩擦风速，由梯度风速和地表摩擦风速叠加得到新的地表风速，经过多次迭代，直到收敛；

(5)台风极值风速预测及灾害风险评估

通过数值模拟得到台风极值风速序列后，需要用极值概率分布模型对该序列进行拟合，进而预测重现期极值风速；常用的极值概率分布有极值I型分布、极值II型分布和极值III 型分布；当模拟的台风样本足够多时，应当首选经验分布。

为了更加清晰地说明上述步骤，现以深圳地区的台风灾害风险评估为例来适当介绍，具体如下：

(1)建立起始点模型

图2所示为基于CMA热带气旋最佳路径数据集统计得到的1949年-2019年历史热带气旋年发生次数年际变化图。热带气旋在西北太平洋区域的年均发生次数为33.6，均方根6.99。由图2可以看出，在西北太平洋区域，热带气旋年频次在20世纪70年代前呈现一定上升趋势，而在70年代后下降趋势明显。

图3所示为采用负二项分布拟合得到的年发生次数概率分布图。基于拟合得到的分布模型参数，可以通过随机抽样获得模拟台风的年发生次数。

西北太平洋区域的历史热带气旋主要生成于北纬5度至30度之间，大约有17％的热带气旋生成于中国的南海。基于历史的起始点信息计算全路径模拟所需的信息：移动速度、朝向、所在位置的洋面温度、相对强度等，进而通过随机抽样可以确定每一条模拟台风的起始点信息。起始点信息的抽样样本需要剔除陆地起始点信息。

(2)建立行进模型、强度模型并进行终止判别

采用随机森林行进模型、强度模型模拟得到热带气旋下一点的经度、纬度和中心气压，并判定是否终止。针对不同的模型整理相应的训练数据集，确定样本总量及样本特征数，设置决策树数量和节点特征数即可进行训练。随机森林训练过程可基于MATLAB、python等平台提供的通用算法。随机森林由大量的单个决策树组成，训练时对数据进行有放回地抽样，并交由每棵“树”进行学习训练，最后集成所有“树”的预测得到最终结果。每棵“树”为一个弱学习机，只能考虑单方面的因素进行学习预测，但随机森林算法对大量单个“树”预测结果的集成分析能使最终结果在准确性上有大幅度提升。由于单个决策树为弱学习机，因此并不需要对其参数进行过多地调整以强化效能。随机森林主要设置参数为决策树的数量以及每个节点训练时采用的特征数。决策树的数量越多越精确但会影响计算效率，且在一定数量后准确性提升不明显。每次节点训练时抽取的特征数并无明确值，需要设置在一定范围内以满足相关性与分类能力的平衡，数量可以为总特征数的1/2、开方或开方的2倍等。

选取50棵树进行训练，特征抽取数为总特征的1/3。模拟了71年的热带气旋路径随机样本，与CMA的71年历史热带气旋路径(1949年-2019年)进行对比。从图4a-4b热带气旋密度图可以看出，整体上模拟气旋在计算区域内的密度与历史气旋的密度具有相近的热度分布规律。可以从图中定性地判断模拟路径与CMA历史路径在空间趋势上基本接近。

(3)路径模拟结果检验

将中国东南沿海海岸线划分成25个站点，相邻站点间距100km。模拟了10000年的热带气旋路径随机样本，图5a-5d所示为25个海岸站点250km范围内的CMA历史路径与模拟路径的关键参数对比。可以看出，模拟路径在海岸线各个站点范围内的台风关键参数与CMA的历史路径所对应的关键参数比较接近，可以验证热带气旋全路径模拟方法的有效性。为进一步验证效果，采用K-S检验方法对历史样本和模拟样本移动速度、朝向以及中心气压差三个关键参数的经验分布进行比较，如图6所示。发现除个别参数点外，其余参数K-S检验的P-value结果均大于0.05，表明历史参数经验分布和模拟参数经验分布具有较好的相似性。

(4)建立Yan Meng台风风场模型

假定当某条热带气旋中心距模拟点的最短距离小于250km时，模拟点受该条热带气旋影响，进而激活风场模型用于风险灾害评估。图7为采用Yan Meng台风风场模型计算得到的深圳地区台风年极值风速序列，图8为对应的概率分布图，由于台风极值风速样本足够多，所以采用经验分布。表1给出了10/50/100年重现期下深圳地区的台风极值风速，并与中国荷载规范对应的设计风速进行了对比。由表1可以看出，在10年和50年重现期下，根据荷载规范给出的设计风速偏于风险。

表1

上述具体实施方式用来解释说明本发明，仅为本发明的优选实施例，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改、等同替换、改进等，都落入本发明的保护范围。

Claims (5)

1.一种基于随机森林算法的热带气旋全路径模拟方法，其特征在于：包括如下步骤：

第一步，建立起始点模型，包括抽样生成模拟热带气旋的年发生次数和生成起始点信息；

第二步，建立行进模型，模拟热带气旋的经纬度形成热带气旋路径；

需要事先对路径进行聚类分析以将其分类，第二步中采用基于机器学习的随机森林路径回归模型以及第四步的随机森林路径终止判别均针对不同路径类型分别学习建立；按热带气旋的分布与行进特征分成四类，分别对各类热带气旋的数据进行回归训练；类型1：从生成至终止均在洋面且主要在东西方向运行分布；类型2：主要行进路线为西进方向，可能途径菲律宾、中国 南海、中国东南沿海诸省和东南亚部分区域；类型3：初始行进路线为西北方向，之后向东北方向行进可能途径中国东部沿海诸省、日本和朝鲜半岛；类型4：从生成至终止均在洋面且与类型1相比，热带气旋更多地向北发展；

第二步中采用基于机器学习的随机森林回归模型预测下一点的经度和纬度；基于随机森林回归算法分别对各类路径数据的经纬度信息特点进行学习训练，生成路径预测回归模型，据此模型基于前三步的经纬度信息对下一步经纬度做出预测：

式中：RFr_lat和RFr_lon分别为用于纬度预测和经度预测的随机森林回归模型；j＝1,2,…,n，表示采用的路径类别，n表示路径类别总数；lat为纬度，lon为经度；i+1、i、i-1、i-2、i-3分别表示下一、当前、前一、前二、前三时间步；ε_lat和ε_lon为符合正态分布的零均值随机扰动项；

第三步，建立强度模型，包括洋面强度发展模型和陆地强度衰减模型，并根据热带气旋路径中下一点的经纬度选择相应的强度模型；

第三步中洋面强度发展模型是基于随机森林回归算法对所有路径数据的位置、强度、洋面温度进行学习训练生成的，以预测洋面上下一模拟点的相对强度，从而确定中心气压：

I_i+1＝RFr_I(lat_i,lon_i,I_i,I_i-1,I_i-2,SST_i+1,SST_i,SST_i-1,SST_i-2,MSLP_i,MSLP_i-1,MSLP_i-2)+ε_I(2)

式中：RFr_I为用于热带气旋强度预测的随机森林回归模型；I为相对强度；SST为洋面温度，单位为K；ε_I为符合正态分布的零均值随机扰动项；MSLP为最小海平面气压，单位为hpa；综合考虑洋面温度、平流层温度、相对湿度影响，将中心气压转换为相对强度；

第三步中陆地强度衰减模型是综合考虑登陆持续时间、热带气旋强度大小以及地表粗糙度影响，建立登陆强度回归模型，模拟热带气旋登陆后的强度衰减：

Δp_t＝Δp₀exp(-mt) (3)

式中：Δp₀为登陆时的中心气压差；Δp_t为登陆t小时之后的中心气压差；m＝m₀+m₁Δp₀+ε_m，系数m₀、m₁和随机扰动项ε_m可通过历史登陆热带气旋的气压差变化回归分析得到；

第四步，终止判别，基于生成的气旋强度信息、位置信息判决该条模拟气旋是否终止；

第四步中采用随机森林分类算法对热带气旋的强度、经纬度、洋面温度、行进时间进行学习训练，建立针对不同路径类别的热带气旋进行终止判别：

end_i+1＝RFc_end,j(lat_i,lon_i,SST_i,MSLP_i,N_i) (4)

式中：end表示热带气旋是否终止的判决结果，属于逻辑标量，“1”为终止，“0”，为不终止；RFc_end表示用于终止判定的随机森林分类模型；N表示热带气旋已行进的时间步；

第五步，重复上述步骤直至生成给定模拟年份的所有路径，并检验路径模拟结果，对随机生成的热带气旋路径及强度大样本数据的关键参数进行统计分析；

第六步，基于模拟生成的气旋路径及强度大样本数据建立台风风场工程模型和边界层模型，包括风场模型关键参数的标定，计算模拟点台风风速，以进行极值分布分析。

2.根据权利要求1所述的基于随机森林算法的热带气旋全路径模拟方法，其特征在于：第一步中采用负二项分布对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数；基于历史的热带气旋起始点信息计算全路径模拟所需的信息：移动速度、朝向、所在位置的洋面温度、相对强度，进而通过随机抽样确定每一条模拟台风的起始点信息。

3.根据权利要求1所述的基于随机森林算法的热带气旋全路径模拟方法，其特征在于：第三步中对热带气旋路径中下一点的经纬度进行判断，如果下一点是位于陆地，则选择陆地强度衰减模型，如果下一点是位于洋面，则选择洋面强度发展模型。

4.根据权利要求1所述的基于随机森林算法的热带气旋全路径模拟方法，其特征在于：第五步中模拟结果检验是对比随机模拟的路径及强度与历史热带气旋数据之间的差异，定量刻画模拟结果的精度；第五步中的统计规律分析是将中国东南沿海海岸线划分为若干海岸站点，再统计每一个海岸点250km范围内模拟的和历史的在热带气旋年均频次、移动速度、朝向以及中心气压差的结果差异，以反映热带气旋全路径模拟的空间分布效果；采用Kolmogorov–Smirnov test对历史样本和模拟样本移动速度、朝向以及中心气压差三个关键参数的经验分布进行比较。

5.根据权利要求1～4中任意一项所述的基于随机森林算法的热带气旋全路径模拟方法，其特征在于：第六步中对于进入模拟点250km范围内的台风路径以及强度样本，采用YanMeng风场模型计算模拟点的台风风速，具体步骤如下：

(1)台风风场模型关键参数的标定

采用Vickery经验模型估计最大风速半径R_max和Holland气压剖面参数B：

式中：Δp为中心气压差；ψ为纬度；f_c＝2×7.273×10^-5sinψ；

和ε_B为符合正态分布的零均值随机扰动项；

(2)建立气压模型

采用Holland气压模型：

p＝MSLP+Δpexp[-(R_max/r)^B] (7)

式中：p为距台风中心径向距离为r处的海平面压强；R_max为最大风速半径；B为气压剖面参数；

根据式(5)，气压梯度表达式为：

(3)建立平衡方程

中性大气条件下的Navier-Stokes方程为：

式中：V为空气微团的运动速度，Yan Meng风场模型把V视为梯度风速V_g和地表摩擦风速V′两部分矢量叠加，即V＝V_g+V′；f为Coriolis参数；F为边界层摩擦力；

由于径向风压梯度在台风边界层内随高度变化很小，再忽略边界层以上的摩擦力，式(9)可以分解为梯度层和边界层两部分，分别如式(10)和式(11)：

结合

以及

式(10)和式(11)可进一步化简为：

在随台风中心移动的二维极坐标中，沿径向和切向分解式(10)得到：

式中：V_θg和V_rg分别是切向和径向的梯度风速；c_r＝-V_Tcos(θ-β′)，c_θ＝-V_Tsin(θ-β′)；ρ是空气密度，取1.2kg/m³；θ是正东方向与模拟点和台风中心连线夹角，以逆时针旋转为正；β′是正东方向与台风移动方向的夹角，以逆时针旋转为正；

在上层大气环境中，考虑到V_rg远远小于V_θg，忽略式(14)中的前两项，可以求出V_θg；由于V_rg值很小，近似取0仍可满足风速模拟的准确性，因而，演化出梯度风速的解析式，如下式所示：

在近地边界层中，地表摩擦风速的切向和径向分量V′和V_r′小于相应的梯度风速分量，因此，这些分量对θ的一阶导数可假定同样小于相应的梯度风速分量的一阶导数；于是式(11)可作线性化分解，如下式所示：

上层大气和地表附近的边界条件分别为式(18)和式(19)：

V′|_z′→∞＝0 (18)

式中，V_s表示近地面风速，计算起始点取为z＝h+10即z′＝0处，考虑边界条件后，式(17)可以求解为：

V″＝Dexp[-(1+i)λz′] (20)

式中：边界层地表复常量D＝D₁+iD₂，引入式(21)所示的无量纲参数ξ和带量纲参数λ，最终可推导出地表摩擦风速的解析式，如式(22)所示：

D₁和D₂可按下式计算：

式中：C_d为阻力系数；k_m为运动黏度，取100m²/s；k为Karman常数，取0.4；平均粗糙单元高度

实测得A＝11.4；零平面位移d＝0.75h；z₁₀设定在平均粗糙单元h高10m处；z₀为考虑地形和地表粗糙度影响而引入的等效粗糙长度；

(4)计算流程

Yan Meng台风风场模型计算地表风速的流程为：首先在二维极坐标内由式(16)计算梯度风速，然后取梯度风速作为地表风速的初始值，依次代入解析公式(23)、(21)和(22)求得地表摩擦风速，由梯度风速和地表摩擦风速叠加得到新的地表风速，经过多次迭代，直到收敛；

(5)台风极值风速预测及灾害风险评估

通过数值模拟得到台风极值风速序列后，需要用极值概率分布模型对该序列进行拟合，进而预测重现期极值风速；常用的极值概率分布有极值I型分布、极值II型分布和极值III型分布；当模拟的台风样本足够多时，应当首选经验分布。

Images