还剩10页未读,继续阅读
本资源只提供10页预览,全部文档请下载后查看!喜欢就下载吧,查找使用更方便
文本内容:
一、电阻率测深法原理电阻率法是通过观测地表的电场来了解地下介质电性分布的,要探测一定深度的地层的存在,必须使其明显的影响到地表的电场分布,也就是要求其对观测点处的电场有明显的扰动,而要做到这一点,就要求流入相应深度的电流份额足够多因为在电阻率法中都是使用点电流源,因此需要考察距离点电流源不同距离的时透入某一给定深度以下供电电流所占比例的变化规律在相距的两个异2L性点电流源之间的中垂面上任意一点上的电流密度为AB_IL1」it L2+y2+z23/2式中,为观测点距连线的水平距离;为深度;为供电电流强度透入给定y ABZ I深度以下的相对电流强度为zI/I=-J°°J00dydZ.=1--arctan1-1z223/2冗一z18,L2+y2+z23/2u L下图所示为透入深度以下空间的电流随变化的情况从图中可以看出,当z lz/l L/z值比较小时,透入深度以下空间的电流比例也小,只能探测到近地表的L/z zIz/I情况;增大电测深的供电电极距时,透入某一给定深度以下的供电电流比例L z将随之增大当较大时,就可以探测到较深的部位L/z在研究地下介质电阻率的垂向变化时,希望尽量减小横向电阻率变化的影响如前所述,移动测量电极对地下介质电阻率的横向变化反映非常明显,而移MN动供电电极对地下介质电阻率的横向变化反映则远没有那么明显为了减少横AB向电阻率变化的影响,应该采用一种测量电极基本保持不动,主要移动供电MN电极的装置在实际工作中,一般采用对称四极测深装置,在施工条件限制时,AB也可采用三极测深装置,其他装置则很少使用
二、对称四极测深装置简介对称四极测深装置野外工作布置如下图所示,供电电极和测量电极AB MN都以测点为中心对称布置在一条直线上最初的供电电极距仅数米,逐步取一系列的递增值,每个数量级距离供电极距改变约次,各供电极距在对数5-6AB/2轴上应均匀分布大致按照相同的倍数增大每一个供电极距与前一个供电极距的比值大约为左右选择供电极距时,要求最小的极
1.2—
1.5Ps=O;for kk=l:20jj=kk-ii1-l:kk-ii18;m=
0.11396*
10.-jj/6;for k=nT:-1:1T=Pk*T+Pk+T-Pk.*exp-2*m*h k./T+P k-T-Pk.*exp-2*m*h k;endPs=Ps+T*C kk;end嬲水平层状大地模型电阻率测深正演嬲load modelfile.txt;M=modelfile;n=Ml;P=M2:n+1;h=Mn+2:2*n;ii=l:18;Ps=myfunctionn,P,h;r=
10.ii/6;loglog r,Ps,
5.J;xlabelAB/2;ylabelPs;save Ps.mat;运行程序所得电测深曲线如下如图所示,为取了个不同的供电极距所得的离散化的视电阻率图形18AB/2)反演2用上面已知模型所包含的模型参数[]经过正演modelf ile.txt35010100510所得的数据(保存在中)作为反演过程中已知的观测视电阻率数据,Ps.mat进行自动反演运算,反演过程中设定的初始模型保存在文件中,modelfilel.txt用行向量表示为[]第项为地层层数,第、、项分别为各层电3333331234阻率,第、项分别为各层厚度编程如下56MATLAB概三层大地模型电测深自动反演概load modelfilel.txt loadPs.matr=
10.ii/6;loglog r,Ps,‘.;xlabelAB/2;ylabelPs;hold on;Ml=modelfilel;nl=Ml1;P1=M12:nl+1;hl=Mlnl+2:2*nl;Psl=myfunctionnl,Pl,hl;err=sqrtsumPs-Psl./Ps.^2/18;while err
0.05for j=l:2*nlTmm=Ml2:2*nl;Pl=mml:nl;hl=mm nl+1:2*nl-l;Ps2=myfunction nl,Pl,hl;for i=l:18求取系数矩阵ai,j=10*Ps2i/Psl i-l;%endendb=Ps./PsiT;[U,W,V]=svda;x=V*pinvW*U*b;for jj=l:2*nl-lif xjj
1.5x jj=l.5;elseif xjj-
0.5xjj=-O.5;else endendmm=l+x.*inni;Ml=[nl,mm];Pl=mm1:nl;hl=mm nl+1:2*nlT;Psl=myfunction nl,Pl,hl;err=sqrtsumPs-Psl./Psi.2/18;endM2=[nl,mm];disp M2n2=M2l;P2=M22:nl+1;h2=M2nl+2:2*nl;Ps2=myfunction n2,P2,h2;r=
10.ii/6;loglog r,Ps2;xlabelAB/2;ylabelPs2;hold off;运行所得图形如下:演结果模型为[
3.
000049.
994710.
546598.
96134.91模型行向量中第项为地层层数,第、、项分
9910.9610]o1234别为各地层电阻率,第、项分别为各地层厚度对比可发现,自动反演所得56到的模型误差在可接受范围内,效果比较好距应能反映地表浅层电阻率,最大的极距则能满足勘探深度要求,并保证测深曲线尾支的完整,不妨碍解释最后一个电性层从勘探深度方面考虑,供电电极距应从最小勘探深度的一半到最大勘探深度的倍左右测量电极开始是固AB/25MN定的,例如取
0.5m;直到随着供电电极距的加大电压过小时,才去另一增大值,例如以此类推,一般的大小大约为的在改变时一般3m,MN AB1/3—1/30MN要求有个供电极距以组极距观测因为增大测量电极距会降低勘探22MN MN深度,因此增大测量电极距时,曲线通常会出现脱节现象另外还有一种特殊ps的对称四极装置,它是始终保持称为装置,西方国家用的较多,MN=AB/3,Wenner这种装置的曲线是光滑的,没有脱节问题ps
4.1多层R水=平Vr2地4-z层2^0T2TIRd1J11z=0=03-3dz.-Pi在地面处电流密度法向分量为零,即2除场源点外,电位处处有限,且无穷远处电位为零,即3U…703-4i=L2,……,n-1Uiz=Hi=Ui+iz=Hi在岩层分界面处电流密度法向分量连续,即5工警一竽i=i,2,……n Pidz z=Hi Pi+1dz z=Hj3-5式可用分离变量法求解,设3-1U=RrZz3-6经分离变量后得到两个二阶常微分方程d2R.1°R,2n n+--+mR=0r dr3-7dr2d2Z-m2Z=03-8dz2在岩层分界面处电位连续,即4其中,式的解为第一类和第二类零阶贝塞尔函数和3-7Jomr Y0mr第二类零阶贝塞尔函数在的轴上趋于无穷,这种特征不符Yomr r=0Z合点场源特征,因此应舍去,这样式的解为第一类零阶贝塞尔函数3-7Jomr式的解为3-8Zz=Ae-mz+Bemz3-9由此可写出各层电位积分形式的通解Uir,z=+B memz]J mrdm3-10i0式中,和为待定函数;Am Bimi=l,2,……,n在第一层介质中的电位,还应附加电源点电位,即Uir,z=黑+5J;[Aime-mz Bmemz]J mrdm3-11+10z SlZI f°°+memz]mj mrdmodz2nr2+z23/2ZTlKZTlu据边界条件式可得3-3,Aim=Bim又根据李普希积分:3-12Uir,z=+Bim]e-mz_|_Bimemz}Jomrdm3-13在第n层介质中么,当ZT8时,电位有限,因此Bnm=O,意=e_m|z|J mrdmK u0Unr,z Ame-mzJ rnrdm3-14n0Uir=TimJomrdm3-15代入边界条件后可得到地面电位的表达式式中,T1为电阻率转换函数注意到Tn=Pn,可得到T1的递推公式为Tn=Pn「Ti+i+pi+Ti+PieTmhiT=3-16一i PiT+-T-e-2mhii+1Pi i+1Pi水平地层上视电阻率表达式及滤波算2式对「微分,得电场强度3-15Er=~3-17dr由此得到时的对称四极装置视电阻率表达式f0=产Ps3-18其中一是供电极距AB/2在电测深视电阻率的表达式中的被积函数可以分为两部分,一是包含地下各层所有信息厚度及电阻率的电阻率转换函数7;加,二是与地层参数无关的贝塞尔函数,虽然没有解析计算结果,但可用线性滤波方法进行计算根据电阻率转换函数的变化规律,对加的抽样取对数等间隔比较合适,因此,首先令P«=J/fe2M3-19,心ex=mr,则电阻率的表达式变为根据采样定理,一个函数可以用它在等间距离散抽样点上的抽样值表达00sin[ix-左Ax/Ax]/x=z3-20〃k=-g x-ZAx/Ax将电阻率装换函数用它在数轴上的离散抽样值表达为x1乙T sin[»x-ZAx/Ax]3-211r占21r〃x—k Ax/Ax3-228kN3-23将预先计算出来,实际上取有限项求和就可以达到足够的精度了从频谱分G析的观点看,当电阻率装换函数用它在%数轴上的离散抽样值表达式时,相当于滤取了频率高于的谐波成分,因此这种计算视电阻率的方法称为滤波计算方法,1/2—称为滤波系数为了提高线性滤波计算的精度,减少滤波系数的数目,需要对%的抽样点进行位移,实际使用的线性滤波计算视电阻率的公式为心AAr+s■⑺=£z—Q3-24k=r式中,为供电极距为时的视电阻率;工/-+$为相爪+,时的电阻率转换0r r=6函数;为第左个滤波系数;为抽样间隔;为位移系数Q Axs实际线性滤波计算常用的抽样间隔有两种一种为六点是抽样间隔,即Ar=lnl0/6=101/6,直流电阻率测深曲线一般都采用这种抽样间隔进行计算另一种为十点式抽样间隔,即—=电磁测深曲线一般都采用这种抽样间隔11110/10=1”1°,进行计算下表为常用的一套六点式抽样间隔的滤波系数,共有个系数其位20k=1〜20,移系数s=-
2.1719,^21719=
0.11396o用上述公式就可计算某个供电极距r的视电阻率,只要计算个对应于这个供电极距厂的不同加值的电阻率转换函数工,将其与20mJ下表中对应的个滤波系数相乘再求和就可以了20K1234567Ck
0.003042-
0.
0011980.
012840.
02350.
086880.
23740.6194K891011121314Ck1,
18170.4248-
3.
45072.7044-
1.
13240.393-
0.1436K151617181920Ck
0.05812-
0.
025210.
011250.
0049780.002072-
0.000318一因为要计算一条视电阻率测深曲线就需要计算多个供电极距〃的视电阻率,为减少计算量,取彳==106,则s=e爪+s/.=Az+i=
0.11396x10fc-,/6,这样,计算r不同供电极距的视电阻率所需要的电阻率转换函数大多数可以共用这样取值时,r线性滤波计算视电阻率的公式可以写为22亿=Z1%TQ3-25k=]下面是用数值滤波法计算视电阻率测深曲线的计算机流程输入层参数,包括层数〃、各层的层厚度”和电阻率1q;输入要计算的供电极距范围,并由此得到中的变化范围2r=d i••1〜1•min max9计算%的变化范围:;3-i k-i,-k-i\/min\/nuix4用电阻率转换函数递推公式,循环计算吗=e加+s=o.i函数xio〃6,『Lx时的工吗卜5用滤波公式2-19,循环计算[i=imin,…,imax时的2⑺
四、自动反演法理论反演是由响应到模型的过程,是正演的逆过程,对水平层状大地对称四极电阻率测深曲线的解释就是反演问题,主要解释方法有四种特征点法、直接反演法、正演拟合法,自动反演法由于计算机的便利,我们通常采用自动反演法电阻率测深数据的自动反演方法可归结为寻找模型使=1,2,…,NM下面的目标函数0趋于极小:NS江|氏-加『
44.1i=\式中,%是第个极距的实测视电阻率,河是由模型/正演计算所得的i/%=1,2,...”第,个极距的理论视电阻率,是模型参数地层的电阻率或厚度,M/=l,2,…,NM是模型参数个数,是供电极距数,是范数,当时即为最小二乘法NM NS1a=2将模型在预测模型处展开,并忽略二次项以上的项,式表达可改为求预测模M4-1型修改量使目标函数趋于极小AMa NSNM亿「时当前第2⑷2式要趋于极小,则对于各供电极距要满足下面的线性方程4-2=1,2,…,NSNM QY.=p.-n.4-3」%£dMj解上述线性方程组,就可得到预测模型〃的修改量△〃,于是可得到新的预测模型,计算新模型的理论视电阻率,与实测视电阻率进行对比,如果精度满足要求,则新的预测模型即为反演结果;否则重新展开,计算模型修改量,直到满足要求下面是对反演过程的三个具体问题的具体处理方法偏导数的计算可用差分方法来计算,取则1出30_Pai…,
1.1M j,…,M Pai…,M j,…,MMI,M’NM-MI,M,NMdMj
0.1M;模型参数的处理模型参数中有不同量纲的电阻率和厚度,尤其对电阻率2参数,变化范围很大,如果直接由上述方法求解,不但会导致方程严重病态,4-3而且电阻率和厚度参数的修改量也不会正确,从而导致反演方法不收敛为了解决这个问题,可对方程无量纲化将式的偏导数代入方程式中,有4-34-44-3NM£10[4(必,心,・・,,1・1/,・・・,丹劫-4必,根,,・,,/,・・,,孙“)]言(4-5)J=1Mj=Psi~Pai(1,加,…,M j,…,M丽),,=L2,…,NS2上式两端同时除以4(陷,2,・・・,//,・.)得.陷,加…,监,…,L1K VM1AM.p.*142,YlO\/-----------」处-1—L=2-\,i=l,2,…,NS4-6夕陷,加…,…,Mp/MWMj Paij=l2,Pai M,m
2.IM,,…,MNM4=ioPaiM1,M,…,Mj,…,M2NM AM.x=Mj8=-1Pai这样方程组可写为矩阵形式4-6AX=B(4-7)解线性方程组4-7可得新模型参数M*(M;,1,2,…,NM)为=M+AM.=(l+x.)A/.(4-8)j7另外为了防止模型参数修改过量,实际过程中又可作如下规定勺>时取
1.5为=
1.5;为<-
0.5时取勺=-
0.5,即每次修改量不超过原有模型参数值的一半,保证收敛稳定
(3)线性方程组的求解用奇异值分解法求解,其基本思想是建立在奇异值分解定理上即任意阶矩阵均可分解为这里为阶正交NSxNM4A=LWy7,U NSxNS阵和为阶正交阵,为的奇异值组成的阶对角阵当非奇异V NMxNMW ANSxNM A时,奇异值较大,方程组有广义逆解当接近奇异时,有的奇异值3-7X=DMTU3,A就较小,此时由于卬一中系数过大,上述解的误差就较大为了解决这个问题,维根斯(Wiggins)提出用最接近A的矩阵R来代替A,R=UWy\其中也是将W中小c的奇异值用零代替,因此有较精确的广义逆解下面是计算机反演的具体流程
(1)输入实测视电阻率曲线夕s);()根据实测视电阻率曲线确定预测模型的层数〃;2
(3)根据曲线特征,初步确定预测模型的层参数7=1,2,…
(4)调用正演程序计算预测模型的视电阻率曲线々(厂);
(5)对比自⑺和金⑺曲线,计算拟合误差s;11vf(4Q)—〉——,err=J、(4-9)2⑺J若拟合误差满足精度要求时,输出预测模型参数,反演结束,若不满足精度要求时,由方程计算出模型修改量,得到新的预测模型,返回第四步,直到满足精度要3-6求
五、在三层大地模型电阻率测深正反演中的MATLAB实现)正演1已知模型参数三层水平层状模型,第一层电阻率为50,地层厚度为5;第二层电阻率为10,地层厚度为10;第三层电阻率为100o在程序中用dat.txt文件给出,文件内容()第一项表示层数,后面项为各层的电阻率、最后35010100510332项为各层厚度正演编程如下MATLAB%定义子程序外function()function Ps=myfunction n,P,h()T=P n;C=[
0.003042-o.
0011980.
012840.
02350.
086880.
23740.
61941.
18170.4248-
3.
45072.7044-
1.
13240.393-
0.
14360.05812-
0.
025210.01125-
0.
0049780.002072-
0.000318];ii=l:18;。