1 引言

逆合成孔径雷达成像技术（Inverse Synthetic Aperture Radar，ISAR）对于环境中目标识别、隐藏目标非侵入性检测、雷达制导与导航等都有十分重要的应用.为了研究ISAR中目标散射特征与散射中心提取，从而重构目标特征性物理参数，ISAR模拟成像已成为一个十分重要的研究方向.但是，ISAR成像模拟需要先获得宽频带宽角度的目标散射场数据^[1~3]，这是复杂多目标ISAR成像模拟中最耗时的工作.

电大尺寸目标的电磁散射通常采用高频方法如物理光学（PO）、几何光学（GO），以及基于GO和PO的射线追踪模型等，往往考虑目标的平面近似模型，不适合求解小尺寸目标和具有曲面结构的绕射、爬行波等散射过程.经典的数值方法[4~6]虽然可以对研究目标模型进行准确的散射数值计算，但是它们对于目标模型的离散剖分往往导致庞大的存储量和运算量需求，尤其是成像技术中需要获得目标对不同频率不同入射角度的电磁波散射的大量计算.另一方面，传统的基于波函数解析理论往往只考虑单个体目标的情况^[7]，对于空间分布的多目标模型，尤其当目标尺寸较大或彼此相距较近时，各目标间的耦合散射对总散射场的贡献是难以计算的.近年来，针对多目标的散射有一些算法研究与应用，如用互易性原理分析细长导体圆柱的散射^[8]，但其散射的面积分计算仍是需要耗费大量的计算量.此外，还有采用广义散射矩阵和聚合T矩阵相结合的等效微波网络法计算偏心介质圆柱的散射^[9]，但它们形式往往相当复杂、计算也显得繁琐.

本文提出一种迭代的解析方法，基于柱面波函数展开和加法定理快速计算多导体圆柱的电磁散射，多导体柱间的耦合散射计算由迭代过程完成.对每个目标自身的一阶散射用各自中心坐标系的柱面波函数展开，而由各圆柱中心偏移产生的初始相位来进一步考虑各圆柱间的干涉作用.将每个圆柱目标的一阶散射看成另一圆柱目标的激励场，由柱面波函数加法定理和边界条件获得各目标间的二次耦合散射.类似地，可处理三阶和更高阶的散射计算.对场景中多导体圆柱目标完成计算宽频带和宽角度的散射场数据，然后经过频率向的平滑处理和f-θ平面的插值运算后，用FFT方法进行ISAR成像模拟，很好地重构了各圆柱体的位置和大小尺寸.

2 散射场的解析计算

考虑图 1中多个随机位置和尺寸的理想导体圆柱的散射，TM极化的平面电磁波沿=（cosφ_i，sinφ_i）方向入射，入射电场只有z极化分量，表示为

(1)

其中k=2 πf/c为入射电磁波的波数，f为电磁波频率，c为自由空间中的光速.图中=（cosφ_s，sinφ_s）表示散射方向.

由于二维导体柱不存在交叉极化效应，空间任意一点的散射场也只有z方向分量（设为），它包括每个圆柱目标自身独立的散射（一阶散射）和各圆柱之间的相互耦合散射（高阶散射），本文先阐述解析方法计算各阶散射场.

2.1 一阶散射计算

一阶散射即是不考虑各目标间的相互耦合作用而仅考虑各个圆柱自身的独立散射.设圆柱m的中心ο_m在主坐标系xoy中的位置为（x_m^c，y_m^c）（上标c表示圆柱目标的中心center），对应的圆柱坐标为（ρ_m^c，φ_m^c）.建立圆柱目标m的局部坐标系x_mο_my_m，在该坐标系中的任意位置矢量r_m表示为（x'_m，y '_m）或（ρ'_m，φ'_m）（上标'表示局部坐标系中的坐标）.

将主坐标系中的任意位置矢量r（x，y）转换到局部坐标系x_mο_my_m中，得到式的平面电磁波在该局部坐标系中的表达式为

(2)

与主坐标系中的激励场形式相比，局部坐标系中的激励场引入了一个初始相位，即图 1中ο点和ο_m点所在等相位面的相位延迟（图中粗实线表示入射电磁波的等相位面）.

在局部坐标系x_mο_my_m中，令E ₀=exp[-jkρ_m^c×cos（φ_i-φ_m^c）]，则圆柱m在式的平面波作用下的散射场可由柱面波函数展开方法计算为^[10]

(3)

其中a_m表示圆柱m的半径.对于ka_m < 1，上式求和的截断项数取L=3已有足够的精度^[11].可见它们是以ο_m为中心向外传播的柱面波.所有目标的一阶散射总场表示为各个离散圆柱散射场的叠加，即

(4)

式中各局部坐标分量与主坐标分量的关系为

(5)

对远区散射场，由坐标变量的远区近似式

(6)

以及Hankel函数的大宗量近似式^[12]

(7)

将式（6）和（7）代入式（4），得到远区散射场为

(8)

进一步交换求和顺序，得到

(9)

对后向散射，令φ_s=φ_i+π得到

(10)

可见对于一阶散射而言，各圆柱目标散射的相互作用主要体现在由于各自中心偏移而产生的相位干涉.

2.2 二阶散射计算

由Hankel函数的加法定理^[12]

(11)

其中φ ″为矢量ρ-ρ'与x轴的夹角.

考虑图 2中圆柱m对圆柱n作用的二次散射，即圆柱m的一阶散射场（即公式）对圆柱n的激励场，它可看成位于中心ο_m处不同种类偶极子对柱n的作用.当圆柱体尺寸与其间距离相比很小时，圆柱n对ο_m所张的角度很小，因而在圆柱n上任意点在圆柱m的局部坐标系中的相角φ'_m≈φ_mn+π（其中φ_mn表示圆柱m与圆柱n中心连线在主坐标系中的相角）.因此，激励偶极子可近似为

(12)

此时相位因子变成了常数，下面的二阶散射分析中暂且不考虑.

考虑偶极子E '_zi=H_l²（kρ '_m）对圆柱n的激励，在其局部坐标系x_nο_ny_n中展开为

(13)

设其散射场可展开为如下形式

(14)

由导体圆柱n表面的边界条件E'_zi+E'_zs|_ρ'n=an=0，得到展开系数为

(15)

对所有H_l²（kρ'_m）的激励求和，代入式得到圆柱m对圆柱n的二次散射为

(16)

式中各局部坐标分量见图 2.

圆柱n的二次散射场作用于另一圆柱的三次甚至更高阶的散射，可用类似的方法计算.后文的数值分析得知，当圆柱分布比较稀疏时，只需考虑到二阶散射即可.因此，空间任意一点的总散射场为一阶和二阶散射之和，双站雷达散射截面表示为

(17)

当φ_s=φ_i+π时即为后向散射σ_c.

3 散射的数值分析

考虑三个导体圆柱目标的散射，分别位于ο₁（-0.1，0）、ο₂（0，0.07）、ο₃（0.1，0），其半径相同均为r=0.01 m，模型见图 3中附图.平面电磁波频率f=2 GHz，入射角φ_i=0 °，散射角φ_s∈[0，360 °]共180个点.迭代解析计算中的级数求和截断取L=P=3.迭代解析计算中需要针对每个散射角进行柱面波级数求和运算，计算180个角度的双站散射总共耗时11.5s.数值MoM计算中每个圆柱表面围线等间隔地离散成30小段，共60个未知量.数值MoM方法的主要运算在于阻抗矩阵的填充和矩阵方程的求解，对于双站散射计算仅考虑一个入射角激励，只需执行一次MoM过程和求解一次矩阵方程，耗时1.98s.

图 3给出了双站雷达散射截面σ _bi随散射角φ_s的变化，图中虚线仅包含一阶散射，细实线表示一阶和二阶散射之和，粗实线表示MoM数值结果.可见二阶解析散射基本能与MoM数值结果吻合，而一阶散射与精确的MoM数值的计算结果相差较大，尤其是在前向散射方向.这是因为圆柱的电小尺寸导致电磁波在其上主要以爬行波绕射为主，圆柱之间的距离较近导致它们之间耦合的耦合散射不容忽略.

考虑后向散射，解析计算还是只需执行式（4）和式（16）的级数求和操作，计算180个方向的后向散射仍然只需耗时11.5s.而数值计算需要对每个入射角进行一次MoM运算，每次耗时1.98s.计算180个角度的后向散射需要求解180次矩阵方程，共耗时375.3s，是本文解析迭代方法的30多倍.图 4给出了后向雷达散射截面随入射角的变化关系，可以看出二阶散射能基本与MoM数值结果相吻合，三阶以上的散射贡献很小可忽略.

现考虑位置随机分布的10个导体圆柱目标，中心位置的x和y坐标在-0.5~0.5 m之间均匀分布，半径相同r=0.02 m，模型如图 5a所示.入射电磁波频率f=2 GHz.图 5b 给出了后向散射截面随入射角的变化关系，可以看出：本文的二阶散射解与MoM数值结果吻合相当好，自然一阶散射解有较大差别.本文的快速计算保证了能有效地进行宽频带宽角度散射计算，完成ISAR成像.

4 ISAR成像模拟

如图 6所示，接收雷达位置A点位于目标远区，则目标区域任一点（x，y）到接收A点的径向距离近似为

(18)

设在位置（x，y）处放置散射强度为1的点目标，则A点接收到的该点目标的散射场为

(19)

上式计及了单站雷达电磁波到目标的双程距离，波数k=2 πf/c.

设目标的反射率密度函数为g（x，y），则A点接收的总散射场为^[13]

(20)

式中的常系数exp（-j 4 πfρ₀/c）对成像结果的影响可不考虑（后文忽略该因子）.令

(21)

则（20）式可写成

(22)

可以看出它是二维Fourier逆变换形式.由Fourier变换关系得到

(23)

为了获得目标的重构图像g（x，y）的精确解，式（21）需在X-Y平面上进行二维无穷积分.由式（21）的对应关系可知，需要获得频率范围f ∈（-∞，∞）和相角范围θ ∈（0，2 π）的散射场数据G（X，Y），这自然在实际测量中是无法达到的.现实中往往考虑一定频带f ₁~f ₂范围内的散射场数据，即图 7中半径f ₁~f ₂之间的圆环内的数据，其他地方（图中白色区域）的散射场数据为零.

为了抑制散射场数据截断所引起的Fourier变换的振铃效应（ringing）^[13]，需对空间的散射场数据进行平滑窗处理，这里采用Hamming窗平滑函数，其形状如图 8所示，表达式为

(24)

其中f _1，2表示宽频带的两端点，Δf表示平滑宽度，它控制数据平滑过渡的快慢程度.经过平滑窗处理后，图 7中深灰色区域为原始散射场值（如图中A点），浅灰色区域为平滑数据（如图中B ₁、B ₂点），白色区域为零场值（如图中C ₁、C ₂点）.

为了充分利用快速Fourier变换（FFT）的高效计算方法，需要获得X-Y平面上均匀离散点（X_m，Y_n）上的散射场数据.若考虑对f-θ平面上均匀离散样点（f_i，θ_j）的散射场，则由（X_m，Y_n）对应的频率和相位θ _mn=arg（X_m，Y_n）往往就不在（f_i，θ_j）离散样点上，本文考虑由其周围四点的二维插值计算获得（这里采用三次样条插值），如图 7中A点的散射场由周围P₁，P₂，P₃，P₄四点的数据插值计算.

经平滑和插值后得到X-Y平面上均匀离散点（X_m，Y_n）的散射场数据G（X_m，Y_n），由二维的离散Fourier变换关系得到目标的近似重构图像为

(25)

上式可直接用FFT实现.

前面关于散射场解析计算的推导是考虑入射电磁场在主坐标系原点的相位为零，而ISAR成像原理则是考虑电磁波从雷达出发经目标上各点散射后回到雷达的双程传播过程，为此需将式中的入射电磁波作如下的相位变换：

(26)

经（26）式变换后的入射场在发射雷达位置的相位为零，式中距离ρ为雷达相对于目标区域基准点（全局坐标系原点）的距离.可以看出零相位面的改变只是对入射场引入一个常指数因子，对上面第2节推导的散射计算公式没有本质的影响，只需在散射场公式前乘上常系数exp（-jk ρ）即可.

对图 3中的圆柱目标模型，考虑频率f=1.5~2.5 GHz范围内41个离散频点，方位角θ=0~360 °范围内120个相位方向，用上面第2节描述的方法解析计算f-θ平面上均匀离散样点（f_i，θ_j）上的后向散射电场.对这些散射场数据用Hamming窗进行平滑处理，选择平滑宽度为8个离散频率点即Δf=0.2 GHz.然后经由二维的三次样条插值得到X-Y平面上均匀离散样点（X_m，Y_n）上的场值，最后用FFT算法实现式（25），即可得到目标的重构图像，如图 9所示.与图 3中的附图相比，重构图像相当准确地定位了各圆柱目标的位置和形状廓线.

同样地，对图 5 a中10个随机分布导体圆柱目标的散射场数据进行ISAR成像，如图 10所示，目标成像位置和尺寸与原模型（图 5a）基本吻合.

5 结论

本文提出迭代的解析方法快速完成多个导体圆柱目标的散射计算.首先考虑每个导体柱自身的一阶散射，在各圆柱的局部中心坐标系内用柱面波函数展开，然后经由坐标转换得到多目标总的一阶散射.分析表明：每个圆柱中心偏移产生的初始相位体现了各圆柱干涉作用对散射场的贡献.将各圆柱自身的一阶散射作为其他圆柱的激励场，利用柱面波函数的加法定理和边界条件推导出任意两个圆柱间的二次耦合散射表达式，经由坐标变换得到多目标总的二阶散射场.类似地，将二次散射作为激励场，经适当运算和变换得到三次甚至更高次的散射.由于每个圆柱目标的散射场最终由波函数展开解析计算，相对于MoM等数值方法而言，这种解析方法十分明显地提高了计算速度，尤其是对于后向散射计算.此外，高次散射是由低一次的散射作为激励场计算得到的，这种由低次到高次的一步步迭代计算，物理上十分清晰地计算了各圆柱间的各次耦合散射.

进一步利用这些宽频带宽角度的散射场数据，经由平滑和插值处理后，用FFT快速实现多目标ISAR成像模拟，很好地重构了各圆柱体目标的位置和尺寸.