基于槽波地震和克里金插值的煤层厚度预测方法及装置与流程

本发明涉及工作面煤层厚度预测方法及装置，属于煤田地质领域，具体涉及一种基于槽波地震和克里金插值的工作面煤层厚度预测方法及装置。
背景技术：
：在现代化长壁综合机械化采煤过程中，采煤工作面煤层厚度及其变化是一项极为重要的地质参数，其重要作用可以概括为以下三个方面：(1)采煤工作面煤厚分布直接影响了开拓、采准巷道的布置及采煤方法的选择，煤厚探查不清可能导致大量的资金浪费，并严重影响开采计划；(2)煤厚变化较大的区域，一般更倾向于发生压力集中及瓦斯聚集现象，如果事先没有做好煤厚勘探及安全隐患处理工作，则容易发生冲击地压和煤与瓦斯突出等安全事故；(3)如果能够超前准确预测工作面的煤厚变化，则可以通过采煤机摇臂的自动控制实现智能化、无人化开采，从而避免截齿切割煤层顶板、保护采矿设备、提高开采效率、降低矸石率、规避瓦斯事故等。当前，基于采准巷道中的煤厚采样点的煤厚插值方法并不能保证采煤工作面中间位置的煤厚预测精度；槽波地震勘探技术反演得到的工作面煤厚分布结果在煤厚采样点处一般不等于实测煤厚值，只能给出煤厚趋势性变化的规律。因此，如何构建一种利用多源数据定量预测采煤工作面煤层厚度方法是综合机械化采煤地质保障工作亟待解决的问题之一。技术实现要素：本发明所要解决的技术问题在于针对上述现有技术中的不足，提供了一种基于槽波地震和克里金插值的煤层厚度预测方法及装置，该方法及装置步骤简单、所需原始资料容易获取且使用方面，所预测的采煤工作面煤厚分布较为精确。为解决上述技术问题，本发明提供以下技术方案：一种基于槽波地震和克里金插值的煤层厚度预测方法，具体包括以下步骤：初级厚度分布图获取步骤，根据煤厚采样数据拟合煤厚与群速度的定量关系，基于所述定量关系将工作面群速度分布图转换为初级厚度分布图；煤厚预测误差图获取步骤，计算所述煤厚采样数据对应采样点处的煤厚预测误差，建立煤厚预测的残差变量；基于所述残差变量，对整个工作面残差变量的分布进行克里金插值，得到煤厚预测误差分布图；工作面煤厚分布图获取步骤，叠加所述初级厚度分布图与所述煤厚预测误差分布图得到工作面煤厚分布图。优选的，所述初级厚度分布图获取步骤具体包括：根据待采煤工作面的平面图建立xoy平面直角坐标系；基于所述平面直角坐标系对待采煤工作面煤层群速度进行成像得到工作面群速度分布图，收集待开采工作面顺槽及切眼中实测的煤厚采样数据，根据煤厚采样数据拟合煤厚与群速度的定量关系，基于所述定量关系将工作面群速度分布图转换为初级厚度分布图。优选的，所述煤厚预测误差图获取步骤中，基于以下子步骤建立煤厚预测的残差变量：以预定步长的网格密度对采煤工作面的初级煤厚分布图进行网格划分，每一个网格点上对应一个初级煤厚预测值；遍历网格，将与采样点平面距离最近的网格点作为该采样点处的煤厚预测值，将所述煤厚预测值与该采样点的真实煤厚值之差作为该采样点处的煤厚预测误差。优选的，所述煤厚预测误差图获取步骤中，对整个工作面残差变量的分布进行克里金插值，得到煤厚预测误差分布图包括：按照预定步长的网格密度对采煤工作面的平面图进行网格划分；基于下式计算两两样本点之间的距离dij与半方差rij；其中，n表示可用于煤厚预测的总煤厚采样点数，e表示求期望；xi、yi和di分别为第i个采样点(i＝1,2,…,n)的x坐标、y坐标和实测煤厚值；以点间距离dij为x轴、以变差rij为y轴，把每一个点对投影到所建立的坐标系中，可以得到变差函数点云图；拟合距离与半方差之间的关系，把x轴分成固定间隔的区间，对落入每个区间段点的变差rij求均值，将求取的均值作为相应区段中点的变差值，然后将这些点连接起来，得到实验变差曲线。对实验变差曲线进行拟合，将拟合结果记为r＝r(d)；对各待预测点执行以下步骤以得到整个工作面的煤厚预测残差网格数据：对于某一待预测网格点，根据拟合的实验变差函数r＝r(d)，计算它与第i个采样点(i＝1,2,…,n)间的半方差ri0；基于下式，求解该网格点处的最优系数λi；基于下式计算待预测网格点处的煤厚预测残差em,n；其中m和n分别为待预测网格点在坐标系中的横、纵轴坐标，ei表示第i个采样点处的煤厚预测误差。优选的，所述工作面煤厚分布图获取步骤包括：按照预定步长的网格密度分别对采煤工作面的初级煤厚分布图和煤厚预测误差分布图进行网格划分，每一个网格点上对应一个初级煤厚预测值和一个煤厚预测误差；对同一网格点上的初级煤厚预测值和煤厚预测误差进行相加，得到该网格点上的最终煤厚预测值。一种基于槽波地震和克里金插值的煤层厚度预测装置，具体包括以下模块：初级厚度分布图获取模块，根据煤厚采样数据拟合煤厚与群速度的定量关系，基于所述定量关系将工作面群速度分布图转换为初级厚度分布图；煤厚预测误差图获取模块，计算所述煤厚采样数据对应采样点处的煤厚预测误差，建立煤厚预测的残差变量；基于所述残差变量，对整个工作面残差变量的分布进行克里金插值，得到煤厚预测误差分布图；工作面煤厚分布图获取模块，叠加所述初级厚度分布图与所述煤厚预测误差分布图得到工作面煤厚分布图。优选的，所述初级厚度分布图获取模块执行以下子步骤：根据待采煤工作面的平面图建立xoy平面直角坐标系；基于所述平面直角坐标系对待采煤工作面煤层群速度进行成像得到工作面群速度分布图，采集待开采工作面煤厚采样数据，根据煤厚采样数据拟合煤厚与群速度的定量关系，基于所述定量关系将工作面群速度分布图转换为初级厚度分布图。优选的，所述煤厚预测误差图获取模块中，基于以下步骤建立煤厚预测的残差变量：以预定步长的网格密度对采煤工作面的初级煤厚分布图进行网格划分，每一个网格点上对应一个初级煤厚预测值；遍历网格，将与采样点平面距离最近的网格点作为该采样点处的煤厚预测值，将所述煤厚预测值与该采样点的真实煤厚值之差作为该采样点处的煤厚预测误差。优选的，所述煤厚预测误差图获取模块中，对整个工作面残差变量的分布进行克里金插值，得到煤厚预测误差分布图包括：按照预定步长的网格密度对采煤工作面的平面图进行网格划分；基于下式计算两两样本点之间的距离dij与半方差rij；其中，n表示可用于煤厚预测的总煤厚采样点数，e表示求期望；xi、yi和di分别为第i个采样点(i＝1,2,…,n)的x坐标、y坐标和实测煤厚值；以点间距离dij为x轴、以变差rij为y轴，把每一个点对投影到所建立的坐标系中，可以得到变差函数点云图；拟合距离与半方差之间的关系，把x轴分成固定间隔的区间，对落入每个区间段点的变差rij求均值，将求取的均值作为相应区段中点的变差值，然后将这些点连接起来，得到实验变差曲线。对实验变差曲线进行拟合，将拟合结果记为r＝r(d)；对各待预测点执行以下步骤以得到整个工作面的煤厚预测残差网格数据：对于某一待预测网格点，根据拟合的实验变差函数r＝r(d)，计算它与第i个采样点(i＝1,2,…,n)间的半方差ri0；基于下式，求解该网格点处的最优系数λi；基于下式计算待预测网格点处的煤厚预测残差em,n；其中m和n分别为待预测网格点在坐标系中的横、纵轴坐标，ei表示第i个采样点处的煤厚预测误差。优选的，所述工作面煤厚分布图获取模块执行以下步骤：按照预定步长的网格密度分别对采煤工作面的初级煤厚分布图和煤厚预测误差分布图进行网格划分，每一个网格点上对应一个初级煤厚预测值和一个煤厚预测误差；对同一网格点上的初级煤厚预测值和煤厚预测误差进行相加，得到该网格点上的最终煤厚预测值。因此，本发明具备以下有益效果：本发明提供的利用槽波地震和克里金插值定量预测工作面煤层厚度的方法，相对于现有技术，具有如下优势：1、所提出的煤厚预测方法，综合利用了采煤工作面的煤厚采样数据和槽波地震技术反演得到的煤厚趋势性分布数据，相较于其他方法，所利用的含煤厚信息数据更多，综合预测结果更加接近采煤工作面真实煤厚分布情况；2、该方法继承了克里金插值方法无偏、最优等优点及槽波地震技术能较为准确反演整个工作面煤厚变化趋势的优点，在靠近煤厚采样点的区域，煤厚预测结果主要取决于煤厚采样点，而在远离煤厚采样点的采煤工作面中间区域，煤厚预测结果主要取决于槽波地震技术的煤厚反演结果；3、所需资料少且常见，实现简单，该方法只需提供采煤工作面平面图、采准巷道写实图和槽波地震探勘资料就能较为精确地预测整个工作面煤厚分布，绝大多数矿山都具备上述三种资料，无需其它配套工作。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1是本发明方法的实施流程图；图2是实施例采煤工作面煤厚采样点分布图，横轴是采煤工作面长度方向，纵轴是采煤工作面宽度方向，图中圆圈是采准巷道中煤厚采样点，实心圆点表示落在工作面范围内的穿煤钻孔，“▲”表示工作面回采过程中实测的煤厚点，用于后续预测误差计算及对比；图3是实施例采煤工作面槽波地震勘探观测系统；图4是实施例不同频率下群速度与煤厚理论数学关系，横轴表示群速度(单位：m/s)，纵轴表示煤层厚度(单位：m)，(a)图对应频率65、75、85、95hz，(b)图对应频率100、200、300、400hz；图5是实施例采煤工作面群速度分布图，横轴是采煤工作面长度方向，纵轴是采煤工作面宽度方向；图6是实施例煤厚与群速度拟合曲线，横轴表示群速度(单位：m/s)，纵轴表示煤层厚度(单位：m)，图中“×”表示煤厚采样点；图7是实施例槽波地震勘探得到的初级煤厚分布图，横轴是采煤工作面长度方向，纵轴是采煤工作面宽度方向；图8是实施例煤厚预测残差变量克里金插值结果，横轴是采煤工作面长度方向，纵轴是采煤工作面宽度方向；图9是实施例最终煤厚预测结果，横轴是采煤工作面长度方向，纵轴是采煤工作面宽度方向；图10是实施例槽波地震勘探(a)和本方法(b)的煤厚预测误差直方图，横轴表示煤厚预测误差(单位：m)，纵轴表示频数。具体实施方式下面将结合本发明实施例及其附图，对本发明的技术方案、工作流程进行详细地描述。本发明的保护范围不受以下实例的限制。本发明基于采煤工作面平面图、采准巷道写实图和槽波地震探勘资料，首先对采煤工作面新建xoy直角坐标系，然后对槽波地震勘探资料进行处理及成像，得到初级煤厚分布图，然后基于采煤工作面现有的煤厚采样点，构建煤厚预测的残差变量，然后采用克里金插值方法对整个工作面范围内的残差变量进行插值，得到煤厚预测误差分布图，最后将初级厚度分布图与预测误差分布图叠加，得到最终的工作面煤厚分布图。如图1所示，本发明采用如下技术方案：步骤1：根据待采煤工作面的平面图建立xoy平面直角坐标系，未指明情况下后续所述坐标均指该新建的xoy平面直角坐标系下的坐标。所述待采煤工作面的平面图一般近似为标准长方形，长度方向与工作面运输巷和工作面回风巷平行，宽度方向与工作面切眼平行；所述新建xoy平面直角坐标系，以待采煤工作面左下角为坐标原点，以工作面长度方向为x轴，以工作面宽度方向为y轴。步骤2：根据步骤1中的新建xoy平面直角坐标系及槽波地震资料的观测系统和旅行时资料，选择合适频率采用联合迭代重建层析成像(sirt)技术对待采煤工作面煤层群速度进行成像。选择合适频率：(1)基于公式(1)、(2)和(3)，绘制不同频率f下的“煤厚d-群速度u”理论曲线。其中，ω为角频率，数值上等于2πf；γi(i＝1,2,3)分别为顶板岩层、煤层、底板岩层的传播因子，数值上等于其中βi(i＝1,2,3)分别为顶板岩层、煤层、底板岩层的剪切波速度，c为相速度；d为待采煤工作面煤厚；μi(i＝1,2,3)分别为顶板岩层、煤层、底板岩层的刚度系数，数值上等于ρiβi2(i＝1,2,3)，其中ρi(i＝1,2,3)分别为顶板岩层、煤层、底板岩层的密度；n为阶数；u为群速度；i为虚数单位。为了简化，一般认为底板岩层和顶板岩层具有相同的物理力学参数(即μ1＝μ3，γ1＝γ3)，且n取等于0。(2)选择合适频率的原则。综合考虑不同频率f下的“煤厚d-群速度u”理论曲线簇及采煤工作面采准巷道掘进过程中实测的煤厚分布区间，找出对所预测煤层厚度区间分辨率最高的频率，并将之用于后续的群速度成像。也就是，合适频率f下“煤厚d-群速度u”理论曲线的负相关段所对应的煤厚区间应该与实测的煤厚分布区间一致性最强。步骤3：收集待开采工作面煤厚采样数据，具体包括工作面运输巷、工作面回风巷、工作面联络巷和切眼等巷道在掘进过程中的煤厚测量数据及位于工作面范围内的地表穿煤勘探钻孔数据，根据煤厚采样数据拟合煤厚与群速度的定量关系。所述煤厚采样数据，按照格式(xi,yi,di)进行整理，其中xi、yi和di分别为第i个采样点(i＝1,2,…,n)的x坐标、y坐标和实测煤厚值，n为煤厚采样点总数。根据煤厚采样数据拟合煤厚与群速度的定量关系：根据采煤工作面煤层群速度成像结果，建立各煤厚采样点处群速度于煤厚关系表，为后续拟合两者关系做准备。表1群速度与煤厚关系表煤厚采样点编号采样点平面坐标煤厚值(m)群速度(m/s)1(x1,y1)d1u12(x2,y2)d2u2…………根据表1中的煤厚值列及群速度列，选择函数类别拟合“煤厚d-群速度u”函数关系，推荐选用幂函数模型，即d＝a·ub，其中a和b为待拟合系数。步骤4：根据步骤3中的拟合函数，将工作面群速度分布图转换为初级厚度分布图。具体转换方法为：步骤401：按照1m×1m的网格密度对采煤工作面的群速度分布图进行网格划分，每一个网格点上对应一个群速度值；步骤402：对于某一网格点，按照拟合后的公式d＝a·ub及该网格点所对应的群速度值u，计算该点的煤厚初级预测值；步骤403：对所有网格点重复402步骤，最终得到初级煤厚网格点数据，稍加处理即可得到煤厚初级煤厚等值线图，或初级煤厚云图。步骤5：根据步骤4中的初级煤厚分布图，提取各采样点处的煤厚预测值，并根据实际煤厚值计算对应采样点处的煤厚预测误差，建立煤厚预测的残差变量。具体建立残差变量方法为：步骤501：按照1m×1m的网格密度对采煤工作面的初级煤厚分布图进行网格划分，每一个网格点上对应一个初级煤厚预测值；步骤502：对于某一煤厚采样点，遍历所有网格点，找到与采样点平面距离(xoy平面)最近的网格点作为该采样点处的初级煤厚预测值，然后基于公式(4)用该预测值减去该采样点真实煤厚值，得到该点处的煤厚预测误差。ei＝d1i-dii＝1,2,...,n，式(4)其中ei、d1i和di分别为第i个采样点处煤厚预测误差、初级煤厚预测值和实测煤厚值。步骤503：对所有煤厚采样点重复步骤502，最终得到煤厚预测的残差变量。步骤6：基于步骤5中的煤厚预测残差变量，利用克里金插值法对整个工作面残差变量的分布进行插值，得到煤厚预测误差分布图。具体插值步骤为：步骤601：按照1m×1m的网格密度对采煤工作面的平面图进行网格划分；步骤602：对于观测数据，根据公式(5)、(6)计算两两样本点之间的距离dij与半方差rij；其中，n表示可用于煤厚预测的总煤厚采样点数，e表示求期望。计算完毕后即可得到n2个(dij,rij)数据对。以点间距离dij为x轴、以变差rij为y轴，把每一个点对投影到所建立的坐标系中，可以得到变差函数点云图。步骤603：选择数学模型，拟合距离与半方差之间的关系；把x轴分成固定间隔的区间，对落入每个区间段点的变差rij求均值，将求取的均值作为相应区段中点的变差值，然后将这些点连接起来，得到实验变差曲线。根据上述点集分布规律，选择合适的数学模型，对实验变差曲线进行拟合，将拟合结果记为r＝r(d)。步骤604：对于某一待预测网格点，根据拟合的实验变差函数r＝r(d)，计算它与第i个采样点(i＝1,2,…,n)间的半方差ri0；步骤605：根据公式(7)，求解该网格点处的最优系数λi；步骤606：按照公式(8)计算待预测网格点处的煤厚预测残差em,n；其中m和n分别为待预测网格点在坐标系中的横、纵轴坐标。步骤607：对其他所有网格点重复步骤604～606，得到整个工作面的煤厚预测残差网格数据，稍加处理即可得到预测残差等值线图或预测残差分布图。步骤7：将步骤4中的初级厚度分布图与步骤5中的预测残差分布图叠加，得到最终的工作面煤厚分布图。具体叠加方法为：步骤701：按照1m×1m的网格密度对采煤工作面的初级煤厚分布图进行网格划分，每一个网格点上对应一个初级煤厚预测值；步骤702：按照1m×1m的网格密度对采煤工作面的煤厚预测误差分布图进行网格划分，每一个网格点上对应一个煤厚预测误差；步骤703：将步骤701和步骤702中的同一网格点(横、纵坐标相等分别相等的两个网格点)上的初级煤厚预测值和煤厚预测误差进行相加，得到该网格点上的最终煤厚预测值；步骤704：对所有网格点按照步骤703进行运算，每个网格点上均对应一个最终煤厚预测值，对该网格数据绘制等值线图得到最终的煤厚分布图。实施例：一、工作面概况玉华煤矿位于陕西省铜川市，其2408采煤工作面位于2号盘区。2408工作面长度1775m，宽度240m，采准巷道掘进过程中实测的煤厚采样点显示煤厚分布区间为6.0～13.5m，平均煤厚7.20m.如图2所示：2408采煤工作面运输巷、回风巷和切眼在掘进过程中分别实测了42、38和5个煤厚采样点，如图2中圆圈所示；另外，还有两个地表穿煤勘探钻孔位于2408采煤工作面范围内，如图2中实心点所示；回采过程中，还实测了94煤厚采样点，如图2中“×”所示。其中前两种可以用于预测煤厚，而这三种数据均可用于评价煤厚预测精度。二、槽波地震观测系统基于正演模型，采用中煤科工集团西安研究院有限公司自主研发的geocoal槽波数据处理软件设计了本次槽波地震勘探的观测系统，设计简图如图3所示。沿着运输巷和回风巷每隔10m布置一台检波器，共计354个，如图3中的实心三角形所示；同样地，沿着运输巷和回风巷每隔30m布置一个炮点，共计118个，如图3中的实心点所示。系统采样间隔为0.25ms，采集时长为2s。表2槽波地震观测系统详细情况地震波射线覆盖了整个采煤工作面。槽波勘探实施过程中，总共记录了115个有效的炮集，记录的地震道信噪比较高。三、槽波地震煤厚反演实测2408采煤工作面顶、底板及煤层物理力学参数如表3所示。表3顶、底板及煤层物理力学参数根据上述参数绘制不同频率下的“煤厚d-群速度u”理论曲线簇，如图4所示。由图4可知，75hz下“煤厚d-群速度u”理论曲线负相关段所对应的煤厚区间与实测煤厚采样点所反映的煤厚区间一致性最好，因此选择75hz对采煤工作面群速度进行成像，结果如图5所示。将采准巷道中实测的煤厚采样点及两个地表穿煤勘探钻孔(即图2中的圆圈与实心点)所对应的煤厚和群速度值按照幂函数关系进行拟合，拟合结果为d＝1305u-0.7313。将群速度分布图按照1m×1m的网格密度进行离散，每个网格点上对应一个群速度值，按照拟合公式d＝1305u-0.7313，计算该网格点对应的初级煤厚预测值。对所有网格点重复该计算过程，得到初级煤厚预测的网格点数据，稍加处理即得到采煤工作面初级煤厚分布图，如图7所示。四、初级煤厚预测残差变量插值根据槽波地震勘探反演得到的初级煤厚分布图，以采准巷道实测煤厚采样点及两个地表穿煤钻孔揭露的煤厚为参考，计算各采样点处的煤厚预测残值，并整理得到初级煤厚预测残差变量表。表4初级煤厚预测残差变量五、残差变量插值及最终煤厚预测基于表4中的各煤厚采样点处的平面坐标信息及煤厚预测残差值，采用克里金方法，对整个采煤工作面范围内的煤厚预测残差进行插值，插值结果如图8所示。将图7与图8进行叠合，得到最终预测的煤厚分布图，结果如图9所示。六、误差评价与验证对于一般槽波地震勘探技术而言，煤厚预测工作在得到本发明中的初级煤厚分布图就结束了，即本发明中的步骤4。因此为了直观反映本发明的优越性，基于2408采煤工作面所有实测煤厚采样点(包括回采过程中实测的94煤厚采样点)，计算各采样点处的煤厚预测误差，并据此评价常规槽波地震煤厚反演方法与本发明的煤厚预测精度。统计结果如图10所示，其中图10(a)为槽波地震勘探的煤厚反演结果(初级煤厚预测分布图)的煤厚预测误差直方图，图10(b)为本发明方法的煤厚预测误差直方图。对比可见，本发明方法的煤厚预测精度明显高于常规槽波地震勘探的煤厚反演结果。当前第1页12