examp.m

这是在网上下的一个东东

%VSP velocity inversion example for Chapter 5
%(originally suggested by Gary Pavlis, Indiana University).
clear
rand('state',0);
randn('state',0);
load vsp

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

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

%tril(ones(N,N))*dd;

mls=(G\t);

[U,s,V]=svd(G);

%higher order Tikh regularization
[L1,W] = get_l(N,1);
[U1,sm1,X1,V1]=cgsvd(G,L1);
figure(2)
bookfonts;
[reg_corner1,rho1,eta1,reg_param]=l_curve(U1,sm1,t,'tsvd');
loglog(rho1,eta1,'ko');
reg_corner1=8;
[m1,rho,eta]=tgsvd(U1,sm1,X1,t,reg_corner1);
loglog(rho1,eta1,'ko');
xlim([1e-4,1e-2]);
xlabel('residual norm || G m - d ||_2');
ylabel('solution semi-norm || L m ||_2');
for i=[1 10 20 30 40],
  H=text(rho1(i),eta1(i),['    ',num2str(i)]);
  set(H,'FontSize',18);
end
%print -deps2 c5flcurve1tsvd.eps

[L2,W] = get_l(N,2);
[U2,sm2,X2,V2]=cgsvd(G,L2);
figure(3)
bookfonts;
[reg_corner2,rho2,eta2,reg_param]=l_curve(U2,sm2,t,'tsvd');
reg_corner2=4;
[m2,rho,eta]=tgsvd(U2,sm2,X2,t,reg_corner2);
loglog(rho2,eta2,'ko');
xlim([1e-4,1e-2]);
xlabel('residual norm || G m - d ||_2');
ylabel('solution semi-norm || L m ||_2');
for i=[1 10 20 30 40],
  HH=text(rho2(i),eta2(i),['    ',num2str(i)]);
  set(HH,'fontsize',18);
end
%print -deps2 c5flcurve2tsvd.eps

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

figure(5)
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 -deps2 c5fmtik2tsvd.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))]);