摘要
为了校正已有的确定性台风强度模型并考虑随机性的影响，在确定性台风强度的常微分方程中引入了包括均值和随机噪声的修正项，并采用多种偏态分布模型作为随机噪声的候选概率分布。利用西北太平洋地区的历史台风数据，采用地理加权方法估计了修正项的地理变化均值、标准差、偏度和超值峰度。进一步采用矩估计方法识别了修正项中随机噪声候选概率分布模型的参数，并基于KS距离确定了最优概率分布模型。通过对比分析历史台风强度演化过程的模拟结果，研究了误差项的均值部分和随机部分对模型性能的影响。结果表明，修正项的引入显著改进了模型对历史台风强度的模拟能力，增强了模型对台风强度随机性的反映能力,验证了该模型在台风极值风速分析中的有效性。
Abstract
To adjust existing deterministic typhoon intensity models and account for stochastic influences, a correction term containing both mean and random noise is introduced into the ordinary differential equations governing deterministic typhoon intensity. Various skewed distribution models serve as candidate probability distributions for the random noise. Using historical typhoon data from the Northwest Pacific, the geographically weighted method estimates geographic variation in mean, standard deviation, skewness, and excess kurtosis of the correction term. Furthermore, the method of moments estimates parameters for candidate probability distribution models of the random noise within the correction term, with the optimal probability distribution model determined by KS distance. By comparing simulated results of historical typhoon intensity evolution, the impact of both mean and stochastic components of the error term on model performance is examined. The results indicate that introducing the correction term significantly improves the model′s ability to simulate historical typhoon intensities and enhances its capacity to capture the stochastic nature of typhoon intensity. Additionally, this paper validates model effectiveness in extreme wind speed analysis of typhoons.
Keywords
中国东南沿海地区常年受东亚季风和西北太平洋台风的影响，每年台风登陆都会给沿海地区带来严重的生命财产损失，其中，大部分与台风地表强风有关。因此，准确评估台风影响地区的台风极值风速对保障工程结构抗风安全具有重要意义[1-2]。由于缺乏可靠的长期台风地表风速观测数据，通常采用模拟方法人工生成台风样本，并进行极值风速统计分析[3]。其中，台风强度模型是台风风速模拟的关键之一。
Vickery等[4]基于自回归随机过程模型较早地提出了台风强度模拟方法，并进行了北大西洋地区的台风极值风速分析。Li等[5]的研究指出Vickery等的模型可能存在过拟合问题，并优化了台风强度的自回归模型。尽管自回归模型在描述台风强度演变的随机特征方面取得了重要进展，但自回归模型本质属于统计模型，难以准确反映环境变量与台风强度演化之间的复杂物理联系。事实上大量研究指出，为了保证沿海地区台风统计特性与实测结果的一致性，基于自回归模型的台风强度模拟研究或多或少需要对洋面温度等物理量进行经验调节[6-7]，从而增加了建模复杂性和主观因素的影响。因此，更合理的方式是基于物理机制建立台风风速全过程模拟方法。
事实上，已有相关研究利用数值天气预报模型分析台风全过程演化及其对工程结构的影响[8]。研究结果表明，数值天气预报模型能够综合考虑全物理要素，具有较高的精度。然而，台风危险性分析涉及到大量确定性样本的数值求解，数值天气预报模型较低的分析效率难以适应这一要求。为此，学者们从台风物理背景和基本控制方程出发，深入研究台风演化的物理机制，逐步发展了物理驱动的台风强度演化高效模拟方法[9-11]。事实上，台风本质上是一个作用在洋面（热源）和高空大气（冷凝器）之间的热机，由此决定了台风可能达到的最大潜在强度[12]。基于此，文献[13]综合天气尺度气流竖向切变、大气海洋耦合作用等因素的影响，提出了针对北大西洋地区的台风强度演化高效动力学模型。文献[14]进一步研究了该模型在西北太平洋地区的适用性。值得指出的是，尽管上述台风强度动力学模型建立了台风强度演化与基本物理量之间的确定性联系，但尚未考虑台风中显著的不确定性，从而无法精确评估台风极值风速。为此，本研究基于台风潜在强度理论，在台风强度动力学模型中引入考虑偏态分布特征的随机噪声，合理描述台风强度演化中的不确定性，进而发展随机台风强度模型。
1 台风随机强度模型
1.1 台风潜在强度理论
如图1所示，台风是一个作用在洋面（热源）和对流层顶（冷凝器）之间的热机：海洋表面的热量和水汽通过蒸发传输到大气中，较高的海表温度导致更多的潜热释放，为台风提供能量来源；海洋表面传输到边界层的热量和湿度在对流层内引发强烈的上升气流，导致大规模的云团和降水形成；这些对流活动将潜热释放到高层大气中，进一步驱动台风的增强。基于此，Emanuel提出了给定环境下台风中心地表最大风速可能达到的最大潜在强度，即
(1)
式中：Vp为台风的潜在强度，m/s；ζ为台风热机的效率，ζ=（Ts-To）/Ts; Ts为海表温度，K；To为对流层顶温度，K；k*0为饱和静态能；k为边界层中的湿静能。
图1台风Carnot热机循环系统
Fig.1Typhoon Carnot heat engine circulation system
由于竖向切变、大气-海洋相互作用下洋面温度降低以及中心气压的下降导致饱和比湿度的增加等因素的影响，实际台风中心最大风速无法达到Vp且处于动态变化之中。为此，Emanuel发展了动力学模型描述台风强度的演化，即
(2)
(3)
式中：Cd为表面阻力系数; h为边界层高度; m为湿度变量; S为25~85 kPa的垂直风切变; a1、a2为模型参数; γ=ζ+ακ，α为考虑台风吸热后洋面温度下降效应的因子，κ为中心气压的下降导致饱和比湿度增加的因子。α和κ可由下式计算：
(4)
(5)
式中：Γ为洋面混合层下的温度递减率，K/100 m； hm为海洋混合层深度，m；ut为台风移动速度，m/s；Lv为水汽潜热，2.25×106 J/kg；q*0为表面饱和比湿度；Rd为气体常数，287 J/（kg·K）。
综上，给定洋面温度、对流层顶温度、海洋混合层深度、洋面混合层下的温度递减率、25~85 kPa的垂直风切变，可以根据式（1）~（5）模拟得到确定性台风强度演化过程。
1.2 台风强度随机动力学模型
尽管上述台风潜在强度模型可以模拟台风强度演化过程，但一方面其无法反映台风中显著随机性，因此，难以精确评估台风危险性；另一方面，作者的前期工作表明，由于原模型是基于北大西洋飓风数据提出的，在模拟西北太平洋地区历史台风强度中存在一定的系统性偏差。鉴于此，在台风强度动力学方程基础上引入修正项，用以考虑随机性并校正模型：
(6)
式中：ε为本文引用的修正项，可以将其假设为均值部分和随机噪声的总和，即
(7)
式中：μ（λ，φ）和σ（λ，φ）分别为随经纬度变化的残差均值和残差标准差，ξ为零均值、标准差为1的随机噪声。
采用显式欧拉方法计算式（6）中的随机微分方程。首先，根据上述的观测值计算出Vp、α和γ，然后利用观测到的台风强度和垂直风切变对式（3）进行积分，其中，湿度变量的初始值由式（6）求解，使用前6 h的观测值，ε=0。推导出湿度变量后，用台风强度变化率减去式（6）右侧第一项的值即可计算出ε。
为了反映随机噪声的偏态分布特征，可以采用概率分布模型对其进行建模。本文前期工作表明，随机噪声ξ表现出显著的偏态特征，因此，台风强度随机源概率建模中，除了采用经典的正态分布模型，还需要考虑偏态概率分布模型。文献[15]的工作表明，偏态正态分布、Johnson SU和Johnson SB分布能够较好地描述随机变量的偏态特征。因此，采用如表1所示的4类概率分布模型作为随机噪声ξ的候选概率分布。
表1本研究采用的概率分布模型
Tab.1Probability distribution model used in this study
2 随机源概率模型
2.1 基于地理加权方法的统计矩估计
为确定基本随机源ε概率分布模型，采用矩估计方法分别识别了表1中ξ候选概率分布的参数，并基于KS距离准则选择了最优概率分布模型[16]。为此，需要计算随机残差实测前四阶统计矩及其随经纬度的变化。
本研究采用了中国气象局台风最佳路径数据[17]、NOAA 20世纪再分析气象数据[18]，针对1980—2013年西北太平洋地区发生的1 010个台风，结合式（1）~（6）计算ε的实测值。其中，采用中国气象局台风最佳路径数据的台风地表中心最大风速估计式（6）中的台风强度；采用NOAA 20世纪再分析气象数据估计式（6）中洋面温度、对流层顶温度、海洋混合层深度、洋面混合层下的温度递减率、25~85 kPa的垂直风切变。
由于台风运动的随机性，历史台风中心的位置在空间中是随机散布的。任意空间位置一般不存在ε的多个历史观测记录。为此，采用地理加权方法[19-20]计算平均值、标准差、偏度和超值峰度随纬度和经度变化的函数。设{x1，x2，···，xN}为观测样本，其中，第i个样本xi观测位置位于（λi，φi）。则局部前四阶统计矩可以采用下述公式计算。
均值
(8)
标准差
(9)
偏度
(10)
超值峰度
(11)
式中wi（λ，φ）为第i个样本的权重，即
(12)
其中，di为样本点（λi，φi）到（λ，φ）的距离，h为带宽，h的作用是通过权重系数筛选目标区域附近的台风样本，从而计算局地的统计特征，因此，带宽取值不宜过大；同时，为保证统计特征计算结果的可靠性，需要有足够的台风样本，带宽取值也不宜过小。通过试算发现，采用h=0.4°能够较好地兼顾上述两个准则。
图2为西北太平洋地区的均值、标准差、偏度和超值峰度等高线。可以看出：1）修正项均值和偏度在区域东南部普遍为负，在西北部普遍为正，这说明原模型可能高估了西南偏东南地区的台风强度; 2）标准差的变化规律与海岸线形状一致，表明地表条件是决定噪声强度的重要因素; 3）超值峰度在区域西北部的值较大，这意味着数据在均值附近集中，但极值出现的概率较高。
图2修正项前四阶矩的估计值
Fig.2Estimated value of the first four moments of the correction term
2.2 概率分布拟合与优选
以经纬度分别为145.5°E、12.6°N的位置为例，由式（8）~（12）计算的该地区ε均值、标准差、偏度和超值峰度分别为-0.5 m/s2、2.6 m/s2、0.1、2.5。根据矩估计，候选概率分布模型参数如表1所示；KS距离分别为0.098 9、0.103 8、0.037 4、0.087 4；各模型的拟合结果如图3所示。可以看出，Johnson SU概率分布模型取得KS距离最小值，因此，取其为该点ε的概率模型。进一步，采用上述方法确定了西北太平洋地区其他位置台风强度随机源ε的概率分布模型，结果如图4所示。可以看出，在中国东南沿海的绝大部分区域，Johnson SU分布为最优的概率分布。
图3（145.5°E、12.6°N）处修正项的概率密度函数（PDF）和累积分布函数（CDF）
Fig.3Probability density function (PDF) and cumulative distribution function (CDF) of the correction term at (145.5°E, 12.6°N) position
图4不同经纬度位置的最优概率分布
Fig.4Optimal probability distributions for locations at different latitudes and longitudes
3 模型对比与应用
3.1 历史台风强度模拟
为了验证随机强度模型的准确性及研究残差项各部分对台风强度模拟的影响，考虑以下3种模型：
1）模型1。强度方程中不考虑修正项，即式（7）中的ε=0。
2）模型2。强度方程中只考虑修正项的均值部分，即式（7）中ε=μ（λ，φ）。
3）模型3。在强度方程中考虑了全部修正项，即式（7）中的ε=μ（λ，φ）+σ（λ，φ）ξ（t）。
比较模型1和模型2可以研究修正项均值部分的影响，比较模型2和模型3可以研究修正项随机部分的影响。
用上述3种模型分别模拟了西北太平洋地区1980—2013年1 010个历史台风的强度演化过程。图5给出了4个典型台风的模拟与实测结果对比。为了表征模型3的随机性，给出了模型3模拟结果的中位值、30%中心区域、50%中心区域、70%中心区域和90%中心区域。由图5可以看出，模型1和模型2之间存在明显的差异，而模型2和模型3中位值之间的差异不显著。如图5（b）所示，2007年9月12日—9月18日，模型1模拟台风Pat的强度从15 m/s增加到61 m/s；模型2和模型3分别增加到57 m/s和50 m/s。结果表明，与模型1相比，模型2和模型3的模拟的结果与观测值更为一致。
图5单个样本的模拟和观测风速演变对比
Fig.5Comparison of simulated and observed wind speed evolution for a single sample
为了验证模型对中国东南沿海地区台风强度概率特征的模拟能力，考虑了如图6所示的6个沿海城市局地台风强度的概率分布。其中，局地台风强度定义为台风中心首次进入目标位置中心附近250 km圆形区域时的强度。在大部分地区，模型1的CDF曲线位于模型2的左侧，表明模型1预测的台风强度小于模型2。在上海、宁波、阳江和湛江，模型1与模型2的差距较大，而在福州和汕头，两者的模拟结果相对接近。上述观测结果与图2中修正项均值的空间分布规律一致。正如图2（a）所示，绝大部分地区的修正项均值为正，表明模型1预测的台风强度变化小于模型2。在115°E至120°E与20°N至27°N区域，修正项均值接近于0，福州和汕头均位于该区域，说明这两座城市中，模型1和模型2的差异并不显著。在3个模型中，模型3模拟的CDF曲线与观测值最为接近，这说明充分考虑修正项能够提高模型再现台风强度概率分布的能力。
图66个城市的台风强度概率分布
Fig.6Probability distributions of typhoon intensity for six cities
表2列出了中国东南沿海地区部分选取站点及所研究城市的地理位置信息，根据选取的站点计算台风强度的均值和标准差，结果如图7所示。可以看出，模型2的预测平均值较模型1大1.70 m/s，模型3的平均值相比模型2大0.38 m/s。尽管模型3的平均值略低于观测值，平均差值约为1.50 m/s，但模型3是对原模型的显著改进。在标准差方面，模型1与模型2非常接近，表明修正项的均值对标准差的影响很小。模型3考虑了随机噪声的影响，因此，其标准差比模型2和模型1分别大0.93 m/s和0.77 m/s。此外，在海岸线里程大于1 200 km的站点，模型3与观测值吻合度更高，进一步验证了模型3的有效性。
表2沿海岸线站点和研究城市的地理位置信息
Tab.2Geographical location information of the coastal stations and study cities
图7沿海岸线站点台风强度均值和标准差的演变
Fig.7Evolution of mean and standard deviation of typhoon intensity at coastline stations
为了计算不同年份台风强度的重现值，采用以下公式计算T a重现值：
(13)
式中：F-1（x）为不同站点台风强度的逆累积分布函数，λ为最佳路径数据集中导出的台风年发生率。由于台风强度的观测样本有限，将这些观测数据拟合到Gumbel分布，以计算相应的分位数。
图8展示了x25、x50和x100随站点的变化情况。可以看出，模型3和模型2之间的差异与模型2和模型1之间的差异大致相当，这表明在台风强度模型中，修正项的确定性部分与不确定性部分都起着重要作用。此外，在整个沿海岸线站点上，模型1和模型2的x25、x50和x100与观测值之间的平均相对差异较大，分别为10.3%、10.1%、10.4%和5.1%、5.4%、6.1%；而模型3的相对差异较小，分别为3.2%、3.6%、3.9%。这表明，引入不确定性因素显著提高了台风强度模型的适应性、准确性和可靠性，为台风极值风速分析提供了更加科学的依据和支持。
图8台风强度的25、50、100 a重现值沿海岸线站点的变化
Fig.8Variation of 25, 50, and 100 a recurrence value of typhoon intensity at coastal stations
3.2 台风危险性分析
进一步检测本文台风强度模型在分析极值风速时的性能，利用该模型估算了中国东南沿海25、50、100 a重现期极值风速。基于上述历史台风强度演化过程模拟结果，对各个目标场址，当台风中心进入半径为250 km的圆形区域时，采用Vickery[21]中的台风风场模型模拟了台风风速，并获得各台风样本在目标场址引起的地表风速极值。进一步统计获得了25、50、100 a重现期的极值风速，如图9所示。可以看出，修正项对台风地表风速的影响与对台风强度的影响相似，模型1估计的极值风速小于模型2，而模型2的极值风速又小于模型3。
图9台风地表风速25、50、100 a重现值沿海岸线站点的变化
Fig.9Variation of 25, 50, and 100 a recurrence value of typhoon surface wind speed at coastal stations
4 结论
1）通过在已有的确定性动力学台风强度控制方程中引入修正项，建立了台风随机强度模型，同时引入偏态分布模型，更合理地反映了台风强度的随机特征。
2）发展了基于地理加权方法的台风随机源概率模型识别方法；对历史台风强度的模拟对比表明，原模型总体上低估了中国东南沿海地区的台风强度，修正项的引入显著改善了台风强度的模拟效果。
3）结合台风地表风场模型，分析了中国东南沿海地区25、50、100 a一遇台风地表风速，验证了本文建立的台风随机强度模型在极值风速分析中的有效性。
