青藏高原东北缘形成于晚新生代，是青藏高原外扩和崛起的最新组成部分，也是青藏高原向北东扩展的前缘部位^[1-2]。区域内部构造变形强烈，形变类型丰富，遍布全区的第四纪褶皱、逆冲和走滑断裂表明，整个区域不仅吸收了青藏高原北东向的挤压力，还对调节阿尔金和海原断裂带和昆仑断裂左旋走滑的形成、发展及演化等具有重要的影响^[1-3]。青藏高原东北缘内的强烈的垂直隆升作用，代表青藏高原最新的及正在形成的组成部分在青藏高原地质历史中扮演着重要的边界约束角色^[4-6]，为研究大陆内部变形及其动力学过程提供了天然的实验室。

20世纪90年代GPS开始应用于地球科学领域^[7]，地壳形变与地球动力学中的GPS大地测量学为青藏高原东北缘的构造形变研究带来了更精细的认识。研究结果显示，青藏高原东北缘主要以压缩应变为主(图 1(a))，具体表现为强烈的地壳缩短，同时在岩石圈模拟与地质滑动速率研究中^[8-9]，东昆仑断裂与阿尔金断裂滑动速率存在明显的东向衰减，这并不是块体物质东向挤出^[10]，而是高原内部岩石圈物质的东向流动^[11-13]。GPS垂直速度显示，青藏高原目前的隆升速率大约为2 mm/a，与青藏高原平均的隆升速度大致相同，但与青藏高原中部和南部的相对稳定或下降^[12]有所不同，可能与高地形及重力势能的作用有关。利用GPS速度场得到的应变量/应变率可以用来反映地球表面的震间弹性加载及地震旋回中震前与震后的运动学特征，结合相关的历史地震资料与地质证据可以得到断层和板块间相互作用的关系及作用方式。国内外地球科学家将利用GPS速度场进行地壳应变分析视为一种研究地震旋回与地震危险性评估的重要手段，并不断探索与更新GPS速度场应变分析的研究^[14-25]。另外，有学者针对GPS地壳应变的求解，根据各自的研究方法编写了程序及软件，为地球科学相关领域的研究带来了很大的便利，如QOCA(http://qoca.jpl.nasa.gov/)、VISR^[19]、Geostrain^[23]、Grid_strain and Grid_strain3^[24]、Surfacevel2strain^[22]、SSPX^[20]等。Wu等^[26]认为，在较小区域的研究中使用Delaunay插值可以得到相对准确的应变场，而在大区域的研究中使用最小二乘配置方法比较合适。但这些方法都基于特定的插值函数与插值方法，没有考虑地壳弹性块体的物理特性，且都是在二维空间上考虑GPS反映的平面应变，并没有考虑GPS速度球面三维空间分布，导致得到的结果误差不尽相同。考虑到以上因素，本文利用一种新的三维GPS应变方法研究青藏高原东北缘的构造变形场，以期得到更为精确的地球表面的三维应变场。首先利用区域离散的GPS观测速度场，参照弹性板片的三维GPS速度插值方法得到最佳的薄片弹性参数，模拟得到连续网格化的GPS运动速度场；再对于每个网格作为球面应变的输入，利用最小二乘得到应变张量；最后通过应变张量的分析与垂直运动速度的比较，联合分析青藏高原东北缘的构造形变特征。

1 数据与方法

对青藏高原东北缘的GPS三维速度场^[12]进行研究，首先将收集到的GPS速度场统一转换到相对于欧亚板块下，以更直观地反映印度板块与欧亚板块碰撞对青藏高原东北缘的运动学影响。由于应变场的计算不受参考框架的限制，利用转化后的GPS速度作为输入进行应变分析，在应变计算中假定地球表面弹性板块内部的体力加载致使块体发生了形变，将这些体力约束到观测的GPS点位上，可得到矢量化的形变场，通过调整弹性板片体力的大小让GPS点位的观测速度与计算得到的形变场之间的误差最小，最后得到模拟的连续GPS运动场。在运算模拟中，Haines等^[27]主要利用有限元模拟计算得到GPS形变场；Sandwell等^[28]主要利用计算解析格林函数得到二维GPS形变场，通过调整弹性板块的泊松比使最佳拟合的速度场与GPS速度观测场的误差最小。由于该插值方法考虑了地球表面弹性块体的物理特性，得到的插值结果更加真实可靠，并且弥补了在GPS点位空白区域的计算。本文参考Uieda^[29]的方法，利用格林函数模拟GPS连续运动场，用于解决弹性半空间的Cerruti问题：

$ \left| {\begin{array}{*{20}{c}} {{\boldsymbol{G}_{xx}}}&{{\boldsymbol{G}_{xy}}}&{{\boldsymbol{G}_{xz}}}\\ {{\boldsymbol{G}_{yx}}}&{{\boldsymbol{G}_{yy}}}&{{\boldsymbol{G}_{yz}}}\\ {{\boldsymbol{G}_{zx}}}&{{\boldsymbol{G}_{zy}}}&{{\boldsymbol{G}_{zz}}} \end{array}} \right|\left| {\begin{array}{*{20}{c}} {\boldsymbol{f}_x}\\ {\boldsymbol{f}_y}\\ {\boldsymbol{f}_z} \end{array}} \right| = \left| {\begin{array}{*{20}{c}} {\boldsymbol{u}_x}\\ {\boldsymbol{u}_y}\\ {\boldsymbol{u}_z} \end{array}} \right| $

(1)

$ \left( {{\boldsymbol{\overline{\overline {G}} }^{\rm T}}\boldsymbol{\overline{\overline {W}} }\boldsymbol{\overline{\overline {G}}} + \mu \boldsymbol{\overline{\overline {I}} }} \right)\hat f = \boldsymbol{\overline{\overline {G}} ^{\rm T}}\boldsymbol{\overline{\overline {W}}} \bar u $

(2)

式中，μ为泊松比，G为格林函数，f为区域内部的体力，W为权重，u为GPS的位移场。参考全球大陆地壳的平均泊松比值，结合青藏高原东北缘的速度模型^[30-31]，选择泊松比0.25作为参考标准值，网格距离为3′，同时利用最小二乘方法实现对于正则化参数、耦合参数及体力的最佳深度搜索。该方法利用交叉验证技巧自动选择最佳值来控制插值的参数，交叉验证通过将数据分为两部分分别拟合和评估插值的性能，尤其是程序算法用K-fold交叉验证，明显优于双调和的Spline函数对于每个方向单独的速度插值^[32]。根据算法得到的插值与误差见图 1(b)和1(c)。再利用模拟得到的连续GPS速度值选用球面求解应变的方法求得每个网格内的应变张量场，其中，挤压应变为正，拉张应变为负。

在得到连续的GPS速度场后，将每个网格的速度矢量作为输入，利用球面应变公式^[33]得到位移/速度的泰勒级数展开式为：

$ {u_\theta } = {U_\theta } + \frac{{\partial {u_\theta }}}{{\partial \varphi }}{\rm{\Delta }}\varphi + \frac{{\partial {u_\theta }}}{{\partial \theta }}{\rm{\Delta }}\theta + {r_o}\frac{{\partial {u_\theta }}}{{\partial r}}\frac{{{\rm{\Delta }}r}}{{{r_o}}} + \cdots $

(3)

$ {u_\varphi } = {U_\varphi } + \frac{{\partial {u_\varphi }}}{{\partial \varphi }}{\rm{\Delta }}\varphi + \frac{{\partial {u_\varphi }}}{{\partial \theta }}{\rm{\Delta }}\theta + {r_o}\frac{{\partial {u_\varphi }}}{{\partial r}}\frac{{{\rm{\Delta }}r}}{{{r_o}}} + \cdots $

(4)

$ {u_r} = {U_r} + \frac{{\partial {u_r}}}{{\partial \varphi }}{\rm{\Delta }}\varphi + \frac{{\partial {u_r}}}{{\partial \theta }}{\rm{\Delta }}\theta + {r_o}\frac{{\partial {u_r}}}{{\partial r}}\frac{{{\rm{\Delta }}r}}{{{r_o}}} + \cdots $

(5)

应变场各个分量的表达式为：

$ \begin{array}{c} {e_{rr}} = \frac{{\partial {u_r}}}{{\partial r}};{e_{\theta \theta }} = \frac{1}{r}\frac{{\partial {u_\theta }}}{{\partial \theta }} + \frac{{{u_r}}}{r}, \\ {e_{\varphi \varphi }} = \frac{1}{{r\sin \theta }}\frac{{\partial {u_\varphi }}}{{\partial \varphi }} + \frac{{{u_\theta }}}{r}\cot \theta + \frac{{{u_r}}}{r} \end{array} $

(6)

$ \begin{array}{c} 2{e_{\theta \varphi }} = \frac{1}{r}\frac{{\partial {u_\varphi }}}{{\partial \theta }} - \frac{{{u_\varphi }}}{r}\cot \theta + \frac{1}{{r\sin \theta }}\frac{{\partial {u_\theta }}}{{\partial \varphi }}, \\ 2{e_{\varphi r}} = \frac{1}{{r\sin \theta }}\frac{{\partial {u_r}}}{{\partial \varphi }} + \frac{{\partial {u_\varphi }}}{{\partial r}} - \frac{{{u_\varphi }}}{r}, \\ 2{e_{\theta r}} = \frac{{\partial {u_\varphi }}}{{\partial r}} - \frac{{{u_\varphi }}}{r} + \frac{1}{r}\frac{{\partial {u_r}}}{{\partial \theta }} \end{array} $

(7)

$ \begin{array}{c} 2{\omega _\theta } = \frac{1}{{r\sin \theta }}\frac{{\partial {u_r}}}{{\partial \varphi }} - \frac{{\partial {u_\varphi }}}{{\partial r}} - \frac{{{u_\varphi }}}{r}, \\ 2{\omega _\varphi } = \frac{{\partial {u_\theta }}}{{\partial r}} + \frac{{{u_\theta }}}{r} - \frac{1}{r}\frac{{\partial {u_r}}}{{\partial \theta }}, \\ 2{\omega _r} = \frac{1}{r}\frac{{\partial {u_\varphi }}}{{\partial \theta }} - \frac{{{u_\varphi }}}{r}\cot \theta - \frac{1}{{r\sin \theta }}\frac{{\partial {u_\theta }}}{{\partial \varphi }} \end{array} $

(8)

将式(6)~(8)代入式(1)~(3)，得到球面位移/速度与应变/应变率的关系式为：

$ \begin{array}{c} {u_\varphi } = {U_\varphi }{U_\theta }\cos {\theta _0}\Delta \varphi - {U_r}\sin {\theta _0}{\rm{\Delta }}\varphi + \\ {e_{\varphi \varphi }}{r_o}\sin {\theta _0}{\rm{\Delta }}\varphi + {r_o}{e_{\theta \varphi }}{\rm{\Delta }}\theta + {e_{\varphi r}}{\rm{\Delta }}r + {r_o}{\omega _r}{\rm{\Delta }}\theta - {\omega _\theta }{\rm{\Delta }}r \end{array} $

(9)

$ \begin{array}{c} {u_\theta } = {U_\theta } + {U_\varphi }\cos {\theta _0}\Delta \varphi - {U_r}{\rm{\Delta }}\theta + {e_{\theta \varphi }}{r_o}\sin {\theta _0}{\rm{\Delta }}\varphi + \\ {r_o}{e_{\theta \theta }}{\rm{\Delta }}\theta + {e_{\theta r}}{\rm{\Delta }}r - {\omega _r}{r_o}\sin {\theta _0}{\rm{\Delta }}\varphi + {\omega _\varphi }{\rm{\Delta }}r \end{array} $

(10)

$ \begin{array}{c} {u_r} = {U_r} + {U_\theta }{\rm{\Delta }}\theta + {U_\varphi }\sin {\theta _0}{\rm{\Delta }}\varphi + {e_{\theta r}}{r_o}\sin {\theta _0}{\rm{\Delta }}\varphi + \\ {r_o}{e_{\theta r}}{\rm{\Delta }}\theta + {e_{rr}}{\rm{\Delta }}r - {r_o}{\omega _\varphi }{\rm{\Delta }}\theta + {\omega _\theta }{r_o}\sin {\theta _0}{\rm{\Delta }}\varphi \end{array} $

(11)

最后根据应变张量，利用最小二乘反演方法^[34]得到主应变、剪切应变、面应变及垂直应变，结果见图 2。

2 构造变形分析与讨论

通过应变计算发现，应变的分布与区域的褶皱、断裂分布是一致的，主应变及走滑剪切基本集中在区域的走滑断裂上(图 2)，如阿尔金断裂、东昆仑断裂、海原断裂，且剪切应变最大值在40~60 nstrain(文中计算的速度为应变率，为方便描述用应变代替，1 nstrain = 10^－9 /a)左右，主应变值在80~100 nstrain左右。祁连山褶皱带主要以NE-SW向的压缩应变为主，且祁连山东边与西边的压缩应变明显不同，西边的压缩应变明显小于东边，说明地壳缩短在东北缘东边地区的挤压明显大于西边，与该区域EW向或NE-SW向的地壳挤压是吻合的。同时，西边区域受阿尔金大型走滑断裂的影响，挤压没有那么明显，而在东边区域，鄂尔多斯块体是一个明显的稳定地块，阻挡了地壳NE-SW向挤出，因此地壳压缩应变非常明显。野外地质调查及遥感影像分析结果与GPS速度场较为一致^[35]，青藏高原喜马拉雅地带的变形主要表现在青藏高原的边界上，主要的大型逆冲断裂和走滑断裂也在这里，而其内部的形变则相对较弱，伴随的是大量正断层与走滑断层，这种形变分布方式不仅适用于青藏高原，且几乎遍布整个亚洲大陆：正断层的走向几乎是NS向分布，逆断层局限在高原的边界。从祁连山区域的震源机制解可以看出，青藏高原东北缘地区主要集中了逆冲型地震与走滑型地震。

柴达木盆地区域GPS点位分布很少，应变结果分布没有明显规律。在共和盆地中，应变并没有那么明显，但是地震的活动相对比较频繁，认为GPS现今的应变仅仅可以描述现今的缩短速率，而这些地震可能是来源于历史地震中没有完全破裂的残余能量的累积，所以GPS可以给出一个详细的地壳形变图像，但在地震危险评估中需要结合古地震及构造地质背景的资料进行分析。从垂直应变图像可以看出，垂直应变最高值主要集中在日月山断裂东部及海原断裂的南边，与该区域明显的隆升速率有关，而在祁连山褶皱逆冲带，相对隆升速率则没有那么高，与平面的应变也是一致的。在断裂上，剪切应变从西往东的衰减可以看出区域明显的剪切走滑与挤压的转换，而鄂尔多斯地块与阿拉善地块的阻挡对于海原断裂南部的地壳增厚具有重要的作用，所以该区域在印度板块向北推挤的过程中，应变跨过东昆仑断裂后在祁连山区域发生明显的地壳缩短，而受到相邻块体的相互作用及大型走滑断裂的走滑剪切的转换作用，该区域相对于周围块体作顺时针旋转运动。同时，地壳缩短的强烈变形分异为区域在海原断裂南边明显的地壳增厚，可能是垂直应变主要集中在日月山断裂以东、海原断裂以南的主要构造形变模式。地震资料显示，在西宁盆地下发现大片的低速带^[36]可能受中下地壳流的控制。这种作用模式也可能导致上地壳的增厚，因此该区域的垂直隆升可能不仅是地壳缩短的结果，还受中下地壳流作用模式的控制。从面应变与垂直应变的分布对比结果得出，扩张应变与垂直应变没有很好地对应，说明在考虑区域地壳变形的情况下，不能仅将区域构造变形当作完全弹性体去解释形变行为(完全弹性体体应变为0)，也就是说，面应变的大小不能直接反映垂直应变的拉伸或收缩。

由此可知，本文利用三维的应变计算方法^[29]，考虑到区域的岩石圈物理性质，将离散的GPS速度插值为连续的GPS速度值，经过误差的对比与相关系数的比较，利用区域的泊松比及地球球面作限定条件，相对准确地反映了岩石圈的真实形变行为，与前人的研究结果基本一致。本文反映的是三维形变特征，为探讨区域的构造形变特征提供了更好的认识，但在GPS站点相对较少的区域，反映的仅为应变平滑的结果，不是区域内真实的应变特征。因此，在GPS应变分析中，精密的原始GPS台站速度场有助于提高应变场的精度，尤其对于一些大型的、伴随蠕滑的断裂，可更加准确精细地反映断裂的运动学特征。

3 结语

本文主要基于弹性薄片的格林函数对GPS三维速度进行插值，不仅考虑了站点的分布情况及地壳的弹性特征，同时考虑了GPS速度场在三维空间的分布情况。通过三维的GPS应变分析得到青藏高原东北缘应变的构造分布特征为：

1) 应变的最大值都集中在断裂处，高达100 nstrain/a左右，与地质构造背景及震源机制解的分布是一致的。区域内大型走滑断裂在调节区域构造变形过程中扮演着重要的角色，应变在大型走滑断裂上的分布呈由西往东衰减的现象，与地质上的滑动速率向东衰减是一致的，也反映了区域剪切走滑与挤压缩短的转换。

2) 从垂直应变可以看出，区域内隆升的最大处不是发生在祁连山山体，而是在日月山断裂以东及海原断裂以南，与祁连山的地壳缩短分异为块体旋转及物质东流有关，而这些区域与低速层的分布情况也有很好的对应，说明该区域的隆升不仅是挤压导致的地壳增厚，也可能是中下地壳流导致的增厚，或是二者共同作用的结果。

3) 由面应变与垂直应变的分布对比结果可知，扩张应变与垂直应变没有很好地对应，说明在考虑区域地壳变形的情况下，不能仅将区域的构造变形当作完全弹性体去解释形变行为(完全弹性体体应变为0)，面应变的大小不能直接地反映垂直应变的拉伸或者收缩。

另外，在计算应变的结果中发现一些结果仍存在偏差，可能是区域站点分布稀疏导致的。因此，更多的GPS站点是区域更加精细运动图像的保证，也有助于研究孕震断层的运动学特征。本文研究对理解青藏高原东北缘的运动学过程与断层块体之间的相互作用及青藏高原在东北缘地区的扩张与隆升具有重要意义。

致谢: 感谢Leonardo分享计算三维GPS应变的程序。