基于多面函数法的青藏高原应变特征分析*

吴啸龙 杨志强 武继峰 胡 洋

(长安大学地质工程与测绘学院，西安 710054)

基于多面函数法的青藏高原应变特征分析*

吴啸龙 杨志强 武继峰 胡 洋

(长安大学地质工程与测绘学院，西安 710054)

将曲线拟合法解算大区域低频应变分布与地球参考椭球面上的应变张量矩阵相结合，提出基于多面函数模型解算在给定参考椭球面上的均匀视应变场方法。并根据GPS实测速度场资料解算出青藏高原的均匀视应变场，解算结果与青藏高原各主要构造活动特征以及地震震源机制解所得P、T轴优势方向分布与变化规律吻合较好。认为基于旋转椭球面解算大范围均匀视应变场方法能更加精确地获取大范围内地壳形变的总体变化趋势。关键词 多面函数;参考椭球;应变率;震源机制解;青藏高原

1 引言

WGS84大地坐标系是GPS观测结果所采用的基本参考框架之一，加之参考椭球几何模型更加接近地球的真实形态，因此，利用 GPS观测结果在WGS84定义的参考椭球面上解算地壳应变率，不仅在理论上更加严密，也可以进一步减小模型解算过程中的系统误差。尤其是在大区域地应变解算时，其解算结果将更加逼近真实情况［1］。

青藏高原是世界上海拔最高，地壳最厚、年代最新、面积最广的高原，具有十分复杂的地质结构与构造特征。自20世纪90年代至今，我国先后在青藏高原及其邻域布设的一系列GPS观测站，已经累积了丰富的观测资料。为了更加精确地获取大范围内地壳形变的总体变化趋势，本文利用现有的GPS观测资料，基于多面函数理论模型和WGS84所定义的参考椭球，计算了青藏高原的应变张量，获得了青藏高原整体在参考椭球面上的地壳应变率场分布，并对应变率场分布特征进行了分析。

2 旋转椭球面上的应变张量

应变张量矩阵是一个对称矩阵，在二维情况下为［2－4］:

式(2)就是基于球面的应变张量解。由此可知，应变张量的球面解实际上是椭球面解的一种特殊情况。本文主要研究椭球面上的应变场解算方法，忽略式(1)中的大地高及高程变化，并令uh=0，h=0，可得:

3 基于多面函数理论模型求解应变参数

3.1 多面函数模型［5］

多面函数基本数学模型为:

式中，n为所选建立模型的结点的个数，R(x，y;xj，yj)为核函数。经验上在对地壳运动速率场的拟合中，选取下述函数作为核函数:

式中，(x，y)为推估点的坐标，(xj，yj)为用于建立多面函数模型所选取的测点资料，δ为光滑因子，β为幂数，通过取不同的值来构建不同的核函数，一般取1/2(正双曲型)，3/2(三次曲面型)，－1/2(负双曲型)。

为了考察光滑因子的不同取值与核函数的选取对模型构建的有效性。本文采取交叉验证法对建立的模型进行评价。并分别尝试利用正双曲面、倒双曲面和三次曲面进行建模。选择50、100、145个数据作为结点，其余观测点作为检核点，对高原地区GPS实测速率场进行拟合分析。每次建模时，令光滑因子δ在区间［0，2］内随机取值。并将东向、北向速率各自的拟合精度求平均数，来评价当前模型的拟合精度。通过反复试验，当δ=1.14时，得到的拟合结果最好(表1);进一步选取双三次曲面(β=)作为核函数，多面函数建模效果最好。

表1 多面函数拟合精度Tab.1 Fitting accuracy of multi-surface function

3.2 应变参数解算方法

对GPS速度场的东向分量Ve与北向分量Vn分别建立多面函数模型为:

对式中两项分别求偏微分得:

4 应变场解算结果及分析

分析数据来源于CMONOC工程，数据处理主要分为四步［6］:1)基于 GAMIT软件处理载波相位以获取测站坐标和卫星轨道的单日松驰解;2)利用GLOBK将IGS产出的全球IGS跟踪站单日松驰解作为约束条件，解算得到包含区域点和IGS站的单日松驰解;3)基于ITRF2000框架下的速度解，采用QOCA软件对全部获得的单日松驰解进行平差，估算研究区域内测站坐标和速度;4)基于ITRF2000求解欧亚板块欧拉矢量，并以此为基础进一步解算研究区域各测站相对于稳定欧亚参考框架下的相对运动速率场。

首先选定青藏高原及其邻域75°E ～105°E，23°N～39°N的范围内建立1°×1°的格网节点;然后利用多面函数模型建立研究区域均匀分布的速率场。如图1(蓝色箭头代表拟合多面函数模型拟合结果，黑色箭头代表实测速度场)所示，图中的实测速度场与拟合速度场二者符合甚好，说明利用多面函数法对青藏高原GPS观测速度场进行拟合是一种十分有效的方法。但是利用多面函数法拟合研究区域速率场类似于低通滤波，过滤掉了局部变化或者偏差较大的数据，进而使得到的速率场的整体趋势性特征更加明显。从图1中看出，青藏高原分布的速度场整体上从喜马拉雅碰撞带的主边界中心区域向北东方向偏转，并呈扇形辐射打开。整体上，高原南部的速率场远大于北部的，同时青藏高原的东部边缘区域发生顺时针偏转。这显示出了青藏高原内部形变的两个基本特征:1)高原内部最为显著的特点是印度板块在喜马拉雅冲撞带与欧亚大陆发生碰撞，NNE向主压应力促成了青藏高原整体北东向运移，在青藏高原的北缘受到一系列稳定块体的阻挡，进而形成高原现今南北向的压缩;2)高原南部在侧向挤压和重力拖拽的作用下，塑性的下地壳物质流拖动了上地壳脆性层发生东向逃逸。

基于上述建立的函数模型，按应变参数解算方法求解区域均匀视应变。得到研究区域面膨胀值分布、最大剪应变分布以及最大、最小主应变率场(图2～4)。

图2显示青藏高原南北边界都处于压缩状态，在青藏高原南缘的主边界冲断带附近，地壳被剧烈压缩。这一近NS向的推挤作用通过高原内部各块体的相互作用向北传递，最终在高原北部边界受到塔里木及阿拉善等稳定块体的阻挡形成了南强北弱的地壳压缩。另外在青藏高原中部存在大范围的引张区。这与该区域正断型地震多发以及一系列的南北走向的拉张地堑分布十分吻合。许志勤等［7］研究认为在这一区域下地壳存在大规模的地幔上升流可能是造成高原中部相对较大区域的膨胀作用的深部原因。青藏高原内部存在另一个规模宏大、以川滇块体的边界走滑断裂为界的面膨胀区。结合该地区大量的火山分布，推测其上地壳扩展的原因可能是下地壳软流圈的上涌以及高原物质的大规模东流共同作用结果。

图2 面膨胀Fig.2 Distribution of surface dilatation rate

图3显示出高原内最大剪应变作用最广泛、强烈的区域同样位于青藏高原南部的喜马拉雅主冲断带上。鲜水河断裂与阿萨姆地区的剪应变率次之。最大剪应变分布与高原内部各深大走滑断裂系分布均可以较好吻合。最大剪应变分布能够清晰地显示出阿尔金断裂带、祁连山断裂带以及鲜水河断裂带的基本轮廓，且最大剪应变几乎全部倾向于于青藏高原内部。说明青藏高原相对于其周围稳定地体，具有明显的侧向挤出趋势。如此也验证了基于椭球面解算大区域地壳均匀视应变场分布的方法不仅在理论上更加严密，也可以得到十分精细的结果。

综合来看，面膨胀值的大范围引张地区与最大剪应变分布的峰值区在川滇地块十分吻合，这可能预示着青藏高原的东向挤出在高原东部受到扬子地台的阻挡之后，形成的以“东构造节”为中心的滇藏漩涡。“漩涡构造”的旋转主体正是青藏高原地区剪应变最为显著区域之一的川滇菱形块体。不难看出，川滇菱形块体边界及内部的一系列断裂均显示出明显的剪切走滑特征。而对应嘉黎断裂、理塘断裂的共轭剪切和鲜水河-小江断裂系的左旋走滑，可能是由于旋转体内外圈旋转速度差异造成的。地震波层析资料表明，这一区域分布的大型活动走滑断裂多为岩石圈-地幔剪切带。这似乎也暗示着青藏高原“东向逃逸”实际上是以川滇菱形地块为通道，以大规模软塑性物质流拖动上地壳脆性岩石层的形式发生的。

由图4可看出，青藏高原均匀分布的最大、最小主应变场具有类似于GPS速度场的空间变化趋势:主应变方向自西向东存在顺时针旋转，且以南北向的挤压和东西向的拉张为其主要特征。最大主压应变速率的峰值出现在喜马拉雅主冲撞带附近，其值约(5～9)×10－8/a;其次，高值分布区域主要出现在阿萨姆地区至川滇一带。对最大、最小主应变轴优势方向的特征分布与空间变化进行归纳统计，结果如图5所示。青藏高原内部主应变率方位分布统计特征及变化趋势与徐纪人等［8］根据区域内发生的中、强型地震震源机制解得出的地壳岩石圈应力场P、T轴分布方位具有较好的一致性(表2)。

整体来看，来自印度次大陆的NNE向的推挤作用产生的构造应力为现今青藏高原地区地壳岩石圈形变的主要影响。高原内部主压应变方位集中于近南北向分布:自南向北，青藏高原内主压应变方位由其南端的密集分布，向北过渡逐渐向近NS方向的均匀打开。图6所示为青藏高原内部自西向东的主张应变方位分布。不难看出，青藏高原内部普遍存在近东西向的拉张。张应变优势方位自高原内部过渡到东部边缘时，发生的明显的顺时针偏转。

表2 震源机制解P、T轴主方向与主应变率轴方位对比Tab.2 Contrast between the main direction of P，T-axis of the focal mechanism solutions and the principal strain rate axis

5 结论

本文利用多面函数法建立了青藏高原及其邻近区域地壳水平形变函数模型，得到了1°×1°的格网节点上建立了拟合速率场。通过对比实验，认为当光滑因子δ取1.14，核函数为三次曲面时，建模效果最优。根据此多面函数模型进一步导出了基于旋转椭球面上的应变张量公式，并解算得到了青藏高原及其临区的均匀视应变场。

对青藏高原应变场分析发现，青藏高原内部存在南北向的挤压与东西向的拉伸。最大主压应变的峰值主要集中于喜马拉雅冲断带附近。高原内部东西向的拉张与其间分布的一系列南北向正断层与地堑较好对应。面膨胀值分布特征清晰地显示出高原的南北挤压边界、侧向挤出通道与顺时针旋转等高原形变的主要特征。最大剪应变与高原南北边界及内部一系列大型走滑断层空间分布非常吻合。计算得到的高原总体以及各细部地壳主应变率方位均与地震震源机制解得到的P、T轴方向可以较好吻合。主压应变方向在高原南部相对集中且近似与NS分布，向北逐渐扇形打开且趋于东向偏转。主张应变优势方位自西向东发生明显的顺时针转动，说明本文方法具有较好的精确度与可靠性。

1 赖锡安，等.中国大陆现今地壳运动［M］.北京:地震出版社，2004.(Lai Xian，et al.Present-day crustal movement in China constrained［M］.Beijing:Seismological Press，2004)

2 刘序俨，黄声明，林岩钊.地形旋转张量探讨［J］.大地测量与地球动力学，2010，(5):57－63.(Liu Xuyan，Huang Shengming and Lin Yanzhao.Research on rotation tensor of crustal deformation［J］.Journal of Geodesy and Geodynamics，2010，(5):57 －63)

3 刘序俨，黄声明，梁全强.旋转椭球面上的应变与转动张量表达［J］.地震学报，2007，29(3):240 －249.(Liu Xuyan，Huang Shengming and Liang Quangqiang.Expression of Strain and rotation tensor in geodetic coordinates［J］.Acta Seismologica Sinica，2007，29(3):240 －249)

4 陈健，晁定波.椭球大地测量学［M］.北京:测绘出版社，1989.(Chen Jian and Chao Dingbo.Ellipsoidal geodesy［M］.Beijing:Surveying and Mapping Press，1989)

5 石耀林，朱守彪.用GPS位移资料计算应变方法的讨论［J］.大地测量与地球动力学，2006，(1):1 －8.(Shi Yaoling and Zhu Shoubiao.Discussion on method of calculating strain with GPS displacement data［J］.Journal of Geodesy and Geodynamics，2006，(1):1 －8)

6 王琪，张培震，马宗晋.中国大陆现今构造变形GPS观测数据与速度场［J］.地学前缘，2002，9(2):415－429.(Wang Qi，Zhang Peizhen and Ma Zongjin.Chinese contemporary tectonic deformation of the GPS data and velocity field［J］.Earth Science Frontiers，2002，9(2):415 －429)

7 许志勤，等.青藏高原与大陆动力学——地体拼合、碰撞造山及高原隆升的深部驱动力［J］.中国地质，2006，33(2):221 －238.(Xu Zhiqin，et al.The Qinghai-Tibet Plateau and continental dynamics:a review on terrain tectonics，collisional orogenesis，and progress and mechanisms for the rise of the plateau［J］.Geodesy in China，2006，33(2):221－238)

8 徐纪人，赵志新，石川有三.中国大陆地壳应力场与构造运动区域特征研究［J］.地球物理学报，2008，51(3):770－781.(Xu Jiren，Zhao Zhixin and Ishikawa Yuzo.Regional characteristics of crustal stress field and tectonic motions in and around Chinese mainland［J］.Chinese J Geophys.，2008，51(3):770 －781)

ANALYSIS ON QINGHAI-TIBET PLATEAU STRAIN FIELD BASED ON MULTI-SURFACE FUCTION

Wu Xiaolong，Yang Zhiqiang，Wu Jifeng and Hu Yang
(Academy of Geological Engineering and Surveying，Chang’an University，Xi’an710054)

Combined the method of curve fitting when calculating a lower-frequency strain field in a large area with the strain tensor matrix on the earth reference ellipsoid，we put forward a method of calculating strain parameters on reference ellipsoid by using multi-surface function.According to the GPS data，we get the Qinghai-Tibet plateau strain field.The result is in line with both characters of the major tectonic activity within this area and the main direction of P，T-axis from focal mechanism solutions.

multi-surface function;reference ellipsoid;strain rate;focal mechanism solution;Qinghai-Tibet plateau

P315.72+5

A

1671-5942(2013)04-0017-05

2013-01-12

吴啸龙，男，1988年生，在读博士生，研究方向为地壳形变监测与动态大地测量学.E－mail:xlong_wu@foxmail.com