CN106772380A - 一种圆周合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开了一种圆周合成孔径雷达成像方法，属于电子信号处理技术领域，涉及空间遥感和空对地观测信息处理技术，特别涉及机载圆周合成孔径雷达成像技术。本发明利用求解系统核函数的逆，将原始斜平面回波映射到地平面，再通过方位角频域与参考函数相乘得到信号在波数域的分布，避免了平面波近似，解决了传统极坐标格式算法成像区域尺寸大小受限的问题，可以实现大成像场景的成像；采用伪极坐标代替直角坐标作为中间插值过渡矩阵，考虑了信号在波数域的分布密度特性，插值精度更高，得到的图像分辨率也更高，可实现高分辨成像；成像过程中仅涉及一维插值，且采用了CZT及IFFT等快速算法，具备很高的计算效率，可实现快速成像。

Description

一种圆周合成孔径雷达成像方法

技术领域

本发明属于电子信号处理技术领域，涉及空间遥感和空对地观测信息处理技术，特别涉及一种机载圆周合成孔径雷达(Circular Synthetic Aperture Radar,CSAR)成像技术。

背景技术

SAR是二十世纪雷达技术发展的重要里程碑，它利用雷达回波信号的相关性，累积雷达运动过程中回波信号的多普勒频移，在雷达的运动方向上合成等效的雷达孔径，实现方位向的高分辨成像。因SAR采用主动式工作模式，对自身发射电磁波的反射回波进行成像处理，不受光照、温度等外界环境的限制，可实现全天时、全天候的区域监测成像，且对植被、沙漠覆盖等介质具有穿透能力而在灾害评估、环境监测、海洋观测、资源勘查、植被监测、测绘等领域得到了广泛的应用。

圆周SAR(CSAR)是一种新兴的体制。不同于传统的SAR，它以圆周运动轨迹采集数据代替传统的直线采集，其系统几何示意图如图1所示，雷达平台在高为H的高度面内做匀速圆周运动，其运动轨迹的半径为R_gc，则雷达平台在空间域中的位置可表示为r_s＝(X,Y,Z)＝(R_gccosθ,R_gcsinθ,H)，其中θ∈[0,2π)表示慢时间方位角。当雷达沿圆周轨迹运动时，其波束中心始终指向场景中心，并覆盖地面上以场景中心为圆心、以R₀为半径的圆形场景区域。雷达与场景中心的距离为雷达俯仰角为θ_z＝arctan(H/R_gc)。

CSAR相对于传统SAR具有如下优势：

首先，CSAR由于其特殊的飞行轨迹，使其在距离向和方位向的频谱展到最宽，因此可以获得方位向和距离向最高的分辨率，并且距离向和方位向的分辨率将保持一致，相对于传统聚束SAR在分辨率方面又有了一个新的提高。

其次，传统的条带SAR还是聚束SAR他们都是以直线运动轨迹来采集数据的，导致雷达只能在一个有限的观测角度内观测目标区域内的目标，这样会丧失很多目标特征。而CSAR的观测角度可以达到360度，可以对目标做全向观测，获得目标的全视角特征。

然而，现有的CSAR成像算法主要分为两类：参数化方法和非参数方法。参数化方法是指利用现代谱估计等参数估计方法提取散射点幅度和位置信息，包括基于RELAX的目标特征提取方法以及基于广义Radon变换的成像算法等，但此类算法对于散射点特征不明显的均匀场景成像效果较差。非参数方法是通过信号时域或频域分布利用聚焦方法重建目标函数，主要包括：后向投影(Back-projection,BP)算法，共焦成像算法，基于格林函数分解的成像算法等。BP算法将成像区域划分为网格(像素点)，计算每一像素点对应的积分路径，然后沿该路径进行积分(需要插值)，即将信号投影到对应的像素点，完成该点的聚焦。共焦成像算法与BP算法类似，也需计算每一像素点的路径，再利用空变滤波得到聚焦结果。以上两种算法均需对图像每一像素点进行逐点计算，因而计算量巨大，不适用于大场景高分辨的实际应用。而基于格林函数分解的成像算法虽较前面两种算法效率更高，但是算法中涉及到极坐标信号到直角坐标的二维插值，导致成像精度受插值算法的影响，且会引入额外的插值误差，进而影响成像效果。

发明内容

本发明的发明目的在于：针对上述存在的问题，公开了一种圆周合成孔径雷达成像方法，针对高分辨的应用条件，在CSAR成像过程中充分考虑信号的极坐标格式分布特性，通过算法改进，解决传统极坐标到直角坐标二维插值带来的成像误差的同时，采用傅里叶变换和快速分数阶傅里叶变换等快速算法提高算法效率，是一种高效高分辨圆周合成孔径雷达成像方法。

本发明的一种圆周合成孔径雷达成像方法包括下列步骤：

步骤S1：对CSAR原始回波进行距离向处理，即对CSAR原始回波沿距离向变换到距离频域；

再将距离向处理结果与距离向参考信号相乘，得到距离向匹配滤波后的斜平面信号S₁(ω,θ)，其中ω为快时间角频率，θ为雷达方位角，斜平面信号S₁(ω,θ)也称快时间频率-方位角域信号；

步骤S2：对斜平面信号S₁(ω,θ)进行地平面转换，得到地平面信号S₂(ω_g,θ)，即基于系统核函数Λ(ω,ω_g)与逆核函数Λ^-1(ω_g,ω)，将信号S₁(ω,θ)转换为地平面信号S₂(ω_g,θ)，其中ω_g为地平面快时间频率；

步骤S3：将地平面信号S₂(ω_g,θ)沿角度向进行傅里叶变换得到地平面快时间频率-方位角频域信号S₃(ω_g,ξ)，其中ξ为对应雷达方位角θ的方位角频率域；

再将S₃(ω_g,ξ)与方位向参考信号S_g0(ω_g,θ)在方位角频域进行匹配滤波，得到极坐标系下的空间频率谱F_p(ρ,θ)，其中ρ为距离向空间频率；

步骤S4：对空间频率谱F_p(ρ,θ)进行频率归一化，在-π＜ρ≤π的限制下，得到归一化后的极坐标系下的空间频谱F_p0(ρ,θ)；

对空间频谱F_p0(ρ,θ)进行角度向插值后，再进行径向插值，得到伪极坐标下的空间频率谱F_pp(ξ_x,ξ_y)，所述F_pp(ξ_x,ξ_y)点数为2N×2N，由竖直子频谱F_v(ξ_x,ξ_y)和水平子频谱F_h(ξ_x,ξ_y)组成：

对于F_v(ξ_x,ξ_y)，频点分布为：

对于F_h(ξ_x,ξ_y)，频点分布为：

其中d表示径向索引，且-N≤d＜N；m表示角度向索引，且

步骤S5：对伪极坐标系下的空间频谱F_pp(ξ_x,ξ_y)的竖直子频谱F_v(ξ_x,ξ_y)、水平子频谱F_h(ξ_x,ξ_y)分别进行伪极坐标成像处理，得到最终的成像结果。

本发明的有益效果是：利用伪极坐标代替直角坐标作为中间插值矩阵，减小误差的同时也减小了计算量，本发明可操作性强，效率高；利用求解系统核函数的逆，将原始斜平面回波映射到地平面，再通过方位角频域与参考函数相乘得到信号在波数域的分布，避免了平面波近似，解决了传统极坐标格式算法成像区域尺寸大小受限的问题，可以实现大成像场景的成像；采用伪极坐标代替直角坐标作为中间插值过渡矩阵，考虑了信号在波数域的分布密度特性，插值精度更高，因而得到的图像分辨率也更高，可实现高分辨成像。

附图说明

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

图2是本发明的处理流程图；

图3是伪极坐标频谱分布示意图；

图4是极坐标系到伪极坐标系转换过程中角度向插值示意图；

图5是极坐标系到伪极坐标系转换过程中径向插值示意图；

图6是伪极坐标傅里叶变换成像仿真结果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

参见图2，通过本发明的一种圆周合成孔径雷达的高效高分辨成像方法，能将输入的输入为CSAR的原始回波，经过成像处理后，获得CSAR二维图像，具体实施步骤如下：

步骤S1：对CSAR原始回波进行距离向处理：沿距离向变换到距离频域，并与距离向参考信号相乘，得到距离向匹配滤波后的斜平面快时间频率-方位角域信号S₁(ω,θ)，其中ω为快时间角频率，θ为雷达方位角，S₁(ω,θ)表达式为：

S₁(ω,θ)＝∫∫f(x,y)g_θ(x,y,ω)dxdy

其中f(x,y)为目标的散射系数，为CSAR斜平面格林函数，k＝ω/c为波数，虚数单位下同；

步骤S2：将斜平面信号S₁(ω,θ)转换到地平面信号，得到地平面信号S₂(ω_g,θ)，其中ω_g为地平面快时间频率；

步骤S21：计算系统核函数

其中，θ_x＝arcsin(R₀/R_gc)，ρ为距离向空间频率；

步骤S22：对核函数求逆，得到逆核函数Λ^-1(ω_g,ω)，再将逆核函数Λ^-1(ω_g,ω)与斜平面信号S₁(ω,θ)关于斜平面快时间频率ω进行卷积(利用FFT实现)，得到地平面信号S₂(ω_g,θ)，公式为：

S₂(ω_g,θ)＝∫_ωΛ^-1(ω_g,ω)S₁(ω,θ)dω

步骤S3：对地平面信号S₂(ω_g,θ)进行极坐标角度向处理，与方位角频域的角度向参考信号S_g0(ω_g,ξ)进行匹配滤波，得到极坐标系下的目标特性函数的空间频率谱F_p(ρ,θ)：

步骤S31：计算方位向参考信号S_g0(ω_g,ξ)，其表达式为：

其中，表示ξ阶第二类Hankel函数。

步骤S32：将地平面信号S₂(ω_g,θ)沿角度向进行傅里叶变换，得到地平面快时间频率-方位角频域的信号S₂(ω_g,ξ)，与角度参考信号S_g0(ω_g,ξ)相乘，得到方位角频域的空间频谱F_p(ρ,ξ)，其表达式为：

其中，表示ξ阶第一类Hankel函数。

步骤S33：对F_p(ρ,ξ)沿角度向进行傅里叶逆变换，得到极坐标系下的目标特性函数的空间频率谱F_p(ρ,θ)；

步骤S4：对极坐标系下的空间频率谱F_p(ρ,θ)进行坐标转换处理，得到伪极坐标下的空间频率谱F_pp(ξ_x,ξ_y)：

其中F_pp(ξ_x,ξ_y)的点数为2N×2N，由竖直子频谱F_v(ξ_x,ξ_y)和水平子频谱F_h(ξ_x,ξ_y)组成，其频谱分布示意图参见图3。为简化叙述，以下均以竖直子频谱F_v(ξ_x,ξ_y)为例进行讨论，对于水平子频谱F_h(ξ_x,ξ_y)，其原理与F_v(ξ_x,ξ_y)相同；

步骤S41：对空间频率谱F_p(ρ,θ)进行频率归一化，在-π＜ρ≤π的限制下，得到归一化后的极坐标系下的空间频谱F_p0(ρ,θ)；

步骤S42：对极坐标系下归一化后的空间频谱F_p0(ρ,θ)先进行角度向插值，得到等斜率间隔的极坐标频谱F(ρ,θ_p)，过程示意图如图4所示，其中θ_p为伪极坐标下的方位角，呈现为等斜率间隔分布。θ_p(m)为θ_p的离散形式，对于F_v(ξ_x,ξ_y)有：

步骤S43：对角度向插值后的频谱F(ρ,θ_p)进行径向插值，得到同心矩形分布的伪极坐标频谱F_pp(ρ_p,θ_p)。ρ_p(d,m)为伪极坐标下第m个角度上第d个频点的半径，对于角度向插值后的F_v(ξ_x,ξ_y)有：

对于单个方位角θ_p(m)，按上式即可得到该角度下沿同心矩形分布的径向频谱值；对每一个方位角索引m，重复以上方法即可得到伪极坐标分布的频谱F_pp(ρ_p,θ_p)＝F_pp(ξ_x,ξ_y)，插值过程示意图如图5所示。

步骤S5：对F_pp(ξ_x,ξ_y)进行伪极坐标成像处理，得到最终的成像结果f(k₁,k₂)＝f(x,y)：

步骤S51：分别对F_pp(ξ_x,ξ_y)的竖直子频谱F_v(ξ_x,ξ_y)(也可表示为F_v(m,d)，下标v用于标识竖直方向)和水平子频谱F_h(ξ_x,ξ_y)(也可表示为F_h(m,d)，下标h用于标识水平方向)的每一径向单元，沿角度向进行α＝-d/N的CZT，得到角度向处理后的子频谱和以子频谱为例，其公式为：

步骤S52：对角度向处理后的子频谱和分别沿径向进行快速傅里叶逆变换(IFFT)，得到子频谱图像f_v(k₁,k₂)和f_h(k₁,k₂)，k₁和k₂为图像坐标，公式如下：

步骤S53：将子频谱图像f_v(k₁,k₂)和f_h(k₁,k₂)相加，得到最终形成的图像f(k₁,k₂)，即：

f(x,y)＝f_v(k₁,k₂)+f_h(k₁,k₂)

实施例

采用机载雷达参数进行成像仿真，载机的标准飞行高度3000m，雷达与场景中心俯仰角为θ＝30°，中心频率9.6GHz，带宽1.2GHz，方位向采样间隔0.05°。设置了10个点目标，散射强度相同，水平和竖直方向间隔均为0.3m。采用本发明得到的成像结果如图6所示。由成像结果可见利用本发明的基于伪极坐标傅里叶变换的成像结果可以实现高效高分辨率成像。

Claims (3)

1.一种圆周合成孔径雷达成像方法，其特征在于，包括下列步骤：

步骤S1：对CSAR原始回波进行距离向处理，即对CSAR原始回波沿距离向变换到距离频域；

再将距离向处理结果与距离向参考信号相乘，得到距离向匹配滤波后的斜平面信号S₁(ω,θ)，其中ω为快时间角频率，θ为雷达方位角；

步骤S2：对斜平面信号S₁(ω,θ)进行地平面转换，得到地平面信号S₂(ω_g,θ)，即基于系统核函数Λ(ω,ω_g)与逆核函数Λ^-1(ω_g,ω)，将信号S₁(ω,θ)转换为地平面信号S₂(ω_g,θ)，其中ω_g为地平面快时间频率；

步骤S3：将地平面信号S₂(ω_g,θ)沿角度向进行傅里叶变换得到地平面快时间频率-方位角频域信号S₃(ω_g,ξ)，其中ξ为对应雷达方位角θ的方位角频率域；

再将S₃(ω_g,ξ)与方位向参考信号S_g0(ω_g,θ)在方位角频域进行匹配滤波，得到极坐标系下的空间频率谱F_p(ρ,θ)，其中ρ为距离向空间频率；

步骤S4：对空间频率谱F_p(ρ,θ)进行频率归一化，在-π＜ρ≤π的限制下，得到归一化后的空间频谱F_p0(ρ,θ)；

对归一化后的空间频谱F_p0(ρ,θ)进行角度向插值后，再进行径向插值，得到伪极坐标下的空间频率谱F_pp(ξ_x,ξ_y)，所述F_pp(ξ_x,ξ_y)点数为2N×2N，由竖直子频谱F_v(ξ_x,ξ_y)和水平子频谱F_h(ξ_x,ξ_y)组成：

对于F_v(ξ_x,ξ_y)，频点分布为：

对于F_h(ξ_x,ξ_y)，频点分布为：

其中d表示径向索引，且-N≤d＜N；m表示角度向索引，且

步骤S5：对空间频谱F_pp(ξ_x,ξ_y)的竖直子频谱F_v(ξ_x,ξ_y)、水平子频谱F_h(ξ_x,ξ_y)分别进行伪极坐标成像处理，得到最终的成像结果。

2.如权利要求1所述的方法，其特征在于，步骤S4中，对归一化后的空间频谱F_p0(ρ,θ)进行角度向插值具体为：

对空间频率谱F_p0(ρ,θ)进行角度向插值，得到等斜率间隔的极坐标频谱F(ρ,θ_p)，其中θ_p(m)为伪极坐标下的方位角θ_p的离散形式，对于竖直子频谱有：

${\mathrm{\&theta;}}_p\left(m\right)=\arctan \left(\frac{{\mathrm{\&xi;}}_y}{{\mathrm{\&xi;}}_x}\right)=\arctan \left(\frac{N}{2m}\right),-\frac{N}{2}\&le;m<\frac{N}{2};$

利用Sinc函数对待插值点信号进行重建，公式如下：

$F(\mathrm{\&rho;},{\mathrm{\&theta;}}_p(m\left)\right)=\underset{d=-M}{\overset{M-1}{\&Sigma;}}{F}_{p0}(\mathrm{\&rho;},\frac{\mathrm{\&pi;}d}{M})\frac{sin\&lsqb;\frac{1}{2}(2M-1)\left({\mathrm{\&theta;}}_p\right(m)-\frac{\mathrm{\&pi;}d}{M})\&rsqb;}{2M\sin \&lsqb;\frac{1}{2}\left({\mathrm{\&theta;}}_p\right(m)-\frac{\mathrm{\&pi;}d}{M})\&rsqb;},$

其中，2M为插值核点数；

再对角度向插值得到的等斜率间隔的极坐标频谱F(ρ,θ_p)采用三次样条插值算法进行径向插值，得到同心矩形分布的伪极坐标频谱F_pp(ρ_p,θ_p)＝F_pp(ξ_x,ξ_y)，ρ_p(d,m)为伪极坐标下第m个角度上第d个频点的半径，对于竖直子频谱有：

${\mathrm{\&rho;}}_p(d,m)=\sqrt{{\mathrm{\&xi;}}_x^2+{\mathrm{\&xi;}}_y^2}=\frac{\mathrm{\&pi;}d}{N}\sqrt{1+{\left(\frac{2m}{N}\right)}^2},-N\&le;d<N.$

3.如权利要求1或2所述的方法，其特征在于，步骤S5中的伪极坐标成像处理具体为：

步骤S51：分别对竖直子频谱F_v(ξ_x,ξ_y)、水平子频谱F_h(ξ_x,ξ_y)的每一径向单元，沿角度向进行α＝-d/N的Chirp-Z变换，得到角度向处理后的子频谱和其中k₁表示图像水平方向索引；

步骤S52：对角度向处理后的子频谱和分别沿径向进行快速傅里叶逆变换，得到子频谱图像f_v(k₁,k₂)和f_h(k₁,k₂)，其中k₂为图像竖直方向索引；

步骤S53：将子频谱图像f_v(k₁,k₂)和f_h(k₁,k₂)相加，得到最终的成像结果。