有耗对称介质层反射系数的快速计算

孙宪钢，魏 兵，石 磊

(1.西安电子科技大学 物理与光电工程学院，陕西 西安 710071；2.西安电子科技大学 信息感知协同创新中心，陕西 西安 710071)

平面层状结构在生活中十分常见，除了美观之外，它的电磁传输特性也十分显著，尤其是周期平面层状结构在吸波材料[1-4]、频率选择表面[5]、地质构造研究[6]和光子晶体[7-8]等方面都有着广泛的应用。分析电磁波在层状介质中的反射和透射，常采用传播矩阵法[9-13]。在20世纪60年代，DE JAGER等[14]用传播矩阵法对非均匀层状结构的反射和透射进行了探究。2000年，郑宏兴等[10]利用传播矩阵法计算了各向异性分层介质的反射和透射。2017年和2019年，王飞等[12-13]用传播矩阵法分别对分层有耗手征介质和含石墨烯分界面有耗分层介质的反射和透射进行了计算。可以看到，传播矩阵法应用较为广泛，但是，该方法以矩阵连乘的形式处理层状结构，当层数较多时，需要考虑矩阵连乘的数量以及计算结果的稳定性问题。特别是传播矩阵法在计算有耗层状体系时，如果材料电导率较大或结构层数较多，则容易产生较大的计算量并出现数值失稳问题，导致计算结果发散[15]。

有耗对称的平面介质层，作为一种特殊的层状结构，在以上各方面也都有应用，可采用传播矩阵法计算其反射系数，但是也会出现相同的问题。为此，采用等效四端网络法计算有耗对称平面介质层的反射系数。1990年，李世智等[16]在研究地表探测时，对此法进行了初步的研究。1998年，段方勇等[17]对李世智等的研究进一步加深，并且给出了一些理论推导。2007年，MEHTA等[18]在近地表层特性研究时也运用此法做了一定的运算。笔者通过对等效四端网络法进行详细研究，发现该方法对界面矩阵进行了简化处理，相对于传播矩阵法，降低了矩阵的数量，避免了传播矩阵法中出现的大量矩阵连乘，在计算多个频点或宽频段时可以快速计算有耗对称介质层的反射系数，并且以光子晶体为例，对比了两种方法在具体计算其反射系数时的计算时间和消耗的内存。

1 理论公式

1.1 有耗介质中拉普拉斯波动方程解

设平面电磁波沿z轴正方向传播，其中电场强度为Ex，磁场强度为Hy，则可以得到其一维有耗介质的Maxwell支配方程为

(1)

其中，ε为介电常数，μ为磁导系数，σ为电导率，σm为导磁率。为了得到频域解析解，对式(1)的时域t进行拉普拉斯变换到s域[19]，化简求解可以得到电磁场满足

(2)

其中，γ(s)=((σm+μs)(σ+εs))1/2，为损耗介质中的电波传播系数；B(s)和A(s)为上下行波幅值，由边界条件确定。令s=jω，即可得到频域解析解，再令σ=0，式(2)退化为一般无耗介质波动方程的解。

1.2 等效四端网络法求反射系数

设平面电磁波从真空(ε0，μ0)垂直入射有N层介质的第1层介质的分界面上，入射平面为xOz平面，如图1所示。各层介质是均匀有耗的，介质分界面都为平行平面，第n层介质的厚度为dn，相应的电磁参数为：介电系数εn，磁导系数μn，电导率σn，导磁率σmn，n=1，2，…，N。

图1 分层介质

设第ln个界面的电场强度和磁场强度分别为En和Hn，其表达式可由式(2)获得。由于在层状结构的输入界面l1和输出界面ln+1各有2个变量，因此，此层状介质系统可以等效于有2个输入量和2个输出量的四端网络，这个四端网络的参数可以通过矩阵[Uij](i=1，2；j=1，2)联系，所以，系统的左右边界的矩阵表达式为

(3)

式(3)中矩阵元素Uij的值与所有介质层的物理参数有关。假设任意的第n层的左右边界面的关系矩阵为[uij](i=1，2；j=1，2)已知，则应用n次递推，式(3)可以改写为

(4)

由式(3)和式(4)可以得到

(5)

式(5)表示出了四端网络的参数由每一层的界面关系矩阵决定，其中，第n层介质中两界面间的关系矩阵[uij](i=1，2；j=1，2)可以表示为

(6)

对于关系矩阵[uij](i=1，2；j=1，2)，可参考文献[17]。由式(6)可知，当第n层介质的厚度为0时，两界面间的关系矩阵[uij](i=1，2；j=1，2)退化为单位矩阵，即介质层面之间的矩阵为单位矩阵。其中γn=((σmn+μns)(σn+εns))1/2，为每一层的损耗介质的电波传播系数，bn=γn(s)/(σn+εns)为计算的中间变量。所以输入界面和输出界面的关系矩阵[Uij](i=1，2；j=1，2)可以通过每层左右边界面的关系矩阵[uij](i=1，2；j=1，2)求出。入射波Ei，反射波Er和透射波Et的表达式可以表示为

(7)

(8)

将式(8)的总入射场和总透射场表达式代入式(3)，可以求出反射系数R为

(9)

可以看出，已知系统四端网络的参数、前后区域的电磁参数以及该区域的电波传播常数，即可用式(9)计算有耗分层介质的反射系数，以上参数的计算可以通过式(5)和式(6)得到。

2 算 例

算例1有耗3层介质的反射。设介质层的厚度为：d1=d3=0.5 m，d2=0.2 m；相应的电磁参数为：ε1=ε3=4ε0，ε2=ε0，μ1=μ2=μ3=μ0，σ1=σ3=0.05 S/m，σ2=0.03 S/m，下标代表层数编号；入射区域和透射区域为空气。设定入射平面波的频率为0.000 1～1 GHz，采用等效四端网络法计算反射系数，并给出传播矩阵法和FDTD方法[20-21]的计算结果。由图2可知，三者符合良好，可以看出，3层有耗介质层的反射呈现单调性变化。

图2 3层有耗介质的反射系数

算例2有耗4层介质的反射。设介质层的厚度为d1=d4=0.5 m，d2=d3=0.2 m；相应的电磁参数为：ε1=ε4=4ε0，ε2=ε3=2ε0，μ1=μ2=μ3=μ4=μ0，σ1=σ4=0.005 S/m，σ2=σ3=0.02 S/m，下标代表层数编号，入射区域和透射区域为空气。设定入射平面波的频率为0.000 1～1 GHz。采用等效四端网络法计算反射系数，并给出传播矩阵法和FDTD方法的计算结果。由图3可知，三者符合良好，可以看出，此时有耗介质层的反射在高频情形呈现出周期性变化，在低频时呈现单调性变化(类似于算例1)。

图3 4层有耗介质的反射系数

图4 光子晶体的反射系数

图5 光子晶体

由图2和图3可知采用等效四端网络法计算有耗对称介质层的反射系数，无论是奇数层或偶数层，所得结果均与传播矩阵法及FDTD方法的计算结果较好吻合。将文中采用的等效四端网络法与传播矩阵法对比，给出算例1、算例2和算例3中具体的计算参数，如表1和表2所示。

表1 算例1和算例2等效四端网络法和传播矩阵法对比

表2 算例3等效四端网络法和传播矩阵法对比

由表1可知，算例1、2中采用等效四端网络法和传播矩阵法计算有耗对称介质层的反射系数时，所需内存相差不大，但是传播矩阵法所需的计算时间约为等效四端网络法的3.4倍。等效四端网络法的计算效率高于传播矩阵法。由表2可知，与传播矩阵法相比，算例3中等效四端网络法在计算时间上减少90%以上，而且传播矩阵法所用时间已超1小时以上。由此可见，当层数越多时，计算频段越宽或采样频点较多时，等效四端网络法的计算效率将远高于传播矩阵法的计算效率。

3 总 结

在计算有耗对称介质层的反射系数时，结合传输线理论和传播矩阵的思想，笔者采用了等效四端网络法计算其反射系数，并且与传播矩阵法在计算效率上进行了比较。相对于传播矩阵法，等效四端网络法对界面矩阵进行了简化处理，使得矩阵计算时赋值次数成倍减少，降低了计算量，不论是周期还是非周期的介质层都可以快速计算其反射系数，掌握其反射特性。从仿真算例可以看出，如果介质层层数越多，或者计算频段越宽，等效四端网络法的计算效率将远高于传播矩阵法，同时FDTD方法的验证表明了计算结果的准确性，为其在实际应用和设计中提供了一种快速计算电磁反射系数的方案。