一种面向灾害风险评估的热带气旋全路径模拟方法与流程

本发明属于工程结构抗风性能设计领域，尤其涉及一种面向灾害风险评估的热带气旋全路径模拟方法。

背景技术：

台风是热带气旋的一个类别。在气象学上，按世界气象组织定义：热带气旋中心持续风速达到12级(即每秒32.7米或以上)称为飓风(hurricane)，飓风的名称使用在北大西洋及东太平洋；而西北太平洋使用的近义字是台风(typhoon)。在中国台湾、日本等地，则将中心持续风速每秒17.2米以上的热带气旋皆称为台风。一个发展成熟的台风，按其结构和带来的天气，分为台风眼、涡旋风雨区、外围大风区三部分，从中心向外呈同心圆状排列。台风眼位于台风中心，直径约5～10公里。台风眼内盛行下沉气流，故天气睛朗，风平浪静。台风眼外侧为涡旋风雨区，这里盛行强烈的辐合上升气流，形成浓厚的云层，出现狂风暴雨，风力常常在12级以上，是台风中天气最恶劣的区域。再向外为外围大风区，风速向外减小，风力通常在6级以上。台风过境常常带来狂风暴雨天气，引起海面巨浪，严重威胁沿海工程结构安全。

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

在工程设计领域，面向灾害风险评估的热带气旋数值模拟方法主要为基于单个站点或是局部区域范围的过境热带气旋实测数据，应用台风风场模型通过数值计算模拟台风风速，具体步骤如下：(1)采用模拟圆法提取一定范围内的过境台风记录，并进行数据预处理；(2)统计得到模拟点处各个台风关键参数的概率模型；(3)依据关键参数的概率模型进行随机抽样，得到一组台风关键参数，结合台风风场模型，计算模拟点的台风风速；(4)重复上步，得到模拟点的台风最大风速序列，用于台风灾害风险评估(推算不同重现期对应的极值风速或确定台风平均风剖面参数等)。上述基于单个站点或是局部区域范围的台风数值模拟方法操作步骤少，计算效率高，但同时也存在着局限性，即需要统计得到模拟点处各个台风关键参数的概率模型。但是，中国东南沿海有些地区比如浙江省的热带气旋事件总体样本非常少，采用模拟圆法筛选得到的过境台风记录不足以建立该区域内的台风关键参数的概率模型。

技术实现要素：

本发明提出了一种能实现精细化评估、结果精确可靠、计算效率高、适用范围广的面向灾害风险评估的热带气旋全路径模拟方法。

本发明采用的技术方案是：

一种面向灾害风险评估的热带气旋全路径模拟方法，其具体步骤包括：

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

第二步，建立行进模型，模拟热带气旋的移动速度和朝向形成热带气旋路径；

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

第四步，路径模拟结果检验，对随机生成的热带气旋路径及强度大样本进行统计结果检验以及统计规律分析；

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

进一步，第一步中采用负二项分布对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。根据西北太平洋区域的历史热带气旋年发生次数数据特点，对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。目前主要有负二项分布和泊松分布两种模型来拟合热带气旋的年发生次数概率分布。georgiou对墨西哥湾和大西洋海岸遭遇的台风进行了研究，比较了负二项分布和泊松分布对台风年发生率模型的拟合程度，发现负二项分布能够更好地模拟海面上的台风序列的发生率，而泊松分布能更好地模拟沿岸发生的台风序列的发生率，这是因为台风的移动路径受气候参数影响较大，而台风的发生不受这些因素的影响。因此，在全路径模拟中，起源于洋面上的台风年发生次数宜用负二项分布模拟。

进一步，第一步中基于历史的热带气旋起始点信息计算全路径模拟所需的信息：移动速度、朝向、所在位置的洋面温度、相对强度，进而通过随机抽样确定每一条模拟台风的起始点信息。

进一步，第二步中采用vickery路径回归模型预测下一点的移动速度和朝向：

式中v为移动速度，△lnv＝lnvi+1-lnvi；θ为朝向，△θ＝θi+1-θi；ψ纬度；λ经度；εv和εθ为符合正态分布的零均值随机扰动项；回归系数ai，bi只跟地理区域有关，将西北太平洋区域划分成若干网格，收集每个网格内的历史热带气旋路径信息，并区分向东朝向和向西朝向，系数ai，bi可通过每个网格的历史热带气旋信息回归分析得到。

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

进一步，第三步中洋面强度发展模型是根据洋面强度回归模型，预测洋面上下一点的相对强度，从而确定中心气压：

式中i为相对强度；t为洋面温度；εi为符合正态分布的零均值随机扰动项，系数ci可以通过每个网格的历史热带气旋信息回归分析得到。本发明综合考虑洋面温度、平流层温度、相对湿度等环境因素影响，将中心气压转换为相对强度。

进一步，第三步中陆地强度衰减模型利用登陆强度回归模型，模拟热带气旋登陆后的强度衰减：

△p(t)＝△p0exp(-mt)(3)

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

进一步，第四步中模拟结果检验是对比随机模拟的路径及强度与历史热带气旋数据之间的差异，定量刻画模拟结果的精度。

进一步，第四步中的统计规律分析是将中国东南沿海海岸线划分为若干海岸站点，再统计每一个海岸点250km范围内模拟的和历史的在热带气旋年均频次、移动速度、朝向以及中心气压差的结果差异，从一定程度上可以反映热带气旋全路径模拟的空间分布效果。

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

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

采用vickery经验模型估计最大风速半径rmax和holland气压剖面参数b：

式中△p为中心气压差；ψ为纬度；fc＝2×7.273×10^-5sinψ；和εb为符合正态分布的零均值随机扰动项；

(2)建立气压模型

采用holland气压模型：

p＝p0+△pexp[-(rmax/r)^b](6)

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

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

(3)建立平衡方程

中性大气条件下的navier-stokes方程为：

式中，v为空气微团的运动速度，yanmeng风场模型把v视为梯度风速vg和地表摩擦风速v’两部分矢量叠加，即v＝vg+v’；f为coriolis参数；f为边界层摩擦力；

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

结合以及式(9)和式(10)可进一步化简为：

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

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

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

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

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

v′|z′→∞＝0(17)

计算起始点取为z＝h+10即z’＝0处，考虑边界条件后，式(16)可以求解为：

v″＝dexp[-(1+i)λz′](19)

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

d1和d2可按下式计算：

式中，cd为阻力系数；km为运动黏度，取100m²/s；k为karman常数，取0.4；平均粗糙单元高度h＝az0^0.86，实测得a＝11.4；零平面位移d＝0.75h；z10设定在平均粗糙单元h高10m处；z0为考虑地形和地表粗糙度影响而引入的“等效粗糙长度”；

(4)计算流程

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

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

通过数值模拟得到台风极值风速序列后，需要用极值概率分布模型对该序列进行拟合，进而预测重现期极值风速。常用的极值概率分布有极值i型分布(即gumbel分布)、极值ii型分布(即fréchet分布)和极值iii型分布。当模拟的台风样本足够多时，应当首选经验分布。

进一步，第五步之前要先对热带气旋路径中下一点的经纬度进行判断，当模拟点距热带气旋中心最短距离小于250km时进行第五步，当模拟点距热带气旋中心最短距离大于250km时判断热带气旋行进路径是否满足终止条件，满足终止条件则结束，不满足则继续第二步。

有益效果：本发明所提供的一种基于整个西北太平洋海域的热带气旋全路径模拟方法能够克服局部区域历史热带气旋样本不足的问题，对整个西北太平洋区域的热带气旋生成、行进路径、行进方向、洋面强度发展以及登陆强度衰减进行随机模拟，生成大量符合历史样本特征的热带气旋路径及强度随机事件样本集，据此进行中国东南沿海地区台风灾害风险的定量精细化评估，必将获得更加精确可靠的结果，从而有效地解决局部区域历史样本不足的问题，为灾害风险的定量精细化评估奠定了基础。该方法计算效率高，可推广到其他区域进行飓风灾害风险评估分析。

附图说明

图1为本发明的技术流程图。

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

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

图4为历史热带气旋起始点空间分布图。

图5a为行进模型回归系数a1空间分布图：(左)向西朝向；(右)向东朝向。

图5b为行进模型回归系数a2空间分布图：(左)向西朝向；(右)向东朝向。

图5c为行进模型回归系数a3空间分布图：(左)向西朝向；(右)向东朝向。

图5d为行进模型回归系数εv空间分布图：(左)向西朝向；(右)向东朝向。

图6a为影响浙江省近海100km范围内的cma历史路径与模拟路径比较。

图6b为影响福建省近海100km范围内的cma历史路径与模拟路径比较。

图6c为影响广东省近海100km范围内的cma历史路径与模拟路径比较。

图7为海岸站点的空间分布图。

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

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

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

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

图9为热带气旋全路径模拟中激活风场模型的示意图。

图10为深圳地区台风年极值风速序列。

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

具体实施方式

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

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

第一步，建立起始点模型，包括抽样生成模拟热带气旋的年发生次数和生成起始点信息；具体的采用负二项分布对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。根据西北太平洋区域的历史热带气旋年发生次数数据特点，对年发生次数进行模型参数拟合，并抽样生成模拟热带气旋的年发生次数。目前主要有负二项分布和泊松分布两种模型来拟合热带气旋的年发生次数概率分布。georgiou对墨西哥湾和大西洋海岸遭遇的台风进行了研究，比较了负二项分布和泊松分布对台风年发生率模型的拟合程度，发现负二项分布能够更好地模拟海面上的台风序列的发生率，而泊松分布能更好地模拟沿岸发生的台风序列的发生率，这是因为台风的移动路径受气候参数影响较大，而台风的发生不受这些因素的影响。因此，在全路径模拟中，起源于洋面上的台风年发生次数宜用负二项分布模拟。

基于历史的热带气旋起始点信息计算全路径模拟所需的信息：移动速度、朝向、所在位置的洋面温度、相对强度，进而通过随机抽样确定每一条模拟台风的起始点信息。

第二步，建立行进模型，模拟热带气旋的移动速度和朝向形成热带气旋路径；采用vickery路径回归模型预测下一点的移动速度和朝向：

式中v为移动速度，△lnv＝lnvi+1-lnvi；θ为朝向，△θ＝θi+1-θi；ψ纬度；λ经度；εv和εθ为符合正态分布的零均值随机扰动项；回归系数ai，bi只跟地理区域有关，将西北太平洋区域划分成若干网格，收集每个网格内的历史热带气旋路径信息，并区分向东朝向和向西朝向，系数ai，bi可通过每个网格的历史热带气旋信息回归分析得到。

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

综合考虑洋面温度、平流层温度、相对湿度等环境因素影响，将中心气压转换为相对强度。洋面强度发展模型根据洋面强度回归模型，预测洋面上下一点的相对强度，从而确定中心气压：

式中i为相对强度；t为洋面温度；εi为符合正态分布的零均值随机扰动项，系数ci可以通过每个网格的历史热带气旋信息回归分析得到。

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

△p(t)＝△p0exp(-mt)(3)

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

第四步，路径模拟结果检验，对随机生成的热带气旋路径及强度大样本进行统计结果检验以及统计规律分析；具体的模拟结果检验是对比随机模拟的路径及强度与历史热带气旋数据之间的差异，定量刻画模拟结果的精度。

统计规律分析是将中国东南沿海海岸线划分为若干海岸站点，再统计每一个海岸点250km范围内模拟的和历史的在热带气旋年均频次、移动速度、朝向以及中心气压差的结果差异，从一定程度上可以反映热带气旋全路径模拟的空间分布效果。

第五步，建立台风风场工程模型和边界层模型，包括风场模型关键参数的标定，计算模拟点台风风速，进行台风风险灾害评估。在这之前要先对热带气旋路径中下一点的经纬度进行判断，当模拟点距热带气旋中心最短距离小于250km时进行第五步，当模拟点距热带气旋中心最短距离大于250km时判断热带气旋行进路径是否满足终止条件，满足终止条件则结束，不满足则继续第二步。

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

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

采用vickery经验模型估计最大风速半径rmax和holland气压剖面参数b：

式中△p为中心气压差；ψ为纬度；fc＝2×7.273×10^-5sinψ；和εb为符合正态分布的零均值随机扰动项；

(2)建立气压模型

采用holland气压模型：

p＝p0+△pexp[-(rmax/r)^b](6)

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

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

(3)建立平衡方程

中性大气条件下的navier-stokes方程为：

式中，v为空气微团的运动速度，yanmeng风场模型把v视为梯度风速vg和地表摩擦风速v’两部分矢量叠加，即v＝vg+v’；f为coriolis参数；f为边界层摩擦力；

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

结合以及式(9)和式(10)可进一步化简为：

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

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

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

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

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

v′|z′→∞＝0(17)

计算起始点取为z＝h+10即z’＝0处，考虑边界条件后，式(16)可以求解为：

v″＝dexp[-(1+i)λz′](19)

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

d1和d2可按下式计算：

式中，cd为阻力系数；km为运动黏度，取100m²/s；k为karman常数，取0.4；平均粗糙单元高度h＝az0^0.86，实测得a＝11.4；零平面位移d＝0.75h；z10设定在平均粗糙单元h高10m处；z0为考虑地形和地表粗糙度影响而引入的“等效粗糙长度”；

(4)计算流程

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

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

通过数值模拟得到台风极值风速序列后，需要用极值概率分布模型对该序列进行拟合，进而预测重现期极值风速。常用的极值概率分布有极值i型分布(即gumbel分布)、极值ii型分布(即fréchet分布)和极值iii型分布。当模拟的台风样本足够多时，应当首选经验分布。

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

(1)建立起始点模型

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

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

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

(2)建立行进和强度模型

采用vickery的行进模型和强度模型模拟热带气旋下一点的移动速度、朝向和中心气压，模型回归系数ai，bi，ci只跟地理区域有关。将西北太平洋区域划分成2度*2度网格，其中在105e-130e,15n-35n区域划分成1度*1度网格。基于cma热带气旋最佳路径数据集收集每个网格内的历史热带气旋路径信息，并区分向东朝向和向西朝向，系数ai，bi，ci可以通过每个网格内的历史热带气旋信息回归分析得到。图5(a-d)分别为模型回归系数a1、a2、a3、εv的空间分布图。

模拟了67年的热带气旋路径随机样本，与cma的67年历史热带气旋路径(1949年-2015年)进行对比，图6(a-c)分别为影响浙江省、福建省、广东省近海100km范围内的cma历史路径与模拟路径空间趋势比较图。可以从图中定性地判断模拟路径与cma历史路径在空间趋势上基本接近。随着纬度增加，热带气旋的朝向从向西朝向逐渐转变为向东朝向。

(3)路径模拟结果检验

将中国东南沿海海岸线划分成25个站点，相邻站点间距100km，图7为海岸站点编号s1-s25的空间分布图。模拟了100000年的热带气旋路径随机样本，图8(a-d)所示为25个海岸站点250km范围内的cma历史路径与模拟路径的关键参数对比。可以看出，模拟路径在海岸线各个站点范围内的台风关键参数与cma的历史路径所对应的关键参数基本吻合，可以验证热带气旋全路径模拟方法的有效性。

(4)建立yanmeng台风风场模型

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

表1深圳地区不同重现期下台风极值风速