CN115095318A - 一种基于粒子滤波的井眼轨迹预测方法及装置 - Google Patents

一种基于粒子滤波的井眼轨迹预测方法及装置 Download PDF

Info

Publication number
CN115095318A
CN115095318A CN202210723052.6A CN202210723052A CN115095318A CN 115095318 A CN115095318 A CN 115095318A CN 202210723052 A CN202210723052 A CN 202210723052A CN 115095318 A CN115095318 A CN 115095318A
Authority
CN
China
Prior art keywords
well
predicted
angle
section
well section
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210723052.6A
Other languages
English (en)
Other versions
CN115095318B (zh
Inventor
曾义金
张洪宝
陶新港
李梦刚
杨顺辉
柏侃侃
周非
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
China University of Petroleum Beijing
Sinopec Research Institute of Petroleum Engineering
Original Assignee
China Petroleum and Chemical Corp
China University of Petroleum Beijing
Sinopec Research Institute of Petroleum Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, China University of Petroleum Beijing, Sinopec Research Institute of Petroleum Engineering filed Critical China Petroleum and Chemical Corp
Priority to CN202210723052.6A priority Critical patent/CN115095318B/zh
Publication of CN115095318A publication Critical patent/CN115095318A/zh
Application granted granted Critical
Publication of CN115095318B publication Critical patent/CN115095318B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B47/00—Survey of boreholes or wells
    • E21B47/02—Determining slope or direction
    • E21B47/022—Determining slope or direction of the borehole, e.g. using geomagnetism
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mining & Mineral Resources (AREA)
  • Geology (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computational Linguistics (AREA)
  • Software Systems (AREA)
  • Biomedical Technology (AREA)
  • Geophysics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Fluid Mechanics (AREA)
  • Artificial Intelligence (AREA)
  • Health & Medical Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geochemistry & Mineralogy (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本说明书提供了一种基于粒子滤波的井眼轨迹预测方法及装置。该方法包括：获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。基于上述方法能够解决现有方法中存在的井眼轨迹预测算法自适应能力差的技术问题，实现井眼轨迹的精确预测和对井眼轨迹预测结果的不确定描述。

Description

一种基于粒子滤波的井眼轨迹预测方法及装置
技术领域
本说明书属于石油天然气勘探开发技术领域，尤其涉及一种基于粒子滤波的井眼轨迹预测方法及装置。
背景技术
在油气开发领域，钻井过程中井眼轨迹精度会影响油井产量。通常情况下，井眼轨迹难以随钻高频实时测量，现场施工人员主要依靠个人经验和间隔测量信息对定向施工参数进行调整，依据调整后的定向施工参数确定井眼轨迹。定向施工参数调整决策不合理，将会导致井眼轨迹预测不准确，进而导致井眼质量和施工效率降低，安全风险增大。
现有的井眼轨迹预测方法包括：几何外推法、力学分析法和统计方法。其中，几何外推法没有考虑钻井参数和地层特性对井眼轨迹的影响，计算准确性不足；力学分析法采用的模型中考虑因素较多，输入参数不易获取，难以在钻井现场广泛应用；统计方法利用邻井或已钻井眼数据建立回归模型，通过回归系数反映待钻井眼的地层特性，但是其自适应能力差。
因此，目前亟需一种解决自适应能力差问题的井眼轨迹预测方法。
发明内容
本说明书提供了一种基于粒子滤波的井眼轨迹预测方法及装置，能够解决现有方法中存在的井眼轨迹预测算法自适应能力差的技术问题，得到待预测井段的井眼轨迹参数和井眼轨迹参数置信区间。
本说明书实施例的目的是提供一种基于粒子滤波的井眼轨迹预测方法，包括：
获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；
根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；
基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
进一步地，所述方法的另一个实施例中，所述获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据，包括：
基于轨迹测量数据对井眼进行第一次分段，得到井眼第一分段结果；
根据所述井眼第一分段结果和录井米数据，得到钻进数据；
根据所述井眼第一分段结果和定向井施工记录数据，得到工具面角平均值；
基于邻井参数数据对井眼进行第二次分段，得到井眼第二分段结果。
进一步地，所述方法的另一个实施例中，所述根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程，包括：
按照以下算式构建井斜角观测方程：
Incb，i＝(θ0，i×cos(tfb，i)+θ1，i×WOBb，i+θ2，i)×Lb，i+θ3，i×Lh，i+θ4，i+Incb，i-1
式中，i和i-1表示待预测的井段的编号，b表示底部，Incb，i为待预测井段底部井斜角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Incb，i-1为待预测井段顶部井斜角，θ0，i、θ1，i、θ2，i、θ3，i、θ4，i为i井段的井斜角状态参数，其中，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，ii井段的为井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差。
进一步地，所述方法的另一个实施例中，所述根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程，包括：
按照以下算式构建井斜角状态方程：
θ0，i＝θ0，i-1+w0，i-1
θ1，i＝θ1，i-1+w1，i-1
θ2，i＝θ2，i-1+w2，i-1
θ3，i＝θ3，i-1+w3，i-1
θ4，i＝θ4，i-1+w4，i-1
式中，w0，i-1、w1，i-1、w2，i-1、w3，i-1、w4，i-1为i-1井段的井斜角高斯白噪声，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，i为i井段的井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差，θ0，i-1为i-1井段的井斜角造斜率修正系数，θ1，i-1为i-1井段的井斜角钻压影响系数，θ2，i-1为i-1井段的井斜角造斜率误差，θ3，i-1为i-1井段的井斜角复合钻进距离斜率变化规律系数，θ4，i-1为i-1井段的井斜角***误差。
进一步地，所述方法的另一个实施例中，所述根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程，包括：
按照以下算式构建方位角观测方程：
Azib，i＝(φ0，i×cos(tfb，i)+φ1，i×WOBb，i+φ2，i)×Lb，i+φ3，i×Lh，i+φ4，i+Azib，i-1
式中，i和i-1表示待预测的井段的编号，b表示底部，Azib，i为待预测井段底部方位角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Azib，i-1为待预测井段顶部方位角，φ0，i、φ1，i、φ2，i、φ3，i、φ1，i为i井段的方位角状态参数，其中，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差。
进一步地，所述方法的另一个实施例中，所述根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程，包括：
按照以下算式构建方位角状态方程：
φ0，i＝φ0，i-1+n0，i-1
φ1，i＝φ1，i-1+n1，i-1
φ2，i＝φ2，i-1+n2，i-1
φ3，i＝φ3，i-1+n3，i-1
φ4，i＝φ4，i-1+n4，i-1
式中，n0，i-1、n1，i-1、n2，i-1、n3，i-1、n4，i-1为i-1井段的方位角高斯白噪声，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差，φ0，i-1为i-1井段的方位角造斜率修正系数，φ1，i-1为i-1井段的方位角钻压影响系数，φ2，i-1为i-1井段的方位角造斜率误差，φ3，i-1为i-1井段的方位角复合钻进距离斜率变化规律系数，φ4，i-1为i-1井段的方位角***误差。
进一步地，所述方法的另一个实施例中，所述基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间，包括：
按照以下方式，通过粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，确定目标井中当前待预测井段的井眼轨迹参数和井眼轨迹参数置信区间：
根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数；
根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角；
根据当前待预测井段的权重归一化后的粒子群，得到当前待预测井段的重采样粒子群；
基于当前待预测井段的重采样粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，获得当前待预测井段的底部最终井斜角、当前待预测井段的底部最终方位角、当前待预测井段底部井斜角预设的预测置信区间、当前待预测井段底部方位角预设的预测置信区间。
进一步地，所述方法的另一个实施例中，所述根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数，包括：
判断当前待预测井段是否为第一个井段；
在当前待预测井段是第一个井段的情况下，依据预处理数据获得待预测目标井的先验高斯分布，并基于先验高斯分布生成初始粒子群；依据初始粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数；
在当前待预测井段不是第一个井段的情况下，获取前一个井段重采样后的粒子群，依据前一个井段重采样后的粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数。
进一步地，所述方法的另一个实施例中，所述根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，包括：
依据当前待预测井段的状态参数，求解井斜角观测方程和方位角观测方程，得到当前待预测井段的底部井斜角和当前待预测井段的底部方位角；
依据当前待预测井段的底部井斜角和当前待预测井段的底部方位角，得到当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差；
根据当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差重新计算当前待预测井段的粒子群中每一个粒子的权重，获得当前待预测井段的重新赋权的粒子群；
依据当前待预测井段的重新赋权的粒子群，获得当前待预测井段的权重归一化后的粒子群。
另一方面，本申请提供了一种基于粒子滤波的井眼轨迹预测装置，包括：
预处理模块，用于获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；
构建模块，用于根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；
粒子滤波模块，用于基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
再一方面，本申请还提供了一种计算机可读存储介质，其上存储有计算机指令，所述计算机可读存储介质执行所述指令时实现上述基于粒子滤波的井眼轨迹预测方法。
本说明书提供的一种基于粒子滤波的井眼轨迹预测方法及装置，通过获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。依据本说明书提供的方法，能够解决现有方法中存在的井眼轨迹预测算法自适应性差的技术问题，实现井眼轨迹预测和对井眼轨迹进行不确定性描述的技术效果。
并且，在基于粒子滤波，确定目标井的井眼轨迹参数和井眼轨迹参数置信区间时，通过待预测的前一个井段的粒子群计算得到待预测井段的状态参数，根据待预测井段的状态参数计算得到待预测井段底部井斜角、待预测井段底部方位角，对粒子群进行归一化、重采样操作后，依据重采样后的粒子群获得待预测井段底部最终井斜角、待预测井段底部最终方位角，以对井眼轨迹进行预测，实现在线自适应校准井眼轨迹预测参数的目的，从而对定向施工参数、随钻数据进行调整。
附图说明
为了更清楚地说明本说明书实施例，下面将对实施例中所需要使用的附图作简单地介绍，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本说明书提供的一种基于粒子滤波的井眼轨迹预测方法一个实施例的流程示意图；
图2是本说明书一个实施例中的预处理数据图；
图3是本说明书一个实施例中的井斜角先验高斯分布求取结果图；
图4是本说明书一个实施例中的方位角先验高斯分布求取结果图；
图5是本说明书一个实施例中的自适应校准流程图；
图6是本说明书一个实施例中的目标井的井斜角预测结果图；
图7是本说明书提供的一种基于粒子滤波的井眼轨迹预测装置一个实施例的模块结构示意图。
具体实施方式
为了使本技术领域的人员更好地理解本说明书中的技术方案，下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明书保护的范围。
考虑到现有的基于统计方法的井眼轨迹预测方法通常利用邻井或已钻井眼数据建立回归模型，通过回归系数反映待钻井眼的地层特性，但是它难以根据邻井数据及时调整自身模型参数，导致其自适应能力较差。
进一步，还考虑到现有的基于几何外推的井眼轨迹预测方法没有考虑地层特性对井眼轨迹的影响，导致其计算结果准确性较差。
针对现有方法存在的上述问题以及产生上述问题的具体原因，本申请考虑可以引入粒子滤波进行井眼轨迹预测以提升方法的自适应能力和预测结果精度。
基于上述思路，本说明书提出一种基于粒子滤波的井眼轨迹预测方法。首先，获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；然后，根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；最后，基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。参阅图1所示，本说明书实施例提供了一种基于粒子滤波的井眼轨迹预测方法。具体实施时，该方法可以包括以下内容。
S101：获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据。
在一些实施例中，所述邻井数据具体可以包括：轨迹测量数据、录井米数据、定向井施工记录数据、邻井参数数据；所述预处理数据具体可以包括：钻进数据、工具面角平均值、井眼第二分段结果。
在一些实施例中，所述轨迹测量数据具体可以包括：轨迹测量深度、井斜角、方位角；所述录井米数据具体可以包括：录井米深度、钻压、转速；所述定向井施工记录数据具体可以包括工具面角；所述邻井参数数据具体可以包括：井眼尺寸、地层、钻具组合；所述钻进数据具体可以包括：滑动钻进距离、滑动钻进平均钻压、复合钻进距离。
在一些实施例中，上述基于邻井数据获得预处理数据，具体实施时，可以包括：
S1：基于轨迹测量数据对井眼进行第一次分段，得到井眼第一分段结果；
S2：根据所述井眼第一分段结果和录井米数据，得到钻进数据；
S3：根据所述井眼第一分段结果和定向井施工记录数据，得到工具面角平均值；
S4：基于邻井参数数据对井眼进行第二次分段，得到井眼第二分段结果。
在一些实施例中，上述基于轨迹测量数据对井眼进行第一次分段，得到井眼第一分段结果，具体实施时，可以包括：依据轨迹测量数据对井眼轨迹分段，将某测点的轨迹测量深度、井斜角、方位角分别作为该井段的底部深度、底部井斜角、底部方位角，将上一个测点的轨迹测量深度、井斜角、方位角分别作为该井段的顶部深度、顶部井斜角、顶部方位角，得到了井眼第一分段结果。
在一些实施例中，上述根据所述井眼第一分段结果和录井米数据，得到钻进数据，具体实施时，可以包括：根据井眼第一分段结果，在录井米数据中提取录井米深度，将每个录井米深度放入对应的第一分段结果中，并提取每个第一分段结果中转速等于0的记录行数，记为滑动钻进距离，并计算转速等于0的行中钻压的平均值，记为滑动钻进平均钻压；提取每个第一分段结果中转速大于0的记录行数，记为复合钻进距离。
在一些实施例中，上述根据所述井眼第一分段结果和定向井施工记录数据，得到工具面角平均值，具体实施时，可以包括：根据定向井施工记录数据获取井眼第一分段结果中每一段的工具面角，并计算其平均值，记为工具面角平均值。
在一些实施例中，上述基于邻井参数数据对井眼进行第二次分段，得到井眼第二分段结果，具体实施时，可以包括：根据井眼尺寸划分井眼第一分段结果，得到第一井段数据，其中，第一井段数据中的每个井段都具有相同的井眼尺寸；根据地层划分第一井段数据，得到第二井段数据，其中，第二井段数据中的每个井段都具有相同的地层；根据钻具组合划分第二井段数据，得到井眼第二分段结果，其中，井眼第二分段结果中的每个井段都具有相同的钻具组合；其中，上述井眼第二分段结果中的最小井段内测点数目不能少于10个；并且，将井眼第二分段结果中每段内最后一个测点的井斜角作为该段的底部井斜角，同时也是下一段的顶部井斜角；将井眼第二分段结果中每段内最后一个测点的方位角作为该段的底部方位角，同时也是下一段的顶部方位角。
在一些实施例中，在计算机的软件交互页面针对上述预处理数据进行可视化显示。
通过上述实施例，同时考虑到了井眼尺寸、地层、钻具组合三种因素对井眼进行分段，可以使分段结果更加准确，提高后续粒子滤波方法计算的精度。
S102：根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程。
在一些实施例中，状态方程用于将***当前时刻的状态预测结果和状态的测量相结合，得到对***的最优估计；观测方程用于对状态进行线性或非线性观测观测，得到***当前时刻的观测值。本说明书中，可以使用前一个井段的井眼轨迹参数的观测值结合预处理数据来更正待预测井段的井眼轨迹参数的值，以得到更加准确的观测值；使用前一个井段的状态参数值结合高斯白噪声推算待预测井段的状态参数值。
在一些实施例中，上述根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程，具体实施时，可以包括：
按照以下算式构建井斜角观测方程：
Incb，i＝(θ0，i×cos(tfb，i)+θ1，i×WOBb，i+θ2，i)×Lb，i+θ3，i×Lh，i+θ4，i+Incb，i-1 (1)
式中，i和i-1表示待预测的井段的编号，Incb，i为待预测井段底部井斜角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Incb，i-1为待预测井段顶部井斜角，θ0，i、θ1，i、θ2，i、θ3，i、θ4，i为i井段的井斜角状态参数，其中，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，ii井段的为井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差。
按照以下算式构建井斜角状态方程：
Figure BDA0003711977730000081
式中，w0，i-1、w1，i-1、w2，i-1、w3，i-1、w4，i-1为i-1井段的井斜角高斯白噪声，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，i为i井段的井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差，θ0，i-1为i-1井段的井斜角造斜率修正系数，θ1，i-1为i-1井段的井斜角钻压影响系数，θ2，i-1为i-1井段的井斜角造斜率误差，θ3，i-1为i-1井段的井斜角复合钻进距离斜率变化规律系数，θ4，i-1为i-1井段的井斜角***误差。
按照以下算式构建方位角观测方程：
Azib，i＝(φ0，i×cos(tfb，i)+φ1，i×WOBb，i+φ2，i)×Lb，i+φ3，i×Lh，i+φ4，i+Azib，i-1(3)
式中，i和i-1表示待预测的井段的编号，Azib，i为待预测井段底部方位角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Azib，i-1为待预测井段顶部方位角，φ0，i、φ1，i、φ2，i、φ3，i、φ1，i为i井段的方位角状态参数，其中，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差。
按照以下算式构建方位角状态方程：
Figure BDA0003711977730000082
Figure BDA0003711977730000091
式中，n0，i-1、n1，i-1、n2，i-1、n3，i-1、n4，i-1为i-1井段的方位角高斯白噪声，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差，φ0，i-1为i-1井段的方位角造斜率修正系数，φ1，i-1为i-1井段的方位角钻压影响系数，φ2，i-1为i-1井段的方位角造斜率误差，φ3，i-1为i-1井段的方位角复合钻进距离斜率变化规律系数，φ4，i-1为i-1井段的方位角***误差。
S103：基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
在一些实施例中，所述目标井的井眼轨迹参数可以包括：待预测井段底部最终井斜角、待预测井段底部最终方位角；所述目标井的井眼轨迹参数置信区间可以包括：待预测井段底部井斜角预设的预测置信区间、待预测井段底部方位角预设的预测置信区间。
在一些实施例中，粒子滤波是一种基于蒙特卡洛思想的非线性、非高斯***滤波方法，它改进了卡尔曼滤波只适用于线性***的缺陷，对***的过程噪声和测量噪声没有任何限制。粒子滤波用一系列带有权值的粒子来表示后验概率密度函数，从而得到待预测的值。由于粒子滤波同时适用于非线性和非高斯的***，所以它的应用范围更广。本说明书中，首先采用粒子滤波更新状态参数，将前一个井段的后验分布作为当前井段的先验分布，利用更新后的状态参数解算井斜角观测方程和方位角观测方程，得到待预测井段底部井斜角和待预测井段底部方位角；然后对粒子群进行权重归一化、重采样操作，得到重采样后的粒子群；最后根据待预测井段底部井斜角、待预测井段底部方位角、重采样后的粒子群，计算待预测井段的井眼轨迹参数和井眼轨迹参数置信区间，以实现井眼轨迹参数的在线自适应校准，然后可以对定向施工参数进行调整；上述后验分布具体指的是前一个井段经过重采样后的粒子群。
在一些实施例中，上述通过粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间，具体实施时，可以包括：
按照以下方式，通过粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，确定目标井中当前待预测井段的井眼轨迹参数和井眼轨迹参数置信区间：
S1：根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数；
S2：根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角；
S3：根据当前待预测井段的权重归一化后的粒子群，得到当前待预测井段的重采样粒子群；
S4：基于当前待预测井段的重采样粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，获得当前待预测井段的底部最终井斜角、当前待预测井段的底部最终方位角、当前待预测井段底部井斜角预设的预测置信区间、当前待预测井段底部方位角预设的预测置信区间。
在一些实施例中，上述根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数，具体实施时，可以包括：
S1：判断当前待预测井段是否为第一个井段；
S2：在当前待预测井段是第一个井段的情况下，依据预处理数据获得待预测目标井的先验高斯分布，并基于先验高斯分布生成初始粒子群；依据初始粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数；
S3：在当前待预测井段不是第一个井段的情况下，获取前一个井段重采样后的粒子群，依据前一个井段重采样后的粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数。
在一些实施例中，上述状态参数具体可以包括：井斜角状态参数和方位角状态参数。
在一些实施例中，上述初始粒子群具体可以包括：井斜角初始粒子群和方位角初始粒子群；所述井斜角初始粒子群和所述方位角初始粒子群均为一组随机产生的具有上述先验高斯分布的数据集合，井斜角初始粒子群中每个粒子的权重相等，方位角初始粒子群中每个粒子的权重相等。
其中，上述井斜角初始粒子群可以记为
Figure BDA0003711977730000101
其为N×5的矩阵，N表示井斜角初始粒子群的粒子数，每个粒子的权重可以记为1/N；上述方位角初始粒子群可以记为
Figure BDA0003711977730000102
其为M×5的矩阵，M表示方位角初始粒子群的粒子数，每个粒子的权重可以记为1/M。
在一些实施例中，上述在当前待预测井段是第一个井段的情况下，依据预处理数据获得待预测目标井的先验高斯分布，具体实施时，可以包括：
S1：根据井眼第二分段结果，基于顶部井斜角、滑动钻进距离、滑动钻进平均钻压、复合钻进平均钻压、底部井斜角，拟合井斜角观测方程，得到每一个井段的井斜角状态参数；
S2：根据井眼第二分段结果，基于顶部方位角、滑动钻进距离、滑动钻进平均钻压、复合钻进平均钻压、底部方位角，拟合方位角观测方程，得到每一个井段的方位角状态参数；
S3：计算每一个井段的井斜角状态参数的期望和方差、每一个井段的方位角状态参数的期望和方差，得到待预测目标井的先验高斯分布。
在一些实施例中，所述拟合井斜角观测方程和所述拟合方位角观测方程的方法具体可以包括：最小二乘法、牛顿迭代法、线性逼近约束优化(COBYLA)法。
在一些实施例中，上述依据初始粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数，具体实施时，可以包括：
S1：将井斜角初始粒子群第1列数据作为井斜角造斜率修正系数，第2列数据作为井斜角钻压影响系数，第3列数据作为井斜角造斜率误差，第4列数据作为井斜角复合钻进距离斜率变化规律系数，第5列数据作为井斜角***误差，求解井斜角状态方程，得到当前待预测井段的井斜角状态参数；
S2：将方位角初始粒子群第1列数据作为方位角造斜率修正系数，第2列数据作为方位角钻压影响系数，第3列数据作为方位角造斜率误差，第4列数据作为方位角复合钻进距离斜率变化规律系数，第5列数据作为方位角***误差，求解方位角状态方程，得到当前待预测井段的方位角状态参数。
在一些实施例中，上述依据前一个井段重采样后的粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数，具体实施时，可以包括：
S1：将前一个井段重采样后的井斜角粒子群第1列数据作为前一个井段的井斜角造斜率修正系数，第2列数据作为前一个井段的井斜角钻压影响系数，第3列数据作为前一个井段的井斜角造斜率误差，第4列数据作为前一个井段的井斜角复合钻进距离斜率变化规律系数，第5列数据作为前一个井段的井斜角***误差，将上述前一个井段重采样后的井斜角粒子群代入井斜角状态方程，求解得到当前待预测井段的井斜角状态参数；
S2：将前一个井段重采样后的方位角粒子群第1列数据作为前一个井段的方位角造斜率修正系数，第2列数据作为前一个井段的方位角钻压影响系数，第3列数据作为前一个井段的方位角造斜率误差，第4列数据作为前一个井段的方位角复合钻进距离斜率变化规律系数，第5列数据作为前一个井段的方位角***误差，将上述前一个井段重采样后的方位角粒子群代入方位角状态方程，求解得到当前待预测井段的方位角状态参数。
在一些实施例中，上述根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，具体实施时，可以包括以下步骤：
S1：依据当前待预测井段的状态参数，求解井斜角观测方程和方位角观测方程，得到当前待预测井段的底部井斜角和当前待预测井段的底部方位角；
S2：依据当前待预测井段的底部井斜角和当前待预测井段的底部方位角，得到当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差；
S3：根据当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差重新计算当前待预测井段的粒子群中每一个粒子的权重，获得当前待预测井段的重新赋权的粒子群；
S4：依据当前待预测井段的重新赋权的粒子群，获得当前待预测井段的权重归一化后的粒子群。
在一些实施例中，上述重新赋权的粒子群，具体包括：井斜角重新赋权的粒子群、方位角重新赋权的粒子群。
在一些实施例中，上述权重归一化后的粒子群，具体包括：井斜角权重归一化后的粒子群、方位角权重归一化后的粒子群。
在一些实施例中，上述依据当前待预测井段的底部井斜角和当前待预测井段的底部方位角，得到当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差，具体实施时，可以包括：
按照以下算式计算当前待预测井段的井斜角绝对误差：
Figure BDA0003711977730000121
式中，
Figure BDA0003711977730000122
为当前待预测井段的井斜角绝对误差，
Figure BDA0003711977730000123
为第k个当前待预测井段底部井斜角，
Figure BDA0003711977730000124
为第k个当前待预测井段底部井斜角的实测值。
按照以下算式计算当前待预测井段的方位角绝对误差：
Figure BDA0003711977730000125
式中，
Figure BDA0003711977730000126
为当前待预测井段的方位角绝对误差，
Figure BDA0003711977730000127
为第l个当前待预测井段底部方位角，
Figure BDA0003711977730000128
为第l个当前待预测井段底部方位角的实测值。
在一些实施例中，上述根据当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差重新计算当前待预测井段的粒子群中每一个粒子的权重，获得当前待预测井段的重新赋权的粒子群，具体实施时，包括：
按照以下算式计算当前待预测井段的井斜角重新赋权的粒子群的权重：
Figure BDA0003711977730000129
式中，
Figure BDA00037119777300001210
为当前待预测井段的井斜角重新赋权的粒子群中第k个粒子的权重，RInc为常数，根据经验确定。
按照以下算式计算当前待预测井段的方位角重新赋权的粒子群的权重：
Figure BDA0003711977730000131
式中，
Figure BDA0003711977730000132
为当前待预测井段的方位角重新赋权的粒子群中第l个粒子的权重，RAzi为常数，根据经验确定。
在一些实施例中，上述依据当前待预测井段的重新赋权的粒子群，获得当前待预测井段的权重归一化后的粒子群，具体实施时，可以包括：
按照以下算式计算当前待预测井段的井斜角权重归一化后的粒子群的权重：
Figure BDA0003711977730000133
式中，
Figure BDA0003711977730000134
为井斜角权重归一化后的粒子群中第k个粒子的权重。
按照以下算式计算当前待预测井段的方位角权重归一化后的粒子群的权重：
Figure BDA0003711977730000135
式中，
Figure BDA0003711977730000136
为方位角权重归一化后的粒子群中第l个粒子的权重。
在一些实施例中，上述根据当前待预测井段的权重归一化后的粒子群，得到当前待预测井段的重采样粒子群，具体实施时，可以包括：对当前待预测井段的权重归一化后的粒子群进行按照权重大小进行重采样，权重越高的粒子重采样次数越多，经过重采样之后，粒子总数不变，权重大的粒子分成了多个粒子，权重小的粒子被抛弃，重采样之后，每个粒子的权重相同；其中，所述重采样的方法具体可以包括：随机重采样、***采样、残差采样、多项式采样；所述重采样粒子群具体可以包括：井斜角重采样粒子群、方位角重采样粒子群。
在一些实施例中，上述基于当前待预测井段的重采样粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，获得当前待预测井段的底部最终井斜角、当前待预测井段的底部最终方位角、当前待预测井段底部井斜角预设的预测置信区间、当前待预测井段底部方位角预设的预测置信区间，具体实施时，可以包括：
按照以下算式计算当前待预测井段的底部最终井斜角：
Figure BDA0003711977730000137
式中，Incfinal,b,i为当前待预测井段的底部最终井斜角，
Figure BDA0003711977730000138
为当前待预测井段的井斜角重采样粒子群中第k个粒子的值。
按照以下算式计算当前待预测井段的底部最终方位角：
Figure BDA0003711977730000139
式中，Azifinal,b,i为当前待预测井段的底部最终方位角，
Figure BDA00037119777300001310
为当前待预测井段的方位角重采样粒子群中第l个粒子的值。
在一些实施例中，在计算机的软件交互页面针对上述待预测井段的底部最终井斜角、上述待预测井段的底部最终方位角进行可视化显示。
在一些实施例中，在得到当前待预预测井段的井眼轨迹参数之后，可以利用该井眼轨迹参数进行井眼轨迹的预测，即利用井眼轨迹参数依据轨迹形状计算相应测点的位置坐标，最后得到井眼轨迹形态的空间数据，现场工作人员依据空间数据进行决策思考，为指挥钻井工作提供理论依据。
在所述当前待预测井段底部井斜角预设的预测置信区间可以为当前待预测井段底部井斜角95％置信区间的情况下，按照以下算式计算当前待预测井段底部井斜角预设的预测置信区间下限：
Incpred,min＝μInc-1.96σInc (13)
式中，μInc为待预测井段底部井斜角的均值，σInc为待预测井段底部井斜角的方差，Incpred,min为待预测井段底部井斜角预设的预测置信区间下限。
在所述当前待预测井段底部井斜角预设的预测置信区间可以为当前待预测井段底部井斜角95％置信区间的情况下，按照以下算式计算当前待预测井段底部井斜角预设的预测置信区间上限：
Incpred,max＝μInc+1.96σInc (14)
Incpred,max为待预测井段底部井斜角预设的预测置信区间上限。
因此，待预测井段底部井斜角预设的预测置信区间可以记为[Incpred,min,Incpred,max]。
在所述当前待预测井段底部方位角预设的预测置信区间可以为当前待预测井段底部方位角95％置信区间的情况下，按照以下算式计算当前待预测井段底部方位角预设的预测置信区间下限：
Azipred,min＝μAzi-1.96σAzi (15)
式中，μAzi为待预测井段底部方位角的均值，σAzi为待预测井段底部方位角的方差，Azipred,min为待预测井段底部方位角预设的预测置信区间下限。
在所述当前待预测井段底部方位角预设的预测置信区间可以为当前待预测井段底部方位角95％置信区间的情况下，按照以下算式计算当前待预测井段底部方位角预设的预测置信区间上限：
Azipred,max＝μAzi+1.96σAzi (16)
式中，Azipred,min为待预测井段底部方位角预设的预测置信区间上限。
因此，待预测井段底部方位角预设的预测置信区间可以记为[Azipred,min,Azipred,max]。
在一些实施例中，根据待预测井段底部井斜角预设的预测置信区间、待预测井段底部方位角预设的预测置信区间就可以对目标井的井眼轨迹进行不确定描述。
在一个具体的场景示例中，可以应用本说明书提供的基于粒子滤波的井眼轨迹预测方法对目标井眼的井眼轨迹参数进行求取。其中，基于邻井数据得到的预处理数据如图2所示；目标井一共被分为4个井段，井斜角状态参数先验高斯分布如图3所示，θi,prior(i＝0,1,2,3,4)表示井斜角状态参数先验高斯分布；方位角状态参数先验高斯分布如图4所示，φi,prior(i＝0,1,2,3,4)表示方位角状态参数先验高斯分布；自适应校准流程如图5所示，依据邻井数据获得目标井的先验高斯分布，并且将前一个井段的重采样粒子群作为前一个井段的后验分布，同时也作为当前待遇测井段的先验分布，依据粒子滤波的计算结果对井眼轨迹预测结果进行自适应校准，根据井眼轨迹预测结果，指导基于随钻数据的钻井作业，在自适应校准过程中，粒子滤波的相关参数可以进行自适应调整；井斜角预测结果如图6所示，由图6可知基于粒子滤波的井眼轨迹预测方法能取得较好的预测效果。
虽然本说明书提供了如下述实施例或附图所示的方法操作步骤或装置结构，但基于常规或者无需创造性的劳动在所述方法或装置中可以包括更多或者部分合并后更少的操作步骤或模块单元。在逻辑性上不存在必要因果关系的步骤或结构中，这些步骤的执行顺序或装置的模块结构不限于本说明书实施例或附图所示的执行顺序或模块结构。所述的方法或模块结构的在实际中的装置、服务器或终端产品应用时，可以按照实施例或者附图所示的方法或模块结构进行顺序执行或者并行执行(例如并行处理器或者多线程处理的环境、甚至包括分布式处理、服务器集群的实施环境)。
基于上述基于粒子滤波的井眼轨迹预测方法，本说明书还提出一种基于粒子滤波的井眼轨迹预测装置的实施例。如图7所示，所述基于粒子滤波的井眼轨迹预测装置具体包括以下模块：
预处理模块701，用于获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；
构建模块702，用于根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；
粒子滤波模块703，用于基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
在一些实施例中，上述预处理模块701具体可以用于基于轨迹测量数据对井眼进行第一次分段，得到井眼第一分段结果；根据所述井眼第一分段结果和录井米数据，得到钻进数据；根据所述井眼第一分段结果和定向井施工记录数据，得到工具面角平均值；基于邻井参数数据对井眼进行第二次分段，得到井眼第二分段结果。
在一些实施例中，上述构建模块702具体可以用于：
按照以下算式构建井斜角观测方程：
Incb，i＝(θ0，i×cos(tfb，i)+θ1，i×WOBb，i+θ2，i)×Lb，i+θ3，i×Lh，i+θ4，i+Incb，i-1 (17)
式中，i和i-1表示待预测的井段的编号，b表示底部，Incb，i为待预测井段底部井斜角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Incb，i-1为待预测井段顶部井斜角，θ0，i、θ1，i、θ2，i、θ3，i、θ4，i为i井段的井斜角状态参数，其中，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，ii井段的为井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差。
按照以下算式构建井斜角状态方程：
Figure BDA0003711977730000161
式中，w0，i-1、w1，i-1、w2，i-1、w3，i-1、w4，i-1为i-1井段的井斜角高斯白噪声，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，i为i井段的井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差，θ0，i-1为i-1井段的井斜角造斜率修正系数，θ1，i-1为i-1井段的井斜角钻压影响系数，θ2，i-1为i-1井段的井斜角造斜率误差，θ3，i-1为i-1井段的井斜角复合钻进距离斜率变化规律系数，θ4，i-1为i-1井段的井斜角***误差。
按照以下算式构建方位角观测方程：
Azib，i＝(φ0，i×cos(tfb，i)+φ1，i×WOBb，i+φ2，i)×Lb，i+φ3，i×Lh，i+φ4，i+Azib，i-1(19)
式中，i和i-1表示待预测的井段的编号，b表示底部，Azib，i为待预测井段底部方位角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Azib，i-1为待预测井段顶部方位角，φ0，i、φ1，i、φ2，i、φ3，i、φ1，i为i井段的方位角状态参数，其中，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差。
按照以下算式构建方位角状态方程：
Figure BDA0003711977730000171
式中，n0，i-1、n1，i-1、n2，i-1、n3，i-1、n4，i-1为i-1井段的方位角高斯白噪声，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差，φ0，i-1为i-1井段的方位角造斜率修正系数，φ1，i-1为i-1井段的方位角钻压影响系数，φ2，i-1为i-1井段的方位角造斜率误差，φ3，i-1为i-1井段的方位角复合钻进距离斜率变化规律系数，φ4，i-1为i-1井段的方位角***误差。
在一些实施例中，上述粒子滤波模块703具体可以用于按照以下方式，通过粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，确定目标井中当前待预测井段的井眼轨迹参数和井眼轨迹参数置信区间：根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数；根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角；根据当前待预测井段的权重归一化后的粒子群，得到当前待预测井段的重采样粒子群；基于当前待预测井段的重采样粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，获得当前待预测井段的底部最终井斜角、当前待预测井段的底部最终方位角、当前待预测井段底部井斜角预设的预测置信区间、当前待预测井段底部方位角预设的预测置信区间。
本说明书实施例还提供了一种基于粒子滤波的井眼轨迹预测方法的计算机存储介质，所述计算机存储介质存储有计算机程序指令，在所述计算机程序指令被执行时实现：获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
在本实施例中，上述存储介质包括但不限于随机存取存储器(Random AccessMemory，RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(Hard DiskDrive,HDD)或者存储卡(Memory Card)。所述存储器可以用于存储计算机程序指令。网络通信单元可以是依照通信协议规定的标准设置的，用于进行网络连接通信的接口。
虽然本说明书提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或客户端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。第一，第二等词语用来表示名称，而并不表示任何特定的顺序。
本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构、类等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
通过以上的实施例的描述可知，本领域的技术人员可以清楚地了解到本说明书可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本说明书的技术方案本质上可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，移动终端，服务器，或者网络设备等)执行本说明书各个实施例或者实施例的某些部分所述的方法。
虽然通过实施例描绘了本说明书，本领域普通技术人员知道，本说明书有许多变形和变化而不脱离本说明书的精神，希望所附的权利要求包括这些变形和变化而不脱离本说明书的精。

Claims (10)

1.一种基于粒子滤波的井眼轨迹预测方法，其特征在于，包括：
获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；
根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；
基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
2.根据权利要求1所述的方法，其特征在于，基于邻井数据获得预处理数据，包括：
基于轨迹测量数据对目标井的井眼进行第一次分段，得到井眼第一分段结果；
根据所述井眼第一分段结果和录井米数据，得到钻进数据；
根据所述井眼第一分段结果和定向井施工记录数据，得到工具面角平均值；
基于邻井参数数据对井眼进行第二次分段，得到井眼第二分段结果。
3.根据权利要求1所述的方法，其特征在于，根据预处理数据，构建井斜角观测方程，包括：
按照以下算式构建井斜角观测方程：
Incb，i＝(θ0，i×cos(tfb，i)+θ1，i×WOBb，i+θ2，i)×Lb，i+θ3，i×Lh，i+θ4，i+Incb，i-1
式中，i和i-1表示待预测的井段的编号，b表示底部，Incb，i为待预测井段底部井斜角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Incb，i-1为待预测井段顶部井斜角，θ0，i、θ1，i、θ2，i、θ3，i、θ4，i为i井段的井斜角状态参数，其中，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，i井段的为井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差。
4.根据权利要求1所述的方法，其特征在于，根据预处理数据，构建井斜角状态方程，包括：
按照以下算式构建井斜角状态方程：
θ0，i＝θ0，i-1+w0，i-1
θ1，i＝θ1，i-1+w1，i-1
θ2，i＝θ2，i-1+w2，i-1
θ3，i＝θ3，i-1+w3，i-1
θ4，i＝θ4，i-1+w4，i-1
式中，w0，i-1、w1，i-1、w2，i-1、w3，i-1、w4，i-1为i-1井段的井斜角高斯白噪声，θ0，i为i井段的井斜角造斜率修正系数，θ1，i为i井段的井斜角钻压影响系数，θ2，i为i井段的井斜角造斜率误差，θ3，i为i井段的井斜角复合钻进距离斜率变化规律系数，θ4，i为i井段的井斜角***误差，θ0，i-1为i-1井段的井斜角造斜率修正系数，θ1，i-1为i-1井段的井斜角钻压影响系数，θ2，i-1为i-1井段的井斜角造斜率误差，θ3，i-1为i-1井段的井斜角复合钻进距离斜率变化规律系数，θ4，i-1为i-1井段的井斜角***误差。
5.根据权利要求1所述的方法，其特征在于，根据预处理数据，构建方位角观测方程，包括：
按照以下算式构建方位角观测方程：
Azib，i＝(φ0，i×cos(tfb，i)+φ1，i×WOBb，i+φ2，i)×Lb，i+φ3，i×Lh，i+φ4，i+Azib，i-1
式中，i和i-1表示待预测的井段的编号，b表示底部，Azib，i为待预测井段底部方位角，tfb，i为工具面角平均值，WOBb，i为滑动钻进平均钻压，Lb，i为滑动钻进距离，Lh，i为复合钻进距离，Azib，i-1为待预测井段顶部方位角，φ0，i、φ1，i、φ2，i、φ3，i、φ1，i为i井段的方位角状态参数，其中，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差。
6.根据权利要求1所述的方法，其特征在于，根据预处理数据，构建方位角状态方程，包括：
按照以下算式构建方位角状态方程：
φ0，i＝φ0，i-1+n0，i-1
φ1，i＝φ1，i-1+n1，i-1
φ2，i＝φ2，i-1+n2，i-1
φ3，i＝φ3，i-1+n3，i-1
φ4，i＝φ4，i-1+n4，i-1
式中，n0，i-1、n1，i-1、n2，i-1、n3，i-1、n4，i-1为i-1井段的方位角高斯白噪声，φ0，i为i井段的方位角造斜率修正系数，φ1，i为i井段的方位角钻压影响系数，φ2，i为i井段的方位角造斜率误差，φ3，i为i井段的方位角复合钻进距离斜率变化规律系数，φ4，i为i井段的方位角***误差，φ0，i-1为i-1井段的方位角造斜率修正系数，φ1，i-1为i-1井段的方位角钻压影响系数，φ2，i-1为i-1井段的方位角造斜率误差，φ3，i-1为i-1井段的方位角复合钻进距离斜率变化规律系数，φ4，i-1为i-1井段的方位角***误差。
7.根据权利要求1所述的方法，其特征在于，通过粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间，包括：
按照以下方式，通过粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，确定目标井中当前待预测井段的井眼轨迹参数和井眼轨迹参数置信区间：
根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数；
根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角；
根据当前待预测井段的权重归一化后的粒子群，得到当前待预测井段的重采样粒子群；
基于当前待预测井段的重采样粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，获得当前待预测井段的底部最终井斜角、当前待预测井段的底部最终方位角、当前待预测井段底部井斜角预设的预测置信区间、当前待预测井段底部方位角预设的预测置信区间。
8.根据权利要求7所述的方法，其特征在于，根据井斜角状态方程以及方位角状态方程，求取当前待预测井段的状态参数，包括：
判断当前待预测井段是否为第一个井段；
在当前待预测井段是第一个井段的情况下，依据预处理数据获得待预测目标井的先验高斯分布，并基于先验高斯分布生成初始粒子群；依据初始粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数；
在当前待预测井段不是第一个井段的情况下，获取前一个井段重采样后的粒子群，依据前一个井段重采样后的粒子群求解井斜角状态方程以及方位角状态方程，得到当前待预测井段的状态参数。
9.根据权利要求7所述的方法，其特征在于，根据当前待预测井段的状态参数，求解井斜角观测方程以及方位角观测方程，得到当前待预测井段的权重归一化后的粒子群、当前待预测井段的底部井斜角、当前待预测井段的底部方位角，包括：
依据当前待预测井段的状态参数，求解井斜角观测方程和方位角观测方程，得到当前待预测的井段底部井斜角和当前待预测井段的底部方位角；
依据当前待预测井段的底部井斜角和当前待预测井段的底部方位角，得到当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差；
根据当前待预测井段的井斜角绝对误差和当前待预测井段的方位角绝对误差重新计算当前待预测井段的粒子群中每一个粒子的权重，获得当前待预测井段的重新赋权的粒子群；
依据当前待预测井段的重新赋权的粒子群，获得当前待预测井段的权重归一化后的粒子群。
10.一种基于粒子滤波的井眼轨迹预测装置，其特征在于，包括：
预处理模块，用于获取待预测目标井周边的邻井数据，并基于邻井数据获得预处理数据；
构建模块，用于根据预处理数据，构建井斜角观测方程和状态方程，以及方位角观测方程和状态方程；
粒子滤波模块，用于基于粒子滤波，求解井斜角观测方程和状态方程，以及方位角观测方程和状态方程，以确定目标井的井眼轨迹参数和井眼轨迹参数置信区间。
CN202210723052.6A 2022-06-24 2022-06-24 一种基于粒子滤波的井眼轨迹预测方法及装置 Active CN115095318B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210723052.6A CN115095318B (zh) 2022-06-24 2022-06-24 一种基于粒子滤波的井眼轨迹预测方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210723052.6A CN115095318B (zh) 2022-06-24 2022-06-24 一种基于粒子滤波的井眼轨迹预测方法及装置

Publications (2)

Publication Number Publication Date
CN115095318A true CN115095318A (zh) 2022-09-23
CN115095318B CN115095318B (zh) 2023-05-23

Family

ID=83292309

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210723052.6A Active CN115095318B (zh) 2022-06-24 2022-06-24 一种基于粒子滤波的井眼轨迹预测方法及装置

Country Status (1)

Country Link
CN (1) CN115095318B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002129882A (ja) * 2000-10-25 2002-05-09 Nippon Telegr & Teleph Corp <Ntt> 埋設管推進機の方向制御方法および装置
CA2636564A1 (en) * 2008-06-25 2009-12-25 Uti Limited Partnership In-drilling alignment
US20100100329A1 (en) * 2008-10-22 2010-04-22 Gyrodata, Incorporated Downhole surveying utilizing multiple measurements
US20100126770A1 (en) * 2008-11-24 2010-05-27 Pathfinder Energy Services, Inc. Non-Azimuthal and Azimuthal Formation Evaluation Measurement in a Slowly Rotating Housing
CN109915018A (zh) * 2017-12-07 2019-06-21 中国石油化工股份有限公司华北油气分公司石油工程技术研究院 一种三维井眼轨迹控制中工具面角的确定方法
US20190257189A1 (en) * 2018-02-19 2019-08-22 Gyrodata, Incorporated Determining Direct Hit or Unintentional Crossing Probabilities for Wellbores
CN111119855A (zh) * 2018-11-01 2020-05-08 中国石油化工股份有限公司 基于地层和钻具特性的井眼轨迹预测方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002129882A (ja) * 2000-10-25 2002-05-09 Nippon Telegr & Teleph Corp <Ntt> 埋設管推進機の方向制御方法および装置
CA2636564A1 (en) * 2008-06-25 2009-12-25 Uti Limited Partnership In-drilling alignment
US20100100329A1 (en) * 2008-10-22 2010-04-22 Gyrodata, Incorporated Downhole surveying utilizing multiple measurements
US20100126770A1 (en) * 2008-11-24 2010-05-27 Pathfinder Energy Services, Inc. Non-Azimuthal and Azimuthal Formation Evaluation Measurement in a Slowly Rotating Housing
CN109915018A (zh) * 2017-12-07 2019-06-21 中国石油化工股份有限公司华北油气分公司石油工程技术研究院 一种三维井眼轨迹控制中工具面角的确定方法
US20190257189A1 (en) * 2018-02-19 2019-08-22 Gyrodata, Incorporated Determining Direct Hit or Unintentional Crossing Probabilities for Wellbores
CN111119855A (zh) * 2018-11-01 2020-05-08 中国石油化工股份有限公司 基于地层和钻具特性的井眼轨迹预测方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
MIAO YINSEN 等: "-Nonlinear state-space modeling approaches to real-time autonomous geosteering" *
丁华,王秀坤,孙焘: "基于粒子群算法的井眼轨迹优化研究", 大连理工大学学报 *
刘修善;王超;: "空间圆弧轨迹的解析描述技术", 石油学报 *
宋晓健等: "基于数据融合的近钻头井眼轨迹参数动态测量方法" *
沙林秀;潘仲奇;: "基于NPSO的三维复杂井眼轨迹控制转矩的优选", 石油机械 *
王越之,王用军: "用灰色动态GM(1,1)模型预测定向井井眼轨迹", 石油天然气学报 *

Also Published As

Publication number Publication date
CN115095318B (zh) 2023-05-23

Similar Documents

Publication Publication Date Title
AU2017208259B2 (en) Controlling wellbore operations
EP2929141B1 (en) Weighting function for inclination and azimuth computation
US10359541B2 (en) Creating virtual production logging tool profiles for improved history matching
BR112020022959A2 (pt) Taxa de penetração de ferramenta de perfuração terrestre e sistema de previsão dedesgaste e métodos relacionados
CN111989456B (zh) 用于井筒操作中定向钻井的基于模型的参数估计
WO2020080973A1 (ru) Способ и система комбинированного сопровождения процесса бурения скважины
CN112162316B (zh) 一种avo波形数据驱动的高分辨率井震融合叠前反演方法
CN115329657B (zh) 钻井参数优化方法及装置
US20230195952A1 (en) Well planning based on hazard predictive models
EP3186478B1 (en) Directional driller quality reporting system and method
RU2720115C1 (ru) Способ автоматизированного процесса геологической проводки скважин и система для его осуществления
CN110737977A (zh) 一种隧道围岩变形预测方法及预测装置
CN107273659A (zh) 一种基于ransac算法改进的用于空间碎片光电跟踪的轨迹预测方法
CN115095318A (zh) 一种基于粒子滤波的井眼轨迹预测方法及装置
CN107780926B (zh) 井眼截面形状测量方法、装置和终端
RU2745152C1 (ru) Способ объединения модели геологической проводки скважины с оперативной петрофизической интерпретацией данных гис в режиме реального времени и система, реализующая способ
CN111856573B (zh) 一种古应力反演方法
CN114488297A (zh) 一种断层识别方法及装置
EP3449290A1 (en) Method for wellbore survey instrument fault detection
CN115162958B (zh) 一种定向钻孔控制方法及***
US11156740B2 (en) Electrofacies determination
CN113027427B (zh) 一种基于贝叶斯理论的含可信度地层压力随钻修正方法
WO2022071898A1 (en) Method of determining boring points using geostatistical interpolation method
CN110967756B (zh) 基于正态分布的微地震定位精度评价方法及***
US20210396127A1 (en) Estimating borehole shape between stationary survey locations

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant