SEISMIC FORWARD MODELLING FOR SYNTHETIC P-P & P-S REFLECTION (REVIEWED & IMPLEMENTED BY MATLAB)

Figure.1. Exploration Seismology

On this article, I just tried to present some aspect about exploration seismology and visualized by MATLAB, especially a seismic modelling by 2D raytracing in a layered medium and AVO modelling based on synthetic data. It is very pleasant for us of being able to share, to discuss, and to improve our knowledge about geophysical method for oil and gas exploration. Hope at least this little bit will interest us to propose our kind ideas.

Seismic Forward Modelling

The most important technique for exploring oil and gas is to perform the subsurface model of earth’s interior by seismic imaging. “Imaging”, we mean the representation of the earth’s subsurface model. Due to this ultimate goal, geophysicist has imagined seismic forward modelling method by using both numerical study in geoscience and computation technology. Seismic forward modelling will assist us to understand the aspect of seismic traveltime, arrival time, elastic impedance, reflectivity, amplitude that can be generated by seismic energy and the other aspect and it also can be accomplished by ray tracing and synthetic data

.

Figure.2. (a). Geology Model of Fault and Anticline (b). Seismic Section Expressed Geology Model After Post-Stack Kirchhoff Depth Migration

Seismic Raytracing In a Layered Medium

In a layered medium, the slowness parallel to the interfaces is conserved. Let us first consider a model of homogeneous, plane layers, in the X-Z plane. If the velocity in a i-th layer with thickness ∆Z_i , offset ∆X_i is V_i, in this layer we have.

Figure.3. (a). Schematic of raytracing principle through the i-th layer. At each layer the ray is either reflected or transmitted (b). Rays propagate from sources to receivers through flat layers. Red asterisk denotes seismic source, blue box denotes receiver.

Normally we have to be satisfied a unique constant parameter (p) for any particular ray, is commonly called the ray parameter. θ_i is the angles involved in Snell’s law, can be considered as the angles between the raypaths and the normal to the interface. Up to the critical angle, the raypath intersect at the interface satisfying Snell’s law.

From the geometry of the ray segment in the layer, we can easily calculate the traveltime (T) and how far it goes horizontally to the X axis.

Anyway, Snell’s law is incorporated by replacing the trigonometric functions, so this is similar equation.

Because the medium is a discrete layer, so the summation of ∆X_i refers to the offset of ray will be traveled. In practise, we have to fix it iteratively by using shooting method until any particular ray refers to the receivers position. The summation of ∆T_i refers to total traveltime, also commonly called as two way traveltime (twt).

The offsets are examined by fan of rays and hopefully p values will be found that bracket the desired offset. An alternate approach to finding the intersection is computed iteratively by using bisection technique. If the layer is not flat, we can intersect ray fit to the layer by using spline function, polynomial derivative, interpolation, and line-search method.

Figure.4. (a). Ray diagram for wireline/horizontal survey in a layered medium(b). Seismic traveltime are computed from raytracing

Figure.5. (a). Ray diagram for borehole/vertical seismic profile in a layered medium (b). Seismic traveltime are computed from raytracing

Bisection technique can be applied for solving nonlinear equations like f (x) = 0, only in the case where we know some interval [a, b] on which f (x) is continuous and the solution uniquely exists and, most importantly, f (a) and f (b) have the opposite signs. The procedure toward the solution of f (x) = 0 is described as follows and is cast into MATLAB program.

Amplitude Versus Offset (AVO) Modelling

A variation of reflection and transmission coefficients with incidence angles, and thus offset, is commonly known as amplitude versus offset (AVO). Aki and Richards said that plane-wave reflection and transmission coefficients as a function of incidence angle and elastic properties of the medium. So, to understand AVO is done in order to determine the variations in P-wave (Vp), S-wave (Vs), and density of the medium. It means that changes in sesismic amplitudes with offset are due to contrast in the physical rock properties.

Castagna provided a relationship between Vp and Vs for water-saturated clastic silicate rocks, such as mudstone.

Gardner derived a relationship between density and Vp using empirical observations.

When a plane wave reaches a i-th interface with the particular angle, some of it’s energy is reflected by the boundary and some of its energy is transmitted through the boundary. The portioning of seismic energy is described by the plane wave amplitudes of reflected and transmitted as a function of incidence angle at a i-th interface.

Figure.6. Seismic energy proportion are generated by an incident P-wave through subsurface boundary

Several seismologist has made approximation for plane-wave amplitudes reflection and transmission coefficient. They are Zoeppritz, Shuey, and Thomsen equation. Below, I enclose the PP reflection coefficient which only in a short equation, the complete equations are written in MATLAB function programs.

Figure.7. Raytracing for P-P mode by an incident P-wave in a layered medium

Figure.8. Raytracing for P-S mode by an incident P-wave in a layered medium

Zoeppritz equations are quite complicated and thus their physical meanings are not easily seen. PP and PS reflection coefficient is given below.

Shuey re-wrote the approximation which has been introduced by Aki and Richards in below form. Shuey approximation is fully linearized and more applicable to analyze gradient AVO.

Thomsen approximation is commonly used for various computation to describre the magnitude and the type of rock. His approximation has become the most effective way to compute the portioning of seismic energy at an interface. PP reflection coefficient is given below.

Synthetic Test For P-P & P-S Reflection

For proper analyze, we also need to examine the variation of amplitude as changes of incidence angle. The AVO as seismic response generate from the synthetic data, such as elastic parameter, geology model, and spatial geometry of the source-receiver offsets configuration. Below, I tried to make a simple lithology of 13 layers, including model of P wave velocity, S wave velocity, and density. Seismic reflection geometry is designed from 0 to 2000 metres and covers 2000 metres of maximum depth. Receiver interval 20 metres, 3 shot positions at 0 m, 1000 m, and 2000 m along the surface of geometry setting.

Figure.9. Lithology model of P wave velocity, S wave velocity, and density

Lithology model, geometry experiment, source-receiver configuration, and cross plots are figured by Matlab program (model.m), raytracing process are computed by shooting.m (modified from CREWES program), Reflection coefficient Rpp and Rps are calculated by zoeppritz.m, shuey.m, and thomsen.m, synthetic log for P sonic, S sonic, and density log are estimated by logsyn.m. For study purpose, I enclose Matlab program and it is recomended to modify all source codes if there were wrong algorithms.

Figure.10. Raytracing P-P mode by an incident P-wave in a layered medium of seismic reflection experiment.

Figure.11. Raytracing P-S mode by an incident P-wave in a layered medium of seismic reflection experiment

Figure.12. Stacked section of seimogram synthetic for P-P reflection coefficient after convolving with 20 Hz ricker wavelet.

Figure.13. Stacked section of seimogram synthetic for P-S reflection coefficient after convolving with 20 Hz ricker wavelet

Figure.14. Stacked section of seimogram synthetic for both P-P and P-S reflection coefficient after convolving with 20 Hz ricker wavelet

Figure.15. Comparison of approximation between exact reflection coefficients versus absolut incidence angle from lithology model

Figure.16.(a). Cross plot of intercept vs gradient AVO from the lithology model, (b). Cross plot of reflection coefficient vs sine square of angles.

Synthetic Log & Seismogram

Figure.17 P sonic, S Sonic, and density log (synthetic) from the lithology model

Figure.18 Cross plot of physical rock properties from the lithology model

Figure.19. (a). Poisson ratio, (b). Synthetic seismogram after convolving with 20 Hz ricker wavelet, (c). Stacked section of PP seimic traces from AVO modelling

Further, the results are useful in applying seismic inversion, to get better approximation of actual geology, and also to analyze the physical rock properties in petroleum systems. The majority of oil and gas industry have initiated and still develop their seismic imaging methodology. Be able to characterize accurately the reservoirs that are associated with paradox formation or complex geology is not pure easily applied. As prior, we have to do low-cost investigation before drilling by using sophisticated tools. Usually an extensive and successful drilling program has been conducted by forwarding the best seismic imaging.

Marine seismic data acquisition (From Schlumberger, 1990)

However, either exploration or exploitation of hydrocarbon is not just a journey to employ people, upgrading tools, allocating budget to support acquiring the data, extensive drilling, and the other supporting activities. Naturally, resources and disaster or accident are illustrated as two sides of coin. Always becoming most important notice to consider the aspect of environmental and social sensitivity, safety, energy saving, reduction poverty and prosperity of life.

That’s all folks….

References :

Carcuz, J.R., 2001, A combined AVO analysis of P-P and P-S reflection data: SEG Expanded Abstracts 20,323-325.

R. James Brown and Alexandru Vant. Relative polarity of PP and PS events in the registration process and approximations to the PS reflection coefficient. CREWES Research Report — Volume 17 (2005).

MATLAB SOURCE CODES (model.m)

% model.m

%==========================================================================================

% MAIN PROGRAM

%==========================================================================================

clc; clear all; close all;

format short g

%============================================================================================

% Define Geometry

%============================================================================================

fprintf('---> Defining the interfaces ...\n');

% Input geometry

xmin = 0; xmax = 2000;

zmin = 0; zmax = 2000;

% Make strata layer

zlayer = [0;80;180;300;450;630;800;1050;1200;1400;1620;1820;2000];

nlayer = length(zlayer);

layer = 1:1:nlayer;

thick = abs(diff(zlayer));

x = [xmin xmax]; z = [zlayer zlayer];

dg = 10;

xx = xmin:dg:xmax; nx = length(xx);

zz = repmat(zlayer,1,nx);

fprintf(' Geometry has been defined...[OK]\n');

%============================================================================================

% Source-Receiver Groups

%============================================================================================

fprintf('---> Setting source-receiver configuration ...\n');

% Receiver Interval

dr = 20;

% Mix Configuration

% Source

% xs = [55 1350];

% zs = [0 690];

% ns = length(xs);

%

% % Receiver

% xr = [200 500 850 1600 1900];

% zr = [300 0 220 150 0];

% nr = length(xr);

% --------------------------------------

% Front Configuration

% Source

% xs = xmin;

% zs = 0;

% ns = length(xs);

%

% % Receiver

% xr = xmin+dr:dr:xmax;

% zr = zeros(length(xr));

% nr = length(xr);

% --------------------------------------

% End Configuration

% Source

% xs = xmax;

% zs = 0;

% ns = length(xs);

%

% % Receiver

% xr = xmax-dr:-dr:xmin;

% zr = zeros(length(xr));

% nr = length(xr);

% --------------------------------------

% Split Spread Configuration

% Source

% xs = fix((xmax-xmin)/2);

% zs = 0;

% ns = length(xs);

%

% % Receiver

% xr = [fix((xmax-xmin)/2)-dr:-dr:xmin fix((xmax-xmin)/2)+dr:dr:xmax];

% zr = zeros(length(xr));

% nr = length(xr);

% --------------------------------------

% Off-End Configuration

% Source

% xs = [xmin xmax];

% zs = [0 0];

% ns = length(xs);

%

% % Receiver

% xr = xmin+dr:dr:xmax-dr;

% zr = zeros(length(xr));

% nr = length(xr);

% Front Spread Configuration

% Source

xs = [xmin fix((xmax-xmin)/2) xmax];

zs = [0 0 0];

ns = length(xs);

% Receiver

xr = [xmin+dr:dr:fix((xmax-xmin)/2)-dr fix((xmax-xmin)/2)+dr:dr:xmax-dr];

zr = zeros(length(xr));

nr = length(xr);

% VSP Configuration

% Source

% xs = 1950;

% zs = 0;

% ns = length(xs);

%

% % Receiver

% zr = 0:100:1600;

% xr = 50.*ones(length(zr));

% nr = length(xr);

nray = ns*nr;

fprintf(' Source-Receiver Groups have been setted...[OK]\n');

%============================================================================================

% Create Elastic Parameter

%============================================================================================

fprintf('---> Creating elastic parameter ...\n');

% Create synthetic Vp, Vs, and Density

vlayer = [1800;2000;2200;2500;2400;2700;3150;2950;3440;3750;4000;4350;4600]; % Velocity

vlayer = vlayer(1:length(zlayer));

vel = [vlayer vlayer];

vp = vlayer; % P wave velocity

vs = (vp-1360)./1.16; % S wave velocity based on Castagna's rule

ro = 0.31.*(vp).^0.25; % Density based on Gardner's rule

pois = (vs.^2-0.5*vp.^2)./(vs.^2-vp.^2); % Poisson ratio

fprintf(' No.Layer Depth(m) Vp(m/s) Vs(m/s) Density(kg/m^3) Poisson Ratio \n');

disp([layer' zlayer vp vs ro pois]);

fprintf(' Elastic parameters have been created...[OK]\n');

figure;

set(gcf,'color','white');

% Plot P wave velocity

subplot(1,3,1);

stairs(vp,zlayer,'LineWidth',2.5,'Color',[0.07843 0.1686 0.549]);

ylabel('Depth (m)','FontWeight','bold','Color','black');

title('V_{p} (m/s)','FontWeight','bold');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','top','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle','-','FontWeight',...

'demi','FontAngle','italic');

% Plot S wave velocity

subplot(1,3,2);

stairs(vs,zlayer,'LineWidth',2.5,'Color',[1 0 0]);

title('V_{s} (m/s)','FontWeight','bold');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','top','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle','-','FontWeight',...

'demi','FontAngle','italic');

% Plot Density

subplot(1,3,3);

stairs(ro,zlayer,'LineWidth',2.5,'Color',[0.07059 0.6392 0.07059]);

title('Density (kg/m^3)','FontWeight','bold');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','top','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle','-','FontWeight',...

'demi','FontAngle','italic');

% Plot Geology Model

figure;

set(gcf,'color','white');

color = load('mycmap.txt');

pcolor(xx,zz,repmat(vlayer,1,nx)); shading flat; hold on

colormap(color); colorbar ('horz');

axis([0 xmax zmin-0.03*zmax zmax])

set(gca,'YDir','reverse','XaxisLocation','bottom',....

'Ytick',zlayer,'FontWeight','demi','PlotBoxAspectRatioMode','Manual',...

'PlotBoxAspectRatio',[2.4 1.2 1],'Position',[0.04 0.30 0.90 0.60]);

% Plot Source-Receiver Group

plot(xs,zs,'r*','markersize',12); hold on

plot(xr,zr,'sk','markersize',4,'markerfacecolor','c'); hold on

xlabel('Velocity (m/s)','FontWeight','bold','Color','black');

ylabel('Depth (m)','FontWeight','bold','Color','black');

%============================================================================================

% Run Ray Tracing

%============================================================================================

fprintf('---> Starting ray tracing ...\n');

wat = waitbar(0,'Raytracing is being processed, please wait...');

%============================================================================================

xoff = [];

% Loop over for number of souce

tic

for i=1:ns

% Loop over for number of receiver

for j=1:nr

%====================================================================================

% Loop over for number of layer

for k=1:nlayer

%=====================================================================================

% Declare reflection boundary

if and(zr(j) < zlayer(k),zs(i) < zlayer(k))

zm = zz(k,:); zf = min(zm)- 100000*eps;

% Downgoing path

d = find(zlayer > zs(i));

if(isempty(d)); sdown = length(zlayer); else sdown = d(1)-1; end

d = find(zlayer > zf);

if(isempty(d)); edown = length(zlayer); else edown = d(1)-1; end

zd = [zs(i);zlayer(sdown+1:edown)]; nd = length(zd);

% Upgoing path

u = find(zlayer > zr(j));

if(isempty(u)); sup = length(zlayer); else sup = u(1)-1; end

u = find(zlayer > zf);

if(isempty(u)); eup = length(zlayer); else eup = u(1)-1; end

zu = [zr(j);zlayer(sup+1:eup+1)]; nu = length(zu);

zn = [zd;(flipud(zu))]; nrefl = length(zn)-1;

%=================================================================================

% Declare elastic parameter

% Downgoing elastic parameter

vpd = [vp(sdown:edown);vp(edown)];

vsd = [vs(sdown:edown);vs(edown)];

rod = [ro(sdown:edown);ro(edown)];

% Upgoing elastic parameter

vpu = [vp(sup:eup);vp(eup)];

vsu = [vs(sup:eup);vs(eup)];

rou = [ro(sup:eup);ro(eup)];

% Combine model elastic parameter

vpp = [vpd(1:end-1);flipud(vpu(1:end-1))];

vss = [vsd(1:end-1);flipud(vsu(1:end-1))];

vps = [vpd(1:end-1);flipud(vsu(1:end-1))];

rho = [rod(1:end-1);flipud(rou(1:end-1))];

%=================================================================================

% Start Raytracing (P-P, S-S, or P-S mode)

ops = 1; % ops=1 for PP mode; ops=2 for PS mode

[xh,zh,vh,pp,teta,time] = shooting(vpp,vps,zn,xx,xs(i),xr(j),ops);

theta = abs(teta); twt(k,j,i) = time;

% Plot Ray

if ops == 1

plot(xh,zh,'k-');

title('Seismic Raytracing (P-P mode)','FontWeight','bold');

elseif ops == 2

xd = xh(1:nd+1); xu = xh(nd+1:end);

zd = zh(1:nd+1); zu = zh(nd+1:end);

plot(xd,zd,'k-',xu,zu,'r-');

title('Seismic Raytracing (P-S mode)','FontWeight','bold');

end

%=================================================================================

% Compute Reflection Coefficient (Downgoing-Upgoing)

for c=1:nrefl-1

% Reflection Coefficient of Zoeppritz Approximation

[rc1,teta] = zoeppritz(pp,theta(c),vpp(c),vss(c),rho(c),theta(c+1),...

vpp(c+1),vss(c+1),rho(c+1),ops);

rcz(c,j,i) = rc1;

% Reflection Coefficient of Shuey Approximation

[A,B,rc2,teta] = shuey(theta(c),vpp(c),vss(c),rho(c),theta(c+1),...

vpp(c+1),vss(c+1),rho(c+1),ops);

rcs(c,j,i) = rc2; AA(c,j,i) = A; BB(c,j,i) = B;

% Reflection Coefficient of Thomsen Approximation

[rc3,teta] = thomsen(theta(c),vpp(c),vss(c),rho(c),theta(c+1),...

vpp(c+1),vss(c+1),rho(c+1));

rct(c,j,i) = rc3;

angle(c,j,i) = teta.*(180/pi);

end

%====================================================================================

end

end % for horizon/reflector end

xoff = [xoff xr(j)];

waitbar(j/nray,wat)

end % for receiver end

% Save Data

save(['time_shot',num2str(i),'.mat'],'twt');

save(['reflz_shot',num2str(i),'.mat'],'rcz');

save(['refls_shot',num2str(i),'.mat'],'rcs');

save(['reflt_shot',num2str(i),'.mat'],'rct');

save(['teta_shot',num2str(i),'.mat'],'angle');

save(['intercept_shot',num2str(i),'.mat'],'AA');

save(['gradient_shot',num2str(i),'.mat'],'BB');

waitbar(i/ns,wat)

end % for sources end

%=============================================================================================

close(wat);

toc

fprintf(' Ray tracing has succesfully finished...[OK]\n');

xx = repmat(xr,nlayer,1);

% Plot Traveltime

figure;

for i=1:ns

tm = load(['time_shot',num2str(i),'.mat']); tt = tm.twt;

plot(xx(2:nlayer,:),tt(2:nlayer,:,i),'b.'); hold on

xlabel('Horizontal Position (m)','FontWeight','bold','FontAngle','normal','Color','black');

ylabel('Time (s)','FontWeight','bold','FontAngle','normal','Color','black');

title('Traveltime','FontWeight','bold');

set(gca,'Ydir','reverse');

set(gcf,'color','white');

end

% Plot Reflection Coefficient

figure;

for i=1:ns

reflz = load(['reflz_shot',num2str(i),'.mat']); reflz = reflz.rcz;

refls = load(['refls_shot',num2str(i),'.mat']); refls = refls.rcs;

reflt = load(['reflt_shot',num2str(i),'.mat']); reflt = reflt.rct;

theta = load(['teta_shot',num2str(i),'.mat']); angle = theta.angle;

%Reflection Coefficient of Zoeppritz

subplot(1,3,1)

plot(abs(angle(1:nlayer,:,i)),reflz(1:nlayer,:,i),'r.'); hold on

xlabel('Incidence Angles (degree)','FontWeight','bold','Color','black');

ylabel('Reflection Coefficient','FontWeight','bold','Color','black');

grid on; title('Rpp Zoeppritz','FontWeight','bold','Color','black')

set(gca,'YColor',[0.04314 0.5176 0.7804],'XColor',[0.04314 0.5176 0.7804]); hold on

%Reflection Coefficient of Shuey

subplot(1,3,2)

plot(abs(angle(1:nlayer,:,i)),refls(1:nlayer,:,i),'g.'); hold on

xlabel('Incidence Angles (degree)','FontWeight','bold','Color','black');

grid on; title('Rpp Shuey','FontWeight','bold','Color','black')

set(gca,'YColor',[0.04314 0.5176 0.7804],'XColor',[0.04314 0.5176 0.7804]); hold on

%Reflection Coefficient of Thomsen

subplot(1,3,3)

plot(abs(angle(1:nlayer,:,i)),reflt(1:nlayer,:,i),'b.'); hold on

xlabel('Incidence Angles (degree)','FontWeight','bold','Color','black');

grid on; title('Rpp Thomsen','FontWeight','bold','Color','black')

set(gca,'YColor',[0.04314 0.5176 0.7804],'XColor',[0.04314 0.5176 0.7804]); hold on

set(gcf,'color','white');

end

for i=1:ns

refls = load(['refls_shot',num2str(i),'.mat']); refls = refls.rcs;

theta = load(['teta_shot',num2str(i),'.mat']); angle = theta.angle;

Rt = load(['intercept_shot',num2str(i),'.mat']); Rt = Rt.AA;

Gt = load(['gradient_shot',num2str(i),'.mat']); Gt = Gt.BB;

Ro(1:nlayer,:,i) = Rt(1:nlayer,:,i);

Go(1:nlayer,:,i) = Gt(1:nlayer,:,i);

Rc(1:nlayer,:,i) = refls(1:nlayer,:,i);

inc(1:nlayer,:,i) = angle(1:nlayer,:,i);

end

Rp = reshape(Ro,nr*ns*nlayer,1);

G = reshape(Go,nr*ns*nlayer,1);

R = reshape(Rc,nr*ns*nlayer,1);

teta = reshape(inc,nr*ns*nlayer,1);

% Plot Attributes

figure;

% Rp-G Cross Plot

subplot(1,2,1)

plot(Rp,G,'r.'); hold on

plot(Rl,G,'m-'); hold on

xlabel('R','FontWeight','bold','Color','black');

ylabel('G','FontWeight','bold','Color','black');

grid on; title('R-G Cross Plot','FontWeight','bold','Color','black')

set(gca,'YColor',[0.04314 0.5176 0.7804],'XColor',[0.04314 0.5176 0.7804]); hold on

% R-sin^2(teta) Cross Plot

subplot(1,2,2)

plot((sin(teta).^2),R,'g.'); hold on

xlabel('sin^2(teta)','FontWeight','bold','Color','black');

ylabel('R(teta)','FontWeight','bold','Color','black');

grid on; title('R-sin^2(teta) Cross Plot','FontWeight','bold','Color','black')

set(gca,'YColor',[0.04314 0.5176 0.7804],'XColor',[0.04314 0.5176 0.7804]); hold on

set(gcf,'color','white');

%============================================================================================

% AVO Modelling

%============================================================================================

fprintf('---> Starting AVO Modelling ...\n');

% Make wavelet ricker (f = 2 Hz)

dt = 0.004; f = 20;

[w,tw] = ricker(dt,f);

% Create zeros matrix for spike's location

tmax = max(tt(:));

tr = 0:dt:tmax; nt = length(tr);

t = tt(2:nlayer,:);

% Take reflectivity into spike location

spikes = zeros(nt,nray);

rcz = real(reflz(2:nlayer,:));

for k=1:nlayer-1

for j=1:nray

ir(k,j) = fix(t(k,j)/dt+0.1)+1;

spikes(ir(k,j),j) = spikes(ir(k,j),j) + rcz(k,j);

end

end

% Convolve spikes with wavelet ricker

for j=1:nray

seisz(:,j) = convz(spikes(:,j),w);

end

ampz=reshape(seisz,length(tr),nr,ns);

for i=1:ns

ampz_shot(:,:,i) = ampz(:,:,i);

save(['ampz_shot',num2str(i),'.mat'],'ampz_shot');

end

% CMP Stacked

Stacked = zeros(nt,nr);

for i=1:ns

files = load(['ampz_shot',num2str(i),'.mat']);

files = files.ampz_shot;

Stacked = Stacked + files(:,:,i);

end

% Save Seismic Data

save(['seis','.mat'],'Stacked');

% Plot AVO

figure;

wig(xr,tr,Stacked,'black'); hold on

xlabel('Offset (m)','FontWeight','bold','Color','black');

ylabel('Time (s)','FontWeight','bold','Color','black');

title('Synthetic Seismogram','FontWeight','bold','Color','black');

axis tight

set(gca,'YColor',[0.04314 0.5176 0.7804],'XColor',[0.04314 0.5176 0.7804]); hold on

set(gcf,'color','white');

fprintf(' AVO Modelling has been computed...[OK]\n');

%============================================================================================

% Seismogram Synthetic

%============================================================================================

% Create Smooth Synthetic Log From Elastic Data Set

[pson,sson,spson,rlog,zlog] = logsyn(vp,vs,ro,zlayer,500);

vplog = 1.0e6./pson; vslog = 1.0e6./sson;

rholog = 0.31.*(vplog).^0.25; vpvs = vplog./vslog;

poisson = (vslog.^2-0.5*vplog.^2)./(vslog.^2-vplog.^2);

% Plot Synthetic Data Log

figure;

set(gcf,'color','white');

% Plot P Sonic Log

subplot(1,3,1);

stairs(pson,zlog,'LineWidth',1,'Color',[0.07843 0.1686 0.549]);

xlabel('P sonic (US/m)','FontWeight','bold','Color','black');

ylabel('Depth (m)','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','XDir','reverse','YColor',...[0.04314 0.5176 0.7804],'XAxisLocation','top','XColor',...

[0.04314 0.5176 0.7804],'MinorGridLineStyle','-','FontWeight','demi','FontAngle','italic');

% Plot S Sonic Log

subplot(1,3,2);

stairs(sson,zlog,'LineWidth',1,'Color',[1 0 0]);

xlabel('S sonic (US/m)','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','XDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','top','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle','-','FontWeight',...

'demi','FontAngle','italic');

% Plot Density Log

subplot(1,3,3);

stairs(rlog,zlog,'LineWidth',1,'Color',[0.07059 0.6392 0.07059]);

xlabel('Density (kg/m^3)','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','top','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle','-','FontWeight',...

'demi','FontAngle','italic');

% Plot Rock Properties

figure;

% Cross Plot P and S wave velocity

subplot(2,2,1);

plot(vplog,vslog,'.','Color',[0.07843 0.1686 0.549]);

xlabel('Vp (m/s)','FontWeight','bold','Color','black');

ylabel('Vs (m/s)','FontWeight','bold','Color','black');

title('Cross Plot of Vp and Vs','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle',...

'-','FontWeight','demi','FontAngle','italic');

% Plot S Sonic Log

subplot(2,2,2);

plot(vplog,vpvs,'.','Color',[1 0 0]);

xlabel('Vp (m/s)','FontWeight','bold','Color','black');

ylabel('Vp/Vs ratio','FontWeight','bold','Color','black');

title('Cross Plot of Vp and Vp/Vs Ratio','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],...

'MinorGridLineStyle','','FontWeight',...

'demi','FontAngle','italic');

% Plot Density Log

subplot(2,2,3);

plot(rholog,vpvs,'.','Color',[0.07059 0.6392 0.07059]);

xlabel('Density (kg/m^3)','FontWeight','bold','Color','black');

ylabel('Vp/Vs ratio','FontWeight','bold','Color','black');

title('Cross Plot of Vp and Vp/Vs Ratio','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle',...

'-','FontWeight','demi','FontAngle','italic');

subplot(2,2,4);

plot(poisson,vpvs,'m.');

xlabel('Poisson Ratio','FontWeight','bold','Color','black');

ylabel('Vp/Vs ratio','FontWeight','bold','Color','black');

title('Cross Plot of Vp and Vp/Vs Ratio','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle',...

'-','FontWeight','demi','FontAngle','italic');

set(gcf,'color','white');

% Integrate Depth and Velocity To Get Two-Way Time

tstart = 0; tracenum = 10;

[tz,zt] = sonic2tz(spson,zlog,-100,tstart);

tzobj = [zt tz];

vins = vplog; dens = rholog; z = zlog;

% Determine Ricker Wavelet

dt = 0.004; f = 20;

[w,tw] = ricker(dt,f);

[theo,tlog,rcs,pm,p] = seismogram(vins,dens,z',w,tw,tzobj,0,0,1);

% Convolves Reflectivity With Wavelet

nint = length(tlog);

RC = zeros(length(tlog),1);

for i=1:nint

ir(i,:) = fix(tlog(i,:)/dt+0.1)+1;

RC(ir(i,:),:)= RC(ir(i,:),:) + rcs(i,:);

end

synth = convz(RC,w);

% Compare AVO Data Estimation

rseis = load(['seis','.mat'],'Stacked'); num = 5:15;

rseis = rseis.Stacked; seis = rseis(:,num);

% Converts surface seismic data from time to depth

[seist,zp,dzp,tz,zt]=seis2z(synth,tlog,vp,zlayer,vs,1,min(zlayer),max(zlayer));

[seisr,zp,dzp,tz,zt]=seis2z(seis,tr,vp,zlayer,vs,1,min(zlayer),max(zlayer));

figure;

subplot(131);

stairs(poisson,zlog,'LineWidth',1,'Color',[0.07843 0.1686 0.549]);

xlabel('Poisson Ratio','FontWeight','bold','Color','black');

ylabel('Depth (m)','FontWeight','bold','Color','black');

title('Poisson Ratio','FontWeight','bold','Color','black');

set(gca,'XMinorGrid','on','YMinorGrid','on','YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle',...

'-','FontWeight','demi','FontAngle','italic');

subplot(132);

wig(tracenum,zp',seist,'red',0.0002); hold on

xlabel('Trace Number','FontWeight','bold','Color','black');

title('Synthetic Seismogram','FontWeight','bold','Color','black');

set(gca,'YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle',...

'-','FontWeight','demi','FontAngle','italic');

subplot(133);

wig(num,zp',seisr,'blue'); hold on

xlabel('Trace Number','FontWeight','bold','Color','black');

title('AVO Estimation','FontWeight','bold','Color','black');

axis tight

set(gca,'YDir','reverse','YColor',[0.04314 0.5176 0.7804],...

'XAxisLocation','bottom','XColor',[0.04314 0.5176 0.7804],'MinorGridLineStyle',...

'-','FontWeight','demi','FontAngle','italic');

set(gcf,'color','white');

SUBFUNCTION

shooting.m

function [xh,zh,vh,pp,teta,time] = shooting(vpp,vps,zn,xx,xs,xr,ops)

%============================================================================================

% Some Constants

%============================================================================================

itermax = 20;

offset = abs(xs-xr); xc = 10;

%============================================================================================

% Determine Option

%============================================================================================

if ops == 1;

vh = vpp;

elseif ops == 2

vh = vps;

end

% Initial guess of the depth & time

zh = zn - 100000*eps;

t = inf*ones(size(offset));

p = inf*ones(size(offset));

%============================================================================================

% Start Raytracing

%============================================================================================

% Trial shooting

pmax = 1/min(vh);

pp = linspace(0,1/max(vh),length(xx));

sln = vh(1:length(zh)-1)*pp;

vel = vh(1:length(zh)-1)*ones(1,length(pp));

dz = abs(diff(zh))*ones(1,length(pp));

if(size(sln,1)>1)

xn = sum((dz.*sln)./sqrt(1-sln.^2));

tt = sum(dz./(vel.*sqrt(1-sln.^2)));

else

xn = (dz.*sln)./sqrt(1-sln.^2);

tt = dz./(vel.*sqrt(1-sln.^2));

end

xmax = max(xn);

%============================================================================================

% Bisection Method

%============================================================================================

% Start Bisection Method

for k=1:length(offset)

% Analyze the radius of target

n = length(xn);

xa = xn(1:n-1); xb = xn(2:n);

opt1 = xa <= offset(k) & xb > offset(k);

opt2 = xa >= offset(k) & xb < offset(k);

ind = find((opt1) | (opt2));

if (isempty(ind))

if (offset(k) >=xmax);

a = n; b = [];

else

a = []; b = 1;

end

else

a = ind; b = ind + 1;

end

x1 = xn(a); x2 = xn(b);

t1 = tt(a); t2 = tt(b);

p1 = pp(a); p2 = pp(b);

iter = 0; err= (b-a)/2;

% Minimize the error & intersect the reflector

while and(iter < itermax,abs(err) < 1)

iter = iter + 1;

xt1 = abs(offset(k) - x1); xt2 = abs(offset(k) - x2);

if and(xt1 < xc,xt1 <= xt2)

% Linear interpolation

t(k) = t1 + (offset(k)-x1)*(t2-t1)/(x2-x1);

p(k) = p1 + (offset(k)-x1)*(p2-p1)/(x2-x1);

elseif and(xt2 < xc,xt2<=xt1)

% Linear interpolation

t(k) = t2 + (offset(k)-x2)*(t1-t2)/(x1-x2);

p(k) = p2 + (offset(k)-x2)*(p1-p2)/(x1-x2);

end

% Set new ray parameter

if (isempty(a));

p2 = p1; p1 = 0;

elseif (isempty(b))

p1 = p2; p2 = pmax;

end

pnew = linspace(min([p1 p2]),max([p1 p2]),3);

% Do shooting by new ray parameter

sln = vh(1:length(zh)-1)*pnew(2);

vel = vh(1:length(zh)-1)*ones(1,length(pnew(2)));

dz = abs(diff(zh))*ones(1,length(pnew(2)));

if (size(sln,1)>1)

xtemp = sum((dz.*sln)./sqrt(1-sln.^2));

ttemp = sum(dz./(vel.*sqrt(1-sln.^2)));

else

xtemp = (dz.*sln)./sqrt(1-sln.^2);

ttemp = dz./(vel.*sqrt(1-sln.^2));

end

xnew = [x1 xtemp x2]; tnew = [t1 ttemp t2]; xmax = max(xnew);

% Analyze the radius of target

n = length(xnew);

xa = xnew(1:n-1);xb = xnew(2:n);

opt1 = xa <= offset(k) & xb > offset(k);

opt2 = xa >= offset(k) & xb < offset(k);

ind = find((opt1) | (opt2));

a = ind; b = ind + 1;

x1 = xnew(a); x2 = xnew(b);

t1 = tnew(a); t2 = tnew(b);

p1 = pnew(a); p2 = pnew(b);

err = (b - a)/2;

% Declare ray parameter

if xr > xs; pp = p; else pp = -p; end

% Compute travel time & angle

dx = real((pp.*vh.*dz)./sqrt(1-pp.*pp.*vh.*vh));

xx = xs + cumsum(dx); xh = [xs;xx];

dz = real(dx.*sqrt(1-pp.*pp.*vh.*vh))./(pp.*vh);

dt = dz./(vh.*sqrt(1-pp.*pp.*vh.*vh));

tt = cumsum(dt); time = tt(end);

teta = real(asin(pp*vh));

end % End the while loop

end % End the offset loop

% zoeppritz.m

function [RC,incidence] = zoeppritz(pp,teta1,vp1,vs1,rho1,teta2,vp2,vs2,rho2,ops)

costeta1 = cos(teta1);

costeta2 = cos(teta2);

j1 = asin(pp*vs1);

j2 = asin(pp*vs2);

cosj1 = cos(j1);

cosj2 = cos(j2);

a = rho2*(1-2*vs2^2*pp^2)-rho1*(1-2*vs1^2*pp^2);

b = rho2*(1-2*vs2^2*pp^2)+2*rho1*vs1^2*pp^2;

c = rho1*(1-2*vs1^2*pp^2)+2*rho2*vs2^2*pp^2;

d = 2*(rho2*vs2^2-rho1*vs1^2);

E = b*costeta1/vp1+c*costeta2/vp2;

F = b*cosj1/vs1+c*cosj2/vs2;

G = a-d*(costeta1/vp1)*(cosj2/vs2);

H = a-d*(costeta2/vp2)*(cosj1/vs1);

D = E*F+G*H*pp^2;

if ops==1

Rpp = ((b*costeta1/vp1-c*costeta2/vp2)*F-(a+d*(costeta1/vp1)*(cosj2/vs2))*H*pp^2)/D;

RC = Rpp;

elseif ops==2

Rps = (-2*(costeta1/vp1)*(a*b+c*d*(costeta2/vp2)*(cosj2/vs2)*(pp*vp1)))/(vs1*D);

RC = Rps;

end

incidence = teta1;

% shuey.m

function [A,B,RC,incidence] = shuey(teta1,Vp1,Vs1,rho1,teta2,Vp2,Vs2,rho2,ops)

i=(teta1+teta2)/2;

Z1 = Vp1*rho1;

Z2 = Vp2*rho2;

Z = (Z1+Z2)/2;

dZ = Z2-Z1;

Vp = (Vp1+Vp2)/2;

dVp = Vp2-Vp1;

Vs = (Vs1+Vs2)/2;

dVs = Vs2-Vs1;

Rb = dVs/(2*Vs);

gam = (Vs1+Vs2)/(Vp1+Vp2);

rho = (rho1+rho2)/2;

Rro = (rho2-rho1)/(2*rho);

po1 = (((Vp1/Vs1)^2)-2)/(2*(((Vp1/Vs1)^2)-1));

po2 = (((Vp2/Vs2)^2)-2)/(2*(((Vp2/Vs2)^2)-1));

po = (po1+po2)/2;

dpo = po2-po1;

A = (1/2)*(dZ/Z);

B = -2*(1-(po/(1-po)))*A-(1/2*((1-3*po)/(1-po))*(dVp/Vp))+(dpo/(1-po)^2);

C = 1/2*(dVp/Vp);

if ops==1

Rpp = A+B*(sin(i))^2+C*(sin(i))^2*(tan(i))^2;

RC = Rpp;

elseif ops==2

Rps = -gam*((Rro+2*gam(2*Rb+Rro))*sin(i));

RC = Rps;

end

incidence = teta1;

% thomsen.m

function [Rpp,incidence] = thomsen(teta1,Vp1,Vs1,rho1,teta2,Vp2,Vs2,rho2)

i = (teta1+teta2)/2;

Z1 = Vp1*rho1;

Z2 = Vp2*rho2;

Z = (Z1+Z2)/2;

dZ = Z2-Z1;

Vp = (Vp1+Vp2)/2;

dVp= Vp2-Vp1;

Vs = (Vs1+Vs2)/2;

G1 = rho1*Vs1^2;

G2 = rho2*Vs2^2;

G = (G1+G2)/2;

dG = G2-G1;

Rpp=((1/2)*(dZ/Z))+1/2*((dVp/Vp)-((2*Vs/Vp)^2*(dG/G)))*((sin(i))^2)+((1/2)*(dVp/Vp)*((sin(i))^2)*((tan(i))^2));

incidence = teta1;

% logsyn.m

function [pson,sson,spson,rlog,zlog] = logsyn(vp,vs,ro,z,n)

curvep = fit(z,vp,'pchipinterp'); csp = curvep.p;

curves = fit(z,vs,'pchipinterp'); css = curves.p;

curver = fit(z,ro,'pchipinterp'); csr = curver.p;

if n==0,

pson = 1.0e6./vp;

sson = 1.0e6./vs;

spson = (pson + sson)/2;

rlog = 1.0e6./ro;

zlog = z;

else

zlog = linspace(min(z),max(z),n);

plog = ppval(csp,zlog);

slog = ppval(css,zlog);

rlog = ppval(csr,zlog);

pson = 1.0e6./plog - (10)*rand(1,length(zlog));

pson = pson(:);

sson = 1.0e6./slog - (10)*rand(1,length(zlog));

sson = sson(:);

spson = (pson + sson)/2;

rlog = rlog - 0.1*rand(1,length(zlog));

rlog = rlog(:);

end;

% seis2z.m

function [seisz,z,dz,tz,zt]=seis2z(seis,t,vp,zv,vs,ops,zmin,zmax)

% converts PP or PS surface seismic data from time to depth.

%

% seis : PP or PS data in time.

% t : time axis of PP or PS data.

% vp : P-wave velocity log.

% zv : depth axis for velocity log.

% vs : S-wave velocity log.

% ops : conversion type, ops=1 for PP data and ops=2 for PS data

% zmin : minimum desired depth.

% zmax : maximum desired depth.

% seisz : the output PP or PS data.

% z : the output depth-axis.

if ops==1; % For PP data

vs = vp;

end;

v = vs; % Get the same sample for both PP and PS data.

dt = t(2)-t(1); % Time interval

[nt,nx] = size(seis); % nt is the number of time steps.

pson=(1.e06./vp); % P-sonic

sson=(1.e06./vs); % S-sonic

spson=(pson + sson)/2; % Fake PS-sonic

% Computes an approximate 2-way time-depth curve from a sonic log for use with depth conversion.

tstart=0;

[tz,zt] = sonic2tz(spson,zv,-10,tstart); % Get time-depth curve.

dz = min(v)*dt/2; % Sample depth

nz = floor((zmax-zmin)/dz)+1; % Number of depth sample

z = (0:nz-1)*dz + 0; % Depth axis

z = z(:); % Vector for depth

t2 = interp1(zt,tz,z);

nt2 = length(t2);

seisz = zeros(nt2,nx); % Data in depth, nx = number of offsets

for k=1:nx,

seisz( : ,k) = sinci(seis(: ,k),t',t2); % Interpolate the amplitude at each depths sample

end

% wig.m

function wig(x,t,Data,style,dmax,showmax,plImage,imageax)

cla;

showmax_def = 500;

style_def = 'wiggle';

if nargin==1,

Data = x;

t = 1:1:size(Data,1);

x = 1:1:size(Data,2);

dmax=max(Data(:));

style=style_def;

showmax=showmax_def;

end

if nargin==2,

Data = x; dmax = t;

t = 1:1:size(Data,1);

x = 1:1:size(Data,2);

style=style_def;

showmax=showmax_def;

end

if nargin==3,

style = style_def;

dmax = max(abs(Data(:)));

showmax = showmax_def;

end

if nargin==4,

dmax=max(abs(Data(:)));

showmax=showmax_def;

end

if nargin==5,

showmax=showmax_def;

end

if nargin<7

plImage=0;

end

if isempty(dmax),

dmax=max(abs(Data(:)));

end

if isempty(showmax),

showmax=100;

end

if nargin==7,

imageax=[-1 1].*dmax;

end

if plImage==1,

imagesc(x,t,Data);

if (length(imageax)==1)

imageax=[-1 1].*abs(imageax);

end

caxis(imageax);

hold on

end

if (showmax>=0)

if length(x)>1

dx = x(2)-x(1);

else

dx = max(Data)*15;

end

ntraces = length(x);

d = ntraces/showmax;

if d<=1; d=1; end

d=round(d);

dmax = dmax/d;

for i=1:d:length(x)

xt=dx*Data(:,i)'./dmax;

if (strmatch('black',style)==1)

xt1 = xt; xt1(find(xt1>0))=0;

f1=fill(x(i)+[xt,fliplr(xt1)],[t,fliplr(t)],[0 0 0]); hold on

set(f1,'LineWidth',0.0001)

set(f1,'EdgeColor',[0 0 0])

plot(xt+x(i),t,'k-','linewidth',.05); hold on

elseif (strmatch('red',style)==1)

xt2 = xt; xt2(find(xt2>0))=0;

f2=fill(x(i)+[xt,fliplr(xt2)],[t,fliplr(t)],[1 0 0]); hold on

set(f2,'LineWidth',0.0001)

set(f2,'EdgeColor',[1 0 0])

plot(xt+x(i),t,'r-','linewidth',.05); hold on

elseif (strmatch('green',style)==1)

xt3 = xt; xt3(find(xt3>0))=0;

f3=fill(x(i)+[xt,fliplr(xt3)],[t,fliplr(t)],[0 1 0]); hold on

set(f3,'LineWidth',0.0001)

set(f3,'EdgeColor',[0 1 0])

plot(xt+x(i),t,'g-','linewidth',.05); hold on

elseif (strmatch('blue',style)==1)

xt4 = xt; xt4(find(xt4>0))=0;

f4=fill(x(i)+[xt,fliplr(xt4)],[t,fliplr(t)],[0 0 1]); hold on

set(f4,'LineWidth',0.0001)

set(f4,'EdgeColor',[0 0 1])

plot(xt+x(i),t,'b-','linewidth',.05);hold on

else

plot(xt+x(i),t,'k-','linewidth',.05); hold on

end

if i==1, hold on;end

end

end

if length(x)>1

axis([min(x)-1.2*x(1) max(x)+1.2*x(1) min(t) max(t)])

else

axis([min(x)-1.2*min(x) max(x)+1.2*max(x) min(t) max(t)])

end

set(gca,'Ydir','reverse')

Posted in: