一种适用于煤层底板复杂陷落柱的全波形反演方法

1.本发明涉及一种适用于煤层底板复杂陷落柱的全波形反演方法，属于地震勘探技术领域。


背景技术：

2.陷落柱是煤系下伏灰岩中溶洞顶部覆岩发生塌陷所形成的陷落体,是我国华北型煤田中一种典型的隐伏垂向构造。因具有随机性及不确定性，其诱发的突水隐蔽性强、速度快、水量大、破坏性强，对煤矿安全生产的危害极大。地震三维勘探在绿色煤炭资源丰富的中西部地区受地表地形以及浅部上组煤的开采和多层采空区及地表沉陷的影响，陷落柱的探测精度受到了很大限制，矿井地震勘探具有距离目标体近、地震波能量和高频成分衰减少、地震波传播路径干扰少、分辨率高等特点。因此，开展井下针对陷落柱探测具有重要意义。
3.目前工作面回采过程中针对煤层陷落柱的探测以槽波最为典型。当陷落柱发育至煤层时，会破坏煤层的连续性，利用槽波的地震波属性进行探测，典型的探测方式是利用工作面其中一条巷道布置检波器，另外一条巷道进行震源激发进行ct探测。自evison在新西兰首次于煤层中激发并且记录到了槽波，随后众多学者开展了大量研究，目前，槽波勘探已成为煤层工作面回采面内构造探测的常规工作。
4.在槽波发育的地区，槽波勘探作为主要的煤层勘探方法，但是，若煤层相对较薄或者煤层槽波不发育时，槽波则不能作为探测手段，此时发现沿着顶底板传播的折射波具有陷落柱探测可能性。因此，当单一类型的地震波成像精度受限，开展全波场信息的综合使用可以提升成像分辨率。为了获得煤田高精度的速度模型，全波形反演方法在煤田中的应用已经有了一些初探。有学者将全波形反演方法应用到煤田多参数反演中，并且全波形反演结果与测井曲线的吻合度比较高。还有一些学者在煤矿井下针对煤层内地质构造高精度成像问题，发展了槽波全波形反演方法，相比于常规的槽波成像方法，利用全波形反演可以获得更高精度的反演结果。但是目前研究仅限于合成数据的参数反演，同时现场陷落柱条件复杂，目前的全波形反演方法无法对条件复杂多个陷落柱的精确反演，因此如何能对条件复杂的陷落柱进行精确的全波形反演是亟需解决的问题。


技术实现要素：

5.针对上述现有技术存在的问题，本发明提供一种适用于煤层底板复杂陷落柱的全波形反演方法，能对煤层底板下方条件复杂的陷落柱进行精确定位及形态探测，从而为后续开采提供数据支撑。
6.为了实现上述目的，本发明采用的技术方案是：一种适用于煤层底板复杂陷落柱的全波形反演方法，该方法的具体步骤为：
7.a、在煤层工作面一侧的巷道壁相互等间距布设多个震源；在煤层工作面另一侧的巷道壁相互等间距布设多个检波器，多个检波器均通过连接主线与采集基站及采集主机连
接，使检波器、采集基站及采集主机形成地震观测系统；激发多个震源，通过地震观测系统接收观测数据并获得各个地震子波；
8.b、先构建基于卷积型的目标函数，具体定义如下：
[0009][0010]
式中，d是观测数据，u是合成数据，x
r
是检波器的位置参数，*是时间卷积算子，x
ref
表示提取参考道的位置参数；表示欧式范数，m＝(v
p
，v
s
)为速度场模型参数，密度为常数；
[0011]
其中合成数据中的正传波场通过弹性波动方程合成，具体公式如下：
[0012]
l[m(x)]u(x,t；m)＝s(x
s
,t)
ꢀꢀ
(2)
[0013]
其中，x
s
是震源坐标，l[]是弹性波正演模拟算子，u(x,t；m)是正传波场，s(x
s
,t)为震源；
[0014]
然后将步骤a获得地震子波及震源坐标代入公式(2)计算得出震源的正传波场；
[0015]
c、利用拉格朗日伴随方法获得模型参数的梯度，弹性波动方程的伴随方程为：
[0016][0017]
其中表示伴随算子，u
′
(x,t；m)表示反转波场，伴随源s
′
(x
r
,t)为：
[0018][0019]
其中表示互相关运算；
[0020]
根据公式(4)计算得出基于卷积型目标函数对应的伴随源，然后将求取的伴随源代入公式(3)中进行数据模拟，最终获得伴随源的反转波场；
[0021]
d、先确定速度场模型参数m＝(v
p
，v
s
)的梯度表达式如下：
[0022][0023]
其中表示正传波场u(x,t；m)的应力波场变量，表示表示反传波场u
′
(x,t；m)的应力波场变量；v
p
为纵波速度，v
s
为横波速度；
[0024]
将步骤b计算出的正传波场和步骤c计算出的反传波场代入到公式(5)中，最终计算得到速度场模型参数的梯度；
[0025]
e、采用共轭梯度方法求取更新方向，根据得到的梯度对速度场模型参数进行迭代更新，直至数据残差达到设定阈值或迭代更新次数达到预设值，停止迭代更新确定速度场模型参数，此时采用该参数的速度场模型进行全波形反演。
[0026]
进一步，所述步骤e的具体过程为：
[0027]
采用共轭梯度方法求取更新方向δm
k+1
，
[0028][0029]
其中下标k表示迭代次数，s
k
是用于计算共轭方向的一个中间变量，上标t表示矩阵转置，更新后的速度场模型为：
[0030][0031]
其中，t
k
表示选取的更新步长，t
k
δm
k+1
即为此次迭代过程中速度场模型参数的更新量；
[0032]
通过对速度场模型参数进行迭代更新，直至数据残差达到设定阈值或迭代更新次数达到预设值，停止迭代更新确定速度场模型参数，此时采用该参数的速度场模型进行全波形反演。
[0033]
进一步，所述步骤a中在煤层工作面的切眼处等间距设有一排检波器，使其与煤层工作面另一侧的巷道壁内设置的检波器形成l型观测系统。在切眼处增设检波器能提高观测系统采集地震数据的精确性，从而使后续反演的精度更高。
[0034]
与现有技术相比，本发明先获取震源的观测数据，然后经过计算获得震源的正传波场和伴随源的反传波场，然后确定速度场模型的梯度表达式，接着根据正传波场和反传波场及梯度计算表达式，计算得到速度场模型参数的梯度，有利于获得准确的速度场模型参数的更新方向，最终通过对速度场模型参数进行迭代更新，直至数据残差达到设定阈值或迭代更新次数达到预设值，停止迭代更新确定速度场模型参数，此时采用该参数的速度场模型进行全波形反演。因此本发明能对煤层底板下方条件复杂的陷落柱进行精确定位及形态探测，从而为后续开采提供数据支撑。
附图说明
[0035]
图1是某煤矿工作面已知的陷落柱位置及形态示意图；
[0036]
图2是采用本发明实施例对图1煤矿工作面进行纵波反演成像图；
[0037]
图3是采用本发明实施例对图1煤矿工作面进行横波反演成像图；
[0038]
图4是采用现有的纵波层析反演对图1煤矿工作面进行反演成像图；
[0039]
图5是本发明实施例中l型地震观测系统的位置示意图；
[0040]
其中，星号表示震源；三角形表示检波器。
具体实施方式
[0041]
下面将结合本文实施例中的附图，对本文实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本文一部分实施例，而不是全部的实施例。基于本文中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文保护的范围。
[0042]
本发明实施例的反演方法，具体步骤为：
[0043]
a、在煤层工作面一侧的巷道壁相互等间距布设多个震源；在煤层工作面另一侧的
巷道壁相互等间距5m布设多个检波器，在煤层工作面的切眼处等间距5m设有一排检波器，各个检波器均通过连接主线与采集基站及采集主机连接，如图5所示，使检波器、采集基站及采集主机形成l型地震观测系统；激发多个震源，通过地震观测系统接收观测数据并获得各个地震子波；
[0044]
b、先构建基于卷积型的目标函数，具体定义如下：
[0045][0046]
式中，d是观测数据，u是合成数据，x
r
是检波器的位置参数，*是时间卷积算子，x
ref
表示提取参考道的位置参数；表示欧式范数，m＝(v
p
，v
s
)为速度场模型参数，密度为常数；
[0047]
其中合成数据中的正传波场通过弹性波动方程合成，具体公式如下：
[0048]
l[m(x)]u(x,t；m)＝s(x
s
,t)
ꢀꢀ
(2)
[0049]
其中，x
s
是震源坐标，l[]是弹性波正演模拟算子，u(x,t；m)是正传波场，s(x
s
,t)为震源；
[0050]
然后将步骤a获得地震子波及震源坐标代入公式(2)计算得出震源的正传波场；
[0051]
c、利用拉格朗日伴随方法获得模型参数的梯度，弹性波动方程的伴随方程为：
[0052][0053]
其中表示伴随算子，u
′
(x,t；m)表示反转波场，伴随源s
′
(x
r
,t)为：
[0054][0055]
其中表示互相关运算；
[0056]
根据公式(4)计算得出基于卷积型目标函数对应的伴随源，然后将求取的伴随源代入公式(3)中进行数据模拟，最终获得伴随源的反转波场；
[0057]
d、先确定速度场模型参数m＝(v
p
，v
s
)的梯度表达式如下：
[0058][0059]
其中表示正传波场u(x,t；m)的应力波场变量，表示表示反传波场u
′
(x,t；m)的应力波场变量；v
p
为纵波速度，v
s
为横波速度；
[0060]
将步骤b计算出的正传波场和步骤c计算出的反传波场代入到公式(5)中，最终计算得到速度场模型参数的梯度；
[0061]
e、采用共轭梯度方法求取更新方向δm
k+1
，
[0062][0063]
其中下标k表示迭代次数，s
k
是用于计算共轭方向的一个中间变量，上标t表示矩阵转置，更新后的速度场模型为：
[0064][0065]
其中，t
k
表示选取的更新步长，t
k
δm
k+1
即为此次迭代过程中速度场模型参数的更新量；
[0066]
通过对速度场模型参数进行迭代更新，直至数据残差达到设定阈值或迭代更新次数达到预设值，停止迭代更新确定速度场模型参数，此时采用该参数的速度场模型进行全波形反演。
[0067]
效果验证：
[0068]
在某煤矿其工作面走向长度为500m，切眼长200m，如图1所示，通过现有的多种方法先确定其工作面内存在五个大小形态不一的陷落柱，陷落柱的位置及速度场模型参数如表1所示，五个陷落柱的尺寸差异较大，最小的陷落柱长轴约15m，短轴约10m，最大的陷落柱长轴约56m，短轴约26.6m。此外，陷落柱呈现不同的产状，长轴与水平方向的夹角在
‑
37.5
°
到37
°
之间变化。由于位于同一工作面的陷落柱波速特征差异不大，确定五个陷落柱的纵波速度为2600m/s，横波速度为1501m/s，底板围岩的纵波速度为3200m/s，横波速度为1847m/s；
[0069]
表1
[0070][0071]
然后在相同的观测数据和观测系统条件下，分别采用本发明实施例和现有的纵波层析反演方法对该工作面进行反演成像，其中本发明实施例获得的反演成像如图2和3所示，现有的纵波层析反演成像如图4所示；
[0072]
将两者成像图分别与图1已知的现有位置及形态示意图进行比较可知，现有的层析成像结果分辨率明显低，并且无法准确地指示出尺寸陷落柱的形态以及物性参数，同时，对小尺度的陷落柱，不能有效地指示。而采用本发明实施例的方法，可以得到高分辨率的成像结果，并且显示的陷落柱位置及形态与实际情况吻合度较高，因此本发明具有较高的精度。