大地電磁與地震反射波走時(shí)數(shù)據(jù)二維聯(lián)合反演

陳　霜，譚捍東，高敬語(yǔ)

(1.核工業(yè)二〇八大隊(duì)，內(nèi)蒙古 包頭　014010；2.中國(guó)地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京　100083)

?

大地電磁與地震反射波走時(shí)數(shù)據(jù)二維聯(lián)合反演

陳霜1,2，譚捍東2，高敬語(yǔ)2

(1.核工業(yè)二〇八大隊(duì)，內(nèi)蒙古 包頭014010；2.中國(guó)地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京100083)

由于地球物理反演存在固有的多解性問(wèn)題，如何減小多解性，提高地球物理方法反演結(jié)果的可靠性和準(zhǔn)確性是當(dāng)今地球物理領(lǐng)域研究的熱點(diǎn)。為了減小大地電磁與地震反射波走時(shí)資料反演的多解性，在Gallardo和Meju研究的交叉梯度理論的基礎(chǔ)上，實(shí)現(xiàn)大地電磁與地震反射波走時(shí)數(shù)據(jù)的二維聯(lián)合反演算法；設(shè)計(jì)理論地電模型和速度模型進(jìn)行合成數(shù)據(jù)的單獨(dú)反演和聯(lián)合反演試算，并驗(yàn)證大地電磁和地震反射波走時(shí)數(shù)據(jù)二維聯(lián)合反演算法的正確性和有效性。反演結(jié)果表明，利用交叉梯度方法耦合大地電磁電阻率參數(shù)和地震波速度參數(shù)能使聯(lián)合反演過(guò)程中的這兩個(gè)模型在更新的時(shí)候相互約束，聯(lián)合反演結(jié)果比單一反演結(jié)果更接近真實(shí)模型，從而使反演結(jié)果更可靠，減少了反演解的多解性。

大地電磁法；地震反射走時(shí)；交叉梯度；二維聯(lián)合反演

0　引　言

地球物理方法可劃分為電法、磁法、地震、重力等分支，根據(jù)地質(zhì)體的物性參數(shù)不同，可以采用不同的地球物理方法對(duì)同一地質(zhì)體問(wèn)題進(jìn)行剖析。但是由于地球物理反演問(wèn)題存在多解性和不適定性，并且單一的地球物理方法獲取的數(shù)據(jù)往往不能彌補(bǔ)其他物性參數(shù)的特點(diǎn)，從而降低了反演結(jié)果的可信度，影響反演的可靠性。為了解決這一問(wèn)題，國(guó)內(nèi)外專家在21世紀(jì)初對(duì)同一地質(zhì)體的不同物性參數(shù)的聯(lián)合反演進(jìn)行了深入的研究[1-13]，其中，頗具代表性的是Gallardo和Meju提出的交叉梯度(cross-gradient)法[1-3]。該方法通過(guò)在反演的目標(biāo)函數(shù)中加入交叉梯度函數(shù)來(lái)增加反演過(guò)程中對(duì)模型的約束，比單一的反演效果更好。近年來(lái)，國(guó)內(nèi)學(xué)者對(duì)聯(lián)合反演研究也有一定的進(jìn)展[14-19]，如陳曉等在聯(lián)合反演中引入Tikhonov正則化思想，實(shí)現(xiàn)了速度參數(shù)和電阻率參數(shù)在隨機(jī)分布情況下的地震資料與大地電磁數(shù)據(jù)的聯(lián)合反演研究[19]，彭淼等實(shí)現(xiàn)了大地電磁和天然地震的三維聯(lián)合反演算法[17]。聯(lián)合反演是反演技術(shù)發(fā)展的必然趨勢(shì)，許多聯(lián)合反演實(shí)例證明了聯(lián)合反演效果的優(yōu)越性[21]。

由于大地電磁探測(cè)和地震勘探各自存在著優(yōu)勢(shì)和局限,因此聯(lián)合反演比單一反演更優(yōu)越[18]。本文通過(guò)對(duì)地質(zhì)體電性參數(shù)和速度參數(shù)的深入了解，開(kāi)展大地電磁與地震反射波走時(shí)數(shù)據(jù)的聯(lián)合反演研究，在其反演目標(biāo)函數(shù)中引入模型光滑項(xiàng)、正則化因子和交叉梯度函數(shù)，通過(guò)設(shè)計(jì)理論地電模型和速度模型對(duì)算法進(jìn)行有效性試驗(yàn)，來(lái)討論該聯(lián)合反演方法的可靠性。

1　大地電磁和地震反射波走時(shí)二維介質(zhì)正演理論

1.1大地電磁法二維正演

有限差分(FDTD)法是一種流行的電磁場(chǎng)數(shù)值計(jì)算方法，已被廣泛地應(yīng)用于求解麥克斯韋方程，具有簡(jiǎn)潔、直觀、計(jì)算速度快等優(yōu)點(diǎn)[26]。大地電磁場(chǎng)滿足麥克斯韋方程組：

▽×E=iωμ0H

(1)

▽×H=σE

(2)

▽×B=0

(3)

▽×D=0

(4)

其中：μ0表示地下介質(zhì)體的磁導(dǎo)率。

在地下介質(zhì)電阻率是二維分布情況下，大地電磁場(chǎng)可分解為T(mén)E和TM兩種模式。根據(jù)大地電磁響應(yīng)數(shù)據(jù)分析，與TE極化模式有關(guān)的場(chǎng)分量滿足以下場(chǎng)方程[25]：

(5)

(6)

(7)

與TM極化模式有關(guān)的場(chǎng)分量滿足以下場(chǎng)方程：

(8)

(9)

(10)

給定特定的邊界條件，對(duì)上述方程組采用有限差分算法進(jìn)行數(shù)值模擬可得到離散節(jié)點(diǎn)的場(chǎng)值，按照下面的公式可進(jìn)一步整理得到視電阻率和相位。

TE極化模式計(jì)算公式：

(11)

(12)

TM極化模式計(jì)算公式：

(13)

(14)

1.2地震反射波走時(shí)二維正演

本文中地震反射波走時(shí)計(jì)算的正演方法采用三點(diǎn)擾動(dòng)的射線追蹤算法[4-7](圖1)，也叫偽彎曲法[23]。該方法運(yùn)用Fermat原理，要求對(duì)一個(gè)初始假設(shè)的射線路徑不斷地進(jìn)行擾動(dòng)，直到其滿足平穩(wěn)狀態(tài)為止。實(shí)際上，也就是將沿著射線路徑段的走時(shí)進(jìn)行最小化計(jì)算，使其達(dá)到穩(wěn)態(tài)效果。即利用射線曲率的性質(zhì)，通過(guò)對(duì)每一個(gè)分段進(jìn)行最小化計(jì)算[22]，使其最終的擾動(dòng)值達(dá)到最小。

如圖1所示，假設(shè)Xk-1和Xk+1是該射線上固定的兩個(gè)端點(diǎn)，已知Xk-1、Xk、Xk+1是相鄰的三個(gè)點(diǎn)，從Xk-1到Xk+1，Xk是其一個(gè)路徑上的初始點(diǎn)，如果有一個(gè)新的點(diǎn)Xk′取代點(diǎn)Xk使走時(shí)最小化，那么Xk的位置將會(huì)被改變。因此，為了尋找這個(gè)新的點(diǎn)Xk′，需要得到兩個(gè)變量，就是從Xk-1、Xk+1中點(diǎn)Xmid出發(fā)的偏移量R以及其方向矢量n。方向矢量n可以通過(guò)獲得反射波最小走時(shí)所走的射線路徑的彎曲方向得出，然后可以確定R的估計(jì)值[20]。

圖1　三點(diǎn)擾動(dòng)方案圖解Fig.1　Chart of three-point perturbation scheme

由于新的Xk′點(diǎn)的射線的切線方向與Xk+1-Xk-1的方向近似平行，因此就能得到近似滿足方程(15)的點(diǎn)Xk′的偏移方向n′，該方向可表示如下：

n′=(gradV)-[(gradV)·(Xk+1-Xk-1)]。

(15)

通過(guò)計(jì)算，可以得到該方向的單位矢量n=n′/|n′ |。

新點(diǎn)Vk′的速度可以近似表達(dá)式如下：Xk′=Xmid+[n·(gradX)mid]R。

因此，沿著方向n的擾動(dòng)值滿足方程：

Rc=-(cVmid+1)/{4cn·(gradV)mid} +{(cVmid+1)2/

[4cn·(gradV)mid]2+L2/(2cVmid)}1/2

(16)

Xk′=Xk+n′ Rc

(17)

2　交叉梯度聯(lián)合反演

2.1交叉梯度函數(shù)

2003年，Gallardo和Meju首次提出交叉梯度函數(shù)的概念[15]，本文以地震波速度參數(shù)和大地電磁電阻率參數(shù)為例，根據(jù)交叉梯度原理，得到其三維計(jì)算公式如下：

t(x,y,z)=▽mr(x,y,z)×▽ms(x,y,z)

(18)

圖2　兩種不同物性模型間的交叉梯度計(jì)算示意圖Fig.2　Charts of cross gradient value for two different property models

式中：▽ms和▽mr分別代表取對(duì)數(shù)的地震波速度和取對(duì)數(shù)的電阻率。交叉梯度函數(shù)表示的是兩種不同物性模型單元參量梯度的叉乘[16]。

交叉梯度具有如下兩條特性，如圖2所示：

(1)當(dāng)兩種物性參數(shù)的其中一種參數(shù)為常數(shù)或二者的變化方向平行時(shí)，交叉梯度值恒等于零；

(2)當(dāng)兩種物性參數(shù)朝不相同的方向變化時(shí)，交叉梯度值不等于零[14]。

設(shè)二維地質(zhì)體的走向沿y軸正方向延伸，則交叉梯度t的z分量和x分量是零，y分量表達(dá)式如下：

(19)

采用五點(diǎn)中心差分網(wǎng)格，將上式用中心差分法離散，得到：

(20)

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

大地電磁反演的方法較多[24]，本文中地震反射波走時(shí)與大地電磁數(shù)據(jù)聯(lián)合反演算法所采用的目標(biāo)函數(shù)為：

(21)

(22)

(23)

(24)

(25)

i=1, 2, ...,N

(26)

在進(jìn)行計(jì)算的時(shí)候，寫(xiě)成矩陣的形式如下：

(27)

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

對(duì)電阻率模型求一階偏導(dǎo)：

(28)

目標(biāo)函數(shù)對(duì)速度模型求一階偏導(dǎo)：

(29)

對(duì)(29)式再求二次偏導(dǎo)：

2λ[BMT]TBMT

(30)

因此，得到大地電磁高斯牛頓反演方法的最終迭代方程，即可求得電阻率迭代更新表達(dá)式：

(31)

地震波速度最小二乘法反演的最終迭代方程可表示為：

(32)

對(duì)應(yīng)的，也可求出地震波速度更新表達(dá)式。

2.3聯(lián)合反演流程

本文實(shí)現(xiàn)大地電磁和地震反射波走時(shí)數(shù)據(jù)聯(lián)合反演計(jì)算的流程如下(圖3)：

(1)首先設(shè)置初始模型。輸入初始速度模型v0和初始電阻率模型m0。

(2)最小二乘法實(shí)現(xiàn)地震波速度模型更新，計(jì)算目標(biāo)函數(shù)梯度gS；求解速度模型擾動(dòng)量Δvi。

(3)高斯牛頓法實(shí)現(xiàn)大地電磁電阻率模型更新[16-17]，計(jì)算目標(biāo)函數(shù)梯度，求一階偏導(dǎo)；計(jì)算近似Hessian矩陣，求二階偏導(dǎo)；計(jì)算電阻率模型擾動(dòng)量Δmi。

(4)計(jì)算速度模型更新vi+1=vi+Δv和電阻率模型更新ρi+1=ρi+Δρ。

(5)計(jì)算地震和大地電磁數(shù)據(jù)擬合差之和。

(6)判斷擬合差之和是否滿足目標(biāo)函數(shù)設(shè)置的收斂條件。如果不滿足，迭代繼續(xù)，返回步驟(2)重新計(jì)算；如果滿足收斂條件，完成聯(lián)合反演程序，跳出迭代，輸出反演結(jié)果(圖3)。

圖4　反演理論模型Fig.4　Model for inversion

3　理論模型合成數(shù)據(jù)算例

在LINUX系統(tǒng)下，利用FORTRAN語(yǔ)言實(shí)現(xiàn)大地電磁與地震反射波走時(shí)數(shù)據(jù)的二維聯(lián)合反演，設(shè)計(jì)了三個(gè)不同的二維模型算例，來(lái)驗(yàn)證該算法的正確性和有效性。

3.1單一高阻高速棱柱體模型

設(shè)計(jì)的第一個(gè)模型，電阻率模型與速度模型空間分布一致。大地電磁網(wǎng)格為119×53，網(wǎng)格間距為100m，地震網(wǎng)格為103×53，網(wǎng)格間距也為100m。在此需要說(shuō)明的是，大地電磁的網(wǎng)格是在地震網(wǎng)格的基礎(chǔ)上，水平方向向各自方向分別擴(kuò)展8個(gè)不等間距的網(wǎng)格。

大地電磁模型為高阻體，如圖4(a)所示，異常體電阻率值為100Ω·m，背景值為10Ω·m，頂板埋深為2 000m，底板埋深3 000m，寬度1 000m。選取不等間距的25個(gè)測(cè)點(diǎn)，均分布于地表，頻率選取9個(gè)，分別為0.010 00Hz、0.050 00Hz、0.100 00Hz、0.500 00Hz、1.00 00Hz、5.00 00Hz、10.000Hz、50.000Hz、100.00Hz。地震反射波走時(shí)反演模型為高速體，如圖4(b)所示，異常體速度值為3.5km/s，背景值為3.0km/s，頂板埋深為2 000m，底板埋深為3 000m，寬度為1 000m，選取地表的60個(gè)激發(fā)點(diǎn)，100個(gè)接收點(diǎn)。

圖5　大地電磁反演結(jié)果(電阻率單位：Ω·m)Fig.5　Inversion results for MT data

大地電磁反演的初始模型選取背景值為10Ω·m的均勻半空間，地震反射波走時(shí)反演模型選取速度背景值為3.0km/s的均勻半空間，反演迭代20次后，單一反演的結(jié)果圖分別如圖5(a)、圖6(a)所示，聯(lián)合反演結(jié)果如圖5(b)、圖6(b)所示。

從圖5(a)、 (b)和圖6(a)、 (b)可以看出，聯(lián)合反演結(jié)果比單一反演結(jié)果要好。大地電磁聯(lián)合反演高值恢復(fù)較單一反演更接近真實(shí)模型，收斂較快，且模型邊界光滑度較好；地震聯(lián)合反演速度高值恢復(fù)較好，對(duì)異常體外的速度值也有明顯的改善。

3.2單一低阻低速棱柱體模型

本文設(shè)計(jì)的第二個(gè)模型，網(wǎng)格剖分同第一個(gè)模型。速度異常體與電阻率異常體空間分布一致。

圖6　地震波速度反演結(jié)果(速度單位：km/s)Fig.6　Inversion results for seismic data

大地電磁模型為低阻體，如圖4(c)所示，頂板埋深為2 000m，底板埋深為3 000m，寬度為1 000m，異常體電阻率值為10Ω·m，背景值為100Ω·m，選取不等間距的30個(gè)測(cè)點(diǎn)，均分布于地表，頻率選取9個(gè)，分別為0.010 00Hz、0.050 00Hz、0.100 00Hz、0.500 00Hz、1.000 0Hz、5.000 0Hz、10.000Hz、50.000Hz、100.00Hz。地震模型為低速體，如圖4(d)所示，異常體速度值為2.5km/s，背景值為3.0km/s，頂板埋深為2 000m，底板埋深為3 000m，寬度為1 000m，選取地表的60個(gè)激發(fā)點(diǎn)，100個(gè)接收點(diǎn)。

大地電磁反演的初始模型選取背景值為100Ω·m的均勻半空間，地震反射波走時(shí)反演初始模型選取速度背景值為3.0km/s的均勻半空間，反演迭代30次后，單一反演的結(jié)果分別如圖5(c)、圖6(c)所示，交叉梯度聯(lián)合反演得到反演結(jié)果如圖5(d)、圖6(d)所示。

從圖5(c)、 (d)和圖6(c)、 (d)中可以看出，聯(lián)合反演結(jié)果優(yōu)于單一反演結(jié)果，低值異常體都能較好地恢復(fù)。大地電磁聯(lián)合反演結(jié)果不僅能恢復(fù)電阻率低值，還能對(duì)異常體的模型邊界有較好的恢復(fù)。地震聯(lián)合反演速度低值恢復(fù)地較好，對(duì)異常體邊界也有明顯的改善。

3.3高阻低速、低阻高速棱柱體組合模型

設(shè)計(jì)的組合模型為分開(kāi)的低阻高速體和高阻低速體，網(wǎng)格剖分同第一種模型。大地電磁模型左邊為低阻體，如圖4(e)，異常體電阻率值為1Ω·m；右邊為高阻體，異常體電阻率值為100Ω·m；背景電阻率值為10Ω·m；頂板埋深為2 000m，底板埋深為3 000m，寬度分別為500m。選取分布于地表的25個(gè)測(cè)點(diǎn)，頻率選取9個(gè)，分別為0.010 00Hz、0.050 00Hz、0.100 00Hz、0.500 00Hz、1.000 0Hz、5.000 0Hz、10.000Hz、50.000Hz、100.00Hz。地震模型左邊為高速異常體，如圖4(f)，速度值為3.5km/s；右邊為低速異常體，速度值為2.5km/s;頂板埋深為2 000m，底板埋深為3 000m，寬度分別為500m;背景速度值為3.0km/s。選取地表的60個(gè)激發(fā)點(diǎn)，100個(gè)接收點(diǎn)。

大地電磁反演的初始模型選取背景值為10Ω·m的均勻半空間，地震反射波走時(shí)的反演模型選取速度背景值為3.0km/s的均勻半空間，反演迭代20次后，單一反演的結(jié)果圖分別如圖5(e)、圖6(e)所示。經(jīng)過(guò)聯(lián)合反演計(jì)算，得到的反演結(jié)果如圖5(f)、圖6(f)所示。

從圖5(e)、(f)和圖6(e)、(f)可以看出，聯(lián)合反演效果均優(yōu)于單一反演效果。模型異常能較好地恢復(fù)，且反演邊界也與模型邊界吻合較好。聯(lián)合反演大地電磁模型中高阻異常體的電阻率與實(shí)際值偏差稍大，低阻異常體的電阻率值擬合較好，模型邊界恢復(fù)也更接近真實(shí)模型。聯(lián)合反演地震模型異常體速度值擬合較好，高速異常體外部高值異常消除，低速異常體低值恢復(fù)更好，二者的模型邊界較單一反演也有較好的改善。

4　結(jié)　論

本文對(duì)地震反射波走時(shí)和大地電磁數(shù)據(jù)的二維聯(lián)合反演進(jìn)行了較系統(tǒng)的研究，實(shí)現(xiàn)大地電磁與地震反射波走時(shí)數(shù)據(jù)的二維聯(lián)合反演算法。地震方面，基于射線追蹤正演算法和最小二乘法反演方法，實(shí)現(xiàn)正則化的加入；大地電磁方面，基于有限差分正演算法和高斯牛頓反演方法，實(shí)現(xiàn)正反演研究。加入交叉梯度函數(shù)后，設(shè)計(jì)理論地電模型和速度模型進(jìn)行了合成數(shù)據(jù)的單獨(dú)反演和聯(lián)合反演試算，成功實(shí)現(xiàn)了交叉梯度聯(lián)合反演計(jì)算。反演結(jié)果表明，聯(lián)合反演方法優(yōu)于單一反演方法，利用交叉梯度方法耦合大地電磁電阻率參數(shù)和地震波速度參數(shù)能使聯(lián)合反演過(guò)程中的這兩個(gè)模型在更新的時(shí)候相互約束、共同影響，使反演結(jié)果更接近真實(shí)模型，減少了反演解的多解性。本文成果可以為開(kāi)發(fā)其他地球物理方法聯(lián)合反演算法提供參考。

[1]GALLARDOL,MEJUM.Jointtwo-dimensionalDCresistivityandseismictraveltimeinversionwithcross-gradientsconstraints[J].JournalofGeophysicalResearchSolidEarth, 2010, 109:223-229.

[2]FREGOSOE,GALLARDOL.Cross-gradientsjoint3Dinversionwithapplicationstogravityandmagneticdata[J].Geophysics, 2009, 74(4) :32-41.

[3]HEINCKEB,HOBBSR.JointinversionofMT,gravityandseismicdataappliedtosub-basaltimaging[J].SEGExpandedAbstract, 2006, 25(1) : 784-786.

[4]JUNHOUM,CLIFFORDHThurber.Afastalgorithmfortwo-pointseismicraytracing[J].BulletinofTheSeismologicalSocietyofAmerica, 1987,77(3) :972-986.

[5]CLIFFORDHThurber.Improvedthree-dimensionalvelocitymodelsandearthquakelocationsforCalifornia[J].JournalofGeophysicalResearch, 2003,108(2):1029-2002.

[6]CLIFFORDThurber,DONNAEberhart-Phillips.Localearthquaketomographywithflexiblegridding[J].Computers&Geosciences, 1999,25(7) : 809-818.

[7]ROECKERC,THURBERD,MCGHEED.Jointinversionofgravityandarrivaltimedatafromparkfield:newconstraintsonstructureandhypocenterlocationsneartheSAFODdrillsite[J].GeophysicalResearchLetters, 2004, 31(12):399-420.

[8]于鵬,戴明剛,王家林,等. 電阻率和速度隨機(jī)分布的MT與地震聯(lián)合反演[J]. 地球物理學(xué)報(bào), 2009, 52(4):1089-1097.

[9]楊輝,王家林,吳健生,等. 大地電磁與地震資料仿真退火約束聯(lián)合反演[J]. 地球物理學(xué)報(bào), 2002,45(5):723-734.

[10]MEJUMA,GALLARDOLA,MOHAMEDAK.Evidenceforcorrelationofelectricalresistivityandseismicvelocityinheterogeneousnear—surfacematerials[J].GeophysicalResearchLetters, 2003,30(7) :1373-1376.

[11]陳曉,于鵬,張羅磊,等. 地震與大地電磁測(cè)深數(shù)據(jù)的自適應(yīng)正則化同步聯(lián)合反演[J]. 地球物理學(xué)報(bào), 2011,54(10):2673-2681.

[12]林昌洪，譚捍東，佟拓.大地電磁面積性資料和稀疏測(cè)線資料的三維反演解釋[J].現(xiàn)代地質(zhì)，2012，26(6):1185-1192.

[13]JEGENMD,RICHARDWH,PASCALT,etal.Jointinversionofmarinemagnetotelluricandgravitydataincorporatingseismicconstraints-preliminaryresultsofsub-basaltimagingofftheFaroeShelf[J].EarthandPlanetaryScienceLetters, 2009, 282(1):48-53.

[14]彭淼,譚捍東,姜枚,等. 基于交叉梯度耦合的大地電磁與地震走時(shí)資料三維聯(lián)合反演[J]. 地球物理學(xué)報(bào), 2013,56(8):2728-2738.

[15]GALLARDOLA,MEJUMA.Jointtwodimensionalcrossgradientimagingofmagnetotelluricandseismictraveltimedataforstructuralandlithologiealclassification[J].Geophysicaljournalinternational, 2007,169(3):1261-1272.

[16]周麗芬. 大地電磁與地震數(shù)據(jù)二維聯(lián)合反演研究[D].北京:中國(guó)地質(zhì)大學(xué)(北京),2012:7-36.

[17]彭淼.大地電磁與天然地震數(shù)據(jù)聯(lián)合反演研究[D] .北京:中國(guó)地質(zhì)大學(xué)(北京), 2012:25-47.

[18]劉彥,呂慶田,孟貴祥,等. 大地電磁與地震聯(lián)合反演研究現(xiàn)狀與展望[M]. 地球物理學(xué)進(jìn)展, 2012,27(6):2444-2451.

[19]陳曉,于鵬,張羅磊,等. 大地電磁與地震正則化同步聯(lián)合反演[J] . 地震地質(zhì), 2010,32(3): 402-408.

[20]張美根,程冰潔,李小凡,等. 一種最短路徑射線追蹤的快速算法[J] . 地球物理學(xué)報(bào), 2006,49(5):1467-1474.

[21]王俊,孟小紅,陳召曦,等. 交叉梯度理論及其在地球物理聯(lián)合反演中的應(yīng)用[J]. 地球物理學(xué)進(jìn)展, 2013,28(4):2094-2013.

[22]CLIFFORDHThurber.Earthquakelocationsandthree-dimensionalcrustalstructureintheCoyoteLakeArea,centralCalifornia[J].JournalofGeophysicalResearch, 1983, 88(10):8226-8236.

[23]郭繼茹,馮晅,王俊祥,等. 最佳路徑射線追蹤算法研究[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),2008,38(1):72-76.

[24]譚捍東, 陳樂(lè)壽,魏文博,等. 先進(jìn)的大地電磁資料處理和反演方法在INDEPTH-MT中的應(yīng)用研究[J]. 現(xiàn)代地質(zhì), 1997,11(3):393-400.

[25]譚捍東, 魏文博,陳樂(lè)壽,等. 大地電磁標(biāo)定曲線的畸變及校正[J]. 現(xiàn)代地質(zhì), 1998,12(4):603-606.

[26]林樹(shù)海，王偉利. 電磁波場(chǎng)數(shù)值計(jì)算中時(shí)間域有限差分法與時(shí)間域偽譜法的對(duì)比研究[J]. 現(xiàn)代地質(zhì), 2012, 26(6):1193-1198.

附公式中變量說(shuō)明：

μ0，地下介質(zhì)體的磁導(dǎo)率；E，電場(chǎng)強(qiáng)度；D，電位移矢量；B，磁感應(yīng)強(qiáng)度；H，磁場(chǎng)強(qiáng)度；▽ms，取對(duì)數(shù)的地震波速度；▽mr，取對(duì)數(shù)的電阻率；σ，介質(zhì)的電導(dǎo)率；ρ，電阻率；φ，相位；ω，角速度。

Ex，x方向電場(chǎng)強(qiáng)度；Ey，y方向電場(chǎng)強(qiáng)度；Ez，z方向電場(chǎng)強(qiáng)度；Hx，x方向磁場(chǎng)強(qiáng)度；Hy，y方向磁場(chǎng)強(qiáng)度；Hz，z方向磁場(chǎng)強(qiáng)度；mr，電阻率值；ms，速度值。

Two-dimension Joint Inversion of Magnetotelluric and Seismic Reflection Travel Time Data

CHEN Shuang1,2, TAN Handong2, GAO Jingyu2

(1.GeologicalExplorationPartyNo.208,CNNC,Baotou,InnerMongolia014010,China;2.SchoolofGeophysicsandInformationTechnology,ChinaUniversityofGeosciences,Beijing100083,China)

Because of the problem of multiple solutions for geophysical inversion, how to reduce the multiple solutions to improve the reliability and accuracy of the results of geophysical inversion is a hot geophysical research topic today. To reduce the multiple solutions of magnetotelluric and seismic data inversion, this study develops a 2-D joint inversion algorithm of magnetotelluric data and seismic data on the base of cross-gradient theory which is proposed by Gallardo and Meju, and the theoretical geoelectric model and velocity model are designed for single inversion and joint inversion, then the 2-D joint inversion algorithm of magnetotelluric and seismic data is verified to be correct and valid. The inversion results show that the above inversion algorithm with cross-gradient constraints can cause the resistivity model and velocity model constrain each other so that the joint inversion result is closer to the reality than the single inversion result. Therefore, the joint inversion is more reliable, which can reduce the multiple solution of geophysical inversion.

magnetotelluric method；seismic reflection travel time；cross-gradient method；2-D joint inversion

2015-09-14；改回日期：2016-05-05；責(zé)任編輯：潘令枝。

國(guó)家自然科學(xué)基金項(xiàng)目(41374078)；國(guó)土資源部地質(zhì)調(diào)查項(xiàng)目(12120113086100,12120113101300)。

陳霜，女，碩士研究生，工程師，1985年出生，地球物理學(xué)專業(yè)，主要從事地球物理地電學(xué)算法研究工作。Email：chens130@163.com。

譚捍東，男，1966年生，教授，地球探測(cè)與信息技術(shù)專業(yè)，主要從事電法勘探理論及應(yīng)用研究。

Email：thd@cugb.edu.cn。

P631

A

1000-8527(2016)03-0597-09