matlab实习作业

             Matlab语言及应用结课报告

             题目：利用matlab语言建立速度扫描程序

         学院：地空学院

             年级：061114

             姓名：张宁  学号：20111000486

             指导老师： 王贤敏

2015年6月25日

                    目  录

第一章 模型背景与意义.................................1

第二章 建立模型的思路与步骤...........................2

第三章 模型演示及结果................................ 4   

第四章 结论与评价....................................14

第五章 参考文献......................................15

附录.................................................16

         第一章  模型背景与意义

              1.1 MATLAB产生的历史背景   

   MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。 

  MATLAB是矩阵实验室（Matrix Laboratory）的简称，和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 

  MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且mathwork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。 

MATLAB的应用：利用MATLAB可以进行以下各种工作   

①数值计算   

②数值和符号计算   

③工程和科学绘图   

④控制系统的设计与仿真   

⑤数字图像处理   

⑥数字信号处理   

⑦通讯系统设计与仿真   

⑧财务与金融工程 

1．2 MATLAB的语言特点和开发环境

   MATLAB作为一种科学计算的高级语言之所以受欢迎，就是因为它有丰富的函数资源和工具箱资源，编程人员可以根据自己的需要选择函数，而无需再去编写大量繁琐的程序代码，从而减轻了编程人员的工作负担。被称为编程语言的MATLAB最大的特点就是简洁开放的程序代码和直观实用的开发环境。具体地说MATLAB主要有以下特点：   

●库函数资源丰富  

●语言精炼，代码灵活   

●运算符多而灵活 

●面向对象控制功能优良，使界面编程更方便、自由。  

●程序设计自由 

●图形功能强大  

●程序的兼容性好  

●源代码开放  

●形形色色的工具箱 

第二章.建立模型的思路与步骤

              2.1.建立模型的思路

地震勘探利用地下介质弹性和密度的差异，通过观测和分析人工地震产生的地震波在地下的传播规律，推断地下岩层的性质和形态的地球物理勘探方法叫做地震勘探。地震勘探是钻探前勘测石油与天然气资源的重要手段，在煤田和工程地质勘查、区域地质研究和地壳研究等方面，也得到广泛应用。

    地震波的传播速度这一参数参与到了地震数据的采集处理与解释全部三个过程中，是地震勘探过程中最重要的参数之一。而地震勘探的最终目的则是获得与地震有关的物性参数，因此正确的提取参数是各项处理工作得以顺利进行的前提和保证。

当地区的地震地质条件复杂、得不到高质量的速度谱时，便可以应用第四种速度分析方法，即速度扫描。速度扫描是最直观、最简单、最可靠的速度分析方法，因为他是直接从动校正记录或叠加道上提取速度。

利用MATLAB编程，进行速度扫描这一方法的数值实验，一方面熟悉了MATLAB语言的使用，学习并掌握了程序的编写和调试过程。另一方面加深了对本专业课程的认识，有助于巩固并深透的理解所学的理论知识，并练习了将理论知识化为实际应用的方法，是非常有意义的一次学习过程。

    速度扫描方法本质上还是动校正的应用，动校正是将共中心点道集中炮检距不同的各道上来自同一界面同一点的反射波到达时ti校正为共中心点处回声时间的一种处理方法，动校正的目的是使各道反射波到达时相同，从而在叠加时实现同相叠加，得到反射波能量突出的叠加道。

正常时差公式：

           Δti=ti-t0=sqrt（t0^2+xi^2/v^2)-t0

在上式中，t0=2h/v为共中心点处回声时间，即自激自收时间，xi为炮检距。

     本次速度扫描的程序就是基于动校正原理进行编写。

               2.2 建立模型的步骤

   Ⅰ.先进行问题分析，明确问题

     本文的问题是实现地震勘探中速度扫描这一速度分析方法，并能成功利用编写出的程序对原始的地震记录进行速度扫描并得出成果。

     Ⅱ.建立模型，并利用matlab函数求解

     本文的matlab程序包含三个部分，第一部分是给各个变量赋值，并导入地震记录数据；第二部分是计算过程，包含了一个三重循环语句，原理是动校正公式，即对上文给出的不同的V值进行动校正然后得出结果；第三部分是画图，这里调用了wigb画图函数，wigb是根据地震数据画出二维剖面图的一个地震处理专用的作图程序。

     Ⅲ.对结果进行分析说明

     Ⅳ.评价并提出改进意见     

        第三章 模型演示及结果

    本章模型设计是利用的Tesseral 2-D软件进行的，Tesseral 2-D是第一个基于PC的商业                  化的全波场模拟软包，用它可以非常容易地建立复杂的地质模型剖面，并且模拟不同的地震观测系统。

    Tesseral 2-D – 采用人机交互的方式建立、检查速度-深度模型，并与常规地震数据处理和解释直接集成。系统中包含一个典型的地质环境数据库，包括P波和S波的速度及密度和岩性符号等。这样可以很方便地输入岩石物理特性来构建地质模型。同时，也可以 把不同网格距和不同文本格式的外部模型输入到该系统进行模拟计算。

    所建立的地质模型一共有三层介质，第一层是白色部分，第二层是绿色部分，第三层是紫色部分，第一层介质设定的层速度（P波速度）为2000m/s, 第二层介质设定的层速度（P波速度）为2500m/s，第三层介质设定的层速度（P波速度）为3000m/s。

这个地质模型包含了一个垂直台阶和一个正断层，垂直台阶部分的深度是2000m和2600m，正断层断面之上与水平方向成20°夹角，断面之下与水平方向成40°夹角。

     如下图所示是建立的地震地质模型

正演模拟得到了102张的共炮点原始地震记录。以下分别是第1、21、41、61、81、101炮的原始地震记录。

                            图3-2 共炮点地震记录

                              图3-3 共炮点地震记录

                                   图3-4 共炮点地震记录

                                  图3-5 共炮点地震记录

                                 图3-6 共炮点地震记录

                               图3-7  共炮点地震记录

利用segy软件将原始的炮集记录转换为segy格式，再选取道间隔、起始道、终止道等数据，为了便于分析反射波同相轴，也可以通过选取采样点进行切割。如图所示为抽道集方法图。

以上4张CDP道集记录都是基于48次覆盖的观测系统抽取的，分别从第3炮、第7炮、第10炮、第25炮的第105道记录开始抽取，道间隔是102道。

扫描速度是依据地质模型设定的层速度选择的，从1900m/s到3100m/s，100m/s的递增，一共13个扫描速度。

为了方便同相轴的识别与比较，将每个扫描速度动校正之后的道集记录压缩并合并。所得结果与分析如后面所示。

从图中可以很直观的看到第一层反射波同相轴的校正情况，第一层在选取2000m/s的扫描速度下校正情况最好，同相轴比较平直，在选取2000m/s以上的扫描速度的情况下，反射波同相轴上凸程度递增，这表示校正不足，而在选取1900m/s的扫描速度时，反射波同相轴下凹，校正过量了。这张图只得出了第一层的最佳动校正速度，而在下面那两层的反射波能量很弱，基本无法分辨出同相轴的变化。

图中显示的的第一层反射波的动校正情况与图4-9是一致的，但是图4-10已经能够模糊的辨认出第三层的反射波同相轴的校正情况，约为2900-3100m/s，而第二层的反射波同相轴仅仅很模糊的辨认。但是因为能量比较弱，情况还是不明显，但是比图4-9要好，我认为这是因为图4-10选择的道集中第二层与第三层的反射界面较图4-9中的要浅一些。从这一里我们也可以得出二层与三层的速度分界面是倾斜的，而且是从左到右向上倾斜。

图中对于这三层的反射波同相轴已经能够很好的辨认了，分析同相轴的校正情况也可以得出第二层和第三层大致的动校正速度范围。

                      图4-12 速度扫描成果图

在这个图上，第二层的反射波同相轴要比第一层清晰，我认为是波在向下传播时碰到了断层改变了频谱和能量分布特征，表现为了能量的会聚。

            第四章 结论与评价

从成果图可以看到做出的matlab程序是比较成功的，但是在成图效果上较差，主要问题是地震波的振幅过大，并且经过压缩后十二张图合成一张使得振幅太大导致图片很黑。

    Matlab是一个强大的数学工具,它的应用广泛,涉及到各个领域.它使用起来十分方便,不用麻烦去定义变量.它的绘图能力很强,甚至可以模拟出三维视图.矩阵是它应用的核心,许多工程繁琐的运算都需要靠矩阵来化简,这正是它的生命力所在.但是,他的函数很多,开始学时记的比较痛苦,我已经深深感觉到了.它的语法简单,它的数组定义十分符合自然,是从1开始的,数组元素的调用也很接近数学的表达.此外,函数的名字也很符合英文规则。

    通过学习matlab,我又一次锻炼了自己的思维.也让我明白了学习一门语言对学习其他方面的帮助指导作用.同时,它也加强了我理论联系实际的能力.这是一个专业课的基础工具,学好它是必要的.

                第五章 参考文献

[1].杨淑卿,陈云峰等.STseis叠前成像系统偏移速度分析技术[J].油气地球物

    理,2011,01(06):11-14.

[2].张恒超.几种简单实用的多次波去除方法及其应用[J].内蒙古石油化 

    工,2011,21(03):21-24.

[3].刘超颖,王成祥.CFP层速度扫描建模方法[J].石油物探,2003,03(04):294-297.

[4].潘祖福,吴彦伟等.单次相遇时距曲线折射静校正方法[J].天然气工

    业,1992,05(04):104-105.

[5].俞寿朋.地震勘探资料的数字处理[J].石油地球物理勘探,1977,01(04):65-71.

[6].王建华.一种求取地层层速度与空间归位的新方法[J].石油物探,1978,04(01):38-52.

[7].肖奼莉,李保田.提高速度分析精度的方法[J].石油地球物理勘

    探,1983,02(04):154-166.

[8].陈天胜,刘洋.纵波和转换波联合AVO反演方法研究[J].中国石油大学学

    报,2006,01(03):33-37.

[9].潘宏勋,方伍宝.地震速度分析方法综述[J].勘探地球物理进展,2006,05(07):305-311.

[10].陈雨红,李强.分偏移距动校正方法的改进[J].石油物

     探,2007,02(07):107-111+125+13.

[14].王辉,丁志峰.浅层地震勘探资料处理中的速度分析参数选取[J].地震地

     质,2006,04(02):597-603.

[11].孙祥娥,詹毅等.双参数速度扫描方法研究[J].石油天然气学

     报,2006,06(05):65-68+5+4.

[12].包吉山,明.二维倾斜叠加速度分析[J].石油地球物理勘

     探,1988,01(04):34-54+130.

[13].Hans J.Tieman,张伯君.非扫描叠加速度分析法[J].石油物探译

     丛,1994,02(05):1-14.

[14].王辉,丁志峰.浅层地震勘探资料处理中的速度分析参数选取[J].地震地

     质,2006,04(02):597-603.

[15].魏修成,刘洋.高精度转换波速度分析[J].石油勘探与开发,2000,02(04):57-59+10-3.

[16].宋建国.井间地震偏移速度扫描和反射波超级叠加技术[J].石油大学学

     报,1999,01(06):36-39+6.

[17].王立新.自动拾取叠加速度及提取偏移速度方法[J].石油物探,1999,04(03):107-113.

[18].王棣,匡斌,孟祥宾等.三维叠前深度偏移的CRP道集速度扫描分析技术及应用[J].石

     油地球物理勘探,2002,06(02):615-619.

[19].张凯,李振春,曹景忠.速度扫描法共聚焦点道集偏移速度分析[J].同济大学学

     报,2005,09(02):1265-1269.

[20].陈云峰,杨淑卿,王华忠.STseis地震成像处理系统偏移速度分析技术[J].海洋地质前

     沿,2012,08(04):-69.

附录Ⅰ：

load ap;

A0 = ap;

[w,l]=size(A0);

A = A0';

dett=0.004;

x = 1:1:53;

z = 0:0.004:2482*0.004;

wigb(A,1,x,z);

v = 1900:100:3100;

lv = length(v);

[l,lx] = size(A);

delta_t = zeros(lx,lv); 

t0 = 2.05;

for i0 = 1:1:lv

    for i1 = 1:1:lx

        delta_t(i1,i0) = ( abs( x(i1) - 27 )*2*40 )^2 / ( 2 * t0 * v(i0)^2 );

    end

end

%NMO

A1 = zeros(l,lx,lv); 

for i2 = 1:1:lv

    for i3 = 1:1:lx

        m0 = ceil( delta_t(i3,i2)/0.004 );   %取整

        n1 = m0;

        n2 = l-n1;

        for i4 =  n1:1:n2-1

            A1(i4-n1+1,i3,i2) = A(i4+1,i3);

        end

    end

end

figure

set(gcf,'outerposition',get(0,'screensize'));

for i5 = 1:lv

    subplot(1,lv,i5)

    wigb(A1(:,:,i5),1,x,z)

    str = num2str(v(i5));

    xlabel(str)

End

附录Ⅱ：

if nargin == 0, nx=10;nz=10; a = rand(nz,nx)-0.5; end;

[nz,nx]=size(a);

trmx= max(abs(a));

if (nargin <= 4); amx=mean(trmx); end;

if (nargin <= 2); x=[1:nx]; z=[1:nz]; end;

if (nargin <= 1); scal =1; end;

if nx <= 1; disp(' ERR:PlotWig: nx has to be more than 1');return;end;

 % take the average as dx

    dx1 = abs(x(2:nx)-x(1:nx-1));

    dx = median(dx1);

 dz=z(2)-z(1);

 xmx=max(max(a)); xmn=min(min(a)); 

 if scal == 0; scal=1; end;

 a = a * dx /amx; 

 a = a * scal;

 fprintf(' PlotWig: data range [%f, %f], plotted max %f \\n',xmn,xmx,amx);

% set display range 

x1=min(x)-2.0*dx; x2=max(x)+2.0*dx;

z1=min(z)-dz; z2=max(z)+dz;

set(gca,'NextPlot','add','Box','on', ...

  'XLim', [x1 x2], ...

  'YDir','reverse', ...

  'YLim',[z1 z2]);

    fillcolor = [0 0 0];

    linecolor = [0 0 0];

    linewidth = 0.1;

    z=z';   % input as row vector

    zstart=z(1);

    zend  =z(nz);

for i=1:nx,

  if trmx(i) ~= 0;    % skip the zero traces

    tr=a(:,i);  % --- one scale for all section

    s = sign(tr) ;

    i1= find( s(1:nz-1) ~= s(2:nz) );   % zero crossing points

    npos = length(i1);

    %12/7/97 

    zadd = i1 + tr(i1) ./ (tr(i1) - tr(i1+1)); %locations with 0 amplitudes

    aadd = zeros(size(zadd));

[zpos,vpos] = find(tr >0);

    [zz,iz] = sort([zpos; zadd]);   % indices of zero point plus positives

    aa = [tr(zpos); aadd];

    aa = aa(iz);

    % be careful at the ends

        if tr(1)>0, a0=0; z0=1.00;

        else,       a0=0; z0=zadd(1);

        end;

        if tr(nz)>0, a1=0; z1=nz;

        else,       a1=0; z1=max(zadd);

        end;

    zz = [z0; zz; z1; z0];

    aa = [a0; aa; a1; a0];

    zzz = zstart + zz*dz -dz;

    patch( aa+x(i) , zzz,  fillcolor);

    line( 'Color',[1 1 1],'EraseMode','background',  ...

         'Xdata', x(i)+[0 0], 'Ydata',[zstart zend]); % remove zero line

%'LineWidth',linewidth, ...

%12/7/97 'Xdata', x(i)+[0 0], 'Ydata',[z0 z1]*dz);  % remove zero line

    line( 'Color',linecolor,'EraseMode','background',  ...

     'LineWidth',linewidth, ...

     'Xdata', tr+x(i), 'Ydata',z);  % negatives line

   else % zeros trace

    line( 'Color',linecolor,'EraseMode','background',  ...

     'LineWidth',linewidth, ...

         'Xdata', [x(i) x(i)], 'Ydata',[zstart zend]);

   end;

end;