基于GPS数据的几种地壳形变分析方法

基于GPS数据的几种地壳形变分析方法

张永奇

长安大学地质工程和测绘学院，西安（710054）

摘要：主要介绍了几种在地壳形变分析中的方法，介绍了他们的原理，以及在地壳形变监测分析中的应用，同时也说明了几种方法的优缺点。

关键字：DDA；位错模型；刚体旋转加均匀应变模型；经验正交函数分解方法模型；数值流形法

引言

地震的孕育和发生在本质上是地壳内部的应变能逐渐积累并突然释放的结果。伴随着大震孕育或应变能的显著积累, 岩石圈表层必然会表现出某种形式和量级的地壳形变。基于这样的认识, 地壳形变监测一向是地震监测或地震危险性分析的最重要手段之一[1]。

自20 世纪80 年代以来,VLBI、GPS、SL R 等空间大地测量技术的精度得到大幅度提高,其观测结果已被应用于板块运动、冰后期反弹、火山、地震、地球自转和地球系统内部物质再分布等地球动力学过程的研究。尤其是GPS 观测技术,由于GPS 接收机价格不是特别昂贵,观测比较容易实施,在地壳形变监测中得到了广泛应用。我国自20 世纪80 年代开始,在青藏高原、南北地震带、川滇和华北首都圈等地壳运动活跃区和地震监测重点区布设了一系列GPS 观测网,用于研究这些区域的地壳形变和地震以及判断板块内部块体构造运动模型。这些GPS 区域网中的部分站点已经并入国家连续GPS 观测网或定期复测GPS 观测网,列入长期观测计划。随着连续和分期GPS 网观测的持续,积累的地表运动观测数据越来越多。这是一类观测精度非常高(上百km 范围内可以测量1mm/ a 的水平方向变化) 、分布在监测区域地表上的站点坐标时间序列。利用这类时空观测数据提取精确可靠的地壳形变信息,反演地球动力学机制是现代大地形变测量数据处理的科学目标,为此,很多学者开展了一系列的理论和方法研究,也取得了不少应用成果[3]。

进一步的基于GPS观测数据的地壳运动和形变的分析方法，尚处于研究阶段。目前初步的应用方法大体可以分为两类：1.纯数学方法，包括多面函数拟合法，多项式拟合法，双三次样条函数法，协方差经验函数法等；2.以物理力学为基础的方法，如欧拉矢量法，非连续形变分析法（discontinuous deformation analysis ，DDA）等。前者尽管在数学上是自洽的，但物理力学含有不清，且都是以地壳连续性假设为基础的；后者物理力学意义清楚，且考虑了地壳的非连续性，但描述块体形变不够精细和有效。如DDA对块体形变的描述仅用3个参数，过于简单只适合与小尺度块体系统的变形分析。

本文主要介绍现今在地壳形变分析的一些常用方法，如非连续变形分析方法，位错运动模型，刚体旋转加均匀应变模型，经验正交函数分解方法模型，数值流形法。以及对这些方法的优劣进行比较，在此基础上对一些可能应用在地壳形变分析的新方法进行展望，如小波理论，灰色理论，神经网络等。

1.非连续变形分析(DDA)方法

模拟介质不连续缝的历史可追溯到30年前的Goodman、Taylor和Brekke等教授发展的节理单元。对岩土裂缝的数值计算发展很快,并已在岩石工程中得到广泛应用。Cundall介绍的离散元法现在被广泛应用于节理或块状岩石。两者是用虚拟力来调整滑动和阻止块体重叠的一种方法,有时候可达到稳定。1993年,由石根华博士发展的非连续变形分析是全新的数值

方法。这种方法用位移作为未知数,解平衡方程时则是用与有限单元法中结构矩阵分析相同的方法。非连续变形分析严格遵循经典力学规则,它可用来分析块体系统的力和位移的相互作用,对各块体允许有位移、变形和应变;对整个块体系统,允许滑动和块体界面间张开或闭合。如果知道每个块体的几何形状、荷载及材料特性常数,以及块体接触的摩擦角、粘着力和阻尼特征,DDA 即可计算应力、应变、滑动、块体接触力和块体位移。DDA 方法自提出以后已被应用到许多工程领域。在国内,长江科学院与石根华博士合作研究了三峡工程中的岩石稳定性问题以及清江隔河岩水电工程的岩体稳定性问题。在国际上,DDA 方法已被美国陆军工程师团应用于军事工程的安全坚固性分析。在日本,陈光齐博士用DDA 数值分析成功地解释了日本Toyohama 隧道坍塌的机制。在DDA 相关理论研究方面,石根华博士把数学中的“数值流形”概念引入DDA,进一步发展成为能满足更多工程要求的“数值流形方法”。陈光齐博士在此基础上又进一步发展了“高阶流形方法”。在内部变形和块体边界处理的数值实现方面,很多学者也做过改进工作。

1.1 DDA 方法中块体的位移和变形

在DDA 方法中,块体系统的大位移和大变形是通过分步迭代计算的小位移和小变形累加来实现的。由于每一步都是小位移,因此可以设每一块体在每一步过程中具有常应力和常应变,块体任一点),(y x 的位移),(v u 可用6个位移不变量来表示,即

),,,,,(000xy y x r v u γεε,其中,),(00v u 是块体内特殊点(如块体的重心)),(00y x 的刚体位移;0r 是块体绕转动中心),(00y x 的转动角(以rad 为单位);xy y x γεε,,是块体的法向应变和切向应变。由此则可推导出:

?????????

?????????????????????-------=??????xy y x r v u x x y y x x y y x x y y v u γεε00000000020)(1020)

()(01 （1.1） 可以证明上述位移函数是块体变形的全一阶近似。

非连续变形分析的平衡方程式由总势能最小化原理来建立,即由各种力和应力产生的总势能min =∏来推导，则得到平衡方程式为:

0=?∏?ri d （1.2） 式中,i 代表第i 个块体;ri d 是块体i 的位移变量;r=1、2、3、4、5、6,对应于上述6个位移不变量。方程式0,000=?∏?=?∏?v u ,分别代表作用于块体x 、y 方向上所有荷载和接触力平衡。方程式00=?∏?r 代表作用于块体i 上的所有荷载和接触力的力矩平衡。方程式0=?∏?x ε,0=?∏?y ε , 0=?∏?xy γ代表沿x 、y 方向块体i 上的所有外力和应力的

平衡。式(1.2)的系数由下列微分求得:

6,5,4,3,2,1,,02==??∏?s r d d sj

ri （1.3） 上式中,当i=j 时,则为块体i 的系数元素,由块体i 的材料特性和相关块体的相互作用决定,构成6×6阶对称阵。当i ≠j 时,则为块体系统中块体i 的相关联元素,即由块体i 和块体j 之间的接触所决定,也构成6×6阶阵。把块体系统中所有自身的系数子阵和块体间的相关联子阵叠加,则构成总体平衡方程式为:

?????

???????=??????????????????????n n nn n n F F F D D D K K K K 21211111 （1.4）

式中,每个系数元素ij K 都是6×6阶子矩阵;i i F D ,是6×1阶子矩阵,其中i D 代表块体i 的变形变量),,,,,(654321i i i i i i d d d d d d ;i F 是块体i 上分配给6个变形变量的荷载,它可由下式求出,

6,5,4,3,2,1,0)0(==?∏?r d ri

1.2 DDA 方法中块体系统运动的处理

非连续变形分析发展了一种完善的运动学理论,在块体没有被别的块体嵌入的加卸载条件下,它可以计算出许多块体的大变形结果。非连续变形方法用加约束弹簧的方式,克服了块体系统中块体间的相互侵入问题。数学处理上通过引入“弹簧子矩阵”,使严格的不等式成为线性方程系统。

平衡方程式通过用加刚度很大的弹簧在一个或两个方向上锁住块体的运动。如两个块体在它们之间有接触拉力,在锁定撤离后它们将被分开。总方程式在选择锁定位置时必须反复求解。用此方法,块体间即使有拉力和嵌入,仍可在极短时间内正确选定锁定位置。此外,块体的接触力在非连续变形分析中可精确计算。

DDA 的数学模型采用隐式算法,既考虑到静力学又考虑到动力学。时间步用于静力学和动力学两者,惟一的差别是静力学计算假定每一时间步开始速度为零,而动力学计算则保持前一时间步的速度。在大变形情况下,DDA 方法的每一步由前一步形成的变形了的块体形状和位置开始,而建立和求解的平衡方程式是针对当前块体的几何形状和受力状态。

1.3 DDA 法与有限单元法的比较

DDA 是与有限单元方法相平行的一种数值分析方法。它解算的网格与有限单元类型相似,但所有单元是被事先存在的不连续缝所包围的实际隔离块体,这是非连续变形分析显著超过有限元分析的优点。非连续变形分析在块体边界不是连续体,形成网格单元的块体可以是任意条边的,凸状形或非凸状形的,甚至可以是带孔的多接点的多边形。块体网格不要求块体顶点与另一块体顶点相接触,而有限单元法限定只能用标准形状的单元,且顶点相接。此外,在DDA 方法中,当块体接触时,库仑定律可用于接触面,而联立平衡方程式是对每一荷载或时

间增量来选择和求解的。在有限单元法的情况下,未知数是所有节点的自由度之和;在DDA 法的情况下,未知数是所有块体的自由度之和。从理论观点看,DDA 是有限单元法的广义化。

虽然非连续变形分析类似于离散单元法,但它是更接近于与有限单元法相平行的一种方法,因为它是:①采用总势能最小化来建立平衡方程式;②选择位移作为联立方程式的未知数;③把刚度、质量和荷载的子矩阵采用叠加的方式加到联立方程的系数矩阵中去。

非连续变形分析的块体刚度矩阵比有限单元分析的单元刚度矩阵更加简单。该方法把接触块体的位移锁定以及附加的杆单元应用到有限元分析中去非连续变形分析的特点是:完全的运动学及其数值可靠性、完全一阶位移近似、严格的平衡要求、正确的能量守恒和高计算效率。这一方法可靠的原因是分析非常接近实际,力学现象的数学和数值描述与块体运动相一致。

在不连续情况中,大位移和大变形是更为复杂和重要的问题。当块体移动或变形时,新的块体形状和位置将比在连续力学更为敏感的方式下产生不同的块体接触和影响破坏模式。

1.4 DDA 在现代地壳运动中的应用研究

1.4.1 理论模型研究

DDA 方法起源于岩石力学分析,如何科学有效地把DDA 方法用于地壳运动研究,这既是一个理论问题,也是一个应用问题,要实现这一点必须首先研究DDA 的适用性问题。目前,DDA 方法还仅限于二维平面块体运动的方程,它能解决的问题也仅限于较小范围的二维平面问题,而地壳运动的范围可能大到整个板块乃至全球。因此,首先必须把DDA 方法从平面方程发展到适应于大范围地壳运动的球面方程。通过长时间的探索发现,可由下列途径实现:

1)建立球面上块体任一点的位移和位移不变量之间的函数关系。从球坐标系位移与变形的基本关系式出发,通过设定块体的全一阶近似表达式,可以建立球面上的位移),(v u 与位移不变量),,,,,(000φθθφεεεr v u 的函数关系。

2)在势能最小原理下建立球面上块体的联立方程式。分别求出弹性应力、初始应力、点荷载、体积荷载、惯性力、设定点位移等的势能表达式,再通过微分解出系数子矩阵,叠加到联立方程式中。

3)建立球面上由多个块体组成的块体系统运动学方程。在联立方程式建立后,块体系统中块体之间的相互作用问题必须给以精确的数学描述。首先是块体间不允许侵入,通过在边界上加刚性弹簧可以在很高的精度上满足不侵入的问题。刚性弹簧的添加需要通过块体边界的嵌入距离来判定,由于在球面坐标系中计算侵入距离将很复杂,为了使计算简便、高效,可选用局部平面投影坐标系,把球面问题转化为与平面相类似的问题。

2.地壳块体位错运动模型及其反演

2.1块体位错运动模型

Matsuura 率先提出了一种用块体运动与断层位错运动合成来表述块体边界的地面运动的理论模型,并且运用于解释美国加州大地测量得到的地面移动数据。随后,这种模型广泛运用于板块边界、断层带大地形变观测数据的解释。当断层发生地震时,该模型也可以用于描述地震断层滑动分布造成的地表形变场。图1为断层位错运动模型。图中,XOY 为地面,断层面长为L 、宽为W 、倾角为δ,断层位于(Z ≤0)弹性半无限空间中,其下边缘埋深为d,断层面

上的(负)位错分布记为U1、U2和U3,依次称为走滑位错、倾滑位错和张裂位错分量。根据弹性均匀各项同性半无限空间中的位错理论,只要知道上述断层位错模型的4个几何位置参数和3个断层面上的位错分量,就可以计算地表任意点(X,Y)的位移。

具体运用上述位错模型时,一般采用所谓的块体位错运动模型,它实际上是把块体边界的运动表述成块体间的相对运动叠加上断层(负)位错引起的形变。如图2所示,假定A 地块和B 地块间存在着长趋势的稳定的相对运动AB V ,A 地块与B 地块的边界由断层面F 分开,在断层面的上部Σ的阻碍作用以负位错分布模型来表示。则地面任一点的运动y 表示为:

∑+=f V y AB （2.1）

式中, ∑f 表示断层面上负位错分布引起的地壳形变。式(2.1)中,如果y 是两期观测期间监测点的位移,则AB V 表示对应于观测期间的A 地块与B 地块的相对位移, ∑f 表示负位错引起的监测点位移。如果y 是监测点的平均速度,则AB V 表示A 地块与B 地块的相对速度, ∑f 表示负位错引起的监测点速度。

为了描述断层面上位错的非均匀分布,一般将断层面分成矩形格网,地面点位运动的位错贡献是所有断层面上格网单元位错贡献的总和。上述两个块体的运动可以推广到多个块体及多条断层边界。很显然,这种模型也可以表述地震引起的共震位移,这时只要在式(2.1)中令0=AB V 即可。

图1断层位错模型

图2 块体（负）位错运动模型

2.2 大地测量数据反演方法

反演是相对于正演而言的。在式(2.1)中,利用地表观测位移或速度确定断层位错模型包括块体运动参数就是一个反演问题。通常情况下,式(1)是非线性函数,反演问题的解可能不惟

一。一般通过泰勒展开的方法取一次项,将式(1)转换成线性函数。对于线性反演问题,根据观测数据和模型参数的个数确定采用求解的方法。当观测数据个数多于未知参数个数时,称反演问题超定,通常采用最小二乘法;而当观测数据个数少于未知参数个数时,称反演问题欠定,通常要采用顾及参数先验信息的贝叶斯方法,其本质上是顾及参数先验信息的最小二乘法。由于即使在超定的情况下也需要利用先验值(初值)来线性化观测方程;另外,对模型参数解的一些物理属性如光滑、连续等条件也可以在数学上表达为参数的先验信息,所以后一种顾及参数先验信息的最小二乘法或顾及参数先验信息的迭代最小二乘法是现在大地测量数据反演方法的主流。下面以上面的块体断层运动模型反演为例,简要介绍这种反演方法。

在上面的块体断层位错模型式(2.1)中, ∑f 为每一条断层长度、宽度、深度、倾角及负位错分量的函数。假定把式(1)右端所有待求的模型参数记为x ,则有:

)(x f y y =+ε （2.2） 式中,x 为1?m 阶待求未知数向量,包括所有块体相对运动速度分量及每个断层的几何参数和负位错分量;y 为1?n 阶观测数据向量,它由所有观测点上的相对运动速度分量构成;y ε为1?n 阶观测误差向量。

假设x0表示模型参数的先验信息,即x x x =+ε0 （2.3） 式中,x ε是m 维随机误差向量。通常假定随机向量y ε、x ε均服从高斯分布,且它们的期望为零,方差阵分别为∑∑y x ,。迭代求解过程可以表述为:

k k k k k r M b x x 11-++= （2.4）

这里,

∑∑--+=11x

y k T k A A M （2.5） []∑∑---+-=101)()(x k y k T k x x x f y A r （2.6）

其中,k A 的元素按下式计算: k x x j i ij x f A =??= （2.7）

k b 是一个控制迭代步长的正数,0≤k b ≤1;0x 是初始迭代向量,也就是先验值。取一个小的γ作为迭代阈值,控制迭代终止。于是,迭代问题的解可以定义为:

1+∧=k x x ，如果γγ

111

111)(---+++∑∑+=y x k T k k A A C (2.9) 很显然,对于上述非线性反演问题的解法,其反演结果对初值的依赖性较大,如何给定比较好的初值直接影响解的正确性。另一方面,有关模型参数的先验信息或条件不仅仅只有式(2.3)的表达形式,还可以有其他的表达方式,要根据反演问题的具体情况而定。上面的反演是针对两期观测得到的位移或速度来进行的,把模型参数当作是固定的值。而实际上对于连续

GPS 站观测或多期连续观测,不断地有新的观测数据产生,部分模型参数本身也可能随时间变化,这时就有必要引入一种序贯最小二乘方法或卡尔曼滤波方法来反演时变的模型参数。这类时变的反演方法是目前研究的一个趋势。

3.用刚体旋转加均匀应变模型研究应力应变分布

在研究全球板块运动特征建立板块运动模型时,通常假定板块是刚性的,许多学者做过这方面研究,取得了一定的成果,这也证明板块的

实际上板块不是纯粹的刚体,由于板块间的相互作用,在板块边界上基本上都存在着宽度不等的形变带,在内部一些薄弱地区,板块内部的应力积累会引起这些地区的形变。沿板块边界和内部频繁发生的地震,就是这种形变存在的极好例证。因此,建立了块体的刚体运动和内部的弹塑性形变模型。欧拉矢量对于在球面上描述块体刚性运动方面比较好,用欧拉矢量表示块体的整体旋转运动,并假定块体内部发生均匀应变。模型为:

???????????????????

???????-------=??????en n e z y x n e w w w r r r r r r r r r V V εεεφλλφφγγφφφλλφλφλφcos )()(00cos sin )(0cos )(cos sin sin cos sin 0000 （3.1） 其中, x w ，y w ，z w 为绝对欧拉旋转矢量; e ε,n ε，en ε为应变参数。通过解算可得块体旋转参数及中误差值,块体线应变及中误差值,块体速率残差及运动速率值。其 中,e V , n V ,en V 分别为模型求出的各块体运动速率;e V ? n V ? , en V ? 分别为各块体速率残差。再解出各块体的面膨胀率、最大剪应变率及主应变率大小和方向值,求出主应力的大小及方向值。

4.经验正交函数分解方法模型

4.1 经验正交函数分解的基本原理与计算方法

经验正交函数分解的基本原理与主成分分析基本相同, 如果观测要素是时间、空间变化极为显著的变量, 那么由其构成的要素场的空间分布和时间变化特征就可分解为空间函数和时间函数的线性组合。

设某一区域观测要素场, 有n 个测点, 进行了m 次观测, 其观测值构成矩阵X 为

????????????????????mn mj m m in ij i i n j n j x x x x x x x x x x x x x x x x

2121

222221

111211 (4.1) 经验正交函数分解是把要素场观测序列分解成正交的时间函数和空间函数, 即

L T X ?= (4.2) 其中

??????????=?mm m m n m t t t t T 1111，????

??????=?mn m n n m l l l l L 1111 式中, T 、L 分别为时间函数和空间函数矩阵。把式(4.2) 写成和的形式

∑==m h hj ih ij l t x 1

),,2,1;,2,1(n j m i == （4.3）

式中,hj l j 表示序号为h 的空间函数, 它只依赖于空间变化, 不随时间变化, 也称为序号h 的典型场;ih t 表示序号为h 的空间函数在第i 时刻的权重系数, 它只依赖于时间变化, 不随空间变化, 称为时间函数。在实际应用中, 在式(4.3) 中取有限的几项, 如H 项(H < m ) 就可达到对原观测要素场有较好的拟合。

∑==H

h hj ih ij l t x

1?),2,1;,,2,1(n j m i == （4.4） 式中, ij x

? 为ij x 的估值, H 为典型场的个数。 计算方法如下:

(1) 求协方差矩阵为

X X m

S T 1= （4.5） S 的特征值方程为0=-I S λ （4.6）

(2) 求解特征值方程的特征值λ和对应的特征值向量L , 特征值按降序排列, 即 n λλλ≥≥≥ 21

对应的特征值向量为

????????????=n l l l L 112111 ????????????=n l l l L 222212 ????

?

???????=nn n n n l l l L 21 （4.7） （3）求时间权系数T 为

∑∑===n

j kj n j kj

ij ki l l x t 121 （4.8）

(4) 典型场个数的确定。一般在式(4.3) 右端取前几项就可以得到较好的近似结果, 当k =

H 项时, 整个序列式(4.4) 展开的精度用下式来估计, 即

∑∑===n j jj H k k

S R 11λ （4.9）

根据所要求的精度R 来确定H 值, 一般R 在90% 以上就可确定H 值, 即式(4.3) 应取头几项的项数。特征值向量k L 又称为典型场,L L 依特征值λ按降序排列, 把排在第一位特征值向量1L 称为第一典型场,2L 称为第二典型场, 依次类推。

值得注意的是, 解方程式(4.6) 的方法很多, 一般有迭代法、雅可比法、ivens-Householder 法、QR 法等, 用不同的解法可能得出不同的解, 所以向量L 的解不一定是唯一的, 但是。它们之间往往只相差一个非零的乘常数, 因此L 的值只有相对意义。对不同的数据只有用同一种解法才有可比性。本文用的是Given s-Householder 法。另外, 若矩阵S 是秩亏的, 其秩亏数是r, S 是实对称方阵, 所以S 有r 个零特征值。由式(4.9) 可知, 忽略0=λ 的λ项对拟合精度没有影响, 所以零特征值所对应的特征值向量就可以忽略。

4.2 地壳垂直形变场经验正交函数分解

众所周知, 地壳垂直形变是时间和空间的函数, 那么, 根据经验正交函数分解的基本原理, 某地区的地壳形变场的时间序列构成的要素场就可以分解为时间函数和空间函数。在式 (4.1) 中, 若ij x 为地面上第i 点第j 时段的垂直形变速率(单位为a mm ) , 经验正交函数分解计算方法就变成地壳垂直形变场经验正交函数分解的方法, 此时定义式(4.7) 中的1L 为地壳垂直形变第一典型场,2L 为地壳垂直形变第二典型场。地壳垂直形变典型场与地壳垂直形变场的共同点是它们都包含地壳垂直形变空间变化信息; 不同点是, 前者仅包含地壳垂直形变空间变化信息, 不同的典型场表示由不同的因素引起的地壳形变空间变化特征, 同时地壳垂直形变典型场还能区分不同的因素在地壳形变中所起作用的大小, 而后者既包含地壳垂直形变空间变化信息, 还包含地壳垂直形变时间变化信息。

5.数值流形法

5.1 数值流形方法的思想

数值流形方法是基于有限覆盖技术建立的,有限覆盖技术的基本思想就是在求解区域上

构造一组函数,称为覆盖函数,这种覆盖函数要具有两个基本性质:1) 局部非零性,只在一个局部区域范围内不为0;2) 这组覆盖函数之和在求解区域内恒为1。由于这组函数具有性质

2) ,有时亦称为单位分解函数。数值流形中的覆盖有数学覆盖和物理覆盖,数学覆盖可以任意选择,它由占整个材料体的许多有限重叠覆盖所组成。常规的网格和域,诸如规则的格子、有限元的网格或级数的收敛域,能转换为有限数学覆盖。物理网格包括材料的边界、节理、裂缝及不同材料区域的交界面,不能人为选择。物理覆盖由数学覆盖和物理网格组成。 在研究区域上形成数学覆盖和物理覆盖后,就可以建立有限覆盖位移函数了。对于二维流形方法的计算,设定义在物理覆盖i U 上的位移函数),(y x u i 和),(y x v i ,其中,它们可以是常量、线性的、高阶多项式或局域级数。这些覆盖位移函数可用权函数),(y x w i 连接在一起,并且),(y x w i 应满足：0),(≥y x w i ，i U y x ∈),(，0),(=y x w i ，i U y x ?),(及

1),(),(=∑∈i U y x y x

整个物理覆盖系统上的位移函数可用覆盖位移函数的加权平均表示如下:

[]

{}[]{}

∑∑∑∑======??????=??????n i i i n i ij m j ij n

i i i i D y x T D y x T y x v y x u y x w y x v y x u 1111),(),(),(),(),(),(),(

其中，[]

??????=---),(),(),(),(),(12,212,22,112,1y x t y x t y x t y x t y x T j j j j ij {}??

????????=-j i j i ij d d D 2,12, []()im i i i T T T y x T 21

),(= {}{}T im i i i D D D y x D 21),(=

位移矩阵[]??

????=

=--),(),(),(),(),(),(),(),()(),(2,112,22,21,22,112,12,11,121

y x t y x t y x t y x t y x t y x t y x t y x t T T T y x T m m m m im i i i

其中, i 表示第i 个物理覆盖。

物理覆盖的每一个交集定义为单元,假设单元e 为物理覆盖

)(),2(),1(q V V V e e e 的交集。对任一e y x ∈),(,则0),(≠y x T i 0),(),(=∈y x T r U e i e ， )(r U e e ∈；

[]{}[]

{}∑∑∑∑======??????????=??????q r r e r e q r m

j j r e j r e r e r e q r r e D y x T

D y x T y x v y x u y x w y x v y x u 1)()(11)()()()(1)(),(),(),(),(),(),(),(

有了待定系数的位移,就可以利用和有限元中相同的方法进行求解,即有位移函数可以求出应力、应变,不过它们也都含有待定的常数,求出整个待分析域的总势能,然后根据最小势能原理,列出总体平衡方程式,然后即可由此方程求出其中待定的常数,最终得到问题的解答,这部分和有限元中的分析方法基本相同。

数值流形方法用于地壳形变场分析中的可行性

5.2 地壳变形的力学特征

地壳是由复杂的岩石结构组成,从几何特征看,地壳是由厚度不均匀,且含有节理、裂隙、夹层、断层等不连续面的球形壳体,从物理特性看,岩石是一种分布不均匀、各向异性的非线性介质。从荷载环境看,地壳受到从宇宙形体的吸引力、地幔的浮力、重力和构造应力等的作用,受力情况十分复杂。地壳的运动与变形过程不但与外界的作用力有关,更主要的是与地壳中介质的物理性质和地壳中各种不规则的构造分布有关。如地壳的各向异性、非弹性介质以及地壳中的不连续界面分布对地壳运动与变形过程和地震活动等物理过程都产生非常重要的影响。因此,研究地壳的形变必须要考虑到地壳介质的各向异性,非弹性性质及不连

续面的变化。

5.3 数值流形方法在地壳形变分析中的可行性

数值流形方法建立的基本思路就是建立能解决连续和非连续部分共存特殊介质在外力作用下的力学特征,所以不但能解决模拟连续部分在外力作用下的应力和应变的情况,而且能够模拟裂纹等不连续缝在外力作用下对地壳运动的影响。相对于目前模拟连续介质的有限元和非连续介质的离散元及非连续变形分析(DDA) 相比,数值流形方法具有与实际地壳介质更合理的特点。所以,数值流形方法从一出现就吸引了许多学者的注意,目前数值流形在岩土工程、金属成型等大变形中得到了有效的应用,取得了较好的成果,但是在地壳形变的分析中还没有得到实际的应用。由上面对地壳变形的力学特征的分析中,可知数值流形方法能充分考虑到地壳运动的力学特征,把连续和不连续面更好地融合在一起解决。

有限元方法和非连续变形分析在地壳形变的分析中,已经取得了有效的成果,而有限元和非连续变形是数值流形方法的特例,可以设想数值流形方法也可以用在地壳形变的分析中去,去分析地壳运动的应力应变过程,从而对地壳变形给出科学合理的物理解释。

参考文献

[1] 甘卫军，张锐,等.中国地壳运动观测网络的建设及应用[J].国际地震动态，2007（7）：

43～44.

[2] 彭华,张冬菊.刚体旋转加均匀应变理论在地壳形变中的应用[J].山西建筑,

2007,33(35):121～122.

[3] 陈永奇,伍吉仓.利用GPS监测区域地壳形变的理论与方法[J].武汉大学学报·信息科

学版,2007,32(11):962～963.

[4] 王泽民.非连续变形分析与现代地壳运动研究.武汉大学学报·信息科学

版,2001,26(2):122～124

[5] 吴云,申重阳,周硕愚,等.基于边界元的非连续( 块体系统) 形变反分析法[J]. 武汉

大学学报·信息科学版,2003,28(3):345～349.

[6] 崔笃信, 刘文义.经验正交函数分解在地壳垂直形变场分析中的应用[J].地震,

2000,20(3):82～84.

[7] 张冬菊.数值流形法在地壳形变场分析中的可行性研究[J].山西建

筑,2005,31(17):57～58.

[8] 王敏,沈正康,牛之俊,等.现今中国大陆地壳运动与活动块体模型[J].中国科学（D

辑）,2003,33(增刊)：21～30.

[9] 蒋锋云，王双绪，张希,等.用二维有限元研究青藏块体东北缘地壳的水平运动[J].大

地测量与地球动力学,2007,27(8):23～29.

[10] 许才军, 晃定波.有限元分析中断裂区域单元的划分[J].地壳形变与地

震,1997,17(3):34～38.

[11] 张燕,吴云,刘永启,等.小波分析在地壳形变资料处理中的应用[J].地震学

报,2004,26(增刊)：104～109.

浅谈地形变观测的发展与地震预报

浅谈地形变观测的发展与地震预报 地震预报是一个系统工程，需要长期的科学探索和科学积累，地形变观测系统已成为地震预报学科领域不可取代的力学型的基础观测系统。本文简述了地形变观测的发展历程以及其在地震预报中的运用。 标签：地震预报地形变测量 我国地域辽阔，地质构造复杂，是世界上地质灾害较严重的国家之一。地震就是对人们生活和生命、财产安全影响较大经常发生的一种灾害。地震的前兆是可以预见的。如果在地震前能获知一些前兆特征，便能提醒人们引起注意，从而减少地震造成的灾害。 地震预报是对未来破坏性地震发生的时间、地点和震级及地震影响的预测。 地震预报目前或今后很长一段时间都将是以观测为主的试验性科学。物理统计分析方法是研究地震孕育动力学过程或进行预测的最具有潜力的途径之一。地形变测量就是物理统计分析方法的其中一种。尽管2008年汶川和2011年东日本大地震预报失败，但这两次和其他大地震前后GPS和其他观测得到的地壳形变表明，大地震是有前兆的，是可以预报的. 地震孕育的能量来源于地壳运动过程中产生的应变能积累，地形变观测是监测地壳运动与变形、认识地震孕育过程与开展地震预测的重要手段。如2008年汶川地震后，通过对其地形变化监测，监测显示，汶川地震引起震中区域监测点的水平位移量达238厘米，沉降量达到70厘米，隆起量达30厘米。通过数据观测，得出结果：地震造成灾区地形发生重大变化地形发生重大变化。根据这个结论，开展地形变分析后可以有效预报大地震后的余震及其震级。 地形变测量是指对一个地区地面的相对变化进行的重复地形变测量或连续观测。一般是在某个区域布设相当数量的形变观测点，在区域外设立基准点，利用常规大地测量仪器和工程测量中监测变形的方法。通过定期观测，并通过适当的数据处理和物理解释准确确定区域内各测点的空间随时间变化的“绝对”位移与方向，以达到预测地震之目的。地形变观测目前已具有多种手段并形成一定规模。地震地形变观测系统目前所涵盖的主要观测手段有：垂直形变测量网、水平形变测量网即GPS区域复测、重力测量网、跨断层形变测量网、重力与固体潮观测台网、地倾斜与固体潮观测台网、地应变与固体潮观测台网和连续GPS监测网也叫GPS基准网等。所采用的手段和仪器主要有：水准测量、三角测量、倾斜仪、伸缩仪、电阻丝应变仪、激光测距仪、测潮仪；近些年来还发展到应用一系列空间对地观测技术，如GPS（全球空间定位系统）等。 地形变观测系统所测定的物理量包括：位移、旋转、速度，加速度、应变（应力）、蠕变、位错、重力、固体潮汐、地下介质物性参量（密度、勒夫数）、电离层与对流层介质物性参量（电子浓度、湿度…）等的空间分布及其随时间变化。

地壳形变

大地测量联合反演理论及应用 联合反演模型的研究：赵少荣(1991)系统地研究了“基于固体力学的大地测量反演问题”,给出了基于固体力学模式的多类大地测量反演的解算模型；许才军(1994)给出了大地测量联合反演构造应力场的解算模型; 晁定波等(1997)进一步提出了四维整体大地测量有限单元法,把四维整体大地测量模型(包括GPS,水准,GPS水准和重力监测数据)与固体力学基本方程结合起来,对大地测量和地球物理数据进行有效的整体处理,强化了边界的大地测量约束,具有改善刚度矩阵方程数学性质和降低其阶数的优点,从而提高反演解的稳定性和可靠性。P. Segall等(1997)根据永久GPS网的扩建在时空上为地壳形变量提供了足够的数据,建立了一种网络反演滤波模型,它可以综合频繁收集到的各类大地测量网络数据来估计断层滑动的时空分布,反演获得各种参数,包括观测误差、局部移动、瞬时和空间平滑参数。大地测量反演模型的建立涉及到正演模型和观测数据的类型、反演模型选择及反演参数辨识。正演模型是反演模型的基础,只有清楚正演问题建立正演模型,才能给出反演模型。大地测量反演已由单一的观测数据的反演,发展到多种观测数据的大地测量联合反演问题,这种多种观测数据的大地测量联合反演问题不仅仅指大地测量的多种观测数据,而且也包括不同类型的(地震、地质和大地测量等)多种观测数据。 联合变形、重力和地震资料反演研究地球内部介质参数：地球内部介质的弹性参数μ和λ是利用地表的形变观测资料研究地球的应力场及内部的物质运动所必需的参数。观测资料表明,在地球内部,介质的弹性参数μ和λ存在横向和纵向上的不均匀性。（1）利用变形资料可以研究地球介质弹性参数μ和λ；2）利用地震和重力观测资料可以研究参数μ和λ的不均匀性。由于参数μ、λ与介质的密度、横波速度和纵波速度参数有关,而密度、横波速度和纵波速度参数又是地震观测资料的场源,密度参数还是重力数据的场源,因此,进一步探求参数μ和λ结构的问题,可以转化为利用地震和重力观测资料反演密度、横波速度和纵波速度参数结构的问题。联合变形、重力和地震资料研究地球内部介质参数的联合反演问题可以表示为: BX = UL NX = g SX = d （1）式(1)可以通过下列目标函数取极小值求解:Φ= ( UL-BX)TWL( UL-BX)+(g-NX)T

地壳形变知识点整理

?全球问题：①地球动力现象引起的地震海啸火山等自然灾害，给人类生命财产带来损失②全球气候变暖、海平面上升、局部地层沉降和海上溢油公害等是随着工业发展引起的环境问题③由于人口增加和陆地资源枯竭，需开拓生存空间和寻找新矿产资源。 ?地壳运动：指在地球内部构造应力作用下所引起的地壳一些元素的相对运动。它们可以是⊥运动、水平运动或地倾斜运动，综合表现为大面积的地壳形变。广：地质旋回、狭：构造旋回。分类：水平运动指组成地壳的岩层，沿∥地球表面方向的运动，又称造山运动或褶皱运动；⊥运动又称升降运动或造陆运动，表现为岩层部分区域的隆起和相邻区域的下降。长期运动是在地质时间尺度内的运动，由几千年到几百万年，它与板块运动有关；瞬变运动与地震和火山等活动相联系的。范围：全球板块运动和区域及局部地壳运动。地壳运动监测：测定板块运动参数、大陆板块和海洋板块的内部形变、板块边界与大地震有关的区域形变和局部形变、其他地震活动区的区域形变和局部形变，主要通过建立全球测定板块运动监测网和区域、局部地壳运动监测网实现。 ?地壳形变：指在地球内力和外力作用下，地球的地壳表面产生的升降、倾斜、错动等现象及其相应的变化量。成因：①人类活动，地表形变：离散性、短暂性和局部性； ②地球自转和极移，形变：全球规模特性，可理论计算；③日月等天体对固体地球在引力作用下，形变：固体潮,理论上可严格计算；④大地构造运动（地球内部的构造原因），地壳构造形变：连续性、长期性、区域性、复杂性。地壳形变测量：对一个地区的地壳表面的相对变化进行重复或连续的观测称为地壳形变测量。特点：①以动态观测替代大地测量只以静态观测方式来测定地面点变化，并分析研究其物理意思；②主要在构造带、多震区和具有潜在地震危险的重点地区及在大坝等要害部位进行，而大地测量未考虑这些；③测点设置要求稳定可靠，布网边长短、测量精度高、复测周期密。任务：监测地壳形变运动，具体观测元素是地表点位置变化。种类:①全球板块运动监测：主要用来测定板块运动参数，测定大陆板块和海洋板块的内部形变，VLBI、SLR 、GPS。②全国及区域地壳形变测量：精密水准、高精度流动重力、高精度空间大地。③断层形变测量：短水准、短基线、短边GPS网等。④定点形变测量：有效地监测地壳的连续变动，地倾斜、地应变、重力台站。 全球板块运动监测：意义：①寻找矿产资源，板块运动的边界是生成矿源的地点；②防灾减灾，地震一般发生在板块的边界。方法：VLBI测量、SLR测量、GPS测量等。 ?VLBI：可提供整体运动和地壳运动的丰富信息，短时间测量即可得极高精度。特点：目标源和参考源（点源）距离近；两种源时间间隔小。①纯几何方法，不涉及地球重力场；②不受气候限制，有长期稳定性；③提供以河外射电源为参考的最佳准惯性参考系。?SLR：测定台站的位置变化或站间基线长度的变化率来确定板块运动参数。在一列选定的时间间隔内联合求解测站的坐标、卫星轨道和地球定向参数EOP→测站间基线的时间序列→线性拟合→基线长度变化率。不足：①由于采用可见光，受天气影响大，不能全天候观测；②SLR台站建立和维护费用高。精度：取决于在各所选时间间隔内测得的站间基线精度、基线的观测误差；EOP误差以及卫星的定轨误差均有影响。?GNSS WGS84：世界大地坐标系统，地心系，GPS广播星历。WGS84（G1150）：美国对WGS84第3次精化获得的框架。G框架用GPS资料确定，1150开始用的时间。 ?空间大地测量方法：建立现代板块运动模型。观测值：站坐标和站速度（与参考框架有关）→绝对运动参数，站间基线长度变化率（无关）→相对~。原理：板块构造学说认为相邻两板块之间的相对运动实际上是围绕通过地球中心的一个轴的旋转运动，通常用欧拉定理来表述V=ω×r。相对运动参数：，k、l板块上i、 j测站间基线长度的变化 率，kl ω是k对l的相对运动角速度，R i 、R j：i、j的坐标矢量，测两板块间若干ij B →由上式加权的最小二乘平差得kl ω。绝对运动参数：Ω(ω x ，ωy，ωz)角速度矢量，r(x,y,z)位置矢量，λ、φ经纬度。欧拉定理V=ω×r→地心系→地球近似为球体→站心系→角速度ω、纬度Φ、经度Λ。 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - = ? ? ? ? ? ? ? ? ? ? z y x G z y x x y x z y z V V V ω ω ω ，? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - = ? ? ? ? ? ? ? ? ? ? z y x G z y x r r r r r r V V V ω ω ω λ ? λ ? λ ? ? λ ? ? cos cos sin cos cos cos sin sin cos sin ， ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - = ? ? ? ? ? ? z y x r n e r r r r r V V ω ω ω ? λ ? λ λ ? λ cos cos sin sin sin sin cos ，。 区域地壳形变监测：区域空间尺度由几百到1000km的瞬变运动，GPS/INSAR/GRAVITY+VLBI/SLR；局部几百米到几十公里+LEVELING。 ?2000国家GPS大地控制网：2000国家重力、GPS大地网。缺陷：密度<全国天文大地网；所提供的低密度的3维地心坐标框架不能完整实现中国的3维地心系。解决：~GPS 网+全国天文大地网联合平差，将后者纳入3维地心系，并提高其精度和现势性，使我国的大地坐标框架在密度和分布方面实现我国3维地心大地坐标系的跨越式进步。CGCS2000通过~大地控制网的坐标和速度具体实现。参考历元2000.0。 ?GPS用于地壳⊥形变监测：①有可能利用GPS观测直接得到毫米级的大地高数据。重复GPS观测可以求定大地高的变化（或站心系U分量的变化）。②由于椭球体法线与该点夹角很小，大地高和正常高方向基本重合，∴可用大地高的变化代替正常高的变化，即可利用重复GPS观测取代精密水准以监测地面的升降变化。u分量与正常高变化：?h=?u/cos αAB。∵αAB=1°(S AB=110km)→cos αAB ≈1∴S AB ≈100km，可用?u代替?h。 ?InSAR：合成孔径雷达干涉测量一种以合成孔径雷达天线记录到的回波信号为信息源，利用干涉测量获取地球表面的三维地形、地表形变和地物特征变化等信息的测量技术。处理流程：配准、生成干涉图、平地相位消除、去地形影响、干涉图滤波、相位解缠。星载INSAR优点：①全球范围②成像测绘带宽③全天候④轨道稳定，多次成像重复性好⑤飞行平稳，无运动补偿问题⑥下视角小，图像阴影小。局限：①轨道固定，侦察

地壳形变课后题答案

第一章绪论 1. 地壳运动与变形的基本内涵是什么？ 地壳运动：是在地球内部构造应力的作用下，所引起的地壳一些元素的相对运动，它可以是垂直运动，水平运动或地倾泻运动。综合表现为大面积的地壳变形。 地壳形变：是指在地球内力和外力作用下，地球的地壳表面产生的升降、倾斜、错动等现象及其相应的变化量。 2. 你对地壳运动与变形的具体表现形式有哪些认识？ 其表现形式主要有：断裂，地震，火山，滑坡，泥石流，地面裂缝与沉降，冰川运动。 3. 现代大地测量学在地球动力学研究中可以发挥哪些作用？ 4. 简述地球物理大地测量学的发展。 第二章地壳形变测量 1. 地壳形变测量有哪几类？各有什么特点？ （1）全球板块运动监测 主要用来测定板块运动参数，测定大陆板块和海洋板块的内部形变，其观测手段主要采用VLBI、SLR 和GPS等空间测量技术。 （2）全国及区域地壳形变测量 测定亚板块及构造块体的地壳形变，给出全国大陆动力学的边界条件，以及全国大陆应力场、形变场变化过程的总体和分区特征；区域地壳形变测量主要测定块体边界与大地震有关的区域形变，它可以给出大陆内部地形变的时空演变图象。 （3）断层形变测量 在各活动构造块体边界上进行的近场构造变形测量。能够直接测定块体边界断裂及其不同段落的现今活动方式、相对位移速率以及它们随时间变化的过程，提供震间、震前、同震与震后滑动等构造活动的微动态信息。 目前以短水准、短基线、短边GPS网以及由水管倾斜仪、伸缩仪、蠕变仪、短边激光测距仪或重力仪组成的台阵等为主要手段。 （4）定点形变测量 主要包括地倾斜、地应变和重力(固体潮汐)台站观测。这种方法可以有效地监测地壳的连续变动，可以通过不同时间间隔的采样，在相当宽的频带范围内对地壳动力学现象进行观测。 2. GPS连续跟踪站测量属于定点形变测量吗?为什么? 属于 第三章：地球参考系统与地球参考框架 1. 什么是惯性参考系？ 其时间是均匀流逝的，空间是均匀和各向同性的；在这样的参考系内，描述运动的方程有着最简单的形式。这样的参考系就是惯性系 牛顿运动定律成立的参考系，称为惯性参考系，简称惯性系。 2. 我国的大地测量坐标基准有哪些？ 1954年北京坐标系；1980年西安坐标系（1980年国家大地坐标系）；2000国家大地坐标系。 3. 大地坐标系建立的基本原理是什么？ 4. 全球最优的协议地球参考架CTRF是如何建立的？ 由IERS组织，遵循以下原则①原点位于包括海洋和大气在内的整个地球的质心；②长度国际单位制为米，尺度为广义相对论下的局域地球框架；③定向为最初由国际时间局（BIH）所给出1984.0定向；④定向的时变通过一个关于全球的水平构造运动的非净旋转条件（地壳无整体旋转）来保证。

GNSS地壳形变异常检测理论与方法

一一 第４５卷一第６期测一绘一学一报 V o l ．４５,N o ．６ 一２０１６年６月A c t aG e o d a e t i c ae tC a r t o g r a p h i c aS i n i c a J u n e ,２０１６ 引文格式:X U K e k e ．T h e o r y a n d M e t h o d f o rD e t e c t i n g A b n o r m a l C r u s t a l D e f o r m a t i o nU s i n g G N S S [J ]．A c t aG e o d a e t i c a e t C a r t o g r a p h i c a S i n i c a ,２０１６,４５(６):７５６．(徐克科．G N S S 地壳形变异常检测理论与方法[J ]．测绘学报,２０１６,４５(６):７５６．)D O I :１０．１１９４７/j ．A G C S ．２０１６．２０１６０１５２ G N S S 地壳形变异常检测理论与方法 徐克科１, ２１．河南理工大学测绘与国土信息工程学院,河南焦作４５４０００;２．同济大学测绘与地理信息学院,上海２０００９２ T h e o r y a n dM e t h o d f o r D e t e c t i n g A b n o r m a l C r u s t a l D e f o r m a t i o nU s i n g G N S S X UK e k e １, ２１．S c h o o l o f S u r v e y i n g a n dL a n d I n f o r m a t i o nE n g i n e e r i n g o fH e n a nP o l y t e c h n i cU n i v e r s i t y ,J i a o z u o ４５４０００,C h i n a ;２．C o l l e g eo f S u r v e y i n g a n dG e o ＧI n f o r m a t i o no f T o n g j i U n i v e r s i t y ,S h a n g h a i ,２０００９２,C h i n a 一一地震的孕育过程是一个长期缓慢的过程, 地震发生时断层破裂所释放的能量只是一部分,其中有很大一部分能量在地震前后以无震蠕滑的形式释放.因此,检测这些形变异常信息对于地震危险性评估至关重要.随着G N S S 监测网络数据的持续积累,时空分辨率越来越高,就有可能从中挖掘出更有价值的形变信息.本文利用G N S S 时空数据, 从高精度数据处理二全球板块运动模型构建二多尺度速度场与应变场估计二地壳微动态形变异常检测和断层滑移时空反演等方面,研究了G N S S 地壳形变异常检测理论与方法.主要研究内容与结果如下: (１)重建了G N S S 测站速度估计模型:法方程重构模型二基线向量最小二乘模型二卡尔曼滤波模型和坐标时序模型.前３种模型将坐标二速度二年二半年周期项作为参数,在进行网平差的同时一并求解速度值.坐标时序模型引入白噪声和幂律噪声组合,分析了区域周期项振幅与相位空间分布的一致性.基于４种模型,估计了川滇陆态网G N S S 测站速度.结果差异在１mm /a 内,精度在亚毫米级.从而验证了这些速度估计模型的一致性. (２ )提出了基于统计假设检验和稳健估计的全球板块运动模型构建方法.利用I T R F ２００８V E L 求取的全球板块欧拉参数与先前模型具有较好的一致性.建立了相对欧亚板块背景场的中国大陆速度场,分析了中国大陆现今地壳运动特征. (３)构建了G N S S 多尺度速度场和应变场估计模型. 利用负位错模拟数据作了大量试验,结果表明,不同尺度的应变场具有检测不同空间范围的地壳形变的能力.利用２００９ ２０１１年陆态网数据估计并分析了中国大陆３ ８尺度的应变场及形变特征. (４ )联合卡尔曼滤波和主成分时空分析,集时空滤波与形变检测于一体,构建了瞬态无震蠕滑时空检测模型,进一步提高了数据时空信噪比.利用滇西南陆态网 G N S S 数据,检测其时空分布特征与２０１１年缅甸Mw ７．２级地震相对应.得出了滇西南区域的断层活动可能会受到缅甸地震带的影响. (５)提出了基于G N S S 基线面应变和G N S S 网形的 震前形变异常检测方法.得出２０１３年芦山M s ７．０地震和日本近年来的４次地震在震前数月的时间内均有不同程度的异常偏离.尤其G N S S 基线夹角二基线方位角和第一剪应变的异常变化更为突出.推断震前可能产生了强烈的左旋剪切构造应力变化,加速了芦山地震的孕育发生. (６ )以地壳形变检测与断层滑移反演为一体,构建了G N S S 主成分和卡尔曼滤波时空反演模型.通过模拟试验,得出了正确反演断层滑移时空分布所需要的最低信噪比和最优的台站分布密度.以２００５年苏门答腊Mw ８．７地震震后余滑和２００６年墨西哥慢滑移为例,反演了断层蠕滑时空分布. 中图分类号:P ２２８一一一一文献标识码:D 文章编号:１００１Ｇ１５９５(２０１６)０６Ｇ０７５６Ｇ０１基金项目:国家９７３项目(２０１３C B ７３３３０４ );国家自然科学基金(４１４０４０２３ )收稿日期:２０１６Ｇ０４Ｇ０５ 作者简介:徐克科(１９７９ ),男,副教授,２０１５年７月毕业于同济大学,获工学博士学位(指导教师:伍吉仓教授),研究方向为G N S S 数据处理与地壳形变分析. A u t h o r :X U K e k e (１９７９ ),m a l e ,r e c e i v e dh i sd o c t o r a l d e g r e e f r o mT o n g j iU n i v e r s i t y o nJ u l y ２０１５(P h Da d v i s o r :P r o f ．W U J i c a n g ),m a j o r si n G N S S d a t a p r o c e s s i n g a n d c r u s t a l d e f o r m a t i o na n a l y s i s E Ｇm a i l :１２x k k ＠t o n g j i ．e d u ．c n

GPS地壳形变

GPS地壳形变与地震预测 自20世纪80年代以来,VLBI、GPS、SLR等空间大地测量技术的精度得到大幅度提高,其观测结果已被应用于板块运动、冰后期反弹、火山、地震、地球自转和地球系统内部物质再分布等地球动力学过程的研究[1,2]。尤其是GPS观测技术,由于GPS接收机价格不是特别昂贵,观测比较容易实施,在地壳形变监测中得到了广泛应用。如为了监测地壳形变,日本在境内布设了1 000多个连续GPS观测站点,平均每隔30 km就有一个GPS观测站点[3];美国为了监测南加州圣安德烈斯断层带的运动和形变,布设了超过250个站的连续GPS监测网[4]。我国自20世纪80年代开始,在青藏高原、南北地震带、川滇和华北首都圈等地壳运动活跃区和地震监测重点区布设了一系列GPS观测网,用于研究这些区域的地壳形变和地震以及判断板块内部块体构造运动模型[5]。这些GPS区域网中的部分站点已经并入国家连续GPS观测网或定期复测GPS观测网,列入长期观测计划。随着连续和分期GPS网观测的持续,积累的地表运动观测数据越来越多。这是一类观测精度非常高(上百km范围内可以测量1mm/a的水平方向变化)、分布在监测区域地表上的站点坐标时间序列。利用这类时空观测数据提取精确可靠的地壳形变信息,反演地球动力学机制是现代大地形变测量数据处理的科学目标,为此,很多学者开展了一系列的理论和方法研究,也取得了不少应用成果。 从现有的研究成果来看,通过全球连续GPS网观测数据已经成功地解算了全球板块运动模型参数,其结果与基于地质和地球物理调查数据(百万年平均)得到的板块运动模型参数(NUVEL-1A)非常接近，连续或分期的GPS观测网成功地观测到了同震位移和震后位移[6],并且用于反演地震断层面上的位错分布，或结合地震波数据联合反演,增加人们对地震破裂过程的认识，通过GPS观测结果结合地质调查和地球物理勘探结果,研究划分板块内部次一级的构造运动块体,以及确定各块体之间的相对运动等[7]。目前,国内外研究关注的焦点和难点是如何通过GPS观测得到的站点坐标时间序列提取块体边界或断层带的地壳形变信息,进而研究如地震孕震形变等地球动力学过程[8]。

昆仑山大地震震后形变反映的地壳岩石流变特性

第30卷第9期 岩 土 力 学 V ol.30 No.9 2009年9月 Rock and Soil Mechanics Sept. 2009 收稿日期：2008-03-13 基金项目：中国科学院知识创新工程重要方向项目资助（No.KZCX2-VW-134）；国家自然科学重大研究项目重点基金（No.90814014）资助；973课题（2004CB418405）资助；地震基本科研业务项目（020*******）资助。 第一作者简介：张晁军，男，1965年生，博士，高级工程师，主要从事慢地震及计算地球动力学的研究。E-mail: zhangchaojun@https://www.doczj.com/doc/ff3857416.html, 文章编号：1000－7598 (2009) 09－2552－07 昆仑山大地震震后形变反映的地壳岩石流变特性 张晁军1， 2，石耀霖2，马 丽3 （1. 中国地震台网中心 信息网络部，北京 100045；2. 中国科学院研究生院地球科学学院 计算地球动力学实验室，北京 100049； 3. 中国地震局地震预测研究所，北京 100036） 摘 要：2001年昆仑山地震是我国近50年来最大的地震。用3种模型对震后中国地震局跨断层布设的4个GPS 站点记录到显著的震后形变进行了模拟研究。结果表明：单纯的上地壳10 km 为滞弹性的模型不能解释震后形变的幅度；30 km 的弹性上地壳覆盖在 40 km 的柔性下地壳上的松弛模型可以解释变形速率指数衰减的主要特征；而两个模型的结合不但能解释整 体指数衰减的特征，而且还能更好地拟合震后最初几周的较高形变速率。结果表明，在昆仑山断层两侧存在着流变性质的 差异。 关 键 词：震后形变；黏滞系数；流变模型；滞弹性；标准线性固体；Maxwell 体 中图分类号：P 554 文献标识码：A Numerical simulation of crust rheological property reflected by post-seiemic deformations of Kunlun large earthquake ZHANG Chao-jun 1, 2 , SHI Yao-lin 2, MA Li 3 (1. China Earthquake Network Center, Beijing 100045, China; 2. Laboratory of Computational Geodynamics, Graduate School, Chinese Academy of Sciences, Beijing 100049, China; 3. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China) Abstract: The Mw7.8 Kunlun earthquake of 14 November, 2001, in the northern Tibetan Plateau of China, is the largest event in the Chinese continent in the latest 50 years. In this paper, layered viscoelastic models are calculated by PSGRN/PSCMP code; and the results are used to fit post-seismic deformation obtained from four GPS stations of CEA across the earthquake fault. The results show that model of a single anelastic layer of 10 km thick upper crust overlying an elastic lower crust cannot explain the amplitude of deformation; the relaxation model of 30 km elastic upper crust overriding 40 km ductile lower crust can explain the main features of exponential attenuation of post-seismic deformation. Combination of the two models, however, can fit the character better with deformation rate being higher in the first few weeks and slower rate thereafter. The viscous layer of lower crust provide the main control of the post-seismic deformation of long term decay of deformation in months, while the upper anelastic layer may contribute to the high rate of initial few weeks after the main earthquake. The result also suggests that rheological differences exist at two sides of the Kunlun Mountains fault. Key words: post-seismic deformation; viscosity rheological model; anelasticity standard line solid; maxwell body 1 前 言 2001年11月14日17时26分，在昆仑山口地区发生了7.8级强烈地震（美国国家地震信息中心NEIC 地震目录）。震后中国地震局利用GPS 开展了相对密集的地壳形变观测。穿越震区沿青藏公路跨断层布设了4个固定GPS 站点：格尔木（JB30）、昆仑山（KLGD ）、不冻泉（BDGD ）和五道梁 （WDGD ），其中不冻泉和五道梁位于东昆仑断裂南侧，距地震地表破裂带距离分别为30 km 和 130 km ；昆仑山和格尔木位于东昆仑断裂北侧，分别距地震地表破裂28 km 和100 km （图1）。从 JB30、KLGD 、BDGD 和WDGD 站点震后地表蠕动位移东向（E ）运动的时间序列[1]可以看出，在断层南侧虽然WDGD 和BDGD 两者与断层距离不同，东向运动幅度与时间变化特征仍很相似，变形速率

基于GPS数据的几种地壳形变分析方法

基于GPS数据的几种地壳形变分析方法 张永奇 长安大学地质工程和测绘学院，西安（710054） 摘要：主要介绍了几种在地壳形变分析中的方法，介绍了他们的原理，以及在地壳形变监测分析中的应用，同时也说明了几种方法的优缺点。 关键字：DDA；位错模型；刚体旋转加均匀应变模型；经验正交函数分解方法模型；数值流形法 引言 地震的孕育和发生在本质上是地壳内部的应变能逐渐积累并突然释放的结果。伴随着大震孕育或应变能的显著积累, 岩石圈表层必然会表现出某种形式和量级的地壳形变。基于这样的认识, 地壳形变监测一向是地震监测或地震危险性分析的最重要手段之一[1]。 自20 世纪80 年代以来,VLBI、GPS、SL R 等空间大地测量技术的精度得到大幅度提高,其观测结果已被应用于板块运动、冰后期反弹、火山、地震、地球自转和地球系统内部物质再分布等地球动力学过程的研究。尤其是GPS 观测技术,由于GPS 接收机价格不是特别昂贵,观测比较容易实施,在地壳形变监测中得到了广泛应用。我国自20 世纪80 年代开始,在青藏高原、南北地震带、川滇和华北首都圈等地壳运动活跃区和地震监测重点区布设了一系列GPS 观测网,用于研究这些区域的地壳形变和地震以及判断板块内部块体构造运动模型。这些GPS 区域网中的部分站点已经并入国家连续GPS 观测网或定期复测GPS 观测网,列入长期观测计划。随着连续和分期GPS 网观测的持续,积累的地表运动观测数据越来越多。这是一类观测精度非常高(上百km 范围内可以测量1mm/ a 的水平方向变化) 、分布在监测区域地表上的站点坐标时间序列。利用这类时空观测数据提取精确可靠的地壳形变信息,反演地球动力学机制是现代大地形变测量数据处理的科学目标,为此,很多学者开展了一系列的理论和方法研究,也取得了不少应用成果[3]。 进一步的基于GPS观测数据的地壳运动和形变的分析方法，尚处于研究阶段。目前初步的应用方法大体可以分为两类：1.纯数学方法，包括多面函数拟合法，多项式拟合法，双三次样条函数法，协方差经验函数法等；2.以物理力学为基础的方法，如欧拉矢量法，非连续形变分析法（discontinuous deformation analysis ，DDA）等。前者尽管在数学上是自洽的，但物理力学含有不清，且都是以地壳连续性假设为基础的；后者物理力学意义清楚，且考虑了地壳的非连续性，但描述块体形变不够精细和有效。如DDA对块体形变的描述仅用3个参数，过于简单只适合与小尺度块体系统的变形分析。 本文主要介绍现今在地壳形变分析的一些常用方法，如非连续变形分析方法，位错运动模型，刚体旋转加均匀应变模型，经验正交函数分解方法模型，数值流形法。以及对这些方法的优劣进行比较，在此基础上对一些可能应用在地壳形变分析的新方法进行展望，如小波理论，灰色理论，神经网络等。 1.非连续变形分析(DDA)方法 模拟介质不连续缝的历史可追溯到30年前的Goodman、Taylor和Brekke等教授发展的节理单元。对岩土裂缝的数值计算发展很快,并已在岩石工程中得到广泛应用。Cundall介绍的离散元法现在被广泛应用于节理或块状岩石。两者是用虚拟力来调整滑动和阻止块体重叠的一种方法,有时候可达到稳定。1993年,由石根华博士发展的非连续变形分析是全新的数值