基于分裂Bregman迭代的混合正则化重力场反演

曹书锦，朱自强，鲁光银



基于分裂Bregman迭代的混合正则化重力场反演

曹书锦，朱自强，鲁光银

(中南大学地球科学与信息物理学院，湖南长沙，410083)

针对重力场反演中深度分辨率过低的问题，提出基于混合稀疏正则化联合反演。引入L1多尺度小波和全变差技术，使在剔除干扰异常及规避反演过度聚焦的同时，尽可能地保留深度方向上的分辨率；针对在混合正则化反演中预条件矩阵左乘(求逆)计算量过大的问题，构建新的深度权项实现基于最小分裂Bregman迭代算法快速三维密度反演成像。研究结果表明：本文反演方法可行、有效，具有较强的适应性。

混合正则化；分裂Bregman迭代算法；深度加权矩阵；多尺度小波；全变差正则化

由于反演着重点或目的不同，所需要建立的反演模型也不同，同时也决定了反演的总体设计、反演过程、反演特征及反演效果等不同。本质上，每种正则化方法都有一定的侧重性，而在面对较复杂的地质构造时，这种优势即转变成过度正则化效应。同时，由于重力异常的幅值随着地下物性网格质心深度的增加而迅速衰减，导致反演深度分辨率过低[1]，如何得到更清晰的地质界面仍是当前地球物理反演研究的一个重点问题。在过去的十几年中，众多研究者在边界识别/边界保存的位场反演技术如全变差正则化[2]、聚焦反演[3]和尖锐边界/陡边界反演[4]等方面进行了大量的研究。Tikhonov正则化虽然具有成像速度快、求解稳定、定位准确等优点，但由于该算法导出的是最小二乘问题近似解，导致目标区域与背景区域间的边界模糊，成像效果不好，分辨能力较低，多个物体之间存在伪迹。全变差正则化算法能在对问题进行反演过程中得到跳跃性较大的参数部分，有效地避免了物性边界的过度光滑，使得目标区域与背景区域间的边界分明，具有很好的保边缘性[5]。仅在先验地质信息充足的情况下，尖边界反演方法才能取得很好的反演效果，但其对初始模型的设置和先验地质信息的准确性要求较高，难于满足起伏地形条件下的反演需求，且存在反演深部构造不够准确以及加入约束条件较困难等不足[6]。聚焦反演具有成像快速、对单一异常聚焦效果明显等优点，但由于存在过度正则化效应，在地质构造复杂且各构造权重不均一的情况下，其易导致构造形态发生畸变及对目标区域的定位不准确[7]。为此，本文作者引入一类多尺度小波稀疏算子，以期在边界保持的同时借助其具有的一定光滑效应，在此基础上引入具有边界保持特性的全变差算子，将两者在分裂Bregman迭代算法的框架下有机地实现混合正则化反演[8]。针对在分裂Bregman迭代算法框架下的预条件算子构建难度大(即当预条件矩阵为左除矩阵时，其为稀疏矩阵，由于每次反演迭代都需要计算其临时中间逆矩阵，但其逆矩阵为非稀疏矩阵，因而导致不能有效地存储预条件矩阵的逆矩阵)，引入一类与场源衰减特性无关的预条件矩阵，以实现快速反演。使用2组不同模型(同一深度与不同深度)和实际观测数据，对比研究多类位场反演方法在不同深度异常的适应性，并通过数值模型和实际观测数据的反演验证本文算法的正确性，即对不同深度异常的适应性。

1 混合正则化反演

将反演的地下空间划分为大小相等、紧密排列的立方体格网单元，且任一立方体内部的密度均匀[9]。重力正演可以表示如下：

其中：表示观测重力，可具体为g；为示模型空间密度分布的重力正演核函数；物性参数向量。重力场物性反演待求的参数就是各网格单元内的密度。传统位场反演等价于求解线性方程(1)的反问题，利用正则化理论构建三维反问题[10]：

其中：()通常为L1或L2罚项，为模型泛函；为数据泛函，为正则化参数。

单一正则化方法都有一定的侧重性，也意味着对相应的问题具有非常重要的单一优势，而在面对较复杂的地质构造时，这种优势则转变成过度正则化效应。为确保反演的稳定性和健全性，又不损失重构结果的分辨率，引入L1多尺度小波[11−12]与全变差正则化算子[13]以构建重力场联合稀疏约束反演目标函数，以期在保留位场多尺度特征的同时，增强边界区分的能力，从而在更好地得到清晰构造边界的同时快速而稳定地求解。联合稀疏约束反演的可通过修改式(2)实现[8]：

使用L1范数，则有[15]

其导数为

其中，导数定义为：

；

；

；

。

其中：h，h和h分别为物性网格在，和方向的边长。函数(,)的定义为

其中：sign为符号函数，

当与的方向相反时，(,)则等于0。对于在物性不连续接触界面(如断层接触面、地层分界面等)，

通常采用重加权迭代最小二乘法序贯求解式(3)[16]，但采用传统的序贯方法解混合正则化优化问题的效率较低。在解混合正则化问题上，Wright[17]提出了固定点迭代算法，但其仅为线性收敛，存在收敛速度慢的问题；Goldstein等[10]提出了分裂迭代正则化算法，以降低模型的计算复杂性；Jia等[18]在其基础上，对其收敛性进行了进一步分析。为此，对式(3)引入收敛更快的最小分裂Bregman迭代算法框架以实现重力场联合稀疏约束反演。利用1和2分别替换式(3)中的和等2项，得一约束优化问题[10]：

将式(5)转化为无约束优化问题：

对式中求导，并令求导式等于0，则迭代形式为

其中：

2 深度加权

用式(6)进行物性反演，其目标是取得地下物性的分布。由于重力没有或较少含有深度方向的信息，且在深度方向的分辨率较低，其确定异常边缘上的能力较强，主要源于重力异常的幅值与场源到观测点的距离的幂次方成反比，而使重力异常的幅值随着物性网格质心深度的增加而迅速衰减产生的趋肤效应，导致重构结果的密度分布趋向地表附近，而不是按照地质体的真实深度分布，因此，可以通过引入深度加权函数近似地补偿核系数矩阵随着深度的衰减。仅考虑深度方向的上距离而忽略观测点与密度体水平方向上的差异，可以使用近似反映核函数随着深度的衰减。Li等[19]在三维磁化率反演中，引入深度加权函数：

式中：为物性网格质心埋深；0为常数，取决于物性网格的边长以及观察面高度；为衰减因子，取决于地球物理位场的场源特性。通过调整0和，以近似表达反演核函数的趋肤效应。

混合正则化反演其主要受式(6)的求解速度控制，尤其是进行海量面数据的三维反演时，这种瓶颈效应更加明显。为了加快反演速度，应该使反演收敛速度更快、迭代次数尽可能小并减少反演核矩阵参与运算的次数。式(6)的收敛受系数矩阵的条件数控制，对于三维密度成像反演而言，其系数矩阵的条件数很大，一般可以通过左乘1个预条件矩阵优化目标函数。Pilkington[20]引入深度加权矩阵式(7)到三维磁化率反演目标函数中，得出简单可靠的预条件矩阵(其中=3.0，为单位矩阵)，并得到了比较理想的结果。但在重力反演中，一般为1.5~2.0[21]；在重力梯度反演中，一般为2.0~3.0；在磁化率反演中，=3.0[20]。这致使在重力场多分量联合反演中难以确定。对于混合正则化密度反演，其预条件矩阵为

对于小模型密度反演成像问题，直接对式(6)左乘预条件矩阵2其计算效率将很快。但对大尺度模型的三维物性反演，引入预条件共轭梯度法可快速近似地计算式(6)[22]；当混合正则化反演中两正则项分别为多尺度小波算子和全变差算子时，的预条件算子可近似地为其自身的对角矩阵：

2与3均为稀疏带状矩阵，因而，其左乘(求逆)的计算和存储成本远远比1的高，且对于相对较大模型的三维密度反演成像时，3的逆矩阵将极难计算或存储。但由于1应用于多阶重力场联合反演并不准确，通过以上3类预条件算子1，2与3相比，选择一个更简单的基于反演核矩阵的先验深度权[23]：

因此，在不考虑物性网格在水平方向上的差异及网格边界效应的影响时，将随着物性网格深度的增大而减小。若应用于(2)式，则将优先约束在浅层网格上，从而使得在作用深部的物性网格权重更大[24]，此时，便可以作为求解式(6)的预条件矩阵。

引入深度加权项式(8)，对式(2)改写为

类似于式(3)，引入2个稀疏正则化项，并将数据泛函的正则化参数缩并到两混合正则化参数内：

深度加权的分裂Bregman迭代无约束优化算法的具体实施步骤如下。

1) 选定误差限t、正则化参数和。

4) 鉴于在地球物理重磁场反演时，通常使用非线性共轭梯度法或预条件共轭梯度法易于获得光滑边界，而难于获得尖锐边界。因此，利用重加权正则化共轭梯度法求解式(10)以实现重力场联合稀疏约束反演。多正则化参数的重加权正则化共轭梯度法算法流程见文献[3, 25]。

5) 利用soft−shrinkage公式[26]对式(5)中约束项更新：

7) 重复步骤3)~6)，直至满足反演终止条件为止。

3 模型实验

3.1 实验模拟Ⅰ

如图1(a)所示，构建1个含有2个同深度异常体的模型对上述方法原理进行验证。2个异常体的长×宽×高均为200 m×200 m×200 m，异常体顶板埋深均为200 m，底板埋深为为400 m，剩余密度均为1.0 g/cm3。将场源空间划分为20×20×10=4 000个单元格，每个单元格沿，和轴方向的长度均为50 m；观测网测点高度为25 m，测点间距为50 m，共有20×20=400个采集数据。

(a) 正演模型；(b) 光滑反演结果；(c) Marquardt 反演结果[27]；(d) Occam反演结果[28]；(e) 聚焦反演结果[29]；(f) 混合正则化反演结果

反演数据为添加3%的高斯白噪声的正演模拟结果。在光滑反演、Marquardt反演及Occam反演中，=4，a=1，a=1，a=1，a=0.000 5，=1.8，如图1(a)~(d)所示[27−29]；在图1(e)中，聚焦算子为，=1.2×10−6；图1(f)中，混合正则化反演的参数为=0.01，=0.1，=1。在图1(b)，(c)和(d)中三者效果相当，三者在深度上的分辨率均较低，尤以光滑反演(图1(b))在深度上的分辨率最低，异常底板分界面非常模糊，难以辨识；而聚焦反演取得了较清晰的边界，由于聚焦作用的影响，存在过度聚焦的现象(图1(e))，即反演得出异常的规模远比实际的要小，而异常幅值要比实际正演模型的大得多。相比于聚焦反演(如图1(f)所示)，本文混合正则化反演方法的聚焦作用更明显，但其反演轮廓与正演模型的相当。

3.2 实验模拟Ⅱ

构建1个含2个不同深度异常体的模型对上述方法原理进行验证，如图2(a)所示。2个异常体长×宽×高均为200 m×200 m×200 m；右侧异常体顶板埋深为200 m，底板为400 m，左侧异常体顶板埋深均为300 m，底板为500 m；剩余密度均为1.0 g/cm3。将场源空间划分为20×20×16=6 400个单元格，每个单元格沿，和轴方向的长度均为50 m；观测网测点高度为25 m，观测点间距为50 m，总共有20×20=400个采集数据。

(a)正演模型；(b) 光滑反演结果；(c) Marquardt 反演结果[27]；(d) Occam反演结果[28]；(e) 聚焦反演结果[29]；(f) 混合正则化反演结果

反演数据为添加3%的高斯白噪声正演模拟结果，在光滑反演、Marquardt反演及Occam反演中，=4，a=1，a=1，a=1，a=0.000 5，=1.8，如图2(a)~(d)所示；在混合正则化反演中，=0.01，=0.1，=1，如图2(f)所示。在图2(b)，(c)和(d)中三者效果相当，异常底板分界面非常模糊，难以辨识；而聚焦反演取得了较为清晰的边界(图2(e))，将两者区分开来，但其反演推断异常空间位置与正演模型的存在较大偏差，尤其在深度方向存在差异较大。这主要受到聚焦算子的适应性和聚焦阀值的难以选取的影响，此处聚焦算子为，=0.001，其中与前次模型(即图1(e))相差很大，其原因为：在限制反演密度模型上限的同时(异常幅值上限为1.2 g/cm3)，进一步提高了深度方向上的聚焦效应。相比于聚焦反演，本文反演方法的聚焦作用更明显，其聚焦作用介于Occam反演、Smoothness反演和Marquardt反演这三者与聚焦反演之间，其反演结果轮廓与正演模型的相当，且深度与实际正演模型的相当(图2(f))。与其他反演方法相比，混合正则化反演具有更强的适应性。

4 实例

加拿大卑诗省在2007年设施1个旨在为吸引矿业投资者向地质构造复杂、勘探程度较低、潜在成矿潜力巨大的British Columbia地区投资的地球物理及地球化学综合勘探项目[30]。本文实例采用该项目组中的1个名为Quest—South的测区，见图3。为便于处理，选取其中1个矩形测区(见图4)。在测区内地势从丘陵到连绵陡峭的山脉，海拔高程为−150~370 m。

图3 Quest South测区示意图[30]

图4 重力异常

为对重力异常进行分离，用向上延拓到3 km平面的重力场作为区域场，将区域异常从重力异常中减去从而得到剩余异常场，见图5。

图5 局部重力异常

从图5可以发现反演算法在浅部的正异常取得了较清晰的结果。反演密度切片见图6。图6中的(1)，(3)和(4)号异常均具有较清晰的边界；而图6中的(2)和(5)号异常所示在深部反演的分辨率有所下降，这主要受到区域场和剩余异常场划分的影响。

图6 混合正则化重力场反演结果组合切片图

5 结论

1) 分析了混合正则化框架下预条件算子的实现及计算性能，给出了适合多阶重力场分量的深度权公式，并在此基础上给出了混合正则化反演算法。

2) 充分对比了多类单一正则化反演方法，引入基于L1多尺度小波正则项和全变差正则的混合正则化，在剔除无效异常的同时，尽可能地保留深度方向上的分辨率及在边界上的识别能力。

3) 与其他反演方法相比，混合正则化反演对不同深度异常源具有更强的适应性，并有效地避免了聚焦反演存在的过度聚焦现象。

[1] 姚长利, 郝天珧, 管志宁, 等. 重磁遗传算法三维反演中高速计算及有效存储方法技术[J]. 地球物理学报, 2003, 46(2): 2252−2258. YAO Changli, HAO Tianyan, GUAN Zhining, et al. High speed computation and efficient storage in 3D gravity and magnetic inversion based on genetic algorithms[J]. Chinese Journal of Geophysics, 2003, 46(2): 2252−2258.

[2] 李星秀, 韦志辉. 形态成分正则化约束的图像恢复方法[J]. 计算机工程与应用, 2010,46(17): 27−29. LI Xingxiu, WEI Zhihui. Image restoration via regularization constraints of morphological components[J]. Computer Engineering and Applications, 2010, 46(17): 27−29.

[3] Zhdanov S. Geophysical inverse theory and regularization problems[M]. Boston: Elsevier Science Ltd, 2002: 155−165.

[4] 张罗磊, 于鹏, 王家林, 等. 光滑模型与尖锐边界结合的MT二维反演方法[J]. 地球物理学报, 2009, 52(6): 1625−1632. ZHANG Luolei, YU Peng, WANG Jialin, et al. Smoothest model and sharp boundary based two-dimensional magnetotelluric inversion[J]. Chinese Journal of Geophysics, 2009, 52(6): 1625−1632.

[5] 岳建惠. 电阻率成像反问题的混合正则化方法研究[D]. 大连: 大连海事大学应用数学系, 2012: 4−5. YUE Jianhui. Research of mixture regularization methods for EIT inverse problem[D]. Dalian: Dalian Maritime University. Department of Applied Mathematics, 2012: 4−5.

[6] 谭捍东, 余钦范, John B, 等. 大地电磁法三维快速松弛反演[J]. 地球物理学报, 2003, 46(6): 1218−1226. TAN Handong, YU Qinfan, John B, et al. Three-dimensional rapid relaxation inversion for the magnetotelluric method[J]. Chinese Journal of Geophysics, 2003, 46(6): 1218−1226.

[7] 黄嵩. 电阻抗静态成像中正则化算法研究[D]. 重庆: 重庆大学电气工程学院, 2005: 27−50. HUANG Song. Regularization algorithm research of static electrical impedance tomography[D]. Chongqing: Chongqing University. College of Electrical Engineering, 2005: 27−50.

[8] 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.

[9] LI Xiong, Chouteau M. Three-dimensional gravity modeling in all space[J]. Surveys in Geophysics, 1998, 19(4): 339−368.

[10] Goldstein T, Osher S. The split Bregman method for L1- regularized problems[J]. SIAM Journal on Imaging Sciences, 2009,2(2): 323−343.

[11] Davis K, Li Y. Efficient 3D inversion of magnetic data via octree mesh discretization, space-filling curves, and wavelets[J]. SEG Technical Program Expanded Abstracts, 2010, 29(1): 1172−1177.

[12] Davis K, LI Yaoguo. Fast solution of geophysical inversion using adaptive mesh, space-filling curves and wavelet compression[J]. Geophysical Journal International, 2011, 185: 157−166.

[13] Routh S, Qu L, Sen K, et al. Inversion for non-smooth models with physical bounds[J]. SEG Technical Program Expanded Abstracts, 2007, 26(1): 1795−1799.

[14] Bertete-Aguirre H, Cherkaev E, Oristaglio M. Non-smooth gravity problem with total variation penalization functional[J]. Geophysical Journal International, 2002,149(2): 499−507.

[15] Loris I, Verhoeven C. Iterative algorithms for total variation-like reconstructions in seismic tomography[J]. International Journal on Geomathematics, 2012, 3(2): 179−208.

[16] Abbasbandy S. Improving Newton-Raphson method for nonlinear equations by modified a domian decomposition method[J]. Applied Mathematics and Computation, 2003, 145(2): 887−893.

[17] Wright J. Primal-dual interior-point methods[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1997: 127−156.

[18] Jia R 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.

[19] LI Yaoguo, Oldenburg W. 3-D inversion of magnetic data[J]. Geophysics, 1996, 61(2): 394−408.

[20] Pilkington M. 3-D magnetic imaging using conjugate gradients[J]. Geophysics, 1997,62(4): 1132−1142.

[21] LI Yaoguo, Oldenburg W. 3-D inversion of gravity data[J]. Geophysics, 1998, 63(1): 109−119.

[22] Koh K, Kim S 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.

[23] Portniaguine O, Zhdanov S. 3-D magnetic inversion with data compression and image focusing[J]. Geophysics, 2002, 67(5): 1532−1541.

[24] Pilkington M. 3D magnetic data-space inversion with sparseness constraints[J]. Geophysics, 2008,74(1): L7−L15.

[25] ZHANG Luolei, YU Peng, WANG Jialin, et al. Smoothest model and sharp boundary based two-dimensional magnetotelluric inversion[J]. Chinese Journal of Geophysics, 2009, 52(6): 1360−1368.

[26] WANG Yilun, YIN Wotao, ZHANG Yin. A fast algorithm for image deblurring with total variation regularization[R]. Texas: Technical Report TR07-10 of Rice University Computational and Applied Mathematics, 2007: 1−19.

[27] Thanassoulas C, Tselentis G A, Dimitriadis K. Gravity inversion of a fault by marquardt’s method[J]. Computers & Geosciences, 1987,13(4): 399−404.

[28] 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.

[29] Zhdanov M, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004, 69(4): 925−937.

[30] Reichheld S. Documentation and assessment of exploration activities generated by geoscience BC data publications[R]. British Columbia, CANADA: Geoscience BC, 2013: 125−130.

Gravity inversion based on hyper-parameter regularization inversion via iteration splitting bregman algorithm

CAO Shujin, ZHU Ziqiang, LU Guangyin

(School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Considering the depth resolution of gravity inversion is too low, a hyper-parameter regularization algorithm based on multi-scale L1 wavelet operator and total variation operator was proposed. Regarding the difficulty to calculate the inverse matrix of preconditioning matrix, a new preconditioner was introduced to carry out the hyper-parameter regularization inversion, which was solved by the minimum splitting Bregman iterative algorithm. The results show that the proposed method is vilid and effective, and has strong adaptability.

hyper-parameter regularization; splitting Bregman iterative algorithm; depth weighting matrix; multi-scale wavelet; total variation operator

10.11817/j.issn.1672-7207.2015.05.018

P631

A

1672−7207(2015)05−1699−08

2014−05−28；

2014−07−21

国家自然科学基金资助项目(41174061, 41374120) (Projects(41174061, 41374120) supported by the National Natural Science Foundation of China)

曹书锦，博士，从事地球物理勘探数据处理与反演等研究；E-mail: shujin.cao@163.com

(编辑 陈灿华)