地壳运动是一个极其复杂的动力学过程，许多函数模型被用来定量描述、研究地壳运动，早期的大部分模型是在板块构造理论框架下，利用长期积累的地球物理资料基于欧拉刚体定点旋转定理建立的刚性地壳运动模型(RRM)^[1-2]，认为板块运动时不发生变形或变形只发生在边界一定范围内的带状区域。事实上，地壳在运动时任何一点都是弹塑性的，板块内部是会变形的^[3]。李延兴等^[4]根据弹性应变理论，考虑板内形变提出了块体运动的整体旋转与均匀应变模型(REHSM)和整体旋转与线性应变模型(RELSM)两种弹性地壳运动模型，对中国大陆及周边地区的地壳运动及应变的分析结果也证实了这两种模型优于RRM模型对实测速度场的拟合结果^[4-6]。但REHSM和RELSM模型仅对块体内部形变作“均匀”和“线性”假设，具有一定的局限性。本文利用半参数模型中“非参数”补偿因对板内“均匀”和“线性”应变假设不合理引起的不规则运动部分进行补偿，称为半参数弹塑性地壳运动模型(SEMIREHSM)，并对中国大陆构造环境监测网络在环渤海区域的近几期GPS速度场进行拟合分析。

1 半参数弹塑性运动模型的构建 1.1 半参数模型及其补偿最小二乘估计

半参数模型一般表示为^[7]：

$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{S}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} $

(1)

式中，L为n维观测向量，其真值为${\mathit{\boldsymbol{\tilde L}}}$；B为n×m维系数矩阵；X为m维参数向量，是确定性未知量，其估值为${\mathit{\boldsymbol{\hat X}}}$，真值为${\mathit{\boldsymbol{\tilde X}}}$；S为n维非参数向量(信号)，其估值为Ŝ，真值为${\mathit{\boldsymbol{\tilde S}}}$；Δ为观测误差(噪声)，Δ~N(0，Σ_Δ)；Σ_Δ为Δ的方差阵。

式(1)的误差方程为：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{B\hat X}} + \mathit{\boldsymbol{\hat S}} - \mathit{\boldsymbol{L}} $

(2)

式中，V表示n维残差向量，通常采用如下补偿最小二乘估计准则求解：

$ {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{PV}} + a{{\mathit{\boldsymbol{\hat S}}}^{\rm{T}}} + \mathit{\boldsymbol{R\hat S}} = \min $

(3)

式中，P为n×n方阵，是观测值的权阵，R为n×n的适当选择的正定或半正定矩阵(本文仅讨论正定矩阵情形)，二次型${{\mathit{\boldsymbol{\hat S}}}^{\rm{T}}}\mathit{\boldsymbol{R\hat S}}$表示对S的某种度量，a∈[0, +∞)称为平滑参数，在极小化过程中在V与Ŝ之间起平衡作用。

顾及式(2)~(3)得式(1)的补偿最小二乘估计为:

$ \mathit{\boldsymbol{\hat S}} = {\left( {a\mathit{\boldsymbol{R}} + \mathit{\boldsymbol{M}}} \right)^{ - 1}}\mathit{\boldsymbol{ML}} $

(4)

$ \mathit{\boldsymbol{\hat X}} = {\mathit{\boldsymbol{N}}^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{\hat S}}} \right) $

(5)

式中，N=B^TPB，M=P-PBN^-1B^TP。需要说明的是，式(4)、(5)给出的参数解和系统误差估计均为有偏估计。关于补偿最小二乘解的相关统计性质及均方误差、单位权方差以及正则参数等的计算公式和方法可参阅文献[7]和[8]。

1.2 REHSM和RELSM两种弹性运动模型

RRM模型在球面坐标系下表示为：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{g}}} = {{\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right]}_{\rm{g}}} = }\\ {\left[ {\begin{array}{*{20}{c}} { - 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}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right]} \end{array} $

(6)

式中，下标g为刚性运动，V_e和V_n为一点水平运动速率的东西分量和南北分量，r为地球半径，λ和φ为块体一点的经度和纬度，ω_x、ω_y、ω_z分别为欧拉矢量在三个坐标轴上的分量，相应欧拉极计算公式为：

$ \left\{ \begin{array}{l} \omega = \sqrt {\omega _x^2 + \omega _y^2 + \omega _z^2} \\ l = a\tan \left( {\frac{{{\omega _y}}}{{{\omega _x}}}} \right)\\ \varphi = a\tan \left( {\frac{{{\omega _z}}}{{\sqrt {\omega _x^2 + \omega _y^2} }}} \right) \end{array} \right. $

(7)

式中，ω为欧拉旋转角速度，l为欧拉极经度，φ为欧拉极纬度。

假设块体内部的力学响应应变是均匀的，则REHSM模型表示为：

$ {\mathit{\boldsymbol{V}}_{\rm{h}}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right]_{\rm{h}}} = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right]_{\rm{g}}} + {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right]_{\rm{s}}} $

(8)

式中，h为一点的综合运动，s为塑性应变(形变)，V_h由块体整体刚性运动和弹塑性运动(形变)复合而成，式(8)可进一步表示为^[4]：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{h}}} = {{\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right]}_{\rm{h}}} = }\\ {\left[ {\begin{array}{*{20}{c}} { - 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}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{\varepsilon _{ee}}}&{{\varepsilon _{en}}}\\ {{\varepsilon _{ne}}}&{{\varepsilon _{nm}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {r\left( {\lambda - {\lambda _0}} \right)\cos \varphi }\\ {r\left( {\varphi - {\varphi _0}} \right)} \end{array}} \right]} \end{array} $

(9)

式中，ε_ee、ε_en、ε_ne、ε_nn为4个应变张量矩阵元素，λ₀和φ₀分别为块体几何中心的经度和纬度。比较式(6)和(9)可知，式(9)中的第1项对应式(6)中的刚性整体旋转运动部分，第2项则是伴随刚性整体旋转的均匀应变部分。在式(9)中，进一步假设4个应变参数是位置的线性函数，即RELSM模型可以表示为^[4]:

$ \begin{array}{*{20}{c}} {{{\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right]}_{\rm{h}}} = }\\ {\left[ {\begin{array}{*{20}{c}} { - 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}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{A_0}}&{{B_0}}\\ {{B_0}}&{{C_0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y \end{array}} \right] + \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {{A_1}}&{{B_2}}\\ {{B_1}}&{{C_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x^2}}\\ {{y^2}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{A_2}}&{{B_1}}\\ {{B_2}}&{{C_1}} \end{array}} \right]xy} \end{array} $

(10)

式中，A₀~A₂、B₀~B₂、C₀~C₂均为应变参数，x=r(λ-λ₀)cosφ，y=r(φ-φ₀)。

1.3 半参数弹塑性应变运动模型(SEMIREHSM)

从式(9)和(10)给出的两种弹性运动模型来看，如何处理块体内部各点偏离所在块体整体运动的趋势部分是其关键。假设REHSM模型可以确定块体运动的主要部分，偏差部分仅由对块体内部应变作“均匀”假设与实际不符引起的，则可以采用数学模型进行补偿。由于这部分偏差是模型表达不完善引起的，不同于观测噪声，因此不宜将其简单地归入随机噪声，而半参数模型可用于解决此类问题。在式(9)中增加表示偏离块体整体旋转与均匀应变的“非参数”S即可得到半参数弹塑性应变运动模型(SEMIREHSM)：

$ \begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{V}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{V}}_{\rm{n}}}} \end{array}} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} { - 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}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right] + }\\ {\left[ {\begin{array}{*{20}{c}} {{\varepsilon _{ee}}}&{{\varepsilon _{en}}}\\ {{\varepsilon _{ne}}}&{{\varepsilon _{nm}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {r\left( {\lambda - {\lambda _0}} \right)\cos \varphi }\\ {r\left( {\varphi - {\varphi _0}} \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{S}}_{\rm{n}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{e}}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_{\rm{n}}}} \end{array}} \right]} \end{array} $

(11)

式中，S_e和S_n分别表示非参数在东西向和南北向的速度分量，Δ_e和Δ_n分别表示观测速度噪声在东西向和南北向的噪声分量。

2 不同模型在GPS速度场拟合中的比较及分析

使用中国大陆构造环境监测网络在环渤海地区的GPS观测资料，研究区内(30.96°~42.54°N，111.48°~128.11°E)共有482个监测站，包括32个CORS基准站(2010~2014年)和450个区域站(一期和二期分别建设200、250个)，区域网共3期观测资料(2009、2011、2014年)，每期每站连续观测3 d，采样间隔为30 s。GPS站速率均由各观测站的时间序列通过st_filter软件进行处理获得，并进一步处理得到ITRF2008下速度场，其精度优于2 mm/a。

鉴于环渤海地区的地壳运动具有一定的长期稳定性，其运动在整体上应以刚性运动为主，即研究区内各站应具有大体一致的运动趋势。在数据处理过程中，首先根据式(6)利用最小二乘法估计板块刚体运动模型的欧拉运动参数，利用估计参数计算各站运动速率，根据站速度分量标准差大于1.5倍中误差的原则进行剔除，重复进行该步骤，最后共剔除24个站。利用剩余的458个站结果分别利用式(6)~(11)给出的4种方法进行建模，结果见表 1，表中的速率残差标准差计算方法为：

$ {S_{\Delta V}} = {\left[ {\frac{1}{{2n - R}}\left( {\sum\limits_{i = 1}^n {\left( {\Delta {V_{\rm{e}}}} \right)_i^2} + \sum\limits_{i = 1}^n {\left( {\Delta {V_{\rm{n}}}} \right)_i^2} } \right)} \right]^{\frac{1}{2}}} $

式中，ΔV_e和ΔV_n分别表示东西向和南北向速率残差分量，R为模型未知参数的个数，n为观测站数。4种模型与实际速度场相比的速率残差如图 1所示。

由图 1可以看出，4种模型计算的最大速率残差在东西向和南北向均未超过6 mm/a和5 mm/a，说明4种模型在分析和描述地壳长期运动趋势方面都是有效的；从速率残差分布的总体情况来看，RRM模型较REHSM和RELSM两种模型稍大，而后两者差别不大；而SEMIREHSM模型无论在东西分量还是南北分量上，其速率残差都明显小于上述3种模型，在-3~3 mm之间。从表 1可以看出，RRM、REHSM、RELSM 3种模型的速率残差标准差依次减小，分别为1.7 mm/a、1.5 mm/a、1.4 mm/a，而SEMIREHSM模型残差最小，为0.7 mm/a，明显小于上述3种模型，这说明SEMIREHSM较REHSM、RELSM、RRM模型更为合理。鉴于环渤海区域块体主要由坚硬的陆地岩石和海水组成，且被多种地质构造所围限，其内部构造运动及形变相当复杂，在这种情况下，对其内部运动规律作任何简单假设都将难以符合实际，即REHSM和RELSM模型简单假设其内部形变呈“均匀”和“线性”变化存在一定的局限性，而本文提出的SEMIREHSM模型假定地壳内部形变主体上呈“均匀”变化，即假定每点除了作整体运动外，还具有大致相同的均匀形变运动，但同时并不认为板内形变完全遵循“均匀”变化，而是还具有偏离“均匀”应变的形变存在，巧妙地利用“非参数”对这部分偏离“均匀”应变的形变部分加以描述。显然，假定板内形变应变“不均匀”也包括了“线性”变化的情况，因此，本文方法实际上是对REHSM和RELSM两种模型的精化，已被上述结果所证实。

3 结语

本文从模型精化的角度出发，利用半参数模型对两种弹性运动模型进行精化，并取得了满意的结果。相对于两种弹性地壳运动模型，应用于GPS速度场，新模型不仅具有更高的拟合精度，还具有良好的适应性和模型解释能力；将半参数模型引入地壳形变分析，涉及到正则矩阵的构造问题，且建模时采用了时间序列法构造正则矩阵，这等同于假定板内形变连续且具有一定的周期性。由于地震的发生和地壳的构造活动是紧密联系的，因此，若将几年或稍长尺度的地壳运动看作是处于某段地震酝酿平静期的周期运动，则采用时间序列法构造地形变分析的正则矩阵就具有一定的合理性。尽管如此，利用半参数模型进行地壳形变分析，尚需结合地壳运动的动力学机制，构造具有地学物理意义的正则矩阵，以加强半参数模型在地学数据处理领域应用的深度和广度。