台风浪集合预报方法研究

刘安琪，李铖，郭文云，葛建忠*

（1.华东师范大学 河口海岸学国家重点实验室，上海200241；2.上海市海洋监测预报中心，上海200062；3.上海海事大学 海洋科学与工程学院，上海201306）

0 引言

台风期间，受台风影响的海域的环流、水位、波浪等要素通常会发生剧烈变化，其中水位的异常变化会导致风暴潮现象的发生，而剧烈的海表风场会激发产生巨大的波浪，即台风浪现象[1]。与风暴潮相比，台风浪作用周期更短，但产生的瞬时破坏作用更大，会破坏近岸工程，并威胁人们的生命以及财产安全。2019年，我国近海共发生有效波高4.0 m（含）以上的灾害性海浪过程39次，因灾直接经济损失达0.34 亿元，死亡（含失踪）22 人。台风浪在向近岸传播的过程中，会与复杂海底地形以及不规则岸线发生非线性相互作用[2]，使得波浪预测难度增大，另外近岸波浪观测资料的缺乏也限制了区域波浪预报产品预报精度的提升。

对台风浪的预报，国内外很多学者从不同的方向做了很多努力，取得了丰硕的成果。CHAWLA等[3]利用了美国国家环境预报中心（National Centers for Environmental Prediction，NCEP）的CFSR（Climate Forecast System Reanalysis）再分析数据，使用第三代风浪模型（WAVEWATCH III）进行了31 a 的后向预报，分析发现，在南半球，CFSR对风速的过度预测对波浪预测的准确性有重要影响。ALVES 等[4]基于NCEP 和美国海军数字气象和海洋学研究中心（Fleet Numerical Meteorology and Oceanography Center，FNMOC）的数值模式对有效波高等要素进行综合概率预报，结果表明多模式联合预报方式的准确性比单一数值模式有了很大的提升。FAN 等[5]使用多个全球大气模式数据作为波浪模型输入得到的有效波高，与空间上不同近岸浮标的观测值的吻合程度各不相同。PAN 等[6]提出了一种优化模型-集合波高的算法，计算了不同风模式输入波浪模型产生的波高，该算法在空间上给不同风模式赋不同的权重，调整后的波高在近岸范围内与浮标观测值的吻合较好。ZIEGER 等[7]使用欧洲中期天气预报中心集合预报模式（The European Centre for Medium-Range Weather Forecasts-Ensemble Prediction System，ECMWF-EPS）中经过矫正的大气风场生成的集合风场进行波浪模拟，并对比了近岸区域的浮标观测数据，部分风场的结果和观测的有效波高峰值吻合较好。

综上可见，近岸波浪预报结果的准确性在很大程度上依赖于风场的准确性。近岸波浪由于各物理过程的非线性作用导致不确定性很强，很难找到一种传统模式或者参数化方案在整个预报时段内都表现较好。在实际预报中，既要关注预报结果的准确性，也要关注预报结果的概率区间预测。本文旨在提出一种新的台风路径集合方法开展波浪预报实验，并给出台风浪的极值范围作为预报参考。

集合风场目前主要有以下几种生成方法：一种是使用单一业务预报中心的数据作为数据源（如应用较多的增长繁殖法[8]）以及特定业务中心集合预报系统中的集合成员[7]，这种方法的优势在于数据一致性较好，不足之处是会引入单一模式的固有误差；另一种是使用多业务预报中心数据成员，基于“共识”的方式挑选出其中较好的成员[9]；第三种是基于不同业务中心历史预报误差分配权重，使用“概率圆”的方法进行集合概率预报[10]，这种方法能避免单一模式的偶然误差，不足在于各预报源的“好坏”一般由其历史预报误差确定，无法在预报过程中实时变动；还有一种是使用研究区域内历史台风运动数据作为数据源，并基于台风的历史最佳路径集数据确定所有状态的概率转移情况，从而在实际预报过程中对台风移动的路径进行模拟[11]，这种方法的不足之处在于对实时预报数据的利用较少，且状态转移概率情况受限于研究区域。这些方法都只聚焦于历史数据或实时预报数据，没有同时兼顾到两者。

本文使用实况台风数据融合背景分析风场，驱动近岸波浪模式以及海洋模式进行浪-流耦合模拟，使用后报方式验证数值模式；在数值模式精度良好的基础上，提出一种新的台风预报路径集合生成方法，考虑研究区域内历史台风数据和业务中心实时预报数据，构建集合路径并形成集合风场，开展近岸区域波浪要素计算的数值实验，并考虑集合路径情形下有效波高计算结果的覆盖率，以此提供一种新的台风浪概率数值预报方法。

1 研究区域和研究方法

1.1 研究区域与数据概况

长江口区域处于亚热带季风气候区域，最大风速大多发生在夏秋季的台风时期。该区域内的波浪以风浪为主[12]，由于入海口分叉较多，地形和岸线较复杂，区域内不同位置的波高差异性很大，最大波高大多出现在台风期间。据统计，平均每年有2～3 个台风对该区域产生直接影响[13]。长江口区域存在诸多重要的近岸站点（见图1）。本文采用H1与H2两浮标平台的实测波高数据，采样频率为1 h，其中H1位于杭州湾北岸近岸区域（121°33′12.05″E，30°44′52.84″N），H2 位于长江口外30 m 水深处（122°48′29.29″E，30°58′12.19″N），两个站位兼顾了外海与近岸的波浪差异及变化。

图1 研究区域及主要岸基及浮标站点Fig.1 Study area and coastal and buoy observation sites

1.2 数值模式

本文选择数值模式作为研究方法，波浪模式选用第三代海浪数值模式（Simulating Waves Nearshore，SWAN），考虑风能输入、波浪折射、三波和四波相互作用、深度诱导波破碎、底摩擦等因素；风场输入采用JANSSEN 风增长方式，白冠耗散使用KOMEN 模式，底摩擦采用JONSWAP 模式进行参数设置。考虑到潮汐作用，除为波浪模式输入风场外，还必须为其输入潮流和潮位作为驱动要素，其由海洋模式进行计算，从而进行高精度的浪流耦合计算（见图3）。

海洋模式选用有限体积海洋数值模式（Finite-Volume Community Ocean Model，FVCOM）。它是一种无结构三角形网格架构、有限体积自由表面的三维海洋数值模型，其原始方程包括动量、质量连续方程，以及温度、盐度和密度方程。和结构网格海洋模型相比，三角形网格在岸线相接处有局部加密效果，且过渡更加平滑，在长江口杭州湾区域已广泛被使用[14-15]。

为提高预报时效并保证计算精度，SWAN 模式采用大小区域嵌套计算的方式。大区域覆盖了研究区域内台风从产生到对近岸造成影响的整个过程，能给小区域提供合理的计算边界（见图2a）。为兼顾计算效率，采用正交曲线网格，覆盖范围为23°～42°N，115°～130°E，包括了东海、黄海、渤海和部分西北太平洋洋面，开边界处分辨率为3 km 左右。近岸区域，由于岸线和地形的复杂性很高，对波浪要素的影响十分显著。为更加精细地刻画各种局地变化，模拟近岸波浪，保证计算结果的准确性，小区域（见图2b）采用高分辨率无结构三角形网格，主要包括长江口以及杭州湾，河流上边界延伸至大通站，在河道以及沿岸处进行加密，分辨率为100 m 左右。在预报计算时，大区域计算较快，将其波浪谱结果作为小区域的计算边界，并在小区域内进行精细化计算，从而在保证计算结果准确性的同时兼顾计算效率。浪-流耦合计算框架见图3。

图2 研究区域内SWAN模式计算网格Fig.2 SWAN model grid in the study area

图3 浪-流耦合计算框架Fig.3 Wave-flow coupled computing framework

2 集合风场的构建

2.1 集合路径的形成

本文中以24 h、48 h、72 h 3个预报时刻为代表，构建集合路径并形成集合风场。集合路径的构建采用以下步骤：首先形成某一时刻的预报位置集合，并为它们计算分配相应的权重，再将3个代表性预报时刻各自位置集合进行组合，形成所有的集合路径，并计算出每条路径的权重。整个流程的重点是对某一时刻各预报位置的确定和权重分配，本文采用如下方法：统计业务预报中心（中央气象台）历史台风位置预报误差，并形成当前预报时刻位置偏差的概率分布情况，选取两个特征概率从而确定两个特征偏差距离；同时，使用台风历史最佳路径集数据，基于面积相似法[16]优选出若干相似路径并综合距离相似及时间相似搜寻出预报时刻相似路径上的对应位置，统计相似位置相对当前时刻预报位置的方向，这些方向代表着台风移动方向以及移动速度的不同趋势，使用聚类分析的方法确定若干个主要聚类方向以及各自的权重；将业务预报位置作为中心，偏差距离作为特征半径，聚类方向结果作为特征方向，形成某一预报时刻的特征预报位置集合。

2.1.1 代表性位置偏差计算

本文搜集了中央气象台2008—2020 年东海区域共318场台风的历史预报数据以及对应的实际台风过程数据（网址：http://typhoon.zjwater.gov.cn）；以实时起始时刻和预报时刻的位置形成的区域筛选出此区域内相同预报时效的历史预报数据，统计各预报时效下历史位置预报偏差的累计概率情况。预报偏差EDis为根据两点经纬度计算出的距离误差，公式为：

式中：R为地球半径;a为两点纬度之差；b为两点经度之差；φ1和φ2为两点纬度。

各预报时效位置偏差的累计概率曲线见图4。从图中可以看出，随着累计概率的增大，各预报时刻位置偏差增大的速率都越来越大。本文选取70%和95%概率为典型代表作为特征概率，确定各预报时刻两个特征偏差距离作为各预报时刻的两个特征半径，结果见表1。

表1 特征偏差距离（单位：km）Tab.1 Feature deviation distance(unit：km)

图4 各预报时刻位置偏差的累计概率Fig.4 Cumulative probability of position deviation at each forecast moment

2.1.2 主要聚类方向计算

搜集了1949—2020 年西太平洋区域台风最佳路径集资料（网址：http：//tcdata.typhoon.org.cn）[17-18]，对实际台风过程中的各预报时刻，基于面积相似法[16]，以历史台风轨迹和实时台风轨迹上所有路径点为分割点将整个面积分成L段，采用积分方式计算历史台风路径和当前路径围成的面积，并作为相似度依据。计算公式为：

式中：S为总面积；Si为面积微元。

采用这种方式，当预报台风、预报起始时刻以及预报时效发生变化时都会导致相似路径成员发生变化。获取所有相似路径后，选择时间和距离最近的对应预报时刻的相似位置，将业务预报位置和所有的相似预报位置投影为二维平面上的点，计算出所有相似位置相对于业务预报位置的方向角（0°～359°）。将这些方向角作为数据，使用一维均值聚类方法形成特征聚类方向。计算方法为：

式中：样本共划分k簇；Ci为第i簇；ui为簇Ci的均值向量，也称为质心；x为样本点；E为平方误差。

式（3）中由于k值不同，计算出的E也不同。根据经验，一般将k从2～10 赋值试算E，随着k增大，E呈下降趋势，以E在下降过程中梯度最大处的k值作为最优的聚类个数（簇个数）。每个簇质心的方向角代表一个特征方向，该簇包含的方向角数量在所有方向角数量中的占比为该特征方向的权重。

2.1.3 形成集合路径

获得各预报时刻的特征半径和特征方向后，对某一特定预报时刻而言，以业务预报位置为中心可以形成特定时刻所有的特征成员位置；从每个预报时刻中各取一个成员位置，再连接所有选出的成员位置作为一条路径，集合路径的数量为各预报时刻成员位置进行组合的所有情况，每条路径的概率也随之确定（见图5）。具体流程如下：

图5 台风集合路径（0～48 h）Fig.5 Typhoon track ensembles(0～48 h)

①对特定预报时刻形成特征半径及特征方向，获得各特征方向的权重。以24 h 为例，基于历史业务预报数据确定指定概率下特征半径为Rt1；基于相似路径数据，使用聚类分析方法收敛到最佳聚类个数Kt1（Kt1为式（3）中最优k取值，不同数据集最优值不同），每个聚类的中心作为一个特征方向，每个聚类样本量占总样本的比例作为特征方向的权重（Wt11,Wt12…,Wt1Kt1），对应的特征点空间位置为Pt11,Pt12…,Pt1Kt1。

②对24 h、48 h和72 h预报时刻，从每一个时刻选择一个特征位置形成一条成员路径，例如Pt1Kt1→Pt2Kt2→Pt3Kt3，相应地，这条路径的权重为Wt1Kt1·Wt2Kt2·Wt3Kt3，3 个预报时刻一共可以形成Kt1·Kt2·Kt3条路径，计算每条路径的权重并进行集中构成集合路径。

2.2 集合风场的形成

在具有确定性的台风实况路径的后报验证中，将台风风场模型融合背景风场作为风场输入，而预报实验中仅考虑台风模型风场。

台风风眼处风速较小且向外迅速增加，超过眼壁后风速又急剧减小直至消失。本文采用对称风场模型刻画台风风场，其气压场由藤田公式[19]给出：

式中：r为距离台风中心的距离；P为空气质点的气压，Pe为环境气压，Pc为台风中心气压；R为最大风速半径。通过藤田公式可以计算出整个台风风域内的气压分布情况。

台风风速根据梯度风关系结合气压分布给出[20]，计算公式如下：

式中：f= 2Ωsinφ为科氏参数，Ω 为地球自转角速度，φ为纬度；ρa为空气密度。

在后报验证中，本文使用CFSR 作为背景风场。CFSR 是第三代再分析产品，它是一个全球性的、高分辨率的、大气-海洋-陆地-表面-海冰耦合系统，旨在提供这一时期这些耦合域状态的最佳估计。该数据经纬度分辨率可以达到0.5°，时间分辨率可以达到1 h。两种风场融合公式为[21]：

式中：WI为台风模型风场计算风速；WQ为CFSR 风场风速；r为网格点与台风中心的距离；e、c为系数；R为最大风速半径，计算公式为[21]：

式中：Pc为台风中心气压。

3 应用示例

本文以2019 年9 号台风“利奇马”为例，使用实况台风路径对数值模式进行后报验证；选择2019年8 月9 日00 时（北京时，下同）作为起始预报时刻进行72 h的虚拟集合预报实验。

台风“利奇马”于8月7日05时被中央气象台升格为台风，23 时升级为超强台风，此后台风一直向西北方向移动，靠近浙江沿岸地区，8 月10 日01 时45分左右在浙江省温岭市城南镇沿海登陆，此时中心附近最大风力达到16级，中心最低气压为930 hPa，随后台风继续移动至黄海海域，8月11日20时50分左右在青岛市沿海再次登陆，此时中心最大风力仍有9 级，中心最低气压为980 hPa，此后台风移动至渤海海面并不断减弱直至消亡。

自起始预报时刻，构建24 h、48 h、72 h 的集合路径。选取从台风自洋面生成位置—各时刻业务预报位置区段作为实时轨迹，并与历史轨迹做相似度计算，计算出24 h、48 h、72 h 预报时刻情形下相似度最高的各50 条路径（见图6）。从图中可以看出，不同预报时刻的实时轨迹不同，相似路径的选择结果也不同。获取相似路径后，采用聚类分析方法确定各预报时刻的特征方向以及权重（见表2）。

表2 各预报时刻特征方向及权重Tab.2 Feature directions and weights for each forecast moment

图6 各预报时刻的相似路径及相似位置Fig.6 Similar tracks and similar positions for each forecast time

形成各预报时刻预报位置以及集合成员位置后（见图7），采用2.1节中的方法形成48条（Kt1·Kt2·Kt3=4×3×4）集合路径，并计算出每一条成员路径的权重。

图7 各预报时刻特征位置及实况路径（0～72 h）Fig.7 Location of features at each forecast time and the real track(0～72 h)

考虑到模式冷启动计算的不稳定性，FVCOM模式和SWAN 模式都提前7 d 开始运行。计算完成后，验证部分站点有效波高的分布情况，给出各预报时刻下部分重点近岸站点以及整个研究区域有效波高的概率分布预测作为参考。

3.1 后报检验与数值预报实验

从与H1、H2 站点实测波高数据的对比结果来看（见图8），后报检验（实况路径）的效果较好，自起始预报时刻～72 h 预报时段，H1 和H2 站点波高的均方根误差平均值分别为0.31和0.60。使用集合路径进行数值预报实验的结果符合实际波高的基本发展趋势，在台风离岸较近时间段内（24～72 h），集合结果基本能够覆盖实际波高的波动范围，且在一定程度上能够反映出波高的极值。值得注意的是，H1 站点位于杭州湾靠内区域，在24～48 h 时段内，集合预报实验结果都小于实际值（见图8a），且两种概率路径（70%和95%）区间值的差异性小于H2 站点，而在48～72 h 时段内，波高集合的结果基本覆盖了实际值；H2 站点位于长江口外，集合预报实验在24 h 附近的波高预测偏大（见图8b），在48 h 附近的波高集合结果基本覆盖了实际值，在此时段内，95%概率路径结果的覆盖效果更好。

图8 H1（a）与H2（b）站点有效波高的验证Fig.8 Validation of effective wave height at H1(a)and H2(b)sites

为量化集合预报在两特征概率情形下各预报时段内对实际有效波高的预测情况，我们对每个实测站位计算其覆盖指数，将指定时段内处于集合预报预测区间范围内的实测波高点数量占该时段内总实测波高点数量的比例作为指定时段的覆盖指数。计算公式为：

式中：CR 表示覆盖指数；Oin表示预报时段内观测点含于集合预报实验结果上下限范围内；Oall为预报时段内所有观测点的数量。作为集合预报结果覆盖范围的参考，CR 的波动范围为0～100%，其值越大越好。CR 的计算结果见表3。从表中可以看出，在0～24 h，台风离岸较远，集合成员离业务预报位置较近，集合预报实验的准确性非常依赖业务预报的准确性，集合结果的发散程度小，因此，该时段由于业务预报效果较差，导致两站的概率路径的CR都为0，其中H1站位集合结果区间范围小；在24～48 h预报时段的CR 较高，两概率路径下的CR 分别为39%和52%，其中H2 站位的波高峰值结果偏大，两概率路径的CR分别为22%和54%；在48～72 h预报时段的两概率路径的CR有所上升，分别为74%和94%。

表3 各站位集合预报的覆盖指数（单位：%）Tab.3 Cover rate of ensemble forecasts at each station(unit：%)

3.2 重要沿岸站点各预报时刻有效波高预测的概率分布参考

在实测站点数据验证良好的基础上，分析重要沿岸以及近岸浮标站点各预报时刻实况路径有效波高的分布情况，并给出集合预报路径相对实况路径有效波高预测的概率分布参考，结果见图9，其中站点自江苏沿岸开始，经沿北支、南支到杭州湾北岸区域，外加部分口外站点（见图1）。近岸站点波高预测概率分布在70%和95%概率路径下表现出的特征有所不同。95%路径在24 h 预报时刻显示出较大的概率分布（见图9b），这一现象在沿岸及近岸站点如H5、堡镇、赵家沟、芦潮港有所体现，24 h台风距离站点较远，95%路径的部分集合成员距离上述站点更近；而70%路径在48 h 预报时刻显示出较大的概率分布（见图9a），该现象在H2、绿华山等部分口外站点有所体现，48 h 台风位于站点附近，70%路径的部分集合成员距离上述站点更近；两种概率路径在72 h预报时刻给出的概率分布情况基本一致。从空间分布来看，H2、嵊山等离岸稍远的站点在整个时间序列上表现出相对较大的波高预测值，而徐六径、崇西、堡镇等站点波高预测相对较小。

图9 重要沿岸站点有效波高预测概率分布Fig.9 Predicted probability distribution of significant wave height at important coastal stations

3.3 有效波高空间预测的概率分布参考

为分析研究区域空间上有效波高的变化情况，图10 和图11 分别给出70%与95%概率路径集合风场情形下各预报时刻有效波高的空间分布情况。图中为分别计算所有集合成员有效波高并减去实况后报的结果，并根据权重加权分别累计差值大于0和小于0的增减量分布。

图10 70%集合风场情形各预报时刻有效波高空间分布Fig.10 The spatial distribution of significant wave height at each forecast moment for the 70%ensemble wind field case

图11 95%集合风场情形各预报时刻有效波高空间分布Fig.11 The spatial distribution of significant wave height at each forecast moment for the 95%ensemble wind field case

在台风中心移向研究区域的时段，与实况后报相比，集合预报的波高结果在台风中心附近的远岸区域表现为正增量，而在近岸处表现为负增量（见图10b、10c、11b、11c），这表明集合预报对波高空间场有重塑作用，导致空间上波高局地分布的差异。两概率路径情形的正增量最高都超过了4 m（见图10b、11b），两者负增量在近岸区域最大为-1 m 左右（见图10c、11c）。当台风中心渐远时，较大概率路径情形的正增量区域（见图11e）相对较小概率路径情形（见图10e）偏移近岸研究区域更远，最大值均为1 m 左右，而较小概率路径情形的负增量区域总体相对偏移近岸研究区域更远（见图10f、11f）。当台风中心离研究区域更远时，两者都不再有显著的正增量贡献（见图10h、11h），且负增量空间分布的差异不大（见图10i、11i）。从整个时序来看，空间中波高正增量区域基本和台风中心位置保持一致，负增量区域滞后于台风中心位置，该特性使得集合预报对部分口外站位波高峰值后的下降时段刻画更好（见图8b，24～72 h）。对比两种路径概率情形可以看出，范围偏小的概率路径（70%）集合预报波高加权结果的空间分布局地差异更大，波高随空间分布的梯度更大。

4 结论

针对近岸区域波高受到很多因素的影响导致的预测难度较大的问题，本文提出了一种新型的台风集合路径的方法。基于历史台风数据和业务中心实时预报数据，采用统计学的方法给出任意实时情况下台风各预测移动方向的概率，使用动态权重生成台风预报路径集合并进行数值预报实验，得出实验结果的概率分布。以长江口及杭州湾附近为研究区域，结合数值模式对台风过境期间近岸区域的有效波高进行后报验证，在验证基础上应用该新型集合路径生成方法进行预报试验，其中在48～72 h时段内，近岸区域H1 浮标集合预报结果的覆盖指数能达到52%左右，H2浮标的覆盖指数能达到94%左右。该方法可适用于台风期间风暴潮、台风浪等要素的预报，可为近岸海域的复杂动力要素预报工作提供新的集合预报技术。