基于Matlab實(shí)現(xiàn)的地震波場(chǎng)邊界處理軟件

1、基于Matlab實(shí)現(xiàn)的地震波場(chǎng)邊界處理軟件姓名：姚嘉德 學(xué)號(hào)：2015301130007院系：資源與環(huán)境科學(xué)學(xué)院摘要：用有限差分法模擬地震波場(chǎng)是研究地震波在地球介質(zhì)中傳播的有效方法。但我們?cè)趯?shí)驗(yàn)室進(jìn) 行波場(chǎng)數(shù)值模擬時(shí)有限差分網(wǎng)格是限制在人工邊界里面，即引入了人工邊界條件。本文采用Clayton_Engquist_Majda二階吸收邊界條件，通過(guò)MATLAB編程實(shí)現(xiàn)了這一算法。依靠MATLAB具有 更加直觀的、符合大眾思維習(xí)慣的代碼，為用戶提供了友好、簡(jiǎn)潔的程序開發(fā)環(huán)境，方便同行們交流。利用Matlab本身所具有可視化功能以及像素識(shí)別功能，可以將生成的動(dòng)畫電影進(jìn)行識(shí)別，用于地震局實(shí)時(shí)分析有著深遠(yuǎn)

2、意義。關(guān)鍵詞：有限差分法，地震波場(chǎng)，吸收邊界條件，MATLAB矢量幀,像素識(shí)別Abstract： Modeling seismic wave field with the Finite Difference Method (FDM) is an effective method to study theseismic wave propagation in the earth medium. When we model seismic wave field in the laboratory, the finitedifference grids are restricted in the a

3、rtificial boundary. So it should introduce the artificial boundary conditions.This paper adopts Clayton_Engquist_Majda second absorbing boundary conditions and realizes the arithmetic with MATLAB. The MATLAB codes are direct and accord with our thinking custom. So it can provide the friendlyand succ

4、inct programming environment and is easy to communicate with other.Using the functions of Matlab that make visualization come true and identify the pixel,we can identify the earthquake wave field.Key words: finite difference method, seismic wave field, numerical modeling, absorbing boundary conditio

5、ns,MATLAB一、引言用有限差分法模擬地震波場(chǎng)是研究地震波在地球介質(zhì)中傳播的有效方法 。但我們?cè)趯?shí)驗(yàn)室進(jìn) 行波場(chǎng)數(shù)值模擬時(shí)，只能在有限的空間進(jìn)行，所以有限差分網(wǎng)格是限制在人工邊界里面，即引入了 人為的邊界條件。這種人為邊界條件的引入將對(duì)有限區(qū)域內(nèi)的波場(chǎng)值的計(jì)算帶來(lái)嚴(yán)重影響，所以必 須進(jìn)行特殊的邊界處理。邊界條件處理的好壞直接影響地震正演模擬的最終效果。本文中我們采用Clayton_Engquist_Majda二階吸收邊界條件2。被稱作是第四代計(jì)算機(jī)語(yǔ)言的MATLAB語(yǔ)言，利用其豐富的函數(shù)資源把編程工作者從繁瑣的程 序代碼中解放出來(lái)。MATLAB用更加直觀的、符合大眾思維習(xí)慣的代碼，為用戶提

6、供了友好、簡(jiǎn)潔的 程序開發(fā)環(huán)境。本文介紹運(yùn)用MATLAB實(shí)現(xiàn)帶有吸收邊界條件的地震波場(chǎng)數(shù)值模擬方法和步驟，便于 同行們交流，亦可用于本科地震理論的教學(xué)中，讓學(xué)生們?cè)诔绦蜓菔局欣斫獾卣鸩ǖ膫鞑ヒ?guī)律。二.、Clayton_Engquist_Majda二階吸收邊界條件我們給定二維標(biāo)量聲波波動(dòng)方程（含震源）：（1）式中：是聲波波場(chǎng)，是聲波速度， 是震源。 對(duì)（1）式進(jìn)行時(shí)間和空間2階精度有限差分離散（見圖1），整理后可得（2）式中， 為別為空間、時(shí)間離散步長(zhǎng)， ， ，為震源函數(shù)。 震源函數(shù)：（3）Clayton_Engquist_Majda二階吸收邊界條件的微分表達(dá)式可參見文獻(xiàn)2，其左、右、上、下 邊

7、界的差分格式分別為：三、基本算法步驟從圖1可以看出，k+1時(shí)刻的波場(chǎng)值由k時(shí)刻和k-1時(shí)刻的波場(chǎng)值決定。所以在MATLAB里實(shí)現(xiàn)的基本算法步驟如下：(1) 初始時(shí)刻的全波場(chǎng)值均為零，P(i, j, dt)=0（在MATLAB中初始從dt開始，不能從0開始）；(2) 時(shí)刻2dt時(shí)，在炮點(diǎn)S (m, n)給出一個(gè)脈沖震源Src(t)（見式（3），其它點(diǎn)波場(chǎng)P =0；(3) 時(shí)刻t大于或等于3dt時(shí)，P(i, j, k+1)根據(jù)式（2）計(jì)算，其它點(diǎn)波場(chǎng)P=0；(4) 在波傳播到四周邊界時(shí)，左、右、上和下邊界的波場(chǎng)值分別由式（4-1）、（4-2）、（4-3）和（4-4）計(jì)算出來(lái)。四、數(shù)值模擬由于是計(jì)算

8、機(jī)模擬，為了能說(shuō)明問(wèn)題且便于計(jì)算，我們?cè)O(shè)地質(zhì)模型邊界為1，具體詳細(xì)參數(shù)如下見表1：V22 第一步：速度文件的載入及相關(guān)整理（替換）clc; clear; %清除工作空間及顯示屏幕load vm_0.mat; % 載入速度文件，里面包含v(j, i)Nx=101; Nz=101; Nt=800;hx=0.01;hz=0.01;dt=0.002; % 模擬參數(shù)見表1for i=1:Nxfor j=1:Nzr(j,i)=v(j,i)*dt/hx;r2(j,i)=(v(j,i)*dt/hx)2;s(j,i)=2-4*(v(j,i)*dt/hx)2;endend% 簡(jiǎn)縮“常量”u=zeros(Nz,Nx

9、,Nt); % 分布空間，預(yù)值充零第二步：賦初值for k=1:2for j=1:Nzfor i=1:Nxu(j,i,k)=0;endendend第三步：邊界條件處理及7點(diǎn)差分計(jì)算波場(chǎng)延拓for k=2:Nt-1for j=1:Nzfor i=2:Nx-1if j=1u(j,i,k+1)=(2-2*r(j,i)-r2(j,i)*u(j,i,k)+2*r(j,i)*(1+r(j,i)*u(j+1,i,k)-r2(j,i)*u(j+2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j+1,i,k-1);% 上邊界吸收elseif j=Nzu(j,i,k+1)=(2

10、-2*r(j,i)-r2(j,i)*u(j,i,k)+2*r(j,i)*(1+r(j,i)*u(j-1,i,k)-r2(j,i)*u(j-2,i,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j-1,i,k-1);% 下邊界吸收elseif j=30&i=51&(k-1)*dt <= 6*pi/100 %炮點(diǎn)S(30,51)，子波持續(xù)時(shí)間條件u(j,i,k+1)=c(j,i)2*dt2*sin(50*(k-1)*dt)*exp(-188*(k-1)*dt-3*pi/100)2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u

11、(j+1,i,k)+u(j-1,i,k)+s(j,i)*u(j,i,k)-u(j,i,k-1);elseu(j,i,k+1)=r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k)+s(j,i)*u(j,i,k)-u(j,i,k-1);endendend；for i=1:Nx4.1 地質(zhì)模型的構(gòu)造本文選取了兩個(gè)較為簡(jiǎn)單的地質(zhì)模型，分別是均勻介質(zhì)模型（見圖2）和層狀均勻介質(zhì)模型（見圖3）。4.2 程序及相關(guān)說(shuō)明根據(jù)上述我們建立的地質(zhì)模型，生成相應(yīng)的速度文件，再聯(lián)合表 1中的模擬參數(shù)和吸收邊界條件，就可以編程實(shí)現(xiàn)波場(chǎng)模擬。平時(shí)表示波場(chǎng)的習(xí)慣u(x,

12、z,t)對(duì)應(yīng)在matlab程序中，為了便于成圖則被表示為u(z,x,t)，即u(i,j,k)變?yōu)閡(j,i,k)。詳細(xì)過(guò)程如下： % 波場(chǎng)快照?qǐng)D顯示for k=1:Nt % 表示所有時(shí)刻，也可以是等間隔時(shí)間如 k=1:10:Ntfigure(k)imagesc(u(:,:,k) ; colormap(gray); % 同理可用surf(u(:,:,k)xlabel('x'); ylabel('z');% zlabel('Amplitude') 在surf(u(:,:,k)中用到title('Wave Field Snapshot'

13、);axis squareend% 二維電影動(dòng)畫放映，除了imagesc,還有contour、surf等繪圖命令clf; shg;M=moviein(Nt); %預(yù)設(shè)畫面矩陣for k=1:Ntimagesc(u(:,:,k)colormap(gray)xlabel('x');ylabel('z');title('Wave Field Snapshot');axis(0 Nx 0 Nz);M(:,k)=getframe; % 捕獲畫面end% movie(M,5,10)%電影回放，每秒10幀，重復(fù)播放5次xin_xi_ti_shi=('*

14、該程序已經(jīng)運(yùn)行結(jié)束*')4.3 運(yùn)行結(jié)果顯示運(yùn)行以上程序，可以得到所有時(shí)刻的波場(chǎng)快照?qǐng)D，在這里只選擇部分圖件進(jìn)行展示。for j=2:Nz-1if i=1u(j,i,k+1)=(2-2*r(j,i)-r2(j,i)*u(j,i,k)+2*r(j,i)*(1+r(j,i)*u(j,i+1,k)-r2(j,i)*u(j,i+2,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j,i+1,k-1); % 左邊界吸收elseif i=Nxu(j,i,k+1)=(2-2*r(j,i)-r2(j,i)*u(j,i,k)+2*r(j,i)*(1+r(j,i)*u(j,i-

15、1,k)-r2(j,i)*u(j,i-2,k)+(2*r(j,i)-1)*u(j,i,k-1)-2*r(j,i)*u(j,i-1,k-1); % 右邊界吸收elseif j=30&i=51&(k-1)*dt <= 6*pi/100 % 同上u(j,i,k+1)=c(j,i)2*dt2*sin(50*(k-1)*dt)*exp(-188*(k-1)*dt-3*pi/100)2)+r2(j,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k)+s(j,i)*u(j,i,k)-u(j,i,k-1);elseu(j,i,k+1)=r2(j

16、,i)*(u(j,i+1,k)+u(j,i-1,k)+u(j+1,i,k)+u(j-1,i,k)+s(j,i)*u(j,i,k)-u(j,i,k-1);endendendend第四步：四個(gè)角點(diǎn)的處理for k=1:Ntu(1,1,k)=1/2*(u(1,2,k)+u(2,1,k);u(1,Nx,k)=1/2*(u(1,Nx-1,k)+u(2,Nx,k);u(Nz,1,k)=1/2*(u(Nz-1,1,k)+u(Nz,2,k);u(Nz,Nx,k)=1/2*(u(Nz,Nx-1,k)+u(Nz-1,Nx,k);end第五步：保存結(jié)果及成圖u % 在主窗口顯示數(shù)值結(jié)果，按照時(shí)間切片save U u

17、 ; % 保存數(shù)值結(jié)果 圖4是同一速度均勻介質(zhì)模型的波場(chǎng)快照?qǐng)D。a、b、c和d是分別200ms、400ms、600ms和800ms時(shí)刻的波場(chǎng)快照，震源位于模型中心，波速為V1。從四個(gè)小圖上可以看出在波逐步向外傳播的過(guò)程中，當(dāng)波傳播到邊界時(shí)逐漸被吸收從而沒有人為的反射波。圖5是二層均勻介質(zhì)模型的波場(chǎng)快照?qǐng)D。a、b、c和d是分別200ms、400ms、600ms和800ms時(shí)刻的波場(chǎng)快照，震源位于模型上中部，波速分別為 V1和V2。從四個(gè)小圖上可以看出在波逐步向外傳播的過(guò)程中，當(dāng)波傳播到地層邊界時(shí)有很強(qiáng)的反射，而當(dāng)波傳播到人工邊界時(shí)逐漸被吸收沒有人為的反射波。五、成像總結(jié)本文側(cè)重于通過(guò)數(shù)值模擬介紹

18、運(yùn)用MATLAB實(shí)現(xiàn)帶吸收邊界條件的波場(chǎng)模擬的簡(jiǎn)要方法和基本步驟，因而文中數(shù)值模型和其計(jì)算差分節(jié)點(diǎn)設(shè)置都僅考慮十分簡(jiǎn)單的情形。在實(shí)際應(yīng)用中，除計(jì)算節(jié)點(diǎn)數(shù)增多外，人們不斷通過(guò)引入差分格式，差分精度等一系列措施改進(jìn)正演模擬，從而提高計(jì)算效率，減少計(jì)算誤差。以上數(shù)值模擬程序在MATLAB7.0中運(yùn)行通過(guò)，可供同行們交流。對(duì)于二層均勻介質(zhì)模型，運(yùn)行的結(jié)果以P(z, x)形式保存成dat文件，然后用地震成圖軟六、Matlab的動(dòng)畫實(shí)現(xiàn)原理MATLAB中，創(chuàng)建電影動(dòng)畫的過(guò)程分為以下四步：step1：調(diào)用moviein函數(shù)對(duì)內(nèi)存進(jìn)行初始化(該步驟在Matlab5.3以上均可省略)，創(chuàng)建一個(gè)足夠大的矩陣，使之能夠容納基于當(dāng)前坐標(biāo)軸大小的一系列指定的圖形（此處稱為幀）。step2：調(diào)用getframe函數(shù)生成每個(gè)幀。該函數(shù)返回一個(gè)列矢量，利用這個(gè)矢量，就可以創(chuàng)建一個(gè)電影動(dòng)畫矩陣。getframe函數(shù)可以捕捉動(dòng)畫幀，并保存到矩陣中。一般將該函數(shù)放到for循環(huán)中得到一系列的動(dòng)畫幀。該函數(shù)格式有：(1)F=gefframe，從當(dāng)前圖形框中得到動(dòng)畫幀(2)F=gefframe(h)，從圖形句柄h中得到動(dòng)畫幀(3)F=getframe(h，rect)，從圖形句柄h的指定區(qū)域rec中得到動(dòng)畫幀step3：調(diào)用mo