企业官网建设流程全解析

简介：这份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 时的特例，这时候只需保留 n=1 的项，散射光强随角度呈 (1+cos²θ) 分布，强度还强烈依赖 λ^-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（从 n=1 开始）。 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, derivative=True) yn_x = spherical_yn(idx, x) ynp_x = spherical_yn(idx, x, derivative=True) jn_mx = spherical_jn(idx, m * x) jnp_mx = spherical_jn(idx, m * x, derivative=True) # Riccati-Bessel 函数 psi = z*j_n，xi = 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, dtype=complex) bn = np.zeros(n_stop, dtype=complex) 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, derivative=True)一次性返回所有阶的导数，省去手写递推，代价是 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, I0=1.0, r=1.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), dtype=complex) s2 = np.zeros(len(theta_deg), dtype=complex) 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, I0=1.0, r=1000.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_by="number"): """ 计算对数正态分布颗粒群的平均远场散射光强。 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, I0=1.0, r=1000.0 ) # 对粒径做数值积分，并归一化到概率质量 numerator = np.trapz(I_matrix * pdf[:, None], d_grid, axis=0) 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 递推，并把两条路径在 x=200 附近对比 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 相关代码里都保留了这三条校验，靠它们拦下过无数次因为符号约定和截断造成的翻车。希望帮到你。

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