基于寬范圍巖石物性約束的大地電磁和地震聯(lián)合反演

陳曉, 于鵬， 鄧居智， 張磊, 葛坤朋， 王彥國

1 東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室， 南昌 330013 2 中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074 3 同濟(jì)大學(xué)海洋地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室， 上海 200092

?

基于寬范圍巖石物性約束的大地電磁和地震聯(lián)合反演

陳曉1，2, 于鵬3*， 鄧居智1， 張磊1, 葛坤朋1， 王彥國1

1 東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室， 南昌 330013 2 中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074 3 同濟(jì)大學(xué)海洋地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室， 上海 200092

巖石物性約束是地球物理聯(lián)合反演重要的實(shí)現(xiàn)方式.在以往聯(lián)合反演研究中，這種約束通常以點(diǎn)對點(diǎn)的映射方式出現(xiàn)，要求經(jīng)驗(yàn)性的巖石物性關(guān)聯(lián)是相對精確的，無形中提高了巖石物性約束應(yīng)用于聯(lián)合反演的條件.基于此，本文結(jié)合全局尋優(yōu)的模擬退火算法，提出了寬范圍巖石物性約束技術(shù).該技術(shù)將巖石物性約束與模擬退火模型擾動、模型接受準(zhǔn)則相結(jié)合，實(shí)現(xiàn)了空間對空間的映射，將電阻率和速度耦合在一起.文章以MT和地震正則化同步聯(lián)合反演為例，檢驗(yàn)了不同的先驗(yàn)信息條件下的寬范圍巖石物性約束技術(shù)的適用性.模型試驗(yàn)表明：新技術(shù)容易實(shí)現(xiàn)，容錯性高，可降低“不精確”巖石物性關(guān)聯(lián)信息引入帶來的風(fēng)險，提高巖石物性約束應(yīng)用的靈活性.

聯(lián)合反演； 巖石物性約束； 大地電磁； 地震； 模擬退火

1 引言

在當(dāng)前地球物理聯(lián)合反演研究中，不同巖石物性參數(shù)之間的耦合方式主要有兩種：一種為構(gòu)造幾何約束，要求不同的巖石物性參數(shù)的分布或變化特征在空間上保持一定的一致性；另一種是巖石物性約束，要求不同的巖石物性參數(shù)可以通過某種經(jīng)驗(yàn)性、統(tǒng)計性的巖石物性關(guān)聯(lián)耦合起來.

作者曾在多篇文獻(xiàn)(Chen et al.,2010,2015;陳曉等，2011a; 陳曉，2013)中從不同角度分析過MT和地震聯(lián)合反演的發(fā)展概況和趨勢.本文僅從巖石物性約束角度，回顧一下電阻率和速度耦合在MT和地震聯(lián)合反演中的應(yīng)用.Dell′Aversana(2001)利用測井資料確定了電阻率、速度和密度的經(jīng)驗(yàn)關(guān)系，運(yùn)用順序反演的方法研究了地質(zhì)情況復(fù)雜的逆沖斷層帶.Heincke等(2006,2010,2014)利用巖石物性約束，實(shí)現(xiàn)了MT、重力以及地震數(shù)據(jù)的聯(lián)合反演，并指出聯(lián)合反演可以確定出單一方法無法確定的玄武巖下伏地層的空間分布.Colombo和De Stefano(2007)通過測井資料確定了電阻率和速度的經(jīng)驗(yàn)關(guān)系，利用地震資料聯(lián)合重力和MT數(shù)據(jù)重建速度分布，進(jìn)而進(jìn)行疊前深度成像,并指出聯(lián)合反演可以減少解的多解性并提高非地震方法的分辨率，非地震方法可以為低信噪比的地震資料提供有益的補(bǔ)充.Jegen等(2009)利用大洋鉆探資料，以電阻率、速度和密度的巖石物性關(guān)聯(lián)為出發(fā)點(diǎn)，實(shí)現(xiàn)了MT、重力和地震數(shù)據(jù)的聯(lián)合反演，獲得了玄武巖下伏地層的分布.在國內(nèi)，陳曉等(2011b)將電阻率和速度的巖石物性關(guān)聯(lián)引入MT和地震的正則化同步聯(lián)合反演中，檢驗(yàn)了巖石物性約束在MT和地震聯(lián)合反演中的可行性.構(gòu)造幾何約束靈活方便，對先驗(yàn)信息要求相對較少.交叉梯度(Gallardo and Meju,2004; Gallardo et al., 2012) 是一種典型的構(gòu)造幾何約束方式.有學(xué)者(Moorkamp et al.,2011; Lelièvre et al.,2012)研究了多種聯(lián)合反演耦合方式，指出與交叉梯度約束相比，巖石物性約束可更好地促進(jìn)巖石物性耦合,但這種約束方式的適用性較弱.

作者認(rèn)為開展基于巖石物性約束的聯(lián)合反演研究存在以下幾點(diǎn)制約：第一，巖石物性關(guān)聯(lián)不易建立，不是所有的勘探區(qū)域都有詳細(xì)的測井或巖石物性統(tǒng)計資料；第二，傳統(tǒng)的點(diǎn)對點(diǎn)映射式的巖石物性關(guān)聯(lián)不夠精確，一條擬合曲線不能完全反映巖石物性參數(shù)之間的空間對應(yīng)關(guān)系；第三，巖石物性關(guān)聯(lián)區(qū)域特征明顯，適應(yīng)性有限，一或幾“孔”所得的巖石物性關(guān)聯(lián)容錯率低，用于整條剖面乃至整個研究區(qū)是不嚴(yán)謹(jǐn)?shù)?基于此，本文以MT和地震聯(lián)合反演為例，嘗試探索一種容易實(shí)現(xiàn)、空間對空間、容錯率高的巖石物性關(guān)聯(lián)技術(shù)，進(jìn)而提高巖石物性約束方式應(yīng)用的靈活性.

2 電阻率和速度巖石物性關(guān)聯(lián)

2.1 傳統(tǒng)的巖石物性約束

電阻率和速度之間的巖石物性關(guān)聯(lián)很難用一個明確的、通用的表達(dá)式來確定.雖然Faust公式(格勞爾,1987)可以將速度、電阻率和埋深聯(lián)系在一起，但不斷有學(xué)者(徐群洲等，2000；嚴(yán)良俊等，2004)對Faust公式提出疑問：Faust公式是在深層地層中得到的關(guān)系，是否具有廣泛的適用性？因此鮮見文獻(xiàn)利用Faust公式開展MT和地震聯(lián)合反演研究.從作者所查閱的文獻(xiàn)來看：根據(jù)測井鉆孔資料擬合出巖石物性關(guān)聯(lián)曲線是耦合電阻率和速度的重要途徑，且速度或速度的對數(shù)，和電阻率或電阻率的對數(shù)之間滿足某種線性或非線性的經(jīng)驗(yàn)性關(guān)聯(lián).

需要指出，測井鉆孔資料所展示的巖石物性關(guān)聯(lián)特征并不是一一對應(yīng)的，通常一種巖石物性參數(shù)(比如：速度)的取值對應(yīng)了另外一種巖石物性參數(shù)(比如：電阻率)一個可能的分布空間，即真實(shí)的巖石物性關(guān)聯(lián)特征應(yīng)該是空間對空間的.但在以往應(yīng)用中，學(xué)者們大多利用的是一個線性或非線性的關(guān)聯(lián)公式將一種巖石物性參數(shù)直接映射到另外一種巖石物性參數(shù)，即利用的是一種點(diǎn)對點(diǎn)的映射.雖然這種映射在一定程度上反映了巖石物性參數(shù)的分布特征，但無法完整地反映巖石物性參數(shù)之間的空間對應(yīng)關(guān)系.因此利用點(diǎn)對點(diǎn)映射式的巖石物性關(guān)聯(lián)來耦合兩種巖石物性參數(shù)，進(jìn)而縮小聯(lián)合反演的解空間，是以拋棄部分解空間為代價的，當(dāng)經(jīng)驗(yàn)性的巖石物性關(guān)聯(lián)“不精確”時，被拋棄的解空間反而可能包含真實(shí)解.

2.2 寬范圍巖石物性約束技術(shù)

針對傳統(tǒng)的巖石物性約束技術(shù)的不足，作者結(jié)合非?？焖俚哪M退火算法(VFSA)，提出了寬范圍巖石物性約束技術(shù).VFSA是一種常見的全局尋優(yōu)的算法，其基本原理、技術(shù)細(xì)節(jié)請參見文獻(xiàn)(楊輝等，2002；于鵬等，2009；Chen et al.,2011).下面以速度映射電阻率為例來介紹新技術(shù)的具體流程：

(1) 在某一溫度，利用(1)式(Ingber,1989)隨機(jī)產(chǎn)生與溫度有關(guān)的物性模型：

(1)

(2) 根據(jù)已有的經(jīng)驗(yàn)性的速度-電阻率巖石物性關(guān)聯(lián)，由速度映射出電阻率，這一步和傳統(tǒng)的巖石物性約束方式相同.

(3) 不直接將映射所得的電阻率值進(jìn)行正、反演運(yùn)算，而是將映射出的電阻率進(jìn)行上下浮動，獲得一個新的電阻率分布區(qū)間[Ai,Bi].

(4) 根據(jù)新的電阻率擾動區(qū)間，再根據(jù)VFSA的模型擾動方式(公式(1))，擾動出新的電阻率.

(5) 將擾動所得速度、電阻率模型進(jìn)行正演，然后根據(jù)VFSA的模型接受準(zhǔn)則來判斷是否接受新的速度、電阻率模型.

(6) 假如新的速度、電阻率模型被接受，則成功地完成一次速度到電阻率的映射，反之，當(dāng)前溫度當(dāng)次模型擾動失敗，進(jìn)行下一次擾動.重復(fù)步驟(1)—(6)直至達(dá)到當(dāng)前溫度下的最大擾動次數(shù).

(7) 降低溫度，重復(fù)步驟(1)—(7)直至滿足模擬退火迭代終止條件.

寬范圍巖石物性約束技術(shù)的核心是不再簡單地將點(diǎn)對點(diǎn)映射所得的巖石物性參數(shù)直接代入正、反演運(yùn)算，而是對所獲得巖石物性參數(shù)進(jìn)行再次加工，獲得新的巖石物性分布空間，并結(jié)合全局尋優(yōu)的VFSA算法進(jìn)行重新搜索，進(jìn)而實(shí)現(xiàn)電阻率和速度參數(shù)的耦合.

3 MT和地震聯(lián)合反演的實(shí)現(xiàn)

MT和地震正則化同步聯(lián)合反演的目標(biāo)函數(shù)由數(shù)據(jù)擬合和模型約束兩部分組成(Zhdanov, 2002；Chen et al., 2010)：

Pα(m) =φ(m)+αS(m)=WMTφMT+WSφS+αS(m)

=min,

(2)

其中，m為待解模型，Pα(m)為目標(biāo)函數(shù)，φ(m)為數(shù)據(jù)擬合函數(shù)，φMT是視電阻率或相位的觀測值和計算值的均方根誤差，φS是地震走時的觀測值和計算值的均方根誤差，WMT、WS為兩方法的權(quán)重，S(m)為模型穩(wěn)定函數(shù)，本文選擇最小模型約束，α為正則化因子.

MT二維正演采用通用的有限單元法(陳樂壽和王光鍔，1990)，地震正演采用簡單而快速精確的二維速度隨機(jī)分布的射線追蹤方法(高爾根和徐果明，1996)，具體正演公式，自適應(yīng)正則化算法、物性參數(shù)隨機(jī)分布的建模方法和VFSA技術(shù)細(xì)節(jié)，在相關(guān)文獻(xiàn)中(于鵬等，2009；陳曉等，2011a；Chen et al.,2011；陳曉，2013)有詳細(xì)描述，本文不再贅述.

4 模型試驗(yàn)

為了模擬經(jīng)驗(yàn)性的巖石物性關(guān)聯(lián)不一定存在、存在也不一定精確這些在實(shí)際資料處理中經(jīng)常遇到的情況，本文設(shè)計了三個模型試驗(yàn)來檢驗(yàn)新技術(shù)的適用性.

4.1 相對精確的巖石物性關(guān)聯(lián)條件

設(shè)計了如圖1所示的地球物理模型.四套地層自上而下的電阻率分別為：5、10、300和1000 Ωm，異常塊體A和B的電阻率分別為30和15 Ωm；四套地層自上而下的速度分別為：3.2、3.6、5.0和5.5 km·s-1，異常塊體A和B的速度分別為4.3和4.0 km·s-1.反演目標(biāo)是異常體A和B，即電阻率差異僅為15 Ωm的兩個低阻異常.MT模擬測點(diǎn)共33個，測點(diǎn)間距為1 km.MT計算頻點(diǎn)在320～0.00055 Hz之間取40個.

作者借鑒文獻(xiàn)(Dell′Aversana P,2001)中的電阻率和速度之間的關(guān)聯(lián)表達(dá)式，結(jié)合上述地球物理模型，利用線性回歸方法，重新計算了相關(guān)待定系數(shù)(公式(3))，得到如圖2所示的相對精確的巖石物性關(guān)聯(lián).

ln(lnρ)= 0.646×(V-2.405),

(3)

圖1 模型試驗(yàn)1設(shè)計的模型(色標(biāo)同圖3)(a) 電阻率模型； (b) 速度模型.Fig.1 Designed models for model test 1 (color scale is the same as the Fig.3)(a) Resistivity model; (b) Velocity model.

圖2 模型試驗(yàn)1的速度和電阻率物性關(guān)聯(lián)Fig.2 Relationship between velocity and resistivity for model test 1

其中，V為速度，ρ為電阻率.

模型試驗(yàn)中，作者設(shè)定異常塊體A、B存在相同的電阻率擾動空間(10～40 Ωm)，并將兩個異常塊體的頂界面和物性分布予以松約束，即放開反演，其余反演參數(shù)予以緊約束.異常塊體A、B的速度變化空間分別為：4.1～4.5 km·s-1，3.8～4.2 km·s-1.異常塊體A、B頂界面的變化空間為其真實(shí)數(shù)值的±30%.

本文設(shè)計了三種聯(lián)合反演方案：方案一，無巖石物性約束的聯(lián)合反演；方案二，基于點(diǎn)對點(diǎn)映射的巖石物性約束的聯(lián)合反演；方案三，基于寬范圍巖石物性約束的聯(lián)合反演，其中速度映射出電阻率重新擾動±20%作為電阻率再次搜索所需要的物性范圍.模擬退火初始擾動結(jié)果見圖3(a，b)，三種方案的反演結(jié)果見圖3(c—h)，迭代誤差曲線見圖3i,聯(lián)合反演結(jié)果的耦合情況見圖3j.聯(lián)合反演結(jié)果的定量統(tǒng)計見表1.

分析圖3和表1，可見無巖石物性約束的聯(lián)合反演(方案一)對兩套低阻層識別不清，其余兩種基于巖石物性約束的聯(lián)合反演(方案二、三)都可以較好地區(qū)分這兩套低阻層.分析圖3i，可見與方案一相比，方案二、三迭代效率更高.分析圖3j，可見通過巖石物性約束電阻率和速度耦合在一起，可以有效地縮小聯(lián)合反演的解空間，反演結(jié)果更接近真實(shí)值.根據(jù)表1，可以看出，在巖石物性經(jīng)驗(yàn)關(guān)系相對精確的條件下，傳統(tǒng)的點(diǎn)對點(diǎn)映射約束和寬范圍巖石物性約束的聯(lián)合反演效果相當(dāng)，均可以提高解的分辨能力，優(yōu)于無巖石物性約束的聯(lián)合反演效果.

4.2 “不精確”的巖石物性關(guān)聯(lián)條件

為了檢驗(yàn)新的巖石物性約束技術(shù)在先驗(yàn)的巖石物性關(guān)聯(lián)信息“不精確”時的適用性，作者對4.1節(jié)中的模型進(jìn)行修改(圖4)，增加了一個尖滅層，且使目標(biāo)地層分布得更深.此外作者調(diào)整了目標(biāo)體的巖石物性取值使其與公式(3)存在一定偏差，具體如下：異常塊體A、B電阻率分別為50和15 Ωm,尖滅層電阻率為80 Ωm，尖滅層速度為4.5 km·s-1.

圖3 模型試驗(yàn)1不同方案聯(lián)合反演結(jié)果(色標(biāo)同圖1)(a) 電阻率初始擾動模型； (b) 速度初始擾動模型； (c) 方案一電阻率反演結(jié)果； (d) 方案一速度反演結(jié)果； (e) 方案二電阻率反演結(jié)果； (f) 方案二速度反演結(jié)果； (g) 方案三電阻率反演結(jié)果； (h) 方案三速度反演結(jié)果； (i) 迭代誤差曲線圖； (j) 聯(lián)合反演結(jié)果耦合圖(黑線為公式(3)).Fig.3 Joint inversion results by different schemes for model test 1 (color scale is the same as the Fig.1)(a) Initial perturbed resistivity model; (b) Initial perturbed velocity model; (c) Resistivity inversion by scheme 1; (d) Velocity inversion by scheme 1; (e) Resistivity inversion by scheme 2; (f) Velocity inversion by scheme 2; (g) Resistivity inversion by scheme 3; (h) Velocity inversion by scheme 3; (i) Curves of iteration errors; (j) Coupling diagram of joint inversion (black line represents Equation (3)).

圖4 模型試驗(yàn)2設(shè)計的模型(色標(biāo)同圖6和8)(a) 電阻率模型； (b) 速度模型.Fig.4 Designed models for model test 2 (color scale is same as Figs.6 and 8)(a) Resistivity model; (b) Velocity model.

圖5 模型試驗(yàn)2的速度和電阻率物性關(guān)聯(lián)Fig.5 Relationship between velocity and resistivity for model test 2

如圖5所示，公式(3)此時不是嚴(yán)格精確的.模型試驗(yàn)中，作者設(shè)定異常塊體A、B存在相同的電阻率擾動空間(5～60 Ωm).此外，將兩個異常塊體的頂界面和物性分布，以及尖滅層的物性分布予以松約束，即放開反演.其余反演參數(shù)予以緊約束.異常塊體A、B以及尖滅層的速度變化空間分別為：4.1～4.5 km·s-1，3.8～4.2 km·s-1，4.3～4.7 km·s-1，電阻率變化空間分別為：5～60 Ωm，5～60 Ωm，40～120 Ωm.異常塊體A、B頂界面的變化空間為其真實(shí)數(shù)值的±30%.

本文設(shè)計了三種聯(lián)合反演方案：方案一，無巖石物性約束的聯(lián)合反演；方案二，基于點(diǎn)對點(diǎn)映射的巖石物性約束的聯(lián)合反演；方案三，基于寬范圍巖石物性約束的聯(lián)合反演，并將由速度映射出電阻率重新擾動±80%作為電阻率再次搜索所需要的物性范圍.模擬退火初始擾動結(jié)果見圖6(a，b)，三種方案的反演結(jié)果見6(c—h)，迭代誤差曲線見圖6i,聯(lián)合反演結(jié)果的耦合情況見圖6j.聯(lián)合反演結(jié)果的定量統(tǒng)計見表2.

對比圖6(c—h)，可見無巖石物性約束的聯(lián)合反演(方案一)對兩套低阻層的區(qū)分效果不佳，電阻率高低相間分布，不能還原真實(shí)的電阻率模型，考慮巖石物性的兩種聯(lián)合反演(方案二、三)均可以將兩套低阻層區(qū)分開來；分析圖6i，可見基于巖石物性約束的聯(lián)合反演迭代效率更高，對比圖3i,可見當(dāng)先驗(yàn)的巖石物性關(guān)聯(lián)“不精確”時，方案二迭代效率要低于方案三.

分析表2，可見方案二、三對兩套低阻層的刻畫明顯優(yōu)于方案一，但方案二在尖滅層處的速度反演結(jié)果偏高，電阻率反演結(jié)果偏低，這是引用“不精確”巖石物性關(guān)聯(lián)所帶來的不穩(wěn)定因素造成的，由于寬范圍巖石物性約束技術(shù)對“不精確”巖石物性關(guān)聯(lián)的映射結(jié)果進(jìn)行了重新擾動、搜索和判斷，降低了“不精確”先驗(yàn)信息的束縛.

表2 模型試驗(yàn)2反演結(jié)果定量統(tǒng)計表

圖6 模型試驗(yàn)2不同方案聯(lián)合反演結(jié)果(色標(biāo)同圖4和8)(a) 電阻率的初始擾動模型； (b) 速度的初始擾動模型； (c) 方案一電阻率反演結(jié)果； (d) 方案一速度反演結(jié)果； (e) 方案二電阻率反演結(jié)果； (f) 方案二速度反演結(jié)果； (g) 方案三電阻率反演結(jié)果； (h) 方案三速度反演結(jié)果； (i) 迭代誤差曲線圖； (j) 聯(lián)合反演結(jié)果的物性耦合圖(黑線為公式(3)).Fig.6 Joint inversion results by different schemes for model test 2 (color scale is the same as the Figs.4 and 8)(a) Initial perturbed resistivity model; (b) Initial perturbed velocity model; (c) Resistivity inversion by scheme 1; (d) Velocity inversion by scheme 1; (e) Resistivity inversion by scheme 2; (f) Velocity inversion by scheme 2; (g) Resistivity inversion by scheme 3; (h) Velocity inversion by scheme 3; (i) Curves of iteration errors; (j) Coupling diagram of joint inversion (black line represents Equation (3)).

圖7 不同迭代次數(shù)對應(yīng)的反演結(jié)果的物性耦合(黑線為公式(3))(a) 1次迭代； (b) 50次迭代； (c) 100次迭代； (d) 250次迭代.Fig.7 Coupling relationships of inversion results to different iterative numbers (black line represents Equation (3))(a) 1 iteration; (b) 50 iterations; (c) 100 iterations; (d) 250 iterations.

圖7展示了方案三第1、50、100和250次迭代所接受的電阻率和速度模型的巖石物性耦合情況.可以看出在寬范圍巖石物性約束技術(shù)約束下，電阻率和速度呈現(xiàn)一種空間對空間的映射關(guān)系.迭代初期，模擬退火有較高的概率接受高能量狀態(tài)，接受的模型與真實(shí)的模型偏差較大，模型參數(shù)分布在“不精確”的物性關(guān)聯(lián)附近，隨著反演迭代的進(jìn)行，目標(biāo)函數(shù)的下降，模擬退火算法不斷篩選新模型，模型參數(shù)分布逐漸偏離“不精確”的物性關(guān)聯(lián)，接受的模型逐漸接近真實(shí)值.圖4清晰地展示出：當(dāng)先驗(yàn)的巖石物性關(guān)聯(lián)條件“不精確”時，新的巖石約束技術(shù)結(jié)合VFSA算法的模型擾動和接受準(zhǔn)則，確定巖石物性參數(shù)的空間分布，進(jìn)而促進(jìn)巖石物性耦合是可行的.

通過該模型試驗(yàn)，可以看出與以往點(diǎn)對點(diǎn)映射式的巖石物性關(guān)聯(lián)方式相比，寬范圍巖石物性約束技術(shù)具有較高的容錯率，可適用于先驗(yàn)的巖石物性關(guān)聯(lián)“不精確”的條件.

4.3 不存在巖石物性關(guān)聯(lián)條件

在4.1節(jié)和4.2節(jié)中，寬范圍巖石物性約束技術(shù)的應(yīng)用是以存在先驗(yàn)的巖石物性經(jīng)驗(yàn)關(guān)系為前提的.但假如連這個條件都不存在或者不易獲得呢？即電阻率和速度不存在明確的巖石物性經(jīng)驗(yàn)關(guān)系.基于這種條件，對上文提出的巖石物性約束技術(shù)進(jìn)行拓展：首先根據(jù)各套地層的電阻率和速度的可能分布范圍，建立簡單靈活的線性關(guān)系，以此來代替經(jīng)驗(yàn)性的巖石物性關(guān)聯(lián)，然后套用寬范圍巖石物性約束技術(shù)的相關(guān)流程.

選擇4.2節(jié)中的模型，反演過程中各種參數(shù)的選擇(包括物性和界面的擾動空間)和4.2節(jié)完全相同.但假設(shè)這次模型試驗(yàn)不存在先驗(yàn)的巖石物性經(jīng)驗(yàn)關(guān)系，僅有以下先驗(yàn)信息：目標(biāo)層的電阻率變化大體范圍分別為：5～60 Ωm，5～60 Ωm，40～120 Ωm；目標(biāo)層速度變化大體范圍分別為：4.1～4.5 km·s-1，3.8～4.2 km·s-1，4.3～4.7 km·s-1；目標(biāo)層的電阻率和速度大體為正比例關(guān)聯(lián).

首先，根據(jù)目標(biāo)層物性參數(shù)可能的分布范圍，運(yùn)用兩點(diǎn)公式，建立線性的巖石物性關(guān)聯(lián)式.考慮到異常塊體A、B和尖滅層的空間位置，建立兩套巖石物性關(guān)聯(lián).對于異常塊體A和B，巖石物性參數(shù)的上限為lnln60 Ωm,4.5 km·s-1,下限為lnln5 Ωm,3.8 km·s-1，套用兩點(diǎn)公式，電阻率和速度滿足：

ln(lnρ)=1.334×(V-3.8)+0.476.

(4)

對于尖滅層，巖石物性參數(shù)的上限為lnln120 Ωm,4.7 km·s-1,下限為lnln40 Ωm,4.3 km·s-1，套用兩點(diǎn)公式，電阻率和速度滿足：

ln(lnρ)=0.6525×(V-4.3)+1.305.

(5)

然后，進(jìn)入寬范圍巖石物性約束流程，將由速度映射出電阻率重新擾動±80%，作為電阻率再次搜索所需要的物性范圍.電阻率和速度聯(lián)合反演結(jié)果見圖8(a，b)，聯(lián)合反演結(jié)果的耦合情況見圖8c.聯(lián)合反演結(jié)果的定量統(tǒng)計見表3.

本次聯(lián)合反演結(jié)果要遜于4.2節(jié)中寬范圍巖石物性約束的聯(lián)合反演結(jié)果，這是因?yàn)楸敬文Ｐ驮囼?yàn)只采用了一個根據(jù)目標(biāo)層大體物性分布范圍快速建立的、簡單的巖石物性關(guān)聯(lián)，這樣的“經(jīng)驗(yàn)關(guān)系”勢必是不精確的，但與4.2節(jié)中無巖石物性約束的聯(lián)合反演結(jié)果(圖6c，6d)相比，本次聯(lián)合反演依然可以清晰地揭示目標(biāo)層電阻率的高低分布.

通過建立一個大體的物性關(guān)聯(lián)，進(jìn)而套用寬范圍巖石物性約束技術(shù)，是寬范圍巖石物性約束技術(shù)的拓寬，不必要求“精確”或“不精確”的巖石物性經(jīng)驗(yàn)關(guān)系的存在，只需知道目標(biāo)層的巖石物性參數(shù)的大體分布范圍.

圖8 模型試驗(yàn)3聯(lián)合反演結(jié)果(色標(biāo)同圖4和圖6)(a) 電阻率反演結(jié)果； (b) 速度反演結(jié)果； (c) 聯(lián)合反演結(jié)果的物性耦合圖.Fig.8 Joint inversion results for model test 3 (color scale is the same as the Figs.4 and 6)(a) Resistivity inversion; (b) Velocity inversion; (c) Coupling of physical properties in joint inversion results.

反演目標(biāo)反演參數(shù)反演參數(shù)的均值和標(biāo)準(zhǔn)差名稱真實(shí)值平均值標(biāo)準(zhǔn)差異常塊體A電阻率(Ωm)5028.5146.270速度(km·s-1)4.34.3010.032頂界面深度(km)-2.5-2.5400.041異常塊體B電阻率(Ωm)1510.1712.173速度(km·s-1)4.04.0120.048頂界面深度(km)-3.25-3.4130.109尖滅體電阻率(Ωm)8073.41418.358速度(km·s-1)4.54.4950.035

簡而言之，本文提出的新技術(shù)，是利用先驗(yàn)的巖石物性參數(shù)分布范圍構(gòu)建某種巖石物性關(guān)聯(lián)，再將通過巖石物性關(guān)聯(lián)映射所得巖石物性參數(shù)進(jìn)行再擾動，結(jié)合模擬退火算法的模型產(chǎn)生和模型接受準(zhǔn)則進(jìn)而實(shí)現(xiàn)寬范圍巖石物性約束.根據(jù)先驗(yàn)信息確定的巖石物性關(guān)聯(lián)的精確程度、巖石物性參數(shù)的再擾動量均會影響模擬退火算法的解空間，進(jìn)而影響聯(lián)合反演目標(biāo)函數(shù)的收斂速度和聯(lián)合反演的精度.重新確定出的模擬退火的解空間越接近真實(shí)解，模擬退火算法便越容易實(shí)現(xiàn)全局尋優(yōu).本次模型試驗(yàn)中，作者以模擬退火算法的解空間需要包含真實(shí)解為基本要求，設(shè)計了模型的電阻率和速度的分布以及物性參數(shù)的再擾動量.因?yàn)閺睦碚撋戏治?，只要真?shí)解分布在再次擾動出的區(qū)間，聯(lián)合反演便可能反演出真實(shí)的物性參數(shù)分布，而不必要求先驗(yàn)的巖石物性關(guān)聯(lián)是嚴(yán)格精確的.

5 結(jié)論

通過模型試驗(yàn)，可以看出：與傳統(tǒng)的點(diǎn)對點(diǎn)映射式巖石物性關(guān)聯(lián)相比，寬范圍巖石物性約束技術(shù)對先驗(yàn)信息要求更低，容易實(shí)現(xiàn)；借助模擬退火可以實(shí)現(xiàn)空間對空間的映射關(guān)聯(lián)；并可以降低由于“不精確”巖石物性經(jīng)驗(yàn)關(guān)聯(lián)帶來的風(fēng)險，容錯性高.

總結(jié)一下該技術(shù)的應(yīng)用條件：(1)當(dāng)先驗(yàn)的巖石物性經(jīng)驗(yàn)關(guān)系可信度比較高時，選用點(diǎn)對點(diǎn)映射或?qū)挿秶鷰r石物性約束技術(shù)進(jìn)行聯(lián)合反演都是可行的；(2)當(dāng)先驗(yàn)的巖石物性經(jīng)驗(yàn)關(guān)系可信度不高時，建議選用寬范圍巖石物性約束技術(shù)以降低“不精確”信息引入帶來的風(fēng)險.(3)不存在明確的巖石物性經(jīng)驗(yàn)關(guān)系時，可選擇一個大體的巖石物性關(guān)聯(lián)，進(jìn)而套用寬范圍巖石物性約束技術(shù).

需要指出的有兩點(diǎn)，第一，新技術(shù)不局限于電阻率和速度的物性耦合約束，其適用于所有兩兩可能的巖石物性參數(shù)之間物性耦合；此外，新技術(shù)不局限于模擬退火這一種算法，其具有運(yùn)用到其他優(yōu)化算法中的可能.第二，如何在“寬范圍”內(nèi)實(shí)現(xiàn)高效便捷的巖石物性參數(shù)耦合，是該技術(shù)進(jìn)一步改進(jìn)和完善的重點(diǎn).

致謝 感謝同濟(jì)大學(xué)張羅磊副教授對本文撰寫提出的有益意見和建議.感謝兩位審稿專家的寶貴意見.

Chen L S, Wang G E. 1990. Magnetotelluric Sounding (in Chinese). Beijing: Geological Publishing House.Chen X, Yu P, Zhang L L, et al. 2010. Regularized synchronous joint inversion of magnetotelluric and seismic data.∥ SEG Technical Program Expanded Abstracts. SEG, 910-914.

Chen X, Yu P, Zhang L L, et al. 2011a. Adaptive regularized synchronous joint inversion of MT and seismic data.ChineseJ.Geophys. (in Chinese), 54(10): 2673-2681, doi: 10.3969/j.issn.0001-5733.2011.10.024.

Chen X, Yu P, Zhang L L, et al. 2011b. Joint inversion of magnetotelluric and seismic based on petrophysical property.∥ 12th China International Geo-electromagnetic induction workshop (in Chinese). Nanchang: Chinese Geophysical Society, 181-183.Chen X, Yu P, Zhang L L, et al. 2011. Applying improved very fast simulated annealing in regularized geophysical joint inversion.∥ Proceedings of the 2011 7th International Conference on Natural Computation. Shanghai, China: IEEE, 1998-2001.Chen X. 2013. Regularized joint inversion of magnetotelluric and seismic data [Ph. D. thesis] (in Chinese). Shanghai: Tongji University.

Chen X, Deng J Z, Yu P, et al. 2015. Bayesian joint inversion of MT and seismic data based on simulated annealing method.∥ International Workshop on Gravity, Electrical & Magnetic Methods and their Applications. Chengdu, China, 470-473. Colombo D, De Stefano M. 2007. Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data: Application to prestack depth imaging.TheLeadingEdge, 26(3): 326-331.Dell′Aversana P. 2001. Integration of seismic, MT and gravity data in a thrust belt interpretation.FirstBreak, 19(6): 335-341.

Gallardo L A, Meju M A. 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints.JournalofGeophysicalResearch, 109:111-121.

Gallardo L A, Fontes S L, Meju M A, et al. 2012. Robust geophysical integration through structure-coupled joint inversion and multispectral fusion of seismic reflection, magnetotelluric, magnetic, and gravity images: Example from Santos Basin, offshore Brazil.Geophysics, 77(5): B237-B251.Gao E G, Xu G M. 1996. A new kind of step by step iterative ray-tracing method.ChineseJ.Geophys. (in Chinese), 39(Suppl.): 302-308.Graul M. 1987. Seismic Lithology (in Chinese). Beijing: Petroleum Industry Press, 135-138.

Moorkamp M, Bjorn H, Jegen M, et al. 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data.Geophys.J.Int., 184(1): 477-493.

Heincke B, Jegen M, Hobbs R. 2006. Joint inversion of MT, gravity and seismic data applied to sub-basalt imaging.∥ SEG Technical Program Expanded Abstracts. SEG, 784-789.

Heincke B, Jegen M, Moorkamp M, et al. 2010. Adaptive coupling strategy for simultaneous joint inversions that use petrophysical information as constraints.∥ SEG Technical Program Expanded Abstracts 2010. SEG, 2805-2809.Heincke B, Jegen M, Moorkamp M. 2014. Joint-inversion of magnetotelluric, gravity and seismic data to image sub-basalt sediments offshore the Faroe-Islands.∥ SEG Technical Program Expanded Abstracts. SEG, 770-775.Ingber L. 1989. Very fast simulated re-annealing.MathematicalandComputerModelling, 12(8): 967-973.

Jegen M D, Hobbs R W, Tarits P, et al. 2009. Joint inversion of marine magnetotelluric and gravity data incorporating seismic constraints: Preliminary results of sub-basalt imaging off the Faroe Shelf.EarthandPlanetaryScienceLetters, 282(1-4): 47-55.

Lelièvre P G, Farquharson C G, Hurich C A. 2012. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration.Geophysics, 77(1): K1-K15.

Xu Q Z, Wu J H, Zhang L F. 2000. Rock wave velocity and resistivity in fluid flooding process under strata condition.GeophysicalProspectingforPetroleum(in Chinese), 39(1): 124-126.

Yan L J, Su Z L, Hu W B, et al. 2004. Near surface velocity imaging with the resistivity inverted from shallow transient electromagnetic sounding data.GeophysicalProspectingforPetroleum(in Chinese), 43(5): 487-491.

Yang H, Wang J L, Wu J S, et al. 2002. Constrained joint inversion of magnetotelluric and seismic data using simulated annealing algorithm.ChineseJ.Geophys. (in Chinese), 45(5): 723-734.Yu P, Dai M G, Wang J L, et al. 2009. Joint inversion of magnetotelluric and seismic data based on random resistivity and velocity distributions.ChineseJ.Geophys. (in Chinese), 52(4): 1089-1097, doi: 10.3969/j.issn.0001-5733.2009.04.026.

Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. Netherlands: Elsevier Science.

附中文參考文獻(xiàn)

陳樂壽, 王光鍔. 1990. 大地電磁測深法. 北京: 地質(zhì)出版社.

陳曉, 于鵬, 張羅磊等. 2011a. 地震與大地電磁測深數(shù)據(jù)的自適應(yīng)正則化同步聯(lián)合反演. 地球物理學(xué)報, 54(10): 2673-2681, doi: 10.3969/j.issn.0001-5733.2011.10.024.

陳曉, 于鵬, 張羅磊等. 2011b. 基于巖石物性關(guān)聯(lián)的MT與地震聯(lián)合反演.∥ 第十屆中國國際地球電磁學(xué)術(shù)討論會. 南昌: 中國地球物理學(xué)會, 181-183.

陳曉. 2013. 大地電磁與地震正則化聯(lián)合反演研究[博士論文]. 上海: 同濟(jì)大學(xué).

高爾根, 徐果明. 1996. 二維速度隨機(jī)分布逐步迭代射線追蹤方法. 地球物理學(xué)報, 39(Suppl.): 302-308.

格勞爾 M. 1987. 地震巖性學(xué). 北京: 石油工業(yè)出版社, 135-138.

徐群洲, 伍菁華, 張露菲. 2000. 地層條件下流體驅(qū)排時巖石波速和電阻率的實(shí)驗(yàn)研究. 石油物探, 39(1): 124-126.

嚴(yán)良俊, 蘇朱劉, 胡文寶等. 2004. 淺層瞬變反演電阻率資料在表層速度成像中的應(yīng)用. 石油物探, 43(5): 487-491.

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

于鵬, 戴明剛, 王家林等. 2009. 電阻率和速度隨機(jī)分布的MT與地震聯(lián)合反演. 地球物理學(xué)報, 52(4): 1089-1097, doi: 10.3969/j.issn.0001-5733.2009.04.026.

(本文編輯 何燕)

Joint inversion of MT and seismic data based on wide-range petrophysical constraints

CHEN Xiao1,2, YU Peng3*, DENG Ju-Zhi1, ZHANG Lei1, GE Kun-Peng1, WANG Yan-Guo1

1FundamentalScienceonRadioactiveGeologyandExplorationTechnologyLaboratory,EastChinaInstituteofTechnology,Nanchang330013,China2HubeiSubsurfaceMulti-scaleImagingKeyLaboratory(SMIL),ChinaUniversityofGeosciences,Wuhan430074,China3StateKeyLaboratoryofMarineGeology,TongjiUniversity,Shanghai200092,China

Petrophysical constraints play an important role in a geophysical joint inversion. A relationship of different physical properties is often used as a mapping from one physical property to another one (point to point) in a general joint inversion, which requires the empirical relationship of different physical properties is accurate. This requirement potentially reduces the flexibility of the application to the joint inversion with petrophysical constraints. In this paper, we develop a wide-range petrophysical constraints technology based on a global optimal simulated annealing (SA) algorithm. By combining model perturbation，model acceptance criterion of the SA method and petrophysical constraints, a space to space mapping is fulfilled in this method. Therefore, the physical properties of resistivity and velocity can be coupled together. And we apply this technology to a regularized synchronous joint inversion of magnetotelluric and seismic data and test the applicability of the technology based on wide-range petrophysical constraints under different priori conditions. Model tests show that the new technology is robust and easy to use. Meanwhile, the risk introduced by the incorrect relationship of physical properties is reduced, and the flexibility of the application of the petrophysical constraints is also improved.

Joint inversion; Petrophysical constraints; Magnetotelluric (MT); Seismic; Simulated annealing

10.6038/cjg20161228.

同濟(jì)大學(xué)海洋地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室開放基金(MGK1506),中國地質(zhì)大學(xué)(武漢)地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室開放基金(SMIL-2015-09)，國家自然科學(xué)基金(41604104,41504054,41504098，41664006)，江西省教育廳科學(xué)技術(shù)研究項(xiàng)目(GJJ150581)，東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室開放基金(RGET1405)，東華理工大學(xué)博士科研啟動基金(DHBK201403)和國防科工局項(xiàng)目(科工技[2013]969號)聯(lián)合資助.

陳曉，男，1986年生，東華理工大學(xué)講師，博士，主要研究方向?yàn)槎喾N地球物理數(shù)據(jù)的聯(lián)合反演.E-mail:dwjhtj@hotmail.com

*通訊作者 于鵬,男,1969年生,教授,博士生導(dǎo)師，主要從事綜合地球物理正反演方法的研究工作.E-mail:yupeng@#edu.cn

10.6038/cjg20161228

P631

2015-09-01，2016-10-03收修定稿

陳曉, 于鵬， 鄧居智等. 2016. 基于寬范圍巖石物性約束的大地電磁和地震聯(lián)合反演. 地球物理學(xué)報，59(12):4690-4700,

Chen X, Yu P, Deng J Z, et al. 2016. Joint inversion of MT and seismic data based on wide-range petrophysical constraints.ChineseJ.Geophys. (in Chinese),59(12):4690-4700,doi:10.6038/cjg20161228.