利用U-Net网络增强骨密度全波形反演

谢辉武 杨 艳

(武汉大学物理科学与技术学院 武汉 430072)

0 引言

全波形反演(Full waveform inversion,FWI)是基于声波方程约束的优化问题。全波形反演最初由Tarantola 提出，主要是时域反演。当然作为波动方程约束的问题，也被推广到了频域。全波形反演最初是用在地质勘探中，使用地震波来进行反演，之后也用在了工程和医学等领域。为了解决反演中的病态问题，目前最常用的是Tikhonnv正则化方法。使用更凸的目标函数可能减小反演对准确初始模型的需求，如Warner等[1]提出的自适应反演，Yang[2]提出了使用新的目标函数可以减小反演的病态问题。Bernard 等[3]尝试将全波形反演应用在骨密度测量中。但是己有的研究无法解决骨模型中数值差异大、更容易陷入局部极值的问题，反演结果中皮质骨的声速值只有2800 m/s，距离皮质骨的实际声速3500～4200 m/s还有很大的差距。

1 方法

本节将介绍如何将神经网络用于增强全波形反演。

1.1 逆时偏移成像

逆时偏移成像是一种基于波动理论的成像方法。成像的步骤为：首先通过模拟零时刻到最大时刻的发射点波场传播得到模型的发射点波场。然后将接收点接收到的波场从最大时刻逆时传播到零时刻，得到接收点波场。最后按照一定的成像条件使用两种波场进行成像。该方法最早是由Whitemore[4]在1983年的SEG 会议上提出的。本文使用的成像条件是简单且常用的非时空移相关成像法，将震源波场与接收波场进行零延时互相关就可得到偏移成像结果。

正向传播使用有限差分法数值模拟波动方程：

其中，ρ是密度，K=ρv2是弹性模量，p是波场值，s是发射点波场。

正向传播可以得到接收数据。将接收数据逆时延拓，在合适的初始模型上传播可以得到接收点波场R，将其与发射点波场执行互相关运算可以得到偏移成像。在地质成像的过程中因为模型较大，常常需要使用一些时间换空间的算法，如使用边界存储等方法减小存储量[5]。与地质成像不同，骨骼的模型大小较小，储存波场所需的空间也不是很大。在偏移成像时可以直接存储正向和反向传播的波场，所以不需要使用边界储存等用时间换空间的算法。

因为长骨模型较为简单，在进行逆时偏移时选取的初始背景模型可以是匀质模型，速度假定为软组织平均速度1500 m/s，密度为1000 kg/m3。

1.2 神经网络设计

Ronneberger 等[6]提出了一种网络结构UNet。这种网络结构由下采样部分和上采样部分构成，下采样获取信息，上采样重构图像。为了在重构时可以获取原始的位置信息，采用了跳跃连接层。U-Net 的输出可以用来进行图像分割，之前主要应用于医学图像分割，近年来也有将其应用到反演中。如Yang 等[7]提出的通过地震波形来生成地下速度模型，Zhu 等[8]使用U-Net 来提取走时。本文使用U-Net 提取逆时偏移成像的信息，重建出骨密度测量中皮质骨，骨髓和软组织的分布。

1.2.1 数据集

人体皮质骨声速为3500～4200 m/s，密度为1800 kg/m3左右。皮质骨的厚度为4 mm 左右。可以根据这些数据生成大量有着不同组织厚度的模型，并给相同的组织在合理范围内提供不同的速度与密度值，使用有限差分法来正向模拟声波的传播，将接收点相应位置的数据收集起来进行逆时互相关偏移成像。将不同模型生成的偏移成像作为神经网络的输入。本文使用1500 m/s、1000 kg/m3的背景值，使用225 种不同厚度和位置的简单平值骨骼模型，将其速度和密度赋值为不同的3 组，共得到675 种简单模型。再使用随机方法生成背景上的675 种不同形状和物性的模型，最后使用的样本一共有1350 种。因为模拟使用是波动方程的正向推理，当差分精度合适时可以得到足够精度的接收数据。超声检测装置使用线性阵探头，可以在手臂长骨的上下面各放置一个。假设探头的间距是2 mm，在61 mm×51 mm 的模型中可以有62 个模拟源。表1介绍了一些正向模拟时的参数。

表1 超声正向模拟时参数设置Table 1 Parameter setting during ultrasonic forward simulation

使用的超声源为雷克子波，它可以看作是高斯函数的二次导数。它的表达式为公式(3)。波形和频谱图如图1所示。

图1 雷克子波波形与频谱图Fig.1 Ricker wavlet wavefomr and spectrum

1.2.2 网络结构

本文采用的U-Net 网络结构如图2所示。输入是1通道的偏移成像结果，下采样过程每个单元由2个卷积层和其后的ReLU 激活层构成，如果需要也可以在每个单元的输入前加上一个参数标准化层(BN 层)和在卷积层之间添加Dropout 层。本文因为网络任务简单，训练过程也没有出现过拟合的现象，所以并没有使用Dropout 层。每个单元之后经过2×2 的最大池化层，张量尺寸减小。上采样使用反卷积来完成。在反卷积的过程中，由于输入张量的大小会出现奇数，反卷积后的张量大小与下采样过程中的大小不匹配。而进行跳跃连接需要两端大小相同，这时可以使用简单的裁剪将上采样过程中的张量大小裁剪为与下采样过程中对应深度的张量大小相同。经过网络后输出为3 通道的分类结果，最后采用交叉熵作为损失函数。将之前用模拟方法生成的数据集送入网络进行训练。该网络参数量为31, 054, 275。其中有31, 042, 499 个参数是需要训练的。该网络的训练量并不大，甚至可以在一般的CPU上完成训练。

图2 U-Net 神经网络结构Fig.2 Neural network structure diagram

该网络训练的正确率和损失变化如图3所示。可以看到随着训练次数增加，正确率有了很大的提高，可以达到99%以上。

图3 神经网络训练效果Fig.3 Neural network training results

表2展示了一些训练时的网络超参数，U-Net的卷积核滤波器基础数目为64，每一层的数目为64×2depth，depth表示网络层在U-Net中的深度。训练时一共用了100个Epoch，学习率设为1.0×10−5。

表2 U-Net 超参数设置Table 2 U-Net hyperparameters

1.3 全波形反演的目标函数

全波形反演的优化目标是让检测波场和通过预测模型模拟得到的接收波场均方差最小，因为该优化问题在数学上是病态的，为了结果的稳定需要添加上正则化，常用的有Tikhnonv 正则化[9]。本文提出可以将通过神经网络预测的模型m0加入目标函数作为约束，看作是另外一种正则化。

式(4)中，pobs和pcal分别代表观测的波场和模拟计算的波场，m是反演模型。而不同的α代表相应的正则系数，用来控制该项对结果的影响力。因为m0的数值并不一定准确，所以需要将其系数调整为合适的值。为了增强分布约束，可以在反演开始的几个迭代上使用较大的αcons，然后再将其调小，使波形匹配占据主要地位。

全波形反演可以使用各种优化方法，如梯度下降法、高斯-牛顿法等。因为反演对速度和密度的敏感性不同，实际反演中还需要一定的反演策略，比如说对密度使用更小的更新步长。本文综合考虑反演准确度和计算速度后决定使用高斯-牛顿法，其中观测波场对模型的导数可以使用伴随状态法(格林函数法)求解[10]。

2 模拟实验

本文主要验证了以下几种情况的对比。

2.1 正常数值的模型

如图4所示，使用的模型长度为6 cm，宽度为5 cm，网络间距为1 mm。皮质骨声速为4000 m/s，密度为1800 kg/m3。经过偏移成像和神经网络后的预测结果如图5所示，可以看到模型被准确地分为软组织、皮质骨、骨骼3类。可知对于这种简单模型，神经网络的预测结果非常准确。当不加约束，只使用Tikhnonv 正则化时，结果如图6(a)和图6(b)所示；而添加约束后的结果如图6(c)和图6(d)所示。根据对比可以明显看出当约束与模型相同时对反演结果有很大的提升，中间的骨髓区域结果完整。接下来讨论对于骨质疏松，声速和密度都有所下降的结果。

图4 正常的骨骼模型Fig.4 Normal bono model

图5 神经网络对图3模型的分类预测Fig.5 The neural network predict result for Fig.3

图6 正常模型两种反演方法的结果对比Fig.6 The normal model results comparison between two methods

2.2 低密度模型

如图7所示，使用的模型大小与之前相同，皮质骨声速为3500 m/s，密度为1700 kg/m3。经过偏移成像和神经网络后的分类结果与图5相同。然后执行全波形反演，最后的结果如图8所示。因为此时模型中数值的变化没有之前的大，没有加约束的结果也还可以，不过在中间的骨髓区域加入了约束的结果形状和数值都优于未加约束的结果。对比图8(a)与图8(c)可以看出速度图像的形状更好，在骨髓区域的数值也优于未加约束的结果。图8(b)和图8(c)的对比更加明显，密度的结果图形状更好。图9展示了速度和密度的皮质骨横向剖面和反演区域中间的纵向剖面结果对比。经过图片的对比可以看出约束后骨髓区域的数值更加优秀。密度区域与假设的理论模型更加接近。而且进行了约束的反演对于正则系数的健壮性有很大的提高。

图7 低密度模型Fig.7 Low density bone model

图8 低密度模型两种反演方法的结果对比Fig.8 Results comparison between two methods for low density model

图9 低密度模型两种反演方法剖面对比Fig.9 Comparison of profile results of two inversion methods for low-density model

3 结论

与骨密度检测的“金标准”双能X 射线比较来说，声学法仪器便宜，无辐射，有可能小型化后方便家庭进行提前检查。而定量声学检测法的准确度还需要实践的检验。如果使用本文的方法可以对全波形骨密度测量的准确度有一定的提高，还可以减小反演对初值的依赖。不过现有的骨密度检测的标准是骨矿物质密度，使用的是骨头对光子的吸收能力相对值。使用声学检测需要在未来进一步研究声学的相对标准或它与光学法的结果对应关系。