一维电测深正反演程序

电阻率测深法原理电阻率法是通过观测地表的电场来了解地下介质电性分布的，要探测一定深度的地层的存在，必须使其明显的影响到地表的电场分布，也就是要求其对观测点处的电场有明显的扰动，而要做到这一点，就要求流入相应深度的电流份额足够多。因为在电阻率法中都是使用点电流源，因此需要考察距离点电流源不同距离的时透入某一给定深度以下供电电流所占比例的变化规律。在相距2L的两个异性点电流源AB之间的中垂面上任意一点上的电流密度为：J=式中，y为观测点距AB连线的水平距离；z为深度；I为供电电流强度。透入给定深度z以下的相对电流强度为：IzI=Lπ下图所示为透入深度z以下空间的电流Iz/I随L/z变化的情况。从图中可以看出，当L/z值比较小时，透入深度z以下空间的电流Iz/I比例也小，只能探测到近地表的情况；增大电测深的供电电极距L时，透入某一给定深度z以下的供电电流比例将随之增大。当L/z较大时，就可以探测到较深的部位。在研究地下介质电阻率的垂向变化时，希望尽量减小横向电阻率变化的影响。如前所述，移动测量电极MN对地下介质电阻率的横向变化反映非常明显，而移动供电电极AB对地下介质电阻率的横向变化反映则远没有那么明显。为了减少横向电阻率变化的影响，应该采用一种测量电极MN基本保持不动，主要移动供电电极AB的装置。在实际工作中，一般采用对称四极测深装置，在施工条件限制时，也可采用三极测深装置，其他装置则很少使用。对称四极测深装置简介对称四极测深装置野外工作布置如下图所示，供电电极AB和测量电极MN都以测点O为中心对称布置在一条直线上。最初的供电电极距仅数米，逐步取一系列的递增值，每个数量级距离供电极距改变约5—6次，各供电极距AB/2在对数轴上应均匀分布（大致按照相同的倍数增大）。每一个供电极距与前一个供电极距的比值大约为1.2—1.5左右。选择供电极距时，要求最小的极距应能反映地表浅层电阻率，最大的极距则能满足勘探深度要求，并保证测深曲线尾支的完整，不妨碍解释最后一个电性层。从勘探深度方面考虑，供电电极距AB/2应从最小勘探深度的一半到最大勘探深度的5倍左右。测量电极MN开始是固定的，例如取0.5m；直到（随着供电电极距的加大）电压过小时，才去另一增大值，例如3m，以此类推，一般MN的大小大约为AB的1/3—1/30。在改变MN时一般要求有2个供电极距以2组MN极距观测。因为增大测量电极距MN会降低勘探深度，因此增大测量电极距时，ρs曲线通常会出现脱节现象。另外还有一种特殊的对称四极装置，它是始终保持MN=AB/3，称为Wenner装置，西方国家用的较多，这种装置的ρs曲线是光滑的，没有脱节问题。多层水平地层上视电阻率表达式水平地层上地面点电流源的电场如下图所示，假定地面水平,在地下有n层水平层状地层，各层的电阻率分别为ρ1，ρ2，……ρn；厚度分别为h1,h2,……hn，其中hn→∞；每层地面到地面的距离为H1，H2，……Hn,其中Hn→∞。在A点有一电流源供电，电流强度为I。用柱坐标系，将原点设在A点，z轴垂直向下。由于问题的解具有轴对称性，与ψ无关,因此电位分布满足以下形式的拉普拉斯方程∂2U∂r2+以及如下边界条件电源点附近，趋于地面点电流源的正常电位，即U1R=r2+z2在地面处电流密度法向分量为零，即1ρ1∂U1∂zz=0（3）除场源点外，电位处处有限，且无穷远处电位为零，即U在岩层分界面处电位连续，即Uiz=Hi=U在岩层分界面处电流密度法向分量连续，即1ρi∂Ui∂zz=式（3-1）可用分离变量法求解，设U=RrZz（经分离变量后得到两个二阶常微分方程∂2R∂r2+1∂2Z∂z2-其中，式（3-7）的解为第一类和第二类零阶贝塞尔函数J0(mr)和Y0(mr)。第二类零阶贝塞尔函数Y0(mr)在r=0的Z轴上趋于无穷，这种特征不符合点场源特征，因此应舍去，这样式（3-7）的解为第一类零阶贝塞Zz=Ae-mz+Bemz由此可写出各层电位积分形式的通解Uir,z=I2π式中，Ai(m)和Bi(m)为待定函数；i在第一层介质中的电位，还应附加电源点电位，即U1r,z=Iρ由∂据边界条件式（3-3），可得A又根据李普希积分：1R=0∞e-m得U1r,z=I2π0在第n层介质中么，当z→∞时，电位有限，因此BnmUnr,z=I2π0∞代入边界条件后可得到地面电位的表达式U1r=I2π式中，T1为电阻率转换函数。注意到Tn=ρTTi=ρiT2）水平地层上视电阻率表达式及滤波算法式（3-15）对微分，得电场强度（3-17）由此得到时的对称四极装置视电阻率表达式（3-18）其中是供电极距（AB/2）在电测深视电阻率的表达式中的被积函数可以分为两部分，一是包含地下各层所有信息（厚度及电阻率）的电阻率转换函数，二是与地层参数无关的贝塞尔函数，虽然没有解析计算结果，但可用线性滤波方法进行计算。根据电阻率转换函数的变化规律，对的抽样取对数等间隔比较合适，因此，首先令，则电阻率的表达式变为（3-19）根据采样定理，一个函数可以用它在等间距离散抽样点上的抽样值表达（3-20）将电阻率装换函数用它在数轴上的离散抽样值表达为（3-21）记（3-22）则（3-23）将预先计算出来，实际上取有限项求和就可以达到足够的精度了。从频谱分析的观点看，当电阻率装换函数用它在数轴上的离散抽样值表达式时，相当于滤取了频率高于的谐波成分，因此这种计算视电阻率的方法称为滤波计算方法，称为滤波系数。为了提高线性滤波计算的精度，减少滤波系数的数目，需要对的抽样点进行位移，实际使用的线性滤波计算视电阻率的公式为（3-24）式中，为供电极距为时的视电阻率；为时的电阻率转换函数；为第个滤波系数；为抽样间隔；为位移系数。实际线性滤波计算常用的抽样间隔有两种。一种为六点是抽样间隔，即，直流电阻率测深曲线一般都采用这种抽样间隔进行计算。另一种为十点式抽样间隔，即，电磁测深曲线一般都采用这种抽样间隔进行计算。下表为常用的一套六点式抽样间隔的滤波系数，共有20个系数：，其位移系数，。用上述公式就可计算某个供电极距的视电阻率，只要计算20个对应于这个供电极距的不同值的电阻率转换函数，将其与下表中对应的20个滤波系数相乘再求和就可以了。因为要计算一条视电阻率测深曲线就需要计算多个供电极距的视电阻率，为减少计算量，取，则，这样，计算不同供电极距的视电阻率所需要的电阻率转换函数大多数可以共用。这样取值时，线性滤波计算视电阻率的公式可以写为（3-25）下面是用数值滤波法计算视电阻率测深曲线的计算机流程：输入层参数，包括层数、各层的层厚度和电阻率；输入要计算的供电极距范围，并由此得到中的变化范围：；计算的变化范围：；用电阻率转换函数递推公式，循环计算，时的；（5）用滤波公式2-19，循环计算，时的。四、自动反演法理论反演是由响应到模型的过程，是正演的逆过程，对水平层状大地对称四极电阻率测深曲线的解释就是反演问题，主要解释方法有四种：特征点法、直接反演法、正演拟合法，自动反演法。由于计算机的便利，我们通常采用自动反演法。电阻率测深数据的自动反演方法可归结为寻找模型使下面的目标函数趋于极小：（4-1）式中，是第个极距的实测视电阻率，是由模型正演计算所得的第个极距的理论视电阻率，是模型参数（地层的电阻率或厚度），是模型参数个数，是供电极距数，是范数，当时即为最小二乘法。将模型在预测模型处展开，并忽略二次项以上的项，式4-1表达可改为求预测模型修改量使目标函数趋于极小：（4-2）式4-2要趋于极小，则对于各供电极距要满足下面的线性方程：（4-3）解上述线性方程组，就可得到预测模型的修改量，于是可得到新的预测模型，计算新模型的理论视电阻率，与实测视电阻率进行对比，如果精度满足要求，则新的预测模型即为反演结果；否则重新展开，计算模型修改量，直到满足要求。下面是对反演过程的三个具体问题的具体处理方法。（1）偏导数的计算。可用差分方法来计算，取，则（4-4）（2）模型参数的处理。模型参数中有不同量纲的电阻率和厚度，尤其对电阻率参数，变化范围很大，如果直接由上述方法求解，不但会导致方程4-3严重病态，而且电阻率和厚度参数的修改量也不会正确，从而导致反演方法不收敛。为了解决这个问题，可对方程4-3无量纲化。将式4-4的偏导数代入方程式4-3中，有（4-5）上式两端同时除以得（4-6）令这样方程组4-6可写为矩阵形式（4-7）解线性方程组4-7可得新模型参数为（4-8）另外为了防止模型参数修改过量，实际过程中又可作如下规定：时取；时取，即每次修改量不超过原有模型参数值的一半，保证收敛稳定。（3）线性方程组的求解。用奇异值分解法求解，其基本思想是建立在奇异值分解定理上即任意阶矩阵均可分解为，这里为阶正交阵和为阶正交阵，为的奇异值组成的阶对角阵。当非奇异时，奇异值较大，方程组3-7有广义逆解，当接近奇异时，有的奇异值就较小，此时由于中系数过大，上述解的误差就较大。为了解决这个问题，维根斯（Wiggins）提出用最接近的矩阵来代替，，其中是将中小的奇异值用零代替，因此有较精确的广义逆解。下面是计算机反演的具体流程：输入实测视电阻率曲线；根据实测视电阻率曲线确定预测模型的层数；根据曲线特征，初步确定预测模型的层参数；调用正演程序计算预测模型的视电阻率曲线；对比和曲线，计算拟合误差；（4-9）若拟合误差满足精度要求时，输出预测模型参数，反演结束，若不满足精度要求时，由方程3-6计算出模型修改量，得到新的预测模型，返回第四步，直到满足精度要求。MATLAB在三层大地模型电阻率测深正反演中的实现1）正演已知模型参数：三层水平层状模型，第一层电阻率为50，地层厚度为5；第二层电阻率为10，地层厚度为10；第三层电阻率为100。在程序中用dat.txt文件给出，文件内容（35010100510）第一项表示3层数，后面3项为各层的电阻率、最后2项为各层厚度。MATLAB正演编程如下：%定义function子程序%functionPs=myfunction(n,P,h)T=P(n);C=[0.003042-0.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=1:18;Ps=0;forkk=1:20jj=(kk-ii(1)):-1:(kk-ii(18));m=0.11396*10.^(jj/6);fork=(n-1):-1:1T=P(k)*((T+P(k))+(T-P(k)).*exp(-2*m*h(k)))./((T+P(k))-(T-P(k)).*exp(-2*m*h(k)));endPs=Ps+T*C(kk);end%%水平层状大地模型电阻率测深正演%%loadmodelfile.txt;M=modelfile;n=M(1);P=M(2:(n+1));h=M((n+2):(2*n));ii=1:18;Ps=myfunction(n,P,h);r=10.^(ii/6);loglog(r,Ps,'.');xlabel('AB/2');ylabel('Ps');savePs.mat;运行程序所得电测深曲线如下：如图所示，为取了18个不同的供电极距AB/2所得的离散化的视电阻率图形。反演用上面已知模型modelfile.txt所包含的模型参数[35010100510]经过正演所得的数据（保存在Ps.mat中）作为反演过程中已知的观测视电阻率数据，进行自动反演运算，反演过程中设定的初始模型保存在modelfile1.txt文件中，用行向量表示为[333333]第1项为地层层数，第2、3、4项分别为各层电阻率，第5、6项分别为各层厚度。MATLAB编程如下：%%三层大地模型电测深自动反演%%loadmodelfile1.txtloadPs.matr=10.^(ii/6);loglog(r,Ps,'.');xlabel('AB/2');ylabel('Ps');holdon;M1=modelfile1;n1=M1(1);P1=M1(2:(n1+1));h1=M1((n1+2):(2*n1));Ps1=myfunction(n1,P1,h1);err=sqrt(sum(((Ps-Ps1)./Ps).^2)/18);whileerr>0.05forj=1:2*n1-1mm=M1(2:2*n1);mm(j)=1.1*mm(j);P1=mm(1:n1);h1=mm((n1+1):(2*n1-1));Ps2=myfunction(n1,P1,h1);fori=1:18a(i,j)=10*(Ps2(i)/Ps1(i)-1);%求取系数矩阵endendb=(Ps./Ps1)'-1;[U,W,V]=svd(a);x=V*pinv(W)*U'*b;forjj=1:2*n1-1ifx