当前位置：文档之家 › 大地电磁非线性共轭梯度拟三维反演

大地电磁非线性共轭梯度拟三维反演

  • 格式：pdf
  • 大小：966.41 KB
  • 文档页数：9

胡祖志,胡祥云,何展翔.大地电磁非线性共轭梯度拟三维反演.地球物理学报,2006,49(4):1226～1234Hu Z Z ,Hu X Y,He Z X.Pseudo 2Three 2Dimensional magnetotelluric inversion using nonlinear conjugate gradients.Chinese J .G eophys .(in Chinese ),2006,49(4):1226～1234大地电磁非线性共轭梯度拟三维反演胡祖志1,2,胡祥云1,何展翔21中国地质大学地球物理与空间信息学院,武汉　4300742东方地球物理公司综合物化探事业部,河北涿州　072751摘　要　提出了非线性共轭梯度法大地电磁拟三维反演.该方法选取共轭梯度反演算法为拟三维反演的核心.在计算灵敏度(Jacobian )矩阵时,吸取近似灵敏度矩阵思想,采用一维灵敏度矩阵来代替三维灵敏度矩阵,并对非测点的灵敏度元素提出一种近似方法.在第一次反演之后,采用拟牛顿法更新灵敏度矩阵.拟三维反演法在很大程度上节省了计算时间,并且理论模型和实际资料的反演试算结果表明大地电磁拟三维反演法具有一定的实用价值.关键词　大地电磁,拟三维,灵敏度,非线性共轭梯度,拟牛顿法文章编号　0001-5733(2006)04-1226-09中图分类号　P631收稿日期　2005-05-12,2006-03-16收修定稿作者简介　胡祖志,男,1981年生,硕士,主要从事电磁法勘探与正反演研究.E 2mail :airfoxh @Pseudo 2Three 2Dimensional m agnetotelluric inversion usingnonlinear conjugate gradientsHU Zu 2Zhi 1,2,HU X iang 2Y un 1,HE Zhan 2X iang 21Institute o f G eophysics &G eomatics ,China Univer sity o f G eo sciences ,Wuhan 430074,China2Department o f Non 2seismic Exploration ,Bureau o f G eophysical Pro specting ,H ebei Zhuo zhou 072751,ChinaAbstract We propose nonlinear conjugate gradient pseudo 232D MT inversion alg orithm.The conjugate gradient method was selected as the inversion kernel.When calculating sensitivities matrix ,instead of com puting 32D sensitivities ,we attem pt to com pute 12D sensitivities.Since the sensitivity elements of non 2stations are zero ,which should be non 2negligible values ,we adopt a method of approximation.We update sensitivities matrix by the quasi 2Newton method after the first inversion.Psudeo 232D MT inversion can save com putational time greatly.Inversion results of a synthetic data set and real data set show that the psudeo 232D MT inversion has s ome practical applications.K eyw ords Magnetotelluric ,Pseudo 2three 2dimensional ,Sensitivity ,Nonlinear conjugate gradient ,Quasi 2Newton method1　引　言 大地电磁(MT )是前苏联学者T ikhonov [1]和法国学者Cagniard [2]在20世纪50年代提出来的利用天然交变电磁场研究地球电性结构的一种地球物理勘探方法,该方法在国内外的地热田的调查、矿产普查和勘探、地壳和上地幔电性结构的研究[3]、海洋地球物理、环境地球物理研究和地质工程中,尤其在石油天然气勘探方面,发挥着举足轻重的作用[4].当前MT 资料处理使用最多的还是一维和二维正反演.一维大地电磁反演理论比较成熟,方法众多.虽然在原理上二维和三维反演方法与一维方法在本质上是一致的,但是二维和三维反演方法在20世纪80年代前期研究进展相对缓慢.随着计算机技术的飞速发展,高维MT 反演方法也得到了快速发第49卷第4期2006年7月地　球　物　理　学　报CHI NESE JOURNA L OF GE OPHY SICSV ol.49,N o.4Jul.,2006展.目前国际上比较常用的二维反演方法有Occam 反演法[5]、RRI反演法[6]、共轭梯度反演法[7],国内有视模反演法[8]、连续介质快速反演法[9]、拟波动方程反演法[10]、模拟退火法[11]等.自20世纪70年代开始,三维大地电磁正演模拟算法得到了很大的发展.最初很多是采用积分方程法[12,13].在层状模型中仅有几个异常体时,这种方法计算的精度高,速度比较快,但随着模型复杂度的增加,花费的计算时间也随之剧增.有限差分法[14～16]、有限元法[17,18]、边界元法[19]、有限体积法[20]、整体积分局部微分法[21]、基于传统方法上进行一些改进的混合算法[22]和解大型稀疏矩阵方法的进步,使得它们在三维模拟中得到广泛的应用.随着三维正演的发展,MT的三维反演研究也日趋升温,近几年来的SEG年会专门开辟了三维电磁模拟与反演专题.国内的中国地质大学、中国石油大学、中国地震局地质研究所以及东方地球物理公司等单位也开展了三维模拟的应用研究[23].目前国内外提出的三维反演方法众多[24],有以共轭梯度法极大似然反演[25]、非线性共轭梯度反演[26]和人工神经网络反演[27]为代表的真三维反演方法,还有以拟线性近似反演[28]、快速松弛反演[29]等为代表的近似三维反演方法.虽然Madden et al.[14]提出了使用共轭梯度松弛法进行大地电磁三维反演,但共轭梯度法仍以其良好的稳定性和内存需求不高的特点,赢得了很多学者的关注,如S pitzer[30]研究了利用共轭梯度法进行电阻率法的三维正演问题;Zhang et al.[31]、吴小平等[32]分别研究了利用共轭梯度法进行电阻率法的三维反演问题;R odi et al.[7]在共轭梯度法的基础上又提出非线性共轭梯度(N LCG)法大地电磁二维反演,取得显著效果;而Newman et al.[26]则使用非线性共轭梯度法解决三维大地电磁反演问题.他们使用N LCG三维反演方法在串行机和并行机上反演合成模型数据,也说明了它的有效性,但还未见到应用该方法对实际三维资料进行处理的报道.地球物理反演中灵敏度矩阵的求取是反演中一个重要的分支.许多学者提出了很多计算灵敏度矩阵的方法,如扰动法、灵敏度方程法、伴随方程法和近似法等.在大地电磁近似反演方法中,Smith et al.[6]、Ellis et al.[33]曾经用近似灵敏度矩阵反演大地电磁二维数据,Farquhars on et al.[34]在伴随方程法基础上提出了一种灵敏度矩阵的近似方法,戴世坤等[9]提出连续介质反演则采用等效一维模型的方法计算灵敏度矩阵,谭捍东等[29]在RRI反演法的基础上进行了改进,在求取灵敏度矩阵时也采用近似法,发展了快速大地电磁三维反演方法.引入了近似计算灵敏度矩阵的方法使得反演的时间大大缩短,并且反演的结果比较令人满意.基于此,我们提出非线性共轭梯度法大地电磁拟三维反演法.三维正演采用交错采样有限差分方法[15],在计算三维灵敏度矩阵时,用一维灵敏度矩阵代替,对非测点的灵敏度元素按某一规则进行插值构建,从而形成拟三维灵敏度矩阵,并且在迭代反演过程中,采用拟牛顿法对拟灵敏度矩阵进行更新.2　拟三维反演方法211　目标泛函的构建假定离散化的三维地电模型向量m=[m1, m2,…,m M]T,观测数据向量d=[d1,d2,…,d N]T,误差向量e=[e1,e2,…,e N]T,经典的地球物理反演问题一般表示为d=F(m)+e,(1)其中F是正演模拟算子.在T ikhonov正则化理论基础上构建三维反演的目标泛函Φ:Φ(m)=d-F(m)T C-1e d-F(m)+λm T W T Wm,(2) C-1e是正定矩阵,为误差向量e的协方差矩阵.(2)式的第二项定义了模型空间的稳定化泛函.在拟三维反演中,选取矩阵W为简单的二阶差分算子来近似逼近Laplacian算子Δ2.正则化参数λ为一个正数,用来控制模型拟合的光滑度,在极小化数据拟合差范数和模型范数之间充当“折中”的角色.在选取λ参数时需谨慎,λ越大模型越光滑,当λ→0,反演问题变成接近于病态的最小二乘反演方法,由此会得到一个不稳定模型[35].212　非线性共轭梯度法最速下降法是最简单的梯度方法,但是在实际应用中经常收敛很慢.Fletcher and Reeves[37]提出了一种较好的非线性共轭梯度法(N LCG),之后P olak[38]对此进行了改进.这种方法与Hestenes and Stiefel[39]提出的线性共轭梯度法(CG)很密切,当目标泛函是正定二次并且N LCG中使用了精确的线性搜索,则N LCG就与CG方法等价.非线性共轭梯度法已经广泛地应用于地球物理反演问题[7,26,40]当中.我们在大地电磁拟三维反演中7221　4期胡祖志等:大地电磁非线性共轭梯度拟三维反演采用P olak2Ribiere共轭梯度法去极小化目标泛函(2)式.N LCG法沿搜索方向线性极小化确定一个模型序列:m0=mΦ(ml +αl p l)=minαΦ(ml+αp l)(3)m l+1=m l+αl p l l=0,1,2,…. m为给出的已知初始模型向量序列,在每次迭代极小化目标泛函时需要沿搜索方向pl计算搜索步长αl,搜索的迭代方向为p0=-C0g0p l=-C l g l+βl p l-1 l=1,2,…,(4)其中,gl代表Φ在m=ml处的梯度,C l是预条件因子,βl有多种形式,不同的形式代表不同的共轭梯度方法.本文选取PR(P olak2Ribiere)共轭梯度法βl =g T l C l(g l-g l-1)g T l-1C l g l-1.(5)213　预条件预条件因子Cl对共轭梯度法有很大的影响, C l越近似Hessian矩阵(目标函数对模型参数的二次导数),共轭梯度法执行的效率越高.R odi et al.[7]提出如下的方法构造预条件因子:C l=(κl I+λW T W)-1.(6) Newman[41]则采用另一种近似Hessian矩阵的方法:C l=[diag(A T C-1e A)+λW T W]-1.(7)式中κl是一个特定的标量,A为灵敏度矩阵,diag表示提取矩阵A T C-1eA的对角数值.前一种方法构造简单,后一种方法需要对灵敏度矩阵进行矩阵计算,需要很大的存储空间.在计算方程(4)时,可以通过下式应用线性共轭梯度法进行求解:C-1l q=g.(8) 214　近似灵敏度矩阵思想目标泛函(1)式的梯度为ΔΦ(m)=-2A(m)T C-1e d-F(m)+2λW T Wm,(9)其中A ij(m)=9j F i(m),i=1,2,…,N;j=1,2,…,M(10)代表Jacobian矩阵或者灵敏度矩阵.21411　对非测点的灵敏度矩阵元素的插值在拟三维反演中,可以选取一维灵敏度矩阵代替真三维灵敏度矩阵,除了与测点位置相对应的灵敏度元素不为零外,其他的元素皆为零,A=0Aj1N10Aj2N2……0AjNsNs,(11)其中,jk代表第k个测点在灵敏度矩阵中第j k列的位置,Ns为测点个数.AjNiNi=9F19m1…9F19m k…9F19m Nz9φ19m1…9φ19m k…9φ19m Nz……………9F Nf9m1…9F Nf9m k…9F Nf9m Nz9φNf9m1…9φNf9m k…9φNf9m Nz,(12)为第Ni个测点的一维灵敏度矩阵,N f为频点数,φ为相位.对非测点的灵敏度元素,根据电磁波的传播特性,并类比地震波的吸收[42]规律,我们提出如下方法进行近似:A～jNi=AjNiNie-ks j,j=1,2,…,M,(13)其中,k为与频率有关的系数,sj=Δx2j+Δy2j+Δz2j,Δxj,Δy j和Δz j分别为第j个非测点与测点下面长方体单元的x、y和z方向的距离.21412　拟牛顿更新灵敏度矩阵在每次反演迭代中,灵敏度矩阵的计算要占用很长的时间.为了节省计算时间,在第一次反演之后,我们采用拟牛顿法更新灵敏度矩阵[43]:A n+1=A n+[Δd n-A nΔm n]Δm TnΔm TnΔmn,(14)其中,Δdn=d n-d n-1,Δm n=m n-m n-1.在试验中发现,采用拟牛顿法更新灵敏度矩阵所占用的计算机CPU时间比直接计算灵敏度矩阵所用的时间减少了1Π2～1Π3.215　线性搜索在极小化目标泛函Φ(m+αp)时,需要确定搜索步长α,因为只包含一个变量α,所以可以用全局最优化的一维线性搜索方法,如精确线性搜索法、黄金分割法、插值法和不精确线性搜索法等求取.但是全局最优化方法有可能在每次反演迭代时导致多次8221地球物理学报(Chinese J.G eophys.)49卷　三维正演计算,会极大地降低反演的效率.R odi and Mackie [7]提出了一种改进的单变量高斯-牛顿线性搜索方法,Ψl 表示极小化的单变量泛函,Ψ～l 表示它的高斯-牛顿近似:Ψl =Φ(m l +αp l ),(15)Ψ～l (α;m lin =Φ～(m l +αp l ;m lin ),(16)其中m lin 表示一线性化的模型.线性搜索产生一系列模型:m l ,k =m l +αl ,k p l ,k =0,1,2,…,n k ,(17)其中αl ,0=0,Ψ～l (αl ,k +1;m l ,k )=min αΦ～(α;m l ,k ),k =0,1,2,…,n k .(18)因为Φ～(m ;m l ,k )是关于m 的二次函数,Ψ～l (α;m l ,k )是关于α的二次函数,则很容易求解(18)式极小化方程αl ,k +1=αl ,k-g Tl ,k p lp T l H ～l ,k p l,(19)其中,H ～是Hessian 矩阵H 的近似矩阵.在线性搜索中,我们采用如下的方案:(1)保持在线性搜索中得到的最优模型(最小的Φ),记做m l ,best =m l +αl ,best p l ;(2)在迭代中如果Ψl (αl ,k )>Ψl (αl ,k -1),则通过二等分方法计算下步的步长αl ,k +1=(αl ,k +αl ,best )2;(20) (3)在线性搜索的第2步或者其后续的步骤里,如果当前的模型和先前的最优模型有一个极小最优模型,也就是说当Ψ′l (αl ,k )Ψ′l (αl ,best )<0时,则计算αl ,k +1,以获得用三次近似Ψl (α)的局部极小.在α=αl ,k 和α=αl ,best 两点处,三次近似值与Ψl 及三次近似的导数值与Ψ′l 值相等.当计算的目标泛函和插值近似的目标泛函满足下面的不等式时,就可以认为线性搜索是成功的:Ψl (αl ,k +1)-Ψ～l (αl ,k +1;m l ,k )≤τΨl (αl ,k +1),(21)式中τν1;如果在允许最大线性搜索次数内Ψl (αl ,k +1)>115Ψl (αl ,best ),则线性搜索失败.最终经l 次线性搜索后的结果被确定为最优模型m l +1=m l ,best .(22) 如果线性搜索成功,新的搜索方向p l +1将由(4)式计算;如果线性搜索失败,p l +1将由(4)式的右边第一项计算得到.216　收敛标准第k 次迭代单一模式视电阻率和相位联合反演数据拟合差可由下式计算:χ2XY =6Nsj =16Nfi =1ln ρXY obsa ij -ln ρXY cala ij εXY obsij 2+φXY obs ij ΠφXY calij εφXY obs ij 2,(23)χ2YX =6Nsj =16N fi =1ln ρYX obsa ij-ln ρYX cala ijεYX obs ij2+φYX obs ij ΠφYX cal ij εφYX obsij 2,(24)式中ρXY obs a ij 、φXY obs a ij 、ρYX obs a ij 和φYX obsij 分别是XY 模式和YX 模式的视电阻率与相位的观测值,ρXY cal a ij 、φXY cal a ij 、ρYX cala ij 和φYX calij 分别是XY 模式和YX 模式的视电阻率与相位的计算值,εXY obsij、εφXY obs ij 、εXY obs ij和εφXY obsij 分别是XY 模式和YX 模式的视电阻率与相位的方差.对于三维大地电磁的单一模式反演,判断其收敛标准为χ—2XY =χ2XY 2・N s ・N f≤1,(25)χ—2YX=χ2YX2・N s ・N f≤1.(26) 对于联合反演,判断反演的收敛标准为χ—2both =χ—2XY +χ—2YX 2≤1.(27)3　理论模型与实例 为了对前面建立的拟三维反演算法的有效性进行检验,我们对理论合成模型数据和实测资料进行了反演计算.下面的程序均在主频为214GH z ,内存为512M 的奔腾Ⅳ处理器的PC 机上运行.311　模型一:二维长方体二维长方体模型如图1a 所示.二维长方体的走向为x 方向,尺寸为1km ×2km ,电阻率为5Ωm ,顶面埋深为700m ,背景电阻率为100Ωm.有限差分三维正演计算网格为N x =54,N y =54,N z =43(N z 包含9221　4期胡祖志等:大地电磁非线性共轭梯度拟三维反演图1　二维长方体及三维长方体模型(a )二维长方体Π三维长方体在x =0处的垂直断面图;(b )三维长方体在z =700～1700m 间的平面俯视图.Fig.1　Geometry of 22D cuboid m odel and 32D cuboid m odel(a )Cross section of a 22D Π32D cuboid m odel at x =0;(b )Vertial view of a 32D cuboid m odel at z =700～1700m.图2　模型一的XY 模式拟三维反演结果的水平切片图和垂直断面(a )在z =1000m 处的水平切片图;(b )在x =0处的垂直断面图.图中黑线代表真实模型边界(下同).Fig.2　Plan view and cross section of the inverted m odel from pseudo 232D XY m ode(a )Plan view of z =1000m ;(b )Cross section of x =0.Black line denotes the boundary of the m odel.图3　模型二的YX 模式拟三维反演结果的水平切片图和垂直断面(a )在z =1000m 处的水平切片图;(b )在x =0处的垂直断面图.Fig.3　Plan view and cross section of inverted m odel from pseudo 23D YX m ode(a )Plan view of z =1000m ;(b )Cross section of x =0.0321地球物理学报(Chinese J.G eophys.)49卷　图4　某地区大地电磁XY模式反演模型在不同深度的水平切片“3”代表测点位置.Fig.4　Plan view at different depths of the inverted m odel from pseudo232D XY m ode“3”denotes MT site.空气层).对所示模型进行正演计算,视电阻率和相位值分别加入5%和215%的随机误差,用来模拟实测数据.取初始模型为100Ωm,XY模式拟三维反演经两次迭代反演后退出,耗时54min,拟合差为0192%.图2是模型一的XY模式拟三维反演第二次迭代反演结果.从图2中可以看出,反演结果不但很好地反映了二维体的形态和走向,还得到了与真实模型相近的电阻率值,在真实模型边界范围之内,反1321　4期胡祖志等:大地电磁非线性共轭梯度拟三维反演图5　测线x=2000m(A)和x=4000m(B)的XY模式实测资料与反演模型响应拟断面对比(a),(c)为实测视电阻率和相位拟断面图;(b),(d)为反演的视电阻率和相位拟断面图.Fig.5　C omparis on of apparent resistivity and phase(XY m ode)pseudosections between the observed and predicted data along profile x=2000m(A)and x=4000m(B)(a)and(c)are observed apparent resistivity pseudosection and phase pseudosection,respectively,(b)and(d)are predicted apparent resistivity pseudosection and phase pseudosection,respectively.演的电阻率值最大不超过16Ωm.312　模型二:三维长方体三维长方体模型如图1(a,b)所示,为一个低阻三维长方体模型.长方体电阻率为5Ωm,尺寸为2km ×2km×1km,埋深h=700m,中心位于(0,0, 1200m),背景电阻率为100Ωm.有限差分三维正演计算网格同模型一.选取初始模型为100Ωm,YX模式拟三维反演经3次迭代反演后退出,耗时106min,拟合差为2141%.图3为模型二的YX模式拟三维反演第3次迭代反演结果,它能较好地反映异常体的空间展布和电性特征,在真实模型边界范围之内,反演的电阻率值要比实际值偏大,但最大值不超过40Ωm.从图2和图3中可以看出,XY模式和YX模式的拟三维反演得到的地电模型能够较好地反映真实的理论模型.313　实测资料的反演我们用拟三维反演法对某地区的实测资料进行了反演.选取该地区测线8条,线距1000m,每条测线各有10个测点,点距400m,采用10个频点,范围从0175～320H z.选取该地区的XY资料进行反演.初始模型为10Ωm,经过4次迭代反演,耗时122min,拟合差为1158%.图4是大地电磁XY模式拟三维反演得到的模型在不同深度的水平切片图.图5(A,B)是其中两条测线的XY模式实测资料与反演模型响应拟断面对比图.可以看出,实测与反演的视电阻率拟断面图吻合得比较好,反演的相位拟断面图拟合程度相对要差,但大体的趋势一致.2321地球物理学报(Chinese J.G eophys.)49卷　因此可以说得到的地电模型是能够代表实测资料所反映的真实地电构造所具有的特征.4　结　论 通过对两个理论模型合成数据进行反演试算,结果表明,YX模式、XY模式的拟三维反演都能够很好地重建地电模型,这也证明了我们的拟三维反演算法的可行性与正确性.对实际的大地电磁XY 模式资料进行了反演,得到的地电模型理论响应与实测资料拟合得比较好.反演试算的结果表明大地电磁拟三维反演法具有一定的实用价值.致　谢　本项研究中的三维正演程序得到Randy Mackie教授的授权许可,成文过程中得到了王家映教授的指导,在此表示感谢.参考文献(References)[1]　T ikhonov A N.On determining electrical characteristics of the deeplayers of the earth’s crust.Dokl.Akad.Nauk S.S.S.S.R.,1950,73(2):295～297[2]　Cagniard L.Basic theory of the magnetotelluric method of geophysicsprospecting.G eophysics,1953,18(3):605～635[3]　WEI W enbo,Unsw orth,M Jones A,et al.Detection of widespreadfluids in the T ibetan crust by magnetotelluric studies.Science,2001,292:716～718[4]　王家映.石油电法勘探.北京:石油工业出版社,1992 W ang J Y.Oil E lectrical Exploration(in Chinese).Beijing: Petroleum Industry Press,1992[5]　deG root2Hedlin C,C onstable S.Occam’s inversion to generatesm ooth,tw o2dimensional m odels from magnetotelluric data.G eophysics,1990,55(12):1613～1624[6]　Smith J T,Booker J R.Rapid inversion of tw o2and three2dimensionalmagnetotelluric data.J.G eophys.Res.,1991,96:3905～3922 [7]　R odi W L,M ackie R L.N onlinear conjugate gradients alg orithm for22D magnetotelluric inversion.G eophysics,2001,66(1):174～187 [8]　戴世坤.用视模型空间对比进行地球物理反演.地球物理学报,1994,37(增刊Ⅱ):524～533 Dai S K.G eophysical inversion using apparent m odel space contrast method.Chinese J.G eophys.(in Chinese),1994,37(Suppl.Ⅱ):524～533[9]　戴世坤,徐世浙.MT二维和三维连续介质快速反演.石油地球物理勘探,1997,32(3):305～317 Dai S K,Xu S Z.Rapid inversion of magnetotelluric data for22D and32D continuous media.Oil G eophysical Prospecting(inChinese),1997,32(3):305～317[10]　师学明.二维大地电磁拟波动方程反演法研究[博士论文].武汉:中国地质大学,1999 Shi X M.A pseudo wave equation inversion for tw o2dimensional magnetotelluric data[Ph.D.thesis](in Chinese).Wuhan:ChinaUniversity of G eosciences,1999[11]　杨　辉,王永涛,戴世坤等.带地形的MT多参量二维快速模拟退火约束反演.石油地球物理勘探,2003,38(2):213～217 Y ang H,W ang Y T,Dai S K,et al.22D inversion of magnetotelluric multi2parameters with topography using fast simulatedannealing.Oil G eophysical Prospecting(in Chinese),2003,38(2):213～217[12]　H ohmann G W.Three2dimensional induced polarization andelectromagnetic m odeling.G eophysics,1975,40(2):309～324 [13]　W annamaker P E,H ohmann G W,W ard S H.M agnetotelluricresponses of three2dimensional bodies in layered earths.G eophysics,1984,49(9):1517～1533[14]　M adden T R,M ackie R L.Three2dimensional magnetotelluricm odeling and inversion.Proc.IEEE.,1989,77:318～333 [15]　M ackie R L,Smith J T,M adden T R.Three2dimensionalelectromagnetic m odeling using difference equations:Themagnetotelluric exam ple.Radio Science,1994,29(4):923～935 [16]　谭捍东,余钦范,John Booker等.大地电磁三维交错采样有限差分数值模拟.地球物理学报,2003,46(5):705～711 T an H D,Y u Q F,Booker J,et al.M agnetotelluric three2 dimensional m odeling using the staggered2grid finite differencemethod.Chinese J.G eophys.(in Chinese),2003,46(5):705～711[17]　Pridm ore D F,H ohmann G W,W ard S H,et al.An investigation offinite element m odeling for electrical and electromagnetic data in threedimensions.G eophysics,1981,46(7):1009～1024[18]　M itsuhata Y,Uchida T.32D magnetotelluric m odeling using the T2Ωfinite2element method.G eophysics,2004,69(1):108～119 [19]　徐世浙,阮百尧,周　辉等.大地电磁场三维地形影响的数值模拟.中国科学(D辑),1997,27(1):15～20 Xu S Z,Ruan B Y,Zhou H,et al.Numeric m odeling of the three2 dimensional topographic effects in magnetotelluric.Science in China(Series D)(in Chinese),1997,27(1):15～20[20]　Haber E,Ascher U M.Fast finite v olume simulation3Delectromagnetic problems with highly discontinuous coeffcients.SIAMput,2001,22(6):1943～1961[21]　X ie G,Li J,M ajer E,et al.32D electromagnetic m odeling andnonlinear inversion.G eophysics,2000,65(3):804～822[22]　Lee K H,Pridm ore D F,M orris on H F.A hybrid three2dimensionalelectromagnetic m odeling scheme.G eophysics,1981,46(5):796～805.[23]　何展翔,贺振华,王绪本等.油气非地震勘探技术的发展趋势.地球物理学进展,2002,17(3):473～479 He Z X,He Z H,W ang X B,et al.T endency of non2seismic techniques of hydrocarbon prospecting.Progress in G eophysics(inChinese),2002,17(3):473～479[24]　胡祖志,胡祥云.三维大地电磁反演综述.地球物理学进展,2005,20(1):214～220 Hu Z Z,Hu X Y.Review of three dimensional magnetotelluric inversion methods.Progress in G eophysics(in Chinese),2005,20(1):214～220[25]　M ackie R L,M adden T R.Three2dimensional magnetotelluric3321　4期胡祖志等:大地电磁非线性共轭梯度拟三维反演inversion using conjugate gradients.G eophys.J.Int.,1993,115:215～229[26]　Newman G A,Alumbaugh D L.Three2dimensional magnetotelluricinversion using non2linear conjugate gradients.G eophys.J.Int.,2000,140:410～424[27]　S pichak V,P opova.Artificial neural netw ork inversion ofmagnetotelluric data in terms of three2deimensional earthmacroparameters.G eophys.J.Int.,2000,142:15～26 [28]　Zhdanov M S,Fang S,Hursan G.E lectromagnetic inversion usingquasi2linear approximation.G eophysics,2000,65(5):1501～1513 [29]　谭捍东,余钦范,Booker J等.大地电磁三维快速松弛反演.地球物理学报,2003,46(6):850～855 T an H D,Y u Q F,Booker J,et al.Three2dimensional rapid relaxation inversion for the magnetotelluric method.Chinese J.G eophys.(in Chinese),2003,46(6):850～855[30]　S pitzer K.A32D finite difference alg orithm for DC resistivitym odeling using conjugate gradient methods.G eophys.J.Int.,1995,123:903～914[31]　Zhang J,M ackie R L,M adden T R.32D resistivity forward m odelingand inversion using conjugate gradients.G eophysics,1995,60(5):1313～1325[32]　吴小平,徐果明.利用共轭梯度方法电阻率三维反演研究.地球物理学报,2000,43(3):420～427 Wu X P,Xu GM.S tudy on32D resistivity inversion using conjugate gradient method.Chinese J.G eophys.(in Chinese),2000,43(3):420～427[33]　E llis R G,Farquhars on C G,Oldenburg D W.Approximate inversemapping inversion of COPROD22Data.J.G eomag.G eoelectr,1993,45:1001～1012[34]　Farquhars on C G,Oldenburg D W.Approximate sensitivities for theelectromagnetic inverse problem.G eophys.J.Int.,1996,126:235～252[35]　R odi W L,M ackie R L.N onlinear conjugate gradients alg orithm for22D magnetotelluric inversion.G eophysics,2001,66(1):174～187 [36]　Parker R L.The inverse problem of electromagnetic induction:existence and construction of s olutions based upon incom plete data.J.G eophys.Res.,1980,85:4421～4425[37]　Fletcher R,Reeves C M.Function minimization by conjugatepute J.,1964,7:149～154[38]　P olak E.C om putational M ethods in Optimization:A UnifiedApproach.Academic Press,1971[39]　Hestenes M R,S tiefel E.M ethods of conjugate gradients for s olvinglinear systems.Journal o f Research o f the National Bureau o fStandards,1952,49:409～436[40]　E llis R G,Oldenburg D W.The pole2pole32D Dc2resistivity inverseproblem:A conjugate gradient approach.G eophys.J.Int.,1994,119:187～194[41]　Newman GA.S olution accelerators for large2scale3D electromagneticinverse problems.74th Ann.Internat.M tg.Expl.G eophys.Expanded Abstracts,2004[42]　张胜业,潘玉玲.应用地球物理学原理.武汉:中国地质大学出版社,2000 Zhang S Y,Pan Y L.Fundamentals of Applied G eophysics(in Chinese).Wuhan:China University of G eosciences Press,2000 [43]　Christiansen A V,Auken E.Optimizing a layered and laterallyconstrained2D inversion of resistivity data using Broyden’s updateand1D derivatives.J.Appl.G eophys.,2004,56:247～261(本文编辑　汪海英)4321地球物理学报(Chinese J.G eophys.)49卷　。

合集下载

大地电磁三维反演方法综述

大地电磁三维反演方法综述

万方数据　万方数据　万方数据　万方数据　万方数据　万方数据　万方数据大地电磁三维反演方法综述作者：胡祖志， 胡祥云， HU Zu-zhi， HU Xiang-yun作者单位：胡祖志,HU Zu-zhi(中国地质大学地球物理与空间信息学院,武汉,430074)， 胡祥云,HU Xiang-yun(中国地质大学地球物理与空间信息学院,武汉,430074;中国科学院地质与地球物理研究所,北京,100029)刊名：地球物理学进展英文刊名：PROGRESS IN GEOPHYSICS年，卷(期)：2005,20(1)被引用次数：17次1.Raiche A A pattern recognition approach to geophysical inversion using neural nets[外文期刊] 19912.Smith J T;BookerJR Rapid inversion of two- and three- dimensional magnetotelluric data 19913.deGroot-Hedlin C;ConstableS Occam's inversion to generate smooth,two-dimensional models from magnetotelluric data[外文期刊] 1990(12)4.Rodi W L;MackieRL Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[外文期刊] 2001(01)5.Wannamaker P E;StodtJA;Rijo L Two-dimensional topographic responses in magnetotellurics modelled using finite elements 1986(11)6.陈乐寿有限元法在大地电磁正演计算中的应用及改进 1981(03)7.Rodi W L A technique for improve the accuracy of finite element solution for magnetotelluric data 19768.Smith J T Conservative modeling of 3-D electromagnetic fields,Part Ⅱ:Biconjugate gradientsolution and an accelerator[外文期刊] 1996(05)9.王家映地球物理反演理论 199810.Nagihara S;StuartAH Three-dimensional gravity inversion using simulated annealing:Constraints on the diapiric roots of allochthonous salt structures[外文期刊] 2001(05)11.Mackie R L;SmithJT;Madden T R There-dimensional electromagnetic modeling using difference equations:The magnetotelluric example 1994(04)12.刘延年;冯纯伯神经网络在控制领域中的应用 1994(01)13.Mitsuhata Y;UchidaT3D magnetotelluric modeling using the T-Ω finite-element method[外文期刊] 2004(01)14.徐世浙;阮百饶;周辉大地电磁场三维地形影响的数值模拟 1997(01)15.徐世浙;王庆乙;王军用边界单元法模拟二维地形对大地电磁场的影响[期刊论文]-地球物理学报 1992(03)16.Xiong Z Electromagnetic modeling of three-dimensional structures by the method of systemiteration using integral equations 1992(12)17.Wannamaker P E Advances in three-dimensional magnetotelluric modeling using integral equations 1991(11)18.Wannamaker P E;HohmannGW;SanFilipo W A Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations 1984(01)19.陈超;刘江平;余丰求解位场反演问题的混合编码遗传算法[期刊论文]-地球物理学报 2004(01)20.Smith J T Conservative modeling of 3 D electromagnetic fields,Part Ⅰ:Properties and error analysis[外文期刊] 1996(05)21.Tarantola A Inverse Problem Theory:Methods of Data Fitt ting and Model Parameter EstimationElsevier Publishing Co 198722.杨文采非线性地球物理方法:回顾与展望[期刊论文]-地球物理学进展 2002(02)23.杨辉;王永涛;戴世坤带地形的MT多参量二维快速模拟退火约束反演[期刊论文]-石油地球物理勘探 2003(02)24.师学明;王家映;张胜业多尺度逐次逼近遗传算法反演大地电磁资料[期刊论文]-地球物理学报 2000(01)25.徐义贤;王家映大地电磁的多尺度反演[期刊论文]-地球物理学报 1998(05)26.吴小平;汪彤彤电阻率三维反演方法研究进展[期刊论文]-地球物理学进展 2002(01)27.徐世浙;李予国;刘斌大地电磁Hx型波的二维地形改正的方法与效果[期刊论文]-地球物理学报 1997(06)28.Chouteau M;BouchardK Two-dimensional terrain corrections in magnetotellurics surveys[外文期刊] 1988(06)29.潘和平;刘国强应用BP神经网络预测煤质参数及含气量 1997(02)30.徐迪;马大军;李元熹基于神经元网络的股票市场预测 1997(06)31.魏文博我国大地电磁测深新进展及瞻望[期刊论文]-地球物理学进展 2002(02)32.吴小平;徐果明利用共轭梯度方法电阻率三维反演研究[期刊论文]-地球物理学报 2000(03)33.Zhang J;MackieRL;Madden T R3-D resistivity forward modeling and inversion using conjugate gradients[外文期刊] 1995(05)34.Spitzer K A 3-D finite difference algorithm for DC resistivity modeling using conjugate gradient methods 199535.Madden T R;MackieRL There-dimensional magnetotelluric modeling and inversion[外文期刊] 1989(2)36.Spichak V;Popova Artificial neural network inversion of magnetotelluric data in terms of three deimensional earth macroparameters 200037.Hohmann G W There-dimensional induced polarization and electromagnetic modeling 1975(02)38.Spichak V;MenvilleM;Roussignol M Three-dimensional inversion of the magnetotelluric fields using Bayesian statistics 199539.谭捍东;余钦范;Booker J大地电磁三维快速松弛反演[期刊论文]-地球物理学报 2003(06)40.Zhdanov M S;FangS;Hursan G Electromagnetic inversion using quasi-linear approximation[外文期刊] 2000(05)41.Newman G A;AlumbaughDL Three-dimensional magnetotelluric inversion using non-linear conjugate gradients[外文期刊] 200042.Mackie R L;MaddenTR Three-dimensional magnetotelluric inversion using conjugate gradients[外文期刊] 19931.叶益信.胡祥云.金钢燮.陈玉萍.何梅兴.周晓晨.YE Yi-xin.HU Xiang-yun.KimKangsop.CHEN Yu-ping.HE Mei-xing.ZHOU Xiao-chen大地电磁二维陡边界反演应用效果分析[期刊论文]-地球物理学进展2009,24(2)2.胡祖志.胡祥云.何展翔.HU Zu-Zhi.HU Xiang-Yun.HE Zhan-Xiang大地电磁非线性共轭梯度拟三维反演[期刊论文]-地球物理学报2006,49(4)3.师学明.王家映.张胜业.胡祥云.SHI XUE-MING.WANG JIA-YING.ZHANG SHENG-YE.HU XIANG-YUN多尺度逐次逼近遗传算法反演大地电磁资料[期刊论文]-地球物理学报2000,43(1)4.王若.王妙月.底青云二维大地电磁数据的整体反演[期刊论文]-地球物理学进展2001,16(4)5.吴成平.王卫平.胡祥云.方迎尧.WU Cheng-ping.WANG Wei-ping.HU Xiang-yun.FANG Ying-yao频率域直升机航空电磁法视电阻率转换及应用[期刊论文]-物探与化探2009,33(4)6.胡祖志.胡祥云.吴文鹂.桑卓.HU Zu-zhi.HU Xiang-yun.WU Wen-li.SANG Zhuo大地电磁二维反演方法对比研究[期刊论文]-煤田地质与勘探2005,33(1)7.徐凯军.李桐林.刘展.XU Kai-jun.LI Tong-lin.LIU Zhan基于积分方程法的大地电磁三维反演[期刊论文]-物探化探计算技术2009,31(6)8.胡祖志.何展翔.王财富.陈英.罗卫峰大地电磁拟三维反演及其应用[会议论文]-20059.杨迪琨.胡祥云.YANG Di-Kun.HU Xiang-Yun含噪声数据反演的概率描述[期刊论文]-地球物理学报2008,51(3)10.欧东新.王家林.Ou Dongxin.Wang Jialin二维块状结构大地电磁快速反演[期刊论文]-石油物探2005,44(5)1.仇根根.钟清.刘君平.白大为.袁永真大地电磁测深视电阻率和相位曲线之间近似互算方法及程序实现[期刊论文]-物探化探计算技术 2012(4)2.张春贺.乔德武.李世臻.张颖.杨辉.胡来东.尚应军.徐雷良.柴继堂.谭捍东.刘劲松复杂地区油气地球物理勘探技术集成[期刊论文]-地球物理学报 2011(2)3.卓贤军.陆建勋.赵国泽.底青云极低频探地(WEM)工程[期刊论文]-中国工程科学 2011(9)4.陈向斌.吕庆田.张昆大地电磁测深反演方法现状与评述[期刊论文]-地球物理学进展 2011(5)5.卓贤军.陆建勋"极低频探地工程"在资源探测和地震预测中的应用与展望[期刊论文]-舰船科学技术 2010(6)6.白俊雨.王绪本.周军大地电磁实码广义遗传反演算法研究[期刊论文]-物探化探计算技术 2009(4)7.叶益信.胡祥云.金钢燮.陈玉萍.何梅兴.周晓晨大地电磁二维陡边界反演应用效果分析[期刊论文]-地球物理学进展 2009(2)8.徐凯军.李桐林.刘展基于积分方程法的大地电磁三维反演[期刊论文]-物探化探计算技术 2009(6)9.张罗磊.于鹏.王家林.吴健生光滑模型与尖锐边界结合的MT二维反演方法[期刊论文]-地球物理学报 2009(6)10.王猛.邓明.张启升.陈凯.崔金岭控制海底电磁激发脉冲发射的时间同步技术[期刊论文]-地球物理学进展2009(4)11.林昌洪.潭捍东.佟拓大地电磁法三维共轭梯度反演研究[期刊论文]-应用地球物理（英文版） 2008(4)12.朱传庆.杨书江.李同彬.刘海樱最优化电性分层技术在柯克亚山前带的应用[期刊论文]-地球物理学进展2007(2)13.余年.王绪本.阚瑷珂.朱迎堂倾子和视倾子的研究及在断裂解释中的应用[期刊论文]-工程地球物理学报2007(4)14.屈有恒.张贵宾.赵连锋.文战久井地有限线源三维电阻率反演研究[期刊论文]-地球物理学进展 2007(5)15.刘国栋矿产资源调查的物探方法和仪器设备[期刊论文]-物探与化探 2007(z1)16.曹中林.何展翔.昌彦君MT激电效应的模拟研究及在油气检测中的应用[期刊论文]-地球物理学进展 2006(4)17.胡祖志.胡祥云.何展翔大地电磁非线性共轭梯度拟三维反演[期刊论文]-地球物理学报 2006(4)本文链接：/Periodical_dqwlxjz200501037.aspx。

大地电磁全信息资料三维共轭梯度反演研究_英文_

大地电磁全信息资料三维共轭梯度反演研究_英文_

1Manuscript received by the Editor September 31, 2010; revised manuscript received January 11, 2011.*This work is jointly supported by the National Hi-tech Research and Development Program of China (863 Program) (No. 2007AA09Z310), National Natural Science Foundation of China (Grant No.40774029, 40374024), the Fundamental Research Funds for the Central Universities (Grant No. 2010ZY53), and the Program for New Century Excellent Talents in University (NCET).1. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Beijing, 100083, China.2. Key Laboratory of Geo-detection (China University of Geosciences), Ministry of Education, Beijing, 100083, China.3. School of Geophysics and Information Technology, China University of Geosciences, Beijing, 100083, China.Three-dimensional conjugate gradient inversion ofmagnetotelluric full information data*APPLIED GEOPHYSICS, Vol.8, No.1 (March 2011), P. 1 - 10, 3 Figures.DOI: 10.1007/s11770-010-0266-9Lin Chang-Hong 1,2,3, Tan Han-Dong 1,2,3, and Tong Tuo 1,2,3Abstract : Based on the analysis of impedance tensor data, tipper data, and the conjugate gradient algorithm, we develop a three-dimensional (3D) conjugate gradient algorithm for inverting magnetotelluric full information data determined from five electric and magnetic field components and discuss the method to use the full information data for quantitative interpretation of 3D inversion results. Results from the 3D inversion of synthetic data indicate that the results from inverting full information data which combine the impedance tensor and tipper data are better than results from inverting only the impedance tensor data (or tipper data) in improving resolution and reliability. The synthetic examples also demonstrate the validity and stability of this 3D inversion algorithm.Keywords : Magnetotelluric, full information data, 3D inversion, conjugate gradientIntroductionThe magnetotelluric (MT) method can provide more than 20 parameters for geological interpretation after the observation data are processed, such as apparent resistivity, phase, impedance, tipper, skew, and etc. Apparent resistivity, phase, and impedance are often applied for quantitative interpretation, while others are used for qualitative interpretation. Previously, the most frequently used parameters are the apparent resistivity and phase, especially the 2D apparent resistivity and phase. Some 2D inversion algorithms have been developed to invert 2D apparent resistivity and phase data for the E- and H-polarization modes (Jupp and Vozoff, 1977; DeGroot-Hedlin and Constable, 1990; Smithand Booker, 1991; Siripunvaraporn and Egbert, 2000; Rodi and Mackie, 2001). In recent years, 3D inversion has been a hot topic in electromagnetics. The apparent resistivity and phase data can be used for 3D inversion (Mackie and Madden, 1993; Spichak and Popova, 2000; Zhdanov et al., 2000; Tan et al., 2003b; Hu et al., 2006; Lin et al, 2009) and scholars also developed algorithms which use the impedance data for 3D inversion (Newman and Alumbaugh, 2000; Siripunvaraporn et al., 2005; Avdeev and Avdeeva, 2006; Lin et al., 2008).The magnetotelluric method usually measures five electric and magnetic field components: two horizontal electric components, two horizontal magnetic components, and one vertical magnetic component. Both the impedance tensor data and the apparent resistivity-phase data which can be calculated from impedance2Three-dimensional conjugate gradient inversionare determined from the relationship between the two electric and magnetic fields horizontal components,x x y y E H Z E H ªºªº «»«»¬¼¬¼, 21sZ U ZP , and Arg Z I . Thus, inverting only the impedance data or apparent resistivity-phase data means that only four electric and magneticfield components are used for interpretation. If we want to use the full information data (five electric and magnetic field components), another parameter which is related to the horizontal and vertical magnetic fields, tipper, should be considered in the inversion. Inmagnetotelluric theory, tipper is defined as x z y H H T H ªº«»¬¼.The tipper data were seldom obtained because past magnetotelluric field work usually did not acquire the vertical magnetic fields. It was considered that the tipper data signal-to-noise was low and it was difficult to apply the tipper data for quantitative interpretation. Therefore, the tipper data were only employed for recognition of geological structure strike. However, modern magnetotelluric field practice typically includes measurement of the vertical magnetic fields, particularly at long periods where a tri-axial magnetometer is used. Berdichevsky et al. (2003) studied the tipper data generated from 2D geological models and concluded that the tipper data is affected by magnetotelluric static shift smaller than the impedance tensor data (or the apparent resistivity-phase data). If the tipper data is included in the inversion, the reliability of geoelectrical models can be substantially improved. Some 2D inversion algorithms (Siripunvaraporn and Egbert, 2000; Rodi and Mackie, 2001) allowed inversion of the tipper data, which became included along with the TM and TE apparent resistivity-phase data in 2D profile interpretation. From the 2D inversion of field MT, Tuncer et al. (2006) concluded that the tipper data add useful information to the inversion and enhance its resolution. From the 2D inversion example of the synthetic data, Becken et al. (2008) found that the tipper data are most important for sensing deep anomalies.Based on the analysis of impedance tensor data, tipper data, and the conjugate gradient algorithm, we develop a 3D conjugate gradient inversion algorithm for inverting magnetotelluric full information data determined from five electric and magnetic field components and discuss the importance of full information data for 3D inversion. In this paper, we introduce a 3D conjugate gradient inversion algorithm for magnetotelluric full information data. Forward modeling of the impedance tensor andtipper response will be discussed first. Second, we will discuss the theory of the inversion, including the objective function and the pseudo-forward modeling problem. Finally, the results of 3D inversion with synthetic data will be shown.Forward problemBecause impedance and tipper data are used in 3D inversion, the 3D forward modeling problem needs to consider how to correctly simulate the impedancetensor response xx xy yx yy Z Z Z Z Z ªº«»¬¼ and the tipper response zx zy T T T ªº ¬¼ on a surface generated from complex realistic 3D geology.The magnetotelluric source S can always be divided into two orthogonal sources, SX and SY (Tan et al., 2004). Under the influence of SX, the electric field components on the surface can be denoted by E x 1 and E y 1 and the magnetic field components by H x 1, H y 1, and H z 1. The corresponding values under the influence of SY can be denoted by E x 2, E y 2, H x 2, H y 2, and H z 2. Then, the expressions used to compute the impedance tensor response are:12211221211212211221122121121221,,,.x y x y xx x y x y x x x x xyx y x y y y y y yx x y x y y x y x yy x y x y E H E H Z H H H H E H E H Z H H H H E H E H Z H H H H E H E H Z H H H H ­ ° °° ° °®° ° °°° ¯(1)The expressions used to compute the tipper responseare:2112122112211221,.y z y z zx x y x y x z x Z zy x y x y H H H H T H H H H H H H H T ­ ° °®° ° ¯(2)From equations (1) and (2), we know that the electricand magnetic fields on the surface are needed first if we want to obtain the impedance tensor or tipper responses. For a given polarization mode, the magnetic field H j at3Lin et al.the jth site can be expressed by a vector H (Newman and Alumbaugh, 2000):.h T j j H g H(3)where H j denotes the magnetic field at the center of jth site and H is a vector composed of the three magnetic field components at each grid subsurface cell, which can be transformed into each other using the interpolator vector, h g T j . For three directions x , y , and z , equation (3) can be correspondingly written as:,,.h Tjx jx h Tjy jy h Tjz jz H g H H g H H g H ­ °° ®° °¯ (4)Similarly, the electric field E j for a given polarization mode can also be expressed by a vector H :.e T j j E g H(5)where E j denotes the electric field at the center of jthsite and the vector H denotes the magnetic field onthe model boundaries calculated by the staggered-gridfinite difference method, which can be transformed into each other by using the interpolator vector, e g T j . For components at the two direction x and y of the electric field E j , equation (5) can be written as:,.e Tjx jx e T jy jy E g H E g H ­ °® °¯(6)From equations (4) and (6), we know that the vector H is needed first if we want to obtain the electric and magnetic field on the surface. Vector H can be calculated by the staggered-grid finite difference method (Newman and Alumbaugh, 2000; Tan et al., 2003a). With a staggered grid, a discretization of the integral form of the Maxwell equations: 0H l J S E S,E l H S.d d d d i d V P Z ­ °® °¯³³³³³³³³v v(7) (7)have the forward modeling equation of the magnetic field H at the sampling point of each grid cell as: ,KH s (8)where K is a sparse complex-symmetric matrix and sis the source vector related to the boundary conditions. Two forward modeling equations can be obtained, KH 1 = s 1 and KH 2 = s 2, if the magnetic field values to be solvedunder the action of the source SX and vectors at the right end of column of the equations are denoted by H 1 and s 1 and the magnetic field values to be solved under the action of the source SY and vectors at the right end of the equations are denoted by H 2 and s 2. These two forward modeling equations are solved to give magnetic fields H 1 and H 2 of the two polarization modes at sampling locations of each grid cell.Vector H is obtained by solving equation (8). Substituting the vector H into equations (3) and (5), the electric and magnetic fields on the surface are calculated. Then, the impedance tensor and tipper response of the 3D model is computed by equations (1) and (2).Inversion problemThe objective functionBased on the regularized solution method, the objective function of the 3D conjugate gradient inversionalgorithm for inverting magnetotelluric full informationdata can be defined as:d m \100obs T obs TT \\ 100(())(())()().obs T obs T T D F m V D F m m m L L m m O \\ 100(())(())()(obs T obs T T D F m V D F m m m L L m O (9)where 12,,...,T obs obs obs obsN D D D D ªº ¬¼ is a vector composed of the observed impedance or tipper data with D obsidenoting the real or imaginary part of impedance (or tipper) for a particular observation site and frequency, N is the total number of observations, F (m ) denotes the forward modeling function to calculate the impedance or tipper responses, V is the variance of error vector, λ is the regularization parameter, and L is the simple second-difference Laplacian. We have m = [m 1, m 2, …,m M ] where M is the total number of model blocks andm j is the logarithm of resistivity for a unique block (log ρ), and m 0 is the priori model. The objective functionhas two terms: the first term is data fit ψd and the second term is the model fit ψm being constraints for smoothingthe model. The regularization parameter λ balances the influence of ψd and ψm on the image solution. When boththe impedance data and the tipper data are used in theinversion, data fit is expressed as: .d z t \\Z\ (10)where ψz is the misfit of the impedance data, ψt is the4Three-dimensional conjugate gradient inversion misfit of the tipper data, and the weight factor ω is used to balance the ratio of ψz and ψt . Generally, the tipper data are small and easily influenced by noise in field work and it is difficult to acquire good-quality tipper data. Therefore, the value of ω can be decided by the data quality. The tipper data weights can increase (or decrease) according to the quality of the tipper data.The gradient of the objective function can be expressed as: 1022().T Tg A V e L L m m O (11)where A is the Jacobian matrix and the error vector ise = D obs-F (m ). Calculation of the Jacobian matrix – the pseudo-forward modeling problemThe most important part in the inversion process is to determine the Jacobian matrix. Once the value of theJacobian matrix is calculated, the model update at each iteration can be easily computed. The 3D conjugate gradient inversion algorithm flowchart for inverting magnetotelluric full information data is similar withthat for impedance data. The main computations in the complete inversion process are calculations of the valueof A (m ref )T V -1e and A (m ref )p . We will describe how to compute the values of A (m ref )T V -1e and A (m ref )p when the tipper data are inverted in this algorithm, which is quitesimilar with the process when the impedance data are inverted (Lin et al., 2008).Differentiating with respect to m k , equation (8) yields: 1/(/).k k H m K K m H w w w w (12)Substituting equation (4) into equation (12) and with ∂h g Tj /∂m = 0, we get:111/(/),/(/),/(/).h Tjx k jxk h T jy k jy k h T jz k jz k H m g K K m H H m g K K m HH m g K K m H ­w w w w °°w w w w ®°w w w w °¯ (13) Differentiating with respect to m k , equation (2) yields: k y z k z y y x y x k zx m H H m H H H H H H m T w w w w w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) >k y z k z y y x y x k zx m H H m H H H H H H m T w w w w w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y m H H H H º w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) >k y z k z y y x y x k zx m H H m H H H H H H m T w w w w w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y m H H m H H w w //)(21121)()/211212z y z y k y z H H H H m H H k x y k y x H H m H H w /212121221/(),x y x y H H H H k x z k z x m H H m H H //)(12211)()/122121z x z x k x z H H H H m H H k x y k y x m H H m H H w w //2121k y z k z y y m H H m H H w w w //)(21121)()/211212z y z y k y z H H H H m H H k x y k y x m H H m H H w w w //212121221/(),x y x y H H H H º ¼k x z k z x y m H H m H H //)(12211)()122121z x z x k x z H H H H H k x y k y x m H H m H H w w w //2121ky z k z y y x y x k zx m H H m H H H H H H m T w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w w k x y k y x k x y m H H m H H m H H w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H º w w ¼ (14) k y z k z y y x y x m H H m H H H H H //)((21121221)()//21121221z y z y k y z k z y H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w ///(2121122211221/)/(),x y k x y x yH H m H H H H º w w ¼k x z k z x y x y x m H H m H H H H H H //)((12211221)()//12212112z x z x k x z k z x H H H H m H H m H w w w w k x y k y x k x y m H H m H H m H H w w w w w w ///(212112k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H m H H w w w w k x y k y x k x y m H m H H m H H w w w w w w u ///(212112 k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) z k z y y x y x k zx H m H H H H H H m T w w /)((/121221()//21221y k y z k z y H m H H m H H w w y k y x k x y H m H H m H H w w u //(21122211221/)/(),x y k x y x y H H m H H H H º w w ¼z k z x y x y x k xy H m H H H H H H m T w w /)((/211221()//12112x k x z k z x H m H H m H H w w w w y k y x k x y H m H H m H H w w w w u //(21122211221/)/().x y k x y x y H H m H H H H º w w ¼ y z k z y y x y x k zx H H m H H H H H H m T w w //)((/21121221()//121221z y k y z k z y H H H m H H m H H w w x y k y x k x y H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>x z k z x y x y x k xy H H m H H H H H H m T w w w w w w //)((/12211221()//212112z x k x z k z x H H H m H H m H H w w w w x y k y x k x y H H m H H m H H w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ ky z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼>k x z k z x y x y x k xy m H H m H H H H H H m T w w w w w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) k y z k z y y x y x k zx m H H m H H H H H H m T w w //)((/21121221)()//21121221z y z y k y z k z y H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w u ///(2121122211221/)/(),x y k x y x y H H m H H H H º w w ¼k x z k z x y x y x k xy m H H m H H H H H H m T w w //)((/12211221)()//12212112z x z x k x z k z x H H H H m H H m H H w w k x y k y x k x y m H H m H H m H H w w w w w w u ///(2121122211221/)/().x y k x y x y H H m H H H H º w w ¼ (14) (14)Substituting equation (13) into equation (14) anddifferentiating each tipper element at the jth observationsite (∂T zxj /∂m k , ∂T zyj /∂m k ) with respect to the kth modelparameter m k yields:1112/(/)(/),T T zxj k jtx k jtx k T m g K K m H g K K m Hw w w w w w 2112/(/)(/),T T zxj k jtx k jtx k T m g K K m H g K K m H w w w w w w1112/(/)(/).TT zyj k jty k jty k T m g K K m H g K K m H w w w w w w 2112/(/)(/).T T zyj k jty k jty k T m g K K m H g K K m H w w w w w w >))((2212211Tjy h z T jz h y y x y x T jtx g H g H H H H H g 211222122()()/(h T h Ty z y z x jy y jx x y x H H H H H g H g H H H º ¼))((2212211T jy h z T jz h y y x y x T jtx g H g H H H H H g 22112221221()()/(),h T h T y z y z x jy y jx x y x y H H H H H g H g H H H H º ¼>))((1112212T jz h y T jy h z y x y x T jtx g H g H H H H H g 211211122()()/(h T h T y z y z y jx x jy x y x H H H H H g H g H H H º ¼>))((1112212Tjz h y T jy h z y x y x T jtx g H g H H H H H g 211211122()()/(h T h T y z y z y jx x jy x y x H H H H H g H g H H H º ¼))((1112212Tjz h y T jy h z y x y x T jtx g H g H H H H H g 22112111221()()/(),h T h T y z y z y jx x jy x y x y H H H H H g H g H H H H º ¼))((2212211T jzh x T jx h z y x y x T jty g H g H H H H H g 122122122()()/(h T h T x z x z x jy y jx x y x H H H H H g H g H H H º ¼>))((2212211Tjz h x T jx h z y x y x T jty g H g H H H H H g 122122122()()/(h T h T x z x z x jy y jx x y x H H H H H g H g H H H º ¼))((2212211T jz h x Tjxhz y x y x T jty g H g H H H H H g 21221221221()()/(),h T h T x z x z x jy y jx x y x y H H H H H g H g H H H H º ¼and>))((1112212T jx h z T jz h x y x y x T jty g H g H H H H H g 122111122()()/(h T h T x z x z y jx x jy x y x H H H H H g H g H H H Hº ¼>))((1112212Tjx h z T jz h x y x y x T jty g H g H H H H H g 122111122()()/(h T h T x z x z y jx x jy x y x H H H H H g H g H H H H º ¼))((1112212T jx h z T jz h x y x y x T jty g H g H H H H H g 21221111221()/().h T h T x z x z y jx x jy x y x y H H H H H g H g H H H H º ¼ , (15)where h g T jx is the interpolator vector between H jx and H , h g T jyis the interpolator vector between H jy and H , and h g T jz is the interpolator vector between H jz and H . H jx , H jy and H jz denote the x , y , and z components of the magnetic field at the jth observation site, respectively.5Lin et al.Since the tipper data are used in the inversion and the elements of the Jacobian matrix are the partial derivatives of the data with respect to the model parameters, ∂T zxj /∂m k and ∂T zyj /∂m k in equation (15) are the elements of the Jacobian matrix A . Then, A (m ref )T V -1e can be expressed as:111()(/)().NTT ref n k n n A m V e D m V e w w ¦,(16)where D n is the tipper data, N is the total number of tipperdata (number of observation sites × frequency × 4), and k = 1, 2…, M denotes model parameters.Substituting equation (15) into equation (16) and since matrix K is symmetric, we get:11111211()/()/().NNT T T ref k n n k n n n n A m V e H K m K g V e H K m K g V e w w w w ¦¦1211211/()/().NNT Tk n n k n n n n H K m K g V e H K m K g V ew w w w ¦¦ , (17)where 1g n and 2g n are described in equation (15). Forexample, 1g nT = 1g Tjty and 2g nT = 2g Tjty , if T zyj is the nth datavalue.If we set11111()Nn n n v K g V e ¦,(18)and12121(),Nn n n v K g V e ¦(19)then1111(),N n n n Kv g V e ¦(20) 2121(),N n n n Kv g V e ¦(21)11122()//.T T T ref k k A m V e H K m v H K m v w w w w (22)If we take v 1 as the vector to be solved and n Nn n e V g )(111 ¦as the vector on the right side of the equation, equation (20) is similar to the forward modeling equation (8). Therefore, we call equation (20) a pseudo-forward modeling equation. The value of v 1 can be obtained by solving the pseudo-forward modeling problem one time. Similarly, the value of v 2 can be acquired by solving equation (21). Finally, A (m ref )T V -1e in equation (22) can be obtained with the values of v 1 and v 2.Similarly, derivation of A (m ref )p can be expressed as:11111()(/)(/MMT T ref nk k nk k k A m p g K K m H p g K K m Hw w w w ¦¦211211()(/)(/).M MT Tref n k k n k k k k A m p g K K m H p g K K m H pww w w ¦¦ (23)If we set 1111(/),Mkkk t K K m H p w w ¦(24) 1221(/),Mkk k t K K m Hp w w ¦(25)Then111(/),M k k k Kt K m H p w w ¦ (26)and221(/),Mk k k Kt K m H p w w ¦(27) 1212().T Tref n n A m p g t g t(28)If we take t 1 as the vector to be solved and ¦ w w Mk kk p H m K 11)/(as the vector on the right side of the equation (26), equation (26) also is similar to the forward modeling equation (8). Therefore, we call equation (26) a pseudo-forward modeling equation. The value of t 1 can be obtained by solving the pseudo-forward modeling problem one time. Similarly, the value of t 2 can be acquired by solving equation (27). Finally, A (m ref )p in equation (28) can be obtained by substituting the values of t 1 and t 2 into equation (28).Synthetic exampleBased on the theory and equations, we have developed a procedure to realize the 3D conjugate gradient inversion algorithm for inverting full information data. To verify the correctness of this algorithm, we designed two test models composed of two prisms.Model 1This model is composed of two conductive rectangular prisms which are embedded in a 100 Ω•m half-space6Three-dimensional conjugate gradient inversion(see Figure 1). The 10 Ω•m prism (ρ1) has a dimensions of 1.25 km × 1.25 km × 0.4 km with the top located 0.2 km below the earth’s surface. The 1 Ω•m prism (ρ2) has a dimensions of 1.75 km × 1.75 km × 1.55 km with the top located 1.05 km below the earth’s surface. The size and position of the two prisms are shown in Table 1. The size of the 3D grid is 36 × 36 × 37 (including 10 air layers). The impedance tensor and the tipper data atfour frequencies (10 Hz, 3.3 Hz, 1 Hz, and 0.33 Hz) at the grid center on the earth’s surface (289 sites) were calculated by staggered-grid finite difference modeling. The impedance tensor and tipper data were inverted using the 3D conjugate gradient inversion algorithm after Gaussian random noise was added to the data. All results were computed on a PC with a 2.2 GHz Intel (R) core (TM) 2 Duo CPU and 2 GB of physical memory.Fig. 1 View of the two-prism model 1.(a) Side view along negative Y direction. (b) Top view from negative to positive Z direction.Table 2 Statistical data for three scenariosData used in the inversionIteration numbersRuntimes (minutes)Data fit Z xx + Z xy + Z yx + Z yy662212 1.047T zx + Z zy4755380.997Z xx + Z xy + Z yx + Z yy + T zx + Z zy8378011.071The 3D inversion was carried out in three scenarios: 1)inverting the impedance tensor data (Z xx , Z xy , Z yx and Z yy ); 2) inverting the tipper data (T zx and T zy ); and 3) jointly inverting both the impedance tensor and the tipper data (Z xx , Z xy , Z yx , Z yy , T zx and T zy ). One percent Gaussian random noise was added to the impedance tensor data and three percent Gaussian random noise was added to the tipper data. Table 2 shows the total number of iterations, the CPU runtime, and the data fit for the three scenarios.The inversion results are shown in Figure 2. The top row shows the test model, the second row shows the results from the inversion of the impedance data, the third row shows the results from the inversion of the tipper data, and the fourth row shows the results from the inversion of both the impedance and the tipper data. The black dashed lines show the prism margins. The first column shows the horizontal slices at 0.35 km depth, the second column shows the horizontal slices at 1.8 km depth, the third column shows the vertical slices at Y = -0.375 km along the X axis, and the fourth column shows the vertical slices at X = -0.375 km along the Y axis. From this figure we see that the results of inverting the impedance tensor data show the shallow prism better. For instance, the size, shape, and resistivity (the lowest resistivity value is the same as the real value, 10 Ω•m) of the shallow prism are well recovered. Only its center is a little bit deeper than the real model, the center of the deep prism is much shallower than the real model, and the lowest resistivity value is only recovered to 23 Ω•m. Relatively, the results of inverting the tipper data recover. the deep prism much better. The size and shape of theTable 1 The size and position of the two prismsPrism Resistivity In the x direction In the y direction In the z direction Depth of top ρ110 Ω•m -1.25 Km, 0 Km -1.25 Km, 0 Km -0.2 Km, -0.6 Km 0.2 Km ρ21 Ω•m-0.75 Km, 1 Km-0.75 Km, 1 Km-1.05 Km, -2.6 Km1.05 Kmρ1ρ21.75 km1.25 km1.75 km1.25 kmXYSurface0.2 km0.4 km1.05 km1.25 km 1.55 kmρ2ρ11.75 kmXZ(a) (b)。

大地电磁测深反演方法现状与评述

大地电磁测深反演方法现状与评述


犚 犲 狏 犻 犲 狑狅 犳犿 犪 狀 犲 狋 狅 狋 犲 犾 犾 狌 狉 犻 犮犱 犪 狋 犪 犻 狀 狏 犲 狉 狊 犻 狅 狀犿 犲 狋 犺 狅 犱 狊 犵
， ＣＨＥ ＮＸ ｉ ａ ｎ ｂ ｉ ｎ Ｑ ｉ ｎ ｔ ｉ ａ ｎ， ＨＡＮＧ Ｋ ｕ ｎ Ｌ Ｚ ｇ ｇ
（ 犐 狀 狊 狋 犻 狋 狌 狋 犲 狅 犻 狀 犲 狉 犪 犾犚 犲 狊 狅 狌 狉 犮 犲 狊， 犆 犺 犻 狀 犲 狊 犲犃 犮 犪 犱 犲 犿 犲 狅 犾 狅 犻 犮 犪 犾犛 犮 犻 犲 狀 犮 犲 狊， 犅 犲 犻 犻 狀 ０ ０ ０ ３ ７， 犆 犺 犻 狀 犪） 犳犕 狔狅 犳犌 犵 犼 犵１
犃 犫 狊 狋 狉 犪 犮 狋 ｍ ｒ ｏ ｖ ｉ ｎ ｈ ｅａ ｃ ｃ ｕ ｒ ａ ｃ ｆｄ ａ ｔ ａｉ ｎ ｔ ｅ ｒ ｒ ｅ ｔ ａ ｔ ｉ ｏ ｎｆ ｏ ｒｍ ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃｓ ｏ ｕ ｎ ｄ ｉ ｎ ＭＴ） ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄｈ ａ ｓ Ｉ ｐ ｇｔ ｙｏ ｐ ｇ ｇ（ ，ｗ ａ ｌ ｗ ａ ｓｂ ｅ ｅ ｎａ ｎｉ ｍ ｏ ｒ ｔ ａ ｎ ｔｓ ｕ ｂ ｅ ｃ ｔｏ ｆｔ ｈ ｅｍ ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃｓ ｏ ｕ ｎ ｄ ｉ ｎ ｉ ｅ ｌ ｄ ． Ｉ ｎｔ ｈ ｉ ｓｐ ａ ｅ ｒ ｉ ｔ ｈｔ ｈ ｅｍ ａ ｉ ｎｌ ｉ ｎ ｅｏ ｆｔ ｈ ｅ ｙ ｐ ｊ ｇ ｇｆ ｐ ， ， ｍ ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄｄ ｅ ｖ ｅ ｌ ｏ ｉ ｎ ｉ ｓ ｔ ｏ ｒ ｔ ｈ ｅｈ ｏ ｍ ｅａ ｎ ｄａ ｂ ｒ ｏ ａ ｄｍ ａ ｏ ｒ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄ ｓ ｕ ｃ ｈａ ｓｂ ｏ ｓ ｔ ｉ ｃ ｋ， ｇ ｐ ｇｈ ｙ ｊ ，ｏ ，ｃ ，ｒ ，ｑ ， ｕ ａ ｓ ｉ ｌ ｉ ｎ ｅ ａ ｒａ ｒ ｏ ｘ ｉ ｍ ａ ｔ ｉ ｏ ｎ ｓ ｉ ｍ ｕ ｌ ａ ｔ ｉ ｎ ｅ ｉ ｓ ｍ ｉ ｃ ｃ ｃ ａ ｍ，ｃ ｕ ｒ ｖ ｅｃ ｏ ｍ ａ ｒ ｉ ｓ ｏ ｎ ｏ ｎ ｕ ａ ｔ ｅｇ ｒ ａ ｄ ｉ ｅ ｎ ｔ ａ ｉ ｄｒ ｅ ｌ ａ ｘ ａ ｔ ｉ ｏ ｎ ｐ ｐ ｇｓ ｐ ｊ ｇ ｐ ，ｍ ，ｎ ｓ ｉ ｍ ｕ ｌ ａ ｔ ｅ ｄａ ｎ ｎ ｅ ａ ｌ ｉ ｎ ｕ ｌ ｔ ｉ ｓ ｃ ａ ｌ ｅｇ ｅ ｎ ｅ ｔ ｉ ｃａ ｌ ｏ ｒ ｉ ｔ ｈ ｍ，ｃ ｏ ｎ ｕ ａ ｔ ｅｇ ｒ ａ ｄ ｉ ｅ ｎ ｔｍ ａ ｘ ｉ ｍ ｕ ｍｌ ｉ ｋ ｅ ｌ ｉ ｈ ｏ ｏ ｄ ｏ ｎ ｌ ｉ ｎ ｅ ａ ｒｃ ｏ ｎ ｕ ａ ｔ ｅ ｇ ｇ ｊ ｇ ｊ ｇ ，ｂ ，ａ ，ｑ ｕ ａ ｎ ｔ ｕ ｍ ｐ ａ ｔ ｈｉ ｎ ｔ ｅ ｒ ａ ｌａ ｌ ｏ ｒ ｉ ｔ ｈ ｍ ａ ｎ ｄｄ ａ ｍ ｉ ｎ ａ ｒ ｔ ｉ ｃ ｌ ｅｓ ｗ ａ ｒ ｍ ｒ ａ ｄ ｉ ｅ ｎ ｔ ａ ｅ ｓ ｉ ａ ｎ ｒ ｔ ｉ ｆ ｉ ｃ ｉ ａ ｌｎ ｅ ｕ ｒ ａ ｌｎ ｅ ｔ ｗ ｏ ｒ ｋ ｇ ｇ ｐ ｇｐ ｇ ｙ ｏ ｔ ｉ ｍ ｉ ｚ ａ ｔ ｉ ｏ ｎｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎ ｍ ｅ ｔ ｈ ｏ ｄａ ｒ ｅｒ ｅ ｖ ｉ ｅ ｗ ｅ ｄ ｍ ａ ｉ ｎ ｌ ａ ｓ ｅ ｄｏ ｎｔ ｈ ｅｗ ｏ ｒ ｋａ ｃ ｈ ｉ ｅ ｖ ｅ ｄｉ ｎｔ ｈ ｅａ ｕ ｔ ｈ ｏ ｒ ｓ ′ｇ ｒ ｏ ｕ ．Ｔ ｈ ｅ ｐ ｙｂ ｐ ， ｍ ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄｉ ｓｃ ｌ ａ ｓ ｓ ｉ ｆ ｉ ｅ ｄｉ ｎ ｔ ｏｑ ｕ ａ ｌ ｉ ｔ ａ ｔ ｉ ｖ ｅａ ｒ ｏ ｘ ｉ ｍ ａ ｔ ｅｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄ ｔ ｈ ｅｌ ｉ ｎ ｅ ａ ｒｏ ｒｎ ｏ ｎ ｌ ｉ ｎ ｅ ａ ｒ ｇ ｐ ｐ ｉ ｔ ｅ ｒ ａ ｔ ｉ ｖ ｅｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎ ｍ ｅ ｔ ｈ ｏ ｄｂ ａ ｓ ｅ ｄｏ ｎｔ ｈ ｅｏ ｂ ｅ ｃ ｔ ｉ ｖ ｅｆ ｕ ｎ ｃ ｔ ｉ ｏ ｎａ ｎ ｄｔ ｈ ｅｇ ｌ ｏ ｂ ａ ｌｏ ｔ ｉ ｍ ｉ ｚ ａ ｔ ｉ ｏ ｎｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎ ｍ ｅ ｔ ｈ ｏ ｄｂ ｊ ｐ ｙ ｓ ｕ ｍｍ ａ ｒ ｉ ｚ ｉ ｎ ｎ ｄｃ ｏ ｍ ａ ｒ ｅ ｄａ ｎ ｄｒ ｅ ｖ ｉ ｅ ｗ ｅ ｄｖ ａ ｒ ｉ ｏ ｕ ｓｍ ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄ ｓｆ ｒ ｏ ｍｔ ｈ ｅｏ ｂ ｅ ｃ ｔ ｉ ｖ ｅｆ ｕ ｎ ｃ ｔ ｉ ｏ ｎ ｇａ ｐ ｇ ｊ ， ，ａ ｔ ｈ ｅｓ ｅ ｎ ｓ ｉ ｔ ｉ ｖ ｉ ｔ ａ ｔ ｒ ｉ ｘｃ ａ ｌ ｃ ｕ ｌ ａ ｔ ｉ ｏ ｎａ ｎ ｄｔ ｈ ｅｃ ｏ ｎ ｖ ｅ ｒ ｅ ｎ ｃ ｅｓ ｅ ｅ ｄ ｌ ｓ ｏｐ ｏ ｉ ｎ ｔ ｅ ｄｏ ｕ ｔｔ ｈ ｅｐ ｒ ｏ ｂ ｌ ｅ ｍ ｓ ．Ｔ ｈ ｅ ｃ ｏ ｎ ｓ ｔ ｒ ｕ ｃ ｔ ｉ ｏ ｎ ｙｍ ｇ ｐ ， ｄ ｅ ｖ ｅ ｌ ｏ ｉ ｎ ｒ ｏ ｓ ｅ ｃ ｔ ｓｏ ｆ ｔ ｈ ｅｍ ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄａ ｒ ｅ ａ ｌ ｓ ｏ ｆ ｏ ｒ ｅ ｃ ａ ｓ ｔ ｅ ｄ ． Ｉ ｎ ｔ ｈ ｅ ｆ ｕ ｔ ｕ ｒ ｅ ｔ ｈ ｅ ｒ ｅ ｓ ｅ ａ ｒ ｃ ｈｗ ｏ ｕ ｌ ｄｂ ｅ ｐ ｇｐ ｐ ｇ ， ｆ ｏ ｃ ｕ ｓ ｅ ｄｏ ｎｈ ｅ ｔ ｅ ｒ ｏ ｅ ｎ ｅ ｏ ｕ ｓｇ ｅ ｏ ｅ ｌ ｅ ｃ ｔ ｒ ｉ ｃｓ ｔ ｒ ｕ ｃ ｔ ｕ ｒ ｅ ｓｅ ｌ ｅ ｃ ｔ ｒ ｏ ｍ ａ ｎ ｅ ｔ ｉ ｃｆ ｉ ｅ ｌ ｄｒ ｅ ｓ ｏ ｎ ｓ ｅｗ ｉ ｔ ｈｌ ａ ｎ ｄ ｆ ｏ ｒ ｍｅ ｆ ｆ ｅ ｃ ｔ ｔ ｈ ｅｐ ｒ ａ ｃ ｔ ｉ ｃ ａ ｌｏ ｆ ｇ ｇ ｐ ， ｎ ｏ ｎ ｌ ｉ ｎ ｅ ａ ｒ ｔ ｈ ｒ ｅ ｅ ｄ ｉ ｍ ｅ ｎ ｓ ｉ ｏ ｎ ａ ｌ ｌ ｏ ｂ ａ ｌ ｏ ｔ ｉ ｍ ｉ ｚ ａ ｔ ｉ ｏ ｎ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎａ ｎ ｄｈ ｂ ｒ ｉ ｄ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄ ｄ ｅ ｔ ａ ｉ ｌ ｅ ｄ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄａ ｎ ｄ ｇ ｐ ｙ ｈ ｏ ｗｔ ｈ ｅｐ ａ ｒ ａ ｍ ｅ ｔ ｅ ｒｓ ｅ ｔ ｔ ｉ ｎ ｓａ ｆ ｆ ｅ ｃ ｔｔ ｈ ｅ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｒ ｅ ｓ ｕ ｌ ｔ ｓａ ｎ ｄｊ ｏ ｉ ｎ ｔ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｏ ｆｏ ｔ ｈ ｅ ｒｇ ｅ ｏ ｈ ｓ ｉ ｃ ａ ｌｄ ａ ｔ ａ， ｏ ｒｃ ｏ ｎ ｓ ｔ ｒ ａ ｉ ｎ ｅ ｄ ｇ ｐ ｙ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎ ． ， ， 犓 犲 狑 狅 狉 犱 狊 ａ ｎ ｅ ｔ ｏ ｔ ｅ ｌ ｌ ｕ ｒ ｉ ｃｓ ｏ ｕ ｎ ｄ ｉ ｎ ｉ ｎ ｖ ｅ ｒ ｓ ｉ ｏ ｎｍ ｅ ｔ ｈ ｏ ｄ ｓ ｒ ｅ ｖ ｉ ｅ ｗ ｍ ｇ ｇ 狔
基于数据空间算法的大地电磁各向异性三维反演

基于数据空间算法的大地电磁各向异性三维反演

基于数据空间算法的大地电磁各向异性三维反演大地电磁（EM）法是一种非侵入式的地球物理探测方法，广泛应用于矿产勘探、环境地质和地下水资源调查等领域。

大地电磁法通过测量地面上的电磁场进行数据采集，然后利用数学模型和反演算法，推断出地下的电导率分布。

大地电磁各向异性反演是一项复杂的任务。

各向异性是指地下岩石或土壤的电导率在不同方向上有所差异。

对于大地电磁法来说，各向异性的存在会导致数据解释和模型反演的困难。

因此，开发一种有效的算法来处理大地电磁各向异性问题具有重要意义。

基于数据空间的反演算法是一种可以用于处理反演问题的有效方法。

该方法将数据空间中的电磁场响应看作是一个变量，而不是通过模型计算得出的固定值。

这样做的目的是将反演问题转化为一个优化问题，通过调整模型参数，使模拟的数据与观测数据之间的误差最小化。

大地电磁各向异性三维反演中，关键是在反演过程中同时估计各向异性参数和电导率分布。

一种常用的方法是使用多个相对应的模型，每个模型都对应一个不同的各向异性参数。

然后，将这些模型与观测数据进行匹配，找到与观测数据最匹配的模型。

具体的反演算法包括以下几个步骤：1.初始化模型参数和各向异性参数的范围；2.生成初始模型，并计算相应的电磁场响应；3.对于每个模型，根据观测数据和模拟数据之间的差异，计算一个误差函数；4.使用优化算法（如梯度下降法或牛顿法），调整模型参数，直到误差函数最小化；5.重复步骤3和4，直到得到最佳模型。

在实际应用中，还可以引入先验信息来增加反演的稳定性和准确性。

例如，可以使用地质或地球物理学的先验知识来限制模型参数的范围，或者通过加权处理来调整观测数据和模拟数据之间的差异。

总之，基于数据空间的算法对于大地电磁各向异性三维反演有着重要的意义。

虽然该方法存在一定的计算复杂性，但可以通过合适的优化算法和先验信息引入来提高反演的稳定性和准确性。

随着计算能力的不断提高，基于数据空间算法将在大地电磁反演中发挥更加重要的作用。