摘要:天线的方向系数是天线的核心性能指标之一, 准确计算方向系数是高性能天线应用的核心要求. 本文基于平面近场测试理论、实测数据和快速傅里叶变换算法, 系统阐述基于近场测试来数值计算天线方向系数的原理, 并进行深入的误差分析. 本文选择一种方向图函数和方向系数已知的被测天线, 来检验所讨论的误差评估方案. 评估分两步实现, 第一步, 针对这一天线, 采用标准的近场测试配置, 仿真模拟出(相当于实际测量出)一套平面近场数据. 第二步, 基于这套近场数据, 利用数值积分计算出天线方向系数. 本文使用或提出了四种数值算法, 分析了提出的后三种算法本身的误差来源, 并开发出程序搜索方案, 确定出后两种算法的最小误差界. 随后, 利用这四种数值算法分别得出天线的方向系数. 结果表明, 计算所得的近场方向系数都比真实方向系数大, 但误差不超过0.6 dB. 这一结果对实际应用中正确评估基于近场测试的天线方向系数准确性有重要参考价值.
关键词: 天线/
方向系数/
数值积分/
逼近

English Abstract

全文HTML

--> --> -->

天线是实现电磁波辐射与接收的基本设备, 广泛应用于广播、通信、雷达、导航等应用系统中. 一支好的天线, 首先要能够高效率辐射和接收电磁波, 即输入天线的电磁场功率能够大部分被辐射出去或接收回来. 其次, 要根据实际需求, 实现方向性. 简单来说, 方向性是指天线辐射或接收的电磁波功率密度跟空间方向紧密相关: 某些方向电磁波功率密度很大, 而某些方向电磁波功率密度很小直至忽略不计. 典型如跟踪雷达天线, 它的电磁波功率只聚焦于空间很狭窄的方向区域, 以实现辨别和测定特定目标位置. 电磁波天线已经用于国家经济活动的不同层面和千家万户. 这一事实表明, 无论是从应用物理研究角度还是适应实际发展需求角度, 高效天线设计与制造不会是一个过时的课题.
天线辐射的电磁场空间方向分布特性是天线的主要性能, 天线的方向系数是性能核心指标之一. 虽然天线的形式多种多样, 但物理上, 任一天线的辐射性能可通过天线的归一化方向图函数F(θ, ?)来描述. 天线在空间球坐标(θ, ?)方向的方向系数D(θ, ?)定义为该方向辐射强度与整个空间所有方向的平均辐射强度之比. 方向系数因此可理解为天线聚集电磁波能量程度大小的度量. 天线方向系数D(θ, ?)与归一化方向图函数F(θ, ?)的关系如下^[1]:

$D\left( {\theta,\phi } \right) = \frac{{4{\text{π}}{F^2}\left( {\theta,\phi } \right)}}{{\displaystyle\int_0^{2{\text{π}}} {\displaystyle\int_0^{\text{π}} {{F^2}\left( {\theta,\phi } \right)\sin \theta {\rm{d}}\theta {\rm{d}}\phi } } }},$

(1)式是计算任意天线方向系数的基本公式. 工程计算中, 通常所说的天线方向系数是指D(θ, ?)的最大值D₀. 根据(1)式, D₀公式为

${D_0} = \frac{{4{\text{π}}}}{{\displaystyle\int_0^{2{\text{π}}} {\displaystyle\int_0^{\text{π}} {{F^2}\left( {\theta,\phi } \right)\sin \theta {\rm{d}}\theta {\rm{d}}\phi } } }}.$

如上所述, 方向系数D₀是一支天线最重要的性能指标之一, 其准确可靠计算是天线设计与应用需要关注的核心问题之一. 天线另一个重要性能指标“增益”G(θ, ?)与方向系数D(θ, ?)紧密相关, 本文着重于方向系数的计算, 对天线增益不做进一步讨论, 在功率损耗可忽略的情况下, 天线方向系数与天线增益大小相同.
根据(2)式, 计算方向系数必须首先获得天线归一化方向图函数F(θ, ?)的值, 但是除了几种几何形状简单的天线外, 实际上很难获得F(θ, ?)的解析表达式. 工程上, 通过实际天线的平面近场测试数据来计算F(θ, ?)是一种可行方法, 基于天线平面近场测试数据实现F(θ, ?)和D₀的高精度计算成为电磁场计算的一个重要课题.
天线平面近场测试已形成一套比较成熟的理论^[2], 其中探头在天线辐射方向端测试采样的所有点构成一个横向、纵向间距均匀的矩形栅格平面, 该平面与被测天线相距一般为2λ—10λ, 这里λ为被测电磁波波长, 获得数据为探头在测试采样点感应到的复电压或复电流, 该数据与被测天线激发的电磁场成正比. 显然, 平面近场测试获得数据为一个复数矩阵, 矩阵大小由测试选择的横向、纵向采样点数决定.
平面近场测试理论是一个有一定误差的工程理论, 这里列举平面近场测试的主要误差源: 第一, 天线远场与平面波谱的简单正比关系式是平面近场波数谱积分利用驻相法获得的近似计算结果; 第二, 平面近场离散采样带来的高波数谱贡献被省略的误差; 第三, 平面近场有限大小测试平面导致的截断误差; 第四, 一般近场测试数据处理未考虑探头方向性带来的误差. 当然, 也可从电磁场等效原理来分析平面近场测试技术的误差源. 另外, 平面近场测试只考虑被测天线靠近探头侧即主要辐射方向半空间的电磁场分布, 实际天线一般满足此条件, 背瓣方向电磁场相对较弱, 在计算时一般不予考虑, 认为天线偏离探头侧的背面半空间电磁场为零.
如上所述, 尽管平面近场测试理论和技术有一定误差, 但是实际工作中仍然被广泛使用, 实践证明平面近场测试误差在根据被测天线特点选择对应的测试条件后可以容忍直至忽略不计. 因此, 基于天线平面近场测试数据实现F(θ, ?)和D₀的高精度计算仍具有重要的工程意义.
实际应用中, 一般是在获得了被测天线近场复数矩阵数据后, 通过近场-远场快速傅里叶变换(fast Fourier transform, FFT)快速算法^[3]来计算天线的归一化方向图函数F(θ, ?), 再由(2)式结合数值积分算法来计算天线的方向系数D₀, 从而形成一整套完备的天线性能数据, 特别是较为准确可靠的方向系数D₀的值, 以用于雷达等实际探测应用. 这就引出了本文要关注的关键问题: 根据(2)式, 结合数值积分算法获得的方向系数准确可靠吗? 与被测天线真实方向系数比较有多大误差? 分析这种误差并给出评估结论, 成为本文主要的研究内容之一.
事实上, 本文作者所在单位已经有专业科技人员开发出数值计算程序—处理近场 (processing near field, PNF), 以计算、处理平面近场测试获得的被测天线近场数据, 包括计算被测天线方向系数D₀. 这一PNF程序在本单位被广泛使用于各种天线的平面近场测试数据分析与处理, 经过实践验证, 算法本身具有公认的基本准确度, 其优点是计算速度很快, 但是本文作者对其算法的具体细节未能获知, 其严格准确度需要通过测试数据计算进行评估.
本文作者也曾利用美国矩阵实验室(Matrix Laboratory, MATLAB)软件开发了相关程序计算D₀值. 结果显示, 针对不同的近场幅相分布数据, 与PNF对比, 两者计算所得的D₀结果在大多数情况下存在一定偏离. 本文的仿真算例中表1—表7给出了共7个算例的结果比对, 其中方法1)乃PNF程序计算方法, 而方法2)—4)乃本文将要叙述的方法. 鉴于此, 迫切需要对天线方向系数的算法进行全面分析评估, 以求既准确又快速地逼近天线的真实方向系数D₀. 平面近场测试环境下的被测天线本身具有真实的方向系数(真值D₀), 但是由于平面近场测试方法具有一定误差, 根据该方法获得的测试数据矩阵计算获得的方向系数仍与真值D₀有一定偏离, 这是近场测试方法具有的系统误差, 这种偏离误差成为本文研究评估的基本误差之一. 基本误差之二为基于平面近场测试数据矩阵的数值积分算法本身具有的积分算法误差. 在整个过程中, 被测天线方向系数的真值D₀是评估的根本参考标准.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(20)式理论真值Dpa	40.9512	0; 应用方法4)对被测天线阵列本身进行数值积分, 结果为40.9512, INVuv误差界1.247607 × 10–8; 积分真值区间: [40.95119, 40.95129], 方法4)计算结果与被测天线理论真值近乎完全吻合, 就积分本身结果而言, 方法4)可作为积分真值参考基准	21.51	解析解, 公式应用范围 受限
1) PNF	41.13	比真值偏大约0.18 dB	1.11	速度最快
2) 本文算法基础积分求和估算	40.9522	比真值偏大约0.001 dB	20.56	速度较快, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	40.9522	比真值偏大约0.001 dB, INVuv误差界1.247085 × 10–8, 积分真值区间: [40.9521, 40.9523]	211.78	速度较慢, 计算时间依赖于FFT点数和测试数据矩阵大小
4) 累加求和被积函数 解析值	40.95218	比真值偏大约0.001 dB, quad2d()算法本身误差, 不存在被积函数值误差, INVuv误差界9.120334 × 10–9, 积分真值区间: [40.952139, 40.952217]	1895.52	速度最慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表1激励均匀分布四种方法方向系数计算比较评估
Table1.Comparison and evaluation of four computation methods of directivity in uniform distribution of electromagnetic excitation.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(21)式理论真值Dpa	13.5284	0; 应用方法4)对被测天线阵列本身进行数值积分, 结果为13.5299, INVuv误差界1.248662 × 10–4; 积分真值区间: [13.52897, 13.53091], 方法4)计算结果与被测天线理论真值接近完全符合, 就积分本身结果而言, 方法4)可作为积分真值参考基准.	24.84	解析解, 公式应用范围 受限
1) PNF	13.73	比真值大约0.2 dB	1.03	速度最快
2) 本文算法基础积分求和估算	13.8459	比真值偏大约0.33 dB	20.34	速度较快, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	13.8598	比真值偏大约0.33 dB, INVuv误差界1.249746 × 10–3, 真值区间: [13.8493, 13.8703]	310.81	速度较慢, 计算时间依赖于FFT点数和测试数据矩阵大小
4) 累加求和被积函数解析值	13.8184	比真值偏大约0.29 dB, quad2d()算法本身误差, INVuv误差界3.546120 × 10–4, 真值区间: [13.8154, 13.8213]	1611.23	速度最慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表7激励均匀分布幅相误差随机正态分布四种方法方向系数计算比较评估
Table7.Comparison and evaluation of four computation methods of directivity in uniform distribution of electromagnetic excitation with normal random errors for amplitude & phase.

本文采用的数值积分数理计算方法本质上是逼近积分真值, 总具有一定的数值误差, 不过积分计算程序函数会给出误差界(error bound), 误差界表示积分逼近值与积分真值的误差绝对值上界, 对同一被积函数, 不同方法的误差界一般不同, 误差界越小, 对应方法获得的积分值越可信, 越逼近积分真值, 也即收敛性越好; 但是数值积分算法无法给出积分真值的确切数值, 只能给出由误差界和积分逼近值决定的一个真值所在区间. 这是数值积分的基本特性.
基于上述两类基本误差, 本文的研究思路为: 建立一个可以计算方向系数真值的理论阵列模型, 先求其严格的方向系数; 再按照实际近场扫描系统的扫描范围与采样间隔, 用计算的方法得到采样点的场强数据(模拟探头测量), 再用模拟探头测量的有限离散数据计算方向系数, 与理论真值比较, 这样在控制基本误差之二最小的情况下可评估算法的基本误差之一的误差程度, 这是对基本误差之一进行评估. 对基本误差之二进行评估主要依赖于MATLAB函数quad2d( )最小积分误差界的搜索和相关计算.
目前对方向系数的计算算法研究可分为两类: 一类为数值方法, 另一类为解析方法. 数值方法中一部分文献在(2)式的基础上引入简化假定^[4-10], 如认为方向图函数只与θ有关, 对积分区间进行截断处理等; 一部分文献采用计算电磁学的数值方法^[11-14]来计算特定属性天线的电磁场, 然后获得方向系数的仿真解. 解析方法中一类文献聚焦于将方向系数的积分公式进行重新解析分析, 如文献[15]利用电场强度的球面波级数展开和傅里叶级数技术重新分析获得方向系数的级数表达式, 文献[16]则利用远场的空间带宽有限特性, 引入采样技术和插值展开获得了方向系数的闭式表达式, 文献[17]对一般空间阵列的方向系数进行了解析分析获得了通用公式. 文献[18]分析了环形微带贴片天线的方向系数逼近解析公式, 也即该公式是近似解析解. 文献[19]专题分析了扫描共线阵列天线的辐射功率和方向系数. 这些文献有一个共同特点: 分别计算辐射强度和辐射功率, 再得出方向系数, 而不是只关注归一化方向图函数F(θ, ?), 这点是与本文不同的地方.
但是, 基于天线平面近场测试理论、数据、FFT算法和数值积分的通用的方向系数计算求解方法研究文献报道很少, 针对基本误差二, 本文试图研究一类基于前述理论和技术的逼近数值积分真值的数学物理计算方法, 对算法计算误差源进行深入分析, 提出尽可能减小误差、逼近积分真值的较准确计算方案, 然后将积分计算结果与近场被测天线真实方向系数进行比较, 以评估这类算法的系统误差, 为理论和实际应用提供参考. 本文研究方法属于数值方法, 不是解析方法, 与上述参考文献研究方向合在一起构成方向系数计算算法的两类方向.

根据平面近场测试理论^[2], 任一天线的方向图函数和方向系数评估都可以基于平面近场测试探头扫描采样点矩阵数据进行, 对天线方向系数的求解转变为基于对应扫描采样点矩阵的近场测试数据矩阵的计算求解.
针对基本误差之一分析, 需要从被测天线模拟出探头扫描有限采样平面处的探头输出值, 根据文献[2], 采用图1坐标系, 探头输出值为
图 1 矩形阵列天线方向图函数求解坐标系
Figure1. The solution coordinate for the pattern of planar array antennas.

$\begin{split} {P_{\rm{B}}}\left( {{{{r}}_{{0}}}} \right) =\, & \frac{{8{{\text{π}}^2}}}{{\omega \mu }}\int_{ - \infty }^\infty \int_{ - \infty }^\infty {k_z}{{G}}\left( {{k_x}, - {k_y}} \right) \\ & \times {{F}}\left( {{k_x},{k_y}} \right)\exp \left( { - {\rm{j}}{{k}} \cdot {{{r}}_{{0}}}} \right){\rm{d}}{k_x} {\rm{d}}{k_y}, \end{split}$

式中k_x和k_y是二重积分的自变量, 有如下关系式:

$\left\{ \begin{aligned} & {{k}} = \hat x{k_x} + \hat y{k_y} + \hat z{k_z}, \\ & {k_z} \!=\! \left\{ \begin{aligned} & \sqrt {{k^2} - {k_x}^2 - {k_y}^2},\quad~ {k_x}^2 + {k_y}^2 \leqslant {k^2}, \\ & - {\rm{j}}\sqrt {{k_x}^2 + {k_y}^2 - {k^2}},~{k_x}^2 + {k_y}^2 > {k^2},\end{aligned} \right. \\& k = {{2{\text{π}}}/ {\lambda,}}\end{aligned} \right.$

G(k_x, –k_y)为探头的平面波谱矢量, F(k_x, k_y)为被测天线的平面波谱矢量, r₀为扫描平面的探头位置矢量, 又因为任一天线的平面波谱矢量A(k_x, k_y)和其远场E(r)有如下关系^[2]:

${{E}}\left( {{r}} \right) = {\rm{j}}\frac{{2{\text{π}}}}{r}{k_{0z}}{{A}}\left( {{k_{0x}},{k_{0y}}} \right)\exp \left( { - {\rm{j}}kr} \right),$

式中k₀为对应远场点r(r, θ, φ)方向(θ, φ)的矢量, 有如下关系式:

$\left\{ \begin{aligned}& {{{k}}_0} = k\hat r = k\left( {\hat x\sin \theta \cos \phi + \hat y\sin \theta \sin \phi + \hat z\cos \theta } \right), \\& k = {{2{\text{π}}} / {\lambda,}}\end{aligned} \right.$

式中k _0z为k ₀的z向分量, k _0z = k cosθ, 其他分量含义类似. 根据(5)式, 任一天线远场矢量方向图函数P(θ, φ)与其平面波谱A(k _x, k _y)有如下关系:

${{P}}\left( {\theta,\phi } \right) = C\cos \theta {{A}}\left( {k\sin \theta \cos \phi,k\sin \theta \sin \phi } \right),$

式中C为与θ, φ无关的部分, k为波数大小.
根据(7)式, 只要知道了被测天线和探头的方向图函数, 即可获知它们对应的位于半径为k的圆内区域的平面波谱函数, 考虑到平面波谱的凋落波成分, 在平面近场测试条件下根据(3)式可获得探头输出值. 为了简化计算和提高测量及计算的精度, 在探头与被测天线同为线极化且方向都平行于z轴的情形下(参见图1, 这也是实际近场测试常见的配置, 此时被测天线和探头主极化都为e_θ方向), 为了(3)式有效采样离散化来模拟扫描平面处的探头响应值, 可设探头矢量方向图函数为cosθe_θ, 被测天线方向图函数为f(θ, φ)e_θ, e_θ为球面坐标系对应θ的基矢, 根据(7)式, (3)式成为

$\left\{ \begin{aligned} & {P_{\rm{B}}}\left( {{{{r}}_{{0}}}} \right) \\ = \,& {C_{\rm{p}}}\int_{ - \infty }^\infty \int_{ - \infty }^\infty \cos \theta {{{e}}_\theta } \cdot \frac{{f\left( {\theta,\phi } \right)}}{{\cos \theta }}{{{e}}_\theta }\\ & \times\exp \left( { - {\rm{j}}{{k}} \cdot {{{r}}_{{0}}}} \right){\rm{d}}{k_x}{\rm{d}}{k_y} \\= \,&{C_{\rm{p}}}\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {f\left( {\theta,\phi } \right)\exp \left( { - {\rm{j}}{{k}} \cdot {{{r}}_{{0}}}} \right){\rm{d}}{k_x}} {\rm{d}}{k_y}} \\= \,& {C_{\rm{p}}}\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {f\left( {u,v} \right)\exp \left( { - {\rm{j}}{{k}} \cdot {{{r}}_{{0}}}} \right){\rm{d}}{k_x}} {\rm{d}}{k_y},} \\u = \,& \sin \theta \cos \phi,\;v = \sin \theta \sin \phi, \\[-10pt]\end{aligned} \right.$

式中C_p为与θ, φ无关的部分, k _x, k _y与θ, φ关系如下:

$\left\{ \begin{aligned} & {k_x} = k\sin \theta \cos \phi = ku, \\& {k_y} = k\sin \theta \sin \phi = kv, \\ & k = {{2{\text{π}}} / {\lambda,}} \end{aligned} \right.$

(8)式在平面近场测试条件下进行采样离散和截断^[20], 省掉C_p, 有如下关系式:

$\left\{ \begin{aligned}& {P_{\rm{B}}}\left( {{{{r}}_{{0}}}} \right) = \sum\limits_{\nu = - {M_x}}^{{M_x} - 1}\sum\limits_{\mu = - {M_y}}^{{M_y} - 1}f\left( {\frac{{{\text{π}}\nu }}{{{L_x}k}},\frac{{{\text{π}}\mu }}{{{L_y}k}}} \right)\\ &\quad\quad\quad\quad\times\exp \left( { - {\rm{j}}{{{k}}_{\nu \mu }} \cdot {{{r}}_{{0}}}} \right),\\& {{{k}}_{\nu \mu }} = \frac{{{\text{π}}\nu }}{{{L_x}}}\hat x + \frac{{{\text{π}}\mu }}{{{L_y}}}\hat y + {\gamma _{\nu \mu }}\hat z, \\& {\gamma _{\nu \mu }} = \sqrt {{k^2} - {{\left( {\frac{{{\text{π}}\nu }}{{{L_x}}}} \right)}^2} - {{\left( {\frac{{{\text{π}}\mu }}{{{L_y}}}} \right)}^2}}, \\& {M_x} \geqslant \frac{{{L_x}}}{{{\lambda / 2}}},~~~~{M_y} \geqslant \frac{{{L_y}}}{{{\lambda / 2}}},\end{aligned} \right.$

式中L_x为扫描平面x向长度的一半, L_y为扫描平面y向长度的一半, M_x为扫描平面x向采样点数的一半, M_y为扫描平面y向采样点数的一半, 为包括所有的传输波谱, M_x和M_y在L_x和L_y确定后需满足(10)式中的不等式.
给定一个已知准确方向图函数f(θ, φ) (即f(u, v))和准确方向系数D₀的被测天线(本文数值算例部分选定一个平面阵列天线作为评估对象, 该天线阵列根据文献[17]可计算出本文条件下的准确方向系数), 根据(10)式求解出探头输出值后, 再根据平面近场测试理论和前述被测天线和探头极化假定后有如下关系式^[2]:

$\left\{ \begin{aligned} & u = \sin \theta \cos \phi,\;v = \sin \theta \sin \phi,\;k = {{2{\text{π}}}/ {\lambda,}} \\ & {f_t}\left( {u,v} \right) = \\ & \Bigg| \sum\limits_{i = 1}^{2{M_x}}\sum\limits_{l = 1}^{2{M_y}}{U_{il}} \exp \left\{ {{\rm{j}}k\left[ {\left( {i - 1} \right){D_x}u+ \left( {l - 1} \right){D_y}v} \right]} \right\} \Bigg|,\end{aligned} \right.\tag{11a}$

式中f_t(u, v)为由被测天线平面近场测试数据计算出的方向图函数, D_x为扫描平面x向探头采样间距, D_y为扫描平面y向探头采样间距. f_t(u, v)与f(u, v)之间存在着一定的系统误差, 本文将在后节数值算例部分对由此导致的方向系数误差给出评估.
(11a)式中U_il为P_B(r₀)在采样点处的离散值, 这一近场测试扫描采样平面矩阵相当于一个平面矩形阵列天线(因为其方向图表达式与对应的单元位于采样点并且单元电压激励与采样点探头输出值相同的平面矩形阵列天线方向图表达式阵因子完全一致). 一般矩形阵列天线的示意图如图1所示, 选取矩形阵列平面为xy平面, 阵列平面每个交点数据代表一个天线单元(单元结构一样, 只是每个单元激发的电压或电流幅相可能不一样). 阵列单元数为N_x × N_y个, 沿x和y轴的单元间距分别为d_x和d_y, 而(θ, ?)为空间电磁场场点r的球坐标方向角.
类似地, 针对图1的结构参数, 单元各向同方向性的该平面天线阵列的方向图函数为^[1]

$\left\{ \begin{aligned} & u = \sin \theta \cos \phi,\;v = \sin \theta \sin \phi,\;k = {{2{\text{π}}} / \lambda }, \\& f\left( {u,v} \right) =\\ & \Bigg| \sum\limits_{i = 1}^{{N_x}}\sum\limits_{l = 1}^{{N_y}} {U_{il}} \exp \left\{ {{\rm{j}}k\left[{( i \!-\! 1){d_x}u \!+\! (l-1){d_y}v} \right]} \right\} \Bigg|.\end{aligned} \right.\tag{11b}$

这里f(u, v)无下标t, 表示这是平面天线阵列的理论方向图函数.
近场测试中, 利用FFT从平面近场数据矩阵得到的方向图函数通常是直角坐标形式F_t(u, v), 这里使用下标t表示该函数是测试数据产生的方向图函数, 与(11a)式的未归一化方向图函数f_t(u, v)相对应. 直角坐标变量(u, v)对应于球坐标中单位球面上某点在xy平面上的投影坐标, 它们的绝对值都不大于1. 注意(u, v)实际上只代表上半球面或下半球面上的空间方向, 满足如下关系式^[3]:

$ \left\{ \begin{aligned} & \left\{ \begin{aligned} & u = \sin \theta \cos \phi, \\ & v = \sin \theta \sin \phi, \end{aligned} \right. \\& {F_t}\left( {\theta,\phi } \right) = {F_t}\left( {u,v} \right).\end{aligned} \right.$

如前所述, 根据平面近场测试理论结果的半空间范围特性, 结合天线实际设计与应用, 可以认为平面天线阵只具有前向辐射, 无需考虑后向辐射, 故F_t(θ, ?)的θ只取区间[0, π/2], 方向系数写为D_pa, 表达式(2)式修正为

$\left\{ \begin{aligned} & {D_{{\rm{pa}}}} = \frac{{4{\text{π}}}}{{\displaystyle\int_0^{2{\text{π}}} {\displaystyle\int_0^{{{\text{π}} / 2}} {{F_t}^2\left( {\theta,\phi } \right)\sin \theta {\rm{d}}\theta {\rm{d}}\phi } } }}, \\ & {D_{{\rm{pa}}}}({\rm{dB}}) = 10\lg \left( {{D_{{\rm{pa}}}}} \right).\end{aligned} \right.$

设有一被积函数H(u, v), 满足

$\begin{split} &H\left( {u,v} \right)\left| {\frac{{D\left( {u,v} \right)}}{{D\left( {\theta,\phi } \right)}}} \right| \\= & {F_t}^2\left( {\theta,\phi } \right)\sin \theta = {F_t}^2\left( {u,v} \right)\sin \theta, \end{split}$

式中$ \left| \cdots \right| $为雅可比行列式, 是被积函数积分元变换的基本因子, 可以写为

$\frac{{D\left( {u,v} \right)}}{{D\left( {\theta,\phi } \right)}} = \left| {\begin{array}{*{20}{c}} {\dfrac{{\partial u}}{{\partial \theta }}}&{\dfrac{{\partial u}}{{\partial \phi }}} \\ {\dfrac{{\partial v}}{{\partial \theta }}}&{\dfrac{{\partial v}}{{\partial \phi }}} \end{array}} \right| = \sin \theta \cos \theta. $

根据(12)式, 有

$\begin{split}& \int_{ - 1}^1 {\int_{ - \sqrt {1 - {u^2}} }^{\sqrt {1 - {u^2}} } {H\left( {u,v} \right){\rm{d}}v} } {\rm{d}}u \\ =\, & \int_0^{2{\text{π}}} {\int_0^{{{\text{π}} / 2}} {{F_t}^2\left( {\theta,\phi } \right)\sin \theta \frac{1}{{\left| {\dfrac{{D\left( {u,v} \right)}}{{D\left( {\theta,\phi } \right)}}} \right|}}\left| {\frac{{D\left( {u,v} \right)}}{{D\left( {\theta,\phi } \right)}}} \right|{\rm{d}}\theta {\rm{d}}\phi } } \\=\, & \int_0^{2{\text{π}}} {\int_0^{{{\text{π}} / 2}} {{F_t}^2\left( {\theta,\phi } \right)\sin \theta {\rm{d}}\theta {\rm{d}}\phi } } = \frac{{4{\text{π}}}}{{{D_{{\rm{pa}}}}}}.\\[-18pt]\end{split} $

根据(14)式、(15)式和(12)式, 有

$\begin{split} H\left( {u,v} \right)\, & = \frac{{{F_t}^2\left( {u,v} \right)\sin \theta }}{{\left| {\dfrac{{D\left( {u,v} \right)}}{{D\left( {\theta,\phi } \right)}}} \right|}} = \frac{{{F_t}^2\left( {u,v} \right)}}{{\left| {\cos \theta } \right|}} \\ &= \frac{{{F_t}^2\left( {u,v} \right)}}{{\sqrt {1 - {u^2} - {v^2}} }},\end{split}$

很显然, (13)式可重写为

$\left\{ \begin{aligned}& {D_{{\rm{pa}}}} = \frac{{4{\text{π}}}}{{{\rm{INVuv}}}}, \\& {\rm{INVuv}} = \int_{ - 1}^1 {\int_{ - \sqrt {1 - {u^2}} }^{\sqrt {1 - {u^2}} } {\frac{{{F_t}^2\left( {u,v} \right)}}{{\sqrt {1 - {u^2} - {v^2}} }}{\rm{d}}v{\rm{d}}u,} } \\& {D_{{\rm{pa}}}}({\rm{dB}}) = 10\lg \left( {{D_{{\rm{pa}}}}} \right).\end{aligned} \right.$

本文基本误差之二问题最后变成逼近计算近场测试数据矩阵方向系数(18)式的问题. 基于PNF程序的算法可以快速计算平面近场测试数据的方向系数, 但所得系数D_pa没有精度控制, 无法估计计算结果的可靠性. 本文将在(18)式基础上引入数值积分精度控制, 实现可靠逼近INVuv的真值, 也就是实现可靠逼近近场测试数据矩阵方向系数积分式真值.

为研究各种计算INVuv的方法, 本文主要采用体现各种算法内核的程序数值仿真形式, 分析误差来源, 给出积分逼近值、误差界和积分真值所在区间, 进行仿真验证、评估. 共讨论四种逼近方法.
1) PNF: 作者所在单位专业科技人员开发的程序, 作为一种近场数据处理程序得到广泛使用, 本文将其作为参考评估基础. 它的方向系数算法类型肯定属于数值积分, 但是因为PNF是编译后可执行程序, 具体算法未知, 只能通过与其他算法的结果进行比较来分析、评价其准确性和快速性.
2) 本文所述的FFT算法^[3]基础积分估算程序, 无精度控制. 其基本点是在已知矩形阵列幅相分布基础上, 利用FFT算法快速算出阵面方向图函数F_t(u, v). 该方法提供了有限的、离散的(u, v)处函数值, 然后基于对函数F_t(u, v)特性的分析来数值求积求和逼近INVuv. 此方法由于积分数据受FFT算法限制, 计算所覆盖的u值范围未完全覆盖积分定义域[–1, 1], 因此INVuv在(u, v)定义域边缘处的积分值被忽略掉, 导致INVuv偏小. 另外, 近似数值求积求和也引入了一定计算误差. 此方法的结果误差界要通过与以下两种方法尤其是第4)种方法的比较才能有一定的估计, 正是由于它没有误差控制, 因此计算时间在第2), 3), 4)种方法中是最短的, 也即速度在三种方法中最快.
3) 在第2)种方法基础上, 引入MATLAB二维插值函数interp2()和二维数值积分基本函数quad2d()的积分程序. 积分基本函数quad2d()提供误差容差控制, 积分方法是平铺(tiled), 不是两次一重积分迭代(iterated), 平铺的基本方法是将非矩形平面积分区域通过内部变换函数变换为一批小矩形积分区域再进行数值求积, 而且quad2d()函数能处理积分区域边界处的被积函数暇点, 这种能力正好适应(18)式中被积函数在积分边界处的暇点特性, 具体信息可参阅文献[21]. 一般来说, 这一方法的积分结果比第2)种方法的结果更逼近积分真值, 它在第2)种方法得到的F_t(u, v)离散值基础上进行二维插值, 来估算所需插值点处的函数值. 这些插值点由求积函数quad2d()根据积分算法需要产生, 从而使本方法获得了更多的被积函数值. 这些插值点处的被积函数值相对于函数值真值来说仍存在一定插值误差, 虽然经过MATLAB的三次Hermite或样条插值方法可尽量逼近函数值真值, 但插值误差不可能完全消除.
4) 基于累加求和的F_t(u, v)积分计算, 称为累加求和被积函数解析值方法. 注意到被积函数F_t(u, v)的“未归一化”累加求和解析表达式f _t(u, v)即为(11a)式:

$\left\{ \begin{aligned}& u = \sin \theta \cos \phi,\;v = \sin \theta \sin \phi,\;k = {{2{\text{π}}} /{\lambda,}} \\& {f_t}\left( {u,v} \right) = \\ & \Bigg| \sum\limits_{i = 1}^{2{M_x}}\sum\limits_{l = 1}^{2{M_y}}{U_{il}} \exp \left\{ {{\rm{j}}k\left[ {( i \!-\! 1 ){D_x}u + ( {l \!-\! 1} ){D_y}v} \right]} \right\} \Bigg|,\end{aligned} \right.$

式中U_il为P_B(r₀)在采样点处的离散值, U_il也可理解为探头采样平面对应的平面矩形阵列天线离散单元处激发的电压或电流的幅相复数值. F_t(u, v)与f_t (u, v)关系为

${F_t}\left( {u,v} \right) = \frac{{{f_t}\left( {u,v} \right)}}{{\max \left( {{f_t}\left( {u,v} \right)} \right)}}.$

很显然, 要得到F_t(u, v)的值, 必须知道f_t(u, v)的最大值. 有两种方法: 一是由第2)种方法得到的主瓣位置(u₀, v₀)(即f_t(u, v)最大值所在方向)代入(19)式求得. 这里, 得到的最大值因为是离散数据, 在第2)种方法二维FFT计算点数(一般取为2¹⁰ × 2¹⁰)较多的情形下只有一定的微小误差; 二是理论分析得到f_t(u, v)函数最大值, 比如图1相控阵的最大值方向(u₀, v₀)在某些情形下可有理论分析值, 代入(11a)式即可获得f_t(u, v)函数最大值, 这一方法没有函数值误差. 虽然理论分析法精度更高, 但是理论分析值可应用情形不多, 第一种方法是(11a)式获得最大值的通用方法.
根据(11)式、(18)式和(19)式, 本方法引入MATLAB二维数值积分基本函数quad2d()的积分程序来计算积分值. 该方法与第3)种方法不同之处在于被积函数即为根据近场测试理论得出的F_t(u, v)解析表达式, 而第3)种方法中, 被积函数F_t(u, v)的值来源于被插值(u, v)处的估计值, 引入了插值误差. 本方法避免了数值计算F_t(u, v)这一影响积分精度的关键因子时引入的误差, 直接引入F_t(u, v)累加求和解析表达式. 该方法根本误差源只有一项: 来源于MATLAB二维数值求积函数quad2d()的计算精度误差, 因此在同样容差精度条件下其结果是最逼近积分真值的, 虽然计算耗时最长.
本文基于MATLAB和其包含的二维数值积分基本函数quad2d()开发了数值积分计算天线方向系数的程序, 程序中根据quad2d()函数返回的信息不断调节计算容差, 实现最小容差定位和最佳积分逼近值, 同时给出误差界和积分真值区间. 误差界越小, 积分值越逼近积分真值. 需要说明的是, 此处积分真值与被测天线方向系数真值是不同的概念, 两者之间如前所述有一定的系统偏离.
鉴于以上分析和论证, 方法4)可作为对于同一被积函数不同数值积分算法关于积分真值的参考基准, 数值算例和结果可以予以验证, 算例结果应体现在同样被积函数条件下方法4)的误差界最小, 方法4)是积分真值的最佳逼近.
综上所述, 对于第1)种方法, 只能将其计算结果与其他三种方法进行误差比较来给出评估; 第2), 3)种方法可以用第4)种方法的结果作为评估积分真值本身的参考基准, 这点将由下节的数值算例和结果验证.

为了比较这四种方法的结果, 本文选取某个矩形阵列天线作为被测天线, 矩形阵列天线示意图如图1所示. 矩形阵列结构参数为: N_x = 60, N_y = 40, d_x = 0.65λ, d_y = 0.65λ, λ为计算辐射电磁波波长. 结合实际精度和速度需要, 二维FFT计算点数一般设为2¹⁰ × 2¹⁰. 这一阵列将作为以下7个算例的平面近场测试被测天线(单元各向同方向性), 它的方向图函数和方向系数真值都可获得. 文献[17]给出了任意立体阵列的方向系数解析解, 自然本被测天线也包含在内, 结合平面近场测试的测试条件和假定, 单元各向同方向性的相位线性递增分布的平面天线阵列方向系数有如下表达式^[17]:

$\begin{split}& {D_{{\rm{pa}}}} = 10 \times \\& \lg \!\left[ \frac{{2{{\left( {\displaystyle\sum\limits_{n = 1}^N {{I_n}} } \right)}^2}}}{{\rm real} \left( \displaystyle\sum\limits_{m = 1}^N \displaystyle\sum\limits_{n = 1}^N {{I_m}{I_n}} \exp ( {{\rm{j}}( {\beta _m} \!-\! {\beta _m} )} )\dfrac{{\sin ( {k{r_{mn}}} )}}{{k{r_{mn}}}} \right)} \! \right],\end{split}$

式中N表示天线阵列的单元总数, I表示单元激励电流或电压振幅值, β表示单元激励电流或电压相位值(弧度), k为波数, m和n是单元下标, r_mn表示单元m和单元n之间的空间距离, real表示取复数实部. 下文算例1至算例4天线阵列属于这一类别.
类似地, 单元各向同方向性的相位任意分布的天线阵列的方向系数有如下表达式^[17]:

${D_{{\rm{pa}}}} \!=\! 10\lg \left[{\frac{{\operatorname{real} \left( {2\displaystyle\sum\limits_{m = 1}^N {\displaystyle\sum\limits_{n = 1}^N {{I_m}{I_n}} } \exp \left\{ {{\rm{j}}k\sin {\theta _0}\left[ {\left( {{x_m} - {x_n}} \right)\cos {\phi _0} + \left( {{y_m} - {y_n}} \right)\sin {\phi _0}} \right]} \right\} \cdot \exp \left( {{\rm{j}}\left( {{\beta _m} - {\beta _m}} \right)} \right)} \right)}}{{\operatorname{real} \left( {\displaystyle\sum\limits_{m = 1}^N {\displaystyle\sum\limits_{n = 1}^N {{I_m}{I_n}} } \exp \left( {{\rm{j}}\left( {{\beta _m} - {\beta _m}} \right)} \right)\dfrac{{\sin \left( {k{r_{mn}}} \right)}}{{k{r_{mn}}}}} \right)}}} \right],$

式中x和y是单元的直角坐标值, (θ₀, ?₀)为根据第3节第2)种方法求解出的方向图函数最大值球坐标方向角, 其他物理量含义同(20)式. 下文算例5至算例7天线阵列属于这一类别.
被测天线与探头扫描采样平面距离d₀ = 3λ, (10)式中有关结构参数如下: L_x = N_x·d_x = 39λ, L_y = N_y·d_y = 26λ, 采样间距D_x = 0.5λ, D_y = 0.5λ, M_x, M_y由(10)式给出为M_x =78, M_y = 52.
本文共针对7个算例展开数值仿真计算, 7个算例给定的电压激励模式都是针对矩形阵列被测天线, 电压激励模式选择遵循如下原则: 尽量覆盖方向系数高、低的情形, 同时结合天线电讯设计工程实际中常用的电压激励幅相分布, 以全面评估、验证、判断这四种方法哪一种计算的方向系数D_pa的准确性最好. 计算过程简介如下: 先由(10)式模拟算出被测天线在探头采样平面的激励电压值, 再由(11a)式和(18)式等对模拟探头测试数据矩阵进行四种方法的数值积分计算, 并将结果与被测天线的方向系数理论真值进行比较, 从而获得近场测试数据方法计算方向系数的误差评估结果, 具体如以下7个算例所示.
本文的计算机软硬件环境配置如下: 操作系统为Windows 7 SP 1 X64旗舰版, 内存4 GB, 处理器Intel Core2 Duo CPU P9400@2.4GHz, MATLAB版本为8.1.0.604(R2013a) 32-bit(Win32). 第1)种方法为Windows可执行程序, 计时采用秒表; 第2), 3), 4)种方法为MATLAB程序, 可在程序中精确计时, quad2d()函数初始绝对误差控制为10^–7, 该值不同, 第3), 4)种方法耗时会不同. 四种方法都在上述同一计算平台进行验证评估.
算例1　激励模式满足幅度均匀分布、相位同相的条件. 理论分析可知主瓣波束指向(即方向图函数最大值方向)满足(θ₀, ?₀) = (0°, 0°), 即指向阵面法向. 此时, 数值积分时不存在f_t(u, v)函数最大值误差. 这种幅相分布称为幅相均匀分布, 在工程实际应用中, 雷达天线发射通常采用这种幅相分布. 其特点是方向图波束宽度很窄, 副瓣电平为–13.3 dB左右, 方向系数较大. 对此分布, 利用(20)式被测天线方向系数理论真值公式和所列四种方法对近场测试矩阵的计算和评估结果如表1所列. 计算结果表明, 天线阵列方向系数达到40 dB以上, 实现了高方向性, 与均匀分布的规律一致.
算例2　激励模式满足幅度均匀分布条件, 但相位则依x轴、y轴可分离线性递增分布(广义递增, 相邻单元相位增量为负也叫递增). 这种相位分布是典型的相控阵扫描相位分布, 在工程实际中广泛应用, 能实现相控阵波束扫描功能. 本例设定波束扫描指向为(θ₀, ?₀) = (40°, –60°), 即在算例1的基础上实现相控阵波束扫描, 也就是主瓣波束扫向空间方向(θ₀, ?₀) = (40°, –60°). 相比于(θ₀, ?₀) = (0°, 0°)的情况, 这一相控阵扫描阵列方向系数要小, 且扫描角偏离(θ₀, ?₀) = (0°, 0°)越大, 方向系数就越小. 针对此分布的四种方法计算与评估结果如表2所列. 计算结果表明, 本例相控阵扫描阵列方向系数比算例1的情形下降2.0 dB左右.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(20)式理论真值Dpa	38.9473	0; 应用方法(4)对被测天线阵列本身进行数值积分, 结果为38.9474, INVuv误差界2.257101 × 10–8; 积分真值区间: [38.94733, 38.94746], 方法4)计算结果与被测天线理论真值近乎完全符合, 就积分本身结果而言, 方法4)可作为积分真值参考基准	20.92	解析解, 公式应用范围 受限
1) PNF	39.33	比真值偏大约0.38 dB	1.12	速度最快
2) 本文算法基础积分求和估算	39.2175	比真值偏大约0.27 dB	20.92	速度较快, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	39.2142	比真值偏大约0.27 dB, INVuv误差界1.248032 × 10–7, 积分真值区间: [39.21385, 39.21457]	201.61	速度较慢, 计算时间依赖于FFT点数和测试数据矩阵大小
4) 累加求和被积函数解析值	39.2023	比真值偏大约0.25 dB, quad2d()算法本身误差, INVuv误差界2. 873319 × 10–8, 积分真值区间: [39.20225, 39.20241]	1276.35	速度最慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表2激励幅度均匀相位扫描分布四种方法方向系数计算比较评估
Table2.Comparison and evaluation of four computation methods of directivity in uniform amplitude & linear scanning phase distribution of electromagnetic excitation.

算例3　激励的幅相遵循常用的Taylor分布. 在相同副瓣电平设计要求的条件下, Taylor分布易于在实际工程中实现. 工程实际中, 此分布一般用于天线接收时的幅相分布设计. 本例假定幅度在x方向为–35 dB副瓣电平的Taylor分布, 在y方向为–25 dB副瓣电平的Taylor分布, 相位扫描波束指向为(θ₀, ?₀) = (40°, –60°). 针对此分布的四种方法计算和评估结果如表3所列.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(20)式理论真值Dpa	37.8093	0; 应用方法4)对被测天线阵列本身进行数值积分, 结果为37.809329, INVuv误差界6.831371 × 10–10; 积分真值区间: [37.809327, 37.809330], 方法4)计算结果与被测天线理论真值近乎完全符合, 就积分本身结果而言, 方法4)可作为积分真值参考基准	21.47	解析解, 公式应用范围 受限
1) PNF	38.40	比真值偏大约0.59 dB	1.04	速度最快
2) 本文算法基础积分求和估算	38.1953	比真值偏大约0.39 dB	20.65	速度较快, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	38.2152	比真值偏大约0.4 dB, INVuv误差界1.249154 × 10–8, 真值区间: [38.21517, 38.21523]	189.59	速度较慢, 计算时间依赖于FFT点数和测试数据矩阵大小
4) 累加求和被积函数解析值	38.190513	比真值偏大约0.38 dB, quad2d()算法本身误差, INVuv误差界1.248422 × 10–9, 真值区间: [38.19051, 38.190516]	1465.8	速度最慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表3激励Taylor分布四种方法方向系数计算比较评估
Table3.Comparison and evaluation of four computation methods of directivity in Taylor distribution of electromagnetic excitation.

算例4　FFT计算精度的影响. 在算例3基础上, 将FFT算法分辨率进一步提高, 一维FFT点数由2¹⁰变为2¹¹, 注意本文所指FFT实际上是二维FFT运算. 因为F_t(u, v)的真实离散值数量同步增大, 数值积分计算结果一般会更加逼近真值. 但是, FFT点数不宜太多, 否则计算逼近时间太长, 不便实用, 如本例第3)种方法耗时太多. 在此情况下, 考察积分结果的变化, 四种方法计算和评估结果如表4所列.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(20)式理论真值Dpa	37.8093	0; 应用方法4)对被测天线阵列本身进行数值积分, 结果为37.809329, INVuv误差界6.831371 × 10–10; 积分真值区间: [37.809327, 37.809330], 方法4)计算结果与被测天线理论真值近乎完全符合, 就积分本身结果而言, 方法4)可作为积分真值参考基准.	21.18	解析解, 公式应用范围 受限
1) PNF	38.36	比真值偏大0.56 dB, 结果因FFT点数变化与表3相比略有变化	1.09	速度最快, 本例FFT点数为211 × 211, 其他算例FFT点数都为210 × 210
2) 本文算法基础积分求和估算	38.1533	比真值偏大0.34 dB	82.18	计算时间加长, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	38.1535	比真值偏大0.34 dB, INVuv误差界1.805049 × 10–8, 真值区间: [38.153468, 38.153549]	900.46	计算时间长, 计算时间依赖于FFT点数和测试数据矩阵大小, 实用性减弱
4) 累加求和被积函数解析值	38.13822	比真值偏大0.33 dB, quad2d()算法本身误差, INVuv误差界1.247355 × 10–9, 真值区间: [38.138217, 38.138222], 结果因FFT点数变化造成最大值略有变化, 最终结果与表3相比略有变化	1468.83	速度较慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表4激励Taylor分布FFT精度提高四种方法方向系数计算比较评估
Table4.Comparison and evaluation of four computation methods of directivity in Taylor distribution of electromagnetic excitation with raised precision of FFT algorithm.

算例5　幅度、相位随机“均匀”分布误差的影响. 以上算例为高方向性的计算实例, 在其基础上, 可以研究低方向性情况下的性能. 一种可能性为在上述幅相分布的基础上引入随机幅相分布误差, 因此本算例在算例3的基础上引入幅度、相位随机“均匀”分布误差. 四种方法计算和评估结果如表5所列. 结果表明, 随机均匀分布幅相误差对矩形阵列天线的方向系数数值影响很大(方向系数降低24 dB左右), 实现了低方向性矩形阵列幅相分布. 注意到PNF更加逼近方向系数真值.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(21)式理论真值Dpa	13.3604	0; 应用方法4)对被测天线阵列本身进行数值积分, 结果为13.3604, INVuv误差界3.473624 × 10–5; 积分真值区间: [13.3602, 13.3607], 方法4)计算结果与被测天线理论真值近乎完全符合, 就积分本身结果而言, 方法4)可作为积分真值参考基准.	25.09	解析解, 公式应用范围 受限
1) PNF	13.63	比真值大约0.27 dB	1.03	速度最快
2) 本文算法基础积分求和估算	13.6866	比真值偏大约0.33 dB	20.25	速度较快, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	13.6806	比真值偏大约0.32 dB, INVuv误差界7.892290 × 10–4, 真值区间: [13.6742, 13.6870]	279.98	速度较慢, 计算时间依赖于FFT点数和测试数据矩阵大小
4) 累加求和被积函数解析值	13.6580	比真值偏大约0.30 dB, quad2d()算法本身误差, INVuv误差界1. 382844 × 10–4, 真值区间: [13.6569, 13.6592]	1259.79	速度最慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表5激励Taylor分布幅相误差随机均匀分布四种方法方向系数计算比较评估
Table5.Comparison and evaluation of four computation methods of directivity in Taylor distribution of electromagnetic excitation with uniform random errors for amplitude & phase.

算例6 幅度、相位随机“正态”分布误差对Taylor分布的影响. 在算例3的基础上, 引入幅度、相位随机“正态”分布误差. 四种方法计算和评估结果如表6所列. 结果表明, 随机正态分布幅相误差对方向系数数值影响显著(方向系数减小25 dB左右), 实现了低方向性矩形阵列幅相分布. 注意到PNF更加逼近方向系数真值、速度最快, 其他方法逼近积分真值规律和计算速度上与高方向系数的情形类似.

方法	方向系数计算值/dB	误差分析、估计	计算时间/s	方法特点和局限
(21)式理论真值Dpa	12.6206	0; 应用方法4)对被测天线阵列本身进行数值积分, 结果为12.6213, INVuv误差界1.249209 × 10–4; 积分真值区间: [12.6205, 12.6220], 方法4)计算结果与被测天线理论真值近乎完全符合, 就积分本身结果而言, 方法4)可作为积分真值参考基准	24.84	解析解, 公式应用范围 受限
1) PNF	12.90	比真值偏大约0.28 dB	1	速度最快
2) 本文算法基础积分求和估算	13.0038	比真值偏大约0.38 dB	20.40	速度较快, 计算时间依赖于FFT点数和测试数据矩阵大小
3) 二维插值估计被积函数	13.0030	比真值偏大约0.38 dB, INVuv误差界1.249323 × 10–3, 真值区间: [12.9943, 13.0116]	297.41	速度较慢, 计算时间依赖于FFT点数和测试数据矩阵大小
4) 累加求和被积函数解析值	12.9737	比真值偏大约0.35 dB, quad2d()算法本身误差, INVuv误差界2.116856 × 10–4, 真值区间: [12.9722, 12.9751]	1561.10	速度最慢, 计算时间依赖于近场测试数据矩阵大小, 本例为104 × 156, 最佳逼近积分真值

表6激励Taylor分布幅相误差随机正态分布四种方法方向系数计算比较评估
Table6.Comparison and evaluation of four computation methods of directivity in Taylor distribution of electromagnetic excitation with normal random errors for amplitude & phase.

算例7 幅度、相位随机“正态”分布误差对均匀分布的影响. 在算例1基础上引入幅度、相位随机“正态”分布误差, 进一步考察较小方向系数的矩形阵列天线. 四种方法计算和评估结果如表7所列. 四种方法逼近积分真值规律和计算速度上与算例6的情形类似, 正态分布误差对均匀幅相分布的影响比对Taylor分布的影响大, 本例方向系数下降27 dB左右. 注意到PNF更加逼近方向系数真值.
上述7个算例的计算仿真和评估验证表明, 计算天线方向系数的数值积分四种方法各有特点和局限, 如下所述.
第1)种方法PNF对高方向系数的计算误差相对较大, 误差较大值接近0.6 dB, 对低方向系数的计算逼近方向系数真值, 速度总是最快. 这表明PNF计算高方向系数天线时, 算法对应INVuv的积分结果比INVuv积分真值偏小, 造成计算的方向系数值比第4)种方法计算的方向系数值偏大, 因为第4)种方法对应INVuv的积分结果(如算例1—算例7误差分析、估计部分结果证明)总是INVuv积分本身真值的最佳逼近(近乎完全符合). 但是所有方法(包括第4)种方法)计算出的方向系数总是比方向系数真值偏大, 这就造成PNF计算的方向系数更加远离方向系数真值, 如算例1—算例4所示. PNF计算低方向系数天线时, 算法对应INVuv的积分结果比INVuv积分真值偏大, 造成计算的方向系数值比第4)种方法计算的方向系数值偏小, 同样因为第4)种方法对应INVuv的积分结果总是INVuv积分本身真值的最佳逼近(近乎完全符合). 但是所有方法(包括第4)种方法)计算出的方向系数总是比方向系数真值偏大, 这样就造成PNF计算的方向系数更加逼近方向系数真值, 如算例5—算例7所示.
第2)种方法速度稍慢, 计算时间依赖于FFT点数和测试数据矩阵大小, 对于测试矩阵数据总数在16000 (104 × 156 = 16224)规模, 一般计算时间在20 s左右; 第3)种方法速度更慢, 基于第(2)种方法已知的方向图函数值, 引入MATLAB二维插值函数interp2()和二维数值积分基本函数quad2d(), 采用三次Hermite插值或样条插值和二维数值积分精度控制, 一般来说使计算结果更加逼近积分真值, 第4)种方法根据前面小节的理论分析论证和本小节的数值算例仿真结果验证, 误差界一般情况下最小, 因此可作为积分真值最佳逼近, 也可作为其他三种方法积分计算结果本身准确性评估的基础, 该方法根本误差源只有quad2d()的计算误差, 但是耗时最长.
所有方法计算结果都比被测天线方向系数真值偏大. 其中PNF对于高方向性天线计算误差最大值达0.6 dB; 对于低方向性天线, 计算结果更加逼近方向系数真值. 第2), 3), 4)种方法计算误差最大值为0.4 dB, 在这三种方法中, 第4)种方法误差最小. 在四种方法中, 第4)种方法总是积分本身真值的最佳逼近, 但是耗时最长.

基于平面近场测试理论获得的矩形阵列幅相数据, 利用MATLAB插值函数和数值积分函数对被测天线的方向系数进行了逼近求解. 提出了三种新的近场测试数据数值积分计算方法, 并与已有的PNF计算程序进行了仿真数据对比, 给出了这些方法的特点和局限. 结果显示: 所有方法计算结果都比被测天线方向系数真值偏大. 其中第1)种方法PNF对于高方向性天线计算误差最大值达0.6 dB; 对于低方向性天线, 计算结果更加逼近方向系数真值. 第2)、第3)、第4)种方法计算误差最大值为0.4 dB, 在这三种方法中, 第4)种方法误差最小. 在四种方法中, 第4)种方法总是积分本身真值的最佳逼近, 但是耗时最长. 第1)种方法PNF速度总是最快. 第2)种方法FFT算法基础积分估算程序速度稍慢, 但是计算误差对于高方向系数情形较小. 第3)种方法基于MATLAB插值函数interp2()和数值积分函数quad2d()的积分程序一般来说比第2)种方法计算精度更高, 速度较慢, 积分结果一般来说相对于第2)种方法更加逼近方向系数真值. 第4)种方法累加求和被积函数解析值方法可作为方向系数积分式本身真值的最佳逼近, 相比于第2)种、第3)种方法离被测天线方向系数真值更近, 但是耗时最长, 不大能满足快速响应的需求. 本工作是对基于平面近场测试数据的一类天线方向图和方向系数计算技术的系统评估, 对天线方向系数的计算数理方法研究具有理论和应用价值.