CN103776524B - 一种测量强聚焦超声非线性声场分布的方法 - Google Patents

Abstract

本发明公开了一种测量强聚焦超声非线性声场分布的系统及其方法，属于超声声场测量技术领域。本发明步骤为：一、获得球壳换能器轴向及其焦平面径向的声压分布曲线；二、根据瑞利积分计算球壳换能器轴向及焦平面径向的声压分布曲线，将计算结果与测量结果比较，获得有效参数；三、将有效参数代入非线性声传播模型，计算R₂₁，再改变P₀，计算不同P₀对应的R₂₁，得R₂₁与P₀的关系曲线；四、对球壳换能器施加激励，测量焦点波形，计算得R₂₁′，插值得到R₂₁′对应的P₀′；五、将P₀′带入非线性声传播模型，得到非线性声场分布。本发明有效解决了水听器在测量强非线性声场带宽不足的问题，使用相对值R₂₁作为判断依据，消除了水听器灵敏度不高带来的测量隐患。

Description

一种测量强聚焦超声非线性声场分布的方法

技术领域

本发明涉及超声声场测量技术领域，更具体地说，涉及一种测量强聚焦超声非线性声场分布的系统及其方法。

背景技术

近十多年来，高强度聚焦超声（HIFU）作为一种新兴的非侵入式治疗肿瘤方法得到人们极大的关注。HIFU治疗肿瘤的主要原理是将体外超声能量聚焦于体内肿瘤组织内，靶区组织吸收声能量，温度在短时间内上升至65℃以上，发生不可逆凝固性坏死。为确保治疗的安全性和有效性，治疗前必须准确描述HIFU声场分布。由于HIFU治疗时声能量很高，焦点声压甚至能达到100MPa量级；强非线性导致声波波形严重畸变，甚至形成冲击波，频谱展宽。强声压和宽频带是HIFU声场测量中的两个挑战，要求声场测量系统既要能承受强声压又要有很宽的频带，而目前一直没有理想的声场测量系统能够满足上述两点。

传统的测量声场方法（例如使用PVDF水听器测量）无法满足HIFU声场的极端条件，目前比较常用的方法是使用光纤水听器，利用光纤端口的声-光效应（如压力变化导致光反射系数变化）测量HIFU声场。虽然光纤水听器灵敏度不如传统PVDF水听器高且噪声较大，但是其频带宽且鲁棒性好。然而，在实际测量中发现，当冲击波产生后，光纤水听器带宽仍然不能满足测量要求，测得的峰值正压常常小于实际值。此外，研究者们还提出了一些非侵入式的通过监测HIFU声场中质点位移的测量声场技术，如通过磁共振MRI（magneticresonanceimaging）、粒子速度成像PIV(particleimagevelocimetry)等，但这些方法还很不成熟，且造价昂贵。

另一方面，HIFU非线性声传播模型目前已经较为成熟。Westervelt方程描述了精确到二阶的非线性声场。对于平面和弱聚焦换能器（半张角小于16°），在近轴近似下，Westervelt方程可以化简为著名的Khokhlov–Zabolotskaya–Kuznetsov(KZK)方程；对于强聚焦换能器，在椭球坐标系下，Westervelt方程可以化简为Spheroidalbeamequation(SBE)模型。在已知边界条件，即换能器表面声压分布后，非线性声场模型可以准确的预测HIFU声场分布。

目前，有研究者提出了结合实验测量与数值模型来确定高强度聚焦超声非线性声场的方法。Canney等人于2008年提出了该方法：首先根据线性声场测量和数值计算结果得到换能器有效参数；其次根据测量的焦点声压和换能器聚焦增益得到换能器表面声压P₀；假设激励换能器的电信号V₀与换能器表面声压P₀成线性关系，那么在任意电信号V激励下，换取表面声压p=P₀(V/V₀)；最后将表面声压p作为边界条件，可以计算得到在该激励下的声场分布。然而，此方法的缺点在于当激励电信号很强时，电信号与换能器表面声压的线性关系无法得到保证。

发明内容

1.发明要解决的技术问题

本发明的目的在于克服现有高强度聚焦超声非线性声场的测量方法，在强非线性条件下，测量系统带宽不够，导致声场测量准确度不高的问题，提供了一种测量强聚焦超声非线性声场分布的系统及其方法，本发明提供的技术方案测量准确度高、容易实施，测量过程中对水听器的灵敏度要求较低，避免了因水听器灵敏度低导致检测结果不准确的隐患。

2.技术方案

为达到上述目的，本发明提供的技术方案为：

本发明的一种测量强聚焦超声非线性声场分布的方法，其测量步骤为：

步骤一、在线性声场条件下，测量球壳换能器轴向以及该球壳换能器焦平面径向的声压，获得球壳换能器轴向及其焦平面径向的声压分布曲线；

步骤二、根据瑞利积分计算球壳换能器轴向及该球壳换能器焦平面径向的声压分布曲线，将计算结果与步骤一获得的测量结果相比较，调节球壳换能器的张口半径a和几何焦距d，使得计算获得的球壳换能器轴向及其焦平面径向声压分布曲线与步骤一的测量结果相吻合，获得球壳换能器有效参数a和d；

步骤三、将步骤二获得的有效参数代入非线性声传播模型，计算焦点的二次谐波与基波的比值R₂₁，再改变球壳换能器表面声压幅值P₀，计算不同表面声压幅值P₀对应的R₂₁，得R₂₁与P₀的关系曲线；

步骤四、对球壳换能器施加单频电信号激励，测量焦点波形，对该焦点波形进行傅里叶变换，计算得R₂₁′，根据步骤三得到的R₂₁与P₀的关系曲线进行插值，得到R₂₁′对应的球壳换能器表面声压P₀′；

步骤五、将步骤四获得的表面声压P₀′带入步骤三所述的非线性声传播模型，计算得到步骤四施加的电信号激励对应的非线性声场分布。

更进一步地，步骤一测量声压的条件为：对球壳换能器施加低能量激励，使得球壳换能器焦点处基波幅度与二次谐波幅度的差值不低于40dB，忽略声传播过程中的非线性效应，将声场近似为线性声场。

更进一步地，步骤二涉及的瑞利积分公式为：

$p\left(M\right)=-\frac{\mathrm{jk}}{2\mathrm{\π}}{\∫\∫}_{\mathrm{Source}}{p}_0\left(N\right)\frac{\exp (-\mathrm{jkr})}{r}\mathrm{dS}$

式中，p(M)为声场空间任意一点M的声压，p₀(N)为球壳换能器表面任意一点N处声压；k＝ω₀/c₀为波数，其中，ω₀＝2πf₀为声波的角频率，f₀为声波的频率，c₀为声速；r为M、N两点间距离，Source表示球壳换能器表面区域。

更进一步地，步骤二要求计算获得的球壳换能器轴向及其焦平面径向声压分布曲线与步骤一测量结果的吻合程度为：测量结果与理论计算结果的主瓣偏差低于5%。

更进一步地，步骤三在椭球坐标系下计算声场分布，椭球坐标系中一点(σ,η,φ)与直角坐标系(x,y,z)的变换关系为：

z＝bση

式中，-∞<σ<+∞，0≤η≤1，2b为椭球焦距；为简化声传播方程，作变换η=cosθ,0≤θ≤π/2；

所述的非线性声传播模型为SBE模型，即：

$\begin{array}{c}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_s\&PartialD;\mathrm{\&sigma;}}+\frac{1}{2}\frac{\sin 2\mathrm{\&theta;}}{\mathrm{\&sigma;}(1+{\mathrm{\&sigma;}}^2)}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_s\&PartialD;\mathrm{\&theta;}}+\frac{\mathrm{\&epsiv;}\sqrt{{\mathrm{\&sigma;}}^2+{\sin }^2\mathrm{\&theta;}}}{\mathrm{\&sigma;}(1+{\mathrm{\&sigma;}}^2)}(\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&theta;}}^2}+\cot \mathrm{\&theta;}\frac{\&PartialD;P}{\&PartialD;\mathrm{\&theta;}})+\frac{E}{\mathrm{\&sigma;}}\frac{\&PartialD;P}{\&PartialD;{\mathrm{\&tau;}}_s}\\ =-\frac{\sqrt{{\mathrm{\&sigma;}}^2+{\sin }^2\mathrm{\&theta;}}}{\mathrm{\&sigma;}}(\mathrm{\&alpha;b}\frac{{\&PartialD;}^{3}P}{\&PartialD;{\mathrm{\&tau;}}_s^3}+\frac{b}{2l_D}{\frac{{\&PartialD;}^2{P}^2}{\&PartialD;{\mathrm{\&tau;}}_s^2}}^{})E,(\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_0<0)\end{array}$

$\begin{array}{c}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_p\&PartialD;\mathrm{\&sigma;}}-\frac{\mathrm{\&sigma;}\sin 2\mathrm{\&theta;}}{1+{\mathrm{\&sigma;}}^2}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_p\&PartialD;\mathrm{\&theta;}}-\frac{\mathrm{\&epsiv;}(2-\cos \mathrm{\&theta;})}{1+{\mathrm{\&sigma;}}^2}(\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&theta;}}^2}+\cot \mathrm{\&theta;}\frac{\&PartialD;P}{\&PartialD;\mathrm{\&theta;}})\\ =(\mathrm{\&alpha;b}\frac{{\&PartialD;}^{3}P}{\&PartialD;{\mathrm{\&tau;}}_p^3}+\frac{b}{2l_D}{\frac{{\&PartialD;}^2{P}^2}{\&PartialD;{\mathrm{\&tau;}}_p^2}}^{})E,(\mathrm{\&sigma;}\&GreaterEqual;{\mathrm{\&sigma;}}_0)\end{array}$

其中，σ＜σ₀＜0部分为声场靠近声源表面的球面波近似部分，σ≥σ₀部分为声场远离声源表面的平面波近似部分，σ₀为球面波区域与平面波区域分界面；式中，P=p/P₀为归一化声压，p为点处声压，P₀为球壳换能器表面声压；τ_s＝ω₀[t+b(σ²+sin2θ)^1/2/c₀]为球面波弛豫时间，τ_p＝ω₀[t-bσcosθ/c₀]为平面波弛豫时间，t为球壳换能器激发声波在传播介质中的传播时间；E＝(σ²+cos2θ)/(1+σ²)；α为声传播介质声衰减系数；l_D＝ρ₀c₀ ³/(βP₀ω₀)为冲击波形成距离，ρ₀为声传播介质密度，β为声传播介质非线性系数，ρ₀，α，β均为常数。

更进一步地，步骤三所用的SBE模型的边界条件为：

$P(\mathrm{\&sigma;}=-{\mathrm{\&sigma;}}_{\max },\mathrm{\&theta;}{\mathrm{\&tau;}}_s)=\left\{\begin{array}{cc}\sin {\mathrm{\&tau;}}_s& \mathrm{\&theta;}\&le;{\mathrm{\&theta;}}_0\\ 0& \mathrm{\&theta;}>{\mathrm{\&theta;}}_0\end{array}\right.$

式中，σ_max=d/b， ${\mathrm{\&theta;}}_0={\tan }^{-1}\left(\frac{\tan {\mathrm{\&alpha;}}_0}{\sqrt{1+1/{\mathrm{\&sigma;}}_{\max }^2}}\right),{\mathrm{\&alpha;}}_0={\sin }^{-1}\left(\frac{a}{d}\right).$

本发明的一种测量强聚焦超声非线性声场分布的系统，包括示波器、计算机、信号发生器、功率放大器、三维运动平台、水听器、球壳换能器和水槽，所述的示波器、计算机、信号发生器、功率放大器和球壳换能器依次相连，功率放大器将信号发生器的输出信号放大，驱动球壳换能器激发声场；所述的球壳换能器固定于水槽内部，水槽内盛装有水；所述的三维运动平台控制水听器移动，测量球壳换能器轴向及该球壳换能器焦平面径向声压分布；所述的示波器采集水听器的输出信号；所述的计算机用于调控整个测量系统，计算机完成的功能包括：控制信号发生器输出，控制三维运动平台运动以及存储数据。

3.有益效果

采用本发明提供的技术方案，与已有的公知技术相比，具有如下显著效果：

（1）本发明的一种测量强聚焦超声非线性声场分布的方法，通过引入焦点处二次谐波与基波比值R₂₁反映换能器的边界条件，再根据非线性声传播模型，得到对应激励下的声场分布；由于仅需要测量二次谐波和基波，有效的解决了强非线性条件下，测量系统带宽不够的问题；此外由于R₂₁为相对值，测量过程中对水听器的灵敏度要求较低，避免了因水听器灵敏度低导致检测结果不准确的隐患；

（2）本发明的一种测量强聚焦超声非线性声场分布的系统，结构简单、便于控制，测量准确度高、容易实施。

附图说明

图1是本发明中强聚焦超声非线性声场分布的测量流程图；

图2是本发明的一种测量强聚焦超声非线性声场分布的系统结构示意图；

图3是本发明进行实例测量使用的强聚焦球壳换能器的示意图；

图4是本发明通过测量得到的球壳换能器线性声场轴向分布和焦平面径向分布以及数值仿真结果图；

图5是本发明根据SBE模型计算得到的R₂₁与球壳换能器表面声压P₀的关系曲线图；

图6是本发明球壳换能器表面声压与信号发生器输出信号的关系曲线图；

图7中的（a）～（f）分别是本发明在三种激励下测量得到的以及使用本发明的测量方法得到的焦点波形、相应频谱对比图；

图8中的（a）～（c）分别是本发明在三种激励下球壳换能器焦平面上径向声压分布图。

示意图中的标号说明：

1、示波器；2、计算机；3、信号发生器；4、功率放大器；5、三维运动平台；6、水听器；7、球壳换能器；8、水槽。

具体实施方式

为进一步了解本发明的内容，结合附图和实施例对本发明作详细描述。

实施例1

结合图2，本实施例测量强聚焦超声非线性声场分布的系统，包括示波器（Agilent54830B）1、计算机2、信号发生器（Agilent33250A）3、功率放大器（A150，ENI，Rochester）4、三维运动平台（NewportESP7000）5、水听器（HNA-0400，ONDA或者HFO-660，ONDA）6、球壳换能器7和水槽8，所述的示波器1、计算机2、信号发生器3、功率放大器4和球壳换能器7依次相连，功率放大器4将信号发生器3的输出信号放大，驱动球壳换能器7激发声场；所述的球壳换能器7固定于水槽8内部，水槽8内盛装有水；所述的三维运动平台5控制水听器6在三维空间（x，y，z）内移动，测量球壳换能器7轴向及该球壳换能器7焦平面径向声压分布。所述的示波器1采集水听器6的输出信号。所述的计算机2用于调控整个测量系统，完成功能包括：控制信号发生器3输出，控制三维运动平台5运动以及存储数据。

本实施例的一种测量强聚焦超声非线性声场分布的方法（测量流程参见图1），其测量步骤为：

步骤一、给球壳换能器7施加低能量的激励，使得焦点处二次谐波幅度比基波幅度至少低40dB，此时高次谐波幅度远小于基波幅度，忽略声传播中的非线性效应，将声场看作线性声场。将球壳换能器7固定在水槽8中，控制三维运动平台5带动水听器6逐点测量球壳换能器7轴向以及该球壳换能器7焦平面（焦平面即指垂直于球壳换能器7对称轴且经过焦点的平面）径向的声压，获得球壳换能器7轴向及其焦平面径向的声压分布曲线。

步骤二、根据瑞利积分计算球壳换能器7轴向及该球壳换能器7焦平面径向的声压分布曲线，涉及的瑞利积分公式为：

$p\left(M\right)=-\frac{\mathrm{jk}}{2\mathrm{\&pi;}}{\&Integral;\&Integral;}_{\mathrm{Source}}{p}_0\left(N\right)\frac{\exp (-\mathrm{jkr})}{r}\mathrm{dS}$

式中，j为单位虚数，p(M)为声场空间任意一点M的声压，p₀(N)为球壳换能器7表面任意一点N处声压；k＝ω₀/c₀为波数，其中，ω₀＝2πf₀为声波的角频率，f₀为声波的频率，c₀为声速；r为M、N两点间距离，dS为球壳换能器7表面面积元，Source表示球壳换能器7表面区域，该表面区域级指球壳换能器7用于电-声转化的的压电陶瓷部分。

将计算结果与步骤一获得的测量结果相比较，一般情况下，由于球壳换能器7封装等原因，球壳换能器7表面振动并不是理想的同相位同振幅振动，这导致根据球壳换能器7标称参数计算得到的声场与实际测量声场有差异。微调球壳换能器7的张口半径a和几何焦距d（参看图3），使得计算获得的球壳换能器7轴向及其焦平面径向声压分布曲线与步骤一的测量结果相吻合，吻合程度要求主瓣的偏差低于5%，此时的球壳换能器7参数a和d为球壳换能器7的有效参数。

结合图3，本实施例中球壳换能器7的标称几何尺寸张口半径a为5.0cm，几何焦距d为10.0cm。球壳换能器7中心频率为1.1MHz。利用本实施例球壳换能器7有效参数的确定方法，信号发生器3输出10mv信号，经功率放大器4放大后驱动球壳换能器7，此时由于驱动信号很小，球壳换能器7激发的声场满足线性声场条件，测量球壳换能器7轴向和焦平面径向声压分布。根据球壳换能器7标称参数计算轴向及焦平面径向声场分布，与测量结果比较，发现焦点位置与测量结果不符合，调整球壳换能器7几何焦距为10.3cm，再次计算轴向和焦平面径向声压分布，此时计算结果与实验结果吻合很好。因此球壳换能器7有效参数为：张口半径5.0cm，几何焦距10.3cm。如图4所示，实际测量的球壳换能器7轴向声压和焦平面径向声压分布与计算结果吻合很好。

步骤三、将步骤二获得的有效参数代入非线性声传播模型，计算焦点的二次谐波与基波的比值R₂₁，再改变球壳换能器7表面声压幅值P₀，计算不同表面声压幅值P₀对应的R₂₁，得R₂₁与P₀的关系曲线。具体为：对于本实施例使用的大张角强聚焦球壳换能器7，在椭球坐标系下计算声场分布。椭球坐标系为以正交的椭球和双曲面为坐标轴的坐标系，椭球坐标系中一点(σ,η,φ)与直角坐标系(x,y,z)的变换关系为：

z＝bση

式中，-∞<σ<+∞，0≤η≤1，σ代表椭球面坐标，2b为椭球焦距，η代表双曲面坐标，代表原点与点(x,y,z)的连线在x-y平面的投影与x轴的夹角；为简化声传播方程，作变换η=cosθ,0≤θ≤π/2。大张角强聚焦球壳换能器7的声场将分为两个部分：靠近球壳换能器7表面的球面波近似部分和远离球壳换能器7表面的平面波近似部分。椭球坐标系下非线性声传播方程，即Spheroidalbeamequation(SBE)模型表示为：

$\begin{array}{c}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_s\&PartialD;\mathrm{\&sigma;}}+\frac{1}{2}\frac{\sin 2\mathrm{\&theta;}}{\mathrm{\&sigma;}(1+{\mathrm{\&sigma;}}^2)}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_s\&PartialD;\mathrm{\&theta;}}+\frac{\mathrm{\&epsiv;}\sqrt{{\mathrm{\&sigma;}}^2+{\sin }^2\mathrm{\&theta;}}}{\mathrm{\&sigma;}(1+{\mathrm{\&sigma;}}^2)}(\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&theta;}}^2}+\cot \mathrm{\&theta;}\frac{\&PartialD;P}{\&PartialD;\mathrm{\&theta;}})+\frac{E}{\mathrm{\&sigma;}}\frac{\&PartialD;P}{\&PartialD;{\mathrm{\&tau;}}_s}\\ =-\frac{\sqrt{{\mathrm{\&sigma;}}^2+{\sin }^2\mathrm{\&theta;}}}{\mathrm{\&sigma;}}(\mathrm{\&alpha;b}\frac{{\&PartialD;}^{3}P}{\&PartialD;{\mathrm{\&tau;}}_s^3}+\frac{b}{2l_D}{\frac{{\&PartialD;}^2{P}^2}{\&PartialD;{\mathrm{\&tau;}}_s^2}}^{})E,(\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_0<0)\end{array}$

$\begin{array}{c}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_p\&PartialD;\mathrm{\&sigma;}}-\frac{\mathrm{\&sigma;}\sin 2\mathrm{\&theta;}}{1+{\mathrm{\&sigma;}}^2}\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&tau;}}_p\&PartialD;\mathrm{\&theta;}}-\frac{\mathrm{\&epsiv;}(2-\cos \mathrm{\&theta;})}{1+{\mathrm{\&sigma;}}^2}(\frac{{\&PartialD;}^{2}P}{\&PartialD;{\mathrm{\&theta;}}^2}+\cot \mathrm{\&theta;}\frac{\&PartialD;P}{\&PartialD;\mathrm{\&theta;}})\\ =(\mathrm{\&alpha;b}\frac{{\&PartialD;}^{3}P}{\&PartialD;{\mathrm{\&tau;}}_p^3}+\frac{b}{2l_D}{\frac{{\&PartialD;}^2{P}^2}{\&PartialD;{\mathrm{\&tau;}}_p^2}}^{})E,(\mathrm{\&sigma;}\&GreaterEqual;{\mathrm{\&sigma;}}_0)\end{array}$

上式中，声场分为靠近球壳换能器7表面（即σ＜σ₀＜0）的球面波近似部分，和远离球壳换能器7表面（即σ≥σ₀）的平面波近似部分，σ0为球面波区域与平面波区域分界面。对于轴对称的球壳换能器7，声场分布与无关；P=p/P₀为归一化声压，p为点处声压，P0为球壳换能器7表面声压；τ_s＝ω₀[t+b(σ²+sin2θ)^1/2/c₀]为球面波弛豫时间，τ_p＝ω₀[t-bσcosθ/c₀]为平面波弛豫时间，ω0为声波的角频率，t为球壳换能器激发声波在传播介质中的传播时间，c₀为声传播速度；ε＝1/(2kb)，k=ω₀/c₀为波数；E＝(σ²+cos2θ)/(1+σ²)；α为声传播介质的声衰减系数；l_D＝ρ₀c₀ ³/(βP₀ω₀)为冲击波形成距离，ρ₀为声传播介质的密度，β为声传播介质的非线性系数。ρ₀，c₀，α，β为声传播介质的常数，对于本实施例使用的纯水，ρ₀=1000kg/m³，c₀=1495m/s，α=0.025Np/m/MHz²，β=3.5。

SBE模型的边界条件为：

$P(\mathrm{\&sigma;}=-{\mathrm{\&sigma;}}_{\max },\mathrm{\&theta;}{\mathrm{\&tau;}}_s)=\left\{\begin{array}{cc}\sin {\mathrm{\&tau;}}_s& \mathrm{\&theta;}\&le;{\mathrm{\&theta;}}_0\\ 0& \mathrm{\&theta;}>{\mathrm{\&theta;}}_0\end{array}\right.;$

式中，σ_max=d/b，d为球壳换能器7几何焦距； ${\mathrm{\&theta;}}_0={\tan }^{-1}\left(\frac{\tan {\mathrm{\&alpha;}}_0}{\sqrt{1+1/{\mathrm{\&sigma;}}_{\max }^2}}\right),{\mathrm{\&alpha;}}_0={\sin }^{-1}\left(\frac{a}{d}\right),$ a为球壳换能器7张口半径。

将步骤二得到的球壳换能器7有效参数代入边界条件，数值求解SBE模型。本实施例通过有限差分算法计算SBE模型，具体过程为：第一步、通过傅里叶级数展开，将时域声压P变换到频域其中C_n为n次谐波幅度；第二步、将傅里叶级数展开式代入SBE模型，可以得到关于C_n的抛物线型偏微分方程；第三步、对于各次谐波，利用隐式后向差分算法，将偏微分方程离散化，再进行求解得到C_n。由于声场的该计算方法已经较为成熟，其具体计算过程就不再赘述，可参见相关文献。计算焦点处二次谐波与基波的比值R₂₁：

R₂₁＝C₂/C₁；

改变球壳换能器7表面声压幅值P₀，计算不同球壳换能器7表面声压幅值P₀对应的R₂₁，得到R₂₁与P₀的曲线。高强度聚焦超声的传播非线性与球壳换能器7表面声压幅度、传播距离以及传播介质的声学参数有关，当确定了传播介质，在焦点的波形非线性畸变程度（用焦点处的二次谐波C2与基波C1衡量）仅仅与球壳换能器7表面声压幅度有关，且表面声压幅度越强，焦点波形的非线性畸变越明显，即R₂₁越大，因此R₂₁随P₀增加而单调增加（参见图5）。

步骤四、对球壳换能器7施加任意幅度的单频电信号激励，激发声场，利用光纤水听器测量焦点波形p(t)，对该焦点波形进行傅里叶变换：

$S\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\&infin;}^{\&infin;}p\left(t\right)e^{-\mathrm{j\&omega;t}}\mathrm{dt}$

计算得R₂₁′：

R₂₁′＝S(2ω₀)/S(ω₀)

利用焦点二次谐波与基波比值R₂₁与球壳换能器7表面声压P₀的单调关系，测定R₂₁即可获得P₀；根据步骤三得到的R₂₁与P₀的关系曲线进行插值，得到R₂₁′对应的球壳换能器7表面声压P₀′。

步骤五、将步骤四获得的表面声压P₀′带入步骤三所述的SBE模型，计算得到步骤四施加的电信号激励对应的非线性声场分布。如图6所示，信号发生器3输出信号幅值与球壳换能器7表面声压幅值的关系曲线。依据图6，在任意激励下的声场分布均可由SBE模型计算得到。

如图7所示，在表面声压分别为81kPa、259kPa、350kPa时，测量和计算得到的焦点波形，以及各自相应频谱。在表面声压为81kPa和259kPa时，实验测量的波形和理论计算的结果吻合很好；在表面声压为350kPa时，冲击波产生，实验测量结果和计算结果存在较大差异。对焦点波形进行傅里叶变换，比较各次谐波。示波器采样频率为100MHz，分析50MHz内的各次谐波幅值。由图7可见，在表面声压350kPa时，50MHz内的各次谐波幅值吻合的很好，这说明波形的差异主要由于测量系统带宽不够。另外，图7中在表面声压81kPa时，频谱在高于10MHz后，存在明显差异；主要原因是此时高次谐波太小，淹没在系统噪声中。

如图8所示，在表面声压分别为81kPa，259kPa，350kPa时，焦平面径向声压分布。计算结果和实验结果吻合较好。

综上所述，本实施例提供的强聚焦超声非线性声场分布测量方法，通过引入焦点二次谐波与基波比值R₂₁反映球壳换能器7边界条件，再根据非线性声传播模型，得到对应激励下的声场分布。由于仅需要测量二次谐波和基波，有效的解决了强非线性下，测量系统带宽不够的问题；此外由于R₂₁为相对值，测量中不需要使用水听器6的灵敏度，因此避免了水听器6灵敏度不高的隐患；本实施例的测量方法操作简便，实用性强，使用的测量系统结构简单、便于控制，测量准确度高、容易实施。

Claims (6)

1.一种测量强聚焦超声非线性声场分布的方法，其测量步骤为：

步骤一、在线性声场条件下，测量球壳换能器轴向以及该球壳换能器焦平面径向的声压，获得球壳换能器轴向及其焦平面径向的声压分布曲线；

步骤二、根据瑞利积分计算球壳换能器轴向及该球壳换能器焦平面径向的声压分布曲线，将计算结果与步骤一获得的测量结果相比较，调节球壳换能器的张口半径a和几何焦距d，使得计算获得的球壳换能器轴向及其焦平面径向声压分布曲线与步骤一的测量结果相吻合，获得球壳换能器有效参数a和d；

步骤三、将步骤二获得的有效参数代入非线性声传播模型，计算焦点的二次谐波与基波的比值R₂₁，再改变球壳换能器表面声压幅值P₀，计算不同表面声压幅值P₀对应的R₂₁，得R₂₁与P₀的关系曲线；

步骤四、对球壳换能器施加单频电信号激励，测量焦点波形，对该焦点波形进行傅里叶变换，计算得R₂₁′，根据步骤三得到的R₂₁与P₀的关系曲线进行插值，得到R₂₁′对应的球壳换能器表面声压P₀′；

步骤五、将步骤四获得的表面声压P₀′带入步骤三所述的非线性声传播模型，计算得到步骤四施加的电信号激励对应的非线性声场分布。

2.根据权利要求1所述的一种测量强聚焦超声非线性声场分布的方法，其特征在于：步骤一测量声压的条件为：对球壳换能器施加低能量激励，使得球壳换能器焦点处基波幅度与二次谐波幅度的差值不低于40dB，忽略声传播过程中的非线性效应，将声场近似为线性声场。

3.根据权利要求2所述的一种测量强聚焦超声非线性声场分布的方法，其特征在于：步骤二涉及的瑞利积分公式为：

$p\left(M\right)=-\frac{jk}{2\mathrm{\&pi;}}\&Integral;{\&Integral;}_{Source}{p}_0\left(N\right)\frac{\exp (-jkr)}{r}dS$

式中，p(M)为声场空间任意一点M的声压，p₀(N)为球壳换能器表面任意一点N处声压；k＝ω₀/c₀为波数，其中，ω₀＝2πf₀为声波的角频率，f₀为声波的频率，c₀为声速；r为M、N两点间距离，Source表示球壳换能器表面区域。

4.根据权利要求3所述的一种测量强聚焦超声非线性声场分布的方法，其特征在于：步骤二要求计算获得的球壳换能器轴向及其焦平面径向声压分布曲线与步骤一测量结果的吻合程度为：测量结果与理论计算结果的主瓣偏差低于5％。

5.根据权利要求2或3所述的一种测量强聚焦超声非线性声场分布的方法，其特征在于：步骤三在椭球坐标系下计算声场分布，椭球坐标系中一点与直角坐标系(x,y,z)的变换关系为：

z＝bση

式中，-∞<σ<+∞，0≤η≤1，2b为椭球焦距；为简化声传播方程，作变换η＝cosθ,0≤θ≤π/2；

所述的非线性声传播模型为SBE模型，即：

$\begin{array}{c}\frac{{\&part;}^{2}P}{\&part;{\mathrm{\&tau;}}_s\&part;\mathrm{\&sigma;}}+\frac{1}{2}\frac{\sin 2\mathrm{\&theta;}}{\mathrm{\&sigma;}(1+{\mathrm{\&sigma;}}^2)}\frac{{\&part;}^{2}P}{\&part;{\mathrm{\&tau;}}_s\&part;\mathrm{\&theta;}}-\frac{\mathrm{\&epsiv;}\sqrt{{\mathrm{\&sigma;}}^2+{\sin }^2\mathrm{\&theta;}}}{\mathrm{\&sigma;}(1+{\mathrm{\&sigma;}}^2)}(\frac{{\&part;}^{2}P}{\&part;{\mathrm{\&theta;}}^2}+\cot \mathrm{\&theta;}\frac{\&part;P}{\&part;\mathrm{\&theta;}})+\frac{E}{\mathrm{\&sigma;}}\frac{\&part;P}{\&part;{\mathrm{\&tau;}}_s}\\ \begin{array}{cc}=-\frac{\sqrt{{\mathrm{\&sigma;}}^2+{\sin }^2\mathrm{\&theta;}}}{\mathrm{\&sigma;}}(\mathrm{\&alpha;}b\frac{{\&part;}^{3}P}{\&part;{\mathrm{\&tau;}}_s^3}+\frac{b}{2l_D}\frac{{\&part;}^2{P}^2}{\&part;{\mathrm{\&tau;}}_s^2})E& (\mathrm{\&sigma;}<{\mathrm{\&sigma;}}_0<0)\end{array}\end{array}$

$\begin{array}{c}\frac{{\&part;}^{2}P}{\&part;{\mathrm{\&tau;}}_p\&part;\mathrm{\&sigma;}}-\frac{\mathrm{\&sigma;}\sin \mathrm{\&theta;}}{1+{\mathrm{\&sigma;}}^2}\frac{{\&part;}^{2}P}{\&part;{\mathrm{\&tau;}}_p\&part;\mathrm{\&theta;}}-\frac{\mathrm{\&epsiv;}(2-cos\mathrm{\&theta;})}{1+{\mathrm{\&sigma;}}^2}(\frac{{\&part;}^{2}P}{\&part;{\mathrm{\&theta;}}^2}+\cot \mathrm{\&theta;}\frac{\&part;P}{\&part;\mathrm{\&theta;}})\\ \begin{array}{cc}=(\mathrm{\&alpha;}b\frac{{\&part;}^{3}P}{\&part;{\mathrm{\&tau;}}_p^3}+\frac{b}{2l_D}\frac{{\&part;}^2{P}^2}{\&part;{\mathrm{\&tau;}}_p^2})E& (\mathrm{\&sigma;}\&GreaterEqual;{\mathrm{\&sigma;}}_0)\end{array}\end{array}$

其中，σ<σ₀<0部分为声场靠近声源表面的球面波近似部分，σ≥σ₀部分为声场远离声源表面的平面波近似部分，σ₀为球面波区域与平面波区域分界面；式中，P＝p/P₀为归一化声压，p为点处声压，P₀为球壳换能器表面声压；τ_s＝ω₀[t+b(σ²+sin2θ)^1/2/c₀]为球面波弛豫时间，τ_p＝ω₀[t-bσcosθ/c₀]为平面波弛豫时间，t为球壳换能器激发声波在传播介质中的传播时间；E＝(σ²+cos2θ)/(1+σ²)；α为声传播介质声衰减系数；l_D＝ρ₀c₀ ³/(βP₀ω₀)为冲击波形成距离，ρ₀为声传播介质密度，β为声传播介质非线性系数，ρ₀，α，β均为常数。

6.根据权利要求5所述的一种测量强聚焦超声非线性声场分布的方法，其特征在于：步骤三所用的SBE模型的边界条件为：

$P(\mathrm{\&sigma;}=-{\mathrm{\&sigma;}}_{max},\mathrm{\&theta;},{\mathrm{\&tau;}}_s)=\left\{\begin{array}{cc}{\mathrm{sin\&tau;}}_s& \mathrm{\&theta;}\&le;{\mathrm{\&theta;}}_0\\ 0& \mathrm{\&theta;}>{\mathrm{\&theta;}}_0\end{array}\right.$

式中，σ_max＝d/b， ${\mathrm{\&theta;}}_0={\tan }^{-1}\left(\frac{{\mathrm{tan\&alpha;}}_0}{\sqrt{1+1/{\mathrm{\&sigma;}}_{max}^2}}\right),$ ${\mathrm{\&alpha;}}_0={\sin }^{-1}\left(\frac{a}{d}\right).$