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

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

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

2 中國(guó)地質(zhì)大學(xué)地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室，北京 100083

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

1 引 言

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

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

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

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

2 三維正演

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

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

2.1 有限差分法求取二次電磁場(chǎng)

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

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

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

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

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

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

求出二次磁場(chǎng)值列向量H2后，地表某觀測(cè)點(diǎn)j處的二次磁場(chǎng)值H2j可以用磁場(chǎng)值列向量H2來(lái)表示：

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

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

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

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

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

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

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

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

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

3 三維反演

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

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

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

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

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

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

3.2 反演流程

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

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

（2）一維正演，計(jì)算一次電場(chǎng)E1和一次磁場(chǎng)H1；

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

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

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

（6）通過(guò)“擬正演”問(wèn)題計(jì)算ATV-1e；

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

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

（9）通過(guò)“擬正演”問(wèn)題計(jì)算f=Api；

（10）計(jì)算模型的更新步長(zhǎng)

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

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

3.3 雅可比矩陣的計(jì)算-“擬正演”問(wèn)題

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

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

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

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

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

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

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

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

其中，

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

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

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

則：

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

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

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

令：

則有：

把t1視為解向量，視為方程右端向量，式（28）也是一個(gè)和式（7）類(lèi)似的正演方程，我們稱(chēng)之為“擬正演”方程.通過(guò)解一次“擬正演”方程，可以得到t1.把t1代入式（29），可以計(jì)算出Ap.這樣通過(guò)解一次“擬正演”問(wèn)題，也可以直接得到Ap.

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

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

4.1 模型一

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

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

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

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

4.2 模型二

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

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

圖5 模型一的三維反演結(jié)果圖中第一行表示真實(shí)模型，第二行為場(chǎng)源位于x=0km，y=-7km，z=0km處時(shí)的三維反演結(jié)果，第三行為場(chǎng)源位于x=0km，y=-1km，z=0km處時(shí)的三維反演結(jié)果，第四行為場(chǎng)源位于x=0km，y=-0.3km，z=0km處時(shí)的三維反演結(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é) 論

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

圖6 模型二的三維反演結(jié)果圖中第一行表示真實(shí)模型，第二行為三維反演結(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]羅延鐘，周玉水，萬(wàn)樂(lè).一種新的CSAMT資料近場(chǎng)校正方法.勘查地球物理勘查地球化學(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]湯井田，何繼善.水平電偶源頻率測(cè)深中全區(qū)視電阻率定義的新方法.地球物理學(xué)報(bào)，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ū)頻域電磁測(cè)深研究.中國(guó)科學(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ì)上水平電偶源頻率電磁測(cè)深的阻抗實(shí)部等效電阻率.物探與化探，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ū)反演對(duì)比研究.吉林大學(xué)學(xué)報(bào)（地球科學(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]李帝銓?zhuān)豕饨?，底青云?基于遺傳算法的CSAMT最小構(gòu)造反演.地球物理學(xué)報(bào)，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é)報(bào)，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]湯井田，張林成，肖嘵等.基于頻點(diǎn)CSAMT一維最小構(gòu)造反演.物探化探計(jì)算技術(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]雷達(dá).起伏地形下CSAMT二維正反演研究與應(yīng)用.地球物理學(xué)報(bào)，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é)報(bào)，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等.大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬.地球物理學(xué)報(bào)，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é)報(bào)，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é)報(bào)，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.