基于波形包络的各向同性介质超声波速度自动计算方法与流程

本发明涉及一种基于波形包络的各向同性介质超声波速度自动计算方法。技术背景随着岩石物理研究的兴起，岩石的速度信息是最受关注的参数之一，获得岩石速度的方法有很多，超声波波速测试便是一种快速、便捷、无损的波速测试方法，已经得到广泛应用。超声波测量岩石波速的原理是：在测试样品两端分别放上压电换能器探头，一端发射特定频率的超声脉冲信号，另一端则接收信号，通过分析超声波在样品中传播的时间，便可求得岩石特定频率的超声波速度，这个速度也称相速度。实际测试过程中，由于边界条件等的作用，穿过岩石的波很少是单频的，相反，它们具有有限频宽，这些不同频率的单频波将相互叠加并以波包的形式向前传播，由此产生了群速度，即具有有限频宽的某一调制载波包络面的传播速度。相速度与群速度不是一个概念，但是，在各向同性、均匀、弹性介质中，岩石的群速度和相速度在大小、方向上是相同的。现有的超声波波速测试分析中，主要侧重于通过人工拾取接收波形的起跳时间作为波在岩石中传播的时间，进而计算岩石样品的超声波速度。实验室条件下岩石中传播的超声波频率并不单一且压电换能器探头尺寸有限，使得其超声波初至拾取存在人为误差且效率低。因此，本发明通过自动拾取超声波波包初至时间，获得超声波在岩石中的传播时间，进而达到计算各向同性介质超声波速度的目的。目前还未见到通过拾取波形包络初至计算岩石超声波速度的方法，这里建立一种基于波形包络的各向同性介质超声波速度自动计算方法，通过本法明方法，能简单、快速的获得各向同性介质超声波速度。技术实现要素：本发明主要目的是计算各向同性介质超声波速度，其特征在于通过求取波形包络，通过本发明方法自动拾取波包初至时间，进而计算各向同性介质超声波速度。它给出了超声波波形包络的求取方法、波包初至时间自动拾取方法、给出了各向同性介质超声波速度求取方法，具有简单、精确的特点。本发明解决其技术问题所采取的技术方案如下:(1)数据导入：导入超声波波形数据，获得原始超声波波形数据；(2)数据滤波：构建带通滤波器，获得噪声压制后的波形数据；(3)波形包络计算：基于滤波后的数据，利用本发明方法求取波形的上包络和下包络；(4)波包初至拾取：利用本发明方法自动拾取波包初至时间；(5)超声波速度计算：利用本发明方法计算样品超声波速度。基于波形包络的各向同性介质超声波速度自动计算方法，通过求取波形包络，通过自动拾取超声波波包初至时间，达到计算样品超声波速度的目的，所述样品超声波速度计算方法具体实现步骤如下：步骤一：导入波形数据x1(t(i)),(i＝1,2…N)；步骤二：通过构建带通滤波器，获得去除噪声的超声波波形数据x2(t(i)),(i＝1,2…N)；步骤三：求取x2(t(i)),(i＝1,2…N)的所有极大值以及极小值点；步骤四：利用三次样条插值算法对上述求得得极大值与极小值进行插值处理，计算超声波波形的上包络E_top(t(i)),(i＝1,2…N)与下包络E_bot(t(i)),(i＝1,2…N)；步骤五：取上下包络之差ΔE(t(i))＝E_top(t(i))-E_bot(t(i))其中，i＝1,2…N；步骤六：利用步骤三方法求取ΔE(t(i)),(i＝1,2…N)所有极小值点；步骤七：给定待测样品的高度信息(L)、速度区间(Lower_velocity<Velocity<Higher_velocity)，计算波包初至时间上下界限(Lower_time<t<Higher_time)；步骤八：通过构建超声波波包初至时间应满足的条件，对ΔE(t(i)),(i＝1,2…N)的极小值点做筛选，拾取波包初至时间FA(FirstArrival)；步骤九：计算待测样品的超声波速度Velocity＝L/FA。进一步地，所述步骤二中，通过对原始波形x1(t(i)),(i＝1,2…N)进行快速傅里叶变换，获得x1(t(i)),(i＝1,2…N)的傅里叶变换X1(f(i)),(i＝1,2…N)，构建带通滤波器如下式：令X2(f(i)),(i＝1,2…N)为滤波处理后的频率域信号，则有：X2(f(i))＝X1(f(i))×H(f(i)),(i＝1,2…N)(2)对X2(f(i)),(i＝1,2…N)做傅里叶反变换，获得去除噪声的波形数据x2(t(i)),(i＝1,2…N)。进一步地，所述步骤三中，做x2(t(i)),(i＝1,2…N)的差分diff(t(i)),(i＝1,2…N-1)，即：diff1(t(i))＝x2(t(i+1))-x2(t(i))i＝1,2…N-1(3)令同理，再做diff1(t(i)),(i＝1,2…N-1)的差分diff2(t(i)),(i＝1,2…N-2)，即：diff2(t(i))＝diff1(t(i+1))-diff1(t(i))i＝1,2…N-2(5)令indmax、indmin为x2(t(i)),(i＝1,2…N)的极大值、极小值横坐标(时间)逻辑索引，则有：按照极大值与极小值的定义，数据的第一个数据x2(t(1))和最后一个数据x2(t(N))不可能是极值点，那么x2(t(i)),(i＝1,2…N)最终的极大值、极小值横坐标(时间)逻辑索引值：进而可求出x2(t(i)),(i＝1,2…N)的极大值(t(indmax1(i)),x2(t(indmax1(i))),(i＝1,2…N)与极小值点(t(indmin1(i)),x2(t(indmin1(i))),(i＝1,2…N)。进一步地，所述步骤四中，对上述极值点在时间序列t(i),(i＝1,2…N)上做三次样条插值，根据极大值点获得x2(t(i)),(i＝1,2…N)的上包络E_top(t(i))(i＝1,2…N)，根据极小值点获得x2(t(i)),(i＝1,2…N)的下包络E_bot(t(i))(i＝1,2…N)。进一步地，所述步骤五中，做上下包络线之差ΔE(t(i))(i＝1,2…N),如下式：ΔE(t(i))＝E_top(t(i))-E_bot(t(i))i＝1,2…N(10)进一步地，所述步骤六中，利用步骤三方法可求得ΔE(t(i))(i＝1,2…N)的所有极小值点(t(indmin2(i)),(x2(t(indmin2(i))),(i＝1,2…N)。进一步地，所述步骤七中，输入待测样品的速度上下限，以及样品高度信息，根据下式子可求得超声波波包初至时间的上下限区间范围：(Lower_time<FirstArrival<Higher_time)进一步地，所述步骤八中，当极小值点(t(indmin2(i0)),(x2(t(indmin2(i0)))满足如下条件时：①极小值点横坐标t(indmin2(i0))(时间)落在波包初至时间区间内(即Lower_time<t(indmin2(i0))<Higher_time)；②给定阈值Num(足够大)，使得ΔE(t(indmin2(i+1)))≥ΔE(t(indmin2(i))),(i＝i0,i0+1…i0+Num-1)成立；③第一个同时满足上述两个条件的极小值点t(indmin2(i0))。认为该极小值点横坐标t(indmin2(i0))(时间)为波包初至时间FA。进一步地，所述步骤九中，利用本发明方法自动拾取的超声波波包初至时间，通过公式计算获得样品的超声波速度Velocity。本发明的有益效果是：可以自动拾取超声波波包初至时间，通过本方法能简单、快速精确的计算各向同性介质超声波速度。附图说明下面结合附图和实施例对本发明进一步说明。图1超声波速度计算流程图图2实例应用示意图具体实施方式以下将结合附图对本发明进行详细说明，但本发明并不限于此。采用美国GCTS超声测试系统获得样品超声波形(仪器探头固有延时3.2μs)，获得14块样品的超声波波形，挑选1-1号样品作为重点说明实例。实例具体实施步骤如下：步骤一：导入1-1号样品原始超声波形数据x1(t(i)),(i＝1,2…683)；步骤二：通过对原始波形x1(t(i)),(i＝1,2…683)进行快速傅里叶变换，获得x1(t(i)),(i＝1,2…683)的傅里叶变换X1(f(i)),(i＝1,2…683)，构建带通滤波器如下式：令X2(f(i)),(i＝1,2…683)为滤波处理后的频率域信号，则有：X2(f(i))＝X1(f(i))×H(f(i)),(i＝1,2…683)(14)对X2(f(i)),(i＝1,2…683)做傅里叶反变换，获得去除噪声的波形数据x2(t(i)),(i＝1,2…683)(见图2)。步骤三：做x2(t(i)),(i＝1,2…683)的差分diff(t(i)),(i＝1,2…682)，即：diff1(t(i))＝x2(t(i+1))-x2(t(i))i＝1,2…682(15)令同理，再做diff1(t(i)),(i＝1,2…682)的差分diff2(t(i)),(i＝1,2…681)，即：diff2(t(i))＝diff1(t(i+1))-diff1(t(i))i＝1,2…681(17)令indmax、indmin为x2(t(i)),(i＝1,2…683)的极大值、极小值横坐标(时间)逻辑索引，则有：按照极大值极小值的定义，数据的第一个数据x2(t(1))和最后一个数据x2(t(683))不可能是极值点，那么x2(t(i)),(i＝1,2…683)最终的极大值、极小值横坐标(时间)逻辑索引值：进而可求出x2(t(i)),(i＝1,2…683)的极大值(t(indmax1(i)),x2(t(indmax1(i))),(i＝1,2…683)与极小值点(t(indmin1(i)),x2(t(indmin1(i))),(i＝1,2…683)。步骤四：对上述极值点在时间序列t(i),(i＝1,2…683)上做三次样条插值，根据极大值点获得x2(t(i)),(i＝1,2…683)的上包络E_top(t(i))(i＝1,2…683)，根据极小值点获得x2(t(i)),(i＝1,2…683)的下包络E_bot(t(i))(i＝1,2…683)(见图2)。步骤五：做上下包络线之差ΔE(t(i))(i＝1,2…683),如下式：ΔE(t(i))＝E_top(t(i))-E_bot(t(i))i＝1,2…683(22)步骤六：按照步骤三的方法，求得ΔE(t(i))(i＝1,2…683)的所有极小值点。步骤七：输入待测样品的速度上下限，以及样品高度信息，根据下式子可求得超声波波包初至时间的上下限区间范围：那么波包初至时间区间为：2.56e-05s<FA<7.03e-05s(见图2)。步骤八：筛选满足如下条件的极小值点：①极小值点横坐标(时间)落在波包初至时间区间内(即2.56e-05s<t<7.03e-05s)，可知，只有满足条件的(t(83),(x2(t(83))以及(t(139),(x2(t(139))，即(3.28e-05,-1.06e-04)以及(5.52e-05,0.0176)满足条件；②给定阈值Num＝30，使得ΔE(t(i)+1)≥ΔE(t(i)),(i＝i0,i0+1…i0+29)成立；③第一个同时满足上述两条件的极小值点。最终求得极小值点(3.28e-05,-1.06e-04)，即波包初至时间FA＝3.28e-05s(见图2)。步骤九：根据如下公式，计算出样品超声波速度(见图2)：实施例1：山西晋城寺河矿区煤系地层岩石超声波速度计算，于2016年7月10日进行。对研究区内煤系地层岩石(煤、砂岩、泥岩)进行加工制作。根据实验室仪器测量规格，为尽可能降低岩石各向异性对波速测试的影响，对所采样品在垂直或者平行层理方向进行钻、切、磨等加工，制成直径38mm的标准圆柱体(共计14个样品)，磨光两端面使其互相平行且垂直圆柱体轴线。在室温常压下，采用美国GCTS超声测速系统进行测试，获得超声波波形，根据本发明方法计算研究区煤系地层岩石超声波速度值，如表1。表1研究区煤系地层岩石超声波速度计算表编号高度(mm)阈值(Num)速度区间(m/s)速度(m/s)备注1-167.06301000-30002265.54无烟煤，垂直层理方向1-268.21301000-30001937.78无烟煤，垂直层理方向1-367.04301000-30002394.29无烟煤，垂直层理方向1-464.33301000-30002233.68无烟煤，垂直层理方向1-560.21301000-30002403.06无烟煤，垂直层理方向1-660.00301000-30001948.05无烟煤，垂直层理方向2-167.841001000-60001856.09粗砂岩，垂直层理方向2-266.551001000-60001670.01粗砂岩，垂直层理方向2-364.09301000-60003355.50中砂岩，垂直层理方向2-465.31301000-60003899.10中砂岩，垂直层理方向2-563.64301000-60004821.21细砂岩，垂直层理方向2-663.33301000-60005066.40细砂岩，垂直层理方向2-763.72301000-60004192.10粉砂岩，垂直层理方向2-863.22301000-60003209.14粉砂岩，垂直层理方向天然岩石并不完全满足各向同性、均质性等的理想条件，但是针对于满足横向上(HTI)或垂向上(VTI)各向同性的岩石，其速度计算同样适用于本发明方法。因此，本发明一般要求待测样品通过垂直于或者平行层理方向钻取。综上，基于波形包络的各向同性介质超声波速度自动计算方法是一种有效的超声波速度计算方法，以上所述仅为实现本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。当前第1页1 2 3