-
近年来,磁力勘探主要应用在军事侦察、水文及工程地质勘探等领域。运用磁探测技术可以对地下或者水下小尺度磁性目标进行定位与识别[1],为下一步的工作提供数据支持。反演方法是解释实测数据的重要一步,其主要目的是利用地表磁测数据来评估地下未知磁异常体的轮廓形态[2]。在磁测数据线性反演方法中,假设观测面以下的实测区域可以离散化为具有固定物性参数值的规则立方体模型,由于地下待求未知参数远远多于观测数据点数[3],因此会导致线性方程求解的模糊性,难以得到方程的最优解。测量数据的噪声和固有的场源不唯一性,导致磁测数据的反演是病态多解的。为了降低测量误差等因素的影响,常常使用L2或者Tikhonov正则化思想[4, 5]。正则化思想的两个重要的方面是模型约束器和正则化参数的选择,前者可以使数据获得稳定解,后者是数据项和正则项之间的折中系数。在反演问题中正则化参数的选择方法是多年来研究的热点问题之一。确定正则化参数的主要方法有LC曲线(L-curve, LC)、广义交叉验证(generalized cross validation, GCV)及莫罗佐夫偏差原则(Morozov discrepancy principle, MDP)[6-8]。其中LC曲线法计算量较大,需要计算多个正则化参数来确定最优解; GCV方法难以保证在在每次迭代过程中都能获得GCV功能矩阵的极小值;相比前两种方法,MDP目前应用较为广泛,但是在求解过程中会产生过高的正则化参数值,导致结果过于光滑[9-11]。
针对上述反演过程存在的问题,本文提出了基于χ2准则的磁梯度张量3D聚焦反演方法。通过向目标函数中加入深度加权矩阵和最小支撑矩阵,避免了反演解的不稳定性;采用奇异值分解法将反演解转化到合理的值域范围内;利用χ2准则对正则化参数进行自适应选择,提高了迭代速度和求解精度。
-
磁梯度张量正演公式如下:
$$ \boldsymbol{d} = \boldsymbol{Am},p \ll n $$ (1) 式中,核函数矩阵A的元素Aij为第i(i=1, 2, 3…p)个观测点观测的第j(j=1, 2, 3…n)个单元的响应;向量d∈Rp为观测数据矩阵;向量m∈Rn为待测模型参数矩阵。在磁梯度张量数据实际测量中,由于操作误差等影响导致实测数据中包含一定的噪声。设e∈Rp代表实测数据中的误差,则式(1)表示为:
$$ {\boldsymbol{d}_{\rm obs}} = \boldsymbol{Am} + \boldsymbol{e} $$ (2) 根据式(2)可知,磁梯度张量反演的过程就是利用包含噪声的实际测量数据dobs求解模型参数m。式(2)表示的病态反问题可以表示为Tikhonov最小二乘公式:
$$\begin{align} \boldsymbol{m}\left( \alpha \right) = &\arg \mathop {\min }\limits_m \{ {\left\| {{\boldsymbol{W}_d}\left( {\boldsymbol{Am} - {\boldsymbol{d}_{\rm obs}}} \right)} \right\|}^2_2 + \\ & {\alpha ^2}\left\| {\boldsymbol{D}(\boldsymbol{m} - {\boldsymbol{m}_{\rm apr}})} \right\|^2_2 \}\end{align} $$ (3) 式中,‖Wd(Am-dobs)‖22为加权数据保真度,‖D(m-mapr)‖22为正则项,Am为预测数据向量,Wd=diag(1/η1, …, 1/ηp)为数据加权矩阵,ηi为第i个基准噪声的标准偏差,D为正则化矩阵,mapr为模型参数先验信息向量,α为正则化参数。设Ã=WdA,$\boldsymbol{\tilde d} $obs=Wddobs,Z=m-mapr。转换后式(3)表示为:
$$ \begin{array}{c} \boldsymbol{Z}\left( \alpha \right) = \arg \mathop {\min }\limits_y \left\{ {{\left\| {\boldsymbol{\tilde A}\boldsymbol{m} - \boldsymbol{\tilde r}} \right\|}^2_2 + {\alpha ^2}{\left\| {\boldsymbol{DZ}} \right\|}^2_2} \right\},\\ \boldsymbol{\tilde r }= ({{\boldsymbol{\tilde d}}_{\rm obs}} - \boldsymbol{\tilde A}{\boldsymbol{m}_{\rm apr}}) \end{array} $$ (4) 此时m(α)仅由正则化参数α决定,用矩阵Ã(α)表示为:
$$ \begin{align} \boldsymbol{Z}\left( \alpha \right) &= {({\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A} + {\alpha ^2}{\boldsymbol{D}^{\rm T}}\boldsymbol{D})^{ - 1}}{\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde r} = \boldsymbol{\tilde A}\left( \alpha \right)\boldsymbol{\tilde r},\\ &\boldsymbol{\tilde A}\left( \alpha \right) = {({\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A} + {\alpha ^2}{\boldsymbol{D}^{\rm T}}\boldsymbol{D})^{ - 1}}{\boldsymbol{\tilde A}^{\rm T}} \end{align} $$ (5) $$ \boldsymbol{m}\left( \alpha \right) = {\boldsymbol{m}_{\rm apr}} + \boldsymbol{Z}\left( \alpha \right) = {\boldsymbol{m}_{\rm apr}} + \boldsymbol{\tilde A}\left( \alpha \right)\boldsymbol{\tilde r} $$ (6) Hansen指出矩阵D为可逆矩阵[10],因此:
$$ \begin{array}{c} ({\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A} + {\alpha ^2}{\boldsymbol{D}^T}\boldsymbol{D}) = \\ {\boldsymbol{D}^T}({({\boldsymbol{D}^T})^{ - 1}}{\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A}{\boldsymbol{D}^{ - 1}} + {\alpha ^2}{I_n})\boldsymbol{D} \end{array} $$ (7) 定义右预处理矩阵$ \boldsymbol{\mathop A}\limits^ \approx$=ÃD-1,正则化反演矩阵$ \boldsymbol{\mathop A}\limits^ \approx(\alpha)=({\boldsymbol{\mathop A}\limits^ \approx }^{\rm T}\boldsymbol{\mathop A}\limits^ \approx+\alpha^2I_n)^{-1}{\boldsymbol{\mathop A}\limits^ \approx}^{\rm T}$,此时J(α)=DZ(α),则式(4)表示为:
$$ \boldsymbol{J}\left( \alpha \right) = \arg \mathop {\min }\limits_J \left\{ {{\left\| {\boldsymbol{\mathop A}\limits^ \approx \boldsymbol{J} - \boldsymbol{\tilde r}} \right\|}^2_2 + {\alpha ^2}{\left\| \boldsymbol{J} \right\|}^2_2} \right\} $$ (8) $$ \boldsymbol{m}\left( \alpha \right) = {\boldsymbol{m}_{\rm apr}} + {\boldsymbol{D}^{ - 1}}\boldsymbol{J}\left( \alpha \right) $$ (9) 正则化矩阵D(k)∈Rn×n可以表示为D(k)=We(k)Wdepth。文献[6]给出深度加权矩阵Wdepth=diag(1/(zj+ξ)β)来抵消核函数随深度增加所带来的衰减,其中zj为地下第j个剖分单元的中心埋深,ξ>0是为了避免出现奇点而引入的数值,β是可以控制加权函数大小的正数。Vatankhah指出,利用最小支撑矩阵We(k)=diag((m(k)-m(k-1))+ε2)-1/2获得的反演结果具有较好的非光滑特征,能够提高块状异常体的分辨率,k>0,W(0)=I, m(0)=mapr,ε>0为聚焦参数。通过正则化矩阵对正则项进行约束,减小模型解的多解性,提高了反演解的稳定性。在反演过程中,为了得到与实际情况更加吻合的物性参数分布,利用惩罚函数在反演迭代过程中去除病态解,在迭代过程中设置物性参数上下限约束函数[mmin, mmax],对物性参数进行约束使反演结果控制在合理的范围内:
$$ \left\{ \begin{array}{c} m\left( j \right) = {m_{\max}},若m\left( j \right) > {m_{\max}}\\ m\left( j \right) = {m_{\min}},若m\left( j \right) < {m_{\min}}\\ j = 1,2,3, \ldots ,n; \end{array} \right.{\rm{ }} $$ (10) 设置终止迭代条件:当χComputed2=‖(dobs)i-(dpre)i/ηi‖22≤p+$\sqrt {2p} $或达到最大迭代次数k时迭代终止。迭代过程为:
$$ \boldsymbol{J}({\alpha ^{(k + 1)}}) = ({\boldsymbol{\mathop A}\limits^ \approx}^{\rm T}\boldsymbol{\mathop A}\limits^ \approx + {({\alpha ^{(k)}})^2}{\boldsymbol{I}_n})^{ - 1}\boldsymbol{\mathop A}\limits^ \approx{\boldsymbol{\tilde r}^{(k)}} $$ (11) $$ {\boldsymbol{m}^{(k + 1)}} = {\boldsymbol{m}^{(k)}} + {({\boldsymbol{D}^{\left( {k + 1} \right)}})^{ - 1}}\boldsymbol{J}({\alpha ^{(k + 1)}}) $$ (12) 在求解J(α)时,对右处理矩阵$ \boldsymbol{\mathop A}\limits^ \approx$∈Rp×n进行奇异值分解(singular value decomposition,SVD)。矩阵$ \boldsymbol{\mathop A}\limits^ \approx$可以表示为$ \boldsymbol{\mathop A}\limits^ \approx$=U∑VT,其中,U、V分别为p×p和n×n的正交矩阵。∑=diag([δ1, δ2, …δp]),其中奇异值δ1≥δ2≥…≥δp>0。为了防止高频观测误差被放大,引入正则化滤波因子fi(α)=$\frac{{\delta ^2_i}}{{\delta ^2_i + {\alpha ^2}}} $, 1≤i≤p。因此,对式(11)进行奇异值分解得:
$$ \boldsymbol{J}\left( \alpha \right) = \sum\limits_{i = 1}^p {\frac{{\delta ^2_i}}{{\delta ^2_i + {\alpha ^2}}}\frac{{\boldsymbol{u}^{\rm T}_i\boldsymbol{\tilde r}}}{{{\delta _i}}}} {\boldsymbol{v}_i} = \sum\limits_{i = 1}^p {{f_i}\left( \alpha \right)\frac{{\boldsymbol{u}^{\rm T}_i\boldsymbol{\tilde r}}}{{{\delta _i}}}{\boldsymbol{v}_i}} $$ (13) 式中,ui和vi分别为矩阵U和V的第i列。
-
徐新禹等讨论了在Tikhonov正则化方法中如何有效的选择正则化参数[12, 13]。Zhdanov提出一种自适应正则化参数的选择方法,定义α(k)=α(1)q(k),q为衰减因子[14]。但是这种方法不能保证每次迭代得到的正则化参数为最佳解,影响求解精度。还有一种常用的方法是MDP法,但Kilmer指出该方法会导致反演结果过度光滑[15]。因此,本文提出基于χ2准则的正则化参数选择方法。
χ2准则实际上是MDP方法的一种扩展算法。Mead指出MDP是基于后验策略的偏差原理[11],通过利用先验误差水平信息,使得实测数据与预测数据的残差为最小值的正则化参数即为最优正则化参数,MDP实际上就是遵循自由度为p-n的χ2分布,可以表示为:
$$ \left\| {\boldsymbol{\mathop A}\limits^ \approx \boldsymbol{J}\left( \alpha \right) - \boldsymbol{\tilde r}} \right\|^2_2 = p - n,p \ge n $$ (14) 由式(14)可知,当p<n时,上述方法是不可用的。而基于χ2准则的正则化参数确定方法是通过确定Tikhonov正则化目标函数的极小值来确定最优正则化参数,且无需误差先验信息。该方程遵循自由度为p的χ2分布,可以表示为:
$$ \left\| {{\boldsymbol{W}_d}\left( {\boldsymbol{Am} - {\boldsymbol{d}_{\rm obs}}} \right)} \right\|^2_2 + {\alpha ^2}\left\| {\boldsymbol{D}\left( {\boldsymbol{m} - {\boldsymbol{m}_{\rm apr}}} \right)} \right\|^2_2 = p $$ (15) 定义一个单直立长方体理论模型,纵横向切片图如图 1(a)、(b)所示。将地下待测空间划分为22×22×10=4 840个单元格,每个单元格均为边长为0.1 m的正方体。在仅考虑感应磁化条件下,异常体的磁化强度为20 A/m,磁倾角为55°,磁偏角为-18°。根据d=Am计算右端项的理论值,并加入0.01×std(d, 1)×rand(484, 1)的误差作为实测数据。定义初始正则化参数α(1)=(n/p)γ(max(δi)/mean(δi)),其中0≤γ≤2,δi为对应核函数矩阵$\boldsymbol{\mathop A}\limits^ \approx $的奇异值,定义模型重构相对误差rms=‖mexact-m(k)‖2/‖mexact‖2,当χComputed2≤506或最高迭代次数K=100时迭代终止。本文仅选择垂直向和水平向分辨率相对较高的Bzz分量进行方法的有效性验证[16-17], 反演结果如图 1(c)~1(f)所示。
从反演结果可以看出,两种方法都使得反演结果具有较好的横向和纵向分辨率,但是从表 1可知,基于χ2准则的正则化参数确定方法比MDP迭代次数少,求得的反演结果精度更高。由图 2可以看出,正则化参数初始值很大,表明在迭代初期正则项起到较大的约束作用,随后迅速衰减,最后随迭代次数增加趋于稳定。
表 1 不同迭代原则下各部分参数值
Table 1. Parameter Values at Different Iterations Principle
迭代原则 α(1), γ=2 相对误差 迭代次数 χ2准则 112.6 0.7018 6 MDP 112.6 0.7096 20 -
将地下待测空间划分为22×22×10=4 840个单元格,每个单元格均为边长为0.1 m的正方体,地下正演模型如图 3所示。在仅考虑感应磁化的条件下,磁倾角为55°,磁偏角为-16°,3个不同位置、不同形态的异常体的具体位置及参数如表 2。如图 4所示,在反演时加入0.001×std(d, 1)×rand(484, 1)的误差作为实测数据,并利用Bzz分量进行反演。反演过程中设置最高迭代次数为15次,设置物性参数阈值为m(k)>8A/m,所得反演结果的空间分布如图 5所示。从图中结果可以看出:该反演方法能够较好的反映出不同位置、不同物性参数模型的轮廓形态,具有较好的横向分辨率,纵向上随深度增加模型分辨率有所下降,存在边缘效应,这是由于重构模型的分辨率不可避免的随深度增加而下降,但是能够大致分布在相应的位置。
表 2 异常体参数
Table 2. Abnormal Body Parameters
异常体 X方向位置/m Y方向位置/m Z方向位置/m 磁化强度M(A·m-1) 异常体1 1.55~1.85 0.35~1.75 0.1~0.4 20 异常体2 0.25~0.55 0.75~1.75 0.3~0.6 30 异常体3 1.35~1.75 0.25~0.95 0.2~0.4 40 -
测区位于石家庄某地,测区内放置两个长方体铁块。在地理坐标系下铁块1位于测区东南方向,中心坐标为(1.55, 0.5, 0.16),块体尺寸为0.2 m×0.4 m×0.05 m,铁块2位于测区西北方向,中心坐标为(0.65, 1.2, 0.27),块体尺寸为0.2 m×0.4 m×0.1 m。实验装置为十字形传感器阵列结构,利用Bartington公司的三轴磁通门传感器搭建磁梯度张量探头。图 6(a)展示了测量装置及铁块的实际位置分布,图 6(b)为Bzz分量实际测量数据。将待测空间划分为22×22×10=4 840个单元格,每个单元格均为边长为0.1 m×0.1 m×0.1 m的长方体,地面的观测点为22×22=484个网格。分别采用MDP准则和χ2准则对磁异常进行反演,得到的反演结果如图 7所示。图 7(a)、7(c)为模型体反演结果平面图,图 7(b)、7(d)为物性参数阈值m(k)>5 A/m时的三维反演结果。其中利用χ2准则进行计算需要迭代12次,而利用MDP方法进行迭代需要17次。根据反演结果可以看出,相比MDP准则,χ2准则的反演结果的横向和纵向模型分辨率较好,能够更好地确定异常体的轮廓及位置,进一步证明了本文方法的有效性。
-
1) 利用深度加权矩阵和最小支撑矩阵对经典Tikhonov正则化理论框架下的反演模型进行约束,提高了反演解的稳定性;
2) 采用奇异值分解法(SVD)对核函数矩阵进行分解,利用物性参数上下限函数将反演解转换到合理的值域范围内;
3) 利用χ2准则对迭代过程中的正则化参数进行选择,实现了正则化参数选择的自适应功能,提高了迭代速度和求解精度;
4) 虽然本文方法能够使获得的物性参数值分布在相应的位置,但是重构模型的分辨率不可避免地随着深度的增加而下降,存在一定的边缘效应。
Three-dimensional Focusing Inversion of Magnetic Gradient Tensor Data Based on the χ2 Principle
-
摘要: 针对铁磁性物质反演中正则化参数自适应选择的问题,提出了基于χ2准则的磁梯度张量3D聚焦反演方法。利用深度加权矩阵和最小支撑矩阵对经典Tikhonov正则化理论框架下的反演模型进行约束得到目标函数,避免了由于反演参数多于采集点数而导致反演解的多解性,并有效解决了核函数随深度增大而快速衰减的问题。通过对目标函数进行迭代奇异值分解获得最佳物性参数,并根据χ2准则自适应地确定目标函数在迭代过程中的正则化参数,提高了迭代速度和求解精度。仿真和实验结果表明:该方法能准确还原磁性异常体的轮廓形态,具有较好的模型分辨率。Abstract: The method of magnetic gradient tensor three-dimensional pose inversion is researched. Three-dimensional focusing inversion of magnetic gradient tensor data based on the χ2 principle was proposed for regularization parameter adaptively selected issues of small target ferromagnetic materials in the inversion process. Depth weighting matrix and minmal support matrix were added in the inversion model of Tikhonov regularization theoretical framework to obtain the objective function to avoid the multiple solutions of inversion problem caused by the number of inverse parameter far than observational points and compensate kernel function avoid diminshing rapidly with depth. This function was calculated the optimal physical parameters by singular value decomposition. and according the χ2 principle to adaptively determine the regularization parameter in the iterative process, thereby this principle increased the iterative speed and the solver accuracy. Three-dimensional focusing inversion of magnetic gradient tensor data produces models with non-smooth properties, for which typical implementations in this field use the iterative minmum support stabilizer and both regularizer and regularizing parameter are updated each iteration. The χ2 principle generalizes the Morozov discrepancy principle to the augmented residual of the Tikhonov regularized least squares problem. For weighting of the data fidelity by a known Gaussian noise distribution on the measured data and, when the stabilizing, or regularization, term is considered to be weighted by unknown inverse covariance information on the model parameters, the minmum of the Tikhonov functional becomes a random variable that follows a χ2-distribution. Simulation and experimental results show that:this inversion method can reflect the outline shade of the magnetic anomaly and has good model resolution, which provide the theoretical foundation and practical experience for the application of the inversion method of magnetic gradient tensor.
-
表 1 不同迭代原则下各部分参数值
Table 1. Parameter Values at Different Iterations Principle
迭代原则 α(1), γ=2 相对误差 迭代次数 χ2准则 112.6 0.7018 6 MDP 112.6 0.7096 20 表 2 异常体参数
Table 2. Abnormal Body Parameters
异常体 X方向位置/m Y方向位置/m Z方向位置/m 磁化强度M(A·m-1) 异常体1 1.55~1.85 0.35~1.75 0.1~0.4 20 异常体2 0.25~0.55 0.75~1.75 0.3~0.6 30 异常体3 1.35~1.75 0.25~0.95 0.2~0.4 40 -
[1] 卞光浪, 翟国军, 黄谟涛.顾及地磁背景场的多目标磁异常分量换算方法[J].武汉大学学报·信息科学版, 2011, 36(8):914-918 http://ch.whu.edu.cn/CN/abstract/abstract614.shtml Bian Guanglang, Zhai Guojun, Huang Motao. Transformation of Multiobjects Magnetic Anomaly Components with Geomagnetic Effect[J]. Geomatics and Information Science of Wuhan University, 2011, 36(8):914-918 http://ch.whu.edu.cn/CN/abstract/abstract614.shtml [2] 陈召曦, 孟小红, 郭良辉.重磁数据三维物性反演方法进展[J].地球物理学进展, 2012, 27(2):503-511 http://www.doc88.com/p-5146290251314.html Chen Zhaoxi, Meng Xiaohong, Guo Lianghui. Review of 3D Property Inversion of Gravity and Magnetic Data[J]. Progress in Geophysics, 2012, 27(2):503-511 http://www.doc88.com/p-5146290251314.html [3] Boulanger O, Chouteau M. Constraint in 3D Gravity Inversion[J]. Geophysical Prospecting, 2001, 49(2):265-280 doi: 10.1046/j.1365-2478.2001.00254.x [4] Aster R C, Borchers B, Thurber C H. Parameter Estimation and Inverse Problems[M]. Amsterdam:Second Edition Elsevier Inc, 2011 [5] Karaoulis M, Revil A, Minsley B, et al. Timelapse Gravity Inversion with an Active time Constraint[J]. Geophys.J.Int, 2013, 196(2):748-759 https://academic.oup.com/gji/article/196/2/748/572578/Time-lapse-gravity-inversion-with-an-active-time [6] Li Y, Oldenburg D W. 3-D Inversion of Magnetic data[J]. Geophysics, 1996, 61(2):394-408 doi: 10.1190/1.1443968 [7] Farquharson C G, Oldenburg D W. A Comparison of Automatic Techniques for Estimating the Regularization Parameter in Non-linear Inverse Problems[J]. Geophys.J.Int. 2004, 156(3):411-425 doi: 10.1111/gji.2004.156.issue-3 [8] Vatankhah S, Ardestani V E, Renaut R A. Automatic Estimation of the Regularization Parameter in 2-D Focusing Gravity Inversion:Application of the Method to the Safo Manganese Mine in Northwest of Iran[J]. Journal of Geophysics and Engineering, 2014, 11(4):045001 doi: 10.1088/1742-2132/11/4/045001 [9] Wahba G. Practical Approximate Solutions to Linear Operator Equations when the Data are Noisy[J]. SIAM Journal on Number Anal, 1977, 14(4):651-667 doi: 10.1137/0714044 [10] Hansen P C. Analysis of Discrete Ill-posed Problems by Means of the L-curve[J]. Society for Industrial and Applied Mathematics, 1992, 34(4):561-579 https://mathscinet.ams.org/mathscinet-getitem?mr=1193012 [11] Mead J L. Discontinuous Parameter Estimates with least Squares Estimators[J]. Applied Mathematics and Computation, 2013, 219(10):5210-5223 doi: 10.1016/j.amc.2012.11.067 [12] 徐新禹, 李建成, 王正涛, 等. Tikhonov正则化方法在GOCE重力场求解中的模拟研究[J].测绘学报, 2010, 39(5):465-470 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=chxb201005007&dbname=CJFD&dbcode=CJFQ Xu Xinyu, Li Jiancheng, Wang Zhengtao, et al. The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Msission[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(5):465-470 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=chxb201005007&dbname=CJFD&dbcode=CJFQ [13] 吴小平.利用共轭梯度法的电阻率三维反演若干问题研究[J].地震地质, 2001, 23(12):321-327 http://www.cqvip.com/qk/95728X/200102/5259337.html Wu Xiaoping. Study on Some Problems for 3D Resistivity Inversion Using Conjugate Gradient[J]. Seismology and Geology, 2001, 23(12):321-327 http://www.cqvip.com/qk/95728X/200102/5259337.html [14] Zhdanov M S, Tartaras E.Three-dimensional Inversion of Multitransmitter Electromagnetic data Based on the Localized Quasi-linear Approximation[J]. Geophys. J.Int, 2002, 148(3):506-519 https://www.researchgate.net/publication/229735175_Three-dimensional_inversion_of_multitransmitter_electromagnetic_data_based_on_the_localized_quasi-linear_approximation [15] Kilmer M E, O'Leary D P. Choosing Regularization Parameters in Iterative Methods for Ill-posed Problems[J]. SIAM Journal on Matrix Analysis and Applications, 2001, 22(4):1204-1221 doi: 10.1137/S0895479899345960 [16] 陈少华, 朱自强, 鲁光银, 等.重力梯度张量的预条件共轭梯度法反演[J].中南大学学报(自然科学版), 2013, 44(2):619-625 http://www.doc88.com/p-0387192660080.html Chen Shaohua, Zhu Ziqiang, Lu Guangyin, et al. Inversion of Gravity Gradient Tensor Based on Preconditioned Conjugate Gradient[J]. Journal of Central South University(Science and Technology), 2013, 44(2):619-625 http://www.doc88.com/p-0387192660080.html [17] 杨娇娇, 刘展, 陈晓红, 等.基于深度加权的重力梯度张量数据的3D聚焦反演[J].世界地质, 2015, 34(1):210-218 doi: 10.3969/j.issn.1004-5589.2015.01.025 Yang Jiaojiao, Liu Zhan, Chen Xiaohong, et al. Three-dimensional Focusing Inversion of Gravity Gradient Tensor Data Based on Depth Weighting[J]. Global Geology, 2015, 34(1):210-218 doi: 10.3969/j.issn.1004-5589.2015.01.025 -