重力变化与地壳和地幔内部的物质运移有关，而大地震的孕育和发生也伴随着地球内部物质的运移，所以提取及分析重力变化信息有助于研究大地震的孕育和发生机理。早在1964年日本新潟地震、1968年新西兰因南格华地震及1971年美国圣费尔南多地震前后都观测到重力异常变化，我国的海城地震和唐山地震前后也出现了与地震有关的重力变化^[1]。对于近年来中国大陆发生的一些强震，也有震前重力异常的研究报告。如2013年的四川芦山地震，基于流动重力观测的研究结果表明，震前康定和小金出现重力负异常区，四川盆地出现重力正异常区，震中位于正、负异常区交汇的高梯度带上^[2]。又如2014年云南鲁甸地震，根据西昌测网的流动观测数据，地震发生在重力变化梯级带方向由NE向转为近EW向的转折部位，这种重力异常变化与地壳内部介质孔隙流体的运移相关^[3]。关于2013年岷县-漳县M_S6.6地震也有重力方面的研究。祝意青等^[4]认为，震中位置处于NE向的重力变化高梯度带，且这些异常变化是深部地下流体变化的响应。本研究采用距震中较近的2个台(兰州、高台)的重力观测资料，试图提取岷县-漳县M_S6.6地震的震前异常信息，并对重力异常的形成机理进行初步探讨。

1 资料与资料处理

去趋势波动分析(deterended fluctuation analysis，DFA)方法最早应用于医学领域，对类似心电图信号进行分析处理，用于诊断人体的健康状态。和一般的分析方法相比较，其优点是可以在非稳定时序信号中检测出嵌入的长程交互作用，并避免一些外在非稳定信号的干扰^[5]。也有学者将去趋势波动分析方法应用于地电、地磁资料处理，结果表明，分形标度异常与地震火山活动具有较好的关联性^[6-8]。荣扬名等^[9]用该方法处理2011年日本东北大地震前3个近台的地磁资料发现，特低频信号分形标度特征指标在震前25~30 d左右出现同步的异常增加，认为是与这次地震有关联的异常信息。刘立炜等^[10]将该方法应用于祁连山断裂和海原断裂带的跨断层水准资料分析，结果表明，去趋势波动分析对周围4级以上地震显示出折点异常。和其他领域相比，去趋势波动分析在地学领域应用的报道较少，本研究尝试将该方法应用于重力观测数据中。

上世纪90年代，为了分析一些生理信号的分形标度特征，Peng等^[11]提出去趋势波动分析方法。假设得到的时间序列的数据长度为N, 将其分为k段(或k个盒子)，数据序列为B(1), B(2), B(3), …, B(i)，则累计时间序列可表示为：

$ y(k) = \sum\limits_{i = 1}^k {\left[ {B(i) - {B_{{\rm{ave}}}}} \right]} $

(1)

式中，B_ave为该段的均值。

$ F(n) = \sqrt {\frac{1}{N}\sum\limits_{k = 1}^N {{{\left[ {y(k) - {y_n}(k)} \right]}^2}} } $

(2)

式中，n为每段数据个数，各段的平均波动用F(n)表示。logF(n)和logn对应的斜率为α，称为标度指数(分形指数)，用于描述波动性特征。当0 < α < 0.5时，时间序列负相关，是一种反持久性状态；当0.5 < α < 1时，未来变化趋势与其历史变化趋势呈正相关，具有记忆效应和分形特征。

选取兰州、高台2个台2011~2015年的重力观测资料进行分析，这2个台分别距震中约190 km和710 km, 2个台相距约490 km。数据采样率为每秒1个样点。首先对个别缺数进行线性内插，保证数据在时序上的完整性。将每天数据形成一个数据文件，采用去趋势波动分析进行处理。主要计算步骤为：1)本项计算确定数据最小段数为4，由于重点关注100~200 s的周期范围，每段最小数据点数为100，最大数据点数为200，并以此得到数据分段数k；2)计算每段的去均值累加和得到一个相应序列；3)对上述序列进行多项式拟合，得到拟合曲线值；4)据此求得所有段的平均波动值F(n)；5)对F(n)和n取双对数再相除，每天得到一个分形指数(α值)。分形指数在时序上经过光滑处理后，结果如图 1所示。

2 结果与讨论

如图 1所示，2个台的重力分形指数的时序变化整体较为相似。本文还不能以确定的物理模型定量给出异常阈值，考虑到台站之间背景值的差异性，这里经验性地给出9.4倍的标准差值作为重力分形指数的异常阈值，高台台和兰州台的异常阈值分别是0.68和0.53，将该台分形指数最小值和阈值之间视为常态变化, 则这2个台的常态变化分别是0.36~0.68、0.3~0.53。结果显示，在5~8月间都存在一个异常升高的过程，兰州台在2013年5月初重力分形指数上升超过阈值，大约在7月初达到峰值0.59。高台台重力分形指数则是从5月下旬上升超过阈值，在7月中旬达到峰值0.71，峰值出现时间略微滞后于兰州台。2013-07-22距离兰州约190 km的岷县与漳县交界处发生M_S6.6地震。高台台和兰州台之间的距离约为490 km，高台台距离震中约710 km。在这样一个空间跨度上，2个台的重力分形指数的时序变化显示大致相似的特征，意味着这种主体变化不应是来自地壳浅部，更多地反映了地壳深部物质密度的变化。

关于重力变化的机理，国内外学者提出了许多模式进行解释。这些模式可概括为2种观点：1)有关区域的地壳物质受到挤压或拉伸，密度发生变化；2)外部流体的迁移引起重力变化^[12-14]。近年来，关于青藏高原及周缘板块运动和动力学问题成为世界地学领域的热点问题之一。虽然关于下地壳物质是否流动尚有争议^[15]，但不论是否存在热物质的流动，不同学科关于这一区域地下结构的研究结果还是给我们提供了许多新的认识。滕吉文等^[16]认为，汶川地震是地壳深部物质与能量强烈交换和物质相向(上、下)运动的结果。对于甘肃东南部(岷县-漳县地震的发震区域)，一条经过岷县震区、西南-北东向展布的320 km测线的电磁测深结果表明，震源处于上部高阻体和下部低阻体的陡变位置，在研究区域内中下地壳广泛存在着低阻层，深度约为10~40 km, 根据剖面图推测，走向应该是北西方向。尽管这次地震发生机理尚不清晰，但能量的聚集和地震的发生一定与低阻层有关^[17]。在兰州以西这一区域，一条经过古浪、北东19°方向展布的340 km测线的电磁测深结果表明，约10~40 km深度存在北西向展布、成层性较好的低阻层，几条大断裂都收敛于低阻层^[18]。由此可知，兰州以东及以西的河西地区在中下地壳存在低阻层。祝意青等^[19]基于GPS测量和流动重力观测给出2011~2014年青藏高原东北缘的面膨胀率分布及重力分布，结果显示，高台台处于面收缩率较大区域，兰州台处于面收缩率较小的区域，震中位置处于2个面收缩率峰值区之间的过渡带，也是在重力变化高梯度带上，研究范围内的重力变化是该地区深部壳、幔物质运移引起的重力变化效应。GPS测量的是地表形变程度，这3个位置之间的地表形变状态存在很大的差异性, 而兰州台和高台台的重力分形指数时序图却显示出相似的变化特征。由此认为，本文研究得到的结果更多的是反映了地壳中下层(或者低阻层)乃至地幔的物质(或能量)运移变化。2个台在震前出现的重力分形指数升高意味着物质(或能量)运移引起重力持续记忆性的变化，具有越来越强的正相关特征，是对能量汇集过程的动态描述。这种异常的过程，应该与2013-07-22发生的岷县-漳县地震有关联。

3 结语

通过对兰州、高台2个台2011~2015年的重力观测资料进行去趋势波动性分析，得到这2个台的重力分形指数时序图。结合其他学者的研究结果，初步得到以下认识：

1) 2个台的重力分形指数的时序变化整体较为相似。结合其他学者的研究报告认为，这种重力变化更多的是反映了地壳中下层(或者低阻层)乃至地幔的物质运移变化。

2) 从2013年5月开始到7月下旬，2个台的重力分形指数异常上升，应该与岷县-漳县M_S6.6地震有关联。本项研究得到的异常阈值为9.4倍的重力分形指数标准差是一个考虑了台站背景差异的经验性指标，是否具有普适性尚需进一步研究。

本文是去趋势波动性分析应用于大地震前异常信息提取的初步性、探索性研究，结果表明，该方法是可靠、有效的。期望后续得到大量数据及震例进行更深入、更广泛的研究，分形指数有望成为一个研究地震孕育和发生的标志性参量。