可控源音頻大地電磁三維共軛梯度反演研究

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

1 中國地質大學（北京）地球物理與信息技術學院，北京 100083

2 中國地質大學地下信息探測技術與儀器教育部重點實驗室，北京 100083

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

1 引 言

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

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

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

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

2 三維正演

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

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

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

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

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

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

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

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

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

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

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

根據二次電場與二次磁場之間的差分關系式

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

2.2 三維正演結果的檢驗

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

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

2.3 三維視電阻率和相位響應

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

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

3 三維反演

可控源音頻大地電磁數據的三維反演是一項復雜的工作，由于實測數據量大，參加反演的模型參數多，對計算機硬件資源的要求高，需要花費大量的計算時間.如何快速而又可靠地得出三維反演的結果是可控源音頻大地電磁三維反演問題向實用化方向發展的關鍵.針對這個問題，我們采用共軛梯度反演方法求解可控源音頻大地電磁資料的三維反演問題.

3.1 目標函數

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

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

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

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

3.2 反演流程

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

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

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

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

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

（5）計算數據擬合差eTV-1e，如果擬合差達到設定的精度，退出循環，結束程序，否則繼續；

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

（7）計算目標函數ψi，目標函數的梯度gi；

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

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

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

（11）更新模型參數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是固定值，在反演中我們選取剩余電阻率Δρ為反演模型參數，正演方程（7）兩端同時對模型參數Δρ求偏導數，則有：

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

與三維正演相對應，反演過程中的電磁場總場對模型參數的偏導數也可分解成一次場和二次場來求對模型參數的偏導數：

其中，一次場對剩余電阻率Δρ的偏導數為零（一次場與剩余電阻率Δρ無關）：

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

把方程（12）中視電阻率和相位構成的復視電阻率對第k個模型參數Δρ求偏導數，可得：

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

其中，

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

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

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

則：

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

同理，對于Ap，經過推導可得：

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

令：

則有：

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

4 理論模型合成數據反演算例

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

4.1 模型一

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

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

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

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

4.2 模型二

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

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

圖5 模型一的三維反演結果圖中第一行表示真實模型，第二行為場源位于x=0km，y=-7km，z=0km處時的三維反演結果，第三行為場源位于x=0km，y=-1km，z=0km處時的三維反演結果，第四行為場源位于x=0km，y=-0.3km，z=0km處時的三維反演結果.黑色虛線表示棱柱體的邊界.第一列為深度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 結 論

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

圖6 模型二的三維反演結果圖中第一行表示真實模型，第二行為三維反演結果.黑色虛線表示棱柱體的邊界.第一列為深度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資料近場校正方法.勘查地球物理勘查地球化學文集第20集.北京：地質出版社，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]湯井田，何繼善.水平電偶源頻率測深中全區視電阻率定義的新方法.地球物理學報，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]蘇發，何繼善.組合波近區頻域電磁測深研究.中國科學（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]湯井田，何繼善.水平多層介質上水平電偶源頻率電磁測深的阻抗實部等效電阻率.物探與化探，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一維全區反演對比研究.吉林大學學報（地球科學版），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最小構造反演.地球物理學報，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]何梅興，胡祥云，陳玉萍等.奧克姆一維反演的應用.工程地球物理學報，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一維最小構造反演.物探化探計算技術，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數值反演.石油地球物理勘探，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二維正反演研究與應用.地球物理學報，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]胡祖志，胡祥云，何展翔.三維大地電磁數據的二維反演解釋.石油地球物理勘探，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]林昌洪，譚捍東，佟拓.利用大地電磁三維反演方法獲得二維剖面附近三維電阻率結構的可行性.地球物理學報，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等.大地電磁法三維交錯采樣有限差分數值模擬.地球物理學報，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]胡祖志，胡祥云，何展翔.大地電磁非線性共軛梯度擬三維反演.地球物理學報，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]林昌洪，譚捍東，佟拓.傾子資料三維共軛梯度反演研究.地球物理學報，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.