实验测量证实：生物组织、大地、岩石、水、液晶、聚合物材料等众多媒质的电参数展现了频率相关性(这些媒质称作电色散媒质)；不同媒质还可能展现不同的电色散特性^{[1, 2]}。至今已提出好几种经验模型以描述不同的电色散特性，一般可概括为两大类型^[3]：第一大类主要有Debye模型^[4](如人体组织、等离子体、土壤)、Drude模型^[5](如等离子体、金属)和Lorentz模型^[2](如光学材料、人工介质)，其共性在于：媒质的(复)介电常数是jω(j表示虚数单位，ω表示角频率)整数次幂的函数^[3]；第二大类主要包含Cole-Cole(C-C)模型^{[6, 7, 8]}(如生物组织、高分子材料)、Davidson-Cole(D-C)模型^[3](如甘油、丙二醇)和Havriliak-Negami(H-N)模型^[9](如聚合物)，其共性在于：媒质的(复)介电常数是jω分数次幂的函数^[3]。后文将以不同媒质对应的色散模型(或所属大类)命名以区分不同的媒质类型(如称色散特性满足C-C模型的媒质为C-C媒质)。

近几年来，由于媒质的色散特性引起了较多关注^{[10, 11]}，相应地，FDTD法的解算对象已由常规的非色散媒质^[12]拓展到色散媒质^{[6, 13]}。

FDTD建模第一大类(Debye、Drude、Lorentz)媒质，已发展了多种方法，主要有3种思路：1) 引入ADE^{[2, 7, 8, 9]}；2) 利用递归卷积(recursive convolution, RC)^{[2, 13]}；3) 定义移位算子(shift operator, SO)^{[2, 14]}。

然而，FDTD建模第二大类(C-C、D-C、H-N)媒质，仍是当前的难点之一，其主要困难是差分实现分数阶导数^[3]。对此，文献[15]进行了一系列开创性的基础工作：利用帕德(Padé)近似法，导出了一组整数阶ADEs，克服了分数阶导数困难，分别提出了处理(无磁、线性、各向同性的)单极、多极C-C媒质^{[7, 8]}、单极D-C媒质^[3]和单极H-N媒质^[9]的FDTD方案；一维的系列仿真结果初步证实了这些方案的可行性^{[7, 8, 9]}。文献[6]给出了建模多极C-C媒质的FDTD改进方案和三维的仿真测试。

为了增加通用性和降低复杂度，方便统一编程处理不同类型色散媒质的电波传播问题，本文主要以文献[7, 8, 9]为基础，提出了一种FDTD改进方案，创新(改进)之处主要体现在：1) 从单极色散媒质推广到多极情形；2) 从无磁媒质推广到有磁情形；3) 增加了静态电导率项；4) 增加了三维数值算例检验该方案的性能和精度。

1 问题描述

首先给出后文满足的假设条件和必要的符号说明：1) 假设所研究的色散媒质是各向同性、线性、有磁和电色散的，该假设条件经常应用于地球物理勘探、微波医学成像、红外与毫米波技术等领域^{[1, 2, 3, 4, 5]}；2) 符号顶部加点表示其值为复数。

若问题空间r处存在一般的(同时存在电极化损耗和欧姆损耗)多极H-N媒质，则频域本构关系为：

$ \mathbf{\dot{D}}(\mathbf{r},\omega )={{\varepsilon }_{0}}{{\dot{\varepsilon }}_{\text{r}}}(\mathbf{r},\omega )\mathbf{\dot{E}}(\mathbf{r},\omega ) $

(1)

式中， $\mathbf{\dot{D}} $和 $ \mathbf{\dot{E}} $分别表示电位移矢量和电场强度矢量；ω = 2πf，f为频率；ε₀为真空介电常数；多极H-N媒质的相对介电常数定义为^[16]：

$ {{\dot{\varepsilon }}_{\text{r}}}(\mathbf{r},\omega )\triangleq {{\varepsilon }_{\infty }}(\mathbf{r})+\frac{{{\sigma }_{\text{s}}}(\mathbf{r})}{\text{j}\omega {{\varepsilon }_{\text{0}}}}+\sum\limits_{w=1}^{W}{{{{\dot{f}}}_{w}}(\mathbf{r},\omega )\Delta {{\varepsilon }_{w}}(\mathbf{r})} $

(2a)

式中，ε_∞为光学相对介电常数；σ_s为静态电导率；Δε_w = ε_s,w−ε_∞；ε_s,w为第w极的静态相对介电常数；W为极总数；第w极的复值因子 $ {{\dot{f}}_{w}} $为：

$ {{\dot{f}}_{w}}(\mathbf{r},\omega )={{\left\{ 1+{{[{{{\dot{z}}}_{w}}(\mathbf{r},\omega )]}^{{{\alpha }_{w}}(\mathbf{r})}} \right\}}^{-{{\beta }_{w}}(\mathbf{r})}} $

(2b)

式中，α_w(0≤α_w≤1)；β_w(0≤β_w≤1)为第w极色散参数； $ {{\dot{z}}_{w}} $= jωτ_w(r)；τ_w表示第w极的弛豫时间。

本文采用的色散模型是较为通用的形式：1) 当W = 1时，简化为单极情形；2) 当W = 1且α_s = 0时，简化为文献[9]的式(1)；3) 如果W = 1、α_w = 1且σ_s = 0时，退化为文献[3]的式(1)；4) 如果β_w = 1且σ_s = 0时，退化为文献[8]的式(1)；5) 如果W = 1、β_w = 1且σ_s = 0时，变为文献[7]的式(1)；6) 如果α_w = 1且β_w = 1时，变为Debye模型^[4]；7) 假如α_w = 1、β_w = 1且Δε_w = 0时，转化为非色散情形。

2 辅助微分方程(ADEs)

为了克服FDTD方案中分数(α_w、β_w一般为分数)阶导数的困难，首先，借鉴文献[7, 8, 9]的研究思路，利用帕德(Padé)近似法^[15]，并由 $ {{\dot{f}}_{w}}$在 $ {{\dot{z}}_{w}} $= 1处解析，将 $ {{\dot{f}}_{w}} $近似为：

$ {{\dot{f}}_{w}}(\mathbf{r},\omega )\approx \frac{\sum\limits_{n=0}^{N}{{{p}_{n}}{{[{{{\dot{z}}}_{w}}(\mathbf{r},\omega )-1]}^{n}}}}{\sum\limits_{m=0}^{M}{{{q}_{m}}{{[{{{\dot{z}}}_{w}}(\mathbf{r},\omega )-1]}^{m}}}} $

(3a)

式中，N、M分别表示分子、分母多项式的阶数；p_n、q_m分别表示分子、分母多项式的系数，可通过求解线性方程组获得，相关细节参见文献[7, 8, 9]。

将式(3a)代入式(2a)，易得相对介电常数的帕德近似值为：

$ \begin{matrix} {{{\dot{\hat{\varepsilon }}}}_{\text{r}}}(\mathbf{r},\omega )\approx {{\varepsilon }_{\infty }}(\mathbf{r})+\frac{{{\sigma }_{\text{s}}}(\mathbf{r})}{\text{j}\omega {{\varepsilon }_{\text{0}}}}+ \\ \sum\limits_{w=1}^{W}{\Delta {{\varepsilon }_{w}}(\mathbf{r})\frac{\sum\limits_{n=0}^{N}{{{p}_{n}}{{[\text{j}\omega {{\tau }_{w}}(\mathbf{r})-1]}^{n}}}}{\sum\limits_{m=0}^{M}{{{q}_{m}}{{[\text{j}\omega {{\tau }_{w}}(\mathbf{r})-1]}^{m}}}}} \\ \end{matrix} $

(3b)

这样，相对介电常数已全部转化为jω整数次幂的函数，有利于后文FDTD方案的讨论。

其次，定义第w极(辅助的)频域电极化强度矢量为：

$ {{\mathbf{\dot{P}}}_{w}}(\mathbf{r},\omega )\triangleq {{\varepsilon }_{0}}\Delta {{\varepsilon }_{w}}(\mathbf{r})\frac{\sum\limits_{n=0}^{N}{{{p}_{n}}{{[\text{j}\omega {{\tau }_{w}}(\mathbf{r})-1]}^{n}}}}{\sum\limits_{m=0}^{M}{{{q}_{m}}{{[\text{j}\omega {{\tau }_{w}}(\mathbf{r})-1]}^{m}}}}\mathbf{\dot{E}}(\mathbf{r},\omega ) $

(4)

最后，应用傅里叶逆变换(inverse Fourier transform, IFT)，可以得到一组(W个)关于时域电极化强度矢量P_w的ADEs为：

$ \begin{align} & \sum\limits_{m=0}^{M}{{{Z}_{m}}\tau _{w}^{m}(\mathbf{r}){{{\partial }^{m}}{{\mathbf{P}}_{w}}(\mathbf{r},t)}/{\partial {{t}^{m}}}\;}={{\varepsilon }_{0}}\Delta {{\varepsilon }_{w}}(\mathbf{r})\times \\ & \ \ \ \ \ \ \ \ \ \ \ \sum\limits_{n=0}^{N}{{{Y}_{n}}\tau _{w}^{n}(\mathbf{r}){{{\partial }^{n}}\mathbf{E}(\mathbf{r},t)}/{\partial {{t}^{n}}}\;} \\ \end{align} $

(5)

其中：

$ {{Y}_{n}}=\sum\limits_{k=n}^{N}{{{(-1)}^{k-n}}{{p}_{k}}\frac{k!}{n!(k-n)!}} $

(6)

$ {{Z}_{m}}=\sum\limits_{k=m}^{M}{{{(-1)}^{k-m}}{{q}_{k}}\frac{k!}{m!(k-m)!}} $

(7)

式中，E表示时域电场强度矢量；t表示时间。至此，这组辅助微分方程由于包含了整数阶的时间导数，可以方便地通过有限差分法实现时间的离散化。

3 时域有限差分方案

前述电色散媒质中，时域(微分形式)的麦克斯韦(Maxwell)方程组为^[8]：

$ \nabla \times \mathbf{H}={{\varepsilon }_{0}}{{\varepsilon }_{\infty }}\partial \mathbf{E}/\partial t+{{\sigma }_{\text{s}}}\mathbf{E}+\partial \sum\limits_{w=1}^{W}{{{\mathbf{P}}_{w}}}/\partial t $

(8)

$ \nabla \times \mathbf{E}=-{{\mu }_{0}}{{\mu }_{\text{r}}}\partial \mathbf{H}/\partial t $

(9)

式中，H(为了表达简洁，省略了自变量，下文同)表示磁场强度矢量H(r,t)；μ_r表示相对磁导率，当μ_r = 1时，媒质简化为文献[7, 8, 9]的无磁特例。

在t=iΔt(i和Δt分别表示时间步的指标和步长)时刻的电场强度矢量Eⁱ，其n阶的时间偏导数可通过中心差分近似为^{[7, 8, 9]}：

$ {{{\partial }^{n}}{{\mathbf{E}}^{i}}}/{\partial {{t}^{n}}}\;\approx \sum\limits_{m=-{{l}_{n}}}^{{{l}_{n}}}{Q_{m}^{n}{{\mathbf{E}}^{i+m}}}/\Delta {{t}^{n}} $

(10)

式中， $ {{l}_{n}}=\left\lceil n/2 \right\rceil $表示对n/2上取整；Qn m为中心差分参数^{[6, 7, 8, 9]}；P_w的中心差分近似可类似地获得。至此，就实现了式(5)的FDTD时间离散化。

文献[9]显示，欲使FDTD方案稳定，需满足N≤M。为了简化描述，假定M = N(其原因后文给出)，则有l_M=l_N。进一步假定已知iΔt时刻的P_w、E和(i+0.5)Δt时刻的H，分别差分式(8)～式(9)，并引入辅助矢量Φ_w^[8]，经过整理，多极H-N媒质的FDTD时间迭代步进方案可简述如下：

1) 计算辅助矢量为：

$ {{\Phi }_{w}}={{({{\zeta }_{{{l}_{M}},w}})}^{-1}}\sum\limits_{k=-{{l}_{M}}}^{{{l}_{M}}-1}{({{\theta }_{k,w}}{{\mathbf{E}}^{i-{{l}_{M}}+1+k}}-{{\zeta }_{k,w}}\mathbf{P}_{w}^{i-{{l}_{M}}+1+k})} $

(11)

其中：

$ {{\theta }_{k,w}}={{\varepsilon }_{0}}\Delta {{\varepsilon }_{w}}\sum\limits_{n=\upsilon }^{N}{({{{Y}_{n}}\tau _{w}^{n}Q_{k}^{n}}/{\Delta {{t}^{n}}}\;)} $

(12)

$ {{\zeta }_{k,w}}=\sum\limits_{m=\upsilon }^{M}{({{{Z}_{m}}\tau _{w}^{m}Q_{k}^{m}}/{\Delta {{t}^{m}}}\;)} $

(13)

式中，k = −l_M,−l_M+1,…,l_M；ν = max{2|k|−1,0}。

2) 更新(i+1)Δt时刻的电场强度矢量为：

$ {{\mathbf{E}}^{i+1}}=\frac{\left( {{\varepsilon }_{0}}{{\varepsilon }_{\infty }}-\frac{\Delta t{{\sigma }_{\text{s}}}}{2} \right){{\mathbf{E}}^{i}}+\sum\limits_{w=1}^{W}{\mathbf{P}_{w}^{i}-\sum\limits_{w=1}^{W}{{{\Phi }_{w}}}+\Delta t\nabla \times {{\mathbf{H}}^{i+\frac{1}{2}}}}}{{{\varepsilon }_{0}}{{\varepsilon }_{\infty }}+\sum\limits_{w=1}^{W}{{{({{\zeta }_{{{l}_{M}},w}})}^{-1}}{{\theta }_{{{l}_{M}},w}}+\frac{\Delta t{{\sigma }_{\text{s}}}}{2}}} $

(14)

3) 更新(i+1)Δt时刻的电极化强度矢量为：

$ \mathbf{P}_{w}^{i+1}={{\Phi }_{w}}+{{({{\zeta }_{{{l}_{M}},w}})}^{-1}}{{\theta }_{{{l}_{M}},w}}{{\mathbf{E}}^{i+1}} $

(15)

4) 更新(i+1.5)Δt时刻的磁场强度矢量为：

$ {{\mathbf{H}}^{i+1.5}}={{\mathbf{H}}^{i+0.5}}-{\Delta t\nabla \times {{\mathbf{E}}^{i+1}}}/{({{\mu }_{\text{r}}}{{\mu }_{0}})}\; $

(16)

前述FDTD方案中，对于每个离散的网格，每个场分量的存储量为M(W+1)+W，磁场以外的计算时间(乘法运算次数)为(2M+1)W+3，当W=1时与文献[9]的结果相同。

该FDTD方案的时间迭代步进流程(技术路线)如图 1所示^[6]。

4 结果与讨论

由于前述FDTD方案中采用了帕德(Padé)近似，下面首先检验帕德近似的精度。

4.1 帕德近似的精度

定义帕德近似的相对均方(relative mean square, RMS)误差为^{[3, 7]}：

$ e({{\dot{\hat{\varepsilon }}}_{\text{r}}},{{\dot{\varepsilon }}_{\text{r}}})\triangleq \sqrt{\frac{\int_{{{f}_{\text{L}}}}^{{{f}_{\text{H}}}}{{{\left| {{{\dot{\hat{\varepsilon }}}}_{\text{r}}}-{{{\dot{\varepsilon }}}_{\text{r}}} \right|}^{\text{2}}}\text{d}f}}{\int_{{{f}_{\text{L}}}}^{{{f}_{\text{H}}}}{{{\left| {{{\dot{\varepsilon }}}_{\text{r}}} \right|}^{\text{2}}}\text{d}f}}} $

(17)

式中，媒质相对介电常数的解析值 ${{\dot{\varepsilon }}_{\text{r}}} $和Padé近似值 ${{\dot{\hat{\varepsilon }}}_{\text{r}}} $可分别通过式(2a)和式(3b)求得；f_H、f_L分别表示上、下限频率。

算例 1 以文献[9]中的均匀单极H-N媒质(聚丙烯酸正辛酯，聚丙烯酸的一种衍生物)为测试目标，各模型参数分别为W = 1、ε_∞ = 2.61、σ_s = 0.001 S/m、ε_s,1 = 3.88、τ₁ = 39.98 ps、α₁ = 0.73、β₁ = 0.66。上、下限频率分别取为f_H = 100 GHz、f_L = 10 MHz，分母多项式的阶数取为M = 3,4,…,10，分子多项式阶数分别取为N = M、N = M−1和N = M−2时，帕德近似的相对均方误差e随M的变化关系如图 2所示(误差轴为对数坐标，后文同)，其中当M=3、4、5、6时，本文的误差与文献[9]的误差对比如表 1所示。

分析图 2发现：1) 当M相同时，在N = M、N = M−1、N = M−2的3种情形下，第一种情形误差最小，这正是前文讨论FDTD时间步进方案时选取N = M的依据；2) 在N = M、N = M−1、N = M−2的3种情形下，随着M的增加，误差均逐渐减小(然而，实现FDTD方案所需的存储量M(W+1)+W和计算时间(2M+1)W+3也同时增加，建议折中选择M、N)。

结合表 1可知：1) 在表中所列举的12种情形下，本文的误差均略大于文献[9]的误差(适当增加离散频点数可使二者相当)，因此，文献[9]的单极情形可视为本文改进方案的一个特例；2) 当N = M = 4时，本文的误差约为0.003 6，不仅计算精度能被一般的工程应用所接受，而且计算速度比较适中。

算例2 (为不失一般性)选取一种2极(即W=2)的均匀H-N媒质为测试目标，其余模型参数分别为ε_∞ = 2.61、σ_s = 0.01 S/m、ε_s,1 = 3.88、ε_s,2 = 3.50、τ₁ = 39.98 ps、τ₂ = 29.98 ns、α₁ = 0.73、α₂ = 0.80、β₁ = 0.66、β₂ = 0.70。频率上、下限及多项式的阶数均与算例1相同，相对均方误差e随M的变化关系如图 3所示。

由图 3可见：对于多极H-N媒质，也有类似于算例1(单极情形)的结果，只是误差略微增加。当N = M = 4时，相对均方误差约为0.004 3，还能满足一般的工程应用需求，因此，在后文的FDTD方案中，多项式的阶数也默认这组取值。

4.2 FDTD方案的可行性

为了检验前述FDTD改进方案的可行性，再给出两个算例，选择算例基于以下考虑：1) 问题空间涵盖不同物理维度；2) 色散媒质涵盖不同类型；3) 媒质涵盖不同极数；4) 色散模型参数涵盖文献作者自拟和真实测量；5) 测试采用不同类型目标参数。

算例 3 出自文献[3]的一维问题：单极均匀无磁D-C媒质(媒质参数分别为W = 1、μ_r = 1、ε_∞ = 2、σ_s = 0 S/m、ε_s,1 = 50、τ₁ = 153 ps、β₁ = 0.9和α₁ = 1.0)充满z≥0半空间，其余空间为真空；一列入射平面波(电场沿x方向极化)沿+z方向传播，平面波由下述的超宽带(ultra-wideband, UWB)调制高斯脉冲源s激励产生^[3]：

$ s(t)={{\text{e}}^{-{{a}^{2}}{{(t-4/a)}^{2}}}}\sin (2\pi {{f}_{\text{c}}}(t-4/a)) $

(18)

式中，参数a =1.26×10¹⁰ s^-1，中心频率f_c = 6 GHz，该脉冲的频谱主要涵盖0.1～10 GHz频率范围。

FDTD离散化的空间、时间步长分别取为Δz = 1.1 mm、Δt = 1.5 ps，离散网格数、时间步数分别取为n_z = 200、n_t = 3 000。仿真区域−z一侧采用简单的吸收边界：E_x(z_max, t+Δt) = E_x(z_max−Δz, t)，其中z_max表示对应的边界位置坐标。

记z、z+d分别表示D-C媒质中的两个不同位置，FDTD仿真并且存储这两处的时域电场E_x(z, t)、E_x(z+d, t)，后处理(傅里叶变换)获得相应的频域电场 ${{\dot{E}}_{x}}$(z, ω)、 $ {{\dot{E}}_{x}} $(z+d, ω)，则色散媒质相对介电常数的FDTD估计值为^{[3, 7]}：

$ {{\dot{\tilde{\varepsilon }}}_{\text{r}}}(d,\omega )=-{{({{{c}_{0}}({{\gamma }_{\text{R}}}+\text{j}{{\gamma }_{\text{I}}})}/{\omega }\;)}^{2}} $

(19)

式中，c₀表示真空光速；参数γ_R、γ_I分别为：

$ {{\gamma }_{\text{R}}}(d,\omega )={\ln (|{{{{\dot{E}}}_{x}}(z+d,\omega )}/{{{{\dot{E}}}_{x}}(z,\omega )}\;|)}/{d}\; $

(20)

$ {{\gamma }_{\text{I}}}(d,\omega )={(\arg \{{{{\dot{E}}}_{x}}(z+d,\omega )\}-\arg \{{{{\dot{E}}}_{x}}(z,\omega )\})}/{d}\; $

(21)

选取d = 20Δz时，复值相对介电常数实部和负虚部的解析值、Padé近似值、FDTD估计值随频率的变化关系如图 4所示(频率轴为对数坐标，后文同)。

从图 4可以看出，在UWB条件下，以一维问题的相对介电常数为测试目标，其FDTD估计值和Padé近似值或解析值之间差异不大(实际结果显示相对均方误差均小于1%)；再与文献[3]中图 4a的结果对比发现，二者也基本保持一致，因此，文献[3]的单极情形可视为本文改进方案的又一个特例。

算例 4 选取一个更具代表性的三维问题：计算空气(视为真空)中均匀4极C-C弱磁性介质球的后向雷达散射截面(radar cross section, RCS)^[2]，球半径为50 mm，由浸润型脂类组织构成，其模型参数(源于实验测量^[17])分别为ε_∞ = 2.5、σ_s = 0.035 S/m、Δε₁ = 9.0、Δσ₂ = 35、Δε₃ = 3.3×10⁴、Δε₄ = 1.0×10⁷、τ₁ = 7.96 ps、τ₂ = 15.92 ns、τ₃ = 159.15 μs、τ₄ = 15.915 ms、β₁ = 0.80、β₂ = 0.90、β₃ = 0.95、β₄ = 0.99；其余模型参数分别为μ_r = 1.2、W = 4、α₁ = α₂ = α₃ = α₄ = 1.0；采用的入射平面波与算例3相同。

FDTD离散化的空间网格尺寸、时间步长分别取为Δ = Δ_x = Δ_y = Δ_z = 3 mm、Δt = Δ/(2.c₀)，离散获得的网格数、时间步数分别为n = 50×50×50、n_t = 5 000。主FDTD区域周围用10层卷积完全匹配层(convolution perfectly matched layer, CPML)^[2]吸收边界截断。

作为比较，RCS的解析值由Mie级数法^[2]获得，两种方法的计算结果对比如图 5所示。

通过图 5发现，在较宽泛的频率范围，以三维吸波介质球的后向RCS为测试目标，其FDTD估计值和解析值仍然符合较好(结果显示，相对均方误差较一维情形有所增加，原因可能是网格剖分粗糙所致，有望通过先进的网格剖分技术进一步减小误差^[2])。

综合以上4个算例可看出，利用前述FDTD改进方案模拟UWB电磁波在C-C、D-C、H-N等几类电色散媒质中传播，估计相对介电常数、后向雷达散射截面所产生的误差，一般的工程应用可以接受。这些仿真结果初步证实，利用本文发展的FDTD方案，模拟一般电色散媒质中的电波传播是可行的。

5 结 束 语

本文改进了文献[7, 8, 9]提出的FDTD方案，改进后的方案具有适用范围广、实现复杂度低等优势，适用于统一处理射频、微波、红外和毫米波等工程中一般电色散媒质的电波传播问题，媒质可以是一维或多维的、单极或多极的、无磁或有磁的。算例的初步结果显示，在超宽带频谱范围，该方案具有较高的数值精度。该方案的鲁棒性还需经历实验测量、噪声影响等因素的检验。