王月 張捷

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

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

中圖分類號：P313.23，P315.2 文獻標識碼：A 文章編號：1000-0666（2018）01-0022-10

0 引言

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

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

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

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

1 研究方法

1.1 地震走時正演模擬

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

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

1.2 航空電磁正演模擬

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

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

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

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

地震走時層析成像可以解決近地表地面起伏大的問題，得到地下速度結構，但不能反演出隱藏的低速層。電磁數(shù)據(jù)反演可以得到高電阻率區(qū)和低電阻率區(qū)的分布，但反演結果的垂直分辨率較低。如果將地震走時數(shù)據(jù)和電磁數(shù)據(jù)聯(lián)合起來，同時反演地下的速度和電阻率，能夠為地球物理解釋提供更加豐富的信息，降低結果的非唯一性。

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

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

2 棋盤模型測試

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

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

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

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

3 實際數(shù)據(jù)應用

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

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

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

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

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

4 結論與討論

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

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

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

參考文獻：

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

周輝，楊寶俊，劉財.1994.地震—大地電磁測深資料綜合反演[J].地球物理學報，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.