简介一套基于Mie散射理论、在MATLAB环境下运行的光学参数计算程序包适用于大气气溶胶、云滴、纳米颗粒等球形粒子的散射、吸收与消光特性分析。程序依据H.A. Mie于1908年提出的经典理论通过输入颗粒半径、复折射率和入射光波长即可快速求解Mie散射系数、消光系数、吸收系数及不同角度的散射光强分布从而支持对颗粒光学特性的快速评估。资源包共15个文件以14个MATLAB脚本.m文件为主体另含1个zip压缩包整体仅13KB。脚本功能模块覆盖吸收、消光/散射效率、散射振幅函数、复折射率、雾衰减计算并附有可直接运行的示例脚本使用者可基于这些模块快速搭建自己的散射计算流程。已有510人学习查看适合大气科学、光学工程、环境监测等专业的学生与研究者作为入门及进阶的实用工具。 上周整理旧工程时我又看到了那个叫 broken2t1 的文件夹。它记录了我调试米氏散射Mie scattering时最狼狈的一段时间明明折射率只是加了很小一点虚部算出来的散射系数和消光系数却直接变成了 NaN。后来才知道那不是物理出了问题而是递推算法在强吸收条件下踩了数值地雷。如果你也在做球形颗粒的光吸收、散射或消光计算相信迟早会碰到类似情况。这篇就把我当时的排查过程、几个系数的物理含义、稳定数值算法的改造方法以及从系数换算到实际样品透过率的套路一次性讲清楚。文章适合正在写光学仿真、做纳米材料表征、算大气颗粒物散射或只是想搞懂 Mie 公式到底怎么落地的人。1. broken2t1 这个版本标签背后是一场数值事故1.1 从旧工程里翻出来的坏结果那个工程里我想要算的是某批球形颗粒的消光光谱。颗粒折射率实部大概在 2.0 附近吸收不能忽略所以折射率要写成复数的形式[ m n i k ]我当时顺手把工作分支命名为 broken2t1。含义很简单实部取 2.0同时第一次把吸收项 t1 加进去结果代码就 break 了。崩掉的位置很有意思。颗粒粒径不算大尺寸参数不过个位数按理说 Mie 级数收敛很快。但程序算到大约第 40 多项时an和bn的值开始抖动随后直接变成nan。第一反应是公式抄错了、或者索引越界了但把无吸收情况拿回来一跑又完全正常。只要 k 不为 0同样一套代码就崩。这个问题非常典型吸收被引入后折射率变成复数用于计算内部场的球 Bessel 函数不再“温和”递推过程会出现灾难性的数值放大。Mie 理论本身没有错错的是实现方式。1.2 “2t1”背后的物理场景先别急着改代码得搞清楚这个 m2ik 到底对应什么情况。真空或空气中的普通介质折射率实部通常在 1.0 到 1.5 之间。但是很多实际颗粒并非如此半导体纳米颗粒、高折射率陶瓷粉体、某些聚合物微球实部可以到 1.8 甚至 2.5再叠加吸收k 可能是 0.01、0.1甚至更高。在这个区间里Rayleigh 近似已经不够用了但颗粒又没有大到能直接用几何光学所以必须完整求解 Mie 理论与 Maxwell 方程。于是你既要处理高阶项又要面对复折射率带来的病态计算。2. 吸收、散射、消光三个系数到底在算什么2.1 Qext、Qsca、Qabs 的关系Mie 散射计算的核心输出有三个无量纲效率因子[ Q_{\rm ext} \frac{2}{x^2}\sum_{n1}^{\infty}(2n1){\rm Re}(a_nb_n) ][ Q_{\rm sca} \frac{2}{x^2}\sum_{n1}^{\infty}(2n1)(|a_n|^2|b_n|^2) ][ Q_{\rm abs} Q_{\rm ext}-Q_{\rm sca} ]很多工程文档习惯把 Qext 叫消光系数把 Qsca 叫散射系数把 Qabs 叫吸收系数。严格说它们不是“系数”而是单个颗粒的消光、散射、吸收效率以颗粒的几何投影面积为基准。最终的实际消光截面是[ C_{\rm ext} Q_{\rm ext}\cdot \pi\left(\frac{D}{2}\right)^2 ]这三个量的物理逻辑非常直观一束光打到颗粒上一部分被颗粒吸收转成热量一部分被重新散射到其他方向加在一起就是入射光被“消掉”的总量。所以 Qext 永远等于 Qsca 加 Qabs这不是近似而是能量守恒的直接结果。2.2 an 和 bn 代表什么级数里的 an 和 bn 是米氏散射系数。an 对应电多极项bn 对应磁多极项。n1 是偶极项n2 是四极项n3 是八极项依此类推。当颗粒尺寸远小于波长也就是尺寸参数 (x2\pi r/\lambda\ll 1) 时只需要保留 n1 这一项这就是 Rayleigh 散射极限。颗粒变大后高阶项逐渐不可忽略散射光的前向分量增强吸收和散射的相对占比也会跟着改变。我个人的理解方式是an、bn 本质上描述了颗粒内部场和外部入射场的“共振匹配程度”。实部决定振荡相位虚部代表吸收损耗。折射率虚部一旦变大内部场的衰减就快计算时那些带 m 指数的函数会剧烈变化数值上稍不小心就失真。3. 强吸收下的失稳复折射率如何击穿向上递推3.1 灾难性抵消才是元凶Mie 计算里散射系数需要用到 Riccati-Bessel 函数[ \psi_n(z)z j_n(z), \qquad \xi_n(z)z h_n^{(1)}(z) ]对实数 x 而言(\psi_n(x)) 可以用简单的向上递推[ \psi_{n1}(x)\frac{2n1}{x}\psi_n(x)-\psi_{n-1}(x) ]这对实数参数非常稳定。问题出在把 zmx 带入复平面之后。当 m 有虚部时(\psi_n(mx)) 的振荡幅度会随 n 指数级增长但计算机里的双精度浮点数有上限。继续用同一套向上递推很快就会溢出。更隐蔽的问题是灾难性抵消即使没溢出如果计算过程中出现两个非常大的数相减得到一个很小的数那么有效位数会被大量吃掉。这就像用两个亿级数字相减去查一份几块钱的零头结果当然不可靠。我在 broken2t1 分支里遇到的 NaN正是这么来的。3.2 稳定化思路不直接算函数改算比值解决思路并不神秘既然 (\psi_n(mx)) 本身会爆炸那就不算它改算更温和的对数导数[ D_n(z) \frac{\psi_n(z)}{\psi_n(z)} ]这个比值在复平面里相对可控。然后用它重建 an 和 bn。要实现 D_n必须采用向下递推而不是向上递推。向下递推的公式是[ D_{n-1}(z)\frac{n}{z}-\frac{1}{D_n(z)n/z} ]从足够大的 n 处设一个初值逐步回推。这样每步都在做除法而不是大数相减复折射率带来的指数增长就被抑制住了。这是 Bohren 和 Huffman 在《Absorption and Scattering of Light by Small Particles》里给出的经典稳定方案今天绝大多数现代 Mie 代码还在沿用。4. 稳定版算法实现与 m20.1i 案例验证4.1 截断阶数 nmax 怎么选Mie 级数理论上要加到无穷大实际计算肯定要截断。截断阶数 nmax 的选择直接影响精度和速度。经验公式是[ n_{\max} \approx x 4.05 x^{1/3} 2 ]这个公式在尺寸参数 x 从很小到几百都很可靠。但对于复折射率尤其是虚部较大的强吸收颗粒我会建议再往后多推几十项。向下递推也需要一个“起跑距离”否则初始值误差还没衰减就进入了目标区间。下面是完整可跑的稳定版计算函数只依赖 numpyimport numpy as np def logder_D(nmax, z): 对数导数 D_n(z)psi_n(z)/psi_n(z)向下递推。 返回下标从 0 到 nmax 的完整数组。 D np.zeros(nmax 1, dtypecomplex) for n in range(nmax, 0, -1): D[n - 1] n / z - 1.0 / (D[n] n / z) return D def riccati_psi_chi(nmax, x): 实数参数 x 下 Riccati-Bessel 函数 psi_n(x) 和 chi_n(x) 使用稳定的向上递推。 psi np.zeros(nmax 1) chi np.zeros(nmax 1) psi[0] np.sin(x) chi[0] -np.cos(x) if nmax 1: psi[1] np.sin(x) / x - np.cos(x) chi[1] -np.cos(x) / x - np.sin(x) for n in range(1, nmax): psi[n 1] (2 * n 1) / x * psi[n] - psi[n - 1] chi[n 1] (2 * n 1) / x * chi[n] - chi[n - 1] return psi, chi def mie_q(m, x): 返回 (Qext, Qsca, Qabs) m: 颗粒折射率 / 介质折射率复数 x: 尺寸参数 2*pi*r/lambda nstop int(x 4.05 * x ** (1 / 3) 2) nmax nstop 30 D_mx logder_D(nmax, m * x) psi, chi riccati_psi_chi(nmax, x) xi psi 1j * chi qext_sum 0.0 qsca_sum 0.0 for n in range(1, nstop 1): A D_mx[n] / m n / x B m * D_mx[n] n / x an (A * psi[n] - psi[n - 1]) / (A * xi[n] - xi[n - 1]) bn (B * psi[n] - psi[n - 1]) / (B * xi[n] - xi[n - 1]) qext_sum (2 * n 1) * (an bn).real qsca_sum (2 * n 1) * (abs(an) ** 2 abs(bn) ** 2) Qext 2.0 / (x ** 2) * qext_sum Qsca 2.0 / (x ** 2) * qsca_sum Qabs Qext - Qsca return Qext, Qsca, Qabs代码核心就是用向下递推算出 D_n(mx)用它替换掉 an、bn 原始表达式中的危险项。psi 和 chi 仍然按实数递推因为这里的宗量是实数尺寸参数 x没有溢出问题。4.2 用 m20.1i 跑出来的趋势以波长 550nm 为例介质是空气颗粒折射率为 (2.00.1i)分别算几个典型粒径。尺寸参数和主要趋势如下表颗粒直径 D尺寸参数 x主要表现50nm0.29Rayleigh 区吸收占绝对主导200nm1.14吸收仍然高于散射高阶项开始出现500nm2.86散射显著增强前向散射峰形成1000nm5.71散射接近甚至超过吸收衍射效应突出这个趋势是判断结果是否合理的很好的校准线。如果代码在 x 很小的时候给不出“吸收主导”或者在 x 很大的时候给不出“散射上升”大概率是某个细节写错了。一个常用的校验方法是拿 Rayleigh 极限公式对比[ Q_{\rm abs} \approx 8x,{\rm Im}\left(\frac{m^2-1}{m^22}\right) ]当 x 很小时这个近似解和完整 Mie 代码应当非常接近。我在实际项目中都是先跑这个极限确认通过后再上千纳米级的完整计算。5. 消光系数换算到实际样品从单个颗粒到宏观透过率5.1 Lambert-Beer 关系怎么用代码算出来的是每个颗粒的效率因子但实验测量通常拿到的是悬浮液或粉末层的透过率 T。两者之间通过颗粒数浓度 N 和光程 L 连接[ T \exp(-N L C_{\rm ext}) ]其中消光截面 (C_{\rm ext}Q_{\rm ext}\pi(D/2)^2)。如果样品里颗粒尺寸有分布还要对粒径分布做积分。很多卖仪器的测试报告只给“消光值”如果你自己要做波长扫描最好把 Qext 光谱自己算一遍。这个公式用起来有一个常见的坑颗粒数浓度 N 的单位。质量浓度 1g/L 的球形颗粒数浓度要先用密度换算成体积再除以单个颗粒体积。这一步算错后面所有绝对量都会偏好几个数量级。5.2 从消光光谱能看出什么我的经验是别只看 Qext 单条曲线要把吸收、散射、消光三条曲线同时放在一张图里。这样能立刻判断某个波段的消光到底是吸收贡献还是散射贡献。举个例子如果你在可见光区看到一个宽谱消光峰同时 Qabs 明显大于 Qsca说明材料以吸收为主这通常是带隙吸收或等离激元吸收的特征。反过来如果 Qsca 在短波段快速上升而 Qabs 很弱则是典型的散射主导常见于高折射率透明颗粒。另外单位体积内的散射强度还和颗粒浓度有关。我曾经遇到过“消光峰位置看着没错但峰强总是偏低”的情况最后发现是粒径分布没考虑进去。Mie 计算对不同粒径非常敏感粒径相差 20%消光光谱的形状就会明显变化。做粒径分布反演时至少要覆盖 3 到 5 个粒径点再插值。5.3 实际工程里的几条经验最后分享几条从 broken2t1 这个坑里总结出来的经验。第一折射率 m 是颗粒相对介质的不是颗粒绝对折射率。颗粒在水中和颗粒在空气中的 m 完全不同m20.1i 这个值是相对空气的。如果你把介质折射率也代成 1.5一定搞混。第二双精度浮点不是万能的。遇到折射率虚部特别大比如 k1 的强吸收颗粒常规实现还是可能出问题。这时可以先试试增加 nmax 余量。如果还是不稳定就需要用更高精度的递推策略或者调用成熟开源库比如 miepython而不是自己硬造轮子。第三结果里如果出现 Qext 小于 Qsca或者 Qabs 为负那就是数值误差已经大到不能用了。先检查 x 是否用了直径而不是半径再检查折射率是否写成 2.0 而不是 2.00.0j。这些低级错误我在实际项目里都见过不止一次。米氏散射计算看起来是一组标准公式但真正落地时数值稳定性、截断阶数、物理量纲这些细节都会决定结果是否可信。希望这篇能帮你少踩一次 broken2t1 式的坑。本文还有配套的精品资源点击获取