-
摘要:
为了解决现有水平定向钻进导向方法测量精度低，定位深度浅，测量数据不连续以及受磁场干扰影响大的问题，提出一种基于磁信标的水平定向钻进导向定位方法，实现水平定向钻进的位置和姿态测量；设计了基于三阶磁梯度张量的测量阵列，通过仿真验证该方法可以消除地磁场的影响，实现位置测量；依据惯性导航理论，利用MEMS传感器、测量阵列位置信息和零速信息进行多源数据融合，对惯性导航解算出的姿态误差进行估计和修正，提高姿态测量精度。通过地面位置测量实验对基于三阶磁梯度张量的位置测量算法进行验证，位置测量结果的平均相对误差为4.57%；通过地面姿态测量实验对多源数据融合算法进行验证，结果表明，倾角误差在0.3°以内，工具面角误差在0.7°以内，方位角的误差值在1°以内，证明了提出的姿态测量算法的准确性和有效性。研究成果可以提高非开挖水平定向钻进适用范围、实时性和定位精度，为水平定向钻进导向技术自动化提供了理论基础。
Abstract:In order to solve the problems of existing horizontal directional drilling (HDD) guidance methods, including low measurement accuracy, shallow positioning distance, discontinuous measurement data and great interference by magnetic field, a magnetic beacon-based horizontal directional drilling guidance and positioning method was proposed to achieve the position and attitude measurement of HDD. Meanwhile, a measurement array based on the third-order magnetic gradient tensor was designed, which verifies through simulation that the method could eliminate the influence of the geomagnetic fields and achieve the position measurement. In order to improve the attitude measurement accuracy, the multisource data fusion was performed with MEMS sensors, measurement array position information and zero velocity information based on the inertial navigation theory to estimate and correct the attitude errors solved by inertial navigation. Besides, the position measurement algorithm based on the third-order magnetic gradient tensor was verified by ground-based position measurement experiments, with the average relative error of 4.57% in the position measurement results. In addition, the multi-source data fusion algorithm was verified by ground attitude measurement experiments. The results show that the error is within 0.3° for the inclination angle, within 0.7° for the tool face angle, and within 1° for the azimuth angle. This indicates the accuracy and effectiveness of the position measurement algorithm based on the third-order magnetic gradient tensor and the attitude measurement algorithm with multi-source data fusion proposed in this paper. The research results could improve the applicable range, real-time performance and positioning accuracy of trenchless horizontal directional drilling, thereby providing the theoretical basis for the automation of HDD guidance technology.
-
非开挖水平定向钻进技术是市政管线铺设、煤矿地质异常体探测、瓦斯抽采防治以及长距离山体水平勘查的有效手段[1]。水平定向钻进施工流程是利用钻机在入口处开始导向钻进，钻进过程中监控钻头姿态和空间位置参数，通过控制手段使钻孔按设计轨迹钻进成孔，再利用扩孔设备实现分级扩孔至设计直径，最后实施管道回拉敷设[2]。实现水平定向钻进精确导向，对提高钻进效率，避免孔壁坍塌，降低钻进成本具有重要意义[3]。
目前，常用的水平定向钻进导向测量方法分为格栅磁场导向方法、手持式导向仪和随钻测量系统三种。刘金祯[4]等人提出的格栅磁场导向方法通过在钻孔轨迹上方的地面布置多个通电线框并标定几何形状，产生固定磁场，利用底部钻具中的信号棒测量磁场强度进行钻头位置计算，该方法可以有效屏蔽环境磁场对信号棒的干扰，测量距离深，但对于河流、沼泽等无法布置通电线框的施工场地，该方法无法应用。手持式导向仪目前是美国Ditch Witch公司的Subsite系统、Digital Control Inc.公司的猎鹰系统较为先进，该方法通过在底部钻具上安装低频电磁场信号发射源，利用手持地面雷达进行信号强度检测并进行位置计算，该方法设备简单，无线测量操作方便，但存在探测深度浅，受环境磁场干扰严重的问题[5]。随钻测量系统是由三轴磁强计和三轴加速度计进行姿态测量，根据钻探进尺递推钻头位置坐标，该方法可以对钻头姿态和钻孔轨迹进行实时测量，操作简单，但存在受环境磁场干扰严重，累计误差大，测量精度不足的问题[6]。针对上述问题，提出了一种基于磁信标的水平定向钻进位置测量和姿态测量方法。
磁信标定位技术，已经在军用和水下定位领域展开应用[7-8]，该技术通过设置通电螺线圈或旋转永磁体发射具有强穿透力的低频磁场信号，利用磁强计对磁信号进行接收并解算，实现对目标载体的定位，具有强鲁棒性、高稳定性以及高抗干扰性的优点[9]。因此，针对水平定向钻进施工特点，采用磁信标定位技术、惯性导航技术和多源数据融合技术[10]，以期实现一种对施工场地条件要求低，受环境磁场干扰小，适用范围广，测量精度高的水平定向钻进导向方法。通过开展磁信标和测量阵列功能与结构设计，对采集数据进行数据融合，以达到水平定向钻进底部钻具姿态和位置的实时高精度测量的目的。为提高非开挖水平定向钻进适用范围和定位精度，水平定向钻进导向技术向自动化和智能化发展提供了基础。
1 基于磁信标测量阵列的位置测量方法
基于磁场梯度张量不变量定位方法，可以实现对目标的单点实时定位[10]，测量的构型采用十字形的较多[11]。但由于地磁场的存在，很难准确测量出地面磁信标在测点处的真实值。因此，为了消除水平定向钻进中地磁场等干扰磁场的影响，提出一种基于磁信标的三阶磁梯度张量目标定位方法，通过差分处理对磁场梯度张量定位方法进行改进[12-13]。
图1中，对测量阵列中的三轴磁强计进行编号，分别为0、1、2、3、4、5，共计6个三轴磁强计，按图1所示排列，下文采用
$ {F}_{0}、{F}_{1}、{F}_{2}、{F}_{3}、{F}_{4}、{F}_{5} $ 代表磁强计，其中$ x、y $ 轴处于同一平面，z轴指向大地，间隔距离为L。分别对$x、y、{\textit{z}}$ 求导。![]()
图 1 基于三阶磁梯度张量的测量阵列Figure 1. Measurement array based on the third-order magnetic gradient tensor$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}_{xxx}}}&{{{\boldsymbol{B}}_{xyx}}}&{{{\boldsymbol{B}}_{x{\textit{z}}x}}} \\ {{{\boldsymbol{B}}_{yxy}}}&{{{\boldsymbol{B}}_{yyy}}}&{{{\boldsymbol{B}}_{y{\textit{z}}y}}} \\ {{{\boldsymbol{B}}_{{\textit{z}}x{\textit{z}}}}}&{{{\boldsymbol{B}}_{{\textit{z}}y{\textit{z}}}}}&{{{\boldsymbol{B}}_{{\textit{zzz}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x \\ y \\ {\textit{z}} \end{array}} \right] = - 4\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}_{xx}}} \\ {{{\boldsymbol{B}}_{yy}}} \\ {{{\boldsymbol{B}}_{{\textit{zz}}}}} \end{array}} \right] $$ (1) 式中：
${{\boldsymbol{B}}_{xi}}$ 、${{\boldsymbol{B}}_{yi}}$ 、${{\boldsymbol{B}}_{{\textit{z}}i}}$ 为$ i $ 号磁强计磁场强度分量值；${{\boldsymbol{B}}_{xxx}}$ 、${{\boldsymbol{B}}_{xyx}}$ 、${{\boldsymbol{B}}_{x\textit{z}x}}$ 、${{\boldsymbol{B}}_{yxy}}$ 、${{\boldsymbol{B}}_{yyy}}$ 、${{\boldsymbol{B}}_{y\textit{z}y}}$ 、${{\boldsymbol{B}}_{\textit{z}x\textit{z}}}$ 、${{\boldsymbol{B}}_{\textit{z}y\textit{z}}}$ 、${{\boldsymbol{B}}_{\textit{zzz}}}$ 为磁强计磁场强度的三阶导数值；${{\boldsymbol{B}}_{xxxi}}$ 、${{\boldsymbol{B}}_{xyxi}}$ 、${{\boldsymbol{B}}_{x{\textit{z}}xi}}$ 、${{\boldsymbol{B}}_{yxyi}}$ 、${{\boldsymbol{B}}_{yyyi}}$ 、${{\boldsymbol{B}}_{y{\textit{z}}yi}}$ 、${{\boldsymbol{B}}_{{\textit{z}}x{\textit{z}}i}}$ 、${{\boldsymbol{B}}_{{\textit{z}}y{\textit{z}}i}}$ 、${{\boldsymbol{B}}_{{\textit{zzz}}i}}$ 为$ i $ 号磁强计磁场强度的三阶导数值；${{\boldsymbol{B}}_{xx}}$ 、${{\boldsymbol{B}}_{yy}}$ 、${{\boldsymbol{B}}_{{\textit{zz}}}}$ 为$ i $ 号磁强计磁梯度张量值，($i = 0,1,2, 3,4,5$ )。由麦克斯韦方程组可知，磁场的散度及旋度为零，张量矩阵具有无迹性与对称性，因此${{\boldsymbol{B}}_{xyx}} = {{\boldsymbol{B}}_{yxx}}$ ，${{\boldsymbol{B}}_{x{\textit{z}}x}} = {{\boldsymbol{B}}_{{\textit{z}}xx}}$ ，${{\boldsymbol{B}}_{yxy}} = {{\boldsymbol{B}}_{xyy}}$ ，${{\boldsymbol{B}}_{y{\textit{z}}y}} = {{\boldsymbol{B}}_{{\textit{z}}yy}}$ 。以磁强计${F}_{0}、{F}_{2}、 {F}_{4}$ 为例，计算${{\boldsymbol{B}}}_{yxy0}、{{\boldsymbol{B}}}_{yyy0}、$ ${{\boldsymbol{B}}}_{y{\textit{z}}y0}、{{\boldsymbol{B}}}_{yy0}$ 。磁强计${F}_{0}、 {F}_{2}、{F}_{4}$ 的$ y $ 轴分量为${{\boldsymbol{B}}}_{y0}、{{\boldsymbol{B}}}_{y2}、{{\boldsymbol{B}}}_{y4}$ 。其中，${{\boldsymbol{B}}_y}_0 = {{\boldsymbol{B}}_y}({y_0})$ ，${{\boldsymbol{B}}_{y2}} = {{\boldsymbol{B}}_y}({y_0} - L)$ ，${{\boldsymbol{B}}_{y4}} = {{\boldsymbol{B}}_y}({y_0} + L)$ 。$ {y_0} $ 为磁强计$ {F_0} $ 的$ y $ 轴位置坐标，在磁偶极子模型中，误差余项较小，忽略此项后有：$$ {{\boldsymbol{B}}_{yyy0}} = {{\boldsymbol{B}}_{yyy}}({y_0}) = \frac{{{{\boldsymbol{B}}_y}({y_0} + L) + {{\boldsymbol{B}}_y}({y_0} - L) - 2{{\boldsymbol{B}}_y}({y_0})}}{{{L^2}}} $$ (2) $$ {{\boldsymbol{B}}_{yy0}} = {{\boldsymbol{B}}_{yy}}({y_0}) = \frac{{{{\boldsymbol{B}}_y}({y_0} + L) - {{\boldsymbol{B}}_y}({y_0} - L)}}{{{L^2}}} $$ (3) 根据式(1)和
${{\boldsymbol{B}}_{xx}} + {{\boldsymbol{B}}_{yy}} + {{\boldsymbol{B}}_{{\textit{zz}}}} = 0$ 可知：$$\left\{\begin{aligned} & {{\boldsymbol{B}}_{xxx0}}x + {{\boldsymbol{B}}_{xyx0}}y + {{\boldsymbol{B}}_{x{\textit{z}}x0}}{\textit{z}} = - 4{{\boldsymbol{B}}_{xx0}} \\ & {{\boldsymbol{B}}_{yxy0}}x + {{\boldsymbol{B}}_{yyy0}}y + {{\boldsymbol{B}}_{y{\textit{z}}y0}}{\textit{z}} = - 4{{\boldsymbol{B}}_{yy0}} \\ & {{\boldsymbol{B}}_{xxx1}}x + {{\boldsymbol{B}}_{xyx1}}y + {{\boldsymbol{B}}_{x{\textit{z}}x1}}{\textit{z}} = - 4{{\boldsymbol{B}}_{xx1}} + {{\boldsymbol{B}}_{xxx1}}L \end{aligned}\right. $$ (4) 由以上各式可得测量阵列测得的位置信息：
$$ \left[ {\begin{array}{*{20}{c}} x \\ y \\ {\textit{z}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}_{xxx0}}}&{{{\boldsymbol{B}}_{xyx0}}}&{{{\boldsymbol{B}}_{x{\textit{z}}x0}}} \\ {{{\boldsymbol{B}}_{yxy0}}}&{{{\boldsymbol{B}}_{yyy0}}}&{{{\boldsymbol{B}}_{y{\textit{z}}y0}}} \\ {{{\boldsymbol{B}}_{{\textit{z}}x{\textit{z}}1}}}&{{{\boldsymbol{B}}_{{\textit{z}}y{\textit{z}}1}}}&{{{\boldsymbol{B}}_{{\textit{zzz}}1}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} { - 4{{\boldsymbol{B}}_{xx0}}} \\ { - 4{{\boldsymbol{B}}_{yy0}}} \\ { - 4{{\boldsymbol{B}}_{xx1}} + {{\boldsymbol{B}}_{xxx1}}L} \end{array}} \right] $$ (5) 其中：
$$ \left\{ \begin{split} & {{{\boldsymbol{B}}_{xxx0}} = \frac{{{{\boldsymbol{B}}_{x1}} + {{\boldsymbol{B}}_{x3}} - 2{{\boldsymbol{B}}_{x0}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{xyx0}} = \frac{{{{\boldsymbol{B}}_{y1}} + {{\boldsymbol{B}}_{y3}} - 2{{\boldsymbol{B}}_{y0}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{x{\textit{z}}x0}} = \frac{{{{\boldsymbol{B}}_{{\textit{z}}1}} + {{\boldsymbol{B}}_{{\textit{z}}3}} - 2{{\boldsymbol{B}}_{{\textit{z}}0}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{yxy0}} = \frac{{{{\boldsymbol{B}}_{x2}} + {{\boldsymbol{B}}_{x4}} - 2{{\boldsymbol{B}}_{x0}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{yyy0}} = \frac{{{{\boldsymbol{B}}_{y2}} + {{\boldsymbol{B}}_{y4}} - 2{{\boldsymbol{B}}_{y0}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{y{\textit{z}}y0}} = \frac{{{{\boldsymbol{B}}_{{\textit{z}}2}} + {{\boldsymbol{B}}_{{\textit{z}}4}} - 2{{\boldsymbol{B}}_{{\textit{z}}0}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{xxx1}} = \frac{{{{\boldsymbol{B}}_{x0}} + {{\boldsymbol{B}}_{x5}} - 2{{\boldsymbol{B}}_{x1}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{xyx1}} = \frac{{{{\boldsymbol{B}}_{y0}} + {{\boldsymbol{B}}_{y5}} - 2{{\boldsymbol{B}}_{y1}}}}{{{L^2}}}} \\ & {{{\boldsymbol{B}}_{x{\textit{z}}x1}} = \frac{{{{\boldsymbol{B}}_{{\textit{z}}5}} + {{\boldsymbol{B}}_{\textit{z}}}_0 - 2{{\boldsymbol{B}}_{{\textit{z}}1}}}}{{{L^2}}}} \end{split} \right. $$ (6) 式(6)是利用三阶求导的位置测量方法将所有量都转成了磁场分量，消除了矢量项，在很大程度上消除了地磁场和其他磁场的干扰。
2 基于卡尔曼滤波的姿态测量算法
利用磁信标和测量阵列得到底部钻具位置后，再通过陀螺仪和加速度计构成惯性导航系统进行姿态测量。由于导航系统中通过积分递推计算，陀螺仪的误差会随着时间累积；加速度计处于振动等动态环境下，测量结果也会有误差，但不会随时间累积；同时，磁场的变化会引起磁强计测量值的误差。因此，为了修正导航系统中的误差，采用卡尔曼滤波进行多源数据融合，利用外部量测信息来修正导航系统的误差[14]。首先将本文导航系统中涉及到的坐标系定义如图2所示。
1)导航坐标系
本文选取“北–东–地”地理坐标系作为导航参考坐标系，用
${o_n}{x_n}{y_n}{{\textit{z}}_n}$ 表示，${o_n}{x_n} $ 轴指向地理北向，${o_n}{y_n} $ 轴指向地理东向，${o_n}{\textit{z}_n}$ 轴垂直于当地椭球面指向地心；是惯导系统在求解导航参数时所采用的参考坐标系。2)载体坐标系
载体坐标系用
${o_b}{x_b}{y_b}{{\textit{z}}_b}$ 表示，选择MEMS传感器的中点为载体坐标系原点，$ {o_b}{x_b} $ 轴指向底部钻具的正左方，$ {o_b}{y_b} $ 指向底部钻具的正前方，${o_b}{{\textit{z}}_b}$ 轴指向底部钻具的正下方，构成标准右手系，记为$ b $ 系。根据水平定向钻进过程中的运动特征和工作环境[15]，设计了卡尔曼滤波模型的状态方程和量测方程。
2.1 状态方程的建立
基于惯性导航系统误差方程建立了15维的Kalman滤波模型，状态方程如下：
$$ {\boldsymbol{\dot X}}{\text{ = }}F{\boldsymbol{X}} + GW $$ (7) 式中：X为状态量；F为状态转移矩阵；W为系统噪声；
${\boldsymbol{G}}$ 为15维单位矩阵，是系统噪声转移矩阵，其中状态量选取北东地导航坐标系中的姿态误差、速度误差、位置误差、陀螺仪常值零偏和加速度计常值零偏：$$ {\boldsymbol{X}} = {\left[ {\delta \theta ,\delta \alpha ,\delta \varphi ,\delta {v_{\rm{N}}},\delta {v_{\rm{E}}},\delta {v_{\rm{D}}},\delta L,\delta \lambda ,\delta h,\varepsilon _x^b,\varepsilon _y^b,\varepsilon _{\textit{z}}^b,\nabla _x^b,\nabla _y^b,\nabla _{\textit{z}}^b} \right]^{\rm{T}}} $$ (8) 式中：
$ \delta \theta ,\delta \alpha ,\delta \varphi $ 分别为倾角、工具面角和方位角误差；$\delta {v_{\rm{N}}},\delta {v_{\rm{E}}},\delta {v_{\rm{D}}}$ 分别为北、东、地方向上的速度误差；$ \delta L,\delta \lambda ,\delta h $ 分别为纬度、经度、深度的位置误差；$\varepsilon _x^b,\varepsilon _y^b,\varepsilon _{\textit{z}}^b,\nabla _x^b,\nabla _y^b,\nabla _{\textit{z}}^b$ 分别为陀螺仪和加速度计沿载体坐标系三轴方向上的零偏。2.2 量测方程的建立
对于水平定向钻进的姿态测量包括钻进过程中的动态测量和停钻的静态测量，根据水平定向钻进的施工特点，本文在动态情况下利用
$ {x_b} $ 轴速度为零的准零速信息，在静态情况下利用三轴速度为零的零速信息以及磁信标测量阵列测得位置信息、磁强计和加速度计得到的姿态信息，对惯性测量解算的结果进行修正。建立量测方程如下：$$ {{\boldsymbol{Z}}_k}{\text{ = }}{{\boldsymbol{H}}_k}{\boldsymbol{X}} + {\boldsymbol{\upsilon }} $$ (9) 式中：
$ {{\boldsymbol{Z}}_k} $ 为量测量；$ {{\boldsymbol{H}}_k} $ 为量测转移矩阵；$ k $ 为当前状态($ k = 1 $ 为顶进状态，$ k = 2 $ 为静止状态)；$ {\boldsymbol{\upsilon }} $ 为量测噪声。1) 准零速修正
准零速是指当载体坐标系相对于导航坐标系某一方向上的运动速度始终为零，根据动态零速修正非完整约束条件，理想情况下，水平定向钻进坐标系
$ {x_b} $ 轴速度为0，但实际上由于惯性导航存在姿态漂移和振动干扰[16-17]，其输出载体坐标系下$ {x_b} $ 轴速度${\boldsymbol{v}}_{{\rm{INS}}x}^b$ 不为0。北东地坐标系n下的速度$ {{\boldsymbol{v}}^n} $ 到载体坐标系下的速度$ {{\boldsymbol{v}}^b} $ 转换表达式为：$$ {{\boldsymbol{v}}^b} = {\boldsymbol{C}}_n^b{{\boldsymbol{v}}^n} $$ (10) 式中：
$ {\boldsymbol{C}}_n^b $ 为准确的姿态转移矩阵，$ {\boldsymbol{C}}_n^b = {\left( {{\boldsymbol{C}}_b^n} \right)^{\rm{T}}} $ ，因此：$$ {\boldsymbol{v}}_{{\rm{INS}}}^b = ({\boldsymbol{C}}{_n^{b'}}){{\boldsymbol{v}}^n} $$ (11) 式中：
${\boldsymbol{C}}{_b^{n'}}$ 为惯性导航输出姿态矩阵，当3个平台失准角较小时，有如下转换关系[1]：$$ {\boldsymbol{C}}_n^b = [{\boldsymbol{E}} + \left( {{\text{φ}} \times } \right)]{\boldsymbol{C}}{_n^{b'}} $$ (12) 其中，
$ {\boldsymbol{E}} $ 为三阶单位矩阵，$$ \left( {{\text{φ}} \times } \right) = \left[ {\begin{array}{*{20}{c}} 0&{ - \delta \varphi }&{\delta \gamma } \\ {\delta \varphi }&0&{ - \delta \theta } \\ { - \delta \gamma }&{\delta \theta }&0 \end{array}} \right] $$ 因此，
$ {x_b} $ 轴准零速修正的速度误差可表示为：$$ \delta {{\boldsymbol{v}}^b_x} = {\boldsymbol{v}}_{{\rm{INS}}}^b - {{\boldsymbol{v}}^b} $$ (13) 2) 零速修正
水平定向钻进施工期间，会有更换钻杆的静止时间，这时水平定向钻进底部钻具的3个方向速度为0，由于导航解算中高度通道不稳定，不考虑惯导解算出的高度方向速度
${V_{\rm{D}}}$ [18]。选取北东地导航坐标系下东、北方向速度误差$\delta {V}_{{\rm{E}}}^{n}、\delta {V}_{{\rm{N}}}^{n}$ 作为量测量：$$ \left\{ \begin{split} & {\delta V_{\rm{E}}^n = V_{\rm{E}}^n - 0} \\ & {\delta V_{\rm{N}}^n = V_{\rm{N}}^n - 0} \end{split} \right. $$ (14) 3) 磁信标测量阵列位置测量信息
根据底部钻具上的测量阵列可以在静止时测量到当前底部钻具的位置
$\left[ {x,y,{\textit{z}}} \right]$ ，钻进的入口处测量阵列的测得位置为$\left[ {{x_0},{y_0},{{\textit{z}}_0}} \right]$ 。因此，底部钻具的位置信息为：$$ \left\{\begin{split} & {L_{\rm{M}}} = {L_0} + {{(x - {x_0})} \mathord{\left/ {\vphantom {{(x - {x_0})} {\left( {{R_1} + {\textit{z}}} \right)}}} \right. } {\left( {{R_1} + {\textit{z}}} \right)}} \\ & {\lambda _{\rm{M}}} = {\lambda _0} + {{(y - {y_0})} \mathord{\left/ {\vphantom {{(y - {y_0})} {\left( {\left( {{R_2} + {\textit{z}}} \right) \cdot \cos {L_D}} \right)}}} \right. } {\left( {\left( {{R_2} + {\textit{z}}} \right) \cdot \cos {L_D}} \right)}} \\ & {h_{\rm{M}}} = {\textit{z}} \end{split} \right. $$ (15) 式中：
$ {L_0} $ 、$ {\lambda _0} $ 分别为起始时刻测量阵列安装的纬度、经度信息；$ {R_1} $ 和$ {R_2} $ 分别为当地地球的子午圈和卯酉圈半径；$h_{\rm{M}}$ 为当前高度。4) 基于加速度计和磁强计的姿态测量信息
在静止状态，还可以利用安装在底部钻具的测量阵列中的三轴磁强计和三轴加速度计可以进行姿态测量(图3)。当底部钻具处于静止状态时，MEMS传感器只受到重力加速度的作用，加速度计的测量值是重力加速度的3个分量值为
${\left[ {\begin{array}{*{20}{c}} {{g_x}}&{{g_y}}&{{g_\textit{z}}} \end{array}} \right]^{\rm{T}}}$ [19-20]。北东地导航坐标系中，只有垂直方向的重力场，可以表示为${\left[ {\begin{array}{*{20}{c}} 0&0&g \end{array}} \right]^{\rm{T}}}$ 。其中，地磁场可以表示为${\left[ {\begin{array}{*{20}{c}} 0&{{{\boldsymbol{B}}_{\boldsymbol{e}}}\cos \beta }&{{{\boldsymbol{B}}_{\boldsymbol{e}}}\sin \beta } \end{array}} \right]^{\rm{T}}}$ ，当钻具静止时，可以根据三轴加速度计计算得到倾角和工具面角，通过三轴磁强计可以计算得到方位角，再从导航坐标系转换到钻具所处的载体坐标系${o_b}{x_b}{y_b}{\textit{z}_b}$ 后可得钻具的倾角${\theta _{\rm{M}}}$ 、工具面角${\alpha _{\rm{M}}}$ 和方位角${\varphi _{\rm{M}}}$ 的表达式：![]()
图 3 姿态测量传感器安装方法及坐标系建立Figure 3. Installation method of attitude measurement sensor and establishment of coordinate system$$ \begin{aligned} & {\alpha _M} = \arctan \left( - \frac{{{g_x}}}{{\sqrt {{g_x}^2 + {g_y}^2} }}\right) \\ & {\theta _M} = \arctan \frac{{{g_y}}}{{{g_{\textit{z}}}}} \\ & {\varphi _M} = \arctan \left(\frac{{g({B_{\textit{z}}}{g_y} - {B_y}{g_{\textit{z}}})}}{{{B_x}(g_y^2 + g_{\textit{z}}^2) - {g_x}({B_y}{g_y} - {B_{\textit{z}}}{g_{\textit{z}}})}}\right)\end{aligned} $$ (16) 式中：
$g = \sqrt { {{g^2_x}}+{{g^2_y}}+{g^2_\textit{z}} }$ ，为测量的重力加速度值。将本节得到4项量测信息：准零速信息、零速信息、磁信标测量阵列位置测量信息以及加速度计和磁强计的姿态测量信息，输入到卡尔曼滤波器中进行数据融合，算法流程如图4所示。图中，MEMS-IMU指微机电惯性测量元件，DZUPT为准零速修正，ZUPT为零速修正，
$[x,y,{\textit{z}}]$ 为磁信标测量阵列测得的位置信息。![]()
图 4 基于卡尔曼滤波的姿态测量算法流程Figure 4. Process flow of Kalman filter-based attitude measurement algorithm动态时采用的量测量为
$ {{\boldsymbol{Z}}_1} $ ：$$ {{\boldsymbol{Z}}_1}{\text{ = }}\left[ {\delta {{\boldsymbol{v}}^b_x}} \right] $$ (17) 静态采用的量测量为
$ {{\boldsymbol{Z}}_2} $ ，式中下标IMU指MEMS-IMU惯性解算得到的姿态、位置和速度。$$ {{\boldsymbol{Z}}_2}{\text{ = }}\left[ {\begin{array}{*{20}{c}} {{\alpha _{{\rm{IMU}} - }}{\alpha _M}} \\ {{\theta _{{\rm{IMU}}}} - {\theta _M}} \\ {{\varphi _{{\rm{IMU}}}} - {\varphi _M}} \\ {{L_{{\rm{IMU}}}} - {L_M}} \\ {{\lambda _{{\rm{IMU}}}} - {\lambda _M}} \\ {{h_{{\rm{IMU}}}} - {h_M}} \\ {\delta V_E^n} \\ {\delta V_N^n} \end{array}} \right] $$ (18) 3 模拟实验系统设计与结果分析
3.1 磁信标仿真与设计
利用COMSOL物理场仿真软件模拟磁信标磁场分布，仿真参数为单层密绕，磁信标磁芯长度为800 mm，直径为70 mm，密绕漆包线长度为600 mm。
漆包线线径相对磁信标横截面直径小得多，若进行多层密绕，可认为磁信标产生磁场强度是叠加的，若采用
$ n $ 层密绕，磁信标产生磁场强度是单层密绕磁信标产生磁场强度的$ n $ 倍。本文仿真中$ n $ 值取1，相对磁导率为10 000，根据麦克斯韦电磁场理论，仿真结果如图5所示。![]()
图 5 加入磁芯的磁信标磁场分布Figure 5. Magnetic field distribution of a magnetic beacon with a magnetic core仿真结果显示，磁信标中心产生的磁场强度到达到
$2 \times {10^8}\;{\rm{n}}{\rm{T}}$ ，仿真磁信标的空间磁场大小见表1。磁感应强度在距离磁信标原点10 m处轴向磁场强度为54 nT，可以被磁强计测量到。同时考虑到一般钻进距离、安全性和便携性，本文设计的磁信标的主要参数为4层密绕漆包线，漆包线直径为2.2 mm，密绕长度为600 mm，单层为250匝，磁芯长度为800 mm、磁芯直径为70 mm，相对磁导率为10000。制作的磁信标经过现场实际测量，磁场强度满足要求。表 1 磁信标空间位置磁场强度Table 1. Magnetic field strength at the spatial position of the magnetic beacon距离中心点位置/m 径向场强/nT 轴向场强/nT 10 125.77 54 8 423.97 118.81 6 1204.5 459.3 4 3340.4 1000.1 2 93917 24923 3.2 场景模拟实验
为验证本文提出算法和设计磁信标的精度和可行性，进行了位置和姿态的场景模拟实验。
3.2.1 位置测量实验方案与结果分析
位置测量实验场地选择在变化磁场相对较少的空旷操场上，如图6所示。磁信标架设在距离地面2.2 m处，以磁信标中心为原点建立参考坐标系
$oxy{\textit{z}}$ ，轴向为$ x $ 轴，径向为$ y $ 轴。测量方案如下：沿
$ y $ 轴方向共设置101个测点，等距排列，间隔距离为20 cm，测量范围为$ y = $ ±10 m，测点均位于距离$ x $ 轴3 m处。采用手持式激光测距仪测量每个测点位置作为真实值。根据场地磁信标布置位置，采用米尺划定$ y $ 轴方向测量范围，并标记101个测点的位置，方便后续测量，按照测点依次排放测量阵列，按照三阶磁梯度张量定位算法，正反接入电源进行测量，得到位置的测量值，再利用手持激光仪测量真实值。实验中为保证位置真实值测量准确，采用的手持激光仪的参数见表2。表 2 手持式激光测距仪技术参数Table 2. Technical parameters of the handheld laser range finder参数 数值 测量距离/m 0~100.0 30 m内测量精度/mm ±1.5 最小显示单位/mm 1 激光等级 不低于Ⅱ 激光点直径精度/(mm·m−1) 0.6 根据设计的磁信标在10 m处产生的磁场强度，可以被测量阵列所用磁强计敏感到，本文所用磁强计型号为RM3100，参数见表3。
表 3 测量阵列中磁强计技术参数Table 3. Technical parameters of the magnetometer in the measurement array参数 参数值 动态范围/
Gauss±13 比例系数
稳定性/% F·S0.5 非线性/
% F·S1.0 偏移量/
mGauss1.0 根据位置测量结果，测量值和真实值的差值为测量误差，相对误差为测量误差相对于各测量轴方向总长度所得长度乘以100%。以
$ x $ 轴为例，本次实验的$ x $ 轴测量距离为3 m，所测得的最大误差为$0.278\;{\rm{m}}$ ，最大相对误差为$0.278 \div 3 \times 100\text{%} = 9.27\text{%}$ ，101个测点的平均误差为$0.135\;{\rm{m}}$ ，平均相对误差为$0.135 \div 3 \times 100\text{%} = 4.52\text{%}$ ，其余各轴同理计算，x、y、z轴误差见表4。表 4 各轴误差Table 4. Error of various axes位置测量 最大误差/
m相对误差/
%平均误差/
m平均相对误差/
%$ x $轴 0.278 9.27 0.135 4.52 y轴 0.702 7.02 0.306 3.06 $ {\textit{z}} $轴 0.261 11.86 0.135 6.14 经位置实验验证，提出的基于三阶磁梯度张量的位置测量方法的三轴平均相对误差为
$(4.52 + 3.06 + 6.14) \div 3 \times 100\text{%} = 4.57\text{%}$ ，满足工程的需要。位置测量实验表明距离磁信标位置越近，定位误差越小，可以通过提升测量阵列磁强计和磁信标测量精度和磁场强度，提升测量距离。3.2.2 姿态测量实验方案与结果分析
姿态测量实验采用的MEMS传感器型号为XNA100C，加速度计和陀螺仪技术相关参数见表5。
表 5 测量阵列中加速度计和陀螺仪技术参数Table 5. Technical parameters of accelerometer and gyroscope in the measurement array陀螺仪参数 数值 加速度计参数 数值 带宽/Hz 100 带宽/Hz 100 零偏稳定性/[(°)·hr−1] 12 零偏稳定性/10−6g 30 偏移量/[(°)·s−1] 0.003 偏移量/10−3g 0.1 随机游走/[(°)·hr−1/2] 0.2 随机游走/(10−6g·hr−1/2) 340 动态范围/[(°)·s−1] ±300 动态范围/g ±8 姿态测量实验将MEMS传感器固定在Wi-Fi履带机器人上模拟钻进过程中的姿态变化(图10)。实验设定姿态测量轨道长度为12.6 m，考虑到水平定向钻进主要进行竖直方向的造斜，以及实验场地限制，实验主要模拟了倾角和方位角变化。
设置模拟钻进速度为0.07 m/s，总行进时长180 s，两端坡度为16°。具体钻进实验参数见表6。
表 6 实验中轨迹参数设置Table 6. Trajectory parameter settings in the experiments时间/s 方位角/(°) 工具面角/(°) 倾角/(°) 0~40 0 0 −16 40~70 0 0 0 70~100 15 0 0 100~120 −5 0 0 120~140 −20 0 0 140~165 0 0 0 165~180 0 0 16 传感器开机稳定后按照设定轨迹模拟钻进。对采集数据进行多源数据融合，得到倾角、工具面角和方位角3个姿态角的最优估计值如图11—图13。图中，红色曲线是陀螺仪解算出的姿态角，绿色曲线是加速度计和磁强计算得到的姿态角，蓝色曲线是经过多源数据融合之后得到姿态角。
图11为倾角的测量结果。从图中可以看出，陀螺仪解算的倾角随时间发生了明显的漂移；磁强计和加速度计解算的倾角受环境影响，测量噪声大；采用本文提出的多源数据融合算法解算的倾角可以看出，倾角由16°变为0°再变为−16°，与表6中设置的倾角参数一致，且测量误差小于0.3°说明了本方法的有效性。
图12为工具面角的测量结果。从图中可以看出，陀螺仪解算的工具面角随时间发生了明显的漂移，并且存在10°左右的波动；磁强计和加速度计解算的工具面角没有发生明显漂移，但与真实值相差较多，采用多源数据融合算法解算的工具面角，与真实值0°基本一致，测量误差小于0.7°，但测量场地不是完全水平的，会造成工具面角在误差允许范围±0.5°的小范围内波动。
图13为方位角的测量结果。从图中可以看出，陀螺仪解算的方位角有实验设定值的对应变化趋势，但随时间发生了明显的漂移；磁强计和加速度计解算的方位角有很大的波动，噪声完全淹没了有效信息；采用多源数据融合算法解算的方位角，在0~70 s、140~180 s在0°附近波动，75~100 s在15°附近波动，105~120 s在−5°附近波动，125~140 s在−20°附近波动，与实验设定基本一致，测量误差小于1°，验证了本文方法的正确性和有效性，但前70 s的部分与实验设定相比偏移较大，主要原因是在操控小车沿设定轨迹行驶时进行了小的方向调整。三个姿态角中方位角相对于其他两个角的测量误差较大，这主要由于在导航解算系统中，融合数据对方位角修正效果相比其他两个角较弱，但测量精度满足水平定向钻进工程实际应用的需求。
4 结 论
a. 设计了一种底部钻具的测量阵列，提出了基于三阶张量的磁梯度张量目标定位方法，来消除水平定向钻进位置算法中地磁场的干扰，并通过仿真确定了磁信标的参数，制作了满足磁场强度的磁信标。
b. 通过提出采用多源数据融合的姿态测量方法，将磁信标测量阵列位置信息、加速度计和磁强计的姿态信息以及准零速、零速信息进行卡尔曼滤波融合，来修正导航系统中的姿态误差，得到最优估计角度。
c. 通过位置测量实验，在20 m范围内取101个测点进行基于三阶磁梯度张量算法的位置计算。结果表明，位置测量三轴的平均相对误差为4.57%，在距离磁信标越远处测量误差越大，为保证远距离的测量精度，可以提升磁信标磁场强度或使用高精度的磁强计，还可以通过增设磁信标延长测量距离。
d. 姿态测量实验利用履带机器人作为载体，进行连续钻进姿态角解算实验，倾角、工具面角和方位角的误差值在1°以内，证明了姿态算法的准确性和卡尔曼滤波效果的有效性。测量精度满足水平定向钻进的要求。本研究暂未进行现场试验，今后将进一步将测量阵列结构化，设计适用于水平定向钻进的测量短节，为水平定向钻进提供一种高精度导向测量技术。
-
![]()
图 1 基于三阶磁梯度张量的测量阵列
Fig. 1 Measurement array based on the third-order magnetic gradient tensor
![]()
图 3 姿态测量传感器安装方法及坐标系建立
Fig. 3 Installation method of attitude measurement sensor and establishment of coordinate system
![]()
图 4 基于卡尔曼滤波的姿态测量算法流程
Fig. 4 Process flow of Kalman filter-based attitude measurement algorithm
![]()
图 5 加入磁芯的磁信标磁场分布
Fig. 5 Magnetic field distribution of a magnetic beacon with a magnetic core
表 1 磁信标空间位置磁场强度
Table 1 Magnetic field strength at the spatial position of the magnetic beacon
距离中心点位置/m 径向场强/nT 轴向场强/nT 10 125.77 54 8 423.97 118.81 6 1204.5 459.3 4 3340.4 1000.1 2 93917 24923 
下载: 导出CSV
表 2 手持式激光测距仪技术参数
Table 2 Technical parameters of the handheld laser range finder
参数 数值 测量距离/m 0~100.0 30 m内测量精度/mm ±1.5 最小显示单位/mm 1 激光等级 不低于Ⅱ 激光点直径精度/(mm·m−1) 0.6
下载: 导出CSV
表 3 测量阵列中磁强计技术参数
Table 3 Technical parameters of the magnetometer in the measurement array
参数 参数值 动态范围/
Gauss±13 比例系数
稳定性/% F·S0.5 非线性/
% F·S1.0 偏移量/
mGauss1.0
下载: 导出CSV
表 4 各轴误差
Table 4 Error of various axes
位置测量 最大误差/
m相对误差/
%平均误差/
m平均相对误差/
%$ x $轴 0.278 9.27 0.135 4.52 y轴 0.702 7.02 0.306 3.06 $ {\textit{z}} $轴 0.261 11.86 0.135 6.14
下载: 导出CSV
表 5 测量阵列中加速度计和陀螺仪技术参数
Table 5 Technical parameters of accelerometer and gyroscope in the measurement array
陀螺仪参数 数值 加速度计参数 数值 带宽/Hz 100 带宽/Hz 100 零偏稳定性/[(°)·hr−1] 12 零偏稳定性/10−6g 30 偏移量/[(°)·s−1] 0.003 偏移量/10−3g 0.1 随机游走/[(°)·hr−1/2] 0.2 随机游走/(10−6g·hr−1/2) 340 动态范围/[(°)·s−1] ±300 动态范围/g ±8
下载: 导出CSV
表 6 实验中轨迹参数设置
Table 6 Trajectory parameter settings in the experiments
时间/s 方位角/(°) 工具面角/(°) 倾角/(°) 0~40 0 0 −16 40~70 0 0 0 70~100 15 0 0 100~120 −5 0 0 120~140 −20 0 0 140~165 0 0 0 165~180 0 0 16
下载: 导出CSV
-
[1] 姚宁平,王力,张金宝,等. 煤矿井下连续管钻进管柱分析及射流钻进实验[J]. 煤田地质与勘探,2023,51(1):298−308. YAO Ningping,WANG Li,ZHANG Jinbao,et al. Analysis and jet drilling test of coiled tubing drilling string for underground coal mine[J]. Coal Geology & Exploration,2023,51(1):298−308.
[2] PATINO–RAMIREZ F,LAYHEE C,ARSON C. Horizontal directional drilling (HDD) alignment optimization using ant colony optimization[J]. Tunnelling and Underground Space Technology,2020,103:103450. DOI: 10.1016/j.tust.2020.103450
[3] 王小波. 无线电磁波随钻测量系统姿态精度的影响因素分析[J]. 煤田地质与勘探,2021,49(6):258−264. WANG Xiaobo. Factors affecting the attitude accuracy of wireless electromagnetic wave MWD system[J]. Coal Geology & Exploration,2021,49(6):258−264.
[4] 刘金祯,金键,陈雪华,等. 地面信标控向系统在海河穿越中的应用[J]. 石油工程建设,2003,29(6):58−60. LIU Jinzhen,JIN Jian,CHEN Xuehua,et al. Application of ground beacon directional control system in pipeline crossing through Haihe River[J]. Petroleum Engineering Construction,2003,29(6):58−60.
[5] 孙平贺,刘伟胜,杨涵涵,等. 中国非开挖水平定向钻进装备与技术研究应用进展[J]. 工程科学学报,2022,44(1):122−130. SUN Pinghe,LIU Weisheng,YANG Hanhan,et al. Progress in research and applications of trenchless horizontal directional drilling equipment and technology in China[J]. Chinese Journal of Engineering,2022,44(1):122−130.
[6] LIU Tao,WANG Boxiong,CUI Yuanyuan,et al. Direction and position measurement in HDD using two magnetic fields[J]. Sensors and Actuators A:Physical,2012,185:168−172. DOI: 10.1016/j.sna.2012.08.005
[7] ZHANG Yuexin. A fusion methodology to bridge GPS outages for INS/GPS integrated navigation system[J]. IEEE Access,2019,7:61296−61306. DOI: 10.1109/ACCESS.2019.2911025
[8] 常帅,付晓梅,张翠翠,等. 基于磁信标的水下SLAM方法[J]. 水下无人系统学报,2019,27(3):277−283. CHANG Shuai,FU Xiaomei,ZHANG Cuicui,et al. An underwater SLAM approach using magnetic beacons[J]. Journal of Unmanned Undersea Systems,2019,27(3):277−283.
[9] ZHENG Yuanxun,LI Qinghua,WANG Changhong,et al. Magnetic–based positioning system for moving target with feature vector[J]. IEEE Access,2020,8:105472−105483. DOI: 10.1109/ACCESS.2020.3000305
[10] 石智军,许超,李泉新. 煤矿井下近水平随钻测量定向钻孔轨迹设计与计算方法[J]. 煤田地质与勘探,2015,43(4):112−116. SHI Zhijun,XU Chao,LI Quanxin. Trajectory design and calculation of nearly horizontal MWD directional borehole in underground coal mine[J]. Coal Geology & Exploration,2015,43(4):112−116.
[11] 吕俊伟,迟铖,于振涛,等. 磁梯度张量不变量的椭圆误差消除方法研究[J]. 物理学报,2015,64(19):190701. DOI: 10.7498/aps.64.190701 LYU Junwei,CHI Cheng,YU Zhentao,et al. Research on the asphericity error elimination of the invariant of magnetic gradient tensor[J]. Acta Physica Sinica,2015,64(19):190701. DOI: 10.7498/aps.64.190701
[12] LEE K M,LI Min. Magnetic tensor sensor for gradient–based localization of ferrous object in geomagnetic field[J]. IEEE Transactions on Magnetics,2016,52(8):4002610.
[13] 李光,随阳轶,刘丽敏,等. 基于差分的磁偶极子单点张量定位方法[J]. 探测与控制学报,2012,34(5):50−54. LI Guang,SUI Yangyi,LIU Limin,et al. Magnetic dipole single–point tensor positioning based on the difference method[J]. Journal of Detection & Control,2012,34(5):50−54.
[14] LI You,ZHUANG Yuan,ZHANG Peng,et al. An improved inertial/wifi/magnetic fusion structure for indoor navigation[J]. Information Fusion,2017,34:101−119. DOI: 10.1016/j.inffus.2016.06.004
[15] 张浩. 非开挖管道定向穿越技术[J]. 中国石油和化工标准与质量,2012(8):115. ZHANG Hao. Trenchless pipeline directional crossing technology[J]. China Petroleum and Chemical Standard and Quality,2012(8):115.
[16] 赵玉,赵忠,范毅. 零速修正技术在车载惯性导航中的应用研究[J]. 压电与声光,2012,34(6):843−847. ZHAO Yu,ZHAO Zhong,FAN Yi. Study on application of zero velocity update technology to inertial navigation system[J]. Piezoelectrics & Acoustooptics,2012,34(6):843−847.
[17] LIU Xixiang,WANG Songbing,ZHANG Tongwei,et al. A zero–velocity detection method with transformation on generalized likelihood ratio statistical curve[J]. Measurement,2018,127:463−471. DOI: 10.1016/j.measurement.2018.05.113
[18] ZHANG Peng,LI You,ZHUANG Yuan,et al. Multi–level information fusion with motion constraints:Key to achieve high–precision gait analysis using low–cost inertial sensors[J]. Information Fusion,2023,89:603−618. DOI: 10.1016/j.inffus.2022.09.009
[19] 黄欣,熊智,许建新,等. 基于零速/航向自观测/地磁匹配的行人导航算法研究[J]. 兵工学报,2017,38(10):2031−2040. DOI: 10.3969/j.issn.1000-1093.2017.10.020 HUANG Xin,XIONG Zhi,XU Jianxin,et al. Research on pedestrian navigation algorithm based on zero velocity update/heading error self–observation/geomagnetic matching[J]. Acta Armamentarii,2017,38(10):2031−2040. DOI: 10.3969/j.issn.1000-1093.2017.10.020
[20] 燕斌,程建远,蔡远利,等. 基于MEMS陀螺的矿用新型钻机开孔定向仪研制[J]. 煤田地质与勘探,2018,46(6):187−192. YAN Bin,CHENG Jianyuan,CAI Yuanli,et al. Development of a new type of drilling orientation device based on MEMS gyroscope[J]. Coal Geology & Exploration,2018,46(6):187−192.
-
期刊类型引用(2)
1. 李浩，姚宁平，陆承达，张金宝，魏宏超，吴加俊，李金宇. 煤矿井下定向钻孔轨迹模型预测控制方法. 煤田地质与勘探. 2025(02): 205-212 . 
本站查看
2. 刘东，郑媛媛. 水平定向钻钻进施工难题及合格率提升对策研究. 工程技术研究. 2024(22): 152-154 . 百度学术
其他类型引用(1)
计量
- 文章访问数: 209
- HTML全文浏览量: 12
- PDF下载量: 31
- 被引次数: 3
