﻿ 基于阻尼最小二乘法的位场欧拉反演方法
 地球物理学报  2019, Vol. 62 Issue (10): 3710-3722 PDF

Euler deconvolution of potential field based on damped least square method
LIU Qiang, YAO ChangLi, ZHENG YuanMan
Key Laboratory of Geo-detection(China University of Geosciences, Beijing), Ministry of Education, Beijing 100083, China
Abstract: Euler deconvolution is a popular method to estimate the source location. But the problem that how to obtain stable and reliable solutions has not been solved yet. The interference from adjacent sources which may be a main interference factor could cause the normal equation matrix (GTG) singular (GTG does not exist) and result in generating spurious solutions. And then in our study, the damped least square method which is often used to deal with singular matrices is applied to obtain Euler solution and we find that the damping factor, the eigenvalue of the normal equations matrix and even the relative position of the origin of coordinate could influence the accuracy of the Euler solutions. Specifically, the larger the damping factor, the closer the solution gets to zero, and the larger the size of the eigenvalue of the normal equations matrix, the smaller the impact of damping factor on solution, and also the matrix with small eigenvalue is sensitive to the change of the relative position of the coordinate origin. Accordingly, we propose a new method, called variable coordinate system selection method, to select best solutions. And in order to obtain more reliable solutions, an effective strategy based on the spatial relationship between the window center and the solution is used. Model test and real data application show that our selection methods are easy to perform and could dramatically reduce the number of spurious solutions.
Keywords: Euler deconvolution    Matrix singularity    Damped least square method    Change of coordinate system
0 引言

1 方法原理

 (1)

 (2)

 (3)

 (4)

 (5)

 (6)

(6) 式表明，对于任意的k>0, ‖P‖≠0总有‖P(k)‖ < ‖P(0)‖(Hoerl and Kennard, 1970)，因此相对于P而言，P(k)有向原点压缩的趋势，向原点压缩程度可由公式k(GTG+kI)-1‖表示.

 图 1 系数矩阵大小与反演结果向原点压缩程度分析图 (a)垂直磁化的球体模型，图中的1-6为滑动窗口中心点的位置; (b)滑动窗口中的系数矩阵最大特征值; (c)滑动窗口中的反演结果分布图(图中的数字与a图中的数字对应). Fig. 1 The figures used to analyse the relationship between the normal equations matrix and the accuracy of Euler solutions (a) Synthetic magnetic generated by 7 spheres, vertical magnetization, 1-6 represent the center of the sliding window; (b) The maximum eigenvalues obtained from six windows; (c) The inversion results obtained from the six windows.

 图 2 阻尼最小二乘反演结果随坐标系变化时的分布图 (a)坐标原点位于测区内时的反演结果; (b)坐标原点位于测区外时的反演结果. Fig. 2 The influence of the change of the relative position of the coordinate origin on Euler solutions obtained by damped least square method (DLS) (a) The distribution of the solutions when the coordinate origin is in the study area; (b) The distribution of the solutions when the coordinate origin is out of the study area.

 图 3 坐标系变化时的常规最小二乘反演结果分布图 (a)坐标原点位于测区内时的反演结果; (b)坐标原点位于测区以外时的反演结果. Fig. 3 The influence of the change of the relative position of the coordinate origin on Euler solutions obtained by ordinary least square method (OLS) (a) The distribution of the solutions when the coordinate origin is in the study area; (b) The distribution of the solutions when the coordinate origin is out of the study area.

 (7)

(7) 式表明，在某一次数据反演的窗口中，随着坐标原点远离测区，ΔP变大，P(k)偏离P的距离增大，P(k)越容易靠近坐标原点.又由于λ越小，反演结果向原点压缩程度越大.因此随着坐标原点远离测区，远离场源的窗口中产生的反演结果越容易靠近原点而偏离测区，此时测区内虚假反演结果的数量会急剧减少(图 4)，反演结果的发散性减小.

 图 4 坐标原点远离测区时，测区内反演结果的数量变化曲线图 Fig. 4 The number of solutions in study area with the coordinate origin away from the study area

 图 5 变坐标尺度下的反演结果分布情况对比图 (a)微小坐标尺度下的常规最小二乘反演结果；(b)微小坐标尺度下的阻尼最小二乘反演结果；(c)极大坐标尺度下的常规最小二乘反演结果；(d)极大坐标尺度下的阻尼最小二乘反演结果. Fig. 5 The distribution of the inversion results when the coordinate-scale changed (a) The results obtained by OLS when the coordinate-scale is tiny; (b) The inversion results obtained by DLS when the coordinate-scale is tiny; (c) The results obtained by OLS when the coordinate-scale is large; (d) The inversion results obtained by DLS when the coordinate-scale is large.

2 解的筛选

 (8)

 图 6 变坐标系筛选方法确定场源位置的过程及效果图 (a)坐标原点位于测区右上角时反演结果分布图; (b) a图中的反演结果对应的滑动窗口中心点位置; (c)坐标原点位于测区左下角时反演结果分布图; (d) c图中反演结果对应的滑动窗口中心点的位置; (e)将a、c图中相同滑动窗口中的两个反演结果取平均值得到的反演结果分布图; (f) e图中的反演结果对应的滑动窗口中心点位置. Fig. 6 The application of variable coordinate system selection method (a) The Euler solutions obtained by DLS when the origin of the coordinates is in the upper right corner; (b) The position of the sliding window corresponding to the results in (a); (c) The Euler solutions obtained by DLS when the origin of the coordinates is in the lower left corner; (d) The position of the sliding window corresponding to the results in (c); (e) Averaging the inversion results in the same window in (a) and (c); (f) The position of the sliding window corresponding to the results in (e).
3 模型实验 3.1 理想球体模型实验

 图 7 理想模型实验 (a)阻尼最小二乘反演结果分布图; (b)坐标原点位于测区右上角时的反演结果分布图; (c)坐标原点位于测区左下角时的反演结果分布图; (d)变坐标系筛选后的反演结果分布图. Fig. 7 An ideal model test (a) The inversion results generated by DLS; (b) The inversion results when the coordinate origin is in the upper right corner; (c) The inversion results when the coordinate origin is in the lower left corner; (d) The inversion results after selection.
3.2 复杂模型测试

 图 8 磁异常及变坐标系筛选方法确定的反演结果分布图 (a)组合模型磁异常分布图; (b)常规最小二乘反演结果分布图; (c)阻尼最小二乘反演结果分布图; (d)坐标原点位于测区左下角时反演结果分布图；(e)坐标原点位于测区右上角时反演结果分布图; (f)将c、d图中相同滑动窗口中的两个反演结果取平均值后得到的反演结果分布图; (g) d图中的反演结果对应的滑动窗口中心点; (h) e图中的反演结果对应的滑动窗口中心点; (i) f图中的反演结果对应的滑动窗口中心点. Fig. 8 Synthetic data produced by different geometries and the distribution of the Euler solutions (a) Model anomalous magnetic field; (b) The inversion results obtained by OLS; (c) The inversion results obtained by DLS; (d) The distribution of DLS Euler solutions when the coordinate origin is in the lower left corner; (e) The distribution of OLS Euler solutions when the coordinate origin is in the upper right corner; (f) Averaging the Euler solutions in the same sliding window in (c) and (d); (g) The position of the sliding window center corresponding to the results in (d); (h) The position of the sliding window center corresponding to the results in (e); (i) The position of the sliding window center corresponding to the results in (f).

 图 9 估算的场源位置分布图 (a)移动速率筛选后的反演结果水平分布图；(b) X-Z方向的反演结果分布图，红色方框为真实场源位置；(c)文献(Gerovska and Araúzo-Bravo, 2003)中的X-Z方向反演结果分布图，红色方框为真实场源位置；(d) Y-Z方向的反演结果分布图; (e)文献(Gerovska and Araúzo-Bravo, 2003)中的Y-Z方向反演结果分布图. Fig. 9 The distribution of the estimated depths (a) The estimated horizontal location of the source after selecting by moving-rate selection method; (b) The projection of the Euler solutions on X-Z plane by our method, the red rectangles respresent the position of the source; (c) The projection of the Euler solutions on X-Z plane by Gerovska and Araúzo-Bravo (2003), the red rectangles respresent the position of the source; (d) The projection of the Euler solutions on Y-Z plane by our method; (e) The projection of the Euler solutions on Y-Z plane by Gerovska and Araúzo-Bravo (2003).
4 实际应用

 图 10 湖南某实际铜镍矿勘探区地质图 Fig. 10 Geologic map of the copper-nickel deposit at Hunan Province
 图 11 实测数据等值线图及欧拉反演结果分布图 (a)该区域的矿区磁异常分布图; (b)常规最小二乘欧拉反演结果分布; (c)坐标原点位于测区左下角时反演结果分布; (d) c图中的反演结果对应的滑动窗口中心点位置; (e)坐标原点位于测区右上角时反演结果分布; (f) e图中的反演结果对应的滑动窗口中心点位置. Fig. 11 Magnetic total field data from Hunan province and the distribution of Euler solutions (a) Magnetic total field data; (b) The Euler solutions obtained by OLS; (c) The Euler solutions obtained by DLS when the coordinate origin is in the lower left corner; (d) The position of the sliding window center corresponding to the results in (c); (e) The Euler solutions obtained by DLS when the coordinate origin is in the upper right corner; (f) The position of the sliding window center corresponding to the results in (e).

 图 12 估算的场源位置分布图 (a)欧拉反演结果分布图, 颜色代表高程; (b) lp范数反演结果分布图，p=0; (c)图b中三条剖面的反演结果分布图，图中白色圆点为欧拉反演结果. Fig. 12 The distribution of the estimated depths (a) The distributon of the Euler solutions after selection, and the color represents the depth value; (b) Inversion results of the magnetic using lp-like-norm inversion, p=0; (c) The distribution of the estimated depths from three cross sections in (b), the white dots are the position of Euler solution.
5 结论

References
 Barbosa V C F, Silva J B C, Medeiros W E. 1999. Stability analysis and improvement of structural index estimation in Euler deconvolution. Geophysics, 64(1): 48-60. DOI:10.1190/1.1444529 Beiki M. 2013. TSVD analysis of Euler deconvolution to improve estimating magnetic source parameters: An example from the Åsele area, Sweden. Journal of Applied Geophysics, 90: 82-91. DOI:10.1016/j.jappgeo.2013.01.002 Davis K, Li Y G, Nabighian M. 2010. Automatic detection of UXO magnetic anomalies using extended Euler deconvolution. Geophysics, 75(3): G13-G20. DOI:10.1190/1.3375235 Dewangan P, Ramprasad T, Ramana M V, et al. 2007. Automatic interpretation of magnetic data using Euler deconvolution with nonlinear background. Pure and Applied Geophysics, 164(11): 2359-2372. DOI:10.1007/s00024-007-0264-x Fan M N. 2006. The study and application of Euler deconvolution [Ph. D. thesis] (in Chinese). Changchun: Jilin University. Fedi M. 2016. An unambiguous definition of the structural index. // Society of Exploration Geophysicists. SEG Technical Program Expanded Abstracts 2016. 1537-1541. Florio G, Fedi M, Pateka R. 2014. On the estimation of the structural index from low-pass filtered magnetic data. Geophysics, 79(6): J67-J80. DOI:10.1190/geo2013-0421.1 Gerovska D, Araúzo-Bravo M J. 2003. Automatic interpretation of magnetic data based on Euler deconvolution with unprescribed structural index. Computers & Geosciences, 29(8): 949-960. Gerovska D, Stavrev P, Araúzo-Bravo M J. 2005. Finite-difference Euler deconvolution algorithm applied to the interpretation of magnetic data from northern Bulgaria. Pure and Applied Geophysics, 162(3): 591-608. DOI:10.1007/s00024-004-2623-1 Hansen R O, Suciu L. 2002. Multiple-source Euler deconvolution. Geophysics, 67(67): 525-535. Hoerl A E, Kennard R W. 1970. Ridge regression: biased estimation for nonorthogonal problems. Technometrics, 12(1): 55-67. DOI:10.1080/00401706.1970.10488634 Li Z, Yao C, Zheng Y, et al. 2018. 3D magnetic sparse inversion using an interior-point method. Geophysics, 83(3): J15-J32. DOI:10.1190/geo2016-0652.1 Liu Q, Yao C, Zheng Y. 2019. A new method to select the best Euler solutions based on the angle variation. 81st EAGE Conference and Exhibition 2019. Marquardt D W. 1963. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics, 11(2): 431-441. DOI:10.1137/0111030 Martelet G, Sailhac P, Moreau F, et al. 2001. Characterization of geological boundaries using 1-D wavelet transform on gravity data: Theory and application to the Himalayas. Geophysics, 66: 1116-1129. DOI:10.1190/1.1487060 Melo F F, Barbosa V C. 2018. Correct structural index in Euler deconvolution via base-level estimates. Geophysics, 83(6): J87-J98. DOI:10.1190/geo2017-0774.1 Melo F F, Barbosa V C, Uieda L, et al. 2013. Estimating the nature and the horizontal and vertical positions of 3D magnetic sources using Euler deconvolution. Geophysics, 78(6): J87-J98. DOI:10.1190/geo2012-0515.1 Mushayandebvu M F, Lesur V, Reid A B, et al. 2004. Grid Euler deconvolution with constraints for 2D structures. Geophysics, 69(2): 489-496. DOI:10.1190/1.1707069 Nabighian M N, Hansen R O. 2001. Unification of Euler and Werner deconvolution in three dimensions via the generalized Hilbert transform. Geophysics, 66(6): 1805-1810. DOI:10.1190/1.1487122 Reid A B, Allsop J M, Granser H, et al. 1990. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics, 55(1): 80-91. DOI:10.1190/1.1442774 Reid A B, Ebbing J, Webb S J. 2014. Avoidable Euler errors-the use and abuse of Euler deconvolution applied to potential fields. Geophysical Prospecting, 62(5): 1162-1168. DOI:10.1111/1365-2478.12119 Reid A B, Thurston J B. 2014. The structural index in gravity and magnetic interpretation: Errors, uses, and abuses. Geophysics, 79(4): J61-J66. DOI:10.1190/geo2013-0235.1 Salem A, Williams S, Fairhead D, et al. 2007. Interpretation of magnetic data using tilt-angle derivatives. Geophysics, 73(1): L1-L10. Silva J B, Barbosa V C. 2003. 3D Euler deconvolution: Theoretical basis for automatically selecting good solutions. Geophysics, 68(6): 1962-1968. DOI:10.1190/1.1635050 Stavrev P Y. 1997. Euler deconvolution using differential similarity transformations of gravity or magnetic anomalies. Geophysical Prospecting, 45(2): 207-246. DOI:10.1046/j.1365-2478.1997.00331.x Stavrev P, Reid A. 2006. Degrees of homogeneity of potential fields and structural indices of Euler deconvolution. Geophysics, 72(1): L1-L12. Thompson D T. 1982. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data. Geophysics, 47(1): 31. DOI:10.1190/1.1441278 Wang Y G. 2011. Two-dimensional forward calculation of magnetic terrain fluctuation anomaly and its application in Changjie copper and nickel deposit exploration in Tongdao County of Hunan Province [Master's thesis] (in Chinese with English abstract). Beijing: China University of Geosciences. Yao C L, Guan Z N, Wu Q B, et al. 2004. An analysis of Euler deconvolution and its improvement. Geophysical and Geochemical Exploration (in Chinese), 28(2): 150-155. 范美宁. 2006.欧拉反褶积方法的研究及应用[博士论文].长春: 吉林大学. http://cdmd.cnki.com.cn/article/cdmd-10183-2006109772.htm 王彦国. 2011.二维起伏地形磁异常正演及其在湖南通道长界铜镍矿床勘探中的应用[硕士论文].北京: 中国地质大学(北京). http://cdmd.cnki.com.cn/article/cdmd-11415-1011078029.htm 姚长利, 管志宁, 吴其斌, 等. 2004. 欧拉反演方法分析及实用技术改进. 物探与化探, 28(2): 150-155. DOI:10.3969/j.issn.1000-8918.2004.02.017