refsgymake.m

这是matlab在地球物理数据处理方面的源码

clear
t = [0:1526]*0.001;
load receivers
rec = receivers;
clear receivers

geom = [100:198]';  % Field REcord file
geom(:,2) = zeros(size(geom));  % SP stations
geom(1,2) = 100.5;  geom(2:99,2) = rec(1:98,4);
geom(:,3) = zeros(99,1);  % First channel station
geom(1,3) = 121;  geom(2:99,3) = (120 + [1:98])';
geom(:,4) = zeros(99,1);  % Last channel station
geom(1,4) = 168; geom(2:99,4) = (167 + [1:98])';
geom(:,5) = zeros(99,1);  % rec receivers reference number
geom(1,5) = 21; geom(2:99,5) = (20 + [1:98])';
geom(:,6) = zeros(99,1);  % SP elevetion
geom(1,6) = rec(1,3);  geom(2:99,6) = rec(1:98,3);
geom(:,7) = zeros(99,1); % SP Northing
geom(1,7) = rec(1,1)-2;  geom(2:99,7) = rec(1:98,1)
geom(:,8) = zeros(99,1);  % SP easting
geom(1,8) = rec(1,2);  geom(2:99,8) = rec(1:98,2);




for i = 1:106;
   if floor(log10(i)) == 0
      filename = ['file000',int2str(i),'.dat']
      [junk,traces] = sg2read(filename);
   elseif  floor(log10(i)) == 1
      filename = ['file00',int2str(i),'.dat']
      [junk,traces] = sg2read(filename);
   else  floor(log10(i)) == 2
      filename = ['file0',int2str(i),'.dat']
      [junk,traces] = sg2read(filename);
   end
   
   traces = traces(:,1:48); % remove aux.
   [m,n] = size(traces);
   trahead = zeros(101,48);
   
trahead(1,:) = 1:n;  	% Trace Sequence Number *
trahead(2,:) = 1:n;  	% Reel trace sequence number *
trahead(3,:) = geom(i,1)*ones(1,n);    	% Orignal Field Record Number *
trahead(4,:) = 1:n;  	% Trace Number within original field record *
trahead(5,:) = geom(i,2)*ones(1,n);	% Energy Source Number (Shot Point Number)
trahead(6,:) = [geom(i,3):geom(i,4)]+geom(i,2);	% CDP Ensemble Number
trahead(7,:) = 1:n;		% Trace number within CDP ensemble
trahead(8,:) = ones(1,n); 	% Trace Identification * see tech info
trahead(9,:) = 8*ones(1,n);	% Number of vertically summed traces (shot stack)
trahead(10,:) = ones(1,n);	% Number of horizontally summed traces	
trahead(11,:) = ones(1,n);	%Data use, 1 = production, 2 = test
start = geom(i,3); finish = geom(i,4);
trahead(12,:) = (10*([start:finish]-geom(i,2))*4); 	% Trace-Shot Point Offset in dmeters
trahead(13,:) = 10*rec([geom(i,5):geom(i,5)+47],3)';	% Receiver Elevation (properly from sea level)
trahead(14,:) = 10*geom(i,6)*ones(1,n);	% Source Elevation at surface
trahead(15,:) = zeros(1,n);	% Source Depth (a positive number)
trahead(16,:) = zeros(1,n);	% Receiver Datum elevation
trahead(17,:) = zeros(1,n); 	% Source Datum elevation
trahead(18,:) = zeros(1,n);	% Water Depth at Source
trahead(19,:) = zeros(1,n);	% Water Depth at Group
trahead(20,:) = -10*ones(1,n);	% Scaler to multiply elevations and depths above
trahead(21,:) = -10*ones(1,n);	% Scaler to multiply co-ordinates below
				% Note, scalars are plus/minus 1, 10, 100, ...
				% Positive means multiply, negative means divide
trahead(22,:) = 10*geom(i,8)*ones(1,n);	% X source co-ordinate
trahead(23,:) = 10*geom(i,7)*ones(1,n);	% Y Source co-ordinate
trahead(24,:) = 10*rec([geom(i,5):geom(i,5)+47],2)';		% X receiver co-ordinate
trahead(25,:) = 10*rec([geom(i,5):geom(i,5)+47],1)';;	% Y receiver co-ordinate
trahead(26,:) = ones(1,n);	% Co-ordinate units, 1(meter or feet),2(seconds of arc)
trahead(27,:) = zeros(1,n);	% Weatering velocity
trahead(28,:) = zeros(1,n);	% Subweathering velocity
trahead(29,:) = zeros(1,n);	% Uphole time at source
trahead(30,:) = zeros(1,n);	% Uphole time at receiver
trahead(31,:) = zeros(1,n);	% Source static correction
trahead(32,:) = zeros(1,n);	% Reciever static correction
trahead(33,:) = zeros(1,n);	% Total static applied (zero if none)
trahead(34,:) = zeros(1,n);	% Lag time A in ms
trahead(35,:) = zeros(1,n);	% Lag time B in ms
trahead(36,:) = zeros(1,n);	% Delay Recording time in ms
trahead(37,:) = zeros(1,n); 	% Mute time Start
trahead(38,:) = zeros(1,n);	% Mute time end
trahead(39,:) = m*ones(1,n);	% Number of samples in trace *
trahead(40,:) = (t(2)-t(1))*1e6*ones(1,n); % Sample interval in microseconds for this trace
trahead(41,:) = ones(1,n);	% Gain type of field instruments, see tech standards
trahead(42,:) = ones(1,n);	% Instrument gain constant
trahead(43,:) = ones(1,n);	% Instrument early gain
trahead(44,:) = ones(1,n);	% Correlated 1 = no, 2 = yes
trahead(45,:) = zeros(1,n);	% Sweep freq at start
trahead(46,:) = zeros(1,n);	% Sweep freq at end
trahead(47,:) = zeros(1,n);	% Sweep length
trahead(48,:) = zeros(1,n);	% Sweep type, see tech standards
trahead(49,:) = zeros(1,n);	% Sweep taper at start in ms
trahead(50,:) = zeros(1,n);	% Sweep taper at end in ms
trahead(51,:) = ones(1,n);	% Sweep taper type
trahead(52,:) = zeros(1,n);	% Alias filter freqency if used
trahead(53,:) = zeros(1,n);	% Alias filter slope if used
trahead(54,:) = zeros(1,n);	% Notch filter frequency if used
trahead(55,:) = zeros(1,n);	% Notch filter slope
trahead(56,:) = zeros(1,n);	% Low cut if used
trahead(57,:) = zeros(1,n);	% High Cut if used
trahead(58,:) = zeros(1,n);	% Low cut slope
trahead(59,:) = zeros(1,n);	% High cut slope 
trahead(60,:) = 2000*ones(1,n); % Data recorded year
trahead(61,:) = 831*ones(1,n);	% Day of year 

filename(10:12) = 'sgy'
fid = sgywrit2(filename,traces,t,trahead);
   
end