李亞萍,孫付平,朱新慧,丁　赫,陳競(jìng)男

(1.信息工程大學(xué) 導(dǎo)航與空天目標(biāo)工程學(xué)院,河南 鄭州 450000；2.信息工程大學(xué) 地理空間信息學(xué)院，河南 鄭州 450000)

?

抗差主成分估計(jì)在板塊運(yùn)動(dòng)參數(shù)解算中的應(yīng)用

李亞萍1,孫付平1,朱新慧1,丁赫1,陳競(jìng)男2

(1.信息工程大學(xué) 導(dǎo)航與空天目標(biāo)工程學(xué)院,河南 鄭州 450000；2.信息工程大學(xué) 地理空間信息學(xué)院，河南 鄭州 450000)

摘要：為解決在采用空間大地測(cè)量數(shù)據(jù)求解板塊運(yùn)動(dòng)歐拉參數(shù)過程中可能出現(xiàn)的矩陣病態(tài)問題，以及消除粗差對(duì)解算結(jié)果的影響，利用抗差主成分估計(jì)方法，采用IGG3方案的選權(quán)迭代方案對(duì)北美板塊的2 577個(gè)臺(tái)站進(jìn)行檢測(cè)，剔除包含粗差的異常站數(shù)據(jù)，求出歐拉參數(shù)最優(yōu)解,并進(jìn)行誤差估計(jì)，建立北美板塊的運(yùn)動(dòng)模型。計(jì)算結(jié)果表明，剔除異常臺(tái)站前后精度有一定的提高，參數(shù)結(jié)果與其他模型的解算結(jié)果具有較強(qiáng)的一致性。

關(guān)鍵詞：抗差主成分估計(jì)；歐拉參數(shù)；誤差估計(jì)；北美板塊

在板塊運(yùn)動(dòng)機(jī)理研究中，通常假定板塊是剛體，而事實(shí)上板塊并非嚴(yán)格意義上的剛體，板塊內(nèi)局部地區(qū)和邊界地區(qū)都存在著形變。因此，在建立板塊運(yùn)動(dòng)模型時(shí)，必須選擇板塊穩(wěn)定主體內(nèi)的臺(tái)站，才能完全反映板塊的主體運(yùn)動(dòng)。許多研究基于最小二乘求解的，當(dāng)觀測(cè)誤差服從正態(tài)分布且觀測(cè)方程狀態(tài)良好時(shí)，采用最小二乘估計(jì)具有較好的優(yōu)越性，但對(duì)偏離主體分布的粗差不夠敏感。若觀測(cè)臺(tái)站所處的環(huán)境發(fā)生變化，臺(tái)站就有可能會(huì)出現(xiàn)異常，有粗差的存在。本文嘗試以北美板塊臺(tái)站數(shù)據(jù)為實(shí)例利用抗差主成分估計(jì)的方法對(duì)臺(tái)站進(jìn)行篩選處理，以期剔除異常臺(tái)站，選擇更合理的臺(tái)站解算歐拉參數(shù)的最優(yōu)解。

1歐拉參數(shù)求解原理

通常用歐拉矢量參數(shù)反映板塊運(yùn)動(dòng)的狀態(tài)，根據(jù)歐拉定律，板塊運(yùn)動(dòng)歐拉矢量Ω與測(cè)站地心速度V之間關(guān)系

(1)

式中:r為測(cè)站地心位置矢量。已知V和r可估計(jì)出Ω。地心速度與歐拉矢量的關(guān)系式：

(2)

板塊運(yùn)動(dòng)主要反映在水平站速度之中，VLBI、SLR、GPS 3種技術(shù)確定垂向站速度受到的干擾因素還不能進(jìn)行規(guī)范的建模，直接采用地心站速度V可能會(huì)使板塊運(yùn)動(dòng)參數(shù)估計(jì)受到垂向站速度誤差的影響。因此，在實(shí)用中常把測(cè)站地心速度轉(zhuǎn)換為站心速度，采用其水平分量來求解板塊運(yùn)動(dòng)歐拉矢量Ω。

1)測(cè)站地心速度轉(zhuǎn)換為站心速度:

(3)

將式(2)代入式(3)得到站心坐標(biāo)速度與歐拉矢量的關(guān)系式

(4)

(5)

2)精度評(píng)定。速度單位權(quán)中誤差σ0

(6)

速度協(xié)方差矩陣Qvv,

(7)

式中:A′為剔除后的系數(shù)矩陣。

(8)

(9)

2抗差主成分估計(jì)模型及計(jì)算

2.1抗差主成分估計(jì)模型

設(shè)參數(shù)平差的函數(shù)模型為

(10)

式中：L為觀測(cè)值向量，X為待求參數(shù)向量，A為X系數(shù)矩陣, Δ為誤差向量。

X的LS估計(jì)為

(11)

設(shè)N=AΤPA；?1,?2,?3,…,?t為N的t個(gè)特征值；λ1,λ2,…,λt所對(duì)應(yīng)標(biāo)準(zhǔn)化特征向量。

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

記

(12)

式(10)改寫為

(13)

如果矩陣A呈病態(tài)，假定有λl+1,λl+2,…,λt接近于零，并將Λ,Λ1,Λ2記為

(14)

若矩陣A呈病態(tài)，則式(12)中的Φ1為前l(fā)個(gè)標(biāo)準(zhǔn)化非零特征向量組成的矩陣，Φ2為λl+1,λl+2,…,λt對(duì)應(yīng)的標(biāo)準(zhǔn)化特征向量。

則有定義推出主成分估計(jì)

(15)

根據(jù)抗差M估計(jì)原理，可得參數(shù)X的抗差估計(jì)為

(16)

可推導(dǎo)得到抗差主成分估計(jì)

(17)

(18)

2.2參數(shù)計(jì)算

建立誤差方程

V=AX+L.

(19)

其中

(20)

整個(gè)迭代過程用流程圖1表示。需要指出的是，初步處理的數(shù)據(jù)，在此次計(jì)算中作為原始數(shù)據(jù)。求解(AΤPA)特征根以及相應(yīng)的特征向量，并按照第二種方法，嘗試進(jìn)行主成分個(gè)數(shù)的確定，進(jìn)而確定Λ1，且根據(jù)實(shí)驗(yàn)精度要求以及以往的經(jīng)驗(yàn)值給定值C=0.95，也即通常所說的置信度為95%。

圖1　計(jì)算流程

3結(jié)果分析

3.1數(shù)據(jù)說明

本文采用的數(shù)據(jù)從美國(guó)國(guó)家大地測(cè)量局(NGS)網(wǎng)站下載獲取，包括站點(diǎn)坐標(biāo)以及速度矢量，為了充分更精確地研究板塊的運(yùn)動(dòng)機(jī)制，還將布設(shè)在SLR、VLBI的臺(tái)站數(shù)據(jù)進(jìn)行統(tǒng)一處理，納入到板塊運(yùn)動(dòng)參數(shù)求解中，分布展點(diǎn)見圖2。為避免

圖2　臺(tái)站展點(diǎn)圖

原始數(shù)據(jù)存在的較明顯的粗差帶來不必要的干擾，先采用文獻(xiàn)[1]中提到的3個(gè)原則進(jìn)行初步篩選。對(duì)數(shù)據(jù)進(jìn)行歷元統(tǒng)一，統(tǒng)一到歷元2005。經(jīng)過初步處理，選取2 557個(gè)站參與抗差估計(jì)計(jì)算。

3.2結(jié)果分析

從計(jì)算過程可知，法方程系數(shù)矩陣的條件數(shù)N=256，可知存在中等程度病態(tài)問題。在迭代過程中經(jīng)過多次試驗(yàn)，發(fā)現(xiàn)降權(quán)中的k1=3.0更為合適，因?yàn)樵跐M足篩選臺(tái)站條件下，盡可能多的臺(tái)站有利于提高解算精度。水平速度殘差百分比見圖3，剔除異常臺(tái)站前后水平速度均方根誤差。歐拉參數(shù)及其中誤差分別見表1、表2。

圖3　水平速度殘差百分比

對(duì)歐拉參數(shù)進(jìn)行檢驗(yàn)，選擇其它模型進(jìn)行比較，結(jié)果如表3所示。

表1　剔除異常臺(tái)站前后水平速度均方根誤差(RMSE)

表2　歐拉參數(shù)中誤差

表3　幾種模型的比較

注：Ⅰ.NNR-NUVEL1A[10]估計(jì)結(jié)果；

Ⅱ.朱文耀(2000)ITRF96[11]估計(jì)結(jié)果；

Ⅲ.柴洪洲 ITRF2005[12]估計(jì)結(jié)果；

Ⅳ.朱新慧 ITRF2008VEL[13]估計(jì)結(jié)果;

V.本文嶺估計(jì)結(jié)果

由表1可知，剔除異常臺(tái)站前，臺(tái)站的水平方向的速度東向和北向的均方根誤差(E-RMSE 、N-RMSE)分別達(dá)到12.78,20.46,而利用抗差模型剔除異常站后,二者分別減小到0.46,0.73,說明剔除異常站后，說明利用抗差主成分模型計(jì)算的精度明顯提高。

表2給出了利用抗差方法后求得的歐拉參數(shù)的中誤差，結(jié)果表明其精度達(dá)到10-7量級(jí)，說明求解的參數(shù)具有較好的精度，可靠性強(qiáng)。

水平速度殘差在一定程度上反映數(shù)據(jù)的質(zhì)量，由圖3可知，殘差在1 mm以內(nèi)的百分比均能達(dá)到80%以上，說明數(shù)據(jù)精度質(zhì)量很好。由表3可以看出，本文求得的結(jié)果與其他學(xué)者的結(jié)果基本上在同一個(gè)量級(jí)，但是也有一定的差別，原因是其他學(xué)者的研究大多基于地學(xué)資料、ITRF系列的研究，采用的數(shù)據(jù)量相對(duì)本文較少，臺(tái)站的選擇不同，求解的結(jié)果也會(huì)有不同。這些結(jié)果相差在一定的范圍內(nèi)，說明板塊的整體運(yùn)動(dòng)在一定時(shí)期內(nèi)是相對(duì)穩(wěn)定的。

4結(jié)束語

利用抗差主成分估計(jì)方法，采用IGG3方案的選權(quán)迭代方案對(duì)北美板塊的2 577個(gè)臺(tái)站數(shù)據(jù)進(jìn)行檢測(cè)，剔除包含粗差的異常站數(shù)據(jù)，求出歐拉參數(shù)的最優(yōu)解,并進(jìn)行誤差估計(jì)，建立北美板塊的運(yùn)動(dòng)模型；計(jì)算結(jié)果表明，臺(tái)站的水平方向的速度東向和北向的均方根誤差明顯降低，殘差在1 mm以內(nèi)能達(dá)到80%以上，說明抗差的效果較好，求解的參數(shù)具有較好的精度，可靠性強(qiáng)。不足之處，僅對(duì)北美板塊的運(yùn)動(dòng)參數(shù)進(jìn)行解算與分析，需利用該方法進(jìn)行選站，以期解算全球范圍內(nèi)的板塊參數(shù)，進(jìn)一步利用最優(yōu)歐拉參數(shù)進(jìn)行理論速度的計(jì)算，進(jìn)而可以與實(shí)測(cè)速度進(jìn)行比較，分析各個(gè)塊體內(nèi)部變形情況。

參考文獻(xiàn)：

[1]金雙根.確定板塊運(yùn)動(dòng)學(xué)模型的臺(tái)站選取[J].大地測(cè)量與地球動(dòng)力學(xué)，2003，23(3)：56-60.

[2]黃立人.用于相對(duì)穩(wěn)定點(diǎn)組判別的QUAD法[J].大地測(cè)量與地球動(dòng)力學(xué)，2002，22(2)：10-15.

[3]徐克科.利用抗差估計(jì)進(jìn)行剛性板塊異常臺(tái)站探測(cè)及其運(yùn)動(dòng)參數(shù)求取[J].大地測(cè)量與地球動(dòng)力學(xué)，2014，34(2)：95-99.

[4]柴洪洲.最小二乘配置方法確定中國(guó)大陸主要塊體運(yùn)動(dòng)模型[J].測(cè)繪學(xué)報(bào)，2009，58(1)：61-65.

[5]明鋒，柴洪洲.利用半?yún)?shù)模型求解板塊運(yùn)動(dòng)參數(shù)[J].大地測(cè)量與地球動(dòng)力學(xué)，2008，28(3)：104-107.

[6]隋立芬，張超.抗差主成分估計(jì)及其應(yīng)用[J].測(cè)繪學(xué)院學(xué)報(bào)，1996，13(3)：195-199.

[7]楊元喜.抗差估計(jì)理論及其在大地測(cè)量中的應(yīng)用[D].北京：中國(guó)科學(xué)院，1991.

[8]周江文.經(jīng)典誤差理論與抗差估計(jì)[J].測(cè)繪學(xué)報(bào)，1989，18(2)，115-120.

[9]黃維彬.近代平差理論及其應(yīng)用[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板塊運(yùn)動(dòng)模型[J].天文學(xué)報(bào)，2000，41(3)：312-319.

[12] 柴洪洲.全球板塊運(yùn)動(dòng)背景場(chǎng)及ITRF2005VEL的建立[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2009，26 (2)：79-85.

[13] 朱新慧，孫付平，王刃.運(yùn)用空間技術(shù)建立絕對(duì)板塊運(yùn)動(dòng)模型[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2012，37(3)：282-285.

[14] 張洪文，趙忠海，朱李忠，等.利用GPS基線長(zhǎng)度研究南極板塊穩(wěn)定性[J].測(cè)繪與空間地理信息，2015,38(5)：18-20.

[責(zé)任編輯：李銘娜]

DOI:10.19349/j.cnki.issn1006-7949.2016.08.009

收稿日期：2015-05-11

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(41374027)

作者簡(jiǎn)介：李亞萍(1989-)，女，碩士研究生.

中圖分類號(hào)：P226.3

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：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