基于最小二乘配置模型探讨青藏高原巴彦喀拉块体的运动特征

张贵钢,杨志强,王庆良,刘同文

(1.长安大学地质工程与测绘学院,陕西西安710054;2.中国地震局第二监测中心,陕西西安710054)

0 引言

巴彦喀拉块体位于青藏高原东缘,北部以东昆仑断裂带为界,南部为鲜水河断裂,东侧以龙门山断裂带为界与四川盆地相接。作为青藏板块中部的次级块体,其南北控制断裂活动都较活跃。其北部东西走向的东昆仑断裂平均滑动速率为12～13 mm/a,晚第四纪活动十分强烈[1],2001年昆仑山发生8.1级地震;南部为鲜水河断裂,走向为北西—南东向,滑动速率约14 mm/a[2],20世纪先后发生炉霍7.2级、道孚6.9级地震,属于地震多发带。作为巴彦喀拉块体东向边界的龙门山断裂带为北东走向,与上述两条断裂带不同,其滑动速率小于3 mm/a[3],但仍发生了汶川8.0级地震,显示出较为强烈的逆冲构造特征[3-4]。

该区域控制断裂活跃的构造特性,说明巴彦喀拉块体内部在青藏高原北东向挤压及物质东移的影响下,板块的整体运动及内部构造活动都比较剧烈。笔者利用最小二乘配置模型,对巴彦喀拉块体的整体运动趋势及内部应变场特征进行分析,以期为进一步揭示青藏板块在四川盆地的阻挡下地壳的形变特征及深部动力机制提供参考。

1 地壳形变模型

1.1 板块运动模型

过去对板块运动的研究,一直都把板块和板内地块视为刚体,在其内部是不产生变形的[5]。刚体板块或地块在球面上的运动可描述为[6]

式中:Ve、Vn为地块上任一点(λ,φ)的东向与北向速率;r为地球的平均半径;wx、wy、wz分别为地块在三维空间直角坐标系下绕x、y、z轴的旋转角速率。

实际上,板块或地块更接近于弹性体或粘弹性体,由于区域应力及地壳介质的差异,其内部形变会呈现不同的变形特征。通常利用板块的刚性旋转及内部的应变场来系统描述板块运动及内部形变特征[7]。

随着GPS技术的发展,大量观测资料保证了块体运动及应变场计算的可行性。许多学者对板块运动特征及其内部形变进行了深入研究[8-9],在研究过程中需要解决两个问题:①如何有效区分监测点运动速率中的块体刚性旋转部分和内部形变部分;②针对监测点分布疏密程度的差异,如何适当对监测点数据进行加权分配,以避免监测点密集区权重较大,造成分析结果的偏差。

1.2 最小二乘配置模型

最小二乘配置也称为拟合推估,由 Trarup于1969年提出,最初是用来研究地球形状与重力场的一种数学方法[10]。其基础方程为式中:L为观测向量;A为系数矩阵;X为倾向参数向量;t为观测点的信号向量;ε为观测值的噪声向量。

将最小二乘配置模型应用于块体运动及内部形变分析,把块体刚性旋转的欧拉矢量视为模型中的倾向参数,观测速度扣除整体刚性旋转后的值则看作信号向量与噪声的叠加。基础方程可转化为[11]

式中:V0为监测点的运动速度;Ω为欧拉矢量;A为系数矩阵(GPS观测站坐标信息);Vs为观测速度中扣除整体刚性旋转后剩余部分中的有效信号部分。

式(3)中欧拉矢量Ω,有效信号部分Vs的解Ω、Vs及自协方差 CΩΩ、CVsVs可表示为[12]

式中:C-为观测速度向量协方差与观测噪声协方差之和;Coo为观测速度向量的协方差;Cnn为观测误差的自协方差;Cuo为推估点与观测点信号向量的协方差矩阵。

Coo、Cuo通常采用高斯型函数 F(d)构建。高斯型函数主要是通过建立距离 d与协方差σ(d)的映射关系,拟合函数 F(d)[12]。但在求取σ(d)时经常会出现负值,这显然违反高斯函数的正定特性。因此,通过距离与协方差拟合的方式求取模型系数 K的方法并不完全适用[12-14]。

笔者通过建立点位自协方差与相距最远点之间互协方差的比较关系,削弱距离较远点的影响权值,从而求取较为合理的模型参数 K,构建高斯型函数。

有效信号部分Vs的自协方差(0)为

式中:q为参与建模点的数量。在研究范围内,假设GPS监测点 Pi与 Pj间距离为所有监测点两间距的最大值(dmax),则这两点的协方差γ(dmax)可以认为远远小于信号向量的自协方差(0);因此,设定γ(dmax)=(0),系数 R远小于1。在确定R的基础上,可以相应得到 K,从而建立高斯型函数,确定点位之间的协方差矩阵。

2 巴彦喀拉块体运动特征

基于“中国地壳运动观测网络”1999—2007年共32个GPS监测点的观测资料,采用在欧亚框架下的形变速率。在确定高斯型函数的过程中,取R=10-7,得 K=0.041;构建协方差矩阵时,由于东向速率与北向速率相关性很小,假设其为相互独立观测量,因此,东向和北向速率的协方差σen=0。

基于最小二乘配置模型,得到巴彦喀拉块体的欧拉矢量为((0.45°±0.14°)/Ma,(-17.5°± 5.93°)N,(-77.91°±12.9°)E),其速度场模型如图1。

图1 巴彦喀拉块体速度场模型Fig.1 Velocity Field Model of Bayanhar Block

板块的整体运动在 GPS监测点观测资料中主要反映了运动量的一致性,而内部的应变趋势则体现板块内部各点的形变差异,这与最小二乘配置模型中的倾向参数和随机参数的设置一致。此外,利用观测点间的距离构建协方差矩阵,可以有效避免观测点分布密度较大区域对板块整体运动计算影响过大的问题。因此,将最小二乘配置模型应用于区域地壳形变场的建立,将板块的整体刚性运动视作观测速度的倾向参数,将内部弹性形变视作观测速度的信号参数,从而有效区分板块运动中的整体旋转及内部变形量。根据点位之间的距离确定高斯型函数,继而确定待估点与观测点之间的协方差矩阵,可以有效避免观测点密集区所占权值过大的问题。

图1表明,通过最小二乘配置模型推估的块体内部运动模型与监测点运动的趋势基本吻合,说明通过该方法建立的速率场完全可信。在建立区域地壳速度场的基础上,首先扣除形变场中整体旋转及倾向参数部分,利用扣除后的信号部分(设定噪声的期望值为零,因此不作考虑),通过微分方法可以求得区域应变场[13-14]。

结果显示,在巴彦喀拉块体西部,速度方向主要为北东向,这是印度板块—欧亚板块碰撞下北北东向主压应力影响的结果。巴彦喀拉块体东侧,北向速率逐渐减小,东向速率逐渐加大,特别在康定附近,速度逐渐变为南东东向。在印度板块向欧亚板块俯冲作用下,巴彦喀拉块体内部主压应变呈北东向;青藏高原内部物质东移使得在靠近四川盆地处,主压应变逐渐变为近东西向,这表明在该区域内青藏高原内部两种作用形式对地壳的变形都产生了重要影响。

值得注意的是,图2、3中的应变率场及面膨胀率显示在靠近四川盆地龙门山断裂带附近,无论是应变率还是面膨胀率都较小。其压缩量最大区域处于鲜水河断裂北段,龙门山断裂带前面膨胀率接近于0,这与前期地壳缩短量较小的结论基本吻合。由此可以判断,汶川地震显示出的较强逆冲特性,并非由上地壳挤压逆冲引起,而应该是与下地壳流变的耦合。

图2 巴彦喀拉块体应变率场Fig.2 Strain Field of Bayanhar Block

3 结语

(1)利用最小二乘配置方法建立了巴彦喀拉块体内部的形变场模型,结果表明,最小二乘配置模型通过建立推估点与观测点的协方差矩阵,能够较好地处理点位疏密的定权和块体整体刚性旋转与内部弹性形变的区分这两个问题,建立的形变场模型与观测结果基本吻合。

图3 巴彦喀拉块体面膨胀率Fig.3 Superficial Expansivity of Bayanhar Block

(2)通过对巴彦喀拉块体的速度场、应变率、面膨胀率的特征分析可以看出,作为青藏板块内部的次级块体,该区域同时受到印度板块—欧亚板块碰撞产生的北东向挤压及青藏高原物质东移两方面的作用。应变率与面膨胀结果显示,龙门山断裂前地壳挤压特征并不明显,龙门山断裂带逆冲特性并非通过上地壳脆性挤压缩短来实现,而应该是与下地壳低速高导物质的活动有关。

成文中得到长安大学张永志教授、赵丽华博士以及中国地震局地质研究所武艳强博士的指导,谨致谢忱。

[1] 张培震,王 敏,甘卫军,等.GPS观测的活动断裂滑动速率及其对现今大陆动力作用的制约[J].地学前缘,2003,10 (增):81-92.

[2] 徐锡伟,闻学泽,郑荣章,等.川滇地区活动块体最新构造变动样式及其动力来源[J].中国科学:D辑,2003,33(增):151-162.

[3] 张培震,徐锡伟,闻学泽,等.2008年汶川8.0级地震发震断裂的滑动速率、复发周期和构造成因[J].地球物理学报, 2008,51(4):1066-1073.

[4] 郭祥云,陈学忠,李艳娥.2008年5月12日四川汶川8.0级地震与部分余震的震源机制解[J].地震,2010,30(1):50-60.

[5] 李延兴,杨国华,李 智,等.中国大陆活动地块的运动与应变状态[J].中国科学:D辑,2003,33(增):65-81.

[6] 孙付平,赵 铭.全球五大板块的运动和变形:用卫星激光测距数据导出的站速度分析[C]∥赵 铭,郑大伟.天文地球动力学论文集.上海:上海科技教育出版社,1997:50-58.

[7] Gan WJ,Zhang P Z,Shen Z K,et al.Present-day Crustal Motion Within the Tibetan Plateau Inferred from GPS Measurements[J].Journalof GeophysicalResearch,2007,112 (B08416).DOI:10.1029/2005JB004120.

[8] 魏子卿.现今绝对板块运动[J].地球科学与环境学报,2009, 31(4):331-343.

[9] 刘经南,姚宜斌,施 闯.中国地壳运动整体速度场模型的建立方法研究[J].武汉大学学报:信息科学版,2002,27(4): 331-336.

[10] 崔希璋,於宗俦,陶本藻.广义测量平差[M].武汉:武汉测绘科技大学出版社,2000.

[11] El-Fiky G S,Kato T.Continuous Distribution of the Horizontal Strain in the Tohoku District,Japan,Predicted by Leastsquares Collocation[J].Journal of Geodynamics,1998,27(2): 213-236.

[12] 江在森,刘经南.应用最小二乘配置建立地壳运动速度场与应变场的方法[J].地球物理学报,2010,53(5):1109-1117.

[13] 江在森,马宗晋,张 希,等.GPS初步结果揭示的中国大陆水平应变场与构造变形[J].地球物理学报,2003,46(3):352-358.

[14] 江在森,张 希,陈 兵,等.华北地区近期地壳水平运动与应力应变场特征[J].地球物理学报,2000,43(5):657-665.