基于泰勒级数展开的位场视延拓成像方法

王彦国，　张　瑾，　聂逢君

(1.东华理工大学核工程与地球物理学院，江西 南昌　330013；2.东华理工大学地球科学学院，江西 南昌　330013)



基于泰勒级数展开的位场视延拓成像方法

王彦国1，张瑾1，聂逢君2

(1.东华理工大学核工程与地球物理学院，江西 南昌330013；2.东华理工大学地球科学学院，江西 南昌330013)

摘要：位场快速反演算法-位场视延拓成像法是利用向下延拓因子的泰勒级数展开代替理论向下延拓来进行近似延拓，然后再将近似延拓场向上延拓到某一深度上，最后将该深度上延拓场再次利用泰勒级数展开进行近似向下延拓实现的。该方法可以在无任何地质约束情况下快速有效地完成场源成像工作。滤波特性表明该方法属于中低通滤波，对高频干扰具有较强的压制能力。模型试验表明泰勒级数展开项数越少，成像异常的分辨率越高；且位场梯度数据的成像准确度更高；二维及三维实例应用表明视延拓成像方法可有效地反映出地质体的空间赋存状态，具有良好的纵向和横向分辨能力。

关键词：位场；泰勒级数展开；视延拓成像；快速反演

王彦国，张瑾，聂逢君.2016. 基于泰勒级数展开的位场视延拓成像方法[J].东华理工大学学报：自然科学版，39(1)：72-79.

Wang Yan-guo, Zhang Jin, Nie Feng-jun.2016. Apparent continuation imaging for potential-field based on Taylor series expansion [J].Journal of East China University of Technology (Natural Science), 39(1):72-79.

利用重磁数据反演地下地质体的空间赋存状态是位场进行地质解释的重要手段。常规反演算法是基于反演理论，在最小二乘法下使目标函数达到极小的线性或非线性反演(郭良辉等，2009；张维宸等，2008)。但这类方法需要引入先验信息约束，计算维数和计算量较大，制约了其在三维实际资料中的应用。近年来，Portniaguine等(2002)、Li等(2003)、姚长利等(2003，2007)提出了一系列改进措施，在一定程度上改善了反演的计算效率。重力归一化总梯度(肖一鸣等，1984；曾华霖等，1999；肖鹏飞等，2006)不需附加任何地质条件下能半定量地确定场源位置，已被广泛应用于寻找油气构造中，但此法受谐波数、向下延拓的影响，使场源深度和水平位置的估计受到了一定限制，且该方法不能同时反映不同深度上的场源位置。Euler反褶积法(Thompson, 1982; Reid, et al, 1990; Hansen et al., 2002)也可以用于场源几何参数反演，该方法是利用位场一阶导数和计算点点位坐标实现的，导数计算结果或计算点位稍有偏差,便会在反演中产生较大的偏离，从而影响反演精度。重磁相关成像法是由郭良辉等(2009，2010)引入的，该方法同样无需先验信息约束，计算方便稳定，对地下场源具有良好的成像效果，但需事先确定地下场源的几何形状；马国庆等(2013)针对相关成像的不足进行了改进，使得精确反演地质体位置的同时还可以检测出场源的几何类型。

位场向下延拓是一个典型的不适定问题，目前解决这一问题行之有效的措施是迭代方法。王彦国等(2011)提出了位场向下延拓的泰勒级数迭代法，该迭代法相对于常用的积分迭代法(徐世浙，2006)，具有收敛速度快和保幅能力强的优点，但方法更易受干扰影响。为此，曾小牛等(2013)对泰勒级数迭代法进行了改进，提出了正则-积分迭代法，有效地解决了高频干扰带来的不稳定性问题。本文在前人工作基础之上，提出了基于泰勒级数展开的位场视延拓成像方法。该算法不仅适合重力及其各阶导数，而且还可以对磁异常(直接)解析信号进行处理。文中利用二维、三维及含噪模型重力数据和二、三维实际资料进行了试验和验证。

1视延拓成像的基本原理及其滤波特性

1.1基本原理

在波数域中，位场向上延拓算子φ和向下延拓算子ψ分别为：

φ=e-ωh，ψ=eωh

(1)

图1　视延拓滤波特性曲线Fig.1　Filtering characteristic of apparent continuationa.N的影响；b.k的影响

其中，ω是波数，h是延拓深度。王彦国等(2011)给出了下延算子ψ在h=0处Taylor展开的前N项和φ0→h：

(2)

本文同样基于迭代算法的回返思想，提出了一种通过极值反映场源质心位置的方法——视延拓成像法，其原理如下：

采用φ0→h进行近似向下延拓，然后采用φ将向下延拓的近似波谱向上延拓到某一深度：

φ0→kh=e-ω(1-k)h

(3)

其中，k为延拓返回系数，则深度为kh的位场波谱近似延拓算子可以表示为：

(4)

由Φkh作用于原始波谱得到的是深度为kh面上的似波谱，为进一步得到h深度的似波谱，需对kh深度上的波谱进行向下延拓，在此仍采有上述的Taylor级数前N项和代替FFT下延算子，即对下延算子在深度kh处进行Taylor级数展开：

(5)

由此得到h深度的近似延拓算子Φh：

(6)

及h深度的近似波谱Uh：

Uh=U0·Φh

(7)

其中，U0是已知平面的原始波谱。对Uh进行反变换即可得到h深度的近似异常值uh。

由于Φh是近似延拓算子的组合，故称Φh为视延拓算子，将得到h深度上的似向下延拓场称为视延拓异常。

1.2滤波特性

由公式(6)可以看出，视延拓算子Φh与展开项数N和延拓系数k有关。图1给出了Φh与N、k的关系，可以看出：(1)Φh是中低通滤波，放大中低频信号的同时压制了高频干扰，故滤波结果具有一定的稳定性；(2)当k固定时，Φh的振幅随着N的增加而增大，而且主峰逐渐向高频靠近；(3)当N固定时，Φh的振幅随着k的增加而增大同样向高频靠近。由此可知，选择较小的展开项数N和较大的回返延拓系数k的组合与选择较大展开项数和较小的回返延拓系数相当。

通过大量模型试算发现，视延拓成像中的回返系数k与地质体埋深关系不大，仅与Taylor级数展开项数N和求导阶次n有关，且展开项数N越大，回返系数k越小；求导阶次越高，回返系数k越大，具体对应关系见表1。需要指出的是，磁异常与重力异常的一阶导数相当，因此对磁异常进行视延拓成像则应选取表1中的阶次n=1。

图2　不同展开项数N的视延拓成像Fig.2　Apparent continuation imaging of different number of expansion terms(a) 重力异常；(b)N=1；(c) N=2; (d) N=3

项数N123二维三维二维三维二维三维00.360.470.170.250.090.1510.470.570.250.320.150.2020.570.630.320.380.200.25

2模型试算

2.1二维单体模型

为了说明位场视延拓成像的有效性，这里设计了一个半径为0.5 km，埋深为1 km，剩余密度为1.0 g/cm3的无限长水平圆柱体构成的单体模型，模型产生的重力异常见图2a。利用Taylor级数展开项分别为1、2、3的视延拓滤波(公式7)对重力异常(图2a)进行相应处理，计算结果分别见图2b-d。从图中可以看出，不同展开项数下视延拓场的极值所对应的位置恰为模型体的质心所在，等值线在纵向上表现为椭圆形，在横向上关于模型体呈现左右对称分布；但从等值线的疏密程度和数值变化区间可以看出，N=1时的视延拓成像结果应具有更高的分辨率，为此下文中的视延拓成像中均采用N=1进行。

2.2二维含噪叠加模型

为了说明本文方法处理复杂问题的能力，采用了三个水平圆柱体组成的叠加模型进行试算。从左到右三个模型体参数分别是：半径0.5 km，埋深1 km；半径1 km，埋深2 km；半径2 km，埋深3 km，其模型体剩余密度均为1.0 g/cm3。为了解本文方法对噪声的敏感度，对叠加模型体的重力异常添加了1%的随机噪声，含噪重力异常如图3a所示。含噪重力异常的视延拓成像结果见图3b，可以看出，视延拓成像结果稳定性较强，极大值对地质体质心位置有一定指示，但受异常叠加影响，视延拓场极大值点与地质体1、2的质心位置存在较大偏差。为进一步提高对地质体显示的准确性，对含噪重力异常的垂向一阶导数进行了视延拓成像处理，结果见图3c，易看出，随机干扰只是分布在浅层近地表，而并不影响实际地质体所反映出的异常形态，且相对于图3b，垂向导数视延拓场可以更加清晰、准确地反映出了地质体1、2的质心位置。

图3　含噪模型重力异常及其垂向一阶导数的视延拓成像Fig.3　Apparent continuation imaging of gravity anomaly including random noise and its first vertical derivative (a) 含1%噪声的重力异常；(b)重力异常的视延拓成像；(c)垂向一阶导数的视延拓成像

2.3三维模型

为说明视延拓成像可以处理三维数据，设计了两个参数不同的球体组成的模型，浅部球体质心空间位置为(2.5 km，2.5 km，1 km)，半径为0.5 km；深部地质体质心在(7.5 km，7.5 km，3 km)处，半径为1 km，两个模型体的剩余密度均为1.0 g/cm3，模型体产生的重力异常如图4所示。重力异常视延拓成像在不同深度上的切片见图5。可以看出，成像结果在空间中存在的两个极大值点同样与模型体的质心位置相对应，这说明本文方法同样可以处理三维位场数据。

图4　三维模型体重力异常Fig.4　Three dimensional forward gravity anomaly

图5　三维重力异常不同深度切片的视延拓成像Fig.5　Apparent continuation imaging of different depths based on the data of Fig. 4

3实际资料应用

3.1二维重力剖面

这里选取了二连盆地某地区的一条重力剖面进行处理与解释，该剖面所在地区地表出露的均为第三纪、第四纪沉积地层，基底钻探验证为高密度高电阻的花岗岩。图6a是该剖面的实测重力异常以及采用迭代滤波法(王彦国，2013)分离获得的区域重力异常和剩余重力异常。可以看出，分离出的剩余重力异常正负相间出现，正异常与零值线之间覆盖面积与负异常的相当，这说明了场分离的正确性。利用视延拓成像对剩余重力异常的处理结果见图6b。图中看出，该剖面的地下构造在重力上反映为“三凸夹两凹”的特征，这与电法反演结果(图6c)完全一致，同时视延拓成像结果在局部细节上的展现也与电法结果具有一定的相似性(如15 km和35 km处的异常特征)，这充分说明了本文方法的实用性。

图6　二维实测重力异常视延拓成像Fig.6　Apparent continuation imaging of two-dimensional real gravity data(a) 实测重力异常及其分离；(b) 剩余重力异常视延拓成像；(c) 同剖面上CSAMT视电阻率反演断面图

图7　内蒙古某矿区实测磁异常Fig.7　The measured magnetic anomaly of a certain mining area in inner Mongolia

3.2三维磁力数据

为进一步验证算法实用性，对内蒙古某矿区实测磁力数据进行相应处理与解释，其原始磁力异常见图7。为减小剩磁和磁化角度的影响，首先对磁异常进行了直接解析信号处理，然后对解析信号进行视延拓成像，不同深度上的视延拓成像切片见图8。可以看出，视延拓成像结果展示出大部分矿体埋藏深度在200～300 m之间，其西南角的磁性体为马蹄形，且磁化强度较大，中西部的磁性体主要呈现为条带状，并以NE走向为主。

图8　实测磁异常直接解析信号不同深度上的视延拓成像Fig.8　Apparent continuation imaging at the difference depths based on direct analytic signal of real gravity data

4结论

本文在泰勒级数迭代法的研究基础上，提出了一种可以反演地质体质心位置的视延拓成像方法。该算法无需地质先验信息约束，方法简单，易于实现，计算稳定快速，是一种在深部矿产资源勘查和地质调查中有着良好应用前景的方法。

参考文献

郭良辉，孟小红，石磊，等. 2009. 重力和重力梯度数据三维相关成像[J]. 地球物理学报，52(4):1098-1106.

郭良辉，孟小红，石磊. 2010. 磁异常△T三维相关成像[J]. 地球物理学报，53(2):435-441.

马国庆，杜晓娟，李丽丽. 2013. 改进的位场相关成像方法[J]. 地球科学-中国地质大学学报，38(5):1121-1127.

王彦国，张凤旭，王祝文，等. 2011. 位场向下延拓的泰勒级数迭代法[J]. 石油地球物理勘探， 46(4):657-662.

王彦国. 2013. 位场数据处理的高精度方法研究及应用[D]. 吉林大学, 31-35.

肖鹏飞，李明，徐世浙，等. 2006. 重力归一化总梯度的稳定算法[J]. 石油地球物理勘探，41(5):596-560.

肖一鸣，张林详. 1984. 重力归一化总梯度法在寻找油气中的应用[J]. 石油地球物理勘探，(3):247-254.

徐世浙. 2006. 位场延拓的积分-迭代法[J]. 地球物理学报，49(4):1176-1182.

姚长利，郝天珧，管志宁，等. 2003. 重磁遗传算法三维反演中高速计算及有效存储方法技术[J]. 地球物理学报，46(2):252-258.

姚长利，郑元满，张聿文. 2007. 重磁异常三维物性反演随机子域法方法技术[J]. 地球物理学报，50(5):1576-1583.

曾华霖，李小孟，姚长利，等. 1999. 改进的重力归一化总梯度法及其在胜利油区油气藏探测中的应用效果[J]. 石油勘探与开发，26(6):1-6.

曾小牛，李夕海，牛超，等. 2013. 位场向下延拓的波数域正则-积分迭代法[J]. 石油地球物理勘探，48(4):643-650.

张维宸，刘建芬，谢连文，等. 2008. 利用航磁数据推断隐伏(半隐伏)岩体[J]. 东华理工大学学报：自然科学版，31(4):349-356.

Hansen R O, Laura S. 2002. Multiple-source Euler deconvolution[J]. Geophysics, 67(2):525-535.

Li Y G, Oldenburg D W. 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method[J]. Geophys. J. Int., 152:251-265.

Portniaguine O, Zhdanov M S. 2002. 3-D magnetic inversion with data compression and image focusing[J]. Geophysics, 67:1532-1541.

Reid A B, Allsop J M, Granser H, et al. 1990. Magnetic interpretation in three dimensions using Euler deconvolution[J]. Geophysics, 55(1):80-91.

Thompson D T.1982. EULDPH:A new technique for making computerassisted depth estimates from magnetic data[J]. Geophysics, 47(1):31-37.

Apparent Continuation Imaging for Potential-field Based on Taylor Series Expansion

WANG Yan-guo1,ZHANG Jin1,NIE Feng-jun2

(1. College of Nuclear Engineering and Geophysics, East China University of Technology, Nanchang,JX 330013, China; 2 College of Earth Sciences, East China University of Technology, Nanchang, JX 330013, China)

Abstract:This paper presents a fast inversion method called apparent continuation imaging method in potential-field. This method can make imaging of field source without any priori geological information, and can be achieved by three steps. The first step is using Taylor series expansion of downward continuation operator to obtain approximate continuation field; The second step is taking the approximate continuation field above to be upward continuation of a certain depth; The last step is once again using Taylor series expansion of downward continuation operator to be downward continuation from the certain depth to the given depth. Filtering characteristic shows that the method presented by this paper is low-middle pass filtering and can reduce the high-frequency interference. The single model test indicates that field imaging resolution is higher when the number of Taylor series terms is smaller, and combined model including random noise shows that the accuracy of inverted result using gradient data is higher. Two and three real applications show that apparent continuation imaging method can reflect the spatial extension states of geological bodies with high resolution.

Key Words:potential field; Taylor series expansion; apparent continuation imaging; fast inversion

中图分类号：P631

文献标识码：A

文章编号：1674-3504(2016)01-0072-08

doi:10.3969/j.issn.1674-3504.2016.01.012

作者简介：王彦国(1985—)，男，博士，讲师，主要从事重磁数据处理与解释方法的研究. E-mail:wangyg8503@126.com

基金项目：国家自然科学基金(41504098)；国防科工委科研项目(科工技〔2013〕969号)；东华理工大学博士科研启动基金(DHBK2013213)。

收稿日期：2015-09-01