做颗粒光散射计算的人大概率逃不开一个问题我要算的颗粒是任意形状不是标准球Mie理论给不了答案应该怎么办如果你也卡在这那我强烈建议你亲手写一遍离散偶极子近似DDA的MATLAB实现。我记得第一次跑通自己的DDA代码把球形颗粒的消光效率和Mie理论叠在一张图上对比时那种“这玩意儿真的能算出来”的兴奋感比直接用现成工具强烈得多。这套方法的核心思路非常朴素把任意形状的颗粒切成一个个小体积单元每个单元等效成一个偶极子然后自洽求解所有偶极子在入射光照射下的响应最后再汇总成全空间的散射场。它对形状没有要求球、椭球、立方体、团簇、甚至随机粗糙表面颗粒都能算这也是它在纳米光学、大气辐射、生物医学光散射领域里被广泛使用的原因。这篇整理会从物理原理一路讲到MATLAB工程落地包括网格生成、极化率赋值、耦合方程构建、迭代求解、截面和远场角分布的后处理最后再补上我调试过程中踩过的一堆坑。适合刚接触DDA、想自己从零搭一套可运行代码的人参考也适合已经会用现成开源包、但想知道“里面到底干了什么”的朋友。1. DDA在解一个什么样的“散射问题”1.1 为什么不用Mie理论形状自由是核心卖点Mie理论是光散射分析的基石但它只对均匀球体有解析解。一旦颗粒变成椭球、立方体、圆柱、甚至两个球粘在一起的二聚体Mie理论就无能为力了。虽然也有针对椭球的扩展解析方法但适用范围依然极其有限更别说随机取向、表面粗糙、多层结构这些更复杂的场景。DDA的思路是“以离散换形状自由”。它把颗粒占用的空间离散成规则格点每个格点是一个小立方体体积是d³。只要d远小于波长颗粒内部和外部的电磁场在这个尺度上变化足够缓慢那么每个小立方体对光场的响应就可以用一个点偶极子来等效描述。这样再复杂的形状在程序代码里就只体现为一个“索引数组”——坐标为格点、在颗粒内部则标记为1在外部则标记为0。只要这个掩模数组能生成出来任意形状都能参与计算。相比有限时域差分FDTD和有限元FEM这类通用电磁场求解器DDA的离散方式非常“节约”它只在颗粒体积内布点颗粒外部不需要额外网格。颗粒在空气中时这个优势特别明显不浪费任何一个格点在入射波传播路径上。1.2 偶极子离散与自洽方程设入射场是平面波在颗粒内部第i个格点位置r_i处的入射场为E_inc(r_i) E0 exp(i k · r_i)颗粒的折射率是m_p周围介质折射率是m_med通常用相对折射率m m_p / m_med。每个偶极子的极化率为α_i偶极矩为P_i。问题是P_i不仅受入射场驱动还受颗粒内部其它所有偶极子辐射场的驱动。所以每一个偶极矩都要满足一个耦合方程P_i α_i [ E_inc,i E_scat,i ]其中E_scat,i是其他所有偶极子在第i个偶极子位置产生的散射场之和。这是典型的“多体自洽”问题本质上是一个3N维的线性方程组N是偶极子总数。整理成矩阵形式就是(1/α_1, 0, 0, ...0, 1/α_2, 0, ...... ) P A P E_inc这里的A矩阵包含了偶极子两两之间的相互作用对角线块是极化率的倒数非对角线块是点偶极子辐射场的格林张量。有个特别容易困惑的地方为什么方程里ij那一项用的是1/α而不是格林函数本体因为点偶极子在自己所在位置的辐射场是发散的物理上偶极子自身的“退极化场”已经包含在极化率α的定义里了。换句话说α是考虑了自场修正之后的有效极化率。所以在构建方程时对角线不能直接套用非对角线的格林函数公式而必须用1/α来表征每个偶极子的自响应。很多初写代码的人在这一步翻车把对角线也填成某个格林函数表达式算出来结果完全失真。2. 从零搭建MATLAB求解器五个关键模块2.1 网格生成从三维坐标到颗粒掩模第一步是在三维空间中生成规则格点。我习惯用ndgrid而不是meshgrid。因为meshgrid在三维下会把第一维和第二维的顺序搞错导致坐标和索引对不上。ndgrid生成的X、Y、Z数组维度是[ nx, ny, nz ]和空间坐标的直觉一致。生成一个球形颗粒的掩模代码如下% 物理参数长度单位统一用微米 lambda 0.5; % 入射光波长 d 0.010; % 偶极子间距单位微米即10 nm R 0.050; % 颗粒半径单位微米即50 nm % 生成覆盖颗粒的最小矩形网格 half ceil(R / d); xx (-half : half) * d; [X, Y, Z] ndgrid(xx, xx, xx); % 球体掩模 mask (sqrt(X.^2 Y.^2 Z.^2) R); % 提取颗粒内部格点坐标 pos [X(mask), Y(mask), Z(mask)]; N size(pos, 1);这里有几个细节值得注意。一是边界条件用“形心是否落在球内”来判断这简单但会让阶梯边界产生约半个格点的尺寸误差。如果颗粒很小、格点数很少这个误差会直接影响截面结果的精度。更精细的做法是对球边界做体积加权修正但对于刚上手的第一版代码直接用形心判断即可。二是格点数量N和颗粒体积的关系颗粒真实体积是4/3πR³离散后的体积是N*d³两者之间会有偏差。在调试收敛性时我会输出这两个体积做对比如果偏差超过百分之几说明网格太粗。2.2 极化率把材料参数映射到每个偶极子每个偶极子的极化率α由材料折射率和格点体积一起决定。最简单的方案是Clausius-Mossotti关系这也是DDA里最基础的极化率模型m 1.5 0.02i; % 相对折射率 V d^3; % 格点体积 epsilon m^2; alpha_CM 3 * V * (epsilon - 1) / (epsilon 2); alpha alpha_CM * ones(N, 1);这个公式的物理含义是把一个介质小球放在外场中其内部场会因为退极化效应而修正极化率自然和介电常数、体积强相关。CM关系形式简单对折射率不太高、格点相对波长不太大的场景已经足够。但如果你希望和成熟开源包的精度对齐CM关系在高折射率或较大格点时误差会明显偏大因为真实格点上的偶极子受到的晶格修正不只是简单的退极化场。业界常用的做法是Lattice Dispersion RelationLDR修正它会根据入射波的波矢、格点间距和材料折射率对极化率做进一步校准。LDR实现的细节涉及多项展开系数首次动手不建议自己推导我的建议是先用CM版本跑通全流程确认DDA整体架构没问题之后再去查阅文献或参考成熟开源包源码把极化率替换成LDR版本。这样分两步走能最大程度减少把“物理问题”和“工程bug”混在一起排查的痛苦。2.3 耦合矩阵的显式构建与内存边界对于第一个可用版本我建议直接显式构建整个耦合矩阵M然后用MATLAB的“\”直接求解或者交给bicgstab迭代。这样逻辑最清晰也容易调试。相互作用格林张量的构建是整段代码的核心。在MATLAB里先把偶极子坐标pos拆成三者分量然后双重循环遍历所有偶极子对。对于第i个偶极子和第j个偶极子点偶极子的散射场作用矩阵有一个标准闭式公式。我当时手写的时候喜欢把公式拆成“径向点积部分”和“垂直部分”来避免出错。M zeros(3*N, 3*N); E0vec [1, 0, 0]; kdir [0, 0, 1]; % 入射方向沿 z k 2 * pi / lambda; for i 1:N row0 (i-1)*3; % 对角线块自作用 M(row01:row03, row01:row03) eye(3) / alpha(i); for j 1:N if i j continue; end rvec pos(i,:) - pos(j,:); r norm(rvec); rhat rvec / r; g 1/r^2 - 1i*k/r; % 点偶极子格林张量 Aij exp(1i*k*r)/r * ( ... k^2 * (rhat.*rhat - eye(3)) ... (3*rhat.*rhat - eye(3)) * g ); % 方程左端 M * P E_inc因此非对角块是 -Aij col0 (j-1)*3; M(row01:row03, col01:col03) -Aij; end end这段代码可以跑但循环是O(N²)N在几百时没有问题N到两三千时已经很吃力到五千以上基本没法等。再看内存一个3N乘3N的复数稠密矩阵每个元素占16字节N1000时是144MBN2000时接近580MBN5000时直接超过3.6GB。所以显式矩阵方案的上限大约在N1000到1500左右超过这个量级必须换成“矩阵向量算子迭代求解”的进阶模式。这部分我在第五章展开。需要特别说明的是坐标和索引的分量排布约定。我习惯把未知向量P打包成[px(1:N); py(1:N); pz(1:N)]的结构这样第i个偶极子的三个分量分别在第i、Ni、2Ni的位置。矩阵构建时二维块的行列索引要严格对应这个排布。另一种常见排布是[px1, py1, pz1, px2, py2, pz2, ...]的交错式两种都可以但后处理提取分量时容易把自己绕晕。建议从一而终。2.4 求解器直接解还是迭代解小N场景直接用M\b简单可靠。当N变大、矩阵无法显式存储时必须走迭代路线。两个推荐gmres和bicgstab。MATLAB内置的gmres只需要提供矩阵向量乘的函数句柄不需要转置矩阵向量乘这对DDA特别友好因为格林张量矩阵不是对称矩阵手写转置算子容易出错。bicgstab收敛通常更快但它同时需要A函数和A函数第一版代码不建议上来就搞这个。所以我的顺序是N小于1500用M\bN更大用gmres配一个Afun函数句柄。迭代求解的一大痛点是收敛判据。DDA矩阵的条件数会随着颗粒尺寸、折射率虚部的变化剧烈波动。金属颗粒在等离子体共振波长附近极化率很大矩阵条件数明显恶化迭代很容易“看起来收敛了结果却是错的”。我自己遇到过两次这种伪收敛后来都会额外检查残差对应的物理量比如用求出的P向量重算一次C_ext看它和迭代中间值是否稳定而不是单纯相信求解器返回的flag。2.5 后处理截面与远场方向图解出P向量后最常见的物理量是消光截面、吸收截面和远场角分布。消光截面C_ext可以用光学定理计算% 入射场各偶极子处的值 Einc zeros(N, 3); Einc(:,1) exp(1i * k * (pos * kdir.)); % 这里的kdir是行向量...细节略 for i 1:N Einc(i,:) E0vec * exp(1i * k * dot(pos(i,:), kdir)); end % 偶极矩P向量拆分 px Pvec(1:N); py Pvec(N1:2*N); pz Pvec(2*N1:3*N); % C_ext 4*pi*k/|E0|^2 * sum( imag( P_i . conj(Einc_i) ) ) Cext 4 * pi * k * sum(imag(px .* conj(Einc(:,1)) ... py .* conj(Einc(:,2)) ... pz .* conj(Einc(:,3))));注意这里已经假设入射场振幅|E0|1并采用了空间相位因子exp(i k·r)的约定。如果你换用了其他约定公式中的符号会变化这是跨程序对比时最容易出问题的地方。远场散射方向图则是把每个偶极子作为辐射源对方向r̂_s做相参叠加E_s(r̂_s) ∝ k² * exp(ikr)/r * Σ_i exp(-i k r̂_s · r_i) * [ r̂_s × (r̂_s × P_i) ]要对一组给定的方位角(θ, φ)计算角分辨散射强度只需要再写一个循环遍历方向采样点。散射截面C_sca等于远场强度对全立体角积分再除以入射强度。吸收截面则用C_abs C_ext - C_sca得到这样比直接用近似公式更稳健也避开了“吸收公式里要不要含偶极子相干项”的坑。我第一次写吸收截面时就踩过这个坑直接对每个偶极子单独算辐射功率再求和忽略了偶极子之间的相干项算出来的C_sca和远场精确积分差了一大截。从那以后我后处理一律走远场积分的路子。3. 上手跑通示例50nm介质球3.1 参数选择与网格密度估算一个稳妥的起步算例是半径R50nm、相对折射率m1.50.02i的介质球入射波长λ500nm。这样颗粒尺寸远小于波长散射处于瑞利-米氏过渡区Mie理论可以给出精确参考值便于验证。网格间距d的选择经验是|m| k d 1。这里k 2π/λd 10nm时代入得到m约1.5、λ500nm|m|kd约0.19满足条件。N约为419。这个规模用显式矩阵求解在普通笔记本上几秒钟就能完成非常适合当第一个测试用例。如果你希望更高精度可以把d减小到5nm。但那样N会涨到接近4200显式矩阵内存达到约2.5GB普通机器已经很难受必须切到算子模式和迭代求解。一种折中的中间档是d6nm左右N约1200到1500显式矩阵约250MB还能一战。这个趋势说明一个现实DDA的精度是靠网格密度换来的而网格密度直接决定内存和计算量。所以好的DDA实现必须在“显式便利”和“算子高效”之间选择合适路径。3.2 与Mie理论对比验证跑完代码下一步一定是拿着结果和Mie理论做对比。这一步建议不要只对比一个波长而是扫一个波长范围把C_ext随λ变化的曲线叠在一起看。我的对比流程是这样的固定颗粒半径和折射率模型在λ400nm到900nm之间取20到30个点每个点都调用DDA计算C_ext同时调用一个Mie理论函数算出C_ext_ref。然后把两条曲线画在一起计算相对误差。你会发现大部分区域的误差都在可接受范围内但在某些特定波长可能会出现系统性偏高或偏低。出现偏差先不要急着怀疑代码优先检查三件事网格d是否满足|m|kd1的条件颗粒边界阶梯效应对体积的误差是否超过2%折射率虚部符号和相位约定是否和Mie参考程序一致。我印象中至少有一次查了一整天最后发现是Mie参考函数使用的折射率虚部是负号约定而DDA代码是正号约定导致吸收项符号反了误差大得离谱。3.3 从球体扩展到任意形状DDA真正的价值在非球形。把网格生成部分的mask表达式换一下就能立刻从球变成其他形状。比如椭球% 椭球半轴 a, b, c mask (X.^2/a^2 Y.^2/b^2 Z.^2/c^2 1);立方体mask (abs(X) Lx/2) (abs(Y) Ly/2) (abs(Z) Lz/2);两个球组成的二聚体可以先用两个mask判断再做逻辑或。我做过一个随机取向椭球的代码只是把pos重新随机旋转一下再生成掩模整个DDA核心逻辑一行没改。这就是DDA对初学者的最大友好之处形状只决定掩模掩模在物理计算里只表现为格点列表和数量。你写好了球体就写好了所有形状。4. 常见问题与避坑实录4.1 误差偏大到底该调哪个参数这是后台被问得最多的问题。误差来源可以粗暴分成两类空间离散误差和边界逼近误差。空间离散误差由d过大引起。颗粒内部场的空间变化由波长和折射率决定折射率越高电磁场在颗粒内部变化越剧烈需要更小的d。经验判据|m|kd1是底线实际工程我一般取d λ / (15 |m|)这样C_ext误差通常能控制在百分之几以内。边界逼近误差则由球形边界上格点“阶梯化”引起。对同一个d颗粒越小、边界附近的格点占比越高相对误差越大。比如50nm半径的球用10nm格点直径方向只有10层格点边界格点比例很高误差明显如果换成200nm半径的球同样的d边界占比降低误差就会小很多。所以小颗粒想要同样精度必须付出更细网格的代价。我的习惯是做一个网格收敛测试把d从20nm逐步减到8nm观察C_ext是否单调逼近某个稳定值。如果结果还在明显漂移那就继续加密。做这个测试时不要忘了记录每个d对应的N和内存占用因为你会很快撞到显式矩阵内存上限。4.2 迭代求解器卡死或“伪收敛”用gmres或bicgstab时最常见的问题是迭代几百步不收敛或者flag为0但物理结果明显错误。先说卡死。高折射率金属颗粒、或者颗粒尺寸较大时矩阵对角占优性下降迭代收敛会很慢。一个有效的预处理是对角预处理把整个矩阵按对角块做归一化。对角块恰好是各偶极子的1/α物理上相当于把偶极子的自响应剥离剩余的相互作用是短程辐射项条件数会改善很多。MATLAB中可以用gmres(Afun, b, restart, tol, maxit, L, U)传入预处理矩阵L和U。简单的Jacobi预处理就够用。再说伪收敛。gmres返回的残差是基于矩阵方程M P E_inc的但物理上我们更关心截面量。有时候残差降到1e-8C_ext还在随迭代步数缓慢漂移这是因为C_ext是偶极矩和入射场的点积某些高频误差分量可能在残差范数里权重很小却在物理量里权重很大。我处理的方法是每隔一定迭代步数把当前的P向量取出来算一次C_ext画一条C_ext随迭代步数的曲线看它是否平台。平台震荡剧烈说明迭代参数没调好平台单向漂移说明收敛判据太低。这个习惯后来救了我好几次。4.3 内存不足的工程化解法显式矩阵方案在N到1500以上就开始吃力2万甚至10万偶极子的大规模运算必须彻底放弃A矩阵。这里给出的工程路径是一是在不改变物理模型的前提下把数据类精度从double切到single。但DDA的复平面运算里精度损失可能让消光截面误差放大所以只对极度内存受限的场景用且必须和double版本做对比验证。二是把“显式矩阵”换成“矩阵向量算子”。写一个函数输入P向量输出P在方程左端的像。算子内部不需要构建完整矩阵而是按偶极子对关系循环累加。N在几千时MATLAB的向量化循环还能扛住。三是更进一步利用DDA相互作用矩阵的平移不变性用FFT加速矩阵向量乘。这一步是大型DDA性能的关键也是成熟开源包的核心做法。4.4 单位、符号约定与跨程序对比单位问题几乎是每个写DDA的人都会踩的坑。我的建议是代码内统一用微米作为长度单位所有涉及频率、波矢、截面的量都据此推导。不要混用纳米和微米否则一个10^-3的系数错误会让你怀疑整个物理模型。符号约定上DDA有两种主流选择空间因子exp(-i k·r)和exp(i k·r)。它们会显著改变插值、格林函数和吸收项的符号。我使用exp(i k·r)约定对应的折射率虚部为正代表材料吸收。在和外部Mie程序对比前先确认双方约定一致。一个小技巧是算一个已知解析解的小球在真空中的消光效率虚部设为零且折射率近1看结果是否趋近于零。这样就可以快速锁定符号方向。还有折射率应该是相对折射率还是绝对折射率。颗粒在空气里两者相等但颗粒在水里就完全不同。DDA里的极化率和格林函数都假设背景介质是均匀且无源的所以入射场、波矢、极化率都必须在背景介质折射率下重新定义。我接触过几个初学者把颗粒绝对折射率直接当成相对折射率用导致共振波长整体偏移排查时非常迷惑。5. 性能优化进阶不做FFT就慢到没朋友5.1 矩阵向量算子的实现思路当N超过两三千显式矩阵开始不可行第一步还是实现Afun算子即不存储A矩阵而是在函数内动态计算A*x。这样内存占用从O(N²)降到O(N)时间则是每次矩阵向量乘都是O(N²)循环。Afun的基本结构是输入x向量先拆成三个分量px、py、pz然后对每个偶极子i循环所有j≠i累加格林张量作用在P_j上的结果。虽然还是O(N²)但因为没有保存3N乘3N的复数矩阵N到5000时内存不再是瓶颈只是时间涨得厉害。MATLAB里写的第一个Afun很容易写成三重循环外层i、内层j、每层里面再算3乘3矩阵乘。这样在N3000时一次矩阵向量乘可能就要好几分钟迭代几百步根本无法忍受。必须改成向量化对每个i一次性用数组运算算出所有j的位移和格林项。核心是避免在j循环里使用norm函数改为对整个j向量作减法、求模、广播点积。我建议第一次写Afun时先用显式矩阵的结果做验证随机生成P分别用A*P和Afun(P)的结果对比最大误差应小于1e-10。通过这一步再进入迭代求解能省掉后面好几个小时的排查。5.2 平移不变结构与FFT加速DDA的相互作用矩阵M只依赖于两个格点之间的相对位移Δr r_i - r_j。这意味着如果所有偶极子都排在一个完整、规则的矩形网格上那么格林张量矩阵就具有Toeplitz结构矩阵向量乘可以用三维FFT加速到O(N log N)。注意这里有个容易忽略的细节颗粒形状通常是不规则的偶极子只占据矩形网格的一部分。FFT加速时不能在稀疏的子集上直接做卷积而必须先把偶极子响应扩展到完整矩形网格上非颗粒区域的偶极矩设为0再做三维卷积最后从结果中提取颗粒区域内格点的部分。这个过程相当于把空外部区域的贡献显式算成了零物理上完全一致但FFT必须做一次三维尺寸的卷积内存占用取决于包围盒尺寸而不是颗粒体积本身。实际实现FFT版Afun时要处理的三个格林张量分量分别对三个偶极矩分量做卷积输出三个散射场分量总共是9次三维卷积。我建议分阶段完成先实现一个朴素FFT版本在颗粒近似填满矩形网格的场景下验证再用细长颗粒验证稀疏场景下的加速度和内存代价。FFT加速后N到数万甚至数十万都成为可能。颗粒出现凹槽或尖锐结构时包围盒体积可能远大于颗粒体积FFT的内存优势会被削弱这也是新型DDA变体试图解决的方向之一但对大多数应用场景标准FFT加速已经足够用。5.3 给新手的优化优先级如果你想把代码从“能跑”升级到“能算大颗粒”我建议按照以下优先级推进不要一上来就攻坚FFT。第一步把显式矩阵改成Afun算子配合gmres跑通N3000到5000的算力边界。这一步能帮助你理解迭代求解和内存瓶颈。第二步优化Afun的向量化把O(N²)循环做到尽可能快。这里通常能做到一个量级的加速并且Afun本身可以作为后续FFT版本的验证基准。第三步实现FFT矩阵向量乘。这也是工作量最大的一步建议把格林张量公式拆开写清楚分块测试卷积结果不要一步到位。第四步加LDR极化率修正。当你已经确信DDA算法本身没问题再用LDR去抠那百分之几的精度边际收益最大。优化要随时回到物理验证每改一个模块都要用上一章的50nm球在某个波长下重算一遍和Mie理论对比。不然你可能在加速的同时悄悄引入了一个符号错误直到算一个大颗粒时才暴露那才是最难受的调试场景。最后再分享一个我个人的小习惯每次改完网格生成逻辑我都顺手打印一行体积信息比较离散体积N*d³和理论体积。这个检查几乎不花时间却能提前预警很多边界上的隐蔽问题。DDA这个方向看起来公式多、矩阵大但真正把它跑通之后你会明显感觉到所谓“任意形状散射计算”并没有想象中那么高不可攀。如果你也被某个形状奇怪的颗粒困扰着强烈建议照着上面的代码框架从最小算例开始一步步把它算出来。