一种油井测斜数据抗差处理算法研究
李毅，张春熹**
基金项目：国家863 高技术研究发展计划项目“高精度光纤惯性组合测斜仪技术研究”(2007AA06Z101)
作者简介：李毅，（1988-），男，硕士生，光纤陀螺测斜仪测斜数据事后处理算法。
通信联系人：张春熹，（1965-），男，教授，博士生导师，主要研究方向：微弱信号检测. E-mail:
zhangchunxi8231@126.com
（北京航空航天大学光电技术研究所，北京 100083）
5 摘要：为了提高井眼轨迹工程参数的测量精度，实现井眼轨迹科学可靠的还原，对四种基于
选权迭代法的稳健估计方法进行了详细分析和比较，讨论了四种权函数抗差性与临界值的关
系。通过对几种算法的理论分析，并利用VC 环境的MFC 平台完成对四种算法的仿真，用实
测数据详细对比了四种算法的抗差效果和运算效率，得到了针对井眼轨迹工程参数测量数据
的有效抗差处理方案，使完全正确处理500 个测斜数据的时间缩短至469 毫秒，有效提高了
10 工作效率，同时也为其他测量数据的预处理提供了参考。
关键词：数据处理；抗差算法；稳健估计；运算效率
 0 引言
30 在油井测斜过程中，井眼轨迹位置参数的确定与描述为复现井眼轨迹提供了最直接有效
的方法。为了能够精确地描述井眼轨迹的空间形状，必须要得到足够多的合理数据。但是，
在测斜过程中, 由于井下环境的干扰和传感器的偶然跳动等不确定因素的影响, 实测数据中
可能出现不合理的跳点，称之为粗差或野值。含有粗差的测斜数据不仅会降低数据的可信度，
而且在其后续处理中会严重影响处理结果的质量，甚至得出错误的结论。因此，需要对测斜
35 数据粗差进行识别、剔除与补正。由于测斜数据是时间序列，可以采用广义卡尔曼滤波、偏
最小二乘法、建立神经网络模型等方法预测下一个观测值[1,2,3]，并求其残差来识别粗差, 也
可以利用小波分析等方法来处理[4]。但上述方法需要预先建立动力学方程组或计算复杂度过
大，不适合在测斜数据预处理阶段使用。针对以上情况，本文以基于光纤陀螺的惯性技术测
斜仪为工程背景，对所测得的井眼轨迹数据进行分析，运用稳健估计的思想，通过探讨和对
40 比测量界常用的Huber 法、Hampel 法、丹麦法以及IGG 法的处理效果和运行效率，并在
VC 环境下进行仿真，据此确立一种适用于测井数据抗粗差的稳健估计方法，该方法不仅能
 有效降低粗差对数据的干扰，还能实现误差补偿，从而正确体现测斜仪器的精度。
1 稳健估计及选权迭代原理
1.1 稳健估计的基本思想
45 稳健估计（Robust Estimation），在测量中也称为抗差估计，它是在针对最小二乘法抗
粗差的干扰差这一缺陷提出来的，在粗差不可避免的情况下，选择适当的估计方法使未知量
估计值尽可能减免粗差的影响，得出最佳估值[5]。抗差估计的原则是要充分地利用有效信息，
限制利用不可用信息，排除有害信息。但是往往由于事先不大可能准确知道观测数据中有效
信息和有害信息所占的比例以及它们具体包含在哪些观测值中，从抗差估计的主要目标着眼
50 要冒险损失一些效率的风险，去获得较可靠的、具有实际意义的和较有效的估值。
设观测样本L1，L2，…，Ln，X 为待估参数，观测值i l 的分布密度为( , ˆ) i f l X ，由极大似
然估计准则为
 选择能够满足稳健化要求的函数ρ （或ϕ ），利用（2）式或（3）式估算参数X，这
就是M 估计[6]。
1.2 基于等权独立观测的选权迭代法
M 估计的估计方法有许多种，在测量平差中应用最广泛，易于程序实现的是选权迭代
法。设独立观测值为L，未知参数向量为ˆX 60 ，误差方程及权阵为
 式中： i b 为1×t 系数向量。
考虑误差方程，M 估计的函数
( , ˆ) iρ l X 可表述为
( , ˆ) ( ) i i ρl X =ρv （5）
由于待处理的每组测量数据都是同一仪器在某段特定时间所测得的数据，所以本文采用
等权观测的方法。按M 估计极大似然估计准则并取ρ 函数为式（5） ，则有
 称为稳健权矩阵，其元素i ω 称为稳健权因子，是相应残差i v 的 函数。将误差方程（7）
75 带入所得M 估计的法方程式为
BTWBXˆ = BTWl （8）
当选定ρ 函数后，稳健权矩阵W 可以确定，但是i ω 是i v 的函数，故需要对权进行迭
代求解。其迭代计算步骤为：
1）列误差方程，令权因子初值为1，即W＝I。
2）解算法方程（8），得出参数ˆX 和残差V 的第一次估值：
80 Xˆ (1) =(BTWB)−1BTWL
V(1) =BXˆ (1) −L
3）由
V (1) 按( i)
i
i
v
v
ϕ
=ω 确定各观测值新的权因子
Xˆ (2) =(BTW(1)B)−1BTW(1)L
V(2) =BXˆ (2) −L
4）由
V (2) 构造新的等价权
85 W(2) ，在解算法方程，类似迭代计算，直至结果收敛。即得：
Xˆ(k) =(BTW(k−1)B)−1BTW(k−1)L （9）
V(k) =BXˆ (k−1) −L （10）
从迭代过程可以看出，随着ρ 函数选取的不同，构成了权函数的多种不同形式，但权
函数总在一个平差过程中随改正数变化，其中i ω 与i v 大小成反比，因此经过多次迭代，可
以使含有粗差的观测值的权函数为零（或接近零），使其在平差中不起作用，从而在变权过
程中实现参数估计的稳健性[7]。
  90 1.3 几种典型的选权迭代法
(1) Huber 函数
1,
( )
,
u k
u k u k
u
ω
⎧ ≤
⎪ = ⎨ > ⎪⎩
（11）
式中：u=v/s；s 在u≤k区间取标准差，u 在u >k 区间s 取中位绝对差。
(2) 丹麦函数
2
1 ,
( )
exp(1 ( / ) ) ,
u c
u
u c u c
ω
= ⎧⎨⎪ ≤
⎩⎪ − >
（12）
(3) IGG 函数
0
0 0 1
1
1 ,
( ) / ,
0,
u k
u k u k u k
k u
ω
⎧ ≤
⎪
= ≤ < ⎨⎪
⎩ ≤
（13）
95 (4) Hampel 函数
1 ,
,
( )
,
( )
0 ,
u a
a a u b
u
u
c u
a b u c
c b u
c u
ω
⎧ ≤
⎪⎪
< ≤
⎪⎪
= ⎨⎪ − < ≤
⎪ −
⎪
< ⎪⎩
i
（14）
式中：u=v/MAD，a、b、c 为调制系数
2 权函数的理论分析
由（11）~（14）式可以看出，四种权函数的共同特点是将权函数的取值分为若干个区
间，根据自变量不同的取值决定其所占权重。不同的区间划分方式必将影响到权函数的适用
100 范围以及数据处理效果。
从式（11）的Huber 权函数可以看出，其权函数是单调的，包含正常域和可疑域，对正
常域的观测值采用最小二乘估计，对其它观测值一律降权，正常域临界值为k，当改正数大
于k 时, 其所占权重与其成反比, 误差愈大, 对应的权重愈小, 与此相应该观测值对参数估
计的影响也愈小。
105 式（12）的丹麦权函数和Huber 权函数类似，也包含正常域和可疑域，对正常域的观测
值采用最小二乘估计，对其它观测值一律降权，不同的是其降权函数呈指数变化。
式（13）的IGG 方案是基于测量误差的有界性提出来的，它对测量抗差估计比较有效[8]。
相比前两种方法增加了淘汰域，对位于淘汰域的数据直接将其舍去。根据不同的工程数据的
特点，采用不同的阈值范围，通过将超出误差范围的异常值直接剔除，能显著提高计算速度。
110 与前述三者不同的是，式（14）的Hampel 权函数分为四段，不仅有正常域、可疑域和
淘汰域，而且将可疑域分成了两段，通过对可疑数据的权重做多范围的调整，其稳健性是四
种函数中最好的，但是，由于迭代次数的增加牺牲了计算效率。
 3 算例
3.1 仿真数据
115 现以苏73-12x 井的上测过程部分数据为例（如图1 所示），图中显示该井部分井斜数
据，该组数据满足实验数据的全面性，7 处粗差包括了仪器数据中可能出现的所有情况，可
以有效评判算法的正确性。在粗差出现数量上包含了单值型粗差，成片型粗差；在粗差数值
上包括了正向型粗差、负向型粗差以及正负交替型粗差。
120 图1 带粗差的井斜数据
Fig. 1 inclination date with outliers
3.2 仿真流程及环境
由于抗差估计的选权迭代法中，临界值的选取在平差过程中显得相当重要，它的合理选
125 取直接关系到平差的结果和平差是否能顺利进行[9]，故每种方法都选定了几组各具特点的临
界值。
将基于稳健估计的测斜数据抗差算法在VC 环境下的MFC 平台下编制，并在Windows
XP Service Pack 3 操作系统中进行仿真，过程如下，仿真流程图见图2。
1) 确定数据迭代步长及更新步长。由于本系统数据更新间隔为100ms，在井下运行速
130 度为1～2m/s，根据相关石油标准[10]同时考虑到数据处理的稳健性，将数据处理容
量定为25 个，每次更新10 个数据，即用本次得到校正的15 个数据加上10 个待处
理数据作为下次处理的数据。同时起始端先处理50 个数以保证初始信息的准确性。
2) 选定具体的权函数，根据权函数不同每次计算不同的量，如均值、标准差等。每次
迭代更新相应数据的权重。
135 3) 设定迭代收敛规则。根据井斜角的实际工程情况，当前后两次迭代的加权平均值之
差的绝对值小于0.3°时，即停止迭代。将校正数据输出至最终文件。循环(1)至(3)，
将所有待校正数据进行抗差处理。
 140 图2 算法流程图
Fig. 2 algorithm flow graph
3.3 仿真结果
对于Huber 法和丹麦法， c 分别取2. 5σˆ 、 2. 0σˆ 、1.5σˆ 、σˆ 、0. 8σˆ 、0. 4σˆ ， 试
145 验数据结果如表1 和表2 所示。对于IGG 法，取了5 组不同的k0 和k1，比较抗差效果，
如表3 所示。对于Hampel 法，分别取5 组数分析a，b，c 的取值对估计结果的影响，结果
如表4 所示。此外，表5 总结了将四种方法抗差处理结果达到100%准确率时临界取值及运
行效率的对比情况。图3 显示了经过抗差处理后的数据情况，证明了上述四种算法抗差的有
效性。
 图 3 抗差结果
Fig. 3 anti-outliers result
学术论文网Tag：代写论文 代写代发论文 代发论文 职称论文发表