（大地测量学与测量工程专业论文）有限元及数值方法在大地测量反演中的应用.pdf

摘要 摘要 在最近的几十年中，有限单元法发展迅速，不仅应用于动态分析，流体动 力学以及航空航天，机械，土木工程等领域，在地球物理学领域也应用广泛。 有限元法不仅能够解决各种复杂问题，而且还具有相当高的精度。 本文主要介绍有限元反演方法及其在大地测量中的应用，研究了分离节点 技术表述断层及其在大地测量反演中的应用，对二维平面问题进行算例分析发 现，在引入断层位移情况下，断层节点及其周围的应力场和应变场发生明显改 变，包括主应力和主应变大小及主方向的改变。分析还表明，在弹性体区域中， 远离断层周围的区域主应力分布影响不大，这与圣维南原理也是相符的。 利用矩阵奇异值分解方法研究了断层位错模型参数反演观测方程系数矩阵 的奇异值分布规律。探讨了位移观测网相对于断层的位置对反演断层位错模型 参数结果的影响。模拟算例表明，位移观测网应布设在以断层为中心的周围， 以提高反演解的稳定性。 建立了南加州地区的有限元模型，运用有限元法进行了板块边界力的反演 的模拟实验和实际计算。通过数值模拟，探讨和分析了几种因素对反演结果的 影响，结果表明有限元法进行大地测量反演有较好的效果。 利用南加州地震中心的速度场数据，在将北美板块固定时，反演其比较近 似的边界力，在断层位置作用2 m m 向西方向水平位移和l m m 向北方向水平位 移，此时模拟结果较好。并以此研究其主应变场和主应力场得出，圣安德列斯 断层周围具有一定的张拉应力和应变，断层呈现出右旋运动的特征，这与其它文 献资料相吻合。 关键词：反演，有限元法，分离结点法，数值方法，断层 a b s t r a c t a b s t r a c t i nr e s e n ts e v e r a ld e c a d e s ，t h ef i n i t ee l e m e n tm e t h o d ( f e m ) h a sd e v e l o p e d p r o m p t l y ，i ta p p l i e sn o to n l yi nt h ea r e ao fd y n a m i ca n a l y s i s ，f l u e n td y n a m i c s ， a e r o s p a c ee n g i n e e r i n g , m a c h i n e ，c i v i le n g i n e e r i n ge t c ，a l s oi nt h e a r e al i k e g e o p h y s i c s u s i n gt h i sm e t h o d ，al o to fd i f f i c u l tp r o b l e m sa r es o l v e dw i t hh i g h p r e c i s i o n i nt h i st h e s i s ，t h ei n v e r s i o no ff e mw h i c ha p p l i e si ng e o d e t i ci n v e r s i o ni s i n t r o d u c e d , a n dt h es p l i tn o d et e c h n i q u ei ss t u d i e dw h e naf a u l ti si n t r o d u c e d a t w o - d i m e n s i o n a lp r o b l e mi sa n a l y z e da n df r o mt h i sac o n c l u s i o ni sd r o w nt h a t ，w h e n i n t r o d u c i n gt h ef a u l td i s p l a c e m e n t s ，t h es t r e s sa n ds t r a i nf i e l d sc h a n g e so b v i o u s l yo n o ra r o u n dt h ef a u l t - n o d e s ，i n c l u d i n gt h es i z ea n dd i r e c t i o n a l s o ，t h ee f f e c t so ft h e p r i n c i p l es t r e s sa n ds t r a i nw h i c ha r ea w a yf r o mt h ef a u l ta r e aa r en o ts om a g n i f i c e n t , t h i si si nl i n ew i t l lt h es a i n t - v e n a n tp r i n c i p l e m a t r i xs v d ( s i n g u l a rv a l u ed e c o m p o s i t i o n ) i su s e df o r s t u d y i n g t h e d i s t r i b u t i o no fs i n g u l a rv a l u e so ft h eo b s e r v a t i o ne q u a t i o nc o e f f i c i e n tm a t r i xi n i n v e r s i n go ff a u l td i s l o c a t i o np a r a m e t e r sb yd i s p l a c e m e n tm e a s u r e d t h ee f f e c to f r e l a t i v ep o s i t i o nb e t w e e ns u r v e y i n gn e t w o r ka n dt h ef a u l tt ot h ei n v e r s i o nr e s u l t sh a s b e e ns t u d i e d t h er e s u l t so fs i m u l a t i v ee x a m p l e ss h o wt h a tt h es u r v e y i n gn e t w o r k s h o u l db ed i s t r i b u t e da r o u n dt h ec e n t e ro ft h ef a u l tm o n i t o r e ds oa s t oi m p r o v et h e s t a b i l i t yo fi n v e r s i o n t h ef e mm o d e lo fs o u t h - c a l i f o r n i ah a ss e tu p t h ei n v e r s i o nc o m p u t a t i o nf o r p l a t ed r i v i n gf o r c e sa c t i n ga tt h eb o u n d a r i e si sm a d eu s i n gf i n i t ee l e m e n tm e t h o db y s i m u l a t i o ne x a m p l ea n da c t u a le x a m p l e t h es i m u l a t i o ne x p e r i m e n ti l l u s t r a t e st h e f e a s i b i l i t ya n dv a l i d i t yo f t l 蚯sm e t h o d 。 t h ed a t ao fv e l o c i t yf i e l do fs o u t h e r nc a l i f o r n i ae a r t h q u a k ec e n t e ri su s e dt o i n v e r s et h ea p p r o x i m a t eb o u n d a r yf o r c e w h e nt h en o r t h e r na m e r i c a np l a t ei sf i x e d t h ee x p e r i m e n t sd i s p l a yt h a tt h es i m u l a t i o nr e s u l tw i l lb eb e t t e r ，w h e nah o r i z o n t a l d i s p l a c e m e n to f2m i l l i m e t e r sa n da v e r t i c a ld i s p l a c e m e n to f1m i l l i m e t e ra r ea p p l i e d n a b s l r a c t a c c o r d i n gt ot h ep r e v i o u sr e s u l tt h ep r i n c i p l es t r e s sa n ds 仃a i nm o u n ds a na n d r e a s f a u l ta r es t u d i e df u r t h e rt of m dt h a tt h e r ea r et e n s i l es t r e s sa n ds 仃a i na r o u n dt h ef a u l t w h i c hs h o w sar i g h tl a t e r a ls u i k e - s l i pm o t i o nc h a r a c t e la n dt h i si sc h e c k e db ys o m e d o c u m e n t s k e yw o r d s ：i n v e r s i o n ，f e m ，s p l i tn o d em e t h o d ，n u m e r i c a lm e t h o d ，f a u l t 学位论文版权使用授权书 本人完全了解同济大学关于收集、保存、使用学位论文的规 定，同意如下各项内容：按照学校要求提交学位论文的印刷本和 电子版本；学校有权保存学位论文的印刷本和电子版，并采用影 印、缩印、扫描、数字化或其它手段保存论文；学校有权提供目 录检索以及提供本学位论文全文或者部分的阅览服务；学校有权 按有关规定向国家有关部门或者机构送交论文的复印件和电子 版；在不以赢利为目的的前提下，学校可以适当复制论文的部分 或全部内容用于学术活动。 学位论文作者签名： 年月 日 同济大学学位论文原创性声明 本人郑重声明：所里交的学位论文，是本人在导师指导下，进行 研究工作所取得的成果。除文中已经注明引用的内容外，本学位论文 的研究成果不包含任何他人创作的、已公开发表或者没有公开发表的 作品的内容。对本论文所涉及的研究工作做出贡献的其他个人和集 体，均已在文中以明确方式标明。本学位论文原创性声明的法律责任 由本人承担。 学位论文作者签名： 年月 日 第l 章引言 第1 章引言 任何事物都是具有联系的。人与人之间，人与自然之间，人与社会之间， 无不包含千丝万缕的联系。从整体到局部，从简单到复杂，又从复杂到简单。 仔细一想，有限元方法与之是如此相似，将一个连续复杂的物体离散为一个个 简单明了的单元来研究，然后通过单元之间的“局部 关系，又将其一个个单 元组装成为“整体”，从而达到其近似计算的目的。 1 1 有限元法进展和展望 1 1 1 有限单元法进晨 有限元法的基本思想可以追溯到2 0 世纪5 0 年代，1 9 4 3 年，c o u r a n t 把整 个求解区域分成很多的三角形子区域，然后假定待求函数为线性函数，把三角 形顶点处的函数值作为未知数，并放宽边界条件，只要求在有限节点上满足就 可以【1 1 。1 9 4 5 至1 9 5 5 年，a r g y r i s 等人在结构矩阵分析方面取得重大进展【2 】【3 1 。 1 9 5 6 年，t u r n e r ，c l o u g h 等人把刚架位移法的思路，推广应用于平面弹性力学 问题：他们把连续体划为三角形和矩形单元，单元中的位移函数采用近似表达 式，推导单元的刚度矩阵，建立结点位移与结点力之间的单元刚度方程【4 】【5 1 。 1 9 6 0 年，c l o u g h 进一步求解了平面弹性问题，并第一次提出了“有限单元法 的名称【6 】。初期的有限元法是建立在虚功原理的基础上。1 9 6 3 1 9 6 4 年m e l o s h 和j o n e s 等人证明了有限元法是基于变分原理的里兹法的另一种形式，确认了 有限元法是处理连续介质问题的一种普遍方法，扩大了有限元法的应用范围 【7 】【1 0 】 o 从2 0 世纪6 0 年代后期开始，进一步利用加权余量法，主要是伽辽金法【n 】， 来确定单元特性和建立有限元求解方程，使之应用于已知问题的微分方程和边 界条件、但变分的泛函尚未找到或者根本不存在的情况，进一步扩大了有限元 法的应用领域。 1 9 6 7 年首次出版有限元专著 结构与连续力学的有限元法，由 第l 章引言 z i e n k i e w i c z 与y k c h e u n g ( 张佑启) 合作。z i e n k i e w i c z 等把复杂结构的计算 问题转化为简单单元的分析和集合问题，许多经典的数学近似方法以及工程中 所用的各直接近似方法都属于这一范畴【1 2 】【。有限单元这一术语的出现，意味 着直接应用可用于离散系统的标准研究方法。在概念上，使我们对方法的理解 得到改善；在计算上，可对各种问题应用统一的方法，并研制出标准的计算程 序。 近3 0 多年来，伴随着电子计算机科学和技术的快速发展，有限元法发展迅 速兴起。从静态弹性力学分析开始，对不均匀介质和复杂形状物体进行应力和 变形分析并得到很精确的结果，因此很快发展到应用于连续介质物体的其他分 支：动态分析，应力波，热传导渗流，流体动力学，电场磁场以及航空航天， 机械，土木工程等领域【1 5 】【1 8 】。在地球物理学中，由于地质构造形态的复杂性， 地球介质的不均匀性，在有限单元法出现前，只能进行一些简单的情况分析。 6 0 年代后期，有限元法被引进地质构造分析以来，得到迅速的发展，对地球科 学做定量分析做出重要贡献。已出现多种新单元( 先后有等参元、高次元、不 协调元、拟协调元、杂交元、样条元、边界元、罚单元，还有半解析的有限条 等不同单元) 和求解方法( 如半带宽与变带宽消去法、超矩阵法、波前法、子 结构法、子空间迭代法等) 。能解决各种复杂耦合问题的软件和软件系统不断涌 现。对网格自动剖分和网格自适应过程的研究，大大加强了有限元法的解题能 力，使有限单元法逐渐趋于成熟。有限元法作为一种离散化的数值解法，也已 成为应用数学的一个新的分支【1 9 】 f 髓方法即有限单元法，作为一种技术更多的与f 髓软件的发展紧密的结 合起来。某种主软件的f 翻方法必然会一直朝该f 髓方法的方向发展，只有当 新的f 跚方法比现有的f e m 方法更加优越时才会放弃现有的f 尉方法，从而使 f 脚方法有较大的发展。因此f e m 方法仍将在当今世界广泛应用。 当今主流的f 脒软件有德国的a s k a 、英国的p a f e c 、法国的s y s t u s 、美国的 a b q u s 、a d i n a 、a n s y s 、b e r s a f e 、b o s o r 、c o s m o s 、e l a s 、m a r c 和s t a r d y n e 等公 司的产品。这些软件所代表的方法有tc o s m o s 软件使用的快速有限元算法( f f e ) 。 在传统有限元分析的数值计算方法之中，有直接计算法( d i r e c ts o l v e r ) 与迭 代法( i t e r a t i v e ) 两种由于在过去的经验中。迭代法一直无法直接而有效 的保证数值计算的收敛性，快速有限元法是一种可以保证收敛性的迭代法该 方法计算速度也很快。m a r c 软件以l a g r a n g e 算法为主，兼有a l e 和e u l e r 算法； 2 第1 章引育 以显式求解为主，兼有隐式求解功能。a n s y s 软件有直接求解器，如波前求解器， 可计算出线性联立方程组的精确解。a n s y s 程序还提供了一个有效的稀疏矩阵求 解器，它既可用于线性分析，也可用于非线性分析。即要求求解精度又要求求 解时间的静态及瞬态分析中，该求解器可代替迭代求解器。稀疏矩阵求解器只 能用于真正的对称矩阵，与波前及其它直接求解器相比，稀疏矩阵求解器能显 著加速求解速度闭j 。 i i 2 有限单元法发展趋势 随着有限元技术的不断发展和成熟，有限元法越来越多的应用于实际工程 和技术领域，当今国际上有限元法的发展趋势： 1 从单纯的结构力学计算发展到求解许多物理场问题 有限元分析方法最早是从结构化矩阵分析发展而来，逐步推广到板、壳和实 体等连续体固体力学分析，实践证明这是一种非常有效的数值分析方法。而且从 理论上也已经证明，只要用于离散求解对象的单元足够小，所得到的解就可足够 逼近于精确值。所以近年来有限元方法已发展到流体力学、温度场、电传导、 磁场、渗流和声场，构造应力场等问题的求解计算，最近又发展到求解几个交 叉学科的问题。 2 由求解线性工程问题进展到分析非线性问题 随着科学技术的发展，线性理论已经远远不能满足设计的要求。例如建筑行 业中的高层建筑和大跨度悬索桥的出现，就要求考虑结构的大位移和大应变等 几何非线性问题：航天和动力工程的高温部件存在热变形和热应力，也要考虑材 料的非线性问题：诸如塑料、橡胶和复合材料等各种新材料的出现，仅靠线性计 算理论就不足以解决遇到的问题，只有采用非线性有限元算法才能解决。非线性 的数值计算是很复杂的，它涉及到很多专门的数学问题和运算技巧，很难为一般 工程技术人员所掌握。为此近年来国外一些公司花费了大量的人力和投资开发 诸如m a r c 、a b q u s 和a d i n a 等专长于求解非线性问题的有限元分析软件，并广 泛应用于工程实践。这些软件的共同特点是具有高效的非线性求解器以及丰富 和实用的非线性材料库【2 l 】。 3 更为强大的网格处理能力有限元法求解问题的基本过程主要包括： 分析对象的离散化、有限元求解、计算结果的后处理三部分。由于结构离 第l 章引言 散后的网格质量直接影响到求解时间及求解结果的正确性与否，近年来各软件 开发商都加大了其在网格处理方面的投入，使网格生成的质量和效率都有了很 大的提高，但在有些方面却一直没有得到改进，如对三维实体模型进行自动六 面体网格划分和根据求解结果对模型进行自适应网格划分，除了个别商业软件 做得较好外，大多数分析软件仍然没有此功能。自动六面体网格划分是指对三 维实体模型程序能自动的划分出六面体网格单元，现在大多数软件都能采用映 射、拖拉、扫略等功能生成六面体单元，但这些功能都只能对简单规则模型适 用，对于复杂的三维模型则只能采用自动四面体网格划分技术生成四面体单元。 对于四面体单元，如果不使用中间节点，在很多问题中将会产生不正确的结果， 如果使用中间节点将会引起求解时间、收敛速度等方面的一系列问题，因此人 们迫切的希望自动六面体网格功能的出现。自适应性网格划分是指在现有网格 基础上，根据有限元计算结果估计计算误差、重新划分网格和再计算的一个循 环过程。对于许多工程实际问题，在整个求解过程中，模型的某些区域将会产 生很大的应变，引起单元畸变，从而导致求解不能进行下去或求解结果不正确， 因此必须进行网格自动重划分。自适应网格往往是许多工程问题如裂纹扩展、 薄板成形等大应变分析的必要条件。 1 2 大地测反演理论进展和展望 1 2 1 反演理论进展 正演问题是按事物的一般原理( 或模型) 以及相关的条件( 初始条件、边界条 件) 来预测事物的结果。与正演问题相反，反演问题则是由结果及某些一般原理 ( 或模型) 出发去确定表征问题特征的参数( 或称模型参数) 2 2 】大地测量反演理 论是利用大地测量变形观测数据，研究地球表面客观形变的演化特征和规律， 推求地球内部的物性参数和特征，揭示地球内部动力过程的- - n 边缘学科是 大地测量学的重要组成部分1 2 3 。虽然大地测量反演是一个较为年轻的学科分支， 但最近几十年国内外在该领域的研究还是取得了令人欣慰的成绩。1 9 6 7 1 9 7 0 年间，美国地球物理学家b a c k u s 和应用数学家g i l b e r t 发表了一系列反演方面的 文章，这些文章奠定了反演理论的基础，不仅为以后的反演研究提供了重要的 参考，而且引入了很多重要概念( 如分辨率矩阵，信息矩阵，协方差矩阵等) ， 第1 章引言 这对理解反演问题的实质和评价反演效果有重要意义。这一理论被称为 b a c k u s - g i l b e r t ( 或b g ) 理论，它的反演模型为连续模型彤】。1 9 8 7 年，t a t a n t o l a 从概率的角度提出了广义反演方法，该方法的优点是把模型协方差和数据协方 差引入到反演模型中，从而有效地减小了解的非唯一性，方便了反演工作的实 施。这些理论为反演工作奠定了基础，形成较为完备的理论系统，成为反演理 论基础【2 5 1 。 1 模型构制：大地测量反演问题的模型构制是一项非常重要的工作，它直 接影响其它工作的作用和效率。不合理的模型构制可导致大地测量反演问题出 现错误的解甚至无解。根据构制模型的参数数量和观测数据数量之间的关系， 大地测量反演问题的求解可分为适定问题，超定问题，欠定问题和混定问题。 2 解的适定性问题：这个问题有三个内容，即解的存在性、唯一性和稳定 性。解的存在性是解决对于给定的一组大地测量数据是否一定存在能拟合这组 数据的解或者模型。唯一性是表示能拟合这组数据的模型是否唯一。解的稳定 性则是表示反演得到的模型参数对这组数据的敏感程度，如果观测数据的稍微 扰动所引起模型参数的扰动量也是微小的，则称该模型是稳定的。从重要性来 讲，对于一个问题的提出，首先应该满足解的存在性条件，其次是唯一性条件 和稳定性条件。如果不满足存在性条件，证明该模型是错误的，应更换模型。 但在大地测量反演的实践中，解的存在性已被大量的事实所证型2 3 1 。问题的提 出同样应该满足唯一性条件和稳定性条件，如果一个问题的解不唯一证明我们 对该问题的认识程度还不够，存在偏差，而当问题的解表现不稳定时，降低了 该解的可信度和实用性。 3 求解方法：大地测量反演问题的求解方法是人们最为关心、工作量最大 的二项内容f 2 3 】。一个有效的反演方法不仅能得到合理的模型解释，而且可以节 省大量的人力、物力和财力，因此广大地学工作者才坚持不懈地探索、研究新 的更加有效的反演方法，特别是非线性反演方法。 4 解的评价：无论线性反演还是非线性反演，大地测量反演问题的解都是 非唯一的，怎样从这些解中提取有用的信息就是解的评价所要解决的问题。由 于大地测量所研究的问题较复杂，对于反演解的评价也相对困难。目前常使用 的方法有比较先验信息与反演值之差法、求观测值与反演模型理论值拟合差法 和多种算法检核法。构制好模型后，在处理观测数据和模型参数时存在三种不 同的观点和方法，分别为【刈： 第l 章引言 1 ) 把观测数据和模型参数( 既可以是离散模型，也可以是能用有限个参数 表征的连续模型) 都看成是随机变量，通过研究它们所遵循的概率分布对观测资 料进行反演。 2 ) 把观测数据视为随机变量，把模型看成是由一些确定参数所决定的。 反演的任务就是估算这些待定的参数以及它们可能具有的误差。 3 ) 视观测数据为随机变量模型为连续函数，形成一套连续介质的反演 理论和方法，这就是著名的b u c k u s - g ii b e r t 反演法。 随着科学技术的迅速发展，最近十几年各种观测资料在时空分布范围和测 量精度上都有很大的改进和提高。这为大地测量反演的研究工作提供了必需的 第一手资料，也为大地测量反演的研究方向拓宽了道路。特别是计算机科学的 飞快发展，为我们可能进行的大量的计算和设计工作提供了保障。在新形式下， 大地测量反演的研究工作应该充分利用现有的有利条件，发挥各种优势，在提 高反演精度的同时尽可能拓宽反演的应用面。 1 2 2 反演方法进展 一个有效的反演方法不仅能节省大量的人力、物力和财力，而且可以得到 合理的模型解释，所以大地测量反演方法的选取是十分重要的。大地测量反演 方法按分类的依据不一样，可以有多种分类方式，下面简单介绍几种常用的分 类方式：按照大地测量反演所使用的观测量种类可分为单一数据反演和多种数 据联合反演：按照大地测量反演过程不同可分为逆解法反演和直接法反演：按照 大地测量反演所使用的优化算法不同可分为局部优化法反演和全局优化法反演： 按照大地测量反演使用的计算方法不同可以分为解析法和数值法。传统的大地 测量反演方法存在很多弊端。由于大地测量反演所研究的问题比较复杂，一般 都是非线性问题，所以怎样选择一种合理有效的、能够反演非线性问题的大地 测量反演方法是至关重要的传统的处理方法是用参数置换法或者泰勒级数展 开法把非线性问题转化为线性问题不管采用哪种线性化方法，都必然或多或 少地存在线性化误差，这就影响了反演模型的精度。对于线性反演问题，当反 演的模型参数较多时，传统方法除计算量大的缺点外也可能由于参数的相关 性而出现病态的情况这样会使得反演的解不稳定甚至出现较大的偏差。而 且传统的反演方法一般为局部最优的优化法，初值的选择对这种方法的反演效 6 第l 章引言 果十分重要，如果选择的初值与真实的值相差比较大，很容易陷入局部最优。 而实际的操作中选择这样合适的初值是很困难的在大地测量反演方法中，对 于不需要计算导数的方法其效率一般较低，收敛速度也较慢，而对于需要计算 导数的方法，虽然其效率一般较高，收敛速度也较快，但如果传统方法的目标 函数不存在导数或者导数难以计算时，这些传统方法是不便使用的【2 7 】。如上所 述，寻求优秀的大地测量反演问题的求解方法是很多国内外学者的追求目标。 虽然大地测量反演的理论和方法远没有地球物理反演的理论和方法那样成熟丰 富，但国内外大地测量工作者还是克服重重困难，取得了令人骄傲的成绩。利 用大地测量数据反演地震的位错模式【2 8 】，是前期大地测量反演研究成果的杰出 代表。张祖胜把这一反演方法进行系统的归纳，并用于地震断层【2 9 】。 1 3 有限单元法在大地测量反演理论中的进展 有限单元法在大地测量反演中的应用理论已取得丰硕的成果：赵少荣运用 有限元法和正则化方法建立了动态大地测量数据反演的基本理论，并在此基础 上进一步建立了基于力学模式的动态大地测量反演的一般理论【3 们。许才军从固 体力学的基本方程出发，结合数值方法，提出用大地测量数据反演构造应力场 理论，研究了由大地测量数据反演线弹性构造应力场的方法1 3 1 】。党亚民把分割 节点技术应用到有限元数值解算中，解决了有限元解算中节点位移不连续的问 题【3 2 l ；安美建把具有良好的全局搜索能力的随机优化算法即遗传算法和有限单 元法结合应用于大地测量反演中，并取得较好的效果【3 引。许忠淮等提出了两种 在作二维线弹性问题有限单元分析时，根据区域内部应力方向观测结果反演区 域边界作用力的方法，即约束反演法和应力张量拟合法p 4 】。有限单元法在实际 的大地测量反演中也有很好的应用。许才军根据弹性理论的叠加原理，研究了有 限元分析中的单元划分问题。通过有限元结点位移反演边界力系数与理论边界 力系数进行比较来判别断裂区单元划分是否合理，应用模拟数据定量分析了二 维和三维有限元法的单元划分，为有限元法在地学研究中的实际应用提供参考 依据【3 1 1 。安美建利用遗传有限单元反演法对东亚部分地区构造应力场进行二维 平面应力模型拟合反演运算，得到了影响该地区现今构造应力场分布的边界作 用力相对值和方向p 引。许忠淮根据中国东部及附近地区的主压应力轴方向观测 结果，采用二维线弹性有限单元模型，反演出了周围板块边界作用力的相对大 7 第1 章引言 小【列。焦明若等利用粘弹性有限单元模型，采用边界位移速率结果，对中国大 陆及其邻区的基本构造应力场进行了数值模拟，解释了中国大陆及其邻区地震 分布的特点口5 1 。黄立人等根据研究得到的华北部分地区地壳构造分布的空间特 征和在深度方向上的地壳介质分层特征，设计了三维线弹性有限单元介质模型， 并在顾及和不顾及区域本底应力场作用方式的情况下，计算模拟了构造块体间 的相对运动【3 6 1 ，王凯英等利用有限元模拟川滇地区断层相互作用，分析显示断 层的相互作用是块体非均匀运动过程中应力场调整的反应，是块体运动的结果 1 3 7 j 。为将数值流形方法从平面扩展到球面，用球面三角形网格作为数学网格去 覆盖研究区域，得到球面上的数值流形单元，推导了球面上数值流形方法的有限 单元覆盖的应变矩阵及刚度矩阵、初应力矩阵、点荷载矩阵、体荷载矩阵、惯 性力矩阵和固定点矩阵，从而为将数值流形方法应用于大范围地壳运动的数值 模拟打下了基础【3 引。蒋锋云利用青藏块体区域内部相对运动的水平速度场，通过 构建该区二维有限元模型，结合最优化原理，反演区域分段边界位移。将反演结 果作为边界约束并将变换参考框架后的水平速度场作为有限元节点载荷，模拟 区域内部水平运动场并计算其主应变率场、最大剪应变率场，以分析这不同时段 中两种场的变化特征及其与地震孕育之间的关系【3 们。根据要解决的问题的要求， 采用适当的数值处理方法，一般都可以满足精度的要求。但由于大地测量反演 问题的复杂性，求其精确解的可能性不大。而有限元法的精度可以满足大地测 量反演问题的要求，而且可以方便地编程实现具体问题的计算，所以近些年有 限单元法在大地测量反演中得以广泛采用。作为一个较为年轻的学科分支，大 地测量反演理论正以良好势头不断发展和完善，这个过程中有限单元法起着非 常重要的作用。 1 4 本文研究的主要内容 本文首先介绍有限元法基本原理分割结点技术原理及其应用，以及反演 常用的几种数值方法，通过模拟算例分析，编制适当的程序进行位移场，主应 变场，主应力场的分析，对比其结果的可靠性，并利用已观测的g p s 数据用 有限元软件a n s y s 建立相应的模型，通过适当的网格划分和边界条件的选取， 对美国南加州的圣安德列斯断裂带进行模拟算例分析考虑太平洋板块对加州 海岸的挤压作用，进行数值模拟其位移场及应力应变场并与观测值进行对比， 第1 章引言 得到反演结果并对周围构造应力场进行相应的解释，以及对可能发生的地震灾 害进行适当的预测。 研究内容主要有： 1 本文第一章为引言，介绍有限元法及大地测量反演理论进展，有限元法 在大地测量反演中的应用的文献综述。 2 接下来第二章讲述有限单元法的基本原理及在大地测量反演中的实现。 3 第三章主要利用分离结点技术原理及在大地测量中的应用。 4 第四章研究了基于数值分析的断层位错模型的参数反演理论及算例分 析，得出相应的结论。 5 第五章利用南加州地区的g p s 位移场数据与用有限元软件加渺s 对其 断层带划分后的位移场进行反演分析，得出相应的结论。 6 第六章为结论分析和展望，由前面的公式推导及算例分析最后得出结论 并提出展望。 9 第2 章有限元基本原理及其在大地测量反演中的应用 第2 章有限元基本原理及其在大地测量反演中的实现 2 1 有限单元法原理 有限单元法的基本原理是将连续的求解区域离散为一组有限个、且按一定 方式相互联结在一起的单元的组合体。由于单元能按不同的联结方式进行组合， 且单元本身又可以有不同的形状，因此可以模型化几何形状复杂的求解域【4 0 1 。 它广泛的应用于各种类型的偏微分方程问题的求解，并通过变分，使原偏微分 方程变为相应的泛函，这样问题就变成求相应的泛函的驻值。用有限单元法求 解问题时是利用在每一单元内假设的近似函数来分片地表示全求解域上待求的 未知场函数。单元内的近似函数通常由未知场函数或其导数在单元的各个结点 的数值和其插值函数来表达。这样一来，一个问题的有限元分析中，未知场函 数或其导数在各个结点上的数值就成为新的未知量( 也即自由度) ，从而使一个 连续的无限自由度问题变成离散的有限自由度问题1 4 0 i 。使用有限单元法，从数 学角度考虑包括以下几点【4 1 1 ： 1 从古典变分方法原理去定义微分方程边值问题的广义解以及在古典变分 方法的框架对有限元进行理论分析。 2 保证偏微分方程边值问题的提法正确，即要求解存在、唯一和稳定，即 保证数值解法是可靠的。 3 有限元中重要的一点是采用了分块多项式插值函数，因此，有限元的误 差估计转化为插值逼近的误差估计问题。 4 有限元的收敛性和误差估计。 2 2 有限单元法实现过程 由于本文将用到四节点等参单元，特别是在下章分割节点技术在二维平面 问题中的应用时要用到，故在2 2 1 节先讲述四节点等参单元的基本推导。 2 2 1 四节点等参单元 四节点等参单元刚，采用比常应变的三角形单元次数更高的位移模式，可 l o 第2 章有限元基本原理及其在大地测量反演中的应用 以更好地反映弹性体中的位移状态和应力状态。 首先建立规整形状的母单元，如取边长为2 的正方形单元( 如图所示) ，其 形心处设局部坐标，得到如下的形函数，作坐标变换： ( a ) 母单元 亏 4 y - - ( f ，n ) y ， ( 2 1 ) ，- i y ，v ( b ) 等参元 图2 1 四结点平面等参元 使得图4 a 中的孝，7 平面上的4 个角点分别映射成图4 b 中叫平面上的4 个 点，其坐标为而，乃( i = l ，2 ，3 ，4 ) 。由于形函数m 是双线性的，善，7 平面上的 正方形被映射到砂平面上以而，y l 为角点的四边形。所以坐标变换式起了把砂 平面上的所有四边形单元( 称子单元) 都映射到勃平面上的正方形单元( 称母 单元) 的作用。同时，还可以把和刁坐标看成为子单元的局部坐标( 图4 b ) ， 该局部坐标系是用一组不超过1 的无量纲数来定出单元中的点，单元各边的方 程分别是f = 1 和刁= 1 。 位移模式为； 材= 口i + 口2 善+ 口3 叩+ 口4 f ，7l ( 2 2 ) ，= 口5 + 口6 善+ 口7 r l + ao 勿j 将结点的局部坐标值代入上式，列出四个结点处八个位移分量表达式，可 获两组四元联立方程，解出未知参数q 。q 。，回代后得到用结点位移表示的位 移模式 x 、j玎 芦j i ，- ， m = x 第2 章有限元基本原理及其在大地测量反演中的应用 甜= 嗽1 ，= ( 2 3 ) 闫鲥 式中 n i = ( 1 + 彘x i + r o ) 4 ( 2 4 ) 其中，磊= ，r o = 刁协，( ，= 1 ，2 ，3 ，4 ) 即： l = 丢( 1 吲( 1 刊= 丢( 1 一善一，7 + 矧 2 = 丢( 1 瑚( 1 刊= 丢( 1 + 善一刁一c v ) 3 = 丢( 1 + 似l + 功= 丢( 1 + 善+ ，7 + 矧 = 丢( 1 吲( 1 训= 丢( 1 一f + ，7 一俐 矩形单元的位移模式比常应变三角形单元中采用的线性位移模式增加翻项( 即 相当于砂项) ，故称之为双线性模式位移模式选定以后，即可按确定的公式 来推导单元刚度矩阵。 2 2 2 弹性力学基本方程与虚功原理【4 2 】 1 弹性力学基本方程 弹性力学基本方程包括平衡方程，几何方程，物理方程。对于一般的 各向同性弹性体其基本方程用张量形式表示如下： 1 ) 平衡方程 嘞，+ e = 0 ( 2 5 ) 2 ) 几何方程 = 寺( 甜“+ “) ( 2 6 ) 3 ) 物理方程 勺= 去【( 1 + p ) 一p 仃肚岛】 ( 2 7 ) 上式对应于平面应力问题，对于平面应变问题用e ii 告，。= 代替 即可 1 2 第2 章有限元基本原理及其在大地测量反演中的应用 2 边界条件 边界条件包括应力边界条件和位移边界条件和混合边界条件。 性力学的第一类、第二类和第三类边值问题。 1 ) 应力边界条件 o u r t i2 f 2 ) 位移边界条件 材，= 瓦 3 ) 混合边界条件 刀，= z ( 在砖上) 约= 玩( 在s 。上) 分别称为弹 ( 2 8 ) ( 2 9 ) ( 2 1 0 ) ( 2 1 1 ) 3 弹性体虚功原理 虚功原理可以表述为：在弹性体上，外力在任意一组几何可能位移上所做 的功，等于任意一组静力可能的应力在与上述几何可能位移所对应的应变上所 做的功。对于任意可能的静力位移材? i 和静力可能应力仃；及所对应的应变占：虚 功原理可以用下式表示： f f , # d v + u , d s = l 菇d y ( 2 1 2 ) 2 2 3 具体实现过程 上述的三角形单元和矩形单元采用的线性或双线性模式，是对实际位移场 的最低级逼近，精度受局限，而且矩形单元难以适应不规则边界。因此有必要 构造出具有较高精度的、能适应不规则边界的曲边或直边的四边形单元。 如果选择单元的插值函数是协调而完备的，有限单元法的解是收敛的，应 用有限单元法分析问题的过程【2 0 1 ，主要有以下几个步骤： 1 结构的离散化 它是有限元法分析的第一步，也是有限元法的基础。这个过程的主要任务 是选择网格划分的形式，确定较优的网格划分方案。本文网格划分时使用的是 a n s y s 软件，进行的平面四节点划分形式。 2 选择位移函数 为了能用节点位移来表示单元内任意一点的位移、应变和应力，首先假定单 第2 章有限元基本原理及其在大地测量反演中的应用 元内任意一点的位移是坐标的某种简单函数，称之为位移函数。即式( 2 4 ) 形式。 3 分析单元的力学特征 这个过程主要分析节点位移表示的单元应变和单元应力，即单元应变与结点 位移的关系矩阵，称应变矩阵 b 。 坩= 陋弦广 ( 2 1 3 ) 再由材料的本构关系( 即物理方程) ，得到单元弹性矩阵 d ，从而推出用结 点位移表示单元应力表达式： p ) | = 【d ) = 【d ) 。= p 广 ( 2 1 4 ) 其中， s = d l b 。 然后考虑结点平衡求得单元结点力与结点位移的关系，由矩阵 k 。表示， 称单元刚度矩阵。根据虚功原理或最小势能原理( 平衡条件) ，也可导出用结点 位移表示结点力的表达式： 扩) = 陋r p p 妣p 广= 陆r p 广 ( 2 1 5 ) 其中，单元刚度矩阵： 陆r = 陋r 【d 】陋k 哟脑蛊陋r 【d 】陋】矿 ( 2 1 6 ) 4 计算等效节点载荷 用有限单元法分析问题时是假定力，通过节点从一个单元传递到另一个单 元的，而实际上力是通过单元的公共边来传递的，所以作用于单元上的各种力 就必须等效地作用到节点上。 5 整体分析 由单元刚度矩阵求整体刚度矩阵，并和等效节点载荷一起构成整体平衡方 程，即l k 8 篁f ( 2 1 7 ) 式中k 为结构的总体刚度矩阵，艿为结构的节点位移列矩阵，为结构的 等效节点载荷列矩阵。 6 应用位移边界条件 利用位移边界条件( 式2 9 ) 消除整体刚度矩阵的奇异性，使整体平衡方程可 解 7 求解结构平衡方程 结构的平衡方程是以总体刚度矩阵为系数的线性代数方程组，解方程组可 1 4 第2 章有限元基本原理及其在大地测量反演中的应用 以求得未知的节点位移。 8 计算单元应力 由节点位移和弹性力学的几何方程可以求出单元应变，然后由单元应变和 弹性力学的物理方程可以计算出单元应力 2 3 有限单元法反演边界力原理 对于线弹性或者粘弹性的材料，在小变形的情况下，根据叠加原理，力可 以叠加，位移模型是线性的，位移也可以叠加。利用这一原理可采用叠加的方 法进行构造应力场的反演。假设有一构造框架，有m 个受力边界，每个受力边 界在个方向受力，则用有限单元法反演构造应力场的做法如下： 在第一个受力边界的第一个方向作用单位力，利用有限单元法正演，可以 得到一个位移场和一个应力场，同样对其它的受力边界的同一方向均分别施加 单位力，这样可以得到m 个位移场和肘个应力场。对于地面的一个测量点在该 方向的位移值为s ，则由叠加原理可知有下式成立： 材 嘞岛= s j ( 2 1 8 ) ，1 1 上式中口。为测量点，如果在第f 个方向要产生位移s ，在第f 个受力边界的第 一个方向所需加载力的大小，为待求未知量：s , j 为在第，个受力边界的第个方 向加载单位力时该方向上的位移值。由于测量数据为含有误差的量，上式可改 写为： 吖 1 ，= y 口。s 。一s ， ( 2 1 9 ) 这样，”1 由一个测量点在个方向的位移可以得到个方程。未知数的数量 为m n 个，所以至少需要知道m 个点在各个方向的位移，才能确定各个边的 边界力在各个方向的分量。如果测量点数量多于m 个，则用最小二乘平差法求 未知量，这样就得到了边界力在各个方向的分量f 4 ”。由边界力在各个方向的分 量就可以唯一确定各个边界力的大小和方向。如果用测量点的应力值进行边界 力的反演，其过程和使用的公式与用位移数据进行反演一样，只是使用的数据 不同。 第3 章分割结点技术及其在大地测量反演中的廑旦 第3 章分割结点技术及其在大地测量反演中的应用 3 1 分割结点技术原理 利用有限元方法解决大地测量或地球物理的反演问题时，首先将一个连续 体用假想的剖面分成有限个单元，通过构造插值函数，利用结点的平衡条件求得 结点上位移的近似值，进而求出单元的应变和应力的近似值。因此，对于一个有 限元问题，其解算步骤概括起来是由离散化、单元分析以及总体分析三个主要步 骤构成。一般情况下，由不同的单元所共有的结点其位移是相等的。但被“分割 的结点，则属于这样一些特殊的结点，即与这些结点的位移和这些结点所在的单 元有关。例如，对于某