李 輝,王華忠,張 兵,3

(1.同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與反演成像研究組,上海200092;2.青鳳致遠(yuǎn)應(yīng)用地球物理研究所,上海200093;3.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

層析反演中的正則化方法研究

李 輝1,2,王華忠1,張 兵1,3

(1.同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與反演成像研究組,上海200092;2.青鳳致遠(yuǎn)應(yīng)用地球物理研究所,上海200093;3.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

正則化可顯著降低層析反演解的非唯一性，提高層析反演結(jié)果的質(zhì)量。主要研究了模型參數(shù)正則化和數(shù)據(jù)正則化。地下介質(zhì)參數(shù)之間的關(guān)聯(lián)性如何加入模型正則化是討論的問題之一；觀測數(shù)據(jù)之間的關(guān)聯(lián)性加入數(shù)據(jù)正則化的方法則是另一個(gè)主要議題。此外，討論了Tikhonov正則化和預(yù)條件兩種模型正則化實(shí)現(xiàn)策略，指出前者理論比較直觀，后者計(jì)算效率更高，并證明了兩者在理論上的等價(jià)性。模型正則化通過構(gòu)造各向異性光滑算子加入地質(zhì)構(gòu)造特征，數(shù)據(jù)正則化則通過在層析矩陣中加入預(yù)先構(gòu)造的數(shù)據(jù)預(yù)條件矩陣來實(shí)現(xiàn)。通過層析偏移速度分析給出了模型正則化和數(shù)據(jù)正則化的具體實(shí)現(xiàn)策略。理論分析和層析偏移速度分析的數(shù)值實(shí)驗(yàn)說明本文的模型正則化和數(shù)據(jù)正則化可顯著提高層析反演的質(zhì)量。

層析偏移速度分析；模型正則化；數(shù)據(jù)正則化；預(yù)條件；地質(zhì)構(gòu)造約束

隨著勘探地震技術(shù)的發(fā)展以及石油工業(yè)需求的提高,疊前深度偏移逐漸成為工業(yè)應(yīng)用中偏移技術(shù)的主流。在觀測數(shù)據(jù)品質(zhì)有保障的前提下,疊前深度偏移在實(shí)際應(yīng)用中成功與否最重要的前提是速度模型的準(zhǔn)確程度。深度域速度建模的方法包括分析類方法(如掃描速度分析、剩余曲率分析)、層析方法、波動方程類反演(如全波形反演)方法等。其中,分析類方法存在較強(qiáng)的理論假設(shè),導(dǎo)致估計(jì)的速度精度及分辨率非常低;波動方程類反演方法的效率非常低,且嚴(yán)重依賴于初始速度模型;基于射線理論的層析方法是目前工業(yè)界應(yīng)用最廣泛的深度域速度估計(jì)方法[1-3],但存在明顯的缺點(diǎn),如:只能反演光滑的背景速度、建立的反問題病態(tài)性較明顯等,所幸的是,正則化可顯著解決射線層析中的上述問題[4-5]。盡管正則化的本質(zhì)思想十分明確,但具體實(shí)現(xiàn)方法及種類繁多,本文針對模型參數(shù)之間的關(guān)聯(lián)性和數(shù)據(jù)之間的關(guān)聯(lián)性討論層析反演中加入正則化的策略。

正則化是地球物理反問題中非常重要的環(huán)節(jié),可以顯著降低層析反演解的非唯一性,提高層析反演結(jié)果的質(zhì)量。地震層析反演中的模型參數(shù)是地下介質(zhì)參數(shù)(通常是地震波速度),數(shù)據(jù)是觀測到的波場信息(通常是旅行時(shí))。地下不同空間區(qū)域的介質(zhì)參數(shù)之間有一定關(guān)聯(lián)性;同樣,不同觀測點(diǎn)的數(shù)據(jù)也相互關(guān)聯(lián)。此外,測量的地震波數(shù)據(jù)良莠不齊,可靠性存在明顯差異。如何把這些信息加入層析過程以提高反演質(zhì)量是一個(gè)值得研究的課題。

在地震層析反演中,利用地下構(gòu)造信息約束模型參數(shù)的空間分布特征,把模型參數(shù)在空間中的相關(guān)特征通過構(gòu)造信息提取出來,結(jié)合Tikhonov模型正則化或預(yù)條件思想把此信息加入層析反演中,可顯著改善估計(jì)的模型參數(shù),這里稱之為模型正則化。Zhou等[6]、Zhou[7]和Zdraveva等[8]把光滑算子加入靈敏度核函數(shù)中對反演的模型直接進(jìn)行正則化約束;Clapp[9]和Clapp等[10]利用預(yù)條件思想實(shí)現(xiàn)對反演模型的預(yù)條件約束。事實(shí)上,Tikhonov模型正則化和預(yù)條件模型正則化等價(jià),雖然兩者有不同的表達(dá)形式及理論基礎(chǔ)。類似于模型參數(shù),利用協(xié)方差矩陣可提取出觀測數(shù)據(jù)之間的相關(guān)性,將數(shù)據(jù)相關(guān)性作為數(shù)據(jù)先驗(yàn)信息加入反演過程是貝葉斯反演思想[4]的一個(gè)關(guān)鍵環(huán)節(jié)。此外,測量的數(shù)據(jù)有優(yōu)劣之分,反演中如何通過正則化突出質(zhì)量較高的數(shù)據(jù)也是本文的研究內(nèi)容。層析中加入數(shù)據(jù)關(guān)聯(lián)性和突出高質(zhì)量的數(shù)據(jù)稱為數(shù)據(jù)正則化。

1 模型正則化理論及實(shí)現(xiàn)策略

層析反演中針對模型施加正則化是把我們對模型參數(shù)的認(rèn)識加入層析反演中。地下介質(zhì)的模型參數(shù)化后,不同參數(shù)之間存在一定的關(guān)聯(lián)性,因?yàn)樗袇?shù)組成的地下介質(zhì)遵循一定的地質(zhì)規(guī)律。這里研究的模型正則化試圖在層析過程中加入地質(zhì)構(gòu)造約束,我們通過在層析矩陣中加入光滑矩陣來實(shí)現(xiàn)。光滑矩陣可約束模型中不同方向的突變成分,令垂直于反射界面與平行于反射界面的方向光滑強(qiáng)度不同來實(shí)現(xiàn)對模型的各向異性光滑約束。

模型正則化的具體實(shí)現(xiàn)策略包括Tikhonov模型正則化和預(yù)條件模型正則化兩種。無論是Tikhonov正則化還是預(yù)條件模型正則化,光滑矩陣的構(gòu)建是重要的一環(huán)。

1.1 Tikhonov模型正則化

貝葉斯框架下層析反演等價(jià)于優(yōu)化如下的目標(biāo)泛函[4]:

(1)

(2)

式中:矩陣A是A(m)的線性化算子;Δm=m-m0是當(dāng)前迭代中模型的更新量;m0是參考模型;Δd=Am0-d是當(dāng)前迭代中的數(shù)據(jù)殘差。

目標(biāo)泛函((1)式)和方程(2)中的協(xié)方差矩陣及阻尼因子實(shí)現(xiàn)了對層析反演施加正則化。不考慮正則化,僅用數(shù)據(jù)殘差進(jìn)行層析反演,即協(xié)方差矩陣為單位矩陣且ε=0時(shí),層析方程(2)退化成一般的層析方程,即:

(3)

令方程(2)中的數(shù)據(jù)協(xié)方差矩陣為單位矩陣,只考慮模型正則化時(shí)的層析方程為:

(4)

模型協(xié)方差矩陣的逆矩陣不能直接構(gòu)造,一般利用其它矩陣代替,從而實(shí)現(xiàn)對模型的正則化約束。Xu等[11]和Zhou[7]約束了模型的突變成分,試圖得到模型的光滑解,方程(4)可寫成:

(5)

式中:R是粗糙化算子,用于提取模型的突變部分。方程(5)即為層析中加入光滑正則化后的方程,其對應(yīng)的解為光滑解。事實(shí)上,方程(5)所表達(dá)的模型正則化即為Tikhonov模型正則化[12]。

模型正則化的關(guān)鍵是如何在方程(5)中引入構(gòu)造信息。提取模型突變部分可對模型向量做微分,典型的方式為拉普拉斯算子作用到模型向量上。已知模型的構(gòu)造特征時(shí),可令R垂直于地質(zhì)界面方向，提取Δm的突變信息時(shí),空間展布范圍較窄,平行于界面方向時(shí)展布范圍較寬,從而構(gòu)造出包含地質(zhì)構(gòu)造特征的各向異性粗糙化算子R。等價(jià)于約束平行于界面方向較垂直界面方向更加光滑,實(shí)現(xiàn)對模型的各向異性光滑約束。此時(shí)，方程(5)的解即為引入地質(zhì)構(gòu)造信息的模型正則化反演結(jié)果。

粗糙化算子R的逆矩陣是一個(gè)光滑算子,相對于粗糙化矩陣,光滑矩陣更容易構(gòu)造,且物理意義更為明確。所以,我們進(jìn)一步改造方程(5),利用光滑算子對層析問題施加正則化。在方程(5)兩端左乘矩陣(RTR)-1,得到:

(6)

即:

SSTATAΔm+εΔm=SSTATΔd

(7)

式中:S是R的逆矩陣,是一個(gè)光滑算子。光滑矩陣S中地質(zhì)構(gòu)造特征體現(xiàn)為沿構(gòu)造界面方向的平滑范圍較大,垂直構(gòu)造界面方向的平滑范圍較小。

1.2 預(yù)條件模型正則化

方程(3)所示的層析問題在模型預(yù)條件[11]思想下可表示為:

ASu=Δd

(8)

(9)

式中:S是預(yù)條件算子;u為層析方程的預(yù)條件解。先求解方程(8),將求解結(jié)果u代入公式(9)即可得到最終的解Δm。預(yù)條件層析方程(8)的阻尼最小二乘方程為:

STATASu+εu=STATΔd

(10)

方程(10)對應(yīng)的解是方程(8)的阻尼最小二乘解[13-14]。

下面討論模型預(yù)條件方程((8)式和(9)式)與Tikhonov模型正則化策略((5)式)之間的關(guān)系,從而探討方程(10)中矩陣S與地質(zhì)構(gòu)造正則化之間的關(guān)系。把預(yù)條件方程(9)代入方程(10),預(yù)條件層析的阻尼最小二乘方程變成:

STATAΔm+εS-1Δm=STATΔd

(11)

在方程(11)兩端左乘矩陣(ST)-1,得到:

(12)

即:

(13)

對比方程(5)和方程(13),當(dāng)方程(13)中矩陣S為方程(5)中矩陣R的逆矩陣時(shí)兩方程等價(jià)。

由上述分析可知,理論上,預(yù)條件層析問題方程(8)和方程(9)的阻尼最小二乘解等價(jià)于方程(5)所示Tikhonov模型正則化的解,前提是預(yù)條件算子S與正則化算子R互為逆矩陣。與Tikhonov模型正則化相同,當(dāng)S是光滑算子時(shí),預(yù)條件方程(8) 的最小二乘阻尼解u利用(9)式計(jì)算出的最終解Δm等價(jià)于對層析施加光滑正則化后得到的解。

對比方程(5),方程(7)和方程(8)可知,構(gòu)建方程(5)需要計(jì)算2次矩陣乘,構(gòu)建方程(7)需計(jì)算3次矩陣乘,構(gòu)造方程(8)僅需1次矩陣乘。從構(gòu)建層析方程的計(jì)算效率上看,方程(8)最優(yōu)。Clapp等[10]指出方程(8)所示的預(yù)條件正則化方式收斂性優(yōu)于方程(5)所示的模型約束的Tikhonov正則化。

1.3 光滑矩陣的構(gòu)建策略

如何構(gòu)建地質(zhì)構(gòu)造約束的光滑矩陣是模型正則化的關(guān)鍵點(diǎn)。光滑矩陣的每一行都是地下介質(zhì)空間中的一個(gè)光滑函數(shù)[6]。我們令此光滑函數(shù)為高斯光滑函數(shù),在不考慮地質(zhì)構(gòu)造時(shí),光滑矩陣中第i行第j列的元素為:

(14)

式中:xi,yi,zi是第i行高斯函數(shù)中心點(diǎn)空間坐標(biāo)的3個(gè)分量,xj,yj,zj是第j列對應(yīng)的空間位置的3個(gè)坐標(biāo)分量,σx,σy,σz是高斯函數(shù)在3個(gè)坐標(biāo)方向的標(biāo)準(zhǔn)差。將地質(zhì)構(gòu)造特征加入光滑矩陣中,(14)式表達(dá)的高斯函數(shù)在三維空間平移加旋轉(zhuǎn)至一個(gè)局部笛卡爾坐標(biāo)系中,該局部坐標(biāo)系的坐標(biāo)方向標(biāo)記為u,v,w,其中u與地質(zhì)界面的走向一致,w坐標(biāo)軸與地質(zhì)界面垂直,u,v,w組成一個(gè)右手系,稱此坐標(biāo)系為“局部地質(zhì)坐標(biāo)系”。在此局部坐標(biāo)系中,定義光滑矩陣中的高斯函數(shù)為:

(15)

式中:σu i,σv i,σw i是高斯函數(shù)在局部坐標(biāo)系中的標(biāo)準(zhǔn)差。在實(shí)際應(yīng)用中,σu i和σv i比σw i大,即在平行于地質(zhì)界面方向的光滑范圍較大,垂直于地質(zhì)界面方向的光滑范圍小。這樣,層析反演模型參數(shù)的空間分布特征被已知的地質(zhì)特征約束。

局部地質(zhì)坐標(biāo)系由原始笛卡爾坐標(biāo)系的平移和旋轉(zhuǎn)得到。原始坐標(biāo)系的原點(diǎn)位于(0,0,0)T,局部地質(zhì)坐標(biāo)系的原點(diǎn)在原始坐標(biāo)系中位于(xi,yi,zi)T,所以坐標(biāo)變換中的平移量等于(xi,yi,zi)T。平移后的坐標(biāo)系通過三維旋轉(zhuǎn)可得局部地質(zhì)坐標(biāo)系。如圖1所示,首先把坐標(biāo)系以z坐標(biāo)軸為中心旋轉(zhuǎn)角度φ,x坐標(biāo)軸旋轉(zhuǎn)至與地質(zhì)走向一致,此時(shí)x,y坐標(biāo)分別寫成x′,y′;然后以x′坐標(biāo)軸為中心旋轉(zhuǎn)角度θ,z坐標(biāo)軸旋轉(zhuǎn)至垂直于地質(zhì)界面的方向。旋轉(zhuǎn)后的坐標(biāo)系即為局部地質(zhì)坐標(biāo)系。旋轉(zhuǎn)角度的定義為面對中心坐標(biāo)軸觀察時(shí)沿逆時(shí)針方向旋轉(zhuǎn)的角度,如圖1所示。以z軸與x軸為中心旋轉(zhuǎn)的旋轉(zhuǎn)矩陣分別為:

(16)

(17)

局部地質(zhì)坐標(biāo)系中的坐標(biāo)(u,v,w)T與原始坐標(biāo)系中坐標(biāo)(x,y,z)T的關(guān)系為:

(18)

圖1 三維坐標(biāo)旋轉(zhuǎn)示意圖解

從上述過程可知,構(gòu)建含構(gòu)造特征的光滑矩陣時(shí)需要地質(zhì)構(gòu)造的傾角信息以及高斯光滑函數(shù)的標(biāo)準(zhǔn)差σu i,σv i,σw i。數(shù)值實(shí)驗(yàn)部分將結(jié)合具體數(shù)據(jù)給出傾角信息和標(biāo)準(zhǔn)差的提取過程。

2 數(shù)據(jù)正則化理論及實(shí)現(xiàn)策略

層析反演中針對數(shù)據(jù)施加正則化是把我們對數(shù)據(jù)的認(rèn)識加入層析反演中。層析反演中正則化關(guān)注的數(shù)據(jù)特點(diǎn)主要有兩種:①類似于不同模型參數(shù)間的關(guān)系,不同數(shù)據(jù)之間也存在一定的關(guān)聯(lián)性,例如單炮數(shù)據(jù)中同一波前的信號特征(如旅行時(shí)、振幅、相位)在不同檢波點(diǎn)的變化應(yīng)連續(xù);②不同觀測數(shù)據(jù)在測量過程中有不同的準(zhǔn)確性,即不同數(shù)據(jù)的可信程度不同,例如層析反演中檢測波前到達(dá)時(shí),不同數(shù)據(jù)的檢測精度不同。如何在層析過程中考慮數(shù)據(jù)的特點(diǎn)是數(shù)據(jù)正則化的主要任務(wù)。我們在基于角度域共成像點(diǎn)道集的層析偏移速度分析中考慮數(shù)據(jù)正則化問題。

數(shù)據(jù)之間的關(guān)聯(lián)性同樣可以通過對數(shù)據(jù)進(jìn)行光滑加入層析反演中。此時(shí)需對相關(guān)聯(lián)的數(shù)據(jù)平滑,消除它們之間由人為因素引入的突變。層析方程(3)在考慮數(shù)據(jù)關(guān)聯(lián)性后可表達(dá)為:

AΔm=CΔd

(19)

其中,矩陣C對數(shù)據(jù)向量Δd改造,使得相關(guān)聯(lián)的數(shù)據(jù)之間變化平緩,符合物理規(guī)律。層析偏移速度分析在成像域(偏移剖面和共成像點(diǎn)道集)測量數(shù)據(jù),實(shí)現(xiàn)反射波旅行時(shí)層析,具體實(shí)現(xiàn)方法不再贅述,可參考文獻(xiàn)[1]和文獻(xiàn)[15]。層析偏移速度分析中考慮數(shù)據(jù)關(guān)聯(lián)性主要包括兩方面:①反射點(diǎn)鄰近、反射角相同的射線對應(yīng)的數(shù)據(jù)殘差不應(yīng)存在突變(圖2a);②同一反射點(diǎn)反射角鄰近的反射射線對應(yīng)的數(shù)據(jù)殘差不應(yīng)存在突變(圖2b)。矩陣C修改滿足上述兩個(gè)條件的數(shù)據(jù),使數(shù)據(jù)間滿足上述物理規(guī)律。

圖2 存在關(guān)聯(lián)性的反射射線

不同于其它層析方法,層析偏移速度分析中的數(shù)據(jù)測量過程比較復(fù)雜[15],首先在偏移剖面中掃描反射面傾角,同時(shí)在成像道集中拾取不同反射角度的反射點(diǎn)成像深度,最后把不同角度之間的成像深度差轉(zhuǎn)換成反射射線的旅行時(shí)殘差。復(fù)雜的數(shù)據(jù)測量過程可能導(dǎo)致測量數(shù)據(jù)的誤差較大,例如反射面傾角誤差、不同反射角成像深度誤差等都將轉(zhuǎn)移至最終的數(shù)據(jù)誤差中,所以,在層析偏移速度分析中考慮測量數(shù)據(jù)的準(zhǔn)確性尤為重要。我們對此問題的解決方案是壓制精度較低數(shù)據(jù)同時(shí)提升精度較高數(shù)據(jù)在層析反演過程中所起的作用,具體策略是在層析方程(3)中加入權(quán)系數(shù)矩陣,即:

WAΔm=WΔd

(20)

圖3 偏移剖面中的傾角與相似系數(shù)示意

圖4 角度域共成像點(diǎn)道集中的拾取曲線與相似系數(shù)示意

式中:W是對角矩陣,其對角元素是數(shù)據(jù)對應(yīng)的權(quán)系數(shù),數(shù)據(jù)精度越高系數(shù)越大,數(shù)據(jù)精度越低系數(shù)越小。關(guān)于數(shù)據(jù)質(zhì)量的正則化的關(guān)鍵點(diǎn)是如何考察不同數(shù)據(jù)的測量精度。層析偏移速度分析中數(shù)據(jù)測量精度取決于在偏移剖面上掃描反射面傾角和在共成像點(diǎn)道集中拾取成像深度的精度,所以只需考察反射面傾角和成像深度的精度即可。反射面傾角準(zhǔn)確時(shí),偏移剖面上反射同相軸沿傾角方向連續(xù)性較好,否則同相軸不連續(xù)。如圖3所示,在偏移剖面中通過沿傾角方向信號的相似系數(shù)考察所測傾角的精度,其中紅色實(shí)線是正確反射面傾角方向,虛線是不正確的傾角參考方向,傾角方向準(zhǔn)確時(shí)計(jì)算的相似系數(shù)接近1,其它方向的相似系數(shù)相對非常小。圖3中定義豎直向上方向的角度為0,逆時(shí)針方向旋轉(zhuǎn)為正角度。同樣,在共成像點(diǎn)道集中拾取不同角度成像深度的精度也用相似系數(shù)考察。圖4a顯示了在角度域共成像點(diǎn)道集中自動拾取的成像深度,圖4b是不同層位拾取結(jié)果的相似系數(shù)對比,圖中5個(gè)紅色圓點(diǎn)是圖4a中5個(gè)同相軸拾取位置的相似系數(shù),可見,拾取質(zhì)量比較好的同相軸對應(yīng)的相似系數(shù)接近1,拾取質(zhì)量較差的同相軸對應(yīng)的相似系數(shù)較小。最終,方程(20)中矩陣W中的權(quán)系數(shù)可通過上述兩種相似系數(shù)來確定:

(21)

綜合模型正則化和數(shù)據(jù)正則化,可得最終的層析方程。

Tikhonov模型正則化策略((7)式)結(jié)合數(shù)據(jù)正則化策略((19)式和(20)式)可得層析方程為:

SSTATWTWAΔm+εΔm=SSTATWTCΔd

(22)

模型預(yù)條件正則化策略((8)式和(9)式)結(jié)合數(shù)據(jù)正則化策略((19)式與(20)式)可得層析方程為:

WASu=WCΔd

(23)

Δm=Su

(24)

上文分析了層析問題關(guān)于模型正則化與數(shù)據(jù)正則化的理論方法。模型正則化包括Tikhonov模型正則化和預(yù)條件模型正則化,前者理論直觀,后者計(jì)算效率較高。理論分析結(jié)果表明,Tikhonov模型正則化和預(yù)條件模型正則化在理論上等價(jià)。數(shù)據(jù)正則化考慮了數(shù)據(jù)參數(shù)之間的關(guān)聯(lián)性以及不同的數(shù)據(jù)質(zhì)量。結(jié)合模型正則化和數(shù)據(jù)正則化,可得最終的正則化方程。

3 數(shù)值實(shí)驗(yàn)

本文理論部分地質(zhì)構(gòu)造約束的模型正則化與數(shù)據(jù)正則化可應(yīng)用在不同種類的層析成像中,例如初至波層析、井間層析、傾斜層析以及層析偏移速度分析等。層析偏移速度分析[1,16]是一種經(jīng)典的深度域速度建模方法,其缺點(diǎn)是射線角度覆蓋范圍較小,所以相應(yīng)的層析問題病態(tài)性嚴(yán)重。這里將地質(zhì)約束的正則化策略應(yīng)用于層析偏移速度分析,以提高層析的穩(wěn)定性,同時(shí)改善反演的質(zhì)量。相對于偏移距域共成像點(diǎn)道集(ODCIG),角度域共成像點(diǎn)道集(ADCIG)更加符合物理規(guī)律,并且在ADCIG中實(shí)現(xiàn)層析偏移速度分析時(shí)層析核函數(shù)的構(gòu)建只需進(jìn)行初值射線追蹤,而ODCIG中的層析偏移速度分析構(gòu)建核函數(shù)時(shí)必須做耗時(shí)且不穩(wěn)定的兩點(diǎn)射線追蹤,所以,本文在角度域共成像點(diǎn)道集中實(shí)現(xiàn)層析偏移速度分析。

層析偏移速度分析的實(shí)現(xiàn)流程如圖5所示?；谝阎某跏妓俣饶Ｐ瓦M(jìn)行疊前深度偏移并生成ADCIG,以道集中同相軸是否拉平為準(zhǔn)則考察速度場的優(yōu)劣。ADCIG中同相軸不平時(shí)在偏移剖面中選取反射點(diǎn)用于反射層析,并計(jì)算反射點(diǎn)處的構(gòu)造傾角,在ADCIG中拾取不同角度的成像深度,同時(shí)把不同角度的成像深度轉(zhuǎn)化為旅行時(shí)差,即為層析反演中的數(shù)據(jù)殘差。根據(jù)反射點(diǎn)的地下位置、反射角度和偏移剖面中的界面傾角,利用初值射線追蹤模擬出反射射線,建立無約束的層析矩陣,即層析方程(3)中的矩陣A。同時(shí),利用偏移剖面和ADCIG建立正則化矩陣方程。求解層析方程更新初始速度,得到更新后的速度場。如此,則實(shí)現(xiàn)了層析偏移速度分析的一次迭代。

圖5 層析偏移速度分析流程

上文討論了兩種構(gòu)造約束的模型正則化策略,這里,通過預(yù)條件正則化測試本文地質(zhì)構(gòu)造約束的正則化方法。此外,我們也對兩種數(shù)據(jù)正則化方式進(jìn)行測試,以考察數(shù)據(jù)正則化對層析的影響。

利用二維理論模型和實(shí)際數(shù)據(jù)的二維測線對本文正則化方法進(jìn)行測試。其中,理論模型測試分為兩部分:①模型正則化測試,詳細(xì)表述模型正則化中構(gòu)造約束的光滑矩陣如何構(gòu)建;②數(shù)據(jù)正則化測試,給出數(shù)據(jù)正則化影響層析結(jié)果的直觀認(rèn)識。

3.1 理論模型測試

理論模型參數(shù)及觀測系統(tǒng)如表1所示。理論模型(圖6a)中有4個(gè)沉積層,其中第3層內(nèi)有一個(gè)高速體，其偏移結(jié)果如圖6b所示。初始速度模型和相應(yīng)的偏移結(jié)果如圖7 所示。從偏移結(jié)果中提取地質(zhì)構(gòu)造傾角和高斯光滑函數(shù)的標(biāo)準(zhǔn)差是構(gòu)建光滑矩陣的關(guān)鍵點(diǎn)。下面將從初始速度偏移剖面出發(fā)獲取地質(zhì)構(gòu)造傾角和高斯函數(shù)標(biāo)準(zhǔn)差。模型的構(gòu)造信息可從圖7b所示的偏移剖面中提取,具體策略與偏移速度分析中的反射點(diǎn)傾角提取策略[15]相同,從圖7a中得到的反射界面傾角如圖8所示,其中無構(gòu)造區(qū)域的傾角置為0。利用偏移剖面和反射界面傾角計(jì)算出高斯光滑函數(shù)在不同空間點(diǎn)的標(biāo)準(zhǔn)差σu和σw,計(jì)算公式為:

(25)

(26)

式中:s是偏移剖面上沿構(gòu)造方向的相似系數(shù);sthresh是相似系數(shù)門檻值,滿足s

表1 理論模型參數(shù)與觀測系統(tǒng)參數(shù)

圖6 理論速度模型及相應(yīng)偏移結(jié)果

圖7 初始速度模型及相應(yīng)偏移結(jié)果

為測試數(shù)據(jù)正則化對層析結(jié)果的影響,人為修改層析偏移速度分析中的數(shù)據(jù)使之污染,引入部分錯(cuò)誤數(shù)據(jù)。針對此被污染的數(shù)據(jù),分別利用無數(shù)據(jù)正則化(圖12)的層析方法和加入數(shù)據(jù)正則化(圖13) 的層析方法更新圖7a所示的初始速度模型,更新后的速度模型分別如圖12a和圖13a所示。對比圖12a和圖13a 可以看出,圓圈區(qū)域由于數(shù)據(jù)問題,無數(shù)據(jù)正則化的層析速度明顯偏大,加入數(shù)據(jù)正則化后的層析速度偏小;但就速度更新的特點(diǎn)來看,兩者的速度更新方向均正確,前者速度更新量過大,后者速度更新量偏小。前者影響層析的收斂性,可能導(dǎo)致反演結(jié)果不收斂;后者影響層析的收斂速度,導(dǎo)致反演收斂慢。從上述分析可知數(shù)據(jù)正則化的作用是以犧牲收斂效率為代價(jià)保證正確收斂。相應(yīng)的偏移剖面如圖12b和13b所示,ADCIG如圖12c,圖12d,圖13c和圖13d所示。

圖8 初始速度偏移剖面中提取的界面傾角

圖9 不同空間位置的標(biāo)準(zhǔn)差

圖10 構(gòu)造約束的反演速度模型及相應(yīng)的偏移結(jié)果

圖11 無構(gòu)造約束反演的速度模型

3.2 實(shí)際數(shù)據(jù)測試

用于測試高斯束層析的數(shù)據(jù)來源于中國某探區(qū)的一條二維測線。該測線共232炮,每炮最大偏移距3200m,最小偏移距100m,左側(cè)單邊接收。道間距50m,炮間距大部分為100m,小部分為50m。模型共952個(gè)CDP點(diǎn),CDP間隔25m,深度采樣點(diǎn)數(shù)為501,深度采樣間隔為10m。初始速度模型通過速度分析獲得。初始速度模型和相應(yīng)偏移剖面如圖14a和圖14b所示。在橫向位置為2500,15000,22500m的3個(gè)CDP點(diǎn)抽取出3個(gè)角度域成像道集如圖15a,圖15b和圖15c所示。從成像道集中可以看出,x=2500m附近初始速度比較準(zhǔn)確,同相軸基本水平;x=15000m附近同相軸有微小的下拉現(xiàn)象,速度略微偏大;x=22500m附近的同相軸下拉更加明顯,速度偏大。在此初始速度模型的基礎(chǔ)上,利用本文模型正則化和數(shù)據(jù)正則化約束的角度域?qū)游銎扑俣确治鼋?jīng)過兩次迭代得到的速度模型如圖16a所示,其中的構(gòu)造特征不是很明顯,這是因?yàn)槌跏妓俣饶Ｐ捅容^接近真實(shí)速度的低波數(shù)成分,從而速度更新量比較小。圖16b是層析反演的速度模型與初始速度模型之差,從圖16b中可以看出,速度模型的更新量體現(xiàn)了地下介質(zhì)的地質(zhì)構(gòu)造特征。利用圖16a 中的反演速度模型偏移出的成像剖面如圖16c所示,從圖16c中可以看出,模型中部小斷塊區(qū)域成像更清楚,同時(shí)右側(cè)的反射同相軸更加連續(xù),相對于原始剖面成像效果有明顯改善。相應(yīng)的成像道集如圖17a,圖17b和圖17c所示,相比于圖15a,圖15b和圖15c可以看出,原始成像道集中同相軸的下拉現(xiàn)象已經(jīng)被完全校正,成像效果更好。

圖12 無數(shù)據(jù)約束的反演速度模型及相應(yīng)的偏移結(jié)果

圖13 數(shù)據(jù)約束的反演速度模型及相應(yīng)的偏移結(jié)果

圖14 初始速度模型及相應(yīng)的成像剖面

圖15 初始速度模型偏移出的ADCIG

圖16 反演結(jié)果及相應(yīng)的成像剖面

圖17 反演速度模型偏移出的ADCIG

4 結(jié)論與討論

正則化是層析反演技術(shù)的主要發(fā)展方向之一,本文討論了關(guān)于模型和數(shù)據(jù)的正則化方法及其在層析偏移速度分析中的應(yīng)用。

模型正則化考慮了不同模型參數(shù)之間的關(guān)聯(lián)性。我們討論了模型正則化的兩種具體實(shí)現(xiàn)方式——Tikhonov模型正則化與預(yù)條件模型正則化。Tikhonov模型正則化從貝葉斯反演理論出發(fā),導(dǎo)出在層析過程中利用光滑矩陣實(shí)現(xiàn)模型正則化的公式。預(yù)條件層析問題的阻尼最小二乘解轉(zhuǎn)化成原模型參數(shù)實(shí)現(xiàn)了預(yù)條件模型正則化。分析了模型正則化的兩種實(shí)現(xiàn)策略之間的關(guān)系,證明了在一定條件下兩者的等價(jià)性。但預(yù)條件正則化的實(shí)現(xiàn)效率明顯比Tikhonov正則化高。

在上述模型正則化的框架下,把地質(zhì)構(gòu)造特征加入層析反演中,其中包含地質(zhì)構(gòu)造特征的光滑矩陣是關(guān)鍵點(diǎn)。模型正則化中的光滑矩陣可利用各向異性光滑算子構(gòu)建,詳細(xì)討論了如何構(gòu)建融合地質(zhì)構(gòu)造特征的各向異性光滑算子。

類似于模型參數(shù),不同的觀測數(shù)據(jù)之間也存在關(guān)聯(lián)性。討論了反射數(shù)據(jù)之間的兩種關(guān)聯(lián)特征——反射點(diǎn)鄰近且反射角相同的射線對應(yīng)的數(shù)據(jù)之間差別較小,反射點(diǎn)相同且反射角鄰近的射線對應(yīng)的數(shù)據(jù)之間差別也較小。此外,不同測量數(shù)據(jù)的測量精度有高低之分,尤其在層析偏移速度分析中由于數(shù)據(jù)的測量過程比較復(fù)雜,人為引入觀測誤差的可能性更高。

利用一個(gè)二維理論模型對模型正則化和數(shù)據(jù)正則化的具體實(shí)現(xiàn)做了詳細(xì)討論,主要包括模型正則化中光滑算子的標(biāo)準(zhǔn)差的計(jì)算策略,數(shù)據(jù)正則化中權(quán)系數(shù)的計(jì)算方式。加入正則化后層析反演得到的速度模型表現(xiàn)出了地質(zhì)構(gòu)造特征,偏移剖面和成像道集的質(zhì)量均可與準(zhǔn)確速度的偏移結(jié)果媲美。二維測線實(shí)際數(shù)據(jù)的應(yīng)用也得到了滿意的速度模型以及相應(yīng)的偏移結(jié)果。

今后的工作方向之一是在精細(xì)的偏移結(jié)果中提取出反射界面和斷層,結(jié)合測井?dāng)?shù)據(jù)作為先驗(yàn)信息,將正則化有機(jī)地融入層析反演中,以期得到更高波數(shù)的速度成分,進(jìn)一步提高層析速度反演的精度和分辨率。

致謝:感謝中國石油化工集團(tuán)公司科技部及中國石油化工股份有限公司石油物探技術(shù)研究院對WPI研究組及本項(xiàng)研究的支持。

[1] Woodward M J,Nichols D,Zdraveva O,et al.A decade of tomography[J].Geophysics,2008,73(5):VE5-VE11

[2] Hardy P.Ongong R&D in ray based tomography:well worth the effort[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,4806-4810

[3] Lambare G,Guillaume P,Montel J P.Recent advances in ray-based tomography[J].Expanded Abstracts of 76thEAGE Annual Conference,2014,We-G103-01

[4] Tarantola A.Inverse problem theory and methods for model parameter estimation[M].Philadelphia:Society for Industrial and Applied Mathematics,2004:1-352

[5] Zhdanov M S.Geophysical inverse theory and regularization problems[M].Amsterdam:Elsevier Science,2002:1-628

[6] Zhou C,Brandsberg-Dahl S,Jiao J.A continuation approach to regularize the reflection tomography with a 3D Gaussian filter[J].Expanded Abstracts of 71stEAGE Annual Conference,2009,U031

[7] Zhou C.Incorporating geologic information into reflection tomography with a dip oriented Gaussian filter[J].Expanded Abstracts of 75thEAGE Annual Conference,2013,Th-04-02

[8] Zdraveva O,Hydal S,Woodward M.Tomography with geological constraint-an alternative solution for resolving of carbonates[J].Expanded Abstracts of 75thEAGE Annual Conference,2013,Th-04-04

[9] Clapp R G.Geologically constrained migration velocity analysis[D].San Francisco:Stanford University,2001

[10] Clapp R G,Biondi B,Claernout J F.Incorporating geologic information into reflection tomography[J].Geophysics,2004,69(2):533-546

[11] Xu S,Zhang Y,Huang T.Enhanced tomography resolution by a fat ray technique[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg, 2006,3354-3358

[12] Tikhonov A N,Arsenin V Y.Solution of ill-posed problems[M].Washington:Winston & Sons,1977:1-258

[13] Backus G E,Gilbert F.The resolving power of gross earth data[J].Geophysical Journal of the Royal Astronomical Society,1968,16(1):169-205

[14] 楊文采.地球物理反演的理論與方法[M].北京:地質(zhì)出版社,1997:1-271 Yang W C.Theory and methods of geophysical inversion[M].Beijing:Geological Publishing House,1997:1-271

[15] 張兵.基于深度域共成像點(diǎn)道集的層析速度反演與建模方法研究[D].上海:同濟(jì)大學(xué),2011 Zhang B.The study of tomography and velocity model building based on common imaging gathers[D].Shanghai:Tongji University,2011

[16] Stork C.Reflection tomography in the postmigrated domain[J].Geophysics,1992,57(5):680-692

(編輯：顧石慶)

The study of regularization in tomography

Li Hui1,Wang Huazhong1,Zhang Bing1,2

(1.WavePhenomenaandInversionImagingGroup(WPI),TongjiUniversity,Shanghai200092,China; 2.QingfengzhiyuanAppliedGeophysicsInstitute,Shanghai200093,China;3.SinopecGeophysicalResearchInstitute,Nanjing211103,China)

Regularization in tomography is able to weaken the non-uniqueness of tomography to improve the inversion resultThe discussion of regularization in this paper includes model-regularization and data-regularizationModel parameters are not isolated,how to add the relationship of these parameters into tomography is one of the missions hereSimilarly,considering datum relationship in tomography is another problemThe so-called “straightforward regularization” and the “precondition regularization” are focused,and we achieve that the former is intuitionistic and the latter is more efficiencyAlso,we point out that the above two algorithms are equivalent to each other,and this will be shown in this paperThe geological structure characteristics of the medium can be integrated into the tomography using the model-regularization with anisotropic smooth matrix.The data-regularization is realized with another smooth operator which will be integrated into the tomographic matrix.The model-regularization and data-regularization are tested with tomographic migration velocity analysis (MVA) algorithm.The results of theory and numerical experiments with tomographic MVA show that the proposed model-regularization and the data-regularization are both able to improve the quality of tomography obviously.

tomographic MVA,model-regularization,data-regularization,precondition,geological structure constraint

2014-11-24;改回日期：2015-02-26。

李輝(1985—),男,博士,現(xiàn)從事射線類偏移與反演的研究工作。

國家自然科學(xué)基金(41374117)、國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(2011CB201002)、國家科技重大專項(xiàng)項(xiàng)目(2011ZX05003-003,2011ZX05005-005-008HZ,2011ZX05006-002)和中國石化地球物理重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(33550006-14-FW2099-0026)共同資助。

P631

A

1000-1441(2015)05-0569-13

10.3969/j.issn.1000-1441.2015.05.010