利用震后GPS数据反演汶川地区有效黏滞系数

朱良玉 王双绪 蒋锋云

（中国西安710054中国地震局第二监测中心）

引言

近年来，随着大地测量技术的发展，许多研究人员开始利用GPS观测资料研究地球介质的黏性结构.孙荀英等（1994）采用三维黏弹性有限元方法拟合1976年唐山MS7.8地震同震和震后形变，在国内首次开展了利用震后形变资料研究地球黏滞系数的工作，其结果认为华北板块下方软流层的黏滞系数为7×1018Pa·s.王庆良和巩守文（1997）采用1990年共和MS6.9地震后的水准资料，反演得到共和震区下地壳有效黏滞系数为1018Pa·s.朱守彪和蔡永恩（2006）采用有限元软件，根据1999年集集MS7.6地震GPS震后时间序列反演了该区域下地壳黏滞系数约为3.6×1019Pa·s.张晁军等（2008）采用3种不同时段的方法计算了青藏高原下地壳的黏滞系数，量级在1017—1021Pa·s之间不等，并认为下地壳的黏滞系数随着空间和时间的变化均有相应的变化.郝明等（2011）采用最小二乘法对1976年唐山MS7.8地震、1966年邢台MS7.2地震、1970年通海MS7.8地震和1973年炉霍MS7.3地震的震后垂直位移进行拟合，计算得到这4次地震的有效黏滞系数均在1017—1019Pa·s量级.这些方法主要采用震后观测到的位移时间序列，通过震后形变模型来达到反演地壳黏滞系数的目的，以期能够获得研究区站点的介质信息.但是由于地球介质的复杂性，仅仅通过少数连续台站观测数据来获取整个研究区域的黏滞系数，显然是不够充分的.另外在实际工作中，由于观测能力的限制，往往很难在震区布设足够密集的连续GPS观测台站.因此，如何采用震区大量流动观测台站数据反演研究区的黏滞系数，对于厘清震后形变模型和震区介质信息具有重要意义.

2008年汶川地震发生之后，许多学者对汶川地震的发震机理（王庆良等，2008；张培震等，2009；滕吉文等，2010；蒋锋云等，2011）以及断层参数和滑动分布（陈运泰，2008；Ji，Hayes，2008；Shen et al，2009；谈洪波等，2009；王敏，2009；许才军等，2009；徐锡伟等，2010；Wang et al，2011）进行细致深入的研究，取得了显著成果.这为研究汶川地震后的地壳运动模式和介质物性参数奠定了基础.因此，基于上述考虑，本文拟采用汶川震后2009—2011年区域GPS速度场数据，结合三维半空间准动态黏弹性松弛模型（Wang et al，2006）反演该地区地壳介质的黏滞系数.

1 数据源及正演模型

本文采用中国地壳运动观测网络2009和2011年两期GPS观测数据，通过GAMIT／GLOBK软件解算出研究区域的GPS速度场（相对于欧亚块体），如图1所示.由于选取的震后时间较短，震后松弛范围较小，结合汶川震区的断层分布特点，本文选取了如图1所示的龙门山断裂区域作为研究区域.该区域包括35个GPS区域观测站，每个观测站均有2009和2011年两期观测数据.

本文采用的正演模型是在Lorenzo-Martin等（2006）提出的震后形变模型基础上，经过修正之后的模型.Lorenzo-Martin等在研究1966年智利MS9.5地震时，认为震后观测的地壳形变包含了板块汇聚产生的形变、震后黏弹性松弛及断层面上一些无震滑动产生的形变.其模型为

图1 观测的GPS水平形变场（黑色箭头）与模拟的GPS形变场（红色箭头）.红色和黑色圆圈表示相应的误差椭圆Fig.1 The measured（black arrows）and modeled（red arrows）GPS velocities where black and red circles denote corresponding errors

式中，Vmod为模型计算的地壳形变，Vpla为活动块体汇聚所产生的形变，Vrel为汶川地震引起的震后黏弹性松弛所产生的形变，Vss为震后无震滑动产生的形变，Ucos为同震滑动分布，Upos为震后滑动分布，χ为耦合因子，η为黏滞系数.该模型认为在断层一侧稳定不动的情况下，运动板块汇聚所产生的应变能会产生两个作用：一是板块运动块体分量；二是断层无震滑动分量.因此需要一个耦合因子来调和震后形变中块体运动分量与无震滑动分量.而本文研究的汶川地震是板内地震，其观测数据遍布断层两侧，不同于Lorenzo-Martin等（2006）研究的智利地震.因此本文对式（1）所示的震后形变模型进行了修正，将块体运动分量与无震滑动分量看作相互独立的两个分量，采用的模型如下：

1.1 构造背景设定

龙门山断裂位于巴颜喀拉地块东边界.由于受到该地块持续向东的构造加载，该地区观测的GPS水平形变场必然包含了该块体的一个整体性运动.根据龙门山断裂以西距离断裂200km范围内的23个GPS观测站，采用整体无旋转基准计算得到该区域的刚性块体运动，表示如下：

式中，R为地球半径，本文取为6 371km；λ，φ分别为数据点的经、纬度；wx，wy，wz为采用龙门山断裂以西所选数据通过最小二乘法计算得到的欧拉极三分量.

对于龙门山断裂以东区域，GPS观测点相对于华南地块运动较小，基本上可看成一个整体.因此，将解算得到的GPS观测转换到相对于华南块体的运动框架下，经过下式得到块体的运动速度：

式中，R为地球半径；λ，φ为数据点的经、纬度；λ0，φ0为华南块体质心位置的经、纬度；u，α，w分别为块体质心运动的大小、方向及旋转角速度.此处采用王敏（2009）给出的华南地块运动参数，即λ0＝110.71°，φ0＝27.66°，u＝6.6mm，α＝108.3°，w＝0.2×10－9rad／a.

1.2 断层破裂参数设定

由式（2）可知，震后黏弹性松弛和震后非震滑动均与同震破裂参数密切相关，因此，断层同震破裂参数对于计算汶川地震的震后松弛至关重要.本文的同震破裂参数引自Ji和Hayes（2008）一文.根据赵珠等（1997）给出的龙门山断裂带P波和S波的地壳平均速度模型，王椿镛等（2002）反演获得的川滇及邻区三维地壳速度结构，以及Wang（2001）对中国大陆热流结构和强度的研究结果，同时参考张晁军等（2008）给出的炉霍地区地壳模型，建立了该地区的地壳结构模型，如表1所示.研究区介质分为上下两层，弹性层覆盖于麦克斯威尔体上，黏弹性松弛产生的震后形变为

其中，X2009，X2011分别表示2009和2011年发生的震后形变，η为黏滞系数，Ucos为汶川地震的同震滑动分布.

1.3 震后余滑参数的设定

董运洪等（2012）根据汶川震后水准资料研究发现，该区存在快速滑动背景，因此在模拟汶川地震震后滑动时必须考虑余滑效应所产生的震后形变.震后余滑作为与同震紧密相关的一种震后运动形式，与同震滑动分布密切相关.一般认为，在发震断层面上，同震滑动较大且断层摩擦系数较小的区域，更容易发生余滑.本文采用许才军等（2009）的同震滑动分布，将汶川地震余滑断层粗略地分为东北段和西南段，通过下式计算得到龙门山断裂带断层面上2009—2011年的3个滑动分量：

表1 汶川地震的地壳介质模型Table 1 Parameters of crustal material model for Wenchuan earthquake

式中，α为断层的方位角；δ为倾角；Us，Ud，Uo分别为断层面滑动三分量，此处仅考虑断层面上的剪切分量Us和Ud，未考虑拉张分量Uo，计算结果详见表2；ΔVN，ΔVE分别为断层上下盘北向、东向速度的平均值之差；ΔVV为龙门山断裂带上2009—2011年北西盘相对于南东盘的垂直形变，由董运洪等（2012）给出的水准数据计算得到.

表2 汶川地震震后龙门山断裂的余滑参数Table 2 After-slip parameters of Longmenshan fault after Wenchuan earthquake

这样，得到断层的震后滑动速率分布后，采用Sun等（1996）给出的球面位错公式（式（7））估算震后断层运动在断层左右两盘产生的弹性形变.

式中，Upos为震后断层位移； 为断层位移与产生形变之间的格林函数，具体形式请参考Sun等（1996）一文.

2 反演结果

将式（3）、（5）、（7）带入式（2）得到龙门山断裂以西的模拟震后水平形变.同理，将式（4）、（5）、（7）带入式（2）可得到龙门山断裂以东四川盆地的震后水平形变.对于给定的断层参数，设定不同的黏滞系数，即可得到不同的震后形变场.为了反演断层区域的黏滞系数，本文采用统计量Δ为

式中，N为观测点数目；VmodE和VmodN分别为模拟计算得到的震后形变场的东西分量和南北分量；VmodE，VmodN，σmodE，σmodN分别为观测所得的震后东西向、南北向速度及其相应误差.当Δ最小时，即认为取得最佳的黏滞系数.由于黏滞系数η与统计量Δ之间为非线性关系，本文采用区间法，通过先验信息确定中下地壳黏滞系数η的搜索区间一般为1018—1025Pa·s，然后通过搜索找出使Δ最小的η值.

图2 残差加权平方和VTPV（Δ）与中下地壳黏滞系数η之间关系Fig.2 Residuals VTPV （Δ）versus effective viscosityηof middle to lower crust

如图2所示，通过搜索发现，Δ随介质黏滞系数η的变化没有规律，仅有两个底值区域；随着黏滞系数的增加，残差平方和不再变化，无法得到合理的反演结果.值得注意的是，通过式（5）计算黏弹性松弛产生的震后形变时，没有考虑介质的横向非均匀性.实际上，断层两侧的介质特性严重影响着地表观测的分布形态.因此，为了得到更佳的模拟结果，对断层两侧采用不同的黏滞系数值来计算该区域震后松弛产生的形变.龙门山断裂以西的形变场表示如下：

而对于龙门山断裂以东的四川盆地，则采用

分别将式（3），（9），（7）与式（4），（10），（7）带入式（2），以式（8）为目标函数，采用格网搜索法所得结果如图3与图4所示.

从图3可以看出，残差随着黏滞系数的增大而逐渐减小.当黏滞系数大于7×1022Pa·s时，ΔE不再减小而保持稳定值.这说明增加黏滞系数确实能有效改善拟合效果.但黏滞系数增加到一定程度之后，黏弹性松弛效应对数据拟合效果已不能进一步改进.因此，根据本文的正演模型，反演得到四川盆地16—40km地壳的黏滞系数为7×1022Pa·s.龙门断裂以东取得最佳反演结果时的模拟结果详见表3.

表3 龙门山断裂以东与观测GPS速度场最佳拟合结果Table 3 The best-fitting-observed results of GPS velocity inversion in eastern side of Longmenshan fault 单位：mm·a－1

图4给出了龙门山断裂以西川西块体残差加权平方和VTPV（ΔW）与中下地壳黏滞系数ηW之间的关系.从该图可以看出，该区域16—40km地壳的黏滞系数在（7.0—9.5）×1019Pa·s范围内时，可使得残差最小.因此我们认为龙门山断裂以西川西高原中下地壳的黏滞系数大致在（7.0—9.5）×1019Pa·s范围内.但仍需注意的是，在图4所示的区间内，偏差值依然存在波动，并未收敛于一个极值点上.这说明川西地区中下地壳黏滞系数的横向变化依然很强烈，单一黏滞系数值仅仅体现的是该区域黏滞系数的平均值.从图1的GPS速度场分布也可以发现，龙门山断裂以西岷山地块附近拟合效果欠佳，这与图4所反映的多黏滞系数值现象是一致的.表4给出了龙门山断裂以西最佳拟合结果的反演数据.

表4 龙门山断裂以西与观测GPS速度场最佳拟合结果Table 4 The best-fitting observed results of GPS velocity inversion to the western side of Longmenshan fault 单位：mm·a－1

3 讨论与结论

由于本文的计算模型比较复杂，为了粗略地评估反演结果，忽略了计算过程中对误差的放大，将观测数据与模拟数据之差的绝对值与观测数据的原始误差之和作为模拟数据的误差.将模拟结果加上误差椭圆之后显示在图1上.从图1中可以看出，GPS形变场的拟合效果总体上较好，能够反映该区域的运动特征，但是在一些区域存在较大差异.

1）从图1中的红色误差椭圆上可以看出，岷江断裂以东的观测区域内拟合误差较大，拟合结果比实际观测量级大很多，尤其如H021，H005等测点模拟观测与实际观测基本无法匹配.这无疑说明，川西地区中下地壳黏滞系数的横向变化很强烈，单一黏滞系数值并不能在岷江断裂以东、龙门山断裂以北的区域取得较好的反演效果.实际上，上述区域的地壳水平运动比龙门山山前区域小，但比龙门山断裂以东的四川盆地大.这与本文得到的龙门山断裂以西岷江地块中下地壳的黏滞系数比龙门山山前区域大而比四川盆地小的结论相一致.从图4可以推测岷江地块中下地壳黏滞系数在1020Pa·s左右，而川西地区中下地壳黏滞系数约为7×1019Pa·s.

2）龙门山主断裂灌县—江油段附近的观测点H045、H049、H050、H058及H384模拟效果最差，最大误差达20mm（图1）.其主要原因在于，汶川震中附近断裂由于震前闭锁，在同震破裂后存在一个显著的滑脱面（见张培震等（2009）中图2），因而在震后短时间内发震断层上余滑效应显著（表2）.本文在计算余滑分量时选取了相对平均的断层滑动，不能完全消去余滑产生的效果.因此，震中附近观测点的拟合效果往往比较差.

本文根据汶川地区的地质背景按断层两侧将研究区域分开反演.其中龙门山断裂以西16—40km中下地壳的有效黏滞系数为（7.0—9.5）×1019Pa·s，其松弛时间约为31—41年.这与朱守彪和蔡永恩（2006）、张晁军等（2008，2009）、杨强和党亚明（2010）等反演的青藏高原中下地壳黏滞系数大体相当，均为1017—1020Pa·s量级；龙门山断裂以东四川盆地16—40km中下地壳的有效黏滞系数为7×1022Pa·s，其松弛时间约为3万年，呈现极强的整体性运动.石耀霖和曹建玲（2008）通过实验室流变结果估算了中国大陆岩石圈的等效黏滞系数，其中50km地壳深度青藏高原等效黏滞系数为1020Pa·s，而其它地区由于地壳较薄，等效黏滞系数为1021—1022Pa·s.这与本文结论基本一致，从侧面验证了本文研究结果的可靠性.

当前采用大地测量资料反演地壳的有效黏滞系数，主要是利用GPS时间序列资料.但由于观测成本的限制，连续站观测资料台站数目较少，对地壳黏滞系数的空间约束条件较差.本文利用大量的区域流动站观测资料，大大增加了反演地壳黏滞系数的约束条件，使得反演结果更能体现区域的黏滞系数特性.通过反演计算，我们认为龙门山断裂以西川西高原中下地壳（16—40km）的黏滞系数初步估算为（7.0—9.5）×1019Pa·s，而四川盆地16—40km中下地壳的黏滞系数约为7×1022Pa·s.

感谢德国波兹坦地球科学中心汪荣江教授提供PSGRN／PSCMP程序，感谢审稿专家对本文提出中肯的修改意见.

陈运泰.2008.汶川特大地震的震级和断层长度［J］.科技导报，26（10）：26－－27.

Chen Y T.2008.On the magnitude and the fault length of the great Wenchuan earthquake［J］.Science and Technology Review，26（10）：26－－27（in Chinese）.

董运洪，罗三明，韩月萍，陈长云.2012.汶川MS8.0地震北川附近同震和震后垂直位移特征［J］.地震学报，34（5）：611－－620.

Dong Y H，Luo S M，Han Y P，Chen C Y.2012.Co-and post-seismic vertical displacements of Wenchuan MS8.0earthquake near Beichuan［J］.Acta Seismologica Sinica，34（5）：611－－620（in Chinese）.

郝明，王庆良，占伟，李煜航.2011.中国东西部地区震后形变有效松弛时间研究［J］.大地测量与地球动力学，31（6）：49－－53.

Hao M，Wang Q L，Zhan W，Li Y H.2011.On effective relaxation time of post-seismic deformation in eastern and western China［J］.Journal of Geodesy and Geodynamics，31（6）：49－－53（in Chinese）.

蒋锋云，王庆良，朱良玉.2011.汶川地震震前龙门山以西地壳垂直隆升机理分析［J］.大地测量与地球动力学，31（5）：26－－29.

Jiang F Y，Wang Q L，Zhu L Y.2011.Mechanism analysis of vertical uplift of crust in the west of Longmenshan before Wenchuan MS8.0earthquake［J］.Journal of Geodesy and Geodynamics，31（5）：26－－29（in Chinese）.

石耀霖，曹建玲.2008.中国大陆岩石圈等效粘滞系数的计算和讨论［J］.地学前缘，15（3）：82－－95.

Shi Y L，Cao J L.2008.Effective viscosity of China continental lithosphere［J］.Earth Science Frontiers，15（3）：82－－95（in Chinese）.

孙荀英，刘激扬，王仁.1994.1976年唐山地震震时和震后变形的模拟［J］.地球物理学报，37（1）：45－－55.

Sun X Y，Liu J Y，Wang R.1994.Simulation of coseismic and postseismic deformation in Tangshan earthquake of 1976［J］.Chinese J Geophys，37（1）：45－－55（in Chinese）.

谈洪波，申重阳，李辉，李冀，玄松柏，邢乐林.2009.汶川大地震震后重力变化和形变的黏弹分层模拟［J］.地震学报，31（5）：491－－505.

Tan H B，Shen C Y，Li H，Li J，Xuan S B，Xing L L.2009.Simulation of post-seismic gravity change and deformation of the Wenchuan earthquake based on viscoelastic layered half-space model［J］.Acta Seismologica Sinica，31（5）：491－－505（in Chinese）.

滕吉文，张永谦，阮小敏.2010.汶川8.0级大地震的孕育和发生的深层过程与动力成因探讨［J］.国际地震动态，（6）：2.

Teng J W，Zhang Y Q，Ruan X M.2010.Deep processes and dynamic responses associated with the Wenchuan MS8.0 earthquake of 2008［J］.Recent Developments in World Seismology，（6）：2（in Chinese）.

王敏.2009.基于GPS同震位移场约束反演2008年5·12汶川大地震破裂空间分布［J］.地球物理学报，52（10）：2519－－2526.

Wang M.2009.Coseismic slip distribution of the 2008Wenchuan great earthquake constrained using GPS coseismic displacement field［J］.Chinese J Geophys，52（10）：2519－－2526（in Chinese）.

王椿镛，Mooney W D，王溪莉，吴建平，楼海，王飞.2002.川滇地区地壳上地幔三维速度结构研究［J］.地震学报，24（1）：1－－16.

Wang C Y，Mooney W D，Wang X L，Wu J P，Lou H，Wang F.2002.3Dvelocity structure of lower crust and upper mantle in Chuan－－Dian area［J］.Acta Seismologica Sinica，24（1）：1－－16（in Chinese）.

王庆良，巩守文.1997.由1990年共和7.0级地震震后垂直形变求得的地球介质有效弛豫时间和粘滞系数［J］.地震学报，19（5）：480－－486.

Wang Q L，Gong S W.1997.Estimating the effective viscosity from the vertical deformation in Gonghe earthquake of 1990［J］.Acta Seismologica Sinica，19（5）：480－－486.

王庆良，崔笃信，王文萍，张四新，刘锦文，史旗.2008.川西地区现今垂直地壳运动研究［J］.中国科学：D辑，38（5）：598－－610.

Wang Q L，Cui D X，Wang W P，Zhang S X，Liu J W，Shi Q.2008.Investigating vertical crustal deformation of western Sichuan area present［J］.Science in China：Series D，38（5）：598－－610（in Chinese）.

徐锡伟，陈桂华，于贵华，孙鑫喆，谭锡斌，陈立春，孙建宝，陈于高，陈文山，张淑萍，李康.2010.5·12汶川地震地表破裂基本参数的再论证及其构造内涵分析［J］.地球物理学报，53（10）：2321－－2336.

Xu X W，Chen G H，Yu G H，Sun X Z，Tan X B，Chen L C，Sun J B，Chen Y G，Chen W S，Zhang S P，Li K.2010.Reevaluation of surface rupture parameters of the 5·12Wenchuan earthquake and its tectonic implication for Tibetan uplift［J］.Chinese J Geophys，53（10）：2321－－2336（in Chinese）.

许才军，刘洋，温阳茂.2009.利用 GPS资料反演汶川 MW7.9级地震滑动分布［J］.测绘学报，38（3）：195－－201，215.

Xun C J，Liu Y，Wen Y M.2009.MW7.9Wenchuan earthquake slip distribution inversion from GPS measurements［J］.Acta Geodaetica et Cartographica Sinica，38（3）：195－－201，215（in Chinese）.

杨强，党亚明.2010.利用GPS速度场估算青藏高原地壳韧性层等效粘滞系数分布的研究［J］.测绘学报，39（5）：497－－502.

Yang Q，Dang Y M.2010.A research about effective viscosity of Tibetan Plateau lithosphere viscoelastic ductile layer using GPS velocity fields［J］.Acta Geodaetica et Cartographica Sinica，39（5）：497－－502（in Chinese）.

张晁军，曹建玲，石耀霖.2008.从震后形变探讨青藏高原下地壳粘滞系数［J］.中国科学：D辑，38（10）：1250－－1257.

Zhang C J，Cao J L，Shi Y L.2008.Exploring viscosity of lower crust from postseismic deformation［J］.Science in China：Series D，38（10）：1250－－1257（in Chinese）.

张晁军，石耀霖，马丽.2009.昆仑山大地震震后形变反映的地壳岩石流变特性［J］.岩土力学，30（9）：2552－－2558.

Zhang C J，Shi Y L，Ma L.2009.Numerical simulation of crust rheological property reflected by post-seismic deformations of Kunlun large earthquake［J］.Rock and Soil Mechanics，30（9）：2552－－2558（in Chinese）.

张培震，闻学泽，徐锡伟，甘卫军，王敏，沈正康，王庆良，黄媛，郑勇，李小军，张竹琪，马胜利，冉勇康，刘启元，丁志峰，吴建平.2009.2008年汶川8.0级特大地震孕育和发生的多单元组合模式［J］.科学通报，54（7）：944－－953.

Zhang P Z，Wen X Z，Xu X W，Gan W J，Wang M，Shen Z K，Wang Q L，Huang Y，Zheng Y，Li X J，Zhang Z Q，Ma S L，Ran Y K，Liu Q Y，Ding Z F，Wu J P.2009.2008Wenchuan 8.0earthquake preparation and occurrence of multi-unit combined mode［J］.Chinese Sci Bull，54（7）：944－－953（in Chinese）.

赵珠，范军，郑斯华，长谷川昭，堀内茂木.1997.龙门山断裂带地壳速度结构和震源位置的精确修订［J］.地震学报，19（6）：615－－622.

Zhao Z，Fan J，Zheng S H，Hasegawa A，Horiuchi S.1997.Crustal structure and accurate hypocenter determination along the Longmenshan fault［J］.Acta Seismologica Sinica，19（6）：615－－622（in Chinese）.

朱守彪，蔡永恩.2006.利用GPS观测的时间序列资料反演地壳地幔粘性结构［J］.地球物理学报，49（3）：771－－777.

Zhu S B，Cai Y E.2006.Inversion of viscous properties of crust and mantle from the GPS temporal series measurements［J］.Chinese J Geophys，49（3）：771－－777（in Chinese）.

Ji C，Hayes G.2008.Preliminary result of the May 12，2008 MW7.9eastern Sichuan，China earthquake［EB／OL］.［2008-10-12］.http：∥earthquake.usgs.gov／earthquakes／eqinthenews／2008／us2008ryan／finite＿fault.php.

Lorenzo-Martin F，Roth F，Wang R J.2006.Inversion for rheological parameters from post-seismic surface deformation associated with the 1960Valdivia earthquake，Chile［J］.Geophys J Int，164（1）：75－－87.

Shen Z K，Sun J B，Zhang P Z，Wan Y G，Wang M，Bürgmann R，Zeng Y H，Gan W J，Liao H，Wang Q L.2009.Slip maxima at fault junctions and rupturing of barriers during the 2008Wenchuan earthquake［J］.Nature Geoscience，2（10）：718－－724.

Sun W K，Okubo S，Peter V.1996.Global displacements caused by point dislocations in a realistic Earth model［J］.J Geophys Res，101（B4）：8561－－8577.

Wang Q，Qiao X J，Lan Q G，Jeffrey F，Yang S M，Xu C J，Yang Y L，You X Z，Tan K，Chen G.2011.Rupture of deep faults in the 2008Wenchuan earthquake and uplift of the Longmen Shan［J］.Nature Geoscience，4（9）：634－－640.

Wang R J，Lorenzo-Martin F，Roth F.2006.PSGRN／PSCMP：A new code for calculating co-and post-seismic deformation，geoid and gravity changes based on the viscoelastic-gravitational dislocation theory［J］.Computers ＆ Geosciences，32：527－－541.

Wang Y.2001.Heart flow pattern and lateral variations of lithosphere strength in China Mainland：Constraints on active deformation［J］.Phys Earth Planet Inter，126（3／4）：121－－146.