2. 中国人民解放军 93011部队气象台, 延吉 133000
2. Meteorological Observatory, Unit No. 93011 of PLA, Yanji 133000
雷暴[1]是一种典型的中小尺度强对流天气系统,常伴有闪电和雷击。雷暴产生于变化剧烈的积雨云中,常常伴有强降水,有时伴有冰雹和龙卷风,具有极强的破坏力和杀伤力,给人们生命和财产安全带来巨大威胁[2]。
目前,国内外采用多种手段和预报方法进行雷暴的预报研究。2008年,Schmeits等[3]通过逻辑回归的方法建立了基于模式输出产品的雷暴预报系统,该系统能够实现12 h雷暴强度和区域预警。2011年,Drira等[4]利用观测资料,采用人工神经网络方法成功实现了印度Paradeep地区和澳大利亚新南威尔士州的雷暴预报,有效地提高了这两个区域的雷暴预报和分类的准确程度。吴蓁等[5]基于构成要素预报方法, 结合常规观测资料、卫星云图和雷达观测资料,对河南地区雷暴天气进行要素分析,可提前几小时进行雷暴短时预报。张强等[6]利用站点观测资料,采用改进遗传算法优化小波神经网络的方法进行雷暴预报,结果表明与BP神经网络等方法相比,该方法预报准确率更高,具有更好的非线性处理能力和泛化性。安洁等[7]利用WRF预报产品和观测资料,采用综合指标叠套的方法,研究不同海域雷暴的影响因子及其预报模型,结果表明,在东海和南海区域试报时准确率可达70%以上,说明综合指标叠套法用于海上雷暴预报是可行的。
以数值预报为基础的天气预报技术路线包含两个方面的内容,一是提升数值预报模式的能力,二是提高以优选预报因子为主要手段的数值预报产品解释应用能力[8]。贝叶斯分类法是气象统计预报与分析的重要方法,是目前公认的一种简单而有效的数据挖掘分类算法[9],贝叶斯判别准则(Bayes Discriminatory Criterion,BDC)通过计算预报样本属于所有类的概率,具有最大概率的类便是预报样本所属的类[10]。采用BDC建立雷暴预报方程,因子选择时受显著性检验标准的影响较大,并且不能剔除那些通过检验但是预报效果不佳的因子。而利用二进制粒子群(Binary Parti⁃ cle Swarm Optimization, BPSO)这种启发式算法,能够从大量待选因子中,筛选出相应适应度函数下的最优因子组合,即适应度函数对应的最优子集,提升BDC的预报效果。因此本文拟采用BPSO-BDC方法探索和研究雷暴预报的最优模型。
1 基本原理 1.1 BDC的基本原理BDC的基本原理[11]是:统计样本类别的历史信息,总结出客观事物分类的规律性,建立判别函数和判别准则,依据总结出的判别函数和判别准则,判断新样本的所属类别。
如果把k个因子的样品与预报量的N个级别对应地分为N组,并把它们看成N个总体中抽取的样本。然后假设由k个因子x1, x2, …, xk组成的每一个体都来自于N个总体A1, A2, …, AN中的一个。现在任给一个体,记为:
$ {X_0} = {\left( {{x_{10}}, {x_{20}}, \cdots , {x_{k0}}} \right)^\prime } $ | (1) |
要判断该个体来自于哪一个总体,即需要按一定规则把X0划归到N个总体中的某一个An(n=1, 2, …, N)。
BDC可以描述为,假设各总体An(n=1, 2, …, N)的概率分布密度函数fn(X)为已知,又知其先验概率为pn, 在错分损失相等的前提下,可建立判别函数:
$ {p_n}{f_n}\left( x \right)\;\;\left( {n = 1, 2, \cdots , N} \right) $ | (2) |
如果:
$ {p_l}{f_l}\left( X \right)\; = \mathop {\max }\limits_{1 \le n \le N} \left[ {{p_n}{f_n}\left( X \right)\;} \right] $ | (3) |
则把个体X0划归到总体中的Al。
在进行因子筛选时,试验采用t检验对所有强对流指数进行显著性检验,选出在有无雷暴两类总体中具有显著性差异的强对流指数,作为初选因子。假设某一强对流指数xi满足二级判别函数显著性检验的条件,其中A、B两类总体的样本容量分别为nA, nB,可建立对应统计量t:
$ t = \sqrt {\frac{{{n_A}{n_B}\left( {{n_A} + {n_B} - 2} \right)}}{{{n_A} + {n_B}}}} \cdot \frac{{\overline {x_l^A} - \overline {x_l^B} }}{{\sqrt {{n_A}{{\left( {s_i^A} \right)}^2} + {n_B}{{\left( {s_i^B} \right)}^2}} }} $ | (4) |
统计量服从自由度为(nA+nB-2)的t分布。其中:
$ \left\{ \begin{array}{l} {\left( {s_i^A} \right)^2} = \frac{1}{{{n_A}}}\sum\limits_{k = 1}^{{n_A}} {{{\left( {x_{ki}^A - \overline {x_l^A} } \right)}^2}} \\ {\left( {s_i^B} \right)^2} = \frac{1}{{{n_B}}}\sum\limits_{k = 1}^{{n_B}} {{{\left( {x_{ki}^B - \overline {x_l^B} } \right)}^2}} \end{array} \right. $ | (5) |
试验发现显著性标准对初选因子个数的影响较大。当满足显著性检验的因子个数较多时,由于判别因子间的相互影响,会导致整体判别效果的下降。这时采用Fisher准则对初选因子进行排序,选取其中量值较大的因子作为预报因子。Fisher准则的表达式为:
$ \lambda = \frac{{{{\left( {\overline {{y^A}} - \overline {{y^B}} } \right)}^2}}}{{\sum\nolimits_{k = 1}^{{n_A}} {{{\left( {y_k^A - \overline {{y^A}} } \right)}^2} + \sum\nolimits_{k = 1}^{{n_B}} {{{\left( {y_k^B - \overline {{y^B}} } \right)}^2}} } }} = \mathit{max} $ | (6) |
其中ykA(k=1, 2, …, nA)表示A类的判别函数值,ykB(k= 1, 2, …, nB)表示B类的判别函数值,nA和nB分别是A和B类的样本数,yA和yB分别为A和B类的数学期望。
1.2 PSO算法粒子群优化(Particle Swarm Optimization,PSO)算法,又称微粒群算法,是1995年Kennedy和Ebehrart通过研究类比鸟类和鱼类的群体行为,而在此基础上提出的一种群智能算法[12]。PSO算法的思想来自人工生命和演化计算理论,利用计算机模仿鸟群飞行觅食中的集体协作行为,通过个体间信息的传递实现集体内部相互影响,引导整个群体向可能解的方向移动,以便快速找出最优解。PSO算法中的粒子是一种在向量空间中以一定速度运行的个体,任何一个粒子n在t时刻所处的位置为xnt=[xn1t, xn2t, …, xndt],具有速度vnt=[vn1t, vn2t, …, vndt],其中d为求解问题的维数。每一个粒子在解空间中所占据的位置代表了当前粒子所对应的一个解,而粒子的速度矢量代表了解变化的方向和程度。算法在迭代的同时会记录每个粒子曾经的最优位置pBest,以及群体曾经的最优位置gBest。
PSO算法中粒子在t+1时刻速度和位置更新公式如下:
$ \begin{array}{l} v_{nd}^{t{\rm{ + }}1}{\rm{ = }}v_{nd}^t + {c_1}rand\left( {} \right) \times \left( {pBeset_{nd}^t - x_{nd}^t} \right) + {c_2}\times rand\left( {} \right) \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {gBest - x_{nd}^t} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l} v_{nd}^t = {v_{\max }}, v_{nd}^t > {v_{\max }}\\ v_{nd}^t = - {v_{\max }}, v_{nd}^t < - {v_{\max }} \end{array} \right. \end{array} $ | (7) |
$ x_{nd}^{t + 1} = x_{nd}^t + v_{nd}^{t + 1} $ | (8) |
在公式(7)中,c1和c2是学习因子,rand()为随机数。在迭代过程中,粒子会根据自身现状、最优位置pBest和群体曾经的最优位置gBest来调整速度,进而更新自身位置,即解空间的解。在适应度函数标准下,粒子通过不断进化,最终得到全局最优解。
PSO算法有易于陷入局部最优解的缺点,如果当前最优解为局部最优,那么一旦所有粒子都收敛于该位置,粒子将很难跳出局部最优。Y.Shi和R.C.Eberhart发现, vndt由于具有随机性,并且本身缺乏记忆能力,导致粒子有扩大搜索空间、探索新的搜索区域的趋势[13]。因此提出,引入惯性权重ω来控制历史速度对当前速度的影响,不再需要vmax,同样达到维护全局和局部搜索平衡的目的。得到粒子速度更新的新公式如下:
$ \begin{array}{l} v_{nd}^{t + 1} = \omega \times v_{nd}^t + {c_1} \times rand\left( {} \right) \times \left( {pBeset_{nd}^t - x_{nd}^t} \right) + {c_2} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;rand\left( {} \right) \times \left( {gBest - x_{nd}^t} \right) \end{array} $ | (9) |
有实验表明,较大的惯性权重有利于展开全局寻优,较小惯性权重有利于局部寻优[14]。
1.3 BPSO搜寻最优子集为了解决工程实际中的组合优化问题,Kennedy和Ebehrart提出了BPSO算法[15]。BPSO算法的粒子位二进制串,位置为离散型的0和1,速度表示一个位串转变成1的概率。首先粒子位置进行二进制编码,每个二进制位利用公式(9)产生速度,再将速度值转换为变换概率,即位变量值取1的概率。这样就不再需要位置更新公式(8)。
速度值表示二进制位取1的概率,一般是采用(10)式Sigmoid函数,将速度值映射到区间[0, 1]上。
$ S\left( {v_{nd}^t} \right) = \frac{1}{{1 + \exp \left( { - v_{nd}^t} \right)}} $ | (10) |
这里S(vndt)表示位置vndt+1取1的概率,粒子通过(11)式来改变位值:
$ x_{nd}^{t + 1} = \left\{ \begin{array}{l} 1, rand\left( {} \right) \le S\left( {v_{nd}^t} \right)\\ 0, rand\left( {} \right) > S\left( {v_{nd}^t} \right) \end{array} \right. $ | (11) |
其中rand()为一个随机数,从区间[0, 1]中随机产生。为了避免S(vndt)过于靠近0或1,定义参数vmax为最大速度值,这里取4。用式(12)来限制vndt,以限制位vndt+1取0或1的概率:
$ x_{nd}^{t + 1} = \left\{ \begin{array}{l} v_{nd}^{t + 1}, \;\; - {v_{\max }} \le v_{nd}^{t + 1} \le v_{\max }^{}\\ {v_{\max }}, \;\;v_{nd}^{t + 1} > v_{\max }^{}\\ - {v_{\max }}, \;\;v_{nd}^{t + 1} < - v_{\max }^{} \end{array} \right. $ | (12) |
考虑到实际情况中雷暴只有两种可能,因此对粒子的位置采用二进制编码,粒子位置向量xnt=[xn1t, xn2t, …, xndt]中的d代表预报因子的个数,二进制位所对应的粒子位置代表着解空间中的每一维,分别与进行建模的雷暴预报因子相对应。若某个二进制位取值为1,则表示对应的预报因子被选为预报因子。这样每一个粒子位置,就对应着一组预报因子,也就是一个普通子集。粒子运动速度用式(13)表示:
$ \begin{array}{l} v_{nd}^{t + 1} = {\omega ^t} \times v_{nd}^{t + 1} + {c_1} \times r \times \left( {pBeset_{nd}^t - x_{nd}^t} \right)/\Delta t + {c_2} \times r \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {gBeset_{nd}^t - x_{nd}^t} \right)/\Delta t \end{array} $ | (13) |
其中,ωt为惯性权重;r为[0, 1]上的随机数;Δt通常取单位时间;c1和c2是学习因子,这里都取2。
为了维护全局和局部搜索平衡,采用自适应惯性权重系数进行优化,惯性权重ωt的迭代公式如下:
$ {\omega ^t} = {\omega _{\max }} - t\left( {{\omega _{\max }} - {\omega _{\min }}} \right)/{t_{\max }} $ | (14) |
其中,t为迭代次数,tmax为最大迭代次数,ωmax和ωmin分别代表最大和最小惯性权重系数,分别取0.9和0.4。
1.4 评分标准及适应度函数为了对模型的预报效果进行预报质量评定,首先定义预报和观测列联表,见下表。
采用TS评分、空报率、漏报率三种评分标准,对漳州站、义乌站、乐东站雷暴预报的结果进行评价,同时也是预报因子的适应度函数,它们的表达式分别为:
$ TS{\rm{评分}}{I_{TS}} = \frac{{{n_{11}}}}{{{n_{11}} + {n_{01}} + {n_{10}}}} $ | (15) |
$ {\rm{空报率}}{I_{FAR}} = \frac{{{n_{10}}}}{{{n_{11}} + {n_{10}}}} $ | (16) |
$ {\rm{漏报率}}{I_{MAR}} = \frac{{{n_{01}}}}{{{n_{11}} + {n_{01}}}} $ | (17) |
其中,TS评分、空报率、漏报率在排除大量预报对象不出现数据的情况下,从三个角度对模型的预报能力进行评价,能够较好地反映出小概率事件的预报水平。因此选择三种评分标准综合考虑,对雷暴预报模型的结果进行评价。
适应度函数是预报因子组合的评价标准,也就是粒子解空间位置所对应解的评价标准,决定了粒子最优位置pBest和群体最优位置gBest的更新。适应度函数选择的好坏将直接影响模型预报效果,为了避免因评价标准不同而影响预报结果,将模型的三种评分标准定为因子组合的适应度函数。
1.5 BDC最优子集的实现BPSO算法筛选BDC最优子集的具体流程如下:
(1) 数据处理。首先,将强对流指数数据进行标准化处理。其次,当某一时次有雷暴数或无雷暴数过少时,如果使用BPSO算法进行因子筛选,不仅计算量大而且效果不理想,这时采用Fisher判别准则进行因子筛选。
(2) 初选因子。与雷暴相关并且物理意义明确的因子有很多,如果从全部因子中筛选出最优子集,明显工作量过大,因此需要从所有相关因子中剔除一部分,主要为对该时次雷暴预报贡献小的因子。根据样本数据,确定初选因子的阈值,以适应度函数为标准,初步筛选出一定数量的因子,这里控制初步筛选的因子数量为12个。
(3) 二进制编码。种群规模越大,则个体越具有多样性,陷入局部收敛的可能性就越小,但是计算量越大。而种群规模太小,则搜索空间受限,难以跳出局部收敛,因此构造BPSO算法的初始种群包含30个粒子。
随机对种群中的每个粒子进行二进制编码,构建粒子位二进制串,其中“1”代表该位置所对应的因子被选入模型,“0”代表该位置所对应的因子没有被选入模型。比如有12个初选预报因子,其中第n个粒子第t步的位置向量xnt=[1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0],表示第n个粒子第t步选取第1、3和9个预报因子,组成预报因子组合,构建BDC雷暴预报模型。考虑到BDC建立的雷暴预报模型,既要包含一定数量物理意义明确的因子,又要避免因选取因子过多而导致误差增加,把模型因子数量控制在3~5个。
(4) 计算初始最优位置。计算初始粒子的后验概率,评价粒子的适应度函数值,记录种群中每个粒子的初始位置及适应度函数值。将粒子的初始位置作为粒子的历史最优位置,比较种群中所有粒子历史最优位置所对应的适应度函数值,选取其中适应度函数值最大的粒子,将该粒子的历史最优位置作为种群的历史最优位置。
(5) 迭代计算。根据粒子的位置向量计算粒子的速度向量,再更新粒子的位置向量。因为要控制模型因子数量,所以就要控制位置向量中“1”的个数。但是位置向量中的一个分量代表着一个预报因子,某一预报因子是否组成预报因子组合,构建BDC雷暴预报模型,不应受到其他预报因子的影响。所以如果直接在位置迭代公式中限制位置向量中分量得“1”的个数,将会直接阻断BPSO算法的择优操作,很有可能无法获得最优解,有悖于BPSO算法的基本原理。
为了在遵守BPSO算法原理的前提下,实现限定预报因子个数的目的,在位置迭代公式运行完之后,如果位置向量中“1”的个数超出了3~5个因子的范围,则逐个改变位置向量中的分量,直到找到满足要求的一系列的新位置向量。通过比较所有新位置向量的适应度函数值大小,找出适应度函数值最大的新位置向量替换旧的位置向量。
比如第n个粒子第t步迭代的位置向量为xnt=[1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0],表示第n个粒子第t步选取第1、3、5、6、9、11个预报因子,一共有6个因子,超过了范围,接下来介绍限定预报因子的具体方法:
步骤1,将位置向量xnt=[1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0]中的“1”依次变为“0”,则得到6个新的位置向量xn1t=[0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0]、xn2t=[1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0]、xn3t=[1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0]、xn4t=[1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0]、xn5t= [1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0]、xn6t=[1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0],这6个新的位置向量中“1”的个数都为5,满足因子的范围要求。
步骤2,计算6个新位置向量的适应度函数值,并且比较6个适应度函数值的大小。
步骤3,假设第3个位置向量的适应度函数值最大,用第3个位置向量xn3t=[1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0]来替代原来的位置向量xnt=[1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0],则最终得到的第n个粒子第t步的位置向量为xnt=[1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0]。
同理,如果第n个粒子第t步迭代的位置向量中“1”的个数小于3,则依次把位置向量中的“0”变成“1”,得到一系列“1”的个数为3的新位置向量,再比较新位置向量的适应度函数值,选取其中适应度函数值最大的新位置向量,代替第n个粒子第t步迭代得到的位置向量。
计算所有粒子的后验概率,评价粒子的适应度函数值,记录种群中每个粒子的位置及适应度函数值。比较粒子现有的位置向量对应的适应度函数值,以及粒子历史最优位置所对应的适应度函数值,选取其中最大的适应度函数值,用其所对应的位置向量作为粒子历史最优位置。比较种群中所有粒子历史最优位置对应的适应度函数值,以及种群的历史最优位置所对应的适应度函数值,选取其中适应度函数值最大的粒子位置向量,将该位置作为种群的历史最优位置。
迭代算法,当种群的历史最优位置满足条件时结束算法,此时种群历史最优位置的粒子位二进制串中“1”所对应的因子,即为筛选出的雷暴预报因子,所有因子的组合即为适应度函数下的最优子集。
2 资料及模型因子的选取 2.1 资料我国雷暴天气多出现于夏季[16],选择雷暴数量较多的8月进行研究。我国的T511数值天气预报系统[17]采用4DVar同化方法,能够提供0~72 h的逐3 h气象要素、物理量场以及时段降水产品,具有较高的准确率,选取T511L60数值预报产品进行雷暴3~48 h的逐3 h预报。使用的资料包括T511L60数值预报产品和单站实况资料。
(1) T511L60数值预报产品。T511L60数值预报产品范围为70°—140°E,10°—55°N;分辨率为1°×1°预报时效为0~72 h,间隔3 h;数值预报产品内容包括单层场和多层场基本要素场;资料长度为2010—2014年逐年8月(北京时间,下同);起报时间20时。
(2) 实况资料。实况资料来源于绘图报。选用的单站包括漳州站、义乌站、乐东站;选用的要素为雷暴;资料长度为同数值预报产品。
2.2 强对流指数考虑到雷暴产生所需要的充沛的水汽、足够的对流冲击力、大量的不稳定能量三个条件[18],选取与雷暴发生密切相关的27个强对流指数场,比如高空温度,高空露点温度,Adedokun1指数(θw-850-θs-850,θw-850是850 hPa湿球位温,θs-850是500 hPa饱和湿球位温),LIsfc指数(T500-T2m→500hPa',其中T500为500 hPa等压面上的环境温度,T2m→500hPa'定义为平均气块从地面2 m高度沿干绝热线上升,到凝结高度后再沿湿绝热线上升至500 hPa时所具有的温度),K指数((T850-T500)+ Td850-(T700-Td700)),Showalter指数(T500-T850hPa→500hPa',T850hPa→500hPa'定义为850 hPa等压面上的湿空气块沿干绝热线抬升,达到抬升凝结高度后再沿湿绝热线上升至500 hPa时所具有的温度),TEI925指数(h500-h925, h=(Cpd-rC1)T+Lv·r+(1+r)gZ, Cpd为干空气的定容比热,C1为湿空气的定容比热,r为比湿,Lv为相变潜热, T为环境温度,Z为位势高度,g=9.8)。
2.3 样本消空对漳州站、义乌站、乐东站的实况资料进行统计,得出在2010—2013年8月逐3 h的资料中发生雷暴的次数分别为85、57、148次,分别占样本的0.086,0.057,0.149。相对于所有的样本数据而言,雷暴发生的次数是比较少的,因此要尽可能多的剔除无雷暴样本同时保留有雷暴样本。利用强对流指数历史最小值,对建模样本进行“消空”。将强对流指数场中,存在两个及两个以上因子等于对应历史最小值的样本剔除,得到消空后的样本中雷暴数的比例分别为0.231,0.192,0.326,实现了尽可能多的剔除无雷暴样本,同时保留有雷暴样本的目的。
2.4 模型因子的选取将消空后样本集中的建模样本根据雷暴有无分成两类,根据适应度函数标准,采用BPSO-BDC方法选取雷暴预报因子的最优子集,可确定单站各预报时次参与建模的因子。表 2给出了漳州站、义乌站、乐东站以TS评分为适应度函数,确立的8月份24 h雷暴预报因子的信息。
雷暴预报模型的步骤可概括为以下四步:
第一步:计算漳州站、义乌站、乐东站2010—2013年8月强对流指数场, 将有两个及以上指数等于历史最小值的样本剔除,对建模样本进行“消空”。
第二步:以适应度函数值为标准,采用BPSO-BDC方法,从强对流指数场中选取雷暴预报的最优子集。
第三步:通过BPSO-BDC建立雷暴预报模型,进行稳定性检验,评估稳定性检验结果。
第四步:利用2010—2014年8月资料进行拟合和预报,使用TS评分、空报率、漏报率进行评价,并与BDC和逐步判别的雷暴预报模型进行对比。
4 稳定性检验从消空样本中依次随机选出5%、10%、15%、20%、25%、30%的数据作为检验样本,其余样本作为建模样本,可得到BPSO-BDC建立雷暴预报方程的稳定性检验结果。图 1所示为以TS评分为适应度函数,对漳州站、义乌站、乐东站雷暴预报方程进行稳定性检验,并以TS评分为标准进行评估的结果。
BDC最优子集雷暴预报方程,其稳定性检验结果的TS评分有高低起伏,但是漳州站、义乌站、乐东站,随机选出5%、10%、15%、20%、25%、30%的样本数据进行检验,TS评分都在0.21~0.35之间变化幅度较小,且处于较高水平,因此可以认为BPSO-BDC方法是稳定的,能够进行雷暴预报。
5 预报结果分析 5.1 最优子集分析使用BPSO算法筛选BDC雷暴预报的最优子集,在迭代过程中,因子组合的适应度函数值是不断变化的,初始最优位置所对应的适应度函数值比较低,随着种群历史最优位置的变化,适应度函数值逐渐增高,到最后最优子集选出后,适应度函数值保持不变。图 2给出了以TS评分作为适应度函数,24 h漳州站、义乌站、乐东站BDC最优子集雷暴预报模型在因子筛选过程中适应度函数值的变化曲线。
从图中可以看出,漳州站和乐东站在前10步和前8步迭代中,TS评分均快速提高,平衡一段后,从第20步和21步又开始了逐渐的提升,最终分别在第39步和第30步进入收敛状态,最优子集所对应的TS评分分别为0.526和0.527。而义乌站则在第8步开始陷入了局部最优位置,因子组合的TS评分达到了0.418,在第37步终于跳出了局部最优位置,后迅速增高,在第42步的时候进入收敛状态,最优子集对应的TS评分为0.477。
5.2 拟合结果分析考虑到实况资料中记录的是观测前3 h到观测时的情况,因此在进行预报结果检验时,如果具体时次预报有雷暴,地面观测资料记录无雷暴,但是3 h前的记录或者3 h后的记录显示有雷暴,那么可以认为该时次实况为有雷暴,雷暴预报成功。
将2010—2013年8月数据带入雷暴预报方程,评价BPSO-BDC建立的BDC最优子集雷暴预报方程的拟合效果,并与BDC和逐步判别模型进行对比。图 3给出了义乌站24 h雷暴预报方程的拟合结果。其中,图中“BDC”对应BDC雷暴预报模型,“逐步判别”对应逐步判别雷暴预报模型,“BDC-MAR、BDC-FAR、BDC-TS”分别对应以漏报率、空报率、TS评分为适应度函数进行因子选取,得到相应三种最优子集后再采用BPSO-BDC方法建立的雷暴预报模型拟合结果。
从义乌站拟合效果可以看出,BPSO-BDC方法建立雷暴预报模型的预报结果较BDC和逐步判别方法有明显的提高。分别使用TS评分、空报率、漏报率作为适应度函数,利用预报方程的TS评分结果都明显好于BDC和逐步判别方法。以TS评分来看,使用TS评分作为适应度函数的模型预报结果最好,然后依次是使用漏报率和空报率作为适应度函数的模型、逐步判别模型、BDC模型。另外,逐步判别方法和使用漏报率作为适应度函数的BPSO-BDC方法,都实现了无漏报,但是逐步判别方法的空报率是所有方法中最高的,因此逐步判别方法的TS评分并不高。使用TS评分作为适应度函数的BPSO-BDC方法不仅实现了空报率为0,而且漏报率也较低,因此其TS评分在所有方法中是最高的。
5.3 预报结果分析通过拟合结果分析,可以得出使用TS评分作为适应度函数,可以综合降低空报率和漏报率,使得TS评分最高。因此在进行预报结果分析时,仅将使用TS评分作为适应度函数的BPSO-BDC方法与BDC和逐步判别方法进行对比。
表 3给出了漳州站、义乌站、乐东站利用2014年8月份T511L60数值预报产品和实况资料,采用BPSO-BDC、BDC和逐步判别方法建立雷暴预报模型进行24 h和48 h雷暴预报的结果。从漳州、义乌、乐东三站TS评分的平均结果来看,24 h的雷暴预报结果中,BPSO-BDC模型的TS评分达到了0.697,其次是BDC模型的0.259,逐步判别模型最低(0.133);48 h的雷暴预报结果中,BPSO-BDC模型的TS评分达0.418,其次为逐步判别模型(0.281),BDC模型仅0.206。无论是24 h还是48 h预报,BPSO-BDC模型的TS评分都明显高于BDC和逐步判别方法建立雷暴预报模型。从空报率的三站平均来看,BPSO-BDC模型的空报率明显低于BDC和逐步判别模型,且48 h空报率较24 h有明显下降。从漏报率的三站平均来看,BPSO-BDC模型24 h预报的漏报率达到0.048,处于极好的水平,48 h预报的漏报率虽然比BDC和逐步判别模型较高,但是相差不大。
对比三种建模方法可以发现,在建模因子的选取时,BDC方法仅仅使用显著性检验标准来筛选因子,因子的选择受显著性检验标准的影响较大,没有考虑到因子间的相互影响,并且不能剔除那些通过检验但是预报效果不佳的因子;逐步判别方法考虑到各因子间的相互影响,综合分析因子组合的判别能力,预报效果有所提高,但是以F检验为标准,虽然筛选出了判别能力显著的因子组合,但是并不能确保预报效果;BPSO-BDC方法则能够立足于试验目的,通过调整适应度函数,达到较好的预报效果。
6 结论与讨论利用T511L60数值预报产品和单站实况资料,采用PSO-BDC方法对漳州站、义乌站、乐东站进行雷暴数值预报研究,并与使用BDC和逐步判别方法进行雷暴预报的结果进行对比,探索雷暴预报最优方法,结论如下:
(1) 采用BPSO算法,通过选取适应度函数,从几十个待选因子中,为雷暴预报模型智能地选取出雷暴预报因子的最优子集。从预报模型因子选用的角度,进一步提升了贝叶斯分类法的预报效果,实现了贝叶斯分类法的有益改进。
(2) 从稳定性检验结果来看,虽然BPSO-BDC建立的BDC最优子集雷暴预报方程的TS评分有高低起伏,但是变化幅度较小且处于较高水平,因此可以认为BPSO-BDC方法是稳定的,能够进行雷暴预报。
(3) 从预报结果来看,BPSO-BDC模型的TS评分和空报率明显优于BDC和逐步判别模型,24 h的平均漏报率也明显低于BDC和逐步判别模型,48 h的平均漏报率与BDC和逐步判别模型大致处于同一水平。可以认为,BPSO-BDC模型的预报效果明显好于BDC和逐步判别模型。
本文从选因子的方法上对BDC雷暴预报模型进行了改进,得到了较好的雷暴预报结果。但是没有考虑不同因子个数对模型预报结果的影响。在后续工作中,将考虑因子个数变化时,模型的预报结果变化情况。另外,考虑预报结果时,没有具体分析预报结果好与坏的原因。在后续工作中,将从建模因子、环流形势等方面考虑预报结果好与坏的原因。
[1] |
卓鸿, 王冀, 霍苗, 等. 不同类型大尺度环流背景下首都国际机场的雷暴特征分析[J]. 暴雨灾害, 2016, 35(04): 371-377. DOI:10.3969/j.issn.1004-9045.2016.04.009 |
[2] |
张文龙, 王迎春, 崔晓鹏, 等. 北京地区干湿雷暴数值试验对比研究[J]. 暴雨灾害, 2011, 30(03): 202-209. DOI:10.3969/j.issn.1004-9045.2011.03.002 |
[3] |
Schmeits M J, Kok K J, Vogelezang D H P, et al. Probabilistic Forecasts of (Severe) Thunderstorms for the Purpose of Issuing a Weather Alarm in the Netherlands[J]. Weather & Forecasting, 2008, 23(6): 1253-1267. |
[4] |
Drira A, Derbel N. Prediction and classification of thunderstorms using artificial neural network[J]. International Journal of Engineering Science & Technology, 2011, 3(5): 1-6. |
[5] |
吴蓁, 俞小鼎, 席世平, 等. 基于配料法的"08.6.3"河南强对流天气分析和短时预报[J]. 气象, 2011, 37(1): 48-58. DOI:10.7519/j.issn.1000-0526.2011.01.006 |
[6] |
张强, 行鸿彦, 徐伟. 基于改进遗传小波神经网络的雷暴预报方法[J]. 南京信息工程大学学报, 2015, 7(3): 221-226. |
[7] |
安洁, 齐琳琳, 王东明. 综合指标叠套方法在海上雷暴预报中的应用[J]. 海洋预报, 2015, 32(1): 58-62. DOI:10.11737/j.issn.1003-0239.2015.01.009 |
[8] |
段旭, 丁圣, 许美玲, 等. 预报因子选取及方程建立人机交互平台[J]. 气象, 2010, 36(11): 120-125. DOI:10.7519/j.issn.1000-0526.2010.11.019 |
[9] |
邓桂骞, 赵跃龙, 刘霖, 等. 一种优化的贝叶斯分类算法[J]. 计算机测量与控制, 2012, 20(1): 199-201. |
[10] |
徐冰霖, 李战怀. 基于贝叶斯准则的航天器故障判决方法[J]. 计算机仿真, 2013, 30(9): 117-120. |
[11] |
袁野, 胡邦辉, 刘丹军, 等. 基于贝叶斯分类判别方法的雷暴预报研究[C]//中国气象学会年会灾害天气事件的预警、预报及防灾减灾论文集. 2009: 1613-1621
|
[12] |
Kennedy J, Eberhart R. Particle swarm optimization[C]//IEEE International Conference on Neural Networks, Proceedings. 1995: 1942-1948
|
[13] |
Eberhart R C, Shi Y. Comparing inertia weights and constriction factors in particle swarm optimization[C]//Proc. of the 2000 Congress on Evolutionary Computation. 2000: 84-88
|
[14] |
Zhang L P, Yu H J, Hu S X. A new approach to improve particle swarm optimization[C]. Lecture Notes in Computer Science. Chicago: Springer2 Verlag, 2006: 134-139
|
[15] |
Kennedy J, Eberhart R C. A discrete binary version of the particle swarm algorithm[C]//IEEE International Conference on Systems, Man, and Cybernetics, 1997. Computational Cybernetics and Simulation. IEEE, 2002: 4104-4108
|
[16] |
王学良, 余田野, 汪姿荷, 等. 1961-2013年中国雷暴气候特征及东亚夏季风影响研究[J]. 暴雨灾害, 2016, 35(05): 471-481. DOI:10.3969/j.issn.1004-9045.2016.05.009 |
[17] |
甘少华, 单军辉, 郭卫东. T511L60全球中期数值预报性能统计评估[J]. 解放军理工大学学报(自然科学版), 2013, 14(1): 107-111. |
[18] |
李新峰, 梅珏, 陈志豪, 等. 2013年9月13日上海地区雷暴天气过程分析[J]. 暴雨灾害, 2015, 34(04): 343-351. DOI:10.3969/j.issn.1004-9045.2015.04.007 |