高斯光束入射椭球形气泡的光散射特性研究

Study on Light Scattering Characteristics of Spheroidal Bubble with Gaussian Light Incidence

摘要: 气泡的光散射研究是颗粒测量的重要组成部分。目前的光散射技术在计算大气泡时耗时较长、效率较低。本文提出了一个球体模型来分析椭球形气泡的光散射特性。用几何光学近似方法模拟了轴向高斯光束照射椭球形气泡时散射光强的角分布。计算中考虑了衍射、散度因子和相移。研究了高斯光束束腰半径对大气泡散射特性的影响。

Abstract: The study of light scattering of bubbles is an important part of particle measurement. Current light scattering techniques are less efficient for calculating large bubbles. In this paper, aspheroidal model is proposed to analyze the light scattering characteristics of bubbles. The angular distribution of scattered light intensity of a spheroidal bubble illuminated by an axial Gaussian beam is simulated by geometric optical approximation. Diffraction, divergence factor and phase shift are considered in the calculation. The influence of the waist radius of a Gaussian beam on the scattering characteristics of a large spheroidal bubble is studied.

1. 介绍

舰船尾流、鱼群游动、海水的冲击等是导致海水产生气泡的原因，而气泡的存在会影响光在水中的传播。而海水中的气泡在湍流的作用下往往不是规则的球形气泡，而是类似于椭球形的气泡，因此研究椭球形气泡的光散射特性具有重要的作用。

在光散射的研究历程中，通过麦克斯韦方程组和边界条件求解出的洛伦兹–米理论(简称米理论)是研究光散射特性的经典理论[1]。但是由于Mie理论的算法中含多种复杂函数组合，计算也十分耗时。在光散射问题求解过程中，研究者们也做出了许多的努力，提出了许多的计算方法，如Debye级数[2]、分离变量法、Watson变换[3]、T矩阵法和驻相法[4]等。其中分离变量法和T矩阵法是通过求解麦克斯韦方程组推导出了球面粒子散射平面波的严谨解释，但在计算中未能描述入射光线经过任意次数的内反射是的散射光强度。

在这种研究背景下，几何光学近似(GOA)方法得到了发展[5]-[8]。在气泡散射特性的研究中几何光学近似方法被认为是一种更好的计算方法，它清晰地解释了气泡的散射机制并且对求大尺寸气泡的计算有更高的计算效率。关于几何光学近似方法的发展最早在1957年Van de Hulst给出了光束入射均匀球形粒子散射强度的几何光学近似方法详细推导过程。在此理论基Hodkinson和Greenleaves在计算中考虑了消光系数，得出了球形粒子光散射的强度分布曲线。对于球形粒子的光散射，后续Xiaobing Zhou等对其不断完善，研究了大尺寸粒子的近场和远场散射并编写了对应的GOA计算程序。随着学者的深入研究，GOA方法得到了一定程度上的发展，在研究中将球形粒子散射问题转化为了非球形粒子是散射问题，从介质分布均匀的粒子扩展到了非均匀分布的粒子。F. Xu等人在GOA方法的理论基础上发展了双层球状粒子前向散射的GOA散射模式，并将其研究成果与Mie理论结果进行比较，获得了高度的一致性，验证了研究成果的准确性。Zhibo Zhang等人利用GOA方法探究了规则多面体粒子的散射特性。相较于平面波与粒子散射的GOA方法，F. Xu等人也研究了轴向高斯波束对均匀球形和均匀椭球形粒子散射的GOA方法[9] [10]。

2. 方法描述

湍流气泡的粒径往往会远大于可见光波长，假设一个椭球形气泡的中心是笛卡尔坐标系的原点，记作O_p。我们考虑一个TEM₀₀高斯光束，设在轴高斯波束束腰半径为ω₀，波长为λ，电矢量沿z轴方向极化，沿x轴方向传播入射到椭球形气泡上，高斯光束束腰中心O_G位于(d，0，0)处(如图1所示)，定义了椭球形气泡的椭球度c = a/b为形状参数，且当c > 1时为长椭球，当c < 1时为扁椭球。a和b分别表示沿x轴方向的半轴和沿y轴方向上的半轴，因此b为椭球形气泡在投影截面(y-z平面)的半径。

Figure 1. Gaussian beam incidence ellipsoidal bubble scattering model

图1. 高斯波束入射椭球形气泡散射模型

高斯光束照射到气泡表面上$A\left(x,y,z\right)$ 处的电场强度复振幅S_G可表示为：

$S_G\left(x,y,z\right)=\left|S_G\right|\exp \left(i{\varphi }_i\right)$ (1)

式中的S_G和Φ_i分别表示电场强度振幅值和高斯波束相位，解析式可表达为：

$\left|S_G\right|=\frac{{\omega }_0}{\omega \left(x\right)}\cdot \exp \left(-\frac{z^2+y^2}{{\omega }^2\left(x\right)}\right)$ (2)

${\varphi }_i\left(x,y,z\right)=-k\left[\left(x-d\right)+\frac{y^2+z^2}{2R}\right]+\arctan \left(\frac{x-d}{x_R}\right)$ (3)

其中ω₀为束腰半径，λ为真空中光波波长。波数k，x平面处的束宽半径$\omega \left(x\right)$ ，瑞利长x_R以及等相面曲率半径R可分别表示为：

$\{\begin{array}{l}k=2\pi /\lambda \\ \omega \left(x\right)={\omega }_0\sqrt{1+{\left(\frac{x}{f}\right)}^2}\\ x_R=\frac{\pi {\omega }_0^2}{\lambda }\\ R=\left(x-d\right)\{1+\left[x_R/\left(x-d\right)\right]{}^2\}\end{array}$ (4)

根据几何光学近似的观点，光束入射到椭球形气泡表面后会与气泡发生相互作用，如图1所示，其中一部分光线会在外表面直接反射(p = 0)，一部分光线会折射进入气泡内部。除了这部分光线之外，还有一部分光线会发生衍射。粒子的散射光就是由反射、折射和衍射光线叠加而成，其复振幅可表达为：

$s_j=s_d+{\displaystyle \sum _{p=0}^{\infty }{s}_{j,p}}j=1,2$ (5)

式中S_j为散射光振幅函数，S_d为衍射光振幅函数，S_j,p为第p阶出射光的振幅函数。j = 1表示垂直极化(振动方向垂直于散射面)分量，j = 2表示平行极化(振动方向平行于散射面)分量。在距离粒子中心距离为r的观测点远场散射强度为：

$I_j=\frac{I_0}{{\left(kr\right)}^2}{i}_i\left(\theta \right)=\frac{I_0}{{\left(kr\right)}^2}{\left|s_j\left(\theta \right)\right|}^{{}_{{}^2}}$ (6)

式中I₀是高斯光束中心的强度。

2.1 散射角

在一阶近似中，高斯波束在任意入射点$A\left(x,y,z\right)$ 处的相位可由公式(3)表示。在几何光学近似中，当高斯波束入射椭球形气泡时，其方向矢量可表示为[11]

$L_o=L_{X}i+L_{y}j+L_{Z}k$ (7)

其中

$L_x=-\frac{\partial \varphi }{\partial x}=k+\frac{k\left(y^2+z^2\right)\left[x^R-\left(x-d\right)\right]}{\left[{\left(x-d\right)}^2+x_R^2\right]}-\frac{x_R}{{\left(x-d\right)}^2+x_R^2}$ (8)

$L_y=-\frac{\partial \varphi }{\partial y}=\frac{2x_R\left(x-d\right)y}{{\omega }_0^2\left[{\left(x-d\right)}^2+x_R^2\right]}$ (9)

$L_z=-\frac{\partial \varphi }{\partial z}=\frac{2x_R\left(x-d\right)z}{{\omega }_0^2\left[{\left(x-d\right)}^2+x_R^2\right]}$ (10)

将入射点坐标代入上述公式并对矢量进行归一化便可得到入射点A处的入射单位向量，再通过斯涅尔定理可以得到折射光线和反射光线，对光线进行矢量追踪，从而得到各阶出射光线的偏转角${{\theta }^{\prime }}_p$ 和散射角${\theta }_p$ 。

${{\theta }^{\prime }}_p=2\pi l+q{\theta }_p$ (11)

2.2. 相位

如图2所示，我们考虑一个被高斯光束照射的椭球形气泡，首先，与球形颗粒类似，我们主要考虑三个方面的相移：由光程差引起的相移Φ_p,pH、由高斯光束波前曲率引起的相移Φ_G和椭球内外焦线引起的相移Φ_P,FL。对于球形颗粒，如果入射角已知，那么球内光线入射角始终保持不变，则球内连续两个内反射点的长度差保持不变，即球内路程差有明确的解析解。而对于椭球形颗粒，当光线进入粒子后，每次光线到达椭球形粒子的边界时入射角就会发生变化，因此椭球形气泡内部连续内反射点之间的长度没有明确的解析表达式，因此我们需要找出内反射之间的路径关系，并逐步进行计算。

如图2所示，我们以O_p为圆心，以半长轴a为半径做参考圆，其中R、R₁为半径与参考圆的交点，则反射光线(p = 0)的相移为

${\varphi }_{0,pH}=k\left(\overline{R{O}_p}-\overline{A{S}_o}+\overline{O_p{R}_0}-\overline{S_0{A}_0}\right)$ (12)

对于直接透射穿过气泡无内反射的光线(p = 1)，其透射光线相移可表示为：

${\varphi }_{1,pH}=k\left(\overline{R{O}_p}-\overline{A{S}_o}\right)-k{m}_r\overline{S_0{S}_1}+k\left(\overline{O_p{R}_1}-\overline{S_1{A}_1}\right)$ (13)

因此，对于经过$p-1$ 次内反射的光线，其相移可表示为：

${\varphi }_{p,pH}=k\left(\overline{R{O}_p}-\overline{A{S}_o}\right)-k{m}_r{L}_p+k\left(\overline{O_p{R}_p}-\overline{S_p{A}_p}\right)$ (14)

其中，L_p表示光线在椭球形气泡内部的总路径。

Figure 2. Phase calculation diagram

图2. 相位计算示意图

综合，p阶出射光线的相移可表示为：

${\varphi }_{p,pH}=\{\begin{array}{l}k\sqrt{x_0^2+y_0^2-\frac{{\left(y_0-x_0\tan \beta \right)}^2}{{\tan }^2\beta +1}}+k\left[\sqrt{x_0^2+y_0^2-\frac{{\left(y_0-x_0\tan {\theta }_0\right)}^2}{{\tan }^2{\theta }_0+1}}\right]p=0\\ k\sqrt{x_0^2+y_0^2-\frac{{\left(y_0-x_0\tan \beta \right)}^2}{{\tan }^2\beta +1}}-k{m}_r{L}_p+k\left[\sqrt{x_p^2+y_p^2-\frac{{\left(y_p-x_p\tan {\theta }_p\right)}^2}{{\tan }^2{\theta }_p+1}}\right]p\ge 1\end{array}$ (15)

其中β为入射光线与x轴的夹角，θ_p为p阶光线的散射角。关于β、θ_p、(x_p, y_p)的确定都是由矢量射线追踪程序确定的。对于椭球截面，当椭球度c接近于1时有：

$L_p=2pa\cos {\theta }_r$ (16)

$k\sqrt{z_0^2+y_0^2-\frac{{\left(y_0-z_0\tan \beta \right)}^2}{{\tan }^2\beta +1}}=ka\cos {\theta }_i$ (17)

$k\sqrt{z_p^2+y_p^2-\frac{{\left(y_p-z_p\tan {\theta }_p\right)}^2}{{\tan }^2{\theta }_p+1}}=ka\cos {\theta }_i$ (18)

通过上式不难发现方程已恢复至球形粒子散射的解析式，其中θ_i为入射角，θ_r为折射角。

对于高斯波束与气泡的散射光线相位，还需要考虑入射波束的波阵面引起的相移Φ_G：

${\varphi }_G={\varphi }_{A{S}_0}-{\varphi }_{PQ}=k\left[a-\sqrt{x_0^2+y_0^2-\frac{{\left(y_0-x_0\tan \beta \right)}^2}{{\tan }^2\beta +1}}\right]-\left[k\left(d+a\right)-{\varphi }_i\right]$ (19)

对于一个球形颗粒，光线出现的焦点线的数目只由相对折射率m决定，而对于椭球体颗粒，则由折射率m和椭球度c共同决定。Van de Hulst提出两个相邻交叉的射线形成一个焦点，相位会推进π/2，因此有：

${\varphi }_{p.FL}=\frac{\pi }{2}\left(1+p-2l-\frac{1}{2}s-\frac{1}{2}q\right)$ (20)

由于椭球形气泡内部空气的折射率小于外部介质的折射率，因此，可能会发生全反射现象，全反射产生的相移为Φ_j,r，其中j = 1或2分别表示垂直或平行分量。

${\varphi }_{1,r}=2{\tan }^{-1}\left[\frac{{\left({\sin }^2{\theta }_{i,p}-1/m_r^2\right)}^{1/2}}{\cos {\theta }_{i,p}}\right]$ (21)

${\varphi }_{2,r}=2{\tan }^{-1}\left[\frac{m_r{}^2{\left({\sin }^2{\theta }_{i,p}-1/m_r^2\right)}^{1/2}}{\cos {\theta }_{i,p}}\right]$ (22)

至此，对于高斯波束入射的情况，散射光线总相位可以表示为

${\varphi }_p=\frac{\pi }{2}+{\varphi }_{p,pH}+{\varphi }_{p.FL}+{\varphi }_G+{\varphi }_{j,r}$ (23)

2.3. 衍射

高斯波束在衍射部分的描述与球形粒子相同，将三维球体转化为半径R = a的二维圆盘，根据Chevaillier等给出的衍射场计算方法[12] [13]，高斯波束对半径为R的圆盘衍射的远场电场强度复振幅公式为：

$s_d=\frac{x^2}{2\pi }\left(\frac{{\omega }_0}{\omega }\right){\displaystyle {\int }_0^{2\pi }d{\phi }^{\prime }}{\displaystyle {\int }_0^1\exp \left(-A{t}^2\right)}\exp \left[i\left(Bt\right)+c{t}^2\right]t\text{d}t$ (24)

其中x为粒子尺寸参量，x = kR，其他参量可表示为$A={\left(a/\omega \right)}^2,B=-x\tan \theta \cos \left(\phi -{\phi }^{\prime }\right),C=\frac{x^2}{2kR}$ 。

2.4. 散射光复振幅

当波束入射椭球形气泡表面时，会有一部分能量被反射，另一部分能量被折射，${\epsilon }_{j,p}^2$ 表示第p阶出射光线所占的强度比例，则光线经过$p-1$ 次内反射后，射线的振幅${\epsilon }_{j,p}$ 可表示为：

${\epsilon }_{j,p}=\{\begin{array}{l}{r}_{j,p}p=0\\ {\left(1-r_{j,0}^2\right)}^{1/2}{\left(1-r_{j,p}^2\right)}^{1/2}\underset{n=1}{\overset{p-1}{\Pi }}\left(-r_{j.n}\right)p\ge 1\end{array}$ (25)

$r_{j,p}$ 为菲涅尔反射系数，j = 1，2分别表示垂直和平行极化

$\begin{array}{l}{r}_{1,p}=\frac{\cos {\theta }_{i,p}-m_r\cos {\theta }_{r,p}}{\cos {\theta }_{i,p}+m_r\cos {\theta }_{r,p}}\\ r_{2,p}=\frac{m_r\cos {\theta }_{i,p}-\cos {\theta }_{r,p}}{m_r\cos {\theta }_{i,p}+\cos {\theta }_{r,p}}\end{array}$ (26)

当粒子入射气泡表面时，会根据颗粒表面的局部曲率收敛或发散，由于气泡局部曲率引起的散度因子D可表示为：

$D_G=\frac{\cos {\theta }_i\sin \tau }{\sin {\theta }_p\left|\frac{d{{\theta }^{\prime }}_p}{d\tau }\right|}$ (27)

与球形颗粒不同，对于椭球形粒子$\text{d}{{\theta }^{\prime }}_p/\text{d}\tau$ 与${{\theta }^{\prime }}_p$ 、τ没有明确的关系，很难得出确定的解析解，则

$\frac{\text{d}{{\theta }^{\prime }}_p^{}}{\text{d}\tau }=\frac{\left|{{\theta }^{\prime }}_{p,t+1}^{}-{{\theta }^{\prime }}_{p,t-1}^{}\right|}{\sqrt{R_{t-1}^2+R_{t+1}^2-2R_{t-1}{R}_{t+1}\cos \left({\tau }_{t+1}-{\tau }_{t-1}\right)}/R_t}$ (28)

综上所述，当高斯波束入射到吸收性气泡时，第p阶光线远场散射强度复振幅函数可表示为[14] [15]：

$S_{j,p}=kR\left|S_G\right|{\epsilon }_{j,p}{D}_G^{1/2}\exp \left(i{\varphi }_p\right)$ (29)

其中S_G为高斯光在入射点的振幅。

3. 数值模拟结果与讨论

本文中，我们考虑高斯光束束腰中心OG与粒子中心Op在同一位置处，在任何情况下，光线经过几次内反射后就会完全消失，因为光线每经过一次反射后光强就会显著降低，所以本文中我们考虑了pmax = 6来计算散射光强。设气泡在x轴上的长度$a=100\text{ μm}$ ，高斯波束波长$\lambda =632.8\text{ nm}$ ，束腰半径${\omega }_0=1000\text{ μm}$ ，气泡内部介质均匀且气泡表面光滑，气泡的光吸收影响忽略不计。

对于椭球形颗粒，当光线进入粒时到达椭球形粒子的边界时入射角会发生变化，所以我们使用矢量射线追踪(VRT)的方法获得散射角与入射角之间的关系，为验证MATLAB程序结果的正确性，我们将椭球形气泡的椭球度$c\to 1$ ，将得到的散射角结果与球形气泡在基于几何光学近似原理下得出的结果进行比较。

由图3可知，当椭球度c = 1.01时，对于p = 0，1，2，3阶光线，使用矢量射线追踪法模拟的结果与几何光学近似方法下的球形气泡的结果具有很好的一致性。验证了矢量射线追踪法模拟计算的散射角的可行性。

为了分析形状参数的影响，我们还给出了散射角和入射角的关系，如图4所示，对于不同的椭球度c值，$p=0$ 阶的曲线几乎完全一致，因此，占有大部分入射光强的$p=0$ 阶的散射光强与球形气泡相似，然而$p=1$ 阶曲线开始出现了部分偏差，对于$p\ge 2$ 阶的曲线出现了更大程度的偏差。对于$c<1$ 的椭球形气泡，在入射角较小的区域散射角图斜率变化较大，在较大的入射角部分，散射角曲线变化较为缓慢。对于$c>1$ 的气泡，散射角曲线变化先是比较缓慢，而后迅速变化。

(a) (b)

(c) (d)

Figure 3. Relation curve of scattering Angle and incidence Angle when c = 1.01. (a) p = 0; (b) p = 1; (c) p = 2; (d) p = 3

图3. c = 1.01时散射角和入射角关系曲线。(a) p = 0；(b) p = 1；(c) p = 2；(d) p = 3

(a) (b)

(c) (d)

Figure 4. Curve of incidence Angle and scattering Angle under different c values. (a) p = 0; (b) p = 1; (c) p = 2; (d) p = 3

图4. 不同c值下入射角与散射角曲线。(a) p = 0；(b) p = 1；(c) p = 2；(d) p = 3

由图5(a)图可以看出气泡的相对折射率对0阶出射光线的散射角几乎没有影响，在不同m值的情况下三条分布图线高度重合，占有大部分入射光强的p = 0阶散射强度几乎不受气泡相对折射率的影响。根据图5(b)~(d)可以看出对于$p\ge 1$ 阶的散射光线其散射角分布受相对折射率影响较大，随着气泡相对折射率的增大，气泡光散射对应的最大散射角在逐渐减小。若GOA散射计算中仅考虑p = 0和p = 1阶光线，根据几何光学近似方法的推导过程，当散射角度大于一阶出射光线的最大散射角时，GOA方法将不能准确计算出其前向散射强度，且从图5中可以看出对于相对折射率小于0的粒子，其相对折射率数值越接近1，p = 1阶散射光线的最大散射角越小，前向散射强度计算角度范围越小。

(a) (b)

(c) (d)

Figure 5. Curve of incidence Angle and scattering Angle under different m values when c=1.01. (a) p = 0; (b) p = 1; (c) p = 2; (d) p = 3

图5. C = 1.01时不同m值下入射角与散射角曲线。(a) p = 0；(b) p = 1；(c) p = 2；(d) p = 3

如图6所示，我们模拟了椭球度c = 1.05，$a=100\text{ μm}$ ，$\lambda =632.8\text{ nm}$ ，m = 0.75的气泡在不同高斯光束束腰半径时椭球形气泡的光强分布。当高斯光束束腰半径无限增大时，照射到气泡上的高斯光束可近似看作平行光。观察图5(b)图可以发现，在散射角接近80˚时，平行极化强度会急剧下降，此时光线的入射角接近全反射临界角。观察图5(a)、图5(b)，我们可以明显的发现随着高斯光束束腰半径的增大，气泡的垂直极化方向和平行极化方向的光强都在增强，且前向散射角区间光强的变化范围也随着高斯光束束腰半径的增大而增大。

(a) (b)

Figure 6. Intensity distribution ((a) Vertical polarization (b) Parallel polarization) of diffused bubble light under different Gaussian beam waist radius

图6. 不同高斯光束束腰半径下气泡散射光强分布。(a) 垂直极化；(b) 平行极化

4. 总结

本文在平面波入射球形气泡的几何光学近似方法的基础上，推导出了轴向高斯光入射椭球形气泡的几何光学近似方法。在文章中我们分析了椭球度对散射角的影响以及高斯光束束腰半径对椭球形气泡散射光强的影响，得出了0˚~180˚散射角之间的平行极化与垂直极化光强分布图。研究发现随着高斯光束束腰半径的增大，气泡的垂直极化方向和平行极化方向的光强逐渐增强，且随着高斯光束的束腰半径对椭球形气泡的前向散射光强的影响较大，而后向散射相对比较稳定，因此可通过其分析气泡对高斯光束散射的光强曲线来反演入射高斯光束的束腰半径。在理论推导和算法实现的过程中可以看出几何光学近似方法可以快速且有效地计算在轴高斯光束照射椭球形大气泡的散射强度分布，研究结果可为优化气泡的测量技术提供理论依据。