2020, Vol. 41
2. 中国科学院 测量与地球物理研究所 大地测量与地球动力学;
3. 中国地质大学(武汉)复杂系统先进控制与智能自动化湖北省重点实验室，湖北 武汉 430074;
4. 中国地质大学(武汉) 智能地学信息处理湖北省重点实验室，湖北 武汉 430074;
5. 中国地质大学(武汉) 地球科学学院，湖北 武汉 430074
2. State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics (CAS), Wuhan 430077, China;
3. Hubei Key Laboratory of Advanced Control and Intelligent Automation for Complex Systems, China University of Geosciences (Wuhan), Wuhan 430074, China;
4. Hubei Key Laboratory of Intelligent Geo-Information Processing, China University of Geosciences (Wuhan), Wuhan 430074, China;
5. School of Earth Sciences, China University of Geosciences (Wuhan), Wuhan 430074, China
地球天然脉冲电磁场(the Earth’s natural pulse electromagnetic field, ENPEMF)是指可在地球表面接收的由天然场源产生的综合电磁总场[1-2].地震、滑坡等地质灾害现象可在地表产生甚低频(very low frequency, VLF)信号脉冲波动，“微破裂机-电转换”机制和“地壳波导”是上述电磁现象的机理之一[3].ENPEMF信号中的电磁异常信息具有潜在变化趋势及典型的非平稳特征，可反映地质活动的孕育发展趋势，可用于震前电磁异常监测分析[4].
随着信息技术、人工智能及机器学习理论的不断发展，电磁信息预测模型逐渐成为前沿热点.由于影响ENPEMF信号的场源较多，具有非平稳信号的特点，很难根据脉冲波形归纳其规律.非线性模型中的神经网络方法是基于经验风险最大化原则的机器学习算法，非线性拟合能力较强，可以对采集到的震前ENPEMF信号进行建模，并拟合其强度趋势的变化特点.
径向基函数(radial basis function, RBF)神经网络可逼近任意的非线性函数，其优点是学习速度快、非线性逼近能力强，具有良好泛化能力，可为ENPEMF信号构建非线性预测模型.本文提出基于混沌参数优化RBF算法的预测模型对其强度趋势进行预测，为数据分析和灾害监测提供支持.首先采用混沌理论对实测ENPEMF数据进行分析，其中假邻近法(false nearest neighbor, FNN)及自相关函数法分别求得嵌入维数和延迟时间等混沌特征参数，并优化RBF神经网络；采用训练完成的混沌参数优化RBF神经网络模型对ENPEMF数据进行预测，并与传统RBF神经网络进行比较.结果表明，混沌参数RBF算法可基于混沌理论确定动态系统的混沌特性，从采集到的震前ENPEMF信号强度数据中找到其变化态势，并预测14 d(4月7~20日)的ENPEMF数据强度趋势，且预测效果及精度均优于传统RBF预测模型.
1 混沌理论由于震前ENPEMF信号的产生机理复杂、孕育过程非线性，其信号强度数据具有非平稳特点和混沌特性.因此，本文引入混沌理论对其数据内部特征进行挖掘，找到隐藏的混沌特点.假设一段时间内采集的震前ENPEMF信号数据为{x(tj), j=1, 2, …, n}，其中n表示采集的数据点数，通过混沌理论中的假邻近法及自相关函数法对数据进行处理，得到数据变化形式：
|
(1) |
式中，τ和m分别表示震前ENPEMF信号数据的延迟时间和嵌入维数，用于描述该信号隐藏的混沌特征并为RBF神经网络输入节点提供判断依据.
1.1 假邻近法时间序列的本质是将系统高维空间坐标的运动轨迹投影到低维空间.当嵌入维数较小时，系统空间轨道中本来相距很远的相点相互挤压折叠，未能充分展开，这些点为假邻近点[5].随着嵌入空间维数的增加，轨道逐渐展开，投影到低维空间的假邻近点随之分离.当所有的假邻近点消失时所对应的最小嵌入空间维数即为最佳嵌入维数.给定正整数m，可构造m维重构向量：
|
(2) |
在m维重构空间中，采用欧氏度量来计算ym(n)的紧邻点y′m(n)：
|
(3) |
将维数从1维增加到m+1维.m+1维空间重构向量如式(4)所示：
|
(4) |
紧邻点y′m+1(n)如式(5)所示：
|
(5) |
在m+1维空间中，考察紧邻点y′m+1(n)是否与m维空间的紧邻点y′m(n)一致.如果相同，则m为嵌入维数.在实际计算中，欧氏距离为
|
(6) |
给定参数Rτ，如果满足以下条件：
|
(7) |
在m+1维空间中，计算重构向量ym+1(n)与其紧邻点y′m+1(n)的距离，若满足：Rm+1(n)/Rm(n)≥Rτ，其中10≤Rτ≤50时结果稳定，则y′m+1(n)为重构向量y(m+1)(n)的假邻近点，即判据2为
|
(8) |
对所有的重构向量，利用判据找出邻近点中的假邻近点，并记录下所有假邻近点的数目FN(m).继续增加维数，当找到一个整数mε使得FN(mε)=0时，mε即为所求嵌入维数.当假邻近点所占比率即假邻近率随着嵌入维数的增加趋于平稳不再降低时，所对应的嵌入维数m为最佳嵌入维数.本文通过统计假邻近点数的比率随嵌入维数升高逐渐减小，最后维持不变的情况，确定最优嵌入维数.
1.2 自相关函数法由于实际时间序列长度有限且存在噪声，选取合适的延迟时间至关重要.延迟时间τ过小，将使重构的系统由于相关性较强造成相空间的挤压，不能充分展示系统的动力特征；延迟时间τ太大，会造成相邻两时刻的动力学形态剧烈变化，使构造的相空间比实际空间复杂.
自相关函数法可在降低相关性的同时保证原动力学的系统信息不丢失，使重构相空间能充分展现系统拓扑性质和几何性质.首先写出时间序列的自相关函数，然后做出自相关函数随时间变化的函数图，找到自相关函数首次达到零点时对应的时间即为延迟时间τ.自相关函数定义为
|
(9) |
自相关函数值随时间变化逐步下降，当其下降到初始值的(1-1/e)时对应的时间为所求延迟时间.
2 RBF神经网络算法径向基函数(RBF)神经网络是一种包含输入层、隐含层和输出层的三层神经网络[6].其中输入层与隐含层之间为非线性变换，从隐含层到输出层为线性变换.
在RBF神经网络中，输入层仅作为通道传输信号.隐含层中神经元的变换函数为径向基函数，通过非线性变换可将信号从输入层传递到隐含层.输出层是对输入信号的响应.RBF神经网络结构可根据具体问题在训练阶段进行自适应调整.
在网络中，X=(x1, x2, …, xn)T为输入样本，Y=(y1, y2, …, yn)T为输出响应.RBF神经网络算法需要求解基函数的中心Ci，基函数的宽度Di以及隐含层到输出层的权值Wi三个参数.
RBF神经网络的训练过程分为两步：首先进行无监督学习，计算输入层与隐含层之间的和，得到隐含层输出：
|
(10) |
其中：i=1, 2, …, N; ‖X-Ci‖为欧氏范数.
然后采用最小二乘法求隐含层与输出层之间的权值ωi.最终得到RBF神经网络的输出y：
|
(11) |
将经过参数训练的RBF神经网络用于预测混沌时间序列.可利用混沌时间序列数据校正上述权值参数，提高神经网络的非线性泛化能力.其中，当输入层节点个数为混沌时间序列的嵌入维数m时，预测结果较好；隐含层节点数目根据实验实时调整确定：预先设定RBF神经网络的精度值，隐含层节点个数递增，当神经网络达到预设精度时，该节点个数即为神经网络的隐含层节点数.
本文提出的基于混沌参数优化RBF神经网络算法的预测模型工作流程如图 1所示.
|
图 1 本文算法预测模型工作流程 Fig.1 Workflow of proposed forecasting model |
首先对实测ENPEMF数据进行混沌分析，然后分别采用假邻近法和自相关函数法求得最优嵌入维数和延迟时间等混沌特征参数，将得到的参数作为确定输入层节点个数的依据并优化RBF神经网络；对参数优化的RBF神经网络进行训练，学习其内部混沌特征，最后用训练完成的混沌参数优化RBF算法预测ENPEMF信号的强度趋势.
3 实验仿真隐藏在混沌时间序列内部的某些特征信息可以通过延迟时间及嵌入维数表现出来.本文用经典的Rossler混沌时间序列验证假邻近法及自相关函数法获得嵌入维数和延迟时间等参数的可行性.
Rossler系统可用式(12)微分方程组进行描述：
|
(12) |
选取参数a=b=0.2，c=5，初值x(0)=-1，y(0)=0，z(0)=1，积分时间步长h=0.05，生成长度为3 000的连续混沌时间序列如图 2所示.
|
图 2 Rossler混沌时间序列 Fig.2 Rossler chaotic time series |
Rossler混沌时间序列在有限区域内运动时趋向于一个稳定的点，完全展开系统内部的混沌特性.嵌入维数的取值范围为[1, 8]，阈值的判别门限范围为[2, 15]，Rossler时间序列长度为3 000，采用假邻近法计算时间序列x分量的嵌入维数m，结果如图 3所示.
|
图 3 假邻近法求嵌入维数 Fig.3 Using FNN method to find the embedded dimension |
当嵌入维数从1增加到4时假邻近率急速下降；当嵌入维数达到5时，假邻近率趋于平缓，此时的嵌入维数达到理想值，为所求最佳嵌入维数，即Rossler混沌时间序列的嵌入维数为5.
采用自相关函数法求延迟时间τ，选择x分量序列计算延迟时间τ，仿真结果如图 4所示.
|
图 4 自相关函数法求延迟时间 Fig.4 Using autocorrelation function method to calculate the delay time |
图 4中，直线为初始值的(1-1/e)，曲线为自相关函数曲线，选取值为自相关函数下降到初始值的(1-1/e)时所对应的时间，即Rossler混沌时间序列的延迟时间τ为16 s.
取迭代后长度为2 000的时间序列，其中前500个数据进行训练，后1 500个数据进行预测.以嵌入维数m=5作为RBF神经网络的输入节点判断依据，采用训练完成的混沌参数优化RBF算法对数据进行预测.图 5为预测模型的仿真结果.“×”表示真实值，“-”表示预测值.图 6为预测结果的绝对误差.
|
图 5 Rossler混沌时间序列的预测结果 Fig.5 Forecasting result of Rossler chaotic time series |
|
图 6 预测结果的绝对误差 Fig.6 The absolute error of forecasting result |
图 5和图 6中，混沌参数优化RBF算法可以较好拟合经典的Rossler混沌时间序列，且预测误差较小.
4 混沌参数优化RBF算法的预测研究本文使用俄罗斯科学院托木斯克分院GR-01型设备接收ENPEMF信号.在武汉九峰地震台放置了3台设备，接收方向为W-E和N-S，3个通道为CN1，CN2，CN3.设备记录了ENPEMF信号的AH数据(超过设定阈值的脉冲幅度)及NH数据(超过设定阈值的脉冲个数)，可表征地表天然磁场的强弱[7].设备的工作频率为甚低频段：5~25 kHz，在武汉设置的接收频率为14.5 kHz.
ENPEMF信号为非周期、非平稳信号，具有明显混沌特性，设备输出为数字量化后的信号，数据存储格式为：时间-幅度-脉冲数(t-AH-NH)[8].幅度单位仅为信号包络大小变化的参照量，已不具有原来量纲的直接意义.
对采集的ENPEMF数据进行混沌特性分析，找到其内部隐藏的混沌特征及趋势变化特点，结合RBF神经网络算法对其信号强度进行预测，识别孕震信息[9].本文采用Grassberger和Procaccia提出的关联维数算法(G-P算法)求解ENPEMF信号关联维数，可判断信号是否具有混沌特性[10].
对于时间序列x(1), x(2), …, x(t)，其长度为M，对其进行相空间重构，得到向量：X(t)=[x(t), x(t+τ), …, x(t+(m-1)τ)], 其中t=1, 2, …, N, N=M-(m-1)τ.给定正数ε足够小，当空间向量间的距离小于ε时，向量关联.关联向量的关联积分表达式为
|
(13) |
其中θ(·)为Heaviside阶跃函数，满足式(14)：
|
(14) |
当时间序列的长度N→∞，半径ε→0时，关联积分与半径的关系为
|
(15) |
其中，D为所求关联维数，变形后得到式(16)：
|
(16) |
给定一系列半径ε和嵌入维数m，作半径随嵌入维数变化的关联积分图组，用最小二乘法对图中lbC(N, ε)~lbε最接近直线的一段拟合最佳直线，该直线斜率即所求关联维数D.
本文数据于2013年4月20日在中国四川省芦山地震期间收集，地震的位置和地震台的位置如图 7所示[11].
|
图 7 芦山地震的位置和地震台的位置 Fig.7 Location of the Lushan earthquake and seismic station |
图 8为4月10~20日通道2的AH数据，箭头指向中国四川省7.0级芦山地震发生时间.
|
图 8 4月10~20日通道2的AH数据 Fig.8 The AH data of the channel 2 from April 10th to 20th |
地震发生前的11 d内，ENPEMF数据从14日到15日有较大的峰值变化，在16日回落至正常水平.在17日和18日观测到显著的峰值变化，在19日信号又跌落至正常.因此，在地震发生前ENPEMF信号脉冲强度会发生剧烈变化.
本文选择经过平滑及归一化的20 d ENPEMF数据作为实验数据(4月1~20日).其中前6 d数据(4月1~6日)为模型的训练样本，后14 d数据(4月7~20日)作为模型的预测样本.利用G-P算法计算4月1~6日ENPEMF信号的关联维数.信号的关联积分组在一定范围内呈近似直线分布；随着嵌入维数的增加，直线斜率增大，且最后关联维数趋于稳定，说明ENPEMF信号具有混沌特性.
对选取的ENPEMF信号中前6 d数据(4月1~6日)进行预处理，采用假邻近法计算嵌入维数m，自相关函数法计算延迟时间τ，得到图 9和10，从中了解震前ENPEMF信号数据的混沌特性.
|
图 9 假邻近法求嵌入维数 Fig.9 Using FNN method to find the embedded dimension |
|
图 10 自相关函数法计算延迟时间 Fig.10 Using the autocorrelation function method to calculate the delay time |
图 9中，随着嵌入维数的增加，当嵌入维数为4时，假邻近率趋于平稳，此时的嵌入维数即ENPEMF信号数据的最佳嵌入维数.图 10中，信号自相关函数达到初始值的(1-1/e)时τ=5 s，即ENPEMF信号数据的延迟时间为5 s.用得到的参数确定RBF神经网络的输入节点个数为4，进而将训练完成的混沌参数优化RBF算法用于数据预测.
选择4月1~6日数据作为训练样本，训练混沌参数优化RBF神经网络，其中输入层有4个节点，输出层1个节点，隐含层有6个节点.选择径向基高斯函数作为隐含层神经元传递函数，输出为线性函数.最后，利用训练完成的混沌参数优化RBF神经网络预测模型和传统的RBF神经网络预测模型分别实现对4月7~20日ENPEMF数据的单步预测.图 11为所提混沌参数优化RBF神经网络预测模型的结果，图 12为传统RBF神经网络预测模型结果.
|
图 11 混沌参数优化RBF预测结果 Fig.11 The results of the chaotic parameters optimized RBF forecasting model |
|
图 12 传统RBF预测结果 Fig.12 The results of the traditional RBF forecasting model |
图 11和12中，分别采用混沌参数优化RBF预测模型与传统RBF预测模型对4月7~20日ENPEMF信号强度数据进行预测，两种模型均可以模拟出采集到的地震前14 d(7~20日)实际ENPEMF信号强度的波动.对于整体数据范围，传统RBF神经网络模型不能较好地跟踪实际值的变化，而本文所提预测模型对ENPEMF信号具有较好的跟踪拟合性能；对于17日的数据剧烈波动时刻，混沌参数优化RBF预测模型相较于传统的RBF预测模型拟合效果更好，具有较好的预测结果，误差较小，预测优势明显.为更精确评估所提预测模型的预测效果，选取绝对误差作为ENPEMF数据预测精度评价指标，结果如图 13所示.
|
图 13 两种算法的绝对误差值对比 Fig.13 Comparison of absolute error between two algorithms |
混沌参数RBF神经网络预测模型仅在17日剧烈波动时段存在预测误差，整体上的预测误差均小于传统RBF神经网络预测模型.为验证混沌参数优化RBF算法预测结果的稳健性和可靠性，本文采用互相关系数对两种算法预测值与实际值之间进行量化测量，如式(17)所示：
|
(17) |
得到混沌参数优化RBF神经网络模型的互相关系数结果为r1=0.800 4，略大于传统RBF神经网络模型的互相关系数r2=0.792 6.因此，本文所提优化算法的预测效果优于传统RBF神经网络算法.
综上，混沌参数优化RBF预测模型能够较好地反映采集到的强震前14 d(4月7~20日)ENPEMF信号强度变化的趋势和规律，可以满足对强震前ENPEMF信号强度趋势的预测需要，期望为地震和地质灾害前的电磁预测发挥积极的作用.
5 结论1) 本文提出了一种混沌参数优化RBF神经网络预测模型.通过混沌理论确定的特征参数优化RBF神经网络，进而对震前ENPEMF强度数据进行预测，并与传统RBF神经网络算法进行比较.
2) 本文混沌参数优化RBF预测模型可有效地预测强震前14 d(4月7~20日)ENPEMF信号的强度趋势，且预测效果优于传统RBF神经网络算法，期望为地震及地质灾害前的电磁监测分析提供依据.
| [1] |
Vorobyov A A. On probability of electrical discharges in the Earth's interior[J]. Geologiya i Geophysika, 1970, 12: 3-13. |
| [2] |
Sibgatulin V G, Peretokin S A, Khlebopros R G. Fundamental peculiarities of the entropy model of energy processes in seismic areas[J]. Earth Science Frontiers, 2007, 14(6): 222-225. DOI:10.1016/S1872-5791(08)60013-5 |
| [3] |
Malyshkov Y P, Malyshkov S Y. Periodicity of geophysical fields and seismicity:possible links with core motion[J]. Russian Geology and Geophysics, 2009, 50: 115-130. DOI:10.1016/j.rgg.2008.06.019 |
| [4] |
郝国成, 白雨晓, 吴敏, 等. 基于BSWT-DDTFA方法的地球天然脉冲电磁场震前信号时频分析研究[J]. 地球物理学报, 2018, 61(10): 4063-4074. (Hao Guo-cheng, Bai Yu-xiao, Wu Min, et al. Time-frequency analysis of the Earth's natural pulse electromagnetic field before earthquake based on BSWT-DDTFA method[J]. Chinese Journal of Geophysics, 2018, 61(10): 4063-4074. DOI:10.6038/cjg2018L0431) |
| [5] |
Kugiumtzis D. State space reconstruction parameters in the analysis of chaotic time series:the role of the time window length[J]. Physica D, 1996, 95(1): 13-28. DOI:10.1016/0167-2789(96)00054-1 |
| [6] |
Gencay R, Liu T. Nonlinear modeling and prediction with feed forward and recurrent networks[J]. Physica D, 1997, 108(1/2): 119-134. |
| [7] |
Hao G C, Bai Y X, Liu H, et al. The Earth's natural pulse electromagnetic fields for earthquake time-frequency characteristics:insights from the EEMD-WVD method[J]. Island Arc, 2017, 27: e12256. |
| [8] |
郝国成, 龚婷, 董浩斌, 等. 基于聚类经验模态分解的地球天然脉冲电磁场时频与能量谱分析:以芦山Ms7.0地震为例[J]. 地学前缘, 2015, 22(5): 231-238. (Hao Guo-cheng, Gong Ting, Dong Hao-bin, et al. Time-frequency characteristics and energy analysis of the Earth's natural pulse electromagnetic fields based on ensemble empirical mode decomposition:The Lushan Ms7.0 earthquake as an example[J]. Earth Science Frontiers, 2015, 22(5): 231-238.) |
| [9] |
郝国成, 陈忠昌, 赵娟, 等. 基于NSTFT-WVD变换的芦山Ms7.0级地震前后地球天然脉冲电磁场信号时频分析[J]. 地学前缘, 2016, 23(1): 276-286. (Hao Guo-cheng, Chen Zhong-chang, Zhao Juan, et al. Time-frequency analysis of the Earth's natural pulse electromagnetic field signal before and after the Lushan Ms7.0 earthquake based on NSTFT-WVD transform[J]. Earth Science Frontiers, 2016, 23(1): 276-286.) |
| [10] |
Lorenz E N. Atmospheric predictability as revealed by naturally occurring analogues[J]. Journal of the Atmospheric Sciences, 1969, 26: 636-643. DOI:10.1175/1520-0469(1969)26<636:APARBN>2.0.CO;2 |
| [11] |
Hao G C, Tan F, Hu X Y, et al. A matching pursuit-based method for cross-term suppression in WVD and its application to the ENPEMF[J]. IEEE Geoscience and Remote Sensing Letters, 2019, 16(8): 1304-1308. DOI:10.1109/LGRS.2019.2894223 |