在医用超声生理组织成像、工业材料和设备无损评估、水声参量阵等领域，含气泡介质都有着广泛的应用。其中，气泡尺寸分布的声学测量与反演是获取气泡群物理参数，进而对应用效果进行定量化分析评估的重要环节。目前，气泡群尺寸分布参数的声学测量与反演方法按测量方式分为被动测量法和主动测量法。Leighton^[1]提出被动测量法，即通过测量海水中破碎波产生的气泡向外辐射的声压信号，再通过时频分析获得气泡半径分布及对应的气泡数量。该方法受信噪比的影响较大，低信噪比时经常无法检测到信号，导致得到的分布结果不够准确。主动测量法利用发射声信号，经过气泡层后，接收到的声学数据，再利用物理模型进行反演，进而得到气泡的分布情况。其中测量方法主要是基于气泡共振响应获取相关参数，如Terril^[2-3]和Matthieu Cavaro^[4]通过测量声速进而反演气泡参数，Medwin等^[5-9]提出的脉冲回声法和连续波法，Vagle等^[10]提出的回声测深方法结合声逆散射方法等。此外还有基于衰减量反演气泡分布的经典方法，如Clay和Medwin^[11-12]共同提出的共振估计理论，其假设外部激励仅受与其发生频率共振的半径的气泡的影响，不考虑液体的粘滞和热传导并认为阻尼系数为一个常数，进而化简了声衰减的积分公式。Commander^[13-14]提出的有限元法化简了声衰减积分公式。Caruthers^[15]提出的迭代法对共振估计法进行模拟和修正等。在气泡振动与声学特性的理论研究方面，Minnaert^[16]忽略了气泡表面张力和介质中粘滞阻力的影响，最早推导出气泡的共振频率公式。国内的钱祖文等^[17-18]对人工产生的气幕及海洋中的气泡进行了大量研究，运用声学测量手段对气泡幕和尾迹气泡群的尺寸分布和气体体积分数进行了反演，并建立了相关的参数模型。综合已有的研究方法可以看出，用声学方法研究气泡分布主要依托于含气泡介质中气泡群的声学参数，即声散射特性，如气泡对声速频散、声衰减、声散射截面^[19]等的影响特性。将气泡对水介质物理特性影响的规律加以修正，即可用来粗略地反演含气泡海水介质中相关的气泡参数。线性方法基于共振估计法，计算线性散射截面或衰减截面进而反演气泡相关参数^[20-21]，对共振态气泡较为敏感，而对非共振态气泡分析的精度不高，常会出现频率模糊问题，使分析结果出现偏差甚至错误。非共振气泡对声散射特性的影响程度、进而对气泡尺寸分布反演结果的影响程度，是本文讨论的重点。

本文综合研究和分析非共振气泡对含气泡介质中的声频散、声衰减和声散射等声学特性的影响，指出仅考虑共振气泡影响是传统气泡分布参数线性声学反演方法的误差来源；建立气泡群声散射数理模型，利用B样条差值求解不适定的逆散射问题方程，获得了比共振反演方法更准确的气泡分布参数反演结果，抑制了经典线性反演方法中的伪峰问题。

1 非共振气泡对含气泡介质中声散射特性的影响

根据Kargl^[20]的研究，含气泡介质中的等效波数k_m满足:

$ {k_m} = \frac{\omega }{c} - i\alpha $

(1)

式中：c为声速; ω为声激励角频率；α为衰减。且根据Medwin的经典散射理论：

$ \frac{{{\sigma _e}}}{{{\sigma _s}}} = \frac{\delta }{{ka}} $

(2)

式中：σ_e为气泡的衰减截面；σ_s为散射截面；δ为阻尼系数；k为含气泡介质中的等效波数；a为气泡半径。从式(2)可以看出，通过研究非共振气泡与气泡群的等效波数和散射截面的关系，能够分析非共振气泡对气泡群声散射特性的影响。

1.1 非共振气泡对声频散、声衰减及等效波数的影响

Van Wijingaarden最先提出了以混合介质为基础的模型理论，导出了混合介质的线性声波动方程。Commander和Prosperetti在Van Wijingaarden混合介质模型的基础上推导出了一个等效波数的最低阶近似表达式，并讨论了混合介质内的声速和声衰减等声特性。在不考虑多重散射的情况下，假设声波在含气泡海水介质中线性传播，气泡在声波激励下作周期振荡，得到的介质的等效波数k_m为：

$ k_m^2 = {k^2} + 4\pi {\omega ^2}\int\limits_0^\infty {\frac{{an(a){\rm{d}}a}}{{\omega _0^2 - {\omega ^2} + ib}}} $

(3)

式中ω为角频率，共振频率和阻尼表达式为：

$ \omega _0^2 = \frac{{{p_0}}}{{\rho {a^2}}}\left( {{\rm{Re}} (\phi ) - \frac{{2\sigma }}{{a{p_0}}}} \right) $

(4)

$ b = {b_v} + {b_t} + {b_a} = \frac{{2\mu }}{{\rho {a^2}}} + \frac{{{p_0}}}{{2\rho {a^2}\omega }} {\rm{Im}} (\phi ) + \frac{{{\omega ^2}a}}{{2c}} $

(5)

阻尼b包括粘滞阻尼b_v、热阻尼b_t和再辐射b_a。其中参数解释和推导过程见参考文献[22]，此处不再赘述。根据等效波数可以得到特定分布下的声速和衰减参数。

1.2 非共振气泡对气泡的声散射截面的影响

气泡散射模型的研究主要分为3个体系，分别为Wildt体系、Devin体系以及非线性体系。这3种体系的假设条件和推导依据有所差异，因此对非共振气泡的影响逐一进行讨论。

1.2.1 Wildt理论

Wildt通过精确地表示入射波p_i(幅度为A的平面波)和散射波p_s(幅度为B/r的球面波)，再通过散射功率4πr²|p_s|²/(2ρ_liqc)与入射强度|p_i|²/(2ρ_liqc)的比值，其中r为声压与气泡中心的距离，ρ_liq为水的密度，得到散射截面与散射波的归一化振幅|B/A|的关系式为：

$ {\sigma _s} = 4\pi |B/A{|^2} = \frac{{4\pi R_0^2}}{{{{(\omega _{{\rm{res}}}^2/{\omega ^2} - 1)}^2} + \delta _{{\rm{Wildt}}}^2}} $

(6)

式中：ω_res为气泡共振频率，满足如下表达式：

$ \gamma \frac{{{\omega ^2}}}{{\omega _{res}^2}} = \left( {1 + \frac{{{R_{{\rm{ Laplace }}}}}}{{{R_0}}}} \right) {\rm{Re}} \varGamma (\omega ) - \frac{{{R_{{\rm{ Laplace }}}}}}{{3{R_0}}} $

(7)

式中：γ为比热比；R₀为气泡平衡半径；R_Laplace为拉普拉斯半径；Γ是复多变指数。

1.2.2 Devin理论

Devin理论在Wildt的基础上，根据将气泡径向位移与强迫压力相关联的线性微分方程。给出了刚度参数K近似为共振角频率平方ω₀²的函数。Medwin将Devin的阻尼因子表达式代入式(6)，得到了Devin体系的散射截面为：

$ {{\sigma _s} = \frac{{4\pi R_0^2}}{{{{(\omega _{res}^2/{\omega ^2} - 1)}^2} + {\delta ^2}}}} $

(8)

且:

$ {{\sigma _e} = {\sigma _s}\frac{\delta }{\varepsilon }} $

(9)

其中：

$ \delta = \frac{{4{\eta _s}}}{{{\rho _{liq}}R_0^2\omega }} + \frac{{3(\gamma - 1)}}{X}\zeta (X)\frac{K}{{{\omega ^2}}} + \varepsilon $

(10)

式中：ε₀是无量纲气泡半径，定义为ε₀≡ε(ω₀)=ω₀R₀/c; X为扩散比系数。

1.2.3 非线性理论

Church在以上2种理论的基础上，考虑了气泡的高阶振动，且仅考虑粘滞阻尼系数δ_vis，通过对非线性方程Rayleigh-Plesset进行线性化，得到简化的自由气泡的散射截面表达式：

$ {\sigma _s} = \frac{{4\pi R_0^2}}{{{{(\omega _o^2/{\omega ^2} - 1)}^2} + {{(2{\delta _{vis}}/\omega )}^2}}} $

(11)

式中ω₀为与入射声频率相关的函数：

$ \omega _0^2 = 3 {\rm{Re}} \varGamma \frac{{{P_{{\rm{gas}}}}}}{{{\rho _{{\rm{liq}}}}R_0^2}} - \frac{\tau }{{{\rho _{{\rm{liq}}}}R_0^3}} $

(12)

式中P_gas算是气泡内部的声压。

1.2.4 非共振气泡对散射截面的影响分析

对于服从特定分布的气泡群，假设气泡群孔隙率较低，不考虑气泡间的粘连和相互作用，可以通过积分∫σ_sN(a)da的形式求得气泡群的总散射截面；其中σ_s为单个气泡的散射截面，气泡尺寸分布函数N(a)为单位体积的含气泡介质中，半径为a的气泡数量。从上文总结的3种散射截面理论，可以看出单个气泡的散射截面均与激励频率和气泡半径有关。当忽略非共振气泡对散射截面的影响时，也可以理解为，仅当激励频率为某半径气泡的共振频率时，该半径的气泡才对气泡群的总散射截面有贡献。

1.3 非共振对声散射特性影响的仿真与分析

通过综合在不同气泡分布下，是否考虑非共振气泡影响下的声速频散、声衰减和声散射截面，就可得到非共振气泡对声学特性的响应。

本文列举4种假设分布的情况。4种分布的主峰都设在100 μm处，峰值处的气泡数量为5.5×10⁹。图 1(a)为单一分布(分布函数为δ函数)，其余3种情况为高斯分布，标准差σ分别为0.15(图 2(a))、0.5(图 3(a))和1(图 4(a))。

在气泡单一分布情况时(图 1)，不存在非共振气泡，因此是否考虑非共振气泡影响的声速曲线和衰减曲线基本重合(图 1(b)、1(c))。图 2中，在激励为100 μm气泡的共振频率时，非共振气泡比单一分布时增加，声速曲线(图 2(b))发生畸变。衰减曲线(图 2(c))普遍升高，峰值向大于共振频率方向偏移。对于图 3和图 4的情况，随着非共振气泡进一步增多，声速曲线(图 3(b)、4(b))畸变更加明显，声衰减曲线(图 3(c)和4(c))与考虑共振气泡的声衰减曲线差异越来越大。

声速频散和衰减的变化可以用来综合反映非共振气泡对等效波数的影响程度。可见随着非共振气泡比重的增大，对声速频散和衰减的影响越来越大。在实际的海洋环境中，气泡的尺寸分布具有较宽的尺度分布范围。这说明当使用声速和衰减对气泡尺寸进行反演时，如果忽略了非共振气泡的影响一定会带来较大的误差。

在图 1(d)、2(d)、3(d)、4(d)中，呈现了考虑非共振气泡后在3种体系下的散射截面对比，分别描述了前述的Wildt、Devin和非线性理论中是否考虑非共振气泡的散射截面结果，实线代表忽略非共振的散射截面仿真结果，虚线代表考虑非共振的结果。由图可见3种散射截面体系下的规律是基本一致的。当气泡群为单一分布，不存在非共振气泡，因此图 1(d)中是否考虑非共振气泡影响的总散射截面曲线形状吻合位置接近。当服从图 2(a)中的分布时，若激励频率为100 μm气泡的共振频率，非共振气泡增多，散射截面曲线等效于将分布函数进行不同的加权，仅考虑共振气泡时加权较为均匀，体现为图 2(d)中3条实线基本重合。在考虑非共振气泡(虚线)时，出现了一些的明显峰值，且后者曲线的峰值向大于共振频率方向偏移。对于图 3(d)和图 4(d)，随着非共振气泡进一步增多(尺寸分布函数标准差变大)，是否考虑非共振气泡对散射截的分析的影响差异则更加明显，可以认为非共振气泡对散射截面的影响越来越大。

由以上分析可知，非共振气泡对等效波数和散射截面的影响是不能被忽略的，尤其是当等效波数和散射截面作为反演的目标函数时，对反演结果的准确程度尤其重要。综上说明传统的共振估计法在物理原理上是有较大缺陷的，需要探求可充分考虑非共振气泡对声散射特性影响的反演方法。

2 考虑非共振气泡影响的气泡群尺寸分布反演

通过进行非共振气泡对声散射特性的影响分析，并为了能基于声学手段获得准确的气泡尺寸分布，必须考虑非共振对散射强度和衰减的影响。这样做的过程中就需要求解不适定的逆散射问题方程(核函数呈病态的第1类Fredholm积分方程):

$ \alpha (f) = \int_0^\infty K (f,a)N(a){\rm{d}}a $

(13)

式中：α(f)是单位体积气泡群对不同入射频率f引起的衰减；K(f, a)是核函数；N(a)是待求的单位体积的含气泡水介质中半径为a的数量(即气泡群尺寸分布函数)。为方便数值形式求解，将方程(13)进行离散化：

$ {\alpha _i} = \sum\limits_j {{K_{ij}}} {N_j},i = 1,2, \cdots ,m,j = 1,2, \cdots ,n $

(14)

方程以矩阵形式可写为：

$ \left[ {\begin{array}{*{20}{c}} {{\alpha _1}}\\ {{\alpha _2}}\\ {{\alpha _3}}\\ \vdots \\ {{\alpha _m}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{K_{11}}}&{{K_{12}}}&{{K_{13}}}& \cdots &{{K_{1n}}}\\ {{K_{21}}}&{{K_{22}}}&{{K_{23}}}& \cdots &{{K_{2n}}}\\ {{K_{31}}}&{{K_{32}}}&{{K_{33}}}& \cdots &{{K_{3n}}}\\ \vdots & \vdots & \vdots &{}& \vdots \\ {{K_{m1}}}&{{K_{m2}}}&{{K_{m3}}}& \cdots &{{K_{mn}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{N_1}}\\ {{N_2}}\\ {{N_3}}\\ \vdots \\ {{N_n}} \end{array}} \right] $

(15)

式中：α是m×1矢量；N是n×1矢量；K是m×n矩阵。将积分方程转换为离散和时，为了得到精确的积分结果，积分步长通常较小，因此常遇到的条件是m < n，这对应于一个不适定的问题。对于一个不适定的病态矩阵方程，本文采用B样条插值方法来获得方程的解，它是对正则化方法的一种改进。正则化方法首先由Phillips提出，然后由Twomey改进，它根据所拥有的先验信息对解施加约束。通常使用的一个约束是平滑度约束，首先假设要寻找的解是平滑的，并定义其度量为平滑度。在定义了这个度量之后，就可以从所有可能的解中挑选出最平滑的解。此处的平滑度s(N)可定义为:

$ s(\mathit{\boldsymbol{N}}) = \int_0^\infty {{{\left( {\frac{{{\partial ^2}N(a)}}{{\partial {a^2}}}} \right)}^2}} {\rm{d}}a $

(16)

在离散形式下，可以写为:

$ S(\mathit{\boldsymbol{N}}) = \sum\limits_i {{{({N_{i + 1}} - 2{N_i} + {N_{i - 1}})}^2}} $

(17)

这是一个二次项，矩阵形式为:

$ S(\mathit{\boldsymbol{N}}) = {\mathit{\boldsymbol{N}}^{\rm{T}}}\mathit{\boldsymbol{HN}} $

(18)

式中H为:

$ \mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{r}} 1&{ - 2}&1&0& \cdots &0&0&0\\ { - 2}&5&{ - 4}&1&{}&{}&{}&{}\\ 1&{ - 4}&6&{ - 4}&{}&{}&{}&{}\\ \vdots &{}&{}&{}&{}&{}&{}& \vdots \\ 0&0&0&0& \cdots &1&{ - 2}&1 \end{array}} \right] $

(19)

为实现约束最小化s(N)，现寻找一种能够同时满足N=K^-1α的解，其中K^-1为K的广义逆。利用拉格朗日算子λ，通过确定N(a)，来使总约束l最小化，其表达式为：

$ l = \lambda {\mathit{\boldsymbol{N}}^{\rm{T}}}{\bf{HN}} + {(\mathit{\boldsymbol{\alpha }} - \mathit{\boldsymbol{KN}})^{\rm{T}}}(\mathit{\boldsymbol{\alpha }} - \mathit{\boldsymbol{KN}}) $

(20)

对N(a)的每一个元素进行微分并令其等于零，可得到：

$ \mathit{\boldsymbol{N}} = {({\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} + \lambda \mathit{\boldsymbol{H}})^{ - 1}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{\alpha }} $

(21)

通过选择合适的拉格朗日算子λ的值，能够使得定义|α-KN|²的残差小于预定值。在这种方法中，通过施加平滑约束来克服解的不唯一性。而将矩阵λH添加到K^TK中会导致所得矩阵的特征值增加，从而避免了不稳定性。

对于气泡群的衰减截面，有定义式：

$ \alpha (f) = \int_0^\infty {\frac{{4\pi a\delta n(a)/k}}{{{{[{{({f_a}/f)}^2} - 1]}^2} + {{(ka)}^2}}}} {\rm{d}}a $

(22)

式中f_a为半径为a气泡的共振频率。而若只考虑共振气泡对声衰减的影响，则有：

$ \alpha (f) = \frac{{4\pi a\delta n(a)/k}}{{{\delta _{0r}}}}\int\limits_0^\infty {\frac{{{\rm{d}}a}}{{{{[{{({f_a}/f)}^2} - 1]}^2} + {\delta ^2}}}} $

(23)

式中δ_{0 r}是共振半径的阻尼系数。这样可以得到2个有着不同核函数的第1类Fredholm方程。假设气泡的分布分别服从幂指数分布和高斯分布(最常见的2种近似气泡分布)，首先用高斯数值积分计算得到声衰减作为目标函数，再使用B样条插值的方式进行反演，就可解算出气泡分布。反演的结果对比如图 5所示。

通过对常见气泡分布进行的数值仿真，可以看到，是否考虑非共振气泡对尺寸分布的反演存在着影响，尤其是在大半径处可能会出现本不存在的峰值起伏。在考虑非共振气泡的影响后，这些“伪峰值”的情况就可以得到改善。

为方便对比2种反演效果，可分别选择观测2种反演结果与假设分布数量的最大峰值偏差、次大峰值偏差和第3峰值偏差进行比较；而曲线整体的差异程度则是使用和基准的相关系数来度量。计算得出的结果如表 1。

由表可知，2种假设分布下，考虑了非共振影响的反演结果与假设分布都有着更大的相关系数。峰值偏差也更小。可以得出：在不考虑干扰的前提下，对于常见的气泡分布(幂指数分布和高斯分布见图 6、7)，考虑非共振影响能够得到更好的反演效果。

下面对该方法下的反演过程进行误差干扰分析。在同上假设气泡尺寸分布为幂指数分布和高斯分布的条件下，给衰减量分别加入1%、3%、5%和10%倍原衰减量均值的高斯噪声干扰。再经过100次蒙特卡洛实验后，得到的反演效果对比如下：对比发现，噪声干扰的引入加剧了忽略非共振气泡带来的“伪峰值”，误差影响较大。尤其在大半径位置得到的反演结果已不具备参考价值。而考虑非共振的结果更接近假设的分布，抗干扰能力更好。

干扰条件下的反演效果如表 2。从表中也可以看出，相关系数越接近1，说明反演效果越好。假设幂指数分布时，当引入误差干扰，忽略非共振的反演结果影响很大，而对考虑了非共振的反演结果几乎没有影响。这个结论在假设高斯分布时尤为显著，在较大误差5%、10%时，忽略非共振的反演结果与假设分布的相关系数仅有0.33甚至更小，反演几乎失效，这在图 7(c)(d)中也可以看出。以上这些都说明了在考虑非共振影响时其拥有更好的抗干扰能力。

3 结论

1) 在是否考虑非共振气泡两种情况下，等效波数(声速频散、衰减)和散射衰减截面都有显著的差异。非共振气泡对声散射特性的影响，导致在反演阶段必然会产生较大误差，为后续的讨论提供了正向的理论依据；

2) 以B样条插值的方式对分布参数进行了反演，是否考虑非共振时的反演效果有较大区别；

3) 反演结果表明，考虑非共振气泡时，反演结果能够对“伪峰”进行抑制，具有更好的准确性和抵抗干扰的能力。因此为了保证反演效果，考虑非共振气泡的影响是有必要的。