基于泛克里格方法的位场扩边处理

李 盼， 戴前伟，2， 吕宏安

(1.中南大学 地球科学与信息物理学院，长沙 410083；2.中南大学 有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083；3.湖南省永吉高速公路建设开发有限公司，吉首 416000)

0 引言

一些常见的位场数据处理与转换可以在空间域或者频率域内进行[1-2](位场导数[3]、位场延拓[4]等)。

由于在频率域的处理相对于空间域的处理更为简洁[5]，所以目前的位场数据处理与转换更多的是在频率域内进行的。

在频率域内对位场数据进行快速傅立叶变换时要求数据为规则网格数据，为了减弱边界效应以及Gibbs效应[6]，要求网格数据在x、y上的点数均为2的幂次方，并且扩边后的边界值均等于同一数值[7]。所以当原始的位场网格数据不满足上述要求时，就必须采用适当的扩边算法对原始数据进行扩边处理以满足频率域快速傅立叶变换的要求。目前比较常用的扩边算法是余弦扩边，该方法利用余弦公式对数据区域进行扩展，扩展后的边界值均为零，从而满足了频率域快速傅立叶变换的要求。

虽然余弦扩边方法有着简单易实现的特点，但是作为一种优秀的扩边算法，应该使得扩展后的数据区域能够与原始数据区域平滑、连续地衔接，并且扩展后的数据能尽可能地接近真实位场值，很明显余弦扩边并不能保证扩边区域与原始区域的平滑连接，更不能保证扩展区域的数据接近真实的位场值，因此,寻找一种更加合理的扩边方法是有必要的。扩边算法的实质是利用空间插值算法估计出未知区域的位场值，前人具体研究过各种插值方法在地球物理数据网格化中的应用[8-10]，其中克里格方法是公认的效果最好的空间插值方法之一，作为普通克里格的改进，泛克里格方法考虑到了非平稳随机变量的漂移与波动，充分利用了数据点的空间相关性，能够尽可能准确地对未知点进行估计。其他常见的空间插值方法有反距离加权法、最小曲率法、径向基函数法、自然邻点插值法、线性插值三角网法等方法。其中反距离加权方法容易出现“牛眼”现象[11]，不宜用来做扩边处理；而最小曲率法虽然能很好地保证数据的圆滑，但是在数据稀疏区域插值结果可能会出现严重的“震荡”现象[12]，严重影响扩边效果，需要对其稳定性做进一步改进。笔者设计了理论模型对泛克里格扩边方法的效果进行了验证，同时为了说明该方法的效果，还将其与余弦扩边方法以及其他几种效果较好的径向基函数法、自然邻点插值法、线性插值三角网法的扩边效果进行了比较，对原始数据进行扩边，然后利用涉及到频率域处理的ISVD算法[13]换算位场垂向一阶导数，通过比较换算结果与理论值的误差来说明各个扩边方法的优劣性。

1 泛克里格插值方法原理[14]

1.1 非平稳随机变量的漂移与波动

在普通克里金方法中，区域变量Z(x)的数学期望是假设不变的，而在泛克里格方法中随机变量Z(x)在研究区域内的数学期望E[Z(x)]是变化的，即：

E[Z(x)]=m(x)

(1)

非平稳随机变量包括两个部分：①漂移；②波动。其中漂移是指随机变量在点x处的数学期望，即式(1)中的m(x)，而波动R(x)是指点x处随机变量与漂移的差。漂移与波动的关系为式(2)。

Z(x)=m(x)+R(x)

(2)

从实际意义上来说，漂移可以理解为变量在研究区域的某种明确的变化趋势，而波动则是随机变量在m(x)附近的随机摆动，波动的数学期望为零，即：

E[R(x)]=E[Z(x)]-E[m(x)]=0

(3)

漂移一般用一次或二次多项式表示，即：

(4)

式中：fl(x)为已知函数；μl为未知系数。

在曲面插值中，漂移通常采用式(3)表示。

m(x,y)=μ0+u1x+u2y+u3xy+u4x2+u5y2

(5)

1.2 泛克里格方程组

区域变量Z(x)在待插值点x0处的真实值Z(x0)是未知的，因此需要用x0影响范围内的测量值的加权平均来估计待插值点的值：

(6)

而待插值点x0处的真实值Z(x0)与估计值Z*(x0)之差应该满足无偏条件，即：

E[Z*(x0)-Z(x0)]≡0

(7)

结合式(6)、式(7)及式(1)、式(4)可得：

(8)

如果要使式(8)对任何一组系数都成立，则必须满足：

(9)

式(9)即为用k+1个方程表示的无偏条件。

文献[9]推导了用变差函数表示的估计方差表达式：

(10)

采用拉格朗日乘数法，将无偏条件式(9)代入式(10)得目标函数式(11)。

(11)

将式(11)对λi及μl求一阶偏导，令其为“0”，

便可以得到泛克里格方程组(12)。

图1 理论模型的平面图及其重力场值的灰度图Fig.1 The model plane position and its gravity anomaly gray image

表1 理论模型的参数

(12)

2 理论模型实验

笔者设计了图1所示的理论模型对本文提出的方法效果进行验证，为了更加明显地突出边界处理的效果，特意在数据区域边界放置了四个长方形异常体2、3、4、5，模型参数如表1所示。

用文献[15]介绍的式(13)、式(14)计算了模型的重力场及其垂向一阶导数理论值，网格密度为101×101。

(13)

(14)

分别用余弦法、泛克里格法、径向基函数法、自然邻点插值法、线性插值三角网法对重力场数据进行扩边，扩边后的边界值统一设定为零，然后利用离散的网格数据基于ISVD算法[13]换算出其垂向一阶导数。图2、图3是基于各种扩边方法通过频率域处理换算垂向一阶导数与理论值的对比及误差绝对值对比。表2统计了基于各种扩边方法换算出的垂向一阶导数与理论值的误差。

图2 基于各种扩边方法通过频率域处理换算垂向一阶导数与理论值的对比Fig.2 The compare between transformed vertical first derivative and the theoretical value based on different edge extending methods

图3 基于各种扩边方法换算出的垂向一阶导数与理论值的误差绝对值Fig.3 The absolute value of the error of transformed vertical first derivative and the theoretical value based on different edge extending methods

图2是基于各种扩边方法对原始数据进行了扩边之后，利用ISVD算法换算出的位场垂向一阶导数，并基于式(14)计算出的理论值作为对比。通过对比可以看出来基于离散值换算出的结果与理论值是有一定误差的，误差主要出现于曲线的“峰顶”和“谷底”处，基于各种扩边方法换算的结果都非常接近，主要区别在于数据边界位置，即图2两个红色椭圆所圈出的区域，在其他地方，这五种扩边方法换算结果曲线都比较吻合，而在边界出则出现了比较明显的分离，其中青色虚线所代表的余弦扩边法出现了比较明显的偏差，其余四种方法比较接近。图3给出了各种方法的换算结果与理论值的误差绝对值的曲线图，由图3可以看出，在靠中间的区域，所有方法的误差曲线都非常接近，只存在细微的差距，主要区别在边界处，余弦法的表现是比较差的，出现了明显的误差陡增的现象，这主要是余弦法在扩边时，对于原始数据与扩边数据的衔接处处理得不够光滑，使得连接处出现“沟壑”现象，使得数据严重失真，这必然造成换算结果的误差增大。而其他四种方法相对来说表现比较接近，在对应图2椭圆所圈出的边界处只出现了较小幅度的误差增加，整体来说，没有出现明显的误差陡增，比较平稳。仔细对比的话，泛克里格方法对应的误差曲线处于最下方，也就是说误差最小。表2计算了各种扩边方法换算结果与理论值的最大误差与均方根误差，由表2可以看出，泛克里格方法无论是最大误差还是均方根误差都是最小的，而余弦法则是表现最差的。

表2 基于各种扩边方法换算出的垂向一阶导数与理论值的误差Tab.2 The error of transformed vertical first derivative and the theoretical value based on different edge extending methods

3 实际应用

我们选取了Gooper于2008年在SEG网站上提供的南非Witwatersrand Basin某区域的开源重力数据，利用泛克里格扩边方法对原始数据做了扩边处理，然后用ISVD算法换算了其一、二、三阶垂向导数，如图3所示。由图3可以看出，一阶垂向导数的灰度图相对原始重力异常而言能够观察到非常明显的地质体主体框架，右上角的环状地质体边界清晰可辨，其左侧和下方的地质结构体相对原始重力异常图有明显的突显；而二阶和三阶垂向导数则突出了更多的地质细节，地质体与周围围岩区分明显，能够观测到明显的褶皱和断层的分布和走向。

图4 实测重力异常数据及其一、二、三阶垂向导数Fig.4 The measured gravity anomaly data and its first, second and third order vertical derivative

4 总结

鉴于传统的余弦位场扩边方法不能保证扩边后位场数据的连续且平滑的要求，笔者旨在从目前常用的地球物理数据网格化方法中寻找一种合适的扩边方法，通过参考多篇相关文献，筛除了效果不佳的方法，从中挑出了效果较好的泛克里格法、径向基函数法、自然邻点插值法、线性插值三角网法这几种方法与余弦扩边法对位场数据做了扩边处理，通过换算其垂向一阶导数对比了各自的精度。结果表明，泛克里格扩边方法的精度明显高于余弦扩边法，而与其他扩边方法相比，其精度也是最高的。通过实测数据的应用结果表明，基于泛克里格扩边的位场数据的高阶导数边值处理是比较成功的。