一种热带气旋路径相似性快速计算方法研究
2023-09-16王青颜韩屹刘丽
王青颜，韩屹，刘丽
（1.海南省海洋监测预报中心，海南海口570206；2.国家海洋环境预报中心，北京100081；3.航天宏图信息技术股份有限公司，北京100195）
0 引言
热带气旋是一个天气尺度的非锋面低压系统，在温暖水域上发展，有组织的对流，最大平均风速超过17.2 m/s且其区域在中心附近延伸超过一半并持续至少6 h。热带气旋路径被认为是由多种因素综合影响作用的结果。两条路径相似的热带气旋在一定程度上能反映热带气旋的环境条件和演变过程的相似，也反应了诸多影响因子综合作用的等效[1]。
在台风风暴潮预报工作中，会选择一条热带气旋路径作为目标。对历史热带气旋进行筛选，首先得到若干与目标路径相似的候选结果，再逐一分析其与目标热带气旋的强度、形势场的差异，得到与目标热带气旋“相似”的历史热带气旋过程，再考察历史过程中验潮站的增水情况，作为开展当前台风风暴潮预报的重要参考依据。因此，通过路径相似寻找与当前热带气旋相近的历史过程是开展台风风暴潮预报的必要基础工作之一。
在已有热带气旋路径相似性的研究中，研究人员提出了多种相似性度量方法。例如，刘勇等[2]通过引入相似离度方法对热带气旋的相似程度进行评估；孔令娜[3]提出了一种热带气旋路径相似程度的度量方法，即通过豪斯多夫距离法进行判别；郑霞[4]应用主成分分析方法确定了台风各影响因子的权重，再基于加权的动态时间规整法开展热带气旋相似程度评估；程源清[5]使用弗雷歇距离法及改进方法来计算热带气旋的相似性；CHEN等[6]和DI等[7]以动态时间规整法为基础开展了热带气旋相似性计算和分析。此外，相关研究也提出了不同的热带气旋路径相似性的计算方法及实现步骤[8-10]。
目前，热带气旋路径相似计算方法主要可分为两种：一是基于地理信息系统（Geographic Information System，GIS）的空间分析功能，通过以样本热带气旋做缓冲区，分析其与目标热带气旋路径的空间关系，进而确定其相似程度[11-12]；二是计算样本热带气旋与目标热带气旋路径之间的差异性指标（通常使用距离指标衡量），如欧氏距离、弗雷歇距离、豪斯多夫距离等。
上述方法也存在一定的不足，如在GIS 空间分析方法中缓冲区大小与空间分析结果强相关，如何科学合理地确定缓冲区半径尚无较为明确的方法；在通过距离法判别热带气旋路径相似的方法中，欧式距离法需要两个路径中离散点数量相同，豪斯多夫距离法可理解为取样本热带气旋到目标热带气旋路径点集的最短距离的最大值，如路径点集中出现极端点，将导致该距离出现较大偏差；同样地，弗雷歇距离法对于热带气旋路径点集的位置也比较敏感，较为相近的路径会因为点集位置不同得出差异较大的计算结果。此外，弗雷歇距离法、动态时间规整法的计算耗时较大。
本文提出了一种热带气旋路径相似的快速计算方法，通过评估目标热带气旋与样本热带气旋路径的相近度、重叠度和连续度，确定两者的相似程度。为验证算法的有效性，设计了对比试验，结果表明与现有算法相比，本方法在保持较好相似性匹配的基础上，计算效率较高。
1 数据源及数据处理
数据源为中国气象局（China Meteorological Administration，CMA）热带气旋最佳路径数据集[13-14]。该数据集提供1949 年以来西北太平洋（赤道以北，东经180°以西）海域热带气旋每6 h 的位置和强度信息，包括时间（世界时，下同）、强度、纬度、经度、中心最低气压等。
本文使用的CMA 热带气旋最佳路径数据的时间范围是1949—2020 年。为后续开展热带气旋路径相似性计算分析，我们对数据集作如下处理：
①删除热带气旋副中心情况，形成包括热带气旋序号、时间、经度、纬度、强度、中心最低气压的路径点数据集；
②按热带气旋序号逐条遍历路径点数据，当一条路径点中出现前后相邻点空间位置重合时，做删除处理，相同位置保留一个点，避免存在路径点重叠的情况；
③将经上步处理后的热带气旋路径点集合数据按热带气旋序号排序，形成路径（点）集合，用于后续计算分析。
在实际工作中，一般是在某个区域内进行热带气旋路径相似匹配评估，如针对热带气旋72 h 预报集合路径覆盖范围或者预报员指定的海域等。在这种情况下，需额外进行数据处理，即对每条热带气旋路径查找满足空间范围条件的路径点并单独存储。特殊情况是一条热带气旋路径中在某个时间段的路径点满足空间范围约束条件，但之后的某段时间的某个或者某几个路径点不满足，之后的时间段的路径点又满足条件，在该情况下，需获取第一个和最后一个满足空间约束的路径点，然后依次取出中间路径点（包括中间不满足需求的路径点），确保不出现在时间序列上割裂热带气旋路径的情况。经上述处理后，生成满足空间范围条件的热带气旋路径点集合。
2 热带气旋路径相似算法
本文提出了一种热带气旋路径相似的评估算法。通过构建轨迹点集合的椭圆序列，设计以相近度、重叠度和连续度作为指标，综合评价两条热带气旋路径的相似程度。
2.1 定义
分析热带气旋路径相似时，先选定一条热带气旋路径作为目标热带气旋路径，定义为Cp，然后在历史热带气旋路径集合中逐条计算相近度，选取的热带气旋路径称为样本热带气旋路径，定义为Cq。
目标热带气旋路径由若干个路径点组成，故Cp=
以目标热带气旋路径点集合中顺序排列的两个路径点为焦点，可生成一个椭圆，将逐个路径点依次处理，可生成一系列的椭圆即椭圆序列，为目标热带气旋路径椭圆序列，定义为Ep=
为了便于理解上述定义，下图给出了一个由4个路径点构成的热点气旋路径
图1 热带气旋路径点及椭圆序列定义示意图Fig.1 The diagram of tropical cyclone track points and elliptical sequence
2.2 相近度
相近度可以考察目标热带气旋路径和样本热带气旋路径在地理空间上的接近程度。通过计算Cp中的路径点落入Sq中椭圆的数量以及Cq中的路径点落入Cp中椭圆的数量，经归一化处理后得到相近度。具体的计算步骤如下：
①椭圆参数设置
椭圆在其焦点确定的情况下，其长半轴和短半轴参数直接决定了椭圆的形状，推荐椭圆长半轴a应用如下公式：
式中：K为尺度参数，取值为；Deuc为椭圆焦点的距离，采用欧氏距离。长半轴确定后，短半轴b计算方法为：
式中：a为长半轴；e为椭圆焦距。
②路径点与对应椭圆序列的匹配情况计算
分别计算Cp与Sq、Cq与Ep的匹配情况。以前者为例，匹配程度用m(Pi,Sq)表达，计算方法如下：
情况1：当Pi与Sq中任意椭圆处于包含关系时；
情况2：除去情况1的其他情况。
类似地，用同样方法计算后者的M( )Qi,Ep。
③相近度计算
目标热带气旋路径与样本热带气旋路径的相近度M(Cp,Cq)的计算公式如下：
式中：M(Cp,Cq)的取值范围为0～1，越接近1表明两条路径更相近。
2.3 重叠度
当目标热带气旋和样本热带气旋的路径点均处于对应椭圆中时，虽然其空间位置不同，但相近度的评分是相同的，这种情况下仅依靠相近度无法区分。为解决这种情况，设计了通过目标热带气旋与样本热带气旋路径椭圆序列重叠程序的指标来区分其相似程度。显然地，重叠程度高的椭圆序列所对应的两条路径相似性更高。
本文通过路径点与椭圆的归一化距离（以目标热带气旋路径点与对应的样本热带气旋椭圆为例，反之则是样本热带气旋路径点与对应的目标热带气旋椭圆），作为评价椭圆序列重叠情况的指标。
①目标热带气旋路径点与样本热带气旋椭圆的归一化距离计算
设Pi为目标热带气旋路径中的点，Pi∈Cp，Si为样本热带气旋路径椭圆序列Sq中的椭圆，Si∈Sq。计算公式为：
式 中：K为Si椭 圆 的 焦 点Sf1、Sf2或 椭 圆 原 点Sfc；deuc(P,K)为点Pi到点K的欧式距离。
Α取值公式为：
式中：b为椭圆短半轴；a为椭圆长半轴。Β取值公式为：
式中：e为所对应的椭圆焦距。
②样本热带气旋路径点与目标热带气旋椭圆的归一化距离计算
设Qi为样本热带气旋路径中的点，Qi∈Cq，Ei为目标热带气旋路径椭圆序列Ep中的椭圆，Ei∈Ep。dr(Q,E)的计算公式同dr(P,S)。
③重叠度计算
设目标热带气旋路径与样本热带气旋路径的共享度为O(Cp,Cq)。计算公式为：
式中：O(Cp,Cq)的取值范围在0～1，越接近1表明两条路径更为接近。
2.4 连续度
将Cp与Sq匹配的顺序序列定义为Seque(Pi,Sq)，如Seque(P1,Sq)=，将首次匹配位置定义为f(Pi,Sq)。计算方法为：
①若i=1 且Seque(Pi,Sq)不为空的情况下，f(P1,Sq)= minSeque(P1,Sq)；
②若i＞1 且Seque(Pi,Sq)不为空的情况下，f(Pi,Sq) =min[Seque(Pi,Sq)≥f(Pi-1,Sq)] ；
③其他情况下，f(Pi,Sq) =-1。
目标热带气旋路径首匹配位置序列Sequep用
连续度计算公式为：
式中：C(Cp,Cq)的取值范围在0～1，越接近1表明两条路径更为连续。grade函数计算规则如下:
①当i=1 且Sequep[i] 的 值≠-1 时，grade(Sequep[i] )=1；
②当i＞1 且Sequep[i] 的值≠-1 且Sequep[i] ≥Sequep[i] - 1 时，grade(Sequep[i] )= 1；
③其他情况下，grade(Sequep[i] )= 0。
2.5 综合评分
综合评分采用相近度、重叠度和连续度的均值，计算公式如下：
式中:M(Cp,Cq)、O(Cp,Cq)、C(Cp,Cq)分别表示评分之前的相近度、重叠度和连续度。
3 试验
为验证本文方法的有效性，选取影响南中国海区域的若干历史热带气旋路径作为目标热带气旋路径，分别使用本文算法与前人研究中效果较好的弗雷歇距离法和动态时间规整法开展相似性计算，同时比较不同算法的计算时间。
3.1 试验环境
本文试验选用的硬件是联想P43s笔记本电脑，硬件配置为中央处理器i7-8565U、16G 内存、1T 硬盘，操作系统为64 位Windows 10 操作系统，本文算法的开发语言是Python，版本为3.9，开发及调试使用的软件为PyCharm社区版。
3.2 形态相似结果
选择南中国海作为研究区域，选择该区域内西移路径和特殊路径的8 个热带气旋为试验对象（即目标热带气旋路径），其信息和空间范围条件见表1。
表1 样本热带气旋及空间范围Tab.1 Experimental tropical cyclones and spatial extent
图2—9 绘制了使用不同方法（本文方法、动态时间规整法、弗雷歇距离法）计算的目标热带气旋路径在热带气旋路径数据集中的路径相似的结果。
图2 热带气旋“康森”（201002）路径相似计算结果对比Fig.2 Comparison of similar paths of tropical cyclone"Conson"（201002）searched by different methods
图3 热带气旋“纳沙”（201121）路径相似计算结果对比Fig.3 Comparison of similar paths of tropical cyclone"Nesat"（201121）searched by different methods
图4 热带气旋“尤特”（201311）路径相似计算结果对比Fig.4 Comparison of similar paths of tropical cyclone"Utor"（201311）searched by different methods
图5 热带气旋“威马逊”（201410）路径相似计算结果对比Fig.5 Comparison of similar paths of tropical cyclone"Rammasun"（201410）searched by different methods
图6 热带气旋“达维”（201726）路径相似计算结果对比Fig.6 Comparison of similar paths of tropical cyclone"Damrey"（201726）searched by different methods
从试验结果看，3 种方法均能够实现相似热带气旋路径的分析计算，并取得较好的结果。本文方法与其他两种方法的结果具有较高的一致性，但根据目标热带气旋路径的不同，表现有所区别。具体为：
①在西移路径中，热带气旋在菲律宾以东洋面生成，经过南海中北部海域后向偏西方向移动，在粤西或海南岛登陆（如热带气旋“康森”、“纳沙”、“尤特”、“威马逊”、“山竹”）或经过南海中南部海域在越南登陆（如热带气旋“达维”）。在指定空间范围中，本文方法比弗雷歇距离法给出的相似路径（图中绘制的排序前4的路径）与样本路径更为近似（见图2、3、4、5、6、8），主要体现在筛选出的相似路径与目标路径在整体空间形态上一致，且筛选出的相似路径上的路径点与目标路径点的距离更近。本文方法与动态时间规整法给出的相似路径结果也基本一致，在空间形态上表现较好。
以热带气旋“康森”为例，经3 种方法计算后的路径相似排序前10的结果见表2。本文方法与动态时间规整法的结果有8条相同，与弗雷歇距离法有6条相同，证明路径相似性计算结果具有较高的一致性。
表2 3种方法下热带气旋“康森”（201002）路径相似性计算结果排序表Tab.2 Similar paths of tropical cyclone"Conson"(201002)searched by three methods
②出现特殊路径时，选择了两种情况来具体考察3种方法的计算结果。
一是热带气旋“艾云尼”，其在南沙西北部海域生成后，按偏北路径移动并在广东徐闻登陆后南下，在海南海口二次登陆后进入粤西海域，并在广东阳江三次登陆。从图7可以看出本文方法的结果与动态时间规整法、弗雷歇距离法的结果不同，3种方法都给出了200802 号和197233 号热带气旋是最为相近的路径（在本文方法和动态时间规整法的路径相似性结果排序中均为1、2，在弗雷歇距离法中排序为1、4）。图7 中本文方法计算出的其余路径（198611 号、196426 号）与其余两种方法差别较大，其原因是本文算法是建立在路径点的椭圆序列匹配基础上的，与其余两种算法基于距离的原理不同。本文方法路径匹配的相似程度由椭圆序列的匹配程度决定，当缺少整体匹配情况时，局部匹配度高的结果会排序靠前。此外，动态时间规整法结果（197523 号、195040 号）、弗雷歇距离法结果（196535 号、197523 号）也表明其与目标路径在空间上的差别很大。经3 种方法计算后的热带气旋“艾云尼”路径相似性排序前5的结果见表3。
表3 3种方法下热带气旋“艾云尼”（201805）路径相似性计算结果排序表Tab.3 Similar paths of tropical cyclone"Ewiniar"（201805）searched by three methods
图7 热带气旋“艾云尼”（201805）路径相似计算结果对比Fig.7 Comparison of similar paths of tropical cyclone"Ewiniar"（201805）searched by different methods
图8 热带气旋“山竹”（201826）路径相似计算结果对比Fig.8 Comparison of similar paths of tropical cyclone"Mangkhut"（201826）searched by different methods
二是热带气旋“娜基莉”，其在南海中部生成，先向东移动再转向掉头向西，经海南中南部海域后在越南登陆。3 种方法都没有匹配到与目标路径整体相似的结果，动态时间规整法和弗雷歇距离法给出的排序第一的结果为199034 号，从图9 可以看出该路径与目标路径在空间形态上匹配度并不高。主要原因是上述两种方法的计算原理是通过构建序列点集的距离矩阵，寻找从左下—右上的最短路径，并以这条路径值作为相似性评分依据，当目标路径出现转向移动，而样本路径正好位于其转向的中间区域时，目标路径点集和样本路径点集的距离虽然在形态不一致，但是两者距离较小，就会出现相似性评分较高情况。本文方法给出的排序第一的结果为201726 号，虽然整体上与目标路径不一致，但与转向后的目标路径相似性较高。出现这种情况的原因与热带气旋“艾云尼”类似，在整体匹配不到的情况下，给出了局部匹配度最好的结果。经3 种方法计算后的热带气旋“娜基莉”路径相似性序前5的结果见表4。
表4 3种方法下热带气旋“娜基莉”（201928）路径相似性计算结果排序表Tab.4 Similar paths of tropical cyclone "Nakri"（201928）searched by three methods
图9 热带气旋“娜基莉”（201928）路径相似计算结果对比Fig.9 Comparison of similar paths of tropical cyclone“Nakri”(201928)searched by different methods
综合以上两种情况，在出现热带气旋特殊路径时，本文方法可以在缺少整体相似的情况下给出局部匹配度最好的结果，这种处理方法是较为合理的。
3.3 计算时间对比
分别使用本文算法、动态时间规整法和弗雷歇距离法对样本热带气旋路径开展相似性计算，计算结果见表5。需要说明的是，热带气旋路径条数为在空间范围条件限制下参与相似性计算的路径条数，3 种方法的计算时间均为调试模式的计时结果。
表5 样本热带气旋相似性计算时间Tab.5 The duration time of similar paths of experimental tropical cyclone computed by three methods
从表5 分析可知，在热带气旋路径相似的分析中，使用相同的路径数据集，本文方法计算的所需时间最短（平均约为0.49 s）、弗雷歇距离法次之（平均约为0.85 s），动态时间规整法耗时最长（平均约为1.76 s）。相较于其他两种方法，本文算法在计算效率上优势明显，其原因是本文方法在评价指标计算时均可直接进行浮点运算，不涉及递归处理，故计算速度快，耗时低。
4 结论及讨论
4.1 结论
本文针对寻找相似热带气旋路径的问题，提出了一种快速计算方法并设计了试验方案对其有效性进行验证。选取南中国海区域的8个热带气旋作为样本，通过与动态时间规整法和弗雷歇距离法的对比，证明本文方法是可行的。本文方法具备如下特点：
①对区域内常规的经菲律宾东部进入南中国海向西或西北移动的热带气旋路径，本文方法的相似计算结果与动态时间规整法具有较高的一致性，排序靠前的相似路径与目标路径在空间形态上相似度较高，两种算法均优于弗雷歇距离法。
②对区域内热带气旋特殊路径的情况，本文方法的结果排序是优化给出整体空间形态接近的结果，次之给出局部相似的结果。在存在整体相似结果的情况下，本文算法与动态时间规整法具有较好的一致性。
③本文方法在计算效率上比其他两种方法具有较明显的优势，在相同条件下，计算耗时仅为动态时间规整法的27%、弗雷歇距离法的57%。
综上，本文提出的方法在热带气旋路径相似匹配上具有较好效果，与现有的动态时间规整法、弗雷歇距离法得到的结果的一致性高，针对特殊路径可以给出局部相似度高的结果作为参考，且本文方法的计算效率高。
4.2 讨论
从试验结果来看，本文方法在筛选热带气旋主要路径的样本时有较好的效果，但在处理特殊路径的热带气旋时，可能会出现由于无法匹配到整体形态相似路径进而给出局部路径相似的结果。在这种情况下，在设计算法时如何定量化的定义“相似”值得进一步讨论，是倾向于局部相似，还是寻找空间距离远但整体形态更接近的样本，这就是下一步的工作方向。此外，路径相似只是热带气旋相似判定的第一步，在路径相似的基础上，还需要对热带气旋的强度、背景场等开展综合性分析，才能判断热带气旋是否“相似”。在下一步的工作中，除了开展针对特殊路径情况下相似计算的算法优化外，还将考虑加入热带气旋强度变化以及相似背景场的相关指标。
