重力及其垂直梯度交叉约束联合密度结构反演方法

马国庆, 李欣蔚, 王泰涵*, 熊盛青, 高桐

1 吉林大学地球探测科学与技术学院, 长春 130026 2 中国自然资源航空物探遥感中心, 北京 100083

0 引言

重力梯度数据具有较高水平分辨率,相对原始重力场可凸显较浅地质体.重力梯度测量可通过不同高度测量重力值的方式来获得,其与重力联合测量可以探测不同深度的空洞、人防工程(Cordell, 1979; Schlinger, 1990; Bate, 2005; Butler, 2012). 通过重力异常反演地下密度分布可以圈定场源体的分布,主要的反演方法包括光滑约束反演、粗糙约束反演、聚焦反演、协克里金反演、模糊C均值聚类反演等(Li and Oldenburg, 1998; Portniaguine and Zhdanov, 1999; Zhdanov et al., 2004; Pilkington, 2009; Geng et al., 2014; Sun and Li, 2015; 王泰涵等, 2020; 李芳等, 2021).

重力和梯度的联合反演方法现今主要采用将不同分量数据组合为一个矩阵,通过数据约束来提升反演结果的精度(Tikhonov and Arsenin, 1977; Qin et al., 2016; Wang et al., 2022).Wu等(2013)提出加权最小二乘估计方法,通过自适应权值实现同源数据的联合.Capriotti和Li(2014)基于灵敏度矩阵平衡核函数衰减速率,并定义了一种重力和张量梯度数据联合的公式.Qin等(2016)以非线性牛顿法和最小梯度支撑进行重力和重力梯度反演计算,获得了较清晰的异常边界分布特征.高秀鹤等(2019)利用密度阈值不断更新协方差矩阵来实现重力和重力梯度数据协克里金反演.重力与重力梯度联合反演还可以采用交叉梯度方法,其是从结构上对不同物性参数进行耦合的联合反演方法,主要用于不同物性源地球物理数据的联合反演(Gallardo and Meju,2003; Fregoso and Gallardo, 2009; Zhou et al., 2015).重力和重力梯度数据是同源数据,但由于不同参量数据对于深浅场源反演效果的差异,为此采用交叉梯度的结构约束特征可提升反演的分辨率和精度(Zhao et al., 2018; Entezar-Saadat et al., 2020).

为了充分利用重力与重力梯度对地下密度结构的敏感性差异，我们提出重力及其垂直梯度交叉约束联合密度结构反演方法,以重力与重力梯度单独反演结果作为权函数来约束交叉梯度联合反演过程,有效恢复地下不同深度场源的密度结构分布特征.最后通过模型实验和实际数据证明该方法较常规数据联合密度反演方法在反演分辨率和精度的优势，以及在实际应用中的适用性.

1 基本原理

现今重力数据(Vz)和垂直梯度数据(Vzz)密度反演通常将不同数据放置于同一矩阵,采用Tikhonov正则化方法进行反演,其目标函数的表达式为(Tikhonov and Arsenin, 1977; Qin et al., 2016)：

(1)

式中为Φd数据拟合项；Φm为模型约束项；α为正则化参数，用于权衡数据拟合项和模型约束项对正则化方程的影响；Wd为数据加权矩阵；F为重力及其垂直梯度的组合核函数矩阵[Fz,Fzz]T；m为物性(密度)参数；d为重力及其垂直梯度观测数据的组合矩阵[Vz,Vzz]T；Wz为深度加权矩阵,可以有效抵消核矩阵随深度的固有衰减,以降低趋肤效应的影响(Li and Oldenburg, 2000).

采用不同大小、不同埋深的异常体A和异常体B来分析重力和重力梯度联合的Tikhonov正则化反演方法效果,其具体参数如表1所示,分布如图1a所示.图1b为重力异常更好地体现较深场源信息,但两个异常体并未分开.图1c为垂直梯度异常更好地反映较浅场源的信息,并且两个异常体的边界较清晰.图1d、e分别为利用Tikhonov正则化方法反演得到的重力和重力梯度的密度反演结果在y=1000 m位置的切片图,可以看出重力结果在深部的分辨率高于重力梯度,重力梯度对于浅部地质目标有更高的反演精度.重力及其垂直梯度异常虽反应地下同一密度结构,但是由于数据敏感度差异，使恢复密度结构易于集中在变化较大区域,模型分辨率产生了差异性.图1f为重力和重力梯度组合矩阵的Tikhonov正则化方法反演结果,对于浅部地质体的反演效果有一定的提升,但是对于深部地质体依旧较差.

图1 两个棱柱体模型的正演异常及不同反演方法结果在y=1000 m处剖面图(a) 两个棱柱体组合模型; (b) 重力异常; (c) 重力垂直梯度异常; (d) Tikhonov正则化方法Vz数据反演; (e) Tikhonov正则化方法Vzz数据反演; (f) Tikhonov正则化方法[Vz, Vzz]T数据反演; (g) 交叉梯度Vz数据反演; (h) 交叉梯度Vzz数据反演; (i) 交叉梯度Vzz数据反演结果三维透视图.Fig.1 Forward response and Section at y=1000 m of inversion results based on different methods of two-prism model(a) The model combined by two prisms; (b) Gravity anomaly; (c) Gravity vertical gradient anomaly; (d) Tikhonov regularization inversion result of Vz data; (e) Tikhonov regularization inversion result of Vzz data; (f) Tikhonov regularization inversion result of [Vz, Vzz]T; (g) Cross-gradient inversion result of Vz data; (h) Cross-gradient inversion result of Vzz data; (i) 3D perspective view of cross-gradient inversion result of Vzz data.

表1 两个棱柱体模型Table 1 Parameters oftwo-prism model

基于不同重力参量数据对应同源的密度结构,但是不同数据反演结果存在相似性和差异性,采用交叉梯度引入的结构约束可有效提升结果的精度和分辨率.交叉梯度函数(Gallardo and Meju, 2003)的表达式为：

其中，m1和m2分别为不同的物性参数，当两个梯度方向平行或者不存在时t为0.其在x、y、z三个方向上的分量可以表示为：

(3)

采用交叉梯度的方法进行重力数据和重力梯度联合反演,其目标函数为：

(5)

式中，λz和λzz为交叉约束项系数.式中参数的下角标z和zz分别代表重力和重力梯度数据相关参数.Wz1和Wz2分别为重力及其梯度数据对应的深度加权矩阵.采用图1b、c所示的重力及其梯度数据进行交叉梯度反演，得到图1g、h的两个反演结果，深部和浅部分辨率有提高,但与实际模型依旧存在较大偏差.

为了提升重力和重力梯度联合反演的分辨率来获取不同深度场源的分布,本文提出重力及其垂直梯度交叉约束联合密度结构反演方法,其首先对重力异常与重力梯度进行单独Tikhonov正则化反演,然后将各自反演结果形成的归一化权函数进行交叉互换使用,并采用交叉梯度进行结构约束来实现重力与其梯度的联合反演,以期达到不同深度场源分辨率的提升,最终取双方结果的平均值作为最终的密度分布,具体的计算流程如图2所示.

图2 交叉约束联合密度结构反演技术流程Fig.2 The flow of the cross-constrained density joint inversion based on gravity and vertical gradient data

图2中，Tikhonov正则化方法重力数据反演结果mz1和重力垂直梯度数据反演结果mzz1形成的交叉约束权函数Wρz和Wρzz的表达式为：

(6)

其中(mz1)max、(mzz1)max为矩阵中的最大的元素；(mz1)min、(mzz1)min为矩阵中最小的元素.为了避免产生奇点,令εz和εzz分别等于0.1倍的(mz1)min和0.1倍的(mzz1)min.μ是一个控制交叉约束权在反演中约束影响的常数,通常取值为1.不同参量结果互为权函数的交叉梯度反演表达式为：

最后，取交叉约束权联合密度反演结果mz2和mzz2的平均值作为最终密度分布mFin.

采用表1中的模型来试验交叉约束联合密度结构反演方法的效果,图3a、b为采用该方法进行反演获得的重力和重力梯度的反演结果,可以看出异常的位置准确,水平分辨率较高.图3c为重力和重力梯度的反演结果取平均值融合后得到的最终结果,可以看出相对放置同一矩阵反演结果的精度大,且密度值更接近理论值.图3d为密度大于0.65 g·cm-3的三维图,其中蓝色边框为实际模型范围,本文方法结果形态较前文的其他反演方法更符合实际密度分布,且密度幅值更接近模型的真实密度.为了验证所提出交叉约束联合密度结构反演方法的准确性,将不同反演方法结果在y=1000 m处的重力与正演重力异常进行对比,从图4可以看出在异常变化平缓区域全都拟合较好,而异常变化率较大的位置,本文方法的反演结果的拟合差更小.

图3 交叉约束联合密度结构反演结果(a) Vz数据反演结果在y=1000 m处剖面图; (b) Vzz数据反演结果在y=1000 m处剖面图; (c) 反演最终结果; (d) 密度大于0.65 g·cm-3反演结果三维透视图.Fig.3 Inversion results of the cross-constrained joint inversion method(a) Section at y=1000 m of inversion result of Vz data; (b) Section at y=1000 m of inversion result of Vzz data; (c) Section at y=1000 m of final inversion result; (d) 3D perspective view of inversion result with density larger than 0.65 g·cm-3.

图4 不同方法反演结果在y=1000m处的与观测数据拟合曲线Fig.4 The fitting curve map of inversion results of different methodsand observed data at section of y=1000 m

2 模型试验

为了检验本文方法的抗噪能力,对图5a所示的模型的重力和梯度数据加入最大值为10%的高斯噪声,含噪声的数据正演异常如图5b、c所示.表2统计了重力数据进行正则化反演、重力及其梯度数据正则化联合反演、重力及梯度数据交叉梯度联合反演以及本文方法的预测数据与观测数据的拟合差进行对比，可以看出本文方法较其他反演方法拟合差更小,证明了本文方法抗噪性良好.

图5 含噪数据正演异常和不同方法反演结果在y=1000 m处剖面图(a) 两个棱柱体组合模型; (b) 重力异常; (c) 重力垂直梯度异常; (d) Tikhonov正则化方法Vz数据反演结果; (e) Tikhonov正则化方法Vzz数据反演结果; (f) Tikhonov正则化方法[Vz, Vzz]T反演结果; (g) 交叉梯度Vz数据反演结果; (h) 交叉梯度Vzz数据反演结果; (i) 交叉约束Vz数据反演结果; (j) 交叉约束Vzz数据反演结果; (k) 交叉约束方法反演最终结果; (l) 交叉约束方法结果三维透视图.Fig.5 Forward anomaly and sections at y=1000 m of inversion results of different methods of data with noise(a) The model combined by two prisms; (b) Gravity anomaly; (c) Gravity vertical gradient anomaly; (d) Tikhonov regularization inversion result of Vz data; (e) Tikhonov regularization inversion result of Vzz data; (f) Tikhonov regularization inversion result of [Vz, Vzz]T; (g) Cross-gradient inversion result of Vz data; (h) Cross-gradient inversion result of Vzz data; (i) Cross-constrained inversion result of Vz data; (j) Cross-constrained inversion result of Vzz data; (k) Final inversion result of cross-constrained method; (l) 3D perspective view of cross-constrained inversion result.

表2 含噪数据不同反演方法拟合差Table 2 The misfit of the different inversion methods for data with noise

对比重力异常图5b和重力梯度异常图5c,重力梯度数据水平分辨率高于重力数据.重力数据较梯度数据对于深部异常更为敏感,重力梯度数据对于浅部异常更敏感.含噪声模型实验的Tikhonov正则化反演方法梯度数据结果图5e浅部异常的分辨率和水平分辨率优于重力数据反演结果图5d,图5d深部异常的分辨率优于图5e，验证了重力和梯度数据的数据特征.并且梯度数据对于高频信号较为敏感,噪声属于高频信号,可以看到有梯度数据参与反演的Tikhonov正则化方法Vzz数据反演结果图5e、Tikhonov正则化方法[Vz,Vzz]T反演结果图5f和交叉梯度Vzz数据反演结果图5h在浅部都一定程度上受噪声干扰.采用交叉约束联合密度结构反演方法,可以能够降低噪声的影响,使反演结果更集中于密度高值分布的区域.交叉约束联合密度结构反演方法的结果图5i、j在两个异常体定位和分辨率方面明显优于图5d—f所示的Tikhonov正则化方法反演结果以及图5g、h交叉梯度的反演结果.图5l为密度值大于0.65 g·cm-3的部分的三维密度结果,比较符合真实密度分布.

采用不同密度的异常体组合模型来验证交叉约束密度结构反演方法对密度差异区间较大的复杂地质体的适用性,模型参数见表3,图6a为异常体的密度真实分布.

图6 三个棱柱体模型及其正演异常(a) 三个棱柱体组成的模型; (b) 重力异常; (c) 重力垂直梯度.Fig.6 Three-prism model and forward anomaly(a) The model combined by three prisms; (b) Gravity anomaly; (c) Gravity vertical gradient anomaly.

表3 三个棱柱体模型参数Table 3 Parameters of three-prism model

重力异常图6b和重力垂直梯度异常6c中异常体③的分辨率较其他两个异常体高,梯度数据原始异常的异常体②和③的边界较重力异常的边界更为清晰.取x=1400 m,y=1000 m和z=400 m处的反演结果图对比不同方法的反演结果.由于Tikhonov和交叉梯度法的反演结果幅值较低,为了方便对比不同反演方法的结果,将两种反演方法的结果统一色标柱以便对比.

图7d—f梯度数据的Tikhonov正则化反演结果较重力Tikhonov正则化反演结果图7a—c整体分辨率较高.Tikhonov正则化方法[Vz,Vzz]T反演结果图7h以及交叉梯度方法的图7k和图7n中异常体②水平方向密度分布较Tikhonov正则化法的反演结果图7b、e更为收敛.Tikhonov正则化方法和交叉梯度方法的在z=400 m 反演结果图7c、f、i、l、o在水平方向上的分辨率都表现良好,三个异常体的边界都比较明晰.由本文方法的水平方向剖面图7r可以看出密度高值集中于异常体中心,其垂向剖面图图7p、q较Tikhonov正则化方法和交叉梯度方法的垂向剖面图7a、b、d、e、g、h、j、k、m、n密度分布更符合实际情况,说明本文方法在提升纵向分辨率上优于Tikhonov正则化方法和交叉梯度方法.由图7r本文方法的反演结果得到的三个异常体中心的密度值分别为0.11 g·cm-3、0.46 g·cm-3、0.73g·cm-3，幅值较接近理论模型的值,证明了本文方法在反演情况复杂的异常体时,对异常体位置定位准确,密度分布更加接近实际情况,三个异常体的边界都很清晰,并且有较高的纵向分辨率.

图7 不同反演方法在x=1400 m, y=1000 m, z=400 m反演结果图(a),(b)及(c)为Tikhonov正则化方法Vz数据反演; (d), (e)及(f) Tikhonov正则化方法Vzz 数据反演; (g), (h)及(i) Tikhonov正则化方法[Vz, Vzz]T反演; (j), (k)及(l) 交叉梯度Vz数据反演; (m),(n)及(o) 交叉梯度Vzz数据反演; (p), (q)及(r) 交叉约束方法反演.Fig.7 Slices at x=1400 m, y=1000 m and z=400 m of inversion results of diffrent methods(a), (b) and (c) Tikhonov regularization inversion of Vz data; (d), (e) and (f) Tikhonov regularization inversion of Vzz data; (g), (h) and (i) Tikhonov regularization inversion of [Vz, Vzz]T; (j), (k) and (l) Cross-gradient inversion of Vz data; (m), (n) and (o) Cross-gradient inversion of Vzz data; (p), (q) and (r) Cross-constrained method.

3 实际应用

3.1 辽源矿区重力及垂直梯度联合反演

辽源煤矿位于长白山余脉与松嫩平原交界,矿产资源丰富.由于长期的无序开采,地面塌陷给人民生活和城市建设带来安全隐患 (武欣等,2022).图8为辽源矿区图,图中标注了开采的矿点、已废弃矿点以及地表塌陷区域,黑色线框内区域为研究区域,面积为3 km×2 km.我们将目标区域划分为30×20×10个棱柱体单元,每个单元的大小为100 m×100 m×30 m.通过重力及其梯度交叉约束联合密度结构反演方法,确定采空区的位置,以圈定测区内采空区范围.

图8 辽源矿区Fig.8 Liaoyuan mining area

采空区引起的重力异常相对较小,为了寻找采空区采用CG-5重力仪对测区进行高精度的微重力测量(图9a、c).仪器读数分辨率1 μGal,测量的均方误差为0.087 mGal.对实测重力异常去除背景场,并对重力数据和垂直梯度数据做了地形改正(陈善, 1986; 叶周润等, 2011),得到的剩余异常和重力梯度异常如图9b、d所示.围岩密度为2.1 g·cm-3,采空区密度为1.0 g·cm-3,采空区引起的异常密度差接近-1.1 g·cm-3.

图9 野外重力测量及重力异常(a) 重力梯度测量; (b) 重力异常; (c) CG-5重力仪; (d) 重力垂直梯度.Fig.9 Gravity measurement of field exploration and gravity anomaly(a) Gravity gradiometry; (b) Gravity anomaly; (c) CG-5 Gravimeter; (d) Gravity vertical gradient anomaly.

重力及其梯度交叉约束联合密度结构反演方法结果如图10a所示,采空区密度为负值,通过本文方法计算得到的采空区密度最低值为-0.88 g·cm-3与推测的剩余密度值较为符合.在反演结果150 m的切片结果上根据异常边界圈定了Ⅰ-Ⅴ号采空区,Ⅵ号区域附近为塌陷区,所以也将其认定为一个采空区(图10a).根据图8辽源矿区图,将研究区域内的已知矿点画在150 m的反演结果上,10个矿点的位置,都位于圈定的采空区内,表明本文方法的反演结果与以往已知矿点位置对应良好.采空区的深度范围为0～250 m,平均深度200 m.图10c为测线LL′测线的高密度电阻率剖面,推测在该剖面上采空区深度大约为170 m.图10b为反演结果在测线同一位置的剖面图,该剖面三处区域的形状及位置与高密度电阻率剖面上的范围有较好的位置对应,验证了反演结果的准确性.

图10 反演密度分布结果与高密度电阻率法结果对比(a) 本文交叉约束方法反演结果; (b) 交叉约束反演结果在LL′测线的剖面图; (c) LL′测线的高密度电阻率法剖面图.Fig.10 Comparison of inversion results of density distribution and high density resistivity method(a) Inversion result of cross-constrained method; (b) Slice of Line LL′ of the cross-constrained method; (c) Slice of Line LL′ of the high-density resistivity method.

3.2 长白山火山区重力及垂直梯度联合反演

长白山火山区是新生代多成因复合火山,地球物理勘探发现长白山天池火山下方有岩浆囊的存在,火山有再次喷发的风险,因此研究长白山火山岩浆囊的位置和范围对火山的监测具有重大价值 (郭履灿等, 1996; 吕政等, 2007).

本文选择430 km×430 km 的长白山火山区卫星重力数据,去除沉积层和莫霍面起伏产生的影响,得到图11a剩余重力异常(Blakely,1995).对重力异常求导得到图11b垂直梯度异常来进行重力及其垂直梯度联合反演.

图11 长白山研究区域重力异常(a) 重力异常; (b) 重力垂直梯度异常.Fig.11 Gravity anomaly of Changbai Mountain survey area(a) Gravity anomaly; (b) Gravity vertical gradient anomaly.

将地下网格剖分为43×43×10个单元,每个单元大小为10 km×10 km×3 km.岩浆囊在高温熔融状态较围岩呈低密度,重力异常显示为低值.重力及其梯度交叉约束联合密度结构反演方法的三维密度结果如图12a.图12b、c为穿过岩浆囊的两条剖面.根据通过岩浆囊中心的L1和L2剖面图,岩浆囊呈拱形,位于北向250～370 km,东向190～260 km,呈东北-西南方向分布.本文方法估计的岩浆囊的中心埋深约15 km,深度范围约为7～20 km.汤吉等(2001)根据大地电磁测深反演剖面结果显示10～30 km的深度范围内存在低阻体,认为其可能为活动的岩浆囊; 管彦武等(2020)利用重力剖面数据及地质信息进行人机交互建立的密度模型推断深度7～15 km处存在岩浆囊,本文推断的岩浆囊深度整体范围及中心埋深与这些前人研究比较吻合,进一步验证了本文方法的适用性.

图12 长白山区域反演结果(a) 交叉约束方法结果三维透视图; (b) 反演结果在L1处的剖面图; (c) 反演结果在L2处的剖面图.Fig.12 Inversion result of Changbai Moutain area(a) 3D perspective view of the cross-constrained inversion result; (b) Slice of Line L1 of the inversion result; (c) Slice of Line L2 of the inversion result.

4 结论

本文提出了重力及其垂直梯度联合的交叉约束联合密度结构反演方法,通过建立交叉权函数的进行交叉梯度联合反演,有效提升了结果的分辨率和精度.通过模型试验证明本文方法能高分辨率地获取不同深度场源的密度变化特征,尤其是相对以往方法对于较深场源的反演精度大大提升.此外,本文方法的抗噪性较强,对复杂地质体适用性良好.将本文方法应用于辽源矿区的实际数据反演,圈定了6个采空区,其位置与电法资料以及已知矿点完全吻合,且描绘出了大致范围.通过本文方法对长白山火山区的卫星数据的应用,我们估计了岩浆囊的位置以及范围,其结果与前人研究结果基本一致.实际数据应用证明了本文方法的实用性,为区域规划建设及火山灾害监测提供了准确的依据.

致谢感谢审稿专家提出的修改意见.