數據驅動的塊排列正則化多道疊前地震反演

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

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

0 引言

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

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

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

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

1 基本原理

1.1 塊排列正則化

從疊后地震數據中提取一個塊排列矩陣P，期望將反演參數重新排序為分段序列.與傳統的直接平滑方法不同，該方法首先根據觀測數據的局部相似性，將特征相近的采樣點重新排列在一起，然后再采用約束的方法去除噪聲，最后重排回原始位置.為了得到塊排序矩陣，構造目標函數為：

(1)

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

基于提取的塊排列矩陣，可以構建具有多道構造約束的塊排列正則項：

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

(3)

可以看到，當相鄰采樣點存在邊界時，地震記錄振幅的二階差分相對較大，即分母比較大，使相應的權重系數變小，達到減弱邊界平滑的作用.

圖1 地震記錄塊的提取(a) 二維實際資料； (b) 采用滑動窗從地震記錄剖面提取的地震記錄塊.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方程可以描述振幅隨偏移距變化的現象，但它具有嚴重的非線性，在正演模擬和反演時計算復雜(Zhi et al., 2016).本文采用Fatti近似方程(Fatti et al., 1994)進行疊前地震反演:

(4)

圖2 塊排列矩陣的作用(a) 合成疊后地震記錄； (b) 真實縱波速度模型； (c) 真實模型所有道連接的一維速度序列； (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近似三項系數隨入射角度變化；(b)Fatti近似三項隨入射角度變化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道不同入射角的觀測記錄，W表示地震子波矩陣，A表示方程(4)中Fatti近似公式的系數，mj表示第j道反演的縱波阻抗和橫波阻抗，其可以展開為：

(6)

考慮到橫向連續性，將疊前地震反演目標函數寫成多道的形式：

(7)

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

(8)

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

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

(2)采用井震標定的方法，將測井數據進行時深轉換，并估計用于反演的地震子波；

(3)將井數據進行低通濾波，沿層位進行插值，獲得低頻的縱波阻抗，橫波阻抗和密度的低頻初始模型；

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

2 合成數據測試

為了說明基于塊排列正則化的多道疊前地震反演方法的可靠性，將該方法應用于Marmousi模型的疊前反演中，并與常規的基于模型約束的單道反演方法對比.在合成數據測試中，采用均方根誤差衡量反演結果的精度：

(9)

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

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

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

圖4 觀測地震記錄(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 真實模型和初始模型真實(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 反演結果傳統基于模型約束的(a)縱波阻抗和(b)橫波阻抗反演結果；塊排列正則化的(c)縱波阻抗和(d)橫波阻抗反演結果.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 局部放大的傳統(a)縱波阻抗和(b)橫波阻抗反演結果，塊排列正則化(c)縱波阻抗和(d)橫波阻抗反演結果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隨超參數λP和λS的變化； (b) 目標函數值隨迭代次數的變化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)橫波阻抗的初始模型；傳統基于模型約束的(c)縱波阻抗和(d)橫波阻抗反演結果；基于塊排列正則化的(e)縱波阻抗和(f)橫波阻抗的反演結果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) 傳統縱波阻抗反演結果局部放大； (c)和(d) 塊排列正則化縱波阻抗反演結果局部放大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) 傳統橫波阻抗反演結果局部放大； (c)和(d) 塊排列正則化橫波阻抗反演結果局部放大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)橫波阻抗的反演結果與測井數據比較Fig.13 Comparison between borehole data and inversion results of (a) IP and (b) IS

3 實際數據測試

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

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

4 結論

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