0   引　言

多尺度目标电磁散射特性的快速计算与分析，在军事和民用领域都有很大的需求与应用场景，对雷达探测与成像、目标识别、隐身和反隐身技术、微波遥感等技术的实现具有关键作用。在实际应用场景中，由于制作工艺、人为因素、环境因素等影响导致电磁系统的不确定性是普遍存在的^[1]。同时，大多数非合作目标经常由于种种不可抗原因，导致能够获取的目标网格信息^[2-3]、目标材料信息^[4-5]较为有限，造成信息参数不完备的情况。因此，由于目标外形、材料属性等高维度不确定性因素^[6-7]使得目标电磁散射特性也充满了多重不确定性。

目前，对含有高维度不确定性因素的情况，多个不确定性因素间的分级关系尚不明确，不确定因素对散射特性的影响关系缺乏研究。在过去几十年中，蒙特卡洛(Monte Carlo，MC)方法^[8]和准MC法^[9]得到了广泛的应用，其通过多次随机采样，逐次计算随机外形、随机材料参数的目标电磁散射特性，来获取不确定性参数目标的电磁统计特性，该方法的本质是将不确定性问题转化为确定性问题进行计算。因此，通常将MC方法的结果作为标准值，用于验证其他数值方法、近似方法计算结果的准确性。然而，基于MC的方法非常耗时，尤其针对多维问题，计算时间会显著增加。

随后，针对简单的不确定性问题和多参数问题，研究人员提出了广义多项式混沌(generalized polynomial Chaos，gPC)法^[10-13]。发展到2018年，gPC的理论框架才逐渐完善，主要包括随机伽辽金方法(stochastic Galerkin method，SGM)和随机点配置方法(stochastic collocation method，SCM)。其中SGM将随机方程转化为一个耦合且增广的确定性方程，再利用确定性方法求解，通过有限阶广义多项式使得计算误差较小，但是该方法属于侵入性方法，在计算过程中会改变求解器。SCM能够在不改变求解器的情况下给出不确定性问题的计算结果，但该方法的计算精度低于SGM。这两种方法各有利弊，并且对于复杂目标的多维随机问题，gPC方法存在耦合度高和计算资源消耗大的问题^[14-15]。

近年来，基于有限差分时域法^[16-17]和矩量法(method of moments, MoM)^[18]，研究人员提出了几种可以计算参数不确定性的数值估算方法，这类方法统称为扰动法^[19-22]。该类方法利用非均匀有理B样条(nonuniform rational B-spline surface，NURBS)建模技术^[23-25]实现目标外形的构建，并通过引入随机变量任意改变目标的外形。在扰动法的求解过程中，由于不需要逐次填充阻抗矩阵和逐次求逆，因此所需的计算时间大大减少。

2019年，王珂琛等人^[19]通过将MoM中的阻抗矩阵、右边向量、等效电流等展开为关于外形向量的一阶泰勒级数(1-st Taylor expansion)，再通过数学公式求导的方式，计算阻抗矩阵、右边向量关于外形向量的导数，进而得到具有外形不确定的理想电导体(perfect electrical conductor，PEC)目标的电流变化量。2021年，何姿等人^[20]提出了基于混合场积分方程(combined field intergral equation, CFIE)的二阶泰勒级数展开 (2-nd Taylor expansion) 方法，相较于一阶泰勒级数，该方法能够准确计算更大的外形扰动量，主要解决金属目标的外形不确定的问题。2022年，李世玺等人^[21]基于体面积分方程(volume-surface integral equation, VSIE)，提出了外形、介电参数不确定的电磁散射分析方法，其能够解决目标外形、介电参数同时存在不确定性时的目标散射特性统计分析。2023年，何姿等人^[22]又提出了基于旋转对称体(body of revolution, BoR)非均质介质体的几何外形和介电参数的高维不确定性分析方法，当引入微小的扰动量时，计算结果能够与MC方法保持良好的一致性。但是，随着引入的随机扰动量的增加，上述方法的统计误差逐渐增大。另一方面，扰动法在使用过程中存在较多难点。例如，阻抗矩阵、右边向量关于外形向量的数学公式求导过程非常繁琐；编写程序的过程中也极易出错。因此，如何寻求更简洁高效的求导方法，并进一步扩大扰动区间，为后续的研究提供了方向。

在频率和角度等标量域的研究过程中，发现渐近波形估计(asymptotic waveform evaluation, AWE)技术^[26-28] 以 Padé 有理函数取代泰勒级数，面对相同的扰动量其计算精度高于泰勒技术，因此AWE技术也被认为是实现更好性能的可行方法。近年来 AWE 技术得到不断改进，一些改进方法被提出，如 GAWE、MGAWE 和 WCAWE ^[29-32]。这些方法主要应用于频率、角度等标量参数的外推^[33-35]。在一些研究中^[36-37]，AWE也被用于结合快速多极子的混合有限元边界积分方法(finite element-boundary integral the multilevel fast multipole algorithm, FE-BI-MLFMA)，实现金属和介质目标的宽带散射特性的加速计算。

综上所述，目前现有的扰动方法无法对较大扰动量的不确定性问题进行准确、高效的计算。本文首先介绍了NURBS建模技术，用于实现任意模型的目标构建，以及外形、介电参数不确定性目标的参数提取技术。然后，介绍了基于泰勒级数展开的不确定性目标电磁散射特性的计算方法。进一步地，为了扩大扰动法能够准确计算的扰动量，并且节省求导时间，将AWE技术和伪谱法(pseudospectral derivative method, PSDM)联合应用至外形、介电常数的不确定性分析中。AWE技术使用高阶Padé有理函数代替低阶泰勒级数展开，可使得算法适应目标外形动态变化范围更大的不确定性，而使用PSDM求导可实现更高的计算效率。最后，利用反向传播(back propagation, BP)神经网络实现对不确定性目标散射特性的置信度和不确定度的智能化计算。

1   不确定性目标的NURBS建模技术

理想情况下分析电磁散射问题时，目标的几何外形、介电常数等参数信息不会随时间、外界环境的变化发生改变。但实际的生产生活中，外界环境中存在诸多不可控因素，比如制作工艺、日常损耗、人为因素、特殊的自然环境等，导致目标的形状、涂敷材料等参数发生未知的变化，通常将这类问题归为不确定性问题。为模拟现实生活中出现的不确定性问题，以目标外形不确定性为例，通过给目标外形引入服从均匀分布的随机变化，等效真实的目标外形不确定性场景。

对于具有不确定性目标的几何建模，需要在建模后的分析过程中对目标模型进行方便的控制，并且能够描述目标的结构。对于导弹、飞机等复杂目标，其几何外形具有不规则性，采用传统的建模方法难以实现对目标几何外形的不确定性建模。NURBS建模技术是一种基于B样条基函数的几何建模技术，可以通过一组参数曲面片组合逼近目标模型，其使用的单元很少而模拟精度很高，并且根据理论公式，模型的形状可以通过少量的控制点坐标变化实现相应改变。以立方体目标为例，如图1所示，通过改变1号控制点的坐标，实现目标外形的随机变化。因此，本文采用NURBS建模技术，即采用NURBS实现对模型的几何建模，使用的NURBS建模软件为犀牛软件。

图  1  外形不确定的立方体模型

Fig.  1  Cubic model with shape uncertainty

下载: 全尺寸图片 幻灯片

NURBS的形状可以由若干个控制点通过一个双变量分段有理函数控制，如图2所示，曲面上任意一点的坐标可以表示为

图  2  含有6个控制点的NURBS

Fig.  2  NURBS surface with 6 control points

下载: 全尺寸图片 幻灯片

式中：U和V分别为u，v方向上控制点的个数；$w_{i,j}$为相应的权重；${N_{i,p}}(u)$和${N_{j,q}}(v)$分别为由节点矢量${\boldsymbol{U}} = [{u_0},{u_1},\cdots ,{u_{n + k + 1}}]$和${\boldsymbol{V}} = [{v_0},{v_1},\cdots ,{v_{n + k + 1}}]$根据Cox-DeBoor递推公式而得到的p阶、q阶规范B样条基函数；${{{P}}_{ij}} \left( {{P_{ijx}},{P_{ijy}},{P_{ij{\textit{z}}}}} \right)$为x, y, z方向上控制点的坐标。

${N_{i,p}}(u)$的定义如下式所示：

令分段有理基函数为

则NURBS任意一点坐标表达式可写为

根据式(4)发现，通过调整控制点${{{P}}_{ij}}$的坐标能够改变NURBS的外形。而且采用NURBS建模技术建立目标几何模型，当目标的局部几何外形发生变化时，仅须改变相关的控制点坐标，并且不影响其他部分的几何形状。

2   不确定性目标参数提取技术

2.1   外形不确定目标参数提取技术

将控制点位置坐标设置为外形向量${\boldsymbol{\alpha}} = \left[ {{\alpha _1},{\alpha _2}, \cdots ,{\alpha _t}} \right]$(t为外形向量元素的个数)，即通过外形向量${\boldsymbol{ \alpha }}$来描述目标几何外形的不确定性。

NURBS上任意一个点的x, y, z方向坐标可以表示为

为实现NURBS与MoM程序相结合，在NURBS上放置RWG(Rao-Wilton-Glisson)基函数，每个三角形顶点的坐标用相应的NURBS的控制点坐标表示。进一步地，RWG基函数${{\boldsymbol{f}}_n}({\boldsymbol{r}})$的几何信息均可以用外形向量${\boldsymbol{\alpha }}$表示：

式中：${l_n}({\boldsymbol{\alpha}} )$为相邻三角形公共边的长度；$A_n^ + ({\boldsymbol{\alpha}} )$和$A_n^ - ({\boldsymbol{\alpha }})$为三角形$T_n^ +$,$T_n^ -$的面积；${\boldsymbol{\rho }}_n^ +$为三角形$T_n^+$的顶点指向该三角形上的场点的矢量；${\boldsymbol{\rho}} _n^ -$为三角形${T_n^ -}$上的场点指向该三角形的顶点的矢量。结合MoM原理，用带有外形向量${\boldsymbol{\alpha }}$的RWG基函数离散积分方程后，MoM的矩阵方程可以看作是关于外形向量${\boldsymbol{\alpha }}$的函数等式，如公式(7)所示：

式中：${\boldsymbol{Z}}({\boldsymbol{\alpha}} )$和${\boldsymbol{b}}({\boldsymbol{\alpha}} )$分别为关于外形向量${\boldsymbol{\alpha}}$的阻抗矩阵和右边向量；${\boldsymbol{x}}({\boldsymbol{\alpha}} )$为关于外形向量${\boldsymbol{\alpha}}$的电流。

2.2   介质不确定目标参数提取技术

相较于2.1节中外形向量${\boldsymbol{ \alpha }}$的参数提取过程，介电参数的提取过程较为简单。首先，构建介电参数向量${\boldsymbol{\varepsilon}} = \left[ {{\varepsilon _1},{\varepsilon _2}, \cdots ,{\varepsilon _l}} \right]$(l为介电参数种类的个数)，即通过改变介电参数向量${\boldsymbol{\varepsilon}}$来描述目标介电参数的不确定性。结合MoM原理，右边向量与介电参数无关，将介电参数向量${\boldsymbol{\varepsilon }}$带入矩阵方程后，可以得到

式中，${\boldsymbol{Z}}({\boldsymbol{\varepsilon}} )$和${\boldsymbol{x}}({\boldsymbol{\varepsilon}} )$分别为关于介电参数${\boldsymbol{\varepsilon}}$的阻抗矩阵和电流。

3   基于泰勒级数展开的不确定性目标电磁散射特性分析方法

3.1   外形不确定目标电磁散射特性分析技术

把相应的外形向量元素${\alpha _i} \in [{\overline \alpha _i},{\underline {\alpha } _i}]$$(i = 1, \cdots ,t)$在一个区间内变化，然后通过三角基函数引入到积分方程公式当中，可以相应地表示目标外形不确定变化，从而建立带有外形向量${\boldsymbol{\alpha }}$的矩阵方程，$\alpha _i^{\mathrm{c}}$和$\Delta {\alpha _i}$分别定义为取值范围的中值和随机变化量，其表达式如下：

因此$\left[ {{{\overline \alpha }_i},{{\underline {\alpha } }_i}} \right]$可表示为$[\alpha _i^{\mathrm{c}} - \Delta {\alpha _i},\alpha _i^{\mathrm{c}} + \Delta {\alpha _i}]$。对于所有的随机变量，其中值向量${{\boldsymbol{\alpha}} ^{\mathrm{c}}} = [\alpha _1^{\mathrm{c}}, \alpha _2^{\mathrm{c}}, \cdots ,\alpha _t^{\mathrm{c}}]$和变化量向量$\Delta {\boldsymbol{\alpha}} = \left[ {\Delta {\alpha _1},\Delta {\alpha _2} ,\cdots ,\Delta {\alpha _t}} \right]$。

3.1.1   一阶泰勒级数展开

MoM积分方程中阻抗矩阵以及右边向量在中值${{\boldsymbol{\alpha}} ^{\mathrm{c}}}$处的一阶泰勒级数展开^[36-38]如下所示：

3.1.2   二阶泰勒级数展开

MoM积分方程中阻抗矩阵以及右边向量在中值${{\boldsymbol{\alpha}} ^{\mathrm{c}}}$处的二阶泰勒级数展开如下所示：

式中，阻抗矩阵、右边向量关于外形向量${\boldsymbol{\alpha}}$的导数计算过程见附录。因此，公式(7)的矩阵方程可以表示为

式中：${\boldsymbol{Z}}\left( {{{\boldsymbol{\alpha}} ^{\mathrm{c}}}} \right)$、${\boldsymbol{b}}\left( {{{\boldsymbol{\alpha}} ^{\mathrm{c}}}} \right)$和${{\boldsymbol{x}}^{\mathrm{c}}}$对应着外形向量${\boldsymbol{\alpha}}$取中值时的阻抗矩阵、右边向量以及感应电流系数，满足

$\Delta {\boldsymbol{x}}$为待求电流系数的扰动半径，通过推导公式(15)可以得到

对于一阶泰勒级数展开，$\Delta {{\boldsymbol{x}}_\alpha }$的表达式为

对于二阶泰勒级数展开，$\Delta {{\boldsymbol{x}}_\alpha }$的表达式为

通过电流系数的扰动半径$\Delta {{\boldsymbol{x}}_\alpha }$可以得到因目标外形变化导致的电流系数变化，进而可以计算出因目标外形变化对应的电磁散射特性，有

3.2   介电参数不确定目标电磁散射特性技术

将介电参数向量${\boldsymbol{\varepsilon}} = \left[ {{\varepsilon _1},{\varepsilon _2}, \cdots ,{\varepsilon _l}} \right]$引入到MoM矩阵方程中。对于目标材料介电参数的不确定性，可以表示为相应的随机变量${\varepsilon _i} (i = 1, \cdots , l)$在一个区间$[ {{{\overline \varepsilon }_i},{{\underline{\varepsilon} }_i}} ]$内的随机取值。$\varepsilon _i^{\mathrm{c}}$和$\Delta \varepsilon _i^{}$分别定义为取值范围的中值及其随机变化量，其表达式如下所示：

因此，$[ {{{\overline \varepsilon }_i},{{\underline{\varepsilon } }_i}} ]$可表示为$\left[ {\varepsilon _i^{\mathrm{c}} - \Delta {\varepsilon _i},\varepsilon _i^{\mathrm{c}} + \Delta {\varepsilon _i}} \right]$。对于所有的随机变量，其中值向量${{\boldsymbol{\varepsilon}} ^{\mathrm{c}}} = \left[ {\varepsilon _1^{\mathrm{c}},\varepsilon _2^{\mathrm{c}}, \cdots ,\varepsilon _l^{\mathrm{c}}} \right]$和变化量向量$\Delta {\boldsymbol{\varepsilon}} = \left[ {\Delta {\varepsilon _1},\Delta {\varepsilon _2}, \cdots ,\Delta {\varepsilon _l}} \right]$。

3.2.1   一阶泰勒级数展开

根据扰动法的原理，将阻抗矩阵在${{\boldsymbol{\varepsilon}} ^{\mathrm{c}}}$处用一阶泰勒级数展开，如下所示：

3.2.2   二阶泰勒级数展开

MoM积分方程中阻抗矩阵在${{\boldsymbol{\varepsilon}} ^{\mathrm{c}}}$用二阶泰勒级数展开，如下所示：

则矩阵方程可以表示为如下形式：

式中，${\boldsymbol{Z}}\left( {{{\boldsymbol{\varepsilon}} ^{\mathrm{c}}}} \right)$、${\boldsymbol{b}}$和${{\boldsymbol{x}}^{\mathrm{c}}}$为介电参数取中值时所对应的阻抗矩阵、右边向量以及电流系数。通过式(25)可以得到因目标材料属性变化带来的电流系数的变化量$\Delta {{\boldsymbol{x}}_\varepsilon }$：

对于一阶泰勒级数展开，$\Delta {{\boldsymbol{x}}_\varepsilon }$可以表示为

对于二阶泰勒级数展开，$\Delta {{\boldsymbol{x}}_\varepsilon }$可以表示为

通过电流系数的扰动半径$\Delta {{\boldsymbol{x}}_\varepsilon }$可以得到因材料介电参数不确定性导致的电流系数的变化，进而可以计算出因材料介电参数变化对应的电磁散射特性，有

叠加由于目标外形变化导致的电流系数的变化量$\Delta {{\boldsymbol{x}}_\alpha }$及由于目标材料属性变化产生的电流系数的变化量$\Delta {{\boldsymbol{x}}_\varepsilon }$，就可以计算出目标外形和介电参数同时变化对应的电磁散射特性，有

另外，基于MoM的计算平台，泰勒级数展开扰动法的计算复杂度和内存复杂度均为$O{(}{{N}^\dagger }^{2}{)}$(${{N}^\dagger }$表示未知量的个数)，应用MLFMA的MoM平台，其内存复杂度和计算复杂度均为$O\text{(}{{N}}^{\text{†}}\cdot \text{lg}{{N}}^{\text{†}}\text{)}$。

4   基于AWE技术的新型不确定性目标电磁散射特性分析方法

不同于AWE技术在角度、频率等标量场中的应用，针对具有可变几何形状的不可穿透非合作目标，本节将AWE技术扩展到矢量场，将雷达散射截面积(radar cross section, RCS)的高效计算求解器与AWE技术相结合，并将CFIE视为形状矢量的函数。为了简化积分方程中外形向量的导数计算，首先将PSDM推广到AWE中。随后，得到Padé 有理函数的系数，整个过程中，Padé 有理函数的系数计算只需进行一次。最后，在不同外形向量下分别计算双站RCS。对于电大规模目标，MLFMA可加快计算速度。结合AWE技术提出的新型不确定方法通过数值实例进行验证，证明了其准确性。与基于泰勒级数展开的传统扰动方法相比，该方法可以获得更大的扰动半径，且比传统的MC方法效率更高。

4.1   外形不确定目标电磁散射特性分析技术

4.1.1   外形不确定目标散射特性计算方法

首先，给出CFIE的表达式：

式中：$\alpha$为CFIE的组合因子，其值为[0,1]；${\boldsymbol{\alpha}}$为任意情况下的外形向量，${\boldsymbol{\alpha}} = [{\alpha _{\text{1}}},\cdots ,{\alpha _s},\cdots ,{\alpha _n}]$。将CFIE中的阻抗矩阵、右边向量展开为泰勒级数的形式，有

式中：L和M分别为Padé近似分子、分母的阶次，$L + M$为泰勒级数的阶次；${{{\boldsymbol{Z}}}^{(n)}}({{{\boldsymbol{\alpha}} }^{\mathrm{c}}})$为阻抗矩阵关于外形向量的n阶导数；${{{\boldsymbol{b}}}^{(n)}}({{{\boldsymbol{\alpha}} }^{\mathrm{c}}})$为右边向量关于外形向量的n阶导数。相应地，电流可以表示为矩矢量${{\boldsymbol{m}}_n}$与外形变化量的乘积形式：

将公式(32)～(34)带入公式(31)中，通过匹配等号两侧的同阶次项，有

能够得到矩矢量${{\boldsymbol{m}}_n}$的各阶表达式为：

实际编程计算过程中，电流${I}({{\boldsymbol{\alpha }}})$和矩向量${{\boldsymbol{m}}_{{n}}}$均为一维数组的形式。其中电流可表示为$\left[ {{{I}}_1}, {{{I}}_2},\cdots , {{{I}}_{{e}}},\cdots, {{{I}}_{{{N}^\dagger }}} \right]$，矩向量可表示为$\left[ {{m_{n1}},{m_{n2}},\cdots ,{m_{n{{e}}}},\cdots ,{m_{n{{N}^\dagger }}}} \right]$。为了进一步提高扰动法准确计算的扰动区间，利用Padé近似替换泰勒级数，有

式中：${a_i}_{{e}}$和${b_j}_{{e}}$分别为系数向量${{\boldsymbol{a}}_i}$和${{\boldsymbol{b}}_j}$的元素，即$\left[ {{a_i}_1,{a_i}_2,\cdots ,{a_i}_{{e}},\cdots ,{a_{i{{N}^\dagger }}}} \right]$，$[ {{b_j}_1,{b_j}_2,\cdots ,{b_j}_{{e}},\cdots ,{b_{j{{N}^\dagger }}}} ]$。此处不再赘述${a_i}_{{e}}$和${b_j}_{{e}}$的推导过程，其思路依旧为匹配等式两侧的同幂次项，对公式(38)进行化简得到：

通过上述推导，能够求解出公式(38)中的最终电流，进而获得外形不确定目标的电磁散射特性。

4.1.2   外形向量的各阶导数求解

不难发现，阻抗矩阵和右侧向量中包含大量与外形向量相关的量。若使用数学公式直接计算外形向量的高阶导数会导致计算负担随着控制点的增加而显著增加，为了克服这一问题，我们采用基于插值概念的PSDM进行计算。通常将外形向量的变化量表示为$\Delta {{\boldsymbol{\alpha}} } = [\Delta {\alpha _1},\cdots ,\Delta {\alpha _s},\cdots ,\Delta {\alpha _n}]$。$[{{{\boldsymbol{\alpha}} }^{\mathrm{c}}} - \Delta {{\boldsymbol{\alpha}} },{{{\boldsymbol{\alpha}} }^{\mathrm{c}}} + \Delta {{\boldsymbol{\alpha}} }]$表示外形向量的实际变化区间，在该区间中选取$N + 1$个插值向量${{{\boldsymbol{\alpha}} }_j}$，每个向量${{{\boldsymbol{\alpha}} }_j}$表示一个模型的外形控制点坐标信息。通过CGL(Chebyshev-Gauss-Lobatto)方法选取$N + 1$个插值向量，有

${{{\boldsymbol{\alpha}} }_j}$的元素${\alpha _{js}}$可以通过下式计算：

在PSDM中，通过构造微分矩阵${\boldsymbol{D}}$实现对外形向量的求导，其元素表达式为

为了统一微分矩阵${{\boldsymbol{D}}}$中的下标，对${\alpha _{js}}$进行转化，即对${\alpha _{is}}$而言，$- {N \mathord{\left/ {\vphantom {N {\text{2}}}} \right. } {\text{2}}} \leqslant i \leqslant {N \mathord{\left/ {\vphantom {N {\text{2}}}} \right. } {\text{2}}}$，而对于$\tilde{\alpha}_{js}$而言，${\text{0}} \leqslant j \leqslant N$，表达式为

以右边向量为例，计算其关于外形向量的一阶导数${{{\boldsymbol{b}}}^{({\text{1}})}}({{{\boldsymbol{\alpha}} }^{\mathrm{c}}})$，有

微分矩阵${{\boldsymbol{D}}}$的维度是$(N + {\text{1}}) \times (N + {\text{1}})$，$\left[ {{{\boldsymbol{b}}}({{{\boldsymbol{\alpha}} }_j})} \right]$是维度为$(N + {\text{1}}) \times {{N}^\dagger }$的矩阵。求得$\left[ {{{{\boldsymbol{b}}}^{({\text{1}})}}({{{\boldsymbol{\alpha}} }_j})} \right]$之后，$\left[ {{{{\boldsymbol{b}}}^{({\text{2}})}}({{{\boldsymbol{\alpha}} }_j})} \right]$可以进一步计算得到：

进而右边向量关于外形向量的n阶导数可以通过下式获得：

同理，矩阵矢量乘关于外形向量的各阶导数，也可以通过上述方式获得：

通过式(51)发现，并未直接求取阻抗矩阵关于外形向量的导数，原因有二：其一，存储阻抗矩阵的各阶导数会消耗大量内存，而存储矩阵矢量乘则能节省大量内存空间；其二，在使用PSDM求导的过程中，加法定理依旧适用，因此可以通过MLFMA方法计算。但是运用MLFMA计算过程中，阻抗矩阵的获取较为繁琐，而计算矩阵矢量乘是更为简洁的方式。因此，基于AWE技术的外形不确定性分析方法能够高效、准确地分析目标外形变化时的散射特性。

4.2   介电参数不确定目标电磁散射特性分析技术

4.2.1   介电参数不确定目标散射特性计算方法

首先给出关于介电参数的矩阵方程：

以多介质PMCHWT平台为例，其阻抗矩阵、电流可表示为如下关于介电参数的形式：

式中：${{{\boldsymbol{J}}}^i}({\boldsymbol{\varepsilon}} )$为等效电流；${{{\boldsymbol{M}}}^i}({\boldsymbol{\varepsilon}} )$为等效磁流。

通过将阻抗矩阵、电流展开为关于介电参数的泰勒级数形式，并带入公式(52)中，可得到各阶矩向量的表达式：

利用Padé近似有理多项式代替泰勒级数展开，逐元素计算电流系数：

通过上述推导，能够求解出公式(57)中的最终电流，进而获得介电参数不确定目标的电磁散射特性。

4.2.2   介电参数的各阶导数求解

相较于矩阵方程对外形向量的求导，介质目标求解平台中，右边向量与介电参数无关，因此对介电参数的求导相对简单。根据CGL插值方法选取$N + 1$个介电参数向量${{\boldsymbol{\varepsilon}} _j}$如下：

式中，${\boldsymbol{\varepsilon}} = \left[ {{{\boldsymbol{\varepsilon}} ^{\mathrm{c}}} - \Delta {\boldsymbol{\varepsilon}} ,{{\boldsymbol{\varepsilon}} ^{\mathrm{c}}} + \Delta {\boldsymbol{\varepsilon}} } \right]$。此时，构造微分矩阵${\boldsymbol{D}}$如下：

类似地，计算矩阵矢量乘关于介电参数的导数：

进而，计算出矩阵矢量乘的高阶导数：

因此，基于AWE技术的介电参数不确定性分析方法能够高效、准确地分析目标介电参数变化时的散射特性。基于MoM的计算平台，应用AWE技术扰动法的计算复杂度和内存复杂度均为$O{(}{{N}^\dagger }^{2}{)}$，应用MLFMA之后，其内存复杂度和计算复杂度均为$O\text{(}{{N}}^{\text{†}}\cdot \text{lg}{{N}}^{\text{†}}\text{)}$。

5   智能化不确定度计算评估系统

5.1   置信度与不确定度评估方法

5.1.1   基于Kullbuck-Leibler距离加权的置信度计算方法

评估置信度水平的方法众多，如基于Birchfield距离、基于Kullbuck-Leibler距离、基于Hausdorff距离等^[5]。结合与不确定性问题的匹配程度，本文选择基于Kullbuck-Leibler距离加权的置信度水平评估模型^[39]。首先，给定扰动法的概率密度函数$p$和MC法的概率密度函数$q$，则对称Kullbuck-Leibler距离的计算式如下：

假设本文中两类相互独立模型各个角度下的RCS数据服从高斯分布，其概率密度函数$p = N({\mu _p},{\sigma _p})$和$q = N({\mu _q},{\sigma _q})$的Kullbuck-Leibler距离可以计算如下：

离散度${d_i}$定义如下：

尺度参数$\rho$选取为全体分布的平均离散度，即各分布之间的平均距离表达式为

将上述两式代入求权函数的式子中，并对其进行归一化，可以得到最终的权函数如下：

最后将每个角度下的相似度权函数和该角度下概率值相乘以后累加得到所有角度下的一个联合概率值，即两组数据之间的相似度表达式如下：

对于目标的雷达散射特性，依据扰动法置信度水平，当其满足一定的阀值以后再去计算其不确定度。不确定度表示RCS在一定范围内的抖动程度，满足一定置信度水平情况下计算出来的不确定度具有可靠性。

5.1.2   基于RMSE加权的不确定度评估方法

不确定度的评估方法包括基于方差、基于均值方差加权、基于均方根误差(root mean square error, RMSE)加权的三种方法。综合三种方法的性能，选择基于RMSE加权的方法。首先，构造权函数如下：

式中：$\alpha < 1$、$\; \beta < 0$为引入的尺度变换参数；${R_i}$为第$i$个角度的RMSE，给定测试集$({X^v},{Y^v})$有

式中：${N_v}$为样本集中数据样本个数；${\hat y _i}$为扰动法单个角度所有样本的值；$\overline y$为MC法单个角度所有样本的值。每个角度所有样本下的RMSE与其在加权模型中所占的比重呈反比，即RMSE越小，该角度所占的权重越大。

由于上述每个角度的RCS数据互不相关，合成的标准不确定度计算如下：

结合RMSE权重的不确定度模型能更真实地反映和真实模型之间的偏差，该种方法比不加权重直接计算方差的方法更加准确。

5.2   基于BP神经网络的智能化不确定度评估系统

结合计算置信度和不确定度的计算方法，利用BP神经网络构建了一个智能化RCS预测和评估系统。以目标的外形存在不确定性为例，首先，分别通过扰动法和MC的程序计算出模型改变1000次时的RCS结果，再结合引入随机变量的外形信息，分别构建训练样本集。其次，根据样本数据训练BP神经网络，合理设置训练步数和权重系数，当loss满足要求后结束训练。选取部分数据进行测试，并计算真实值与预测值的RMSE，当RMSE达到一定阈值时，得到满足要求的测试结果。最后，计算出预测值和真实值的均值、方差等数据，并用于单独计算二者的置信度和不确定度。通过对比预测值和真实值的置信度与不确定度，说明智能化RCS预测与评估系统的有效性。相较于数值方法，该系统最大的优势在于保证精度的同时大大减少了计算时间。程序框架如图3所示。

具体地，在训练BP神经网络时，分别针对扰动法和MC的样本数据优化各自的权重系数和偏置值，得到各自的RCS预测值用于计算均值和方差，如图4所示。BP网络包括4层隐藏层，每层有100个节点，训练时间为11.5 min，将置信度和不确定评估算法集成至该智能化系统中。根据全波方法生成的样本数据，先计算出基于扰动法的真实值的均值和方差以及基于MC方法的真实值的均值和方差共4组数据；再根据网络测试数据，计算出基于扰动法的预测值以及基于MC方法的预测值的均值和方差共4组数据；最后利用这8组数据计算出基于真实值的扰动法与MC方法的置信度以及基于预测值的扰动法与MC方法的置信度。当置信度满足要求之后，计算不确定度。

图  3  智能化RCS预测和评估系统

Fig.  3  Intelligent RCS prediction and evaluation system

下载: 全尺寸图片 幻灯片

图  4  BP神经网络示意图

Fig.  4  Schematic diagram of BP neural network

下载: 全尺寸图片 幻灯片

6   算例分析

为验证本文提出方法的准确性，分别利用基于泰勒级数展开的扰动法、基于AWE技术的扰动法分析外形不确定目标、介电参数不确定目标的电磁散射特性，并通过实测数据验证了方法的有效性。另外，利用立方体算例，对搭建的智能网络进行验证。以下算例分析中，均将MC方法作为标准值供其他方法对比。其中，MC方法随机采样1 000次，须进行1 000次的阻抗矩阵填充，1000次迭代求解。基于泰勒级数展开的扰动法，阻抗矩阵的求导次数以及迭代求解次数与引入随机变量的个数一致。基于AWE技术扰动法须进行$(N + 2) \times (L + M) + 1$次阻抗矩阵填充，$L + M + 1$次迭代求解，在本文中，$N$等于$L + M + 2$，其中$L = M = 3$^[36]。另外，计算了不同角度下的方差线，来表明因目标外形变化导致的RCS波动情况。如果方差的绝对值大，说明在该角度下外形变化对RCS影响较大，反之RCS影响较小。引入随机变量时，随机数在扰动区间中满足均匀分布。

6.1   外形不确定目标电磁散射特性计算

基于CFIE平台分析外形不确定的金属立方体目标，利用NURBS技术构建立方体模型，如图5所示，该模型由8个控制点、6个NURBS组成，模型边长为2 m，未知量为7 200，本算例未使用MLFMA算法。入射波频率为300 MHz，入射角${\theta }_{\text{i}}=0{\text{°}}、 {\phi }_{\text{i}}= 0{\text{°}}$，散射角${\varphi _{\mathrm{s}}} = 0{\text{°}}$, $- {\text{90}}{\text{°}} \leqslant {\theta _{\mathrm{s}}} \leqslant {\text{90}}{\text{°}}$。将2、4、5、6号控制点设置为随机变量，扰动半径为0.5 m。

图  5  外形不确定金属立方体目标

Fig.  5  Metallic cube target with uncertain shape

下载: 全尺寸图片 幻灯片

由图6(a)发现，AWE技术扰动法与MC方法吻合良好，而二阶泰勒级数展开法出现了较大误差，说明本文提出的基于AWE技术的新型不确定性目标分析方法有较高的精度。通过图6(a)中的方差线分布还发现，在30°、150°附近以及65°～115°范围内，外形变化对RCS影响较大。图6(b)RCS的概率密度函数(probability density function, PDF)分布体现了图6(a)中RCS均值在不同区间的分布情况，在区间 [−5, −2] 和 [7, 10] 中，PDF 分别为 13.18% 和 23.07%。这两个区间的权重明显大于其他区间的权重，表明模型的 RCS 值更多地分布在这两个区间，从图6(c)累积概率密度函数(cumulative probability density function, CDF)分布也可以得出同样的结论。在该算例中，PDF和CDF是目标外形随机变化 1000 次时 RCS 统计分布的特征，能反映外形不确定立方体目标的RCS 真实分布。这种统计分布分析对目标识别和隐身设计也有帮助。

图  6  外形不确定立方体目标散射特性

Fig.  6  Scattering characteristics of target with uncertain shape cube

下载: 全尺寸图片 幻灯片

将上述三种方法的内存需求和求解时间进行比较，结果如表1。可以看出：基于AWE技术的扰动法在内存消耗方面和MC方法相近，但是速度比MC方法提高了19.6倍；而基于二阶泰勒级数展开的方法不仅内存消耗大，求解时间也远大于AWE方法。

表  1  立方体目标内存需求与求解时间情况

Tab.  1  Cube target-memory and time consumption

方法	内存消耗/MB	求解时间/s
AWE-扰动法	888.0	2138.5
泰勒级数	5547.1	7250.7
MC方法	877.6	42055.8

下载: 导出CSV 

| 显示表格

6.2   介电参数不确定目标电磁散射特性计算

6.2.1   基于泰勒级数展开的介电参数不确定目标电磁散射特性分析

分析一个介电参数不确定金属-介质BoR导弹模型，尺寸参数如图7所示。金属部分采用线段剖分，介质部分采用三角面元网格剖分，未知量数量为75，未使用MLFMA。入射波频率为300 MHz，入射波为${\theta }_{\text{i}}=0{\text{°}}、{\varphi }_{\text{i}}=0{\text{°}}$，散射角为${\varphi _{\text{s}}} = 0{\text{°}}$,$0{\text{°}} \leqslant {\theta _{{\mathrm{s}}}} \leqslant 180{\text{°}}$。介质#1相对介电常数为(3，−0.6)，并引入0.1的随机变量。导弹模型RCS的均值和方差，如图8所示。

图  7  介电参数不确定导弹模型及其控制点信息

Fig.  7  Medium uncertain missile model and control point information

下载: 全尺寸图片 幻灯片

图  8  介电参数不确定导弹目标双站RCS均值和方差

Fig.  8  Bistatic RCS mean and variance of the missile model with medium uncertainty

下载: 全尺寸图片 幻灯片

通过图8发现：基于泰勒级数展开的计算结果与MC方法吻合良好，说明当引入的随机变量较小时，该方法能够有效分析介质参数不确定目标的电磁散射特性；在0°～180°范围内，方差线绝对值均较小，说明此时介电参数的变化对RCS影响较小。

如表2所示，基于泰勒级数展开方法相比于MC方法内存会稍高，这是因为基于泰勒级数展开方法除了需要存储中值矩阵${\boldsymbol{Z}}({{\boldsymbol{\varepsilon}} ^{\mathrm{c}}})$外，还需要存储随机变量所对应的导数矩阵，但是相比于MC方法将会提高近10倍的效率，证明了该方法的高效性。

表  2  导弹目标内存需求和求解时间

Tab.  2  Missile target-memory requirements and solution times

方法	内存消耗/MB	求解时间/s
泰勒级数	109.1	27.7
MC方法	83.0	303.0

下载: 导出CSV 

| 显示表格

6.2.2   基于AWE技术的介电参数不确定目标电磁散射特性分析

分析一个介电参数不确定非均匀介质的圆柱目标，模型如图9所示。

图  9  非均匀介质圆柱模型

Fig.  9  Cylindrical model of inhomogeneous medium

下载: 全尺寸图片 幻灯片

该圆柱模型由两种不同介质材料的圆柱#1和#2上下摆放构成，两个圆柱半径均为0.5 m，高度均为0.5 m。未知量为6528，未使用MLFMA。在此算例中，将两种介质的介电参数设置为随机变量。入射波频率为300 MHz，入射角为${\varphi _{\mathrm{i}}} = 0{\text{°}}$、${\theta _{\mathrm{i}}} = 0{\text{°}}$，散射角为${\varphi _{\mathrm{s}}} = 0{\text{°}}$、$0{\text{°}} \leqslant {\theta _{\mathrm{s}}} \leqslant 180{\text{°}}$。圆柱#1的相对介电参数为(3，0)，相对磁导率为(1，0)；圆柱#2的相对介电参数为(4，0)，相对磁导率为(1，0)，其相对介电常数的扰动半径为0.4。圆柱目标双站RCS的均值及统计变化结果如图10所示。在45°～180°范围内，随着介电参数的变化RCS出现波动，尤其是45°～90°角度范围，说明此时介电参数变化对RCS影响较大。

图  10  圆柱的双站RCS均值和方差

Fig.  10  Mean and variance of bistatic RCS for the cylindrical

下载: 全尺寸图片 幻灯片

由表3可以看出，基于AWE技术的扰动法内存需求稍微大于MC方法，求解时间比采样1000次的MC方法提高了近30倍的效率，体现了本文方法相比于MC方法的高效性。

表  3  圆柱目标内存需求和求解时间

Tab.  3  Cylindrical target-memory requirements and solution times

方法	内存消耗/MB	求解时间/s
AWE-扰动法	1692.7	1428
MC方法	1692.4	43446

下载: 导出CSV 

| 显示表格

6.3   实测数据对比

实际加工不同外形尺寸的金属圆柱目标，并测试其在8～18 GHz范围内的单站RCS，其中入射角度为${\theta _{{\mathrm{i}}}} = 90{\text{°}}$、${\varphi _{{\mathrm{i}}}} = 0{\text{°}}$，散射角度为${\theta _{\mathrm{s}}} = 90{\text{°}}$、${\varphi _{{\mathrm{s}}}} = - 30{\text{°}} \sim 30{\text{°}}$，扫描间隔为0.2°，极化方式包括HH极化和VV极化。不确定外形的圆柱模型包括改变圆柱的高度、改变圆柱的直径两种情况。圆柱目标的实测场景如图11所示，其中圆柱高度改变时，直径为120 mm，高度在270～290 mm变化；圆柱直径改变时，高度为300 mm，直径在108～116 mm变化。测试结果如图12和13所示。

图  11  外形不确定金属圆柱目标及实测场景

Fig.  11  Metallic cylindrical target with uncertain shape and the measured scenario

下载: 全尺寸图片 幻灯片

图  12  改变圆柱高度时RCS的均值与方差

Fig.  12  Mean and variance of RCS when varying cylindrical height

下载: 全尺寸图片 幻灯片

图  13  改变圆柱直径时RCS均值与方差

Fig.  13  Mean and variance of RCS when varying cylindrical diameter

下载: 全尺寸图片 幻灯片

可以看出：1)改变圆柱高度时，HH极化下，在−29°、−9°、11°等角度RCS出现较大波动；VV极化下，RCS整体趋势较为稳定。2)改变圆柱直径时，HH极化下，在−23°、−11°、−10°、12°、13°、14°、20°、25°等角度RCS出现较大波动；VV极化下，除在±12°略有波动外，RCS整体趋势较为稳定。

另外，针对上述均值方差信息计算了各种情况时的置信度以及重合度，如表4所示，置信度均大于95%，重合度均大于75%。

表  4  不确定外形圆柱目标置信度与重合度分析

Tab.  4  Confidence and coincidence analysis of uncertain shape cylindricity

类别	置信度/%	重合度/%
变高度-HH	98.17	86.56
变高度-VV	97.84	76.52
变直径-HH	98.10	81.74
变直径-VV	98.38	75.57

下载: 导出CSV 

| 显示表格

6.4   基于AGX边缘计算盒子的智能化不确定度评估系统

分析一个具有不确定外形的1 m金属立方体目标，利用NURBS建模方法建立该目标模型，模型如图14所示，该模型由6个NURBS、8个控制点组成。电磁波垂直入射，频率为300 MHz，入射角度为${\theta _{\mathrm{i}}} = 0{\text{°}}、{\varphi _{\mathrm{i}}} = 0{\text{°}}$，观察角度为${\varphi _0} = 0{\text{°}}$，$0{\text{°}} \leqslant {\theta _0} \leqslant 180{\text{°}}$。将1号控制点的x、y、z坐标引入随机变量，其扰动区间为[−0.2$\lambda$,0.2$\lambda$]。AGX集成硬件系统实现如图15所示。

图  14  边长1 m的金属立方体

Fig.  14  Metal cube with a side length of 1 m

下载: 全尺寸图片 幻灯片

图  15  AGX边缘计算盒子系统实物

Fig.  15  Physical picture of AGX edge computing box system

下载: 全尺寸图片 幻灯片

基于Kullbuck-Leibler距离加权的置信度评估方法分别计算了预测值和真实值的置信度，二者分别为82.9%和83.4%。此时，认为置信度水平是可靠性的，计算基于RMSE加权的不确定度，如表5所示。可以看出，不确定度的预测值和真实值基本一致，说明了智能化RCS评估系统的有效性。

表  5  智能化不确定度评估系统的预测值和真实值

Tab.  5  Intelligent uncertainty assessment: predicted and real value

方法	预测值/dBsm	真实值/dBsm
AWE-扰动法	0.32	0.32
MC方法	0.28	0.27

下载: 导出CSV 

| 显示表格

最后，将智能化RCS预测与评估系统的计算时间与分别使用MC方法和扰动法等数值方法的计算时间进行对比，相较于扰动法速度提升了147.3倍，相较于MC方法速度提升了1381.5倍。

7   结　论

本文所提方法能够实现不确定参数目标的散射特性准确、高效计算，并结合评估方法对均值、方差进行了统计分析，利用神经网络对不确定度的智能计算进行了初步探索。所提方法对研究信息参数不完备的非合作目标，尤其是计算外形、介电参数不确定目标的电磁散射特性具有推进作用。不同于传统的求导方法，本文基于PSDM和AWE技术实现各阶导数电流的快速计算，节省大量公式求导的时间。另外，通过实测数据对比，说明本文方法能够有效分析目标的不确定性问题，并大幅提高了计算效率。在接下来的工作中，针对扰动法，我们将不断探索更高效的计算方法，并将扰动法和神经网络的结合作为研究重点。

附 录A

1)关于外形的导数计算公式

对于电场积分方程(electrical field integral equation，EFIE)，${{\partial {\boldsymbol{Z}}\left( {{{\boldsymbol{\alpha}} ^{\rm{c}}}} \right)} \mathord{\left/ {\vphantom {{\partial Z\left( {{\alpha ^c}} \right)} {\partial {\alpha _i}}}} \right. } {\partial {\alpha _i}}}$和${{\partial {\boldsymbol{b}}\left( {{{\boldsymbol{\alpha}} ^{\mathrm{c}}}} \right)} \mathord{\left/ {\vphantom {{\partial b\left( {{\alpha ^c}} \right)} {\partial {\alpha _i}}}} \right. } {\partial {\alpha _i}}}$分别表示为：

式中：k为波数；$\eta$为自由空间的波阻抗；$o$和$s$分别为观察基函数和源基函数的高斯点个数；${W_o}\left( {{K_o}} \right)$和${W_s}\left( {{K_s}} \right)$为相应高斯点的权重。对于磁场积分方程(magnetic field integral equation，MFIE)，${{\partial {\boldsymbol{Z}}\left( {{{\boldsymbol{\alpha}} ^{\mathrm{c}}}} \right)} \mathord{\left/ {\vphantom {{\partial Z\left( {{\alpha ^c}} \right)} {\partial {\alpha _i}}}} \right. } {\partial {\alpha _i}}}$和${{\partial {\boldsymbol{b}}\left( {{{\boldsymbol{\alpha}} ^{\mathrm{c}}}} \right)} \mathord{\left/ {\vphantom {{\partial b\left( {{\alpha ^c}} \right)} {\partial {\alpha _i}}}} \right. } {\partial {\alpha _i}}}$分别表示为：

2)关于介电参数的导数计算公式

以VSIE为例推导阻抗矩阵关于介电参数的导数，首先给出VSIE的阻抗矩阵表达式：

式中：${{\boldsymbol{Z}}}_{mn}^{{\mathrm{DD}}}$和${{\boldsymbol{Z}}}_{mn}^{{\mathrm{DM}}}$为关于介电参数的函数，分别表示介质与介质之间的相互作用，以及介质对金属的作用；${{\boldsymbol{Z}}}_{mn}^{{\mathrm{MD}}}$和${{\boldsymbol{Z}}}_{mn}^{{\mathrm{MM}}}$与介电参数无关，分别表示金属对介质的作用以及金属与金属之间的作用。${{\boldsymbol{Z}}}_{mn}^{{\mathrm{DM}}}$和${{\boldsymbol{Z}}}_{mn}^{{\mathrm{DD}}}$关于介电参数的导数表达式如下：

式中：$\omega$为角频率；${\varepsilon _0}$为真空中的介电参数；$\mu$为磁导率；${{\boldsymbol{f}}}_m^V({{\boldsymbol{r}}})$和${{\boldsymbol{f}}}_n^V({{\boldsymbol{r}}}')$为SWG基函数；$G({{\boldsymbol{r}},{\boldsymbol{r}}'})$为格林函数；${K_n}({{{\boldsymbol{\varepsilon}} }^{\mathrm{c}}}) = {{\left( {{{{\boldsymbol{\varepsilon}} }^{\mathrm{c}}}({\boldsymbol{r}}) - {\varepsilon _0}} \right)} \mathord{\left/ {\vphantom {{\left( {{{\varepsilon }^c}(r) - {\varepsilon _0}} \right)} {{{\varepsilon }^c}(r)}}} \right. } {{{{\boldsymbol{\varepsilon}} }^{\mathrm{c}}}({\boldsymbol{r}})}}$。

3)Padé近似系数推导过程

首先，将公式(38)进行化简得到：

进一步地，公式(A.8)可以展开为

通过匹配公式(A.9)两侧关于$({\left. {{\alpha _s}} \right|_{\boldsymbol{\alpha}} } - {\left. {{\alpha _s}} \right|_{{{\boldsymbol{\alpha}} _0}}})$的同幂次项，可得：

运用递归的思想，可推导得公式(A.11)，即正文的公式(39) 。需注意的是在公式(A.11)中${{\boldsymbol{b}}_{\text{0}}} \equiv {\text{1}}$，

考虑到泰勒级数的最高展开阶次为$L + M$，公式(A.9)可重写为

匹配公式(A.12)中的同幂次项后，发现(A.12)左侧存在许多零项。整理如下：

为便于表达，我们将公式(A.13)统一为

进而，${b_{{{je}}}}$可以化简为

至此，正文部分公式(38)～(40)的关系推导完毕。