采用三维空间数据计算汶川地区地壳应变

刘明学 段虎荣

1 陕西铁路工程职业技术学院，渭南市站北街东1号，714000

2 西安科技大学测绘科学与技术学院，西安市雁塔路58号，710054

传统地壳应变计算方法是三角形法［1－3］，该方法主要考虑X、Y两个方向的主应变和XY平面的剪应力，其计算结果反映区域平面应变趋势［4］。也有学者对四边形单元法开展相关研究［5－7］。石耀霖对球面坐标系计算应变进行研究［8］。武艳强采用球谐函数整体解算了GPS应变场［9］。但以上都是研究平面的力学问题。本文在前人的基础上，利用2005～2008的GPS观测资料，采用四边形计算法对汶川地区地壳运动和区域应变场特征进行分析。应变计算过程中引入高程分量，提出四边形网络构造，推导X、Y、Z3 个方向及XY、XZ、YZ3个面的地壳应变计算公式，分析了汶川地区应变特征，不仅能反映该区域平面的应变结果，还能反映垂直应变变化趋势。

图1 四边形单元示意图Fig.1 Schematic diagram of quadrilateral element

1 应变计算方法

假定两个GPS观测点之间是线性变化的，根据观测网建立四边形图形网（空间直角坐标系），通过同一个四边形不同方向上的变化计算出各边长方向的线应变，再用不同方向的线应变计算X、Y、Z3个坐标轴方向的主应变和XY、XZ、YZ3个面的剪应变，最后得出应变分布的计算公式。在引入高程的基础上，对四边形进行划分（图1），利用位移场直接计算应变。

图1为GPS观测点在XY平面上的投影，其中，X轴指向正北方向，Y轴指向正东方向，1、2、3…8为四边形单元的顶点，x、y、z为各点的三维空间位置，u、v、w为各点的三维位移量。四边形单元4个角点的三维坐标分别为（x1，y1，z1）、（x2，y2，z2）、（x3，y3，z3）、（x4，y4，z4），4个角点的位移观测值分别为（u1，v1，w1）、（u2，v2，w2）、（u3，v3，w3）、（u4，v4，w4）。该四边形内任意一点（x，y，z）的位移场（u，v，w）是坐标的线性函数，可表示为：

将四边形的4个角点x轴方向的位移和坐标值代入式（1）中第1式，得：

通过解式（2），可求出系数a1、a2、a3、a4，同理可解得系数b1、b2、b3、b4、c1、c2、c3、c4。

主应变计算公式为：

剪应变计算公式为：

2 观测数据

国家测绘地理信息局、国家重点基础研究发展计划项目“活动地块边界带的动力过程与强震预测”项目组、国家重大科学工程“中国地壳运动观测网络”在龙门山断裂带两侧布设了一系列数据监测点，采用GPS流动观测为主，测量了监测点的水平位移。中国地震局在2007－04～07进行了多期观测，观测周期为3～4d。汶川地震发生后，又对这些GPS 点进行了复测，观测周期为2～3d。国家测绘地理信息局在2005～2008－05－12按等级的不同，采用1～4d的观测周期进行观测，震后又对其进行周期为1d的复测。四川省地震局和重庆市建立以服务为主的GPS 连续观测网络，提供了垂向同震位移资料。本文采用双差模式，数据采样间隔为30s，24h为一时段，用GAMIT／GLOBK 软件对上述单位的观测资料进行处理。

3 汶川地区应变特征计算分析

分析研究区域GPS 观测点的分布位置（图2），对观测点进行四边形网格化，构建出72个四边形，如图3所示。

根据式（3）、（4）解算每个四边形应变参数，求取X、Y、Z方向的主应变和XY、XZ、YZ面的剪应变，并作出四边形对应等值线图（图4～9）。

图4中粗黑实线表示断层（下同）。由图4可看出，X方向主应变值处于［－0.3，0.75］，以龙门山断裂带为界，龙门山断裂带以西为负应变，以东为正应变，在成都以南和映秀以西附近分别出现应变峰值，映秀、汶川、北川均处于应变高梯度带区域，说明映秀以西附近地区在地震中有向东移动的现象。该结果与实际位移图中映秀以西区域向东移动，成都、都江堰区域向西北移动相对应。

图5显示青川及附近区域负应变比较大，成都附近和黑水西北方向分别出现正应变峰值。对比实际位移图像，青川在Y轴沿负方向移动量相当大，而映秀附近则和断裂带上该区域大多数移动方向一致，向Y轴的正方向移动。汶川、北川均处于应变高梯度带区域。映秀、汶川、北川处于挤压状态，Y方向主应变值处于［－28，12］，远大于X方向。

图2 汶川地区断层及GPS观测点分布图Fig.2 The distribution of faults and GPS

图3 研究区域四边形网格划分Fig.3 The quadrilateral mesh of regional

图4 X 方向主应变等值线图Fig.4 The contour map of principal strain in Xdirection

图5 Y 方向主应变等值线图Fig.5 The contour map of principal strain in Ydirection

从图6可以看出，Z方向的主应变值处于［－0.002，0.002 8］，以龙门山断裂带为界，龙门山断裂带以西为负应变，以东为正应变，在成都以南和茂县以西附近分别出现应变峰值，映秀、汶川、北川均处于应变高梯度带区域。

图6 Z 方向主应变等值线图Fig.6 The contour map of principal strain in Zdirection

图7中XY面剪应变的值处于［－20，12］，高值区出现在映秀－北川－青川断裂带附近，表明这片区域的剪应变值偏大，处于挤压状态。该状态与震后的地质结果一致，这一带发生右旋剪切应变，全长超过300km，地表破裂带沿着映秀－北川走向断断续续分布，破裂带切割了多种地貌单元，包括河流、山脉、冲洪积扇、桥梁和公路等，使道路发生曲拱，桥梁垮塌、位移。

图7 XY 面剪应变等值线图Fig.7 The contour map of shear strain in XYdirection

由图8可以看出，XZ面剪应变值处于［－0.3，0.45］，负高值区出现在北川附近，小金、宝兴附近出现较大正剪应变，汶川、映秀处于梯度带区域。

图9中YZ面的剪应变值为［－5，9］，可以看出，平武、梓潼附近出现高正剪应变，宝兴、邛崃、成都一带出现负剪应变，北川、汶川、映秀处于剪应变梯度带区域。

图8 XZ 面剪应变等值线图Fig.8 The contour map of shear strain in XZdirection

图9 YZ 面剪应变等值线图Fig.9 The contour map of shear strain in YZ direction

从图4～6可以看出，Y方向主应变最大，X方向次之，Z方向最小，北川、汶川均处于主应变的高梯度带区域。从图7～9可以看出，XY面剪应变最大，YZ面次之，XZ面最小，汶川、映秀均处于剪应变梯度带区域。汶川、北川地区处于水平挤压应变状态，西北地区垂直分量处于向上拉张状态，东南地区垂直分量处于向下挤压状态，与文献［10－11］结果一致。

4 结 语

本文基于汶川震前的GPS定期观测数据与震后第一时间的复测数据，计算汶川地区地壳应变问题，分析汶川地区地壳变形的基本特征，得到如下结论：1）平面应变结果显示汶川、映秀地区处于东西方向挤压状态，高程结果显示在主断裂带东南方向应变为正值，西北方向应变为负值，即东南地区处于向下挤压状态，西北地区处于向上拉张状态；2）Y方向主应变最大，X方向次之，Z方向最小，北川、汶川均处于主应变的高梯度带区域；3）XY面剪应变最大，YZ面次之，XZ面 最小，汶川、映秀均处于剪应变梯度带区域，汶川地震恰好发生在应变高梯度带。

致谢：感谢国家测绘地理信息局、国家重点基础研究发展计划项目“活动地块边界带的动力过程与强震预测”项目组为本文研究提供观测数据。

［1］孟国杰，申旭辉，伍吉仓，等.GPS揭示的贝加尔湖地区现今地壳形变特征［J］.大地测量与地球动力学，2006，26（1）：15－20（Meng Guojie，Shen Xuhui，Wu Jicang，et al.Present－Day Crustal Deformation Revealed by GPS in Baikal Lake Region，Russia［J］.Journal of Geodesy and Geodynamics，2006，26（1）：15－20）

［2］肖根如，甘卫军，陈为涛.地应变计算Delaunay三角网在MATLAB与GMT 环境下的相互转换［J］.大地测量与地球动力学，2010，30（3）：122－126（Xiao Genru，Gan Weijun，Chen Weitao.Exchange of Delaunay Tin between MATLAB and GMT in Crustal Strain Calculation［J］.Journal of Geodesy and Geodynamics，2010，30（3）：122－126）

［3］吕志鹏，伍吉仓，孟国杰.条件数对三角形法地应变计算精度的影响［J］.大地测量与地球动力学，2013，33（6）：127－129（LüZhipeng，Wu Jicang，Meng Guojie.Effects of Condition Number on Precision of Ground Strain Calculation for Triangle Method［J］.Journal of Geodesy and Geodynamics，2013，33（6）：127－129）

［4］伍吉仓，邓康伟，陈永奇.三角形形状因子对地壳形变计算精度的影响［J］.大地测量与地球动力学，2003，23（3）：26－30（Wu Jicang，Deng Kangwei，Chen Yongqi.Effects of Triangle Shape Factor on Precision of Crustal Deformation Calculated［J］.Journal of Geodesy and Geodynamics，2003，23（3）：26－30）

［5］邓继华，蔡松柏，钟勇，等.高精度非线性平面应力单元的研究与应用［J］.重庆建筑大学学报，2007，29（3）：83－87（Deng Jihua，Cai Songbai，Zhong Yong，et al.Study and Application of High Precision Nonlinear Plane Stress Element［J］.Journal of Chongqing Jianzhu University，2007，29（3）：83－87）

［6］邓继华，蔡松柏.用四边形平面应力单元进行平面梁的几何非线性分析［J］.长沙理工大学学报：自然科学版，2007，4（2）：32－35（Deng Jihua，Cai Songbai.The Geometrical Non－Linearity Analysis for Plane Beam with Quadrilateral Plane Stress Element［J］.Journal of Changsha University of Science and Technology：Natural Science，2007，4（2）：32－35）

［7］孙聪，李春光，郑宏，等.基于四边形网格的边坡上限有限元法［J］.岩石力学与工程学报，2015，34（1）：114－120（Sun Cong，Li Chunguang，Zheng Hong.Upper Bound Limit Analysis of Slopes Using Quadrilateral Finite Elements［J］.Chinese Journal of Rock Mechanics and Engineering，2015，34（1）：114－120）

［8］石耀霖，朱守彪.用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）

［9］武艳强，江在森，杨国华，等.用球谐函数整体解算GPS应变场方法研究［J］.大地测量与地球动力学，2009，29（6）：68－73（Wu Yanqiang，Jiang Zaisen，Yang Guohua，et al.Research on Method for Entire Calculation of GPS Strain Field by Using Spherical Harmonic Function［J］.Journal of Geodesy and Geodynamics，2009，29（6）：68－73）

［10］万永革.汶川地震产生的应变场和位移场［J］.防灾科技学院学报，2009，11（1）：5－9（Wan Yongge.Displacement and Strain Field Caused by Wenchuan Earthquake［J］.Journal of Institute of Disaster－Prevention Science and Technology，2009，11（1）：5－9）

［11］孙东生，Lin Weiren，崔军文，等.非弹性应变恢复法三维地应力测量——汶川地震科学钻孔中的应用［J］.中国科学：地球科学，2014，44（3）：510－518（Sun Dongsheng，Lin Weiren，Cui Junwen，et al.Three－Dimensional in Situ Stress Determination by Anelastic Strain Recovery and Its Application at the Wenchuan Earthquake Fault Scientific Drilling Hole［J］.Science China Earth Sciences，2014，44（3）：510－518）