一种利用联合反演方法检测构造煤发育的方法与流程

本发明涉及地震勘探技术领域，具体涉及一种利用多分量地震数据模型的联合反演方法检测构造煤发育的方法。

背景技术：

煤作为一种断裂、易碎的介质，在构造应力作用下，煤很容易破碎变形。在世界各地的煤田中，构造变形煤(tdc)是导致煤矿灾难的因素之一。tdc被认为与煤炭和天然气的突出有关。通常情况下，由于较大的吸附表面，tdc比未变形的煤更容易吸收更多的气体。在中国，很大一部分煤矿是瓦斯高突出的煤矿。在隧道施工前确定tdc的开发和分布将有助于避免煤和瓦斯突出的灾难。地震方法经常用于工程安全评估。在煤田地震勘探中，从地震资料中提取与tdc发育和分布有关的信息已成为一个热点和难点的研究课题。

许多学者根据地震勘探数据来探测tdc进行了大量的研究。对于一些浅煤层，当煤区变形时，可以通过对地震叠加剖面的解释来找到相应的指示。这些预测包括相位异常、振幅变化和同相轴不连续。然而，在煤田地震学中，构造变形被认为是煤层的固有属性，这是很难直接识别出来的。因此，一些学者试图通过地震属性分析来确定煤层中发育异常的地质构造。ge等人用一种层间地震法来找出煤质中的异常现象。王等人应用极限学习机和主成分分析来预测tdc厚度。其他学者试图通过将地震数据转换成地层岩性的属性数据来确定煤块内发育的tdc。李等人只使用p波叠前反演，推导出了描述tdc的弹性阻抗、声阻抗和孔隙度等评价参数。

当目标煤层随着深度变得越来越薄时，在地震勘探中只使用p波会导致对煤层属性的解释有很大的不确定性。为了弥补传统的p波勘探的不足，煤田多分量地震勘探已经进行了试验。多组分地震勘探是可行的，可为储层预测提供可靠的p-s波速度比(α/β值)。stewart等人指出，多组分地震的目的是记录和利用反射的p和s波(pp和ps波)，以获得与岩石属性相关的更多信息。在多分量地震方法中，联合pp和ps波的振幅偏移距变化(avo)反演是一个岩石性质特征的高精度分析方法。与传统的pp波avo反演相比，在反演中加入ps波将减少多解问题。veire和landro应用最小二乘理论联合pp和ps反演。kurt提出了一种通用的线性算法，用于联合pp和ps的avo反演。du和yan应用联合pp和ps的avo反演来检测流体。lu等人在联合pp和ps的avo反演中使用了泰勒展开方法。在煤田勘探领域，lu等人论证了煤系地层预测的成功案例，其中联合pp和ps的avo反演方法发挥了重要作用。

然而，可靠的地震反演依赖于超前的岩石物理知识。近年来，对煤的测井曲线和岩石物理特性进行了深入的研究。李等人通过对测井资料的分析，对储层渗透率进行了评价。徐等人系统地研究了tdcs的测井曲线特征。对煤岩物理试验的研究主要集中在微观结构，力学性能和弹性性质。然而，对于地震反演来说，与地震波速度有关的经验关系是最重要的。王等人的岩石物理试验表明，在相似的地层条件下，tdcs的声速和体积密度低于未变形的煤层。然而，对于在不同位置收集的煤炭样品，不管p或s波的速度是很难找到确定的阈值来区分tdc和未变形的煤。p-s波速度和密度之间的关系很差，不变形和变形的煤之间的密度差并不明显。morcote等人关于动态弹性特性的测试数据表明，与单p或s波速度相比，煤的纵横波速度比(α/β)值受到压力的影响较小。因此，对于tdc检测来说，α/β值是一个有效且稳定的因素。从王等人的测试中可以发现，大部分变形的煤样的α/β值小于1.7。从陈等人的实验数据，对tdc的超声速度进行了分析，得出了相同的经验关系。该文提出了一种基于模型的多分量地震数据联合反演的tdc检测方法。基于模型的反演可以将高频岩石物理模型与地震反演结果联系起来。mallick应用基于模型的反演将高频测井信息添加到avo反演结果中。spikes和dvorkin使用基于模型的反演来推导岩性和物理参数。然而，关于基于模型的联合pp和ps的avo反演理论的研究却很少。

技术实现要素：

针对现有技术的不足，本发明旨在提供一种利用多分量地震数据模型的联合反演方法检测构造煤发育的方法，应用最小二乘法理论，把pp和ps波的avo道集联合反演成相应的α/β值，另外利用p-s波速度和密度测井产生的先验模型在联合反演中得到应用，以提高反演精度。

为了实现上述目的，本发明采用如下技术方案：

一种利用联合反演方法检测构造煤发育的方法，包括如下步骤：

s1、从采集的三分量地震数据中分离出pp波与ps波炮集数据；

s2、对pp波与ps波炮集数据进行叠前预处理；

s3、对叠前预处理之后的pp波与ps波炮集数据进行速度分析，并抽成pp波共成像点道集和ps波共成像点道集；

s4、将pp波共成像点道集与ps波共成像点道集分别转换为pp波ava道集和ps波ava道集；

s5、用p波测井曲线与pp波井旁道子波进行合成得到合成的pp波ava道集，并用合成的pp波ava道集标定p波测井曲线与步骤s4中得到的pp波ava道集；

s6、用p波与s波测井曲线与ps波井旁道子波进行合成得到合成的ps波ava道集，并用来标定p波与s波测井曲线与步骤s4中得到的ps波ava道集；

s7、根据标定后的pp波与ps波反射时间关系，将步骤s4中得到的ps波ava道集压缩至pp波零偏移距反射时间上；

s8、将测井数据沿着地震层位外推，插值形成p波速度、s波速度以及密度的初始模型；

s9、分别估计步骤s4中得到的pp波ava道集与步骤s7中压缩至pp波零偏移距反射时间上的ps波ava道集的地震子波；

s10、建立联合反演的雅各比矩阵j，并建立目标函数；

s11、高斯-牛顿反演并更新初始模型；若更新后的模型能够满意则输出模型；若更新后的模型不满意，则将该模型作为初始模型，并回到步骤s10；

s12、将更新后输出的模型作为反演模型，计算纵横波速度比模型；

s13、基于纵横波速度比模型预测构造煤发育与分布。

本发明的有益效果在于：经过模型试验结果表明，本发明方法具有较高的抗噪声能力和良好的tdcs识别能力。

附图说明

图1为当p波入射到地下弹性界面时，反射波和透射波的示意图。

图2为pp和ps子波的时间样本图。

图3为本发明实施例提供的用基于模型的pp和ps道集联合反演tdc预测的实施流程图；

图4为含煤地层的理论模型示意图；

图5为合成ava道集示意图，其中(a)为纯pp波合成ava道集，(b)为纯ps波合成ava道集，(c)为加入了10％随机噪声的pp波合成ava道集，(d)为加入

图6为用纯合成的pp波和ps波ava道集得出的基于模型的反演结果示意图，其中(a)表示p波速度(α)，(b)表示s波速度(β)，(c)表示密度(ρ)，(d)表示p、s波速度比(α/β)。

图7为加入了10％随机噪声的合成pp波和ps波道集得出的基于模型的反演结果示意图，其中(a)代表p波速度(α)，(b)代表s波速度(β)，(c)代表密度(ρ)，(d)代表p、s波速度比(α/β)。

具体实施方式

以下将结合附图对本发明作进一步的描述，需要说明的是，本实施例以本技术方案的前提，给出了详细的实施方式和具体的操作过程，但本发明的保护范围并不限于本实施例。

以下先对基于模型的联合pp和ps的avo反演理论作进一步的描述。

1、基于模型的联合反演

在多分量地震勘探中，如图1所示，p波源常用于引发地震波，并在地面上布置多分量检波器，以记录反射波。通常，多分量检波器有三个分量：一个垂直分量(z)和两个水平分量(x和y)，反射的pp和ps波是由多分量检波器在地面上同时记录的。由于p波的偏振方向与传播方向平行，所以大部分的pp波的能量都是由z分量记录的。然而，由于s波的偏振方向与传播方向正交，所以大部分的ps波能量被x和y分量接收。

对于地下弹性界面(图1)，使用α1，β1，ρ1来表示纵波速度、横波速度和上地层的密度，α2，β2，ρ2表示纵波速度、横波速度、下地层密度，pp和ps反射系数(rpp和rps)由zoeppritz方程分别给出：

其中i1表示入射(或反射)p波的角度，j1表示反射的s波角度，i2和j2分别表示透射的p波和s波角度，另外：

其中p为射线参数^[41]；而其他参数a，b，c，d，f，h，d没有具体的物理意义。它们只是用来简化rpp和rps的表达式。在反演时间窗内的所有样本中都使用了p波和s波速度以及密度参数，模型参数矩阵e可以表示为：

e＝(αβρ)^t,(4)

其中：

α＝(α1α2lαm),(5)

β＝(β1β2lβm),(6)

ρ＝(ρ1ρ2lρm).(7)

矩阵α，β，ρ是由纵波速度、横波速度和密度在时间样品1，2，…m的时间窗内联合反演。基于方程(1)(3)，考虑每个时间样本作为一个界面，分别计算pp和ps波反射系数矩阵rpp和rps。然后，合成pp波和ps波的振幅随入射角变化(ava)集合可以推出：

spp＝wpprpp,(8)

sps＝wpsrps,(9)

其中：

是在时间窗口中时间样本1、2、…m-1的pp波和ps波的反射系数。因为一个样品的反射系数的计算需要连续两个样品的弹性参数，矩阵中的元素数量rpp或rps是比矩阵α，β或者ρ要少。wpp和wps分别是pp和ps子波矩阵，它们的形式是：

其中,分别是pp和ps子波事件样品的振幅(图2)。n是子波时间样本的最大数量，应该不到m。在实际应用程序中，pp和ps子波可以通过在频率域中构造振幅和相位谱提取出来。

在提出的反演方法中，第一步是直接计算模型参数矩阵e，根据阻尼最小二乘理论，高斯牛顿公式化简如下：

e^(k)＝e^(k-1)+δe^(k),(13)

h＝j^tj,(15)

k表示迭代次数，e^(k)是第k次迭代之后最新的模型参数矩阵。δe^(k)是第k次模型的修改矩阵：

δe^(k)＝(δα^(k)δβ^(k)δρ^(k))^t.(16)

j和h分别是jacobian和hessian矩阵^[45，46]。在联合avo反演中，如果pp和ps波ava道集的入射角(在等式(1)和(2)中的i1)是θ1，θ2，…θx，则jacobian矩阵j形式如下：

在等式(12)中，λ是一个阻尼参数，i是一个单位矩阵。这两个参数用在稀疏最小二乘法问题中以确保总是存在一个唯一解。一开始给出先验模型矩阵e⁽⁰⁾，则第k次迭代之后的矩阵e^(k)和δe^(k)就可以导出来。通过控制迭代次数，建立了联合avo反演的目标函数如下：

其中和分别是输入的pp和ps波ava道集；和分别是第k次迭代之后合成的pp和ps波ava道集。权重因子ω大小从0到1。如果ω大于0.5，地震数据更多的涉及到约束，否则，先验模型对约束的控制更大。当q^(k)达到合适的最小值或者几乎不变时，迭代就可以停止了，与此相对应得e^(k)则是输出的反演模型矩阵。

基于上述基于模型的联合pp和ps的avo反演理论，本实施例提供了一种利用多分量地震数据模型的联合反演方法检测构造煤发育的方法,包括如下步骤(如图3所示)：

s1、从采集的三分量地震数据中分离出pp波与ps波炮集数据；

s2、对pp波与ps波炮集数据进行叠前预处理；

s3、对叠前预处理之后的pp波与ps波炮集数据进行速度分析，并抽成pp波共成像点道集和ps波共成像点道集；

s4、将pp波共成像点道集与ps波共成像点道集分别转换为pp波ava道集和ps波ava道集；

s5、用p波测井曲线与pp波井旁道子波进行合成得到合成的pp波ava道集，并用合成的pp波ava道集标定p波测井曲线与步骤s4中得到的pp波ava道集；

s6、用p波与s波测井曲线与ps波井旁道子波进行合成得到合成的ps波ava道集，并用来标定p波与s波测井曲线与步骤s4中得到的ps波ava道集；

s7、根据标定后的pp波与ps波反射时间关系，将步骤s4中得到的ps波ava道集压缩至pp波零偏移距反射时间上；

s8、将测井数据沿着地震层位外推，插值形成p波速度、s波速度以及密度的初始模型；

s9、分别估计步骤s4中得到的pp波ava道集与步骤s7中压缩至pp波零偏移距反射时间上的ps波ava道集的地震子波；

s10、建立联合反演的雅各比矩阵j，并建立目标函数；

s11、高斯-牛顿反演并更新初始模型；若更新后的模型能够满意则输出模型；若更新后的模型不满意，则将该模型作为初始模型，并回到步骤s10；

s12、将更新后输出的模型作为反演模型，计算纵横波速度比模型；

s13、基于纵横波速度比模型预测构造煤发育与分布。

以下通过模型试验对本实施例方法的性能作进一步的说明。

本试验中使用了一个合成数据集来证明本实施例所提出的tdc检测方法的有效性。含煤地层的理论模型如图4所示，其中3-5层为煤层。从表1中可以看到层5是tdc煤层，其弹性参数远远低于未变形的煤层(4和6层)。5层的α/β值小于周围的岩石。在pp反射时域内，分别由地震子波(主频率30hz)矩阵和pp和ps反射系数矩阵的乘积合成纯pp波ava道集(图5a)和ps波ava道集(图5b)。对应的入射角度为5°到30°，间隔为5°。然后用10％的随机噪声(图5c，d)分别加入纯ava道集。

表1：含煤地层的理论模型参数

如图6所示，虽然煤层很薄，初始模型与实际模型相差很大，但是基于模型的联合反演在无噪声条件下也能得到准确的结果。变形煤在任何反演模型上都很容易识别。将10％的随机噪声添加到合成的ava道集中，如图7所示，在反演模型中有很多异常。然而，反演出的α/β模型(图7d)仍比其他反演模型稳定(图7c)，变形煤层仍是清晰的。

本实施例提出了一种基于模型的多分量地震数据联合反演的tdc检测方法。α/β值是用来区分tdcs和未变形的煤。在反演中同时使用了声波和测井资料。因此，测井数据揭示的tdc发育的先验知识可以添加到反演结果中。模型试验结果表明，基于模型的联合反演能有效地识别出在噪声条件下的煤层内的tdc。

专业术语：

tdc构造变形煤

p-wave压缩波/纵波

s-wave剪切波/横波

pp-wave转换纵波

ps-wave转换横波

α地层中纵波速度

α1上部地层纵波速度

α2下部地层纵波速度

β地层中横波速度

β1上部地层横波速度

β2下部地层横波速度

α/β纵横波速度比

ρ地层密度

ρ1上部地层的密度

ρ2下部地层的密度

i1纵波入射角

j1横波反射角

i2纵波透射角

j2横波透射角

p射线参数

rpp转换纵波的反射系数

rps转换横波的反射系数

rpp转换纵波的反射系数矩阵

rps转换横波的反射系数矩阵

α地层的纵波速度矩阵

β地层的横波速度矩阵

ρ地层的密度矩阵

e地层的模型参数矩阵

δe地层的模型矫正矩阵

e⁽⁰⁾地层的先验模型矩阵

e^(k)地层第k次迭代之后的矫正模型矩阵

wpp转换纵波子波矩阵

wps转换横波子波矩阵

j雅克比矩阵

h海塞矩阵

输入的转换纵波振幅随角度变化道集

输入的转换横波振幅随角度变化道集

第k次迭代之后合成的转换纵波ava道集

第k次迭代之后合成的转换横波ava道集

q^(k)第k次迭代的目标函数

ω权重系数

对于本领域的技术人员来说，可以根据以上的技术方案和构思，作出各种相应的改变和变形，而所有的这些改变和变形都应该包括在本发明权利要求的保护范围之内。