基于Mie理论的散射光强计算：Python实现与避坑经验

简介这份MATLAB代码包围绕Mie理论实现散射光强计算与分析适用于大气物理、环境光学、生物医学检测等方向的研究者与学生帮助解决微小颗粒散射建模与可视化问题。包内共10个文件以.m脚本为主另含一个.mat数据文件压缩包仅34KB结构紧凑。脚本覆盖了颗粒尺寸参数计算、消光系数求解、折射率定义、散射振幅与特定角度光强计算等核心环节并提供可视化绘图与测试脚本便于直接运行和验证。数据文件则提供了特定波长如1.06微米的光谱信息用户只需输入粒径、折射率和波长即可快速获得散射光强分布图并分析消光特性适用于污染物监测、光学材料评估等场景。已有575人浏览/学习适合具备MATLAB基础、希望绕开复杂级数推导并直接应用Mie理论建模的读者也可作为教学演示与科研参考脚本使用。1. 基于Mie理论的散射光强为什么球形粒子的散射角分布让人反复翻车做颗粒物光学测量的人早晚会撞上这个题目一束激光打到直径和波长同量级的球上散射光强在每个方向都不一样甚至会在前向出现尖峰、后向出现零碎振荡。用瑞利散射公式估一下可以但一旦尺度参数超过 1角度分布就不是简单的余弦形状误差可以到几个数量级。基于Mie理论的散射光强计算就是通过严格求解电磁场在球形边界上的边值问题把任意粒径、折射率、波长组合下的散射系数和角分布算出来。它要回答的是“已知颗粒和光某个方向散射光强是多少”反过来也能从多点光强反演粒径分布。适合粒度仪、气溶胶消光、乳浊液浓度检测、生物组织模型这些方向。接下来给出能直接跑的 Python 实现并把最容易踩的参数符号、远场条件和级数截断拉出来说清楚。2. 从 Maxwell 方程到散射光强先算 an/bn再把 S1/S2 换成瓦每球面度2.1 Mie 理论到底在算什么严格解和瑞利近似的边界Mie 解的核心思想并不复杂把入射平面波、颗粒内部场和外部散射场都按球谐函数展开然后在球面上匹配电场和磁场的切向分量。匹配结果把所有信息压缩成两组散射系数 an 和 bn其中 n 从 1 取到某个足够大的阶数。an 对应电多极子贡献bn 对应磁多极子贡献。它们只和两个无量纲量有关尺度参数 x π·d / λ_medium以及颗粒相对周围介质的复折射率 m。瑞利近似是 Mie 理论在 x 远小于 1 时的特例这时候只需保留 n1 的项散射光强随角度呈 (1cos²θ) 分布强度还强烈依赖 λ^-4。可是当直径涨到和波长同一量级高阶多极子开始起作用前向散射明显变强后向会出现干涉振荡。这种“大颗粒不是变大版瑞利”的现象是许多做尘埃散射或气泡检测的人第一次翻车的根源。Mie 理论能统一覆盖从分子散射到毫米波雨滴的尺度范围这也是它在雷达气象和激光粒度仪里长期没被替代的原因。2.2 一份能用的 Mie 散射系数代码截断阶数和贝塞尔函数的选择常见做法是直接调 scipy 的球贝塞尔函数把 an/bn 用 Riccati-Bessel 函数写出。这个版本不是最快的但方便逐行检查适合 x 在几百以内的情况。import numpy as np from scipy.special import spherical_jn, spherical_yn def mie_an_bn(m, x): 返回散射系数 an, bn从 n1 开始。 m : 相对复折射率例如 (n_particle 1j*kappa_particle)/n_medium x : 尺度参数x pi * 直径 / 介质内波长 # Wiscombe 经验截断后面再补一项安全余量 n_stop int(x 4.0 * x**(1.0 / 3.0) 2.0) 1 idx np.arange(n_stop 1) jn_x spherical_jn(idx, x) jnp_x spherical_jn(idx, x, derivativeTrue) yn_x spherical_yn(idx, x) ynp_x spherical_yn(idx, x, derivativeTrue) jn_mx spherical_jn(idx, m * x) jnp_mx spherical_jn(idx, m * x, derivativeTrue) # Riccati-Bessel 函数 psi z*j_nxi z*(j_n i*y_n) # 导数用 d(z*j_n)/dz j_n z*j_n 直接算 psi_x x * jn_x psi_x_d jn_x x * jnp_x xi_x x * (jn_x 1j * yn_x) xi_x_d jn_x x * jnp_x 1j * (yn_x x * ynp_x) psi_mx (m * x) * jn_mx psi_mx_d jn_mx (m * x) * jnp_mx an np.zeros(n_stop, dtypecomplex) bn np.zeros(n_stop, dtypecomplex) for n in range(1, n_stop 1): a_num m * psi_mx[n] * psi_x_d[n] - psi_mx_d[n] * psi_x[n] a_den m * psi_mx[n] * xi_x_d[n] - psi_mx_d[n] * xi_x[n] b_num psi_mx[n] * psi_x_d[n] - m * psi_mx_d[n] * psi_x[n] b_den psi_mx[n] * xi_x_d[n] - m * psi_mx_d[n] * xi_x[n] an[n - 1] a_num / a_den bn[n - 1] b_num / b_den return an, bn这段代码里最值得注意的调参点是n_stop。截断阶数取x 4*x**(1/3) 2是 Bohren-Huffman 书上的经验式对中等折射率颗粒通常能把散射截面的误差压到 1e-8 以下。这里的m必须是复数即使材料无吸收也要写1.33 0j否则后面算 S1/S2 时有些分支会出实数错误。spherical_jn(idx, x, derivativeTrue)一次性返回所有阶的导数省去手写递推代价是 x 几千后精度变差这个留给避坑章说。2.3 散射幅度 S1/S2 到光强的换算远场项别漏 r²有了 an/bn 还不够散射角分布需要先组合出幅度函数 S1(θ) 和 S2(θ)。这里要用到角函数的递推def mie_s1_s2(an, bn, mu): 由 an/bn 和 cos(theta) 计算散射幅度 S1, S2。 n_max len(an) pi_n np.zeros(n_max 1) tau_n np.zeros(n_max 1) pi_n[1] 1.0 tau_n[1] mu for n in range(2, n_max 1): pi_n[n] ((2 * n - 1) * mu * pi_n[n - 1] - n * pi_n[n - 2]) / (n - 1) tau_n[n] n * mu * pi_n[n] - (n 1) * pi_n[n - 1] s1 0j s2 0j for n in range(1, n_max 1): factor (2 * n 1) / (n * (n 1)) a an[n - 1] b bn[n - 1] s1 factor * (a * pi_n[n] b * tau_n[n]) s2 factor * (a * tau_n[n] b * pi_n[n]) return s1, s2然后把 S1/S2 转成远场光强。对非偏振入射光角分布要取两个偏振态的平均def mie_intensity_farfield(m, wavelength_medium, diameter, theta_deg, I01.0, r1.0): 计算单个球形颗粒的远场散射光强。 wavelength_medium : 介质内波长单位与 diameter 一致 theta_deg : 散射角可以是数组 r : 观测点到颗粒中心的距离单位与 wavelength_medium 一致 theta_deg np.atleast_1d(theta_deg) x np.pi * diameter / wavelength_medium an, bn mie_an_bn(m, x) mu np.cos(np.deg2rad(theta_deg)) s1 np.zeros(len(theta_deg), dtypecomplex) s2 np.zeros(len(theta_deg), dtypecomplex) for i, mu_i in enumerate(mu): s1[i], s2[i] mie_s1_s2(an, bn, mu_i) k 2.0 * np.pi / wavelength_medium # 非偏振|S1|^2 |S2|^2 后除以 2 intensity I0 / (k * r) ** 2 * (np.abs(s1) ** 2 np.abs(s2) ** 2) / 2.0 return intensity, x, s1, s2这里最容易被忽略的是r。Mie 理论给的是远场渐进解光强必须带 (1/(k·r))² 的衰减。如果你把 r 取成毫米级、颗粒直径十几微米这条公式没问题可一旦 r 小到和颗粒直径可比计算出来的“散射光强”就不再是远场而是近场干涉拿去和实验对比会系统性偏离。另外如果入射光是线偏振且偏振方向相对散射面夹角为 φ强度公式要换成 |S1|²sin²φ |S2|²cos²φ不能继续用非偏振平均式。3. 用 Python 算单粒子散射光强角分布复折射率、粒径和波长要统一写对3.1 复折射率的符号约定和单位换算复折射率通常写成 m n i·κ其中 κ 是吸收指数。这里的第一大坑是时间因子约定。按标准光学约定 exp(-iωt)κ0 表示吸收但不少老代码用 exp(iωt)对应 κ0。你从别的源码抄 an/bn 公式时如果发现吸收材料算出来的消光截面为负十有八九是符号没对齐。我一般会在所有接口注释里固定写死“本模块一律用 exp(-iωt)m 虚部为正表示吸收。”单位方面波长和直径必须同单位。如果你在真空或空气中工作直接用波长 λ0 当 wavelength_medium 就行如果在水中测试颗粒折射率是相对于水的介质内波长要取 λ0 / n_water。一个常见错误是拿真空波长算 x又把颗粒折射率按真空值填导致尺度参数偏大一圈角分布整体错位。3.2 计算角分布光强的脚本和三个调参位置把上一章的函数串起来跑一个直径 2 微米水珠在 532 纳米光下的角分布# 参数颗粒相对空气折射率水在 532nm 附近取实部 1.33 m 1.33 1e-8j wavelength_medium 0.532 # 微米空气/真空 diameter 2.0 # 微米 theta np.linspace(0, 180, 181) I, x, s1, s2 mie_intensity_farfield( m, wavelength_medium, diameter, theta, I01.0, r1000.0 ) # 打印尺度参数和 90° 方向光强供自检 print(x , x) print(I(90°) , I[90])这段脚本的三个调参位置分别是m、wavelength_medium和r。m的实部决定颗粒和介质的光学对比度虚部控制吸收wavelength_medium影响所有尺度相关的振荡周期r决定了光强的绝对量级如果只想看归一化角分布把r固定成 1000 微米即可因为归一化时它会约掉。跑完后你把角度数据画出来会看到前向不再服从瑞利余弦。直径 2 微米、波长 532 纳米对应 x 约 11.8前向峰已经很明显后向出现若干极小值点。这些振荡不是数值噪声是干涉项在角方向的真实表现。3.3 用光学定理检验你的散射效率拿到 an/bn 后第一步自检永远是用光学定理验证消光截面。Mie 理论里消光效率 Q_ext 很容易直接从系数累加def mie_extinction_and_scattering(an, bn, x): 由 an/bn 返回消光截面效率 Q_ext 和散射截面效率 Q_sca。 n_range np.arange(1, len(an) 1) coeff_sum (2 * n_range 1) * (an bn).sum() q_ext coeff_sum.real * 2.0 / x ** 2 q_sca 2.0 / x ** 2 * ((2 * n_range 1) * (np.abs(an) ** 2 np.abs(bn) ** 2)).sum() return q_ext, q_sca对于非吸收颗粒Q_ext 应等于 Q_sca因为能量没有被吃掉对于吸收颗粒Q_abs Q_ext - Q_sca必须为正。光学定理还告诉你消光截面正比于前向散射幅度 S(0) 的实部C_ext (4π/k²)·Re(S1(0))。你在 0° 方向算出的 S1 如果和系数累加对不上说明 π/tau 递推或 an/bn 公式里至少有一处抄错了。这个自检能筛掉九成实现错误。4. 从单颗粒到颗粒群粒径分布加权与平均散射光强4.1 为什么工程里总要对粒径分布做加权真实颗粒体系极少是单分散的激光粒度仪里看到的角分布是成千上万颗粒的叠加。把单颗粒散射光强当成“核函数”把粒径分布函数和它卷积得到的就是群散射角分布。反过来从测到的光强分布反推粒径分布就是粒度仪的反演问题。这里有个容易犯糊涂的点用数量分布加权还是体积分布加权。如果你关心的是颗粒数目浓度比如空气里的细颗粒物计数用数量分布如果你关心的是质量浓度比如乳液浊度通常要乘 d³因为大颗粒虽少但贡献了大部分散射面积和质量。两种加权下的平均散射角分布会差很多尤其在宽分布体系里。4.2 用对数正态分布加权实现平均角散射强度颗粒粒径分布常用对数正态分布描述参数是几何中位径 d_med 和几何标准差 sigma_g。下面这段代码算一组粒径下的加权平均角分布def lognormal_pdf(d, d_med, sigma_g): if d 0: return 0.0 return 1.0 / (np.sqrt(2.0 * np.pi) * d * sigma_g) * \ np.exp(-(np.log(d / d_med) ** 2) / (2.0 * sigma_g ** 2)) def ensemble_intensity(m, wavelength_medium, d_med, sigma_g, theta_deg, weigh_bynumber): 计算对数正态分布颗粒群的平均远场散射光强。 weigh_by: number 数量加权volume 体积加权乘 d^3 theta_deg np.atleast_1d(theta_deg) # 采样范围取中位径的 exp(±3*sigma_g)覆盖 99.7% 以上质量 d_grid np.linspace(d_med * np.exp(-3.0 * sigma_g), d_med * np.exp(3.0 * sigma_g), 128) pdf np.array([lognormal_pdf(d, d_med, sigma_g) for d in d_grid]) if weigh_by volume: pdf pdf * d_grid ** 3 # 每一行是一个粒径下的角分布 I_matrix np.empty((len(d_grid), len(theta_deg))) for i, d in enumerate(d_grid): I_matrix[i, :], _, _, _ mie_intensity_farfield( m, wavelength_medium, d, theta_deg, I01.0, r1000.0 ) # 对粒径做数值积分并归一化到概率质量 numerator np.trapz(I_matrix * pdf[:, None], d_grid, axis0) denominator np.trapz(pdf, d_grid) return numerator / denominator这段代码把 128 个粒径分别算一遍 Mie 散射再按概率密度加权积分。sigma_g是几何标准差无量纲典型值在 1.1 到 2.0 之间取 1 意味着单分散代码会因为 d 范围收缩到一点而退化实际使用时至少给 1.05。d_grid的边界按 exp(±3·sigma_g) 截断避免积分上限拖到几乎贡献为 0 的粒径。如果粒径范围很宽或者需要对更多角度做实时计算128 次完整 Mie 循环会慢。常见做法是先对粒径网格算一次散射截面和角分布再在时间维度缓存。工程上没人想对 1 万个粒径逐个跑递归换成一个稀疏网格加插值精度下降不超过几个百分点速度能快几十倍。4.3 从角分布光强反推粒径的几个坑反演不是直接把群散射光强对角求导。因为不同粒径产生的角分布高度相似反演是一个病态问题需要用 Tikhonov 正则化或截断奇异值分解。我在实际项目里的习惯是先把前向模型写得尽量准确再叠加零均值高斯噪声跑完反演看恢复的 d_med 和 sigma_g 是否在误差范围内。如果前向模型里折射率虚部填错反演出来的粒径分布会整体偏移。更隐蔽的问题是探测器只覆盖有限角度范围比如只能测 10° 到 170°此时前向大颗粒的部分信息丢失反演结果对 20 微米以上的大粒子几乎不敏感。这时候需要联合消光系数或浊度一起反演单靠光强角分布不够。5. 基于 Mie 理论的散射光强计算避坑五条血泪经验5.1 现象粒径一大散射系数开始出现 NaN现象同样的代码x 小于 50 时算得顺滑x 超过 200 时 an/bn 开始吐 NaN或者级数中间出现 1e30 量级的超大值。原因Riccati-Bessel 函数 ψ_n 和 ξ_n 在阶数超过 x 之后一个指数增长一个指数衰减直接由球贝塞尔函数乘积再相减有效数字会全部损失。这是浮点精度问题不是物理问题。解决把 an/bn 换成基于对数导数 D_n(z) 的递推D_n 用从高阶向低阶的向后递推。这也是 Bohren-Huffman 书里经典算法的核心。平时我保留两条路径x 小于 200 用上面的 scipy 版本x 大于等于 200 走 D_n 递推并把两条路径在 x200 附近对比 Q_ext误差小于 1e-8 才认为并线正确。5.2 现象折射率虚部符号不对吸收变“增益”现象算碳颗粒或金属颗粒时消光截面出现负值或者吸收截面 Q_abs 为负。原因复折射率虚部符号完全取决于时间因子约定。按本模块的 exp(-iωt)m n iκκ0 是吸收但有些旧代码从 Fortran 移植过来时使用 exp(iωt)存储的虚部是负值。如果你不检查时间因子直接换上负数虚部颗粒的“吸收”会变成“增益”能量不守恒。解决在参数入口统一做一次兼容注释写死时间因子并且如果检测到传入虚部小于 0给出警告而不是静默接受。自检方法是对吸收颗粒算 Q_abs如果 Q_ext - Q_sca 略微为负立刻怀疑符号。5.3 现象近场光强和远场光强对不上现象把观测距离 r 从 1 米减小到 0.1 毫米角分布的相对形状开始发生变化甚至出现前向强度下降这种不合直觉的结果。原因mie_intensity_farfield里的强度公式只保留远场渐进项。r 太小时颗粒内部和表面附近的消逝场、高阶多极子近场项还没衰减完远场公式天然失效。解决仿真前先验算 far-field 条件通常要求 r 远大于 d²/λ 且 r 远大于 λ。做实验室光路设计时我会把 r 至少取到 max(1000λ, d²/λ) 的十倍再和光学仿真软件的结果对比确认角分布不随 r 变化。5.4 现象用非偏振光公式去套线偏振实验现象实验用的是偏振激光探测器测到的角分布和后向散射比和计算值差一大截尤其是侧向 90° 附近。原因非偏振光的强度公式是 (|S1|² |S2|²)/2它等效于对偏振角 φ 做平均。但实际线偏振光散射强度与入射偏振方向有关I ∝ |S1|²sin²φ |S2|²cos²φφ 是偏振方向相对散射面的夹角。解决在功能里显式加入偏振角参数默认 φ0探测器只测同一偏振分量时按完整公式。如果不确定实验光路消偏程度就用手动旋转波片在探测器前测两个正交偏振分量再做平均比盲改折射率强。5.5 现象级数截断项数按经验取后向振荡对不上文献现象固定 n_stop 100算直径 10 微米银颗粒时前向匹配得不错后向却持续抖动和论文数据对不上。原因截断阶数必须随 x 增加而且高折射率或强吸收颗粒的收敛更慢。Wiscombe 的 x 4x^{1/3} 2 是对多数情况的保守估计但金属颗粒经常要额外加 20 到 30 项。解决把 n_stop 调成上一个经验值的上限并在系数累加时检查最后一项贡献。我一般会在循环里记录连续两个阶数对 S1 的相对增量当增量小于 1e-8 时提前退出如果到 n_stop 还没达到阈值自动把 n_stop 翻倍重跑。避免“看起来收敛、其实后向还没收”的假象。6. 进阶验证用光学定理、能量守恒和相函数归一化校准你的 Mie 计算6.1 光学定理消光截面和 S(0) 的前向关系前向散射幅度 S1(0) 的实部按光学定理直接决定消光截面。这个关系不依赖材料是麦克斯韦方程组的整体结果用来验证 an/bn 的实现非常合适。# 在 0° 方向计算 S1 s1_0, _ mie_s1_s2(an, bn, 1.0) k 2.0 * np.pi / wavelength_medium c_ext_from_s0 4.0 * np.pi / (k * k) * s1_0.real q_ext_from_s0 c_ext_from_s0 / (np.pi * (diameter / 2.0) ** 2) q_ext, q_sca mie_extinction_and_scattering(an, bn, x) print(Q_ext from coefficients:, q_ext) print(Q_ext from S(0): , q_ext_from_s0)如果两个 Q_ext 不能在小数后 6 位对上先检查 pi/tau 递推和 an/bn 公式的符号。这个自检比任何外部文献比对都有效因为它直接检查你代码内部的电磁场自洽性。6.2 无吸收颗粒的能量守恒无吸收颗粒的 Q_ext 和 Q_sca 必须精确相等。原因是散射本身不消耗能量只是把能量重新分配到各个方向。如果算出来 Q_ext 比 Q_sca 大一个可观量说明级数截断不够如果 Q_sca 比 Q_ext 大说明复折射率虚部符号错了。我把这一条放在 CI 脚本里只要当天改了递推代码就跑一次无吸收验证失败直接报错。6.3 散射相函数的归一化与后向散射比单颗粒 Mie 计算的最终交付物通常是相函数 P(θ)它满足 2π ∫ P(θ) sinθ dθ 1。用 s1/s2 构造相函数时分子是 (|S1|² |S2|²)/2分母是散射截面乘 k²。做完归一化后可以提取后向散射比 P(180°)/P(0°)这是激光雷达和气溶胶反演的重要参数。我的习惯是每次算完都把 x、Q_ext、Q_sca、g 因子和后向比打出来存成一个哈希表以后再跑同粒径参数时先比对旧值变化超过 1e-4 就停下来找原因。从入行到现在我所有 Mie 相关代码里都保留了这三条校验靠它们拦下过无数次因为符号约定和截断造成的翻车。希望帮到你。本文还有配套的精品资源点击获取