一种基于光纤陀螺测斜仪的测斜数据的事后处理系统

摘要

本发明公开了一种基于光纤陀螺测斜仪的测斜数据的事后处理系统，该系统首先用稳健方法的抗差估计准则将测斜数据中的粗差剔除，然后采用预测滤波原理将测斜数据中的毛刺去除，最后用三次B样条数据拟合方法将测斜数据进行平滑处理。本系统能够克服现有光纤陀螺测斜仪所测数据和相关石油标准在井眼轨迹还原过程中测斜数据不匹配的现状，为最终井眼轨迹的还原数据提供了可靠的测斜数据的事后处理。

说明书

技术领域

本发明涉及一种油井井眼轨迹还原的方法，更特别地说，是指一种 利用光纤陀螺测斜仪的测量得到的测斜数据进行井眼轨迹还原的事后处 理系统。

背景技术

《中华人民共和国石油天然气行业标准》，SY/T5435－2000，“定向 井轨道设计与轨迹绘图”。行业标准中规定了井眼轨迹还原以每25米～ 30米采集一个测斜数据点（简称为测点）作为还原出井眼轨迹的单位长 度。但是，光纤陀螺测斜仪实时测量得到的测斜数据更新时间较短（每 100毫秒更新一次），仪器井下运动速度为1～2米/秒，所以每30米有 300～600个测点。

光纤陀螺测斜仪能够测量的参数包括有：

井斜角（单位：度），井眼轴线上某点作井眼轴线的切线，该切线向 井眼前进方向延伸的部分称为井眼方向线，井眼方向线与重力线之间的 夹角就是井斜角。井斜角是钻井井眼轨道设计和实际的井眼轨迹计算中 的一个十分重要的参数。

方位角（单位：度），某测点处的井眼方向线投影到水平面上，称为 井眼方位线，或井斜方位线。以正北方向线为始边，顺时针旋转至井斜 方位线所转过的角度即方位角。

测深（单位：米），就是井深，指井下钻具的长度。

垂深（单位：米），就是测点的垂直深度。从测点做井口平面的垂线 就是垂深。

水平投影长度（单位：米），就是东北偏移，是指以井口为原点，井 眼轨迹在东北平面的投影，通常设定北为正向，北东轴之间为第一象限， 东南之间为第二象限，三四象限以此类推。

井眼曲率（单位：度/30米），也叫狗腿度与全角变化率。SY/T5313-93 《钻井工程术语》标准10.43条款中，对全角变化率定义为“单位井段 长度井眼轴线在三维空间的角度变化”，而单位井段长度取决于生产实 际中测斜需要。它既包含了井斜角的变化又包含着方位角的变化。井的 全角变化率过大会增加钻井和采油作业的困难，易导致井下事故。井眼 全角变化率引起的管杆弯曲,降低了管杆的自由度,将发生应力的集中释 放,使固定阀在下冲程中有开启过程,造成漏失,降低了泵效。

在井眼轨迹还原过程中，由于行业标准是以长度为计量单位，而光 纤陀螺测斜仪是以时间为单位，这必将造成实测测斜数据与需要还原的 井眼轨迹数据存在差异。如果抽取部分的实测测斜数据完成井眼轨迹还 原，将会引起还原的井眼轨迹存在很大偏差，对后续的地质资料解释带 来误导。

发明内容

为了解决光纤陀螺测斜仪的实测测斜数据与需要还原的井眼轨迹数 据存在的计量单位差异，本发明提出一种基于光纤陀螺测斜仪的测斜数 据的事后处理系统，该系统首先用稳健方法的抗差估计准则将测斜数据 中的粗差剔除，然后采用预测滤波原理将测斜数据中的毛刺去除，最后 用三次B样条数据拟合方法将测斜数据进行平滑处理。本发明采用300～ 600个测点对应于一个井眼轨迹还原点（简称为还原点），能够有效地抗 粗差，滤波及平滑测点。事后处理系统充分利用了大量的实时测斜数据， 不仅提高了仪器测量精度，同时使还原井眼轨迹的数据源更科学，合符 地质资料。

本发明的一种基于光纤陀螺测斜仪的测斜数据的事后处理系统，该 系统是为了解决光纤陀螺测斜仪的实测测斜数据与需要还原的井眼轨迹 数据存在的计量单位差异，该差异将造成还原的井眼轨迹存在很大偏差， 对后续的地质资料解释带来误导；其特征在于：该系统包括有稳健估计 的IGG迭代模块（10）、滤波毛刺去除模块（20）和三次样条平滑模块（30）；

在稳健估计的IGG迭代模块（10）中包括有下列处理步骤：

步骤101：稳健估计的IGG迭代模块10接收光纤陀螺测斜仪实时测 量输出的完整测斜数据DD；

步骤102：对所述完整测斜数据DD采用M估计的估计方法进行处理， 得到估计值DDA；

步骤103：对估计值DDA采用选权迭代法中的IGG方法按权函数 $\omega \left(u\right)=\left(\begin{array}{cc}1& \left|u\right|\le k_0\\ k_0/\left|u\right|& k_0\le \left|u\right|<k_1\\ 0& k_1\le \left|u\right|\end{array}\right)$进行迭代，得到抗差测斜数据DDB；

ω表示权重，u表示完整测斜数据与加权平均值的差，k₀表示第一淘 汰点，k₁表示第二淘汰点；

步骤104：统计抗差测斜数据DDB并进行加权处理，得到加权测斜数 据DDC；求加权测斜数据DDC的当前加权平均值A_n；当前加权平均值A_n中的n表示迭代次数；

步骤105：求当前加权平均值A_n与前一次加权平均值A_n-1的差值记为 ΔA_n，ΔA_n＝A_n-A_n-1；

步骤106：判断ΔA_n＜10^-4，即当前后两次迭代的加权平均值之差是否 小于10^-4，小于则停止迭代，输出M值估计数据DM；

若ΔA_n≥10^-4，重复步骤101至步骤105进行重复的抗差处理，直至 满足前后两次迭代的加权平均值之差小于10^-4；

在滤波毛刺去除模块（20）对接收到的的M值估计数据DM进行平滑 滤波处理，得到去毛刺后的测斜数据ADM；

在三次样条平滑模块（30）中对去毛刺后的测斜数据用三次B样条 进行平滑处理的步骤有：

步骤301：对接收到的所述去毛刺后测斜数据ADM进行滤波窗宽和 更新步长设置，得到输入数据样本；

所述的滤波窗宽设置为25，更新步长设置为10。

步骤302：将输入数据样本进行三次B样条拟合，输出拟合测斜数据 SDM；

步骤303：提取出拟合测斜数据SDM中的前10个数作为样条平滑后 的测斜数据估计参数；

提取出拟合测斜数据SDM中的后15个数更新输入数据样本；

步骤304：判断测斜数据长度DDL是否至结束，若到达测斜数据长度 的尾部，则结束对去毛刺后测斜数据的三次样条平滑处理，输出拟合测 斜数据SDM；若未到达测斜数据长度的尾部，则重复步骤302和步骤303， 直至到达测斜数据长度的尾部。

附图说明

图1是光纤陀螺测斜仪的结构图。

图2是光纤陀螺测斜仪的测量数据结构图。

图3A是本发明测斜数据事后处理系统的结构框图。

图3B是本发明测斜数据事后处理系统的流程图。

图3C是应用本发明系统后的轨迹复现示意图。

图4A是某型号光纤陀螺测斜仪测得的加入七处粗差的井斜数据波形 图。

图4B是对图4A进行了稳健估计的IGG迭代模块处理后的测斜数据 波形图。

图4C是对图4B部分进行放大的测斜数据波形图。

图5A为滤波之后两组数据对比图。

图5B为细节对比图。

图5C为平滑滤波处理之后的井斜数据图。

图6为经过样条平滑后的井斜数据图。

具体实施方式

下面将结合附图对本发明做进一步的详细说明。

参见图1所示，光纤陀螺测斜仪包括有中心处理器1、信号处理器2、 计数器3、光缆绞盘4、张力传感器5、井架6、光缆8和下井探管10组 成，中心处理器1与信号处理器2通过RS232接口连接，光缆绞盘4与 信号处理器2通过电缆线11连接，计数器3安装在光缆绞盘4上，光缆 绞盘4上缠绕有光缆8，井架6放置在被测井口7上方，井架6上安装有 张力传感器5，光缆8一端通过张力传感器5与下井探管10一端连接。

中心处理器1的控制流程为：系统上电后，首先进行自检，获得正 确的运行状态信息后进行测量参数设置；测量参数设置结束后测斜仪进 入初始对准状态，获得下井探管10机体坐标系12相对于当地地理坐标 系13的转换矩阵对准结束后进行测斜仪初始状态装订、输入当地经 纬度后进行下井测量状态；下井探管10到达指定深度后停止下行，并进 行速度误差归零处理，下井探管10进入上行测量状态；下井探管10到 达被测井口7后进行井口二次对准处理，实现井眼轨迹闭合测量，并对 测量数据进行存储、分析、绘图处理，并将结果输出。

中心处理器1由一公知计算机101和安装在所述计算机内的适用于 光纤陀螺测斜仪的捷联惯性组合测量单元，以及适用于下井探管10位置 设定的初始对准单元，以及测斜数据的事后处理系统。

所述捷联惯性组合测量单元是指对捷联惯性测量输出的姿态信息、 速度信息分别与磁测量输出的磁姿态信息、光缆绞盘运动的速度信息和 下井探管静止时的零速信息进行信息比较后，并对所述信息比较输出的 差值进行卡尔曼滤波处理后，用于估计捷联惯性组合的相关误差并进行 修正和补偿。经所述修正和补偿后的捷联惯性测量将输出的姿态信息给 计算机进行精确的井眼轨迹绘图及显示。

参见图2所示，所述捷联惯性组合测量单元包括捷联惯性测量102、 磁通门测量103、光缆运行速度测量104、卡尔曼滤波处理105、三个信 息比较模块（信息比较A106、信息比较B107和信息比较C108）和零 速修正模块（下井探管10静止时信息），下面将从各部分能够实现的功 能进行说明。

一、捷联惯性测量102

惯性测量分为平台式惯性测量和捷联式惯性测量，捷联式惯性测量 是将陀螺和加速度计直接固连在机体上的一种惯性测量方式。捷联式惯 性测量是用数学平台取代平台式惯性测量中的物理平台而进行的一种测 量模式。该捷联式惯性测量具有机械结构简单、尺寸小、相对造价较低， 但由于陀螺和加速度计直接固连在机体上，导致陀螺和加速度计的输入 动态范围较大，则在对使用合适的陀螺和加速度计提出了更高的要求， 此外选用数学平台代替物理平台对数据处理用的计算机要求较高。

在本发明中，捷联惯性测量102首先接收A）光纤陀螺输出的下井探 管机体坐标系O_bX_bY_bZ_b下的角速度信息；和B）加速度计输出的下井探管 机体坐标系O_bX_bY_bZ_b下的比力信息；然后对所述角速度信息和所述比力信 息采用航迹推算处理后，输出C）下井探管的速度信息V_x、V_y、V_z；和D） 下井探管姿态信息的方位角ψ、横滚角φ、倾斜角θ；最后，将E）所述 速度信息V_x、V_y、V_z输出给信息比较B单元和信息比较C单元；将F）所 述姿态信息输出给信息比较A单元。

卡尔曼滤波器105输出的所述状态变量X用于在线补偿所述捷联惯 性测量102的误差，其补偿后的表示下井探管姿态信息的输出为方位角 ψ₀、横滚角φ₀、倾斜角θ₀。并将这些姿态信息输出给计算机显示界面上 绘图显示（井眼轨迹绘制），以方便操作者实时监测。

垂直井深计算单元用于计算井眼的垂直井深，并将井深参数输出给 计算机显示界面显示，以方便操作者实时监测。垂直井深计算单元的井 深计算是通过将接收的缆长信息与下井探管10倾斜角θ₀的余弦相乘得到 被测点的垂直井深。

二、磁通门测量103

磁通门现象是一种普遍存在的电磁感应现象，它是利用高导磁材料 在传感器线圈系统中的磁饱和特性而测量磁场的仪器，是应用最广泛的 弱磁测量仪器之一，其测量范围为10^-12～10^-3T，这种仪器结构简单、体 积小、重量轻、功耗低、抗震性好、测量范围宽、灵敏度高、分辨率高、 稳定性好、适于测量恒定磁场或缓慢变化的磁场。

在本发明中，磁通门组件包括三个正交安置的磁通门，即X轴磁通 门、Y轴磁通门和Z轴磁通门。

在本发明中，磁通门测量103首先接收A）磁通门组件输出的下井 探管机体坐标系O_bX_bY_bZ_b下的地磁分量信息；和B）加速度计输出的下 井探管机体坐标系O_bX_bY_bZ_b下的比力信息；然后对所述地磁分量信息和 所述比力信息采用磁测量计算，并将计算出的相对于磁北的方位角转换 为相对于真北的方位角ψ_c后，输出下井探管姿态信息的方位角ψ_c、横 滚角φ_c、倾斜角θ_c；最后将C）所述姿态信息输出给信息比较A单元。

三、光缆运行速度104

光缆8的一端连接在光缆绞盘4上，另一端固定在下井探管10上。 在本发明中，对下井探管10在井中运行的速度可通过计数器3测量得到 的光缆长度在单位时间内光缆长度增量得到下井探管运行时在其机体坐 标系O_bX_bY_bZ_b下的速度信息。

全光纤数字测斜仪在测量过程中涉及两个坐标系，一个是当地地理 坐标系O_tX_tY_tZ_t13，一个是下井探管机体坐标系O_bX_bY_bZ_b12。这两个坐 标系存在的坐标转换关系为$\left(\begin{array}{c}{x}_b\\ y_b\\ z_b\end{array}\right)=C_t^b\left(\begin{array}{c}{x}_t\\ y_t\\ z_t\end{array}\right),$表示坐标转换矩阵。通过这 个坐标转换矩阵可将下井探管机体坐标系O_bX_bY_bZ_b下的光缆运行速度 信息转换为当地地理坐标系O_tX_tY_tZ_t下的速度信息V_xc、V_yc、V_zc，所述速 度信息V_xc、V_yc、V_zc输出给信息比较B单元进行速度差值比较。

下井探管机体坐标系O_bX_bY_bZ_b与当地地理坐标系O_tX_tY_tZ_t之间的转换 矩阵为 $C_t^b=\left(\begin{array}{ccc}\cos \psi \cos \phi +\sin \psi \sin \theta \sin \phi & -\sin \psi \cos \phi +\cos \psi \sin \theta \sin \phi & -\cos \theta \sin \theta \\ \sin \psi \cos \theta & \cos \psi \cos \theta & \sin \theta \\ \cos \psi \sin \theta -\sin \psi \sin \theta \cos \phi & -\sin \psi \sin \phi -\cos \psi \sin \theta \cos \phi & \cos \theta \cos \phi \end{array}\right),$这个 矩阵中的角度分别有下井探管姿态信息的方位角ψ、横滚角φ、倾斜角θ。

四、卡尔曼滤波器105

卡尔曼滤波器105接收

A）信息比较A单元输出的姿态差值$Z_{\mathrm{KA}}=\left(\begin{array}{c}\psi -{\psi }_c\\ \theta -{\theta }_c\\ \phi -{\phi }_c\end{array}\right),$和

B）信息比较B单元输出的速度差值$Z_{\mathrm{KB}}=\left(\begin{array}{c}{V}_x-V_{\mathrm{xc}}\\ V_y-V_{\mathrm{yc}}\\ V_z-V_{\mathrm{zc}}\end{array}\right)$和

C）信息比较C单元输出的速度差值$Z_{\mathrm{KC}}=\left(\begin{array}{c}{V}_x-0\\ V_y-0\\ V_z-0\end{array}\right),$利用离散型卡尔曼滤波 器对接收信息进行数据融合实现对状态变量 X＝[δr_x,δr_y,δr_z,δv_x,δv_y,δv_z，η_x，η_y，η_z,δf_x,δf_y,δf_z,δω_x,δω_y,δω_z] 的进行最优估计，根据估计的状态变量对所述捷联惯性测量进行在线误 差补偿，输出补偿后的下井探管姿态信息$\left(\begin{array}{c}{\psi }_0\\ {\theta }_0\\ {\phi }_0\end{array}\right)=\left(\begin{array}{c}\psi -{\eta }_x\\ \theta -{\eta }_y\\ \phi -{\eta }_z\end{array}\right)$给所述计算机显示 界面101。

式中，δr_x表示当地地理坐标系O_tX_tY_tZ_t下的位置误差X轴分量，δr_y表 示当地地理坐标系O_tX_tY_tZ_t下的位置误差Y轴分量，δr_z表示当地地理坐标 系O_tX_tY_tZ_t下的位置误差Z轴分量，δv_x表示当地地理坐标系O_tX_tY_tZ_t下的 速度误差X轴分量，δv_y表示当地地理坐标系O_tX_tY_tZ_t下的速度误差Y轴分 量，δv_z表示当地地理坐标系O_tX_tY_tZ_t下的速度误差Z轴分量，η_x表示方位 角ψ误差，η_y表示倾斜角θ误差，η_z表示横滚角φ误差，δf_x表示X轴加速 度计误差，δf_y表示Y轴加速度计误差，δf_z表示Z轴加速度计误差，δω_x表 示X轴光纤陀螺误差，δω_y表示Y轴光纤陀螺误差，δω_z表示Z轴光纤陀 螺误差。

所述的卡尔曼滤波器105的滤波处理采用离散型卡尔曼滤波基本方 程，其各状态变量初值均设为零。

五、信息比较

在本发明中，信息比较由信息比较A单元106、信息比较B单元107 和信息比较C单元108组成，其中，

信息比较A单元106，用于完成对所述捷联惯性测量102输出的下井 探管姿态信息的方位角ψ、横滚角φ、倾斜角θ与所述磁测量103输出的 下井探管姿态信息的方位角ψ_c、横滚角φ_c、倾斜角θ_c相减输出姿态差值 $Z_{\mathrm{KA}}=\left(\begin{array}{c}\psi -{\psi }_c\\ \theta -{\theta }_c\\ \phi -{\phi }_c\end{array}\right)$给卡尔曼滤波器105；

信息比较B单元107，用于完成对所述捷联惯性测量102输出的下井 探管的速度信息V_x、V_y、V_z与所述光缆运动速度测量104输出的速度信 息V_xc、V_yc、V_zc相减输出速度差值$Z_{\mathrm{KB}}=\left(\begin{array}{c}{V}_x-V_{\mathrm{xc}}\\ V_y-V_{\mathrm{yc}}\\ V_z-V_{\mathrm{zc}}\end{array}\right)$给卡尔曼滤波器105；

信息比较C单元108，用于完成对所述捷联惯性测量102输出的下井 探管的速度信息V_x、V_y、V_z与所述下井探管10静止量输出的零速度信息 相减输出速度差值$Z_{\mathrm{KC}}=\left(\begin{array}{c}{V}_x-0\\ V_y-0\\ V_z-0\end{array}\right)$给卡尔曼滤波器105。

六、零速修正

零速修正是指全光纤数字测斜仪中下井探管处于静止位置时，其速 度为零的现象来修正捷联惯性测量输出的速度误差。该“零速修正”通 过设置的一按键实现，参见图2A所示。采用所述零速修正模式时，下井 探管实际的运行速度为零作为输出给信息比较C单元的比较信息。

本发明捷联惯性组合测量单元可以提供四种组合测量工作模式，各 个模式通过在界面中按下相应的键实现。其中，

模式一：捷联惯性测量＋缆长  （测量模式中按键设定为①）

该组合测量模式用于常规有线测井的快速连续测量，尤其适用于套 管井或有磁性物质干扰的井眼轨迹测量。

组合测量所需的信息流程为：第一路信息，捷联惯性测量首先接收A） 光纤陀螺输出的下井探管机体坐标系O_bX_bY_bZ_b下的角速度信息；和B）加 速度计输出的下井探管机体坐标系O_bX_bY_bZ_b下的比力信息；然后对所述角 速度信息和所述比力信息采用航迹推算处理后，输出C）下井探管的速度 信息V_x、V_y、V_z；最后将所述速度信息V_x、V_y、V_z输出给信息比较B单元。 第二路信息，通过计数器对光缆缆长的记录，再经坐标转换获得的表示 井下探管的运行速度信息V_xc、V_yc、V_zc输出给信息比较B单元。在信息比 较B单元中第一路信息与第二路信息进行速度差值比较后输出给卡尔曼 滤波器处理。第三路信息，经卡尔曼滤波输出的所述状态变量X用于在 线补偿所述捷联惯性测量的误差，其补偿后表示下井探管姿态信息的输 出为方位角ψ₀、横滚角φ₀、倾斜角θ₀。所述姿态信息输出给计算机显示 界面绘图显示，以方便操作者实时监测。第四路信息，所述姿态信息输 出给垂直井深计算单元进行垂直井深计算，并输出井深参数给计算机显 示界面显示，以方便操作者实时监测。

模式二：捷联惯性测量＋零速修正（测量模式中按键设定为②）

该组合测量模式用于不要求连续测量的测井过程中，可用在随钻测 井（MWD）过程中，实时监测钻头的钻进方向。

组合测量所需的信息流程为：第一路信息，捷联惯性测量首先接收A） 光纤陀螺输出的下井探管机体坐标系O_bX_bY_bZ_b下的角速度信息；和B）加 速度计输出的下井探管机体坐标系O_bX_bY_bZ_b下的比力信息；然后对所述角 速度信息和所述比力信息采用航迹推算处理后，输出C）下井探管的速度 信息V_x、V_y、V_z；最后将所述速度信息V_x、V_y、V_z输出给信息比较C单元。 第二路信息，下井探管实际的运行速度为零作为输出给信息比较C单元 的比较信息。在信息比较C单元中第一路信息与第二路信息进行速度差 值比较后输出给卡尔曼滤波器处理。第三路信息，经卡尔曼滤波输出的 所述状态变量X用于在线补偿所述捷联惯性测量的误差，其补偿后表示 下井探管姿态信息的输出为方位角ψ₀、横滚角φ₀、倾斜角θ₀。所述姿态 信息输出给计算机显示界面绘图显示，以方便操作者实时监测。第四路 信息，所述姿态信息输出给垂直井深计算单元进行垂直井深计算，并输 出井深参数给计算机显示界面显示，以方便操作者实时监测。

模式三：捷联惯性测量＋磁测量（测量模式中按键设定为③）

该组合测量模式用于无磁环境干扰的测井作业，实现高精度连续测 量。

组合测量所需的信息流程为：第一路信息，捷联惯性测量首先接收A） 光纤陀螺输出的下井探管机体坐标系O_bX_bY_bZ_b下的角速度信息；和B）加 速度计输出的下井探管机体坐标系O_bX_bY_bZ_b下的比力信息；然后对所述角 速度信息和所述比力信息采用航迹推算处理后，输出C）表示下井探管姿 态信息的方位角ψ、横滚角φ、倾斜角θ；最后将所述姿态信息输出给信 息比较A单元。第二路信息，磁通门测量首先接收A）磁通门组件输出的 下井探管机体坐标系O_bX_bY_bZ_b下的地磁分量信息；和B）加速度计输出的 下井探管机体坐标系O_bX_bY_bZ_b下的比力信息；然后对所述地磁分量信息和 所述比力信息采用磁测量计算，并将计算出的相对于磁北的方位角转换 为相对于真北的方位角ψ_c后，输出表示下井探管姿态信息的方位角ψ_c、 横滚角φ_c、倾斜角θ_c；最后将所述姿态信息输出给信息比较A单元。在信 息比较A单元中第一路信息与第二路信息进行姿态差值比较后输出给卡 尔曼滤波器处理。第三路信息，经卡尔曼滤波输出的所述状态变量X用 于在线补偿所述捷联惯性测量的误差，其补偿后表示下井探管姿态信息 的输出为方位角ψ₀、横滚角φ₀、倾斜角θ₀。所述姿态信息输出给计算机 显示界面绘图显示，以方便操作者实时监测。第四路信息，所述姿态信 息输出给垂直井深计算单元进行垂直井深计算，并输出井深参数给计算 机显示界面显示，以方便操作者实时监测。

模式四：捷联惯性测量＋缆长＋零速修正＋磁测量（测量模式中 按键设定为④）

该组合测量模式用于有缆长辅助且无连续测量要求的各类测井环 境，通过多种数据融合技术的应用，可提供精确估计和补偿捷联惯性的 误差，能提供高精度的测量结果。

在本发明中，所述卡尔曼滤波器对状态变量进行实时在线估计。四 种不同的相关测量模式在卡尔曼滤波处理中输出的所述状态变量X是相 同的，但使用的观测量是不同的。

当选取模式一时，其相应的量测阵 $H_1=\left(\begin{array}{ccccccccccccccc}0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\end{array}\right),$则补偿后表示下井探 管姿态信息的输出为$\left(\begin{array}{c}{\psi }_0\\ {\theta }_0\\ {\phi }_0\end{array}\right)=\left(\begin{array}{c}\psi -{\eta }_x\\ \theta -{\eta }_y\\ \phi -{\eta }_z\end{array}\right).$

当选取模式二时，其相应的量测阵 $H_{\mathrm{2,3}}=\left(\begin{array}{ccccccccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right),$则补偿后表示下井探 管姿态信息的输出为$\left(\begin{array}{c}{\psi }_0\\ {\theta }_0\\ {\phi }_0\end{array}\right)=\left(\begin{array}{c}\psi -{\eta }_x\\ \theta -{\eta }_y\\ \phi -{\eta }_z\end{array}\right).$本发明中，模式三与模式二的量测阵 是相同的，则输出也相同，即模式三与模式二都是以速度差值作为观测 量。

当选取模式四时，其相应的量测阵 $H_4=\left(\begin{array}{cccccccccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\end{array}\right),$则补偿后表示下井 探管姿态信息的输出为$\left(\begin{array}{c}{\psi }_0\\ {\theta }_0\\ {\phi }_0\end{array}\right)=\left(\begin{array}{c}\psi -{\eta }_x\\ \theta -{\eta }_y\\ \phi -{\eta }_z\end{array}\right).$

在本发明中，信号处理器2用于实现井下数据的接收的解码和井上 下传控制指令及光缆缆长数据的编码和发送。信号处理器2由FPGA逻辑 处理器、功率放大电路和光电转换电路组成，FPGA逻辑处理器接收由中 心处理器1输出的光缆缆长L（缆长是指从被测井口7处到下井探管10 所处当前位置时光缆通过光缆绞盘4的长度，其长度由计数器3读取， 且输出给信号处理器2）和下井探管10运行状态指令，并将所述信息按 照曼彻斯特编码格式进行编码处理后输出给功率放大电路；功率放大电 路对所述编码信息进行放大、滤波处理输出给光电转换电路，在光电转 换电路中将电信号转换为光强信号，所述光强信号通过光缆耦合器耦合 进入光缆8（七芯铠装光缆）传递给下井探管10。

参见图3A、图3B所示，本发明为了解决光纤陀螺测斜仪的实测测斜 数据与需要还原的井眼轨迹数据存在的计量单位差异，该差异将造成还 原的井眼轨迹存在很大偏差，对后续的地质资料解释带来误导，因此。 提出一种内嵌在中心处理器1（或者称上位机）中的测斜数据的事后处理 系统（采用计算机编程语言实现）。所述测斜数据的事后处理系统包括有 稳健估计的IGG迭代模块10、滤波毛刺去除模块20和三次样条平滑模块 30。

光纤陀螺测斜仪用于测量表征井眼轨迹的参数（至少包括井斜角、 方位角、工具面角等），多组这样的井眼轨迹参数形成井眼轨迹的完整测 斜数据记为DD，井眼轨迹参数的组数个数称为测斜数据长度记为DDL。

本测斜数据处理系统首先接收A)光纤陀螺油井连续测斜仪的输出也 就是原始测斜数据DD，代表所有原始测斜数据组成的集合，包含了7种 数据：L斜深,θ井斜角,ψ方位角,γ为工具面角，H垂深，V速度，Nl 北向偏移，S1东向偏移，ρ造斜率；

将原始测斜数据DD经过稳健估计思想中的M估计方法IGG为权函数 的选权迭代法进行抗差处理之后输出B)D_M为抗差处理后的数据，即将数 据内的粗差剔除。粗差就是在测斜过程中，由于井下环境的干扰和传感 器的偶然跳动等不确定因素的影响，实测数据中可能出现不合理的跳点。 经过抗差处理后的数据在细节图上还有很多杂散数据，输入至预测滤波 系统，输出C)D_F为预测滤波后的数据。

所述杂散数据就是毛刺，定义为在采样时偏离正常数据的任何小的 跳变。

然后将上述数据经过三次B样条拟合输出D)D_S为用三次B样条数据 拟合方法将预处理之后的数据，此时，经过处理的测斜数据仍然是随时 间改变而改变的曲线。

对上述数据根据测斜仪速度与位移的关系进行测斜数据采样输出E) AL为基于深度的采样数据。

所述测斜仪速度与位移的关系即S＝V×ΔT。经过抗差与平滑之后的 测斜数据仍然是时间的函数，为了和国家标准匹配25～30m采集一个数 据复现井眼轨迹，那么S=25～30m，测斜仪的速度V是可测量，从而计算 出采集数据时间间隔△T。如果测斜仪以1m/s的速度测量数据，那么我 们的采样的时间间隔为25～30s。采集到的样本深度间隔为25～30m。

本系统能够克服了现有仪器所测数据和相关石油标准在井眼轨迹还 原过程中数据不匹配的情况，为最终井眼轨迹的还原提供可靠的数据源。

从图3C中可以看出：对原始测斜数据经过抗差处理、预测滤波、三 次B样条拟合之后，得到与真实井眼轨迹相似的轨迹曲线。在这条曲线 上可以根据距离进行采样，每25～30m采集一个数据，对应时间约为25～ 30s，从而得到符合SY/T5435－2000的中华人民共和国石油天然气行业 标准的“定向井轨道设计与轨迹绘图”。在工程应用中，当实时采集到测 斜数据时，通过标准中公开的表中的关系将采集数据填入，所需数据通 过查表方式能够得到。

（一）稳健估计的IGG迭代模块10

步骤101：稳健估计的IGG迭代模块10接收光纤陀螺测斜仪实时测 量输出的完整测斜数据DD；

步骤102：对所述完整测斜数据DD采用M估计的估计方法进行处理， 得到估计值DDA；

在本发明中，所述M估计的估计方法参考《测量数据处理理论与方 法》，邱卫宁等人编著，武汉大学出版社，2008年6月。

在本发明中，所述M估计的估计方法也可以是以n个观测样本 L₁,L₂,…,L_i,…,L_n，X为待估参数。任意一个观测样本L_i的分布密度记为 表示采集到的样本值，由极大似然估计准则 $\underset{i=1}{\overset{n}{\Sigma }}\ln f(L_i,\hat{X})=\max ,$求导得到其中， ρ(x)＝-lnf(x)。在这里当选择出能够满足稳 健化要求的函数ρ，利用估算参数X。

步骤103：对估计值DDA采用选权迭代法中的IGG方法按权函数 $\omega \left(u\right)=\left(\begin{array}{cc}1& \left|u\right|\le k_0\\ k_0/\left|u\right|& k_0\le \left|u\right|<k_1\\ 0& k_1\le \left|u\right|\end{array}\right)$进行迭代，得到抗差测斜数据DDB；

ω表示权重，u表示完整测斜数据与加权平均值的差，k₀表示第一淘 汰点，k₁表示第二淘汰点。

在本发明中，由于本系统数据更新间隔为100ms，在井下运行速度为 1～2m/s，根据相关石油标准同时考虑到数据处理的稳健性，将数据处理 容量定为25个，但每次更新10个数据，即用本次得到校正的15个数据 加上10个待处理数据作为下次处理的数据。同时起始端先处理50个数 以保证初始信息的准确性。设定迭代步长为25及更新步长为10。

步骤104：统计抗差测斜数据DDB并进行加权处理，得到加权测斜数 据DDC；求加权测斜数据DDC的当前加权平均值A_n；当前加权平均值A_n中的n表示迭代次数；

步骤105：求当前加权平均值A_n与前一次加权平均值A_n-1的差值记为 ΔA_n，ΔA_n＝A_n-A_n-1；

步骤106：判断ΔA_n＜10^-4，即当前后两次迭代的加权平均值之差是否 小于10^-4，小于则停止迭代，输出M值估计数据DM；

若ΔA_n≥10^-4，重复步骤101至步骤105进行重复的抗差处理，直至 满足前后两次迭代的加权平均值之差小于10^-4。

图4A是以苏73－12x油井测斜仪的上测过程部分数据。横坐标为数 据个数（500个），纵坐标为井斜角α。为保证实验数据的全面性，增加7 处粗差，在粗差出现数量上包含了单值型粗差，成片型粗差；在粗差数 值上包括了正向型粗差、负向型粗差以及正负交替型粗差。

图4B是经稳健估计的IGG迭代模块10处理后的测斜数据。数据经 过抗差处理后，毛刺少了，但从细节看还是有很多毛刺，如图4C所示。

（二）滤波毛刺去除模块20

工程实际中，采集10～15个的测斜数据对应时间为1～2秒，对应 光纤陀螺测斜仪则运动了1～2米，在如此短的距离，井斜角和方位角的 变化不仅没有实际的物理意义，而且对后期平滑带来影响，所以应将此 类数据作为噪声除去。

滤波毛刺去除模块20对接收到的的M值估计数据DM（即加权平均 值）进行平滑滤波处理，得到去毛刺后的测斜数据ADM。

图5A为滤波之后两组数据对比图，图5B为细节对比图（为方便观 察，将原始数据错开一个极小量）。图5C为平滑滤波处理之后的井斜数 据图。

在本发明中，平滑滤波处理可以参考《最优估计的计算方法》的固 定区间最优平滑估计章节的内容，史忠科著，科学出版社出版，2001年 5月。

（三）三次样条平滑模块30

由于井眼轨迹在空间里是三维连续曲线，可用三次参数样条函数描 述。三次参数样条函数是根据一些有效特征点进行曲线拟合以达到对未 知特征进行插值预测。在本发明中，采用三次参数样条函数是根据一些 井眼轨迹的有效特征点进行曲线拟合，以达到对未知特征进行插值预测， 实现对缺失数据修补的目的。

在本发明中，对去毛刺后的测斜数据用三次B样条进行平滑处理， 得到样条测斜数据。

步骤301：对接收到的所述去毛刺后测斜数据ADM进行滤波窗宽和 更新步长设置，得到输入数据样本；

在本发明中，所述的滤波窗宽设置为25，更新步长设置为10。

由于光纤陀螺测斜仪数据更新间隔为100ms，井下运行速度为1～ 2m/s，根据相关石油标准，考虑到数据处理的可靠性，将数据处理容量 定为25个，但每次更新15个数据，即用本次得到校正的15个数据加上 10个待处理数据作为下次处理的数据。

步骤302：将输入数据样本进行三次B样条拟合，输出拟合测斜数据 SDM；

步骤303：提取出拟合测斜数据SDM中的前10个数作为样条平滑后 的测斜数据估计参数；

提取出拟合测斜数据SDM中的后15个数更新输入数据样本；

步骤304：判断测斜数据长度DDL是否至结束，若到达测斜数据长度 的尾部，则结束对去毛刺后测斜数据的三次样条平滑处理，输出拟合测 斜数据SDM；若未到达测斜数据长度的尾部，则重复步骤302和步骤303， 直至到达测斜数据长度的尾部。

在本发明中，采用三次参数样条函数是根据一些井眼轨迹的有效特 征点进行曲线拟合，以达到对未知特征进行插值预测，实现对缺失数据 修补的目的。将经过前面两步处理后的测斜数据用三次B样条平滑，使 得井眼轨迹在空间里是三维一个连续曲线。而在样条节点处，三次样条 连续，并有连续的一、二阶导数，这与工程景象吻合。如图6所示的经 过样条平滑后的井斜数据图。