基于交叉梯度耦合的大地電磁與地震走時資料三維聯(lián)合反演

彭 淼，譚捍東，姜 枚，錢 輝，譚嘉言

1 中國地質(zhì)科學(xué)院地質(zhì)研究所，大陸構(gòu)造與動力學(xué)國家重點實驗室，北京 100037

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

1 引 言

隨著地球深部探測工作的大力開展和地球物理各種單方法正反演技術(shù)的日趨成熟，利用多種方法對同一地質(zhì)體進行探測并做綜合解釋已經(jīng)得到廣泛應(yīng)用［1－8］.但由于受制于方法本身的單一性和局限性，利用單方法反演結(jié)果進行綜合解釋往往很難獲得較為統(tǒng)一的地質(zhì)—地球物理模型.為了突破這一限制，國內(nèi)外學(xué)者對聯(lián)合反演方法進行了深入研究［9－30］.在大地電磁與地震走時資料的聯(lián)合反演方面，研究主要集中在依靠先驗資料所獲得的參數(shù)換算關(guān)系來解決電阻率與速度模型的耦合問題［13－21］.但事實證明，只有在一定地質(zhì)條件下電阻率與速度之間才存在直接的函數(shù)關(guān)系，例如Meju和Gallardo就曾對英國Quorn地區(qū)某沉積地層進行了巖石物性標(biāo)本的測定，發(fā)現(xiàn)電阻率的對數(shù)與速度的對數(shù)之間呈某種線性關(guān)系［24］.不過這一參數(shù)耦合方法對先驗物性關(guān)系的依賴性較大，具有不確定性，而且一旦這種關(guān)系得不到正確描述，聯(lián)合反演結(jié)果會適得其反，因而其實用價值十分有限.

近年來，采用模型結(jié)構(gòu)耦合的聯(lián)合反演方法在國際上較為流行［25］，比較有代表性的成果是由Gallardo和Meju提出的交叉梯度（cross－gradient）耦合方法［26］.他們進行了大地電磁和折射地震二維聯(lián)合反演算法的研究［27］，通過在聯(lián)合反演的目標(biāo)函數(shù)中引入交叉梯度函數(shù)來約束不同參數(shù).這種方法相比傳統(tǒng)的直接關(guān)系的耦合更具普遍適用性，于是一些學(xué)者沿用這一耦合方法進行了各類地球物理數(shù)據(jù)聯(lián)合反演的研究，也都取得了良好的效果［28－29］.最近，已經(jīng)有研究者分別采用直接關(guān)系耦合和交叉梯度耦合研究了大地電磁、重力與地震走時資料的三維聯(lián)合反演的算法，提出了整個算法的架構(gòu)，并對相容和不相容模型做了模型測試［30］.相比國外，國內(nèi)的已有成果更多地關(guān)注數(shù)據(jù)擬合，較少考慮不同物性參數(shù)間的耦合問題［23］.

在研究地球深部殼幔構(gòu)造方面，大地電磁法和天然地震方法一起被視為兩大支柱方法.大地電磁提供的非常低的電阻率一般是地下連通的流體的反應(yīng)，這通常與寬頻地震所獲得的低速體相對應(yīng)，在綜合解釋中相互驗證，互為補充.本文在大地電磁和地震走時資料的三維聯(lián)合反演研究中首次引入交叉梯度，通過設(shè)計單棱柱體模型和雙棱柱體模型對聯(lián)合反演算法進行模型測試，探討該聯(lián)合反演方法的可行性和適用性.

2 大地電磁和地震走時資料的三維正反演方法

2.1 大地電磁三維正反演

求解大地電磁正演問題所常用的數(shù)值模擬方法包括有限單元法、有限差分法、積分方程法和邊界單元法，本文聯(lián)合反演中大地電磁的正演采用三維交錯采樣有限差分算法［31］.

大地電磁反演方法主要包括OCCAM反演，快速松弛反演（RRI），數(shù)據(jù)空間反演（Data Space）和共軛梯度反演法（CG）等.為了使OCCAM反演理論適用于大地電磁三維反演，Siripunvarapor等人在OCCAM法基礎(chǔ)上進行了改進：通過一系列數(shù)學(xué)變換將M行M 列的模型空間的叉積矩陣轉(zhuǎn)化為N行N列的數(shù)據(jù)空間的叉積矩陣，演變成數(shù)據(jù)空間反演法.由于模型參數(shù)M比數(shù)據(jù)個數(shù)N 大得多，反演中直接求取模型協(xié)方差矩陣Cm，因而大大減小了計算量并減低了計算機內(nèi)存需求，是一種較為實用的大地電磁三維反演方法.考慮到上述優(yōu)點，本文聯(lián)合反演算法中的大地電磁模塊采用數(shù)據(jù)空間反演.讀者可參考Siripunvarapor等人的相關(guān)文獻［32，33］詳細(xì)了解該反演方法，本文不做贅述.

2.2 地震走時資料三維正反演

通過記錄地表觀測到的首波初至，獲取從震源（人工震源和天然地震）到地震臺站之間的走時，可以反演觀測臺站下方的速度結(jié)構(gòu).地震走時數(shù)據(jù)的正演問題通過求解程函方程來實現(xiàn)：

式中，T為走時，x為某點的位置向量，v代表該點地震波速度，這一方程描述了不同時刻波前的幾何形態(tài)的變化.

地震走時資料的正演算法主要包括三類：基于斯奈爾定律的兩點射線追蹤法（如打靶法和彎曲法），基于圖論、費馬原理、惠更斯原理的最短路徑算法以及基于程函方程波前擴散類的有限差分?jǐn)?shù)值求解算法（如快速行進法）.本文采用快速行進三維有限差分算法（Fast Marching Method，簡稱為FMM）作為聯(lián)合反演中地震模塊的正演算法.這一算法由Sethian等人［34］首次提出，算法采用迎風(fēng)差分格式，以窄帶技術(shù)為波前擴展方式，通過窄帶來近似代表波前，并用窄帶的擴展推進來近似模擬波前的傳播過程.此后，De Kool和 Rawlinson［35］將該方法做了改進，能同時計算折射波和反射波的射線路徑.

聯(lián)合反演中地震走時資料的反演則采用子空間（Subspace）反演法［36－37］.這一方法同時構(gòu)建多個方向進行最小化迭代，擴展成在模型空間中的n維子空間.在每次迭代過程中，子空間反演的模型更新量可以寫成n個M 維基向量｛aj｝的加權(quán)和：

式中，A＝［aj］為M行n列的投影矩陣，加權(quán)系數(shù)為μj相應(yīng)基向量aj的長度.

子空間反演的關(guān)鍵問題取決于基向量的構(gòu)造，一般以最速下降向量作為構(gòu)造基礎(chǔ).該方法的優(yōu)勢在于模型更新的計算量小，只需要求解n×n階的線性方程.一般n取值較小，例如Rawlinson等［38］對澳大利亞Tasmania的地震走時反演所采用的子空間維數(shù)n小于20.

3 基于交叉梯度耦合的三維聯(lián)合反演算法

3.1 交叉梯度函數(shù)

對大地電磁和地震走時資料做聯(lián)合反演首先要解決電阻率與速度之間的參數(shù)耦合問題.考慮到巖石物性關(guān)系耦合方法不適用于復(fù)雜地質(zhì)條件，本文采用交叉梯度來對電阻率和速度進行耦合［25－27］，其函數(shù)表達(dá)式可寫成Δ Δ

式中，t代表交叉梯度，ms和mr分別代表地震波速度和電阻率的梯度.

交叉梯度函數(shù)表示為速度和電阻率兩種模型單元參量的梯度的叉積，在聯(lián)合反演中通過直接加入到目標(biāo)函數(shù)來實現(xiàn)大地電磁電阻率模型和地震波速度模型的有效連接.交叉梯度具有以下主要性質(zhì)：

（1）兩種模型的梯度向同一方向變化或者其中一個模型的值為常數(shù)時，交叉梯度函數(shù)為零；

（2）模型參數(shù)改變的大小不影響交叉梯度函數(shù)值的大小；

（3）只有當(dāng)兩個模型參數(shù)朝不同方向變化時交叉梯度函數(shù)的值才不為零.這一耦合方法通過增強不同參數(shù)模型間結(jié)構(gòu)的相似性，達(dá)到結(jié)構(gòu)化的參數(shù)耦合的目的.

對于三維聯(lián)合反演而言，交叉梯度為矢量，可以展開為沿x，y，z三個方向的分量：Δ Δ

其中，

對于各交叉梯度分量，分別采用差分算法進行二維空間的離散化，求取各個分量的長度，再對交叉梯度三維矢量作求模計算.

3.2 聯(lián)合反演目標(biāo)函數(shù)

本文聯(lián)合反演算法包括大地電磁和地震兩大模塊，其核心思想是：每次聯(lián)合反演迭代過程中，電阻率模型的更新不僅考慮數(shù)據(jù)擬合和模型約束，而且還兼顧上次迭代更新的地震波速模型，通過計算它與電阻率模型的交叉梯度來連接兩者，與此同時，迭代中地震波速模型的更新也兼顧上次迭代更新的電阻率模型.因此，聯(lián)合反演算法所采用的目標(biāo)函數(shù)包含以下兩部分：

式中，ms和mr分別為地震波速和電阻率模型向量；ms0和mr0分別為地震波速和電阻率先驗?zāi)Ｐ拖蛄?；Cm為模型協(xié)方差矩陣，Cd為數(shù)據(jù)協(xié)方差矩陣，在（8）和（9）式中具有不同Cm和Cd；ds和dr分別為地震走時和大地電磁觀測數(shù)據(jù)向量，g（·）和f（·）分別為地震和大地電磁的正演算子，α、β、λ和η為權(quán)重因子，本文選取的權(quán)重因子如下：α ＝0.1，β＝1.0～3.0，λ＝1.0，0.5（初始值為1.0，對數(shù)比率上的步長為0.5），η＝1.0～3.0；t（·）為交叉梯度算子，其表達(dá)式見式（3）.

由此可見，當(dāng)未加入交叉梯度項時，大地電磁模塊和地震模塊之間沒有數(shù)據(jù)關(guān)聯(lián)，所獲得的結(jié)果就是大地電磁或地震單獨反演的結(jié)果.與單一反演不同，聯(lián)合反演目標(biāo)函數(shù)中都加入了交叉梯度函數(shù)，每迭代一次所修改生成的電阻率（或速度）模型參數(shù)都會傳遞給地震（或大地電磁）模塊，通過計算交叉梯度來實現(xiàn)兩者的耦合，增強兩者的相似性，以達(dá)到聯(lián)合反演的目的.

3.3 網(wǎng)格匹配方式

為了獲得穩(wěn)定而且足夠精確的正演解，不同方法對正演模擬的網(wǎng)格剖分有著不同要求.考慮大地電磁方法的分辨能力、場的衰減特點以及正演計算效率，大地電磁正演網(wǎng)格一般按照不等間隔剖分，且隨深度的增加以及在研究區(qū)域的外圍網(wǎng)格變得稀疏.地震走時反演的正演網(wǎng)格可采用等間隔剖分，也可采用不等間隔剖分.在聯(lián)合反演中，為了對兩種模型的物性參數(shù)進行耦合，需要計算電阻率和速度模型的交叉梯度值，但網(wǎng)格的差異勢必給參數(shù)耦合帶來困難.為解決這一問題，我們在大地電磁地下網(wǎng)格中定義一塊公共反演網(wǎng)格（圖1），并采用如下網(wǎng)格匹配方法：

（1）公共網(wǎng)格為大地電磁和地震方法探測能力相當(dāng)?shù)姆秶?，在此公共網(wǎng)格范圍內(nèi)地震反演網(wǎng)格和大地電磁網(wǎng)格一樣，交叉梯度值也只在這一公共范圍內(nèi)計算，即電阻率和速度值只在這一公共范圍內(nèi)相互約束.

圖1 聯(lián)合反演三維網(wǎng)格剖分示意圖（a）沿水平方向網(wǎng)格剖分；（b）三維有限差分網(wǎng)格單元；（c）沿垂直方向網(wǎng)格剖分，（a）、（c）圖中粗線框內(nèi)為大地電磁與地震走時反演的公共反演網(wǎng)格，該網(wǎng)格以等間隔剖分.Fig.1 Diagram of 3－D grid for joint inversion（a）Horizontal mesh grid；（b）A 3－D element for forward difference calculation；（c）Vertical mesh grid，in the figure（a）and（c）.Dlue bound represents common inversion grid of magnetotelluric grid and seismic grid，which is divided by equal interval.

（2）對地震方法分別定義正演網(wǎng)格和反演網(wǎng)格，地震反演網(wǎng)格的范圍可大于所定義的公共反演網(wǎng)格，但公共部分的網(wǎng)格和大地電磁是一樣的；公共網(wǎng)格外為大地電磁的不等間隔剖分網(wǎng)格.

（3）本文的地震正演網(wǎng)格采用等間隔剖分方式，且正演網(wǎng)格要比反演網(wǎng)格精細(xì)，反演初始速度值和反演獲得的新模型速度值通過插值在正演網(wǎng)格和反演網(wǎng)格之間傳遞.

3.4 聯(lián)合反演算法流程

大地電磁與地震走時資料的三維聯(lián)合反演的算法流程圖如下（圖2）.

聯(lián)合反演的算法框架可劃分為兩大模塊：地震走時反演和大地電磁反演.兩大程序模塊之間不是彼此獨立的，它們之間通過交叉梯度連接.對于地震走時反演模塊，首先輸入地震初始速度模型，對該模型進行FMM射線追蹤，求取理論走時，并利用梯度類反演方法（本文采用子空間反演）計算速度模型的擾動量，獲得新速度模型.再對新速度模型做正演以計算地震數(shù)據(jù)的擬合差（rmsseis）.另一方面，大地電磁的反演與地震反演同時進行，在模型迭代環(huán)節(jié)中與速度模型交替更新；模型更新后以交叉梯度來耦合電阻率模型和速度模型.大地電磁的正演采用三維交錯采樣有限差分算法，反演采用三維MT梯度類反演法（如數(shù)據(jù)空間反演法）.當(dāng)每次電阻率模型和速度模型都更新完成后，比較地震擬合差與大地電磁擬合差（rmsMT）之和與期望值（ε0）的大小，如果rmsseis＋rmsMT≥ε0，那么繼續(xù)迭代，直至第k次迭代的聯(lián)合反演滿足rmsseis＋rmsMT＜ε0為止.

另需強調(diào)的是，為了使地震速度模型參數(shù)與大地電磁的反演模型參數(shù)相匹配，在反演計算中（包括交叉梯度計算）需對電阻率取對數(shù)，當(dāng)下次迭代進行大地電磁的正演計算時再將對數(shù)電阻率轉(zhuǎn)換回來.另外，該算法結(jié)構(gòu)與傳統(tǒng)的順序反演類的聯(lián)合反演算法有顯著不同.后者首先對某一方法（如地震）進行單獨反演，生成新速度模型，然后通過巖石物性關(guān)系將它轉(zhuǎn)化為另一方法的模型參數(shù)（如電阻率），使其作為該方法反演的初始模型進行反演，生成的模型作為聯(lián)合反演模型.

4 算例分析

基于以上理論和算法，在Linux操作系統(tǒng)下采用Fortran語言并結(jié)合Shell腳本成功開發(fā)了三維聯(lián)合反演程序.為了檢驗算法的有效性和正確性，同時對比聯(lián)合反演結(jié)果與單方法反演的差別，本文進行了理論模型合成數(shù)據(jù)的反演研究，分別設(shè)計了單棱柱體模型和雙棱柱體模型，下面依次作分析討論.

4.1 單棱柱體模型

圖2 聯(lián)合反演流程示意圖Fig.2 The flow chart of the joint inversion algorithm

電阻率為1000Ωm，P波速度大小為7.4km／s，大小為12km×12km×3km的單棱柱體頂面埋深為500m，埋藏于電阻率為100Ωm的地下均勻半空間，圍巖的速度為遞增的一維模型，第一層速度大小為5.87km／s，每層遞增0.13km／s，直到第12層為7.33km／s.P波速度模型的地下三維網(wǎng)格剖分為12×12×12，x，y方向分別以單元網(wǎng)格2km進行等間隔剖分，z方向以每個網(wǎng)格0.5km進行等間隔剖分.電阻率模型的網(wǎng)格剖分在速度模型網(wǎng)格的基礎(chǔ)上向外進行了延展，延展網(wǎng)格以不等間隔剖分：沿x，y方向的延展網(wǎng)格間距在兩側(cè)均依次為5km和7km，沿z方向的延展網(wǎng)格間距向下依次為0.8、1、1.5和2km，總體網(wǎng)格剖分為16×16×16.如圖3a、d所示，地面設(shè)計了16個大地電磁測點（紅色三角），25個震源（紅色五角星）和100個地震觀測臺站（黑色倒三角）.

大地電磁采用全部波阻抗張量元素，正演采用0.01、0.1、1、10Hz和100Hz五個頻率計算出合成數(shù)據(jù)；利用地震FMM三維正演程序計算走時，再分別加入1%的高斯隨機誤差作為觀測數(shù)據(jù).對模型的合成數(shù)據(jù)分別進行了單獨反演和聯(lián)合反演，初始模型采用電阻率為100Ωm的均勻半空間和P波速度遞增的一維模型（與棱柱體圍巖結(jié)構(gòu)相同，前文已做描述）.圖3b和e為單棱柱體模型的單獨反演結(jié)果圖.

圖3 單棱柱體模型聯(lián)合反演與單一反演結(jié)果對比圖（a）電阻率理論模型；（b）大地電磁單一反演結(jié)果；（c）聯(lián)合反演電阻率模型；（d）P波速度理論模型；（e）地震走時資料單一反演結(jié)果；（f）聯(lián)合反演P波速度模型.各子圖中上圖為沿z方向俯視2km切片圖，下圖為沿x方向y＝0km位置剖面圖，圖中黑色虛線所示為棱柱體邊界范圍，（a）、（b）、（c）中藍(lán)色實線為P波速度模型的網(wǎng)格邊界.Fig.3 The comparison of joint inversion and separate inversion by a single－prism model（a）Theoretical resistivity model；（b）Separate MT inversion results；（c）Resistivity model from joint inversion；（d）Theoretical P－wave velocity model；（e）Separate seismic inversion results；（f）P－wave velocity model from joint inversion，in each figure，the one above shows the horizontal slices at 2km depth and the one below shows the vertical slices at y＝0km along the x axis，the black dashed lines show the prism margins and the blue solid lines show seismic grid margins in the figure（a）、（b）and（c）.

結(jié)果顯示：電阻率的反演結(jié)果對單棱柱體的空間形態(tài)有較好的恢復(fù)，但對棱柱體真實電阻率（1000Ωm）卻反演得不夠好，圍巖電阻率出現(xiàn)一定范圍的假異常；P波速度的反演結(jié)果只對單棱柱體淺部有較好恢復(fù)，對棱柱體真實波速（7.4km／s）也恢復(fù)得不好.相比單獨反演，圖3c和f的聯(lián)合反演結(jié)果對單棱柱體的空間形態(tài)有更好的恢復(fù)，棱柱體的反演電阻率更接近于真實電阻率；聯(lián)合反演不僅對單棱柱體淺部速度結(jié)構(gòu)有較好恢復(fù)，而且單棱柱體深部的速度特征也得到了一定程度的恢復(fù).

4.2 雙棱柱體模型

首先設(shè)計了一個含高阻高速和低阻低速異常體的雙棱柱體模型：左側(cè)為高阻高速棱柱體（電阻率為1000Ωm，P波速度大小為7.2km／s），右側(cè)為低阻低速棱柱體（電阻率為10Ωm，P波速度大小為4.8km／s），大小均為16km×8km×4.5km，頂面埋深500m.雙棱柱體埋藏于電阻率為100Ωm的地下均勻半空間，圍巖為速度遞增的一維模型，第一層速度大小為5.36km／s，每層遞增速度為0.14km／s，直到第14層為7.14km／s.P波速度模型的地下三維網(wǎng)格剖分為14×14×14，x，y方向分別以單元網(wǎng)格2km進行等間隔剖分，z方向以每個網(wǎng)格0.5km進行等間隔剖分.電阻率模型的網(wǎng)格向外進行了延展，延展網(wǎng)格以不等間隔剖分：沿x，y方向的延展網(wǎng)格間距在兩側(cè)均依次為3km、5km和7km，沿z方向的延展網(wǎng)格間距向下依次為0.8、1、1.5、2、3km和5km，總體網(wǎng)格剖分為20×20×20.如圖4a和d所示，地面設(shè)計了16個大地電磁測點（黑色三角），25個震源（紅色五角星）和100個地震觀測臺站（黑色倒三角）.

大地電磁正演采用0.01、0.1、1、10Hz和100Hz五個頻率計算出全部阻抗張量的合成數(shù)據(jù)，利用地震FMM三維正演程序計算走時，分別加入1%的高斯隨機誤差作為觀測數(shù)據(jù).對模型的合成數(shù)據(jù)分別進行了單獨反演和聯(lián)合反演，圖4b和e所示的單獨反演結(jié)果顯示：電阻率的反演結(jié)果對雙棱柱體的空間形態(tài)有較好恢復(fù)，但圍巖電阻率出現(xiàn)一定范圍的假異常，特別是在右側(cè)低阻棱柱體深部出現(xiàn)明顯的假低阻異常區(qū)，該假異常應(yīng)是由于三維網(wǎng)格剖分比較稀疏所致（剖分過密會大大增加內(nèi)存需求和計算量）；P波速度的反演結(jié)果則只對淺部有較好恢復(fù)，而且左側(cè)高速棱柱體明顯不如右側(cè)低速棱柱體恢復(fù)得好.聯(lián)合反演在同等網(wǎng)格剖分條件下進行.相比單獨反演，圖4c和f的聯(lián)合反演結(jié)果對異常體的空間形態(tài)和位置都恢復(fù)更好，特別是對速度模型異常體的深部特征有更好恢復(fù)，而且部分消除了電阻率模型中雙棱柱體圍巖的假異常，使反演結(jié)果得到明顯改善.

為了討論交叉梯度與巖石物性關(guān)系耦合方法的區(qū)別，考慮異常體的電阻率和速度不是同向變化的情形（高阻低速或低阻高速），本文進一步設(shè)計了含高阻高速與低阻高速的雙棱柱體模型：由兩個相隔8km的高速的（P波速度大小為7.2km／s）單棱柱體組成，左側(cè)為高阻棱柱體（電阻率為1000Ωm），右側(cè)為低阻棱柱體（電阻率為10Ωm），大小均為12km×6km×3km，頂面埋深500m.圍巖為速度遞增的一維模型，第一層速度大小為5.96km／s，每層遞增速度為0.04km／s，直到第16層為6.54km／s.網(wǎng)格剖分方式、設(shè)計的震源和測點位置如圖5（a）和（d）.

同樣以電阻率為100Ωm的均勻半空間和P波速度遞增的一維模型作為初始模型，對模型的合成數(shù)據(jù)分別進行了單獨反演和聯(lián)合反演.大地電磁的單獨反演使異常體得到一定恢復(fù)，但左側(cè)高阻體的反演值普遍不足200Ωm，與真實電阻率（1000Ωm）有一定差距；對右側(cè)低阻體的反演雖然更接近真實電阻率（10Ωm），但與真實空間位置出現(xiàn)一定偏離.地震波速的反演結(jié)果同樣只對淺部（1.5km以上）有較好恢復(fù)，深度2km切片已顯示所恢復(fù)得到的P波速度值普遍不到6.5km／s，而高速異常體的真實P波速度則為7.2km／s，說明對于這一模型地震走時的單獨反演很難恢復(fù)深度在1.5km以下的高速異常體.而圖5c和f所示的聯(lián)合反演結(jié)果不僅使左側(cè)高阻體和右側(cè)低阻體的反演值都更加接近于真實電阻率，而且對右側(cè)低阻體的空間形態(tài)和位置也有更好的恢復(fù).深度2km切片顯示所恢復(fù)得到的P波速度值可達(dá)到7km／s，很接近真實速度值；而且深度在1.5km至2.5km間的高速異常體得到了一定程度的恢復(fù).

通過對比雙棱柱體模型合成數(shù)據(jù)反演算例，我們認(rèn)為對于電阻率和速度不是同向變化的異常體，通過聯(lián)合反演依然可以同時改善大地電磁和地震走時的單獨反演結(jié)果，這也進一步說明了交叉梯度耦合的特點和優(yōu)勢：并不依賴于不同物性間的特定關(guān)系，而是依靠地下結(jié)構(gòu)的相似性進行反演，因而具有更普遍的適用性.

圖4 高阻高速—低阻低速雙棱柱體模型聯(lián)合反演與單一反演結(jié)果對比圖（a）電阻率理論模型；（b）大地電磁單一反演結(jié)果；（c）聯(lián)合反演電阻率模型；（d）P波速度理論模型；（e）地震走時資料單一反演結(jié)果；（f）聯(lián)合反演P波速度模型.各子圖中由上至下依次為沿z方向俯視3km切片圖，沿x方向y＝0km位置剖面圖和沿y方向x＝2km位置剖面圖，黑色虛線所示為棱柱體邊界范圍，（a）、（b）、（c）中藍(lán)色實線為P波速度模型的網(wǎng)格邊界.Fig.4 The comparison of joint inversion and separate inversion by a resistive－h(huán)igh velocity and conductive－low velocity double－prism model（a）Theoretical resistivity model；（b）Separate MT inversion results；（c）Resistivity model from joint inversion；（d）Theoretical P－wave velocity model；（e）Separate seismic inversion results；（f）P－wave velocity model from joint inversion，each figure from the top down shows the horizontal slices at 3km depth，the vertical slices at y＝0km along the xaxis and the vertical slices at x＝2km along the yaxis，the black dashed lines show the prism margins and the blue solid lines show seismic grid margins in figure（a）、（b）and（c）.

圖5 高阻高速—低阻高速雙棱柱體模型聯(lián)合反演與單一反演結(jié)果對比圖（a）電阻率理論模型；（b）大地電磁單一反演結(jié)果；（c）聯(lián)合反演電阻率模型；（d）P波速度理論模型；（e）地震走時資料單一反演結(jié)果；（f）聯(lián)合反演P波速度模型.各子圖中由上至下依次為沿z方向俯視2km切片圖，沿x方向y＝0km位置剖面圖，沿y方向x＝－6km位置剖面圖，黑色虛線所示為棱柱體邊界范圍，（a）、（b）、（c）中藍(lán)色實線為P波速度模型的網(wǎng)格邊界.Fig.5 The comparison of joint inversion and separate inversion by a resistive－h(huán)igh velocity and conductive－low velocity double－prism model（a）Theoretical resistivity model；（b）Separate MT inversion results；（c）Resistivity model from joint inversion；（d）Theoretical P－wave velocity model；（e）Separate seismic inversion results；（f）P－wave velocity model from joint inversion，each figure from the top down shows the horizontal slices at 2km depth，the vertical slices at y＝0km along the xaxis，the vertical slices at x＝－6km along the yaxis，the black dashed lines show the prism margins and the blue solid lines show seismic grid margins in the figure（a）、（b）and（c）.

5 結(jié) 論

聯(lián)合反演能夠克服單一方法的局限性，減小地球物理反演的多解性，是今后地球物理反演發(fā)展的方向.基于大地電磁數(shù)據(jù)空間三維反演和地震走時資料的子空間三維反演方法，本文實現(xiàn)了在交叉梯度耦合下的三維聯(lián)合反演算法.設(shè)計了單棱柱體模型和雙棱柱體模型對聯(lián)合反演程序進行了理論模型測試.通過對比合成數(shù)據(jù)的聯(lián)合反演和單一反演結(jié)果，分析表明：

（1）無論是單棱柱體模型還是雙棱柱體模型，大地電磁聯(lián)合反演結(jié)果比單獨反演對異常體的空間形態(tài)都有更好的恢復(fù)，其中單棱柱體模型反演的異常體電阻率更接近于真實電阻率，雙棱柱體模型的聯(lián)合反演消除了圍巖的部分假異常.

（2）速度模型的單一反演很難恢復(fù)深部的高速異常體，而聯(lián)合反演不僅對異常體淺部速度結(jié)構(gòu)有更好恢復(fù)，而且使異常體深部的速度特征也能得到一定程度的恢復(fù).

（3）高阻高速－低阻高速的雙棱柱體組合模型的合成數(shù)據(jù)反演算例表明：聯(lián)合反演能同時改善電阻率和速度反向變化異常體的單獨反演結(jié)果，進一步說明叉梯度耦合并不依賴于不同物性間的特定關(guān)系，而是依靠地下結(jié)構(gòu)的相似性，因而具有更普遍的適用性.

（References）

［1］Aric K，Adam A，Smythe D K.Combined seismic and magnetotelluric imaging of upper crystalline crust in the Southern Bohemian Massif.First Break，1997，15（8）：265－271.

［2］Favetto A，Pomposiello C，Booker J，et al.Magnetotelluric inversion constrained by seismic data in the Tucumán Basin（Andean Foothills，27°S，NW Argentina）.J.Geophys.Res.，2007，112（B9），B09104.

［3］Mu?oz，Bauer K，Moeck I，et al.Exploring the Gro? Sch?nebeck（Germany）geothermal site using a statistical joint interpretation of magnetotelluric and seismic tomography models.Geothermics，2010，39（1）：35－45.

［4］高德章，趙金海，薄玉玲等.東海重磁地震綜合探測剖面研究.地球物理學(xué)報，2004，47（5）：853－861.Gao D Z，Zhao J H，Bo Y L，et al.A profile study of gravitative－magnetic and seismic comprehensive survey in the East China Sea.Chinese J.Geophys.（in Chinese），2004，47（5）：853－861.

［5］胡平，劉保金，白立新等.奧林匹克公園地區(qū)隱伏斷裂綜合探測.地球物理學(xué)報，2010，53（6）：1486－1494.Hu P，Liu B J，Bai L X，et al.Synthetic exploration of the buried faults in Olympic park area.Chinese J.Geophys.（in Chinese），2010，53（6）：1486－1494.

［6］李德春，楊書江，胡祖志等.三維重磁電震資料的聯(lián)合解釋——以庫車大北地區(qū)山前礫石層為例.石油地球物理勘探，2012，47（2）：353－359.Li D C，Yang S J，Hu Z Z，et al.Integrated interpretation of 3Dgravity，magnetic，electromagnetic and seismic data：a case study of conglomerate mass investigation in piedmont area of Kuche Depression.Oil Geophysical Prospecting （in Chinese），2012，47（2）：353－359.

［7］姜枚，譚捍東，錢輝等.金川銅鎳礦床的地球物理深部結(jié)構(gòu)與成因模式.礦床地質(zhì)，2012，31（2）：207－215.Jiang M，Tan H D，Qian H，et al.Geophysical deep structure and genetic model of Jinchuan copper－nickel deposit.Mineral Deposits （in Chinese），2012，31（2）：207－215.

［8］姜枚，彭淼，王有學(xué)等.喜馬拉雅東構(gòu)造結(jié)巖石圈板片深俯沖的地球物理證據(jù).巖石學(xué)報，2012，28（6）：1755－1764.Jiang M，Peng M，Wang Y X，et al.Geophysical evidence for deep subduction of Indian lithospheric plate beneath Eastern Himalayan Syntaxis.Acta Petrologica Sinica （in Chinese），2012，28（6）：1755－1764.

［9］Moorkamp M，Jones A G，Eaton D W.Joint inversion of teleseismic receiver functions and magnetotelluric data using a genetic algorithm： Are seismic velocities and electrical conductivities compatible？Geophys.Res.Lett.，2007，34（16）：L16311.

［10］Moorkamp M，Jones A G，F(xiàn)ishwick S.Joint inversion of receiver functions，surface wave dispersion，and magnetotelluric data.J.Geophys.Res.，2010，115（B4）：B04318.

［11］Roux E，Moorkamp M，Jones A G，et al.Joint inversion of long－period magnetotelluric data and surface－wave dispersion curves for anisotropic structure：Application to data from Central Germany.Geophys.Res.Lett.，2011，38（5）：L05304.

［12］彭淼，譚捍東，姜枚等.利用接收函數(shù)和大地電磁數(shù)據(jù)聯(lián)合反演南迦巴瓦構(gòu)造結(jié)中部地區(qū)殼幔結(jié)構(gòu).地球物理學(xué)報，2012，55（7）：2281－2291.Peng M，Tan H D，Jiang M，et al.Joint inversion of receiver functions and magnetotelluric data：Application to crustal and mantle structure beneath central Namche Barwa，eastern Himalayan syntaxis.Chinese J.Geophys.（in Chinese），2012，55（7）：2281－2291.

［13］Kozlovskaya E，Hjelt S E.Modeling of elastic and electrical properties of solid－liquid rock system with fractal microstructure.Physics and Chemistry of the Earth，Part A：Solid Earth and Geodesy，2000，25（2）：195－200.

［14］Heincke B，Jegen M，Hobbs R.Joint inversion of MT，gravity and seismic data applied to sub－basalt imaging.／／SEG Technical Program Expanded Abstracts，2006，25（1）：784－789.

［15］Colombo D，Stefano M D. Geophysical modeling via simultaneous joint inversion of seismic， gravity， and electromagnetic data：application to prestack depth imaging.The Leading Edge，2007，26（3）：326－331.

［16］Vermeesch P M，Morgan J V，Christeson G L，et al.Threedimensional joint inversion of traveltime and gravity data across the Chicxulub impact crater.J.Geophys.Res.，2009，114（6）：B02105.

［17］楊輝，王家林，吳健生等.大地電磁與地震資料仿真退火約束聯(lián)合反演.地球物理學(xué)報，2002，45（5）：723－734.Yang H，Wang J L，Wu J S，et al.Constrained joint inversion of magnetotelluric and seismic data using simulated annealing algorithm.Chinese J.Geophys.（in Chinese），2002，45（5）：723－734.

［18］于鵬，戴明剛，王家林等.電阻率和速度隨機分布的MT與地震聯(lián)合反演.地球物理學(xué)報，2009，52（4）：1089－1097.Yu P，Dai M G， Wang J L，et al.Joint inversion of magnetotelluric and seismic data based on random resistivity and velocity distributions.Chinese J.Geophys.（in Chinese），2009，52（4）：1089－1097.

［19］陳高，于鵬，陳曉等.改進的大地電磁與地震資料聯(lián)合反演方法在黔中隆起區(qū)的適用性研究.石油物探，2010，49（2）：158－163，197.Chen G，Yu P，Chen X，et al.Adaptability research on improved joint inversion method of magnetotelluric and seismic data in middle Guizhou Uplift.Geophys.Pros.for Petroleum.（in Chinese），2010，49（2）：158－163，197.

［20］陳曉，于鵬，張羅磊等.地震與大地電磁測深數(shù)據(jù)的自適應(yīng)正則化同步聯(lián)合反演.地球物理學(xué)報，2011，54（10）：2673－2681.Chen X，Yu P，Zhang L L，et al.Adaptive regularized synchronous joint inversion of MT and seismic data.Chinese J.Geophys.（in Chinese），2011，54（10）：2673－2681.

［21］陳華根，李嘉虓，吳健生等.MT－重力模擬退火聯(lián)合反演研究.地球物理學(xué)報，2012，55（2）：663－670.Chen H G，Li J X，Wu J S，et al.Study on simulatedannealing MT－gravity joint inversion.Chinese J.Geophys.（in Chinese），2012，55（2）：663－670.

［22］何委徽，王家林，于鵬.地球物理聯(lián)合反演研究的現(xiàn)狀與趨勢分析.地球物理學(xué)進展，2009，24（2）：530－540.He W H，Wang J L，Yu P.Overview of the status and prospect of geophysical joint inversion.Progress in Geophysics.（in Chinese），2009，24（2）：530－540.

［23］劉彥，呂慶田，孟貴祥等.大地電磁與地震聯(lián)合反演研究現(xiàn)狀與展望.地球物理學(xué)進展，2012，27（6）：2444－2451.Liu Y，Lu Q T，Meng G X，et al.Joint electromagnetic and seismic inversion survey：status and prospect.Progress in Geophysics （in Chinese），2012，27（6）：2444－2451.

［24］Meju M A，Gallardo L A，Mohamed A K.Evidence for correlation of electrical resistivity and seismic velocity in heterogeneous near－surface materials.Geophys.Res.Lett.，2003，30（7）：1373－1376.

［25］Haber E， Oldenburg D.Joint inversion：a structural approach.Inverse Problems，1997，13（1）：63－77.

［26］Gallardo L A，Meju M A.Characterization of heterogeneous near－surface materials by joint 2Dinversion of DC resistivity and seismic data.Geophys.Res.Lett.，2003，30（13）：1658.

［27］Gallardo L A，Meju M A.Joint two－dimensional crossgradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification.Geophys.J.Int.，2007，169（3）：1261－1272.

［28］Linde N，Binley A， Tryggvason A，et al.Improved hydrogeophysical characterization using joint inversion of cross－h(huán)ole electrical resistance and ground－penetrating radar traveltime data.Water Resources Research，2006，42（12）：W04410.

［29］Hu W Y， Abubakar A， Habasshy T M. Joint electromagnetic and seismic inversion using structural constraints.Geophysics，2009，74（6）：R99－R109.

［30］Moorkamp M，Heincke B，Jegen M，et al.A framework for 3－D joint inversion of MT，gravity and seismic refraction data.Geophys.J.Int.，2011，184（1）：477－493.

［31］譚捍東，余欽范，Booker J等.大地電磁法三維交錯采樣有限差分?jǐn)?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.

［32］Siripunvaraporn W，Egbert G.An efficient data－subspace inversion method for 2－D magnetotelluric data.Geophysics，2000，65（3）：791－803.

［33］Siripunvaraporn W，Egbert G，Lenbury Y，et al.Threedimensional magnetotelluric inversion：data－space method.Phys.Earth Planet.Interiors，2005，150（1－3）：3－14.

［34］Sethian J A，Popovici A M.3－D traveltime computation using the fast marching method.Geophysics，1999，64（2）：516－523.

［35］De Kool M，Rawlinson N，Sambridge M.A practical gridbased method for tracking multiple refraction and reflection phases in three－dimensional heterogeneous media.Geophys.J.Int.，2006，167（1）：253－270.

［36］Kennett B L N，Sambridge M S，Williamson P R.Subspace methods for large inverse problems with multiple parameter classes.Geophys.J.Int.，1988，94（2）：237－247.

［37］Sambridge M S. Non－linear arrival time inversion：constraining velocity anomalies by seeking smooth models in 3－D.Geophys.J.Int.，1990，102（3）：653－677.

［38］Rawlinson N，Reading A M，Kennett B L N.Lithospheric structure of Tasmania from a novel form of teleseismic tomography.J.Geophys.Res.，2006，111（B2）：B02301.