區(qū)域地殼應(yīng)變的簡單分析

1、課程編號：1300894（課頭號：20161024137） 課程性質(zhì)：必修大地形變測量學(xué)課程設(shè)計及實習(xí)實 習(xí) 報 告 學(xué)院： 測繪學(xué)院 專業(yè)： 測繪工程 地點： 武漢大學(xué) 班級： 2 0 1 3 0 5 組號： A方向 姓名： * 學(xué)號： 2013301610* 教師： 許才軍 賈建鋼 汪建軍 2016年9月5日 至 2016年 9月23日課程設(shè)計與實習(xí)概述一、課程設(shè)計與實習(xí)目的 1.對課程中所學(xué)到的知識有更深的理解，做到理論聯(lián)系實際，通過課程設(shè)計與實習(xí)鞏固所學(xué)知識。 2. 培養(yǎng)學(xué)生的編程能力、實際動手能力、分析問題以及解決復(fù)雜工程問題的能力。 3.&#

2、160;為畢業(yè)設(shè)計打下基礎(chǔ)，使學(xué)生對于論文的寫作和科研的過程有一定認(rèn)識。  二、課程設(shè)計題目名稱 區(qū)域地殼應(yīng)變分析三、課程設(shè)計要求 1. 查閱相關(guān)文獻, 深入理解區(qū)域地殼應(yīng)變分析基本原理與方法，推導(dǎo)給出相關(guān)數(shù)學(xué)模型2. 通過編程實現(xiàn)不同尺度應(yīng)變值歸化，并通過計算分析，討論不同圖形尺度應(yīng)變值的特點及歸化的原理與方法。3. 完成一篇課程設(shè)計論文區(qū)域地殼應(yīng)變的簡單分析*（中國武漢 武漢大學(xué)）摘要 主要介紹了進行區(qū)域地殼應(yīng)變分析的相關(guān)方法及相關(guān)分析.研究利用有限單元應(yīng)變分析法求解應(yīng)變值的尺度性對性問題，探討了根據(jù)疊加原理進行有限單元法中區(qū)域單

3、元合理劃分的問題.從統(tǒng)計分析入手發(fā)現(xiàn)了應(yīng)變值與計算圖形尺度存在顯著的相關(guān)性，據(jù)此提出進行應(yīng)變值尺度歸化的必要性，給出了歸化方法.關(guān)鍵詞 區(qū)域地殼應(yīng)變 有限單元法 尺度相對性 歸化引言應(yīng)變分析中所用的基本數(shù)據(jù)有兩種：一是新、舊測量的原始資料，二是根據(jù)新、舊測量結(jié)果的比較得出的位移場。利用原始觀測資料時，要求新、舊測量的網(wǎng)形和觀測量都相同，因而其適用性受到一些限制。但它的優(yōu)點是，不依賴其他的觀測量，避免了監(jiān)測網(wǎng)平差中因基準(zhǔn)點設(shè)定不當(dāng)?shù)仍驇淼挠绊?。利用位移場時，由于網(wǎng)中各點的位移的向量是根據(jù)新、舊平差結(jié)果的坐標(biāo)之差得出的，為了使位移場能反映實際地殼應(yīng)變，把殘余的誤差影響化為最小，必須采用特殊的平

4、差方法，例如自由網(wǎng)平差和擬穩(wěn)平差。利用位移場的優(yōu)點是：所有的觀測量都可用于應(yīng)變分析，并不要求新、舊測量的觀測量都相同，只是要求它們屬于同一大地基準(zhǔn)。此外，在求定位移場的平差過程中，可以濾掉觀測數(shù)據(jù)的粗差和估計觀測質(zhì)量，而且由各點位移向量的圖解，可以看出各點位移的趨勢。由于這些原因，人們側(cè)重于利用位移場（許才軍等，2006）。1、區(qū)域地殼應(yīng)變分析的有限單元法1.1有限單元法的基本思想有限單元法的基本思想是將連續(xù)的求解區(qū)域離散成一組有限個、且按一定方式相互聯(lián)結(jié)在一起的單元的組合體。由于單元能按不同的聯(lián)結(jié)方式進行組合，且單元本身又有不同的形狀，因此可以最大程度上模型化幾何形狀復(fù)雜的求解域。有限單元法

5、作為數(shù)值分析方法的另一個重要特點是利用在每一個單元內(nèi)假設(shè)的近似函數(shù)來分片地表示全求解域上的待求未知場函數(shù)。單元內(nèi)的近似函數(shù)通常由未知場函數(shù)及其導(dǎo)數(shù)在單元各個結(jié)點的數(shù)值與其插值函數(shù)來表達。這樣一來，一個問題的有限元分析中，未知場函數(shù)及其導(dǎo)數(shù)在各個結(jié)點上的數(shù)值就成為新的未知量，從而使一個連續(xù)的無限自由度問題變成離散的有限自由度問題。一經(jīng)求解出這些未知量，就可以通過插值函數(shù)計算出各個單元內(nèi)場函數(shù)的近似值，從而得到整個求解域上的近似解。顯然，隨著單元數(shù)目的增加，也即單元尺寸的縮小，或者隨著單元自由度的增加及插值函數(shù)精度的提高，解的近似程度將不斷改進。如果單元是滿足收斂要求的，近似解將收斂于精確解。1.

6、2有限單元應(yīng)變分析法采用有限單元分析法進行應(yīng)變分析，就是把一個大的區(qū)域分割成一些有限的小區(qū)域。三角網(wǎng)和三邊網(wǎng)的基本圖形是三角形把它作為有限單元比較方便。分別對各三角形進行應(yīng)變分析，就可以得到接近于真實的地殼應(yīng)變情況。韋爾施（Welsch）曾用這種單元劃分法分析了美國西部圣安德列斯斷層上霍利斯特附近的局部和區(qū)域應(yīng)變形式（Welsch，1982）。黃立人（1999）在研究華北部分地區(qū)水平變形時，根據(jù)地區(qū)介質(zhì)構(gòu)造以及研究區(qū)內(nèi)斷裂優(yōu)勢走向分布和復(fù)雜程度，采用8節(jié)點6面體單元用北西向的24各不等間距剖面剖分整個區(qū)域。但是有限單元法是區(qū)域型方法,必須在全區(qū)域進行剖分,如果剖分后的單元和節(jié)點數(shù)目多,則得到的

7、線性代數(shù)方程組很大,特別是三維問題和地球物理中常遇的無界區(qū)域問題,一般需要在中、大型計算機上才能完成有限元法的計算。大地測量反演構(gòu)造應(yīng)力場,也主要采用有限元進行數(shù)值模擬,反演結(jié)果的可靠性和準(zhǔn)確度取決于多種因素,其中在構(gòu)造應(yīng)力場分布為非線性的區(qū)域,如斷裂帶、邊界轉(zhuǎn)折點處的單元劃分也是主要因素之一。如果這些特殊區(qū)域單元劃分過大,將使反演結(jié)果不可靠;單元劃分過小,又引起不必要的計算機內(nèi)存和機時的浪費。因此,有限單元分析法中,合理地劃分單元是一個十分重要的問題。1.3有限元的劃分根據(jù)彈性理論的疊加原理,同時作用于物體的兩組或兩組以上外力的總和在物體內(nèi)所產(chǎn)生的效果(應(yīng)力,應(yīng)變及位移)等于各組外力分別作用

8、而產(chǎn)生的效果的總和。因此,對于一個確定了構(gòu)造格架的地塊來說,設(shè)有m個受力邊界,并假設(shè)各邊界所受邊界力和地塊的介質(zhì)參數(shù)均已知,則可以正演計算該地塊在總邊界力(外力)作用下的位移場和應(yīng)力場。同時,又可以分別在m個受力邊界的各個方向(平面問題為x、y 2個方向,三維問題有x、y、z 3個方向)作用單位邊界力,在地塊介質(zhì)參數(shù)不變的情況下,仍可正演計算地塊在各個單位邊界力作用下的位稱場和應(yīng)力場。設(shè)為地塊單元結(jié)點k的總位移值,為地塊單元結(jié)點k在第i邊界第j方向作用單位力后的位移值,為第i邊界第j方向的邊界力系數(shù),由疊加原理可得: (1.1)其中n為空間維數(shù),對于平面問題,n=2,立體問題n=3。如果,分別

9、代表k結(jié)點在總邊界(外)力和單位邊界(外)力作用下獲得的某一應(yīng)力分量或主應(yīng)力值,(1)式也成立,它們是等價的。當(dāng)單元劃分不合理時,可認(rèn)為單元結(jié)點位移值有誤差,由(1)式可得的誤差方程式: (1.2)寫成矩陣形式為:V = A- L (1.3)其中V為改正數(shù)向量,A為系數(shù)矩陣,為邊界力系數(shù)列陣,L為常數(shù)項列陣。由(3)可得邊界力系數(shù)的最小二乘解: (1.4)由此可得邊界力F: (1.5)如果有限單元格網(wǎng)劃分特別是斷裂區(qū)單元劃分合適,則有限元法計算結(jié)果精度高,利用單元結(jié)點位移值,不論是斷裂區(qū)還是非斷裂區(qū)單元結(jié)點位移,都能較好地反演出邊界力系數(shù),使之與其理論值相符。反之,若反演的邊界力系數(shù)相差較大且

10、都與理論值有差別,則說明結(jié)點位移值計算精度不高,或者說,斷裂區(qū)單元結(jié)點位移精度與非斷裂區(qū)單元結(jié)點位移精度不匹配,這可認(rèn)為主要是由斷裂區(qū)單元劃分不合適而引起的。需要重新劃分?jǐn)嗔褏^(qū)單元,直至由斷裂區(qū)單元結(jié)點位移和非斷裂區(qū)單元結(jié)點位移反演的邊界力系數(shù)都與理論值趨于一致。二、應(yīng)變值的尺度相對性問題2.1不同尺度應(yīng)變值的分布特征利用華北GPS監(jiān)測網(wǎng)1993,1995年觀測資料來簡單分析圖形單元應(yīng)變值與圖形尺度之間的關(guān)系. 以華北GPS監(jiān)測網(wǎng)1992, 1995年觀測資料為例, 分析圖形單元應(yīng)變值與圖形單元尺度之間的關(guān)系. 由GPS測點直接構(gòu)成53個最小圖形單元(張希等, 1998) ; 為了更好地反映統(tǒng)

11、計特征, 我們又取更大一些點距構(gòu)成27個圖形單元, 以這80個不同圖形尺度的應(yīng)變 作為樣本, 分別以最大剪應(yīng)變值與單元圖形尺度為縱、橫坐標(biāo), 得到了最大剪應(yīng)變-圖形尺度分布圖如下：. 圖2.1 剪應(yīng)變-圖形尺度分布圖從圖中可以看出兩個特點: 一是一些高值點基本上分布于圖形尺度較小的區(qū)間; 二是圖形尺度小的區(qū)間內(nèi)的應(yīng)變值分布更離散(量值變化幅度更大) . 利用冪指函數(shù)模型進行非線性回歸, 可求解得回歸方程： (1.6)式中 ,應(yīng)變值數(shù)值單位取( mm /km,即微應(yīng)變),d為計算圖形尺度(單位為km),取d = S(S 為計算圖形面積).計算出的復(fù)相關(guān)系數(shù)R= 0.575, 遠大于臨界值,表明應(yīng)

12、變與尺度之間的相關(guān)性顯著.若以圖形尺度相近的每8個圖形的最大剪應(yīng)變值取平均,得到最大剪應(yīng)變均值分布圖(圖2).以這組應(yīng)變均值為采樣點進行回歸計算,得到冪指數(shù)回歸方程 (1.7)相關(guān)系數(shù)R= 0.957如圖： 圖2.2 8個圖形尺度相近的取平均后回歸方程為了反映不同尺度應(yīng)變值離散程度的差異, 我們又以式( 1)來確定不同圖形尺度應(yīng)變值的統(tǒng)計平均值, 并分組計算了均方差, 得到應(yīng)變值均方差-圖形尺度分布圖如下： 圖2.3 應(yīng)變值均方差-圖形尺度分布圖由圖2.3可見，應(yīng)變值的均方差也是呈冪指數(shù)衰減趨勢的.2.2應(yīng)變值與圖形尺度相關(guān)的必然性在應(yīng)用三角形直接在橢球上計算應(yīng)變時，根據(jù)公式我們可以直接看到計

13、算的結(jié)果與三角形的大小有明顯關(guān)系。另外，統(tǒng)計結(jié)果也反映了這樣的趨勢：一是高值點基本分布于圖形尺度較小的（定于圖形尺度為 ，單位為km，s為所在三角形的面積）的區(qū)間；二是圖形尺度較小的區(qū)間應(yīng)變值分布更離散，量值變化幅度更大（江在森等，2000）。根據(jù)彈性力學(xué)理論，用大地形變資料計算應(yīng)變時，實際上假定觀測資料構(gòu)成的最小圖形內(nèi)的平均應(yīng)變值。因此大地形變資料計算的應(yīng)變只能是相對測點分布的近似平均值。由于地殼應(yīng)變的空間分布不均勻，在測點密集跨度小的情況下，對實際應(yīng)變場才具有較高的分辨率，而跨度較大的圖形就反映不出細部的劇烈變化。另外，從地殼運動特征來看，也會使小尺度應(yīng)變值具有分布更離散的特點。因此，有必

14、要分析尺度的相對性并且進行歸化。三、由大地形變資料求解應(yīng)變值的尺度歸化3.1歸化方法介紹把不同尺度的應(yīng)變看成具有不同數(shù)學(xué)期望和不同方差母體隨機變量，就可以采用將一般正態(tài)分布化為標(biāo)準(zhǔn)正態(tài)分布的方法對不同尺度應(yīng)變值的分布，我們分別表示為: (3.1)令： (3.2)對上式擴充一個加常數(shù)和一個乘常數(shù)則得到： (3.3)其中加常數(shù)u0表示歸化某一標(biāo)準(zhǔn)尺度的應(yīng)變值的數(shù)學(xué)期望，表示某一標(biāo)準(zhǔn)尺度的均方差，這樣就可以把不同尺度的應(yīng)變值從統(tǒng)計分布特征上聯(lián)系起來了。從而歸化處理的計算公式可寫為：(3.4)其中：為歸化后的應(yīng)變值；為原應(yīng)變值；為某一標(biāo)準(zhǔn)尺度應(yīng)變值的統(tǒng)計平均值；為計算尺度（與相對應(yīng)的尺度）應(yīng)變值的統(tǒng)計

15、平均值；m0為標(biāo)準(zhǔn)尺度應(yīng)變值的統(tǒng)計均方差；為計算尺度應(yīng)變值的統(tǒng)計均方差（許才軍等，2006）。4.2歸化步驟總結(jié)由4.1可總結(jié)出歸化需要以下幾個步驟：（1）、根據(jù)實際資料計算不同尺度應(yīng)變的統(tǒng)計平均值，亦可用擬合得到的經(jīng)驗函數(shù)值確定；（2）分別計算不同尺度的應(yīng)變值的均方差（3）選定某一歸化尺度，這一尺度的應(yīng)變統(tǒng)計平均值為，均方差為（4）根據(jù)給定公式進行歸化。對于同一個地區(qū)相同的觀測資料，當(dāng)觀測樣本較多時，可以由應(yīng)變觀測值用冪函數(shù)模型進行回歸，得出統(tǒng)一回歸方程，這樣其應(yīng)變值的統(tǒng)計平均值和均方差都可由回歸方程來計算。參 考 文 獻1許才軍，申文斌，晁定波.地球物理大地測量學(xué)原理與方法M.武漢：武漢大

16、學(xué)出版社.2006，208-209.2許才軍，張朝玉.地殼形變測量與數(shù)據(jù)處理M.武漢：武漢大學(xué)出版社，2009,167-169.3黃立人，馬青.華北部分地區(qū)水平變形的力學(xué)機制J.地震學(xué)報,1999,21(1):50-56.4江在森.地形變資料求解應(yīng)變值的尺度相對性研究J.地震學(xué)報,2000,4(22),352-359.5許才軍，晁定波.有限元分析中斷裂區(qū)域單元的劃分J.地殼形變與地震,1997,17(3).6劉孫君.GPS約束下青藏高原地殼運動位移場模擬及應(yīng)力應(yīng)變分析D.中國科學(xué)院，2003.四、附錄：程序源代碼（采用Fortran語言編寫）1、 有限單元劃分不合理時，進行平差求邊界力系數(shù)的最

17、小二乘解：PROGRAM LGCIMPLICIT NONEreal,allocatable:a(:,:),b(:,:),R(:,:),R0(:,:),L(:,:),C(:,:),D(:,:)integer m,n,i,jprint*,'請輸入矩陣行數(shù)、列數(shù)m,n(m>=n):'read*,mread*,nallocate(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='formatted')do i=1,m read(1,*)a(

18、i,:)enddoclose(1)open(2,file='changshuxiang.txt',form='formatted')do i=1,m read(2,*)L(i,1)enddoclose(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='result.txt',form='formatted')do i=1,n write(

19、3,*)D(i,:)enddoclose(3)print*,R0print*,Ddeallocate(a,b,R0,R,L,C,D)end program lgc!-subroutine zhuanzhi(V,Vt,m,n) !求矩陣的轉(zhuǎn)置integer m,n,i,jreal*8:V(m,n),Vt(n,m)do i=1,m do j=1,n Vt(j,i)=V(i,j) enddoenddoend subroutine zhuanzhi!-subroutine cheng(P,Q,R,m1,n,m2) !矩陣相乘integer n,m1,m2,i,jreal*8:P(m1,n),Q(n,m

20、2),R(m1,m2)do i=1,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) enddo enddoenddoend subroutine cheng !-subroutine inv(A,invA,N) !計算逆矩陣implicit noneinteger ninteger ireal*8:A(n,n),invA(n,n),E(n,n)E=0 do i=1,n !設(shè)置E為單位矩陣 E(i,i)=1end docall mateq(A,E,invA,N,N)end subroutine invsubroutine ma

21、teq(A,B,X,N,M) !高斯列主元消去法計算矩陣方程 AX=Binteger N,M,ireal*8:A(N,N),B(N,M),X(N,M)real*8:btemp(N),xtemp(N)do i=1,M btemp=B(:,i) call elgauss(A,btemp,xtemp,N) X(:,i)=xtempend doend subroutine mateqsubroutine elgauss(A,b,x,N) !高斯列主元消去法 Ax=binteger i,k,Ninteger id_max !主元素標(biāo)號real*8:A(N,N),b(N),x(N)real*8:Aup(N

22、,N),bup(N)real*8:Ab(N,N+1) !Ab為增廣矩陣 Abreal*8:vtemp1(N+1),vtemp2(N+1)Ab(1:N,1:N)=AAb(:,N+1)=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(

23、i,k)/Ab(k,k) Ab(i,:)=Ab(i,:)-temp*Ab(k,:) end doend doAup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)call uptri(Aup,bup,x,n)end subroutine elgausssubroutine uptri(A,b,x,N) ! 上三角方程組的回帶方法Ax=binteger i,j,Nreal*8:A(N,N),b(N),x(N)x(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) end do

24、x(i)=x(i)/A(i,i)end doend subroutine uptri!-截圖： 圖4.1 系數(shù)矩陣截圖 圖4.2 平差求邊界力系數(shù)的最小二乘解運行截圖2、應(yīng)變值的尺度歸化PROGRAM GUIHUAIMPLICIT NONEinteger u0,sigma0integer nreal,allocatable:T(:),u(:),sigma(:),T1(:)integer i,j,kprint*,'請輸入標(biāo)準(zhǔn)尺度的數(shù)學(xué)期望和均方差'read*,u0,sigma0print*,'請輸入待歸化的一般正態(tài)分布的個數(shù):'read*,nif(n<2)

25、pauseallocate(T(n),u(n),sigma(n),T1(n)print*,'請依次輸入各組應(yīng)變值、期望、方差'do i=1,nread*,T(i),u(i),sigma(i)enddodo j=1,nT1(j)=(T(j)-u(j)/sigma(j)T1(j)=u0+(T(j)-u(j)*sigma0/sigma(j)enddoprint*,'原應(yīng)變值',' 歸化后的應(yīng)變值'do k=1,nprint*,T(k),T1(k)enddodeallocate(T,u,sigma,T1)end截圖： 圖4.3 應(yīng)變值的尺度歸化程序運行截圖課程設(shè)計感想一轉(zhuǎn)眼長達三周的實習(xí)就要結(jié)束了，雖然中間經(jīng)歷了頗多挫折，但臨結(jié)束卻還是有幾分不舍，或許是因為已經(jīng)大四了。言歸正傳，下面談一下這次實習(xí)的大概經(jīng)理與自己的一些感受吧。最初知道這個實習(xí)長達三周后，其實我的內(nèi)心是拒絕的，因為作為一名考研黨，在這三周里要分出大量的時間去機房，必然減少了復(fù)習(xí)的時間，再加上幾門選修課穿插其間，心里那個苦的呦！再回顧這三周的實習(xí)工作，感受很深，收獲頗豐，雖然大部分時間都是在機房度過，中間難免在教學(xué)樓、圖書館之間來回奔波，但每天都過得很充實，而且