曹诗嘉,侯静惟,方伟华,林 伟
(1.北京师范大学 地理科学学部 环境演变与自然灾害教育部重点实验室,北京100875;2.北京师范大学 地理科学学部 减灾与应急管理研究院,北京100875)
中国是遭受台风灾害最严重的国家之一,常因台风造成严重的人员伤亡和巨大的经济损失。例如,2006年第8号台风“桑美”在浙江省苍南县登陆,引起了狂风巨浪,共造成483人死亡,千余条船沉没,直接经济损失达196.5亿元人民币[1],2014年10号台风“威马逊”登陆风速达到17级,造成海南省216个乡镇(街道)受灾,受灾人口325.8万人,直接经济损失108.28亿元[2]。近年来我国沿海台风易发区经济快速发展,对于台风致灾因子的暴露程度不断增加;而沿海桥梁、铁路、石油石化以及核电站等重大工程的选址及设计需要对局地风速进行评估。因此,面向台风灾害风险管理及工程设计需求,采取有效手段对中国台风大风危险性进行评估具有重要意义。
目前学界对于台风大风危险性评估方法主要有三类,其中第一类一般通过极值分布模型对历史台风风速序列进行模拟从而计算单点的风速重现期;第二类通常针对某个区域,通过统计该区域历史台风参数概率分布,并结合风场模型和极值分布模型计算该小区域的风速重现期[3];第三类一般采用全路径模拟的方法生成大样本台风随机事件,然后结合风场模型和极值分布模型计算台风影响区域的风速重现期[3]。与第二类和第三类方法相比,第一类方法在历史风速观测时间序列较长时,可充分利用历史样本信息得到较为准确的风速重现期计算结果,其缺点在于当历史数据不足时,重现期计算结果存在一定不确定性。国外利用第一种方法即历史观测数据进行台风大风危险性评估采用的数据类型主要包括台风近中心历史风速观测数据及气象站点历史风速观测数据两类。一些学者利用台风近中心历史风速观测数据采用韦伯分布(Weibull Distribution)、广义帕累托分布(General Pareto Distribution,GPD)等极值分布函数对美国区域最大风速进行拟合,并计算了一定空间尺度的年发生超越概率[4,5],此类方法可以满足台风综合危险性评估需要,但其不足之处在于台风近中心历史风速数据无法反映整个区域的空间差异性,因此不适宜用于大区域范围台风灾害危险性评估;另外一些学者采用气象站点历史风速观测数据,同样利用极值模型对强风危险性进行评估[3],此类方法可以较好地反应风速分布的空间差异性,缺点在于气象站点历史风速观测数据往往存在数据量不足的问题。
中国大风危险性研究多是针对小区域或台站进行的不同分布模型对比研究或大风重现期推算实证研究[6-9],如综合岛屿站观测、台风记录、船舶报告以及数值模式计算风速,组成大风年最大值序列,利用矩估计法拟合大风序列,通过检验选定最适风速序列及分布模型,绘制中国近海50年一遇和100年一遇大风极值等值线[10],此类研究的缺点在于面向区域较小,且缺少对评估结果不确定性的讨论。利用气象站点历史风速观测数据计算重现期,数据不足、数据一致性问题、测量误差是结果不确定性的3个最主要来源。为了提高结果可靠性,需要对不确定性进行量化,目前国内已有研究采用广义极值分布模型,利用极大似然法及渐进分布理论推导出不同年遇水平设计风速和一定置信度下的置信区间,并通过气象站点年最大风速资料进行实证研究[11],该研究对单个站点重现期计算结果的不确定性进行了刻画,对于工程的选址设计具有一定意义,而为了反映整个中国东部沿海台风大风危险性分布特征及不确定性大小,还需对站点计算结果进行插值,从而得到整个空间连续的重现期评估结果。
本文的研究目标为利用气象站点历史风速观测资料,评估不同重现期下中国东部台风风速空间分布并给出一定置信度下的风速置信区间。具体过程为,首先,基于中国地面气象站点历史风速观测数据,提取每个站点历史台风过程影响风速的年最大值;其次,针对每个站点利用极大似然法估计3种经典极值分布模型的参数估计值,分别选择最适模型进行站点台风年最大值风速拟合;然后,利用bootstrap重采样方法,得到90%置信度下的的风速置信区间;最后,从4种空间插值方法中选择最优者插值得到不同重现期及90%置信度下的中国台风大风空间分布并对中国台风大风危险性区域分布及不确定性特征进行分析。
本文采用的风速资料来自中国756个一般气象观测站1951-2014年日值风速观测值,指标为日极大风速(3 s瞬时风速的日最大值)。由于中国东部沿海地区是受台风影响的主要地区,因此选取中国东部、中部、南部13个省、直辖市、自治区作为本文的研究区,考虑到研究区边缘插值风速的准确性及连续性,故保留中国东部、中部、南部21个省、直辖市、自治区的359个一般气象观测站的气象观测数据,359个一般气象观测站站点空间分布如图1所示。
图1 研究区地面一般气象观测站点分布图
本文采用中国气象局上海台风研究所(CMA-STI)整编的西北太平洋热带气旋最佳路径数据[12],包括1949-2014年共2182场台风,其中登陆中国台风共计611场。数据具体指标包括:国际热带气旋编号、中国热带气旋编号、英文名称、每个台风路径点的年、月、日、时、等级、经度、纬度、近中心最低气压以及最大持续风速等,记录时间间隔为6 h。
台风大风重现期计算的前提是提取每个气象站点受台风影响风速,其提取方法简述如下:①基于站点经纬度信息,提取一定空间范围内历史台风的路径点;②获得历史台风对该站点影响的起止时间,并提取该时间范围对应的站点观测风速,作为历史台风对该站点的影响风速。这种方法的关键在于台风影响范围的界定,由于台风是由比较均匀的热带海洋气团发展起来,因此台风气压场、风场分布具有一定的对称性,可近似将台风看作圆对称的涡旋,其半径变化范围小到上百公里,大到上千公里[13]。可基于台风中心位置,经验地设定台风影响半径从而划定影响范围。若半径设定太小,则无法获取台风外围大风信息,可能导致获取的格点风速样本不足;若半径设定过大,虽可获得完整的台风大风序列,但可能将其他天气系统引起的局地大风误认为台风大风[14],在本研究中,经验性地将台风影响范围设定为距台风中心500 km。
利用上述气象站点台风影响风速提取方法,对359个气象站点的历史台风影响风速进行提取,对部分历史样本小于15个的站点予以剔除,最终得到风速数据提取结果。图2为站点历史台风过程影响风速最大值统计图,从图中可以看出,台风风速样本数超过15个的站点共有275个,且空间分布相对较广泛,其最大值出现在浙江大陈岛站,为59.5 m/s,站点台风风速高值区主要分布在东南沿海地区,部分内陆地区也出现了风速高值,如山东泰山站及安徽黄山站,体现了地形对于风速的影响作用。
图2 359个一般气象观测站历史台风风速过程最大值
空间插值方法主要包括点插值和面插值两种,已知某点数据推求空间区域内任一点数值应采用点插值方法进行插值[15],点插值方法又包括了克里金法、多项式回归法、最近邻法、反距离权重法等。为确定最适用中国东部台风大风风速插值的方法,以2.1中得到的站点历史台风风速最大值作为检验样本,选用克里金法、反距离权重法、自然邻域法以及最近邻法4种点插值方法对站点风速进行空间插值,然后通过交叉验证法确定最适插值方法。各插值模型主要参数设置如下:最大搜索半径设为500 km,搜索范围内的最小样本数设为15,其中反距离权重法采用固定搜索半径,幂参数设为3,克里金插值法的半变异函数采用指数模型。
具体计算方法为对站点i,利用其他站点的记录,选用一种方法计算该站点的插值风速,并将其与该站点的实测值进行比较,循环计算所有站点的均方根误差(RMSE)以及决定系数(R2),以所有站点RMSE及R2的均值作为插值效果比较的判定标准,均方根误差及决定系数的公式分别如式(1)
(1)
(2)
利用达标站点台风过程风速的最大值的统计结果作为检验样本,比较4种空间插值方法模拟结果的交叉检验误差,如表1所示。可以看出,对于台风风速的插值,克里金插值法效果最好,站点平均RMSE为6.45,R2达到了0.41。
表1 4种站点风速空间插值方法交叉验证比较
气候统计学中,通常采用经典极值理论对气象要素极值进行拟合[16],经典极值理论包括Gumbel模型, Frechet模型以及Weibull模型3种,其累积概率分布函数分别如式(3)~式(5)所示,其中μ为位置参数,σ为尺度参数,α为形状参数,x为连续型随机变量,超过定值x则表示极端事件发生。根据历史样本年最大值分布计算极端事件的重现期(Return Period,RP)是极值统计最重要的应用之一,由式(6)可知每一个重现期对应一个极值分位数,表示极端事件的极值变量的数值大小。同时,对于给定重现期,极值分位数越大说明超越概率越小,则极端事件发生的可能性也就越小。
图3 (No.59758)海口站3种极值函数拟合下的台风年极值风速累积概率及90%置信区间
表2 (No.59758)海口站不同极值分布拟合参数及拟合效果检验比较
对于每一个气象站点,利用年极值抽样方法,选择该站点历史台风影响风速的极大值作为极值拟合的样本,然后从3种极值分布模型中选择最适者进行极值拟合,得到模拟风速的期望值(即位于所有风速估计值中50%分位数的风速值)并利用bootstrap方法得到90%置信度下的风速区间。
(3)
(4)
(5)
(6)
为定量刻画不同极值模型的拟合效果,①采用Kolmogorov-Smirnov(KS)检验方法评价拟合结果与样本概率分布的一致性。KS检验根据样本得到经验分布与假设分布函数之间的偏差,描述两个独立统计样本相似程度[17-18]。统计量D及其对应的P值可用于表征模拟分布与样本分布的相似性,P值越高,则相似性越大;②针对每个站点,计算每个模型模拟CDF与经验CDF的均方根误差RMSE,用于协助判断拟合效果。RMSE越小模拟平均误差越小,若出现不同模型的KS检验结果相同的情况,则选择RMSE最小者作为最适模型。
以海口站点为例,利用上述站点大风重现期及不确定性估计方法,得到3种极值模型模拟下的台风年极值风速累积概率分布90%置信区间与观测值的对比图。由图 3可知海口站点采用Weibull分布函数的拟合效果较好,观测值均位于模拟风速的90%置信区间内,特别是对于样本累积概率分布尾部的拟合效果较好。Gumbel函数及Frechet函数的拟合效果较差,样本累积概率分布尾部超出了模拟风速的90%置信区间。
图4 台风极大风速数据达标站点及其最适极值分布类型
图5 不同重现期及90%置信度的中国台风极大风速空间分布
表2给出了利用三种极值函数对海口站台风风速年最大值进行拟合得到的模型参数[19],以及对应的KS检验及RMSE结果。可以看出,对于海口站年极值风速拟合,采用Weibull分布的D统计量最小、P统计量最大、RMSE结果与Gumbel模型结果相差较小,因此对海口站点台风年极大风速适合采用Weibull函数进行拟合。其他站点最适极值模型的确定方法与海口站点相同。
对359个气象站点,利用2.3中的站点大风重现期估计方法,提取台风影响期间站点风速的年最大值作为极值拟合的样本,利用Gumbel、Frechet以及Weibull函数进行拟合。计算过程中,为了减小因数据量不足导致的参数估计的不确定性及满足极值分布函数对拟合样本数量的要求,对极值拟合样本数少于15个的站点予以剔除,最终得到台风极大风速数据达标站点及其对应的最适极值分布(如图4所示)。可以看出: ① 359个站点中共计133个站点满足拟合样本数量需求,最适极值分布确定为Gumbel、Frechet以及Weibull的站点个数分别为60个,27个以及46个,即对于不同站点的风速极值,应采用不同极值分布函数进行拟合;② 359个站点中位于研究区外符合样本数量要求的站点较少,这是因为这些区域的站点受台风影响频次较少,影响时长较短;③ 359个站点中位于研究区内江西省、福建省的站点满足拟合样本数量要求的数量偏少,江西18个站点中仅3个站点满足要求,福建省22个站点中仅 5个站点满足要求,造成这一结果的主要原因是尽管此区域台风影响频次较多,影响时间较长,但未达标站点历史台风观测数据较少。
利用各达标站点最适极值分布模型,分别计算各站点20年一遇及50年一遇的台风风速期望值及90%置信区间;然后,利用克里金法插值得到中国1 km网格分辨率的20年一遇及50年一遇重现期风速期望值及90%置信区间,计算结果如图5所示。从图5可以看出:① 台风极值大风高值区主要分布在我国东部地区,尤其是东南沿海地区,该区域是我国历史上受台风灾害影响最严重的地区,发生台风大风极端事件的可能性也最大;② 部分内陆地区也出现了风速高值区,主要原因是由于历史台风影响的观测风速样本不足;③ 典型重现期台风风速估计的不确定性较大,20年一遇的台风风速90%置信区间约49 m/s,50年一遇的台风风速90%置信区间近61 m/s,置信区间的高值区主要分布在中国东南沿海,置信区间整体由东南沿海向内陆递减。
本文基于中国地面气象站点历史风速观测数据,提取了每个站点历史台风过程影响风速,在此基础上优选出最适极值分布以及空间插值模型,计算了不同重现期下,中国东部地区空间上连续分布的台风风速及90%置信度下的风速区间,并对其区域分布特征及不确定性进行了分析。主要结论如下:
(1)采用交叉检验的方法比较了4种空间插值方法对于站点极大风速的插值总误差,结果表明,对于站点台风极大风速插值,克里金法插值效果最优。
(2)由于历史样本数量不足、数据一致性问题及测量误差等原因,基于站点历史风速观测数据的台风大风重现期估计结果不确定性较大,为了保证评估结果的可靠性,需要对其进行不确定性量化。
(3)结合极值分布模型以及空间插值方法评估中国台风大风危险性,相较于传统小区域或台站尺度的评估有一定改进,但仍存在一些局限性,包括:① 气象站点历史风速观测样本不足可能影响极值风速评估精度及区域插值效果;② 本文经验地将台风影响范围确定为500公里,主观性较强,今后的研究可通过敏感性分析方法对台风影响最适范围进行评估[14];③ 由于台风登陆后下垫面情况比较复杂,局地地形以及地表粗糙度的变化对于近地表风速影响较大,仅依靠空间插值的方法难以精确刻画该影响机制。后续的台风大风危险性研究可基于大样本台风路径事件集,结合台风风场模型得到空间上连续分布并且样本充足的台风风速序列[7-9],在此基础上进行不同重现期极值风速评估将更具可靠性[19]。