李桐林，張镕哲, 樸英哲，2

1.吉林大學地球探測科學與技術學院，長春 130026 2.金策工業(yè)綜合大學資源探測工學系，朝鮮 平壤 999093

?

大地電磁測深與地震初至波走時交叉梯度反演

李桐林1，張镕哲1, 樸英哲1，2

1.吉林大學地球探測科學與技術學院，長春 130026 2.金策工業(yè)綜合大學資源探測工學系，朝鮮 平壤 999093

在地球物理勘探過程中，觀測數(shù)據(jù)有限且存在誤差，導致反演結果存在非唯一解。為了解決非唯一性問題，筆者研究了交叉梯度數(shù)學理論，實現(xiàn)了用交叉梯度約束2D 大地電磁和地震初至波走時的聯(lián)合反演，并建立了2個模型。通過各自模型實驗和不同模型間的對比，證明了交叉梯度聯(lián)合反演比單獨反演修復幾何形態(tài)和物性參數(shù)的效果好，且交叉梯度與巖石物性無關;利用得到的交叉梯度圖進一步證明了聯(lián)合反演模型相似度高，也從交叉繪圖中得到聯(lián)合反演的物性關系更好的結論，反向證明了聯(lián)合反演的正確性。

大地電磁測深;地震初至波走時;交叉梯度;聯(lián)合反演

0 引言

單一的地球物理方法僅反應相應的物性分布，且每種地球物理方法的物理模型修復特征不同。聯(lián)合反演是結合多種地球物理場數(shù)據(jù)，通過反演地下的物性參數(shù)和幾何形態(tài)得到相同地質體[1-2]。聯(lián)合反演的參數(shù)耦合基本分為2類：1)基于物性的耦合。像直流電阻率法與大地電磁法同一物性的聯(lián)合反演，或利用經(jīng)驗及分析巖石物理關系的聯(lián)合反演，如文獻[3-4]實現(xiàn)了大地電磁測深法和重力以及地震數(shù)據(jù)間的二維聯(lián)合反演。2)空間分布結構性耦合，與物性梯度方向無關的結構約束聯(lián)合解釋是以曲率差為約束條件進行反演[5-6]。如文獻[7]實現(xiàn)了直流電阻率和地震數(shù)據(jù)的聯(lián)合反演，將交叉梯度加入到目標函數(shù)中用來約束不同參數(shù)之間的關系，該方法不依賴先驗信息即可耦合本質上不同的地球物理資料并進行解釋；文獻[8]完成了對magnetotelluric(MT)和地震間的三維交叉梯度聯(lián)合反演。

筆者先從交叉梯度數(shù)學的算法出發(fā)，對其基本理論進行剖析；再編寫大地電磁測深和地震初至波走時聯(lián)合反演程序，模擬研究大地電磁測深和地震初至波走時二維聯(lián)合反演，并將聯(lián)合反演結果和單獨反演結果進行對比。

1 基本理論和算法

1.1 交叉梯度

將交叉梯度函數(shù)應用于電阻率和速度間的聯(lián)合反演中，定義函數(shù)為

(1)

1.2 交叉梯度的離散化

交叉梯度函數(shù)是二階非線性方程，不存在不連續(xù)點和奇異點。令地質體走向沿y方向。此時模型參數(shù)對y求偏導數(shù)均為0，因此T僅在y分量上有函數(shù)值，Tx=Tz=0，則式(1)可表示為

(2)

式中，Tx、Ty、Tz分別為T在x,y,z三個方向上的分量。

地下2D網(wǎng)格化模型結構如圖1所示：Δx和Δz分別為單元的寬度和高度；mrc，mrr，mrb，msc，msr，msb分別為中心塊、右邊塊、下邊塊的電阻率的對數(shù)和慢度。采用前向差分法[10]對公式(2)進行離散化處理，可得

(3)

圖1 地下2D網(wǎng)格化模型結構Fig.1 Underground 2D grid model structure

1.3 大地電磁和地震走時的正演和雅克比矩陣

1.3.1 地震走時二維正演和雅克比矩陣

地震波走時可通過對Eikonal方程

(4)

數(shù)值求解得到。式中：t為時間；v為速度。本文采用有限差分法[11-12]求解方程(4)。

雅克比矩陣是由走時對模型慢度求導所得，對于一個給定的射線，走時的表達式為

(5)

其中：lj為第j個單元內的射線長度；msj為第j個單元內的慢度。由式(5)可見，雅克比矩陣可通過計算射線路徑求得。本文地震走時采用FAST(first arrival seismic tomography)程序，該程序可計算二、三維的初至波走時和雅克比矩陣。

1.3.2 大地電磁法二維正演和雅克比矩陣

由電磁理論可知，在非均勻地電結構中,大地電磁場由一次場和二次場組成。其中一次場為背景產(chǎn)生的場，二次場為異常體產(chǎn)生的場。由二次場的麥克斯韋方程可得兩種模式transverse electric polarization(TE)和transverse magnetic polarization(TM)的表達式，通過推導麥克斯韋方程可得兩種模式下二次場的赫姆霍茲方程。再用有限元法求解偏微分方程[13]即可實現(xiàn)大地電磁二維正演。文獻[14]給出了大地電磁二維有限元正演的算法并編程實現(xiàn)，由于正演的計算需要有限單元網(wǎng)格，因此對地下模型進行了網(wǎng)格化，且每個網(wǎng)格內的電導率是均勻的。

利用互換原理[15]計算地下電導率的大地電磁測深響應的雅克比矩陣，對赫姆霍茲方程中每個反演單元的電導率求導，結果與分布在網(wǎng)格節(jié)點新場源下的亥姆霍茲方程結果相似。根據(jù)電磁互換原理，將原正演和觀測點作為單位新場源正演求得的場值用在雅克比矩陣計算中；輔助場的雅克比矩陣是沿走向方向的場通過麥克斯韋方程的近似計算求得；視電阻率和視相位的雅克比項需要沿走向的場和輔助場的雅克比矩陣計算求得。文獻[16]給出了MT反演中求取拉格朗日乘子的方法。

在大地電磁正演和雅克比矩陣計算采用Occam2DMT(v3.0)程序[17-18]進行。

1.4 聯(lián)合反演目標函數(shù)

傳統(tǒng)的聯(lián)合反演目標函數(shù)包含對最小二乘數(shù)據(jù)的擬合和平滑約束，平滑約束幫助克服了非唯一解問題和不穩(wěn)定問題。不同的是本文增加了交叉梯度算法，該算法在每次反演迭代中，地震波速度的更新需使目標函數(shù)的數(shù)據(jù)擬合和正則化項最小化，且下次反演迭代中電阻率模型的更新需兼顧上次反演迭代中生成的地震波速度模型，直到滿足期望的擬合差。目標函數(shù)由數(shù)據(jù)擬合差和模型方差組成。定義目標函數(shù)為

(6)

約束條件為

首先，通過泰勒級數(shù)展開法處理模型正演響應和約束條件可得

(7)

其中：k為迭代次數(shù)；Ar和As分別為fr和fs的雅克比矩陣；B為交叉梯度函數(shù)的導數(shù)，在二維情況下，可通過離散化求得：

(8)

然后，將式(7)代入到式(6)可得

(9)

約束條件為

其中：

其次，利用拉格朗日乘子方法，將約束條件加入到目標函數(shù)中可得

(10)

其中，Λ為拉格朗日算子。

對式(10)求導乘以1/2得

(11)

式中，mi為第i個單元的模型參數(shù)。

最后，在約束條件下由式(11)=0，可得

(12)

將式(12)代入式(11)中可得

(13)

通過上述過程可求得目標函數(shù)結果。

1.5 程序流程

交叉梯度聯(lián)合反演主要通過兩個循環(huán)完成:內循環(huán)和外循環(huán)。外循環(huán)主要進行地震初至波走時和MT的正演、雅克比矩陣的計算，并在給定的一個較大的β初始值下逐漸減小，進行數(shù)據(jù)擬合，得到最小二乘模型。內部循環(huán)進行交叉梯度約束的聯(lián)合反演，主要計算拉格朗日算子，得到模型更新量，并對反演后的模型相似度進行判定，判斷模型是否達到擬合程度，最終得到反演模型(圖2)。

圖2 交叉梯度聯(lián)合反演算法Fig.2 Algorithm of the crossgradient joint inversion

2 交叉梯度聯(lián)合反演模型試驗

R為電阻率；vP為縱波速度。圖3 MT電阻率(a)和地震速度(b)理論模型一Fig.3 MT resistivity (a)and seismic velocity (b)theory model 1

圖4 模型一單獨反演(a、b)和聯(lián)合反演(c、d)結果Fig.4 Separate inversion (a、b) and joint inversion (c、d) in model 1

本文設計了2個理論模型，模型一如圖3所示。在MT電阻率模型(圖3a)中：左側低阻異常體電阻率為1×101.5Ω·m，其頂部和底部分別位于地下0.3和0.7 km；右側高阻高速異常體電阻率為1×104Ω·m，其頂部和底部分別位于地下0.3和1.2 km；均勻半空間背景場電阻率為1×102.5Ω·m。在地震速度模型(圖3b)中：左側為低速異常體，速度為3 000 m/s；右側為高速異常體，速度為6 000 m/s ；均勻半空間背景場速度為4 000 m/s。MT電阻率模型和地震速度模型采用相同的均勻網(wǎng)格化，水平節(jié)點125個，垂直節(jié)點61個，網(wǎng)格節(jié)點間距為50 m。MT電阻率模型觀測點為7個，間距為1 km，分布于地表；頻率為1～1 000 Hz，共10個頻率(對數(shù)間隔)。地震速度模型為地表放炮，共7個炮點，炮間距為1 km，檢波器間距為0.2 km，位于兩口井中，水平位置分別為1.5，4.5 km, 每口井內有10個檢波器。

模型一反演結果如圖4所示，其中MT和地震兩者的初始模型均為圖3理論模型中沒有異常體的背景場。由圖4a、b可見，單獨反演結果可確定異常體大概的位置，但反演異常體的準確值效果明顯一般，并且?guī)缀涡螒B(tài)也與真實的異常體(圖中白框所示)差別很大。而且，圖4b中高速體周圍產(chǎn)生的異常明顯，且深部異常體明顯差于淺層異常體的反演效果，深部反演出的異常體真實值較差。這是由射線均集中到模型的上邊中心部分，走時數(shù)據(jù)對模型的兩邊和底部不敏感所致(圖5)。由圖4c、d可見，交叉梯度聯(lián)合反演明顯優(yōu)于單獨反演，反演出的異常體幾何形態(tài)更接近真實異常體，反演出的異常體值與真實值很相似，并且地震深部的反演效果得到了提高。模型一聯(lián)合反演和單獨反演的速度和電阻率交叉梯度值如圖6所示。由圖6可見，聯(lián)合反演遠小于單獨反演的交叉梯度值，表明聯(lián)合反演后速度模型和電阻率模型的結構相似度較高。模型一電阻率和速度的交繪圖如圖7所示。由圖7可見，聯(lián)合反演推斷的物性間相互關系比單獨反演明顯。

圖5 模型一地震速度射線分布Fig.5 Distribution of rays in seismic velocity in model 1

圖6 模型一單獨反演(a)和聯(lián)合反演(b)的電阻率和速度交叉梯度值Fig.6 Cross-gradients values of separate inversion (a) and joint inversion (b) of resistivity and seismic velocity in model 1

圖8 MT電阻率(a)和地震速度(b)理論模型二 Fig.8 MT resistivity (a) and seismic velocity (b) theory model 2

圖9 模型二單獨反演(a、b)和聯(lián)合反演(c、d)結果Fig.9 Separate inversion (a、b) and joint inversion (c、d) in model 2

模型二中的異常體具有相同的速度和不同的電阻率，如圖8所示。在MT電阻率模型(圖8a)中，異常體為2個大小相同的長方形板塊體，電阻率分別為1×105Ω·m和10 Ω·m，背景場電阻率為1×103.5Ω· m，頂部和底部距地表分別為2.5 km和9 km；地表分布了間距為9 km的9個觀測測點；頻率為0.01～300 Hz，共8個頻率(對數(shù)間隔)。在地震速度模型(圖8b)中，背景場速度為4 000～5 700 m/s，從地表逐層遞增，異常體為2個相同的長方形板塊體，速度均為6 000 m/s；共10個地震炮點置于地下0.5 km處，49個間距為2 km的檢波器位于地表。模型二中2個模型采用相同的均勻網(wǎng)格化，水平節(jié)點205個，垂直節(jié)點51個，網(wǎng)格節(jié)點間距500 m。

圖10 模型二地震速度射線分布Fig.10 Distribution of rays in seismic velocity in model 2

圖11 模型二單獨反演(a)和聯(lián)合反演(b)的電阻率和速度交叉梯度值Fig.11 Cross-gradients values of separate inversion (a) and joint inversion (b) of resistivity and seismic velocity in model 2

圖12 模型二單獨反演(a)和聯(lián)合反演(b)電阻率、速度交繪圖Fig.12 Crossplots of separate inversion (a) and joint inversion (b) resistivity and seismic velocity in model 2

模型二反演結果如圖9所示，其中MT、地震初始模型均為沒有異常體的背景場。由圖9a可見，MT電阻率單獨反演結果在物性參數(shù)和幾何形態(tài)上均與實際結果有偏差，且在異常體之間和周圍圍巖中出現(xiàn)假異常體的現(xiàn)象較嚴重。圖9b可見，深部異常體明顯差于淺層異常體的反演效果，在速度單獨深部反演出的異常體真實值較差。這是由于射線均集中至模型的上邊中心部分，走時數(shù)據(jù)對模型的底部不敏感所致，如圖10所示。由圖9c和d可知，交叉梯度聯(lián)合反演明顯優(yōu)于單獨反演效果，如恢復異常體的幾何形態(tài)和物性參數(shù)，并對地震深部的異常進行恢復。

模型二中單獨反演和聯(lián)合反演的電阻率和速度交叉梯度值如圖11所示。由圖11可見，聯(lián)合反演遠小于單獨反演的交叉梯度值，表明聯(lián)合反演后不同物性模型的結構相似度更高。模型二中電阻率和速度的交叉繪圖如圖12所示。由圖12可見，聯(lián)合反演比單獨反演推斷物性間的相互關系更明顯。

3 結語

綜上，筆者研究了交叉梯度的基本理論，建立了2個模型，通過模擬實驗分別得到單獨反演和聯(lián)合反演結果，從而證明了交叉梯度適合應用于聯(lián)合反演中，可較好地改善單獨反演的幾何形態(tài)和物性參數(shù)。由2個模型實驗對比可知，電阻率和速度不是同向變化的異常體，但通過交叉梯度聯(lián)合反演可同時改善大地電磁和地震走時的單獨反演結果；表明交叉梯度聯(lián)合反演與巖石物性無關，僅取決于地下結構相似性。我們還可以通過聯(lián)合反演的交叉繪圖得到不同物性參數(shù)間關系，進一步證明了地下地質體中，不同物性具有一定的關系。將交叉梯度聯(lián)合反演方法應用于地球物理勘探中，可以更好地解決非唯一解問題和更準確地尋找地下目標體。

不足之處在于本文采用的是規(guī)則網(wǎng)格，在反演效果上差于不規(guī)則網(wǎng)格，希望在以后研究中解決上述問題。

[1] 楊輝，戴世坤，宋海斌，等.綜合地球物理聯(lián)合反演綜述[J].地球物理學進展，2002,17(2)：262-271. Yang Hui, Dai Shikun, Song Haibin, et al. Integrated Geophysical Joint Inversion Review[J]. Progress in Geophysics, 2002, 17(2): 262-271.

[2] 于鵬，王家林，吳健生，等.地球物理聯(lián)合反演的研究現(xiàn)狀和分析[J].勘探地球物理進展，2006,29(2)：87-93. Yu Peng, Wang Jialin, Wu Jiansheng, et al. Research Status and Analysis of Joint Inversion[J]. Progress in Exploration Geophysics，2006,29(2)：87-93.

[3] Heincke B, Hobbs R. Joint Inversion of MT, Gravity and Seismic Data Applied to Sub-Basalt Imaging[C]//Annual Meeting. New Orleans: SEG, 2006: 784-787.

[4] Colombo D, Stefano M D. Geophysical Modeling via Simultaneous Joint Inversion of Seismic Gravity, and Electromagnetic Data: Application to Prestack Depth Imaging[J]. The Leading Edge, 2007, 26(3): 326-331.

[5] Haber E, Oldenburg D. Joint Inversion: A Structural Approach Inverse Problems[J]. IOP Science, 1997, 13: 63-77.

[6] Zhang J, Morgan F D. Joint Seismic and Electrical Tomography[C]// Symposium on the Application of Geophysics to Engineering and Environmental Problems. [S.l.]: SEG,1997: 391-396.

[7] Gallardo L A, Meju M A. Characterization of Heterogeneous Near-Surface Materials by Joint 2D Inversion of DC Resistivity and Seismic Data[J]. Geophysical Research Letters, 2003, 30(13): 1658-1661.

[8] 彭淼，譚捍東，姜枚，等.基于交叉梯度耦合的大地電磁與地震走時資料三維聯(lián)合反演[J].地球物理學報，2013,56(8)：2728-2738. Peng Miao, Tan Handong, Jiang Mei , et al. Three-Dimension Joint Inversion Magnetotelluric and Seismic Travel Time Data with Cross-Gradient Constraints[J]. Chinese Journal of Geophysics, 2013, 56(8)：2728-2738.

[9] 王俊，孟小紅，陳朝曦，等.交叉梯度理論及其在地球物理聯(lián)合反演中的應用[J].地球物理學進展，2013, 28(4)：2094-2103. Wang Jun, Meng Xiaohong, Chen Zhaoxi, et al. The Theory of Cross-Gradient and Its Application in Geophysical Joint Inversion[J]. Progress in Geophysics, 2013, 28(4)：2094-2103.

[10] Gallardo L A, Meju M A. Joint Two-Dimensional DC Resistivity and Seismic Travel Time Inversion with Cross-Gradients Constraints[J]. Jouranl of Geophsical Research, 2004, 109: 3311.

[11] Vidale J. Finite-Difference Calculation of Travel-Ti-mes in Three Dimenisons[J]. Geophysics, 1990, 55: 521-526.

[12] Zelt C A, Barton P J. Three-Dimensional Seismic Refraction Tomography: A Comparison of Two Methods Applied to Data from the Faeroe Basin[J]. Geophys Res, 1998, 103: 7187-7210.

[13] 劉小軍，王家林，于鵬，等.基于二次場的二維大地電磁有限元法數(shù)值模擬[J].同濟大學學報：自然科學版，2007,35(8)：1113-1117. Liu Xiaojun, Wang Jialin, Yu Peng, et al. Numerical Simulation of Two-Dimensional Magnetotelluric Secondary Field Based on the Finite Element Method[J]. Journal of Tongji University：Natural Science, 2007, 35 (8)：1113-1117.

[14] Wannamaker P E, Stodt J A, Rijo L. A Stable Finite Element Solution for Two-Dimensional Magnetotelluric Modeling[J]. Geophys J R,1987, 88: 277-296.

[15] de Lugao P P, Wannamaker P E. Calculating the Two-Dimensional Magnetotlluric Jacobian in Finite Elements Using ReciProcity[J]. Geophys J,1996, 127: 806-810.

[16] 樸英哲，李桐林，劉永亮，等.在大地電磁二維Occam反演中求取拉格朗日算子方法的改進[J].吉林大學學報：地球科學版，2014，44(2)：660-667. Pak Yongchol, Li Tonglin, Liu Yongliang, et al. Improvement of Choosing Lagrange Multiplier on MT 2D Occam Inversion[J]. Journal of Jilin University : Earth Science Edition, 2014, 44(2): 660-667.

[17] de Groot-Hedlin C D, Constable S C. Occam’s Inversion to Generate Smooth, Two-Dimensional Models from Magnetotelluric Data[J]. Geophysics, 1990, 93: 1613-1624.

[18] Constable S C, Parker K L, Constable C G. Occam’s Inversion a Practical Algorithm for Generating Smooth Models from EM Sounding Data[J]. Geophysics,1987, 52: 289-300.

Joint Inversion of Magnetotelluric and First-Arrival Wave Seismic Traveltime with Cross-Gradient Constraints

Li Tonglin1,Zhang Rongzhe1, Pak Yongchol1，2

1.CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130026,China2.ResourceExplorationandEngineeringDepartment,KimchaekUniversityofTechnology,Pyongyang999093,Korea

In the process of geophysical exploration, observation data are limited and errors exist, leading to the presence of non-unique inversion results. In order to solve the problem of non-uniqueness, the authors study the mathematical theory of cross-gradient, and achieve the joint inversion of 2D magnetotelluric and first-arrival wave seismic traveltime with cross-gradient. At the same time, we establish two models. Through test of each model and comparison between the two models, it is proven not only that the result of the cross-gradient joint inversion is better than the single inversion in remediation on geometry and physical property parameters, but also that the cross-gradient is independent from the petro-physical properties. Further, the cross-gradient figure evidences the high similarity of the joint inversion models, and the cross-plot of joint inversion presents a better correlation with physical properties. These, in turn, prove the validity of the joint inversion.

magnetotelluric;first-arrival wave seismic traveltime;cross-gradient ;joint inversion

10.13278/j.cnki.jjuese.201503304.

2014-09-13

中國地質調查局項目([2013]01-060-004);國家重大科學儀器設備開發(fā)專項項目(2011YQ05006009)

李桐林(1962--)，男，教授，博士生導師，主要從事電磁法理論和應用研究，E-mail:litl@jlu.edu.cn。

10.13278/j.cnki.jjuese.201503304

P631.3

A

李桐林，張镕哲，樸英哲. 大地電磁測深與地震初至波走時交叉梯度反演.吉林大學學報:地球科學版,2015,45(3):952-961.

Li Tonglin, Zhang Rongzhe, Pak Yongchol.Joint Inversion of Magnetotelluric and First-Arrival Wave Seismic Traveltime with Cross-Gradient Constraints.Journal of Jilin University:Earth Science Edition,2015,45(3):952-961.doi:10.13278/j.cnki.jjuese.201503304.