1、第 37 卷第 1 期2023 年 2 月现代地质GEOSCIENCEVol.37No.1Feb.，2023DOI:10.19657/j.geoscience.1000 8527.2022.074基于乘积型目标函数的电阻率法三维反演谢海军，谭捍东，马逢群(中国地质大学(北京)地球物理与信息技术学院，北京100083)收稿日期:2022-03-20;改回日期:2022-10-20。基金项目:国家自然科学基金项目(41830429);山西省重点研发计划项目(202102080301001)。作者简介:谢海军，男，硕士研究生，1998 年出生，地球探测与信息技术专业，主要从事地球物理算法研究。Ema
2、il:2010200055 。通信作者:谭捍东，男，教授，博士生导师，1966 年出生，地球探测与信息技术专业，主要从事地球物理算法研究。Email:thd 。摘要:为了减小地球物理反演的多解性，通常采用累加型目标函数，即数据拟合差项加上正则化项。基于累加型目标函数的反演存在如何优选正则化因子的问题，这个过程通常需要做大量的反演计算。本研究系统分析和实现了基于乘积型目标函数(数据拟合差项和正则化项相乘)的电阻率法三维反演，乘积型目标函数的反演不存在优选正则化因子的问题。三维正演采用非结构化有限单元法，三维反演采用有限内存拟牛顿方法。使用理论模型合成数据进行了三维反演，检验了基于乘积型目标函数的
3、电阻率法三维反演的可行性与有效性。反演算例结果表明:基于乘积型目标函数的反演方法能够可靠恢复异常体的电阻率值、形态和位置。关键词:电阻率法;正则化;乘积型目标函数;三维反演中图分类号:P631文献标志码:A文章编号:1000 8527(2023)01 0067 07Three-dimensional Inversion for esistivity Method Basedon Multiplicative egularizationXIE Haijun，TAN Handong，MA Fengqun(School of Geophysics and Information Technology
4、，China University of Geosciences，Beijing100083，China)Abstract:To reduce the multiplicity of geophysical inversion solution，additive objective functions are usuallyused，i.e.，the data term and regularization term Inversion based on additive objective function has the prob-lem of optimizing the regular
5、ization factor，which usually requires much inversion calculations This study sys-tematically analyzed and implemented the 3D resistivity inversion method，based on multiplicative objectivefunction(i.e.，multiplication of the data term and regularization term)The problem of optimizing regulariza-tion f
6、actor is absent in the inversion of multiplicative objective function We used the unstructured finite elementmethod for resolving the 3D forward problem，and the L-BFGS(Limited-memory Broyden-Flecher-Goldfarb-Shanno)method for 3D inversion The feasibility and effectiveness of our approach is verified
7、 by using the syn-thetic data of the theoretical model The results show that the inversion method(based on multiplicative objec-tive function)can reliably recover the resistivity value，shape and position of anomaliesKey words:resistivity method;regularization;multiplicative objective function;3D inv
8、ersion0引言电阻率法三维反演问题属于大规模非线性问题，同其他反演问题一样，具有固有的多解性，因此需要寻找一些方法减少反演问题的多解性。正则化方法1 2 正是解决上述问题的常用方法之一。正则化反演方法在目标函数中增加了正则化项进行约束，一般采用累加型目标函数，即数据拟合差项和正则化项相加，用正则化因子(拉格朗日因子)控制数据拟合差项和正则化项之间的权重。正则化因子取值越大，则正则化项起主导作用，反演结果更偏重于光滑模型，反之数据拟合差项起主导作用，则更偏重于数据拟合。只有选择适当的正则化因子才能得到可靠的反演结果，但适当的正则化因子需要选择不同的正则化因子进行反演计算，对比分析反演结果来确
9、定。这种选择正则化因子的做法相当耗时，是累加型目标函数的缺点。为了解决累加型目标函数优选正则化因子耗时的问题，基于 Bakushinskii3、Blaschke 等4 和Jin5 的研究成果，Van den berg 等人6 8 首先提出了基于乘积型目标函数的对比源反演(ContrastSource Inversion)方法。在 Van den berg 等人的研究中，全变分正则化项9 11 作为目标函数的乘法约束，这种做法的优点是基于乘积型目标函数的反演方法在反演过程中本身会以迭代方式自动调整数据拟合差项和正则化项之间的权重，而不需要优选正则化因子。乘积型目标函数在过去二十多年中受到了广泛关
10、 注，并 成 功 地 应 用 于 物 理 学 的 不 同 分支12 20。在地球物理反演领域，乘积型目标函数应用较少，代表性的成果有:Habashy 等21 基于乘积型目标函数开发了电磁测井三维反演算法;Zaslavsky 等22 基于乘积型目标函数实现了井间电磁法和可控源电磁(CSEM)法的三维反演算法;Abubakar23 基 于 乘 积 型 目 标 函 数 分 别 实 现 了CSEM 数据和大地电磁数据的二维和三维反演;Alpak 等24 基于乘积型目标函数实现了二维感应测井和井旁压力测量数据的联合反演。以上研究成果都是基于结构化网格剖分，目前还未见有基于乘积型目标函数的电阻率法三维反演
11、的成果。本文基于非结构化网格和乘积型目标函数，实现了电阻率法三维有限内存拟牛顿(L BFGS)反演，设计双异常体模型进行了合成数据反演，检验了该反演方法的有效性和可行性。1电阻率法三维正演如图 1 所示，将电导率 分解为背景电导率p和异常电导率s，无穷远边界为以及地表与空气的界面为s，点电源 A 的三维地电场电势满足的微分方程25 为:(x，y，z)(x，y，z)=I(x x0)(y y0)(z z0)(1)其中:I 为电流，为狄拉克函数，为电势，点电源位于r0=(x0，y0，z0)。图 1电阻率法三维模型示意图Fig.1Three-dimensional model of resistivi
12、ty method在地表与空气分界面采用第二类边界条件，在无穷远边界采用混合边界条件，则电阻率法边值问题表示为:n=0，s(2)n+rcos=0，(3)其中:n 为边界外法线方向的单位矢量，为 n 和r 的夹角。采用非结构化四面体网格(图 2)有限单元法26 27 对以上方程进行离散，得到待求解大型线性方程，其矩阵形式为:K=F(4)其中:K 为与网格剖分相关的总体系数矩阵，为节点上异常场电势，F 为场源项。求解此线性方程组就可得到异常场电势，将其由解析解计算的一次场电势相加得到所有节点的总场电势，即所有测点电势，从而求出所有测点处的视电阻率值。图 2非结构化网格示意图Fig.2Schemat
13、ic diagram of unstructured grid86现代地质2023 年2电阻率法三维反演2.1目标函数电阻率法乘积型目标函数定义为:(m)=d(m)m(m)(5)目标函数梯度为:g=2JTWTdWd dobs F(m)m(m)+d(m)Tmm(m)(6)目标函数的数据拟合差项 d(m)为:d(m)=dobs F(m)TV1 dobs F(m)(7)其中:m=m1，m2，mMT为模型参数向量，M是模型参数个数，dobs=d1，d2，dNT是观测数据向量，di是某一观测点上的电阻率值，N 表示观测数据总个数，F(m)为正演响应，V1为数据协方差矩阵，JT是雅可比矩阵。目标函数的正则
14、化项 m(m)为:m(m)=(m)2+(8)式中:为表示模型之间的梯度光滑约束的粗糙度算子28 29;为常数。图 3 为二维非结构化网格模型示意图，中心深色网格代表当前需要光滑的模型单元，黑点表示各个网格单元中心点，各浅色单元中心点到深色单元的距离用实线表示，深色网格单元与周边共顶点的浅色网格单元的距离成反比，与周边网格的体积成正比，此概念可类比到三维。据此在非结构化有限元网格上导出 L2范数定义的粗糙度算子，为使正则化项函数离散化，假设每个单元模型参数为常数，则粗糙度算子 的离散形式为:(m)2=mi=1ViN(i)i=1j(mij rij)2(9)mij=mi mj(10)rij=(xi
15、xj)2+(yi yj)2+(zi zj)2(11)j=VjN(i)k=1Vk(12)其中:mi为 第 i 个单元的模型参数，Vi为第 i 个单元的体积，wj为体现相邻元素间导电性差异的体积加权函数，N(i)为第 i 个单元所有相邻单元的集合，rij为模型参数的空间距离。此外，公式(8)中的 为常数，其作用是防止正则化项在迭代过程中为 0，导致迭代失败，由以下公式确定6 7，23:=NV(13)其中:N 表示观测数据总个数，V 表示剖分单元的尺寸。在非结构化网格中，每个四面体单元的尺寸无法具体算出，所以采取 V 为模型总尺寸除去剖分单元总个数得出的单元平均体积这种形式。需要注意的是，与正则化因
16、子(拉格朗日因子)的选择相比，正则化项对 的变化不敏感21，24，所以采用平均体积这种做法对反演影响很小。图 3二维非结构化网格模型向量相互作用关系(中心深色网格代表当前需要光滑的模型单元，黑点表示各个网格单元中心点)Fig.3Vector interaction relationship of 2D unstructuredmesh model2.2L BFGS 反演流程本文采用 L BFGS 实现电阻率法三维反演。L BFGS 来源于传统的牛顿法，它改进了高斯牛顿法，使反演更快速高效。L BFGS 中不再存储完整的海森矩阵，只存储计算中的向量序列，而且仅存储它们最新的 n 个数值，当需要海森矩阵时就用最新的 n 个值计算。海森矩阵的逆和迭代步长的求取，以及“拟正演”的详细求解过程可参考王堃鹏等30 的成果。L BFGS 反演流程如下:(1)选择反演初始模型、最小数据拟合差、最大迭代次数、常数 ;(2)计算目标函数(m)与梯度 g;(3)计算拟海森矩阵的逆 Hk1与拟牛顿方向 pk=Hk1g;(4)设定初始步长，寻找一个当前迭代的合适步长 a 来最小化目标函数;(5)更新模型 mk+