2017, Vol. 60
地震是应变能累积到一定程度后的释放过程，其中一种释放形式是以地震波向外辐射，这可通过地表的仪器观测到.而除了最终将衰减的地震波的传播，地球还会产生永久的变形，该变形会进一步产生重力位能的永久变化.重力位能的变化要比地震波能量大得多(Chao et al., 1995).实际上，在地震发生之前重力位能与流体静力平衡状态下的弹性应变能保持平衡状态，它们也是用位移表示的一阶能量；而地震发生后产生的弹性应变能、地震波能量都属于量级更小的二阶能量(Dahlen，1977).尽管我们无法察觉，但由于地球的永久变形，地球内部确实发生了重力位能与弹性应变能之间的能量转化.Chao等(1995)指出，地震有使全球重力位能减小的趋势.另一方面，Tanimoto和Okamoto(2000)的研究表明，拉张的构造活动使地壳的重力位能减小，而挤压的构造活动使地壳的重力位能增大.因此他们指出地壳的重力位能的同震变化能够很好地揭示构造活动的挤压或拉张状态.
重力位能的变化是由于地球物质在永久变形的位移下克服重力作用所作的功.Chao等(1995)应用Gilbert(1970)提出的由地球简正模计算地表位移的理论，发展了计算重力位能同震变化的简正模方法.该方法后来被Tanimoto等(2002)运用到计算地壳的重力位能变化中.然而，他们的结果表明：对于非常浅的地震，尽管简正模方法截断到低阶球谐项时对于全球重力位能同震变化的计算是收敛的，但在计算地壳的重力位能变化时，该方法是不收敛的.因此计算时需要截断到很高的球谐阶数，从而增大了很多的计算量.为了解决这个问题，Okamoto和Tanimoto(2002)提出了直接求解微分方程组边值问题的数值计算方法，即首先解边值问题获得垂直位移，然后对垂直位移进行积分运算从而得到重力位能的同震变化.实际上，只要我们在求解微分方程组的边值问题时增加一个变量就可以省略计算垂直位移再进行积分运算的中间变量和过程，从而将上述两个步骤合而为一.该变量在地表的值就是全球重力位能的变化；而该变量在地壳上下边界处值的差就是地壳的重力位能变化，具体内容详见本文的计算理论部分.
印度板块和欧亚板块的碰撞产生了世界上最高的青藏高原.青藏高原是全球构造活动最为活跃的地区之一.青藏高原的形成、演化和隆升机制一直是地球科学的研究热点(傅容珊等，2000；张健和石耀霖，2002；邓宾等，2016).由于印度板块向北的俯冲、碰撞，青藏高原的南北向呈明显的挤压状态，而东西向呈现出拉张状态.Tanimoto等(2002)研究了青藏高原地区1977—1998年的地震对地壳重力位能的影响.结果表明，高海拔地区的地震使地壳的重力位能减小，而低海拔地区的地震使重力位能增大.由此很容易看出地壳重力位能的这种变化与构造活动状态的对应关系.除此之外，他们也从地壳重力位能的变化角度给出了青藏高原东南缘重力坍塌的证据.
本文首先给出直接计算重力位能变化的理论；然后给出理论结果，并与Okamoto和Tanimoto(2002)的结果进行了比较；最后将该理论运用于青藏高原地区，为了了解最近十多年来该地区地震对地壳重力位能同震变化的影响，我们采用了1976—2013年的地震.
2 计算理论重力位能的同震变化可以表示为(Chao et al., 1995)
|
(1) |
其中ρ是密度，r是地球内部任一点的位置，g是地球重力场(由于不考虑地球的自转，也即引力场)，u是地震导致的同震位移场.对于SNREI地球，重力场是球对称的，因此(1) 式又可以写为
|
(2) |
其中，r是地球内部任一点的地心距，g是重力值(向下为正)，ur是垂直方向的同震位移，θ和φ分别是任一点的余纬和经度.由于对称性和球函数的正交关系，只有零阶的垂直位移对(2) 式中的积分才有贡献.Chao等(1995)发展了计算ΔGPE的简正模方法.在该方法中，垂直位移用简正模和地震矩表示，因此，(2) 式变为
|
(3) |
其中，
|
(4) |
R是地球半径，Mrr=M0sin2δsinβ是地震矩张量的一个分量，其中M0是标量地震矩，δ和β分别为断层的倾角和滑动角，rs是震源的地心距, ωn是本征频率，Un是与本征频率对应的垂直位移本征函数(即简正模)，In是规格化因子，
|
(5) |
显然，(3) 式的积分区间是地心到地表，那么所得即为全球重力位能变化；而如果积分区间只是整个地壳，即下地壳到地表，则所得为地壳重力位能变化.
必须注意的是，(4) 式中不可能对所有的无穷多的简正模求和，所以总要截断到某一球谐阶数.然而正如Tanimoto等(2002)指出的那样，对于震源深度较浅的地震，计算地壳重力位能时，(4) 式中的截断阶数需要很大，从而增大了计算量.因此，Okamoto和Tanimoto(2002)提出了直接的数值方法，即解下面的边值问题：
|
(6a) |
|
(6b) |
其边界条件为
|
(7) |
式中，y1表示垂直位移，y2表示垂向的正应力，μ和λ为弹性拉梅常数.需要注意的是，实际上(6b)式等号右边省略了震源函数项(Sun, 1992).由于该项包含Dirac函数，直接求解比较困难，因此为了简化计算，通过数学手段可将该项转移到震源处的边界条件中，从而对(6) 式求解两次得到该边值问题的解，即为一个特解和一个通解的组合.由该组合满足地表边界条件即可以获得地球内和表面任一点的垂直位移，然后利用(2) 式即可获得重力位能的变化.
尽管上述数值方法克服了简正模的收敛问题，但该方法包含了两个步骤.实际上，只要定义一个新变量就可以将这两步合并成一步.该变量定义为
|
(8) |
可以发现，该变量在地表的值即为全球重力位能变化，在上下地壳边界处的值之间的差即为地壳重力位能变化.利用零阶的微分方程组(Sun, 1992)，并顾及(8) 式，可以得到如下的边值问题：
|
(9a) |
|
(9b) |
|
(9c) |
|
(9d) |
其边界条件为
|
(10) |
(9c)式的加入是为了提供足够的边界条件.同样地，(9b)式等号右边也省略了震源函数项，通过相同的求解过程即可获得最终的结果.
根据位错理论(Sun, 1992), 对于SNREI地球，任意的地震都可以表示为四个独立的断层模式的线性组合，而其中只有水平拉张和垂直拉张型断层有零阶解.我们分别用上标数字22和33表示这两种断层模式.这两种断层的通解是一样的，可由(9) 式直接计算；它们的特解按如下的方法计算：取震源函数为初值，从震源积分到地表.上述两种断层模式的震源函数分别为
|
(11) |
其中，UdS是滑动量与断层面积的乘积，它与剪切模量的乘积即为标量地震矩M0.最后根据边界条件分别确定组合系数，从而得到这两种断层模式的解.对于任意的地震其解就按下面的方式组合：
对于剪切型地震：
|
(12) |
对于拉张型地震：
|
(13) |
通过比较(11) 式和(12) 式和Okamoto和Tanimoto(2002)一文中的(8) 式和(9) 式，我们发现他们采用了两种震源函数的差为特解的积分初值，所以他们的结果仅考虑了剪切型地震.Vavryuk(2011)研究了拉张地震的理论，通过增加一个滑动的倾向角(即滑动方向与断层面法向不垂直)，可以将任一断层分解为剪切和拉张分量.我们的理论不仅适合于剪切型地震，而且适合于拉张型地震，因此具有更大的应用价值.
3 数值结果与讨论 3.1 两种拉张断层产生的重力位能同震变化理论结果首先计算了两种独立的拉张断层引起的重力位能的同震变化，也就是y6.由于目前还没有观测到地震发生在深度700 km以下，因此仅计算了震源深度为1~700 km的结果.计算中，采用了PREM地球模型(Dziewonski and Anderson, 1981).为了消除PREM模型中的海洋层，我们将海洋层和最上层的地壳参数进行了平均.计算结果如图 1所示.图 1a和图 1b分别给出了不同深度的地震对全球和地壳重力位能的影响，图中的不连续性(突跳)是由于地球内部边界处密度的不连续性造成的.尽管垂直位移在震源处是不连续的，但它不会导致重力位能的不连续性.
|
图 1 不同震源深度的地震产生的全球和地壳重力位能变化(规格化系数M0) Fig. 1 Changes of global and crustal gravitational potential energy caused by earthquakes of varied focal depths (with normalization coefficient M0) |
由图 1a，我们发现水平拉张断层产生的全球重力位能变化要比垂直拉张断层产生的全球重力位能变化大.需要注意的是计算理论值时，取标量地震矩M0=1.如果考虑具体的地震，只需要将M0乘到结果中即可.由(12) 式中关于倾角和滑动角的表达式可以看出：如果滑动角在0°~180°之间变化(逆冲断层)，则该地震将使全球重力位能减小；如果滑动角在-180°~0°之间变化(正断层)，则该地震将使全球重力位能增大.通常大地震都发生在俯冲带，且多表现为逆冲型，也就是滑动角在0°~180°之间变化.因此这也就是为什么全球地震有使全球重力位能减小的趋势，正如Chao等(1995)所指出的那样.从图中也可以看出，因为两条曲线的差异随深度增大而减小，因此(震级相同的地震)震源越深，该地震引起的重力位能变化越小.
与上述结果正好相反，从图 1b可以看出：只要不是太浅的地震，垂直拉张断层产生的地壳重力位能变化要比水平拉张断层产生的地壳重力位能大.但如果震源深度小于4.7 km(所采用的地球模型不同，该值有少许的变化)，则结果正好相反.图中只给出了震源深度在1~50 km的结果，这是因为大地震的震源深度通常在这个范围之内.随着震源深度的增大，图中的曲线慢慢向上直到转换带，然后向下，但一直到700 km仍是大于0的.因此，当震源深度大于4.7 km时，正断层导致地壳重力位能减小，而逆冲断层导致地壳重力位能增大.通常正断层对应于拉张构造，而逆冲断层对应于挤压构造.所以Tanimoto和Okamoto(2000)指出地壳重力位能的同震变化能够很好地反映构造活动是拉张或是挤压的状态.尽管这种状态实际上完全决定于断层滑动角，但是采用地壳重力位能的同震变化相较于采用震源机制(震源小球)能够更加直观、并且能够定量地描述构造活动的状态.也必须看到，采用全球重力位能的同震变化同样也可以反映出构造活动的状态，只是其反映的构造活动的挤压或拉张状态与地壳重力位能变化正好相反.
另外，我们也可以发现图 1b中粗实线(标注为33-22) 与Okamoto和Tonimoto(2002)中的图 7除了有个比例因子2之外，完全一致.这是因为(12) 式前有一个系数1/2.
由(13) 式，以及基于图 1中标注为22和33的曲线的值都大于0的结果，可以得出如下的结论：所有的拉张型地震，不管地震发生在何处，断层的倾角是多大，它总是使全球和地壳重力位能变大.
3.2 青藏高原地区的结果我们计算了1976—2013年间青藏高原地区(15°N—50°N，50°E—115°E)的地震产生的地壳重力位能的同震变化.地震的震源参数搜集自GCMT(www.globalcmt.org)，在这一时间段内共记录了该区域发生的1891次地震.
首先将该区域划分为1°×1°的网格，计算在每个小网格内发生的地震对地壳重力位能变化的累积影响.用空心圆代表能量增大，实心圆代表能量减小，圆的半径与增大或减小的能量大小对应.最后将圆圈画在该小网格内.结果如图 2所示.需要说明的是：每个小网格内的小圆圈代表该网格内的地震对全球的地壳重力位能的影响，而不是对该小网格内的地壳重力位能的影响.
|
图 2 地震对地壳重力位能的影响 (a) 1976—2013年；(b)1976—1998年；(c)1999—2013年.图中也画出了青藏高原地区的5个主要活动断裂带：① 阿尔金；② 海原；③ 昆仑；④ 玛尼—玉树—鲜水河；⑤ 嘉黎(Gan et al., 2007). Fig. 2 Influence of earthquakes on crustal gravitational potential energy during different time periods (a) 1976—2013; (b) 1976—1998; (c) 1999—2013. Five major active fault zones in the Tibetan plateau: ① Altyn Tagh; ② Haiyuan; ③ Kunlun; ④ Mani-Yushu-Xianshuihe; ⑤ Jiali (Gan et al., 2007). |
由图中结果可以看出：在青藏高原的边缘地壳重力位能增大；而在其内部，除了东北部以外，地壳重力位能都是减小的.这说明边缘处于挤压状态，而内部处于拉张状态.对于这样的状态也有许多的证据予以证实，如震源机制结果(Molnar and Chen, 1983)和GPS结果(Wang et al., 2001；朱守彪等, 2005；Gan et al., 2007).Gan等(2007)不仅给出了该地区相对于欧亚板块的速度场，而且给出了其相对于青藏高原局部参考框架的速度场.这样人们能够更加直观地看出青藏高原的变形情况.从他们的结果中可以明显看出青藏高原中部的速度场从羌塘块体由西向东递增；速度场由中部向东南缘递增.这些都能够很好地解释图 2中反映出的拉张状态(实心圆).另外，地壳重力位能减小的地方基本在高程比较大的地方，这从图 2中绘出的高程阴影图也可直观看出.比较图 2b和2c可以看出，相比于1976—1998年期间，在1999—2013年期间青藏高原中部的拉张活动更加明显.
图中也画出了青藏高原地区5条主要的断裂带(以带圆圈的数字标出).从图中可以看出：阿尔金、海原和昆仑断裂带所处区域主要表现为挤压状态，而玛尼—玉树—鲜水河和嘉黎断裂带之间的区域表现为拉张状态.玛尼—玉树—鲜水河断裂带基本上将青藏高原分成东北和西南两块，其中东北部主要表现为挤压状态，而西南部主要表现拉张状态.
在青藏高原，另一个典型的构造活动是重力坍塌.一般认为它是因为高原物质由于隆升被抬升到很高的地方，然而底下没有足够的支撑，因此在重力作用下，物质向下坍塌.为此，我们给出了由于正断层(此处定义其滑动角范围为-135°~-45°)导致的重力位能变化，如图 3所示.结果表明正断层地震集中在青藏高原的西南和东南部，最近的GPS结果表明，尽管其他地区由于板块的碰撞继续隆升，然而这两个区域部分地区不是隆升，而是下沉(Liang et al., 2013).这说明重力坍塌是正断层地震发生的重要原因.青藏高原的东南部是物质向东运动受阻而转向东南流出的通道，因此物质从高处流出也是重力坍塌的原因.对比图 2b和2c，发现近十几年来青藏高原中部和西部地震使得更多的地壳重力位能减小.这也许意味着这些区域更多的重力坍塌.关于这一点，还没有足够的观测予以证实.然而通过比较也可发现，青藏高原东部的地震造成了更多的地壳重力位能增大，而且最近的青藏高原东部边缘的地震，如2008年的汶川地震，使得地壳重力位能明显增大.另外，图 2b和Tanimoto等(2002)中的图 3具有非常好的一致性.
|
图 3 正断层地震导致的地壳重力位能减小 Fig. 3 Sketch showing decrease of crustal gravitational potential energy caused by a normal-faulting earthquake |
通过定义一个新的变量，我们根据点源位错理论，提出了计算地球重力位能同震变化的新方法.该方法整合了以前计算采用的两步法，避免了利用垂直位移这一中间变量和计算积分过程，而是在求解边值问题时一并求出.并且该方法不仅适用于剪切型地震，而且适用于拉张型地震.在理论上改进了以前的数值解法或简正模方法只考虑剪切型地震的不足.
逆冲型地震使全球重力位能减小，而正断层地震使全球重力位能增大.通常大地震属于逆冲型，因此这种地震主导着全球重力位能的同震变化，所以地震有使全球重力位能减小的趋势.然而对于地壳重力位能来说，这种情况正好相反，除非震源深度小于4.7 km.但是事实上，这么浅的地震通常相对来讲震级都比较小，因此这些非常浅的地震对本文的结论没有太大的影响.
青藏高原是世界上构造活动最为活跃的地区之一，了解该地区的构造活动状态尤为重要.地壳重力位能的同震变化为我们提供了一种定量且非常直观的手段，通过它可以很直观地知道该地区或其部分地区的构造活动状态是挤压或拉张.在青藏高原的边界地区和东北部地区，具有明显的挤压运动；在青藏高原的中部和东南部，具有明显的拉张运动.我们也能清楚地识别青藏高原东南部的重力坍塌现象，且该现象已为GPS资料所证实.
| Chao B F, Gross R S, Dong D N. 1995. Changes in global gravitational energy induced by earthquakes. Geophys. J. Int., 122(3): 784-789. DOI:10.1111/gji.1995.122.issue-3 | |
| Dahlen F A. 1977. The balance of energy in earthquake faulting. Geophys. J. Int., 48(2): 239-261. DOI:10.1111/j.1365-246X.1977.tb01298.x | |
| Deng B, Yong Z Q, Liu S G, et al. 2016. Cenozoic mountain-building processes in the Daliangshan, southeastern margin of the Tibetan Plateau: Evidence from low-temperature thermochronology and thermal modeling. Chinese J. Geophys. , 59(6): 2162-2175. DOI:10.6038/cjg20160621 | |
| Dziewonski A M, Anderson D L. 1981. Preliminary reference earth model. Physics of the Earth and Planetary Interiors, 25(4): 297-356. DOI:10.1016/0031-9201(81)90046-7 | |
| Fu R S, Xu Y M, Huang J H, et al. 2000. Numerical simulation of the compression uplift of the Qinghai-Xizang plateau. Chinese J. Geophys. , 43(3): 346-355. DOI:10.3321/j.issn:0001-5733.2000.03.008 | |
| Gan W J, Zhang P Z, Shen Z K, et al. 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. J. Geophys. Res., 112(B8): B08416. DOI:10.1029/2005JB004120 | |
| Gilbert F. 1970. Excitation of the normal modes of the Earth by earthquake sources. Geophys. J. Int., 22(2): 223-226. | |
| Liang S M, Gan W J, Shen C Z, et al. 2013. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. J. Geophys. Res., 118(10): 5722-5732. DOI:10.1002/2013JB010503 | |
| Molnar P, Chen W P. 1983. Focal depths and fault plane solutions of earthquakes under the Tibetan plateau. J. Geophys. Res., 88(B2): 1180-1196. DOI:10.1029/JB088iB02p01180 | |
| Okamoto T, Tanimoto T. 2002. Crustal gravitational energy change caused by earthquakes in the western United States and Japan. Earth and Planetary Science Letters, 195(1-2): 17-27. DOI:10.1016/S0012-821X(01)00576-3 | |
| Sun W K. 1992. Potential and gravity changes raised by dislocations in spherically symmetric earth models [Ph. D. thesis]. Tokyo: The University of Tokyo. | |
| Tanimoto T, Okamoto T. 2000. Change of crustal potential energy by earthquakes: An indicator for extensional and compressional tectonics. Geophys. Res. Lett., 27(15): 2313-2316. DOI:10.1029/1999GL011122 | |
| Tanimoto T, Okamoto T, Terra F. 2002. Tectonic signatures in coseismic gravitational energy change. Geophys. J. Int., 149(2): 490-498. DOI:10.1046/j.1365-246X.2002.01662.x | |
| Vavryuk V. 2011. Tensile earthquakes: Theory, modeling, and inversion. J. Geophys. Res., 116(B12): B12320. DOI:10.1029/2011JB008770 | |
| Wang Q, Zhang P Z, Freymueller J T, et al. 2001. Present-day crustal deformation in China constrained by Global Positioning System measurements. Science, 294(5542): 574-577. DOI:10.1126/science.1063647 | |
| Zhang J, Shi Y L. 2002. The role of gravitational potential energy in raising and spreading of Qinghai-Xizang plateau. Chinese J. Geophys. , 45(2): 226-232. DOI:10.3321/j.issn:0001-5733.2002.02.009 | |
| Zhu S B, Cai Y E, Shi Y L. 2005. Computation of the present-day strain rate field of the Qinghai-Tibetan plateau and its geodynamic implications. Chinese J. Geophys. , 48(5): 1053-1061. DOI:10.3321/j.issn:0001-5733.2005.05.011 | |
| 邓宾, 雍自权, 刘树根, 等. 2016. 青藏高原东南缘大凉山新生代隆升建造过程——多封闭系统低温热年代学与热模型限制. 地球物理学报, 59(6): 2162–2175. DOI:10.6038/cjg20160621 | |
| 傅容珊, 徐耀民, 黄建华, 等. 2000. 青藏高原挤压隆升过程的数值模拟. 地球物理学报, 43(3): 346–355. DOI:10.3321/j.issn:0001-5733.2000.03.008 | |
| 张健, 石耀霖. 2002. 青藏高原隆升及伸展变形中的重力位能. 地球物理学报, 45(2): 226–232. DOI:10.3321/j.issn:0001-5733.2002.02.009 | |
| 朱守彪, 蔡永恩, 石耀霖. 2005. 青藏高原及邻区现今地应变率场的计算及其结果的地球动力学意义. 地球物理学报, 48(5): 1053–1061. DOI:10.3321/j.issn:0001-5733.2005.05.011 | |