本发明公开了一种基于反射槽波频散相似度的断层成像方法，通过以下步骤实现：步骤一：时频分解；步骤二：求直达槽波频散曲线；步骤三：消除直达槽波；步骤四：成像。本发明充分利用了槽波的频散性质，依据频散性减少了其他噪声对成像的影响，减少了画弧现象，大大提高了断层成像的精度。
1.一种基于反射槽波频散相似度的断层成像方法，其特征在于：通过以下步骤实现：步骤一：时频分解，用广义S变换求出各道地震数据的时频谱其中sn(t)为第n道地震数据的时域形式，Sn(τ,f)为其时频谱，n∈[1,N]，N为总道数；
步骤二：求直达槽波频散曲线，依据各道的炮点位置(xsn,ysn)与检波点位置(xrn,yrn)，求出直达槽波的传播距离Ln再由Ln对Sn(τ,f)进行坐标转换，可求出各道直达槽波的频散曲线，叠加后获得总体数据中直达槽波的频散曲线C(v,f)步骤三：消除直达槽波，根据Ln将C(v,f)从“速度-频率”域转换回“时间-频率”域后，与Sn(τ,f)相减，就可以得到消除了直达槽波后的各道时频谱S’n(τ,f)步骤四：成像，首先选定一个成像区域，区域的一个点P，坐标为(xp,yp)，假设P点是一个反射点，则槽波从第n道数据的炮点(xsn,ysn)出发，在P点反射后，被(xrn,yrn)处的检波器接收到，则传播路径的总长度为如果P点是反射点，则根据Lpn求出的频散曲线Cpn(v,f)应该与C(v,f)一致，因此，将所有道Cpn(v,f)与C(v,f)的相关系数的和Rp作为该成像点的值，当该点确实为反射点时，Rp较大，而当它不是反射点时，Rp较小，对成像区域内所有成像点进行上述处理，得到整个区域的成像结果。
2.根据权利要求1所述的基于反射槽波频散相似度的断层成像方法，其特征在于，所述步骤一具体为：
2)用广义S变换求sn(t)对应的时频谱Sn(τ,f)；
3.根据权利要求1或2所述的基于反射槽波频散相似度的断层成像方法，其特征在于，所述步骤二具体为：
2)从第n道道头中读取炮点位置(xsn,ysn)与检波点位置(xrn,yrn)；
4)时频谱Sn(τ,f)转换为速度-频率谱Cn(v,f)；
6)将所有道的Cn(v,f)叠加得到直达槽波的频散曲线C(v,f)。
4.根据权利要求3所述的基于反射槽波频散相似度的断层成像方法，其特征在于，所述步骤三具体为：
2)令n＝1，读取第n道对应的直达槽波传播距离Ln；
4)由Ln将C(v,f)转换到“时间-频率”域后与Sn(τ,f)相减，得到消除了直达槽波后时频谱S’n(τ,f)；
5.根据权利要求4所述的基于反射槽波频散相似度的断层成像方法，其特征在于，所述步骤四具体为：
2)令x＝1，y＝1，从第1个点开始对成像范围内的每一个点逐个处理；
3)令n＝1，从道头中读取炮点位置(xsn,ysn)与检波点位置(xrn,yrn)；
4)计算槽波从炮点到反射点，再到检波点的传播距离Lpn；
6)求当以(x,y)为坐标的点P为反射点时，第n道反射槽波的频散曲线Cpn(v,f)；
7)求Cpn(v,f)与C(v,f)的相关系数Rpn；
9)将所有道的Rpn叠加后的值Rp作为P点对应的成像值；
10)令x＝x+1或y＝y+1，求下一个点的Rp直到成像区域内所有点成像完毕；
基于反射槽波频散相似度的断层成像方法
技术领域
[0001] 本发明属于煤矿井下物探技术领域，具体涉及一种基于反射槽波频散相似度的断层成像方法。
背景技术
[0002] 理论的反射槽波成像方法为延时求和偏移法。其思想是依据理论频散曲线对各道信号进行相位补偿，然而在计算出延迟时间后将各道信号延时时刻的振幅叠加起来。实际资料处理中，由于频散曲线的理论值与实际值有一定差异，相位补偿的效果并不理想。因此常用方法是依据中心频率进行带通滤波后求包络，再对包络进行延时叠加。当反射槽波较强能形成单独包络时，采用此方法能够获得一定的成像效果。然而，由于井下噪声较强且直达槽波波列较长，反射槽波波列经常被淹没在噪声或直达槽波波列中，无法形成包络，成像效果较差。已有的反射槽波法对断层成像的效果并不理想，难以满足实际需要。
[0003] 事实上，频散性是槽波区别于其他弹性波的重要特性之一。已有方法无论是相位补偿还是求包络，都没有利用频散性，甚至将其视为有害信息而加以消除。针对这一问题，我们研究发现反射对槽波频散的影响并不大，发生反射前后槽波在“速度-频率”域内具有很强的相似性。因此可以利用这一性质，首先提取直达槽波的“速度-频率”谱作为模板，在成像过程中，根据成像点的位置将数据转换到“速度-频率”域内以后计算其与直达槽波“速度-频率”谱的相似性。如果成像点是反射点，则相似性就高。通过计算每一个成像点处反射槽波与直达槽波“速度-频率”谱的相似度，就可以对断层成像。
发明内容
[0004] 有鉴于此，本发明的主要目的在于提供一种基于反射槽波频散相似度的断层成像方法。
[0005] 为达到上述目的，本发明的技术方案是这样实现的：
[0006] 本发明实施例提供一种基于反射槽波频散相似度的断层成像方法，通过以下步骤实现：
[0007] 步骤一：时频分解，用广义S变换求出各道地震数据的时频谱
[0008]
[0009] 其中sn(t)为第n道地震数据的时域形式，Sn(τ,f)为其时频谱，n∈[1,N]，N为总道数；步骤二：求直达槽波频散曲线，依据各道的炮点位置(xsn,ysn)与检波点位置(xrn,yrn)，求出直达槽波的传播距离Ln
[0010]
[0011] 再由Ln对Sn(τ,f)进行坐标转换，可求出各道直达槽波的频散曲线，叠加后获得总体数据中直达槽波的频散曲线C(v,f)
[0012]
[0013] 步骤三：消除直达槽波，根据Ln将C(v,f)从“速度-频率”域转换回“时间-频率”域后，与Sn(τ,f)相减，就可以得到消除了直达槽波后的各道时频谱S’n(τ,f)[0014]
[0015] 步骤四：成像，首先选定一个成像区域，区域的一个点P，坐标为(xp,yp)，假设P点是一个反射点，则槽波从第n道数据的炮点(xsn,ysn)出发，在P点反射后，被(xrn,yrn)处的检波器接收到，则传播路径的总长度为
[0016]
[0017] 如果P点是反射点，则根据Lpn求出的频散曲线Cpn(v,f)应该与C(v,f)一致，因此，将所有道Cpn(v,f)与C(v,f)的相关系数的和Rp作为该成像点的值，
[0018]
[0019] 当该点确实为反射点时，Rp较大，而当它不是反射点时，Rp较小，对成像区域内所有成像点进行上述处理，得到整个区域的成像结果。
[0020] 上述方案中，所述步骤一具体为：
[0021] 1)令n＝1，读取第n道地震数据sn(t)；
[0022] 2)用广义S变换求sn(t)对应的时频谱Sn(τ,f)；
[0023] 3)保存Sn(τ,f)；
[0024] 4)令n＝n+1，处理下一道直到所有道处理完毕。
[0025] 上述方案中，所述步骤二具体为：
[0026] 1)令n＝1，读取第n道时频谱Sn(τ,f)；
[0027] 2)从第n道道头中读取炮点位置(xsn,ysn)与检波点位置(xrn,yrn)；
[0028] 3)计算直达槽波的传播距离Ln；
[0029] 4)时频谱Sn(τ,f)转换为速度-频率谱Cn(v,f)；
[0030] 5)令n＝n+1，处理下一道直到所有道处理完毕；
[0031] 6)将所有道的Cn(v,f)叠加得到直达槽波的频散曲线C(v,f)。
[0032] 上述方案中，所述步骤三具体为：
[0033] 1)读取直达槽波频散曲线C(v,f)；
[0034] 2)令n＝1，读取第n道对应的直达槽波传播距离Ln；
[0035] 3)读取第n道时频谱Sn(τ,f)；由Ln将C(v,f)转换到“时间-频率”域后与Sn(τ,f)相减，得到消除了直达槽波后时频谱S’n(τ,f)；
[0036] 4)保存S’n(τ,f)；
[0037] 5)令n＝n+1，处理下一道直到所有道处理完毕。
[0038] 上述方案中，所述步骤四具体为：
[0039] 1)选定成像范围；
[0040] 2)令x＝1，y＝1，从第1个点开始对成像范围内的每一个点逐个处理；
[0041] 3)令n＝1，从道头中读取炮点位置(xsn,ysn)与检波点位置(xrn,yrn)；
[0042] 4)计算槽波从炮点到反射点，再到检波点的传播距离Lpn；
[0043] 5)读取S’n(τ,f)；
[0044] 6)求当以(x,y)为坐标的点P为反射点时，第n道反射槽波的频散曲线Cpn(v,f)；
[0045] 7)求Cpn(v,f)与C(v,f)的相关系数Rpn；
[0046] 8)令n＝n+1，处理下一道直到所有道处理完毕；
[0047] 9)将所有道的Rpn叠加后的值Rp作为P点对应的成像值；
[0048] 10)令x＝x+1或y＝y+1，求下一个点的Rp直到成像区域内所有点成像完毕；
[0049] 11)最终显示成像结果。
[0050] 与现有技术相比，本发明的有益效果：
[0051] 本发明充分利用了槽波的频散性质，依据频散性减少了其他噪声对成像的影响，减少了画弧现象，大大提高了断层成像的精度。
附图说明
[0052] 图1为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的总流程图；
[0053] 图2为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的时频分解流程图；
[0054] 图3为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的求直达槽波频散曲线流程图；
[0055] 图4为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的直达槽波压制流程图；
[0056] 图5为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的反射槽波成像流程图；
[0057] 图6为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的断层位置与观测系统布置图；
[0058] 图7为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的第20炮模型数据；
[0059] 图8为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的直达槽波频散曲线；
[0060] 图9为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的按直达波路径求出第30道的速度-频率谱；
[0061] 图10为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的压制直达槽波后结果；
[0062] 图11为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的按反射路径求出的速度-频率谱；
[0063] 图12为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的已有方法的断层成像结果；
[0064] 图13为本发明实施例提供一种基于反射槽波频散相似度的断层成像方法的本发明的断层成像结果。
具体实施方式
[0065] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0066] 本发明实施例提供一种基于反射槽波频散相似度的断层成像方法，通过以下步骤实现：
[0067] 步骤一：用广义S变换求出各道地震数据的时频谱
[0068]
[0069] 其中sn(t)为第n道地震数据的时域形式，Sn(τ,f)为其时频谱。n∈[1,N]，N为总道数；
[0070] 具体地，如图2所示，令n＝1，读取第n道地震数据sn(t)；用广义S变换求sn(t)对应的时频谱Sn(τ,f)；保存Sn(τ,f)；令n＝n+1，处理下一道直到所有道处理完毕。
[0071] 步骤二：依据各道的炮点位置(xsn,ysn)与检波点位置(xrn,yrn)，求出直达槽波的传播距离Ln
[0072]
[0073] 具体地，如图3所示，令n＝1，读取第n道时频谱Sn(τ,f)；从第n道道头中读取炮点位置(xsn,ysn)与检波点位置(xrn,yrn)；计算直达槽波的传播距离Ln；时频谱Sn(τ,f)转换为速度-频率谱Cn(v,f)；令n＝n+1，处理下一道直到所有道处理完毕；将所有道的Cn(v,f)叠加得到直达槽波的频散曲线C(v,f)。
[0074] 步骤三：由Ln对Sn(τ,f)进行坐标转换，可求出各道直达槽波的频散曲线，叠加后获得总体数据中直达槽波的频散曲线C(v,f)
[0075]
[0076] 再根据Ln将C(v,f)从“速度-频率”域转换回“时间-频率”域后，与Sn(τ,f)相减，就可以得到消除了直达槽波后的各道时频谱S’n(τ,f)
[0077]
[0078] 具体地，如图4所示，读取直达槽波频散曲线C(v,f)；令n＝1，读取第n道对应的直达槽波传播距离Ln；读取第n道时频谱Sn(τ,f)；由Ln将C(v,f)转换到“时间-频率”域后与Sn(τ,f)相减，得到消除了直达槽波后时频谱S’n(τ,f)；保存S’n(τ,f)；令n＝n+1，处理下一道直到所有道处理完毕。步骤四：成像，首先选定一个成
[0079] 像区域，区域的一个点P，坐标为(xp,yp)。
[0080] 假设P点是一个反射点，则槽波从第n道数据的炮点(xsn,ysn)出发，在P点反射后，被(xrn,yrn)处的检波器接收到，则传播路径的总长度为
[0081]
[0082] 如果P点是反射点，则根据Lpn求出的频散曲线Cpn(v,f)应该与C(v,f)一致，因此，将所有道Cpn(v,f)与C(v,f)的相关系数的和Rp作为该成像点的值，
[0083]
[0084] 当该点确实为反射点时，Rp较大，而当它不是反射点时，Rp较小，对成像区域内所有成像点进行上述处理，得到整个区域的成像结果。
[0085] 具体地，如图5所示，读取直达槽波频散曲线C(v,f)；选定成像范围；令x＝1，y＝1，从第1个点开始对成像范围内的每一个点逐个处理；令n＝1，从道头中读取炮点位置(xsn,ysn)与检波点位置(xrn,yrn)；计算槽波从炮点到反射点，再到检波点的传播距离Lpn；读取S’n(τ,f)；求当以(x,y)为坐标的点P为反射点时，第n道反射槽波的频散曲线Cpn(v,f)；求Cpn(v,f)与C(v,f)的相关系数Rpn；令n＝n+1，处理下一道直到所有道处理完毕；将所有道的Rpn叠加后的值Rp作为P点对应的成像值；令x＝x+1或y＝y+1，求下一个点的Rp直到成像区域内所有点成像完毕；最终显示成像结果。
[0086] 下面以理论模型为例说明本发明效果：模型由顶板、煤层、底板三层[0087] 组成。围岩纵波速度3800m/s，横波速度2000m/s，密度2.4g/cm3；煤层厚度5m，纵波速度1800m/s，横波速度1100m/s，密度1.4g/cm3。模型中包含一条斜向断层，落差5m。断层位置与观测系统布置如图6所示。图7为根据三维弹性波正演模拟方法得到的第20炮正演模拟的结果，从中可以分辨出直达槽波与断层的反射槽波。图8由所有炮叠加得到的直达槽波“速度-频率”谱，其中黑色线条为根据物性参数计算出的理论频散曲线。图9是第20炮的第30道数据时频谱转换到“速度-频率谱”的结果，压缩其中的直达槽波，结果如图10所示。根据反射的几何关系，第20炮第20道数据的反射槽波应来自于断层上的点(x＝429,y＝137)。
如果将这一点作为反射点，求第20炮第30道数据的“速度-频率谱”，结果如图11所示。可以看出反射槽波与直达槽波的频散曲线具有很强的相似性。本发明正是利用这一相似性对断层成像的。图12是用已有方法成像的结果。图13是用本发明方法成像的结果，可以看出与已有方法相比，新方法的成像结果减少了画弧效应，具有更高的分辨率，能够更准确反映断层的发育情况。断层的像在距离模型右边界20m范围内有能量减弱的现象，这是因为能够在这段范围形成有效反射的炮点已经超出了模型范围，所以这段断层无法成像是正常现象。
[0088] 以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。
