龙门山断裂周边区域在汶川Ms8.0地震和芦山Ms7.0级地震前的地壳形变特征对比研究*

周德敏 甘卫军 李金平 陈为涛 丁晓光 梁诗明

1)中国地震局地质研究所，北京 100029

2)云南师范大学旅游与地理科学学院，昆明 650500

3)中国地震局地壳运动监测工程研究中心，北京100036

2013-04-20 四川芦山发生Ms7.0级强震，震中位于30.314°N、102.934°E，震源深度12.3 ～15 km，震源断层走向N40°E，倾角约35°［1－2］。推测发震断层为龙门山断裂系前山断裂的双石-大川断裂，但也不排除是大邑断裂的可能性［3］。无论是双石-大川断裂还是大邑断裂，其破裂发震的构造动力学背景都是青藏高原中东部巴颜喀拉块体的东向逃逸对华南块体西北端四川盆地的强烈挤压，属碰撞带内部的逆断层型地震。本次地震与2008年汶川Ms8.0地震相距不到100 km，两个地震是否为主-余震关系尚未取得共识。

汶川Ms8.0地震后，许多学者从地震应力触发的角度研究了龙门山断裂带周边区域的地震危险性，认为龙门山断裂的西南段处于应力加载状况［4－5］。因而，本次地震的发生可认为是该区域地壳应力场和龙门山断裂带应力分布状况在汶川地震5 a 后的又一次调整，而这种应力调整，必然伴随着一定区域地壳形变或地壳应变率的变化。本文基于龙门山地区多年来由“中国地壳运动观测网络”、“中国大陆构造环境监测网络”GPS 观测站和汶川Ms8.0地震后局域连续GPS 观测网所积累的GPS地壳运动观测资料，分析对比龙门山断裂周边地区在两次大震之前的地壳形变特征，以便更好地认识大震孕育区的震前地壳形变场变化规律和方式。

1 GPS 资料处理

本文选取龙门山断裂周边74 个GPS 观测站(图1)，站点类型和观测状况见表1。采用JPL 开发的高精度GPS 处理软件GIPSY-OASIS 处理上述资料。

表1 本文所用龙门山断裂周边GPS 站点观测信息Tab.1 Information of GPS stations aroundin Longmenshan fault zone

图1 龙门山断裂带及周边地区GPS 站点分布Fig.1 Distribution of GPS stations around Longmenshan fault zone

第一，首先从JPL 网站获取相应观测日的IGS 精密卫星轨道、精密卫星时钟改正、卫星掩蔽信息等信息，采用精密单点定位方法［6］，以消电离层影响的线性组合作为观测量，卫星截止高度角15°，数据采样率300 s。固定JPL 精密星历和时钟，采用先验对流层干延迟模型和对流层投影函数分别为全球大气压和温度模型GPT(global pressure and temperature)及GMF［7］进行处理;海洋潮汐影响改正由FES2004 海洋潮汐模型基于格林函数在线计算获得(http://holt.oso.chalmers.se/loading)。处理过程中考虑星历、钟差、仪器偏差(DCB)、接收机天线相位中心和卫星相位中心等产品的自洽性。估计参数除基准站坐标和接收机钟差外，还包括对流层延迟参数，从而解算每个本地站和部分IGS 全球核心参考站的单日坐标松弛约束解。其中，部分IGS 全球核心站的加入是为了建立本地网与全球参考框架之间的联络纽带。

第二，将已获得精密单点定位结果的所有站点组网，进一步采用载波相位整周模糊度解的Ambizap 算法［8］，利用固定点法则确定观测网参数的各种线性组合，生成唯一、自洽且消除整周模糊度的单日解，并通过7 参数变换将单日解转换到ITRF2008框架下［9］。对网中所有独立的站间基线进行整周模糊度影响改正，进一步提高单日坐标松弛约束解的精度。

第三，根据JPL 提供的每日参考框架转换参数，将单日松弛约束解转换到指定的全球参考框架下(如ITRF2008)，得到各站点的坐标变化时间序列。

第四，利用QOCA 软件［10］对上述龙门山地区所有GPS 观测站和全球80 个IGS 核心参考站不同时期、不同网络的单日松弛约束解进行联合平差。在假定各站点三维坐标随时间呈线性变化的约束下，确定ITRF2008 参考框架下各站点坐标和水平速度的估值及其误差。

2 两个地震前的GPS 速度场差异

选取龙门山断裂东侧四川盆地内部的9 个相对较为稳定的GPS 站点，建立一个“稳定四川盆地”的参考框架，将两套GPS 速度场均归化到这一参考框架下。首先通过这9 个点在ITRF2008 参考框架下的速度矢量，求解欧拉旋转参数;然后，将区域内所有的GPS 速度矢量扣除这一欧拉刚性旋转，从而使这9 个点相对于稳定的四川盆地不再有整体的平移或旋转。这样，其他点旋转后的速度矢量，即被归化至“稳定四川盆地”参考框架。图2 和图3 分别为两个地震之前的GPS 速度场。

由图2 可见，在汶川地震前，位于整个龙门山断裂带西侧的H051、H050、H045、H061 等GPS 观测点相对于稳定的四川盆地，水平向缩短量小于2 mm/a，即相对于四川盆地并没有明显的水平向缩短变形。此外，介于玛尼-玉树-鲜水河断裂带与东昆仑断裂带之间的巴颜喀拉块体东向挤出受到华南块体的阻挡，在高塑性上地壳物质的挤压和重力共同作用下，龙门山断裂带中段、北段向东北方向逃逸，量值为1 ～2 mm/a。

在汶川地震后(见图3)，由于震后形变弛豫效应，位于龙门山断裂带中北段的H045、H050、H051等GPS 观测点相对于稳定的四川盆地，水平向缩短量不小于10 mm/a;而位于龙门山断裂带南段的H061、H377、H384 等GPS 观测点与汶川地震前相同，其相对于稳定的四川盆地，水平向缩短量小于2 mm/a，并没有明显的水平向挤压缩短。

图2 汶川Ms8.0地震前GPS 水平运动速度场(“稳定四川盆地”参考框架)Fig.2 Horizontal velocity field before Wenchuan Ms8.0 earthquake(In Sichuan basin fixed reference frame)

图3 芦山Ms7.0地震前GPS 水平运动速度场(“稳定四川盆地”参考框架)Fig.3 Horizontal velocity field before Lushan Ms7.0 earthquake(In Sichuan basin fixed reference frame)

需要说明的是，对比汶川Ms8.0地震震前的速度场，芦山Ms7.0地震前GPS 站点的速度场误差相对较大，主要是由于GPS 观测时段跨度较短，但这种相对较大的误差不会对基本认识有实质性的影响。

3 震前应变率场

利用位移求解应变，分为平面坐标系求解和球面坐标系求解两大类。通常情况下，如果位移场覆盖范围较小，则可用忽略地球曲率影响的平面坐标系公式计算相应区域的均匀应变、整体平移和旋转［11］。但更严密的做法应是无论区域大小，均尽可能采用球面坐标系下的计算方法，以克服地球曲率或投影变形等因素产生的影响［12］。采用下式进行应变计算［13］:

式中，φ、θ 和r 为地球坐标系下的经度、余纬度和地心距;uφ、uθ和ur表示各观测点的东向、北向和垂向位移;φ0、θ0和r0表示观测网中心点的坐标;Uφ、Uθ和Ur表示观测网中心点的三方向位移;Δφ、θ0和r0表示观测点到观测网中心点的距离;和表示各应变分量;和表示旋转量。只要知道3 个或以上观测点的东向、北向和垂向三维速度值，就可用最小二乘拟合，求出观测网的9 个参数Uφ、和

为了由离散的GPS 速度矢量获得空间上连续变化的应变率场，首先采用二维“高张力样条”(τ=0.95)内插算法对非均匀的GPS 速度场进行0.5° ×0.5°的均匀预测加密，然后根据每个1° ×1°网格边界及内部的9 个内插速度矢量，计算出相应1° ×1°网格内的均匀应变率。这种基于9 个加密网格点的单元区域均匀应变率计算，既能控制应变计算单元网格的密度，又能充分反映那些高密度观测区实际观测结果对应变计算的贡献。

这里所采用的二维张力样条是一种基于所有已知点的通用内插算法。其中的张力参数τ(0≤τ ＜1)表示曲面膜中张力应变能所占总应变能的份额。当τ=0 时，这种内插算法等同于曲率平方和最小的双调和样条。随着张力的增加，整个曲面膜曲率的平方和也随之增加。而当τ→1 时，这种算法将等同于基于所有已知点的线性内插。Gan［14］通过对数种方法的比较发现，“高张力样条”法对速度场的内插效果较优。基于上述方法，获取龙门山地区在汶川Ms8.0地震和芦山Ms7.0地震前的地壳应变率场(图4、图5)。

在2008年汶川Ms8.0地震前(图4)，龙门山断裂带中段、北段的典型应变率值小于10 nanostrain/a，而南段的应变率值约为30 nanostrain/a，反映了巴颜喀拉地块东部物质向东顺时针扭转运动并未引起龙门山断裂带中段、北段的压缩应变。

图4 汶川Ms8.0地震震前应变率场Fig.4 Strain rate filed in Longmenshan fault zone before Wenchuan Ms8.0 earthquake

图5 芦山Ms7.0地震前地壳应变率场Fig.5 Strain rate filed in Longmenshan fault zone before Lushan Ms7.0 earthquake

在2008年汶川Ms8.0地震后(图5)，由于地震弛豫效应对周围区域的影响，龙门山断裂带周边地区的应变率值明显增加，其中龙门山断裂带中段、北段的水平应变率值达到200 nanostrain/a 以上，四川盆地、巴颜喀拉地块及鲜水河断裂带均有明显增大;但芦山Ms7.0地震的发震断层附近，应变率值在汶川地震前、后未有明显变化，挤压缩短的典型应变率仍为30 nanostrain/a，其结果与Wu 等［15］获得的龙门山断裂带南段水平应变率明显增大的结论不同。

4 结 论

本研究初步结论如下:

1)受汶川Ms8.0地震震后形变弛豫效应影响，龙门山断裂带中北段水平向缩短量不小于10 mm/a，而龙门山断裂带南段与汶川地震震前相比，其水平向缩短量小于2 mm/a，并没有明显的水平向挤压缩短。

2)由于汶川Ms8.0地震对周边区域震后地壳变形的调整作用，龙门山断裂带南段的水平应变率在震后有明显增大，但对比周边区域，芦山Ms7.0地震震中附近区域由于孕震断层及其附近区域的地壳形变因弹性变形已趋于极限状态，其地壳形变趋缓或保持低的形变状态而应变率值为一个相对的“低值区”，其挤压缩短的典型应变率为30 nanostrain/a，而周边挤压缩短的典型应变率量值至少在60 nanostrain/a 以上。

3)两次强震都发生在挤压缩短的典型应变率“低值区”断层上，其原因可能是即将发生大地震的孕震区域，其应力已处于岩石破裂的临界状态，岩石的弹性变形积累也已趋于极限状态，因此，其地壳形变的增量或应变率将趋缓或保持低值。

1 陈运泰，杨智娴，张勇.浅谈庐山地震［J］.地震学报，2013，35(3):285－295.(Chen Yuntai，Yang Zhixian，Zhang Yong.A brief talk on the 20 April 2013 Lushan Mw6.7 earthquake［J］.Acta Seismologica Sinica，2008，35(3):285－295)

2 张勇，许立生，陈运泰.庐山4.20地震破裂过程及其致灾特征初步分析［J］.地球物理学报，2013，56(4):1 408－1 411.(Zhang Yong，Xu Lisheng，Chen Yuntai.Rupture process of the Lushan 4·20 earthquake and preliminary analysis on the disaster-causing mechanism［J］.Chinese J Geophys，2013，56(4):1 408－1 411)

3 徐锡伟，闻学泽，韩竹军，等.四川芦山7.0级强震:一次典型的盲逆断层型地震［J］.科学通报，2013(20):1－7.(Xu X W，Wen X Z，Han Z J，et al.Lushan Ms7.0 earthquake:A blind reserver-fault earthquake［J］.Chin Sci Bull，2013(20):1－7)

4 Parsons T，Ji C.Stress changes from the 2008 Wenchuan earthquake and increased hazard in the Sichuan basin［J］.Nature，2008，454(7 203):509－510.

5 万永革，沈正康，盛书中，等.2008年汶川大地震对周围断层的影响［J］.地震学报，2009，31(2):128－139.(Wan Yongge，Shen Zhengkang，Sheng Shuzhong，et al.The influence of 2008 Wenchuan earthquake on surrounding faults［J］.Acta Seismologica Sinica，2009，31(2):128－139)

6 Zumberge J F，Helfin M B，Jefferson D C，et a1.Precise point positioning for the efficient and robust analysis of GPS data from large networks［J］.J Geophys Res，1997，102(B3):5 005－5 017.

7 Boehm J，Werl B，Schuh H.Troposphere mapping functions for GPS and very long baseline interferometry from European center for midium-range weather forcasts operational analysis data［J］.Journal of Geophysical Research，2006，111:B02406.

8 Blewitt G.Fixed point theorems of GPS carrier phase ambiguity resolution and their application to massive network processing:ambizap［J］.Journal of Geophysical Research，2008，113:B12410.

9 Altamimi Z，Collilieux X，Metivier L.ITRF2008:an improved solution of the international terrestrial reference frame［J］.Journal of Geodesy，2011，85:457－473.

10 Dong D N，Fang P，Bock Y，et al.Anatomy of apparent seasonal variations from GPS derived site position time series［J］.Journal of Geophysical Research，2002，107:901－916.

11 Jaeger J C.Elasticity，fracture and flow［M］.New York:Methuen，1964.

12 石耀霖，朱守彪.用GPS 位移资料计算应变方法的讨论［J］.大地测量与地球动力学，2006，26(1):1－8.(Shi Yaolin，Zhu Shoubiao.Discussion on method of calculating strain with GPS displacement data［J］.Journal of Geodesy and Geodynamics，2006，26(1):1－8)

13 Savage J C，Gan Weijun，Svarc J L.Strain accumulation and rotation in eastern California shear zone［J］.Journal Geophysical Research，2001，106(B10):21 995－22 007.

14 Gan W，Zhang P.Present-day crustal motion within the Tibetan plateau inferred from GPS measurements［J］.Journal of Geophysical Research，2007，112(B8):B08416.

15 Wu Y Q，Jiang Z S，Wang M，et al.Preliminary results pertaining to co-seismic displacement and pre-seismic strain accumulation of the Lushan Ms7.0 earthquake，as reflected by GPS surveying［J］.Chin Sci Bull，2013，58(20):3 460－3 446.