王寧， 周輝， 王玲謙， 趙海波

1 東北石油大學(xué)地球科學(xué)學(xué)院， 大慶 163318 2 中國(guó)石油大學(xué)(北京)， 北京 102249 3 大慶油田勘探開發(fā)研究院， 大慶 163318

0 引言

地震反演作為地震勘探中一種重要的解釋方法，將地震觀測(cè)數(shù)據(jù)轉(zhuǎn)化為彈性參數(shù)，可以有效去除地震觀測(cè)數(shù)據(jù)中帶限子波的影響，補(bǔ)充低頻和高頻信息(Tarantola, 1987; Shuey, 1985; 高剛等, 2013; 曹俊興等, 2011).根據(jù)反演得到的彈性參數(shù)和不同的巖石物理模型，可以獲得孔隙度、泥質(zhì)含量等各種巖石物理參數(shù)(Karimi et al., 2010; Zong et al., 2012, 2013; 印興耀等, 2014；白俊雨，2020).傳統(tǒng)的疊后地震反演只能提供波阻抗剖面，不足以描述流體類型和地下精細(xì)構(gòu)造(Sharma and Chopra, 2015；張巖等，2021；周單等，2022).相比之下，疊前地震數(shù)據(jù)包含了振幅隨偏移距變化(AVO)的信息，且與縱波速度、橫波速度和密度的反射系數(shù)有關(guān)(Aki and Richards, 1980).因此，疊前地震反演為流體預(yù)測(cè)和巖性識(shí)別提供了更豐富、可靠的信息(宗兆云等, 2012; Chen et al., 2020; Liu et al., 2020a, b; Chen et al., 2021；吳奎等，2021).

由于地震觀測(cè)數(shù)據(jù)的帶限性和噪聲，地震反演存在嚴(yán)重的不適定性(Varela et al., 2006).為了獲得穩(wěn)定、準(zhǔn)確的反演結(jié)果，正則化方法可以對(duì)反演結(jié)果施加特定的約束，從而獲得具有期望特征的反演結(jié)果.目前有兩種典型的正則化地震反演方法，包括平滑約束和塊約束(VanDecar and Snieder, 1994; Ulrych and Sacchi, 2005；曹丹平等, 2010；印興耀等, 2020a；楊俊等，2020).Tikhonov(1963)提出了L2范數(shù)來(lái)約束反演參數(shù)，可以根據(jù)初始模型的低頻趨勢(shì)得到平滑的反演結(jié)果.Buland和Omre(2003)將貝葉斯框架與地震反演相結(jié)合，假設(shè)反演參數(shù)服從高斯先驗(yàn)分布，使反演結(jié)果平滑.然而，這些平滑約束會(huì)模糊邊界和斷層，并不能提高反演結(jié)果的分辨率.為了加強(qiáng)反演結(jié)果的不連續(xù)性，目前提出了多種形式的塊約束.Alemie和Sacchi(2011)引入貝葉斯框架下的三元柯西概率分布，基于稀疏的假設(shè)條件，重構(gòu)地下反射系數(shù)序列，提高反演結(jié)果的分辨率.Zhang等(2013)將L1范數(shù)引入疊前反射系數(shù)重構(gòu)中，通過(guò)提供準(zhǔn)確的地下反射系數(shù)序列，增強(qiáng)了反演彈性參數(shù)的塊狀特征，更清晰地刻畫了反演結(jié)果中的邊緣和斷層.塊狀約束的主要挑戰(zhàn)和局限性是地下薄層和弱反射層的預(yù)測(cè)(Wang et al., 2019).此外，這些方法都是逐道進(jìn)行的，忽略了地下構(gòu)造的空間連續(xù)性.單道處理方法降低了反演剖面的信噪比，不可避免地干擾了進(jìn)一步的解釋和儲(chǔ)層預(yù)測(cè)(Ma et al., 2017, 2019).

目前，有許多研究致力于多道同時(shí)反演和反射系數(shù)重構(gòu)(Gholami and Sacchi, 2013; Hamid et al., 2018; Huang et al., 2021; 霍國(guó)棟等, 2017; 印興耀等, 2020b).Gholami(2015, 2016)利用全變分約束進(jìn)行多道非線性波阻抗反演，該方法可以使反演結(jié)果具有塊狀特征，增強(qiáng)其空間連續(xù)性.Hamid和Pidlisecky(2015)將所有道連接成一道，并對(duì)相鄰道反演結(jié)果施加水平平滑約束.當(dāng)真實(shí)情況違背事先的假設(shè)條件，這些多道反演方案不能提供準(zhǔn)確的反演結(jié)果.因此，后續(xù)發(fā)展出了數(shù)據(jù)驅(qū)動(dòng)的多道地震反演方法.Karimi(2015)利用平面波分解方法提取地下構(gòu)造傾角，然后利用光滑算子沿局部斜率約束反演結(jié)果.Hamid和Pidlisecky(2016)利用觀測(cè)數(shù)據(jù)計(jì)算地震記錄同相軸的局部?jī)A角，并將傾角信息和測(cè)井?dāng)?shù)據(jù)引入到地震反演中，提高反演結(jié)果的橫向連續(xù)性.Huang等(2022)從地震數(shù)據(jù)中估算了地震斜率屬性，并將其納入全變分正則化中，獲得高分辨率反演結(jié)果.

受數(shù)據(jù)驅(qū)動(dòng)多道反演方法優(yōu)勢(shì)的啟發(fā)，本文將Wang等(2021)提出的多道疊后地震波阻抗反演方法擴(kuò)展到疊前地震反演.與以往的工作不同的是，該方法沒有直接從疊前地震數(shù)據(jù)中提取構(gòu)造信息.考慮到疊后地震數(shù)據(jù)具有更高的信噪比，本文利用疊后地震記錄的局部相似性，提取了記錄構(gòu)造延展的塊排列矩陣.準(zhǔn)確的構(gòu)造信息可以有效降低疊前反演的病態(tài)性，獲得穩(wěn)定可靠的反演結(jié)果.此外，與多道疊后地震反演不同，內(nèi)存和計(jì)算效率是多道疊前地震反演的主要限制因素.本文將逐道地震記錄不匹配項(xiàng)與多道構(gòu)造約束相結(jié)合，對(duì)原目標(biāo)函數(shù)進(jìn)行優(yōu)化，可以節(jié)省大型正演矩陣的內(nèi)存，提高計(jì)算效率.

1 基本原理

1.1 塊排列正則化

從疊后地震數(shù)據(jù)中提取一個(gè)塊排列矩陣P，期望將反演參數(shù)重新排序?yàn)榉侄涡蛄?與傳統(tǒng)的直接平滑方法不同，該方法首先根據(jù)觀測(cè)數(shù)據(jù)的局部相似性，將特征相近的采樣點(diǎn)重新排列在一起，然后再采用約束的方法去除噪聲，最后重排回原始位置.為了得到塊排序矩陣，構(gòu)造目標(biāo)函數(shù)為：

(1)

其中D表示一階差分矩陣，x表示待反演參數(shù).該目標(biāo)函數(shù)的求解等價(jià)于傳統(tǒng)的TSM(Travelling Salesman)最優(yōu)化問題(Cormen et al., 2001).每個(gè)地震記錄塊中各個(gè)采樣點(diǎn)的振幅可以等價(jià)于TSM問題里的城市坐標(biāo)，商人需要一次性通過(guò)各個(gè)城市并選取最短的路徑.隨機(jī)最近鄰啟發(fā)式算法是求解該優(yōu)化問題的一種常用方法，本文采用該方法獲取塊排列矩陣.圖2展示了塊排列矩陣應(yīng)用到推覆體模型后的結(jié)果.可以觀察到，二維速度剖面經(jīng)過(guò)塊排列矩陣處理后可以轉(zhuǎn)換為一維塊狀序列.

基于提取的塊排列矩陣，可以構(gòu)建具有多道構(gòu)造約束的塊排列正則項(xiàng)：

r(m)=‖QLPm‖1,(2)

(3)

可以看到，當(dāng)相鄰采樣點(diǎn)存在邊界時(shí)，地震記錄振幅的二階差分相對(duì)較大，即分母比較大，使相應(yīng)的權(quán)重系數(shù)變小，達(dá)到減弱邊界平滑的作用.

圖1 地震記錄塊的提取(a) 二維實(shí)際資料； (b) 采用滑動(dòng)窗從地震記錄剖面提取的地震記錄塊.Fig.1 The seismic patch extraction from seismic profile(a) The 2D field seismic data; (b) The seismic patches decomposed from the seismic profile with a sliding window.

1.2 多道疊前地震反演

Zoeppritz方程可以描述振幅隨偏移距變化的現(xiàn)象，但它具有嚴(yán)重的非線性，在正演模擬和反演時(shí)計(jì)算復(fù)雜(Zhi et al., 2016).本文采用Fatti近似方程(Fatti et al., 1994)進(jìn)行疊前地震反演:

(4)

圖2 塊排列矩陣的作用(a) 合成疊后地震記錄； (b) 真實(shí)縱波速度模型； (c) 真實(shí)模型所有道連接的一維速度序列； (d) 塊排列矩陣作用后的一維速度序列.Fig.2 The function of patch-ordering matrix(a) The synthetic seismic data; (b) True P-wave velocity model; (c) The contaminated 1D sequence of the true model; (d) The reordered 1D velocity sequence via patch-ordering scheme.

圖3 (a)Fatti近似三項(xiàng)系數(shù)隨入射角度變化；(b)Fatti近似三項(xiàng)隨入射角度變化Fig.3 (a) The amplitude variation with incident angle of three coefficients inFatti approximation; (b) The amplitude variation with incident angle of three terms in Fatti approximation

基于褶積模型，正演模擬可以寫成矩陣的形式：

dj=WADmj,(5)

其中，dj表示第j道不同入射角的觀測(cè)記錄，W表示地震子波矩陣，A表示方程(4)中Fatti近似公式的系數(shù)，mj表示第j道反演的縱波阻抗和橫波阻抗，其可以展開為：

(6)

考慮到橫向連續(xù)性，將疊前地震反演目標(biāo)函數(shù)寫成多道的形式：

(7)

為了直接使用L-BFGS迭代類算法求解目標(biāo)函數(shù)(Nocedal and Wright, 2006)，將目標(biāo)函數(shù)中L1范數(shù)寫成平滑的形式：

(8)

基于塊排列正則化的多道疊前地震反演，其工作流程總結(jié)如下：

(1)利用高信噪比疊后地震記錄，提取塊排列矩陣，自適應(yīng)地記錄地下構(gòu)造信息；

(2)采用井震標(biāo)定的方法，將測(cè)井?dāng)?shù)據(jù)進(jìn)行時(shí)深轉(zhuǎn)換，并估計(jì)用于反演的地震子波；

(3)將井?dāng)?shù)據(jù)進(jìn)行低通濾波，沿層位進(jìn)行插值，獲得低頻的縱波阻抗，橫波阻抗和密度的低頻初始模型；

(4)使用L-BFGS算法求解多道疊前地震反演的目標(biāo)函數(shù).

2 合成數(shù)據(jù)測(cè)試

為了說(shuō)明基于塊排列正則化的多道疊前地震反演方法的可靠性，將該方法應(yīng)用于Marmousi模型的疊前反演中，并與常規(guī)的基于模型約束的單道反演方法對(duì)比.在合成數(shù)據(jù)測(cè)試中，采用均方根誤差衡量反演結(jié)果的精度：

(9)

圖4a為Marmousi模型合成的疊后地震記錄，使用主頻為30 Hz，采樣間隔為1 ms的雷克子波.圖4b—d展示了添加20%隨機(jī)噪聲的疊前地震角道集.真實(shí)的Marmousi模型縱波阻抗和橫波阻抗如圖5a、b所示.初始模型通過(guò)采用81×81的平滑窗得到，如圖5c、d所示.

為了便于比較，本文也提供了傳統(tǒng)的基于模型約束的反演結(jié)果，如圖6a、b所示.可以看到，傳統(tǒng)的反演結(jié)果存在嚴(yán)重的噪聲干擾，很難刻畫地下構(gòu)造信息.相比之下，提出的方法可以提供更干凈準(zhǔn)確的反演剖面.為了更直觀地比較傳統(tǒng)方法和提出方法的反演結(jié)果，對(duì)反演結(jié)果進(jìn)行局部放大，如圖7所示.傳統(tǒng)反演方法得到的縱波阻抗和橫波阻抗的均方根誤差(Root Mean Square Error，RMSE)分別為747.0和565.4.提出方法的反演結(jié)果，其對(duì)應(yīng)的RMSE分別為568.56和506.10.

在基于塊排列正則化的多道疊前地震反演中，目標(biāo)函數(shù)包含超參數(shù)λP和λS，本文采用試錯(cuò)法確定超參數(shù)，得到最優(yōu)的反演結(jié)果，雖然目前有很多超參數(shù)確定的方法，但對(duì)不同的目標(biāo)函數(shù)，其方法的適用性依然存在問題(張宏兵和楊長(zhǎng)春, 2003; 張宏兵等, 2005).RMSE隨超參數(shù)λP和λS的變化如圖8a所示，其中叉號(hào)為縱波阻抗的RMSE，圓圈為橫波阻抗的RMSE.理論上，較大的超參數(shù)會(huì)導(dǎo)致反演結(jié)果平滑甚至模糊，很難刻畫邊界和斷層等細(xì)節(jié)信息.過(guò)小的超參數(shù)不能有效壓制反演結(jié)果中的噪聲.為了進(jìn)一步驗(yàn)證L-BFGS算法在求解目標(biāo)函數(shù)過(guò)程中的收斂性，記錄了目標(biāo)函數(shù)值隨迭代次數(shù)的變化，如圖8b所示，可以看出目標(biāo)函數(shù)隨著迭代次數(shù)的增加，逐漸收斂到一個(gè)穩(wěn)定的值，因此求解多道反演目標(biāo)函數(shù)的策略是有效的.

圖4 觀測(cè)地震記錄(a) 疊后地震記錄；角度為(b)0°； (c) 15°； (d) 30°的含噪角道集.Fig.4 The observed seismic record(a) The post-stack seismic data; The noisy angle gathers of (b) 0°, (c) 15°, and (d) 30°.

圖5 真實(shí)模型和初始模型真實(shí)(a)縱波阻抗和(b)橫波阻抗模型； (c) 縱波阻抗和(d)橫波阻抗的初始模型.Fig.5 The true and initial modelsThe true (a) IP and (b) IS models; The initial (c) IP and (d) IS models.

圖6 反演結(jié)果傳統(tǒng)基于模型約束的(a)縱波阻抗和(b)橫波阻抗反演結(jié)果；塊排列正則化的(c)縱波阻抗和(d)橫波阻抗反演結(jié)果.Fig.6 The inversion resultsThe conventional model-based (a) IP and (b) IS inversion results; The patch-ordering regularized (c) IP and (d) IS inversion results.

圖7 局部放大的傳統(tǒng)(a)縱波阻抗和(b)橫波阻抗反演結(jié)果，塊排列正則化(c)縱波阻抗和(d)橫波阻抗反演結(jié)果Fig.7 Zoomed-in view of conventional (a) IP and (b) IS inversion results,patch-ordering regularized (c) IP and (d) IS inversion results

圖8 (a) RMSE隨超參數(shù)λP和λS的變化； (b) 目標(biāo)函數(shù)值隨迭代次數(shù)的變化Fig.8 (a) The RMSE variation with hyper-parameters λP and λS; (b) The value of objective function versus iterations

圖9 (a) 疊后地震記錄； (b) 小角度、(c)中角度和(d)大角度部分疊加角道集Fig.9 (a) The post-stack seismic data; Partially stacked angle gathers of (b) small-angle range,(c) middle-angle range, and (d) large-angle range

圖10 (a)縱波阻抗和(b)橫波阻抗的初始模型；傳統(tǒng)基于模型約束的(c)縱波阻抗和(d)橫波阻抗反演結(jié)果；基于塊排列正則化的(e)縱波阻抗和(f)橫波阻抗的反演結(jié)果Fig.10 Initial model of (a) IP and (b) IS; Model-based inversion results of (c) IP and (d) IS; Patch-ordering regularized inversion results of (e) IP and (f) IS

圖11 (a)和(b) 傳統(tǒng)縱波阻抗反演結(jié)果局部放大； (c)和(d) 塊排列正則化縱波阻抗反演結(jié)果局部放大Fig.11 Zoomed-in view of (a) and (b) conventional IP inversion results,(c) and (d) patch-ordering regularized IP inversion results

圖12 (a)和(b) 傳統(tǒng)橫波阻抗反演結(jié)果局部放大； (c)和(d) 塊排列正則化橫波阻抗反演結(jié)果局部放大Fig.12 Zoomed-in view of (a) and (b) conventional IS inversion results,(c) and (d) patch-ordering regularized IS inversion results

圖13 (a)縱波阻抗和(b)橫波阻抗的反演結(jié)果與測(cè)井?dāng)?shù)據(jù)比較Fig.13 Comparison between borehole data and inversion results of (a) IP and (b) IS

3 實(shí)際數(shù)據(jù)測(cè)試

將提出的方法應(yīng)用于實(shí)際的地震數(shù)據(jù)中，部分疊加角道集和疊后地震記錄如圖9所示.小角度，中角度和大角度部分疊加角道集角度范圍為0°～15°、10°～25°和20°～35°.縱波阻抗和橫波阻抗的初始模型通過(guò)測(cè)井?dāng)?shù)據(jù)10～15 Hz低通濾波和插值得到，如圖10a、b所示.塊排列矩陣通過(guò)疊后地震記錄提取，然后，利用提取的塊排列矩陣構(gòu)建多道反演的目標(biāo)函數(shù).基于塊排列正則化的反演結(jié)果如圖10e、f所示.可以看到在反演剖面上有一個(gè)大的阻抗異常體，反演結(jié)果清晰地刻畫出研究區(qū)域潛山構(gòu)造.為了進(jìn)一步說(shuō)明提出方法的優(yōu)勢(shì)，傳統(tǒng)的基于模型約束的反演結(jié)果如圖10c、d所示.可以看到反演的剖面信噪比較低，干擾了地下層位和構(gòu)造的刻畫.為了更直觀地比較，將圖10中反演結(jié)果的淺層和深層進(jìn)行局部放大.圖11a、c展示了淺層基于模型和塊排列正則化反演的縱波阻抗，圖11b、d展示了深部縱波阻抗反演結(jié)果.可以看到，提出的多道反演策略可以在提高分辨率的同時(shí)加強(qiáng)反演結(jié)果的橫向連續(xù)性.此外，基于模型約束和塊排列正則化的橫波阻抗反演結(jié)果如圖12所示，可以得到相同的結(jié)論.

為了定量評(píng)價(jià)塊排列正則化反演策略的可靠性，反演結(jié)果與位于CDP 34處的測(cè)井?dāng)?shù)據(jù)進(jìn)行比較，如圖13所示.塊排列正則化反演得到的縱波阻抗和橫波阻抗(藍(lán)色曲線)和測(cè)井曲線(紅色曲線)匹配較好.它們的相關(guān)系數(shù)為0.85和0.84.黑色曲線表示基于模型約束的反演結(jié)果，可以發(fā)現(xiàn)，基于模型約束的反演結(jié)果雖然低頻趨勢(shì)與井?dāng)?shù)據(jù)基本一致，但信噪比較低，出現(xiàn)很多不規(guī)則抖動(dòng).基于模型約束的反演結(jié)果與測(cè)井曲線的相關(guān)系數(shù)分別為0.81和0.79.因此，提出的方法可以提供可靠的反演結(jié)果，并且能更清晰地刻畫地下構(gòu)造的空間展布.

4 結(jié)論

本文提出了一種基于塊排列正則化的多道疊前地震反演方法.該方法利用疊后地震數(shù)據(jù)提取出一個(gè)可以記錄構(gòu)造方向的塊排列矩陣.求取塊排列矩陣的目標(biāo)函數(shù)等價(jià)于傳統(tǒng)的TSM問題，該目標(biāo)函數(shù)可以用隨機(jī)最近鄰啟發(fā)式方法有效地求解.然后利用提取的塊排列矩陣構(gòu)造正則項(xiàng)，考慮到重新排序的彈性參數(shù)序列的局部空間連續(xù)性，采用變權(quán)重系數(shù)的策略，在斷層或分界面處采用較小的權(quán)重系數(shù)，在連續(xù)構(gòu)造處采用較大權(quán)重系數(shù)，增強(qiáng)連續(xù)性.因此，這種正則化方法可以在保護(hù)邊緣的同時(shí)降低反演結(jié)果中的噪聲.目標(biāo)函數(shù)由觀測(cè)數(shù)據(jù)不匹配項(xiàng)、縱波阻抗和橫波阻抗塊排列正則項(xiàng)構(gòu)成.在合成數(shù)據(jù)測(cè)試中，與傳統(tǒng)的基于模型約束的反演結(jié)果相比，塊排列正則化反演結(jié)果中的縱波阻抗和橫波阻抗的變化更清晰、更準(zhǔn)確.在實(shí)際資料測(cè)試中，提出方法的反演結(jié)果可以穩(wěn)定、準(zhǔn)確地描述地下地層的延展.因此，該方法可以作為今后儲(chǔ)層預(yù)測(cè)的有效工具.此外，該方法有望應(yīng)用于地震數(shù)據(jù)處理和全波形反演中，在增強(qiáng)信號(hào)或反演結(jié)果空間連續(xù)性的同時(shí)，保護(hù)邊緣.