CN114856537A - 一种脉冲中子氧活化测井水流速度最优化计算方法 - Google Patents

Abstract

本发明公开了一种脉冲中子氧活化测井水流速度最优化计算方法，主要实现了根据采集的活化时间谱计算井下水流速度。本发明采用蒙特卡洛数值模拟技术证明了可以采用洛伦兹分布来表征脉冲中子氧活化测井的氧活化反应率分布函数和探测器响应函数。据此获得了活化时间谱中各时刻计数率的理论值。结合实测活化时间谱各时刻计数率，建立求解水流速度的最优化模型。通过寻优获得使各时刻理论计数率与实测计数率间误差最小的水流速度，该水流速度与实际水流速度最为接近。本发明从原理上描述了水流速度与活化时间谱计数率之间的关系，可以直接从活化时间谱中提取水流速度，避免了估计活化水流到达探测器时间，且能够计算复杂管柱下的多水流速度。

Description

一种脉冲中子氧活化测井水流速度最优化计算方法

技术领域

本发明涉及一种脉冲中子氧活化测井水流速度最优化计算方法，主要用于根据脉冲中子氧活化测井采集的活化时间谱计算井下水流速度。同时也解决了复杂管柱结构下存在多水流时的速度计算问题。

背景技术

注水开发是实现油田中后期稳产、增产的关键措施之一，准确的监测注水剖面对评价注水效果具有重要意义。在注入剖面监测中最常用的方法是采用涡轮流量计，但涡轮流量计在较长的油管中测量有遇阻遇卡的风险，且涡轮流量计对启动排量有一定要求，不适用于对低流量的测量。井温测井，放射性示踪测井以及频谱噪声测井也是近年来用于评价注入剖面的生产测井方法。然而，井温测井难以获得高质量的定量结果；当注水井管柱较为粗糙或注入水质欠佳时，放射性示踪测井必须进行沾污校正；频谱噪声测井虽然可以直接监测近井区的水流流动，但却很难区分油管和环空中的水流。除此之外，像分布式温度传感器(DTS)和分布式声波传感器(DAS)这类分布式光纤传感器技术目前也已用于注入剖面的监测中。但DTS需要在理想的测量环境下才能进行较好的定量分析，而DAS的定量解释方法目前仍在研究。相比而言，现阶段脉冲中子氧活化测井技术在油田注入剖面的监测中具有明显的优势。脉冲中子氧活化测井仪器在井中进行定点测量，仪器D-T中子源发射的14MeV的快中子能够穿透套管和油管与井内水流中的氧原子核发生活化核反应，因此脉冲中子氧活化测井不会受管柱结构的影响。此外，脉冲中子氧活化测井还能够同时监测油管和油套环空中的水流。在复杂管柱结构下，它还能区分不同测量井段的水流方向。

当D-T中子源发射的快中子与井内水流中的氧原子核发生活化核反应时，被活化的氧原子核会释放出活化伽马射线，当活化水流靠近并通过仪器伽马探测器时，其释放出的活化伽马射线会被伽马探测器采集，并在记录的活化时间谱上形成一个明显的信号峰。现有方法根据计时原理，通过估算活化水流到达伽马探测器的时间，结合已知源距即可计算水流速度。

为了准确估算活化水流到达伽马探测器的时间，一些研究人员提出活化时间谱中信号峰的峰位所对应的时间即为活化水流到达伽马探测器的时间，在解释时采用加权平均法对峰位进行近似计算。一些学者也提出了采用函数拟合法求取峰位时间，比如采用高斯函数加线性函数残式最小二乘法拟合确定峰值对应时间；或利用高斯函数拟合信号峰，同时采用双高斯函数来对多水流下重叠的信号峰进行分离；以及将活化时间谱的形态分为对称型和非对称型两种，分别采用高斯函数与泊松函数来进行拟合计算等方法。然而由于活化水流在流动时水中的活化氧原子核会发生衰变，因此活化时间谱中的信号峰并非是严格对称的。这使得采用加权平均法或高斯函数拟合法的结果存在一定的误差。此外，通过固定的拟合函数来求解活化水流到达探测器的时间只能得到一个估计解，并无法从原理上来解释活化时间谱中信号峰与水流速度之间的关系。

发明内容

本发明主要克服现有技术中的不足之处，目的在于提供一种基于脉冲中子氧活化测井采集的活化时间谱来准确计算水流速度的方法，同时也适用于复杂管柱结构下井中存在多水流时的速度计算。为了达到以上技术目的，本发明采用以下技术方案：

脉冲中子氧活化测井仪器采用D-T中子源向井中发射能量为14MeV的快中子，快中子与井眼以及地层中的氧原子核发生活化核反应产生放射性同位素¹⁶N，¹⁶N以7.13s的半衰期发生β衰变产生处于激发态的¹⁶O原子核，当处于激发态的氧原子发生退激回到基态时会释放出能量为6.13MeV和7.12MeV的活化伽马射线。仪器伽马探测器会记录不同时刻的活化伽马射线计数率，即活化时间谱。因此可以将¹⁶N看作一种短时效的示踪剂。脉冲中子氧活化测井采集的活化时间谱伽马计数率随时间的变化关系可表示为：

式中λ为¹⁶N的衰减系数；C(t)为活化时间谱t时刻计数率；S(t)为t时刻中子源强度；t₀为中子脉冲起始时间；t_a为中子脉冲持续时间；A为流动空间横截面积；z为到探测器的距离；函数D(z)为t时刻的氧活化反应率分布；R(z)为伽马探测器响应函数；L为伽马探测器到源的距离；v为水流速度。该模型描述了由活化水流引起的活化时间谱中不同时刻的计数率。然而，在实际测量中R(z)和D(z)难以确定，且当井内存在多个水流时，活化时间谱会更为复杂，因此直接利用此模型来通过活化时间谱的计数率计算水流速度是不可行的。

为了获得脉冲中子氧活化测井中氧活化反应率分布函数D(z)与仪器响应函数R(z)的表征形式，采用蒙特卡洛数值模拟技术建立两个计算模型对其进行模拟研究。第一个模型如图1所示，用来确定D(z)的分布规律。地层为直径90cm，高120cm的纯砂岩圆柱体，钢制套管直径13.97cm(5.5in)，水泥环厚度为3cm，氧活化测井仪器直径38mm包含一个D-T中子源和两个不同源距的闪烁体伽马探测器，其中长源距80cm，短源距30cm，仪器在井内居中测量。为了更有效的进行模拟，该计算模型移除了油管，而且井眼内充满了清水。为了获得氧活化反应率随与中子源距离的变化关系，井内流体在模型中为静止状态。

采用上述计算模型模拟得到的氧活化反应率分布如图2所示。可以看出井中的流体距离中子源越近则越易被活化并产生更多的活化氧原子核。氧活化反应率分布可以近似采用洛伦兹分布来表示，如式(2)所示，两者的相关系数(R²)可达到0.9962。

式中，h为洛伦兹峰高，p为洛伦兹峰位，w为洛伦兹峰宽。图2所展示的是井内流体静止条件下的氧活化反应率分布。但在实际井中，井内流体是以特定的速度流动的，并且脉冲中子源发射中子也会持续一定的时间，一般为1-15sec。因此在中子脉冲停止后的氧活化反应率分布会和图2所展示的有所不同。实际上根据脉冲中子氧活化测井原理，氧活化反应率分布可以看作是水流速度与脉冲时间的函数，如式(3)所示。因此不同脉冲时间与水流速度下的氧活化反应率分布可以根据该式计算出来。

式中，D′(z)为中子脉冲结束后的氧活化反应率分布函数，z为到中子源的距离。

第二个计算模型用来确定仪器探测器的响应函数，如图3所示。该计算模型的模拟地层和井筒参数与第一个模型一致。为了确定探测器的响应函数，沿着井筒设置多个探测器，并且把活化水流看作是井筒内的一个放射性点源。该计算模型可以等效模拟活化水流沿井筒流动时的探测器响应。通过统计各探测器的计数率以及放射源到各探测器的距离可以得到仪器探测器响应函数，如图4所示。为了方便分析，源强随着时间的衰减暂时先不考虑。显然放射源(活化水)距离探测器越近，则探测器响应强度越大而且探测器的响应函数同样可以采用洛伦兹分布来近似描述，相关系数可达0.9999。

基于以上分析，在考虑¹⁶N随时间的衰减后，脉冲中子氧活化测井测量响应的理论模型可以表示为式(4)。

式中，n为井中的水流数。该式描述了随着活化水流动，仪器伽马探测器的响应规律。可以用此模型来预测活化时间谱中各道的计数率。基于上述结论，即可根据实测活化时间谱计数率，采用最优化方式实现井下水流速度的计算。

一种脉冲中子氧活化测井水流速度最优化计算方法其特征在于所述方法包括以下步骤：

S1：对活化时间谱进行质量控制。质量控制包括异常点修正以及降噪处理两部分，用于降低仪器噪声以及放射性统计涨落对活化时间谱的影响，提高活化时间谱的信噪比，使活化时间谱中的水流特征更为明显，便于后续的处理。

S2：对活化时间谱进行静态活化氧校正。由于氧原子核在地层和井下水泥环中均有存在并且可被活化，这部分静态的活化氧原子核释放的活化伽马射线同样会被探测器采集。消除活化时间谱中静态活化氧产生的计数率，提取出来自活化水流的流动活化氧计数率，同样能够使水流特征更为明显。谱线中的静态活化氧计数可以视为一个指数衰减函数，因此可以通过函数拟合的方式来估算静态活化氧计数率并对其进行校正。静态活化氧通常仅对距离中子源最近的探测器产生影响，而当探测器距离中子源较远时，则受静态活化氧影响较小。

S3：构造目标函数。水流速度最优化计算方法的关键在于建立能够评价水流速度准确性的目标函数。式(4)描述了不同水流数以及不同水流速度下所产生的活化时间谱各时刻理论计数率。根据实测活化时间谱的各时刻计数率与理论计数率之间的误差可以判断计算水流速度与实际水流速度之间的接近程度，当误差最小时计算水流速度与实际水流速度最为接近。

基于此原理可以构造出计算井下水流速度的目标函数，如式(5)所示。

式中，m为中子脉冲结束后的采样点数，M(t_j)为t_j时刻实测活化时间谱计数率，C(t_j)为t_j时刻理论活化时间谱计数率，e为实测活化时间谱计数率与理论计数率的累计误差。

S4：水流速度优化计算。式(5)所表示的目标函数为一个典型的非线性函数，常规方法难以进行求解，因此本发明引入群体智能最优化计算方法来对目标函数进行求解。其核心思想是将若干个简单个体构成一个群体，通过合作、竞争、交互与学习等机制表现出高级和复杂的功能，完成复杂问题的求解。本发明中采用粒子群优化(PSO)来对该目标函数进行求解。需指出的是，本发明并不对求解目标函数所采用的智能计算方法进行限定，任何同类型的计算方法均可用于对目标函数进行求解，在此仅对PSO算法求解该目标函数的方案进行说明。

PSO是受鸟类觅食行为启发而提出的一种多维进化计算方法。PSO首先在可行域内初始化一组粒子，每一个粒子具有速度、位置与适应性三个属性并可表示为可行域内的一个解，其中适应性由式(5)表示的目标函数确定，其值越小则适应性越好。PSO更新种群的方式如式(6)和式(7)所示：

式中，V_i ^k为第k次迭代中第i个粒子的速度向量，

为第k次迭代中第i个粒子的位置向量，P_i ^k为第k次迭代中第i个粒子的个体最优位置向量，

为第k次迭代中粒子群的全局最优位置向量，ω为惯性权重，c₁，c₂为加速因子，r₁，r₂为区间[0，1]上的随机数。

在每次迭代中，PSO通过评价每一个粒子的适用性来更新其个体最优位置以及群体最优位置，直至群体最优的适应性满足要求。对于脉冲中子氧活化测井的水流速度计算这一优化问题，其解向量可以写为：

S＝[h_Di,w_Di,v_i,h_R,w_R] i＝1,...,n (8)

式中，h_Di为第i个水流产生的氧活化反应率分布的洛伦兹峰高，w_Di为第i个水流产生的氧活化反应率分布的洛伦兹峰宽，h_R为仪器响应函数的洛伦兹峰高，w_R为仪器响应函数的洛伦兹峰宽，v_i为第i个水流的速度。

本发明所提出的脉冲中子氧活化测井最优化水流速度计算方法的计算流程如图5所示。异常点处理与活化时间谱降噪是两个必要的预处理步骤，可以降低放射性统计涨落与仪器噪声产生的不利影响。静态活化氧校正对于距离中子源较近的探头采集的活化时间谱是必要的处理步骤，而当探头距离中子源较远，受静态活化氧影响不大时，可依据实际情况选择性处理。之后，根据实测活化时间谱计数率构造出目标函数。采用PSO对目标函数进行寻优处理，当满足理论计数率与实测计数率间误差最小时，即可获得实际水流速度。

和现有方法相比，本发明所提出的最优化水流速度计算方法可以直接从活化时间谱中提取水流速度，同时从原理上描述了水流速度与活化时间谱计数率之间的关系。避免了现有方法对活化水流到达探测器时间的估计，并且能够用于复杂管柱下多水流活化时间谱的水流速度计算。

附图说明

图1为确定脉冲中子氧活化测井氧活化反应率分布函数的蒙特卡洛计算模型示意图。

图2为归一化氧活化反应率与距中子源距离的函数关系图。

图3为确定脉冲中子氧活化测井探测器响应函数的蒙特卡洛计算模型示意图。

图4为归一化探测器响应强度与距探测器距离的函数关系图。

图5为脉冲中子氧活化测井水流速度最优化计算方法流程图。

图6为实例井管柱结构及射孔层位示意图。

图7为脉冲中子氧活化测井仪器结构示意图。

图8为各测点水流速度优化计算结果图。

图9为层流与紊流速度剖面对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合应用实例，对本发明进行进一步详细说明。应当理解，此处所描述的具体应用实例仅用以解释本发明，并不用于限定本发明。

应用实例：

本实例以某一实测注水井中脉冲中子氧活化测井在各测点采集的活化时间谱来说明本发明的应用效果。

该实例井采用油套分注的注水方式进行注水开发，即从油管与油套环空中同时注水，管柱结构如图6所示。该井套管直径139.7mm，油管直径62mm。油管与环空的平均注水量分别为25.86m³/d和33.31m³/d。测量井段为2125m到2220m，其中包含两个射孔层且两个层均在喇叭口以上。层1起止深度为2186m到2194m，层2起止深度为2206m到2210m。从油管中注入的水会从该井的底部喇叭口进入油套环空向上流动。显然，在油管中仅存在向下流动的水流，而在油套环空内上水流和下水流均存在。

为了评价两个解释层的注入情况，利用脉冲中子氧活化测井仪对不同深度的8个测点进行定点测量，分别是2180m，2190m，2193m，2202m，2210m和2218m。解释层位处的流量变化可以通过位于解释层两端的测点所测量的流速以及对应的流动截面积来确定。该井所使用的测量仪器直径为38mm，测井时置于油管中。仪器包含8个伽马探测器分别对称的位于D-T中子源的两侧。自上而下可以依次标记为U4，U3，U2，U1，D1，D2，D3，D4，对应的源距为240cm，150cm，80cm和30cm，仪器结构如图7所示。这种探测器分布可以在一次测量中实现油套管中上水流和下水流的同时测量并且提高了仪器对流速的测量范围。

仪器在不同测点采集到的氧活化时间谱以及采用本发明方法优化得到理论谱如图8所示。可以看出油管中的水流信号峰和环空中的水流信号峰相比峰值更高，峰宽更窄，并且油管中的活化水流信号峰具有更好的信噪比，这是因为油管中的水流更靠近中子源的原因。

采用本发明计算的各测点水流速度同样显示在图8中。本发明不需要计算活化水流从中子源流动到探测器的时间，这也降低了放射性统计涨落以及¹⁶N同位素衰减对结果的影响，并使多水流信号的分离成为可能。根据地面统计的注水量以及管柱截面积可以计算出油管和环空顶部和底部的实际水流速度。油管中的水流平均速度为15.88cm/s，2180m和2210m处测点的环空水流平均速度分别为3.13cm/s和2.43cm/s。本发明根据该位置处的测点采集的活化时间谱计算得到的水流速度以及与实际水流速度间的对比结果如表1所示。

可以看出，油管中的水流速度计算结果与实际结果十分接近，相对误差最低能达到0.76％。而油套环空内的水流速度计算结果与实际水流速度相比误差较大。这主要是由于井下流体流动速度分布的不均性所导致的。

表1计算流速与实际流速对比结果

流体的流动类型可以根据雷诺数(Re)来划分。Re是一个无因次数，可以用来表征流体的流动情况。它可以根据式(9)计算：

式中，ρ为流体密度；v为流体速度；μ为流体粘度；D为等效直径。

流体的流动类型分为层流和紊流两类，其流动速度剖面如图9所示。大量实验证明当Re<2000时，流动类型为层流；而当Re>4000时，流动类型为紊流。表2展示了上述已知平均流速的测点的流动类型以及雷诺数。可以看出环空中的流体流动类型为层流，对于层流来说，中心处的流体速度要大于周边的流动速度，因此中间处流动速度较快的流体会对最终的计算结果造成影响。而对于油管中的紊流来说，其速度剖面更为平坦。这意味着速度剖面上各处的流体速度均接近于平均流速，这也使得油管速度的计算结果更为准确。

表2已知流速测点的雷诺数与流动类型

在实际脉冲中子氧活化测井中，流体流动速度分布的不均匀性无法避免。因此，最常用的处理方法是利用已知流速的测点来对所有计算结果进行刻度，计算对应的速度剖面校正系数，如式(10)所示。

式中，C_v为速度剖面校正系数，v_m为已知测点平均流速，v_a为已知测点计算流速。将其他测点的计算结果乘以速度剖面校正系数C_v即可得到最终的计算结果。

表3展示了综合流动空间截面积和本发明计算得到的流体流动速度，最终获得的各射孔层精细注入剖面定量解释结果，包含层内各小层的吸水量与吸水指数。

表3吸水剖面定量解释结果

Claims (6)

1.一种脉冲中子氧活化测井水流速度最优化计算方法其特征在于所述方法包括以下步骤：

S1：对活化时间谱进行质量控制，质量控制包括异常点修正以及降噪处理两部分，用于降低仪器噪声以及放射性统计涨落对活化时间谱的影响，提高活化时间谱的信噪比，使活化时间谱中的水流特征更为明显，便于后续的处理；

S2：对活化时间谱进行静态氧校正，静态活化氧原子核主要来自地层和井下水泥环中的氧元素，这部分静态的活化氧原子核释放的活化伽马射线会影响距离中子源最近的探测器，静态活化氧计数可视为一个指数衰减函数，可以通过函数拟合的方式来估算静态氧计数并对其进行校正；

S3：构造目标函数，根据井内流体静止时的氧活化反应率分布函数D(z)以及仪器探测器响应函数R(z)，建立计算活化时间谱各时刻理论计数率的模型，结合实测活化时间谱的各时刻计数率，根据活化时间谱理论计数率与实测计数率之间的累计误差构造目标函数；

S4：水流速度优化计算，采用群体智能最优化计算方法来对目标函数进行求解，通过寻优处理来最小化目标函数，求解目标函数解向量，从而获得使各时刻理论计数率与实测计数率间误差最小的水流速度，该水流速度为与实际水流速度最接近的最优化计算结果。

2.如权利要求1中所述的一种脉冲中子氧活化测井水流速度最优化计算方法其特征在于所述步骤3中的井内流体静止时的氧活化反应率分布函数D(z)与探测器响应函数R(z)通过蒙特卡洛数值模拟技术获得，均可采用洛伦兹分布来表征：

式中，h为洛伦兹峰高，p为洛伦兹峰位，w为洛伦兹峰宽。

3.如权利要求1中所述的一种脉冲中子氧活化测井水流速度最优化计算方法其特征在于所述步骤3中计算活化时间谱各时刻理论计数率的模型为：

式中，n为井中的水流数，λ为¹⁶N的衰减系数，D′(z)为中子脉冲结束后的氧活化反应率分布函数，R(z)为伽马探测器响应函数，其中自变量z为到探测器的距离，L为伽马探测器到中子源的距离，v_i为第i个水流速度，C(t)为活化时间谱t时刻计数率。

4.如权利要求3中所述的中子脉冲结束后的氧活化反应率分布函数D′(z)可以根据井内流体静止时的氧活化反应率分布函数D(z)来计算，计算公式为：

式中，A为流动空间截面积，z为到中子源的距离，t₀为中子脉冲起始时间，t_a为中子脉冲持续时间。

5.如权利要求1中所述的一种脉冲中子氧活化测井水流速度最优化计算方法其特征在于所述步骤3根据活化时间谱理论计数率与实测计数率之间的累计误差构造的目标函数为：

式中，m为中子脉冲结束后的采样点数，M(t_j)为t_j时刻实测活化时间谱计数率，C(t_j)为t_j时刻理论活化时间谱计数率，e为实测活化时间谱计数率与理论计数率的累计误差。

6.如权利要求1中所述的一种脉冲中子氧活化测井水流速度最优化计算方法其特征在于所述步骤4采用群体智能最优化计算方法来对目标函数进行求解的解向量表示为：

S＝[h_Di,w_Di,v_i,h_R,w_R]i＝1,...,n (5)

式中，h_Di为第i个水流产生的氧活化反应率分布的洛伦兹峰高，w_Di为第i个水流产生的氧活化反应率分布的洛伦兹峰宽，h_R为仪器响应函数的洛伦兹峰高，w_R为仪器响应函数的洛伦兹峰宽，v_i为第i个水流的速度。

Images