20世纪90年代以来，基于欧拉定理的全球板块刚性运动模型得到快速发展^[1]。这些模型假定块体在运动时其内部不发生形变，但实际上板块运动时其边缘具有发散特点和一定的弹塑性^[2]。李延兴等^[3-4]提出块体的弹性运动模型，主要包括整体旋转与均匀应变模型(REHSM)和整体旋转与线性应变模型(RELSM)，是对地壳刚性运动模型(RRM)的发展。但上述2种模型分别假设块体内部形变在空间的分布为均匀和线性变化，对于地质构造复杂的区域仍具有一定局限性。考虑到块体运动整体上一般以刚性运动为主，具有明显的整体运动趋势和长期稳定性特点，本文借鉴选权拟合思想^[5]，提出对整体旋转和线性应变模型中的欧拉参数附加不等式约束，将欧拉参数约束到由块体刚性运动模型所确定的长期运动的平均欧拉参数值附近，即建立附加不等式约束的地壳整体旋转和线性应变分析模型(IC-RELSM)，从而抑制隐伏断层、断裂带等局部不规则构造运动及异常测站对整个模型参数求解的影响。对中国大陆构造环境监测网络在环渤海区域的近几期GPS观测速度场数据进行拟合分析，结果表明，该方法具有一定的优越性。

1 附加不等式约束的整体旋转和线性应变模型构建 1.1 地壳运动的整体旋转与线性应变模型

传统的块体运动模型不考虑块体内部形变，将块体水平运动等效为绕欧拉轴的有限旋转运动，其模型为：

$ \left[ {\begin{array}{*{20}{c}} {{V_E}} \\ {{V_N}} \end{array}} \right]{\text{ = }}r\left[ {\begin{array}{*{20}{c}} {{\text{ - cos}}\lambda {\text{sin}}\varphi }&{{\text{ - sin}}\lambda {\text{sin}}\varphi }&{{\text{cos}}\varphi } \\ {{\text{sin}}\lambda }&{{\text{ - 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分别为板块任意一点在EW方向和NS方向的运动速率；r为地球平均半径；λ、φ为经纬度；ω_x、ω_y、ω_z为欧拉旋转参数。

由于地壳刚性运动模型会忽略块体内部的形变，与实际块体运动情况不符。一般情况下当块体受外力作用时，块体局部将发生一定程度的变形，某点的运动应由块体整体刚性旋转运动和局部弹塑性运动复合而成。若假设块体内部的形变具有均匀性的特点，则整体旋转与均匀应变模型可表示为^[3]：

$ \begin{gathered} \left[ {\begin{array}{*{20}{c}} {{V_E}} \\ {{V_N}} \end{array}} \right]{\text{ = }}r\left[ {\begin{array}{*{20}{c}} {{\text{ - cos}}\lambda {\text{sin}}\varphi }&{{\text{ - sin}}\lambda {\text{sin}}\varphi }&{{\text{cos}}\varphi } \\ {{\text{sin}}\lambda }&{{\text{ - cos}}\varphi }&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _x}} \\ {{\omega _y}} \\ {{\omega _z}} \end{array}} \right] + \hfill \\ {\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} {\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 _E}}&{{\varepsilon _{EN}}} \\ {{\varepsilon _{NE}}}&{{\varepsilon _N}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x \\ y \end{array}} \right] \hfill \\ \end{gathered} $

(2)

式中，ε_E、ε_N、ε_EN、ε_NE为应变参数；x、y为球面坐标系的坐标值，可表示为：

$ \left\{ \begin{array}{l} x = r(\lambda - {\lambda _0}){\rm{cos}}\varphi \\ y = r(\varphi - {\varphi _0}) \end{array} \right. $

(3)

式中，(λ₀, φ₀)为区域的几何中心经纬度，其他参数含义同式(1)。若进一步假设块体内部形变呈线性变化，则地壳运动的整体旋转与线性应变模型为^[4]：

$ \begin{gathered} {\left[ {\begin{array}{*{20}{c}} {{V_E}} \\ {{V_N}} \end{array}} \right]_h} = \left[ {\begin{array}{*{20}{c}} { - r{\text{cos}}\lambda {\text{sin}}\varphi }&{ - r{\text{sin}}\lambda {\text{sin}}\varphi }&{r{\text{cos}}\varphi } \\ {r{\text{sin}}\lambda }&{ - r{\text{cos}}\lambda }&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\omega _x}} \\ {{\omega _y}} \\ {{\omega _z}} \end{array}} \right] + \hfill \\ {\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}} {{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 \hfill \\ \end{gathered} $

(4)

式中，A₀~A₂、B₀~B₂、C₀~C₂均为应变参数。

1.2 附加不等式约束的整体旋转与线性应变模型(IC-RELSM)构建

附加不等式约束的模型一般可表示为^[6]：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{BX}} - \mathit{\boldsymbol{V}}\\ \mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{W}} \end{array} \right. $

(5)

式中，B为误差方程的系数矩阵；L为观测向量；V为观测向量的改正数；G为表达参数间约束关系的系数矩阵；W为约束条件中的常数向量。在最小二乘准则下，不等式约束平差问题可转化为：

$ \left\{ \begin{array}{l} \min f(\mathit{\boldsymbol{X}}) = {(\mathit{\boldsymbol{BX}} - \mathit{\boldsymbol{L}})^{\rm{T}}}\mathit{\boldsymbol{P}}(\mathit{\boldsymbol{BX}} - \mathit{\boldsymbol{L}})\\ \mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{W}} \end{array} \right. $

(6)

附加不等式约束的模型一般可通过迭代求解获得参数最佳估值，具体流程见文献[6]。

在建立模型解决实际问题时，应当根据具体问题对参数作具体分析，找出合理的权阵或参数约束矩阵，附加约束条件应尽可能符合实际^[5]。研究表明，区域地壳运动总体上主要遵循刚体运动规律，其在运动中表现出明显的整体运动趋势，这种趋势性运动一般具有长期稳定性特点，其定量形式可由欧拉运动矢量决定。为完整地刻画地壳运动，还必须考虑块体内部的不规则形变，某点的地壳运动应为趋势性刚体运动和内部不规则形变的复合形式。在建立区域地壳运动模型时，应充分考虑研究区所在块体整体运动趋势的长期稳定性(即应当考虑各点运动的平滑性)^[7]。为此，本文提出在研究区域地壳形变时，可借鉴选权拟合思想^[5]，仅对部分参数施加约束，即将模型参数中所包含的欧拉参数约束到由本区域长期积累的观测资料所确定的平均运动欧拉参数值附近，从而使形变参数既能顾及整体运动的平滑性，又能更真实地反映偏离块体整体运动趋势的精确性。

按式(5)进行地壳形变分析时，只需在式(4)中增加针对欧拉矢量的不等式约束模型，即可获得附加不等式约束的整体旋转与线性应变模型：

$ \left\{ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{V_E}}\\ {{V_N}} \end{array}} \right]{\rm{ = }}r\left[ {\begin{array}{*{20}{c}} {{\rm{ - cos}}\lambda {\rm{sin}}\varphi }&{{\rm{ - sin}}\lambda {\rm{sin}}\varphi }&{{\rm{cos}}\varphi }\\ {{\rm{sin}}\lambda }&{{\rm{ - cos}}\varphi }&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\\ \mathit{\boldsymbol{GW}} \le \mathit{\boldsymbol{H}} \end{array} \right. $

(7)

式中，W为研究区内各点在坐标轴上的3个欧拉旋转参数分量及9个应变参数组成的待求参数向量；G为约束矩阵，由于待求参数向量都为单一参数，此处G取单位阵；H为参照块体历史平均运动的欧拉参数值确定的区域地壳形变分析模型中针对欧拉参数不等式约束的临界向量。H取下界时，G和H相应部分皆作反符号处理。

2 环渤海区域地壳运动-形变分析

本文选用2009~2014年中国大陆构造环境监测网络在环渤海及其邻区的GPS对地观测速度场数据，共计482个监测台站，其中基准站32个，区域站450个。基准站采集数据的时间为2010~2014年，区域站采集数据的时间为2009年、2011年和2014年。通过对监测站观测资料的处理得到的环渤海及其邻域在ITRF2008框架下的统一速度场数据的精度优于2 mm/a。

根据文献[8]中环渤海区域大致范围，将区域外的测站点剔除，最终保留环渤海区域430个测点速度场数据。图 1(a)为2009~2014年环渤海区域相对于ITRF2008框架的速度场图像，可反映本区域地壳运动的整体趋势。鉴于研究区地壳运动主要以刚性运动为主，运动趋势具有一定的长期稳定性，研究区内各监测台站运动趋势大体一致，因此认为此区域在无偏估计下的旋转参数为其真值，即将李延兴等^[9]针对研究区解算的欧拉旋转参数作为具有先验信息的参数，将参数的约束范围扩大1/3得到约束值的上下界。

为验证本文方法，分别采用RELSM和IC-RELSM两种模型对环渤海区域的速度场数据进行拟合分析，结果见表 1，同时给出2种模型解算的速率残差(图 2)。

对本文结果进行如下分析：

1) 将文献[9]中欧拉运动参数作为对整体旋转与线性应变模型中相应欧拉参数附加不等式约束的先验参考值，考虑到区域地壳运动虽具有相对长期稳定的特征，但区域内部由于块体之间的相互挤压及各种地质构造运动的存在，使得块体内部形变相对于长期平均运动趋势仍会存在较大偏离。因此本文采用宽松约束，在施加约束时上下限均允许相对于李延兴等^[9]给出的欧拉参数存在1/3的偏差，即-2.050×10^-9≤ω_x≤-1.025×10^-9、-2.167×10^-9≤ω_y≤-1.083×10^-9、3.463×10^-9≤ω_z≤6.925×10^-9。从图 2(a)、2(b)可以看出，RELSM模型的速率拟合残差值在EW方向上绝大多数集中在5 mm/a以内，IC-RELSM模型则集中在4 mm/a以内，2种模型的拟合残差在NS方向整体情况相当，大多数均在3 mm/a以内。从表 1中均方根误差和残差标准差来看，RELSM模型均方根误差为0.673×10^-10 rad/a，残差标准差为1.753 mm/a，IC-RELSM模型分别为0.497×10^-10 rad/a和1.493 mm/a，IC-RELSM模型的结果明显优于RELSM模型。

2) 根据表 2中RELSM模型和IC-RELSM模型的应变参数，分别计算环渤海区域的主应变场参数，计算结果为：RELSM模型主压应变轴方向约为NE46.74°~110.82°，IC-RELSM模型主压应变轴方向约为NE54.85°~105.86°，与李延兴等^[9]研究渤海盆地时给出的主压应变轴的方向范围NE65.5°~89.3°相比，IC-RELSM模型结果与之更为一致。图 1(b)为由IC-RELSM模型计算的环渤海区域主应变场的空间分布，从图中可以看出，环渤海区域总体上呈现NEE-SWW方向压缩与NNW-SSE方向拉张的基本特征，这与地质研究结果大体一致。

3) 本文在构建GPS速度场IC-RELSM模型时未进行异常测站筛选，但从速度场拟合残差和残差标准差来看，仍可取得较好的结果，表明将整体旋转与线性应变模型中的欧拉参数约束到长期运动的平均欧拉参数值附近，可以在一定程度上抑制异常测站对速度场整体拟合结果的影响。上述现象的原因主要为IC-RELSM模型通过不等式约束可加强块体平均整体运动趋势对各点的平滑作用，异常测站的影响在迭代求解过程中被削弱，从而可提高拟合精度。

3 结语

本文顾及地壳运动的物理实际，将块体长期运动的平均欧拉参数值作为先验信息，建立针对RELSM模型中欧拉参数的不等式约束模型，从而得到更为合理、接近实际的速度场拟合模型IC-RELSM模型。IC-RELSM模型通过对欧拉参数附加不等式约束，强调块体内部各点运动的平滑性(与块体整体刚性运动趋势的一致性)和内部构造运动引起的不规则形变的平衡，即兼顾地壳运动的整体性和局部差异性。对环渤海区域GPS速度场的拟合结果表明，IC-RELSM模型不仅可保持RELSM模型的优良性，通过不等式约束还可以在一定程度上抑制异常测站对拟合结果的影响。但如何充分利用先验信息建立符合实际的不等式约束模型，不同区域可能存在差异。今后研究中将根据区域具体实际情况，进一步挖掘并合理利用更多的先验信息，构建符合客观实际的地壳形变分析模型。