好。仿真结果和解析结果的误差主要是由于数值计 算方法造成的，加密离散网格可有效降低数值误 差，然而这将极大地增加数值计算的复杂度。
10
5
结论
8
sample points AFS_SB_MLFMA MIE
RCS(dBsm)
6
4
2
针对插值模型比较法和二分法自适应采样策略 的不足，本文提出了一种新的基于 Stoer-Bulirsch 插 值的自适应采样方法。通过比较 n 和 n+1 个采样点 的插值结果，以及对插值区域二分加密后继续进行 收敛性判断，同时实现了采样数目的自适应增加和 插值区间的自适应加密，较好地避免了伪收敛情况 和增强了算法的稳健性，达到了以最少的采样点来 插值计算未知函数响应的目的。
0 100 150 200 250 300
Frequency(MHz)
图 2 半径 1m 金属球后向 RCS 对比
参考文献
[1] [2] [3] [4] Yang Ding, Ke-Li Wu, and Da Gang Fang, A broad-band adaptive frequency-sampling approach for microwave-circuit EM simulation exploiting Stoer-Bulirsch algorithm, IEEE Trans. Microwave Theory Tech., vol.51, pp. 928-934, March 2003. Karwowski, Interpolation broadband shielding behavior of wire-grid cages from full-wave electromagnetic simulation, IEEE International Symposium on Electromagnetic Compatibility, EMC Europe 2012, pp. 1-4, Sep. 2012. J. Reddy, M. D. Deshpande, F. B. Beck. Fast RCS computation over a frequency band using method of moments in conjunction with asymptotic waveform evaluation technique, IEEE Trans. Antenna and Prop., vol. 46, pp.1229-1233. August 1998. Stoer, J., Bulirsch, R. Introduction to Numerical Analysis, New York: Springer-Verlag, pp. 58-72, 1980.
作者简介： 何十全，男，讲师，主要研究领域为计算电磁学，电磁兼容仿真与应用等；韩奎，男，硕士研究生， 主要研究领域为计算电磁学；聂在平，男，教授、博士生导师，主要研究领域为计算电磁学，非均匀介质 中的场与波等。
·577·
基于 OpenMP 与 VALU 加速技术的多层快速多极子混合并行算法
刘金波 何 芒
(北京理工大学信息与电子学院，北京 100081) hemang@https://www.doczj.com/doc/925597785.html,
摘 要：研究了在共享内存式计算机架构下，基于 OpenMP 及 VALU 硬件加速的多层快速多极子的混合并 行算法。着重讨论了多层快速多极子并行程序设计中的几个关键问题，并详细分析了混合并行计算中的关 键点。针对一些典型散射目标，给出了由混合并行程序得到的雷达散射截面数值结果。通过对计算数据的 对比分析，得到具有实际指导意义的相关结论。 关键词：OpenMP，VALU，多层快速多极子，雷达散射截面
A Hybrid Parallelization Technique Based on OpenMP and VALU Acceleration for the Multilevel Fast Multipole Algorithm LIU Jinbo, HE Mang
(School of Information and Electronics, Beijing Institute of Technology, Beijing 100081) Abstract: In the framework of the shared-memory computer system, a hybrid parallelization technique based on OpenMP and VALU acceleration for the Multilevel Fast Multipole Algorithm (MLFMA) is proposed. In addition, several key issues in the hybrid parallel MoM code design and the major factors influencing the parallel efficiency are addressed in detail. Numerical results of the radar cross section of typical objects are shown, through which comparative analyses are performed to give some concluding remarks. Keywords: OpenMP; VALU; MLFMA; RCS
1
引言
近年来，计算机技术取得了迅猛发展，一台普 通的计算机即可达到 4 核甚至更高的配置水平。从 Intel Pentium IV 开始，每个 CPU 核心都具有一个 浮 点 运 算 单 元 （ FPU ） 和 矢 量 算 术 逻 辑 单 元 （VALU） 。与 FPU 每次只能对一个数据进行操作不 同， VALU 可以同时对 4 个数据进行操作 [1] 。然 而，在计算电磁学中，传统的串行程序仅利用了 CPU 某一个核心中的 FPU 进行数值计算，而其余 计算机资源都处于闲置状态，这就造成了极大浪 费。因此，在现代计算机架构下，如何充分利用计 算资源进一步挖掘计算潜力，成为了计算电磁学领 ·578·
域中新的热点与挑战之一。本文研究了在共享内存 式计算机架构下，基于 OpenMP 及 VALU 硬件加速 的多层快速多极子（ MLFMA ）的混合并行算法， 充分利用了现代计算机的多核结构与矢量计算能 力，大大提高了程序的计算效率。
2
2.1
并行程序的设计
OpenMP 程序设计中的关键问题
由于近场区矩阵填充的并行设计十分简单，故 此着重讨论 MLFMA 中矩阵乘向量（MVP）的远相 互作用。主要通过三个步骤实现：聚集、转移、发 散。其中，聚集和发散过程的程序代码主要包含三
、对子盒子的循环 重循环：对父盒子的循环（ FL ） （SL ）以及对平面波的循环。并行程序设计中，若 将 OpenMP 的指令直接应用到最外层的循环体 FL 上，随着计算机线程的增加，并行 效率将会急剧下降。为了取得更高的并行效 率，我们需要调整循环次序。最佳的循环结构是将 SL 设置为最外层循环。此时，由于子盒子的数目 要远远多于父盒子的数目，和原始程序相比，改进 后的程序更容易实现负载均衡。 转移过程的程序代码同样包含三重循环：盒子 循环及每个盒子的次相邻盒子循环、平面波循环。 由于盒子及其次相邻盒子处于同一层中，因此应用 到聚集和发散过程中的循环策略并不适用于转移过 程。实验数据表明，在低层盒子中，平面波循环作 为最外层循环时可以获得更高的并行效率以及更短 的并行时间，而在高层盒子中保持原始程序的循环 结构更为合适。此时需要选择一个过渡层，并在高 层和低层采用不同的循环策略。然而，过渡层的选 择取决于层数以及最高层盒子的尺寸，若选择不当 将会大大影响并行计算的效率。大量的实验数据表 明，当盒子的数目满足分配到每个核中的循环体数 目大于 60 时，这一层可以作为过渡层。
TFQMR 求解器进行迭代求解，收敛误差设置为 0.001。以理想导体球的双站雷达散射截面（RCS） 的计算为例，分析混合并行程序的效率和精度。导 体 球 采 用 CFIE 进 行 求 解 ， 直 径 为 100 个波 长 （λ） ，共包含 11178357 个未知量。在计算中，分别 应用传统串行程序(Serial)、仅应用 OpenMP 的并行 程序(OpenMP)、仅应用 VALU 硬件加速技术的并 行程序(VALU)以及同时应用了 OpenMP 和 VALU 的混合并行程序 (OMP&VALU) 这四种程序进行计 算。为了分析程序的并行加速效率，分别定义加速 比 S 和混合并行效率 η 如下：
S=T
Tv
(1) (2)
η = S p × 100%
在公式（ 1 ）中，当下角标“ v ”用“ f ”代替 时，表示的是填充近场矩阵时的计算时间；当其用 “m”代替时，表示 MVP 中每一步迭代所用的计算 时间。Tv 和 T 分别表示相应的并行程序和串行程序 的计算时间。在公式（2）中，p 表示的是参与并行 计算的计算机的核数。表 1 列出了应用四类程序分 别计算导体球的 RCS 时的详细计算时间及加速性 能。
90
2.2
VALU 硬件加速技术
80
VALU 硬件加速技术对程序的性能有着相当大 的改进，理论上可以将速度提升 4 倍，但实际加速 能力还取决于问题的种类。然而在任何情况下，应 用了 VALU 硬件 加速 技术的 程序 一定比 不使 用 VALU 时速度要快[4]。同 OpenMP 相比，VALU 加 速技术并行粒度更细，对程序细节的要求也更高。 一般而言，只有以下几种情况可以使用此技术：(1) 纯粹的赋值语句；(2) 多重循环的最内重循环；(3) 循环体内的每一个语句能够独立执行； (4) 对于那 些需要退出条件的循环，必须保证只有一个入口和 一个出口； (5) 循环中除编译器的标准内部函数 外，不能有任何的函数调用。
70
Bistatic RCS (dBsw)
Mie series Serial OMP&VALU
60
50
40
30
20 160
165
170
175
180
θ(deg)
图 1 导体球的双站 RCS 计算结果
表 1 四类程序的计算时间及加速性能
Serial OpenMP VALU OMP&VALU Serial OpenMP VALU OMP&VALU
Tf (s)
6639.5 601.5 1990.3
S
— 11.0 3.3
η(%)
— 92.0 —
3
数值实验
Tm (s)
465.6 42.6 146.0
S
— 10.9 3.2
η(%)
— 90.8 —
为了验证上述关于 OpenMP 与 VALU 混合并行 编 程 理 论 的 正 确 性 ， 编 写 了 基 于 Intel Visual FORTRAN 的 MLFMA 混合并行程序。计算中采用
·579·
图 1 给出了导体球的双站 RCS 的计算结果及 Mie 级数的解析解。从图中可以看到，并行程序的 计算结果同串行程序相比十分吻合。图 2 给出了 VALU 的加速比同并行计算中所用核数的关系。从 图中可以看到，VALU 的加速比相当稳定，并不随 着核数的增加而改变。图 3 给出了 OpenMP 及 OMP&VALU 程序在填充近场区矩阵及迭代计算时 的并行效率。从图中可以看到，VALU 硬件加速技 术的应用并没有对 OpenMP 的效率产生太大影响。
4.0 4.0
100 95 90
100 95 90
Parallel efficiency (%)
85 80 75 70 65 60 55 50
OpenMP NIM filling OMP&VALU NIM filling OpenMP MVP OMP&VALU MVP
85 80 75 70 65 60 55 50 12
1
2
3
4
5
6
7
8
9
10
11
Number of threads
图 3 不同程序的并行效率
3.5
3.5
Speedup of VALU
3.0
3.0
2.5
NIM filling MVP
2.5
4
结论
2.0
2.0
1.5
1.5
1.0
1
2
3
4
5
6
7
8
9
10
11
1.0 12
Number of threads
图2
VALU 的加速比
本文成功地实现了基于 OpenMP 及 VALU 硬件 加速技术的 MLFMA 混合并行程序，充分利用了计 算机的多核资源。为了提高并行效率、降低计算时 间，我们应用了调整循环次序、寻找合适的过渡层 等并行计算的策略。从文中给出的数值实验结果可 以看出，此混合并行程序具有很高的计算效率和精 度，适用于任意形状导体目标的电磁散射计算，并 且具有良好的可扩展性。
参考文献
[1] [2] [3] [4] [5] [6] [7]
余文华，李文兴. 高等时域有限差分法[M], 哈尔滨工程大学出版社, 2011: 23-39
Sanjay Velamparambil, Weng Cho Chew. Analysis and Performance of a Distributed Memory Multilevel Fast Multipole Algorithm[J]. IEEE Antennas and Propagation Magazine, August 2005, 53(8): 2719–2727. Kalyan C. Donepudi, Jian-Ming Jin, Sanjay Velamparambil, et al. A Higher Order Parallelized Multilevel Fast Multipole Algorithm for 3-D Scattering[J], IEEE Antennas and Propagation Magazine, July 2001, 49(7): 1069–1078. Wenhua Yu, Xiaoling Yang, Yongjun Liu, et al. New Development of Parallel Conformal FDTD Method in Computational Electromagnetics Engineering[J], IEEE Antennas and Propagation Magazine, June 2011, 53(3): 15–41.
余文华, 杨小玲, 刘永俊等,并行 FDTD 和 IBM BlueGene/L 巨型计算机结合求解电大尺寸的电磁问题[J], 电波科学学报,
August 2006, 21(4): 562–566. R. F. Harrington. Field Computation by Moment Methods [M]. New York: MacMillan, 1968: 1–25 X.M.Pan, W.C.Pi, and X.Q.Sheng, “On OpenMP parallelization of the multilevel fast multipole algorithm”, Progress In Electromagnetics Research, vol. 112, pp. 199-213, Jan. 2011.
作者简介： 刘金波，男，博士生，主要研究领域为计算电磁学的并行计算；何芒，男，副教授、博士生导师，主 要研究领域为天线理论与设计、计算电磁学等。
·580·
RWG-SED 算法改进及其在辐射计算中的应用
陈 涛 张志衡
(北京无线电测量研究所，北京 100854) chentao---111@https://www.doczj.com/doc/925597785.html,
摘 要：为解决传统 RWG 基方法在求解大型周期阵列天线矩量法方程中出现的计算资源困难，本文运用
子域全域基方法 SED，降低了阻抗矩阵维数，较大程度地缓解了矩阵方程求解压力。为进一步提高阵中各 单元输入阻抗的计算精度，本文提出一种多尺度子域全域基方法 MSED。数值计算表明，这种方法具有提 高 SED 方法计算精度、加速阵列天线电流分布计算的作用，为分析大型有限周期阵列天线提供了途径。 关键词：子域全域基 SED，有限周期阵列天线，RWG 基函数
Analysis of large-scale periodic arrays of crossed dipoles based on RWG-MSED basis function Chen tao, Zhang zhiheng
(Beijing Institute of Radar Measurement, Beijing 100854) Abstract: The conventional approaches of MoM based on RWG basis functions suffer a lot in solving large-scale periodic antennas. In this article, the sub-entire domain basis function is defined on each single cell of the periodic array, considerably reducing matrix size, CPU time and memory cost. For further improving the calculation accuracy of input impedance, a multi-scale sub-Entire domain basis function has been proposed. Numerical results from the new method prove that compared with SED basis function, the proposed method has excellent agreements with the conventional RWG method of moments, and obtains the electric current rapidly, providing the approaches to finite large-scale periodic antennas. Keywords: Sub-entire domain basis function, Finite periodic array，RWG basis function
1
引言
近年来，人们提出了许多用于快速精确求解大
解决计算与存储的困难。但在处理子阵耦合关系中 采用迭代方式，计算较耗时。 CBF[2] 求解了一个规 则排布的微带天线阵，但采用了一种新的高阶基函 数计算单元耦合，效果较好。论文[3]提出一种以单 元为子域、构建有限阵列上的电流分布的子域全域 法 SED，在计算周期结构的频率选择表面 RCS 问 题上得到了较好结果。论文[4
，5]
型阵列天线的方法。其中一类是基于阵列天线周期 结构来减少问题未知数的全域基方法。它们旨在将 一个有限阵列天线上的电流分解为一系列有限的电 流分布形式。这些电流分布形式可以通过精确计算 较小阵列得到。这些方法中，比较有代表性的是宏 基函数 MBF 与特征基函数法 CBF。MBF
[1]
扩展的 SED 基函数
作用域，进一步减小了阵列未知数个数。为提高全 域基函数构建阻抗矩阵时矩阵单元的填充效率，论 文[6]在论文[3]SED 方法的基础上结合 FFT，将阻抗 ·581·
在计算
大型平面周期天线阵中，将大阵划分成多个子阵以
矩阵的填充速度从 O(n ) 提高到 O(n) ， n 为 SED 基函数个数。论文 [7] 改进建立 SED 基函数中的
2
总的基函数个数为 MN。对于阵中第 n 个单元，定 义其上的子域全域基函数 ?∑ M I mn f mn (r ), r ∈ Sn ? g n (r ) = ? m=1 0 , r ? Sn ? ?
Green 函数积分项，求解无限大地平面上半空间中
周期结构的 RCS，同样得到了较好结果。 本文对 SED 这类全域基函数方法进行了详细 验证，将计算平面周期结构 RCS 的算法思想运用 到求解三维有限周期振子阵的辐射问题中，较大程 度地缓解了矩阵方程求解压力。为进一步提高阵中 各单元输入阻抗的计算精度，提出一种多尺度子域 全域基函数 MSED 方法。数值计算结果表明，当单 个天线全域细化分块数并不显著时，天线各端口输 入阻抗计算精度得到较好改善。这种方法具有降低 矩阵维数、加速计算电流分布的作用，为分析大型 周期排列阵列天线提供了思路。
（1）
其中 I mn 为各天线单元全域电流展开系数，
f mn (r ) 为 RWG 基函数。这样，通过 SED 基函数
建立导体表面上电流分布式
J ≈ ∑ i =1α i gi (r )
N
（2）
求解 I mn 的方法是通过引入阵中特定位置单元 附近的 L 个伴随单元，精确计算这（L+1）单元小 阵的电流系数 I1 , I 2 I ( L +1) ，并用有限的电流系数
I1 , I 2
I ( L +1) 建立阵列对应位置上单元的电流初级
分布 I mn 。数据实验以及论文[3]均表明，当 L 取
2
2.1
算法理论
子域全域基 SED 方法
传统矩量法采用 RWG 基函数构建阵列天线各
8，即采用 9 单元小阵就可以模拟周期排布的大
阵，具体对应关系如图 2 所示。
单元电流分布，是一种全波子域基函数方法，计算 精确但不便于构造大型天线阵列。 子域全域基函数法 SED 是一种定义在阵中各 天线单元上的全域基函数；同时，对于整个阵列来 说，每个天线单元又可看作构成整个阵列的一个子 域。SED 利用近似的表面电流分布来定义单个天线 单元电流展开系数。具体如下：
图 2 九单元小阵模拟大阵示意图
同样采用 SED 基函数，建立对电场积分方程 式 3 的 N 次抽样检验，得到 SED 矩量法方程
Zα = β
右端项
（3) （4）
β n = ? ∫ g n ( r ) E in ( r )dr
Sn
Zpq = jwμ∫ +
Sp Sq
∫
gp (r )gq (r ')G(r ?r ')dr ' dr
1 ?S ? gp (r )?'S ? gq (r )G(r ?r ')dr ' dr jwε ∫Sp ∫Sq
M M m=1 m'=1 Spm Sqm' qm'
= jwμ∑∑I pmIqm' ∫ +
∫
f (r ') f pm(r )G(r ?r ')dr ' dr
图 1 周期排布交叉阵子矩形阵列
1 M M ∑∑IpmIqm' ∫Spm ∫Sqm' ?'S ? fqm' (r ')?S ? fpm(r )G(r ?r ')dr ' dr jwε m=1 m'=1
N N M M
构建如图 1 中二维平面中周期排布交叉阵子矩 形 阵 列 ， 单 元 总 数 N = N x N y ， N x、N y 分 别 为
= ∑p=1∑q=1∑m=1∑m'=1 I pmIqm'Z( pm,qm')
x、y 方向的单元数。
当单个天线剖分为 M 个 RWG 基函数时，阵列 ·582·
其中 n = 1, 2
N ， p, q = 1, 2
'
（5）
N ， Z ( pm , qm ')
是阵中第 p, q 号天线单元上第 m, m 号 RWG 基函
数建立的阻抗元素。将式 4、5 带入式 2，即可求得 各单元电流分布。与 RWG 基函数相比，SED 法矩 阵方程未知数大幅度减少。
具有的 RWG 基函数个数为 M=70，阵列总未知数 为 70×225=15750。 分别建立不同分块方式下的 MSED(n)，n 为单 个振子上的全域分块数。计算得阵列中心一行各单 元输入阻抗分布如图 4 所示
2.2
多尺度子域全域基函数 MSED
SED 在计算周期阵列结构的方向图或 RCS 上
较精确[3-7]，而在阵列天线分析中，同样关心的是各 单元端口匹配特性。数值计算表面，SED 方法在计 算周期天线各单元端口输入阻抗上存在一定误差， 尤其是当波束扫描时，这种误差将被放大。 为提高 SED 方法输入阻抗计算精度，下面采 用缩小单个天线上全域基范围的思路，将各天线单 元进行不等密度的全域划分，建立多尺度划分全域 基函数 MSED。以单个交叉振子为研究对象，将其 划分成多个全域区域以建立各区域全域基 MSED， 如图 3 所示
图4
225 元阵列中心行 y 向振子端口输入阻抗比较
表 1 15×15 阵列计算时耗分析
方法 RWG SED MSED(5) MSED(7)
未知数 15750 450 2250 3150
计算 小阵 无 21s 21s 21s
建立全域 基矩阵 无 9.5s 11.9s 13.5s
解方程 182s 0.03s 0.54s 1.31s
总耗 时 182s 30.53s 33.44s 35.81s
图 3 MSED 全域划分方法
以单个振子臂为例，划分全域范围成 5 块（图
以传统 RWG 矩量法为阻抗计算值比较基准， 图 5 说明，SED 方法在计算端口输入阻抗上存在一 定误差，尤其是当波束扫描时，这种误差将被放 大。采用本节 MSED 方法对各振子上的全域范围进 行细化，当细化分块数并不大时，输入阻抗的计算 精度即得到了明显的改善。 计算中，在增加单个振子上全域分块数时，亦 增大了 g n 基函数个数，这也增加了 SED 矩量法方 程的求解时间。极端情况即当分块数选为单个振子 上 RWG 基函数个数时，MSED 即退化为传统的矩 量法。因此，相对于 RWG 以及 SED 方法，MSED 是对计算精度与矩阵方程求解难度之间的权衡。图
3 左振子）与 7 块（图 3 右振子） ，分别记为 MSED(5)与 MSED(7)。当原单个振子上 RWG 基函
数数目为 M 时，分块后的 MSED 基函数个数分别 降为 M/7 与 M/5。
3
仿真算法比较
为检验 MSED 方法的合理性，下面首先建立如
图 2 所示的 225 单元周期天线阵列，其中天线单元 选用如图 4 中振子结构：两振子交叉对称放置，振 子臂长均为 L = 0.5m ，振子横截面等效直径 5mm 。两振子中心采用 δ (t ) 电压源馈电，幅度分 别取 1 V 和-1 V ，频率 f = 300 MHz ；单元间距
5 说明，当分块数达到 7 后，采用 MSED 得到的输
入阻抗已经较好地逼近于 RWG 精确解，而矩阵未 知数降低到原来的 1/ 5 。具体时耗分析见表 1 。 ·583·
0.51λ 。用三角面片对其进行剖分，得各交叉振子
（注：计算数据均选用 DELL-T7500 服务器计算得 到，机器配置为 CPU 3.07G*12，RAM 96G） 为进一步检验 MSED 方法计算更大规模阵列天 线时的准确性与有效性，分别建立了 361 元、729 元 周 期 阵 列 ， 阵 列 总 未 知 数 分 别 为 25270 和
计算速度上有较大地提高，同时，对计算机内存要 求也大大降低。因此，当传统矩量法达到计算机速 度与容量限制时，依然可以采用 MSED 法计算更大 的阵列。
表 2 361 阵列计算时耗分析
方法 RWG MSED(5) MSED(7) 未知 数 25270 3610 5054 计算 小阵 无 21s 21s 建立全域 基矩阵 无 30.1s 34.5s 解方程 7 min 1.8s 4.2s 总耗 时 7 min 52.9s 59.7s
51030 。分别比较了不同方法计算阵列天线输入阻
抗，如图 5、6 所示。
表 3 729 阵列计算时耗分析
方法 RWG MSED(5) MSED(7) 未知数 51030 7290 10206 计算 小阵 无 21s 21s 建立全域 基矩阵 无 115.3s 121.7s 解方程 48 min 11.1s 27.8s 总耗时 48 min 147.4s 170.5s
图 5 361 元阵列中心行振子输入阻抗比较
方法 RWG MSED(5) MSED(7)
表 4 阵列计算内存要求
361 元阵列 18.3G 0.39G 0.77G 729 元阵列 78.2G 16.0G 3.1G
4
结论
图 6 729 元阵列中心行振子输入阻抗比较
表 2、3 分别给出了不同方法计算阵列的时 耗。表 4 给出了不同方法计算阵列的内存要求。相 对于求解 M 元周期阵列天线对应的 RWG 矩量法方 程，MSED(n)方法可以将方程未知数降为 M × n ，
本文分析周期排布的交叉振子阵列输入阻抗 时，为提高计算精度，在 SED 算法思路的基础 上，进一步细化了单个振子上的子域全域区域，提 出多尺度子域全域基 MSED 方法。计算结果与矩量 法精确求解数据具有较好吻合性。通过大量数值实 验表明，这种方法在计算周期阵列交叉振子天线中 具有快速、精确、节省内存空间的优点，为分析大 规模周期排布阵列天线提供了有效途径。
·584·
参考文献
[1] [2] [3] [4] [5] [6] [7] E. Suter and J. Mosig, A subdomain multilevel approach for the MoM analysis of large planar antennas, Microwave Opt. Technol. Lett., 2000. Junho Yeo and Mittra, R., Numerically efficient analysis of microstrip antennas using the characteristic basis function method (CBFM)，IEEE Trans. Antennas Propagat., 2003. T. J. Cui, W. B. Lu, Z. G. Qian, and W. Hong, Accurate analysis of large-scale periodic structures using an efficient sub-entiredomain basis function method, in Proc. IEEE Int. Symp. Antennas Propagation, 2004 H. Li, B. Z. Wang, and P. Du, An Extended Sub-Entire Domain Basis Function Method for Finite Periodic Structures, IEEE Trans. Antennas and Wireless Propagat, 2008 H. Li, B. Z. Wang, and P. Du, An Dimension reduction-based sub-entire domain basis function method for finite periodic structures, IEEE Microwave Conference, 2009 H. Zhao, G. H. Li, W. B. Lu, J. F. Zhang and T. J. Cui, An Efficient Algorithm for Accurate Simulation of Split-Ring Resonators Using a Modified Sub-Entire Domain Basis Function Method, IEEE Microwave Conference Proceedings, 2005 T. J. Cui, W. B. Lu, J. L. Cai, X. B. Wang and B. Zhao, Analysis of periodic array structures in half space with infinite electric plane using sub-entire-domain basis function method, IEEE International Workshop on Metamaterials，2008
作者简介： 陈涛，男，助理工程师，主要研究领域为电磁数值计算；张志衡，男，研究员，主要研究领域为电磁 数值计算、反射面天线等。
·585·
微纳尺寸结构的精确电磁模拟
孟 敏 聂在平
(电子科技大学电子工程学院，成都 611731) mengmin@https://www.doczj.com/doc/925597785.html,
摘 要：本文研究微纳尺寸结构的精确电磁模拟。通过分离电荷和电流建立的增广型电场积分方程，可以 克服传统电场积分方程的低频崩溃问题。但我们通过研究增广电场积分方程的矩阵性态，发现其在微纳尺 寸模拟时的不稳定性。本文提出一种基于几何尺度归一化的改进型增广电场积分方程。数值结果验证了本 文方法的正确性和通用性。 关键词：电场积分方程，低频崩溃，增广电场积分方程，归一化
Accurate Electromagnetic Modeling of Micro-Scale or Nano-Scale Structures MENG Min, NIE Zaiping
(School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu 611731) Abstract: We study the accurate electromagnetic modeling of micro-scale or nano-scale structures in this paper. The augmented electric field integral equation (A-EFIE) obtained by separating the charges and currents can overcome the low-frequency breakdown problem of the traditional electric field integral equation (EFIE). However, after carefully studying the matrix property of the A-EFIE, we found that the A-EFIE is not stable when modeling the micro-sale or nano-scale structures. We propose an improved A-EFIE (IA-EFIE) by introducing a normalization factor based on geometrical characterizations. Numerical results are presented to demonstrate the accuracy and generality of this method. Keywords: EFIE; low-frequency breakdown; A-EFIE; normalization or loop-star）[3, 4]可通过准亥姆霍兹分解，将电流源 分解成无散度源和有散度源，分离矢量位和标量位 贡献，并通过频率归一化使二者的贡献均被有效计 入。该方法可有效克服低频崩溃问题。然而，对于 复杂互联结构，全局环基函数的搜索将极为费时。 近年来，有学者提出增广型电场积分方程方法（AEFIE）[5, 6]，作为一种更为简单的方案，一定程度 上替代环树基分解方法。该方法的思路是分离电流 和电荷，使电荷也成为独立的未知量进入计算，通 过追加电流连续性方程，使得矩阵系统有唯一解。 由于可以避免复杂的环基函数搜索， A-EFIE 在实 现上更为简单，计算上也更为高效。 然而，通过矩阵性态分析发现， A-EFIE 方法
1
引言
随着半导体工艺和技术的进步，微纳结构的精 确电磁模拟变得极为迫切。表面积分方程方法是常 用的一种分析方法。然而，由于这类结构通常工作 在 GHz 频段，使得其电尺寸（对比波长）极小。 此时，积分方程存在低频崩溃问题[1]。以电场积分 方程（EFIE）为例，由于奇异项（矢量位）和超奇 异项（标量位）的同时存在，用传统 RWG 基函数[2] 离散该方程时，奇异项部分由于数值误差被噪声淹 没，从而使矩阵病态。从算子方程看，忽略奇异项 使得算子存在零空间。环树基分解方法（ loop-tree ·586·
在处理工作在 GHz 频段下的微纳结构时，存在不 稳定性问题。问题的原因在于 A-EFIE 矩阵的构建 依赖于实际工作频率，而非目标本身的电尺寸。本 文通过提出基于几何尺寸的归一化方法，修正原 AEFIE 系统，使其能稳定描述工作在任意频率下的 电小尺寸问题，如工作在 GHz 频段下的微纳尺寸 结构。
一个鞍点系统，数学上可以被精确求解[7]。
3
改进型 A-EFIE
上述 A-EFIE 系统成为鞍点系统可被精确求解 的前提假设是频率极低，即 k0 → 0 。对于工作在
2
增广电场积分方程（A-EFIE）
GHz 频段的微纳尺寸结构，该假设显然无法满足。 为此，我们提出一种基于几何尺度归一化的改进方 案，表达式如下[8]：
考虑理想导电体（PEC）受到电磁源的激励， 表面将产生感应电流。该电流可以通过下述电场积 分方程（EFIE）求解：
?A ? Lc ? ? FiD
? DT iPi Lc iB ? ?ik0 Lc J ? ?η0 1V ? ? i ?c ρ ? = ? ? (9) ?? 0 r ? ? 0 ? 2 2 k0 Lc Ir ? ? Lc ?
其中， Lc 是一个待求结构的特性长度因子。 显然，该系统对工作在任意频段的电小尺寸问题均 是一个鞍点系统，各矩阵块完全由结构的电尺寸描 述，而不依赖于特定的工作频率。
(1)
其中， g (r, r ′) 为标量格林函数。通过 RWG 基 函数离散电流，并利用伽略金方法，可以将上述积 分方程离散为如下一个矩阵方程：
4
数值算例
Z iI = V
阻抗矩阵 Z 包含矢量位和标量位贡献
(2)
Z = Z A + Z S = ik0η0 A +
其中，
η0
ik0
S
(3)
本节给出数值算例。由于球体散射具有解析解 可供对比，我们计算了一个半径为 1μm 的金属球在 f=20GHz 平面波照射下的散射结果，如图 1 所示：
[ A] ji = ∫ drf j (r )i ∫ dr ' g (r,r ') fi (r ') (4) S S'
[S ] ji = ∫ dr?i f j (r ) ∫ dr ' g (r,r ')? 'i fi (r ') (5)
S S'
若分离电荷作为独立的未知量，则标量位部分 贡献可以采用定义在单个三角形上的脉冲基函数 p(r) 展开及测试：
[P] ji = ∫ drp j (r) ∫ dr ' g (r,r ') pi (r ') (6) S S'
通过分析，上述 S 和 P 存在一个简单的关系：
图 1 半径 1μm 金属球的双站 RCS（f=20GHz）
(7) 其中 D 是一个联系 RWG 基函数和脉冲基函数定义 域的稀疏映 射矩阵。因 此， A-EFIE 可以如下 建 立：
? ?A DT iP iB ? ?ik0 J ? ?η0 1V ? i? =? (8) ? ? ? ? 2 ? ? ? FiD k0 Ir ? ?c0 ρ r ? ? 0 ? 其中的 B 和 F 是两个变换矩阵，引入后可以去
S = D iPiD
T
从图中可以发现， EFIE 由于低频崩溃问题， 结果错误，而 A-EFIE 由于描述该类问题时的不稳 定性，单精度结果也不正确，需要采用双精度计算 方可得到与解析解吻合的结果。但本文提出的修正 型 A-EFIE 方法，即使采用单精度计算，也能得到 理想的结果。
掉该系统的线性相关性[6]。在极低频时，该系统是 ·587·
5
结论
本文通过研究增广电场积分方程的矩阵性态， 发现其在微纳尺寸模拟时的不稳定性，并提出一种
基于几何尺度归一化的改进型增广电场积分方程。 该方法适用于工作在任意频段的电小尺寸问题的分 析。数值结果验证了本文方法的正确性和通用性。
参考文献
[1] [2] [3] [4] [5] [6] [7] [8] D. R. Wilton and A. W. Glisson, “On improving the electric field integral equation at low frequencies,” in URSI Radio Science
Meeting Dig., Los Angeles, CA, Jun. 1981, pp. 24–24
S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surface of arbitrary shape,” IEEE Trans. Antennas
Propag., vol. 30, pp. 409–418, May 1982
G. Vecchi, “Loop-star decomposition of basis functions in the discretization of the EFIE,” IEEE Trans. Antennas Propag., vol. 47, pp. 339-346, Feb. 1999 J.S. Zhao and W. C. Chew, “Integral equation solution of Maxwell’s equations from zero frequency to microwave frequencies,”
IEEE Trans. Antennas Propagat. vol. 48, pp. 1635-1645, Oct. 2000
Z. G. Qian and W. C. Chew, “An augmented EFIE for high speed interconnect analysis,” Microw. Opt. Tech. Lett., vol. 50, no. 10, pp. 2658–2662, Oct. 2008 Z. G. Qian and W. C. Chew, “Fast full-wave surface integral equation solver for multiscale structure modeling,” IEEE Trans.
Antennas Propag., vol. 57, pp. 3594–3601, Nov. 2009
M. Benzi, G. H. Golub, and J. Liesen, “Numerical solution of saddle point problems,” Acta Numer., vol. 14, pp. 1–137, Apr. 2005 M. Meng, and Z. P. Nie, “Improved augmented electric field integral equation method for modeling micro-scale structures,” Electromagnetics, vol. 33, pp. 144-152, 2013
作者简介： 孟敏，女，讲师，主要研究领域为计算电磁学，电磁散射与辐射，电磁理论及应用等；聂在平，男， 教授、博士生导师，主要研究领域为电磁散射与逆散射，计算电磁学，天线工程等。
·588·
衍射高斯波束分析方法的三维化建模方法
焦天栋 1 陆泽健 1 刘小明 1 俞俊生 1 陈晓东 2 姚 远 1
(北京邮电大学电子工程学院，北京 100876)1 (伦敦大学玛丽女王学院电子工程与计算机科学学院，伦敦)2 dongdong_1800@https://www.doczj.com/doc/925597785.html,
摘 要：衍射高斯波束分析方法由于采用几何光学计算极大提高了计算速度，同时利用衍射来补充结果保 证了计算精度。但是传统 DGBA 算法只能对二维对称结构进行建模，这是相较 GRASP 分析的一大缺陷。 本文提出一种新的 DGBA 建模方法，使其能够实现三维衍射高斯波束的分析。文章末尾验证了这种建模方 法的精度。 关键词：衍射高斯波束分析方法，三维化，建模
Diffracted Gaussian Beam Analysis for Three-dimensional Structure Modeling Jiao Tiandong1, Lu Zejian1,Liu Xiaoming1, Yu Junsheng1, Chen Xiaodong2, Yao Yuan1
(School of Electronic Engineering, Beijing University of Posts and Telecommunications, Beijing 100876)1
(School of Electronic Engineering and Computer Science, Queen Mary University of London, London)2
Abstract: Diffracted Gaussian beam analysis(DGBA) takes the advantage of computational efficiency of geometric optics, and introduces diffraction to improve the calculation accuracy. But previous DGBA is only capable of modeling two-dimensional symmetrical structured, which is a major flaw in contrast with the commercial software GRASP. This article put forward a new method of modeling so that it can be utilized for 3D diffracted Gaussian beam analysis. The precision of the analytical methods has been verified in this paper. Keywords: Diffracted Gaussian beam analysis; Three-dimensional; modeling
1
引言
度可保证的优点，但是缺点是不适用于三维结构。 本文介绍了它的三维化建模方法。
毫米波太赫兹频段介于微波与红外波段之间， 是电磁波频谱中很有研究价值的一个频段，也是比 较新兴的研究领域。准光技术是毫米及亚毫米波系 统分析与设计的一个重要手段。准光器件具有加工 难度相对低，承受功率高，可多波束多极化工作， 损耗小等优点，在毫米波亚毫米波的传输分析设计 中得到广泛的应用。衍射高斯波束分析方法作为分 析准光系统的重要方法，具有计算速度快、计算精
基金项目：本研究由十二五民用航天项目和中央高校基本 科研业务费专项资金资助。
2
衍射高斯波束分析方法的介绍
基于几何光学的衍射高斯波束分析方法 （ DGAB ） 首 先 由 伦 敦 大 学 玛 丽 女 王 学 院 的 C.Riechmann 博士等人一起研发的旨在快速分析电 大尺寸反射镜系统的电磁场分布的一种有效的方法 及实现。其主要分析方法如图 1-A： 对入射场（可以是来自馈源的近场或者前一个 镜子的辐射场）在入射平面上进行过采样，然后将 其分解为一系列的三维的高斯波束（如图 1-B） ，并 ·589·
用相应的高斯波束系数表示。每一个高斯波束都传 向待分析的镜子。 对每一束传播到镜面且能量大于阈值的波束进 行分析。对于投射到镜面边缘的所有高斯波束，要 根据相关理论分别对衍射和反射进行计算。对于投 射到镜面中心附近的波束只考虑反射的影响。反射 只需要根据简单的反射定律根据镜面的曲率半径将 改变波束的复参数，而衍射则利用高斯波束入射到 基尔霍夫半平面的正则问题的渐进解来计算。
反射镜的入射波束和反射波束也均在此平面内。对 于不满足此条件的准光系统，原始的 DGBA 无法对 其进行分析。这就大大限制了 DGBA 的应用。
图 2 抛物面反射镜设计示意图
be ana ly
}
di ffr act ed be am s re fle cte d be am s
sed
o u t put p la ne
GB
}
ex
di ffr act ed be am s
in pu t p lan e
pre v i ou s ref lec to r
e x la n p L0
n pa
si
on e
新的建模方法则突破了这个限制，对于多反射 镜系统，反射镜不必沿着光线传播的方向摆放，而 是可以在三维空间中自由摆放。对于每个反射镜， 我们都会建立一个坐标系，各个坐标系之间通过坐 标变换来进行转换。 当我们假设反射镜足够大的时候，我们可以暂 时忽略镜面的边缘衍射，这时候高斯波束入射到镜 面之后以反射为主。曲面上高斯波束反射示意图如 图-4 所示：
sr ξ
r
GB expansion
ref lec
tor to
pan
si o
n
ηr θ θ
n ξi si ηi
η ξ
图 1 DGBA 分析一个反射镜的模块化过程以及入射 波在入射平面内展开成高斯模的示意图
图 3 曲面上的高斯波束反射
最后，出射平面上的场分布由反射部分和衍射 部分叠加而成。如果是最后一个反射镜，近场就由 此得到。如果不是最后一个反射镜所得到的场将作 为下一个反射镜的入射场。这使得衍射高斯波束分 析方法能够做到高度的模块化。
其中反射光线的曲率半径 R1 和 R2 可以通过 如下公式求解：
Σr ? ? ? cos( 2 α r ) Δ r ? sin( 2 α ) Δ r r ? Σi ? ? ? ? ? ? ? = ? cos( 2 α i ) Δ i ? ? ? ? sin( 2 α ) Δ ? i i ? ? ? ?1 / cos θ + cos θ ? 1 / cos θ ? cos θ + ? ? 0 ?
r
r
1 / cos θ ? cos θ 1 / cos θ + cos θ 0
0?? Σc ? 0 ? ? cos( 2 α c ) Δ c ? 2 ? ? sin( 2 α c ) Δ c ??
? ? ? ? ?
(1)
3
三维 DGBA 的建模方法介绍
i 式 (1) 中 ， Σ i ,c,r = (1 / R 2,c ,r + 1 / R1i ,c ,r ) / 2 ，
i Δ i ,c , r = (1 / R 2,c , r ? 1 / R 1i ,c , r ) / 2 。而镜面的两个主曲
如图-3 所示，是原始 DGBA 中抛物面反射镜 的设计示意图。从图中可以看出，反射镜是用半径 为 R 的光柱来截取相应的二次曲面得到的。反射镜 以 XoZ 平面为对称面，馈源中心以及入射光线和反 射光线均在此平面内，对于多反射镜的准光系统， 也必须满足所有的反射镜关于 XoZ 平面对称，各个 ·590·
c c 率半径 R1 , R2 和旋转角度 α c 的求解是一个难点。
在此，我们利用二次曲面的解析方程根据微分几何 的知识来进行求解。 设高斯波束与标准反射面的交点为
P0 ( x0 , y0 , z0 ) ,则点 P0 ( x0 , y0 , z0 ) 的主曲率半径
C R1,2
为如下方程的两个根 (rt ? s2 )R2 + h[2 pqs ? (1+ p2 )t ? (1+ q2 )r]R + h4 = 0
2 2 2
1 0.9 0.8
展展展展展展展展展展展展展
展展展 展展展
归归归归归 (d B )
(2) 其中， r =? z / ?x |x=x0,y=y0 ， s =? z / ?x?y |x=x0 , y=y0 ，
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -80 -60 -40 -20 0 20 40 60 80
t = ?2 z / ?y 2 |x= x0 , y = y0 ， p = ?z / ?x |x=x0 , y= y0 ，
q = ?z / ?y |x = x0 , y = y0 和 h = 1 + p + q 。
2 2
以入射点
P0
处反射面的切线建立的坐标系
1
r(lamda)
? ? ξ ,η ?i （ c c ）是由入射高斯波束的方向 s 与入射点 P
? 处反射面的法线 n 所确定的，即
? ? ? ξc = n × si
入射点 P 0 切平面
d y d x
传传传传D后展展展展展展展展展 展展展 展展展
0
0.9 0.8 0.7
归归归归归 (dB)
(3) 处反射面的主法截线 C1 和 C 对应的主
2
? ? η c = n × ξ?c
0.6 0.5 0.4 0.3 0.2 0.1 0 -80 -60 -40 -20 0 r(lamda) 20 40 60 80
为如下方程的两个根
dy dy (GM ? FN )( )2 + (GL ? EN ) + (FL ? EM ) = 0 dx dx (4)
图 4 展开平面和传播距离 D 之后重构场和原始场分布对比
展展后入散展散展 1 0.9 0.8 0.7 0.6 0.5
其中 E、F 和 G 为曲面的第一基本量，L、M 和 N 为曲面的第二基本量。 根据主切平面易得到主法截线 C 1 和 C 的切线方 ?c 向 u? 和 v ，即两个主曲率方向，那么旋转角度 α c 为
2
c
? α c = ± cos (ξ?c ? u c )
?1
? 其中，当 ξ c
与 u?
c
? 关于法线 n
(5) 为逆时针转动时
0.4 0.3 0.2 0.1 0 -0.5 0 0.5
取+，反之取-。
4
仿真结果分析 5 结论
图 5 展开后反射场和入射场对比
下面是对三维 DGBA 算法仿真结果的验证。图 4—图 5 是高斯波束斜入射一个抛物面反射镜的相 。这个系统的频率是 关仿真结果（入射角为 30 度） 30G Hz，反射镜的直径是 1m，束腰半径是 5λ， 采样间距是 50λ。 图 4 是展开平面和传播了距离 D=0.5m 之后重 构场和原始场分布的对比图，从图中可以看出，展 开平面和传播了距离 D=0.5m 之后的重构场和原始 场非常吻合。图 6 是展开后反射场和入射场分布的 对比图，展开平面的反射场和入射场分布也非常吻 合，验证了算法的正确性。
三维 DGBA 建模方法解决了原始 DGBA 算法 无法处理三维结构的问题，同时在计算时间相比物 理光学法有优势的情况下精度也得到了保证。由于 本文中假设反射镜足够大，所以没有考虑边界衍射 的情况，所以在今后工作中，可以加入对边界衍射 处理的算法，以进一步提高算法的精度。
参考文献
[1] [2] [3] [4] C. Rieckmann, M. R. Rayner, C. G. Parini, D. H.Martin and R.S.Donnan, "Novel modular approach based on Gaussian beam diffraction for analyzing quasi-optical multi-reflector antennas", lEE Proc.-Microw. Antennas Propag., 149, No. 3, June 2002. DAUBECHIES, I. “The wavelet transform, time-frequency localization and signal analysis”, IEEE Trans. Inf. Theory, 1990, 36, (5), pp. 961-1005 TAKENAKA, T., and FUKUMITSU, O. “Asymptotic representation of the boundary-diffraction wave for a three-dimensional Gaussian beam incident upon a Kirchhoff half-screen”, J. Opt. Soc. Am., 1982, 72, (3), pp. 331-336 P. F. Goldsmith, "Quasi-optical Systems: Gaussian Beam Quasi-optical Propagation and Applications," IEEE Press, 1998
·591·
基于双线性变换的扩展 ADI-FDTD 算法
夏 芬 褚庆昕
(华南理工大学电子与信息学院，广州 510640) xia.fen@https://www.doczj.com/doc/925597785.html,, qxchu@https://www.doczj.com/doc/925597785.html,
摘 要：本文提出了分裂时间步下的双线性变换技术。这种新的双线性变换可通过 Z 变换技术，用于加载 集总网络的扩展无条件稳定算法的研究，为 FDTD 算法分析复杂的集总网络提供了另一个指导思想。为了 验证该双线性变换技术能用于扩展无条件文件稳定算法的可行性，文章基于该变换分析了扩展 ADI-FDTD 算法，并用这种分析方法计算了具体的实例。从算法的计算结果可知，分裂时间步下的双线性变换技术可 用于扩展 ADI-FDTD 算法中。 关键词：分裂时间步，双线性变换，Z 变换技术，扩展 ADI-FDTD 算法
Extended ADI-FDTD Method Based on Bilinear Transformation XIA FEN , CHU QING-XIN
(School of Electronic and Information Engineering, South China University of Technology, Guangzhou 510640) Abstract: Bilinear transformation based on the time split scheme is proposed in this paper. Through Z-tranform, this new bilinear transformation can be used to research the extension of the unconditionally stable algorithm, providing a guiding ideology for the finite difference time domain (FDTD) algorithm to analyze complex lumped network. For the sake of showing the validity of the proposed scheme, a specific example is computed and it shows the new bilinear transformation can be used to the extended ADI-FDTD method. Keywords: Time split scheme; Bilinear transformation; Z-transform; Extended ADI-FDTD method
1
引言
时域有限差分(FDTD)算法可以直观地反映电磁 场的时间变化过程，且计算程序容易编写，被广泛 用于分析电磁问题的各个方面[1]。用 FDTD 算法分 析微波有源电路也是众多学者研究的一个热点[2]， 若加载器件是比较复杂的集总网络时，通常采用 PLRC 分析方法[3]或者 Z 变换分析法[4]。若一个完 整的时间步分裂成两个子时间步时，子时间步下的
PLRC 技术已做推导[5]，但是 Z 变换技术的分析推 导至今尚未见报。若给定集总负载模型，可得到模 型在 s 域的导纳形式，这种导纳形式要进行 Z 变 换，首先要将导纳从 s 域变换到 z 域，通常采用双 线性变换[4]。 本文基于完整时间步下的双线性变换，提出了 分裂时间步的双线性变换技术，并将这种技术用于 分析扩展 ADI-FDTD 算法。
2
基金项目：国家自然科学基金资助项目 (No. 61171029) ，
分裂时间步下的双线性变换
(No. 61201036) 和 广 东 省 自 然 科 学 基 金 资 助 项 目
（No.S2012040006893） 。
要分析 Z 变换技术，先介绍分裂时间步下的双 线性变换。线性系统的系统函数为：
·592·
H (s) =
b s+a
(1)
这个系统也可以写成微分方程的形式： dy ( t ) + ay ( t ) = bx ( t ) dt 方程的解用导数的近似积分形式来表示：
y ( t ) = ∫ y ' (τ )dτ + y ( t0 )
t0 t
2 ? ( Δt ) B[i + 1] = ? 1 + ? 2εμ ( Δx )2 ?
? ? ? ?
(11)
(2)
1 ? Δt ? n ? 1 1? 1 1 ?? ? n? X [i + 1] = Ezn ? i, j, k + ? + H y ? i + , j, k + ? ? H y ? i ? , j, k + ? ? 2 ? 2εΔx ? ? 2 2? 2 ?? ? ? 2 ? ? 1 1? 1 1 ?? Δt ? n ? ? H x ? i, j + , k + ? ? H xn ? i, j ? , k + ? ? 2εΔy ? ? 2 2? 2 2 ?? ? ?
(3)
? n? 1 1 ? ? ? n? ? E x ? i + 2 , j , k + 1? ? Ex ? i + 2 , j , k ? ? ? ? ? ? ? ? ? 4εμΔxΔz ? n ? 1 1 ? ?? n? ? Ex ? i ? , j , k + 1? + E x ? i ? , j , k ? ? ? ? 2 ? ? 2 ?? ?
(12)
( Δt )
2
当 t = nT 且 t0 = nT-T/2 时, 公式(3)近似为：
y (t ) = T? T ?? T? ? ? y ' ( nT ) + y ' ? nT ? ? ? + y ? nT ? ? (4) 4? 2 ?? 2? ? ? ?
加载集总器件后，扩展算法的表达式变为：
X [i + 1]? = Δt n + 1 ? 1? J z 4 ? i, j , k + ? 2ε 2? ?
(13)
当 t = nT 时，微分方程(2)可以写为：
y ' ( nT ) = ? ay ( nT ) + bx ( nT )
(5) (6)
将式子(5)代入(4)中，并应用 Z 变换，可得：
1 1 ? ? bT ? ? aT ? ? aT ? ? 2 2 ? X ( z ) + z X ( z )? ?1 + ? Y ( z ) ? ?1 ? ? z Y (z) = 4 ? 4 ? 4 ? ? ? ?
一旦给定集总网络，其导纳 Y(s)可得到，再通 过双线性变换公式(8)，可得到 Z 域下的导纳 Y(z)。 Y(z)与 E(z)、J(z)的关系为：
E (z) J (z) = ΔzY ( z ) ΔxΔy
(14)
在 Z 域，系统函数可以表示为：
H (z) = Y (z) =
1 ? ? bT ? ?1 + z 2 ? 4 ? ?
假定 Y(z)的表达式为：
b 4 1? z +a 1 ? T 1+ z 2
1 ? 2
X ( z)
? aT ? ? aT ? ?1 + ? ? ?1 ? ?z 4 ? ? 4 ? ?
1 ? 2
=
(7)
Y (z) =
ΔxΔy M ∑ ck z 2 Δz k =0
k
1+∑ d k z
k =0
M
-
k 2
(15)
将(15)代入到(14)中， J(z) 与 E(z)的关系式：
Jz
n+
对比公式(1)和(7)，可得到在分裂时间步下 s 域 到 z 域的双线性变换为：
4 1? z 2 s= 1 ? T 1+ z 2
?
1
1 2
+∑ d k J z
k =0
M
1 k n+ 2 2
= c0 Ez
n+
1 2
+∑ ck Ez
k =0
M
1 k n+ 2 2
(16)
(8)
公式(8)为分裂时间步下的双线性变换公式，该 方法可以为扩展 FDTD 算法分析复杂的集总网络提 供了另一个指导思想。
上式可以等效为下面的迭代形式: 1 n+ ? n+ 1 2 2 ? J z = c0 Ez +w1 1 k 1 k ? n+ n+ ? ? wk = wk +1 +ck Ez 2 2 ? d k J z 2 2 ? 1 M 1 M n+ n+ ? wM = cM Ez 2 2 ? d M J z 2 2 ? ?
(17)
3
扩展 ADI-FDTD 算法
假设分析区域是理想媒质，集总元件加载在 Z 方向上，传统 ADI-FDTD 算法的 Ez 在第一子时间 步下的迭代公式可以写为： n→n+1/2：
A[i +1] ? Ez
n+ 1 2 1 n+ ? 1? 1? ? 2 ? i ?1, j, k + ? + B[i +1] ? Ez ? i, j, k + ? 2? 2? ? ?
将等式 (17) 代入等式 (16) 中，扩展 ADI-FDTD 算法变为： Δt ? c0 B[i + 1]+ = (18) 4ε
X [i + 1]? = ? 1? Δt ? n + 1 ? ? J z 2 ? i, j , k + ? + w1 ? 4ε ? 2? ? ?
(19)
(9)
等式(18)和等式(19)构成了扩展 ADI-FDTD 算 法的关键部分。
1 n+ ? 1? +C[i +1]? Ez 2 ? i +1, j, k + ? = X [i +1] 2? ?
其中：
4
(10)
数值分析
( Δt ) A[i + 1] = C[i + 1] = ? 2 4εμ ( Δx )
2
用包含集总网络的微带线结构来验证上文提出 的设计方法，如图 1 所示。计算区域为 12 × 18 × ·593·
2mm, 网格尺寸为 Δx = Δy = 0.2 mm, Δz = 0.1 mm。 介质板厚度为 3Δz，介电常数为 εr = 2.2，微带尺寸 为 6Δx × 80Δy。导带一端接阻值为 100?、振幅为 1V、频率为 20GHz 的电阻性正弦电压源，另外一 端(30Δx, 80Δy)处，导带与地板之间接集总网络，该 集总网络结构如图 2 所示。
1
Z算算算算算算
0.5 电电（V）
0
-0.5
－Z-FDTD -Z-ADI-FDTD 200 400 时时时 600 800
-1 0
ε = 2.2
图 3 集总负载两端的电压图
图 3 给出了扩展 ADI-FDTD 算法和扩展 FDTD 算法的计算结果。从该图中可以看出，两者的计算 结果基本吻合。进而说明了分析方法的可行性。
图 1 微带线端接集总网络结构图
5
结论
图 2 集总网络电路图
本文提出了分裂时间步下的双线性变换技术。 一旦给定任意线性的集总网络，就可以得到导纳 Y(s)，通过分裂时间步下的双线性变换技术可得到 z 域的导纳 Y(z)，导纳 Y(z)可以采用 Z 变换技术得 到电流和电场之间的关系。数值计算结果表明：基 于双线性变换的扩展 ADI-FDTD 算法的计算结果与 扩展传统的 FDTD 算法吻合，说明了上述的分析方 法的可行性。
参考文献
[1] [2] [3] [4] [5] Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 2nd ed. Boston, MA: Artech House, 2000 J. A. Pereda, F. Alimenti, P. Mezzanotte, L. Roselli, and R. Sorrentino, “A new algorithm for the incorporation of arbitrary linear lumped networks into FDTD simulators,” IEEE Trans. Microw. Theory Tech., vol. 47, no. 6, pp. 943-949, Jun. 1999 J. Y. Lee, J. H. Lee, and H. K. Jung, “Linear lumped loads in the FDTD method using piecewise linear recursive convolution method,” IEEE Microwave Wireless Compon. Lett., vol. 16, no. 4, pp. 158-160, Apr. 2006 H. E. Abd EI-Raouf, W. Yu, and R. Mittra, “Application of the Z-transform technique to modeling linear lumped loads in the FDTD,” Proc. Inst. Elect. Eng.-Microw. Antennas Propag., vol. 151, no. 1, pp. 67-70, Feb. 2004 Fen Xia, Qing-Xin Chu, Yong-Dan Kong, and Zhi-Yong Kang, “The ADI-FDTD Method Including Lumped Networks Using Piecewise Linear Recursive Convolution Technique,” Progress in Electromagnetics Research M, vol. 30, pp. 66-77, 2013
作者简介： 夏芬，女，硕士，主要研究领域为计算计算电磁学；褚庆昕，男，教授、博士生导师，研究方向包括 射频与微波电路、有源集成天线和时域计算电磁学等。
·594·
陷波超宽带天线同极化电场分量的时域信号分析
杨 光 1 褚庆昕 2
(广州慧睿思通信息科技有限公司，广州 511442)1 (华南理工大学电子与信息学院，广州 510641)2 yang.guang@https://www.doczj.com/doc/925597785.html,
摘 要：超宽带（UWB）天线在引入陷波结构之后，会对信号传输造成一定程度的失真，导致振铃现象产 生。本文从 UWB 天线使用不同技术引入陷波结构的角度出发，研究和分析不同陷波技术对 UWB 信号瞬时 响应的同极化方向电场分量的影响，结果表明陷波 UWB 天线仍然能够很好地进行 UWB 脉冲信号传输。 关键词：超宽带天线（UWB） ，时域分析，同极化电场分量
Time Domain Analysis of Co-polarized Component of the Transient Response on Band-Notched UWB Antennas YANG Guang1, CHU Qing Xin2
(Guangzhou HuiRuiSiTong Information Technology LTD, Guangzhou 511442)1 (School of Electronic and Information Engineering, South China University of Technology, Guangzhou 510641)2 Abstract: After the introduction of band-notched structures into the Ultra-wide-band (UWB) antennas, they cause certain degrees of distortion of the signal transmission and lead to the ringing effects. In this paper, from the views of the different band-notched structures, the influences on time domain analysis of co-polarized component of the transient response have been researched and analyzed. The results show that the band-notched UWB antennas are still good candidates for UWB pulse signal transmission. Keywords: Ultra-wide-band (UWB), Time Domain Analysis, Co-polarized Component
1
引言
天线传统的电性能参数可以通过跟频率有关的 频域参数描述，例如效率、输入阻抗、增益和辐射 方向图[1]。对于 UWB 系统，UWB 天线占用非常宽 的频带，天线参数是用频率的函数来实现的，仅用 这些电性能参数并不能完整的描述 UWB 系统。采 用传统频域参数来评估天线性能具有局限性，因而 需要对 UWB 天线展开时域特性的评估[2]。 在 UWB 天线中引入陷波结构，可以让天线在 频域中具有良好的陷波特性，能够有效地抑制窄带
基金项目：国家自然科学基金资助项目(61171029)
系统与超宽带系统之间的相互干扰。但是，陷波结 构对天线的时域传输特性有何影响，是否会引发超 宽带信号畸变，是值得深入研究的问题。因此，对 于 UWB 天线，除了频域特性外，时域特性也是评 估 UWB 天线的重要标准。
2
超宽带天线的时域瞬时响应
为了评估 UWB 天线的时域性能，除了传递函 数、群延时、相关系数和脉冲宽度拉伸比[3]之外， 需要对天线的瞬时响应进行研究，可以从频域的传 输方程 H(jω)入手。为了简化，只考虑同极化方向 电场[4]，亦即： ·595·
?2 H ( jω ) ? + H同极化 ( jω ) = ? 同极化 ?0 ?
ω >0 ω≤0
(1）
瞬态响应的 H 同极化(t)同极化方向电场 H+同极化(jω) 是传输函数解析的逆傅里叶变换 H+同极化(t)的实部， 表达式是：
+ H同极化 (t) = Re[ H同极化 (t)]
（2）
(c)天线 3
(d)天线 4
+ H同极化 (t,θ ,ψ ) 是天线接收信号瞬态响应的包
图 1 时域分析讨论的四种超宽带天线
络函数。 首先，天线瞬时响应（同极化方向电场）的最 大值 P 定义为： + P (θ ,ψ ) = max H同极化 (t,θ ,ψ ) （3）
t
4
同极化电场方向时域信号分析
其次，天线瞬时响应的持续时间，可以由最大 值一半（Full Width at Half Maximum）的全宽时间 跨度得到，定义为：
四种 UWB 天线同极化电场分量方向的瞬时响 应和其包络曲线如图 2，图 3，图 4 和图 5 所示：
τ FWHM (θ ,ψ ) = t1 ||H
?t2 |t
+ 2 < t1 ,| H同极化 (t 2 ,θ ,ψ
+ 同极化 ( t1 ,θ ,ψ
)|= | p /2
(4）
)|= | p /2
最后，讨论振铃现象的本质，需要研究振铃信号 的最大值和持续时间，定义为包络曲线从最大值滑落 到某特定值（通常定义为最大值的 α 分式） ，如：
τ r ,α (θ ,ψ ) = tα ||H
?t p |t
+ 同极化 ( tα
,θ ,ψ )|=α p
图 2 天线 1 瞬时响应和其包络曲线
(5）
+ α < t p ,| H同极化 (t p ,θ ,ψ
)|= p
3
时域分析讨论的陷波超宽带天线
在天线上开槽实现陷波的技术主要是分为两 类，一类是在辐射体上开槽，一类是在馈线端开槽 [5] 。以图 1 中四种超宽带天线为示例来完成超宽带 天线时域分析的讨论，图 1（a）为无陷波超宽带天 线，图 1 中（b） c） d）陷波天线则是在（a）的 （ （ 基础上分别在辐射体开槽、馈线端开槽以及结合这 两种槽实现陷波的超宽带天线[6]。
图 3 天线 2 瞬时响应和其包络曲线
(a)天线 1
(b)天线 2
图 4 天线 3 瞬时响应和其包络曲线
·596·