电法报告-直流电测深正演曲线

关于多层水平地层地电断面电测深曲线的正演的读书报告姓名：林俊班级：061084-27学号：20081003195指导老师：师学明日期：二〇一一年五月1目录前言…………………………………………………………2目的……………………………………………………………2任务要求………………………………………………………2工作过程………………………………………………………2成果……………………………………………………………2原理…………………………………………………………3§1-多层水平地层上的对称四极电测深视电阻率表示式……31.多层水平地层地面点电流源的电场……………………………32.多层水平地层上电测深的ρs表示式和电阻率转换函数……53.电阻率转换函数的递推公式……………………………………6§2-水平地层上视电阻率的滤波算法…………………………6§3-多层水平地层的电测深曲线类型…………………………9A二层情况…………………………………………………………9B三层情况…………………………………………………………9C四层及多层情况…………………………………………………9编程…………………………………………………………10感想…………………………………………………………182关于多层水平地层地电断面电测深曲线的正演的读书报告前言目的：熟悉并掌握多层水平地层地电断面直流电测深曲线的正演任务要求：编制适用于n层地电断面的正演电测深程序（编程环境不限制，可用C语言，C++，VC，VB，matlab，推荐用matlab）。每个同学计算两个标准地电模型的正演计算第一个模型：二层G型地电模型第一层地层电阻率10欧姆米，第一层厚度10米；第二层地层电阻率100欧姆米第二个模型：三层H型地电模型第一层地层电阻率：班号（4）×100欧姆米，第一层厚度15米；第二层地层电阻率：序号（27）×1欧姆米，第二层厚度20米；第三层地层电阻率：1000欧姆米AB/2为13个：2,3,4.5,6,9,12,15,20,30,45,60,90,120（米）。工作过程：先进行原理分析，再用matlab进行编程，最后小结。成果：用matlab实现了n层地电断面的直流电测深正演。3原理§-多层水平地层上的对称四极电测深视电阻率表示式1.多层水平地层地面点电流源的电场如图所示，水平地面下有层水平地层，各层电阻率分别为ρ1、ρ2…ρn；各层厚度分别为h1、h2…hn-1；各层底面到地表的距离分别为H1、H2…Hn-1，Hn→∞。求解思路：U→E→ρs→Ti(λr)①用分离变量法求方程②边界条件③通解④各层的电位表达式⑤利用边界条件求待定系数A2、A3，…,AnB1、B2，…,Bn-1设地面点电流源A的强度为I。为求各层中的电位表达式，将柱坐标系的原点设在A点，Z轴垂直向下。在所设条件下，电位与角无关，满足如下形式的拉普拉斯方程：022122zUrUrrU4边界条件：（1）0|221rzU（2）0|01zzU（3）RCRIUrz11012|22（22rzR）（4）IIHziHziUU||111（5）IIHziiHzizUzU|1|11111i由分离变量)()(),(zZrRzrU得到零阶贝赛尔方程，其通解为：001)()()(dmmrJemBemAUmzmz式中：A(m),B(m)为待定的积分变量m的函数；J0(mr)为零阶贝赛尔函数。利用边界条件，可以得到第一层电位公式：00111)())((22dmmrJeemBeIUmzmzmz第二层以下至n-1层，第i层的电位为：00)()()(dmmrJemBemAUmzimzii第n层内的电位表达式，由0zU得：()100()()nmznUAmeJmrdm电测深只在地面工作，即z＝0，故只需求出B1，111002()()2IUBmrJmrdm式中：J０(mr)为零阶第一类贝赛尔函数；B１(m)为积分变量m的函数。5对于层数确定的水平地层，根据地层界面上电位和电流密度法向分量连续的边界条件，可具体求出B1(m)的表示式。例如，最简单的二层水平地层，利用ρ1和ρ2岩层分界面的相应边界条件可具体求出112122121)2(112)(mhmheKeKImB①n=3时，112112222((3)1223122()212231223)()1mhmhhmhmhhmhKeKeBmKeKeKKe其中232323K2.多层水平地层上电测深的ρs表示式和电阻率转换函数令)(2)(11mBImB②则地面上电位公式为：dmmrJmBIU)()(212001若采用→0的装置测量，相应的ρs表达式为：rUIrIErS2222mdmmrJmBrS)()(211021令)](21[)(11mBmT③则多层水平地层上的电测深ρs公式简写成：mdmmrJmTrS)()(1012式中，T１(m)定义为电阻率转换函数又称核函数。可见，电阻率转换函数与各层的层参数(厚度和电阻率)及积分量m有关。63.电阻率转换函数的递推公式对于二层水平地层情况，若将①式先后代入②式和③式，便得到二层水平地层的电阻率转换函数:112122121)2(111)(mhmheKeKmT归纳每一层的电阻率转换函数，就可导出电阻率转换函数的递推公式：)1)(()1()1)(()1()(212212iiiimhimhimhimhiiiemTeemTemT1i1()()()1()()/iiiiiithmhTmTmTmthmh④nnmT)(式中11)(22iimhmhieemhth为双曲函数i1()()()1()()/iiiiiiTmthmhTmTmthmh由此可知：当给定n，ρ１，ρ２．．．ρn-1,ρn和h1，h2，．．．hn-1。可递推出上层或下层的电阻率转换函数。电阻率转换函数递推公式④的导出，免去应用边界条件解方程组求系数B１(m)的计算。§-2水平地层上视电阻率的滤波算法在电测深视电阻率的表达式中的被积函数可以分为两部分：一是包含地下各层所有信息（厚度及电阻率）的电阻率转换函数T1（m）；此外是与地层参数无关的贝塞尔函数。其虽没有解析计算结果，但可用线性滤波方法进行计算。根据电阻率转换函数的变化规律，对m的抽样取对数等间隔比较合适，因此，首先令ex=mr，则电阻率的表达式变为dxeJereTrxxxs)()()(1217根据采样定理，一个函数可以用它在等间距离散抽样点上的抽样值表达：kxxkxxxkxxkfxf/)(]/)(sin[)()(将电阻率装换函数用它在x数轴上的离散抽样值表达为：kxkxxxkxxxkxreTreT/)(]/)(sin[)()(11记dxeJexxkxxxkxCxxk)(/)(]/)(sin[12则kkxksCreTr)()(1将Ck=预先计算出来，实际上取有限项求和就可以达到足够的精度了。从频谱分析的观点看，当电阻率装换函数用它在x数轴上的离散抽样值表达式时，相当于滤取了频率高于1/2Δx的谐波成分，因此这种计算视电阻率的方法称为滤波计算方法，Ck称为滤波系数。为了提高线性滤波计算的精度，减少滤波系数的数目，需要对x的抽样点进行位移，实际使用的线性滤波计算视电阻率的公式为21)()(1kkkksxksCreTr式中，ρs（r）为供电极距为r时的视电阻率；T1(ekΔx+s/r)为mk=ekΔx+s/r时的电阻率转换函数；Ck为第k个滤波系数；Δx为抽样间隔；s为位移系数。实际线性滤波计算常用的抽样间隔有两种。一种为六点是抽样间隔，即Δx=（ln10）/6=101/6，直直流电阻率测深曲线一般都采用这种抽样间隔进行计算。另一种为十点式抽样间隔，即Δx=（ln10）/10=101/10，电磁测深曲线一般都采用这种抽样间隔进行计算。8下表为常用的一套六点式抽样间隔的滤波系数，共有20个系数：k=1~20，其位移系数s=-2.1719，e-2.1719=0.11396。用上述公式就可计算某个供电极距r的视电阻率，只要计算20个对应于这个供电极距r的不同m值的电阻率转换函数T1（mk），将其与下表中对应的20个滤波系数相乘再求和就可以了。因为要计算一条视电阻率测深曲线就需要计算多个供电极距r的视电阻率，为了减少计算工作量，取ri=eiΔx=10k/6，对应的mk=ekΔx+s/ri=e(k-i)Δx+s，这样，计算不同供电极距的视电阻率所需要的电阻率转换函数大多数可以共用。下面是用数值滤波法计算视电阻率测深曲线的计算机流程：（1）输入层参数，包括层数n、各层的层厚度和电阻率；（2）输入要计算的供电极距范围，并由此得到r=eiΔx中i的变化范围：imin~imax；（3）计算k-i的变化范围：（k-i）min~（k-i）max；（4）用电阻率转换函数递推公式，循环计算mj=ejΔx+s，j=（k-i）min，…，（k-i）max时的T1（mj）；（5）用滤波公式，循环计算ri=eiΔ，i=imin，…，imax是的ρs（ri）。9§-3多层水平地层的电测深曲线类型A二层情况二层地点断面上层电阻率为ρ１，下层电阻率为ρ２。有如下两种电测深曲线类型：（如右图）（1）D型：ρ１ρ２，基底为低阻；（2）G型：ρ１ρ２，基底为高阻。B三层情况三层断面共有以下4种类型：（如右图）（1）A型：ρ1ρ2ρ3，电阻率递增；（2）K型：ρ1ρ2ρ3，中间高阻层；（3）H型：ρ1ρ2ρ3，中间低阻层；（4）Q型：ρ1ρ2ρ3，电阻率递减。C四层及多层情况四层及多层断面电阻率测深曲线类型用三层断面类型的组合表示。四层断面及其电阻率测深曲线类型共有8种，例如ρ1ρ2ρ3ρ4的四层断面及其电阻率测深曲线类型为HA型。五层断面及其电阻率测深曲线类型共有16中，例如ρ1ρ2ρ3ρ4ρ5的五层断面及其电阻率测深曲线类型为KHK型。每多一层，曲线类型增加一倍。N层地电断面的电阻率测深曲线类型数为2n-1。10编程要求计算两个标准地电模型的正演：第一个模型：二层G型地电模型ρ1=10Ω·mh1=10mρ2=100Ω·m第二个模型：三层H型地电模型ρ1=班号（4）×100=400Ω·mh1=15mρ2=序号（27）×1=27Ω·mh2=20mρ3=1000Ω·mAB/2为13个：2,3,4.5,6,9,12,15,20,30,45,60,90,120程序设计流程图11程序如下：function[ps]=dcszhengyan(ab2,n,p,h)c=[0.003042-0.0011980.012840.02350.086880.23740.61941.18170.4248-3.45072.7044-1.13240.393-0.14360.05812-0.02520.01125-0.0049780.002072-0.000318];%滤波系数C，共20个rn=length(ab2);%极距的个数d=log(10)/6;%取样间距s=-2.1719;%位移系数ps=zeros(rn,1);%视电阻率fori=1:rnfork=1:20m=exp(k*d+s)/ab2(i);T2=p(n);%电阻率转换函数forj=n:-1:2T1=p(j-1)*(p(j-1)*tanh(m*h(j-1))+T2)/(p(j-1)+T2*tanh(m*h(j-1)));%电阻率转换函数T2=T1;12endps(i)=ps(i)+T1*c(k);endendloglog(ab2,ps,'