基于三维反演的重力异常分离方法研究

2020-12-16陈煜烽邱姜歆邓居智

东华理工大学学报(自然科学版) 2020年5期

陈辉，陈煜烽，邱姜歆，邓居智

(东华理工大学放射性地质与勘探技术国防重点学科实验室，江西南昌 330013)

重力布格异常包含诸多信息，是地下半空间所有密度不均匀体所引起的重力异常综合响应(曾华霖，2005)。在实际重力数据处理中，人们往往只关注某些目标地质体(如断裂构造、基底起伏或成矿岩体等)的重力异常，由此必须从综合布格重力异常中提取单个或若干目标体引起的重力异常进行反演解释，即异常分离。然而，由于实际地质体间存在错综复杂的组合关系，提取某个地质体或几个目标体的重力异常相对较困难(管志宁等，2002)。因此，重力异常分离成为当前重力数据处理关键环节和难点问题。

经过几十年来的深入研究，目前具有代表性的重力异常分离方法有多项式拟合(Martín et al.，2011;Ojo et al.，1997)、滤波技术(Guo et al.，2013;郭良辉等，2012;Keating et al.，2011)、延拓技术(骆遥，2011;Meng et al.，2009)、趋势分析(羊春华，2005;Abdelrahman et al.，1991;王彦国等，2017)、小波分解(Singh et al.，2015;Xu et al.，2009;杨文采等，2001)、经验模态分解(曾琴琴等，2010)、最小曲率法(纪晓琳等，2015)、主成分分析(Pei et al.，2012;Zhang et al.，2009)和有限元法(Agarwal et al.，2010)等。其中，多项式拟合、滤波技术和延拓技术是较为成熟的方法。这些方法能够利用特定理论模型和实测数据验证其有效性，但需要通过大量试验选择合理参数组合，如小波分析中母函数、滤波技术中的滤波因子选择，延拓技术的延拓高度选择等。在实际数据处理中，有时可以采用多种方法在不同参数条件下进行对比试验，并结合地质资料进行优化组合(Martínez-Moreno et al.，2015;刘彦等，2012;Sedighi et al.，2009)。近年来，随着重力三维反演技术逐渐走向成熟，在地质和地球物理等相关资料约束下可获取地下三维密度分布信息(Lelièvre et al.，2009;Williams，2008;Shin et al.，2006;陈辉等，2015;Li et al.，1998a)。Li等(1998b)提出基于三维反演的磁测数据异常分离技术，通过理论模型和实际资料验证有效性。该方法能够比较直观地提取不同深度的磁异常。

笔者介绍重力三维反演的基本原理，给出利用三维反演技术实现重力异常分离的方法和技术流程;进而通过理论模型试验探讨方法的有效性;最后将该方法应用于江西九瑞矿集区实测布格重力异常处理中，提取了不同深度的局部重力异常，并结合地质资料进行分析和解释。

1 理论与方法

1.1 三维重力反演基本原理

笔者采用加拿大哥伦比亚大学UBC-GIF小组研发的三维重力反演软件Grav3D。其基本原理是基于Li等(1998a)提出的反演算法。该反演方法构建的目标函数为

φ=φd+μφm

(1)

式中，φd为数据目标函数，φm为模型目标函数，μ是平衡模型适用性和数据拟合误差的一个权重参数。数据目标函数定义为观测数据与模型正演响应权重差的L-2范数，即

φd=‖Wd(dfwd-dobs)‖2

(2)

式中，Wd为对角阵，由观测数据标准误差的倒数构成，dfwd为正演模拟数据，dobs为观测数据。模型目标函数φm定义为反演密度模型和先验信息参考模型的权重差，即

(3)

式中，mref为先验信息参考模型，m为每次迭代的反演模型，ws、wx、wy和wz为权重函数。其中，空间权重函数ws控制各个网格单元和参考模型网格单元的物性特征，而权重函数wx、wy、wz控制单元模型各个方向的平滑度，wr(z)是深度或距离的权重函数，αs、αx、αy和αz为各方向的相关性系数。

为了克服重力数据三维反演垂向分辨率差及多解性问题，通常将已知物性、钻孔和地质信息作为约束条件融合到反演中。由式(3)可知，通过预先设置参考模型mref、边界范围以及权重函数三种方式可以实现约束反演。本研究只采用无约束三维反演。

1.2 基于三维反演的重力异常分离

为了实现布格重力异常分离，首先对重力数据区域进行扩边处理，并利用三维反演得到地下三维密度模型，然后根据需要提取密度模型中的目标区块进行正演模拟，得到局部重力异常，最后通过与布格重力异常相减得到区域重力异常。

(1)数据准备。首先将实测数据通过常规数据改正和处理后得到布格重力异常△g;为了消除人为噪声影响、提高三维反演的稳定性，将布格重力异常△g向上延拓0.5～1倍点距。另外，三维反演效果也会受到地形起伏的影响，因此需要足够密度的地形数据加以约束。

(2)数据扩边。为了削弱边界效应对三维反演的不利影响，获得深部密度分布信息，通常需要对反演区域进行扩边，为此研究区域的观测布格重力数据也需进行相应的扩边处理。通常利用研究区域周边大面积小比例尺的数据与观测布格重力数据进行融合拼接，或者利用重磁处理软件(如Oasis Montaj)中的FFT算法对观测布格重力数据进行扩边得到。扩边范围通常需要达到数据观测区域的3倍以上。

(3)初始模型建立。合理的初始模型是重力三维反演的关键。初始模型由研究区域、扩边区域和延拓区域三部分组成。其中，研究区域是为了提取浅部局部异常，需要获得较为可靠的密度模型，横向网格大小通常为布格重力数据观测点距或网格大小的0.5～2倍，垂直方向采用随深度指数增加的网格，第一层网格厚度一般取横向网格的0.5～1倍。扩边区域是为了提取中深部的局部异常，横向网格尺寸为研究区域网格大小的1～3倍，而垂向网格尺寸与研究区域一致。延拓区域是为了减小反演的边界效应，提高研究区域和扩边区域内反演密度模型的可靠性，初始横向网格大小与扩边区域一致，然后呈指数递增，垂向网格与研究区域一致。另外，初始模型的剩余密度一般取为零。当研究区地形起伏较大时，初始模型需要加入地形数据。

(4)重力三维反演通常采用三维无约束反演。在布格重力异常反演数据中需给定数据误差(数据不确定性)，通常依据数据质量确定。当数据质量难于确定时，可根据多次反演试验的结果进行设定。

(5)局部异常和区域异常分离。根据研究目标需求在反演三维密度模型中将对应区块的密度体单独提取，并利用三维重力正演软件进行正演模拟，得到局部重力异常数据，最后利用布格重力异常数据减去局部重力异常数据获得区域重力异常数据。

2 理论模型试验

为了验证本算法的有效性，可设计一个区域重力异常和局部重力异常叠加密度模型(图1a)。其中，区域密度模型为一个半径3.0 km的深部球体，顶部埋深为10 km，质心位于(0，0)，密度为1.8 g/cm3。局部密度模型为5个不同大小、位置、埋深和密度的异常体。其中，块体1和2密度为1.5 g/cm3、块体3和4为1 g/cm3、块体5为1.2 g/cm3(表1)。每个块体位置和埋深如图1a所示。重力异常数据观测范围为7 km×7 km，采样间隔0.25 km。图1b给出重力异常叠加分布特征。由图1可见，5个异常体所产生的局部重力异常混杂在埋藏较深的大异常体所产生的区域重力异常之中，异常体边界难于识别。

表1 理论模型参数Table 1 The parameter of synthetic models

对叠加重力场数据进行三维反演获取地下三维密度模型，然后提取研究区范围内3.5 km以浅密度体进行正演模拟得到局部异常，并与叠加异常相减得到区域重力异常。同时采用二阶趋势分析和解析延拓(向上延拓3 500 m)进行重力异常分离。图2引用三种分离方法得到的区域重力异常与区域密度模型正演响应的对比结果。趋势分析法(图2b)分离的区域异常形态与正演响应(图2a)基本一致，但最大异常幅值相差最大(0.595 mGal，约为7.44%)。向上延拓分离(图2c)的区域重力异常形态与正演响应形态有一定差别，在近北东方向有一定拉伸，呈现为近椭圆形，最大异常幅值相差0.567 mGal，约为7.09%。利用三维反演方法(图2d)分离得到的区域重力异常形态与正演响应形态基本一致，最大异常幅值相差最小，为0.416 mGal，约为5.20%。

图3引用三种分离方法得到的局部重力异常及局部密度模型正演响应的对比结果。三种方法提取的局部重力异常形态(图3b，c，d)与正演响应(图3a)基本一致，但从细节上看，趋势分析和向上延拓分离出的区域异常边界不清晰。从局部重力异常最大幅值来看，趋势分析相差最大(为0.396 mGal，约为4.44%);向上延拓相差0.377mGal，约为4.43%，而三维反演相差最小(0.294 mGal，约为3.44%)。

综上所述，无论从形态还是从分布特征上，趋势分解、解析延拓和三维反演方法分离的区域和局部重力异常与正演响应均呈现基本一致特征，均取得较好效果。然而，从细节上利用重力三维反演分离区域和局部重力异常与理论值更接近;从数值上利用三维反演提取局部和区域重力异常与正演响应最大幅值相差最小，具有更高准确度。由此说明基于三维反演的重力异常分离算法较传统的异常分离方法更有效。

3 实例分析

九瑞矿集区位于扬子地块中部，下扬子断裂拗陷带的西缘，北临大别山隆起，南靠江南褶皱造山带，是长江中下游铜铁金硫多金属成矿带的重要组成部分。区内发育“多位一体”(斑岩型、矽卡岩型、块状硫化物)的城门山铜多金属矿床，以矽卡岩型、块状硫化物矿床为主的武山铜矿，以及与次火山岩作用有关的洋鸡山金矿等重要金属矿床(邓居智等，2015)，如图4所示。该区位于九瑞—彭泽复式向斜的西段，主要褶皱自北向南有:邓家山—通江岭向斜、界首—大桥背斜、横立山—黄桥向斜、大冲—丁家山背斜、乌石街—赛湖向斜、长山—城门湖背斜、新塘向斜。区内断裂主要为北东向、北西向、北东东向。北东向断裂属郯庐—赣江构造带中庐山穹断束西部边缘断裂带的一部分，对九瑞地区构造格局起着重要的控制作用。该区盖层出露的地层比较齐全，从震旦系至第四系，除缺失志留系部分地层外，其余均有出露。区内岩浆活动较强烈，小岩体星罗棋布，岩浆岩活动有燕山期和喜山期两个岩浆旋回。该区发现已知矿床近30个，受成矿地质、断裂构造和岩体共同控制。

该区已开展1∶5万比例尺重力测量。研究区位置见图4，而图5给出九瑞矿集区1∶5万布格重力异常图。布格重力异常整体呈现由东北向西南逐渐降低，中部出现明显北东向密集梯度带特征，与丰山洞-城门山深部断裂相吻合。该深部断裂带对九瑞矿集区构造格局起着重要的控制作用。布格重力异常无法揭示地层和构造分布特征，需要进行异常分离，提取局部重力异常。为此，对研究区数据通过FFT算法进行扩边处理，然后实现无约束条件带地形三维反演得到地下三维密度分布，并根据研究目标提取不同深度的重力异常进行分析和解释。

为了进一步说明本算法的有效性和可行性，利用三阶趋势分析、低通滤波(波长为70 km)和向上延拓(上延6 km)等传统方法进行异常分离。从图6可见，4种方法提取的区域重力异常均表现为西南低，沿东北方向逐渐增大，在东北角呈现为高值的特征。其中，趋势分析等值线表现为近北西向的弧形，数值变化激烈，三维反演、低通滤波和向上延拓重力异常等值线表现为近北西方向直线，幅值变化较平滑且基本接近。通过与该区区域地质及深部构造对比发现，基于三维反演的重力异常分离方法更加符合实际地质情况。

从图7可见，4种方法提取的区域重力异常整体呈西南向北东延伸，北西为高值，东南为低值。另外，区内存在三条明显高异常区，分别为中南部近东西向分布的大规模带状高异常区、中部近西北向延伸的高异常带和东北部近西北向高异常带，与该区背斜的分布特征吻合。西北部、东南部及东北部的低异常区则与向斜分布相对应。同时，已知矿床基本出现在局部重力高值异常的边界处。三维反应提取的局部重力异常幅值变化较小，高值异常边界清晰。

为进一步了解不同深度重力异常分布特征，分别提取了三维密度模型中<0.5 km、0.5～1 km、1～5 km、5～10 km的密度体进行正演模拟，得到不同深度的局部重力异常(图8)。在浅部(图8a)可以看出研究区存在东西向和近东西向三条重力高异常带，分别对应大桥-宝山复式背斜、丁家山-大冲背斜和城门山-长山背斜;四条东西向和近东西向的重力高异常带分别对应于江岭-邓家山复式向斜、黄桥复式向斜、赛城湖-乌石街复式向斜和新圹向斜;背斜和向斜交替出现为研究区地层的主要特征。随着深度增加(图8b)高值异常范围逐渐增大，低值异常区逐渐变小。在深部(图8c，d)表现为东北重力异常值高、西南低的特征。由此表明背斜在深部逐渐消失，而向斜逐渐融合为一体。结合矿床和矿点分布特征，发现向斜和背斜接触带(五通组和黄龙组交接处)为成矿有利区。