Matlab编程实现流体替换前后合成地震记录对比.docVIP

Matlab编程实现流体替换前后合成地震记录对比 我采用三层模型对流体替换进行试算： 所用模型用FORTRAN编程建模如下图所示： 所用参数具体如下： 第一层深度为500m,纵波速度为2000m/s，横波速度为1200m/s，密度2.4g/cm3。 第一层深度为650m,纵波速度为2350m/s，横波速度为1200m/s，密度2.1g/cm3。 第一层深度为110m,纵波速度为2400m/s，横波速度为1400m/s，密度2.2g/cm3。 下面对第二层进行流体替换，并做简要分析：由Gassmann方程求解，易得到给定孔隙度情况下，饱和含气或含水之后的纵横波速度值和密度值，由计算得到的速度和密度进行褶积得到合成地震记录。 以孔隙度为10%和30%为例进行计算，可得如下结果； φ=10%时，流体替换后，含水情况，Vp=3224m/s,Vs=1172,ρ=2.2g/cm3. φ=10%时，流体替换后，含气情况，Vp=2356m/s,Vs=1196.4,ρ=2.114g/cm3. φ=30%时，流体替换后，含水情况，Vp=2724m/s,Vs=1123,ρ=2.4g/cm3. φ=30%时，流体替换后，含气情况，Vp=2333m/s,Vs=1188,ρ=2.142g/cm3. 由以上数据可进行编程合成流体替换前后的地震记录。其实现程序如下： %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %流体替换前后合成地震记录对比程序 depth1=500; depth2=650; depth3=1100; velocity1=2000; velocity2=2724; velocity3=2400; velocity4=2100; sample_inter=1; trace_num=49; trace_length=1500; density1=2.4; density2=2.4; density3=2.2; density4=2.0; wave_t=160; fp=45; t1=2*depth1/velocity1; t2=t1+2*(depth2-depth1)/velocity2; t3=t2+2*(depth3-depth2)/velocity3; ref3=(density4*velocity4-density3*velocity3)/(density4*velocity4+density3*velocity3); ref2=(density3*velocity3-density2*velocity2)/(density3*velocity3+density2*velocity2); ref1=(density2*velocity2-density1*velocity1)/(density2*velocity2+density1*velocity1); sample_num=trace_length/sample_inter; for i=1:trace_num for j=1:sample_num R(j,i)=0; end k1=int32(t1*1000/sample_inter); k2=int32(t2*1000/sample_inter); k3=int32(t3*1000/sample_inter); R(k1,i)=ref1; R(k2,i)=ref2; R(k3,i)=ref3; end w_t=-wave_t/(2*1000):0.001:wave_t/(2*1000); rick=(1-2*(pi*fp*w_t).^2).*exp(-(pi*fp*w_t).^2); for i=1:trace_num syn(:,i)=conv(R(:,i),rick); end wave_l=wave_t/sample_inter; syn_l=wave_l+sample_num-1; for i=1:trace_num for j=(wave_l/2):1:(syn_l-wave_l/2) syn_n((j+1-wave_l/2),i)=syn(j,i); end end figure(2); wigb(syn_n); %title(替换前的原始地震记录) title(第二层流体替换--含水（孔隙度30%）) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 作图如下： 图1 原始的未替换前的地震记录 由上图可知，未进行流体替换之前一二层界面的反射振幅能量很弱，第三层的能量较强。 图2 第二层流体替换后（含水，φ=10%） 图3 第二层流体替换后（含气，φ=10%） 由上两图对比可知，当对第二