地球物理学报  2011, Vol. 54 Issue (6): 1634-1641   PDF    
欧拉反褶积与解析信号相结合的位场反演方法
张季生1,2 , 高锐1,2 , 李秋生1,2 , 管烨1,2 , 彭聪3 , 王海燕1,2     
1. 中国地质科学院地质研究所, 北京 100037;
2. 中国地质科学院深部探测与地球动力学开放实验室, 北京 100037;
3. 中国地质科学院矿产资源研究所, 北京 100037
摘要: 由于解析信号具有不受(二维)或少受磁化方向影响, 能够较好反映磁性体边界的特性, 因此受到人们的重视.欧拉反褶积法可以确定场源的位置和深度以及形状因子, 具有较强的适应性.因此前人提出将二者相结合的方法.针对前人提出的方法中存在受高频干扰严重的问题, 本文提出低阶的欧拉反褶积与解析信号相结合的位场反演方法.本方法在反演中只需计算磁异常的一、二阶导数, 这样能够将高频噪声的干扰减少到人们可以控制的水平.文中详细推导了反演计算公式.模型计算证明了方法的正确性, 同时探讨了噪声和邻近磁性体干扰影响的问题.对安徽庐江-枞阳火山岩盆地磁异常反演计算的结果证实了方法的实用性.本方法也可以应用于重力异常的反演计算中.
关键词: 磁异常      重力异常      快速反演      欧拉反褶积      庐江-枞阳盆地     
A combined Euler and analytic signal method for an inversion calculation of potential data
ZHANG Ji-Sheng1,2, GAO Rui1,2, LI Qiu-Sheng1,2, GUAN Ye1,2, PENG Cong3, WANG Hai-Yan1,2     
1. Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China;
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
Abstract: The advantage of using the analytic signal method to determine source boundary parameters is the independence from the magnetization direction (2D) or is only slightly dependent on the magnetization direction. The position, depth and structural index of source can be obtained by the EULDPH method. The combined Euler and analytic signal method with high order derivative has been presented. But one has to calculate the third derivative of potential field anomaly. The method, therefore, is sensitive to noise in the data. Especially, in the case of complicated feature of magnetic anomaly, the noise distorts the higher order analytic signal so badly that it might spoil computed results. In order to circumvent this difficulty we develop a combined Euler and analytic signal method for an inversion calculation of potential data, with which the order of derivative calculated is not more than two. Synthetic data and real data demonstrate the effectiveness of this improved method. In addition, this method can be used for inversion of gravity anomaly..
Key words: Magnetic anomaly      Gravity anomaly      Automatic inversion      EULDPH      Lujiang-Zongyang basin     
1 引言

目前人们已经研究开发了多种基于利用磁异常导数的快速反演方法,用这些方法可以确定场源的参数,例如边界的位置和埋藏深度.解析信号法(AnalyticSignalMethod, 也称总梯度模法)由于解析信号具有不受(二维)或少受磁化方向影响,能够较好地反映磁性体边界,利用解析信号模确定磁源深度和产状,将其应用于地质构造“填图"已受到人们的重视[1-6].欧拉反褶积法[7-10]以Euler齐次方程为理论基础,利用总场及其梯度,依据欧拉齐次方程组确定场源的位置和深度,并可以得出形状因子(StructuralIndex, 又译为构造指数)对地质体进行识别,方法具有较强的适应性.但由于欧拉反褶积法采用窗口滑动计算,可以得到一系列深度点,在众多的深度点中如何分辨和确定有效的深度点一直是困扰人们的一个难题.

Salem 等(2003,2002)将Euler方程与解析信号表达式联立,推导出ANEUL 反演方法[1112],该方法可以同时计算磁性体的几何形状、水平位置和埋藏深度.但是用该方法需要计算磁异常的三阶导数,我们知道,计算导数相当于高通滤波.计算三阶导数不可避免地引进高频噪声干扰,而这种高频噪声干扰在实际资料反演计算中难以压制.为了改进上述方法,Salem 等(2008)又提出Tiltangle 方法[13],该方法虽然只需要计算一、二阶导数,但是需要确定磁性体走向方向,并沿垂直磁性体走向方向进行插值,然后通过求解超定线性方程组,获得场源的相应参数.这种方法利用多个测点上的场值进行反演计算,较利用单点反演压制随机噪声能力强,但是需要逐个确定磁性体走向,并要求解超定线性方程组,计算量较大.

本文提出同时计算磁性体形状和埋深等参数的低阶欧拉反褶积与解析信号相结合的反演方法.该方法在反演中只需计算磁异常的一、二阶导数,这样将高频噪声的干扰减少到人们可以控制的水平.

2 原理

Reid等[8]给出Euler齐次方程如下:

(1)

其中,T为位于点(x0y0z0)的磁性体在点(x,y,z)处产生的磁异常,n为形状因子,B为正常场存在的偏差.

计算式(1)对x的偏导数,经整理得

(2)

x= x0y= y0z=0,埋深z0=h,方程(2)可以化简为

(3)

同理,有

(4)

(5)

hTzz(0)= (n+1)Tz(0),(5)对(3),(4),(5)两边平方,加后开方,得

(6)

其中,

向上延拓高度Δz时,即令x= x0y= y0z=-Δ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 水平圆柱磁性体反演计算结果(a)水平圆柱磁性引起的磁异常,水平圆柱体模型绘在图的最下方.中心埋深为3km, 磁化倾角为30°,磁偏角为0°;(b)计算的解析信号曲线;(c)计算的形状因子曲线;(d)计算的埋藏深度曲线. Fig. 1 Total magnetic field response of the horizontalc ylinder shown at the bottom of the figure (a) The depth to the center of the cylinder isequal to 3 km.The inclination is 30°; (b) The analytic signal calculated;(c)The structural index calculated;(d) The depth calculated.
图 2 无限延深薄板磁性体反演计算结果(a)无限延深薄板磁性体引起的磁异常,无限延深薄板模型绘在图的最下方.顶面埋深为1km, 倾角为50°,磁化倾角为70°;(b)计算的解析信号曲线;(c)计算的形状因子曲线;(d)计算的埋藏深度曲线. Fig. 2 Total magnetic field response of the sheet shown at the bottom of the figure (a) The depth to the top of the sheet is equal to 1 km.The dipping is 50°.The inclination 70°;(b) The analytic signalc alculated;(c) The structural index calculated;(d) The depth calculated.
图 3 无限延深台阶磁性体反演计算结果(a)无限延深台阶磁性体引起的磁异常,无限延深台阶模型绘 在图的最下方.顶面埋深为1km,倾角为50°,磁化倾角为70°;(b)计算的解析信号曲线;(c)计算的形状因子曲线; (d)计算的埋藏深度曲线. Fig. 3 Total magnetic field response of the step shown at the bottom of the figure (a) The depth to the top of the step is equal to 1 km.The dipping is 50°.The inclination 70°;(b) The analytic signal calculated;(c)The structural index calculated; (d) The depth calculated.
表 1 二度磁性体模型的计算结果 Table 1 Calculation result of two dimensional magnetic models

图 1~图 3可以看出解析信号的极大值位于磁性体在地面的投影位置.离该极大值点越远,计算形状因子和埋深的误差越大.这是由于本方法计算的假设前提是解析信号的极大值位于磁性体中心的上方.

3.2 随机噪声对计算精度的影响

像所有的梯度定量解释方法一样,随机噪声对欧拉反褶积与解析信号相结合方法的计算精度有一定的影响.由于随机噪声位于高频段,因此可以用向上延拓的方法减少噪声的影响.本文对上面讨论的水平圆柱体的理论异常加上一组幅值为磁异常的2%~8%的随机噪声.通过向上延拓不同高度,延拓高度由0.2km 增加到1.8km, 计算相应的形状因子和埋藏深度(见图 4图 5).由图 4图 5 可以看出,延拓高度较小时,计算的误差大.随着延拓高度的增大,计算的误差减小,当延拓超过一定高度时(在本例中,大于1.4km)计算结果是准确的.

图 4 对水平圆柱体(见图 1)加上一组随机噪声,向上延拓不同高度后的形状因子计算结果. Fig. 4 Calculation result of the structural index for the cylinder (seen Fig.1)after adding different random noise and continuing upward
图 5 对水平圆柱体(见图 1)加上一组随机噪声,向上延拓不同高度后的埋藏深度计算结果 Fig. 5 Calculation result of the depth for the cylinder (seen Fig.1)after adding random noise and continuing upward
3.3 邻近磁性体干扰对计算精度的影响

为了研究邻近磁性体干扰的影响,我们计算了相距一定水平距离的垂直磁化无限延深薄板和水平圆柱体模型的埋深和形状因子(见图 6).用本文提出的方法计算了不同距离间隔的埋深h和形状因子n,并计算了其相对误差.为了便于分析,绘制了磁性体的几何参数埋深h和形状因子n的计算误差与间隔L与埋深h的比值之间的关系曲线(见图 7).从图 7可以看出,相邻磁性体的干扰所产生的计算相对误差具有随它们之间的间隔增大而减小的趋势.当间隔与深度的比值大于11时,几何参数的计算误差小于4%.需要说明的是,几何参数的计算误差与磁性体的相对大小规模及形状有关.Hsu等[4]曾用所提出的算法研究垂直磁化无限延深的台阶和厚板叠加模型的参数计算误差关系证实了这一点.

图 6 垂直磁化水平叠加的无限延伸薄板和水平圆柱体引起的磁异常图中L表示无限延伸薄板和水平圆柱体模型之间的间隔. Fig. 6 Magnetic anomaly caused by the sheet and cylinder magnetic models magnetized vertically L indicates the gap distance.
图 7 磁性体的埋深h和形状因子SI的计算误差与间隔L与埋深h的比值之间的关系 Fig. 7 Relationship between the calculation error of depth h and SI and the ratio of gap distance L to depth h
3.4 直立长方体磁性模型的计算

与上述模型不同,对于球体,只有在垂直磁化的情况下,解析信号的极大值位于磁性体中心的上方[11].而长方体等非球体,理论上不能严格符合欧拉方程条件.斜磁化长方体磁模型解析信号的形态不具对称性[6].实际中非理想模型多见,但是在一定的条件下是可以用本文提出的方法进行反演计算的.本文选用一系列埋深不同的垂直磁化直立长方体磁性模型,用过模型体中心剖面的计算结果进行分析和讨论.直立长方体磁性模型的顶面的宽度为20km, 长度为20km.延深为50km, 埋藏深度为1~5km, 间隔为0.5km.图 8 为埋藏深度为2km的直立长方体的磁异常平面图.图 9 为过直立长方体(图 8)的中心剖面磁异常和解析信号图.由图 9可以看出,在直立长方体边缘上方出现类似于无限延深厚板的解析信号极大值.我们用前面提出的方法计算磁性体的参数.

图 8 垂直磁化直立长方体的磁异常图 图中,灰色正方形表示直立长方体磁性模型的顶面在地面的投影,其宽度为20km, 长度为20km.模型的顶面埋深为2km, 延深为50km.磁异常的单位为nT. Fig. 8 Magnetic anomaly map of the cuboids magnetized vertically The gray square indicates the projection of the top plane on the surface, which width is 20 km and length 20 km.Depth of the top plane is 2 km and depth deepened 50 km.Unit of magnetic anomaly is nT.
图 9 过垂直磁化直立长方体(见图 8)中心的剖面磁异常和解析信号图 Fig. 9 Magnetic anomaly and amplitude of analytic signal profile across the center of cuboids magnetized vertically (seen Fig.8)

不同深度垂直磁化直立长方体的计算结果(见图 10图 11)表明,随着埋深增大,计算结果的误差增加.当埋深不大于2km 时,埋深计算误差小于10.5%,解析信号极大值点相对于直立长方体边缘位移了一倍点距的距离,形状因子为0.028~0.013,近于无限延深台阶的特征.

图 10 埋深不同的垂直磁化直立长方体的埋藏深度计算误差 Fig. 10 Calculation error of the depth of the cuboids magnetized vertically with the various depth
图 11 埋深不同的垂直磁化直立长方体的解析信号极大值点相对于直立长方体边缘的位移误差 Fig. 11 Calculation error of the horizontal shift of the edges of the cuboids magnetized vertically with the various depth
4 实例

庐江-枞阳陆相火山岩盆地位于长江中下游前陆区,即大别-苏鲁超高压碰撞造山带前陆缩短带.地处郯-庐断裂带的南段[14-16].庐枞盆地为一“继承式"或“基底坳陷型"火山盆地,呈耳状,长轴约56km, 短轴约24km, 面积约1032km2.庐枞盆地磁异常大致呈NE 向展布,北东宽,南西窄(见图 12).前人对该地区岩石磁性的研究表明:区内上元古界新生界岩层磁性很弱,属微磁或无磁性岩石.中碱性侵入岩及潜火山岩具中等到强磁性.庐枞盆地磁异常主要由这些中碱性岩类引起[1718].

图 12 庐江-枞阳地区磁异常(ΔT)图 Fig. 12 Magnetic anomaly(ΔT)map of Lujiang-Zongyang area

由于斜磁化球体的解析信号极大值点偏离球心在地面的投影位置[11],因此我们首先对庐枞盆地的磁异常进行化极计算(见图 13).由于庐枞盆地的磁异常深部为规模大的正长岩基[18],因此我们对化极磁异常进行磁场分离.用向上延拓5km 的磁场作为区域场,用化极磁异常减去区域场得到局部磁异常(见图 14).先对分离后的中碱性侵入岩体及潜火山岩引起的局部磁异常向上延拓1km 和1.5km;再计算上延后磁异常的解析信号和一阶增强解析信号的模,搜索出它们的极大值;最后按公式(11),(12)计算各极大值点处的形状因子和埋深.反演计算结果(见图 15)表明:庐枞盆地内部磁性体的形状因子大小范围为0.35~2.91,即中碱性侵入岩体及潜火山岩的形状近似为无限延深的板状体和水平圆柱体、球体.板状体顶面埋深为0.3~1.2km.水平圆柱体和球体的中心埋深为1.6~2.9km.

图 13 庐江-枞阳地区化极磁异常(ΔT)图 Fig. 13 Magnetic anomaly(ΔT)reduced to pole of Lujiang-Zongyang area
图 14 庐江-枞阳地区局部化极磁异常(ΔT)图 Fig. 14 Magnetic anomaly(ΔT)reduced to pole and continued 5 km upward of Lujiang-Zongyang area
图 15 庐江-枞阳地区磁性体深度图图中背景为解析信号模图,蓝色圆表示磁性体中心埋深,绿色圆表示磁性体顶面埋深. Fig. 15 Depth of magnetic bodies in the Lujiang-Zongyang area Background is the modulus of analytic signal, blue circle represents the central depth of magnetic bodies, green circle represents the topdepth of magnetic bodies.
5 结论与讨论

与传统的位场反演方法不同,本文提出的欧拉反褶积与解析信号相结合的位场反演方法可以同时计算磁性体的形状因子和埋深.即用该方法可以确定磁性体的形状,并判断磁性体在地下的产出状态.这个方法在计算中仅仅需要计算磁异常的一、二阶导数,这样将高频噪声的干扰减少到人们可以控制的水平.随机噪声干扰的实验证实了这一点.相邻磁性体的干扰所产生的计算误差随它们之间的间隔增大而减小,干扰误差水平与磁性体的相对大小规模及形状有关.实际中非理想模型是多见的,对于垂直磁化直立长方体磁模型的计算结果表明,在一定的条件下是可以用本文提出的方法进行反演计算的.对安徽庐枞火山岩盆地磁异常反演计算的结果证实了该方法的实用性.

本方法也可以应用于重力异常的反演计算中,计算时需要注意形状因子的不同.根据泊松公式,磁异常与重力异常存在一阶导数的关系,因此前者的形状因子比后者的大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