地球物理学报  2020, Vol. 63 Issue (7): 2737-2750   PDF    
空-地-井重力异常正则化协同密度反演方法
王泰涵1,2, 马国庆1,2, 熊盛青3, 刘燕戌3, 曾昭发1, 李丽丽1,2, 孟庆发1,2     
1. 吉林大学地球探测科学与技术学院, 长春 130026;
2. 吉林大学国家发展与安全研究院, 长春 130026;
3. 中国地质调查局自然资源航空物探遥感中心, 北京 100083
摘要:围绕综合利用空-地-井多源重力异常联合反演提升精度的目标,提出正则化协同密度反演方法,从而有效利用多维数据的横向和纵向变化特征提高反演精度,且无需先验信息的约束.此外,还提出利用奇异值谱和深度分辨率工具评价航空数据观测高度和层数对反演分辨率的影响.通过理论模型试验本文方法在空-地、地-井等不同情况下的应用效果,结果表明空-地-井重力异常正则化协同密度反演方法能获得更高分辨率、高精度的反演结果,且证明不同的钻孔位置对反演结果有较大的影响,从而可指导实际的空-地-井联合勘探.最后,利用文顿盐丘实际勘探数据进行空地和空地井重力数据协同反演,反演结果垂向分辨率明显优于地面观测数据反演结果,盐盖的顶面及中心埋深与前人研究和地质资料解释相吻合,验证了方法正确性及实用性,为推进我国空-地-井立体重磁勘探提供了重要的技术手段.
关键词: 空-地-井重力      协同密度反演      分辨率      文顿盐丘     
Joint regularized density inversion method of airborne, surface and borehole gravity anomaly data
WANG TaiHan1,2, MA GuoQing1,2, XIONG ShengQing3, LIU YanXu3, ZENG ZhaoFa1, LI LiLi1,2, MENG QingFa1,2     
1. College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China;
2. Institute of National Development and Security Studies, Jilin University, Changchun 130026, China;
3. China Aero Geophysical Survey and Remote Sensing Center for Natural Resources, Beijing 100083, China
Abstract: In order to comprehensively utilize the joint inversion of airborne-surface-borehole multi-source gravity data anomalies to improve the accuracy, a regularized density inversion method was proposed, so as to effectively utilize the horizontal and vertical variation characteristics of multi-dimensional data to improve the inversion accuracy without the constraint of prior information. In addition, singular value spectrum and depth resolution plot tool were used to evaluate the effect of observation height and flight layer numbers on inversion resolution. Theoretical model was used to test the application effect of the method in airborne-surface, surface-borehole and other conditions. The results showed that joint regularized density inversion of airborne-surface-borehole gravity anomaly can obtain higher resolution and precision, and proved that different drilling positions have a greater influence on inversion results, which can guide the actual airborne-surface-borehole joint exploration. Finally, using actual field data of Vinton salt dome to carry out airborne-surface and airborne-surface-borehole joint inversion, the vertical resolution of inversion result was obviously better than the inversion result of surface data. The top depth and center depth of cap rock was consistent with predecessors' research and geological interpretation, which verifies the correctness and practicability of the method, and provides an important technical means for promoting three-dimensional airborne-surface-borehole gravity and magnetic exploration in China.
Keywords: Airborne-surface-borehole gravity    Joint density inversion    Resolution    Vinton salt dome    
0 引言

重力法在矿产及油气资源勘探和生产中一直发挥着重要作用, 其能够直接感知地下质量的不均匀分布所产生的响应.三维密度反演(Li and Oldenburg, 1998Portniaguine and Zhdanov, 1999Zhdanov et al., 2004; Zhdanov, 2009; Geng et al., 2014秦朋波和黄大年, 2016)是常用的重力资料定量解释方法, 在矿产和油气资源勘探中扮演着重要的角色.然而, 传统的地面重力异常的反演方法存在着垂向分辨率不足的问题, 为实际资料解释带来多解性和不确定性.由于重力场以距场源距离的二次方衰减, 这种对近地表略敏感的异常特征会导致反演对较深部敏感度的丢失, 即深部成像结果分辨率较差.因此, 一种常用策略是在反演中加入深度权函数约束来抵消这种固有衰减, 使场源成像在更合理的深度上(Li and Oldenburg, 1996Zhdanov, 2002).另外还可联合多种数据协同反演来达到这一目的, 例如基于重力仪器测量物理量不同获得的重力与重力梯度数据和重力张量梯度多分量数据协同反演, 以及由于重力仪测量环境与空间位置不同获得的航空、地面和井中观测数据协同反演等.简单的说协同反演是指用不同数据体, 围绕同一研究目标, 协同一致完成目标体场源计算的过程, 由于多元数据带来比单一地面数据更加丰富的信息, 有助于减少多解性, 同时提高解的精度及准确性.

井中重力测量因其靠近地下目标体的优势, 可以在地下空间更接近于解释人员感兴趣和关注的目标体, 获得其随深度分布的异常数据, 所以长期以来在井中重力仪器和解释方法方面得到了相应的发展(Lautzenhiser, 1983LaFehr, 1983MacQueen, 2007).近些年随着适用于较窄钻孔的小井眼重力仪被开发出来(Nind et al., 2007Seigel et al., 2010), 在矿产和油气资源勘查领域的应用越来越广泛, 同时许多解释技术也与仪器和测量技术并行发展.Chasseriau和Chouteau(2003)将最小构造法(Li and Oldenburg, 1998)应用于井中与地面重力数据反演, 模型试验清楚地表明加入井中数据后, 棱柱体的垂向分布范围得到了更为真实的再现.Sun和Li(2010)利用地面和井中重力数据对两条平行的倾斜脉状模型进行二维反演, 在迭代过程中不断抽取倾角信息联合井中的岩石物理属性作为约束, 能够明显提高倾斜构造反演的分辨率, 并表明在反演过程中井中数据对于提取倾角信息至关重要.Shamsipour等(2010)提出了一种基于协同克里金的反演方法对地面和井中重力异常进行反演, 穿过模型中心的钻孔数据对提高反演的深度分辨率具有重要作用.Mosher和Farquharson(2013)同样联合地面和井中重力数据进行最小构造法反演, 并指出在同时反演多个距离较近的井中重力数据时, 使用l1范数有助于准确地重建井中或井之间的精确密度分布.Geng等(2015)在重力梯度数据协同克里金反演中引入井重力数据获得了更为准确的密度模型.Rim和Li(2015)介绍了井中矢量重力正演模拟方法, 并对井中垂直和矢量重力测量数据对应包含的信息量进行了研究, 数值模拟表明, 矢量重力在异常质量定位、异常质量表征、油藏注液运动跟踪等方面优于垂直重力数据.Han等(2018)将地面和井中重力张量梯度数据应用于位场偏移成像中, 能够快速地恢复地下地质体信息, 提高常规偏移方法解释精度.

航空重力测量具有测量速度快、测量控制范围广和测区基本不受限制等优点(郭志宏等, 2008熊盛青, 2009孙中苗等, 2015), 能够反映更深的地质体及构造特征.另外, 改变不同飞行高度可以获得不同高度的测网数据, 能够近似反映地下不同深度场源分布特征.基于此原理, 张永明和张贵宾(2009)研究了二维重力正则化“等维”反演, 提高了重磁位场反演的垂向分辨率.Pilkington(2012)依据最佳设计思想, 对反演核函数特征值分解, 比较分析了不同高度观测情况下重力张量梯度分量组合所包含的信息量.王逸宸(2014)在聚焦反演中引入重力等维测网, 构建的等维聚焦反演方法明显抑制聚焦反演结果上漂现象.

航空和地面重力测量数据具有较高的水平分辨能力, 但缺少垂向分辨率, 与地面数据相比, 航空重力数据能够表征更深的场源信息, 而井中重力数据具有较高的垂向分辨能力, 但实际情况下受钻孔数量和方位限制而缺少水平分辨率.针对空地井立体勘探技术体系, 如何充分地发挥空地井测量数据的优势, 更好地协调不同数据之间协同反演是进一步提高重力方法探测能力和反演分辨率的关键.因此本文基于立体数据的空间分布特征, 建立空地井重力数据三维协同反演技术, 有效地利用其包含纵向变化的信息, 以减少多解性.利用模型计算讨论航空多飞行高度观测重力数据和井中重力数据的加入对反演结果分辨率的影响, 定性、定量地评价航空飞行高度和钻孔位置, 对实际勘探给予合理化指导和建议.最后通过实际案例对本文空地井协同反演方法进行了验证.

1 空地井重力数据联合反演方法 1.1 正演问题

首先, 在反演之前, 我们将地下三维模型空间划分为有限多的单位密度长方体, 三维网格中每个单元的相对密度即是反演所要求取的参数.根据不同的探测方式, 观测点可以位于地表以上、地表面或井中(钻孔内), 假设地下三维网格剖分为m个长方体单元, 给定任一种方式获得的n个观测数据, 重力异常与密度分布呈线性关系:

(1)

式中:g是观测重力异常;ρ是地下空间密度分布;G是重力正演核函数.

在计算正演核函数时定义笛卡尔坐标系中的x轴, y轴和z轴分别代表北向、东向和深度方向.根据Haáz(1953)提出的积分解, 由密度为ρ的单一棱柱体单元C在观测点O=(x0, y0, z0)处产生的重力异常为

(2)

式中:xi=x0ξi, yj=y0ηj, zk=z0ζk, i, j, k=1, 2;ξ, η, ζ定义长方体单元的八个角点;rijk=γ为万有引力常数, 一般取γ=6.67×10-11m3/(kg·s2).因此, 由密度分布在观测点处产生的重力异常可以近似为所有长方体单元产生响应的叠加

(3)

需要注意的是井中核函数正演时, 位于长方体单元边界位置的观测数据会导致解析正演公式不适用, 采取的措施是将其点沿z方向略微移动很小一段距离, 以避免公式(2)奇异性.

1.2 反演问题

三维重力反演问题往往是欠定的, 通常需要引入Tikhonov正则化约束来解决反演过程的非唯一性和不确定性, 恢复出较为合理的解, 反演问题表述如下:

(4)

式中:φ (ρ)是数据拟合函数, s(ρ)是稳定函数, 在本文中我们仅讨论一种近光滑的线性稳定泛函‖ρ22, α为正则化参数, 以自适应的方法在反演迭代中选取(Zhdanov, 2002Wang et al., 2017).Wd是数据误差加权矩阵, 是由观测数据噪声标准差倒数组成的对角矩阵;Wz为模型加权函数, 如果为单位矩阵I, 则相当于对解给予无加权约束.另外为改善位场反演深度分辨率, Wz被定义为一种基于深度加权方式的函数能够抵消重力核函数矩阵G随深度的固有衰减, 克服趋肤效应.此类加权方式的形式有很多(Li and Oldenburg, 1996Portniaguine and Zhdanov, 2002Zhdanov, 2002), 针对本文空地井数据联合反演的特点可选用广义的深度加权方式(Li and Oldenburg, 2000), 分别为基于核函数矩阵的加权方式:

(5)

或基于由网格单元与观测点之间距离定义的深度加权方式:

(6)

其中, β为施加权重的强度, 取值在0.5到1.5之间, 通常取1;Gij为第j个单位密度的网格单元对第i个观测点的正演响应;ΔVj为第j个网格单元的体积;rij为第i个观测点到第j个网格单元中心的距离;r0为常数, 通常选为反演网格间距的一半.以上两种加权函数仅与观测数据本身有关, 不受观测方向的变化且无需先验信息, 因此适用于空地井重力数据联合反演.

2 空、地测网数据协同反演

与地面重力测量相比, 航空重力测量具有采集效率高、测网设计限制少等优点.飞机在不同飞行高度下还可以获取多层的空中测网数据, 低高度观测数据反映浅部场源信息较明显, 较高的观测面数据能一定程度上反映较深部场源信息.这样, 通过航空测网数据的引入, 反演的数据体相当于反映地下不同深度的场源信息, 增加的量会提高常规重力反演结果的分辨率.本文通过利用奇异值谱和深度分辨率图像(Depth Resolution Plot, DRP)工具(Fedi et al., 2005)对空地联合反演的优势进行理论量化分析, 首先建立模型和网格剖分如图 1a, 地下网格单元为20×20×10, 改变不同的观测高度z值即可得到各高度观测面的核函数矩阵.以重力数据为例, 求取z=0, -100, -200, -300, -400, -500 m的六个观测面的重力正演核函数矩阵G, 对其进行奇异值分解如下:

图 1 (a) 模型网格剖分, (b)奇异值谱分析 Fig. 1 (a) Model grid meshing, (b) Spectral analysis of singular values

(7)

式中, U是左奇异向量, VT为右奇异向量, Σ是奇异值矩阵, 对角线元素为包含的奇异值σi, 其余元素为零.

离散反演的解可以表示成不同奇异值和奇异向量加权后作用于观测数据的结果.Hansen(1998)指出奇异值σi具有逐渐衰减为零的特征, 其谱是连续的, 系数矩阵维数(阶数)大小会影响奇异值的分布.因此, 将不同组合情况的奇异值谱成图, 若奇异值幅值变大, 则表征有更多能量较大的奇异值(特征值)参与了反演的重建, 使反演计算更为准确.为了便于比较, 所有组合情况的奇异值谱只展示前400个(地面核函数包含奇异值个数)序号, 小序号表征能量较大的奇异值(图 1b).随着航空观测数据的增多, 参与反演重建的奇异值增多, 奇异值幅值有所增加.从奇异值谱曲线形态来看, 增加多高度观测数据能改变奇异值谱的衰减程度, 当地面数据联合两层空中测网数据形成三层观测数据时, 再增加观测面数据对于曲线形态改变不大, 考虑其可能对反演结果的提升不大.

要恢复出具有有用深度信息的解, 必须有足够的右奇异向量(Fedi et al., 2005).右奇异向量包含了反演重建中的重要信息, 分析右奇异值向量可以理论分析反演包含深度分辨率的相关信息.对于每一个右奇异向量vi, 分配一个长度为nz(地下纵向剖分层数)的向量si, 定义Vi表示vi在三维空间(nx×ny×nz)的展开, 这种可视化的深度分辨率图像(Fedi et al., 2005Paoletti et al., 2014)计算方法为:

(8)

其中每个si幅值介于0和1之间, 没有物理单位, 其大小可以用来表征第i个奇异向量中携带的深度信息, 矩阵(si, …, sQ)表示所有可能参与重建的奇异向量包含的深度信息.

当联合不同高度观测数据的奇异值进行反演重建时, DRP分析可以得到其理论的深度分辨率(图 2), 横坐标代表参与反演的奇异值数量, 纵坐标代表地下网格深度, 图像幅值强弱刻画其对应的深度分辨率信息大小.结果清晰表明随着更多观测面数据的增多, 反演可利用的奇异值数量改变, 参与数据重建的有效信息增多, 理论深度分辨率得到较为明显的增加.然而, 与奇异值谱分析的结果类似, 地面和两层空中测网即三个高度面数据理论上即可获得较高的深度分辨率, 再添加观测数据获得的深度分辨率增加较少且较弱, 如图 2.虽然图像上显示再增加高度观测异常可利用更多大序号(低幅值)的奇异值和奇异向量参与反演, 理论上提升一定分辨率, 但在实际情况下一方面增加的低幅值奇异值易受噪声干扰, 甚至被噪声淹没, 另一方面过多的航空观测数据会增加实际工作量和反演计算量, 成本可能远大于效果的提升, 得不偿失.

图 2 不同高度组合的深度分辨率图 Fig. 2 Depth resolution plot of different height combinations

为了验证空-地协同反演方法的有效性, 我们利用设计的模型(图 1)进行试验, 理论计算的各高度观测数据如图 3, 由于航空和地面观测系统不同, 在模型试验里我们不考虑噪声的干扰.分别利用地面观测数据(z=0)和地面与多高度数据组合的6种方式进行三维反演, 地下剖分方式与图 1一致, 为20×20×10个立方体, 反演结果的中心切片如图 4所示, 图中黑框代表真实模型的位置.图 4a为地面数据单独反演结果, 由于单一观测数据缺乏深度分辨能力, 反演结果底部分辨率较低.图 4bf为空地不同高度数据协同反演结果, 清晰地发现随空中测量数据的加入, 结果的底部上收, 深部分辨率得以提升.与图 4c相比, 图 4df改变的效果不大, 表明三个高度面即可获取较为真实的模型分布, 这一点与奇异值和DRP分析的结果一致, 验证了理论分析的有效性.图 5为各个反演结果的三维可视图, 可以提供更直观的效果展示.

图 3 双立方体模型不同高度观测重力异常 Fig. 3 Gravity anomaly of two cube model observed at different heights
图 4 不同高度组合数据反演结果x=1000 m处切片 Fig. 4 Inversion results of different height combinations sliced at x=1000 m
图 5 不同高度组合数据反演结果三维可视图(0.15 g·cm-3以下不显示) Fig. 5 3D perspective view of data inversion results using different height combinations (density contrast below 0.15 g·cm3 removed)
3 井中异常数据引入对反演结果影响

与航空和地面重磁测量一样, 井中测量获得重磁异常也是地下所有场源叠加产生的结果.对于埋藏较深的地质体, 其在地面或者航空观测面产生的异常相对较弱, 而井中测量能够更为靠近场源, 且获取的是随深度变化的观测数据列, 其对场源纵向上分布变化较为敏感.因而在地面重力数据反演中引入井中异常数据理论上能够对纵向上给予约束, 提高传统地面数据结果垂向分辨率和解释精度.由于井中重力测量获取的是垂向剖面数据, 其异常表征受井位限制, 且缺少水平方向性, 水平分辨率随距井距离增加而急剧下降, 使得井中测量控制区域有限.因此, 本文设计两种模型针对于异常旁井和穿异常井两种情况来讨论不同井位对反演结果的影响.

3.1 埋深相同的双模型

图 6a模拟两个同深度分布相距300 m的立方体, 周围分布五口旁井(A, B, C, D, E), 异常体空间三方向展布均为300 m.地面观测数据范围都为0到2000 m, xy方向点距为100 m, 这样地面观测点数为400个, 地面观测重力异常如图 6b.图 7为计算的五口井中观测重力异常, 每个井中观测20个点, 深度方向范围为50到1000 m, 垂向间距为50 m.

图 6 双立方体模型及其地面观测重力异常 Fig. 6 Two cube model and its observed surface gravity anomaly
图 7 井中观测重力异常 Fig. 7 Borehole observed gravity anomalies

本文将地下剖分为20×20×10=4000个块体, 分别利用地面重力异常数据、多个井中重力异常数据以及地面与井中重力数据进行反演, 反演模型的中心剖面(x=1000 m)如图 8所示, 黑框为模型体真实位置, 红色直线代表联合的井位.图 8a为地面重力数据反演结果, 预测模型的水平分辨率较高, 随深度增加分辨率逐渐下降.图 8b为ACD三口井中重力数据联合反演结果, 模型的垂向分辨率较好, 但密度值较小且在C井外侧产生虚假异常.图 8c为ABE三口井中重力数据联合反演结果, 可以发现结果分辨率在靠近井的一侧较高, 随离井距离急剧下降.图 8d为BED三口井中重力数据联合反演结果, 已经得不到可靠的场源分布范围.地面数据反演的结果具有较高的水平分辨率, 但缺少垂向分辨率;相反, 井中重力数据单独反演结果在靠近井附近垂向分辨率较好, 但由于井中重力数据没有水平方向性, 反演时会在无异常体的另一侧产生虚假异常.实验表明只使用多口井中重力数据联合反演的结果会存在较大的局限性, 准确性甚至不及原始地面数据反演结果.图 8eg为地面和各方位旁井重力数据协同反演结果, 既提高了地面反演结果的垂向分辨率又弥补井中重力反演缺少水平分辨率的不足.图 9给出了各种反演结果的三维可视化图, 为便于比较设置密度值小于0.2 g·cm-3不显示.旁井观测数据的引入可有效提高目标体的分辨能力, 这种提升随距井距离衰减的较为明显, 如仅利用旁井重力数据来获得较为准确的解释结果, 需在异常体附近布设多个钻孔.在下文第二个模型中, 我们考虑穿异常井数据对反演结果的改善程度.

图 8 地面与异常旁井数据协同反演结果切片图(x=1000 m处) Fig. 8 Slice map of joint inversion results of surface and anomaly-side borehole gravity data (at x=1000 m)
图 9 反演结果三维可视图(0.2 g·cm-3以下不显示) Fig. 9 3D perspective view of inversion results (density contrast below 0.2 g·cm-3 removed)
3.2 埋深不同的双模型

第二个模型采用不同埋深的两个棱柱体模型如图 10a, 两个异常体x, y, z三方向展布分别为500 m, 300 m, 300 m, 相距为300 m, 剩余密度为1 g·cm-3.区域内设计八口井, 其中存在穿地质体井(A、C、D、E)和旁井(B、P1、P2、P3), 钻孔A和C为完全穿透, D为贴近地质体上表面, E为进入地质体内部但未穿透的井.地面观测数据范围和观测点数与上个模型一致, 地面观测重力异常如图 9b.图 9c为计算的八口井中观测重力异常, 垂向观测间距为50 m, 除D、E井其他观测都为20个点, D井中有7个观测点, E井中有10个观测点.

图 10 深浅模型体及其地面观测重力异常 Fig. 10 Deep and shallow distributed model and its observed surface gravity anomaly

同样, 对地面重力数据和地面与井中观测数据组合进行三维反演, 对地下空间采用相同的剖分网格, 反演结果中心切片(x=1000 m)如图 12所示, 黑框为模型体真实位置, 白色直线代表协同反演联合的井位, 其中实线为穿异常井, 虚线为异常旁井.由于井中观测数据较少, 该模型中不考虑多口井数据单独反演, 均考虑地面与井中数据协同反演的情况.图 12a为单一地面重力数据反演结果, 结果对于浅部棱柱体有所表征, 深部棱柱体效果较差, 模型密度较低且底部分辨率较差.地面重力数据联合旁井P1、P2和P3重力数据协同反演时(图 12b), 结果垂向分辨率明显增加, 但较深棱柱体仍得不到表征, 与实际情况不符.当有穿过较深棱柱体的井中重力数据时(图 12c), 深部地质体成像结果得以增强, 但浅部地质体表征较弱, 说明穿异常井位能明显提升引起该异常地质体的分辨率.在此基础上增加旁井B可以进一步控制深部棱柱体成像结果的水平分布范围(图 12e), 图 12d12f的反演结果对比可以进一步证实这一说法, 联合叠加异常之间的B井重力数据进行协同反演使深浅分布的棱柱体水平分辨率更高.图 12fg为检验钻孔穿透深部地质体与否对该地质体反演结果的影响, 可以清晰地发现钻孔达到较深地质体上表面时反演结果即可分辨出其密度分布范围, 随着钻孔穿入越深, 反演密度体越接近真实模型, 当完全穿透地质体时反演结果分辨率最高(图 12f), 能够减少重力数据对于叠加场源解释的多解性.因此在实际资料处理时, 选择地面(或航空)重力异常极值上方或旁侧的井中重力数据进行协同反演能有助于获取更精细的目标体密度分布结果, 减少反演多解性, 提高解释准确性.图 13提供更为直观的反演结果三维可视化图.

图 11 各方位井中观测重力异常 Fig. 11 Borehole gravity anomaly observed in all positions
图 12 地面与穿异常井数据协同反演结果切片图(x=1000 m处) Fig. 12 Slice map of joint inversion results of surface and anomaly-through borehole gravity data (at x=1000 m)
图 13 反演结果三维可视图(0.2 g·cm-3以下不显示) Fig. 13 3D perspective view of inversion results (density contrast below 0.2 g·cm-3 removed)
4 实际案例

为了进一步验证空地井三维立体数据协同反演的有效性和适用性, 本文选取美国墨西哥湾地区的文顿盐丘实际勘探数据进行协同反演, 以恢复地下密度分布.文顿盐丘位于路易斯安那州西南部靠近德克萨斯州, 属穿刺性盐丘(Thompson and Eichelberger, 1928), 也是最早在盐丘区两翼找到石油的盐丘之一.盐丘主要由盐岩组成的盐核和一个较为明确的盐盖组成(Marfurt et al., 2004), 盐体的东西向分布较南北向分布宽, 盐盖主要成分从上至下由石灰岩逐渐变为石膏和硬石膏.Bell Geospace公司于2008年7月对该区进行了航空重力及张量梯度数据采集, 飞行高度平均85 m, 测线平均线距150 m.许多学者对该区进行了研究(Ennen and Hall, 2011; Oliveira Jr and Barbosa, 2013; Geng et al., 2014), 大部分结果表明该区重力异常主要由顶部高密度的盐盖引起, 平均剩余密度约为0.55 g·cm-3, 埋深范围约为200到600 m.在地面和航空实测数据中心截取5.4 km×4.5 km的子区域进行反演研究, 反演区域地下空间被离散为37×31×20=22940个小长方体, 每个长方体的大小为150 m×150 m×50 m.图 14ab分别表示反演所用地面和航空重力异常, xy方向分别指示北方向和东方向, 高异常指示盐盖分布位置.由于该地区尚无井中重力数据, 我们根据已有解释结果, 建立盐盖密度模型和四口不同方位钻孔(图 15), 井中观测间距为20 m, AC为穿盐盖井, BD为盐盖旁井.模拟计算井中观测重力异常并加入均值为0标准差为0.01 mGal的高斯随机噪声, 该噪声水平相当于现今用于矿产勘探的井中重力仪灵敏度的两倍(Seigel et al., 2010).

图 14 实测布格重力异常(a地面, b航空) Fig. 14 Real Bouguer gravity anomaly (a is surface, b is airborne)
图 15 盐盖模型和井中观测重力异常 Fig. 15 Salt cap rock model and observed borehole gravity anomalies

利用航空、地面和井中重力数据进行协同三维密度反演, 截取过盐盖中心的三方向切片结果(图 16), 其中16a—c为单一地面重力数据反演结果, 16d—f为航空-地面重力数据协同反演结果, 16g—i为航空-地面-井中重力数据协同反演结果.所有反演结果均指示推测出的高密度盐盖东西向最大长度是1100 m, 南北向最大长度是900 m, 并且盐盖顶面东南部埋深较浅约200 m, 向西北部迅速变深, 具有从东北—西南向的构造延伸方向, 这一结果与Coker等(2007)的构造解释吻合, 验证了反演结果准确性.单一地面重力数据反演结果(图 16ac)能比较准确地指示盐盖水平位置和分布, 但地质体深部分辨率略低, 反演的剩余密度比真实密度略小.通过引入航空重力数据进行协同反演, 地质体信息得到增强(图 16df), 反演结果剩余密度得到提高.当综合航空-地面-井中重力数据进行协同反演时, 基于井中重力数据的高垂向分辨能力, 反演结果的垂向分辨率明显提高(图 16gi), 地质体上下边界清晰, 尤其是深部无冗余构造, 且剩余密度与真实密度更为接近, 进一步降低反演的多解性, 提高解释精度.图 17展示了空地井立体重力数据协同反演预测数据与实际观测数据的拟合情况, 经过盐盖中心的两个剖面和井中的预测数据(图中红色圈)与实际地面和井中观测重力异常(图中蓝色星号)基本吻合, 进一步说明了协同反演的正确性.因此, 在实际地区工作中, 在开展区域地面重力测量同时应该适当开展航空重力测量, 最好获得两个飞行高度的空中测网数据, 空地重力协同反演可减少地面数据反演多解性.条件允许的情况下, 在航空和地面数据圈定的异常靶区中选取异常内部及异常旁侧的钻孔进行井中重力测量, 联合穿异常井和异常旁井数据进行空地井立体数据协同反演, 能够对反演结果在井周围及垂向分布上给予约束, 进一步提高分辨率和深部分辨能力.

图 16 反演结果三方向切片图(a, d, g为y=2700 m处切片;b, e, h为x=2550 m处切片;c, f, i为z=300 m处切片) Fig. 16 Three-direction slices of inversion results (a, d, g is sliced at y=2700 m; b, e, h is sliced at x=2550 m; c, f, i is sliced at z=300 m)
图 17 反演预测数据与实测数据拟合图(蓝色星号为实际观测的, 红色圈为反演预测的) Fig. 17 Fitting results of inversion predicted and real observed data (blue asterisk for observed, red circle for predicted)
5 结论

传统地面重力异常的反演和解释存在垂向分辨率不足的问题, 给实际资料解释带来多解性和不确定性.针对该问题, 我们利用空地井立体数据进行协同反演, 在不引入额外先验地质信息约束情况下仍能获得较高精度和分辨率的结果.理论和模型试验表明, 空地协同反演能够提升地质体深部分辨率, 地面联合两层空中测网数据即可获得较高分辨率的结果.井中重力数据表征异常体随深度变化的异常, 但由于其缺乏水平方向性且井位控制范围有限, 仅利用单独或几个钻孔内观测数据反演在实际资料解释中仍会存在较强的多解性.因而结合地面和航空资料, 进行空地井协同反演能够提高结果垂向分辨率和解释精度.穿异常体井作用大于异常体旁井, 穿异常体井配合异常体旁井可以更好地实现场源的精细定位.此外, 美国文顿盐丘实际数据处理结果验证空地井协同反演方法的正确性及适用性.因此, 在实际地区工作中, 可以适当开展航空及井中重力测量技术, 综合利用多维数据的横向和纵向变化特征来实现空地、地井或空地井立体数据的三维协同反演, 提高结果精度和分辨率, 尤其是提高深部分辨能力, 帮助解释人员获得更为可靠的解释模型.

References
Chasseriau P, Chouteau M. 2003. 3D gravity inversion using a model of parameter covariance. Journal of Applied Geophysics, 52(1): 59-74. DOI:10.1016/S0926-9851(02)00240-9
Coker M O, Bhattacharya J P, Marfurt K J. 2007. Fracture patterns within mudstones on the flanks of a salt dome:Syneresis or slumping?. Gulf Coast Association of Geological Societies Transactions, 57: 125-137.
Ennen C, Hall S. 2011. Structural mapping of the Vinton salt dome, Louisiana, using gravity gradiometry data.//81st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 30(1): 830-835.
Fedi M, Hansen P C, Paoletti V. 2005. Analysis of depth resolution in potential-field inversion. Geophysics, 70(6): A1-A11. DOI:10.1190/1.2122408
Geng M X, Huang D N, Yang Q J, et al. 2014. 3D inversion of airborne gravity-gradiometry data using cokriging. Geophysics, 70(6): G37-G47.
Geng M X, Yang Q J, Huang D N. 2015. 3D joint inversion of gravity-gradient and borehole gravity data. Exploration Geophysics, 48(2): 151-165.
Guo Z H, Xiong S Q, Zhou J X, et al. 2008. The research on quality evaluation method of test repeat lines in airborne gravity survey. Chinese Journal of Geophysics (in Chinese), 51(5): 1538-1543.
Haáz I B. 1953. Relationship between the potential of the attraction of the mass contained in a finite rectangular prism and its first and second derivatives. Geophysical Transactions, 2: 57-66.
Han M R, Wan L, Zhdanov M. 2018. Joint iterative migration of surface and borehole gravity gradiometry data.//88th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1399-1403.
Hansen P C. 1998. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Denmark: Society for Industrial and Applied Mathematics.
LaFehr T R. 1983. Rock density from borehole gravity surveys. Geophysics, 48(3): 341-356.
Lautzenhiser T V. 1983. Amoco/Lacoste and Romberg automated borehole gravity meter.//53rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 20-21.
Li Y G, Oldenburg D W. 1996. 3-D inversion of magnetic data. Geophysics, 61(2): 394-408. DOI:10.1190/1.1443968
Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data. Geophysics, 63(1): 109-119. DOI:10.1190/1.1444302
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
MacQueen J D. 2007. High-resolution density from borehole gravity data.//77th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 741-744.
Marfurt K J, Zhou H W, Sullivan E C. 2004. Development and Calibration of New 3-D Vector VSP Imaging Technology: Vinton Salt Dome. LA. Houston, TX: University of Houston, 169.
Mosher C R W, Farquharson C G. 2013. Minimum-structure borehole gravity inversion for mineral exploration:A synthetic modeling study. Geophysics, 78(2): G25-G39.
Nind C J M, Seigel H O, Chouteau M, et al. 2007. Development of a borehole gravimeter for mining applications. First Break, 25(7): 71-77.
Oliveira V C Jr, Barbosa V C F. 2013. 3-D radial gravity gradient inversion. Geophysical Journal International, 195(2): 883-902. DOI:10.1093/gji/ggt307
Paoletti V, Hansen P C, Hansen M F, et al. 2014. A computationally efficient tool for assessing the depth resolution in large-scale potential-field inversion. Geophysics, 79(4): A33-A38. DOI:10.1190/geo2014-0017.1
Pilkington M. 2012. Analysis of gravity gradiometer inverse problems using optimal design measures. Geophysics, 77(2): G25-G31.
Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596
Portniaguine O, Zhdanov M S. 2002. 3-D magnetic inversion with data compression and image focusing. Geophysics, 67(5): 1532-1541. DOI:10.1190/1.1512749
Qin P B, Huang D N. 2016. Integrated gravity and gravity gradient data focusing inversion. Chinese Journal of Geophysics (in Chinese), 59(6): 2203-2224. DOI:10.6038/cjg20160624
Rim H, Li Y G. 2015. Advantages of borehole vector gravity in density imaging. Geophysics, 80(1): G1-G13.
Seigel H O, Nind C J M, Milanovic A, et al. 2010. Results from the initial field trials of a borehole gravity meter for mining and geotechnical applications.//80th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 92-96.
Shamsipour P, Marcotte D, Chouteau M, et al. 2010. 3D stochastic inversion of gravity data using cokriging and cosimulation. Geophysics, 75(1): I1-I10.
Sun J J, Li Y G. 2010. Inversion of surface and borehole gravity with thresholding and density constraints.//80th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 1798-1803.
Sun Z M, Zhai Z H, Wu F M, et al. 2015. Algorithm comparison for strapdown airborne gravimetry. Chinese Journal of Geophysics (in Chinese), 58(5): 1547-1554. DOI:10.6038/cjg20150508
Thompson S A, Eichelberger O H. 1928. Vinton salt dome, Calcasieu parish, Louisiana. AAPG Bulletin, 12(4): 385-394.
Wang T H, Huang D N, Ma G Q, et al. 2017. Improved preconditioned conjugate gradient algorithm and application in 3D inversion of gravity-gradiometry data. Applied Geophysics, 14(2): 301-313. DOI:10.1007/s11770-017-0625-x
Wang Y C. 2014. Study of equidimension focusing inversion of gravity data[Master's thesis] (in Chinese). Beijing: China University of Geosciences.
Xiong S Q. 2009. The present situation and development of airborne gravity and magnetic survey techniques in China. Progress in Geophysics (in Chinese), 24(1): 113-117.
Zhang Y M, Zhang G B. 2009. Research on the equivalent dimensions of 2-D gravity. Chinese Journal of Engineering Geophysics (in Chinese), 6(2): 137-142.
Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. Salt Lake City: Elsevier Science.
Zhdanov M S, Ellis R, Mukherjee S. 2004. Three-dimensional regularized focusing inversion of gravity gradient tensor component data. Geophysics, 69(4): 925-937. DOI:10.1190/1.1778236
Zhdanov M S. 2009. New advances in regularized inversion of gravity and electromagnetic data. Geophysical Prospecting, 57(4): 463-478. DOI:10.1111/j.1365-2478.2008.00763.x
郭志宏, 熊盛青, 周坚鑫, 等. 2008. 航空重力重复线测试数据质量评价方法研究. 地球物理学报, 51(5): 1538-1543. DOI:10.3321/j.issn:0001-5733.2008.05.028
秦朋波, 黄大年. 2016. 重力和重力梯度数据联合聚焦反演方法. 地球物理学报, 59(6): 2203-2224. DOI:10.6038/cjg20160624
孙中苗, 翟振和, 吴富梅, 等. 2015. 捷联式航空重力测量算法比较. 地球物理学报, 58(5): 1547-1554. DOI:10.6038/cjg20150508
王逸宸. 2014.重力等维聚焦反演研究[硕士论文].北京: 中国地质大学(北京).
熊盛青. 2009. 我国航空重磁勘探技术现状与发展趋势. 地球物理学进展, 24(1): 113-117.
张永明, 张贵宾. 2009. 二维重力"等维"反演研究. 工程地球物理学报, 6(2): 137-142. DOI:10.3969/j.issn.1672-7940.2009.02.002