-
Notifications
You must be signed in to change notification settings - Fork 10
/
TV_tomo.m
180 lines (146 loc) · 5.48 KB
/
TV_tomo.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
function [sTV] = TV_tomo(inputs)
%============================================================
% [ss,sg, normError,D] = TV_tomo(inputs)
%
% Implementation of total variation (TV)-regularized tomography, based on
% MTV regularization from Lin et al. 2015. and Chambolle 2004 (please see
% below for citations).
%
% Inputs:
% inputs.lam1: regularization param. 1
% inputs.lamTV: TV regularization parameter
% inputs.tomoMatrix: tomography matrix
% inputs.refSlowness: reference slowness
% inputs.travelTime: travel times
% inputs.validBounds: valid boundaries for TV-regularized inversion
% inputs.normNoise: Euclidian norm of noise vector
% inputs.sTrue: true slowness
% inputs.noiseRealiz: noise realization number (for display purposes
% only)
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The methods implemented here are the same as described in:
% 1. M.J. Bianco and P. Gerstoft, "Travel time tomography with adaptive
% dictionaries," IEEE Trans. on Computational Imaging, Vol. 4, No. 4, 2018.
%
% Implemented here is a TV-regulared tomography method based on MTV
% regularization, developed in:
% 2. Y. Lin and L. Huang, "Quantifying subsurface geophysical properties
% changes using double-difference seismic-waveform inversion with a
% modificed total-variation regularization scheme," Geophysical Journal
% Int., vol. 203, no. 2, 2015.
%
% This implementation of MTV is uses TV image denoising code from
% Chambolle 2004:
% 3. A. Chambolle, "An algorithm for total variation minimization and
% applications," J. Mathematical Imaging and Vision, vol. 20, no. 1--2,
% 2004.
%
% Further implemented is the LSQR sparse least squares algorithm, developed
% in the reference:
% 4. C.C. Paige and M.A. Saunders, "LSQR: And algorithm for sparse linear
% equations and sparse least squares," ACM Trans. on Mathematical Software,
% vol., no. 1, 1982.
%
% If you find this code useful in your research, please cite the above
% references (1-4).
%
% LST version 1.0
% Michael J. Bianco, December 11th, 2018.
% email: [email protected]
%============================================================
lam1=inputs.lam1;
lamTV=inputs.lamTV;
A=inputs.tomoMatrix;
s0=inputs.refSlowness;
Tarr=inputs.travelTime;
vb=inputs.validBounds;
normNoise=inputs.normNoise;
sTrue=inputs.sTrue;
gg=inputs.noiseRealiz;
sTV=zeros(size(sTrue(:))); % initializating reference to be 0
normError = norm(Tarr-A*s0); % initial error (of reference solution)
A = sparse(A);
[nrays,npix] = size(A);
for nn = 1:inputs.solIter
% global slowness
disp(['TV: Realization #',num2str(gg),', Iteration #',num2str(nn)])
dt = Tarr-A*(sTV+s0);
ds = lsqrSOL(nrays,npix,A,dt,lam1,1.0e-6,1.0e-6,1.0e+5,1e3,0);
sg = ds+sTV;
valb = double(vb(:));
n=length(sTrue);
f=sg;
alpha = 0.25; % optimal per chambolle 2004
NIT=5000;
GapTol=1e-2;
verbose=0;%true;
lamCham=2/lamTV; % per Zhu 2010
[u, ~,~, ~,~, ~] ...
= TV_Chambolle(zeros(n),zeros(n),reshape(f,size(sTrue)),lamCham,alpha,NIT,GapTol,verbose);
% ss=u(:);
sTV=u(:).*valb;
s=sTV+s0;
Tref = A*(s); % new reference travel time for global estimate
% calculating errors
normError_new = norm(Tref-Tarr); % error in travel time
normError = [normError,normError_new];
if inputs.plots ==true
[nrow,ncol]=size(sTrue);
Splot = reshape(sg+s0,nrow,ncol);
uplot = reshape(s,nrow,ncol);
figure(inputs.figNo)
clf;
subplot(3,2,1)
imagesc(Splot,inputs.lims)
colormap(gca,'default')
title('$\widehat{\mathbf{s}}_\mathrm{g}$','interpreter','latex','fontsize',16)
bigTitle='Total variation (TV) inversion example';
text(0,-20,bigTitle,'fontsize',15,'interpreter','latex')
h= colorbar;
ylabel(h,'Slowness (s/km)')
xlabel('Range (km)')
ylabel('Range (km)')
subplot(3,2,2)
colormap(gca,'default')
imagesc(uplot,inputs.lims)
h= colorbar;
ylabel(h,'Slowness (s/km)')
xlabel('Range (km)')
ylabel('Range (km)')
title('$\widehat{\mathbf{s}}_\mathrm{TV}$','interpreter','latex','fontsize',16)
% plotting slices
% horizontal slice
sliceLoc = 48;
uSlice = uplot(sliceLoc,:);
sTrue_slice = sTrue(sliceLoc,:);
subplot(3,2,3)
plot(1:100,sTrue_slice,'k',1:100,uSlice,'r')
legend('True','Estimated')
ylabel('Slowness (s/km)')
xlabel('Range (km)')
title('$\widehat{\mathbf{s}}_\mathrm{TV}$, horizontal slice',...
'interpreter','latex','fontsize',16)
% vertical slice
sliceLoc = 30;
uSlice = uplot(:,sliceLoc);
sTrue_slice = sTrue(:,sliceLoc);
subplot(3,2,4)
plot(1:100,sTrue_slice,'k',1:100,uSlice,'r')
legend('True','Estimated')
ylabel('Slowness (s/km)')
xlabel('Range (km)')
title('$\widehat{\mathbf{s}}_\mathrm{TV}$, vertical slice',...
'interpreter','latex','fontsize',16)
subplot(3,2,5)
plot(0:nn,normError,'b-o',0:nn,ones(1,nn+1)*normNoise,'r')
ylabel('Travel time error norm (s)')
xlabel('iteration #')
legend('error','noise')
title('Travel time error vs. iteration',...
'interpreter','latex','fontsize',16)
drawnow
end
end
end