基于速度場模型的區(qū)域CGCS2000坐標實時轉換

嚴劍鋒 王英剛

(南京市不動產登記中心, 江蘇 南京 210005)

0 引言

將WGS84(World Geodetic System 1984)坐標系(實時歷元和最新框架)下的定位結果轉換到2000國家大地坐標系(CGCS2000)系統(tǒng)下,需要經過將測量結果進行坐標框架和歷元轉換兩個步驟。坐標框架間的轉換參數(shù),在ITRF官方網站可以直接獲取[1],而歷元間的轉換,則需要站點速度信息。因此建立高精度的速度場模型對實現(xiàn)全球導航衛(wèi)星系統(tǒng)(Global Navigation Satellite System，GNSS)測量成果向CGCS2000的轉換和各時期測量成果的統(tǒng)一有著重要的作用。

為獲得高精度、實時動態(tài)的速度場模型,國內很多學者先后利用中國地殼運動觀測網絡工程數(shù)據(jù)開展了關于中國地殼運動速度場的研究。劉經南等對我國多個高精度全球定位系統(tǒng)(Global Positioning System，GPS)區(qū)域網進行了整體解算,得到了統(tǒng)一參考框架下的中國地殼運動水平速度場[2]。楊元喜等利用最小二乘配置法將歐拉矢量模型作為地殼運動的趨勢項模型,提高了中國大陸地殼水平運動場的計算精度[3]。江在森采用最小二乘配置法建立了中國大陸地殼運動模型,并獲取了水平應變場結果[4-5]。武艷強等在利用多面函數(shù)擬合方法建立速度場模型的基礎上首次將其應用于大區(qū)域應變場的解算[6]。蔣志浩[7]等采用有限元插值方法建立了中國大陸區(qū)域地殼運動速度場模型。本文基于美國大地測量所(National Geodetic Survey,NGS)196個連續(xù)運行參考站(Continuously Operating Reference Stations，CORS)近8年的觀測數(shù)據(jù),經過基線解算和時間序列分析獲取CORS站點的點位速度,然后分析對比了幾種建立速度場模型的方法,最后基于最優(yōu)速度場模型和坐標框架轉換理論,建立融合速度場改正信息的七參數(shù)轉換模型,并應用于WGS84與CGCS2000坐標轉換。程鵬飛,成英燕[8]以板塊運動模型劃分板塊邊界為基礎,在同一個板塊上,利用監(jiān)督分類七參數(shù)法選取基準站,并采用速度場歸算的方法進行高精度的GNSS測量坐標歸算到CGCS2000坐標。

1 速度場模型建立

常用的地殼運動速度場模型的建立方法包括歐拉矢量法和擬合法。歐拉矢量法的特點是把整個建模區(qū)域當成剛性塊體,認為板塊運動是整體性的,內部形變一致;擬合方法是基于空間大地測量數(shù)據(jù),通過數(shù)學擬合建立速度場模型,對研究區(qū)域沒有剛彈性的要求。

1.1 局部歐拉矢量法

局部歐拉矢量法就是應用歐拉矢量基本原理,把待求點周圍局部范圍看成是剛體,并建立模型來獲取待求點的速度。

1.1.1歐拉矢量基本原理

假定一個剛體,繞某一固定點旋轉,那么這個剛體中所有點位的偏移可以看作是該剛體繞固定點的固定軸旋轉的結果。因此,如果把地心當作固定點,研究區(qū)域當作剛體,那么根據(jù)歐拉定理研究區(qū)域內某一點的點位變化速度可表示為:

(1)

式(1)中的歐拉矢量可以根據(jù)研究區(qū)域3個或3個以上速度信息已知的點位求得。如果超過3個已知點,可以根據(jù)最小二乘準則求得歐拉矢量最佳估值[8]。

1.1.2局部歐拉矢量法

以歐拉矢量法基本原理作為基礎,建立局部歐拉矢量模型的關鍵在于如何選擇待求點周圍區(qū)域的已知速度點,本文采用移動窗口法進行點位選取。

移動窗口法一般應遵循的原則為:在盡量小的區(qū)域內選擇已知點來建立模型,而且選擇的已知點應盡量均勻分布在待求點的四周,且數(shù)量在4個以上。但是,有的待求點在區(qū)域的邊界或周圍點位比較稀疏,很難滿足以上要求,此時應擴大搜索范圍或用其他方式獲取待求點速度信息[9]。

1.2 數(shù)學擬合法

1.2.1反距離加權法

反距離加權法是一種基于鄰近的多個點的加權平均值來估計未知點的一種插值方法,待插值點的數(shù)值和它到已知點的距離呈反相關,即待求點離已知點的距離越近,它的數(shù)值就越接近已知點[10]。計算待插值點的速度公式為:

(2)

式中，V(x0)為待插值點x0的速度;V(xi)為建模點xi的速度參考值;n為所需建模點個數(shù);d(xi-x0)為建模點到待插值點之間的距離;D(xi)是由d(xi-x0)確定的權系數(shù),k為距離d(xi-x0)的冪,對插值結果影響較大,一般取k=2。如果待插值點周圍的建模點數(shù)少且速度變化快時,插值精度不高。

1.2.2多面函數(shù)擬合法

多面函數(shù)擬合的理論基礎為:任一光滑數(shù)學表面,可以用一系列規(guī)則數(shù)學表面逼近。式(3)給出了該方法的數(shù)學模型:

(3)

式中,V(x,y)為點(x,y)處的點位速度值;Q(x,y;xi,yi)為和函數(shù),一般形式為Q(x,y;xi,yi)=((x-xi)2+(y-yi)2+(z-zi)2+σ2)β;σ為光滑因子,取值一般在0.001～1.00;β為指數(shù)參數(shù),一般取0.5、1、1.5或-0.5等;u為所取節(jié)點個數(shù);αi(i=1,2,…,u)為待估參數(shù)。根據(jù)公式(3),利用間接平差理論可以方便地求得待定參數(shù)αi,并且進行精度評定[11]。

1.2.3最小二乘配置法

地殼運動的點位坐標速度可描述成兩部分:板塊長期整體運動和短期局部運動[12]。最小二乘配置法既能夠估計非隨機參數(shù)(長期的運動趨勢),又能夠估計隨機信號和誤差(局部的小范圍運動),如果將地殼運動整體特征的歐拉矢量模型作為非隨機參數(shù)項,將局部小范圍運動信息作為隨機信號項,根據(jù)最小二乘配置原理,點位坐標速度的觀測函數(shù)模型為[13]:

Vn×1=An×3Ω3×1+Un×mSm×1+nq×1

(4)

式中,V為站點三維速度向量;Ω為歐拉參數(shù)向量;A為其系數(shù)矩陣(由點位坐標信息組成)；AΩ表示板塊運動中整體的、剛性的變化部分;US表示從點位速度V中扣除某種剛性運動后剩余部分中的有效信號部分的剩余速度。按最小二乘配置模型,式(4)中非隨機參數(shù)Ω和信號S的解為:

Ω=(ATC-1A)-1ATC-1V

S=CuoC-1(V-AΩ)

(5)

其中,

C=Coo+Cnn

(6)

式中,Coo為站點速度向量V的協(xié)方差；Cnn為V觀測誤差自協(xié)方差；Cuo為推估(點)信號與已測(點)信號的協(xié)方差矩陣。

上面共用到三個協(xié)方差陣Coo、Cnn和Cuo,其中Cnn為已知,由點位速度解算結果給出,而Coo和Cuo是未知的,計算公式如下:

F(d)=C(0)exp(-k2d2)

(7)

式(7)為高斯型函數(shù)。式中，d為觀測點間距離；C(0)表示d=0的方差。假設空間上兩個站點的速度值與兩站間的距離呈反相關,即距離越近,點位速度相關性越高,則高斯型函數(shù)能夠較好地表達不同站點間的協(xié)方差F(d)隨著點間距離d的增大而衰減。由此可以根據(jù)已有的點位速度數(shù)據(jù)建立協(xié)方差函數(shù),從而構建協(xié)方差陣的Coo和Cuo的值。

1.3 速度場模型精度對比

根據(jù)前面提到的速度場模型建立方法,編程實現(xiàn)模型并進行精度評定,最終選擇最優(yōu)速度場模型參與CGCS2000坐標轉換。

基于NGS的196個建模點,分別利用局部歐拉矢量法、反距離加權法、多面函數(shù)法和最小二乘配置法建立速度場模型,然后選擇56個點位速度已知點對模型進行精度評定,圖1給出了幾種方法得到的點位速度誤差曲線圖。

不同速度場模型計算得到外符合精度如表1所示。

表1 幾種速度場模型外符合精度 單位:mm

由圖1可知,幾種速度場模型得到的待求點誤差都在毫米級,從誤差分布來看,局部歐拉矢量和最小二乘配置法計算誤差大都在3 mm以內,其他兩種方法個別點誤差較大。從表1可知,最小二乘配置法外符合精度最高,各個方向都在1 mm左右,且X方向速度值誤差僅為0.69 mm。

2 CGCS2000坐標轉換

2.1 融合速度場改正信息的七參數(shù)實時轉換模型

實時坐標轉換參數(shù)的求取是基于速度場模型、坐標框架轉換理論實現(xiàn)的。一般情況下,CGCS2000坐標實時轉換是結合CNSS測量手簿中的七參數(shù)轉換模塊來實現(xiàn)的,因此實現(xiàn)實時轉換的關鍵部分是求取融合速度場的實時轉換參數(shù)。

圖1 速度場模型點位誤差圖

2.1.1 融合速度場信息的實時轉換參數(shù)求取

(8)

同理,若某個點位的變化速度為V,某坐標系統(tǒng)定義的歷元為t0,那么實時歷元下獲取的坐標相對于t1歷元時的坐標的位置變化值S可表示為:

S=(t1-t0)×V

(9)

2.1.2CGCG2000坐標實時轉換

WGS84向CGCG2000坐標轉換即ITRF2008框架下當前歷元的坐標向目標框架ITRF97下的2 000.0歷元坐標進行轉換。假設測站點的空間點位速度為ΔX、ΔY、ΔZ,單位為mm/年,轉換公式可表示為:

(10)

結合GNSS RTK手簿要想實現(xiàn)坐標的實時轉換,這里有一個假設,就是在保證測量精度的情況下,將一定測區(qū)范圍內的點位速度視為一致,然后將速度場改正信息以平移參數(shù)的形式融合到七參數(shù)轉換模型,進而實現(xiàn)GNSS測量成果向CGCS2000坐標系統(tǒng)的實時轉換。

3 實例分析

實現(xiàn)坐標轉換選擇最小二乘配置法建立速度場模型,以經緯度為單位,分別以1°×1°、5°×5°、10°×10°和整個NGS的CORS網為試驗對象,在每個區(qū)域內選擇20～30個CORS站點,獲得CGCS2000系統(tǒng)下的空間直角坐標,并對轉換結果進行精度評定。基線解算時選擇2015年第126天的數(shù)據(jù),并從SOPAC獲得點位的精確坐標。轉換結果在X、Y、Z三個方向上的點位誤差如圖2所示。

通過精度統(tǒng)計,表2給出了轉換后各個待求點的外符合精度。

表2 基于速度場轉換的待求點外符合精度 單位:m

從圖2可以看出,研究區(qū)域的大小對融合速度場信息的七參數(shù)轉換獲得的CGCS2000坐標結果的精度影響不大,各區(qū)域在X、Y、Z幾個方向上的誤差都在4 cm以內,且大部分分布在2 cm內。由表2可知，各個研究區(qū)域所有點的外符合精度都優(yōu)于2 cm。

由于在七參數(shù)改正模型中加入了速度改正信息,使得地面點位的相對位置不再因為地殼運動隨時間的變化而變化,因此在坐標轉換時不再受到轉換區(qū)域大小的影響,也有利于各區(qū)域坐標轉換參數(shù)的統(tǒng)一。

4 結束語

(1)通過對比分析幾種常用的速度場模型建立方法,考慮點位長期、整體運動趨勢及小范圍局部運動情況,采用最小二乘配置法建立速度場模型,并為CGCS2000坐標轉換提供點位速度信息。

(2)通過對不同范圍的研究區(qū)域坐標轉換精度分析,融合速度場改正信息的七參數(shù)轉換精度與區(qū)域大小無關,且可以獲取高精度的轉換結果。