近年来,由于空间测地技术的飞速发展,利用各种卫星资料得到的高精度地壳运动数据为研究地壳运动提供了前所未有的广阔空间。其中,最有效、最直接、最方便的GPS观测技术获得了广泛的应用。作为地壳运动的典型区域之一,青藏高原受到许多研究部门的重视。研究人员布设了大量的GPS站点,并进行了多个时段的观测,获得了大量高精度地壳水平运动资料。许多专家、学者利用这些资料,借助各种数学方法,做了大量的研究(王庆良等,2002; 张培震等,2003; 甘卫军等,2004; 李延兴等,2004; 张晓亮等,2004; 江在森等,2006; 蒋锋云等,2007)。在众多的数值模拟方法中,有限元方法能够细化模型,研究局部的块体运动变形特征,更可能解决用解析方法难以处理的几何问题及材料非均匀性问题。本文中,我们结合GPS大地水平观测资料来研究和探讨用二维有限元技术提取地震前兆的相关方法和问题。

在利用GPS速度场,借助有限元技术研究地壳水平运动时,边界条件往往很难确定。为使建立的有限元模型在数学上保持稳定,通过扣除区域块体的刚性运动,得到反映区域内部相对运动的GPS速度场,然后用这个速度场研究区域地壳水平运动,这样模型就不会产生平移和旋转,我们就能够相对方便地确定模型边界约束。

基于力学模式的大地测量反演是在选择力学模式和一定的边界条件下,利用有限元等数值计算方法,以位移等观测值与相应点计算值之差的平方和趋于最小为准则,来反求待定参数的方法(独知行,卢秀山,2003)。这个反演系统主要包含两个部分:① 力学模型求解部分。从力学角度看,力学模型求解问题实际是边值问题。大地测量反演问题所考虑的力学模型多为非线性的,预测数据与模型参数具有非线性关系。② 目标函数拟合部分。目标函数采用最小二乘拟合,属于非线性优化问题。这些特点决定了基于力学模式的大地测量反演属于非线性反演问题,同时也决定了它的研究内容和方法的特殊性。

笔者结合以上两个方面,利用高精度GPS地壳水平运动速度场成果,通过消除区域刚性运动获得反映青藏高原内部相对运动的速度场。在理想假设的基础上,借助二维有限元技术构建该区平面二维线弹性有限元模型。根据最优化原理,将消除刚性运动后的GPS速度场作为内部约束条件,反演青藏块体边界分段位移,重构统一、均匀的内部变形场,揭示其与构造及强震孕育之间的关系。

本文所用资料为由赖锡安等(2004)主编的《中国大陆现今地壳运动》一书中发表的1998～2003年中国大陆20°～41°N、74°～107°E范围内所有GPS测站速度场(图1)。

如图1所示,我们共选取了青藏高原及其周围471个GPS站点的速度资料。这些站点包括塔里木地块的站点,阿拉善南部的一些站点,鄂尔多斯西缘的部分站点,华南地块西缘的部分站点,川滇地块外围云南西南部的一些站点。可以看出:① GPS站点在空间分布上具有明显的不均匀性,绝大部分都分布在南北地震带和祁连山地块上,其他地方站点则相对稀少,但站点分布基本上还是能反映整个青藏高原的整体运动形态。② 90E°经线以东的GPS站点以向北运动为主,以西的站点向北东—北东东—东西—东东南方向做顺时针旋转运动。③ 站点运动速率由南向北、从西向东呈明显的减小趋势。

图1 青藏高原及周围相对欧亚板块GPS速度场(1998～2003年)
Fig.1 GPS velocity field of the Qinghai-Tibet plateau relative to Eurasia plate(1998-2003)

在分析GPS区域水平位移资料时,参考框架的选取理论上不影响对区域内部水平差异运动的分析结果,但在实际应用中会在一定程度上掩盖和模糊差异构造变形的许多细节,从而影响人们对问题的直观分析(甘卫军等,2004)。因此,选取适当的参考框架,消除区域整体平移和旋转,有利于突出区域相对运动,反映区域内部块体相对变形特征。因此我们将研究区域看作刚体,利用欧拉定理求出块体在球面上的欧拉旋转矢量,然后用相对于欧亚板块的GPS速度场减去利用欧拉矢量反算得到的整个区域的块体刚性水平运动,从而得到反映块体内部相对变形的地壳水平运动速度场。

用测站地心直角坐标系计算块体欧拉旋转矢量的公式为

[V_n

V_e]=[Rsinλ -Rcosλ 0

-Rsinφ -Rsinφsinλ Rcosφ][ω_x

ω_y

ω_z].(1)

式中,V_n,V_e为测站南北向和东西向速率,R为地球半径,λ为经度,φ为纬度。ω_x、ω_y、ω_z为块体欧拉旋转矢量。知道块体欧拉旋转矢量,就可以求出块体上任意一点的刚性运动。用GPS观测得到的速度场减去刚性运动就得到反映块体内部相对运动的速度场(图2)。

图2 扣除刚性整体运动后青藏高原及周围GPS速度场(1998～2003年)
Fig.2 Relative horizontal movement GPS velocity field inside the Qinghai-Tibet block(1998-2003)

(1)青藏地块应变率大小和方向存在明显的分区性:喜马拉雅带及附近地区主应变率量值高; 以90°E经线为界,以西的主压应变率方向基本上为SN、NNE向,以东的主压应变率方向为NNE—EW—EES,呈顺时针旋转; 青藏块体中部以北地区,主压应变方向与主干断裂方向的夹角在30°～60°之间,为左旋剪切变形区。从羌塘地块东部到川滇地块,主压应变率方向是随着主干断裂方向逐渐偏转的,并且与主干断裂的夹角在30°～60°之间,拉张变形明显,为走滑变形区; 印度板块的东、西两个“触角”出现较大量值的应变率,西部的帕米尔“触角”主要以NEE向压应变为主,东部的阿萨姆“触角”由于物质顺时针流出,出现较大的张应变。

(2)整个青藏高原的内部变形以南北向的地壳缩短为主。

(3)研究区M_S≥7.0地震大都发生在相对速度方向变化明显或速率量值衰减较快的部位。

(4)与区域主干断裂的构造活动背景相一致的最大剪应变率的高值区是强震比较容易发生的部位。

(5)与区域主干断裂的构造活动背景相一致的面积膨胀率梯度变化剧烈部位也是强震危险区。

在用GPS水平速度场成果研究地壳水平运动时,由于测点分布不均匀且稀疏,存在着许多困难。利用GPS高精度水平速度场成果,借助二维有限元技术,通过理想的假设,构建区域二维有限元模型,结合最优化原理来反演区域边界位移,重建均匀的、统一的变形场,来研究区域地壳运动变形与强震之间的关系是一条行之有效的途径。