矿区重力三维物性反演的参考模型构建方法及应用

2020-07-26周俊杰陈聪喻翔高玲举陈涛

物探与化探 2020年4期

周俊杰，陈聪，喻翔，高玲举，陈涛

(核工业北京地质研究院，北京 100029)

0 引言

地球物理三维反演是获取深部地质结构信息的重要技术之一，是近年来矿产勘查工作攻深找盲的热点研究问题之一[1-2]。重力异常对横向密度变化敏感度高，在矿产勘查中常用于推断地层起伏以及局部岩体分布等，但垂向分辨能力较差，通常需要引入先验信息来约束，以期改善重力反演结果的质量。因此，约束的构建和引入是重力三维反演中的关键问题。Li 和Oldenburg[3]、Williams[4]、姚长利等[5]较早的介绍了重力三维约束反演的方法原理，通过理论模型和典型实测案例说明了重力三维反演中约束的重要性。随后，关于重力三维约束反演的相关研究不断涌现。祁光等[6]、兰学毅等[7]使用交互方式研究了位场反演中的先验信息约束问题，并在矿集区得到成功应用。郭冬等[8]对三维密度填图进行了研究，认为地表地质信息和钻孔信息对提升重力反演的垂向分辨能力具有帮助。罗凡等[9]在深部磁铁矿勘查中利用先验信息约束的重磁反演圈定了矿体范围。陈长敬等[10]采用电法资料构建约束进行重力三维反演，利用电法垂向分辨率高的特点改善了重力反演结果。

重力反演的约束方式有多种，如参考模型约束[11]、物性范围约束[12]、方向约束[13]、聚焦约束[14]以及结构约束[15]等，可适用不同种类、不同性质的先验信息的引入。相较于其他方式，参考模型具有理论简单，稳定可靠等特点，但其构建流程较为繁琐。使用参考模型进行约束的主要步骤为：广泛收集已知资料并转化为地下空间物性信息，构建参考模型并衡量相应权重，最后代入反演计算得到结果。在实际应用中，常见的可利用先验信息包括地表地质图、物性统计数据、钻孔资料以及其他地球物理成果等。笔者针对参考模型的构建及约束的引入问题进行分析，通过理论算例和实测资料分析探讨重力三维约束反演的效果，结果显示参考模型法可有效提升反演质量，值得在实践中推广应用。

1 重力三维约束反演方法

重力三维约束反演利用重力数据及相关地质信息约束，计算地下密度三维分布特征。Li 和Oldenburg[16]提出的重力三维约束反演方法在重力勘探领域得到广泛应用，其基本思想是将地下三维空间剖分为包含若干直立长方体单元的规则网格，通过求解每个单元物性值来获取地下密度分布情况。本文采用以此理论为基础的UBC-GIF GRAV3D程序进行反演计算。为便于先验信息的引入，反演算法在满足数据拟合的同时也引入了模型拟合，其目标函数可简写为：

φ(m)=‖Wd[d0-g(m)]‖2+

β[α1‖Wm(m-mref)‖2+α2‖D(m-mref)‖2]，

(1)

其中：第一项为重力数据拟合目标函数；第二项为模型目标函数，内含最小模型项和光滑模型项，β为正则化参数。两项目标函数均用L2范数定义。Wd为观测数据标准差对角阵，d0为重力数据向量，m为待求模型向量，g为正演算子。模型项中，mref为参考模型向量，D为三维光滑度矩阵，α1,α2分别为最小项和光滑项的权重因子。在数据项和模型项目标函数的共同作用下，反演迭代所得结果将在满足数据拟合的基础上，同时符合参考模型的特征，并具有光滑特性。参考模型作为已知向量，已经包含在目标函数中，需要在反演前预先构建。

2 参考模型构建方法

参考模型是指和待反演模型具有相同网格架构，且内含地下空间已知物性信息的三维数据体。根据式(1)，反演迭代解将尽可能的贴近参考模型特征，以此来达到约束的效果。为表明参考模型约束的作用，使用典型理论模型进行分析测试。图1所示为测试采用的理论模型，其东向、北向和垂向网格规模为100×1×30，其中剖面方向单元长度为10 m，走向方向单元长度1 000 m，为降低边界效应影响，在模型体边缘部位均进行了扩边处理。模型由浅地表低密度体和起伏基底构成，各地质单元密度参数如表1所示。取平均密度为2.6×103kg/m3，对剩余密度进行三维反演。首次反演采用无约束方式，即不考虑先验信息，参考模型为0的情形，如图2a所示；假设已知基底的大致深度信息，将参考模型设置为层状模型，代入反演所得结果如图2b所示；假设除基底信息外，还已知浅地表密度信息，将其纳入参考模型，代入反演所得结果如图2c所示。3种情况均使用GRAV3D程序的默认反演参数，仅对比使用不同参考模型时反演效果。

图2 不同参考模型下反演结果对比Fig.2 Comparison of inverted results under different reference model constraints

表1 理论密度模型参数设置

图1 理论密度模型Fig.1 Synthetic density model

从理论模型反演试验可以看出，参考模型约束对于改善反演结果质量具有成效。当不使用参考模型(即参考模型为0)时，反演仅在正则化条件下对观测数据进行拟合，受制于重力反演体积效应的影响，其结果仅在横向上较好地分辨出基底隆起位置，但其层位难以进行解释。当已知测区内基底平均深度时，可根据盖层和基底密度特点构建一个层状模型，反演结果则会在满足数据拟合的同时尽可能符合参考模型特征，相较于无约束反演，在分辨出基底隆起横向位置的同时还能反演出基底和盖层的分界面，但基底隆起处密度界面模糊，不利于隆起的进一步解释。当已知浅地表信息时，可在层状参考模型基础上再引入局部低密度体。值得注意的是，两者的可靠程度不同，通常浅地表信息更加准确，可赋予高权重值，反演结果会更加贴近高权重位置的参考模型密度特征。在本例中，引入地表信息后，模型权重取地表30 m内从100衰减至10，其余部位权重均为1。反演结果可看出，地表密度分布更符合设定信息，而层状基底的隆起部位更加明显，埋深也与已知数值相近。可见，将不同的地质信息以参考模型方式引入到重力反演中，可有效提高反演质量，得到更加可靠的反演结果。

在实际应用中，重力三维反演的关键在于参考模型的构建。一般可先在不引入先验信息的情况下开展无约束反演，获得地下密度分布初步认识；之后结合测区地质图、物性资料和相关勘探成果，根据岩体和地层的大致分布特点和密度统计规律构建参考模型。对于较可靠的地表信息可赋予高权重，而无法确认是否准确的深部推测资料可赋予低权重。此外，钻孔信息也可赋予较高权重加以利用。实验表明，通常情况下高权重设置为低权重的100倍左右较为适宜，反演结果既能具有参考模型特征，又能接受迭代扰动，满足正常的观测数据拟合要求。

3 张麻井矿区应用实例

3.1 矿区概况

笔者结合张麻井矿区实例来分析参考模型在重力三维反演中的构建及应用。张麻井矿区(图3a)位于沽源火山断陷盆地西南缘，区域面积约64 km2，盆地基底为太古宇红旗营子群中深变质岩系，研究区中未出露；盖层为侏罗系上统及白垩系下统火山熔岩、火山碎屑岩系。其中，区内西北部大面积出露张家口组二段(J3z2)粗面岩，中部及南部多为张家口组三段(J3z3)流纹岩、角粒熔岩，西南部出露张家口组三段(J3z3)凝灰岩，东北部出露花吉营组(K1h1)安山岩等。区内地表第四系覆盖广阔，为下赋地层界线的划定带来一定困难。矿床主要受构造、岩浆活动、热液蚀变等因素综合控制，摸清成矿地质结构对深部找矿具有较大价值[17-18]。

1—第四系；2—花吉营组安山岩；3—张家口组二段粗面岩；4—张家口组三段第8层流纹岩；5—张家口组三段第7层角粒熔岩；6—张家口组三段第6层凝灰岩；7—张家口组三段第5层板状流纹岩；8—张家口组三段第1～4层凝灰岩；9—次流纹斑岩；10—花岗斑岩；11—硅化带；12—地质界线；13—实测断层；14—推测断层；15—剖面；16—钻孔

在测区内开展了1∶25 000大比例尺重力测量，对布格重力异常进行分离后得到剩余布格重力异常，其形态可基本反映2 km深度范围内地层走势和岩体分布[19]。研究区西北部呈现高值异常，表征地下高密度地层较厚，而中部低值异常则表征低密度地层较厚。在研究区开展了岩石密度测量统计，并结合前人研究成果，获得本区岩石物性主要特征(表2)。盖层侏罗系上统及白垩系下统岩石密度较低，其中张家口组流纹岩(J3z3)密度为2.34×103kg/m3，凝灰岩(J3z3)密度为2.41×103kg/m3，花吉营组安山岩(K1h)密度为2.41×103kg/m3。张家口组角粒熔岩(J3z3-7)和粗面岩(J3z2)密度较高，分别为2.51×103kg/m3和2.52×103kg/m3。由于区内基底资料较为缺乏，据以往资料太古宇混合变质基底(Arh)密度为2.69×103kg/m3，远高于火山岩密度。此外，次流纹斑岩(Cλπ)也显示出密度高的特征，为2.59×103kg/m3。

3.2 无参考模型约束反演

根据测区基本情况，将地下介质剖分为138×152×35的规则网格，单元长度为50 m。为削弱边界效应，各方位向外各扩展10个网格。使用2.57×103kg/m3作为背景密度值，结合物性统计数据，将剩余密度限定在(-0.5 ～0.5)×103kg/m3之间。反演使用GRAV3D程序，除物性范围和参考模型外均使用默认参数。图4为无参考模型约束的反演结果。从密度切片可见，地下密度分布在西北较大，东南较小，与剩余布格重力异常特征相符，但密度分层不明显，深部密度较低，与地质认识不符。此外，浅地表密度形态较为杂乱，和岩体露头、出露地层特征差异很大。这说明无约束反演结果分辨能力较弱，可靠性低，不利于开展进一步的地质解释工作。

图4 无参考模型约束反演结果外观(a)和切片(b)Fig.4 Perspective(a) and sliced view(b) of inverted result without reference model constraint

3.3 参考模型构建

3.3.1 物性填充法构建地表约束

根据岩石密度特征和地质图相关层位几何信息可构建地表物性约束。利用自主研发的地表信息数字化软件将岩层密度值填充到对应位置，其网格大小和前述三维网格的平面设置相同。首先根据地质图将第一层地表物性进行数字化填充，之后顺层推断。推断完毕后，将平面网格按照深度关系合为一体，组装为三维网格模型。图5所示为根据测区地质图对地表岩性进行数字化编码，进而整合为三维地表参考模型结果，结合表2中统计的各类岩性密度值进行物性填充，整合后即可用作参考模型。

图5 基于地质图的物性数字化和三维模型构建Fig.5 Geological-map-based density digitization and 3D model construction

表2 研究区岩石密度特征

3.3.2 分层赋值法构建地层约束

在测区内收集了密度界面反演结果和二维人机交互正反演结果[19]，可将其作为深部地层信息引入到反演中。首先将密度界面按照三维网格参数进行重采样，随后转化为网格模型，在对应地层单元上赋予相应的物性代号值(图6a、b)，之后结合表2中物性统计值进行填充。对于多层地层，需注意地层整合及穿插关系，断层处的地层错断也需注意，以免发生物性填充错误。由于深部地层信息通常为推断结果，在反演时应赋予不高于地表信息的权重。如有钻孔信息，则应在钻孔处赋予较高的权重。

3.3.3 模型融合

地表信息和地层信息构建为三维模型后，将两者进行拓扑融合(图6c、d)，主要思路为权重大的物性单元将权重小的代替。此外还需注意地表和露头的连续性，可采用剖面扫描方式逐个进行检查。之后还需将相应的权重模型也进行融合，确保准确度高的物性单元配以高权重。

图6 基于地层信息的模型构建及其与地表模型的融合Fig.6 Strata-based model construction and fusion with surface model

3.4 参考模型约束反演结果

图7为对张麻井铀矿研究区的重力三维约束反演结果， GRAV3D程序参数设置除上述参考模型外，其他均与图4相同。反演初始均方差为1.3×107，经12次迭代反演，最终均方差为2.8×103，收敛曲线呈对数下降并趋于稳定，数据拟合非常理想。此外，该模型也具有参考模型提供的地表地质和地层起伏基本特征，同时密度的空间分布信息更为丰富，对中浅层密度不均匀体也有很好的体现。图8a、b分别为使用和不使用参考模型得到的反演结果对比。尽管无约束结果也能描述基底的隆起，但在深部无法反映出火山盖层，浅部物性分布杂乱；和无约束情况相比较，约束反演结果则有了很大改善：反演结果反映的地层露头形状和地质图相匹配，密度幅值符合统计规律，呈现出随深度增加，密度逐渐增大的特点，且地层具有明显的分层特征，呈现西北隆起，东南凹陷的特点，与重力异常特征相符。此外，中浅层隐伏岩体在反演剖面上也反映较为明显，呈现出次高密度异常。

综合分析反演结果可知，地表第四系密度较低，在研究区广泛分布，但其埋深有限，多在50 m以内。北侧出露张家口组二段粗面岩(J3z2)，呈现出密度中等偏高的特点，且密度随深度加大而逐渐变大。南侧出露多为流纹岩(J3z3)，仅西南侧为凝灰岩(J3z3)，其密度比粗面岩较低，角粒熔岩(J3z3-7)在西大山以西部分出露，和围岩相比密度值略高。东北部出露安山岩(K1h)，其密度相对较低。从图7b剖面可知，基底和火山盖层均呈现北侧隆起、中南部凹陷，南侧隆起的特征，其中南侧隆起比参考值更高，但低于北侧隆起。三维密度模型对岩体分布也有反映，在张麻井矿床处显示局部高密度特征，表征了次流纹斑岩的赋存；西大山西侧角粒熔岩也显示出局部次高密度特征。西大山隆起处也出现次高密度特征。根据SN向密度剖面AA′结果(图8)，测区北部为地层隆起，地表有张家口组二段粗面岩出露；而中部和南部主要为张家口组三段流纹岩。结合岩层密度特征和地质特征，可推断F45下盘张家口组二段粗面岩顶界面埋深约为海拔600～1 000 m之间；太古界混合变质基底顶界面埋深在海拔100 m以下。在F45上盘，变质基底顶界面埋深在海拔100～1 000 m之间。据钻孔资料，张麻井ZKB-1在海拔810 m处见张家口组二段粗面岩，密度剖面解释结果与钻孔揭露位置相符。综上所述，参考模型约束反演由于先验信息的引入而变得更加可靠。相比无约束反演的分辨能力更好、质量更高；相较于收集的界面反演资料，反演结果信息更为丰富，在体现地层错断和刻画中浅层岩体赋存方面优势更为明显。