水平运动场滤波试验与应变计算＊

杨 博 张风霜 占 伟 韩月萍

(中国地震局第一监测中心,天津 300180)

水平运动场滤波试验与应变计算＊

杨 博 张风霜 占 伟 韩月萍

(中国地震局第一监测中心,天津 300180)

GNSS观测资料为揭示地壳形变提供了必要依据,但由于噪声的存在使得地壳形变信息难以恰当地被利用或获取。从连续变化的角度出发,提出利用多核函数进行滤波与解析。该方法有如下特点：1)可实现对水平运动场进行不同层次的滤波与信息分离;2)可获得水平运动场的数学解析;3)可便利地计算连续应变场的各类应变参量。利用华北地区 GNSS结果,给出利用该方法进行滤波与应变计算的实际算例。

GNSS;水平运动;多核函数;滤波;应变计算

1 引言

在地球动力学研究领域,以 GNSS观测资料为依托的研究已是一个重要分支;就地壳形变研究而言,利用 GNSS资料人们有望从空间与时间的角度观察地壳的运动及形变的详细过程。因此,多角度的研究也相应展开[1-12]。在早期研究中,有的学者提出将形变场分块研究的方法[13],实质上是假定各块内部介质各向同性,这使我们对块体的形变特征有一个总体的印象。然而,这种结果毕竟不是空间上的细部描述,真实的地壳介质也非各向同性,故而显得比较粗略。随着研究的深入,近年来人们希望能够更加客观地反映形变场连续的细部变形,故而提出了一些新的方法,如李延兴[14]给出了可描述地壳趋势性连续应变的模型,杨国华[7]在此基础上进一步提出了非线性应变连续模型,但从实质上讲都未离开低频应变的范畴;江在森[6]提出的用最小二乘配置法解求应变的方法,使应变场的描述更加可观,但从频谱域的角度来看,所提取的信息仍属于中低频段内的信息;Savage[16]也给出了在球面上计算应变得方法。然而,它们均存在着 10-9/a量级的偏差[17];此外,从频谱域的角度来看信息也不够完整。石耀霖、刘序俨[18,19]提出了在球面或椭球面计算水平应变场的算式,在理论上解决了上述问题。从技术层面来看,这些方法均是以运动“观测值”等为依托,直接进行应变场的计算。但是,由于不同作者资料选用的不同或模型选用的不同等,使计算结果存在一定的差别,并在一定程度上也导致了应变分析的混乱;此外,在分析时又往往忽略对应变场相应的运动场的计算与分析,故而不免有些缺憾。因此,本文试图从运动“观测值”入手,通过数理方法,结合物理含义获取较高信噪比的运动场,实现无偏应变场的解析算法。

2 水平运动场的解析描述与滤波计算

2.1 多核函数与水平运动场的解析表达

多核函数法的基本思想是,任何数学表面和任何不规则的圆滑表面,总可以通过一系列有规则的数学表面的合成,以任意精度逼近。其数学表达式为：

式中,sj(x,y,xj,yj)为核函数,(xj,yj)为核点的坐标,cj为待定系数。

一般说来核函数是一个非常简单的曲面形态,通常为“钵”形或“钟”形曲面函数。尽管哪一种形态的核函数更优尚无非常明确的定论,但在此后的应用研究中,该法已被广泛地用于地壳形变研究[20,21]、DEM内插[22]和地球重力场元内插[23],并被证明是非常成功和有效的。

然而,多核函数在具体应用中的效果如何仍存在选择问题。对于不同的描述对象可能会有不同的选择,即使对相同的对象也并非是唯一选择,甚至很难区分它们之间的优劣。所以,在实际应用中主要取决于使用者对多核函数及描述对象的理解。一般说来它取决于 3个方面:1)核函数的选择;2)核函数中的参数选择;3)核函数的核点选择。由于我们所描述的对象是较大尺度的地壳应变场,而地壳介质由于其各向非同性、活动构造带的存在及应力大小与方向的空间变化等,使得应变场具有广谱性。这就是说,我们所描述的对象是很复杂的。据此,这里所选用的核函数为“钵”形函数,具体为:

式中:dj为球面上两点间的大地线长度 (以 km为单位);ST=(s1,s2,…,snx);C=(c1,c2,…,cnx)。

在 GNSS计算中,通常是以 ITRF作为参考框架获取速度场,而在实际分析时为了突出相对变化,往往从中剔出欧拉运动,即相对于区域整体无旋转基准[24]。实现假定 (λj,φj)为测站 j的位置坐标,ve(λj,φj)、vn(λj,φj)为其相应的东向和北向相对运动,Voe(λj,φj)、Von(λj,φj)为相应的欧拉运动。由于欧拉运动是通过解析式计算,故只需对相对运动进行解析即可。在用多核函数进行数值逼近时应将所有测站位置作为核函数的核点位置,因此对于任一方向均可列出 n个方程 (测站的个数),故可求解n个未知数(C)。这时可获得东、北向相对运动的解析式 fe(λ,φ)、fn(λ,φ)。由此即可计算研究区任意位置的相对水平运动结果。即

同理,也可获得相应误差的解析式:

2.2 水平运动的滤波计算

由于上述数值逼近所得到的是水平运动观测值的解析结果,它不仅包含运动信息,同时也包含误差干扰,所以由此得到的运动结果实际上是各种成分的综合结果。因此,进行滤波与信息分离是必要的。多核函数不但可以进行数值逼近,而且也具有滤波功能,故在进行空间滤波时同样利用多核函数法。滤波计算时,滤波效果或滤波层次取决于核点的选用与分布。从性质上看,核点的间隔越小 (不得小于测站间的平均间隔),由此所获得的运动信息所包含的频谱带宽越宽,可涵盖从低频至高频的范围;反之带宽越窄,则越趋于低频。在进行滤波计算时,由于运动场是由东、北两个分量共同描述,所以为了保持运动场描述的协调性和彼此之间的相关性不被破坏,进行两分量的滤波计算时应采用相同的滤波计算准则。故东、北向运动分量的滤波函数均表述为:

式中：AT=(a1,a2,…,anx)、BT=(b1,b2,…,bnx),均为待定系数;ST=(s1,s2,…,snx)为核函数阵;si(λ, φ,λi,φi)=exp(-(di/780)2.8)(di为球面上两点间的大地线长度,以 km为单位),(λi,φi)为核点位置坐标。

在实际计算时,首先以相同的网格化形式 (等步长)获取两个分量的数据,并以此作为“观测值”。然后,根据最小二乘法求解上述待定系数,即对每个分量通过误差方程

建立相应的法方程

并由此求解 X(待定系数)：

式中,N=DTPL、W=DTPD、L为相应分量观测值的矩阵,P为相应分量观测值的权矩阵 (pi=1/m2i,mi由式(5)获得)。当 X求得后,式 (6)变为经滤波后运动计算的函数解析式,具体为:

据此可获得任意位置水平运动滤波值。

根据误差传播律,可进行滤波值的误差估计。由式(7)可计算相应运动分量的单位权方差估值,东向为:

北向为 :

式中,nx为核点的个数,n为测站数量,且大于 nx。式(11)、(12)中的残差均是测站点位置的相应残差。令

式中 j=ee,nn。东向运动方差为:

北向运动方差为：

由式(14)、(15)及测站东、北分量间的相关系数即可计算滤波运动的误差椭圆。

2.3 应变场的计算

在球面上进行应变计算应参照石耀霖、刘序俨[18,19]所给出的应变算式。但由于其中存在非微分项的运动,所以其大小必然与运动的参考基准有关。关于该问题杨国华[17]指出,只要依据 ITRF参考框架所得到的水平运动 (位移)即可获得无偏应变计算结果。由于文献[18,19]所给的算式是依据余纬度而建立的算式,此外其坐标系的使用也与现行的不一致,故使用起来不太方便。为此,进行变换是必要的。假定东西向应变为εe(λ,φ),南北向应变为εn(λ,φ),它们之间的剪应变为εen(λ,φ)。滤波后在 ITRF框架下运动的解析式为:

其中

在现行球面坐标系统下球面应变算式为:

式中R为地球的平均半径。

在椭球面上为:

式中,Rn为卯酉圈曲率半径,Rm为子午圈曲率半径,λ为经度,φ为纬度。

最大主应变为:

最小主应变为:

最大剪应变为:

面应变为:

最大主应变方向为:

3 算例

以 2003—2007年华北地区 GPS流动测站在ITRF参考框架的平均水平运动速度结果为例进行分析。为便于观察,首先剔除其中的欧拉运动,获得相对于区域无旋转基准的运动结果 (图 1),然后再依据式(3)、(4)对研究区相对运动速度及误差结果进行数据逼近获得解析式。在此基础上结合式 (6)及以经、纬步长分别为 0.2°计算网格点上的运动结果,再进行空间滤波计算。作为算例,这里选用的核点经、纬步长分别为 1.0°和 2.0°,依据式 (6)～(10)进行计算,结果见图 2～3。

图 1 未经滤波处理的相对区域无旋转基准的水平运动结果(误差椭圆置信度为 95%)Fig.1 Horizontal movements without filtering and relative to no-rotation datum (confidence of error ellipse is 95%)

图 2 步长为 1.0°核点间隔的空间滤波(a)(误差椭圆置信度为 95%)与残差结果(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是未经处理的相对区域无旋转基准的相对水平运动。众所周知,华北地区是一个形变信息较弱、干扰较强的地区。从单站运动误差的角度来看,大部分测站的运动信息都在限差范围内。如果由此判断运动场没有比较显著的信息可能会犯错误。其一,根据 GNSS计算软件的处理过程及后处理来看,测站间误差的关联性并非完全恰当(后处理时常常难以顾及前期处理的协方差矩阵);其二,如果完全由误差所致,则运动场结果不太可能存在某种有序的变化。所以对观测资料作进一步处理,以削弱噪声干扰、凸显有序的运动信息是必要的。图 2(a)是核点间隔步长为 1.0°滤波处理结果及相应的误差椭圆。从误差的角度来看,图 2(a)所展现的运动误差明显减小,从运动的有序性来看,运动的趋势性变得较为明朗,因此可使我们比较容易地把握其运动的基本特征。图 2(b)是与此滤波相应的残差结果,不难发现它具有典型的噪声特点。仅此已初步表明,本文的方法是可行的,得到的结果也是有效的。为了获得更多的认识,给出核点间隔步长为 2.0°滤波处理结果、相应的误差椭圆及残差 (图 3)。由图3(a)可知其信噪比更高,有序变化的特征更加突出。从残差的分布也可看出,噪声的特性依然突出(其中或多或少地包含着一些高频成分)。然而,在实际分析中什么样的结果为最佳?难以给出确定的结论,因为它与研究的具体对象或目标有关,故不能一概而论。从滤波的角度来看,核点间隔越长所体现的运动信息越趋于低频。

图 3 步长为 2.0°核点间隔的空间滤波(a)(误差椭圆置信度为 95%)与残差结果(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)核点间隔的东向应变率结果Fig.4 Space filtering results of the east strain rate by kernel pointwith step-size 1.0°(a)and 2.0°(b)

由于滤波结果是通过模型解析式的计算获得,所以根据此解析式即可非常便利地计算与其相应的应变场。图 4所展现的是在上述两种核点间隔情况下的东向应变率结果(图中黑点为 GNSS测站,蓝点表示一些城市的位置,灰粗线为该区主要活动断层,黑粗线为海岸线)。相比之下,图 4(a)所包含的频谱范围较宽,较高频成分相对突出,曲面虽光滑但起伏相对较大;其数值较大,一般为 10-8/a,梯度也较大。图 4(b)的趋势性变化特征则显现得更加清晰,但变化较为宽缓,数值也较小,一般为 10-9/a。由此可进一步看出,核点间的间隔越小,包含不同频谱形变信息就越全面;核点间隔越大,较高频成分被滤掉的就越多。一般地说,进行滤波通常主要考虑两个方面的问题,即噪声削弱与不同层次的形变信息分离。从这层意义讲,图 4(a)主要体现的是对噪声削弱后的处理结果,因此不同频谱形变的成分也得以较全面地保留。但这样的结果有时难以把握趋势的特征性,或与研究目标所需的某种信息显得不够突出。图 4(b)主要显示较低频的形变信息,趋势性变化较为清晰,较高频的信息已被削弱。这说明可通过设定核点的不同间隔来满足不同的需求。若以研究趋势性变化为目标,则图 4(b)的结果为最佳。图5是核点步长为2.0°的主应变计算结果,由此看出主应变方向在空间上不是恒常不变的,这不但表明了介质并非各向同性,同时也展现了主应变场在空间上的变化趋向,可使人们更加细致地了解其变化特征。

图 5 步长 2.0°核点间隔的主应变率结果Fig.5 Space filtering results of the main strain rate by kernel pointwith step-size 2.0°

4 结语

综上所述,本文通过多核函数法实现了运动场的解析和滤波计算,同时也给出了运动场空间滤波结果严密精度评定算式。算例表明：本文方法具有降噪与信息分离的功能。此外,由于滤波所依据的是可解析模型,所以又可非常方便地获得相应的应变场结果。但是,本方法更加完美的使用还需做进一步的工作,如不同情况下核函数的最优选用,网格的划分、核函数点的选择与滤波和信息分离的关系等。

致谢 感谢杨国华研究员的指导!

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 顾国华,张晶.中国地壳运动观测网络基准站 GPS观测的时间序列结果 [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连续观测资料进行强震危险性预测的探讨[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 黄立人,王敏.中国大陆构造块体的现今活动和变形[J].地震地质,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观测结果推导中国大陆现今水平应变场[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初步结果揭示的中国大陆水平应变场与构造变形[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 杨国华,江在森,王敏.印尼地震对我国川滇地区地壳水平活动的影响[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 张希,等.用最小二乘配置获得地形变应变场动态图像的几个问题研究[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 杨国华,等.汶川地震对华北地区水平形变场影响及有关含义的讨论 [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 杨国华,等.华北地区的水平运动场与昆仑山 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 王敏,等.现今中国大陆地壳运动与活动块体模型[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 王琪,等.天山现今地壳快速缩短与南北地块的相对运动[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 李延兴,等.中国大陆活动块体的运动与应变状态[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 李延兴,等.中国大陆及周边地区的水平应变场[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 杨博,等.断裂区带变形分析方法及应用[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 刘序俨,黄声明,梁全强.旋转托球面上的应变与转动张量表达 [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 黄立人,刘天奎.用速度面拟合法研究华北部分地区的现今地壳垂直运动[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 杨国华,黄立人.速度面拟合法中多面函数几个特性的初步数值研究[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 周光文,黄莜容.一种改进的多面函数法[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 杨国华,等.中国大陆整体无净旋转基准及其应用[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级地震中长期危险性预测研究”;“强化华北地区强震监视跟踪”;“十一五”国家科技支撑课题(2006BAC01B03-01-01)

杨博,男,1985年生,主要从事 GPS应用与数据处理研究.E-mail：boyyang1234@qq.com

P207

A