CN115629428A - 一种泥浆侵入下电阻率及核测井联合反演储层参数的方法 - Google Patents

Abstract

本发明是一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，可用于模拟泥浆侵入条件下任意时刻的测井曲线及储层物理参数反演。首先构建井筒数值模型，然后根据地层储层特征对井筒数值模型进行赋值，再对赋值后的数值模型进行数值模拟仿真得到储层的电阻率响应、中子响应和密度响应曲线，将实测测井曲线与数值模拟响应曲线进行对比，不断调整数值模型参数的赋值直至数值模拟曲线与实测测井曲线基本一致，最后得到地层条件下的真实储层参数。该方法基于岩石物理参数结合井筒信息及相渗特征曲线，从仪器原理出发综合渗流特征、电磁场及核测井原理，通过多时刻、多物理场多维度测井响应信息反演储层实际储层参数，提高了储层流体性质识别及参数评价精度。

Description

一种泥浆侵入下电阻率及核测井联合反演储层参数的方法

技术领域

本发明涉及石油天然气勘探开发储层测井评价技术领域，具体地，涉及一种泥浆侵入下电阻率及核测井联合反演储层参数的方法。

背景技术

目前国内低渗油气田储量规模较大，但是该类油气藏往往由于局部泥质重、灰质重或是孔隙结构复杂，导致储层物性差，因此自然产能低。人工开采时，钻井在受泥浆侵入后，测井曲线响应受到多种因素影响，导致油气、水层电阻率差异小，给储层评价带来了较大的困难。

目前主要是采用常规流体识别及储层参数测量针对不同的目标采用不同的测井项目进行模拟，如现有的一种基于并行计算技术的随钻电阻率测井联合反演方法，包括以下步骤：A.载入随钻电阻率测井数据；B.对测井曲线光滑滤波、自动分层；C.根据曲线分离程度选取地层模型侵入深度、侵入带电阻率、地层电阻率初始值；D.采用有限元素法，基于多核并行或GPU并行计算技术，计算实测曲线相应深度点的随钻电阻率测井响应；E.计算正演仿真响应值与实测数据相对误差，若相对误差小于设定阈值，则输出该模型值作为反演结果，若否，则根据马奎特迭代算法计算模型改变量，重置地层模型参数，返回步骤D，直至输出反演结果。

但是，采用常规流体识别及储层参数测量针对不同的目标采用不同的测井项目进行模拟得到的流体性质与孔隙结构参数精度低。

发明内容

本发明为解决上述技术方案的问题提供了一种泥浆侵入下电阻率及核测井联合反演储层参数的方法。本方案能够开展基于岩石物理参数结合井筒信息及相渗特征曲线的电阻率及核测井联合模拟，并通过多时刻、多物理场等多维度测井响应信息联合反演得到可靠的储层参数，提高了储层流体性质识别及参数评价精度，为油气勘探开发提供技术支持。

本发明采用的技术方案是：一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，包括以下步骤：

步骤一：构建井筒数值模型；根据储层伽玛、电阻率、中子、密度曲线特征对地下储层实测测井曲线模型进行分层，然后按层进行所述地下储层实测测井曲线方波化，每一层代表一个岩石物理模型；

步骤二：井筒数值模型初始值设定；将岩石物理模型参数、渗流信息和井眼信息参数输入所述步骤一中的所述岩石物理模型，其中岩石物理模型参数和渗流信息参数是按层赋值，每一层按地层储层特征将参数输入所述步骤一中的所述的岩石物理模型，井眼信息参数按照全井段处理，即输入统一的数值；

步骤三：动态数值模拟；根据相关多孔介质渗流理论和扩散理论对所述步骤二获得的所述井筒初始数值模型进行钻井液动态侵入模拟，获得任意时刻的储层动态流体分布剖面，包括储层动态饱和度和矿化度剖面，通过相关公式对所获得的储层动态流体分布剖面进行转换，获得相应的电阻率、迁移长度、密度在不同时刻的径向分布剖面；

步骤四：电阻率、中子、密度测量响应数值模拟；根据电磁场有限元仿真方法结合仪器结构参数及工作原理获得数值模拟的电阻率测井响应，根据仪器结构参数及核测井快速模拟方法获得数值模拟的中子、密度测量曲线；

步骤五：反演储层参数；用所述步骤一中的所述地下储层实测测井曲线与所述步骤四获得的数值模拟曲线进行对比，通过模型更新及迭代，使所述数值模拟曲线与所述实测测井曲线一致，最终得到反演的岩石物理结构的储层参数。

在上述技术方案中，首先构建井筒数值模型，然后根据地层储层特征对井筒数值模型进行赋值及数值模拟，然后将实测测井曲线与数值模拟曲线进行对比，不断调整数值模型参数的赋值直至数值模拟曲线与实测测井曲线基本一致，进而得到地层条件下真实的储层参数。本方案的处理方法能够开展电阻率及核测井联合动态模拟及反演，得到储层准确的流体性质与孔隙结构参数。

优选的，在步骤一中，测井曲线是对地层物理性质随井深变化的记录，以自然伽玛半幅点，参考电阻率、中子、密度曲线划分各小层顶、底界面。方波化是是指在图形显示方式下对测井曲线的数值分段进行平均，形成锯齿形态，采用每一小层的测井曲线平均值代表该小层。

优选的，在步骤二中，输入模型参数包括岩石物理模型参数、渗流信息和井眼信息参数，岩石物理模型参数和渗流参数是按层处理，井眼信息参数按照全井段处理，岩石物理模型具体参数包括岩石骨架参数、地层物性参数、地层孔隙流体参数、泥质含量。岩石骨架参数包括矿物组分、含量和密度；地层物性参数包括饱和度、渗透率、孔隙度和孔隙结构；地层孔隙流体参数包括流体组分、含量、密度、矿化度和粘度；渗流信息包括相渗曲线和毛管压力曲线；井眼信息参数包括井眼直径、泥浆信息、钻柱压力和温度。岩石骨架参数通过X衍射资料获取，包括流体组分、密度、矿化度、粘度的孔隙流体信息通过水分析实验材料获得，渗流信息通过地区驱替实验获得，井眼信息通过钻井日志、地质日报获得，地层孔隙结构参数通过电阻率岩心实验获得，泥质含量、孔隙度、渗透率、饱和度通过每一小层的测井曲线平均值采用如下的测井解释公式获得。泥质含量初始值获取公式为

式中,I_sh表示泥质含量指数，小数；V_sh表示地层泥质含量；小数；GR表示地层自然伽马测量值，gAPI；GR_max表示本井纯泥岩的自然伽马最大值，gAPI；GR_min本井纯砂岩的自然伽马最小值，gAPI。

孔隙度初始值获取公式为

式中：

表示计算的密度孔隙度，v/v；DEN表示密度曲线测井读值，g/cm3；DEN_ma表示骨架的密度响应参数，g/cm3；DEN_f表示流体的密度响应参数，g/cm3；V_sh表示地层泥质含量，v/v；DEN_sh表示泥岩的密度响应参数，g/cm3；

表示计算的中子孔隙度，v/v；CNC表示中子曲线测井读值，％；CNC_ma表示骨架的中子响应参数，％；CNC_f表示流体的中子响应参数，％；CNC_sh表示泥岩的中子响应参数，％；

表示地层孔隙度，v/v。

渗透率计算可采用区域岩心孔渗回归关系模型，例如，研究区模型公式为：

式中：K表示地层渗透率，mD；

表示地层孔隙度，v/v。

含水饱和度初始值获取公式为：

式中：S_w表示地层含水饱和度；R_t表示地层真电阻率，Ω.m；V_sh表示地层泥质含量，小数；R_sh表示纯泥岩电阻率，Ω.m；R_w表示地层水电阻率，Ω.m；

表示地层孔隙度，小数；a表示岩性系数；m表示胶结指数；n表示饱和度指数，a，m，n采用区块经验值或者其他测井解释方法获得。

优选的，在步骤三中，相渗曲线及毛管压力曲线通过岩心实验获得。根据包括相渗曲线及毛管压力曲线的相关渗透率理论和离子扩散理论模拟获得任意时刻包括地层饱和度及矿化度的储层流体动态分布剖面，通过电阻率及核测井相关公式转换获得相应的电阻率、迁移长度、密度在不同时刻的径向分布剖面。电阻率径向剖面转化公式为：

式中：R_t表示地层电阻率；R_w表示地层水电阻率；S_w表示地层含水饱和度；

表示地层孔隙度；R_sh表示泥质电阻率；V_sh表示泥质含量；a表示岩性系数；m表示胶结指数；n表示饱和度指数；c表示泥质含量指数系数，一般计算公式为

a，m，n采用区块经验值或者其他测井解释方法获得。

密度径向剖面转换公式为：

式中：ρ_b表示地层密度；ρ_w表示地层水密度；ρ_h表示油气流体密度；

表示地层孔隙度；ρ_sh表示泥质密度；C_sh表示泥质孔隙度；S_w表示地层含水饱和度；ρ_i骨架中第i种组分的密度；C_i表示骨架中第i种组分的含量；n_c表示骨架中组分种类总数。

中子迁移长度径向剖面计算公式为：

式中：ξ表示地层迁移长度；S_w表示地层含水饱和度；ξ_w表示地层水迁移长度；ξ_h表示油气流体迁移长度；

表示地层孔隙度；ξ_m表示骨架迁移长度；η表示指数系数。

优选的，在步骤四中，根据步骤三转换的电阻率径向分布剖面获取电阻率测井响应，主要采用电磁场有限元仿真方法结合仪器结构参数及工作原理数值模拟获得对应的电阻率测井响应。其中采用三维矢量有限元素法，对测井物理模型进行空间离散，建立单元电磁场离散方程，对全部单元安装、消元、求解，进而得到电测井仪器响应。电阻率测井仿真为求解给定边界条件下的麦克斯韦方程问题。从麦克斯韦方程出发，电磁场满足下面方程：

式中：E表示电场强度；H表示磁场强度；J为源电流密度；ω为角频率；σ为电导率；μ为磁导率。

从上述两个方程出发，可以推导出电场所满足的矢量波动方程为：

为复介电常数；ε＝ε_rε₀，其中ε₀为真空介电常数；ε_r为相对介电常数。

令

E＝E_p+E_s

式中：背景场E_p是当全部空间被电导率为σ₀的介质填充时的电场，它满足方程：

其中，

由上述三式变化可得到矢量波动方程：

背景场E_p通过解析方法计算得到，二次场E_s则由有限元素法计算。上式的解变化平缓，可以利用稀疏一些的网格来进行求解，减少了计算工作量。选取足够大区域，使边界上的电场衰减到近似为0，则上式只需满足边界条件方程：

式中：

为求解区ω的边界，n为其法线方向。

考虑边界条件方程，将矢量波动方程转化为其弱积形式：

式中：N为矢量基函数，V为单个求解元素，Ω为整个求解空间，E_p表示背景场电场，E_s表示二次场电场,ε_c0表示背景场的介电常数，ε表示二次场介电常数，ω为角频率，μ为磁导率。

优选的，在步骤四中，根据步骤三转换的地层密度径向分布剖面和迁移长度径向分布剖面，结合仪器结构参数及核测井快速模拟方法数值模拟获得中子和密度测井响应。

其基本原理是在一定体积内中子、伽马光子随时间变化率等于产生率减去泄漏率和吸收率，可以用下式来描述这一过程。探测器处通量可以表示为每个从源发射的粒子对探测器贡献的总和。

对于密度测井，伽马源发射伽马光子，与原子核发生康普顿散射、光电散射并产生散射伽马射线，探测器主要探测散射伽马射线的计数率。对于中子孔隙度测井，脉冲中子源发射快中子，与原子核发生非弹性散射、弹性散射反应，中子能量降低变为热中子，探测器主要探测热中子的通量。根据上面理论，背景区域通量可以表示为：

N_B(r_R)＝∫dr∫dE∫dΩψ_B(r_s,r,E,Ω)S(r_R,r,E,Ω)

式中N_B(r_R)表示背景区域探测器的通量，n；r_s表示空间中源的位置，无量纲；r_R表示空间中探测器的位置，无量纲；r表示空间中某一点与源的距离，cm；Ω表示粒子的碰撞角度，°；E表示粒子的能量，MeV；ψ_B表示位于r_s处的放射源发射的粒子在r处的通量，n/(cm²·s)。S表示r_R处的探测器响应函数，可以表示为：

S(r_R,r,E,Ω)＝∫dE′∫dΩ′∑_s(r,E,Ω→E′Ω′)Ψ⁺(r_R,r,E′,Ω′)

式中∑_S表示粒子的宏观散射截面，1/cm；Ω′表示粒子碰撞后角度，°；E′表示粒子碰撞后能量，MeV。Ψ⁺为重要性方程。定义背景通量灵敏度方程为：

变化区域对探测器通量的贡献可以表示为：

式中Δ∑为变化区域粒子截面，1/cm；∑_B为背景区域粒子截面，1/cm。探测器通量为背景通量与变化区域通量的叠加，表示为：

对于密度测井，变量设置为密度值如下式所示，C_p为密度系数，无量纲；ρ表示密度，g/cm³；ρ_B表示背景密度，g/cm³。

对于中子测井，变量设置为迁移长度导数，如下式所示，C_N为中子系数；Lm为迁移长度，无量纲；LmB为背景迁移长度，无量纲。

优选的，在步骤五中，用地下储层实际测井曲线与数值模拟曲线进行对比，通过非线性迭代反演及模型更新，不断调用步骤三和步骤四正演仿真过程使数值模拟曲线与实测测井曲线一致，最终得到地层条件下包括孔隙度、渗透率、饱和度的储层参数。其中非线性迭代反演是通过求使目标方程达到最小值的一组模型参数，目标方程为：

目标方程受如下公式约束

式中：W_d为数据权重；d_i是模拟得到的某一测井响应；d_m是测量得到的测井响应；α为阻尼因子；n_c为预先设定的矿物组分数量；x为岩石物理模型参数向量；可以表示为：

S_w为地层含水饱和度，

表示地层孔隙度，C_i第i种矿物组分含量，C_sh表示泥质孔隙度，模拟得到的测井响应d可以表示为：

d＝[CNC,DEN,1/R]^T

CNC、DEN、R分别为中子测井，密度测井与电阻率测井值。

与现有技术相比，本发明的有益效果在于：该方法基于岩石物理参数结合井筒信息及相渗特征曲线，从仪器原理出发综合渗流特征、电磁场及核测井原理，通过多时刻、多物理场多维度测井响应信息反演储层实际储层参数，提高了储层流体性质识别及参数评价精度，为油气勘探开发提供技术支持。

附图说明

图1是本发明的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法的流程示意图。

图2是本发明的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法的工作流程示意图。

图3是本发明的泥浆侵入动态剖面图。

图4是本发明的电阻率及核测井响应图。

图5是本发明的实际井基于数值模拟仿真处理结果示意图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部品会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本专利的限制。

本发明实施例的附图中相同或相似的标号对应相同或相似的部品；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”“长”“短”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本专利的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

下面通过具体实施例，并结合附图，对本发明的技术方案作进一步的具体描述：

实施例1

如图1-图5所示为一种泥浆侵入下电阻率及核测井联合反演储层参数的方法的实施例1，

步骤一：构建井筒数值模型；根据储层伽玛、电阻率、中子、密度曲线特征对地下储层实测测井曲线模型进行分层，然后按层进行所述地下储层实测测井曲线方波化，每一层代表一个岩石物理模型；

步骤二：井筒数值模型初始值设定；将岩石物理模型参数、渗流信息和井眼信息参数输入所述步骤一中的所述岩石物理模型，其中岩石物理模型参数和渗流信息参数是按层赋值，每一层按地层储层特征将参数输入所述步骤一中的所述的岩石物理模型，井眼信息参数按照全井段处理，即输入统一的数值；

步骤三：动态数值模拟；根据相关多孔介质渗流理论和扩散理论对所述步骤二获得的所述井筒初始数值模型进行钻井液动态侵入模拟，获得任意时刻的储层动态流体分布剖面，包括储层动态饱和度和矿化度剖面，通过相关公式对所获得的储层动态流体分布剖面进行转换，获得相应的电阻率、迁移长度、密度在不同时刻的径向分布剖面；

步骤四：电阻率、中子、密度测量响应数值模拟；根据电磁场有限元仿真方法结合仪器结构参数及工作原理获得数值模拟的电阻率测井响应，根据仪器结构参数及核测井快速模拟方法获得数值模拟的中子、密度测量曲线；

步骤五：反演储层参数；用所述步骤一中的所述地下储层实测测井曲线与所述步骤四获得的数值模拟曲线进行对比，通过模型更新及迭代，使所述数值模拟曲线与所述实测测井曲线一致，最终得到反演的岩石物理结构的储层参数。

本实施例的有益效果：本方案的处理方法能够开展电阻率及核测井联合动态模拟及反演，得到储层准确的流体性质与孔隙结构参数。

实施例2

一种泥浆侵入下电阻率及核测井联合反演储层参数的方法的实施例2，在实施例1的基础上，如图1-图5所示，对步骤一至步骤四进一步限定。

具体的，在步骤一中，测井曲线是对地层物理性质随井深变化的记录，以自然伽玛半幅点，参考电阻率、中子、密度曲线划分各小层顶、底界面。方波化是是指在图形显示方式下对测井曲线的数值分段进行平均，形成锯齿形态，采用每一小层的测井曲线平均值代表该小层。

具体的，在步骤二中，输入模型参数包括岩石物理模型参数、渗流信息和井眼信息参数，岩石物理模型参数和渗流参数是按层处理，井眼信息参数按照全井段处理，岩石物理模型具体参数包括岩石骨架参数、地层物性参数、地层孔隙流体参数、泥质含量。岩石骨架参数包括矿物组分、含量和密度；地层物性参数包括饱和度、渗透率、孔隙度和孔隙结构；地层孔隙流体参数包括流体组分、含量、密度、矿化度和粘度；渗流信息包括相渗曲线和毛管压力曲线；井眼信息参数包括井眼直径、泥浆信息、钻柱压力和温度。岩石骨架参数通过X衍射资料获取，包括流体组分、密度、矿化度、粘度的孔隙流体信息通过水分析实验材料获得，渗流信息通过地区驱替实验获得，井眼信息通过钻井日志、地质日报获得，地层孔隙结构参数通过电阻率岩心实验获得，泥质含量、孔隙度、渗透率、饱和度通过每一小层的测井曲线平均值采用如下的测井解释公式获得。泥质含量初始值获取公式为

式中,I_sh表示泥质含量指数,V_sh表示地层泥质含量；小数；GR表示地层自然伽马测量值，gAPI；GR_max表示本井纯泥岩地层自然伽马最大值，gAPI；GR_min表示本井纯砂岩地层自然伽马最小值，gAPI。

孔隙度初始值获取公式为

式中：

表示计算的密度孔隙度，v/v；DEN表示密度曲线测井读值，g/cm3；DEN_ma表示骨架的密度响应参数，g/cm3；DEN_f表示流体的密度响应参数，g/cm3；V_sh表示地层泥质含量，v/v；DEN_sh表示泥岩的密度响应参数，g/cm3；

表示计算的中子孔隙度，v/v；CNC表示中子曲线测井读值，％；CNC_ma表示骨架的中子响应参数，％；CNC_f表示流体的中子响应参数，％；CNC_sh表示泥岩的中子响应参数，％；

表示地层孔隙度，v/v。

渗透率计算可采用区域岩心孔渗回归关系模型，例如，研究区模型公式为：

式中：K表示地层渗透率，mD；

表示地层孔隙度，v/v。

含水饱和度初始值获取公式为：

式中：S_w表示地层含水饱和度；R_t表示地层真电阻率，Ω.m；V_sh表示泥质含量，小数；R_sh表示纯泥岩电阻率，Ω.m；R_w表示地层水电阻率，Ω.m；

表示孔隙度，小数；a表示岩性系数；m表示胶结指数；n表示饱和度指数，a，m，n采用区块经验值或者其他测井解释方法获得。

具体的，在步骤三中，相渗曲线及毛管压力曲线通过岩心实验获得。根据包括相渗曲线及毛管压力曲线的相关渗透率理论和离子扩散理论模拟获得任意时刻包括地层饱和度及矿化度的储层流体动态分布剖面，通过电阻率及核测井相关公式转换获得相应的电阻率、迁移长度、密度在不同时刻的径向分布剖面。电阻率径向剖面转化公式为：

式中：R_t表示地层电阻率；R_w表示地层水电阻率；S_w表示地层含水饱和度；

表示地层孔隙度；R_sh表示泥质电阻率；V_sh表示泥质含量；a表示岩性系数；m表示胶结指数；n表示饱和度指数；c表示泥质含量指数系数，一般计算公式为

a，m，n采用区块经验值或者其他测井解释方法获得。

密度径向剖面转换公式为：

式中：ρ_b表示地层密度；ρ_w表示地层水密度；ρ_h表示油气流体密度；

表示地层孔隙度；ρ_sh表示泥质密度；C_sh表示泥质孔隙度；S_w表示地层含水饱和度；ρ_i骨架中第i种组分的密度；C_i表示骨架中第i种组分的含量；n_c表示骨架中组分种类总数。

中子迁移长度径向剖面计算公式为：

式中：ξ表示地层迁移长度；S_w表示地层含水饱和度；ξ_w表示地层水迁移长度；ξ_h表示油气流体迁移长度；

表示地层孔隙度；ξ_m表示骨架迁移长度；η表示指数系数。

具体的，在步骤四中，根据步骤三转换的电阻率径向分布剖面获取电阻率测井响应，主要采用电磁场有限元仿真方法结合仪器结构参数及工作原理数值模拟获得对应的电阻率测井响应。其中采用三维矢量有限元素法，对测井物理模型进行空间离散，建立单元电磁场离散方程，对全部单元安装、消元、求解，进而得到电测井仪器响应。电阻率测井仿真为求解给定边界条件下的麦克斯韦方程问题。从麦克斯韦方程出发，电磁场满足下面方程：

式中：E表示电场强度；H表示磁场强度；J为源电流密度；ω为角频率；σ为电导率；μ为磁导率。

从上述两个方程出发，可以推导出电场所满足的矢量波动方程为：

为复介电常数，ε＝ε_rε₀其中ε₀为真空介电常数，ε_r为相对介电常数。

令

E＝E_p+E_s

式中，背景场E_p是当全部空间被电导率为σ₀的介质填充时的电场，它满足方程：

其中，

由上述三式变化可得到矢量波动方程：

背景场通过解析方法计算得到，二次场则由有限元素法计算。上式的解变化平缓，可以利用稀疏一些的网格来进行求解，减少了计算工作量。选取足够大区域，使边界上的电场衰减到近似为0，则上式只需满足边界条件方程：

式中：

为求解区ω的边界，n为其法线方向。

考虑边界条件方程，将矢量波动方程转化为其弱积形式：

式中：N为矢量基函数，V为单个求解元素，Ω为整个求解空间，E_p表示背景场电场，E_s表示二次场电场,ε_c0表示背景场的介电常数，ε表示二次场介电常数，ω为角频率，μ为磁导率。

具体的，在步骤四中，根据步骤三转换的地层密度径向分布剖面和迁移长度径向分布剖面，结合仪器结构参数及核测井快速模拟方法数值模拟获得中子和密度测井响应。

其基本原理是在一定体积内中子、伽马光子随时间变化率等于产生率减去泄漏率和吸收率，可以用下式来描述这一过程。探测器处通量可以表示为每个从源发射的粒子对探测器贡献的总和。

对于密度测井，伽马源发射伽马光子，与原子核发生康普顿散射、光电散射并产生散射伽马射线，探测器主要探测散射伽马射线的计数率。对于中子孔隙度测井，脉冲中子源发射快中子，与原子核发生非弹性散射、弹性散射反应，中子能量降低变为热中子，探测器主要探测热中子的通量。根据上面理论，背景区域通量可以表示为：

N_B(r_R)＝∫dr∫dE∫dΩψ_B(r_s,r,E,Ω)S(r_R,r,E,Ω)

式中N_B(r_R)表示背景区域探测器的通量，n；r_s表示空间中源的位置，无量纲；r_R表示空间中探测器的位置，无量纲；r表示空间中某一点与源的距离，cm；Ω表示粒子的碰撞角度，°；E表示粒子的能量，MeV；ψ_B表示位于r_s处的放射源发射的粒子在r处的通量，n/(cm²·s)。S表示r_R处的探测器响应函数，可以表示为：

S(r_R,r,E,Ω)＝∫dE′∫dΩ′∑_s(r,E,Ω→E′Ω′)Ψ⁺(r_R,r,E′,Ω′)

式中∑_S表示粒子的宏观散射截面，1/cm；Ω′表示粒子碰撞后角度，°；E′表示粒子碰撞后能量，MeV。Ψ⁺为重要性方程。定义背景通量灵敏度方程为：

变化区域对探测器通量的贡献可以表示为：

式中Δ∑为变化区域粒子截面，1/cm；∑_B为背景区域粒子截面，1/cm。探测器通量为背景通量与变化区域通量的叠加，表示为：

对于密度测井，变量设置为密度值如下式所示，C_p为密度系数，无量纲；ρ表示密度，g/cm³；ρ_B表示背景密度，g/cm³。

对于中子测井，变量设置为迁移长度导数，如下式所示，C_N为中子系数；Lm为迁移长度，无量纲；LmB为背景迁移长度，无量纲。

本实施例的有益效果：在步骤一中，对测井曲线方波化，对测井曲线的数值分段进行平均分层，可以使每小层的测井曲线更具有代表性。在步骤四中，采用三维矢量有限元素法，对测井物理模型进行空间离散，建立单元电磁场离散方程，对全部单元安装、消元、求解，可以准确的得到电测井仪器响应。

实施例3

一种泥浆侵入下电阻率及核测井联合反演储层参数的方法的实施例3，在实施例1或2的基础上，如图1-图5所示，对步骤五进一步限定。

具体的，在步骤五中，用地下储层测井曲线与数值模拟曲线进行对比，通过非线性迭代反演及模型更新，不断调用步骤三和步骤四正演仿真过程使数值模拟曲线与实测测井曲线一致，最终得到地层条件下包括孔隙度、渗透率、饱和度的储层参数。其中非线性迭代反演是通过求使目标方程达到最小值的一组模型参数，目标方程为：

目标方程受如下公式约束

式中：W_d为数据权重；d_i是模拟得到的某一测井响应；d_m是测量得到的测井响应；α为阻尼因子；n_c为预先设定的矿物组分数量；x为岩石物理模型参数向量；可以表示为：

S_w为地层含水饱和度，

表示地层孔隙度，C_i第i种矿物组分含量，C_sh表示泥质孔隙度，模拟得到的测井响应d可以表示为：

d＝[CNC,DEN,1/R]^T

CNC、DEN、R分别为中子测井，密度测井与电阻率测井值。

本实施例的有益效果：在步骤五中，通过非线性迭代反演及模型更新，可以快速的调用步骤三和步骤四正演仿真过程使数值模拟曲线与实测测井曲线一致，最终得到地层条件下包括孔隙度、渗透率、饱和度的储层参数。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (10)

1.一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，包括如下步骤：

步骤一：构建井筒数值模型；根据储层伽玛、电阻率、中子、密度曲线特征对地下储层实测测井曲线模型进行分层，然后按层进行所述地下储层实测测井曲线方波化，每一层代表一个岩石物理模型；

步骤二：井筒数值模型初始值设定；将岩石物理模型参数、渗流信息和井眼信息参数输入所述步骤一中的所述岩石物理模型，其中岩石物理模型参数和渗流信息参数是按层赋值，每一层按地层储层特征将参数输入所述步骤一中的所述的岩石物理模型，井眼信息参数按照全井段处理，即输入统一的数值；

步骤三：动态数值模拟；根据相关多孔介质渗流理论和扩散理论对所述步骤二获得的所述井筒初始数值模型进行钻井液动态侵入模拟，获得任意时刻的储层动态流体分布剖面，包括储层动态饱和度和矿化度剖面，通过相关公式对所获得的储层动态流体分布剖面进行转换，获得相应的电阻率、迁移长度、密度在不同时刻的径向分布剖面；

步骤四：电阻率、中子、密度测量响应数值模拟；根据电磁场有限元仿真方法结合仪器结构参数及工作原理获得数值模拟的电阻率测井响应，根据仪器结构参数及核测井快速模拟方法获得数值模拟的中子、密度测量曲线；

步骤五：反演储层参数；用所述步骤一中的所述地下储层实测测井曲线与所述步骤四获得的数值模拟曲线进行对比，通过模型更新及迭代，使所述数值模拟曲线与所述实测测井曲线一致，最终得到反演的岩石物理结构的储层参数。

2.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤二中，所述岩石物理模型参数包括岩石骨架参数、地层物性参数、地层孔隙流体参数、泥质含量；所述渗流信息包括相渗曲线和毛管压力曲线；所述井眼信息参数包括井眼直径、泥浆信息、钻柱压力和温度。

3.根据权利要求2所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，所述岩石骨架参数包括矿物组分、含量和密度；所述地层物性参数包括饱和度、渗透率、孔隙度和孔隙结构；所述地层孔隙流体参数包括流体组分、含量、密度、矿化度和粘度，所述泥浆信息包括钻井液矿化度、泥浆固含、钻井液粘度。

4.根据权利要求3所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤二中初始值设定，所述岩石骨架参数通过X衍射资料获取，包括流体组分、密度、矿化度、粘度等，孔隙流体信息通过水分析实验材料获得，所述渗流信息通过岩心驱替实验获得，所述井眼信息通过钻井日志、地质日报获得，所述地层孔隙结构参数通过电阻率岩心实验获得，所述泥质含量、所述孔隙度、所述渗透率、所述饱和度通过每一小层的所述测井曲线平均值采用常规测井解释公式获得。

5.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤三中，电阻率径向分布剖面基于饱和度和矿化度动态剖面采用印尼公式转换获得，密度径向分布剖面基于饱和度动态剖面、孔隙度、孔隙流体组分和岩石骨架组分密度值，通过岩石物理体积模型计算获得。

6.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤三中，迁移长度径向分布剖面基于饱和度动态剖面、孔隙度、孔隙流体和岩石骨架组分迁移长度值，通过岩石物理体积模型计算获得。

7.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤四中，所述仪器为电阻率测井仪、中子测井仪、密度测井仪。

8.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤四中，根据所述步骤三转换的地层电阻率径向分布剖面，采用三维矢量有限元素法，对测井物理模型进行空间离散，建立单元电磁场离散方程，对全部单元安装、消元、求解，进而得到所述电阻率测井仪响应曲线。

9.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤四中，根据所述步骤三转换的地层密度径向分布剖面和迁移长度径向分布剖面，再结合仪器结构参数及核测井快速模拟方法得到所述密度、中子响应曲线。

10.根据权利要求1所述的一种泥浆侵入下电阻率及核测井联合反演储层参数的方法，其特征在于，在所述步骤五中，用地下储层测井曲线与数值模拟曲线进行对比，通过非线性迭代反演及模型更新，不断调用所述步骤二和所述步骤四正演仿真过程使数值模拟曲线与实测测井曲线一致，最终得到地层条件下的储层参数。

Images