摘要
与传统近场逆合成孔径雷达(inverse synthetic aperture radar,ISAR)成像算法相比,基于传统Omega-K近场成像算法的ISAR成像不会产生几何形变。但是由于雷达平台所需载频高、成像积累时间有限,这使得采用传统Omega-K近场成像算法直接在空间域进行聚焦时会出现图像混叠,且需要进行大量补零操作,无法实现快速高效成像。针对上述问题提出了一种波数域近场ISAR成像及形变校正算法。通过分析斜视成像二维波数域支撑区的斜拉特性,结合坐标旋转和方位谱重采样的方法,实现了波数谱的正侧化,再将子孔径数据聚焦在方位波数域,实现了方位向聚焦,最后通过Chirp-Z变换进行形变校正。经仿真验证,所提算法在解决形变校正问题的同时还能大幅减少运算量,显著提升了成像效率。
Abstract
Inverse synthetic aperture radar (ISAR) imaging does not produce geometric deformation based on the traditional Omega-K near-field imaging algorithm compared with the traditional near-field ISAR imaging algorithm. However, due to the high carrier frequency required by the radar platform and the limited imaging accumulation time, the imaging aliasing occurs when the traditional Omega-K near-field imaging algorithm is used to focus directly in the spatial domain, and a large number of zero-padding operations is required, which makes it impossible to achieve fast and efficient imaging. Aiming at the above problems, a near-field ISAR imaging and deformation correction algorithm in wavenumber domain was proposed. By analyzing the cable-stayed characteristics of the two-dimensional wavenumber domain support area of squint imaging, combined with coordinate rotation and azimuth spectrum resampling, the positive lateralization of the wavenumber spectrum was realized. The sub-aperture data was focused in the azimuth wavenumber domain, and the azimuth focusing was realized. Finally, the deformation correction was performed by Chirp-Z transform. The simulation results show that the proposed algorithm can greatly reduce the computational complexity and improve the imaging efficiency while solving the deformation correction problem.
0 引言
近年来,近场逆合成孔径雷达(inverse synthetic aperture radar,ISAR)成像的应用领域越来越广泛,特别是在隐蔽目标检测和目标散射特性测量等方面,近场ISAR成像展示出了重要的应用价值。在隐蔽目标检测方面,利用电磁波对非金属的穿透探测能力,毫米波近场ISAR成像技术逐渐应用于机场、高铁等的安检领域,并能够快速检测出隐蔽的危险物品[1-4]。在目标散射特性测量方面,基于ISAR成像的近场目标散射特性测量系统具有可快速部署、制造成本低、可测量大尺寸目标等优势,为目标散射特性测量提供了新的技术途径[5-6]。
受限于器件的发展水平,国内在雷达近场成像方面的研究近年来才开始兴起。张彪等[7]采用基于格林函数分解的成像算法实现了太赫兹近场成像;杨啸宇[8]对毫米波近场阵列成像方法以及基于图形处理器(GPU)的并行加速方法开展了研究;符吉祥[9]提出了一种基于子孔径的空变误差补偿近场ISAR成像方法。随着近场成像技术研究的不断深入,对近场成像性能及图像处理的要求也在逐渐提高。一方面,小尺寸隐蔽目标安检、高精度目标散射特性测量等应用对近场ISAR成像的效率和分辨率提出了更高要求;另一方面,基于线性调频(LFM)信号和宽带合成信号的高分辨ISAR成像算法存在回波处理数据量大、计算效率低等问题[10]。这些使得高分辨近场ISAR的快速成像成为亟待解决的问题。
与传统近场ISAR成像算法相比,采用传统Omega-K空间域近场成像算法的ISAR图像不会产生几何形变。但由于雷达平台所需载频非常高、成像积累时间十分有限,使得直接采用传统Omega-K近场成像算法在空间域进行聚焦时会出现图像混叠,这导致在ISAR成像时需要进行大量补零操作,无法实现快速高效成像。针对上述问题,本文首先对传统近场ISAR成像算法的几何形变问题进行分析;然后对传统Omega-K空间域近场成像算法及其运算量进行分析;最后提出波数域近场ISAR成像及形变校正算法,并对其形变校正性能和运算量进行仿真验证。
1 传统近场ISAR成像分析
ISAR匀速转动平台斜距模型如图1所示。
图1ISAR匀速转动平台的斜距模型
建立二维直角坐标系oxy,目标旋转中心位于原点o,(x,y)为目标散射点坐标;雷达在y轴上,R0为雷达与目标旋转中心间的中心斜距,R(ta)为雷达与目标散射点间的瞬时斜距,其中ta为方位向时间;θ为目标散射点绕中心旋转的初始角;ρ为初始长度,即旋转半径;ωta为旋转增量,其中ω为旋转角速度。
根据图1,瞬时斜距R(ta)的表达式为
(1)
设目标散射点(x,y)的初始坐标x=ρcosθ,y=ρsinθ,则式(1)可以表示为
(2)
当成像场景中旋转增量较小,即ωta→0时,式(1)存在近似。将式(1)在ωta=0处展开,可得
(3)
根据式(3),散射点的聚焦位置坐标(x′,y′)的表达式为
(4)
(5)
可见散射点的聚焦位置与真实位置之间存在偏移,ISAR图像发生形变,且随着R0变小,形变会越来越严重。为验证R0变化对图像形变的影响,采用传统的包络对齐自聚焦ISAR成像算法,分别进行远场(中心斜距为10km)和近场(中心斜距为50m)ISAR成像仿真,其仿真条件如表1所示。
表1近远场成像仿真条件
基于传统成像算法的ISAR近远场成像结果对比如图2所示。能够明显看出,与图2(a)相比,图2(c)所示的近场成像结果不再是矩形,目标散射点聚焦位置发生偏移,图像发生了形变。
图2基于传统成像算法的ISAR近远场成像结果
2 传统Omega-K算法空间域近场成像
2.1 传统Omega-K算法
在高速运动平台上进行ISAR成像时,目标与雷达平台之间的运动成像模型可以等效为静止目标与雷达平台间的斜视成像模型。Omega-K算法的空间域成像模型如图3所示。建立三维直角坐标系OXYZ,雷达平台沿正X轴方向以速度v飞行,点B为场景目标中心点,θ0为雷达波束中心斜视角,θBW为雷达方位向波束宽度,H为雷达平台飞行高度,RB为雷达平台与场景目标中心点B间的最小斜距,Xn为目标散射点n与场景目标中心点B间的方位向距离。
图3Omega-K算法的空间域成像模型
设X=vta为雷达平台的方位向采样点坐标,则目标散射点瞬时斜距
(6)
假设雷达发射信号为线性调频信号,雷达的基带回波可以表示为
(7)
式中:Rr为雷达平台的距离向采样点坐标;wr(Rr),wa(X)分别为距离向和方位向的空间域窗函数;γ为调频斜率;c为光速;λ为波长。
对式(7)沿距离向进行快速傅里叶变换(FFT),得到频域基带回波
(8)
式中:Kr=Krc+ΔKr为距离向波数,其中Krc=4πfc/c为距离向中心波数,fc为载波中心频率,ΔKr∈[-2πγTp/c,2πγTp/c]为距离向波数范围,Tp为发射脉冲宽度;wr(Kr)为距离向波数域窗函数。
设Kx,c0为方位向波数谱的偏移量,则构造的中心校正函数
(9)
将式(8)乘以式(9)进行中心校正,然后沿方位向作FFT,得到中心校正后的基带回波二维波数谱
(10)
其中
(11)
式中:Kx为方位向波数。构造距离向匹配函数,其表达式为
(12)
将式(10)乘以式(12)进行距离向匹配滤波,得到距离校正后的基带回波二维波数谱
(13)
与式(10)相比,式(13)仅保留了相位项。设Rs为场景中心斜距,则构造的相位补偿函数
(14)
将式(10)乘以式(14),得到相位补偿后的基带回波二维波数谱
(15)
这里采用Stolt插值进行解耦合,引入新的距离向波数Ky,则波数域的映射关系表达式为
(16)
Stolt插值后的二维波数谱表达式为
(17)
对式(17)进行距离向快速傅里叶逆变换(IFFT),得到的距离向聚焦的频域基带回波
(18)
式中:Cr为距离向带宽常数。根据不同的距离单元构造方位向偏移校正函数,其表达式为
(19)
将式(18)与式(19)相乘,最后进行方位向IFFT,即可实现方位向聚焦。方位向聚焦后的时域基带回波
(20)
式中:Ca为方位向带宽常数。由式(20)可以看出,目标散射点最终在空间域聚焦,且方位向目标散射点的聚焦位置与实际位置相同,图像未发生形变。
虽然采用Omega-K算法进行空间域成像能够实现良好聚焦且不会产生图像形变,但是由于雷达载频高、成像积累时间有限、斜距短等原因,导致与雷达的运动距离相比目标场景明显较大,这使得直接在空间域进行聚焦成像时会出现图像混叠问题。
采用Omega-K算法进行斜视子孔径成像时,图像的混叠问题是方位向空间域支撑区较短导致的。图4给出了子孔径成像的方位向空间域支撑区示意图,其中Xsub为子孔径数据录取长度,L为实际合成孔径长度。
图4子孔径成像的方位向空间域支撑区示意图
可以看出,在雷达录取数据时,其子孔径数据录取长度为Xsub,而雷达照射范围内的数据录取长度为Xsub+L。在实际应用中,Xsub一般远小于L。由于传统的Omega-K类算法均为二维空间域成像算法,数据录取长度即为最终成像场景的空间域支撑区长度。而在子孔径成像时,由于数据录取长度较短,即空间域支撑区较短,而实际照射场景区域远大于空间域支撑区,最终在进行方位向IFFT时,会导致超出空间域支撑区的场景目标图像发生混叠。如图4中的粗实线所示,仅处于Xsub范围内的目标散射点在方位向聚焦在正确位置上,而处于混叠区的目标散射点会在Xsub范围内发生混叠。
为了解决图像混叠问题,必须要进行补零操作。但补零操作明显增加了成像算法的运算量,这不利于快速成像处理。
2.2 运算量分析
下面对补零操作后成像算法的运算量进行分析。
首先对方位向数据进行补零,扩展方位向空间域支撑区,扩展后的支撑区应覆盖方位向的场景宽度。数据补零后的方位向点数
(21)
式中:Naz为初始方位向点数;N0为补零数。
整个算法包括两次距离向FFT操作、两次方位向FFT操作、四次复乘操作以及一次插值操作。设初始数据大小为Nrg×Naz,其中Nrg为初始距离向点数。经过方位向补零后,数据的方位向点数为Na。每次距离向或补零后方位向FFT操作的浮点运算量为5NrgNalog2Nrg或5NrgNalog2Na,而每次复乘操作则需要6NrgNa次浮点运算。相比较而言,插值操作的浮点运算量最大,每次为2(2Mker-1)NrgNa,其中Mker为插值核长度。则传统Omega-K算法总运算量
(22)
3 波数域近场ISAR成像及形变校正算法
为了实现精确快速且无形变的聚焦处理,本文提出了波数域近场ISAR成像及形变校正算法。首先,通过分析斜视成像二维波数域支撑区的斜拉特性,结合坐标旋转和方位谱重采样方法,实现了波数谱的正侧化,以实现支撑区的最大化利用并提高成像质量,从而消除斜视带来的影响。在方位向聚焦方面,该算法采用了类似频谱分析(spectral analysis)的处理方法,将子孔径数据聚焦在方位向波数域,最后通过Chirp-Z变换(CZT)进行形变校正。这种方法避免了对子孔径方位向空间域支撑区进行大量补零操作,从而有效减少了成像算法的运算量。
3.1 二维波数谱正侧化处理
由于波数域支撑区倾斜较大,导致可利用的波数域支撑区较小,这会使成像分辨率恶化,影响最终的成像效果。传统的处理方法[11]是沿视线方向进行插值。该方法需要先将二维波数域进行坐标旋转,再进行两次插值处理,将倾斜的支撑区旋转为正侧视的支撑区。上述操作会极大增加算法的运算量,不适合快速成像。采用传统斜视处理的距离走动时域校正方法可以实现波数域的坐标旋转[12-14],这里考虑构造一个波数域旋转校正函数来实现。坐标旋转校正前后波数域支撑区示意图如图5所示,图中θ2为初始倾斜角,K′x,K′y为波数域坐标旋转后的新波数。可见,旋转校正后波数域支撑区实现了最大程度利用。
图5坐标旋转校正前后波数域支撑区示意图
对式(8)进行距离向匹配滤波,得到频域基带回波,其表达式为
(23)
构造波数域旋转校正函数,其表达式为
(24)
将式(23)与式(24)相乘,得到坐标旋转变换后的频域基带回波,其表达式为
(25)
用Φ(Kr,X)表示式(25)中的指数项相位,对Φ(Kr,X)进行泰勒展开,可得
(26)
式中:O(·)为高阶无穷小函数;k为高阶无穷小函数的阶数。
针对空变问题[15]进一步分析回波的波数谱。对式(25)进行方位向傅里叶变换,得到回波二维波数谱
(27)
受距离走动线性校正的影响,雷达与目标旋转中心的斜距由R0变为R′0=R0+Xnsin θ0,将其代入式(27)并进行整理,可得
(28)
由式(28)的指数项可以看出,Xn与方位向波数Kx的二次项相耦合,会产生调频斜率的空变。本文采用波数域方位谱重采样方法,通过方位向插值引入新的波数K′x,使得
(29)
即只存在波数K′x的一次项与Xn的耦合项(该耦合项反映了目标的方位向位置),从而消除了波数Kx的二次项以及高次项与Xn的耦合项,这样就能够进行方位向的聚焦处理。将式(29)代入式(27),整理得到的正侧化处理后的二维波数谱
(30)
式(30)表明,对多普勒中心和距离走动进行校正,并采用方位向重采样插值处理可以消除方位空变。此时,原本斜视点目标的二维波数谱已经完全等效于正侧视成像的情况。
3.2 波数域近场成像
结合前面提到的正侧化思想提出了一种波数域近场ISAR成像算法。为了解决空间域成像支撑区不足的问题,该算法选择在波数域实现方位向的聚焦成像。构造一致聚焦函数
(31)
将式(30)与式(31)相乘,进行统一相位补偿,将距离向目标散射点位置调整到以场景中心为参考中心,得到
(32)
接下来进行扩展Stolt插值。不同于传统Stolt插值,扩展Stolt插值保留了方位向的波数项,便于方位向的波数域聚焦。
令
(33)
经过扩展Stolt插值的回波二维波数谱可以表示为
(34)
对式(34)作距离向IFFT,即可实现距离向的聚焦。距离脉压后的频域回波信号
(35)
构造统一的“去斜”校正函数,其表达式为
(36)
将式(35)与式(36)相乘,得到“去斜”处理后的频域回波信号
(37)
对式(37)作方位向IFFT,得到时域回波信号
(38)
构造方位空间域“扳直”校正函数,其表达式为
(39)
则经方位向“扳直”校正处理后的时域回波信号
(40)
最后对式(40)作方位向FFT,得到方位向波数域聚焦的回波信号
(41)
实现方位向波数域聚焦,最大无模糊成像的方位幅宽Xwid=fpvcos θ0/γa,其中fp为发射脉冲重复频率,γa为方位向调频斜率[16]。在通常情况下,无需进行方位向补零或仅需少量补零即可最终实现无混叠的成像。因此,相比于传统的Omega-K算法,本文所提算法降低了运算量,提升了成像效率。
3.3 形变校正
由前面推导可知,在方位向波数域将目标散射点n聚焦在4πcos2θ0Xn/(λR′0)位置处时,由于该位置与R′0相关,因此最终成像结果将会出现一定的形变。该形变可以在最终方位向聚焦时,通过使用CZT替代FFT来进行方位向校正。基于CTZ的形变校正流程如图6所示,图中x(m,n)表示成像完成后任一距离单元的数据,x(m,n′)表示形变校正后的数据。
图6基于CTZ的形变校正流程
以场景中心为参考中心,定义尺度因子ε,则CZT中使用的两个复数A0,W的表达式为
(42)
定义线性相位函数
(43)
线性相位函数用以保证CZT后场景中心不会发生偏移,其中定义尺度因子ε则是为了保证变换前后目标散射点回波的幅度基本不变。
至此就完成了对目标散射点的距离向空间域和方位向波数域的聚焦处理,同时也完成了形变校正。波数域近场ISAR成像及形变校正算法处理流程如图7所示。
图7波数域近场ISAR成像及形变校正算法处理流程
3.4 运算量分析
波数域近场ISAR成像及形变校正算法包括了两次距离向FFT操作、三次方位向FFT操作、五次复乘操作以及两次插值操作。设ISAR图像有Nrg×Naz个像素单元,每次距离向或方位向FFT操作的浮点运算量为5NrgNazlog2Nrg或5NrgNazlog2Naz,每次复乘操作的浮点运算量为6NrgNaz,每次插值操作的浮点运算量为2(2Mker-1)NrgNaz,则波数域近场ISAR成像及形变校正算法总运算量
(44)
4 仿真实验与结果分析
为了验证所提算法的有效性,将分别从形变校正和运算量两方面进行仿真实验与结果的对比分析。
4.1 点目标模型成像形变校正
采用如图2(a)所示的目标散射点分布,仿真验证所提算法的有效性。设场景宽度为30m,雷达系统参数如表2所示。
表2雷达系统参数
形变校正前后的目标散射点成像结果如图8所示,图中图框上方为远端,下方为近端,红色直线为垂线。
图8形变校正前后目标散射点成像结果
由图8可见,近场条件下各个散射点的聚焦位置在距离向和方位向都出现了明显偏移,目标散射点图像整体形状发生改变,而形变校正后目标散射点图像整体形状恢复为矩形。
单点插值后的目标散射点成像结果如图9所示。从单点插值结果可以看出各散射点聚焦良好。
图9单点插值后的目标散射点成像结果
4.2 F16飞机模型成像形变校正
设场景宽度为20m,雷达系统参数如表2所示。F16飞机模型散射点示意图如图10所示。形变校正前后F16飞机模型成像结果如图11所示。
图10F16飞机模型散射点示意图
图11形变校正前后F16飞机模型成像结果
从图 10 和图 11 中能够明显看到:成像形变校正前尾翼部分不对称;通过本文所提算法进行成像形变校正后,尾翼形变回复良好,尾翼对称且散射点聚焦良好。
4.3 运算量对比分析
2.2 节和3.4节已经分别对两种算法作了理论运算量分析,现将通过具体参数进行运算量对比。运算量对比所需参数如表3所示。
表3运算量对比所需参数
传统Omega-K成像算法首先需要补零。经计算,补零数N0=6 000,则方位向点数Na=6 512。将Na代入式(22),可得其总运算量Ncom1约为1.0×109。波数域近场ISAR成像及形变校正算法无需补零,其总运算量Ncom2约为9.9×107。
对比两种算法的实际运算量可以看出,本文所提的波数域近场ISAR成像及形变校正算法的运算量的数量级明显小于传统Omega-K成像算法的。算法运算量的大幅度减小能够有效提升成像效率。
5 结束语
相比于传统的基于Omega-K的空间域近场ISAR成像算法,本文所提波数域近场ISAR成像及形变校正算法无需进行方位向补零或仅需少量补零操作,即可实现无混叠、无形变的近场ISAR成像。该算法运算量小,有效提升了成像效率。