基于UNet++卷积神经网络的重力异常三维密度反演

李柏森,鲁宝亮,2,3*,安国强,巨鹏,2,3,朱武,3,4,苏子旺

1 长安大学地质工程与测绘学院,西安 710054 2 海洋油气勘探国家工程研究中心,北京 100028 3 长安大学西部矿产资源与地质工程教育部重点实验室,西安 710054 4 自然资源部生态地质与灾害防控重点实验室,西安 710054

0 引言

三维物性结构是近年来地球物理反演的热点.而重磁方法作为勘探地球物理研究中的重要方向,是研究地壳结构与密度、磁性特征的主要方法,也是地球深部探测向全三维、大深度发展的重要方向.故通过建立合理的地球内部结构及物性赋存状态的地质模型,进行重磁三维物性反演问题研究,对解决地球内部物理性质、结构、构造以及能源矿产分布等问题具有重要的科学研究意义.

利用重力异常进行三维物性反演主要采用基于直立长方体单元剖分的正则化反演(Aghaee et al.,2021; Chasseriau and Chouteau,2003; Codd et al.,2021; Commer,2011; Dramsch et al.,2021; Geng et al.,2020; Guo et al.,2021; Hou et al.,2016; Li and Oldenburg,1996,1998,2003; Meng et al.,2016; Nagihara and Hall,2001; Pilkington,2006; Sun and Li,2014; Sun et al.,2018; Wang et al.,2017; Zhdanov et al.,2004; 刘银萍等,2013; 秦朋波和黄大年,2016; 姚长利等,2007).该方法首先将地下场源区域规则剖分为长方体单元,其次建立目标函数(含引入模型约束项),通过求解目标函数来反演确定剖分单元的密度.其中,模型约束有深度加权约束(Commer,2011; Li and Oldenburg,1996,1998,2003)、模型特征约束,其中分为光滑约束(Li and Oldenburg,1998,2003)、聚焦约束(Commer,2011)、稀疏约束(Pilkington,2006)、物性范围约束(Commer,2011; Li and Oldenburg,2003)等,这些模型约束条件的加入对正则化问题的解有很大的改善,但是模型约束中正则化参数以及经验选取的约束因子等对反演结果带来了较大的影响和不确定性.同时在求解最优解的过程中也经常容易陷入局部极值状态,无法得到全局最优解.

被证明可以拟合任意函数以及可以进行大数据处理的深度学习逐渐被应用于地球物理领域(胡琪鑫和徐亚,2022;刘俊等,2022;欧炳霖等,2023;王玥天等,2023;朱振宇等,2023).而大数据时代的到来以及近几年硬件技术的发展,使得神经网络迎来了新一轮研究和应用热潮.对地球物理来说,深度学习可以不依赖先验模型以及约束项、不需要存储模型矩阵、极强的非线性复杂函数逼近性能正是实际处理当中所需要的(Yu and Ma,2021; 张志厚等,2021a,2022).

目前深度学习神经网络在地震、测井、电法或者岩性分类方面的应用较为广泛,而这些成果也为深度学习在重磁领域应用提供了思路.国内外有很多学者使用不同的网络结构,例如BP神经网络、卷积神经网络(Convolutional Neural Networks,CNN)、全卷积神经网络(Fully Convolutional Neural Networks,FCN)等,对数据网格化(马国庆等,2021,2022)、二维(耿美霞和杨庆节,2013)或三维物性(Codd et al.,2021; Elawadi et al.,2001; Guo et al.,2021; Huang et al.,2021,2022; Lv et al.,2023; Wang et al.,2021; Yang et al.,2022; Zhang L Z et al.,2022; Zhang S et al.,2022; 管志宁等,1998; 张志厚等,2021a,b)、界面(He et al.,2021; Santos et al.,2015)以及特征位置(Poulton et al.,1992)等进行了反演方法研究,均取得了良好的效果.

但是之前使用的卷积神经网络最后通常为全连接层,并且无法提取融合不同层次的数据特征,使得应用场景及准确率受限; 其次,对重力数据来说,通常在正负半区均有分布,且绝对值相对较大,所以之前使用的tanh、sigmoid或者ReLU会存在梯度消失问题; 此外学习率通常需要随着训练进行动态调整,而传统使用的学习率衰减方法会使得陷入局部最优无法跳出,无法进行全局寻优.因此,不同的网络结构、激活函数以及学习率对反演具有重要影响,在应用过程中如何选取仍需进一步开展研究.

在前人的研究基础上,本文研究了基于改进的UNet++网络的重力异常三维密度深度学习反演方法,使用了对输入数据绝对值较大时梯度更为稳定、输入值为负时输出不为0的LeakyReLU激活函数,加入Bacth Normalization层对每层数据进行归一化操作,使其分布与输入数据相同,以增加网络训练的稳定性以及收敛速度,并采用了基于余弦退火的学习率变更策略以提高网络的全局最优化能力,避免陷入局部极小值.反演流程主要包括模型数据集及标签集的生成、深度神经网络的构建及训练以及模型预测.通过组合模型的试验,与UNet、FCN以及正则化方法对比,验证了本文方法具有较强的抗噪性及泛化能力,并通过对南澳大利亚Olympic Dam多金属矿区的实测重力数据反演,进一步证明了本文方法的正确性以及有效性.

1 方法原理

1.1 UNet++网络结构

神经网络通过大量数据集和对应标签的学习,来逼近数据集与标签之间的非线性关系.对数据来说,通常需要提取不同层次的特征,才能更全面地学习逼近.对于传统卷积神经网络来说,最后通常为全连接层,就使得应用范围受限.而随着全卷积神经网络的提出,通过卷积层和池化层对输入数据进行处理和下采样,学习不同层次的特征,而后通过反卷积层将其恢复至原始数据大小,这一系列也被称为编码-解码架构.

UNet是最常用的编码器-解码器网络,也是一种被使用较多的“端到端”网络,已广泛应用于各种场景(Ronneberger et al.,2015).过去大多数研究都将其作为主干,并对不同的任务进行了一些更改.UNet++是最具代表性的基于UNet的分割架构之一(Zhou et al.,2018).为了减少来自编码器和解码器子网络的特征图之间的语义差距,UNet++使用了一系列嵌套和密集的跳过路径(图1),而不仅仅是编码器和解码器网络之间的连接.与UNet降采样到最底层才返回不同,UNet++每进行一次降采样,就将其上采样回去并进行结果预测,最后再将不同尺度特征的预测结果进行融合输出预测结果,从而使得可以提取不同大小层次的特征.相较于UNet,UNet++填满了原本长链接部分的空白,通过抓取不同层次的特征进行整合叠加,从而实现对不同大小的目标的识别,提高网络精度.

本文基于UNet++网络架构搭建了如图1所示网络结构,网络共有37层,其中包含18个卷积层(卷积核尺寸3×3,跨度1)、4个最大池化层(宽度2,跨度2)、9个反卷积层(卷积核尺寸2×2,跨度2)、4个跳跃链接层以及输入输出层,相对于CNN来说未使用全连接层,既极大减少了节点参数,也减少了二维数据一维化导致的数据特征丢失,提高了计算效率以及反演精度.

为了提高训练网络的稳定性、增加神经网络的泛化能力,通常还会加入Dropout层或Batch Normalization层.Dropout层是指在训练过程中以某个概率暂时地使某些神经元停止工作,使得每次训练都类似于训练不同的网络结构,进而提高网络的泛化能力.但全卷积神经网络在执行过程中深层本身链接较少,使用Dropout放弃部分链接可能使得有效信息无法提取到.在训练过程中,每个batch的训练都会使得网络的参数进行更新.在神经网络较深时,后面数据的学习会受到前面数据分布的影响,而随着深度增加,受影响程度会越来越大,此时就会使得每层学习的数据分布不同,网络收敛速度变慢,泛化能力变差.本文通过加入Batch Normalization层对每层数据进行归一化操作,使其分布与输入数据相同,通过防止梯度爆炸或消失增加网络训练的稳定性、提高收敛速度(Ioffe and Szegedy,2015).

1.2 激活函数的选择

激活函数是通过作用在神经元上,将输入数据非线性地映射到输出端,进而提高网络的非线性拟合能力,使得网络能力更强,可以拟合更复杂的数据,从而表征输入输出之间更复杂的函数关系.目前常用的激活函数有sigmoid(乙状函数)、tanh(双曲正切函数)、ReLU(Rectified Linear Unit,整流线性单元)以及LeakyReLU(Leaky Rectified Linear Unit,渗漏整流线性单元),其具体函数图像见图2.

图2 激活函数曲线示意图

神经网络中参数的训练更新主要依赖反向传播算法,其核心内容是结合损失函数及最优化方法,而梯度类算法是目前应用最为广泛的一种,在训练过程中需要计算每一次残差的梯度值进行更新.对重力数据来说,数值在正负半区均有分布,且变化较大,绝对值一般分布在几十到上百.从函数图像可以看出,在输入绝对值大于2时,sigmoid、tanh变化平缓,梯度趋近于0,ReLU在输入值为负时输出恒为0,梯度为0,从而使得反向传播时无法更新或更新较小,增加了收敛难度; 由于LeakyReLU通过解决输入值为负时输出恒为0,减少了梯度消失的可能性,所以其具有更好的收敛能力以及更小的计算成本(Xu et al,2015).本文使用表现更好的LeakyReLU作为网络中的激活函数,并将在下文第2.2节中开展tanh、sigmoid、LeakyReLU对训练结果影响的实验研究.

1.3 基于余弦退火的学习率衰减

对神经网络训练来说,每轮次的更新步长会对最终训练结果以及收敛速度产生较大影响.当使用梯度下降算法来优化目标函数的时候,在训练初期学习率大一些,使得网络收敛迅速;在训练后期学习率小一些,使得网络更好的收敛到最优解;当越来越接近Loss值的全局最小值时,学习率应该变得更小来使得模型尽可能接近这一点.

对于传统衰减策略来说,通常是指定衰减方法,如线性衰减(Exponential Decay)、指数衰减(Linear Decay)等,使得学习率随训练过程逐渐减小从而使得结果收敛.因为目标优化函数可能不是严格的凸函数,所以可能存在多个局部最优解,而由于初始权重的随机性,在训练时梯度下降算法可能陷入局部最小值.此时可以通过突然提高学习率,来“跳出”局部最小值并找到通向全局最小值的路径,所以本文还进行了余弦退火(Cosine Annealing)学习率的应用.通过模拟余弦函数,先快速下降,再线性上升,通过热重启的办法使得寻优有可能跳出局部极小值(图3),进而提高网络训练收敛的全局搜索能力(Loshchilov and Hutter,2017),计算公式为

(1)

图3 不同学习率衰减策略及其应用示意图

其中,i为退火循环次数,α为退火率,ηt为当前批次学习率,ηmin,ηmax为学习率的最大值与最小值,Tcur,Ti分别为上一次热重启后训练轮次与给定热重启间隔轮次.

从图3a可以看出,线性衰减和指数衰减均为单调减小至某个最小值后保持恒定,使得当陷入局部最小值时无法寻找到全局最优(图3b),而余弦退火学习率衰减,可以在某一步时增大学习率,进而通过增大步长跳出局部最优,增强了全局寻优能力,训练结束后保存所有轮次中的最优解.为进一步证明余弦退火的全局寻优能力,将在下文第2.2节中开展线性衰减、指数衰减、余弦退火对训练结果影响的实验研究.

1.4 损失函数及学习算法的选择

损失函数与传统地球物理反演方法最小化的目标函数相似,通过计算输出结果与标签集的误差,利用反向传播算法来修改网络权值等参数,进而达到拟合效果.本文使用平均绝对误差(MAE)与均方差(MSE)进行衡量:

(2)

强大的学习算法也是收敛至全局最优解必不可少的一环.核心算法之一是随机梯度下降法(Stochastic Gradient Descent,SGD)(Hinton and Salakhutdinov,2006),但是其收敛速度慢,全局最优化能力较差,所以后续通过加入一阶动量或二阶动量的方法改进出了SGDM(SGD with Momentum)(Qian,1999)、AdaGrad(Adaptive Gradient)(Duchi et al.,2011).而为了解决AdaGrad二阶动量持续累积导致的训练提前结束问题,不计算全部历史梯度而只关注一段时间内的梯度平均值,改进得到了RMSprop.而本文采用Kingma和Ba(2017)提出的Adam算法,通过同时使用一阶动量以及二阶动量,使得收敛速度变快,计算效率提升.

2 数据集生成与网络训练

2.1 数据集生成

对深度学习来说,通常需要大量的数据集.本文将地下半空间模型均匀剖分为直立六面体,在其中随机生成长宽高不等的直立六面体组合及台阶状、饼状模型(图4).通过给定剩余密度值以及六面体顶点坐标,即可利用正演公式(3)、(4)计算其引起的重力异常及其垂向一阶导数(Nagy et al.,2000).由于只计算非零密度值的异常体,所以计算速度较快.再将总数据集分为三部分:训练集、验证集、测试集,训练集是指神经网络用于学习训练网络的数据集;验证集用于确定网络结构以及调整网络中的超参数;测试集是用于检验网络最终性能与泛化能力.公式(3)、(4)为

图4 数据集模型示意图

gz=Gρ{-(x-ξ)ln[r+(y-η)]-(y-η)ln[r+(x-ξ)]

(3)

(4)

数据集生成是将地下半空间模型均匀剖分为64×64×20个立方体,立方体边长为100 m,在其中随机生成长宽高不等的直立六面体组合及饼状、不同厚度的倾斜板状体模型,训练集剩余密度在-1.0、-0.8、-0.6、0.6、0.8、1.0 g·cm-3中随机取值,共生成模型16000组,利用正演公式计算其引起的重力异常生成数据集,其中训练集、验证集和测试集比例为80∶15∶5,即分别包含12800、2400、800组数据,每轮次训练后进行一次验证.

2.2 激活函数、学习率衰减策略及不同网络结构的确定

由于不同学习参数会对训练结果产生一定影响,本文选取了对训练结果影响较大的三个参数进行实验,分别是激活函数、学习率更新策略以及网络结构.同样使用前文生成的数据集,batchsize固定为64,给定最大训练轮次200,同时若验证集精度连续50轮次提升小于1×10-5则终止,训练结束后保存所有轮次中最优的网络参数,即以验证集精度最高的网络参数作为训练结果并用于模型实验.本文实验硬件平台为CPU:i9-10900X,GPU:RTX2080Ti×2,RAM:128GB,软件平台为Python3.9.12、TensorFlow2.9.0_GPU版.

首先是选取了三种激活函数,分别为LeakyReLU、sigmoid、tanh,使用UNet++网络以及余弦退火学习率更新策略,不同激活函数的收敛曲线见图5.从训练集MSE收敛情况可以看出LeakyReLU同样训练轮次下收敛更快,收敛结果也优于sigmoid和tanh;从验证集曲线来看使用LeakyReLU波动更小,稳定性更高,过拟合可能性要小于其他两种.

图5 不同激活函数均方差收敛曲线

其次选取了线性衰减、指数衰减以及余弦退火学习率更新进行实验,使用UNet++网络结构,LeakyReLU作为激活函数,最小学习率为1×10-7.MSE收敛曲线见图6.在训练刚开始阶段三种方法收敛速度接近,但随着训练进行指数衰减和线性衰减下降趋于平缓,而由于余弦退火具有热重启策略,收敛曲线具有起伏,但收敛情况优于线性衰减和指数衰减.

对于网络结构实验,选取了常用的三种网络模型,分别是UNet++、UNet以及FCN,使用LeakyReLU作为激活函数,以及余弦退火作为学习率更新方法,收敛曲线见图7.从收敛曲线可以看出UNet++在迭代相同轮次的情况下训练集和验证集MSE收敛情况都要优于UNet以及FCN.

图7 不同网络结构均方差收敛曲线

综上可以得出,UNet++网络模型结合LeakyReLU作为激活函数,在使用余弦退火学习率更新策略的情况下具有较强的全局寻优及防止过拟合能力.

3 模型实验

3.1 直立六面体组合

模型实验一为两个不同大小、密度、位置的直立六面体模型,具体参数及三维示意图如表1所示,其重力异常如图8a所示.为了测试噪声对反演结果的影响,对重力异常添加了其原始数据绝对值最大值的3%、5%、10%的正态分布随机噪声,噪声均值和方差见表2.

表1 直立六面体模型参数

表2 不同噪声水平及反演方法结果对比

图8 不同方法在不同噪声水平数据的反演结果

利用UNet++、FCN、UNet网络结构及传统方法进行反演对比,不同的网络结构均选取LeakyReLU激活函数、余弦退火学习率,即图7训练的网络,传统方法反演使用正则化方法,加入深度加权、物性约束以及最小长度约束,使用L2范数作为目标函数.反演结果的均方根误差见表2,三维结果见图8e—t.结果表明,在不加噪的情况下UNet++与FCN反演结果接近,优于UNet,能够非常准确地刻画真实模型的地下空间位置以及相对剩余密度大小,边界及异常体形态清晰,且明显优于传统正则化方法.不同噪声水平下的反演结果的均方根误差(表2)及三维几何形态(图8),可以看出UNet++均优于其他两种网络,也仍然优于传统方法,说明UNet++受噪声影响较小,更为稳定,抗噪性较好.

3.2 复杂组合模型

该实验为三种模型的组合,分别为两个直立六面体、两个饼状模型、一个台阶状模型,其具体参数及三维示意图如表3所示,不同噪声水平的重力异常见图9a—d,噪声的均值及方差见表4.

表3 复杂模型中各模型参数

表4 不同噪声水平及反演方法结果对比

图9 不同方法在不同噪声水平数据的反演结果

不同网络结构及传统正则化方法反演结果见图9e—t,均方根误差见表4.由于该种模型组合未经训练,所以反演结果均方根误差均大于3.1节中直立六面体模型.在无噪声情况下,UNet++网络反演结果仍然可以较为准确的刻画异常体形态特征及剩余密度情况,与FCN和UNet结果接近,均优于传统正则化方法.随着噪声水平增加,三种网络结构对饼状模型的反演逐渐变差,但UNet++网络对边界刻画明显优于其他两种网络,FCN和传统正则化方法对饼状模型较难识别,UNet受噪声干扰较大,虚假异常较多.

综合模型实验结果,从无噪声的反演结果来看,本文提出的方法模型可以应用于重力异常的三维密度反演,且即使未经过训练的组合模型也可以具有一定的反演效果,说明了该网络具有较强的泛化能力,但反演精度会低于经过训练的模型组合;从添加噪声后的反演结果来看,即使在添加10%的噪声后,也可以反演得到异常体模型的空间位置及相对剩余密度情况,且反演结果均优于FCN、UNet以及传统正则化方法,说明本文所搭建的网络具有一定的有效性及稳定性.

4 实测数据测试

为进一步检验本文所搭建网络的可行性及实用性,利用Olympic Dam多金属矿区布格重力数据进行了反演实验.Olympic Dam多金属矿区位于澳大利亚南部高勒克拉通省(图10),是20世纪70年代发现的一个超大型Cu-U-Au-Ag矿床,同时也是典型的IOCG矿床(Iron Oxide-Copper-Gold矿床).矿床产于一个复杂的角砾岩群中,主要矿石为富含赤铁矿和磁铁矿的角砾岩(图11).Olympic Dam多金属矿床具有明显的重力异常特征.为了消除背景场,获得由矿体引起的剩余重力异常,使用二维小波多尺度分解方法(杨文采等,2001)对该区域的布格重力异常数据进行异常的分离,如图12a所示为重力异常,图12b为区域重力异常,图12c为剩余重力异常.

图10 高勒克拉通省的基底地质图(Direen and Lyons,2007)

图11 超巨型Olympic Dam角砾岩群的简化地质图及剖面图(Kirchenbaur et al.,2016)

图12 重力异常及其分离结果

利用本文所提方法及传统正则化方法对该区剩余重力异常进行了三维密度反演,两种方法反演得到的高密度岩体平面位置分布重合度较高,本文方法反演的岩体更为聚焦,边缘更清晰,而正则化反演更为光滑.从本文方法反演结果(图13a)来看,高密度岩体平面位置分布与剩余重力异常分布形态相似,纵向上主要分布在300～1400 m的深度范围,并与Fe元素含量5%等值线(图11a)吻合度较高.4口钻井岩心实测密度曲线(图14)显示高密度岩体(实测密度大于2.9 g·cm-3)主要分布集中在400～1400 m范围内;浅部上覆盖层为低密度,小于2.6 g·cm-3,深部为一套变质的火山岩、沉积岩及侵入岩,密度相对矿体较小,在2.6～2.9 g·cm-3之间(Direen 和 Lyons,2007),反演结果与钻井RD16W1井密度变化较为一致(图15),其中该井在455～1316 m深度部分的密度为实测值,其余深度部分密度根据Direen和Lyons(2007)所给数据推测绘制.另外,本文方法反演结果与前人反演结果(Yang et al.,2019; Zhang S et al.,2022)及正则化反演结果(图13c)的平面形态相似,但最大分布深度略有差别.以上分析表明本文方法反演结果水平位置及垂向深度与目标岩体位置重合度较高.另外,分别对深度学习及正则化方法反演结果进行了重力数据重建(图16a、d),计算了重建数据与剩余重力异常的均方根误差(图16b、e)及残差(图16c、f),结果显示两种方法重建结果与剩余重力异常形态相似,重建效果均较好,UNet++重建结果均方根误差为0.6433 mGal,略小于正则化方法重建结果.综上,说明本文所提方法在实际数据处理中也表现出良好的可行性和实用性.

图13 反演结果三维透视图及切片图

图14 钻井岩心实测密度曲线

图15 过RD16W1井密度剖面(y=6063.68 km)

图16 UNet++与传统方法反演结果重建数据及残差

5 结论

本文提出了基于UNet++网络的重力异常三维物性反演,采用LeakyReLU激活函数以及加入Batch Normalization层,提高了网络更新收敛速度和稳定性;引入基于余弦退火的学习率更新策略,明显提升了网络的全局最优化能力,进而实现了重力异常数据地下异常体模型的三维物性反演.将本文方法与UNet、FCN以及传统正则化方法应用于不同噪声水平模型实验,结果表明本文搭建的网络可以准确地反映地下异常体模型的三维物性,且未经训练的模型也可以较为准确的反映异常体特征,说明网络具有良好的泛化能力及较强的抗噪性.并将该方法应用于澳大利亚Olympic Dam多金属矿区的实际重力数据,反演结果与地质情况、钻井资料以及正则化反演结果吻合度较高,进一步表明本文所提出的方法具有良好的可行性及实用性.

但是深度学习方法通常需要足量、多样且具有代表性的样本,而增加样本量会使得计算数据量及复杂度大为上升,如何通过小批量样本训练网络达到相近结果是本文方法的一个局限性.本文反演未使用先验模型信息及约束项,仅仅使用了一种网络架构,而不同的网络具有不同的优势,也就会导致反演结果有所不同.在下一步工作中将优化网络结构及样本特征性,尝试不同方法,加入约束项及先验信息,进而增加网络收敛能力及应用范围.另外,在本文研究中发现采用深度学习方法的反演结果出现了垂向分层现象,在其他深度神经网络重力反演结果中也出现了该现象(Yang et al.,2022; Zhou et al.,2023),具体原因仍需深入研究.

致谢感谢两位匿名评审专家及编辑部提出的宝贵修改意见和建议,感谢南澳大利亚资源信息网站提供的Olympic Dam多金属矿区布格重力数据及钻井资料.