采用三維空間數(shù)據(jù)計算汶川地區(qū)地殼應(yīng)變

劉明學(xué) 段虎榮

1 陜西鐵路工程職業(yè)技術(shù)學(xué)院，渭南市站北街東1號，714000

2 西安科技大學(xué)測繪科學(xué)與技術(shù)學(xué)院，西安市雁塔路58號，710054

傳統(tǒng)地殼應(yīng)變計算方法是三角形法［1－3］，該方法主要考慮X、Y兩個方向的主應(yīng)變和XY平面的剪應(yīng)力，其計算結(jié)果反映區(qū)域平面應(yīng)變趨勢［4］。也有學(xué)者對四邊形單元法開展相關(guān)研究［5－7］。石耀霖對球面坐標系計算應(yīng)變進行研究［8］。武艷強采用球諧函數(shù)整體解算了GPS應(yīng)變場［9］。但以上都是研究平面的力學(xué)問題。本文在前人的基礎(chǔ)上，利用2005～2008的GPS觀測資料，采用四邊形計算法對汶川地區(qū)地殼運動和區(qū)域應(yīng)變場特征進行分析。應(yīng)變計算過程中引入高程分量，提出四邊形網(wǎng)絡(luò)構(gòu)造，推導(dǎo)X、Y、Z3 個方向及XY、XZ、YZ3個面的地殼應(yīng)變計算公式，分析了汶川地區(qū)應(yīng)變特征，不僅能反映該區(qū)域平面的應(yīng)變結(jié)果，還能反映垂直應(yīng)變變化趨勢。

圖1 四邊形單元示意圖Fig.1 Schematic diagram of quadrilateral element

1 應(yīng)變計算方法

假定兩個GPS觀測點之間是線性變化的，根據(jù)觀測網(wǎng)建立四邊形圖形網(wǎng)（空間直角坐標系），通過同一個四邊形不同方向上的變化計算出各邊長方向的線應(yīng)變，再用不同方向的線應(yīng)變計算X、Y、Z3個坐標軸方向的主應(yīng)變和XY、XZ、YZ3個面的剪應(yīng)變，最后得出應(yīng)變分布的計算公式。在引入高程的基礎(chǔ)上，對四邊形進行劃分（圖1），利用位移場直接計算應(yīng)變。

圖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）。該四邊形內(nèi)任意一點（x，y，z）的位移場（u，v，w）是坐標的線性函數(shù)，可表示為：

將四邊形的4個角點x軸方向的位移和坐標值代入式（1）中第1式，得：

通過解式（2），可求出系數(shù)a1、a2、a3、a4，同理可解得系數(shù)b1、b2、b3、b4、c1、c2、c3、c4。

主應(yīng)變計算公式為：

剪應(yīng)變計算公式為：

2 觀測數(shù)據(jù)

國家測繪地理信息局、國家重點基礎(chǔ)研究發(fā)展計劃項目“活動地塊邊界帶的動力過程與強震預(yù)測”項目組、國家重大科學(xué)工程“中國地殼運動觀測網(wǎng)絡(luò)”在龍門山斷裂帶兩側(cè)布設(shè)了一系列數(shù)據(jù)監(jiān)測點，采用GPS流動觀測為主，測量了監(jiān)測點的水平位移。中國地震局在2007－04～07進行了多期觀測，觀測周期為3～4d。汶川地震發(fā)生后，又對這些GPS 點進行了復(fù)測，觀測周期為2～3d。國家測繪地理信息局在2005～2008－05－12按等級的不同，采用1～4d的觀測周期進行觀測，震后又對其進行周期為1d的復(fù)測。四川省地震局和重慶市建立以服務(wù)為主的GPS 連續(xù)觀測網(wǎng)絡(luò)，提供了垂向同震位移資料。本文采用雙差模式，數(shù)據(jù)采樣間隔為30s，24h為一時段，用GAMIT／GLOBK 軟件對上述單位的觀測資料進行處理。

3 汶川地區(qū)應(yīng)變特征計算分析

分析研究區(qū)域GPS 觀測點的分布位置（圖2），對觀測點進行四邊形網(wǎng)格化，構(gòu)建出72個四邊形，如圖3所示。

根據(jù)式（3）、（4）解算每個四邊形應(yīng)變參數(shù)，求取X、Y、Z方向的主應(yīng)變和XY、XZ、YZ面的剪應(yīng)變，并作出四邊形對應(yīng)等值線圖（圖4～9）。

圖4中粗黑實線表示斷層（下同）。由圖4可看出，X方向主應(yīng)變值處于［－0.3，0.75］，以龍門山斷裂帶為界，龍門山斷裂帶以西為負應(yīng)變，以東為正應(yīng)變，在成都以南和映秀以西附近分別出現(xiàn)應(yīng)變峰值，映秀、汶川、北川均處于應(yīng)變高梯度帶區(qū)域，說明映秀以西附近地區(qū)在地震中有向東移動的現(xiàn)象。該結(jié)果與實際位移圖中映秀以西區(qū)域向東移動，成都、都江堰區(qū)域向西北移動相對應(yīng)。

圖5顯示青川及附近區(qū)域負應(yīng)變比較大，成都附近和黑水西北方向分別出現(xiàn)正應(yīng)變峰值。對比實際位移圖像，青川在Y軸沿負方向移動量相當大，而映秀附近則和斷裂帶上該區(qū)域大多數(shù)移動方向一致，向Y軸的正方向移動。汶川、北川均處于應(yīng)變高梯度帶區(qū)域。映秀、汶川、北川處于擠壓狀態(tài)，Y方向主應(yīng)變值處于［－28，12］，遠大于X方向。

圖2 汶川地區(qū)斷層及GPS觀測點分布圖Fig.2 The distribution of faults and GPS

圖3 研究區(qū)域四邊形網(wǎng)格劃分Fig.3 The quadrilateral mesh of regional

圖4 X 方向主應(yīng)變等值線圖Fig.4 The contour map of principal strain in Xdirection

圖5 Y 方向主應(yīng)變等值線圖Fig.5 The contour map of principal strain in Ydirection

從圖6可以看出，Z方向的主應(yīng)變值處于［－0.002，0.002 8］，以龍門山斷裂帶為界，龍門山斷裂帶以西為負應(yīng)變，以東為正應(yīng)變，在成都以南和茂縣以西附近分別出現(xiàn)應(yīng)變峰值，映秀、汶川、北川均處于應(yīng)變高梯度帶區(qū)域。

圖6 Z 方向主應(yīng)變等值線圖Fig.6 The contour map of principal strain in Zdirection

圖7中XY面剪應(yīng)變的值處于［－20，12］，高值區(qū)出現(xiàn)在映秀－北川－青川斷裂帶附近，表明這片區(qū)域的剪應(yīng)變值偏大，處于擠壓狀態(tài)。該狀態(tài)與震后的地質(zhì)結(jié)果一致，這一帶發(fā)生右旋剪切應(yīng)變，全長超過300km，地表破裂帶沿著映秀－北川走向斷斷續(xù)續(xù)分布，破裂帶切割了多種地貌單元，包括河流、山脈、沖洪積扇、橋梁和公路等，使道路發(fā)生曲拱，橋梁垮塌、位移。

圖7 XY 面剪應(yīng)變等值線圖Fig.7 The contour map of shear strain in XYdirection

由圖8可以看出，XZ面剪應(yīng)變值處于［－0.3，0.45］，負高值區(qū)出現(xiàn)在北川附近，小金、寶興附近出現(xiàn)較大正剪應(yīng)變，汶川、映秀處于梯度帶區(qū)域。

圖9中YZ面的剪應(yīng)變值為［－5，9］，可以看出，平武、梓潼附近出現(xiàn)高正剪應(yīng)變，寶興、邛崍、成都一帶出現(xiàn)負剪應(yīng)變，北川、汶川、映秀處于剪應(yīng)變梯度帶區(qū)域。

圖8 XZ 面剪應(yīng)變等值線圖Fig.8 The contour map of shear strain in XZdirection

圖9 YZ 面剪應(yīng)變等值線圖Fig.9 The contour map of shear strain in YZ direction

從圖4～6可以看出，Y方向主應(yīng)變最大，X方向次之，Z方向最小，北川、汶川均處于主應(yīng)變的高梯度帶區(qū)域。從圖7～9可以看出，XY面剪應(yīng)變最大，YZ面次之，XZ面最小，汶川、映秀均處于剪應(yīng)變梯度帶區(qū)域。汶川、北川地區(qū)處于水平擠壓應(yīng)變狀態(tài)，西北地區(qū)垂直分量處于向上拉張狀態(tài)，東南地區(qū)垂直分量處于向下擠壓狀態(tài)，與文獻［10－11］結(jié)果一致。

4 結(jié) 語

本文基于汶川震前的GPS定期觀測數(shù)據(jù)與震后第一時間的復(fù)測數(shù)據(jù)，計算汶川地區(qū)地殼應(yīng)變問題，分析汶川地區(qū)地殼變形的基本特征，得到如下結(jié)論：1）平面應(yīng)變結(jié)果顯示汶川、映秀地區(qū)處于東西方向擠壓狀態(tài)，高程結(jié)果顯示在主斷裂帶東南方向應(yīng)變?yōu)檎担鞅狈较驊?yīng)變?yōu)樨撝担礀|南地區(qū)處于向下擠壓狀態(tài)，西北地區(qū)處于向上拉張狀態(tài)；2）Y方向主應(yīng)變最大，X方向次之，Z方向最小，北川、汶川均處于主應(yīng)變的高梯度帶區(qū)域；3）XY面剪應(yīng)變最大，YZ面次之，XZ面 最小，汶川、映秀均處于剪應(yīng)變梯度帶區(qū)域，汶川地震恰好發(fā)生在應(yīng)變高梯度帶。

致謝：感謝國家測繪地理信息局、國家重點基礎(chǔ)研究發(fā)展計劃項目“活動地塊邊界帶的動力過程與強震預(yù)測”項目組為本文研究提供觀測數(shù)據(jù)。

［1］孟國杰，申旭輝，伍吉倉，等.GPS揭示的貝加爾湖地區(qū)現(xiàn)今地殼形變特征［J］.大地測量與地球動力學(xué)，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］肖根如，甘衛(wèi)軍，陳為濤.地應(yīng)變計算Delaunay三角網(wǎng)在MATLAB與GMT 環(huán)境下的相互轉(zhuǎn)換［J］.大地測量與地球動力學(xué)，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］呂志鵬，伍吉倉，孟國杰.條件數(shù)對三角形法地應(yīng)變計算精度的影響［J］.大地測量與地球動力學(xué)，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］.大地測量與地球動力學(xué)，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］鄧繼華，蔡松柏，鐘勇，等.高精度非線性平面應(yīng)力單元的研究與應(yīng)用［J］.重慶建筑大學(xué)學(xué)報，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］鄧繼華，蔡松柏.用四邊形平面應(yīng)力單元進行平面梁的幾何非線性分析［J］.長沙理工大學(xué)學(xué)報：自然科學(xué)版，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］孫聰，李春光，鄭宏，等.基于四邊形網(wǎng)格的邊坡上限有限元法［J］.巖石力學(xué)與工程學(xué)報，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 位移資料計算應(yīng)變方法的討論［J］.大地測量與地球動力學(xué)，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］武艷強，江在森，楊國華，等.用球諧函數(shù)整體解算GPS應(yīng)變場方法研究［J］.大地測量與地球動力學(xué)，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］萬永革.汶川地震產(chǎn)生的應(yīng)變場和位移場［J］.防災(zāi)科技學(xué)院學(xué)報，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，崔軍文，等.非彈性應(yīng)變恢復(fù)法三維地應(yīng)力測量——汶川地震科學(xué)鉆孔中的應(yīng)用［J］.中國科學(xué)：地球科學(xué)，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）