CN106324689A - 一种水平井地层环境下电阻率各向异性识别方法 - Google Patents

Abstract

本发明公开了一种水平井地层环境下电阻率各向异性识别方法，首先根据测井获得方位自然GR和方位电阻率测井并提取地层倾角和地层方位角；从测井获得水平井中多个源距和不同频率下的仪器方位电阻率测井响应值，通过环境校正和自动分层、假层剔除，得到标准化的测井数据和地层分层界面；然后根据倾角提取结果、分层界面、视电阻率测量值和邻井数据，构建水平井地层模型；最后反演计算地层模型，输出各向异性系数和地层电阻率反演结果。该方法应用随钻方位电阻率测井数据提取水平井环境下地层真实电阻率信息和各向异性系数信息，完善了现有水平井测井评价资料，在水平井测井资料综合解释评价中，通过地层真实电阻率和各向异性系数反演结果可精确计算储层含油饱和度等，为储层评价提供可靠参数。

Description

一种水平井地层环境下电阻率各向异性识别方法

技术领域

本发明涉及油田开发技术领域，尤其涉及一种水平井地层环境下电阻率各向异性识别方法。

背景技术

随着复杂油气田勘探开发的不断深入，大斜度井、水平井等复杂工艺井的广泛应用，随钻测井技术研究与随钻测井仪器研发得到足够的重视，并得到了快速的发展。传统的随钻电磁波测井仪器发射线圈与接收线圈共轴，测量得到的地层信号为地层信息的平均值，不具备方位特性，不能够准确获得地层电阻率各向异性等信息。方位随钻电磁波测井仪器与传统仪器有很大的不同，方位电磁波测量仪器均采用轴向倾斜或横向线圈混合，能够更好地提供地层方位信息，指示地层的各向异性，并识别地层边界。三大石油测井服务公司相继推出了具有方位探测能力的随钻方位电磁波测井仪器，2005年斯伦贝谢公司推出了PeriScope方位电阻率测量仪，2006年贝克休斯公司推出了随钻方位电磁波测井仪APR，2007年哈里伯顿公司推出了方位深电阻率测量仪ADR。方位电磁波测量仪器均采用轴向，倾斜或横向线圈混合，能够更好地提供有关地层方位的信息，指示地层的各向异性，识别地层边界。

电测井响应的地层电阻率各向异性反演属于非线性问题，单一测井信息的片面性与反演的多解性，增加了测井反演、测井解释与评价的难度。联合反演为解决这一问题提供了有效途径。联合反演是指利用不同物理机制的两种或两种以上测井数据进行地质模型参数反演。联合反演在本质上是通过增加特定探测目标的有效信息量(增加约束)，来达到更准确反映地质目标体的目的。

物性同源是联合反演增加该源有效信息量的基本条件。同一口井不同测井系列的随钻电阻率测井与电缆电阻率测井，针对相同原状地层具有相同的目标物性测量项目，使其联合反演成为可能。随钻测井与电缆测井由于时间的推移具有不同的泥浆侵入深度(侵入半径)及侵入带电阻率，且随钻测井仪器与电缆测井仪器两者间探测特性及各自不同测井曲线间的探测特性存在差异，因此随钻电缆测井联合反演较单一反演具有相当丰富的测井信息。

随钻电缆联合反演以阻尼最小二乘法为基础，将随钻与电缆测井响应、随钻与电缆测井时刻地层模型参数有机统一起来进行反演。但是，传统电阻率联合反演主要是针对旋转对称性地层进行一维和(或)二维反演，既不适用非对称地层结构的水平井和大斜度井，也不能提供层边界距离和各向异性等信息。因此，需要一种能提供层边界距离和各向异性等信息的基于三维空间地层模型的地层电阻率各向异性识别方法。

发明内容

(一)要解决的技术问题是提供一种能提供各向异性等信息的基于三维空间水平井地层模型的随钻方位电阻率测井各向异性反演和识别方法。该方法利用水平井中方位电阻率测井数据对各向异性的很好指示作用，并联合视电阻率曲线、自然GR曲线和中子密度孔隙度曲线，在三维空间水平井地层模型上实现了地层各向异性反演和识别方法。

(二)技术方案

本发明的目的是通过以下技术方案来实现的：

本发明提供的水平井地层环境下电阻率各向异性识别方法，包括以下步骤：

步骤1，根据水平井测井获得方位自然GR和方位电阻率测井并提取地层倾角和地层方位角；

步骤2，水平井测井获得多个源距和不同频率下的仪器方位电阻率测井响应值，测井获得中子密度和孔隙度曲线，并利用地层倾角和地层方位角，对方位电阻率测井和方位电阻率测井数据进行方位校正、环境校正得到标准化的测井数据；

步骤3，将标准化的水平井测井数据进行预处理，并进行相关对比和活度法进行地层自动分层得到地层界面；

步骤4，筛选地层界面并剔除自动分层产生的假层得到真实层界面；

步骤5，根据倾角和真实层界面构建地层模型，并初始化地层模型；

步骤6，利用有限元方法计算水平井地层初始模型的随钻方位仪器正演响应，对比水平井实测数据与地层初始模型的正演响应，建立目标函数，用最小二乘方法进行多参数迭代反演求解目标函数最优解的模型改变量得到反演结果；

步骤7，格式化输出反演结果。

进一步，在步骤1中是通过利用水平井测井曲线高程差和非线性拟合方法来提取地层倾角和方位角，具体步骤如下：

步骤11：利用ABG测井仪器测量的若干方位自然GR曲线进行倾角提取，对每条GR曲线进行相关对比分析，得到对应地层的高程差，获取地层倾斜层面上的测量点；并通过以下方程来表示井壁与倾斜层相交的展开图：

y＝Asin(ωx-β)+y₀

＝y₀+Acosβsin(ωx)+(-Asinβ)cos(ωx) (1)；

式中，y表示井壁与倾斜层相交的函数值；A表示函数值变化振幅；y0表示函数值y的均值；β表示函数值y所满足正弦函数的初始相位；ω表示函数值y所满足正弦函数的周期；

步骤12：通过非线性拟合求解方程(1)中的待定系数，利用最小二乘法得到的矩阵方程：

其中：

a₀、a₁、a₂表示拟合的多项式系数；表示拟合函数时自变量取值；

步骤13：求解步骤12式，得正弦曲线参数表达式和和倾角倾向表达式：

$\left\{\begin{array}{c}{y}_0=a_0\\ \mathrm{\β}=\arctan (-a_2/a_1)\\ A=a_1/cos\mathrm{\β}\end{array}\right.---\left(3\right);$

倾角倾向表达式：

$\left\{\begin{array}{c}Dip=\arctan \left(\frac{2A}{D_{el}}\right)\\ Dir=x|y=y_{min}\end{array}\right.---\left(4\right);$

式中，Dip表示倾角；Dir表示倾向；D_el表示电直径即探测深度；y_min表示表示井壁与倾斜层相交的一周上最小值点。

进一步，在步骤2中，利用倾角和方位角对方位信息中的方位电阻率测井和方位电阻率测井数据进行方位校正，具体步骤如下：

步骤21：通过以下方程表示所述方位信息：

式中，y(x)表示方位曲线函数值；B表示方位曲线y的平均值；A表示位曲线满足正弦函数的振幅；

步骤22：按照以下解析解来求取所述方程(5)获得A、B和的值：

步骤23：将A、B值带人所述方程(5)形成已知系数方程；

步骤24：将已知系数方程所对应的函数值作为方位校正值。

进一步，在步骤3中，利用相关对比和活度法进行地层自动分层，具体步骤如下：

步骤31：对测井数据进行光滑滤波处理；

步骤32：对测井数据进行归一化处理；

步骤33：计算水平井测井数据的相关性系数并进行相关对比；计算水平井测井数据的曲线活度并进行活度分层；所述相关对比是指不同探测深度测井曲线归一化后数据在固定窗长下的相关性对比；所述活度分层是按以下方式进行的：首先计算曲线活度，判断曲线活度是否大于预设的活度阈值，如果是，则将大于活度阈值的极值点的位置定义为分层界面；

所述相关性系数按以下公式来计算：

$Z_i=\frac{{\mathrm{\Σ}}_{i=1}^n(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{{(x_i-\stackrel{\‾}{x})}^2+{(y_i-\stackrel{\‾}{y})}^2}}$

式中，i表示采样点位置，n表示窗长，表示窗长内对比条曲线平均值，表示被对比条曲线平均值，x_i表示对比曲线第i个采样点值，y_i表示被对比曲线第i个采样点曲线值，Z_i表示相关系数。

进一步，在步骤31中对标准化的测井数据进行预处理采用最小二乘滑动平均法或加权滑动平均法进行平滑滤波，所述平滑滤波采用以下的方式：

线性函数平滑：

二次函数平滑：

钟形函数平滑法：

汉明函数平滑法：

式中；T_i表示第i个采样点值，m表示滤波窗长，表示滤波后第i各采样点取值。

进一步，在步骤32中对对测井数据进行归一化处理采用极值归一化方法、孔隙度归一化方法或密度归一化方法进行，归一化方法具体如下：

对电阻率采样极值归一化方法：

对中子孔隙度采样线性归一化方法：

对中子密度采样线性归一化方法：

式中；x_ij表示第j条曲线的第i个采样点值，x_minj表示第j条曲线的最小值，x_maxj表示第j条曲线的最大值，Y_ij表示归一化后第j条曲线的第i个值；TNPL_j表示孔隙度曲线第j个值，用归一化后的值代替原始值；ALCD_j表示中子密度第j个值，用归一化后的值代替原始值。

进一步，在步骤33中，对测井数据进行相关对比和活度分层，具体步骤如下：

步骤331：所述活度分层中的活度按照以下公式进行定义：

$E\left(d\right)={\mathrm{\Σ}}_{i=d-n}^{i=d+n}{\[x\left(i\right)-\stackrel{\‾}{x}\left(d\right)\]}^2---\left(12\right)$

$\stackrel{\‾}{x}\left(d\right)=\frac{1}{2n+1}{\mathrm{\Σ}}_{i=d-n}^{i=d+n}x\left(i\right)---\left(13\right)$

式中，E(d)表示d的活度函数值，x(t)表示测井曲线测量值，表示测井曲线在区间[d-n，d+n]内的平均值；d表示采样点位置；n表示采样窗长；i表示在一个采样窗长中第几个采样点；

步骤332：根据活度函数值与预设活度阈值进行比较来确定层界面，所述活度阈值为所有活度极值从大到小排列的前5％处作为活度阈值，大于活度阈值的极值点作为层界面，小于阈值的点不做层界面。

进一步，在步骤4中的剔除自动分层产生的假层，具体步骤如下：

步骤41：根据层界面位置，用两个相邻层界面相减的绝对值，即为相邻层界面间垂直距离；

步骤42：判断两相邻界面垂直距离是否小于预设距离值，如果是，则该层界面为无效界面，并视为假层；

步骤43：如果否，则返回步骤41重复直至结束；

其中，在计算出相邻层界面垂直距离后，根据层界面垂直距离，把层界面距离小于预设定值的层厚看作是假层，并剔除假层。

进一步，在步骤6中利用有限元方法计算建立的地层初始模型正演响应，具体步骤如下：

步骤61)求解给定边界条件下麦克斯韦Maxwell方程的问题，将Maxwell方程转化为波动方程：

$\▿\×(\▿\×\frac{\stackrel{\→}{E}}{\mathrm{\μ}})-{\mathrm{\ω}}^2\mathrm{\ϵ}\stackrel{\→}{E}=-j\mathrm{\ω}{\stackrel{\→}{J}}_s---\left(14\right);$

其中，E表示电场强度；μ表示振幅；ω表示角频率；ε表示介电常数；

步骤62)结合边界条件波动方程归结为场能量泛函：

$F\left(\stackrel{\→}{E}\right)=\frac{1}{2}\∫\∫{\∫}_V\[\frac{1}{{\mathrm{\μ}}_r}\left(\▿\×\stackrel{\→}{E}\right)\·\left(\▿\×\stackrel{\→}{E}\right)-{\mathrm{\ω}}^2\mathrm{\ϵ}\stackrel{\→}{E}\×\stackrel{\→}{E}\]dV+{\mathrm{j\ω\μ}}_0\∫\∫{\∫}_V\stackrel{\→}{E}\×\stackrel{\→}{J}dV---\left(15\right);$

其中，表示目标函数；μ₀表示初始振幅；V表示求解区域；ω表示角频率；ε表示介电常数；J表示电流密度；

步骤63)应用有限单元剖分场域，并选取相应的插值基函数，对能量泛函进行空间离散，得到如下离散化泛函形式：

$F=\frac{1}{2}{\mathrm{\Σ}}_{e=1}^M({\left\{E^e\right\}}^T\[A^e\]\{E^e\}-{\mathrm{\ω}}^2{\mathrm{\ϵ\μ}}_0{\left\{E^e\right\}}^T\[B^e\]\{E^e\}+{\mathrm{j\ω\μ}}_0{\left\{E^e\right\}}^T\{C^e\left\}\right)---\left(16\right)$

其中，A^e表示函数离散形成系数矩阵；B^e表示等式右端项离散矩阵；C^e表示函数值离散矩阵；E^e表示自变量矩阵；M表示自变量个数；T表示矩阵的转置；

步骤64)求解线性方程组得到所需的参数。

进一步，在步骤6中利用最小二乘方法进行多参数迭代反演，具体步骤如下：

步骤621)通过对比模型正演响应和实测数据，建立目标函数，通过最小二乘法求解实际测量值与模拟值的残差：

$\underset{x\∈R^n}{min}f\left(x\right)=\frac{1}{2}r{\left(x\right)}^{T}r\left(x\right)=\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{\[r_i\left(x\right)\]}^2,m\≥n---\left(17\right)$

其中，r:Rⁿ→R^m是x的非线性函数；R表示自变量向量；m表示未知量个数；n表示正演响应曲线个数；

步骤622)通过梯度方法求解目标函数雅可比矩阵，形成雅可比线性方程组，求解方程组计算目标函数最速下降方向；

按照以下公式来计算目标函数的梯度：

$g\left(x\right)={\mathrm{\Σ}}_{i=1}^m{r}_i\left(x\right)\▿r_i\left(x\right)=J{\left(x\right)}^{T}r\left(x\right)---\left(19\right)$

其中，J(x)是r(x)的雅可比矩阵：

求解目标函数最速下降方向；

步骤623)利用黄金分割方法，在最速下降方向上确定步长；

步骤624)根据求得模型改变方向和改变步长，确定模型改变量，改变模型，完成一次迭代反演；

步骤625)设置迭代终止条件，循环调用迭代反演，直至满足迭代终止条件，输出结果。

(三)有益效果

与现有技术和产品相比，本发明有如下优点：

本发明采用方位随钻电阻率测井识别三维水平井地层中电阻率各向异性，该方法应用方位电阻率测井数据反演地层各向异性能够获得井眼-地层相对倾角、泥浆侵入深度、侵入带电阻率、水平电阻率以及各向异性系数等，完善了现有测井评价资料，可得到全面地层信息。在后期解释评价中，通过每层地层电阻率各向异性反演结果可精确计算储层含油饱和度、可动油饱和度等，为储层评价提供可靠参数。

附图说明

图1为本发明的水平井地层环境下电阻率各向异性识别方法具体实施例的流程图。

图2为本发明的水平井地层建模示意图。

图3为本发明的水平井下视电阻率随地层各向异性系数变化示意图。

图4为本发明的水平井环境下视电阻率的相位差电阻率与幅度比电阻率差值随各向异性系数变化示意图。

图5为本发明的水平井中各向异性地层中视电阻率随倾角变化示意图。

图6为本发明的实施例中对某段实测水平井资料的地层各向异性反演结果图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实施方式对本发明作进一步的详细描述。

实施例1

本实施例提供的水平井地层环境下电阻率各向异性识别方法，包括以下步骤：

步骤1，根据水平井测井获取方位自然GR和方位电阻率测井响应值，利用高程差和曲线拟合提取测井仪器-地层相对倾斜角度和地层方位角；

步骤2，根据步骤1获得的倾角和方位角，对方位电阻率测井和方位电阻率测井数据进行方位校正；

步骤3，对进行方位校正后的标准化数据，通过相关对比方法求取曲线活度，用活度法进行地层界面划分，地层界面划分时考虑自然GR曲线、方位电阻率测井曲线、方位电阻率曲线、中子密度曲线和孔隙度曲线；

本实施例利用相关对比和活度法进行地层自动分层，包括以下步骤：

测井数据光滑滤波和无效数据剔除；

数据归一化处理；

相关对比和活度分层。

本实施例在进行相关对比和自动分层前，需对测量的数据进行预处理，预处理的内容包括曲线光滑滤波，去除无效数据并对缺失数据进行合理范围内最大限度的补充。常用的平滑滤波方法有最小二乘滑动平均法、加权滑动平均法等，其中上述两种方法又包含了多种具体的实现类型，最小二乘滑动平均法、加权滑动平均法：

$\stackrel{\‾}{T_i}=\frac{1}{2m+1}{\mathrm{\Σ}}_{k=-m}^{+m}{T}_{i+k}---\left(5\right)$

$\stackrel{\‾}{T_i}=\frac{1}{35}(-3(T_{i-2}+T_{i+2})+12(T_{i-1}+T_{i+1})+17T_i---\left(6\right)$

$\stackrel{\‾}{T_i}=0.11(T_{i-2}+T_{i+2})+0.24(T_{i-1}+T_{i+1})+0.3T_i---\left(7\right)$

$\stackrel{\‾}{T_i}=0.04(T_{i-2}+T_{i+2})+0.24(T_{i-1}+T_{i+1})+0.44T_i---\left(8\right)$

本实施例对滤波和无效数据处理后的数据进行归一化处理，归一化方法为：

极值归一化方法：

$Y_{ij}=\frac{{x_{ij}}^{0.26}-{x_{\min j}}^{0.26}}{{x_{\max j}}^{0.26}-{x_{\min j}}^{0.26}}---\left(9\right)$

孔隙度归一化方法：

${\mathrm{TNPL}}_j=\frac{{\mathrm{TNPL}}_j+0.15}{0.45-0.15}---\left(10\right)$

密度归一化方法：

${\mathrm{ALCD}}_j=\frac{{\mathrm{ALCD}}_j+1.95}{2.95-1.95}---\left(11\right)$

本实施例对数据进行相关对比和活度分层，包括以下步骤：

1)对归一化后的数据进行相关对比和活度分层，活度的定义为：

$E\left(d\right)={\mathrm{\Σ}}_{i=d-n}^{i=d+n}{\[x\left(i\right)-\stackrel{\‾}{x}\left(d\right)\]}^2---\left(12\right)$

$\stackrel{\‾}{x}\left(d\right)=\frac{1}{2n+1}{\mathrm{\Σ}}_{i=d-n}^{i=d+n}x\left(i\right)---\left(13\right)$

式中，E(d)表示d的活度函数值，x(t)表示测井曲线测量值，表示测井曲线在区间[d-n，d+n]内的平均值。

2)跟据3)中求得的活度值，设置阈值，根据阈值确定层界面。

步骤4，对步骤3划分的层界面进行筛选和甄别，剔除假层，获得真实层界面；

本实施例结合水平井层界面垂直深度、地层倾角信息，计算相邻层界面见垂直距离，通过距离的大小进行有效界面判断，两相邻界面距离过小时视其中一个位无效界面视为假层，根据活度值大小，对界面进行筛选和甄别，剔除假层，获得真实层界面。

步骤5，结合倾角、界面、视电阻率、井眼轨迹垂深-斜深-水平位移变化，建立初步三维水平井多层地质模型；

本实施例结合倾角、界面、视电阻率、井眼轨迹垂深-斜深-水平位移变化，建立初步三维多层地质模型。根据视电阻率分离特征，初始化地层背景电阻率、泥浆侵入、各向异性等参数；跟据井径、泥浆等已知信息，初始化地层模型中井眼参数；根据倾角、层界面、垂深初始化地层模型中井眼轨迹与层界面位置参数；

步骤6，利用有限元方法计算三维地层模型正演响应，通过模型正演响应和实测数据对比，建立目标函数，用最小二乘方法求解目标函数最优解的模型改变量，进行迭代反演；

本实施例利用有限元方法计算三维水平井地层模型正演响应；利用最小二乘方法进行迭代反演

其中，利用有限元方法计算建立的水平井地层模型响应的步骤，包含以下几个步骤：

1)有限元方法计算建立的地层模型响应实质归结为求解给定边界条件下麦克斯韦(Maxwell)方程的问题，Maxwell方程转化为波动方程：

$\▿\×(\▿\×\frac{\stackrel{\→}{E}}{\mathrm{\μ}})-{\mathrm{\ω}}^2\mathrm{\ϵ}\stackrel{\→}{E}=-j\mathrm{\ω}{\stackrel{\→}{J}}_s---\left(14\right)$

2)结合边界条件波动方程归结为场能量泛函：

$F\left(\stackrel{\→}{E}\right)=\frac{1}{2}\∫\∫{\∫}_V\[\frac{1}{{\mathrm{\μ}}_r}\left(\▿\×\stackrel{\→}{E}\right)\·\left(\▿\×\stackrel{\→}{E}\right)-{\mathrm{\ω}}^2\mathrm{\ϵ}\stackrel{\→}{E}\×\stackrel{\→}{E}\]dV+{\mathrm{j\ω}}_0\∫\∫{\∫}_V\stackrel{\→}{E}\×\stackrel{\→}{J}dV---\left(15\right)$

3)应用有限单元剖分场域，并选取相应的插值基函数，对能量泛函进行空间离散，可以得到如下离散化泛函形式：

$F=\frac{1}{2}{\mathrm{\Σ}}_{e=1}^M({\left\{E^e\right\}}^T\[A^e\]\{E^e\}-{\mathrm{\ω}}^2{\mathrm{\ϵ\μ}}_0{\left\{E^e\right\}}^T\[B^e\]\{E^e\}+{\mathrm{j\ω}}_0{\left\{E^e\right\}}^T\{C^e\left\}\right)---\left(16\right)$

4)求解这个线性方程组得到所需的参数。

其中，利用最小二乘方法进行多参数反演的步骤，包含以下几个步骤：

1)通过对比模型正演响应和水平井实测数据，建立目标函数，通过最小二乘法求解实际测量值与模拟值的残差：

$\underset{x\∈R^n}{min}f\left(x\right)=\frac{1}{2}r{\left(x\right)}^{T}r\left(x\right)=\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{\[r_i\left(x\right)\]}^2,m\≥n---\left(17\right)$

其中r:Rⁿ→R^m是x的非线性函数,把非线性最小二乘问题看作为无约束极小化的特殊情形

2)通过梯度方法求解目标函数雅可比矩阵，形成雅可比线性方程组，求解方程组计算目标函数最速下降方向。设J(x)是r(x)的雅可比矩阵：

则目标函数的梯度为：

$g\left(x\right)={\mathrm{\Σ}}_{i=1}^m{r}_i\left(x\right)\▿r_i\left(x\right)=J{\left(x\right)}^{T}r\left(x\right)---\left(19\right)$

令方程组(19)等0，即可求解目标函数最速下降方向；

3)利用黄金分割方法，在最速下降方向上确定步长；

4)根据求得水平井模型改变方向和改变步长，确定模型改变量，改变模型，完成一次迭代反演

5)设置迭代终止条件，循环调用迭代反演，直至满足迭代终止条件，输出结果。

步骤7，格式化输出反演结果；

本实施例根据反演结果，对测量点对应位置按测井曲线间隔对反演结果重采样，同时对逐层反演的方波形结果进行非线性插值，获得光滑的结果曲线，便于绘制图件。

实施例2

本实施例提供一种可以进在水平井和大斜度井等复杂井眼环境下利用方位电阻率测井进行地层各向异性识别的方法；通过如下步骤来实现：

步骤1利用测井曲线高程差和非线性拟合方法，提取地层倾角和方位角，包含以下几个步骤：

1)利用ABG测井仪器测量的四条方位自然GR曲线进行倾角提取，对四条GR曲线进行相关对比分析，可以得到对应层的高程差，即为地层倾斜层面上的六个点。由于井壁与倾斜层相交的展开图在图像上表现为单周期的正弦函数，满足方程：

y＝Asin(ωx-β)+y₀ (1)

2)非线性拟合，求解方程(1)中的待定系数，利用最小二乘法得到的矩阵方程：

其中：a₀＝y₀，a₁＝Acosβ，a₂＝-Asinβ，φ₀(x)＝1，φ₁(x)＝sin(ωx)，φ₂(x)＝cos(ωx)

3)求解(2)式，可得正弦曲线参数表达式：

$\left\{\begin{array}{c}{y}_0=a_0\\ \mathrm{\β}=\arctan (-a_2/a_1)\\ A=a_1/cos\mathrm{\β}\end{array}\right.---\left(3\right)$

和倾角倾向表达式：

$\left\{\begin{array}{c}Dip=\arctan \left(\frac{2A}{D_{el}}\right)\\ Dir=x|y=y_{min}\end{array}\right.---\left(4\right)$

或者利用水平井测井曲线高程差和非线性拟合方法，提取地层倾角和方位角，包含以下几个步骤：

1)利用其它水平井测井仪器测量的三条或四条带方位信息曲线(方位自然GR、方位电阻率测井、方位电阻率、方位密度、方位孔隙度等)进行倾角提取，对方位曲线进行相关对比分析，可以得到对应层的高程差，即为地层倾斜层面上的三个和三个以上点。由于井壁与倾斜层相交的展开图在图像上表现为单周期的正弦函数，满足方程：

y＝Asin(ωx-β)+y₀或y＝Acos(ωx-β)+y₀ (5)

2)用解析解或者非线性拟合方法，求解方程(1)中的待定系数，利用最小二乘法得到的矩阵方程：

3)求解(2)式，可得正弦曲线参数表达式：

和倾角倾向表达式：

步骤2中，利用骤1获得的倾角和方位角，对方位电阻率测井和方位电阻率测井数据进行方位校正，方位信息满足正弦(或余弦)规律，且周期为2π，可设曲线方程为：

已知在x＝0，1/2π，π，3/2π的值，通过解析解求A、B和解析解如下：

求得后，A、B值，取函数y(x)＝Asin(x)+B或y(x)＝Acos(x)+B在x＝0，1/2π，π，3/2π处的值作为方位校正后值。

步骤3中，利用相关对比和活度法进行地层自动分层，包括的测井数据光滑滤波和无效数据剔除，采用大于3个点的曲线拟合方法进行滤波，或采用插值法包括线性插值、非线性插值、面积插值；数据归一化处理采用指数归一化方法、线性归一化方法；相关对比和活度分层或者拐点方法分层。

步骤4，对薄层数目进行确定，剔除自动分层产生的假层；

步骤5，根据步骤1提取的倾角，步骤3获得的分层数据，结合视电阻率构建地层初始模型；

步骤6中，计算构建的初始地层模型的方位电阻率测井正演响应，与实测曲线对比，用最小二乘方法进行多参数迭代反演；利用有限元方法计算建立的地层模型响应。

利用有限元方法进行方位随钻电磁波仪器的数值仿真，研究了仪器的方位电阻率测井响应同水平井地层各向异性的关系，研究结果表明方位电阻率测井对各向异性有很好的指示作用，在各向异性地层，幅度比电阻率小于相位差电阻率，且随各向异性系数增大，二者差异增大；随着倾角变化，方位随钻测井仪器测井响应分离，倾角越大曲线分离越大。此时，综合使用地层电阻率和方位电阻率测井进行地层各向异性识别是有效的。因此，本发明提出利用仪器测量得到的方位电阻率测井和电阻率曲线进行地层各向异性识别的方法，并给出侵入剖面图，在地质导向和钻后评价中均具有重要意义。本发明在不增加现有测井技术成本的前提下完成对水平井地层各向异性的提取，更为行之有效的反应地层真实信息，完善了现有水平井测井评价资料；本发明使用同时三维反演得到的水平井地层环境下各向异性地层电阻率和侵入剖面，更加贴近真实地层状况，得到全面地层参数信息，进行更为真实有效的地层评价和储量计算。

实施例3

如图所示，图1为本发明的水平井地层环境下电阻率各向异性识别方法具体实施例的流程图；图2为本发明的水平井地层建模示意图，图中TVD表示层边界位置，RM表示泥浆电阻率，RI表示侵入深度，RXO表示侵入带电阻率，RS表示围岩电阻率，Rh-Rv表示目的水平电阻率-垂直电阻率，DH表示井眼直径；图3为本发明的水平井下视电阻率随地层各向异性系数变化；图4为本发明的视电阻率的相位差电阻率与幅度比电阻率差值随各向异性系数变化；图5为本发明的各向异性地层中视电阻率随倾角变化；图6为本发明的实施例中对某段实测井资料的地层各向异性反演结果图，其中RI为侵入深度，TVD为井眼垂深，DTB为层边界距离，GR为自然GR，ARM48P为48in中频相位差电阻率，ARH48P为48in高频相位差电阻率，ARM32P为32in中频相位差电阻率，ARM16P为16in中频相位差电阻率，RXO为侵入带电阻率，RH为水平电阻率，RV为垂直电阻率，TNPL为中子孔隙度，ALCDLC为中子密度。

如图1所示，图1为本发明的水平井地层环境下电阻率各向异性识别方法的流程图。

步骤101，水平井地层电阻率各向异性识别是基于三维空间进行的，地层倾角和和方位角等非对称性，对方位电阻率测井和电阻率曲线有较大影响。因此首要的是结合三维实际情况，联合其它测井资料进行倾角提取和地层方位角提取；

步骤102，在进行水平井地层各向异性识别时，计算出的层边界距离是采样点到地层界面间的距离，点到面的距离，在地层建模时为简化模型，层界面倾角默认为0度，因此需对方位测井曲线进行倾角校正；

步骤103，在对水平井中方位测井数据进行方位校正后得到标准化的测井数据，0度方位对应高边或底边，180度方位与之相反。应用相关对比和活度方法进行分层，分层前因实测数据可能有异常值，需进行光滑滤波处理和对异常数据重采样，消除异常值对分层结果的影响；

步骤104，结合层界面垂直深度、地层倾角信息，计算相邻层界面见垂直距离，通过距离的大小进行有效界面判断，两相邻界面距离过小时其中一个视为无效界面视为假层，根据活度值大小，对界面进行筛选和甄别，剔除假层，获得真实层界面。

步骤105，结合倾角、界面、视电阻率、井眼轨迹垂深-斜深-水平位移变化，建立初步三维多层地质模型。根据视电阻率分离特征，初始化地层背景电阻率、泥浆侵入、各向异性等参数；跟据井径、泥浆等已知信息，初始化地层模型中井眼参数；根据倾角、层界面、垂深初始化地层模型中井眼轨迹与层界面位置参数；

步骤106，水平井地层电阻率各向异性识别是基于正演仿真和迭代反演完成的。因此首要的是建立正演仿真的三维水平井地层模型，其建立步骤如下：

联合其它测井资料，确定地层倾角和方位角；

结合电阻率曲线、自然GR曲线和中子密度孔隙度曲线对地层分层，确定层界面个数并对每层背景电阻率赋初值；

步骤107，根据反演结果，对测量点对应位置按测井曲线间隔对反演结果重采样，同时对逐层反演的方波形结果进行非线性插值，获得光滑的结果曲线，便于绘制图件。

本实施例在计算时，充分考虑到水平井、大斜度井三维空间复杂环境中井眼、泥浆、倾角、侵入、各向异性、层厚及各个层中背景电阻率等多因素的影响，并且联合其它测井数据进行地层建模和界面划分、倾角提取，计算过程中充分考虑到各个影响因素耦合作用，更贴近于实际测井环境情况，得到的结果更为贴近真实值。同时一方面解决了传统反演方法不能进行非对称地层和大倾角地层的反演；另一方面，可提供地层各向异性的识别信息，获得水平电阻率和垂直电阻率，在实际生产应用中具有非常重要的应用价值

以上实施例仅为本发明的一种实施方式，其描述较为具体和详细，但不能因此而理解为对本发明专利范围的限制。其具体结构和尺寸可根据实际需要进行相应的调整。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (10)

1.一种水平井地层环境下电阻率各向异性识别方法，其特征在于，包括以下步骤：

步骤1，根据水平井测井获得方位自然GR和方位电阻率测井并提取地层倾角和地层方位角；

步骤2，水平井测井获得多个源距和不同频率下的仪器方位电阻率测井响应值，水平井测井获得中子密度和孔隙度曲线，并利用地层倾角和地层方位角，对方位电阻率测井和方位电阻率测井数据进行方位校正、环境校正得到标准化的测井数据；

步骤3，将标准化的测井数据进行预处理，并进行相关对比和活度法进行地层自动分层得到地层界面；

步骤4，筛选地层界面并剔除自动分层产生的假层得到真实层界面；

步骤5，根据倾角和真实层界面构建地层模型，并初始化地层模型；

步骤6，利用有限元方法计算三维水平井地层初始模型的随钻方位仪器正演响应，对比实测数据与地层初始模型的正演响应，建立目标函数，用最小二乘方法进行多参数迭代反演求解目标函数最优解的模型改变量得到反演结果；

步骤7，格式化输出反演结果。

2.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤1中是通过利用测井曲线高程差和非线性拟合方法来提取地层倾角和方位角，具体步骤如下：

步骤11：利用ABG测井仪器测量的若干方位自然GR曲线进行倾角提取，对每条GR曲线进行相关对比分析，得到对应地层的高程差，获取地层倾斜层面上的测量点；并通过以下方程来表示井壁与倾斜层相交的展开图：

y＝Asin(ωx-β)+y₀

＝y₀+Acosβsin(ωx)+(-Asinβ)cos(ωx) (1)；

式中，y表示井壁与倾斜层相交的函数值；A表示函数值变化振幅；y₀表示函数值y的均值；β表示函数值y所满足正弦函数的初始相位；ω表示函数值y所满足正弦函数的周期；

步骤12：通过非线性拟合求解方程(1)中的待定系数，利用最小二乘法得到的矩阵方程：

其中：

a₀＝y₀，a₁＝Acosβ，a₂＝-Asinβ，

a₀、a₁、a₂表示拟合的多项式系数；表示拟合函数时自变量取值；

步骤13：求解步骤12式，得正弦曲线参数表达式和和倾角倾向表达式：

$\left\{\begin{array}{c}{y}_0=a_0\\ \mathrm{\β}=\arctan (-a_2/a_1)\\ A=a_1/cos\mathrm{\β}\end{array}\right.---\left(3\right);$

倾角倾向表达式：

$\left\{\begin{array}{c}Dip=\arctan \left(\frac{2A}{D_{el}}\right)\\ Dir=x|y=y_{min}\end{array}\right.---\left(4\right);$

式中，Dip表示倾角；Dir表示倾向；D_el表示电直径即探测深度；y_min表示表示井壁与倾斜层相交的一周上最小值点。

3.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤2中，利用倾角和方位角对方位信息中的方位电阻率测井和方位电阻率测井数据进行方位校正，具体步骤如下：

步骤21：通过以下方程表示所述方位信息：

式中，y(x)表示方位曲线函数值；B表示方位曲线y的平均值；A表示位曲线满足正弦函数的振幅；

步骤22：按照以下解析解来求取所述方程(5)获得A、B和的值：

步骤23：将A、B值带人所述方程(5)形成已知系数方程；

步骤24：将已知系数方程所对应的函数值作为方位校正值。

4.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤3中，利用相关对比和活度法进行地层自动分层，具体步骤如下：

步骤31：对水平井测井数据进行光滑滤波处理；

步骤32：对水平井测井数据进行归一化处理；

步骤33：计算水平井测井数据的相关性系数并进行相关对比；计算水平井测井数据的曲线活度并进行活度分层；所述相关对比是指不同探测深度测井曲线归一化后数据在固定窗长下的相关性对比；所述活度分层是按以下方式进行的：首先计算曲线活度，判断曲线活度是否大于预设的活度阈值，如果是，则将大于活度阈值的极值点的位置定义为分层界面；

所述相关性系数按以下公式来计算：

$Z_i=\frac{{\mathrm{\Σ}}_{i=1}^n(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{{(x_i-\stackrel{\‾}{x})}^2+{(y_i-\stackrel{\‾}{y})}^2}}$

式中，i表示采样点位置，n表示窗长，表示窗长内对比条曲线平均值，表示被对比条曲线平均值，x_i表示对比曲线第i个采样点值，y_i表示被对比曲线第i个采样点曲线值，Z_i表示相关系数。

5.根据权利要求4所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤31中对标准化的水平井测井数据进行预处理采用最小二乘滑动平均法或加权滑动平均法进行平滑滤波，所述平滑滤波采用以下的方式：

线性函数平滑：

二次函数平滑：

钟形函数平滑法：

汉明函数平滑法：

式中；T_i表示第i个采样点值，m表示滤波窗长，表示滤波后第i各采样点取值。

6.根据权利要求4所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤32中对对测井数据进行归一化处理采用极值归一化方法、孔隙度归一化方法或密度归一化方法进行，归一化方法具体如下：

对电阻率采样极值归一化方法：

对中子孔隙度采样线性归一化方法：

对中子密度采样线性归一化方法：

式中；x_ij表示第j条曲线的第i个采样点值，x_minj表示第j条曲线的最小值，x_maxj表示第j条曲线的最大值，Y_ij表示归一化后第j条曲线的第i个值；TNPL_j表示孔隙度曲线第j个值，用归一化后的值代替原始值；ALCD_j表示中子密度第j个值，用归一化后的值代替原始值。

7.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤33中，对测井数据进行相关对比和活度分层，具体步骤如下：

步骤331：所述活度分层中的活度按照以下公式进行定义：

$E\left(d\right)={\mathrm{\Σ}}_{i=d-n}^{i=d+n}{\[x\left(i\right)-\stackrel{\‾}{x}\left(d\right)\]}^2---\left(12\right)$

$\stackrel{\‾}{x}\left(d\right)=\frac{1}{2n+1}{\mathrm{\Σ}}_{i=d-n}^{i=d+n}x\left(i\right)---\left(13\right)$

式中，E(d)表示d的活度函数值，x(t)表示测井曲线测量值，表示测井曲线在区间[d-n，d+n]内的平均值；d表示采样点位置；n表示采样窗长；i表示在一个采样窗长中第几个采样点；

步骤332：根据活度函数值与预设活度阈值进行比较来确定层界面，所述活度阈值为所有活度极值从大到小排列的前5％处作为活度阈值，大于活度阈值的极值点作为层界面，小于阈值的点不做层界面。

8.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤4中的剔除自动分层产生的假层，具体步骤如下：

步骤41：根据层界面位置，用两个相邻层界面相减的绝对值，即为相邻层界面间垂直距离；

步骤42：判断两相邻界面垂直距离是否小于预设距离值，如果是，则该层界面为无效界面，并视为假层；

步骤43：如果否，则返回步骤41重复直至结束；

其中，在计算出相邻层界面垂直距离后，根据层界面垂直距离，把层界面距离小于预设定值的层厚看作是假层，并剔除假层。

9.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤6中利用有限元方法计算建立的地层初始模型正演响应，具体步骤如下：

步骤61)求解给定边界条件下麦克斯韦Maxwell方程的问题，将Maxwell方程转化为波动方程：

$\▿\×(\▿\×\frac{\stackrel{\→}{E}}{\mathrm{\μ}})-{\mathrm{\ω}}^2\mathrm{\ϵ}\stackrel{\→}{E}=-j\mathrm{\ω}{\stackrel{\→}{J}}_s---\left(14\right);$

其中，E表示电场强度；μ表示振幅；ω表示角频率；ε表示介电常数；

步骤62)结合边界条件波动方程归结为场能量泛函：

$F\left(\stackrel{\→}{E}\right)=\frac{1}{2}\∫\∫{\∫}_V\[\frac{1}{{\mathrm{\μ}}_r}\left(\▿\×\stackrel{\→}{E}\right)\·\left(\▿\×\stackrel{\→}{E}\right)-{\mathrm{\ω}}^2\mathrm{\ϵ}\stackrel{\→}{E}\×\stackrel{\→}{E}\]dV+{\mathrm{j\ω\μ}}_0\∫\∫{\∫}_V\stackrel{\→}{E}\×\stackrel{\→}{J}dV---\left(15\right);$

其中，表示目标函数；μ₀表示初始振幅；V表示求解区域；ω表示角频率；ε表示介电常数；J表示电流密度；

步骤63)应用有限单元剖分场域，并选取相应的插值基函数，对能量泛函进行空间离散，得到如下离散化泛函形式：

$F=\frac{1}{2}{\mathrm{\Σ}}_{e=1}^M({\left\{E^e\right\}}^T\[A^e\]\{E^e\}-{\mathrm{\ω}}^2{\mathrm{\ϵ\μ}}_0{\left\{E^e\right\}}^T\[B^e\]\{E^e\}+{\mathrm{j\ω}}_0{\left\{E^e\right\}}^T\{C^e\left\}\right)---\left(16\right)$

其中，A^e表示函数离散形成系数矩阵；B^e表示等式右端项离散矩阵；C^e表示函数值离散矩阵；E^e表示自变量矩阵；M表示自变量个数；T表示矩阵的转置；

步骤64)求解线性方程组得到所需的参数。

10.根据权利要求1所述的水平井地层环境下电阻率各向异性识别方法，其特征在于，在步骤6中利用最小二乘方法进行多参数迭代反演，具体步骤如下：

步骤621)通过对比模型正演响应和实测数据，建立目标函数，通过最小二乘法求解实际测量值与模拟值的残差：

$\underset{x\∈R^n}{min}f\left(x\right)=\frac{1}{2}r{\left(x\right)}^{T}r\left(x\right)=\frac{1}{2}\underset{i=1}{\overset{m}{\Σ}}{\[r_i\left(x\right)\]}^2,m\≥n---\left(17\right)$

其中，r:Rⁿ→R^m是x的非线性函数；R表示自变量向量；m表示未知量个数；n表示正演响应曲线个数；

步骤622)通过梯度方法求解目标函数雅可比矩阵，形成雅可比线性方程组，求解方程组计算目标函数最速下降方向；

按照以下公式来计算目标函数的梯度：

$g\left(x\right)={\mathrm{\Σ}}_{i=1}^m{r}_i\left(x\right)\▿r_i\left(x\right)=J{\left(x\right)}^{T}r\left(x\right)---\left(19\right)$

其中，J(x)是r(x)的雅可比矩阵：

求解目标函数最速下降方向；

步骤623)利用黄金分割方法，在最速下降方向上确定步长；

步骤624)根据求得模型改变方向和改变步长，确定模型改变量，改变模型，完成一次迭代反演；

步骤625)设置迭代终止条件，循环调用迭代反演，直至满足迭代终止条件，输出结果。