摘 要

MATLAB是一款历史悠久的专业数学软件，现在已经发展成了一种集数值运算、数据可视化、图形界面设计、仿真等功能于一体的综合性软件。

离散偶极子近似（Discrete Dipole Approximation DDA）是一种用来计算任意尺寸和形状的小粒子电磁散射场的方法。

本文利用了MATLAB和DDA方法，对在有入射光的影响下，小粒子周围的散射场进行了数值计算和仿真。第一部分根据电磁场理论，建立了DDA的模型，然后对散射场的计算公式进行了推导，并讨论了使DDA方法有效的判据和其准确性，此外还介绍了解线性方程组的方法：双共轭梯度法。第二部分开始编写MATLAB程序，详细地讲解了程序中用到的所有函数，及仿真需要的物理量；然后以图片的形式，从多个角度直观地显示出散射电场的场强大小，并作出了分析和总结。结果表明：DDA方法具有优良的性能，可以快速、准确地模拟出小粒子的散射场。

关键词：MATLAB；离散偶极近似；DDA；双共轭梯度法

Abstract

MATLAB is a professional mathematical software with a long history, which has developed into a comprehensive software integrating numerical operation, data visualization, graphical interface design, simulation and other functions.

Discrete Dipole Approximation is a method used to calculate the electromagnetic scattering field of small particles of arbitrary size and shape.

In this paper, MATLAB and DDA methods are used to numerically calculate and simulate the scattering field around small particles under the influence of incident light. In the first part, the model of DDA is established according to the theory of electromagnetic field, and then the calculation formula of scattering field is deduced, and the criterion of DDA method and its accuracy are discussed. The second part began to write MATLAB program, explained in detail all the functions used in the program, and the physical quantities needed for simulation; Then, in the form of pictures, the intensity of scattering electric field is displayed intuitively from many angles, and the results are analyzed and summarized. The results show that the DDA method has good performance and can simulate the scattering field of small particles quickly and accurately.

Key Words: MATLAB; Discrete Dipole Approximation; DDA; Biconjugate gradients method

目录

第1章 绪论 1

1.1 目的及意义 1

1.2 国内外研究现状 1

1.3 论文的主要内容 2

1.4 论文的章节安排 2

第2章 离散偶极子近似理论 4

2.1 DDA基本理论 4

2.2 DDA建模 4

2.3 散射场 5

2.4 DDA有效性的判据 6

2.4.1 判据一 6

2.4.2 判据二 7

2.4.3 判据三 7

2.5 双共轭梯度法 7

2.5.1 算法A 8

2.5.2 算法B 8

2.6 DDA方法的准确性 9

第3章 程序及仿真结果 11

3.1 MATLAB简介 11

3.2 程序说明 11

3.2.1 生成偶极子坐标 11

3.2.2 生成DDA系数矩阵 12

3.2.3 利用双共轭梯度法解线性方程组 13

3.2.4 计算场强 13

3.2.5 pcolor函数简介 14

3.3 仿真结果 14

3.3.1 参数设置 15

3.3.2 垂直于传播方向的散射截面 15

3.3.3 平行于传播方向的散射截面 17

第4章 总结与展望 19

4.1 工作总结 19

4.2 未来展望 20

参考文献 22

附录 23

致谢 27

第1章 绪论

1.1 目的及意义

离散偶极子近似（Discrete Dipole Approximation, DDA）是一种用来计算任意尺寸、形状的粒子的散射场的方法。

在电动力学中，计算物体对电磁波的散射和吸收实际上就是求解物体内部和周围的电磁场分布。麦克斯韦方程组从理论上可以解出介质存在时电场的分布。然而，由于其复杂性，解析解只能在对称系统中求得。因此，对于一般物体，通常采用数值近似法求出其周围的电磁场分布。离散偶极近似就是这样的方法。

后来证明，DDA 也可以从电场的积分方程中推导出来，通过将散射体离散成小的立方体。这个推导显然首先由Goedecke和O'Brien完成，然后由其他人进一步发展^[1]。需要注意的是，由DDA的两种推导得到的最终方程本质上是相同的。唯一的不同之处在于，基于积分方程的推导使我们对近似有了更多的数学理解，从而为改进方法指明了方向，而基于点偶极子的模型在物理上更加清晰。

小粒子的电磁散射在天体物理学、大气科学、海洋学等许多工程和科学领域有着重要的应用。由于该问题的重要性和复杂性,已经成为电磁散射理论的热门。

目前，航天遥感正朝着高空间和时间分辨率的方向飞快地发展着。然而，大气气溶胶和高空中的卷云会对遥感产生剧烈的影响。为了进一步发展遥感技术，大气气溶胶和卷云中的散射传输是亟待解决的关键问题。为了解决这一系列的重大科学问题，有必要对小粒子的电磁散射理论及其计算和仿真进行深入的研究探讨。

1.2 国内外研究现状

1964年，Howard DeVoe最早地建立了DDA的体系结构，并发展成了一种经典的物理模型。该模型可根据分子等单体的光学性质来求解分子晶体等聚集体的光学性质。DeVoe认为，如果找到聚合体的总偶极矩，就可以推导出它的消光值、吸收值和散射值。分子偶极矩包含入射电场和其他偶极电场，通过这些就可以建立并求解包含所有偶极矩的线性方程组。但是，在他的论文中，描述偶极子间电场的作用时使用的是静电场的表达式，而不是振荡偶极子的电场。

1973年，Purcell确立了DDA的基本原理和一些重要的概念，并求解了任意形状小粒子的光散射特性。他是在关于星际尘埃散射星光的研究过程中发展出这一方法的，因而并没有引用DeVoe的研究成果。相对于DeVoe建立的模型，Purcell明确了DDA的基本概念和在计算小粒子散射场方面的用途。Purcell第一次发表了如何利用三维空间下的偶极子模型仿真小粒子的散射，得到的结果比DeVoe的模型更加准确。Purcell通过Clausius-Mossotti关系式确定了偶极阵列中偶极子的极化率(polarizability) α{\displaystyle \alpha }与颗粒材料的介电常数ε{\displaystyle \epsilon }的关系（DeVoe是使用的分子消光系数导出的分子偶极极化率）^[2]。在此基础上，求出极化率就解决了DDA方法中的基本问题：偶极矩与局部电场之间的数值关系。这样就可以构建出描述偶极矩之间关系的线性方程组{\displaystyle {\boldsymbol {P}}_{i}=\alpha {\boldsymbol {E}}_{i}}。

1991年，Goodman指出，当偶极子组成的阵列在空间上具有周期性时，复数的共轭梯度算法中的矩阵乘法本质上是一种卷积运算，那么久可以利用快速傅里叶变换(Fast-Fourier Transform FFT)对原有算法进行改进，大大的减小了其时间复杂度和空间复杂度。

1993年，Draine与Goodma指出，使用Lattice Dispersion Relation (LDR)关系描述偶极子的极化率能够使DDA方法的解与米氏散射的结果更为接近，因而取代了之前常用的Clausius-Mossotti plus radiative reaction (CMRR)关系。

2008年，Maltsev和Yurkin从电磁场的散射理论出发，系统地研究了DDA方法中产生的误差与离散的子单元尺寸d的关系，推导出了DDA方法中离散化误差的表达式。并对正方体形状的偶极子和散射体实际形状之间的差别导致的误差进行了分析。

1.3 论文的主要内容

本文利用MATLAB计算了小粒子的散射场，对仿真的结果做出了分析。并讨论了偶极子数量对结果收敛性的影响，和迭代算法解线性方程组的计算效率。具体包括：

1. 计算偶极子周围的电磁场；
2. 推导小粒子的散射场公式；
3. 基于MATLAB编程，计算球形粒子的散射场；
4. 对程序中的主要部分进行讲解；
5. 根据显示的散射场图片，做出结论；
6. 讨论偶极子数目对结果收敛性的影响；
7. 讨论使DDA方法有效的判据；
8. 讨论DDA的准确性；
9. 总结离散偶极近似方法过去的成果，展望其未来。

1.4 论文的章节安排

本文共设4章，从实际应用出发，到理论推导和MATLAB的程序实现，系统地阐述了用离散偶极子近似方法模拟小粒子的散射场的过程。具体的章节内容安排如下：

1. 绪论：介绍了DDA方法的历史背景，说明了为何要采用这一方法；并简要总结了DDA近年的发展。最后对论文的总体结构进行了说明。
2. 离散偶极子近似理论：先对DDA的基本理论进行了概括，然后说明如何使用DDA方法建立粒子的模型，综合论述了散射场的计算方式、DDA有效性的判据及双共轭梯度法等。
3. 程序及仿真结果：根据第2章中的结论，利用MATLAB里相关的函数，编写程序，以图片的形式显示出散射场，并得出结论。
4. 总结与展望：与第1章相呼应，对DDA方法的成果进行总结，与其它方法进行对比，指出DDA在电磁散射仿真方法中的地位，并对未来该方法可能继续发展下去的研究方向进行了探讨。

第2章 离散偶极子近似理论

2.1 DDA基本理论

离散偶极子近似(DDA)，也被称为耦合偶极子法(Coupled Dipole Method, CMD)，是计算任意几何形状和组成的粒子对电磁波的散射和吸收的一般性方法。这种方法将小粒子离散成多个可极化单元，若单元的尺寸远小于入射波的波长，则可将其视为偶极子，从而使用足够多的偶极子阵列来模拟散射粒子^[3]。

在入射电磁波的诱导下，偶极子将发生极化。这些偶极子和入射场之间相互作用， 产生了一个线性方程组，求解可以得到偶极子的极化强度，所有需进行计算的散射量都可以从这些极化强度中得出^[4]。

DDA方法将在入射电磁波影响下小粒子的散射情况的求解转换成求解3N维线性方程组（其中N为偶极子数目），从而可以利用成熟的迭代算法，使整个运算过程具有良好的速度和精度。

2.2 DDA建模

本文的公式采用高斯单位制(Gaussian units)。

DDA方法建模的思路是利用尺寸为d的一系列立方体来离散化任意尺寸和形状的物体，每个立方格点上都有一个偶极子，给这些偶极子赋予极化率，那么它们便有了具体的物理意义。偶极子可以理解为一个极化的球体。

令偶极子的半径为a，因为偶极子的体积要等于立方单元的体积，所以：

于是可以推出：

可见偶极子的直径大于立方单元的边长，所以相邻的极化球体即偶极子之间具有重叠的部分。从表象上看，这体现了球形偶极子模型中的多极子效应，可以通过利用多种不同的模型来减小建模所带来的的误差。

所有立方单元的体积之和应等于被建模的物体的体积，所以：

其中，N为偶极子的总数，V为物体的体积。于是有：

其中，R_eq为物体的等效球体半径，令球体的体积等于物体的体积可得：

当物体为球形粒子时，等效球体半径R_eq就等于粒子的半径。

可见，在给定了物体的体积V（当物体为球体时给定半径）和偶极子总数N（考虑到结果的准确性，N应越大越好；但是另一方面，N越大对计算机的要求就越高，程序的时间复杂度和空间复杂度都将很大，计算效率大打折扣。N应满足的条件见2.4节）的条件下，立方单元的尺寸和偶极子半径都由上述公式确定^[5]。

2.3 散射场

假设中心位置为r_i的偶极子的极化率为α_i，P_i为该偶极子在局域场E_loc(r_i)处的极化强度，由电介质中极化强度和电场强度的关系得：

上式中E_loc(r_i)为r_i处入射场E_inc和其他所有偶极子在ri处场强的叠加：

E_inc为入射场，其振幅为E₀，波矢为k，是一个均匀平面波：

-A_ij·P_j为r_j处的偶极子在r_i处的场强，对-A_ij·P_j从j ＝ 1到N（且j ≠ i）求和，就得到了除r_i处偶极子外所有其他的偶极子在r_i处的场强的叠加^[6]：

A_ij·P_j的表达式如下：

上式中r_ij = | r_i - r_j |，为r_i处偶极子和r_j处偶极子间的距离。若补充定义A_ii = α_i^-1，那么全部偶极子的极化强度可以列成一个N维的线性方程组^[7]：

其中A为DDA方法的系数矩阵，A=Λ-T。Λ为极化率张量的逆矩阵，T为场极化张量，T可表示为：

上式中，为r_j指向r_i的单位矢量。式(2.11)是一个矢量方程组，因为是在3维空间中，每个矢量都有3个互相垂直的分量，所以式(2.11)本质上是一个3N维的线性方程组。求解这个线性方程组便能计算出全部偶极子的极化强度。再根据电场强度和极化强度的关系便可得到每一处的电场强度^[8]。

2.4 DDA有效性的判据

使用DDA方法需要将散射体离散为偶极子阵列，易知，偶极子数目越少所得到的结果就越不准确，但是较少的偶极子数目有利于计算机进行运算；相反，偶极子数目越多，得到的结果将越准确，不过这将对计算机硬件有较高要求。

可见，偶极子数目是一个很重要的参数，不能随意设置，应当有理论依据作为支撑点。这样才能兼顾准确性和计算效率，得到令人满意的结果。

在DDA的众多文献中，主要有以下几种DDA有效性的判定方法，我们可以根据实际情况从中选择一种合适的。

2.4.1 判据一

因为偶极子模型是离散的，由N个立方单元组成。显然，N的值越大，通过建模引起的误差就越小。但是，因为N只能取有限的值，不可能取无穷大，所以模型误差就一定会存在。

例如，对于球形粒子来说，表面应该是光滑的球面，但实际建立的模型表面形状会和粒子有一定的误差，即模型的表面具有一定的粒度。这是离散偶极近似方法中不能忽视的离散化误差。

为了使DDA方法能投入实际的使用当中，就需要使误差控制在一定的范围内，Draine根据理论推导和实际结果相结合，得出了使DDA方法有效的一个判据：若要结果的百分比误差小于，那么就要偶极子的总数大于某一个值^[9]：

其中，m为粒子的折射率，m可以为复数（吸收材料）。上式表明，当粒子的折射率m越大时，偶极子总数也应该越大，否则将会导致较大的误差。

2.4.2 判据二

从另一方面来看，要DDA方法得到的结果越准确，就需要每个偶极子内部的电场比较稳定，不能有太大的波动。即立方单元的尺寸d要远小于入射电磁波的波长，本文中使用波矢k来描述电磁波的性质。所以有如下判定公式^[10]：

上式与判据一类似，都说明了当粒子的折射率越大时，若要结果保持在一定的误差内，就需要更大的偶极子总数。