大地电磁阻抗张量旋转方法和曲线圆滑方法的比较

大地电磁阻抗张量旋转方法和曲线圆滑方法的比较赵维俊;孙中任【摘要】随着三维深部地质填图和含页岩气盆地潜力的评价项目展开,大地电磁方法的重要性越来越突显.大地电磁资料处理和解释中,至今还存在几个挑战性难题,例如阻抗张量旋转、TE和TM模式识别、静态矫正、TE和TM模式选择二维反演问题•从WinGlink软件大地电磁模块出发，使用5种阻抗旋转方法和3种数据圆滑方法处理在内蒙古扎鲁特旗采集的大地电磁数据.以数据质量不同的3个测点的处理结果为依据,对各种旋转方法和圆滑方法进行详细评价,数据极化图显示地下结构模型是复杂的三维结构•对这些方法的比较研究，对三维结构下大地电磁数据处理，特别是对使用WinGlink软件处理大地电磁数据具有重要现实意义.期刊名称】《物探与化探》年(卷),期】2013(037)006【总页数】8页(P1125-1132)【关键词】大地电磁;阻抗张量旋转;曲线圆滑；WinGlink软件【作者】赵维俊;孙中任【作者单位】沈阳地质矿产研究所,辽宁沈阳110034;沈阳地质矿产研究所,辽宁沈阳110034【正文语种】中文【中图分类】P631大地电磁测深［1-2］是一种被动源电磁测深方法。它利用在地表测量随着时间变化的电场和磁场来推断地下电阻率分布的电磁方法。大地电磁(MT)方法在石油天然气勘探［3］、矿产勘查［4］和地热勘查［5］中得到了广泛应用。近年来，中国地质调查局开展了三维地质填图和深部地质调查试点工作，MT这种球物理工作的重要作用越来越突出，而随着对能源需求的日益增加，非常规能源页岩气的勘探与开发已经提上日程，MT在页岩气勘探中也必将起到重要作用。但是，在MT资料处理与反演中，还存在一些难题没有很好地解决，例如大地电磁数据二维反演的数据旋转方向选择［6］、反演数据极化模式选择［7-8］等。在中国地质调查局试点项目“松辽盆地外围深部地质调查”支持下，2012年在内蒙古扎鲁特地区采集了5条剖面的大地电磁数据。笔者利用实测大地电磁数据研究阻抗张量旋转和数据圆滑问题。从WinGlink(WG)软件大地电磁模块出发，使用5种阻抗旋转方法和3种数据圆滑方法处理采集大地电磁数据。选用数据质量不同3点处理结果为依据，对各种旋转方法和圆滑方法进行详细评价。又选择某条剖面连续5点用5种阻抗旋转方法来处理，对处理结果进行详细分析和评价。最后，实验选择的3点数据极化图显示地下结构模型是复杂三维结构。5条MT剖面总长495km，测量仪器采用PhoenixGeophysics公司的V8多功能电磁仪，测量点距为1km，有效频率范围0.01~320Hz,Ex电场方向为测线方向。时间序列MT数据经过SSMT2000和MTEDITOR预处理后输出EDI文件。在实测MT数据中选择质量不同的3个测点作方法对比研究，它们在所有数据中具有一定代表性。这三点数据质量仅仅从视电阻率曲线和阻抗相位的连续性和趋势上就可作判断，然而它们的维数(dimensionality)或者说它们的地下模型是怎样的呢。知道地下模型的维数对数据处理结果判断更有利，对于如何获得测点地下模型的维数，通常是从离散度(swiftskew)、扁率(ellipticity)、倾子振幅(tippermagnitude)和极化图(polardiagram)来判断分析。下面选择最能体现维数的极化图来讨论。图1给出了3个测点的极化图。110点Zxy的极化图为花生状，Zxx极化图为四瓣花状。从整个频率范围看，110点地下为三维结构的，极化形状一致，阻抗走向基本相同。60点情况与110截然不同，79~320Hz是二维的，6.9~57Hz是一维，0.59-4.7Hz是二维的，0.073~0.43Hz是复杂的三维，Zxy极化图显示4个花瓣形。0.017-0.054Hz是三维的。所以60点地下模型随着频率变化而变化，从二维变到一维，再变到二维，然后从复杂三维到三维变化。40点情况也比较复杂，57-320Hz是二维的，1.17-40Hz是三维的，0.037-0.86Hz是二维的，0.0017-0.0269Hz是三维的。总之，110点数据是一致的三维的，60点是一、二、三维变换的，40点数据是二、三维变换的。2.1阻抗张量在Fourier频率域里，复数阻抗张量Z是把地球表面测量的天然磁场H映射到天然电场E。它的计算公式为视电阻率pxy、pyx和阻抗相位申xy、申yx,通过阻抗张量，由计算。式中，3是角频率，p0是真空磁导率。阻抗张量准确估计是非常重要的，因为视电阻率和相位都是由它导出的。阻抗张量包含地质结构重要信息，是一个重要研究对象。研究阻抗张量的目的是确定地质模型维数和方向性及其随周期增大的变化特征。阻抗张量是地质模型的响应函数［9］。地质模型是一维均匀层状模型，阻抗张量Z1D为xy与yx视电阻率重合在一起；xy与yx阻抗相位相差180°。地质模型是二维时，观测的Z2D为式中，a是地质构造走向顺时针旋转到观测电场Ex方向的角度。后面讨论阻抗旋转，就是把观测电场Ex顺时针旋转到地质构造走向上。二维阻抗旋转方程产生几个旋转不变量［1-2］。地质模型是三维时，阻抗张量是普通二维复数矩阵。当地质模型是二维区域构造叠加局部三维畸变时，Bahr分解［10-12］、GB分解［13-14］方法可以用来判断构造维数。Martin等利用Bahr旋转不变量和WAL旋转不变量编写出程序WALDIM来分析构造维数［15-16］。Utada和Munekane提出没有区域构造二维的假设新阻抗分解方法［17］。Caldwell提出在纯三维地质模型下phasetensor分解方法［18］。蔡军涛等基于phasetensor方法提出共轭变换阻抗法［19］。在MT商业软件中没有这些阻抗分解模块，可能这些方法并没有适应大部分实际测量数据的三维复杂地质模型，然而这些方法确实对数据处理和解释有一定指导作用。2.2阻抗张量旋转在数据输入到WG软件MT模块中，或者MTSounding模块中，都有阻抗张量旋转功能。为什么阻抗张量需要旋转？二维MT反演程序是针对二维地质模型的。反演之前，TE和TM极化模式需要被识别出来。而TE和TM模式要被识别出来，阻抗张量就必须旋转到地质构造走向或者说地电结构的走向。如果地质模型是二维，阻抗张量直接旋转到地质构造走向上就可以。WG软件提供5种旋转方法。第一种方法是旋转到主轴方法(PrincipleAxis),阻抗主轴基本就是阻抗走向的含义。阻抗走向(impedancestrike)有4种定义［20］:①最大化|Zxy+Zyx|2:②最大化|Zxy|2+|Zyx|2:③最小化|Zxx|2+|Zyy|2:④最小化|Zxx+Zyy|2。对于理想二维地质模型，这四种方法得到同样阻抗走向角度，在三维地质模型下，它们得到的角度有一些差别，WG软件里没有明确指出采用哪种计算方法。第二种是旋转到用户定义角度方法(UserDefined)，就是旋转到用户规定的一个角度。第三种是旋转到最大Zxy方法(ZxyMax)，这种方法旋转阻抗张量到最大Zxy方向上，最大Zxy的一组曲线为TE视电阻率曲线。第四种是旋转到最小Zxy方法(ZxyMin)，这种方法旋转阻抗张量到最小Zxy方向上，最小Zxy的一组曲线为TE视电阻率曲线。第五种是LaTorraca分解方法（LaTorracaDec）[21],它不是真正的旋转，而是一种分解方法。以下图2~图11中，红色曲线代表xy分量，软件默认TE模式；蓝色曲线代表yx分量，软件默认TM模式。选择数据质量较好的110点数据，用WG提供的5种阻抗张量旋转方法进行处理。处理结果显示（图2）,最小Zxy方法与最大Zxy方法相反除外，其余4个方法得到的结果基本一致。TE模式视电阻率曲线在TM视电阻率之上，随着周期增大，TM视电阻率加速下降。然而用户定义方法的TE和TM模式视电阻率之差最小，特别是在小周期。PrincipleAxis,ZxyMax和LaTorracaDec方法的结果最相似。40点的数据在0.1Hz以下出现跳跃，没有趋势，数据质量较差。处理结果显示（图3）,PrincipleAxis、UserDefined和ZxyMin方法的处理结果类似，对TM模式的处理结果要优于TE模式，ZxyMax和LaTorracaDec的处理结果与此相反。60点的数据质量中等，其视电阻率曲线在10s周期左右突然向上跳跃,100s开始跳跃严重。处理结果显示（图4）,除ZxyMin方法，其余4种方法处理结果基本类似，TE模式的处理结果好于TM模式的。然而仔细观察发现，PrincipleAxis和UserDefined方法结果非常相似，ZxyMax和LaTorracaDec方法结果非常相似,它们对处理结果的差异表现在两条视电阻率曲线相交的前部分。质量不同的3点数据处理结果显示数据质量好坏严重影响阻抗张量旋转结果,影响TE和TM模式视电阻率的识别。数据质量越好，ZxyMin之外的4种阻抗张量旋转方法处理结果越一致，TE和TM模式视电祖率曲线越容易识别。ZxyMin方法与ZxyMax方法相反，在实际处理中可以单独使用ZxyMax方法，只有在特别情况下才需要用到ZxyMin方法。这3点数据不能说明用哪个旋转方法好，但是从处理结果观察，它们都分成处理结果相反的两组。单点MT数据旋转不能说明更多问题，因为在相同地质结构下，连续几个测点应该具有类似的曲线形态，旋转后结果应该类似，所以选用一条测线上连续5点来作阻抗张量旋转，采用除ZxyMin外4种旋转方法来比较它们的旋转结果。图5~图8显示这连续5点不是在同一个地质结构下，把30、40点归为A组，50、60、70点归为B组，A组和B组处于两个不同地质结构下。A组里PrincipleAxis和UserDefined方法处理结果类似，而ZxyMax和LaTorracaDec方法处理结果相反。B组里PrincipleAxis方法处理结果显示，50点和60点结果类似，而70点结果相反，所以PrincipleAxis处理识别出错。为校正PrincipleAxis旋转结果，可以使用UserDefined90°来处理，或者通过模式交换(ModeSwap)功能实现交换。UserDefined,ZxyMax,LaTorracaDec处理结果基本一致，TE模式视电阻率曲线绝大部分在TM模式视电阻率曲线之上。然而在50点上，只有UserDefined方法旋转后两条视电阻率曲线相交。B组的60、70点在这4种旋转方法下，两条视电阻率曲线都相交。在WG大地电磁模块中，有三种视电阻率和相位曲线圆滑方法。它们是Sutarno相位一致性方法［22］(Sutarnophaseconsistentsmoothing),D+函数圆滑方法［23-24］(D+functionsmoothing)和数值圆滑方法［25］(Numericalsmoothing)。这些圆滑过程实际是一种滤波。Sutarno相位一致性方法是通过假设阻抗函数是最小相位，对相位做Hilbert变换得到阻抗振幅。它的基本思想是：标准分析技术对相位函数能作一致性估计，而不能对振幅函数作一致性估计。这种方法主要用来确认视电阻率和相位是一致的。MT响应的最小相位性质在1D和2D结构中能预测出来，在3D结构中没有被证实［26］。D+函数圆滑方法是把同一分量(xy或yx)的视电阻率和相位通过D+函数联系起来，本质上是找到一个一维D+函数地质模型来最佳匹配视电阻率和相位，这个函数适用于大部分2D模型和部分3D模型。数值圆滑方法是以快速傅里叶变换为基础的低通滤波，这个滤波独立计算4个分量中每一个分量的光滑曲线。采用三种圆滑方法对数据质量不同的3个点作处理分析。110点的圆滑结果（图9）显示，Sutarno相位一致性方法对TM视电阻率曲线圆滑效果不好，圆滑曲线和实测曲线不匹配，但对相位的圆滑效果很好；D+函数和数值圆滑对视电阻率和相位都有很好的圆滑效果。数据质量较差的60点圆滑结果（图10）显示，3种方法对相位曲线的圆滑效果好，但Sutarno方法在2s周期后，视电阻率圆滑曲线和原始曲线相差很大；D+函数圆滑与数值圆滑结果类似，但是数值圆滑的结果更逼近实测曲线。数据质量最差的40点处理结果（图11）中，Sutarno方法相位圆滑效果好，视电阻率曲线圆滑在1s后出现畸变；D+函数圆滑方法对TM模式相位圆滑结果产生很大误差，对视电阻率曲线圆滑效果比较好；数值圆滑方法对视电阻率曲线和相位均有较好的圆滑效果。上述结果显示，数值圆滑方法最逼近原始曲线，D+函数方法次之，Sutarn。相位—致性方法最差。实际处理MT数据经验也显示Sutarno方法圆滑的结果较差，原因可能为Sutarno方法假设阻抗函数是最小相位，实际数据很多并不满足这个要求。D+函数圆滑方法假设是严格的一维地质模型，用一些深度delta函数来逼近实测曲线。数值圆滑方法处理结果最逼近原始曲线，因为它没有任何物理假设前提，直接做Fourier低通滤波。D+函数方法在实验的3点数据还是没有数值光滑方法好。对WG大地电磁模块中的5种阻抗张量旋转方法作旋转对比。ZxyMin方法排除，因为它只在少数情况下使用，在数据质量较好的情况下，其余4种旋转方法处理结果类似，若数据质量较差，则这4种旋转方法处理结果通常分为相反的两组。LaTorracaDec方法实质不是旋转方法，而是一种分解方法，如果没有特殊原因，谨慎使用。UserDefined方法具有灵活性，可以按照输入的任意角度进行旋转。不论选用哪一种旋转方法，也需参照同一条剖面连续测点的旋转结果。通常WG软件默认使用PrincipleAxis方法，因为在理想二维地质模型下，实测数据必须旋转到主轴上来识别TE和TM模式，然后才能进行TE和TM反演。在三维地质模型情况下，使用哪种旋转方法要具体情况具体分析。对于WG大地电磁模块中的3种圆滑方法，数据质量越好，3种圆滑方法处理结果越好。Sutarno相位一致性圆滑方法在数据质量较好的情况下，其圆滑效果不好，因此在三维复杂地质情况下，强烈建议不使用该方法；D+函数圆滑方法在数据质量较差的情况下不建议使用，此情况下建议使用数值圆滑方法。【相关文献】SimpsonF，BahrK.Practicalmagnetotellurics[M].CambridgeU-niversityPress，2005.BerdichevskyMN，DmitrievVI.Modelsandmethodsofmagnetotellurics[M].SpringerPress，2008.VozoffK.Themagnetotelluricmethodintheexplorationofsedimentarybasins[J].Geophysics，1972，37（1）：98-141.LivelybrooksD，MareschalM，BlaisE，etal.MagnetotelluricdelineationoftheTrillabellemassivesulfidebodyinSudbury，Ontario[J].Geophysics，1996，61（4）：971-986.StanleyWD，BoehlJE，BostickFX，etal.GeothermalsignificanceofmagnetotelluricsoundingintheEasternSnakeRiverPlain-YellowstoneRegion[J].JournalofGeophysicalResearch，1977，82（17）：2501-2514.[6]陈小斌，赵国泽，马宵•关于MT二维反演中数据旋转方向的选择问题初探[J].石油地球物理勘探，2008，43（1）：113-128.[7]蔡军涛，陈小斌.大地电磁资料精细处理和二维反演解释技术研究（二）：反演数据极化模式选择[J』•地球物理学报，2010,53（11）:2703-2714.[8]贺春艳，郭秋峰•应用Winglink进行大地电磁测深极化模式识别的研究[G]//山东地球物理六十年.青岛:中国海洋大学出版社，2010:557-564.ChristophersonKR，JonesA，MackieR.Magnetotelluricsfornaturalresourcesfromacquisitionthroughinterpretation[M].SEGpublishing，2002.BahrK.Interpretationofthemagnetotelluricimpedancetensor：regionalinductionandlocaltelluricdistortion[J].JournalofGeophysics，1988，62：119-127.BahrK.Geologicalnoiseinmagnetotelluricdata：aclassificationofdistortiontype[J].PhysicsofEarthplanetinterior，1991，66：24-38.[12]PracserE,SzarkaL.AcorrectiontoBahr's“phasedeviation“methodfortensordecomposition[J].EarthPlanetsSpace，1999，51：1019-1022.GroomRW，BaileyRC.Decompositionofmagnetotelluricimpedancetensorinthepresenceoflocalthree-dimensionalgalvanicdistortion[J].JournalofGeophysicalResearch，1989，94（B2）：1913-1925.GroomRW，BaileyRC.Analyticinvestigationsoftheeffectsofnear-surfacethree-dimensionalgalvanicscatterersonMTtensordecompositions[J].Geophysics，1991，56（4）：496-518.MatiA，QueraltP，JoesAG，etal.ImprovingBahr'sinvariantparametersusingtheWALapproach[J].GeophysicalJournalInternational，2005，163（1）：38-41.MartiA，QueraltP，LedoJ.WALDIM：AcodeforthedimensionalityanalysisofmagnetotelluricdatausingtherotationalinvariantoftheMagnetotellurics[J].ComputersandGeosciences，2009，35：2295-2303.UtadaH，MunekaneH.Ongalvanicdistortionofregionalthree-dimensionalmagnetotelluricimpedance[J].GeophysicalJournalInternational，2000，140（2）：385-398.CaldwellTG，BibbyHM，BrownC.Themagnetotelluricphasetensor[J].GeophysicalJournalInternational，2004，158（2）：457-469.[19]蔡军涛，陈小斌，赵国泽.大地电磁资料精细处理和二维反演解释技术研究（一）：阻抗张量分解与构造维性分析[J]•地球物理学报，2010,53（10）:2516-2526.MT/EMAPDataInterchange