王月 張捷

摘要：地震走時(shí)數(shù)據(jù)可以反演近地表的速度，但不能反演出隱蔽層和低速層。航空電磁數(shù)據(jù)可以反演近地表的高電阻率和低電阻率，但是對(duì)垂直方向的分辨率低。聯(lián)合三維地震走時(shí)數(shù)據(jù)和航空電磁數(shù)據(jù)反演了近地表的速度和電阻率結(jié)構(gòu)，并采用棋盤(pán)模型測(cè)試和實(shí)際數(shù)據(jù)應(yīng)用展示了聯(lián)合反演出的速度模型優(yōu)于單獨(dú)的走時(shí)反演出的速度模型。結(jié)果表明：該算法可以應(yīng)用到處理數(shù)據(jù)量大的勘探地震成像中，能提供優(yōu)化的速度結(jié)構(gòu)。

關(guān)鍵詞：聯(lián)合反演；航空電磁；地震走時(shí)；近地表

中圖分類號(hào)：P313.23，P315.2 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1000-0666（2018）01-0022-10

0 引言

地球物理反演是指利用接收到的信號(hào)推測(cè)地下介質(zhì)的物理性質(zhì)及結(jié)構(gòu)分布，不同信號(hào)能夠反演出不同的物理參數(shù)信息，為地球物理解釋提供證據(jù)，有助于尋找礦產(chǎn)、油藏等資源。近地表地質(zhì)結(jié)構(gòu)復(fù)雜，分布著沼澤、沙漠、戈壁、風(fēng)化層等松軟結(jié)構(gòu)，而且存在山地等起伏地形，給地下介質(zhì)反演帶來(lái)很大困難。因此，利用有限的數(shù)據(jù)資料反演出近地表的地質(zhì)結(jié)構(gòu)，有利于地球物理深部成像。

聯(lián)合反演在地球物理學(xué)的發(fā)展始于由Vozoff和Jupp（1975）提出的利用直流電磁測(cè)深數(shù)據(jù)和大地電磁測(cè)深數(shù)據(jù)聯(lián)合反演一維層狀電阻率模型，采用的反演算法為迭代二階馬奎特阻尼最小二乘法。隨后，聯(lián)合不同數(shù)據(jù)的反演逐漸發(fā)展，Gomez-Tre-vino和Edwards（1983）聯(lián)合可控源電磁數(shù)據(jù)和直流電測(cè)深數(shù)據(jù)利用奇異值分解法（SVD）反演了地層的電阻率和厚度，Roth和Holliger（1998）利用高分辨率的瑞雷波和導(dǎo)波聯(lián)合反演的P波速度和S波速度，Bosch和McGaughey（2001）在先驗(yàn)條件下利用重力數(shù)據(jù)和電磁數(shù)據(jù)聯(lián)合反演二維密度結(jié)構(gòu)和磁導(dǎo)率分布，Tryggvason等（2002）利用P波走時(shí)和S波走時(shí)信息聯(lián)合反演縱橫波速度和地震定位，Pieta和Bala（2008）聯(lián)合直流電測(cè)深數(shù)據(jù)和電磁數(shù)據(jù)應(yīng)用蒙特卡羅全局優(yōu)化算法實(shí)現(xiàn)了并行計(jì)算反演視電阻率等。

地震數(shù)據(jù)和非地震數(shù)據(jù)之間的聯(lián)合反演，也得到廣泛發(fā)展。周輝等（1994）利用廣義線性反演方法及非線性反演的預(yù)條件最速下降法開(kāi)展了一維地震一大地電磁測(cè)深資料反演方法研究，得出順序反演效果優(yōu)于地震、電磁單獨(dú)反演效果，而同時(shí)反演的效果最優(yōu)。王西文（1997）采用剝離法對(duì)重力、地震資料聯(lián)合反演目的層密度值，進(jìn)而預(yù)測(cè)油氣藏。Hering等（1995）利用一維直流電測(cè)深數(shù)據(jù)和地震面波數(shù)據(jù)聯(lián)合反演得到近地表的電阻率模型和面波速度模型。Misiek等（1997）利用電法數(shù)據(jù)和勒夫波數(shù)據(jù)聯(lián)合反演了Borsod地區(qū)近地表的速度和電阻率結(jié)構(gòu)。Meju和Gallardo（2003）驗(yàn)證了近地表地區(qū)電阻率和縱波速度存在結(jié)構(gòu)一致的現(xiàn)象，提出了交叉梯度算法，實(shí)現(xiàn)了二維直流電測(cè)深數(shù)據(jù)和地震走時(shí)數(shù)據(jù)的聯(lián)合反演（Gallardo，Meju，2003，2004），以及二維大地電磁和地震走時(shí)數(shù)據(jù)的聯(lián)合反演（Gallardo，Meju，2007）。交叉梯度在聯(lián)合反演中得到極大地推廣（Linde et al，2006；Fregoso，Gallardo，2009；Ogunbo，Zhang，2014）。

本文利用地震走時(shí)數(shù)據(jù)和航空電磁數(shù)據(jù)，應(yīng)用交叉梯度對(duì)模型的結(jié)構(gòu)限制，聯(lián)合反演了近地表的速度和電阻率模型。地震走時(shí)正演是三維模擬，航空電磁正演是一維模擬，通過(guò)三維Tik-honov正則化（Tikhonov，Arsenin，1977），聯(lián)合反演出三維速度模型和三維電阻率模型，再對(duì)合成棋盤(pán)模型進(jìn)行了測(cè)試，并應(yīng)用到一組實(shí)際數(shù)據(jù)中，檢驗(yàn)反演結(jié)果。

1 研究方法

1.1 地震走時(shí)正演模擬

地震波走時(shí)正演模擬的目的是計(jì)算波從源到接收器的時(shí)間以及射線路徑，主要有基于射線的方法和基于網(wǎng)格的方法?；谏渚€的方法簡(jiǎn)單高效，能解決各向同性和各向異性介質(zhì)，模擬的到時(shí)相對(duì)準(zhǔn)確，比如打靶法（Anderson，Kak，1982；Wil-lianMSon，1990）和彎曲法（Julian，Gubbins，1977；Graeber，Asch，1999）?；诰W(wǎng)格的方法種類繁多，應(yīng)用原理各不相同，可以計(jì)算模型中每個(gè)網(wǎng)格點(diǎn)的到時(shí)，能夠穩(wěn)定處理各向同性介質(zhì)，并且能同時(shí)高效計(jì)算出到時(shí)和射線路徑，如求解程函方程的有限差分法（Vidale，1988，1990）和應(yīng)用費(fèi)馬原理的最短路徑法（Nakanishi，Yamaguchi，1986；Moser，1991）。

求解程函方程和最短路徑的方法可以提供初至波到時(shí)，而打靶法和彎曲法只能得到一個(gè)不能確定是初至波的走時(shí)，因此，本文采用了求解程函方程的方法正演模擬地震走時(shí)和射線路徑。射線追蹤的程函方程為（Vidale，1988，1990）：式中：s為慢度；t為走時(shí)。此方法是采用有限差分向外擴(kuò)展的方式，利用已知點(diǎn)的到時(shí)尋找到達(dá)下一個(gè)網(wǎng)格點(diǎn)的用時(shí)最短的路徑。同時(shí)，可以根據(jù)最終走時(shí)的梯度從接收點(diǎn)回溯到源，得到射線路徑。這種方法得到的地震走時(shí)準(zhǔn)確而且計(jì)算效率高，計(jì)算速度與定義的網(wǎng)格點(diǎn)個(gè)數(shù)成正比。由于這種立方體向外擴(kuò)展計(jì)算波前的方法不總是遵循波在介質(zhì)中傳播的方向，存在波前不連續(xù)的問(wèn)題，比如，初至波可能在傳回立方體之前需要對(duì)立方體外的區(qū)域再進(jìn)行擴(kuò)展計(jì)算。但是隨著后續(xù)的發(fā)展，Vidale提出的方法在地震波走時(shí)計(jì)算中得到廣泛的應(yīng)用（VanTrier，Symes，1991；Qin et al，1992；Hole，Zelt，1995）。

1.2 航空電磁正演模擬

航空電磁包括頻率域系統(tǒng)和時(shí)間域系統(tǒng)，時(shí)間域電磁場(chǎng)通?？商峁┍阮l率域地球更深部的信息（Holladay et al，1997），而頻率域電磁場(chǎng)能得到高分辨率的淺層信息，因此，在近地表研究中，本文采用頻率域電磁場(chǎng)數(shù)據(jù)反演淺層的電阻率結(jié)構(gòu)。

時(shí)空域電磁場(chǎng)的正演計(jì)算根據(jù)麥克斯韋定律：式中：D=εE；B=μH；D為位移電流；E為電場(chǎng)；H為磁場(chǎng)；B為磁通量；σ為介質(zhì)的電導(dǎo)率；ε為介質(zhì)的電容率；μ為介質(zhì)的磁導(dǎo)率；Jm為磁流；Je為電流；r和t為空間坐標(biāo)和時(shí)間坐標(biāo)。

橫斷電場(chǎng)是由垂直磁偶極子源產(chǎn)生的（Ward，Hohmann，1988）。本文在正演模擬中采用水平圓環(huán)發(fā)射器，利用橫斷電場(chǎng)計(jì)算瞬時(shí)電壓。一維電磁場(chǎng)正演模擬計(jì)算效率高，得到了廣泛發(fā)展（Seng-piel，2010；Beard，2000；Siemon，2009）。由于近地表區(qū)域結(jié)構(gòu)復(fù)雜，橫向變化大，一維電阻率結(jié)構(gòu)不能很好地反映地下結(jié)構(gòu)，因此需要二維或者三維反演。三維電磁數(shù)據(jù)可以提供豐富的地下結(jié)構(gòu)信息，但計(jì)算量大且耗時(shí)（Cox et al，2010；Haber etal，2004）。Auken和Christiansen（2004）首次提出用一維正演模擬，加入水平方向約束反演出偽二維電阻率模型。Schultz和Ruppel（2005）通過(guò)多組一維數(shù)據(jù)，在反演目標(biāo)函數(shù)中引入二維正則化項(xiàng)，同時(shí)反演出二維電阻率模型，并應(yīng)用到實(shí)際數(shù)據(jù)中。為了提高計(jì)算效率，本文根據(jù)一維頻率域航空電磁數(shù)據(jù)的正演模擬理論，提出了應(yīng)用Tikhonov正則化矩陣加強(qiáng)相鄰網(wǎng)格點(diǎn)之間的約束，實(shí)現(xiàn)反演三維電阻率的偽三維算法。

1.3 聯(lián)合反演方法

地震走時(shí)層析成像可以解決近地表地面起伏大的問(wèn)題，得到地下速度結(jié)構(gòu)，但不能反演出隱藏的低速層。電磁數(shù)據(jù)反演可以得到高電阻率區(qū)和低電阻率區(qū)的分布，但反演結(jié)果的垂直分辨率較低。如果將地震走時(shí)數(shù)據(jù)和電磁數(shù)據(jù)聯(lián)合起來(lái)，同時(shí)反演地下的速度和電阻率，能夠?yàn)榈厍蛭锢斫忉屘峁└迂S富的信息，降低結(jié)果的非唯一性。

本文聯(lián)合三維地震走時(shí)計(jì)算和一維航空電磁模擬，通過(guò)三維Tikhonov正則化的應(yīng)用，在三維交叉梯度的結(jié)構(gòu)約束作用下，同時(shí)反演出近地表三維速度模型和電阻率結(jié)構(gòu)。聯(lián)合反演的目標(biāo)函數(shù)為：式中：mr為電阻率模型；ωr表示電磁數(shù)據(jù)的權(quán)重；dr表示觀測(cè)到的數(shù)據(jù)；Gr表示正演模擬的數(shù)據(jù)；L為Tikhonov正則化矩陣；τr為正則化項(xiàng)的平滑系數(shù)；MS為慢度模型；Gs為正演模擬的走時(shí)數(shù)據(jù)；ds為觀測(cè)到的走時(shí)數(shù)據(jù)；ωs為數(shù)據(jù)的權(quán)重；τs為模型的平滑系數(shù)；t是三維二參量的交叉梯度；τ為交叉梯度的權(quán)重。交叉梯度因子根據(jù)下列公式計(jì)算（Gallardo，Meju，2007）：

t（x，y，z）=▽log（mr（x，y，z））×▽MS（x，y，z）=txi+tyj+tzk（4）為了得到迭代反演的模型更新量，對(duì)公式（3）進(jìn)行泰勒展開(kāi)，得到：式中：β為加在聯(lián)合反演上平衡mr和MS模型更新量的參數(shù)，定義為：式中，mws和mwr分別代表模型的權(quán)重，由于電阻率的對(duì)數(shù)比慢度近似高3個(gè)量級(jí)，因此，模型權(quán)重比值乘以系數(shù)1000。三維模型反演的參數(shù)多，計(jì)算數(shù)據(jù)量大，因此利用共扼梯度算法求解公式（5），得到模型更新量。

2 棋盤(pán)模型測(cè)試

為了對(duì)比聯(lián)合反演和單獨(dú)的走時(shí)反演、電磁反演的結(jié)果，本文對(duì)棋盤(pán)模型進(jìn)行測(cè)試，真實(shí)模型如圖1所示。圖1a～d展示了不同深度速度模型的4個(gè)切面，可以看出速度模型內(nèi)部共有9個(gè)異常塊體，5個(gè)高速體和4個(gè)低速體交錯(cuò)放置。圖1e～h展示了不同深度的電阻率模型的4個(gè)切面，同樣可以看出，在電阻率模型內(nèi)部，高速體位置對(duì)應(yīng)高電阻率異常，低速體位置對(duì)應(yīng)低電阻率異常。速度模型的觀測(cè)系統(tǒng)，炮點(diǎn)在地下10m，檢波器在地表，炮間距為250m，檢波器的間距為50m，共有169個(gè)炮點(diǎn)，5929個(gè)檢波器。電阻率模型的觀測(cè)系統(tǒng)，發(fā)射接收裝置離地表50m，在X方向和Y方向的間距均是200m，每個(gè)發(fā)射器伴隨著一個(gè)接收器，二者之間的水平距離為7.9m。

首先，分別進(jìn)行單獨(dú)的地震初至波走時(shí)層析成像和偽三維航空電磁反演，將沒(méi)有異常的背景速度結(jié)構(gòu)和背景電阻率結(jié)構(gòu)作為反演的初始模型，得到三維速度模型和三維電阻率模型（圖2）。然后，將單獨(dú)反演的結(jié)果作為聯(lián)合反演的初始輸入模型。進(jìn)行聯(lián)合反演之前，我們首先需要討論在目標(biāo)函數(shù)中應(yīng)用的權(quán)重參數(shù)（τr，τs，ωs，β和γ）的選擇。根據(jù)數(shù)據(jù)和模型的二階范數(shù)的L曲線，選定速度和電阻率模型的正則化項(xiàng)平滑因子τr和τs分別為0.001和0.5。設(shè)定參數(shù)ωr=1，ωs=1，β=1，然后測(cè)試了不同的交叉梯度的權(quán)重下的聯(lián)合反演。圖3為不同γ值反演的地震走時(shí)數(shù)據(jù)不匹配度、頻率域航空電磁數(shù)據(jù)（FrequencyAirborne Electromagnetic，簡(jiǎn)稱FAEM）不匹配度、歸一化的交叉梯度值以及目標(biāo)函數(shù)值隨迭代次數(shù)的變化，交叉梯度在第三次反演迭代時(shí)引入。圖3c顯示隨著γ的增大，交叉梯度對(duì)反演的作用越來(lái)越大，電阻率模型和速度模型在結(jié)構(gòu)上越來(lái)越相似，作為補(bǔ)償，走時(shí)數(shù)據(jù)的不匹配度逐漸增加。在4組測(cè)試中，我們選擇較合理的γ=3時(shí)的反演結(jié)果展示在圖4中。

從圖4a～d中反演的速度模型中可以看出，Z為250m，300m和400m的切面中有明顯的低速體，說(shuō)明聯(lián)合反演中交叉梯度項(xiàng)起作用，使得電阻率模型的信息傳遞給了速度模型。Z為450m的切面中低速體不明顯，是由于電阻率反演在深度上的分辨率低，且航空磁測(cè)的探測(cè)深度有限，深部的電阻率反演結(jié)果和真實(shí)模型有差距，導(dǎo)致傳遞給速度模型的信息有差異，因此深部的速度模型也沒(méi)有明顯的低速結(jié)構(gòu)。

通過(guò)以上合成模型測(cè)試，可知地震走時(shí)數(shù)據(jù)和航空電磁數(shù)據(jù)聯(lián)合反演的方法能夠得到比單獨(dú)的地震走時(shí)反演近地表速度和單獨(dú)的電磁數(shù)據(jù)反演近地表電阻率更準(zhǔn)確的結(jié)果。通過(guò)交叉梯度的作用，速度模型和電阻率模型互相傳遞信息，使同一地區(qū)的2個(gè)物理特征量保持結(jié)構(gòu)上的一致。

3 實(shí)際數(shù)據(jù)應(yīng)用

我們將聯(lián)合地震走時(shí)數(shù)據(jù)和航空電磁數(shù)據(jù)反演三維速度和電阻率結(jié)構(gòu)的方法應(yīng)用到野外采集的實(shí)際數(shù)據(jù)中。加拿大阿爾伯塔地區(qū)的地表起伏嚴(yán)重，走時(shí)數(shù)據(jù)采集有難度，因此走時(shí)數(shù)據(jù)的觀測(cè)系統(tǒng)比電磁數(shù)據(jù)的觀測(cè)系統(tǒng)較稀疏。地震數(shù)據(jù)觀測(cè)系統(tǒng)在X方向覆蓋10000m，Y方向覆蓋7000m，有3542個(gè)炮點(diǎn)，4621個(gè)接收器，共接收了6355116條地震記錄，炮間距和接收器間距都為50m。電磁數(shù)據(jù)采集儀器為Fugro Airborne RESOLVE Sys-tem，觀測(cè)系統(tǒng)共有6933個(gè)發(fā)射接收裝置，每條測(cè)線的間距為100m，發(fā)射器和接收器水平距離為7.9m，高程隨地表變化，平均在地表以上35m的位置，接收到數(shù)據(jù)包括二次場(chǎng)與一次場(chǎng)的比值的實(shí)部和虛部，共有5個(gè)頻率（0.390kHz，1.787kHz，8.391kHz，40.740kHz，131.630kHz）。

首先根據(jù)地震記錄在共炮域拾取初至波走時(shí)，利用拾取到的走時(shí)文件進(jìn)行初至波走時(shí)反演。從圖5a～d反演的速度模型4個(gè)Y方向切面（4151m，5851m，7351m和8551m處）可以看出，該研究地區(qū)地表起伏大，速度分層明顯，從地表到400m處大致為一層，400～600m大致為第二層。在淺層分布著一些離散的、小的低速體，速度橫向不均勻。利用航空電磁數(shù)據(jù)進(jìn)行單獨(dú)的偽三維反演出電阻率模型，從圖5e～h可以看出，在淺層存在明顯的低電阻率層，在Y=8551m的切面中，X從60000m到63000m之間凸起地表處有明顯的高電阻率異常區(qū)。

進(jìn)一步將單獨(dú)的反演結(jié)果作為聯(lián)合反演的初始速度模型和初始電阻率模型輸入，設(shè)置不同的權(quán)重因子，比較反演結(jié)果。固定wr=1，ωs=1，β=1，測(cè)試交叉梯度權(quán)重因子γ分別為0.1、0.5、1.0和10.0作用下電磁數(shù)據(jù)不匹配度、走時(shí)數(shù)據(jù)不匹配度、歸一化的交叉梯度和目標(biāo)函數(shù)隨迭代次數(shù)的變化，如圖6所示。由圖6a看出，隨著交叉梯度權(quán)重的增大，電磁數(shù)據(jù)的不匹配度沒(méi)有明顯的變化，仍然能夠降到很小的值。由圖6b看出，當(dāng)γ<1時(shí)，走時(shí)數(shù)據(jù)的不匹配度仍然會(huì)隨著迭代次數(shù)的增加而降低，但當(dāng)γ=10時(shí)，在交叉梯度引入的第三次迭代后，走時(shí)不匹配度陡增，而且隨著反演迭代次數(shù)增加而增加，這是由于交叉梯度的權(quán)重過(guò)大，強(qiáng)烈限制了速度模型和電阻率模型在結(jié)構(gòu)上逐漸一致，而電阻率模型的分辨率較低，傳遞給速度模型的結(jié)構(gòu)信息不準(zhǔn)，所以在走時(shí)數(shù)據(jù)上會(huì)出現(xiàn)不匹配度增加的現(xiàn)象。由圖6c看出，隨著γ的增加，速度模型和電阻率模型在結(jié)構(gòu)上越來(lái)越相似，交叉梯度的值越來(lái)越小，這也驗(yàn)證了交叉梯度在聯(lián)合反演中對(duì)模型結(jié)構(gòu)上約束的作用。當(dāng)γ<1時(shí)，圖6d中目標(biāo)函數(shù)的值仍然隨著迭代次數(shù)的增加而降低，γ>1時(shí)，其值隨著迭代次數(shù)增加而增加，與走時(shí)數(shù)據(jù)的變化相同。進(jìn)一步對(duì)比γ分別為 0.1和0.5的情況，在第四次迭代之前二者數(shù)值幾乎重疊，從第五次迭代開(kāi)始，γ=0.5的曲線整體在γ=0.1曲線的下方，取值更小。也就是說(shuō)，當(dāng)γ=1時(shí)，交叉梯度的作用變大，聯(lián)合反演中模型的匹配度開(kāi)始以犧牲數(shù)據(jù)的擬合度來(lái)實(shí)現(xiàn)，速度模型和電阻率模型在結(jié)構(gòu)上的相似度較高。因此，本文認(rèn)為對(duì)這組地震走時(shí)數(shù)據(jù)和航空電磁數(shù)據(jù)的聯(lián)合反演，γ=1時(shí)反演的速度模型和電阻率模型較為可靠。

聯(lián)合反演的速度模型和電阻率模型如圖7所示。反演的速度模型中，圖7a中的淺層部分出現(xiàn)了明顯的低速塊體，在圖7d中低速體的邊緣更加清晰。反演的電阻率模型中，圖7e中淺層出現(xiàn)橫向不均勻現(xiàn)象，與速度模型有相似結(jié)構(gòu)，圖7f中在薄的低電阻率結(jié)構(gòu)下分布著不均勻的高電阻率結(jié)構(gòu)，圖7g、7h中在X為60000～62000m，凸起的地表處淺層的電阻率值偏低，這是由于在速度結(jié)構(gòu)中此區(qū)域?yàn)榈退賲^(qū)，交叉梯度的作用使得速度結(jié)構(gòu)信息傳遞給了電阻率結(jié)構(gòu)，使二者在結(jié)構(gòu)上趨于一致。

為了進(jìn)一步對(duì)比單獨(dú)反演算法和聯(lián)合反演算法得到的速度結(jié)構(gòu)的差異，我們分別在2個(gè)速度模型的基礎(chǔ)上對(duì)地震數(shù)據(jù)進(jìn)行靜校正，然后實(shí)施疊加處理，如圖8所示，箭頭指示了疊加剖面變化明顯的位置，尤其在圖8b和d中，基于聯(lián)合反演得到的速度結(jié)構(gòu)計(jì)算出來(lái)的靜校正處理得到的疊加剖面圖顯示了反射層較好的連續(xù)性和較少的起伏變化。從對(duì)比的結(jié)果可以看出，整體上聯(lián)合反演得到的速度結(jié)構(gòu)產(chǎn)生的疊加剖面要優(yōu)于單獨(dú)的走時(shí)反演得到的速度結(jié)構(gòu)產(chǎn)生的疊加剖面。因此，聯(lián)合航空電磁數(shù)據(jù)和地震走時(shí)數(shù)據(jù)反演近地表的電阻率和速度結(jié)構(gòu)可以提供更準(zhǔn)確的速度值來(lái)計(jì)算靜校正，有利于我們進(jìn)行后續(xù)的地震數(shù)據(jù)處理和地下深部成像。

4 結(jié)論與討論

本文提出了利用地震初至波走時(shí)數(shù)據(jù)和航空電磁數(shù)據(jù)聯(lián)合反演近地表三維速度和電阻率結(jié)構(gòu)的方法，該方法應(yīng)用了三維地震走時(shí)反演，一維電磁正演，結(jié)合三維Tikhonov正則化，并利用三維交叉梯度算子對(duì)結(jié)構(gòu)的約束，實(shí)現(xiàn)了偽三維反演，得到近地表的三維速度結(jié)構(gòu)和電阻率結(jié)構(gòu)。這種偽三維反演算法大大提高了計(jì)算效率，節(jié)省了存儲(chǔ)空間。電磁反演能對(duì)高電阻率區(qū)和低電阻率區(qū)成像，初至波走時(shí)反演對(duì)低速區(qū)成像有困難，而電磁數(shù)據(jù)和走時(shí)數(shù)據(jù)可以互相補(bǔ)充，聯(lián)合反演同一區(qū)域的速度和電阻率結(jié)構(gòu)，提高速度成像質(zhì)量。通過(guò)對(duì)理論模型測(cè)試和實(shí)際數(shù)據(jù)應(yīng)用，驗(yàn)證了這種聯(lián)合反演方法的可行性和結(jié)果的可靠性。

本文仔細(xì)討論了不同的權(quán)重因子在聯(lián)合反演中的作用。交叉梯度的權(quán)重過(guò)小，聯(lián)合反演的模型之間傳遞的結(jié)構(gòu)信息不足以影響聯(lián)合反演結(jié)果；交叉梯度的權(quán)重過(guò)大，數(shù)據(jù)的不匹配度增加，反演結(jié)果出現(xiàn)假象。因此，在聯(lián)合反演中，需要根據(jù)走時(shí)數(shù)據(jù)、電磁數(shù)據(jù)、交叉梯度以及目標(biāo)函數(shù)隨反演迭代次數(shù)變化的曲線判斷權(quán)重因子的合理取值。在實(shí)際數(shù)據(jù)的應(yīng)用中，本文對(duì)比了聯(lián)合反演的速度結(jié)果和單獨(dú)的走時(shí)反演的速度結(jié)果，然后分別基于這2種速度結(jié)構(gòu)對(duì)數(shù)據(jù)做了靜校正和疊加處理。疊加剖面顯示，聯(lián)合反演的速度結(jié)構(gòu)優(yōu)于單獨(dú)反演結(jié)果。因此，電磁數(shù)據(jù)和走時(shí)數(shù)據(jù)的聯(lián)合反演有助于得到更接近地下真實(shí)結(jié)構(gòu)的模型。

航空電磁數(shù)據(jù)具有采集方便、數(shù)據(jù)量大的特點(diǎn)，能夠反映地下電阻率的分布。航空電磁數(shù)據(jù)和地震走時(shí)數(shù)據(jù)聯(lián)合反演近地表的電阻率結(jié)構(gòu)和速度結(jié)構(gòu)，為地球物理解釋提供了更豐富的信息，也降低了反演的非唯一性問(wèn)題。非地震數(shù)據(jù)和地震數(shù)據(jù)的結(jié)合有效地?fù)P長(zhǎng)避短，在勘探地球物理成像領(lǐng)域具有廣泛的應(yīng)用前景。因此，在未來(lái)聯(lián)合更多的數(shù)據(jù)，比如重力數(shù)據(jù)、水文監(jiān)測(cè)數(shù)據(jù)、波形數(shù)據(jù)等，對(duì)研究地球物理成像有極大的推動(dòng)作用。

參考文獻(xiàn)：

王西文.1997.用重力、地震資料聯(lián)合反演直接預(yù)測(cè)油氣藏的方法[J].石油地球物理勘探，32（2）；221-228.

周輝，楊寶俊，劉財(cái).1994.地震—大地電磁測(cè)深資料綜合反演[J].地球物理學(xué)報(bào)，37（增刊1），471-485.

Auken E，Christiansen A V.2004.Layered and laterally constrained 2Dinversion of resistivity data[J].Geophysics，69（3）：752-761.

Anderson A H，Kak A C.1982.Digital ray tracing in two-dimensional re-fractive fields[J].Journal of the A coustical Society of America，72（5）：1593-1606.

Beard L P.2000.Comparison of methods for estimating earth resistivityfrom airborne electromagnetic measurements[J].Journal of AppliedGeophysics，45（4）：239-259.

Bosch M，McGaughey J.2001.Joint inversion of gravity and magnetic dataunder lithologic constraints[J].The Leading Edge，20（8）：877-881.

Cox L H，Wilson G A，Zhdanov M S.2010.3D inversion of airborne elec-tromagnetic data using a moving footprint[J].Exploration Geophys-ics，41（4）：250-259.

Fregoso E，Gallardo L A.2009.Cross-gradients joint 3D inversion withapplications to gravity and magnetic data[J].Geophysics，74（4）：L31-L42.

Gallardo L A，Meju M A.2003.Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismicdata[J].Geophysical Research Letters，30（13）：1658.

Gallardo L A，Meju M A.2004.Joint two-dimensional DC resistivity andseismic travel time inversion with cross-gradients constraints[J].Journal of Geophysical Research，109（B3），B03311，1-11.

Gallardo L A，Meju M A.2007.Joint 2D cross-gradients imaging of mag-netotelluric and seismic travel-time data for structural and litholog-ical classification[J].Geophysical Journal International，169（3）：1261-1272.

Graeber F M，Asch G.1999.Three-dimensional models of P wave veloci-ty and P-to-S velocity ratio in the southern central Andes by sim-ultaneous inversion of local earthquake data[J].Journal of Geophys-ical Research，104（B9）：20237-20256.

Gomez-Trevino E，Edwards R N.1983.Electromagnetic soundings in thesedimentary basin of southern Ontario-A case history[J].Geophys-ics，48（3）：311-330.

Haber E，Ascher U，Oldenburg D.2004.Inversion of 3D electromagneticdata in frequency and time domain using an inexact all-at-onceapproach[J].Geophysics，69（5）：1216-1228.

Hering A，Misiek R，Gyulai A，et al.1995.A joint inversion algorithm toprocess geoelectric and surface wave seismic data.Part Ⅰ：Basis ideas[J].Geophysical Prospecting，43（2）：135-156.

Hole J A，Zelt B C.1995.3-D finite-difference reflection travel times[J].Geophysical Journal International，121（2）：427-434.

Holladay J S，Lo B，Prinsenberg S J.1997.Bird orientation effects inquantitative airborne electromagnetic interpretation of pack ice thick-ness sounding[C].Oceans Conference，Marine Technology Society/Institute of Electrical and Electronics Engineers，Proceedings，2：1114-1119.

Julian B R，Gubbins D.1977.Three-dimensional seismic ray tracing[J].Journal of Geophysics，43：95-113.

Linde N，Binley A，Tryggvason A，et al.2006.Improved hydrogeophysicalcharacterization using joint inversion of cross-hole electrical resist-ance and ground-penetrating radar traveltime data[J].Water Re-sources Research，42（12）：W12404.

Meju M A，Gallardo L A.2003.Evidence for correlation of electrical resis-tivity and seismic velocity in heterogeneous near-surface materials[J].Geophysical Research Letters，30（7）：1373.

Misiek R，Liebig A，Gyulai A，et al.1997.A joint inversion algorithm toprocess geoelectric and surface wave seismic data.Part Ⅱ：applica-tions[J].Geophysical Prospecting，45（1）：65-85.

Moser T J.1991.Shortest path calculation of seismic rays[J].Geophys-ics，56（1）：59-67.

Nakanishi I，Yamaguchi K.1986.A numerical experiment on nonlinearimage reconstruction from first-arrival times for two-dimensionalisland arc structure[J].Journal of Physics of the Earth，34（2）：195-201.

Ogunbo J N，Zhang J.2014.Joint seismic traveltime and TEM inversionfor near surface imaging[C].SEG：84th Annual International Meet-ing，2104-2108.

Pieta A，Bala J.2008.Joint inversion of direct current（DC）and electro-magnetic（EM）measurements in parallel computing environment[C].Barcelona；14th EAGE European Meeting of Environmentaland Engineering Geophysics.

Qin F，Luo Y，Olsen K B，et al.1992.Finite-difference solution of theeikonal equation along expanding wavefronts[J].Geophysics，57（3）.478-487.

Roth M，Holliger K.1998.Joint inversion of Rayleigh and guided waves inhigh-resolution seismic data using a genetic algorithm[C].SEG：67th Annual International Meeting，1570-1573.

Schultz G，Ruppel C.2005.Inversion of inductive electromagnetic data inhighly conductive terrains[J].Geophysics，70（1）：G16-G28.

Sengpiel K-P.2010.Approximate inversion of airborne EM data from amulti-layered ground[J].Geophysical Prospecting，36（4）：446-459.

Siemon B.2009.Levelling of helicopter-borne frequency-domain elec-tromagnetic data[J].Journal of Applied Geophysics，67（3）：206-218.

Tikhonov A N，Arsenin V Y.1977.Solutions of Ill-Posed Problems[M].New York：Wiley.

Tryggvason A，Rognvaldsson S Th，F(xiàn)lovenz G.2002.Three-dimensionalimaging of the P-and S-wave velocity structure and earthquake lo-cations beneath Southwest Iceland[J].Geophysical Journal Interna-tional，151（3）：848-866.

Van Trier J，Symes W W.1991.Unwind finite-difference calculation oftraveltimes[J].Geophysics，56（6）：812-821.

Vidae J E.1990.Finite-difference calculations of traveltimes in three di-mensions[J].Geophysics，55（5）：521-516.

Vidale J E.1988.Finite-difference calculations of traveltimes[J].Bulle-tin of the Seismological America，78（6）：2062-2076.

Vozoff K，Jupp D L B.1975.Joint inversion of geophysical data[J].Geo-physical Journal International，42（3）：977-991.

Ward S H，Hohmann G W.1988.Electromagnetic theory for geophysicalapplications[M].SEG，131-311.

William P R.1990.Tomographic inversion in reflection seismology[J].Geophysical Journal International，100（2）：255-274.