基于Delauany三角網(wǎng)與反距離內(nèi)插實現(xiàn)GPS區(qū)域網(wǎng)中噪聲剔除

紀(jì)海源 張偉佳 梁 磊 何遠(yuǎn)梅 郝 明

(1. 陜西工業(yè)職業(yè)技術(shù)學(xué)院 土木工程學(xué)院, 陜西 咸陽 712000; 2. 陜西新旅程測繪科技有限公司,陜西 西安 712000; 3. 中國地震局第二監(jiān)測中心, 陜西 西安 710054)

0 引言

中國陸態(tài)網(wǎng)絡(luò)工程的區(qū)域站已經(jīng)有兩千多個了,這些區(qū)域站的觀測方式不是連續(xù)的,而是采取流動觀測的方式,其優(yōu)點是可以節(jié)約人力、財力、物力,觀測效率高,具有很強(qiáng)的靈活性;缺點是得到的坐標(biāo)時間序列不連續(xù),而在全球定位系統(tǒng)(Global Positioning System，GPS)觀測結(jié)果中存在非構(gòu)造噪聲是有目共睹的,時間序列的不連續(xù)無法使用QOCA(Quasi-Observation Combination Analysis)地殼形變分析軟件進(jìn)行非構(gòu)造噪聲的剔除。因此區(qū)域站速度所反映的地殼形變中包含有大量的誤差,剔除區(qū)域站時間序列中的非構(gòu)造噪聲成為國內(nèi)外學(xué)者研究的熱點[1]。田云峰等認(rèn)為GPS坐標(biāo)中非構(gòu)造噪聲存在區(qū)域相關(guān)性[2-3],于是設(shè)想到區(qū)域站分布和連續(xù)站分布的關(guān)系,是否能尋找某種方法將連續(xù)站和區(qū)域站關(guān)聯(lián)起來。 Delaunay三角形能把三角形頂點和內(nèi)部的點聯(lián)系起來,用連續(xù)站構(gòu)建Delaunay三角形成本文新的探索方法;其次采用反距離權(quán)內(nèi)插將區(qū)域站和Delaunay三角形的頂點聯(lián)系起來,用連續(xù)站的周年項和半周年項參數(shù)(QOCA軟件結(jié)果)去擬合區(qū)域站的周年項和半周年項參數(shù),最終剔除區(qū)域站中的部分非構(gòu)造噪聲,提高GPS流動觀測區(qū)域站三維運動速度場的精度,為深入認(rèn)識地塊動力學(xué)機(jī)理提供重要基礎(chǔ)資料。

1 連續(xù)站時間序列參數(shù)的提取

本文研究區(qū)域為中國地震監(jiān)測陸態(tài)網(wǎng)在青藏高原東北緣地區(qū)(101°～107°E,32°～40°N)分布的GPS連續(xù)站和流動觀測區(qū)域站,分布如圖1所示;采用GAMIT/GLOBK軟件處理13個連續(xù)站和153個流動觀測區(qū)域站,得到連續(xù)站高精度位移序列,再使用QOCA地殼分析軟件[2]的analyze_tseri和pca模塊剔除GPS連續(xù)站中的非構(gòu)造噪聲,在QOCA軟件的結(jié)果文件中,將三個方向分析擬合出來的周年項、半周年項參數(shù)等提取出來(表1～3)。s1、s2和s3、s4分別為年、半年周期項系數(shù)。

圖1 汾渭地區(qū)GPS連續(xù)站和區(qū)域站分布圖注：五角形為GPS連續(xù)站；圓形為區(qū)域站。

表1 北向參數(shù)提取結(jié)果

表2 東向參數(shù)提取結(jié)果

表3 垂向參數(shù)提取結(jié)果

2 Delaunay三角網(wǎng)的構(gòu)建

2.1 VB實現(xiàn)Delaunay三角網(wǎng)

狄洛尼(Delaunay)三角網(wǎng)是Voronoi圖的對偶圖,由對應(yīng)的Voronoi多邊形共邊的點連接而成。狄洛尼三角形由三個相鄰點連接而成,這三個相鄰點對應(yīng)的Voronoi多邊形有一個公共頂點,此頂點同時也是狄洛尼三角形外接圓的圓心[4-5]。本文采用VB(Visual Basic)編程實現(xiàn)了13個連續(xù)站狄洛尼(Delaunay)三角網(wǎng)的生成(圖2)。

圖2 Delaunay三角網(wǎng)

2.2 區(qū)域站在Delaunay三角網(wǎng)中的識別

在生成三角網(wǎng)之后,需要判斷每個區(qū)域站隸屬于哪個三角形。判斷點在三角形內(nèi)的算法有好多種,在這里選擇面積法:如果區(qū)域站與三角形三個頂點的連線構(gòu)成三個小三角形的面積之和與大三角形的面積相等,則確定點落在這個三角形里面(圖3)。

圖3 部分區(qū)域站在Delaunay三角網(wǎng)中的分布

3 根據(jù)連續(xù)站擬合區(qū)域站參數(shù)

3.1 反距離權(quán)算法簡介

“反距離加權(quán)算法”又稱為倒數(shù)距離加權(quán)插值或“Shepard方法”[6]。設(shè)有n個點,平面坐標(biāo)為(xi,yi),垂直高度為zi,i=1,2,…,n,倒數(shù)距離加權(quán)插值的插值函數(shù)為

f(x,y)=

3.2 VB實現(xiàn)區(qū)域站參數(shù)擬合

將提取的連續(xù)站周年項、半周年項參數(shù)準(zhǔn)備好,VB編寫程序?qū)崿F(xiàn)反距離權(quán)內(nèi)插的程序,用三個頂點的參數(shù)擬合出區(qū)域站的周年項、半周年項參數(shù)。表4～6是提取的部分三角形三個方向擬合的區(qū)域站參數(shù)結(jié)果。

表4 北向參數(shù)擬合結(jié)果

表5 東向參數(shù)擬合結(jié)果

表6 垂向參數(shù)提取結(jié)果

4 區(qū)域站非構(gòu)造噪聲的剔除

在上一節(jié)已經(jīng)將區(qū)域站的周年項、半周年項參數(shù)擬合出來了,可以應(yīng)用數(shù)學(xué)模型計算其中包含的非構(gòu)造噪聲,然后從位移序列中剔除;對位移序列分量每日解觀測序列建立參數(shù)模型:

y(ti)=a+bti+s1×sin(2πti)+s2×

cos(2πti)+s3×sin(4πti)+s4×cos(4πti)(2)

式中,ti(i=1…N)為以年為單位的時間,在這里我們以每日的單日解構(gòu)成時間序列;a為初始位置;b為速率;s1、s2和s3、s4分別為年、半年周期項系數(shù)[7]。

將表5中擬合的參數(shù)按照測站代入位移序列的數(shù)學(xué)模型,把計算出來的非構(gòu)造噪聲從位移序列y(ti)中剔除,這樣可以得到相對精確的地殼運動形變,下面舉例G120站的剔除結(jié)果。

從表7的噪聲列可以看到北南方向誤差最大可以剔除掉1.2 mm,最小為0.5 mm,噪聲的均方根為0.89 mm;東西方向誤差最大可以剔除掉1.0 mm,最小為0.2 mm,噪聲的均方根為0.66 mm;垂直向最大可以剔除掉3.5 mm,最小1.2 mm,噪聲的均方根為2.44 mm。噪聲在水平上的均方根為1.11 mm。垂直方向的噪聲均方根是水平方向的兩倍多,說明這種方法在垂直方向上影響比水平向影響更為明顯;這也能說明GPS在垂直方向監(jiān)測中包含的誤差要比水平方向包含的誤差大。

表7 G120站剔除非構(gòu)造噪聲結(jié)果 單位：mm

5 剔除前后結(jié)果對比分析

將剔除非構(gòu)造噪聲前后的153個區(qū)域站位移序列分別準(zhǔn)備在兩個文件里面,用Matlab編寫程序?qū)崿F(xiàn)批處理擬合區(qū)域站的線性速度[8-9],結(jié)果見表8。

表8 部分區(qū)域站剔除非構(gòu)造噪聲前后速度對比表 單位：mm/年

對剔除后的時間序列用加權(quán)最小二乘擬合速度[10],經(jīng)過對比153個站發(fā)現(xiàn)(表10只列出部分):在南北方向上速度差最大為0.8 mm最小為0 mm,速度差均方根0.33 mm;在東西方向上速度差最大為1.0 mm,最小為0.1 mm,速度差均方根為0.55 mm;在垂直方向上速度差最大為-3.0 mm,最小為 0 mm,速度差均方根為1.71 mm;水平方向上的速度差均方根為0.58 mm;垂直方向上的速度差均方根是水平方向上的2.13倍;這個結(jié)論也說明在GPS在高程方向上包含的非構(gòu)造噪聲要大于水平方向上的。

通過Delaunay三角網(wǎng)和反距離內(nèi)插算法相結(jié)合的處理,發(fā)現(xiàn)此方法能夠從區(qū)域站中剔除部分非構(gòu)造誤差;同時,與QOCA剔除連續(xù)站的結(jié)果對比發(fā)現(xiàn):前者在水平速度差均方根上為后者的2.0倍,在垂直速度差均方根上為后者的2.7倍,說明Delaunay三角網(wǎng)和反距離內(nèi)插方法剔除區(qū)域站非構(gòu)造噪聲明顯比QOCA軟件剔除連續(xù)站的非構(gòu)造噪聲大兩倍多。至于這兩種方法存在差距的原因和這種新方法在區(qū)域站非構(gòu)造噪聲剔除中應(yīng)用的可靠性還有待于深入研究,筆者將繼續(xù)改變網(wǎng)形以及內(nèi)插的方法,對比分析區(qū)域站和連續(xù)站位移時間序列中存在的關(guān)系。

6 結(jié)束語

使用QOCA軟件能夠提取GPS連續(xù)站的時間序列,借助連續(xù)站時間序列參數(shù)構(gòu)建Delaunay三角網(wǎng),結(jié)合反距離內(nèi)插算法能夠?qū)崿F(xiàn)GPS流動觀測區(qū)域站中非構(gòu)造噪聲的剔除;GPS區(qū)域站中非構(gòu)造噪聲的剔除,提高了地殼三維運動速度的精度,為地殼運動研究以及地震預(yù)報提供了基礎(chǔ)數(shù)據(jù)。