自板块构造运动理论提出以来，许多学者基于欧拉定理和地球物理资料，建立了一批全球板块运动模型^[1-3]。这些模型大多将板块运动看成为一种刚性运动，认为板块运动时其内部不发生变形，或形变只发生在板块边缘一定条带范围内。基于该认识，在建立块体运动模型时，许多学者倾向于将偏离整体运动趋势的点视为异常点并加以剔除。事实上，板块运动时其边缘具有发散性和一定弹塑性^[4-5]。李延兴等^[6-8]提出块体的弹性运动模型，主要包括整体旋转与均匀应变模型(REHSM)和整体旋转与线性应变模型(RELSM)，这2种模型都已考虑板内形变，是对地壳刚性运动模型(RRM)的发展，对中国大陆及周边地区的地壳运动分析结果也证实这2种模型优于以往刚性运动模型对实测速度场的拟合结果。由此可见，块体运动整体上以刚性运动为主，具有明显的整体运动趋势和长期稳定性特点，但由于板内构造活动，板内存在与整体运动趋势不同的形变运动。因此，当存在足够数目的观测点时，异常测站的剔除应适当考虑块体内部形变，同时为防止粗差对分析结果的影响，还应采取抗差措施，以使速度场模型更趋于实际。

本文利用整体旋转与线性应变模型(RELSM)描述块体内部形变，采用高崩溃污染率抗差估计方法对中国大陆构造环境监测网络在环渤海区域的GPS观测速度场数据进行测站筛选及拟合分析，获得优于基于刚性模型进行测站筛选与拟合的结果。

1 块体的刚性运动模型和整体旋转与线性应变模型 1.1 块体的刚性运动模型(RRM)

根据欧拉定点旋转定理，板块刚性运动模型(RRM)可表示为^[1-2]：

$ \left[ {\begin{array}{*{20}{c}} {{V_E}}\\ {{V_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - r{\rm{cos}}\lambda {\rm{sin}}\varphi }&{ - r{\rm{sin}}\lambda {\rm{sin}}\varphi }&{r{\rm{cos}}\varphi }\\ {r{\rm{sin}}\varphi }&{ - r{\rm{cos}}\varphi }&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right] $

(1)

式中，V_E和V_N分别为测站E向速度分量和N向速度分量，r为地球半径，λ和φ分别为测站大地经纬度，ω_x、ω_y、ω_z为板块欧拉旋转矢量在3个坐标轴上的分量。根据欧拉旋转矢量中3个分量可计算板块运动的欧拉极所对应的经度l、纬度φ及旋转角速度ω，计算公式为：

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

(2)

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

研究表明，在几十年时间尺度内块体运动具有弹性运动性质^[9]，如果将块体看成是由具有连续介质的固态岩石圈组成，则某点的运动应由块体整体刚性旋转运动和局部弹性运动复合而成。假设块体内部形变为均匀变化，则可得到块体水平运动的整体旋转与均匀应变模型(RELSM)为^[6]：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{V_E}}\\ {{V_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - r{\rm{cos}}\lambda {\rm{sin}}\varphi }&{ - r{\rm{sin}}\lambda {\rm{sin}}\varphi }&{r{\rm{cos}}\varphi }\\ {r{\rm{sin}}\lambda }&{ - r{\rm{cos}}\lambda }&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {{\varepsilon _{ee}}}&{{\varepsilon _{en}}}\\ {{\varepsilon _{ne}}}&{{\varepsilon _{nn}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {r(\lambda - {\lambda _0}){\rm{cos}}\varphi }\\ {r(\varphi - {\varphi _0})} \end{array}} \right] \end{array} $

(3)

式中，ε_ee、ε_ne、ε_nn为3个待求应变参数，λ₀、φ₀为研究块体所包含区域的几何中心的大地经纬度。进一步假设块体内部呈线性变化，则可得到块体水平运动的整体旋转与线性应变模型(RELSM)为^[7]：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{V_E}}\\ {{V_N}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - r{\rm{cos}}\lambda {\rm{sin}}\varphi }&{ - r\sin \lambda {\rm{sin}}\varphi }&{r{\rm{cos}}\varphi }\\ {r{\rm{sin}}\lambda }&{ - r{\rm{cos}}\lambda }&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _x}}\\ {{\omega _y}}\\ {{\omega _z}} \end{array}} \right] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {{A_{\rm{0}}}}&{{B_{\rm{0}}}}\\ {{B_{\rm{0}}}}&{{C_{\rm{0}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y \end{array}} \right] + \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {{A_{\rm{1}}}}&{{B_{\rm{2}}}}\\ {{B_{\rm{1}}}}&{{C_{\rm{2}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x^2}}\\ {{y^2}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{A_{\rm{2}}}}&{{B_{\rm{1}}}}\\ {{B_{\rm{2}}}}&{{C_{\rm{1}}}} \end{array}} \right]xy \end{array} $

(4)

式中，x=r(λ－λ₀)cosφ，y=r(φ－φ₀)，A₀~A₂、B₀~B₂、C₀~C₂为9个待求应变参数。

研究表明^[7]，REHSM与RELSM两种模型可同时顾及到板块的整体运动趋势和内部形变，模型的拟合精度优于RRM。但REHSM过于绝对地假设块体内部的均匀形变，RELSM则假设块体内部形变呈线性变化，可顾及板内形变的空间差异性，更接近实际地壳运动模型。

2 基于RELSM的抗差估计测站筛选及速度场拟合方法

研究表明，区域内某点的地壳运动应为趋势性刚体运动与内部不规则形变的复合形式，在建立区域地壳运动模型时，既要充分考虑研究区所在块体的整体运动趋势的长期稳定性(块体内部各点运动与块体整体运动的一致性)，又要考虑块体所在区域的构造活动所引起的不规则形变运动(各点运动的局部差异性)。因此，在建模分析测站筛选过程中，可采用能兼顾整体性和局部差异性的异常测站剔除机制来平衡某点运动的整体性和局部差异性。研究表明^[10]，在同一运动块体中，相对于由全体观测站组成的整体速度场，尽管异常测站表现出与其他站点不一致的运动特征，但一般认为各点运动在空间的分布具有连续性，该点的运动出现异常应具有过渡性。因此，在对异常测站进行筛选时，可采取类似于粗差处理的方式，按基于等价权理论的抗差估计进行处理^[11]。另外，考虑到引起块体内部形变的原因和动力学机制及各点的异常程度都不明确，本文采用高崩溃污染率抗差估计方法进行测站筛选，该方法的优点为只要异常点不超过50%都能进行判别^[10]。具体步骤为^[12]：

1) 取所有测站速率值的中位数med(L)，将各测站速率值与中位数相减，求得初始残差为：

$ \Delta {L_i} = {L_i} - {\rm{med}}(L) $

(5)

一般情况下，E向和N向的运动具有明显差异，因此需在2个方向分别进行计算。

2) 计算Δ L的初始均方差因子抗差解：

$ {\sigma _{\Delta \mathit{\boldsymbol{L}}}} = {\rm{med}}(\left| {\Delta \mathit{\boldsymbol{L}}} \right|)/0.674{\kern 1pt} {\kern 1pt} {\kern 1pt} 5 $

(6)

3) 取强淘汰函数：

$ \bar P_i^0 = \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1&{\left| {\Delta {L_i}/{\sigma _{\Delta L}}} \right| < {k_0}} \end{array}}\\ {\begin{array}{*{20}{c}} 0&{\left| {\Delta {L_i}/{\sigma _{\Delta L}}} \right| \ge {k_0}} \end{array}} \end{array}} \right. $

(7)

式中，k₀一般取1.0或1.5，本文取k₀=1.5。

4) 以${\bar P_i^0}$为初始权，运用RELSM进行加权最小二乘估计，求得参数估计值${\mathit{\boldsymbol{\hat X}}}$、观测值残差V及方差因子σ：

$ \left\{ \begin{array}{l} \begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat X}} = {{({\mathit{\boldsymbol{A}}^{\rm{T}}}\bar P_i^0\mathit{\boldsymbol{A}})}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\bar P_i^0\mathit{\boldsymbol{L}}}\\ {\mathit{\boldsymbol{V}} = \mathit{\boldsymbol{A\hat X}} - \mathit{\boldsymbol{L}}} \end{array}\\ \sigma = {\rm{med(}}\left| \mathit{\boldsymbol{V}} \right|{\rm{)/0}}{\rm{.674}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{5}} \end{array} \right. $

(8)

5) 采用IGGⅢ权函数构建等价权：

$ {{\bar P}_i} = \left\{ \begin{array}{l} {P_i}, {{\bar V}_i} = \left| {\frac{{{V_i}}}{\sigma }} \right| \le {k_1}\\ {P_i}\frac{{{k_0}}}{{\left| {{{\bar V}_i}} \right|}}{\left( {\frac{{{k_1} - \left| {{{\bar V}_i}} \right|}}{{{k_2} - {k_{\rm{1}}}}}} \right)^2}, {k_1}{\rm{ < }}{{\bar V}_i} \le {k_2}\\ 0, {{\bar V}_i} > {k_2} \end{array} \right. $

(9)

式中，k₁=1.0~1.5，k₂=2.5~3，本文取k₁=1.5，k₂=3。将式(9)代入步骤4)，重新求解相关参数，重复步骤4)和5)，直至满足$\left| {{{\mathit{\boldsymbol{\hat X}}}^{k + 1}} - {{\mathit{\boldsymbol{\hat X}}}^k}} \right| < \delta $(δ为大于0的极小值，本文取δ=10^－10)时停止，将等价权为0所对应的测站剔除，再用剩余测站数据进行模型构建和参数求解。

3 算例及分析

本文数据来源于2009~2014年中国大陆构造环境监测网络在环渤海及其邻区的GPS对地观测速度场数据，共计482个监测台站，其中基准站32个、区域站450个，GPS站速率通过st_filter软件处理获得，相对于ITRF2008其精度优于2 mm/a。在34°~43°N、113°~126°E范围内截取环渤海区域共401个测站的观测数据。图 1为2009~2014年环渤海区域在ITRF2008框架下的速度场图像，该图像可反映研究区测站空间分布和地壳运动的整体趋势。

为验证本文方法及效果，设计3种方案对环渤海区域GPS速度场数据进行测站筛选及拟合分析。

1) 方案1：利用RRM进行抗差迭代计算，完成测站筛选后利用RRM进行拟合分析。

2) 方案2：利用RELSM进行抗差迭代计算，完成测站筛选后利用RELSM进行拟合分析。

3) 方案3：不经过测站筛选，对全部测站数据直接利用RELSM进行速度场拟合分析。

将速率拟合的残差标准差作为精度衡量指标对3种方案的拟合结果进行评价，计算公式为：

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

(10)

式中，ΔV_E和ΔV_N分别为模型解算的E向和N向速率残差值，t为模型的未知参数个数。表 1为3种方案解算的欧拉运动参数及速率残差标准差，图 2~4为3种方案的速率残差。为说明本文剔除异常测站的合理性，根据中国地震局提供的近20 a数据绘制环渤海区域主要地震分布图(图 5)，图 6为本文剔除的异常测站分布图。

从图 2~4可以看出，方案1和方案2计算的速率残差最大值在EW向和NS向均未超过5 mm/a和3 mm/a；而方案3除个别测站的最大速率残差达到10 mm/a，绝大多数测站的拟合残差在EW向和NS向分别集中在6 mm/a和5 mm/a以内，说明3种方案在分析和描述地壳长期运动趋势方面都有效。由表 1可以看出，方案1、方案2和方案3的拟合残差标准差分别为1.44 mm/a、1.39 mm/a和1.77 mm/a，表明即使采用RRM进行合理测站筛选后仍可得到优于RELSM的拟合结果；尽管方案3已考虑板内形变，采用RELSM对数据进行拟合分析，但未经过测站筛选，其结果最差；方案2既考虑板内形变，又采用高崩溃污染率抗差估计方法抑制粗差，同时进行严格的测站筛选，结果最优。

经统计，方案1检测出异常测站15个，方案2检测出异常测站24个，且几乎包含方案1中检测出的全部异常测站(仅1个测站不重合)。这表明方案2不仅可以顾及块体刚性运动的整体趋势，还可进一步检测出扣除局部构造形变(偏离线性应变)后的异常测站。从图 5、6可以看出，本文剔除的异常测站多分布在以郯庐和蓬莱为主的断裂带上以及近20 a来地震发生点附近，说明异常测站所在区域现今构造运动比较活跃，偏离块体刚性运动程度相对较大。从地壳运动的地球物理实际来看，方案2在顾及块体内部构造运动后，剔除异常测站进行速率拟合分析的方法更加合理。

4 结语

本文提出的顾及板内形变的抗差估计测站筛选及速度场拟合方法能够顾及地壳运动的物理实际，同时兼顾地壳运动的整体性和局部差异性，但由于RELSM假定块体内部形变呈线性变化，具有一定的局限性，在复杂地质构造区域可能难以适用。今后将着力研究地壳形变机制，构建更加精细和符合实际的地壳形变模型。