﻿ 基于L<sub><i>p</i></sub>范数稀疏优化算法的重力三维反演
 地球物理学报  2019, Vol. 62 Issue (10): 3699-3709 PDF

3D inversion of gravity data using Lp-norm sparse optimization
LI ZeLin, YAO ChangLi, ZHENG YuanMan
Key Laboratory of Geo-detection(China University of Geosciences, Beijing), Ministry of Education, Beijing 100083, China
Abstract: Three-dimensional inversion for density distribution has become a common tool for quantitative interpretation of gravity data in recent years. However, as gravity data lacks depth resolution, such inversion has severe non-uniqueness. To reduce this problem, the most common approach is to introduce extra prior information. To further perfect this tool, we propose a 3D sparse gravity inversion method which minimizes an Lp-norm (0 ≤ p ≤ 1) of the density model subject to bound constraints. Compared with the traditional L2-norm inversion approach, our method permits to use known physical property information more effectively and produce solutions characterized by sharp boundaries and high depth resolution. To better understand our method, we analyze the equivalence between our method and binary or ternary inversion. Moreover, we also point out several issues of this method that need to be noticed. The validity of our method is tested by both synthetic-and real-data examples.
Keywords: Three-dimensional gravity inversion    Sparse    Lp-norm    Physical property information
0 引言

Li等(2018)经过分析后认为，在磁数据反演中，当模型范数中仅包含模型长度项(即不包含模型导数相关项)时，第二类算法和第一类算法中的二值反演是等价的.Li等(2018)通过磁数据模型试验和实际数据试验验证了两类算法的等价性.在本文中，我们将Li等(2018)提出的Lp范数稀疏优化算法扩展到重力数据反演中，分析了稀疏反演算法与二值、三值反演的等价性，指出了其在实际应用时需要注意的问题，最后通过模型试验和矿区数据反演验证了稀疏反演算法的有效性.

1 重力Lp范数稀疏反演方法 1.1 正演问题

 (1)

1.2 反演问题

 (2)

 (3)

 (4)

 (5)

 (6)

mj=0时，(6)式无定义，为了防止(6)式奇异，将(6)式修改为

 (7)

 (8)

 (9)

(8) 式属于不等式约束线性反演问题，对于该问题的求解，当前已有多种算法成功应用到了重磁领域，如非线性投影法(Li and Oldenburg, 1996Lelièvre and Oldenburg, 2006)，内点法(Li and Oldenburg, 2003)，以及梯度投影法(Lelièvre et al., 2009).本文选择内点法求解该问题，详细的算法见附录A.

(1) 利用附录A的算法求解(8)式得到m(n).

(2) nn+1，更新μ(n)，按照(7)式由m(n-1)生成W(n).

1.3 稀疏反演方法与二值和三值反演的等价性分析

Li等(2018)经过分析后认为，在磁测数据反演中，当模型范数中仅包含模型长度项(即不包含模型导数相关项)时，稀疏反演和二值反演是等价的.我们认为，该等价性也可以扩展到重力反演中.

 图 1 不同物性上界约束条件下球体模型的反演结果中心切片图 (a—b)、(c—d)以及(e—f)分别为物性上界设置为2、1以及0.5 g·cm-3时的反演结果.物性下界均设置为0 g·cm-3.黑框代表了真实模型的位置. Fig. 1 Inversion results of the anomaly produced by the sphere model by using different upper bounds (a—b) The inversion result with an upper bound of 2 g·cm-3. (c—d) The inversion result with an upper bound of 1 g·cm-3. (e—f) The inversion result with an upper bound of 0.5 g·cm-3. The lower bounds of all inversions are set to 0 g·cm-3. The position of the true model is indicated by the black lines.

2 模型试验

2.1 倾斜薄板模型

 图 2 薄板模型所产生的重力数据 Fig. 2 The gravity anomaly produced by the thin dipping slab model

 图 3 倾斜薄板模型不同方法的反演结果中心切片图 (a—b)无约束L2范数反演结果；(c—d)物性上下界约束L2范数反演结果；(e—f)物性上下界约束Lp范数反演结果，其中p设置为0.在物性约束反演中，物性上、下界分别设置为1和0 g·cm-3.黑框代表了真实模型的位置. Fig. 3 Inversion results of the gravity anomaly produced by the thin dipping slab model by using different methods (a—b) The result of unconstrained L2-norm inversion. (c—d) The result of bound-constrained L2-norm inversion. (e—f) The result of bound-constrained Lp-norm inversion with p=0. The upper and lower bounds of bound-constrained inversions are set to 1 and 0 g·cm-3, respectively. The position of the true model is indicated by the black lines.
2.2 组合块体模型

 图 4 组合长方体模型所产生的重力数据 Fig. 4 The gravity anomaly produced by the model composed of two rectangular prisms

 图 5 组合长方体模型重力数据不同方法的反演结果切片图 (a—b)物性上下界约束L2范数反演结果；(c—d)物性上下界约束Lp范数反演结果，其中p设置为0.物性上、下界分别设置为1和-1 g·cm-3.黑框代表了模型的真实位置. Fig. 5 Inversion results of the gravity anomaly produced by the model composed of two rectangular prisms by using different methods (a—b) The result of bound-constrained L2-norm inversion. (c—d) The result of bound-constrained Lp-norm inversion with p=0. The upper and lower bounds of both inversions are set to 1 and-1 g·cm3, respectively. The position of the true model is indicated by the black lines.
2.3 组合倾斜板状体模型

 图 6 组合倾斜板状体模型所产生的重力数据 Fig. 6 The gravity anomaly produced by the model composed of two dipping slabs

 图 7 组合倾斜板状体模型不同方法的反演结果中心切片图 (a—b)无约束L2范数反演结果；(c—d)物性上下界约束L2范数反演结果；(e—f)物性上下界约束Lp范数反演结果，其中p设置为0.在物性约束反演中，物性上、下界分别设置为1和0 g·cm-3.黑框代表了真实模型的位置. Fig. 7 Inversion results of the gravity anomaly produced by the model composed of two dipping slabs by using different methods (a—b) The result of unconstrained L2-norm inversion; (c—d) The result of bound-constrained L2-norm inversion; (e—f) The result of bound-constrained Lp-norm inversion with p=0. The upper and lower bounds of bound-constrained inversions are set to 1 and 0 g·cm-3, respectively. The position of the true model is indicated by the black lines.

2.4 不准确的物性信息对反演结果的影响

 图 8 倾斜薄板模型不同物性上界的反演结果中心切片图 (a—b)物性上界设置为0.5 g·cm-3时的反演结果；(c—d)物性上界设置为2 g·cm-3时的反演结果.黑框代表了真实模型的位置. Fig. 8 Inversion results of the gravity anomaly produced by the thin dipping slab model by using different upper bounds (a—b) The inversion result with an upper bound of 0.5 g·cm-3; (c—d) The inversion result with an upper bound of 2 g·cm-3. The position of the true model is indicated by the black lines.
2.5 不同类型的干扰场源对反演结果的影响

 图 9 不同组合块体模型重力数据的反演结果切片图 (a—b)两个块体密度均为正值时的反演结果; (c—d)两个块体密度一正一负时的反演结果.黑框代表了真实模型的位置. Fig. 9 Inversion results of gravity anomalies produced by different models composed of two rectangular prisms (a—b) The inversion result of the gravity anomaly produced by the model composed of two positive prisms; (c—d) The inversion result of the gravity anomaly produced by the model composed of a positive prism and a negative prism. The position of the true model is indicated by the black lines.
3 实际资料试验

 图 10 圣尼古拉斯硫化物铜锌矿区剩余重力数据(两条互相垂直的黑线代表了已知地质剖面的位置) Fig. 10 The gravity anomaly of San Nicolas massive sulfide copper-zinc deposit (The two orthogonal lines indicate the positions of known geological sections)

 图 11 圣尼古拉斯硫化物铜锌矿区重力数据不同方法的反演结果切片图 (a—b)物性上下界约束L2范数反演结果；(c—d)物性上下界约束全变分反演结果；(e—f)物性上下界约束Lp范数反演结果，其中p设置为0.物性上、下界分别设置为1.4和-0.2 g·cm-3.黑框代表了已知矿体的位置. Fig. 11 Inversion results of the gravity anomaly of San Nicolas massive sulfide copper-zinc deposit by using different methods (a—b) The result of bound-constrained L2-norm inversion; (c—d) The result of bound-constrained total varation inversion; (e—f) The result of bound-constrained Lp-norm inversion with p=0. The upper and lower bounds of both inversions are set to 1.4 and-0.2 g·cm-3, respectively. The position of the ore body is indicated by the black lines.
4 结论

 (A1)

 (A2)

 (A3)

 (A4)

 (A5)

(1) 由mk-1生成XkYk.

(2) 采用预处理共轭梯度法求解(A5)式得到Δm.

(3) mkmk-1+γβΔm，其中γ=0.925，

(4) λk←[1-min(β, γ)]λk-1.

(5) kk+1.

References
 Aster R C, Borchers B, Thurber C H. 2005. Parameter Estimation and Inverse Problems. Boston, MA: Academic Press. Blakely R J. 1996. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press. Boulanger O, Chouteau M. 2001. Constraints in 3D gravity inversion. Geophysical Prospecting, 49(2): 265-280. DOI:10.1046/j.1365-2478.2001.00254.x Camacho A G, Montesinos F G, Vieira R. 2000. Gravity inversion by means of growing bodies. Geophysics, 65(1): 95-101. DOI:10.1190/1.1444729 Daubechies I, DeVore R, Fornasier M, et al. 2010. Iteratively reweighted least squares minimization for sparse recovery. Communications on Pure and Applied Mathematics, 63(1): 1-38. Farquharson C G. 2008. Constructing piecewise-constant models in multidimensional minimum-structure inversions. Geophysics, 73(1): K1-K9. DOI:10.1190/1.2816650 Guan Z N, Hou J S, Huang L P, et al. 1998. Inversion of gravity and magnetic anomalies using pseduo-BP neural network method and its application. Acta Geophysica Sinica (in Chinese), 41(2): 242-251. Holland P W, Welsch R E. 1977. Robust regression using iteratively reweighted least-squares. Communications in Statistics-Theory and Methods, 6(9): 813-827. DOI:10.1080/03610927708827533 Krahenbuhl R A, Li Y G. 2006. Inversion of gravity data using a binary formulation. Geophysical Journal International, 167(2): 543-556. DOI:10.1111/j.1365-246X.2006.03179.x Krahenbuhl R A, Li Y G. 2009. Hybrid optimization for lithologic inversion and time-lapse monitoring using a binary formulation. Geophysics, 74(6): I55-I65. DOI:10.1190/1.3242271 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. 2006. Magnetic forward modelling and inversion for high susceptibility. Geophysical Journal International, 166(1): 76-90. DOI:10.1111/j.1365-246X.2006.02964.x 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 Lelièvre P G, Oldenburg D W, Williams N C. 2009. Integrating geological and geophysical data through advanced constrained inversions. Exploration Geophysics, 40(4): 334-341. DOI:10.1071/EG09012 Li W B, Lu W T, Qian J L, et al. 2017. A multiple level-set method for 3D inversion of magnetic data. Geophysics, 82(5): J61-J81. DOI:10.1190/geo2016-0530.1 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. 2003. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method. Geophysical Journal International, 152(2): 251-265. DOI:10.1046/j.1365-246X.2003.01766.x Li Z L, Yao C L, Zheng Y M, et al. 2018. 3D magnetic sparse inversion using an interior-point method. Geophysics, 83(3): 1-56. Liu S, Hu X Y, Liu T Y, et al. 2013. Magnetization vector imaging for borehole magnetic data based on magnitude magnetic anomaly. Geophysics, 78(6): D429-D444. DOI:10.1190/geo2012-0454.1 Liu S, Hu X Y, Xi Y F, et al. 2015. 2D sequential inversion of total magnitude and total magnetic anomaly data affected by remanent magnetization. Geophysics, 80(3): K1-K12. DOI:10.1190/geo2014-0019.1 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. DOI:10.6038/cjg20130522 Lu W T, Qian J L. 2015. A local level-set method for 3D inversion of gravity-gradient data. Geophysics, 80(1): G35-G51. DOI:10.1190/geo2014-0188.1 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 Menke W. 2012. Geophysical Data Analysis: Discrete Inverse Theory: MATLAB Edition. New York: Academic Press. Nocedal J, Wright S. 2006. Numerical Optimization. New York: Springer. Phillips N D. 2001. Geophysical inversion in an integrated exploration program: Examples from the San Nicolas deposit [Master's thesis]. Vancouver: University of British Columbia. 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. 3D magnetic inversion with data compression and image focusing. Geophysics, 67(5): 1532-1541. DOI:10.1190/1.1512749 Rao B D, Kreutz-Delgado K. 1999. An affine scaling methodology for best basis selection. IEEE Transactions on signal processing, 47(1): 187-200. DOI:10.1109/78.738251 Uieda L, Barbosa V C F. 2012. Robust 3D gravity gradient inversion by planting anomalous densities. Geophysics, 77(4): G55-G66. DOI:10.1190/geo2011-0388.1 Van Zon T, Roy-Chowdhury K. 2006. Structural inversion of gravity data using linear programming. Geophysics, 71(3): J41-J50. DOI:10.1190/1.2197491 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. Yao C L, Zheng Y M, Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces. Chinese Journal of Geophysics (in Chinese), 50(5): 1576-1583. 管志宁, 侯俊胜, 黄临平, 等. 1998. 重磁异常反演的拟BP神经网络方法及其应用. 地球物理学报, 41(2): 242-251. DOI:10.3321/j.issn:0001-5733.1998.02.013 刘银萍, 王祝文, 杜晓娟, 等. 2013. 基于Extrapolation Tikhonov正则化算法的重力数据三维约束反演. 地球物理学报, 56(5): 1650-1659. DOI:10.6038/cjg20130522 孟小红, 刘国峰, 陈召曦, 等. 2012. 基于剩余异常相关成像的重磁物性反演方法. 地球物理学报, 55(1): 304-309. DOI:10.6038/j.issn.0001-5733.2012.01.030 姚长利, 郝天珧, 管志宁, 等. 2003. 重磁遗传算法三维反演中高速计算及有效存储方法技术. 地球物理学报, 46(2): 252-258. DOI:10.3321/j.issn:0001-5733.2003.02.020 姚长利, 郑元满, 张聿文. 2007. 重磁异常三维物性反演随机子域法方法技术. 地球物理学报, 50(5): 1576-1583. DOI:10.3321/j.issn:0001-5733.2007.05.035