2017, Vol. 60
2. 合肥工业大学矿床成因与勘查技术研究中心, 合肥 230009;
3. 安徽省地质调查院, 合肥 230001
2. Ore Deposit and Exploration Centre (ODEC), Hefei University of Technology, Hefei 230009, China;
3. Geological Survey of Anhui Province, Hefei 230001, China
基于位场观测数据的边界识别技术已广泛应用于地球物理勘探工作中，运用该方法可以有效识别隐伏断层和地质体的边界，对于揭露深部地质结构、圈定隐伏岩体和预测深部找矿靶区有着积极的意义.
依据原始重磁异常难以划分地质体的分布范围，但重磁异常水平导数的极大值、垂直导数的零值与地质体的边界相对应，垂直导数(Evjen, 1936)是最早被用来描述密度体或磁性体边界的方法，早期的边界识别方法(Cordell, 1979; Nabighian, 1984; Blakely and Simpson, 1986; Roest et al., 1992; Hsu et al., 1996)均不能同时显示出浅部和深部异常体的边界，这是由于异常的导数随深度的衰减速度较快，无法清晰地识别出较深地质体的边界.为了解决这一问题，同时清晰地识别出深、浅部的地质体边界，相关研究开始致力于均衡边界识别滤波器.第一个均衡边界识别滤波器倾斜角法(Miller and Singh, 1994)能很好地均衡不同强度异常的振幅，但不能清晰地识别地质体的边界；Rajagopalan和Milligan (1994)使用自动增益控制技术均衡不同大小的振幅，该方法的结果依赖于窗口大小的选择；基于tilt angle的总水平导数THDR方法(Verduzco et al., 2004)受噪声影响较大，也不能很好地给出较深地质体的边界；利用解析信号振幅(Roest et al., 1992)均衡总水平导数的Theta图(Wijns et al., 2005)取得了很好的效果，但该方法识别的边界存在一定的发散；Cooper和Cowan (2006)提出了改进的倾斜角TDX方法进行异常体边界的识别；Cooper在2009年提出采用解析信号的希尔伯特变换进行地质体边界的识别，该方法能清晰地显示磁源体的边界；马国庆等(2012, 2013)和Ma (2015)提出增强型均衡滤波器和增强型局部相位滤波器，利用不同阶水平导数之间的组合进行边界识别；李丽丽等(2014)使用归一化总水平导数法获得异常体的水平位置和深度.近年来基于重力梯度张量(GGT)数据的特征值进行边界识别也屡见报道(Sertcelik and Kafadar, 2012；Zhou et al., 2013; Yuan et al., 2014；Wang et al., 2015；Yuan and Yu, 2015; 袁园等，2015)，此外，基于小波分析的多尺度边缘检测也得到应用(Hornby et al., 1999; 严加永等，2015).以上工作使得均衡滤波器的研究不断深入，更清晰的识别深、浅部地质体的边界.
本文基于Theta图法定义了新的边界识别滤波器，该滤波器显著压制了Theta图法识别深部地质体边界时的发散现象，较好地平衡了深部和浅部边界，能清晰地识别出地质体的边界.本文设计了两组模型，其中一组模型同时包含了正、负异常并加入高斯噪声，来检验方法的稳定性和边界识别效果.文中以长江中下游成矿带的庐枞(庐江-枞阳)矿集区实测重磁数据为例，进行了实例分析.
2 基于Theta法改进的均衡滤波器 2.1 Theta图滤波器Roest等(1992)证明了解析信号的最大振幅可以用来识别地质体的边界.解析信号的表达式为
|
(1) |
式中f为位场的总强度，
从式(1)可以得到解析信号的最大振幅为
|
(2) |
Wijns等(2005)基于解析信号振幅提出了Theta图法，该方法的表达式为
|
(3) |
即：
|
(4) |
Theta图方法对深部地质体边界的识别十分有效，但深、浅部地质体的边界在平面图上的展现效果正好相反，浅部的地质体边界清晰而收敛，深部的地质体边界清晰但显著发散.本文针对Theta图法的不足，进行了改进，首先对位场的总强度f求取垂向一阶偏导数fz，在垂向一阶偏导数的基础上求取新的Theta值，本文定义为Theta1，其表达式为
|
(5) |
式(5)中的垂向二阶偏导数fzz使得Theta1的计算结果细节十分丰富，通过实际应用我们发现过多的信息可能会造成地质上的解释困难，考虑到fzx、fzy和fzz是相对独立的物理量，fzx、fzy用于识别x、y水平方向的边界，对上式稍作改动，用垂向一阶导数fz代替高阶的fzz，并除以采样间距的倍数，使物理量单位一致，定义为Theta2，其表达式为
|
(6) |
式(6)中，p为无量纲的量，经测试p取值0.05~5效果较好，较大的p值会降低均衡效果，较小的p值能使边界更收敛，h为网格化的数据采样间距，如本文使用的庐枞矿集区1：5万重磁数据网格化后的间距为0.5 km.
2.3 Theta图滤波器的主要计算公式在波数域中kx、ky分别是x、y方向的波数，波数k的表达式为
|
(7) |
则其垂向一阶偏导数为
|
(8) |
fzx、fzy和fzz的表达式分别为：
|
(9) |
|
(10) |
|
(11) |
由上述一系列计算公式可知，分别求取kx、ky和k，根据偏导数的需求将观测平面位场f作Fourier变换(使用FFT计算)并乘以相应波数，再作Fourier逆变换(使用IFFT计算)(刘东甲等, 2009, 2012).
为了避免边部第一类间断点引起的吉布斯效应及边部计算点数据的损失，需要对原始数据做扩边处理，本文使用了余弦扩边法；为消除频率域中高频干扰，使用最小二乘二十五点圆滑.
除了波数域计算程序，本文还分别编写了相应的空间域计算程序，并考虑到fzz在空间域直接使用高阶导数公式计算可能会增加噪声干扰，使用了水平导数和进行变换，公式为
|
(12) |
空间域编程计算相对繁杂，计算耗时更长一些；而波数域编程计算程序相对简洁，执行效率更高些.通过模型数据对波数域和空间域计算结果进行了对比，Theta1、Theta2方法以及本文使用的其他方法两者较为接近，异常体边界都能准确、清楚的识别，相对而言空间域计算结果噪声要低些.
2.4 本文用于对比的其他滤波器为了比较本文滤波器的边界识别效果，我们分别选取了Theta图(Wijns et al., 2005)、THDR方法(Verduzco et al., 2004)、TDX方法(Cooper and Cowan, 2006)和基于三维梯度张量的滤波器(Yuan and Yu, 2015)进行对比.
2.4.1 THDR和TDX方法Miller和Singh (1994)定义了倾斜角法(tilt angle)滤波器，表达式为
|
(13) |
Verduzco等(2004)建议对T值再计算一次水平导数，定义为THDR，表达式为
|
(14) |
Cooper和Cowan (2006)建议对T值的分子、分母进行对调，定义为TDX，表达式为
|
(15) |
Zhou等(2013)和Yuan和Yu (2015)分别对位场数据的三维梯度张量矩阵进行了研究.公式为
|
(16) |
据Yuan和Yu (2015)的研究成果，本文选取效果较好的N_ED_A进行对比，其表达式为
|
(17) |
为了检验边界识别方法的效果，建立了如图 1所示的一组模型.假设在100 km×100 km的平面范围内分布有四个长方体场源，其编号依次为1、2、3、4，各异常体的埋深逐步加深，顶、底面深度依次为0.5~1 km、1.5~2 km、2.5~3 km和3.5~4 km.
|
图 1 理论模型的平面图和3D图 (a)平面图；(b) 3D图. Fig. 1 Planar and 3D views of synthetic model (a) Planar view; (b) 3D view. |
以重力模型为例，本文考虑到两种情况：(1)模型中长方体均为正的剩余密度，各长方体剩余密度均设为0.3 g·cm-3，图 2a为由此模型合成的重力异常.(2)模型中同时含有正、负剩余密度的长方体，设1号、3号长方体的剩余密度为0.3 g·cm-3，2号、4号长方体的剩余密度为-0.3 g·cm-3，并且在重力异常中叠加异常峰值3%的正态分布随机噪声，图 3a为由此模型合成的重力异常.两个模型的采样间距h均取0.5 km.
|
图 2 不同边界识别滤波器的效果对比 (a)合成的重力异常；(b) Theta图法；(c) THDR方法；(d) TDX方法；(e)基于三维梯度张量的N_ED_A方法；(f) Theta1方法；(g) Theta2方法，p取2. Fig. 2 Comparison of results using different filters (a) Synthetic gravity anomalies; (b) Theta map method; (c) THDR method; (d) TDX method; (e) N_ED_A of method based on 3D gradient tensor; (f) Theta1 method; (g) Theta2 method with p of 2. |
|
图 3 不同边界识别滤波器的效果对比 (a)合成的重力异常；(b) Theta图法；(c) THDR方法；(d) TDX方法；(e)基于三维梯度张量的N_ED_A方法；(f) Theta1方法；(g) Theta2方法，p取2. Fig. 3 Same as Fig.2 but for gravity data of synthetic model added by 3% Gaussian noise |
对模型中四个长方体的剩余密度均为正的情况，以合成的理论重力异常数据为基础，分别使用Theta图、THDR、TDX和基于三维梯度张量的N_ED_A方法，以及本文的Theta1和Theta2方法进行边界识别，图 2a为模型合成的重力异常，图 2b-g为不同方法的边界识别结果，对计算结果进行了归一化处理.
图 2b是Theta图法边界识别结果，可以看出随着1~4号长方体的深度逐渐加深，Theta图法识别的深部异常体的边界也变得更加宽大.图 2c为THDR方法的识别结果，由图可见该方法对浅部的边界识别清晰，随着长方体的深度逐渐加深，识别的边界信号越来越模糊.图 2dTDX方法识别的边界比较均衡，也比较清晰，特别是识别的深部边界比Theta法要收敛.图 2e是基于三维梯度张量的N_ED_A方法的边界识别结果，该方法较好的均衡了深、浅部异常，但由于运用了高阶偏导数fzz，产生了噪声干扰.图 2f为本文改进的Theta1方法的结果，相对于Theta图法，Theta1方法识别的边界更加均衡，深部边界更收敛，由于也运用了高阶偏导数fzz，图 2f同图 2e产生了类似的噪声干扰，虽然这两种方法产生了噪声干扰，但更清晰均衡的显示了边界位置.针对Theta1方法使用高阶偏导数fzz产生的噪声干扰，Theta2方法进行了改进在垂向上用一阶导数fz代替高阶导数fzz，图 2g是Theta2方法的边界识别结果，由图 2b-g不同方法的结果对比可知，Theta2方法的边界识别效果相对最好，深、浅部边界均衡，边界清晰且更为收敛，不产生干扰.
为了测试算法的稳定性.对模型中四个长方体的剩余密度包含正、负，且叠加了异常峰值3%高斯噪声的情况，以合成的理论重力异常数据为基础，分别使用前文所述的多种方法进行边界识别，图 3a为模型合成的重力异常，图 3b-g为不同方法的边界识别结果，对计算结果进行了归一化处理.
由图 3边界识别结果可以看出，由于正、负异常同时存在的原因，Theta图法、THDR方法、TDX方法和Theta2方法，都产生了额外的虚假边界，不同的是THDR方法识别的长方体边界明显被噪声所掩盖，而Theta图法、TDX方法和Theta2方法识别的边界在噪声背景下，仍然能清晰的识别出模型的边界，Theta法识别的边界最为宽大，TDX法识别的边界比Theta法要收敛些，本文提出的Theta2方法对模型周边的噪声压制范围最大，在低背景噪声下显著的突出模型的边界位置，并且识别的边界比Theta法和TDX方法更为收敛.基于三维梯度张量的N_ED_A方法和Theta1方法运用了更高阶的垂向偏导数，未产生额外的边界，但整体噪声水平偏高，这两种方法在放大噪声干扰的同时，也突出显示出了边界位置，相对而言，Theta1方法比N_ED_A方法能在放大噪声干扰的同时更加均衡突出的识别出模型边界.
4 庐枞矿集区重、磁边界识别 4.1 地质背景庐枞矿集区位于长江中下游断陷带内，扬子板块北缘，西邻秦岭-大别造山带，是长江中下游成矿带七个重要的矿集区之一(常印佛等，1991；翟裕生等，1992；唐永成等，1998；Pan and Dong, 1999；袁峰等，2008；周涛发等, 2008, 2010；董树文等，2010；高锐等，2010；吕庆田等，2014).矿集区主体为庐枞火山岩盆地，盆地总体呈NE向展布，长轴方向长约60 km，面积约1032 km2，形似“耳状”(任启江等，1993).
庐枞矿集区内出露的地层主要有志留系、泥盆系、石炭系、二叠系、三叠系、侏罗系、白垩系及第四系.寒武-奥陶纪碳酸盐岩及碎屑岩主要出露于庐枞矿集区北部盛桥-东顾山地区，志留纪至中三叠世地层主要出露于庐枞火山岩盆地周边.盆地内部地表主要出露早白垩世陆相火山岩，盆地的直接基底为中三叠统-中侏罗统地层.庐枞盆地的火山岩由老至新划分为龙门院旋回、砖桥旋回、双庙旋回和浮山旋回.四个旋回火山岩的形成时代分别为：龙门院旋回134.8±1.8 Ma；砖桥旋回134.1±1.6 Ma；双庙旋回130.5±0.8 Ma；浮山旋回127.1±1.2 Ma，均为早白垩世岩浆活动的产物(周涛发等，2008).
|
图 4 庐枞矿集区地质矿产略图(据周涛发等，2010) Fig. 4 Geological map of Luzong ore district (after Zhou et al., 2010) |
庐枞矿集区内侵入岩分布广泛，侵入岩体规模大小悬殊，计有各类侵入体40余个.侵入岩岩石类型为闪长岩、二长岩、正长岩和A型花岗岩.闪长岩体主要出露在庐枞盆地的北西侧沙溪地区，产出有沙溪斑岩型铜矿；A型花岗岩体主要分布在庐枞盆地的东南侧，形成NE走向的岩浆岩带；二长岩体主要沿盆地中心分布；正长岩体在多个矿床的深部都有发现，如龙桥、罗河、泥河、小岭等铁矿，这些正长岩多以岩盖或岩基产出，侵位于中三叠世东马鞍山组，或早白垩世的龙门院组和砖桥组.侵入岩可划分成早、晚两期，早期侵入岩主要为二长岩和闪长岩类，以沙溪岩体、焦冲岩体、龙桥岩体和拔茅山岩体等为代表，成岩时代为134-130 Ma，主要分布在庐枞盆地北部；晚期侵入岩分为两类，第一类主要为正长岩类，以巴坛岩体、大缸窑岩体、罗岭岩体等为代表，成岩时代为129-123 Ma，主要分布在庐枞矿集区南部；第二类主要为A型花岗岩，以枞阳岩体和黄梅尖岩体等为代表，其形成时代约为126-123 Ma (周涛发等，2010；范裕等，2014).
庐枞火山岩盆地受控于4条边界断裂，即盆地西侧的罗河-缺口断裂，北侧的庐江-黄姑闸-铜陵拆离断层，东侧的长江逆冲断裂和陶家湾-施家湾断裂(任启江等，1993；吕庆田等，2014).
4.2 物性分析物性是连接地球物理与地质的纽带和桥梁.据以往庐枞矿集区的物性工作(张季生等，2010)，特别是深部探测技术与实验研究专项SinoProbe-03项目在庐枞矿集区实测了大量岩石、矿石的密度和磁化率(刘彦等，2012)并开展了三维岩性填图和地质建模实验(郭冬等，2014；祁光等，2014；严加永等，2014).庐枞矿集区的岩石物性研究显示(图 5)，主要岩性之间密度和磁化强度差异明显，火山沉积岩、正长岩和陆相红层沉积一般具有较低的密度(2.4~2.5 g·cm-3)，二长岩、闪长岩、次火山岩和砂岩、灰岩密度较高(2.6~2.7 g·cm-3).祁光等(2014)对庐枞矿集区基于先验密度信息约束的三维地质建模工作对物性进行了进一步细分，如白垩系浮山组(K1f)粗面岩的密度范围在2.40~2.51 g·cm-3，白垩系双庙组-龙门院组(K1sh-K1l)粗安岩、碎屑岩的密度范围在2.56~2.66 g·cm-3.火山岩及次火山岩的磁性与周边老地层及红层沉积差异更加明显.
|
图 5 庐枞矿集区岩(矿)石磁化率与密度对应关系图(据严加永等，2014) Fig. 5 Cross plot of rock′s density and susceptibility relationship in Luzong ore district (after Yan et al., 2014) |
庐枞矿集区东西结构由“两坳一隆”组成，西侧为潜山-孔城坳陷，东侧为庐枞火山岩盆地，二者之间以一基底隆起相隔，基底隆起带大致沿义津-罗河-缺口分布，由古生界-中生界地层组成.两个坳陷的物性特征完全不同，潜山-孔城坳陷无磁、低密度、低阻，浅部低速，反映典型的沉积盆地特征，与陆相碎屑沉积的物性特征吻合；庐枞火山岩盆地总体上表现为强磁、高阻和浅部高速特征，重力异常高低分明，界线清晰，反映侵入岩体的近垂直接触关系，这种物性特征与火山岩、侵入岩区相吻合.庐枞矿集区的“地垒”构造主要出现在矿集区的北部和东部，北部以庐江-黄姑闸-铜陵拆离断层(LHTD)为界大面积出露古生界地层，东部紧邻火山岩盆地为沿江隆起带，由逆冲的古生界-中生界地层组成，产生了局部高重力异常(吕庆田等，2014).
4.3 重、磁资料边界识别结果本文使用的重、磁资料来自SinoProbe-03项目和安徽庐枞地区隐伏及深部矿资源潜力调查项目，为1：5万经网格化的重、磁数据(点线距为500 m×500 m).
4.3.1 磁力数据边界识别庐枞矿集区位于中低纬度，需要通过化极将斜磁化转换为垂直磁化，据国际标准地磁参考场IGRF12模型，矿集区中心的地磁倾角为46.29°，磁偏角为-4.58°.如图 6a所示，庐枞矿集区航磁化极异常总体呈NE向展布，北东宽、南西窄，呈“蝌蚪”状(刘彦等，2012).沿异常走向长约95 km，北东宽约27 km，南西宽约18 km，火山岩盆地孙家坂、将军庙附近为主要的异常区，一些“孤立”异常主要沿火山岩盆地周边分布，重要异常有罗河、龙桥、焦冲、钱家铺、杨家市异常等，这些孤立异常基本对应了已发现的玢岩型铁硫矿床.盆地周边还分布有若干独立异常带或异常.盆地西北部沙溪地区存在航磁异常带，呈NNE向带状展布，强度达2000 nT，对应了沙溪斑岩型铜矿床.盆地北部边缘塘串河处分布一孤立航磁异常，同时存在局部重力异常，大致呈椭圆状，强度达1400 nT.盆地东部边缘施家湾存在一弱磁异常.
|
图 6 庐枞矿集区磁力异常边界识别效果对比 (a) 1：5万化极后磁力异常图；(b) Theta图法；(c) THDR方法；(d) TDX方法；(e)基于三维梯度张量的N_ED_A方法；(f) Theta1方法；(g) Theta2方法，p取1.8. Fig. 6 Comparison of the results using different filters on magnetic data in Luzong ore district (a) After reduction to pole magnetic anomalies; (b) Theta map method; (c) THDR method; (d) TDX method; (e) N_ED_A method based on 3D gradient tensor; (f) Theta1 method; (g) Theta2 method with p of 1.8. |
庐枞火山岩盆地-1 km以浅分布大量的矿点，在进行边界识别前，首先将航磁数据进行向上延拓，消除浅表干扰，对比了向上延拓0.5、1、2、3 km和5 km的数据结果，向上延拓1 km后的数据边界识别结果既能消除干扰，也保留了地质体的异常信息.由图 6b-g的边界识别结果可知：THDR方法边界比较模糊；Theta图法、TDX方法和Theta2方法的边界与实测资料显示的边界比较吻合，Theta图法的识别边界相对最为宽大；TDX方法相对于Theta图法，边界识别效果显著改善，Theta图中盆地边缘宽大的边界在TDX方法识别结果中明显变得更加细腻；Theta2方法的边界识别结果相对Theta图法和TDX方法，又进一步得到改善，不仅识别的主要地质体的边界更加细腻、收敛，而且在空间上计算区域的整体背景值低，在细节方面，Theta图法和TDX方法中识别的一些较小的呈半封闭的边界，在Theta2方法中能够识别出完整的呈封闭形态的边界.基于三维梯度张量的N_ED_A方法和Theta1方法结果近似一致，通过与实测资料对比可以看出这两种方法都识别了主要地质体边界，同时产生了丰富的细节边界，Theta1方法识别的边界更清晰些.
4.3.2 重力数据边界识别庐枞矿集区布格重力异常如图 7a所示，总体呈北东向展布，西北和东南侧边缘重力低，中部整体为重力高值背景区，在高背景重力场上又呈现为正、负高低异常相间.盆地中心的孙家坂、杨家市一带为重力高，盆地边缘的沙溪、塘串河、杨家桥等地均呈现重力高异常.结合矿床的分布位置，已知的矿床主要分布在重力高值梯度带上，如罗河铁矿、沙溪铜矿、龙桥铁矿等，少量分布在重力低值背景区的相对高值梯度带上，如泥河铁矿.
|
图 7 庐枞矿集区重力异常边界识别效果对比 (a) 1：5万布格重力异常图；(b) Theta图法；(c) THDR方法；(d) TDX方法; (e)基于三维梯度张量的N_ED_A方法；(f) Theta1方法；(g) Theta2方法，p取1.8. Fig. 7 Same as Fig.6 but for 1:50000 Bouguer gravity anomalies |
庐枞矿集区的重力变化相对较平缓，直接使用布格重力异常数据进行边界识别.由图 7b-g的边界识别结果可知：THDR方法边界比较模糊；Theta图法的边界线条较发散；TDX方法的边界比Theta图法收敛，细节更多；基于三维构造张量的N_ED_A方法和Theta1方法结果近似一致，Theta1方法的结果相对更清晰；Theta2方法识别的主要边界与Theta图法和TDX方法结果类似，边界线条更加收敛，在细节方面，Theta图法和TDX方法中识别的一些较小的呈半封闭的边界，在Theta2方法中呈完整的封闭形态.庐枞矿集区的重力实测数据变化范围幅度不大，不存在很显著的突变，因此不同的边界识别方法对重力数据的计算结果较为相似，相对而言Theta2方法的边界识别结果在保持了清晰的同时边界更加收敛.
4.4 地质解释从庐枞矿集区边界识别结果与实测资料对比可知，Theta2方法效果相对最好，准确刻画了重、磁实测数据的主要边界，没有产生明显的额外边界，识别的边界更为收敛.本文将Theta2方法的重、磁实测数据边界识别结果与矿集区地质图相叠合，如图 8所示.
|
图 8 庐枞矿集区地质和磁力(a)、重力(b)边界识别结果叠合图(地质图据吕庆田等，2014) Fig. 8 Superosed maps of the Luzong ore district showing geology and edge detection results of magnetic (a) and gravity (b) data (geologic map after Lü et al., 2014) |
前人的研究已经证实：庐枞矿集区的负磁异常区与白垩系红盆沉积和老地层分布区吻合，正磁异常区与火山岩和侵入岩分布区吻合，边界清晰.重力局部高异常与古生界-中生界地层、闪长岩等对应，火山岩区中的局部高重力异常或与基底老地层隆起有关，或与中基性侵入体有关；火山岩区之外的局部重力低与白垩系红盆沉积或第四系沉积对应，火山岩区内部的局部重力低是由于正长岩和二长岩岩基、岩盖引起(刘彦等，2012；吕庆田等，2014；祁光等，2014；严加永等，2014).
郯庐断裂带在边界识别结果中得到非常清晰的揭露，识别的磁力边界位置与郯庐断裂带的地质位置十分吻合，在重力边界识别中也有一条近似平行于郯庐断裂的NE向线条，但位置偏东南.郯庐断裂带两侧物性存在较大差异，郯庐断裂带以西表现为极高阻、低密度、高磁性的特征，与大别超高压岩层(UHP)变质带的中酸性岩浆活动有关(吕庆田等，2014).据庐枞矿集区NW向跨郯庐断裂带的人工地震剖面Lz09-02解译结果(吕庆田等，2014)，在CDP点801东南侧为郯庐断裂带主干断裂，与磁力边界识别的断裂位置接近；而重力识别的边界与CDP点1201附近的小型断裂位置接近，推测该位置可能是郯庐断裂带的一条次级断层.
在庐枞矿集区东南侧，识别出NE-NNE走向的一条重力边界，该位置与深部探测解译的CTF (长江逆冲断裂)(吕庆田等，2014)毗邻，推测是沿(长)江断裂的一部分，在磁力边界识别中只有局部地段有所显示，如黄梅尖岩体的东部边界，重、磁边界在位置上正好能相互沟通衔接.从物性方面分析，庐枞矿集区的东部紧邻火山岩盆地为沿江隆起带，由逆冲的古生界(Pz)-中生界(Mz)地层组成，密度范围在2.61~2.79 g·cm-3(据祁光等，2014)，相对盆地内部的二长岩、粗安岩、正长岩等密度范围在2.47~2.62 g·cm-3(据刘彦等，2012)，密度平均要高约0.155 g·cm-3，产生了局部高重力异常，两者之间的边界在重力数据边界识别结果中非常明显.
目前，对于庐枞盆地的边界分布，有不同的认识(任启江等，1993；董树文等，2010；高锐等，2010；吕庆田等，2014).深反射地震剖面揭示庐枞火山岩盆地是一个沿着北东向罗河-缺口边界断裂向东发育的“耳状”非对称盆地(董树文等，2010；高锐等，2010).深部探测项目在庐枞矿集区实施的五条综合地球物理探测剖面解译结果认为，庐枞火山岩盆地呈不对称的“箕状”，四周由向盆地倾斜的边界断裂围限(吕庆田等，2014).庐枞盆地是一个以中基性火山岩为主的盆地，据庐枞矿集区立体探测成果报告(中国地质科学院矿产资源研究所，20141))岩、矿石磁化率变化范围0~0.2 SI, 统计钻孔实测磁化率存在两个峰值，第一个多在0~400 m范围内，磁化率平均值约为0.02254 SI，最高值在0.12 SI左右，磁性主要由安山岩引起；第二个峰值在大多集中在600~1100 m范围内，磁化率平均值为0.04234 SI，最高值超过0.2 SI，磁性主要由矿体引起，矿体附近的围岩剩磁和磁化率也相对较大，说明存在一定程度的蚀变和矿化.结合图 5的统计结果，盆地内部与正磁异常有关的火山岩和侵入岩的磁化率约为0.02~0.04 SI，盆地外围与负磁异常有关的白垩系红盆沉积和老地层的磁化率约为0.00001~0.00002 SI，两者磁化率相差103数量级，因此其边界在磁力实测数据的识别结果上非常明显(图 8a)，构成一个近似封闭的环形.边界所围限的区域呈NE走向，西界位于罗河-缺口断裂附近，北界处于龙桥、焦冲一线，东界由黄梅尖岩体东侧经过钱家铺、杨家市、陶家巷延伸至枞阳，边界整体与地质单元一致，局部有出入.对比庐枞盆地的地质边界和磁力数据经上延1 km后识别的磁力边界，可见边界总体呈现向盆地内部收敛的特征，与综合地球物理探测的结果(吕庆田等，2014)相吻合.
1)中国地质科学院矿产资源研究所.2014.《庐枞矿集区立体探测技术与深部成矿预测示范》项目研究报告.
重力识别边界与庐枞盆地的地质边界相关性不太强，局部相吻合，如盆地东北侧的黄梅尖岩体.黄梅尖岩体为正长岩，密度范围在2.53~2.63 g·cm-3；黄梅尖岩体东侧为逆冲的古生代-中生代灰岩地层组成，密度范围在2.61~2.79 g·cm-3；黄梅尖岩体北侧为早、中侏罗系碎屑岩，密度范围在2.62~2.72 g·cm-3(据祁光等，2014)；黄梅尖岩体密度平均低约0.1~0.12 g·cm-3，产生了局部低重力异常，地质体的边界在重力数据边界识别结果中非常明显.
在盆地内部，位于阳家墩、孙家坂、白柳之间的重力识别封闭边界，在地质图上与白垩系浮山组(K1f)位置一致，周边主要为白垩系双庙组(K1sh)、西北侧为白垩系砖桥组(K1zh).据物性资料，白垩系浮山组(K1f)粗面岩的密度范围在2.40~2.51 g·cm-3，而白垩系双庙组-龙门院组(K1sh-K1l)粗安岩、碎屑岩的密度范围在2.56~2.66 g·cm-3(祁光等，2014).即白垩系双庙组(K1sh)和白垩系砖桥组(K1zh)物性差异不显著，在上述物性统计中合并归入白垩系双庙组-龙门院组(K1sh-K1l)，但围限在其中的白垩系浮山组(K1f)密度平均低约0.155 g·cm-3，产生了局部低重力异常，地质体的边界在重力数据边界识别结果中非常明显.有意义的是在浮山组(K1f)重力边界围限区域也识别出一组磁力封闭边界，较大的封闭椭圆位于阳家墩、店桥、白柳之间，较小的封闭椭圆位于钱家铺的西北角；由航磁异常图可知，该处处于盆地内部的低磁异常区.地质工作已经证实该处是一个破火山口(任启江等，1993)，推测后期充填了中酸性侵入岩(A型花岗岩或正长岩)，在物性上具有低磁、低密度的特点，据安徽省地质调查院了解该处未施工钻孔，没有准确的岩性和物性资料.
在庐枞矿集区的外围还分布有一些环形或近似呈封闭形态的识别边界，这些封闭区域可能与深部的隐伏岩体相对应，如磁力数据在沙溪地区所识别的封闭区域处于地表出露的沙溪岩体西南侧，可能代表了沙溪岩体在深部向西南侧延伸；磁力数据在塘串河地区识别出较好的封闭区域，与钻孔揭露的隐伏基性-超基性岩体相对应；磁力与重力数据均在裕丰村西北侧识别出封闭区域，最近的勘探成果也表明该区域深部存在隐伏的闪长玢岩，且发育有矽卡岩型铁铜矿化.
4.5 小结通过物性分析来解释边界识别结果，是将地球物理方法和地质认识相结合的纽带.通过上述工作，我们发现庐枞盆地内部火山岩和侵入岩(磁化率约为0.02~0.04 SI)较盆地外围白垩系红盆沉积和老地层(磁化率约为0.00001~0.00002 SI)的磁化率高约103数量级，因此盆地边界在磁力实测数据的识别结果上非常明显.对庐枞矿集区的重力识别边界封闭区域进行物性分析，可知存在密度异常的地质体相对于周边地层或岩体有≥0.1 g·cm-3的密度差即可以产生很好的重力边界识别结果，如文中的黄梅尖岩体相对密度低约0.1~0.12 g·cm-3，白垩系浮山组(K1f)相对密度低约0.155 g·cm-3等.上述关系是在庐枞矿集区1：5万尺度重磁数据分析中得到，当使用具体矿床更大比例尺(如1：1万)数据时，可能会有差异.通过本文设置的理论重力模型(图 1)试验，当地质体具有≥0.01 g·cm-3的密度差时理论上可以识别出边界，但这与地质体的形状、规模和埋深有关，仅供理论参考，不具有实际地质意义.
5 结论(1)本文基于Theta图法定义了新的边界识别滤波器，首先对位场数据求取垂向一阶导数，对垂向一阶导数再使用Theta图公式进行边界识别，本文的Theta2方法对Theta1方法的垂向高阶导数项进行了变换，更加清晰的描述地质体的主要边界.通过模型验证，该滤波器显著压制了Theta图法对深部地质体边界的放大作用，较好地平衡了深部和浅部边界.通过与传统的Theta图法、THDR方法和TDX等方法对比，本文定义的边界识别滤波器能够更清晰的圈定出地质体的水平边界位置.
(2)以长江中下游成矿带庐枞矿集区为例，开展了1：5万重磁数据的处理分析，结果表明：重磁数据的检测结果精确刻画了郯庐断裂带的位置，磁力边界位于郯庐断裂带主干断裂附近；区域重力在长江北岸识别的重力边界，推测是沿(长)江断裂的一部分；庐枞盆地的磁力数据检测边界整体与盆地的地质边缘一致，明确了边界断裂在深部倾向盆地内部；识别出庐枞盆地外围一系列环形边界，这些边界封闭区域与最新勘探发现的深部岩体及铁铜矿化体相对应，对于指导区域深部找矿工作有着重要的意义.
(3)物性是连接地球物理与地质的纽带和桥梁，通过对庐枞矿集区重、磁数据边界识别异常区的物性对比分析，我们认为只有了解更多的物性数据，多学科综合分析，才能合理、正确的解释边界识别结果蕴含的地质意义.
致谢感谢SinoProbe-03项目组吕庆田研究员、严加永副研究员的帮助！文中图件由Wessel and Smith开发的免费作图软件GMT绘制.感谢审稿专家及编辑提出的宝贵修改意见！
| Blakely R J, Simpson R W. 1986. Approximating edges of source bodies from magnetic or gravity anomalies. Geophysics, 51(7): 1494-1498. DOI:10.1190/1.1442197 | |
| Chang Y F, Liu X P, Wu Y C. The Copper-Iron Belt of the Lower and Middle Reaches of the Changjiang River (in Chinese).Beijing: Geological Publishing House, 1991: 127-308. | |
| 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. | |
| Cordell L. 1979. Gravimetric expression of graben faulting in Santa Fe Country and the Espanola Basin, New Mexico.//New Mexico Geological Society Guidebook, 30th Field Conference, 59-64. | |
| Dong S W, Xiang H S, Gao R, et al. 2010. 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), 26(9): 2529-2542. | |
| Evjen H M. 1936. The place of the vertical gradient in gravitational interpretations. Geophysics, 1(1): 127-136. DOI:10.1190/1.1437067 | |
| Fan Y, Qiu H, Zhou T F, et al. 2014. LA-ICP MS zircon U-Pb dating for the hidden intrusions in the Lu-Zong basin and its tectonic significance. Acta Geologica Sinica (in Chinese), 88(4): 532-546. | |
| Gao R, Lu Z W, Liu J K, et al. 2010. A result of interpreting from deep seismic reflection profile:Revealing fine structure of the crust and tracing deep process of the mineralization in Luzong deposit area. Acta Petrologica Sinica (in Chinese), 26(9): 2543-2552. | |
| Guo D, Yan J Y, Lü Q T, et al. 2014. 3D density mapping constrained by geological information:model study and application. Acta Geologica Sinica (in Chinese), 88(4): 763-776. | |
| Hornby P, Boschetti F, Horowitz F G. 1999. Analysis of potential field data in the wavelet domain. Geophysical Journal International, 137(1): 175-196. DOI:10.1046/j.1365-246x.1999.00788.x | |
| Hsu S H, Sibuet J C, Shyu C T. 1996. High-resolution detection of geologic boundaries from potential-field anomalies:An enhanced analytic signal technique. Geophysics, 61(2): 373-386. DOI:10.1190/1.1443966 | |
| Li L L, Huang D N, Han L G. 2014. Application of the normalized total horizontal derivative (NTHD) in the interpretation of potential field data. Chinese J. Geophys. (in Chinese), 57(12): 4123-4131. DOI:10.6038/cjg20141223 | |
| Liu D J, Hong T Q, Jia Z H, et al. 2009. Wave number domain iteration method for downward continuation of potential fields and its convergence. Chinese J. Geophys. (in Chinese), 52(6): 1599-1605. DOI:10.3969/j.issn.0001-5733.2009.06.022 | |
| Liu D J, Hong T Q, Liao X T, et al. 2012. Iterative solution of integral equation for potential field continuation from an irregular surface to a horizontal plane. Chinese J. Geophys. (in Chinese), 55(10): 3467-3476. DOI:10.6038/j.issn.0001-5733.2012.10.030 | |
| Liu Y, Lü Q T, Yan J Y, et al. 2012. The structure of Luzong ore district and its metallogenic indication from gravity and magnetic information. Acta Petrologica Sinica (in Chinese), 28(10): 3125-3138. | |
| Lü Q T, Liu Z D, Tang J T, et al. 2014. Upper crustal structure and deformation of Lu-Zong ore district:Constraints from integrated geophysical data. Acta Geologica Sinica (in Chinese), 88(4): 447-465. | |
| Ma G Q, Du X J, Li L L. 2013. New edge detection method of potential field data-enhanced horizontal derivative method. Progress in Geophys. (in Chinese), 28(1): 402-408. DOI:10.6038/pg20130145 | |
| Ma G Q, Huang D N, Yu P, et al. 2012. Application of improved balancing filters to edge identification of potential field data. Chinese J. Geophys. (in Chinese), 55(12): 4288-4295. DOI:10.6038/j.issn.0001-5733.2012.12.040 | |
| Ma G Q, Liu C, Huang D N. 2015. The removal of additional edges in the edge detection of potential field data. Journal of Applied Geophysics, 114: 168-173. DOI:10.1016/j.jappgeo.2015.01.007 | |
| 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. DOI:10.1016/0926-9851(94)90022-1 | |
| Nabighian M N. 1984. Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms:Fundamental relations. Geophysics, 49(6): 780-786. DOI:10.1190/1.1441706 | |
| Pan Y M, Dong P. 1999. The Lower Changjiang (Yangzi/Yangtze River) metallogenic belt, east central China:Intrusion-and wall rock-hosted Cu-Fe-Au, Mo, Zn, Pb, Ag deposits. Ore Geology Reviews, 15(4): 177-242. DOI:10.1016/S0169-1368(99)00022-0 | |
| Qi G, Lü Q T, Yan J Y, et al. 2014. 3D geological modeling of Luzong ore district based on priori information constrained. Acta Geologica Sinica (in Chinese), 88(4): 466-477. | |
| Rajagopalan S, Milligan P. 1994. Image enhancement of aeromagnetic data using automatic gain control. Exploration Geophysics, 25(4): 173-178. DOI:10.1071/EG994173 | |
| Ren Q J, Wang D Z, Xu Z W, et al. 1993. Formation and development of the Mesozoic Lujiang-Zongyang volcanic-structural depression in Anhui Province and their relation to mineralization. Acta Geologica Sinica (in Chinese), 67(2): 131-145. | |
| Roest W R, Verhoef J, Pilkington M. 1992. Magnetic interpretation using the 3-D analytic signal. Geophysics, 57(1): 116-125. DOI:10.1190/1.1443174 | |
| Sertcelik I, Kafadar O. 2012. Application of edge detection to potential field data using eigenvalue analysis of structure tensor. Journal of Applied Geophysics, 84: 86-94. DOI:10.1016/j.jappgeo.2012.06.005 | |
| Tang Y C, Wu Y C, Chu G Z, et al. Geology of Copper-Gold Polymetallic Deposits in the along Changjiang Area of Anhui Province (in Chinese).Beijing: Geological Publishing House, 1998: 10-21. | |
| 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-125. DOI:10.1190/1.1651454 | |
| Wang J, Meng X H, Li F. 2015. Improved curvature gravity gradient tensor with principal component analysis and its application in edge detection of gravity data. Journal of Applied Geophysics, 118: 106-114. DOI:10.1016/j.jappgeo.2015.04.013 | |
| Wijns C, Perez C, Kowalczyk P. 2005. Theta map:Edge detection in magnetic data. Geophysics, 70(4): L39-L43. DOI:10.1190/1.1988184 | |
| Yan J Y, Lü Q T, Chen M C, et al. 2015. Identification and extraction of geological structure information based on multi-scale edge detection of gravity and magnetic fields:An example of the Tongling ore concentration area. Chinese J. Geophys. (in Chinese), 58(12): 4450-4464. DOI:10.6038/cjg20151210 | |
| Yan J Y, Lü Q T, Chen X B, et al. 2014. 3D lithologic mapping test based on 3D inversion of gravity and magnetic data:A case study in Lu-Zong ore concentration district, Anhui province. Acta Petrologica Sinica (in Chinese), 30(4): 1041-1053. | |
| Yuan F, Zhou T F, Fan Y, et al. 2008. Source, evolution and tectonic setting of Mesozoic volcanic rocks in Luzong basin, Anhui Province. Acta Petrologica Sinica (in Chinese), 24(8): 1691-1702. | |
| Yuan Y, Huang D N, Yu Q L, et al. 2014. Edge detection of potential field data with improved structure tensor methods. Journal of Applied Geophysics, 108: 35-42. DOI:10.1016/j.jappgeo.2014.06.013 | |
| Yuan Y, Huang D N, Yu Q L. 2015. Using enhanced directional total horizontal derivatives to detect the edges of potential-field full tensor data. Chinese J. Geophys. (in Chinese), 58(7): 2556-2565. DOI:10.6038/cjg20150730 | |
| Yuan Y, Yu Q L. 2015. Edge detection in potential-field gradient tensor data by use of improved horizontal analytical signal methods. Pure and Applied Geophysics, 172(2): 461-472. DOI:10.1007/s00024-014-0880-1 | |
| Zhai Y S, Yao S Z, Ling X D. Regularities of Metallogenesis for Copper (Gold) Deposits in the Middle and Lower Reaches of the Yangtze River Area (in Chinese).Beijing: Geological Publishing House, 1992: 1-20. | |
| Zhang J S, Gao R, Li Q S, et al. 2010. Characteristics of gravity and magnetic field of Luzong volcano basin and its periphery. Acta Petrologica Sinica (in Chinese), 26(9): 2613-2622. | |
| Zhou T F, Fan Y, Yuan F. 2008. Advances on petrogensis and metallogeny study of the mineralization belt of the Middle and Lower Reaches of the Yangtze River area. Acta Petrologica Sinica (in Chinese), 24(8): 1665-1678. | |
| Zhou T F, Fan Y, Yuan F, et al. 2010. Temporal-spatial framework of magmatic intrusions in Luzong volcanic basin in East China and their constrain to mineralizations. Acta Petrologica Sinica (in Chinese), 26(9): 2694-2714. | |
| Zhou W N, Du X J, Li J Y. 2013. The limitation of curvature gravity gradient tensor for edge detection and a method for overcoming it. Journal of Applied Geophysics, 98: 237-242. DOI:10.1016/j.jappgeo.2013.09.008 | |
| 常印佛, 刘湘培, 吴言昌. 长江中下游铜铁成矿带.北京: 地质出版社, 1991: 127-308. | |
| 董树文, 项怀顺, 高锐, 等. 2010. 长江中下游庐江-枞阳火山岩矿集区深部结构与成矿作用. 岩石学报, 26(9): 2529–2542. | |
| 范裕, 邱宏, 周涛发, 等. 2014. 安徽庐枞盆地隐伏侵入岩的LA-ICP MS定年及其构造意义. 地质学报, 88(4): 532–546. | |
| 高锐, 卢占武, 刘金凯, 等. 2010. 庐-枞金属矿集区深地震反射剖面解释结果-揭露地壳精细结构, 追踪成矿深部过程. 岩石学报, 26(9): 2543–2552. | |
| 郭冬, 严加永, 吕庆田, 等. 2014. 地质信息约束下的三维密度填图技术研究及应用. 地质学报, 88(4): 763–776. | |
| 李丽丽, 黄大年, 韩立国. 2014. 归一化总水平导数法在位场数据解释中的应用. 地球物理学报, 57(12): 4123–4131. DOI:10.6038/cjg20141223 | |
| 刘东甲, 洪天求, 贾志海, 等. 2009. 位场向下延拓的波数域迭代法及其收敛性. 地球物理学报, 52(6): 1599–1605. DOI:10.3969/j.issn.0001-5733.2009.06.022 | |
| 刘东甲, 洪天求, 廖旭涛, 等. 2012. 位场曲化平积分方程的迭代解. 地球物理学报, 55(10): 3467–3476. DOI:10.6038/j.issn.0001-5733.2012.10.030 | |
| 刘彦, 吕庆田, 严加永, 等. 2012. 庐枞矿集区结构特征重磁研究及其成矿指示. 岩石学报, 28(10): 3125–3138. | |
| 吕庆田, 刘振东, 汤井田, 等. 2014. 庐枞矿集区上地壳结构与变形:综合地球物理探测结果. 地质学报, 88(4): 447–465. | |
| 马国庆, 杜晓娟, 李丽丽. 2013. 位场数据边界识别的新方法-增强型水平导数法. 地球物理学进展, 28(1): 402–408. DOI:10.6038/pg20130145 | |
| 马国庆, 黄大年, 于平, 等. 2012. 改进的均衡滤波器在位场数据边界识别中的应用. 地球物理学报, 55(12): 4288–4295. DOI:10.6038/j.issn.0001-5733.2012.12.040 | |
| 祁光, 吕庆田, 严加永, 等. 2014. 基于先验信息约束的三维地质建模:以庐枞矿集区为例. 地质学报, 88(4): 466–477. | |
| 任启江, 王德滋, 徐兆文, 等. 1993. 安徽庐枞火山-构造洼地的形成、演化及成矿. 地质学报, 67(2): 131–145. | |
| 唐永成, 吴言昌, 储国正, 等. 安徽沿江地区铜金多金属矿床地质.北京: 地质出版社, 1998: 10-21. | |
| 严加永, 吕庆田, 陈明春, 等. 2015. 基于重磁场多尺度边缘检测的地质构造信息识别与提取-以铜陵矿集区为例. 地球物理学报, 58(12): 4450–4464. DOI:10.6038/cjg20151210 | |
| 严加永, 吕庆田, 陈向斌, 等. 2014. 基于重磁反演的三维岩性填图试验-以安徽庐枞矿集区为例. 岩石学报, 30(4): 1041–1053. | |
| 袁峰, 周涛发, 范裕, 等. 2008. 庐枞盆地中生代火山岩的起源、演化及形成背景. 岩石学报, 24(8): 1691–1702. | |
| 袁园, 黄大年, 余青露. 2015. 利用加强水平方向总水平导数对位场全张量数据进行边界识别. 地球物理学报, 58(7): 2556–2565. DOI:10.6038/cjg20150730 | |
| 翟裕生, 姚书振, 林新多. 长江中下游地区铁铜矿床.北京: 地质出版社, 1992: 1-120. | |
| 张季生, 高锐, 李秋生, 等. 2010. 庐枞火山岩盆地及其外围重、磁场特征. 岩石学报, 26(9): 2613–2622. | |
| 周涛发, 范裕, 袁峰. 2008. 长江中下游成矿带成岩成矿作用研究进展. 岩石学报, 24(8): 1665–1678. | |
| 周涛发, 范裕, 袁峰, 等. 2010. 庐枞盆地侵入岩的时空格架及其对成矿的制约. 岩石学报, 26(9): 2694–2714. | |