地球物理学报  2016, Vol. 59 Issue (5): 1585-1595   PDF    
基于多尺度球面小波解算GPS应变场的方法及应用
苏小宁1,2, 孟国杰1, 王振3    
1. 中国地震局地震预测研究所, 北京 100036;
2. 同济大学测绘与地理信息学院, 上海 200092;
3. Earthquake Research Institute, University of Tokyo, Tokyo Japan
摘要: 提出了利用多尺度球面小波解算GPS应变场的方法,该方法通过建立不规则分布的小波基函数以凸显测站非均匀分布的特征,获取不同空间尺度的应变场.通过对模拟数据进行对比分析验证该方法的有效性,模拟结果显示,多尺度球面小波应变计算方法稳健性较好.在位移场中加入高斯随机误差(均值为0.0 mm,标准差为0.5 mm)对应变计算结果的影响可忽略不计.从原有位移场随机抽取的60%的数据解算应变场时,仍然能够获取比较可靠的结果,但如果出现大范围的数据空缺,则会对应变结果产生明显影响.基于华北地区GPS测站原始观测数据,通过高精度数据处理方法解算了该区域GPS速度场,根据多尺度球面小波方法,解算了该区域的应变率场及其误差,并分析了其空间分布特征.
关键词: 球面小波     多尺度     GPS应变场     跨断裂GPS剖面     模拟分析    
Methodology and application of GPS strain field estimation based on multi-scale spherical wavelet
SU Xiao-Ning1,2, MENG Guo-Jie1, WANG Zhen3    
1. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
2. College of Surveying and Geo-informatics, Tongji University, Shanghai 200092, China;
3. Earthquake Research Institute, University of Tokyo, Tokyo Japan
Abstract: A methodology was proposed to estimate spherical GPS strain based on multi-scale spherical wavelet. By constructing an irregular distribution of wavelet basis functions to highlight non-uniform distribution of stations, this methodology can be used to calculate the strain fields at different spatial scales. The methodology is verified and evaluated using synthetic displacements. The simulated results show that the proposed methodology is robust. If Gaussian errors, with a mean value of 0.0 mm and a standard deviation of 0.5 mm are introduced into displacements, its effect to derived strain field is negligible. If 60% data are extracted from the original displacement field to calculate its strains, the obtained strain field is still reliable. Nonetheless, if large scale of data is missed, it will bring noticeable effect to the derived strain field. Based on the GPS raw data in North China, we obtain a displacement field utilizing a precise data processing strategy. Subsequently, we calculate strain rate field and its associated uncertainty, which spatial characteristics are finally analyzed.
Key words: Spherical wavelet     Multi-scale     GPS strain fields     Cross-fault GPS profile     Simulation analysis    
1 引言

随着全球范围大量GPS测站的布设与加密,GPS观测资料在研究现今地壳形变及其动力学机制等方面得到广泛应用(Gan et al.,2007; Meng et al.,2008; Liang et al.,2013Wang et al.,2011).由于GPS速度场建立于某个参考基准,例如,全球参考基准,区域参考基准,选择不同的参考基准会导致速度场存在较大差异,而GPS应变场与基准选择无关.因此,地壳应变场的解算在地壳形变和断层活动性研究方面具有重要地位(Poyraz,2015; Pearson and Snay,2014; Craig and Calais,2014).

GPS应变场解算主要有两种模型:一种是根据断层几何和运动特征建立块体运动学模型,计算块体内部应变场,称为物理模型;另一种是直接利用GPS速度场进行计算,称为数学模型(Tape et al.,2009; 葛伟鹏等,2013Wu et al.,2011),目前已经提出的数学模型主要有最小二乘配置、距离加权、多面函数、张力样条、Kriging差值等(武艳强等,2009江在森和刘经南,2010Shen et al.,1996Gan et al.,2007刘晓霞等,2014).由于解算方法不同,即使利用同一数据也难以得到相同的应变结果(武艳强等,2009江在森和刘经南,2010),因此如何把实际资料反映的不同空间尺度、不同时间尺度的地壳运动与变形信息客观地表达出来,仍然是值得探究的问题.不同于其他数学模型,通过建立不规则分布的小波基函数,多尺度球面小波方法能够充分体现GPS测站不均匀分布的特征;同时也可通过选择不同的尺度因子,研究不同空间尺度的GPS应变场.本文分析了多尺度球面小波原理,探讨了利用多尺度球面小波解算GPS应变场的方法,通过模拟断层闭锁积累的应变,验证了该方法的有效性,并利用本文提出的应变计算方法研究了华北地区不同空间尺度的GPS应变场分布特征.

2 多尺度球面小波解算GPS应变场 2.1 多尺度球面小波基本原理

球面小波是小波理论从无限的平面空间到有限的球面空间的延伸,是小波分析理论的前沿研究分支,已在地磁场和重力场研究中得到了广泛应用(马涛等,2003Holschneider et al.,2003; Chambodut et al.,2005),同时也可用于建立GPS地壳运动速度场和应变场模型.首先定义在单位球上的小波母函数,然后通过作用在母函数上的仿射变换(平移、伸缩)得到小波基函数.任何二维平面上的小波母函数通过逆映射变换都可以得到其球面上的表达式(Bogdanova et al.,2005),所以平面上的小波母函数都有其球面的逆映射形式.本文选择用高斯函数的差(Difference of Gaussians,DOG)(包伯成等,2009)形成的球面DOG小波作为基函数,其表达式如式(1)所示:

式中a=1/2qq为尺度因子,γ是观测点与小波基函数所在位置之间的夹角,单位为弧度,α是1~2之间的常数,取1.25,.由式(1)可以看出:(1)DOG小波基函数满足在有效支撑区域之外快速衰减的要求;(2)不同尺度下DOG小波的衰减特征存在差异,尺度因子越大则衰减越快,其影响范围越小;(3)不同尺度因子的基函数对同一观测站的影响不同.

小波母函数的伸缩使得小波方法具有不同的分辨率.对于单个测站而言,不同时刻的数据代表时间分辨率,而同一时刻不同位置的测站则代表空间分辨率.将待处理信号用球面小波的方法在不同尺度上进行分解,进而得到对应的平滑信息和细节信息.

2.2 分解尺度的选择

通常GPS测站在空间上分布不均匀,直接用GPS速度场解算应变率场时不同测站的权重不同,导致在观测点密集的地方权重较大,故权重的选择对解算结果产生很大的影响(江在森和刘经南,2010).分布不均匀的GPS测站也可视为不同间隔的空间序列,测站的疏密程度代表了不同的空间尺度.通过选择建立不规则分布的小波基函数,多尺度球面小波方法充分考虑了测站不均匀分布的特征.选择分解尺度的原则是:在测站较密集的区域分解尺度较大,在测站较稀疏的区域分解尺度较小.如果选择等边三角形作为图形单元表达整个球面,那么不同尺度下三角形的平均边长如表 1所示(Wang and Dahlen,1995),设立基函数需要满足的条件是在图形单元覆盖的区域至少包含3个GPS测站.当只包含3个GPS测站时,最小图形单元对应的尺度因子为最大分解尺度,而最小分解尺度可根据观测网的范围确定.

表 1 球面三角形网格的几何特征 Table 1 Geometric property of the spherical-triangular grids

图 1为实际分布的GPS测站所确定的最大分解尺度.从该图可以看出,只有两个小区域能够达到测站间距为6.89 km的10级分解,测站间距为13.78 km的9级分解较为零散地分布于整个区域,测站间距为27.55 km的8级分解覆盖了大部分区域,整个区域最小分解尺度可达到6级,所以小波分解的最大分解尺度可以选择6~9级,从而在不同区域获取不同空间分辨率的应变信息.

图 1 由GPS站分布确定的球面小波最大分解尺度 Fig. 1 Maximum resolved scale of spherical wavelet by the locations of GPS sites
2.3 模型参数的解算

在球面上GPS速度场可表示为

式中,依次表示垂直向、南北向和东西向测站速度的3个方向,θ为测站纬度,Φ为测站经度.利用地球表面不均匀分布的GPS站速度v(θΦ)构建小波基函数,进而描述不用空间尺度的地球物理现象和不同噪声源的特征.

已知任何平方可积空间的标量函数f(fL2(S2))可以写成

因此,利用式(3)将式(2)的分解问题转换为3个标量函数的求解问题.把式(3)带入式(2)可得

式中gk(θφ)为DOG小波基函数,由式(1)确定,v(θφ)为实际观测的GPS测站速度,âkĉk为待估计的模型参数.由于gk(θφ)为非正交基,使得法方程系数矩阵秩亏,导致待估参数的解算具有不唯一性,本文采用最小二乘岭估计和误差传播定律解算模型参数及其误差.

2.4 应变场及其误差解算

Savage等(2001)孟国杰等(2008)均给出了球坐标系下三维位移与应变的微分表达式.考虑到地球自由表面的应变计算实质仍然是二维问题,而GPS测站不在同一球面,因此在利用三维应变计算公式解算二维球面应变时需要将不在同一个球面的测站速度投影到同一球面,但由于参与应变计算的测站很难满足同量值的垂向运动,因此如何投影是非常复杂的问题;此外,相对水平形变产生的应变,同量级的垂向形变对水平应变的影响可忽略不计.基于此两点,本文暂不考虑垂向速度对水平向应变率的影响,此时应变张量计算如式(5)所示:

式中为速度梯度张量,r是地球半径,θ是余纬,φ是经度.值得注意的是,式(5)定义的坐标系以南向为正,而GPS速度场通常以北向为正,在计算时需要对不同坐标系的速度和应变张量进行坐标旋转.

依据模型参数,根据应变基函数可以计算速度梯度张量,计算方法如式(6)所示,该式中G(θφ)=,.再由式(5)可以解算出任意目标点的应变张量,进而推导出需要的应变参数.已知模型参数的误差,则可由误差传播定律推导出应变张量的误差.

3 应变计算方法有效性验证

利用位错模型分别正演走滑、逆冲和倾滑三类断层在地表产生的位移场和应变场,然后采用球面小波方法,由位移场数据计算应变张量.Baxter等(2011)认为走滑断层在断层附近能够引起与计算方法无关的应变失真现象,为检验本文应变计算方法的有效性,仅对逆冲和倾滑两类断层进行分析.限于篇幅,本文仅展示对倾滑断层的分析结果.

由位错模型正演的应变称为“理论应变”.由正演的位移场,利用本文的应变计算方法解算的应变为“计算应变”.通过“理论应变”和“计算应变”对比,分析本文应变计算方法的有效性.对于逆冲和倾滑断层,均考虑位移场加入随机误差和不加入随机误差、测站空间均匀分布和不均匀分布等不同情况,其目的是研究本文应变计算方法对误差的敏感程度,以及数据分布不均匀性对应变计算结果的影响.其中增加误差的方法是在正演的地表位移场中加入均值为0.0 mm,标准差为0.5 mm的高斯随机白噪声.产生测站空间不均匀分布的位移场的方法有两种,一是从正演的均匀分布的位移场中随机抽取其中的60%,二是在该区域有GPS测站的位置正演位移场(图 2).

图 2 由位错模型正演的位移场
(a) 原始位移场; (b) 加入误差的位移场; (c) 加入误差并随机提取其中60%的位移场; (d) 实际GPS测站的位移场.
Fig. 2 Displacement fields from dislocation model
(a) Original displacements; (b) Displacements with stochastic errors added; (c) Randomly extractd 40 percent of displacements, with stochastic errors introduced; (d) Displacements at real GPS sites.
3.1 基于位错模型正演的位移场和应变场

设定一条倾角为45°的倾滑断层,断层起点为(115°E,41°N),终点为(120°E,40°N),宽度为40 km,假定地表至地下15 km处完全锁定,15 km以下存在20 mm的逆冲和20 mm左旋的滑动量.利用均匀各向同性半无限空间的位错模型(Okada,1985),计算由于断层闭锁在地表产生的采样间隔为 0.25°×0.25°位移场和应变场.图 3为倾滑断层位错在地表产生的主应变/面应变和最大剪应变.为显示清晰,对主应变进行了重采样,从图中可以看出梯度较大的应变场主要分布在断层附近.

图 3 由位错模型正演的应变场(单位:10-8)
(a) 主应变和面应变; (b) 最大剪应变.
Fig. 3 Displacement and strain fields from dislocation model (unit: 10-8)
(a) Principal strain and dilatation strain fields; (b) Maximum shear strain fields.
3.2 利用多尺度球面小波方法解算的应变场

在构建DOG小波基函数时,需要通过对尺度因子增加比例系数,以调节基函数的强度和覆盖范围,使得应变场既能够分布光滑,也能够突出小区域的细节,该系数类似于高斯距离加权法的平滑半径和张力样条的张力参数.Tape等(2009)认为该系数的选择范围应为1~2,通过多次重复试算,本文最终选择1.25.

利用正演的位移场,由本文的应变计算方法解算应变场,结果如图 4所示,小波尺度因子选择3~9.与图 3对比可以看出,该结果与理论应变在应变梯度比较大的区域几乎一致,在离断层较远的区域存在微小差异,因此,可以认为二者在空间分布特征和量值上都非常接近.

图 4 利用多尺度球面小波方法解算的应变场(原始位移场数据,单位:10-8)
(a) 主应变和面应变; (b)最大剪应变.
Fig. 4 Strain distribution using multi-scale spherical wavelet (pure displacement fields from dislocation model, unit: 10-8)
(a) Principal and dilatation strain; (b) Maximum shear strain.
3.3 数据误差对解算应变的影响

考虑到GPS测站速度是含有误差的,因此在正演的位移场中加入均值为0.0 mm、标准差为0.5 mm的高斯随机白噪声,以研究本文应变计算方法对于误差的敏感程度.加入误差的位移场如图 2b所示,从该图可以看出,相对于原始位移场,加入误差之后的位移场一致性稍差,但分布形态并未改变.以加入误差的位移场作为输入数据,由本文应变计算方法计算应变场,结果如图 5所示.对比图 3图 5可以看出,二者的主应变/面应变和最大剪应变分布特征均十分相似,量值也基本一致,可以认为二者非常接近.说明本文应变计算方法对于高斯随机分布的偶然误差不敏感,方法具有较好的稳健性.

图 5 基于多尺度球面小波方法解算的应变场(位移场加入误差,单位:10-8)
(a) 主应变和面应变; (b)最大剪应变.
Fig. 5 Strain and displacement distribution using multi-scale spherical wavelet (displacement fields with stochastic errors added, unit: 10-8)
(a) Principal and dilatation strains; (b) Maximum shear strain.
3.4 测站空间分布对解算应变的影响

GPS观测网的测站分布通常是不规则的,因此有必要分析测站空间分布的不均匀性对本文应变计算方法有效性的影响.从图 2b位移场中随机剔除其中40%的数据,剔除后的位移场如图 2c所示.可以看出剔除之后位移场并未出现区域性的空缺.由图 2c所示的位移场,利用本文应变计算方法计算的应变场,结果如图 6所示.对比图 3图 6可以看出,二者的主应变/面应变和最大剪应变分布形态仍然十分相似,虽然图 6中应变的量值仅为图 3的85%,但仍然可以认为二者非常接近.说明利用60%的位移场数据能够恢复真实的应变场,表明本文应变计算方法在数据比较稀疏时仍然能够获得好的应变解算结果.

图 6 基于多尺度球面小波方法解算的应变场(随机提取60%的位移场并加入误差,单位:10-8)
(a) 主应变和面应变; (b)最大剪应变.
Fig. 6 Displacement distribution by multi-scale spherical wavelet approach (displacement fields added by stochastic errors and deleted to 60%, randomly, unit: 10-8)
(a) Principal and dilatation strains; (b) Maximum shear strain.

为进一步分析数据空间分布对应变计算结果的影响,在该条断层周边区域存在GPS测站的位置正演了位移场,并加入均值为0.0 mm、标准差为0.5 mm的高斯随机白噪声,该位移场如图 2d所示,由此图可以看出,在断层北侧测站比较稀疏,平均间距约为120 km;断层东东南侧存在一个400 km×400 km的空区,其他区域分布较为均匀.由图 2d所示的位移场,利用本文应变计算方法计算的应变场,结果如图 7所示.对比图 7图 3可以看出,二者的应变分布特征有相似性,基本反映了断层的活动特征,但在量值上图 7中的应变仅为图 3的63%,说明二者存在较大差异.该结果似乎与图 7的结论存在矛盾,但考虑到应变本质反映的是形变梯度,因此当测站间距跨越了变形区域,形变就可能会被忽略,其计算的应变也会与真实结果存在较大差异.

图 7 多尺度球面小波方法解算的应变(真实测站位移场并加入误差,单位:10-8)
(a) 主应变和面应变; (b) 最大剪应变.
Fig. 7 Strain and displacement distribution using multi-scale spherical wavelet (displacements with stochastic errors introduced in real GPS sites, unit: 10-8)
(a) Principal and dilatation strains; (b) Maximum shear strain.
4 华北地区应变率场分布特征

将多尺度球面小波应变计算方法应用于华北地区(其范围如图 8所示)实测GPS站速度数据分析.测站由两部分组成:(1)跨张-渤断裂和郯-庐断裂布设的6个GPS连续观测站及2条区域站剖面,其中连续站如表 2所示,区域站每年观测两期,在2008-2012年期间进行了9期观测;(2)中国构造环境监测网络(简称“陆态网”)GPS基准站和复测区域站.以上共计28个连续站,308个复测区域站.考虑到华北地区受到2011年日本东北大地震同震(王敏等,2011Cheng et al.,2014),以及震后余滑和松弛的持续性影响,为保证速度场的精度,在分析时未选择日本东北大地震之后的观测资料,同时剔除观测周期不足两年的测站,最终得到220个站的速度,扣除区域整体运动后,速度场结果如图 8所示,从该图可以看出,华北地区的运动特征在不同区域存在差异性.在张-渤断裂带以北,GPS速度场一致性比较好,在其以南比较杂乱.

图 8 华北地区GPS速度场(无整体旋转基准) Fig. 8 GPS velocity fields in North China with respect to No-Net-Rotation reference frame

表 2 跨断裂GPS连续观测站 Table 2 Continues GPS sites of cross fault

图 9给出了不同空间尺度的面应变率,约定压缩为正,膨胀为负.从图 9a9d空间解析度依次增加,对应的尺度因子分别为6~9,其中图 9a要求最小站点距离小于110.2 km,最大面应变率为1.7×10-8/a;图 9d要求最小站点距离小于12.8 km,最大面应变率为5.3×10-8/a.可以看出不同空间尺度应变率在一些区域存在差异,小的分解尺度得到的应变率场反映出大尺度应变积累特征,大的分解尺度得到的应变率场能够反映出小区域的应变积累特征,说明本文应变计算方法具有依据站点密度计算不同空间尺度应变的优点.

图 9 基于多尺度球面小波方法解算的不同尺度GPS面应变率场(单位:10-8/a) Fig. 9 GPS strain rate fields in different scale calculated by multi-scale spherical wavelet (units: 10-8/a)

图 10给出了最大分解尺度为8时的主应变/面应变率和最大剪应变率,从该图可以看出,主应变率显示该区域的形变场主要受控于近东西向的挤压;面应变率呈现压缩和拉张交替出现的特征,其中廊坊—天津一带呈面压缩,张-渤断裂带面压缩和面膨胀交替出现.最大剪应变率的分布显示,该区域存在两个最大剪应变率较大的区域,分别分布于沿着张-渤断裂带,保定与沧州之间.该结果与Wu等(2006)杨国华等(2013)的应变率场结果在大空间尺度上比较一致,但在局部区域存在差异.可以看出,利用本文应变计算方法可以获取小空间尺度的应变积累特征.

图 10 华北地区GPS应变率场(单位:10-8/a)
(a) 主应变率和面应变率; (b) 最大剪应变率.
Fig. 10 GPS strain rate fields in North China (units:10-8/a)
(a) Principal and dilatation strain rates; (b) Maximum shear strain rate.

图 11为多尺度球面小波方法计算的面应变率中误差,从该图可以看出,面应变率的中误差均小于1×10-8/a,如果取2倍中误差作为面应变率可以分辨的阈值,那么图 10中几个典型区域应变结果均是可信的.

图 11 华北地区GPS面应变率场中误差 Fig. 11 Standard deviation of GPS strain rate field in North China
5 结论

如何由地表大地测量得到的位移、基线变化等观测量获取地壳应变一直是地学研究中的热点问题,本文将多尺度球面小波方法应用于GPS应变场解算,通过对模拟数据的分析,验证了该方法的有效性.多尺度球面小波应变计算方法稳健性较好,对数据中加入的高斯随机误差反映不敏感,随机误差对计算结果的影响基本可以忽略.从规则分布的位移场数据随机抽取其中的60%解算应变场,结果显示,即使测站分布比较稀疏,该方法仍然能够获取比较可靠的结果;但如果位移场数据出现大区域的空缺,则会对应变结果产生显著影响.

基于华北地区新布设的跨断层GPS连续观测站和区域剖面站,以及“陆态网”GPS基准站和复测区域站数据,解算了该区域的GPS速度场,并利用本文的应变计算方法计算了华北地区的应变率场及其误差.结果显示,廊坊—天津一带面应变呈压缩特征,张-渤断裂带呈压缩和膨胀交替特征,该区域存在两个最大剪应变率较大的区域,分别分布于张-渤断裂带和保定与沧州之间.此外,不同空间尺度的应变率存在差异性,利用本文提出的应变计算方法可以获得小空间尺度的应变积累特征.

致谢    中国构造环境监测网络国家数据中心提供了部分原始GPS观测数据.部分图形采用GMT软件绘制.伍吉仓教授和Carl Tape博士对本研究提出了许多宝贵意见,在此一并表示感谢.

参考文献
[1] Bao B C, Hu W, Liu Z, et al. 2009. Dynamical analysis of DOG wavelet mapping with dilation and translation. Acta Physica Sinica (in Chinese), 58(4): 2240-2247.
[2] Baxter S C, Kedar S, Parker J W, et al. 2011. Limitations of strain estimation techniques from discrete deformation observations. Geophys. Res. Lett., 38(1): L01305, doi: 10.1029/2010GL046028.
[3] Bogdanova I, Vandergheynst P, Antoine J P, et al. 2005. Stereographic wavelet frames on the sphere. Appl. Comput. Harmonic Anal., 19(2): 223-252.
[4] Chambodut A, Panet I, Mandea M, et al. 2005. Wavelet frames: an alternative to spherical harmonic representation of potential fields. Geophys. J. Int., 163(3): 875-899.
[5] Cheng J, Liu M, Gan W J, et al. 2014. Seismic impact of the MW9.0 Tohoku earthquake in Eastern China. Bull. Seismol. Soc. Am., 104(3): 1258-1267, doi: 10.1785/0120130274.
[6] Craig T J, Calais E. 2014. Strain accumulation in the New Madrid and Wabash Valley Seismic Zones from 14 years of continuous GPS observation. J. Geophys. Res., 119(12): 9110-9129, doi: 10.1002/2014JB011498.
[7] Gan W J, Zhang P Z, Shen Z K, et al. 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements. J. Geophys. Res., 112(B8): B08416, doi: 10.1029/2005JB004120.
[8] Ge W P, Wang M, Shen Z K, et al. 2013. Intersiesmic kinematics and defromation patterns on the upper crust of Qaidam-Qilianshan block. Chinese J. Geophys. (in Chinese), 56(9): 2994-3010, doi: 10.6038/cjg20130913.
[9] Holschneider M, Chambodut A, Mandea M. 2003. From global to regional analysis of the magnetic field on the sphere using wavelet frames. Phys. Earth Planet. Inter., 135(2-3): 107-124.
[10] Jiang Z S, Liu J N. 2010. The method in establishing strain field and velocity field of crustal movement using least squares collocation. Chinese J. Geophys. (in Chinese), 53(5): 1109-1117, doi: 10.3969/j.issn.0001-5733.2010.05.011 .
[11] Liang S M, Gan W J, Shen C Z, et al. 2013. Three-dimensional velocity field of present-day crustal motion of the Tibetan Plateau derived from GPS measurements. J. Geophys. Res., 118(10): 5722-5732, doi: 10.1002/2013JB010503.
[12] Liu X X, Jiang Z S, Wu Y Q. 2014. The applicability of Kriging interpolation method in GPS velocity gridding and strain calculating. Geomatics and Information Science of Wuhan University (in Chinese), 39(4): 457-461.
[13] Ma T, Li F, Yue J L. 2003. Application of wavelet analysis in geophysics and geodesy. Progress in Geophysics (in Chinese), 18(1): 49-52.
[14] Meng G J, Ren J W, Wang M, et al. 2008. Crustal deformation in western Sichuan region and implications for 12 May 2008 MS8.0 earthquake. Geochemistry, Geophysics, Geosystems, 9(11): Q11007, doi: 10.1029/2008GC002144.
[15] Meng G J, Ren J W, Wu J C, et al.2008. Computation of strain and rotation tensor and their uncertainty for small arrays in spherical coordinate system. Acta Seismologica Sinica,30(1):67-75.
[16] Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am., 75(4): 1135-1154.
[17] Pearson C F, Snay R A. 2014. Strain partitioning along the western margin of North America. Journal of Structural Geology, 64: 67-78.
[18] Poyraz F. 2015. Determining the strain upon the eastern section of the North Anatolian fault zone (NAFZ). Arab. J. Geosci., 8(3): 1787-1799, doi: 10.1007/s12517-014-1301-6.
[19] Savage J C, Gan W J, Svarc J L. 2001. Strain accumulation and rotation in the Eastern California Shear Zone. J. Geophys. Res., 106(B10): 21995-22007.
[20] Shen Z K, Jackson D D, Ge B Z. 1996. Crustal deformation across and beyond the Los Angeles basin from geodetic measurements. J. Geophys. Res., 101(B12): 27957-27980, doi: 10.1029/96JB02544.
[21] Tape C, Musé P, Simons M, et al. 2009. Multiscale estimation of GPS velocity fields. Geophys. J. Int., 179(2): 945-971.
[22] Wang M, Li Q, Wang F, et al. 2011. Far-field coseismic displacements associated with the 2011 Tohoku-oki earthquake in Japan observed by Global Positioning System. Chinese Sci. Bull., 56(23): 2419-2424, doi: 10.1007/s11434-011-4588-7.
[23] Wang Q, Qiao X J, Lan Q G, et al. 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan. Nature Geoscience, 4(9): 634-640, doi: 10.1038/ngeo1210.
[24] Wang Z, Dahlen F A. 1995. Spherical-spline parameterization of three-dimensional earth models. Geophys. Res. Lett., 22(22): 3099-3102.
[25] Wu J C, Tang H W, Chen Y Q, et al. 2006. The current strain distribution in the North China Basin of eastern China by least-squares collocation. Journal of Geodynamics, 41(5): 462-470, doi: 10.1016/j.jog.2006.01.003.
[26] Wu Y Q, Jiang Z S, Yang G H, et al. 2009. The application and method of GPS strain calculation in whole mode using least square collocation in sphere surface. Chinese J. Geophys. (in Chinese), 52(7): 1707-1714, doi: 10.3969/j.issn.0001-5733.2009.07.005.
[27] Wu Y Q, Jiang Z S, Yang G H, et al. 2011. Comparison of GPS strain rate computing methods and their reliability. J. Geophys. Int., 185(2): 703-717.
[28] Yang G H, Yang B, Chen X, et al. 2013. The basic characteristics of spatial variation of 3D deformation field in current North China. Geomatics and Information Science of Wuhan University (in Chinese), 38(1): 31-35.
[29] 包伯成, 胡文, 刘中等. 2009. DOG小波映射伸缩和平移的动力学分析. 物理学报, 58(4): 2240-2247.
[30] 葛伟鹏, 王敏, 沈正康等. 2013. 柴达木—祁连山地块内部震间上地壳块体运动特征与变形模式研究. 地球物理学报, 56(9): 2994-3010, doi: 10.6038/cjg20130913.
[31] 江在森, 刘经南. 2010. 应用最小二乘配置建立地壳运动速度场与应变场的方法. 地球物理学报, 53(5): 1109-1117, doi: 10.3969/j.issn.0001-5733.2010.05.011.
[32] 刘晓霞, 江在森, 武艳强. 2014. Kriging方法在GPS速度场网格化和应变率场计算中的适用性. 武汉大学学报·信息科学版, 39(4): 457-461.
[33] 马涛, 李斐, 岳建利. 2003. 小波分析在地球物理及大地测量中的应用. 地球物理学进展, 18(1): 49-52.
[34] 孟国杰,任金卫,伍吉仓等. 2008. 球坐标系中图形单元应变与旋转张及其误差解算. 地震学报,30(1):67-75.
[35] 王敏, 李强, 王凡等. 2011. 全球定位系统测定的2011年日本宫城MW9.0级地震远场同震位移. 科学通报, 56(20): 1593-1596.
[36] 武艳强, 江在森, 杨国华等. 2009. 利用最小二乘配置在球面上整体解算GPS应变场的方法及应用. 地球物理学报, 52(7): 1707-1714, doi: 10.3969/j.issn.0001-5733.2009.07.005.
[37] 杨国华, 杨博, 陈欣等. 2013. 华北现今三维形变场空间变化的基本特征. 武汉大学学报, 38(1): 31-35.