地壳内部的应力一般由静岩应力和构造应力两部分组成.静岩应力(也称正压力)主要由重力决定，其大小是深度的函数.地壳内部压力与上覆物重力之间的平衡状态称为静岩应力状态，与静止流体中的应力状态完全等价.静止流体中任何一点的压力在各个方向相等，并随着深度线性增加.然而，在地球内部，尤其在地球浅层，岩石并非全部处于静岩应力状态，而多为多种应力状态，当三个主应力不等时，即存在构造应力.其实，构造应力就是附加在正压力之上的一种应力，又称偏应力，是由相邻地块地形高低、密度、力的大小和方向、温度场等的不同所产生的负载应力差(党亚民等，2009).

与构造应力状态相关的研究是一个重要的地学研究方向，不少学者对重力异常、密度异常与构造应力的关系进行过探讨，或者对典型区域的构造应力场分布进行过研究.Turcotte和Schubert(1982)推导了使用地壳平均密度计算构造应力的计算公式；Bowin等(1986)利用点质量元法建立起场源深度与卫星重力位系数阶数之间的关系；Hardebeck和Michael(2004)探讨了圣安德斯断裂带上的应力方向.我国的构造应力研究始于20世纪70年代，之后很多学者在此方面取得了一系列成果.张赤军等(1998)给出了利用地面、卫星重力资料研究岩石层密度的方法；方剑(1994)对这种方法进行了推广，给出了用重力场模型计算不同深度密度异常的公式；张永庆等(2009)利用1996年丽江地震序列反演了震区应力状态；邱泽华等(2010)研究了汶川地震前姑咱地震台观测到的异常应变变化.万永革等(2011)探讨了不同应力状态和摩擦系数对综合P波辐射花样的影响.总体来说，目前关于构造应力状态的研究途径主要包括地质资料分析、原地应力测量、震源机制解、断层擦痕反演、实验模拟以及数值模拟等.

地震是岩石圈应力场动态演化的结果，一般认为是由地球内部的岩体受到构造应力作用，导致岩体突然断裂错动而产生，因此研究地壳构造应力场分布具有重要的科学意义.鉴于此，本文提出一种基于流动重力/GPS联合观测数据计算构造应力的新方法，该算法可望进一步挖掘目前越来越丰富的流动重力/GPS联合观测数据的应用潜力，并为地震中长期预测研究提供支持.

重力/GPS结合是一种新的测量方法(宁津生和王正涛，2010).近年来，基于重力/GPS数据的研究取得了不少进展.王谦身等(2013)对陕西榆林—重庆综合地球物理大断面中从陕西咸阳至中秦岭北侧测段特异重力场进行了分析和探讨.陈石等(2014)研究了云南鲁甸M_S6.5地震震源区和周边三维密度结构及重力场变化.Zhang等(2014)利用两条垂直于龙门山断裂带的布格重力异常剖面数据，通过莫霍面和均衡面的对比发现龙门山以东地区地壳基本处于均衡状态，而龙门山以西地区则处于非均衡状态.Fu等(2014)发现龙泉山脉是四川盆地地壳均衡状态的分界线，并据此推测龙泉山断裂带是一条深埋的高角度断裂带.Fu和Zhang(2014)发现龙门山中南段地壳均衡负异常显著，具有发生7级以上强震的深部动力条件.杨光亮等(2015)利用重力/GPS观测数据反演了金川—芦山—犍为观测剖面的密度构造，分析了芦山地震及龙门山地区地壳构造背景和深部动力环境特征.

综上可知，利用重力/GPS联测数据进行研究是国内地学界近年来的一个研究热点，并取得了不错的成绩，但迄今为止尚未看到利用重力/GPS联测数据进行地壳构造应力计算方面的研究.

地球是一个具有圈层结构的球体，大体可分为地壳、地幔和地核三层.从漫长的地质年代来看，相对轻且坚硬的地壳漂浮在重而柔软的地幔之上，如同固态冰漂浮在液态水之上.如果没有外力作用，即处在一种理想的静水平衡状态，地壳自然地漂浮在地幔之上(如图 1a)，此时地壳处于均衡状态，处于地壳底部的莫霍面与均衡面重合.换句话说，当不存 在垂向外力作用时，莫霍面与均衡面一致.但是，当 有垂向外力作用于地壳时，均衡面将偏离莫霍面，地壳就不均衡(图 1b).加在地壳上的垂向外力越大，莫霍面与均衡面之间的差异也越大，地壳越不均衡(图 1c).

图 1

Fig. 1

图 1 莫霍面、均衡面的分布与地壳承载外力之间的关系 Fig. 1 The relationship between the force the crust bears and the distribution of the Moho and the isostatic crustal thickness

对于实际地球而言，地壳具有弹性，可以承载各个方向的构造应力，包括垂向构造应力.局部地区均衡面与莫霍面之间的剩余物质(地壳地幔物质密度差)所产生的附加浮力，一般是通过地壳来承载并平衡.理论上讲，均衡面与莫霍面之间的剩余物质所产生的附加浮力，基本上与地壳所承载的垂向外力大小相等、方向相反.因此，我们可以通过莫霍面与均衡面之间的剩余物质承受的附加浮力来定量计算地壳所承受的垂向构造应力.根据Fu等(2014)可知，基于地表进行的流动重力/GPS联合观测数据，可推算莫霍面与均衡面之间的差异.若考虑其间剩余物质密度，即可进一步计算区域地壳所承受的垂向构造应力.

依据流动重力/GPS联测数据计算地壳承载的垂向构造应力的具体步骤如下：(1)依据流动重力/GPS联合观测数据，经过固体潮改正、中间层改正、高程改正和地形改正，获得测区布格重力异常场；(2)参考区域岩石圈密度构造模型，依据布格重力异常数据推算莫霍面深度；(3)依据GPS观测获得的高程数据，利用均衡理论计算均衡面深度；(4)依据均衡面与莫霍面之间的差异，以及其间的剩余物质密度，计算剩余物质产生的附加浮力，该附加浮力与地壳所承受的垂向构造应力大小相等、方向相反.

当地幔和地壳的密度分别简化为ρ_M和ρ_C时(图 2a)，根据浮力原理，地壳承载的垂向构造应力的计算公式如下：

图 2

Fig. 2

图 2 Moho面与均衡面之间剩余物质产生的附加浮力计算模型，该附加浮力与地壳承载的垂向构造应力大小 相等、方向相反 (a) 均一密度地壳模型； (b) 层状密度地壳模型. Fig. 2 Computation model of additional buoyancy generated by the remaining materials between the Moho and the isostatic crustal thickness. This additional buoyancy is equal to the vertical tectonic stress of the Crust in magnitude, but in opposite direction (a) A uniform crustal model; (b) A layered crustal model.

(1)

公式(1)中，d_Moho为Moho面深度，d_iso为均衡面深度，g为地表重力加速度.当地壳密度为分层构造时(图 2b)，假设均衡面处在密度层ρ_m的内部，则相应的垂向构造应力计算公式为

(2)

公式(2)中，ρ_i和d_i(i=1，2，…，n)分别代表从地表到地壳深部各密度层的密度和厚度，ρ_m为第m层地壳的密度，ρ_M为地幔的密度.

需要特别指出的是，本文计算的地壳垂向构造应力，不是指构造应力随深度的详细分布，而是指地壳整体承受的构造应力，即垂向构造应力沿着地壳垂向积分的结果.此外，莫霍面不是固态与液态物质的分界线，而是一个化学界面，其两侧物质的物性明显不同，但实际固液分界线的位置不影响上述附加浮力以及垂向构造应力的计算结果.

2008年汶川M_W7.9地震以来，不少学者在巴颜喀拉块体东边界及周边地区开展了高精度流动重力/GPS联合观测(图 3).中国科学院地质与地球物理研究所跨龙门山断裂带中北段建立了两条流动重力/GPS联合观测剖面(Zhang et al.，2014)；湖北省地震局与中国地震局地球物理研究所共同在芦山M_W6.6地震周围地区展开流动重力/GPS联合剖面观测(杨光亮等，2015).此外，中国地震局地震预测研究所于2012年围绕龙泉山脉建立了一个重力观测网，并进行了流动重力/GPS联合观测(Fu et al.，2014).上述观测相互补充，形成了一个较为理想的重力/GPS联合观测网(图 3)，为巴颜喀拉块体东边界垂向构造应力的计算提供了坚实的数据基础.图 3中，黑色粗线条为流动重力/GPS联合观测路线，红色线条为区域主要断层.汶川—茂汶断裂带、映秀—北川断裂带和灌县—安县断裂带组合成著名的龙门山断裂带，该断裂带是巴颜喀拉地块和华南地块的边界带，两侧的地形地貌差异巨大.巴颜喀拉地块上地壳朝南东的水平运动在四川盆地西缘受到华南地块的阻挡，转换成龙门山断裂带的逆冲运动(杜方等，2009).本文统一整理了上述三个来源的布格重力异常数据，获得巴颜喀拉块体东边界及周边地区布格重力异常场.

图 3

Fig. 3

图 3 巴颜喀拉块体东边界流动重力/GPS联合观测网 Fig. 3 The hybrid gravity and GPS observation network at the east boundary of Bayanhar block

Moho面是地球内部接近地表的一个比较大的间断面，是地壳与地幔之间的一个化学界面，同时也是一个地震波速度间断面.Moho面的变形通常可用来衡量上地幔岩石圈的变形.由于缩短效应，青藏高原的地壳厚度大约是地球大陆地壳厚度平均值的两倍.从四川盆地到青藏高原内部，Moho面的巨大变化是观测到的布格重力异常的最主要贡献(Jiang et al.，2004).布格重力异常反应的是地球内部物质密度的横向变化，应用布格重力异常数据，通过重力正反演分析，可得到相应的深部地壳构造，特别是莫霍界面(Moho)深度(地壳厚度)的起伏变化(王谦身等，2009).

为简化计算，我们假设研究区域的地壳密度均一(模型2a)，则该地区地壳厚度的变化可近似表示为密度层的缓慢变化，该密度层的密度为地幔与地壳的密度差.假设研究区域地壳密度为2.67 g·cm^-3(海平面以上物质的平均密度)，则地表观测的布格重力异常变化可简化为由地下Moho面附近的密度层引起，因此可以直接利用布格重力异常数据计算Moho面的起伏变化.Sun和Zhou(2012)把大地震引起的海水再分布导致的重力变化转化为密度层的变化，其基本研究思路与本文一致.具体地说，若地表重力观测站下方H米深度处一个较大的局部区域内存在一个厚度为h的密度层，该密度层的密度为Δρ，则该密度层在地表产生的重力变化可表示为(郭俊义，2001)

(3)

公式(3)表示的是观测点下方局部区域内厚度为h而且可被近似当作平面的质量层在地表引起的重力异常值的近似值.其中，G为万有引力常数.Wang等(2010)对四川地区进行过详细研究，认为该地区 地壳地幔的密度差异Δρ为0.649±0.036 g·cm^-3. 取Δρ=0.649 g·cm^-3，代入公式(3)，可获得Δg与h之间的定量关系.进一步，如果我们已知一个测站的地壳厚度h₀与布格重力异常值Δg₀，则整个重力观测网内任意测站对应的地壳厚度可根据该测站的布格重力异常值推测出来(Fu et al.，2014).与我们研究小组的前期工作(Fu et al.，2014)类似，我们以四川盆地的资阳站为基础进行分析.资阳站所在地区的均衡地壳厚度为39.5287 km，观测到的布格重力异常为-113.9032 mGal.假设该处地壳厚度处于均衡状态(该测站距离青藏高原较远，处于均衡状态的可能性较大)，则图 3所示的重力网所有测站对应的Moho面深度都可以推测出来.

图 4给出了巴颜喀拉块体东边界及周边地区Moho面深度分布.测线附近的Moho面深度计算结果精度较高，用彩色标示.距离测线较远地区的Moho面结果是基于插值方法获得，精度相对较低，用灰色等值线标示.整体上，研究区域Moho面由东南向西北逐渐增厚，变化趋势比较平缓.在华南块体的四川盆地区域的Moho面深度在40 km左右变动，但进入巴颜喀拉块体的松潘高原以后，Moho面深度逐步增加到50 km以上.

图 4

Fig. 4

图 4 巴颜喀拉块体东边界及周边地区Moho面深度分布 Fig. 4 The distribution of Moho depth at the east boundary of Bayanhar block and around area

由GPS观测数据，依据均衡理论可以确定不同高程对应的均衡补偿厚度，然后根据区域平均地壳厚度可进一步计算该处的均衡地壳厚度(王谦身等，2009).具体地说，根据Airy均衡理论，处于均衡状态的地壳厚度H与地形高度h满足如下关系(郭俊义，2001)：

(4)

其中，r是Airy均衡模型的山根，H₀是海平面以下平均地壳厚度，ρ_t是海平面以上物质的密度(2.67 g·cm^-3)，Δρ是地壳和地幔物质的密度差异.依据Wang等(2010)，研究区域壳幔密度差Δρ与地壳平均厚度分别取为0.649 g·cm^-3和37.9 km.依据图 3给出的GPS高程数据，本文利用公式(4)估算了巴颜喀 拉块体东边界及周边地区均衡地壳厚度分布(图 5).

图 5

Fig. 5

图 5 巴颜喀拉块体东边界及周边地区均衡面深度分布 Fig. 5 The distribution of isostasy thickness at the east boundary of Bayanhar block and around area

图 5给出了巴颜喀拉块体东边界及周边地区地壳均衡厚度分布形态.同图 4，精度较高的结果用彩色标示，精度相对较低的结果用灰色等值线标示.由图 5可知，研究区域内部的均衡地壳厚度基本上在38~56 km之间变化，且受地形影响较大，总体分布形态不是特别规则，其中山区短波长信号较多，呈现出明显的不规则特性.

将由布格重力异常数据计算出的Moho面深度与由Airy均衡理论计算得到的均衡地壳厚度做对比，给出二者的差异.当某处二者的差异很小或者接近于零时，即表明该地处于均衡状态.如差异较大，则处于不均衡状态.二者相差越大，表明该区地壳越不均衡.当均衡面深度大于Moho面深度时，根据大陆均衡原理，该地区应该“下降”，使地形高度减小，才能趋于均衡；反之，当均衡面深度小于Moho面深度时，地壳亦不均衡，该区应该“上升”(王谦身等，2009).总之，通过均衡面与Moho面深度的比较，可确定某一地区(地带)是否达到均衡，且可借此计算地壳承载的垂向构造应力的大小和方向.

均衡面与Moho面的不一致一般会造成附加浮力，该浮力主要由地壳来承载.因此，根据均衡面与莫霍面之间的深度差，以及地壳与地幔之间的密度差异(即剩余密度)，可计算区域地壳承载的垂向构造应力.取地壳地幔密度差为Δρ=649 kg·m^-3，则巴颜喀拉块体东边界及周边地区垂向构造应力分布如图 6所示.同图 4，精度较高的结果用彩色标示，精度相对较低的结果用灰色等值线标示.

图 6

Fig. 6

图 6 巴颜喀拉块体东边界及周边地区地壳垂向构造应力分布 Fig. 6 The distribution of vertical tectonic stress of the Crust at the east boundary of Bayanhar block and around area

由图 6可知，龙泉山断裂带以东地区，地壳承载的垂向构造应力几乎为零.龙泉山断裂带与龙门山断裂带之间地区，地壳承载的垂向构造应力为正值(向上为正).向西穿越龙门山断裂带，地壳承载的垂向构造应力再度接近为零.到隶属于巴颜喀拉地块的松潘高原等地，地壳承载的垂向构造应力转为负值.总体上，研究区域从东南到西北形成“零-正-零-负”模式的地壳垂向构造应力总体分布形态.研究区域西北部的松潘高原地区，地壳承载的垂向构造应力大约为10~20 MPa，该地区的垂向构造应力释放的较为完整，近期再次发生大地震的可能性不大.研究区域西南部的鲜水河断裂带东南段附近地区，地壳承载的垂向构造应力较大，大约为40~50 MPa，该地区垂向构造应力积累较高，且梯度变化剧烈.鲜水河断裂带东南段附近明显的高梯度垂向构造应力分布，与该地区较大的地形落差有关.

巴颜喀拉块体东边界及周边地区的上述垂向构造应力分布状态可用龙门山挤压增厚模型解释.Wang等(2011)指出，龙门山是由于四川盆地地壳被动插入到龙门山下部，引起地壳增厚造成.在这种龙门山挤压增厚模式下，四川盆地西部地壳由于龙门山的压力作用，引起下沉，致使Moho面深度大于均衡面深度，该效应对应向上的垂向构造应力；反之，松潘高原地区在其底部被动插入的四川盆地地壳向上的弯曲应力作用下，Moho面深度小于均衡面深度.为保持平衡，松潘高原地区的地壳产生向下的垂向构造应力.

值得注意的是，汶川M_W7.9地震和芦山M_W6.6 地震都不是发生在垂向构造应力很高的地区，而是发生在由正向负变化的高梯度带附近.祝意青等(2013)发现，强震易发生在重力变化正、负异常区过 渡的高梯度带上，在汶川M_W7.9地震和芦山M_W6.6 地震前均有重力场正负异常高梯度变化.垂向构造应力场可能和重力场时空变化一样，对强震地点预测具有指示意义，相关研究有待深入进行.

本文基于流动重力/GPS联合观测数据，提出了一种计算地壳垂向构造应力的新方法，拓展了流动重力/GPS联合观测数据的应用潜力.该方法的计算步骤如下：首先利用重力/GPS联合观测数据，经布格改正获取区域布格重力异常场，然后结合区域地震波层析成像结果或密度模型，反演Moho面分布；接着依据GPS高程数据，利用地壳均衡模型计算区域均衡面分布；最后，利用均衡面与Moho面之间剩余物质所受附加浮力与地壳承载垂向应力大小相等、方向相反的基本思想，依据均衡面与Moho面之间的差异及其间的剩余密度，计算地壳承载的垂向构造应力.

本文收集了三个研究小组(Zhang et al.，2014; Fu et al.，2014; 杨光亮等，2015)给出的重力/GPS联合观测数据，利用上述构造应力计算方法，计算了巴颜喀拉块体东边界及周边地区的地壳垂向构造应力分布形态，发现龙泉山断裂带以东地区的地壳垂向构造应力基本为零，龙泉山断裂带与龙门山断裂带之间地区的地壳垂向构造应力为正值，然后向西逐渐转化为负值.鲜水河断裂带东南段周边蓄积了40~50 MPa的垂向构造应力，且应力梯度变化剧烈；松潘高原垂向应力蓄积相对较小，大约只积蓄了10~20 MPa的垂向构造应力.汶川M_W7.9地震和芦山M_W6.6地震都发生在垂向构造应力正、负过渡区的高梯度带上，发生在垂向构造应力零值线附近.最后，利用龙门山挤压增厚模型解释了巴颜喀拉块体东边界从东南到西北的“零-正-零-负”模式的地壳垂向构造应力总体分布形态.

感谢审稿专家给出的建设性意见.