一种反射声波测井的波场分离方法及装置与流程

本发明涉及反射声波测井技术领域，尤其涉及一种反射声波测井的波场分离方法及装置。

背景技术：

随着油气勘探开发的不断深入，人们对地球物理资料的分辨率与精度要求越来越高。常规测井技术虽然分辨率与精度都很高，但受探测范围与探测尺度的限制，无法探测井壁附近以外的地区。

在井中通过声波测井对井周围数米到数十米范围内的地层构造及地质体进行反射波成像的技术称为反射声波测井技术，它把常规测井技术的测量范围从距井壁一米左右提高到几十米。该技术应用前景十分巨大，可以用来显示与井相交的地质界面；探测井旁的倾斜地层界面、裂缝、断层等；在水平井还可以用来追踪油储边界等。

利用常规的地震处理方法对反射波数据进行偏移处理，可以得到地质构造的图像。然而，不同于地震测量方法，反射声波测井的声源和接收器均布置在同一井中。因此，由于工作模式与数据采集方式的不同，给数据采集、信号处理以及后期的测井解释工作都带来了一定的困难。用于偏移成像的反射信号相对于高幅度井中直达波来说较为微弱，在准确记录反射信号的前提下，利用波场分离技术将反射波从全波形中提取出来，才能对反射信号进行成像处理。

反射声波测井的数据处理思路一般都借鉴了地震资料处理中较为成熟的滤波方法，或者直接使用地震处理软件，Tang在2007年，针对反射声波测井模型的特殊性，分析了反射波的传播特征，在参数估计法的基础上提出一种针对反射声波测井数据的线性预测滤波法。线性预测滤波法的思路为：首先利用已知的走时特征，根据各个井中直达波阵型的传播速度，构建传播矩阵，利用参数估计法估计得到直达波，随后从测井波形中减去井中直达波，剩余的波形即为反射波。

常规线性预测法中，以偶极反射声波测井为例，通常认为只有弯曲波一种直达波振型，构建传播矩阵时，使用的速度仅选取对应地层的横波速度。然而，采用上述方法构建传播矩阵时，存在以下问题：弯曲波是一种频散波，不同频率的弯曲波对应的传播速度是不同的，在低频时弯曲波的传播速度趋向于横波速度，高频时趋向于流体速度。因此，偶极反射声波测井中，如果假设只有一种直达波阵型，并且设定这种阵型的传播速度为横波速度，势必只能滤掉对应的弯曲波中以横波速度传播的部分，而以低于横波速度传播的弯曲波部分将留在剩余信号中，对井外微弱反射波的识别与提取形成严重干扰。因此，目前的反射声波测井成像结果中普遍存在成像精度低、可靠性较差、多解性强等问题。

技术实现要素：

本发明的目的在于，解决现有的反射声波测井技术中普遍存在成像精度低、可靠性较差、多解性强等问题，提供了一种反射声波测井的波场分离方法及装置，利用慢度时间相干法求出更准确的传播矩阵，基于频散效应考虑了直达波的多个阵型，能够稳定有效地从实测波形中提取出反射波，从而为反射声波测井提供更加准确有效的信息。

为了实现上述目的，一方面，本发明提供了一种反射声波测井的波场分离方法。该方法包括步骤：对实测波形进行预处理，并对预处理后的实测波形中的直达波的慢度范围进行计算；根据计算出的慢度范围构建传播矩阵；根据传播矩阵与预处理后的实测波形，反演得到直达波的波形；将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形，并对反射波波形进行偏移成像，获得最终的反射波波形。

优选地，慢度范围通过以下步骤获取：慢度时间相干法通过对阵列波形的相关函数的计算，得到阵列波形的到时与慢度。

优选地，相关函数的计算公式为：

其中，ρ(S,T)为相关函数，S为慢度，T为计算时窗起始位置，m为接收器的总个数，T_w为对比窗口长度，w_i(t)为第i道波形，δ为各接收器之间的间距。

优选地，根据计算出的慢度范围构建传播矩阵步骤具体为：根据计算出的慢度范围求出速度范围，根据速度范围构建传播矩阵。

优选地，传播矩阵为：

其中，N是接收器的总个数，n为不大于N的自然数，L为直达波阵型的总个数，l为不大于L的自然数，ω为角频率，d为相邻两个接收器之间的距离，V_l为直达波中第l个阵型的速度。

优选地，直达波的波形计算公式为：

其中，～表示复共轭，T表示转置，E为传播矩阵，A为直达波的阵型组成的矩阵，W为预处理后的实测波形组成的矩阵。

优选地，将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形步骤包括：将第n个接收器对应的直达波频谱从其对应全波形的频谱中移除，得到第n个接收器的剩余信号对应的频谱，

其中，W_n(ω)是第n个接收器接收到的实测波形，N为接收器的总个数，L为直达波阵型的总个数，A_l(ω)为直达波的第l个阵型，l为不大于L的自然数，R_n(ω)为第n个接收器的剩的信号对应的频谱。

优选地，将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形，并对反射波波形进行偏移成像，获得最终的反射波波形具体包括：将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形，提高反射波波形的信噪比，并对提高信噪比后的反射波波形进行偏移成像，获得最终的反射波波形。

另一方面，本发明提供了一种反射声波测井的波场分离装置，其特征在于，该装置包括：计算模块，用于对实测波形进行预处理，并对预处理后的实测波形中的直达波的慢度范围进行计算；构建模块，用于根据计算出的慢度范围构建传播矩阵；反演模块，根据传播矩阵与预处理后的实测波形，反演得到直达波的波形；偏移模块，用于将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形，并对反射波波形进行偏移成像，获得最终的反射波波形。

优选地，计算模块具体用于对实测波形进行预处理，并利用慢度时间相干法对预处理后的实测波形中的直达波的慢度范围进行计算。

本发明提供的一种反射声波测井的波场分离方法及装置，利用慢度时间相干法求出更准确的传播矩阵，基于频散效应考虑了直达波的多个阵型，能够稳定有效地从实测波形中提取出反射波，从而为反射声波测井提供更加准确有效的信息。

附图说明

图1为本发明实施例提供的一种反射声波测井的波场分离方法的流程示意图；

图2为为本发明实施例提供的一种反射声波测井的波场分离装置的结构示意图；

图3为本发明实施例提供的对X582m处的实测波形进行预处理后的波形图；

图4为图3中的波形对应的时间慢度相干图；

图5为对图3中的波形选取L＝1进行线性预测滤波后的波形图；

图6为图5中的波形对应的时间慢度相干图；

图7为对图3中的波形选取L＝5进行线性预测滤波后的波形图；

图8为对图3中的波形选取L＝4进行线性预测滤波后的波形图；

图9为对图3中的波形选取L＝3进行线性预测滤波后的波形图；

图10为图9中的波形对应的时间慢度相干图；

图11为本发明实施例提供的对X720m-X460m深度范围的实测波形预处理后的全波波形图；

图12为对图11中的波形选取L＝1进行线性预测滤波，沿时距曲线叠加增强和FK滤波后的波形图；

图13为对图11中的波形选取L＝4进行线性预测滤波，沿时距曲线叠加增强和FK滤波后的波形图。

具体实施方式

下面通过附图和实施例，对本发明的技术方案作进一步的详细描述。

图1为本发明实施例提供的一种反射声波测井的波场分离方法的流程示意图。如图1所示，该方法包括步骤201-204：

步骤201，对实测波形进行预处理，并对预处理后的实测波形中的直达波的慢度范围进行计算。

具体地，慢度范围通过以下步骤获取：慢度时间相干法通过对阵列波形的相关函数的计算，得到阵列波形的到时与慢度。

在本实施例中，采用Xmac-F1正交偶极子阵列声波测井仪器获取实测波形。该测井仪器源距为3.2004m，有八个四分量接收器，分别记录XX、YY、XY、YX的分量信号，八个接收器之间的间距为0.1524m。数据接收的深度范围为X720m到X460m，共有1707个深度点。波形的采样间隔为36us，采样时间点共400个，对应记录时间长度为14.4ms。

选取该测井仪器在深度X582m处测得的全波波形进行预处理，包括去增益和带通滤波。图3为本发明实施例提供的预处理后的全波波形图。如图3所示，图中有幅度较强的弯曲波。经过速度分析，该深度点对应的横波速度约为3000m/s，根据该速度绘出的时距曲线与弯曲波的到时有较高的匹配度。

慢度时间相干法是一种二维网格搜索法，可以找出阵列波形的相关函数极值区域对应的波的到时与慢度，是目前广泛使用的较为成熟的计算慢度的方法。相关函数的计算公式为：

其中，ρ(S,T)为相关函数，S为慢度，T为计算时窗起始位置，m为接收器的总个数，T_w为对比窗口长度，w_i(t)为第i道波形，δ为各接收器之间的间距。

利用慢度时间相干法对图3中的波形进行处理，处理结果如图4所示，图中有两个极值区域，在一个极值区域中，慢度S1＝330us/m，对应速度为3000m/s；在另一个极值区域中，慢度S2＝430us/m，对应速度为2300m/s。慢度为S1的极值区域对应横波速度的部分，慢度为S2的极值区域则对应由于频散造成的弯曲波中以低于横波速度传播的部分。由于频散，弯曲波对应的慢度为一个范围。

步骤202，根据计算出的慢度范围构建传播矩阵。

具体地，根据慢度范围S1-S2，求出对应的速度范围，构建传播矩阵E。与未考虑频散采用地层横波速度构建的传统传播矩阵相比，根据慢度范围构建的传播矩阵E更加接近实际情况。

步骤203，根据所述传播矩阵与所述预处理后的实测波形，反演得到所述直达波的波形。

具体地，利用参数估计法，根据所述传播矩阵与所述预处理后的实测波形，反演得到所述直达波的波形。

在参数估计法中，将直达波看作是以已知速度传播的模式，直达波的阵型表达式为A_l(ω)exp(-iωz/v_l)(l＝1,....,L)，其中ω为角频率，A_l(ω)为直达波的第l个阵型,v_l是第l种直达波振型对应的传播速度，z代表第n个接收器的所处位置,L是直达波的振型个数。例如，对单极数据，L＝3，直达波的阵型包括P波、S波、斯通利波；偶极数据L＝1，对应一种阵型弯曲波。假设该直达波从第n个接收器传至第m个接收器，在数学上可以把这一过程表示为A_l(ω)exp(-iω(m-n)d/v_l)，其中d是相邻两个接收器之间的距离。指数m与n之间大小关系的不同，对应该直达波的传播方向不同，m大于n意味着该直达波向前传播，反之m小于n意味着该直达波向后传播。基于上述理论，未知的直达波频谱A与阵列声波数据W之间的关系可用以下关系式表示：

其中，E_l＝exp(-iωd/v_l)，N是阵列接收器的个数；A_n(ω)(n＝1,…,L)为接收器n接收到的井中直达波的L个振型，为待求量；W_n(ω)(n＝1,…,N)是测井仪器接收到的实测波形。上式的矩阵形式为E·A＝W。

E矩阵为传播矩阵，该矩阵由井中直达波各个阵型的速度v_l决定。当L<N的时，通过最小二乘法求得A，计算公式为：

其中，～表示复共轭，T表示转置，E为所述传播矩阵，A为所述直达波的阵型组成的矩阵，W为所述预处理后的实测波形组成的矩阵。将估算得到的各个阵型的频谱进行相加，则得到了对应该接收器位置n处的直达波的完整频谱。

步骤204，将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形，并对反射波波形进行偏移成像，获得最终的反射波波形。

具体地，将反演得到的直达波的波形从预处理后的实测波形中移除，得到反射波波形，提高反射波波形的信噪比，并对提高信噪比后的反射波波形进行偏移成像，获得最终的反射波波形。

将第n个接收器对应的直达波频谱从其对应全波形的频谱中减去，得到的剩余的第n个接收器信号对应的频谱，

其中，W_n(ω)是第n个接收器接收到的实测波形，N为接收器的总个数，L为直达波阵型的总个数，A_l(ω)为所述直达波的第l个阵型，l为不大于L的自然数，R_n(ω)为所述第n个接收器的剩的信号对应的频谱。

将上式变换回时间域，得到对应于第n个接收器的、时域的、去除了直达波的剩余信号。对n＝1,…,N重复，则得到完整的阵列反射信号。

本实施例中，假设井中直达波中只有一种振型，为以横波速度传播的弯曲波。取横波速度为3000m/s，对预处理后的实测波形进行线性预测滤波，即步骤202-205中的操作。滤波结果如图5所示，以横波速度滤波后，时距曲线，即图中虚线标定的部分被过滤掉了，但是在时距曲线之后到达的波形，过滤效果不佳。如图6所示，仍存在慢度为S2的极值区域。

如图3所示，两个极值区域整体的上下边界分别为：上边界对应的慢度为250us/m，速度为4000m/s；下边界对应的慢度为500us/m，速度为2000m/s。将此范围分为5等份，即取L1＝5，即假设井中直达波有五种不同速度传播的阵型，五种速度分别为V1₁＝4000m/s，V1₂＝3500m/s,V1₃＝3000m/s,V1₄＝2500m/s,V1₅＝2000m/s。根据此速度范围，构建传播矩阵E1，对预处理后的实测波形进行线性预测滤波。如图7所示，与图3和图6相比，图7中的波形过滤效果得到明显的提高。另外，选取L2＝4，即假设井中直达波有四种不同速度传播的阵型，四种速度分别为V2₁＝4000m/s，V2₂＝3300m/s,V2₃＝2700m/s,V2₄＝2000m/s，构建传播矩阵E2，对预处理后的实测波形进行线性预测滤波，滤波结果如图8所示；选取L3＝3，即假设井中直达波有三种不同速度传播的阵型，三种速度分别为V3₁＝4000m/s，V3₂＝3000m/s,V3₃＝2000m/s，滤波结果如图9所示。

如图7、图8和图9所示，图3中的波形经过滤波后，弯曲波均受到较大的压制，且残余的弯曲波的幅度相近。由此可知，在线性滤波中，速度的选取以及直达波振型个数L的选取对滤波效果的影响较小，此方法具有较高的稳定性。如图10所示，图中已不存在极值区域，说明弯曲波被有效地压制。

利用倾角叠加技术与FK滤波技术，提高反射波波形的信噪比。如图7、图8和图9所示，在线性滤波后，由于残余的弯曲波，井壁不规则所产生的散射波和实测波形采集时产生的不可避免的数据噪声，仍有的残余的波形存在。对线性滤波后的波形沿时距曲线进行叠加增强处理、以及FK滤波，以进一步提升反射信号信噪比。

对提高信噪比后的反射波波形进行偏移成像。本实施例中，对所测井的整个深度范围，即X720m到X460m，重复步骤201-206。图11为本发明实施例提供的对X720m-X460m深度范围的实测波形预处理后的全波波形图；图12为对图11中的波形选取L＝1进行线性预测滤波，沿时距曲线叠加增强和FK滤波后的波形图；图13为对图11中的波形选取L＝4进行线性预测滤波，沿时距曲线叠加增强和FK滤波后的波形图。

如图11所示，图中弯曲波的幅度较大，井外反射波微弱，难以识别。如图12和图13所示，因考虑到频散效应采用了线性预测滤波方法，井中弯曲波被进一步的压制，井外反射信号更加清晰，图中可清晰地看见第一反射信号1，第二反射信号2，第三反射信号3。对于第二反射信号2，图12中由于与残余弯曲波幅度差别较大而显得较微弱，图13中，则在X660深度附近清晰可见下行反射信号2。这进一步说明了采用本发明提出的技术方案，可以获得更好的波场分离效果，进而提升反射声波测井的有效性。

本发明实施例提供的一种反射声波测井的波场分离方法，通过引入慢度时间相干法，从实测波形中求得直达波的慢度范围，进而求得更接近实际情况的传播矩阵。本发明实施例不仅考虑以横波速度传播的弯曲波，还考虑了直达波以整个慢度范围传播的各个振型。基于频散效应，结合慢度时间相干法的对预处理后的实测声波进行处理，较常规方法而言，实测波形中的弯曲波被进一步压制，使得微弱的反射波可以被有效提取。并且，本发明实施例中振型个数以及对应速度的选取，对滤波效果的影响较小。由此可知，该方法具有较好的稳定性，适于用在实际处理中。

图2为本实施例提供的一种反射声波测井的波场分离装置的结构示意图。如图2所述，该装置400包括：

计算模块401用于对实测波形进行预处理，并对预处理后的实测波形中的直达波的慢度范围进行计算。计算模块具体用于对实测波形进行预处理，并利用慢度时间相干法对预处理后的实测波形中的直达波的慢度范围进行计算。

构建模块402用于根据计算出的慢度范围构建传播矩阵。

反演模块403根据所述传播矩阵与所述预处理后的实测波形，反演得到所述直达波的波形。

偏移模块404用于将所述反演得到的直达波的波形从所述预处理后的实测波形中移除，得到反射波波形，并对所述反射波波形进行偏移成像，获得最终的反射波波形。

本实施例提供的一种反射声波测井的波场分离装置，通过计算模块利用慢度时间相干法求出直达波的慢度范围，再用构建模块根据慢度范围构建更准确的传播矩阵，反演模块基于频散效应考虑了直达波的多个阵型从而使偏移模块稳定有效地从实测波形中提取出反射波，为反射声波测井提供更加准确有效的信息。

以上的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。