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

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

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

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

圖2 Delaunay三角網(wǎng)

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

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

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

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

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

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

f(x,y)=

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

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

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

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

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

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

在上一節(jié)已經(jīng)將區(qū)域站的周年項(xiàng)、半周年項(xiàng)參數(shù)擬合出來了,可以應(yīng)用數(shù)學(xué)模型計(jì)算其中包含的非構(gòu)造噪聲,然后從位移序列中剔除;對(duì)位移序列分量每日解觀測序列建立參數(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)為以年為單位的時(shí)間,在這里我們以每日的單日解構(gòu)成時(shí)間序列;a為初始位置;b為速率;s1、s2和s3、s4分別為年、半年周期項(xiàng)系數(shù)[7]。

將表5中擬合的參數(shù)按照測站代入位移序列的數(shù)學(xué)模型,把計(jì)算出來的非構(gòu)造噪聲從位移序列y(ti)中剔除,這樣可以得到相對(duì)精確的地殼運(yùn)動(dòng)形變,下面舉例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é)果對(duì)比分析

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

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

對(duì)剔除后的時(shí)間序列用加權(quán)最小二乘擬合速度[10],經(jīng)過對(duì)比153個(gè)站發(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倍;這個(gè)結(jié)論也說明在GPS在高程方向上包含的非構(gòu)造噪聲要大于水平方向上的。

通過Delaunay三角網(wǎng)和反距離內(nèi)插算法相結(jié)合的處理,發(fā)現(xiàn)此方法能夠從區(qū)域站中剔除部分非構(gòu)造誤差;同時(shí),與QOCA剔除連續(xù)站的結(jié)果對(duì)比發(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)插的方法,對(duì)比分析區(qū)域站和連續(xù)站位移時(shí)間序列中存在的關(guān)系。

6 結(jié)束語

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