梁宏达
【摘 要】大地电磁作为一种重要的地球物理勘探方法,近些年来被广泛的地应用于勘察实践中.在野外大地电磁勘察资料解释的重要手段是大地电磁反演,反演效果的好坏将直接影响解释结果的正确与否.本文研究对比4种目前常用的大地电磁反演方法,即BOSTICK反演、RRI反演、OCCAM反演和NLCG反演,分析4种反演方法各自的优缺点,并对大地电磁正演模型数据进行反演,分析4种反演方法的反演效果,直观上对反演结果进行评价对比. 【期刊名称】《工程地球物理学报》
【年(卷),期】2012(009)005
【总页数】7页(P537-543)
【关键词】大地电磁;反演方法;分析;对比
【作 者】梁宏达
【作者单位】中南大学地球科学与信息物理学院,湖南长沙410083
【正文语种】中 文
【中图分类】P631.3
1 引 言
大地电磁测深法MT(Magnetotelluric)自20世纪50年代初由A.N.Tikhonov和L.Cagniard分别提出以来[1，2]，由于其勘探成本低廉,具有较大的勘探深度并且不受高阻层屏蔽影响, 对低阻层也有较高的分辨能力，目前已得到广泛的应用，成为矿区矿、地震预报、油气田普查勘探、岩石圈结构探测、地下水和地热资源寻、海洋地质调查等的重要手段[3～9]。
MT的反演研究一直是MT方法研究的一个重要部分，目前大地电磁常用的反演方法主要有Bostick(Bostick，1977)、Occam(Constable等人,1987, 1990)、RRI(Smith J T,Booker J R,1991)、贝叶斯统计(Spichak V,Menville M,1995)、零空间反演法(Michael M D,Guust,1996)、SBI(Smith T,Hoversten S C,1999)、NLCG(Newman G A , Alumbaugh D L,2000)等。随着理论不断完善、方法不断更新，很多新的算法也在不断涌现[10～12]。
2 反演方法
目前大地电磁的反演方法可分为两类：一类是直接反演方法，即直接从观测数据出发求得地电模型；另一类称为间接反演，即首先给出初始模型，经过正演计算，然后用模型的地表电磁响应对观测数据进行反演拟合，拟合过程中采用各种优化方法和数学手段，寻求最佳拟合的地电模型，间接反演又分为线性反演和非线性反演[13～15]。
2.1 BOSTICK反演
Bostick反演方法是基于大地电磁测深曲线低频渐进线的性质，将视电阻率随周期变化的曲线变换成为电阻率随深度变化的曲线[16]。
BOSTICK反演公式为：
(1)
它可以提供与实测视电阻率曲线相对应的电阻率随深度变化的地电模型，并且可以根据电阻率-深度曲线的拐点确定电性界面的位置。该方法的特点是理论基础简单、容易计算，但
由于其反演结果受数据误差影响大，反演精度不高，很难满足实际应用的需要，但经常被用作现场实时处理和构建间接反演的初始模型。
2.2 OCCAM反演法
由S.C.Constable 等提出的Occam 反演是一种由电磁测深数据产生光滑模型的实用算法，并且是一种带平滑约束的非线性最小二乘解决方案[17，18]。该方法的原理是：在保证电性分布连续或光滑的条件下，寻求与实测数据拟合最好的地电模型，直到达到指定的拟合精度。Occam反演要寻求的解是：能尽可能的与实际观测数据相吻合，同时又具有最小粗糙度的地电结构。
正演问题的解可以表示为：dj=Fj[m],j=1,2L,M。包含M个实测值的数据可表示为d∈EM，模型可记为m∈EN,Fj 是与第j个数据相联系的正演涵数，用向量可表示为：d=F[m]。数据拟合差可写成：X2=‖wd-wF[m]‖2，其中：
w=diag{1/σ1,1/σ2,L,1/σM}
(2)
最优化过程是对由拉格朗日(Lagrange)乘子构造的一个无约束的目标函数U最小化：
U=R1+μ-1{‖wd-wF[m]‖2-X*2}
(3)
式中R1为粗糙度，μ-1为拉格朗日乘子，X*2为反演所要求达到的拟合误差。
在反演迭代过程中目标函数U将趋于极小值。因此可令mU(U在m处的梯度)等于0，可得模型向量m满足：
μ-1(WJ)TWJm-μ-1(WJ)Twd+∂T∂m=0
或：
μ-1(WJ)TWJm-μ-1(WJ)Twd
(4)
式中的M×N的矩阵J是雅可比矩阵：
J=mF
(5)
Occam 反演结果对初始模型依赖程度较小，当然给定的初始模型越好，同样会更容易得到好的反演结果，而且会加速收敛，减少反演过程的迭代次数。Occam 反演的迭代过程稳定，具有收敛速度快的特点，一般只需要几次迭代就能得到较理想的反演结果。在反演的每一次迭代过程中都求解一次完整的关于偏导数矩阵(N×M)，它虽然可能得到较为详细的模型参数，但同时反演过程也会随着测点数的增加以及初始模型网格的加密而变得缓慢，另外随着迭代次数的增加，容易陷入冗余的构造。
2.3 RRI反演法
为了避免像Occam等法直接线性搜索带来的繁重计算，Smith等人提出，可以通过解与一维相近的反演问题，来计算在每个测量位置下面的电阻率扰动，把二维反演问题转化为一系列一维反演问题[19]。RRI关键是用前一次迭代二维正演产生的电、磁场来近似计算场的横向偏导数项，得到一个近似线性泛函式。据此可用类似一维的反演方法求得各测点下的模型值，经插值形成新的二维电性断面重新迭代。
在RRI中目标函数离散形式为：
W=(Rm-c)T(Rm-c)+β eTe
(6)
式中，R就是所谓的测点粗糙度矩阵，m是尝试模型，c是偏差向量，β为拉格朗日常数，上标T表示矩阵转置。数据残差e由下式确定：
(d-e)-d0=Fm-Fm0
(7)
m0和m分别表示初始模型和迭代产生的新模型，d0和d分别为正演计算数据和测量数据，F为Frechet偏导数矩阵。通过给定初始模型m0计算出d0和各测点下的积分核函数及数据残差e，然后对目标函数求极小得到模型改正量及新模型m，再以新模型为初始模型重复上述的过程逐次迭代直至结果满足精度要求。
RRI对初始模型有较大的依赖性，当给定的初始模型条件较好时，反演效果才能得到保证。
反演的迭代过程基本是稳定收敛的，这与在算法中引入了水平最平缓约束有关。RRI 避免了求解雅可比矩阵，减少了内存需求，极大地提高了计算速度。但由于测点间的模型参数完全采用插值求取，因此在遇到大区域、稀疏测点的二维反演时，模型的不确定性较大。
2.4 NLCG反演方法
Rodi和Mackie在2001年提出了非线性共轭梯度法(NLCG)[20]。将反演问题表示成：
d=F(m)+e
(8)
式中，d为数据矢量，m为模型矢量，F为正演响应。把m设成模型网格共M个元素，每个MR代表唯一网格的对数视电阻率。定义ψ：
ψ(m)=(d-F(m))TV-1(d-F(m))
+λmTLTLm
(9)
给定λ、V和L的值。正则化因子λ是一个正数。正定矩阵V表示探测数据的误差。当模型网格唯一时，Lm近似为logρ的拉普拉斯算子。
定义目标函数的梯度和赫赛函数分别为：
g(1×M)和H(H×M)。
令A表示正演函数F的雅可比矩阵：
Aij(m)=∂jFi(m),i=1,…，N，j=1,…，M带入式(9)得到：
g(m)=-2A(m)TV-1(d-F(m))
+2λLTLm
(10)
H(m)=2A(m)TV-1A(m)+2λLTL
(11)
式中Bi和Fi的赫赛函数，q=V-1(d-F(m))。
式(10)中的A(m)TV-1(d-F(m))以及式(11)中的A(m)TV-1A(m)可以不进行雅可比矩阵A的计算而直接求得。
使用非线性共轭梯度直接解决最小化问题，算法引用了Polak-Ribiere的非线性共轭梯度方法来求解式(9)极小值。其求解过程如下：
m0=given
ψ(mk+αkpk)=minψ(mk+αpk)
mk+1=mk+αkpk,k=0,1,2……
(12)
搜索方向迭代如下：
p0=-C0g0
pk=-Ckgk+βkpk-1,k=1,2,……
(13)
(14)
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议