CN105938503A - 一种方向信号多层界面识别方法 - Google Patents

Abstract

本发明公开了一种使用方向信号识别多个地层层界面的方法，首先根据测井获得方位自然GR和方向信号并提取地层倾角和地层方位角；然后测井获得多个源距和不同频率下的仪器方向信号测井响应值，得到标准化的测井数据和真实层界面；构建地层模型；最后计算地层模型输出反演结果绘制井眼轨迹‑地层界面的关系图；该方法应用方向信号识别多个层界面能够获得井眼‑地层全面地层信息，完善了现有测井评价资料，在地质导向中，可提前获得未钻遇层电阻率剖面信息和储层位置关系，为实时钻进过程中钻头调整提供数据支撑；在后期评价中，通过层边界距和每层电阻率反演结果可精确计算储集层厚度、储层含油饱和度等，为储层评价提供可靠参数。

Description

一种方向信号多层界面识别方法

技术领域

本发明涉及油田开发技术领域，尤其涉及一种方向信号多层界面识别方法。

背景技术

随着复杂油气田勘探开发的不断深入，大斜度井、水平井等复杂工艺井的广泛应用，随钻测井技术研究与随钻测井仪器研发得到足够的重视，并得到了快速的发展。随钻电磁波测井仪器由于其较深的探测深度，能够提前预测层界面的存在，在层界面的预测以及地质导向中有着重要的作用。但是由于传统的随钻电磁波测井仪器发射线圈与接收线圈共轴，测量得到的地层信号为地层信息的平均值，不具备方位特性，当仪器接近或离开储层界面时，难以确定储层界面位于仪器的上方还是下方，给钻头调整带来了不确定性，不利于钻井过程中的地质导向。为了提供更完善的钻前预测，提供实时的地质导向信息，三大石油测井服务公司相继推出了具有方位探测能力的随钻方位电磁波测井仪器。2005年斯伦贝谢公司推出了PeriScope方位电阻率测量仪，2006年贝克休斯公司推出了随钻方位电磁波测井仪APR，2007年哈里伯顿公司推出了方位深电阻率测量仪ADR。方位电磁波测量仪器均采用轴向，倾斜或横向线圈混合，能够更好地提供有关地层方位的信息，指示地层的各向异性，识别地层边界。

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

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

随钻电缆联合反演以阻尼最小二乘法为基础，将随钻与电缆测井响应、随钻与电缆测井时刻地层模型参数有机统一起来进行反演。

首先建立二维快速反演所需的数据库。随钻和电缆的反演数据库的建立，都是以它们的正演程序为基础，计算不同井径、泥浆电阻率、围岩电阻率、目的层电阻率、侵入带电阻率、层厚和侵入半径等因素下的视电阻率值，并将各个模型下计算得到的结果保存成数据库(但是在此过程中，所有模型均是基于二维旋转对称模型，无法考虑水平井和大斜度井下非对称地层，无法考虑层边界影响，也不能进行多边界识别)；数据库建立好后，可进行二维快速反演为联合反演赋初值。二维快速反演过程中，由于库中数据量非常大，若先读取整个数据然后进行查询，花费的时间和内存将非常大，所以一般使用定位读取某段数据库的方法，不仅使得读取的库变小，也加快了查询的速度。利用二维快速反演得到联合反演初值(得到径向侵入剖面模型，包括侵入半径、侵入带电阻率和地层电阻率)后即可进行联合反演。流程如下：

(1)将随钻与电缆测井数据进行深度对齐(实际的随钻与电缆测井数据往往存在深度上的差异)，读入随钻与电缆测井数据文件，加载随钻与电缆测井数据；

(2)对测井曲线进行自动分层，导入自动分层结果文件，定义层界面位置；

(3)按照二维快速反演程序对测井数据进行查库反演，导入查库反演结果文件，给联合反演初始地层模型参数赋初值；

(4)分别调用随钻、电缆正演程序，计算设置的地层模型下的随钻、电缆仪器响应；

(5)判断计算出的响应值是否充分逼近实测数据，若是，则输出该模型值作为反演结果，若否，则使用最小二乘迭代法计算模型参数改变量，重置地层模型参数，重复步骤(4)、(5)，直至输出反演结果(地层电阻率等)。

但是，传统电阻率联合反演主要是针对旋转对称性地层进行一维和(或)二维反演，既不适用非对称地层结构的水平井和大斜度井，也不能提供层边界距离和各向异性等信息。

因此，需要一种能提供层边界距离和各向异性等信息的基于三维空间地层模型的方向信号多层边界识别方法。

发明内容

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

(二)技术方案

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

本发明提供的一种方向信号多层界面识别方法，包括以下步骤：

步骤1，根据测井获得方位自然GR和方向信号并提取地层倾角和地层方位角；

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

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

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

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

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

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

步骤8，根据反演结果绘制井眼轨迹-地层界面的关系图，并输出关系图。

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

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

$\begin{array}{c}y=A\sin (\mathrm{\ω}x-\mathrm{\β})+y_0\\ =y_0+A\cos \mathrm{\β}\sin \left(\mathrm{\ω}x\right)+(-A\sin \mathrm{\β})\cos \left(\mathrm{\ω}x\right)\end{array}---\left(1\right);$

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

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

其中：

α₀＝y₀，α₁＝Acosβ，α₂＝-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);$

倾角倾向表达式：

$\{\begin{array}{c}Dip=\arctan \left(\frac{2A}{D_{el}}\right)\\ Dir=x|y=y_{\min }\end{array}---\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：对归一化测井数据进行相关对比和活度分层。

进一步，在步骤3中对标准化的测井数据进行预处理采用最小二乘滑动平均法或加权滑动平均法进行平滑滤波；对测井数据进行归一化处理采用极值归一化方法、孔隙度归一化方法或密度归一化方法进行。

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

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

$E\left(d\right)={\Σ}_{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}{\Σ}_{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表示在一个采样窗长中第几个采样点；

步骤322：根据活度函数值与预设活度阈值进行比较来确定层界面。

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

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

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

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

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

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

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

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

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

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

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

$F=\frac{1}{2}{\Σ}_{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中利用最小二乘方法进行多参数迭代反演，具体步骤如下：

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的非线性函数；R表示自变量向量；m表示未知量个数；n表示正演响应曲线个数；

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

则目标函数的梯度为：

$g\left(x\right)={\Σ}_{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)的雅可比矩阵：

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

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

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

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

进一步，在步骤7中，根据反演结果，对测量点对应位置按测井曲线间隔对反演结果重采样，同时对逐层反演的方波形结果进行非线性插值，获得光滑的结果曲线；以及

在步骤8中，根据反演结果输出文件，绘制原始曲线测井响应、反演结果曲线、井眼轨迹-地层关系图并对地层进行色标填充，填充颜色根据地层岩性和电阻率值，按国际通用色标填充。

(三)有益效果

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

本发明采用方向信号来识别多个地层层界面，该方法应用方向信号识别多个层界面能够获得井眼-地层相对倾角、储层厚度、测量点到多个层界面间距离以及每层对应的电阻率信息等，完善了现有测井评价资料，可得到全面地层信息。在地质导向中，可提前获得未钻遇层电阻率剖面信息和储层位置关系，为实时钻进过程中钻头调整提供数据支撑；在后期评价中，通过层边界距和每层电阻率反演结果可精确计算储集层厚度、储层含油饱和度等，为储层评价提供可靠参数。

附图说明

图1为本发明的方向信号多层界面识别方法具体实施例的流程图；

图2为本发明的方向信号在层界面处响应曲线图；

图3为本发明的方向信号随地层倾角变化图；

图4为本发明的有限元对地层模型正演仿真模拟视电阻率图；

图5为本发明的大斜度井地层建模示意图。

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

图7为本发明的实施例中对某段实测井资料的多边界反演结果图。

具体实施方式

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

实施例1

如图所示，图1为本发明的方向信号多层界面识别方法具体实施例的流程图；本实施例提供一种方向信号多层界面识别方法，具体包括以下步骤：

步骤1，根据测井获得方位自然GR和方向信号并提取地层倾角和地层方位角；

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

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

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

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

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

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

步骤8，根据反演结果绘制井眼轨迹-地层界面的关系图，并输出关系图。

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

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

$\begin{array}{c}y=A\sin (\mathrm{\ω}x-\mathrm{\β})+y_0\\ =y_0+A\cos \mathrm{\β}\sin \left(\mathrm{\ω}x\right)+(-A\sin \mathrm{\β})\cos \left(\mathrm{\ω}x\right)\end{array}---\left(1\right);$

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

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

其中：

α₀＝y₀，α₁＝Acosβ，α₂＝-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);$

倾角倾向表达式：

$\{\begin{array}{c}Dip=\arctan \left(\frac{2A}{D_{el}}\right)\\ Dir=x|y=y_{\min }\end{array}---\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：对测井数据进行光滑滤波处理，剔除测井数据中无效的数据，所述无效数据分为两类，一类是明显标记为的无效数据，比如值为-999.25和等于仪器测量上限2000.00的值，对这类无效数据直接删除，然后用临近有效值进行插值，对删除的无效值位置重采样；第二类无效数据是因为位置原因的跳变，个别采样点异常，该类数据并非完全无效，只是有不同程度失真，对该类数据光滑滤波处理，滤波的方式有以下几种可供选择：

线性函数平滑：

二次函数平滑：

钟形函数平滑法：

汉明函数平滑法：

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

式中，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：进行相关对比和活度分层，相关对比是不同探测深度(不同源距或者不同频率)测井曲线归一化后数据在固定窗长下相关性对比，计算相关性系数，相关性系数越大表示相关性越好；活度分层的方法是先计算曲线活度(见步骤32)，设置活度阈值，大于活度阈值的极值点位置定义为分层界面，活度阈值可调，通常设置为所有活度极值从大到小排列的前5％处作为活度阈值。

计算相关性公式：

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

在步骤3中对标准化的测井数据进行预处理采用最小二乘滑动平均法或加权滑动平均法进行平滑滤波；对测井数据进行归一化处理采用极值归一化方法、孔隙度归一化方法或密度归一化方法进行。

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

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

$E\left(d\right)={\Σ}_{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}{\Σ}_{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表示在一个采样窗长中第几个采样点；

步骤322：根据活度函数值与预设活度阈值进行比较来确定层界面，活度阈值并不固定，这与仪器的信噪比有关，对于不同的仪器有不同的阈值，通常设置为所有活度极值从大到小排列的前5％处作为活度阈值，大于阈值的极值点作为层界面，小于阈值的点不做层界面。

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

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

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

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

在计算出相邻层界面距离后，根据层界面距离，把层界面距离小于一定值(通常设置为仪器分辨率)的层厚看作是假层，去掉其中的一个(活度值较小的那个)层界面来剔除假层。

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

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

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

其中，E表示电场强度；μ表示振幅；ω表示角频率；ε表示介电常数；这是电磁波在地层中传播的理论基础，地层模型不同，电磁波所感应出的电磁场也不同，通过电磁场的不同反过来也可以推导出地层参数。影响电磁波在地层中传播的因素主要有：地层电阻率、层厚、各向异性系数、井径、泥浆电阻率、收发线圈距离、倾角等。

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

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

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

$F=\frac{1}{2}{\Σ}_{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中利用最小二乘方法进行多参数迭代反演，具体步骤如下：

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的非线性函数；R表示自变量向量；m表示未知量个数；n表示正演响应曲线个数；

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

则目标函数的梯度为：

$g\left(x\right)={\Σ}_{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)的雅可比矩阵：

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

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

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

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

在步骤7中，根据反演结果，对测量点对应位置按测井曲线间隔对反演结果重采样，同时对逐层反演的方波形结果进行非线性插值，获得光滑的结果曲线；以及

在步骤8中，根据反演结果输出文件，绘制原始曲线测井响应、反演结果曲线、井眼轨迹-地层关系图并对地层进行色标填充，填充颜色根据地层岩性和电阻率值，按国际通用色标填充。

实施例2

本发明提供一种使用方向信号识别多个地层层界面的方法，该方法包括：

步骤1，测井获得方位自然GR和方向信号提取地层倾角和方位角；

步骤2，测井获得多个源距和不同频率下的仪器方向信号测井响应值，测井获得中子密度和孔隙度曲线，对未进行方位校正的方向信号曲线进行方位校正；

步骤3，将步骤1、步骤2获得的测井响应曲线，进行相关对比，应用活度方法自动分层；

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

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

步骤6，计算步骤5构建的初始地层模型的方向信号正演响应，与实测曲线对比，用最小二乘方法进行多参数迭代反演；

步骤7，格式化输出反演结果，反演结果中根据井眼与地层界面关系，输出测量点分别到多个边界的边界距；

步骤8，根据输出结果，绘制井眼轨迹-地层关系图。

本实施例提供的方法应用方向信号识别多个层界面，能够获得井眼-地层相对倾角、储层厚度、测量点到多个层界面间距离以及每层对应的电阻率信息等，完善了现有测井评价资料，可得到全面地层信息。

实施例3

本实施例提供的方向信号多层界面识别方法，包括以下步骤：

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

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

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

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

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

其中：α₀＝y₀，α₁＝Acosβ，α₂＝-Asinβ，

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)$

和倾角倾向表达式：

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

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

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

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

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

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

本实施例利用相关对比和活度法进行地层自动分层，首先对测井数据光滑滤波和无效数据剔除；然后数据归一化处理；最后再进行相关对比和活度分层。

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

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

$\stackrel{\‾}{T_i}=\frac{1}{35}(-3(T_{i-2}+T_{i+1})+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)$

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

极值归一化方法：

孔隙度归一化方法：

密度归一化方法：

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

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

$E\left(d\right)={\Σ}_{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}{\Σ}_{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，结合倾角、界面、视电阻率、井眼轨迹垂深-斜深-水平位移变化，建立初步三维多层地质模型；

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

其中，图5为本发明的大斜度井地层建模示意图，B1表示第一个边界位置，B2表示第二个边界位置，Rm表示泥浆电阻率，Ri表示侵入深度，Rxo表示侵入带电阻率，Rs1表示上围岩电阻率，Rs2表示下围岩电阻率，Rt(Rh-Rv)表示目的层电阻率(水平电阻率-垂直电阻率)；图6为本发明的水平井地层建模示意图，图中DTB1表示到上界面距离，Ri表示侵入深度，Rxo表示侵入带电阻率，Rm表示泥浆电阻率，Rt(Rh-Rv)表示目的层电阻率(水平电阻率-垂直电阻率)。

其中，图7为本发明的实施例中对某段实测井资料的多边界反演结果图。

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

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

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

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

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

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

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

$F=\frac{1}{2}{\Σ}_{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)={\Σ}_{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，格式化输出反演结果；

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

步骤8，根据步骤7处理得到的反演结果，绘制井眼轨迹-地层关系图，还原地层界面位置关系，进而实现多边界识别；

本实施例根据反演结果输出文件，绘制原始曲线测井响应、反演结果曲线、井眼轨迹-地层关系图并对地层进行色标填充，填充颜色根据地层岩性和电阻率值，按国际通用色标填充。

实施例4

本实施例提供一种可以进在水平井和大斜度井等复杂井眼环境下利用方向信号进行多边界识别的方法；通过如下步骤来实现：

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

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

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

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

其中：α₀＝y₀，α₁＝Acosβ，α₂＝-Asinβ，

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)$

和倾角倾向表达式：

$\{\begin{array}{c}Dip=\arctan \left(\frac{2A}{D_{el}}\right)\\ Dir=x|y=y_{\min }\end{array}---\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中，计算构建的初始地层模型的方向信号正演响应，与实测曲线对比，用最小二乘方法进行多参数迭代反演；利用有限元方法计算建立的地层模型响应的步骤，包括以下几个步骤：

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

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

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

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

$F=\frac{1}{2}{\Σ}_{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)求解这个线性方程组得到所需的参数。

本实施例中的利用有磁偶极子解析解方法计算建立的地层模型响应，步骤如下：

设单位磁偶极子源随时间的变化关系为exp(iωt)，其中ω为角频率，并假设在地层直角坐标系中(水平面为xy面)源点的位置坐标为r_t＝(x_t，y_t，z_t)，场点的位置坐标为r＝(x，y，z)，则方向单位磁偶极子在均匀各向异性介质中产生的Hertz势可表示为:

${\Π}^{\left(3\right)}(r,r_t)=\left(\exp \right(-{\mathrm{ik}}_{h}R)/4\mathrm{\π}R){\hat{e}}_z---\left(17\right)$

式中，μ_b为均匀介质磁导率，σ_hb为均匀各向异性介质的水平复电导率。经过转化处理，(17)式可以表示为如下Sommerfeld积分形式：

${\Π}^{\left(3\right)}(r,r_t)={\hat{e}}_z\frac{1}{4\mathrm{\π}}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\λ}}{{\mathrm{\Λ}}_h}\exp (-{\mathrm{\Λ}}_h|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}---\left(18\right)$

式中，l_v为v阶Bessel函数，λ为积分变量。方向单位磁偶极子在均匀各向异性介质中产生的Hertz势可以表示为：

${\Π}^{\left(1\right)}(r,r_t)=\frac{\exp (-{\mathrm{ik}}_{v}S)}{4\mathrm{\π}KS}{\hat{e}}_x+\frac{(x-x_t)(z-z_t)}{4\mathrm{\π}{\tilde{r}}^2}(K\frac{\exp (-{\mathrm{ik}}_{v}S)}{S}-\frac{\exp (-{\mathrm{ik}}_{h}R)}{R}){\hat{e}}_z---\left(19\right)$

${\Π}^{\left(2\right)}(r,r_t)=\frac{\exp (-{\mathrm{ik}}_{v}S)}{4\mathrm{\π}KS}{\hat{e}}_y+\frac{(y-y_t)(z-z_t)}{4\mathrm{\π}{\tilde{r}}^2}(K\frac{\exp (-{\mathrm{ik}}_{v}S)}{S}-\frac{\exp (-{\mathrm{ik}}_{h}R)}{R}){\hat{e}}_z---\left(20\right)$

式中，为各向异性系数，σ_vb为均匀各向异性介质的垂向复电导率。经推导(19)、(20)式可以分别表示为如下Sommerfeld积分形式：

$\begin{array}{c}{\Π}^{\left(1\right)}(r,r_t)={\hat{e}}_x\frac{1}{4\mathrm{\π}K}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\λ}}{{\mathrm{\Λ}}_v}\exp (-{\mathrm{\Λ}}_{v}K|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}\\ +{\hat{e}}_z\frac{1}{4\mathrm{\π}}\frac{z-z_t}{\left|z-z_t\right|}\frac{\∂}{\∂x}{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\λ}}\left(\exp \right(-{\mathrm{\Λ}}_{v}K\left|z-z_t\right|)-\exp (-{\mathrm{\Λ}}_h\left|z-z_t\right|\left)\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}\end{array}---\left(21\right)$

$\begin{array}{c}{\Π}^{\left(2\right)}(r,r_t)={\hat{e}}_y\frac{1}{4\mathrm{\π}K}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\λ}}{{\mathrm{\Λ}}_v}\exp (-{\mathrm{\Λ}}_{v}K|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}\\ +{\hat{e}}_z\frac{1}{4\mathrm{\π}}\frac{z-z_t}{\left|z-z_t\right|}\frac{\∂}{\∂x}{\∫}_0^{\mathrm{\∞}}\frac{1}{\mathrm{\λ}}\left(\exp \right(-{\mathrm{\Λ}}_{v}K\left|z-z_t\right|)-\exp (-{\mathrm{\Λ}}_h\left|z-z_t\right|\left)\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}\end{array}---\left(22\right)$

式中，

由电磁场与Hertz势之间的关系式

及

其中，为均匀各向异性介质的电导率张量。将(17)、(19)、(20)式带入(23)式，可以得到沿三个方向单位磁偶极子产生的电场和磁场各分量的解析式及Sommerfeld积分形式，其中电场和磁场z分量的Sommerfeld积分形式分别表示为：

$E_z^{\left(1\right)}(r,r_t)=\frac{{\mathrm{i\ω\μ}}_{b}K}{4\mathrm{\π}}\frac{\∂}{\∂y}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\λ}}{{\mathrm{\Λ}}_v}\exp (-{\mathrm{\Λ}}_{v}K|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}---\left(20a\right)$

$E_z^{\left(2\right)}(r,r_t)=-\frac{{\mathrm{i\ω\μ}}_{b}K}{4\mathrm{\π}}\frac{\∂}{\∂x}{\∫}_0^{\mathrm{\∞}}\frac{\mathrm{\λ}}{{\mathrm{\Λ}}_v}\exp (-{\mathrm{\Λ}}_{v}K|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}---\left(20b\right)$

$E_z^{\left(3\right)}(r,r_t)=0---\left(20c\right)$

$E_z^{\left(1\right)}(r,r_t)=-\frac{1}{4\mathrm{\π}}\frac{z-z_t}{\left|z-z_t\right|}\frac{\∂}{\∂x}{\∫}_0^{\mathrm{\∞}}\mathrm{\λ}\exp (-{\mathrm{\Λ}}_h|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}----\left(20d\right)$

$E_z^{\left(2\right)}(r,r_t)=-\frac{1}{4\mathrm{\π}}\frac{z-z_t}{\left|z-z_t\right|}\frac{\∂}{{\∂}_y}{\∫}_0^{\mathrm{\∞}}\mathrm{\λ}\exp (-{\mathrm{\Λ}}_h|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}---\left(20e\right)$

$H_z^{\left(3\right)}(r,r_t)=\frac{1}{4\mathrm{\π}}{\∫}_0^{\mathrm{\∞}}\frac{{\mathrm{\λ}}^3}{{\mathrm{\Λ}}_h}\exp (-{\mathrm{\Λ}}_h|z-z_t\left|\right)J_0\left(\mathrm{\λ}\tilde{r}\right)d\mathrm{\λ}---\left(20f\right)$

上述各分量均被表示成了波模积分的形式。例如

$H_z^{\left(1\right)}(r,r_t)={\∫}_0^{\mathrm{\∞}}{\tilde{H}}_z^{\left(1\right)}\left(\mathrm{\λ}\right)d\mathrm{\λ}$

式中，为某一λ对应的波模。

根据Maxwell方程组，电场和磁场波模的切向分量与纵向分量之间的关系可表示为：

${\tilde{H}}_x=-\frac{1}{K^2{\mathrm{\λ}}^2}{\mathrm{i\ω\ϵ}}_h\frac{\∂{\tilde{E}}_z}{\∂y}+\frac{1}{{\mathrm{\λ}}^2}\frac{\∂}{\∂z}\frac{\∂{\tilde{H}}_z}{\∂x}---\left(21a\right)$

${\tilde{H}}_y=\frac{1}{K^2{\mathrm{\λ}}^2}{\mathrm{i\ω\ϵ}}_h\frac{\∂{\tilde{E}}_z}{\∂x}+\frac{1}{{\mathrm{\λ}}^2}\frac{\∂}{\∂z}\frac{\∂{\tilde{H}}_z}{\∂y}---\left(21b\right)$

${\tilde{E}}_x=\frac{1}{{\mathrm{\λ}}^2}i\mathrm{\ω}\mathrm{\μ}\frac{\∂{\tilde{H}}_z}{\∂y}+\frac{1}{K^2{\mathrm{\λ}}^2}\frac{\∂}{\∂z}\frac{\∂{\tilde{E}}_z}{\∂x}---\left(21c\right)$

${\tilde{E}}_y=-\frac{1}{{\mathrm{\λ}}^2}i\mathrm{\ω}\mathrm{\μ}\frac{\∂{\tilde{H}}_z}{\∂x}+\frac{1}{K^2{\mathrm{\λ}}^2}\frac{\∂}{\∂z}\frac{\∂{\tilde{E}}_z}{\∂y}---\left(21d\right)$

将(20)式代入(21)式即可分别计算得到电场和磁场x、y分量的Sommerfeld积分表达式。

本实施例中的利用最小二乘方法进行多参数反演的步骤，包括以下几个步骤：

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)={\Σ}_{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)设置迭代终止条件，循环调用迭代反演，直至满足迭代终止条件，输出结果。

利用有限元方法进行方位随钻电磁波仪器的数值仿真，研究了仪器的方向信号同地层界面位置和地层电阻率的关系，研究结果表明方向信号对层界面由很好的指示作用，在层界面处，方向信号响应幅度值达到极大值，随着仪器远离层界面，方向信号值减弱，在足够远的地方，测量点与层边界距离大于仪器探测深度，方向信号响应幅度等于0。此时，综合使用地层电阻率和方向信号进行多边界识别是有效的。因此，本发明提出利用仪器测量得到的方向信号和电阻率曲线进行多边界识别的方法，并给出电阻率剖面和井眼轨迹-地层位置关系图，在地质导向和钻后评价中均具有重要意义。本发明在不增加现有测井技术成本的前提下完成对地层多边界的边界距与地层电阻率的提取，更为行之有效的反应地层真实信息，完善了现有测井评价资料；本发明使用同时三维反演得到的多层边界距与地层电阻率，更加贴近真实地层状况，得到全面地层参数信息，进行更为真实有效的地层评价和储量计算。

步骤7，格式化输出反演结果，反演结果中根据井眼与地层界面关系，输出测量点分别到多个边界的边界距；

步骤8，根据输出结果，绘制井眼轨迹-地层关系图。

实施例5

如图1所示，图1为本发明的方向信号多层界面识别方法的流程图。

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

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

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

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

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

步骤106，方向信号多层界面识别是基于正演仿真和迭代反演完成的。因此首要的是建立正演仿真的三维地层模型，其建立步骤如下：

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

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

(3)利用有限元方法计算上述建立的地层模型响应。在一实施例中，利用有限元方法计算正演模型响应，包含以下几个步骤：

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

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

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

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

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

$F=\frac{1}{2}{\Σ}_{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)={\Σ}_{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)设置迭代终止条件，循环调用迭代反演，直至满足迭代终止条件，输出结果。

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

步骤108，根据反演结果输出文件，绘制原始曲线测井响应、反演结果曲线、井眼轨迹-地层关系图并对地层进行色标填充，填充颜色根据地层岩性和电阻率值，按国际通用色标填充。流程结束。

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

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

Claims (10)

1.一种方向信号多层界面识别方法，其特征在于，包括以下步骤：

步骤1，根据测井获得方位自然GR和方向信号并提取地层倾角和地层方位角；

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

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

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

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

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

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

步骤8，根据反演结果绘制井眼轨迹-地层界面的关系图，并输出关系图。

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：对归一化测井数据进行相关对比和活度分层。

5.根据权利要求4所述的方向信号多层界面识别方法，其特征在于，在步骤3中对标准化的测井数据进行预处理采用最小二乘滑动平均法或加权滑动平均法进行平滑滤波；对测井数据进行归一化处理采用极值归一化方法、孔隙度归一化方法或密度归一化方法进行。

6.根据权利要求4所述的方向信号多层界面识别方法，其特征在于，在步骤3中，对数据进行相关对比和活度分层，具体步骤如下：

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

$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表示在一个采样窗长中第几个采样点；

步骤322：根据活度函数值与预设活度阈值进行比较来确定层界面。

7.根据权利要求1所述的方向信号多层界面识别方法，其特征在于，在步骤4中的剔除自动分层产生的假层，具体步骤如下：

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

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

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

8.根据权利要求10所述的方向信号多层界面识别方法，其特征在于，在步骤6中利用有限元方法计算建立的地层初始模型正演响应，具体步骤如下：

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

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

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

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

其中，表示目标函数；μ₀表示初始振幅；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)求解这个线性方程组得到所需的参数。

9.根据权利要求10所述的方向信号多层界面识别方法，其特征在于，在步骤6中利用最小二乘方法进行多参数迭代反演，具体步骤如下：

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的非线性函数；R表示自变量向量；m表示未知量个数；n表示正演响应曲线个数；

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

则目标函数的梯度为：

$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)的雅可比矩阵：

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

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

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

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

10.根据权利要求1所述的方向信号多层界面识别方法，其特征在于，在步骤7中，根据反演结果，对测量点对应位置按测井曲线间隔对反演结果重采样，同时对逐层反演的方波形结果进行非线性插值，获得光滑的结果曲线；以及

在步骤8中，根据反演结果输出文件，绘制原始曲线测井响应、反演结果曲线、井眼轨迹-地层关系图并对地层进行色标填充，填充颜色根据地层岩性和电阻率值，按国际通用色标填充。