基于贝叶斯理论的反射率法煤田地震波阻抗反演_杨真.pdf
第 48 卷 第 3 期 煤田地质与勘探 Vol. 48 No.3 2020 年 6 月 COAL GEOLOGY 2. College of Geol- ogy and Environment, Xi’an University of Science and Technology, Xi’an 710054, China Abstract Seismic inversion can reflect the morphology of subsurface ation and facies ination more intui- tively. Seismic impedance inversion can directly produce elastic ination of underground media, providing re- liable ination for subsequent coal identification, water channel prediction and collapse column distinguishing. A Bayesian-based inversion of reflectivity 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 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 . Compared with traditional inversion , the resolution and accuracy are improved. It can effectively identify the distribution of thin coal and deep coal. It can provide valuable ination for identifying the coal and the collapse column by using seismic inversion technology and predicting the water distribution in the roof and the floor. Keywords seismic inversion; reflectivity ; Bayesian theory; wave propagation effect; coal exploration 地震波阻抗反演是将地震资料转换为波阻抗剖 面的一种技术,它将地震资料、测井数据、地质解 释相结合,充分利用了测井资料具有较高的垂向分 辨率和地震资料具有较高的横向分辨率的特点,反 ChaoXing 第 3 期 杨真等 基于贝叶斯理论的反射率法煤田地震波阻抗反演 205 演获得的波阻抗剖面,不仅便于解释人员将地震资 料与测井资料连接对比,而且对地层相关参数的变 化研究提供可靠资料,从而指导资源的勘探开发。 煤炭资源的勘探开发离不开地震勘探技术,其 中波阻抗反演可以预测煤层位置,指示煤系气分布 情况,分析煤层厚度的变化,估算资源储量,查明 陷落柱发育情况,同时对富水区也能有显著的指示 作用, 并对后续的进一步研究提供直观可靠的资料[1-4]。 因而,在煤炭及其相关资源的勘探开发中越来越受 到重视。但是,目前的波阻抗反演算法都是基于简 单的褶积模型进行的,通常需要一定的假设条件。 当不满足这些假设条件时,严重制约反演精度,同 时地震资料的带限性质会影响反演分辨率,从而影 响勘探效果,增加勘探开发风险[5-6]。此外,在处理 薄煤层问题时,由于非均质性强,转换波及透射波 损失严重,这些方法没有考虑地震波的损失,从而 制约反演精度。 弹性波动方程是一种较为精确地表达地震波传 播的方法,但受计算效率、反演稳定性和对地震资 料要求高等因素的影响,实际应用面临诸多挑战。 而反射率法是一种以一维递推方式求解弹性波动方 程解析解的方法,能对全波场模拟。除了一次反射, 反射率法还能够模拟透射波、 多种转换波和多次波, 考虑了地震波相位变化以及透射损失的影响,且反 射率法计算精度较高,计算成本适中,具有较强的 实用价值[7-9]。 反射率法由 K. Fuchs 等[10]提出,经过不断的完 善[11-14], 已经较为成熟并被成功应用于地震反演中。 Zhao H.等[15]基于Kennett递归矩阵, 通过Levenberg- Marquardt 算法实现了频率波数反演。W. P. Gouveia 等[16]在文献[10]的技术基础上求解波动方程,并联 合贝叶斯定理实现了反演, 降低了不确定性。 M. Sen 等[17]在 Kennett 递归算法的基础上,采用高斯–牛顿 法求解反演的目标函数。 但这些方法都是利用基本 反射率法进行的,计算过程复杂,计算时间长。Liu Hongxing 等[18]运用矢量化的反射率法借助高斯分 布,实现了基于反射率法的叠前反演,提高了计算 效率,但忽视了参数极值,精度有待进一步提高。 陈莉[19]利用矢量化的反射率法进行正演, 引入微分 拉普拉斯分布,并结合贝叶斯理论,利用高斯–牛 顿法求解目标函数, 在保证计算效率的同时提高了 反演精度。 由于地下情况的复杂性和获取资料的不完备 性,反演通常具有多解性,由于地震数据带限性质 以及反演中的病态条件,使得同一地震数据的反演 存在多个不同的结果;并且反演过程是不稳定的, 若地震资料中的噪声强烈或者其他方面的干扰因素 严重,会导致较大的误差。针对这些问题,基于贝 叶斯理论的反演方法是一种有效的解决方法,不仅 能引入先验信息约束反演过程, 还能够提高稳定性, 保证反演效果[18-19]。因此,将贝叶斯理论引入波阻 抗反演中,引入柯西分布作为先验约束,有助于进 一步提高反演结果的精度;结合反射率法,考虑地 震转换波和透射损失等影响,提高正演精度;最后 给出波阻抗反演的求解算法。 通过实际数据的应用, 证明该方法的有效性,为煤炭资源的勘探开发提供 一种行之有效的高精度高分辨率的反演方法,可以 用于识别薄煤层及小构造。 1 原理与方法 1.1 反射率法正演原理 反射率法是在一维假设下求解弹性波方程的一 种方法,因而能够模拟多次波、转换波,并考虑地 震波在传播过程中的损耗等影响。煤层与周围地层 的波阻抗差异通常较大,传统方法的误差较大,不 满足精细勘探开发的要求。反射率法则适用于波阻 抗差异大的情况下的反演问题,然而考虑到传统的 基于递归矩阵的反射率法的计算过程复杂且耗时, 研究基于矢量化的反射率法[18-20]展开,该方法简化 了计算过程,有效减小了运算耗时。其求解过程是 在频率–慢度域进行的,利用数值积分方法将频率– 慢度域的反射率变换到时空域或者截距时间–射线 参数域。 根据反射率法,假设在水平层状介质的第一层 顶部放置一个炮点和多个检波器, 则频率–慢度p 域的总反射系数 , pR可以通过频率–慢度域中的 六元矢量 [1 ,,6 ]www获得[18] pp00 ps00 ,41 ,51 p p Rww Rww 1 式中 0 w为地面接收到的总反射响应;R 为反射 系数,R 的两个下标 pp 分别表示入射波类型和反 射波类型,ps 同; n w为地下第n层界面以下的总 反射响应。 T spsspppsn RRRRR w 2 式中n1,,N;为缩放因子;为计算行 列式。 假设地下介质在第 N 个地层界面之下为半空间 弹性介质,无反射波,因此,第 N 层介质以下的响 应可表示为 ChaoXing 206 煤田地质与勘探 第 48 卷 T100000 N w 3 反射率法是从底界面开始通过传递矩阵 n Q递 归计算获得顶界面之下的总反射响应,从而计算得 到频率–慢度域的反射系数, 然后对反射系数与子波 的乘积进行反变换得到–p域的合成地震记录[20], 递归过程可以表示为 0011NN wQ QQw 4 1, nnnnnnn wQ wQT E T 5 研究引入变量 n E、 n T和 n T来表示 n Q,其中 矩阵 n E表达了地震波经过第n层介质时其相位的 变化, n T和 n T表述了第n层介质对地震波振幅的 影响,包含 16 个独立分量,其表达式分别为 i1/1/i1/1/i1/1/i1/1/ pspspsps diag e1ee1e hvvhvvhvvhvv n E 6 111213131511 21212121 313233333531 314233334531 51515151 616263636561 00 00 n tttttt tttt tttttt tttttt tttt tttttt T 7 615131312111 65453515 635133331213 635133332113 62423212 615131312111 00 00 n tttttt tttt tttttt tttttt tttt tttttt T 8 式 中 p v、 s v分 别 表 示 纵 横 波 速 度 ; p q 22 p 1/vp、 22 ss 1/qvp分别表示纵横波的垂 直慢度, p sinpv表示水平慢度。 令 222 ss 21/,pvv,h为第 n 层介质厚 度,则 16 个独立分量分别表示为 2 11p s12p 2 13p s15s 2 21s31p s 2 32p33p s 35s42p 22 45s51p 22 61p s62p 22 63p s /2/ /2/ i/i2 4ii2 2i2i 4ii/ 44 4 tpq qtpq tpq qtpq tqtpq q tp qtpq q tqtq tp qtq tp q qtpq tp q qt ,, ,, ,, ,, ,, ,, ,, , 65s 4pq 9 通过递归方法根据式4、式5计算得到 0 w, 代入式1即可得到频率–慢度域的反射率 , pR, 对其与地震子波乘积后反变换得到–p的合成地 震记录 i 1 , , ed 2π t t pSRp d 10 由于波阻抗反演中使用的地震记录是叠加剖 面,因此,先按照斯奈尔定律 p sinpv将频率– 慢度域的反射系数 , pR重采样,获得频率–角度 域的反射系数 , R,与地震子波乘积后,进行傅 里叶反变换,经过叠加得到合成地震记录 i 1 d ,d 2π t t SRe 11 1.2 基于贝叶斯理论的反演方法 贝叶斯理论的基本表达式为 ||PPPm dd mm 12 式中|P m d为后验概率分布;|P d m为似然 函数; 0 exp m PPmF m为先验分布, 其中 0m P 为常数项。 F m是反演正则化约束 2 F m Q mu m 13 式中m为模型参数,假设有 m N个模型参数待反 演; Q为正则化因子;u为模型背景。 假设地震数据中的噪声是相互独立的并且服从 高斯分布,则噪声可以表示为 T1 0d exp0.5PPCnn 14 式中n为噪声; 2 dd CI为噪声的协方差矩阵, 其中 2 d 为噪声方差;I为单位矩阵。 联合dG mn和式13可得到如下的似然 函数 T 1 0d |exp0.5PP d mdG mCdG m15 其中,dG mn,表示观测获得的实际地震 数据;G为正演算子,在文中表示反射率法,G m 为合成地震剖面; 0 P可用式16表示。 1/2/2 0d 12π Nd P C 16 为了提高反演结果的分辨率,研究将柯西分布 引入到地震反演中。高斯分布由于引入了模型参数 的先验信息,从而提高了反演过程的稳定性,但是 高斯分布是一种光滑分布,对模型的极大极小值有 压制作用, 从而使得反演结果具有一定的平滑效应, 即降低了反演的分辨率[17]。柯西分布则是一种相对 ChaoXing 第3期 杨真等 基于贝叶斯理论的反射率法煤田地震波阻抗反演 207 稀疏的分布,在反演过程中有利于保护模型参数中 的高值,从而有效提高反演分辨率。柯西分布的正 则化项可以表示为 T 1 2ln 1 Nm i i F mmuΦ mu 17 式中 T1 ii i ΦDΣD;矩阵Σ为模型参数的协方 差矩阵;矩阵 i D表示分选矩阵。 将先验分布Pm与似然函数|Pd m代入式12 中即可获得后验概率分布,因此,可以表示为 T 2 d 1 |exp 2 P m ddG mdG mF m 18 提取后验概率分布最大解即可获得相应的反演 波阻抗信息,即对后验概率分布函数求偏导数,并 令其等于零, 求解获得的参数就是最大后验概率解, 上式等价于求解下述目标函数的极小值 T 0 2 d 1 2 SmdG mdG mF m d 19 采用普遍适用的高斯–牛顿迭代来求解上式, 为 了保证方法的计算效率在求解过程中忽略高阶项, 则高斯–牛顿的迭代更新公式可以表示为 1kkk mmm 20 式中为反演迭代步长;为 0 Sm对模型参数的 一阶偏导。 T 20 d SG m mG mdQ mu mm 21 假设子波矩阵为S,反射系数为R,则合成地 震记录可表示为GSR,从而根据反射率法,可 得其偏导数为 S G m R mm 22 其 中 , i ** , , 1 ed 2π ij it nn p t p S R R mm , 00 00 ** *2 0 41 14 , 1 i nn n p ww ww Rmm m w , 0 011 ** n NN nn wQ Q QQw mw , **** nnnn nnnnnn nnnn QTET E TTTT E mmmm 。 获得m后,代入式21更新模型参数,直到满 足终止条件输出反演结果。综上所述,基于贝叶斯理 论的反射率法波阻抗反演的流程如图1所示。 图 1 基于反射率法波阻抗反演流程 Fig.1 Flowchart of the proposed 2 方法测试 2.1 模型数据分析 模型由六层砂泥介质组成,其中砂岩层较薄, 相应的纵横波与密度信息如图2所示。选择主频 为30 Hz的雷克子波,基于反射率法实现正演, 得到该模型的道集图3a,从角道集中可以看出, 由于岩层厚度较薄,同相轴叠加在一起,同时可 以发现,由于反射率法可以模拟多次波,地震记 录中存在一些能量较弱的同相轴;且在不同界面 处,同相轴的能量强度并不相同图3d,这是由 于地震波的各种传播效应所引起的振幅变化,在 实 际 地 层 中 这 些 影 响 是 真 实 存 在 的 。 而 利 用 Zoeppritz方程正演获得的地震记录如图3b所示, 可 以 看 出 在 上 下 层 介 质 性 质 相 同 的 情 况 下 , Zoeppritz方程正演的地震记录在两个薄层底界面 反射同相轴能量完全相同图3e,没有考虑传播 效应对振幅引起的变化,且无多次波同相轴存在; 从图3c可以看出,两种方法合成的地震记录存在 差异。在薄地层中,多次波反射能量常常叠加在 一次波反射之上,严重降低了地震记录的分辨率 且难以完全去除,反演结果出现偏差;而透射损 失也引起地震记录的振幅变化,若不考虑该影响, 则正演的地震记录与真实的地震记录存在差异, 影响反演的精度。Zoeppritz方程忽略了波的传播 效应,不能将多次波和透射损失考虑在内,正演 模拟能力有限,而反射率法可以模拟多次波,并 考虑各种波的传播效应,更接近真实的地震记录。 ChaoXing 208 煤田地质与勘探 第48卷 因此,基于反射率法的反演方法不要求对多次波 进行处理,简化处理流程,降低对地震波振幅处 理的要求,正演结果更加符合实际情况。 图 2 测试模型 1 Fig.2 The test model 1 2.2 实际数据应用 将提出的方法应用于布尔台矿区三维工区,识 别薄层煤和深层煤炭资源, 指导进一步的勘探开发。 工区内含煤地层为侏罗系中下统延安组,其煤系的 沉积基底为三叠系上统延长组。地层保存完整,未 遭受后期剥蚀。含煤地层由陆源碎屑岩组成,其岩 性组合为各级粒度的砂岩、粉砂岩、泥岩及煤层, 呈规律性交替出现。沉积相由河流相、湖泊三角洲 相、湖相组成,为一套大型内陆盆地含煤建造,具 有良好的经济价值和开发潜力。 工区面元大小为10 m5 m,工区内有30口井 数据包括纵横波速度与密度曲线可供利用,但井 深度较浅。图4为该工区随机抽取的某测线的地震 叠加剖面,从图4可以看出,在0.250.40 s位置具 有强振幅反射特征,根据实际勘探开采结果,该位 置处发育多层煤层,且分布广泛,如图中箭头指示; 而在深部地层,没有明显的强震幅反射特征。首先 从三维地震数据中提取出地震子波,然后进行井震 标定,便于测井数据与地震资料的对应,经过井震 标定后,两类数据可以较好的匹配,便于反演过程 中实现井的约束。在反演中,提取E井作为验证井 未参与反演过程。 该测线的反演结果如图5a所示, 可以看出, 反演结果能够获得地下介质的弹性信息, 图 3 正演地震记录三个负极性反射是三个砂泥岩底界面的反射情况 Fig.3 Forward gathers ChaoXing 第3期 杨真等 基于贝叶斯理论的反射率法煤田地震波阻抗反演 209 图 4 研究区的二维测线叠加剖面 Fig.4 Stacked section of 2D survey line in the study area 有效地识别不同的地层,分辨能力强。通过弹性信 息的分布情况,可以直接推断出在0.250.40 s之间 分布有交替呈现的煤层,其厚度分布也可从图中直 观地估算获得,煤层阻抗值较低图中蓝色表示阻抗 值小,指示煤层。同时,从反演结果中,可以清晰 地判别在1.31.4 s分布有连续性较好的较厚煤层, 有利于开发过程中的相关决策,而深部煤层在原始 的地震叠加剖面上体现不明显,表明了方法在煤层 预测及煤厚识别方面的有效性。 图5b为常规方法反 演获得的波阻抗分布,虽然也能指示煤层位置,但 其对于浅层薄煤层反演的效果较差,分辨率较低 如图中箭头指示位置。 为了更直观地说明方法的 优势,图6展示了过验证井剖面的反演结果,图 中黑色直线表示验证E井的位置,图6a为本文提 出方法的反演结果,图6b为常规方法的反演结果, 对比发现,常规方法对于较薄的煤层及其夹层不能 较好地反映出来,分辨率有待进一步提高,而提出 方法能较好地反演出不同厚度煤层的分布情况。 图6c是验证井与井旁道反演结果的对比, 可以看出, 虽然该井并未参与反演, 但反演结果能与测井数据匹 配, 提出方法及传统方法的反演结果与井数据的相关 系数分别为0.84与0.68,再次说明了反射率法波阻 抗法在煤田勘探中应用的可行性。 3 结 论 a. 提出的基于贝叶斯理论的反射率法波阻抗 反演方法能够模拟全波场和波的传播效应,考虑了 透射损失、多次波等对地震响应的影响,有效简化 前期的数据处理流程并提高反演精度。 b. 在贝叶斯框架下,引入柯西分布先验信息, 有效地保护了地下介质弹性参数的极值,从而提高 反演的分辨率,不仅能有效识别煤层,还能估算煤 层厚度,可以为顶底板水及水通道的预测提供有力 的方法支撑;该方法也可推广到反演弹性阻抗等信 息,在煤层的勘探开发中具有显著的优势和良好的 应用前景。 图 5 选取的二维测线的反演结果 Fig.5 The inversion result of the corresponding line ChaoXing 210 煤田地质与勘探 第48卷 图 6 不同方法反演结果 Fig.6 Inversion results by different c. 该方法要求具有横波测井数据,在未能提供 横波数据的地区可通过模式识别等方法预测横波曲 线;此外该方法未考虑横向相关性,是单道反演方 法,也是基于各向同性理论的反射率法,因此,在 后续的进一步研究中,将考虑横向相关性,并研究 各向异性介质反射率法,以更真实地表达地震波的 传播特征,提升反演能力和效果。 请听作者语音介绍创新技术成果 等信息,欢迎与作者进行交流 参考文献References OSID 码 [1] 黄亚平,束荣华,董守华,等. 波阻抗反演技术在煤田地震岩 性勘探中的应用[J]. 煤田地质与勘探,2008,36562–64. HUANG Yaping,SHU Ronghua,DONG Shouhua,et al. Application of wave impedance inversion technique to the seis- mic lithological exploration in coalfield[J]. Coal Geology Ex- ploration,2008,36562–64. [2] 任川,潘冬明,彭刘亚,等. 利用弹性波阻抗反演预测构造煤 发育[J]. 煤田地质与勘探,2014,42392–95. REN Chuan,PAN Dongming,PENG Liuya,et al. Prediction of tectonic coal development using elastic impedance inversion[J]. Coal Geology Exploration,2014,42392–95. [3] 王千遥, 单蕊. 波阻抗反演技术在煤田勘探中的应用[J]. 煤炭 技术, 2018,37897–100. WANG Qianyao,SHAN Rui. Application of wave impedance inversion technique in coalfield exploration[J]. Coal Technol- ogy,2018,37897–100. [4] LIU Xingye,CHEN Xiaohong,LI Jingye,et al. Facies iden- tification based on multikernel relevance vector machine[J]. IEEE Transactions on Geoscience and Remote Sensing,2020, 58101–14. [5] 刘兴业,陈小宏,李景叶,等. 基于核贝叶斯判别法的储层物 性参数预测[J]. 石油学报,2016,377878–886. LIU Xingye,CHEN Xiaohong,LI Jingye,et al. Reservoir physical property prediction based on kernel-Bayes discriminant [J]. Acta Petrolei Sinica,2016,377878–886. [6] 刘兴业,李景叶,陈小宏,等. 联合多点地质统计学与序贯高 斯模拟的随机反演方法[J]. 地球物理学报,2018,617 364–373. LIU Xingye,LI Jingye,CHEN Xiaohong,et al. A stochastic inversion integrating multi-point geostatistics and se- quential Gaussian simulation[J]. Chinese Journal of Geophysics, 2018,6172998–3007. [7] FERTIG J G, MLLER G. Computations of synthetic seismo- grams for coal seams with the reflectivity [J]. Geophysi- cal Prospecting,1978,264868–883. [8] LIU Xingye,LI Jingye,CHEN Xiaohong,et al. Bayesian discriminant analysis of lithofacies integrate the Fisher transfor- mation and the kernel function estimation[J]. Interpretation, 2017,52SE1–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, 1562455–2468. [10] FUCHS K,MLLER G. Computation of synthetic seismo- grams with the reflectivity and comparison with ob- servations[J]. Geophysical Journal International, 1971, 234417–433. [11] KENNETT B L N. Theoretical reflection seismograms for elastic media[J]. Geophysical Prospecting,1979,272301–321. 下转第 218 页 ChaoXing 218 煤田地质与勘探 第 48 卷 version of on ground GPR data[J]. Geophysics,2012,776 79–91. [26] 张彬. 探地雷达中的逆时偏移及速度估计[D]. 长沙中南大 学,2010. ZHANG Bin. Reverse time migration and velocity estimation of ground penetrating radar[D]. ChangshaCentral South Univer- sity,2010. [27] BRADFORD J H. Measuring water content heterogeneity using multifold GPR with reflection tomography[J]. Vadose Zone Journal,2008,71184–193. [28] 薛桂霞,邓世坤,刘秀娟. 逆时偏移在探地雷达信号处理中的 应用[J]. 煤田地质与勘探,2004,32155–57. XUE Guixia,DENG Shikun,LIU Xiujuan. An application of reverse-time migration in the ground-penetrating radar data processing[J]. Coal Geology Exploration, 2004, 321 55–57. [29] PARSEKIAN A D,SLATER L,NTARLAGIANNIS D,et al. Uncertainty in peat volume and soil carbon estimated using ground-penetrating radar and probing[J]. Soil Science Society of America Journal,2012,7651911–1918. [30] 高妍. 几种速度分析方法的对比分析与应用[J]. 工程地球物 理学报,2014,114436–440. GAO Yan. Comparison of several velocity s and applica- tion[J]. Chinese Journal of Engineering Geophysics,2014, 114436–440. [31] TOLDI J. Velocity analysis without picking[J]. Geophysics, 1989,542191–199. [32] ANNAN A P. Ground Penetrating Radar,in near surface geophysics[M]. TulsaSociety of Exploration Geophysicists, 2005. [33] HUSMAN J A, SPERL C, BOUTEN W, et al. Soil water content measurements at different scalesAccuracy of time domain re- flectometry and ground-penetrating radar[J]. Journal of Hydrol- ogy,2001,2451/2/3/448–58. [34] BERTEUSSEN K A, URSIN B. Approximate computation of the acoustic impedance from seismic data[J]. Geophysics,1983, 48101351–1358. [35] ZENG Zhaofa,CHEN Xiong,LI Jing,et al. Recursive im- pedance inversion of ground-penetrating radar data in stochastic media[J]. Applied Geophysics,2015,124615–625. [36] 冯德山, 戴前伟, 何继善. 探地雷达的正演模拟及有限差分波 动方程偏移处理[J]. 中南大学学报自然科学版, 2006, 372 361–365. FENG Deshan, DAI Qianwei, HE Jishan. Forward simulation of ground penetrating radar and its finite difference wave equation migration processing[J]. Journal of Central South Uni- versityScience and Technology,2006,372361–365. [37] 冯德山,戴前伟. 探地雷达时域多分辨法MRTD三维正演模 拟[J]. 地球物理学进展,2008,2351621–1625. FENG Deshan,DAI Qianwei. Application of the multi- resolution time domain in three dimensional forward simulation of ground penetrating radar[J]. Progress in Geo- physics,2008,2351621–1625. 责任编辑 聂爱兰 上接第 210 页 [12] FRYER G J. A slowness approach to the reflectivity of seismogram synthesis[J]. Geophysical Journal International, 1980,633747–758. [13] MULLER G. The reflectivity A tutorial[J]. Journal of Geophysics,1985,581/2/3153–174. [14] MALLICK S,FRAZER L N. Practical aspects of reflectivity modeling[J]. Geophysics,1987,52101355–1364. [15] ZHAO H, BJRN U, AMUNDSEN L. Frequency-wave number