2. 中色杰泰地球物理科技(北京)有限公司, 北京 100012;
3. 河北工程大学地球科学与工程学院, 河北邯郸 056038
2. Sino-GT Geophysical Technology(Beijing) Company Limited, Beijing 100012, China;
3. College of Geosciences and Engineering, Hebei University of Engineering, Handan Hebei 056038, China
基于重力数据的地质体解释主要包括其空间位置及物性分布属性.边界识别、深度反演等方法可快速的获取地质体的空间位置参数信息(Fedi and Abbas, 2013;周帅等,2016;张琦等,2018),而物性反演方法通过求解物性参数与观测异常间方程的最优化问题实现场源的密度或磁性参数的计算.重力数据物性反演面临严重的多解性问题,为了获得唯一、可靠的物性参数反演结果,通常采用正则化方法加入模型约束函数降低反演多解性,常用的约束方法有最小质量约束、最小体积约束、最大光滑约束等(Last and Kubik, 1983;Guillen and Menichetti, 1984;Li and Oldenburg, 1996).正则化方法可以获得数学意义上的最优解,但获得的解释结果可能与现有的地质或其他地球物理资料差别较大.加入更多的先验信息可有效地提高反演模型的可靠性.Silva等(2001)讨论了现有物性反演方法中的约束条件,并建议约束条件的选择应根据实际的地质问题确定.Lelièvre和Oldenburg(2009)在Li和Oldenburg(1996)最大光滑模型约束函数的基础上将地质构造走向信息约束到重磁数据的物性参数反演中.为了结合地质目标的先验深度信息,Commer(2011)提出基于地质体深度参数的加权函数,提高了反演的深度分辨率.Zhang等(2015)通过拉格朗日因子和松弛变量将已有地质构造信息加入到目标函数中,更好地揭示了上地幔的密度分布情况.Carter-McAuslan等(2015)对重力数据和地震初至走时数据联合反演,并采用模糊C-均值聚类法获得了更加准确的密度模型和速度模型.郭曼等(2018)和殷长春等(2018)通过大地电磁数据与重力数据的联合反演降低了位场数据反演的多解性.Giraud等(2017)结合了已有地质和岩石物理信息进行重磁联合物性参数成像,有效的降低了反演的不确定性.
通过加入已有地质或其他地球物理的先验信息,物性反演有效地改善了多解性严重的问题,并提高了反演结果的可靠性.对于先验信息缺乏的情况,充分挖掘数据本身包含的约束信息并用于反演计算中,可获得更可靠的物性分布结果.Cella和Fedi(2012)指出物性反演中的深度加权函数的衰减系数与构造指数有关,而后Ialongo等(2014)将欧拉反褶积获得的构造指数信息作为衰减系数参与到物性反演中,在深度方向上获得了更准确的物性分布.Paoletti等(2013)将多尺度快速成像法计算的地质体空间位置参数信息约束到光滑反演中,获得了更高分辨率的物性空间分布结果.Fregoso等(2015)在重磁交叉梯度联合反演计算过程中结合了欧拉反褶积输出的深度参数.Lo Re等(2016)利用快速几何参数反演方法获取的空间位置参数约束到微重数据的密度成像中,有效地反演出深浅断层构造的密度分布结构.Sun和Chen(2016)将相关系数成像结果形成约束条件参与到反演计算过程中,在磁法数据的解释中获得更可靠的反演结果.
本文通过基于重力梯度数据的几何参数反演方法获取地质目标的空间位置信息,并将获取的空间位置先验信息转化为优化约束条件参与到物性反演中,通过对航空重力梯度张量数据进行处理解释,进一步提高了重力梯度张量数据物性反演的可靠性和精度.
1 原理与方法 1.1 正演计算在重力数据反演中,一般采用长方体单元作为正演计算的基本单元,每个单元内密度均匀且空间位置已知,进而建立剩余密度与原始异常数据的定量计算关系式,通过求解各个单元内的物性参数完成地质体三维物性反演计算.
剩余密度为ρ的地质体在其外部空间引起的引力位为
(1) |
其中,r为观测点(x, y, z)到场源点(ξ, η, ζ)的距离,G是万有引力常量.(ξ, η, ζ)为场源内任何一点的坐标位置,观测点坐标为(x, y, z).对引力位求垂向一阶导数获得原始重力异常,对x, y, z三个方向求解二阶导数获得梯度张量矩阵,公式为
(2) |
在地球外部引力位满足拉普拉斯方程Vxx+Vyy+Vzz=0,并且重力场具有无旋性,因此,重力梯度张量矩阵中只有五个独立分量.
在反演计算时需要将地下空间进行离散化处理,本文将地下空间划分为大小相等的长方体单元,剖分长方体基本单元的边长与点距相等.对于观测平面N个观测点获得的重力梯度分量的异常值为
(3) |
其中,ρj表示第j个基本单元的剩余密度,Gij代表梯度分量中第j个单元的几何架构.i观测点异常gαβi是地下M个基本单元引起的异常叠加之和,公式为
(4) |
即:
(5) |
梯度张量数据正演矩阵形式为
(6) |
其中,d是N维观测异常向量;m是M维剩余密度向量;A为N×M维灵敏度矩阵.
1.2 聚焦反演反演问题为最优化问题,即求解方程:
(7) |
由于反演问题的多解性和不稳定性,直接求解公式(7)的最优化问题无法获得稳定的唯一解,需要加入约束条件来获得合理的反演结果.吉洪诺夫正则化思想引入到反演计算中可有效地解决这一问题,正则化方程为
(8) |
其中φd(m)为数据拟合项;s(m)为模型约束项,也称作稳定函数;λ表示正则化参数.根据不同的先验约束信息,选择不同的稳定函数构建模型约束函数,常用的稳定函数是基于模型参数及梯度的不同阶范数构建的.
为了获得陡峭边界地质体的空间物性分布,Portniaguine和Zhdanov(1999)提出最小梯度支撑方程,公式为
(9) |
其中,e为极小值,e可避免|▽m|=0时奇点出现.e值选择时应保证|▽m|趋近于0时,
重力梯度全张量矩阵包括五个独立梯度分量,在进行航空重力梯度数据反演时选择五个梯度独立分量数据进行联合反演,能够更好地刻画地质体的空间特征.本文选择最小梯度支撑方程作为稳定函数的聚焦反演方法用于陡峭边界地质体的物性反演.聚焦反演正则化方程表达式为
(10) |
重力全张量矩阵五个独立分量联合反演的目标函数表达式为
(11) |
其中,dgxz、dgyz、dgzz、dgxx、dgxy是梯度张量数据的五个独立分量引起的重力异常;Axz、Ayz、Azz、Axx、Axy是对应梯度分量的灵敏度矩阵.为了获得公式(11)的反演解,采用非线性共轭梯度算法求解反演中的最优化问题(Newman and Alumbaugh, 2000).
1.3 深度加权函数由于趋肤效应的原因,随着反演深度的增加位场数据的核函数逐渐衰减,反演的物性结果集中在浅部.在实际数据处理中,一般引入深度加权函数抵消核函数的衰减效应,进而提高位场数据的深度分辨率.Li和Oldenburg(1998)引入了深度加权函数进行约束反演,其表达式为
(12) |
其中,z为网格单元中心点深度;z0为常数,取决于单元尺寸及观测面高度;β为常数,对于重力梯度数据β等于3.对目标函数模型项加入深度约束函数可有效的抵消核函数的衰减.Cella和Fedi(2012)指出深度衰减系数β与目标体的形态(构造指数)有关,在重力梯度数据反演中选择为确定值会在深度方向上产生较大的物性反演误差.
本文为了提高位场数据反演的深度分辨率,引入基于目标体中心埋深的深度加权函数(刘银萍等,2013),其表达式为
(13) |
其中θ是一个常数,其值越大对近地表压制越大,反之反演结果越发散,一般取θ=0.001;z为深度变量;d为地下空间纵向尺寸;r为常数,一般取1;zc表示异常体的中心埋深;α取值范围为0.8~0.9;β取值范围为24~26.通过深度估计或成像方法可以获取目标体的中心埋深位置,本文将获取的目标体深度信息作为先验约束信息参与到反演计算中,通过公式(13)对目标函数施加约束,可有效地提高反演结果的深度分辨率.
1.4 空间梯度加权函数Commer等(2011)提出空间梯度加权函数概念,并用于电磁数据的物性反演中.在进行非线性共轭梯度计算时,通过对目标函数中的数据拟合函数梯度施加空间梯度加权实现约束反演,表达式为
(14) |
其中f(x, y, z)为空间梯度加权函数,包括x、y、z三方向的空间约束:
(15) |
其中f(z)为深度加权函数(公式(13)),f(x)与f(y)分别为x和y方向的水平梯度加权函数,控制反演结果横向的聚焦范围.f(x)与f(y)具有相同的表达式,以f(x)为例,表达式为
(16) |
其中,x1、x2分别是x方向目标体边界位置;dx为x方向研究区域范围;α为常数,一般取0.001;r控制曲线的陡度.f(x)水平梯度加权函数对[x1, x2]区域给予较大权值,对区域外部分赋予较小权值.f(y)采用相同的定义方式并对y方向目标体[y1, y2]区域给予较大权值,通过f(x)与f(y)乘积得到的水平梯度加权函数对反演方程的约束可提高反演的水平分辨率.对于水平梯度加权函数中的目标体水平位置,本文选择几何参数反演方法快速的对边界位置参数进行求取.
空间梯度加权函数f(x, y, z)包含了基于深度参数信息的深度加权以及基于水平位置信息的水平梯度加权函数,可有效的结合几何参数反演获取的目标体空间位置信息,通过对目标函数中的数据拟合函数的梯度施加约束,进而控制了反演结果的物性参数空间分布.基于目标体参数先验信息的优化约束反演方法可在临近干扰异常、深浅叠加异常等复杂地质条件下获得更可靠、高分辨率的反演结果,有效地的提高了反演结果的空间分辨率.
此外,为了获得更合理的反演物性分布,本文采用稳定的不等式约束条件控制反演的物性参数在合理的分布范围内(Barbosa et al., 1999),在物性反演过程中需要加入地质体物性上下限约束条件,强制将反演密度值限制在一定的范围内.
2 理论模型试验 2.1 单一模型试验基于最小梯度支撑方程的聚焦反演方法可有效地反演出具有陡峭边界的地质体物性参数,建立如图 1所示的棱柱体模型说明传统聚焦反演方法的有效性以及两种深度加权函数的应用效果.建立棱柱体模型的空间位置如图 1所示,剩余密度为1000 kg·m-3,图 2g表示棱柱体引起的加入5%高斯噪声的重力异常图.图 2a—f为张量矩阵六个分量.
采用Li和Oldenburg(1998)提出的深度加权函数(公式12)对张量矩阵五个独立分量gxx、gxy、gxz、gyz、gzz进行联合反演,如图 3为反演的物性参数的水平及垂直方向剖面图.深度加权函数可以改善反演的趋肤效应,但反演异常体深度位置与实际埋深存在一定误差.
采用基于异常体中心埋深的深度加权函数可有效地提高深度方向的反演精度,为了获取公式(13)中地质体的中心埋深参数,本文采用归一化总梯度算法对梯度分量gzz进行成像(Zhou et al., 2017),成像结果的最大值对应于地质体的中心埋深位置.图 4为过棱柱体中心位置的归一化总梯度方法成像结果,最大值对应的中心埋深为145 m,将其作为先验信息约束到聚焦反演中,获得如图 5所示的反演结果.反演结果的深度分辨率得到明显的提升,重构异常体的水平位置及深度埋深与理论模型更加吻合.采用基于深度先验信息的深度加权函数进行反演计算,可有效地提高反演结果的深度分辨率.
建立如图 6所示的组合棱柱体模型,两个深浅棱柱体的中心埋深分别为110 m和150 m;剩余密度均为1000 kg·m-3;较浅棱柱体x和y方向水平边界位置范围分别为220~320 m和100 ~200 m,较深棱柱体x和y方向水平边界位置范围分别为220~320 m和360~460 m.图 6d为模型引起的重力异常图.图 7a、b为采用传统聚焦反演方法获得两个模型体物性反演结果的水平及垂直剖面,结果表明反演的较深棱柱体的物性分布在深度方向上具有较大的误差.
为了提高反演结果的水平及纵向分辨率,采用Zhou等(2017)提出的归一化边界识别方法EDTN对棱柱体模型异常三维成像.图 8a为不同深度位置EDTN成像结果切片图,成像结果的局部极大值与目标体的水平边界位置及中心埋深相对应.较浅棱柱体在110 m时出现极大值,较深棱柱体在150 m位置出现局部极大值.提取归一化边界识别方法局部极大值得到如图 8b的棱柱体空间位置,x方向棱柱体水平边界位置坐标范围分别为218~322 m和216~324 m,y方向棱柱体水平边界位置坐标范围分别为98~200 m和358~464 m,与实际模型体边界位置基本吻合.将快速成像方法计算的不同棱柱体水平边界位置及中心埋深应用到公式(16)空间梯度加权函数f(x, y, z)计算中,进而获得如图 7c—d的自约束物性反演的横向及垂向物性切片图.相比传统聚焦反演方法,加入约束条件后的反演结果可更好地反映深部棱柱体的物性分布特征,具有更高的水平及深度分辨率.
美国文顿岩丘位于路易斯安那州与德克萨斯州交界位置,该地与墨西哥湾盆地区相邻,发育盐丘构造.石灰岩、石膏和硬石膏组分的岩盖与巨大的盐核组成了盐丘构造.2008年,Bell Geospace公司对盐丘所在区域进行航空重力梯度全张量测量,测量范围自北纬30.07°至北纬30.23°,自西经93.53°至西经93.66°,飞行总测线长度1087.5 km,测量面积为192.6 km2,飞行高度大约为75 m.对实测数据进行向心力校正、自梯度校正、加速度计补偿、地形改正等数据改正后得到如图 9所示的全张量梯度数据,本文选择该区域的航空重力梯度数据进行优化约束物性反演.
针对航空重力梯度测量的五个独立梯度分量进行联合反演,选择基于目标体中心埋深的深度加权函数进行优化约束反演,其中的深度参数先验信息采用已有的多维斜导数深度估计方法获得,具体计算过程参考Zhou和Huang(2015).图 10为采用多维斜导数深度估计方法得到的两个高度的梯度比率R异常图,选择R=2等值线用于水平距离计算及深度反演.采用向上延拓100 m和150 m得到如图 10a、b的梯度比率异常图,虚线对应盐丘构造位置R=2的等值线,计算得到盐丘的中心埋深为345.5 m.
将基于多维数据斜导数深度估计方法获得的盐丘中心埋深作为先验信息约束到深度加权函数中,物性反演结果的水平及垂向剖面图如图 11所示.水平剖面图显示盐丘在东西方向延伸大约1.5 km,而南北方向延伸1 km左右.垂直剖面结果显示盐丘在深度方向有一定延深,绘制不同深度位置的物性三维切片图,如图 12所示.结果显示盐丘的高密度值在深度方向上大致从250 m延深至550 m,Oliveira和Barbosa(2013)给出岩丘的底面深度位置为460 m,Geng等(2014)基于协同克里格算法对全张量数据进行反演,获取的岩丘空间分布范围在200 ~650 m之间.本文方法反演结果与已有反演方法的解释结果基本吻合,验证了方法的实用性.
重磁数据物性反演存在多解性和不确定性问题,为了提高位场数据反演的水平及深度分辨率,本文充分挖掘原始观测数据中包含的目标体信息,采用几何参数快速反演方法获取目标体空间位置信息,并将其应用到物性反演的约束条件建立中.本文选择聚焦反演方法进行陡峭边界地质体的物性反演,并在反演之前首先采用快速反演方法提取出目标体空间位置参数,然后将其作为先验条件参与到深度加权函数及水平梯度加权函数的计算中,并通过空间梯度加权函数对数据拟合函数的梯度施加约束实现物性参数计算.单一目标体反演计算结果表明,加入深度先验信息的深度加权函数相比传统加权函数具有更高的深度分辨率.组合目标体反演结果表明,边界位置信息可有效地提高不同埋深目标体的水平分辨率.最后,将优化约束反演方法应用到美国文顿岩丘的实测航空重力全张量数据中,采用多维斜导数法反演的深度参数进行约束反演有效地提高了岩丘的垂向分辨率,进一步验证了本文建立方法的有效性和实用性.
Barbosa V C F, Silva J B C, Medeiros W E. 1999. Stable inversion of gravity anomalies of sedimentary basins with nonsmooth basement reliefs and arbitrary density contrast variations. Geophysics, 64(3): 754-764. DOI:10.1190/1.1444585 |
Carter-McAuslan A, Lelièvre P G, Farquharson C G. 2015. A study of fuzzy c-means coupling for joint inversion, using seismic tomography and gravity data test scenarios. Geophysics, 80(1): W1-W15. |
Cella F, Fedi M. 2012. Inversion of potential field data using the structural index as weighting function rate decay. Geophysical Prospecting, 60(2): 313-336. DOI:10.1111/gpr.2012.60.issue-2 |
Commer M. 2011. Three-dimensional gravity modelling and focusing inversion using rectangular meshes. Geophysical Prospecting, 59(5): 966-979. |
Commer M, Newman G A, Williams K H, et al. 2011. 3D induced-polarization data inversion for complex resistivity. Geophysics, 76(3): F157-F171. DOI:10.1190/1.3560156 |
Fedi M, Abbas M A. 2013. A fast interpretation of self-potential data using the depth from extreme points method. Geophysics, 78(2): E107-E116. DOI:10.1190/geo2012-0074.1 |
Fregoso E, Gallardo L A, García-Abdeslem J. 2015. Structural joint inversion coupled with Euler deconvolution of isolated gravity and magnetic anomalies. Geophysics, 80(2): G67-G79. DOI:10.1190/geo2014-0194.1 |
Geng M X, Huang D N, Yang Q J, et al. 2014. 3D inversion of airborne gravity-gradiometry data using cokriging. Geophysics, 79(4): G37-G47. DOI:10.1190/geo2013-0393.1 |
Giraud J, Pakyuz-Charrier E, Jessell M, et al. 2017. Uncertainty reduction through geologically conditioned petrophysical constraints in joint inversion. Geophysics, 82(6): ID19-ID34. DOI:10.1190/geo2016-0615.1 |
Guillen A, Menichetti V. 1984. Gravity and magnetic inversion with minimization of a specific functional. Geophysics, 49(8): 1354-1360. DOI:10.1190/1.1441761 |
Guo M, Deng J Z, Chen X, et al. 2018. Bayesian joint inversion of magnetotelluric and gravity based on petrophysical Constraints. Progress in Geophysics (in Chinese), 33(5): 1897-1902. DOI:10.6038/pg2018BB0461 |
Ialongo S, Fedi M, Florio G. 2014. Invariant models in the inversion of gravity and magnetic fields and their derivatives. Journal of Applied Geophysics, 110: 51-62. DOI:10.1016/j.jappgeo.2014.07.023 |
Last B J, Kubik K. 1983. Compact gravity inversion. Geophysics, 48(6): 713-721. DOI:10.1190/1.1441501 |
Lelièvre P G, Oldenburg D W. 2009. A comprehensive study of including structural orientation information in geophysical inversions. Geophysical Journal International, 178(2): 623-637. DOI:10.1111/gji.2009.178.issue-2 |
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 |
Liu Y P, Wang Z W, Du X J, et al. 2013. 3D constrained inversion of gravity data based on the Extrapolation Tikhonov regularization algorithm. Chinese Journal of Geophysics (in Chinese), 56(5): 1650-1659. |
Lo Re D, Florio G, Ferranti L, et al. 2016. Self-constrained inversion of microgravity data along a segment of the Irpinia fault. Journal of Applied Geophysics, 124: 148-154. DOI:10.1016/j.jappgeo.2015.12.002 |
Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International, 140(2): 410-424. DOI:10.1046/j.1365-246x.2000.00007.x |
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, Ialongo S, Florio G, et al. 2013. Self-constrained inversion of potential fields. Geophysical Journal International, 195(2): 854-869. DOI:10.1093/gji/ggt313 |
Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596 |
Silva J B C, Medeiros W E, Barbosa V C F. 2001. Potential-field inversion:choosing the appropriate technique to solve a geologic problem. Geophysics, 66(2): 511-520. DOI:10.1190/1.1444941 |
Sun S D, Chen C. 2016. A self-constrained inversion of magnetic data based on correlation method. Journal of Applied Geophysics, 135: 8-16. DOI:10.1016/j.jappgeo.2016.09.022 |
Yin C C, Sun S Y, Gao X H, et al. 2018. 3D joint inversion of magnetotelluric and gravity data based on local correlation constraints. Chinese Journal of Geophysics, 61(1): 358-367. DOI:10.6038/cjg2018K0765 |
Zhang Q, Yu P, Zhang D L, et al. 2018. Preliminary inversion of the potential field model based on DEXP. Progress in Geophysics (in Chinese), 33(5): 2076-2082. DOI:10.6038/pg2018BB0558 |
Zhang Y, Yan J G, Li F, et al. 2015. A new bound constraints method for 3-D potential field data inversion using Lagrangian multipliers. Geophysical Journal International, 201(1): 267-275. DOI:10.1093/gji/ggv016 |
Zhou S, Huang D N. 2015. Determining the depth of certain gravity sources without a priori specification of their structural index. Journal of Applied Geophysics, 122: 172-180. DOI:10.1016/j.jappgeo.2015.09.021 |
Zhou S, Huang D N, Jiao J. 2016. A filter to detect edge of potential field data based on three-dimensional structural tensors. Chinese Journal of Geophysics (in Chinese), 59(10): 3847-3858. DOI:10.6038/cjg20161028 |
Zhou S, Jiang J B, Lu P Y. 2017. Source parameters estimation by the normalized downward continuation of gravity gradient data. Journal of Applied Geophysics, 147: 81-90. DOI:10.1016/j.jappgeo.2017.10.011 |
郭曼, 邓居智, 陈晓, 等. 2018. 基于岩石物性约束的大地电磁与重力贝叶斯联合反演. 地球物理学进展, 33(5): 1897-1902. DOI:10.6038/pg2018BB0461 |
刘银萍, 王祝文, 杜晓娟, 等. 2013. 基于Extrapolation Tikhonov正则化算法的重力数据三维约束反演. 地球物理学报, 56(5): 1650-1659. |
殷长春, 孙思源, 高秀鹤, 等. 2018. 基于局部相关性约束的三维大地电磁数据和重力数据的联合反演. 地球物理学报, 61(1): 358-367. DOI:10.6038/cjg2018K0765 |
张琦, 于平, 张代磊, 等. 2018. 基于DEXP的位场模型初步反演. 地球物理学进展, 33(5): 2076-2082. DOI:10.6038/pg2018BB0558 |
周帅, 黄大年, 焦健. 2016. 基于三维构造张量的位场边界识别滤波器. 地球物理学报, 59(10): 3847-3858. DOI:10.6038/cjg20161028 |