CN101900812A - 一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法 - Google Patents

Abstract

本发明公开了一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，将用于聚束合成孔径雷达的传统极坐标格式成像方法发展为适用于圆迹合成孔径雷达的三维成像方法，传统极坐标格式成像方法只校正了一次越距单元离徙动项，其平面波假设造成距离弯曲，因此只能对远场小场景进行成像，本发明在原有极坐标格式成像方法的基础上，还通过方位频域的二次滤波以及图像域的空变滤波这两次校正去除了二次越距离单元徙动项，较传统极坐标格式成像方法，可以对更大的场景进行精确的三维成像。

Description

一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法

技术领域

本发明涉及三维成像技术领域，特别是一种将大场景极坐标格式算法应用于圆迹合成孔径雷达的成像方法。

背景技术

近年来，圆迹合成孔径雷达(circular synthetic aperture radar，CSAR)因其具有全方位观测和高分辨三维成像能力开始受到人们的广泛关注。2004年，法国宇航实验室(ONERA)在瑞典开展了第一次的L波段和P波段的机载CSAR实验，2006年该实验室又在法国和德国进行了X波段的机载CSAR实验，得到了分辨率为0.15m的全方位观测图像。

CSAR通过在方位向形成的圆形合成孔径，对目标区域进行360°全方位观测，其波束始终指向同一目标区域。它独具的大合成孔径角使其能够达到工作波长量级的平面分辨率，这是传统直线合成孔径雷达(以下简称为线性SAR)：条带模式SAR和聚束模式SAR所达不到的；发射信号为宽带信号的条件下，二维合成孔径又决定了CSAR系统具有三维分辨能力。单次圆周飞行即可得到高分辨三维图像的特点使CSAR系统在军事应用领域中非常具有实际意义。

现有的CSAR成像算法主要分两类：参数化方法和非参数化方法。参数化方法是指应用现代普估计的方法提取散射点的幅度信息和位置信息，基于RELEX的三维目标特征提取方法是目前比较常用的。非参数化方法是通过聚焦的方法重建目标函数，主要有：后向投影(Back projection，BP)算法，共焦成像算法，基于格林函数傅里叶分析的三维成像算法等。BP算法首先计算图像中每个象素点对应的积分路径，然后对该路径插值并积分，完成该象素点的聚焦。共焦成像算法与BP算法原理类似，也要先计算每个象素点的路径，然后通过空变滤波得到聚焦图像。这两种算法都需要对图像中每个象素点逐一操作，计算量比较大。基于格林函数傅里叶分析的三维成像算法，采用逐高度平面成像的方法，通过格林函数的分解，对回波进行统一操作即可重建目标函数，较前两种CSAR成像算法效率更高，但是该算法中的求逆系统核以及反卷积等逆运算增大了算法的复杂性，而且算法对雷达平台的运动稳定性要求较高。

在现有圆迹合成孔径的暗室或机载实验中，由于雷达平台运动控制的难度，多采用BP算法以得到较精确的成像结果，但是在高分辨或大场景的实际应用条件下，BP算法的效率将成为主要问题。随着CSAR研究工作的不断开展和深入，迫切需要一种适用于CSAR的适应性强、高效、简洁、准确的成像算法。

极坐标格式算法(Polar Format Algorithm，PFA)是聚束模式信号处理的经典算法，最早是作为一种有效的旋转目标成像方法提出的，但很快就被成功用于聚束模式SAR成像中，PFA算法以其简单、高效、计算量小和适用性好等优点，一直被广泛应用于聚束模式成像领域。

极坐标格式算法利用解线频调(Dechirp)技术以场景中心的回波信号作为参考信号，对回波信号解调和相位补偿。解调后场景中心的越距离单元徙动(range cell migration，RCM)项被完全消除，场景中其余点还存在剩余RCM。然后通过极坐标系到直角坐标系的插值操作，去除残余RCM中的一次项。

PFA算法存的主要问题是它的平面波假设，而实际发射微波为球面波，平面波假设将造成距离弯曲，他对图像质量的影响主要有两个方面：一次相位误差造成图像的几何失真，二次两位误差造成散焦。这限制了PFA算法在大场景成像中的应用。2004年由Walter G.Carrara等人提出的对信号相位历史的沿航迹方向校正和格式化算法(along-track alignment andformatting system，ATAFS)就很好的解决了聚束SAR中的这一问题。将ATAFS算法与PFA算法相结合的大场景极坐标格式算法(Widefield PolarFormat Algorithm，WPFA)可以获得较高的图像质量，且对场景大小没有限制。

CSAR波束始终指向同一场景区域，可以被认为是一种特殊的聚束模式SAR，因此PFA算法在CSAR成像领域应该具有一定的可移植性，但目前并未见相关报道。

一般的CSAR的成像几何关系如图1所示，以场景中心作为坐标系原点建立直角坐标系OXYZ，雷达平台在高度H，围绕原点作平行于X-Y平面，半径为R的匀速圆周运动，波束始终指向场景区，场景被照射区横截面半径为X₀，高度为2h₀，场景中心与雷达的仰角设为α，场景中心与雷达的距离设为R₀。θ为方位向角度，每隔角度间隔Δθ，雷达发射一次脉冲。则雷达在极坐标系下的坐标为(R₀，α，θ)，其中

场景中的任意点目标设为P，在极坐标系下的坐标设为

其中，ρ为点目标P位置向量的长度，η为点目标P位置向量与X-Y平面的夹角，

为点目标P位置向量在X-Y平面的投影与X轴的夹角，设点目标P在直角坐标系下的坐标设为(x_p，y_p，z_p)，其散射系数设为σ_p。

发明内容

本发明的目的是公开一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，针对高分辨和大场景的应用条件，将用于聚束SAR的二维极坐标格式成像方法发展为适用于CSAR的三维成像算法，并通过算法改进，解决传统PFA算法中由平面波假设带来的距离弯曲问题。是一种适用于高分辨大场景的高效、简洁、准确的圆迹合成孔径雷达三维成像方法。

为达到上述目的，本发明的技术解决方案是：

一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其包括步骤如下：

步骤S1：将圆迹合成孔径雷达(Circular SAR，CSAR)的雷达原始回波信号送入距离向处理单元，沿原始回波信号的距离向变换到距离压缩频域，并通过参考函数共轭相乘处理，使场景中心处的相位历程为0，生成第1信号S₁(θ，k)，其中，θ为方位向雷达平台旋转角，k为斜距向波数；

步骤S2：将第1信号送入高度向处理单元，乘以高度滤波函数，生成第2信号S₂(θ，k)，在正交直角坐标系OXYZ下，以场景中心为原点，高度滤波函数为关于成像平面高度z₀的函数H₂(z₀；k)，其中，z₀＝-h₀+i×Δz，i＝0，1，2…，Δz为成像高度间隔，h₀为场景区高度的一半，由i＝0开始，至z₀＞h₀结束；更新参数cosα和R₀，α为Z轴上高度z₀点相对雷达的仰角，R₀为Z轴上高度z₀点与雷达的距离；

步骤S3：将第2信号送入方位向处理单元，通过方位向傅里叶变换变换至方位频域，然后作方位向二次滤波，最后经傅里叶逆变换变换至方位空间域，生成第3信号S₃(θ，k)，第3信号已去除了部分二次越距离单元徙动(range cell migration-RCM)项；

步骤S4：将第3信号送入坐标变换单元，由极坐标格式数据经二维插值得到直角坐标系下的数据，生成第4信号S₄(k_x，k_y)，其中，k_x为X方向波数，k_y为Y方向波数，第4信号已去除了一次RCM项；

步骤S5：将第4信号送入二维逆傅里叶变换单元，生成第5信号S₅(x，y)，第5信号为成像高度z＝z₀，其中，x为X方向空间域，y为Y方向空间域，z为Z方向空间域，S₅(x，y)残留部分二次及高次RCM的二维图像；

步骤S6：将第5信号送入空变滤波单元，将第5信号表示的二维图像按不同半径分区，各分区分别作二维傅里叶变换，然后与各分区对应的滤波器H₆(l；k_x，k_y)相乘，其中l＝1，2，...为分区序号，最后将各分区分别作二维傅里叶逆变换，变换至二维图像域，各分区拼接后生成第6信号S₆(z₀；x，y)。

步骤S7：判断单元，若z₀+Δz＜h₀，则令i＝i+1，并转到步骤S2，否则直接进入步骤S8；

步骤S8：将第6信号送入三维图像输出单元，对各高度层的二维图像按高度顺序排列，获得最终的CSAR三维图像：第7信号S₇(x，y，z)。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其所述步骤S1中距离向处理单元的具体步骤为：

步骤S11：如果回波信号为步进频连续波信号，则将回波信号与参考函数H₁(θ，k)共轭相乘，生成第1信号S₁(θ，k)，否则转入步骤S12。

步骤S12，如果发射信号为chirp信号，且回波已经二维解线频调(dechirp)处理，则进行距离去扭(range deskew)操作，生成第1信号S₁(θ，k)，否则转入步骤S13；

步骤S13：如果发射信号为线性调频(chirp)或其他宽带信号，且回波未作其他处理，则先进行距离压缩，然后在距离频域-方位空间域，与参考函数H₁(θ，k)作共轭相乘，生成第1信号S₁(θ，k)，否则退出程序。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其所述步骤S3中方位向处理单元的具体步骤为：

步骤S31：将第2信号沿方位向作傅里叶变换，得到二维频域信号：第31信号S₃₁(k_θ，k)，其中k_θ为方位波数；

步骤S32：将第31信号与方位二次滤波函数H₃(k_θ，k)相乘，得到第32信号S₃₂(k_θ，k)，这一操作去处了部分二次RCM项；

步骤S33：将第32信号作方位向傅立叶逆变换，变换到方位空间域，生成第3信号S₃(θ，k)。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其所述步骤S4中坐标变换单元的具体步骤为：

步骤S41：通过直角坐标系中均匀分布的X方向波数域值k_x、Y方向波数域值k_y计算对应的极坐标系下的斜距向波数k，及角度θ；

步骤S42：通过二维插值将极坐标系下的数据转换到直角坐标系，生成第4信号S₄(k_x，k_y)。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其所述步骤S6中空变滤波单元的具体步骤为：

步骤S61：将第5信号表示的二维图像按不同半径分成互不重叠的若干分区，分区准则为分区内外圈上的目标，残余二次RCM引起的相位误差差值不超过π，表达式为：

$k_{\max }{\sin }^2\mathrm{\&alpha;}({r_{l+1}}^2-r_l^2)/\left(2R_0\right)<\mathrm{\&pi;}$

其中，k_max为最大斜距向波数，r_l+1，r_l分别为第1个分区的外圈半径和内圈半经，π为圆周率。

步骤S62：各分区分别作二维傅里叶变换，得到各分区二维频谱；

步骤S63：各分区频谱与各自的相位补偿滤波器H₆(l；k_x，k_y)相乘，生成第63信号S₆₃(l；k_x，k_y)；

步骤S64：将第63信号各自作傅立叶逆变换，回到图像域，并拼接为第6信号S₆(z₀；x，y)；

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，在CSAR的成像几何中，雷达发射脉冲的角度间隔Δθ应小于或等于

其中R为雷达平台轨道半径，H为雷达平台高度，X₀为场景照射区半径。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，所述步骤S2中成像高度间隔Δz应小于或等于

其中k_min为最小斜距向波数。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其所述步骤S32中方位二次滤波函数H₃(k_θ，k)的表达式为：

H₃(k_θ，k)＝exp[jk_θ ²/(2kR₀)]。

所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其所述步骤S63中相位补偿滤波器H₆(l；k_x，k_y)的表达式为： $H_6(l;k_x,k_y)=\exp [j\sqrt{{k_x}^2+{k_y}^2}\sin \mathrm{\&alpha;}\tan \mathrm{\&alpha;}(r_l^2+{r_{l+1}}^2)/4/R_0].$

本发明的有益效果：本发明提供一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，将传统极坐标格式算法发展为适用于圆迹合成孔径雷达的三维成像方法，本发明可操作性强，效率高，且对雷达平台不稳定运动有很强的适应能力，因此非常具有现实意义。

本发明通过预处理单元不仅能够处理发射为chirp信号的回波信号，以及发射信号为chirp信号且回波经二维dechirp后的信号，还能处理发射信号为步进频连续波信号的回波信号，因而对于发射信号形式的要求降低，便于实际应用。

本发明利用方位向频域二次滤波以及分区域相位校正的方法分两次去除了二次RCM带来的相位误差，解决了传统极坐标格式算法对成像区域尺寸的限制问题，可以实现大场景精确成像。

本发明具有可操作性强，高效，精确，适应性强等特点，有很大的实际意义。

附图说明

图1是圆迹合成孔径雷达三维成像几何示意图；

图2是本发明的计算机或专用设备处理流程图；

图3是本发明的距离向处理单元处理流程图；

图4是圆迹合成孔径雷达中雷达与目标几何关系示意图；

图5是本发明的方位向单元处理流程图；

图6是本发明的坐标变换单元处理流程图；

图7是本发明的空变滤波单元处理流程图；

图8是本发明的空变滤波处理单元中的分区示意图；

图9是传统极坐标格式算法的成像结果；

图10是本发明实施例1的17个点目标的成像结果；

图11是本发明实施例的图9中点目标A的放大图；

图12是本发明实施例的图9中点目标B的放大图；

图13是本发明实施例2的23个点目标的成像结果侧视图；

图14是本发明实施例2的23个点目标的成像结果正视图；

图15是本发明实施例2的23个点目标的成像结果俯视图。

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

如图2本发明提出的圆迹合成孔径雷达的大场景极坐标格式三维成像方法的流程图所示，输入为CSAR的原始回波，经过三维成像处理后，可获得CSAR三维图像，具体实施步骤如下：

步骤S1：将CSAR的雷达原始回波信号送入距离向处理单元，沿原始回波信号的距离向变换到距离压缩频域，并通过参考函数共轭相乘处理，使场景中心处的相位历程为0，生成第1信号S₁(θ，k)，其步骤如图3所示，具体实施过程为：

步骤S11：如果回波信号为步进频连续波信号，则将回波信号与参考函数H₁(θ，k)共轭相乘，生成第1信号S₁(θ，k)，否则转入步骤S12。

如果回波信号为步进频连续波信号，则来自点目标P的回波表达式为：

s_r(θ，k)＝σ_pexp{-jkR_p}    (1)

其中

为点目标P与雷达之间的距离，其中，(x_p，y_p，z_p)为点目标P的直角坐标，为描述简便，省略公式前的常系数，该简化对本方法描述不造成影响。参考函数H₁(θ，k)的表达式为：

H₁(θ，k)＝exp(-jkR₀)    (2)

将回波信号与H₁(θ，k)共轭相乘，生成第1信号S₁(θ，k)，表达式为：

S₁(θ，k)＝exp[-jk(R_p-R₀)]    (3)

如果发射信号不是步进频连续波信号，转入步骤S12；

步骤S12，如果发射信号为chirp信号，且回波已经二维dechirp处理，则进行距离去扭操作，生成第1信号S₁(θ，k)，否则转入步骤S13；

如果发射信号为chirp信号，且回波已经二维dechirp处理，则回波信号表达式为：

$s_r(\mathrm{\&theta;},t)=\mathrm{rect}\left(\frac{t-{2R}_p/C}{T_p}\right)e^{-j\frac{4\mathrm{\&pi;}}{C}[\mathrm{\&gamma;}(t-\frac{{2R}_0}{C})+f_c]R_{\mathrm{\&Delta;}}+j\frac{4\mathrm{\&pi;\&gamma;}}{C^2}{R}_{\mathrm{\&Delta;}}^2}---\left(4\right)$

其中，t为快时间，γ为脉冲信号调频率，T_p为发射脉冲宽度，rect[t/T]表示|t|≤T，R_Δ＝R_p-R₀。经距离去扭，并转换到距离频域后，生成第1信号S₁(θ，k)，表达式如式(3)。

如果发射信号不是chirp信号，或回波未经二维dechirp处理，则转入步骤S13；

步骤S13：如果发射信号为chirp或其他宽带信号，且回波未作其他处理，则先进行距离压缩，然后在距离频域-方位空间域，与参考函数H₁(θ，k)作共轭相乘，生成第1信号S₁(θ，k)，否则退出程序。

如果发射信号为chirp或其他宽带信号，如果发射信号为chirp或其他宽带信号，且回波未作其他处理，其回波表达式为：

$s_r(\mathrm{\&theta;},t)=p(t-\frac{{2R}_p}{C})---\left(5\right)$

其中，p(t)为发射信号，对回波作距离向傅里叶变换并匹配滤波后得：

S_r(θ，k)＝|P(k)|²exp(-jkR_p)    (6)

将S_r(θ，k)与参考函数H₁(θ，k)作共轭相乘，生成第1信号S₁(θ，k)：

S₁(θ，k)＝|P(k)|²exp[-jk(R_p-R₀)]    (7)

为描述简便，省略公式前的|P(k)|²，该简化对本方法描述不造成影响，则式(7)形同式(3)，如果发射信号为其他信号，退出程序。

步骤S2：将第1信号送入高度向处理单元。乘以高度滤波函数H₂(z₀；k)：

$H_2(z_0,k)=\exp \{-\mathrm{jk}[\sqrt{R^2+H^2}-\sqrt{R^2+{(H-z_0)}^2}\left]\right\}---\left(8\right)$

其中z₀＝-h₀+i×Δz，由i＝0开始，至z₀＞h₀结束。更新参数，α为Z轴上z＝z₀点相对雷达的仰角，则

R₀为Z轴上z＝z₀点与雷达的距离，则

通过高度滤波生成第2信号S₂(θ，k)：

$S_2(\mathrm{\&theta;},k)=\exp \{-\mathrm{jk}[R_p-\sqrt{R^2+{(H-z_0)}^2}\left]\right\}---\left(9\right)$

如图4所示，β为雷达位置向量与点目标P位置向量之间的夹角，则由余弦定理可知P与雷达之间的距离R_p为：

$R_p=\sqrt{R_0^2+{\mathrm{\&rho;}}^2-{2R}_0\mathrm{\&rho;}\cos \mathrm{\&beta;}}=\sqrt{{(R_0-\mathrm{\&rho;}\cos \mathrm{\&beta;})}^2+{\left(\mathrm{\&rho;}\sin \mathrm{\&beta;}\right)}^2}---\left(10\right)$

当满足R₀-ρcosβ＞＞ρsinβ，即在远场条件下时，对ρsinβ作泰勒展开至二次项，可上式可近似为：

R_p≈(R₀-ρcosβ)+(ρsinβ)²/2(R₀-ρcosβ)    (11)

其中，传统PFA算法只校正了式(11)中的一次项，却忽略了二次项，因此限制了其在大场景成像中的应用，本发明的目的是通过相位补偿的方法校正距离弯曲的二次项，以适应大场景精确成像。将式(11)代入式(9)得第2信号S₂(θ，k)的近似表达式为：

S₂(θ，k)≈exp[jkρcosβ-jk(ρsinβ)²/2(R₀-ρcosβ)]    (12)

只考虑z＝z₀高度层的点目标的RCM，即η＝0，为了更清晰地描述RCM矫正的过程，以下步骤都只给出了η＝0时的表达式。在远场条件下，R₀-ρcosβ≈R₀，则式(12)右边相位的表达式可近似为：

步骤S3：将第2信号送入方位向处理单元，通过方位向傅里叶变换变换至方位频域，然后作方位向二次滤波，最后经傅里叶逆变换变换至方位空间域，生成第3信号S₃(θ，k)，第3信号已去除式(13)的第二项，其步骤如图5所示，具体实施过程为：

步骤S31：将第2信号沿方位向作傅里叶变换，得到二维频域信号，第31信号S₃₁(k_θ，k)：

S₃₁(k_θ，k)＝FT_θ{S₂(θ，k)}    (14)

根据驻定相位原理，方位向波数kθ的表达式可通过式(13)中的相位γ对θ求导数求得，在远场条件下，式(13)等号右侧的第一项远大于后两项，因此k_θ的表达式近似为

式(13)等号右侧的第三项相对于前两项是缓变的，则式(14)可近似为：

步骤S32：将第31信号与方位二次滤波函数H₃(k_θ，k)相乘，H₃(k_θ，k)的表达式为：

H₃(k_θ，k)＝exp[jk_θ ²/(2kR₀)]    (16)

得到第32信号S₃₂(k_θ，k)：

这一操作去处了部分二次RCM项；

步骤S33：将第32信号作方位向傅立叶逆变换，变换到方位空间域，生成第3信号S₃(θ，k)：

其中，

分别为点目标P在X方向和Y方向的直角坐标，k_x＝-kcosαcosθ，k_y＝-kcosαsinθ分别为X方向和Y方向的波数，(θ，kcosα)为极坐标，(k_x，k_y)为相应的直角坐标，(θ，k)到(k_x，k_y)的变换需通过步骤S4实现；

步骤S4：将第3信号送入坐标变换单元，由极坐标格式数据经二维插值得到直角坐标系下的数据，生成第4信号S₄(k_x，k_y)，第4信号已去除了一次RCM项，其步骤如图6所示，具体实施过程为：

步骤S41：通过直角坐标系中均匀分布的X方向波数值k_x、Y方向波数值k_y计算对应的极坐标系下的斜距向波数k’，及角度θ’为：

$\left\{\begin{array}{c}{k}^{\&prime;}=\sqrt{{k_x}^2+{k_y}^2}/\cos \mathrm{\&alpha;}\\ {\mathrm{\&theta;}}^{\&prime;}=\arg (-k_x-{\mathrm{jk}}_y)\end{array}\right.---\left(19\right)$

步骤S42：通过二维插值将极坐标系下的数据转换到直角坐标系，生成第4信号S₄(k_x，k_y)，sinc插值函数的表达式为：

$S_4(k_x,k_y)=S_3({\mathrm{\&theta;}}^{\&prime;},k^{\&prime;})$

$=\underset{m=0}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}{S}_3(\mathrm{m\&Delta;\&theta;},k_{\min }+\mathrm{n\&Delta;k})\sin c[({\mathrm{\&theta;}}^{\&prime;}-\mathrm{m\&Delta;\&theta;})/\mathrm{\&Delta;\&theta;}]\sin c[(k^{\&prime;}-k_{\min }-\mathrm{n\&Delta;k})/\mathrm{\&Delta;k}]---\left(20\right)$

其中，S₅(mΔθ，k_min+nΔk)为第3信号S₃(θ，k)的离散表达式，m＝0，1，2，...，M，n＝0，1，2，...，N，为离散序号，M为序列m的最后一个值，N为序列n的最后一个值。上式是理想二维插值的表达式，通常选择4点或8点插值函数即可满足精度要求，经坐标变换后，第4信号S₄(k_x，k_y)的近似解析表达式为：

$S_4(k_x,k_y)\&ap;\exp [-{\mathrm{jk}}_x{x}_p-{\mathrm{jk}}_y{y}_p-j\sqrt{k_x^2+k_y^2}\frac{{\left(\mathrm{\&rho;}\sin \mathrm{\&alpha;}\right)}^2}{2R_0\cos \mathrm{\&alpha;}}]$

$=\exp [-{\mathrm{jk}}_x{x}_p-{\mathrm{jk}}_y{y}_p-j\sqrt{k_x^2+k_y^2}\sin \mathrm{\&alpha;}\tan \mathrm{\&beta;}\frac{{\mathrm{\&rho;}}^2}{{2R}_0}]---\left(21\right)$

步骤S5：将第4信号送入二维逆傅里叶变换单元，生成第5信号S₅(x，y)：

S₅(x，y)＝IFT2_x，y[S₄(k_x，k_y)]    (22)

第5信号成像高度为z＝z₀，且残留部分二次及高次RCM的二维图像；

步骤S6：将第5信号送入空变滤波单元，将第5信号表示的二维图像按不同半径分区，各分区分别作二维傅里叶变换，然后与各分区对应的滤波器H₆(l；k_x，k_y)相乘，其中l∈N为分区序号，最后将各分区作二维傅里叶逆变换，变换至二维图像域，生成第6信号S₆(z₀；x，y)，其步骤如图7所示，具体实施过程为：

步骤S61：将第5信号表示的二维图像按不同半径分成互不重叠的若干分区，如图8的分区示意图所示，分区准则为分区内外圈上的目标，残余二次RCM，即式(13)等号右侧第二项，引起的相位误差差值不超过π，表达式为：

$k_{\max }{\sin }^2\mathrm{\&alpha;}({r_{l+1}}^2-r_l^2)/\left({2R}_0\right)<\mathrm{\&pi;}---\left(23\right)$

其中r_l+1，r_l分别为第1个分区的外圈半径和内圈半经，分区后生成第61信号S₆₁(l；x，y)：

S₆₁(l；x，y)＝win(l)S₅(x，y)

$\mathrm{win}\left(l\right)=\left\{\begin{array}{cc}1,& r_l\&le;\sqrt{x^2+y^2}<r_{l+1}\\ 0,& \mathrm{else}\end{array}\right.---\left(24\right)$

其中win(l)为分区1的窗函数。

步骤S62：对上述各分区S₆₁(l；x，y)分别作二维傅里叶变换，得到各分区二维频谱，即第62信号：S₆₂(l；k_x，k_y)＝FT_x，y{S₆₁(l；x，y)}；

步骤S63：各分区频谱与各自的相位补偿滤波器H₆(l；k_x，k_y)相乘，H₆(l；k_x，k_y)的表达式为：

$H_6(l;k_x,k_y)=\exp \left(j\sqrt{{k_x}^2+{k_y}^2}\sin \mathrm{\&alpha;}\tan \mathrm{\&alpha;}\frac{r_l^2+{r_{l+1}}^2}{{4R}_0}\right)---\left(25\right)$

生成第63信号S₆₃(l；k_x，k_y)：

$S_{63}(l,k_x,k_y)\&ap;{\mathrm{FT}2}_{x,y}\left[\mathrm{win}\left(l\right)\right]\&CircleTimes;\exp [-{\mathrm{jk}}_x{x}_p-{\mathrm{jk}}_y{y}_p]---\left(26\right)$

步骤S64：将第63信号各自作傅立叶逆变换，回到图像域，并拼接为第6信号S₆(z₀；x，y)：

$S_6(z_0;x,y)\&ap;\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\mathrm{FT}2\left\{S_{63}(l,k_x,k_y)\right\}---\left(27\right)$

步骤S7：判断单元，若z₀+Δz＜h₀，则令i＝i+1，并转到步骤S2，否则直接进入步骤S8；

步骤S8：将第6信号送入三维图像输出单元，对各高度层的二维图像按高度顺序排列，获得最终的CSAR三维图像：第7信号S₇(x，y，z)。

以位于三维空间中的点目标为例，雷达中心频率0.5GHz，信号带宽为0.25GHz，步进频率间隔为2.31MHz，雷达飞行轨迹半径800m，高度2000m，成像区域半径为100m，方位向采样间隔Δθ＝0.42°，成像区域共有17个点目标，且假设都位于同一平面：z₀＝0平面，其中1个点位于场景中心，8个点均匀分布在以场景中心为圆心，半径为50m的圆上，另8个点位于场景边缘，且分布在以场景中心为圆心，半径为100m的圆上。则通过上述大场景极坐标格式三维成像方法得到的17个点目标的成像结果如图10-12所示，图10为17个点目标在高度z₀＝0的二维图像，图像作了归一化处理，动态显示范围为【-40dB 0dB】，从图中看到，通过上述处理步骤后，能够精确地重建出17个点目标的空间位置和幅度信息，图11、图12分别是场景中心目标A与场景边缘目标B的放大图，动态显示范围为【-20dB 0dB】，可以看到场景边缘的点目标聚焦性能与场景中心的点目标聚性能一样好。图9是将传统极坐标格式算法应用于本例的成像结果，由于未作任何相位误差校正，目标的聚焦性能随着其与场景中心距离的增大而变差，直至完全散焦。

为了能更清晰地显示三维成像的效果，其他参数不变，缩小场景区范围，场景区的参数变化为：成像区半径为30m，成像区高度为20m，共设置了23个点目标，17个点目标位于z₀＝0平面，其中1个点目标位于原点处，8个点目标均匀分布在原点为圆心，半径为15m的圆上，8个点均匀分布在以圆点为圆心，半径为15m的圆上，另6个点均匀分布在Z轴上，间隔为3m。图13-15为通过本大场景极坐标格式三维成像方法获得的三维成像结果，图13-15分别为成像结果的侧视，正视和俯视的-3dB轮廓图。可以看到，23个点目标的幅度和位置信息都得到了很好的重建。

本发明应用用于聚束模式SAR的二维极坐标格式成像方法的原理，发明了适用于圆迹合成孔径雷达的极坐标格式三维成像方法。

本发明通过距离向处理单元不仅能够处理发射步进频连续波信号的CSAR，以及发射调频脉冲信号，且回波经二维解线频调处理的CSAR，还能处理发射信号为调频脉冲信号或其他宽带信号，且回波未经其他处理的CSAR，因而对于发射信号形式的要求降低，便于实际应用。

本发明通过方位向频率域的二次滤波，去除了部分二次RCM项带来的相位误差，该校正为统一校正，与目标位置无关。

本发明的空变滤波，通过对场景各分区进行分别操作，去除了残余二次RCM项带来的相位误差，该校正与目标与成像区中心的距离有关。

本发明通过两次相位误差校正，校正了二次RCM项带来的相位误差，而传统极坐标格式成像方法只去除了一次RCM的影响，因此，本发明较传统极坐标格式方法，能够对更大的场景进行精确聚焦成像。

本发明上述的方法，已经在计算机上应用MATLAB软件得到验证，图10-15是在MATLAB软件环境下使用本发明所述方法获得的结果。该方法在计算机或专用设备上主要通过8个程序模块来实现，如图2所示，所述距离向处理单元、高度向处理单元、方位向处理单元、坐标变换单元、二维逆傅里叶变换单元、空变滤波单元、判断单元和三维图像输出单元，分别完成步骤S1、步骤S2、步骤S3、步骤S4、步骤S5、步骤S6、步骤S7和步骤S8的功能。

距离向处理单元，输入为CSAR原始回波数据，输出为处于距离向频域、方位向时域的第1信号；

高度向处理单元，输入为第1信号，输出为处于高度z0＝0平面的距离向频域、方位向时域的第2信号；

方位向处理单元，输入为第2信号，输出为已去除部分二次RCM引起的相位误差的距离向频域、方位向时域的的第3信号；

坐标变换单元，输入为第3信号，输出为由极坐标数据格式变换为直角标数据格式的距离向频域、方位向时域的第4信号；

二维逆傅里叶变换单元，输入为第4信号，输出为残余二次及高次RCM引起的相位误差的图像域的第5信号；

空变滤波单元，输入为第5信号，输出为去除了残余二次RCM引起的相位误差的图像域的第6信号；

判断单元，输入为第6信号，如果已处理完所有高度层，则进入步骤S8，否则转入步骤S2；

三维图像输出单元，输入为第6信号，输出为CSAR的三维成像结果的第7信号。

本发明上述的方法，已得到验证，该圆迹合成孔径雷达的大场景极坐标格式成像方法已用于仿真和数据处理，本发明方法的有效性得到了验证。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的范围之内，因此，本发明的保护范围应以权利要求书的保护范围为准。

Claims (9)

1.一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，包括步骤如下：

步骤S1：将圆迹合成孔径雷达的雷达原始回波信号送入距离向处理单元，沿原始回波信号的距离向变换到距离压缩频域，并通过参考函数共轭相乘处理，使场景中心处的相位历程为0，生成第1信号S₁(θ，k)，其中，θ为方位向雷达平台旋转角，k为斜距向波数；

步骤S2：将第1信号送入高度向处理单元，乘以高度滤波函数，生成第2信号S₂(θ，k)，在正交直角坐标系OXYZ下，以场景中心为原点，高度滤波函数为关于成像平面高度z₀的函数H₂(z₀；k)，其中，z₀＝-h₀+i×Δz，i＝0，1，2...，Δz为成像高度间隔，h₀为场景区高度的一半，由i＝0开始，至z₀＞h₀结束；更新参数cosα和R₀，α为Z轴上高度z₀点相对雷达的仰角，R₀为Z轴上高度z₀点与雷达的距离；

步骤S3：将第2信号送入方位向处理单元，通过方位向傅里叶变换变换至方位频域，然后作方位向二次滤波，最后经傅里叶逆变换变换至方位空间域，生成第3信号S₃(θ，k)，第3信号已去除了部分二次越距离单元徙动(range cell migration，RCM)项；

步骤S4：将第3信号送入坐标变换单元，由极坐标格式数据经二维插值得到直角坐标系下的数据，生成第4信号S₄(k_x，k_y)，其中，k_x为X方向波数，k_y为Y方向波数，第4信号已去除了一次越距离单元徙动项；

步骤S5：将第4信号送入二维逆傅里叶变换单元，生成第5信号S₅(x，y)，第5信号为成像高度z＝z₀，其中，x为X方向空间域，y为Y方向空间域，z为Z方向空间域，S₅(x，y)残留部分二次及高次越距离单元徙动的二维图像；

步骤S6：将第5信号送入空变滤波单元，将第5信号表示的二维图像按不同半径分区，各分区分别作二维傅里叶变换，然后与各分区对应的滤波器H₆(l；k_x，k_y)相乘，其中l＝1，2，...为分区序号，最后将各分区分别作二维傅里叶逆变换，变换至二维图像域，各分区拼接后生成第6信号S₆(z₀；x，y)。

步骤S7：判断单元，若z₀+Δz＜h₀，则令i＝i+1，并转到步骤S2，否则直接进入步骤S8；

步骤S8：将第6信号送入三维图像输出单元，对各高度层的二维图像按高度顺序排列，获得最终的圆迹合成孔径雷达三维图像：第7信号S₇(x，y，z)。

2.根据权利要求1所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S1中距离向处理单元的具体步骤为：

步骤S11：如果回波信号为步进频连续波信号，则将回波信号与参考函数H₁(θ，k)共轭相乘，生成第1信号S₁(θ，k)，否则转入步骤S12。

步骤S12，如果发射信号为线性调频(chirp)信号，且回波已经二维解线频调(dechirp)处理，则进行距离去扭(range deskew)操作，生成第1信号S₁(θ，k)，否则转入步骤S13；

步骤S13：如果发射信号为chirp或其他宽带信号，且回波未作其他处理，则先进行距离压缩，然后在距离频域-方位空间域，与参考函数H₁(θ，k)作共轭相乘，生成第1信号S₁(θ，k)，否则退出程序。

3.根据权利要求1所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S3中方位向处理单元的具体步骤为：

步骤S31：将第2信号沿方位向作傅里叶变换，得到二维频域信号：第31信号S₃₁(k_θ，k)，其中k_θ为方位波数；

步骤S32：将第31信号与方位二次滤波函数H₃(k_θ，k)相乘，得到第32信号S₃₂(k_θ，k)，这一操作去处了部分二次越距离单元徙动项；

步骤S33：将第32信号作方位向傅立叶逆变换，变换到方位空间域，生成第3信号S₃(θ，k)。

4.根据权利要求1所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S4中坐标变换单元的具体步骤为：

步骤S41：通过直角坐标系中均匀分布的X方向波数域值k_x、T方向波数域值k_y计算对应的极坐标系下的斜距向波数k，及角度θ；

步骤S42：通过二维插值将极坐标系下的数据转换到直角坐标系，生成第4信号S₄(k_x，k_y)。

5.根据权利要求1所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S6中空变滤波单元的具体步骤为：

步骤S61：将第5信号表示的二维图像按不同半径分成互不重叠的若干分区，分区准则为分区内外圈上的目标，残余二次越距离单元徙动引起的相位误差差值不超过π，表达式为：

$k_{\max }{\sin }^2\mathrm{\&alpha;}({r_{l+1}}^2-{r_l}^2)/\left({2R}_0\right)<\mathrm{\&pi;}$

其中，k_max为最大斜距向波数，r_l+1，r_l分别为第l个分区的外圈半径和内圈半经，π为圆周率。

步骤S62：各分区分别作二维傅里叶变换，得到各分区二维频谱；

步骤S63：各分区频谱与各自的相位补偿滤波器H₆(l；k_x，k_y)相乘，生成第63信号S₆₃(l；k_x，k_y)；

步骤S64：将第63信号各自作傅立叶逆变换，回到图像域，并拼接为第6信号S₆(z₀；x，y)；

6.根据权利要求1所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，圆迹合成孔径雷达的成像几何中，雷达发射脉冲的角度间隔Δθ应小于或等于

其中R为雷达平台轨道半径，H为雷达平台高度，X₀为场景照射区半径。

7.根据权利要求1所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S2中成像高度间隔Δz应小于或等于

其中k_min为最小斜距向波数。

8.根据权利要求3所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S32中方位二次滤波函数H₃(k_θ，k)的表达式为：H₃(k_θ，k)＝exp[jk_θ ²/(2kR₀)]。

9.根据权利要求5所述的一种圆迹合成孔径雷达的大场景极坐标格式三维成像方法，其特征在于，所述步骤S63中相位补偿滤波器H₆(l；k_x，k_y)的表达式为：

Images