考慮電流型畸變的大地電磁三維反演

李鑫， 白登海， 閆永利， 馬曉冰， 孔祥儒

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京　100029　2 中國(guó)科學(xué)院大學(xué)， 北京　100039

?

考慮電流型畸變的大地電磁三維反演

李鑫1,2， 白登海1， 閆永利1， 馬曉冰1， 孔祥儒1

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 北京1000292 中國(guó)科學(xué)院大學(xué)， 北京100039

摘要大地電磁測(cè)深(MT)的觀測(cè)數(shù)據(jù)易受到由近地表小尺度非均勻體或地形起伏引起的電流型畸變干擾，消除或壓制這種干擾對(duì)獲取可靠的深部電性結(jié)構(gòu)至關(guān)重要.當(dāng)區(qū)域結(jié)構(gòu)為二維時(shí)，電流型畸變可采用張量分解等方法予以消除或壓制.當(dāng)區(qū)域結(jié)構(gòu)為三維時(shí)，畸變問(wèn)題更加復(fù)雜和嚴(yán)重，傳統(tǒng)張量分解方法往往效果不佳或無(wú)效，嚴(yán)重地制約了MT三維反演技術(shù)的實(shí)用性.對(duì)此，本文提出一種考慮電流型畸變的MT三維反演算法，將完整的電流型畸變參數(shù)引入到目標(biāo)函數(shù)，并采用非線性共軛梯度法與電阻率參數(shù)同時(shí)反演，從而達(dá)到壓制畸變的目的.該算法有兩個(gè)關(guān)鍵點(diǎn)：一是通過(guò)分析實(shí)測(cè)數(shù)據(jù)所遭受畸變的分布特征，在目標(biāo)函數(shù)中對(duì)其進(jìn)行有效約束；二是在迭代過(guò)程中，通過(guò)自適應(yīng)地調(diào)整雙正則化因子保障算法的穩(wěn)定和效率.理論模型測(cè)試結(jié)果顯示，常規(guī)三維反演算法不能合理解釋數(shù)據(jù)中的畸變成分，而只能通過(guò)引入虛假異常體強(qiáng)制地?cái)M合受畸變數(shù)據(jù)，從而造成電阻率模型嚴(yán)重失真.與之相比，本文算法能夠在反演中自動(dòng)求解各測(cè)點(diǎn)所受到的畸變，獲得更接近真實(shí)的電阻率模型.

關(guān)鍵詞電流型畸變； 靜態(tài)效應(yīng)； 三維反演； 大地電磁測(cè)深

1引言

大地電磁測(cè)深(MT)是一種探測(cè)地球深部電性結(jié)構(gòu)的有效方法，它通過(guò)在地球表面同步觀測(cè)天然電、磁場(chǎng)分量，并在頻率域根據(jù)各分量之間的線性關(guān)系獲得多種響應(yīng)函數(shù)，如阻抗張量(視電阻率/相位)、傾子等，進(jìn)而通過(guò)反演手段尋找能夠滿足觀測(cè)數(shù)據(jù)并符合物理規(guī)律的電阻率模型(Simpson and Bahr, 2005).然而，實(shí)際觀測(cè)的MT數(shù)據(jù)極易受到電流型畸變(Galvanic distortion)的干擾.若不對(duì)該畸變進(jìn)行妥善處理而貿(mào)然開(kāi)展反演，將難以獲得真實(shí)的地電模型，甚至因此導(dǎo)致錯(cuò)誤的解釋(Jiracek, 1990; 晉光文和孔祥儒, 2006).

過(guò)去數(shù)十年間，國(guó)內(nèi)外同行針對(duì)如何消除或壓制電流型畸變已開(kāi)展了大量的理論分析和技術(shù)實(shí)踐，提出了諸多解決方案，其中被廣泛使用的是阻抗張量分解方法(Swift, 1967; Bahr, 1988; Groom and Bailey, 1989,1991; Groom and Bahr, 1992; McNeice and Jones,2001).這類方法通常假設(shè)區(qū)域結(jié)構(gòu)可簡(jiǎn)化為一、二維結(jié)構(gòu)，進(jìn)而根據(jù)阻抗張量/畸變張量的數(shù)學(xué)/物理特征將畸變從阻抗張量中分解出來(lái)，獲得基本不受畸變干擾的正常阻抗張量.然而，在實(shí)際運(yùn)用中，這類方法受到兩方面的嚴(yán)重制約：首先是其僅適用于較簡(jiǎn)單區(qū)域結(jié)構(gòu)情況(一、二維)，其次是其僅能求解部分畸變參數(shù)，而對(duì)于另外一部分關(guān)鍵的畸變參數(shù)(如靜位移因子)則無(wú)法求取.靜態(tài)位移效應(yīng)會(huì)造成觀測(cè)視電阻率曲線發(fā)生整體性偏移，一般可采用曲線平移校正、標(biāo)志層校正、空間濾波、輔助瞬變電磁或大極距直流電測(cè)深等方法予以應(yīng)對(duì)(Jones, 1988; Pellerin and Hohmann, 1990; Sternberg et al., 1988).這些方法在特定環(huán)境下具有一定效果，但是均存在不足.比如，很難找到一個(gè)普遍適用準(zhǔn)則來(lái)進(jìn)行曲線平移，而標(biāo)志層校正只適用于盆地等特定的地質(zhì)構(gòu)造環(huán)境，其他空間濾波、瞬變電磁或大極距直流電測(cè)深等均需要增加額外的觀測(cè).為此，DeGroot-Hedlin (1995)創(chuàng)造性地將靜位移和構(gòu)造走向參數(shù)同時(shí)引入到二維OCCAM反演中，與電阻率模型一并求解以消除其干擾，取得了較好的效果.

近年來(lái)MT三維正反演技術(shù)發(fā)展迅速，逐漸開(kāi)始被用于一些實(shí)際資料解釋中(Zhang et al., 2015;Xiao et al., 2015).在三維區(qū)域結(jié)構(gòu)下，電流型畸變對(duì)MT觀測(cè)數(shù)據(jù)所造成的影響更加復(fù)雜和嚴(yán)重，上述現(xiàn)有的張量分解或靜校正方法都不再適用，因此發(fā)展一套能夠配合MT三維反演的畸變消除技術(shù)很有必要.已有一些學(xué)者對(duì)此開(kāi)展了研究，如Utada和Munekane(2000)利用水平電場(chǎng)與垂直磁場(chǎng)的空間導(dǎo)數(shù)關(guān)系壓制三維區(qū)域結(jié)構(gòu)下的電流型畸變干擾.Garcia和Jones (2002)在GB分解的基礎(chǔ)上(Groom and Bailey, 1989)，假設(shè)相鄰測(cè)點(diǎn)遭受不同的畸變影響但反映相同的區(qū)域電性結(jié)構(gòu)，利用牛頓法求解三維區(qū)域結(jié)構(gòu)中的畸變參數(shù).Bibby等 (2005)利用不受電流型畸變影響的相位張量獲得地下維度信息，進(jìn)而采用低維度(一，二維)數(shù)據(jù)求解畸變參數(shù).Patro等(2013)直接三維反演了相位張量數(shù)據(jù)，有效壓制了畸變干擾，然而由于相位張量數(shù)據(jù)中不包含幅值信息導(dǎo)致初始電阻率模型的選取對(duì)反演結(jié)果影響很大.Sasaki和Meju (2006)參考DeGroot-Hedlin (1995)的思路，試圖將靜位移參數(shù)與三維電阻率結(jié)果同時(shí)反演.然而，在三維區(qū)域結(jié)構(gòu)中，由于阻抗張量各分量的模值和相位均受到畸變干擾，僅僅采用靜位移參數(shù)來(lái)描述電流畸變是不夠的 (Jones,2011).針對(duì)這一缺陷，Avdeeva等(2015)將完整的電流型畸變參數(shù)引入到三維反演中，采用擬牛頓法將其與電阻率參數(shù)同步反演，取得了較好的效果.

本文首先給出了三維區(qū)域結(jié)構(gòu)下電流型畸變作用于區(qū)域阻抗張量的數(shù)學(xué)表達(dá)形式，并依據(jù)大量實(shí)測(cè)數(shù)據(jù)分析了電流畸變各參數(shù)的大致分布特征.根據(jù)上述信息，建立了新的目標(biāo)函數(shù)，進(jìn)一步采用非線性共軛梯度算法同時(shí)反演了三維電阻率結(jié)構(gòu)和電流型畸變參數(shù).三維理論模型測(cè)試表明，該算法可以有效壓制電流型畸變的影響，獲得更接近真實(shí)的電阻率模型.

2電流型畸變

2.1電流型畸變的物理涵義

電流型畸變是由近地表小尺度非均勻體電性界面處的電荷積累所產(chǎn)生的散射場(chǎng)與區(qū)域結(jié)構(gòu)所產(chǎn)生的正常場(chǎng)矢量疊加所造成的(Berdichevsky and Dmitriev,1976).假設(shè)三維區(qū)域背景的電導(dǎo)率分布為σ0(r)，而近地表分布著多個(gè)局部三維散射體，其中第k個(gè)散射體的電導(dǎo)率為σk(r)，它與區(qū)域電導(dǎo)率的差異為δσk(r)=σk(r)-σ0(r).Chave和Smith(1994)推導(dǎo)了第k個(gè)散射體內(nèi)外任意位置r處的電場(chǎng)E(r)的積分方程解為

=E0(r)+eI(r)+eG(r).

(1)

式中r′為該散射體內(nèi)任意點(diǎn)的空間位置，E(r′)為r′處的電場(chǎng)矢量，Vk是第k個(gè)散射體的體積，g(r,r′)為標(biāo)量全空間格林函數(shù)，積分對(duì)近地表所有局部散射體所占據(jù)的區(qū)域進(jìn)行.右側(cè)第一項(xiàng)E0(r)通常被稱為正常場(chǎng)或區(qū)域場(chǎng)，是由我們感興趣的區(qū)域結(jié)構(gòu)所產(chǎn)生的電場(chǎng).右側(cè)剩余的兩項(xiàng)分別為近地表局部非均勻體產(chǎn)生的感應(yīng)型畸變電場(chǎng)eI(r)和電流型畸變電場(chǎng)eG(r).這兩個(gè)畸變場(chǎng)的強(qiáng)度取決于近地表局部異常體的體積、形狀及其與區(qū)域電導(dǎo)率的差異情況.采用波恩近似，可將(1)式積分內(nèi)的近地表散射體內(nèi)部的電場(chǎng)E(r′)替換為區(qū)域電場(chǎng)E0(r)，得到觀測(cè)電場(chǎng)E(r)與區(qū)域電場(chǎng)E0(r)之間的關(guān)系：

=c(r)E0(r).

(2)

上式中G(r)為格林并矢張量函數(shù).對(duì)于常規(guī)大地電磁的觀測(cè)頻段(尤其在低頻段)，感應(yīng)型畸變迅速減小，通常忽略不計(jì).因此，畸變張量c(r)一般可以簡(jiǎn)化為一個(gè)與頻率無(wú)關(guān)的二階實(shí)數(shù)張量.對(duì)于磁場(chǎng)，一般認(rèn)為在低頻段局部電荷積累對(duì)觀測(cè)磁場(chǎng)B(r)的影響一般較小(Chave and Smith, 1994)，不予考慮.將(2)式代入阻抗張量的一般表達(dá)式，可以得到觀測(cè)阻抗張量D(r)與區(qū)域阻抗張量Z(r)二者的關(guān)系式如下：

D(r)=E(r)B(r)-1=c(r)E0(r)B0(r)-1

=c(r)Z(r),

(3)

上式中D(r)和Z(r)均為二階復(fù)數(shù)張量，隨頻率而變化.當(dāng)觀測(cè)數(shù)據(jù)不受電流型畸變影響時(shí)，c(r)退化為單位矩陣I(晉光文和孔祥儒, 2006).

2.2三維區(qū)域結(jié)構(gòu)中的電流型畸變

當(dāng)區(qū)域電性結(jié)構(gòu)是三維時(shí)，區(qū)域阻抗張量Z的對(duì)角、非對(duì)角分量均不為零，將(3)式展開(kāi)有

(4)

很明顯，觀測(cè)阻抗張量D各分量的幅值和相位均受到了畸變張量c的影響.上式中包含畸變張量c中4個(gè)和區(qū)域阻抗張量Z中8個(gè)(實(shí)、虛部各4個(gè))共12個(gè)待求解參數(shù)，在不作假設(shè)的情況下，無(wú)法根據(jù)觀測(cè)阻抗張量D中8個(gè)(實(shí)、虛部各4個(gè))已知量建立的方程組直接求解這個(gè)非線性的欠定問(wèn)題.

2.3電流型畸變參數(shù)的分布特征

2.1節(jié)中已證實(shí)MT數(shù)據(jù)受到的電流型畸變通常只與觀測(cè)位置附近的淺層環(huán)境相關(guān)，因此具有很強(qiáng)的隨機(jī)性.Degroot-Hedlin (1995)在將靜位移因子引入到二維反演中時(shí)，采用了對(duì)數(shù)域內(nèi)所有測(cè)點(diǎn)所受靜態(tài)位移之和為零作為約束，測(cè)試結(jié)果表明當(dāng)測(cè)點(diǎn)間距較小且數(shù)量較多時(shí)，這種約束是大致合理的.Sasaki和Meju (2006)依據(jù)不同地質(zhì)環(huán)境中采集的MT和瞬變電磁(TEM)數(shù)據(jù)分析了靜態(tài)位移參數(shù)的分布特點(diǎn)，認(rèn)為其大體呈高斯分布，在高度風(fēng)化或沉積低阻蓋層地區(qū)表現(xiàn)為負(fù)靜態(tài)位移(視電阻率曲線下移)，而在火山巖分布地區(qū)表現(xiàn)為正靜態(tài)位移(視電阻率曲線上移).Avdeev等(2015)根據(jù)經(jīng)驗(yàn)假設(shè)電流型畸變張量的對(duì)角分量大體以1為中心呈對(duì)稱分布，非對(duì)角分量大體以0為中心呈對(duì)稱分布.需要指出的是，上述假設(shè)均缺乏充分的實(shí)測(cè)數(shù)據(jù)的證實(shí).

基于美國(guó)臺(tái)陣項(xiàng)目USARRAY(www.usarray.org)所采集的587個(gè)測(cè)點(diǎn)的高質(zhì)量MT數(shù)據(jù)，我們采用相位張量分解方法(Caldwell et al., 2004;Bibby et al., 2005)估算了所有測(cè)點(diǎn)遭受的畸變情況，統(tǒng)計(jì)結(jié)果如圖1所示.這些觀測(cè)數(shù)據(jù)覆蓋范圍廣，地質(zhì)構(gòu)造環(huán)境豐富，具有較好的代表性.分析結(jié)果表明，畸變張量的非對(duì)角分量(c12,c21)取值以0為中心呈大致對(duì)稱；對(duì)角分量(c11,c22)取值主要集中在1的附近，但在1的右側(cè)具有較長(zhǎng)的尾部，在1的左側(cè)出現(xiàn)的次數(shù)迅速下降.

圖1　USARRAY實(shí)測(cè)數(shù)據(jù)所受的電流型畸變參數(shù)值分布非對(duì)角分量(c12,c21)取值以0為中心呈近對(duì)稱分布；對(duì)角分量(c11,c22)在1附近集中出現(xiàn), 但在1的右側(cè)具有更長(zhǎng)的尾部.Fig.1　Histograms of the galvanic distortion parameters of the USARRAY MT dataThe value of off-diagonal elements symmetrically distributed around zero, while the diagonal elements are concentrated around one but have a longer tail on the right side.

3考慮電流型畸變的三維反演算法

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

與其他地球物理反演問(wèn)題一樣，能夠擬合一組MT觀測(cè)數(shù)據(jù)的電阻率模型和電流型畸變參數(shù)都不是唯一的，因此需要根據(jù)這兩類參數(shù)各自的特征在反演過(guò)程中進(jìn)行有效的約束.基于常規(guī)MT反演的目標(biāo)函數(shù)，我們重新定義了一個(gè)考慮電流型畸變的目標(biāo)函數(shù)：

(5)

右側(cè)第一項(xiàng)Φd為觀測(cè)數(shù)據(jù)擬合項(xiàng)，其具體形式為

(6)

右側(cè)第二項(xiàng)Φm為電阻率模型約束項(xiàng)，定義為

(7)

右側(cè)第三項(xiàng)Φc為電流型畸變參數(shù)約束項(xiàng)，c與cprior分別為預(yù)測(cè)和先驗(yàn)畸變參數(shù)，在反演開(kāi)始前同樣需要為c設(shè)置一個(gè)初始狀態(tài)，即初始畸變參數(shù)c0；λ2為該項(xiàng)的正則化因子.從統(tǒng)計(jì)學(xué)的觀點(diǎn)來(lái)看，對(duì)于大致滿足高斯分布的數(shù)據(jù)采用L2范數(shù)解比較合理(王家映, 1998).上文2.3節(jié)分析已表明畸變張量的對(duì)角、非對(duì)角分量大致分別以1和0呈近似對(duì)稱高斯分布，因此我們采用了預(yù)測(cè)和先驗(yàn)畸變參數(shù)二者之差的L2范數(shù)這種約束形式，即

(8)

在沒(méi)有精確的先驗(yàn)信息情況下，同樣可根據(jù)畸變張量的分布特征將cprior設(shè)為單位矩陣I.

3.2非線性共軛梯度反演算法

共軛梯度法(CG)最早由Stiefel(1951)和Hestenes(1951)分別獨(dú)立提出，Hestenes和Stiefel(1952)給出了用CG方法求解正定系數(shù)矩陣的線性方程組的詳細(xì)過(guò)程.經(jīng)過(guò)Fletcher和Reeves(1964)發(fā)展后，CG被用于解決非線性最優(yōu)化問(wèn)題.Mackie和Madden(1993)首先在三維MT反演中采用了CG方法.Rodi和Mackie (2001)發(fā)展了非線性共軛梯度法(NLCG)，并將其成功用于MT二維反演中.由于無(wú)需直接計(jì)算和存儲(chǔ)雅克比矩陣且收斂速度較快，NLCG很快得到了廣泛應(yīng)用.Newman和Alumbaugh(2000，注：該工作實(shí)際在Rodi和Mackie (2001)之后)首先把NLCG方法用于三維MT反演中，有效降低了三維反演的計(jì)算成本.

本文算法以目標(biāo)函數(shù)(5)式關(guān)于電阻率模型m和電流型畸變參數(shù)c的負(fù)梯度方向分別作為二組參數(shù)的初始搜索方向，之后在每一次迭代中以當(dāng)前的兩組Polak-Ribiere共軛梯度方向作為搜索前進(jìn)方向(Polak and Ribière, 1969).在每一次迭代中，通過(guò)三次多項(xiàng)式插值確定兩組最優(yōu)搜索步長(zhǎng)，按照兩組搜索方向和步長(zhǎng)更新兩組參數(shù)，并通過(guò)Armijo條件判斷當(dāng)前迭代中目標(biāo)函數(shù)的下降是否充分，若不充分則需要調(diào)整正則化因子，重新計(jì)算目標(biāo)函數(shù)值并將搜索方向置為當(dāng)前的負(fù)梯度方向，重新開(kāi)始迭代直至收斂.算法的詳細(xì)步驟見(jiàn)圖2.

兩組共軛梯度方向的計(jì)算是上述算法的一個(gè)關(guān)鍵點(diǎn).將(6)—(8)式代入(5)式可得到目標(biāo)函數(shù)的完整表達(dá)形式，將其對(duì)電阻率模型m及畸變參數(shù)c分別求偏導(dǎo)，可以得到當(dāng)前的兩組梯度方向：

(9)

式中Jm為受到電流型畸變c影響的預(yù)測(cè)阻抗張量dpred關(guān)于電阻率模型m的敏感度矩陣(雅克比矩陣)，表達(dá)式如下

(10)

式中N=Nsite×Nfreq為數(shù)據(jù)總量(Nsite為測(cè)點(diǎn)數(shù)，Nfreq為頻點(diǎn)數(shù))，M為電阻率模型m的電阻率參數(shù)總量.對(duì)于第i個(gè)測(cè)點(diǎn)的第j個(gè)頻點(diǎn)而言，

(11)

注意第i個(gè)測(cè)點(diǎn)所受電流型畸變ci與頻率無(wú)關(guān)，同等作用于該測(cè)點(diǎn)的所有頻點(diǎn).由于引起電流型畸變的淺層異常體的規(guī)模遠(yuǎn)小于MT信號(hào)波長(zhǎng)，MT數(shù)據(jù)的橫、縱向分辨率不足以分辨這些異常體，因此一般認(rèn)為電阻率模型m與電流型畸變c無(wú)關(guān)，可得

圖2　自適應(yīng)壓制電流型畸變的三維MT反演算法流程圖Fig.2　Diagram of the 3D MT inversion algorithm with adaptive galvanic distortion suppression

(12)

與常規(guī)的NLCG方法相比，本算法計(jì)算Jm僅需要將各測(cè)點(diǎn)的預(yù)測(cè)畸變參數(shù)與常規(guī)的敏感度矩陣相乘即可.對(duì)于如何高效計(jì)算和存儲(chǔ)常規(guī)的敏感度矩陣，前人已開(kāi)展了大量有意義的工作(Mackie and Madden, 1993; Farquharson and Oldenburg, 1996; 胡祖志等, 2006; Egbert and Kelbert, 2012).本文采用了Egbert和Kelbert (2012)給出的電磁反演問(wèn)題中敏感度矩陣的一般表達(dá)形式，為了縮小計(jì)算量和存儲(chǔ)的規(guī)模，將敏感度矩陣分解為多個(gè)單元分別予以求解.

(9)式中的Jc是電流型畸變的敏感度矩陣，同樣對(duì)于第i個(gè)測(cè)點(diǎn)的第j個(gè)頻點(diǎn)而言,有

(13)

(14)

3.3自適應(yīng)雙正則化因子策略

正則化反演的一個(gè)難點(diǎn)是如何確定合理的正則化因子.(5)式中的λ1若取值過(guò)大將會(huì)對(duì)模型約束過(guò)于嚴(yán)厲，難以擬合觀測(cè)數(shù)據(jù)，而取值過(guò)小又容易過(guò)度擬合，導(dǎo)致模型的穩(wěn)定性欠佳.本文在引入電流型畸變參數(shù)約束項(xiàng)的同時(shí)引入了一個(gè)新的正則化因子λ2，在反演進(jìn)程中，需要不斷調(diào)整這兩個(gè)正則化因子，保持?jǐn)?shù)據(jù)擬合項(xiàng)(Φd)、電阻率模型約束項(xiàng)(Φm)及畸變參數(shù)約束項(xiàng)(Φc)三者的平衡.

正則化因子通?？刹捎脽o(wú)偏風(fēng)險(xiǎn)預(yù)測(cè)法、廣義交叉檢驗(yàn)法、L-曲線法等方法確定，但這些方法或需要已知數(shù)據(jù)的誤差水平，或需要進(jìn)行大量的實(shí)驗(yàn)性計(jì)算，在三維反演中難以實(shí)現(xiàn).近年來(lái)，一種自適應(yīng)正則化因子的方法在MT反演中得到了廣泛應(yīng)用(Zhdanov, 2002; 陳小斌等, 2005; 吳小平和徐果明, 2000).我們將該方法改進(jìn)后拓展到雙正則化因子情況.正則化因子λ1、λ2的初始值仍需在反演前根據(jù)具體情況預(yù)設(shè)，通?？稍O(shè)置為較大的值以保障在反演前期模型的穩(wěn)定.在反演迭代過(guò)程中，當(dāng)目標(biāo)函數(shù)下降到小于某一閾值時(shí)(通常設(shè)置為一極小值)，將正則化因子乘以某一衰減因子獲得到新的正則化因子，重新計(jì)算目標(biāo)函數(shù)值和梯度方向，并將搜索方向置為當(dāng)前的負(fù)梯度方向，以此放寬對(duì)模型的約束，進(jìn)一步擬合觀測(cè)數(shù)據(jù).衰減因子的取值并沒(méi)有明確的標(biāo)準(zhǔn)，Zhdanov(2002)認(rèn)為衰減因子的經(jīng)驗(yàn)取值區(qū)間為[0.5,0.9].此外，我們?cè)谶@個(gè)過(guò)程中引入了一組更新效率因子δm和δc:

(15)

式中δm和δc分別為電阻率模型和畸變參數(shù)的更新效率因子.當(dāng)算法判斷需要衰減正則化因子進(jìn)一步擬合觀測(cè)數(shù)據(jù)時(shí)，通過(guò)比較當(dāng)前的兩組更新效率因子自適應(yīng)地判斷是衰減λ1還是衰減λ2.如果電阻率模型更新效率優(yōu)于畸變參數(shù)(即δm>δc)，衰減λ2以放寬對(duì)畸變參數(shù)的約束，反之衰減λ1以放寬對(duì)模型參數(shù)的約束.通過(guò)采用上述自適應(yīng)調(diào)整正則化因子策略，算法既能滿足對(duì)觀測(cè)數(shù)據(jù)的高效擬合，又能保障電阻率模型和畸變參數(shù)的穩(wěn)定.

4三維理論模型測(cè)試

4.1模型設(shè)置與響應(yīng)數(shù)據(jù)

圖3所示為一個(gè)理論三維電阻率模型(mtrue).在100 Ωm的背景半空間中分布著三個(gè)電阻率分別為10 Ωm、500 Ωm和5000 Ωm的異常塊體.采用交錯(cuò)網(wǎng)格有限差分算法(Staggered-Grid Finite Difference method)作正演計(jì)算(Mackie and Madden, 1993)，共得到225個(gè)測(cè)點(diǎn)的全阻抗張量數(shù)據(jù)(周期范圍為1～1000 s，周期數(shù)為8).該數(shù)據(jù)不含電流畸變成分，下文中稱為正常數(shù)據(jù).由于理論模型測(cè)試的主要目的是探討新方法在壓制電流型畸變上的有效性，所以該數(shù)據(jù)中未加入任何人為噪聲.

4.2反演結(jié)果對(duì)比和分析

4.2.1電阻率模型比較分析

為了便于比較分析，所有反演方案均采用如下設(shè)置：

圖3　三維理論模型示意圖(a) 平面視圖； (b) 沿測(cè)線P1的剖面視圖； 白色小圓點(diǎn)為225個(gè)MT測(cè)點(diǎn)； 黑色虛線為選取的穿過(guò)模型的測(cè)線P1，白色大圓點(diǎn)S1為選取的一個(gè)代表性測(cè)點(diǎn).Fig.3　Sketch map of the 3-D synthetic model(a) Plan view; (b) Cross-section along the profile P1. Small circles are the location of 225 MT sites; the black dashed line represents a selected profile P1, the bigger circle S1 marks a selected site.

(1) 正常數(shù)據(jù)：由4.1節(jié)所示模型的正演計(jì)算所得，不含電流畸變成分和任何噪音；

(2) 畸變數(shù)據(jù)：將隨機(jī)生成的理論畸變參數(shù)(ctrue)作用于正常數(shù)據(jù)產(chǎn)生的受畸變的數(shù)據(jù)(稱為畸變數(shù)據(jù)).

兩類反演方法：

(1) 不考慮電流型畸變因素的常規(guī)三維反演方法(Kelbert et al., 2014)(稱為常規(guī)方法)；

(2) 考慮電流型畸變因素的三維反演算法(本文方法).

測(cè)試中共設(shè)計(jì)了4種測(cè)試方案，方案1和方案2分別為常規(guī)方法和本文方法對(duì)正常數(shù)據(jù)的反演，目的是檢測(cè)在無(wú)畸變條件下(即均采用正常數(shù)據(jù))本文方法與常規(guī)方法所得結(jié)果的一致性.常規(guī)方法的可靠性已經(jīng)過(guò)了大量的測(cè)試，若方案1和方案2的結(jié)果相同或接近，則表明本文方法在面對(duì)正常數(shù)據(jù)時(shí)是可靠的.方案3和方案4分別為常規(guī)方法和本文方法對(duì)畸變數(shù)據(jù)的反演，目的是表明在面對(duì)畸變數(shù)據(jù)時(shí)的常規(guī)方法的局限性以及本文方法的有效性.四種測(cè)試方案的具體情況見(jiàn)表1.

方案1：常規(guī)方法反演正常數(shù)據(jù)

方案1的目的是檢測(cè)常規(guī)方法在數(shù)據(jù)未受到電流型畸變干擾情況下的反演效果.將正常數(shù)據(jù)和先驗(yàn)電阻率模型(mprior)的響應(yīng)數(shù)據(jù)代入(6)式，可以得到數(shù)據(jù)擬合項(xiàng)(Φd)的初值為6.589.經(jīng)過(guò)17次迭代，隨著預(yù)測(cè)電阻率模型逐漸接近理論模型，Φd的值由6.589減小為1.003(圖4a中的黑色實(shí)線)，達(dá)到理想擬合值1(圖4a中的黑色虛線)，滿足擬合條件，迭代正常終止.電阻率模型約束項(xiàng)(Φm)的值反映了電阻率模型的粗糙度，它的初值為0(初始和先驗(yàn)電阻率模型相同，均為均勻半空間)，在迭代過(guò)程中隨著電阻率模型逐漸偏離先驗(yàn)電阻率模型，Φm增大至0.455(圖4b中黑色實(shí)線)，趨于其理論值2.761(圖4b中黑色虛線).圖5b、6b分別展示了方案1反演所得電阻率模型沿P1測(cè)線的剖面圖及不同深度的切片圖，與理論模型(圖5a、6a)對(duì)比，除了三個(gè)異常體底部由于MT方法的特性無(wú)法得到精確約束外，電阻率大小和尺度都得到了較好的還原，尤其是各塊體之間的橫向電性分界面清晰可辨，表明常規(guī)方法對(duì)正常數(shù)據(jù)的反演是正確有效的.

表1　四種反演方案下目標(biāo)函數(shù)中各項(xiàng)的變化情況

方案2：本文方法反演正常數(shù)據(jù)

采用本文方法同時(shí)反演電阻率模型參數(shù)與電流型畸變參數(shù)時(shí)，所面臨的一種風(fēng)險(xiǎn)是可能會(huì)將一些觀測(cè)數(shù)據(jù)中包含地下電性結(jié)構(gòu)的信息誤判為畸變作用，從而導(dǎo)致真實(shí)的電性結(jié)構(gòu)被掩蓋或者扭曲，這相當(dāng)于方法本身帶入了某種不確定的噪聲.因此，有必要驗(yàn)證當(dāng)數(shù)據(jù)未受畸變干擾時(shí)本文方法是否能夠同樣穩(wěn)定地還原真實(shí)電性結(jié)構(gòu).將λ2的初值設(shè)為600，初始畸變張量(c0)和先驗(yàn)畸變張量(cprior)均設(shè)為單位矩陣I，其余參數(shù)與方案1保持一致.經(jīng)過(guò)30次迭代后，該項(xiàng)由6.589減小為1.036(圖4a中藍(lán)色實(shí)線)，同樣達(dá)到預(yù)期目標(biāo)，反演終止.電阻率模型約束項(xiàng)(Φm)的終值為0.508(圖4b中藍(lán)色實(shí)線)，接近方案1結(jié)果.由于數(shù)據(jù)未受到任何畸變影響，整個(gè)反演過(guò)程中畸變參數(shù)約束項(xiàng)(Φc)保持為一接近于0的極小值(圖4c中的藍(lán)色實(shí)線)，符合實(shí)際.電阻率模型反演結(jié)果如圖5c、6c所示，與方案1反演結(jié)果(圖5b、6b)基本一致，三個(gè)異常體均得以成功還原，說(shuō)明在沒(méi)有畸變時(shí)本文方法與常規(guī)方法是一致的，沒(méi)有帶入附加噪聲，也沒(méi)有丟失有效信息.

圖4　四種反演方案下Φd，Φm和Φc三項(xiàng)隨迭代次數(shù)的變化情況(a) 數(shù)據(jù)擬合項(xiàng)(Φd)； (b) 電阻率模型約束項(xiàng)(Φm)；(c) 電流型畸變參數(shù)約束項(xiàng)(Φc).黑色虛線為各項(xiàng)的理想值.Fig.4　Variation of Φd, Φm and Φc versus iteration number for the four different inversion strategies(a) Data misfit term (Φd); (b) Resistivity model constraint term (Φm); (c) Galvanic distortion constraint term (Φc). The black dashed lines represent the theoretical values of these terms.

圖5　比較四種反演方案所得模型沿測(cè)線P1的電阻率剖面圖(a) 理論模型； (b) 方案1所得反演結(jié)果，較好還原了三個(gè)異常體； (c) 方案2所得反演結(jié)果，與方案1基本一致； (d) 方案3所得反演結(jié)果，嚴(yán)重偏離理論模型(如圖中的虛假低阻異常A)； (e) 方案4所得反演結(jié)果， 與方案3相比更接近理論模型， 畸變得到有效壓制.Fig.5　Comparison of resistivity cross-sections along the profile P1 derived from the four inversion strategies(a) Synthetic model; (b) Inversion result derived from strategy 1 (conventional inversion of the undistorted data); (c) Inversion result derived from strategy 2 (new algorithm inversion of the undistorted data); (d) Inversion result derived from strategy 3 (conventional inversion of the distorted data); (e) Inversion result derived from strategy 4 (new algorithm inversion of the distorted data).

方案3：常規(guī)方法反演畸變數(shù)據(jù)

將畸變數(shù)據(jù)和先驗(yàn)電阻率模型(100 Ωm的均勻半空間)的響應(yīng)數(shù)據(jù)代人(6)式，得到數(shù)據(jù)擬合項(xiàng)(Φd)的初值為10.061.其他參數(shù)與方案1相同.經(jīng)過(guò)84次迭代后陷入局部最優(yōu)，反演被強(qiáng)行終止.數(shù)據(jù)擬合項(xiàng)(Φd)最終減小至6.475(圖4a中紅色實(shí)線)，遠(yuǎn)大于理想值1(圖4a中黑色虛線)，數(shù)據(jù)擬合質(zhì)量差；電阻率模型約束項(xiàng)(Φm)的終值為4.651(圖4b中紅色實(shí)線)，遠(yuǎn)大于理論值2.761(圖4b中黑色虛線)，表明反演過(guò)程中常規(guī)算法對(duì)電阻率模型的約束基本失效.同時(shí)，三個(gè)異常塊體的形態(tài)和電阻率幅值均遭受到嚴(yán)重的畸變，深部背景空間中出現(xiàn)了大范圍的虛假低阻區(qū)域(如圖5d中的A)，并且模型淺層出現(xiàn)了大量虛假的小異常體(圖6d).這說(shuō)明，當(dāng)觀測(cè)數(shù)據(jù)含有電流型畸變時(shí)，常規(guī)方法的反演結(jié)果偏離真實(shí)結(jié)構(gòu).

圖6　比較四種反演方案所得電阻率模型在不同深度的水平切片圖(a)—(e)與圖5相同.Fig.6　Comparison of horizontal resistivity slices at different depths derived from the four inversion strategies(a)—(e) same as in Fig.5.

圖7　比較所有測(cè)點(diǎn)受到的電流畸變理論值與方案4所得預(yù)測(cè)值(a) 電流型畸變參數(shù)理論值(ctrue)； (b) 本文方法預(yù)測(cè)的電流型畸變參數(shù)(cpred).Fig.7　Comparison of theoretical and predicted galvanic distortion values(a) Theoretical galvanic distortion values (ctrue); (b) Galvanic distortion values predicted by our new inversion algorithm (cpred).

方案4：本文方法反演畸變數(shù)據(jù)

反演參數(shù)與方案3完全相同，采用本文方法經(jīng)過(guò)77次迭代后，數(shù)據(jù)擬合項(xiàng)(Φd)由10.061降低為2.210(圖4a中綠色實(shí)線)，雖然未達(dá)到理論目標(biāo)值1，但相較于方案3所得最終值(6.475)，數(shù)據(jù)的擬合程度得到了顯著改善；電阻率模型約束項(xiàng)(Φm)(圖4b中綠色實(shí)線)最終值為0.855，畸變參數(shù)模型約束項(xiàng)(Φc)(圖4c中綠色實(shí)線)的最終值為8.635×10-2，均趨近于理論值，表明電阻率模型和電流型畸變參數(shù)都得到了有效地約束.從反演的電阻率模型來(lái)看，盡管數(shù)據(jù)受到了畸變的干擾，本文方法仍較好地還原了理論模型中的三個(gè)異常體.與常規(guī)方法(方案3)的結(jié)果(圖5d, 圖6d)比較，本文方法的結(jié)果(圖5e，圖6e)更接近理論模型，表明反演質(zhì)量得到了明顯改善.

4.2.2電流型畸變參數(shù)比較分析

圖7比較了所有測(cè)點(diǎn)所遭受的理論電流型畸變參數(shù)(ctrue)與本文方法預(yù)測(cè)的電流型畸變參數(shù)(cpred)，二者整體上具有較好的一致性，說(shuō)明本文方法有效地還原了電流型畸變參數(shù).這一點(diǎn)對(duì)于評(píng)價(jià)本文方法的有效性相當(dāng)重要，在反演過(guò)程中，電流型畸變成分如果未能得到有效的識(shí)別和還原，勢(shì)必將參與到電阻率模型的還原中，從而導(dǎo)致電阻率模型的畸變.

4.3數(shù)據(jù)擬合分析

僅根據(jù)數(shù)據(jù)擬合項(xiàng)(Φd)的取值情況并不能全面、真實(shí)地反映數(shù)據(jù)的擬合程度.圖8、圖9分別給出了所有測(cè)點(diǎn)阻抗張量的非對(duì)角和對(duì)角分量在周期1.389 s處的視電阻率和相位平面圖.對(duì)比正常數(shù)據(jù)(圖8a，圖9a)和畸變數(shù)據(jù)(圖8b，圖9b)，可以看出當(dāng)三維區(qū)域阻抗張量受到畸變時(shí)，對(duì)角和非對(duì)角分量的視電阻率和相位均受到不同程度的影響，特別是對(duì)角分量受到的影響更嚴(yán)重.比較圖8b和圖8a, 非對(duì)角分量的視電阻率受到畸變的干擾比較明顯，而相位受到的影響較小，這說(shuō)明即使在三維結(jié)構(gòu)中，非對(duì)角分量仍主要表現(xiàn)為靜態(tài)效應(yīng)特征(即幅值變化大，而相位不變或變化很小).比較圖9b和圖9a, 對(duì)角分量的模值和相位同時(shí)受到嚴(yán)重的畸變干擾，這是由于對(duì)角分量的模值較小，更容易受到電流型畸變的干擾.采用常規(guī)算法反演畸變數(shù)據(jù)時(shí)(方案3)，數(shù)據(jù)的擬合程度很差(圖8c對(duì)比圖8a;圖9c對(duì)比圖9a).面對(duì)同樣的畸變數(shù)據(jù)，本文方法顯著改善了數(shù)據(jù)的擬合質(zhì)量(對(duì)比圖8d和圖8a).即使對(duì)于受畸變影響嚴(yán)重的對(duì)角分量，本文方法的擬合質(zhì)量也比常規(guī)方法好很多(對(duì)比圖9d和圖9a).

為了分析數(shù)據(jù)各頻點(diǎn)的擬合情況，圖10給出了S1點(diǎn)的全阻抗張量數(shù)據(jù)隨周期的變化曲線.從圖中可以看出，在計(jì)算的頻段范圍內(nèi)，與常規(guī)方法(方案3)(圖10中的紅色圓點(diǎn))比較，本文方法(方案4)計(jì)算的阻抗張量(圖10中的綠色三角)更接近理論值(圖10中的黑色實(shí)線).進(jìn)一步表明在全頻段范圍內(nèi)，本文方法的擬合質(zhì)量明顯優(yōu)于常規(guī)方法.

圖8　比較周期1.389 s處方案3和方案4預(yù)測(cè)模型響應(yīng)的非對(duì)角分量擬合情況(a) 正常數(shù)據(jù)； (b) 畸變數(shù)據(jù)； (c) 方案3預(yù)測(cè)的模型響應(yīng)； (d) 方案4預(yù)測(cè)的模型響應(yīng).Fig.8　Comparison of the off-diagonal data fits for strategy 3 and 4 at period 1.389 s(a) Synthetic data; (b) Distorted data; (c) Predicted responses by strategy 3; (d) Predicted responses by strategy 4.

5結(jié)論與討論

在不考慮外界人為噪聲的情況下，對(duì)大地電磁觀測(cè)信號(hào)影響最嚴(yán)重的地質(zhì)噪音當(dāng)屬電流型畸變(包含靜位移效應(yīng)).數(shù)學(xué)上來(lái)說(shuō)，解決該問(wèn)題的關(guān)鍵在于求解一個(gè)非線性的欠定問(wèn)題.當(dāng)區(qū)域結(jié)構(gòu)可簡(jiǎn)化為一、二維時(shí)，通過(guò)阻抗張量分解等方法可以在反演前將部分畸變成分分解出來(lái)，在一定程度上緩解畸變帶來(lái)的危害.然而，當(dāng)區(qū)域結(jié)構(gòu)具有較強(qiáng)的三維性時(shí)，目前尚無(wú)有效措施應(yīng)對(duì)該畸變.本文試圖發(fā)展一種考慮電流型畸變的三維反演算法，在反演過(guò)程中自適應(yīng)地同步擬合電阻率模型和電流型畸變參數(shù)，從而達(dá)到壓制畸變的影響、提升電阻率模型的可信度的目的.

與常規(guī)反演電阻率模型一樣，在反演電流型畸變參數(shù)時(shí)同樣面臨非唯一性問(wèn)題.如何在反演過(guò)程中對(duì)其進(jìn)行有效約束是解決該問(wèn)題的關(guān)鍵.基于統(tǒng)計(jì)分析大量實(shí)測(cè)MT數(shù)據(jù)所遭受的畸變分布特征，我們認(rèn)為，在無(wú)法通過(guò)其他手段獲得更精確的先驗(yàn)畸變信息的情況下，采用單位矩陣作為先驗(yàn)畸變參數(shù)并且根據(jù)預(yù)測(cè)畸變參數(shù)和先驗(yàn)畸變參數(shù)之差的L2范數(shù)來(lái)衡量畸變模型的復(fù)雜程度(粗糙度)是合理、有效的.

為了保障反演算法的穩(wěn)定和效率，我們將二維反演中常用的自適應(yīng)正則化因子策略改進(jìn)后推廣到雙正則化因子情況中，并且通過(guò)引入新的更新效率因子，依據(jù)電阻率模型和畸變兩組參數(shù)的收斂情況，在兩個(gè)正則化因子之間自主選擇需要衰減的因子.

圖9　比較周期1.389 s處方案3和方案4預(yù)測(cè)模型響應(yīng)的對(duì)角分量擬合情況(a)—(d)與圖8相同.Fig.9　Comparison of the diagonal data fits for strategy 3 and 4 at period 1.389 s.(a)—(d) same as in Fig.8.

圖10　比較測(cè)點(diǎn)S1處模型理論響應(yīng)和方案3、4所得預(yù)測(cè)響應(yīng)曲線黑色實(shí)線為理論模型響應(yīng)，紅色圓點(diǎn)為方案3預(yù)測(cè)模型響應(yīng)，綠色正方形為方案4預(yù)測(cè)模型響應(yīng).Fig.10　Comparison of the theoretical and predicted response curves for strategy 3 and 4 at site S1.Solid black line: synthetic model responses; red circle: predicted responses by strategy 3;green square: predicted responses by strategy 4.

在反演初期，偏重于電阻率模型和畸變兩組參數(shù)的穩(wěn)定性；在反演后期，當(dāng)兩組參數(shù)趨于穩(wěn)健后，適度放寬約束，進(jìn)一步擬合觀測(cè)數(shù)據(jù).需要注意的是兩個(gè)正則化因子的初值仍需要在反演前根據(jù)具體情況設(shè)定并會(huì)對(duì)反演結(jié)果產(chǎn)生一定影響.

理論模型測(cè)試表明，對(duì)于受到電流型畸變的觀測(cè)數(shù)據(jù)，常規(guī)的三維反演方法難以還原真實(shí)的電性結(jié)構(gòu).面對(duì)同等強(qiáng)度的畸變干擾，采用本文方法能夠同時(shí)還原真實(shí)的電性結(jié)構(gòu)和電流型畸變參數(shù).最后需要指出的是，雖然理論模型測(cè)試驗(yàn)證了本文方法在壓制電流型畸變干擾上的有效性，但地球深部的實(shí)際情況要遠(yuǎn)比理論模型復(fù)雜，面對(duì)受到各類噪聲干擾的實(shí)測(cè)數(shù)據(jù)，該算法要取得理想的效果仍需要進(jìn)一步的研究和改進(jìn).

致謝文中所用MT實(shí)測(cè)數(shù)據(jù)來(lái)自USARRAY項(xiàng)目.感謝俄勒岡大學(xué)Gary Egbert教授和Anna Kelbert博士提供MT三維正反演代碼和熱心幫助.匿名審稿專家提出的寶貴意見(jiàn)和建議對(duì)于本文的改進(jìn)和提高非常重要.

ReferencesAvdeeva A, Moorkamp M, Avdeev D, et al. 2015. Three-dimensional inversion of magnetotelluric impedance tensor data and full distortion matrix.GeophysicalJournalInternational, 202(1): 464-481. Bahr K. 1988. Interpretation of the magnetotelluric impedance tensor: regional induction and local telluric distortion.JournalofGeophysics, 62(2): 119-127. Berdichevsky M N, Dmitriev V I. 1976. Distortion of magnetic and electrical fields by near-surface lateral inhomogeneities.ActaGeodaetica,GeophysicaetMontanistica, 11(3-4): 447-483.

Bibby H M, Caldwell T G, Brown C. 2005. Determinable and non-determinable parameters of galvanic distortion in magnetotellurics.GeophysicalJournalInternational, 163(3): 915-930. Caldwell T G, Bibby H M, Brown C. 2004. The magnetotelluric phase tensor.GeophysicalJournalInternational, 158(2): 457-469.

Chave A D, Smith J T. 1994. On electric and magnetic galvanic distortion tensor decompositions.JournalofGeophysicalResearch:SolidEarth, 99(B3): 4669-4682.Chen X B, Zhao G Z, Tang J, et al. 2005. An adaptive regularized inversion algorithm for magnetotelluric data.ChineseJournalofGeophysics(in Chinese), 48(4): 937-946, doi: 10.3321/j.issn:0001-5733.2005.04.029.

DeGroot-Hedlin C. 1995. Inversion for regional 2-D resistivity structure in the presence of galvanic scatterers.GeophysicalJournalInternational, 122(3): 877-888. Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems.GeophysicalJournalInternational, 189(1): 251-267. Farquharson C G, Oldenburg D W. 1996. Approximate sensitivities for the electromagnetic inverse problem.GeophysicalJournalInternational, 126(1): 235-252.

Fletcher R, Reeves C M. 1964. Function minimization by conjugate gradients.TheComputerJournal, 7(2): 149-154.

Garcia X, Jones A G. 2002. Decomposition of three-dimensional magnetotelluric data. ∥Proceedings of the Second International Symposium.Amsterdam:Elsevier, 35: 235-250.

Groom R W, Bailey R C. 1989. Decomposition of magnetotelluric impedance tensors in the presence of local three-dimensional galvanic distortion.JournalofGeophysicalResearch:SolidEarthandPlanets, 94(B2): 1913-1925.

Groom R W, Bailey R C. 1991. Analytic investigations of the effects of near-surface three-dimensional galvanic scatterers on MT tensor decompositions.Geophysics, 56(4): 496-518.

Groom R W, Bahr K. 1992. Corrections for nearsurface effects: Decomposition of the magnetotelluric impedance tensor and scaling corrections for regional resistivities: Atutorial.SurveysinGeophysics, 13(4-5): 341-379.

Hestenes M R.1951. Iterative methods for solving linear equations. NAML Report 52-9.Los Angeles, CA: National Bureau of Standards.

Hestenes M R, Stiefel E. 1952. Methods of conjugate gradients for solving linear systems.JournalofResearchoftheNationalBureauofStandards, 49(6): 409-436.

Hu Z Z, Hu X Y, He Z X. 2006. Pseudo-Three-Dimensional magnetotelluric inversion using nonlinear conjugate gradients.ChineseJournalofGeophysics(in Chinese), 49(4): 1226-1234, doi: 10.3321/j.issn:0001-5733.2006.04.039.

Jin G W, Kong X R.2006. The Distortion and Decomposition of Magnetotelluric Impedance Tensor (in Chinese). Beijing: Seismological Press.

Jiracek G R. 1990. Near-surface and topographic distortions in electromagnetic induction.SurveysinGeophysics, 11(2-3): 163-203.

Jones A G. 1988. Static shift of magnetotelluric data and its removal in a sedimentary basin environment.Geophysics, 53(7): 967-978.

Jones A G. 2011. Three-dimensional galvanic distortion of three-dimensional regional conductivity structures: Comment on “Three-dimensional joint inversion for magnetotelluric resistivity and static shift distributions in complex media” by Yutaka Sasaki and Max A. Meju.JournalofGeophysicalResearch:SolidEarth, 116(B12), doi: 10.1029/2011JB008665.

Kelbert A, Meqbel N, Egbert G D, et al. 2014. ModEM: A modular system for inversion of electromagnetic geophysical data.Computers&Geosciences, 66: 40-53.

Mackie R L, Madden T R. 1993. Three-dimensional magnetotelluric inversion using conjugate gradients.GeophysicalJournalInternational, 115(1): 215-229. McNeice G W, Jones A G. 2001. Multisite, multifrequency tensor decomposition of magnetotelluric data.Geophysics, 66(1): 158-173. Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients.GeophysicalJournalInternational, 140(2): 410-424. Patro P K, Uyeshima M, Siripunvaraporn W. 2013. Three-dimensional inversion of magnetotelluric phase tensor data.GeophysicalJournalInternational, 192(1): 58-66. Pellerin L, Hohmann G W. 1990. Transient electromagnetic inversion: Aremedy for magnetotelluric static shifts.Geophysics, 55(9): 1242-1250. Polak E, Ribière G. 1969. Note sur la convergence de méthodes de directions conjuguées. Revue Fran?aise d′Informatique et de Recherche Opérationnelle, 16: 35-43.

Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion.Geophysics, 66(1): 174-187.

Sasaki Y, Meju M A. 2006. Three-dimensional joint inversion for magnetotelluric resistivity and static shift distributions in complex media.JournalofGeophysicalResearch:SolidEarth, 111(B5): B05101.Simpson F, Bahr K. 2005. Practical Magnetotellurics. Cambridge:Cambridge University Press.

Sternberg B K, Washburne J C, Pellerin L. 1988. Correction for the static shift in magnetotellurics using transient electromagnetic soundings.Geophysics, 53(11): 1459-1468.

Stiefel E.1951. Some special methods of relaxation technique.∥Peoceedings of the symposium on simultaneous linear equations and the determination of eigenvalues.Washington, D.C.:U.S. Government Printing Office.

Swift C M Jr.1967. A magnetotelluric investigation of an electrical conductivity anomaly in the Southwestern United States[Ph. D. thesis].Massachusetts:Massachusetts Institute of Technology.

Utada H, Munekane H. 2000. On galvanic distortion of regional

three-dimensional magnetotelluric impedances.GeophysicalJournalInternational, 140(2): 385-398. Wang J Y. 1998. Inverse Theory in Geophysics. Wuhan: China University of Geosciences press.

Wu X P, Xu G M. 2000. Study on 3-D resistivity inversion using conjugate gradient method.ChineseJournalofGeophysics(in Chinese), 43(3): 420-427, doi: 10.3321/j.issn:0001-5733.2000.03.016.

Xiao Q B, Shao G H, Jing L Z, et al. 2015. Eastern termination of the altyn tagh fault, western China: Constraints from a magnetotelluric survey.JournalofGeophysicalResearch:SolidEarth, 120(5): 2838-2858.

Zhang L T, Unsworth M, Jin S, et al. 2015. Structure of the Central Altyn Tagh Fault revealed by magnetotelluric data: New insights into the structure of the northern margin of the India-Asia collision.EarthandPlanetaryScienceLetters, 415: 67-79.

Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. Amsterdam: Elsevier Science.

附中文參考文獻(xiàn)

陳小斌, 趙國(guó)澤, 湯吉等. 2005. 大地電磁自適應(yīng)正則化反演算法. 地球物理學(xué)報(bào),48(4): 937-946, doi: 10.3321/j.issn:0001-5733.2005.04.029.

胡祖志, 胡祥云, 何展翔. 2006. 大地電磁非線性共軛梯度擬三維反演. 地球物理學(xué)報(bào),49(4): 1226-1234, doi: 10.3321/j.issn:0001-5733.2006.04.039.

晉光文, 孔祥儒. 2006. 大地電磁阻抗張量的畸變與分解. 北京: 地震出版社.

王家映. 1998. 地球物理反演理論. 北京: 中國(guó)地質(zhì)大學(xué)出版社.

吳小平, 徐果明. 2000. 利用共軛梯度法的電阻率三維反演研究. 地球物理學(xué)報(bào), 43(3): 420-427, doi: 10.3321/j.issn:0001-5733.2000.03.016.

(本文編輯胡素芳)

基金項(xiàng)目國(guó)家自然科學(xué)基金(41330212，41274102，41474056)資助.

作者簡(jiǎn)介李鑫，男，1986年生，在讀博士研究生，主要研究方向?yàn)榇蟮仉姶?MT)三維反演技術(shù)和青藏高原深部電性結(jié)構(gòu)和動(dòng)力學(xué). E-mail:lixin342@mail.iggcas.ac.cn

doi:10.6038/cjg20160632 中圖分類號(hào)P319

收稿日期2015-12-17，2016-03-14收修定稿

Three-dimensional inversion of magnetotelluric resistivity model with galvanic distortion

LI Xin1,2, BAI Deng-Hai1, YAN Yong-Li1, MA Xiao-Bing1, KONG Xiang-Ru1

1InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China

AbstractGalvanic distortion, which is caused by shallow small-scale inhomogeneities or topography, has long been recognized as an obstacle in magnetotelluric (MT) sounding studies. Removal of these distortions from the field-observed MT data is essential to obtain a reliable model of the subsurface electrical structure. In two-dimensional (2-D) situations, a number of impedance decomposition methods have been widely used for distortion removal. For data from three-dimensional (3-D) structures, however, the effects of galvanic distortion become far more complex and cannot be solved by these traditional methods. This paper proposes a method that employs the Nonlinear Conjugate Gradient algorithm (NLCG) to solve both the regional resistivity structures and the parameters of full galvanic distortion during iterations. The statistical distribution of galvanic distortion is firstly estimated by using the phase tensor method and suggests that the off-diagonal elements of the distortion tensor are quasi-symmetrically distributed around zero, and the diagonal elements are centralized around one. Using the prior information contained in the field data, we build a new objective function that takes the galvanic effects into consideration and adds an extra constraint on the distortion parameters. Furthermore, to ensure stability and efficiency of the algorithm, an adaptive regularized strategy is used to determine the trade-off factors between the data misfit term and the roughness terms of the resistivity model and galvanic distortion respectively. The synthetic model study suggests that conventional inversion approach cannot reasonably explain the distorted part of the synthetic data and will inevitably bring erroneous structures to the resistivity model. In contrast, our new inversion algorithm can reliably retrieve both the regional resistivity structures and galvanic distortion, irrespective of whether galvanic distortion is present or not.

KeywordsGalvanic distortion； Static shift； 3-D inversion； Magnetotelluric

李鑫， 白登海， 閆永利. 2016. 考慮電流型畸變的大地電磁三維反演. 地球物理學(xué)報(bào)，59(6):2302-2315,doi:10.6038/cjg20160632.

Li X, Bai D H, Yan Y L. 2016. Three-dimensional inversion of magnetotelluric resistivity model with galvanic distortion.ChineseJ.Geophys. (in Chinese),59(6):2302-2315,doi:10.6038/cjg20160632.