examp.m

来自「这是在网上下的一个东东」· M 代码 · 共 199 行

%VSP velocity inversion example for Chapter 5
clear
rand('state',0);
randn('state',0);

noise = 0.0002;
maxdepth=1000;
depth=0:1:maxdepth;
%true 1-d velocity structure
strue=slowfunctionsmooth(depth);
v=1./strue;

%
% dobs will be a vector containing the depths of the different observations.
% dd is difference in depth between successive observations.
% t contains the noisy observations.
%
%borehole arrival time data
i=0;
dd=20;
for d=dd:dd:1000
  i=i+1;
  dobs(i,1)=d;
  rndstore(i)=randn;
  t(i,1)=quadl(@slowfunctionsmooth,0,d)+noise*rndstore(i);
end
N=i;

figure(1);
bookfonts;
plot(dobs,t,'k');
xlabel('Depth (m)');
ylabel('Arrival Time (s)');
%print -deps2 c5fvspdata.eps

%Number of model parameters
k=1;
M=k*N;

%
% Save this for use in later examples.
%
%save vsp.mat

%System matrix
G=zeros(N,M);
for i=1:N
  G(i,1:k*i) = ones(1,k*i)*dd;
end

%
% Compute the least squares solution.
%
mls=G\t;

%
% Compute the SVD of G.
%
[U,s,V]=svd(G);

%higher order Tikh regularization
[L1,W] = get_l(N,1);
[U1,sm1,X1,V1]=cgsvd(G,L1);
figure(3);
bookfonts;
[reg_corner1,rho,eta,reg_param]=l_curve(U1,sm1,t,'Tikh');
[m1,rho,eta]=tikhonov(U1,sm1,X1,t,reg_corner1);
%print -deps c5flcurve1.eps

[L2,W] = get_l(N,2);
[U2,sm2,X2,V2]=cgsvd(G,L2);
figure(4);
bookfonts;
[reg_corner2,rho,eta,reg_param]=l_curve(U2,sm2,t,'Tikh2');
[m2,rho,eta]=tikhonov(U2,sm2,X2,t,reg_corner2);
%print -deps c5flcurve2.eps

%plot/compare model results
figure(14);
bookfonts;
plot(depth,1000*strue,'k--')
ylabel('True Slowness (s/km)')
xlabel('Depth (m)')
%print -deps c5fvspmod.eps

figure(15);
bookfonts;
plotconst(mls*1000,0,maxdepth);
hold on
plotconstc(mls*1000-1.96*1000*0.0002*sqrt(diag(inv(G'*G))),0,maxdepth,'k--');
plotconstc(mls*1000+1.96*1000*0.0002*sqrt(diag(inv(G'*G))),0,maxdepth,'k--');
ylabel('Slowness (s/km)');
xlabel('Depth (m)');
%print -deps c5fml2.eps


figure(18);
bookfonts;
plotconst(m1*1000,0,maxdepth);
hold on
plot(depth,strue*1000,'k--')
hold off
ylabel('Slowness (s/km)');
xlabel('Depth (m)');
disp(['1st order reg corner is:  ',num2str(reg_corner1)]);
%print -deps c5fmtik1.eps

figure(19);
bookfonts;
plotconst(m2*1000,0,maxdepth);
hold on
plot(depth,strue*1000,'k--')
hold off
ylabel('Slowness (s/km)')
xlabel('Depth (m)')
disp(['2nd order reg corner is:  ',num2str(reg_corner2)]);
%print -deps c5fmtik2.eps

%examine resolution
f2=(diag(s).^2)./(diag(s).^2+reg_corner2^2);

%examine resolution
sdag=inv(s)';

gamma1=sm1(:,1)./sm1(:,2);
f1=gamma1.^2./(gamma1.^2+reg_corner1^2);
F1=diag([f1;ones(N-length(f1),1)]);
R1=X1*F1*inv(X1);

gamma2=sm2(:,1)./sm2(:,2);
f2=gamma2.^2./(gamma2.^2+reg_corner2^2);
F2=diag([f2;ones(N-length(f2),1)]);
R2=X2*F2*inv(X2);

figure(8);
bookfonts;
imagesc(R1);
ylabel('i')
xlabel('j')
H=colorbar;
set(H,'FontSize',18);
colormap('gray');

figure(9);
bookfonts;
imagesc(R2);
ylabel('i')
xlabel('j')
H=colorbar;
set(H,'FontSize',18);
colormap('gray');

%spike resolution example figure
spike=zeros(N,1);
spike(N/2)=1;
r1=R1*spike;
r2=R2*spike;
figure(10);
subplot(2,1,1)
bookfonts
plotconst(r1,0,maxdepth);
ylim([0 0.72])
xlim([0 1000])
H=text(700,.5,'1^{st} Order');
set(H,'FontSize',18);
subplot(2,1,2)
bookfonts
plotconst(r2,0,maxdepth);
ylim([0 0.72])
xlim([0 1000])
xlabel('Depth (m)')
H=text(700,.5,'2^{nd} Order');
set(H,'FontSize',18);
%print -deps c5spikeres.eps

gamma1=sm1(:,1)./sm1(:,2);
gamma2=sm2(:,1)./sm2(:,2);
figure(11);
bookfonts;
semilogy(1:length(gamma1),gamma1.^2./(gamma1.^2+reg_corner1^2),'ko');
hold on
semilogy(1:length(gamma2),gamma2.^2./(gamma2.^2+reg_corner2^2),'k*');
hold off
H=text(15,0.01,'1^{st} order');
set(H,'FontSize',18);
H=text(20,0.00005,'2^{nd} order');
set(H,'FontSize',18);
ylabel('Filter Factors');
xlabel('index, i');
%print -deps c5tikhfilt.eps

dinterp=10:20:990;
strueresamp=spline(depth,strue,dinterp)';

disp(['1st order model 2-norm misfit is ',num2str(norm(m1-strueresamp))]);
disp(['2nd order model 2-norm misfit is ',num2str(norm(m2-strueresamp))]);