1.   引　言

得益于良好的穿透性与指向性, 超声波广泛应用于医学成像与无损检测中, 并被认为具有实现非侵入性治疗的潜质, 自1942年Lynn首次在牛肝组织中造成凝固性坏死中使用高强度聚焦超声以来, 该技术已经得到了长足的发展, 并被应用于前列腺癌^[1-3]、子宫肌瘤^[4-6]、肝癌^[7,8]等多类疾病的治疗中. 由于肋骨具有显著区别于软组织的声学性质, 当使用高强度聚焦超声(HIFU)对肋后的病灶进行治疗时, 肋后区域的声场分布与无肋骨情形有明显的差异. 为研究肋骨对声场分布造成的影响, 2007年Li等^[9]提出了一种完美吸声体模型, 并通过在三维笛卡尔系中求解Khokhlov-Zabolotskaya-Kuznetsov (KZK)方程研究了肋后非线性声场, 该模型基于声束在刚进入肋骨介质时便被完全吸收, 肋骨内部的声压为零的假设, 并忽略了肋骨的声参数同软组织的差异. Yuldashev等^[10]在完美吸声体假设的基础上忽略了肋骨的厚度, 并在三维柱坐标系中使用KZK方程研究了多根肋骨遮挡下的非线性声场. 对于半张角大于16.6°的自聚焦型换能 器, Lin等^[11]基于有限厚度的完美吸声体假设, 通过在椭球系中求解椭球声束方程(SBE)^[12]研究了肋后的非线性声场, 并进一步通过Pennes方程研究了离体组织内的温度场^[13].

与此同时, 在不使用完美吸声体模型的前提下, Zhao等^[14]采用有限元法对二维近似后的肋后声场与温度场进行了计算, de Greef等^[15]采用角谱法在频域上计算了凹球面相控阵换能器在肋后激励的线性声场, Liu等^[16,17]分别使用Rayleigh积分与伪谱法在频域和时域上计算了二维球面换能器阵列在肋后激励的线性声场, 此外还研究了热剂量^[18]的空间分布模式. 这些工作将肋骨视为了强吸声体, 因而能够避免完美吸声体模型导致的误差, 但是并没有体现聚焦超声传播过程中的非线性效应.

为弥补以往工作的不足, 本文对现有的数值计算方法进行了改良, 提出了一种强吸声体模型, 该模型不再假设肋骨内部的声压为零, 而是在考虑肋骨与软组织声参数差异的基础上, 对大张角换能器在肋后区域激励的非线性声场进行了数值求解, 并使用声参数同肋骨接近的ABS塑料^[19]开展了相关实验, 与数值计算的结果进行了比较. 此外, 本文还分别将肋骨视作完美吸声体与强吸声体, 在多层介质模型内研究了肋骨位置对非线性声场的影响, 并进一步比较了所得计算结果的差异.

2.   理论模型

在本文的工作中, 一个中心开口的大张角换能器被用于激励声场, 此时可以在如图1所示的椭球系中研究声束的传播过程, 其基底坐标$\left(\sigma, \eta, \varphi \right)$同$\left (x, y, z \right)$的对应关系为$x = b \sqrt{\left(1+\sigma^2\right)\left(1-\eta^2\right)} \cos\varphi$, $y = b\sqrt{\left(1+\sigma^2\right)\left(1-\eta^2\right)} \sin\varphi$, $z = b \sigma \eta$, 其中, $\sigma \in (-\infty, \infty)$, $\eta \in [0, 1]$, $\varphi \in [0, 2\pi]$, $2 b$为椭球系中两个焦点的间距.

Kamakura等^[12]提出, 椭球系中的聚焦声场可以被分为如图2所示的球面波区域和平面波区域两部分, 二者的分界面为$\sigma = \sigma_0$, 其中换能器的焦距为d, 内口径为$2 a_0$, 外口径为$2 a_1$, 对应的半张角分别为$\gamma_0$与$\gamma_1$. 在椭球系中, 设定换能器表面的几何中心位于$\left(-\sigma_{{\mathrm{max}}}, 1, 0 \right)$处, 几何焦点则位于坐标原点, 此时如图2所示, 可知$b = d/\sigma_{{\mathrm{max}}}$.

记$p_0$为换能器表面声压的振幅, $\overline{p}=p/p_0$为归一化声压, 引入$\theta = {\mathrm{arccos}}\eta \in [0, \pi/2]$并在两个区域分别使用球面波近似和平面波近似, 则可以得到如(1)式所示的归一化SBE方程:

$$\begin{split}& \frac{\partial ^2\overline{p}}{\partial \sigma \partial \tau} + A \cdot \frac{\partial ^2\overline{p}}{\partial \theta \partial \tau} + B \cdot \left ( \frac{\partial ^2\overline{p}}{\partial \theta ^2} + \cot \theta \frac{\partial \overline{p}}{\partial \theta} \right )\\ &+ D \cdot \frac{\partial \overline{p}}{\partial \tau} + E \cdot \frac{\partial ^2\overline{p}}{\partial \varphi ^2} = F \cdot \frac{\partial ^3\overline{p}}{\partial \tau ^3} + G \cdot \frac{\partial ^2\overline{p}}{\partial \tau ^2}. \end{split}$$

(1)

(1)式中各变量的定义为

当$\sigma < \sigma_0$时:

$$\begin{split} &\left\{\begin{aligned} \tau &= \omega \left(t + \frac{b\sqrt{\sigma^2 + \sin^2\theta}}{c_0} \right),\\A &= \frac{\sin 2\theta}{2\sigma (1 + \sigma^2)}, \\ B &= \frac{\varepsilon \sqrt{\sin^2 \theta + \sigma^2}}{\sigma \left (1 + \sigma^2 \right )} ,\\ D &= \frac{\sigma^2 + \cos^2\theta}{\sigma \left (1 + \sigma^2 \right )}, \\ E &= \frac{\sigma^2 + \cos^2\theta}{\left (1 + \sigma^2 \right )\sin^2\theta} \cdot B,\\ F &= \frac{-\alpha b \left(\sigma^2 + \cos^2\theta \right)\sqrt{\sin^2 \theta + \sigma^2}}{\sigma \left (1 + \sigma^2 \right )}, \\ G &= \frac{-b \left(\sigma^2 + \cos^2\theta \right)\sqrt{\sin^2 \theta + \sigma^2}}{2l_D\sigma \left (1 + \sigma^2 \right )}. \end{aligned} \right. \end{split}$$

当$\sigma \geqslant \sigma_0$时:

$$\begin{split} & \left\{\begin{aligned} \tau &= \omega \left(t - \frac{b\sigma\cos\theta}{c_0} \right),\\ A&= \frac{-\sigma\sin\theta}{\left (1 + \sigma^2 \right )\cos\theta}, \\ B& = \frac{-\varepsilon}{\left (1 + \sigma^2 \right )\cos\theta}, \\ D&= 0, \\ E& = \frac{\sigma^2 + \cos^2\theta}{\left (1 + \sigma^2 \right )\sin^2\theta} \cdot B,\\ F &= \frac{\alpha b \left(\sigma^2 + \cos^2\theta \right)}{\left (1 + \sigma^2 \right )\cos\theta},\\ G &= \frac{-b \left(\sigma^2 + \cos^2\theta \right)}{2l_D\left (1 + \sigma^2 \right )\cos\theta}. \end{aligned} \right. \end{split}$$

(2)

这里α为声衰减系数, ω为换能器工作信号的角频率, $\rho_0$为介质密度, $c_0$为介质声速, β为介质的非线性系数, $l_{\mathrm{D}} = \rho_0 c_0^3/(p_0\beta\omega)$为平面波的冲击波形成距离, $\varepsilon = c_0/(2\omega b)$. 在球面波区域和平面波区域分别对归一化声压$\overline{p}$进行傅里叶级数展开, 得

$$\overline{p} = \left\{\begin{aligned} &\sum\limits_{n = -\infty}\limits^{\infty} C_n^{(s)} \cdot {\mathrm{e}}^{{\mathrm{j}}n\tau} ,&\sigma < \sigma_0 ,\\&\sum\limits_{n = -\infty}\limits^{\infty} C_n^{(p)} \cdot {\mathrm{e}}^{{\mathrm{j}}n\tau},& \sigma \geqslant \sigma_0 ,\end{aligned} \right.$$

(3)

其中$C_{-n} = C^{*}_{n}$, 在曲面$\sigma = \sigma_0$上, 有

$$C_n^{(p)} = C_n^{(s)} {\cdot} \exp \left[{\mathrm{j}}nkb \left( \sqrt{\sigma_0^2 + \sin^2\theta} + \sigma_0 \cos\theta \right) \right].$$

(4)

隐去(3)式中的上标, 并将其代入(1)式, 则可得到如下方程:

$$\begin{split} \frac{\partial C_n}{\partial \sigma} \;&= {\mathrm{j}}nG \cdot \sum\limits_{m = -\infty}\limits^{\infty} C_m C_n - n^2F \cdot C_n\\&- \Bigg[ A \cdot \frac{\partial C_n}{\partial \theta} + \frac{B}{jn} \cdot \left( \frac{\partial^2 C_n}{\partial \theta^2} + \cot\theta \frac{\partial C_n}{\partial \theta} \right) \\ &+ D \cdot C_n + \frac{E}{{\mathrm{j}}n} \cdot \frac{\partial^2 C_n}{\partial \varphi^2}\Bigg].\\[-1pt] \end{split}$$

(5)

通过求解(5)式可以在频域上得到非线性声压的分布情况, 声场的边界条件为

$$\begin{aligned} C_n \left (-\sigma_{{\mathrm{max}}},\theta,\varphi \right)= \begin{cases} 0.5{\mathrm{j}},& \theta_0 \leqslant \theta \leqslant \theta_1 \text{且} n = 1, \\ 0,& \text{其他}, \end{cases} \end{aligned}$$

(6)

其中,

$$\begin{split} &\theta_0 = \arctan \left(\tan\gamma_0/\sqrt{1+\sigma_{{\mathrm{max}}}^{-2}} \;\right),\\&\theta_1 = \arctan \left(\tan\gamma_1/\sqrt{1+\sigma_{{\mathrm{max}}}^{-2}} \;\right). \end{split}$$

3.   数值算法与组织模型

3.1   数值算法

设定换能器的几何焦点为坐标原点, 则图3(a)所示的组织模型在椭球系中关于平面$\varphi = 0$对称, 故而可以将计算区域缩小到$\varphi \in [0, \pi]$. 由(2)式可知$E\left(\sigma, \theta, \varphi \right)$在$\theta = 0$处不收敛, 故而在分割空间时可以取$\theta_s = (s-1/2) \cdot {\mathrm{d}}\theta$以避免数值振荡, 再令$\sigma_r = -\sigma_{{\mathrm{max}}} + (r-1) \cdot {\mathrm{d}}\sigma$, $\varphi_p = (p-1) \cdot {\mathrm{d}}\varphi$, 则可以通过$\left (\sigma_r, \theta_s, \varphi_p \right)$对(5)式进行离散化处理, 其中$r = 1, 2, \cdots, N_\sigma$, $s = 1, 2, \cdots, N_\theta$, $p = 1, 2, \cdots, N_\varphi$.

为避免使用有限差分法求解衰减效应时可能出现的数值振荡问题, 可以在数值计算的过程中使用算子分离法^[20], 即在微小步长${\mathrm{d}}\sigma$上依次求解非线性项、衍射项与衰减项的贡献. 首先计算非线性效应的贡献, 对应的方程为

$$\begin{aligned} \frac{\partial C_n}{\partial \sigma} = {\mathrm{j}}nG \cdot \sum\limits_{m = -\infty}\limits^{\infty} C_m C_n. \end{aligned}$$

(7)

由于数值计算中仅能考虑有限数量的谐波, 且非线性效应存在累积性, 在曲面$\sigma = \sigma_r$上可以仅考虑前$M_r$个谐波, 此时由(7)式可得:

$$\begin{aligned} C_{r+1,s,p}^{(n)} = \left\{\begin{aligned} &C_{r,s,p}^{(n)} + {\mathrm{j}}nG_{r,s,p} \cdot \sum\limits_{m = -M_r}\limits^{M_r} C_{r,s,p}^{m}C_{r,s,p}^{n-m},\\&\qquad n \leqslant M_r\\& {\mathrm{j}}nG_{r,s,p} \cdot \sum\limits_{m = -M_r}\limits^{M_r} C_{r,s,p}^{m}C_{r,s,p}^{n-m},\\& \qquad n > M_r \end{aligned}\right. \end{aligned}$$

(8)

其中$C_{r, s, p}^{(n)} = C_n\left (\sigma_r, \theta_s, \varphi_p \right)$, $G_{r, s, p} = G\left (\sigma_r, \theta_s, \varphi_p \right)$, 后文中出现的离散变量同理. 完成非线性项的计算后, 将曲面$\sigma = \sigma_{r+1}$的所有谐波的振幅同小量$\varepsilon$比较, 记最大振幅不小于$\varepsilon$的谐波阶数为N, 并取$M_{r+1} = \max \left(M_r, N \right)$, 则可以避免截止频率过低导致的能量反向传播, 在本文的工作中, $\varepsilon = \left(2 p_0\right)^{-1}$, $M_1 = 8$. 完成非线性项的计算后考虑衍射项的贡献, 其对应的方程为

$$\begin{split} \frac{\partial C_n}{\partial \sigma} = \;&- \Bigg[ A \cdot \frac{\partial C_n}{\partial \theta} + \frac{B}{{\mathrm{j}}n} \cdot \left( \frac{\partial^2 C_n}{\partial \theta^2} + \cot\theta \frac{\partial C_n}{\partial \theta} \right) \\ &+ D \cdot C_n + \frac{E}{{\mathrm{j}}n} \cdot \frac{\partial^2 C_n}{\partial \varphi^2}\Bigg].\\[-1pt] \end{split}$$

(9)

记(9)式中等号右侧的项为$-H_n$, 并使用中心差分格式对其进行计算, 得

$$\begin{split} H_{r,s,p}^{(n)} = \;&A_{r,s,p} \cdot \frac{C_{r,s+1,p}^{(n)}-C_{r,s-1,p}^{(n)}}{2{\mathrm{d}}\theta}\\&+ \frac{B_{r,s,p}}{{\mathrm{j}}n} \cdot \Bigg[ \frac{C_{r,s+1,p}^{(n)} - 2C_{r,s,p}^{(n)} + C_{r,s-1,p}^{(n)}}{\left( {\mathrm{d}}\theta \right)^2}\\ &+ \cot\theta_s \frac{C_{r,s+1,p}^{(n)}-C_{r,s-1,p}^{(n)}}{2{\mathrm{d}}\theta} \Bigg]\\&+ D_{r,s,p} \cdot C_{r,s,p}^{(n)} + \frac{E_{r,s,p}}{{\mathrm{j}}n} \\ &\times \frac{C_{r,s,p+1}^{(n)} - 2C_{r,s,p}^{(n)} + C_{r,s,p-1}^{(n)}}{\left( {\mathrm{d}}\varphi \right)^2}.\\[-1pt] \end{split}$$

(10)

定义矩阵${\boldsymbol{C}}^{(n)}$满足${\boldsymbol{C}}_{s, p}^{(n)} = C_{r, s, p}^{(n)}$, 则可以设出满足如下关系的矩阵${\boldsymbol{J}}$, ${\boldsymbol{K}}$与${\boldsymbol{R}}$:

$$\begin{aligned} \left\{\begin{aligned} \left({\boldsymbol{J}}\cdot {\boldsymbol{C}}^{(n)} \right)_{s,p} &= \frac{C_{r,s+1,p}^{(n)}-C_{r,s-1,p}^{(n)}}{2{\mathrm{d}}\theta},\\\left({\boldsymbol{K}}\cdot {\boldsymbol{C}}^{(n)} \right)_{s,p} &= \frac{C_{r,s+1,p}^{(n)} - 2C_{r,s,p}^{(n)} + C_{r,s-1,p}^{(n)}}{\left( {\mathrm{d}}\theta \right)^2} \\ &+ \cot\theta_s \frac{C_{r,s+1,p}^{(n)}-C_{r,s-1,p}^{(n)}}{2{\mathrm{d}}\theta},\\\left({\boldsymbol{C}}^{(n)}\cdot {\boldsymbol{R}}\right)_{s,p} &= \frac{C_{r,s,p+1}^{(n)} - 2C_{r,s,p}^{(n)} + C_{r,s,p-1}^{(n)}}{\left( {\mathrm{d}}\varphi \right)^2}. \end{aligned} \right. \end{aligned}$$

(11)

进一步地, 定义${\boldsymbol{H}}^{(n)}$满足${\boldsymbol{H}}_{s, p}^{(n)} = H_{r, s, p}^{(n)}$, ${\boldsymbol{A}}$满足${\boldsymbol{A}}_{s, p} = A_{r, s, p}$, ${\boldsymbol{B}}$满足${\boldsymbol{B}}_{s, p} = B_{r, s, p}$, ${\boldsymbol{D}}$满足${\boldsymbol{D}}_{s, p} = D_{r, s, p}$, ${\boldsymbol{E}}$满足${\boldsymbol{E}}_{s, p} = E_{r, s, p}$, 则(10)式可以被改写为如下的矩阵形式:

$$\begin{split} {\boldsymbol{H}}^{(n)} = \;&{\boldsymbol{A}} \circ \left({\boldsymbol{J}}\cdot {\boldsymbol{C}}^{(n)} \right) + \frac{{\boldsymbol{B}}}{{\mathrm{j}}n} \circ \left({\boldsymbol{K}}\cdot {\boldsymbol{C}}^{(n)} \right)\\ &+ {\boldsymbol{D}} \circ {\boldsymbol{C}}^{(n)} + \frac{{\boldsymbol{E}}}{{\mathrm{j}}n} \circ \left( {\boldsymbol{C}}^{(n)} \cdot {\boldsymbol{R}} \right), \end{split}$$

(12)

其中$\circ$表示矩阵的哈达玛积(Hadamard product), ${\boldsymbol{J}}$与${\boldsymbol{K}}$的尺寸为$N_\theta \times N_\theta$, ${\boldsymbol{R}}$的尺寸为$N_\varphi \times N_\varphi$, 其余矩阵的尺寸均为$N_\theta \times N_\varphi$. 由于(12)式在结构上耦合了哈达玛积与一般积, 传统的矩阵运算方式无法在后向隐式差分格式下通过该方程实现对(9)式的求解. 在文献[9, 11]中, 声参数空间分布的非轴对称性被忽略, 因而${\boldsymbol{A}}$, ${\boldsymbol{B}}$, ${\boldsymbol{D}}$, ${\boldsymbol{E}}$可以被改写为大小为$N_\theta \times N_\theta$的对角矩阵, 哈达玛积也相应被改写为一般积.

在改写哈达玛积为一般积的过程中, 为充分考虑声参数在空间分布上的非轴对称性, 可对(12)式进行向量化处理, 记$\vec{{\boldsymbol{A}}} = {\mathrm{vec}}({\boldsymbol{A}})$, $\vec{{\boldsymbol{B}}} = {\mathrm{vec}}({\boldsymbol{B}})$, $\vec{{\boldsymbol{D}}} = {\mathrm{vec}}({\boldsymbol{D}})$, $\vec{{\boldsymbol{E}}} = {\mathrm{vec}}({\boldsymbol{E}})$, $\vec{{\boldsymbol{C}}}^{(n)} = {\mathrm{vec}}({\boldsymbol{C}}^{(n)})$, $\vec{{\boldsymbol{H}}} = {\mathrm{vec}}({\boldsymbol{H}}^{(n)})$, 则有

$$\begin{split} \vec{{\boldsymbol{H}}}^{(n)} =\;& {\boldsymbol{A}}^{({\mathrm{d}})} \cdot {\mathrm{vec}}\left({\boldsymbol{J}}\cdot {\boldsymbol{C}}^{(n)} \right) + \frac{{\boldsymbol{B}}^{({\mathrm{d}})}}{{\mathrm{j}}n} \\ &\cdot {\mathrm{vec}}\left({\boldsymbol{K}}\cdot {\boldsymbol{C}}^{(n)} \right) + {\boldsymbol{D}}^{({\mathrm{d}})} \cdot \vec{{\boldsymbol{C}}}^{(n)}\\ &+ \frac{{\boldsymbol{E}}^{({\mathrm{d}})}}{{\mathrm{j}}n} \cdot {\mathrm{vec}}\left( {\boldsymbol{C}}^{(n)} \cdot {\boldsymbol{R}} \right), \end{split}$$

(13)

其中${\boldsymbol{A}}^{({\mathrm{d}})}$是主对角线为$\vec{{\boldsymbol{A}}}$的对角矩阵, ${\boldsymbol{B}}^{({\mathrm{d}})}$是主对角线为$\vec{{\boldsymbol{B}}}$的对角矩阵, ${\boldsymbol{D}}^{({\mathrm{d}})}$与${\boldsymbol{E}}^{({\mathrm{d}})}$同理. 考虑到${\mathrm{vec}}\left({\boldsymbol{L}}{\boldsymbol{X}}{\boldsymbol{R}} \right) = \left({\boldsymbol{R}}^{\mathrm{T}} \otimes {\boldsymbol{X}} \right) \cdot {\boldsymbol{L}}$, (13)式可以被进一步改写为如下形式:

$$\begin{split} \vec{{\boldsymbol{H}}}^{(n)} = \;& {\boldsymbol{A}}^{({\mathrm{d}})} {\cdot} \left( I_{N_\varphi} \otimes {\boldsymbol{J}} \right) \cdot \vec{{\boldsymbol{C}}}^{(n)} + \frac{{\boldsymbol{B}}^{({\mathrm{d}})}}{{\mathrm{j}}n} {\cdot} \left( I_{N_\varphi} \otimes {\boldsymbol{K}} \right) \cdot \vec{{\boldsymbol{C}}}^{(n)} \\ & + {\boldsymbol{D}}^{({\mathrm{d}})} \cdot \vec{{\boldsymbol{C}}}^{(n)} + \frac{{\boldsymbol{E}}^{({\mathrm{d}})}}{{\mathrm{j}}n} {\cdot} \left({\boldsymbol{R}}^T \otimes I_{N_\theta} \right) {\cdot} \vec{{\boldsymbol{C}}}^{(n)}, \\[-8pt]\end{split}$$

(14)

其中$\otimes$表示矩阵的克罗内克积(Kronecker product), 通过将(14)式代入(9)式中, 即可在后向隐式差分格式下通过矩阵运算实现对衍射项的求解, 所得结果将作为初始条件被用于求解衰减项的贡献, 其中同衰减效应对应的方程为

$$\begin{aligned} \frac{\partial C_n}{\partial \sigma} = - n^2F \cdot C_n. \end{aligned}$$

(15)

考虑到(15)式存在形如$C_{r+1, s, p}^{n} = C_{r, s, p}^{n} \cdot {\mathrm{exp}}(-n^2 F_{r, s, p}\cdot {\mathrm{d}}\sigma )$的解析解, 仅需将初始条件代入解析解中即可求出衰减项的贡献, 从而完成声场的计算. 需要注意的是由于多层介质模型中存在折射与透射现象, 当声束自一类介质入射到另一类介质时, 需要对声场进行相应的修正^[21].

3.2   组织模型与算法验证

使用如图3所示的几何模型描述包含肋骨的生物组织, 超声换能器位于生物组织下方, 其几何焦距、外口径与内口径分别为180 mm, 215 mm和40 mm, 工作频率为0.57 MHz. 声束在依次通过水层和包含多根肋骨的脂肪层后进入肝脏, 其中肋骨的结构被近似为长方体^[9]. 图3(b)为组织模型的$xz$截面图, 其中脂肪层的厚度为H, 肋骨的厚度为$h_\mathrm{B}$, 肋宽与肋隙分别为$w_\mathrm{B}$与$w_\mathrm{S}$, 参考目前开展的活体实验数据, 取$H =\mathrm{ 18\; mm}$, $h_\mathrm{B}= 8 \;{\mathrm{mm}}$, $w_\mathrm{B } = 15\; {\mathrm{mm}}$, $w_\mathrm{S} {= 15\; {\mathrm{mm}}}$.

针对图3所示的模型, 取$\sigma_{{\mathrm{max}}} = 10$, $\sigma_0=-0.5$, 并在3个方向上分别设定步长为${\mathrm{d}}\sigma = 2.5 \times 10^{-3}$, ${\mathrm{d}}\theta = \pi/1000$以及${\mathrm{d}}\varphi = \pi/150$, 考虑到肋骨在y方向上的长度显著大于换能器在对应层面的投影面直径, 数值计算中可以近似地将肋骨在该方向上的长度视作无穷, 记$H_\mathrm{{FS}}$为换能器几何焦点在模型内的深度, $x_\mathrm{c}$为换能器轴线同肋骨阵列中心的距离, 则换能器相对于肋骨阵列的位置能够通过这两个变量进行描述.

为验证数值算法的合理性, 使用ABS塑料制作仿肋模型在脱气水中进行了声场的测量, 图4(a)为声场扫描系统的示意图, 其中水听器与超声换能器分别被固定在水箱底部与步进电机上, 代替肋骨的ABS塑料则按照图4(b)所示的方式被固定在超声换能器下方, 在0.57 MHz下分别测定ABS塑料的声速与声衰减系数为2189 m/s与38 Np/m, 其密度为1200 kg/m³, 仿肋模型的几何参数同图3(b)中的肋骨一致.

在图4所示的模型中$H = h_\mathrm{B}$, 设定$x_\mathrm{c} {= 0\; {\mathrm{mm}}}$并改变$H_\mathrm{{FS}}$的大小, 记录焦平面内的声压峰峰值, 所得结果如图5所示. 研究数据表明, 当完美吸声体模型应用于预测仿肋模型后的声场时, 其预测的声压值普遍低于实际测量值, 且随着$H_\mathrm{{FS}}$值的降低, 其低估程度越发显著, 这表明该模型在描述低$H_\mathrm{{FS}}$情形下声场时存在明显的局限性. 与此同时, 采用强吸声体模型得到的计算结果在$H_\mathrm{{FS}}$不同的多类情形下都和测量结果吻合较好, 这表明强吸声体模型相较于完美吸声体模型能够更好地描述声场的分布.

4.   仿真结果与测温实验

4.1   非线性声场

使用表1^[11,13,22]中的声参数进行数值计算, 其中介质对不同频率声波的声衰减系数为$\alpha\left(f \right) = \alpha\left(f_0 \right)\cdot \left(f/f_0 \right)^\mu$. 对于图3所示的模型, 当$x_\mathrm{c}$从0 mm逐步增大至15 mm时, 换能器轴线从肋骨的中心逐步移动至肋隙的中心, 取$p_0 = \mathrm{0.4\; MPa}$并计算$H\mathrm{_{FS} = 70\; mm}$时的非线性声场, 图6描述了基波振幅在z轴上的分布情况, 其中Field Ref表示肋骨不存在时的计算结果. 图6表明, 换能器轴线相对于肋骨的位置会影响基波振幅沿z轴的分布模式, 其中$x\mathrm{_c = 0\; mm}$与$x\mathrm{_c = 5 \;mm}$时的分布模式相近, $x\mathrm{_c = 10\; mm}$与$x\mathrm{_c }= 15\; {\mathrm{mm}}$时的分布模式相近, 为进一步研究肋骨遮挡对谐波的影响, 记第n阶谐波在z轴上的最大振幅与对应的位置为${A_n}$和${z_n}$, 所得结果列于表2. 可知, 尽管谐波与基波都能会聚在几何焦点附近, 但基波因肋骨遮挡而出现的焦点偏移现象更加显著, 且更容易受到$x\mathrm{_c}$的影响, 此外由于基波振幅的下降导致非线性效应减弱, 谐波振幅的下降程度显著地强于基波. 鉴于声焦点与几何焦点的偏差不超过1 mm, 近似地使用声场在平面$\sigma = 0$上的分布情况来描述焦平面上的声场, 在该平面内, 分别记第n阶谐波在x轴与y轴上的–6 dB宽度为$\mathrm{WX}_{n, -6\;{\mathrm{dB}}}$与$\mathrm{WY}_{{n, {-6\; {\mathrm{dB}}}}}$, 具体数据见表3.

由表3可知, 肋骨的遮挡在不同程度上改变了声压沿x轴与y轴的–6 dB宽度, 在换能器轴线由肋骨中心移动到肋隙中心的过程中, $\mathrm{WX}_{n, -6 \;{\mathrm{dB}}}$呈现出递减的趋势, $\mathrm{WY}_{n, -6\; {\mathrm{dB}}}$则总体上呈现出递增趋势, 当换能器轴线恰好位于肋隙中心时, $\mathrm{WX}_{n, -6\; {\mathrm{dB}}}$显著小于$\mathrm{WY}_{n, -6 \;{\mathrm{dB}}}$.

4.2   热沉积速率

超声波在介质内传播的过程中, 一部分声能会被吸收并导致介质温度的升高, 这一现象被长期作为HIFU治疗的主要作用机制. 通常可以使用热沉积速率Q^[23]描述聚焦超声引发的热效应, 对于本文研究的非线性声场, 存在如下关系:

$$\begin{aligned} Q = \sum_{n > = 1}2\alpha_n I_n = \sum_{n > = 1}4\alpha_n \cdot \frac{p_0^2}{\rho c} \cdot \left|C_n\right|^2, \end{aligned}$$

(16)

其中$\alpha_n$表示介质对第n阶谐波的声衰减系数, 一般情况下高次谐波对应于更大的声衰减系数, 在本文设计的计算模型中, Q在曲面$\sigma = \sigma_r$上的取值总共需要考虑$M_r$个谐波的贡献, 图7描述了$H\mathrm{_{FS}} = 70 \;{\mathrm{mm}}$时热沉积速率在截面$\varphi = 0$内的分布情况.

由图7可知肋后出现了明显的声影区域, 而焦点附近则出现了焦域裂化现象^[24], 记热焦域的最大热沉积速率为$Q_{\mathrm{m}}$, 并进一步研究$H\mathrm{_{FS}}$对$Q_{\mathrm{m}}$的影响以及完美吸声体模型带来的误差, 所得结果如图8所示, 其中图8(a)为采用完美吸声体模型的计算结果, 图8(b)为采用强吸声体模型的计算结果. 如图8(a), (b)所示, 无肋骨情形下$Q_{\mathrm{m}}$随着$H\mathrm{_{FS}}$的增大而单调下降, 这是由于声束在传播过程中的损耗因$H\mathrm{_{FS}}$的增大而加剧, 最终使得到达焦域的能量减小. 当肋骨存在时, $Q_{\mathrm{m}}$受到$H\mathrm{_{FS}}$与$x\mathrm{_c}$两个因素的综合影响, 对于本文使用的中心开孔式自聚焦换能器, 不论使用完美吸声体模型还是强吸声体模型进行计算, 换能器轴线对准肋骨中心的情形都在${45\; {\mathrm{mm}} \leqslant H_{{\mathrm{FS}}} \leqslant 65\; {\mathrm{mm}}}$的范围内具有最大的$Q_{\mathrm{m}}$, 但采用完美吸声体模型得到的$Q_{\mathrm{m}}$仅有采用强吸声体模型所得结果的50%—70%.

4.3   测温实验

考虑到超声照射的过程中, 焦点附近的温度升高主要由热沉积速率Q导致, 可以通过测定该区域的温度变化趋势分析肋骨遮挡对热沉积速率的影响. 对于聚焦超声引发的温度升高现象, 可以通过以下Pennes方程^[23]对离体组织内焦域的温度变化进行研究:

$$\begin{aligned} \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + Q. \end{aligned}$$

(17)

式中, $C_p$与k分别为比热容与导热系数, Q为根据(16)式计算得到的热沉积速率, T为温度. 为进一步验证强吸声体模型的合理性, 本文设计了如图9所示的测温组织模型, 其中图9(a)为由仿肋塑料柱、透声膜与有机玻璃盒三部分组成的定制容器, 图9(b)为模型实物图, 模型中使用的仿肋塑料柱由ABS塑料制成, 其宽度与间隙均被设定为15 mm.

在超声照射开始前, 嵌有热电偶的仿体凝胶被放置在有机玻璃盒的底部, 经脱气处理的新鲜猪肝则被放置在凝胶上方, 其中热电偶自下而上地穿出凝胶并插入猪肝组织. 为确保实验过程中凝胶与猪肝组织的位置不变, 本文使用橡皮筋将透声膜固定在了有机玻璃盒的上方, 并根据实验需求将仿肋塑料柱固定在了透声膜的上方, 此时热电偶到塑料柱上端的距离约为70 mm, 数值计算中对猪肝设定的介质声参数同表1一致, 比热容与导热系数则分别为3700$\; \mathrm{J/(kg\cdot K)}$与0.5$\;\mathrm{W/(m\cdot K)}$^[13].

完成测温组织模型的组装后, 将换能器的几何焦点置于热电偶附近进行10 s的超声照射, 并记录超声照射前后焦域的温度变化情况. 为研究肋骨遮挡造成的影响, 首先在声通道不存在仿肋塑料柱的自由场模型中进行超声照射, 然后在热电偶温度恢复至超声照射前的状态后, 将仿肋塑料柱固定在透声膜上方重复实验, 并根据自由场模型中的最大温升进行归一化处理, 所得结果如图10所示. 研究数据表明, 强吸声体模型能更好地描述肋骨遮挡对焦域温升的影响, 而完美吸声体模型则会导致对焦域温升的低估.

5.   结　论

本文提出了一种将肋骨视作强吸声体而非完美吸声体的计算模型, 并相应改进了现有的SBE方程求解方法, 以避免数值计算过程中潜在的数值振荡风险. 为了比较两类模型的优劣, 进一步使用ABS塑料构建了仿肋模型并进行了声场测量. 研究数据表明, 完美吸声体模型无法合理预测肋骨靠近焦点情况下的声场, 而强吸声体模型则能够有效弥补这一不足. 在此基础上, 本文利用改进的模型计算了大张角换能器在生物组织内激励的非线性声场. 计算结果显示, 尽管肋骨的遮挡会导致声压振幅的下降, 声束仍能较好地聚焦在换能器的几何焦点附近. 其中, 高次谐波振幅因肋骨遮挡而出现的下降程度显著高于低次谐波和基波. 此外, 声场在几何焦平面上的–6 dB宽度因肋骨的遮挡在不同方向上呈现出不同的变化趋势.

进一步地, 本文计算了热沉积速率的空间分布模式, 并比较了完美吸声体模型与强吸声体模型所得计算结果的差异. 计算结果表明, 肋骨的遮挡显著降低了热焦域的最大热沉积速率, 而采用完美吸声体假设则可能导致对下降程度的高估, 从而在临床中引发健康脏器的热损伤, 这一能量低估问题可以通过已开展的测温实验得到验证. 考虑到肝脏区域血供丰富, 在研究HIFU的热效应时, 还需考虑血流灌注的影响. 后续工作将围绕血流灌注的情形进行温度场的计算, 以指导临床治疗参数的选取.