引言

雷暴常伴有短时强降水、破坏性大风和冰雹等灾害性天气，其每年造成的国民经济损失和人员伤亡较重，一直是党和政府关注的焦点之一^[1-3]。因此，加强对雷暴天气的分析和研究，提高雷暴监测、预警和防御能力，对防灾减灾和工农业服务有重要意义。国内许多学者展开了雷暴及其预警方法的研究，俞小鼎等^[4]系统阐述了雷暴的临近预报方法，指出以主观预报为主，结合客观算法，更有利于促进雷暴预报的进一步发展和提高。郝莹等^[5]对比分析了指标叠加法和判别法制作的雷暴潜势预报结果，指出指标叠加法的预报效果较好，并验证了对流有效位能对雷暴发生有较好的指示作用。胡晓等^[6]选取对流参数作为预报因子，在分析并确定阈值的基础上运用权重法建立了雷暴潜势预报模型，检验表明其对雷暴有较好的预报作用。王秀明等^[7]澄清了雷暴发生三要素容易混淆的概念，阐明了雷暴预报中的几个基本问题。田琨等^[8]基于贝叶斯分类法和Logistic回归分析法建立雷暴预报模型，模型预报结果表明Logistic回归分析法对雷暴具有一定的预报能力。上述研究对提高雷暴预报效果有很大帮助，但不足之处在于因子的选取和阈值的确定具有一定的主观性，整体而言识别率不高或识别率高的方法复杂、计算量大。

由于传统雷暴预报方法的主观性和信息不完整性，业务应用上常呈现预报准确率不高、预警时效短、有效防御措施低等弊端。作为自驱动的技术，投影寻踪聚类^[9-10]方法在上世纪70年代被提出，现已在水文、气象上得到广泛应用并取得了较好效果。倪长健和崔鹏^[11]综合利用动态聚类方法和投影寻踪思想，构造投影指标并建立了投影寻踪动态聚类新模型。模型具有计算过程中不需要人为给定参数和聚类结果客观、明确的优点, 同时又具有稳定性好、操作简便等特点^[12-13]。本文从达州市2015年10月—2016年9月雷暴个例的层结条件、水汽条件出发构建指标体系，以850 hPa与500 hPa温差、大气可降水量、K指数、对流有效位能、对流抑制和抬升指数作为雷暴预警因子，构建雷暴预警投影寻踪动态聚类模型，以期为雷暴潜势预报提供一种新思路，进而为雷暴预警提供更加科学的参考依据。

1 资料与方法 1.1 资料说明

使用的资料为2015年10月—2016年9月达州市闪电定位和探空观测资料，闪电定位资料时间间隔为1 h，探空观测资料时间间隔为12 h，个例探空资料统一选取08时(北京时，下同)。闪电现象是雷暴处于发展旺盛成熟状态的明确指标^[14]，云中的放电强度和频次取决于雷暴云垂直伸展高度及强度。利用闪电定位资料来确定达州站(107.50°E，31.2°N)周围半径50 km范围内的雷暴活动^[15]。若观测中无地闪记录则认为无雷暴发生^[8]。随机选取时间段内的50个样本作为雷暴预警研究样本，其中发生雷暴的样本26个，未发生雷暴的样本24个。

1.2 投影寻踪动态聚类模型简介

为最大程度地挖掘数据信息和反映数据特征，从不同角度分析数据用以查找最优投影方向，这就是投影。将高维数据信息投影到低维空间，有利于常规方法处理分析和直观形象分析数据结果。设第i个样本第j个指标为x_ij⁰ (i = 1, …, n；j = 1, …, m样本容量用n表示，选取的指标数目用m表示)，可按照下述步骤构建投影动态聚类模型：

第一步，数据无量纲化，目的在于使各评价指标的量纲差异尽可能的消除，对于越大越优的指标，有

$ {x_{ij}} = \frac{{x_{ij}^0-x_{j\min }^0}}{{x_{j\max }^0-x_{j\min }^0}} $

(1)

若指标是越小越优，则有

$ {x_{ij}} = \frac{{x_{j\max }^0-x_{ij}^0}}{{x_{j\max }^0-x_{j\min }^0}} $

(2)

其中，x_jmin⁰和x_jmax⁰分别为第j个指标的样本最小值和最大值。

第二步，线性投影，即在线性空间实现高维数据的投影。设$\vec a$为m维投影方向向量，其分量为a₁，a₂，…，a_m，分别对应于m个指标的权重，并有$\sum\limits_{j = 1}^m {{a_j} = 1} $，则x_ij投影特征值z_i可表示为

$ {z_i} = \sum\limits_{j = 1}^m {{a_j}{x_{ij}}} \left( {i = 1, 2 \cdots, n} \right) $

(3)

第三步，构建投影指标，这是把高维数据投影聚类到低维空间遵循的关键所在，也有利于寻找最优投影方向。全部样本的投影特征值序列所构成的集合表示为Ω={z₁，z₂，...，z_n}，以动态聚类法为依托，将其聚为p (p ≤ n)类，即：

（1）随机选取p个点，并作为p个聚核，记为L⁰ =(A₁⁰，A₂⁰，...，A_P⁰)；

（2）以L⁰把Ω中的点分为p类，记做Θ⁰ = (Θ₁⁰，Θ₁⁰，...，Θ_p⁰)。其中Θ_i⁰ ={z∈Ω|d(A_i⁰ - z) ≤ d(A_j⁰ - z)，∀j = 1，2，...，p，j≠i}，d(A_i⁰ - z)为点A_i⁰和集合Ω中任一点的绝对值距离。

（3）由Θ出发，对新的聚核L¹，L¹ =(A₁¹，A₂¹，...，A_P¹)加以计算。其中，$A_i^1 = \frac{1}{{{n_i}}}\sum\limits_{{z_i} \in \Theta _i^0}^p {{z_i}} $，类Θ_i⁰中有n_i个点。

（4）重复以上步骤，由此得到一个分类结果序列${ \vee ^k} = \left( {{L^k}, {\Theta ^k}} \right), k = 1, 2, \cdots $。$D\left( {A_i^k, \Theta _i^k} \right) = \sum\limits_{{z_i} \in \Theta _i^k} {\left| {{z_i}-A_i^k} \right|} $，${U_k} = \sum\limits_{i = 1}^p {D\left( {A_i^k, \Theta _i^k} \right)} $，若$\frac{{\left| {{U_{k + 1}}-{U_k}} \right|}}{{{U_{k + 1}}}} \le \varepsilon $，则终止算法。其中，允许误差范围的充分小量用ε表示。这种算法理论证明是收敛的。

属于第h类的所有样本投影特征值构成的集合用Θ_h(h = 1，2，...，p)表示，定义所有两投影特征值间的绝对值距离为d(z_i，z_j)。$dd\left( a \right) = \sum\limits_{h = 1}^p {{D_h}\left( a \right)} $表示类内样本的临近程度，其中${D_h}\left( a \right) = \sum\limits_{{z_i}, {z_j} \in {\Theta _n}} {d\left( {{z_i}, {z_j}} \right)} $为类内聚集度，dd(a)越小，表示类内样本越高的聚集程度。另记投影分散度为$ss\left( a \right) = \sum\limits_{{z_i}, {z_j} \in \Omega } {d\left( {{z_i}, {z_j}} \right)} $，它表示所有样本序列投影特征值的离散程度，所有样本序列投影特征值随着ss(a)愈大，其离散程度则意味着越高。

以ss(a)和dd(a)的定义为依据，投影寻踪动态聚类模型的投影指标QQ(a)为

$ QQ\left( a \right) = ss\left( a \right)-dd\left( a \right) $

(4)

由此可见，投影指标的值随着投影分散的程度越大或者类内聚集的程度越小而越大。如果QQ(a)值取最大时，类间样本尽量散开和类内样本尽量集中的聚类要求就达到了。

第四步，模型求解，其关键在于寻找满足达到最大值时的QQ(a)对应的最优投影方向向量$\overrightarrow {{a^*}} $，所以，下面优化相关问题能够通过投影寻踪动态聚类模型加以描述

$ \begin{array}{l} {\rm{Min}}\;\;\;\;QQ\left( {\vec a} \right)\\ s.t.\;\;\;\;\sum\limits_{j = 1}^m {{a_j} = 1} \end{array} $

(5)

其中，s.t.表示受约束，本文应用免疫进化算法^[16]求解上述优化问题。

2 指标选取

基于雷暴发生的要素，探空资料850 hPa与500 hPa温差（ΔT_850-500）、大气可降水量（PW）、K指数（KI）、对流有效位能（CAPE）、对流抑制（CIN）和抬升指数(LI)共6个因子对雷暴潜势有较好的指示意义, 其中CAPE经过虚温订正。ΔT_850-500表示大气层结静力稳定度。PW是单位面积上整层大气中所有水汽全部凝结并降落到地面的降水量。K反映大气层结稳定情况, 其值越大表示大气层结越不稳定。CAPE为大气浮力对流有效位能的大小。CIN反映处于大气底部的气块若要能自由地参与对流至少要从其他途径获得的能量下限，雷暴的发生往往要求CIN有一个较为合适的值。LI指数表示气块先干绝热上升然后湿绝热线上升到500 hPa处的温度与该处实际大气温度的差值，差值的绝对值越大预示着出现对流天气的可能性越大。

3 应用实例

雷暴是伴有雷击和闪电的局地对流天气，它产生于强烈的积雨云中，常伴有强烈阵雨或暴雨，有时伴有冰雹和龙卷风。雷暴预警的实质就是如何把各雷暴样本的多维分类指标综合成一维或二维指标，然后进行聚类分析，由于该问题涉及到层结条件、水汽条件等诸多不确定因素，而各单因素指标的结果往往是不相容的，现有方法对此类问题的处理在客观性、可操作性方面尚存在一定的局限。选取2015年10月— 2016年9月的50个样本作为研究实例，雷暴分类指标体系由ΔT_850—500、PW、KI、CAPE、CIN和LI共6个指标组成，将其聚为2类，按投影寻踪动态聚类模型的算法流程，其中m=6，n=50，p=2，得到投影方向向量为$\overrightarrow {{a^*}} $=（0.125，0.046 7，0.208 3，0.295，0.141 7，0.183 3）。

表 1（见上页）给出川东北雷暴样本的预报因子及对应的特征值，根据投影特征值z^*(i)的大小和雷暴样本的聚类结果分析可知，(1)其投影特征值在[0.506 8，1]范围，可判为有雷暴发生；若其特征值在[0，0.461 0]范围为第Ⅱ类，其值小于等于0.461 0，可判为无雷暴发生。(2)对模型预报结果进行定量评估^[17]，命中率(POD)为92.31%，漏报率(MR)7.69%, 虚警率(FAR)23.40%，临界成功指数(CSI)为72.00%, 表明该模型对雷暴有一定的预报能力。(3) 6—9月是川东北雷暴的高发期，雷暴潜势较大，一般来说ΔT_850-500越大、KI越大、CAPE越大、LI负值越小、一定的PW和CIN，雷暴发生的潜势越大。6个指标中CAPE所占比重最大，表明雷暴发生时对流有效位能的重要性；其次是KI和LI，它反映了大气层结稳定度和出现对流天气的可能性；然后是ΔT_850—500和CIN，它反映了大气上下层温差和对流抑制；PW占比重最小，但反映了雷暴发生需要一定的水汽条件。

图 1给出雷暴样本投影特征值的散点图，从中可见，样本实测中有雷暴发生26个，没有雷暴发生24个，样本容量和选取较为合理。模型计算结果有无两类的分隔区间为(0.461 0，0.506 8), 靠近分隔区间范围的点相对较少，表明模型有较好的类内聚合性和类间分离性。当模型预报结果和实测在参考线同侧时，预报正确，否则为错误；当实测出现而模型未预报出，则为漏报，当实测未出现而模型预报出，则为空报。通过分析得出，模型虚警率相比同类预报结果偏高，命中率和漏报率优于同类模型结果，临界成功指数与同类模型预报结果相当。

综上所述，投影寻踪聚类模型在川东北雷暴潜势预报的应用表明，该方法所得结果明确、客观，具有良好的类内聚合性和类间分离性，能合理地反映雷暴潜势的层结条件、水汽条件等特征，为雷暴预警提供相应的理论依据。

4 结论与讨论

本文利用川东北闪电定位和高空观测资料，从层结条件、水汽条件出发选取850 hPa与500 hPa温差、大气可降水量、K指数、对流有效位能、对流抑制和抬升指数作为雷暴预警因子构建指标体系，使用投影寻踪聚类模型对雷暴进行预警，得出如下结论：

(1) 模型预报结果定量评估，命中率POD为92.31%，虚警率FAR为23.40%，临界成功指数CSI为72.00%，表明该模型对川东北雷暴潜势预报是有效、合理的。

(2) 与常规预警方法相比，该模型具有识别率高，计算简便，客观性强等特点，为雷暴预警提供一种新的客观方法。但由于雷暴预报的复杂性，上述模型有待进一步完善，如增加物理量预报因子，提高CSI评分。另外，模型虽然给出了有、无雷暴的特征值结果，然而在预报和实测中还存在雷暴强度和等级划分，鉴于此，可考虑进行雷暴强度的分级预报模型研究。