基于最小二乘配置的板块边界地壳形变应变特征分析

王 栋，赵旭坤，刘望明

（陕西铁路工程职业技术学院，陕西 渭南 714000）

关键字：GNSS；最小二乘配置；地壳应变；动力学机制

0 前言

地震、地裂缝及其他地质灾害都伴随着地壳形变的发生，利用有效的监测技术对大规模地质灾害易发区域进行连续高精度监测，可以有效地发现区域的地壳形变特征，结合地质情况总结地质灾害孕育的成因机制对人类生活工作安全具有重要的意义[1]。地表在地球内部应力的作用下发生大规模不均匀的微小变形，这种长期不均匀的运动导致地质灾害的缓慢孕育，GNSS卫星定位技术可以监测这种缓慢的运动现象，能从宏观进行定性分析也能从微观进行定量分析，进而获得地质灾害前兆信息，服务于地震、地裂缝等地质灾害的理论研究[2]。

太平洋板块与北美板块交界地带断裂带纵横，地壳形变剧烈，且多年频繁发生地震，PBO（板块边界观测计划）观测网络拥有密集的GNSS连续运行站，方便获得北美板块西边界的海量GNSS监测数据，利用GAMIT/GLOBK等软件可以获取该区域的地壳速度场，通过速度场获取应变场可以明确区域形变特征，再结合历史地质灾害信息可以预测评估区域的整体地震地质灾害危险性。

大规模GNSS监测网络被广泛应用到地球动力学、地震学、地球物理学等多个地学领域[3]。基于“中国地壳运动观测网络”等GNSS连续运行站网络，我国的学者已获取海量GNSS监测数据，研究了多个区域的地壳形变特征，得出了相应的结论。崔笃信等详细研究了鄂尔多斯周边的形变及应变特征，并与断裂带运动特征进行了结合[4]。高立新等利用多年的GNSS数据研究了鄂尔多斯周边的地壳形变及断裂带的运动特征[5]。唐红涛等利用GNSS速度场剖面进行研究获取了云南部分断裂带为扩张、压缩、扩张特性[6]。徐声鑫分析了云南地区GPS面膨胀与地震间的分布关系[7]。

国外学者研究地壳形变及应变特征也有非常多的成果，世界各区域都建立相应的GNSS形变监测网络，数据量巨大，网络运行稳定。Mohanmmad等详细讨论了最小二乘配置方法的特点及算法，基于模拟数据验证了最小二乘配置解算应变场的合理性[8]。Thomas A.Herring等讨论了PBO网络数据处理的方法[9]。A.V.Vilayev等解算了哈萨克斯坦东南部地震区2009—2014年的地壳运动特征并细致地分析了变形过程[10]。Y.Hamiel等通过GNSS数据分析了河谷周围的地壳形变信息[11]。

1 研究方法

1.1 区域概况及数据来源

北美板块西边界布设了著名的PBO网络，在该区域建立了数千个GNSS连续监测站，并附有其他形变监测设备，如钻孔应变计等。PBO网络提供了详细的GNSS时间序列信息，并做到数据公开，全球大量学者均可进行地壳形变及地震研究。本文利 用Carl Tape等《Multiscale estimation of GPS velocity fields》文章中的PBO速度场资料[12]。其速度场资料来自于NASA Earth Science Research，Education and Applications Solution Network（REASoN）项目，其项目选择了PBO网络中1997—2008年间长达11年的314个连续观测站的速度场数据，其速度场的参考框架为ITRF2005，区域速度场整体朝向NW运动，GNSS西向速度区间为9.1～44.1 mm/yr，北向速度区间为-10.7～25.5 mm/yr，其速度场精度优于1 mm/yr，见图1。

图1 北美板块边界GNSS速度场图Fig.1 GNSSvelocity field diagramof the North American plate boundary

1.2 球面最小二乘配置

最小二乘配置的主体思想是将已知数据点的观测数据分解为整体趋势项、特征信号量及不可避免的随机误差项。通过数据集的统计性质建立协方差函数确定数据集间点与点的相互关系，进而利用已知数据点推求未知数据点。最小二乘配置算法为[13-15]：

式（1）中：L是观测向量，AX是整体趋势性向量，S是已测点的信号向量，S′是未测点的信号向量，S′与观测量L无关，Δ是观测量的误差向量。

用DΔ,QΔΔ,PΔ分别表示误差Δ的方差阵、协因数阵及其权阵。用DY,QYY,PY分别表示随机信号Y的方差阵，协因数阵及其权阵。观测误差和随机信号相互独立，即DΔY=0。取单位权方差σ20=1，则有：

基于广义平差：

可得GNSS观测站的信号值与未进行观测的位置处的信号值估值为：

式（4）中：DS是各GNSS观测站速度信号的方差阵，与经速度值拟合所得的协方差函数有关。DSS′是GNSS观测站的速度信号和未知位置处的速度信号之间的协方差，将通过拟合所得的协方差计算所得。

应变场是通过速度场进行偏导计算所获得的，所以整理速度场和应变场的关系为：

式（5）中：ελ,εφ,ελφ为应变分量，λ,φ为经纬度，R为地球半径，u为形变速度，h为大地高，利用应变分量求解应变特征参数：

式（6）中：ε1为最大主应变，ε2为最小主应变，εarea为面膨胀，rmax为最大剪应变。

2 计算结果及分析

2.1 GNSS相对速度场特征分析

在利用GNSS速度场信息分析区域地壳形变特征时，如果速度场在ITRF框架中，则可以分析区域速度场的整体运动趋势，从北美板块边界的GNSS数据可以得出其整体为WN向运动，为了使区域各个部分的形变差异明显得表现出来，则需要通过欧拉矢量将区域板块的整体运动趋势扣除，才能获取其板块内部各部分之间的差异性运动特征[16-17]。

欧拉旋转参数与GNSS站点速度的关系为：

式（7）中：ve、v n分别表示东向及北向速度分量，R表示地球半径，λ、φ表示经纬度。ωx,ωy,ωz表示欧拉旋转参数，令矩阵A为：

基于最小二乘平差原理，可得欧拉矢量为：

式（9）中：为GNSS速度场误差方差阵。

利用研究区域所有GNSS站点速度信息获取板块整体欧拉旋转参数，进而获取整体运动速度场。整体运动趋势一致性较好，与原始实测速度场对比如下：

北美板块西边界ITRF2005框架实测速度场为蓝色，区域整体运动GNSS速度场为红色，见图2。则实测速度场扣除板块整体运动速度场可得区域GNSS站点明显的差异性运动特征。

图2 站点欧拉速度与实测速度对比Fig.2 Comparison of the Euler speed and the measured speed

从区域差异性运动特征判断出北美板块西边界西向速度值最小为-21.1 mm/yr，最大为21.2 mm/yr，北向速度值最小为-22.0 mm/yr，最大为16.0 mm/yr，见图3。在区域断裂带两边的GNSS站点速度值明显有改变，其相对运动明显，表明该地区地壳形变复杂无序。

图3 区域相对速度场Fig.3 Regional relative velocity field

2.2 GNSS协方差函数拟合

要利用最小二乘配置方法获取区域的地壳形变应变特征，比较重要的是先获取区域的形变协方差函数，利用形变协方差函数才能计算最小二乘配置算法中重要的协方差矩阵，区域整体的形变特征受协方差函数拟合的影响较大，是区域形变特征的体现。利用最小二乘配置算法的前提是获取区域内相邻位置的形变关系，此形变关系与距离呈相关性，根据两点间距离变化与形变特征的对应关系基于统计学原理，拟合协方差函数，就可以获得任意点间的相对形变关系，即求得对应的协方差值。在拟合地壳形变特征关系中指数函数的衰减过程比较符合地壳形变特征的变化，所以将指数函数用来描述地壳形变相关性随距离增加而减小的特征。其函数形式为：

式（10）中：d表示不同GNSS测站的球面距离，C(0)表示速度方差，C(d)表示根据协方差函数获取的两点相距dkm时的关系值，即协方差，k为待定系数，测区不同k值不同，是整体数据统计的结果，C(0)和C(d)公式为：

式（11）中：其中：li、lj代表上述GNSS东或北向相对速度值。

协方差函数拟合的过程为：先计算北美板块西边界各个GNSS站点之间的球面距离，将所有的距离值分成若干段进行统计，分段依次为0～30 km，30～60 km，60～90 km......直至所有GNSS站点间的最大距离，计算各个段内各GNSS站点间的协方差值，即表征了各站点间的相关性，将各段内的协方差值运用于协方差函数的拟合中，值得注意的是，协方差值代表了两点间的相关性，地壳形变相关性计算中，当协方差值为负值，则舍去不用。

利用各段协方差值获取了相应的函数拟合结果，协方差函数分别从北向速度和东向速度进行拟合，表达式为：

利用东向和北向速度值进行相关性的协方差拟合，其东向RMSE为4.61 mm2/yr2，其北向RMSE为6.85 mm2/yr2，其表征相关性强弱的相关系数平方分别为0.95及0.94，代表了指数函数适合运用于该区域的协方差拟合中，根据拟合函数发现当距离值接近350 km时，协方差值接近0，认为超过这个距离时两GNSS站点间无相关性，即相关距离值限值为350 km，则依次类推，北向的相关距离接近360 km，见图4。

图4 指数函数东向及北向协方差拟合示意图Fig.4 Fitting diagram of covariance of exponential function in east and north direction

2.3 区域最小二乘配置推估速度场

获取了区域的协方差函数拟合结果，将各个GNSS站点间的距离值代入，便可获取两GNSS站点间的距离所对应的协方差值，即代表了速度相关性，最小二乘配置算法就利用此相关性来推估区域任意点位处的形变速度值。

将研究区域划分为0.5°×0.5°的经纬度格网，格网点处一般未设立GNSS连续观测站点，可作为推估点利用最小二乘配置算法获取其点位处的速度推估值，从而获取整个北美板块的推估速度场，从推估速度场可以看出整个区域速度场较为均衡，但在北美板块边界处，圣安德列斯断裂周边速度场产生了较大的突变，也代表了区域地壳形变的剧烈性，见图5。

图5 推估速度场Fig.5 Estimated velocity field

2.4 最小二乘配置应变场

利用最小二乘配置算法获取了相应区域的速度场信息，这个速度场在数学上是连续的速度场，即每一点推估的速度值是通过连续函数而导出的，则可以利用连续的速度场对GNSS站点的位置求偏导计算，便可以获取整个区域速度场的变化率，即得到整个北美板块西边界的地壳形变应变场信息。

利用推导的球面应变计算方法，计算相应北美板块区域面膨胀系数，从面膨胀的分布图可以看出：整个北美板块区域面膨胀系数分布极不规则。其中，北美板块西边界带的加利福尼亚州，不规则的面膨胀分布说明其内部产生着复杂的地壳形变，圣弗朗西斯科周边面膨胀约-3×10-7/yr，加利福尼亚中部面膨胀约0.5×10-7/yr，加利福尼亚南部面膨胀约-4×10-7/yr。板块边界地壳形变复杂，压缩和膨胀交替出现，圣马利亚北部圣安德列斯断裂周缘压缩率最大约-11×10-7/yr，洛杉矶沉积盆地西北部膨胀率最大约2×10-7/yr，见图6。

图6 北美板块边界面膨胀示意图Fig.6 Schematic diagram of the expansion of theboundary surfaceof the North American plate

剪应变明显分布于北美板块边界的断裂带周缘，剪应变的值可以直接说明区域地壳形变的剧烈程度，在圣弗朗西斯科东南部剪应变量值约6×10-7/yr，代表该区域地壳形变剧烈。圣安德列斯大断裂的最南端剪应变率发生突变，地壳形变复杂，从1×10-7/yr增加至3×10-7/yr，见图7。

图7 北美板块边界剪应变分布图Fig.7 Distribution of shear strain on the North American plate boundary

最大主应变分布也基本沿着北美板块西边界带的断裂带，其量值约2×10-7/yr。贯穿板块边界的圣安德列斯断裂环绕着明显的主应变值，约1×10-7/yr至3×10-7/yr不等，断裂周缘地壳形变剧烈，见图8。

图8 北美板块边界最大主应变示意图Fig.8 Schematic diagramof the maximumprincipal strain on the North American plate boundary

区域应变矢量分布杂乱无序，在加利福尼亚西北部，N-S向主压应变以及E-W向主压应变占据主导地位，不同区域的应变矢量方向及大小有明显的区别，说明区域地壳形变剧烈且没有规律可循。圣弗朗西斯科东南部是应变张量的极值地带。圣安德列斯断裂周缘的SW-NE向约7×10-7/yr的主压应变与其走向垂直，且同时存在约1.5×10-7/yr的主张应变与断层走向呈一定角度。北美板块地壳形变剧烈，区域断裂带较多，整个地壳被断裂带所分割，每个部分的地壳形变都不尽相同，圣安德列斯断裂作为区域最大断裂带，活动剧烈，其周缘应变值较为突出，断裂两侧应变值差异明显，代表其周缘地壳形变差异性大，见图9。

图9 北美板块边界应变矢量分布图Fig.9 Strain vector distribution map of North American plateboundary

2.5 北美板块边界地壳形变动力学

由于板块常年缓慢漂移，北美板块边界与太平洋板块边界相互作用，导致该区域断裂发育明显，地壳结构呈碎部化，圣安德列斯大断裂为北美板块边界活动最剧烈的断裂带，断裂带活跃，相互挤压碰撞导致北美板块边界带地震活动较为频繁，见图10。

图10 北美板块边界动力学Fig.10 Boundary dynamicsof North American plate

在历史上，北美板块边界发生了多次较为著名的地震，如：1992年的Landers 7.3级地震、Joshua Tree 6.2级地震，1993年的Eureka Valley 6.1级地震，1994年Northridge 6.7级地震，1999年的Hector Mine 7.1级地震，2003年的San Simeon 6.6级地震，2004年的Parkfield 6.0级地震等[12]，见图11。

图11 北美板块边界地震活动分布Fig.11 Distribution of seismic activity on the North American plate boundary

北美板块边界城市及区域从南向北依次为：LOS洛杉矶盆地，SFE文图拉，SM圣马利亚，SSJ圣化金河南部，ST索尔顿海槽，SA圣安德列斯断层，SF圣弗朗西斯科，SAC萨克拉门托。

北美板块边界的圣安德列斯大断裂是边界带构造运动的主导断层，长期以来，太平洋板块东部的洋中脊扩张运动导致太平洋板块向东对北美板块进行挤压，圣安德列斯断层为右旋走滑，断裂两侧相对走滑运动约25～35 mm/yr，圣马利亚北部的圣安德列斯断层附近地壳运动转至东南向，致使断裂附近地壳表现为强烈压缩性质，其量值约-11×10-7/yr。2017年相关学者基于圣安德烈斯断层周缘2010～2015年五期高精度GPS监测资料分析了区域动力学特征，与本文结果基本一致[18]。Smith-Konter&Sandwell在2009年利用610个左右GPS测站的速度矢量建立模型[19]，获取了圣安德烈斯断层地壳形变特征及应变场，应变率的高值区延圣安德烈斯断层分布，特别是帕克菲尔德（Parkfield region）和索尔顿湖（Salton Sea region），其结果与本文基本一致。李英杰等人针对圣安德烈斯断层周缘的地震序列进行分析，表明圣安德列斯断裂周边易发生强震，且强震震中空间上呈现明显的"西北-东南"对称回旋迁移规律[20]。

3 结论

通过北美板块边界高精度GNSS速度场结合球面最小二乘配置方法获取了区域应变场信息，结合地质构造解释了区域地壳形变动力学机制。

（1）利用球面最小二乘配置算法求解区域应变场过程中，基于欧拉矢量扣除板块背景速度场后，准确显示了区域内复杂的相对运动速度场，利用指数函数进行协方差函数拟合后，其东向RMSE为4.61 mm2/yr2，其北向RMSE为6.85 mm2/yr2，相关系数平方分别为0.95及0.94，代表了协方差函数拟合的合理性，利用合理的协方差函数可获取区域较为准确的连续速度场信息。

（2）区域速度场整体朝向西北方向，相对运动速度有明显差异，通过对球面最小二乘配置应变场进行分析：圣马利亚北部圣安德列斯断裂周缘压缩率最大约-11×10-7/yr。圣弗朗西斯科东南部、圣马利亚西北部为剪应变高值区，量值约6×10-7/yr；圣安德列斯断裂两侧是最大主应变较高，沿大断裂分布的主应变率量值在1×10-7/yr至3×10-7/yr之间。

（3）北美板块边界地壳形变的主导断层为圣安德列斯大断裂，太平洋板块向东运动，对北美板块边界进行挤压，其板块边界圣安德列斯断裂为右旋走滑断裂，断裂周缘地壳形变剧烈，其他纵横断裂将板块边界分割成多个运动块体，地壳运动复杂且剧烈，断裂周缘的剧烈形变与地震的孕育直接相关。