CN110886606B - 一种随钻特征量辅助的惯性测斜方法及装置 - Google Patents

一种随钻特征量辅助的惯性测斜方法及装置 Download PDF

Info

Publication number
CN110886606B
CN110886606B CN201911141926.1A CN201911141926A CN110886606B CN 110886606 B CN110886606 B CN 110886606B CN 201911141926 A CN201911141926 A CN 201911141926A CN 110886606 B CN110886606 B CN 110886606B
Authority
CN
China
Prior art keywords
coordinate system
drilling
attitude
speed
information
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.)
Active
Application number
CN201911141926.1A
Other languages
English (en)
Other versions
CN110886606A (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 University of Geosciences Beijing
Original Assignee
China University of Geosciences Beijing
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 University of Geosciences Beijing filed Critical China University of Geosciences Beijing
Priority to CN201911141926.1A priority Critical patent/CN110886606B/zh
Publication of CN110886606A publication Critical patent/CN110886606A/zh
Application granted granted Critical
Publication of CN110886606B publication Critical patent/CN110886606B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • E—FIXED CONSTRUCTIONS
    • E21—EARTH DRILLING; MINING
    • E21B—EARTH DRILLING, e.g. DEEP 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

Abstract

本发明涉及一种随钻特征量辅助的惯性测斜方法及装置，属于定向钻井领域，解决随钻测量精度问题；方法包括根据随钻测量的惯性数据，测量下井探管运行的钻探位置、钻探速度和井孔姿态信息；通过外部测量的钻杆钻进速度信息和准零速信息获取速度参考信息，三轴地磁数据获取井孔姿态参考信息；将所述钻探速度与所述速度参考信息之差作为速度观测量，井孔姿态信息与井孔姿态参考信息之差作为姿态观测量，输入卡尔曼滤波器进行数据融合和最优估计，输出姿态误差对姿态矩阵进行实时修正；根据修正后的姿态矩阵实时计算井眼轨迹姿态信息。本发明有效解决了测斜误差问题，提高了系统的精度。

Description

一种随钻特征量辅助的惯性测斜方法及装置
技术领域
本发明涉及定向钻井领域，尤其是一种随钻特征量辅助的惯性测斜方法及装置。
背景技术
随钻测量技术(MWD-Measurement while drilling)可在钻井过程中测量所钻井眼的轨迹，表达井眼在地下的空间分布情况。其中，测量井眼轨迹的方法称为测斜方法。
在现有技术中，测斜方法大多采用独立的测斜方法；其中，磁通门式测斜仪器具有结构简单、价格低、性能稳定的优点，但无法实现对有磁性干扰的井孔进行的井眼轨迹测量；机械陀螺与加速度组合的惯性测斜仪器，由于结构复杂、易损坏、抗振性差、漂移大、必须不断用其他信息加以修正等缺陷，无法应用于随钻测斜。组合测斜的方法在钻井工程实践中应用不多，且大多采用事后测井方法，并不适用于随钻测量；并且随钻测量井下仪器并不是完全静止，即使在停钻时刻，由于井下钻井液的流动，保护井壁稳定性，井下钻具依然处于缓慢的蠕动状态，现有的零速修正方法也不适用于随钻测量。
发明内容
鉴于上述的分析，本发明旨在提供一种随钻特征量辅助的惯性测斜方法及装置，解决在长时间钻井作业过程中，随钻对井眼轨迹进行高精度测量的问题。
本发明的目的主要是通过以下技术方案实现的：
本发明公开了一种随钻特征量辅助的惯性测斜方法，包括以下步骤：
根据随钻测量的惯性数据，计算下井探管运行的惯性测斜参数；所述惯性测斜参数包括钻探位置、钻探速度和井孔姿态信息；
通过外部测量获取速度参考信息和井孔姿态参考信息，所述速度参考信息包括钻杆钻进速度信息和零速修正信息；所述井孔姿态参考信息通过随钻测量的惯性数据计算获得；
将所述钻探速度与所述速度参考信息之差作为速度观测量，井孔姿态信息与井孔姿态参考信息之差作为姿态观测量，输入卡尔曼滤波器进行数据融合和最优估计，输出姿态误差对姿态矩阵进行实时修正；
根据修正后的姿态矩阵实时计算井眼轨迹姿态信息。
进一步地，所述卡尔曼滤波器包括15维的状态变量；其中，
δθx,δθy为真实地理坐标系相对于计算地理坐标系之间的误差角；
δh为真实地理坐标系中的深度误差；
Figure BDA0002281169160000021
为测斜计算坐标系中的速度误差；
ψx,ψy,ψz为载体坐标系相对计算地理坐标系之间的误差角；
Figure BDA0002281169160000022
为载体坐标系中的陀螺零偏；
Figure BDA0002281169160000023
为载体坐标系中的加速度零偏。
进一步地，零速修正信息为根据井下钻具运动空间约束条件，并考虑实际测量中的振动干扰计算的准零速信息；
所述空间约束条件为，在载体坐标系下钻具沿着轴向Y运动，而与Y垂直的横截面上X、Z轴方向的运动速度在钻具载体系中为0。
进一步地，将所述振动干扰等效为白噪声；
所述准零速信息为载体坐标系下的钻具轴向X、Z轴的准零速信息：
Figure BDA0002281169160000024
式中，υx、υz为载体坐标系下，钻具X、Z轴向由于井下振动产生的白噪声；
Figure BDA0002281169160000031
为载体坐标系下，钻具X、Z轴向叠加了振动白噪声的准零速。
进一步地，在载体坐标系下，外部测量获取的速度参考信息
Figure BDA0002281169160000032
式中；ΔL为钻杆长度增量；t为钻进时间；υy为钻杆钻进速度的测量噪声。
进一步地，所述卡尔曼滤波器中根据速度观测量建立的速度量测方程为：
Figure BDA0002281169160000033
式中，Vi c、
Figure BDA0002281169160000034
为测斜计算坐标系下的钻探速度和速度参考信息，HV为量测矩阵，VV为速度量测噪声；
其中，
Figure BDA0002281169160000035
式中，I为3×3单位矩阵，δVi 1为测斜计算坐标系中速度误差；其中，δθi×为真实地理坐标系相对于计算地理坐标系之间的误差角δθi组成的反对称矩阵；Vi n为真实地理坐标系中的钻探速度；
Figure BDA0002281169160000036
为外部测量获取的速度参考信息；
测斜计算坐标系中的姿态矩阵
Figure BDA0002281169160000037
式中，I为3×3单位矩阵，φi×为载体坐标系相对于当地真实地理坐标系的误差角φi的反对称矩阵；
Figure BDA0002281169160000038
为真实地理坐标系中的姿态矩阵；
量测矩阵HV＝[03×3 I3×3 H1 03×6]；其中
Figure BDA0002281169160000039
Figure BDA0002281169160000045
速度量测噪声VV＝[υx，υy，υz]T。
进一步地，所述卡尔曼滤波器中根据姿态观测量建立的姿态量测方程为：ZA＝HAX+VA；
其中，ZA为姿态观测量；
Figure BDA0002281169160000041
VA为姿态角测量的量测噪声。
进一步地，所述卡尔曼滤波器的状态转移矩阵，
Figure BDA0002281169160000042
其中：
Figure BDA0002281169160000043
Figure BDA0002281169160000044
Figure BDA0002281169160000051
Figure BDA0002281169160000052
Figure BDA0002281169160000053
Figure BDA0002281169160000054
为测斜计算坐标系中各轴向速度；g为重力加速度；R为地球半径；h为井眼轨迹的垂深，
Figure BDA0002281169160000055
为地球自转角速率；
Figure BDA0002281169160000056
为地理坐标系相对于地球坐标系的角速度；
Figure BDA0002281169160000057
为两角速度矢量和在测斜计算坐标系中各轴方向上的分量，i＝x,y,z；
Figure BDA0002281169160000058
分别为加速度计在测斜计算坐标系中各轴向测量值；1/βεx、1/βεy、1/βεz分别为陀螺各轴向随机过程的相关时间；1/β▽x、1/β▽y、1/β▽z分别为加速度计各轴向随机噪声的相关时间；测斜计算坐标系中的姿态矩阵
Figure BDA0002281169160000059
进一步地，根据修正后的姿态矩阵，获得井眼轨迹姿态信息，采用校正平均角法计算井眼轨迹垂深位移、水平位移、东向位移、北向位移。
本发明还公开了一种随钻特征量辅助的惯性测斜装置，包括MEMS-IMU、三轴MEMS磁通门、钻杆长度计算器和数据处理电路；
所述MEMS-IMU包括三轴MEMS陀螺和三轴MEMS加速度计，用于随钻测量三轴陀螺数据和三轴加速度数据；
所述三轴MEMS磁通门，用于随钻测量三轴地磁数据；
钻杆长度计算器，用于测量在单位时间内钻杆的长度增量；
数据处理电路接收所述MEMS-IMU、三轴MEMS磁通门和钻杆长度计算器的测量数据，用于执行如上所述的随钻特征量辅助的惯性测斜方法，利用卡尔曼滤波修正后的姿态矩阵实时计算井眼轨迹姿态信息。
本发明有益效果如下：
1、将MEMS-IMU、磁通门和钻杆长度信息进行融合，利用卡尔曼滤波技术对井斜角误差、方位角误差、工具面角误差等井眼轨迹参数进行最优估计，用来补偿MEMS-IMU测斜系统随时间快速发散的高度通道以及随积分计算不断积累的速度参数、位置参数，修正姿态矩阵，并根据修正后的姿态矩阵实时计算井眼轨迹姿态信息，从而提高随钻测斜的测量精度；
2、针对井下没有GPS等外部信息，利用井下钻具的运动特征，充分利用钻杆长度信息，井下钻具运动空间约束带来的准零速等作为组合测斜算法的量测信息，从而实现组合测斜；
3、在磁通门和钻杆信息不理想或发生故障情况下，MEMS-IMU可以暂时独立的提供井眼轨迹信息，并用MEMS-IMU数据完成自我修正。因而，提高组合测斜系统的可靠性；
4、在满足同样的精度要求的情况下，加入磁通门、钻杆长度信息、准零速降低了测斜算法对MEMS-IMU的精度要求，大大降低系统成本，并且能够实现长时间高动态钻井作业环境下的实时、高精度测斜功能。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图1为本发明实施例一中的三个坐标系关系示意图；
图2为本发明实施例一中的惯性测斜方法流程图；
图3为本发明实施例一中的卡尔曼滤波器原理示意图；
图4为本发明实施例一中的随钻测斜装置安装示意图；
图5为本发明实施例一中的随钻测斜装置的测量原理示意图。
具体实施方式
下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理。
实施例一、
首先，对本实施例中涉及到的三个坐标系进行说明：
三个坐标系分别为：
真实地理坐标系n，即当地真实的东北天地理坐标系，也是导航坐标系；
测斜计算坐标系c，是航迹推算获得的计算地理坐标系；
载体坐标系b，是测斜系统载体的右前上坐标系。三个坐标系关系如图1所示：
将东北天地理坐标系oxnynzn作为捷联惯性测斜仪的导航坐标系时，由于测斜系统存在误差，捷联惯性测斜仪计算的经纬度并不等于载体所在实际位置的经纬度。因此为了分析惯性测斜的误差，引入测斜计算坐标系oxcyczc，作为计算机认为的当地地理坐标系，和真实地理坐标系之间有小角度位置误差δθi(i＝x，y，z)。载体坐标系oxbybzb与真实地理坐标系oxnynzn之间的角度误差为
Figure BDA0002281169160000081
也是井眼轨迹的姿态角
Figure BDA0002281169160000082
测斜计算坐标系与真实地理坐标系之间的角度关系为ψi(i＝x，y，z)，三个角度之间的关系为ψi＝φi-δθi。
本实施例公开了一种随钻特征量辅助的惯性测斜方法，如图2所示，包括以下步骤：
步骤S1、根据随钻测量的包括三轴陀螺数据和三轴加速度数据的惯性数据，进行惯导机械编排，计算下井探管运行的惯性测斜参数，所述惯性测斜参数包括钻探位置、速度和井孔姿态信息；
具体的，所述惯性数据由随钻测量的MEMS-IMU传感器输出，MEMS-IMU传感器依靠三轴MEMS陀螺测量角速度，三轴MEMS加速度计测量加速度；经过惯导机械编排，计算下井探管运行的惯性测斜参数，所述惯性测斜参数包括钻探位置、速度和井孔姿态信息等测斜参数。
由于，MEMS-IMU传感器存在累积误差，必须不断用其他信息加以修正，因此，需要引入外部测量的参考信息。
步骤S2、通过外部测量获取速度参考信息和井孔姿态参考信息，所述速度参考信息包括钻杆钻进速度信息和零速修正信息；所述井孔姿态参考信息通过随钻测量的三轴地磁数据和三轴加速度数据计算获得；
具体的，钻杆钻进速度信息为钻井作业时，钻杆向下的钻进速度即为测斜仪载体坐标系沿轴向Y的速度
Figure BDA0002281169160000083
通过由钻杆长度计算器可测量在时间t内的钻杆钻进的长度增量ΔL；
则，钻杆钻进速度可以由井上钻杆钻进的长度增量ΔL除以时间t获得
Figure BDA0002281169160000084
式中；ΔL为钻杆长度增量；t为钻进时间；υy为钻杆钻进速度的测量噪声。
本实施例适用于井下随钻测量系统，钻具在井下的空间约束条件为，在载体坐标系下钻具沿着轴向Y运动，而与Y垂直的横截面上X、Z方向的运动速度在钻具载体系中为0。但在实际测量中，由于有振动干扰存在，X轴和Z轴的速度不是绝对为零；因此，本实施例在零速修正中采用的零速信息为根据井下钻具运动空间约束条件，并考虑实际测量中的振动干扰计算的准零速信息。
具体的，将所述振动干扰等效为白噪声；则准零速信息为载体坐标系下的钻具轴向X、Z轴的准零速信息
Figure BDA0002281169160000091
式中，υx、υz为载体坐标系下，钻具轴向X、Z轴由于井下振动产生的白噪声；
Figure BDA0002281169160000092
Figure BDA0002281169160000093
为钻具轴向X、Z轴叠加了振动白噪声的准零速。
因此，我们可以获得两种外部速度信息作为速度参考信息：准零速信息以及钻杆钻进速度信息，但是都是仪器载体系中的速度信息，由此，在载体坐标系下，外部测量获取的速度参考信息：
Figure BDA0002281169160000094
式中；ΔL为钻杆长度增量；t为钻进时间；υy为钻杆钻进速度的测量噪声，υx、υz振动产生的白噪声。
具体的，随钻测量的三轴地磁数据由三轴MEMS磁通门测量；
三轴MEMS磁通门测量井底钻具的井斜方位角时，大多采用静态测量方案，由三轴加速度数据结合三轴地磁数据计算钻探井孔姿态的参考信息井斜角α磁、方位角A磁、工具面角γ磁的公式为：
Figure BDA0002281169160000095
Figure BDA0002281169160000096
Figure BDA0002281169160000101
式中，
Figure BDA0002281169160000102
为三轴磁通门在载体坐标系的X、Y、Z轴向测量值，
Figure BDA0002281169160000103
分别为加速度计在载体坐标系的X、Y、Z轴向测量值，g为重力加速度。则有磁通门测量的姿态信息如下：
Figure BDA0002281169160000104
步骤S3、将所述钻探速度与所述速度参考信息之差作为速度观测量，井孔姿态信息与井孔姿态参考信息之差作为姿态观测量，输入卡尔曼滤波器进行数据融合和最优估计，输出姿态误差对姿态矩阵进行实时修正。
具体的包括以下子步骤：
步骤S3-1、建立卡尔曼滤波器的滤波方程
离散线性系统的卡尔曼滤波状态方程和观测方程可表示为
Figure BDA0002281169160000105
Z(t)＝H(t)X(t)+V(t) (5)
式中：X(t)为状态量，Z(t)为量测量，F(t)为状态转移矩阵，H(t)为观测矩阵，G(t)过程噪声转移矩阵，W(t)为过程噪声，V(t)为观测噪声。
进而将状态方程(6)和量测方程(7)离散化可得：
Xk＝Φk,k-1Xk-1+Γk-1Wk-1 (6)
Zk＝HkXk+Vk (7)
式中，Xk为k时刻的n维状态向量(被估计量)，Zk为k时刻的m维量测向量，Φk,k-1为k-1到k时刻的系统一步状态转移矩阵(n×n阶)，Hk为k时刻系统量测矩阵(m×n阶)，Γk-1为系统噪声矩阵(n×r阶)，Wk-1为k-1时刻的系统噪声(r维)，Vk为k时刻m维量测噪声。
状态预测估计方程为：
Figure BDA0002281169160000111
式中，
Figure BDA0002281169160000112
为k-1时刻Xk-1的估计值，
Figure BDA0002281169160000113
为从k-1时刻到k时刻的预测值。
方差预测方程为：
Figure BDA0002281169160000114
式中，Pk-1为估计的协方差矩阵。Qk-1为系统噪声的方差矩阵。
状态预测估计方程为：
Figure BDA0002281169160000115
方差迭代方程：
Figure BDA0002281169160000116
式中，Kk为滤波增益，Rk为量测噪声的方差矩阵。
滤波增益方程为：
Figure BDA0002281169160000117
初始条件为：
Figure BDA0002281169160000118
验前统计量为：
E[Wk]＝0,Cov[Wk,Wj]＝E[WkWj T]＝Qkδkj
E[Vk]＝0,Cov[Vk,Vj]＝E[VkVj T]＝Rkδkj
Cov[Wk,Vj]＝E[WkVj T]＝0
Figure BDA0002281169160000121
如图3所示，可以看出，卡尔曼滤波器包含两次更新回路与两个滤波回路：时间更新与量测更新，滤波计算回路与增益计算回路。卡尔曼滤波的基本方程：状态一步预测方程、方差预测方程、状态估计方程、估计误差方差方程以及滤波增益方程。卡尔曼滤波器通过迭代运算，实现预测与校正更新估计。
图3中符号说明如下：
Figure BDA0002281169160000122
为tk-1时刻的状态估计值；
Figure BDA0002281169160000123
是状态
Figure BDA0002281169160000124
的卡尔曼滤波估计；Φk,k-1为tk-1时刻至tk时刻的一步转移阵；Γk-1为系统噪声驱动阵；Hk为量测阵；Rk为量测噪声方差；Qk为系统噪声方差阵；Kk为滤波增益；Zk量测值；Pk-1为估计均方差阵；Pk,k-1一步预测误差方差阵。
步骤S3-2、建立MEMS-IMU的误差模型；
a.MEMS加速度计误差方程
Figure BDA0002281169160000125
MEMS加速度计零偏的误差模型可用一阶马尔科夫过程模型方程表示：
Figure BDA0002281169160000126
上述方程中，
Figure BDA0002281169160000127
分别为加速度计X、Y、Z轴向理想值；
Figure BDA0002281169160000131
分别为加速度计X、Y、Z轴向在载体坐标系下的零偏；
Figure BDA0002281169160000132
分别为加速度计X、Y、Z轴向测量值；KAx、KAy、KAz分别为加速度计X、Y、Z轴向刻度因子；δKAx、δKAy、δKAz分别为加速度计X、Y、Z轴向刻度因子误差；Aij为加速度计j轴向偏向i轴向的失准角；W▽x、W▽y、W▽z分别为加速度计X、Y、Z轴向随机噪声；1/β▽x、1/β▽y、1/β▽z分别为加速度计X、Y、Z轴向随机噪声的相关时间；
a.MEMS陀螺误差模型
Figure BDA0002281169160000133
MEMS陀螺零偏的误差模型可用一阶马尔科夫过程模型方程表示：
Figure BDA0002281169160000134
上述方程中，
Figure BDA0002281169160000135
分别为X、Y、Z轴向陀螺理想值；
Figure BDA0002281169160000136
分别为陀螺X、Y、Z轴向在载体坐标系下的零偏；
Figure BDA0002281169160000137
Figure BDA0002281169160000138
分别为X、Y、Z轴向陀螺测量值；Kgx、Kgy、Kgz分别为陀螺X、Y、Z轴向刻度因子；δKgx、δKgy、δKgz分别为陀螺X、Y、Z轴向刻度因子误差；Mij为陀螺j轴向偏向i轴向的失准角；Wεx、Wεy、Wεz分别为陀螺X、Y、Z轴向随机噪声；1/βεx、1/βεy、1/βεz分别为陀螺X、Y、Z轴向随机过程的相关时间。
步骤S3-3、建立速度、位置、姿态参数的误差方程。
a.姿态误差方程
Figure BDA0002281169160000141
式中：
Figure BDA0002281169160000142
其中，
ψi为载体坐标系相对于测斜计算坐标系之间的误差角；
δθi为真实地理坐标系相对于测斜计算坐标系之间的误差角；
φi为载体坐标系相对于真实地理坐标系的误差角；即φx为井斜角α的误差；φy为工具面角γ的误差；φz为方位角A的误差；
Figure BDA0002281169160000143
分别为各轴向陀螺在计算坐标系下的零偏；
Figure BDA0002281169160000144
为地球自转角速率；
Figure BDA0002281169160000145
为地理坐标系相对于地球坐标系的角速度；
Figure BDA0002281169160000146
为两角速度矢量和在计算坐标系x，y，z方向上的分量。
Figure BDA0002281169160000147
为载体系与计算坐标系之间的姿态转移矩阵，为3×3姿态矩阵，
Figure BDA0002281169160000148
代表矩阵中第j行、第k列元素。
b.位置误差方程
Figure BDA0002281169160000149
其中δθi(i＝x、y)为当地真实地理坐标系相对于计算地理坐标系之间的误差角，即为角位置误差；h为井眼轨迹的垂深；δh为测斜计算坐标系中的垂深误差；R为地球半径；
Vi c(i＝x、y、z)为测斜计算坐标系中各轴向的速度；
δVi 1(i＝x、y、z)为测斜计算坐标系中个轴向的速度误差；
速度误差满足一下关系式：
Figure BDA0002281169160000151
其中，Vi n为真实地理坐标系中各轴的速度值；I为3×3的单位矩阵；
Figure BDA0002281169160000152
为δθi(i＝x、y、z)组成的反对称矩阵。
c.速度误差方程
Figure BDA0002281169160000153
式中：
Figure BDA0002281169160000154
其中，
Vi c(i＝x、y、z)为测斜计算坐标系中各轴向速度；
δVi 1(i＝x、y、z)为测斜计算坐标系中各轴向速度误差；
Figure BDA0002281169160000155
分别为各轴向加速度计在计算坐标系下的零偏；
fi c(i＝x、y、z)分别为各轴向加速度计在计算坐标系下的测量值；
g为重力加速度。
步骤S3-4、确定状态变量以及状态方程
根据卡尔曼滤波状态方程，组合测斜算法的状态方程可以表达为，重写公式
Figure BDA0002281169160000161
a.确定状态方程中的状态量。
由步骤S3-2中的器件误差方程和步骤S3-3中的速度误差方程、位置误差方程、姿态误差方程，确定组合测斜的状态量为：
Figure BDA0002281169160000162
其中，
δθx,δθy为真实地理坐标系相对于计算地理坐标系之间的误差角；
δh为真实地理坐标系中的深度误差；
Figure BDA0002281169160000163
为测斜计算坐标系中的速度误差；
ψx,ψy,ψz为载体坐标系相对计算地理坐标系之间的误差角；
Figure BDA0002281169160000164
为载体坐标系中的陀螺零偏；
Figure BDA0002281169160000165
为载体坐标系中的加速度零偏。
W为系统噪声，根据步骤S3-2中MEMS-IMU器件的误差模型，可以确定为
W＝[0,0,0,0,0,0,0,0,0,Wεx,Wεy,Wεz,W▽x,W▽y,W▽z]T (24)
G为系统噪声转移矩阵，G＝I15×15。
b.确定状态方程中的状态转移矩阵
根据陀螺误差模型、加速度计误差模型、姿态误差方程、位置误差方程、速度误差方程可以推导得出误差转移矩阵F，具体表示为：
Figure BDA0002281169160000166
其中：
Figure BDA0002281169160000171
Figure BDA0002281169160000172
Figure BDA0002281169160000173
Figure BDA0002281169160000174
Figure BDA0002281169160000175
Figure BDA0002281169160000176
为测斜计算坐标系中各轴向速度；g为重力加速度；R为地球半径；h为井眼轨迹的垂深，
Figure BDA0002281169160000177
为地球自转角速率；
Figure BDA0002281169160000178
为地理坐标系相对于地球坐标系的角速度；
Figure BDA0002281169160000181
i＝x,y,z为两角速度矢量和在测斜计算坐标系中各轴方向上的分量；
Figure BDA0002281169160000182
分别为加速度计在测斜计算坐标系中各轴向测量值；1/βεx、1/βεy、1/βεz分别为陀螺各轴向随机过程的相关时间；1/β▽x、1/β▽y、1/β▽z分别为加速度计各轴向随机噪声的相关时间；测斜计算坐标系中的姿态矩阵
Figure BDA0002281169160000183
步骤S3-5、建立量测量以及量测方程
重写卡尔曼滤波的量测方程Z(t)＝H(t)X(t)+V(t)，确定量测信息Z，量测转移矩阵H，以及量测噪声V。
a.在卡尔曼滤波器中根据速度观测量建立的速度量测方程
所述卡尔曼滤波器中根据速度观测量建立的速度量测方程为：
Figure BDA0002281169160000184
式中，Vi c、
Figure BDA0002281169160000185
为测斜计算坐标系中下的钻探速度和速度参考信息，HV为量测矩阵，VV为速度量测噪声；
其中，
Figure BDA0002281169160000186
式中，I为3×3单位矩阵，δVi1为测斜计算坐标系中速度误差；其中，δθi×为真实地理坐标系相对于计算地理坐标系之间的误差角δθi组成的反对称矩阵；Vi n为真实地理坐标系中的根据惯导数据计算的钻探速度；
Figure BDA0002281169160000187
为外部测量获取的速度参考信息；
测斜计算坐标系中的姿态矩阵
Figure BDA0002281169160000188
式中，I为3×3单位矩阵，φi×为载体坐标系相对于当地真实地理坐标系的误差角φi的反对称矩阵；
Figure BDA0002281169160000191
为真实地理坐标系中的姿态矩阵；
由速度观测量
Figure BDA0002281169160000192
推到可得：
Figure BDA0002281169160000193
由于
Figure BDA0002281169160000194
ψi×＝(φi×)-(δθi×)所以
Figure BDA0002281169160000195
则捷联惯导测斜系统和外部速度信息形成的量测矩阵：
HV＝[03×3 I3×3 H1 03×6]；其中
Figure BDA0002281169160000196
Figure BDA0002281169160000197
速度量测噪声VV＝[υx，υy，υz]T。
b.在卡尔曼滤波器中根据姿态观测量建立的姿态量测方程
所述卡尔曼滤波器中根据姿态观测量建立的姿态量测方程为：
Figure BDA0002281169160000198
其中，ZA为姿态观测量；
Figure BDA0002281169160000201
VA为姿态角测量的量测噪声。
综上所述，
由随钻特征量辅助的惯性测斜系统的量测方程，重写量测方程为：
Z(t)＝H(t)X(t)+V(t)；
Figure BDA0002281169160000202
其中：Z为量测量，H为量测矩阵，V为量测噪声。
步骤S3-6、将建立的卡尔曼滤波连续系统进行离散化；
具体的，根据
Figure BDA0002281169160000203
Figure BDA0002281169160000204
将连续系统离散化：其中M1＝Q(t),Mi+1＝FMi+Mi TFT，在离散化中的参数t为时间，Δt为时间增量。
步骤S3-7、将观测到的量测信息输入卡尔曼滤波器进行数据融合，对状态变量进行最优估计，根据估计的状态变量对所述惯性数据及惯性5测斜参数进行在线误差补偿。
具体包括，
用卡尔曼滤波估计的陀螺和加速度计沿载体坐标系三轴方向上的等效漂移和零位进行测量数据补偿；
对加速度计数据补偿：
Figure BDA0002281169160000205
对陀螺数据补偿：
Figure BDA0002281169160000211
用卡尔曼滤波估计的速度误差对速度及时修正，得到精度较高的速度信息：
Vi n＝[I+(δθi×)](Vi c-δVi 1) (29)
用卡尔曼滤波估计的姿态误差对计算地理坐标系中的姿态矩阵
Figure BDA0002281169160000212
及时修正补偿，获得真实地理坐标系中的姿态矩阵
Figure BDA0002281169160000213
Figure BDA0002281169160000214
Figure BDA0002281169160000215
为3×3姿态矩阵，Tjk代表矩阵中第j行、第k列元素。
用卡尔曼滤波估计的位置误差对测斜计算坐标系中的位置矩阵及时修正补偿：
Figure BDA0002281169160000216
Figure BDA0002281169160000217
为测斜计算坐标系中的位置矩阵，
Figure BDA0002281169160000218
为真实地理坐标系的中的位置矩阵。
步骤S4、根据修正后的姿态矩阵实时计算井眼轨迹姿态信息。
具体的，根据姿态矩阵
Figure BDA0002281169160000219
计算井眼轨迹姿态信息包括：
倾斜角α：α＝sin-1(T23)(32)
工具面角γ：
Figure BDA00022811691600002110
方位角A：
Figure BDA00022811691600002111
优选的，根据校正平均角法获得井眼轨迹信息，垂深增量ΔD、水平位移增量ΔS、东向位移增量ΔE、北向位移增量ΔN。
Figure BDA0002281169160000221
其中Δα＝αk+1-αk，
Figure BDA0002281169160000222
k代表第k个测点，k+1为第k+1个测点，Δα为第k+1与第k个测点之间的井斜角之差值，
Figure BDA0002281169160000223
为第k+1与第k个测点之间的方位角之差值，αc为第k+1与第k个测点之间的井斜角平均值，
Figure BDA0002281169160000224
为第k+1与第k个测点之间的方位角平均值。
则井眼轨迹信息垂深D，水平位移S，东向位移E，北向位移N由以下公式可得：
Figure BDA0002281169160000225
相比于现有技术，本实施例将MEMS-IMU、磁通门和钻杆长度信息进行融合，利用卡尔曼滤波技术，建立完善的系统误差模型，对井斜角误差、方位角误差、工具面角误差等井眼轨迹参数进行最优估计，用来补偿MEMS-IMU测斜系统随时间快速发散的高度通道以及随积分计算不断积累的速度参数、位置参数，修正姿态矩阵，并根据修正后的姿态矩阵实时计算井眼轨迹姿态信息，从而提高随钻测斜的测量精度。
实施例二、
本实施例公开了一种基于多传感器数据融合的随钻测斜装置，包括MEMS-IMU、三轴MEMS磁通门、钻杆长度计算器和数据处理电路；
所述MEMS-IMU包括三轴MEMS陀螺和三轴MEMS加速度计，用于随钻测量三轴陀螺数据和三轴加速度数据；
所述三轴MEMS磁通门，用于随钻测量三轴地磁数据；
钻杆长度计算器，用于测量在单位时间内钻杆的长度增量；
数据处理电路接收所述MEMS-IMU、三轴MEMS磁通门和钻杆长度计算器的测量数据，用于执行实施例一中所述的随钻特征量辅助的惯性测斜方法，利用卡尔曼滤波修正后的姿态矩阵，采用校正平均角法实时计算井眼轨迹位置信息。具体的随钻测斜装置安装示意图如图4所示，随钻测斜装置的测量原理示意图如图5所示；
优选的，在磁通门和钻杆信息不理想或发生故障情况下，MEMS-IMU可以暂时独立的提供井眼轨迹信息，并用MEMS-IMU数据完成自我修正。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (3)

1.一种随钻特征量辅助的惯性测斜方法，其特征在于，包括以下步骤：
根据随钻测量的惯性数据，计算下井探管运行的惯性测斜参数；所述惯性测斜参数包括钻探位置、钻探速度和井孔姿态信息；
通过外部测量获取速度参考信息和井孔姿态参考信息，所述速度参考信息包括钻杆钻进速度信息和零速修正信息；所述井孔姿态参考信息通过随钻测量的惯性数据计算获得；
将所述钻探速度与所述速度参考信息之差作为速度观测量，井孔姿态信息与井孔姿态参考信息之差作为姿态观测量，输入卡尔曼滤波器进行数据融合和最优估计，输出姿态误差对姿态矩阵进行实时修正；
根据修正后的姿态矩阵实时计算井眼轨迹姿态信息；
所述卡尔曼滤波器中根据速度观测量建立的速度量测方程为：
Figure FDA0003184250280000011
式中，Vi c、
Figure FDA0003184250280000012
为测斜计算坐标系下的钻探速度和速度参考信息；HV为量测矩阵，VV为速度量测噪声；其中，引入的测斜计算坐标系oxcyczc为计算机认为的当地地理坐标系，用于分析惯性测斜的误差，与真实地理坐标系之间有小角度位置误差δθi，i代表坐标系的三个轴向x，y，z；
其中，
Figure FDA0003184250280000013
式中，I为3×3单位矩阵，δVi 1为测斜计算坐标系中速度误差；其中，δθi×为真实地理坐标系相对于计算地理坐标系之间的误差角δθi组成的反对称矩阵；Vi n为真实地理坐标系中的钻探速度；
Figure FDA0003184250280000014
为外部测量获取的载体坐标系下的速度参考信息；
测斜计算坐标系中的姿态矩阵
Figure FDA0003184250280000021
式中，I为3×3单位矩阵，φi×为载体坐标系相对于当地真实地理坐标系的误差角φi的反对称矩阵；
Figure FDA0003184250280000022
为真实地理坐标系中的姿态矩阵；
量测矩阵HV＝[03×3 I3×3 H1 03×6]；其中
Figure FDA0003184250280000023
Figure FDA0003184250280000024
速度量测噪声VV＝[υx，υy，υz]T；υx、υz为载体坐标系下，钻具X、Z轴向由于井下振动产生的白噪声，υy为钻杆钻进速度的测量噪声；
X为卡尔曼滤波器的15维的状态变量；包括：
δθx,δθy为真实地理坐标系相对于计算地理坐标系之间的误差角；
δh为真实地理坐标系中的深度误差；
Figure FDA0003184250280000025
为测斜计算坐标系中的速度误差；
ψx,ψy,ψz为载体坐标系相对计算地理坐标系之间的误差角；
Figure FDA0003184250280000026
为载体坐标系中的陀螺零偏；
Figure FDA0003184250280000027
为载体坐标系中的加速度零偏；
零速修正信息为根据井下钻具运动空间约束条件，并考虑实际测量中的振动干扰计算的准零速信息；
所述空间约束条件为，在载体坐标系下钻具沿着轴向Y运动，而与Y垂直的横截面上X、Z轴方向的运动速度在钻具载体系中为0；
将所述振动干扰等效为白噪声；
所述准零速信息为载体坐标系下的钻具轴向X、Z轴的准零速信息：
Figure FDA0003184250280000031
式中，υx、υz为载体坐标系下，钻具X、Z轴向由于井下振动产生的白噪声；
Figure FDA0003184250280000032
为载体坐标系下，钻具X、Z轴向叠加了振动白噪声的准零速；
在载体坐标系下，外部测量获取的速度参考信息
Figure FDA0003184250280000033
式中；ΔL为钻杆长度增量；t为钻进时间；υy为钻杆钻进速度的测量噪声；
所述卡尔曼滤波器中根据姿态观测量建立的姿态量测方程为：ZA＝HAX+VA；
其中，ZA为姿态观测量；
Figure FDA0003184250280000034
VA为姿态角测量的量测噪声；
所述卡尔曼滤波器的状态转移矩阵，
Figure FDA0003184250280000035
其中：
Figure FDA0003184250280000036
Figure FDA0003184250280000041
Figure FDA0003184250280000042
Figure FDA0003184250280000043
Figure FDA0003184250280000044
Figure FDA0003184250280000045
为测斜计算坐标系中各轴向速度；g为重力加速度；R为地球半径；h为井眼轨迹的垂深，
Figure FDA0003184250280000046
为地球自转角速率；
Figure FDA0003184250280000047
为地理坐标系相对于地球坐标系的角速度；
Figure FDA0003184250280000048
为两角速度矢量和在测斜计算坐标系中各轴方向上的分量，i＝x,y,z；
Figure FDA0003184250280000049
分别为加速度计在测斜计算坐标系中各轴向测量值；1/βεx、1/βεy、1/βεz分别为陀螺各轴向随机过程的相关时间；1/β▽x、1/β▽y、1/β▽z分别为加速度计各轴向随机噪声的相关时间；测斜计算坐标系中的姿态矩阵
Figure FDA0003184250280000051
2.根据权利要求1所述的惯性测斜方法，其特征在于，根据修正后的姿态矩阵，获得井眼轨迹姿态信息，采用校正平均角法计算井眼轨迹垂深位移、水平位移、东向位移、北向位移。
3.一种随钻特征量辅助的惯性测斜装置，其特征在于，包括MEMS-IMU、三轴MEMS磁通门、钻杆长度计算器和数据处理电路；
所述MEMS-IMU包括三轴MEMS陀螺和三轴MEMS加速度计，用于随钻测量三轴陀螺数据和三轴加速度数据；
所述三轴MEMS磁通门，用于随钻测量三轴地磁数据；
钻杆长度计算器，用于测量在单位时间内钻杆的长度增量；
数据处理电路接收所述MEMS-IMU、三轴MEMS磁通门和钻杆长度计算器的测量数据，用于执行权利要求1-2任一项所述的随钻特征量辅助的惯性测斜方法，利用卡尔曼滤波修正后的姿态矩阵实时计算井眼轨迹姿态信息。
CN201911141926.1A 2019-11-20 2019-11-20 一种随钻特征量辅助的惯性测斜方法及装置 Active CN110886606B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911141926.1A CN110886606B (zh) 2019-11-20 2019-11-20 一种随钻特征量辅助的惯性测斜方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911141926.1A CN110886606B (zh) 2019-11-20 2019-11-20 一种随钻特征量辅助的惯性测斜方法及装置

Publications (2)

Publication Number Publication Date
CN110886606A CN110886606A (zh) 2020-03-17
CN110886606B true CN110886606B (zh) 2021-09-14

Family

ID=69748114

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911141926.1A Active CN110886606B (zh) 2019-11-20 2019-11-20 一种随钻特征量辅助的惯性测斜方法及装置

Country Status (1)

Country Link
CN (1) CN110886606B (zh)

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111380712B (zh) * 2020-04-05 2021-05-04 新疆正通石油天然气股份有限公司 一种弯曲水平井眼中屈曲管柱钻进的评价方法
CN111502635B (zh) * 2020-04-24 2023-03-14 神华神东煤炭集团有限责任公司 一种煤矿井下防爆型陀螺测斜仪及其使用方法
CN111521178B (zh) * 2020-04-28 2021-01-15 中国人民解放军国防科技大学 基于管长约束的钻探用定位定向仪孔内定位方法
CN111504313B (zh) * 2020-04-28 2021-01-05 中国人民解放军国防科技大学 基于速度信息辅助的钻探用定位定向仪孔内定位方法
CN111521177B (zh) * 2020-04-28 2021-01-05 中国人民解放军国防科技大学 管长信息辅助测速的钻探用定位定向仪孔内定位方法
CN111896026B (zh) * 2020-05-11 2021-05-18 中国科学院地质与地球物理研究所 固态谐振陀螺自校准方法及系统
CN111878056B (zh) * 2020-05-11 2021-04-13 中国科学院地质与地球物理研究所 一种陀螺随钻测量系统及方法
CN111623821B (zh) * 2020-05-15 2022-10-04 天津时空经纬测控技术有限公司 隧道钻孔方向的检测、偏差检测、钻孔位置确定的方法
CN111810112B (zh) * 2020-06-18 2021-12-03 中国地质大学（武汉） 基于粒子滤波和模型预测控制的垂直钻进纠偏控制方法
CN112051569B (zh) * 2020-09-10 2024-04-05 北京经纬恒润科技股份有限公司 雷达目标跟踪速度修正方法及装置
CN113738345A (zh) * 2020-11-25 2021-12-03 中国石油天然气集团有限公司 基于常规随钻工具判断钻井情况的方法
CN112729222A (zh) * 2020-12-14 2021-04-30 北京航空航天大学 一种桩挖转杆位置的实时测量方法
CN113153270A (zh) * 2021-04-27 2021-07-23 西南石油大学 一种近钻头动态井斜角与工具面角的随钻测量方法
CN115163030A (zh) * 2022-06-24 2022-10-11 中国石油天然气集团有限公司 一种井的重钻入处理方法、装置和系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2338075A1 (en) * 2001-01-19 2002-07-19 University Technologies International Inc. Continuous measurement-while-drilling surveying
CN106246168A (zh) * 2016-08-29 2016-12-21 中国科学院地质与地球物理研究所 一种近钻头钻具姿态随钻测量装置及测量方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6631563B2 (en) * 1997-02-07 2003-10-14 James Brosnahan Survey apparatus and methods for directional wellbore surveying
US6883747B2 (en) * 2003-03-28 2005-04-26 Northrop Grumman Corporation Projectile guidance with accelerometers and a GPS receiver
US6918186B2 (en) * 2003-08-01 2005-07-19 The Charles Stark Draper Laboratory, Inc. Compact navigation system and method
WO2006065923A2 (en) * 2004-12-14 2006-06-22 Raytheon Utd Centralizer-based survey and navigation device and method
CN100593689C (zh) * 2006-05-26 2010-03-10 南京航空航天大学 基于捷联惯性导航系统的姿态估计和融合的方法
CN100489459C (zh) * 2006-07-17 2009-05-20 北京航空航天大学 适用于全光纤数字测斜仪的捷联惯性组合测量控制装置
US8439130B2 (en) * 2010-02-22 2013-05-14 Schlumberger Technology Corporation Method and apparatus for seismic data acquisition during drilling operations
CN102337883B (zh) * 2011-09-15 2015-09-16 中煤科工集团重庆研究院有限公司 随钻测量系统钻孔深度实时跟踪测量方法
CN103114846B (zh) * 2013-01-25 2016-05-25 北京航空航天大学 一种基于光纤陀螺测斜仪的测斜数据的事后处理系统
CN105300379B (zh) * 2015-10-13 2017-12-12 上海新纪元机器人有限公司 一种基于加速度的卡尔曼滤波姿态估计方法及系统
NL2016859B1 (en) * 2016-05-30 2017-12-11 Engie Electroproject B V A method of and a device for estimating down hole speed and down hole torque of borehole drilling equipment while drilling, borehole equipment and a computer program product.

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2338075A1 (en) * 2001-01-19 2002-07-19 University Technologies International Inc. Continuous measurement-while-drilling surveying
CN106246168A (zh) * 2016-08-29 2016-12-21 中国科学院地质与地球物理研究所 一种近钻头钻具姿态随钻测量装置及测量方法

Also Published As

Publication number Publication date
CN110886606A (zh) 2020-03-17

Similar Documents

Publication Publication Date Title
CN110886606B (zh) 一种随钻特征量辅助的惯性测斜方法及装置
CN110792430B (zh) 一种基于多传感器数据融合的随钻测斜方法及装置
WO2021227012A1 (zh) 一种姿态测量方法
CN110799727B (zh) 用于生成下向井眼惯性测量单元的输出的系统和方法
US20210348924A1 (en) Attitude measurement method
CN104736963B (zh) 测绘系统和方法
US10047600B2 (en) Attitude reference for tieback/overlap processing
US10550686B2 (en) Tumble gyro surveyor
CN104515527B (zh) 一种无gps信号环境下的抗粗差组合导航方法
US20020133958A1 (en) Continuous measurement-while-drilling surveying
US6480119B1 (en) Surveying a subterranean borehole using accelerometers
EP3221665B1 (en) Inertial carousel positioning
US20040089474A1 (en) Continuous measurement-while-drilling surveying
CN105806364A (zh) 一种矿用回转钻机测斜仪探管的校准方法
CN108387205A (zh) 基于多传感器数据融合的钻具姿态测量系统的测量方法
CN107228664A (zh) 矿用陀螺测斜仪捷联惯导系统姿态解算及零速校正方法
Pecht et al. Observability analysis for INS alignment in horizontal drilling
US4768152A (en) Oil well bore hole surveying by kinematic navigation
EP2800870B1 (en) Navigation device and method for surveying and directing a borehole under drilling conditions
CN111141283A (zh) 一种通过地磁数据判断行进方向的方法
Chao et al. An innovative MEMS-based MWD method for directional drilling
Weston et al. The Combination of Solid-State Gyroscopic and Magnetic Surveys Provides Improved Magnetic-Survey Data and Enhanced Survey Quality Control
ElGizawy et al. Continuous wellbore surveying while drilling utilizing MEMS gyroscopes based on Kalman filtering
Ursenbach et al. Effect of in-drilling alignment with general dynamic error model on azimuth estimation
CA2271156C (en) Method of correcting wellbore magnetometer errors

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