石 磊 陈 石 蒋长胜 徐伟民 卢红艳 郭凤义

（中国北京 100081 中国地震局地球物理研究所）

引言

北京时间2013年4月20日8时2分46秒在四川省雅安市芦山县境内发生了MS7．0地震．截至5月15日该地震共造成196人死亡，1万1 470余人受伤，

震区发生大面积滑坡和建筑物损毁，造成重大经济损失．震源机制和破裂过程等研究结果表明，该地震是继2008年5月12日汶川MS8．0强地震后发生在龙门山断裂带上的又一次大震，为逆冲型地震（张勇等，2013）．

青藏高原东缘的龙门山断裂带北起青川，经北川、映秀、宝兴至泸定附近，全长近500km，宽约40—50km，走向北东（Zhangetal，2010）．龙门造山带位于松潘—甘孜印支造山带与扬子克拉通的中新生代（T3—Q）四川前陆盆地之间，以出露大面积的前震旦纪古老变质杂岩（700—800Ma）、少量的新元古代变质火山岩及震旦纪—早、中三叠世海相沉积岩为特征．龙门山东缘发育一系列的逆冲断裂和飞来峰构造，逆冲作用使山体向东逆冲推覆到四川盆地之上（许志琴等，2007）．该断裂带主要由3条近北东走向的主干断裂组成（图1）：江油—灌县断裂（前山断裂）、北川—映秀断裂（中央断裂）、茂县—汶川断裂（后山断裂）（张培震等，2008）．位于龙门山断裂带西侧的松潘—甘孜地块边缘或内部的断裂带中有少量古生界地层甚至前寒武系变质岩出露，三叠纪末的印支运动使该地区褶皱连山，强烈变形；东部的扬子地块是由前寒武纪变质基底和震旦纪—中生代沉积盖层组成的海相沉积，早、中侏罗世之后为陆相碎屑岩堆积，新生代地层主要分布在成都盆地（Chenetal，2000；楼海等，2008）．

由重力观测得到的重力异常包含了从地表到深部所有密度不均匀体引起的重力效应，信息非常丰富．因此，若要根据重力异常解释某个地质体，必须首先从叠加重力异常中分离出单纯由这个地质体引起的异常，然后对该异常进行分析解释．然而，重力异常分离目前仍是重力资料处理解释中没有很好解决的难题之一（曾华霖，2005）．目前比较具有代表性的重力异常分离方法有：传统向上延拓法、带通滤波法、小波变换多尺度分解法（侯遵泽，杨文采，1997；Fedi，Quarta，1998；杨文采等，2001）、熵滤波法（Nikitinetal，1984）、非线性滤波（Keating，Pinet，2011）、优选延拓法（Pawlowski，1995；Mengetal，2009）和优化滤波法（Guoetal，2013）等．

本文首先简单介绍其中5种重力异常分离方法（带通滤波法、传统向上延拓法、熵滤波法、非线性滤波法和优化滤波法）的基本原理，然后分别利用这5种方法对同一组理论模型数据进行试验，对比分析各种方法的异常分离效果，从中优选出效果较好的方法对芦山地震震区重力异常进行分离，初步分析认识此次地震发生的构造背景．

图1 研究区构造背景简图．图中F1，F2和F3分别代表龙门山后山断裂，中央断裂和前山断裂Fig．1 Tectonic sketch of the research area．F1：Longmenshan back-range fault；F2：Longmenshan central fault；F3：Longmenshan front-range fault

1 重力异常分离方法原理

1．1 带通滤波法

带通滤波法通过对重力异常指定频段的频谱进行滤波通过，对其它频段频谱进行阻止，从而分离得到指定频段的重力异常．所用的带通滤波频率响应函数为

其中，k为波数，k0为低切波数，k1为高切波数．

1．2 向上延拓法

向上延拓法的目的是压制浅层地质因素或干扰引起的异常，突出相对深部地质因素产生的重力异常．在一定范围内向上延拓的高度越大，所反映的地质信息越宏观，也相当于深度越大的．因此，经常通过向上延拓不同高度得到的延拓场，来研究不同深度的场源或构造信息．常用的频率域向上延拓法通过对重力异常频谱进行向上延拓滤波，得到延拓场即区域异常．所用的向上延拓频率响应函数为

其中，h为要延拓的高度．

1．3 熵滤波法

熵滤波法是前苏联学者Nikitin等（1984）在熵的理论基础上提出的一种重力异常分离方法，是一种空间域滑动窗口滤波方法．设网格化后规则重力数据中的任意待滤波点p的异常为fp，熵滤波法所用的滑动窗口的中心为待滤波点，窗口大小自定义，落入滑动窗口内的异常为fi（i＝1，2，…，N），则待滤波点p熵滤波后的区域异常为

其中，ei为各点的熵值，wi为各点的熵权系数．可见，当滑动窗口内某点异常的熵值越小，表明它的变异程度越大，在熵滤波中所起的作用越大，其熵权重值则越大；反之，某点异常的熵值越大，表明它的变异程度越小，在熵滤波中所起的作用越小，其熵权重值则越小．

1．4 非线性滤波法

非线性滤波法是Keating和Pinet（2011）在Naudy和Dreyer（1968）图解法基础上提出的一种重力异常分离方法．在该方法中剖面重力异常被分为两种类型：有一个最大值的简单异常；有一个最大值和两个最小值的复杂异常．

对于网格化后的规则重力数据，定义第二半差值为

其中，C为待滤波点的已知异常值；A，B，D和E分别为待滤波点前后第一节点和第二节点的已知异常值．

异常的形状由T与S的比值R来决定，其中

如果R为负值，说明异常有一个最大值和两个最小值，此时对待滤波点进行非线性滤波，滤波后的异常值为Q＝2（B＋D）／3－（A＋E）／6；如果0≤R≤2，则异常有一个最大值，此时对待滤波点进行非线性滤波，滤波后的异常值为Q＝0．5C＋0．25（B＋D）；而如果R＞2，异常过于平缓，则不对待滤波点进行非线性滤波．

对于面积性的规则数据，其非线性滤波首先沿x方向对各剖面数据进行上述的剖面非线性滤波，得到第一次滤波后的数据；然后再沿y方向对各剖面x方向滤波后的数据重新进行剖面非线性滤波，得到第二次滤波后的数据；最后对x和y方向滤波后的数据做巴特沃斯滤波，滤去高频干扰，得到最终的滤波后的数据，即非线性滤波得到的区域异常．

1．5 优化滤波法

针对传统向上延拓方法在向上延拓时异常分离不彻底问题，意大利学者Pawlowski（1995）根据维纳滤波和格林等效层原理提出了重力异常分离的优选延拓法．应用该方法进行向上延拓时，可以在压制浅源短波长信息的同时，保持深源长波长信息不衰减．之后，许德树和曾华霖（2000）、Meng等（2009）在优选延拓法的基础上给出了基于优选向上延拓的延拓差值法，实现对某一指定频段的局部异常有选择地提取或分离，即类似于带通滤波．郭良辉等（2010）进一步改进优选延拓法而提出优化滤波法，用于分离指定频段的重力异常．它既保留了优选延拓法分离异常较彻底的优点，同时其分离异常不需要已知延拓高度，更能满足实际数据处理需求．

实测重力异常的傅里叶功率谱可通过引入格林等效层概念来模拟，即利用处于不同深度的m＋n层格林等效源薄层来建立实测重力异常的径向平均对数功率谱模型．假设目标层为深部i→j（n＋m≥j≥i≥1）等效层，则优化滤波法的期望是保持目标层场源信息不变，而对其它层场源信息作压制．设维纳滤波器的输入是实测异常g实测（x，y），期望输出是g优化（x，y），则维纳滤波器的频率响应（Wiener，1949；Clarke，1969）为

其中，G代表频率域的频谱，＜＞代表数学期望，＊代表复共轭．式中的分子是滤波器期望输出和输入信号（实测异常）的互功率频谱密度函数，分母是滤波器输入的功率频谱密度函数．

假设目标层与其它深度层场源信息互不相关，则式（6）的维纳滤波转换函数可进一步改写为

我们称式（7）的维纳滤波转换函数为优化滤波算子．当i＝1时，i→j即为1→j，此时的优化滤波是分离1→j层场源信息的低通滤波（图2）；当j＝m＋n时，i→j即为i→m＋n，此时的优化滤波是分离i→m＋n层场源信息的高通滤波；当n＋m＞j≥i＞1时，此时的优化滤波是分离i→j层场源信息的带通滤波．因此，利用式（7）可对实测重力异常进行不同形式的优化滤波，分离出目标层场源在原始观测面上的异常信息．

图2 （a）理论模型各长方体水平面位置分布；（b）理论重力异常等值线图 （单位：10－5 m／s2）；（c）理论含噪重力异常等值线图（单位：10－5 m／s2）Fig．2 （a）Projection of each prism of synthesized model on the xoysurface；（b）Contour map of theoretical gravity anomalies；（c）Contour map of noised gravity anomalies．The unit of gravity anomaly is 10－5 m／s2

2 理论模型数据试验对比

理论模型由处于3个不同深度层、不同大小和剩余密度值的11个直立长方体组合而成．各长方体的几何参数和剩余密度如表1所示．图2a给出了各直立长方体在水平面上的投影分布．假设深层长方体A1，A2和A3产生区域异常，中层长方体B1，B2，B3，B4，B5和浅层长方体C1，C2，C3共同产生局部异常（或剩余异常）．为检验异常分离效果，深层长方体A1，A2和A3距离较近，中、浅层长方体分布分散且部分叠置于深层长方体上，产生叠加异常．数值正演该模型在深度为0km处、测网为101×101、纵横向网格间距均为0．2km的平面网格上的重力异常如图2b所示．对理论重力异常加入标准差为0．2×10－5m／s2的随机高斯噪声，结果见图2c．

图3 5种重力异常分离方法异常分离结果图（单位：10－5 m／s2）Fig．3 Gravity anomaly（in unit of 10－5 m／s2）separation results by five filtering methods

表1 理论模型各长方体的几何参数和剩余密度Table 1 Geometric parameters and residual density of each prism in the synthesized model

下面分别利用带通滤波法、向上延拓法、熵滤波法、非线性滤波法和优化滤波法对理论含噪重力异常进行去除高斯噪声、分离区域异常和分离B，C层局部异常的异常分离试验．对含噪重力异常计算径向对数功率谱，用0—0．274 5cycles／km（频段1）、0．274 5—1．215 7cycles／km（频段2）和1．215 7—5cycles／km（频段3）等3个等效层分段拟合．

在进行去除高斯噪声干扰时，带通滤波法的滤波波长宽度为800m，向上延拓高度为300m；熵滤波法滤波窗口为600m；非线性滤波法的滤波波长宽度为800m；优化滤波法以频段3为浅源信息、频段1和2为深源信息进行滤波，结果见图3a．对比上述各图可见，带通滤波法、熵滤波法、非线性滤波法和优化滤波法均能有效地压制高斯噪声干扰，除了C层有效信号保留得不是很好外（因为其频率与噪声较接近），A和B层的有效信号都得到了较好的保留；而向上延拓法在压制高斯噪声干扰的同时，也部分压制了有效信号．

在进行分离区域异常时，带通滤波法的滤波波长为4 000m，向上延拓高度为1 000m；熵滤波法滤波窗口为3 000m；非线性滤波法的滤波波长为4 000m；优化滤波法以频段2和3为浅源信息、频段1为深源信息进行滤波，所得结果如图3b所示．通过对比分析可知，带通滤波法、向上延拓法和熵滤波法在有效压制高频噪声干扰及B，C层局部异常的同时，也部分压制了A层的区域异常，尤其是向上延拓法的区域异常被压制最多；非线性滤波法较好地保留了区域异常，但B层的局部异常压制不够彻底，部分遗留在区域异常上；优化滤波法不仅有效地压制了高斯噪声干扰及B，C层的局部异常，A层的有效信号也得到了较完整的保留．

分别利用各种方法去除高斯噪声后的异常减去区域异常，即可得到分离B，C层的局部异常，所得结果如图3c所示．通过对比上述各图可知，带通滤波法、向上延拓法和熵滤波法都能有效分离出B，C层的局部异常，但也遗留了A层的部分区域异常，尤其是向上延拓法的区域异常遗留较多；非线性滤波法分离结果没有明显的区域异常，但B层的局部异常被部分压制；优化滤波法则有效压制了高斯噪声及A层的区域异常，B，C层的有效信息都被较好地保留了下来．

通过上述理论模型数据试验对比分析得出，优化滤波法异常分离效果最好，非线性滤波法次之，带通滤波法和熵滤波法再次之，向上延拓法分离效果最差．因此，下面我们将利用优化滤波法对芦山地震震区布格重力异常进行处理，以初步分析认识此次地震发生的构造背景．

3 芦山地震震区重力异常分离与特征分析

本文采用的布格重力异常数据来自中国地质调查局，比例为1∶100万，网格化点、线间距均为5km（图4）．研究区内布格重力异常均为负值，异常值由南东—北西逐渐减小，异常幅度为－465×10－5—－95×10－5m／s2．龙门山断裂带为研究区内一条显著的重力梯度带，近北东走向．其布格重力异常值由成都的－150×10－5m／s2向西到马尔康的－390×10－5m／s2，最大梯度变化达2．5（×10－5m／s2）／km．该重力梯度带在天全附近分为两支，西支与围绕青藏高原北部和东部的巨大重力梯度带重合，向西南方向延伸；东支则向东南方向伸展，至遵义附近．此次芦山MS7．0地震接近研究区内龙门山重力梯度带分叉处，且MS≥3．0余震呈近北东向展布，与梯度带走向基本一致．

图4 研究区布格重力异常图Fig．4 Bouguer gravity anomaly in the research area

龙门山构造带东西两侧，重力异常变化逐渐平缓．扬子地块重力异常等值线以北东向走向为主，在－95×10－5—－170×10－5m／s2范围内变化．内江—江津一带具有一定规模的局部高重力异常，推测与莫霍面隆起有关（Zhangetal，2010）．龙门山西侧松潘—甘孜地块重力异常等值线走向主要为近北北西向，异常幅值为－220×10－5—－460×10－5m／s2，变化较为平缓，且沿北北西向有扭曲现象，显示出青藏高原隆起过程中经历了推覆、挤压和走滑等作用，物质自西向东流动．

通过上述理论模型试验对比分析，本文将采用异常分离效果较理想的优化滤波法对芦山震区及邻区布格重力异常进行分离，并结合其它处理结果进行简单分析．使用优化滤波法进行异常分离，关键是径向对数功率谱的分段拟合．根据实际布格重力异常的径向对数功率谱形状特征（图5），我们采用0—0．0016cycles／km（频段1）、0．001 6—0．005 5cycles／km（频段2）、0．005 5—0．014 5cycles／km（频段3）、0．014 5—0．039 5cycles／km（频段4）、0．039 5—0．059 9cycles／km（频段5）和0．059 9—0．1cycles／km（频段6）等6个等效层分段拟合．将频段1—2对应的深部场源视为区域异常，频段3—6对应的细节之和视为局部异常．

图5 研究区布格重力异常径向对数功率谱及其分段拟合黑色实线为研究区布格重力异常功率谱，蓝色实线为分频段拟合，红色实线为拟合的功率谱Fig．5 Radial logarithm power spectrum of the Bouguer gravity anomaly in the research area and equivalent layer fitting Black curve denotes power spectrum of Bouguer gravity anomaly，blue curve denotes the fitted straight lines in each radial frequency band，and the red curve denotes the fitted radial logarithm power spectrum by using power spectrum density function models

图6给出了采用优化滤波法分离研究区布格重力异常所得到的剩余重力异常，异常幅值在－20×10－5—20×10－5m／s2之间，表明其主要显示地壳中浅部构造特征．该图中龙门山断裂带反映为与地表断裂走向基本一致的高重力异常带，认为该断裂带重力异常主要是由地壳中浅部密度结构不均匀引起的．天全北部附近，沿龙门山断裂带走向的高重力异常带被另一条北北西向的高重力异常带错断，芦山地震震区即位于该异常带连贯性薄弱地段东侧，异常值剧烈变化地区．扬子地块内部在大面积负重力异常上分布着3条走向近似平行的宽缓高重力异常带，异常形式较为单一，表现出稳定块体特征．研究区内松潘—甘孜地块高重力异常带由北部的北北西走向，往南快速变为北西走向，显示青藏高原东缘持续向四川盆地的挤压作用．

图6 研究区剩余重力异常图Fig．6 Residual Bouguer gravity anomaly of the research area

图7给出了优化滤波法分离得到的研究区区域重力异常．其分布形态与布格重力异常的分布形态类似，但主要反映了地壳深部及上地幔顶部场源信息．由图7可见，松潘—甘孜地块与扬子地块分别位于近北东向的巨大重力梯度带两侧，其中红色实线为梯度带变化最为剧烈的地方．同样是在天全北部地区，该重力梯度带分为两支，西支向西南方向延伸，东支则向东南方向伸展，与浅部构造特征存在较大差异．由此推测天全以北附近地区构造复杂，地壳深浅部耦合较差．从图中我们还可以看出，研究区深部块体边界与地表大地构造单元界线位置不同，说明龙门山推覆体下的莫霍面为一较陡的斜坡．

本文还采用优化滤波法提取了研究区基底面深度重力异常．以前人所做的深地震测深结果为约束（崔作舟等，1996；王有学等，2005；朱介寿，2008；熊小松，2010；Wardetal，2013），利用频率域密度界面反演方法得到了芦山地震震区及邻区基底深度分布图（图8）．从图8中可以看出，龙门山断裂带为基底隆起区，基底深度7km左右．该断裂带以东的扬子地块普遍处于基底隆起区，沉积盖层较薄；而断裂带西部的松潘—甘孜地块则位于基底凹陷区，马尔康—康定一线最深达到12km，西北方向逐渐减薄．与汶川地震震中绝大多数分布在基底的隆起区不同（Zhangetal，2010），芦山地震震中主要分布于基底的局部凹陷区，这方面还需要我们进一步深入研究．

图7 研究区区域重力异常图Fig．7 Regional Bouguer gravity anomaly of the research area

4 讨论与结论

本文首先介绍了5种重力异常分离方法（带通滤波法、向上延拓法、熵滤波法、非线性滤波法和优化滤波法）的基本原理，然后通过同一组理论模型数据试验，对比分析5种异常分离方法的结果，从中优选出效果较好的异常分离方法对芦山地震震区布格重力异常进行处理，并初步分析了研究区重力异常特征，得到结论如下：

1）通过对同一组理论模型数据进行异常分离试验对比，得出5种重力异常分离方法中优化滤波法异常分离效果最好，可以有效地提取不同频段的重力异常，并保证有效信号的完整保留；非线性滤波法分离效果次之，可以较好地分离不同深度的重力异常，但部分有效信号稍有损失；带通滤波法和熵滤波法再次之，分离不同深度重力异常的同时部分削减有效信号信息；传统向上延拓法异常分离效果最差，分离重力异常的同时有效信号损失最为严重．

2）研究区内布格重力异常均为负值，异常值由南东—北西逐渐减小．龙门山断裂带为区内一条显著的重力梯度带，近北东走向．该重力梯度带在天全附近分为 两支，西支与围绕青藏高原北部和东部的巨大重力梯度带重合，向西南方向延伸；东支向东南方向伸展至遵义附近．此次芦山地震震区接近研究区内龙门山重力梯度带分叉处，且MS≥3．0余震呈近北东向展布，与梯度带走向基本一致．说明此次强地震的发生与重力梯度带下方深部结构和构造活动密切相关．

图8 研究区基底深度分布图Fig．8 Distribution of basement depths in the research area

3）研究区剩余重力异常中，龙门山断裂带反映为与地表断裂走向基本一致的高重力异常带，认为该断裂带重力异常主要是由地壳中浅部密度结构不均匀引起的．天全北部附近地区，沿龙门山断裂带走向的高重力异常带被另一条北北西向的高重力异常带错断，芦山地震震区即位于该异常带连贯性薄弱地段东侧的异常值剧烈变化地区．同样是在天全以北地区，区域重力异常梯度带分为两支，西支向西南方向延伸，而东支则向东南方向伸展，与浅部构造特征存在较大差异．由此推测芦山地震震区构造复杂，地壳深浅部耦合较差，此次强地震的发生受浅部和深部构造的共同控制．

4）研究区龙门山断裂带为基底隆起区，基底深度为7km左右．断裂带以东的扬子地块普遍处于基底隆起区，沉积盖层较薄；而断裂带西部的松潘—甘孜地块则位于基底凹陷区，马尔康—康定一线最深达到12km．与汶川地震震中绝大多数分布在基底的隆起区不同，芦山地震震中主要分布于基底的局部凹陷区，这方面还需要我们进一步深入研究．

四川省地震局为本研究提供了芦山地震序列速报目录；在论文评审过程中，审稿专家对本文提出了宝贵的建议．作者在此一并表示感谢．

崔作舟，陈纪平，吴苓．1996．花石峡—邵阳深部地壳的结构和构造［M］．北京：地质出版社：1-192．

郭良辉，孟小红，石磊．2010．磁异常ΔT三维相关成像［J］．地球物理学报，53（2）：435-441．

侯遵泽，杨文采．1997．中国重力异常的小波变换与多尺度分析［J］．地球物理学报，40（1）：85-95．

楼海，王椿镛，吕志勇，姚志祥，戴仕贵，尤慧川．2008．2008年汶川MS8．0级地震的深部构造环境：远震P波接收函数和布格重力异常的联合解释［J］．中国科学：D辑，38（10）：1207-1220．

王有学，Mooney W D，韩果花，袁学诚，姜枚．2005．台湾—阿尔泰地学断面阿尔金—龙门山剖面的地壳纵波速度结构［J］．地球物理学报，48（1）：98-106．

熊小松．2010．中国大陆莫霍面深度与变化特征及其地球动力学意义［D］．北京：中国地质科学院：1-200．

许德树，曾华霖．2000．优选延拓技术及其在中国布格重力异常图处理上的应用［J］．现代地质，14（2）：215-222．

许志琴，李化启，侯立炜，付小芳，陈文，曾令森，蔡志慧，陈方远．2007．青藏高原东缘龙门—锦屏造山带的崛起：大型拆离断层和挤出机制［J］．地质通报，26（10）：1262-1267．

杨文采，施志群，侯遵泽．2001．离散小波变换与重力异常多重分界［J］．地球物理学报，44（4）：534-541．

曾华霖．2005．重力场与重力勘探［M］．北京：地质出版社：189-234．

张培震，徐锡伟，闻学泽，冉勇康．2008．2008年汶川8．0级地震发震断裂的滑动速率、复发周期和构造成因［J］．地球物理学报，51（4）：1066-1073．

张勇，许力生，陈运泰．2013．芦山4·20地震破裂过程及其致灾特征初步分析［J］．地球物理学报，56（4）：1408-1411．

朱介寿．2008．汶川地震的岩石圈深部结构与动力学背景［J］．成都理工大学学报：自然科学版，35（4）：348-356．

Clarke G K C．1969．Optimum second-derivative and downward-continuation filters［J］．Geophysics，34：424-437．

Chen Z，Burchfiel B C，Liu Y，King R W，Royden L H，Tang W，Wang E，Zhao J，Zhang X．2000．Global positioning system measurements from eastern Tibet and their implications for India／Eurasia intercontinental deformation［J］．J GeophysRes，105：16215-16227．

Fedi M，Quarta T．1998．Wavelet analysis of the regional-residual and local separation of potential field anomalies［J］．GeophysicalProspecting，46：507-525．

Guo L H，Meng X H，Chen Z X，Li S L，Zheng Y M．2013．Preferential filtering for gravity anomaly separation［J］．Computers＆Geosciences，51：247-254．

Keating P，Pinet N．2011．Use of non-linear filtering for the regional-residual separation of potential field data［J］．JApplGeophys，19：315-322．

Meng X H，Guo L H，Chen Z X，Li S L，Shi L．2009．A method for gravity anomaly separation based on preferential continuation and its application［J］．AppliedGeophysics，6（3）：217-225．

Naudy H，Dreyer H．1968．Essai de filtrage non-linéaire appliquéaux profils aéromagnétiques［J］．GeophysicalProspecting，16（2）：171-178．

Nikitin A A，Vasov O K，Belov A P．1984．Vozmozhnosti kompleksnoy geofizicheskoy interpretatsii na baze entropiynogo fil＇tra．Izvestiya Akademii Nauk Turkmenskoy SSR［J］．SeriyaFiziko-Tekhnicheskikh，KhimicheskikhiGeologicheskikhNauk，2：79-82．

Pawlowski R S．1995．Preferential continuation for potential-field anomaly enhancement［J］．Geophysics，60（2）：390-398．

Ward S，Mikhail K，Fred B，Magdala T，Walter D M，Sierd C．2013．High resolution regional crustal models from irregularly distributed data：Application to Asia and adjacent areas［J］．Tectonophysics，602：55-68．

Wiener N．1949．Extrapolation，Interpolation，andSmoothingofStationaryTimeSeries［M］．New York：John Wiley＆Sons Inc：1-24．

Zhang J S，Gao R，Zeng L S，Li Q S，Guan Y，He R Z，Wang H Y，Lu Z W．2010．Relationship between characteristics of gravity and magnetic anomalies and the earthquakes in the Longmenshan range and adjacent areas［J］．Tectonophysics，491（1／4）：218-229．