林昌洪，譚捍東，舒 晴，佟 拓，譚嘉言

1 中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

2 中國地質(zhì)大學(xué)地下信息探測技術(shù)與儀器教育部重點實驗室，北京 100083

3 中國國土資源航空物探遙感中心，北京 100083

1 引 言

可控源音頻大地電磁法（CSAMT）是在音頻大地電磁法（AMT）基礎(chǔ)上發(fā)展起來的一種人工源頻率域測深方法，在勘探石油、天然氣、地?zé)?、金屬礦產(chǎn)、以及水文和環(huán)境工程中發(fā)揮了重要的作用[1-2］.

然而，與大地電磁法相比，可控源音頻大地電磁資料反演技術(shù)的發(fā)展相對緩慢.其主要原因是源的問題，可控源音頻大地電磁采用人工可控信號源，需要求解有源電磁波波動方程，因此處理方法相對復(fù)雜.目前可控源音頻大地電磁實測資料的反演處理主要采用一維反演技術(shù)和對“被視做遠區(qū)的資料”采用二維大地電磁反演技術(shù).要想獲得“遠區(qū)資料”，需要收發(fā)距足夠大使數(shù)據(jù)觀測區(qū)滿足遠區(qū)場條件.而人工源方法，場源與接收點之間的距離越大，所接收到的電磁場信號的信噪比越低，加上測區(qū)的一些外在噪音，很難采集到高質(zhì)量的可控源音頻大地電磁資料.另外，地下介質(zhì)的電阻率通常是未知的.因此，難以確定滿足遠區(qū)條件的收發(fā)距大小，也很難保證所采集到資料是“遠區(qū)資料”，尤其是低頻段資料.對這種數(shù)據(jù)采用大地電磁二維反演技術(shù)進行處理，容易得到錯誤的地質(zhì)解釋.因而，有些文獻提出了對“近區(qū)資料”和“過渡區(qū)資料”進行校正的方法[3-4］，使得“近區(qū)和過渡區(qū)資料”滿足平面波大地電磁數(shù)據(jù)的特征.但進行校正通常需要事先假設(shè)地電結(jié)構(gòu)，如果假設(shè)的地電構(gòu)造和實際情況相差較大，則容易得到錯誤的校正結(jié)果，并且對數(shù)據(jù)進行校正的過程中人為的影響因素較大.還有一些文獻資料采用對視電阻率的重新定義來達到不做近場校正的目的[5-7］，但還未投入廣泛應(yīng)用.

為了解決這個問題，可以開發(fā)對可控源音頻大地電磁數(shù)據(jù)進行直接反演的算法，有了直接反演算法，則在實際工作中不需要采集“遠區(qū)資料”，可以減小收發(fā)距，增加信噪比，獲得高質(zhì)量的可控源音頻大地電磁資料，并得出較可靠的地質(zhì)解釋.Routh等[1］在1999年提出了可控源音頻大地電磁全資料一維水平層狀介質(zhì)的反演.王若等[8］在2007年采用網(wǎng)格參數(shù)法和剝層法實現(xiàn)了一維全資料可控源音頻大地電磁反演.2008年，朱威等[9］用阻尼最小二乘法完成了可控源音頻大地電磁一維全區(qū)反演.李帝銓等[10］在2008年基于遺傳算法實現(xiàn)了可控源音頻大地電磁一維最小構(gòu)造反演.何梅興等[11］在2008年實現(xiàn)了可控源音頻大地電磁一維occam反演.湯井田等[12］在2011年采用了差商法實現(xiàn)可控源音頻大地電磁一維最小構(gòu)造反演.1999年，Lu等[2］采用快速松弛算法實現(xiàn)了可控源音頻大地電磁三維源二維地質(zhì)結(jié)構(gòu)的反演方法.2006年，底青云等[13］也實現(xiàn)了可控源音頻大地電磁三維源二維地質(zhì)結(jié)構(gòu)的快速松弛反演.2010年，雷達[14］采用occam反演方法實現(xiàn)了起伏地形下可控源音頻大地電磁數(shù)據(jù)的二維反演.

然而，實際的地下地質(zhì)情況比較復(fù)雜，地下電阻率結(jié)構(gòu)可能是三維分布的，若對這樣的實測資料進行一維或二維反演解釋，很可能得到不可靠的地電模型[15-18］.為了充分發(fā)揮可控源音頻大地電磁法的作用，提高可控源音頻大地電磁在石油、礦產(chǎn)等領(lǐng)域的應(yīng)用效果，應(yīng)考慮開發(fā)可控源音頻大地電磁數(shù)據(jù)三維反演算法.在對可控源音頻大地電磁理論和共軛梯度算法深入分析的基礎(chǔ)上，正演采用交錯采樣有限差分法[19-20］，反演采用共軛梯度法[21-28］，我們實現(xiàn)了可控源音頻大地電磁資料三維共軛梯度反演算法.文中，首先討論計算可控源音頻大地電磁響應(yīng)的三維正演問題，其次介紹可控源音頻大地電磁三維共軛梯度反演方法理論，包括目標(biāo)函數(shù)、反演流程和“擬正演”問題的計算，最后給出理論模型合成數(shù)據(jù)的三維反演結(jié)果.

2 三維正演

人工源條件下，場源附近電磁場總場的變化梯度大，直接數(shù)值模擬總場困難.我們采用把電磁場的總場分解成一次場（背景場）和二次場的疊加策略，分別計算出一次場和二次場，再合成總場[29］.如果電阻率、電場總場、磁場總場分別記為ρ、E（ρ）、H（ρ），背景電阻率（可為均勻半空間或?qū)訝罱橘|(zhì)）、一次 電 場、一次磁場總場分別記為ρb、E1（ρb）、H1（ρb），那么剩余電阻率 Δρ=ρ-ρb，二次電場E2（Δρ）=E（ρ）-E1（ρb），二 次 磁 場 H2（Δρ）=H（ρ）-H1（ρb）.

對于一次電場E1和一次磁場H1的計算，可從電磁場的麥克斯韋方程出發(fā)，結(jié)合勢函數(shù)和電磁場理論，推導(dǎo)出水平層狀介質(zhì)在有限長度電偶源激發(fā)條件下，全空間節(jié)點處的電場和磁場的表達式，然后采用漢克爾變換方法實現(xiàn)可控源音頻大地電磁全空間電場和磁場值的計算，為三維正演提供可靠的一次場值.

2.1 有限差分法求取二次電磁場

對于二次磁場H2，采用交錯采樣有限差分方法求解.將電阻率和電磁場總場滿足的麥克斯韋方程：

減去背景電阻率和一次電磁場滿足的麥克斯韋方程：

可得剩余電阻率和二次電磁場所滿足的麥克斯韋方程：

其中，ω表示角頻率，μ0表示真空磁導(dǎo)率，Je表示源電流密度.

在笛卡兒右手坐標(biāo)系中，用交錯采樣剖分網(wǎng)格[19］在研究區(qū)域內(nèi)對剩余電阻率和二次場滿足的麥克斯韋積分方程（5）和（6）進行離散化，可獲得關(guān)于地下各網(wǎng)格單元采樣點處二次磁場的正演方程：KH2=s， （7）其中，K為對稱的大型稀疏系數(shù)矩陣；H2為待求解各網(wǎng)格單元采樣點處的二次磁場三分量組成的列向量；s為與一次場及邊界場值有關(guān)的列向量.求解該線性方程組從而獲得二次磁場值列向量H2.

在求解二次磁場時，采取研究區(qū)域的空中頂邊界、地下底邊界和四個側(cè)邊界處的二次場值為零的方法處理邊界條件.

求出二次磁場值列向量H2后，地表某觀測點j處的二次磁場值H2j可以用磁場值列向量H2來表示：

這里H2j表示地表第j個觀測點處模型塊中心點的二次磁場值，它與H2之間通過內(nèi)插向量hgTj轉(zhuǎn)化.對于H2j的三個方向x、y、z分量，方程（8）分別對應(yīng)為：

根據(jù)二次電場與二次磁場之間的差分關(guān)系式

以及插值關(guān)系，在地表某觀測點j處的二次電場值E2j也可以用磁場值列向量H2來表示：這里，E2j表示的是地表第j個觀測點處模型塊中心點的二次電場值，它和H2之間通過內(nèi)插向量egTj轉(zhuǎn)化，bj表示與背景電阻率以及一次場值相關(guān)的量.對于E2j的兩個水平方向x、y分量，方程（10）分別對應(yīng)

2.2 三維正演結(jié)果的檢驗

為了檢驗三維正演算法的正確性，我們設(shè)計了一個二維低阻棱柱體地電模型.如圖1所示，大小為200m×100m，頂面埋深100m，走向為Y方向，電阻率為10Ωm的二維棱柱體埋藏于100Ωm均勻半空間.棱柱體中心在X方向離坐標(biāo)原點的距離為4500m.電流為1A，方向為Y方向，長度為1m的電偶源位于地表（y=100m，x=0m，z=0m）.

用三維正演算法計算出該棱柱體在地表處的二次電磁場響應(yīng)，再用二維有限單元法對該棱柱體進行數(shù)值模擬，得到二次電磁場響應(yīng).圖2和圖3對比了頻率1Hz時，y=-100m測線處兩種數(shù)值模擬計算結(jié)果得到的二次電場E2y和二次磁場H2y響應(yīng).圖2和圖3中，上圖顯示實部，下圖顯示虛部.由圖可見，兩種計算結(jié)果除了在E2y的虛部的峰值附近存在略微的差別外，其余結(jié)果都擬合得非常好.該結(jié)果表明三維正演算法的計算結(jié)果是準(zhǔn)確可靠的.

2.3 三維視電阻率和相位響應(yīng)

由于在反演中使用視電阻率和相位數(shù)據(jù)，在計算完一次電磁場和二次電磁場再合成電磁場總場后，還需要通過地表電磁場總場計算三維地電模型的視電阻率和相位響應(yīng).參考卡尼亞視電阻率的定義[1-2］，可以得到由地表總電場和總磁場計算地表視電阻率和相位的表達式：

其中，Ex和Ey分別表示地表x和y方向總電場，Hx和Hy分別表示地表x和y方向總磁場.

3 三維反演

可控源音頻大地電磁數(shù)據(jù)的三維反演是一項復(fù)雜的工作，由于實測數(shù)據(jù)量大，參加反演的模型參數(shù)多，對計算機硬件資源的要求高，需要花費大量的計算時間.如何快速而又可靠地得出三維反演的結(jié)果是可控源音頻大地電磁三維反演問題向?qū)嵱没较虬l(fā)展的關(guān)鍵.針對這個問題，我們采用共軛梯度反演方法求解可控源音頻大地電磁資料的三維反演問題.

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

可控源音頻大地電磁三維共軛梯度反演算法的目標(biāo)函數(shù)定義為：

圖4 可控源音頻大地電磁三維共軛梯度反演算法流程圖Fig.4 The flowchart of CSAMT 3Dconjugate gradient inversion

其中，Dobs表示觀測視電阻率或相位數(shù)據(jù)；F（m）為求取可控源音頻大地電磁響應(yīng)的正演函數(shù)；V為數(shù)據(jù)方差；λ為正則化參數(shù)；L為簡單的二次差分拉普拉斯算子，m表示模型參數(shù)，m0為先驗?zāi)Ｐ?目標(biāo)函數(shù)的梯度相應(yīng)可表示為：

這里，A表示雅可比矩陣，數(shù)據(jù)擬合差向量：e=Dobs-F（m）.

3.2 反演流程

可控源音頻大地電磁三維共軛梯度反演算法的大致流程（見圖4）為：

（1）設(shè)置迭代次數(shù)i=0，輸入反演的初始模型、反演數(shù)據(jù)和模型剖分參數(shù)等；

（2）一維正演，計算一次電場E1和一次磁場H1；

（3）三維正演，解正演方程KH2=s得到二次磁場H2；

（4）由二次磁場H2計算二次電場E2，根據(jù)一次電磁場E1、H1和二次電磁場E2、H2合成地表總電磁場，再由地表總電磁場計算視電阻率ρ和相位φ；

（5）計算數(shù)據(jù)擬合差eTV-1e，如果擬合差達到設(shè)定的精度，退出循環(huán)，結(jié)束程序，否則繼續(xù)；

（6）通過“擬正演”問題計算ATV-1e；

（7）計算目標(biāo)函數(shù)ψi，目標(biāo)函數(shù)的梯度gi；

（8）通過hi=Cgi，βi=hTi（gi-gi-1）／hTi-1gi-1更新搜索方向pi=-hi+βipi-1；

（9）通過“擬正演”問題計算f=Api；

（10）計算模型的更新步長

（11）更新模型參數(shù)mi=mi-1+αipi；

（12）i=i+1，回步驟（3）.

3.3 雅可比矩陣的計算-“擬正演”問題

從反演流程可以看出，反演中只要計算出ATq和Ap（其中，q=V-1e，p表示搜索方向），就可以得到每次模型迭代的修正量.可控源音頻大地電磁三維共軛梯度反演算法不通過逐個計算雅可比矩陣A的每個元素來求取ATq和Ap，而是通過解“擬正演”問題來直接計算出ATq和Ap的值，下面我們介紹如何通過解“擬正演”問題來直接計算出ATq和Ap的值.

由于背景電阻率ρb是固定值，在反演中我們選取剩余電阻率Δρ為反演模型參數(shù)，正演方程（7）兩端同時對模型參數(shù)Δρ求偏導(dǎo)數(shù)，則有：

將公式（11）代入方程（15）可得：

與三維正演相對應(yīng)，反演過程中的電磁場總場對模型參數(shù)的偏導(dǎo)數(shù)也可分解成一次場和二次場來求對模型參數(shù)的偏導(dǎo)數(shù)：

其中，一次場對剩余電阻率Δρ的偏導(dǎo)數(shù)為零（一次場與剩余電阻率Δρ無關(guān)）：

則總場對模型參數(shù)的偏導(dǎo)數(shù)只剩下二次磁場對模型參數(shù)的偏導(dǎo)數(shù)，因此，公式（16）和公式（17）可以進一步改寫為：

把方程（12）中視電阻率和相位構(gòu)成的復(fù)視電阻率對第k個模型參數(shù)Δρ求偏導(dǎo)數(shù)，可得：

將方程（18）代入方程（19），則第k個模型參數(shù)Δρ的擾動所引起的第j個觀測點處復(fù)視電阻率的改變量?ρsxy／?Δρ和?ρsyx／?Δρ可以寫成如下形式：

其中，

由于使用視電阻率和相位資料作為反演數(shù)據(jù)，那么雅可比矩陣A的元素是反演數(shù)據(jù)對Δρ求偏導(dǎo)數(shù)，也就是視電阻率或相位數(shù)據(jù)對Δρ求偏導(dǎo)數(shù)（?ρsxy／?Δρ和?ρsyx／?Δρ），即公式（20）中的模型參數(shù)Δρ的擾動所引起的觀測點處復(fù)視電阻率的改變量.因此，ATV-1e可表示為：

其中，n=1，2，…，N表示反演數(shù)據(jù)，N為視電阻率和相位數(shù)據(jù)的總個數(shù)：測點數(shù)×頻率數(shù)×2.根據(jù)公式（20），并由K為對稱陣可進一步整理得：

其中，gn、Cn由式（20）確定.例如：如果ρsxyj為第n個數(shù)據(jù)，那么令：

則：

把ν1視為解向量，視為方程右端向量，式（24）是一個和式（7）類似的正演方程，我們稱之為“擬正演”方程.通過解一次“擬正演”方程，可以得到v1.把v1代入式（25），可以直接計算出ATV-1e.

同理，對于Ap，經(jīng)過推導(dǎo)可得：

其中，k=1，2，…，M 表示模型參數(shù).

令：

則有：

把t1視為解向量，視為方程右端向量，式（28）也是一個和式（7）類似的正演方程，我們稱之為“擬正演”方程.通過解一次“擬正演”方程，可以得到t1.把t1代入式（29），可以計算出Ap.這樣通過解一次“擬正演”問題，也可以直接得到Ap.

4 理論模型合成數(shù)據(jù)反演算例

基于上面的理論和公式，我們編程開發(fā)了可控源音頻大地電磁三維共軛梯度反演程序.為了檢驗三維反演算法的有效性，我們設(shè)計了兩個三維地電模型.

4.1 模型一

設(shè)計的低阻棱柱體模型如圖5第一行所示.大小為200m×200m×100m，頂面埋深為100m，電阻率為10Ωm的低阻棱柱體埋藏于電阻率為100Ωm的均勻半空間.

取棱柱體中心在地表處的投影點為坐標(biāo)原點，在x=0km，y=-7km，z=0km的地表處放置長度為100m的X方向水平電偶源.三維網(wǎng)格剖分為46×46×33（含10個空氣層）.用可控源音頻大地電磁三維共軛梯度反演程序的正演代碼部分計算出單棱柱體模型在地表所有剖分網(wǎng)格單元中心點處產(chǎn)生的9個頻率（4000、2000、1000、500、200、100、10Hz、1和0.1Hz）的視電阻率和相位數(shù)據(jù).

初始模型為100Ωm均勻半空間，正則化因子λ=10-4，對地表900個測點處（測區(qū)范圍X：-300～300m，Y：-300～300m）的9個頻率視電阻率ρsxy和相位數(shù)據(jù)φxy中加入1%高斯隨機誤差后用可控源音頻大地電磁三維共軛梯度反演程序在PC機上進行反演.PC機的配置（下同）為：Intel（R）core（TM）i7處理器，主頻2.93GHz，內(nèi)存4.0G.經(jīng)過33次反演迭代，耗時22小時11分鐘，數(shù)據(jù)的擬合方差從初始值12.06收斂到0.99迭代結(jié)束.反演的結(jié)果見圖5的第二行，從圖中可以看出三維反演結(jié)果基本與理論模型一致（圖5的第一行）.

當(dāng)場源位于x=0km，y=-7km，z=0km時，從場源相對測區(qū)的位置及所使用的頻率分析，測區(qū)可近似看作遠區(qū).為了檢驗三維反演程序是否可用于對過渡區(qū)和近區(qū)數(shù)據(jù)進行三維反演，其他參數(shù)保持不變，把場源位置改為x=0km，y=-1km，z=0km和x=0km，y=-0.3km，z=0km 時，三維反演的結(jié)果見圖5的第三行和第四行.從圖中可以看出，除了當(dāng)場源位于x=0km，y=-0.3km，z=0km時，三維反演得到的低阻體在Y方向稍微有些拉長（分析其原因可能與場源有關(guān)，低阻體拉長的位置正好位于場源位置附近的下方），其余結(jié)果都基本與理論模型相一致.

4.2 模型二

設(shè)計的低阻棱柱體和高阻棱柱體組合模型如圖6第一行所示.大小為200m×200m×100m，頂面埋深為100m，電阻率分別為10Ωm和1000Ωm的兩個棱柱體埋藏于電阻率為100Ωm的均勻半空間.

和模型一類似，在x=0km，y=-7km，z=0km的地表處放置長度為100m的X方向水平電偶源.三維網(wǎng)格剖分為66×46×33（含10個空氣層）.初始模型為100Ωm均勻半空間，正則化因子λ=10-4，對地表1500個測點處（測區(qū)范圍X：-500～500m，Y：-300～300m）的9個頻率視電阻率ρsxy和相位數(shù)據(jù)φxy中加入1%高斯隨機誤差后用可控源音頻大地電磁三維共軛梯度反演程序在PC機上進行反演.經(jīng)過28次反演迭代，耗時35小時30分鐘，數(shù)據(jù)的擬合方差從初始值10.28收斂到0.98迭代結(jié)束.反演的結(jié)果見圖6的第二行.

圖5 模型一的三維反演結(jié)果圖中第一行表示真實模型，第二行為場源位于x=0km，y=-7km，z=0km處時的三維反演結(jié)果，第三行為場源位于x=0km，y=-1km，z=0km處時的三維反演結(jié)果，第四行為場源位于x=0km，y=-0.3km，z=0km處時的三維反演結(jié)果.黑色虛線表示棱柱體的邊界.第一列為深度150m處的水平截面圖，第二列為X=0m處沿Y方向的垂直斷面圖，第三列為Y=0m處沿X方向的垂直斷面圖.Fig.5 The 3Dinversion results of model 1 The top row shows the test model.The second row shows the results from the 3Dinversion when the dipole source locates at the position（x=0km，y=-7km，z=0km）.The third row shows the results from the 3Dinversion when the dipole source locates at the position（x=0km，y=-1km，z=0km）.The fourth row shows the results from the 3Dinversion when the dipole source locates at the position（x=0km，y=-0.3km，z=0km）.The black dashed lines show the prism margins.The first column shows the horizontal slices at 150m depth，the second column shows the vertical slices at X=0malong the Yaxis，and the third column shows the vertical slices at Y=0m along the Xaxis.

5 結(jié) 論

正演采用交錯采樣有限差分數(shù)值模擬方法，反演采用正則化反演方案和共軛梯度反演思路，將反演中的雅可比矩陣計算問題轉(zhuǎn)為求解兩次“擬正演”問題，得到模型參數(shù)的更新步長，我們實現(xiàn)了可控源音頻大地電磁三維共軛梯度反演算法.該反演算法可用于對有限長度電偶源激發(fā)下采集到的可控源音頻大地電磁全區(qū)（近區(qū)、過渡區(qū)和遠區(qū)）視電阻率和相位資料進行三維反演定量解釋，獲得地下三維模型的電阻率結(jié)構(gòu).理論模型合成數(shù)據(jù)的反演算例驗證了所實現(xiàn)的可控源音頻大地電磁三維共軛梯度反演算法的有效性和穩(wěn)定性.

圖6 模型二的三維反演結(jié)果圖中第一行表示真實模型，第二行為三維反演結(jié)果.黑色虛線表示棱柱體的邊界.第一列為深度150m處的水平截面圖，第二列為X=-200m處沿Y方向的垂直斷面圖，第三列為X=200m處沿Y方向的垂直斷面圖，第四列為Y=0m處沿X方向的垂直斷面圖.Fig.6 The 3Dinversion results of model 2 The top row shows the test model.The second row shows the 3Dinversion results.The black dashed lines show the prism margins.The first column shows the horizontal slices at 150mdepth，the second column shows the vertical slices at X=-200malong the Yaxis，the third column shows the vertical slices at X=200malong the Yaxis，and the forth column shows the vertical slices at Y=0malong the Xaxis.

（References）

[1]Routh P S，Oldenburg D W.Inversion of controlled source audio-frequency magnetotellurics data for a horizontally layered earth.Geophysics，1999，64（6）：1689-1697.

[2]Lu X Y，Unsworth J M，Booker J.Rapid relaxation inversion of CSAMT data：Geophysical Journal International，1999，138（2）：381-392.

[3]Yamashita M，et al.CSAMT case histories with a multichannel CSAMT system and discussion of near-field data correction.Phoenix Geophys，Ltd，1985.

[4]羅延鐘，周玉水，萬樂.一種新的CSAMT資料近場校正方法.勘查地球物理勘查地球化學(xué)文集第20集.北京：地質(zhì)出版社，1996.Luo Y Z，Zhou Y S，Wan L.A New Correction Method for CSAMT Near-field Data.The 20thCorpus of Exploration Geophysics and Geochemistry（in Chinese）.Beijing：Geology Press，1996.

[5]湯井田，何繼善.水平電偶源頻率測深中全區(qū)視電阻率定義的新方法.地球物理學(xué)報，1994，（4）：543-552.Tang J T，He J S.A new method of define the full-zone resistivity in horizontal electric dipole frequency soundings on a layered earth.Chinese J.Geophys.（in Chinese），2011，54（1）：245-256.

[6]蘇發(fā)，何繼善.組合波近區(qū)頻域電磁測深研究.中國科學(xué)（E輯），1996 ，26（3）：240-246.Su F，He J S.Study on the combination wave electromagnetic sounding in frequency domain.Science in China （Series E）（in Chinese），1996，26（3）：240-246.

[7]湯井田，何繼善.水平多層介質(zhì)上水平電偶源頻率電磁測深的阻抗實部等效電阻率.物探與化探，1994，18（2）：92-96.Tang J T，He J S.Impedance real part equivalent resistivity in frequency electromagnetic sounding of horizontal electric double source on horizontal multilayer media.Geophysical &Geochemical Exploration.（in Chinese），1994，18（2）：92-96.

[8]王若，王妙月.一維全資料CSAMT反演.石油地球物理勘探，2007，42（1）：107-114.Wang N，Wang M Y.Inversion of 1-D full CSAMT data.Oil Geophysical Prospecting （in Chinese），2007，42（1）：107-114.

[9]朱威，李桐林，尚通曉等.CSAMT一維全區(qū)反演對比研究.吉林大學(xué)學(xué)報（地球科學(xué)版），2008，38（增刊）：12-14.Zhu W，Li T L，Shang X T.Compared study of 1-D fullregion inversion of CSAMT.Journal of Jilin University（in Chinese），2008，38：12-14.

[10]李帝銓，王光杰，底青云等.基于遺傳算法的CSAMT最小構(gòu)造反演.地球物理學(xué)報，2008，51（4）：1234-1245.Li D Q，Wang G J，Di Q Y，et al.The application of Genetic Algorithm to CSAMT inversion for minimum structure.Chinese J.Geophys.（in Chinese），2008，51（4）：1234-1245.

[11]何梅興，胡祥云，陳玉萍等.奧克姆一維反演的應(yīng)用.工程地球物理學(xué)報，2008，5（4）：439-443.He M X，Hu X Y，Chen Y P，et al.Application of 1D Occam′s inversion to CSAMT.Chinese Journal of Engineering Geophysics.（in Chinese），2008，5（4）：439-443.

[12]湯井田，張林成，肖嘵等.基于頻點CSAMT一維最小構(gòu)造反演.物探化探計算技術(shù)，2011，33（6）：577-581.Tang J T，Zhang L C，Xiao X，et al.One dimension CSAMT minimum structure inversion based on the frequency.Computing Techniques for Geophysical &Geochemical Exploration.（in Chinese），2011，33（6）：577-581.

[13]底青云，Martyn Unsworth，王妙月.2.5維有限元法CSAMT數(shù)值反演.石油地球物理勘探，2006，41（1）：100-106.Di Q Y，Unsworth M，Wang M Y.2.5-D finite-element CSAMT numerical inversion.Oil Geophysical Prospecting（in Chinese），2006，41（1）：100-106.

[14]雷達.起伏地形下CSAMT二維正反演研究與應(yīng)用.地球物理學(xué)報，2010，53（4）：982-993.Lei D.Studies and applications of 2-D CSAMT modeling and inversion with a dipole source and topography.Chinese J.Geophys.（in Chinese），2010，53（4）：982-993.

[15]胡祖志，胡祥云，何展翔.三維大地電磁數(shù)據(jù)的二維反演解釋.石油地球物理勘探，2005，40（3）：353-359.Hu Z Z，Hu X Y，He Z X.Using 2-D inversion for interpretation of 3-D MT data.Oil Geophysical Prospecting（in Chinese），2005，40（3）：353-359.

[16]Ledo J.2-D Versus 3-D Magnetotelluric Data Interpretation.Surveys in Geophysics，2005，26：511-543.

[17]林昌洪，譚捍東，佟拓.利用大地電磁三維反演方法獲得二維剖面附近三維電阻率結(jié)構(gòu)的可行性.地球物理學(xué)報，2011，54（1）：245-256.Lin C H，Tan H D，Tong T.The possibility of obtaining nearby 3Dresistivity structure from magnetotelluric 2D profile data using 3Dinversion.Chinese J.Geophys.（in Chinese），2011，54（1）：245-256.

[18]Lin C H，Tan H D，Shu Q，et al.Three-dimensional interpretation of sparse survey lines magnetotelluric data：synthetic examples.Applied Geophysics，2012，9（1）：9-18.

[19]譚捍東，余欽范，Booker John等.大地電磁法三維交錯采樣有限差分數(shù)值模擬.地球物理學(xué)報，2003，46（5）：705-711.Tan H D，Yu Q F，Booker J，et al.Magnetotelluric threedimensional modeling using the staggered-grid finite difference method.Chinese J.Geophys.（in Chinese），2003，46（5）：705-711.

[20]Lin C H，Tan H D，Tong T.Parallel rapid relaxation inversion of 3Dmagnetotelluric data.Applied Geophysics，2009，6（1）：77-83.

[21]Mackie R L，Madden T R.Three-dimensional magnetotelluric inversion using conjugate gradients.Geophys.J.Int.，1993，115：215-229.

[22]Newman G A， Alumbaugh D L.Three-dimensional magnetotelluric inversion using non-linear conjugate gradients.Geophys.J.Int.，2000，140：410-424.

[23]Rodi W， Mackie R L.Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion.Geophysics，2001，66：174-187.

[24]胡祖志，胡祥云，何展翔.大地電磁非線性共軛梯度擬三維反演.地球物理學(xué)報，2006，49（4）：1226-1234.Hu Z Z，Hu X Y，He Z X.Pseudo-three-dimensional magnetotelluric inversion using nonlinear conjugate gradients.Chinese J.Geophys.（in Chinese），2006，49（4）：1226-1234.

[25]Lin C H，Tan H D，Tong T.Three-dimensional conjugate gradient inversion of magnetotelluric sounding data.Applied Geophysics，2008，5（4）：314-321.

[26]Lin C H，Tan H D，Tong T.Three-dimensional conjugate gradient inversion of magnetotelluric impedance tensor data.Journal of Earth Science，2011，22（3）：386-395.

[27]林昌洪，譚捍東，佟拓.傾子資料三維共軛梯度反演研究.地球物理學(xué)報，2011，54（4）：1106-1113.Lin C H，Tan H D，Tong T.Three-dimensional conjugate gradient inversion of tipper data.Chinese J.Geophys.（in Chinese），2011，54（4）：1106-1113.

[28]Lin C H，Tan H D，Tong T.Three-dimensional conjugate gradient inversion of magnetotelluric full information data.Applied Geophysics，2011，8（1）：1-10.

[29]Unsworth J M，Bryan J T，Alan D C.Electromagnetic induction by a finite electric dipole source over a 2-D earth.Geophysics，1993，58（2）：198-214.