周琳，黄江涛，高正红,*

1. 西北工业大学 航空学院，西安 710072 2. 中国空气动力研究与发展中心，绵阳 621000

雷达散射截面(Radar Cross Section, RCS)反映了物体在给定方向上对入射雷达波散射的强弱，是衡量飞机隐身性能的重要指标,考虑隐身的飞行器设计常以减小雷达散射截面作为隐身设计的重要目标。现有飞行器隐身设计中多采用几何光学法(GO)、物理光学法(PO)、几何绕射理论(GTD)、物理绕射理论(PDT)等高频近似算法评估散射体的雷达散射截面[1-2]。高频算法根据高频场的局部性原理，仅根据入射场独立地近似确定表面感应电流[3]，计算速度快，所需内存小。近年来，随着反隐身技术的发展，尤其是全数字化有源相控技术在低空警戒和监视雷达中的应用，对飞行器隐身性能的要求越来越高[4]。高频方法虽然可以预估飞行器目标的镜面散射、边缘绕射等散射特性，但由于算法本身的局限，忽略了一些关键部件间的耦合关系，在弱散射源计算中的精度不足，使其在低RCS构型的评估和优化过程中可信度较低，无法满足低RCS构型精细化设计的需求[5-6]。基于严格理论模型的数值解法如矩量法(MOM)及基于矩量法的快速解法(如多层快速多极子算法(MLFMA))[7-8]从电磁场积分方程出发，可以精确求解三维复杂外形目标的电磁散射。随着高性能计算技术和低可探测飞行器的发展，基于精确建模的积分方程数值解法逐渐成为飞行器隐身设计中重要的电磁分析手段。

飞行器隐身性能与其外形密切相关，设计中常需处理隐身指标与气动指标之间的矛盾[9]。现有的气动隐身一体化设计多采用粒子群算法、遗传算法、神经网络算法等智能搜索算法[1, 6, 9-10]。智能搜索算法具有收敛到全局最优的能力，但优化效率较低，优化所需的雷达散射截面评估次数随设计问题规模的增加迅速增加，与高精度电磁模拟方法结合时计算成本昂贵。基于梯度的优化算法效率较高[11]，但经典基于有限差分的RCS梯度计算代价大，限制了梯度算法在气动/隐身优化中的应用，目前对基于梯度的气动/隐身一体化优化研究较少。基于伴随思想的梯度求解方法可以通过一次原方程求解和一次伴随方程求解得到目标关于所有设计变量的导数，目前已在气动外形优化[12]、气动/结构耦合优化[13]、声爆优化[14]等领域得到广泛应用。基于伴随方法的优化在隐身领域研究较少，Bondeson等[15]推导了二维积分形式的麦克斯韦连续伴随方程，并对二维简单外形(圆形)的RCS进行了优化；Wang和Anderson[16]推导了二维时域麦克斯韦方程的伴随方程，通过给定翼型的目标电流分布，对翼型的几何进行反设计；本文作者[17]在三维矩量法的基础上推导了麦克斯韦方程的离散伴随方程，验证了伴随方法的精度，但仍面临矩量法求解计算成本高、工程应用困难的问题。

针对三维问题中高精度RCS评估梯度计算成本高的问题，本文引入麦克斯韦积分方程的伴随方程，以三维矩量法为例详细推导了麦克斯韦积分方程的离散伴随方程和雷达散射截面的变分表达式。针对矩量法在三维复杂问题中内存、计算量大的问题，采用多层快速多极子算法(MLFMA)进行求解。选取双椎体模型、导弹模型对求解器的计算精度、伴随方法的可靠性进行验证。计算结果表明本文采用的求解器有较高精度，基于矩量法、多层快速多极子算法伴随梯度与差分梯度吻合良好，实现了通过两次线性方程组求解获得雷达散射截面关于所有设计变量的导数。证明了伴随方法的可靠性和高效性，为基于梯度的高精度气动/隐身一体化优化奠定基础。

1 积分方程数值求解方法

本文采用矩量法和多层快速多极子算法求解麦克斯韦积分方程及其离散伴随方程。矩量法是文献[7]提出的一种用于严格计算电磁问题的数值方法。矩量法将积分形式的麦克斯韦方程离散后转化为矩阵方程，通过求解线性方程组来得到目标表面感应电流，进而计算散射场。快速多极子算法(FMM)[8]和多层快速多极子算法[18]是建立在矩量法和线性方程组迭代求解算法基础上的快速算法。

1.1 基于RWG基函数的矩量法

RWG(Rao-Wilton-Glisson)基函数[19]是当前使用较为广泛的一种矩量法基函数，该基函数定义在具有公共边的相邻三角形上，可模拟任意形状物体的表面电、磁流分布。对于RWG基函数矩量法，三角形剖分的网格边长一般在[λ/12，λ/8]之间较为合适[20](λ为入射波长)。根据麦克斯韦方程和导体表面S上切向电场连续条件可得电场积分方程为

(1)

(2)

(3)

(4)

图1 RWG基函数示意图[21]Fig.1 Sketch of RWG basis function[21]

(5)

式中：In为第n个基函数的系数；N为三角形公共边总数。采用RWG基函数作为检验函数(伽略金法)，检测电场积分方程可得

(6)

整理写成矩阵形式：

ZN×NIN×1=VN×1

(7)

式中：Z为矩量法阻抗矩阵；I和V分别为电流系数向量和激励向量。阻抗元素和激励项的表达式为[7, 22]

(8)

(9)

1.2 多层快速多极子算法

矩量法离散得到的阻抗矩阵是复数稠密矩阵，如果用直接法求解矩阵方程，则算法存储量为O(e2)，计算量为O(e3)(相对于矩阵求逆运算量级)；如果用迭代法求解，存储量为O(e2)，每步迭代的计算量为O(e2)(相对于矩阵矢量乘法运算量级)，其中e为未知量的个数。当频率增加时，未知量个数迅速增加，对内存和计算速度提出了很高要求，限制了矩量法在高频散射问题中的应用。为提高矩量法求解问题的规模，在矩量法和线性方程组迭代求解算法的基础上发展了快速多极子算法和多层快速多极子算法，将计算量和存储量进一步降低到O(e1.5)和O(elge)量级，使基于精确建模的积分方程数值解法求解电大尺寸目标的电磁散射问题成为可能[23]。

快速多极子算法把基函数分为G组，将组之间的作用分为近相互作用和远相互作用，即将阻抗矩阵表示为

Z=Znear+Zfar

(10)

对于近相互作用，仍采用矩量法计算；对于远相互作用，则采用快速多极子算法。快速多极子算法利用加法定理和平面波展开定理将格林函数展开为多极子形式，进而将矩阵与向量的乘积运算转化为聚合-转移-配置3个过程。标量格林函数G(r′,r)=e-jkR/R的多极子展开形式为

(11)

(12)

(13)

式中：

(14)

(15)

其中：ji和tj分别为电流J(r)的基函数和测试函数;右端第1项表示邻近组的贡献；右端第2项为各远距离组的贡献，通过式(14)、式(12)和式(15) 表示的聚合、转移和配置3个步骤得到。

多层快速多极子算法是快速多极子算法的多层扩展。多层快速多极子算法基于树形结构，通过在多个层级上分组，组间嵌套，逐层递推来实现FMM[23]，进一步降低了内存需求，提高了矢量矩阵相乘的计算效率。由于FMM和MLFMA本质上加速的是线性方程求解过程中矩阵和向量的乘积运算，在计算过程中并不显式存储远相互作用矩阵。

2 基于伴随方程的梯度计算

隐身优化设计中常需要雷达散射截面关于设计变量的梯度信息。传统的有限差分法计算梯度所需雷达散射截面的评估次数与设计变量个数成正比，当设计变量较多时，有限差分法的计算代价较高，难以在实际问题中应用。基于伴随方程的梯度计算可以通过一次雷达散射截面求解、一次伴随方程求解获得雷达散射截面关于所有设计变量的梯度，有效提高梯度计算效率。本节在1.1节介绍的矩量法基础上对离散伴随方程和雷达散射截面的变分形式进行推导[17]，并对伴随方法所需的计算量和存储量进行简要分析。

2.1 伴随方程推导

典型的雷达散射截面优化设计问题形式为

(16)

需要满足的约束为

R(I(D),D)=V-ZI=0

(17)

式中：σ为雷达散射截面；I为感应电流；D为设计变量。雷达散射截面是几何外形和感应电流的函数，在矩量法求解中阻抗矩阵、激励和解得的感应电流均只由几何形状决定，即

(18)

引入拉格朗日乘子Λ可以构造目标函数：

L=σ+ΛTR

(19)

对式(19)进行求导，有

(20)

从式(20)可以看出，若找到合适的Λ使式(20) 右端第1项为0，可完全消除感应电流对设计变量导数dI/dD的计算，大幅度降低计算量，即

(21)

注意到：

(22)

将式(22)代入式(21)即可得到伴随方程：

(23)

伴随方程的形式与正计算(RCS评估)的形式一致，可以直接采用正计算采用的数值方法(矩量法、多层快速多极子算法)求解。将式(23)代入式(20)，基于伴随方法的目标梯度可以写为

(24)

(25)

2.2 雷达散射截面变分推导

伴随方程式的右端项激励∂σ/∂I为雷达散射截面关于感应电流的导数，以1.1节中介绍的基于RWG基函数的三维矩量法为例对伴随方程的激励项进行推导。

通过矩量法或多层快速多极子算法解得表面感应电流分布后，则可以计算雷达散射截面标量(m2)[24]，即

(26)

根据复数的运算性质将模的平方写成自身与其共轭相乘的形式：

(27)

(28)

式中:上划线表示共轭。由链式求导法则：

(29)

复数求导需分别对复数的实部和虚部求导，即

(30)

则式(29)可以整理为

(31)

式(28)中：g离散后可以写成感应电流和的形式，当采用1.1节介绍的RWG离散感应电流时，具体表达式为

(32)

(33)

将式(33)代入式(31)，有

(34a)

(34b)

将式(34)进一步代入式(30)即可得到伴随方程右端项的表达式为

(35)

2.3 伴随方法计算量和存储量分析

基于伴随方法的雷达散射截面导数求解主要由RCS正计算、伴随方程求解、导数计算3部分构成。

在RCS正计算中，矩量法需要一次阻抗矩阵填充和一次线性方程组求解；多层快速多极子算法需要多次计算ZI，计算次数与所需迭代次数成正比。当采用迭代法求解时，2种算法的计算量和存储量分别为O(e2)和O(elge)。在伴随方程求解中，若正计算的阻抗矩阵为对称阵，则矩量法可以直接使用正计算中的阻抗矩阵，不需要重新填充伴随方程的阻抗矩阵，可通过一次线性方程组求解完成伴随变量计算；多层快速多极子算法不显式存储阻抗矩阵，因此多层快速多极子算法在伴随计算阶段与正计算所需的计算量基本一致。

3 梯度验证

梯度验证采用基于电场积分方程的矩量法和快速多极子算法，线性方程组求解采用基于双共轭梯度算法[25]的迭代法求解。双共轭梯度算法较之普通的共轭梯度方法有较快收敛性，特别是在矩阵为对称阵的情况下，可以有效降低方程组求解的计算量[23]，计算过程中取收敛阈值ε=10-5。本节采用计算电磁学经典算例及实际工程中复杂外形对伴随梯度的可靠性和精度进行探讨。

3.1 双椎体模型

双椎体模型长7.5 in(1 in=25.4 mm)，由2个 椎角分别为22.62°和46.4°的半椎拼接而成[26-27]，是EMCC(ElectroMagnetic Code Consortium)[27]提供用于校验计算电磁学代码的标准算例之一。双椎体模型两端存在尖点，高频方法难以准确计算，需采用高精度数值求解方法。计算采用的雷达波频率f=9 GHz，单站散射水平极化和垂直极化，网格平均尺寸约为λ/10，未知量总数e=11 094。采用矩量法、多层快速多极子算法计算，并用商用软件和实验值(EXP)[26-27]校核计算结果。几种算法的计算结果如图2和图3所示。本文采用的矩量法和多层快速多极子算法均与实验值及商用软件FEKO矩量法的计算结果吻合良好，具有较高精度。

采用基于HMQ全局径向基函数的域元法[28-29]对双椎体模型进行参数化，设计变量分布如图4所示。采用基于矩量法、多层快速多极子算法的伴随方法和基于矩量法的有限差分法求解各设计变量在Z方向的梯度，认为基于矩量法的有限差分(FDD)结果为梯度的精确解。有限差分计算中扰动量Δx=10-3m。设计变量A和B(见图4)在各入射角度的导数计算结果如图5和图6所示。

图2 双椎体水平极化计算结果Fig.2 Numerical results of double-ogive horizontal polarization

基于伴随的矩量法和快速多极子算法计算得到的梯度与有限差分得到的梯度在趋势上吻合较好。其中，基于伴随的矩量法与有限差分的计算结果基本一致，具有较高精度；多层快速多极子算法得到的梯度与有限差分在梯度较小的角域存在一定误差，在梯度峰值附近吻合良好。入射角为0°时的感应电流模值分布和伴随电流模值分布如图7所示。值得注意的是，在流场伴随求解中，流场伴随变量传播方向与流场守恒变量相反[12]，而在电磁场求解中，电磁场伴随变量分布与感应电流分布较为一致。在求解过程中，共计算180个入射角度，其中矩量法计算使用20核并行计算，最大内存占用1.85 G，总计算耗时42 min。多层快速多极子算法计算使用8核并行计算，最大内存使用18.3 MB，计算用时3.7 min。

图3 双椎体垂直极化计算结果Fig.3 Numerical results of double-ogive vertical polarization

图4 双椎体设计变量分布Fig.4 Design variables distribution of double-ogive

图5 双椎体设计变量A导数(水平极化)Fig.5 Gradients of double-ogive design variable A(horizontal polarization)

图6 双椎体设计变量B导数(水平极化)Fig.6 Gradients of double-ogive design variable B(horizontal polarization)

图7 双椎体感应电流及伴随变量云图Fig.7 Contours of surface current and adjoint variable of double-ogive

3.2 导弹模型

导弹模型是宏观上的电大尺寸和细节上的电小尺寸共存的复杂模型，具有较高的实际工程意义，对电磁散射计算和梯度计算提出了较高要求。计算采用的导弹模型全长4 m，弹体直径0.28 m。计算采用雷达波段频率f=1 GHz，单站散射水平极化。网格平均尺寸约为λ/10，未知量e=26 157。采用矩量法、多层快速多极子算法计算，并与商用软件的矩量法进行对比。几种算法的计算结果对比如图8所示，矩量法和多层快速多极子算法的计算结果一致性较好，与商用软件的计算结果区别较小，具有较高的可信度。

采用基于HMQ全局径向基函数的域元法对导弹模型进行参数化，设计变量分布如图9所示。采用基于矩量法、多层快速多极子算法的伴随方法求解各设计变量在z方向的梯度并与基于矩量法的有限差分计算结果对比，计算中差分扰动量Δx=10-2m。其中设计变量C和D(如图9所示)在各入射角度的梯度计算结果如图10和图11 所示。基于伴随方法的矩量法和多层快速多极子算法得到的梯度与有限差分方法得到的梯度在趋势上较为一致。其中，基于伴随的矩量法与有限差分计算得到的梯度吻合较好，误差较小；基于多层快速多极子算法的伴随求解在梯度较小及梯度起伏较为剧烈的区域存在一定误差。由图10 和图11可以看到，雷达散射截面的导数随入射角度变化在大小和正负号上均存在明显起伏。图12为入射角为0°时表面感应电流和伴随变量的分布云图。

图8 导弹单站RCS结果(水平极化)Fig.8 RCS results of missile (horizontal polarization)

图9 导弹模型表面网格Fig.9 Surface mesh of missile model

图10 导弹设计变量C导数(水平极化)Fig.10 Gradients of missile design variable C(horizontal polarization)

图11 导弹设计变量D导数(水平极化)Fig.11 Gradients of missile design variable D(horizontal polarization)

图12 导弹感应电流及伴随变量云图Fig.12 Contours of surface current and adjoint variable of missile

在导数求解过程中，矩量法计算中内存使用10.4 G，计算使用20核并行计算，计算180个入射角度，总计算耗时608 min。多层快速多极子算法内存使用25.4 MB，计算使用8核并行计算，计算用时57 min。

4 结 论

本文提出了一种基于麦克斯韦积分方程离散伴随方程的RCS梯度计算方法，实现了RCS梯度的高效、高精度求解，为基于梯度的电大尺寸目标高精度气动/隐身一体化优化提供了基础。

1) 麦克斯韦积分方程的离散伴随方程形式与原积分方程一致，可直接采用原积分方程的数值解法如矩量法及多层快速多极子算法求解，程序实现容易。

2) 伴随方程的计算量与设计变量个数基本无关，可以实现通过2次线性方程组求解得到雷达散射截面关于所有设计变量的梯度信息。伴随求解的计算量与原方程求解基本一致，存储量在原方程求解的基础上增加不明显。

3) 基于伴随方程的梯度求解方法具有较高精度，多层快速多极子算法求取的梯度在精度上虽然略低于矩量法的计算结果，但在计算时间和内存需求上均明显优于矩量法，更适合于电大尺寸复杂问题的评估和优化。