forked from seismicreservoirmodeling/SeReM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathESSeisInversionDriver.m
150 lines (137 loc) · 4.37 KB
/
ESSeisInversionDriver.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
%% Ensemble Smoother Seismic Inversion Driver %%
% In this script we apply the Ensemble Smoother inversion method
% (Liu and Grana, 2018) to predict the elastic properties (P- and S-wave
% velocity and density) from seismic data.
%% Available data and parameters
% Load data (seismic data, reference elastic properties, and time)
addpath(genpath('../SeReM/'))
load Data/data3.mat
%% Initial parameters
% number of samples (elastic properties)
nm = size(Snear,1)+1;
% number of samples (seismic data)
nd = size(Snear,1);
% number of variables
nv = 3;
% reflection angles
theta = [15, 30, 45];
ntheta = length(theta);
% time sampling
dt = TimeSeis(2)-TimeSeis(1);
% error variance
varerr = 10^-4;
sigmaerr = varerr*eye(ntheta*(nm-1));
% number of realizations
nsim = 500;
%% Wavelet
% wavelet
freq = 45;
ntw = 64;
[wavelet, tw] = RickerWavelet(freq, dt, ntw);
%% Plot seismic data
figure(1)
subplot(131)
plot(Snear, TimeSeis, 'k', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('Near'); ylabel('Time (s)');
subplot(132)
plot(Smid, TimeSeis, 'k', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('Mid');
subplot(133)
plot(Sfar, TimeSeis, 'k', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('Far');
%% Prior model (filtered well logs)
nfilt = 3;
cutofffr = 0.04;
[b, a] = butter(nfilt, cutofffr);
Vpprior = filtfilt(b, a, Vp);
Vsprior = filtfilt(b, a, Vs);
Rhoprior = filtfilt(b, a, Rho);
mprior = [Vpprior Vsprior Rhoprior];
%% Spatial correlation matrix
corrlength = 5*dt;
trow = repmat(0:dt:(nm-1)*dt,nm,1);
tcol = repmat((0:dt:(nm-1)*dt)',1,nm);
tdis = trow-tcol;
sigmatime = exp(-(tdis./corrlength).^2);
sigma0 = cov([Vp,Vs,Rho]);
%% Prior realizations
Vpsim = zeros(nm, nsim);
Vssim = zeros(nm, nsim);
Rhosim = zeros(nm, nsim);
SeisPred = zeros(nd*ntheta, nsim);
for i=1:nsim
msim = CorrelatedSimulation(mprior, sigma0, sigmatime);
Vpsim(:,i) = msim(:, 1);
Vssim(:,i) = msim(:, 2);
Rhosim(:,i) = msim(:, 3);
[SeisPred(:,i), TimeSeis] = SeismicModel (Vpsim(:,i), Vssim(:,i), Rhosim(:,i), Time, theta, wavelet);
end
% plot of prior models
figure(2)
subplot(131)
plot(Vpsim, Time, 'b', 'LineWidth', 1);
hold on;
plot(Vp, Time, 'k', 'LineWidth', 2);
plot(Vpprior, Time, 'r', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('P-wave velocity (km/s)'); ylabel('Time (s)');
subplot(132)
plot(Vssim, Time, 'b', 'LineWidth', 1);
hold on;
plot(Vs, Time, 'k', 'LineWidth', 2);
plot(Vsprior, Time, 'r', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('S-wave velocity (km/s)');
subplot(133)
plot(Rhosim, Time, 'b', 'LineWidth', 1);
hold on;
plot(Rho, Time, 'k', 'LineWidth', 2);
plot(Rhoprior, Time, 'r', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('Density (g/cm^3)');
%% ESMDA seismic inversion
niter = 4;
alpha = niter; % sum 1/alpha = 1
PriorModels = [Vpsim; Vssim; Rhosim];
SeisData = [Snear; Smid; Sfar];
PostModels = PriorModels;
for j=1:niter
[PostModels, KalmanGain] = EnsembleSmootherMDA(PostModels, SeisData, SeisPred, alpha, sigmaerr);
Vppost = PostModels(1:nm,:);
Vspost = PostModels(nm+1:2*nm,:);
Rhopost = PostModels(2*nm+1:end,:);
for i=1:nsim
[SeisPred(:,i), TimeSeis] = SeismicModel (Vppost(:,i), Vspost(:,i), Rhopost(:,i), Time, theta, wavelet);
end
end
% posterior mean models
mpost = mean(PostModels, 2);
Vpmean = mpost(1:nm);
Vsmean = mpost(nm+1:2*nm);
Rhomean = mpost(2*nm+1:end);
%% Plot results
figure(3)
subplot(131)
plot(Vppost, Time, 'b', 'LineWidth', 1);
hold on;
plot(Vp, Time, 'k', 'LineWidth', 2);
plot(Vpmean, Time, 'r', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('P-wave velocity (km/s)'); ylabel('Time (s)');
subplot(132)
plot(Vspost, Time, 'b', 'LineWidth', 1);
hold on;
plot(Vs, Time, 'k', 'LineWidth', 2);
plot(Vsmean, Time, 'r', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('S-wave velocity (km/s)');
subplot(133)
plot(Rhopost, Time, 'b', 'LineWidth', 1);
hold on;
plot(Rho, Time, 'k', 'LineWidth', 2);
plot(Rhomean, Time, 'r', 'LineWidth', 2);
axis tight; grid on; box on; set(gca, 'YDir', 'reverse');
xlabel('Density (g/cm^3)');