摘要 ：重标极差分析法( R/S ， rescaled range analysis method)， 是研究分数布朗运动及自然现象的自仿射分形的有力工具之一， 对于所有时间序列分析都有着广泛的用途. 通过对天水台绝对地磁观测 Z 分量数据序列 R/S 分析， 计算分段数据序列的Hurst指数 H 发现， 昆仑山口西 M _S 8.1地震和汶川 M _S 8.0地震震前2-3个月， H值减小， 分维数D增大， 相关系数C 为负值. 这一动态变化表明， 地震前数据序列由持久性序列发展为反持久性序列， 观测数据的成分发生了改变. 这与地壳内部应力变化和岩石破裂实验结果相符. 所以， 这一动态变化过程可以看作是强震前的一种前兆信号. 关键词 ： R/S 分析 地磁 Z 分量 Hurst指数 分形维数 前兆信号 AA study on fractal Brownian motion of geomagnetic observations before large earthquakes

1. China Earthquake Networks Center, Beijing 100045, China;
2. Tianshui Observation Station, Tianshui， Gansu Province 741020, China;
3. Earthquake Administration of Jiangsu Province, Nanjing 210014, China

Abstract : Rescaled range analysis ( R/S ) method is one of the good tools for studying fractal Brownian motion of time-series data and natural phenomena self-affine fractal, being applied widely for all time-series. In this paper, this method is used to analyze absolute geomagnetic Z component observations at Tianshui Observatory. The Hurst exponent H was calculated for each sub-sequence of data in specific time intervals. Our analysis result indicates that, 2-3 months before Kunlun M _S 8.1 earthquake and Wenchuan M _S 8.0 earthquake, H value decreased and fractal dimension value D increased. Such dynamic variations manifest that the observations varied from persistent behaviour to antipersistent time-series before the earthquakes and the composition of the observed data changed. This variation is consistent with internal crustal stress change and rock failure experiment results. So we may view this dynamic variation as a premonitory symptom of the strong earthquakes. Key words : R/S analysis geomagnetic Z component Hurst exponent fractal dimension premonitory symptom

地震的孕育、 发生是一个极其复杂的过程. 很早以前，人们就注意到地震前后会伴随着地磁异常的变化. 但真正开展震磁效应实验和理论研究是在20世纪50年代以后，而有目的地进行震磁前兆观测与地震预报实践则开始于20世纪60年代( 丁鉴海等，1994 ). 已有震例与野外试验的观测结果显示，在地震孕育、 发展、 发生的各个阶段都不同程度地伴随有地磁异常现象. 目前，震磁效应的研究已取得了一定进展.

长期地震地磁观测表明，观测分量的值是不断变化的. 变化的原因是多方面的，如与太阳、 地球等天体运行有关的昼夜变化、 季节变化、 年变化等周期性变化； 还有在地震的孕育发展过程中，由于地下应力作用，地下岩石的物理、 化学性质就要发生变化，从而导致地下岩石磁性的改变，于是在地面上观测到的地磁数据就会发生局部的微弱变化. 由此可见，地磁场的变化是一个非线性动力学过程，也是一个由规律成分和随机成分叠加而成的时间序列过程. 地震成因的电磁信号淹没在各种干扰信号中，处理这些电磁信号以发现微弱的地震成因的电磁信号是非常重要的，特别是在观测点离震中相对较远的情况下.

20世纪初，英国水文学家Hurst 在对尼罗河潮汐数据进行研究时提出重标极差法( R / S ，rescaled range analysis method). 该方法主要考虑统计意义下的事物发展变化的前后“记忆”特征，用赫斯特指 数H 定量分析时间序列，是研究分数布朗运动及自然现象的自仿射(self-affine)分形的有力工具之一，对于所有时间序列分析都有着广泛的用途.

目前，多维分形分析方法在地震资料处理和异常提取中得到了初步应用. 用来提取地震前兆信号的多维分形方法是基于单维分形分析，而近期破坏过程(裂隙的传播过程)和多维分形结构间定量模型的建立( Kiyashchenko et al ，2003 )对此给以了有力支持.

张继红(1997) 研究表明，连续变化的地磁场具有分维结构. 日本的 Gotoh 等(2003) 用分形方法分析了关岛(Guam)和伊豆(Izu)地震期间超低频(ULF)地磁观测数据，发现两次强震前，ULF地磁观测数据分形维数增加. 希腊 Vallianatos等(2004) 将 R / S 分析方法运用到时间序列的超低频地磁三分量(水平分量X和Y，垂直分量Z)的数据分析，进行单维和多维分形研究，得出了地磁垂直分量比水平分量具有更强的多维分形性. 因此，在后面的应用分析中，我们选取地磁观测垂直分量Z作为研究资料. 1 分数布朗运动与Hurst 指数

Hurst是一个水文工作者，他长期研究尼罗河流量的变化，以决定水库的排水量. 在多年的水文数据中，他发现数据不服从布朗运动及正态分布的特性，而且如果有一年水量较大，次年的水量也往往较大. 在布朗运动中，H=0.5，但一般情形H并不等于0.5. 后来 H 被称为Hurst指数( 彼得斯，1999 ). 此后，Hurst发现一大类自然现象随时间演化的行为都不能用布朗运动刻画，即2 H 都不是整数，因而称为分数布朗运动.

布朗运动是一个特殊的随机过程. 考虑在 x 轴上随机走动的微粒，每隔τ秒跳跃的步长为+ξ或-ξ，设步长服从正态分布

1)当H=0.5时，C=0，即过去和未来增量间的相关系数为0. 表明现在不影响未来. 这说明增量过程是一个独立的随机过程，布朗运动是其特殊情况.

2)当 H ≠0.5时，为分数布朗运动. 这又包括两种情况：

a. 0.5< H <1，有持久性效应. 表明过去增长意味着未来这种趋势将继续下去. 反之，过去的减少趋势就平均而言，意味着未来的连续减少. H 越接近1，趋势越明显； H 越接近0.5，逐渐趋于随机性. 这种长期记忆作用使得随机过程呈现一定的趋势，0< C <1，增量间的相关系数为正，且 C 越接近1，其相关性越强.

b. 0< H <0.5，-0.5< C <0，增量间的相关系数为负，称为反持久性效应(antipersistent). 如果过去是增长的，则下一时刻下降的可能性更大； 反之，过去是下降的，则下一时刻上升的可能性更大. 反持久性效应的强度取决于 H 接近0的程度. H 越接近0，则 C 越接近-0.5，增量间趋于负相关. 因为这种过程有更多的频繁反转，所以比布朗运动波动更加剧烈，这意味着其分形维数大于布朗运动的维数. 2 R / S 分析过程

Hurst 指数可由 R / S 统计法确定( 高红兵等，2001 ). 对一个已知的时间序列

为了比较不同类型的时间序列，赫斯特用原来观测值的标准差去除极差. 极差与标准差之比为R(n)/S(n)，记为R/S，称为重标极差(rescaled range)，这个“重标极差”应该随时间增加.

赫斯特建立以下关系：

且分型维数D=2-H( Smirnova1 et al ，2001 ).

将上述过程用VC++程序实现，可以求出任一时间序列的 H 值. 并将上述分析方法应用到天水台绝对地磁观测 Z 分量数据的分析及地震异常信息的提取. 3 两次强震前天水台绝对地磁观测 Z 分量 H 指数变化分析 3.1 天水台绝对地磁观测情况

天水地震台位于甘肃省，南—北地震构造带中北段( 图 1 ). 地磁台基主要为黄土层，覆盖层厚度约20 m，基岩为泥质砂岩. 台址地磁场磁性分布均匀，梯度为0.2—1.7 nT/m.

绝对地磁 Z 分量观测是“九五”期间配备的高精度G-856质子旋进磁力仪，分辨率为0.1 nT，以及集光电为一体的无磁性地磁绝对观测仪CTM-DI. 该仪器为模拟观测，没有明确的观测频段，输出日均值，数据开始于1994年初，1998—2004年数据毛刺较多，其余时间观测数据平稳. 整个观测过程数据连续性很好. 3.2 昆仑山口西 M _S 8.1地震和汶川 M _S 8.0地震

2001年11月14日17时26分青海省昆仑山口西发生了 M _S 8.1地震，其震中为36.20°N、 90.90°E，震中距天水台为1 350 km. 由于震中位于山中无人区，没有人员死亡( 陈棋福等，2008 ).

2008年5月12日14时28分四川省汶川发生了 M _S 8.0强烈地震，震中位于阿坝州汶川县(31.0°N，103.4°E)，距天水台约450 km. 全国大部分地区均有震感，地震造成几十万人的人员伤亡，严重破坏和巨大损失举世震惊.

图 1 为两次地震震中与天水台相对位置图. 3.3.1 数据的选取和预处理

将天水台绝对地磁观测 Z 分量数据从2000年1月1日—2008年6月16日的3 090个数据作为一个时间序列，作 R / S 分析，得 H=0.896，D=2-H=1.1， 由式 C =2 ^{2 H -1} -1得出 C =0.73. 同时分析发现，观测数据存在不规则的弱的日变化、 季节变化和年变化成分. 这表明天水绝对地磁观测 Z 分量数据序列变化是一个非线性动力学过程，也是一个规律成分和随机成分叠加而成的时间序列过程； 有分维特征，遵循分数布朗运动的性质； 时间序列具有持久性效应，现在的变化趋势是过去变化趋势的延续，且一直持续下去，这种长期记忆作用使得随机过程呈现强的上升趋势，增量间有较强的正相关性.

经分析发现，数据上升趋势对后面研究分析结果有较大的影响. 因此，对数据进行进一步分析前，需要对分析数据按年分段作去趋势预处理. 由于数据中存在弱的不规则周期成分，这些成分为非准周期性的，若作准周期处理，势必会带来不利影响. 另外，后面的 R / S 分析中选取时间段为一个月左右，与上述数据周期存在显著时间差，故不予考虑.

为分析较长时间数据变化情况，选取的两个震例分析数据均在两年以上. 昆仑山口西 M _S 8.1地震，选取数据从2000年1月1日—2001年12月31日，共731个数据点； 汶川 M _S 8.0地震，从2006年1月1日—2008年6月16日，共898个数据点.

图 2 中a ₁ 和b ₁ 分别是选取的两次强震前后天水台绝对地磁 Z 分量原始观测数据变化曲线. 整体上，数据曲线呈上升趋势. 地震前，没有发现明显的异常变化.

天水台绝对地磁观测 Z 分量数据曲线图
(a) 昆仑山口西 M _S 8.1地震前后 Z 分量数据变化；(b) 汶川 M _S 8.0地震前后 Z 分量数据变化(a ₁ )、 (b ₁ )表示天水台绝对地磁观测 Z 分量原始曲线；(a ₂ )、(b ₂ )表示一般多项式分段拟合曲线；(a ₃ )、(b ₃ )表示原始数据与拟合值之差 Fig.2 Variation of absolute geomagnetic Z component at Tianshui observatory
(a) Z component variation before Kunlun 8.1 earthquake; (b) Z component variation before Wenchuan 8.0 earthquake(a ₁ ) and (b ₁ ) represent absolute geomagnetic Z component observation curve; (a ₂ ) and (b ₂ ) represent the general polynomial fitted curve with intervals; (a ₃ ) and (b ₃ ) represent difference between original data and fitted curve

若将选取数据作为两个独立的时间序列，作 R / S 分析，得两序列的 H 值分别为0.73和0.88，表明两个序列都是持久性序列，与原始数据持久上升趋势相符. 原始时间序列按年去掉趋势，作分段曲线拟合，并求与原始数据的差值，得到 图 2 中a ₃ 和b ₃ 曲线. 两个数据序列去趋势处理前后作相关分析，得相关系数分别为0.33和0.08. 这表明，处理后原观测数据的强上升趋势基本消除. 下面的分段 R / S 分析数据是基于预处理后的数据. 3.3.2 地磁数据 R / S 处理与分析

为了进一步分析两次强震前地磁 Z 分量观测数据的变化，对去掉趋势后的时间序列作分段 R / S 分析.

对昆仑山口西8.1级地震去掉趋势的时间序列分段，取30个数据点为一个独立的子时间序列，共得到24个完整的子序列，时间段为2000年1月1日—2001年12月22日. 作分段 R / S 分析，分别计算每个子序列的 H 值，并计算相应的相关系数 C 和分维数 D ，得到 图 3 .

昆仑山口西 M _S 8.1地震天水台绝对地磁观测 Z 分量 R/S 分析图
(a) Hurst指数 H 变化曲线； (b) 相关系数 C 变化曲线； (c) 分维数 D 变化曲线 Fig.3 R/S analysis of Z component data at Tianshui observatory before Kunlun M _S 8.1 earthquake
(a) Hurst exponent H ; (b) correlation coefficient C ; (c) fractal dimension D

由 图 3 可看出，各参数变化大体分3个阶段，出现两个低值点. 第一个低值点出现在2000年4月，其 H =0.56，相应的 C =0.08， D =1.4. 表明虽然 H 值减小，分维数有所增加，数据的成分有改变，但数据序列仍具有持久性，发展变化趋势没有改变，过去的增长意味着这种趋势在未来将继续下去； 第二个低值点出现在2001年6月25日—7月24日，此时， H =0.44＜0.5，表明时间序列由持久性序列变为反持久性序列，相关系数 C =-0.08，过去增量和未来增量不相关. 这一过程有更加频繁的反转，所以比布朗运动的波动更加剧烈. 分形维数 D =1.56，大于布朗运动的分形维数，达到最大值，表明 Z 分量观测数据的成分改变. 随后，数据以持久性序列变化，到2001年11月14日昆仑山口西 M _S 8.1地震发生，有近3个月的时间.

对汶川8.0级地震前一年趋势数据，取40个数据点为一个独立的子时间序列，共得到22个完整的子序列，时间段为2006年1月1日—2008年5月29日. 分别计算每个子序列的 H 值，计算相应的相关系数 C 和分维数 D ，得到 图 4 .

汶川 M _S 8.0地震天水台绝对地磁观测 Z 分量 R/S 分析图
(a) Hurst指数 H 变化曲线； (b) 相关系数 C 变化曲线； (c) 分维数 D 变化曲线 Fig.4 R/S analysis of Z component data at Tianshui observatory before Wenchuan M _S 8.0 earthquake
(a) Hurst exponent H ; (b) correlation coefficient C ; (c) fractal dimension D

由 图 4 可看出，从2006年1月1日—2008年5月29日，天水台绝对地磁观测垂直分量 Z 数据变化大体分为3个阶段.

第一阶段，从2006年1月1日—2006年11月16日. 在2006年1月1日—2006年10月7日这一时间段，Hurst 指数0.5＜ H ＜1， C ＞0， D ＞1. 表明观测数据序列呈非随机的分数布朗运动，增量间的变化不是独立的； 对于每个子时间序列，具有持久性，数据序列的变化发展趋势不会改变，过去的增长意味着这种趋势在未来将继续下去.

2006年10月8日—11月16日， Z 分量数据序列的 H =0.5，数据序列的发展变化出现短暂的随机运动，即布朗运动. 现在与未来增量间的相关系数 C =0，表明现在不影响未来，增量过程是一个独立的随机过程.

第二阶段，从2007年11月17日—2008年4月19日. 在2007年11月17日—2008年3月10日这一时间段，Hurst 指数0.5＜ H ＜1， C ＞0， D ＞1，数据序列呈非随机的分数布朗运动，增量间的变化不是独立的，呈弱正相关性.

2008年3月11日—4月19日，这段时间序列的 H =0.44＜0.5，时间序列呈现反持久性序列. 现在与未来的相关系数 C =-0.08＜0，增量间不相关. 数据变化趋势发生反转，比布朗运动的波动更加剧烈. 分形维数 D =1.56，大于布朗运动的分形维数1.5，达到最大值，这表明 Z 分量观测数据的成分发生改变.

第三阶段，汶川地震发生后，绝对地磁 Z 时间序列的 H 值增加且大于0.5，相应的相关系数 C ＞0，分形维数 D 减小，由反持久性序列发展为持久性序列，恢复序列的增强趋势.

通过以上分析，昆仑山口西地震前3个月左右和汶川地震前2个月左右，天水台绝对地磁观测 Z 分量数据时间序列呈现 H <0.5，分形维数增加. 呈临界态的系统常在空间上呈现分形结构(自相似性)，而在时间上则呈现闪烁噪声(“flicker noise”或“1 / f 噪声”，时间上的标度不变性和长程时间关联). 对地面观测，这一动态过程表现为震前相应分维数增加，和时间序列由随机的、 持久性序列向反持久性序列发展( Gotoh et al ，2003 ). 在自组织临界状态，震源区有发生大量地震的危险，因为任何小尺度的改变(例如，流体填充小的空隙空间)都会牵连到整个系统的改变. Gotoh 等(2003) 研究认为，震前ULF 辐射频谱坡减小(相应的分维数增加)的影响是确实存在的. 这一确定的动态过程，开始于震前的1—2个月，可归因于地震的发展孕育过程. 此时，震源区正处于自组织临界动态演变阶段. 在这种状态下，岩石介质呈蜂窝似的空隙状，其结构越来越具有脆性特征. 这种动态的分维过程的结果是，破裂过程中产生的ULF 电磁信号逐渐被高频的波动所补充.

由于台站使用的仪器没有明确的观测频段，不能确定增加的数据成分一定是高频率波动. 但是可以确定，两次地震前2—3个月，天水台绝对地磁观测 Z 分量时间序列 H 指数均减小到0.44，由持久性增强序列发展成为反持久性序列； 相关系数 C =-0.08，数据序列的增量间不相关； 数据的分维数 D =1.56达到最大值，这表明相对于其它时段，数据成分发生改变. 地震不是发生在 H 值极小的时候，而是发生 H 值转折回返的过程中或者回返以后，亦即地震不是发生在应力强度的极值点上，而是在极值点后. 尤明庆(2000) 研究岩石破裂实验发现，当岩样的细观强度分布比较均匀，岩石破裂的自组织临界点才接近峰值强度； 当岩样强度分布不太均匀、 屈服过程较长时，岩石破裂的自组织临界点则在峰值强度后取得. 这一动态过程与 Vallianatos和Tzanis(1999) 以及 Tzanis和Vallianatos(2002) 提出的，由于破裂扩展引起的裂隙变形过程和电荷的位移引起的地震前兆的物理模型理论相一致( Hayakawa et al ，1999 ).

从地质构造上讲，昆仑山口西 M _S 8.1地震区位于青藏亚板块的东昆仑构造带. 由于印度板块体不断向北推移造成与古亚洲板块多期碰撞闭合，青藏亚板块不仅表现为完整的巨大隆起，而且又是由多个构造块体拼合而成的复杂构造单元. 昆仑山口西 M _S 8.1地震是发生在近东西走向、 几乎直立的断层上的左旋运动，南盘相对于北盘向东运动( 许力生，陈运泰，2004 ). 地表破裂带全长426 km，总体走向90°—110°. 根据Harvard的震源参数，节面Ⅰ走向为94°， P 轴方位是56°. 由 图 1 可看出，天水台位于本次地震的东偏南方向，这表明天水台位于破裂带的延长线上，且近主应力的方向. 这可能是引起震前3个月地磁数据异常变化的原因之一.

天水地磁台位于汶川地震的NE方向，在Ⅵ裂度区内，接近于Ⅶ度区，距震中约450 km，与汶川地震同属于南—北地震构造带，在龙门山断裂带的延长线上. 本次地震发生在我国的南—北地震构造带的中段、 青藏高原东缘的NE向龙门山断裂带上( 图 1 ). 该断裂带显示出较强烈的、 由北西朝南东方向的逆冲运动，并兼有右旋走滑的水平运动分量. 根据汶川地震的震源机制解( 郑勇，2008 )， P 轴走向为102°、 倾角20°，表明主压应力轴为近EW向，压力近水平分布( 张学民等，2009 ). 这种应力变化通过压磁效应、 膨胀磁效应、 感应磁效应、 电动磁效应、 流变磁效应及热磁效应等，在地磁场长趋势变化产生随时间变化的局部前兆异常( 丁鉴海等，1994 ). 根据前人的岩石实验结果，岩石在受压过程中，当EW方向受到巨大压力时，EW，SN，NE和NW四个方向都可观测到不同程度的电磁异常，其中EW向最弱，NE向最为强烈. 这可能是引起天水绝对地磁观测 Z 分量异常的主要原因. 4 讨论与结论

通过上述分析，天水台绝对地磁观测 Z 分量数据序列变化是一个非线性动力学过程，也是一个由规律成分和随机成分叠加而成的时间序列过程. 大部分时间遵循分数布朗运动的性质. 对2001年昆仑山口西 M _S 8.1地震和2008年汶川 M _S 8.0地震，选取时间序列的 R / S 分析结果表明，震前2—3个月，Hurst指数均减小到极小值， H =0.44＜0.5； 分维数增加到极大值， D =1.56； 相关系数 C =-0.08<0. 表明观测数据成分发生了变化，改变了数据发展变化的趋势，由持久性序列发展为反持久性序列. 这可以看作是强震前的一种前兆信号.

两次地震前的2—3个月，天水台绝对地磁观测 Z 分量数据序列Hurst指数 H 减小、 分维数 D 增加. 这表明在强震前，地磁观测数据中有更高频率成分的波动出现，改变了观测数据自身的行为方式，从随机的和持久性的序列发展成为反持久序列( Gotoh et al ，2003 ). 在震中区可能也产生高频段的电磁信号(如，VLF和HF/VHF)，由于在地壳中的快速衰减，它们在地球表面不能被探测到. 所以，这些高频的电磁信号可能是在近地表的二次效应产生的，其产生的机制还不是很清楚.

在原始数据序列分段过程中也发现，每个子序列的数据个数改变时，相应的计算结果相差很大. 昆仑山口西地震数据分析是取30个点为一子序列，而汶川地震取40个点. 这说明，合理的分段，对于将地震前地壳内部应力变化引起的地磁异常显现出来是非常重要的.

到目前为止，还不能说用这一方法可以预报地震. 因为，一方面方法本身没有作处理结果检验和误差分析； 另一方面，选用的两个震例均为8级以上的特大地震，对于震级差异的震例数据分析是否有效还需要在实践中检验. 但是，它毕竟给我们提供了一种地震前兆信息提取的辅助手段.

Vallianatos F, Tzanis A. 1999. A model for the generation of precursory electric and magnetic fields associated with the deformation rate of the earthquake focus[G]//Hayakawa M ed. Atmospheric and Ionospheric Electromagnetic Phenomena Associated with Earthquakes. Tokyo: Terra Sci Publ Co: 287-306. ( 1) 向、 几乎直立的断层上的左旋运动，南盘相对于北盘向东运动( 许力生，陈运泰，2004 ). 地表破裂带全长426 km，总体走向90°—110°. 根据Harvard的震源参数， 了观测数据自身的行为方式，从随机的和持久性的序列发展成为反持久序列( Gotoh et al ，2003 ). 在震中区可能也产生高频段的电磁信号(如，VLF和HF/VHF)，由于在地壳中 过程与 Vallianatos和Tzanis(1999) 以及 Tzanis和Vallianatos(2002) 提出的，由于破裂扩展引起的裂隙变形过程和电荷的位移引起的地震前兆 程较长时，岩石破裂的自组织临界点则在峰值强度后取得. 这一动态过程与 Vallianatos和Tzanis(1999) 以及 Tzanis和Vallianatos(2002) 提出的，由 频(ULF)地磁观测数据，发现两次强震前，ULF地磁观测数据分形维数增加. 希腊 Vallianatos等(2004) 将 R / S 分析方法运用到时间序列的超低频地磁三分量(水平分量X和Y