高阶统计量在形变固体潮汐异常信息提取中的应用①

李希亮，刘希强，张 玲，董晓娜，卢双苓，王 强

（山东省地震局，山东 济南 250014）

0 引言

固体潮汐研究是研究地震成因、地震前兆与地震预测的重要内容之一。近年来通过对全球地震资料的分析，已明确了固体潮汐对地震的发生有着不可忽视的影响。地球固体潮汐是在日、月、地球等天体运动过程中所产生的潮汐形变，通过研究日、月引潮力引起的倾斜、重力值的变化，有可能了解到地壳局部岩石弹性性质的变化。固体潮汐观测到的潮汐畸变、阶变、脉冲等异常变化是震前地壳应力－应变场变化的反映，实际上与地壳介质的应力－应变状态是密不可分的。在“八五”、“九五”、“十五”科技研究期间，利用固体潮观测资料已经对震前地壳形变的中短期异常信息特征有所了解，已取得了一些重要研究成果［1］。

固体潮汐观测数据可以很好地反映地球固体潮汐的变化特征及地震孕育过程的中短期异常特征，而数字化形变资料则提供了该方面研究的可能性。由于采样率的大幅度提高，丰富了潮汐观测的震前变化信息，如地震前短周期的脉动、突跳等短临信息及全球大地震的同震变形等［2］。地球固体潮是一种地球整体受日月引潮力作用所产生的近似粘弹性形变的现象，是一种非平稳过程［3］。对于这种级联影响的非线性非平稳过程，很难从平稳过程的角度直接用固体潮振幅异常（相位异常）进行地震预测。同时大量的固体潮观测实践与研究结果表明，在具有较高灵敏度及精度的潮汐观测中，地震前兆信息对潮汐因子的影响是很微弱的，因此传统的潮汐因子计算方法受到很大的限制。在系统分析和研究了地震预测领域常用的分析方法之后，本文用泰安地震台数字化固体潮汐观测数据，对濮阳ML4．6地震前后异常时间、空间上的特征进行分析，应用高阶统计量方法来获取异常背景信息。有可能对未来发震的时间及地点给出预估计。

1 理论模型与计算方法

1．1 高阶统计量理论［4］

对非高斯信号来说，二阶统计量只是其中一种信息，它不包含相位信息，因此对非最小相位系统的辨识便显得无能为力。随着信号处理应用领域的不断扩大，非高斯、非最小相位、非因果、非平稳信号的处理问题已经成为现代信号处理研究的热点，但基于二阶统计量（如相关函数和功率谱等）的信号分析方法不能从信号中提取非最小相位信息及由于非高斯和非线性带来的信息［5］。

高阶统计量提供了传统的二阶统计量没有的大量丰富信息，能适应于复杂信号或噪声环境。当分析信号为非高斯信号，或者在时间序列分析中对非高斯的AR，MA，ARMA模型进行辨识时，高阶统计量可以有效地抑制高斯或非高斯的有色噪声，并抽取不同于高斯信号的多种信号特征等。理论上，凡是使用功率谱或相关函数等二阶统计量进行过分析和处理，而又未得到满意结果的问题都可以应用高阶统计量方法来进行。高阶统计量方法在地球物理资料数字信号处理和属性参数提取等方面有较强的理论性和实用性。

高阶统计量主要包括高阶矩、高阶累积量、高阶矩谱、高阶累积量谱等四种。

在实际应用中，常用的是随机过程的高阶累积量。设 ｛x（n）｝为平稳随机过程，x1＝x（n），x2＝x（n＋τ1），…，xk＝x（n＋τk－1），则其中k阶累积量和k阶矩为

在k≥3时，k阶累积量和k阶矩统称为高阶统计量。

k阶矩谱定义为k阶矩的k－1维Fourier变换，即

k阶累积谱定义为k阶累积量的k－1维Fourier变换，即

1．2 TOC－AIC模型与计算方法

AIC信息准则由日本学者赤池弘治于1973年提出，其进行震相识别的实质就是求解背景噪声与信号最佳划分点的过程，换句话说就是找出信号时间序列的不稳定点，此点与AIC曲线极小点相对应，即为震相的到时点。AIC的一般计算式为

AIC（n）＝－2log（模型的最大似然函数）＋2n其中n为模型的独立参数的个数。

在震相识别中常见的AIC方法有：自回归AIC方法、基于神经网络的AIC方法、基于小波变换的AIC方法等。AR－AIC 算法（M．LEONARD，1999），即自回归的AIC方法，首先假设可以把地震记录分为两个平稳过程，在每个部分分别建立自回归（AR）模型。这两个AR模型是完全不同的，即一个只包括背景噪声（噪声AR模型）而另一个只包括地震信号（信号AR模型）。当AIC取最小值时，AR模型取最优的阶数，此时AR模型与地震记录这个时间序列最为相符，也就是噪声与地震信号的最优划分点，这一点就是震相的到时点。AIC计算式为

其中为噪声模型的方差为信号模型的方差；mF为噪声AR模型的系数的个数；mB为信号R模型的系数的个数。此做法是在假设地震信号是平稳的，而实际是非平稳的，所以对于低信噪比且初动成隐伏状的地震信号，估计的初至误差较大。

Maeta（1985）提出了不同与AR－AIC的算法，他的做法是直接从地震记录图中计算AIC的值，而没有利用AR系数，其计算式为

在选定的时间窗内对k逐点搜索，AIC的值反映x这个时间序列以k为分界的前后两段方差的最小值，AIC极小值的点即为P波到时。

AIC方法在震相识别中的应用已经很成熟，AR－AIC算法和VAR－AIC算法是比较常用的两种方法，在快速、有效的处理地震记录中起到了很大的作用，推进了地震预警的发展。

震相到时检测技术就是假设震相到时前后的地震记录是两个不同的稳态过程，开始估计的到时把时间序列分成两个局部稳态的时间序列，以检测震相的到时。但是信噪比高低和震级大小对AIC值的影响起决定性作用。同样的在地震前兆数字化观测记录中的时间序列中，如果存在异常信号的介入，将出现两个不同的稳态过程。两个局部稳态的时间序列的划分也很关键，用以检测地震前兆异常信号。

考虑到地震前兆信号的复杂性，不平稳性，以及低信噪比等特点，本文提出了一种直接从前兆数字化观测资料中计算AIC函数的方法。具体思路是用两个时间段数据的三阶累积量代替VAR－AIC方法中的方差，称为TOC－AIC方法（third－order cumulant）。AIC函数表示为

模拟实验证明，当随机高斯噪声中存在异常信号时，此方法可以抑制噪声，精确地提出异常信号的初至信息。在不影响观测分析结果精度的前提下，要对观测数据进行预处理，尽量消除固体潮汐的周期性变化对计算的影响，排除干扰，突出异常信息。如果有新的地震前兆信号到达时，地震前兆数字化记录中就会出现突变和渐变。采用TOC－AIC方法，从背景噪声和地震前兆异常信号中找出信号时间序列中的不稳定点，也就是AIC函数的最小值，来识别地震前兆异常的初至。

2 震例研究

2．1 台址及资料情况

泰安台位于泰山南麓，地处鲁中隆起，泰山山前第四纪活动断裂北侧。台基为太古代花岗片麻岩，岩体完整致密均匀，测量信噪比高，可靠性良好，为全国一类形变基准台。仪器洞室处坐标是东经117°07′22″，北纬36°12″35′。主洞进深约76m，最大覆盖厚度约29m，室温年变幅约0．06℃。布设有水平摆倾斜仪、垂直摆倾斜仪，还平行布设数字化SS－Y伸缩仪、DSQ水管倾斜仪各一套，共计11个测项，均为“九五”数字化前兆台网改造之后新上的数字化仪器，已经积累了近10年的资料，资料可靠性高（图1）。数字化倾斜仪的精度可以达到高于模拟记录的水平，水平摆倾斜仪精度、稳定性都高于垂直摆倾斜仪［6］。这就是在同一洞室里面水平摆的潮汐异常最明显的原因所在。由于台站及其周边地区中强地震较少，相关的形变潮汐短临异常出现的很少，因此濮阳ML4．6地震潮汐短临异常的识别对今后的地震研究工作有一定的借鉴意义。

图1 泰安台形变仪器布设示意图Fig．1 Schematic diagram of deformation instruments at Tai＇an seismic station．

2．2 正常背景场固体潮汐TOC－AIC方法研究

在正常背景场下，固体潮汐观测数据的TOCAIC方法研究结果如图2所示。正常背景场下的处理结果比较平滑，反映不出异常信息的介入。

图2 2006年1月泰安台垂直摆NS向的原始资料和TOC－AIC值Fig．2 The recorded data and TOC － AIC values of NS vertical pendulum at Tai＇an station in January 2006．

2．3 异常背景场固体潮TOC－AIC方法研究

2006年4月9日在与山东莘县相邻的河南濮阳境内发生了ML4．6地震，震源深度12km，震中距泰安地震台165km。在地震前泰安台的潮汐形变观测不同程度的出现了短临异常。将TOC－AIC方法分析引入潮汐形变资料处理中，明显地提取到了我们想要得到的地震前兆信息。图3、4显示结果对异常信息的反映程度虽有不同，但都能明显的反映出异常信息的介入。其中垂直摆的资料分析结果较好的反映了临震信号的介入（图4（a）、（b）），水平摆的分析结果则较好的反映出了异常信号的介入和临震信号的介入（图4（c）、（d））。

图3 濮阳地震前后倾斜仪原始记录Fig．3 The recorded data of pendulum before and after the Puyang earthquake．

3 异常机理分析及其干扰排除

3．1 区域应力场背景特征

山东地区整体平均主压应力方向为近EW向。图5为1970年以来初动符号得到的山东及邻区中小地震的震源机制解得到的主压应力水平投影，主压应力轴优势方位为80°左右，主张应力轴优势方位为350°左右。殷海涛等在山东GPS观测网的基础上，对山东地区地壳运动特征进行了分析，认为山东地区的主应变场呈NEE－SWW向挤压，这是由于鲁西块体主要受来自青藏块体NEE向的挤压［7］。

濮阳ML4．6地震发生在聊考断裂带的西侧。聊考断裂带为一被第四系掩盖的隐伏断裂带，南起河南兰考，北至聊城以北与齐广断裂交会。全长达270多公里。根据物探和钻探资料，断裂走向NNE20°～40°，倾向 NW，倾角35°～60°，为一东升西落的正断层。新第三纪以来，仍表现出强烈活动性，也是山东地震源地之一。濮阳ML4．6地震的震源机制解显示主破裂面走向为312．00°，倾角82°，滑动角13°，主压应力方向为85°，主张应力方向为356°，与该区正常情况下主压应力场方向85°±15°、主张应力场方向350°±15°比较一致，同时符合山东地区地壳运动特征。初步推测此次地震的发震构造为聊考断裂带内规模较小的NW或NNW之次级隐伏断裂，NNE向聊考断裂为控震构造。

3．2 异常机理分析

泰安台震前出现短、临前兆异常，其中水平摆两项短期异常最为显著可靠。结合濮阳ML4．6地震震源机制解和泰安台与震中的位置可见，地震发生前泰安台附近主要受来自NNW方向的压应力和NEE方向的张应力作用，与水平摆出现的NE向倾斜异常相一致。这说明泰安台水平摆出现的异常很可能是震中区应力增强导致的近源区应力变化的结果。随着断层上应力的增强—释放—减弱，异常的表现也随之出现相应的变化，直到最后恢复到原来正常水平［8］。泰安台与濮阳震中之间有多条断层纵横交错，虽然许多断层之问没有直接相连，但它们之间互相交错，使得彼此又具有内在的联系。尤其是NNE向的聊考断裂带与EW、NNW向的多条断裂相交，其活动对与之相关的次级断裂将起到控制作用。因而，泰安台出现的倾斜异常应归为“场兆”异常，是大区域应力场及伴生的多条断层综合效应的结果。

图4 泰安台不同摆不同方向的TOC－AIC值Fig．4 TOC － AIC values of vertical and horizontal pendulum in different directions at Tai＇an station．

3．3 异常干扰排除

形变观测受多种因素影响，其中包含了气象变化、人为作用、仪器故障及未知的诸多情况，使得干扰排除成为识别异常的一项必不可少的工作。但在目前的认识水平下，干扰形态与地震前兆形态极为相似，给干扰与前兆的判定带来很大困难［9－10］。经过与台站工作人员的共同落实，排除了气压、雷电、洞室进人、仪器故障及附加载荷等干扰因素对数字化形变高频记录的影响，提高了濮阳ML4．6地震泰安台形变潮汐短临异常的可靠性。

4 结论与讨论

地震孕育过程中固体潮潮汐伴随着地震孕育产生一定的变化，这种变化在具有较高灵敏度及稳定性的潮汐观测过程中，采用较高分辨率的高阶统计量分析方法能较好的识别。处理结果中的低点就是异常信号的介入点，这对地震异常信息的提取有指导意义，此方法能较好的提取潮汐形变短临异常。

图5 山东及邻区单震震源机制解主压应力轴水平投影Fig．5 The horizontal projection of principal compression stress axis from single earthquake in Shandong province and its adjacent areas．

泰安台水平摆潮汐短临异常形态明显，震后完全恢复到正常年变形态。异常初期表现为地面NEE向倾斜，正好背向震中，与多数倾斜异常统计结果一致。研究表明，在同一洞室水平摆和垂直摆对地震异常的反应程度有差异，虽然它们同是记录倾斜潮汐变化的仪器，但由于原理与结构的不一致，对地震异常的反映程度不一样，高阶统计量方法分析结果中的异常初至也不一致。高灵敏度的水平摆倾斜仪器扑捉地震前兆的能力要好一些，对地壳表面运动反应更敏感，更有利于获取短临异常。泰安台距震中约165km，进一步验证了6．5级以下的中强地震，突变点集中分布在震源附近200km范围内的结论［11］。

本文只是应用某些前兆资料来说明高阶统计量分析方法在提取固体潮潮汐异常的可行性和有效性，对于前兆机理和前兆异常的分析也是比较肤浅的。本文采用高阶统计量方法对地震前兆信号进行分析与特征提取，研究了信号偏离高斯分布的程度，为地震前兆的识别提供了一条新的思路。期望给地震预测提供一些有用的参考。

［1］ 刘序俨，洪星．论潮汐波与地震波［J］．大地测量与地球动力学，2003，23（2）：73－76．

［2］ 张燕，吴云，刘永启，等．小波分析在地壳形变资料处理中的应用［J］．地震学报，2004，26（增）：103－109．

［3］ 李瑞浩．重力学引论［M］．北京：科学出版社，1988．

［4］ 王书明，朱培民，李宏伟，等．地球物理学中的高阶统计量方法［M］．北京：科学出版社，2006．

［5］ 李希亮，刘希强，董晓娜，等．高阶统计量方法在地球物理学中的应用与展望［J］．西北地震学报，2010，32（2）：201－205．

［6］ 王梅，宋治平，李峰，等．由泰安地震台形变资料看固体潮数字化观测运行［J］．地震研究，2004，27：49－56．

［7］ 殷海涛，李 杰，张 玲，等．基于GPS观测网的山东地区地壳运动特征分析［J］．西北地震学报，2008，20（3）：276－281．

［8］ 李杰，李希亮，卢双苓，等．濮阳ML4．6地震前泰安台形变异常特征分析［J］．大地测量与地球动力学，2007，27（4）：100－104．

［9］ 王梅，李峰，孔向阳，等．数字化形变观测干扰识别［J］．大地测量与地球动力学，2004，24（1）：94－98．

［10］ 王梅，宋治平，李峰，等．形变数字化资料综合分析［J］．大地测量与地球动力学，2003，23（4）：60－64．

［11］ 牛安福，张晶，高福旺，等．地壳持续加速变形与地震关系研究［J］．大地测量与地球动力学，2002，22（1）：29－33．