米散射MATLAB实现：从气溶胶到颗粒检测的五个光学参量计算

简介这份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(n1) bn .* conj(bn(n1)))) ... ((2*n 1) ./ (n .* (n 1))) .* real(an .* conj(bn)) ); end逻辑说明Qsca、Qext、Qabs是效率因子乘以几何截面 $\pi r^2$ 就是截面。Qback是后向散射效率注意那个 $(-1)^n$ 相位因子漏了它后向散射会算成前向。g是不对称因子范围 -1 到 10 表示各向同性接近 1 表示强前向散射。参数说明an、bn是 Mie 系数向量x是尺度参数。这段代码里g的求和项容易写错建议用已知结果比如 $x1$$m1.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趋近 2Qabs趋近 0如果 $k0$第三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的求和漏了(2n1)因子。解决确保两个求和用同一个nmax并且Qsca的权重是(2*n1)。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都强制走一遍合成数据验证不然一个符号错误就能让反演结果偏到姥姥家。希望帮到你。本文还有配套的精品资源点击获取