基于ISE-LSC-RELSM的地壳形变分析

Crustal Deformation Analysis Based on ISE-LSC-RELSM

摘要: 将块体运动偏离整体旋转与线性应变(RELSM)的地壳形变看作“信号”，即利用“信号”对偏离RELSM的块体内部不规则形变进行补偿，建立以整体旋转与线性应变为“趋势性”的最小二乘配置模型(LSC-RELSM)，并尝试利用残差的信息扩散函数(ISE)构建方差–协方差矩阵代替以往利用高斯指数函数法确定信号的方差–协方差矩阵来解算配置模型(ISE-LSC-RELSM)。对喜马拉雅–西藏块体的实测GPS速度场数据进行拟合分析结果表明，新方法不仅优化了RELSM，并且还可以取得优于高斯指数法确定方差–协方差的拟合结果，这是因为信息扩散函数可以一定程度的抑制异常测站和粗差的影响。

Abstract: The crustal deformation of blocks deviating from global rotation and linear strain (RELSM) is regarded as “signal”, that is, the irregular deformation of blocks deviating from RELSM is compensated by “signal”, and the least square configuration model (LSC-RELSM) is established with global rotation and linear strain as “trend”. The variance-covariance matrix is constructed by using residual information diffusion function (ISE) to solve the configuration model (ISE-LSC-RELSM) instead of the variance-covariance matrix determined by Gaussian exponential function method. The results of fitting analysis on the measured GPS velocity field data of himalaya-Tibet block show that the new method not only optimizes RELSM, but also achieves better fitting results than the Gaussian exponential method to determine the variance-covariance, because the information diffusion function can suppress the influence of abnormal stations and gross errors to a certain extent.

1. 引言

传统地壳运动模型大都是基于欧拉定理建立的刚体运动模型，这些模型认为板块内部不存在或不考虑内部形变。但事实上，大量研究结果表明，尽管板块运动确实以刚性运动为主，但其内部也存在不同程度的不规则形变 [1] [2] [3]。我国学者李延兴等研究了板块内部形变性质，提出了块体运动的整体旋转与线性应变模型(RELSM) [4]。该模型虽然考虑了块体内部形变，但假设块体内部形变呈线性变化，对于地质构造运动复杂的区域仍有局限性。因此，有学者提出利用最小二乘配置模型改进整体旋转与线性运动模型的方法，即将偏离RELSM的不规则运动部分用信号加以描述，建立利用最小二乘配置优化的整体旋转与线性应变模型 [5]，本文简记为LSC-RELSM。然而，无论是经典刚性最小二乘模型(RRM)还是最小二乘配置模型都易受异常测站或粗差的影响，所以，应考虑建立具有抗差能力的配置模型。柴洪洲等利用等价权迭代方法建立了抗差最小二乘配置模型进行地壳形变分析 [6] [7]，但抗差过程易受最小二乘残差“均摊”效应影响而存在对异常测站影响处理不彻底的问题。本文尝试根据RELSM拟合残差，利用信息扩散估计(ISE)构建方差–协方差矩阵，由于信息扩散函数根据残差分布构建方差协方差矩阵，理论上具有抗差性，故可以解决异常测站对速度场拟合带来的不利影响。对喜马拉雅–西藏块体的地壳形变分析结果表明，该方法可以相对彻底地消除异常测站影响，从而优化LSC-RELSM，本文将该方法简记为ISE-LSC-RELSM。

2. 信息扩散函数构建方差–协方差矩阵原理

2.1. 信息扩散函数原理及其估计

2.1.1. 信息扩散函数定义

设母体K的概率密度函数为 $f_{\left(x\right)}$，W为来自母体K的样本。当由W的数据不能完全精确地解释 $f_{\left(x\right)}$ 时，称W对K来说是非完备的 [8]。当W非完备时，每个样本点 $w_i$ 只是“周围”的代表，不能完全反映出母体的全部特征。此时，如果通过增加样本点数目就能改善不完备的程度，使其趋于或达到完备，即能完全反映母体全部特征。由非完备到完备具有一个过渡阶段，当W非完备时，每一个样本点可以发展成多个样本点。设 $w_i$ 的观测值为 $l_i$，样本W在观测值 $l_i$ 点提供的信息理论上应被其周围点所分享，且越靠近观测值 $l_i$ 点所分享到的信息就越多。我们把从观测值 $l_i$ 所分享得到的这种信息称为从 $l_i$ 扩散来的信息，而将观测值 $l_i$ 点的信息被周围点分享的过程称为信息扩散过程，简称信息扩散 [8]。

2.1.2. 信息扩散估计

设 $W=\left\{w_2,w_2,\cdots ,w_n\right\}$ 是知识样本，U是基础论域，记 $w_i$ 的观测值为 $v_i$，设 $x=h\left(v-v_i\right)$，样本W非完备时存在函数 $\mu \left(x\right)$ 的量值扩散到l，且扩散所得到的原始信息分布

$Q\left(x\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}\mu \left(x\right)}={\displaystyle \underset{j=1}{\overset{n}{\sum }}\mu \left(h\left(v-v_i\right)\right)}$ (1)

能更好地反映出样本W所在总体的规律，我们把这一原理称为信息扩散原理，根据信息扩散原理对母体K概率密度函数的估计称为扩散估计 [8]。

设 $\mu \left(x\right)$ 为定义在 $\left(-\infty ,\infty \right)$ 区域上的一个波雷尔可测函数， ${\Delta }_n>0$ 为常数，n为子样容量，则称

$\hat{f}\left(l\right)=\frac{1}{n{\Delta }_n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}\mu \left[h\left(\frac{l-l_j}{{\Delta }_n}\right)\right]}$ (2)

为母体K的概率密度函数 $f\left(x\right)$ 的一个扩散估计。式中， $\mu \left(x\right)$ 称为扩散函数； ${\Delta }_n$ 为窗宽 [7]。

2.1.3. 信息扩散函数及窗宽确定

由(2)式可知获得母体K概率密度函数的扩散估计关键问题是要确定扩散函数 $k\left(x\right)$ 的具体表达形式，不同的扩散函数所得到的扩散估计也不相同。由于信息扩散的方式与分子扩散非常相似，则类比分子扩散理论确定扩散函数为 [8]：

$\mu \left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{x^2}{2{\sigma }^2}\right)$ (3)

联立(2)式和(3)式得到母体概率密度函数的正态扩散估计为：

$\hat{f}\left(l\right)=\frac{1}{nh\sqrt{2\pi }}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\exp \left(-\frac{{\left(v-v_i\right)}^2}{2h^2}\right)}$ (4)

式中，h是扩散系数，其值的确定式如下 [8] [9]：

$\{\begin{array}{l}h=\frac{\alpha \left(b-a\right)}{n-1}\\ a=\min \left(l_i\right)\\ b=\max \left(l_i\right)\\ i=1,2,\cdots ,n\end{array}$ (5)

式中， $\alpha$ 是n的函数，其取值随n值变化情况见参考文献 [8]，若样本点数目 $n\ge 17$ 时， $\alpha$ 值均为1.420693101。

2.2. 块体的整体旋转与线性应变模型(RELSM)

假设块体内部应变呈线性变化，则其模型为 [2]：

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

式中， $x=r\cos \phi \left(\lambda -{\lambda }_0\right)$， $y=r\left(\phi -{\phi }_0\right)$ ； $V_e,V_n$ 是块体上任意一点的东向与北向速度； $\lambda$ 与 $\phi$ 是该点的经度和纬度；r是地球的半径； ${\lambda }_0$ 与 ${\phi }_0$ 是块体几何中心的经度与纬度； ${\omega }_x,{\omega }_y$ 和 ${\omega }_z$ 是块体的旋转矢量； $A_0,A_1,A_2,B_0,B_1,B_2,C_0,C_1,C_2$ 是块体的应变率参数。

2.3. 基于整体旋转与线性应变的最小二乘配置模型(LSC-RELSM)

最小二乘配置模型的一般形式为 [6]：

$L=AX+Y+\Delta$ (7)

式中，A为倾向参数系数矩阵，X为倾向参数， $Y={\left[\begin{array}{cc}S& S^{\prime }\end{array}\right]}^{\text{T}}$ 为信号， $\Delta$ 为观测噪声。当(7)式的主项取为(6)式时，即为以整体旋转与线性应变为整体趋势项的地壳形变分析的最小二乘配置模型，此时信号Y表示偏离RELSM的不规则形变，也可看作是对整体旋转与线性应变趋势的随机扰动。

已知观测数据为等精度观测，已测点信号S的先验数学期望、未知测站点信号 $S^{\prime }$ 的先验数学期望和观测噪声 $\Delta$ 的先验数学期望均为0。假设信号与观测噪声之间不相关，即

$\{\begin{array}{l}{D}_{S\Delta }=0\\ D_{S^{\prime }\Delta }=0\end{array}$ (8)

信号Y的协方差阵为：

$D_Y=\left[\begin{array}{cc}{D}_{SS}& D_{S{S}^{\prime }}\\ D_{S^{\prime }S}& D_{S^{\prime }{S}^{\prime }}\end{array}\right]$ (9)

通过间接平差得到参数为：

$\hat{X}={\left(A^{\text{T}}{P}_{L}A\right)}^{-1}{A}^{\text{T}}{P}_{L}l$ (10)

$\hat{Y}=D_Y{B}^{\text{T}}P\left(l-A\hat{X}\right)$ (11)

式中， $D_L=D_Y+D_{\Delta }=P_L^{-1}$， $D_Y$ 是观测向量协方差； $D_{\Delta }$ 是观测误差自协方差阵。

已知测站点和未知测站点的信号估值为：

$\begin{array}{l}\hat{S}=D_S{B}^{\text{T}}{P}_L\left(l-A\hat{X}\right)\\ {\hat{S}}^{\prime }=D_{S^{\prime }S}{B}^{\text{T}}{P}_L\left(l-A\hat{X}\right)\end{array}$ (12)

式中， $D_S$ 为已知测站点协方差阵； $D_{S^{\prime }S}$ 为已知测站点和未知测站点之间的协方差阵。在地壳形变分析的配置模型(7)式中，信号的方差–协方差矩阵通常是采用某种随距离衰减的经验协方差函数来确定，很多研究表明高斯指数函数是一种比较适合于地壳形变分析的协方差函数，其具体形式为 [6]：

$D\left(S\right)=D\left(0\right)\exp \left(-K^2{S}^2\right)$ (13)

式中，S表示两个信号间的距离； $D\left(S\right)$ 表示两个信号间的协方差； $D\left(0\right)$ 为两个信号距离为0时的协方差，即信号的方差；K为待求参数，参数具体求解过程详见文献 [2] 和 [6]。

3. ISE-LSC-RELSM原理

由于(13)式中，信号的方差–协方差矩阵实际是根据一点在扣除整体旋转与线性应变后的剩余速率来计算的，则协方差函数参数求解必然受到测站异常速率和粗差的共同影响，为此，需要考虑采取措施抑制这种影响。文献 [6] 利用欧亚板块运动的欧拉矢量对GPS监测速度场进行中心化(扣除欧亚板块刚体运动趋势)后，采用IGGI稳健估计方法抑制粗差影响，取得了较好的效果，但这种方法不仅需要迭代求解，而且迭代过程易受最小二乘固有的“均摊”效应影响而错误判定异常测站，导致结果不精确。本文尝试根据整体旋转与线性应变最小二乘估计残差值，利用式(1)~式(5)给出的信息扩散函数及其估计公式，代入利用式(6)式计算的整体旋转与线性应变最小二乘残差直接计算(不通过高斯指数等协方差函数)并构建信号的方差–协方差矩阵。由于残差的信息扩散函数能事先准确地确定残差分布，故每个方差–协方差矩阵元素本身具有较好的抗差性，从而可以解决高斯指数法计算方差–协方差矩阵对异常测站检测不彻底或对异常测站误判的问题。

4. 算例及分析

4.1. 数据来源及处理

本文采用的算例数据是(王琪等) 2001年发表于《Science》(文献 [10] )的中国大陆地壳运动早期区域变形监测数据中的喜马拉雅山–西藏板块相对于欧亚板块GPS速度场数据，共计54个测点速度场数据。为验证基于信息扩散函数的LSC模型的优势，现采用以下三种方案对喜马拉雅山–西藏板块相对于欧亚板块GPS速度场数据进行拟合分析：

方案I：块体的整体旋转与线性应变模型拟合分析(RELSM)；

方案II：最小二乘配置模型拟合分析(LSC-RELSM)；

方案III：基于信息扩散函数的最小二乘配置模型拟合分析(ISE-LSC-RELSM)。

上述三种方案中，均采用速率拟合残差的均方根误差RMS作为精度评价指标，其表达式为：

$\text{RMS}=\sqrt{\frac{V^{\text{T}}V}{n}}$ (14)

式中，V为速率残差，即拟合值与观测值之差；n为拟合点个数。

三种方案的欧拉运动参数及均方根误差统计结果列入表1，为直观，还同时给出了各测站点在东西向和南北向的速度残差图，如图1、图2所示。

Table 1. Euler motion parameters and root mean square error of each scheme

表1. 各方案的欧拉运动参数和均方根误差

Figure 1. East direction velocity fitting residual

图1. 东方向速率拟合残差

Figure 2. South-North velocity fitting residual

图2. 南北方向速率拟合残差

4.2. 结果分析

(1) 从表1统计的结果可以看出，方案I的均方根误差在东方向和北方向上分别为0.3295 mm/a、0.1025 mm/a；方案II的均方根误差在东方向和北方向上分别为0.2885 mm/a、0.1016 mm/a；方案III的均方根误差在东方向和北方向上分别为0.1745 mm/a、0.0973 mm/a。上述结果表明，方案III的总体拟合精度高于方案I和方案II，且基于信息扩散函数确定信号的方差–协方差矩阵对速度场的拟合结果明显优于利用高斯指数法确定方差–协方差的拟合结果。

(2) 从图1可以反映出：在东方向上方案I和方案II拟合的残差均没有超过1 mm/a，且大部分集中在0.4 mm/a内；方案III拟合的残差均没有超过0.6 mm/a，且大部分集中在0.3 mm/a内。图2可以反映出：在北方向上方案I、方案II和方案III拟合的残差均没有超过0.3 mm/a，且大部分集中在0.2 mm/a内。三种方案的拟合残差均较小，说明这三种方法均能描述喜马拉雅–西藏块体的地壳水平运动。方案III相对于其余两种方案在北方向拟合残差相差不大的原因是地壳形变在北方向上形变量较小，但总体还是优于其他两种方案，说明本文采用的方法结果更精确。

5. 结论

本文将板块内部偏离整体旋转与线性应变的不规则形变视为信号，利用信号对偏离块体的整体旋转与线性应变部分形变进行了补偿从而达到改善RELSM的目的。并且利用信息扩散函数构建方差–协方差矩阵，由于信息扩散函数可以事先确定残差(扣除整体旋转与线性应变趋势运动的剩余速率)分布，因此，这种方法从模型参数解算的角度来看，相当于建立了抗差最小二乘估计模型，最终达到消除异常测站(粗差)影响的目的，其意义在于：

(1) 根据扣除整体趋势运动后的剩余速率的实际分布构建信号的方差–协方差矩阵，可以较好地抑制异常测站或粗差带来的不良影响，为地壳形变分析的抗差最小二乘配置模型解算提供了一种新思路；

(2) 通常，无论哪种模型，例如高斯指数法构建信号的方差–协方差矩阵，一般都假定地壳岩石圈具有连续介质，即假定地壳形变具有连续性，但事实上，由于隐伏断层等复杂构造的存在，有时上述假定难以满足，因此，在形变连续性假设下的一切分析结果都将难以解释。而本文方法，利用信息扩函数构建方差–协方差矩阵，并没有假定介质或形变具有连续性，其结果具有更好的可解释性。

基金项目

贵州大学2020年SRT项目(贵大SRT字(2021)191号)；贵州省科学技术基础研究计划项目(黔科合基础[2017]1054)；贵州大学测绘科学与技术研究生创新实践基地建设项目(贵大研CXJD[2014]002)。

NOTES

^*第一作者。

^#通讯作者。