2. 中国地质科学院深部探测与地球动力学开放实验室, 北京 100037;
3. 中国地质科学院矿产资源研究所, 北京 100037
2. Earthprobe and Geodynamics Open Laboratory, Chinese Academy of Geological Sciences, Beijing 100037, China;
3. Institute of Mineral Resources, Chinese Academy of Geological Sciences, Beijing 100037, China
目前人们已经研究开发了多种基于利用磁异常导数的快速反演方法,用这些方法可以确定场源的参数,例如边界的位置和埋藏深度.解析信号法(AnalyticSignalMethod, 也称总梯度模法)由于解析信号具有不受(二维)或少受磁化方向影响,能够较好地反映磁性体边界,利用解析信号模确定磁源深度和产状,将其应用于地质构造“填图"已受到人们的重视[1-6].欧拉反褶积法[7-10]以Euler齐次方程为理论基础,利用总场及其梯度,依据欧拉齐次方程组确定场源的位置和深度,并可以得出形状因子(StructuralIndex, 又译为构造指数)对地质体进行识别,方法具有较强的适应性.但由于欧拉反褶积法采用窗口滑动计算,可以得到一系列深度点,在众多的深度点中如何分辨和确定有效的深度点一直是困扰人们的一个难题.
Salem 等(2003,2002)将Euler方程与解析信号表达式联立,推导出ANEUL 反演方法[11, 12],该方法可以同时计算磁性体的几何形状、水平位置和埋藏深度.但是用该方法需要计算磁异常的三阶导数,我们知道,计算导数相当于高通滤波.计算三阶导数不可避免地引进高频噪声干扰,而这种高频噪声干扰在实际资料反演计算中难以压制.为了改进上述方法,Salem 等(2008)又提出Tiltangle 方法[13],该方法虽然只需要计算一、二阶导数,但是需要确定磁性体走向方向,并沿垂直磁性体走向方向进行插值,然后通过求解超定线性方程组,获得场源的相应参数.这种方法利用多个测点上的场值进行反演计算,较利用单点反演压制随机噪声能力强,但是需要逐个确定磁性体走向,并要求解超定线性方程组,计算量较大.
本文提出同时计算磁性体形状和埋深等参数的低阶欧拉反褶积与解析信号相结合的反演方法.该方法在反演中只需计算磁异常的一、二阶导数,这样将高频噪声的干扰减少到人们可以控制的水平.
2 原理Reid等[8]给出Euler齐次方程如下:
(1) |
其中,T为位于点(x0,y0,z0)的磁性体在点(x,y,z)处产生的磁异常,n为形状因子,B为正常场存在的偏差.
计算式(1)对x的偏导数,经整理得
(2) |
令
x= x0,y= y0,z=0,埋深z0=h,方程(2)可以化简为
(3) |
同理,有
(4) |
(5) |
hTzz(0)= (n+1)Tz(0),(5)对(3),(4),(5)两边平方,加后开方,得
(6) |
其中,
向上延拓高度Δz时,即令x= x0,y= y0,z=-Δh,埋深h=z0,方程(3)可以化简为
(7) |
同理,有
(8) |
(9) |
对(7),(8),(9)两边平方,相加后开方,整理后得
(10) |
将式(6)与式(10)联立,化简后得
(11) |
将式(11)代入式(6),化简后得
(12) |
其中,
在正常场存在偏差的情况下,我们已知解析信号和一阶增强解析信号的模的极大值在地表和向上延拓高度为z1 的值,利用(11)和(12)式就可以求得磁性体的形状因子和埋深.
由(11)和(12)式可以看出,用该方法需要计算磁异常的一阶和二阶导数,我们知道,计算导数相当于高通滤波,不可避免地引进高频噪声干扰,从而影响计算结果.因此为了获得好的反演结果,在计算中需要压制高频噪声干扰.
3 模型实验 3.1 二度磁性体模型的计算为了验证方法的可靠性,作者分别对水平圆柱体、无限延深薄板和无限延深台阶三种二度磁模型沿垂直磁性体走向剖面上的磁异常进行了反演计算.作者计算了磁异常的导数、解析信号、形状因子和埋藏深度,并绘制了相应的曲线(见图 1~ 图 3).上述三种模型计算结果列于表 1.从表 1 可以看出计算结果的误差很小,说明计算公式(11)和(12)是正确的.
从图 1~图 3可以看出解析信号的极大值位于磁性体在地面的投影位置.离该极大值点越远,计算形状因子和埋深的误差越大.这是由于本方法计算的假设前提是解析信号的极大值位于磁性体中心的上方.
3.2 随机噪声对计算精度的影响像所有的梯度定量解释方法一样,随机噪声对欧拉反褶积与解析信号相结合方法的计算精度有一定的影响.由于随机噪声位于高频段,因此可以用向上延拓的方法减少噪声的影响.本文对上面讨论的水平圆柱体的理论异常加上一组幅值为磁异常的2%~8%的随机噪声.通过向上延拓不同高度,延拓高度由0.2km 增加到1.8km, 计算相应的形状因子和埋藏深度(见图 4 和图 5).由图 4 和图 5 可以看出,延拓高度较小时,计算的误差大.随着延拓高度的增大,计算的误差减小,当延拓超过一定高度时(在本例中,大于1.4km)计算结果是准确的.
为了研究邻近磁性体干扰的影响,我们计算了相距一定水平距离的垂直磁化无限延深薄板和水平圆柱体模型的埋深和形状因子(见图 6).用本文提出的方法计算了不同距离间隔的埋深h和形状因子n,并计算了其相对误差.为了便于分析,绘制了磁性体的几何参数埋深h和形状因子n的计算误差与间隔L与埋深h的比值之间的关系曲线(见图 7).从图 7可以看出,相邻磁性体的干扰所产生的计算相对误差具有随它们之间的间隔增大而减小的趋势.当间隔与深度的比值大于11时,几何参数的计算误差小于4%.需要说明的是,几何参数的计算误差与磁性体的相对大小规模及形状有关.Hsu等[4]曾用所提出的算法研究垂直磁化无限延深的台阶和厚板叠加模型的参数计算误差关系证实了这一点.
与上述模型不同,对于球体,只有在垂直磁化的情况下,解析信号的极大值位于磁性体中心的上方[11].而长方体等非球体,理论上不能严格符合欧拉方程条件.斜磁化长方体磁模型解析信号的形态不具对称性[6].实际中非理想模型多见,但是在一定的条件下是可以用本文提出的方法进行反演计算的.本文选用一系列埋深不同的垂直磁化直立长方体磁性模型,用过模型体中心剖面的计算结果进行分析和讨论.直立长方体磁性模型的顶面的宽度为20km, 长度为20km.延深为50km, 埋藏深度为1~5km, 间隔为0.5km.图 8 为埋藏深度为2km的直立长方体的磁异常平面图.图 9 为过直立长方体(图 8)的中心剖面磁异常和解析信号图.由图 9可以看出,在直立长方体边缘上方出现类似于无限延深厚板的解析信号极大值.我们用前面提出的方法计算磁性体的参数.
不同深度垂直磁化直立长方体的计算结果(见图 10和图 11)表明,随着埋深增大,计算结果的误差增加.当埋深不大于2km 时,埋深计算误差小于10.5%,解析信号极大值点相对于直立长方体边缘位移了一倍点距的距离,形状因子为0.028~0.013,近于无限延深台阶的特征.
庐江-枞阳陆相火山岩盆地位于长江中下游前陆区,即大别-苏鲁超高压碰撞造山带前陆缩短带.地处郯-庐断裂带的南段[14-16].庐枞盆地为一“继承式"或“基底坳陷型"火山盆地,呈耳状,长轴约56km, 短轴约24km, 面积约1032km2.庐枞盆地磁异常大致呈NE 向展布,北东宽,南西窄(见图 12).前人对该地区岩石磁性的研究表明:区内上元古界新生界岩层磁性很弱,属微磁或无磁性岩石.中碱性侵入岩及潜火山岩具中等到强磁性.庐枞盆地磁异常主要由这些中碱性岩类引起[17, 18].
由于斜磁化球体的解析信号极大值点偏离球心在地面的投影位置[11],因此我们首先对庐枞盆地的磁异常进行化极计算(见图 13).由于庐枞盆地的磁异常深部为规模大的正长岩基[18],因此我们对化极磁异常进行磁场分离.用向上延拓5km 的磁场作为区域场,用化极磁异常减去区域场得到局部磁异常(见图 14).先对分离后的中碱性侵入岩体及潜火山岩引起的局部磁异常向上延拓1km 和1.5km;再计算上延后磁异常的解析信号和一阶增强解析信号的模,搜索出它们的极大值;最后按公式(11),(12)计算各极大值点处的形状因子和埋深.反演计算结果(见图 15)表明:庐枞盆地内部磁性体的形状因子大小范围为0.35~2.91,即中碱性侵入岩体及潜火山岩的形状近似为无限延深的板状体和水平圆柱体、球体.板状体顶面埋深为0.3~1.2km.水平圆柱体和球体的中心埋深为1.6~2.9km.
与传统的位场反演方法不同,本文提出的欧拉反褶积与解析信号相结合的位场反演方法可以同时计算磁性体的形状因子和埋深.即用该方法可以确定磁性体的形状,并判断磁性体在地下的产出状态.这个方法在计算中仅仅需要计算磁异常的一、二阶导数,这样将高频噪声的干扰减少到人们可以控制的水平.随机噪声干扰的实验证实了这一点.相邻磁性体的干扰所产生的计算误差随它们之间的间隔增大而减小,干扰误差水平与磁性体的相对大小规模及形状有关.实际中非理想模型是多见的,对于垂直磁化直立长方体磁模型的计算结果表明,在一定的条件下是可以用本文提出的方法进行反演计算的.对安徽庐枞火山岩盆地磁异常反演计算的结果证实了该方法的实用性.
本方法也可以应用于重力异常的反演计算中,计算时需要注意形状因子的不同.根据泊松公式,磁异常与重力异常存在一阶导数的关系,因此前者的形状因子比后者的大1[19].
[1] | Nabighian M N. The analytic signal of two-dimensional magnetic bodies with polygonal cross-section: its properties and use for automated anomaly interpretation. Geophysics , 1972, 37(3): 507-517. DOI:10.1190/1.1440276 |
[2] | Nabighian M N. Additional comments on the analytic signal of two-dimensional magnetic bodies with polygonal cross-section. Geophysics , 1974, 39(1): 85-92. DOI:10.1190/1.1440416 |
[3] | Roest W R, Verhoef J, Pikington M. Magnetic interpretation using the 3-D analytic signal. Geophysics , 1992, 57(1): 116-125. DOI:10.1190/1.1443174 |
[4] | Hsu S K, Sibuet J C, Shyu C T. High-resolution detection of geologic boundaries from potential-field anomalies: An enhanced analytic signal technique. Geophysics , 1996, 61(2): 373-386. DOI:10.1190/1.1443966 |
[5] | Hsu S K, Dorothee C, Shyu C Ti. Depth to magnetic source using analytic signal. Geophysics , 1998, 63(6): 1947-1957. DOI:10.1190/1.1444488 |
[6] | 管志宁, 姚长利. 倾斜板体磁异常总梯度模反演方法. 地球科学 , 1997, 22(1): 81–85. Guan Z N, Yao C L. Inversion of the total gradient modulus of magnetic anomaly due to dipping dike. Earth Science-Journal of China University of Geosciences (in Chinese) (in Chinese) , 1997, 22(1): 81-85. |
[7] | Thompson D T. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data. Geophysics , 1982, 47(1): 31-37. DOI:10.1190/1.1441278 |
[8] | Reid A N, Allsop J M, Granser H, et al. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics , 1990, 55(1): 80-91. DOI:10.1190/1.1442774 |
[9] | Barbosa V C F, Silva J B C, Medeiros W E. Stability analysis and improvement of structural index estimation in Euler deconvolution. Geophysics , 1999, 64: 48-60. DOI:10.1190/1.1444529 |
[10] | 姚长利, 管志宁, 吴其斌, 等. 欧拉反演方法分析及实用技术改进. 物探与化探 , 2004, 28(2): 150–155. Yao C L, Guan Z N, Wu Q B, et al. An analysis of Euler deconvolution and its improvement. Geophysical and Geochemical Exploration (in Chinese) (in Chinese) , 2004, 28(2): 150-155. |
[11] | Salem A, Ravat D. A combined analytic signal and Euler method (AN-EUL) for automatic interpretation of magnetic data. Geophysics , 2003, 68: 1952-1961. DOI:10.1190/1.1635049 |
[12] | Salem A, Ravat D, Gamey T J, et al. Analytic signal approach and its applicability in environmental magnetic investigations. Journal of Applied Geophysics , 2002, 49: 231-244. DOI:10.1016/S0926-9851(02)00125-8 |
[13] | Salem A, Williams S, Fairheard D, et al. interpretation of magnetic data using tilt-angle derivatives. Geophysics , 2008, 71(1): 1-10. |
[14] | 董树文. 长江中下游铁铜矿带成因之构造分析. 中国地质科学院院报 , 1991, 23: 43–56. Dong S W. Tectonic analysis on genesis of metallogenetic belt of middle-lower Yangtze river. Bulletin of the Chinese Academy of Geological Sciences (in Chinese) (in Chinese) , 1991, 23: 43-56. |
[15] | 董树文, 吴宣志, 高锐, 等. 大别造山带地壳速度结构与动力学. 地球物理学报 , 1998, 41(3): 349–361. Dong S W, Wu X Z, Gao R, et al. On the crust velocity levels and dynamics of the Dabieshan orogenic belt. Chinese J| Geophys| (Acta Geophysica Sinica) (in Chinese) (in Chinese) , 1998, 41(3): 349-361. |
[16] | 董树文, 项怀顺, 高锐, 等. 长江中下游庐江-枞阳火山岩矿集区深部结构与成矿作用. 岩石学报 , 2010, 26(9): 2529–2542. Dong S W, Xiang H S, Gao R, et al. Deep structure and ore formation within Lujiang-Zongyang volcanic ore concentrated area in Middle to Lower reaches of Yangtze river. Acta Petrologica Sinica (in Chinese) (in Chinese) , 2010, 26(9): 2529-2542. |
[17] | 戚良刚, 王建伟, 黄传杨, 等. 庐枞地区重磁场特征及其地质认识. 安徽地质 , 2008, 18(4): 277–281. Qi L G, Wang J W, Huang C Y, et al. Characteristic and geological understanding of gravity, aeromagnetic field in Luzong area. Geology of Anhui (in Chinese) (in Chinese) , 2008, 18(4): 277-281. |
[18] | 张季生, 高锐, 李秋生, 等. 庐枞火山岩盆地及其外围重、磁场特征. 岩石学报 , 2010, 26(9): 2613–2612. Zhang J S, Gao R, Li Q S, et al. Characteristics of gravity and magnetic field of Luzong volcano basin and its periphery. Acta Petrologica Sinica (in Chinese) (in Chinese) , 2010, 26(9): 2613-2612. |
[19] | Keating P B. Weighted Euler deconvolution of gravity data. Geophysics , 1998, 63: 1595-1603. DOI:10.1190/1.1444456 |