抗差主成分估计在板块运动参数解算中的应用

李亚萍,孙付平,朱新慧,丁　赫,陈竞男

(1.信息工程大学 导航与空天目标工程学院,河南 郑州 450000；2.信息工程大学 地理空间信息学院，河南 郑州 450000)



抗差主成分估计在板块运动参数解算中的应用

李亚萍1,孙付平1,朱新慧1,丁赫1,陈竞男2

(1.信息工程大学 导航与空天目标工程学院,河南 郑州 450000；2.信息工程大学 地理空间信息学院，河南 郑州 450000)

摘要：为解决在采用空间大地测量数据求解板块运动欧拉参数过程中可能出现的矩阵病态问题，以及消除粗差对解算结果的影响，利用抗差主成分估计方法，采用IGG3方案的选权迭代方案对北美板块的2 577个台站进行检测，剔除包含粗差的异常站数据，求出欧拉参数最优解,并进行误差估计，建立北美板块的运动模型。计算结果表明，剔除异常台站前后精度有一定的提高，参数结果与其他模型的解算结果具有较强的一致性。

关键词：抗差主成分估计；欧拉参数；误差估计；北美板块

在板块运动机理研究中，通常假定板块是刚体，而事实上板块并非严格意义上的刚体，板块内局部地区和边界地区都存在着形变。因此，在建立板块运动模型时，必须选择板块稳定主体内的台站，才能完全反映板块的主体运动。许多研究基于最小二乘求解的，当观测误差服从正态分布且观测方程状态良好时，采用最小二乘估计具有较好的优越性，但对偏离主体分布的粗差不够敏感。若观测台站所处的环境发生变化，台站就有可能会出现异常，有粗差的存在。本文尝试以北美板块台站数据为实例利用抗差主成分估计的方法对台站进行筛选处理，以期剔除异常台站，选择更合理的台站解算欧拉参数的最优解。

1欧拉参数求解原理

通常用欧拉矢量参数反映板块运动的状态，根据欧拉定律，板块运动欧拉矢量Ω与测站地心速度V之间关系

(1)

式中:r为测站地心位置矢量。已知V和r可估计出Ω。地心速度与欧拉矢量的关系式：

(2)

板块运动主要反映在水平站速度之中，VLBI、SLR、GPS 3种技术确定垂向站速度受到的干扰因素还不能进行规范的建模，直接采用地心站速度V可能会使板块运动参数估计受到垂向站速度误差的影响。因此，在实用中常把测站地心速度转换为站心速度，采用其水平分量来求解板块运动欧拉矢量Ω。

1)测站地心速度转换为站心速度:

(3)

将式(2)代入式(3)得到站心坐标速度与欧拉矢量的关系式

(4)

(5)

2)精度评定。速度单位权中误差σ0

(6)

速度协方差矩阵Qvv,

(7)

式中:A′为剔除后的系数矩阵。

(8)

(9)

2抗差主成分估计模型及计算

2.1抗差主成分估计模型

设参数平差的函数模型为

(10)

式中：L为观测值向量，X为待求参数向量，A为X系数矩阵, Δ为误差向量。

X的LS估计为

(11)

设N=AΤPA；ϑ1,ϑ2,ϑ3,…,ϑt为N的t个特征值；λ1,λ2,…,λt所对应标准化特征向量。

为表述方便，假定λ1,λ2,…,λt依次降序排列。

记

(12)

式(10)改写为

(13)

如果矩阵A呈病态，假定有λl+1,λl+2,…,λt接近于零，并将Λ,Λ1,Λ2记为

(14)

若矩阵A呈病态，则式(12)中的Φ1为前l个标准化非零特征向量组成的矩阵，Φ2为λl+1,λl+2,…,λt对应的标准化特征向量。

则有定义推出主成分估计

(15)

根据抗差M估计原理，可得参数X的抗差估计为

(16)

可推导得到抗差主成分估计

(17)

(18)

2.2参数计算

建立误差方程

V=AX+L.

(19)

其中

(20)

整个迭代过程用流程图1表示。需要指出的是，初步处理的数据，在此次计算中作为原始数据。求解(AΤPA)特征根以及相应的特征向量，并按照第二种方法，尝试进行主成分个数的确定，进而确定Λ1，且根据实验精度要求以及以往的经验值给定值C=0.95，也即通常所说的置信度为95%。

图1　计算流程

3结果分析

3.1数据说明

本文采用的数据从美国国家大地测量局(NGS)网站下载获取，包括站点坐标以及速度矢量，为了充分更精确地研究板块的运动机制，还将布设在SLR、VLBI的台站数据进行统一处理，纳入到板块运动参数求解中，分布展点见图2。为避免

图2　台站展点图

原始数据存在的较明显的粗差带来不必要的干扰，先采用文献[1]中提到的3个原则进行初步筛选。对数据进行历元统一，统一到历元2005。经过初步处理，选取2 557个站参与抗差估计计算。

3.2结果分析

从计算过程可知，法方程系数矩阵的条件数N=256，可知存在中等程度病态问题。在迭代过程中经过多次试验，发现降权中的k1=3.0更为合适，因为在满足筛选台站条件下，尽可能多的台站有利于提高解算精度。水平速度残差百分比见图3，剔除异常台站前后水平速度均方根误差。欧拉参数及其中误差分别见表1、表2。

图3　水平速度残差百分比

对欧拉参数进行检验，选择其它模型进行比较，结果如表3所示。

表1　剔除异常台站前后水平速度均方根误差(RMSE)

表2　欧拉参数中误差

表3　几种模型的比较

注：Ⅰ.NNR-NUVEL1A[10]估计结果；

Ⅱ.朱文耀(2000)ITRF96[11]估计结果；

Ⅲ.柴洪洲 ITRF2005[12]估计结果；

Ⅳ.朱新慧 ITRF2008VEL[13]估计结果;

V.本文岭估计结果

由表1可知，剔除异常台站前，台站的水平方向的速度东向和北向的均方根误差(E-RMSE 、N-RMSE)分别达到12.78,20.46,而利用抗差模型剔除异常站后,二者分别减小到0.46,0.73,说明剔除异常站后，说明利用抗差主成分模型计算的精度明显提高。

表2给出了利用抗差方法后求得的欧拉参数的中误差，结果表明其精度达到10-7量级，说明求解的参数具有较好的精度，可靠性强。

水平速度残差在一定程度上反映数据的质量，由图3可知，残差在1 mm以内的百分比均能达到80%以上，说明数据精度质量很好。由表3可以看出，本文求得的结果与其他学者的结果基本上在同一个量级，但是也有一定的差别，原因是其他学者的研究大多基于地学资料、ITRF系列的研究，采用的数据量相对本文较少，台站的选择不同，求解的结果也会有不同。这些结果相差在一定的范围内，说明板块的整体运动在一定时期内是相对稳定的。

4结束语

利用抗差主成分估计方法，采用IGG3方案的选权迭代方案对北美板块的2 577个台站数据进行检测，剔除包含粗差的异常站数据，求出欧拉参数的最优解,并进行误差估计，建立北美板块的运动模型；计算结果表明，台站的水平方向的速度东向和北向的均方根误差明显降低，残差在1 mm以内能达到80%以上，说明抗差的效果较好，求解的参数具有较好的精度，可靠性强。不足之处，仅对北美板块的运动参数进行解算与分析，需利用该方法进行选站，以期解算全球范围内的板块参数，进一步利用最优欧拉参数进行理论速度的计算，进而可以与实测速度进行比较，分析各个块体内部变形情况。

参考文献：

[1]金双根.确定板块运动学模型的台站选取[J].大地测量与地球动力学，2003，23(3)：56-60.

[2]黄立人.用于相对稳定点组判别的QUAD法[J].大地测量与地球动力学，2002，22(2)：10-15.

[3]徐克科.利用抗差估计进行刚性板块异常台站探测及其运动参数求取[J].大地测量与地球动力学，2014，34(2)：95-99.

[4]柴洪洲.最小二乘配置方法确定中国大陆主要块体运动模型[J].测绘学报，2009，58(1)：61-65.

[5]明锋，柴洪洲.利用半参数模型求解板块运动参数[J].大地测量与地球动力学，2008，28(3)：104-107.

[6]隋立芬，张超.抗差主成分估计及其应用[J].测绘学院学报，1996，13(3)：195-199.

[7]杨元喜.抗差估计理论及其在大地测量中的应用[D].北京：中国科学院，1991.

[8]周江文.经典误差理论与抗差估计[J].测绘学报，1989，18(2)，115-120.

[9]黄维彬.近代平差理论及其应用[M].北京：解放军出版社，1992.

[10] ARGUS D F,GORDON R G.No-Net-Rotation Model of Current Plate Velocities Incorporate Motion Model NUVEL-1[J].Geophys.Res.Lett.,1991,18(11):2039-2042.

[11] 朱文耀.基于ITRF96和ITRF97板块运动模型[J].天文学报，2000，41(3)：312-319.

[12] 柴洪洲.全球板块运动背景场及ITRF2005VEL的建立[J].测绘科学技术学报，2009，26 (2)：79-85.

[13] 朱新慧，孙付平，王刃.运用空间技术建立绝对板块运动模型[J].武汉大学学报(信息科学版)，2012，37(3)：282-285.

[14] 张洪文，赵忠海，朱李忠，等.利用GPS基线长度研究南极板块稳定性[J].测绘与空间地理信息，2015,38(5)：18-20.

[责任编辑：李铭娜]

DOI:10.19349/j.cnki.issn1006-7949.2016.08.009

收稿日期：2015-05-11

基金项目：国家自然科学基金资助项目(41374027)

作者简介：李亚萍(1989-)，女，硕士研究生.

中图分类号：P226.3

文献标识码：A

文章编号：1006-7949(2016)08-0038-04

Application of principal components estimate to the calculation of plate motion parameters

LI Yaping1，SUN Fuping1，ZHU Xinhui1，DING He1，CHEN Jingnan2

(1.School of Navigation and Aerospace Engineering,Information Engineering University,Zhengzhou 450000,China；2.School of Geographic Spatial Information,Information Engineering University,Zhengzhou 450000,China)

Abstract：To solve the problem of motion matrix morbid appearing in the process of solving the Euler parameters of plate motion using the space geodetic data,and to eliminate the influence of gross error on the results,this paper detects the 2 577 stations of North American plate by using principal component estimation method with IGG3 iteration scheme.The abnormal station data containing gross error is eliminated.Then the North American plate motion model is established with the optimal solution of Euler parameters and error estimation.The result shows that the accuracy increases to a certain extent before and after excluding abnormal stations and the results have strong consistency with other models.

Key words：principal component estimate；Euler parameter；error estimate;North American plate