摘要:由于皮质骨和软组织间较大的声速差异, 采用固定声速的传统超声波束形成方法无法重建皮质骨图像, 同时皮质骨中较大的衰减也限制了信号信噪比. 为了实现皮质骨超声成像, 本文提出一种采用合成孔径超声提高成像分辨率及信噪比, 利用压缩感知计算延时参数并构建多层声速模型的成像方法. 本文结合时域有限差分仿真方法分析了理想情况下皮质骨成像结果, 并结合软组织覆盖下的离体皮质骨板样本实验, 验证相关方法的可行性. 仿真和实验结果均表明, 本文方法可用于构建多层声速模型并正确重建皮质骨图像. 本研究实现了具有三层声速模型的皮质骨超声成像, 对皮质骨超声成像发展有一定的借鉴意义, 未来将进一步探索在体实验, 以推进骨超声成像的临床应用.
关键词: 合成孔径超声/
皮质骨成像/
声速模型/
压缩感知

English Abstract

全文HTML

--> --> -->

随着人口老龄化加剧, 骨质疏松作为一种常见的代谢性骨病在日益加剧公共卫生负担, 相关骨质检测与成像方法得到了学术界的广泛关注, 成为近年来的研究热点^[1]. 目前常用的骨质诊断技术有双能X线吸收测定法(dual-energy X-ray absorptiometry, DXA)、定量CT法及定量超声法等^[2]. 在骨质评价中, DXA仍是骨质状况评价的金标准. 由于超声评价方法具有无辐射和便携低廉等特点, 更适于临床普查与病患监测, 近年来得到较快发展^[3].
研究表明, 长骨皮质骨能够支持超声导波的传播, 基于轴向传播法所发展起来的超声导波长骨检测技术, 能提供长骨皮质骨的材料特性和内部结构等信息, 已被用于长骨皮质骨评价及骨质疏松诊断^[4,5]. 超声导波可用于长骨中频散和衰减的测量^[6,7], 导波的模式转换规律可用于长骨骨折评估^[8], 其相速度的变化可用于长骨疲劳评估^[9]. 导波信号的Randon变换可用于长骨皮质骨厚度估计^[10]. 临床测试表明, 点医疗超声(point-of-care ultrasound, PoCU)检测长骨骨折的灵敏度可达到64.7%—100%^[11]. 刘洋等^[12]利用超声非线性和时间反转技术对长骨骨裂成像进行了仿真研究, 并分析了裂纹深度对成像结果的影响. 然而, 由于骨与软组织间显著的声速差异, 固定声速的传统超声波束形成用于皮质骨成像仍存在困难. Li等^[13]结合已知声速模型对皮质骨进行了图像重建. Renaud等^[14]采用射线追踪(ray-tracing)法寻找皮质骨最佳成像结果, 并对皮质骨参数进行了估计.
合成孔径技术最早应用于雷达系统, 20世纪70年代起应用于超声成像^[15]. 相比传统定点聚焦, 合成孔径采用动态聚焦, 可提升整体探测区域的成像质量. 已有大量研究将合成孔径超声成像技术用于血管内成像^[16]、无创颈动脉弹性成像^[17]、肝脏病变成像等软组织检测^[18]. 然而骨成像相关研究仍有待完善, 近年来Renaud等^[14]的研究结果表明了合成孔径超声骨成像的可行性. 此外, 考虑皮质骨与软组织间显著的声速差异, 建立正确可靠的声速模型是皮质骨超声成像的关键. 针对多层介质成像建立声速模型最早应用于地球物理学^[19], 近年来广泛应用于无损检测领域^[20]. 相位迁移法(phase shift migration, PSM)是一种可用于多层声速模型的频域重建算法, 最早用于反射地震学中地球内部结构成像^[21]. Olofsson^[22]将PSM用于无损检测领域多层介质成像.
皮质骨可简化为软组织-皮质骨-骨髓三层介质模型, 利用压缩感知方法提取皮质骨上下表面回波, 在无需皮质骨厚度先验信息的情况下建立声速模型. 本文采用合成孔径方法测量超声信号, 通过压缩感知计算延时参数进而建立声速模型, 最后在合成孔径数据和声速模型的基础上, 利用PSM算法对皮质骨在频域进行图像重建. 通过时域有限差分(finite-difference time-domain, FDTD)仿真和离体骨板实验验证上述方法成像的可行性和准确性.

方法流程如图1所示, 主要包括合成孔径测量超声信号, 压缩感知计算延时参数进而估计皮质骨厚度建立声速模型, 及PSM算法重建图像. 下文中将一一展开介绍.
图 1 方法流程图
Figure1. Flow chart of the proposed method.

2

-->

2.1.基于合成孔径方法的超声信号测量

传统波束形成多采用定点聚焦, 即发射和接收时只对一点进行聚焦. 定点聚焦仅在焦点处有较高分辨率, 而焦点以外区域分辨率较低. 此外, 定点聚焦导致成像的方位分辨率受限于换能器参数^[23]. 通常采用大孔径换能器或提高换能器工作频率以得到更高的方位分辨率. 然而过大孔径换能器缺少实用性, 且换能器工作频率与最大可探测深度成反比, 无法同时满足方位分辨率高和探测范围广的要求. 不同于定点聚焦, 图2(a)所示的合成孔径超声通过小孔径换能器的水平移动依次发射接收数据, 研究表明合成孔径超声可有效解决定点聚焦存在的分辨率问题^[23].
图 2 合成孔径超声原理　(a)合成孔径聚焦; (b)发射合成孔径
Figure2. Principle of synthetic aperture ultrasound: (a) Synthetic aperture focusing technique; (b) synthetic transmit aperture.

有关成像方位分辨率的计算如下^[23]:
定点聚焦

${\beta _{\rm{f}}} = {\rm{0}}{\rm{.84}}\frac{\lambda }{d},\;\frac{\lambda }{d} \ll 1,$

$ {\rho _{\rm{a}}} = {\beta _{\rm{f}}}\cdot{R} = {\rm{0}}{\rm{.84}}\frac{\lambda }{d}R; $

合成孔径

${L_{{\rm{sm}}}} = 0.84\frac{\lambda }{d}{z_0},$

$ {\beta _{\rm{s}}} = 0.42\frac{\lambda }{{{L_{{\rm{sm}}}}}} = \frac{d}{{2{z_0}}} ,$

$ {\rho _{\rm{s}}} = {\beta _{\rm{s}}} \cdot {z_0} = \frac{d}{2}. $

(1)和(2)式中, ${\beta _{\rm{f}}}$为定点聚焦的半功率波束角, $\lambda $为工作波长, d为换能器孔径直径, ${\rho _{\rm{a}}}$为方位分辨率, R为成像目标距离超声换能器的斜距. (3)—(5)式中, ${\beta _{\rm{s}}}$为合成孔径的半功率波束角, ${L_{{\rm{sm}}}}$为最大综合长度, ${\rho _{\rm{s}}}$为合成孔径方位分辨率, ${z_0}$为成像目标距离超声换能器的垂直距离. 对比(2)和(5)式可知, 合成孔径换能器尺寸越小方位分辨率越高, 且方位分辨率与距离无关, 保持任意探测区域分辨率的一致性.
仿真与实验均采用图2(b)所示发射合成孔径超声对信号进行测量. 线性阵列共有N个阵元, 每个阵元依次发射, 由全部阵元接收. 发射合成孔径相比于合成孔径聚焦, 可以得到更加完整的数据信息, 实现发射和接收的全动态聚焦, 提高成像质量和帧率^[24]. 发射合成孔径每次发射接收得到1组数据, 通过N次循环可以得到完备集数据, 用于后续图像重建.
2 -->

2.2.压缩感知延时参数计算与声速模型建立

压缩感知广泛应用于具有稀疏变换域的信号^[25]. 假设原信号$x \in {R^N}$为N维实信号, 构造正交基${\psi }= [{\psi _1},{\psi _2}...{\psi _N}] \in {R^{N \times N}}$, 则x在正交基$\psi $下的表示为^[25]

$x = \sum\limits_{i = 1}^N {{\theta _i}} {\psi _i}, $

(6)式的矩阵形式如下:

${x} = {\psi \theta }, $

其中${\theta }= {[{\theta _1},{\theta _2},\cdots,{\theta _N}]^{\rm{T}}}$为x在$ {\psi} $上的系数向量. 假设系数向量${\theta }$中只有K个非零量, 即x可由K个基的线性组合进行重构, 则x是K稀疏的.
考虑皮质骨表面回波信号波包的时域稀疏性, 采用压缩感知计算延时参数, 并建立声速模型. 回波信号可视作发射脉冲经过延时和幅度变化后的线性叠加, 因此设计不同延时发射脉冲作为压缩感知基底, 延时间隔设置为采样间隔. PSM算法要求数据是0延时形式, 首先采用压缩感知提取各个接收通道间的延时, 将完备集数据调整为0延时形式, 以便直接应用PSM算法. 进而提取骨与软组织交界处的两次反射回波所对应基底, 根据基底间延时和已知骨声速估计骨板厚度, 从而建立声速模型.
2 -->

2.3.多层介质PSM原理

PSM是基于爆炸反射模型(exploding reflector model)的一种频域算法. 在爆炸反射模型中, 假设声速为v, 当发射阵元和接收阵元位置固定时, 声波从发射阵元至探测点的路径及探测点至接收阵元的路径一致, 因此反射波可视为由反射表面的一系列点源在$t = 0$时刻产生, 此时声速为${\hat v} = \dfrac{v}{2}$.
首先建立图3(a)所示x-z坐标下的单层介质模型, 根据波动方程有^[22]
图 3 PSM模型　(a)单层介质模型; (b)多层介质模型
Figure3. PSM model: (a) Single layer model; (b) multi-layer model.

$\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {z^2}}} - \frac{1}{{{{\hat v}^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}}} \right)p(t,x,z) = 0, $

其中$p(t,x,z)$表示$(x,z)$处t时刻的声场. (8)式的解为^[22]

$p(t,x,z) = P{{\rm{e}}^{{\rm{i}}({k_x}x + {k_z}z - \varpi t)}}, $

其中, 复数形式P表示幅度, $\varpi $表示角频率, ${k_x}$和${k_z}$分别表示x方向和z方向的波数. 根据(8)和(9)式可以得到^[22]

$\frac{{{\varpi ^2}}}{{{{\hat v}^2}}} = k_{_x}^2 + k_z^2, $

考虑回波传播方向为–z, 解得${k_z}$为^[22]

${k_z}(\varpi ,{k_x}) = - {\rm{sgn}} (\varpi )\sqrt {{{\frac{\varpi }{{{{\hat v}^2}}}}^2} - k_{_x}^2} , $

其中${\rm{sgn}} (\varpi )$为符号函数. 令$P(\varpi ,{k_x})$为特定$(\varpi ,{k_x})$下的幅度, 并设$P(\varpi ,{k_x},z) \equiv $$ P(\varpi ,{k_x}) {{\rm{e}}^{{\rm{i}}{k_z}(\varpi ,{k_x})z}} $. 根据(9)式, 任意声场可表示为^[22]

$p(t,x,z) = \iint {P(\varpi ,{k_x},z){{\rm{e}}^{{\rm{i}}({k_x}x - \varpi t)}}{\rm{d}}{k_x}{\rm{d}}\varpi },$

观察可知(12)式为傅里叶反变换形式, 相应傅里叶变换形式为^[22]

$P(\varpi ,{k_x},z) = \frac{1}{{4{{\text{π}}^2}}}\iint {p(t,x,z){{\rm{e}}^{ - {\rm{i}}({k_x}x - \varpi t)}}{\rm{d}}x{\rm{d}}t}, $

令$P(\varpi ,{k_x},Z)$表示当前深度Z处声场的傅里叶变换形式, 则任意深度z处声场$P(\varpi ,{k_x},z) = $$P(\varpi,{k_x},Z){{\rm{e}}^{{\rm{i}}{k_z}(z - Z)}} $, 设$\Delta z = z - Z$, (12)式可重新写作^[22]:

$\begin{split} & p(t,x,\Delta z + Z)\\ = & \iint {P(\varpi ,{k_x},Z){{\rm{e}}^{{\rm{i}}{k_z}\Delta z}}{{\rm{e}}^{{\rm{i}}({k_x}x - \varpi t)}}{\rm{d}}{k_x}{\rm{d}}\varpi },\end{split}$

根据爆炸反射模型, 爆炸瞬间($t = 0$)声场最为集中, 聚焦程度最大, 因此要求数据为0延时形式.
在单层介质模型基础上设计多层介质模型. 如图3(b)所示, ${Z_0}$为线阵放置位置, ${d_l}$, ${v_l}$, ${Z_l}$分别表示第l ($l = 1,2,3$)层介质的厚度、声速和交界处位置, ${Z_l} = {Z_0} + \sum\limits_{n = 1}^l {{d_n}} $. 在多层介质模型中(11)式改写为^[22]

${k_{zl}}(\varpi ,k_x) = - {\rm{sgn}} (\varpi )\sqrt {{{\frac{\varpi }{{{{\hat v}_l}^2}}}^2} - k_{_x}^2} . $

对于各层介质内$(\Delta z < {d_l})$的声场,

$P(\varpi ,{k_x},\Delta z + {Z_l}) = P(\varpi ,{k_x},{Z_l}){{\rm{e}}^{{\rm{i}}{k_{zl}}\Delta z}}.$

对于交界面处的声场,

$P(\varpi ,{k_x},{Z_l}^ - ) \propto P(\varpi ,{k_x},{Z_l}^ + ).$

因此PSM可对多层介质模型任意深度z成像.
2 -->

2.4.FDTD仿真

目前FDTD方法广泛应用于电磁学、光学、声学等领域^[26]. 仿真中建立二维皮质骨模型, 采用FDTD数值仿真方式模拟超声波动. FDTD方法将计算区间离散化为众多微分单元, 在时间和空间域上进行差分化, 进而得到超声波动控制方程的数值解^[27]. 超声波动控制方程为

$\rho \frac{{{\partial ^2}{ w}}}{{\partial {t^2}}} = (\lambda + \mu )\nabla (\nabla \cdot { w}) + \mu {\nabla ^2}{ w}, $

其中$\rho $为材料密度; w为某点处位移矢量; $\lambda $, $\mu $分别为第一、第二拉梅常数. FDTD仿真中, 皮质骨模型设置为三层介质模型, 上层为3 mm厚软组织, 中间为3.4 mm厚皮质骨, 下层为3.6 mm厚骨髓. 软组织及骨髓的密度设置为1000 kg/m³, 第一拉梅常数设置为2.241 GPa, 第二拉梅常数设置为0 GPa, 声速设置为1500 m/s^[28]. 皮质骨密度设置为1850 kg/m³, 第一拉梅常数设置为9.306 GPa, 第二拉梅常数设置为3.127 GPa, 声速设置为2900 m/s^[29]. 仿真采用128阵元线性阵列, 线阵放置于软组织表面, 阵元间隔0.3 mm, 中心频率6.25 MHz, 网格尺寸设置为0.02 mm × 0.02 mm, 迭代频率250 MHz. 为减少多次回波的干扰, 在模型四周设置4 mm厚完全匹配层(perfectly matched layer). 发射脉冲为高斯包络调制的正弦波. 采用发射合成孔径方法测量信号.
2 -->

2.5.实验设置

实验中选用3.4 mm厚牛胫骨骨板作为实验材料. 测量声速为2900—3100 m/s, 与仿真参数基本一致. 为更好地模拟在体实验, 用3%琼脂制成软组织, 覆盖于骨板四周. 实验装置见图4, 采用Verasonics系统(Vantage 128 or 256, Verasonics Inc, WA, USA)进行信号测量, 线阵紧贴软组织表面, 探头型号为L11-4v, 包含128个阵元, 阵元间距0.3 mm, 中心频率6.25 MHz, 采样频率25 MHz. 发射脉冲为高斯包络调制的正弦波. 为减少成像深度影响, 实验中采用线性时间增益补偿. Verasonics系统控制阵元依次发射脉冲, 每次发射所有128个阵元进行接收, 数据经总线传输到计算机.
图 4 实验装置示意图
Figure4. Experiment setup.

2

-->

3.1.仿真结果

图5为皮质骨FDTD仿真结果, 其中${A_{\rm{N}}}$为归一化幅度. 图5(a)为仿真采用的发射脉冲经过不同延时组成的压缩感知基底. 图5(b)为FDTD仿真的单通道接收信号, 在4.6和7.0 μs附近可观察到两次明显的反射回波, 分别对应软组织和骨板上表面的分界面以及骨板下表面和软组织的分界面. 因超声存在衰减, 第二回波幅度明显小于第一回波. 图5(c)为第一个阵元发射时, 前30个通道的接收信号, 因回波信号到达各个阵元的路径不同, 通道间波形存在延时. 压缩感知计算各通道间接收延时参数, 并对信号进行延时调整, 结果如图5(d)所示. 与图5(c)对比可得, 调整后信号的第一回波延时已对齐, 可进行后续图像重建.
图 5 仿真结果　(a)压缩感知前30个基底; (b)单通道接收信号; (c)单次发射全部通道接收信号; (d)压缩感知调整的单次发射全部通道接收信号
Figure5. Simulated results: (a) The first 30 bases of compressed sensing; (b) received signal of single element; (c) received signals of all elements; (d) compressed sensing based temporally adjusted received signals of all elements.

图6为采用PSM算法对FDTD仿真数据进行图像重建的结果. 未建立声速模型的重建结果如图6(a)所示, 虽可获得清晰的上下骨板界面, 但因声速未经矫正, 骨板厚度存在误差, 平均骨板厚度约为1.8 mm, 相对误差47.1%. 经压缩感知建立的多层声速模型调整后的重建结果如图6(b)所示, 平均骨板厚度约为3.4 mm, 相对误差为0%, 使用多层声速模型可实现骨板厚度的正确成像.
图 6 仿真重建结果　(a)未建立声速模型的仿真重建结果; (b)建立声速模型的仿真重建结果
Figure6. Simulated reconstructed results: (a) Simulated reconstructed result without velocity model; (b) simulated reconstructed result with velocity model.

为验证压缩感知建立声速模型的准确性, 采用FDTD方法对3—5 mm不同厚度的皮质骨进行多次仿真. 压缩感知估计的皮质骨厚度及误差如表1所列.

参量	数值
真实厚度/mm	3.0	3.1	3.2	3.3	3.4	3.5	3.6	3.7	3.8	3.9	4.0	4.1	4.2	4.3	4.4	4.5	4.6	4.7	4.8	4.9
估计厚度/mm	2.7	3.2	3.3	3.0	3.4	3.9	3.9	3.3	3.5	4.0	4.1	4.2	4.2	4.0	4.5	4.5	4.3	4.6	4.6	4.7
相对误差/%	10	3.2	3.1	9.1	0	11	8.3	11	7.9	2.6	2.5	2.4	0	6.9	2.3	0	6.5	2.1	4.2	4.1

表1皮质骨厚度估计及误差
Table1.Estimation and relative error of cortical bone thickness.

压缩感知估计皮质骨厚度的平均相对误差为4.9%, 方差为13.5%.
2 -->

3.2.实验结果

图7为牛胫骨骨板的实验结果. 图7(a)为实验采用的发射脉冲经过不同延时组成的压缩感知基底. 图7(b)为实验的单通道接收信号, 在16.3和18.7 μs附近可观察到两次明显的反射回波, 分别对应软组织和骨板上表面的分界面以及骨板下表面和软组织的分界面. 因超声存在衰减, 第二回波的幅度明显小于第一回波. 图7(c)为第一个阵元发射时, 前30个通道的接收信号, 因回波信号到达各个阵元的路径不同, 通道间波形存在延时. 由压缩感知计算各通道间接收延时参数, 并对相应信号进行延时调整, 结果如图7(d)所示. 与图7(c)对比可得, 调整后信号的第一回波延时已对齐, 从而可进行后续图像重建.
图 7 实验结果　(a)压缩感知前30个基底; (b)单通道接收信号; (c)单次发射全部通道接收信号; (d)压缩感知调整的单次发射全部通道接收信号
Figure7. Experiment results: (a) The first 30 bases of compressed sensing; (b) received signal of single element; (c) received signals of all elements; (d) compressed sensing based temporally adjusted received signals of all elements.

图8为采用PSM算法对实验数据进行图像重建的结果, 图8(a)与图8(b)分别对应未建立声速模型的重建图像以及基于多层声速模型的重建图像. 如图8(a)所示, 虽可看到清晰的上下骨板界面, 但由于没有进行声速矫正, 骨板厚度存在误差, 平均骨板厚度约为1.8 mm, 相对误差47.1%, 形态重建不准确. 而多层声速模型的应用实现了准确的板厚成像. 对同一骨板进行重复实验, 估计得到平均骨板厚度约为3.5 mm, 平均相对误差为3.6%, 方差为5.4%. 可见声速模型可实现较为准确的骨板成像.
图 8 实验重建结果　(a)未建立声速模型的实验重建结果; (b)建立声速模型的实验重建结果
Figure8. Experiment reconstructed results: (a) Experiment reconstructed result without velocity model; (b) experiment reconstructed result with velocity model.

由于皮质骨和软组织间显著的声速差异, 传统超声波束形成无法应用于骨成像, 建立正确声速模型是解决声速差异问题的关键. 目前仿真及实验中声速模型的建立需要事先测量真实骨板尺寸, 声速模型的实际应用受到限制. 因此本文提出了基于压缩感知建立声速模型的方法, 通过压缩感知计算延时参数从而在未测量骨板尺寸的情况下建立较为准确的声速模型.
图9比较了压缩感知和Hilbert变换提取延时的结果. 图9(a)中黑线表示原始信号, 原始信号由两个不同延时的发射脉冲叠加组成, 红线为信噪比20 dB的带噪信号. 图9(b)为压缩感知提取结果, 其中峰值对应发射脉冲起始点. 图9(c)为Hilbert变换提取结果, 其中峰值对应发射脉冲包络顶点. 比较图9(b)和图9(c)可知, 压缩感知的提取结果受噪声影响较小, 延时提取准确, 而Hilbert变换在有噪声条件下提取的信号峰值相比无噪声条件下发生了显著的偏移. 压缩感知是一个求解最优解的过程, 在计算中允许误差项存在, 因此压缩感知法比Hilbert变换更适用于有噪声条件下的延时估计. 且压缩感知提取结果对应脉冲起始点, 可直接用于波形到达时间的计算. 皮质骨仿真和实验结果均验证了压缩感知提取延时的准确性, 厚度估计的相对误差分别为4.9%和3.6%. 如图6(b)和图8(b)所示, 建立声速模型的仿真重建结果和实验重建结果均可以看到清晰的皮质骨上下界面, 且相比图6(a)和图8(a)未建立声速模型的结果, 皮质骨重建厚度准确. 与仿真重建结果比较, 实验重建结果存在明显伪影, 这是实验过程中存在噪声干扰导致的.
图 9 压缩感知与Hilbert变换比较　(a)原始信号及带噪信号; (b)针对原始信号和带噪信号的压缩感知结果; (c)针对原始信号和带噪信号的Hilbert变换结果
Figure9. Comparison of compressed sensing and Hilbert transform: (a) Origin signal and noisy signal; (b) result of compressed sensing; (c) result of Hilbert transform.

在未知实际皮质骨尺寸的情况下, 采用压缩感知方法估计模型各层延时, 结合已知声速可计算模型各层厚度, 正确建立声速模型. 目前本文方法只适用于多层规则介质模型. 射线追踪法可用于多层非规则介质的时域重建. Qin等^[30]研究了多层非规则介质的频域重建. 后续研究将围绕非规则骨板成像展开, 在现有方法上进行改进, 解决水平方向声速变化问题.

本文基于合成孔径法采集超声信号, 通过压缩感知计算延时参数建立声速模型, 采用PSM算法进行图像重建, 实现皮质骨超声成像. 仿真和实验结果表明, 本文采用的压缩感知方法可准确估计皮质骨厚度, 从而建立多层声速模型. 重建图像可观察到清晰的皮质骨上下界面, 且形态正确. 本文方法相比传统超声波束形成, 解决了固定声速重建图像的挑战, 对于皮质骨超声成像的发展有一定借鉴意义. 相比于Hilbert变换, 压缩感知方法更适于计算带噪实验信号的延时参数. 本文方法目前适于规则骨板成像, 考虑到骨的多种形态及复杂结构, 之后的工作应当建立更为实际的骨模型, 寻找更加通用的成像方法. 同时探究在体实验的可行性.