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

中国地质大学（武汉）地空学院 姓 名： 刘 宇 班 级： 061132班 学 号： 20131003811 方法原理一、正演计算 水平层状介质对称四极法电侧深的积分关系式(当MN-0时)可表示为 在上式中, 对理论转换电阻率 T1(λ） 进行积分, 便可以算出理论电阻率ρs 。 但是直接对该式积分较为困难, 所以采用数字滤波方法运算。经过代换、 褶积等一系列运算后, 可得到离散化的数字逮波形式如下 （2）C为滤波系数。I为取点样号。K为序号，s是位移理论转换电阻率的递推公式 由下而上、如下 （3）滤波系数C通过查找资料给出，观测值也能给出，就能根据公式算出视电阻率的值，通过MATLAB程序实现的结果如下：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(:);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 / 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(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);end反演拟合的实现对于任一条电侧深曲线来说, 总会有某一段曲线对某一层层参数的改变反应特别明显, 而对其他层层参数的改变反应较为迟缓, 也就是说, 一个电性层对应着一段测深曲线段。 由此就可以建立电性参数和侧深曲线上的各曲线段一一对应的关系。 在反演拟合中,如果理论曲线和实测曲线的某一段重合不好、 误差大时, 则修改这段理论曲线所对应的层参数,而其他的层参数不变。 如果整条曲线都重合不好,则逐段(层)修改。 对于拟合精度已达到要求的曲线段, 其对应的电性参数则不用修改。 整条曲线所有的曲线段应该修改的都一一修改完毕后,再进行下一步拟合运算, 求出新的拟合精度 即误差, 为下次是否要修改及修改多少作准备。 如此反复修改、 拟合, 就能使整条曲线重合或基本重合。自动化反演可分为如下步骤：①构造初始模型先定性解释地下三层水平结构，即三层视电阻率的值和两层厚度的值，用数组m0来代替这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=