文章编号：1004-0609(2014)10-2601-08
基于混合正则化的重力场约束反演
朱自强，曹书锦，鲁光银
(中南大学 地球科学与信息物理学院，长沙 410083)
摘 要：在实际地下地质构造是一类多尺度的构造(如断层和褶皱等)，而传统的正则化方法多基于最小光滑策略，其反演密度模型一般不易辨识以上构造。为此在分裂Bregman迭代正则化框架下引入混合正则化方法以充分利用非光滑反演和小波多尺度反演算子的特性，引入与衰减系数无关的深度加权矩阵以更好地描述深部异常；针对非光滑反演中异常幅值易于超出现实及理论异常范围，引入密度成像中的约束以确保反演具有物理意义。通过设置两类模型，对比多类正则化反演方法。反演结果显示：混合正则化反演能有效地勾勒异常边界；在处理埋深不同的异常源时，相对于聚焦反演出现的过度聚焦现象而导致的反演深度描述不准确、异常歪斜，混合正则化反演的聚焦效应相对较弱、但深度描述准确。这表明本研究反演确实可行、有效，且具有更强的适应性。
关键词：深度加权；密度约束；聚焦反演；混合正则化；边界识别
中图分类号：P631.1 文献标志码：A
3D gravity inversion with bound constraint based on hyper-parameter regularization
ZHU Zi-qiang, CAO Shu-jin, LU Guang-yin
(School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)
Abstract: Traditional regularization inversions based on minimum smooth cannot distinguish the discontinuity of the underground subsurface (such as faults and folds, etc.). Hyper-parameter regularization method was introduced under split Bregman iterative regularization framework for taking advantage of edge-preserving inversion and wavelet multiscale operator. A new depth matrix was created based on sensitive matrix for giving an exact description of deep density abnormity. Physical bounds were imposed in each inverse iteration to obtain meaningful solutions and to avoid insignificant abnormality in non-smooth inversion. The smoothness inversion, Marquardt inversion, Occam inversion and focusing inversion with two synthetic models were compared. The inversion results show that the hyper-parameter regularization inversion can preserve edge effectively and has relatively weak focusing effect when treating the designed model with different depth sources. Meanwhile, this method can avoid inaccurate depth description and obliqueness caused by over-regularization in focusing inversion. Moreover, the inversion proposed in this study is feasible, effective, and has better adaptability.
Key words: depth matrix; bound constraint; focusing inversion; hyper-parameter regularization; edge-preserving
传统Tikhonov正则化方法是针对线性不适定的问题而提出的，其具有成像速度快、求解稳定、定位准确等优点。然而现实的地质构造可抽象为两部分，其一为地层内物性变化不大的部分，即平滑部分；其二为地层与地层的接触部分，即间断部分[1]。这导致了基于全局平滑策略采用连续函数作罚函数的Tikhonov正则化算法会因L2范数函数平滑效应在一定程度上存在反演质量不高(对比度低，分辨能力较差，多个物体之间存在伪迹等)的弊病[2]。如何得到更清晰的地质界面的反演结果仍是当前地球物理反演研究的一个重点课题[3]。
为此，在过去的十几年中，众多学者在边界识别/边界保存的位场反演技术等方面做出了大量的工作，大体上可以分为如下两类：非光滑反演和多尺度小波域反演。非光滑反演如全变差正则化[4-5]、聚焦反演[6-7]和尖锐边界/陡边界反演[8-9]。尖边界反演在地质先验信息充足的情况下能取得极好的反演效果，但其对初始模型的设置和先验地质信息的准确性要求较高，难适用于起伏地形条件下的反演工作，且存在反演深部构造不够准确以及加入约束条件较困难等不足。聚焦反演虽具有成像快速和对单一异常效果明显等优点，但受聚焦作用的影响，在地质构造复杂且各构造权重不均的情况下，易使得构造形态发生畸变及导致对目标区域的定位不准[2]。与传统正则化方法相比，全变差正则化建立在泛函分析和微分几何的数学模型基础上，所用罚函数的全变差滤波器结构简单，系数容易计算。它对离散函数波动的阻尼作用是离散的，这就使得变差函数很适用于离散分布函数(介质非连续分布区域)的正则化处理[2]，从而消除噪声并保持反演细节，特别是具有较好的保边缘性，且成像速度较快[5, 10]。但是全变差正则化的目标函数收敛程度不如Tikhonov正则化好，存在定位目标区域不准确等缺点[11]。实际上，每种正则化方法都有一定的侧重性，也就意味着对相应的问题有着非常重要的单一优势，而在面对较为复杂的地质构造时，这种优势即转变成过度正则化效应，如Tikhonov正则化在地质构造边界上的过度光滑、全变差正则化只能较好的适用于连续常值区域即地质构造内等[11]。而第二类反演策略，将模型通过基函数投影到不同尺度的空间，并采用一定阈值的罚函数剔除噪声/保留某一特定频率[12-13]，以达到保持不连续边界，获取更强的边界识别能力[14-15]。
为此，本文作者将全变差正则化与多尺度小波进行优化组合，以在利用它们对特定构造信息有某些优点的同时规避其本身固有的缺陷，从而更好地在得到清晰构造边界的同时快速而稳定的求解。在本研究中，引入一类多尺度小波稀疏算子，以借助其在保持边界的同时具有一定光滑效应，在此基础上引入具有边界保持特性的全变差算子，将两者在分裂Bregman迭代算法的框架下有机地实现混合正则化反演[16]；引入与衰减系数
无关的深度加权矩阵，以更好地描述深部异常；针对非光滑反演中，异常幅值易于超出在现实及理论异常范围，引入异常幅值上下限以确保反演具有物理意义。本研究设置两类不同模型(同一深度和一系列多个不同深度)，分别对比多类位场反演及聚焦反演与本研究混合正则化反演在面对不同深度异常的适应性。面对同一深度异常时，由于空间内各异常源权重大小相对平衡，反演深度加权仅表现在所得反演结果是否能与初始正演模型一致。数值结果表明，光滑反演、Marquardt反演[17]和Occam反演[18]均表现为能准确确定异常大体空间位置，但以上三者均难以分辨其异常下底板位置；而面对不同深度异常源时，即使异常幅值相同的异常源，由于深度的不同，导致各异常源权重大小相差巨大。深度加权矩阵在光滑反演、Marquardt反演和Occam反演中的效果均很差；而聚焦反演[19]和混合正则化反演在两类情况下，均表现出较好的聚焦效应(边缘保持效果)。相对于聚焦反演，混合正则化反演并没有出现过度聚焦效应，而表现出一类较弱的聚焦效应，并给出了准确的异常深度。这表明了本研究提出的算法的正确性、即对不同深度异常的适应性。
1 正演
根据万有引力定律，一个体积为V，剩余密度为ρ的地质体，在其外部空间任意一点的引力位为
(1)
式中：
、
分别为依赖于测点r和地质体r′的物理量和核函数；
为密度函数；引力常数
=6.672×10-8 m3/(g-1·s-2)。
在物性反演中，使用长方体对V离散，以式(1)计算各长方体的异常响应。任意长方体在坐标(x，y，z)处的异常值的计算公式可重写为[20]
(2)
式中：
和
分别为长方体体对角线上的两节点坐标；
, 
, 
；
；
。
2 重力场混合正则化反演
重力正演可以表示如下：
d=Gm (3)
式中：d为观测场值所组成的向量且其长度为Nd；G为相对应位场的正演算子，为
大小的矩阵；m为Nm长度的物性参数向量。
传统位场反演等价于求解线性方程(3)的反问题，在正则化理论基础上构建三维反问题
(4)
此处，J(m)通常为L1或L2罚项，
为数据泛函。为进一步提高物性反演成像在深度上的分辨率，引入L1多尺度小波与全变差正则化，在保留位场多尺度特征的同时，增强边界区分的能力，以更好地得到清晰构造边界的同时快速而稳定的求解。联合稀疏约束反演的可通过修改式(4)中J(m)实现，
(5)
式中：W为多尺度正交小波算子，
，I为单位阵；
为全变差算子。W和很相似，可确保找到多个正则化参数使得反演结果稳定。
式(4)可改写为
(6)
式中：γ和λ均为正则化参数。
将
写为泛函形式：
(7)
使用L1范数，则有[21]
(8)
其导数为
(9)
通常采用迭代最小二乘法序贯求解式(6)[22]，即将式(6)中数据泛函分别与两正则项联立，通过如下两步骤迭代求解m，
第一步：
(10)
第二步：
(11)
由于需要通过两步迭代求解式(6)，在解混合正则化问题上，采用传统的序贯方法远比常规单一正则化效率低。WRIGHT[23]提出固定点迭代算法，但其仅为线性收敛，存在收敛速度慢的问题；GOLDSTEIN等[24]提出了分裂迭代正则化算法，以降低模型的计算复杂性；JIA等[25]进一步给出了其收敛性分析。
为此，对式(6)引入收敛更快的最小分裂Bregman迭代算法[16]。利用p1和p2分别替换式(6)中的
和Wm两项，得一约束优化问题
(12)
将式(12)转化为无约束优化问题
(13)
对上式中m求导，并令求导式等于零，那么迭代形式为
(14)
式中：
(15)
，
、
和
为沿坐标系三轴向的差分运算；wk、
、
和
为反演迭代的临时中间值；
为数据加权矩阵。
利用Soft-shrinkage公式[26]对式(12)中约束项求解以更新式(14)中的wk、、和：
(16)
式中：
为软阀值收敛函数[27]，
(17)
和
(18)
混合正则化反演其主要受式(14)的求解速度控制，尤其在做海量面数据的三维反演时，这种瓶颈效应更加明显[27]，那么要减少计算时间，应保持反演收敛速度更快、迭代次数尽可能小和减少反演核矩阵参与运算的次数。式(14)的收敛速度和迭代次数受到系数矩阵的条件数控制，对于三维密度成像反演而言，系数矩阵的条件数相当大。为避免迭代次数太大，可以通过左乘一预条件矩阵，
(19)
对于小模型物性反演成像问题，直接对式(14)左乘Z1其计算效率将大幅提高；但对海量面数据的三维物性反演成像，引入预条件共轭梯度法可快速近似的计算式(14)[28]，当式(5)由多尺度小波算子W和全变差算子组合而成时，GTG的预条件算子可近似地为其自身的对角矩阵
(20)
3 深度加权
位场物性反演的一个目标是取得地下物性的分布，其深度分辨率远逊于重磁位场在确定异常边缘上的分辨率。主要原因在于重力值与场源到观测点的距离的幂次方成反比，从而使重力异常的幅值随着物性网格质心深度的增加而迅速衰减，由此产生一种现象—“趋肤效应”，一般可以通过引入深度加权函数近似补偿核函数随着深度的衰减。仅考虑深度方向上的距离，而忽略观测点与密度体在水平方向上的差异，可以使用
近似表示核函数随着深度的衰减。LI等[29]在三维重力反演中，引入深度加权函数
(21)
式中：z为物性网格中心埋深，z0和κ为常数，z0取决于物性网格的尺寸以及观测点高度。此时式(8)预条件算子可改写为
。通过调整z0和κ的值，以近似表达核函数的趋肤效应[29]。
相对于Z1与Z2，Z3求逆或左乘的计算成本远比前两者要小。在重力反演中κ一般在1.5~2.0之间[29]；而在重力梯度反演中，κ一般在2.0~3.0之间；在磁化率反演中，κ=3[30]。这致使在重磁位场多分量联合反演中难于确定κ。相比前者而言，PORTNIAGUINE等[31]给出了一个更为简单的基于反演核矩阵的先验深度加权Z：
(22)
式中：Gij对应于核矩阵G中第i行第j列的元素。
因此，在不考虑物性网格在水平方向上的差异及网格边界效应的影响时，Z将随着物性网格的深度增加而减小。若应用于式(4)上，Z将优先约束在浅层网格上，从而使得在作用深部的物性网格权重更重[32]。
引入深度加权项(11)，对式(4)改写
(23)
类似于式(6)，引入混合正则化项，则
(24)
式中：
，
；
、
和
的确定可以参考[33]。
通过类似于(12)~(14)推导，式(24)有如下的迭代形式
(25)
式中：
(26)
，
，
，
。
在地球物理重磁场反演中，使用非线性共轭梯度法或预条件共轭梯度法求解反问题，通常难于获得尖锐边界[9]。为在反演结果中保存物性边界，可以利用重加权正则化共轭梯度法求解式(25)[34-35]，并利用式(16)求解更新wk、、和以实现重力场混合正则项联合约束反演。
图1 重力正演模型及不同反演方法效果对比
Fig. 1 Forward model of gravity data and comparison of effect of different inversion methods
4 密度成像中的约束
对于密度成像而言，密度是唯一的反演参数。对于式(4)中的目标函数，可以施加如下密度约束条件
式中：mmin和mmax分别为异常上下限，但这种约束长导致异常倒置现象(如图1(f))。ZHDANOV等[36]给出一约束反演异常的对数变换对，正变换为
(27)
与其对应的反变换为
(28)
重力场混合正则项联合反演是在对数空间m′中进行，无论m′多么大，或多么小，反对数变换式(28)都能确保m在mmin与mmax之间。
而此时，式(24)则修改为
5 模型实验
5.1 实验模拟一
构建一正演模型(如图1(a))对上述方法原理进行验证，将场源空间划分为20×20×16=6400个物性网格，每个物性网格沿x、y和z轴方向的长度均为50 m。如图1(a)所示，地下空间存在两个异常体，两个异常体的大小均为200 m×200 m×200 m，异常体顶板埋深均为200 m，底板为400 m，异常体的剩余密度为1×103 kg/m3。总共有20×20=400个采集数据，测点高度为25 m，测点间距为50 m。
反演数据为添加3%的高斯白噪的正演模拟数据。在光滑反演(图1(b))、Marquardt反演(图1(c))及Occam反演(图1(d))中，λ=4和β=1.8；聚焦反演使用
作为聚焦算子，e=1.2×10-6(图1(e))；在混合正则化反演(图1(f))中，γ=0.01、λ=0.1和μ=1。在图1(b)、(c)和(d)中三者效果相当，三者在深度上的分辨率一般，尤以光滑反演(图1(b))在深度上的分辨率最差，异常底板分界面非常模糊，难以辨识；而聚焦反演取得了较为清晰的边界，由于聚焦作用的影响，存在过度聚焦的现象(图1(e))，即反演得出异常的规模远比实际的要小，而异常幅值要比模型实际要大得多。相比于聚焦反演，本研究反演方法的聚焦作用更为明显，但其反演轮廓与正演模型相当。
5.2 实验模拟二
鉴于光滑反演(图1(b))、Marquardt反演(图1(c))及Occam反演(图1(d))，在深度方向的分辨率较差，限于文章篇幅，便不再对比以上3种算法。为对上述方法原理进行验证，将场源空间划分为40×40×20=32000个物性网格，每个单元格沿x、y和z轴方向的长度均为50 m。总共有40×40=1600个采集数据，测点高度为25 m，测点间距为50 m。如图2所示，地下空间存在两个异常体，其大小均为200 m×200 m×200 m。1号和2号异常体的剩余密度分别为3×103 kg/m3，1×103 kg/m3。以100 m为步长ZΔ沿z轴向下逐一移动1号异常源，ZΔ分别为0、100、200、300和400。
反演数据为添加3%的高斯白噪的正演模拟数据。
图2 正演模型
Fig. 2 Forward model of gravity
在混合正则化反演(图3)中，γ=0.01、λ=0.1和μ=1；聚焦反演使用作为聚焦算子，e=1.2×10-6。异常上下限均为mmin=0和mmax=3。如图3(a)和(b)所示，当ZΔ=0时，聚焦反演与混合正则化反演都取得了较为清晰的边界，并将两者区分开来；混合正则化反演取得了与正演模型较为一致的结果(图3(b))；聚焦反演则在幅值上表现出一定的过度聚焦效应，1号和2号异常其空间范围与实际存在较大偏差，过度聚焦效应明显。当ZΔ=100时，混合正则化反演的2号异常幅值仍然保持在1×103 kg/m3，1号异常幅值保持在2×103 kg/m3，整体上其异常分布与正演模型较为一致(见图3(d))；聚焦反演其2号异常与正演模型较一致，而1号异常由于受过度正则化效应的影响，导致深度不同的异常源出现了一定的歪斜(见图3(c)及(i))。当ZΔ=200时，由于2号异常密度较大且埋深较深，导致其在总场效应上与1号异常相当(见图3(e)和(f))，混合正则化反演和聚焦反演的结果相当，这同时能与图1(e)和(f)相互验证。当ZΔ=300及ZΔ=400时，相较聚焦反演而言，混合正则化反演具有更强的适应性(见图3(g)和(h))，并有效地避免了过度正则化效应(见图3(i)和(j))。
图3 不同反演方法对不同埋深的异常源的反演效果对比
Fig. 3 Comparison of inversion effect of two different inversion methods on density abnormalities with different depth sources
6 结论
1) 充分分析在混合正则化框架下预条件算子的实现及计算性能，给出适合多阶重力分量的预条件算子；并在此基础上给出了适合多阶重力分量的深度加权公式。
2) 引入密度成像中的约束，确保反演异常幅值具有物理意义，数值结果能有效的降低过度聚焦效应。
3) 充分与多类单一正则化反演方法对比，实验结果表明混合正则化反演能有效的勾勒异常边界。
4) 在反演不同埋深的异常源时，与聚焦反演相比，混合正则化反演对复杂地质构造具有更强的适应性，能有效地避免聚焦反演过度聚焦的现象，并较为准确的描述不同深度异常源。
REFERENCES
[1] 王化祥, 唐 磊, 闫 勇. 电容层析成像图像重建的总变差正则化算法[J]. 仪器仪表学报, 2007, 28(11): 2014-2018.
WANG Hua-xiang, TANG Lei, YAN Yong. Total variation regularization algorithm for electrical capacitance tomography [J]. Chinese Journal of Scientific Instrument, 2007, 28(11): 2014-2018.
[2] 黄 嵩. 电阻抗静态成像中正则化算法研究[D]. 重庆: 重庆大学电气工程学院, 2005.
HUANG Song. Regularization algorithm research of static electrical impedance tomography[D]. Chongqing: College of Electrical Engineering Chongqing University, 2005.
[3] 叶益信, 胡祥云, 金钢燮, 陈玉萍, 何梅兴, 周晓晨. 大地电磁二维陡边界反演应用效果分析[J]. 地球物理学进展, 2009, 24(1): 668-674.
YE Yi-xin, HU Xiang-yun, KIM Kang-sop, CHEN Yu-ping, HE Mei-xing, ZHOU Xiao-chen. Application analysis of sharp boundary inversion of magnetotelluric data for 2D structure[J]. Progress in Geophysics, 2009, 24(1): 668-674.
[4] 张宏兵, 尚作萍, 杨长春, 段秋梁. 波阻抗反演正则参数估计[J]. 地球物理学报, 2005, 48(1): 181-188.
ZHANG Hong-bing, SHANG Zuo-ping, YANG Chang-chun, DUAN Qiu-liang. Estimation of regular parameters for the impedance inversion[J]. Chinese Journal Geophysics, 2005, 48(1): 181-188.
[5] 李星秀, 韦志辉. 形态成分正则化约束的图像恢复方法[J]. 计算机工程与应用, 2010, 46(17): 27-29.
LI Xing-xiu, WEI Zhi-hui. Image restoration via regularization constraints of morphological components[J]. Computer Engineering and Applications, 2010, 46(17): 27-29.
[6] WAN L, ZHDANOV M. Focusing inversion of marine full-tensor gradiometry data in offshore geophysical exploration[J]. SEG Technical Program Expanded Abstracts, 2008, 27(1): 751-755.
[7] PORTNIAGUINE O, ZHDANOV M. Focusing Geophysical Inversion Images[J]. Geophysics, 1999, 64(3): 874-887.
[8] 张罗磊, 于 鹏, 王家林, 吴健生, 陈 晓, 李 洋. 二维大地电磁尖锐边界反演研究[J]. 地球物理学报, 2010, 53(3): 631-637.
ZHANG Luo-lei, YU Peng, WANG Jia-lin, WU Jian-sheng, CHEN Xiao, LI-yang. A study on 2D magnetotelluric sharp boundary inversion[J]. Chinese Journal Geophysics, 2010, 53(3): 631-637.
[9] 张罗磊, 于 鹏, 王家林, 吴健生. 光滑模型与尖锐边界结合的MT二维反演方法[J]. 地球物理学报, 2009, 52(6): 1625-1632.
ZHANG Luo-lei, YU Peng, WANG Jia-lin, WU Jian-sheng. Smoothest model and sharp boundary based two-dimensional magnetotelluric inversion[J]. Chinese Journal Geophysics, 2009, 52(6): 1625-1632.
[10] STRONG D, CHAN T. Edge-preserving and scale-dependent properties of total variation regularization[J]. Inverse Problems, 2003, 19(6): S165-S187.
[11] 岳建惠.电阻率成像反问题的混合正则化方法研究[D]. 大连: 大连海事大学, 2012.
YUE Jian-hui. Research of mixture regularization methods for EIT inverse problem[D]. Dalian: Dalian Maritime University, 2012.
[12] MARTINS M, LIMA A, BARBOSA C, SILVA B. Total variation regularization for depth-to-basement estimate: Part 1—mathematical details and applications[J]. Geophysics, 2011, 76(1): I1-I12.
[13] LIMA W A, MARTINS C M, SILVA J B C, BARBOSA V C F. Total variation regularization for depth-to-basement estimate: Part 2—Physicogeologic meaning and comparisons with previous inversion methods[J]. Geophysics, 2011, 76(1): 13-20.
[14] LI Y, OLDENBURG W. 3-D inversion of induced polarization data[J]. Geophysics, 2000, 65(6): 1931-1945.
[15] KANE J, HERRMANN F, TOKSOZ N. Wavelet domain linear inversion with application to well logging[R]. Massachusetts: Massachusetts Institute of Technology Earth Resources Laboratory, 2001.
[16] GHOLAMI A, SIAHKOOHI H. Regularization of linear and non-linear geophysical ill-posed problems with joint sparsity constraints[J]. Geophysical Journal International, 2010, 180(2): 871-882.
[17] THANASSOULAS C, TSELENTIS A, DIMITRIADIS K. Gravity inversion of a fault by Marquardt’s method[J]. Computers & Geosciences, 1987, 13(4): 399-404.
[18] CONSTABLE C, PARKER L, CONSTABLE G. Occam’s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data[J]. Geophysics, 1987, 52(3): 289-300.
[19] ZHDANOV M, ELLIS R, MUKHERJEE S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004, 69(4): 925-937.
[20] LI X, CHOUTEAU M. Three-dimensional gravity modeling in all space[J]. Surveys in Geophysics, 1998, 19(4): 339-368.
[21] BREDIES K, KUNISCH K, POCK T. Total generalized variation[J]. SIAM Journal on Imaging Sciences, 2010, 3(3): 492-526.
[22] ABBASBANDY S. Improving Newton-Raphson method for nonlinear equations by modified Adomian decomposition method[J]. Applied Mathematics and Computation, 2003, 145(2): 887-893.
[23] WRIGHT S J. Primal-dual interior-point methods[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1997.
[24] GOLDSTEIN T, OSHER S. The split Bregman method for L1-regularized problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(2): 323-343.
[25] JIA Q, ZHAO H, ZHAO W. Convergence analysis of the Bregman method for the variational model of image denoising[J]. Applied and Computational Harmonic Analysis, 2009, 27(3): 367-379.
[26] WANG Y, YIN W, ZHANG Y. A fast algorithm for image deblurring with total variation regularization[R]. Texas: Rice University, 2007.
[27] DONOHO D L. Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition[J]. Applied and Computational Harmonic Analysis, 1995, 2(2): 101-126.
[28] KOH K, KIM J, BOYD P. An interior-point method for large-scale L1-regularized logistic regression[J]. Journal of Machine Learning Research, 2007, 8(8): 1519-1555.
[29] LI Y, OLDENBURG W. 3-D inversion of gravity data[J]. Geophysics, 1998, 63(1): 109-119.
[30] VASILEVSKY A, DRUZHININ A, EVANS R, MURPHY C. Feasibility of FTG reservoir monitoring[J]. SEG Technical Program Expanded Abstracts, 2003, 22(1): 1450-1453.
[31] PORTNIAGUINE O, ZHDANOV M. 3-D magnetic inversion with data compression and image focusing[J]. Geophysics, 2002, 67(5): 1532-1541.
[32] PILKINGTON M. 3D magnetic data-space inversion with sparseness constraints[J]. Geophysics, 2008, 74(1): L7-L15.
[33] HE B, YANG H, WANG S. Alternating direction method with self-adaptive penalty parameters for monotone variational inequalities[J]. Journal of Optimization Theory and applications, 2000, 106(2): 337-356.
[34] ZHDANOV M. Geophysical inverse theory and regularization problems[M]. Boston: Elsevier Science Ltd, 2002.
[35] FLETCHER R. Practical methods of optimization[M]. New York: Wiley IEEE Press, 2013.
[36] ZHDANOV M, CAI H, Wilson G. 3-D inversion of squid magnetic tensor data[J]. Journal of Geology And Geosciences, 2012, 1(104): 1-5.
(编辑 李艳红)
基金项目：国家自然科学基金资助项目(41174061，41374120)
收稿日期：2013-08-08；修订日期：2013-12-06
通信作者：朱自强，教授，博士，电话：13507319431；E-mail：13507319431@139.com
