基于貝葉斯理論的反射率法煤田地震波阻抗反演

楊 真，高振宇，劉興業(yè),2

基于貝葉斯理論的反射率法煤田地震波阻抗反演

楊 真1，高振宇1，劉興業(yè)1,2

(1. 神華神東煤炭集團有限責任公司布爾臺煤礦，內蒙古 鄂爾多斯 017000；2. 西安科技大學 地質與環(huán)境學院，陜西 西安 710054)

地震波阻抗反演可以直接提供地下介質的彈性信息，為后續(xù)的煤層識別，水通道預測以及陷落柱識別等提供可靠的資料。在面對復雜地質條件時，現(xiàn)有的波阻抗反演方法精度較低，不能滿足實際勘探的需求。反射率法通過求解彈性波動方程獲得彈性參數(shù)信息，能夠模擬全波場響應以及地震波的各種傳播效應，能更精確地考慮透射損失、多次波等引起的地震波振幅和相位變化，精度較高，且計算效率高于其他波動方程求解方法。因此，基于貝葉斯框架，建立了反射率法的反演目標函數(shù)并通過引入柯西分布作為先驗分布提高反演結果的分辨率，形成了一種基于貝葉斯的反射率法波阻抗反演方法。將該方法應用于布爾臺煤礦地震數(shù)據(jù)中，獲得的波阻抗反演結果與測井資料匹配度較好，相比于常規(guī)方法分辨率更高，精度高，能夠有效識別薄煤層和深部煤層分布，有效驗證了新方法的可行性和有效性。該方法能為利用地震反演技術識別煤層及陷落柱等小構造，預測頂?shù)装逅植继峁┯辛Φ馁Y料。

地震反演；反射率法；貝葉斯理論；波傳播效應；煤田勘探

地震波阻抗反演是將地震資料轉換為波阻抗剖面的一種技術，它將地震資料、測井數(shù)據(jù)、地質解釋相結合，充分利用了測井資料具有較高的垂向分辨率和地震資料具有較高的橫向分辨率的特點，反演獲得的波阻抗剖面，不僅便于解釋人員將地震資料與測井資料連接對比，而且對地層相關參數(shù)的變化研究提供可靠資料，從而指導資源的勘探開發(fā)。

煤炭資源的勘探開發(fā)離不開地震勘探技術，其中波阻抗反演可以預測煤層位置，指示煤系氣分布情況，分析煤層厚度的變化，估算資源儲量，查明陷落柱發(fā)育情況，同時對富水區(qū)也能有顯著的指示作用，并對后續(xù)的進一步研究提供直觀可靠的資料[1-4]。因而，在煤炭及其相關資源的勘探開發(fā)中越來越受到重視。但是，目前的波阻抗反演算法都是基于簡單的褶積模型進行的，通常需要一定的假設條件。當不滿足這些假設條件時，嚴重制約反演精度，同時地震資料的帶限性質會影響反演分辨率，從而影響勘探效果，增加勘探開發(fā)風險[5-6]。此外，在處理薄煤層問題時，由于非均質性強，轉換波及透射波損失嚴重，這些方法沒有考慮地震波的損失，從而制約反演精度。

彈性波動方程是一種較為精確地表達地震波傳播的方法，但受計算效率、反演穩(wěn)定性和對地震資料要求高等因素的影響，實際應用面臨諸多挑戰(zhàn)。而反射率法是一種以一維遞推方式求解彈性波動方程解析解的方法，能對全波場模擬。除了一次反射，反射率法還能夠模擬透射波、多種轉換波和多次波，考慮了地震波相位變化以及透射損失的影響，且反射率法計算精度較高，計算成本適中，具有較強的實用價值[7-9]。

反射率法由K. Fuchs等[10]提出，經(jīng)過不斷的完善[11-14]，已經(jīng)較為成熟并被成功應用于地震反演中。Zhao H.等[15]基于Kennett遞歸矩陣，通過Levenberg- Marquardt算法實現(xiàn)了頻率波數(shù)反演。W. P. Gouveia等[16]在文獻[10]的技術基礎上求解波動方程，并聯(lián)合貝葉斯定理實現(xiàn)了反演，降低了不確定性。M. Sen等[17]在Kennett遞歸算法的基礎上，采用高斯–牛頓法求解反演的目標函數(shù)。但這些方法都是利用基本反射率法進行的，計算過程復雜，計算時間長。Liu Hongxing等[18]運用矢量化的反射率法借助高斯分布，實現(xiàn)了基于反射率法的疊前反演，提高了計算效率，但忽視了參數(shù)極值，精度有待進一步提高。陳莉[19]利用矢量化的反射率法進行正演，引入微分拉普拉斯分布，并結合貝葉斯理論，利用高斯–牛頓法求解目標函數(shù)，在保證計算效率的同時提高了反演精度。

由于地下情況的復雜性和獲取資料的不完備性，反演通常具有多解性，由于地震數(shù)據(jù)帶限性質以及反演中的病態(tài)條件，使得同一地震數(shù)據(jù)的反演存在多個不同的結果；并且反演過程是不穩(wěn)定的，若地震資料中的噪聲強烈或者其他方面的干擾因素嚴重，會導致較大的誤差。針對這些問題，基于貝葉斯理論的反演方法是一種有效的解決方法，不僅能引入先驗信息約束反演過程，還能夠提高穩(wěn)定性，保證反演效果[18-19]。因此，將貝葉斯理論引入波阻抗反演中，引入柯西分布作為先驗約束，有助于進一步提高反演結果的精度；結合反射率法，考慮地震轉換波和透射損失等影響，提高正演精度；最后給出波阻抗反演的求解算法。通過實際數(shù)據(jù)的應用，證明該方法的有效性，為煤炭資源的勘探開發(fā)提供一種行之有效的高精度高分辨率的反演方法，可以用于識別薄煤層及小構造。

1 原理與方法

1.1 反射率法正演原理

反射率法是在一維假設下求解彈性波方程的一種方法，因而能夠模擬多次波、轉換波，并考慮地震波在傳播過程中的損耗等影響。煤層與周圍地層的波阻抗差異通常較大，傳統(tǒng)方法的誤差較大，不滿足精細勘探開發(fā)的要求。反射率法則適用于波阻抗差異大的情況下的反演問題，然而考慮到傳統(tǒng)的基于遞歸矩陣的反射率法的計算過程復雜且耗時，研究基于矢量化的反射率法[18-20]展開，該方法簡化了計算過程，有效減小了運算耗時。其求解過程是在頻率–慢度域進行的，利用數(shù)值積分方法將頻率–慢度域的反射率變換到時空域或者截距時間–射線參數(shù)域。

假設地下介質在第個地層界面之下為半空間彈性介質，無反射波，因此，第層介質以下的響應可表示為：

1.2 基于貝葉斯理論的反演方法

貝葉斯理論的基本表達式為：

假設地震數(shù)據(jù)中的噪聲是相互獨立的并且服從高斯分布，則噪聲可以表示為：

為了提高反演結果的分辨率，研究將柯西分布引入到地震反演中。高斯分布由于引入了模型參數(shù)的先驗信息，從而提高了反演過程的穩(wěn)定性，但是高斯分布是一種光滑分布，對模型的極大極小值有壓制作用，從而使得反演結果具有一定的平滑效應，即降低了反演的分辨率[17]?？挛鞣植紕t是一種相對稀疏的分布，在反演過程中有利于保護模型參數(shù)中的高值，從而有效提高反演分辨率?？挛鞣植嫉恼齽t化項可以表示為：

提取后驗概率分布最大解即可獲得相應的反演波阻抗信息，即對后驗概率分布函數(shù)求偏導數(shù)，并令其等于零，求解獲得的參數(shù)就是最大后驗概率解，上式等價于求解下述目標函數(shù)的極小值：

采用普遍適用的高斯–牛頓迭代來求解上式，為了保證方法的計算效率在求解過程中忽略高階項，則高斯–牛頓的迭代更新公式可以表示為：

圖1 基于反射率法波阻抗反演流程

2 方法測試

2.1 模型數(shù)據(jù)分析

模型由六層砂泥介質組成，其中砂巖層較薄，相應的縱橫波與密度信息如圖2所示。選擇主頻為30 Hz的雷克子波，基于反射率法實現(xiàn)正演，得到該模型的道集(圖3a)，從角道集中可以看出，由于巖層厚度較薄，同相軸疊加在一起，同時可以發(fā)現(xiàn)，由于反射率法可以模擬多次波，地震記錄中存在一些能量較弱的同相軸；且在不同界面處，同相軸的能量強度并不相同(圖3d)，這是由于地震波的各種傳播效應所引起的振幅變化，在實際地層中這些影響是真實存在的。而利用Zoeppritz方程正演獲得的地震記錄如圖3b所示，可以看出在上下層介質性質相同的情況下，Zoeppritz方程正演的地震記錄在兩個薄層底界面反射同相軸能量完全相同(圖3e)，沒有考慮傳播效應對振幅引起的變化，且無多次波同相軸存在；從圖3c可以看出，兩種方法合成的地震記錄存在差異。在薄地層中，多次波反射能量常常疊加在一次波反射之上，嚴重降低了地震記錄的分辨率且難以完全去除，反演結果出現(xiàn)偏差；而透射損失也引起地震記錄的振幅變化，若不考慮該影響，則正演的地震記錄與真實的地震記錄存在差異，影響反演的精度。Zoeppritz方程忽略了波的傳播效應，不能將多次波和透射損失考慮在內，正演模擬能力有限，而反射率法可以模擬多次波，并考慮各種波的傳播效應，更接近真實的地震記錄。因此，基于反射率法的反演方法不要求對多次波進行處理，簡化處理流程，降低對地震波振幅處理的要求，正演結果更加符合實際情況。

圖2 測試模型1

2.2 實際數(shù)據(jù)應用

將提出的方法應用于布爾臺礦區(qū)三維工區(qū)，識別薄層煤和深層煤炭資源，指導進一步的勘探開發(fā)。工區(qū)內含煤地層為侏羅系中下統(tǒng)延安組，其煤系的沉積基底為三疊系上統(tǒng)延長組。地層保存完整，未遭受后期剝蝕。含煤地層由陸源碎屑巖組成，其巖性組合為各級粒度的砂巖、粉砂巖、泥巖及煤層，呈規(guī)律性交替出現(xiàn)。沉積相由河流相、湖泊三角洲相、湖相組成，為一套大型內陸盆地含煤建造，具有良好的經(jīng)濟價值和開發(fā)潛力。

工區(qū)面元大小為10 m×5 m，工區(qū)內有30口井數(shù)據(jù)(包括縱橫波速度與密度曲線)可供利用，但井深度較淺。圖4為該工區(qū)隨機抽取的某測線的地震疊加剖面，從圖4可以看出，在0.25~0.40 s位置具有強振幅反射特征，根據(jù)實際勘探開采結果，該位置處發(fā)育多層煤層，且分布廣泛，如圖中箭頭指示；而在深部地層，沒有明顯的強震幅反射特征。首先從三維地震數(shù)據(jù)中提取出地震子波，然后進行井震標定，便于測井數(shù)據(jù)與地震資料的對應，經(jīng)過井震標定后，兩類數(shù)據(jù)可以較好的匹配，便于反演過程中實現(xiàn)井的約束。在反演中，提取E井作為驗證井未參與反演過程。該測線的反演結果如圖5a所示，可以看出，反演結果能夠獲得地下介質的彈性信息，有效地識別不同的地層，分辨能力強。通過彈性信息的分布情況，可以直接推斷出在0.25~0.40 s之間分布有交替呈現(xiàn)的煤層，其厚度分布也可從圖中直觀地估算獲得，煤層阻抗值較低(圖中藍色表示阻抗值小，指示煤層)。同時，從反演結果中，可以清晰地判別在1.3~1.4 s分布有連續(xù)性較好的較厚煤層，有利于開發(fā)過程中的相關決策，而深部煤層在原始的地震疊加剖面上體現(xiàn)不明顯，表明了方法在煤層預測及煤厚識別方面的有效性。圖5b為常規(guī)方法反演獲得的波阻抗分布，雖然也能指示煤層位置，但其對于淺層薄煤層反演的效果較差，分辨率較低(如圖中箭頭指示位置)。為了更直觀地說明方法的優(yōu)勢，圖6展示了過驗證井剖面的反演結果，圖中黑色直線表示驗證E井的位置，圖6a為本文提出方法的反演結果，圖6b為常規(guī)方法的反演結果，對比發(fā)現(xiàn)，常規(guī)方法對于較薄的煤層及其夾層不能較好地反映出來，分辨率有待進一步提高，而提出方法能較好地反演出不同厚度煤層的分布情況。圖6c是驗證井與井旁道反演結果的對比，可以看出，雖然該井并未參與反演，但反演結果能與測井數(shù)據(jù)匹配，提出方法及傳統(tǒng)方法的反演結果與井數(shù)據(jù)的相關系數(shù)分別為0.84與0.68，再次說明了反射率法波阻抗法在煤田勘探中應用的可行性。

圖3 正演地震記錄(三個負極性反射是三個砂泥巖底界面的反射情況)

圖4 研究區(qū)的二維測線疊加剖面

3 結論

a. 提出的基于貝葉斯理論的反射率法波阻抗反演方法能夠模擬全波場和波的傳播效應，考慮了透射損失、多次波等對地震響應的影響，有效簡化前期的數(shù)據(jù)處理流程并提高反演精度。

b. 在貝葉斯框架下，引入柯西分布先驗信息，有效地保護了地下介質彈性參數(shù)的極值，從而提高反演的分辨率，不僅能有效識別煤層，還能估算煤層厚度，可以為頂?shù)装逅八ǖ赖念A測提供有力的方法支撐；該方法也可推廣到反演彈性阻抗等信息，在煤層的勘探開發(fā)中具有顯著的優(yōu)勢和良好的應用前景。

圖5 選取的二維測線的反演結果

圖6 不同方法反演結果

c. 該方法要求具有橫波測井數(shù)據(jù)，在未能提供橫波數(shù)據(jù)的地區(qū)可通過模式識別等方法預測橫波曲線；此外該方法未考慮橫向相關性，是單道反演方法，也是基于各向同性理論的反射率法，因此，在后續(xù)的進一步研究中，將考慮橫向相關性，并研究各向異性介質反射率法，以更真實地表達地震波的傳播特征，提升反演能力和效果。

請聽作者語音介紹創(chuàng)新技術成果等信息，歡迎與作者進行交流

[1] 黃亞平，束榮華，董守華，等. 波阻抗反演技術在煤田地震巖性勘探中的應用[J]. 煤田地質與勘探，2008，36(5)：62–64. HUANG Yaping，SHU Ronghua，DONG Shouhua，et al. Application of wave impedance inversion technique to the seismic lithological exploration in coalfield[J]. Coal Geology & Exploration，2008，36(5)：62–64.

[2] 任川，潘冬明，彭劉亞，等. 利用彈性波阻抗反演預測構造煤發(fā)育[J]. 煤田地質與勘探，2014，42(3)：92–95. REN Chuan，PAN Dongming，PENG Liuya，et al. Prediction of tectonic coal development using elastic impedance inversion[J]. Coal Geology & Exploration，2014，42(3)：92–95.

[3] 王千遙，單蕊. 波阻抗反演技術在煤田勘探中的應用[J]. 煤炭技術， 2018，37(8)：97–100. WANG Qianyao，SHAN Rui. Application of wave impedance inversion technique in coalfield exploration[J]. Coal Technology，2018，37(8)：97–100.

[4] LIU Xingye，CHEN Xiaohong，LI Jingye，et al. Facies identification based on multikernel relevance vector machine[J]. IEEE Transactions on Geoscience and Remote Sensing，2020，58(10)：1–14.

[5] 劉興業(yè)，陳小宏，李景葉，等. 基于核貝葉斯判別法的儲層物性參數(shù)預測[J]. 石油學報，2016，37(7)：878–886. LIU Xingye，CHEN Xiaohong，LI Jingye，et al. Reservoir physical property prediction based on kernel-Bayes discriminant method[J]. Acta Petrolei Sinica，2016，37(7)：878–886.

[6] 劉興業(yè)，李景葉，陳小宏，等. 聯(lián)合多點地質統(tǒng)計學與序貫高斯模擬的隨機反演方法[J]. 地球物理學報，2018，61(7)：364–373. LIU Xingye，LI Jingye，CHEN Xiaohong，et al. A stochastic inversion method integrating multi-point geostatistics and sequential Gaussian simulation[J]. Chinese Journal of Geophysics，2018，61(7)：2998–3007.

[7] FERTIG J G, MüLLER G. Computations of synthetic seismograms for coal seams with the reflectivity method[J]. Geophysical Prospecting，1978，26(4)：868–883.

[8] LIU Xingye，LI Jingye，CHEN Xiaohong，et al. Bayesian discriminant analysis of lithofacies integrate the Fisher transformation and the kernel function estimation[J]. Interpretation，2017，5(2)：SE1–SE10.

[9] LIU Xingye，LI Jingye，CHEN Xiaohong，et al. Stochastic inversion of facies and reservoir properties based on multi-point geostatistics[J]. Journal of Geophysics and Engineering，2018，15(6)：2455–2468.

[10] FUCHS K，MüLLER G. Computation of synthetic seismograms with the reflectivity method and comparison with observations[J]. Geophysical Journal International，1971，23(4)：417–433.

[11] KENNETT B L N. Theoretical reflection seismograms for elastic media[J]. Geophysical Prospecting，1979，27(2)：301–321.

[12] FRYER G J. A slowness approach to the reflectivity method of seismogram synthesis[J]. Geophysical Journal International， 1980，63(3)：747–758.

[13] MULLER G. The reflectivity method：A tutorial[J]. Journal of Geophysics，1985，58(1/2/3)：153–174.

[14] MALLICK S，F(xiàn)RAZER L N. Practical aspects of reflectivity modeling[J]. Geophysics，1987，52(10)：1355–1364.

[15] ZHAO H，BJ?RN U，AMUNDSEN L. Frequency-wave number elastic inversion of marine seismic data[J]. Geophysics，1994，59(12)：1868–1881.

[16] GOUVEIA W P， SCALES J A. Bayesian seismic waveform inversion：Parameter estimation and uncertainty analysis[J]. Journal of Geophysical Research，1998，103(B2)：2759–2779.

[17] SEN M K，ROY I G. Computation of differential seismograms and iteration adaptive regularization in prestack waveform inversion[J]. Geophysics，2003，68(6)：2026–2039.

[18] LIU Hongxing，LI Jingye，CHEN Xiaohong，et al. Amplitude variation with offset inversion using the reflectivity method[J]. Geophysics，2016，81(4)：R185–R195.

[19] 陳莉. 基于貝葉斯框架的反射率法AVO反演研究[D]. 北京：中國石油大學(北京)，2017. CHEN Li. Prestack AVO inversion based on the reflectivity method using Bayesian theory[D]. Beijing：China University of Petroleum(Beijing)，2017.

[20] PHINNEY R，ODOM R，F(xiàn)RYER G. Rapid generation of synthetic seismograms in layered media by vectorization of the algorithm[J]. Bulletin of the Seismological Society of America，1987，77(6)：2218–2226.

Seismic impedance inversion of coal field with reflectivity method based on Bayesian theory

YANG Zhen1, GAO Zhenyu1, LIU Xingye1, 2

(1. Buertai Coal Mine, Shenhua Shendong Coal Group Corporation Limited, Ordos 017000, China; 2. College of Geology and Environment, Xi’an University of Science and Technology, Xi’an 710054, China)

Seismic inversion can reflect the morphology of subsurface formation and facies information more intuitively. Seismic impedance inversion can directly produce elastic information of underground media, providing reliable information for subsequent coal identification, water channel prediction and collapse column distinguishing. A Bayesian-based inversion of reflectivity method was developed. It can consider the various propagation effects of seismic waves and solve the objective function nonlinearly, which can more accurately calculate the impedance and improve the inversion resolution. The method is applied to the actual data from Buertai coal mine and satisfactory results are output, which effectively verifies the feasibility and effectiveness of the new method. Compared with traditional inversion method, the resolution and accuracy are improved. It can effectively identify the distribution of thin coal and deep coal. It can provide valuable information for identifying the coal and the collapse column by using seismic inversion technology and predicting the water distribution in the roof and the floor.

seismic inversion; reflectivity method; Bayesian theory; wave propagation effect; coal exploration

P631.4

A

10.3969/j.issn.1001-1986.2020.03.029

1001-1986(2020)03-0204-07

2020-01-01；

2020-04-30

國家自然科學基金項目(41904116)；陜西省煤炭綠色開發(fā)地質保障重點實驗室基礎研究計劃項目(MTy2019-20)

National Natural Science Foundation of China(41904116)；Basic Research Project of Shaanxi Provincial Key Laboratory of Geological Support for Coal Green Exploitation(MTy2019-20)

楊真，1970年生，男，陜西橫山人，碩士，高級工程師，從事煤礦管理工作. E-mail：1224252292@qq.com

劉興業(yè)，1991 年生，男，山西晉中人，博士，講師，從事地震勘探工作. E-mail：lwxwyh506673@126.com

楊真，高振宇，劉興業(yè). 基于貝葉斯理論的反射率法煤田地震波阻抗反演[J]. 煤田地質與勘探，2020，48(3)：204–210.

YANG Zhen，GAO Zhenyu，LIU Xingye. Seismic impedance inversion of coal field with reflectivity method based on Bayesian theory[J]. Coal Geology & Exploration，2020，48(3)：204–210.

(責任編輯 聶愛蘭)