MATLAB米散射计算程序：直接输出五个光学系数

简介面向光学、大气科学与遥感方向的学习者与研究人员这份资源提供经典米散射MIE散射理论的MATLAB数值实现用于计算球形粒子在光照下的散射系数、吸收系数、消光系数、后向散射系数及不对称因子适合需要快速复现米散射计算、验证模型或开展颗粒物光学特性分析的中级编程用户。压缩包共3个文件均为m脚本整体约1KB分别承担主计算、散射矩阵系数求解与驱动调用展示的职责结构精简、便于按需修改。目前已有5732人学习下载说明其在教学与科研场景中具备一定认可度。程序侧重宏观散射与吸收性质的计算不含粒子内部电磁场模块且结果已通过图形与文献比对验证读者可据此获得一套可靠的基础计算框架并在此基础上拓展到大气颗粒物、遥感成像散射效应及微纳米材料光学特性等应用。1. 米散射计算工具一份能直接跑出四个系数的 MATLAB 程序做颗粒物光学测量或者大气辐射传输的同行大概都遇到过这个场景实验台上刚测完一组气溶胶样品的透射率导师问你要消光系数和不对称因子你翻出半年前写的散射代码发现里面角度积分步长和折射率虚部符号又得重新推一遍。米散射MIE散射的公式本身不算复杂难的是把贝塞尔函数、递推关系、收敛判据这些细节一次性写对还得保证在粒径参数跨越好几个量级时不发散。这份 MATLAB 程序解决的就是这个痛点——输入复折射率和粒径参数直接吐出散射系数、吸收系数、消光系数、后向散射系数和不对称因子五个量。它适合两类人一是做气溶胶、云物理、水体颗粒物光学特性分析的研究生二是需要在光学检测设备标定环节快速验证理论值的工程开发。整个程序是纯 MATLAB 实现不依赖任何工具箱拿到手改几个参数就能跑。2. 核心算法拆解从 Mie 级数到五个系数的计算链路2.1 为什么选 Mie 级数展开而不是近似公式米散射的严格解来自麦克斯韦方程组在球坐标系下的分离变量最终落到两个散射振幅函数 S1(θ) 和 S2(θ) 的无穷级数上。程序里没有用瑞利近似或者几何光学近似原因很直接当粒径参数 x 2πr/λ 落在 0.1 到 100 这个区间时近似公式的误差能到百分之几十而 Mie 级数在截断项数足够时精度只受浮点运算限制。程序的核心就是构造两个 Riccati-Bessel 函数 ψ_n(x) 和 ξ_n(x)然后通过递推关系算出 Mie 系数 a_n 和 b_n。function [a_n, b_n] mie_coefficients(m, x, n_max) % m: 复折射率 (m n i*k注意虚部符号约定) % x: 粒径参数 x 2*pi*r/lambda % n_max: 级数截断项数 % 初始化递推所需的中间量 nu n_max 1; y x * m; % 用向下递推计算 Riccati-Bessel 函数 psi_n(x) % 避免向上递推在 x 较小时产生的数值溢出 psi zeros(1, nu); psi(nu) 0; psi(nu-1) 1e-30; % 初始种子 for n nu-1:-1:2 psi(n-1) (2*n-1)/x * psi(n) - psi(n1); end % 归一化消除种子值带来的比例因子 psi psi / psi(1); % 计算 xi_n(x) psi_n(x) - i*chi_n(x) % chi_n 用向上递推因为它在 x 小时增长缓慢 chi zeros(1, nu); chi(1) cos(x); chi(2) sin(x) cos(x)/x; for n 2:nu-1 chi(n1) (2*n-1)/x * chi(n) - chi(n-1); end xi psi - 1i * chi; % 计算对数导数 D_n(mx)这是 a_n、b_n 的关键中间量 D zeros(1, nu); D(nu) 0; % 远场截断 for n nu-1:-1:1 D(n) 1 / ((2*n1)/y - D(n1)); end % 组装 Mie 系数 a_n zeros(1, n_max); b_n zeros(1, n_max); for n 1:n_max a_n(n) (D(n)/m n/x) * psi(n) - psi(n-1); a_n(n) a_n(n) / ((D(n)/m n/x) * xi(n) - xi(n-1)); b_n(n) (m*D(n) n/x) * psi(n) - psi(n-1); b_n(n) b_n(n) / ((m*D(n) n/x) * xi(n) - xi(n-1)); end end这段代码里有两个容易翻车的点。第一是 ψ_n 的递推方向当 x 小于 1 时向上递推会让 ψ_n 迅速衰减到浮点下溢所以必须用向下递推加归一化代码里 psi(nu-1) 1e-30 这个种子值就是干这个用的。第二是 D_n(mx) 的递推这里用的是向下递推因为 m 的虚部为正时 y 落在复平面的上半平面向下递推是稳定的。参数 n_max 的选取不能拍脑袋经验公式是 n_max ceil(x 4*x^(1/3) 2)程序里我一般会再乘个 1.2 的余量防止在共振区截断误差偏大。2.2 五个系数的物理含义与计算顺序Mie 系数算出来之后散射系数 Qsca、消光系数 Qext、吸收系数 Qabs 就是级数的加权求和。程序里的计算顺序是 Qsca 和 Qext 先算Qabs 用差值得到这样能保证能量守恒关系 Qext Qsca Qabs 在数值上严格成立。function [Qsca, Qext, Qabs, Qback, g] mie_efficiencies(a_n, b_n, x) % 输入 Mie 系数和粒径参数输出五个无量纲效率因子 n_max length(a_n); n 1:n_max; % 散射效率因子|a_n|^2 |b_n|^2 的加权和 Qsca 2 / x^2 * sum((2*n 1) .* (abs(a_n).^2 abs(b_n).^2)); % 消光效率因子实部 a_n b_n 的加权和 Qext 2 / x^2 * sum((2*n 1) .* real(a_n b_n)); % 吸收效率因子用差值保证能量守恒 Qabs Qext - Qsca; % 后向散射效率因子奇数项和偶数项交替相减 Qback 1 / x^2 * abs(sum((2*n 1) .* (-1).^n .* (a_n - b_n)))^2; % 不对称因子 g散射相函数的一阶矩 % 用 a_n、b_n 的交叉项构造 g_num 0; for n 1:n_max-1 term1 (n*(n2)/(n1)) * real(a_n(n)*conj(a_n(n1)) b_n(n)*conj(b_n(n1))); term2 ((2*n1)/(n*(n1))) * real(a_n(n)*conj(b_n(n))); g_num g_num term1 term2; end g 4 / (x^2 * Qsca) * g_num; endQback 的计算里那个 (-1)^n 交替符号是后向散射特有的因为后向方向对应散射角 θ πLegendre 多项式在 θ π 处的值就是 (-1)^n。不对称因子 g 的公式看起来复杂本质上是 cosθ 在散射相函数上的加权平均g 接近 1 表示强前向散射接近 0 表示各向同性接近 -1 表示强后向散射。程序里 g 的计算用了交叉项求和比直接对相函数做数值积分快一个量级精度也更好因为避开了角度离散化带来的截断误差。2.3 主程序的组织方式与参数入口整个项目的主入口是一个脚本文件把折射率、粒径、波长三个物理量放在最前面后面依次调用系数计算和效率因子计算最后把结果打印到命令行并画一张散射相函数图。这种组织方式的好处是改参数不用翻函数文件坏处是脚本里的变量会污染工作区所以我一般会在脚本开头加 clear 和 close all。% main_mie.m —— 米散射计算主脚本 clear; clc; close all; % 物理参数输入区 lambda 550e-9; % 入射波长单位米550nm 可见光 r 0.5e-6; % 颗粒半径单位米0.5微米 m_real 1.5; % 折射率实部 m_imag 0.01; % 折射率虚部吸收性颗粒取正值 m m_real 1i * m_imag; % 复折射率 % 粒径参数与截断项数 x 2 * pi * r / lambda; n_max ceil(x 4 * x^(1/3) 2); n_max ceil(1.2 * n_max); % 留 20% 余量 % 调用核心计算 [a_n, b_n] mie_coefficients(m, x, n_max); [Qsca, Qext, Qabs, Qback, g] mie_efficiencies(a_n, b_n, x); % 结果输出 fprintf(粒径参数 x %.4f\n, x); fprintf(散射系数 Qsca %.6f\n, Qsca); fprintf(消光系数 Qext %.6f\n, Qext); fprintf(吸收系数 Qabs %.6f\n, Qabs); fprintf(后向散射系数 Qback %.6f\n, Qback); fprintf(不对称因子 g %.6f\n, g);参数入口的设计上波长和半径都用国际单位制这样 x 的计算不会出现单位混用的问题。折射率虚部的符号约定是个高频翻车点物理上吸收性介质的复折射率写成 m n - ik但很多数值程序里为了递推稳定写成 m n ik这份程序用的是后者所以输入吸收性颗粒时虚部取正值。如果拿到的实验数据虚部是负号约定记得在赋值时翻一下否则算出来的 Qabs 会是负数那就闹笑话了。3. 运行环境与实操步骤从零跑通第一组数据3.1 MATLAB 版本要求与文件放置这份程序对 MATLAB 版本的要求很低R2016a 之后的版本都能跑因为用到的语法只有基本的矩阵运算、匿名函数和 fprintf没有用 string 数组或者 arguments 块这些新特性。如果手头是 Octave 也能跑但要注意 Octave 里 1i 虚数单位的写法兼容性没问题fprintf 的格式也一致唯一需要留意的是 Octave 对中文注释的编码偶尔会出问题建议把文件存成 UTF-8 无 BOM 格式。文件放置上整个项目就三个文件主脚本 main_mie.m、系数计算函数 mie_coefficients.m、效率因子计算函数 mie_efficiencies.m。三个文件放在同一个文件夹里在 MATLAB 里把当前工作目录切过去直接运行 main_mie 就行。不需要 addpath不需要安装任何工具箱这一点比那些依赖 Optimization Toolbox 或者 Symbolic Math Toolbox 的散射代码省心不少。3.2 第一组验证数据的选取跑通之后第一件事不是急着换自己的参数而是拿一组已知结果做验证。米散射领域最常用的验证点是 x 1、m 1.5 0.01i 这组参数文献里能查到对应的 Qext 大约在 2.3 附近。程序跑出来如果 Qext 在 2.2 到 2.4 之间说明系数计算和级数截断都没问题。如果偏差很大先检查 n_max 够不够再检查折射率虚部符号。# 在 MATLAB 命令行里直接调用不跑脚本 x 1; m 1.5 0.01i; n_max ceil(1.2 * (x 4*x^(1/3) 2)); [a,b] mie_coefficients(m, x, n_max); [Qs, Qe, Qa, Qb, g] mie_efficiencies(a, b, x); fprintf(Qext %.4f, Qsca %.4f, g %.4f\n, Qe, Qs, g);这组命令的好处是绕开了主脚本里的参数区直接在命令行里改 x 和 m适合快速扫参数。验证通过之后再回到主脚本把波长和半径换成自己实验对应的值。这里有个血泪经验换参数时一定要同步改 n_max程序里虽然有自动计算但如果你手动写死了 n_max 又忘了改大粒径情况下级数截断不够Qext 会偏小而且偏得很隐蔽因为曲线看起来还是光滑的。3.3 批量扫描与结果导出单点计算跑通之后实际工作里更常见的是扫粒径或者扫波长。程序本身没带批量循环但加一层 for 循环就行注意预分配数组不然 MATLAB 的编辑器会给你标黄警告。% 扫粒径半径从 0.05 微米到 5 微米50 个点 r_list linspace(0.05e-6, 5e-6, 50); Qext_list zeros(size(r_list)); g_list zeros(size(r_list)); for k 1:length(r_list) x 2 * pi * r_list(k) / lambda; n_max ceil(1.2 * (x 4*x^(1/3) 2)); [a_n, b_n] mie_coefficients(m, x, n_max); [~, Qext_list(k), ~, ~, g_list(k)] mie_efficiencies(a_n, b_n, x); end % 导出为 csv方便后续在别的工具里画图 T table(r_list, Qext_list, g_list, ... VariableNames, {radius_m, Qext, g}); writetable(T, mie_sweep_result.csv);批量扫描时最容易忽略的是 n_max 随 x 变化这件事。上面代码里每个点都重新算 n_max虽然多花了一点时间但保证了每个点的截断都够。如果为了省事固定一个 n_max小粒径点会过算大粒径点会欠算扫出来的曲线在两端会有莫名其妙的翘尾。导出成 csv 之后不管是拿去 Origin 画图还是丢进 Python 做拟合都方便比在 MATLAB 里截图靠谱。4. 避坑与排查五个让结果对不上的常见问题4.1 消光系数出现负值现象跑出来的 Qext 是负数或者 Qabs 比 Qext 还大。原因几乎总是折射率虚部符号搞反了。物理教材里吸收性介质写成 m n - ik但这份程序的递推公式是按 m n ik 推的输入时虚部要取正值。解决方法是检查赋值那行把 m m_real - 1i*m_imag 改成加号。如果改完还是负的再看一眼波长和半径的单位是不是都用了米有人半径填了微米没乘 1e-6x 大了三个量级级数截断完全不够结果也会乱。4.2 大粒径时 Qext 收敛到 2 但曲线有毛刺现象粒径参数超过 50 之后Qext 理论上应该围绕 2 振荡并逐渐趋近但扫出来的曲线毛刺很多。原因是 n_max 余量不够共振区的 Mie 系数衰减慢截断误差表现为高频振荡。解决办法是把余量系数从 1.2 提到 1.5代价是计算时间增加但曲线会平滑很多。另一个可能是向下递推的种子值 1e-30 在超大 x 时不够小可以改成 1e-50不过这种情况很少见一般 1.5 倍余量就够了。4.3 不对称因子 g 超出 [-1, 1] 范围现象g 算出来是 1.2 或者 -1.5 这种离谱值。原因通常是 Qsca 接近零时除法的数值不稳定比如在某个特定粒径下散射系数恰好很小g 的公式里除以 x^2*Qsca 就会放大误差。解决办法是在计算 g 之前加一个判断Qsca 小于 1e-10 时直接把 g 置零并给个警告。另一种可能是 a_n、b_n 的交叉项求和时索引越界检查循环是不是写成了 n 1:n_max 而不是 n 1:n_max-1因为交叉项要用到 n1。4.4 Octave 下中文注释乱码导致语法错误现象在 Octave 里运行报 parse error指向一行看起来没问题的注释。原因是文件编码和 Octave 的默认编码不一致中文注释被解析成了非法字符。解决办法是用编辑器把文件另存为 UTF-8 无 BOM或者在 Octave 里用 warning(off, Octave:language-extension) 绕过但更彻底的做法是把注释改成英文毕竟这份代码的注释量不大改起来不费事。4.5 批量扫描时内存暴涨现象扫几千个点之后 MATLAB 卡死或者报 out of memory。原因是每次循环都在函数里新建 psi、chi、D 这些数组虽然单个不大但循环次数多了垃圾回收跟不上。解决办法是在循环外预分配所有输出数组并且在 mie_coefficients 里把临时数组用 clear 清掉或者干脆把函数改成接受向量化输入。更简单的做法是分段扫描每 500 个点存一次结果清一次工作区变量。5. 二次开发改造把单点计算改成可复用的光学参数模块拿到这份程序之后大多数人不会止步于跑单点。真正有价值的改造方向是把它封装成一个可复用的函数模块输入可以是向量化的粒径和波长输出直接是一个结构体这样在上层做反演或者拟合时不用每次手写循环。我一般会新建一个 mie_solver.m把参数校验、n_max 自动计算、批量循环全包进去对外只暴露一个接口。function result mie_solver(lambda, r, m) % 向量化米散射求解器 % lambda: 波长标量或向量单位米 % r: 半径标量或向量单位米 % m: 复折射率标量 % 返回结构体字段为 Qsca、Qext、Qabs、Qback、g尺寸与输入广播后一致 [Lambda, R] ndgrid(lambda, r); x_vec 2 * pi * R(:) / Lambda(:); n_pts numel(x_vec); Qsca zeros(n_pts,1); Qext Qsca; Qabs Qsca; Qback Qsca; g Qsca; for k 1:n_pts x x_vec(k); n_max ceil(1.5 * (x 4*x^(1/3) 2)); [a_n, b_n] mie_coefficients(m, x, n_max); [Qsca(k), Qext(k), Qabs(k), Qback(k), g(k)] ... mie_efficiencies(a_n, b_n, x); end sz size(Lambda); result.Qsca reshape(Qsca, sz); result.Qext reshape(Qext, sz); result.Qabs reshape(Qabs, sz); result.Qback reshape(Qback, sz); result.g reshape(g, sz); end这个封装的关键在于 ndgrid 的广播和 reshape 的回填这样输入标量和向量混搭时输出尺寸自动对齐上层调用不用关心维度。余量系数我提到了 1.5因为封装之后通常用于批量反演宁可多算一点也不能在共振区丢精度。验证方法很简单拿改造前的单点结果和改造后的对应元素做减法差值应该在 1e-12 量级如果对不上八成是 reshape 的维度顺序搞反了ndgrid 和 meshgrid 的输出顺序不一样这点我踩过不止一次。另一个值得做的改造是加一层缓存。做反演时同一个折射率会被反复调用几千次每次重算 Mie 系数纯属浪费。可以用 containers.Map 以 x 和 m 的组合做键缓存 a_n、b_n命中率在参数扫描场景下能到百分之九十以上。缓存键的构造要注意浮点数的精度问题我一般用 sprintf(%.6f_%.6f_%.6f, real(m), imag(m), x) 做字符串键精度截到六位小数既避免浮点误差导致的缓存穿透又不会因为精度太高而命中率暴跌。从那以后我每次改完散射代码都强制拿 x1、m1.50.01i 这组基准点跑一遍对比Qext 对不上 2.3 附近就不往下走。这个习惯帮我省掉了好几次把符号错误带进论文图表的麻烦。希望帮到你。本文还有配套的精品资源点击获取