fanyan.asv

用SVD奇异值分解法和阻尼最小二乘算法解决地球物理反问题

%本程序可以实现用svd分解法和阻尼最小二乘法来进行反演，如果choose=0，实现svd分解法，如果choose=1，实现阻尼最小二乘法
function [G,D,DD,M,MM]=fanyan
clear;
x=0:0.02:1;
y=0:1:10;
n=length(x);
m=length(y);
G=zeros(m,n);
choose=input('if you want to use SVD method,choose=0,then damped least-square method,choose=1,so choose=?')
zuni=0.0001;
pi=3.14;
%计算核函数
for i=1:m
    for j=1:n
        G(i,j)=exp(-i*x(j));
    end
end
%给定理论模型
M=zeros(n,1);
for i=1:n
    M(i,1)=1-0.5*cos(2*pi*x(i));
end
D=G*M;
%在理论计算的数据上加噪声得到观测值
N=0.001*randn(m,1);%因为噪声具有随机性，所以得出的结果也可能不一样
DD=D+N; 
%两种反演方法的选择
[U,S,V]=svd(G);
 if   choose==0
     S(11,11)=0; 
     S(10,10)=0;
     S(9,9)=0;
     S(8,8)=0;
     S(7,7)=0;
     SS=pinv(S);
 else
     [U1,S1,V1]=svd(G');
    SS=pinv(S1*S+zuni*eye(n,n))*S1;
 end   
MM=V*SS*U'*DD;%计算反演模型
figure;
plot(x,M,x,MM,'r-*');
legend('理论模型曲线','反演模型曲线');