刘宇-作业5-自动化反演-2016.3.31

1、 中国地质大学（武汉）地空学院 姓 名： 刘 宇 班 级： 061132班 学 号： 20131003811 方法原理一、正演计算 水平层状介质对称四极法电侧深的积分关系式(当MN->0时)可表示为 在上式中, 对理论转换电阻率 T1(） 进行积分, 便可以算出理论电阻率s 。 但是直接对该式积分较为困难, 所以采用数字滤波方法运算。经过代换、 褶积等一系列运算后, 可得到离散化的数字逮波形式如下 （2）C为滤波系数。I为取点样号。K为序号，s是位移理论转换电阻率的递推公式 由下而上、如下 （3）滤波系数C通过查找资料给出，观测值也能给出，就能根据公式算出视电阻率的值，通过MATLAB程

2、序实现的结果如下：clearclose allclcformat short gnlayer=3 %层数res(1)=10res(2)=100res(3)=5thick(1)=2thick(2)=8 ab2=2 3 4.5 6 9 12 15 20 30 45 60 90 120 150 200 300 450 600 900 1200 1500 2000 3000 4500 nab2=length(ab2); for i=1:nab2 r=ab2(i); rhos=sdcs1dford(nlayer,res,thick,r); rhoscal(i)=rhos; end ab2=ab2(:);

3、rhoscal=rhoscal(:);clear aa= ab2 rhoscal;dlmwrite('d01obs.dat',a); %画图 loglog(ab2,rhoscal,'r-o',. 'linewidth',2,. 'MarkerEdgeColor','k',. 'MarkerFaceColor','g',. 'MarkerSize',10)axis equalgrid onxlabel('AB/2 / m')ylabel('RHO

4、 / ohm') 副程序function rhos=sdcs1dford(nlayer,res,thick,r)c=0.003042 -0.001198 0.01284 0.0235 0.08688 0.2374 0.6194 1.1817 0.4248 -3.4507 2.7044 -1.1324 0.393 -0.1436 0.05812 -0.02521 0.01125 -0.004978 0.002072 -0.000318;% rhos=10;% whos;for k=1:20 mk=exp(k*log(10)/6-2.1719)/r; %m的值 T(nlayer)=res(

5、nlayer); for i=nlayer-1:-1:1 t1=T(i+1)+res(i); t2=T(i+1)-res(i); t3=exp(-2*mk*thick(i); t4=t1+t2*t3; t5=t1-t2*t3; T(i)=res(i)*t4/t5; %求T的值 end T1(k)=T(1); %限定值的范围end rhos=0;for k=1:20 rhos=rhos+T1(k)*c(k);end2、 反演拟合的实现对于任一条电侧深曲线来说, 总会有某一段曲线对某一层层参数的改变反应特别明显, 而对其他层层参数的改变反应较为迟缓, 也就是说, 一个电性层对应着一段测深曲线段。

6、由此就可以建立电性参数和侧深曲线上的各曲线段一一对应的关系。 在反演拟合中,如果理论曲线和实测曲线的某一段重合不好、 误差大时, 则修改这段理论曲线所对应的层参数,而其他的层参数不变。 如果整条曲线都重合不好,则逐段(层)修改。 对于拟合精度已达到要求的曲线段, 其对应的电性参数则不用修改。 整条曲线所有的曲线段应该修改的都一一修改完毕后,再进行下一步拟合运算, 求出新的拟合精度 即误差, 为下次是否要修改及修改多少作准备。 如此反复修改、 拟合, 就能使整条曲线重合或基本重合。自动化反演可分为如下步骤：构造初始模型先定性解释地下三层水平结构，即三层视电阻率的值和两层厚度的值，用数组m0来代替

7、这5个参数；初始模型正演即由m0中的5个参数得出视电阻率的值，记为d（cal），即视电阻率推算值；计算误差向量Err=(di(obs)-di(cal)2 误差为视电阻率观测值和视电阻率推算值的差的平方和；误差判断可以给误差一些约束，让误差小于某个值时才算合格；自动修正这一步是自动化反演的关键步骤，即计算出模型的修正值，这要用到雅可比反演方程 （J(T)J+I）m=J（T）d J-雅可比偏导数（可近似为J=d/m）-正则化系数 I -单位矩阵d-d(obs)-d(cal)就可以计算出m的值；模型更新 m1=m0+m回到上述步骤通过MATLAB实现程序如下clearclose allclc cle

8、ar aa = dlmread('d01obs.dat');ab2=a(:,1);rhosobs=a(:,2);loglog(ab2,rhosobs,'r-*'); %构造初始模型nlayer=3; res(1)=8;res(2)=50;res(3)=1;thick(1)=1;thick(2)=5; m0=zeros(2*nlayer-1,1);for i=1:nlayer m0(i)=res(i);endfor i=1:nlayer-1 m0(i+nlayer)=thick(i);endm0=log10(m0); %初始模型正演%nab2=length(ab2

9、);for i=1:nab2 r=ab2(i); rhos=sdcs1dford(nlayer,res,thick,r); rhoscal(i,1)=rhos;endrhoscal;hold onloglog(ab2,rhoscal,'b:o');grid on%计算误差向量for i=1:nab2deltd(i,1)=rhosobs(i)-rhoscal(i);enddeltd; err=0;for i=1:nab2 err=err+deltd(i)*deltd(i);endclear stempstemp=sprintf('err=%g',err)title

10、(stemp)%自动修正m=m0;M=2*nlayer-1;lamda=0.05;for iter=1:10 %第一步：计算雅可比偏导数矩阵J J=scomputeJ(m,ab2);% figure(2)% imagesc(J) A=J'*J+lamda*eye(M,M);% figure(3)% imagesc(A) b=J'*deltd;% figure(4)% plot(b,'r-*') deltm=Ab; m deltm; % %模型约束 % for i=1:M if deltm(i) < -0.5 deltm(i) = -0.5; end if

11、deltm(i) > 0.5 deltm(i) = 0.5; end end % %模型更新m(k+1)=m(k)+deltm % m=m+deltm; m; %重新正演 for i=1:nlayer res(i)=10m(i); end for i=1:nlayer-1 thick(i)=10m(nlayer+i); end res thick for i=1:nab2 r=ab2(i); rhos=sdcs1dford(nlayer,res,thick,r); rhoscal(i)=rhos; deltd(i)=rhosobs(i)-rhoscal(i); end norm(delt

12、d); err=0; for i=1:nab2 err=err+deltd(i)*deltd(i); end %画图 clear stemp stemp=sprintf('err=%g',err) title(stemp) hold on loglog(ab2,rhoscal,'b:o'); legend('rhosobs','rhoscal') hold off pause(1)end副程序就是计算雅可比偏导数的，如下function J=scomputeJ(m,ab2)nab2=length(ab2);N=nab2;nm=le

13、ngth(m);nlayer=floor(nm+1)/2);M=nm;J=zeros(N,M);%计算f0for i=1:nlayer res(i)=10m(i);endfor i=1:nlayer-1 thick(i)=10m(i+nlayer);endfor i=1:nab2 r=ab2(i); rhos=sdcs1dford(nlayer,res,thick,r); f0(i)=rhos;end m0=m;for j=1:M m1=m0; m1(j)=m0(j)+0.05; %给模型一个小扰动，用于计算雅可比偏导数 %分解m1到res,thick for i=1:nlayer res(i)=10m1(i); end for i=1:nlayer-1 thick(i)=10m1(i+nlayer); end %正演 for i=1:nab2 r=ab2(i); rhos=sdcs1dford(nlayer,res,thick,r); f1(i)=rhos; J(i,j)=(f1(i)-f