%十层模型反演两层地形
clc;
tic;
I0=1;
a=100;
load Hdt.mat;%Hdt为理论响应
load dat.txt;%dat为理论模型
load dat0.txt;%dat0为初始模型
M=dat;
M0=dat0;
nt=71;%时间道
n1=1:nt;
%disp(n1);
t=logspace(-6,1,nt)';
c2=dat0(1);
p0=M0(2:(c2+1));
H0=M0((c2+2):(2*c2));
y=zeros(nt,c2);%定义雅克比矩阵大小
b=zeros(nt,1);%定义b的大小
x=zeros(1,c2);
Hdt0=myfunction(p0,H0);
Hdt0=Hdt0';
ff=m0/(4*pi)*(2*pi*I0*a*a*m0/5)^(2/3);
ps0(n1)=ff/t(n1)*(1./(t(n1)*Hdt0(n1)));
err0=(norm((ps-ps0)/ps))^2/nt;
%error(0)=err0;
%load dat1.txt;
%M0=dat1;
%x2=[5,30,80,40,30];%定义模型下限
%x1=[150,110,120,120,120];%定义模型上限
%初始模型视电阻率
%c=c1;
%计算模型变化值x
k=0; %迭代次数while k<5%err0>0.05
k=k+1;
err1=err0;
while k<5%err0>0.05
for j=1:2*c-1
p1(j)=p0(j);
H1(j)=H0(j);
Hdt1=myfunction(p1,H1);
ps1(n1)=ff/t(n1)*(1./(t(n1)*Hdt1(n1)));
err0=norm(abs(ps)-abs(ps1));
err2=norm(abs(ps)-abs(ps2));
%计算导数
z=err2-err0;
y(j)=1001*z./M(j);
[U,W,V]=svd(y);%奇异值分解
end
P=-(y./norm(y));
P=P';
arf=0.5*err0/(y*P);
ll=ll-arf*P'; %更新起始点ll
p1=ll(1:c1);
H1=ll((c1+1):(2*c1-1));
Hdt3=myfunction(p1,H1);
ps1(n1)=ff/t(n1)*(1./(t(n1)*Hdt3(n1)));
err0=norm(abs(ps)-abs(ps1));
error(k)=err0;
M=ll;
while err0>err1
arf=0.5*err0/(y*P);
ll=ll-arf*P'; %更新起始点ll
M=ll;
p1=ll(1:c1);
H1=ll((c1+1):(2*c1-1));
Hdt3=myfunction(p1,H1);
ps1(n1)=ff/t(n1)*(1./(t(n1)*Hdt3(n1)));
err0=norm(abs(ps)-abs(ps1));
error(k)=err0;
end
end
disp(x);
M2=[c,p0,H0];
%显示出反演结果于命令行窗口
disp(M2);
%Hdt4=myfunction(p4,H4);
%Hdt4=Hdt4';
Hdt0=Hdt0';
figure(2);
loglog(t,Hdt0,'-');
figure(3);
loglog(1:k,error);xlabel('k');ylabel('err');title('收敛误差曲线');
disp(M0);
toc;
包括一系列正反演程序_研究正反演_matlab_正反演
版权申诉
5星 · 超过95%的资源 31 浏览量
2022-04-02
10:49:25
上传
评论
收藏 3KB RAR 举报
阿里matlab建模师
- 粉丝: 3233
- 资源: 2782