大地形变测量学课程设计及实习报告区域地壳应变的简单分析大学论文Word下载.docx

1、2. 通过编程实现不同尺度应变值归化，并通过计算分析，讨论不同图形尺度应变值的特点及归化的原理与方法。完成一篇课程设计论文区域地壳应变的简单分析*（中国武汉 武汉大学）摘要 主要介绍了进行区域地壳应变分析的相关方法及相关分析.研究利用有限单元应变分析法求解应变值的尺度性对性问题，探讨了根据叠加原理进行有限单元法中区域单元合理划分的问题.从统计分析入手发现了应变值与计算图形尺度存在显著的相关性，据此提出进行应变值尺度归化的必要性，给出了归化方法.关键词 区域地壳应变 有限单元法 尺度相对性 归化引言应变分析中所用的基本数据有两种：一是新、旧测量的原始资料，二是根据新、旧测量结果的比较得出的位移场

2、。利用原始观测资料时，要求新、旧测量的网形和观测量都相同，因而其适用性受到一些限制。但它的优点是，不依赖其他的观测量，避免了监测网平差中因基准点设定不当等原因带来的影响。利用位移场时，由于网中各点的位移的向量是根据新、旧平差结果的坐标之差得出的，为了使位移场能反映实际地壳应变，把残余的误差影响化为最小，必须采用特殊的平差方法，例如自由网平差和拟稳平差。利用位移场的优点是：所有的观测量都可用于应变分析，并不要求新、旧测量的观测量都相同，只是要求它们属于同一大地基准。此外，在求定位移场的平差过程中，可以滤掉观测数据的粗差和估计观测质量，而且由各点位移向量的图解，可以看出各点位移的趋势。由于这些原因

3、，人们侧重于利用位移场（许才军等，2006）。1、区域地壳应变分析的有限单元法1.1有限单元法的基本思想有限单元法的基本思想是将连续的求解区域离散成一组有限个、且按一定方式相互联结在一起的单元的组合体。由于单元能按不同的联结方式进行组合，且单元本身又有不同的形状，因此可以最大程度上模型化几何形状复杂的求解域。有限单元法作为数值分析方法的另一个重要特点是利用在每一个单元内假设的近似函数来分片地表示全求解域上的待求未知场函数。单元内的近似函数通常由未知场函数及其导数在单元各个结点的数值与其插值函数来表达。这样一来，一个问题的有限元分析中，未知场函数及其导数在各个结点上的数值就成为新的未知量，从而使

4、一个连续的无限自由度问题变成离散的有限自由度问题。一经求解出这些未知量，就可以通过插值函数计算出各个单元内场函数的近似值，从而得到整个求解域上的近似解。显然，随着单元数目的增加，也即单元尺寸的缩小，或者随着单元自由度的增加及插值函数精度的提高，解的近似程度将不断改进。如果单元是满足收敛要求的，近似解将收敛于精确解。1.2有限单元应变分析法采用有限单元分析法进行应变分析，就是把一个大的区域分割成一些有限的小区域。三角网和三边网的基本图形是三角形把它作为有限单元比较方便。分别对各三角形进行应变分析，就可以得到接近于真实的地壳应变情况。韦尔施（Welsch）曾用这种单元划分法分析了美国西部圣安德列斯

5、断层上霍利斯特附近的局部和区域应变形式（Welsch，1982）。黄立人（1999）在研究华北部分地区水平变形时，根据地区介质构造以及研究区内断裂优势走向分布和复杂程度，采用8节点6面体单元用北西向的24各不等间距剖面剖分整个区域。但是有限单元法是区域型方法,必须在全区域进行剖分,如果剖分后的单元和节点数目多,则得到的线性代数方程组很大,特别是三维问题和地球物理中常遇的无界区域问题,一般需要在中、大型计算机上才能完成有限元法的计算。大地测量反演构造应力场,也主要采用有限元进行数值模拟,反演结果的可靠性和准确度取决于多种因素,其中在构造应力场分布为非线性的区域,如断裂带、边界转折点处的单元划分也

6、是主要因素之一。如果这些特殊区域单元划分过大,将使反演结果不可靠;单元划分过小,又引起不必要的计算机内存和机时的浪费。因此,有限单元分析法中,合理地划分单元是一个十分重要的问题。1.3有限元的划分根据弹性理论的叠加原理,同时作用于物体的两组或两组以上外力的总和在物体内所产生的效果（应力,应变及位移）等于各组外力分别作用而产生的效果的总和。因此,对于一个确定了构造格架的地块来说,设有m个受力边界,并假设各边界所受边界力和地块的介质参数均已知,则可以正演计算该地块在总边界力（外力）作用下的位移场和应力场。同时,又可以分别在m个受力边界的各个方向（平面问题为x、y 2个方向,三维问题有x、y、z 3

7、个方向）作用单位边界力,在地块介质参数不变的情况下,仍可正演计算地块在各个单位边界力作用下的位称场和应力场。设为地块单元结点k的总位移值,为地块单元结点k在第i边界第j方向作用单位力后的位移值,为第i边界第j方向的边界力系数,由叠加原理可得: （1.1）其中n为空间维数,对于平面问题,n=2,立体问题n=3。如果,分别代表k结点在总边界（外）力和单位边界（外）力作用下获得的某一应力分量或主应力值,（1）式也成立,它们是等价的。当单元划分不合理时,可认为单元结点位移值有误差,由（1）式可得的误差方程式: （1.2）写成矩阵形式为:V = A- L （1.3）其中V为改正数向量,A为系数矩阵,为边

8、界力系数列阵,L为常数项列阵。由（3）可得边界力系数的最小二乘解: （1.4）由此可得边界力F: （1.5）如果有限单元格网划分特别是断裂区单元划分合适,则有限元法计算结果精度高,利用单元结点位移值,不论是断裂区还是非断裂区单元结点位移,都能较好地反演出边界力系数,使之与其理论值相符。反之,若反演的边界力系数相差较大且都与理论值有差别,则说明结点位移值计算精度不高,或者说,断裂区单元结点位移精度与非断裂区单元结点位移精度不匹配,这可认为主要是由断裂区单元划分不合适而引起的。需要重新划分断裂区单元,直至由断裂区单元结点位移和非断裂区单元结点位移反演的边界力系数都与理论值趋于一致。二、应变值的尺度

9、相对性问题2.1不同尺度应变值的分布特征利用华北GPS监测网1993,1995年观测资料来简单分析图形单元应变值与图形尺度之间的关系. 以华北GPS监测网1992, 1995年观测资料为例, 分析图形单元应变值与图形单元尺度之间的关系. 由GPS测点直接构成53个最小图形单元（张希等, 1998） ; 为了更好地反映统计特征, 我们又取更大一些点距构成27个图形单元, 以这80个不同图形尺度的应变 作为样本, 分别以最大剪应变值与单元图形尺度为纵、横坐标, 得到了最大剪应变-图形尺度分布图如下：. 图2.1 剪应变-图形尺度分布图从图中可以看出两个特点: 一是一些高值点基本上分布于图形尺度较小

10、的区间; 二是图形尺度小的区间内的应变值分布更离散（量值变化幅度更大） . 利用幂指函数模型进行非线性回归, 可求解得回归方程： （1.6）式中 ,应变值数值单位取（ mm /km,即微应变）,d为计算图形尺度（单位为km）,取d = S（S 为计算图形面积）.计算出的复相关系数R= 0.575, 远大于临界值,表明应变与尺度之间的相关性显著.若以图形尺度相近的每8个图形的最大剪应变值取平均,得到最大剪应变均值分布图（图2）.以这组应变均值为采样点进行回归计算,得到幂指数回归方程 （1.7）相关系数R= 0.957如图： 图2.2 8个图形尺度相近的取平均后回归方程为了反映不同尺度应变值离散程

11、度的差异, 我们又以式（ 1）来确定不同图形尺度应变值的统计平均值, 并分组计算了均方差, 得到应变值均方差-图形尺度分布图如下： 图2.3 应变值均方差-图形尺度分布图由图2.3可见，应变值的均方差也是呈幂指数衰减趋势的.2.2应变值与图形尺度相关的必然性在应用三角形直接在椭球上计算应变时，根据公式我们可以直接看到计算的结果与三角形的大小有明显关系。另外，统计结果也反映了这样的趋势：一是高值点基本分布于图形尺度较小的（定于图形尺度为，单位为km，s为所在三角形的面积）的区间；二是图形尺度较小的区间应变值分布更离散，量值变化幅度更大（江在森等，2000）。根据弹性力学理论，用大地形变资料计算应

12、变时，实际上假定观测资料构成的最小图形内的平均应变值。因此大地形变资料计算的应变只能是相对测点分布的近似平均值。由于地壳应变的空间分布不均匀，在测点密集跨度小的情况下，对实际应变场才具有较高的分辨率，而跨度较大的图形就反映不出细部的剧烈变化。另外，从地壳运动特征来看，也会使小尺度应变值具有分布更离散的特点。因此，有必要分析尺度的相对性并且进行归化。三、由大地形变资料求解应变值的尺度归化3.1归化方法介绍把不同尺度的应变看成具有不同数学期望和不同方差母体随机变量，就可以采用将一般正态分布化为标准正态分布的方法对不同尺度应变值的分布，我们分别表示为: （3.1）令： （3.2）对上式扩充一个加常数

13、和一个乘常数则得到： （3.3）其中加常数u0表示归化某一标准尺度的应变值的数学期望，表示某一标准尺度的均方差，这样就可以把不同尺度的应变值从统计分布特征上联系起来了。从而归化处理的计算公式可写为：（3.4）其中：为归化后的应变值；为原应变值；为某一标准尺度应变值的统计平均值；为计算尺度（与相对应的尺度）应变值的统计平均值；m0为标准尺度应变值的统计均方差；为计算尺度应变值的统计均方差（许才军等，2006）。4.2归化步骤总结由4.1可总结出归化需要以下几个步骤：（1）、根据实际资料计算不同尺度应变的统计平均值，亦可用拟合得到的经验函数值确定；（2）分别计算不同尺度的应变值的均方差（3）选定某

14、一归化尺度，这一尺度的应变统计平均值为，均方差为（4）根据给定公式进行归化。对于同一个地区相同的观测资料，当观测样本较多时，可以由应变观测值用幂函数模型进行回归，得出统一回归方程，这样其应变值的统计平均值和均方差都可由回归方程来计算。参 考 文 献1许才军，申文斌，晁定波.地球物理大地测量学原理与方法M.武汉：武汉大学出版社.2006，208-209.2许才军，张朝玉.地壳形变测量与数据处理M.武汉：武汉大学出版社，2009,167-169.3黄立人，马青.华北部分地区水平变形的力学机制J.地震学报,1999,21（1）:50-56.4江在森.地形变资料求解应变值的尺度相对性研究J.地震学报,

15、2000,4（22）,352-359.5许才军，晁定波.有限元分析中断裂区域单元的划分J.地壳形变与地震,1997,17（3）.6刘孙君.GPS约束下青藏高原地壳运动位移场模拟及应力应变分析D.中国科学院，2003.四、附录：程序源代码（采用Fortran语言编写）1、 有限单元划分不合理时，进行平差求边界力系数的最小二乘解：PROGRAM LGCIMPLICIT NONEreal,allocatable:a（:,:）,b（:）,R（:）,R0（:）,L（:）,C（:）,D（:）integer m,n,i,jprint*,请输入矩阵行数、列数m,n（m=n）:read*,mread*,nall

16、ocate（a（m,n）,b（n,m）,R0（n,n）,R（n,n）,L（m,1）,C（n,m）,D（n,1）open（1,file=juzhen.txt,form=formatteddo i=1,m read（1,*）a（i,:enddoclose（1）open（2,file=changshuxiang.txt read（2,*）L（i,1）close（2）call zhuanzhi（a,b,m,n）call cheng（b,a,R0,n,m,n）call inv（R0,R,n）call cheng（R,b,C,n,m,n）call cheng（C,L,D,n,m,1）open（3,file

17、=result.txtdo i=1,n write（3,*）D（i,:close（3）print*,R0print*,Ddeallocate（a,b,R0,R,L,C,D）end program lgc!-subroutine zhuanzhi（V,Vt,m,n） !求矩阵的转置real*8:V（m,n）,Vt（n,m） do j=1,n Vt（j,i）=V（i,j） enddoend subroutine zhuanzhisubroutine cheng（P,Q,R,m1,n,m2） !矩阵相乘integer n,m1,m2,i,jP（m1,n）,Q（n,m2）,R（m1,m2）do i=1

18、,m1 do j=1,m2 R（i,j）=0 do k=1,n R（i,j）=R（i,j）+P（i,k）*Q（k,j）end subroutine cheng subroutine inv（A,invA,N） !计算逆矩阵implicit noneinteger ninteger iA（n,n）,invA（n,n）,E（n,n）E=0 do i=1,n !设置E为单位矩阵 E（i,i）=1end docall mateq（A,E,invA,N,N）end subroutine invsubroutine mateq（A,B,X,N,M） !高斯列主元消去法计算矩阵方程 AX=Binteger

19、N,M,iA（N,N）,B（N,M）,X（N,M）btemp（N）,xtemp（N）do i=1,M btemp=B（:,i） call elgauss（A,btemp,xtemp,N） X（:,i）=xtempend subroutine mateqsubroutine elgauss（A,b,x,N） !高斯列主元消去法 Ax=binteger i,k,Ninteger id_max !主元素标号A（N,N）,b（N）,x（N）Aup（N,N）,bup（N）Ab（N,N+1） !Ab为增广矩阵 Abvtemp1（N+1）,vtemp2（N+1）Ab（1:N,1:N）=AAb（:,N+1）=

20、bdo k=1,N-1 elmax=dabs（Ab（k,k） id_max=k do i=k+1,n if （dabs（Ab（i,k）elmax） then elmax=Ab（i,k） id_max=i end if end do vtemp1=Ab（k,: vtemp2=Ab（id_max,: Ab（k,:）=vtemp2 Ab（id_max,:）=vtemp1 do i=k+1,N temp=Ab（i,k）/Ab（k,k） Ab（i,:）=Ab（i,:）-temp*Ab（k,:Aup（:）=Ab（1:N）bup（:）=Ab（:,N+1）call uptri（Aup,bup,x,n）end

21、subroutine elgausssubroutine uptri（A,b,x,N） ! 上三角方程组的回带方法Ax=binteger i,j,Nx（N）=b（N）/A（N,N） !回带部分do i=n-1,1,-1 x（i）=b（i） do j=i+1,N x（i）=x（i）-a（i,j）*x（j） x（i）=x（i）/A（i,i）end subroutine uptri截图： 图4.1 系数矩阵截图 图4.2 平差求边界力系数的最小二乘解运行截图2、应变值的尺度归化PROGRAM GUIHUAinteger u0,sigma0T（:）,u（:）,sigma（:）,T1（:integer

22、i,j,k请输入标准尺度的数学期望和均方差read*,u0,sigma0请输入待归化的一般正态分布的个数:if（n2） pauseallocate（T（n）,u（n）,sigma（n）,T1（n）请依次输入各组应变值、期望、方差read*,T（i）,u（i）,sigma（i）do j=1,nT1（j）=（T（j）-u（j）/sigma（j）T1（j）=u0+（T（j）-u（j）*sigma0/sigma（j）原应变值, 归化后的应变值do k=1,nprint*,T（k）,T1（k）deallocate（T,u,sigma,T1）end 图4.3 应变值的尺度归化程序运行截图课程设计感想一转眼

23、长达三周的实习就要结束了，虽然中间经历了颇多挫折，但临结束却还是有几分不舍，或许是因为已经大四了。言归正传，下面谈一下这次实习的大概经理与自己的一些感受吧。最初知道这个实习长达三周后，其实我的内心是拒绝的，因为作为一名考研党，在这三周里要分出大量的时间去机房，必然减少了复习的时间，再加上几门选修课穿插其间，心里那个苦的呦！再回顾这三周的实习工作，感受很深，收获颇丰，虽然大部分时间都是在机房度过，中间难免在教学楼、图书馆之间来回奔波，但每天都过得很充实，而且在本次课程设计的过程中踏踏实实的学到了很多东西，这是在课本中找不到的。伴随着大一新生开学的欢声笑语，在9月5日，我们正式开始了为期三周的课程

24、设计，前一天老师已经为每个人分好了题目，我的题目是“区域地壳应变分析”，首先要做的就是查找资料，几乎花费了我一周时间，细细的看，细细的筛选，查找的过程总是枯燥乏味，很多不理解的公式、定理也要慢慢去思索，比如说弹性力学中的叠加原理和最小位能原理，不过我觉得这个过程才是提高自己最多的，首先是寻找有用信息的能力大大提高了，其次查找的资料和论文上的东西让我对以前学习到的东西理解的更为透彻明朗，了解很多以前自己不知道的东西，除此之外，我们还去了珞珈山人防隧道去参观了几个实验室，那天领队的汪老师为我们讲解了很多，在那潮湿阴暗的山洞中同学们对新事物的向往却如同烈火一般熊熊燃烧着，大家还在山洞前合了影，算是本

25、次实习最让人难以忘怀的记忆了。第二周是在浓浓的编程氛围中度过的，不仅让我重新熟练了好久没用的Fortran，提高了自己的编程能力，更让我深刻认识到应用编程去接解决实际问题的困难所在，大大增强了自己动手去解决问题的能力。虽然这次课程只有那么短暂的3周时间，我感觉到这些天我的所学胜过以前很多，更是为不久之后的毕业设计打下一些基础，这次任务原则上是设计，其实就是一次大的作业，是让我对课本知识的巩固和对基本公式的熟悉和应用，编程过程中那些繁琐的数据，使我做事的耐心和仔细程度得以提高。不仅如此，我还深刻地体会到做研究的不易和艰辛，这个过程一丝一毫都马虎不得，正如老师所说：“我们是科班出身”。总之，这次课程设计使我收获很多、学会很多、比以往更有耐心很多。感谢学校及老师给我们这次课程设计的机会，最真挚的感谢各位辅导老师 ，在设计过程中，老师精心的辅导和不厌其烦地的态度才使得我以顺利的完成这次设计，他那无私的奉献的精神照耀着我们对学习的热爱，同时也增加我们对知识的追求和欲望度。实习成绩评定表综合评语：口试成绩所占比例30%报告成绩70%总评成绩 评阅指导教师： 年 月 日