基于超声背散射的骨质评价实验研究

2021-10-27丑幸幸

陕西科技大学学报 2021年5期

汤伟,黄璜,丑幸幸

(陕西科技大学电气与控制工程学院,陕西西安 710021)

0 引言

近年来,人口老龄化发展形势严峻,骨质疏松症[1]发病率逐年升高,现已成为人们生活中常见的代谢性疾病之一.人体骨骼在结构上分为皮质骨和松质骨,而松质骨由大量的骨小梁交错组织而成.骨质疏松症主要表现为松质骨退化,骨密度下降,骨脆性增加[2].根据中国卫生健康委员会的调查结果显示,我国当下存在骨密度检测率低、低骨量人群年轻化、骨质疏松症认知率低等问题[3].临床对于骨质疏松症的建议是预防大于治疗,因此,及时的骨质评价对预防骨质疏松症十分必要.

目前基于射线的骨质检测手段,如DXA、X线、定量CT 等辐射大,不适宜人群定期使用;而基于定量超声的检测方法无辐射、可靠性高,更适用于社区普查、孕妇复检等情景.基于定量超声对松质骨及其骨质疏松的评价方法主要分为透射法和背散射法[4],其中透射法将信号在骨中的衰减和声速作为骨质评价的主要参数,但是很少反映骨微结构信息,如骨小梁厚度、间距等,而这些微结构信息已被证实与骨质疏松有直接关联[5].而背散射法通过骨小梁单元等骨微结构对超声信号的散射作用,使接收的背散射信号包含了丰富的骨密度和骨微结构信息.但是背散射法还需要建立更真实的骨模型和提取骨微结构信息更准确的算法[6].

在背散射信号分析算法发展中,多元线性回归模型没有考虑超声参数与骨参数之间的非线性联系[7],骨质评价准确度低;机器学习算法需要手动提取特征信息,在处理大量数据时操作繁杂[8];而CNN 通过最小化预测误差得到超声参数与骨参数之间的非线性关系[9],可以自主进行特征学习,相比其他算法,评价结果更为准确.

本文基于超声背散射方法,首先对骨样本进行μ-CT 断层扫描,获得大量骨样本原始图像,对原始数据图像进行处理,导入Mimics中重建三维松质骨模型;然后将重建好的不同结构的骨模型导入Wave3000 Plus软件中,模拟超声在骨模型中的传播过程,求解波动方程,获得大量背散射信号;最后将获得的背散射信号分为多组数据集输入至改进的CNN 网络中训练测试,对骨样本进行预测评价.

1 实验部分

1.1 松质骨模型建立

在骨骼结构中,松质骨虽然占骨量的20%,但构成了80%的骨表面.现阶段对松质骨模型的研究,单圆柱模型[10]和多圆柱模型[11]与真实骨骼存在较大误差,如表1所示.使用计算机仿真技术对骨科领域中的相关问题进行分析研究,避免了传统离体和在体试验方法成本高、周期长的问题.Mimics能够对二维医学影像信息进行数字化三维重建,通过骨样本CT 图像数据三维重建出的松质骨模型更近似真实骨骼.

表1 不同骨模型对比

本文对松质骨的重建过程包括原始数据处理、数据导入、轮廓提取、区域增长、模型建立.

(1)原始数据处理

本文原始图像数据来自骨样本断层扫描,扫描仪型号Micro-CT Skyscan1076,仪器工作参数在电源75 KV,电流130μA,功率10 W,松质骨样本取自真实骨骼,在去除皮质骨等预处理之后,根据骨的承重方向切割成表面平整的立方体结构(约为20 mm*20 mm*20 mm).

对于原始数据图像的预处理操作过程如下:首先μ-CT 对骨样本扫描得到松质骨CT 灰度图像,建立松质骨二值图像.然后选取一幅松质骨灰度图像作为原始图像,从任意方向随机对骨小梁表面进行腐蚀,将松质骨二值模型中的l(骨小梁)替换为0(骨髓),模拟骨流失过程,松质骨的孔隙度不断增大.最后得到不同孔隙度结构的松质骨图像模型,如图1所示.

图1 图像预处理

(2)数据导入

Mimics支持DICOM、TIFF、BMP 等多种图像格式,本文导入处理后的图像数据为BMP 格式.进入软件界面后点击File,打开New Project,导入处理后的CT 断层扫描图像,层距为1 mm,在界面窗口中调整影像所示方向,其中冠状位、矢状位、横断位三个位置界面显示的大写字母分别是:A、P和L、R 以及T、B,分别代表重建模型的前面、后面、左侧、右侧、上部、下部,可以通过上述不同视图快速定位.

(3)轮廓提取

在操作窗口选择阈值工具,设置骨组织CT值,阈值单位为HU,最小值是226 HU.在设置阈值时,阈值既不能太低,也不能太高.阈值过低会受到软组织干扰,过高则会造成骨组织信息丢失,本文设置阈值为255 HU 时,分割效果最佳,如图2所示.

图2 轮廓分割

(4)区域增长

选择区域增长按钮,点击感兴趣区域,具有近似灰度值的骨骼周围组织将被选中挑选除去,如图3所示.

图3 感兴趣区域选取

(5)模型建立

选择合适的蒙板,点击Calculate重建按钮,操作界面生成三维立体骨骼模型,可以通过平移、缩放、旋转,直观化展现重建的松质骨模型结构,如图4所示.本文根据正常人群和骨质疏松人群骨骼孔隙度的大小,重建了正常骨模型和异常骨质疏松模型共四组,并导出为STL格式文件.

图4 骨模型三维重建

1.2 超声背散射仿真实验

Wave3000 Plus是一种超声波仿真软件,适用于无损测试材料评估和超声医学等研究领域.该软件可以从STL文件中导入任何形状的3D 对象,釆用时域差分方法[12](FDTD)求解完全的粘弹性波动方程,求解快速,一次模型求解可以得到全场波的传播,有利于大量数据的采集.FDTD通过对超声背散射信号在松质骨中的传播进行数值仿真,在考虑粘弹性材料对声波吸收的基础上,定义吸收边界条件,模拟超声在重建骨模型中的传播过程.

其中,波动方程形式如式(1)所示:

式(1)中:ρ代表材料密度,单位kg/m2;λ表示第一个Lame波常数,单位N/m2;μ表示第二个Lame波常数,单位N/m2;η表示剪切粘度,φ表示堆积粘度,单位N·s/m2;▽表示梯度算子;▽·表示散度算子;∂表示偏微分算子;t表示时间,单位s;w可用式(2)计算得到.

仿真实验步骤为:首先导入STL 格式三维重建好的骨骼模型,3D 模型边界按照空间坐标分为Xmin、Xmax、Ymin、Ymax、Zmin、Zmax 六个平面;模型导入后,对模型材料进行定义.

选取库中材料骨骼和血液,骨模型样本被假定浸在血液中,并假定血液填充所有孔隙空间,根据材料库中材料固定参数如表2所示.

表2 骨模型材料参数

实验边界条件设置为Xmin 平面为吸收层,Ymin、Ymax、Zmin、Zmax 平面为完全匹配吸收层,吸收层的作用是吸收传播过程中的干扰信号;然后定义发射接收源,背散射法主要通过一个超声探头对信号进行发射接收,实验在Xmax平面设置一个1 MHz的正弦脉冲信号源,设置工作条件为Pluse echo方式.最后设置工作参数时间为30μs,步长比例为0.8,分辨波长为0.5 mm,周期为5,仿真结果如图5所示.从仿真结果可以看出,波面被骨小梁严重散射.

图5 背散射信号

1.3 背散射信号处理

图6是一个背散射信号单周期时域波形示例.在信号起点之后,可以观察到一个明显的强反射回波,这是因为超声在传播过程中遇到液体表面产生较大声阻抗差异时出现的强反射,因此在获得背散射信号后,需要对获得的原始信号进行处理.刘成成等[13,14]在选取背散射信号的有效区间时,取了多组T1、T2的值,计算了不同背散射信号有效区间对应的背散射参数,通过分析背散射参数与骨密度之间的相关性,提出了背散射信号提取标准.

图6 背散射信号提取示例

该提取标准如下:从信号的起点位置去除时域部分时间长度为T1的波段,以去除信号的末端做为新起点,重新选取一段时间长度为T2的信号作为有效信号区间(ROI),T2长度包含的信号区间则是声波在松质骨模型中传播的有效信号.第五强强等[15]在该提取标准下取了相关性较高的背散射有效信号区间作为机器学习模型的输入特征,骨质评价实验的准确率达到了82.86%.

本文在该提取标准下,对获得的3 000组背散射信号进行了相同的预处理.为了适应后续CNN的输入,需将预处理后的背散射信号数据进一步进行格式转换,使背散射信号数据转变为二维数字矩阵,最后将转换的数据切分成64*64尺寸大小的图片,并制作成数据集.其次,将数据集中1 000组作为训练集,剩余2 000组划分为四组测试集,训练集中包含有随机数量正常骨模型[16]下仿真得到的背散射有效信号.

1.4 背散射信号分析

(1)背散射信号分析方法

参考人群中骨质疏松的临床判断标准,通过背散射信号分析进行骨质评价时,需要提取骨微结构信息更准确的算法.然而传统的多元线性回归模型并未考虑超声参数与骨参数之间的非线性联系,有监督学习的超声背散射方法和ANN[17]等机器学习算法需要手动提取特征,操作繁杂.使用CNN可以将超声参数映射到高维空间中,通过最小化预测误差得到超声参数与骨参数之间的非线性关系,在骨质评价中优于传统算法.

ResNet作为经典的深度CNN 模型,本身具有很好的非线性表达能力.但是经过预处理和格式转换后的背散射信号数据是二维的数字矩阵,只能从时间和空间维度提取一些特征信息来识别,会导致分类效果变差.因此针对使用数据的特性,对ResNet18模型进行改进,提出一种ResNet16的改进CNN 模型.刘银萍等[18]已经验证Res Net模型并不是所有的层都是必要的,较少的层结构更有利于模型识别分类.在数据集数量不大的情况下,不需要多层次的残差网络,否则易出现过拟合等现象[19].

(2)改进网络结构

ResNet18的模型结构包含了1 个卷积层,4个残差块结构,每个残差块包括4个卷积层,以及1个全连接层,通过平均池化与全连接层相连,池化的结果会使得特征减少和模型参数减少[20].平均池化是对邻域内特征点进行平均运算,而最大池化是对邻域内特征点取最大值.在CNN 模型中,特征提取的误差主要来自两个方面:(1)邻域大小受限造成的估计值方差增大;(2)卷积层参数误差造成估计均值的偏移.因此本文在扩宽网络结构时,使用平均池化和最大池化分别将深度特征信息进行融合,有利于提高模型原始数据的特征信息利用率.若去除其余卷积层,虽然模型的参数会得到减少,但是特征信息却不能得到有效融合,会使特征提取产生误差[21].

本文改进的网络结构参数与未改进的Res-Net18和具有相同层数的VGG16结构参数如表3所示.其中卷积层conv1的主要作用是通过卷积运算减少输入数据大小,残差块conv2_x、conv3_x、conv4_x、conv5_x 由多个卷积层组成,其作用是对缩减后的数据进行特征提取.本文首先用背散射信号数据集从头开始训练;更改输入数据的大小,减少训练时间;缩减ResNet18 网络的深度,降低冗余参数,去除一个残差结构conv5_x,分解为两个并列的卷积层,拓宽整体模型的宽度,所用卷积核都是3*3大小,将conv5_x卷积层分成两个并列的卷积层有两个好处:首先拓宽CNN 模型的宽度有利于在提取深度特征信息时提高信息的利用率,进而提高模型的准确度;其次分解过后,使冗余参数得到了减少,降低了训练时间.分解之后的卷积层用最大池化和平均池化可以将得到的特征信息进行深度融合,使CNN 模型提取的特征信息更加准确[22].

表3 CNN 模型参数

改进CNN 模型结构如图7 所示.由于CNN模型的训练对实验配置要求较高,本文实验的训练环境采用的是Windows10系统,使用显卡的型号为RTX2060,其中计算机最大内存容量为16 GB,计算机CPU 型号为AMD R7-4800H,为八核十六线程处理器.模型的训练和编程软件使用Python3.7版本和Anaconda3,环境后端框架为Tensor Flow.

图7 改进网络结构

2 仿真结果

2.1 训练结果

最后将实现的16层ResNet16在数据集上训练测试,并与未改进的Res Net18以及具有相同层数的VGG16进行简单的性能比较.如图8 所示,改进模型训练效果的误差率对比未改进的Res-Net18和VGG16有明显降低.

图8 不同模型训练对比

由于本实验是为了对骨质状况进行分类,所以使用的评价指标主要为准确率(Accuracy)和特异性(Specificity),其中准确率可以由公式(3)和公式(4)来表示:

其中,TP(True Positive)代表真正例,即被正确预测为正例的个数;FN代表假负例,即被错误预测为负例的个数;FP代表假正例,即被错预测为负例的个数;TN代表真负例,即被正确预测为负例的个数.

2.2 测试结果

本文使用改进的CNN 模型对测试集进行测试,在迭代500次训练后,改进CNN 模型收敛效果趋于稳定.

将重建的不同骨模型下获得的四组数据集作为输入,对CNN 模型进行测试,实验准确率达到了97.3%,其中损失函数变化和准确率变化如图9所示.本文改进的ResNet16模型的预测分类准确率对比文献15 中有监督学习算法和文献17 中ANN 算法分类准确率结果如表4所示.

图9 改进模型准确率和损失函数

表4 不同算法性能对比

通过不同算法之间的对比,实验结果表明ResNet18中的层可以进行改进,而相同层数结构下,ResNet模型基于跳层连接思想比VGG16 网络和ANN 等算法能更好的表征特征信息,有利于深层特征的提取,提升模型分类准确率.本文提出的改进CNN 算法在提取骨微结构信息上鲁棒性更好,准确度更高,对骨样本有更好的区分能力.