一维大地电磁Occam反演拉格朗日乘子的搜索

张君涛，周 军，王绪本，夏时斌，钟红梅

（1.成都理工大学地球探测与信息技术教育部重点实验室，成都 610059；2.四川省核工业地质调查院，成都 610061）

一维大地电磁Occam反演拉格朗日乘子的搜索

张君涛1，周 军1，王绪本1，夏时斌1，钟红梅2

（1.成都理工大学地球探测与信息技术教育部重点实验室，成都 610059；2.四川省核工业地质调查院，成都 610061）

在大地电磁反演中，Occam法因其在反演稳定性和模型分辨率等方面的优势，得到广泛应用。但由于其每次迭代都需要不断地搜索拉格朗日乘子，因而拉格朗日乘子的搜索效率对Occam法反演的运算速度起着至关重要的作用。为提高拉格朗日乘子的搜索效率，这里提出将拉格朗日乘子的搜索从自然数域中转到以10为底的对数域中进行，同时在区间最小值的搜索中采用二次函数极小值搜索法。通过大量一维大地电磁反演验证，该方法将拉格朗日乘子的搜索效率提高了20%～50%，大大提高了Occam反演的运算速度，方便了其在高维反演中的应用。

大地电磁；Occam反演；拉格朗日乘子；对数；二次函数搜索

0 引言

在大地电磁反演中，利用目标函数梯度信息的最优化方法占据着主要的地位［1－3，9］，包含Occam法［4－5］、非线性共轭梯度法（NLCG）［7］、拟牛顿法等反演方法。其中Occam法因其在反演稳定性和模型分辨率等方面的优势，在现今的大地电磁反演中得到了广泛地应用，但由于该方法在每次反演迭代中，都要搜索一个满足数据误差最小或者模型粗糙度最大的拉格朗日乘子，因而需要进行多次正演，增加了反演的计算量，限制了其在高维反演的应用。

很多学者在优化Occam反演的计算量上做了大量的研究，Siripunvaraporn［6－7］等研究了数据空间Occam反演，并且成功的应用到了MT二维、三维反演，极大的加强了Occam法在二维、三维反演中的应用；吴小平［11］等在讨论了Occam反演中数据拟合差随拉格朗日乘子变化的基础上，采用拉格朗日乘子在一定步长下逐次递减的求取方法，每次迭代只需一次正演，极大地提高了计算速度，但初始拉格朗日乘子及减小比例的设定对反演稳定性及结果影响较大。在对拉格朗日乘子搜索方式的改进上，也有很多人提出了改进方式［12－13］，但依然有很多地方值得探讨。

基于上述考虑，作者对Occam法反演中提高拉格朗日乘子的搜索效率上做了进一步研究，结合拉格朗日乘子在Occam方法反演迭代中的变化特性，提出了对该参数应当采用对数变换，并结合二次搜索方法来开展搜索计算。将上述方法应用到了一维Occam反演的拉格朗日乘子搜索中，试图提高拉格朗日乘子搜索的计算速度。这里首先介绍了Occam法反演的原理，然后详细介绍了对数二次法搜索原理，并进行一维模型算例分析，说明了其快速搜索的优势。

1 Occam反演基本原理

一维Occam反演的目标函数包含数据目标函数和模型粗糙度函数：

其中：μ为拉格朗日乘子；大括号内的为数据目标函数：d为实测数据；F为正演算子；W为数据误差加权矩阵；X＊为拟合差期望值；‖∂zm‖2为模型粗糙度函数。

反演的过程就是寻找使目标函数得到最小值的模型。为此Occam法将正演算子F线性化，在初始模型m0处按一阶泰勒级数展开：

其中：J（m0）为正演算子F对模型参数m0的导数，称为雅可比矩阵。将式（2）带入式（1）中，由此得到的目标函数已经简化为一个关于模型参数m的二次函数：

其中：d0＝d－F（m0）＋Jm0，J＝J（m0），再求式（3）关于m的一阶导数，得到梯度函数g（m）：

此时的梯度函数是一个关于模型参数的一次函数，令梯度g（m）＝0，求得的模型参数即为目标函数取得极小值的模型：

由于正演算子是非线性的，目标函数也不是二次函数，故反演过程中需要不断迭代，反复修正初始模型m0，重复式（5）过程，最终求得目标函数的最小值。

在每次迭代中，对于不同的拉格朗日乘子u，利用式（5）可得到不同的模型参数及其相应的绝对拟合误差。相比其他反演方法直接输入一个经验的拉格朗日乘子，Occam法反演提出在每次迭代中都要通过搜索，得到一个满足数据拟合误差最小的μ值（如果数据拟合误差低于目标拟合误差，则搜索满足目标数据拟合差的最小模型粗糙度的μ值）。用此μ值求得模型参数m再作为下一次迭代的初始模型，进行第二次迭代。这样反复迭代，最后找到一个满足目标拟合差且粗糙度最小的模型，或者是达到最大迭代次数后得到的拟合差最小的模型。正是由于Occam法反演对拉格朗日乘子的搜索，使其具有很高的反演稳定性和模型分辨率，但同时也增加了反演的运算量。

2 拉格朗日乘子的搜索方法

考虑到拉格朗日乘子的搜索效率对Occam法反演运算时间的巨大影响，作者对拉格朗日乘子的搜索方式进行改进。

所有患者采用随机数字表法分为舌下含服组和口服组,每组各40例。口服组患者口服酒石酸美托洛尔片(阿斯利康制药有限公司,产品批号1407184)25 mg(1片);舌下含服组患者舌下含服酒石酸美托洛尔片25 mg(1片)。

对于μ值的搜索，一般采用先搜索到一个谷值区间，然后在此区间上搜索数据拟合差的最小值，搜索到最小值后进行下一次迭代；而在搜索过程中，一旦有数据拟合误差小于目标拟合误差就直接进入满足拟合差的最小模型粗糙度的μ值搜索。这里也采用以上搜索策略，但做了两点改进。

2.1 拉格朗日乘子对数化搜索

理论上μ值可以为大于“0”的任意数，但如果直接对μ进行搜索，在区间搜索过程中，μ值的修正步长是随着μ值的变化而变化的，且初始μ值如果给的不合适，区间搜索的次数会相对较大。考虑到将原目标函数中的拉格朗日乘子μ变为10λ后，并不会改变其原有目标函数的极值特性，只要满足μ＝10λ其对应的模型参数及其对应的目标函数与数据拟合差都是一致的，因而作者大胆设想将μ值做对数变化，从而将目标函数中对μ值的搜索转为对10λ中的λ做搜索。将μ＝10λ代入到目标函数式（1）中得式（6）。

式（6）经过式（3）、式（4）、式（5）运算过程可得每次迭代中模型参数的求取公式：

可以看出式（7）只是将式（5）中的μ变为10λ，从而将对μ的搜索转为对λ的搜索，这将更加快速地搜索到谷值区间，并减小初始λ值对搜索效率的影响。

从图1可以看出，图1（a）绝对拟合误差曲线整体上变化平缓，更近似一个二次曲线；如果我们将区间搜索步长在以10为底的对数上设为“1”，一般在左右几个步长之内就能搜索到谷值区间，从而摆脱了初始值的设定对搜索次数的影响。图1（b）绝对拟合误差曲线在最小值左右变化相当剧烈，而在μ值较大时变化相当平缓。这给区间搜索步长的确定带来很大的不便，且初始μ值对区间搜索的次数有很大的影响。正是基于上述因素，这里更加明确将拉格朗日乘子的搜索从自然数域中转到了以10为底的对数域进行，将拉格朗日乘子变化为μ＝10λ，从而将对μ值的搜索改为对λ的搜索，其实质就是在以10为底的对数上进行μ值搜索。为便于进行谷值区间最小值搜索，将搜索到的谷值区间三点按λ值从大到小保存为（λ0，rms0）、（λ1，rms1）、（λ2，rms2）。

图1 模型一首次迭代绝对拟合误差曲线Fig.1 The absolute error curve of the first iteration of the model 1 （a）拉格朗日乘子以10为底；（b）拉格朗日乘子用自然数

2.2 二次函数搜索区间极小λ值

做了上述变化后能快速地搜索到一个谷值区间，但此区间的范围相对自然数搜索的区间较大，如果区间最小值搜索继续使用二分法进行，则搜索次数相对自然数搜索的谷值区间二分法，增加了很多次。同样以图2（a）为例，首次迭代对数搜索区间次数为4次，但用二分法进行区间最小值搜索则用了8次，共计12次；自然数区间搜索用了9次，但最小值搜索只用了5次，共14次；二者相比，对数搜索基本没有优势。因此，提出利用二次函数极小值来搜索区间最小值。

构建二次函数：

本文对于退出搜索需要满足以下要求：

其中ε≤1。当满足式（9）要求后，则取拟合差小的λ值及其相对应的模型进入下一次迭代。

3 算例分析

为了验证本方法的搜索效率，建立了四种地电模型（图2）。特别说明模型一为陈小斌［10］在自适应正则化反演中的八层验证模型，这里利用此模型将本方法与传统拉格朗日乘子二分法的搜索过程进行了详细地比较。四个模型的反演目标拟合差均为0.001。为便于叙述，将作者所描述的搜索方法称为“对数二次法搜索”，将传统拉格朗日乘子二分法的搜索称为“自然二分法搜索”。

表1为八层模型在首次迭代中的μ值变化过程，对数二分法只经过4次区间搜索加3次最小值搜索就完成了迭代，而自然二分法经历了9次区间搜索加5次最小值搜索才完成迭代。而从搜索的结果来看，两种方法搜索到的μ值基本上没有差异，甚至本例中对数二次法搜索得到的绝对拟合误差还更小一点。可见对数二分法在搜索上的优势，这种优势在首次迭代中体现得最为明显。

图2 地电模型及Occam反演结果Fig.2 The geoelectric model and Occam inversion results （a）模型1；（b）模型2；（c）模型3；（d）模型4

表1 模型1首次迭代μ值搜索过程Tab.1 The searching process of Lagrange multiplierin first iteration of model 1

表2、表3分别为模型一Occam反演用对数二次法搜索与自然二分法搜索的详细迭代过程。结果可以看出，对数二次法区间搜索往往只需要3次～4次，而3次搜索已经是区间搜索的极限；二次最小值搜索基本上也只需要1次～3次。整体对比，每次迭代对数二次法搜索的次数都明显少于自然二分法搜索。

表4为四种模型分别用两方搜索法进行Occam反演的μ值变化次数对比。无一例外，对数二次法搜索次数均少于自然二分法搜索，减小比例在20% ～50%之间。

从图2反演结果上看，本方法在提高μ值搜索效率的同时，对Occam反演的结果没有任何影响。

表2 模型1每次迭代对数二次法搜索次数Tab.2 The search process of each iteration of model 1by logarithmic quadratic method

表3 模型1每次迭代自然二分法搜索次数Tab.3 The search process of each iteration of model 1by natural dichotomymethod

表4 四种模型μ值变化总次数对比Tab.4 The comparison of the number of Lagrange multiplier's search

4 结论

作者提出了一种Occam反演中拉格朗日乘子的搜索方法，此方法主要做了以下两方面的改进。

1）利用关于拉格朗日乘子的绝对拟合误差曲线的形态特征，用公式μ＝10λ对拉格朗日乘子做个变化，从而将对μ值的搜索变为对λ值的搜索，即将拉格朗日乘子的搜索从自然数域中转化到了以10为底的对数域中进行。在运用过程中，发现对λ进行区间搜索的确很快，但也变相地加大了区间最小值的搜索范围，在实际算例中，整体的搜索效率并没有提高多少，故而进行了第二次改进。

2）利用区间搜索到的已知三点，构建关于λ值的拟合差二次函数，并求得此二次函数取得极小值时的λ值；用此λ值代入反演中求得拟合差，再得到一个新的谷值区间。这样反复进行二次函数构建，寻找区间最小值。在实际算例中，我们发现往往只需要1次～3次二次函数构建就能搜索到区间的最小值，这大大地提高了区间最小值的搜索效率。

将以上两种改进运用到拉格朗日乘子的搜索中，用模型算例验证发现，此方法将搜索效率提高了20%～50%，大大减少了Occam法反演的运算时间。

本搜索方法仅仅在一维大地电磁反演中进行了验证，有待在二维、三维中运用。

［1］ 石应骏，刘国栋，吴广耀，等.大地电磁测深法教程［M］.北京：地展出版社，1985.

SHI Y J，LIU G D，WU G Y，Tutorial of MT sounding［M］.Beijing：Geological Publishing House，1985.（In Chinese）

［2］ 陈乐寿，刘任，王天生.大地电磁测深资料处理与解释［M］.北京：石油工业出版社，1989.

CHEN L S，LIU R，WANG T S.Data processing and interpretation of MT sounding［M］.Beijing：Petroleum Industry Press，1989.（In Chinese）

［3］ 王家映.地球物理反演理论［M］.北京：高等教育出版社，2002.

WANG J Y.Inverse theory in geophysics［M］.Beijing：Higher Education Press，2002.（In Chinese）

［4］ CONSTABLE S C，PARKER R L，CONSTABLE C G.Occam's inversion：A practical algorithm for generating smooth models from electromagnetic sounding data［J］.Geophysics，1987，52（3）：289－300.

［5］ DEGROOT－HEDLIN C D，CONSTABLE S C.Occam's Inversion to Generate Smooth Two Dimensional Models form MagnetotelluricData［J］.Geophysics，1990，55（12）：1613－1624.

［6］ SIRIPUNVARAPORN W，EGBERTG.An efficient data－subspaceinversion method for 2－D magnetotelluricdata［J］.Geophysics，2000，65（3）：791－803.

［7］ SIRIPUNVARAPORN W，UYESHIMA M，EGBERT G.Three dimensional inversion for network－magnetotelluricdata［J］.Earth Planets Space，2004，56（9）：893－902.

［8］ RODI W L，MACKIE R L.Nonlinear conjugate gradients algorithm for 2－D magnetotelluricinversion［J］.Geophysics，2001，66（3）：174－187.

［9］ 韩波，胡祥云，何展祥，等.大地电磁反演方法的数学分类［J］.石油地球物理勘探，2012，47（1）：177－187.

HAN B，HU X Y，HE ZH X，et al.Mathematical classification of magnetotelluric inversion methods［J］.Oil Geophysical Prospecting，2012，47（1）：177－187.（In Chinese）

［10］陈小斌，赵国泽，汤吉，等.大地电磁自适应正则化反演算法［J］.地球物理学报，2005，48（4）：937－946.

CHEN X B，ZHAO G Z，TANG J，et al.An adaptive regularized inversion algorithm for magnetotelluric data ［J］.Chinese J Geophys，2005，48（4）：937－946.（In Chinese）

［11］吴小平，徐果明.大地电磁数据的OCCAM反演改进［J］.地球物理学报，1998，41（4）：547－554.

WU X P，XU G M.Improvement of OCCAM＇s inversion for MT date［J］.Chinese J Geophys，1998，41 （4）：547－554.（In Chinese）

［12］朴英哲，李桐林，刘永亮.在大地电磁二维Occam反演中求取拉格朗日乘子方法改进［J］.吉林大学学报：地球科学版，2014，44（2）：660－667.

PAK Y C，LI T L，LIU Y L.Improvement of choosing lagrange multiplier on MT 2DOccam inversion ［J］.Journal of Jilin University：Earth Science Edition，2014，44（2）：660－667.（In Chinese）

［13］刘俊峰，邓居智，陈辉，等.一种用于Occam反演中搜索拉格朗日乘子的方法［J］.工程地球物理学报，2013，10（3）：344－350.

LIU J F，DENG J Z，CHEN H，et al.A method used for searching lagrange multiplier in Occam inversion ［J］.Chinese journal of engineering geophysics，2013，10（3）：344－350.（In Chinese）

Lagrange multiplier's search in one－dimensional magnetotelluric Occam's inversion

ZHANG Jun－tao1，ZHOU Jun1，WANG Xu－ben1，XIA Shi－bin1，ZHONG Hong－mei2
（1.Key Lab of Earth Exploration ＆information Techniques of Ministry of Education，Chengdu university of technology，Chengdu 610059，China；2.Sichuan Nuclear Industry Geological Survey Institute，Chengdu 610059，China）

In magnetotelluric inversion，because of its stability in inversion and the advantage of model resolution and so on，Occam method is widely used.But in each iteration，Occam＇s inversion needs to constantly search Lagrange multiplier.The searching efficiency of Lagrange multiplier to Occam＇s inversion speed plays a vital role.To improve the search efficiency of Lagrange multiplier，a new searching method was presented in this paper.First Lagrange multiplier of transfer from the natural number field to search for the 10logs base domain.Second in the search of interval minimum search method using quadratic function.Through a number of one－dimensional magnetotelluric inversions，Lagrange multiplier search efficiency is increased by 20%－50%，and the operation speed of Occam＇s inversion is greatly improved.

magnetotelluric；Occam inversion；Lagrange multiplier；logarithmic；quadratic function

P 631.3

：A

10.3969／j.issn.1001－1749.2015.06.03

1001－1749（2015）06－0687－06

2014－11－18改回日期：2014－12－14

地质调查项目（12120113095100）；重大仪器专项（2011YQ05006007）；四川省科技支撑计划项目（2013FZ0084）；四川省教育厅项目（13ZB0062）

张君涛（1989－），男，硕士，主要研究方向为电磁法数据处理与正反演解释，E－mail：Zjuntao＠foxmail.com。