基于DCT的广义指数阈值衰减凸集投影算法在位场数据补空中的应用

图1 POCS算法原理^[10]

Fig.1 The mechanism of POCS method^[10]

基于DCT的数据补空凸集投影理论公式:

(2)

D_{k} = D_{obs} + M C^{- 1} T_{k} C D_{k - 1}, k = 1,2, 3, \dots, K,

式中:K是总迭代次数;D_obs(x,y)是含缺失数据的原始观测数据;D_k代表第k次迭代补空的数据;C和C^-1代表二维离散余弦正、反变换;M是采样矩阵,其与输入数据维度相同,如果某点值为0,代表该点有数据,无需插值,反之为1,则代表该点需要进行插值计算;T_k是门限值矩阵,可表示为:

(3)

T_{k} (i, j) = \{\begin{array}{l} 1, & | S_{k} (i, j) | \geq p_{k} \\ 0, & | S_{k} (i, j) | < p_{k} \end{array}

式中:p_k是第k次迭代阈值,S_k是第k次的DCT变换谱。

Abma^[13]给出的线性阈值衰减形式为:

(4)

p_{k} = p_{k}^{\max} - \frac{(k - 1) (p_{k}^{ma x} - p_{k}^{\min})}{K - 1},

式中: $p_{k}^{\max}$ =max(abs(S_k)), $p_{k}^{\min}$ =min(abs(S_k)),k=1,2,…,K。

本文给出的广义指数阈值衰减形式为:

(5)

p_{k} = p_{k}^{\max} * \exp [- {(\frac{k - 1}{K - 1})}^{Para} \ln (p_{k}^{\max} / p_{k}^{\min})]

式中:Para为控制衰减速度的参数,当Para=1时为Gao等^[10]提出的阈值衰减情况,当Para=0.5时即为张华等^[11]提出的阈值衰减形式。

凸集投影算法进行数据补空的基本流程为:对原始输入数据采用DCT处理,给定一个阈值,使得变换谱中振幅大于等于这个阈值时保留,振幅低于这个阈值则充零代替,然后对该数据进行反离散余弦变换乘以采样矩阵,最后将其与原始待补空数据叠加,而在下一次迭代计算时,减小阈值重复上述操作,直到达到结束迭代计算的条件。本算法需要实现对不规则数据网格化,以满足DCT变换的需要,算法的关键在于选择合适的稀疏变换和阈值衰减方式,其工作流程图见图2。

图2

图2 基于DCT的POCS数据补空流程

Fig.2 The flow-process diagram of POCS method to interpolate based on DCT

2 模型试验

为了验证指数阈值衰减凸集投影算法在位场中的应用效果及优越性,建立了一个由参数不同的4个模型体产生的叠加重力异常进行试算,使用的模型体参数如表1所示,模拟测区面积为500 m×500 m,点距与线距均为10 m。图3是完整的叠加重力异常,其中黄色虚线所圈部分为地质体边缘在地面的正投影,粉红色锯齿圈出区域为待补空部位。

表1 模型参数

Table 1 Model parameters

模型编号	模型类型	x方向范围/m	y方向范围/m	z方向范围/m	密度/(g·cm^-3)
Model 1	圆柱体	-200~200	-190~-110	10~90	2
Model 2	棱柱体	-50~50	-50~50	20~400	-2
Model 3	棱柱体	-200~0	100~200	30~1000	2
Model 4	球体	100~200	100~200	10~110	2

图3

图3 组合模型产生的重力异常及缺失部位(粉色部分)

Fig.3 Forward gravity anomaly with complete data and the target zone (pink part)

图4展示了Surfer软件中4种常用的插值算法所获得的数据补空结果及其与图3的残差情况。可以看出,所有的常规插值算法都在模型体1处存在等值线的间断,尤其反距离加权平均法的最为明显,最小曲率法的等值线畸变最不明显,但最小曲率法的残差值最大,其次是反距离加权法的。另外,常规插值方法的残差均在原始数据的正异常处为负值,而负异常处为正值,这表明了补空的数据没有达到真实值的幅度。表2给出了4种方法的插值均方根误差(RMS),可以看出,反距离加权法的数据补空精度最差,其次是最小曲率法,二者误差是克里金法和径向基函数法的2倍左右。

图4

图4 常规插值方法的数据补空及与原数据的残差

a、b—克里金法插值及残差;c、d—径向基函数法插值及残差;e、f—反距离加权法插值及残差;g、h—最小曲率法及残差

Fig.4 Results of interpolation using conventional methods and the differences with the original data

a、b—interpolated result of Kriging and its residuo;c、d—interpolated result of radial basis function and its residuo;e、f—interpolated result of inverse distance to a power and its residuo;g、h—interpolated result of minimum curvature and its residuo

表2 不同常规插值方法的数据补空误差

Table 2 RMS errors of interpolating data using different conventional interpolation methods

处理方法	均方根误差/mGal
克里金法	0.123
径向基函数法	0.113
反距离加权法	0.285
最小曲率法	0.217

为了了解线性阈值及指数阈值衰减凸集投影与总迭代次数关系,图5绘制了不同模式下的凸集投影方法数据补空误差与最大迭代次数K关系曲线。可以看出,无论是线性衰减,还是不同参数的指数衰减,数据补空误差均是随着总迭代次数K增加而逐渐减小,不过指数衰减的速率要明显大于线性衰减,其中Para=0.5时衰减速率最快。线性阈值衰减方式POCS算法基本收敛在0.070 mGal,指数阈值衰减基本收敛在0.030 mGal,不过指数阈值衰减的误差随着迭代次数变化呈现出了一定的波动性。对比图5与表2可以看出,线性阈值衰减、指数阈值衰减凸集投影算法的数据补空精度分别是常规方法效果最好的径向基函数法的1.6倍和4.5倍,这表明了采用基于指数阈值衰减的凸集投影算法能够获得更好的数据补空效果。

图5

图5 线性及指数阈值衰减凸集投影算法的理论数据补空误差与总迭代次数K关系曲线

Fig.5 The relationship between the errors of interpolating data by using POCS with different threshold models and the total number of iterations

图6则是基于线性与指数阈值衰减的POCS算法迭代800次数据补空后的重力异常及其与图3的残差,可以看出,无论线性阈值,还是不同参数指数阈值POCS法处理后的重力异常等值线均较圆滑,并不存在明显的等值线畸变情况。从残差图中可以看出,POCS算法的残差在原始数据表现为正异常区域内为正值,而负异常区为负值,这恰于常规插值方法数据补空的残差是相反的,也就是说,基于POCS算法的数据补空的幅值略大于实际数据的。

图6

图6 基于线性、指数阈值衰减的凸集投影算法数据补空结果及与原数据的残差(K=800)

a、b—线性阈值插值及残差; c、d—指数阈值(Para=0.5)及残差;e、f—指数阈值(Para=1)及残差;g、h—指数阈值(Para=2)及残差

Fig.6 Results of interpolation using POCS method with linear and exponential threshold and the difference with the original data (K=800)

a、b—interpolated result using linear threshold and its residuo;c、d—interpolated result using exponential threshold (Para=0.5) and its residuo;e、f—interpolated result using exponential threshold (Para=1) and its residuo;g、h—interpolated result using exponential threshold (Para=2) and its residuo

为了进一步证实本文算法的有效性,对理论重力异常(图3)添加了30%随机噪声(图7),然后再对含噪重力异常进行了数据缺失处理(图7粉色部位)。图8是4种常用插值方法获得的含噪重力数据补空结果及其与无噪声重力异常的残差,表3给出了4种插值方法的误差。对比表3与表2可以看出,相对无噪声时,含噪数据的常规插值方法补空误差有所增加,不过规律基本保持,还是反距离加权法的误差最大,最小曲率法次之。从图8可以看出,4种常规插值方法获得的补空数据仍在在模型体1位置存在等值线的畸变情况,另外还可以看出插值区域的噪声水平与周边数据噪声水平明显不一致,补空区的等值线更加圆滑些。

图7

图7 图3添加30%随机噪声后的重力异常及缺失部位

Fig.7 Gravity anomaly and missing site after adding 30% random noise of Fig. 3

图8

图8 常规插值方法对含噪缺失数据体的补空结果及与原数据的残差

a、b—克里金法插值及残差;c、d—径向基函数法插值及残差;e、f—反距离加权法插值及残差;g、h—最小曲率法及残差

Fig.8 Results of interpolating noise-corrupted data using conventional interpolation methods and the difference with the original noise-free data

表3 不同常规插值方法对含噪重力异常的数据补空误差

Table 3 RMS errors of interpolating noise-corrupted data using different conventional interpolation methods

插值方法	均方根误差/mGal
克里金法	0.191
径向基函数法	0.185
反距离加权法	0.315
最小曲率法	0.280

图9是线性及指数阈值衰减POCS算法的补空的关系曲线。可以看出,线性阈值衰减POCS法的补空数据仍然是随着总迭代次数的增加而减小,基本收敛在0.137 mGal,但不同para时的指数阈值衰减POCS法数据插值误差则是随着迭代增加,先减小后增加,最终基本收敛在0.272 mGal处,推测是误差增加是由于噪声的反复重建导致的。图10可以看出,线性阈值衰减的POCS法(图10a)也在模型体1处产生了明显的等值线间断情况,补空区内的异常同样较为光滑,与残差图(图10b)仍是南、北两侧为大面积的正差值,而中部则整体表现为大面积数据与无噪声理论值之间均方误差与总迭代次数K的负差值,这表明含噪数据进行补空时,线性阈值衰减的POCS法数据填充的误差主要是数据未真实恢复引起的。指数阈值衰减的POCS法(图10c、e、g)补空的数据明显也存在着局部的波动,即含有一定的噪声,残差图(图10d、f、h)则是正负异常间隔出现,即认为主要误差是由噪声干扰引起的。

图9

图9 含噪声时凸集投影算法的数据补空误差与总迭代次数K关系曲线

Fig.9 The relationship between the errors of interpolation data using POCS with different threshold models for the noise-corrupted data and the total number of iterations

图10

图10 图7采用线性与指数阈值衰减POCS算法的数据补空结果及其残差结果

a、b—线性阈值插值及残差; c、d—指数阈值(Para=0.5)及残差;e、f—指数阈值(Para=1)及残差;g、h—指数阈值(Para=2)及残差

Fig.10 Results of interpolation using POCS method with linear and exponential threshold model and the differences with the noise-free data

3 应用实例

为了检验本文算法的实际资料处理能力,对黑龙江嫩北农场布格重力异常进行挖空与数据补空处理。图11是网格密度500 m×250 m的布格重力异常,布格重力异常总误差为0.018 mGal。由于该区中部NWW流向的科洛河在汛期时河床较宽,最宽处可达1 km,汛期时重力数据将无法采集,故这里对科洛河的大致流域区域进行了数据挖空处理(图11粉色部位)。

图11

图11 黑龙江嫩北农场布格重力异常及缺失部位(粉色部分)

Fig.11 Bouguer gravity anomaly of Nenbei Farm in Heilongjiang Province and the target zone (pink part)

图12是理论模型测试所提及的4种常用插值方法处理的结果及其残差分布,可以看到,这4种插值方法数据补空结果基本一致,很难看出彼此差异,且等值线畸变现象也难以直接显示,其主要原因是挖空区域宽度仅为2~3个点距。从残差图上可以看出,补空数据与实际数据的残差仍是大面积的正、负值,即误差具有一定的范围,这显然不利于后期数据处理与解释。表4是的4种常规插值方法的数据补空误差,可以看出,径向基函数法的数据补空误差最小,其次是最小曲率法和克里金法,反距离加权法的误差仍是最大,不过所有常规方法的数据补空误差均远大于布格重力异常的总误差,也就是说补空数据质量较差。

图12

图12 实际数据的常规插值方法数据补空结果

a、b—克里金法插值及残差;c、d—径向基函数法插值及残差;e、f—反距离加权法插值及残差;g、h—最小曲率法及残差

Fig.12 Results of interpolation using conventional methods and the differences with the real data

表4 不同常规方法对实际资料进行数据补空的误差

Table 4 The errors of interpolation results using conventional methods for the read data

处理方法名称	均方根误差/mGal
克里金法	0.105
径向基函数法	0.091
反距离加权法	0.298
最小曲率法	0.105

图13给出了基于线性、不同参数Para指数阈值衰减POCS算法数据补空均方根误差与总迭代次数K的关系曲线,可以看出,无论是线性、还是不同参数的指数阈值衰减,POCS算法数据重建误差均随着K增加而逐渐减小,线性阈值衰减POCS算法的误差基本收敛于0.149 mGal,远高于布格重力异常总误差,即线性阈值衰减的POCS算法重建的数据质量也较低;指数阈值衰减POCS算法的误差收敛于0.008 mGal,则低于布格重力异常总误差,重建数据质量较高,重建数据完全可以替代真实数据。图14给出了Para=0.5时指数阈值衰减POCS法的处理结果以及残差分布,残差图中显示了缺失数据宽度较小的西北侧基本为0,而补空区较宽的东南侧基本是正负异常间隔出现,且数据幅值较小,这体现出了本文方法在数据补空方面的优势。

图13

图13 不同衰减方式POCS补空的误差与总迭代次数K的关系曲线

Fig.13 The relationship between the errors of interpolating data using POCS with different threshold models and the total number of iterations

图14

图14 实际数据指数阈值衰减POCS算法处理结果

a—指数阈值(Para=0.5)数据补空结果;b—与实测数据的残差 (K=6500)

Fig.14 Using exponential threshold-based POCS method to interpolate real field data

a—interpolated result (Para=0.5);b—the residuo (K=6500)

4 结论

针对缺失的重磁数据需要有效恢复的问题,本文采用了基于DCT的广义指数阈值衰减凸集投影算法进行数据补空工作。通过模型分析与实例应用,证实了基于指数阈值衰减的凸集投影算法在数据补空中存在补空精度高及补空痕迹小的优点,另外补空数据还能与实际数据噪声含量保持一致。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

孙洪泉

. 地质统计学及其应用[M]. 北京: 中国矿业大学出版社, 1990.

Sun H

. Geological statistics and its application [M]. Beijing: China University of Mining and Technology Press, 1990.

[2]

Briggs I

Machine contouring using minimum curvature

[J]. Geophysics, 1974, 39(1):39-48.

DOI:10.1190/1.1440410 URL [本文引用: 1]

[3]

Cordell

A scattered equivalent-source method for interpolation and gridding of potential-field data in three dimensions

[J]. Geophysics, 1992, 57(4):629-636.

DOI:10.1190/1.1443275 URL [本文引用: 1]

[4]

郭志宏

一种使用的等值线型数据网格化方法

[J]. 物探与化探, 2001, 25(3):203-208.

Guo Z

A practical contour type data gridding technique

[J]. Geophysical and Geochemical Exploration, 2001, 25(3):203-208.

[5]

闫浩飞, 刘国峰, 薛典军, 等.

基于凸集投影方法的重磁数据规则缺失重建

[J]. 地球物理学进展, 2016, 31(5):2192-2197.

Yan H

, Liu G

, Xue D

, et al.

Reconstruction of gravity/magnetic data with the projection-onto-convex-sets methods

[J]. Progress in Geophysics, 2016, 31(5):2192-2197.

[6]

王万银, 邱之云, 刘金兰, 等.

位场数据处理中的最小曲率扩边和补空方法研究

[J]. 地球物理学进展, 2009, 24(4):1327-1338.

Wang W

, Qiu Z

, Liu J

, et al.

The research to the extending edge and interpolation based on the minimum curvature method in potential field data processing

[J]. Progress in Geophysics, 2009, 24(4):1327-1338.

[7]

王明, 刘前坤, 李芳, 等.

基于热传导模型位场网格数据补空方法研究

[J]. 物探与化探, 2015, 39(S1):144-151.

Wang

, Liu Q

, Li

, et al.

Filling grid dummy values by heat conduction model

[J]. Geophysical and Geochemical Exploration, 2015, 39(S1):144-151.

[8]

曾小牛, 李夕海, 刘继昊, 等.

基于凸集投影的重力数据扩充下延一体化方法

[J]. 石油地球物理勘探, 2019, 54(5):1166-1173.

Zeng X

, Li X

, Liu J

, et al.

An integration of interpolation, edge padding, and downward continuation for gravity data based on the projection onto convex sets

[J]. Oil Geophysical Prospecting, 2019, 54(4):1166-1173.

[9]

曾小牛, 李夕海, 侯维君, 等.

基于凸集投影的重力同时填充扩边和去噪方法

[J]. 石油地球物理勘探, 2020, 55(1):197-205.

Zeng X

, Li X

, Hou W

, et al.

Simultaneous interpolation, edge padding and denoising method for gravity data based on the projection onto convex sets

[J]. Oil Geophysical Prospecting, 2020, 55(1):197-205.

[10]

Gao J

, Chen X

, Li J

, et al.

Irregular seismic data reconstruction based on exponential threshold model of POCS method

[J]. Applied Geophysics, 2010, 7(3):229-238.

DOI:10.1007/s11770-010-0246-5 URL [本文引用: 4]

[11]

张华, 陈小宏.

基于jitter 采样和曲波变换的三维地震数据重建

[J]. 地球物理学报, 2013, 56(5):1637-1649.

Zhang

, Chen X

Seismic data reconstruction based on jittered sampling and curvelet transform

[J]. Chinese J. Geophys., 2013, 56(5):1637-1649.