计算生物医学超声学

-2章着重介绍了生物医学超声非线性声场和热场的经典理论与计算；第3-5章针对不同介质中的单微泡和微泡群，通过建模仿真详细讨论了微泡动力学及声散射问题；第6章详细讨论了波束合成的经典和前沿计算问题；第7-8章针对弹性成像介绍了组织弹性模量重建和分数阶黏弹性参量成像；第9章介绍了稀疏表达模型在超声成像中的应用；0章着重介绍了超声辐射力与微粒声操控物理原理和计算结果；1-12章主要介绍生物医学超声空化、相变与多模空化物理机制和计算仿真。

目 录
第1章 治疗超声的非线性声场1
1.1 引言1
1.2 非线性声场的理论基础与计算模型1
1.2.1 超声波传播的非线性特性1
1.2.2 理想流体中有限振幅声波的线性传播方程3
1.2.3 黏滞流体中有限振幅声波的非线性传播方程6
1.3 声波非线性传播分析14
1.3.1 时域波形畸变15
1.3.2 线性与非线性声压分布对比16
1.3.3 声波传播非线性变化提升组织热效应17
1.4 高强度聚焦超声非线性特性的影响因素18
1.4.1 换能器中心频率对非线性特性的影响18
1.4.2 超声波传播深度对非线性特性的影响19
1.4.3 传播介质的声衰减系数对非线性特性的影响20
1.5 本章小结20
主要参考文献21
第2章 治疗超声中生物组织温度场23
2.1 引言23
2.2 治疗超声中生物组织温度场的理论与计算23
2.2.1 生物组织中超声非线性传播模型23
2.2.2 治疗超声中生物组织内温升与热传导24
2.2.3 治疗超声中生物组织热剂量25
2.2.4 治疗超声中生物组织温度场计算流程26
2.3 治疗超声中软组织温度场27
2.3.1 治疗超声中软组织温度场模型27
2.3.2 肌肉与脂肪组织中温度场特性27
2.3.3 生物组织不同声衰减系数和热导率下的温度场特性29
2.3.4 不同超声换能器参数下生物组织温度场特性32
2.3.5 治疗超声中软组织热损伤特性35
2.4 治疗超声中含血管软组织温度场39
2.4.1 含血管软组织温度场模型39
2.4.2 不同血管直径下含血管软组织温度场特性41
2.4.3 不同血管位置下含血管软组织温度场特性42
2.4.4 不同血流速度时含血管软组织温度场特性43
2.5 治疗超声中骨-组织温度场44
2.5.1 治疗超声中腿骨-软组织界面温度场模型45
2.5.2 不同的腿骨与超声焦点相对位置下腿骨-软组织界面温度场特性46
2.5.3 不同超声作用参数下腿骨-软组织界面温度场特性49
2.5.4 经颅超声治疗中颅骨-脑组织温度场模型54
2.5.5 不同超声换能器参数下经颅治疗超声中颅骨-脑组织温度场特性56
2.6 本章小结59
主要参考文献59
第3章 单微泡振动与声散射62
3.1 引言62
3.2 超声场中单微泡振动理论建模与仿真62
3.2.1 超声场中微泡的线性和非线性特性62
3.2.2 无包膜自由微泡振动的理论模型64
3.2.3 包膜微泡振动的理论模型67
3.2.4 单微泡振动模型求解与散射回波计算分析70
3.3 不同微泡参数对微泡振动特性及回波特性的影响74
3.3.1 微泡初始半径对微泡振动特性及回波特性的影响74
3.3.2 微泡包膜厚度对微泡振动特性及回波特性的影响76
3.3.3 包膜弹性对微泡振动特性及回波特性的影响77
3.3.4 包膜黏性对微泡振动特性及回波特性的影响78
3.4 不同介质参数对微泡振动特性及回波特性的影响79
3.4.1 介质表面张力对微泡振动特性及回波特性的影响79
3.4.2 介质黏滞系数对微泡振动特性及回波特性的影响80
3.5 不同声场参数对微泡振动特性及回波特性的影响83
3.5.1 声驱动幅度对微泡振动特性及回波特性的影响83
3.5.2 声驱动相位对微泡振动特性及回波特性的影响85
3.5.3 激励频率对微泡振动特性及回波特性的影响87
3.5.4 声驱动脉冲个数对微泡振动特性及回波特性的影响90
3.6 本章小结92
主要参考文献93
第4章 黏弹性介质中的微泡动力学95
4.1 引言95
4.2 黏弹性介质中的微泡振动模型研究95
4.2.1 基于Maxwell力学模型的微泡振动研究95
4.2.2 基于Kelvin-Vogit力学模型的微泡振动研究98
4.3 黏弹性介质中微泡振动的Yang-Church模型100
4.3.1 Yang-Church模型的形成100
4.3.2 黏弹性介质中微泡振动的共振频率102
4.4 黏弹性介质中微泡振动特性分析102
4.4.1 激励声压对微泡振动特性的影响102
4.4.2 激励频率对微泡振动特性的影响103
4.4.3 微泡初始半径对微泡振动特性的影响104
4.4.4 介质表面张力对微泡振动特性的影响105
4.4.5 介质黏滞系数对微泡振动特性的影响105
4.5 黏弹性介质中的被动空化检测106
4.5.1 被动空化检测信号及频谱分析106
4.5.2 次谐波信号分析106
4.6 本章小结109
主要参考文献110
第5章 微泡群动力学与声散射111
5.1 引言111
5.2 理论基础及计算原理111
5.2.1 液体中双微泡模型及作用机制111
5.2.2 液体中微泡群模型构建113
5.2.3 声场中微泡群振动模型与声散射计算方法116
5.3 微泡间相互作用的Bjerknes力117
5.3.1 微泡初始半径和微泡间距对Bjerknes力的影响117
5.3.2 Bjerknes系数及其影响因素120
5.4 液体中微泡群的动力学特性123
5.4.1 微泡群模型的理论验证123
5.4.2 基于立方体和球体结构的微泡群动力学124
5.4.3 基于高斯分布的微泡群动力学126
5.5 液体中微泡群的声散射特性127
5.5.1 基于立方体结构的微泡群的声散射特性127
5.5.2 微泡数目、初始半径和微泡间距对微泡群声散射特性的影响129
5.6 本章小结131
主要参考文献132
第6章 波束合成与控制134
6.1 引言.134
6.2 常用超声波束合成技术135
6.3 超声自适应波束合成技术138
6.3.1 信号模型138
6.3.2 最小方差无失真响应法140
6.3.3 空间平滑技术141
6.3.4 稳健的波束合成技术144
6.4 改进的超声自适应波束合成技术144
6.4.1 特征空间自适应波束合成144
6.4.2 维纳自适应波束合成145
6.4.3 广义相干系数自适应波束合成146
6.5 聚焦波成像中的波束合成方法147
6.6 合成孔径成像中的波束合成方法151
6.7 平面波成像中的波束合成方法153
6.8 弱聚焦宽波束成像中的波束合成方法158
6.9 本章小结161
主要参考文献162
第7章 血管壁弹性模量重构164
7.1 引言164
7.2 弹性成像理论基础164
7.2.1 弹性力学基本理论164
7.2.2 弹性模量重构165
7.3 血管壁弹性模量重构有限元仿真170
7.3.1 血管壁弹性模量重构算法170
7.3.2 血管壁斑块模型的弹性模量重构171
7.3.3 弹性模量重构的评估和校正175
7.4 血管腔内超声（IVUS）成像的血管壁弹性重构177
7.4.1 实验材料和方法177
7.4.2 结果178
7.5 本章小结179
主要参考文献179
第8章 生物组织分数阶黏弹性测量与估计方法181
8.1 引言181
8.2 黏弹性理论及其应用181
8.2.1 整数阶黏弹性理论与模型181
8.2.2 分数阶黏弹性理论与模型184
8.3 黏弹性接触理论188
8.3.1 弹性接触理论基础188
8.3.2 黏弹性Boltzmann积分189
8.4 压痕试验192
8.4.1 试验材料制备192
8.4.2 松弛192
8.4.3 蠕变192
8.4.4 松弛和加载-卸载193
8.5 模型参数估计方法及参数物理意义193
8.5.1 模型参数估计193
8.5.2 KVFD模型参数物理意义197
8.6 本章小结198
主要参考文献198
第9章 基于稀疏表达模型的超声成像及GPU并行计算200
9.1 引言200
9.2 稀疏表达模型及GPU 并行计算200
9.2.1 稀疏表达模型及其应用201
9.2.2 稀疏表达的计算202
9.2.3 串行计算及并行计算204
9.2.4 基于GPU 的并行计算204
9.3 基于稀疏表达模型的波束合成205
9.3.1 超声成像的波束合成问题205
9.3.2 频域波束合成206
9.3.3 应用稀疏表达模型恢复波束207
9.4 基于稀疏表达模型的解卷积208
9.4.1 超声成像的解卷积问题208
9.4.2 联合稀疏表达模型提升超声图像的分辨率208
9.5 稀疏表达模型的结果与GPU并行计算的性能分析209
9.5.1 稀疏表达模型的建立及优化问题的求解209
9.5.2 稀疏表达模型在超声波束合成的应用212
9.5.3 稀疏表达模型在超声成像解卷积的应用212
9.5.4 GPU并行计算的实现与性能分析216
9.6 本章小结219
主要参考文献219
第10章 基于声辐射力的微球体操控221
10.1 引言221
10.2 基于声辐射力的微球体操控的理论基础222
10.2.1 声辐射力计算的基本原理222
10.2.2 驻波场中声辐射力的计算公式224
10.2.3 微球体与振动微泡在操控声场中的运动225
10.3 行波声场中微球体声辐射力的计算227
10.4 驻波声场中微球体声辐射力的计算228
10.4.1 微粒、微泡的声辐射力对比228
10.4.2 入射声束频率对声辐射力的影响231
10.4.3 微球体半径对声辐射力的影响233
10.5 超声驻波场中微球体运动规律的计算235
10.5.1 呈180°夹角入射的超声驻波场中微球体的运动过程235
10.5.2 呈120°夹角入射的超声驻波场中微球体的运动过程238
10.6 涡旋声场中的微球体声操控245
10.6.1 涡旋声场的研究历程245
10.6.2 圆周点声源涡旋声场的计算理论247
10.6.3 基于声涡旋的声辐射力和微粒操控249
10.6.4 涡旋声场中声操控过程的计算250
10.7 本章小结255
主要参考文献256
第11章 超声空化与声致液滴相变259
11.1 引言259
11.2 超声空化理论与计算原理259
11.2.1 成核机制259
11.2.2 振动过程264
11.2.3 坍塌过程266
11.3 声致液滴相变理论与计算原理267
11.3.1 声致液滴聚焦高次谐波相变机制268
11.3.2 声致液滴相变过程270
11.4 超声空化过程中的影响因素275
11.4.1 Blake阈值275
11.4.2 异质成核阈值275
11.4.3 空化坍塌阈值276
11.4.4 空化坍塌时间277
11.4.5 空化坍塌压力279
11.5

    章 治疗超声的非线性声场
    1.1 引言
    超声波的医学应用可以分为两种：一种是诊断超声，它是利用脉冲回波获取组织的结构图像；另一种是治疗超声，主要是应用超声的生物学效应。在治疗超声应用中，非线性的有限振幅声波产生的物理效应有声辐射压力、声流和空化效应（程建春，2012）。传统非线性声传播理论认为声速随声压的增加而增加，从而使波峰传播快而波谷传播慢，导致波形产生畸变，此时声波中除了基波外，还会衍生出谐波的成分，介质中传播的声波声压并不能采用线性声传播方程确定，而需要通过建立非线性模型和方程来求解。声波非线性的作用会使得焦点处的能量显著提升、组织吸收等热效应明显增强。由于吸收系数与频率呈现幂律关系，谐波的产生增加了超声生物效应发生的可能性。高强度聚焦超声的治疗手段正是依赖于组织的这些生物效应。由于聚焦超声在不同的强度水平所产生的生物效应不同，相应的治疗模式也不同。传统的治疗模式主要利用超声的热效应，加热的速率取决于局部的声强，组织吸热产生蛋白质热凝固变性。将超声强度进一步提高，除了热效应增强，治疗模式还可以利用伴随的达到惯性空化阈值的微泡产生的机械效应。在超高强度的超声作用下，由于焦区非线性传播的影响，超声在焦点处会形成冲击波，并因此可能产生更加强烈的空化作用，造成细胞层面的机械损伤。在现代超声治疗领域中，非线性的效应不可忽视并且已经成为治疗应用的关键技术（章东等，2014）。
    本章内容着眼于高强度超声的非线性变化及生物热效应，从理想流体质点运动方程逐步推导出适合于大振幅声波在黏滞流体中的非线性声场方程，在此基础上对方程进一步做抛物线近似处理得到目前应用很为广泛的KZK方程。而后，采用频域有限差分的方法求解该方程，仿真出声波的非线性变化、非线性声场分布及热效应结果，之后探讨影响非线性变化的因素，对比阐述线性与非线性传播的差异。这些结果可对超声治疗方案的选择和参数优化提供指导，是治疗超声，尤其是超声手术研究中很基础的内容。
    1.2 非线性声场的理论基础与计算模型
    1.2.1 超声波传播的非线性特性
    单频平面声波在无边界的理想介质中传播时，如果幅值微小，可遵循线性声波的传播理论，不存在边界反射、能量耗散和声频散等现象，波形不会发生变化。但是当初始声压逐渐增大或超声频率增大，声波传播就会偏离线性声学的传播规律，衍生出谐波成分，波形会发生畸变，基波的能量转移到谐波成分上，这就使声波传播产生了非线性变化。由频谱分析可知，圆频率为w的周期有限振幅波在耗散介质中传播会产生波形失真，可以将此时的时域波形展开成频率为w，2w，3w， ，mw的一系列简谐函数，并且随着传播距离的增加，高频分量在总能量中所占比重增大，随时间具有累加效应。因此，畸变的波形就是多个谐波成分叠加形成的结果（章东等，2014）。
    介质的密度及压缩系数决定声速，声速的改变将导致声波的折射。空间密度及压缩系数的变化（即声的特性阻抗）或者声吸收会引起声波的散射和反射。生物组织中，折射、反射、散射及吸收都是引起声能量衰减的因素。由于传播介质的非线性和耗散，有限振幅声波在传播中不断衍生出了谐波，谐波继承基波的能量，基波在传播过程中不仅会因为介质的黏滞吸收引起能量衰减，而且由于谐波衍生转移能量造成了附加衰减，因此，焦区非线性热效应就更加显著。
    正是由于声传播的非线性及非线性热效应，焦区介质吸收转化声能的能力将大大提高，局部组织温度会迅速升高，短时间内热传导现象并不明显，所以焦区局部组织可能很快会达到能量饱和的状态。声饱和效应很终将使传递到介质中指定点上的声功率受到，使接收点的超声波强度大大降低，影响超声治疗的效果。
    高强度聚焦超声（high intensity focused ultrasound，HIFU）是一种近年来发展很好迅猛的治疗方法，以其非侵入、治疗时间短等优势已经成为医学超声领域的研究热点（Kennedy，2005）。HIFU治疗手术是在体外把超声波通过一些聚焦手段（声透镜、凹球面自聚焦、电子聚焦等）汇聚，然后辐照到体内，超声波与人体靶区组织相互作用从而产生热效应、机械效应和空化效应等，瞬间（0.5～2s）灭活肿瘤组织，而且不损伤靶区周围的正常组织。
    聚焦超声波在到达靶组织中的指定区域的过程中，往往需要经过多层组织介质（包括皮肤、肌肉、脂肪等）。在每层组织界面上，由于声阻抗的差异，部分声能量被反射，造成能量的无益损耗。声能量的透射系数主要取决于介质声阻抗的差异及组织的厚度，人体组织中，除了脂肪、气体和骨头，大部分组织的声阻抗和水相似，所以，一般使用水作为耦合介质，使得超声能量尽可能少损失地从换能器传递至人体组织。超声传播通过组织时，声波能量会发生衰减，这主要有两个方面的原因：首先，声压波动引起组织微观层面上的剪切运动，入射声波所携带的部分机械能由于黏滞吸收转化为热能，这种黏滞吸收产生的热效应是超声热疗的主要机制；其次，在各向异性介质中，组织声学特性不均匀，与周围组织不同声学特性的小区域将入射声波散射至各个方向，导致声波的入射方向的能量损失。Kashcheeva等（2000）提出理论：对于大多数组织而言，衰减系数a与超声频率f两者之间呈指数关系：
    （1.1）
    式中，a0为频率f0时声波传播的衰减系数；n为组织常量。对于生物组织，f0=1MHz时，n在1～1.4之间取值，衰减系数和频率近似呈线性关系。
    由于声波传播的衰减系数与频率存在这样的对应关系，因此如何合理选取超声治疗的各种参数成了HIFU治疗研究中一项重要内容。对于多数临床治疗应用，HIFU换能器使用0.5～0.8MHz的单频超声，如果假设超声在组织中的传播遵循线性理论，则组织吸收热量接近依赖于入射超声声强及局部组织的声吸收系数。实际上，高强度超声在组织中传播时，波形会产生非线性畸变，从而衍生出谐波成分，而任何超声声波传播过程中由于非线性机制产生的高次谐波成分，都会改变线性条件下组织对于热量的吸收，从而产生许多与线性传播理论预想的不同的现象。HIFU治疗过程中有两种明显的非线性机制：非线性声传播和空化效应。我们这里主要探讨的是声波传播的非线性机制。
    当大振幅的单频声波在介质中传播时，声波波形由于与介质的相互作用逐渐发生畸变，甚至可能形成冲击波。由于组织的声吸收系数与频率存在依赖关系，因此高次谐波成分更易被组织吸收，进而转化为热能。基频能量向高次谐波转移的程度取决于入射声波的振幅、传播媒介的非线性参量（一般用B/A表示），以及声波在介质中的传播距离。在HIFU治疗的研究中，随着靶组织深度的增加和声强的提高，非线性效应的影响变得越来越重要，波形非线性畸变导致组织能量的额外吸收对HIFU治疗中的热效应有着越来越重要的影响。而根据对HIFU的很新研究，利用分裂焦点（Lu et al.，2014）、组织毁损（Xu et al.，2004）及强非线性引起的冲击波效应（Maxwell et al.，2012）等技术手段在焦区产生更加强烈的非线性生物效应对于治疗模式的改进和治疗效率的提升都大有裨益。
    1.2.2 理想流体中有限振幅声波的线性传播方程
    流体介质在宏观上是连续的，以牛顿质点动力学的观点，流体介质可以看成是由许多紧密相连的微小体积元组成的，每个体积元内的介质集中在一点。在声扰动下，质点在平衡位置附近振动，又引起周围其他介质的振动，从而使振动由近而远传播，这种传播称为声波。存在声波的空间称为声场。理想流体中，不考虑黏性，没有声能转化为热能等能量形态变化，没有总能量的损耗；在流体介质中没有切向力，也就没有横波发生。微小声波时，没有线性波形畸变导致的非线性吸收，没有定向声场引起介质定向流动的问题。在理想流体中，质点只在压力的作用下运动，不存在其他外力。
    （1.2）
    质点加速度由两部分组成：靠前部分称为局部加速度，它是空间固定点的加速度，表示空间速度场的时间变化率；第二部分称为对流加速度，它是质点从t时刻所在空间点运动至t+dt时刻空间点所获得的加速度。对于理想流体，外力为声压，根据牛顿第二定律，可以得到质点的运动方程（程建春，2012）：
    （1.3）
    质点在运动中需要遵守质量守恒定律，即单位时间内体积元质量的增加量等于流入体
    （1.4）
    另外一个方程为状态方程，它描述了介质中热力学量，包括压力P、密度p和熵S三者之间的关系，即（程建春，2012）：
    （1.5）
    在声波传播过程中，体积压缩和膨胀的周期远小于热传导需要的时间，因此在声波传播时，质点来不及与周围介质发生热交换，可以认为声波传播是绝热的，熵保持不变，在物态方程中，声压P为密度p的函数，在0附近将声压P展开成泰勒级数（程春，2012）：
    （1.6）
    式中，A为P对p的一阶偏微分，B为P对p的二阶偏微分，C为P对p的三阶偏微分，B/A称为非线性系数，描述了介质的非线性特性，是声波的非线性传播理论中很重要的参数。
    方程（1.3）～（1.5）是声学的3个基本方程，均是非线性的方程。在线性声学研究范围内，可忽略3个基本方程中的非线性项，消去方程中的p和V，推导出描述线性声场的方程（程建春，2012）：
    （1.7）
    方程（1.7）为声波传播线性模型的基本方程。式中，C0为声速。
    （1.10）
    其解为
    （1.11）
    代入式（1.9）得
    （1.12）
    式中，靠前项是沿+z方向、第二项是沿-z方向行进的平面声波。
    在HIFU治疗研究中，需要将声场能量聚焦在很小的区域内，球壳式聚焦换能器应用广泛。图1.1是基于二维轴对称圆柱坐标下的Helmholtz方程，利用Comsol软件进行的仿真计算，波动方程在二维轴对称圆柱坐标下的Helmholtz方程为（章东等，2014）：
    （1.13）
    式中，r和z分别是径向和轴向坐标；p是声压；w为角频率；密度pc和声速Cc考虑到材料的阻尼特性均是复数。
    使用方程（1.13）假设声波传播是线性的（Canino et al.，1998），而且剪切波在组织区域内的振幅要远小于压力波，因此非线性效应和剪切波被忽略。解法采用有限元方法，具体计算过程不作为本章重点内容，不做过多解释。图1.1（a）、（b）是球壳式聚焦换能器、水、组织仿体的耦合模型，四个圆柱形状的完美匹配层PML（r1～r4）和一个球形PML（c1）用于吸收出射波。图1.1（c）描述了计算中的网格划分，为了准确地解决焦点区域内的强烈的压力梯度，模型在焦区内使用λ/6波长的细网格，在其他区域使用λ/4的粗网格。图1.1（d）为声场的计算结果。仿真结果显示出超声波在组织中的传播路径和声场的分布。进一步根据Pennes（1948）生物传热原理可以得出热场仿真。
    图1.1 HIFU声场模型与仿真