CN102787838B - 一种基于Kalman滤波的改进SAGD算法 - Google Patents

Abstract

本发明涉及磁场定位技术领域。基于Kalman滤波的改进SAGD算法，包括：①选取若干目标点作为相对导航定位的参考点，目标点间隔取为预定双井间距Δl的4倍；②磁传感器放置于进尺为z₂的目标点，钻头进尺z₁由z₂-2Δl变化至z₂+2Δl的位置，磁传感器记录整个作用距离的磁场数据，并输入磁场模型；③由已知的双井进尺z₁和z₂计算得到多个r值；④按照磁场模型解算出多个其双井相对的垂直偏离角度θ；⑤将多个测点的r和θ值输入Kalman滤波器，按一定的状态转移规则实施一步预测；⑥对r和θ进行滤波更新，得到其滤波估计；重复⑤-⑥，直至所有测点输入完毕，此时的滤波估计即为最终的定位结果。本方法可保证钻井过程连续性及抗噪声和抗干扰能力。

Description

一种基于Kalman滤波的改进SAGD算法

技术领域

本发明涉及磁场定位技术领域，尤其涉及基于SAGD技术在低信噪比或突发干扰情况下测距导向的不足提出一种基于Kalman滤波的改进算法。

背景技术

据统计，世界稠油、超稠油和天然沥青的储量约为1000×10⁸t，占有油气总储量较大的比例，对于稠油资源相对丰富的中国、美国和加拿大等国家来说，如何经济有效的开采稠油成为石油界的一个重大课题。传统开采稠油是采用热力采油，包括蒸汽吐、蒸汽驱等技术，其对于普通稠油的开采十分成功。然而，传统方法对于地下原油粘度超过10⁴mPa·s的特稠油开采的经济效益很差，而对于粘度超过5×10⁴mPa·s，甚至高达10⁵mPa·s的超稠油、沥青等，传统方法根本无法开采。20世纪70年代末，加拿大石油工业专家R.M.Butler提出了蒸汽辅助重力泄油(Steam AssistedGravity Drainage,SAGD)技术，采用蒸汽驱开采方式，向注汽井连续注入高温、高干度蒸汽，蒸汽向上超覆在地层中形成蒸汽腔，蒸汽腔向上及侧面扩展，与油层中的原油发生热交换，加热后的原油和蒸汽冷凝水靠重力作用泄到下面的水平生产井中产出。生产水平井一般是在接近油柱底部油水界面上，蒸汽通过该水平井上方与前者相平行的第2口水平井或一系列垂直井持续注入，从而在生产井上方形成蒸汽室。SAGD与水平井技术相结合被认为是20世纪末所建立的最著名的油藏工程理论。

SAGD的布井方式主要有双水平井方式、直井水平井联合方式及单井SAGD。实验表明，双水平井情况下的采油率最高，效果最佳，在相同原油和蒸汽压力的情况下，水平井的采油速度大约是垂直井的2-3倍。在水平双井的钻井过程中，通常先打通一口水平井，然后在距离前者的正上方或正下方3-5m的位置打通第二口平行的水平井。SAGD布设水平双井的关键在于保证两口水平井处于同一垂直平面内。因此，在第二口水平井的钻进过程中应对钻进位置进行导航定位，在钻头偏离预定平面时及时给于修正。

20世纪90年代，Arthur F.Kuckes等人在专利US005589775A中提出一种针对SAGD技术的导航定位系统，并设计相应算法(以下简称SAGD算法)。如图1所示，该系统在钻头后部安装一永磁体并随钻头转动，在钻进过程中产生一个不断移动的时变磁场。而在已有水平井中，放置一个磁传感器，并向钻进方向移动位置，保持与钻头一定的进尺间距。通过对磁场数据的分析，可以得到两口水平井的相对距离及垂直偏离角度等信息，实现对钻进位置的导航定位。然而，该算法在钻头钻进过程中仅通过对特定位置磁场的一次测量进行解算，在低信噪比或突发干扰情况下，其定位结果的准确性将变差。

发明内容

(一)要解决的技术问题

本发明要解决的技术问题是提供一种仅在一次连续钻进过程中，在不同位置实现对同一点的多次测量，保证钻井过程连续性及抗噪声和干扰的基于Kalman滤波的改进SAGD算法。

(二)技术方案

为了解决上述问题，本发明提供了一种基于Kalman滤波的改进SAGD算法，包括以下步骤：

①、在已有井中选取若干目标点作为相对导航定位的参考点，目标点间隔取为预定双井间距Δl的4倍，磁传感器将随钻进过程放置于不同目标点，划分各目标点作用距离；

②、在目标点作用距离内，磁传感器记录钻进过程中的磁场数据，计算总磁场的归一化幅值；

③、选取若干离散的归一化幅值作为测点，各测点分别由SAGD算法计算得到多个相对距离r；

④、由上述测点中提取x、y方向的磁场分量，按照磁场模型解算出多个其双井相对的垂直偏离角度θ；

⑤、将多个测点的r和θ值依次输入Kalman滤波器，对于第m次滤波，r_m和θ_m分别按照状态转移规则实施状态转移，给出一步预测值；

⑥、对r和θ进行滤波更新，得到其滤波估计；

重复⑤-⑥，直至所有测点输入完毕，此时的滤波估计即为最终的定位结果。

进一步的技术方案是，在步骤①中，进尺为z₂的目标点，其作用距离划分为钻头进尺z₁由z₂-2Δl变化至z₂+2Δl的位置。

进一步的技术方案是，步骤②中计算归一化幅值的具体方法是，

磁场模型接收磁传感器的三轴磁场数据及控制设备给出的钻头进尺，记为(B_x(t),B_y(t),B_z(t),z(t))。

将采样数据按磁场变化周期进行分段，如图2所示。其中，第i段数据表示如下

$\begin{array}{c}{B}_{x,i}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ B_{y,i}\left(t\right)=B_y\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ B_{z,i}\left(t\right)=B_z\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\end{array}---\left(1\right)$

对B_x,i(t)B_y,i(t)B_z,i(t)分别取希尔伯特变换H[B_x,i(t)]H[B_y,i(t)]H[B_z,i(t)]，并求取角度

$\begin{array}{c}{\mathrm{\&phi;}}_{x,i}\left(t\right)=\arctan \frac{H\left[B_{x,i}\left(t\right)\right]}{B_{x,i}\left(t\right)}\\ {\mathrm{\&phi;}}_{y,i}\left(t\right)=\arctan \frac{H\left[B_{y,i}\left(t\right)\right]}{B_{y,i}\left(t\right)}\\ {\mathrm{\&phi;}}_{z,i}\left(t\right)=\arctan \frac{H\left[B_{z,i}\left(t\right)\right]}{B_{z,i}\left(t\right)}\end{array}---\left(2\right)$

以表示B_x,i(t)和B_y,i(t)在该周期内的平均幅值，则有

$\begin{array}{c}\left|{\stackrel{\&OverBar;}{B}}_x\left(i\right)\right|=\frac{1}{\mathrm{\&pi;}}\sqrt{{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,i}\left(t\right)\sin {\mathrm{\&phi;}}_{x,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2+{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,i}\left(t\right)\cos {\mathrm{\&phi;}}_{x,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2}\\ \left|{\stackrel{\&OverBar;}{B}}_y\left(i\right)\right|=\frac{1}{\mathrm{\&pi;}}\sqrt{{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,i}\left(t\right)\sin {\mathrm{\&phi;}}_{y,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2+{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,i}\left(t\right)\cos {\mathrm{\&phi;}}_{y,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2}\end{array}---\left(13\right)$

定义

$B_{\mathrm{mag}}\left(i\right)=\sqrt{{\left|{\stackrel{\&OverBar;}{B}}_x\left(i\right)\right|}^2+{\left|{\stackrel{\&OverBar;}{B}}_y\left(i\right)\right|}^2}---\left(14\right)$

搜索B_mag的最大值B_max＝max(B_mag)

对B_mag其进行归一化处理：

$B_{\mathrm{mag}}^{\&prime;}=\frac{B_{\mathrm{mag}}}{B_{\max }}---\left(15\right)$

进一步的技术方案是，步骤③计算双井相对距离具体过程为：

SAGD算法给出的磁场模型有如下归一化幅值与Δz/r关系式

$B_{\mathrm{mag}}^{\&prime;}=\frac{\sqrt{2{(\mathrm{\&Delta;z}/r)}^4-2{(\mathrm{\&Delta;z}/r)}^2+5}}{\sqrt{5}{[1+{(\mathrm{\&Delta;z}/r)}^2]}^{\frac{5}{2}}}---\left(16\right)$

其中，Δz＝z₁-z₂，r为两口水平井实际垂直距离；

绘制上式的关系曲线图，选取若干离散归一化数值{η₁,η₂,...,η_M}，对于每一个η_m，在图3中找到满足公式(6)的两个解a_m和b_m；同时，在步骤②得到的B′_mag中搜索找到两个归一化磁场满足

B′_mag(i_m)＝B′_mag(j_m)＝η_m         (17)

并由脚标提取相应的进尺

$\begin{array}{c}{z}_{m,1}=z\left(t\right){|}_{t=\frac{2\mathrm{\&pi;}}{\mathrm{\&omega;}}{i}_m}\\ z_{m,2}=z\left(t\right){|}_{t=\frac{2\mathrm{\&pi;}}{\mathrm{\&omega;}}{j}_m}\end{array}---\left(18\right)$

进一步解算出相对距离

$r_m=\left|\frac{z_{m,1}-z_{m,2}}{a_m-b_m}\right|---\left(19\right)$

进一步的技术方案是，步骤④计算双井垂直偏离角度的具体过程为：

重新将作用时段距离内的数据按如下形式分段：

$\begin{array}{c}{B}_{x,m}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}{i}_m<\mathrm{\&omega;t}<2\mathrm{\&pi;}{i}_{m+1}},m=\mathrm{1,2},...,M-1\\ B_{x,m}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}{i}_m<\mathrm{\&omega;t}<2\mathrm{\&pi;}{j}_m},m=M\end{array}---\left(20\right)$

B_y,m,φ_x,m,φ_y,m,φ_z,m分别由与公式(2)相似形式进行分段给出；

由磁场模型得到如公式(11)所示的θ值计算公式，并将第m段数据输入磁场模型得到的相对角度记为θ_m

${\mathrm{\&theta;}}_m=a\tan 2[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,m}\left(t\right)\sin {\mathrm{\&phi;}}_{z,m}\left(t\right)d\left(\mathrm{\&omega;t}\right),-{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,m}\left(t\right)\sin {\mathrm{\&phi;}}_{z,m}\left(t\right)d\left(\mathrm{\&omega;t}\right)]---\left(21\right)$

进一步的技术方案是，步骤⑤对r和θ一步预测具体过程为：

定义

r_mθ_m：第m次滤波的测量值

由第m-1次滤波的一步预测值

第m次滤波估计值

p_r,mp_θ,m：第m次滤波方差

对于第m(2≤m≤M)次滤波前的一步预测，其状态转移规则为：

$\begin{array}{c}{\hat{r}}_{m|m-1}=\sqrt{A^2+B^2}\\ {\widehat{\mathrm{\&theta;}}}_{m|m-1}=a\tan 2(B,A)\end{array}---\left(22\right)$

其中，

上式中，

$\mathrm{\&Delta;}{z}_m=\left\{\begin{array}{c}{z}_{m+\mathrm{1,1}}-z_{m,1},m<M\\ z_{m,2}-z_{m,1},m=M\end{array}\right..$

当m＝1时， ${\hat{r}}_1=r_1{\widehat{\mathrm{\&theta;}}}_1={\mathrm{\&theta;}}_1$

进一步的技术方案是，步骤⑥滤波更新具体步骤为：

将一步预测值和测量值输入Kalman滤波方程

$\begin{array}{c}{\hat{r}}_m={\hat{r}}_{m|m-1}+K_{r,m}(r_m-{\hat{r}}_{m|m-1})\\ {\widehat{\mathrm{\&theta;}}}_m={\widehat{\mathrm{\&theta;}}}_{m|m-1}+K_{r,m}({\mathrm{\&theta;}}_m-{\widehat{\mathrm{\&theta;}}}_{m|m-1})\end{array}---\left(23\right)$

其中，

$K_{r,m}=\frac{p_{r,m}}{p_{r,m}+{\mathrm{\&sigma;}}^2},K_{\mathrm{\&theta;},m}=\frac{p_{\mathrm{\&theta;},m}}{p_{\mathrm{\&theta;},m}+{\mathrm{\&sigma;}}^2}$

上式中，σ²为噪声方差，

重复步骤⑤和⑥，共进行M次滤波，输出和作为最终定位结果。

(三)有益效果

本发明采用数字信号处理中的滤波方法对SAGD算法进行改进，改进SAGD算法具有以下特点和优势：

1、仅在一次连续钻进过程中，在不同位置实现对同一点的多次测量，保证了钻井过程的连续性。

2、距离解算过程采用归一化磁场，消除了地下不同区域磁场衰减率不同对定位的影响。

3、通过对多次测量结果进行Kalman滤波，降低了解算位置的方差，有效提高导航定位的精度。

4、有效弥补SAGD算法在低信噪比或突发干扰情况下导航定位能力不足的缺陷，具有较强的抵抗噪声和干扰的能力。

附图说明

图1为本发明基于Kalman滤波的改进SAGD算法结构原理图；

图2为标记点设置示意图；

图3三轴磁场数据按变化周期分段示意图；

图4为相对距离解算示意图；

图5为标记点1采集的磁场数据图；

图6和图7为Kalman滤波过程图。

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本发明提出的改进SAGD算法，具体实施的原理结构如图1所示。包括钻头、旋转磁体、磁传感器、标记点、控制设备、磁场模型和Kalman滤波器。钻头在目标井中不断向前钻进，旋转磁体安装于钻头上，产生周期时变的磁场。控制设备控制钻头进尺并根据定位结果调整钻头姿态，保持目标井与已有井平行且相对距离固定，同时将钻头进尺和姿态信息送入Kalman滤波器。标记点设置于已有井中，作为钻头钻进过程中导航定位的参考点，每当钻头经过标记点时进行一次定位。磁传感器放置于标记点上，采集三轴磁场数据并送入磁场模型。磁场模型从磁场数据中提入若干测点，并依次进行解算，将单个测点的定位结果作为测量值输入Kalman滤波器。Kalman滤波器接收测量值和进尺姿态数据，按照本发明设计的状态转移规则及滤波更新方程对测量结果进行滤波估计，最终输出高精度的定位结果。

本发明提出的改进SAGD算法，具体实施操作流程包括以下步骤：

1、在已有井中选取若干目标点作为相对导航定位的参考点，目标点间隔取为预定双井间距Δl的4倍左右。将磁传感器放置于第一个目标点上，钻头向该目标点方向钻进，在其由目标点后方2Δl位置进尺到前方2Δl位置后，磁传感器移至第二个目标点。依次照此法将磁传感器向第三、四、……个目标点移动，如图2所示。

2、当磁传感器放置于进尺为z₂的目标点，在钻头进尺z₁由z₂-2Δl变化至z₂+2Δl的过程中，磁传感器记录整个作用距离的磁场数据，并输入磁场模型。步骤3～5中磁场模型提取多个测点并给出各测点的定位结果，具体实施方式为：

3、磁场模型接收磁传感器的三轴磁场数据及控制设备给出的钻头进尺，记为(B_x(t),B_y(t),B_z(t),z(t))。

将采样数据按磁场变化周期进行分段，如图3所示。其中，第i段数据表示如下

$\begin{array}{c}{B}_{x,i}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ B_{y,i}\left(t\right)=B_y\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ B_{z,i}\left(t\right)=B_z\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\end{array}---\left(1\right)$

对B_x,i(t)B_y,i(t)B_z,i(t)分别取希尔伯特变换H[B_x,i(t)]H[B_y,i(t)]H[B_z,i(t)]，并求取角度

$\begin{array}{c}{\mathrm{\&phi;}}_{x,i}\left(t\right)=\arctan \frac{H\left[B_{x,i}\left(t\right)\right]}{B_{x,i}\left(t\right)}\\ {\mathrm{\&phi;}}_{y,i}\left(t\right)=\arctan \frac{H\left[B_{y,i}\left(t\right)\right]}{B_{y,i}\left(t\right)}\\ {\mathrm{\&phi;}}_{z,i}\left(t\right)=\arctan \frac{H\left[B_{z,i}\left(t\right)\right]}{B_{z,i}\left(t\right)}\end{array}---\left(2\right)$

以表示B_x,i(t)和B_y,i(t)在该周期内的平均幅值，则有

$\begin{array}{c}\left|{\stackrel{\&OverBar;}{B}}_x\left(i\right)\right|=\frac{1}{\mathrm{\&pi;}}\sqrt{{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,i}\left(t\right)\sin {\mathrm{\&phi;}}_{x,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2+{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,i}\left(t\right)\cos {\mathrm{\&phi;}}_{x,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2}\\ \left|{\stackrel{\&OverBar;}{B}}_y\left(i\right)\right|=\frac{1}{\mathrm{\&pi;}}\sqrt{{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,i}\left(t\right)\sin {\mathrm{\&phi;}}_{y,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2+{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,i}\left(t\right)\cos {\mathrm{\&phi;}}_{y,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2}\end{array}---\left(24\right)$

定义

$B_{\mathrm{mag}}\left(i\right)=\sqrt{{\left|{\stackrel{\&OverBar;}{B}}_x\left(i\right)\right|}^2+{\left|{\stackrel{\&OverBar;}{B}}_y\left(i\right)\right|}^2}---\left(25\right)$

搜索B_mag的最大值B_max＝max(B_mag)

对B_mag其进行归一化处理：

$B_{\mathrm{mag}}^{\&prime;}=\frac{B_{\mathrm{mag}}}{B_{\max }}---\left(26\right)$

4、磁场模型绘制如图4中粗线所示的归一化幅值与Δz/r关系曲线，

$B_{\mathrm{mag}}^{\&prime;}=\frac{\sqrt{2{(\mathrm{\&Delta;z}/r)}^4-2{(\mathrm{\&Delta;z}/r)}^2+5}}{\sqrt{5}{[1+{(\mathrm{\&Delta;z}/r)}^2]}^{\frac{5}{2}}}---\left(27\right)$

在纵坐标中选取若干归一化数值{η₁,η₂,...,η_M}作为测点，对于每一个η_m，在图3中找到满足公式(6)的两个解a_m和b_m；同时，在B′_mag中搜索找到两个归一化磁场幅值满足

B′_mag(i_m)＝B′_mag(j_m)＝η_m          (28)

并由脚标提取相应的进尺

$\begin{array}{c}{z}_{m,1}=z\left(t\right){|}_{t=\frac{2\mathrm{\&pi;}}{\mathrm{\&omega;}}{i}_m}\\ z_{m,2}=z\left(t\right){|}_{t=\frac{2\mathrm{\&pi;}}{\mathrm{\&omega;}}{j}_m}\end{array}---\left(29\right)$

进一步解算出相对距离

$r_m=\left|\frac{z_{m,1}-z_{m,2}}{a_m-b_m}\right|---\left(30\right)$

5、按测点数量对x、y方向的磁场分量按如下形式进行重新分段，

$\begin{array}{c}{B}_{x,m}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}{i}_m<\mathrm{\&omega;t}<2\mathrm{\&pi;}{i}_{m+1}},m=\mathrm{1,2},...,M-1\\ B_{x,m}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}{i}_m<\mathrm{\&omega;t}<2\mathrm{\&pi;}{j}_m},m=M\end{array}---\left(31\right)$

B_y,mB_z,m与B_x,m分段形式相似。同时，φ_x,m,φ_y,m,φ_z,m按公式2给出。

磁场模型给出如公式(11)所示的θ值计算公式，并将第m段数据输入磁场模型得到的相对角度记为θ_m

${\mathrm{\&theta;}}_m=a\tan 2[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,m}\left(t\right)\sin {\mathrm{\&phi;}}_{z,m}\left(t\right)d\left(\mathrm{\&omega;t}\right),-{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,m}\left(t\right)\sin {\mathrm{\&phi;}}_{z,m}\left(t\right)d\left(\mathrm{\&omega;t}\right)]---\left(32\right)$

步骤6～8中Kalman滤波器接收磁场模型数据及控制设备的姿态数据对r和θ进行Kalman滤波，定义

r_mθ_m：第m次滤波的测量值

由第m-1次滤波的一步预测值

第m次滤波估计值

p_r,mp_θ,m：第m次滤波方差

6、磁场模型输入第一组数据r₁θ₁，Kalman滤波器进行初始化： ${\hat{r}}_1=r_1{\widehat{\mathrm{\&theta;}}}_1={\mathrm{\&theta;}}_1$ p_r,1＝1p_θ,1＝1

7、对上一次滤波估计值实施状态转移，给出一步预测值。

$\begin{array}{c}{\hat{r}}_{m|m-1}=\sqrt{A^2+B^2}\\ {\widehat{\mathrm{\&theta;}}}_{m|m-1}=a\tan 2(B,A)\end{array},2\&le;m\&le;M---\left(33\right)$

其中，

上式中，

$\mathrm{\&Delta;}{z}_m=\left\{\begin{array}{c}{z}_{m+\mathrm{1,1}}-z_{m,1},m<M\\ z_{m,2}-z_{m,1},m=M\end{array}\right.$

8、将一步预测值和测量值输入Kalman滤波方程，对r和θ进行滤波更新，得到其滤波估计

$\begin{array}{c}{\hat{r}}_m={\hat{r}}_{m|m-1}+K_{r,m}(r_m-{\hat{r}}_{m|m-1})\\ {\widehat{\mathrm{\&theta;}}}_m={\widehat{\mathrm{\&theta;}}}_{m|m-1}+K_{r,m}({\mathrm{\&theta;}}_m-{\widehat{\mathrm{\&theta;}}}_{m|m-1})\end{array}---\left(34\right)$

其中，

$K_{r,m}=\frac{p_{r,m}}{p_{r,m}+{\mathrm{\&sigma;}}^2},K_{\mathrm{\&theta;},m}=\frac{p_{\mathrm{\&theta;},m}}{p_{\mathrm{\&theta;},m}+{\mathrm{\&sigma;}}^2}$

上式中，σ²为噪声方差。

重复7～8，经过M次滤波，输出和作为最终定位结果。

9、控制设备根据和调整下一个目标点作用距离内的钻头姿态，保持目标井在已有井正下方，并保持预定的相对距离。下面以一个实例介绍本发明方法

假设已有一长为100米的水平井，欲在其正下方5米处打一平行井，按步骤1每隔20米设一个标记点，共设5个标记点，进尺位置分别为10米、30米、50米、70米和90米。钻头钻进速度0.2m/s，转速120转/分，当钻头由0钻至20米距离时，磁传感器放置于第一个标记点，记录三轴磁场变化，如图5所示。

标记点1作用距离内持续约100秒，磁场共变化200个周期，将磁场数据分为200段，并计算200个离散点的总磁场幅值，并对其进行归一化。绘制如图4所示的曲线，然后取一组20个离散归一化值：

0.193,0.210,0.230,0.252,0.277,0.304,0.334,0.367,0.403,0.442,0.483,0.528,0.574,0.622,0.671,0.721,0.769,0.816,0.860,0.900，

并根据图4找到相应的a_m和b_m值。在200个总磁场归一化幅值中搜索20个分别与上述20个离散归一化值最接近的幅值，由幅值的脚标提取相应的进尺z_m,1和z_m,2。

由公式9解算出20个相对距离的测量值，由公式10计算方位角测量值，如图6和7中“□”型符号标示。按步骤6-8进行Kalman滤波，由图6和7所示Kalman滤波过程可以看出，经过若干次滤波后，滤波估计逐渐收敛到较为精确的结果，定位结果的方差随测点增加逐渐减小，经过15次左右的滤波后，方差减小到原来的5％以下。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (8)

1.一种基于Kalman滤波的改进SAGD算法，包括以下步骤：

①、在已有井中按一定间隔选取若干目标点作为相对导航定位的参考点，磁传感器将随钻进过程放置于不同目标点，划分各目标点作用距离；

②、在作用距离内，磁传感器记录钻进过程中的磁场数据，计算总磁场的归一化幅值；

③、选取若干离散的归一化幅值作为测点，各测点分别由SAGD算法计算得到多个相对距离r；

④、由上述测点中提取x、y方向的磁场分量，按照磁场模型解算出多个其双井相对的垂直偏离角度θ；

⑤、将多个测点的r和θ值依次输入Kalman滤波器，各组r和θ值按照状态转移规则进行一步预测；

⑥、对r和θ进行滤波更新，得到其滤波估计；

重复⑤-⑥，直至所有测点输入完毕，此时的滤波估计即为最终的定位结果。

2.根据权利要求1所述的基于Kalman滤波的改进SAGD算法，其特征在于，在步骤①中，目标点间隔取为预定双井间距Δl的4倍。

3.根据权利要求1所述的基于Kalman滤波的改进SAGD算法，其特征在于，在步骤①中，进尺为z₂的目标点，其作用距离为钻头进尺z₁由z₂-2Δl变化至z₂+2Δl的位置。

4.根据权利要求1所述的基于Kalman滤波的改进SAGD算法，其特征在于，步骤②中将连续磁场数据按周期进行离散化，计算各周期平均幅值，并以最大值为1进行归一化，具体方法是

在作用距离内，记录三轴磁场分量及对应的钻头进尺，即(B_x(t),B_y(t),B_z(t),z(t))对B_x(t)B_y(t)B_z(t)分别取希尔伯特变换H[B_x(t)]H[B_y(t)]H[B_z(t)]，并求取角度

$\begin{array}{c}{\mathrm{\&phi;}}_x\left(t\right)=\arctan \frac{H\left[B_x\left(t\right)\right]}{B_x\left(t\right)}\\ {\mathrm{\&phi;}}_y\left(t\right)=\arctan \frac{H\left[B_y\left(t\right)\right]}{B_y\left(t\right)}\\ {\mathrm{\&phi;}}_z\left(t\right)=\arctan \frac{H\left[B_z\left(t\right)\right]}{B_z\left(t\right)}\end{array}---\left(1\right)$

将采样数据按磁场变化周期进行分段，其中，第i段数据表示如下

$\begin{array}{c}{B}_{x,i}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ B_{y,i}\left(t\right)=B_y\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ {\mathrm{\&phi;}}_{x,i}\left(t\right)={\mathrm{\&phi;}}_x\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ {\mathrm{\&phi;}}_{y,i}\left(t\right)={\mathrm{\&phi;}}_y\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\\ {\mathrm{\&phi;}}_{z,i}\left(t\right)={\mathrm{\&phi;}}_z\left(t\right){|}_{2\mathrm{\&pi;}(i-1)<\mathrm{\&omega;t}<2\mathrm{\&pi;i}}\end{array}---\left(2\right)$

以表示B_x,i(t)和B_y,i(t)在该周期内的平均幅值，则有

$\begin{array}{c}\left|{\stackrel{\&OverBar;}{B}}_x\left(i\right)\right|=\frac{1}{\mathrm{\&pi;}}\sqrt{{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,i}\left(t\right)\sin {\mathrm{\&phi;}}_{x,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2+{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,i}\left(t\right)\cos {\mathrm{\&phi;}}_{x,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2}\\ \left|{\stackrel{\&OverBar;}{B}}_y\left(i\right)\right|=\frac{1}{\mathrm{\&pi;}}\sqrt{{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,i}\left(t\right)\sin {\mathrm{\&phi;}}_{y,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2+{\left[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,i}\left(t\right)\cos {\mathrm{\&phi;}}_{y,i}\left(t\right)d\left(\mathrm{\&omega;t}\right)\right]}^2}\end{array}---\left(3\right)$

定义

$B_{\mathrm{mag}}\left(i\right)=\sqrt{{\left|{\stackrel{\&OverBar;}{B}}_x\left(i\right)\right|}^2+{\left|{\stackrel{\&OverBar;}{B}}_y\left(i\right)\right|}^2}---\left(4\right)$

搜索B_mag的最大值B_max＝max(B_mag)

对B_mag其进行归一化处理：

$B_{\mathrm{mag}}^{\&prime;}=\frac{B_{\mathrm{mag}}}{B_{\max }}---\left(5\right).$

5.根据权利要求3所述的基于Kalman滤波的改进SAGD算法，其特征在于，选取一系列归一化数值{η₁,η₂,...,η_M}作为测点，分别由SAGD算法计算得到多个目标点处双井相对距离，具体过程为：

绘制SAGD算法描述的如下归一化磁场与Δz/r关系曲线

$B_{\mathrm{mag}}^{\&prime;}=\frac{\sqrt{2{(\mathrm{\&Delta;z}/r)}^4-2{(\mathrm{\&Delta;z}/r)}^2+5}}{\sqrt{5}{[1+{(\mathrm{\&Delta;z}/r)}^2]}^{\frac{5}{2}}}---\left(6\right)$

其中，Δz＝z₁-z₂，r为两口水平井实际垂直距离；

选取若干归一化数值{η₁,η₂,...,η_M}，对于每一个η_m，找到满足公式(6)的两个解a_m和b_m；同时，在步骤②得到的B′_mag中搜索找到两个归一化磁场满足

B′_mag(i_m)＝B′_mag(j_m)＝η_m         (7)

并由脚标提取相应的进尺

$\begin{array}{c}{z}_{m,1}=z\left(t\right){|}_{t=\frac{2\mathrm{\&pi;}}{\mathrm{\&omega;}}{i}_m}\\ z_{m,2}=z\left(t\right){|}_{t=\frac{2\mathrm{\&pi;}}{\mathrm{\&omega;}}{j}_m}\end{array}---\left(8\right)$

进一步解算出相对距离

$r_m=\left|\frac{z_{m,1}-z_{m,2}}{a_m-b_m}\right|---\left(9\right).$

6.根据权利要求1所述的基于Kalman滤波的改进SAGD算法，其特征在于，步骤④中计算双井垂直偏离角度的具体过程为：

重新将作用时段距离内的数据按如下形式分段：

$\begin{array}{c}{B}_{x,m}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}{i}_m<\mathrm{\&omega;t}<2\mathrm{\&pi;}{i}_{m+1}},m=\mathrm{1,2},...,M-1\\ B_{x,m}\left(t\right)=B_x\left(t\right){|}_{2\mathrm{\&pi;}{i}_m<\mathrm{\&omega;t}<2\mathrm{\&pi;}{j}_m},m=M\end{array}---\left(10\right)$

B_y,m,φ_x,m,φ_y,m,φ_z,m分别由与公式(2)相似形式进行分段给出；

由磁场模型得到如公式(11)所示的θ值计算公式，并将第m段数据输入磁场模型得到的相对角度记为θ_m

${\mathrm{\&theta;}}_m=a\tan 2[{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{x,m}\left(t\right)\sin {\mathrm{\&phi;}}_{z,m}\left(t\right)d\left(\mathrm{\&omega;t}\right),-{\&Integral;}_0^{2\mathrm{\&pi;}}{B}_{y,m}\left(t\right)\sin {\mathrm{\&phi;}}_{z,m}\left(t\right)d\left(\mathrm{\&omega;t}\right)]---\left(11\right).$

7.根据权利要求1所述的基于Kalman滤波的改进SAGD算法，其特征在于，对M组r和θ做Kalman滤波，并以第M次Kalman滤波估计作为最终定位结果输出。

8.根据权利要求1所述的基于Kalman滤波的改进SAGD算法，其特征在于，步骤⑤中所述的状态转移规则的具体方法为

定义

r_mθ_m：第m次滤波的测量值

由第m-1次滤波的一步预测值

第m次滤波估计值

p_r,mp_θ,m：第m次滤波方差

对于第m(2≤m≤M)次滤波前的一步预测，其状态转移规则为：

$\begin{array}{c}{\hat{r}}_{m|m-1}=\sqrt{A^2+B^2}\\ {\widehat{\mathrm{\&theta;}}}_{m|m-1}=a\tan 2(B,A)\end{array}---\left(12\right)$

其中，

上式中，

$\mathrm{\&Delta;}{z}_m=\left\{\begin{array}{c}{z}_{m+\mathrm{1,1}}-z_{m,1},m<M\\ z_{m,2}-z_{m,1},m=M\end{array}\right.$

当m＝1时， ${\hat{r}}_1=r_1{\widehat{\mathrm{\&theta;}}}_1={\mathrm{\&theta;}}_1.$