大地电磁测深一维正反演(附matlab代码)(共15页)

1、精选优质文档-倾情为你奉上版权所有，盗版也无所谓！一切为了财富值！嘻嘻！（Captain Hong）大地电磁测深一维正反演摘 要 本文推导了大地电磁测深的理论计算表达式，并以水平层状介质为例，利用推导的正演计算式在MATLAB软件平台上进行正演，比较了不同层介质参数的视电阻率曲线。简要介绍了阻尼最小二乘法反演的基本原理和反演迭代步骤，并对多种层介质进行了反演。关键词 大地电磁，一维正反演，阻尼最小二乘法1 引 言20世纪50年代初，苏联学者吉洪诺夫和法国学者卡尼亚的经典著作奠定了大地电磁测深法（MT）的基础。它是利用大地仲频率范围很宽（）广泛分布的天然变化的电磁场，进行深部地质构造研究的一种频

2、率域电磁测深法。由于该法不需要人工建立场源，装备轻便、成本低，且具有比人工源频率测深法更大的勘探深度，所以除主要用于研究地壳和上地幔地质构造外，也常被用来进行油气勘查、地热勘探以及地震预报等研究工作。几十年来，由于大地电磁测深法具有以下几个优点：不受高阻屏蔽，对低阻分辨率高；不用人工供电，勘探成本低且工作方便；勘探深度范围大。使大地电磁法在矿产勘探及普查、地壳岩石圈电性结构研究、海洋地球物理勘探、地热勘探、能源勘探、隐伏岩溶水结构、天然地震预测等都扮演着至关重要的角色。大地电磁也存在一些缺点，比如在实际应用的过程中整理后的数据存在分散的情况；频率范围不够宽，特别是缺少高频成分，受噪音影响大信噪

3、比低；所需观察时间长，致使野外工作效率低。随着基础理论、技术手段、仪器设备的不断完善和发展，进一步改进和解决这些问题，才能将大地电磁法更好的应用于生产服务当中。2 视电阻率及水平地层大地电磁测深曲线的理论计算方法2.1大地电磁测深理论的几点假设和论证吉洪诺夫和卡尼亚提出了假设并论证了以下几点：将场源近似地看为平面电磁波垂直入射大地。引入波阻抗的概念（Z=E/H），表征地球电性分布对大地电磁场的响应。利用单点大地电磁场观测研究地球电性分布是可能的。2.2视电阻率及水平地层上的理论计算表达式视电阻率概念是从均匀介质中电阻率和波阻抗关系引申出来的。在均匀介质中有借用这一关系式，把非均匀介质的地面波阻

4、抗代入上式，称相应的电阻率为视电阻率，用表示： （2-1）式中波阻抗的第二个脚码表示层状介质总的层数，第一个脚码表示波阻抗所在层面位置的编号，表示层介质情况下第一层顶面处的波阻抗。通常，视电阻率不是介质的真电阻率，它是介质电阻率的综合反映，并和电磁波的周期（或频率）有关，因为不同周期电磁波的穿透深度不同，当频率很高时，由于趋肤效应，电磁波只能集中在第一层介质中，电磁场不受下伏岩层电阻率的影响，这时视电阻率。随着电磁波信号周期的增大，它的穿透深度也增大，视电阻率值将受到深部介质电阻率分布的影响。显然，视电阻率和地下介质电阻率分布以及电磁波信号周期之间的函数关系，可以由地面波阻抗递推公式给出。但是

5、，我们通常用阻抗比（或称为变换函数）的递推公式来表示。定义变换函数为式中 代入后得到变换函数的递推关系： （2-2）地面的波阻抗为于是，层状介质的视电阻率公式为 （2-3）其中： （2-4）当然，式（2-2）也可以写成双曲线正切的形式，此时式（2-4）将有相应的变换。变换函数还可以用反射系数来表示，这时有 （2-5） （2-6）地球物理工作者通常把野外观测求得的不同周期的地面波阻抗，换算为视电阻率，利用随信号周期变化的视电阻率曲线研究地下介质电阻率的分布。2.3水平地层大地电磁测深曲线的理论计算方法大地电磁测深的理论曲线是指在给定地下介质电阻率分布的情况下，通过计算得出的视电阻率和信号周期之间

6、的函数曲线。当层状一维介质的地电参数，和，给定时，我们根据下列递推公式来讨论在计算机上进行计算的程序设计问题，即 其中为第层的复波数。当波长以千米为单位时，于是和也都是复数，但二者均为量纲一参数。考虑到Matlab软件平台必须把复数分解为实部和虚部在进行运算，为此，令， 求解的实部和虚部：即 （2-7） （2-8）求解的实部和虚部，并将其中复数：做如下变换，令有 （2-9）其中， （2-10） （2-11）可以求得： （2-12） （2-13）对每一个值计算时，递推计算是由下而上逐层进行的。由于，故有，。可以计算出的实部和虚部： 它可以看做是令，即取来求相应的和，接着就可以算出和，这就完成了由

7、计算再计算的一个循环计算。然后使依次递减，再做类似的运算，到并求出和计算相应的视电阻率值：3 水平地层上的正演模拟3.1二层水平地层二层断面的视电阻率函数表达式为视电阻率曲线以的数值为纵坐标，以数值为横坐标绘在双对数坐标系上。（1） G型：指型地点断面曲线，G型曲线的中部存在有的极小值。图1 G型地层视电阻率测深曲线：400,1500， ：100（2） D型：指型地点断面曲线，D型曲线的中部存在有的极大值。图2 D型地层视电阻率测深曲线：1500,400， ：1003.2三层水平地层三层断面的视电阻率函数表达式为视电阻率曲线以的数值为纵坐标，以数值为横坐标绘在双对数坐标系上。（1） H型：指型

8、地电断面的曲线。H型曲线的值先减小后增大再减小，曲线左支趋近于第一层电阻率值，曲线右支趋近于第三层电阻率值。图3 H型地层视电阻率测深曲线：400，1500，400， ：100，200（2） K型：指型地电断面的曲线。K型曲线的值先增大后减小再增大，曲线左支趋近于第一层电阻率值，曲线右支趋近于第三层电阻率值。图4 K型地层视电阻率测深曲线：1500，400，1500， ：100，200（3） A型：指型地电断面的曲线。A型曲线的值先减小再增大，曲线左支趋近于第一层电阻率值，曲线右支趋近于第三层电阻率值。图5 A型地层视电阻率测深曲线：400，1500，3000， ：100，200（4） Q型：

9、指型地电断面的曲线。Q型曲线的值先增大再减小，曲线左支趋近于第一层电阻率值，曲线右支趋近于第三层电阻率值。图6 Q型地层视电阻率测深曲线：3000，1500，400， ：100，2004 水平地层上的阻尼最小二乘法反演4.1阻尼最小二乘法原理MT的反演问题属于离散非线性反演，可利用泰勒级数在给定的初始模型处展开并线性化，考察下面的目标函数 （4-1）其中为电阻率，小标和分别表示观测量和拟合量，为实际观测数据的周期点数。把做泰勒展开，并保留一次项，把上式改写为： （4-2）其中为模型参数的个数。整理后得 （4-3）下面是对反演过程的三个具体问题的具体处理方法。（1） 偏导数的计算。可以采用差分方

10、式进行求解。取，可得 （4-4）（2） 对模型参数的处理模型参数中的电阻率和厚度有不同量纲，尤其是电阻率参数，变化范围很大，如果直接由上述方法求取，不但会导致方程（3-3）严重畸变，而且电阻率和厚度参数的修改也不会正确，从而会导致反演方法不收敛。为了解决这个问题，可对方程（3-4）无量纲化，即进行量纲化归一化，将（3-4）式的偏导数代入方程式（3-3）中，可以得到 （4-5）其中括号内参数为元函数的模型参数列表，把上式写成矩阵形式，便构成下面的方程组 （4-6）式中 （4-7） （4-8） （4-9）由上述方程组可知最小二乘法修改步长。（3） 最小二乘法步长修改中的矩阵是一个对称、正定矩阵。当

11、近乎奇异时，它有一个或数个小的特征值，从而会使得参数修正步长和实际参数差向量之间相差很大，和的方向甚至会相反。对许多地球物理问题来说，当选择的初始参数与实际值相差较大、误差较大或参数间线性相关时，就有可能出现这种情况。当矩阵具有一个或多个接近于零的特征值时，进行逆运算，小的特征值就会对反演解有相当大的影响，导致参数改变向量变得很大，从而使得方程的线性化不再是精确的，这是造成不稳定的原因。为了避免上述情况的出现，可以在矩阵中加上一项，其中是单位矩阵，。其结果是使得该矩阵中的任何一个特征值都增加了一个量，变成为，从而使得矩阵不再是奇异的，亦即使变成，这样对角项就不会变得很大，也就被“阻尼”了。其参

12、数估计值向量的表达式为： （4-10）4.2阻尼最小二乘法反演阻尼最小二乘法反演步骤：（1） 输入初始模型参数和，计算目标函数值；（2） 输入一个初始阻尼因子和缩放系数（一般取，本文给定初始阻尼因子为2）；（3） 求解模型修改量，并求取修改后的模型参数；（4） 由新的模型参数，计算新的目标拟合函数；（5） 对和的大小进行比较，如果<，则证明本次迭代是收敛的，即令阻尼因子减小为，并把修正后的模型参数作为初始模型参数重新进行迭代；如果>，则证明本次迭代是发散的，即令阻尼因子增大为重新迭代。4.3反演结果本文目标函数，初始阻尼因子，缩放系数。（1）两层G型反演结果真值：；初值：。迭代8次

13、，反演参数为（2） 三层H型反演结果真值：；初值：。迭代14次，反演参数为。（3） 四层HA型反演结果真值：初值：迭代14次，反演参数为5 结论反演结果表明，阻尼最小二乘法反演收敛速度和稳定性都较好，但对初值的选取略有要求，初值选取的不适，特别是四层介质的初值选取不适，会导致目标函数取得局部最小值，从而局部收敛偏离真值。可以先进行博斯蒂克反演，进行半定量的解释，求得第一层介质参数，从而约束阻尼最小二乘法的初值选取。参考文献1 张建. CSAMT一维全区约束反演研究与应用D. 桂林理工大学, 2011.2 李斌，郭嵩巍，郑凯，等. 基于MATLAB的大地电磁测深正反演实现以层状一维介质的正演和阻

14、尼最小二乘法反演为例J. 内蒙古石油化工. 2010, 36(9): 35-36.3 周俊杰，强建科，汤井田. 一种优化的CSAMT一维反演方法Z. 桂林: -59.4 刘崧.谱激电法.中国地质大学出版社.1988.5.程序附录：MTforward.m（正演程序）clear allro=200,500;h=1000;T=1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700,800,900,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,20000,3

15、0000,40000,50000,60000,70000,80000,90000,;ros=zeros(1,length(T);Y=zeros(1,length(T);X=zeros(1,length(T);for i=1:1:length(T)X(i)=(sqrt(10*ro(1)*T(i)/h(1);endfor i=1:1:length(T)n=length(ro);u=zeros(1,n);v=zeros(1,n);u(n)=1;v(n)=0;P=zeros(1,n);Q=zeros(1,n);A=zeros(1,n-1);B=zeros(1,n-1);G=zeros(1,n-1);f

16、or j=1:1:length(G)G(j)=(4*pi*h(j)/(sqrt(10*ro(j)*T(i);endwhile(n>1)P(n)=(1-(sqrt(ro(n)/ro(n-1)*u(n)2-(sqrt(ro(n)/ro(n-1)*v(n)2)/(1+sqrt(ro(n)/ro(n-1)*u(n)2+(sqrt(ro(n)/ro(n-1)*v(n)2);Q(n)=(-2)*sqrt(ro(n)/ro(n-1)*v(n)/(1+sqrt(ro(n)/ro(n-1)*u(n)2+(sqrt(ro(n)/ro(n-1)*v(n)2);A(n-1)=exp(-G(n-1)*(P(n)*

17、cos(G(n-1)-Q(n)*sin(G(n-1);B(n-1)=exp(-G(n-1)*(Q(n)*cos(G(n-1)-P(n)*sin(G(n-1);u(n-1)=(1-A(n-1)2-B(n-1)2)/(1+A(n-1)2+B(n-1)2);v(n-1)=(-2*B(n-1)/(1+A(n-1)2+B(n-1)2);n=n-1;endros(i)=ro(1)*(u(1)2+v(1)2);Y(i)=(u(1)2+v(1)2);endloglog(X,Y)xlabel('周期（s）')ylabel('视电阻率/ohm.m')grid on;MTinv.m（

18、反演程序）clear allro=500,100,500,1000;h=300,200,400;T=1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700,800,900,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,20000,30000,40000,50000,60000,70000,80000,90000,;Y=zeros(1,length(T);Yi=zeros(1,length(T);for i=1:1:length(T)Y(i)=MT(r

19、o,h,T(i);endX=zeros(1,length(T);for i=1:1:length(T)X(i)=(sqrt(10*ro(1)*T(i)/h(1);endro1=400,50,400,900;h1=100,100,100;for i=1:1:length(T)Yi(i)=MT(ro1,h1,T(i);endobject=0;for i=1:1:length(T)object=object+(Yi(i)-Y(i)2;endA=zeros(length(T),length(ro)+length(h);B=zeros(length(T),1);zuni=2;sb=0;while(obj

20、ect)>0.;for i=1:1:length(T)for j=1:1:length(ro)ro1(j)=ro1(j)*1.1;Y_temp=MT(ro1,h1,T(i);A(i,j)=10*(Y_temp/Yi(i)-1);ro1(j)=ro1(j)/1.1;endfor k=length(ro)+1:1:length(ro)+length(h)h1(k-length(ro)=h1(k-length(ro)*1.1;Y_temp=MT(ro1,h1,T(i);A(i,k)=10*(Y_temp/Yi(i)-1);h1(k-length(ro)=h1(k-length(ro)/1.1;

21、endI=zuni*eye(length(ro)+length(h);endfor i=1:1:length(T)B(i)=(Y(i)/Yi(i)-1;end% x=pinv(A)*B;% x=lsqnonneg(A,B);% x=AB; x=(A'*A+I)(A'*B);for j=1:1:length(ro)ro1(j)=ro1(j)+ro1(j)*x(j);endro1for i=length(ro)+1:1:length(ro)+length(h)h1(i-length(ro)=h1(i-length(ro)+h1(i-length(ro)*x(i);endh1for i=1:1:length(T)Yi(i)=MT(ro1,h1,T(i);endtemp=object;object=0;for i=1:1:length(T)object=object+(Yi(i)-Y(i)2;endif object<temp zuni=zuni/10;elseif object=temp zuni=zuni*1;elsei