简介：这份MATLAB程序面向光学、大气科学与遥感领域的学习者与科研人员，提供经典米散射（MIE散射）理论的数值实现，用于计算球形粒子在光照下的散射系数、吸收系数、消光系数、后向散射系数及不对称因子等宏观光学参数。资源包共3个文件，均为m脚本，压缩包约1KB，分别承担主计算、散射矩阵系数求解与结果调用展示等职责，结构精简，便于直接运行与二次修改。程序不涉及粒子内部电磁场计算，聚焦散射与吸收等宏观性质，且已通过图形与文献结果比对验证，可靠性有保障。目前已有5713人学习下载，适合需要快速复现MIE理论、研究颗粒物光学特性或搭建遥感散射模型的读者参考使用。

1. 米散射计算工具：从气溶胶到颗粒检测的 MATLAB 实现

做光学颗粒测量或者大气辐射传输的开发者，大概率绕不开米散射这套公式。它描述的是球形粒子对电磁波的散射行为，输入折射率和尺度参数，输出散射系数、吸收系数、消光系数、后向散射系数和不对称因子。听起来像教科书里的东西，但真到项目里，比如激光粒度仪标定、气溶胶光学厚度反演、水体颗粒物后向散射建模，没有一套能跑、能改、结果可信的 MATLAB 程序，光靠查表或者调库会非常被动。这份资源就是一套完整的米散射 MATLAB 实现，把五个核心光学参量一次性算清楚，适合做光学测量、遥感反演和颗粒物检测的开发者直接拿来当计算内核用。下面按“能算什么 → 怎么跑 → 怎么改”的顺序拆一遍。

2. 核心算法拆解：Mie 级数怎么收敛出五个光学参量

2.1 输入参数与物理量的对应关系

米散射的输入其实就两组：粒子的复折射率 $m = n + ik$，以及尺度参数 $x = 2\pi r / \lambda$。折射率实部 $n$ 决定散射强弱，虚部 $k$ 决定吸收。尺度参数把粒径 $r$ 和波长 $\lambda$ 绑在一起，所以同一颗粒子在不同波长下表现完全不同。程序里通常把这两个量作为入口参数，再配合粒子数浓度或者路径长度，才能把“系数”从单粒子量转成宏观量。

常见做法是把单粒子的散射截面 $C_{sca}$、吸收截面 $C_{abs}$、消光截面 $C_{ext}$ 先算出来，再乘以数浓度 $N$ 得到系数。后向散射系数和不对称因子则分别对应角分布积分和相函数的一阶矩。理解这层映射，后面看代码就不会迷。

2.2 Mie 级数的截断与递推

米散射的解析解是一组无穷级数，实际计算必须截断。截断项数 $n_{max}$ 取小了结果不收敛，取大了白白耗算力。工程上常用经验公式：

% 计算 Mie 级数截断项数 % x 为尺度参数，经验公式来自 Bohren & Huffman function nmax = mie_nmax(x) nmax = round(2 + x + 4 * x^(1/3)); % 保证至少 5 项，避免小粒子时截断过少 nmax = max(nmax, 5); end

这个公式在 $x$ 从 0.1 到 100 的范围内都够用。逻辑说明：$x$ 很小时（瑞利区），几项就收敛；$x$ 大时（几何光学区），需要更多项。参数说明：x是尺度参数，nmax返回整数项数。我一般会在算完后检查一下最后两项的贡献，如果还大于总和的千分之一，就把nmax再放大 20% 重算，避免玄学误差。

2.3 五个光学参量的计算链路

核心计算围绕 Mie 系数 $a_n$ 和 $b_n$ 展开。下面这段是典型的实现骨架：

% 计算 Mie 系数 a_n, b_n % m 复折射率, x 尺度参数, nmax 截断项数 function [an, bn] = mie_ab(m, x, nmax) an = zeros(1, nmax); bn = zeros(1, nmax); % 用 Riccati-Bessel 函数的递推计算 psi, xi [psi, xi] = bessel_riccati(x, nmax); for n = 1:nmax % 尺寸参数与折射率组合 mx = m * x; [psi_mx, xi_mx] = bessel_riccati(mx, nmax); % Mie 系数标准公式 an(n) = (m*psi_mx(n)*psi(n) - psi(n)*psi_mx(n)) / ... (m*psi_mx(n)*xi(n) - xi(n)*psi_mx(n)); bn(n) = (psi_mx(n)*psi(n) - m*psi(n)*psi_mx(n)) / ... (psi_mx(n)*xi(n) - m*xi(n)*psi_mx(n)); end end

逻辑说明：bessel_riccati负责生成 Riccati-Bessel 函数 $\psi_n$ 和 $\xi_n$，这是 Mie 级数的基函数。参数说明：m是复折射率，x是尺度参数，nmax是截断项数。注意mx = m * x是复变量，MATLAB 对复数 Bessel 函数的支持没问题，但递推方向要选对，否则大 $x$ 时会数值溢出。我一般用向下递推算 $\psi$，向上递推算 $\xi$，这是血泪经验，反过来必翻车。

有了 $a_n$ 和 $b_n$，五个参量就是求和与积分的事：

% 由 Mie 系数计算五个光学参量 function [Qsca, Qabs, Qext, Qback, g] = mie_params(an, bn, x) n = 1:length(an); % 散射效率因子 Qsca = (2 / x^2) * sum((2*n + 1) .* (abs(an).^2 + abs(bn).^2)); % 消光效率因子 Qext = (2 / x^2) * sum((2*n + 1) .* real(an + bn)); % 吸收效率因子 = 消光 - 散射 Qabs = Qext - Qsca; % 后向散射效率因子 Qback = (1 / x^2) * abs(sum((2*n + 1) .* (-1).^n .* (an - bn))).^2; % 不对称因子 g = (4 / (x^2 * Qsca)) * sum( ... (n .* (n + 2) ./ (n + 1)) .* real(an .* conj(an(n+1) + bn .* conj(bn(n+1)))) ... + ((2*n + 1) ./ (n .* (n + 1))) .* real(an .* conj(bn)) ); end

逻辑说明：Qsca、Qext、Qabs是效率因子，乘以几何截面 $\pi r^2$ 就是截面。Qback是后向散射效率，注意那个 $(-1)^n$ 相位因子，漏了它后向散射会算成前向。g是不对称因子，范围 -1 到 1，0 表示各向同性，接近 1 表示强前向散射。参数说明：an、bn是 Mie 系数向量，x是尺度参数。这段代码里g的求和项容易写错，建议用已知结果（比如 $x=1$，$m=1.5$）对一下，不对称因子大概在 0.6 左右。

2.4 从单粒子到宏观系数

单粒子截面算完后，宏观系数就是乘以数浓度：

% 单粒子截面转宏观系数 % r 粒子半径(m), N 数浓度(m^-3) function [mu_s, mu_a, mu_e, mu_b] = macro_coeff(Qsca, Qabs, Qext, Qback, r, N) geom = pi * r^2; % 几何截面 mu_s = Qsca * geom * N; % 散射系数 mu_a = Qabs * geom * N; % 吸收系数 mu_e = Qext * geom * N; % 消光系数 mu_b = Qback * geom * N; % 后向散射系数 end

逻辑说明：geom是粒子的几何投影面积，N是数浓度。参数说明：r单位米，N单位每立方米，输出系数单位是每米。如果做大气气溶胶，N可能小到 $10^9$，做水体颗粒物可能到 $10^{15}$，量级跨度大，但公式不变。

3. 跑通程序：从环境准备到结果验证的完整步骤

3.1 环境与文件组织

这套程序是纯 MATLAB 实现，不依赖任何工具箱，基础版 MATLAB 就能跑。我一般把文件分成三块：mie_ab.m算系数，mie_params.m算参量，main.m做调用和画图。目录结构大概这样：

mie_scattering/ ├── main.m % 入口脚本 ├── mie_nmax.m % 截断项数 ├── mie_ab.m % Mie 系数 ├── mie_params.m % 五个参量 ├── bessel_riccati.m % Bessel 函数递推 └── macro_coeff.m % 宏观系数转换

如果你的 MATLAB 版本较老（比如 2016 之前），besselj和bessely对复变量的支持可能有问题，常见做法是自己写递推，不调内置函数。我一般会在bessel_riccati.m里用向下递推算 $\psi$，避免大 $x$ 时溢出。

3.2 一个可复现的调用示例

下面这段main.m可以直接抄：

% 米散射主脚本：计算并绘制五个参量随尺度参数的变化 clear; clc; % 固定复折射率，模拟水中的无机颗粒 m = 1.33 + 0.01i; % 尺度参数范围 x_list = logspace(-1, 2, 200); % 预分配 Qsca = zeros(size(x_list)); Qabs = zeros(size(x_list)); Qext = zeros(size(x_list)); Qback = zeros(size(x_list)); g = zeros(size(x_list)); for k = 1:length(x_list) x = x_list(k); nmax = mie_nmax(x); [an, bn] = mie_ab(m, x, nmax); [Qsca(k), Qabs(k), Qext(k), Qback(k), g(k)] = ... mie_params(an, bn, x); end % 画图 figure; loglog(x_list, Qsca, 'b-', x_list, Qabs, 'r--', x_list, Qext, 'k-.'); xlabel('尺度参数 x'); ylabel('效率因子'); legend('Qsca','Qabs','Qext'); grid on; figure; semilogx(x_list, g, 'm-'); xlabel('尺度参数 x'); ylabel('不对称因子 g'); grid on;

逻辑说明：logspace生成对数分布的尺度参数，覆盖瑞利区到几何光学区。参数说明：m是复折射率，实部 1.33 对应水的折射率，虚部 0.01 表示弱吸收。跑完后你会看到Qext在小 $x$ 时按 $x^4$ 增长（瑞利散射），大 $x$ 时趋近 2（消光佯谬），这是验证程序对不对的硬指标。

3.3 结果验证：三个必须对的基准

跑通不等于算对。我一般会做三个检查：第一，小 $x$ 时Qsca应该正比于 $x^4$，在双对数图上是斜率 4 的直线；第二，大 $x$ 时Qext趋近 2，Qabs趋近 0（如果 $k=0$）；第三，g在 $x$ 很小时趋近 0，大 $x$ 时趋近 1。如果这三个都对，基本可以放心用。

提示：如果Qext在大 $x$ 时震荡剧烈且不收敛，八成是nmax不够，把mie_nmax里的系数从 4 调到 6 再试。

3.4 性能与向量化

上面的main.m是循环调用，200 个点跑起来不到一秒。但如果要做反演，比如在优化循环里反复调用，就得向量化。常见做法是把mie_ab里的for n = 1:nmax改成矩阵运算，或者用arrayfun并行。我一般会先跑通循环版，确认结果对，再改向量化，不然调试起来是黑匣子。

4. 避坑与排查：米散射计算里最容易翻车的五个点

4.1 现象：后向散射系数算出来是负的

原因：Qback公式里漏了abs(...).^2，或者(-1).^n写成了(-1)^n在向量化时出错。解决：检查Qback那行，确保先求和再取模平方，相位因子用(-1).^n逐项乘。

4.2 现象：吸收系数大于消光系数

原因：Qabs = Qext - Qsca里Qext和Qsca的截断项数不一致，或者Qsca的求和漏了(2n+1)因子。解决：确保两个求和用同一个nmax，并且Qsca的权重是(2*n+1)。

4.3 现象：大尺度参数时结果震荡不收敛

原因：nmax取小了，或者 Bessel 函数递推方向反了导致数值溢出。解决：用mie_nmax的经验公式，并且 $\psi$ 用向下递推，$\xi$ 用向上递推。如果还不行，把nmax放大 50% 看结果是否稳定。

4.4 现象：复折射率虚部为负时结果异常

原因：物理上虚部应该非负（表示吸收），如果输入负值，Qabs会变负。解决：在入口加一个检查if imag(m) < 0, error('虚部必须非负'); end，避免后面算出鬼结果。

4.5 现象：小粒子时Qsca不按 $x^4$ 增长

原因：nmax最小项设得太大，比如强制至少 10 项，小 $x$ 时高阶项贡献了虚假散射。解决：mie_nmax里max(nmax, 5)改成max(nmax, 3)，或者对小 $x$ 单独用瑞利近似。

5. 二次开发：把这套内核接进反演和测量流程

这套程序最值钱的地方不是算五个参量，而是它能当计算内核嵌到更大的流程里。我一般会做三件事：第一，把mie_params包装成接受向量输入的函数，方便在优化算法里调用；第二，加一个正向模型，输入粒径分布和折射率，输出光谱消光，用来和实测对比；第三，用fminsearch或lsqnonlin反演粒径和折射率。

下面是一个反演骨架：

% 用实测消光光谱反演粒径和折射率实部 % 假设粒径服从对数正态分布 function cost = inverse_cost(params, lambda, Qext_meas) r_med = params(1); % 中值半径 sigma = params(2); % 几何标准差 n_real = params(3); % 折射率实部 m = n_real + 0.001i; % 固定弱吸收 Qext_model = zeros(size(lambda)); for k = 1:length(lambda) x = 2 * pi * r_med / lambda(k); nmax = mie_nmax(x); [an, bn] = mie_ab(m, x, nmax); [~, ~, Qext_model(k), ~, ~] = mie_params(an, bn, x); end cost = sum((Qext_model - Qext_meas).^2); end

逻辑说明：params是待反演参数，lambda是波长向量，Qext_meas是实测消光。参数说明：r_med中值半径，sigma几何标准差，n_real折射率实部。实际用的时候，Qext_meas得先归一化，不然量级差异会让优化跑偏。

改造时注意两点：一是反演里mie_ab会被调用几千次，必须向量化或者用parfor；二是折射率虚部如果也反演，参数相关性会很强，常见做法是固定虚部或者用多波长联合反演。我一般会先用合成数据验证反演流程，确认能收敛到真值，再上实测数据。从那以后我每次改完mie_ab都强制走一遍合成数据验证，不然一个符号错误就能让反演结果偏到姥姥家。希望帮到你。

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