Previous   |    Next

doi:10.2204/iodp.pr.340T.2012

Appendix

Matlab scripts

The Matlab scripts in this section were used to make a preliminary assessment of data quality from the Expedition 340T VSP. The scripts present time-windowed and band-pass filtered versions of all sensor components, shot-by-shot, for a single VSP station. The window is a short time interval based upon time integration of a 1-D velocity model for the hole (see “Principal results”). Time integration of the sonic log data gives essentially the same predictions. The main script is “cull_vsp_data.” It uses “rdsegy_full” to ready the SEGY versions of the VSP data and “butter_filt” to apply a band-pass Butterworth causal filter.

cull_vsp_data

% CULL_VSP_DATA - script to go through 3 component VSP data with traces

% selected by station ID and plot all components for all shots at a given station

% depth one shot at a time.

% The VSP data from Expedition 340T was sufficiently noisy that arrivals could

% only be identified by focusing on a short time window centred on the expected

% arrival time and comparing filtered and unfiltered versions of the traces of

% a shot

%

% Shot-to-shot variability was high enough that displaying all traces for all

% shots at a single station was not feasible.

%

% Based on this display a initial data quality assessment was made for shots

% from the 340T VSP.

%

% Required files

% SEG-Y files with x,y & z data

% Station text file with 3 columns station ID, depth (m), predicted time(s)

% a. j. harding

stationFile = 'vsp_time_pred.txt'; % 3 column table of predicted arrival times

vspDir = '../Expedition_340T/data/vsp'; % Directory with SEG-Y data

xFile = '340T-U1309D_raw_shot_geo_x.segy'; % Actual files

yFile = '340T-U1309D_raw_shot_geo_y.segy';

zFile = '340T-U1309D_raw_shot_geo_z.segy';

% Length of time display window relative to time precition

windowB = 0.1; % before pick

windowT = 0.2; % after pick

time_pred_offset = 0.010; % Adjustment of prediction for gun delay, gun depth etc.

% Read VSP station file with IDs, depths & time prediction

data = load(stationFile);

station_id = data(:,1);

station_z = round(data(:,2)*1000); % convert m

station_tpred = data(:,3);

max_ID = max(station_id); % Maximum station #

if ~exist('seisz','var')

[seisz,header] = rdsegy_full(fullfile(vspDir,zFile));

seisx = rdsegy_full(fullfile(vspDir,xFile));

seisy = rdsegy_full(fullfile(vspDir,yFile));

ixDepth = find(header.lactive == 11); % Row with station elevation wrt sea level

ixScalar = find(header.sactive == 35); % Scalar for elevation

elevScalar = double(header.short(ixScalar,1)); % Assume it doesn't change

if elevScalar < 0; elevScalar = -1/elevScalar; end

ixSampRate = find(header.sactive == 59);

dt = double(header.short(ixSampRate,1))/1000000.; % convert from INT16

end

% Set up parameters for Butterworth low/highpass filtering

fc = [10,60]; % corner frequencies in Hz

[Bl,Al] = butter_filt(4,fc(2),dt,'low'); % 4-pole low

[Bh,Ah] = butter_filt(3,fc(1),dt,'high'); % 3-pole high

% Read in the 3 components of data & header information

% Only do this once per session

nt = size(seisz,1); % samples

ntotal = size(seisz,2); % # of shots/traces

tv = [0:nt-1]*dt; % time vector for x-axis of plot

% Apply filters as two part cascade low then high

filt_seisx = filter(Bl,Al,seisx);

filt_seisy = filter(Bl,Al,seisy);

filt_seisz = filter(Bl,Al,seisz);

filt_seisx = filter(Bh,Ah,filt_seisx);

filt_seisy = filter(Bh,Ah,filt_seisy);

filt_seisz = filter(Bh,Ah,filt_seisz);

% Use depth from SEG-Y header to assign station IDs to all shots

% used to select data for plotting.

% NOTE must convert to double from INT32

zdata = -double(header.long(7,:))/10000.; % Station depths from segy headers

station_data = ones(ntotal,1);

for j = 1:ntotal

ik = find(abs(station_z - zdata(j)) < 0.5);

station_data(j) = station_id(ik);

end

while(1)

curStation = input('Enter current station (<=0 to finish) ');

if (curStation > max_ID); continue; end

if (curStation <= 0); break; end

% Extract station depth & travel time prediction

ixcur = find(station_id == curStation);

curZ = station_z(ixcur);

curT = station_tpred(ixcur) + time_pred_offset;

% Calculate sliding window around time prediction

tbeg = floor((curT-windowB)/0.1)*0.1;

tfin = ceil((curT+windowT)/0.1)*0.1;

ixb = round(tbeg/dt)+1;

ixt = round(tfin/dt)+1;

ixRange = [ixb:ixt];

% Select traces for the current station

ixselect = find(station_data == curStation);

nselect = length(ixselect);

tv_loc = tv(ixRange);

xTraces = seisx(ixRange,ixselect)';

yTraces = seisy(ixRange,ixselect)';

zTraces = seisz(ixRange,ixselect)';

fxTraces = filt_seisx(ixRange,ixselect)';

fyTraces = filt_seisy(ixRange,ixselect)';

fzTraces = filt_seisz(ixRange,ixselect)';

% Create clipped version of traces

clipAmp = 1.;

ixClip = find(abs(fxTraces(:)) > clipAmp);

cxTraces = fxTraces;

cxTraces(ixClip) = NaN;

ixClip = find(abs(fyTraces(:)) > clipAmp);

cyTraces = fyTraces;

cyTraces(ixClip) = NaN;

ixClip = find(abs(fzTraces(:)) > clipAmp);

czTraces = fzTraces;

czTraces(ixClip) = NaN;

% Now loop through traces plotting 3 windows

% Unfiltered data - autoscaled by Matlab

% Filtered - data also autoscaled

% Clipped - y axis will be no larger than clipAmp

for j = 1:nselect

titleString = sprintf('x Comp Station: %d at z: %d Trace %d', ...

curStation,curZ,ixselect(j));

figure(1) % Unfiltered

clf

subplot(3,1,1)

plot(tv_loc, xTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13); % time prediction

xlim([tv_loc(1),tv_loc(end)]);

title(titleString)

subplot(3,1,2)

plot(tv_loc, yTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title('y Component');

subplot(3,1,3)

plot(tv_loc, zTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title('z Component');

figure(2) % Filtered

clf

subplot(3,1,1)

plot(tv_loc,fxTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title(titleString)

subplot(3,1,2)

plot(tv_loc, fyTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title('y Component');

subplot(3,1,3)

plot(tv_loc, fzTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title('z Component');

figure(3) % clipped

clf

subplot(3,1,1)

plot(tv_loc,cxTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title(titleString)

subplot(3,1,2)

plot(tv_loc, cyTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title('y Component');

subplot(3,1,3)

plot(tv_loc, czTraces(j,:),'k');

hold on

plot(curT,0,'*r','markerSize',13);

xlim([tv_loc(1),tv_loc(end)]);

title('z Component');

pause % Wait for user to press key to continue

end % for nselect

end % while

rdsegy_full

function [seis,header,binh] = rdsegy_full(fname)

% RDSEGY_FULL returns both SEG-Y data and a header structure with active entries

%

% [seis,header,binh] = rdsegy_full(fname)

%

%

% Inputs

% fname: input file name

%

% Outputs:

% seis[nsamps, nx]: Seismogram data

% header - structure containing active header entries from the trace

% header (defined as a non-zero value for any trace).

% .sactive - indices in header of active short integer entries

% .lactive - " " " " long " "

% .short(nactiveS,nx) - matrix with short header entries

% .long(nactiveL, nx) - matrix with long header entries

%

% a. j. harding, March 2012

% Indices for elements of the SEGY trace header

insamp = 58; % Header number of samples

% Open file and position past the EBCDIC header

fid = fopen(fname, 'r','ieee-be');

if fid == -1

error('*** RDSEGY_FULL *** Error opening seismogram file')

end

fseek(fid,0,1);

totalBytes = ftell(fid);

% Skip over the EBCDIC Header

fseek(fid,3200,-1);

% Read the Binary header & check data format

% The data format is stored in element 13 of the binary header. By

% convention Sioseis denotes data in host machine floating point format

% with an ID >= 5. This routine assumes floating point format and does

% format conversion

binh = fread(fid,200,'short');

data_format = binh(13);

switch data_format

case 1 % IBMFP

disp('Format is IBM fp')

form_str = 'uint';

form_str = 'float32';

data_size_bytes = 4;

case 2

form_str = 'int32';

data_size_bytes = 4;

case 3

form_str = 'int16';

data_size_bytes = 2;

case 5 % Native floating point

form_str = 'float32';

data_size_bytes = 4;

otherwise

fprintf(2,'RDSEGY2 BINARY HDR Cannot read data format\n');

fprintf(2,'Data format %d\n',data_format);

error('**EXIT**');

end

disp('***WARNING*** forcing format to native floating point');

form_str = 'float32';

data_format = 5;

% Read the 1st header to find the no. of samples/trace

T1 = fread(fid,120,'short')';

nsamps = T1(insamp);

fseek(fid,-240,0);

% Estimate total number of traces from no. of samples

ntrEst = (totalBytes - 3600)/(240+nsamps*data_size_bytes);

% Preassign an initial storage space for trace headers/data.

% (Preassignment significantly increases read speed)

% The header for each trace is read twice and once as short ints and

% once as long ints, thus there are 2 header buffers ibuf & lbuf.

ntr = max(1000,ntrEst); % Size of buffer increments

ibuf = zeros(120,ntr);

lbuf = zeros( 60,ntr);

%bseis = zeros(nsamps,ntr);

seis = zeros(nsamps,ntr);

% >>>>>>>> Main Reading Loop : Read all traces <<<<<<<<

nx = 0; % Number of traces read

tic

while (1)

T1 = fread(fid,120,'int16');

if (feof(fid) == 1) break; end

nsamps = T1(insamp);

fseek(fid,-240,0);

T2 = fread(fid, 60,'int32');

T3 = fread(fid,nsamps,form_str); % read trace data

if data_format == 1;

T3 = ibm2num(uint32(T3));

end

nx = nx + 1;

% Allocate additional storage space for data arrays

if (rem(nx,ntr) == 1 & nx > 2)

seis = [seis,zeros(nsamp,1000)];

ibuf = [ibuf,zeros(120,1000)];

lbuf = [lbuf,zeros( 60,1000)];

end

% keyboard

ibuf(:,nx) = T1;

lbuf(:,nx) = T2;

seis(:,nx) = T3;

end

%fprintf(1,'There are %d seismograms of %d points each\n',nx,nsamps);

fclose(fid);

seis = seis(:,1:nx);

% 27L - used by UTIG as 4-byte lag value - non standard

long_ix = [1:7,10:17,19:22]; % As per SEG-Y standard

% 47 - SEG-Y standard subweathering velocity

% 50 - source static correction

% 51 - group static correction

% 52 - total static correction

% 53 - lag time A -time difference between source header time and time break

% 54 - lag time B - time between time break and initiation of source

% 55 - deep water delay time between energy source and time when recording starts

short_ix = [15:18,35:36,45:59 63:90]; % As per standard

% Last 60 bytes contain optional entries

aux_short = [ 91:120];

aux_long = [ 46:60];

header.sactive = intersect(find(any(ibuf,2)),union(short_ix,aux_short));

header.lactive = intersect(find(any(lbuf,2)),union(long_ix,aux_long));

header.short = int16(ibuf(header.sactive,1:nx));

header.long = int32(lbuf(header.lactive,1:nx));

read_time = toc;

if nx > 1000

fprintf(1,'Read %d traces in %.3f s\n',nx,read_time);

end

butter_filt

function [B,A] = butter_filt(n,fc,dt,type)

% BUTTER_FILT - creates a butterworth lowpass filter coefficients of order n

%

% [B,A] = butter(filt(n,fc,dt,type)

%

% n - order of the filter.

% The drop off of the filter is approximately 6*n dB/octave

%

% fc - corner frequency(s) of filter.

% dt - sample rate of data

%

% type - 'low', 'high','bandpass'

%

% B,A - numerator & denominator of IIR filter. To apply filter use the

% Matlab function

% To apply the filter use the Matlab function filter

% y = filter(B,A,x)

%

% Implementation based on Wikipedia articles on Butterworth filter & bilinear

% z-transform.

% Though note bilinear tranform implied here is

% s <- 2/dt (1-z)/(1+z)

% As I want stable filter in terms of z not z^(-1)

%

% a. j. harding

fcnorm = fc*dt;

if (n <= 0)

error('Order of filter (%d) must be > 0',n);

end

if (any(fcnorm < 0) | any(fcnorm >= 0.5))

error('Corner frequency %.1f outside range [0,fnyq]');

end

if (nargin < 4)

type = 'low';

end

if ~any(strcmp(type,{'low','high','bandpass'}))

error('Unrecognized filter type %s',type);

end

nby2 = floor(n/2);

theta = (2*[1:nby2]+n-1)*pi/2/n; % roots on unit circle

% fa = 1/pi/dt * tan(pi*fcnorm); % pre-warped corner frequency

av = 1./tan(pi*fcnorm); % scaling factor based on corner frequency

if any(strcmp(type,{'low','bandpass'}))

a = av(end);

asq = a*a;

asqp = asq+1;

if (mod(n,2) == 1) % Odd power => pole at -1

scalef = (1+a);

B = [1,1];

A = [1, (1-a)/scalef];

else

scalef = 1.;

B = 1;

A = 1;

end

% Add response due to other poles as conjugate pairs for effieciency

% & to avoid complex coefficients/calculations.

for k = 1:nby2

scalep = asqp - 2*a*cos(theta(k));

B = conv(B,[1,2,1]);

A = conv(A,[1,-2*(a^2-1)/scalep, (asqp+2*a*cos(theta(k)))/scalep]);

scalef = scalef * scalep;

end;

B = B/scalef;

if strcmp(type,'low')

return

else

Bl = B;

Al = A;

end

end

% Now do high pass end

a = av(1);

asq = a*a;

asqp = asq+1;

if (mod(n,2) == 1) % Odd power => pole at -1

scalef = (1+a);

B = a*[1,-1];

A = [1, (1-a)/scalef];

else

scalef = 1.;

B = 1;

A = 1;

end

% Add response due to other poles as conjugate pairs for effieciency

% & to avoid complex coefficients/calculations.

for k = 1:nby2

scalep = asqp - 2*a*cos(theta(k));

B = asq*conv(B,[1,-2,1]);

A = conv(A,[1,-2*(a^2-1)/scalep, (asqp+2*a*cos(theta(k)))/scalep]);

scalef = scalef * scalep;

end;

B = B/scalef;

if strcmp(type,'bandpass')

B = conv(B,Bl);

A = conv (A,Al);

end

Top of page   |    Previous   |    Next