MATLAB实现直流电测深一维正演模拟

5星 · 超过95%的资源 需积分: 44 101 下载量 122 浏览量 更新于2024-10-10 5 收藏 592KB DOC 举报
"直流电测深一维正演程序设计的MATLAB实现" 直流电测深一维正演程序是用于模拟地下电阻率分布的计算工具,主要应用于地球物理勘探领域。通过这个程序,我们可以理解和分析地下的电性结构。本程序基于线性滤波法,利用MATLAB编写,旨在帮助用户掌握核函数的概念及其在电测深正演计算中的应用。 实验的核心在于将视电阻率计算转化为褶积积分问题。对称四极(MN→0)的视电阻率褶积积分公式是电测深正演的基础,通过变量替换,可以得到适合数值计算的形式。具体来说,通过引入对数变量,将原始公式转换为对数变量后的视电阻率表达式,这表达式是一个褶积运算积分形式,适于用空间域数字滤波方法处理。 正演计算模型的关键是正演滤波系数,这些系数决定了滤波过程的结果。实验中,通过递推公式计算核函数值,进而求得每个采样点的视电阻率。递推公式涉及了各层电阻率、厚度以及空间频率的计算,确保了核函数的精确模拟。 实验内容部分,用户需要输入电性层数、各层电阻率和厚度,然后程序会按照设定的采样间隔计算核函数曲线。MATLAB代码展示了如何实现这一过程,包括输入参数、执行递推公式和绘制核函数曲线。例如,当输入N=2,p=[15, 100],h=[20]时,程序会输出相应的核函数值并绘制成图。 通过这样的实验,学习者不仅可以理解核函数的数学意义,还能掌握线性滤波法在实际问题中的应用。此外,MATLAB的编程实践有助于提升编程技能,加深对电测深理论的理解。对于地球物理学家和地质工程师而言,这样的程序设计和应用能力是至关重要的,它能帮助他们在实际工作中解决复杂的问题,如探测地下矿藏、地下水分布等。
2019-06-29 上传
1dforward.for ! ldforward.for主程序 real ps(30),r(20),p(30),h(9),sd(9) ! ps(30)电测深视电阻率;r(30)极距;p(10)层电阻率;h(9)层厚度,sd(9)深度 character*50 char,name1,name2 integer nalyer,max write(*,*)"请输入模型文件名:" read(*,*)name1 open(2,file=name1) write(*,*)"请输入视电阻率文件名:" read(*,*)name2 open(1,file=name2) read(2,*) char read(2,*)nlayer read(2,*)char do i=1,nlayer read(2,*)p(i) enddo read(2,*)char do i=1,nlayer-1 read(2,*)h(i) enddo !!!!!!!!!!!!!!!!!!!!!!!计算电测深极距共20个,每个数量级5个极距 max=20 do i=1,max r(i)=10.* *(i/5.) enddo !调用一维正演子程序 zyld() call zyld(p,h,nlayer,r ,max,ps) !输出电测深曲线到文件 do i=1,max write(1,*)r(i),ps(i) enddo !!!!!!!!!!!!!!!画模型 sd(1)=0. do i=2,nlayer sd(i)=sd(i-1)+h(i-1) enddo sd(nlayer+1)=10000. open(25,file="model.dat")!画模型的文件 write(25,*)1,p(1) write(25,*)h(1),p(1) do i=2,nlayer write(25,*)sd(i),p(i) write(25,*)sd(i+1),p(i) enddo write(*,*)"模型曲线的文件为model.dat" pause end !!!!!!!!!!!一维正演子程序 subroutine zyld(p,h,nl,r,max,ps) real ps(max),t(nl),ck(20),r(max),m real p(nl),h(nl-1) integer nl,max !max极距个数;nl层数;ps(max)电测深视电阻率;r(max)极距;p(nl)层电阻率; !h(nl-1)层厚度;t(nl)电阻率转换函数;m汉克尔变换积分变量;ck(20)滤波系数; data ck/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/ t(nl)=p(nl) do i=1,max ps(i)=0. do k=1,20 m=0.11396*10.* *(k/6.)/r(i) do j=nl-1,1,-1 t(j)=p(j)*((t(j+1)+p(j)+(t(j+1)-p(j))*exp(-2.*m*h(j)))& /((t(j+1)+p(j))-(t(j+1)-p(j))*exp((-2.*m*h(j))) enddo ps(i)=ps(i)+t(1)*ck(k) enddo enddo end