基于多面函数的环渤海区域地壳水平运动速度场模拟

Simulation of Velocity Field of Horizontal Crustal Movement around the Bohai Sea Based on Polyhedral Function

摘要: 利用2009~2014年间中国大陆构造环境监测网在环渤海及其邻区的多期GPS对地观测速度场数据，分别采用多面函数法、刚体模型(RRM)、整体旋转与均匀应变模型(REHSM)及整体旋转与线性应变模型(RELSM)对速度场数据进行拟合和推估。对四种方法的计算结果进行对比分析，结果表明无论拟合精度还是推估精度，多面函数法都要优于其他三种模型。在此基础上，应用多面函数法按0.5˚ × 0.5˚的间隔进行速度场插值，根据插值结果绘制了环渤海区域地壳运动的水平速度场，并计算出环渤海区域相对于ITRF2008，其东向速度为30.2844 mm/yr，北向速度为−11.2871 mm/yr，整体上以32.3194 mm/yr的速度往东南20.4406˚方向运动。

Abstract: Using the Chinese mainland tectonic environment monitoring network in the Bohai Sea and its adjacent areas, the velocity field data of multi-phase GPS observation in the 2009~2014 years were fitted and estimated by the multifaceted function method, rigid body model (RRM), the whole rotation and homogeneous strain model (REHSM), and the overall rotation and linear strain model (RELSM). The results show that the polyhedral function method is better than the other three models in terms of fitting accuracy and estimation accuracy. On this basis, the velocity field is interpolated at intervals by using the polyhedral function method. According to the interpolation results, the horizontal velocity field of crustal movement in the Bohai Rim region is drawn. Compared with ITRF2008, the east velocity is 30.2844 mm/a and the north velocity is −11.2871 mm/a. On the whole, it moves to the southeast at the speed of 32.3194 mm/a.

1. 引言

针对板块刚性运动模型(RRM) [1] [2] 无法顾及板块内部形变的问题，本世纪初，李延兴等学者做了大量研究，先后假设板块内部形变呈均匀变化和线性变化，建立起了板块水平运动的整体旋转与均匀应变模型(REHSM) [3] 和整体旋转与线性应变模型(RELSM) [4]，并通过实例验证了两种弹塑性板块运动模型的拟合结果都要优于刚性模型，并且把板块内部形变假设成线性变化比假设成均匀变化更为合理 [4] [5]。多面函数拟合方法是由美国的Hardy [6] 教授提出，刘经南 [7]、曾安敏 [8] 先后利用多面函数建立了分布均匀、有实用价值的中国地壳水平运动整体速度场。近年来，姚未正等学者利用GNSS速度场探究地壳形变对地质构造的影响 [9] [10] [11] [12]。对多面函数拟合法对板块是否是刚性没有要求，为了比较弹性板块运动模型和多面函数拟合法哪一种方法更适合环渤海区域地壳水平运动速度场模型构建，本文利用2009~2014年间中国大陆构造环境监测网在环渤海及其邻区的多期GPS对地观测速度场数据，分别采用刚体模型(RRM)、多面函数法、整体旋转与均匀应变模型(REHSM)及整体旋转与线性应变模型(RELSM)对速度场数据进行拟合推估，分析结果表明多面函数更适合环渤海区域地壳水平运动模型构建。

2. 板块刚性–弹性运动模型

传统板块运动模型假定板块内部没有变形，且基于欧拉刚体旋转运动模型，容易推导出板块的刚性旋转运动模型(RRM)为 [1] [2]：

$\left[\begin{array}{c}{V}_e\\ V_n\end{array}\right]=\left[\begin{array}{ccc}-r\cos \lambda \sin \varphi & -r\sin \lambda \sin \varphi & r\cos \varphi \\ r\sin \varphi & -r\cos \varphi & 0\end{array}\right]\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]$ (1)

式(1)中， $V_e$ 和 $V_n$ 分别表示东西向和南北向速度分量，r为块体的平均地心距，实际计算时，以地球的平均半径带入即可， $\lambda$ 和 $\phi$ 为测站的经纬度， ${\omega }_x,{\omega }_y,{\omega }_z$ 为板块旋转角速度在三个坐标轴上的分量大小。

式(1)是在假设块体内部没有形变的条件下得出的，但事实上，板块内部并非绝对干刚体，故在做整体趋势性运动的同时，必将伴随不同程度的内部形变。若假设块体内部形变碎位置呈均匀变化，则可得到一种刚性弹塑性运动模型，即块体的整体旋转与均匀应变模型(REHSM)为 [3]：

$\left[\begin{array}{c}{V}_e\\ V_n\end{array}\right]=\left[\begin{array}{ccc}-r\cos \lambda \sin \phi & -r\sin \lambda \sin \phi & r\cos \phi \\ r\sin \lambda & -r\cos \lambda & 0\end{array}\right]\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]+\left[\begin{array}{cc}{\epsilon }_{ee}& {\epsilon }_{en}\\ {\epsilon }_{ne}& {\epsilon }_{nn}\end{array}\right]\left[\begin{array}{c}r\left(\lambda -{\lambda }_0\right)\cos \phi \\ r\left(\phi -{\phi }_0\right)\end{array}\right]$ (2)

式(2)中， ${\epsilon }_{ee},{\epsilon }_{ne},{\epsilon }_{nn}$ 分别为测站点的三个应变参数， ${\lambda }_0,{\phi }_0$ 为监测点所在块体几何中心的经纬度，其他参数意义同式(1)。上述模型假定块体内部存在形变，从理论上优于刚体运动模型，但假定块体内部形变是均匀变化的，在复杂构造区域是不适用的，具有局限性。若进一步假设块体内部形变应变呈线性变化，则易得块体的整体旋转与线性应变模型(REHLM)为 [4] [5]

$\left[\begin{array}{c}{V}_e\\ V_n\end{array}\right]=\left[\begin{array}{ccc}-r\cos \lambda \sin \varphi & -r\sin \lambda \sin \varphi & r\cos \varphi \\ r\sin \lambda & -r\cos \lambda & 0\end{array}\right]\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]+\left[\begin{array}{cc}{A}_0& B_0\\ B_0& C_0\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\frac{1}{2}\left[\begin{array}{cc}{A}_1& B_2\\ B_1& C_2\end{array}\right]\left[\begin{array}{c}{x}^2\\ y^2\end{array}\right]+\left[\begin{array}{cc}{A}_2& B_1\\ B_2& C_1\end{array}\right]xy$ (3)

式中， $x=r\left(\lambda -{\lambda }_0\right)\cos \varphi$ 、 $y=r\left(\phi -{\phi }_0\right)$， $A_0~A_2$ 、 $B_0~B_2$ 、 $C_0~C_2$ 为待求的应变参数，其具体含义见文献 [5]，其余参数意义同(1)式。

3. 多面函数拟合法

多面函数拟合方法是由美国的Hardy教授提出的，其思想是：任何一个连续广滑的数学表面都可以利用一系列规则的数学函数的总和以任意精度逼近。多面函数拟合方法基本原理如下 [6] [7]：

设参与多面函数拟合的测站点有m个， $V\left(\lambda ,\phi \right)$ 为测站点 $\left(\lambda ,\phi \right)$ 的水平运动速度，则多面函数拟合方法可以用模型表示成如下形式

$V\left(\lambda ,\phi \right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\alpha }_{i}Q\left(\lambda ,\phi ;{\lambda }_i,{\phi }_i\right)}$ (4)

式中，n为选取的结点个数； ${\alpha }_i\left(i=1,2,\cdots ,n\right)$ 为待求系数； $Q\left(\lambda ,\phi ;{\lambda }_i,{\phi }_i\right)$ 为核函数。核函数的形式如下

$Q\left(\lambda ,\phi ;{\lambda }_i,{\phi }_i\right)={\left[{\left(\lambda -{\lambda }_i\right)}^2+{\left(\phi -{\phi }_i\right)}^2+{\delta }^2\right]}^{\beta }$ (5)

式中， $\delta$ 为光滑因子，用来对核函数进行调整； $\beta$ 为指数因子，通常可取1/2 (正双曲面)、−1/2 (负双曲面)、3/2 (三次曲面)等。

将(4)式写成矩阵形式，如下

$\underset{m\times 1}{V}=\underset{m\times n}{Q}\underset{n\times 1}{\alpha }+\underset{m\times 1}{v}$ (6)

即

$\left[\begin{array}{c}{V}_1\left({\lambda }_1,{\phi }_1\right)\\ ⋮\\ V_m\left({\lambda }_m,{\phi }_m\right)\end{array}\right]=\left[\begin{array}{ccc}Q\left({\lambda }_1,{\phi }_1;{\lambda }_{i1},{\phi }_{i1}\right)& \cdots & Q\left({\lambda }_1,{\phi }_1;{\lambda }_{n1},{\phi }_{n1}\right)\\ ⋮& \ddots & ⋮\\ Q\left({\lambda }_m,{\phi }_m;{\lambda }_{i1},{\phi }_{i1}\right)& \cdots & Q\left({\lambda }_m,{\phi }_m;{\lambda }_{n1},{\phi }_{n1}\right)\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ ⋮\\ {\alpha }_n\end{array}\right]+\left[\begin{array}{c}{v}_1\\ ⋮\\ v_m\end{array}\right]$ (7)

根据(7)式按照最小二乘法则对待估参数进行求解，求得待估参数后即可对未知点进行推估。运用多面函数进行拟合推估时，其关键在于结点、光滑因子和指数因子的选取，因此需要根据实际情况多次调试后最终确定。

4. 算列及分析

本文数据源自文献 [5]，其精度优于2 mm/yr。

鉴于地壳运动一般都具有长期稳定性，其运动在整体上应以刚性运动为主，为此在构建速度场拟合模型前应先将那些显著偏离块体整体运动趋势的测站予以剔除，以提高速度场拟合模型自身的精度及推估的合理性和可靠性。具体过程为：首先利用式(1)给出的地壳刚性运动模型计算所有测站在东西和南北方向的速率残差，以3倍验后中误差做为取舍准则，将平差后速率残差大于3倍中误差的测站剔除，利用余下测站重复进行上述操作，直到剩余测站全部满足限差要求为止。通过上述方法，本文共剔除测站18个，剩余测站383个。为验证各种模型的速度场拟合和推估效果，首先从383个测站中选取45个测站点作为推估点，其余点作为模型拟合点，分别采用以下四种方案对环渤海区域进行速度场拟合和推估分析：

方案1：采用RRM对环渤海区域地壳水平运动速度场进行拟合和推估；

方案2：采用REHSM对环渤海区域地壳水平运动速度场进行拟合和推估；

方案3：采用RELSM对环渤海区域地壳水平运动速度场进行拟合和推估；

方案4：采用多面函数法对环渤海区域地壳水平运动速度场进行拟合和推估。

在方案4中，首先从383个测站中选取45个测站点作为推估点，用余下的338个测站进行速度场拟合。在多面函数拟合中，根据数据的分布情况，本文较为均匀的选取66个点作为结点，经过多次调试，最终选择 $\delta =0.1$， $\beta =1/2$。

上述方案中，用速率残差的均方根误差RMSE作为精度评价指标，计算公式为：

$\text{RMSE}=\sqrt{\frac{V^{\text{T}}V}{n}}$ (9)

上式中n为参与拟合点的个数，V为速率残差，即拟合值与观测值之差；推估时，n为参与推估点的个数，V为推估值与观测值之差。表1统计了四种方案拟合和推估的均方根误差，图1为研究区域GPS速度场，图2~5为四种方案的拟合残差图，图6、图7为四种方案东方向(E)和北方向(N)的推估残差图。

Figure 1. GPS velocity field in the Bohai Sea rim region

图1. 环渤海区域GPS速度场

Table 1. Statistics of root mean square error RMSE (mm/yr) of four schemes

表1. 四种方案均方根误差RMSE (mm/yr)统计

Figure 2. Fitting residual of scheme 1

图2. 方案1拟合残差

Figure 3. Fitting residual of scheme 2

图3. 方案2拟合残差

Figure 4. Fitting residual of scheme 3

图4. 方案3拟合残差

Figure 5. Fitting residual of scheme 4

图5. 方案4拟合残差

Figure 6. Estimated residuals of the four schemes in the east-west direction

图6. 四种方案东西方向(E)推估残差

Figure 7. Estimated residuals of the four schemes in the north-south direction

图7. 四种方案南北方向(N)推估残差

5. 结果分析

1) 从图1~4可以看出，多面函数法拟合的残差在东方向上未超过4 mm/yr，北方向未超过3 mm/yr，RRM、REHSM和RELSM在东方向上的拟合残差未超过5 mm/yr，北方向上除个别点外，拟合残差都在3 mm/yr以内，说明四种方案拟合的结果都是有效的；从图5、图6推估点残差图来看，四种方案在东方向上的推估残差都在5 mm/yr以内，北方向上的推估残差都在3 mm/yr以内，表明四种方案的推估结果都是合理的。

2) 从表1拟合的均方根误差统计结果来看，多面函数拟合的东方向均方根差为1.3593 mm/yr、北方向为0.8036 mm/yr，RRM的为1.6090 mm/yr和0.9805 mm/yr，REHSM的为1.5911 mm/yr和0.9276 mm/yr，RELSM的为1.5819 mm/yr和0.9102 mm/yr，无论东方向还是北方向，多面函数拟合的均方根误差都要小于其他三种方案。

3) 从表1推估的均方根误差统计结果来看，多面函数推估的东方向均方根误差为1.8659 mm/yr、北方向为1.0406 mm/yr，RRM的为1.9508 mm/yr和1.1204 mm/yr，REHSM的为1.9414 mm/yr和1.0601 mm/yr，RELSM的为1.9334 mm/yr和1.0576 mm/yr，无论是东方向还是北方向，多面函数推估的均方根误差都要小于其他三种方案。

以上分析结果表明，多面函数法更适合用来建立本研究区域的地壳水平运动速度场。根据多面函数的拟合参数，按0.5˚ × 0.5˚的间隔对环渤海区域进行速度场插值，根据插值结果绘制了环渤海区域的均匀速度场图，如图8所示。基于插值速度场，计算出了环渤海区域的整体运动速率，结果见表2。

Figure 8. 0.5˚ × 0.5˚ interpolation velocity field

图8. 0.5˚ × 0.5˚插值速度场

Table 2. Overall movement speed in the Bohai Sea rim region

表2. 环渤海区域整体运动速度

从表2统计的结果来看，环渤海区域相对于ITRF2008，其东向速度为30.2844 mm/yr，北向速度为−11.2871 mm/yr，整体上以32.3194 mm/yr的速度向SE 20.4406˚方向运动。

6. 结束语

本文利用环渤海区域的多期GPS速度场数据，分别采用多面函数拟合法、整体旋转与均匀应变模型和整体旋转与线性应变模型对环渤海区域的地壳水平运动速度场进行了拟合分析。结果表明多面函数拟合法取得了优于其他模型拟合的结果。究其原因，主要是环渤海区域被多种构造所围限，局部运动差异较大。多面函数法实际上采用分片分区域拟合并串联成整个区域的速度场，兼顾了局部差异与整体趋势，因此取得了较好的拟合结果。不仅如此，多面函数模型还可以通过内插和推估给出一定分辨率的均匀速度场图像，这可以很好地缓解监测点数目不多或测站分布不均匀带来的负面影响，尤其是可以为缺乏监测数据的区域研究提供一定具有实用价值的先验信息，这对于研究区域块体的动态变化机制研究具有积极意义。

基金项目

1) 贵州省教学改革项目(2021020)；2) 贵州大学培育项目(贵大培育[2020]57)；3) 国家自然科学基金项目(41701464)；4) 贵州省科学技术基础研究计划项目(黔科合基础[2017]1054)；5) 贵州大学测绘科学与技术研究生创新实践基地建设项目(贵大研CXJD[2014]002)。

NOTES

^*第一作者。

^#通讯作者。