Code covered by the BSD License  

Highlights from
Large Data in MATLAB: A Seismic Data Processing Case Study

image thumbnail

Large Data in MATLAB: A Seismic Data Processing Case Study

by

 

These are the files used in the webinar on Feb. 23, 2011.

saltModelMigrationRTM.m
%% Seismic Migration Example
%% Data source
% 
dataLocation = 'C:\DistCompDemos\LargeDataSeismic\benchmark';
addpath(dataLocation)
addpath('fileReader')
addpath('saltToothModelData')
addpath('migration')
addpath('gpu')

%% Read in velocity model data
% SegYFileReader open the SEG Y file and reads the header information.
% Creates an object used to access trace data.
tic
V = SegYFileReader('vel_z6.25m_x12.5m_exact.segy',true,false);
toc

%% Read in all traces to memory
tic
V = V(:,:);
toc

imagesc(V), axis equal tight
colormap(seismic)

%% Velocity model parameters
nz = 1911;          % depth samples
dz = 6.25;          % depth spacing (m)
nx = 5395;          % surface samples
dx = 12.5;          % surface spacing (m)

z = (0:nz-1)*dz;
x = (0:nx-1)*dx;

% Modify colormap to show ocean
cmap = colormap(seismic);
cmap(1,:) = [0 0 0];
colormap(cmap)

% update plot
imagesc(x,z,V), axis equal tight
colorbar('NorthOutside')
xlabel('Suface Distance (m)')
ylabel('Depth (m)')

%% Pull out "tooth" segment for processing
idx = 2200:3400;
x = x(idx);
nx = length(x);
V = V(:,idx);
imagesc(x,z,V), axis equal tight
colormap seismic
drawnow

%% Resample
idx = 1:10:1200; 
x = x(idx);
nx = length(x);
idz = 1:16:1911;
z = z(idz);
nz = length(z);
V = V(idz,idx);
imagesc(x,z,V), axis equal tight
colormap seismic
drawnow

%% Calculate parameters and plot
[nz,nx] = size(V);
dx = diff(x(1:2));
dz = diff(z(1:2));

subplot(2,2,1)
imagesc(x,z,V)
xlabel('Distance (m)'); ylabel('Depth (m)');
title('Velocity Model');
hold on
hshot = plot(x(1),z(1),'w*');
hold off
colormap(seismic)

%% Create shot gathers
% Use the velocity model to simulate a seismic survey.  The wave equations
% is solved using finite differences for a defined initial wavefield.

% calculate time step dt from stability crierion for finite difference
% solution of the wave equation.
dt = 0.9*min(min(dz./V/sqrt(2)));

% determine time samples nt from wave travelime to depth and back to
% surface
vmin = min(V(:));
nt = round((sqrt((dx*nx)^2 + (dz*nx)^2)*2/vmin/dt + 1));
t  = (0:nt-1).*dt;

% add region around model for applying absorbing boundary conditions (20
% nodes wide)
V = [repmat(V(:,1),1,20) V repmat(V(:,end),1,20)];
V(end+1:end+20,:) = repmat(V(end,:),20,1);

% Define frequency parameter for ricker wavelet
f  = 20;

%% generate shots and save to file and video
%vidObj = VideoWriter('saltToothModelShots.avi');
%open(vidObj);
data = zeros(size(nt,nx));
figure(gcf)
for ixs = 21:nx+20 % shot loop
    % initial wavefield
    rw = ricker(f,nz+40,dt,dt*ixs,0);
    rw = rw(1:nz+20,:);
    
    % plot initial wavefield
    set(hshot,'XData',x(ixs-20),'YData',z(1));
    subplot(2,2,2)
    imagesc(x,z,rw(1:end-20,21:end-20))
    xlabel('Distance (m)'); ylabel('Depth (m)');
    title(['Shot ',num2str(ixs-20),' at ',num2str(x(ixs-20)),' m']);
    colormap(seismic)
    caxis([-0.14 0.2])
    
    % generate shot record
    tic
    [data snapshot] = fm2d_gpu(V,rw,nz,dz,nx,dx,nt,dt);
    tgpu=toc;
    tic
    [data snapshot] = fm2d(V,rw,nz,dz,nx,dx,nt,dt);
    tcpu=toc;
    fprintf('Single Threaded CPU: \t%12.4f\n',tcpu);
    fprintf('GPU CUDA Kernel \t\t%12.4f\n',tgpu);
    fprintf('Speedup CPU/GPU: \t\t%12.4f\n',tcpu/tgpu);
    %save(['saltToothModelData\snapshot',num2str(ixs-20),'.mat'],'snapshot');
    %save(['saltToothModelData\shotfdm',num2str(ixs-20),'.mat'],'data')
    
    data = data(21:end-20,:)';
    
    if ismember(ixs-20,[1 nx/2 nx])
        start = 1;
    else
        start = nt;
    end
    
    for i = start:nt
        % plot shot record evolution
        ds = zeros(nt,nx);
        %if i <= nt-2
            ds(1:i,:) = data(1:i,:);
        %end
        subplot(2,2,3)
        imagesc(x,t,ds)
        xlabel('X Location'), ylabel('Time (s)')
        title('Shot Record')
        caxis([-0.1 0.1])
        
        % plot wave propagation
        subplot(2,2,4)
        imagesc(x,z,snapshot(1:end-20,21:end-20,i))
        xlabel('Distance (m)'), ylabel('Depth (m)')
        title(['Wave Propagation t = ',num2str(t(i),'%10.3f')])
        caxis([-0.14 0.2])
        
        %writeVideo(vidObj,getframe(gcf));
        drawnow;
    end %shot loop
end
close(vidObj);

%% Traveltime by 2D ray-tracing
% Generate the traveltime field for all z = 0 locations
travelTime = zeros(nz,nx,nx);
subplot(2,2,2)
for ixs = 1:nx
    travelTime(:,:,ixs) = ray2d(V(1:end-20,21:end-20),[1 ixs],dx);
    imagesc(x,z,travelTime(:,:,ixs))
    xlabel('Distance (m)'), ylabel('Depth (m)')
    title(['Traveltime for shot ',num2str(ixs)])
    set(hshot,'XData',x(ixs));
    drawnow
end

% save results for later re-use
%save('saltToothModelData\travelTime.mat', 'travelTime')

%% Process Shots - Kirchhoff Migration
load('travelTime.mat');
Stacked = zeros(nz,nx);
figure(1)
vidObj = VideoWriter('saltToothMigrationKirchhoff.avi');
open(vidObj);
for ixs = 1:nx
    fileName = ['shotfdm',num2str(ixs),'.mat'];
    load(fileName)
    shot = data(21:end-20,:)';
    M = migrate(travelTime,shot,dt,nz,ixs,nx);
    Stacked = Stacked + M;
    
    subplot(2,2,2)
    imagesc(x,z,diff(Stacked,2,1))
    xlabel('Distance (m)'); ylabel('Depth (m)');
    title('Stacked Image');
    caxis([-135 135])
    
    subplot(2,2,3)
    imagesc(x,t,shot)
    xlabel('Distance (m)'); ylabel('Time (s)');
    title(['Current Shot ',num2str(ixs)]);
    caxis([-0.1 0.1])
    
    subplot(2,2,4)
    imagesc(x,t,M)
    xlabel('Distance (m)'); ylabel('Time (s)');
    title(['Current Migrated Shot ',num2str(ixs)]);
    caxis([-5 5])
    
    set(hshot,'XData',x(ixs));
    
    drawnow
    writeVideo(vidObj,getframe(gcf));
end
close(vidObj);

%% Process Shots - Reverse Time Migration
Stacked = zeros(nz+20,nx+40);
show = false;
colormap seismic
%vidObj = VideoWriter('saltToothMigrationRTM.avi');
%open(vidObj);
for ixs = 1:nx
    load(['shotfdm',num2str(ixs),'.mat']);
    load(['snapshot',num2str(ixs),'.mat']);
    shot = diff(data,2,1)';
    %data = [zeros(20,nt); data; zeros(20,nt)];
    tic
    [~, rtmsnapshot] = rtm2d_gpu(V,data,nz,dz,nx,dx,nt,dt);
    tgpu=toc;
    tic
    [~, rtmsnapshot] = rtm2d(V,data,nz,dz,nx,dx,nt,dt);
    tcpu=toc;
    fprintf('Single Threaded CPU: \t%12.4f\n',tcpu);
    fprintf('GPU CUDA Kernel \t\t%12.4f\n',tgpu);
    fprintf('Speedup CPU/GPU: \t\t%12.4f\n',tcpu/tgpu);
    %save(['faultModelData\rtmsnapshot',num2str(ixs),'.mat'],'rtmsnapshot');
    %load(['snapshot',num2str(ixs),'.mat']);
    
    M = 0;
    s2 = 0;
    for i = 1:nt
        M = snapshot(:,:,i).*rtmsnapshot(:,:,nt-i+1)+M;
        s2 = snapshot(:,:,i).^2+s2;
        
        if show
            subplot(2,2,3)
            imagesc(snapshot(1:end-20,21:end-20,i))
            
            subplot(2,2,4)
            imagesc(rtmsnapshot(1:end-20,21:end-20,nt-i+1))
            
            subplot(2,2,2)
            imagesc(diff(M(1:end-20,21:end-20)./s2,2,1))
            
            drawnow
        end
    end
    
    Stacked = Stacked + M;
    subplot(2,2,2)
    imagesc(x,z,diff(Stacked(1:end-20,21:end-20),2,1))
    xlabel('X (m)'); ylabel('Depth (m)');
    title('Stacked Image');
    
    subplot(2,2,3)
    imagesc(x,t,shot)
    xlabel('X (m)'); ylabel('Time (s)');
    title(['Current Shot ',num2str(ixs)]);
    
    subplot(2,2,4)
    imagesc(x,t,M(1:end-20,21:end-20))
    xlabel('X (m)'); ylabel('Time (s)');
    title(['Current Migrated Shot ',num2str(ixs)]);
    
    set(hshot,'XData',x(ixs));    
    drawnow
    
    writeVideo(vidObj,getframe(gcf));
end
close(vidObj)

Contact us