2018, Vol. 38
2. 北京劳雷物理探测仪器有限公司(上海), 上海市龙吴路777号, 200232
2017-08-08四川九寨沟发生7.0级地震，震中位置103.82°E、33.2°N，震源深度约20 km。该地震位于巴颜喀拉地块东部，属松潘-甘孜造山带与西秦岭造山带结合部位，区域地质构造主要包括西秦岭地槽褶皱带、龙门山推覆断褶带和川西北断块等构造单元，平均海拔高度在4 km以上。深部探测成果显示，震中位于布格重力异常梯度带，同时也是莫霍面陡变带，在巴颜喀拉块体地壳厚度为57~64 km，而毗邻该块体东部的四川盆地地壳厚度为40~45 km[1]。近10 a来，巴颜喀拉块体边缘地壳运动活跃，已发生包括2008年汶川MS8.0、2013年芦山MS7.0、2017年九寨沟MS7.0等多次强震。目前尚未有关此次地震震中区域物性变化特征的分析研究。
研究表明，动态重力场变化能够综合反映地下物质运移过程[2-4], 而地壳三维密度变化能够给出地震震中区物性变化及物质运移的方向[5]，揭示地壳运动的深层过程，有利于建立地震发震构造模型[6]。多年来，九寨沟地区已经积累了大量重力复测资料，为研究本次地震孕震过程中场源物性和结构的变化机理提供了数据基础。
1 数据与方法 1.1 数据九寨沟地震区域流动重力网(97~111°E，31~41°N)(图 1)每年进行2期观测，实行绝对重力控制下的相对重力联测模式，具备监测区内MS5.0以上地震的能力[7]。对比发现，2014-04~2016-04、2014-04~2017-04两期重力场累计变化数据与本次地震相关性最大，故选取这两期重力动态变化数据进行研究。以绝对重力基准为控制对数据进行整体平差[8]，各期点值平均精度优于15 μGal(图 2，红星为九寨沟地震震中)。
|
图 1 南北地震带区域流动重力网 Fig. 1 The gravity surveying route of the north-south earthquake belt in China |
图 2(a)为九寨沟地震前2014-04~2016-04重力变化的累积结果，其变化范围在-150~90 μGal之间。图中西北方向，沿张掖-西宁一带显示出重力变化梯度带，负值区(约-120 μGal)位于西南部，正值区(约90 μGal)位于东北部，该现象应与印度块体对东构造结挤压增强有关。图中东北方向，沿银川-延安一带显示出重力变化梯度带，负值区(约-60 μGal)位于东北部，正值区(约90 μGal)位于西南部，其他地区基本不变。
|
图 2 累积重力场动态变化 Fig. 2 Dynamic image of gravity field |
图 2(b)为九寨沟地震前2014-04~2017-04重力变化的累积结果，正负重力变化区分较为明显。以震中西北方向为界(沿张掖-西宁-九寨沟方向)，在青藏高原东北缘呈现梯度带，重力变化范围在-90~120 μGal之间，在震中的西北部出现负变化(-90~0 μGal)。
对比两期重力变化可以看出，巴颜喀拉块体东北部物质在逐渐向西南方向移动，使得在巴颜喀拉和柴达木活动块体交界处形成最大超过-90 μGal的负重力差异(2014-04~2016-04)，并在震前转移、集中到震中及邻区附近(2014-04~2017-04)，使玛沁和马尔康一带显示出最大超过-60 μGal的负重力变化。
1.2 反演算法紧凑约束反演算法能够反映出地下密度变化体的连续性和局部性，为强震分析预报提供重要的判断依据[5]。本文采用紧凑约束反演方法[9-10]，并在反演过程中使用深度加权进行约束[11]。
将研究区地下划分为大小相同的块体组合，其中第i个点的重力变化可以表示为：
| $ \begin{array}{l} {g_i} = \sum\limits_{j = 1}^M {{A_{ij}}{m_j} + {e_j}} \\ i = 1, 2, \ldots N, {\rm{ }}\;\;j = 1, 2, \ldots M \end{array} $ | (1) |
式中，mj为第j个块体的密度，ej为第j个数据点的误差，Aij表示核函数。写成矩阵形式：
| $ \mathit{\boldsymbol{g}}{\rm{ }} = {\rm{ }}\mathit{\boldsymbol{Am}}{\rm{ }} + {\rm{ }}\mathit{\boldsymbol{e}} $ | (2) |
综合考虑密度分布与观测误差的最小准则：
| $ \sum\limits_{j = 1}^M {{f_\rho }(\Delta {\rho _j})} + \sum\limits_{{\rm{ }}i = 1}^N {{f_e}({e_i})} {\rm{ }} \to {\rm{min}} $ | (3) |
采用最小二乘约束下的迭代算法，其解为：
| $ \mathit{\boldsymbol{m}}{\rm{ }} = {\rm{ }}\mathit{\boldsymbol{W}}_m^{ - 1}\mathit{\boldsymbol{A}}{^{\rm{T}}}{({\rm{ }}\mathit{\boldsymbol{AW}}_m^{ - 1}\mathit{\boldsymbol{A}}{^{\rm{T}}} + \mu {\rm{ }}\mathit{\boldsymbol{W}}_e^{ - 1}\mathit{\boldsymbol{}})^{ - 1}}\delta {\rm{ }}\mathit{\boldsymbol{g}} $ | (4) |
式中，第k次迭代的表达式为：
| $ \begin{array}{l} \delta {\rm{ }}\mathit{\boldsymbol{m}}{^{k + 1}} = {\rm{ }}\mathit{\boldsymbol{W}}_{m(k)}^{ - 1}\mathit{\boldsymbol{A}}{^{\rm{T}}}{({\rm{ }}\mathit{\boldsymbol{AA}}{^{\rm{T}}} + \mu )^{ - 1}}\delta {\rm{ }}\mathit{\boldsymbol{g}}{^k}\\ \mathit{\boldsymbol{m}}{^{k + 1}} = {\rm{ }}\mathit{\boldsymbol{m}}{^k} + \delta {\rm{ }}\mathit{\boldsymbol{m}}{^{k + 1}} \end{array} $ | (5) |
式中，Wm是模型的紧凑度权重矩阵：
| $ \mathit{\boldsymbol{W}}{_{m(k)}} = {\rm{diag}}({\rm{ }}\mathit{\boldsymbol{W}}_m^{ - 1}) $ | (6) |
式中，Wmj－1=|mj|α+ε，这里选择ε为无穷小量，α=2为不变值[10]。We是第k次迭代的噪声权重矩阵，
在重力场反演过程中，反演结果会随深度而逐渐衰减。为减小这一效应，对其进行深度加权：
| $ {({W_z})_{jj}} = ({\sum\limits_{i = 1}^M {{\rm{ }}({z_j} + {h_i}))} ^\beta } $ | (7) |
式中，zj是第j个块体的平均深度，hi是可调整参数，用来控制核函数随深度的衰减。由于重力效应随指数平方衰减，选择β=2[11]。为缩小参数空间范围，减少反演问题的多解性以及改善迭代过程中的收敛性，对目标函数进行密度变化界限约束，增加密度变化的上下界条件。
2 三维重力反演结果 2.1 反演模型及参数构建引起重力变化的主要因素包括地下物质运移(物质密度变化)、地下水储量变化和地球形变。在研究区域，地球形变以及水量变化引起的重力变化量较小[12]，相对于重力变化的量级可以忽略，可认为重力变化主要是由地下物质运移所致。
将研究区域(97~111°E，31~41°N)地下划分为5层，厚度均为10 km；再将每层划分为0.1°×0.1°的块体，以每层块体的中心所在平面为计算基准。采用小波分析方法[13]，将重力场累积变化的趋势从图 2中扣除，用获得的局部重力场变化来进行反演计算。
2.2 密度变化反演结果由图 2的重力场动态变化，通过紧凑约束反演算法得到的相应密度变化结果如图 3(2014-04~2016-04)和图 4所示(2014-04~2017-04)。
|
图 3 2014-04~2016-04各层密度变化 Fig. 3 The density variation during 2014-04~2016-04 |
图 3(a)为研究区地壳浅部5 km附近的结果，大部分地区主要呈现很小的负密度变化，变化量在-0.03 g/cm3左右。图 3(b)为研究区地壳浅部15 km附近的结果，大部分区域密度变化不明显，在祁连山北侧的张掖附近出现正密度变化(最大变化量约0.06 g/cm3)，在西宁周边及其西部、延安周边呈现小范围的负密度变化(最大变化量约-0.06 g/cm3)。图 3(c)为部分中地壳25 km附近的密度变化，该部分密度变化较为显著，变化范围-0.1~0.1 g/cm3。在张掖地区附近(巴颜喀拉块体与柴达木地块交界区)出现正负密度变化梯度带，其中南部出现负密度变化，显示物质处在迁出的状态，而北部出现正密度变化，显示物质处于汇聚的状态。鄂尔多斯块体附近(银川、延安和西安一带)负变化(-0.1 g/cm3)向西北及西南方向逐渐减小延伸。在35 km深度的密度变化(图 3(d))更加明显，最大正密度变化达到0.14 g/cm3，最大负密度变化达到-0.15 g/cm3，与图 3(c)比较其变化量级较大，显示出深部物质运动更加活跃。在45 km深度的密度变化(图 3(e))最剧烈，在张掖附近已达到最大-0.2 g/cm3的负密度变化，并相对浅部更加向东南方向延伸，表明其物质迁移更为显著。在九寨沟地震震中的西北方向有一定量的物质积累，造成正密度变化(0.06~0.11 g/cm3)，形成一个包围巴颜喀拉块体东北部的趋势。在鄂尔多斯地块(银川、延安和西安一带)出现负的密度变化(-0.1~-0.15 g/cm3)，说明该位置也出现物质外迁。
由2014-04~2016-04深部地壳结果中看出，张掖地区两侧出现正负密度变化(图 3(a)~(c))，显示出差异性运动，由于受到鄂尔多斯块体以及华南块体的阻挡，在此持续挤压，其作用力在中下地壳达到最强，运动方向朝着九寨沟地震震中区域发展。
2.2.2 3 a尺度区域地壳密度变化特征图 4(a)为研究区浅部5 km深度的结果，有少量的负密度变化分布，但基本呈现不变的状态。图 4(b)为研究区地下15 km深度的密度变化，在震中两侧出现正负密度变化对称分布，在震中西北部出现负密度变化(-0.02~0.05 g/cm3)，在龙门山断裂带末端出现正密度变化(0.02~0.05 g/cm3)。图 4(c)为研究区地下25 km处的密度变化，变化范围在-0.08~0.1 g/cm3之间，在广元附近(龙门山断裂带西北部)的正密度变化加强(0.05~0.11 g/cm3)，并向震中方向延伸，在西宁北部出现正密度变化(0.03~0.05 g/cm3)。图 4(d)~(e)为35 km和45 km深度的密度变化，其变化范围和趋势相同，且十分剧烈，密度变化范围在-0.15~0.17 g/cm3之间。在九寨沟地震震中四周，出现类似四象限分布的图像，并呈现出类似左旋走滑的物质运移方向，这与震源机制解相符合。另一方面，对比图 4(d)和(e)，显示出深部物质运移更加活跃，此次地震的动力来源可能来自于深部的构造运动。
|
图 4 2014-04~2017-04各层密度变化 Fig. 4 The density variation during 2014-04~2017-04 |
从整个2014-04~2017-04的地壳内部密度变化来看，在九寨沟震中区域的NW-SE两侧出现明显的正负密度变化梯度带，并围绕震中区域对称分布。从纵向分布来看，两期密度变化趋势在深度上一致，且变化值随深度增加而增大，表明地壳深部活动更加活跃，来自深部的挤压作用力可能就是形成九寨沟地震的驱动力。
3 讨论受坚硬四川盆地的阻挡，巴颜喀拉块体物质向东流动改为向北、向南和向西南方向。其中，柴达木块体和祁连地块一起整体向东北方向挤压，逐渐转为东南向旋转挤压。该区域内部物质分布极其不均，物质流动产生的各个方向的作用力控制着整个地下构造环境的变化，具备发生强震的条件。
对比2014-04~2016-04和2014-04~2017-04两期密度变化结果，从空间分布看，整个研究区两期密度变化分布都不均匀，震前2014-04~2017-04结果更为显著。变化多分布于震中及其邻区，从比较散乱的正负变化状态(2014-04~2016-04)逐渐趋于相对震中九寨沟四周有序集中(2014-04~2017-04)，反映大震前区域构造活动增强和局部应力集中现象。从纵向看，反演得到的密度变化主要集中于地壳深部(35 km和45 km附近)，说明浅部物质运动较慢，主要是由块体整体运动引起，而深部物质运动活跃，主要是由密度变化引起。值得注意的是，2014-04~2016-04在巴颜喀拉块体东北部(最大为-0.16 g/cm3变化)的物质逐渐向东南方向移动，并在震前一段时间(2014-04~2017-04)使得广元附近(龙门山断裂带西北部)出现正密度变化加强(0.05~0.11 g/cm3)、震中西北方向减弱、最大为-0.08 g/cm3的负密度变化。这显示出青藏高原东南缘物质向东和东南方向移动；在震前2014-04~2017-04结果中，震中区及邻区形成密度变化的正负四象限分布(图 4中红色箭头)，巴颜喀拉块体东南部物质向四川盆地(最大为0.09 g/cm3变化)运移，鄂尔多斯块体西南部物质(最大为-0.10 g/cm3变化)向祁连山方向运移，整体的物质运移方向呈左旋走滑式分布，与震源机制解的结果基本一致。
本文基于改进的紧凑约束重力反演方法，在九寨沟震中及邻区获得震前两期(2014-04~2016-04、2014-04~2017-04)地壳密度变化特征。反演结果显示，本次地震震前一段时间里，地壳深部密度变化大，震中及其周边正负密度变化由分散变集中，并形成密度变化的正负四象限分布现象。反映出巴颜喀拉块体深部物质的向东和向东南移动使九寨沟附近物质积累，深部(35 km和45 km)构造运动剧烈，造成上层脆性地壳(本次地震震源深度在20 km左右)应力积累并达到临界状态。这些可能是造成九寨沟区域发生强震的主要原因。本文结果对于认识孕震区深部密度结构及地球物理场时空变化特征具有一定的意义，对识别与局部重力变化有关的地震前兆也有参考价值。
| [1] |
王椿镛, 杨文采, 吴建平, 等. 南北构造带岩石圈结构与地震的研究[J]. 地球物理学报, 2015, 58(11): 3867-3901 (Wang Chunyong, Yang Wencai, Wu Jianping, et al. Study on the Lithospheric Structure and Earthquakes in North-South Tectonic Belt[J]. Chinese Journal of Geophysics, 2015, 58(11): 3867-3901)
(  0) |
| [2] |
申重阳, 李辉, 孙少安, 等. 重力场动态变化与汶川MS8.0地震孕育过程[J]. 地球物理学报, 2009, 52(10): 2547-2557 (Shen Chongyang, Li Hui, Sun Shaoan, et al. Dynamic Variations of Gravity and the Preparation Process of the Wenchuan MS8.0 Earthquake[J]. Chinese Journal of Geophysics, 2009, 52(10): 2547-2557 DOI:10.3969/j.issn.0001-5733.2009.10.013)
( 0) |
| [3] |
祝意青, 梁伟锋, 赵云峰, 等. 2017年四川九寨沟MS7.0地震前区域重力场变化[J]. 地球物理学报, 2017, 60(10): 4124-4131 (Zhu Yiqing, Liang Weifeng, Zhao Yunfeng, et al. Gravity Changes before the Jiuzhaigou, Sichuan, MS7.0 Earthquake of 2017[J]. Chinese Journal of Geophysics, 2017, 60(10): 4124-4131 DOI:10.6038/cjg20171037)
( 0) |
| [4] |
申重阳. 地壳形变与密度变化耦合运动探析[J]. 大地测量与地球动力学, 2005, 25(3): 8-12 (Shen Chongyang. Preliminary Analysis of Coupling Movement between Crustal Deformation and Density Change[J]. Journal of Geodesy and Geodyamics, 2005, 25(3): 8-12)
( 0) |
| [5] |
Xuan S B, Shen C Y, Li H, et al. Structural Interpretation of the Chuan-Dian Block and Surrounding Regions Using Discrete Wavelet Transform[J]. International Journal of Earth Sciences, 2016, 105(5): 1591-1602 DOI:10.1007/s00531-015-1272-1
( 0) |
| [6] |
陈石, 王青华, 王谦身, 等. 云南鲁甸MS6.5地震震源区和周边三维密度结构及重力场变化[J]. 地球物理学报, 2014, 57(9): 3080-3090 (Chen Shi, Wang Qinghua, Wang Qianshen, et al. The 3D Density Structure and Gravity Change of Ludian MS6.5 Yunnan Epicenter and Surrounding Regions[J]. Chinese Journal of Geophysics, 2014, 57(9): 3080-3090 DOI:10.6038/cjg20140933)
( 0) |
| [7] |
胡敏章, 李辉, 刘子维, 等. 川滇地区2010~2013年重力变化及重力网的地震监测能力[J]. 大地测量与地球动力学, 2015, 35(4): 616-620 (Hu Minzhang, Li Hui, Liu Ziwei, et al. The Gravity Change over Sichuan-Yunnan Region in 2010-2013 and the Earthquake Monitoring Ability of the Gravimetric Network[J]. Journal of Geodesy and Geodynamic, 2015, 35(4): 616-620)
( 0) |
| [8] |
邢乐林, 李辉, 李建国, 等. 陆态网络绝对重力基准的建立及应用[J]. 测绘学报, 2016, 45(5): 538-543 (Xing Lelin, Li Hui, Li Jianguo, et al. Establishment of Absolute Gravity Datum in CMONOC and Its Application[J]. Acta Geodaetica et Cargographica Sinica, 2016, 45(5): 538-543 DOI:10.11947/j.AGCS.2016.20140653)
( 0) |
| [9] |
Last B J, Kubik K. Compact Gravity Inversion[J]. Geophysics, 1983, 48(6): 713-721 DOI:10.1190/1.1441501
( 0) |
| [10] |
Guillen A, Menichetti V. Gravity and Magnetic Inversion with Minimization of a Specific Functional[J]. Geophysics, 1984, 49(8): 1354-1360 DOI:10.1190/1.1441761
( 0) |
| [11] |
Li Y, Oldenburg D W. 3-D Inversion of Gravity Data[J]. Geophysics, 1998, 63(1): 109-119 DOI:10.1190/1.1444302
( 0) |
| [12] |
孙少安, 项爱民, 朱平, 等. 三峡水库首次蓄水引起的重力变化及其机制的初步研究[J]. 地震学报, 2006, 28(5): 485-462 (Sun Shaoan, Xiang Aimin, Zhu Ping, et al. Gravity Change and Its Mechanism after the First Water Impoundment in Three Gorges Project[J]. Acta Seismologica Sinica, 2006, 28(5): 485-462)
( 0) |
| [13] |
杨文采, 施志群, 侯遵泽, 等. 离散小波变换与重力异常多重分解[J]. 地球物理学报, 2001, 44(4): 534-541 (Yang Wencai, Shi Zhiqun, Hou Zunze, et al. Discrete Wavelet Transform for Multiple Decomposition of Gravity Anomalies[J]. Chinese Journal of Geophysics, 2001, 44(4): 534-541)
( 0) |
2. Laurel Geophysical Instruments Co Ltd(Shanghai), 777 Longwu Road, Shanghai 200232, China