大地電磁阻抗張量旋轉(zhuǎn)方法和曲線圓滑方法的比較

大地電磁阻抗張量旋轉(zhuǎn)方法和曲線圓滑方法的比較趙維俊;孫中任【摘要】隨著三維深部地質(zhì)填圖和含頁巖氣盆地潛力的評價(jià)項(xiàng)目展開,大地電磁方法的重要性越來越突顯.大地電磁資料處理和解釋中,至今還存在幾個挑戰(zhàn)性難題,例如阻抗張量旋轉(zhuǎn)、TE和TM模式識別、靜態(tài)矯正、TE和TM模式選擇二維反演問題?從WinGlink軟件大地電磁模塊出發(fā)，使用5種阻抗旋轉(zhuǎn)方法和3種數(shù)據(jù)圓滑方法處理在內(nèi)蒙古扎魯特旗采集的大地電磁數(shù)據(jù).以數(shù)據(jù)質(zhì)量不同的3個測點(diǎn)的處理結(jié)果為依據(jù),對各種旋轉(zhuǎn)方法和圓滑方法進(jìn)行詳細(xì)評價(jià),數(shù)據(jù)極化圖顯示地下結(jié)構(gòu)模型是復(fù)雜的三維結(jié)構(gòu)?對這些方法的比較研究，對三維結(jié)構(gòu)下大地電磁數(shù)據(jù)處理，特別是對使用WinGlink軟件處理大地電磁數(shù)據(jù)具有重要現(xiàn)實(shí)意義.期刊名稱】《物探與化探》年(卷),期】2013(037)006【總頁數(shù)】8頁(P1125-1132)【關(guān)鍵詞】大地電磁;阻抗張量旋轉(zhuǎn);曲線圓滑；WinGlink軟件【作者】趙維俊;孫中任【作者單位】沈陽地質(zhì)礦產(chǎn)研究所,遼寧沈陽110034;沈陽地質(zhì)礦產(chǎn)研究所,遼寧沈陽110034【正文語種】中文【中圖分類】P631大地電磁測深［1-2］是一種被動源電磁測深方法。它利用在地表測量隨著時(shí)間變化的電場和磁場來推斷地下電阻率分布的電磁方法。大地電磁(MT)方法在石油天然氣勘探［3］、礦產(chǎn)勘查［4］和地?zé)峥辈椋?］中得到了廣泛應(yīng)用。近年來，中國地質(zhì)調(diào)查局開展了三維地質(zhì)填圖和深部地質(zhì)調(diào)查試點(diǎn)工作，MT這種球物理工作的重要作用越來越突出，而隨著對能源需求的日益增加，非常規(guī)能源頁巖氣的勘探與開發(fā)已經(jīng)提上日程，MT在頁巖氣勘探中也必將起到重要作用。但是，在MT資料處理與反演中，還存在一些難題沒有很好地解決，例如大地電磁數(shù)據(jù)二維反演的數(shù)據(jù)旋轉(zhuǎn)方向選擇［6］、反演數(shù)據(jù)極化模式選擇［7-8］等。在中國地質(zhì)調(diào)查局試點(diǎn)項(xiàng)目“松遼盆地外圍深部地質(zhì)調(diào)查”支持下，2012年在內(nèi)蒙古扎魯特地區(qū)采集了5條剖面的大地電磁數(shù)據(jù)。筆者利用實(shí)測大地電磁數(shù)據(jù)研究阻抗張量旋轉(zhuǎn)和數(shù)據(jù)圓滑問題。從WinGlink(WG)軟件大地電磁模塊出發(fā)，使用5種阻抗旋轉(zhuǎn)方法和3種數(shù)據(jù)圓滑方法處理采集大地電磁數(shù)據(jù)。選用數(shù)據(jù)質(zhì)量不同3點(diǎn)處理結(jié)果為依據(jù)，對各種旋轉(zhuǎn)方法和圓滑方法進(jìn)行詳細(xì)評價(jià)。又選擇某條剖面連續(xù)5點(diǎn)用5種阻抗旋轉(zhuǎn)方法來處理，對處理結(jié)果進(jìn)行詳細(xì)分析和評價(jià)。最后，實(shí)驗(yàn)選擇的3點(diǎn)數(shù)據(jù)極化圖顯示地下結(jié)構(gòu)模型是復(fù)雜三維結(jié)構(gòu)。5條MT剖面總長495km，測量儀器采用PhoenixGeophysics公司的V8多功能電磁儀，測量點(diǎn)距為1km，有效頻率范圍0.01~320Hz,Ex電場方向?yàn)闇y線方向。時(shí)間序列MT數(shù)據(jù)經(jīng)過SSMT2000和MTEDITOR預(yù)處理后輸出EDI文件。在實(shí)測MT數(shù)據(jù)中選擇質(zhì)量不同的3個測點(diǎn)作方法對比研究，它們在所有數(shù)據(jù)中具有一定代表性。這三點(diǎn)數(shù)據(jù)質(zhì)量僅僅從視電阻率曲線和阻抗相位的連續(xù)性和趨勢上就可作判斷，然而它們的維數(shù)(dimensionality)或者說它們的地下模型是怎樣的呢。知道地下模型的維數(shù)對數(shù)據(jù)處理結(jié)果判斷更有利，對于如何獲得測點(diǎn)地下模型的維數(shù)，通常是從離散度(swiftskew)、扁率(ellipticity)、傾子振幅(tippermagnitude)和極化圖(polardiagram)來判斷分析。下面選擇最能體現(xiàn)維數(shù)的極化圖來討論。圖1給出了3個測點(diǎn)的極化圖。110點(diǎn)Zxy的極化圖為花生狀，Zxx極化圖為四瓣花狀。從整個頻率范圍看，110點(diǎn)地下為三維結(jié)構(gòu)的，極化形狀一致，阻抗走向基本相同。60點(diǎn)情況與110截然不同，79~320Hz是二維的，6.9~57Hz是一維，0.59-4.7Hz是二維的，0.073~0.43Hz是復(fù)雜的三維，Zxy極化圖顯示4個花瓣形。0.017-0.054Hz是三維的。所以60點(diǎn)地下模型隨著頻率變化而變化，從二維變到一維，再變到二維，然后從復(fù)雜三維到三維變化。40點(diǎn)情況也比較復(fù)雜，57-320Hz是二維的，1.17-40Hz是三維的，0.037-0.86Hz是二維的，0.0017-0.0269Hz是三維的?？傊?，110點(diǎn)數(shù)據(jù)是一致的三維的，60點(diǎn)是一、二、三維變換的，40點(diǎn)數(shù)據(jù)是二、三維變換的。2.1阻抗張量在Fourier頻率域里，復(fù)數(shù)阻抗張量Z是把地球表面測量的天然磁場H映射到天然電場E。它的計(jì)算公式為視電阻率pxy、pyx和阻抗相位申xy、申yx,通過阻抗張量，由計(jì)算。式中，3是角頻率，p0是真空磁導(dǎo)率。阻抗張量準(zhǔn)確估計(jì)是非常重要的，因?yàn)橐曤娮杪屎拖辔欢际怯伤鼘?dǎo)出的。阻抗張量包含地質(zhì)結(jié)構(gòu)重要信息，是一個重要研究對象。研究阻抗張量的目的是確定地質(zhì)模型維數(shù)和方向性及其隨周期增大的變化特征。阻抗張量是地質(zhì)模型的響應(yīng)函數(shù)［9］。地質(zhì)模型是一維均勻?qū)訝钅Ｐ?，阻抗張量Z1D為xy與yx視電阻率重合在一起；xy與yx阻抗相位相差180°。地質(zhì)模型是二維時(shí)，觀測的Z2D為式中，a是地質(zhì)構(gòu)造走向順時(shí)針旋轉(zhuǎn)到觀測電場Ex方向的角度。后面討論阻抗旋轉(zhuǎn)，就是把觀測電場Ex順時(shí)針旋轉(zhuǎn)到地質(zhì)構(gòu)造走向上。二維阻抗旋轉(zhuǎn)方程產(chǎn)生幾個旋轉(zhuǎn)不變量［1-2］。地質(zhì)模型是三維時(shí)，阻抗張量是普通二維復(fù)數(shù)矩陣。當(dāng)?shù)刭|(zhì)模型是二維區(qū)域構(gòu)造疊加局部三維畸變時(shí)，Bahr分解［10-12］、GB分解［13-14］方法可以用來判斷構(gòu)造維數(shù)。Martin等利用Bahr旋轉(zhuǎn)不變量和WAL旋轉(zhuǎn)不變量編寫出程序WALDIM來分析構(gòu)造維數(shù)［15-16］。Utada和Munekane提出沒有區(qū)域構(gòu)造二維的假設(shè)新阻抗分解方法［17］。Caldwell提出在純?nèi)S地質(zhì)模型下phasetensor分解方法［18］。蔡軍濤等基于phasetensor方法提出共軛變換阻抗法［19］。在MT商業(yè)軟件中沒有這些阻抗分解模塊，可能這些方法并沒有適應(yīng)大部分實(shí)際測量數(shù)據(jù)的三維復(fù)雜地質(zhì)模型，然而這些方法確實(shí)對數(shù)據(jù)處理和解釋有一定指導(dǎo)作用。2.2阻抗張量旋轉(zhuǎn)在數(shù)據(jù)輸入到WG軟件MT模塊中，或者M(jìn)TSounding模塊中，都有阻抗張量旋轉(zhuǎn)功能。為什么阻抗張量需要旋轉(zhuǎn)？二維MT反演程序是針對二維地質(zhì)模型的。反演之前，TE和TM極化模式需要被識別出來。而TE和TM模式要被識別出來，阻抗張量就必須旋轉(zhuǎn)到地質(zhì)構(gòu)造走向或者說地電結(jié)構(gòu)的走向。如果地質(zhì)模型是二維，阻抗張量直接旋轉(zhuǎn)到地質(zhì)構(gòu)造走向上就可以。WG軟件提供5種旋轉(zhuǎn)方法。第一種方法是旋轉(zhuǎn)到主軸方法(PrincipleAxis),阻抗主軸基本就是阻抗走向的含義。阻抗走向(impedancestrike)有4種定義［20］:①最大化|Zxy+Zyx|2:②最大化|Zxy|2+|Zyx|2:③最小化|Zxx|2+|Zyy|2:④最小化|Zxx+Zyy|2。對于理想二維地質(zhì)模型，這四種方法得到同樣阻抗走向角度，在三維地質(zhì)模型下，它們得到的角度有一些差別，WG軟件里沒有明確指出采用哪種計(jì)算方法。第二種是旋轉(zhuǎn)到用戶定義角度方法(UserDefined)，就是旋轉(zhuǎn)到用戶規(guī)定的一個角度。第三種是旋轉(zhuǎn)到最大Zxy方法(ZxyMax)，這種方法旋轉(zhuǎn)阻抗張量到最大Zxy方向上，最大Zxy的一組曲線為TE視電阻率曲線。第四種是旋轉(zhuǎn)到最小Zxy方法(ZxyMin)，這種方法旋轉(zhuǎn)阻抗張量到最小Zxy方向上，最小Zxy的一組曲線為TE視電阻率曲線。第五種是LaTorraca分解方法（LaTorracaDec）[21],它不是真正的旋轉(zhuǎn)，而是一種分解方法。以下圖2~圖11中，紅色曲線代表xy分量，軟件默認(rèn)TE模式；藍(lán)色曲線代表yx分量，軟件默認(rèn)TM模式。選擇數(shù)據(jù)質(zhì)量較好的110點(diǎn)數(shù)據(jù)，用WG提供的5種阻抗張量旋轉(zhuǎn)方法進(jìn)行處理。處理結(jié)果顯示（圖2）,最小Zxy方法與最大Zxy方法相反除外，其余4個方法得到的結(jié)果基本一致。TE模式視電阻率曲線在TM視電阻率之上，隨著周期增大，TM視電阻率加速下降。然而用戶定義方法的TE和TM模式視電阻率之差最小，特別是在小周期。PrincipleAxis,ZxyMax和LaTorracaDec方法的結(jié)果最相似。40點(diǎn)的數(shù)據(jù)在0.1Hz以下出現(xiàn)跳躍，沒有趨勢，數(shù)據(jù)質(zhì)量較差。處理結(jié)果顯示（圖3）,PrincipleAxis、UserDefined和ZxyMin方法的處理結(jié)果類似，對TM模式的處理結(jié)果要優(yōu)于TE模式，ZxyMax和LaTorracaDec的處理結(jié)果與此相反。60點(diǎn)的數(shù)據(jù)質(zhì)量中等，其視電阻率曲線在10s周期左右突然向上跳躍,100s開始跳躍嚴(yán)重。處理結(jié)果顯示（圖4）,除ZxyMin方法，其余4種方法處理結(jié)果基本類似，TE模式的處理結(jié)果好于TM模式的。然而仔細(xì)觀察發(fā)現(xiàn)，PrincipleAxis和UserDefined方法結(jié)果非常相似，ZxyMax和LaTorracaDec方法結(jié)果非常相似,它們對處理結(jié)果的差異表現(xiàn)在兩條視電阻率曲線相交的前部分。質(zhì)量不同的3點(diǎn)數(shù)據(jù)處理結(jié)果顯示數(shù)據(jù)質(zhì)量好壞嚴(yán)重影響阻抗張量旋轉(zhuǎn)結(jié)果,影響TE和TM模式視電阻率的識別。數(shù)據(jù)質(zhì)量越好，ZxyMin之外的4種阻抗張量旋轉(zhuǎn)方法處理結(jié)果越一致，TE和TM模式視電祖率曲線越容易識別。ZxyMin方法與ZxyMax方法相反，在實(shí)際處理中可以單獨(dú)使用ZxyMax方法，只有在特別情況下才需要用到ZxyMin方法。這3點(diǎn)數(shù)據(jù)不能說明用哪個旋轉(zhuǎn)方法好，但是從處理結(jié)果觀察，它們都分成處理結(jié)果相反的兩組。單點(diǎn)MT數(shù)據(jù)旋轉(zhuǎn)不能說明更多問題，因?yàn)樵谙嗤刭|(zhì)結(jié)構(gòu)下，連續(xù)幾個測點(diǎn)應(yīng)該具有類似的曲線形態(tài)，旋轉(zhuǎn)后結(jié)果應(yīng)該類似，所以選用一條測線上連續(xù)5點(diǎn)來作阻抗張量旋轉(zhuǎn)，采用除ZxyMin外4種旋轉(zhuǎn)方法來比較它們的旋轉(zhuǎn)結(jié)果。圖5~圖8顯示這連續(xù)5點(diǎn)不是在同一個地質(zhì)結(jié)構(gòu)下，把30、40點(diǎn)歸為A組，50、60、70點(diǎn)歸為B組，A組和B組處于兩個不同地質(zhì)結(jié)構(gòu)下。A組里PrincipleAxis和UserDefined方法處理結(jié)果類似，而ZxyMax和LaTorracaDec方法處理結(jié)果相反。B組里PrincipleAxis方法處理結(jié)果顯示，50點(diǎn)和60點(diǎn)結(jié)果類似，而70點(diǎn)結(jié)果相反，所以PrincipleAxis處理識別出錯。為校正PrincipleAxis旋轉(zhuǎn)結(jié)果，可以使用UserDefined90°來處理，或者通過模式交換(ModeSwap)功能實(shí)現(xiàn)交換。UserDefined,ZxyMax,LaTorracaDec處理結(jié)果基本一致，TE模式視電阻率曲線絕大部分在TM模式視電阻率曲線之上。然而在50點(diǎn)上，只有UserDefined方法旋轉(zhuǎn)后兩條視電阻率曲線相交。B組的60、70點(diǎn)在這4種旋轉(zhuǎn)方法下，兩條視電阻率曲線都相交。在WG大地電磁模塊中，有三種視電阻率和相位曲線圓滑方法。它們是Sutarno相位一致性方法［22］(Sutarnophaseconsistentsmoothing),D+函數(shù)圓滑方法［23-24］(D+functionsmoothing)和數(shù)值圓滑方法［25］(Numericalsmoothing)。這些圓滑過程實(shí)際是一種濾波。Sutarno相位一致性方法是通過假設(shè)阻抗函數(shù)是最小相位，對相位做Hilbert變換得到阻抗振幅。它的基本思想是：標(biāo)準(zhǔn)分析技術(shù)對相位函數(shù)能作一致性估計(jì)，而不能對振幅函數(shù)作一致性估計(jì)。這種方法主要用來確認(rèn)視電阻率和相位是一致的。MT響應(yīng)的最小相位性質(zhì)在1D和2D結(jié)構(gòu)中能預(yù)測出來，在3D結(jié)構(gòu)中沒有被證實(shí)［26］。D+函數(shù)圓滑方法是把同一分量(xy或yx)的視電阻率和相位通過D+函數(shù)聯(lián)系起來，本質(zhì)上是找到一個一維D+函數(shù)地質(zhì)模型來最佳匹配視電阻率和相位，這個函數(shù)適用于大部分2D模型和部分3D模型。數(shù)值圓滑方法是以快速傅里葉變換為基礎(chǔ)的低通濾波，這個濾波獨(dú)立計(jì)算4個分量中每一個分量的光滑曲線。采用三種圓滑方法對數(shù)據(jù)質(zhì)量不同的3個點(diǎn)作處理分析。110點(diǎn)的圓滑結(jié)果（圖9）顯示，Sutarno相位一致性方法對TM視電阻率曲線圓滑效果不好，圓滑曲線和實(shí)測曲線不匹配，但對相位的圓滑效果很好；D+函數(shù)和數(shù)值圓滑對視電阻率和相位都有很好的圓滑效果。數(shù)據(jù)質(zhì)量較差的60點(diǎn)圓滑結(jié)果（圖10）顯示，3種方法對相位曲線的圓滑效果好，但Sutarno方法在2s周期后，視電阻率圓滑曲線和原始曲線相差很大；D+函數(shù)圓滑與數(shù)值圓滑結(jié)果類似，但是數(shù)值圓滑的結(jié)果更逼近實(shí)測曲線。數(shù)據(jù)質(zhì)量最差的40點(diǎn)處理結(jié)果（圖11）中，Sutarno方法相位圓滑效果好，視電阻率曲線圓滑在1s后出現(xiàn)畸變；D+函數(shù)圓滑方法對TM模式相位圓滑結(jié)果產(chǎn)生很大誤差，對視電阻率曲線圓滑效果比較好；數(shù)值圓滑方法對視電阻率曲線和相位均有較好的圓滑效果。上述結(jié)果顯示，數(shù)值圓滑方法最逼近原始曲線，D+函數(shù)方法次之，Sutarn。相位—致性方法最差。實(shí)際處理MT數(shù)據(jù)經(jīng)驗(yàn)也顯示Sutarno方法圓滑的結(jié)果較差，原因可能為Sutarno方法假設(shè)阻抗函數(shù)是最小相位，實(shí)際數(shù)據(jù)很多并不滿足這個要求。D+函數(shù)圓滑方法假設(shè)是嚴(yán)格的一維地質(zhì)模型，用一些深度delta函數(shù)來逼近實(shí)測曲線。數(shù)值圓滑方法處理結(jié)果最逼近原始曲線，因?yàn)樗鼪]有任何物理假設(shè)前提，直接做Fourier低通濾波。D+函數(shù)方法在實(shí)驗(yàn)的3點(diǎn)數(shù)據(jù)還是沒有數(shù)值光滑方法好。對WG大地電磁模塊中的5種阻抗張量旋轉(zhuǎn)方法作旋轉(zhuǎn)對比。ZxyMin方法排除，因?yàn)樗辉谏贁?shù)情況下使用，在數(shù)據(jù)質(zhì)量較好的情況下，其余4種旋轉(zhuǎn)方法處理結(jié)果類似，若數(shù)據(jù)質(zhì)量較差，則這4種旋轉(zhuǎn)方法處理結(jié)果通常分為相反的兩組。LaTorracaDec方法實(shí)質(zhì)不是旋轉(zhuǎn)方法，而是一種分解方法，如果沒有特殊原因，謹(jǐn)慎使用。UserDefined方法具有靈活性，可以按照輸入的任意角度進(jìn)行旋轉(zhuǎn)。不論選用哪一種旋轉(zhuǎn)方法，也需參照同一條剖面連續(xù)測點(diǎn)的旋轉(zhuǎn)結(jié)果。通常WG軟件默認(rèn)使用PrincipleAxis方法，因?yàn)樵诶硐攵S地質(zhì)模型下，實(shí)測數(shù)據(jù)必須旋轉(zhuǎn)到主軸上來識別TE和TM模式，然后才能進(jìn)行TE和TM反演。在三維地質(zhì)模型情況下，使用哪種旋轉(zhuǎn)方法要具體情況具體分析。對于WG大地電磁模塊中的3種圓滑方法，數(shù)據(jù)質(zhì)量越好，3種圓滑方法處理結(jié)果越好。Sutarno相位一致性圓滑方法在數(shù)據(jù)質(zhì)量較好的情況下，其圓滑效果不好，因此在三維復(fù)雜地質(zhì)情況下，強(qiáng)烈建議不使用該方法；D+函數(shù)圓滑方法在數(shù)據(jù)質(zhì)量較差的情況下不建議使用，此情況下建議使用數(shù)值圓滑方法?！鞠嚓P(guān)文獻(xiàn)】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]陳小斌，趙國澤，馬宵?關(guān)于MT二維反演中數(shù)據(jù)旋轉(zhuǎn)方向的選擇問題初探[J].石油地球物理勘探，2008，43（1）：113-128.[7]蔡軍濤，陳小斌.大地電磁資料精細(xì)處理和二維反演解釋技術(shù)研究（二）：反演數(shù)據(jù)極化模式選擇[J』?地球物理學(xué)報(bào)，2010,53（11）:2703-2714.[8]賀春艷，郭秋峰?應(yīng)用Winglink進(jìn)行大地電磁測深極化模式識別的研究[G]//山東地球物理六十年.青島:中國海洋大學(xué)出版社，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]蔡軍濤，陳小斌，趙國澤.大地電磁資料精細(xì)處理和二維反演解釋技術(shù)研究（一）：阻抗張量分解與構(gòu)造維性分析[J]?地球物理學(xué)報(bào)，2010,53（10）:2516-2526.MT/EMAPDataInterchange