基金项目:地震科学联合基金(A07058)和2009年震情跟踪定向工作任务(2009029901)基金联合资助.
(Earthquake Administration of Beijing Municipality,Beijing 100080,China)
HPPPR model,academic numerical simulation,seismic synthesis forecast
备注
基金项目:地震科学联合基金(A07058)和2009年震情跟踪定向工作任务(2009029901)基金联合资助.
在综合考虑经典投影寻踪算法特点的基础上,针对投影寻踪计算中存在的一些不利因素,给出相应的解决思路。利用数值仿真技术进行基于粒子群优化算法与厄密特多项式构建的投影寻踪回归模型建模能力与计算精度的检验,再将其应用于多维地震时间序列综合建模预测中。计算结果和进一步分析表明,基于粒子群优化算法与厄密特多项式构建的投影寻踪回归投影寻踪模型具有简单、快速、有效的特点,在实际地震综合预测建模中取得了满意的效果,可作为地震综合预测的一种回归分析方法。
Considering its characteristics,we give the solution to the flaws of the conventional Project Pursuit Regression(PPR)algorithm.We use the numerical emulation technique to test the modeling ability and calculating accuracy of PPR that based on particle swarm optimization and Hermite multinomial,and apply HPPPR to the synthetic,seismic forecasting modeling of multi-dimensional seismic time series.The result shows that PPR model that based on particle swarm optimization and Hermite multinomial is simple,fast and effective,and it is satisfying in the real seismic synthesis forecasting modeling and it can serve as a regressing method.
引言
地震综合预测方法在地震预测研究中占据着十分重要的位置,是地震预测的主要手段之一(梅世蓉等,1993; 时振梁等,1997; 张晓东等,2003; 武安绪等,2005; 武安绪等,2008a; 武安绪等,2008b),其中投影寻踪(Projection Pursuit,简称PP)是最重要的综合预测方法之一(朱令人等,1994; 周仕勇等,1995; 赵翠萍,周仕勇,1999; 武安绪等,2003)。投影寻踪的原理是将高维数据点阵投影到低维(K=1~3)子空间上,寻找能反映高维数据的结构特征的投影方向,然后在低维空间上对数据结构或特征进行分析,以达到研究和分析高维数据的目的。投影寻踪回归方法的思想是将输入变量向若干个一维方向进行投影,对每一个方向的投影值,分别用一维函数拟合,然后用此若干个一维拟合函数的和逼近回归函数,其表达式为(朱令人等,1994; 武安绪等,2003; Friedman,Stuetzle,1981; 成平,李国英,1986; 李祚泳,1998)
Y=Y^-+∑Mm=1βmGm(∑Pj=1αTmjX).(1)
式中,(X,Y)为一对随机变量; Y^-表示Y的均值; P为输入空间维数; M为岭函数的个数; amj是第m个投影方向的第j个分量; βm为权值,表示第m个岭函数对模型输出贡献的大小; Gm为第m个光滑岭函数,满足∑Pj=1α2mj=1,E(gm)=0,E(g2m)=1。
在地震综合预测研究与分析中,虽然投影寻踪回归(Projett Pursuit Regression,简称PPR)方法已获得广泛应用,取得了显著效果(朱令人等,1994; 周仕勇等,1995; 赵翠萍,周任勇,1999; 武安绪等,2005),但在使用式(1)时不难发现该方法存在一些不利因素:岭函数G无具体的函数表达式,必须采用一个较为庞大的函数集合给出; 计算导数时需用差分代替微分; 进行预测时,必须采用函数集合进行内插,计算过程较为复杂。α一般采用高斯—牛顿法优化,易于获得局部解,获得最优解则有一定困难; 分组优化参数α、β与岭函数时,需要先固定其中两个,然后优化另一个,完成后再固定两个优化另一个,逐个进行,使得已经优化好的参数会在下一个参数优化时失去其最优性,也有陷于局部解的危险性,不易获得方法的全局最优解。以上问题都会使PPR的学习过程变得复杂,外推技术繁琐,计算效率低下,效果难以确定。针对这些应用中存在的不利因素和问题,可基于粒子群优化算法(Particle Swarm Optimization,简称PSO)(Kennedy,Eberhart,1995; Eberhart,Kennedy,1995)和厄密特多项式(冯康等,1978; Weaver,1991)来重新构建投影寻踪回归方法(PPR that based on Hermite multinomial and particle swarm optimization,简称HPPPR),简化建模的技术与过程,改变建模结构和学习策略,提高建模的效率与能力。为了检验HPPPR方法的效果,我们将其进行多维样本数值模拟,然后用于地震综合预测试验中,探讨在地震综合预测研究中的有效性和实用性。
1 HPPPR地震综合预测回归模型
1.1 HPPPR投影寻踪回归模型中的岭函数拟合方式经典投影寻踪回归模型中岭函数G的具体表达式没有明确给出(朱令人等,1994; Friedman,Stuetzle,1981),而是采用庞大的简单函数集合去光滑逼近,这样会让使用者在实际投影寻踪回归建模与预测中选择的岭函数时造成一定困难。为了避免使用庞大的简单函数集合,且能保证岭函数G拟合的光滑效果和拟合能力,可以选择一调频调幅与正弦波混和信号
x(t)=[1+0.2sin(2π7.5t)]cos[2π30t+
0.5sin(2π15t)]+sin(2π120t)(2)
作为理论模型实验信号,该模型信号由一基频30 Hz、调制频率15 Hz的调频调幅成分和120 Hz的正弦波叠加而成(图1),以此,通过对各类多项式的分析和调频调幅实际检验,认为只有厄密特多项式具有较强的光滑、细节、趋势、动态、非线性多尺度拟合描述能力,其描述问题复杂程度的能力随阶次的增加会不断增强。因此可采用可变阶的正交厄密特多项式h拟合一维岭函数,其表达式(冯康等,1978; Weaver,1991)为
hn=(n!)-1/2π1/42-(n-1)/2Hn(x)φ(x), -∞<x<∞.(3)
式中,x=αTX; n!代表多项式阶数n的阶乘; φ(x)=1/((2π)1/2)e-(x2)/2; Hn(x)为厄密特多项式, H0(x)=1, H1(x)=2x, …, Hn(x)=2[xHn-1(x)-(n-1)Hn-2(x)]。
这样HPPPR中的岭函数可采用的具体形式为拟合多项式
G(x)=∑Nn=1cn.hn(x).(4)
式中,c为多项式系数,N为多项式阶数。
1.2 粒子群优化算法(1)粒子群优化算法代替高斯—牛顿算法的意义
在经典投影寻踪回归学习过程中,一般采用高斯—牛顿算法进行优化,需要利用差分代替微分算法运算,不易获得全局最优解(朱令人等,1994; 武安绪等,2003; Friedman,Stuetzle,1981),而粒子群优化算法是新发展的优化算法,不再需要进行微分运算,具有高效的搜索思路,设置参量也较为简单,因此可采用粒子群优化算法代替高斯—牛顿算法。
(2)粒子群优化算法的基本原理
粒子群优化算法是Kennedy和Eberhart(1995)受鸟群觅食行为的启发于1995年提出来的。该算法是模拟鸟群飞行觅食的行为,通过鸟之间的集体协作使群体达到最优的过程,所以粒子群优化算法是基于群体智能理论的优化算法。在粒子群优化算法中,所有粒子都有一个由适应度函数决定的适应度和一个由速度Vi=(vi1,vi2,…,vid)决定的在搜索空间单位迭代次数的位移。其中,适应度函数由优化目标定义。粒子群优化算法首先初始化为一群随机粒子(随机解),其中第i个粒子在d维解空间的位置表示为Xi=(xi1,xi2,…,xid),然后算法通过迭代搜索产生最优解。每一次迭代,粒子通过两个极值更新自己的速度和位置。一个是粒子从算法迭代初始到当前迭代次数搜索生成的最优解及个体极值pbesti=(pbesti1,pbesti2,…,pbestid),另一个则是整个粒子种群当前的最优解——全局极值。粒子根据以下公式进行速度和位置的更新:
vid=w×vid+c1×rand()×(pid-xid)+c2×
rand()×(g-xid),(5)
xid=xid+vid.(6)
式中,vid为粒子i在第d维的速度,且vd∈[-vdmax,+vdmax],xid是粒子i在第d维的当前位置; w为惯性权重,它使粒子保持运动惯性,起着调整算法全局和局部收缩能力的作用。w取大值可使算法具有较强的全局搜索能力,取小值则算法倾向于局部搜索,一般的做法是将w初始取0.9并使其随迭代次数的增加递减至0.4,这样可先侧重于全局搜索,使搜索空间快速收敛于某一区域,然后采用局部精细搜索以获得高精度解; rand()是均匀分布在(-1,+1)之间的随机数; c1,c2为学习因子,分别调节向个体最好粒子和全局最好粒子方向飞行的最大步长,若太小,则粒子可能远离目标区域,若太大则会突然向目标区域飞去,或飞过目标区域。通常令c1=c2=2。粒子飞行速度由最大值vdmax限制,若vdmax太大,粒子将飞离最好解,太小则会陷入局部最优。假设将搜索空间的第d维定义为区间[-vdmax,+vdmax],则通常vdmax=kxdmax,0.1<k<1.0,每一维都采用相同设置方法。
式(5)中第1部分可以理解为粒子先前的速度或惯性; 第2部分可以理解为粒子的“认知”行为,表示粒子本身的“思考”能力; 第3部分理解为粒子的“社会”行为,表示粒子之间的信息共享与相互合作。粒子在解空间内不断跟踪个体极值与全局极值进行搜索,直到达到规定的最大迭代次数或小于规定的误差标准为止。
1.3 利用粒子群优化算法和厄密特多项式构建投影寻踪回归学习方法基于粒子群优化算法和厄密特多项式构建的投影寻踪回归学习方法是:逐个增加岭函数,直到满足收敛条件。学习过程可分为两个步骤:① 根据公式(1),随机选择m个初始投影方向,对任意一个投影方向,计算X的一维投影值Z,找到与数据对(Z,Y)拟合最好的岭函数,确定多项式的阶数N和系数c。根据拟合结果计算m个方向上的误差平方和; ② 根据粒子群优化算法的优化规则进行参数优化,直到选出误差最小的那个投影方向作为第一个投影方向。确定了第一个投影方向后,再优化参数,计算第一个方向上的拟合残差,如果满足精度要求则输出结果,否则就增加一个岭函数。按照同样的方法进行参数优化,直到增加的岭函数个数满足拟合精度的要求。
2 HPPPR模型应用分析
2.1 理论数据模拟分析假定一数学函数(王士同,1998)
y=(1.0+(x1)1/2+1/(x2)+x3-1.5)2.(7)
根据此公式产生离散的前17个训练样本和后17个测试样本(图2)建立投影寻踪回归模型,检验多维样本建模与内符预测能力。经过粒子群算法优化HPPPR获得的最佳投影寻踪回归模型由3个岭函数组成,其投影方向分别为A1=[-0.853 79,0.358 52,0.246 57],表明参量x1分量在该投影方向上起最主要作用; A2=[-0.077 92,0.145 80,0.986 24],表明参量x3分量在该投影方向上起最主要作用; A3=[-0.385 06,-0.894 09,0.228 76],表明参量x2分量在该投影方向上起最主要作用。从3个岭函数的投影方向可看出,这3个分量都起作用,说明产生理论样本的分布均匀性与实际情况相符。其学习与测试结果为:学习样本的误差为0.730 264,实测值与预测值的相关系数为0.991 541; 测试样本的误差为0.945 434,实测值与预测值的相关系数为0.975 485; 全部样本的误差为0.814 813,实测值与预测值的相关系数为0.986 548(图2)。如果按±1.0误差门限计算,内检结果的正确率为100%,而按相对误差10%计算,正确率为90%。
通过投影寻踪回归模型对理论样本的学习与内符预测检验结果可证明:基于粒子群优化算法
图2 理论数据通过HPPPR学习与检验获得的实测值与预测值序列(实线表示实测值; 虚线表示预测值)
Fig.2 Real value and prediction value series obtained through HPPPR study and test with theoretical data(solid line:real value,dotted line:prediction value)和厄密特多项式的PPR模型拟合能力较强,训练速度较高,误差收敛较快; 系统受训练样本空间分布影响较小,可得到全局最优解,获得的回归模型具有较强的表达能力与外推泛化能力。因此,笔者利用该方法对实际的多维地震样本资料进行综合建模预测,研究地震资料的综合处理能力。
2.2 基于中国大陆典型震例的HPPPR地震综合预测分析例子(1)影响震级水平的前兆因子复杂性与地震综合预测研究的意义
地震综合预测是地震预测工作中的重要方法之一,也是地震工作者不断追求的目标,但地震前兆十分复杂,成百上千条前兆异常各具特色,其与未来震级水平很难找到具有明确物理意义的对应关系特征量。如何从复杂的前兆现象中提取带有共性的特征参量,从而建立地震综合预测的判据,是地震综合预测的难点,也是关键。地震学家利用经典PPR方法,在这方面已进行了大量的有益的探索,取得了明显的效果(朱令人等,1994; 周仕勇等,1995; 赵翠萍,周仕勇,1999; 武安绪等,2003)。本文中,笔者利用由粒子群优化算法和厄密特多项式构建的PPR模型建立多维地震综合预测模型,探索多种前兆异常、地震活动参量多维时间序列和中强地震震级之间的非线性映射关系,这样可以利用某个地区出现异常后或某一个地区的历史地震演化特征,为预测某地区某时间段将要发生的地震震级水平提出综合预测判据。这在一个地区出现多个异常后要判断可能发生的震级水平,或在年度会商中预测某一个地区未来一年的震级水平,都将是十分重要的研究工作,这些工作在地震发生机制还没有解决之前,具有积极的探索意义。下面笔者利用中国大陆典型实际地震震例数据开展地震前兆综合预测模型分析,研究HPPPR地震综合预测模型的有效性。
(2)中国震例中地震前兆异常基础资料的筛选
从20世纪70年代以来中国大陆发生了大量的中强地震,积累了丰富的震例资料。充分利用这些资料对未来的地震预测工作具有重大意义。笔者将利用这些实际震例资料,建立HPPPR地震综合预测模型。地震研究中,震级是最重要的预测参数之一。影响震级水平的因素很多,笔者采用王炜等(1999)整理出的中国震例中的因素量,即地震条带、地震空区、应变释放、地震频次、b值、地震窗、波速比、短水准、地倾斜、地电、水氡、水位、应力、宏观和异常数量共15项前兆指标作为HPPPR模型的输入,实发震级作为系统的输出,建立多维地震综合预测HPPPR模型,根据前兆异常量预测某地区可能发生的地震的震级水平。
(3)建立地震前兆异常持续时间与实发地震震级之间的HPPPR综合预测模型
由于样本的无序性,笔者按照王炜等(1999)整理出的中国震例样本顺序,从中选择前30个震例作为学习样本,最后15个震例作为内检样本,通过HPPPR模型的学习与预测,获得实测震级的预测值(图3)。经过粒子群算法优化HPPPR获得的最佳投影寻踪回归模型由4个岭函数组成,其投影方向分别为A1=[0.094 225,0.221 052,0.008 896,0.282 276,0.190 953,0.103 522,0.029 945,0.386 918,0.474 192,-0.238 101,0.183 297,0.474 192,0.006 500,-0.205 153,-0.287 369],表明地倾斜、水位分量在该投影方向上起最主要作用; A2=[0.155 244,0.091 555,-0.283 544,-0.050 775,-0.630 332,-0.053 172,0.084 225,0.359 576,-0.140 549,-0.443 728,0.138 252,0.235 230,-0.200 028,0.112 516,0.065 121],表明b值分量在该投影方向上起最主要作用; A3=[0.229 026,-0.500 773,0.144 175,-0.353 133,-0.059 426,0.175 419,-0.027 787,-0.236 302,-0.173 341,-0.302 473,0.441 094,0.015 269,-0.377 031,-0.026 328,-0.034 521],表明地震空区分量在该投影方向上起最主要作用; A4=[-0.358 981,-0.430 502,0.364 036,0.439 754,-0.020 925,0.001 574,0.058 240,0.170 358,0.206 046,-0.237 612,0.041 558,0.013 788,-0.146 901,0.440 076,-0.104 703],表明宏观分量在该投影方向上起最主要作用。从4个岭函数的投影方向可看出起作用的4个不同前兆因子分量,说明前兆因子与所预测的地震的对应情况是复杂的。由图3统计可得:学习样本的误差为0.030 114,实测值与预测值的相关系数为0.938 604; 测试样本的误差为0.051 217,实测值与预测值的相关系数为0.963 980; 全部样本的误差为0.034 216,实测值与预测值的相关系数为0.946 914。如果按±0.5级震级误差门限计算,内检结果的正确率为100%,而按相对误差5%计算,正确率为93%。因此,该模型可作为地震前兆综合预测模型。从HPPPR方法对震例样本的学习和内检结果看,模型学习能力较强,具有较高的外推泛化能力。
3 结语
(1)HPPPR模型用厄密特多项式作为岭函数的具体拟合形成,具有光滑内插与外延优点,与经典的投影寻踪回归方法相比,无需记录岭函数的参数表,计算过程较为简单,且可以揭示高维地震数据中更多的非线性信息(武安绪等,2003); 用简便的粒子群算法优化投影方向,不需要求导运算,不进行复杂的数学运算,增强了投影寻踪回归方法的实用性。应用表明,粒子群参量优化能力很强,速度快,精度高,在参数设置得当的情况下,可很快获得全局解。因此,粒子群优化算法可以作为一些地震模型参量反演的主要优化工具。
(2)HPPPR模型中的参数优化不再需要分组,而是同时完成,避免了循环内部拟合抵消的耦合现象,加快了计算速度,获得的解是全局最优解。多维数值仿真与地震综合预测分析例子证明,基于粒子群优化算法和厄密特多项式构建的HPPPR地震综合预测模型是合理的,建模能力和预测精度都较高,说明该方法是有效的,可用于地震综合预测研究。
- 成平,李国英.1986.投影寻踪——一种新兴的统计方法[J].应用概率统计,(3):138-143.
- 冯康,张建中,张绮霞,等.1978.数值计算方法[M].北京:国防工业出版社.
- 李祚泳.1998.投影寻踪的理论及应用进展[J].大自然探索,(1):47-50.
- 梅世蓉,冯德益,张国民,等.1993.中国地震预报概论[M].北京:地震出版社.
- 时振梁,汪良某,傅征祥,等.1997.中国大陆中长期强震危险性预测方法研究[M].北京:海洋出版社.
- 王士同.1998.神经模糊系统及其应用[M].北京:北京航空航天大学出版社.
- 王炜,蒋春曦,张军,等.1999.BP神经网络在地震综合预报中的应用[J].地震,15(2):118-126.
- 武安绪,吴培稚,谷庭瑶,等.2005.Modular模糊神经网络及其在官厅水库地震危险性估计中的应用[J].大地测量与地球动力学,25(增刊):32-36.
- 武安绪,张晓东,张永仙,等.2008a.查表法设计地震活动模糊预测预警系统[J].地震,28(2):101-107.
- 武安绪,张永仙,张晓东,等.2008b.地震前兆综合预测支持向量机模型研究[J].地震,28(3):55-60.
- 武安绪,朱红彬,邢成起,等.2003.应力释放模型在华北地区四个主要地震活动构造带中长期地震活动危险性趋势分析中的初步应用[M]//段尔焕,杜现泽,徐人平.地震研究与工程抗震.北京:原子能出版社.
- 张晓东,傅征祥,张永仙,等.2003.1999~2002年地震预报研究进展[J].地震学报,25(5):479-491.
- 赵翠萍,周仕勇.1999.新疆主要地震区PP回归综合预报模型研究[J].西北地震学报,21(1):37-43.
- 周仕勇,朱令人,邓传玲.1995.PP聚类在震群分析中的应用研究[J].地震学报,17(3):312-321.
- 朱令人,周仕勇,邓传玲.1994.地震综合预报的新方法——投影寻踪回归[J].地震学报,16(增刊):1-9.
- Weaver H J.1991.离散和连续傅里叶分析理论[M].王中德,张辉译.北京:北京邮电学院出版社.
- Eberhart R,Kennedy J.1995.A new optimizer using particle swarm theory[C]//Proceedings of the Sixth International Symposium on Micro Machine and Human Science,IEEE service center,Piscataway,NJ,Nagoya,Japan,39-43.
- Friedman J H,Stuetzle W.1981.Projection Pursuit Regressio[J].Amer Statist Assor,76:817-823.
- Kennedy J,Eberhart R.1995.Particle swarm optimization[J].IEEE International Conference on Neural Networks,4(27):l942-1948.