余才國(guó) 肖 曉*②③④ 湯井田②③④ 皇祥宇 徐錦通

(①中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙 410083；②有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙 410083；③中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙 410083；④自然資源部覆蓋區(qū)深部資源勘查工程技術(shù)創(chuàng)新中心，安徽合肥 230001)

0 引言

大地電磁測(cè)深(MT)通過觀測(cè)地表相互正交的電磁場(chǎng)分量獲取地下電性結(jié)構(gòu)相關(guān)信息，具有不受高阻層屏蔽、對(duì)高導(dǎo)層分辨能力強(qiáng)、橫向電性差異分辨能力較強(qiáng)、勘探深度大、勘探費(fèi)用低、施工方便、資料處理和解釋技術(shù)成熟等優(yōu)勢(shì)，廣泛應(yīng)用于區(qū)域地質(zhì)調(diào)查、資源勘探和工程勘探等領(lǐng)域[1-3]。近年來，隨著研究涉及的地質(zhì)構(gòu)造復(fù)雜程度明顯增加，對(duì)大地電磁測(cè)深資料的處理解釋要求亦更高[4]。當(dāng)近地表存在幾何尺寸遠(yuǎn)遠(yuǎn)小于電磁波趨膚深度的三維電性不均勻小異常體時(shí)，觀測(cè)數(shù)據(jù)會(huì)受到電流型畸變的影響，若不消除這種影響，會(huì)影響后續(xù)的反演解釋[5-10]。

阻抗張量分解方法是目前常用的大地電磁數(shù)據(jù)處理分析方法，可以去除局部畸變帶來的影響，獲得真實(shí)的區(qū)域阻抗、構(gòu)造方位角、構(gòu)造維性等信息。Swift[11]提出了經(jīng)典的阻抗張量分解方法，具體做法是：旋轉(zhuǎn)阻抗張量，尋找阻抗張量反對(duì)角元素最大時(shí)或者主對(duì)角元素最小時(shí)所對(duì)應(yīng)的角度，該角度即為所求的電性主軸方向。但是，此方法獲得的電性主軸方向存在90°的多解性，故需要依據(jù)實(shí)際地質(zhì)信息進(jìn)一步判斷。該方法為后續(xù)的阻抗分解方法奠定了堅(jiān)實(shí)的基礎(chǔ)。Berdichevsky等[12]研究了水平電性非均勻性引起的畸變問題，并將畸變劃分為感應(yīng)型畸變和電流型畸變。Bahr[13]根據(jù)局部畸變的觀測(cè)阻抗的同一列元素相位相等的特性進(jìn)行了畸變校正，提出了Bahr分解方法，但是對(duì)非嚴(yán)格意義上的三維/二維構(gòu)造效果不理想。Bahr[14]在此方法的基礎(chǔ)上提出了Bahr相位偏差法，即阻抗張量同一列元素的阻抗相位存在一個(gè)角度差，這個(gè)角度差稱為相位偏差角，并定義了與維性密切相關(guān)的表征參數(shù)。Groom等[15]提出了經(jīng)典的GB分解法，將局部畸變分為可確定部分(剪切因子和扭曲因子)和不可確定部分(各向異性因子和靜位移因子)，然而GB分解得到的模型是一個(gè)欠定問題：8個(gè)已知參數(shù)(4個(gè)復(fù)數(shù)形式的觀測(cè)阻抗張量元素)，9個(gè)未知參數(shù)(電性主軸角，4個(gè)實(shí)數(shù)形式的畸變矩陣元素，2個(gè)復(fù)數(shù)形式的區(qū)域阻抗張量元素)。需要特別強(qiáng)調(diào)的是，盡管GB分解方法得到了廣泛應(yīng)用，取得了很好的效果，但實(shí)際應(yīng)用中往往忽略了其應(yīng)用前提，即地下介質(zhì)為二維或近似二維的電性結(jié)構(gòu)[16]。Caldwell等[17]提出了相位張量的概念，基于相位張量不受電場(chǎng)畸變的影響這一特點(diǎn)，從觀測(cè)阻抗中提取區(qū)域阻抗的構(gòu)造電性主軸方向以及二維偏度等表征構(gòu)造屬性的參數(shù)，是目前最常用的維性分析方法之一。Bibby等[18]分析了阻抗張量畸變與相位張量的關(guān)系，將相位張量用于局部畸變校正。蔡軍濤等[19]提出一種共軛阻抗變換法(CCZ)，通過對(duì)大地電磁阻抗張量做共軛阻抗變換，求取畸變因子和最佳電性主軸方位角，該方法具有與相位張量分解相似的優(yōu)點(diǎn)，即不受局部電場(chǎng)畸變的影響，且無(wú)需假設(shè)區(qū)域維性。尹曜田等[5]將遺傳算法引入GB分解的求解過程，通過實(shí)際資料處理證明了該算法的有效性。謝成良等[20]使用相位張量分解獲得電性主軸方位角，將旋轉(zhuǎn)后的區(qū)域阻抗作為GB分解區(qū)域阻抗元素的初值，并將各向異性參數(shù)引入目標(biāo)函數(shù)，采用共軛梯度法進(jìn)行模型最優(yōu)化求解，提高了GB算法的穩(wěn)定性及可靠性。

但前人工作大多集中于二維區(qū)域構(gòu)造。對(duì)于三維區(qū)域結(jié)構(gòu)，一些學(xué)者已經(jīng)開展了一些研究。Utada等[21]利用觀測(cè)阻抗張量和觀測(cè)地磁傳遞函數(shù)(即傾子)確定畸變參數(shù)，壓制電流型畸變，獲得三維區(qū)域電性構(gòu)造。Garcia等[22]假設(shè)相鄰測(cè)點(diǎn)的區(qū)域電性結(jié)構(gòu)變化不大但受畸變的影響不同，提出了一種類似于GB分解且適用于三維區(qū)域的畸變校正方法，使用牛頓法求解三維區(qū)域結(jié)構(gòu)中的畸變參數(shù)。Bibby等[18]利用不受電流型畸變影響的相位張量獲得地下維度信息，進(jìn)而采用低維度數(shù)據(jù)求解畸變參數(shù)。這些方法只適用于特定條件，對(duì)于更一般的三維/三維模型不適用。Neukirch等[23]提出了與相位張量互補(bǔ)的振幅張量，利用“相位張量和振幅張量反映同樣的地下介質(zhì)信息，且所有電流型畸變信息都保存在振幅張量中”這一假說，采用遺傳算法實(shí)現(xiàn)了畸變校正，該方法的實(shí)測(cè)效果及理論效果均較好[24]。

本文借鑒Neukirch等[24]的研究成果，基于阻抗振幅張量與相位張量的聯(lián)合分解，采用改進(jìn)粒子群算法求解畸變參數(shù)，對(duì)三維模型正演數(shù)據(jù)和安徽省霍山地區(qū)中生代火山巖盆地實(shí)測(cè)MT數(shù)據(jù)進(jìn)行處理。

1 大地電磁振幅相位張量分解基本理論

1.1 GB分解理論

對(duì)于三維/二維電性結(jié)構(gòu)模型，Groom等[15]只考慮電流型畸變影響，不考慮磁場(chǎng)感應(yīng)型畸變，將區(qū)域構(gòu)造主軸坐標(biāo)系順時(shí)針旋轉(zhuǎn)角度θ到測(cè)量坐標(biāo)系，則觀測(cè)阻抗Zm為

Zm=R(θ)CZ2DRT(θ)

=R(θ)gTSAZ2DRT(θ)

(1)

式中：R(θ)為單位旋轉(zhuǎn)矩陣;C為畸變矩陣；Z2D為區(qū)域阻抗矩陣；g表示增益因子；T表示扭曲角矩陣；S表示剪切角矩陣；A表示各向異性因子矩陣?；兙仃嘋可表示為

(2)

式中：扭曲因子t=tanφt，其中φt∈[-90°,90°]為扭曲角；剪切因子e=tanφe，其中φe∈[-45°,45°]為剪切角；各向異性因子s=tanφs，其中φs∈[-45°,45°]為各向異性因子角。

GB分解把畸變矩陣分解為可確定部分(扭曲角矩陣T、剪切角矩陣S)和不確定部分(各向異性因子矩陣A與增益因子g的乘積)，并將不確定部分并入?yún)^(qū)域阻抗矩陣Z2D，得到

(3)

式中ZTE和ZTM分別表示二維區(qū)域構(gòu)造下TE模式和TM模式的阻抗張量元素。此時(shí)，式(1)變?yōu)?/p>

(4)

上式中右邊8個(gè)參數(shù)(觀測(cè)阻抗的實(shí)部、虛部)為已知，7個(gè)參數(shù)(畸變參數(shù)t、e、θ和區(qū)域阻抗反對(duì)角元素的實(shí)部、虛部)為未知，利用最優(yōu)化求解即可獲得其解?？梢钥闯?，對(duì)于三維介質(zhì)，阻抗張量的區(qū)域阻抗對(duì)角元素中也存在未知數(shù)，未知數(shù)為11個(gè)，故傳統(tǒng)的GB分解無(wú)法適用于三維介質(zhì)的情況。使用GB分解的前提是區(qū)域構(gòu)造為二維，對(duì)于三維區(qū)域構(gòu)造不再適用，因?yàn)槿S區(qū)域構(gòu)造增加了4個(gè)未知數(shù)，此時(shí)的未知數(shù)個(gè)數(shù)大于方程個(gè)數(shù)，傳統(tǒng)的阻抗分解方法失效，需要探索新的方法求解畸變參數(shù)。

1.2 相位張量

大地電磁阻抗張量Z為一個(gè)復(fù)數(shù)張量，既含有電阻率信息，也含有相位信息，可以表示為Z=X+iY，這里X和Y分別表示阻抗張量的實(shí)部和虛部。根據(jù)Caldwell等[17]對(duì)相位張量的定義，相位張量Φ的表達(dá)式為

(5)

式中：Xpq、Ypq分別表示阻抗張量元素的實(shí)部和虛部，其中p=1、2，q=1、2；det(X)=X11X22-X12X21；Φpq表示相位張量Φ的元素。

只考慮電場(chǎng)畸變影響時(shí)，畸變矩陣C是一個(gè)與頻率無(wú)關(guān)的實(shí)數(shù)矩陣，觀測(cè)阻抗Zm可以寫為

Zm=CZ=C(X+iY)

(6)

記

Xm=CX

(7)

Ym=CY

(8)

根據(jù)式(6)，若能求得畸變矩陣C，即可計(jì)算得到未畸變的阻抗張量Z

Z=C-1Zm

(9)

根據(jù)相位張量的定義，此時(shí)的觀測(cè)相位張量Φm可根據(jù)下式求得

(10)

由上式可以看出，相位張量不受局部電場(chǎng)畸變的影響，因此可廣泛應(yīng)用于區(qū)域構(gòu)造走向的研究和維度分析。相位張量概念的提出對(duì)于阻抗張量分解具有里程碑式的意義，為大地電磁數(shù)據(jù)分析提供了新的思路和方法。但同時(shí)也意味著采用相位張量無(wú)法判斷測(cè)點(diǎn)受畸變的程度，因而只采用不含振幅信息的相位張量進(jìn)行反演，會(huì)很大程度增加反演的非唯一性，導(dǎo)致反演不穩(wěn)定。

1.3 振幅張量

通過相位張量的定義發(fā)現(xiàn)，相位張量是阻抗相位的混合，只包括阻抗的相位信息，僅是阻抗張量信息的一半，這對(duì)于測(cè)點(diǎn)本就稀疏、數(shù)據(jù)不足的大地電磁來說，是一項(xiàng)非常大的損失。Neukirch等[23]為了提取振幅信息，創(chuàng)造性地提出了類似于相位張量的振幅張量。

根據(jù)Neukirch等[23]的研究成果，一個(gè)復(fù)數(shù)M可用歐拉公式展開為一個(gè)實(shí)數(shù)L和一個(gè)角度β的函數(shù)

M=Lei β=L(cosβ+isinβ)

(11)

按照式(11)的形式，可將單個(gè)復(fù)數(shù)擴(kuò)展到二階復(fù)數(shù)張量的形式，則阻抗張量Z可以表示為

Z=PeiΦ=P[cos(Φ)+isin(Φ)]

(12)

式中：為了書寫方便，cos(·)和sin(·)表示對(duì)一個(gè)矩陣的各元素分別求余弦和正弦；P為振幅張量。由于采用張量表示形式，下面給出振幅張量的具體求取方法。

張量下的約束條件為

cos(Φ)[cos(Φ)]T+sin(Φ)[sin(Φ)]T=I

(13)

式中I為二階單位矩陣。通過式(13)和式(5)可得

(14)

(15)

eiΦ=cos(Φ)+isin(Φ)

(16)

將式(16)代入式(12)，便可得到振幅張量P的唯一表示形式

(17)

2 振幅—相位張量分解最優(yōu)化畸變校正算法

2.1 目標(biāo)函數(shù)

Booker[25]提出了一種相位張量參數(shù)化分解方法，把相位張量的奇異值分解與區(qū)域構(gòu)造走向、維度等地質(zhì)信息聯(lián)系起來，這種分解方法同樣適用于同為實(shí)數(shù)矩陣的振幅張量。振幅張量和相位張量都是地下介質(zhì)、構(gòu)造信息的反映，在某種程度上其表現(xiàn)形式相似，而且相位張量不受電流型畸變的影響，所有的畸變信息都隱含在振幅張量里，通過尋找與相位張量相似性最高的振幅張量，便可實(shí)現(xiàn)大地電磁的畸變校正。其中，畸變矩陣C為待求參數(shù)，采用式(9)對(duì)觀測(cè)張量阻抗Zm進(jìn)行校正，根據(jù)式(17)獲得振幅張量P，根據(jù)式(5)獲得相位張量Φ，并對(duì)這個(gè)張量進(jìn)行分解。求解畸變參數(shù)C時(shí)，振幅張量P的分解參數(shù)與相位張量Φ的分解參數(shù)最相似。因此，根據(jù)以上思路構(gòu)建了關(guān)于畸變參數(shù)的求解策略，將畸變求解問題轉(zhuǎn)化為一個(gè)最優(yōu)化問題。

相位張量矩陣使用奇異值分解可表示為[25]

(18)

式中VΦ=R(-θΦ)、ΣΦ和WΦ=R(-ψΦ-θΦ)是奇異值分解后的三個(gè)矩陣，其中旋轉(zhuǎn)矩陣R的表達(dá)式為

(19)

θΦ=arccos[Tr(VΦ)/2]表示笛卡爾坐標(biāo)與矩陣坐標(biāo)之間的角度，這里Tr(·)表示求矩陣的跡，在二維情況下，θΦ表示區(qū)域構(gòu)造走向；φ1和φ2是對(duì)角矩陣ΣΦ的奇異值；偏離角ψΦ可表示為

--------------------

(20)

--------------------

同樣地，對(duì)振幅張量進(jìn)行上述分解，可以得到與相位張量分解參數(shù)相對(duì)應(yīng)的各個(gè)分解參數(shù)。

為了使相位張量與振幅張量的奇異值具有可比性，Neukirch等[23]定義了相位相對(duì)各向異性參數(shù)ηφ和振幅相對(duì)各向異性參數(shù)ηP，其表達(dá)式分別為

(21)

(22)

式中P1和P2分別是振幅張量奇異值分解的兩個(gè)奇異值。

假設(shè)某一測(cè)點(diǎn)有N個(gè)觀測(cè)頻率，且其電流型畸變與頻率無(wú)關(guān)?？上惹蟮眠@N個(gè)頻率的振幅張量和相位張量，然后進(jìn)行上述分解，得到各頻點(diǎn)相應(yīng)的偏離角ψ、相位相對(duì)各向異性參數(shù)ηφ和振幅相對(duì)各向異性參數(shù)ηP，所有頻率同一分解參數(shù)皆可組成1×N階矩陣。為了對(duì)振幅張量和相位張量的分解參數(shù)進(jìn)行對(duì)比，定義目標(biāo)函數(shù)

(23)

通過目標(biāo)函數(shù)值的大小可以衡量相位張量與振幅張量的相似程度，從而通過求目標(biāo)函數(shù)的最小值尋找與相位張量相似程度最高的振幅張量，進(jìn)而得到未畸變阻抗張量，實(shí)現(xiàn)畸變校正。

2.2 粒子群優(yōu)化算法

為了求取目標(biāo)函數(shù)(式(23))的最小值，可以使用簡(jiǎn)單有效的粒子群優(yōu)化算法，將式(23)的值作為粒子適應(yīng)度，此問題只有三個(gè)未知數(shù)，使用粒子群算法能夠有效求解。

粒子群算法(Particle Swarm Optimization，PSO)是由Kennedy等[26]提出的一種基于群體智能進(jìn)化的計(jì)算方法。在給定粒子個(gè)數(shù)后，隨機(jī)初始化粒子的位置和速度，然后通過迭代方法優(yōu)化問題。每一次迭代中，粒子通過跟蹤兩個(gè)“極值”更新自身的位置和速度：一個(gè)是粒子本身所找到的最優(yōu)解，即個(gè)體極值Pbest，另一個(gè)是整個(gè)粒子群中所有粒子在歷代搜索過程中所達(dá)到的最優(yōu)解，即全局極值Gbest。每個(gè)粒子都不斷地改變自身在解空間中的速度，當(dāng)找到這兩個(gè)最優(yōu)解后，可進(jìn)一步實(shí)現(xiàn)PSO中最重要的“加速”過程，以盡可能地朝Pbest和Gbest所指向的區(qū)域“飛”去。

假設(shè)對(duì)于M個(gè)粒子，在一個(gè)D維空間進(jìn)行搜索，粒子i的信息可用兩個(gè)D維向量表示，即位置向量xi和速度向量vi

xi=(xi1,xi2,…，xiD)T

(24)

vi=(vi1,vi2,…，viD)T

(25)

通過測(cè)試發(fā)現(xiàn)，使用標(biāo)準(zhǔn)PSO算法容易陷入局部極小，不能很好地求解畸變參數(shù)。為了提高種群多樣性，避免產(chǎn)生早熟收斂現(xiàn)象，本文對(duì)Clerc[27]提出的具有收縮因子的PSO算法進(jìn)行了改進(jìn)。找到Pbest和Gbest后，粒子即可根據(jù)下式更新速度和位置

vid=k[vid+c1r1(Pbest，d-xid)+

c2r2(Gbest，d-xid)]

(26)

xid=xid+vid

i=1,2,…,M;d=1,2,…,D

(27)

c1=c2=2.0+rand(0,0.5)

(28)

通過改進(jìn)學(xué)習(xí)因子，收縮因子k在每一次迭代中也會(huì)相應(yīng)改變，這樣既可以增強(qiáng)粒子的活性，加強(qiáng)算法的全局探索能力，也可以避免因選取的收縮因子過大而錯(cuò)過最優(yōu)值。

數(shù)值實(shí)驗(yàn)結(jié)果表明，當(dāng)粒子群體數(shù)量M為32時(shí)，速度限制為[-3,3]，針對(duì)此場(chǎng)景問題表現(xiàn)出了良好的效果[28]。

2.3 計(jì)算流程

增益因子g的求解非常困難，到目前為止并沒有一個(gè)特別好的求解策略，因此本文將增益因子g并入?yún)^(qū)域阻抗，只考慮扭曲因子t、剪切因子e和各向異性因子s。利用本文方法進(jìn)行畸變校正的流程見圖1，具體求解過程如下。

圖1 本文方法畸變校正流程

(1)給定32個(gè)初始粒子，隨機(jī)初始化每個(gè)粒子的空間位置x(φt,φe,φs)和速度v，空間位置的三維坐標(biāo)與求解的畸變參數(shù)扭曲角φt、剪切角φe和各向異性因子角φs一一對(duì)應(yīng)，需同時(shí)給定畸變參數(shù)的取值范圍；

(2)利用生成粒子的位置信息生成畸變矩陣的逆C-1，利用式(9)校正阻抗張量Zm，然后利用式(5)和式(17)計(jì)算振幅張量P和相位張量Φ；

(3)對(duì)振幅張量P和相位張量Φ進(jìn)行參數(shù)化分解，根據(jù)分解得到的參數(shù)計(jì)算式(23)目標(biāo)函數(shù)的值F(C)，將得到的目標(biāo)函數(shù)值作為粒子的適應(yīng)度；

(4)根據(jù)步驟(3)計(jì)算得到的值不斷更新粒子的位置x(φt,φe,φs)和速度v，直到達(dá)到預(yù)設(shè)條件，迭代達(dá)到最優(yōu)，迭代終止。

本文迭代終止條件設(shè)置為固定的迭代次數(shù)。經(jīng)過數(shù)值實(shí)驗(yàn)，迭代次數(shù)設(shè)置為1000可以滿足此問題的求解。由于求解參數(shù)只有3個(gè)，故無(wú)需考慮求解速度。

3 模型實(shí)驗(yàn)和實(shí)例分析

3.1 理論模型

為了驗(yàn)證方法的可靠性，設(shè)計(jì)了一組三維模型，對(duì)三維阻抗張量進(jìn)行畸變校正。圖2所示為三維電阻率模型及測(cè)點(diǎn)分布圖。在100Ω·m的均勻半空間中分布著三個(gè)電阻率異常體，其電阻率分別為1000、5000、10Ω·m。點(diǎn)距為200m，共25個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)site01～site25從左到右測(cè)點(diǎn)號(hào)依次增大，計(jì)算頻率為0.01～1000Hz，19個(gè)頻點(diǎn)對(duì)數(shù)等間隔分布。使用MT3D軟件進(jìn)行正演計(jì)算[29]。

圖2 三維理論模型及測(cè)點(diǎn)示意圖

對(duì)正演得到的數(shù)據(jù)施加局部電場(chǎng)畸變影響，即人為添加畸變[30-31]。選取8個(gè)點(diǎn)位分布比較典型的測(cè)點(diǎn)(site01、site05、site08、site10、site13、site15、site19、site23)，隨機(jī)生成了各測(cè)點(diǎn)的畸變參數(shù)。使用本文方法對(duì)其中4個(gè)測(cè)點(diǎn)(site10、site13、site15和site19)畸變數(shù)據(jù)進(jìn)行校正，結(jié)果如圖3所示。

由圖3可見，施加局部電場(chǎng)畸變對(duì)相位的影響明顯弱于對(duì)視電阻率的影響，且視電阻率曲線首支出現(xiàn)不同程度的分離；校正后的視電阻率曲線與未施加畸變的視電阻率曲線幾乎重合，可見很好地恢復(fù)了區(qū)域電性結(jié)構(gòu)，證明此方法對(duì)于恢復(fù)除靜位移因子外的其他三個(gè)畸變參數(shù)(t、e和s)是有效的。所選取測(cè)點(diǎn)隨機(jī)生成的畸變參數(shù)及恢復(fù)的參數(shù)如表1所示?？梢姡齭ite19受畸變嚴(yán)重的點(diǎn)恢復(fù)效果有所欠缺外，其他3個(gè)測(cè)點(diǎn)各畸變參數(shù)相對(duì)誤差都小于3%，對(duì)于人為添加畸變數(shù)據(jù)利用本文方法恢復(fù)效果良好。

表1 三維理論模型畸變參數(shù)表

圖3 三維理論模型測(cè)點(diǎn)site10(a)、site13(b)、site15(c)、site19(d)畸變數(shù)據(jù)校正前、后視電阻率(上)和相位(下)對(duì)比

3.2 實(shí)測(cè)數(shù)據(jù)

安徽霍山—舒城盆地位于金寨斷裂南側(cè)，屬山間拗陷盆地(圖4)，呈東西向，基底除見有少量盧鎮(zhèn)關(guān)巖群、佛子嶺巖群變質(zhì)巖外，多為中、晚侏羅世三尖鋪組、鳳凰臺(tái)組山間磨拉石建造，盆地充填物中有毛坦廠組火山巖[32]。對(duì)霍山地區(qū)中生代火山巖盆地的1條測(cè)線上的15個(gè)MT測(cè)點(diǎn)(圖4中紅框區(qū)域)實(shí)測(cè)數(shù)據(jù)應(yīng)用本文方法進(jìn)行校正。

圖4 安徽霍山地區(qū)MT測(cè)點(diǎn)分布圖

MT數(shù)據(jù)的采集使用加拿大鳳凰公司生產(chǎn)的MTU-5A系統(tǒng)，該系統(tǒng)配備MTC-30磁傳感器，實(shí)際測(cè)量中最低頻率為0.00055Hz，最高頻率為320Hz，共40個(gè)頻點(diǎn)。

由于測(cè)區(qū)存在強(qiáng)烈的電磁干擾，需要對(duì)MT數(shù)據(jù)進(jìn)行精細(xì)的整理和處理。首先對(duì)工區(qū)內(nèi)MT01測(cè)線上的15個(gè)MT測(cè)點(diǎn)都進(jìn)行了時(shí)間域及頻率域數(shù)據(jù)預(yù)處理，對(duì)干擾嚴(yán)重的點(diǎn)進(jìn)行了功率譜挑選、飛點(diǎn)剔除、曲線擬合等處理，然后使用相位張量方法對(duì)這15個(gè)測(cè)點(diǎn)進(jìn)行了維性分析，圖5所示為該測(cè)線相位張量分析圖[33]。

一般地，當(dāng)二維偏離度小于3°時(shí)，表明地下構(gòu)造為二維或準(zhǔn)二維，反之則趨向于三維。分析MT01測(cè)線相位張量圖，發(fā)現(xiàn)該剖面10Hz以下頻段的二維偏離度高于其閾值3°，即表現(xiàn)出明顯的三維性，因而傳統(tǒng)的阻抗分解方法不適用于該剖面。在大地電磁勘探中，由于第四系沉積物的存在，一般表現(xiàn)為一維電性結(jié)構(gòu)，理論上其阻抗張量滿足一維地電模型條件下的特征，即TE和TM模式視電阻率曲線的首支重合，實(shí)際上這也是判斷資料是否存在畸變的一個(gè)重要標(biāo)準(zhǔn)[19]，故可以通過對(duì)比校正前、后視電阻率曲線的首支是否重合評(píng)價(jià)校正結(jié)果。

以測(cè)點(diǎn)site01、site02、site07、site15為例分析校正效果。圖6為這4個(gè)測(cè)點(diǎn)經(jīng)本文方法校正前、后視電阻率和相位的對(duì)比?？梢姡Ｕ暗臏y(cè)點(diǎn)site01和測(cè)點(diǎn)site02的XY和YX模式視電阻率偏離較?。恍Ｕ笫字Ы踔睾?，測(cè)點(diǎn)site01視電阻率曲線在10Hz附近分離，測(cè)點(diǎn)site02視電阻率曲線在50Hz附近分離，這與圖5的相位張量分析結(jié)果一致，在相應(yīng)的重合頻段呈現(xiàn)近圓形，分離的頻段相位張量橢圓率增大。測(cè)點(diǎn)site07同樣具有這種特征。測(cè)點(diǎn)site15的原始視電阻率曲線首支分離程度較大，校正后首支重合，這在一定程度上佐證了本文方法的有效性。

圖5 安徽霍山地區(qū)MT01測(cè)線相位張量分析結(jié)果

圖6 測(cè)點(diǎn)site01(a)、site02(b)、site07(c)、site15(d)經(jīng)本文方法校正前、后視電阻率(上)、相位(下)曲線對(duì)比

4 結(jié)論與討論

傳統(tǒng)阻抗分解算法一般是對(duì)區(qū)域構(gòu)造進(jìn)行一維或者二維假設(shè)，然后選擇相應(yīng)的方法進(jìn)行校正，對(duì)于三維區(qū)域構(gòu)造，傳統(tǒng)的阻抗分解方法失效。在前人的研究基礎(chǔ)上，本文提出基于大地電磁振幅—相位張量分解進(jìn)行大地電磁阻抗畸變校正，無(wú)需進(jìn)行區(qū)域構(gòu)造維性假設(shè)，對(duì)三維構(gòu)造MT數(shù)據(jù)的畸變校正適用。此外，本文方法無(wú)需將各向異性參數(shù)并入?yún)^(qū)域阻抗，通過改進(jìn)粒子群算法實(shí)現(xiàn)了畸變參數(shù)的可靠求取。理論模型和實(shí)測(cè)數(shù)據(jù)計(jì)算結(jié)果分析表明本文方法能夠有效進(jìn)行大地電磁阻抗張量畸變校正。

本文算法存在以下不足：畸變校正只針對(duì)單一測(cè)點(diǎn)，若實(shí)際觀測(cè)數(shù)據(jù)受噪聲影響嚴(yán)重，則校正結(jié)果不夠理想，有待在后續(xù)研究中改進(jìn)，進(jìn)一步討論多站點(diǎn)、其他條件約束下的畸變校正方法。