Parker法视磁化强度填图方法研究及应用

秦 熠，杨宇山，刘天佑，朱 丹

(1.中国地质大学 地球物理与空间信息学院，湖北 武汉 430074；2.四川省煤田地质局 地质工程勘察设计研究院，四川 成都 610000)

1 引 言

视磁化强度填图方法是一种利用磁异常反演地下地质体磁性强弱的方法，利用航磁资料进行视磁化强度填图可以有效确定磁性体边界，划分磁性地层，特别对于存在大规模隐伏岩体地区具有良好的应用效果，作为研究区域构造的辅助手段，已被广泛地应用于地质构造解释与研究之中[1]。

1973年Parker[2]将傅里叶变化引入位场正演计算，提出了连续单界面模型快速正演方法，即Parker法。随后，Oldenburg[3](1974)在Parker法基础上推导出了频率域密度界面迭代反演方法，由于该方法解决了过去空间域离散模型界面反演方法计算速度慢的问题，因此被广泛地应用于起伏界面的计算。王万银等[4](1993)在单界面重力场正反演基础之上研究了更加符合实际情况的双界面重力场正反演问题，并提出了保证反演快速收敛的措施。

视磁化强度填图兴起于20世纪80年代，Pilkington[5](1989)根据Parker界面正反演方法推导出了磁化强度反演公式，并应用于Athabasca盆地。早期视磁化强度填图，尤其在化极处理之后存在负值的问题，难以解释。管志宁等[6](1990)研究了视磁化强度填图出现负值的原因，并提出一种采用立体角滑动窗口迭代反演方法的视磁化强度填图，应用于秦巴地区。这种方法适用于地形起伏不大，观测平面水平的航磁资料处理。张贵宾等[7](1996)从理论与模型算例两方面对Parker公式为基础的频率域位场反演方法做了深入研究，提出新的快速正则化反演方法，并将其推广为单界面无限延深视密度和视磁化强度填图方法。Pilkington[8](1998)提出了等效地形改正方法，实现了起伏地形起伏观测面视磁化强度填图，但是并未解决磁化强度存在负值的问题。眭素文[9](2004)结合前人无负值填图与等效地形改正两种方法，实现了小区域中高山区起伏测量面航磁资料无负值视磁化强度填图。此后，大量学者将视磁化强度填图方法应用于岩体与构造研究中，朱晓颖[10](2007)利用视磁化强度填图方法研究了内蒙古北山地区断裂构造，圈定侵入岩体，深化了研究区区域地质、成矿规律和控矿特征的认识。刘天佑等[11](2012)运用密度填图反演推断了西藏朗县秀沟工区蛇纹石化橄榄岩的范围以及铬铁矿与矿化体的重磁远景。王赛昕等[12](2017)针对攀枝花地区进行视磁化强度填图，利用磁化强度与全铁质量分数的经验公式获得了区内全铁质量分数分布图。

尽管前人对起伏界面模型视磁化强度填图方法进行了大量研究，但是受到方法以及资料的限制，在实际应用时仍然主要采用水平界面模型进行反演，对于区域跨度广地形起伏大的地区会带来较大误差。针对于这个问题，本文提出一种功率谱分析约束下的双界面起伏模型Parker法视磁化强度填图新思路，步骤如下：①根据化极磁异常形态，将等轴状或长轴状孤立正磁异常单独分区；②分别截取每一个小分区磁异常以及研究区全区异常，采用功率谱分析方法反演得到每一个小分区的上界面与下界面埋深以及全区平均上下界面深度；③根据反演得到的小分区上下界面深度构建上界面与下界面起伏形态，其中负磁异常等未分区域深度取全区界面平均深度；④以功率谱分析构建的上下起伏界面为约束，采用双界面Parker视磁化强度填图法反演得到更加符合实际的磁性分布。

2 起伏界面视磁化强度填图原理

Parker法频率域磁性界面正反演方法是一种已知磁化强度参数，反演界面深度的方法。由Parker法磁性界面正反演公式可知：

(1)

由(1)式可知，若已知均匀密度分布值M，则由(1)式可以计算界面起伏h，也可以把(1)式推广为已知界面起伏h，求磁化强度分布，这就是Parker法视磁化强度填图的方法原理。假设地下磁性层磁性界面深度已知，需要反演求解界面之间磁化强度的水平分布，其中磁化强度M是随水平位置变化的函数。类似于Parker法磁性界面深度反演方法中将界面起伏看作平均深度与相对于平均深度起伏两部分组成一样，把磁化强度也看作是平均磁化强度与相对平均磁化强度的变化量组成，即：

M=M0+ΔM

(2)

那么磁异常也可以看作是由平均磁化强度产生的磁异常与相对平均磁化强度的变化量产生的磁异常两部分组成：

(3)

其中：

(4)

(5)

(6)

3 理论模型研究

考虑实际磁性体延深不同的情况，本文建立厚层与薄层岩体模型，分析不同上界面与下界面取值对视磁化强度反演的影响；然后，采用功率谱分析法反演理论模型的上下界面深度，以得到的界面深度约束反演磁化强度，比较其与理论模型的误差，验证功率谱分析法结合视磁化强度填图反演的有效性；最后，建立理论模型模拟地下不同磁化强度岩体分布，比较起伏界面模型视磁化强度填图与水平界面模型视磁化强度填图应用效果。

3.1 界面深度对填图的影响分析

通过单个直立长方体模型建立厚层岩体与薄层岩体理论模型，分析上下界面取值对反演视磁化强度的影响，理论模型正演参数选取见表1，正演磁异常如图1所示。视磁化强度填图计算时取磁性界面水平(图2)，分别绘制厚、薄层岩体深度变化率与相对误差曲线图，横坐标为相对界面深度变化率，其中上界面(图3中蓝线)深度变化率为Δh(选取上界面与理论上界面深度差)与理论上界面深度H的比值，下界面(图3中红线)深度变化率为Δl(选取下界面与理论下界面深度差)与理论模型延深长度L的比值，纵坐标为磁性体分布范围内反演的磁化强度均值与理论磁化强度的相对误差。对于上界面深度影响，固定理论下界面深度不变，只改变上界面深度，然后进行视磁化强度填图计算。当选取上界面深度位于理论上界面之上时，其深度变化率为正。对于下界面深度影响，则固定理论上界面深度不变，只改变下界面深度，再进行视磁化强度填图计算。当选取下界面深度位于理论下界面之上时，其深度变化率为负。

综合分析界面深度变化与相对误差曲线图，可以得到以下几点认识：①不论上界面还是下界面，相同界面深度变化下薄层岩体的磁化强度相对误差始终大于厚层岩体，说明上下界面深度取值对薄层岩体影响更大，对厚层岩体影响相对较小；②不论薄层岩体还是厚层岩体，在30 %深度变化范围内，相同界面深度变化下上界面相对误差要大于下界面相对误差，表明上界面深度取值对填图结果影响要大于下界面深度取值；③当上界面位于理论上界面之上或下界面位于理论下界面之下时，曲线斜率逐渐变小，相对误差变化随深度变化率的增大逐渐变小，即深度变化影响逐渐减弱。当上界面位于理论上界面之下或下界面位于理论下界面之上时，曲线斜率逐渐变大，深度变化影响逐渐增强。

理论模型表明，采用实际资料进行视磁化强度填图时，磁性上界面深度选取至关重要，其对反演计算的磁化强度值影响较大。此外，薄层岩体由于其延深较小，其下磁性界面深度选取不准也会造成反演计算较大误差。

表1 厚层与薄层岩体模型参数

图1 理论模型正演磁异常Fig.1 Theoretical model positive magnetic anomaly

图2 界面深度选取示意图Fig.2 Interface depth selection diagram

图3 界面深度对填图结果影响Fig.3 Interface depth affects the mapping results 注：当选取界面较理论界面深时，上界面深度变化率为负，下界面深度变化率为正

3.2 功率谱分析法确定界面深度效果对比

利用航空磁测资料进行视磁化强度填图时，若能较为准确确定磁性体上下界面深度就能取得良好的反演效果。目前常用的反演磁性体深度方法[13]有欧拉齐次方程法和功率谱分析法、矩谱法等。欧拉齐次方程法不需要预先知道反演场源的形状和参数，能够自动快速估算场源埋深，但是很难确定下界面深度，而功率谱分析法利用场源在对数功率谱曲线上表现出不同梯度的线性下降特征反演场源深度，不仅能够反演磁性体上界面深度，也可以反演磁性体底界面深度，结合功率谱分析对磁性界面的约束与视磁化强度填图法可以得到更加符合实际情况的反演结果。此外，对于孤立形体，为了准确地确定孤立形体的下界面深度，可以采用矩谱法求质心，再通过上界面深度与质心深度获得下界面深度。

采用上述厚层与薄层岩体模型，使用功率谱分析法反演磁性体的上下界面深度，低频段反映深部界面信息，中高频段反映浅部界面信息。根据厚层岩体功率谱曲线(图4a)得到模型上界面深度为37.1 m，绝对误差7.1 m，下界面深度140.8 m，绝对误差10.8 m，然后采用反演得到的界面深度计算厚层岩体边界范围内平均磁化强度为8.49 A/m，与理论磁化强度的相对误差为6 %。由薄层岩体功率谱曲线(图4b)得到模型上界面深度为35.1 m，绝对误差5.1 m，下界面深度61.3 m，绝对误差11.3 m，反演得到薄层岩体平均磁化强度为7.16 A/m，与理论磁化强度的相对误差为11 %(见表2)。

比较反演计算值与理论值，功率谱分析法反演的上界面深度要比下界面深度更准确，但是反演得到的界面深度都要比实际界面深度更深一些。对于厚层岩体来说，采用此方法反演的上下界面误差相较其延深来说影响较小，因此反演得到的磁化强度值也更加接近真实值。但是对于薄层岩体来说，根据该方法确定的上下界面误差相对其延深来说影响较大，造成其视磁化强度填图相对误差要大于厚层岩体。总的来说，运用功率谱分析法能够较为有效地确定磁性体上下界面深度，以此为约束得到的视磁化强度填图相对误差不大于15 %，理论模型应用效果较好。

图4 功率谱分析Fig.4 Power spectrum analysis

3.3 起伏界面模型与水平界面模型填图效果对比

表3 理论模型参数

图5 不同界面模式填图效果对比Fig.5 Comparison of different interface mode mapping effects

以不同埋深直立长方体模型模拟实际磁性岩体分布(图5a)，对比水平界面模型与起伏界面模型视磁化强度反演不同磁性岩体的效果。其中1号模型为厚层弱磁性岩体，磁化强度为2 A/m，2号模型为薄层强磁性岩体，磁化强度为5 A/m，3号模型为厚层中等磁性岩体，磁化强度为3.5 A/m，各磁性体上界面与下界面深度不同，模型参数取值见表3。从正演磁异常(图5b)来看，高磁异常范围与模型分布范围相对应，其中1号模型磁异常略大于2号模型，异常幅值在500～550 nT之间，3号模型磁异常最强，异常幅值在650～700 nT之间，正演磁异常大小与理论模型磁化强度并无直接对应关系。

采用水平界面模型计算时，使用功率谱分析法由全区磁异常得到上界面平均深度30 m，下界面平均深度70 m，选取该参数作为水平界面模型约束进行视磁化强度填图，反演结果见图5(e)。当采用起伏界面模型计算时，按照区内高磁异常个数将全区分为3个异常区，功率谱分析法分别反演每一个异常区域内顶底界面埋深，分区之外的区域顶底界面深度则由全区磁异常功率谱分析求得，然后将反演得到的顶界面与底界面深度拟合连接成面得到磁性层上下界面起伏形态(图5c，图5d)，以得到的起伏界面模型深度为约束对正演磁异常做视磁化强度填图结果(图5f)。

由水平界面模型填图结果可知,1号与2 号岩体磁化强度相似，磁化强度值在2.8～3.8 A/m之间，属于中等磁性岩体，3号岩体磁化强度值略高于1、2号岩体，磁化强度值在3.2～4.2 A/m之间，属于中高磁性岩体。显然该结果与理论模型存在较大偏差，未能正确区分岩体磁性的相对强弱，计算的磁化强度值也并无对应关系。对比起伏界面模型反演结果可知1号岩体磁性相对较弱，磁化强度分布在2～2.8 A/m，2号岩体磁性最强，磁化强度在4～5 A /m，3号岩体磁性强于1号岩体，弱于2号岩体，磁化强度在3～4 A/m。反演计算的磁化强度值与理论模型磁化强度值相对应，能够较为明确地划分为弱、强、中等磁性三类不同磁化强度岩体。从理论模型反演结果来看，起伏界面模型相较水平界面模型不论反演计算的磁化强度值还是不同岩体磁性相对强弱关系都与理论模型更相似，说明结合功率谱分析法反演起伏界面深度约束下的视磁化强度填图能够更为准确地圈定岩体分布范围，划分不同磁化强度岩体。

4 小秦岭地区的应用

4.1 研究区地质概况

小秦岭地区是秦岭成矿带地质研究、资源勘查的主要工作区之一，区内已查明十多处大型钼矿、金矿等金属矿床，是我国金、银、钼、铅等多金属矿的重要产地[14]。研究区位于华北陆块南缘东秦岭钼矿带以西，出露地层主要为太古宇太华群变质岩系、中元古界熊耳群火山岩系和中元古界官道口群、新元古界栾川群、下古生界陶湾群等一套碳酸盐岩组合以及中新生代陆相盆地沉积。受板块间深大断裂和秦岭造山带长期活动影响，区内构造复杂，断裂与褶皱发育，区域构造线主要以东西向展布，断裂构造可以分为近东西向、北东向和北西向三组。该区经历元古代-古生代-中生代多期岩浆活动，岩浆岩十分发育(图6)，其中主要以中生代和古生代酸性岩为主，地表有较大面积出露。研究表明，区内岩浆岩与钼矿、金矿等多金属矿床形成有关，已查明金属矿床大多分布在岩体内部或其边部附近[15]。

4.2 地球物理特征

由于不同时期岩浆侵入作用与火山活动强弱差异较大，造成了地层与岩石不同的磁性特征。前人物性统计资料表明(表4)，秦岭地区大部分沉积岩地层磁性较弱，只有太古代及元古代地层作为古老的结晶基底具有较高的磁性，但极不均匀，磁性变化大[16]。区内岩浆岩磁性特征复杂，不同时期存在明显差异(表5)。早古生代时期形成的岩浆岩磁性最高，中生代次之，古生代最低，其中中生代二长花岗岩和闪长花岗岩较其他时期形成的花岗岩磁性更强，随着岩浆岩从酸性到基性的变化，岩石磁性逐渐增强。

1.中生代酸性岩;2.古生代酸性岩;3.中性岩;4.基性岩;5.金矿点;6.钼矿点;7.铁铜铅矿点;8.断裂 图6 小秦岭岩浆岩分布图(修改自《秦岭成矿带地质矿产图》李宗会2013)Fig.6 Map of Xiaoqinling magmatic rocks

表4 秦岭地区地层磁性

表5 秦岭地区岩浆岩磁性

本文采用秦岭东段地区1∶20万航空磁测资料进行处理解释，平均飞行高度500 m。航磁异常(图7a)表现为东西向和北东向连续条带状异常，显示出小秦岭地区受华北陆块与扬子陆块挤压形成的东西向的线性构造特征。经化极处理消除斜磁化影响(图7b)，磁异常极值北移，高磁异常分布范围与地表大面积出露的燕山期花岗岩对应关系较好；印支期花岗岩与沉积性地层磁性较弱，对应负磁异常或低值正磁异常；具有一定规模的中基性岩能够产生较强的磁异常，而小规模条带状中基性岩产生磁异常较弱。

1.中生代酸性岩;2.古生代酸性岩;3.中性岩;4.基性岩;5.金矿点;6.钼矿点;7.铁铜铅矿点 图7 小秦岭地区磁异常Fig.7 Magnetic anomalies in the Xiaoqinling area

4.3 小秦岭地区视磁化强度填图应用

由于区内钼、金、铅等金属矿床的形成与岩体关系密切，寻找隐伏岩体对于寻找金属矿具有重要指导意义。研究区内岩体分布，首先要提取浅部磁性体异常特征，由于小波多尺度分解对高频信号分辨率高的特点[17]，本文采用小波多尺度分解来分离磁异常，提取的剩余异常见图8(a)。

1.中生代酸性岩；2.古生代酸性岩；3.中性岩；4.基性岩；5.金矿点；6.钼矿点；7.铁铜铅矿点； 8.磁异常小分区；9.推测隐伏中生代花岗岩(C1、C2、C3为出露岩体；Y1、Y2、…Y11为隐伏岩体) 图8 小秦岭视磁化强度填图Fig.8 Magnetization map in Xiaoqinling area

由理论模型可知，对孤立磁异常分区，然后采用功率谱分析法确定磁性体顶底界面进行起伏界面视磁化强度填图比传统水平界面填图效果更优。本文将全区剩余磁异常划分为24个区域，分别对每一个孤立磁异常做功率谱分析，计算可以得到其对应的上下界面埋深。而分区之外的区域上下界面深度则由全区磁异常功率谱分析得到，反演平均上界面深度为1 462 m，下界面深度为8 642 m。此外，由于该地区岩浆岩体出露面积较大，为了提高起伏上界面选取精度，可以结合GPS地形高程数据来确定岩体出露地区起伏上界面的深度，GPS地形相对观测面转换的界面深度由航空磁测平均飞行高度来确定(图8b)。即研究区起伏上界面深度由两种方法共同获取，在岩体未出露的分区，上界面深度由功率谱分析法得到；在岩体出露的分区，上界面深度直接由GPS地形高程数据获取，得到的起伏上下界面如图8c、图8d。

以反演计算的起伏上下界面为约束，对剩余磁异常采用双界面模型进行填图，反演的得到磁化强度分布见图8(e)。不同类型岩体之间存在明显差异，古生代酸性岩和沉积地层大多分布在中低磁化强度区域，中性岩与基性岩本身规模较小，主要分布在中等磁化强度区域，中生代酸性岩则分布在中高磁化强度区域。蔡川-三要-官坡中高磁化强度条带区与地表大面积出露的东西向中生代花岗岩对应关系较好，反映了中生代花岗岩磁性相对较高的特征。华阴-潼关-灵宝高磁化强度带也与地表断续出露的中生代花岗岩对应，推测地层之下可能存在隐伏的中生代花岗岩带。五里川北部北西向长轴状高磁化强度区地表出露较大面积的中性闪长岩以及小面积中生代花岗岩，由中性岩填图磁性特征推测，中性岩之下可能隐伏有较大规模中生代花岗岩体。此外，由中生代酸性岩高磁化强度特征，共圈定区内11个隐伏中生代花岗岩体(图8f)。对比分析研究区内目前已查明多金属矿床分布位置，钼、金矿床与中生代花岗岩关系密切，集中分布在北西区域出露的中生代花岗岩体(C-1)边部，其他区域钼、金矿床相对较少，大多分散分布于圈定的隐伏中生代花岗岩体边部(Y-1、Y-2、Y-3、Y-4、Y-10)。隐伏中生代花岗岩体与围岩接触部位既为成矿有利区域，成矿可能性较高，可以作为后期找矿工作开展的重点方向。

5 结 论

本文以Parker法界面正反演公式为基础，推导了起伏界面模型视磁化强度填图迭代公式；针对地形变化剧烈地区提出一种起伏双界面模型快速视磁化强度填图的新思路，即结合功率谱分析法确定磁性起伏界面，然后以计算的起伏界面来约束视磁化强度填图；理论模型验证了该方法的有效性，相较传统水平界面模型提高了反演精度。将该视磁化强度填图方法应用于小秦岭多金属成矿带，圈定了区内与钼、金成矿密切相关的隐伏中生代花岗岩体11个，其中Y-1、Y-2、Y-3、Y-4、Y-10隐伏岩体边部均存在已查明的钼矿或金矿床与之对应，剩余未发现矿床点的隐伏岩体Y-5～Y-9、Y-11边部位置成矿可能性较高，这些岩体边部附近即为成矿有利位置，可以为今后金属矿产勘查工作提供参考。需要指出的是，本文仅采用航空磁测资料以及地质资料等作为推断依据，采用资料较为单一，得到的结论有待更进一步验证。此外，建议尝试将功率谱法求下界面深度改用矩谱法，先求质心再求下界面深度。