水平運動場濾波試驗與應變計算＊

楊 博 張風霜 占 偉 韓月萍

(中國地震局第一監(jiān)測中心,天津 300180)

水平運動場濾波試驗與應變計算＊

楊 博 張風霜 占 偉 韓月萍

(中國地震局第一監(jiān)測中心,天津 300180)

GNSS觀測資料為揭示地殼形變提供了必要依據(jù),但由于噪聲的存在使得地殼形變信息難以恰當?shù)乇焕没颢@取。從連續(xù)變化的角度出發(fā),提出利用多核函數(shù)進行濾波與解析。該方法有如下特點：1)可實現(xiàn)對水平運動場進行不同層次的濾波與信息分離;2)可獲得水平運動場的數(shù)學解析;3)可便利地計算連續(xù)應變場的各類應變參量。利用華北地區(qū) GNSS結(jié)果,給出利用該方法進行濾波與應變計算的實際算例。

GNSS;水平運動;多核函數(shù);濾波;應變計算

1 引言

在地球動力學研究領域,以 GNSS觀測資料為依托的研究已是一個重要分支;就地殼形變研究而言,利用 GNSS資料人們有望從空間與時間的角度觀察地殼的運動及形變的詳細過程。因此,多角度的研究也相應展開[1-12]。在早期研究中,有的學者提出將形變場分塊研究的方法[13],實質(zhì)上是假定各塊內(nèi)部介質(zhì)各向同性,這使我們對塊體的形變特征有一個總體的印象。然而,這種結(jié)果畢竟不是空間上的細部描述,真實的地殼介質(zhì)也非各向同性,故而顯得比較粗略。隨著研究的深入,近年來人們希望能夠更加客觀地反映形變場連續(xù)的細部變形,故而提出了一些新的方法,如李延興[14]給出了可描述地殼趨勢性連續(xù)應變的模型,楊國華[7]在此基礎上進一步提出了非線性應變連續(xù)模型,但從實質(zhì)上講都未離開低頻應變的范疇;江在森[6]提出的用最小二乘配置法解求應變的方法,使應變場的描述更加可觀,但從頻譜域的角度來看,所提取的信息仍屬于中低頻段內(nèi)的信息;Savage[16]也給出了在球面上計算應變得方法。然而,它們均存在著 10-9/a量級的偏差[17];此外,從頻譜域的角度來看信息也不夠完整。石耀霖、劉序儼[18,19]提出了在球面或橢球面計算水平應變場的算式,在理論上解決了上述問題。從技術(shù)層面來看,這些方法均是以運動“觀測值”等為依托,直接進行應變場的計算。但是,由于不同作者資料選用的不同或模型選用的不同等,使計算結(jié)果存在一定的差別,并在一定程度上也導致了應變分析的混亂;此外,在分析時又往往忽略對應變場相應的運動場的計算與分析,故而不免有些缺憾。因此,本文試圖從運動“觀測值”入手,通過數(shù)理方法,結(jié)合物理含義獲取較高信噪比的運動場,實現(xiàn)無偏應變場的解析算法。

2 水平運動場的解析描述與濾波計算

2.1 多核函數(shù)與水平運動場的解析表達

多核函數(shù)法的基本思想是,任何數(shù)學表面和任何不規(guī)則的圓滑表面,總可以通過一系列有規(guī)則的數(shù)學表面的合成,以任意精度逼近。其數(shù)學表達式為：

式中,sj(x,y,xj,yj)為核函數(shù),(xj,yj)為核點的坐標,cj為待定系數(shù)。

一般說來核函數(shù)是一個非常簡單的曲面形態(tài),通常為“缽”形或“鐘”形曲面函數(shù)。盡管哪一種形態(tài)的核函數(shù)更優(yōu)尚無非常明確的定論,但在此后的應用研究中,該法已被廣泛地用于地殼形變研究[20,21]、DEM內(nèi)插[22]和地球重力場元內(nèi)插[23],并被證明是非常成功和有效的。

然而,多核函數(shù)在具體應用中的效果如何仍存在選擇問題。對于不同的描述對象可能會有不同的選擇,即使對相同的對象也并非是唯一選擇,甚至很難區(qū)分它們之間的優(yōu)劣。所以,在實際應用中主要取決于使用者對多核函數(shù)及描述對象的理解。一般說來它取決于 3個方面:1)核函數(shù)的選擇;2)核函數(shù)中的參數(shù)選擇;3)核函數(shù)的核點選擇。由于我們所描述的對象是較大尺度的地殼應變場,而地殼介質(zhì)由于其各向非同性、活動構(gòu)造帶的存在及應力大小與方向的空間變化等,使得應變場具有廣譜性。這就是說,我們所描述的對象是很復雜的。據(jù)此,這里所選用的核函數(shù)為“缽”形函數(shù),具體為:

式中:dj為球面上兩點間的大地線長度 (以 km為單位);ST=(s1,s2,…,snx);C=(c1,c2,…,cnx)。

在 GNSS計算中,通常是以 ITRF作為參考框架獲取速度場,而在實際分析時為了突出相對變化,往往從中剔出歐拉運動,即相對于區(qū)域整體無旋轉(zhuǎn)基準[24]。實現(xiàn)假定 (λj,φj)為測站 j的位置坐標,ve(λj,φj)、vn(λj,φj)為其相應的東向和北向相對運動,Voe(λj,φj)、Von(λj,φj)為相應的歐拉運動。由于歐拉運動是通過解析式計算,故只需對相對運動進行解析即可。在用多核函數(shù)進行數(shù)值逼近時應將所有測站位置作為核函數(shù)的核點位置,因此對于任一方向均可列出 n個方程 (測站的個數(shù)),故可求解n個未知數(shù)(C)。這時可獲得東、北向相對運動的解析式 fe(λ,φ)、fn(λ,φ)。由此即可計算研究區(qū)任意位置的相對水平運動結(jié)果。即

同理,也可獲得相應誤差的解析式:

2.2 水平運動的濾波計算

由于上述數(shù)值逼近所得到的是水平運動觀測值的解析結(jié)果,它不僅包含運動信息,同時也包含誤差干擾,所以由此得到的運動結(jié)果實際上是各種成分的綜合結(jié)果。因此,進行濾波與信息分離是必要的。多核函數(shù)不但可以進行數(shù)值逼近,而且也具有濾波功能,故在進行空間濾波時同樣利用多核函數(shù)法。濾波計算時,濾波效果或濾波層次取決于核點的選用與分布。從性質(zhì)上看,核點的間隔越小 (不得小于測站間的平均間隔),由此所獲得的運動信息所包含的頻譜帶寬越寬,可涵蓋從低頻至高頻的范圍;反之帶寬越窄,則越趨于低頻。在進行濾波計算時,由于運動場是由東、北兩個分量共同描述,所以為了保持運動場描述的協(xié)調(diào)性和彼此之間的相關性不被破壞,進行兩分量的濾波計算時應采用相同的濾波計算準則。故東、北向運動分量的濾波函數(shù)均表述為:

式中：AT=(a1,a2,…,anx)、BT=(b1,b2,…,bnx),均為待定系數(shù);ST=(s1,s2,…,snx)為核函數(shù)陣;si(λ, φ,λi,φi)=exp(-(di/780)2.8)(di為球面上兩點間的大地線長度,以 km為單位),(λi,φi)為核點位置坐標。

在實際計算時,首先以相同的網(wǎng)格化形式 (等步長)獲取兩個分量的數(shù)據(jù),并以此作為“觀測值”。然后,根據(jù)最小二乘法求解上述待定系數(shù),即對每個分量通過誤差方程

建立相應的法方程

并由此求解 X(待定系數(shù))：

式中,N=DTPL、W=DTPD、L為相應分量觀測值的矩陣,P為相應分量觀測值的權(quán)矩陣 (pi=1/m2i,mi由式(5)獲得)。當 X求得后,式 (6)變?yōu)榻?jīng)濾波后運動計算的函數(shù)解析式,具體為:

據(jù)此可獲得任意位置水平運動濾波值。

根據(jù)誤差傳播律,可進行濾波值的誤差估計。由式(7)可計算相應運動分量的單位權(quán)方差估值,東向為:

北向為 :

式中,nx為核點的個數(shù),n為測站數(shù)量,且大于 nx。式(11)、(12)中的殘差均是測站點位置的相應殘差。令

式中 j=ee,nn。東向運動方差為:

北向運動方差為：

由式(14)、(15)及測站東、北分量間的相關系數(shù)即可計算濾波運動的誤差橢圓。

2.3 應變場的計算

在球面上進行應變計算應參照石耀霖、劉序儼[18,19]所給出的應變算式。但由于其中存在非微分項的運動,所以其大小必然與運動的參考基準有關。關于該問題楊國華[17]指出,只要依據(jù) ITRF參考框架所得到的水平運動 (位移)即可獲得無偏應變計算結(jié)果。由于文獻[18,19]所給的算式是依據(jù)余緯度而建立的算式,此外其坐標系的使用也與現(xiàn)行的不一致,故使用起來不太方便。為此,進行變換是必要的。假定東西向應變?yōu)棣舉(λ,φ),南北向應變?yōu)棣舗(λ,φ),它們之間的剪應變?yōu)棣舉n(λ,φ)。濾波后在 ITRF框架下運動的解析式為:

其中

在現(xiàn)行球面坐標系統(tǒng)下球面應變算式為:

式中R為地球的平均半徑。

在橢球面上為:

式中,Rn為卯酉圈曲率半徑,Rm為子午圈曲率半徑,λ為經(jīng)度,φ為緯度。

最大主應變?yōu)?

最小主應變?yōu)?

最大剪應變?yōu)?

面應變?yōu)?

最大主應變方向為:

3 算例

以 2003—2007年華北地區(qū) GPS流動測站在ITRF參考框架的平均水平運動速度結(jié)果為例進行分析。為便于觀察,首先剔除其中的歐拉運動,獲得相對于區(qū)域無旋轉(zhuǎn)基準的運動結(jié)果 (圖 1),然后再依據(jù)式(3)、(4)對研究區(qū)相對運動速度及誤差結(jié)果進行數(shù)據(jù)逼近獲得解析式。在此基礎上結(jié)合式 (6)及以經(jīng)、緯步長分別為 0.2°計算網(wǎng)格點上的運動結(jié)果,再進行空間濾波計算。作為算例,這里選用的核點經(jīng)、緯步長分別為 1.0°和 2.0°,依據(jù)式 (6)～(10)進行計算,結(jié)果見圖 2～3。

圖 1 未經(jīng)濾波處理的相對區(qū)域無旋轉(zhuǎn)基準的水平運動結(jié)果(誤差橢圓置信度為 95%)Fig.1 Horizontal movements without filtering and relative to no-rotation datum (confidence of error ellipse is 95%)

圖 2 步長為 1.0°核點間隔的空間濾波(a)(誤差橢圓置信度為 95%)與殘差結(jié)果(b)Fig.2 Space filtering results of horizontal motion by kernel point with step-size 1.0°(a)(confidence of error ellipse is 95%)and residuals(b)

圖 1是未經(jīng)處理的相對區(qū)域無旋轉(zhuǎn)基準的相對水平運動。眾所周知,華北地區(qū)是一個形變信息較弱、干擾較強的地區(qū)。從單站運動誤差的角度來看,大部分測站的運動信息都在限差范圍內(nèi)。如果由此判斷運動場沒有比較顯著的信息可能會犯錯誤。其一,根據(jù) GNSS計算軟件的處理過程及后處理來看,測站間誤差的關聯(lián)性并非完全恰當(后處理時常常難以顧及前期處理的協(xié)方差矩陣);其二,如果完全由誤差所致,則運動場結(jié)果不太可能存在某種有序的變化。所以對觀測資料作進一步處理,以削弱噪聲干擾、凸顯有序的運動信息是必要的。圖 2(a)是核點間隔步長為 1.0°濾波處理結(jié)果及相應的誤差橢圓。從誤差的角度來看,圖 2(a)所展現(xiàn)的運動誤差明顯減小,從運動的有序性來看,運動的趨勢性變得較為明朗,因此可使我們比較容易地把握其運動的基本特征。圖 2(b)是與此濾波相應的殘差結(jié)果,不難發(fā)現(xiàn)它具有典型的噪聲特點。僅此已初步表明,本文的方法是可行的,得到的結(jié)果也是有效的。為了獲得更多的認識,給出核點間隔步長為 2.0°濾波處理結(jié)果、相應的誤差橢圓及殘差 (圖 3)。由圖3(a)可知其信噪比更高,有序變化的特征更加突出。從殘差的分布也可看出,噪聲的特性依然突出(其中或多或少地包含著一些高頻成分)。然而,在實際分析中什么樣的結(jié)果為最佳?難以給出確定的結(jié)論,因為它與研究的具體對象或目標有關,故不能一概而論。從濾波的角度來看,核點間隔越長所體現(xiàn)的運動信息越趨于低頻。

圖 3 步長為 2.0°核點間隔的空間濾波(a)(誤差橢圓置信度為 95%)與殘差結(jié)果(b)Fig.3 Space filtering results of horizontal motion by kernel point with step-size 2.0°(a)(confidence of error ellipse is 95%)and residuals(b)

圖 4 步長分別為 1.0°(a)和 2.0°(b)核點間隔的東向應變率結(jié)果Fig.4 Space filtering results of the east strain rate by kernel pointwith step-size 1.0°(a)and 2.0°(b)

由于濾波結(jié)果是通過模型解析式的計算獲得,所以根據(jù)此解析式即可非常便利地計算與其相應的應變場。圖 4所展現(xiàn)的是在上述兩種核點間隔情況下的東向應變率結(jié)果(圖中黑點為 GNSS測站,藍點表示一些城市的位置,灰粗線為該區(qū)主要活動斷層,黑粗線為海岸線)。相比之下,圖 4(a)所包含的頻譜范圍較寬,較高頻成分相對突出,曲面雖光滑但起伏相對較大;其數(shù)值較大,一般為 10-8/a,梯度也較大。圖 4(b)的趨勢性變化特征則顯現(xiàn)得更加清晰,但變化較為寬緩,數(shù)值也較小,一般為 10-9/a。由此可進一步看出,核點間的間隔越小,包含不同頻譜形變信息就越全面;核點間隔越大,較高頻成分被濾掉的就越多。一般地說,進行濾波通常主要考慮兩個方面的問題,即噪聲削弱與不同層次的形變信息分離。從這層意義講,圖 4(a)主要體現(xiàn)的是對噪聲削弱后的處理結(jié)果,因此不同頻譜形變的成分也得以較全面地保留。但這樣的結(jié)果有時難以把握趨勢的特征性,或與研究目標所需的某種信息顯得不夠突出。圖 4(b)主要顯示較低頻的形變信息,趨勢性變化較為清晰,較高頻的信息已被削弱。這說明可通過設定核點的不同間隔來滿足不同的需求。若以研究趨勢性變化為目標,則圖 4(b)的結(jié)果為最佳。圖5是核點步長為2.0°的主應變計算結(jié)果,由此看出主應變方向在空間上不是恒常不變的,這不但表明了介質(zhì)并非各向同性,同時也展現(xiàn)了主應變場在空間上的變化趨向,可使人們更加細致地了解其變化特征。

圖 5 步長 2.0°核點間隔的主應變率結(jié)果Fig.5 Space filtering results of the main strain rate by kernel pointwith step-size 2.0°

4 結(jié)語

綜上所述,本文通過多核函數(shù)法實現(xiàn)了運動場的解析和濾波計算,同時也給出了運動場空間濾波結(jié)果嚴密精度評定算式。算例表明：本文方法具有降噪與信息分離的功能。此外,由于濾波所依據(jù)的是可解析模型,所以又可非常方便地獲得相應的應變場結(jié)果。但是,本方法更加完美的使用還需做進一步的工作,如不同情況下核函數(shù)的最優(yōu)選用,網(wǎng)格的劃分、核函數(shù)點的選擇與濾波和信息分離的關系等。

致謝 感謝楊國華研究員的指導!

1 敬少群,等.GPS時間序列及其對昆侖山口西 8.1級地震的響應[J].地震學報,2005,27(4)：394-401.(Jing Shaoqun,et al.GPS Time-series and its response toMs8.1 Kunlunshan earthquake[J]. Acta Seis mologica Sinica, 2005,27(4)：394-401)

2 顧國華,張晶.中國地殼運動觀測網(wǎng)絡基準站 GPS觀測的時間序列結(jié)果 [J].大地測量與地球動力學,2002, (2)：61-67.(Gu Guohua and Zhang Jing.Time series of displacements from GPS observation at fiducial station in the crustalmovement observation network of China[J].Journal of Geodesy and Geodynamics,2002,(2)：61-67)

3 楊國華,等.利用 GPS連續(xù)觀測資料進行強震危險性預測的探討[J].地震,2008,28(1)：33-39.(Yang Guohua, et al.Discussion on the hazard forecastof strong earthquakes by GPS continuous observation data[J].Earthquake,2008, 28(1)：33-39)

4 黃立人,王敏.中國大陸構(gòu)造塊體的現(xiàn)今活動和變形[J].地震地質(zhì),2003,25(1)：23-32.(Huang Liren and WangMin.Pesent-day activity and defor mation of tectonic blocks in China continent[J].Seismology and Geology, 2003,25(1)：23-32)

5 楊國華,等.由 GPS觀測結(jié)果推導中國大陸現(xiàn)今水平應變場[J].地震學報,2002,24(4)：337-347.(Yang Guohua,et al.Current horizontal strain field in Chinese mainland derived from GPS data[J].Acta Seismologica Sinica, 2002,24(4)：337-347)

6 江在森,馬宗晉,張希.GPS初步結(jié)果揭示的中國大陸水平應變場與構(gòu)造變形[J].地球物理學報,2003,46(3)：352-358.(Jiang Zaisen,Ma Zongjin and Zhang Xi.Horizontal strain field and tectonic deformation of Chinamainland revealed by preliminary GPS result[J].Chinese Journal of Geopysics,2003,46(3)：352-358)

7 楊國華,江在森,王敏.印尼地震對我國川滇地區(qū)地殼水平活動的影響[J].大地測量與地球動力學,2006,(1)：9 -14.(Yang Guohua,Jiang Zaisen andWangMin.Effectof Indonesia earthquake on horizontal crustal movement in Sichuan-Yunnan region[J].Journal of Geodesy and Geodynamics,2006,(1)：9-14)

8 張希,等.用最小二乘配置獲得地形變應變場動態(tài)圖像的幾個問題研究[J].地殼形變與地震,1999,(3)：32-39.(Zhang Xi,et al.Study on some questions of dynamic pictures of crustal deformation and strain fields obtained by the least square collocation[J].Crustal Defor mation And Earthquake,1999,(3)：32-39)

9 楊國華,等.汶川地震對華北地區(qū)水平形變場影響及有關含義的討論 [J].地震,2009,29(1)：77-83.(Yang Guohua,et al. Influence of Wenchuan earthquake on the horizontal deformation field in North China and its implications[J].Earthquake,2009,29(1)：77-83)

10 楊國華,等.華北地區(qū)的水平運動場與昆侖山 8.1級地震的可能關系[J].大地測量與地球動力學,2007,(2)：10-15. (Yang Guohua,et al.Possible relation of horizontalmovement field in North China to Kunlun mountain Ms8.1 earthquake[J].Journal of Geodesy and Geodynamics,2007,(2)：10-15)

11 王敏,等.現(xiàn)今中國大陸地殼運動與活動塊體模型[J].中國科學 (D輯),2003,33(增)：21-32.(WangMin,et al.Contemporary crustal deformation of the Chinese continent and tectonic blockmodel[J].Science in China(Series D),2003,46(Supp)：21-32)

12 王琪,等.天山現(xiàn)今地殼快速縮短與南北地塊的相對運動[J].科學通報,2000,45(14)：1 543-1 547.(Wang Qi,et al.Present-day fast crustal constriction of Tianshan and the motion of north-south block[J].Chinese Seince Bullutin,2000,45(14)：1 543-1 547)

13 李延興,等.中國大陸活動塊體的運動與應變狀態(tài)[J].中國科學 (D輯),2003,33(增)：65-81.(Li Yangxin,et al.Movement and strain conditions of active blocks in the Chinese mainland[J].Science in China(SeriesD),2003, 46(Supp)：65-81)

14 李延興,等.中國大陸及周邊地區(qū)的水平應變場[J].地球物理學報,2004,47(2)：222-231.(Li Yangxing,et al.Horizontal strain field in the Chinese mainland and its surrounding areas[J]. Chinese Journal of Geopysics, 2004,47(2)：222-231)

15 楊博,等.斷裂區(qū)帶變形分析方法及應用[J].西北地震學報,(待刊).(Yang Bo,et al.Method and application on defor mation analysis for fault zone[J].Northwestern Seismological Journal,(to be published))

16 Savage J C,GanW and Svare J L.Strain accumulation and rotation in the Eastern California Shear Zone[J].J.Geophys.Res.,2001,106(B10)：21 995-22 007

17 楊國華,等.應變計算與分析的若干問題及有關偏差的修正[J].大地測量與地球動力學,2010,(4)：59-63. (Yang Guohua,et al.Several problems of strain calculation analysis and the correction of related deviation[J]. Journal of Geodesy and Geodynamics,2010,(4)：59-63)

18 石耀霖,朱守彪.用 GPS位移資料計算應變方法的討論[J].大地測量與地球動力學,2006,(1)：1-8.(Shi Yaolin and Zhu Shoubiao.Discussion onmethod of calculating strain with GPS data[J].Journal of Geodesy and Geodynamics,2006,(1)：1-8)

19 劉序儼,黃聲明,梁全強.旋轉(zhuǎn)托球面上的應變與轉(zhuǎn)動張量表達 [J].地震學報,2007,29(3)：240-249.(Liu Xuyan,Huang Shengming and Liang Quanqiang.Expression of strain and rotation tensor in geodetic coordinates [J].Acta Seis mologica Sinica,2007,29(3)：240-249.

20 黃立人,劉天奎.用速度面擬合法研究華北部分地區(qū)的現(xiàn)今地殼垂直運動[J].地殼形變與地震,1989,(3)：6 -12.(HuangLiren and Liu Tiankui.Research on present -day vertical crustalmotion of part of North China by fitting method of velocity surface[J].Crustal Deformation And Earthquake,1989,(3)：6-12)

21 楊國華,黃立人.速度面擬合法中多面函數(shù)幾個特性的初步數(shù)值研究[J].地殼形變與地震,1990,10(4)：70-82.(Yang Guohua and Huang Liren.Primary numerical study of several characteristics about multi-surface function in fitting method of velocity surface[J].Crustal Deformation And Earthquake,1990,10(4)：70-82)

22 周光文,黃莜容.一種改進的多面函數(shù)法[J].測繪通報, 1996,(6)：6-8.(Zhou Guangwen and Huang Yiurong.A kind of improvingmethod aboutmulti-surface function[J]. Bulletin of Surveying andMapping,1996,(6)：6-8)

23 陶本藻.GPS水準似大地水準面擬合和正常高計算[J].測繪通報,1992,(4)：14-18.(Tao Benzao.Calculation of normal height and GPS quasi-geoid fitting[J].Bulletin of Surveying andMapping,1992,(4)：14-18)

24 楊國華,等.中國大陸整體無凈旋轉(zhuǎn)基準及其應用[J].大地測量與地球動力學,2005,(4)：6-10.(Yang Guohua,et al.No-net-rotation on crustal movement of China mainland and its application[J].Journal of Geodesy and Geodynamics,2005,(4)：6-10)

EXPERIM ENT FOR FILTERING HORIZONTAL MOVEM ENT FIELD AND ITS STRA IN CALCULATION

YangBo,Zhang Fengshuang,ZhanWei and Han Yueping
(First CrustM onitoring and Application Center,CEA,Tianjin 300180)

GNSS data provides uswith the necessary foundation for revealing crustal deformation,but it is difficult to utilize or acquire appropriately various crustal defor mation information because of the existence of all kinds of noise.From the aspect of continuous variation,filtering and analyzing by use of the multi-kernel function are proposed.Thismethod some has characteristics as follows.1)The filtering and information separation at every different level can be realized.2)The analytic expression of the horizontalmovement field can be obtained.3)The various strain parameters of the continuous strain field can be conveniently computed.The example of filtering and strain compution with thismethod based on GNSS results in North China is given.

GNSS;horizontalmovement;multi-kernel function;filtering;strain calculation

1671-5942(2010)05-0106-07

2010-04-23

中國地震局“中國大陸 7、8級地震中長期危險性預測研究”;“強化華北地區(qū)強震監(jiān)視跟蹤”;“十一五”國家科技支撐課題(2006BAC01B03-01-01)

楊博,男,1985年生,主要從事 GPS應用與數(shù)據(jù)處理研究.E-mail：boyyang1234@qq.com

P207

A