﻿ 基于密度模型稀疏表征的重力反演方法
 地球物理学报  2021, Vol. 64 Issue (3): 1061-1073 PDF

Gravity inversion based on sparse representation of density model
YU HuiZhen, WANG JinDuo, WANG QianJun
Research Institute of Exploration and Development, Shengli Oilfield, SINOPEC, Shandong Dongying 257000, China
Abstract: Gravity inversion is an essential tool for imaging subsurface density distribution, and the selection of a proper density model constraining method is key for improving gravity inversion resolution and reliability. Conventional constraining methods match inversion target via adjusting smoothing or sparse constrain weight starting from density model in the mesh space. However, it is difficult to guarantee the reasonability and validity of model constrains when there exists a set of multi-style geologic bodies, an inaccurate anomaly separation or an unreasonable mesh-generation solution. To this end, a novel gravity inversion method has been proposed which is based on sparse representation of density model. Firstly, the representation of the density model to be recovered is assumed to be the linear combination of model feature matrix and sparse decomposition coefficient, then objective function of gravity inversion is re-derived and sparse solving process of decomposition coefficient also is provided. Comparing with traditional gravity inversion methods, the feature model for constructing model feature matrix contains prior geometrical information of multi-style geologic bodies, and the sparsity of decomposition coefficient ensures that the inversion target comes from combination of the most typical geological modes. Finally, the test of synthetic models and the application of field data demonstrate the validity of our proposed gravity inversion method which is based on sparse representation of density model.
Keywords: Sparse representation    Gravity inversion    Model feature matrix    Sparse algorithm
0 引言

1 方法原理 1.1 重力反演方法

 (1)

 (2)

 (3)

k次迭代的密度模型Mk满足的目标函数形式如下：

 (4)

 (5)

1.2 模型特征矩阵构建流程

 图 1 一维密度模型分解示意图 Fig. 1 Diagram of 1D density model decomposition

 (6)

 (7)

(1) 根据已知资料分析，假设地质体发育模式包含以下两类矩形块体特征模型，具体参数为：特征模型Ⅰ(图 2a)：宽W1m、高H1m；特征模型Ⅱ(图 2b)：宽W2m、高H2m；

 图 2 构建模型特征矩阵的特征模型 (a) 特征模型Ⅰ；(b) 特征模型Ⅱ. Fig. 2 Feature models for building model feature matrix (a) Feature model Ⅰ; (b) Feature model Ⅱ.

(2) 首先计算特征模型Ⅰ对应的索引向量Index_I.将特征模型I的中心点坐标置于二维剖分网格左上角第一个网格(图 3a)，以此为起始点进行第1次搜索，计算特征模型I所占反演剖分的索引，索引Index_I_1位置的参数记为整数1，其余为0(图 3b)；向右平移一个网格的距离至左上角第2个网格，计算特征模型I所占反演网格的索引，并将其记录为Index_I_2，索引Index_I_2位置的参数记为整数1，其余为0；以此类推，计算完第一层之后，从第二层的左侧第一个网格点为起始点计算特征模型I的所占网格位置，并进行索引标记；之后，逐层计算，设平移至第22个网格(X=4、Z=3)(图 3c)进行第22次搜索，索引Index_I_22标记如图(图 3d)所示；按此步骤直到遍历所有网格；

 图 3 特征模型Ⅰ对应的模型特征矩阵构建过程 (a) 第1次搜索位置; (b) 第1次的索引; (c) 第22次搜索位置; (d) 第22次的索引. Fig. 3 Building process of model feature matrix corresponding to feature model Ⅰ (a) The first search location; (b) The first index; (c) The 22nd search location; (d) The 22nd index.

(3) 将步骤(2)得到的索引Index_I_k(k=1, …, n)，n为剖分网格个数，此次为63，展开为列向量并进行组合得到特征模型I对应的模型特征矩阵DB1

(4) 特征模型Ⅱ的模型特征矩阵构建.重新执行第(2)、(3)步，计算特征模型Ⅱ下的索引向量，并将其进行组合得到特征模型Ⅱ对应的特征矩阵DB2

(5) 背景模型特征矩阵.类似特征模型Ⅰ、Ⅱ的方式，选择具有趋势特征的模型作为特征模型，将其生成背景模型特征矩阵Dbg

(6) 完成上述步骤，便可生成二维反演所需的模型特征矩阵D=[DB1, DB2, Dbg].

1.3 目标函数及求解算法

 (8)

 (9)

 (10)

WΓk-1为第k-1次迭代求解结果Γk-1对应的对角加权矩阵；ϕ′(Γk-1)为模型约束函数ϕ(Γ)关于Γk-1的一阶偏导，其计算公式为：

 (11)

 (12)

 (13)

(1) 根据实际需求剖分原始反演网格，计算重力正演核函数矩阵A

(2) 根据工区实际情况、重力异常形态等信息设置特征模型(可包括块体、倾斜状岩脉、球体等各类典型地质体)，根据1.2节介绍的过程构建模型特征矩阵D

(3) 设置正则化参数λ1λ2及深度加权参数N(一般N=4)，设置最大迭代次数IterMax

(4) 根据模型特征矩阵列向量个数设置稀疏分解系数初始模型向量Γk-1=ΓInit(k=1)，注意ΓInit不要含有0值；

(5) 根据公式(11)，计算稀疏约束迭代加权矩阵WΓk-1-1

(6) 利用公式(10)和(13)计算第k次的分解系数Γk

(7) 利用公式(6)得到第k次迭代对应的密度模型Mk=Wdepthk

(8) 获得公式(9)的计算结果，若满足目标函数收敛条件则执行步骤(9)；若不满足，重复步骤(5)—(7)，修改反演参数；

(9) 评价(8)获得的反演结果，若符合地质认识，则执行步骤(10)；若不满足，重复(2)—(8)，修改特征模型，重新构建模型特征矩阵并求取反演结果；

(10) 最终的密度反演结果输出.

2 理论模型实验分析

2.1 二维模型

(1) 模型一：无背景模型干扰的两个地质体

 图 4 反演结果对比 (a) 重力异常(黑线)及拟合数据(红线); (b) 密度模型; (c) 常规聚焦反演结果; (d) 本文反演结果. Fig. 4 Comparison of inversion results (a) Gravity anomaly (black line) and fitting data (red line); (b) Density model; (c) The results of conventional focus inversion; (d) The results of this paper.

(2) 模型二：存在背景模型干扰的两个地质体

 图 5 反演结果对比 (a) 重力异常(黑线)及拟合数据(红线); (b) 密度模型; (c) 常规聚焦反演结果; (d) 本文反演结果. Fig. 5 Comparison of inversion results (a) Gravity anomaly (black line) and fitting data (red line); (b) Density model; (c) The results of conventional focus inversion; (d) The results of this paper.

(3) 模型三：存在背景模型干扰的倾斜状地质体

 图 6 构建模型特征矩阵的二维特征模型 Fig. 6 2D feature model for building model feature matrix

 图 7 反演结果对比 (a) 重力异常(黑线)及拟合数据(红线); (b) 密度模型; (c) 常规聚焦反演结果; (d) 本文反演结果. Fig. 7 Comparison of inversion results (a) Gravity anomaly (black line) and fitting data (red line); (b) Density model; (c) The results of conventional focus inversion; (d) The results of this paper.
2.2 三维模型

 图 8 三维密度模型及重力异常 (a) 三维密度模型；(b) 重力正演异常(含3%高斯噪声). Fig. 8 3D density model and gravity anomaly (a) 3D density model; (b) Gravity forward anomaly (including 3% Gaussian noise).

 图 9 构建模型特征矩阵的三维特征模型 (a) 三维特征模型a；(b) 三维特征模型b. Fig. 9 3D feature model for building model feature matrix (a) 3D feature model a; (b) 3D feature model b.

 图 10 反演结果对比 (a) 理论模型A-A′垂直剖面；(b) 理论模型B-B′垂直剖面；(c) 聚焦反演A-A′垂直剖面；(d) 聚焦反演B-B′垂直剖面；(e) 本文方法A-A′垂直剖面；(f) 本文方法B-B′垂直剖面. Fig. 10 Comparison of inversion results (a) A-A′ profile of theoretical model; (b) B-B′ profile of theoretical model; (c) A-A′ profile of focus inversion; (d) B-B′ profile of focus inversion; (e) A-A′ profile of proposed method; (f) B-B′ profile of proposed method.

3 实际资料应用测试

 图 11 实际资料反演结果 (a) 剩余重力异常; (b) 密度三维反演结果; (c) A-A′垂直剖面反演结果；(d) B-B′垂直剖面反演结果. (c)—(d)中黑色多边形代表硫化物矿床的位置. Fig. 11 Inversion results of actual data (a) Residual gravity anomaly; (b) 3D inversion result of density; (c) The result of A-A′ vertical profile; (d) The result of B-B′ vertical profile.The sulphide location is indicated by the black polygon in (c)—(d).

4 结论

 (A1)

 (A2)

 (A3)

 (A4)

 (A5)

 (A6)

ΓΓv都为大小为q×1的向量，则根据MM思想：

 (A7)

 (A8)

 (A9)

 (A10)

 (A11)

References
 Bertete-Aguirre H, Cherkaev E, Oristaglio M. 2002. Non-smooth gravity problem with total variation penalization functional. Geophysical Journal International, 149(2): 499-507. DOI:10.1046/j.1365-246X.2002.01664.x Boyd S, Parikh N, Chu E, et al. 2010. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends® in Machine Learning, 3(1): 1-122. DOI:10.1561/2200000016 Chen S S, Donoho D L, Saunders M A. 2001. Atomic decomposition by basis pursuit. SIAM Review, 43(1): 129-159. DOI:10.1137/S003614450037906X Farquharson C G, Oldenburg D W. 1998. Non-linear inversion using general measures of data misfit and model structure. Geophysical Journal International, 134(1): 213-227. DOI:10.1046/j.1365-246x.1998.00555.x Farquharson C G. 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversions. Geophysics, 73(1): K1-K9. DOI:10.1190/1.2816650 Figueiredo M A T, Bioucas-Dias J M, Nowak R D. 2007. Majorization-minimization algorithms for wavelet-based image restoration. IEEE Transactions on Image Processing, 16(12): 2980-2991. DOI:10.1109/TIP.2007.909318 Gao X H, Huang D N. 2017. Research on 3D focusing inversion of gravity gradient tensor data based on a conjugate gradient algorithm. Chinese Journal of Geophysics (in Chinese), 60(4): 1571-1583. DOI:10.6038/cjg20170429 Guan Z N, Hou J S, Huang L P, et al. 1998. Inversion of gravity and magnetic anomalies using pseudo-BP neural network method and its application. Chinese Journal of Geophysics (Acta Geophysica Sinica) (in Chinese), 41(2): 242-251. 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/j.1365-246X.2009.04188.x 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 Z L, Yao C L, Zheng Y M. 2019. 3D inversion of gravity data using Lp-norm sparse optimization. Chinese Journal of Geophysics (in Chinese), 62(10): 3699-3709. DOI:10.6038/cjg2019M0430 Liu Z, Yu H Z, Chen T. 2011. 3D inversion method of density based on double constraint. Journal of China University of Petroleum (Edition of Natural Science) (in Chinese), 35(6): 43-50. Meng X H, Liu G F, Chen Z X, et al. 2012. 3-D gravity and magnetic inversion for physical properties based on residual anomaly correlation. Chinese Journal of Geophysics (in Chinese), 55(1): 304-309. DOI:10.6038/j.issn.0001-5733.2012.01.030 Nguyen H D. 2017. An introduction to Majorization-Minimization algorithms for machine learning and statistical estimation. WIREs Data Mining and Knowledge Discovery, 7(3): e1198. Oldenburg D W, Li Y G, Ellis R G. 1997. Inversion of geophysical data over a copper gold porphyry deposit: A case history for Mt. Milligan. Geophysics, 62(5): 1419-1431. DOI:10.1190/1.1444246 Paterson N R, Reeves C V. 1985. Applications of gravity and magnetic surveys: The state-of-the-art in 1985. Geophysics, 50(12): 2558-2594. DOI:10.1190/1.1441884 Peng G M, Liu Z, Yu H Z, et al. 2018. 3D gravity inversion based on Cauchy distribution constraint and fast proximal objective function optimization. Chinese Journal of Geophysics (in Chinese), 61(12): 4934-4941. DOI:10.6038/cjg2018L0776 Phillips N, Oldenburg D, Chen J P, et al. 2001. Cost effectiveness of geophysical inversions in mineral exploration: Applications at San Nicolas. The Leading Edge, 20(12): 1351-1360. DOI:10.1190/1.1487264 Portniaguine O, Zhdanov M S. 1999. Focusing geophysical inversion images. Geophysics, 64(3): 874-887. DOI:10.1190/1.1444596 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 Sahoo S K, Makur A. 2015. Signal recovery from random measurements via extended orthogonal matching pursuit. IEEE Transactions on Signal Processing, 63(10): 2572-2581. DOI:10.1109/TSP.2015.2413384 Selesnick I W, Bayram İ. 2014. Sparse signal estimation by maximally sparse convex optimization. IEEE Transactions on Signal Processing, 62(5): 1078-1092. DOI:10.1109/TSP.2014.2298839 Sun J J, Li Y G. 2014. Adaptive Lp inversion for simultaneous recovery of both blocky and smooth features in a geophysical model. Geophysical Journal International, 197(2): 882-899. DOI:10.1093/gji/ggu067 Utsugi M. 2019. 3-D inversion of magnetic data based on the L1-L2 norm regularization. Earth, Planets and Space, 71(1): 73. DOI:10.1186/s40623-019-1052-4 Vatankhah S, Renaut R A, Ardestani V E. 2017. 3-D Projected L1 inversion of gravity data using truncated unbiased predictive risk estimator for regularization parameter estimation. Geophysical Journal International, 210(3): 1872-1887. DOI:10.1093/gji/ggx274 Yao C L, Hao T Y, Guan Z N, et al. 2003. High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms. Chinese Journal of Geophysics (in Chinese), 46(2): 252-258. DOI:10.3321/j.issn:0001-5733.2003.02.020 Zhdanov M S. 2002. Geophysical Inverse Theory and Regularization Problems. Methods in Geochemistry and Geophysics. Amsterdam: Elsevier Science, 36: 3-609. 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 Zhu Z Q, Cao S J, Lu G Y. 2014. 3D gravity inversion with bound constraint based on hyper-parameter regularization. The Chinese Journal of Nonferrous Metals (in Chinese), 24(10): 2601-2608. 高秀鹤, 黄大年. 2017. 基于共轭梯度算法的重力梯度数据三维聚焦反演研究. 地球物理学报, 60(4): 1571-1583. DOI:10.6038/cjg20170429 管志宁, 侯俊胜, 黄临平, 等. 1998. 重磁异常反演的拟BP神经网络方法及其应用. 地球物理学报, 41(2): 242-251. DOI:10.3321/j.issn:0001-5733.1998.02.013 李泽林, 姚长利, 郑元满. 2019. 基于Lp范数稀疏优化算法的重力三维反演. 地球物理学报, 62(10): 3699-3709. DOI:10.6038/cjg2019M0430 刘展, 于会臻, 陈挺. 2011. 双重约束下的密度三维反演. 中国石油大学学报(自然科学版), 35(6): 43-50. 孟小红, 刘国峰, 陈召曦, 等. 2012. 基于剩余异常相关成像的重磁物性反演方法. 地球物理学报, 55(1): 304-309. DOI:10.6038/j.issn.0001-5733.2012.01.030 彭国民, 刘展, 于会臻, 等. 2018. 基于柯西分布约束和快速近端目标函数优化的三维重力反演方法. 地球物理学报, 61(12): 4934-4941. DOI:10.6038/cjg2018L0776 秦朋波, 黄大年. 2016. 重力和重力梯度数据联合聚焦反演方法. 地球物理学报, 59(6): 2203-2224. DOI:10.6038/cjg20160624 姚长利, 郝天珧, 管志宁, 等. 2003. 重磁遗传算法三维反演中高速计算及有效存储方法技术. 地球物理学报, 46(2): 252-258. DOI:10.3321/j.issn:0001-5733.2003.02.020 朱自强, 曹书锦, 鲁光银. 2014. 基于混合正则化的重力场约束反演. 中国有色金属学报, 24(10): 2601-2608.