鄧勇， 李芳， 彭海龍， 任婷， 劉國(guó)昌

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 102249 2 中海石油(中國(guó))有限公司湛江分公司， 湛江 524057 3 中海石油(中國(guó))有限公司海南分公司，?？?570000

0 引言

隨著勘探開發(fā)力度的逐步加大，準(zhǔn)確獲取地下地質(zhì)構(gòu)造形態(tài)及展布特征是當(dāng)前油氣藏開采工作的重要基礎(chǔ).尤其在復(fù)雜斷塊區(qū)域，真實(shí)的地震成像是推動(dòng)區(qū)域油氣藏勘探開發(fā)的重要保障(張昌榕等，2018).然而，常規(guī)的成像方法難以獲取該區(qū)域準(zhǔn)確的成像，使復(fù)雜斷塊區(qū)域的勘探開發(fā)陷入瓶頸.為更好的應(yīng)對(duì)上述挑戰(zhàn)，目前采用疊前深度偏移(Prestack Depth Migration，PSDM)方法解決上述問題(秦寧，2020).

在疊前深度偏移方法中，速度模型的精度制約著偏移成像的準(zhǔn)確度，常用的速度建模方法包括速度分析法、層析反演方法、波動(dòng)方程類方法等(王艷香等，2021).分析類方法存在條件假設(shè)，速度分析結(jié)果精度低，一般只用于初始速度建模(劉定進(jìn)等，2019).波動(dòng)方程方法的反演結(jié)果依賴于初始速度模型的精度，且需要多次迭代才能完成速度反演，效率較低(曹呈浩等，2020).而層析反演方法應(yīng)用過程簡(jiǎn)單方便，具有較好的反演精度和反演效率，是當(dāng)前應(yīng)用最為廣泛的速度建模方法(徐嘉亮等，2021).

傳統(tǒng)的射線類走時(shí)層析反演方法一般采用規(guī)則網(wǎng)格單元的模型參數(shù)化方式表征空間坐標(biāo)問題，在復(fù)雜地質(zhì)構(gòu)造速度場(chǎng)反演過程中，需要對(duì)網(wǎng)格進(jìn)行精細(xì)剖分(郭振波等，2019)，但是在不規(guī)則阻抗界面區(qū)域中，該方法容易出現(xiàn)誤差，給反演結(jié)果帶來假象(徐嘉亮等，2020).目前常規(guī)解決策略是采用不規(guī)則網(wǎng)格單元進(jìn)行模型參數(shù)化，可以較好地表征不規(guī)則異常體的波阻抗界面(何雷宇等，2017).但上述方案并沒有完全解決斷層引起的速度差異性問題，限制了斷層區(qū)域速度建模精度的提升.除此之外，應(yīng)用模型正則化方法，將地質(zhì)構(gòu)造特征引入層析進(jìn)行約束，能夠提升層析質(zhì)量，但是在復(fù)雜斷塊區(qū)域中，該方法由于缺少斷層的約束控制，反演結(jié)果中缺少高波數(shù)的速度成分，因此該方法有待進(jìn)一步提升速度反演精度(李輝等，2015).

為解決斷層區(qū)域的速度精度問題，Birdus(2007)以斷層作為條件，在密集的疊前深度偏移道集網(wǎng)格上進(jìn)行非雙曲剩余曲率分析，然后約束速度更新過程.該方法需要高密度道集，且只能在與斷層對(duì)應(yīng)的有限區(qū)域內(nèi)進(jìn)行約束控制.Rodriguez等(2011)將斷層上下盤分為兩個(gè)區(qū)域分別進(jìn)行速度反演，但該方法僅更新斷層附近的速度，且在更新過程中難以實(shí)現(xiàn)有效迭代.Guo等(2015)提出用斷層附近道集的剩余曲率來約束層析正則化過程，能夠更新整個(gè)模型.但是準(zhǔn)確拾取斷層附近反射信息的剩余曲率較為困難，因此，該方法在應(yīng)用中存在一定的挑戰(zhàn).

為更好地改進(jìn)層析成像在復(fù)雜斷塊區(qū)域速度建模的精度問題，本文在模型正則化層析的基礎(chǔ)上，提出一種基于斷層控制的層析方法，該方法對(duì)成像過程中的斷層進(jìn)行識(shí)別判斷，利用斷層數(shù)據(jù)構(gòu)建預(yù)條件矩陣，然后采用TV正則化方法改進(jìn)目標(biāo)泛函，從而避免層析過程中斷層兩側(cè)速度的平滑過渡，提升斷層兩側(cè)的速度場(chǎng)精度，得到分辨率較高的斷塊區(qū)域速度.經(jīng)實(shí)際鉆井?dāng)?shù)據(jù)驗(yàn)證，該方法得到的速度場(chǎng)和實(shí)際速度場(chǎng)匹配度高，具有較好的應(yīng)用效果.

1 方法理論

在地震反演過程中，多個(gè)不同的解都適合作為同一個(gè)地球物理問題的解，使得反演問題具有不確定性.為解決上述問題，層析反演一般采用先驗(yàn)約束降低多解性，常用的方法是正則化(張力起等，2019).正則化是將逆問題的先驗(yàn)信息和逆解估計(jì)的期望加入到反演過程中，從而消除反演過程的病態(tài)問題，提高反問題求解的魯棒性和適定性(代榮獲等，2019).地球物理反問題中，通過增加正則約束來解決病態(tài)反演問題取得了很好的效果.盡管無法獲得描述真實(shí)地下介質(zhì)的完全精確解，但可以得到滿足應(yīng)用需求的近似最大似然解(蘇揚(yáng)等，2021).

正則化技術(shù)是獲得有地質(zhì)意義解的有效方法，按照應(yīng)用條件一般可分為模型正則化和數(shù)據(jù)正則化.模型正則化利用地下構(gòu)造信息約束模型參數(shù)，實(shí)現(xiàn)改善估計(jì)模型參數(shù)的目的(張兵和王華忠，2019)；數(shù)據(jù)正則化是在層析中加入空間介質(zhì)參數(shù)的數(shù)據(jù)關(guān)聯(lián)性(李輝等，2015)，但是收斂效率較低，收斂效果較差，目前應(yīng)用較少.按照正則化函數(shù)的不同，正則化技術(shù)可以分為Tikhonov正則化和總變差(TV，Total Variation)正則化.相比較而言，Tikhonov正則化計(jì)算過程穩(wěn)定，但是難以保持層析過程中的速度邊界，TV正則化能夠較好的保持速度的邊緣信息，但是穩(wěn)定性較差(王永奎等，2022).

在反演過程中引入模型的相關(guān)認(rèn)識(shí)，即構(gòu)成模型正則化.模型正則化的實(shí)現(xiàn)策略主要包括Tikhonov模型正則化和預(yù)條件模型正則化，預(yù)條件模型正則化方式收斂性又優(yōu)于模型約束的Tikhonov模型正則化，且計(jì)算效率較高，因此本文采用預(yù)條件正則化進(jìn)行層析矩陣構(gòu)建，得到精度相對(duì)較高的初始速度.然后以此為基礎(chǔ)，對(duì)成像數(shù)據(jù)中的斷層進(jìn)行識(shí)別判斷，通過識(shí)別的斷層數(shù)據(jù)建立預(yù)條件矩陣并采用TV正則方法改進(jìn)目標(biāo)泛函，實(shí)現(xiàn)速度場(chǎng)的更新，得到復(fù)雜斷塊區(qū)域的高精度速度場(chǎng).

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

層析反演采用匹配射線走時(shí)的方法來得到速度值，應(yīng)用射線追蹤方法擬合地震數(shù)據(jù)中的反射信息走時(shí)，并通過走時(shí)差的反投影得到速度的更新量(范白濤等，2022).對(duì)于實(shí)際地震數(shù)據(jù)而言，拾取的地震走時(shí)可以表示為：

(1)

式中，τobs表示拾取的走時(shí)，s表示觀測(cè)數(shù)據(jù)慢度，dl表示沿著射線路徑ab的射線長(zhǎng)度.如果已知背景慢度，則計(jì)算的走時(shí)可以表示為：

(2)

式中，τcal表示計(jì)算的走時(shí)，s0表示背景慢度，則走時(shí)擾動(dòng)和慢度擾動(dòng)具有如下關(guān)系：

(3)

將式(3)轉(zhuǎn)為矩陣，則基于射線的線性走時(shí)層析反演表達(dá)式為：

Lδs=δτ，

(4)

式中，L表示線性化算子.由(4)式可知，每一條射線都存在一個(gè)對(duì)應(yīng)的公式，所有射線所代表的公式則構(gòu)成一個(gè)射線層析方程組，因此基于貝葉斯反演理論的最大后驗(yàn)概率密度思想，層析反演的目標(biāo)泛函為：

(5)

層析反演在迭代過程中，采用線性化方式對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，則線性迭代方程為：

(6)

式中，矩陣A表示A(m)的線性化算子，Δm表示當(dāng)前迭代次數(shù)下的模型更新量，Δd表示當(dāng)前迭代中的數(shù)據(jù)殘差.

通過公式(5)和公式(6)中的數(shù)據(jù)協(xié)方差矩陣和正則化參數(shù)，能夠?qū)游龇囱葸^程實(shí)現(xiàn)正則化.在不考慮正則化的情況下，即協(xié)方差矩陣為單位矩陣且阻尼因子為0，意味著僅用數(shù)據(jù)殘差進(jìn)行層析反演，則公式(6)退化為一般的層析方程：

Δd=AΔm.

(7)

ASu=Δd，

(8)

Δm=Su，

(9)

式中：S表示預(yù)條件算子，u為預(yù)條件的解.考慮阻尼因子ε，則層析方程的阻尼最小二乘方程可以表示為：

STATASu+εu=STATΔd.

(10)

將預(yù)條件下的解u=S-1A-1Δd代入方程(10)并適當(dāng)變換得到：

SSTATAΔm+εΔm=SSTATΔd，

(11)

對(duì)于公式(11)而言，如果將地質(zhì)構(gòu)造信息引入預(yù)條件算子S，則公式(11)就是構(gòu)造約束正則化的層析方程.所以，提取地質(zhì)構(gòu)造信息，并將其引入預(yù)條件算子所表示的光滑矩陣是上述步驟的關(guān)鍵點(diǎn).

在同一個(gè)空間坐標(biāo)系下，不同的散射點(diǎn)的位置和傾角信息具有較大的差異性，為構(gòu)建一個(gè)空間上的光滑矩陣，可以通過坐標(biāo)變換的方式將上述地質(zhì)信息融入到變換后的光滑矩陣中(萬弘等，2017)，空間坐標(biāo)的轉(zhuǎn)換示意圖如圖1所示.

圖1 坐標(biāo)變換示意圖Fig.1 Schematic diagram of coordinate transformation

圖1中，逆時(shí)針方向表示旋轉(zhuǎn)的正方向，局部地質(zhì)坐標(biāo)系原點(diǎn)為散射點(diǎn)中心，記為T(u,v,w)，u軸的方向就是反射界面的走向，w軸方向是反射界面垂直方向；空間坐標(biāo)系記為T′(x,y,z).然后參照李輝等人提出的坐標(biāo)轉(zhuǎn)換方法，空間坐標(biāo)經(jīng)過一系列的旋轉(zhuǎn)、平移后變換為局部地質(zhì)坐標(biāo)系(李輝等，2015).記散射點(diǎn)的原始空間坐標(biāo)為T′(x0,y0,z0)，則旋轉(zhuǎn)后的局部地質(zhì)坐標(biāo)可以表示T(u0,v0,w0)，兩個(gè)坐標(biāo)系存在以下關(guān)系：

采訪正值“第十九屆國(guó)產(chǎn)高濃度磷復(fù)肥產(chǎn)銷會(huì)”期間，楊三可來到甕福集團(tuán)展臺(tái)，向國(guó)內(nèi)外客商推介甕福集團(tuán)“DAP+”、磷酸二氫鉀、聚磷酸銨等新型增效肥料。

(u0,v0,w0)=T(x-x0,y-y0,z-z0)，

(12)

式中，T為旋轉(zhuǎn)矩陣，滿足條件：

(13)

式中：φ和θ分別表示坐標(biāo)變換中x和z坐標(biāo)軸的旋轉(zhuǎn)角度.兩參數(shù)在地質(zhì)意義中分別表示地層方位角和傾角.

光滑矩陣的某一行元素就是一個(gè)光滑函數(shù)，用于描述不同介質(zhì)點(diǎn)之間的聯(lián)系，令此光滑函數(shù)為高斯光滑函數(shù)，則光滑矩陣中第i行第j列的元素為：

(14)

式中，σui,σvi,σwi分別是高斯函數(shù)在局部地質(zhì)坐標(biāo)系中不同方向上的標(biāo)準(zhǔn)差.實(shí)際應(yīng)用過程中，通常默認(rèn)地下介質(zhì)為均勻?qū)訝罱Y(jié)構(gòu)，認(rèn)為反射界面方向模型參數(shù)基本一致，垂直于反射界面方向的模型參數(shù)變化劇烈，因此水平方向?qū)?yīng)的光滑范圍相對(duì)較大，即σui、σvi要大于σwi，則相應(yīng)的標(biāo)準(zhǔn)差可以通過式(15)、(16)計(jì)算：

σui=σvi=abs(sin(φ))，

(15)

σwi=ξσui，

(16)

式中，ξ表示參數(shù)因子，通常取0.2～0.5的一個(gè)常數(shù)，至此將地質(zhì)構(gòu)造信息引入光滑矩陣.由于高斯函數(shù)的方差求取簡(jiǎn)單快速，可以作為一個(gè)已知量，所以構(gòu)建光滑矩陣的關(guān)鍵點(diǎn)就是求取地質(zhì)構(gòu)造的傾角和方位角信息.

(2)地層傾角方位角計(jì)算

在成像剖面中，可以拾取主要反射同相軸作為某一套地層的頂?shù)捉缑妫源藶榛€進(jìn)行界面之間層位的插值，通過計(jì)算頂?shù)捉缑娴南嚓P(guān)地層信息，然后結(jié)合地層接觸關(guān)系則可以依次得到整套地層的傾角和方位角信息.開展上述工作時(shí)，首先要判斷頂?shù)捉缑嬷g地層的接觸關(guān)系，通過地層關(guān)系實(shí)現(xiàn)整套地層信息的計(jì)算.根據(jù)地層沉積規(guī)律，頂?shù)捉缑嬷g的地層應(yīng)當(dāng)存在連續(xù)沉積、侵入式沉積、剝蝕沉積等地層關(guān)系，具體如圖2所示.

確定好頂?shù)椎貙咏佑|關(guān)系后，頂?shù)椎貙又g的其他地層可以通過插值進(jìn)行恢復(fù)，從而建立整個(gè)工區(qū)的地層架構(gòu).地層架構(gòu)一旦確立，可以根據(jù)幾何定理得到人工解釋層位數(shù)據(jù)的傾角方位角，地層傾角方位角的幾何關(guān)系如圖3所示.圖3中，θx表示x方向的視傾角，θy表示y方向的視傾角，θ表示地層的真傾角，ω表示地層的方位角.利用人工解釋地層數(shù)據(jù)信息，分別在x方向和y方向建立網(wǎng)格，目標(biāo)點(diǎn)位于中心，解釋數(shù)據(jù)位于網(wǎng)格交點(diǎn)上，則數(shù)據(jù)位置分布如圖4所示.

通過圖3和圖4，則目標(biāo)點(diǎn)的視傾角、真傾角和方位角的公式為：

tanθx=Δx/Δz，

(17)

tanθy=Δy/Δz，

(18)

(19)

(20)

式中，Δx、Δy和Δz分別表示目標(biāo)點(diǎn)和四周8個(gè)點(diǎn)不同方向上坐標(biāo)的差值加權(quán)結(jié)果.完成頂?shù)捉忉尩貙拥姆轿唤莾A角計(jì)算后，則解釋地層之間的地層傾角方位角可以采用圖5的原理進(jìn)行計(jì)算.

圖2 不同沉積模式示意圖(a) 等比例連續(xù)沉積； (b) 等厚底面侵入沉積； (c) 等厚頂面剝蝕沉積.Fig.2 Schematic diagram of different sedimentary modes(a) Equal proportion continuous deposition； (b) Equal thickness bottom intrusive sedimentation； (c)Equal thickness top denudation deposition.

圖3 地層傾角方位角幾何信息示意圖Fig.3 Schematic diagram of geometric information of formation dip angle and azimuth angle

圖4 層位網(wǎng)格示意圖Fig.4 Schematic diagram of horizon grid

圖5 不同沉積模式下地層傾角計(jì)算示意圖(a) 等比例連續(xù)沉積地層信息關(guān)系； (b) 等厚底面侵入沉積地層信息關(guān)系； (c) 等厚頂面剝蝕沉積地層信息關(guān)系.Fig.5 Calculation diagram of formation dip angle under different sedimentary modes(a) Information relationship of equal proportion continuous sedimentary strata； (b) Information relationship of sedimentary strata invaded by equal thickness bottom； (c) Information relation between denudation sedimentary strata of equal thickness top surface.

在圖5a表示等比例連續(xù)沉積，設(shè)定頂界面在平面坐標(biāo)某位置處的角度為φ1，同一平面坐標(biāo)位置處，底界面角度為φ3，則同一平面坐標(biāo)位置處的中間某一深度的地層角度φ2為：

(21)

式中，H1表示頂界面的深度值，H3表示底界面的深度值，H2表示目標(biāo)位置處的深度值，在圖5b中，地層接觸關(guān)系為等厚底面侵入關(guān)系，則中間地層的地質(zhì)信息和頂界面地層一致.在圖5c中，展示的是等厚頂面剝蝕地層的接觸關(guān)系，中間地層的地質(zhì)信息和底界面地層一致.結(jié)合圖5所示的不同地層接觸關(guān)系，運(yùn)用公式(21)可以計(jì)算整個(gè)工區(qū)的構(gòu)造信息.

(3)斷層控制正則化

公式(11)引入構(gòu)造信息后，可以得到較為準(zhǔn)確的速度反演結(jié)果，但斷層區(qū)域速度精度依舊有待提升，因此如何構(gòu)建一個(gè)斷層控制的目標(biāo)泛函是關(guān)鍵.

預(yù)條件模型正則化方法會(huì)使得斷層區(qū)域的速度邊界模糊，難以保持?jǐn)鄬觾蓚?cè)的速度差異性，因此有必要采用新方法提高斷層兩側(cè)速度的精度.考慮TV正則化具有保持邊緣信息的特征，可以在預(yù)條件模型正則化的基礎(chǔ)上引入TV正則化.TV正則化在重構(gòu)模型參數(shù)過程中，能夠得到更為準(zhǔn)確和豐富的參數(shù)信息，但是需要異常邊界信息的支持，從而加快收斂速度，提高收斂精度，因此可以在重構(gòu)算法中引入斷層邊界約束控制方法.根據(jù)當(dāng)前迭代成像結(jié)果中斷層邊界位置信息，構(gòu)建權(quán)重矩陣并將其用于下一次迭代，從而形成邊界約束的TV正則化方法，則改進(jìn)后正則化的目標(biāo)泛函可以表示為：

(22)

式中，‖‖1表示L1范數(shù)，D表示三維梯度算子稀疏矩陣，其主要作用是對(duì)速度模型中的每一個(gè)速度值做梯度運(yùn)算，在每次迭代過程中，使探測(cè)得到的參數(shù)邊界避免正則化，從而保持該位置附近的參數(shù)邊界；S表示預(yù)條件光滑矩陣.然而，在實(shí)際計(jì)算過程中，TV正則化存在著不穩(wěn)定性，因此需要將公式(22)轉(zhuǎn)為線性化方式進(jìn)行求解，然后通過多次迭代并不斷修正模型參數(shù)，使最終的反演結(jié)果不斷逼近真實(shí)值.

對(duì)于模型參數(shù)的修正，需要采用迭代加權(quán)方式求解，通過對(duì)一系列權(quán)系數(shù)不斷更新的加權(quán)最小二乘問題進(jìn)行計(jì)算，得到目標(biāo)函數(shù)的最優(yōu)解.采用加權(quán)方式將L1范數(shù)松弛轉(zhuǎn)為L(zhǎng)2范數(shù)，則公式(22)轉(zhuǎn)變?yōu)椋?/p>

(23)

式中，ε和β表示正則化參數(shù)，取值為[0,1]，用于調(diào)節(jié)不同約束項(xiàng)的權(quán)重，w表示對(duì)角矩陣，其公式可以表示為：

(24)

對(duì)于一個(gè)Nx×Ny×Nz的層析網(wǎng)格，當(dāng)i=1,2,…,Nx，j=1,2,…,Ny，k=1,2,…,Nz時(shí)，則梯度算子矩陣可以表示為：

(25)

由實(shí)際地質(zhì)特征可知，斷面左右的速度具有一定的差異性，要想保持?jǐn)鄬觾蓚?cè)速度突變的真實(shí)性，就需要公式(23)在計(jì)算速度值梯度算子時(shí)不能跨越斷層兩側(cè)，也就是位于斷面網(wǎng)格位置處的速度值不能和斷面兩側(cè)的速度值同時(shí)參與梯度計(jì)算，因此可以將上述斷層代表的邊界信息以權(quán)重矩陣的形式參與下一次迭代.即將此邊界與公式(23)正則化過程結(jié)合，基于上述思想，則構(gòu)建的對(duì)角權(quán)重?cái)?shù)組表達(dá)式為：

(26)

式中P(xi,yj,zk)表示xi,yj,zk位置處網(wǎng)格坐標(biāo)，F(xiàn)(Φ)表示解釋的斷層數(shù)據(jù)坐標(biāo)集，W表示三維矩陣.

由公式(26)可知，矩陣的建立不受斷層類型的影響，一旦獲取到斷層的位置坐標(biāo)，就可以建立對(duì)角權(quán)重矩陣.定義斷層控制的三維誘導(dǎo)矩陣Γ為：

(27)

將上述權(quán)重矩陣公式代入公式(23)，則公式變?yōu)椋?/p>

(28)

隨著速度場(chǎng)更新，斷層在空間上會(huì)發(fā)生飄移，斷層更新后的實(shí)際位置已經(jīng)偏離最初解釋的斷層位置，因此需要對(duì)斷層位置進(jìn)行推測(cè)識(shí)別，采用新的斷層信息對(duì)層析過程進(jìn)行約束控制.

為提升斷層的識(shí)別精度和提取效率，本文采用各向異性擴(kuò)散方程中的梯度結(jié)構(gòu)張量來提取地質(zhì)構(gòu)造特征，通過特征向量的情況來識(shí)別斷層，則三維網(wǎng)格中梯度結(jié)構(gòu)張量的計(jì)算公式為：

(29)

(30)

圖6 特征值和空間地質(zhì)結(jié)構(gòu)關(guān)系圖Fig.6 Relationship between eigenvalues and spatial geological structure

為提升流程效率，只需對(duì)目標(biāo)區(qū)域的主斷層區(qū)域進(jìn)行斷層控制模型更新，因此需要應(yīng)用人工解釋的斷層數(shù)據(jù)F(Φ)對(duì)計(jì)算得到的斷層位置集合Fn進(jìn)行預(yù)估，其中F表示斷層數(shù)據(jù)(三維)集合，Φ表示人工解釋斷層，n表示第n次迭代.速度模型的更新對(duì)于斷層位置偏移量影響相對(duì)較小，因此可以定義偏移量來時(shí)預(yù)估新的斷層位置集，如果推測(cè)出的斷層位置集合Fn∈F(Φ+Δf)，則認(rèn)為是目標(biāo)斷層區(qū)域.根據(jù)該處定義可知，即使人工解釋斷層刻畫位置存在誤差，只要在一定范圍內(nèi)，可以通過后續(xù)更新后的斷層來約束斷層位置.其中Δf表示定義的空間位置誤差項(xiàng)，然后將計(jì)算得到的斷層數(shù)據(jù)引入公式(24)和公式(27)，并結(jié)合公式(23)，則建立的斷層控制矩陣的目標(biāo)泛函可以表示為：

(31)

LW(mn)=(TΓn)Tdiag[max(ΓnTmn,γ)p-2]ΓnT，

(32)

+βLW(mn)mn，

(33)

(35)

(36)

通過上述方法，速度更新過程可以很好的保護(hù)斷層區(qū)域速度的差異性，可以解決常規(guī)層析反演方法在該區(qū)域中遇到的難題，得到斷層區(qū)域較為準(zhǔn)確的速度值.

2 理論模型

為驗(yàn)證本文方法，建立一個(gè)含有復(fù)雜斷層構(gòu)造的地質(zhì)模型，該模型含有多個(gè)斷層，然后在此基礎(chǔ)上開展正演，進(jìn)行速度分析后執(zhí)行疊前深度偏移.其中，真實(shí)速度模型剖面以及采用常規(guī)速度層析方法得到的偏移結(jié)果剖面如圖8所示.

圖7 斷控層析速度建模流程算法示意圖Fig.7 Schematic diagram of fault constrained tomography flow algorithm

對(duì)比圖8可知，常規(guī)層析速度建模方法得到的成像剖面，其構(gòu)造形態(tài)和原始速度模型大體一致，但是斷層下方區(qū)域成像質(zhì)量較差.因此該區(qū)域的成像質(zhì)量有待進(jìn)一步改善.在圖8成像剖面的基礎(chǔ)上，解釋層位和斷層數(shù)據(jù)，然后計(jì)算地層的傾角(所有地層方位角為同一數(shù)值)，則計(jì)算的地層傾角數(shù)據(jù)剖面如圖9所示.在上述解釋數(shù)據(jù)、傾角和方位角數(shù)據(jù)的基礎(chǔ)上，建立斷層約束算子，則速度更新前后的剖面對(duì)比如圖10所示.最后將基于斷控層析速度建模方法得到的速度應(yīng)用于PSDM，則兩種不同速度反演方法得到偏移成像結(jié)果如圖11所示.

對(duì)比圖10和11可知，采用斷控層析反演方法得到速度，能夠更好的描述斷層區(qū)域的速度分布情況，斷層兩邊的速度截?cái)嗪蛯?shí)際情況也較為匹配.采用斷控層析速度進(jìn)行PSDM，斷面成像更清晰，斷層下部的地層成像更為連續(xù)，成像質(zhì)量更好.

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

選取南海某工區(qū)的實(shí)際地震數(shù)據(jù)進(jìn)行測(cè)試，該區(qū)域斷裂系統(tǒng)發(fā)育，斷裂結(jié)構(gòu)復(fù)雜，導(dǎo)致很難獲取該工區(qū)準(zhǔn)確的速度場(chǎng)，從而使得偏移成像結(jié)果不可靠，嚴(yán)重影響該區(qū)域的勘探開發(fā)進(jìn)程.為解決上述問題，本文采用斷控層析速度建模方法獲取準(zhǔn)確的速度場(chǎng)，從而有效提高PSDM成像準(zhǔn)確度.

圖8 速度模型及PSDM成像剖面Fig.8 Velocity model and PSDM imaging profile

圖9 計(jì)算的地層傾角數(shù)據(jù)剖面Fig.9 Calculated dip data profile

(1)地層傾角方位角數(shù)據(jù)的計(jì)算

為準(zhǔn)確計(jì)算地層的傾角和方位角，首先需要根據(jù)地質(zhì)沉積規(guī)律來準(zhǔn)確判斷地層之間的接觸關(guān)系，其次分析區(qū)域的實(shí)際地質(zhì)情況，考慮地層和斷層、地層與地層、斷層與斷層之間的接觸關(guān)系，圖12展示一條二維線上主要斷層解釋數(shù)據(jù)分布情況.

由圖12可知，該區(qū)域地層和地層之間都是等比例沉積，所以只需要計(jì)算出區(qū)域主地層的傾角和方位角即可，然后利用公式(17)計(jì)算該區(qū)域所有地層傾角和方位角，計(jì)算的地層傾角和方位角如圖13和圖14所示.

圖10 不同速度反演方法結(jié)果對(duì)比(a) 常規(guī)層析方法反演結(jié)果； (b) 斷控層析反演結(jié)果.Fig.10 Comparison of results of different velocity inversion methods(a) Inversion results of conventional tomography； (b) Fault controlled tomography inversion results.

圖11 不同速度建模方法的到的PSDM剖面(a) 常規(guī)層析反演方法得到PSDM剖面； (b) 斷控層析反演方法得到PSDM剖面.Fig.11 PSDM profiles with different velocity modeling methods(a) PSDM profile obtained by conventional tomography inversion method； (b) PSDM profile obtained by fault controlled tomography inversion method.

圖12 解釋數(shù)據(jù)剖面Fig.12 Interpretation data profile

圖13 計(jì)算的地層傾角剖面示意圖Fig.13 Calculated dip profile

在速度更新過程中，需要采用梯度結(jié)構(gòu)張量對(duì)斷層的位置進(jìn)行重新識(shí)別.隨著速度的變化，部分?jǐn)鄬拥奈恢冒l(fā)生飄移，但是整體偏移位置不大.經(jīng)過速度更新后，采用梯度結(jié)構(gòu)張量來識(shí)別工區(qū)內(nèi)的斷層，則本研究工區(qū)中識(shí)別的斷層位置和原始人工解釋斷層的疊合情況如圖15所示.

由圖15可知，速度更新后，梯度結(jié)構(gòu)張量自動(dòng)識(shí)別的斷層位置和原始人工解釋斷層位置存在一些差異，但是位置的差異量不大.且根據(jù)識(shí)別結(jié)果可知，在部分區(qū)域，即使沒有人工解釋斷層，很多小斷層也被識(shí)別出來.為了盡可能提升計(jì)算效率，目前只需將解釋斷層附近被識(shí)別出斷層數(shù)據(jù)引入層析過程.根據(jù)解釋斷層的位置，選取附近更新后并被識(shí)別出的斷層數(shù)據(jù)，將其用于斷控層析速度建模，從而實(shí)現(xiàn)速度的更新.對(duì)于識(shí)別出來其他非目標(biāo)區(qū)域的斷層，其斷層數(shù)據(jù)則不用建立斷層約束算子.根據(jù)上述方式，可以利用初始解釋斷層數(shù)據(jù)匹配自動(dòng)識(shí)別的斷層，采用新的斷層數(shù)據(jù)建立斷層約束算子，將其引入改進(jìn)的目標(biāo)函數(shù)用于速度場(chǎng)的迭代更新.

圖14 計(jì)算的地層方位角剖面示意圖Fig.14 Calculated stratigraphic azimuth profile

圖15 層析反演前后的斷層邊界示意圖Fig.15 Fault boundary before and after tomographic inversion profile

(2)斷控層析速度建模以及偏移成像

利用提出的斷層控制正則化層析速度反演方法對(duì)南海復(fù)雜斷塊區(qū)進(jìn)行速度建模，在此基礎(chǔ)上進(jìn)行疊前深度偏移.圖16和圖17分別展示了常規(guī)層析速度、斷控層析速度的剖面對(duì)比和水平切片.

斷層區(qū)域由于斷層兩側(cè)對(duì)接巖性的年代差異，會(huì)表現(xiàn)出速度突變現(xiàn)象.然而，在圖16a和圖17a中，斷層兩側(cè)速度是平滑漸變，斷層兩側(cè)速度差異不明顯.斷層約束層析反演速度建模結(jié)果中，斷層兩側(cè)速度突變明顯，斷面等其他物理性質(zhì)影響體現(xiàn)的更為準(zhǔn)確.為進(jìn)一步分析兩種層析結(jié)果，對(duì)其成像結(jié)果再次對(duì)比，結(jié)果如圖18所示.

對(duì)比圖18可以發(fā)現(xiàn)，在常規(guī)方法的成像剖面中，斷層區(qū)域(藍(lán)色線框和黑色線框區(qū)域)的反射同相軸出現(xiàn)扭曲錯(cuò)段的現(xiàn)象，甚至出現(xiàn)一條垂直的虛假斷層，上述現(xiàn)象的存在，嚴(yán)重干擾對(duì)斷層區(qū)域地質(zhì)狀況的認(rèn)識(shí).經(jīng)過分析可知，常規(guī)速度建模方法導(dǎo)致斷層附近的速度不準(zhǔn)確，使得成像剖面出現(xiàn)虛假的垂直斷層和地層扭曲錯(cuò)斷的現(xiàn)象.而在斷控層析反演成像結(jié)果中，虛假的垂直斷層消失，反射同相軸更加清楚和平整，表明斷層附近的速度更為準(zhǔn)確，剖面成像質(zhì)量提升明顯.

圖16 不同速度建模方法的地震-速度疊合剖面(a) 常規(guī)層析速度； (b) 斷控層析速度.Fig.16 Profile of different velocity inversion methods(a) Conventional chromatographic velocity；(b) Fault constrained tomography velocity.

圖17 不同速度建模方法的地震-速度疊合水平切片(a) 常規(guī)層析速度； (b) 斷控層析速度.Fig.17 Time slices of different velocity inversion methods(a) Conventional chromatographic velocity； (b) Fault constrained tomography velocity.

圖18 不同速度建模方法的成像剖面對(duì)比示意圖(a) 常規(guī)層析PSDM剖面； (b) 斷控層析的PSDM剖面.Fig.18 Comparison diagram of imaging sections with different inversion velocities(a) Conventional tomographic PSDM profile；(b) PSDM profile of fault controlled tomography.

為進(jìn)一步分析本文方法對(duì)于成像的影響，提取兩種方法偏移成像結(jié)果的相干體屬性進(jìn)行對(duì)比，結(jié)果如圖19所示.

從相干體切片中，斷控層析反演速度建模方法的成像結(jié)果優(yōu)勢(shì)更加明顯.在圖19a的常規(guī)層析速度成像結(jié)果中，圓圈中斷層的識(shí)別較為模糊，地層的反射也較為紊亂，難以識(shí)別斷層的組合.然而，在圖19b結(jié)果中，地層反射規(guī)律性強(qiáng)，斷層位置識(shí)別清晰.

為檢查速度反演質(zhì)量，進(jìn)一步量化反演速度結(jié)果，利用穿過斷層的某實(shí)際鉆井?dāng)?shù)據(jù)進(jìn)行驗(yàn)證，井軌跡位置如圖20所示.然后采用相同的地震子波和不同的地震速度模擬地震記錄，將模擬結(jié)果和測(cè)井速度的合成地震記錄進(jìn)行相關(guān)，通過相關(guān)性來判斷當(dāng)前速度的精確度.圖21展示的是測(cè)井?dāng)?shù)據(jù)合成記錄同兩種速度的正演結(jié)果對(duì)比，結(jié)果表明，斷控層析速度的合成地震記錄和測(cè)井地震記錄的相關(guān)系數(shù)為0.75，常規(guī)層析速度的合成地震記錄和測(cè)井地震記錄的相關(guān)系數(shù)為0.58，證實(shí)斷控層析速度建模方法的速度結(jié)果和實(shí)際鉆井速度吻合度更高，最符合實(shí)際速度場(chǎng)分布情況.

4 結(jié)論

本文提出了一種基于斷層控制的層析反演速度建模方法，該方法通過地震解釋數(shù)據(jù)計(jì)算相關(guān)模型參數(shù)，建立斷層約束的正則化約束算子，以此為基礎(chǔ)建立斷層控制的層析反演目標(biāo)泛函，利用斷層數(shù)據(jù)作為速度建模的先驗(yàn)信息，采用迭代最小二乘優(yōu)化方法逼近L1范數(shù)的目標(biāo)泛函，提高速度建模的精度.應(yīng)用本文方法獲取的速度模型，斷層區(qū)域的速度分布更加符合地質(zhì)認(rèn)識(shí)，模型與真實(shí)情況吻合度更高.將反演得到的速度場(chǎng)用于PSDM，提高了復(fù)雜斷塊區(qū)域的成像質(zhì)量，為后續(xù)目標(biāo)評(píng)價(jià)提供了可靠的地震成像數(shù)據(jù)，從而有助于提高復(fù)雜區(qū)油氣勘探成功率.由于相關(guān)研究條件限制，本方法的模型參數(shù)表征采用了矩形網(wǎng)格方式，因此在描述模型邊界的精度方面依舊存在提升空間，下一步研究計(jì)劃采用三角網(wǎng)格進(jìn)行模型參數(shù)表征，開展層析反演的正則化研究工作.

圖19 不同反演速度的成像剖面相干體提取結(jié)果對(duì)比(a) 常規(guī)層析成像結(jié)果提取的相干體； (b) 斷控層析成像結(jié)果提取的相干體.Fig.19 Comparison of coherent volume extraction results of imaging profiles with different inversion velocities(a) Coherent volume extracted from conventional tomography results； (b) Coherent volume extracted from fault control tomography results.

圖20 穿越斷層的井軌跡示意圖Fig.20 Schematic diagram of well trajectory crossing fault

圖21 不同速度的合成記錄對(duì)比示意圖(a) 測(cè)井聲波數(shù)據(jù)； (b) 測(cè)井密度數(shù)據(jù)； (c) 測(cè)井合成記錄； (d) 常規(guī)層析反演速度合成記錄； (e) 斷控層析速度合成記錄.Fig.21 Forward comparison diagram of different velocity(a) Logging acoustic data； (b) Logging density data； (c) Logging synthetic record； (d) Conventional tomographic inversion velocity synthesis record； (e) Fault controlled tomography velocity synthesis record.