Mie散射理论与Python实现：从散射光强计算到粒径反演

简介：基于Mie理论的散射光强计算是光学与大气环境研究中的常见需求。面向需要模拟微小粒子散射行为的科研人员与高年级学生，这份MATLAB代码包可用于求解任意尺寸球形颗粒的散射光强、消光系数以及角度分布。压缩包为rar格式，共10个文件，包含9个m函数脚本和1个mat数据文件，整体大小仅34KB。其中核心Mie计算函数实现完整散射计算，配套脚本覆盖散射振幅、折射率定义、尺寸参数、消光系数与角度数据处理，并提供可视化绘图和测试脚本，另有1.06微米波长的mat数据可直接加载运行。目前已有575人学习下载，代码结构清晰、注释直接，用户只需输入颗粒尺寸、折射率与入射波长，即可得到散射光强分布及关键光学参数的可视化结果，对理解Mie理论、开展气溶胶或微粒光学特性研究具有实用价值。

1. 为什么散射光强不是“越细越暗”：Mie理论到底在算什么

一束激光打过去，悬浮液里粒径差零点几微米，散射光强的角分布就可能从“前后对称”变成“前向一柱擎天”。同一个粒子，波长换一下，曲线走势又完全不同。这背后就是标题里的Mie理论：给一个球形粒子的半径、折射率和入射波长，算出任意方向上的散射光强。它不是一条单调衰减曲线，而是一张随角度剧烈振荡的瓣状分布图。这篇笔记从适用边界讲到可复现的Python实现，再讲到参数影响和高频翻车点，最后给两招验证手段。适合正在做粒径反演、光学传感、大气光学仿真，又想把数值算明白的工程师。

2. 先判断要不要上Mie：尺寸参数、折射率与适用边界

Mie级数在理论上几乎适用于任意半径、任意折射率的球形粒子，但在工程上不是每个散射场景都值得跑一套完整级数。先判断问题落在哪个区间，再决定用什么工具，能省掉大量无效计算。

2.1 判断三段区的唯一标准：尺寸参数x

判断依据不是粒径绝对值，而是粒径和波长的比值。尺寸参数定义为 x=2πr/λ，r是粒子半径，λ是入射光在周围介质中的波长。这里的λ用的是介质中波长：粒子在空气中可按真空波长近似；粒子泡在水里，要把真空中波长除以水的折射率。x很小时，入射波“看”不到粒子细节，整个粒子像一个偶极子；x很大时，粒子相当于把波前裁成一个小圆孔，衍射和几何光学效应混在一起。

很多人误以为x大了就可以直接用几何光学。实际上，大粒子的散射光强角分布里仍包含明显的衍射峰，几何光学只描述反射和折射，算不出衍射峰的位置和宽度。所以哪怕是直径几十微米的粒子，要精确描述角分布，还是要回到Mie级数，只是截断阶数变得很大。我一般先算一次x，再看它落在哪个区间：x小于0.1直接用瑞利闭合式，大于0.1才进级数计算，这样能省掉不必要的循环。

2.2 输入三件套：半径、波长、相对折射率

跑通Mie计算只需要三个物理量：

半径。所有公式的尺寸参数x都用半径r。如果手里只有直径D，先除以2再代入。这个看似简单，但第5章会专门展开，因为不少现成代码接收的是直径，混用后曲线会整体错位。

波长。入参写λ0还是介质中波长λ，直接影响x。严格做法是 λ=λ0/n_medium，n_medium是周围介质的折射率。真实场景里，水中粒子的有效λ只有真空中的四分之三左右。

相对折射率。Mie公式里的m是粒子折射率相对周围介质的比值：m = n_particle / n_medium。材质有吸收时写成复数 m = n + iκ，没有吸收也要写成 n + 0j。举个例子：空气中水雾粒子 m=1.33+0j；水中气泡 m=1.0/1.33≈0.752+0j；水中聚苯乙烯微球 m=1.59/1.33≈1.195+0j。介质折射率这一项最常见被漏掉，直接导致曲线完全不能用。

2.3 该不该自己写：现有Python库与自实现的边界

成熟的Python库已经封装好了Mie级数，比如PyMieScatt和miepython，直接调用能快速拿结果。但要不要自己实现，取决于你的目标：

如果只是想快速画一条相函数曲线，直接调库。这个场景没有理由自己碰递推。如果要把Mie计算嵌入实时采集程序、需要在任意角度自由采样、或者系统需要S1和S2的复数振幅而不是归一化强度，建议自己封装一个无依赖版本。库函数为了通用性往往会做很多类型检查，实时处理时反而显得重。

我的做法是：第一次用成熟库对一遍结果，再用自己的实现去对照。自己写的代码不追求快，追求每一步都能打开看。下面这个实现只依赖numpy，逻辑照着可移植到C或嵌入式环境。

3. 从麦克斯韦方程组到散射光强：Python最小实现

散射光强的直接输出不是“强度”，而是两个复数振幅函数S1和S2。S1对应入射光偏振垂直于散射面的分量，S2对应平行分量。任意偏振态下的散射光强，都能由这两个振幅组合出来。所以我下面的代码先拿S1、S2，最后才做模平方，这样偏振信息不会丢。

3.1 先定截断阶数：Wiscombe公式与系数递推

Mie级数是无穷级数，实际只能截断到有限阶Nmax。截断阶数不够，高阶系数没收敛，大粒子曲线会整体掉精度。工程上最常用的经验公式是Wiscombe公式：Nmax = x + 4x^(1/3) + 2，再向下取整。小粒子时取下限15，避免x太小时Nmax算出来过小。下面这段直接按这个策略实现：

import numpy as np def mie_ab(x, m, nmax=None): """ 计算Mie散射系数 a_n, b_n。 x: 尺寸参数 2*pi*r/lambda m: 相对复折射率，必须写成复数 nmax: 截断阶数，默认用Wiscombe经验公式 返回 a, b，下标0对应n=1 """ if nmax is None: nmax = int(x + 4.0 * x ** (1.0 / 3.0) + 2.0) if nmax < 15: nmax = 15 z = m * x # 对数导数 D_n(z)，从高阶向下递推，避免向上递推在复平面发散 D = np.zeros(nmax + 2, dtype=complex) for n in range(nmax + 1, 0, -1): D[n - 1] = n / z - 1.0 / (D[n] + n / z) # psi_n(x) = x * j_n(x)，xi_n(x) = x * h_n^(1)(x) psi_x = np.zeros(nmax + 1, dtype=complex) xi_x = np.zeros(nmax + 1, dtype=complex) psi_x[0] = np.sin(x) xi_x[0] = np.sin(x) + 1j * np.cos(x) psi_x[1] = np.sin(x) / x - np.cos(x) xi_x[1] = (np.sin(x) / x - np.cos(x)) + 1j * (np.cos(x) / x + np.sin(x)) for n in range(2, nmax + 1): psi_x[n] = (2.0 * n - 1.0) / x * psi_x[n - 1] - psi_x[n - 2] xi_x[n] = (2.0 * n - 1.0) / x * xi_x[n - 1] - xi_x[n - 2] a = np.zeros(nmax + 1, dtype=complex) b = np.zeros(nmax + 1, dtype=complex) for n in range(1, nmax + 1): # 经典D函数形式的a_n, b_n ratio = D[n] / m + n / x a[n] = (ratio * psi_x[n] - psi_x[n - 1]) / (ratio * xi_x[n] - xi_x[n - 1]) ratio = m * D[n] + n / x b[n] = (ratio * psi_x[n] - psi_x[n - 1]) / (ratio * xi_x[n] - xi_x[n - 1]) return a[1:], b[1:]

系数递推是这套计算里最容易写错的地方。用D函数形式有两个好处：一是避开了直接在复自变量m*x上求球贝塞尔函数的导数，二是D从高阶向下递推，数值稳定性比向上递推好得多。向上递推在复数平面经常发散，尤其是折射率虚部接近零的时候，曲线会出现无缘无故的毛刺。

psi_x和xi_x的向上递推是稳定的，因为它们是实自变量x上的球贝塞尔函数组合，不会像D那样振荡放大。x到几百时这个递推也扛得住。a、b的下标我特意保留成物理阶次n，返回时去掉下标0，这样后续S1、S2的求和循环可读性更好。

Wiscombe公式里那个x^(1/3)不是拍脑袋来的，它对应的是高阶球贝塞尔函数趋近于零的临界阶数。x=100时Nmax大约是123，这个量级用上面的循环跑一万个角度也就几十毫秒，性能不是问题。真正的问题是有人把Nmax写死成20，后向散射的振荡细节全被截没了，后面第5章会专门说。

3.2 算S1、S2：角函数π_n、τ_n的递推与端点处理

有了a_n、b_n，还差角函数π_n和τ_n。它们由连带勒让德函数派生，工程里几乎没人直接调勒让德函数，都是用递推：

def pi_tau(mu, nmax): """ 计算角函数 pi_n, tau_n。 mu = cos(theta)，theta是散射角。 返回数组形状 (nmax+1, len(mu))，下标0不使用。 """ mu = np.atleast_1d(mu) pi = np.zeros((nmax + 1, mu.size), dtype=float) tau = np.zeros((nmax + 1, mu.size), dtype=float) # 常规递推 pi[1, :] = 1.0 for n in range(2, nmax + 1): pi[n, :] = ((2.0 * n - 1.0) / (n - 1.0)) * mu * pi[n - 1, :] \ - (n / (n - 1.0)) * pi[n - 2, :] for n in range(1, nmax + 1): tau[n, :] = n * mu * pi[n, :] - (n + 1.0) * pi[n - 1, :] # 端点解析值，处理 mu=±1 处的 0/0 退化 eps = 1e-7 for i, mui in enumerate(mu): if abs(mui - 1.0) < eps: for n in range(1, nmax + 1): pi[n, i] = n * (n + 1.0) / 2.0 tau[n, i] = n * (n + 1.0) / 2.0 elif abs(mui + 1.0) < eps: for n in range(1, nmax + 1): pi[n, i] = (-1) ** (n + 1) * n * (n + 1.0) / 2.0 tau[n, i] = (-1) ** n * n * (n + 1.0) / 2.0 return pi, tau

π_n、τ_n的递推在绝大多数角度上是稳定的，唯独θ=0和θ=π会退化。原因是它们本身含有除以sinθ的项，端点处变成0/0。前向和后向又恰好是激光雷达和粒径分析最关心的方向。我在代码里直接把端点解析值覆盖进去，这样输入角度范围包含0°或180°时也不会出NaN。

端点解析值本身不难记：前向θ=0时π_n=τ_n=n(n+1)/2；后向θ=π时多了个符号，奇偶阶相反。这个细节不处理，画出来的曲线两端会莫名其妙断掉，而且不同numpy版本行为还不一样，有的返回inf，有的直接报错。

3.3 从S1、S2到归一化散射光强

S1和S2的表达式是： S1(θ) = Σ_n (2n+1)/(n(n+1)) [a_n π_n + b_n τ_n] S2(θ) = Σ_n (2n+1)/(n(n+1)) [a_n τ_n + b_n π_n]

用代码写就是一层循环：

def mie_s12(a, b, mu): """ 由散射系数a, b计算振幅函数S1、S2。 mu: cos(theta)，可以是数组。 返回 s1, s2，复数数组，与mu同形状。 """ nmax = len(a) pi, tau = pi_tau(mu, nmax) s1 = np.zeros_like(pi, dtype=complex) s2 = np.zeros_like(pi, dtype=complex) for n in range(1, nmax + 1): fact = (2.0 * n + 1.0) / (n * (n + 1.0)) s1 += fact * (a[n - 1] * pi[n] + b[n - 1] * tau[n]) s2 += fact * (a[n - 1] * tau[n] + b[n - 1] * pi[n]) return s1, s2

拿到S1、S2之后，散射光强的完整表达式是：

I(θ) = (λ² / (4π²R²)) · I0 · (|S1|² + |S2|²) / 2

其中R是观测点到粒子的距离。如果只关心角度相对分布，直接画(|S1|²+|S2|²)就行，不用管前面的常数。入射光如果是自然光，总强度就是对两个偏振方向取平均；如果是线偏振光，还要看偏振面相对散射面的夹角，这就不只是平均能解决的了。

不依赖库、能看清每一步的实现，就是这个三件套：mie_ab拿系数，pi_tau拿角度函数，最后用一层循环凑S1、S2。下面这行是完整调用：

x = 2.0 * np.pi * 0.5 / 0.532 # 半径0.5微米水滴，波长532nm，空气中 a, b = mie_ab(x, 1.33 + 0j) mu = np.linspace(-1.0, 1.0, 501) s1, s2 = mie_s12(a, b, mu) I = np.abs(s1) ** 2 + np.abs(s2) ** 2

算出来的I在μ=1处有一个明显尖峰，那就是前向衍射峰，不是噪声。

4. 粒径、波长、折射率：三个参数如何改变光强角分布

同一套代码，把输入参数换一换，I(θ)的形态会有非常大的变化。这个变化不是单调的，理解它才能决定角度网格怎么布置、反演算法里该用哪一段数据。

4.1 尺寸参数x变大：从“椭圆”到“前向刺猬”

x=0.01时，角分布接近一个椭圆，前后向基本对称，这就是瑞利区。x=1时，前向开始略强，后向仍保留一定强度。x=10时，曲线上出现好几个次级峰，后向还有类似彩虹的鼓包结构。x=50以上，前向主峰变得非常窄，能量高度集中。

前向主峰的半宽大致是1/x弧度。x=100时半宽约0.57°，如果角度网格在[-1,1]上均匀撒500个点，整个前向峰只能拿到几个点，峰高会被严重低估。我一般先把x算出来，估算前向峰半宽，再单独加密前向区间的角度网格。

这里的振荡不是数值噪声，而是粒子内部驻波与外部散射波干涉的结果。粒径在微米级、激光单色性好的实验里，这些高频细节必须保留。有人为了“曲线平滑”把振荡抹掉，结果反演出来的粒径分布偏窄偏假。

4.2 折射率实部与虚部分别管什么

相对折射率的实部决定相位延迟，直接影响前向集中程度。实部从1.1加到1.6，前向峰会越来越尖，次级峰位置会移动。虚部则决定吸收，影响的是振荡振幅。虚部很小比如0.001时，角分布几乎和无吸收一样；虚部到0.05，振荡明显被压平；虚部到0.5，角分布基本变成一个光滑的前向鼓包，次级峰全部消失。

实际判断粒子是散射主导还是吸收主导，最有用的是单次散射反照率ω = Qsca / Qext。Qsca是散射效率，Qext是消光效率。ω接近1说明粒子几乎不吸收，光被弹开为主；ω小于0.5说明吸收主导，这时候靠多角度散射光强反演粒径的敏感性会大幅下降，应该转而测消光或吸收光谱。

4.3 光强角分布的实际用途：粒径反演与相函数

多角度散射光强是粒径反演的主要信息来源。前向区域对大粒径敏感，后向区域对小粒径更敏感。反演时建议做角度加权：前向峰区间用加密网格，后向用稀疏网格，而不是整个角度范围均匀取点。均匀取点会让反演权重大头落在前后向之间那些信息量不大的区域。

归一化后的角分布就是散射相函数，它是辐射传输方程的输入。做大气辐射、海洋光学、云雾多次散射仿真时，先用Mie算出单次散射相函数，再进入Monte Carlo或离散纵标法做传输计算。相函数的形状直接决定多次散射的结果，前向峰不准确，模拟出来的反射率会整体偏小。

还有一点值得注意：S1和S2的差异包含偏振信息。如果只记录(|S1|²+|S2|²)/2，退偏振比就丢了。做偏振激光雷达或偏振成像时，必须保留S1、S2的复数振幅，分别算两个偏振方向的强度。

5. 散射光强计算高频翻车点：现象、原因与修复

这套计算看着简单，实际跑起来翻车的地方非常集中。下面五条是我实践里遇到最多的，每一条都按“现象→原因→解决”写清楚。

5.1 半径和直径混用：曲线错位一个2倍

现象：算出的振荡峰位置整体右移，前向峰宽度几乎是文献值的一倍，换波长后偏差更明显。

原因：尺寸参数x=2πr/λ用的是半径r。部分现成实现或论文公式写成πD/λ，D是直径，两种口径差2倍。混用时x翻倍，所有峰位和峰宽都错位。

解决：在函数入口处统一用半径。我给每组参数单独写一行带注释的转换，比如：

radius_um = 0.5 # 半径，微米 diameter_um = 1.0 # 直径，微米，仅用于检查 x = 2.0 * np.pi * radius_um / wavelength_um

不接受“直径模式”的自动切换，显式写出半径来源，比隐式参数开关可靠得多。

5.2 介质折射率没归一化：水中泡算出假吸收

现象：水中气泡的角分布算出来像碳颗粒一样平滑，完全没有气泡该有的强振荡；更离谱的结果是消光效率大于散射效率。

原因：Mie公式里的m是相对折射率。周围介质是水时，m必须写1.0/1.33而不是1.0。如果把真空折射率直接当成周围介质，粒子和介质的对比度被夸大，虚部对比也被算错，看起来就像强吸收粒子。

解决：先写m = n_particle / n_medium，再进入Mie计算。同步改介质中波长λ=λ0/n_medium。这条规则对任何“粒子在液体里”的场景都适用，不只是气泡。

5.3 θ=0和θ=π处的角函数退化

现象：角度数组包含0°或180°时，输出出现NaN或inf，曲线两端断掉；角度数组不包含端点时一切正常。

原因：π_n、τ_n的定义里含有sinθ分母，端点处分子分母同时趋向零。numpy的浮点运算不会自动处理这种极限。

解决：在pi_tau函数里显式判断端点。前向θ=0时π_n=τ_n=n(n+1)/2；后向θ=π时π_n=(-1)^(n+1)·n(n+1)/2，τ_n=(-1)^n·n(n+1)/2。只要角度数组可能包含端点，就补上这段。

5.4 截断阶数写死导致振荡区丢失

现象：大粒径曲线整体数值偏低，前向峰高度不够，甚至出现Qsca > Qext这种违反能量守恒的结果。

原因：Nmax写死成某个小值，比如15或20。x增大后高阶项还没收敛就截断了，能量从高阶项漏掉。漏掉的不只是后向小抖动，连前向峰都会受影响。

解决：用Wiscombe公式Nmax = x + 4x^(1/3) + 2。x<1时取下限15。还要给对数导数D的数组多留几个缓冲位置，比如Nmax+15，防止高阶递推在复平面边缘抖动。我第3章的代码里D数组长度为nmax+2，实际生产环境建议改成nmax+15，多出来的开销可以忽略。

5.5 只取模不看偏振：偏振光路被平均掉

现象：做偏振激光雷达或偏振散射测量时，系统输出的退偏振比和实测对不上；用自然光场景的库函数去算偏振场景，结果完全不可用。

原因：总强度把S1和S2的模平方平均了。退偏振比依赖|S1|²和|S2|²的差异，平均以后这个差异被抹掉。

解决：在代码链路的最后一步之前不要取模。S1和S2始终保留复数，到需要出数的时候再分别算|S1|²、|S2|²。如果入射光是自然光，平均；如果是线偏振光，按偏振方向相对散射面的夹角重新组合。保留复数振幅这一习惯，能避免后期几乎所有的偏振返工。

6. 验证结果靠谱的两招：能量守恒与瑞利极限对齐

代码写完了，怎么确认它算得对？我每次都用两招做交叉验证，速度很快，能挡住绝大多数系数符号和递推方向错误。

第一招是能量守恒。用得到的a、b算消光效率和散射效率：

n_arr = np.arange(1, len(a) + 1, dtype=float) qext = 2.0 / x ** 2 * np.sum((2.0 * n_arr + 1.0) * np.real(a + b)) qsca = 2.0 / x ** 2 * np.sum((2.0 * n_arr + 1.0) * (np.abs(a) ** 2 + np.abs(b) ** 2))

物理上有两个硬指标：qext必须大于等于qsca；当x很大时，qext应趋近于2，这是大粒子消光效率的极限。如果qsca大于qext，说明a或b的实部符号错了，最常见的是D函数递推方向反了。如果qext在大x时明显偏离2，说明截断阶数不够。

第二招是瑞利极限对齐。取x=0.001，用任何折射率跑一遍，角分布应该与(1+cos²θ)成正比：

x_small = 0.001 a, b = mie_ab(x_small, 1.33 + 0j) mu = np.linspace(-1.0, 1.0, 101) s1, s2 = mie_s12(a, b, mu) I_mie = np.abs(s1) ** 2 + np.abs(s2) ** 2 I_ray = 1.0 + mu ** 2

归一化后两条曲线应该重合。如果形状对不上，基本可以断定是π_n递推或S1、S2组合公式有误。这一招对任何Mie实现都管用，因为它不依赖具体参数，只依赖物理极限。

我现在做粒径反演，上来先跑这两个检查：瑞利极限对齐保证角度部分没错，能量守恒保证级数截断没问题。两项都过了才敢把数据交给业务。曾经为了省时间把Nmax写死，结果换到2微米颗粒时整条曲线变形，浪费了一整天才定位到截断问题。后来把Wiscombe公式直接写进默认参数，再没犯过。希望帮到你。

本文还有配套的精品资源，点击获取