地球物理学报  2017, Vol. 60 Issue (12): 4840-4854   PDF    
非均匀物性条件下多尺度窗口修正法换算磁源重力异常及在寻找DSO的应用
颜廷杰1,2, 甘文龙3, 管彦武1     
1. 吉林大学 地球探测科学与技术学院, 长春 130026;
2. 中国地质调查局发展研究中心, 北京 100037;
3. 加拿大拉贝克世纪铁矿公司, 多伦多 M4Y1M7, 加拿大
摘要:DSO(Direct Shipping Ore,直运块矿)是一种高品位的富铁矿.在加拿大拉布拉多(Labrador)地槽谢菲尔威利(Schefferville)铁矿成矿带,含铁建造苏克曼(Sokoman)组地层全铁含量低、具高密度、强磁性,能够引起高重力异常与高磁异常;而风化淋滤后富集的赤铁矿和针铁矿等(也称DSO)全铁含量高,具高密度、无磁性,仅能够引起高重力异常.采用一般的滤波方法不能提取DSO的重、磁异常.本文采用基于泊松(Poisson)公式的磁场换算磁源重力异常(pseudo-gravity anomaly)方法,由磁场换算磁源重力异常,再与实测重力异常对比,得到纯粹由高密度、无磁性的DSO产生的剩余重力异常,对剩余重力异常采用密度成像与2.5D反演方法解释DSO.泊松公式虽然提出时间很长,但迄今为止仅仅用在资料解释中的定性分析,本文推导并实现了密度磁性非均匀条件下经典泊松公式的形式与实现过程,提出了多尺度窗口滑动线性回归修正的磁场换算磁源重力异常方法,使该公式的数学原理能够对重、磁异常的反演解释定量化.最后本文将多尺度窗口滑动线性回归修正的换算磁源重力异常方法用于加拿大拉布拉多地槽谢菲尔威利铁矿成矿带铁矿勘探,较好地解决了寻找高密度、无磁性DSO的问题.
关键词: 非均匀物性      多尺度窗口滑动线性回归修正      磁源重力异常      DSO      谢菲尔威利铁矿成矿带     
Magnetic fields reduction to pseudo-gravity anomalies using multi-scale windows in the condition of non-uniform property and its application to searching for DSO
YAN Ting-Jie1,2, GAN Wen-Long3, GUAN Yan-Wu1     
1. College of Geo-exploration Science and Technology, Jilin University, Changchun 130026, China;
2. Development and Research Center of China Geological Survey, Beijing 100037, China;
3. Century Iron Mines Corporation, Toronto M4Y1M7, Canada
Abstract: The Direct-Shipping Ore (DSO) is a kind of high-grade iron-rich ore. The Sokoman iron formation with low iron content is of high density and strong magnetism, causing high gravity and magnetic anomalies in the Schefferville mining district, central Labrador Trough in Canada. While the leached enrichment of hematite and goethite (also called DSO) with high iron content is featured by high density and weak magnetism, causing only high gravity anomalies. Conventional filter algorithms cannot extract the magnetic and gravity anomalies arising from DSO. Based on the Poisson equation, we propose a magnetic reduction to the pseudo-gravity anomaly. It transforms magnetic fields into gravity anomalies, then compares the results with gravity observations. The residual gravity anomalies arising from only high-density and non-magnetic DSO are derived. The residual anomalies are interpreted by density inversion and 2.5D interactive inversion. Although the Poisson formula was proposed for a long time, it is used only in qualitative analysis of gravity and magnetic data so far. This paper derives classical Poisson formula in the condition of non-uniform density and magnetic properties, then realizes its process, and presents magnetic field reduction to the gravity field using multi-scale linear regression analysis by moving-windows. The principles of mathematics in this formula can quantify the magnetic and gravity inversion and interpretation. Finally, the pseudo-gravity reduction using multi-scale linear regression analysis by moving-windows is adopted to successfully delineate the DSO with high-density and low magnetism from the Sokoman iron formation in the Schefferville mining district, central Labrador Trough in Canada.
Key words: Non-uniform property    Multi-scale linear regression analysis by moving-windows    Pseudo-gravity anomaly    Direct Shipping Ore (DSO)    Schefferville iron-ore mining district    
1 引言

在加拿大拉布拉多(Labrador)地槽谢菲尔威利(Schefferville)铁矿成矿带(图 1),含铁建造苏克曼(Sokoman)组地层受白垩纪热带雨林气候风化,形成由表生弱变质燧石含铁建造组成的软铁矿石,主要包括经淋滤后的富含赤铁矿和针铁矿的碳酸岩和硅质矿物岩石,其全铁含量可达59%,是高品位的富铁矿,也称为DSO(Direct Shipping Ore,直运块矿),而含铁建造苏克曼组地层含铁燧石岩,全铁低,是一种低品位的贫铁矿(图 2).如何在含铁建造苏克曼组地层找DSO是地球物理勘探方法要解决的问题.

图 1 拉布拉多地槽谢菲尔威利铁矿成矿带1) Fig. 1 Metallogenic belt in Labrador and Schefferville, Canada1)
图 2 地层剖面示意图1) UMH-上部致密块状赤铁矿层;RC-红色硅质岩;LMH-下部致密块状赤铁矿层;LRC-下红色硅质岩;RS-露丝页岩. Fig. 2 Schematic diagram of strata1) UMH-upper high density hematite. RC-red silica rock. LMH-lower high density hematite. LRC-lower red silica rock; RS-rose shale.

1) Claude Duplessis P.Eng.2013. NI 43-101 Technical Report Joyce Lake DSO Iron Project New found Land & Labrador.

位场分离与识别常用向上延拓(Fuller,1967)、匹配滤波(Spector and Grant, 1970)与解析信号振幅(Nabighian,1972)等方法.近年来采用小波分析(Fedu and Quarta, 1998)、边界识别(Miller and Singh, 1994)、最小曲率(Mickus and Aiken, 1991, 纪晓琳等,2015)、经验模态分析等方法(Huang et al., 1998, 曾琴琴等,2010),但是对于识别和提取DSO所产生的重、磁异常有效的方法则不多.

1971年,Seguin(1971)在研究拉布拉多地槽的重、磁异常中认为,层位较低的板岩和石英岩引起区域性的负异常,没有蚀变的含铁岩层引起中等强度的异常,而DSO则引起强度大的局部高频异常,于是采用二次滤波的方法提取DSO产生的局部重力异常,未能获得好的地质效果.然而,由重、磁勘探理论可以知道,局部高频异常是场源埋深较浅的反映,采用二次滤波的方法提取的局部重力异常只能反映埋深较浅的DSO,近年来人们在乔伊斯湖和黑鸟湖等地区的勘察中发现了许多DSO都赋存于埋深相对较深的向斜轴部,所引起的异常幅值和频率往往不如位于向斜两翼的含铁岩层高.另外,拉布拉多地槽的DSO是位于苏克曼含铁岩层中,两者埋深往往相近,频谱特征相似.当不同地质体的埋深相近时,就很难用滤波的方法来分离.图 3a是乔伊斯湖1号勘探剖面,乔伊斯湖地区勘探程度高,该勘探剖面内高密度强磁性的舒克曼地层及DSO基本上已被钻孔所控制,根据图 3a地质模型,分别计算地层和DSO所产生的重力异常进行功率谱分析,图 3b是乔伊斯湖1号勘探剖面地层重力异常的功率谱和DSO重力异常的功率谱,表明富铁矿体(DSO)与苏克曼组地层的频谱都具有相似的特征,即高中低频都有分布,很难构造一个滤波器来从苏克曼组地层与富铁矿体(DSO)共同产生的重力异常中分离出富铁矿体(DSO)的局部重力高异常.

图 3 乔伊斯湖1号勘探剖面与重力异常功率谱曲线 (a)中深紫色填充的地质体为DSO;(b)中黑线为地层和DSO的功率谱,蓝线为地层的功率谱,红线为富铁矿体(或DSO)的功率谱. Fig. 3 No.1 profile and gravity anomaly power spectral curves at Joyce lake (a) The purple area is geologic body DSO; (b) Blue line is power spectrum of strata. Red line is power spectrum of rich iron mine (or DSO).

2015年,Joёl Simard在《Report on a ground gravity survey performed on the Sunny Lake project lac le fer property, Blackbird Lake area Labrador trough, Quebec》报告中采用磁化率成像的方法将处于低磁异常区域及磁化率相对较小的部位圈定为DSO的有利部位.然而,当地层构造为向斜模式(图 2),向斜轴部磁性地层下凹,无论是否存在DSO总是会产生低磁异常,另外,板岩和石英岩等磁化率也是相对较小的,因此,单独依靠磁测资料同样无法很好地解决DSO的反演问题.

目前,对于DSO重磁异常的综合解释人们还是停留在通过局部重力高和磁力低的异常组合特征来定性分析和解释的阶段(王炯辉和余宇星,2015).

Garland(1951)提出了采用泊松公式将重磁异常结合起来进行分析解释,并采用重磁异常确定单个地质体的磁化率和密度的比值,并以此来确定地质体的性质.基于泊松公式的重磁异常对应分析方法为进行DSO的解释提供了很好途径.Cordell等(1971)在北大西洋海山地区采用泊松公式确定单个地质体的物性参数以及剩磁的方向.Chandler等(1981)采用泊松公式分析了叠加场的重磁资料,从而将该方法推广到多源的重磁异常研究中.Barano(1957)提出了一种利用磁源重力异常来解释航磁异常的方法.Von Frese等(1982)用计算磁源重力异常方法研究了北美大陆的重磁异常性.Chandler(1991)应用泊松窗口移动的分析方法研究了明尼苏达中东部一佩尼奥克统造山带前寒武纪对应的地质体,取得了较好的效果.曾昭发等(2006)通过具体的实例分析了基于泊松定理来确定地质体总磁化方向及其在分析火山岩活动中的作用.

在国内,刘心铸(1985)率先从国外引入重磁对应分析方法,介绍了针对单个异常体、叠加异常体和组合异常体等理论模型的重磁对应分析及应用.此后,刘沈衡(1992)林长松等(1992)殷秀华等(1999)宋景明(2007)刘彦华等(2008)赵汝敏等(2013)吴兴等(2014)也把该方法应用于区域重磁资料构造解释中.

上述的关于利用泊松公式进行重磁异常对应分析的文献都是用于区域重磁资料定性解释,未见将该方法用于重磁异常的定量分离提取与反演解释.

加拿大拉布拉多地槽谢菲尔威利铁矿成矿带上含铁建造苏克曼组地层具有高密度、强磁性(表 1),能够产生高重力异常与高磁力异常;而DSO高密度、弱磁性或无磁性(表 1),仅能够产生高重力异常.由此可见,加拿大拉布拉多地槽谢菲尔威利铁矿成矿带上的磁异常主要是含铁建造苏克曼组地层的反映;重力高异常则为含铁建造苏克曼组地层与DSO的共同反映;单独依靠重力或磁法是无法可靠地识别和提取DSO异常并进行反演解释.本文利用含铁建造苏克曼组地层与DSO重磁场特征,采用基于泊松公式的磁场换算磁源重力场方法,从实测磁场换算磁源重力异常,再与实测重力异常对比得出剩余重力异常,该剩余重力异常即高密度、无磁性的富铁矿DSO产生,对剩余重力异常采用密度成像与2.5D反演解释DSO.该方法有效解决了加拿大拉布拉多地槽谢菲尔威利铁矿成矿带找富铁矿问题.

表 1 乔伊斯湖地区岩矿石密度和磁性参数2) Table 1 Density and magnetic values in the Joyce lake region2)

2) Joёl Simard.2015.Report on a ground gravity survey performed on the Sunny Lake project lac le fer property, Blackbird Lake area Labrador trough, Quebec.2)

2 方法原理 2.1 均匀物性条件下的磁源重力异常

物体v在任一点P(x, y, z)点产生的磁位V(P)和引力位U(P)分别可表示为

(1)

其中J(M)是物体的磁化强度,ρ(M)是物体的密度,k是万有引力常数,r为磁性体中M(ξ, η, ζ)点到观测点P(x, y, z)的矢径,其模值为

(2)

如果物体为均匀密度且均匀磁化,即J(M)和ρ(M)均为恒量,则式(1)变为:

(3)

由此可以得出对于同一密度均匀且磁化均匀的物体而言,在P点产生的磁位与引力位之间的关系式为

(4)

这就是说,对于均匀磁化且密度均匀的物体所产生的磁位,可以由该物体的质量所产生的引力位来求得,这就是Poisson公式.

于是由Poisson公式(4),可导出垂直磁化时磁位V与磁源重力异常关系为

(5)

垂直磁化时磁位V为垂直磁化时的Za沿着垂直磁化反方向的曲线积分,即:

(6)

于是(5)式中的磁源重力异常可表示为

(7)

由于垂直磁化时的Za等于垂直磁化时的ΔT,因此对于航磁异常ΔT而言,如果已知泊松比,将ΔT化到地磁极,然后利用(7)式就可以计算磁源重力异常.

由于实测重磁异常在整理过程中所选取的背景场值不同,往往转换的磁源重力异常与实测重力异常有一常数背景差,所以公式(7)后往往需要加一个常数C,即:

(8)

由(8)式可以看出当重磁异常同源时,重力异常与化极磁异常的垂向积分呈线性关系,因此可以通过线性回归法估算出泊松比,然后利用(8)式将磁异常换算成磁源重力异常.

目前,地面高精度磁测多为对总磁场强度绝对值的测量,磁源重力异常往往是通过ΔT来换算,因此,利用(8)式计算磁源重力异常时,首先要将ΔT进行频率域化极处理,获得垂直磁化状态下的Za磁异常.(8)式是空间域换算磁源重力异常的方法,计算效率低而且实现困难.目前,主要是利用频率域进行磁源重力异常的换算,ΔT换算磁源重力异常的频率域计算公式为

(9)

其中Sg′(u, v, 0)、ST(u, v, 0)分别为磁源重力异常Δg′和磁异常ΔT的频谱,qt0qt1分别为地磁场方向因子和磁化方向因子,即:

L0M0N0为地磁场方向余弦,L1M1N1为磁化方向的方向余弦.

2.2 非均匀物性条件下的磁源重力异常换算方法

对于非均匀磁化的物体则无法根据(1)式导出磁位与引力位的Poisson公式,进而利用Poisson公式来将磁异常换算磁源重力异常.

将非均匀磁化的物体离散化成许多小单元vi,并假设离散化后的每个单元vi都是密度均匀且磁化均匀,则根据(4)式有:

(10)

其中Vi(P)和Ui(P)分别为物体单元vi在P点产生的磁位和引力位,Jiρi则是物体单元vi的磁化强度和密度.根据(7)式有:

(11)

其中Zai、Δg′i分别为物体单元vi在垂直磁化情况下产生的磁异常和磁源重力异常,σi为物体每个单元vi的泊松比.

整个物体的磁源重力异常为

(12)

因此,对于磁性和密度非均匀的情况下,利用(12)式计算磁源重力异常时需要采用窗口滑动,通过线性回归法估算出泊松比σi.与均匀物性情况一样,非均匀物性情况下的磁源重力异常也是通过频率域方法来进行计算.

3 理论模型分析

根据苏克曼地层含铁建造的向斜模式(图 2),本文建立了简单的向斜理论模型,分析和探讨采用基于泊松公式的磁场换算磁源重力异常方法提取DSO异常的有效性和精度.

图 4a表示苏克曼地层中未被风化淋滤的铁燧岩向斜模型(无DSO),铁燧岩的磁化强度为10 A·m-1,为了使模型分析简便假设垂直磁化(实际斜磁化情况下可通过化极转化为垂直磁化),剩余密度为0.8×103 kg·m-3(围岩密度按2.6×103 kg·m-3).如图 4b表示苏克曼地层中铁燧岩在向斜核部被风化淋滤形成DSO的模型,DSO的垂直磁化强度为0.5 A·m-1,剩余密度为0.5×103 kg·m-3.

图 4 换算磁源重力异常方法的理论模型检验 (a)无DSO的苏克曼地层含铁建造向斜模型;(b)有DSO的苏克曼地层含铁建造向斜模型;(c)无DSO的磁源重力异常与剩余重力异常;(d)有DSO的磁源重力异常与剩余重力异常. Fig. 4 Theoretical model verification based on the method of converting magnetic-source gravity anomalies (a) Syncline model without DSO in the Sokoman iron-bearing formation; (b) Syncline model with DSO in the Sokoman iron-bearing formation; (c) Magnetic-source gravity anomalies and residual gravity anomalies without DSO; (d) Magnetic-source gravity anomalies and residual gravity anomalies with DSO.

当泊松比取14913319 A·kg/N(苏克曼地层的泊松比)时,用图 4a模型产生的磁异常计算磁源重力异常,然后将其从模型产生的重力异常中减去,获得剩余重力异常,从剩余重力异常的结果可以看出,当无DSO时剩余重力异常为接近0的直线.同样方法,利用图 4b模型产生的重磁异常来提取剩余重力异常,如图 4d所示,提取的剩余重力异常与理论模型DSO引起的重力异常十分接近,由此可以看出,利用本文提出的方法可以有效的提取DSO产生的重力异常,并以此来进行DSO的反演与解释.

提取DSO异常的效果好坏取决于磁场换算磁源重力异常的精度,下面本文将探讨剖面长度、测点点距、密度和磁化强度对精度的影响.

3.1 影响换算磁源重力异常精度的因素

(1) 剖面长度对换算磁源重力异常的影响

图 4a中苏克曼组地层向斜模型为例,向斜模型宽1000 m,苏克曼组地层剩余密度为0.59×103 kg·m-3,磁化强度为10 A·m-1.图 5ad是点距固定为25 m和剖面长度分别为1600 m、3200 m、6400 m、12800 m时苏克曼组地层产生的理论重力异常与由苏克曼组地层产生的磁异常换算磁源重力异常的对比曲线图,由对比结果可以看出,随着剖面长度增大,换算的磁源重力异常的误差逐渐减小.当剖面长度与模型宽度的比值为小于3倍时,换算磁源重力异常的误差较大,而当剖面长度与模型宽度的比值大于3倍时,误差迅速变小(见图 5e).

图 5 剖面长度对换算磁源重力异常的影响 (a)剖面长度1600 m; (b)剖面长度3200 m; (c)剖面长度6400 m; (d)剖面长度12800 m; (e)相对误差曲线.图中红线为磁异常换算的磁源重力异常,蓝线为理论模型正演重力异常. Fig. 5 Effect of profile length on conversion of magnetic-source gravity anomalies (a) The profile length is 1600 m; (b) The profile length is 3200 m; (c) The profile length is 6400 m; (d) The profile length is 12800 m; (e) The curve of relative error. Red lines are gravity anomalies converted from magnetic anomalies. Blue lines are gravity anomalies from forward modeling.

由于在加拿大拉布拉多地槽的勘探区施工条件困难,重力剖面较短,因此在解释时要注意剖面长度的影响,适当地将短剖面外延.

(2) 点距对换算磁源重力异常的影响

同样,以图 4a中的苏克曼组地层向斜模型为例,向斜模型宽1000 m,苏克曼组地层剩余密度为0.59×103 kg·m-3,磁化强度为10 A·m-1.图 6ad是剖面长度固定为6400 m和点距分别为25 m、50 m、100 m、200 m时苏克曼组地层产生的理论重力异常与由苏克曼组地层产生的磁异常换算磁源重力异常的对比曲线图,由对比结果可以看出,随着剖面点距增大,换算的磁源重力异常的误差逐渐增大.当点距与模型宽度的比值为小于1/5时,误差较小,而当点距与模型宽度的比值大于1/5时误差就迅速增大(见图 6e).究其原因,误差增大是由于换算磁源重力异常是在频率域进行的,点距太大不满足采样定理所致.因此在施工中应尽量采用样较小的点距.

图 6 点距对换算磁源重力异常的影响 (a)点距25 m; (b)点距50 m; (c)点距100 m; (d)点距200 m; (e)相对误差.图中红线为磁异常换算的磁源重力异常,蓝线为理论模型正演重力异常. Fig. 6 Effect of station distance on the conversion of magnetic-source gravity anomalies (a) The distance of stations is 25 m; (b) The distance of stations is 50 m; (c) The distance of stations is 100 m; (d) The distance of stations is 200 m; (e) The curve of relative error. Red lines are gravity anomalies converted from magnetic anomalies. Blue lines are gravity anomalies from forward modeling.

(3) 密度参数选取不准对换算磁源重力异常的影响

同样,以图 4a中的苏克曼组地层向斜模型为例,向斜模型宽1000 m,苏克曼组地层剩余密度为0.59×103 kg·m-3,磁化强度为10 A·m-1.图 7为剩余密度取0.1×103~1.0×103 kg·m-3、磁化强度固定取10 A·m-1时换算的磁源重力异常的相对误差曲线,由结果可以看出,剩余密度相对于真实值(0.59×103 kg·m-3)无论是取偏大还是取偏小都会使磁源重力异常误差增大,剩余密度偏差越多误差越大.因此,在资料处理解释时要注意密度变化的影响,在野外注意地层岩性与密度变化,在处理中采取必要的措施提高换算磁源重力异常的精度.

图 7 密度参数不准对换算磁源重力异常的影响 (磁化强度固定取为10 A·m-1) Fig. 7 Effect of inaccurate density value on the conversion of magnetic-source gravity anomalies (The magnetization intensity is fixed at 10 A·m-1)

(4) 磁化强度参数不准对换算磁源重力异常的影响

同样,还是以图 4a中的苏克曼组地层向斜模型为例,向斜模型宽1000 m,苏克曼组地层剩余密度为0.59×103 kg·m-3,磁化强度为10 A·m-1.图 8为磁化强度取3~20 A·m-1、剩余密度固定取0.59×103 kg·m-3时换算的磁源重力异常的相对误差曲线,由结果可以看出,磁化强度相对于真实值(10 A·m-1)无论是取偏大还是取偏小都会使磁源重力异常误差增大,磁化强度偏差越多误差越大.与对密度变化影响的措施一样,在野外注意地层岩性与磁性变化,在处理中采取相应的措施.

图 8 磁化强度参数不准对换算磁源重力异常的影响 (密度固定取为0.59×103 kg·m-3) Fig. 8 Effect of inaccurate magnetization intensity value on the conversion of magnetic-source gravity anomalies (The density is fixed at 0.59×103 kg·m-3)
3.2 泊松比的估算

由上述关于密度和磁化强度对磁源重力异常的影响分析可以看出,密度和磁化强度的选取影响了泊松比的精度,从而影响磁源重力异常的计算精度,因此,要获得高精度的磁源重力异常,首先要估算出合适的泊松比.由磁源重力异常的计算公式(8)可以知道,重磁异常同源、泊松比选取准确时,实测重力异常与换算的磁源重力异常相等,如果根据公式(8)取泊松比取为1 A·kg/N时换算的磁源重力异常记为Δg1,那么实测重力异常与Δg1不相等,但成线性关系(图 9所示),因此,可以利用实测重力异常与Δg1通过线性回归方法获得实际的泊松比,从而得到正确的磁源重力异常Δg′,具体做法是:

图 9 泊松比为1 A·kg/N时的磁源重力异常Δg1与重力异常Δg的线性关系 Fig. 9 Linear relationship between magnetic-source gravity anomalies Δg1 and gravity anomalies Δg

(1) 利用磁异常计算泊松比为1 A·kg/N时的磁源重力异常Δg1.

(2) 做实测重力异常Δg与泊松比为1 A·kg/N时的磁源重力异常Δg1的散点图,通过线性回归求取线性回归的斜率,该斜率即为泊松比.

(3) 将该泊松比乘以泊松比为1 A·kg/N时的磁源重力异常Δg1即为磁异常换算的磁源重力异常Δg′.

3.3 密度与磁性横向不均匀情况的窗口滑动线性回归修正法

从加拿大拉布拉多地槽的实际勘探剖面及钻孔磁化率、密度测井的结果可以看出,苏克曼组地层(铁燧岩)和DSO的磁性和密度是随磁铁矿含量、风化、淋滤作用程度的不同而有变化,如果在提取剩余异常时与上述单一均匀模型一样采用统一的泊松比计算磁源重力异常,势必得不到正确的结果.为了说明泊松比的选取对结果的影响,本文建立了由三个埋深相同但具有不同剩余密度与磁化强度的直立薄板模型进行理论计算和分析,三个直立薄板的磁化强度分别为10 A·m-1、0 A·m-1、3 A·m-1,剩余密度分别为0.8×103 kg·m-3、0.5×103 kg·m-3、0.6×103 kg·m-3,泊松比分别为14913319 A·kg/N、0 A·kg/N、5965327 A·kg/N,并假设第二个直立薄板为DSO.当泊松比统一取第一个直立薄板的泊松比时磁源重力异常换算“不足”(图 10),当泊松比统一取第三个直立薄板的泊松比时磁源重力换算“过量”(图 11),而且无论泊松比取多少,统一的泊松比提取的剩余异常都无法同时消除第一个和第三个直立薄板异常,并精确地提取第二个直立薄板的异常.

图 10 泊松比取14913319 A·kg/N时提取的剩余异常 Fig. 10 Residual anomalies when Poisson's ratio is 14913319 A·kg/N
图 11 泊松比取5965327 A·kg/N时提取的剩余异常 Fig. 11 Residual anomalies when Poisson's ratio is 5965327 A·kg/N

根据(12)式,对于磁性和密度横向非均匀的情况,首先需要采用窗口滑动,通过线性回归法估算出泊松比σi的横向变化值,然后对泊松比取1时的磁源重力异常进行修正.图 12为三个直立薄板异常模型窗口滑动线性回归分析法提取的剩余异常(窗口大小取3500 m,相当于单个异常宽度的1~2倍),从图上可以看出基本消除了第一个和第三个直立薄板异常,并较好地提取了第二个直立薄板.

3.4 密度与磁性横向和纵向不均匀情况的多尺度窗口滑动线性回归修正法

同样从实际勘探剖面及钻孔结果可以看出,苏克曼组地层(铁燧岩)和DSO的磁性和密度在纵向上也是有变化的,苏克曼组地层与DSO除了有上述横向叠加情况,还有垂向叠加的情况.当苏克曼组地层与DSO是垂向叠加时,由于埋藏深度不同,产生的异常范围不同,采用固定窗口的滑动线性回归修正法来估算泊松比,进而计算磁源重力异常,势必得不到理想的结果.由重磁勘探理论可以知道,埋藏浅的地质体产生的异常范围小,必须采用小窗口进行分析,而埋藏深的地质体产生的异常范围大,必须采用大窗口进行分析.为此,论文利用小波变换进行多尺度窗口滑动线性回归修正法对磁源重力异常进行修正,具体做法为:

(1) 将重力异常和磁异常分别进行小波多尺度分解.

(2) 由于尺度相同的重力异常和磁异常所反映的场源深度基本相同,因此,相同尺度a的重力异常和磁异常根据异常的范围选择合适的窗口进行滑动线性回归修正法计算该尺度下的磁源重力异常Δg′a.

(3) 将所有尺度下计算的磁源重力异常Δg′a叠加得到总的磁源重力异常Δg′=∑Δg′a.

(4) 实测重力异常减去总的磁源重力异常得到剩余重力异常.

图 13建立了由三个埋深不同且具有不同剩余密度与磁化强度的直立板状体模型进行理论计算和分析,三个直立板状体的中心埋深分别25 m、100 m、400 m,磁化强度分别为3 A·m-1、0 A·m-1、10 A·m-1,剩余密度分别为0.6×103 kg·m-3、0.5×103 kg·m-3、0.8×103 kg·m-3,泊松比分别为5965327 A·kg/N、0 A·kg/N、14913319 A·kg/N,并假设第二个直立薄板为DSO.首先将重力异常和磁异常进行6个尺度的分解,可获得6个细节和1个逼近,阶次越低的小波细节频率越高,异常范围越小,进行窗口滑动线性回归分析时采用的窗口就越小,本文通过分析各尺度异常的波长,以波长的2~3倍作为窗口大小.图 13模型的各阶重磁异常所采用的窗口分别为750 m、1000 m、1300 m、1800 m、2500 m、3500 m、5000 m.采用窗口滑动线性回归分析法提取的剩余异常曲线如图 13的橙色点划线所示,对比理论DSO异常曲线(图 13中绿色实线)可以看出,窗口滑动线性回归分析法无法很好地同时消除第一个和第三个直立薄板异常;采用多尺度窗口滑动线性回归分析法提取的剩余异常曲线如图 13的红色虚线所示,可以看出,多尺度窗口滑动线性回归分析法能较好地消除第一个和第三个直立薄板异常,提取第二个直立薄板的异常.

图 12 窗口滑动线性回归分析法提取的剩余异常(窗口大小3500 m) Fig. 12 Extracted residual anomalies based on the linear regression analysis with sliding windows (window size is 3500 m)
图 13 多尺度窗口滑动线性回归分析法提取的剩余异常 Fig. 13 Extracted residual anomalies based on the linear regression analysis with multiscale windows
4 应用实例

黑鸟湖的地面重力2014年完成,采用加拿大产的CG5重力仪,点距25 m,线距100 m,测线长度约2400 m.地面高精度磁测2014年完成,采用加拿大产的GSM-19T磁力仪,点距5 m,线距50 m,测线长度约2000 m.航磁2010年完成,飞机飞行高度20 m,ΔT磁异常经过IGRF正常场改正后获得.

把多尺度窗口滑动线性回归分析法应用于黑鸟湖矿区重磁资料的解释,以L12剖面为例,具体做法为:

(1) 由于地面磁测剖面长度较短,为了避免由此带来的误差,勘探剖面的主体磁异常为地面磁测异常,把航磁异常通过向下延拓20 m换算到地面与拼接到地面磁异常上来扩展地面磁异常的剖面长度.

(2) 采用小波多尺度分解方法提取苏克曼地层及矿体产生的局部重力异常和局部磁异常,消除背景场对磁源重力异常的影响.

(3) 利用多尺度窗口滑动线性回归分析法把局部磁异常换算为“磁源重力异常”,该“磁源重力异常”是苏克曼组地层产生的.

(4) 再把实测的局部重力异常减去“磁源重力异常”,剩余的局部重力异常是富铁矿体(或DSO)产生的.

(5) 由于黑鸟湖地区的钻孔资料较少,很难直接进行2.5D/3D人机交互的精细反演解释,因此,采用3D密度成像与3D磁化强度成像(Li and Oldenburg, 2000Boulanger and Chouteau, 2001Peter and Lelievre, 2006)的反演结果作为参考建立2.5D/3D人机交互反演的初始模型,然后结合研究区的地质构造及已有的少量钻孔进行约束,利用人机交互精细反演方法(Rasmussen, 1979Wallace, 2007Yang and Li, 2011)反演富铁矿体(DSO).

图 14a中,黑线是实测局部重力异常,红色虚线是局部磁异常,红色点线是根据局部磁异常换算的磁源重力异常.蓝色虚线是实测局部重力异常减去磁源重力异常所得的剩余异常,可以看出在横坐标1200 m、1400~1800 m两处有明显的高剩余重力与低磁异常的组合特征,推测它们是由高密度无磁性的富铁矿体(DSO)产生的.由于目前3D密度和3D磁化强度成像反演的结果很难满足精细解释的要求(如图 14b图 14c所示),因此,本文参照图 14g地表出露地层与向斜特征、图 14b图 14c三维密度成像与磁化强度成像结果,以此作为2.5D/3D交互反演建模的参考.图 14d是2.5D/3D人机交互反演方法解释结果,红色地质体为DSO.经2005年验证,钻孔BB-14-001、BB-14-005、BB-14-006、BB-14-007、BB-14-010均钻遇较厚的DSO.

图 14 加拿大拉布拉多地槽谢菲尔威利铁矿成矿带L12线反演DSO (a) L12线重磁异常曲线;(b)密度成像结果;(c)磁化强度成像;(d) 2.5D/3D人机交互反演结果;(e)局部重力异常;(f)局部磁异常;(g)地质图. Fig. 14 DSO inversion for line L12 in the Schefferville mining district, central Labrador Trough in Canada (a) The profile of gravity and magnetic anomaly for Line 12; (b) The result of density imaging; (c) The result of susceptibility imaging; (d) The result of 2.5D/3D interactive inversion; (e) Local gravity anomaly; (f) Local magnetic anomaly; (g) Geology map.
5 结论

(1) 采用基于泊松公式的磁场换算磁源重力场方法,从高密度强磁性的贫矿与高密度、无磁性的富铁矿共同产生重磁异常中分离出仅由富铁矿产生的剩余重力异常,再由该剩余重力异常解释富铁矿,较好地解决了在贫矿中找富铁矿的问题.

(2) 磁场换算磁源重力异常的精度与剖面长度、点距有关,也受密度、磁化强度参数取值是否正确影响,而且还受选取的背景场及滑动窗口大小的影响.因此,在实际工作中应该使剖面足够长,将短剖面适当外延,并在施工中应尽量采用样较小的点距,密度和磁化强度的选择要以物性测量结果作为最重要的参考依据.为了减小背景场的影响,在计算磁源重力异常之前应该根据工区实际地质情况将区域背景场与局部场进行分离,利用局部场来换算磁源重力异常.滑动窗口的可以选择异常波长2~3倍,也可通过变差函数等方法获取.

(3) 本文采用密度与磁性横向不均匀情况的窗口滑动线性回归修正与密度与磁性纵向不均匀情况的多尺度窗口滑动线性回归修正有效地解决密度与磁性不均匀情况的磁源重力异常换算.

(4) 本文假设DSO为高密度无磁性,但实际中某些DSO存在一定的磁性,因此会造成提取的剩余异常存在误差.DSO如果有磁性,一般也比较弱,但密度高,其泊松比很小,因此为了减小由此造成的误差,在实际中可根据实际工区的岩矿石物性,选取合适的阈值,当估算的泊松比小于该阈值时则认为是DSO,置其泊松比为0,然后再进行磁源重力异常修正处理.

参考文献
Baranov V A. 1957. A new method for interpretation of aeromagnetic maps:pseudo-gravimetric anomalies. Geophysics, 22(2): 359-383. DOI:10.1190/1.1438369
Boulanger O, Chouteau M. 2001. Constraints in 3D gravity inversion. Geophysics Prospecting, 49(2): 265-280. DOI:10.1046/j.1365-2478.2001.00254.x
Chandler V W. 1991. Moving-window Poisson analysis of gravity and magnetic data from the Penokean orogen, east-central Minnesota. Geophysics, 56(1): 123-132. DOI:10.1190/1.1442948
Chandler V W, Koski J S, Hinze W J, et al. 1981. Analysis of multisource gravity and magnetic data sets by moving-window application of Poisson's theorem. Geophysics, 46(1): 30-39. DOI:10.1190/1.1441136
Cordell L, Taylor P T. 1971. Investigation of magnetization and density of a North Atlantic seamount using Poisson's theorem. Geophysics, 36(5): 919-937. DOI:10.1190/1.1440224
Fedu M, Quarta T. 1998. Wavelet analysis for the regional-residual and local separation of potential field anomalies. Geophysical Prospecting, 46(5): 507-525. DOI:10.1046/j.1365-2478.1998.00105.x
Fuller B D. 1967. Two-dimensional frequency analysis and design of grid operators. Mining Geophysics, 42(2): 658-708.
Garland G D. 1951. Combined analysis of gravity and magnetic anomalies. Geophysics, 16(1): 51-62. DOI:10.1190/1.1437650
Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society A Mathematical, Phtsical and Engineering Sciences, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193
Ji X L, Wang W Y, Qiu Z Y. 2015. The research to the minimum curvature technique for potential field data separation. Chinese Journal of Geophysics (in Chinese), 58(3): 1042-1058. DOI:10.6038/cjg20150329
Li Y G, Oldenburg D W. 2000. Joint inversion of surface and three-component borehole magnetic data. Geophysics, 65(2): 540-552. DOI:10.1190/1.1444749
Lin C S, Guan Z N, Wu C J. 1992. Comparing pseudo-gravity anomaly with gravity anomaly in the East China See and the study on deep geological structure. Acta Oceanologica Sinca (in Chinese), 14(5): 75-85.
Liu S H. 1992. Application of correspondence Analysis to Interpretation of Gravity and Magnetic Anomaly in Zalong Area, Songliao Basin. Contributions To Geology and Mineral Resources Research (in Chinese), 7(3): 97-105.
Liu X T. 1985. The method of gravity and magnetic anomaly analysis, one of the introduction of gravity and magnetic mothed learned from America. Foreign Geoexploration Technology(3): 1-8.
Liu Y H, Chen Z G, Ouyang C L. 2008. The application of correspondence analysis of gravity and magnetic anomalies in Xiangshan area. Geophysical and Geochemical Exploration (in Chinese), 32(6): 586-589.
Mickus K L, Aiken C L V, Kennedy W D. 1991. Regional residual gravity anomaly separation using the minimum-curvature technique. Geophysics, 56(2): 279-283. DOI:10.1190/1.1443041
Miller H G, Singh V. 1994. Potential field tilt-a new concept for location of potential field sources. Journal of Applied Geophysics, 32(2): 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(3): 507-517. DOI:10.1190/1.1440276
Peter G, Lelièvre P G, Oldenburg D W. 2006. Magnetic forward modelling and inversion for high susceptibility. Geophysical Journal International, 166(1): 76-90. DOI:10.1111/j.1365-246X.2006.02964.x
Rasmussen R, Pedersen L B. 1979. End corrections in potential field modeling. Geophysical Prospecting, 27(4): 749-760. DOI:10.1111/j.1365-2478.1979.tb00994.x
Seguin M K. 1971. Discovery of direct-shipping iron ore by geophysical methods in the central part of the Labrador Trough. Geophysical Prospecting, 19(3): 459-486. DOI:10.1111/j.1365-2478.1971.tb00611.x
Song J M. 2007. The gravity anomaly separate technology based on correspondence analysis mothed and its application. Natural Gas Industry (in Chinese), 27(S1): 318-319.
Spector A, Grant F S. 1970. Statistical models for interpreting aeromagnetic data. Geophysics, 35(2): 293-302. DOI:10.1190/1.1440092
Von Frese R R B, Hinze W J, Braile L W. 1982. Regional North American gravity and magnetic anomaly correlations. Geophysical Journal International, 69(3): 745-761. DOI:10.1111/gji.1982.69.issue-3
Wallace Y. 2007. 3D Modelling of Banded Iron Formation incorporating demagnetization-a case study at the Musselwhite Mine, Ontario, Canada. Exploration Geophysics, 38(4): 254-259. DOI:10.1071/EG07027
Wang J H, Yu Y X. 2015. Ore-forming processes and exploration methods of DSO iron deposits in Labrador Trough, Canada. Earth Science Frontiers (in Chinese), 22(2): 187-199.
Wu X, Liu C C, Li J. 2014. Pseudo-gravity technology and its application. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 36(4): 410-414.
Yang Y S, Li Y Y, Liu T Y. 2011. Interactive 3D forward modeling of total field surface and three-component borehole magnetic data for the Daye iron-ore deposit (Central China). Journal of Applied Geophysics, 75(2): 254-263. DOI:10.1016/j.jappgeo.2011.07.010
Yin X H, Ni Y S, Liu Z P, et al. 1999. Internal correspondence analysis of gravimetric and magnetic data for the northern part of south-north Tectonic zone. Seismology and Geology (in Chinese), 21(4): 370-376.
Zeng Q Q, Liu T Y. 2010. A potential field separation method based on Empirical Mode Decomposition. Oil Geophysical Prospecting (in Chinese), 45(6): 914-918.
Zeng Z F, Wu Y G, Hao L B, et al. 2006. The Poisson's theorem based analysis method and application of magnetic and gravity anomalies. Journal of Jilin University (Earth Science Edition) (in Chinese), 36(2): 279-283.
Zhao R M, Zhu G H, Bai G J. 2013. The application of correspondence analysis technique of gravitational and magnetic anomalies in the process of factual data, in Myanmar region. Progress in Geophysics (in Chinese), 28(2): 914-919. DOI:10.6038/pg20130244
纪晓琳, 王万银, 邱之云. 2015. 最小曲率位场分离方法研究. 地球物理学报, 58(3): 1042–1058. DOI:10.6038/cjg20150329
林长松, 管志宁, 吴朝钧. 1992. 东海磁源重力异常、重力异常的对比和深部地质构造研究. 海洋学报, 14(5): 75–85.
刘沈衡. 1992. 重磁异常对应分析应用及评述. 地质找矿论丛, 7(3): 97–105.
刘心铸. 1985. 重磁异常对应分析:赴美学习重磁方法介绍之一. 国外地质勘探技术(3): 1–8.
刘彦华, 陈宗刚, 欧阳长亮. 2008. 重磁异常对应分析在相山地区的应用. 物探与化探, 32(6): 586–589.
宋景明. 2007. 基于对应分析的重力异常分离技术及应用. 天然气工业, 27(S1): 318–319.
王炯辉, 余宇星. 2015. 加拿大拉布拉多地槽DSO型铁矿床成矿规律及找矿方法研究. 地学前缘, 22(2): 187–199.
吴兴, 刘春成, 李军. 2014. 磁源重力技术及其应用. 物探化探计算技术, 36(4): 410–414.
殷秀华, 黎益仕, 刘占坡, 等. 1999. 南北构造带北段重磁异常的对应分析. 地震地质, 21(4): 370–376.
曾琴琴, 刘天佑. 2010. 一种基于经验模态分解(EMD)的位场分离方法. 石油地球物理勘探, 45(6): 914–918.
曾昭发, 吴燕冈, 郝立波, 等. 2006. 基于泊松定理的重磁异常分析方法及应用. 吉林大学学报(地球科学版), 36(2): 279–283.
赵汝敏, 朱光辉, 柏冠军. 2013. 重磁对应分析技术在缅甸X区块的应用研究. 地球物理学进展, 28(2): 914–919. DOI:10.6038/pg20130244