基于地震槽波全波形反演的煤矿井下构造成像方法及系统与流程

本发明涉及煤矿井下地质异常体探测的矿井安全技术领域，更具体地说，涉及一种基于地震槽波全波形反演的煤矿井下构造成像方法及系统，用于对煤层中的陷落柱、小断层、夹矸、煤层起伏走势等进行探测，除了能用于大构造和大的地质异常体的成像外，能够精确地反演煤层回采工作面内的小构造形态和物性参数(包括纵波速度、横波速度、密度等)。

背景技术：

许多煤矿的开采条件复杂而又困难。由于浅层地质条件较好的煤层已经不多，目前煤炭开采逐步向更深、地质条件更差的方向发展。煤层中的地质异常体，如陷落柱、断层、夹矸及冲刷带等，容易诱发煤矿冒顶、塌方、突水等灾害事故，严重威胁煤矿安全生产，同时也会严重影响综采机械的正常开采工作。因此，开展煤矿井下地质条件的精细探测至关重要。

槽波勘探是目前煤矿井下工作面异常体探测较常用的方法。槽波勘探有透射波法、反射波法和透射波-反射波联合勘探法(胡国泽等，2013)。透射槽波方法主要基于层析成像原理探测工作面内部异常构造，目前应用较多(王伟等，2012；姬广忠等，2014；任亚平，2015)；反射槽波方法有基于常规的包络叠加(王一,2017)和动态道集叠加(buchananetal.,1981)，还有基于偏移类方法(huetal.，2007；王季，2015；姬广忠，2017)。虽然常规槽波勘探方法在煤矿井下工作面探测较为成功，但对于小尺度异常体，诸如小断层、夹矸层、煤层弯曲等仍然表现不佳或无能为力。究其原因主要有两点：(1)基于槽波的常规层析成像、反射叠加及偏移方法常将煤层简化成二维平面问题，而槽波是在三维空间中形成并传播的，检波器接收的地震记录是三维空间地质异常产生的信息。二维成像方法必然会造成反演精度下降甚至产生反演假象。要想准确获得煤层小异常体的位置及形态大小，就需要进行三维成像；(2)包络叠加与动态道集叠加方法一般基于地层层状结构假设，对于煤层中存在的复杂构造成像精度较低；偏移类方法对于复杂构造区域，其射线追踪容易出现多路径问题，而且对速度模型精度要求较高；层析成像方法只利用了地震数据中的旅行时或振幅信息，其得到的成像结果只包含了模型的宏观速度场，对地质异常细节的刻画无能为力。

地震全波形反演作为一种高精度、高分辨率的地下介质物性参数反演方法，是近年来地震成像领域研究的热点之一。计算机硬件水平的提高推动了全波形反演方法的快速发展。全波形反演基于全波场正演模拟，最小化模拟数据和观测数据之间的残差，通过局域最优化迭代方法寻找最优的地下模型。由于全波形反演利用了地震记录中全部的运动学和动力学信息，因而具有精细刻画复杂地质构造的能力。相比于层析成像方法，全波形反演结果分辨率具有极大的提高。

技术实现要素：

为了解决传统槽波勘探方法对于煤层小构造探测分辨率不足的问题，本发提供乐一种煤矿井下煤层异常体在含有槽波的波形条件下的全波形反演方法，以有效地反演出煤层中陷落柱、小断层、夹矸体、煤层起伏等地质异常体的形态、大小以及纵、横波速度和密度参数值。

本发明解决其技术问题所采用的技术方案是：构造一种基于地震槽波全波形反演的煤矿井下构造成像方法，包含如下步骤：

s1、建立煤矿井下的初始模型，初始模型中的模型参数为纵波速度、横波以及密度；

s2、获取回采工作面的地震槽波的观测数据以及模拟数据；其中，观测数据是由多个检波器接收的数据，每个检波器形成一道数据；

s3、将观测数据和模拟数据相减获得残差记录，然后对每一道数据的所有采样点的残值进行平方再累加求和，从而构建出观测数据与模拟数据的二范数形式的目标函数；

s4、利用震源进行一次正演模拟，得到震源波场，然后利用残差记录作为震源进行一次逆向正演模拟，得到残差波场，将震源波场和残差波场进行互相关，求得关于纵波速度、横波以及密度的目标函数梯度；

s5、采用l-bfgs方法构建各所述目标函数梯度的下降方向；

s6、给定三个测试步长，通过所述下降方向计算出测试步长对应的三个模型，分别对三个模型进行正演获得模拟地震记录，然后利用模拟记录和对应的真实记录分别计算出三个目标函数值，再根据计算出的三个目标函数值，通过抛物线拟合方法估计最优步长；

s7、利用求得的目标函数下降方向和最优步长，对上次的模型进行更新，通过不断迭代，当残差小于时，即得到最终的反演模型；其中第一次迭代前的模型为所述初始模型。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，所述观测数据是进行振幅平衡和恢复，以及滤波处理后的数据。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，步骤s3中：

二范数形式的目标函数e(m)的具体形式如下：

其中，ucal和uobs分别表示模拟和观测得到的质点位移，δu为ucal和uobs之差，m为模型参数，xs和xr分别为震源和检波器的坐标位置，ns和ng分别表示震源和检波器数目，t表示时间变量t的地震波最大传播时间，上标“*”表示共轭转置。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，步骤s4中：

目标函数关于密度的梯度通过拉梅参数梯度转化算出：

目标函数关于速度的梯度利用链式法则通过拉梅参数梯度转化得到：

其中，s表示对所有震源波场求和，σij表示震源波场质点应力，σij表示残差反传波场质点应力，ψi表示残差反传波场质点位移，vi表示震源波场质点速度，其中i＝x、y或z，j＝x、y或z，x、y、z分别表示x轴、y轴、z轴，vp为纵波速度，vs为横波速度，ρ为密度，t表示时间变量t的地震波最大传播时间，δλ、δμ和δρ分别表示目标函数关于拉梅参数和密度的梯度；δvp，δvs和δρvel分别表示目标函数关于速度参数和密度的梯度。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，步骤s6中：

抛物线拟合方法首先搜索三个测试步长α1、α2和α3以及对应的目标函数值e1、e2和e3，然后利用抛物线拟合方法计算出最优步长，测试步长和目标函数值满足以下限制条件

其中，m表示模型参数，mk表示第k次迭代的模型参数，mtest1、mtest2和mtest3分别表示三个用于测试的模型参数；dk表示第k次迭代模型参数修改量。

由于弹性波全波形反演需要同时反演三个模型参数，每个单独的模型参数m的更新步长是不同的，首先反演λ、μ和ρ，因而步长αk包含三个独立的步长每个单独模型参数的测试步长用对应的梯度最大值、模型最大值、以及缩放因子β得到：

由于当前第k次模型参数mk已知，其对应的目标函数值已计算出，可以作为e1，对应α1＝0，寻找α2和α3采用经验方法：首先给定一个初始值β0，然后根据公式βn＝β0/2n)通过不断改变β值计算α2和e2，直到e2<e1得到满足，然后根据公式βn+1＝βn(1+1/2)改变β值计算α3和e3，直到e3>e2得到满足，其中n＝0、1、…、n-1，n表示最大的搜索数目；如果成功地找到α1、α2和α3以及对应e1、e2和e3，通过抛物线拟合求解出二元一次方程系数，以得到最优的步长，如果在最大搜索数n之内没有找到α2，则迭代循环结束，但对于α3则不同，当在最大搜索数n之内条件e3>e2不能得到满足，则最后一次的α3用作最优步长。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，在步骤s7中所述迭代的迭代公式为：

mk+1＝mk-αkpkgk

其中，pk为预条件算子，gk为梯度，αk为步长，下标k表示第k次迭代时的数值。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，步骤s5中：

l-bfgs求取完整的hesse矩阵的逆的一个近似矩阵，从一个初始的hesse矩阵近似h0出发，通过迭代策略求取

式中，i表示单位矩阵，sk＝mk+1-mk，yk＝gk+1–gk，h0为初始近似；

l-bfgs只存储m步向量序列去创建的近似，迭代式为：

其中，(上标t表示转置)，是第k次迭代的初始矩阵，m为预设的正整数。

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，

进一步地，在本发明的基于地震槽波全波形反演的煤矿井下构造成像方法中，h0＝i。

向量序列保存的步数m选择3～20，当迭代步k超过m，则只保存最新的m步向量序列。

根据本发明的另一方面，本发明为解决其技术问题，还提供了一种基于地震槽波全波形反演的煤矿井下构造成像系统，采用上述的基于地震槽波全波形反演的煤矿井下构造成像方法，以进行煤矿井下构造成像。

本发明与现有槽波成像技术相比，本发明的优点有：

1、对原始数据无需较多处理，可直接用于反演；

2、能够精确地反演出煤层工作面小尺度地质异常形态、规模及其物性参数(纵波速度、横波速度、密度等)；

3、能够适应不同观测系统，如透射法、反射法以及透射-反射联合法；

4、能够解决常规勘探方法(如直流电法、无线电坑透、地质雷达、瑞雷波勘探等)无法解决的煤层工作面问题，如小落差断层、夹矸层、小陷落柱、煤层起伏等。

5、相比常规基于槽波的叠加法和层析成像方法，本发明算法成像的分辨率能提升一至二个量级。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明的基于地震槽波全波形反演的煤矿井下构造成像方法一优选实施例的流程图；

图2为煤层异常体模型；

图3为槽波全波形反演需要的煤矿工作面地震观测系统的各种布设方式；

图4为透射观测系统；

图5为反演所用初始模型；

图6为陷落柱模型反演结果；

图7为断层块模型反演结果；

图8为夹矸模型反演结果；

图9为煤层起伏模型反演结果。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

参考图1，其为本发明的基于地震槽波全波形反演的煤矿井下构造成像方法一优选实施例的流程图。本发明提出的全波形反演方法是基于数据拟合的一种局域优化迭代方法，本实施例的基于地震槽波全波形反演的煤矿井下构造成像方法主要包括如下步骤：

s1、首先建立煤矿井下的初始模型，模型参数为纵波速度、横波速度和密度。对应地，初始模型具有纵波速度子模型、横波速度子模型和密度子模型。

s2、获取回采工作面的地震槽波的观测数据以及模拟数据；其中，观测数据是由多个检波器接收的数据，每个检波器形成一道数据，而每一道数据又包含很多时间点上的数据。本实施例中的观测数据优选为经过的振幅平衡和恢复处理，以及滤波处理等的预处理后的地震数据数据。

s3、将观测数据和模拟数据相减获得残差记录，然后对每一道数据的所有采样点的残值进行平方再累加求和，从而构建出观测数据与模拟数据的二范数形式的目标函数；

全波形反演的目标函数用基于位移的l2范数表示，具体形式如下：

其中，ucal和uobs分别表示模拟和观测得到的质点位移，δu为ucal和uobs之差，m为模型参数，xs和xr分别为震源和检波器位置，ns和ng分别表示震源和检波器数目，t表示时间变量t的地震波最大传播时间，上标“*”表示共轭转置。

s4、利用震源进行一次正演模拟，得到震源波场，然后利用残差记录作为震源进行一次逆向正演模拟，得到残差波场，将震源波场和残差波场进行互相关，求得关于纵波速度、横波以及密度的目标函数梯度；

s5、采用l-bfgs方法构建各所述目标函数梯度的下降方向；

s6、给定三个测试步长，通过所述下降方向计算出测试步长对应的三个模型，分别对三个模型进行正演获得模拟地震记录，然后利用模拟记录和对应的真实记录分别计算出三个目标函数值，再根据计算出的三个目标函数值，通过抛物线拟合方法估计最优步长；

s7、利用求得的目标函数下降方向和最优步长，对上次的模型进行更新，通过不断迭代，当残差小于时，即得到最终的反演模型；其中第一次迭代前的模型为所述初始模型。第一次迭代是利用初始模型进行正演，以后的迭代是利用上一次迭代更新的模型进行正演，换一句话说就是，上一次迭代(计算)得到的模型，作为下一次迭代的初始模型用于反演计算，这样不停地迭代，直到目标函数收敛到极小值。对模型的更新包括对各个子模型的更新。

煤层槽波条件下全波形反演基本原理为：

全波形反演是非线性反演问题，常采用局域最优化方法求解。局域最优化方法需要从一个初始模型出发，通过迭代反演逐步逼近全局极小值以获取最优的地下介质模型。迭代公式具有如下形式：

mk+1＝mk-αkpkgk(2)

其中，pk为预条件算子，gk为梯度，αk为步长，k表示第k次迭代时的数据。当pk取不同形式时对应不同的最优化算法(如最速下降法、共轭梯度法、l-bfgs法等)，本发明文中采用的是lbfgs方法，因而

目标函数关于拉梅参数和密度的梯度可以通过下式算出：

式中，s表示对所有震源波场求和，σij表示震源波场质点应力，σij表示残差反传波场质点应力，ψi表示残差反传波场质点位移，vi表示震源波场质点速度，其中i＝x、y或z，j＝x、y或z，x、y、z分别表示x轴、y轴、z轴。

在实际工作中，纵波速度vp，横波速度vs和密度ρ更常用一些，目标函数关于速度和密度参数的梯度可以利用链式法则通过拉梅参数梯度转化得到，具体如下：

式(3)、(4)中，δλ、δμ和δρ分别表示目标函数关于拉梅参数和密度的梯度；δvp，δvs和δρvel分别表示目标函数关于速度参数和密度的梯度。

l-bfgs最优化方法。当公式(2)中预条件算子p取hesse矩阵的逆h-1时，就是bfgs法，它是一种拟牛顿方法。bfgs方法并不计算完整的hesse矩阵的逆，而是求取它的一个近似矩阵。从一个初始的hesse矩阵近似h0出发，通过迭代策略求取

式中，i表示单位矩阵，sk＝mk+1-mk，yk＝gk+1–gk，h0为初始近似，一般选取h0＝i。梯度是目标函数对模型参数(vp,vs和ρ)的导数，公式(3)和公式(4)就是梯度g，是梯度的具体表达，梯度在公式(5)中需要应用，用于计算yk＝gk+1–gk，目的是求取这即是目标函数下降方向，相当于迭代公式mk+1＝mk-αkpkgk中的pkgk。

由于bfgs方法需要保存对于大模型参数问题要求十分大的计算机内存空间，而且需要较大的计算消耗。有限内存bfgs(l-bfgs)方法能够有效减少内存占用和计算消耗。l-bfgs方法是一个bfgs的变种，它只存储m(正整数)步向量序列去创建的近似，其迭代式为

其中，(上标t表示转置)，是第k次迭代的初始矩阵，一般选取加权的单位矩阵

向量序列保存的步数m一般选择3～20，当迭代步k超过m，则只保存最新的m步向量序列。

步长估计对于目标函数的收敛十分重要。步长太小往往导致收敛缓慢，步长过大则容易导致收敛过程不稳定，易陷入局部极值。步长搜索有精确方法和非精确方法。精确方法具有较快的收敛速度，但由于计算量十分巨大，并不适合全波形反演最优化问题。抛物线拟合方法不要求过多的正演计算，能够较稳定地收敛。

抛物线拟合方法首先需要搜索三个测试步长α1、α2和α3以及对应的目标函数值e1、e2和e3，然后利用抛物线拟合方法计算出最优步长。测试步长和目标函数需要满足以下限制条件

e2(mtest2＝mk+α2dk)<e1(mtest1＝mk+α1dk),

e3(mtest3＝mk+α3dk)>e2(mtest2＝mk+α2dk).(8)

其中，m表示模型参数，mk表示第k次迭代的模型参数，mtest1、mtest2和mtest3分别表示三个用于测试的模型参数；dk表示第k次迭代模型参数修改量。

由于弹性波全波形反演需要同时反演多个物性参数(即上述的模型参数)，每个单独的物性参数更新步长是不同的。一般做法是，首先反演λ、μ和ρ，因而步长αk包含三个独立的步长每个单独模型参数的测试步长可以用对应的梯度最大值、模型最大值以及缩放因子β得到：

由于当前第k次模型参数mk已知，其对应的目标函数值已计算出，可以作为e1，对应α1＝0。寻找α2和α3采用经验方法。首先给定一个初始值β0，然后根据公式βn＝β0/2n(n＝0、1、…、n-1)通过不断改变β值计算α2和e2，直到e2<e1得到满足。然后根据公式βn+1＝βn(1+1/2)(n＝0、1、…、n-1)改变β值计算α3和e3，直到e3>e2得到满足，n表示最大的搜索数目。如果成功地找到α1、α2和α3以及对应e1、e2和e3，通过抛物线拟合求解出二元一次方程系数，可以得到最优的步长。如果在最大搜索数n之内没有找到α2，则迭代循环结束。但对于α3则不同，当在最大搜索数n之内条件e3>e2不能得到满足，则最后一次的α3用作最优步长。本发明采用经验值β0＝0.01进行搜索，这个值在步长搜索中能够较好地寻找到最优步长，而且目标函数能够较平滑地下降。整个步长搜索程序将会增加正演模拟的次数。为了减少计算开销，步长计算中可以只用一些代表性的炮数去计算目标函数。

实施例子为四个含煤异常体模型，分别为陷落柱模型、断层块模型、夹矸模型和煤层弯曲模型。这些模型大小均为x×y×z＝100m×50m×40m，网格间距在三个方向均为0.5m。所有模型的煤厚均为8m，煤层xy中心面位于z＝20m。模型1中，陷落柱位于煤层xy中心面x＝50m，y＝25m位置，直径10m。模型2中，煤层中出现一个下陷的断层块，断距3m(z方向下陷3m)。断层块在xy面上为矩形，x方向上宽度为40m，y方向上宽度为10m。模型3中，煤层中出现两个平行于煤层面的层状夹矸，二者相互错开，夹矸的xy中心面分别位于z＝18m和z＝21m，夹矸厚度2m，x方向宽度40m，y方向宽度15m。模型4中，煤层面出现不同程度的起伏，z方向起伏程度最大达6m。模型参数设置如表1所示。

表1

本实例实施可选择利用不同观测系统，图3为四种形式的观测系统，分别为透射观测系统(图3(a))、反射观测系统(图3(b))、u型观测系统(图3(c))和四周观测系统(图3(d))。图中，“*”代表炮点，黑色线代表测线。本实施例采用其中的透射观测系统进行测试(图4)。炮点和检波器均布设在煤层中，一共布设3条炮线和3条测线。其中炮点所在平面为y＝6m，沿x方向布设，深度分别为z＝16m、z＝20m和z＝24m，一共15炮。检波点所在平面为y＝44m，同样沿x方向布设，深度分别为z＝16m、z＝20m和z＝24m，一共483个检波器。

具体实施方式为，首先利用真实模型(图2)进行正演获得真实地震数据，然后利用初始模型(图5)进行反演。反演过程包含目标函数建立，目标函数梯度求取，l-bfgs法计算下降方向，步长求取和模型迭代更新。最终得到陷落柱、断层块、夹矸层和煤层起伏模型的反演结果，分别如图6至图9所示。其中图2包括陷落柱模型(model1:collapsecolumn)、断层块模型(model2:faultblock)、夹矸模型(model3:dirtbands)和起伏煤层模型(model4:undulatedcoalseam)四部分的图。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。