三维目标反演是地球物理、遥感、计算机视觉等领域中的一个重要研究方向. 在地球物理领域，重磁数据被用于反演地球的内部结构，在金属矿床和油气资源勘探中广泛应用. 其中，重力反演是指利用重力数据演算地下物质分布，目前反演方法可以分为基于物理模型的反演方法和基于深度学习的反演方法. 基于物理模型的反演方法依托地球物理理论，利用线性或者非线性优化算法，使目标函数最小化. 由于地球场固有的等效性以及观测数据的离散性和误差，重力反演被视为不适定性问题，反演结果具有多解性^[1]. 针对该问题，引入正则化约束，如最小体积约束^[2]、深度加权^[3]和最小梯度约束^[4]. 引入遗传算法^[5-6]、粒子群算法^[7-8]和蚁群算法^[9]等全局优化算法，解决了陷入局部最小值的问题，但带来了计算复杂、耗时的问题. 基于深度学习的反演方法通过数据驱动的方式建立重力数据与地下目标的映射关系，当输入的重力数据确定时，反演结果唯一. 在地球物理领域，迄今为止已有多种结构的神经网络被用于重力反演. 王逸宸等^[10]使用卷积神经网络，实现了重力等值线图到异常体大小和埋深信息的映射. 阳前果^[11]基于U-Net^[12]重新设计了Grav-Net，端到端地实现了多异常体含噪数据的三维反演，达到了优于多约束正则化反演的效果. 为了提升三维目标的反演效果，Wang等^[13]对3D U-Net进行改进，使用三维卷积实现二维异常数据到三维异常体密度分布的映射. Zhang等^[14]设计GravInv，利用二维卷积和残差结构进行特征提取，维度变换后使用三维反卷积实现模型的重构. 三维卷积的引入使得网络训练和预测所需的内存和时间成本大幅增加. Zhang等^[15]通过预测异常体模型的特征，设计DecNet，间接地实现了三维重力反演.

因为基于深度学习的反演方法采用较多的三维卷积运算，反演网络存在内存占用大、训练耗时久的问题. 针对该问题，本文提出面向重力反演领域的基于多特征重构的三维目标反演算法，旨在降低反演网络的训练成本，提高反演结果的准确性. 主要工作如下. 1)通过特征分解提取目标的水平区域、中心深度、垂直厚度和剩余密度4类特征，实现三维模型的压缩. 2)设计多特征重构反演网络MRNet，实现对各类特征的预测及三维模型的重构. 3)引入梯度联合，通过原始数据和梯度数据的拼接，增强目标水平区域的特征. 4)引入注意力机制分化网络中不同解码器的特征预测功能，优化反演结果.

重力反演是指依据重力异常数据预测地下异常体密度模型，二者的关系满足：

(1)$ {\boldsymbol{a}} = {\boldsymbol{Gm}}. $

式中：${\boldsymbol{a}}$为重力异常数据；$ {\boldsymbol{m}} $为密度模型；${\boldsymbol{G}}$为核矩阵，表示地下密度模型到重力异常数据的映射关系. 反演问题的核心是求解${{\boldsymbol{G}}^{ - 1}}$.

基于深度学习的反演是数据驱动的方法，利用大量的先验数据优化参数化函数(神经网络)，构建重力异常数据与密度模型的映射关系. 神经网络的基本形式可以表示为

(2)$ \begin{split} \hat {\boldsymbol{m}} =&\; {\mathrm{Net}}\;({\boldsymbol{a}};\theta )= \\ & {f_n}({{\boldsymbol{W}}_n}{f_{n - 1}}( \cdots ({f_2}({{\boldsymbol{W}}_2}{f_1}({{\boldsymbol{W}}_1}{\boldsymbol{a}} + {{\boldsymbol{p}}_1}) + {{\boldsymbol{p}}_2})) \cdots ) + {{\boldsymbol{p}}_n}). \end{split}$

式中：$\hat {\boldsymbol{m}}$为密度模型的预测值；$ {\mathrm{Net}} $为非线性函数，表示神经网络；$\theta $为神经网络中包括${{\boldsymbol{W}}_i}$和${{\boldsymbol{p}}_i}$在内的待优化参数，二者分别表示权重和偏置；${f_i}$为激活函数. 实现基于深度学习的反演关键在于优化神经网络，即求解满足以下方程的最优解：

(3)$ \tilde \theta = {\mathrm{arg}}\mathop {\min }\limits_\theta \frac{1}{N}\mathop \sum \nolimits_{i = 1}^N L({{\boldsymbol{m}}_i},{\mathrm{Net}}\;({{\boldsymbol{a}}_i};\theta )). $

式中：$N$为训练样本数，$({{\boldsymbol{a}}_i},{{\boldsymbol{m}}_i})$为第$i$对训练样本，$L$为损失函数. 本文选择L2损失函数，表达式为

(4)$ L= \frac{1}{N}{{\displaystyle \sum }}_{i=1}^{N}{\Vert {{\boldsymbol{m}}}_{i}-{\mathrm{Net}}\;({{\boldsymbol{a}}}_{i};\theta )\Vert }_{2}^{2}. $

常用的神经网络优化算法有随机梯度下降、小批量随机梯度下降、动量法、AdaGrad和Adam等优化算法^[16].

为了使神经网络直接实现二维数据到三维模型的映射，简单增加输出通道数或三维卷积的方式会大大增加网络的参数量和内存占用. 针对该问题，通过特征分解实现对三维模型的压缩，将三维反演问题分解为对多个二维特征的预测. 设异常体的空间密度模型为$ {\boldsymbol{m}} $，形状为$ H \times W \times L $，则水平区域、中心深度、垂直厚度和剩余密度4类二维异常特征如下.

(5)$ {{r}}(x,y) = \left\{ \begin{gathered} 1 , \quad {\max} _{z = 1}^L{{m}}(x,y,z) > 0; \\ 0, \quad 其他. \\ \end{gathered} \right. $

(6)$ {{c}}(x,y) = \left\{ \begin{array}{ll} {0.5({\max} _{z = 1}^L\{ z|{{m}}(x,y,z) \ne 0\} } {+{\min} _{z = 1}^L\{ z|{{m}}(x,y,z) \ne 0\} )} , & {\max} _{z = 1}^L{{m}}(x,y,z) > 0 ; \\ 0 , & 其他. \\ \end{array} \right. $

(7)$ {{t}}(x,y) =\left\{ \begin{array}{ll} {\max _{z = 1}^L\{ z|{{m}}(x,y,z) \ne 0\} } { - \min _{z = 1}^L\{ z|{{m}}(x,y,z) \ne 0\} } , &{\max} _{z = 1}^L{{m}}(x,y,z) > 0; \\ 0 , &其他. \end{array} \right.$

(8)$ {{d}}(x,y) = {\max} _{z = 1}^L{{m}}(x,y,z). $

式中：$ 0 \leqslant x \leqslant H - 1 $，$ 0 \leqslant y \leqslant W - 1 $，$ 0 \leqslant z \leqslant L - 1 $；$ {\boldsymbol{r}} $、${\boldsymbol{c}}$、$ {\boldsymbol{t}} $和$ {\boldsymbol{d}} $分别为水平区域、中心深度、垂直厚度和剩余密度的异常特征，尺寸都为$ H \times W $. 假设以大小为4×4×4的矩阵表示异常体密度模型，矩阵中心大小为2×2×2的区域存在剩余密度为0.6的异常体. 4类异常特征如图1所示. 水平区域特征中值为1表示该区域存在异常体，中心深度和垂直厚度特征中的值分别表示异常体沿Z轴方向最大索引和最小索引的均值和差值. 剩余密度特征中的值表示该区域异常体的剩余密度.

当异常体所处深度较大或者边缘复杂曲折时，在地表边缘处引起的重力异常变化较小，反演网络往往需要增加网络深度，以提高对异常体边缘的敏感度. 针对该问题，分别计算二维重力异常数据行、列方向上的梯度数据，并与原重力异常数据进行通道方向上的拼接，得到3通道的数据矩阵作为反演网络的输入，如图2所示.

根据多特征图重构和重力梯度联合，提出基于U-Net基本架构的多任务反演网络(multi-feature reconstruction network, MRNet)，如图3所示. MRNet由特征提取模块和特征映射模块组成. 特征提取模块包含1个共享编码器，负责提取共享的多层次特征. 特征映射模块包含4个解码器，分别预测三维重力异常体的4类二维特征，重构得到三维模型.

网络中编码器和解码器的细节如图4所示，二者分别通过多次二维卷积和下(上)采样操作实现特征提取(映射)，且在共享编码器和每个解码器之间增加跨层连接. 为了分化不同解码器的功能，增强对应特征的预测能力，在跨层连接处嵌入CA(coordinate attention)模块^[17]，结构如图5所示. CA模块兼顾通道注意力和空间注意力，在训练过程中为不同的通道和位置进行加权，自适应地关注更重要的通道和位置特征. CA模块被用于网络输入端，实现对3通道数据矩阵的自适应加权.

MRNet将三维反演简化为对二维特征的预测，将尺寸为$ {N_x} \times {N_y} \times {N_z} $的三维异常体模型压缩为$ {N_x} \times {N_y} \times 4 $，其中$ {N_x} $、$ {N_y} $和$ {N_z} $分别为反演区域X、Y和Z轴方向上的网格数. 本文中$ {N_z} $为40，利用该方法可以降低90%的内存占用. 在相同的硬件水平下，更低的内存占用将加速网络训练和预测的过程.

在异常体的4类特征中，水平区域属于较浅层的特征，基于U-Net架构的网络较适合学习该特征. 在训练过程中，MRNet使用异常体实际的水平区域特征约束其他3类特征的学习，使特征映射模块只专注于该区域内部特征的学习. MRNet选择L2范数作为损失函数：

(9)$ \begin{split} L =& \frac{1}{N}\mathop \sum \limits_{i = 1}^N [{\lambda _{\text{r}}}\left\| {{D_{\text{r}}}(E({{\boldsymbol{g}}_i};{\theta _{\text{e}}});{\theta _{\text{r}}}) - {{\boldsymbol{r}}_i}} \right\|_2^2+ \\&{\lambda _{\text{c}}}\left\| {{D_{\text{c}}}(E({{\boldsymbol{g}}_i};{\theta _{\text{e}}});{\theta _{\text{c}}}) \cdot {{\boldsymbol{r}}_i} - {{\boldsymbol{c}}_i}} \right\|_2^2+ \\&{\lambda _{\text{t}}}\left\| {{D_{\text{t}}}(E({{\boldsymbol{g}}_i};{\theta _{\text{e}}});{\theta _{\text{t}}}) \cdot {{\boldsymbol{r}}_i} - {{\boldsymbol{t}}_i}} \right\|_2^2 + \\&{\lambda _{\text{d}}}\left\| {{D_{\text{d}}}(E({{\boldsymbol{g}}_i};{\theta _{\text{e}}});{\theta _{\text{d}}}) \cdot {{\boldsymbol{r}}_i} - {{\boldsymbol{d}}_i}} \right\|_2^2].\end{split} $

式中：$E$为编码器模块函数，${D_{\text{r}}}$、${D_{\text{c}}}$、${D_{\text{t}}}$和${D_{\text{d}}}$分别为水平区域、中心深度、垂直厚度和剩余密度支路上的解码器模块函数，${\theta _{\text{e}}}$、${\theta _{\text{r}}}$、${\theta _{\text{c}}}$、${\theta _{\text{t}}}$和${\theta _{\text{d}}}$分别为编码器模块和相应的解码器模块(水平区域、中心深度、垂直厚度和剩余密度)中待优化的参数，${{\boldsymbol{g}}_i}$、${{\boldsymbol{r}}_i}$、${{\boldsymbol{c}}_i}$、${{\boldsymbol{t}}_i}$和${{\boldsymbol{d}}_i}$分别为第i个样本的重力异常及梯度数据、水平区域、中心深度、垂直厚度和剩余密度特征，N为训练样本数，${\lambda _{\text{r}}}$、${\lambda _{\text{c}}}$、${\lambda _{\text{t}}}$和${\lambda _{\text{d}}}$为不同支路上损失项的权重.

在预测过程中，设置水平区域特征的分类阈值为0.5. 若预测值大于0.5，则认为该区域存在异常体；反之则无. 根据预测特征重构三维模型$ \hat {{m}}\left( {x,y,z} \right) $，可以表示为

(10)$ \begin{gathered} \hat {{m}}(x,y,z) = \\ \left\{ {\begin{array}{*{20}{l}} \hat {{d}}(x,y), \quad \hat {{r}}(x,y) \geqslant 0.5 \;且\; \hat {{c}}(x,y) - \\ \qquad \qquad {{{\hat {{t}}(x,y)}}/{2}} \leqslant z \leqslant \hat {{c}}(x,y)+{{{\hat {{t}}(x,y)}}/{2}}; \\ 0 , \quad 其他. \end{array}} \right. \\ \end{gathered} $

式中：$\hat {\boldsymbol{r}}$、$\hat {\boldsymbol{c}}$、$\hat {\boldsymbol{t}}$和$\hat {\boldsymbol{d}}$分别为水平区域、中心深度、垂直厚度和剩余密度的预测值.

利用Geoprobe生成包含约12 000组样本的模拟数据集. 模型所处空间大小为3 200 m×3 200 m×2 000 m，被划分为均匀的64×64×40个立方体网格单元. 异常体共有4种类型，如图6所示，尺寸参数见表1. 异常体在整个空间中以4 grid(网格)的步长游走，剩余密度的取值为{0.4，0.8，1.2，1.6} g/cm³，背景密度为0. 数据集按类以9∶1的比例随机划分为训练集和测试集. 异常体数据经预处理得到梯度联合数据及4类特征. 在输入到网络之前，所有数据进行归一化.

基于Pytorch平台，使用上述模拟数据集训练MRNet. 实验过程均使用NVIDIA GeForce RTX 2080 Ti计算单元. 设置Adam^[18]优化算法和余弦退火学习率衰减算法，初始学习率为10⁻⁵，训练轮数为80，每轮训练耗时约为3 min，损失函数下降曲线如图7所示. 图中，N_i为迭代次数.

为了测试MRNet的反演效果，从测试集中随机选取4类测试样本，分别为单棱柱体模型(样本A)、双棱柱体模型(样本B)、倾斜板状体模型(样本C) 和阶梯模型(样本D). 其中样本B包含2个异常体，分别处于不同的深度，深层异常体引起的重力异常相较于浅层异常体较小，难以直接分辨. 样本C中的倾斜异常体由3层板状体组成，所处深度较大. 样本D中的薄层异常体由5层板状体组成，其厚度较小，边界曲折，难以预测. 如图8所示为各测试样本的重力异常，使用方框标注异常体的水平轮廓.

将重力异常输入MRNet，根据预测特征重构异常体的反演模型，如图9所示为各测试样本的模拟模型和MRNet的反演模型. 计算MRNet所预测4类特征的最大绝对误差e_r、e_c、e_t、e_d，如表2所示.

使用3D U-Net和DecNet验证MRNet的优越性，设置相同的数据集和超参数训练网络. 3D U-Net每轮训练耗时约为40 min，是MRNet的13倍，验证了特征分解在加速网络训练上的有效性，DecNet每轮训练耗时与MRNet相近. 使用相同的4种测试样本进行预测，如图10所示. 3D U-Net(见图10(a)~(d))可以预测出4类样本的大致位置和形状，但出现了较大的剩余密度偏差，且在模型边缘出现较多的异常点突起. DecNet(见图10(e)~(h))相较于3D U-Net消除了大部分的异常点突起，但异常体边缘较模糊，存在剩余密度不均的现象，对于样本A出现了严重的剩余密度偏差. MRNet(见图9(e)~(h))在DecNet的基础上大幅消除了剩余密度误差，反演模型的边缘更加清晰完整.

为了定量分析3种反演网络的优劣，使用平均绝对误差MAE、全局绝对准确率E_acc和局部相对准确率R_acc衡量反演网络在测试集上的效果.

(11)$ {\text{MAE}} = \frac{1}{V}\iiint_V {\left| {{{m}}(x,y,z) - \hat {{m}}(x,y,z)} \right|} {\mathrm{d}}x{\mathrm{d}}y{\mathrm{d}}z. $

(12)$ {{{E}}_{{\text{acc}}}} = \frac{1}{V}\iiint_V {\left| {{{m}}(x,y,z) - \hat {{m}}(x,y,z)} \right| < {t_1}} {\mathrm{d}}x{\mathrm{d}}y{\mathrm{d}}z. $

(13)$ {{{{R}}_{{\text{acc}}}} = \dfrac{1}{{{V_{\mathrm{a}}}}}\displaystyle\iiint_{{V_{\mathrm{a}}}} {\max \left\{\dfrac{{{{m}}(x,y,z)}}{{\hat {{m}}(x,y,z)}},\dfrac{{\hat {{m}}(x,y,z)}}{{{{m}}(x,y,z)}}\right\} < {t_2} {\mathrm{d}}x{\mathrm{d}}y{\mathrm{d}}z}}.\\ $

式中：${\boldsymbol{m}}$和$ \hat {\boldsymbol{m}} $分别为模拟和反演模型；$V$为全局反演空间的大小，${V_{\mathrm{a}}}$为模拟异常体的体积，设置不同的阈值$ {t_1} $、$ {t_2} $. 各反演网络在测试集上的表现如表3所示. 在以上3种评价指标上，DecNet与3D U-Net有着相近的准确度，但是DecNet由于不同解码器的预测功能未能分化，难以从共享特征中学习到不同的异常体特征，准确度略低. 梯度联合和CA模块的引入使该问题得到一定的解决，MRNet相较于3D U-Net和DecNet在3种指标上均有明显的提高，特别是在R_acc上有着高达30%的提升，达到了较高的反演准确度.

由于样本D属于薄层异常体，且具有复杂轮廓，对算法反演能力的要求更高，以样本D为例进一步分析MRNet在反演异常体轮廓上的优势. 将样本D的反演模型与模拟模型作差得到误差模型，将误差模型分别沿X、Y和Z轴进行最大值投影，结果如图11所示. 图中，黑色方框为异常体轮廓. 3D U-Net反演模型的误差主要聚集在异常体的轮廓，较多的误差点包裹着异常体，使其轮廓难以分辨. DecNet的误差点有所减少，且对于异常水平区域的反演完全准确，误差点主要存在于异常体的下表面. MRNet对于水平区域的反演完全准确，在DecNet的基础上大幅减少了误差点，对异常体整体轮廓的反演更准确.

为了验证本文算法在实际反演中的效果，将MRNet应用于Vinton盐丘的重力数据反演实验中. Vinton盐丘位于美国路易斯安那州西南部，主要由大型盐核及顶部的岩盖组成. 其中岩盖主要由石灰岩和无水岩组成，平均密度约为2.75 g/cm³，剩余密度约为0.55 g/cm³. Ennen等^[19-21]的研究结果显示，盖岩南北向距离约为1 100 m，东西向距离约为1 500 m，顶部所处深度为130~160 m，底部埋深为360~460 m. 2008年，Bell Geospace公司在该区域进行了航空重力数据测量，使用2.2 g/cm³的围岩密度进行地形改正，从中选取盖岩主体所在的4 000 m×4 000 m区域进行实验. 将该区域均匀划分为64×64个网格，在深度方向上划分为40个网格，网格棱长为4 000 m/64=62.5 m. 将区域内的数据进行插值获取尺寸为64×64的数据，如图12(a)所示.

MRNet的反演模型如图12(b)所示. 将反演结果分别沿X、Y和Z轴进行投影，如图13(a)~(c)所示. 从反演结果可知，盖岩顶部埋深约为156 m，底部埋深约为468 m，在X轴方向上长约为1 250 m，在Y轴方向上长约为1 500 m，平均剩余密度约为0.48 g/cm³. 反演结果与已有成果基本一致，验证了MRNet的有效性. 此外，反演结果显示盖岩的高密度部分位于其东南部，盖岩东侧存在小型地质体，这与侯振隆等^[22]的研究结果一致. 3D U-Net和DecNet的反演结果如图14所示. 3D U-Net的反演模型边缘存在大量的密度异常点，周围出现了一些散落的异常体. DecNet虽然能够得到盖岩的大致分布，但剩余密度相较于已有资料出现了明显的偏差，未能捕获盖岩东侧的小型地质体. 总体来说，MRNet对盖岩形状、位置和剩余密度的预测更准确，能够捕获主体东侧的小型地质体.

针对基于深度学习的三维目标反演算法，加速反演网络训练过程，提出基于多特征重构的三维目标反演算法. 通过特征分解实现对三维模型的压缩，降低网络计算时的内存占用. 设计多特征重构反演网络MRNet，实现对水平区域、中心深度、垂直厚度和剩余密度4类特征的预测. 利用特征重构，得到目标的三维模型. 在MRNet的输入端，通过梯度联合增强目标的水平边界信息. 在跨层连接处，引入CA模块，使解码器自适应地关注更重要的特征，实现预测功能的分化. 在模拟数据测试集上MRNet的局部相对准确率相较于3D U-Net提升了30%以上，取得了更优的反演效果. 将MRNet应用于Vinton盐丘地区，获得盖岩的分布情况和剩余密度，且与前人的研究成果相符，验证了MRNet具有较好的泛化性能.