綜合地球物理作業(yè)

1、二度體及三度體重力正演姓名：XXX 學(xué)號(hào)：XXX前言模型正演是地球物理處理中的重要手段，同時(shí)也是一些反問題構(gòu)建的基礎(chǔ)，所以正演問題在地球物理問題中具有非常關(guān)鍵的地位。對(duì)一些地球物理問題的認(rèn)識(shí)往往是從正演問題開始的。在重力勘探中，正演模擬的地位就非常重要。由于重力勘探采用的是相對(duì)測(cè)量的手段，所以正問題主要考慮的基本是地下密度異常所引起的相對(duì)重力異常。重力異常（重力加速度異常）是密度異常的一種體積效應(yīng)。如果已知異常體的形狀、位置等參數(shù)，理論上是可以解析計(jì)算出這一異常體引起異常的函數(shù)。不過實(shí)際實(shí)現(xiàn)過程中，這樣做會(huì)耗費(fèi)巨大的計(jì)算代價(jià)。通常的做法是將地下非規(guī)則的異常體近似成相近的規(guī)則異常體，求解規(guī)則異常

2、體的位置、幾何參數(shù)以及密度；或者對(duì)地下空間進(jìn)行網(wǎng)格剖分，求解每個(gè)網(wǎng)格內(nèi)的密度參數(shù)。規(guī)則異常體正演的參數(shù)空間是異常體中心位置、幾何規(guī)模（例如正方體模型的邊長(zhǎng)）、密度差。而基于網(wǎng)格剖分的模型正演的參數(shù)空間與所選擇的網(wǎng)格大小有關(guān)，如果網(wǎng)格確定，那么參數(shù)空間就是每個(gè)網(wǎng)格里密度差。根據(jù)地質(zhì)實(shí)際，地下構(gòu)造往往在沿著走向上有一定的長(zhǎng)度，在走向上構(gòu)造的變化不大，例如向斜、背斜、斷層等構(gòu)造。對(duì)于這樣的在走向上有一定延續(xù)度并且走向延續(xù)度遠(yuǎn)大于我們研究區(qū)域，這時(shí)我們稱這樣的模型為二度體，即我們只需要考慮垂直測(cè)線方向和垂直地面這兩個(gè)空間尺度。但是對(duì)于塊狀地質(zhì)異常體，例如花崗巖侵入體，這樣沒有明顯構(gòu)造走向的異常體，我們

3、必須考慮空間三個(gè)方向的異常體形態(tài)，對(duì)于這樣的異常體，我們稱為三度體。介于二者之間，我們僅考慮構(gòu)造延伸尺度的異常體，我們稱為二度半異常體。這次作業(yè)中，我主要考慮二度體和三度體的正演。對(duì)于二度體的情況，我考慮使用二度臺(tái)階異常的公式模擬一個(gè)褶皺界面所引起的重力異常。對(duì)于三度體，我主要考慮使用長(zhǎng)方體重力正演公式模擬柱狀花崗巖體所引起的重力異常。相關(guān)公式介紹二度臺(tái)階異常體公式：那么二度矩形異常公式：直立長(zhǎng)方體模型重力異常公式：其中，i,j,k=1,2。計(jì)算代碼正確性驗(yàn)證通過理論表達(dá)式計(jì)算坐標(biāo)原點(diǎn)的理論異常值和程序數(shù)值計(jì)算所得結(jié)果，可以證明所獲得的程序是正確的。二度體正演算例考慮灰?guī)r和砂巖交界的起伏界面所

4、引起的密度異常。給定砂巖密度為，為下伏層，灰?guī)r密度為，為上覆層。那么在界面處，相對(duì)水平層狀介質(zhì)有的密度異常。正演結(jié)果如下： 三度體正演算例：考慮花崗巖侵入大理巖所引起的密度異常，花崗巖密度：，大理巖密度：。結(jié)果如下：附錄：二度體：正演函數(shù)forward：function g=forward(x,T) %¶þ¶È°å×´ÌåÒì³£×Ó³ÌÐò%UNTITLED Summary of this fu

5、nction goes here% Detailed explanation goes here x0=T(1);%Öеãºá×ø±ê a=T(2);%°ë¿í¶È H=T(3);%ÉÏµ× h=T(4);%ÏÂµ× Q=T(5);%ÃܶȲî f=6.67*10(-11);%Ò

6、53;Á¦³£Êý for i=1:length(x) %×ǫ́½×Òì³£ X(i)=x(i)-(x0-a); p1(i)=pi*(H-h); p21(i)=(X(i)2+H2)/(X(i)2+h2); p2(i)=X(i)*log(p21(i); p3(i)=2*H*atan(X(i)/H); p4(i)=-2*h*atan(X(i)/h); g1(i)=f*Q*(p1(i)+p2(i)+p3(i)+p4(i);%À&#

7、237;ÂÛ¼ÆË㹫ʽ %ÓǪ̀½×Òì³£ X(i)=x(i)-(x0+a); p1(i)=pi*(H-h); p21(i)=(X(i)2+H2)/(X(i)2+h2); p2(i)=X(i)*log(p21(i); p3(i)=2*H*atan(X(i)/H); p4(i)=-2*h*atan(X(i)/h); g2(i)=f*Q*(p1(i)+p2(i)+p3(i)+p4(i);%&

8、#192;íÂÛ¼ÆË㹫ʽ g(i)=g1(i)-g2(i);%¶þ¶È°å×´ÌåÒì³££¨×ǫ́½×ÓëÓǪ̀½×Òì³£²

9、 g(i)=106*g(i);%½«ÖØÁ¦Òì³£µÄµ¥Î»×ª»¯³É¸ß˹ endend調(diào)用程序s1：clear allclose allx=0:5:400; %¹Û²âµãºá×ø±

10、êz=0; %¹Û²âµã×Ý×ø±êT=200,20,400,10,0.3;%¸ø³öÕýÑÝÄ£ÐͲÎÊýTg=forward(x,T);%µ÷ÓÃÖØÁ¦Òì³£×

11、1;³ÌÐòfigure(1)plot(x,g,'.-r') %»³öÖØÁ¦Òì³£ylabel('Deltag/g.u.','FontSize',12)xlabel('²âµãºá×ø±ê/m','FontSize',12)mo=zeros(length(x),le

12、ngth(x);for j=(200-20)/5+1:(200+20)/5 for i=3:length(x) mo(i,j)=500; endendfigure(2)pcolor(x,x,mo)axis ij調(diào)用程序s2：clear allclose allx=-2000:20:2000; %¹Û²âµãºá×ø±êz=0; %¹Û²âµã×Ý×ø±êh0

13、=2000;A=1000;mx=-1000:1:1000;y=zeros(1,length(x);for i=1:length(mx) H(i)=h0+A*sin(i-1)/1000*pi); T(1)=mx(i); T(2)=0.5; T(3)=H(i); T(4)=h0; T(5)=-0.3; g=forward(x,T); y=y+forward(x,T); %if (mod(i,10)=0) %figure(i) %plot(x,g,'.-g') %ylabel('Deltag/g.u.','FontSize',12) %xlabel(&

14、#39;²âµãºá×ø±ê/m','FontSize',12) %endendfigure(10000)plot(x,y,'.-r') %»³öÖØÁ¦Òì³£ylabel('Deltag/g.u.','FontSize',12)xlabel('²âµãº&

15、#225;×ø±ê/m','FontSize',12)figure(100000)plot(mx,H,'.-r') %»³öÉî¶È½çÃæylabel('½çÃæÉî¶È/m','FontSize',12)xlabel('²âµãº

16、25;×ø±ê/m','FontSize',12)三度體：正演程序chang：function data_c=chang(dm,x1,x2,y1,y2,z1,z2)xx=zeros(2,1);yy=zeros(2,1);zz=zeros(2,1);xx(1)=x1;xx(2)=x2;yy(1)=y1;yy(2)=y2;zz(1)=z1;zz(2)=z2;sum=0;for ii=1:2 for jj=1:2 for kk=1:2 ro=(xx(ii)2+yy(jj)2+zz(kk)2)0.5; sum=sum+(-1)(ii+jj+kk)*(xx(ii)*log(yy(jj)+ro)+yy(jj)*log(xx(ii)+ro)+zz(kk)*(pi+atan2(zz(kk)*ro,xx(ii)*yy(jj); end endenddata_c=-6.67*dm*sum/1000;end調(diào)用程序changf：clearda=20;a=0:da:5000;middle_x=2500;middle_y=2500;dm=0.5;dl=20;h0=100;data=zeros(length(a)