基于对称多项式约束的重力反演

郭一豪，张志勇,2*，陈晓,2，李曼,2，周峰,2，谢尚平,2，曾志文，程三

1 东华理工大学地球物理与测控技术学院，南昌 330013 2 东华理工大学核资源与环境国家重点实验室，南昌 330013

0 引言

由于地球物理场的等效性、观测数据集的有限性、数据观测存在误差等问题,地球物理反演不可避免的具有多解性,稳定、高精度的反演一直是地球物理反演研究追求的重要目标.Tikhonov正则化理论（Tikhonov and Arsenin,1977）自提出至今,被公认为是可以减少反问题不适定性的有效方法之一.地球物理反问题,在考虑数据误差的基础上引入模型稳定函数进行正则化可以有效降低反问题多解性;同时,通过稳定函数将先验信息融入反演,也使反演更符合实际情况、提高了反演的可靠性.

地球物理正则化反演中模型稳定函数的设计与选择直接影响反演的最终效果,其构建是反演的关键环节之一.根据反演模型的特征,模型稳定函数通常可分为光滑约束和聚焦约束两类.早期地球物理反演研究主要使用最小模型约束函数（胡正旺等,2021）.Constable等（1987）构建了最大化光滑模型稳定函数,实现了大地电磁（Magnetotelluric,MT）数据的Occam反演,研究工作表明采用光滑约束可以在一定程度上防止数据过度拟合,获得与真实地下物性分布更吻合的结果.相较光滑反演约束,Rudin等（1992）、Vogel和Oman（1998）在开展图像信息提取研究中提出了一种改进的总变化量支撑模型稳定函数（Total Variation,TV）,用于从噪声数据中恢复图像.Farquharson和Oldenburg（1998）在一维电磁反演中探讨了L-1范数和L-2范数模型约束反演的鲁棒性,研究表明L-1范数具有构造分段或反演“块状”模型能力.Last和Kubik（1983）提出了最小支撑模型函数（Minimum Support,MS）约束并将其应用于重力反演中,结果表明该约束反演结果具有致密性好、密度对比均匀等地质特征.Portniaguine和Zhdanov（2002）将MS聚焦函数引入三维航空磁法反演,实测数据反演表明,使用MS聚焦函数可以获得清晰的、聚焦的反演模型.Zhdanov等（2004）三维重力梯度张量MS聚焦约束反演的研究表明,该方法可以显著提高反演结果的分辨率.Zhdanov等（2004）进一步研究了基于MS约束的三维MT反演,实现了非二次最小支撑稳定器转化为传统的二次最小范数稳定器,使反问题可以采用统一形式优化.Portniaguine和Zhdanov（1999）研究了最小梯度支撑模型函数（Minimum Gradient Support,MGS）并将其应用于重力反演研究中,对比分析了最光滑约束、总变化量支撑约束、最小支撑和最小梯度支撑的反演结果,研究表明MGS与加权函数相结合,有利于得到清晰、集中的反演模型.改进的TV、MS和MGS约束函数中均存在一个小的正实数,对于改进TV约束函数来说,正实数是为了确保TV在零点处可导；对于MS和MGS约束函数来说,引入的正实数被称为聚焦因子,该因子直接影响反演结果,其选择问题一直是聚焦反演研究的一个难点.Zhao等（2016）提出了一种指数型模型聚焦函数,重磁反演研究表明该函数不仅具有聚焦作用,而且更加简便、稳定,可以避免聚焦因子的选取及函数梯度计算困难等问题.

另外，模型参数物性变换函数的使用一方面可以结合先验物性范围对反演物性进行合理的约束;另一方面,从一定程度上也可起到平衡模型误差与数据误差的作用,降低正则化因子选择的难度.Habashy和Abubakar（2004）提出了用于电磁反演的变换函数,确保模型参数在有实际物理意义的取值区间.此外，Loke和Barker（1995）对直流电阻率、Li和Oldenburg（2003）对磁法、Commer和Newman（2008）对三维可控源电磁（Controled Source Electromagnetic Method,CSEM）、Key（2016）对二维海洋电磁反演问题的模型变换函数进行了研究，均取得了理想效果.

随着勘探程度的提高与物性研究的深入,取得的有效先验信息也越来越丰富.丰富精确的先验信息可以有效地减少反演多解性,提高反演解释的正确性.钻孔和测井的开展提供了局部物性范围和空间展布的双重信息,因此融合钻孔约束的反演被应用于地球物理反演中（Bobée et al.,2010;Yan et al.,2017）.Dell′Aversana（2001）将通过测井资料确定的电阻率、速度和密度的统计关系应用于顺序反演中,研究了地质情况复杂的逆冲断层带.周竹生和周熙襄（1993）将地震资料、测井信息和先验地质知识相结合提出了一种宽带约束反演方法,该方法能提供可与测井资料进行对比的、具有较高分辨率的绝对波阻抗参数,为岩性模拟奠定基础.相比测井及钻孔信息约束,基于标本测试等方式统计获取的信息仅具有物性范围约束能力,其往往难以直接应用于反演约束.

为了将任意多种先验岩石物性约束融入地球物理反演,本文提出一种对称多项式模型约束函数,该模型约束通过连乘形式将多种统计或参考物性引入反演.为实现反演目标函数优化求解,深入研究了该模型约束函数的梯度与海森矩阵；开展了模型和实测数据反演,并通过与其他形式正则化模型约束对比分析,验证了该约束函数的正确性和有效性.

1 对称多项式约束

1.1 对称多项式的定义

定义如下的模型稳定函数：

（1）

1.2 对称多项式的函数表示

式（1）中的φi为对称多项式（为方便讨论忽略单元体积项）,

+c1T（mi）k-1+…+ck-1T（mi）1+cn，

（2）

系数c具有通解（推导见附录）:

（3）

i=1,2,3,…,k.

（4）

其中,ci为基本对称多项式,可通过递归算法求解.

1.3 对称多项式函数曲线

为研究对称多项式函数的特点,将其与光滑和最小支撑两种地球物理中常用的模型约束函数进行对比分析.

光滑约束函数

（5）

最小支撑约束函数

（6）

其中,ε为小的正实数,称为聚焦因子.该因子是聚焦反演的关键参数,Zhdanov和Tolstaya（2004）使用一系列的聚焦因子进行了反演,获得了模型误差函数-聚焦因子曲线,选取曲线斜率最大点处的聚焦因子作为反演最佳参数,该方法需要在反演前进行大量的反演试算,计算效率较低.

为了定性比较基于对称多项式和常用的一些模型稳定函数,假设一个反演单元存在三种可能岩性,物性分别为-1、0、1,绘制模型约束值变化曲线如图1,图中红、绿、蓝、黑线分别为最小模型、聚焦因子为0.01的最小支撑、聚焦因子为0.1的最小支撑以及对称多项式约束函数值变化曲线.

图1 三种模型稳定函数曲线Fig.1 The curves of three types model constraint function

分析图1,基于对称多项式的约束函数（黑色曲线）当反演物性处在给定参考值的最大和最小值之间时，函数取值相对较小,当物性超过给定约束范围，则函数值会快速增大;当物性取值与参考物性一致时，该函数取得极小值.对比最小模型（红色曲线）、最小支撑（绿、蓝色曲线），多项式约束具有明显的物性分类与取值区间约束能力.

2 基于对称多项式的反演

2.1 目标函数

根据正则化反演理论,构造如下目标函数：

（7）

α与τ为用于平衡数据误差与模型误差的平衡因子,称为正则化因子.可采用经验和自适应选择两种方法确定.其中,经验法是利用一定的准则确定某个定值作为近似最佳正则化因子,该值一般在整个反演过程中保持不变.主要的定值方法包括:无偏风险估计法（Wahba,1990）,广义交叉验证法（Wahba,1977；Golub et al.,1979）和L曲线法（Hansen and O′Leary,1993；Hansen,1994）.自适应方法是在反演过程中根据一定的准则自适应调整正则化因子取值的方法.Zhdanov（2002）提出了一种利用初始数据拟合泛函和模型稳定泛函的比值为正则化因子初值,在数据拟合效率下降时对正则化因子进行衰减的自适应方法;陈小斌等（2005）研究了最平缓约束MT反演正则化因子的自适应选择算法.为简化研究,本文反演中α与τ的选择均使用经验法.

2.2 对称多项式函数的一阶导数向量

由公式（1）可见,对称多项式约束求和中的每一项仅与一个单元物性相关.因此有

（8）

其中,φi为关于mi的一元k次多项式函数.进一步扩展，得到对称多项式约束的导数向量为

（9）

2.3 对称多项式函数的二阶导数矩阵

结合一阶导数的推导和公式（1）可知，对称多项式约束函数对物性的二阶导数为

（10）

因此,对称多项式约束项的海森矩阵为

（11）

2.4 共轭梯度优化流程

为简化推导,参考Zhdanov（2015）,记m=T（m）.对目标函数（7）的共轭梯度优化过程进行推导.数据项残差

Rn=A[T-1（mn）]-dobs，

（12）

其中,T-1（·）表示物性变换函数的反函数,A[T-1（mn）]和dobs分别为计算数据和观测数据,n为当前迭代次数.目标函数式（7）的梯度

（13）

（14）

（15）

（16）

（17）

模型更新,

（18）

3 模型试验

分别采用光滑约束、最小支撑约束和对称多项式约束进行重力二、三维模型试验,取变换函数为T（m）=m.

3.1 模型试验一

设计如图2所示的剩余密度模型,在剩余密度为0 g·cm-3的背景中存在一个剩余密度为0.1 g·cm-3的块状异常体.

图2 设计剩余密度模型Fig.2 Synthetic residual density model

观测区域为0～10 km,测点距为500 m,光滑模型约束反演正则化因子选择0.5,MS聚焦约束反演正则化因子选择0.05,对称多项式约束反演取光滑模型约束项权重因子0.5、多项式约束项权重因子150,采用矩形网格剖分,模型单元中心位置见示意图3.反演结果见图4,其中图4a、4b、4c分别为光滑模型、最小支撑（ε=0.03）、对称多项式约束（参考物性为m*=[0,0.1]）的反演结果,图4d、4e、4f为与其对应的误差曲线,图4g为响应曲线.

图3 模型离散单元中心节点位置Fig.3 Element center position of discrete mesh

对比图4a和图4c可以发现,光滑反演结果纵向连续性好,但反演的剩余密度最大仅达到0.06 g·cm-3,与真实物性0.1 g·cm-3存在较大差距,而基于对称多项式约束的反演结果物性与真实模型更接近.对比图4b和图4c可以发现,最小支撑函数约束反演物性可以达到真实物性,但存在聚焦作用过强风险,导致反演异常体规模偏小.由图4d、4e和4f可见三种方案的反演误差函数随迭代次数逐渐下降并最终趋于平缓,光滑与对称多项式约束反演数据误差下降速度比最小支撑模型约束快,最后的拟合差也小于最小支撑模型约束.由图4g可见,三种约束方案均取得了较好的拟合效果,且光滑模型约束与对称多项式约束拟合效果优于最小支撑模型约束结果.

3.2 模型试验二

设计如图5所示的剩余密度模型,在剩余密度为0 g·cm-3的背景中存在剩余密度分别为0.1 g·cm-3和-0.1 g·cm-3的两个块状异常体.

观测区域为0～10 km,测点距为500 m,光滑反演正则化因子选择0.3,MS聚焦约束反演正则化因子选择0.05,对称多项式模型约束反演光滑项权重0.3、多项式项权重3000,单元剖分同图3,反演结果见图6,其中,图6a,6b,6c分别为光滑、最小支撑（ε=0.1）、对称多项式约束（参考物性为m*=[0,0.1,-0.1]）的反演结果,图6d,6e，6f为对应的反演误差曲线,图6g为响应曲线.

对比图6a和图6c可以发现,光滑反演结果连续性好,但反演边界模糊,物性仍然小于真实物性,而基于对称多项式约束的反演结果物性与真实模型更接近,且边界更清晰.对比图6b和图6c可以发现,最小支撑函数约束反演物性可以达到真实物性,但反演聚焦效果强,异常体规模偏小.分析图6d、6e和6f可以发现三种方案的反演误差函数随迭代次数逐渐下降并最终趋于平缓.分析图6g可以发现,三种方案的反演响应与实测响应基本吻合.

图4 合成模型反演（a） 光滑模型约束反演密度断面;（b） 最小支撑模型约束反演密度断面;（c） 对称多项式约束反演密度断面;（d） 光滑模型约束反演误差曲线;（e） 最小支撑模型约束反演误差曲线;（f） 对称多项式约束反演误差曲线;（g） 响应拟合曲线.Fig.4 Synthetic model inversion（a） The density section of smooth model constraint inversion;（b） The density section of minimum support model constraint inversion;（c） The density section of symmetric polynomial model constrain inversion;（d） The smooth model constraint inversion error curves;（e） The minimum support model constraint inversion error curves;（f） The symmetric polynomial model constraint inversion error curves;（g） Data fitting curves.

图5 含两个异常体的剩余密度模型Fig.5 Synthetic residual density model with two abnormal blocks

3.3 模型试验三

设计如图7所示的三维剩余密度模型,在剩余密度为0 g·cm-3的背景中存在两个剩余密度分别为0.1 g·cm-3和-0.1 g·cm-3的块状异常体.

观测区域为长宽各10 km的正方形,测点距为500 m,测线距500 m,使用均匀长方体网格剖分,网格长宽高分别为500 m、500 m和200 m,光滑约束反演正则化因子选择5.0,MS聚焦约束反演正则化因子选择0.001,对称多项式约束反演光滑项权重为5.0、多项式项权重为1000,反演结果见图8、9和10.其中,图8为光滑反演结果,图9为最小支撑（ε=0.1）反演结果,图10为对称多项式约束（参考物性为m*=[0,0.1,-0.1]）的反演结果.各图中（a）为反演模型,（b）为反演物性统计图,（c）为误差曲线图.

图6 两个异常体合成模型反演（a） 光滑模型约束反演密度断面;（b） 最小支撑模型约束反演密度断面;（c） 对称多项式约束反演密度断面;（d） 光滑模型约束反演误差曲线;（e） 最小支撑模型约束反演误差曲线;（f） 对称多项式约束反演误差曲线;（g） 响应拟合曲线.Fig.6 Synthetic model inversion with two abnormal blocks（a） The density section of smooth model constraint inversion;（b） The density section of minimum support model constraint inversion;（c） The density section of symmetric polynomial model constrain inversion;（d）The smooth model constraint inversion error curves;（e） The minimum support model constraint inversion error curves;（f） The symmetric polynomial model constraint inversion error curves;（g） Data fitting curves.

对比分析三种方案的反演结果,可以发现:光滑约束反演,在异常体两侧存在假异常,反演边界模糊,物性小于真实物性；最小支撑反演结果异常体更紧凑,异常体规模与真实模型更接近,但物性范围仅为[-0.07,0.07] g·cm-3，与真实物性依旧存在差异；基于对称多项式约束的反演结果物性值与真实模型更接近,且边界更清晰,但纵向分辨率略低,分析可能与光滑约束权重选择有关.由物性统计直方图可见,对称多项式约束形成三个物性统计中心,反映对称多项式约束有较好的聚焦效果.分析误差曲线图,可以发现三种方案的反演误差函数随迭代次数逐渐下降并最终趋于平缓.

图7 两个异常体三维剩余密度模型Fig.7 Synthetic 3D model with two residual density blocks

图8 光滑约束反演结果（a） 反演剩余密度模型;（b） 反演物性统计图;（c） 反演误差曲线.Fig.8 Smooth constraint inversion results（a） Inversion residual density model;（b） Statistical chart of inversion physical property;（c） Inversion error curves.

图9 最小支撑约束反演结果（a） 反演剩余密度模型;（b） 反演物性统计图;（c） 反演误差曲线.Fig.9 Minimum support constraint inversion results（a） Inversion residual density model;（b） Statistical chart of inversion physical property;（c） Inversion error curves.

图10 对称多项式约束反演结果（a） 反演剩余密度模型;（b） 反演物性统计图;（c） 反演误差曲线.Fig.10 Symmetric polynomial constrain inversion results（a） Inversion residual density model;（b） Statistical chart of inversion physical property;（c） Inversion error curves.

综合分析图4、图6、图8、图9和图10,可以发现,基于对称多项式约束的反演可以将任意数量的物性统计信息融入反演,并能获得类似聚焦反演的效果,反演物性处于给定参考值范围内.

4 实测数据反演

为了验证对称多项式约束函数的实用性,选择某测线实测资料进行反演验证.测线长9 km,点距200 m,已知研究深度内存在侏罗系、志留系、奥陶系地层,其中侏罗系以凝灰岩为主,标本测试密度为2.5～2.55 g·cm-3,志留系以砂岩、砂质泥岩及页岩为主，密度为2.65～2.67 g·cm-3,奥陶系以钙质泥岩、灰岩为主，密度为2.68～2.7 g·cm-3.

分别使用光滑约束、最小支撑约束（ε=0.05）和对称多项式约束（参考物性为m*=[2.53,2.66,2.7]）进行反演.计算采用矩形网格剖分,单元中心节点位置示意图见图11.为便于与先验物性对比,采用绝对密度表示反演结果,见图12.其中,图12a、12b、12c分别为光滑、最小支撑、对称多项式约束反演结果,其中黑线是根据反演模型给出的推测地层界面；图12d、12e、12f为对应的反演误差曲线,图12g为响应曲线.

对比图12a、12b和12c可以发现，光滑反演可以获得平滑过渡的反演模型,但反演物性相对已知物性普遍偏小,边界分辨能力很差.最小支撑反演可以获得浅层物性分布信息,但对深部信息反演存在不足.对称多项式约束反演结果具有清晰边界且反演物性接近给定的先验物性信息.分析图12d、12e、12f和12g可以发现,三种反演方案数据误差均稳定下降并趋于平缓,观测数据与计算数据拟合程度很高.

为了直观比较三种方案的反演物性分布,绘制反演物性的频数分布直方图,图13a、13b和13c分别为光滑、最小支撑和对称多项式约束反演的物性统计结果,图中的蓝色虚线为先验的密度统计值.

从反演的频数分布直方图中可以发现,光滑反演物性存在明显的过渡,物性连续性好.MS反演结果物性分布相对孤立,反演结果存在明显的“分块”.但光滑反演和MS反演结果物性分布峰值与先验的物性统计信息存在明显偏差.基于对称多项式约束的反演结果物性分布峰值与先验统计物性更加接近,一定程度证明了对称多项式的约束能力.

实测数据反演表明对称多项式约束函数可以充分利用先验物性信息,通过将任意多个先验物性信息融入目标函数,一定程度上克服了复杂区域先验统计信息应用困难的问题.

5 结论与展望

本文提出了一种基于对称多项式形式的模型约束函数,并将其应用于重力数据正则化反演,通过模型试验和实测数据处理，形成以下结论.

（1）对称多项式可以将任意多个先验物性信息融入目标函数,一定程度上克服了复杂区域先验统计信息应用困难的问题,提高了反演的实际应用效果.

图11 模型离散单元中心节点位置Fig.11 Element center position of discrete mesh

（2）对称多项式模型约束函数极小值取值明确,反演具有较强的聚焦能力.

（3）在给定参考物性时,目标函数在最小与最大参考物性范围外，对称多项式约束函数取值快速增加,表明该约束隐含有物性范围约束能力.

（4）对称多项式模型约束函数的梯度与海森矩阵具有解析形式,反演目标函数优化求解简单.

此外,智能算法在地球物理中的应用受到了越来越多的关注,以聚类算法为代表的反演方法已成功应用于地球物理反演.大多数聚类算法可以有效地获得聚类中心,但没有明确的代价函数,因此当前只有具有明确代价函数的模糊C均值聚类在地球物理反演中获得了广泛的应用.基于对称多项式约束可以有效地结合各种聚类算法获得的聚类中心,实现聚类反演,如何将对称多项式与聚类算法相结合将是今后的研究方向.

致谢感谢审稿专家对本文详细专业的修改及编辑部的支持.

附录

设有一函数为

（A1）

可以发现:

（A2）

因为

1≤i

（A3）

所以Φ为对称多项式,有

（A4）

其中,σi为基本对称多项式,其表达式如下:

（A5）

图12 实测数据反演（a） 光滑模型约束反演密度断面;（b） 最小支撑模型约束反演密度断面;（c） 对称多项式约束反演密度断面;（d） 光滑模型约束反演误差曲线;（e） 最小支撑模型约束反演误差曲线;（f） 对称多项式约束反演误差曲线;（g） 响应拟合曲线.Fig.12 Field data inversion（a） The density section of smooth model constraint inversion;（b） The density section of minimum support model constraint inversion;（c） The density section of symmetric polynomial model constrain inversion;（d） The smooth model constraint inversion error curves;（e） The minimum support model constraint inversion error curves;（f） The symmetric polynomial model constraint inversion error curves;（g） Response fitting curves.

图13 重力反演密度统计图（a） 光滑反演;（b） 最小支撑反演;（c） 对称多项式约束反演.Fig.13 Gravity inversion density statistics（a） Smooth inversion;（b） Minimum support inversion;（c） Symmetric polynomial constrained inversion.

代入φ（A+Bx）有:

φ（A+Bx）=c0xk+c1xk-1+…+ck-1x1+ck.

（A6）

根据函数定义结合基本对称多项式定义可知,c的通解为

（A7）

（A8）

最终,函数f可以表示为如下的多项式函数乘积形式:

f=（c0xk+c1xk-1+…+ck-1x1+ck）

×（c0xk+c1xk-1+…+ck-1x1+ck）.

（A9）