地球物理学报  2012, Vol. 55 Issue (7): 2450-2461   PDF    
解释位场全张量数据的张量局部波数法及其与常规局部波数法的比较
马国庆 , 杜晓娟 , 李丽丽     
吉林大学地球探测科学与技术学院, 长春 130026
摘要: 全张量探测技术以其信息量大、精度高、干扰小等优点在地球物理领域中得到广泛应用.本文提出采用张量局部波数法来进行位场全张量数据的解释,首先给出了张量局部波数的定义,然后推导出利用张量局部波数法进行反演的基本公式.本文方法在进行张量数据反演时无需事先知道场源体的类型(构造指数)即可获得场源体的位置信息,且可根据位置参数对场源体的类型进行估计.通过理论模型证明张量局部波数法可以很好地完成位场全张量数据的反演工作,并将其与常规局部波数法进行对比,证明全张量局部波数法的反演结果更加准确,即使在测点分布不合理的情况下,张量局部波数法仍可以获得准确的结果.最后应用张量局部波数法对美国得克萨斯州实测重力数据进行了反演,其反演结果与已有的研究成果相一致.
关键词: 位场全张量数据      张量局部波数      构造指数      反演     
Comparison of the tensor local wavenumber method with the conventional local wavenumber method for interpretation of total tensor data of potential fields
MA Guo-Qing, DU Xiao-Juan, LI Li-Li     
College of Geoexploration Science and Technology, Jilin University, Changchun 130026, China
Abstract: Measurements of potential field tensor(PFT)data have been widely used in the exploration of geophysics because of the advantages of big capacity, high-precision, and low noise. In this paper, we propose to use the tensor local wavenumber (TLW) method to interpret total tensor data of potential fields. We firstly define the tensor local wavenumber, and then derive the inversion formula based on the tensor local wavenumber. The TLW method does not require any prior information about the nature of the source to estimate the location parameters, and then the nature of the source can be obtained by the source location parameters. The TLW method is demonstrated on theoretical potential field data, and we compare it with the conventional local wavenumber (CLW) method. The results show that the inversion results of the TLW method are more accurate. Especially in the case that observational points are distributed unreasonable, the TLW method still can give accurate results. At last, we apply the TLW method to a gravity data example in Texas, USA, and the inversion results are in consistent with previous results..
Key words: Potential field tensor data      Tensor local wavenumber      Structural index      Inversion     
1 引 言

位场解释的主要目的是确定场源体的位置及其几何形状.为此人们提出了很多方法来解决这一问题[1-2],随着计算机和商业软件的普及,这些方法已经被广泛应用.场源参数法[3]是其中应用范围很广的一种方法,该方法是利用局部波数来进行场源体的反演,局部波数是表示局部相位(倾斜角)在不同方向上的变化率的量.Bracewell(1965)根据分析信号的概念首次定义了局部波数[4].局部波数法最早是用来获得磁异常体的位置参数,但需要事先知道场源体的类型(构造指数),在实际数据处理中,场源体的形状往往是难以获得的.为了避免这一问题,人们对局部波数法进行改进使其在不需要知道构造指数的情况下直接进行异常体深度的计算[5-7].Salem等(2005)利用局部波数及其相位转换形式直接估计场源的位置参数[8];Pilkington 和Keating 在2006年证明了局部波数为分析信号与其导数的比值[9],因此局部波数受磁化方向的影响较小;Salem等(2008)试验了局部波数及其相位转换形式在三维情况下的应用效果[10];Keating(2009)利用不同高度上的局部波数直接计算场源体的深度及类型[11].

全张量测量开始于20世纪90年代,是一种非常有效的地球物理探测手段[12-17],能够更好地描述异常体的特征,且受噪音的干扰小.现今主要有以下几种方法来进行全张量位场数据解释.Zhang等在2000年提出张量欧拉反褶积法进行全张量重力数据的解释工作[18],但该方法需要事先已知场源体的形状;SchmidtandClark(2006)讨论了磁张量的各个分量的性质及特征[19];Valentin(2007)利用张量反褶积法实现了全张量重力数据的解释工作[20];Beiki(2010)利用全张量重力数据的解析信号来完成全张量数据的反演工作[21],取得了很好的效果.BeikiandPedersen (2010)利用全张量重力数据的特征向量来估计理想地质体的位置[22];Qruc(2010)利用张量不变量的最大值去计算理想地质体的深度[23],其应用范围较窄,吴星(2011)采用张量不变量法解释非全张量卫星重力梯度数据[24];Qruc(2010)利用磁张量数据的解析信号来进行简单地质体的反演[25],但该方法需要事先知道场源体的类型;Beiki, PedersenandNazi(2011)利用特征向量法来进行航磁张量数据的解释[26],在应用中限制条件较多.

本文将局部波数的概念推广到三维空间,提出了张量局部波数法,该方法是利用xyz方向的局部波数来进行全张量数据的反演.本文首先给出了张量局部波数的定义,然后推导了张量局部波数与场源体位置参数之间的对应关系,最后通过所获得的位置参数进行场源体类型的计算.通过理论模型试验证明了张量局部波数法的有效性及相对常规局部波数法的优越性,最后将该方法应用于解释美国得克萨斯州实测重力数据,其反演结果与前人研究成果一致.

2 张量局部波数法

全张量数据是指位场U的二阶导数,其具体形式为

(1)

按照常规倾斜角[6]的定义,给出了张量数据在xyz三个方向上的倾斜角:

(2)

(3)

(4)

θxθyθz分别表示xyz方向的倾斜角.

本文定义不同方向的倾斜角在xyz方向的变化率为张量局部波数,其表达式分别为

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

其中,

(14)

(15)

(16)

Zhang等(2000)给出了张量欧拉反褶积法的基本表达式:

(17)

(18)

(19)

其中xyz表示观测点的坐标,x0y0z0 表示场源体的坐标,BxByBz为区域场在xyz方向的分量.根据常规局部波数法的推导过程[10]进行张量局部波数法的推导.分别求取方程(17)在xyz方向上的导数,由于背景场埋深较大,其导数随距离的衰减速率更快,可以忽略掉区域场的导数项,则获得如下方程:

(20)

(21)

(22)

对公式(20)和(22)分别乘以然后相减得

(23)

对公式(21)和(22)分别乘以,然后相减得:

(24)

(23) 式和(24)式相加后得到

(25)

对(25)式除以∂Ux/∂h并将(5)—(7)式代入后得到如下方程:

(26)

同理,根据方程(18)和(19)可以得知如下方程成立:

(27)

(28)

将线性方程(26)—(28)写成矩阵形式:

(29)

利用最小二乘法来求解方程(29)即可获得场源体的位置信息.在利用张量局部波数法进行求解后采用如下三个筛选准则[18]对初始结果进行筛选来获得更加准确的结果.

(1) 利用一个较小窗口滤除一些孤立的解;

(2) 当求解结果与窗口中心点之间的距离大于窗口半径时该解无效,因为当窗口在场源附近时解的准确性较高;

(3) 经过以上的筛选后,解已经比较集中,对于解的个数小于一定值的异常进行滤除,因为当解的个数小于一定数量时不认为是有效异常.经过上述的筛选后就已经获得比较准确的场源位置结果x0y0z0,根据该结果利用方程(20)—(22)来估计地质体的类型(构造指数),其具体计算公式如下:

(30)

通过求解方程(30)即可获得构造指数的大小.

3 理论模型试验

为了试验张量局部波数法的应用效果,建立如下的模型:在水平位置(25,35),(75,35)和(50,75)处分别存在埋深为10、12.5 和15 m, 半径为10 m的球体,当地磁倾角为60°,磁偏角为0°,计算点距为2m, 地质体在地面所引起的磁异常的张量局部波数如图 1所示.

图 1 球体在地面引起的磁异常的张量局部波数 Fig. 1 Tensor local wavenumber of the magnetic anomaly caused by the spheres

分别利用常规局部波数法和张量局部波数法对异常进行反演,窗口大小均为7×7,见图 2.

图 2 不同方法所反演得到的结果 (a)常规局部波数反演得到的异常体的位置;(b)张量局部波数反演得到的异常体的位置;(c)常规局部波数反演得到的深度分布直方图;(d)张量局部波数反演得到的深度分布的直方图;(e)常规局部波数反演得到的构造指数直方图;(f)张量局部波数反演得到的构造指数直方图.rad表示弧度. Fig. 2 The inversion results computed by different methods (a) The inversed locations by the conventional local wavenumber (CLW) method; (b) The inversed locations by the tensor local wavenumber (TLW) method;(c) Histogram of the inversed depths by the CLW method;(d)Histogram of the inversed depths by the TLW method;(e) Histogram of the inversed structural index by the CLW method;(f)Histogram of the inversed structural index by the TLW method.

图 2a图 2b中分别用+标识异常体的真实位置.图 2a为常规局部波数法反演结果的分布图,图 2c2e为反演深度和构造指数分布直方图,横坐标分别为深度和构造指数,纵坐标为一定区间内所包含解的数量.常规局部波数法反演得到的深度分别为9.80±1.1m、12.5±1.0m、14.8±0.8m, 与理论值存在较小的差距,反演得到的构造指数为2.95±0.17.图 2b为张量局部波数法反演结果的分布图,图 2d2f为反演深度和构造指数分布直方图.张量局部波数法所反演得到的球体深度分别为10.0±0.3m、12.5±0.5 m、14.9±0.4 m, 与理论深度相一致,反演得到的构造指数为2.98±0.08.从反演结果中可以看出,两种方法均能很好地完成异常的反演,其构造指数与球体理论构造指数(SI=3)几乎相等,但张量局部波数法的反演结果相对更加集中,与理论值之间的差距更小.从反演结果的分布图中还可以看出,两种方法均不受磁化方向的干扰.

为了验证方法对不同类型异常的应用效果,构建如下的重力模型:存在一个上顶埋深为20 m、边长为20m 的棱柱体和一个埋深为15 m、半径为10m的球体,计算点距为2 m, 地质体位置如图 3a所示,该模型引起重力异常的张量局部波数如图 3所示.

图 3 模型引起重力异常的张量局部波数 Fig. 3 Tensor local wavenumber of the gravity anoma lycaused by the model

分别利用常规局部波数法和张量局部波数法对该异常进行反演,窗口大小均为7×7,见图 4.

图 4 不同局部波数法所反演得到的结果 (a)常规局部波数反演得到的异常体的位置;(b)张量局部波数反演得到的异常体的位置;(c)常规局部波数反演得到的深度分布直方图;(d)张量局部波数反演得到的深度分布的直方图;(e)常规局部波数反演得到的构造指数直方图;(f)张量局部波数反演得到的构造指数直方图. Fig. 4 The inversion results by different local wavenumber method (a)The inversed locations by the CLW method; (b)The inversed locations by the TLW method; (c) Histogram of the inversed depths by the CLW method; (d) Histogramof the inversed depths by the TLW method ; e) Histogramof the inversed structural index by the CLW method; (f) Histogram of the inversed structural index by theTLW method.

图 4a图 4b中的黑框和+分别表示异常体的真实位置.图 4a、4c和4e为常规局部波数法的反演结果,反演得到的异常体的深度分别为14.7±0.4m(球体)、19.5±1.4 m(棱柱体),构造指数分别为0.22±0.60(棱柱体)、2.04±0.14(球体),与理论值存在一定差距,且常规局部波数法的反演结果在两个异常的连线上产生了一个多余的异常,其深度为16.8±1.0m, 这会为进一步的解释工作带来不必要的麻烦.图 4b、4d 和4f为张量局部波数法的反演结果,其反演得到的异常体的深度分别为15.0±0.4m(球体)、19.8±0.8 m(棱柱体),构造指数分别为0.21±0.22(棱柱体)、2.01±0.06(球体),反演结果与理论值相差较小,发散程度较低,且未产生多余的异常,可很好地描述异常体的真实特征.

为了验证方法的适应性,将上述重力模型的计算点距加大到5m, 其所引起的张量局部波数见图 5.

图 5 模型引起的重力异常的张量局部波数(点距5m) Fig. 5 Tensor local wavenumber of the gravity anomaly caused by the model (interval is 5 m)

利用常规局部波数法和张量局部波数法对该异常进行反演,反演窗口为3×3,其反演结果如图 6所示.

图 6 不同局部波数法反演得到的结果 (a)常规局部波数反演得到的异常体的位置;(b)张量局部波数反演得到的异常体的位置;(c)常规局部波数反演得到的深度分布直方图;(d)张量局部波数反演得到的深度分布的直方图;(e)常规局部波数反演得到的构造指数直方图;(f)量局部波数反演得到的构造指数直方图. Fig. 6 The inversion results by different local wavenumber methods (a)The inversed locations by the CLW methods (a)The inversed locations by the CLW method;(b)The inversed locations by the TLW method; (c) Histogram of the inversed depths by the CLW method;(d) Histogram of the inversed depths by the TLW method ; (e) Histogramof the inversed structural index by the CLW method; (f) Histogram of the inversed structural index by theTLW method.

图 6a图 6b中的黑框和+分别表示异常体的真实位置.从反演结果中可以看出,当加大测量点距时,常规局部波数法反演得到的深度与构造指数与理论值之间的误差加大,尤其是所反演得到的异常体的位置与真实位置有较大偏差,无法完成异常的解释;张量局部波数法在该情况下依旧能很好地完成异常的反演工作,其反演结果与理论值之间的误差依然很小,反演得到的异常体的位置与理论位置一致,发散程度较低.从对比中可以看出,张量局部波数法相对常规局部波数法更加稳定和准确.

4 实际重力数据应用

在将张量局部波数法应用于美国得克萨斯州的一个研究构造上所获得的重力数据[27]的处理中,其异常如图 7所示.

图 7 实测重力异常 Fig. 7 Measured gravity anomaly

利用FFT 变换[28]计算得到异常的全张量数据(见图 8).

图 8 重力全张量异常 Fig. 8 Gravity gradient tensor data

利用方程(5)—(13)计算该地区的张量局部波数见图 9.

图 9 该地区重力异常的张量局部波数 Fig. 9 Tensor local wavenumber of the gravity anomaly in the area

分别利用局部波数法和张量欧拉反褶积法对该数据进行处理(图 10).

图 10 张量局部波数法和张量欧拉反褶积法的反演结果 (a)常规局部波数法反演得到的结果;(b)张量局部波数法反演得到的结果;(c)张量欧拉反褶积法反演得到的结果;(d)常规局部波数反 演得到的深度分布直方图;(e)张量局部波数反演得到的深度分布直方图;(f)张量局部波数反演得到的深度的三维图;(g)常规局部波数 反演得到的构造指数直方图;(h)张量局部波数反演得到的构造指数直方图;(i)张量局部波数反演得到的构造指数三维图. Fig. 10 The inversion results by the tensor local wavenumber method and tensor Euler deconvolution (a)The inversed results by the CLW method;(b)The inversed results by the CLW method;(c)The inversed locati ons by the tensor Euler deconvolution method;(d) Histogram of the inversed depths by theCLW method ;(e) Histogram of the inversed depths by the TLW method;(f) 3D view of the inversed depths by the TLW method;(g) Hstogram of the inversed structural index by the CLW method;(h) Hstogram of the inversed structural index by the TLW method;(i) 3D view of the mversed structural index by the TLW method.

从结果中可以看出,常规局部波数法和张量局部波数法所得到的地质体(SI分别为1.92和1.94)近似为球体(SI=2),其平均深度分别为4.15 和4.14km.利用全张量欧拉反褶积法对该异常进行处理,构造指数选取为1.93,其反演得到的结果的地质体的平均深度为4.18km, 因此本文方法得到的深度与张量欧拉反褶积法相一致.将本文得到结果与前人解释结果[2729-33]进行比较见表 1.

表 1 张量局部波数法反演结果与先前计算结果比较 Table 1 comparison of the inversion results computed by the TLW method and previous results

从结果对比中可以看出,本文方法与ShawandAgarwal, Essa和BluentQruc计算得到的结果相互一致.但是本文方法与其它结果还是有一些差距,这可能是由于本文的全张量数据是计算出来的,而不是实际测量得到的.

5 结 论

本文提出利用张量局部波数法来进行全张量位场数据的解释,张量局部波数是指全张量数据在不同方向上的波数,本文推导出了利用张量局部波数进行异常反演的基本公式.该方法的优势是不需要事先知道场源体的类型就可以获得地质体的位置信息,然后可通过反演得到的位置信息对异常体的类型进行计算.通过理论模型试验本文方法可以很好地完成位场全张量数据的解释工作,且通过与常规局部波数法的比较可以看出,张量局部波数法的反演结果更加稳定,更加准确.最后将其应用于美国得克萨斯州盐丘的重力异常的解释,获得结果与前人解释结果相一致.

致谢

特别感谢提供实际资料的同仁以及为评审本文所付出努力的专家们

参考文献
[1] Blakely R J. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press, 1995 .
[2] Nabighian M N, Asten M. Metalliferous mining geophysics: State of the art in the last decade of the 20th century and the beginning of the new millennium. Geophysics , 2002, 67(3): 964-978. DOI:10.1190/1.1484538
[3] Thurston J B, Smith R S. Automatic conversion of magnetic data to depth, dip, and susceptibility contrast using the SPI method. Geophysics , 1997, 62(3): 807-813. DOI:10.1190/1.1444190
[4] Bracewell R N. The Fourier Transform and Its Application. New York: McGraw-Hill Book Co., 1965. http://cn.bing.com/academic/profile?id=1628671557&encoded=0&v=paper_preview&mkt=zh-cn
[5] Smith R S, Thurston J B, Dai T, et al. iSPI-The improved source parameter imaging method. Geophysical Prospecting , 1998, 46: 141-151. DOI:10.1046/j.1365-2478.1998.00084.x
[6] Thurston J B, Smith R S, Guillon J. A multi-model method for depth estimation from magnetic data. Geophysics , 2002, 67(2): 555-561. DOI:10.1190/1.1468616
[7] Salem A, Smith R S. Depth and structural index from the normalized local wavenumber of 2D magnetic anomalies. Geophysical Prospecting , 2005, 51(1): 83-89.
[8] Salem A, Ravat D, Smith R, et al. Interpretation of magnetic data using an enhanced local wavenumber (ELW) method. Geophysics , 2005, 70(2): L7-L12. DOI:10.1190/1.1884828
[9] Pilkington M, Keating P. The relationship between local wavenumber and analytic signal in magnetic interpretation. Geophysics , 2006, 71(1): L1-L3.
[10] Salem A, Williams S, Fairhead D, et al. Interpretation of magnetic data using tilt-angle derivatives. Geophysics , 2008, 73(1): L1-L10. DOI:10.1190/1.2799992
[11] Keating P. Improved use of the local wavenumber in potential-field interpretation. Geophysics , 2009, 74(6): L75-L85. DOI:10.1190/1.3242270
[12] Jekeli C. A review of gravity gradiometer survey system data analyses. Geophysics , 1993, 58(4): 508-514. DOI:10.1190/1.1443433
[13] Bell R E, Anderson R, Pratson L. Gravity gradiometry resurfaces. The Leading Edge , 1997, 16(1): 55-60. DOI:10.1190/1.1437431
[14] Schmidt P W, Clark D A. Advantages of measuring the magnetic gradient tensor. Preview , 2000, 85: 26-30.
[15] Schmidt D V, Bracken R E. Field experiments with the tensor magnetic gradiometer system at Yuma Proving Ground, Arizona. Proceedings of the Symposium on the Application of Geophysics to Engineering and Environmental Problems. (SAGEEP), 2004. http://www.academia.edu/7065885/Location_and_depth_estimation_of_point-dipole_and_line_of_dipoles_using_analytic_signals_of_the_magnetic_gradient_tensor_and_magnitude_of_vector_components
[16] Gamey T J, Doll W E, Beard L P. Initial design and testing of a full-tensor airborne SQUID magnetometer for detection of unexploded ordnance. SEG Expanded Abstracts , 2004, 23: 798.
[17] Doll W E, Gamey T J, Beard L P, et al. Airborne vertical magnetic gradient for near-surface applications. The Leading Edge , 2006, 25(1): 50-53. DOI:10.1190/1.2164755
[18] Zhang C, Mushayandebvu M F, Reid A B, et al. Euler deconvolution of gravity tensor gradient data. Geophysics , 2000, 65(2): 512-520. DOI:10.1190/1.1444745
[19] Schmidt P W, Clark D A. The magnetic gradient tensor: Its properties and uses in source characterization. The Leading Edge , 2006, 25(1): 75-78. DOI:10.1190/1.2164759
[20] Valentin M, Gwendoline P, Michel D, et al. Tensor deconvolution: A method to locate equivalent sources from full tensor gravity data. Geophysics , 2007, 72(5): 161-169.
[21] Beiki M. Analytic signals of gravity gradient tensor and their application to estimate source location. Geophysics , 2010, 75(6): 159-174.
[22] Beiki M, Pedersen L B. Eigenvector analysis of gravity gradient tensor to locate geologic bodies. Geophysics , 2010, 75(6): 137-149. DOI:10.1190/1.3503577
[23] Qruc B. Depth estimation of simple causative sources from gravity gradient tensor invariants and vertical component. Pure and Applied Geophysics , 2010, 167(10): 1259-1272. DOI:10.1007/s00024-009-0021-4
[24] 吴星, 王凯, 冯炜, 等. 基于非全张量卫星重力梯度数据的张量不变量法. 地球物理学报 , 2011, 54(4): 966–976. Wu X, Wang K, Feng W, et al. Method of tensor invariant based on non-full tensor satellite gravity gradients. Geophys. (in Chinese) (in Chinese) , 2011, 54(4): 966-976.
[25] Qruc B. Location and depth estimation of point-dipole and line of dipoles using analytic signals of the magnetic gradient tensor and magnitude of vector components. Journal of Applied Geophysics , 2010, 70(1): 27-37. DOI:10.1016/j.jappgeo.2009.10.002
[26] Beiki M, Pedersen L B, Nazi H. Interpretation of aeromagnetic data using eigenvector analysis of pseudo gravity gradient tensor. Geophysics , 2011, 76(3): L1-L10. DOI:10.1190/1.3555343
[27] Nettleton L L. Gravity and Magnetics in Oil Prospecting. New York: McGraw-Hill Book Co., 1976. http://cn.bing.com/academic/profile?id=95449890&encoded=0&v=paper_preview&mkt=zh-cn
[28] Mickus K L, Hinojosa J H. The complete gravity gradient tensor derived from the vertical component of gravity: A Fourier transform technique. Journal of Applied Geophysics , 2001, 46(3): 159-174. DOI:10.1016/S0926-9851(01)00031-3
[29] Abdelrahman E M, Bayoumi A I, Elaraby H M. Short note: A least-squares minimization approach to invert gravity data. Geophysics , 1991, 56: 115-118. DOI:10.1190/1.1442946
[30] Shaw R K, Agarwal B N P. A generalized concept of resultant gradient to interpret potential field maps. Geophysical Prospecting , 1997, 45(3): 513-520. DOI:10.1046/j.1365-2478.1997.450279.x
[31] Salem A, Ravatz D, Mushayandebvu M F, et al. Linearized least-squares method for interpretation of potential-field data from sources of simple geometry. Geophysics , 2004, 69(3): 783-788. DOI:10.1190/1.1759464
[32] Essa K S. Gravity data interpretation using the s-curves method. Geophys. Eng , 2007, 4(2): 204-213. DOI:10.1088/1742-2132/4/2/009
[33] Mohan N L, Anandababu L, Roa S. Gravity interpretation using Mellin transform. Geophysics , 1986, 52(1): 114-122.