台风的发生与强度分布

发表于 2019-07-28 分类于闭门造车阅读次数： Disqus：

前言

之前曾经介绍过一点台风的内容，包括台风的路径可视化等，此外在之前也曾介绍过极值分析的一些内容，在这里，我们对台风做了进一步的统计分析，分别用泊松分布和广义帕累托分布拟合了台风的年发生次数和巅峰强度。

1	%matplotlib inline

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import shapely.geometry as ageom
import cartopy.crs as ccrs
import glob
import altair as alt
import seaborn as sns
import matplotlib as plt

from datetime import datetime
from matplotlib.font_manager import FontProperties

alt.renderers.set_embed_options(theme='dark')
plt.style.use('ggplot')
mpl.rcParams['figure.figsize']=(16,9)
font = FontProperties(fname=r"c:\windows\fonts\NotoSansCJK.ttc", size=14)

台风数据整理

我们采用了由CEA提供的热带气旋最佳路径集，首先需要对原始的数据进行整理。为了便于标识每个台风，我们对热带气旋进行了重新编号，而没有采用其自带的编号。这是因为：

台风编号是标记和识别台风的方法之一。我国从1959年起开始对每年发生或进入赤道以北、140度经线以西的太平洋、和南海海域的近中心最大风力大于或等于8级的热带气旋(热带风暴及以上强度)每年按其出现的先后顺序进行编号。70年代后期以来，编号范围扩大为东经150度以西。目前编号范围扩大为东经180度以西。

为减少混乱，日本在1981年获委托为每个西北太平洋及南海区域内的达到热带风暴强度的热带气旋编配一个国际编号，但容许其他地区继续自行给予编号。自此，在大部分国际发布中，发布机构会把国际编号放在括号内（JTWC除外）。但是，各气象机构有时对热带气旋的编号会有差别，主要是因为其对热带气旋强度的评估有所不同。

def pre_process(destfile):
    """
    Separate the record information and track data
    :param destfile:
    :return:
    """

    file_data = pd.read_csv(destfile, sep='\s+', header=None,
                            names=['time', 'level', 'lats', 'lons', 'pres', 'wind',
                                   'ty_serial', 'chn_id', 'name'], dtype=object)
    year = ''.join(filter(str.isdigit, destfile))
    temp = file_data[file_data['time'] == '66666']
    for s in temp.index:
        file_data.loc[s:s + int(temp.loc[s, 'lats']), ['ty_serial', 'chn_id', 'name']] = [year + temp.loc[s, 'lons'],
                                                                                        temp.loc[s, 'pres'],
                                                                                        temp.loc[s, 'chn_id']]

    file_data.drop(temp.index, inplace=True)
    # ty_data.drop(ty_data.columns[[6, 7, 8]], axis=1, inplace=True)
    # ty_data['serial'] =
    return file_data

def merge_bst_files():
    file_list = glob.glob(r"./Data/CHBST/*.txt")
    # Detailed info at http://tcdata.typhoon.org.cn/zjljsjj_sm.html

    track_data = pd.DataFrame(columns=['time', 'level', 'lats', 'lons', 'pres', 'wind',
                                       'ty_serial', 'chn_id', 'name'], )
    for file in file_list:
        temp_ty_data = pre_process(file)
        track_data = track_data.append(temp_ty_data)

    track_data['time'] = track_data['time'].apply(pd.to_datetime, format='%Y%m%d%H')
    track_data[['level', 'lats', 'lons', 'pres', 'wind']] = track_data[['level', 'lats', 'lons', 'pres', 'wind']].apply(pd.to_numeric)
    track_data[['lats', 'lons']] = track_data[['lats', 'lons']] * 0.1
    track_data.to_csv('./all_tracks.csv', sep=',', encoding='utf-8', index=False)

# execute this method if necessary
# merge_bst_files()

总体分析

热带气旋与台风

track_data = pd.read_csv('./all_tracks.csv', dtype={'ty_serial': str})
ty_serials_list = np.unique(track_data['ty_serial'])
years_list = [int(serial[:4]) for serial in ty_serials_list]
years, freqs = np.unique(years_list, return_counts=True)
typhoon_info  = pd.DataFrame(zip(years, freqs), columns=['years', 'freqs'])

plt.bar(years, freqs)
plt.show()

D:\ProgramData\Anaconda3\envs\climate\lib\site-packages\IPython\core\interactiveshell.py:3049: DtypeWarning: Columns (7) have mixed types. Specify dtype option on import or set low_memory=False.
  interactivity=interactivity, compiler=compiler, result=result)

output_7_1

台风的的生成海域是大致可以集中到几个区域中的，假设热带气旋在每个区域中发生的次数服从泊松分布，并且在不同区域之间是相互独立的，整体还是服从泊松分布的，这样的话，可以估计局部和整体拟合分布的参数。由于泊松分布的到达率就是期望，所以用样本平均值作为参数的估计值。

此外，我们也可以试一下负二项分布。

import scipy.stats as st
poisson = st.poisson
nbinom = st.nbinom
mu =  np.mean(typhoon_info.freqs)
var = np.var(typhoon_info.freqs)
nbn_p = mu / var
nbn_r = mu*nbn_p/(1-nbn_p)

poi_rvs = poisson.rvs(mu, size=10000)
poi_pmf = poisson.pmf(np.arange(1,61), mu)
nbn_pmf = nbinom.pmf(np.arange(1, 61), nbn_r, nbn_p)
sns.distplot(typhoon_info.freqs)
plt.plot(poi_pmf, label='poisson fitted curve')
plt.plot(nbn_pmf, label='nbinom fitted curve')
plt.legend()
plt.show()

output_10_0

那么，两种分布的拟合优度如何呢？如同之前做的，我们希望可以用类似k-s检验的方式进行评价。但是，由于k-s检验不能用于离散分布，在这里我们采用皮尔森卡方统计量做拟合优度检验。具体可以分为以下几步：

计算卡方检验统计量， $\chi^2 = \sum{\frac{(O-E)^2}{E}}$ ，类似于观测频率与理论频率偏差的归一化平方和；
确定卡方统计量的自由度，对于拟合优度检验而言，df( the degrees of goodness-of-fit) = n(obs) - n(constraints)，例如对于泊松分布的拟合优度检验，H0假设是一个约束，参数为样本均值为另一个约束，也就是有2个约束；
选择需要的置信水平 $\alpha$ ；
将2计算的卡方统计量与给定自由度和置信水平的阈值比较；
如果 $\chi^2$ 超过阈值，则拒绝原假设（样本分布与假设分布不同）；否则，无法拒绝原假设，但是无法拒绝不代表接受。

我们通过泊松分布的卡方统计量拟合优度检验来进行解释：

1
2
3

years_num =len(typhoon_info.freqs)
max_obs = np.max(typhoon_info.freqs)
years_num, max_obs

(70, 53)

chi_squared_df = pd.DataFrame(columns=['rv', 'obs_freq', 'pmf', 'expected_freq'])
chi_squared_df['rv'] = np.arange(0,max_obs+1)
chi_squared_df['obs_freq'] = np.asarray([np.sum(typhoon_info.freqs==i) for i in np.arange(0,max_obs+1)])
pmf = np.zeros(max_obs+1)
pmf[:-1] = [poisson.pmf(i, mu=mu) for i in np.arange(0, max_obs)]
pmf[-1] = 1 - poisson.cdf(max_obs-1, mu=mu)
chi_squared_df['pmf'] = pmf
chi_squared_df['expected_freq'] = pmf * years_num
chi_squared_df.tail()
chi_sq_statistic = np.sum((chi_squared_df['obs_freq']-chi_squared_df['expected_freq'])**2
                          /chi_squared_df['expected_freq'])
chi_squared_df.plot.bar(x='rv', y='obs_freq')
chi_sq_statistic

64.03201336457617

output_13_1

1	st.chi2.isf(q=0.05,df=68)

88.25016442187413

根据上面的分析可知，卡方统计量小于置信水平为0.05时的阈值，所以无法拒绝原假设，即台风发生频次服从泊松分布的假设。

但是，如果我们假设台风频次服从负二项分布，会发现卡方统计量甚至更小，所以也无法拒绝。所以，到底选择哪个假设，一方面需要数据做支撑，另一方面背后的物理机理和模型的目的同样也非常重要。此外，观察上面的图，我们还可能会发现：

泊松分布存在时间上具有不均匀性（inhomogeneous），从时间上看93年前后的水平也有显著不同，而且在30~32之间的样本数明显偏少。
其他气象因素，例如厄尔尼诺可能会影响台风的发生频次。

强度分析

巅峰强度分布

所谓“巅峰强度”其实是我自己起的名字，我不知道专业的术语中文叫什么，英文应该是Lifetime Maximum Intensity (LMI)，顾名思义可知，是在一个气旋的整个生命周期中所能达到的最大强度，例如以近中心2分钟平均最大风速，反映了一个台风的威力大小。

1
2
3

typhoon_details = pd.DataFrame(columns=['ty_serial', 'level', 'max_wind', 'lifetime'])
typhoon_details['ty_serial'] = list(grouped.groups.keys())
typhoon_details['max_wind'] = grouped['wind'].max().values

1
2
3

reduced_track_data = track_data[track_data.level<9]
reduced_grouped = reduced_track_data.groupby('ty_serial')
typhoon_details['level'] = reduced_grouped['level'].max().values

1	print(typhoon_details[['level','max_wind']].describe())

             level     max_wind
count  2304.000000  2304.000000
mean      3.507812    35.471788
std       1.725629    18.074845
min       1.000000    12.000000
25%       2.000000    20.000000
50%       3.000000    30.000000
75%       5.000000    45.000000
max       6.000000   110.000000

从上面的统计结果我们可以看到，在1949年以来的2304个热带气旋（含热带低压）中，最弱的热带气旋其巅峰强度也不过12m/s，也就是6级强风的样子，而最强的热带气旋巅峰最大风速可达110m/s（396km/h)，也就是比复兴号跑得还要快，对应的则是超强台风级别。此外，平均的巅峰最大风速为35.47m/s，也就是台风级别，中位数为30m/s，那么达到台风级别的热带气旋不足半数。

那么，热带气旋的强度有什么分布规律吗？我们可以先通过作图进行直观的判断。

max_wind_summary = typhoon_details.groupby('max_wind')['ty_serial'].count()
fig, axes = plt.subplots(1,2, figsize=(24,8))
sns.distplot( typhoon_details['max_wind'], bins=10, ax=axes[0],hist_kws=dict(edgecolor="w", linewidth=2))
axes[1].bar(max_wind_summary.index, max_wind_summary, label='frequqencies of max wind')
plt.legend()
plt.show()

output_23_0

左图是我们用seaborn画出来的直方图和对应的核密度估计的概率密度函数图，乍一看，规律性还是很强的：最大风速（巅峰强度）越大，数量越少。

但是，由于在CMA最佳路径数据集中，风速均为整数，如果我们看每一个风速的频数，则会惊喜地发现：经常是某一个风速的频数显著高于其两侧的，如此成“波浪”形态。

咦，难道说台风对特定的风速值有偏好？

我倒不觉得大自然会如此挑剔，更偏向于人为的原因，比如计算或者记录的偏差。当然，如果真的是天然这样，那也太有趣了！

分布拟合

查阅相关文献可知，学术界通常使用Weibull分布或者广义帕累托分布（Generalized Pareto distribution, GPD)来对热带气旋的极端风速进行建模。我们以广义帕累托分布为例。

在拟合之前，我们首先对数据进行清洗，由于广义帕累托分布主要考虑极端情况的分布，需要有一个阈值，在这里我们以台风级(32.7m/s)作为阈值。同时，前面观察到风速取在5的整数倍的情况远多于其他情况，为此我们将采样规整为5的整数倍。

cleaned_max_wind = pd.DataFrame(columns=['max_wind', 'freq'])
step_max_wind = np.arange(35, max_wind_summary.index[-1]+5,5)
temp = np.zeros([len(step_max_wind),2])
obs = []

for i in range(len(step_max_wind)):
    temp[i,:] = [step_max_wind[i], max_wind_summary[(max_wind_summary.index<=step_max_wind[i])
                                                   &(max_wind_summary.index>step_max_wind[i]-5)].sum()]
    obs += [temp[i,0]]*int(temp[i,1])
    
cleaned_max_wind = pd.DataFrame(temp, columns=['max_wind', 'freq'])
cleaned_max_wind = cleaned_max_wind.astype(dtype={'max_wind':'float64', 'freq':'int64'})

cleaned_max_wind.head()

	max_wind	freq
0	35.0	202
1	40.0	202
2	45.0	175
3	50.0	148
4	55.0	99

import lmoments3 as lm
from lmoments3 import distr
paras = distr.gpa.lmom_fit(obs)
paras

OrderedDict([('c', -0.29642660282781813),
             ('loc', 33.48293631100662),
             ('scale', 21.983950160585245)])

fitted_gpa = distr.gpa(**paras)
r_x = np.linspace(35,150,100)
r_y = fitted_gpa.pdf(r_x)
sns.distplot(obs, bins=10,hist_kws=dict(edgecolor="w", linewidth=2))
plt.plot(r_x, r_y, label='fitted by GPD')
plt.legend()
plt.show()

output_31_0

1 2	samples = fitted_gpa.rvs(len(obs)) st.kstest(obs, 'genpareto', args=(paras['c'], paras['loc'], paras['scale']))

KstestResult(statistic=0.10876338638750287, pvalue=2.8521788044114453e-12)

我们看到，虽然估计的曲线和频率直方图以及he密度估计都拟合的较好，但是原假设并没有通过ks检验，这不得不说是件悲伤的事情。出现这种问题的原因可能是我们选择的阈值不太合适，也可能假设确实不成立，还可能是因为样本的取值问题，当然可能有什么原因我也不知道。

小结

我们利用CEA的台风最佳路径集，对1949年以来西北太平洋的热带气旋进行了简单的分析，可以发现，受限于数据的质量、实际系统的复杂性等问题，我并没有得到理想的结果或者可靠的结论，居然有一点挫败感和无力感。

文末福利

wind_level

前言