简介Mie理论自1908年提出以来一直是分析球形微粒光散射的重要工具基于该理论的散射光强MATLAB代码包面向光学、大气科学、环境监测及生物医学领域的研究人员与学生用于计算和可视化微米级球形颗粒在光照下的散射、消光与吸收特性。包内共10个文件含9个M脚本与1个MAT数据文件压缩包仅34KB小巧易用。代码采用模块化设计覆盖颗粒尺寸参数、折射率设定、核心Mie散射系数、不同角度散射光强分布以及数据可视化等关键环节并附带1.06微米波长数据与测试脚本可直接运行验证或继续二次开发。使用者输入颗粒粒径、折射率和光源波长即可快速获得散射光强、消光系数和散射系数进而用于雾滴、粉尘、气溶胶等粒子散射特性的定量分析与对比研究。目前已有575人学习下载适合具备基础MATLAB编程和光学知识的读者快速入门与深入实践也可用于课堂教学演示、科研验证与工程初步估算。1. Mie散射光强从一次粒径测量说起做颗粒表征的人多半被一张散射光强随角度的曲线难住过。一颗直径 600 nm 的聚苯乙烯球在 532 nm 绿光照射下前向光强很强但在 90° 方向的光强会出现细密的振荡结构。如果图省事套用瑞利散射公式或者用衍射理论去算得到的是一条平滑衰减曲线和实测数据完全对不上。问题出在粒子尺寸和波长处于同一个量级粒子内部激发的多极电磁场相互干涉只有 Mie 理论能把这个干涉过程描述清楚。基于 Mie 理论的散射光强计算简单说就是给定入射波长、粒子半径和复折射率算出任意散射角上的光强与偏振状态同时还能给出消光、散射和吸收三个截面。这篇文章写给手里有实验数据、想用 Mie 理论解释现象或者做仿真预估的工程师和研究生。读完你不仅能写出一个可运行的 Mie 散射光强程序还会知道哪些参数对结果影响大结果不对劲时应该先检查哪个环节。2. Mie理论的物理边界适用条件、截面定义和计算选型2.1 尺寸参数与复折射率决定你能不能硬套瑞利公式Mie 理论是麦克斯韦方程组对均匀球状粒子的解析解它的核心输入只有两个物理量尺寸参数和复折射率。尺寸参数定义为x 2 * pi * r / lambda其中 r 是粒子半径lambda 是粒子周围介质中的波长。注意不是真空波长而是介质中的波长。如果粒子处在水中入射光波长是 532 nm水的折射率约 1.33那么介质中波长约 400 nm。很多人第一次算就把这个细节漏掉结果整条散射曲线向大角度偏移。第二个输入是复折射率。写法最常见的是m n i*k其中实部 n 是粒子相对周围介质的折射率比虚部 k 描述吸收。k 为正数表示光在粒子内部被吸收k 越大吸收越强。金属粒子如金、银虚部可以到 1 甚至更高干净的水滴、PS 微球k 接近 0。什么时候可以不用 Mie 理论经验判断标准是尺寸参数 x。当 x 远小于 1比如 r 只有波长的几十分之一粒子感受到的入射场近似均匀激发的主要是电偶极子模式这时候瑞利散射公式就能给出足够好的结果。当 x 大于 1 但还没到数百高阶多极模式开始起作用瑞利公式失效这就是 Mie 理论的主战场。这个区间覆盖了大多数实际场景大气气溶胶、胶体颗粒、血细胞、工业粉体粒径在 0.05 到 50 微米之间可见光波段下基本都落在 Mie 区。2.2 消光、散射与吸收截面Mie理论输出的三套关键数字除了散射光强的角分布Mie 理论还给出三个积分量消光效率因子 Qext、散射效率因子 Qsca 和吸收效率因子 Qabs。它们满足能量守恒关系Qext Qsca Qabs三个量都是无量纲效率因子物理含义是粒子的实际消光截面除以几何截面 pirr。如果粒子不吸收Qabs 等于 0Qext 等于 Qsca。这三个量的用处各不相同。做粒径反演时透射法测的是消光对应 Qext做散射光测量时对应 Qsca做光热治疗或热效应评估时对应 Qabs。实际计算中Mie 理论先算散射系数 a_n 和 b_n再用无穷级数求和得到 Qsca 和 Qext最后相减得到 Qabs。后文代码里我会把这三个量一并输出方便直接用于实验数据标定。值得注意Qsca 随粒子半径的变化不是单调的在 x 从 0.1 增长到 10 的过程中Qsca 会先快速上升、后出现振荡。这意味着在粒径反演时不能简单用“散射光越强、粒子越大”来判断必须结合振荡结构甚至多角度信息才能唯一确定粒径。2.3 什么时候必须用Mie什么时候可以偷懒三种近似边界速查工程上最常犯的错误是近似的过度使用。我一般用一个快速判断表来决定是否还要跑 Mie近似方法适用条件误差特征典型场景瑞利散射x 0.1且mx 不特别大Mie 理论0.1 x 1000数值收敛误差可控微米颗粒、细胞、气溶胶夫琅禾费衍射x 1000且折射率接近 1忽略粒子内部吸收和多极细节大颗粒粒度仪、挡光法注意最后一行x 大于 1000 时粒子等效于一个带相位盘的圆孔衍射图样主要取决于几何投影面积不再包含折射率信息。在激光粒度仪里常用这一近似但如果颗粒有强吸收衍射近似会低估消光截面这时候还是回到 Mie 更稳妥。一个辅助判断办法是先算一次 Mie再算一次瑞利或夫琅禾费近似比较两个结果。如果相对偏差小于你需求的误差容限就把近似方案确定下来后续批量扫描可以省下大量计算时间。否则老老实实跑 Mie不要在这个问题上凭感觉选。3. 用Python算散射光强一套可以直接跑的Mie实现3.1 用scipy计算a_n和b_n散射系数核心代码与数值稳定约定Mie 系数的计算在教科书里有标准公式但直接照抄会踩不少数值坑。下面这份代码基于 scipy.special 的球贝塞尔函数实现可读性好在常见尺寸参数范围内数值稳定。先定义系数计算函数import numpy as np from scipy.special import spherical_jn, spherical_yn def mie_coefficients(m, x): 计算 Mie 散射系数 an 和 bn。 参数 ---- m : complex 粒子折射率相对周围介质的比值如 1.33 0.0j x : float 尺寸参数x 2*pi*r/lambdalambda 为介质中波长 返回 ----- an, bn : ndarray 散射系数数组下标从 n1 开始 # 经验截断阶数保证级数收敛 nmax int(np.ceil(x 4.0 * x**(1.0/3.0) 2.0)) nmax max(nmax, 5) n np.arange(1, nmax 1) z m * x # 球贝塞尔函数及其导数 # spherical_jn(n, x) 得到 j_n(x)n 可以是数组 jn_x spherical_jn(n, x) jn_z spherical_jn(n, z) yn_x spherical_yn(n, x) # 下标 n-1 对应的球贝塞尔函数 jn_x_prev spherical_jn(n - 1, x) jn_z_prev spherical_jn(n - 1, z) yn_x_prev spherical_yn(n - 1, x) # psi_n(z) z * j_n(z) psi_x x * jn_x psi_z z * jn_z # xi_n(x) x * (j_n(x) i * y_n(x)) xi_x x * (jn_x 1j * yn_x) # 导数 d psi_n / d x x * j_{n-1}(x) - n * j_n(x) psi_x_deriv x * jn_x_prev - n * jn_x psi_z_deriv z * jn_z_prev - n * jn_z # d xi_n / d x 同理 xi_x_deriv x * (jn_x_prev 1j * yn_x_prev) - n * (jn_x 1j * yn_x) # Mie 系数标准公式 an (psi_z_deriv * psi_x - m * psi_z * psi_x_deriv) / ( psi_z_deriv * xi_x - m * psi_z * xi_x_deriv) bn (m * psi_z_deriv * psi_x - psi_z * psi_x_deriv) / ( m * psi_z_deriv * xi_x - psi_z * xi_x_deriv) return an, bn这段代码的关键在于用 scipy 的球贝塞尔函数直接构造 psi、xi 和它们的导数。psi_n(z) 的标准定义是 z 乘以球贝塞尔函数 j_n(z)它的解析导数可以用递推关系写成 xj_{n-1}(x) - nj_n(x)这样避免了数值微分带来的误差。参数说明n 从 1 开始取因为 n0 项对应单极子在普通非磁性球散射中不贡献。截断阶数 nmax 采用 x 4x^(1/3) 2 的经验公式这是一个广泛使用的估计值但并不是绝对收敛判据。如果你的计算在相邻 nmax 下结果差异明显需要调大 nmax后文避坑章节会展开讲。3.2 从S1/S2到角分布光强角度离散和递推初始化散射系数的下一步是计算角函数 pi_n 和 tau_n然后累加出两个复振幅函数 S1(theta) 和 S2(theta)。S1 对应垂直偏振分量的散射振幅S2 对应平行偏振分量。对非偏振入射光散射光强正比于两者的模平方平均值。角函数用递推计算注意 pi_0 在数学上必须初始化为 0pi_1 初始化为 1这个细节很多人写错def mie_s12(an, bn, mu): 由 Mie 系数计算 S1 和 S2。 参数 ---- an, bn : ndarray 由 mie_coefficients 得到 mu : ndarray cos(theta)theta 为散射角 返回 ----- S1, S2 : ndarray 复振幅函数 nmax len(an) n np.arange(1, nmax 1) # pi 和 tau 的递推theta 可以是一次性传入的数组 pi_arr np.zeros((nmax, mu.size)) tau_arr np.zeros((nmax, mu.size)) # 初始化 pi_1 1, tau_1 mu并保存 pi_0 0 pi_prev2 np.zeros_like(mu) # 相当于 pi_0 pi_prev1 np.ones_like(mu) # 相当于 pi_1 for idx, nn in enumerate(n): if nn 1: pi_arr[idx] pi_prev1 tau_arr[idx] mu * pi_prev1 else: pi_arr[idx] ((2*nn - 1) / (nn - 1)) * mu * pi_prev1 - ( nn / (nn - 1)) * pi_prev2 tau_arr[idx] nn * mu * pi_arr[idx] - (nn 1) * pi_prev1 pi_prev2 pi_prev1 pi_prev1 pi_arr[idx] factor (2*n 1) / (n * (n 1)) # 级数累加an、bn 是复数数组pi、tau 是实数数组 S1 np.sum(factor[:, None] * (an[:, None] * pi_arr bn[:, None] * tau_arr), axis0) S2 np.sum(factor[:, None] * (an[:, None] * tau_arr bn[:, None] * pi_arr), axis0) return S1, S2初始化时 pi_prev2 对应 pi_0必须设为全 0 数组因为在 n2 的递推中会用上它。tau_1 mu 是从定义式推导出来的不要手写成 1。得到 S1、S2 后非偏振入射光的散射光强角分布为def mie_intensity(S1, S2, wavelength, r, dist1.0): 把 S1/S2 转换成散射光强 I_theta。 返回的单位与入射光强 I0 成比例dist 是观测距离。 x 2.0 * np.pi * r / wavelength prefactor wavelength**2 / (4.0 * np.pi**2 * dist**2) intensity prefactor * (np.abs(S1)**2 np.abs(S2)**2) / 2.0 return intensity这里 wavelength 必须与计算 x 时使用的介质波长一致。如果入射光强记为 I0散射光强就是 I0 * intensity。这个公式的推导基于辐射度学距离 dist 取 1 m 时得到的是单位距离上的归一化值实际实验中有明确的探测距离填对应数值即可。3.3 单点调用与参数化扫描快速定位粒径或折射率的影响把上面的函数组合起来一次 Mie 散射光强计算的完整流程如下import numpy as np import matplotlib.pyplot as plt # 粒子参数半径 300 nm相对折射率 1.330i r 300e-9 # 单位米 n_medium 1.0 # 空气 n_particle 1.33 m n_particle / n_medium wavelength_vac 532e-9 wavelength wavelength_vac / n_medium # 介质中波长, 空气中可忽略差异 x 2 * np.pi * r / wavelength an, bn mie_coefficients(m, x) theta_deg np.linspace(0.1, 180, 901) # 0.1 度开始避免数值奇点 mu np.cos(np.deg2rad(theta_deg)) S1, S2 mie_s12(an, bn, mu) I_theta mie_intensity(S1, S2, wavelength, r, dist1.0) # 归一化画图方便看振荡结构 I_norm I_theta / np.max(I_theta) plt.figure(figsize(6, 4)) plt.semilogy(theta_deg, I_norm) plt.xlabel(Scattering angle (deg)) plt.ylabel(Normalized intensity) plt.show()这段代码跑出来的曲线在 0° 方向最强随着角度增大出现多个峰谷后向 180° 附近还有一个不可忽略的局部峰。Mie 散射的这种振荡是干涉的直接结果用实验数据拟合时这些峰谷的位置就是粒径反演的核心特征。参数化扫描也很简单把 r 换成一个数组在循环里重复以上步骤r_list np.linspace(50e-9, 1000e-9, 50) intensity_90 [] for rr in r_list: x_tmp 2 * np.pi * rr / wavelength an_tmp, bn_tmp mie_coefficients(m, x_tmp) S1_tmp, S2_tmp mie_s12(an_tmp, bn_tmp, np.cos(np.deg2rad(90.0))) I_tmp mie_intensity(S1_tmp, S2_tmp, wavelength, rr, dist1.0) intensity_90.append(I_tmp)扫描后你会发现 90° 散射光强随粒径的变化呈现明显振荡而不是单调增长。这是 Mie 区粒子的典型行为所以在实验里只取单一角度的光强做粒径反演很容易得到多个解必须配合多角度或消光谱数据。4. Mie散射光强计算避坑指南五个高频翻车点4.1 坑一复折射率虚部符号写反吸收计算结果直接崩掉现象算出来的 Qabs 是负值或者散射光强在某一粒径区间出现不可能的高斯尖峰。原因Mie 公式推导时隐含了时间因子约定。国内教材常用 exp(-iomegat)这种情况下复折射率虚部取正表示吸收但有些国外代码或数值包采用 exp(iomegat) 约定虚部符号就要反过来。把符号写反本质上相当于把吸收介质当成了增益介质能量不守恒。解决统一把自己代码中的复折射率写成 m n i*k并且 k 0。然后在算完 Qabs 后立刻检查Qabs 必须大于等于 0且 Qext Qsca Qabs 严格成立。如果程序输出 Qabs 为负先改折射率虚部的符号不要怀疑是其他问题。4.2 坑二截断阶数n_max取太小散射效率曲线出现伪振荡现象粒径扫描时 Qsca 曲线出现锯齿状抖动相邻两个半径的计算结果相差很大而且抖动幅度随时间因子或网格变化不稳定。原因Mie 级数是无穷级数实际计算必须截断。锥形截断阶数取太小高阶项的贡献被截掉但更高阶的散射系数并不单调递减反而会在某个 n 附近出现共振所以截断位置不当会引入伪振荡。解决默认按 n_max ceil(x 4 * x^(1/3) 2) 取但不算完。计算完成后把 n_max 加 10 重算一次对比 Qsca 和几个关键角度的光强。如果相对变化超过 0.1%继续增大 n_max 直到结果稳定。金属粒子吸收强时需要的阶数可能比这个经验公式更大这一点我在算金纳米球时吃过亏。4.3 坑三角度网格太粗后向散射细节被抹平现象实验测得的角分布曲线在 120° 到 170° 之间有明显起伏但程序画出来是一条平滑下降的线。原因Mie 角分布中包含的振荡分量阶数很高最细的振荡结构与 n_max 直接相关。角度离散步长如果大于半个振荡周期振荡特征就被混叠滤掉了。常见做法是用 1° 步长但对 2 微米以上的粒子这个步长已经太粗。解决角度步长取始终不大于 0.5°即 0.0087 rad必要时在 90° 到 180° 区间加密网格。也可以用自适应网格先粗算定位峰谷再在峰谷附近细化。角网格越细级数累加的计算量越大但对现代计算机来说901 个角度点和 200 个散射系数阶数的组合耗时不过几十毫秒不值得在精度上省。4.4 坑四把真空波长当介质波长粒子越小误差越大现象实验结果和仿真对不上但把折射率改成不同的水源值偏差仍然存在而且偏差方向固定。原因尺寸参数 x 2pir/lambda 中的 lambda 是粒子周围介质中的波长。如果粒子在水里真空波长 532 nm 必须除以水的折射率约 1.33得到约 400 nm。直接用 532 nm 相当于把粒子尺寸参数缩小了 1.33 倍散射曲线整体偏移尤其在粒径接近波长时偏移非常明显。解决把所有输入统一成介质内波长和相对折射率。定义 wavelength_med wavelength_vac / n_mediumm n_particle / n_medium。这一步看似简单但在多层介质或微流控芯片场景中水、玻璃、空气的折射率混在一起最容易忘记某个界面上的换算。4.5 坑五散射角用度还是弧度递推公式里搞混现象程序报错、输出 NaN或者 0° 和 180° 方向的结果明显不对称。原因Mie 角函数递推基于 mu cos(theta)如果 theta 是角度制数值cos(90) 不等于 0递推的初始条件全部错位级数结果自然不对。另一个变种是把 mu 直接当成角度传给绘图函数导致坐标轴错乱。解决在代码入口统一规定theta 以角度制接收内部立即转弧度并计算 mu。不要在函数里来回转换。写完函数后做一个最小自检theta0 时 S1 和 S2 应该相等theta180 时散射光强应该在物理合理范围。这两条可以通过直接打印比较确认。5. 散射光强的工程应用粒径反演、偏振特征与近似退化边界5.1 从光强振荡反演粒径极值位置的初步估算与最小二乘拟合Mie 散射光强角分布中含有粒径信息最直观的利用方式是提取振荡峰谷的位置。随着粒径增加前向衍射峰变窄侧向和后向的峰谷数量增多。一个粗略的经验是在 0° 到 90° 范围内峰谷总数近似正比于尺寸参数 x。先用这个关系做粗估再用最小二乘拟合精确定位。拟合的目标函数是实测角分布与 Mie 计算值的残差平方和from scipy.optimize import least_squares def residual(params, measured_intensity, theta_deg, wavelength): r_guess params[0] m_guess 1.33 0j # 折射率可固定也可作为拟合参数 x_tmp 2 * np.pi * r_guess / wavelength an, bn mie_coefficients(m_guess, x_tmp) mu np.cos(np.deg2rad(theta_deg)) S1, S2 mie_s12(an, bn, mu) I_calc mie_intensity(S1, S2, wavelength, r_guess, dist1.0) # 归一化后比较避免绝对强度标定误差 I_calc I_calc / np.max(I_calc) measured_norm measured_intensity / np.max(measured_intensity) return I_calc - measured_norm result least_squares(residual, x0[300e-9], args(measured, theta_deg, wavelength))注意这里只拟合粒径一个参数条件是折射率已知。折射率未知时建议先用多角度数据联立拟合 r 和 m但这样容易陷入局部最优。我的做法是先固定折射率粗拟合 r再固定 r 拟合折射率迭代两三轮后就稳定了。拟合时必须使用归一化光强因为实验中的绝对光强受探测效率、激光功率等因子影响很难精确标定。归一化后拟合问题对强度绝对尺度不敏感但要注意实测数据不要包含多次散射的背景否则角分布形态会被污染。5.2 偏振比与S1/S2夹角从单一角度光强中挤出粒子形状信息在 90° 散射方向垂直偏振分量的散射光强与平行偏振分量的比值对粒子尺寸非常敏感。对球形粒子这个偏振比可以直接由 S1 和 S2 算出polarization_ratio np.abs(S1)**2 / np.abs(S2)**2在瑞利极限下S2 在 90° 趋向于 0偏振比趋于无穷大。随着粒径增大S2 不再为零偏振比显著下降。因此测量单一角度的偏振比就能对粒径范围做初步判断这也是很多流式细胞仪“侧向散射 偏振通道”设计的基本依据。如果粒子不是标准球S1 和 S2 的相位差不再遵循 Mie 理论预期偏振比会偏离球形假设。所以偏振比曲线也可以用来检验粒子球形度实测偏振比与球形 Mie 计算值偏差过大说明粒子形貌或内部结构不均匀。这时候不要强行用 Mie 拟合更合适的是 T 矩阵或离散偶极近似。5.3 瑞利散射与夫琅禾费衍射的替换边界一张对照表工程上避免重复跑 Mie 的办法是先判断能否用近似。我在实际项目中总结了一张更保守的替换边界表按精度需求分成两档精度需求瑞利散射边界Mie 必须使用夫琅禾费衍射边界误差 5%x 0.150.15 ≤ x ≤ 800x 800 且 n−1 较小误差 1%x 0.080.08 ≤ x ≤ 1200x 1200 且 n−1 0.05注意折射率接近 1 是夫琅禾费衍射成立的前提。如果粒子是强吸收材料即使 x 很大粒子内部的相位变化仍会影响前向散射衍射近似可能低估消光截面。遇到这种情况最稳妥的做法是在批量计算前先做几次探针计算用 Mie 结果和近似结果对比确认误差可接受后再启用近似。替换边界的另一个工程含义是市售激光粒度仪通常宣称测量范围从 0.01 微米到数千微米低频端用瑞利近似、高频端用夫琅禾费近似中段用 Mie 是完全正常的。但如果你要分析的是强吸收或高折射率颗粒必须确认仪器软件在对应粒径段使用了 Mie 模型否则反演结果会出现系统偏差。6. 一束光强算完怎么验证自洽性检验与我的实操习惯6.1 三种快速自检能量守恒、截断收敛和积分一致计算完成不等于结果正确。我每次算完一组 Mie 散射光强都会先做三个不依赖外部数据的自检。能量守恒是第一条Qabs 必须不小于 0Qext 等于 Qsca 加 Qabs。偏差超过 1% 时检查 Mie 系数公式中的分子分母是否写反或折射率符号是否正确。截断收敛是第二条把 n_max 提高 10 到 20 阶重算 Qsca 和 90° 光强相对变化应小于 0.1%。若变化明显说明截断不够回到避坑章节第 4.2 节的处理方法。数值积分是第三条把算出的角分布光强在 4π 立体角上积分再与散射截面 pirr*Qsca 对比。由于 S1、S2 和强度公式都归一化到入射光强 I0两者应该一致到 1% 以内。这个检验能发现角度网格太粗或漏掉奇点的问题比肉眼观察曲线可靠得多。6.2 一次可信计算的完整动作清单我习惯把一次完整的 Mie 散射光强计算固化成固定流程防止反复改参数时漏掉步骤步骤动作检查点1确认介质折射率与真空波长换算介质内波长wavelength_med 是否正确2确认粒子折射率写成 nik 且 k0Qabs 不应为负3计算 x选定 nmax 并加 10 复算验证Qsca 相对变化 0.1%4离散角度步长不大于 0.5°90° 到 180° 是否保留细节5计算 S1/S2 和光强归一化后与实验对比0° 与 180° 无异常6立体角积分与 Qsca 对比偏差 1%我自己的习惯是把这些检查项写进一个 assert 函数里每次计算完自动跑一遍。曾经有一次拟合结果一直偏大排查了半小时最后发现是输入粒子半径用了直径数值自检函数里的能量守恒项立刻暴露了异常。这个教训让我明白Mie 计算本身不是黑匣子但人对错误的容忍度却是有限的。希望这套流程能帮你少走这些弯路在散射光强计算上更早拿到可信的结果。本文还有配套的精品资源点击获取
