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

版权所有，盗版也无所谓！一切为了财富值！嘻嘻！（Captain Hong） 大地电磁测深一维正反演 摘 要 本文推导了大地电磁测深的理论计算表达式，并以水平层状介质为例，利用推导的正演计算式在MATLAB软件平台上进行正演，比较了不同层介质参数的视电阻率曲线。简要介绍了阻尼最小二乘法反演的基本原理和反演迭代步骤，并对多种层介质进行了反演。 关键词 大地电磁，一维正反演，阻尼最小二乘法 1 引 言 20世纪50年代初，苏联学者吉洪诺夫和法国学者卡尼亚的经典著作奠定了大地电磁测深法（MT）的基础。它是利用大地仲频率范围很宽（ ）广泛分布的天然变化的电磁场，进行深部地质构造研究的一种频率域电磁测深法。由于该法不需要人工建立场源，装备轻便、成本低，且具有比人工源频率测深法更大的勘探深度，所以除主要用于研究地壳和上地幔地质构造外，也常被用来进行油气勘查、地热勘探以及地震预报等研究工作。 几十年来，由于大地电磁测深法具有以下几个优点：不受高阻屏蔽，对低阻分辨率高；不用人工供电，勘探成本低且工作方便；勘探深度范围大。使大地电磁法在矿产勘探及普查、地壳岩石圈电性结构研究、海洋地球物理勘探、地热勘探、能源勘探、隐伏岩溶水结构、天然地震预测等都扮演着至关重要的角色。大地电磁也存在一些缺点，比如在实际应用的过程中整理后的数据存在分散的情况；频率范围不够宽，特别是缺少高频成分，受噪音影响大信噪比低；所需观察时间长，致使野外工作效率低。随着基础理论、技术手段、仪器设备的不断完善和发展，进一步改进和解决这些问题，才能将大地电磁法更好的应用于生产服务当中。 2 视电阻率及水平地层大地电磁测深曲线的理论计算方法 2.1大地电磁测深理论的几点假设和论证 吉洪诺夫和卡尼亚提出了假设并论证了以下几点： 将场源近似地看为平面电磁波垂直入射大地。 引入波阻抗的概念（Z=E/H），表征地球电性分布对大地电磁场的响应。 利用单点大地电磁场观测研究地球电性分布是可能的。 2.2视电阻率及水平地层上的理论计算表达式 视电阻率概念是从均匀介质中电阻率和波阻抗关系引申出来的。在均匀介质中有 借用这一关系式，把非均匀介质的地面波阻抗代入上式，称相应的电阻率为视电阻率，用 表示： （2-1） 式中波阻抗的第二个脚码表示层状介质总的层数，第一个脚码表示波阻抗所在层面位置的编号， 表示 层介质情况下第一层顶面处的波阻抗。通常，视电阻率 不是介质的真电阻率，它是介质电阻率的综合反映，并和电磁波的周期（或频率）有关，因为不同周期电磁波的穿透深度不同，当频率很高时，由于趋肤效应，电磁波只能集中在第一层 介质中，电磁场不受下伏岩层电阻率的影响，这时视电阻率 。随着电磁波信号周期的增大，它的穿透深度也增大，视电阻率值将受到深部介质电阻率分布的影响。 显然，视电阻率和地下介质电阻率分布以及电磁波信号周期之间的函数关系，可以由地面波阻抗递推公式给出。但是，我们通常用阻抗比（或称为变换函数）的递推公式来表示。定义变换函数 为 式中                          代入后得到变换函数的递推关系： （2-2） 地面的波阻抗为 于是，层状介质的视电阻率公式为 （2-3） 其中： （2-4） 当然，式（2-2）也可以写成双曲线正切的形式，此时式（2-4）将有相应的变换。变换函数 还可以用反射系数来表示，这时有 （2-5） （2-6） 地球物理工作者通常把野外观测求得的不同周期的地面波阻抗，换算为视电阻率，利用随信号周期变化的视电阻率曲线研究地下介质电阻率的分布。 2.3水平地层大地电磁测深曲线的理论计算方法 大地电磁测深的理论曲线是指在给定地下介质电阻率分布的情况下，通过计算得出的视电阻率 和信号周期 之间的函数曲线。 当层状一维介质的地电参数 ， ，…， 和 ， ，…， 给定时，我们根据下列递推公式来讨论在计算机上进行计算的程序设计问题，即 其中 为第 层的复波数。当波长以千米为单位时， 于是 和 也都是复数，但二者均为量纲一参数。考虑到Matlab软件平台必须把复数分解为实部和虚部在进行运算，为此，令 ， 求解 的实部和虚部： 即 （2-7） （2-8） 求解 的实部和虚部，并将其中复数： 做如下变换，令 有                                  （2-9） 其中， （2-10） （2-11） 可以求得： （2-12） （2-13） 对每一个 值计算 时，递推计算是由下而上逐层进行的。由于 ，故有 ， 。可以计算出 的实部和虚部： 它可以看做是令 ，即取 来求相应的 和 ，接着就可以算出 和 ，这就完成了由 计算 再计算 的一个循环计算。然后使 依次递减，再做类似的运算，到 并求出 和 计算相应的视电阻率值： 3 水平地层上的正演模拟 3.1二层水平地层 二层断面的视电阻率函数表达式为 视电阻率曲线以 的数值 为纵坐标，以数值 为横坐标绘在双对数坐标系上。 （1）G型：指 型地点断面曲线，G型曲线的中部存在有 的极小值。 图1 G型地层视电阻率测深曲线 ：400,1500， ：100 （2）D型：指 型地点断面曲线，D型曲线的中部存在有 的极大值。 图2 D型地层视电阻率测深曲线 ：1500,400， ：100 3.2三层水平地层 三层断面的视电阻率函数表达式为 视电阻率曲线以 的数值 为纵坐标，以数值 为横坐标绘在双对数坐标系上。 （1）H型：指 型地电断面的曲线。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型：指 型地电断面的曲线。Q型曲线的 值先增大再减小，曲线左支趋近于第一层电阻率值，曲线右支趋近于第三层电阻率值。 图6 Q型地层视电阻率测深曲线 ：3000，1500，400， ：100，200 4 水平地层上的阻尼最小二乘法反演 4.1阻尼最小二乘法原理 MT的反演问题属于离散非线性反演，可利用泰勒级数在给定的初始模型处展开并线性化，考察下面的目标函数 （4-1） 其中 为电阻率，小标 和 分别表示观测量和拟合量， 为实际观测数据的周期点数。 把 做泰勒展开，并保留一次项，把上式改写为： （4-2） 其中 为模型参数的个数。 整理后得 （4-3） 下面是对反演过程的三个具体问题的具体处理方法。 （1）偏导数的计算。可以采用差分方式进行求解。取 ，可得 （4-4） （2）对模型参数的处理 模型参数中的电阻率和厚度有不同量纲，尤其是电阻率参数，变化范围很大，如果直接由上述方法求取，不但会导致方程（3-3）严重畸变，而且电阻率和厚度参数的修改也不会正确，从而会导致反演方法不收敛。为了解决这个问题，可对方程（3-4）无量纲化，即进行量纲化归一化，将（3-4）式的偏导数代入方程式（3-3）中，可以得到 （4-5） 其中括号内参数为 元函数 的模型参数列表，把上式写成矩阵形式，便构成下面的方程组 （4-6） 式中 （4-7） （4-8） （4-9） 由上述方程组可知最小二乘法修改步长 。 （3）最小二乘法步长修改中的矩阵 是一个对称、正定矩阵。当 近乎奇异时，它有一个或数个小的特征值，从而会使得参数修正步长 和实际参数差向量 之间相差很大， 和 的方向甚至会相反。对许多地球物理问题来说，当选择的初始参数与实际值相差较大、误差较大或参数间线性相关时，就有可能出现这种情况。当矩阵 具有一个或多个接近于零的特征值时，进行逆运算，小的特征值就会对反演解有相当大的影响，导致参数改变向量 变得很大，从而使得方程的线性化不再是精确的，这是造成不稳定的原因。为了避免上述情况的出现，可以在矩阵 中加上一项 ，其中 是单位矩阵， 。其结果是使得该矩阵中的任何一个特征值 都增加了一个量 ，变成为 ，从而使得矩阵不再是奇异的，亦即使 变成 ，这样对角项就不会变得很大，也就被“阻尼”了。其参数估计值向量的表达式为：