基于ABCD分块矩阵的Mie散射MATLAB实现：稳定求解大粒径粒子-北京尧图网络

简介一份面向光学与大气科学领域研究者的Matlab资源包围绕Mie散射理论提供计算系数an、bn的m脚本并以txt说明梳理理论要点与参数设置。压缩包共两个文件m文件为可直接调用的函数实现txt文件对理论背景、散射系数定义及代码调用方式作了注释整体仅1KB便于对照学习。该资源已有五百一十七人学习适合正在学习光散射理论、需要快速上手Matlab计算的本科生或科研人员。通过阅读脚本与说明可了解如何根据波长、粒径、折射率计算an和bn进而分析前向和后向散射强度与极化特征也能在此基础上扩展自定义程序用于气溶胶、雾滴等微粒光学性质研究。包内未直接给出cn和dn代码但txt提及可通过求导获得为理解相位结构函数留下线索。1. 先把标题拆开mieabcd 到底在解决 Mie 散射里的哪个痛点拿到这个标题我第一反应是把它拆成 Mie ABCD MATLAB 三块Mie 是球形粒子散射的经典理论ABCD 是把边界条件方程组组织成 2×2 分块矩阵的写法落点是 MATLAB 实现。很多做纳米光学、大气辐射传输、生物组织散射的同行估计都在某个时刻需要这样一段代码——输入粒径、波长、折射率输出消光截面、散射截面和角分布。常见教科书递归写法的痛点是粒径稍大就数值溢出而用 ABCD 分块矩阵组织递推能把粒径上限撑大几个量级多层核壳粒子也好扩展。这篇文章顺着标题把理论、代码和验证串一遍不需要额外下载专用工具箱MATLAB 基础环境就够。2. 从边界条件到 ABCD 矩阵为什么这套组织方式能避开大粒径溢出2.1 教科书递推公式里藏着的两个不稳定因子Mie 散射求解的硬核是计算散射系数 an 和 bn。入射平面波展开成矢量球谐函数后在球形粒子表面施加电场、磁场切向连续条件整理出来的 an、bn 是几类球贝塞尔函数和汉克尔函数的比值。我早年按教材一步步写代码时第一版是这样处理的直接用besselj和besselh逐阶算 psi_n(x) 和 xi_n(x)然后代入教材公式。算直径 100nm 的粒子完全没问题但把直径拉到几个微米、尺寸参数 x 超过 50 之后结果开始出现 NaN。原因很直白。Riccati-Bessel 函数 psi_n(x) x·j_n(x) 在阶数 n 接近或超过 x 时衰减极快而 xi_n(x) x·h_n^(1)(x) 同时包含衰减和振荡分量。分子里的 psi_n(x) 和分母里的 xi_n(x) 数值跨度可能相差几十个数量级直接相除就是灾难。更麻烦的是 D_n(mx) 这个对数导数项它在 mx 的实部很大时会出现指数级增长向上递推完全不可控。所以在 Mie 计算里稳定性的第一原则是psi 和 xi 用向上递推D_n(mx) 必须用向下递推。所谓“向下递推”是从一个足够大的阶数 N 开始按 D_{n-1}(z) n/z - 1/(D_n(z) n/z) 往小阶数推这样能保证数值沿着“衰减方向”走误差不会被放大。这套思路是 Wiscombe 在 1980 年前后系统整理出来的后来几乎所有好用的 Mie 程序都继承了它。2.2 用 2×2 分块矩阵重写边界条件an 和 bn 从哪来把边界条件用矩阵形式写出来落在纸上是这样的入射场展开系数构成一个列向量散射场展开系数构成另一个列向量边界条件把它们联系起来。对每个阶数 n独立形成一个 2×2 的线性关系这就是 ABCD 分块矩阵最朴素的来源。具体到代码实现我只关心四组中间量A_mat 和 B_mat 是两种模式下的对数导数组合再分别与 psi_n、xi_n 以及 psi_{n-1}、xi_{n-1} 组合成 A_blk、B_blk、C_blk、D_blk 四个向量。最终 an A_blk / B_blkbn C_blk / D_blk。你可能要问既然最后落点只是 an、bn为什么非要用 ABCD 这个名字因为在多层球场景里每一层界面都能写出一组 ABCD 块外层的结果由内层矩阵连乘得到写成矩阵形式后扩展逻辑非常自然。此外ABCD 这套组织方式也方便和 T 矩阵方法互相印证——T 矩阵本身就是把入射场与散射场之间的线性映射完整写成一个矩阵ABCD 相当于在做 T 矩阵的逐阶分解。对只算单层球的人来说理解到“ABCD 就是一组中间变量”就够了对想做核壳结构的人这个视角后面能省不少事。2.3 三种实现路线对比什么时候用哪一种我见过很多人在 MATLAB 里实现 Mie 散射路线大致分三类。第一类是最直接的用besselj和besselh从 0 阶推到 nmax 阶再逐项套教材公式。优点是代码短、容易读懂缺点是 x 一大就翻车而且每次都要算出全部贝塞尔函数浪费在无关阶数上。第二类是 Wiscombe 风格的稳定递推也就是上面说的“psi 向上、D 向下”单层球场景下又快又稳绝大多数公开代码属于这一类。第三类是 T 矩阵或 ABCD 分块形式把递推过程显式组织成矩阵连乘单层球场景下计算量和第二类差不多但扩展多层球、非均匀球时优势明显。我个人的选择是项目里只算均匀球就直接用第二类稳定递推代码量最小只要预期后面可能做核壳、多层或粒径分布扫描就按第三类 ABCD 来写省得日后重构。标题里既然写了 mieabcd下面的实现就按第三类来组织同时保证它和 Wiscombe 稳定递推在数学上完全等价。3. 在 MATLAB 里把 mieabcd 跑起来最小实现与三组必备参数3.1 主函数 mieabcd_anbn参数、输出和截断准则代码第一件事是定义输入输出。输入就三个量尺寸参数 x 2πr/λ·n_medium相对复折射率 m n_particle / n_medium以及最大阶数 nmax。nmax 不传也行我用 Bohren Huffman 的经典截断准则自动估nmax ceil(x 4·x^(1/3) 5)这个经验式在绝大多数场景下够用x 到几百也不怕。输出是 an、bn 两个长度 nmax 的复数列向量以及实际使用的 nmax。我习惯把配套函数统一放在 miematlab_utils 目录下主程序叫 mieabcd_anbn.m后处理函数叫 mieabcd_ang.m。这样同事拿到整套代码一眼就能看出哪个文件是求解核心哪个是角分布后处理。3.2 核心代码从 D_n(mx) 向下递推到 ABCD 分块function [an, bn, nmax] mieabcd_anbn(x, m, nmax) % mieabcd_anbn 用 ABCD 分块矩阵组织 Mie 散射系数求解 % 输入: % x : 尺寸参数, x 2*pi*r/lambda*n_medium % m : 相对复折射率, 虚部大于0表示吸收 % nmax : 最大计算阶数, 缺省时自动截断 % 输出: % an, bn : Mie 散射系数, 长度 nmax 的复数列向量 if nargin 3 || isempty(nmax) nmax ceil(x 4 * x^(1/3) 5); end % 防御 x0 或 m0 的退化输入 if x 0 || abs(m) 0 an zeros(nmax, 1); bn zeros(nmax, 1); return; end n (1:nmax).; % ---- 1. 向下递推 D_n(mx) ---- mx m * x; nstart nmax 50; % 多推50阶作为“无穷远”起点 D zeros(nstart, 1); for jj nstart-1 : -1 : 1 D(jj) (jj1) / mx - 1 / (D(jj1) (jj1) / mx); end Dmx D(1:nmax); % D(1) 对应 D_1(mx) % ---- 2. 向上递推 Riccati-Bessel 函数 ---- psi zeros(nmax, 1); xi zeros(nmax, 1); psi_prev sin(x); % psi_0 psi_cur sin(x)/x - cos(x); % psi_1 psi(1) psi_cur; xi_prev -1i * exp(1i*x); % xi_0 xi_cur (1/x - 1i) * exp(1i*x); % xi_1 xi(1) xi_cur; for jj 2:nmax psi_next (2*jj-1)/x * psi_cur - psi_prev; xi_next (2*jj-1)/x * xi_cur - xi_prev; psi(jj) psi_next; xi(jj) xi_next; psi_prev psi_cur; xi_prev xi_cur; psi_cur psi_next; xi_cur xi_next; end % 构造 psi_{n-1} 和 xi_{n-1} 向量 psi_prev_vec [sin(x); psi(1:nmax-1)]; xi_prev_vec [-1i*exp(1i*x); xi(1:nmax-1)]; % ---- 3. ABCD 分块与 an/bn ---- A_mat Dmx ./ m n ./ x; % TE 模式对数导数组合 B_mat Dmx .* m n ./ x; % TM 模式对数导数组合 A_blk A_mat .* psi - psi_prev_vec; B_blk A_mat .* xi - xi_prev_vec; C_blk B_mat .* psi - psi_prev_vec; D_blk B_mat .* xi - xi_prev_vec; an A_blk ./ B_blk; bn C_blk ./ D_blk; end这段代码有三个地方值得单独说。第一向下递推的起点我加了 nstart nmax 50而不是直接从 nmax 开始。原因是 D_n(mx) 的递推需要一小段“预热距离”来消除起始近似值 0 带来的误差多推几十阶成本很低但能明显提升小粒径场景下 an、bn 的尾阶精度。第二psi_cur 和 xi_cur 的初值必须严格对应 n1 的 Riccati-Bessel 函数psi_1 sin(x)/x - cos(x)xi_1 (1/x - i)·exp(i·x)这是整个递推链的地基写错一位后面全错。第三ABCD 四个分块在代码里并不是真的二维矩阵而是四个并列向量因为每个阶数 n 的方程互相独立把它们组织成向量就是对角化后的分块矩阵这样做既省内存也方便后续做矩阵连乘扩展。3.3 远场角分布把 an/bn 变成可画图的结果散射系数算出来之后消光效率 Qext 和散射效率 Qsca 能立刻得到Qext (2/x²)·Σ(2n1)·Re(anbn)Qsca (2/x²)·Σ(2n1)·(|an|²|bn|²)吸收效率 Qabs Qext - Qsca。角分布需要再走一轮递推计算 π_n 和 τ_n 两个角函数。function [S1, S2, I1, I2] mieabcd_ang(an, bn, theta) % mieabcd_ang 由 an/bn 计算远场散射振幅与强度 % 输入: % an, bn : mieabcd_anbn 输出的散射系数 % theta : 0~pi 的列向量, 单位弧度 % 输出: % S1, S2 : 两个正交偏振方向的复散射振幅 % I1, I2 : 对应散射强度 nmax numel(an); pi_minus1 zeros(size(theta)); % pi_0 0 pi_cur ones(size(theta)); % pi_1 1 tau_cur cos(theta); % tau_1 cos(theta) S1 3/2 * (an(1) * pi_cur bn(1) * tau_cur); S2 3/2 * (an(1) * tau_cur bn(1) * pi_cur); for nn 2:nmax pi_next (2*nn-1)/(nn-1) * cos(theta) .* pi_cur ... - nn/(nn-1) * pi_minus1; tau_next nn * cos(theta) .* pi_next - (nn1) * pi_cur; coeff (2*nn1) / (nn*(nn1)); S1 S1 coeff * (an(nn) * pi_next bn(nn) * tau_next); S2 S2 coeff * (an(nn) * tau_next bn(nn) * pi_next); pi_minus1 pi_cur; pi_cur pi_next; end I1 abs(S1).^2; I2 abs(S2).^2; end这里 π_n 和 τ_n 的递推初值不要记错π_1 1τ_1 cosθπ_0 0。很多资料直接用 n 从 2 开始的递推式但初值对不上结果会在小角度区域出现莫名其妙的抖动。顺带一提θ0° 处的 S1 复振幅和消光截面存在解析关系这一步后面验证章节会用到所以代码里 theta 的采样最好把 0 纳入。4. 避坑指南扩散系数、截断阶数与复折射率的四类翻车现场4.1 现象一大尺寸参数直接 NaN问题出在 psi 向上递推把粒径拉到 2μm、波长 532nm、x 接近 12 时程序返回的 an、bn 出现 NaN。起初我以为是递推式写错了逐项排查后才发现是 psi_cur 和 xi_cur 在 x 较大时增长太快中间量超出了 double 的表示范围。原因在于 Riccati-Bessel 函数 psi_n 在 n 小于 x 时随 x 增长而 xi_n 还叠加了指数衰减振荡直接存数值很容易爆。解决方法是保持 psi 向上递推的同时把 xi_n 的指数因子提取出来xi_n exp(i·x)·P_n(x)P_n 是一个多项式型函数数值幅度温和得多或者干脆把整个递推改用对数域格式。我的代码里直接用 exp(1i·x) 构造 xi_0 和 xi_1本质上已经规避了大 x 溢出所以 x 到几百都不怕。如果你看到 NaN先确认 xi 的初值不是用 besselh 直接算的而是用了显式解析式。4.2 现象二Qext 比文献值小 5%多半是截断阶数不够别急着怀疑物理先看 nmax。我曾经算 x20 的粒子为了省时间把 nmax 手动设成 25Qsca 和文献值差了接近 3%。原因是高阶项虽然在幅度上衰减但叠加了 (2n1) 的权重低估高阶贡献对 Qext 影响不大对 Qsca 却很明显。解决方法是按 nmax ceil(x 4·x^(1/3) 5) 截断这个准则几乎不会少算。更稳的做法是把 nmax 翻倍再算一次 Qext两次相对变化小于 1e-8 才认为截断收敛这个习惯我后面专门写一节。4.3 现象三Qabs 出现负值复折射率虚部符号没统一MATLAB 里指数的时间因子约定、介质的增益还是吸收、代码里 m 虚部的正负这三者但凡有一个不一致就会出现奇怪的负吸收。标准约定是时间因子 exp(-iωt)吸收材料折射率 m n i·n虚部为正。代码中 xi 的构造用 exp(1i·x)对应汉克尔函数的出射波约定所以 m 的虚部必须按正值传进去。如果从某个材料库拿到的折射率虚部是负的先把它取正再喂给函数否则不只 Qabs 为负向后散射强度也会整体错掉。这是黑匣子最容易翻车的地方程序本身不报错但结果全盘皆错。4.4 现象四m1 时出现 0/0粒子“消失”反而算崩了还有个边界情况粒子折射率和环境完全一样散射理应严格为零但公式里 m1 直接让 A_mat 和 B_mat 变得一致分子分母同时趋近于 0浮点除法给出 NaN。解决方式是在函数入口加一条判断abs(m-1) 小于某个阈值时直接返回零向量。这个阈值我给 1e-6但注意它不能设太大否则近背景粒子的微弱散射会被误杀。类似的边界还有 x0 和 m0代码里都做了防御。5. 验证 mieabcd 算得对不对三个不依赖第三方库的独立检验5.1 用光学定理验证消光效率 QextMie 散射满足光学定理消光截面正比于前向散射复振幅的实部。对应到效率公式Qext (4/x²)·Re(S1(0°))。这个关系不依赖粒子折射率是散射理论本身的硬约束。我的验证脚本是这样写的% demo_mie_validation.m lambda 0.532e-6; r 0.15e-6; x 2*pi*r/lambda; m 1.33; [an, bn, nmax] mieabcd_anbn(x, m); % 效率 nvec (1:nmax).; Qext 2/x^2 * sum((2*nvec1) .* real(anbn)); Qsca 2/x^2 * sum((2*nvec1) .* (abs(an).^2abs(bn).^2)); Qabs Qext - Qsca; % 光学定理: 用前向散射振幅验证 Qext theta linspace(0, pi, 2001).; [S1, ~] mieabcd_ang(an, bn, theta); Qext_fwd 4/x^2 * real(S1(1)); % theta0 处的前向振幅 fprintf(x %.6f, nmax %d\n, x, nmax); fprintf(Qext %.12f\n, Qext); fprintf(Qext_fwd %.12f\n, Qext_fwd); fprintf(Qsca %.12f\n, Qsca); fprintf(Qabs %.12e\n, Qabs);如果 Qext 和 Qext_fwd 在小数点后 8 位以内对不上说明 an、bn 或 S1 至少有一个环节错了。我用这个测试抓过一次角分布递推的初值错误那种错误在效率上完全看不出只有光学定理能逮到。注意 theta 数组必须从 0 开始取前向振幅才有意义。5.2 用瑞利极限回归小粒径场景尺寸参数 x 远小于 1 时Mie 散射应退化为瑞利散射Qsca ≈ (8/3)·x^4·|(m²-1)/(m²2)|²。这是一个解析极限用来验证程序的低阶行为非常合适。x_small 1e-3; m_ray 1.5; [an_s, bn_s] mieabcd_anbn(x_small, m_ray, 1); % 强制只算 n1 alpha (m_ray^2 - 1) / (m_ray^2 2); Qsca_ray 8/3 * x_small^4 * abs(alpha)^2; nvec_s 1; Qsca_mie 2/x_small^2 * sum((2*nvec_s1) .* (abs(an_s).^2abs(bn_s).^2)); fprintf(Rayleigh Qsca %.12e\n, Qsca_ray); fprintf(Mie Qsca %.12e\n, Qsca_mie);这里有个细节nmax 必须强制设成 1因为高阶项在 x 很小时幅度是 x^(2n) 量级混进来会污染对比。跑出来两者相对误差一般小于 1e-4主要误差来自 x 不为零的高阶残留。我习惯把这段代码放在测试文件最前面每次改动主函数后先跑一遍瑞利极限挂掉意味着 psi 或 Dmx 的初值被碰坏了。5.3 用收敛性检查保护长期改动第三个验证是收敛性做法是把 nmax 翻倍重算 Qext 和 Qsca比较相对变化[an2, bn2, nmax2] mieabcd_anbn(x, m, 2*nmax); Qext2 2/x^2 * sum((2*(1:nmax2).1) .* real(an2bn2)); fprintf(相对变化 %.2e\n, abs(Qext2 - Qext) / abs(Qext));生产环境里我会把阈值设为 1e-8未收敛就把 nmax 往上加。这句检查看着笨实际比任何高级调试手段都管用——它把“截断阶数够不够”这个最隐蔽的变量固化成了自动检查。我现在的所有仿真代码里都保留了这段改完代码跑一遍能挡住 80% 的回归错误。6. 把 ABCD 矩阵变成你的工具箱多层球扩展与批量仿真小技巧6.1 多层核壳粒子ABCD 连乘的思路单层球只是 ABCD 的基层用法它真正的扩展点在多层球。每一层界面的边界条件贡献一个 2×2 传递矩阵整体散射系数由这些矩阵从内到外连乘得到。伪代码结构是初始化为单位阵从最内层开始对每一层计算 A_l、B_l、C_l、D_l 并左乘到累积矩阵上最后散射系数由累积矩阵右上块和右下块的比值给出。老读者可能发现这与 T 矩阵方法完全同构。如果你手里的代码只有 an、bn 递推没有显式 ABCD改多层球时要把内部场展开系数也保留下来然后逐层拼边界条件。这是可预期的工程量少则半天多则两天但换来的能力是核壳粒子、多层介质球甚至粒径分布包络的仿真都归到同一套框架里。6.2 批量仿真的两个工程习惯仿真里经常要对一组粒径或一组波长做扫描。朴素的 for 循环在 x 超过 100、nmax 超过 150 后会很吃力。我用两个习惯解决一是把独立工况丢给 parfor注意每个 worker 里都要能调用 miematlab_utils 路径下的函数二是准备一份“参考结果”存档跑通一组典型参数后把 Qext、Qsca、Qabs 存成文本以后每次改动程序后自动比对。这等于给代码买了后悔药改坏了能第一时间知道。个人教训早年间我给一组 20 个粒径做扫描跑完发现 Qabs 全是负的问题就出在折射率虚部符号没统一。那之后我把验证代码写成了函数每天动代码前先跑一次光学定理和瑞利极限等公司新同事也能跑通这套测试再放他们改主函数。希望你按这套路径走能少走那段我从 NaN 和负吸收里爬出来的弯路也希望这套 mieabcd 的实现思路帮到你。本文还有配套的精品资源点击获取

基于ABCD分块矩阵的Mie散射MATLAB实现：稳定求解大粒径粒子

相关推荐

从先验到后验：贝叶斯概率模型的工程实践

XAPP888 MMCM动态重配置：DRP接口实现FPGA时钟在线切换

AI Agent核心技术解析：从感知到决策的智能闭环

更多资讯

UUV分布式编队一致性控制仿真：应对通信时延、丢包与拓扑干扰

STM32与BMI088姿态融合解算：从四元数到欧拉角的嵌入式实现

软件测试自动化框架设计：接口与UI自动化分层落地的完整实践指南

Q-Learning样本复杂度解析与优化实践

移动应用与AI系统安全防护实战指南

GTC 2026前瞻：NVIDIA驱动与CUDA环境配置实战指南

今日更新

微信个人号API开发案例：图片/语音/视频消息怎么下载解析存到本地？

Python股票自动交易系统源码解析：从行情获取到风控的完整链路

视酷酷信IM源码部署指南：Java 17+MySQL 8.0.33实战配置

本周精选

防幻觉智能体架构：电源硬件设计中的大模型编排与工具化实践

OpenGL面试核心考点：渲染管线、矩阵变换与性能排查全解析

全国计算机二级C语言真题PDF高效刷题指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？