嚴(yán)劍鋒 王英剛

(南京市不動(dòng)產(chǎn)登記中心, 江蘇 南京 210005)

0 引言

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

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

1 速度場(chǎng)模型建立

常用的地殼運(yùn)動(dòng)速度場(chǎng)模型的建立方法包括歐拉矢量法和擬合法。歐拉矢量法的特點(diǎn)是把整個(gè)建模區(qū)域當(dāng)成剛性塊體,認(rèn)為板塊運(yùn)動(dòng)是整體性的,內(nèi)部形變一致;擬合方法是基于空間大地測(cè)量數(shù)據(jù),通過(guò)數(shù)學(xué)擬合建立速度場(chǎng)模型,對(duì)研究區(qū)域沒(méi)有剛彈性的要求。

1.1 局部歐拉矢量法

局部歐拉矢量法就是應(yīng)用歐拉矢量基本原理,把待求點(diǎn)周圍局部范圍看成是剛體,并建立模型來(lái)獲取待求點(diǎn)的速度。

1.1.1歐拉矢量基本原理

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

(1)

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

1.1.2局部歐拉矢量法

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

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

1.2 數(shù)學(xué)擬合法

1.2.1反距離加權(quán)法

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

(2)

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

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

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

(3)

式中,V(x,y)為點(diǎn)(x,y)處的點(diǎn)位速度值;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é)點(diǎn)個(gè)數(shù);αi(i=1,2,…,u)為待估參數(shù)。根據(jù)公式(3),利用間接平差理論可以方便地求得待定參數(shù)αi,并且進(jìn)行精度評(píng)定[11]。

1.2.3最小二乘配置法

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

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

(4)

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

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

S=CuoC-1(V-AΩ)

(5)

其中,

C=Coo+Cnn

(6)

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

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

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

(7)

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

1.3 速度場(chǎng)模型精度對(duì)比

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

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

不同速度場(chǎng)模型計(jì)算得到外符合精度如表1所示。

表1 幾種速度場(chǎng)模型外符合精度 單位:mm

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

2 CGCS2000坐標(biāo)轉(zhuǎn)換

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

實(shí)時(shí)坐標(biāo)轉(zhuǎn)換參數(shù)的求取是基于速度場(chǎng)模型、坐標(biāo)框架轉(zhuǎn)換理論實(shí)現(xiàn)的。一般情況下,CGCS2000坐標(biāo)實(shí)時(shí)轉(zhuǎn)換是結(jié)合CNSS測(cè)量手簿中的七參數(shù)轉(zhuǎn)換模塊來(lái)實(shí)現(xiàn)的,因此實(shí)現(xiàn)實(shí)時(shí)轉(zhuǎn)換的關(guān)鍵部分是求取融合速度場(chǎng)的實(shí)時(shí)轉(zhuǎn)換參數(shù)。

圖1 速度場(chǎng)模型點(diǎn)位誤差圖

2.1.1 融合速度場(chǎng)信息的實(shí)時(shí)轉(zhuǎn)換參數(shù)求取

(8)

同理,若某個(gè)點(diǎn)位的變化速度為V,某坐標(biāo)系統(tǒng)定義的歷元為t0,那么實(shí)時(shí)歷元下獲取的坐標(biāo)相對(duì)于t1歷元時(shí)的坐標(biāo)的位置變化值S可表示為:

S=(t1-t0)×V

(9)

2.1.2CGCG2000坐標(biāo)實(shí)時(shí)轉(zhuǎn)換

WGS84向CGCG2000坐標(biāo)轉(zhuǎn)換即ITRF2008框架下當(dāng)前歷元的坐標(biāo)向目標(biāo)框架ITRF97下的2 000.0歷元坐標(biāo)進(jìn)行轉(zhuǎn)換。假設(shè)測(cè)站點(diǎn)的空間點(diǎn)位速度為ΔX、ΔY、ΔZ,單位為mm/年,轉(zhuǎn)換公式可表示為:

(10)

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

3 實(shí)例分析

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

通過(guò)精度統(tǒng)計(jì),表2給出了轉(zhuǎn)換后各個(gè)待求點(diǎn)的外符合精度。

表2 基于速度場(chǎng)轉(zhuǎn)換的待求點(diǎn)外符合精度 單位:m

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

由于在七參數(shù)改正模型中加入了速度改正信息,使得地面點(diǎn)位的相對(duì)位置不再因?yàn)榈貧み\(yùn)動(dòng)隨時(shí)間的變化而變化,因此在坐標(biāo)轉(zhuǎn)換時(shí)不再受到轉(zhuǎn)換區(qū)域大小的影響,也有利于各區(qū)域坐標(biāo)轉(zhuǎn)換參數(shù)的統(tǒng)一。

4 結(jié)束語(yǔ)

(1)通過(guò)對(duì)比分析幾種常用的速度場(chǎng)模型建立方法,考慮點(diǎn)位長(zhǎng)期、整體運(yùn)動(dòng)趨勢(shì)及小范圍局部運(yùn)動(dòng)情況,采用最小二乘配置法建立速度場(chǎng)模型,并為CGCS2000坐標(biāo)轉(zhuǎn)換提供點(diǎn)位速度信息。

(2)通過(guò)對(duì)不同范圍的研究區(qū)域坐標(biāo)轉(zhuǎn)換精度分析,融合速度場(chǎng)改正信息的七參數(shù)轉(zhuǎn)換精度與區(qū)域大小無(wú)關(guān),且可以獲取高精度的轉(zhuǎn)換結(jié)果。