Geophysics and Seismology Using MATLAB
MATLAB is a powerful tool for geophysics and seismology applications. It provides functions for modeling, data analysis, and visualization, which are critical for interpreting seismic data, solving wave equations, and more. This section covers several aspects of geophysics and seismology using MATLAB, with examples.
Applications of Geophysics and Seismology Using MATLAB
1. Seismic Data Analysis, Reading and Plotting Seismic Data
Seismic data is typically represented as time series or waveforms. MATLAB provides tools to analyze and visualize such data.
% Reading seismic data
time = 0:0.01:10; % Time vector in seconds
amplitude = sin(2*pi*1*time) + 0.5*sin(2*pi*3*time); % Simulated waveform
plot(time, amplitude);
xlabel('Time (s)');
ylabel('Amplitude');
title('Simulated Seismic Waveform');
This code simulates a seismic waveform and plots it:
A sinusoidal waveform combining two frequencies is visualized, resembling a seismic signal.
Example: Spectral Analysis
Fourier Transform is widely used in geophysics to analyze frequency components in seismic data.
% Spectral analysis of seismic data
Fs = 100; % Sampling frequency in Hz
L = length(amplitude);
Y = fft(amplitude);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L; % Frequency vector
plot(f, P1);
xlabel('Frequency (Hz)');
ylabel('|P1(f)|');
title('Single-Sided Amplitude Spectrum of Seismic Data');
2. Wave Propagation
Wave propagation is central to seismology. MATLAB can model wave propagation using finite-difference or analytical methods.
The 1D wave equation is given by:
∂²u/∂t² = c² ∂²u/∂x²
Here is MATLAB code to simulate wave propagation:
% Simulating 1D wave propagation
c = 3000; % Wave speed in m/s
Lx = 1000; % Length in meters
Nx = 100; % Number of spatial points
dx = Lx/Nx; % Spatial step
dt = 0.001; % Time step
Nt = 200; % Number of time steps
x = linspace(0, Lx, Nx); % Spatial grid
u = zeros(Nx, Nt); % Solution matrix
u(Nx/2, 1) = 1; % Initial condition: impulse at center
for t = 2:Nt-1
for i = 2:Nx-1
u(i, t+1) = 2*u(i, t) - u(i, t-1) + c^2*(dt/dx)^2*(u(i+1, t) - 2*u(i, t) + u(i-1, t));
end
end
imagesc(u'); % Display the wave propagation
xlabel('Position');
ylabel('Time');
title('1D Wave Propagation');
This simulation shows the propagation of a wave over time, starting with an initial impulse.
3. Earthquake Modeling, Moment Tensor Visualization
MATLAB can simulate earthquake sources, moment tensors, and fault plane solutions.
% Visualizing earthquake moment tensor
M = [1 -0.5 0; -0.5 1 0; 0 0 -2]; % Moment tensor
[eigenvectors, eigenvalues] = eig(M);
quiver3(0, 0, 0, eigenvectors(:,1), eigenvectors(:,2), eigenvectors(:,3));
title('Moment Tensor Visualization');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Z-axis');
Advance Examples in Geophysics and Seismology Using MATLAB
Example 1: Inversion of Seismic Data
Inversion is a key technique in geophysics for estimating subsurface properties from observed data. The following example demonstrates a basic seismic inversion using least-squares fitting:
% Seismic data inversion using least squares
% Observed data (d) and forward model matrix (G)
G = [1 2; 3 4; 5 6]; % Forward model
d = [10; 20; 30]; % Observed data
% Least squares inversion
m = pinv(G) * d; % Compute model parameters
disp('Estimated Model Parameters:');
disp(m);
% Check forward model prediction
predicted_d = G * m;
disp('Predicted Data:');
disp(predicted_d);
The output shows the estimated model parameters and the predicted data, which can be compared to the observed data for validation.
Example 2: 2D Finite-Difference Wave Propagation
This example simulates wave propagation in 2D using the finite-difference method to solve the acoustic wave equation:
% 2D Wave propagation using finite-difference
nx = 100; ny = 100; % Grid size
dx = 10; dy = 10; % Spatial steps
nt = 200; dt = 0.001; % Time steps and duration
c = 1500; % Wave speed in m/s
u = zeros(nx, ny); % Current wavefield
u_prev = zeros(nx, ny); % Previous wavefield
u_next = zeros(nx, ny); % Next wavefield
% Initial condition: source in the center
src_x = nx/2; src_y = ny/2;
u(src_x, src_y) = 1;
% Time-stepping loop
for t = 1:nt
for i = 2:nx-1
for j = 2:ny-1
u_next(i, j) = 2*u(i, j) - u_prev(i, j) + (c*dt/dx)^2 * (u(i+1, j) + u(i-1, j) + u(i, j+1) + u(i, j-1) - 4*u(i, j));
end
end
u_prev = u;
u = u_next;
end
imagesc(u); colormap('hot');
title('2D Wavefield');
xlabel('X'); ylabel('Y');
The output is a dynamic visualization of the wave propagating through the medium, useful for studying wave behavior in different subsurface structures.
Example 3: Earthquake Fault Slip Simulation
This example simulates fault slip during an earthquake using the Okada model, which calculates surface displacement due to faulting.
% Okada model for fault slip simulation
x = linspace(-100, 100, 100);
y = linspace(-100, 100, 100);
[X, Y] = meshgrid(x, y);
% Fault parameters
depth = 10; % Depth in km
slip = 1; % Slip in meters
strike = 90; % Strike angle in degrees
dip = 45; % Dip angle in degrees
% Compute surface displacements
[Ux, Uy, Uz] = okada_displacement(X, Y, depth, slip, strike, dip);
% Visualization
quiver(X, Y, Ux, Uy);
title('Surface Displacement Due to Fault Slip');
xlabel('X (km)'); ylabel('Y (km)');
The quiver plot shows the displacement vectors on the surface due to fault slip. This model is valuable in understanding fault behavior during earthquakes.
Example 4: Forward Modeling of Gravity Anomalies
Forward modeling helps calculate the expected gravity anomalies caused by subsurface density variations.
% Forward modeling of gravity anomaly
x = linspace(-1000, 1000, 200); % Observation points in meters
z = 500; % Depth of the source in meters
density_contrast = 500; % kg/m^3
G = 6.67430e-11; % Gravitational constant
% Gravity anomaly calculation
anomaly = (2 * G * density_contrast * z^2) ./ ((x.^2 + z^2).^(3/2));
% Plotting the anomaly
plot(x, anomaly);
title('Gravity Anomaly');
xlabel('Distance (m)');
ylabel('Anomaly (m/s^2)');
This plot shows the gravity anomaly profile over a buried mass, useful for identifying subsurface structures like salt domes or mineral deposits.
Example 5: Seismic Tomography
Seismic tomography reconstructs subsurface velocity models from travel-time data. Here is an example:
% Seismic tomography example
% Observed travel times and path matrix
G = [1 0 1; 1 1 0; 0 1 1]; % Ray paths
t_obs = [2.1; 3.2; 1.5]; % Travel times (s)
% Solve for slowness (inverse of velocity)
s = pinv(G) * t_obs; % Slowness vector
v = 1 ./ s; % Velocity vector
disp('Estimated Velocities:');
disp(v);
The estimated velocities represent different regions of the subsurface traversed by seismic waves.
Useful MATLAB Functions for Geophysics
Practice Questions
Test Yourself
1. Simulate a 2D wave propagation using MATLAB.
2. Perform spectral analysis on real seismic data and interpret the results.
3. Compute and visualize eigenvalues and eigenvectors for a given earthquake moment tensor.