利用ITRF2008分析全球板块运动

李春晓

1 中国科学院云南天文台，昆明市官渡区，650216

MORVEL板块运动模型估计了全球25 个板块（占全球面积比为97.2%）的相对运动参数及其不确定度［1］。基于MORVEL 模型，Donald将剩余的31 个小板块（由Bird［2］定义，占全球面积比为2.8%）统一到MORVEL 模型中，从而构成NNR－MORVEL56［3］绝对板块运动模型。NNR－MORVEL56描述了全球56个板块相对于无旋转（no－net－rotation，NNR）参考架的绝对运动参数，每个板块的几何参数、运动参数及其不确定度可参考文献［1，3］。

ITRF2008［4］是对4种空间大测量技术（VLBI，SLR，GPS，DORIS）的观测数据进行重新解算获得的，考虑到它的解较以往国际地球参考架有了很大提高，本文以ITRF2008 参考架下225个GPS测站的速度场作为计算板块运动参数（欧拉矢量）的原始数据。

1 数学模型

若将地球近似为球体，那么总可以用欧拉定理将板块在地球表面的运动描述为整个板块围绕一个穿过地心的旋转轴作转动。定义CM（center of mass）为由激光测距解算的地球质心（包括地球、海洋、大气的质量）位置，CF（center of figure）为固态地球的形心位置［4－5］，板块的运动可用以下方程表述：

其中，ve、vn分别为测站水平速度的东向和北向分量，λ和φ分别为测站的大地测量经纬度，ωx、ωy、ωz分别为ω的3 个分量。用以下公式将求得的欧拉矢量ω用大地测量坐标（λ，φ）和模量ω表示为：

其中，φ′为地心纬度，f为地球扁率，f＝1／298.257 22。若要求解式（2）中的3 个未知量，则至少需要两个测站的速度场。通常情况下一个板块上会含有多个测站，从而式（2）为超定方程。利用加权最小二乘法求解超定方程可以得到板块的欧拉矢量，而欧拉矢量的估计误差可以通过误差传播定律来求得。假设已知测站在地心坐标系中的三维速度及其不确定度，那么可以利用以下方程得到测站的水平速度及其不确定度：

其中，cov（ven）为测站地平坐标系下水平速度的协方差矩阵，cov（vxyz）为ITRF2008参考架下测站三维速度的协方差矩阵，R为∂ven／∂vxyz，即从地心坐标系到测站地平坐标系的旋转矩阵，RT为转置矩阵。本文假定各个测站之间的水平速度场互不影响［7］，则由加权最小二乘法得到的拟合参数为：

其中，Ai为式（2）中的设计矩阵，表示矩阵的转置，vi为测站i的水平速度，n代表板块上总的测站数目，Wi为加权矩阵：

其中，rei为测站i东向速度残差，rni为测站i北向速度残差，n为板块上测站总数。拟合参数的协方差矩阵可由式（9）来计算：

欧拉矢量在球面坐标下的协方差矩阵及其误差椭圆也可以用误差传播定律来推导。首先，已知欧拉矢量的直角坐标（ωx，ωy，ωz）与球面坐标（ω，λ，φ′）之间有如下关系：

然后求出雅可比矩阵及其逆矩阵，也即

最后求出球面坐标下欧拉矢量的协方差矩阵：

对上面3×3协方差矩阵cov（ωωλφ′）提取出欧拉矢量经度和纬度之间的2×2 协方差矩阵cov（ωλφ′），然后利用误差传播定律推导出欧拉矢量的经度和大地测量纬度之间的协方差矩阵为：

误差椭圆的半长轴、半短轴以及长轴方向可通过对角化cov（ωλφ）求得，分别对应矩阵的两个特征值和较大特征值对应的特征向量。

2 测站筛选

测站的筛选主要遵循以下5个准则：1）测站的观测时间段不少于3a，主要用来避免季节的周期性对测站速度的影响；2）测站的速度场精度优于3mm／a，也即测站的东向和北向的速度不确定度都要小于3 mm／a；3）测站的第二不变应变率小于10－14／a，因为第二不变应变率的大小可以间接衡量板块内测站的形变大小，从而判定测站是否位于板块的刚性区域，这样避免了用测站是否位于板块边界地带作为判定依据的不可靠性；4）剔除受冰期回弹［8－10］的影响而产生附加速度较大的测站［6］，本文选取附加水平速度小于0.5mm／a作为限值［5］，也可以用附加垂直速度小于2.5mm／a作为限值；5）保留速度残差小于3mm／a的测站。

2.1 全球应变率模型

由于板块并非严格的刚体，当板块受到外界驱动（挤压、分离、滑动）时，其整体会发生形变，而且形变在空间上的分布也各不相同。板块上空间点的形变表现为应变。为了衡量全球应变的大小，相应的全球应变率模型相继建立。模型建立的基本方法如下［11］：1）获取全球测站的速度场，主要为GPS速度场，以及中亚地区第四纪断层滑动率和全球震源较浅（小于40km）的地震矩张量等数据；2）根据Kostrov理论将断层滑动率和地震矩张量转换为应变率［12］，并利用双三次Bessel插值［13］将旋转矢量函数（实际上该旋转矢量类似于欧拉矢量，不同的是欧拉矢量应用在整个板块上，为常矢量，而旋转矢量函数则为经度和纬度的函数）展开为球面坐标的函数，其中16个控制点作为未知参数；3）通过最小二乘法对基于模型和观测的速度场和应变率进行拟合，得到16个控制点，从而得到全球分布的旋转矢量函数；4）基于欧拉定理，利用得到的旋转矢量函数计算全球速度场，然后利用速度场计算速度场梯度张量场，根据速度场梯度张量场可以得到第二不变应变率。

本文使用的应变率模型为GSRMv1.2，该模型采用了全球4 281个GPS测站的5 170个水平速度场，第二不变应变率的分辨率为0.2°×0.2°，覆盖范围为68°S～80°N。图1为基于全球应变率模型GSRMv1.2 的第二不变应变率等值线图。从图中可以看出，板块边界地区的应变率比较大，同时有些板块的内部也存在较大形变，如青藏高原地区的应变图表现出欧亚板块内部存在一定的形变，说明板块并非严格的刚体。

图1 基于全球应变率模型GSRMv1.2的第二不变应变率等值线图Fig.1 Contour line for second invariant strain rate based on GSRMv1.2

通过对全球应变率模型［14－15］GSRMv1.2 中分辨率为0.2°×0.2°的网格节点第二不变应变率进行双三次二维平面插值得到测站的应变率，并根据测站筛选的第3 条准则选取了310 个GPS测站。

2.2 冰期回弹模型

冰期回弹模型包含冰模型和地球结构模型，冰模型主要包含冰在全球范围内的分布模型和冰的演化模型，地球结构模型主要包含地球的分层结构以及分层结构的物理参数，如密度、粘弹性等。目前最新的冰模型为ICE－5G，广泛采用的地球结构模型为VM2。本文所采用的冰模型为ICE－3G，地球模型为VM2。冰期回弹模型ICE－5G（VM2）预测的地壳垂直运动速率如图2（a）所示，图中显示冰期回弹对欧亚板块的北欧地区、北美板块和南极洲板块的部分区域影响最显著，垂直运动速率可达到1cm／a甚至更大，而且地壳上升的最大运动速率比地壳下降的运动速率要大。地壳水平运动速率东向分量和北向分量分别如图2（b）和2（c）所示，图中显示北美地区相对于其他地区有较大的地壳水平运动速率。对比不同地区的垂直和水平地壳运动速率看出，垂直运动速率一般比水平运动速率大4～5 倍，但是也存在例外，因为部分地区具有很小的垂直运动速率，但却有很大的水平运动速率，这可能与不同地区地壳对荷载的反应能力不同有关。

基于冰期回弹模型［16］ICE3G－VM2，用TABOO［17－18］软件对§2.1中选取的310个测站的垂直和水平位移速率进行估算，根据测站筛选的第4条准则从中选取242个测站。

图2 基于冰期回弹模型ICE－5G（VM2）的地壳垂直和水平运动速率的等值线图Fig.2 Contour line for rate of vertical and horizontal motion of solid earth based on ICE－5G（VM2）

3 计算结果与分析

利用上述4个准则进行筛选后，GPS测站剩余242个（图1中红色圆点所示）。通过将这些测站分派到全球56个板块中发现，只有16个板块含有GPS测站，其中扬子板块、加勒比板块和安纳托利亚板块各只有1 个GPS 测站，分别为WUHN、CRO1和ANKR。由于解算板块的运动参数至少需要2个测站，所以本文通过最小二乘法剔除速度残差较大（residual＞3σ）的测站并不断迭代，直到测站的数目和残差的均方根收敛为止，最后位于所有13个板块上的GPS测站数目为225个。

基于NNR－MORVEL56的板块边界划分情况，估计全球13个板块相对于ITRF2008的运动参数，建立了ITRF2008板块运动模型，这13个板块分别为阿穆尼亚板块、南极洲板块、阿拉伯板块、澳大利亚板块、欧亚板块、印度板块、北美板块、纳兹卡板块、努比亚板块、太平洋板块、南美板块、索马里板块和巽他古陆板块。表1给出了本文和Altamimi［6］计算的ITRF2008 参考架下全球13个板块的运动参数、中误差和水平速度残差的均方根，表2给出了ITRF2008板块运动模型和MORVEL板块运动模型的相对运动参数。

通过对比本文和Altamimi计算的结果看出，两者得到的欧拉矢量基本一致，但有几个板块在欧拉矢量的位置上达到6°的偏差：阿穆尼亚板块在沿纬度方向上存在6.95°的偏差；北美板块沿经度方向存在6.01°的偏差；巽他古陆板块在沿经度方向上存在6.44°的偏差。对比两者欧拉矢量的中误差可以看出，本文拟合的中误差比Altamimi拟合的中误差要小。尤其是对于巽他古陆板块，两者都选择了GETI和NTUS 两个GPS测站，但是拟合的结果却存在很大的差异。本文计算的欧拉矢量转速为0.346 4°／Ma±0.002 0°／Ma，欧拉极的经度和纬度为－89.207°±1.922°、50.880°±4.217°，中误差较小；Altamimi计算的欧拉极转速为0.388°／Ma±0.308°／Ma，欧拉极的经度和纬度为－87.309°±22.186°、44.436°±44.892°，中误差较大。对比东向速度和北向速度的均方根差可以看出，两者在数量级上保持一致，没有大的偏差。

导致两者之间个别板块存在较大差异的主要原因是模型建立的出发点不一样。Altamimi建立模型（参见式（1））的过程中考虑了原点偏差率，而且拟合方法为同时对所有板块进行拟合，以得到所有13个板块的欧拉矢量和一个原点偏差率，也即一次性拟合45个未知参数，那么各个板块之间的欧拉矢量必然存在联系，各个板块欧拉矢量之间的协方差系数不会近似为0，同时这些板块的欧拉矢量与原点偏差率之间也会存在一定的相关性。本文在建立模型的过程中，不考虑整体原点的偏差率，对13个板块的欧拉矢量分别进行拟合，每次拟合3个参数，消除了各个板块之间的相关性。由于板块之间的相关性，导致巽他古陆板块欧拉矢量的中误差很大，因为该板块仅存在两个测站，相对其他板块，这两个测站所占的权重较小，因而该板块的欧拉矢量受其他板块的影响较大。

测站的筛选准则也会对拟合结果造成影响。Altamimi采用以远离板块边界至少100km 作为判断测站是否位于刚性板块的判据之一，但是由于板块边界具有弥散性和选择100km 作为边界距离值具有较大的主观性，并且在精确计算测站到板块边界的过程中，需要计算点到球面多边形（将板块近似为球面多边形）的最短距离，运算量较大，故本文采用全球应变模型作为测站是否位于刚性板块的判据。造成偏差的另外一个原因是在最小二乘拟合的过程中，本文采用了3σ准则，即测站的东向和北向速度残差都要小于3σ，而Altamimi则采用的是3σ－3mm／a准则，即测站水平速度的残差既要小于3 倍水平速度的不确定度，也要小于3 mm／a。由于测站筛选的准则不同，每个板块上所筛选出的测站也会不同，最终拟合的欧拉矢量也会有偏差。

比较本文计算结果与MORVEL板块运动模型发现，努比亚－南极洲板块对的欧拉极在沿经度方向上的偏差比较大，达到15.8°；索马里－南极洲板块对的欧拉极在沿经度和纬度方向的偏差分别达到15.3°和13.5°；欧亚－北美板块对的欧拉极在沿经度和纬度方向的偏差分别达到10.3°和14.4°；努比亚－北美板块对的欧拉极在沿纬度方向上的偏差非常大，高达30.5°；印度－索马里板块对的欧拉极在沿纬度方向上的偏差为9.9°，在转速上的偏差较大，为0.118°／Ma。这些偏差的出现，可能说明努比亚、索马里、欧亚和印度板块在1～3 Ma 期间的运动确实发生了变化。由于ITRF2008板块运动模型与MORVEL 模型建立的时间尺度不一样，前者采用了近十几a的大地测量资料，为短期的结果，而后者则主要采用近几Ma的地质资料，为长期平均的结果，两者之间的偏差项可能意味着驱动板块运动的合力矩在零点发生了波动。另外，测站在空间上的几何分布不均会对计算结果造成一定的偏差，并且本文没有按照传统的测站筛选准则对测站进行筛选，而是以全球应变率模型为基准进行测站筛选，这也可能是造成偏差的原因，尤其是北美板块。

表1 ITRF2008 参考架下全球13个板块的运动参数Tab.1 Kinematic parameters for 13plates in ITRF2008

表2 ITRF2008 参考架下全球13个板块对的相对运动参数Tab.2 Kinematic parameters for 13pairs of plates in ITRF2008

4 结 语

本文采用基于全球应变模型来判断测站是否位于形变区内作为测站筛选的准则之一，把这种方法应用到ITRF2008板块运动模型中，并与Altamimi建立的板块运动模型和MORVEL 板块运动模型相比较，结果表明本文方法是可行的。地球表面的运动及变化本质上是由地球内部作用力驱动的，如果已知板块的运动参数，则从宏观上为研究地球内部的驱动力提供了方便。同时，以板块作为参考基准，可以为研究区域性地壳应变场和应力场提供借鉴。

［1］Demets C，Gordon R G，Argus D F.Geologically Current Plate Motions［J］.Geophysical Journal International，2010，181（1）：1－80

［2］Bird P.An Updated Digital Model of Plate Boundaries［J］.Geochemistry，Geophysics，Geosystems，2003，4（3）：101－112

［3］Argus D F，Gordon R G，Demets C.Geologically Current Motion of 56Plates Relative to the No－Net－Rotation Reference Frame［J］.Geochemistry Geophysics Geosystems，2011，12（11）：75－87

［4］Altamimi Z，Collilieux X，Métivier L.ITRF2008：An Improved Solution of the International Terrestrial Reference Frame［J］.Journal of Geodesy，2011，85（8）：457－473

［5］Argus D F，Gordon R G，Heflin M B，et al.The Angular Velocities of the Plates and the Velocity of Earth’s Centre from Space Geodesy［J］.Geophysical Journal International，2010，180（3）：913－960

［6］Altamimi Z，Métivier L，Collilieux X.ITRF2008 Plate Motion Model［J］.Journal of Geophysical Research：Solid Earth，2012，117（B7）：47－56

［7］Goudarzi M A，Cocard M，Santerre R.EPC：Matlab Software to Estimate Euler Pole Parameters［J］.GPS Solutions，2014，18（1）：153－162

［8］Peltier W R.Global Glacial Isostasy and the Surface of the Ice－Age Earth：The ICE－5G（VM2）Model and Grace［J］.Earth and Planetary Sciences，2004，32：111－149

［9］Argus D F，Peltier W R.Constraining Models of Postglacial Rebound Using Space Geodesy：A Detailed Assessment of Model ICE－5G（VM2）and Its Relatives［J］.Geophysical Journal International，2010，181（2）：697－723

［10］Klemann V，Martinec Z，Ivins E R.Glacial Isostasy and Plate Motion［J］.Journal of Geodynamics，2008，46（3－5）：95－103

［11］Haines A J，Holt W E.A Procedure for Obtaining the Complete Horizontal Motions within Zones of Distributed Deformation from the Inversion of Strain Rate Data［J］.Journal of Geophysical Research：Solid Earth，1993，98（B7）：12 057－12 082

［12］Holt W E，Chamot－Rooke N，Pichon X L，et al.Velocity Field in Asia Inferred from Quaternary Fault Slip Rates and Global Positioning System Observations［J］.Journal of Geophysical Research，2000，105（B8）：19 185－19 209

［13］David Salomon.Curves and Surfaces for Computer Graphics［M］.New York：Springer，2007

［14］Kreemer C，Haines J，Holt W E，et al.On the Determination of a Global Strain Rate Model［J］.Earth Planets ＆Space，2000，52（10）：765－770

［15］Kreemer C，Holt W E，Haines A J.An Integrated Global Model of Present－Day Plate Motions and Plate Boundary Deformation［J］.Geophysical Journal International，2003，154（1）：8－34

［16］Spada G，Galassi G.New Estimates of Secular Sea Level Risefrom Tide Gauge Data and GIA Modelling［J］.Geophysical Journal International，2012，191（3）：1 067－1 094

［17］Giorgio Spada.The Theory behind TABOO［EB／OL］.http：／／samizdat.mines.edu／taboo／teoria.pdf，2003

［18］Giorgio Spada.TABOO User Guide［EB／OL］.http：／／samizdat.mines.edu／taboo／user guide.pdf，2003