地球物理学报  2016, Vol. 59 Issue (6): 2225-2231   PDF    
引用本文
陈国强, 马国庆. 2016. 位场数据解释的Theta-Depth法. 地球物理学报, 59(6): 2225-2231, doi: 10.6038/cjg20160625.  
CHEN G Q, MA G Q. 2016. Theta-Depth method for the interpretation of potential field data. Chinese J. Geophys. (in Chinese), 59(6): 2225-2231, doi: 10.6038/cjg20160625.  
位场数据解释的Theta-Depth法
陈国强1 , 马国庆2     
1. 吉林大学地球科学学院, 长春 130061;
2. 吉林大学地球探测科学与技术学院, 长春 130021
收稿日期 2015-06-02, 2016-03-16 收修定稿
基金项目: 国家重点基础研究发展计划(973计划)(2015CB453000)、国家自然科学基金(41404089)和中国地质调查局地质矿产调查评价专项项目(GZH003-07-03)联合资助.
作者简介: 陈国强,男,1985年生,博士研究生,主要从事矿产资源定量与定位预测研究.E-mail:chenguoqiang0506@163.com
通讯作者: 马国庆,男,1984年生,讲师,主要从事重磁数据解释方法研究.E-mail:maguoqing@jlu.edu.cn
摘要: Theta图是利用位场(重磁)数据识别边界的常用方法,其表达式为重磁异常水平变化与垂直变化的比值函数.该方法计算浅源地质体边界的效果较好,而由于深源位场数据在换算过程中会产生趋同效应,在深源地质体识别应用中计算结果不准确,为此,本文提出Theta-Depth法并进行地质体埋深的计算.首先给出直接利用Theta图像进行场源体深度估算的方法,然后推导出基于Theta导数的线性方程来自动估算场源位置参数,本文方法可有效地利用Theta图像的特征为约束条件来提高反演结果的精度.理论模型试验证明本文提出的Theta-Depth法能有效地计算出场源体位置和深度.将该方法应用于满都拉地区实测磁数据的解释,帮助圈定了矿脉的分布.
关键词: Theta      位场      深度      线性方程     
Theta-Depth method for the interpretation of potential field data
CHEN Guo-Qiang1, MA Guo-Qing2     
1. College of Earth Science, Jilin University, Changchun 130061, China;;
2. College of Geoexploration Science and Technology, Jilin University, Changchun 130021, China
Abstract: The Theta map is a commonly used tool to detect edges of subsurface geological bodies using potential field (gravity and magnetic) data, which is the ratio of function of horizontal derivative to vertical derivative of gravity or magnetic anomalies. Although the effect of calculation of shallow geological body edges by this method is good, it is not successful in probing the deep geological body edges process due to the convergence of the conversion process that leads inaccurate calculation. To solve this problem, we proposed the Theta-Depth method to compute the depth of the buried geological body. Frist, we directly used the Theta map to estimate the depth of the source. Then we derived a linear equation based on the derivatives of Theta to compute the location and depth of the source. The presented method can use the Theta map as a constraint to obtain more accurate results effectively. The tests show that this Theta-Depth method can effectively calculate the location and depth of the source. We applied the Theta-Depth method to interpretation of measured magnetic data in the Mandula area, which helped delineate the distribution of ore veins..
Key words: Theta      Potential field data      Depth      Linear equation     
1 引言

边界识别是位场数据解释的常规任务之一，现今有很多方法来完成这一任务(Evjen，1936; Cordell and Grauch，1982; Blakely，1995).重磁异常水平导数的最大值、垂直导数的零值与地质体的边界相对应(Roest et al.，1992; Miller and Singh，1994; Blakely，1995)，因而边界识别方法大多是水平与垂直导数的函数.前期边界识别方法仅能给出较浅地质体的边界，而不能给出深部地质体的边界.为了改善这一缺点，人们开始致力于均衡边界识别滤波器的研究(Miller and Singh，1994; Verduzco et al.，2004; Wijns et al.，2005; Cooper and Cowan，2006，2008; Cooper，2009; 马国庆等，2012； Ma，2013)，该类方法能同时识别出不同深度地质体的边界.倾斜角法(Miller and Singh，1994)是第一个提出的均衡边界识别方法，但该方法所得到的边界并不十分清晰.Theta图(Wijns et al.，2005)是一种常用的均衡边界识别滤波器，其能清晰地给出场源体的水平位置.Cooper和Cowan(2006)将该方法与倾斜角法进行比较，并将其应用于实际数据的边界识别任务.马国庆等(2012)将该方法应用于四川盆地线性构造的划分，断裂结果与实际地质资料吻合较好.Ma和Li(2012)总结Theta图法的应用特征，并用于中国内蒙古航磁数据的构造划分.

地质体位置反演是重磁数据解释的主要任务，随着勘探数据量的增加，人们更加倾向于选择自动解释方法(Debeglia and Corpel，1997; Thurston and Smith，1997; Hsu et al.，1998; Smith et al.，1998; Kirkham，2001; Salem et al.，2004; Salem and Smith，2005)来进行计算，该类方法在给定异常后可直接计算出地质体的位置参数，主要包括主要包括欧拉反褶积法(Thompson，1982; Reid et al.，1990)、局部波数法(Salem et al.，2008)和解析信号法(Kirkham，2001)，通过原始异常构建相应函数于位置参数的对应方程，进而完成计算，在实际数据解释中应用较为广泛.欧拉反褶积法是由Thompson在1982年提出应用于磁法数据的解释，随后被应用于重力数据的解释，解析信号和局部波数法在进行磁异常解释时具有不受倾斜磁化干扰的优点.众所周之，方程求解会造成结果的发散效应，因而以上方法在求解后需设定结果筛选准则来获得更加准确的结果，由于以上方法的函数本身与地质体水平位置之间不具有直接的对应关系，往往需要借助其他函数来进行，从而造成结果之间缺乏对应性.

本文提出Theta-Depth法计算异常体的深度信息，首先推导出一种直接利用Theta进行深度计算的方法，通过理论推导证明Theta反正弦函数两条45°等值线之间距离的一半等于地质体的埋深，利用此关系进行深度的估算.此外，本文还推导出一种基于Theta的自动解释方法，其利用Theta导数组成的线性方程来完成场源体位置和深度的计算，并根据Theta函数与地质体位置之间的关系设定相应的筛选准则.

2 Theta-Depth法

Theta是由Wijns于2005年提出的一种进行边界识别的方法，其表达式为

(1)

Theta图的最大值对应着地质体的边界，该方法具有能够同时识别出不同深度地质体边界的优势，在实际数据解释中应用较为广泛.

Nabighian(1972)给出水平位置为0埋深为z的接触带所产生磁异常的水平与垂直导数的表达式为

(2)

(3)

其中，K为磁化率差，F为当地的地磁场强度，c=1-cos2isin2A，A为h轴与磁北方向的夹角，i为磁倾角，tanI=tani/cosA，d为地质体的倾斜角度.

在垂直磁化情况下式(2)和式(3)可以改写为

(4)

(5)

将式(4)和式(5)带入到式(1)，整理后可以得到

(6)

通过式(6)可以看出，当h=0时式(6)的反正弦函数取得最大值90°，因此可通过式(6)反正弦的最大值判断出地质体的水平位置；当h=z时式(6)反正弦函数为45°，随着h的增大式(6)的反正弦函数继续减小直至接近最小值0.通过以上的分析可知，利用Theta反正弦函数两条45°等值线之间距离的1/2可估算出地质体的埋藏深度.需要指出的是，利用45°等值线之间距离计算深度只是其中一种计算深度的方法，任何的角度都可以建立起水平距离与埋深的关系.

本文还推导出一种利用Theta导数组成的线性方程进行位场异常反演的方法.常规欧拉反褶积法的表达式为

(7)

其中，x、y和z为观测点坐标，x0、y0和z0为待求的场源体的空间位置坐标，B为背景场值，N为表征场源体类型的构造指数.计算式(7)在x和y方向上的导数，由于背景异常变化较为平缓，其导数较小可忽略，因此式(7)导数的表达式为

(8)

(9)

分别对式(8)和式(9)乘以和，相加后得到

(10)

其中，THD为总水平导数，，A为梯度模，.

Huang等(1995)证明解析信号的欧拉反褶积方程的表达式为

(11)

式(11)乘以，并和式(10)相减后得到

(12)

计算式(1)在x，y和z方向上的导数为

(13)

(14)

(15)

将式(13)、(14)和(15)带入到式(12)可以得到

(16)

根据式(16)可知，利用Theta的导数所组成的线性方程可估算出场源体的位置坐标x0、y0和z0.方程求解往往会造成结果的发散，因此利用Theta函数极大值与场源体之间的对应关系设定如下的结果筛选准则.

(1)距离Theta极大值超过窗口长度的解去掉.

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

(3)经过以上的筛选后，解应该已经比较集中，对于一定范围内解的总数小于给定值的异常进行滤除，因为当解的个数小于一定数量时不认为是有效异常.

3 模型试验

为验证上述算法的精确性与实用性，本文进行了模型试验，主要分为原始异常和含噪异常来进行试验，并分别针对重磁异常进行了试验来验证方法的稳定性和适应性.

由于Theta法的结果会受到倾斜磁化的干扰，因此利用理论模型产生的磁异常和重力异常来试验Theta-Depth方法的应用效果.图 1a垂直磁化条件下埋深分别为5 m和7 m正方体所产生的磁异常.图 1b为磁异常Theta计算结果的反正弦，其中45°等值线采用虚线表示，并用灰色进行填充之间的区域.通过公式(6)的推导证明地质体的埋深等于Theta反正弦结果的两条45°等值线间距离的一般，其区域极大值标识地质体的水平范围.根据图 1b利用两条45°等值线(虚线)的垂直距离可计算出地质体的埋深分别为4.86 m和6.72 m，与地质体理论埋深差距较小，说明该方法具有良好的应用效果.

图 1 (a) 垂直磁化情况下正方体所产生的磁力异常， (b) 异常Theta的反正弦计算结果，(c) 利用式(16)反演得到的场源体的水平位置， (d) 反演结果的三维显示 Fig. 1 (a) Magnetic anomalies of a cub in the vertical magnetization， (b) Arcsine calculation of anomalies using Theta method， (c) Horizontal position of source by inversion， (d) Three-dimensional display of inversion results

利用公式(16)计算地质体的埋深，其窗口大小为7×7，并采用本文给出的结果筛选准则来获得更加准确的结果.图 1c为利用Theta导数组成线性方程所计算得到的异常体的水平位置分布，其结果与第一种方法所得到的异常体的水平位置相一致，且与地质体的理论位置吻合较好.图 1d为Theta导数线性方程反演结果的三维显示，地质体的埋深分别为4.91 m和6.87 m，与理论埋深相接近.通过该试验证明Theta-Depth法能准确地计算出地下异常体的深度.

下面试验一下Theta-Depth方法在含噪重力异常解释中的应用效果.图 2a为埋深分别为7 m和10 m正方体所产生的重力异常，并加入了均值为0，方差为1 mGal的随机噪声.图 2b为重力异常的Theta结果，其最大值能准确地标识出地质体的水平位置.图 2c为Theta的反正弦计算结果，并利用绿色填充45°等值线间的区域，根据45°等值线之间的距离可以得到异常体的埋深分别为6.81 m和9.66 m.

图 2 (a) 正方体产生的重力异常并加入随机噪声后异常; (b) 图2a所示异常的Theta结果； (c)Theta的反正弦结果; (d) Theta-Depth法反演得到的异常体的水平位置; (e) 筛选后的反演结果; (f) 反演结果的三维显示 Fig. 2 (a) Gravity anomalies of the cube with added random noise; (b) Theta results of anomalies in (a);(c)Arcsine results by Theta method; (d) Horizontal position of abnormal body obtained by inversion of the Theta-Depth method; (e) Inversion results after screening; (f) Three-dimensional display of inversion results

图 2d为利用公式(16)计算得到的结果，其窗口大小为7×7.从结果中可以看出由于噪声的加入产生了较多的干扰异常.采用上述筛选准则来获得更加准确的结果，图 2e为方程反演结果经过筛选后得到的结果筛，可以看出该筛选准则很好地去除了随机噪声的干扰.图 2f为筛选后反演结果的三维显示，异常体的深度为6.87 m和9.73 m，与地质体的真实埋深相接近.

4 实际数据应用效果

将该方法应用于满都拉地区实测磁异常数据的解释中，磁数据的采集点距为40 m，该地区主要是为了寻找地下铁矿场的位置，周围矿产的开采深度在100 m左右.图 3a为化极后磁异常，可以看出该地区构造呈东南向展布.图 3b为磁异常的Theta计算结果，其最大值给出了了地下矿体及地层之间的界线.图 3c为Theta的反正弦结果，并采用绿色对45°等值线间进行填充，根据距离可估算矿脉的深度约为140 m.图 3d为利用式(16)反演得到的异常体的分布，并对结果进行筛选.从结果中可以看出异常体的深度大多集中在80~200 m之间，与另一种Theta-Depth方法计算得到的结果相一致，说明本文提出的Theta-Depth法具有较好的实际应用效果，且反演深度与以往矿产深度相接近.

图 3 (a)满都拉地区化极后磁异常; (b) 磁异常的Theta结果; (c) Theta反正弦结果; (d) 式(16)反演得到的结果 Fig. 3 (a) Magnetic anomalies after reduction to the pole in Mandula area; (b) Theta results of magnetic anomalies; (c) Arcsine results by Theta method; (d) Inversion results by Eq.16
5 结论

地质体的埋深确定是重磁数据处理与解释的主要任务，在混成性地质空间分析中，确定地质体的水平边界与垂向边界对于地质体数字特征研究十分重要，尤其是有关含矿地质体埋深推断事关资源预测精度水平，本文提出的Theta-Depth法较好的完成了上述任务.Theta-Depth法首次将Theta函数进行扩展来完成异常体深度的计算，通过推导证明根据Theta反正弦函数可直接给出地质体的埋深，此外还推导出基于Theta导数的线性方法来估算地质体的水平位置和深度，且该方法在进行反演前不需要给定任何的先验信息.通过理论重力和磁数据证明本文提出的Theta-Depth方法能有效地完成重磁异常的解释工作，且反演结果与理论结果相一致.利用噪声重力异常试验可以看出本文方法依旧可以获得较为准确的结果，且对结果进行筛选后所得到的结果与理论值相一致.将本文方法应用于满都拉地区实际磁异常的解释，获得了岩脉的大致深度.

参考文献
Blakely R J. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press, 1995 .
Cooper G R J. 2009. Balancing images of potential-field data. Geophysics , 74(3): L17–L20.
Cooper G R J, Cowan D R. 2006. Enhancing potential field data using filters based on the local phase. Computers & Geosciences , 32(10): 1585–1591.
Cooper G R J, Cowan D R. 2008. Edge enhancement of potential-field data using normalized statistics. Geophysics , 73(3): H1–H4.
Cordell L, Grauch V J S. 1982. Mapping basement magnetization zones from aeromagnetic data in the San Juan basin, New Mexico.//SEG Technical Program Expanded Abstracts. SEG, 246-247.
Debeglia N, Corpel J. 1997. Automatic 3-D interpretation of potential field data using analytic signal derivatives. Geophysics , 62(1): 87–96.
Evjen H M. 1936. The place of the vertical gradient in gravitational interpretations. Geophysics , 1(1): 127–136.
Hsu S K, Coppens D, Shyu C T. 1998. Depth to magnetic source using the generalized analytic signal. Geophysics , 63(6): 1947–1957.
Huang D, Gubbins D, Clark R A, et al. 1995, Combined Study of Euler's Homogeneity Equation for Gravity and Magnetic field. 57th EAGE conference. Glasgow UK Extended Abstracts, 144.
Kirkham K. 2001. Investigations of a high-resolution aeromagnetic survey over the southeastern portion of the Illinois Basin[Master's thesis]. Carbondale: Southern Illinois University of Carbondale.
Ma G Q. 2013. Edge detection of potential field data using improved local phase filter. Exploration Geophysics , 44(1): 36–41.
Ma G Q, Huang D N, Yu P, et al. 2012. Application of improved balancing filters to edge identification of potential field data. Chinese Journal of Geophysics , 55(12): 4288–4295. doi: 10.6038/j.issn.0001-5733.2012.12.040.
Ma G Q, Li L L. 2012. Edge detection in potential fields with the normalized total horizontal derivative. Computers & Geosciences , 41: 83–87.
Miller H G, Singh V. 1994. Potential field tilt-a new concept for location of potential field sources. Journal of Applied Geophysics , 32(2-3): 213–217.
Nabighian M N. 1972. The analytic signal of two-dimensional magnetic bodies with polygonal cross-section; its properties and use for automated anomaly interpretation. Geophysics , 37: 507–517.
Reid A B, Allsop J M, Granser H, et al. 1990. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics , 55(1): 80–91.
Roest W R, Verhoef J, Pilkington M. 1992. Magnetic interpretation using the 3-D analytic signal. Geophysics , 57(1): 116–125.
Salem A, Ravat D, Mushayandebvu M F, et al. 2004. Linearized least-squares method for interpretation of potential-field data from sources of simple geometry. Geophysics , 69(3): 783–788.
Salem A, Smith R S. 2005. Depth and structural index from normalized local wavenumber of 2D magnetic anomalies. Geophysical Prospecting , 53(1): 83–89.
Salem A, Williams S, Fairhead D, et al. 2008. Interpretation of magnetic data using tilt-angle derivatives. Geophysics , 73(2): L1–L10.
Smith R S, Thurston J B, Dai T F, et al. 1998. iSPITM-The improved source parameter imaging method. Geophysical Prospecting , 46(2): 141–151.
Thompson D T. 1982. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data. Geophysics , 47(1): 31–37.
Thurston J B, Smith R S. 1997. Automatic conversion of magnetic data to depth, dip, and susceptibility contrast using the SPI (TM) method. Geophysics , 62(3): 807–813.
Verduzco B, Fairhead J D, Green C M, et al. 2004. New insights into magnetic derivatives for structural mapping. The Leading Edge , 23(2): 116–119.
Wijns C, Perez C, Kowalczyk P. 2005. Theta map: edge detection in magnetic data. Geophysics , 70(4): L39–L43.
马国庆, 黄大年, 于平, 等. 2012. 改进的均衡滤波器在位场数据边界识别中的应用. 地球物理学报 , 55(12): 4288–4295.