CN106547029A - 一种基于测井资料的线性滑动理论参数反演方法及装置 - Google Patents

Abstract

本发明实施例提供一种基于测井资料的线性滑动理论参数反演方法及装置。所述方法包括：获取地质资料及测井资料，测井资料包括实际声波全波列波形；对实际声波全波列波形进行滤波，分离出实际斯通利波波形；建立实际斯通利波波形与理论斯通利波波形误差的目标函数；根据预设最优化方法，利用实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、测井资料、地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。该方法利用测井资料反演得到了线性滑动理论参数，有效地避免了现有技术中利用岩石物理实验求取该理论参数时的困难。

Description

一种基于测井资料的线性滑动理论参数反演方法及装置

技术领域

本发明涉及地球物理反演技术领域，具体而言，涉及一种基于测井资料的线性滑动理论参数反演方法及装置。

背景技术

当地层介质中含有裂缝时，其速度和密度会呈现方向性，也就是各向异性性质。各向异性基本理论主要有三种：Hudson理论、Thomson理论和线性滑动理论。Hudson理论和Thomsen理论对裂缝尺度有明确规定，并且要求裂缝在地震波长尺度内均匀分布。因此Hudson理论和Thomsen理论在应用中往往对应等效介质模型。但在实际地层中，往往会发育单条或多组任意尺度的裂缝。对于这类裂缝，Hudson理论和Thomsen理论的应用受到限制。而线性滑动理论对裂缝尺度没有做要求，可以用来描述单条或多组任意尺度裂缝模型，进一步地还可用来指示裂缝中所含流体性质，因此线性滑动理论参数的求取就具有重要的理论及实际意义。传统的求取线性滑动理论参数的方法是岩石物理实验，然而岩石物理实验的不足之处在于岩心样品少，实验费用成本高，同时地下地层的真实条件很难模拟。

利用测井资料可以通过不同的方法反演得到地层的纵波、快慢横波、斯通利波的速度，岩石的渗透率和裂缝的反射系数等地质参数。然而目前还没有利用测井资料求取线性滑动理论参数的方法。

发明内容

有鉴于此，本发明实施例的目的在于提供一种基于测井资料的线性滑动理论参数反演方法及装置，以解决上述问题。

为了实现上述目的，本发明实施例采用的技术方案如下：

第一方面，本发明实施例提供了一种基于测井资料的线性滑动理论参数反演方法，所述方法包括：获取地质资料及测井资料，所述测井资料包括实际声波全波列波形；对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形；建立实际斯通利波波形与理论斯通利波波形误差的目标函数；根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。

第二方面，本发明实施例提供了一种基于测井资料的线性滑动理论参数反演装置，所述装置包括：接收模块，用于获取地质资料及测井资料，所述测井资料包括实际声波全波列波形；第一处理模块，用于对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形；第二处理模块，用于建立实际斯通利波波形与理论斯通利波波形误差的目标函数；根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。

与现有技术相比，本发明实施例提供的一种基于测井资料的线性滑动理论参数反演方法及装置，通过获取测井及地质资料，从测井资料中的实际声波全波列波形里分离出实际斯通利波波形，并建立实际斯通利波波形与理论斯通利波波形误差的目标函数，最终根据预设最优化方法，利用所述实际斯通利波波形、以及对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的理论斯通利波波形，获取到当目标函数满足预设条件时的线性滑动理论参数。所述方法利用测井资料反演得到了线性滑动理论参数，有效地避免了现有技术中利用岩石物理实验求取该理论参数时岩心样品少，实验费用成本高，同时地下地层的真实条件很难模拟的困难。

本发明的其他特征和优点将在随后的说明书阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明实施例了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1是本发明实施例提供的服务器的结构示意图。

图2是本发明第一实施例提供的一种基于测井资料的线性滑动理论参数反演方法的流程图。

图3是本发明第一实施例提供的一种基于测井资料的线性滑动理论参数反演方法中步骤S320的一种详细流程图。

图4是本发明第一实施例提供的一种基于测井资料的线性滑动理论参数反演方法中步骤S340的一种详细流程图。

图5本发明第一实施例提供的一种基于测井资料的线性滑动理论参数反演方法中获取所述目标函数的梯度近似值的方法的流程图。

图6是本发明第二实施例提供的一种基于测井资料的线性滑动理论参数反演方法的流程图。

图7是本发明第二实施例提供的一种基于测井资料的线性滑动理论参数反演方法中步骤S420的一种详细流程图。

图8是本发明第二实施例提供的一种基于测井资料的线性滑动理论参数反演方法中步骤S440的一种详细流程图。

图9本发明第二实施例提供的一种基于测井资料的线性滑动理论参数反演方法中获取所述目标函数的梯度近似值的方法的流程图。

图10是本发明第三实施例提供的一种基于测井资料的线性滑动理论参数反演装置500的结构框图。

图11是本发明第三实施例提供的一种基于测井资料的线性滑动理论参数反演装置500中第一处理模块520的一种详细结构框图。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

本发明实施例提供的基于测井资料的线性滑动理论参数反演方法可以应用于服务器中。图1示出了服务器100的结构示意图，请参阅图1，所述服务器100包括存储器110、处理器120以及网络模块130。

存储器110可用于存储软件程序以及模块，如本发明实施例中的基于测井资料的线性滑动理论参数反演方法及装置对应的程序指令/模块，处理器120通过运行存储在存储器110内的软件程序以及模块，从而执行各种功能应用以及数据处理，即实现本发明实施例中的基于测井资料的线性滑动理论参数反演方法。存储器110可包括高速随机存储器，还可包括非易失性存储器，如一个或者多个磁性存储装置、闪存、或者其他非易失性固态存储器。进一步地，上述存储器110内的软件程序以及模块还可包括：操作系统111以及服务模块112。其中操作系统111，例如可为LINUX、UNIX、WINDOWS，其可包括各种用于管理系统任务(例如内存管理、存储设备控制、电源管理等)的软件组件和/或驱动，并可与各种硬件或软件组件相互通讯，从而提供其他软件组件的运行环境。服务模块112运行在操作系统111的基础上，并通过操作系统111的网络服务监听来自网络的请求，根据请求完成相应的数据处理，并返回处理结果给客户端。也就是说，服务模块112用于向客户端提供网络服务。

网络模块130用于接收以及发送网络信号。上述网络信号可包括无线信号或者有线信号。

可以理解，图1所示的结构仅为示意，服务器100还可包括比图1中所示更多或者更少的组件，或者具有与图1所示不同的配置。图1中所示的各组件可以采用硬件、软件或其组合实现。另外，本发明实施例中的服务器还可以包括多个具体不同功能的服务器。

图2示出了本发明第一实施例提供的一种基于测井资料的线性滑动理论参数反演方法的流程图，请参阅图2，本实施例描述的是服务器的处理流程，所述方法包括：

步骤S310，获取地质资料及测井资料，所述测井资料包括实际声波全波列波形。

其中，测井资料是将各种仪器放入井内，沿井身测量到的钻井地质剖面上地层的各种物理参数。所述测井资料包括声波测井资料、井径测井数据、自然伽马测井数据、密度测井数据、微电阻率成像测井数据。其中，所述声波测井资料包括所述实际声波全波列波形，所述实际声波全波列波形增益数据及所述实际声波全波列波形记录开始时间数据。声波全波列波形由多种波形成分构成，包括滑行纵波、滑行横波、伪瑞利波、斯通利波等。根据所述自然伽马测井数据、密度测井数据可以确定介质的密度、速度、孔隙度参数，根据所述微电阻率成像测井数据和声波测井资料可以确定裂缝深度位置和裂缝类型。地质资料是指以文字、图表、声像、标本、样品、岩矿心等不同形式存在的地质、矿产信息或实物。所述地质资料包括目标井地层层位信息及岩性信息。

步骤S320，对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形。

作为一种具体的实施方式，请参阅图3，步骤S320可以包括：

步骤S321，根据所述测井资料，采用阈值法或人机交互的方法进行斯通利波时差提取，获取斯通利波时差。

步骤S322，根据所述测井资料及所述斯通利波时差，对所述实际声波全波列波形进行频谱分析，确定实际斯通利波波形的频率范围。

其中，对所述实际声波全波列波形进行频谱分析可以得到所述实际声波全波列波形中各成分波形的频率范围，因此所述实际斯通利波的频率范围也可以确定。

步骤S323，对所述实际声波全波列波形进行滤波，分离出所述频率范围对应的实际斯通利波波形。

作为一种具体的实施方式，可以采用滤波器对所述实际声波全波列波形进行滤波，分离出所述频率范围对应的实际斯通利波波形。例如，可以采用余弦镶边带通滤波器，其频谱为：

其中，f₁和f₂之间是加余弦窗的低通过渡带，f₂为低截频率，f₃为高截频率，f₃和f₄之间是加余弦窗的高通过渡带。根据所述实际斯通利波波形的频率范围确定所述采用余弦镶边带通滤波器的参数f₁、f₂、f₃及f₄，将所述实际声波全波列波形通过所述余弦镶边带通滤波器，即可分离处所述实际斯通利波波形。

优选的，在步骤S320之前，所述方法还可以包括：

对所述测井资料进行预处理，所述预处理包括井眼环境校正、数据滤波、坏道剔除或数据缺失值填补。

通过对所述测井资料进行预处理，可以有效地提高所述测井资料中所述实际声波全波列波形以及其他后续所需的数据的质量。

步骤S330，建立实际斯通利波波形与理论斯通利波波形误差的目标函数。

所述目标函数可以为：

E(m)＝||Af(m)-d||²＝(Af(m)-d)²

其中，E(m)为所述目标函数；m为线性滑动理论参数，m＝m(S_N,S_T)，S_N和S_T为正向屈服度系数以及切向屈服度系数；d为所述实际斯通利波波形；f(m)为理论声波全波列波形，A为与滤波器对应的滤波算子，Af(m)为对所述理论声波全波列波形f(m)进行滤波，分离出的所述理论斯通利波波形。

步骤S340，根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。

其中，最优化方法是指为了达到最优化目的所提出的各种求解方法。从数学意义上说，最优化方法是一种求极值的方法，即在一组约束为等式或不等式的条件下，使系统的目标函数达到极值，即最大值或最小值。所述预设最优化方法可以为共轭梯度最优化方法。当所述最优化方法为共轭梯度最优化方法时，所述预设条件可以为所述目标函数的梯度近似值的绝对值小于或等于预设精度值。所述线性滑动理论参数可以为正向屈服度系数以及切向屈服度系数。所述预设模型可以包括含孔隙的流体介质理论模型及线性滑动理论裂缝模型，优选的，所述含孔隙的流体介质理论模型为Biot理论模型。

作为一种具体的实施方式，请参阅图4，步骤S340可以包括：

步骤S341，根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型及所述实际斯通利波波形，获取所述目标函数的梯度近似值。

作为一种具体的实施方式，请参阅图5，所述获取所述目标函数的梯度近似值的步骤可以包括：

步骤S351，根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型，获得第一理论声波全波列波形及第二理论声波全波列波形。

步骤S352，对所述第一理论声波全波列波形及所述第二理论声波全波列波形分别进行滤波，分离出对应的第一理论斯通利波波形及第二理论斯通利波波形。

步骤S353，根据所述实际斯通利波波形、所述第一理论斯通利波波形、所述第二理论斯通利波波形及预设增量值，获取所述目标函数的梯度近似值。

下面对步骤S351至步骤S353进行详细说明。

利用Biot理论模型可以模拟孔隙地层中的测井声波传播。由Biot给出的饱和流体孔隙介质弹性波传播理论，可以得到运动方程：

其中，i、j表示x,y,z三个不同方向分量，u_j和U_j分别表示固相和流相的位移在j方向的分量，b_ij表示流体相对固体骨架运动时的耗散系数，σ_ij'j表示作用在固相上的应力分量在j方向上的偏导，s表示作用在流体单元侧面上的应力。

若φ表示孔隙度，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂表示视质量，即流相相对固相运动时的质量耦合系数。三者与固体密度ρ_s和流体密度ρ_f之间满足关系式：

ρ₁₁+ρ₁₂＝(1-φ)ρ_s

ρ₁₂+ρ₂₂＝φρ_f

柯西方程：

其中，e_ij表示固相应变张量，H表示流相应变张量。

本构方程：

其中，σ表示固相应力，s表示流相应力，C表示含裂缝介质的弹性系数矩阵。

利用线性滑动理论可以模拟含裂缝地层中的声波传播。在弱各向异性假设下，认为含裂缝介质的弹性系数矩阵C由各向同性背景介质弹性系数矩阵C^iso和各向异性扰动C^ani组成，即:

C＝C^iso-C^ani

其中，C^iso和C^ani的表达式为：

其中，Q_i表示固相和流相之间体积变化的耦合系数，R表示描述孔隙流体的弹性参数，Δ_N和Δ_T表示正向和切向弱度，取值从0至1，且有

其中，S_N和S_T为正向屈服度系数和切向屈服度系数。

结合上述柯西方程、本构方程和运动方程，可以得到双相介质二维三分量一阶速度—应力方程。假设固相介质速度向量为ν＝(υ_x,υ_y,υ_z)^T，流相介质速度向量为V＝(V_x,V_y,V_z)^T，令则速度和应力分量分别表示为：

其中，所述双相介质二维三分量一阶速度—应力方程中的速度方程构成了速度场，根据测井仪器的结构参数设置声源与接收器的位置关系后，即可得到对应的理论声波全波列波形。

以上各式中，ρ₁₁、ρ₁₂、ρ₂₂、b_ij、R、Q_i、μ、λ、φ、ρ_s、ρ_f均为已知参数，其中，μ、λ、φ、ρ₁₁、ρ₁₂、ρ₂₂、ρ_s、ρ_f可以通过所述测井资料获得，b_ij、R、Q_i可以通过所述地质资料获得，S_N和S_T为未知参数。

首先，假设所述预设线性滑动理论参数为m_k∈Rⁿ，m_k＝m(S_N,k,S_T,k)，k＝0,1,…，所述预设增量值为Δm，Δm＝Δm(ΔS_N,ΔS_T)。ρ₁₁、ρ₁₂、ρ₂₂、b_ij、R、Q_i、μ、λ、φ、ρ_s、ρ_f均为已知参数，将预设线性滑动理论参数m_k与m_k+Δm分别代入所述本构方程，分别获得所述双相介质二维三分量一阶速度—应力方程中的速度方程构成的速度场后，根据实际测井时的声源与接收器的位置关系，即可获得第一理论声波全波列波形f(m_k)及第二理论声波全波列波形f(m_k+Δm)。

然后，根据已知的所述实际声波全波列波形d、滤波算子A，将所述第一理论声波全波列波形f(m_k)及所述第二理论声波全波列波形f(m_k+Δm)分别代入所述目标函数E(m)＝||Af(m)-d||²＝(Af(m)-d)²中，获得对应的目标函数值E(m_k)以及E(m_k+Δm)。

最后，采用差分近似微分形式，即获得与所述预设线性滑动理论参数m_k对应的所述目标函数的梯度近似值

步骤S342，若所述目标函数的梯度近似值的绝对值小于或等于预设精度值，则所述预设线性滑动理论参数为获取的线性滑动理论参数；若所述目标函数的梯度近似值的绝对值大于所述预设精度值，则根据预设规则对所述预设线性滑动理论参数进行更新，重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数。

假设所述预设精度值为ε，其中ε＞0。若则所述预设线性滑动理论参数m_k为获取的线性滑动理论参数；若则根据预设规则对所述预设线性滑动理论参数进行更新，通过线性搜索计算步长α_k，令m_k+1＝m_k+α_kt_k，根据得到t_k，其中 根据更新后的所述预设线性滑动理论参数m_k+1重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数。

所述获取的线性滑动理论参数可以用于含裂缝各向异性地层的地震叠前反演技术，可以用来描述单条或多组任意尺度裂缝模型，指示裂缝中流体性质，因此具有重要的理论及实际意义。

本发明第一实施例提供的基于测井资料的线性滑动理论参数反演方法，通过获取测井及地质资料，从测井资料中的实际声波全波列波形里分离出实际斯通利波波形，并建立实际斯通利波波形与理论斯通利波波形误差的目标函数，最终根据预设最优化方法，利用所述实际斯通利波波形、以及对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的理论斯通利波波形，获取到当目标函数满足预设条件时的线性滑动理论参数。所述方法利用测井资料反演得到了线性滑动理论参数，有效地避免了现有技术中利用岩石物理实验求取该理论参数时岩心样品少，实验费用成本高，同时地下地层的真实条件很难模拟的困难，同时节约了经济成本。

图6示出了本发明第二实施例提供的一种基于测井资料的线性滑动理论参数反演方法的流程图，请参阅图6，本实施例描述的是服务器的处理流程，所述方法包括：

步骤S410，获取地质资料及测井资料，所述测井资料包括由多个实际声波全波列波形构成的实际阵列声波全波列波形。

其中，测井资料是将各种仪器放入井内，沿井身测量到的钻井地质剖面上地层的各种物理参数。所述测井资料包括阵列声波测井资料、井径测井数据、自然伽马测井数据、密度测井数据、微电阻率成像测井数据。其中，所述阵列声波测井资料包括所述由多个实际声波全波列波形构成的实际阵列声波全波列波形，所述实际阵列声波全波列波形增益数据及所述实际阵列声波全波列波形记录开始时间数据。声波全波列波形由多种波形成分构成，包括滑行纵波、滑行横波、伪瑞利波、斯通利波等。根据所述自然伽马测井数据、密度测井数据可以确定介质的密度、速度、孔隙度参数，根据所述微电阻率成像测井数据和阵列声波测井资料可以确定裂缝深度位置和裂缝类型。地质资料是指以文字、图表、声像、标本、样品、岩矿心等不同形式存在的地质、矿产信息或实物。所述地质资料包括目标井地层层位信息及岩性信息。

其中，所述多个实际声波全波列可以为在井轴上设置的由N个接收器构成的阵列仪接收到的信号。

步骤S420，对所述实际阵列声波全波列波形中的多个实际声波全波列波形分别进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形构成实际阵列斯通利波波形。

作为一种具体的实施方式，请参阅图7，步骤S420可以包括：

步骤S421，根据所述测井资料及采用慢度-时间相关法，获取斯通利波时差。

声波测井信号处理就是利用各种信号处理方法从声波全波列波形中对各成分波进行识别，并提取各成分波的速度和幅度等信息。声波测井信号处理的方法主要从时域和频域出发，采用数字信号分析技术对声波全波列波形进行分析处理获取各波成分的信息。在时域中比较常用的方法有阀值法、差值检测法、慢度-时间相关法和n次方根法等，在频域中有互功率谱相位法、协方差分析法、最大似然法等。其中，所述慢度—时间相关法相对于阀值法、差值检测法有更高的准确性，而相对于n次方根法、互功率谱相位法、协方差分析法、最大似然法有更快的速度，因此也是目前普遍使用的时差提取方法。

下面对采用慢度-时间相关法获取斯通利波时差进行详细说明。

假设在井轴上有N个接收器的阵列仪，对于阵列声波测井资料相关函数定义为：

其中，ρ表示相关系数，X_m(t)为N个接收器中第m个接收器的波形，d为接收器间距，T_w为时间窗长度，s为慢度，T为时间。

对于某一时间T和慢度值s，相关系数由上式计算，在STC二维谱上手动选择一个合适的矩形窗口，搜索窗口最大值，找到使相关函数ρ取极大值的T和s值，即为某一成分波(包括纵波、横波、斯通利波)的首波到达时间和慢度。通过该方法可以获得斯通利波的首播到达时间和慢度，即斯通利波时差。

由多个实际声波全波列波形构成的实际阵列声波全波列波形是反映井筒内外介质特征的、包含各种波成分的完整波列，通过多个接收器接收到的所述多个实际声波全波列波形可以提取到更准确的斯通利波时差。

步骤S422，根据所述测井资料及所述斯通利波时差，对所述实际阵列声波全波列波形进行频谱分析，确定实际斯通利波波形的频率范围。

步骤S423，对所述实际阵列声波全波列波形中的多个实际声波全波列波形进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形与所述频率范围对应，所述多个实际斯通利波波形构成实际阵列斯通利波波形。

作为一种具体的实施方式，可以采用滤波器对所述实际阵列声波全波列波形中的多个实际声波全波列波形进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形。例如，可以采用余弦镶边带通滤波器，其频谱为：

其中，f₁和f₂之间是加余弦窗的低通过渡带，f₂为低截频率，f₃为高截频率，f₃和f₄之间是加余弦窗的高通过渡带。根据所述实际斯通利波波形的频率范围确定所述采用余弦镶边带通滤波器的参数f₁、f₂、f₃及f₄，将所述实际声波全波列波形通过所述余弦镶边带通滤波器，即可分离处所述实际斯通利波波形。

优选的，在步骤S420之前，所述方法还可以包括：

对所述测井资料进行预处理，所述预处理包括井眼环境校正、数据滤波、坏道剔除或数据缺失值填补。

通过对所述测井资料进行预处理，可以有效地提高所述测井资料中所述实际声波全波列波形以及其他后续所需的数据的质量。

步骤S430，建立所述实际阵列斯通利波波形与理论阵列斯通利波波形误差的目标函数。

所述目标函数可以为：

其中，E(m)为所述目标函数；m为线性滑动理论参数，m＝m(S_N,S_T)，S_N和S_T为正向屈服度系数以及切向屈服度系数；n为接收器的个数；d_i为第i个接收器接收到的所述实际斯通利波波形；f_i(m)为第i个接收器对应的理论声波全波列波形，A为与滤波器对应的滤波算子，Af_i(m)为对第i个接收器对应的理论声波全波列波形f(m)进行滤波，分离出的第i个接收器对应的理论斯通利波波形。

步骤S440，根据预设最优化方法，利用所述实际阵列斯通利波波形及理论阵列斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论阵列斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论阵列声波全波列波形进行滤波后得到的阵列波形。

其中，最优化方法是指为了达到最优化目的所提出的各种求解方法。从数学意义上说，最优化方法是一种求极值的方法，即在一组约束为等式或不等式的条件下，使系统的目标函数达到极值，即最大值或最小值。所述预设最优化方法可以为共轭梯度最优化方法。当所述最优化方法为共轭梯度最优化方法时，所述预设条件可以为所述目标函数的梯度近似值的绝对值小于或等于预设精度值。所述线性滑动理论参数可以为正向屈服度系数以及切向屈服度系数。所述预设模型可以包括含孔隙的流体介质理论模型及线性滑动理论裂缝模型，优选的，所述含孔隙的流体介质理论模型为Biot理论模型。

作为一种具体的实施方式，请参阅图8，步骤S440可以包括：

步骤S441，根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型及所述实际阵列斯通利波波形，获取所述目标函数的梯度近似值。

作为一种具体的实施方式，请参阅图9，所述获取所述目标函数的梯度近似值的步骤可以包括：

步骤S451，根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型，获得第一理论阵列声波全波列波形及第二理论阵列声波全波列波形。

步骤S452，对所述第一理论阵列声波全波列波形中的多个第一理论声波全波列波形进行滤波，分离出与所述多个第一理论声波全波列波形一一对应的多个第一理论斯通利波波形，所述多个第一理论斯通利波波形构成第一理论阵列斯通利波波形；对所述第二理论阵列声波全波列波形中的多个第二理论声波全波列波形进行滤波，分离出与所述多个第二理论声波全波列波形一一对应的多个第二理论斯通利波波形，所述多个第二理论斯通利波波形构成第二理论阵列斯通利波波形。

步骤S453，根据所述实际阵列斯通利波波形、所述第一理论阵列斯通利波波形、所述第二理论阵列斯通利波波形及预设增量值，获取所述目标函数的梯度近似值。

下面对步骤S451至步骤S453进行详细说明。

利用Biot理论模型可以模拟孔隙地层中的测井声波传播。由Biot给出的饱和流体孔隙介质弹性波传播理论，可以得到运动方程：

其中，i、j表示x,y,z三个不同方向分量，u_j和U_j分别表示固相和流相的位移在j方向的分量，b_ij表示流体相对固体骨架运动时的耗散系数，σ_ij'j表示作用在固相上的应力分量在j方向上的偏导，s表示作用在流体单元侧面上的应力。

若φ表示孔隙度，ρ₁₁和ρ₂₂分别表示介质单位体积内固相和流相部分的有效质量，ρ₁₂表示视质量，即流相相对固相运动时的质量耦合系数。三者与固体密度ρ_s和流体密度ρ_f之间满足关系式：

ρ₁₁+ρ₁₂＝(1-φ)ρ_s

ρ₁₂+ρ₂₂＝φρ_f

柯西方程：

其中，e_ij表示固相应变张量，H表示流相应变张量。

本构方程：

其中，σ表示固相应力，s表示流相应力，C表示含裂缝介质的弹性系数矩阵。

利用线性滑动理论可以模拟含裂缝地层中的声波传播。在弱各向异性假设下，认为含裂缝介质的弹性系数矩阵C由各向同性背景介质弹性系数矩阵C^iso和各向异性扰动C^ani组成，即:

C＝C^iso-C^ani

其中，C^iso和C^ani的表达式为：

其中，Q_i表示固相和流相之间体积变化的耦合系数，R表示描述孔隙流体的弹性参数，Δ_N和Δ_T表示正向和切向弱度，取值从0至1，且有

其中，S_N和S_T为正向屈服度系数和切向屈服度系数。

结合上述柯西方程、本构方程和运动方程，可以得到双相介质二维三分量一阶速度—应力方程。假设固相介质速度向量为ν＝(υ_x,υ_y,υ_z)^T，流相介质速度向量为V＝(V_x,V_y,V_z)^T，令则速度和应力分量分别表示为：

其中，所述双相介质二维三分量一阶速度—应力方程中的速度方程构成了速度场，根据测井仪器的结构参数设置声源与接收器的位置关系后，即可得到对应的理论声波全波列波形。

以上各式中，ρ₁₁、ρ₁₂、ρ₂₂、b_ij、R、Q_i、μ、λ、φ、ρ_s、ρ_f均为已知参数，其中，μ、λ、φ、ρ₁₁、ρ₁₂、ρ₂₂、ρ_s、ρ_f可以通过所述测井资料获得，b_ij、R、Q_i可以通过所述地质资料获得，S_N和S_T为未知参数。

首先，假设所述预设线性滑动理论参数为m_k∈Rⁿ，m_k＝m(S_N,k,S_T,k)，k＝0,1,…，所述预设增量值为Δm，Δm＝Δm(ΔS_N,ΔS_T)。ρ₁₁、ρ₁₂、ρ₂₂、b_ij、R、Q_i、μ、λ、φ、ρ_s、ρ_f均为已知参数，将预设线性滑动理论参数m_k与m_k+Δm分别代入所述本构方程，分别获得所述双相介质二维三分量一阶速度—应力方程中的速度方程构成的速度场后，根据实际测井时接收器的个数、声源与接收器的位置关系，即可获得第一理论阵列声波全波列波形f_i(m_k),i＝1,2,…,n及第二理论声波全波列波形f_i(m_k+Δm),i＝1,2,…,n。

然后，根据已知的所述实际阵列声波全波列波形d_i,i＝1,2,…,n、滤波算子A，将所述第一理论阵列声波全波列波形f_i(m_k),i＝1,2,…,n及所述第二理论阵列声波全波列波形f_i(m_k+Δm),i＝1,2,…,n分别代入所述目标函数中，获得对应的目标函数值E(m_k)以及E(m_k+Δm)。

最后，采用差分近似微分形式，即获得与所述预设线性滑动理论参数m_k对应的所述目标函数的梯度近似值

步骤S442，若所述目标函数的梯度近似值的绝对值小于或等于预设精度值，则所述预设线性滑动理论参数为获取的线性滑动理论参数；若所述目标函数的梯度近似值的绝对值大于所述预设精度值，则根据预设规则对所述预设线性滑动理论参数进行更新，重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数。

假设所述预设精度值为ε，其中ε＞0。若则所述预设线性滑动理论参数m_k为获取的线性滑动理论参数；若则根据预设规则对所述预设线性滑动理论参数进行更新，通过线性搜索计算步长α_k，令m_k+1＝m_k+α_kt_k，根据得到t_k，其中 根据更新后的所述预设线性滑动理论参数m_k+1重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数。

所述获取的线性滑动理论参数可以用于含裂缝各向异性地层的地震叠前反演技术，可以用来描述单条或多组任意尺度裂缝模型，指示裂缝中流体性质，因此具有重要的理论及实际意义。

本发明第二实施例提供的基于测井资料的线性滑动理论参数反演方法，通过获取测井及地质资料，从测井资料中的实际阵列声波全波列波形里分离出实际阵列斯通利波波形，并建立实际阵列斯通利波波形与理论阵列斯通利波波形误差的目标函数，最终根据预设最优化方法，利用所述实际阵列斯通利波波形、以及对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论阵列声波全波列波形进行滤波后得到的理论阵列斯通利波波形，获取到当目标函数满足预设条件时的线性滑动理论参数。所述方法利用测井资料反演得到了线性滑动理论参数，无需推导通过阵列声波测井资料求取线性滑动理论参数的公式，同时有效地避免了现有技术中利用岩石物理实验求取该理论参数时岩心样品少、实验费用成本高、地下地层的真实条件很难模拟的困难，节约了经济成本。此外，由于充分利用了阵列声波测井资料的高分辨率和对地层真实响应的优点，使得获取到的线性滑动理论参数更加可靠。

图10是本发明第三实施例提供的一种基于测井资料的线性滑动理论参数反演装置500的结构框图，请参阅图10，所述基于测井资料的线性滑动理论参数反演装置500包括接收模块510、第一处理模块520及第二处理模块530。

所述接收模块510，用于获取地质资料及测井资料，所述测井资料包括实际声波全波列波形。

作为一种具体的实施方式，所述接收模块510，具体用于获取地质资料及测井资料，所述测井资料包括由多个实际声波全波列波形构成的实际阵列声波全波列波形。

所述第一处理模块520，用于对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形。

作为一种具体的实施方式，所述第一处理模块520，具体用于对所述实际阵列声波全波列波形中的多个实际声波全波列波形分别进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形构成实际阵列斯通利波波形。

作为一种具体的实施方式，请参阅图11，所述第一处理模块520可以包括时差提取模块521、频谱分析模块522及滤波模块523。

所述时差提取模块521，用于根据所述测井资料，采用阈值法或人机交互的方法进行斯通利波时差提取，获取斯通利波时差。

作为一种具体的实施方式，所述时差提取模块521，具体用于根据所述测井资料及采用慢度-时间相关法，获取斯通利波时差。

所述频谱分析模块522，用于根据所述测井资料及所述斯通利波时差，对所述实际声波全波列波形进行频谱分析，确定实际斯通利波波形的频率范围。

作为一种具体的实施方式，所述频谱分析模块522，具体用于根据所述测井资料及所述斯通利波时差，对所述实际阵列声波全波列波形进行频谱分析，确定实际斯通利波波形的频率范围。

所述滤波模块523，用于对所述实际声波全波列波形进行滤波，分离出所述频率范围对应的实际斯通利波波形。

作为一种具体的实施方式，所述滤波模块523，具体用于对所述实际阵列声波全波列波形中的多个实际声波全波列波形进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形与所述频率范围对应，所述多个实际斯通利波波形构成实际阵列斯通利波波形。

优选的，所述基于测井资料的线性滑动理论参数反演装置500还可以包括预处理模块540，所述预处理模块540用于在对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形之前，对所述测井资料进行预处理，所述预处理包括井眼环境校正、数据滤波、坏道剔除或数据缺失值填补。

所述第二处理模块530，用于建立实际斯通利波波形与理论斯通利波波形误差的目标函数；根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。

作为一种具体的实施方式，所述第二处理模块530，具体用于根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型，获得第一理论声波全波列波形及第二理论声波全波列波形；对所述第一理论声波全波列波形及所述第二理论声波全波列波形分别进行滤波，分离出对应的第一理论斯通利波波形及第二理论斯通利波波形；根据所述实际斯通利波波形、所述第一理论斯通利波波形、所述第二理论斯通利波波形及预设增量值，获取所述目标函数的梯度近似值。若所述目标函数的梯度近似值的绝对值小于或等于预设精度值，则所述预设线性滑动理论参数为获取的线性滑动理论参数；若所述目标函数的梯度近似值的绝对值大于所述预设精度值，则根据预设规则对所述预设线性滑动理论参数进行更新，重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数。

所述第二处理模块530，具体用于建立所述实际阵列斯通利波波形与理论阵列斯通利波波形误差的目标函数；根据预设最优化方法，利用所述实际阵列斯通利波波形及理论阵列斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论阵列斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论阵列声波全波列波形进行滤波后得到的阵列波形。

作为一种具体的实施方式，所述第二处理模块530，具体用于根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型，获得第一理论阵列声波全波列波形及第二理论阵列声波全波列波形；对所述第一理论阵列声波全波列波形中的多个第一理论声波全波列波形进行滤波，分离出与所述多个第一理论声波全波列波形一一对应的多个第一理论斯通利波波形，所述多个第一理论斯通利波波形构成第一理论阵列斯通利波波形；对所述第二理论阵列声波全波列波形中的多个第二理论声波全波列波形进行滤波，分离出与所述多个第二理论声波全波列波形一一对应的多个第二理论斯通利波波形，所述多个第二理论斯通利波波形构成第二理论阵列斯通利波波形；根据所述实际阵列斯通利波波形、所述第一理论阵列斯通利波波形、所述第二理论阵列斯通利波波形及预设增量值，获取所述目标函数的梯度近似值。若所述目标函数的梯度近似值的绝对值小于或等于预设精度值，则所述预设线性滑动理论参数为获取的线性滑动理论参数；若所述目标函数的梯度近似值的绝对值大于所述预设精度值，则根据预设规则对所述预设线性滑动理论参数进行更新，重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数。

以上各模块可以是由软件代码实现，此时，上述的各模块可存储于服务器100的存储器110内。以上各模块同样可以由硬件例如集成电路芯片实现。

本发明实施例所提供的基于测井资料的线性滑动理论参数反演装置500，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述方法实施例中相应内容。

在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

另外，在本发明各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。

所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

Claims (10)

1.一种基于测井资料的线性滑动理论参数反演方法，其特征在于，所述方法包括：

获取地质资料及测井资料，所述测井资料包括实际声波全波列波形；

对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形；

建立实际斯通利波波形与理论斯通利波波形误差的目标函数；

根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。

2.根据权利要求1所述的方法，其特征在于，所述预设最优化方法为共轭梯度最优化方法，所述根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形，具体包括：

根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型及所述实际斯通利波波形，获取所述目标函数的梯度近似值；

若所述目标函数的梯度近似值的绝对值小于或等于预设精度值，则所述预设线性滑动理论参数为获取的线性滑动理论参数；

若所述目标函数的梯度近似值的绝对值大于所述预设精度值，则根据预设规则对所述预设线性滑动理论参数进行更新，重新获取所述目标函数的梯度近似值，直至所述目标函数的梯度近似值的绝对值小于或等于所述精度值为止，当前的所述预设线性滑动理论参数为获取的线性滑动理论参数；

所述获取所述目标函数的梯度近似值的步骤包括：

根据所述预设线性滑动理论参数、所述预设增量值、所述测井资料、所述地质资料、所述预设模型，获得第一理论声波全波列波形及第二理论声波全波列波形；

对所述第一理论声波全波列波形及所述第二理论声波全波列波形分别进行滤波，分离出对应的第一理论斯通利波波形及第二理论斯通利波波形；

根据所述实际斯通利波波形、所述第一理论斯通利波波形、所述第二理论斯通利波波形及预设增量值，获取所述目标函数的梯度近似值。

3.根据权利要求1所述的方法，其特征在于，所述线性滑动理论参数包括正向屈服度系数及切向屈服度系数。

4.根据权利要求1至3任一项所述的方法，其特征在于，所述预设模型包括含孔隙的流体介质理论模型及线性滑动理论裂缝模型。

5.根据权利要求1所述的方法，其特征在于，所述测井资料包括由多个实际声波全波列波形构成的实际阵列声波全波列波形；

对所述实际阵列声波全波列波形中的多个实际声波全波列波形分别进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形构成实际阵列斯通利波波形；

建立所述实际阵列斯通利波波形与理论阵列斯通利波波形误差的目标函数；

根据预设最优化方法，利用所述实际阵列斯通利波波形及理论阵列斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论阵列斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论阵列声波全波列波形进行滤波后得到的阵列波形。

6.根据权利要求5所述的方法，其特征在于，所述对所述实际阵列声波全波列波形中的多个实际声波全波列波形分别进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形构成实际阵列斯通利波波形，具体包括：

根据所述测井资料及采用慢度-时间相关法，获取斯通利波时差；

根据所述测井资料及所述斯通利波时差，对所述实际阵列声波全波列波形进行频谱分析，确定实际斯通利波波形的频率范围；

对所述实际阵列声波全波列波形中的多个实际声波全波列波形进行滤波，分离出与所述多个实际声波全波列波形一一对应的多个实际斯通利波波形，所述多个实际斯通利波波形与所述频率范围对应，所述多个实际斯通利波波形构成实际阵列斯通利波波形。

7.根据权利要求1所述的方法，其特征在于，所述对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形，具体包括：

根据所述测井资料，采用阈值法或人机交互的方法进行斯通利波时差提取，获取斯通利波时差；

根据所述测井资料及所述斯通利波时差，对所述实际声波全波列波形进行频谱分析，确定实际斯通利波波形的频率范围；

对所述实际声波全波列波形进行滤波，分离出所述频率范围对应的实际斯通利波波形。

8.根据权利要求1所述的方法，其特征在于，所述对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形之前，所述方法还包括：

对所述测井资料进行预处理，所述预处理包括井眼环境校正、数据滤波、坏道剔除或数据缺失值填补。

9.根据权利要求1所述的方法，其特征在于，所述测井资料还包括井径测井数据、自然伽马测井数据、密度测井数据、微电阻率成像测井数据、所述实际声波全波列波形增益数据及所述实际声波全波列波形记录开始时间数据，所述地质资料包括目标井地层层位信息及岩性信息。

10.一种基于测井资料的线性滑动理论参数反演装置，其特征在于，所述装置包括：

接收模块，用于获取地质资料及测井资料，所述测井资料包括实际声波全波列波形；

第一处理模块，用于对所述实际声波全波列波形进行滤波，分离出实际斯通利波波形；

第二处理模块，用于建立实际斯通利波波形与理论斯通利波波形误差的目标函数；根据预设最优化方法，利用所述实际斯通利波波形及理论斯通利波波形，获取当目标函数满足预设条件时的线性滑动理论参数，所述理论斯通利波波形是对根据预设线性滑动理论参数、预设增量值、所述测井资料、所述地质资料及预设模型获得的理论声波全波列波形进行滤波后得到的波形。