随着科技的发展，GPS技术在地壳形变监测中得到广泛应用，如美国为了监测南加州圣安德烈斯断层带的运动和变形，布设了超过250个站的连续GPS监测网(Hudnut et al., 2001).日本在境内布设了1000多个GPS连续观测站点，平均30 km就有1个GPS点(Kobayashi et al., 1996).我们借助于GPS测地技术可以获得地球表面位移场随时间的变化信息，从而可以研究地球内部应力场的变化引起地表变形的响应分布及构造应力作用(孙文科，1989).

在早期的地应变研究中，由于研究的空间尺度较小，通常是将局部球面运动投影到高斯平面，然后在平面上利用位移(运动)与应变的关系求解研究区的应变场.地壳水平应变场计算模型或计算式可归为四类.一类是在已知平面运动解析式的情况下通过微分方程求解应变，第二类是根据离散点的平面运动通过建立应变与运动的关系方程求解应变，第三类是利用离散边长的相对变化量求解两个正交方向上的正应变和它们之间的剪应变，第四类是基于Delaunay的三角单元法计算应变(张永志等，2007).该方法在假定地壳是弹性体，其应变张量在局部是个常数的条件下，通过小区域的任一个三角形的三个顶点的位移或边长的变化量可以解算应变张量的各分量.江在森运用GPS资料分析了水平应变场空间分布特征与强震的关系(江在森等，2003)，黄建平建立刚体模型与弹性体模型研究了断层或断裂带附近应变的方法(黄建平等，2010)，另外还有武艳强基于最小二乘配置法解算GPS应变场(武艳强等，2009)，由于观测数据含有误差，必然会传播到应变的计算量上，而且计算结果与三角形的形状因子有关，利用该方法要求三角形的形状因子大于0.1 (伍吉仓等，2003).白玉柱应用Okada及Steketee位错模型及目前汶川地震研究成果，反演龙门山近断裂区域在汶川地震中从震源深处至地表的变形(应变)场分布(白玉柱等，2011).

随着GNSS的引入，研究区的空间尺度得以扩大，在平面上的研究已显得不适宜，因此球面应变计算模型相继出现，石耀霖研究了在球面坐标系计算应变的方法(石耀霖等，2006).三角形法在研究区划分不同的三角形，研究者就可能给出不同的结果(王敏等，2003；刘启元等，2003).为了克服三角形计算应变的随意性和结果差异性，EI-Fiky提出了一种基于指数函数的最小二乘拟合推估的方法，通过GPS观测得到的位移计算了日本、意大利、中国华北地区的应变分布，其结果与地质地球物理调查结果符合较好(El-Fiky et al., 1998; 陈永奇等，2007).本文利用计算点与观测点之间的距离，构建协方差与自协方差矩阵，推导了最小二乘预估法计算地壳水平应变的公式，其应变值是与计算点位置有关的变量，而不是一个常量.最小二乘预估法是承袭了随机过程的预估理论，该理论避开了泛函分析采用线性代数理论，计算简单、容易实现为更好地发挥区域GPS网监测在地壳形变和地震预报中的作用，做了一点新的尝试工作.

1 理论与计算方法

假设经过处理的GPS位移观测数据的残余量在空间上是均匀各项同性，其数学期望为零， 位移观测值 l₁,l₂,l₃…l_q组成q维向量 L =[l₁,l₂,l₃…l_q]′，信号s₁,s₂,s₃…s_m组成m维向量 S =[s₁,s₂,s₃…s_m]′，通常上标′表示转置.假设观测量的协方差矩阵只与观测点间距离有关，则可以用观测数据的数据本身构建自协方差 C _ll与互协方差 C _sl，所有的协方差是根据同一个协方差函数C(d)求得的，假设C(d)仅与所讨论的两点间的水平距离有关，则依据 L所得到的S 的最小方差无偏估值为

借助(1)从而可以根据x,y方向上的观测位移值u,v，计算内插点P的位移分量的估计值U _P,V _P 则有

式中C_pi=1/d_pi，d_pi为p和i之间的距离，C_ij=1/d_ij，d_ij为i和j之间的距离，按下式计算P点的应变张量

式中ε_xx,ε_yy,ε_xy是应变的三个分量，可以从(6)式看出，应变分量是观测位移的线性组合，其误差估计可由误差传播定律求得.各向同性介质质中的区域应变能可根据文献(张永志等，2007)求得.

2 渭河盆地的GPS资料及应变场分布

渭河盆地北接鄂尔多斯台地，南邻秦岭皱褶带，是中国重要的构造形变和地震活动区域之一，区内发育一系列活动断裂，地质灾害频繁，为了研究该地区活动断裂的动力学机制，我们利用GAMIT软件对2004-2007年渭河盆地GPS测站的观测资料进行统一解算、处理,获得地表水平运动的统一速率场，选取研究区分布较均匀的30个GPS站点速率，描述其水平运动状况如图1所示.图中黑色箭头为GPS站点水平运动速率分布，箭头大小表示GPS位移速率的大小，箭头指向为位移场速率运动方向.渭河盆地地壳运动在2004-2007年期间整体表现出一种向东向南运动特征，且东、西部地壳水平运动存在一定的差异运动，各站点运动速率为0～6 mm／a.黑色线调为该地区主要断裂带，大部分断层呈东西向水平分布，在盆地东部断层分布随地形向北偏的趋势.

图 1

Fig. 1

图 1 渭盆地GPS位移场及主要断裂分布(2004-2007) Fig.1 The distribution of GPS velocity and faults in Weihe Basin

本文结合2004-2007年的GPS观测资料来对该地区地壳的应变场进行计算，首先选取距计算点最近的3个GPS观测点组成(2)、(3)式，然后利用公式(4)、(5)、(6)计算最大主应变与最小主应变结果如图2，由图2来看主应变呈现出空间分布不均匀性，山区和盆地应变性质相反，山区处于挤压状态，盆地处于拉张状态，这一结论与文献(张永志等，2011)一致.渭河盆地西部主要表现以压应变为主，方向为南-北向、东-西向，南北向压应变大于东西向的压应变，研究区西南部，也是压应变为主但强度相对于西部减弱，方向为西南-东北向，盆地内中部主要东西拉张应变为主，与区域内的断层延伸走向一致，研究区东部地区在水平拉张作用下易引起正断层活动，从而形成地裂缝等地质灾害态势，这与多年来的地质调查结果相一致.

图 2

Fig. 2

图 2 渭河盆地2004-2007主应变矢量图(单位：1×10-9) Fig.2 Diagram of the principal strain vectors of Weihe basin from 2004 to 2007(unit:1×10-9)

图3为2004-2007年最大剪应变率等值线图，根据地震资料记载大部分6级以上地震，其震中位于由震前GPS观测资料解算出的剪应变率的高值区或其边缘，尤其是与区域主干断裂的构造活动背景相一致的剪应变率高值区内(张希等，2004；祝意青等，2013).宝鸡西部位于最大剪应变的高梯度区，等值线密集且应变梯度大，压应变起主导作用，宝鸡西部地区的南与北两侧的GPS也显示出了明显的差异性.在2004-2007年间，该地区又呈现出了应变特征值的高值区，具备引发中强地震的条件.图2、图3应变分布特征可归结为渭河盆地处于构造过渡带，同时受青藏块体顺时针旋转运动和鄂尔多斯块体南缘逆时针旋转运动的共同影响，使得渭河盆地总趋势为向北向东拽拖，两侧的地块产生平移和扭转的原因所致.

图 3

Fig. 3

图3 渭河盆地2004-2007年最大剪应变率等值线图(单位：1×10-9) Fig.3 Isolines of the maximum shear strain rate of Weihe basin from 2004 to 2007(unit:1×10-9)

3 结 论

本文采用的最小二乘预估法计算地壳水平应变的方法，是把观测数据的位移看作是随机信号，利用计算点与观测点之间的距离，构建协方差与自协方差矩阵，从而通过求导计算区域应变场的分布.相对于三角形法在单元三角形内的应变是一常量来说，该方法的优势是计算单元在三角形内的应变是连续函数的变量，其值与该点的位置有关，为 GPS监测地壳形变提供一种新的思路.结合2004-2007年的渭河盆地GPS观测资料，渭河盆地主应变呈现出空间分布不均匀性，山区和盆地应变性质相反，山区处于挤压状态，盆地处于拉张状态.宝鸡西部位于最大剪应变的高梯度区，挤压应变起主导作用.

致 谢 感谢两位评审专家提出的宝贵建议！