CN103728619B - 基于变重频技术的机载大斜视条带sar成像方法 - Google Patents

Abstract

本发明公开了一种基于变重频技术的机载大斜视条带SAR成像方法，包括步骤一：读入原始回波数据及相关成像参数；步骤二：读入变重频相关参数；步骤三：距离向位置补偿以及距离走动补偿处理；步骤四：方位向数据频谱搬移处理；步骤五：方位向拉格朗日插值处理；步骤六：方位向数据频带恢复处理；步骤七：二维傅里叶变换处理；步骤八：一致压缩处理；步骤九：斯托尔特插值处理；步骤十：方位向傅里叶逆变换处理；步骤十一：几何校正处理；步骤十二：距离向傅里叶逆变换处理；本发明由于本发明提出的成像算法利用拉格朗日插值法实现方位向非均匀采样重构，相比于非均匀傅里叶变换，能够更加简单容易地完成方位向非均匀信号的有效重构。

Description

基于变重频技术的机载大斜视条带SAR成像方法

技术领域

本发明属于信号处理领域，特别涉及一种基于变重频技术的机载大斜视条带合成孔径雷达(SyntheticApertureRadar，SAR)成像方法。

背景技术

SAR作为一种重要的微波成像技术，能够克服云雾雨雪和暗夜条件的限制进行对地遥感观测，在军事侦察测绘以及诸多民用领域可发挥重要作用，为越来越多的国家所重视。机载斜视SAR由于其具有较高的灵活性和机动性，能够灵活地选择观测区域与快速重访敏感区域，使得斜视SAR系统观测性能优于传统正侧视SAR系统。因此，近些年来机载斜视SAR已经成为一个重要的发展方向。

但是，机载大斜视SAR系统也面临新困难与挑战。由于距离徙动量随斜视角的增大以及分辨率的提高而成倍增长，进而使得回波接收窗长逐渐增大。此外，同等分辨率条件下斜视SAR系统方位有效带宽大于正侧视SAR系统的方位有效带宽，且随斜视角度的增大，方位向有效带宽逐渐增大，这意味着为避免方位向频谱混叠，斜视SAR系统脉冲重复频率需随之增大，即脉冲重复周期变小。因此，随着斜视角度的增大，回波接收窗长增大而脉冲重复周期变小，当斜视角度增加到一定程度时，脉冲重复周期小于回波接收窗长将导致大斜视SAR系统无法正常工作。除此之外，大斜视SAR系统的大距离徙动量使得原始回波数据量的显著增加，加大了数据存储与实时成像的困难。针对上述大斜视SAR系统面临的问题，目前主要采用变脉冲重复频率技术，即变重频技术。变重频技术通过改变脉冲重复频率，调整回波在接收窗中的位置，进而缩短回波接收窗所需长度以及回波数据量，从而使得大斜视SAR系统正常工作。由于距离走动在距离徙动中的起主导作用且距离走动与斜距无关，因此常用的变重频技术准则为通过变脉冲重复频率消除距离走动。

变重频技术的出现解决了大斜视SAR系统回波录取以及回波海量数据的问题，但是基于变重频技术的机载大斜视条带SAR系统数据处理存在以下挑战：方位向与距离向严重耦合、方位向非均匀采样重构以及距离徙动精确校正，上述问题使得传统成像方法不再适用。因此，发明一种基于变重频技术的机载大斜视条带SAR成像方法显得至关重要。

发明内容

本发明的目的是为了解决上述问题，基于变重频技术的特点，结合传统波数域成像方法，提出了一种基于变重频技术的机载大斜视条带SAR成像方法。

一种基于变重频技术的机载大斜视条带SAR成像方法，包括以下几个步骤：

步骤一：读入原始回波数据及相关成像参数；

读入基于变重频技术的机载大斜视条带SAR二维原始回波仿真复数组S_start以及相应的成像参数，具体包括：方位向采样点数N_a，距离向采样点数N_r，信号采样率f_s，信号带宽Bw，脉冲宽度τ，调频斜率b，仿真中心时刻脉冲重复频率PRF_c，参考斜距R_ref，多普勒中心频率fd₀，多普勒调频率f_r0，卫星速度P_v，等效斜视角信号波长λ，信号载频f₀，信号传播速度c；

步骤二：读入变重频相关参数；

读入变重频相关参数，具体包括：一维方位向各采样时刻脉冲重复频率序列PRF_cha，一维方位向各采样时刻时间序列T_cha以及一维方位向各采样时刻瞬时斜距序列R_cha。

步骤三：距离向位置补偿以及距离走动补偿处理；

将二维回波仿真复数组S_start做距离向位置补偿处理：首先计算按仿真中心时刻脉冲重复频率PRF_c进行均匀采样时的一维方位向各采样时刻序列T_uni；其次计算按均匀采样时刻T_uni时的一维瞬时斜距序列R_uni；再次对二维回波仿真复数组S_start做距离向傅里叶变换，即沿方位向（按行）进行快速傅里叶变换（FFT），得到二维复数组S_{loc_1}；再次将二维复数组S_{loc_1}同距离向位置补偿因子H₁以及距离走动补偿因子H₂相乘，得到二维复数组S_{loc_2}；最后将二维复数组S_{loc_2}做距离向傅里叶逆变换，即沿方位向（按行）进行快速傅里叶逆变换（IFFT），得到调整距离向位置后的二维复数组S_loc；

步骤四：方位向数据频谱搬移处理；

将步骤三得到的二维复数组S_loc同频谱搬移因子H₃相乘，得到方位向频谱搬移后二维复数组S_base；

步骤五：方位向拉格朗日插值处理；

对步骤四得到的二维复数组S_base，利用拉格朗日插值法进行方位向插值处理，得到按PRF_c均匀采样的二维复数组S_uni；

步骤六：方位向数据频带恢复处理；

对步骤五得到的二维复数组S_uni同频带恢复因子H₄相乘，得到方位向频带恢复后二维复数组S_new；

步骤七：二维傅里叶变换处理；

对步骤六得到的二维复数组S_new做二维（方位向与距离向）傅里叶变换：首先沿方位向（按行）进行快速傅里叶变换（FFT），得到方位时域距离频域复数组S_sp1；其次沿距离向（按列）进行快速傅里叶变换（FFT），得到二维频谱复数组S_sp；

步骤八：一致压缩处理；

将步骤七得到的二维频谱复数组S_sp同一致压缩因子H₅相乘，得到粗聚焦复数组S_bulk；

步骤九：斯托尔特插值处理；

将步骤八得到的粗聚焦复数组S_bulk，利用辛格插值法进行斯托尔特插值处理，得到精确聚焦的二维复数组S_wave；

步骤十：方位向傅里叶逆变换处理；

将步骤九得到的复数组S_wave沿距离向（按列）进行快速傅里叶逆变换（IFFT），得到方位时域距离频域复数组S_a-t；

步骤十一：几何校正处理；

将步骤十得到的方位时域距离频域复数组S_a-t同几何校正因子H₆相乘，得到经几何校正后的复数组S_geo；

步骤十二：距离向傅里叶逆变换处理；

将步骤十一得到的复数组S_geo沿方位向（按行）进行快速傅里叶逆变换（IFFT），得到最终的成像结果S_end；

本发明优点在于：

（1）本发明提出了一种基于变重频技术的机载大斜视条带SAR成像方法，解决了目前基于变重频技术的机载大斜视条带SAR原始回波数据没有成像算法的现状。

（2）本发明提出了一种基于变重频技术的机载大斜视条带SAR成像方法，具有简单易实现的特点。由于本发明提出的成像算法利用拉格朗日插值法实现方位向非均匀采样重构，相比于非均匀傅里叶变换，能够更加简单容易地完成方位向非均匀信号的有效重构。

（3）本发明提出了一种基于变重频技术的机载大斜视条带SAR成像方法，具有运行效率高的特点。由于本发明提出的成像算法在成像开始补偿变重频技术带来的回波录取位置变化的同时做距离走动补偿，减小了成像时所需距离向长度，减小了成像数据量，进而提高了成像效率。

（4）本发明提出了一种基于变重频技术的机载大斜视条带SAR成像方法，具有高精度聚焦的特点。由于本发明提出的成像方法是一种基于波数域聚焦原理的成像方法，而波数域成像方法只要满足平台速度恒定这一条件（机载SAR恰好满足这一条件），就能够实现高精度聚焦。

（5）本发明提出了一种基于变重频技术的机载大斜视条带SAR成像方法，具有适用性强的特点。由于本发明提出的成像方法是一种基于波数域聚焦原理的成像方法，而波数域成像方法不受斜视角度的限制，因此，在斜视角度很大的条件下，本发明同样能实现场景的精确聚焦。

附图说明

图1是本发明提出的一种基于变重频技术的机载大斜视条带SAR成像方法流程图；

图2是本发明步骤三的方法流程图；

图3是实施例仿真场景示意图；

图4是实施例成像结果；

图5是实施例左上点目标剖面图；

图6是实施例中间点目标剖面图；

图7是实施例右下点目标剖面图；

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明提出了一种基于变重频技术的机载大斜视条带SAR成像方法，处理的对象是基于变重频技术的机载大斜视条带SAR原始回波数据，得到的结果是一幅大斜视SAR图像。

基于消除距离走动的变重频技术通过改变脉冲重复频率，调整了场景回波在接收窗的位置，从而达到消除距离走动，减少回波数据量的目的。但是，脉冲重复频率的改变导致场景回波在方位向存在非均匀采样的现象。此外，场景回波在接收窗位置的改变也导致目标多普勒历程的改变，加之大斜视SAR系统本身存在的方位向/距离向严重耦合，使得现有成像算法无法实现基于变重频技术的机载大斜视条带SAR系统回波精确成像。本发明通过方位向拉格朗日插值法解决了方位向非均匀采样现象，并通过改进传统的波数域成像算法，实现了场景目标的精确聚焦。

本发明是一种基于变重频技术的机载大斜视条带SAR成像方法，具体流程如图1所示，包括以下步骤：

步骤一：读入基于变重频技术的机载大斜视条带SAR二维原始回波仿真复数组S_start以及相应的成像参数。其中，S_start是大小为N_a×N_r二维复数组，而成像参数具体包括：方位向采样点数N_a，距离向采样点数N_r，信号采样率f_s，信号带宽Bw，脉冲宽度τ，调频斜率b，仿真中心时刻脉冲重复频率PRF_c，参考斜距R_ref，多普勒中心频率fd₀，多普勒调频率f_r0，卫星速度P_v，等效斜视角，信号波长λ，信号载频f₀，信号传播速度c；

步骤二：读入变重频相关参数。具体包括：一维方位向各采样时刻脉冲重复频率序列PRF_cha，一维方位向各采样时刻时间序列T_cha以及一维方位向各采样时刻瞬时斜距序列R_cha。其中，一维数组PRF_cha、T_cha和R_cha大小均为N_a×1；

步骤三：将二维原始回波仿真复数组S_start做距离位置补偿以及距离走动补偿处理，操作流程如图2所示，具体步骤为：

（a）构造一维序列M,N，其中M代表行，N代表列；

M＝[1，2，...，N_a]

N＝[1,2,...,N_r]（1）

（b）计算按仿真中心时刻脉冲重复频率PRF_c进行均匀采样时的一维方位向各采样时刻序列T_uni(M)，并计算均匀采样时刻下的一维瞬时斜距序列R_uni(M)；

$T_{\mathrm{uni}}\left(M\right)=\frac{M-N_a/2}{{\mathrm{PRF}}_c}---\left(2\right)$

（c）计算按PRF_cha采样时场景中心目标回波接收时刻与回波窗开启时刻的时延t_c(M)，并计算按PRF_c采样时场景中心目标回波接收时刻与回波窗开启时刻的时延t_u(M)；

其中，为向下取整函数；

（d）对二维回波仿真复数组S_start(M,N)做距离向傅里叶变换，即沿方位向（按行）进行快速傅里叶变换（FFT），得到二维复数组S_{loc_1}(M,N)；

S_{loc_1}(M,:)＝FFT(S_start(M,:))（6）

其中，S_start(M,:)表示S_start的第M行，S_{loc_1}(M,:)表示S_{loc_1}的第M行，FFT(·)表示对一维数组进行快速傅里叶变换。

（e）计算距离频域每列对应的距离频率数组f_τ(N)；

$f_{\mathrm{\τ}}\left(N\right)=\frac{N-N_r/2}{N_r}\·f_s---\left(7\right)$

（f）计算距离向位置补偿因子H₁(M,N)；

$H_1(M,N)=\exp \{-j\·2\mathrm{\π}\·[(\frac{2R_{\mathrm{cha}}\left(M\right)}{c}-\frac{2R_{\mathrm{uni}}\left(M\right)}{c})-(t_c\left(M\right)-t_u\left(M\right))]\·f_{\mathrm{\τ}}\left(N\right)\}---\left(8\right)$

其中j为单位复数；

（g）计算距离走动补偿因子H₂(M,N)；

H₂(M,N)＝exp{j·2π·f_τ(N)·[λfd₀T_cha(M)]}（9）

（h）将二维复数组S_{loc_1}(M,N)同距离向位置补偿因子H₁(M,N)与距离走动补偿因子H₂(M,N)相乘，得到二维复数组S_{loc_2}(M,N)；

S_{loc_2}(M,N)＝S_{loc_1}(M,N)·H₁(M,N)·H₂(M,N)（10）

（i）将二维复数组S_{loc_2}(M,N)做距离向傅里叶逆变换，即沿方位向（按行）进行快速傅里叶逆变换（IFFT），得到调整距离向位置后的二维复数组S_loc(M,N)；

S_loc(M,:)＝IFFT(S_{loc_2}(M,:))（11）

其中，S_loc(M,:)表示S_start的第M行，S_{loc_2}(M,:)表示S_{loc_1}的第M行，IFFT(·)表示对一维数组进行快速傅里叶逆变换。

步骤四：方位向数据频谱搬移处理，具体可以分为以下几个步骤：

（a）构造方位向频谱搬移因子H₃(M)；

H₃(M)＝exp{j·2π·fd₀·T_cha(M)}（12）

（b）将二维复数组S_loc(M,N)同方位向频谱搬移因子H₃(M)相乘，得到方位向频谱搬移后的二维复数组S_base(M,N)；

S_base(M,N)＝S_loc(M,N)·H₃(M)（13）

步骤五：将步骤四得到的方位向频谱搬移后的二维复数组S_base做方位向拉格朗日插值处理，具体可以分为以下几个步骤：

（a）按列N计算均匀采样时刻序列T_uni(M)在不均匀采样时刻序列T_cha(M)中的位置p(M,N)。具体操作为：以计算第一行第一列的位置p(1,1)为例，首先计算绝对差值Δt(n)＝|T_uni(1)-T_cha(n)|,n＝[1,2,...,N_a]，获取最小的绝对差值Δt_min和对应位置的n，若T_uni(1)＜T_cha(n)，p(1,1)＝n，若T_uni(1)≥T_cha(1)，p(1,1)＝n+1，以此类推，得到每一个位置p(M,N)。

（b）计算对二维均匀采样数组S_uni(M,N)进行拉格朗日插值所需的插值基函数L(M,N)；

$\begin{array}{c}L(M,N)=\underset{\underset{M\&NotEqual;p}{p=0}}{\overset{N_l}{\mathrm{\&Pi;}}}\frac{T_{\mathrm{uni}}\left(M\right)-T_{\mathrm{cha}}[p(M,N)-N_l/2+p]}{T_{\mathrm{cha}}[p(M,N)-N_l/2+M]-T_{\mathrm{cha}}[p(M,N)-N_l/2+p]}\\ N_l/2<\mathrm{pos}(M,N)<N_a-N_l/2\end{array}---\left(14\right)$

其中，N_l为拉格朗日插值阶数。

（c）结合拉格朗日插值基函数L(M,N)，分别对S_base(M,N)的实部S_{base_re}(M,N)和虚部S_{base_im}(M,N)分别进行拉格朗日插值法计算，得到方位向均匀采样的二维复数组S_uni(M,N)；

p(M,N)＜N_l/2,S_{uni_re}(M,N)＝S_{base_re}(M,N)

$N_l/2\&le;p(M,N)\&le;N_a-N_l/2,S_{\mathrm{uni}\_\mathrm{re}}(M,N)=\underset{k=0}{\overset{N_l}{\mathrm{\&Sigma;}}}L(p(M,N)-N_l/2+k,N)\&CenterDot;S_{\mathrm{base}\_\mathrm{re}}(p(M,N)-N_l/2+k,N)---\left(15\right)$

p(M,N)＞N_a-N_l/2,S_{uni_re}(M,N)＝S_{base_re}(M,N_a)

p(M,N)＜N_l/2,S_{uni_im}(M,N)＝S_{base_im}(M,N)

$N_l/2\&le;p(M,N)\&le;N_a-N_l/2,S_{\mathrm{uni}\_\mathrm{im}}(M,N)=\underset{k=0}{\overset{N_l}{\mathrm{\&Sigma;}}}L(p(M,N)-N_l/2+k,N)\&CenterDot;S_{\mathrm{base}\_\mathrm{im}}(p(M,N)-N_l/2+k,N)---\left(16\right)$

p(M,N)＞N_a-N_l/2,S_{uni_im}(M,N)＝S_{base_im}(M,N_a)

其中，S_{base_re}(M,N)是指二维复数组S_base第M行第N列数据的实部，S_{base_im}(M,N)是指二维复数组S_base第M行第N列数据的虚部，S_{uni_re}(M,N)是指二维数组S_uni第M行第N列数据的实部，S_{uni_im}(M,N)是指二维复数组S_uni第M行第N列数据的虚部。

步骤六：方位向数据频带恢复处理，具体可以分为以下几个步骤：

（a）构造方位向频带恢复因子H₄(M)；

H₄(M)＝exp{-j·2π·fd₀·T_uni(M)}（17）

（b）将二维复数组S_uni(M,N)同方位向频带恢复因子H₄(M)相乘，得到方位向频带恢复后的二维复数组S_new(M,N)；

S_new(M,N)＝S_base(M,N)·H₄(M)（18）

步骤七：将步骤六得到的复数组S_new做二维傅里叶变换，具体可以分为以下几个步骤：

（a）将步骤六得到的S_new(M,N)做距离向傅里叶变换，即沿方位向（按行）进行快速傅里叶变换（FFT），得到方位时域距离频域复数组S_sp1(M,N)；

S_sp1(M,:)＝FFT(S_new(M,:))（19）

其中，S_sp1(M,:)表示S_sp1的第M行，S_new(M,:)表示S_new的第M行。

（b）将方位时域距离频域复数组S_sp1(M,N)做方位向傅里叶变换，即沿距离向（按列）进行快速傅里叶变换（FFT），得到二维频域复数组S_sp(M,N)；

S_sp(:,N)＝FFT(S_sp1(:,N))（20）

其中，S_sp(:,N)表示S_sp的第N列，S_sp1(:,N)表示S_sp1的第N列。

步骤八：将步骤七得到的二维频域复数组S_sp(M,N)同一致压缩因子H₅(M,N)相乘，得到粗聚焦复数组S_bulk(M,N)，具体步骤如下：

（a）根据参考斜距R_ref计算最近斜距R_min；

$R_{\min }=R_{\mathrm{ref}}-\frac{c}{2f_s}\&CenterDot;\frac{N_r}{2}---\left(21\right)$

（b）计算方位频域每行对应的方位频率一维数组f_a(M)；

$f_a\left(M\right)=\frac{M-N_a/2}{N_a}\&CenterDot;{\mathrm{PRF}}_c---\left(22\right)$

（c）计算二维一致压缩因子H₅(M,N)所需子因子ξ(M,N)；

（d）结合式（21）～式（23）计算大小为N_a×N_r的二维一致压缩因子H₅(M,N)；

（e）结合式（24）计算一致压缩后的粗聚焦数组S_bulk(M,N)；

S_bulk(M,N)＝S_sp(M,N)·H₅(M,N)（25）

步骤九：将步骤八得到的粗聚焦数组S_bulk(M,N)利用辛格插值法进行斯托尔特插值处理，得到精确聚焦二维复数组S_wave(M,N)，具体步骤如下：

（a）计算斯托尔特插值映射距离频率f'_τ(M,N)；

（b）遍历获取斯托尔特插值映射距离频率f'_τ(M,N)的最大值f'_τ,max与最小值f'_τ,min，并计算斯托尔特插值映射距离频率等分间隔Δf'_τ；

${\mathrm{\&Delta;f}}_{\mathrm{\&tau;}}^{\&prime;}=\frac{f_{\mathrm{\&tau;},\max }^{\&prime;}-f_{\mathrm{\&tau;},\min }^{\&prime;}}{N_r}---\left(27\right)$

（c）计算新距离频率均匀一维数组

$f_{\mathrm{\&tau;}}^u\left(N\right)=f_{\mathrm{\&tau;},\min }^{\&prime;}+N\&CenterDot;{\mathrm{\&Delta;f}}_{\mathrm{\&tau;}}^{\&prime;}---\left(28\right)$

（d）计算二维频域复数组新距离频率均匀序列在每行对应的不均匀托尔特插值映射距离频率f'_τ(M,N)中的位置pos(M,N)。方法为按行进行下列操作：以计算第一行第一列的位置pos(1,1)为例，首先计算绝对差值获取最小的绝对差值Δk_min和对应位置的m，若pos(1,1)＝m-1，若pos(1,1)＝m，以此类推，得到每一个位置pos(M,N)。

（e）结合上步得到的位置pos(M,N)，计算辛格插值所需采样点位置q(M,N,n)；

$q(M,N,n)=\frac{f_{\mathrm{\&tau;}}^u\left(N\right)-f_{\mathrm{\&tau;}}^{\&prime;}(M,(\mathrm{pos}(M,N)+n))}{f_s/N_r},n=[-N_s/2,-N_s/2+1,...,N_s/2-1]---\left(29\right)$

其中，N_s是辛格插值核长度。

（f）结合式（29）利用辛格插值法，计算出经斯托尔特插值后的二维复数组S_wave(M,N)，由于S_wave(M,N)是复数组，需要分别对S_wave(M,N)的实部S_{wave_re}(M,N)和虚部S_{wave_im}(M,N)分别进行辛格插值法计算得出。

p(M,N)＜N_s/2,S_{wave_re}(M,N)＝S_{bulk_re}(M,N)

$N_s/2\&le;p(M,N)\&le;N_r-N_s/2,S_{\mathrm{wave}\_\mathrm{re}}(M,N)=\frac{\underset{n=-N_s/2}{\overset{N_s/2-1}{\mathrm{\&Sigma;}}}{S}_{\mathrm{bulk}\_\mathrm{re}}(M,(\mathrm{pos}(M,N)+n))\&CenterDot;\sin c\left(q(M,N,n)\right)}{\underset{n=-N/2}{\overset{N/2-1}{\mathrm{\&Sigma;}}}\sin c\left(q(M,N,n)\right)}---\left(30\right)$

p(M,N)＞N_r-N_s/2,S_{wave_re}(M,N)＝S_{bulk_re}(M,N_r)

p(M,N)＜N_s/2,S_{wave_im}(M,N)＝S_{bulk_im}(M,N)

$N_s/2\&le;p(M,N)\&le;N_r-N_s/2,S_{\mathrm{wave}\_\mathrm{im}}(M,N)=\frac{\underset{n=-N_s/2}{\overset{N_s/2-1}{\mathrm{\&Sigma;}}}{S}_{\mathrm{bulk}\_\mathrm{im}}(M,(\mathrm{pos}(M,N)+n))\&CenterDot;\sin c\left(q(M,N,n)\right)}{\underset{n=-N/2}{\overset{N/2-1}{\mathrm{\&Sigma;}}}\sin c\left(q(M,N,n)\right)}---\left(31\right)$

p(M,N)＞N_r-N_s/2,S_{wave_im}(M,N)＝S_{bulk_im}(M,N_r)

其中，sinc(·)是指插值函数S_{bulk_re}(M,N)是指二维数组S_bulk第M行第N列数据的实部，S_{bulk_im}(M,N)是指二维数组S_bulk第M行第N列数据的虚部，S_{wave_re}(M,N)是指二维数组S_wave第M行第N列数据的实部，S_{wave_im}(M,N)是指二维数组S_wave第M行第N列数据的虚部。

步骤十：将步骤九得到的复数组S_wave沿距离向（按列）进行快速傅里叶逆变换（IFFT），得到方位时域距离频域复数组S_a-t；

S_a-t(:,N)＝IFFT(S_wave(:,N))（32）

其中，S_wave(:,N)表示S_wave的第N列，S_a-t(:,N)表示S_a-t的第N列。

步骤十一：将步骤十得到的方位时域距离频域复数组S_a-t同几何校正因子H₆相乘，得到经几何校正后的复数组S_geo，具体步骤如下：

（a）计算几何校正因子H₆(M)；

$H_6\left(M\right)=\exp \{-j2\mathrm{\&pi;}\&CenterDot;f_a\left(M\right)\&CenterDot;\frac{\mathrm{\&lambda;}\&CenterDot;{\mathrm{fd}}_0\&CenterDot;T_{\mathrm{uni}}\left(M\right)}{c}\}---\left(33\right)$

（b）将方位时域距离频域复数组S_a-t同几何校正因子H₆(M)相乘；

S_geo(M,N)＝S_a-t(M,N)·H₆(M,N)（34）

步骤十二：将步骤十一得到的复数组S_geo沿方位向（按行）进行快速傅里叶逆变换（IFFT），得到最终的成像结果S_end；

S_end(M,:)＝IFFT(S_geo(M,:))（35）

其中，S_geo(M,:)表示S_geo的第M行，S_end(M,:)表示S_end的第M行。

实施例：

本实施例提出了一种基于变重频技术的机载大斜视条带SAR成像方法，仿真场景设置如图3所示，其成像过程中涉及的成像参数如表1所示。

表1实施例参数

本实施例具体包括以下步骤：

步骤一：读入基于变重频技术的机载大斜视条带SAR二维原始回波仿真复数组S_start以及相应的成像参数。其中，S_start是二维复数组，大小为16384×2048，具体成像参数如表1所示；

步骤二：读入变重频相关参数。具体包括：一维方位向各采样时刻脉冲重复频率序列PRF_cha，一维方位向各采样时刻时间序列T_cha以及一维方位向各采样时刻瞬时斜距序列R_cha。其中，一维数组PRF_cha、T_cha和R_cha大小均为16384×1；

步骤三：将二维原始回波仿真复数组S_start距离向位置补偿以及距离走动补偿处理，具体操作步骤为：

（a）构造一维序列，如式（1）所示，M＝[1,2,...,16384],N＝[1,2,...,2048]；

（b）计算按PRF_c进行均匀采样时的一维方位向各采样时刻序列T_uni(M)以及T_uni(M)对应的一维瞬时斜距序列R_uni(M)，具体过程按式（2）与式（3）进行，T_uni(M)和R_uni(M)均是大小为16384×1的一维数组；

（c）计算按PRF_cha采样时场景中心目标回波接收时刻与回波窗开启时刻的时延t_c(M)，并计算按PRF_c采样时场景中心目标回波接收时刻与回波窗开启时刻的时延t_u(M)，具体过程按式（4）与式（5）进行。其中，t_c(M)和t_u(M)均是大小为16384×1的一维数组；

（d）对二维回波仿真数组S_start按式（6）按行进行FFT操作,得到二维复数组S_{loc_1}；

（e）按式（7）计算距离频域每列对应的距离频域f_τ(N)。其中，f_τ(N)是大小为1×2048的一维数组；

（f）结合式（2）～式（5）以及式（7），按式（8）计算距离向位置补偿因子H₁。其中，H₁是大小为16384×2048的二维复数组；

（g）结合式（7），按式（9）计算距离走动补偿因子H₂。其中，H₂是大小为16384×2048的二维复数组；

（h）结合式（8）和式（9），按式（10）计算经距离向位置补偿以及距离走动补偿后的二维复数组S_{loc_2}(M,N)；

（i）对二维复数组S_{loc_2}按式（11）按行进行IFFT操作，得到二维复数组S_loc；

步骤四：将步骤三得到二维复数组S_loc做方位向数据频谱搬移处理，具体操作步骤为：

（a）按式（12）计算方位向频谱搬移因子H₃(M)。其中，H₃(M)是大小为16384×1的二维复数组；

（b）结合式（12），按式（13）计算变换到方位向频谱搬移后的二维复数组S_base；

步骤五：将方位向频谱搬移后的二维复数组S_base做方位向拉格朗日插值处理，具体操作步骤为：

（a）按列N计算均匀采样时刻序列T_uni(M)在不均匀采样时刻序列T_cha(M)中的位置p(M,N)，

方法为：先计算第一行第一列的位置p(1,1)，首先计算绝对差值

Δt(n)＝|T_uni(1)-T_cha(n)|,n＝[1,2,...,16384]，获取最小的绝对差值Δt_min和对应位置的n，若T_uni(1)＜T_cha(n)，p(1,1)＝n，若T_uni(1)≥T_cha(1)，p(1,1)＝n+1，以此类推，得到每一个位置p(M,N)。

（b）按式（14）计算拉格朗日插值所需的插值基函数L(M,N)。其中L(M,N)是大小为16384×2048的二维数组，且选用的拉格朗日插值阶数N_l＝3；

（c）结合式（14）所计算得到的插值基函数，按式（15）和式（16）分别对S_base的实部S_{base_re}和虚部S_{base_im}进行拉格朗日插值，得到方位向均匀采样的二维复数组S_uni；

步骤六：将步骤五得到的二维复数组S_uni做方位向数据频带恢复处理，具体操作步骤为：

（a）按式（17）计算方位向频带恢复因子H₄(M)。其中，H₄(M)是大小为16384×1的二维复数组；

（b）结合式（17），按式（18）计算方位向频带恢复后的二维复数组S_new；

步骤七：将步骤六得到的二维复数组S_new做二维傅里叶变换，具体操作为：

（a）对二维复数组S_new按式（19）按行进行FFT操作,得到方位时域距离频域二维复数组S_sp1；

（b）对方位时域距离频域二维复数组S_sp1按式（20）按列进行FFT操作，得到二维频域复数组S_sp；

步骤八：将步骤七得到的复数组S_sp同对应的一致压缩因子H₅相乘，得到粗聚焦复数组S_bulk，具体操作为：

（a）根据参考斜距R_ref＝47.33km计算最近斜距R_min，具体计算过程按式（21）进行；

（b）按式（22）计算方位频域每行对应的方位频率f_a(M)。其中，f_a(M)是大小为16384×1的一维数组；

（c）结合式（7）与式（22），按式（23）计算二维一致压缩因子H₅(M,N)所需子因子ξ(M,N)。其中，ξ(M,N)是大小为16384×2048的二维数组；；

（d）结合式（21）和式（23），计算大小为16384×2048的二维一致压缩因子H₅(M,N)，具体计算过程按式（24）进行；

（e）结合式（24），按（25）计算一致压缩后的二维复数组S_bulk；

步骤九：对于步骤八得到的粗聚焦数组S_bulk，利用辛格插值法进行斯托尔特插值处理，得到精确聚焦二维复数组S_wave，具体操作为：

（a）按式（26）计算大小为16384×2048的斯托尔特插值映射距离频率f'_τ(M,N)；

（b）遍历获取托尔特插值映射距离频率f'_τ(M,N)的最大值f'_τ,max＝9.01×10⁷Hz与最小值f'_τ，min＝-1.15×10⁸Hz，并按式（27）计算斯托尔特插值映射距离频率等分间隔Δf'_τ；

（c）计算按式（28）计算大小为1×2048的新距离频率均匀序列

（d）计算二维频域复数组新距离频域均匀序列在每行对应的不均匀托尔特插值映射距离频域f'_τ(M,N)中的位置pos(M,N)，其中pos(M,N)是大小为16384×2048的二维数组。具体方法为按行进行下列操作：以计算第一行第一列的位置pos(1,1)为例，首先计算绝对差值获取最小的绝对差值Δk_min和对应位置的m，若，pos(1,1)＝m-1，若pos(1,1)＝m，以此类推，得到每一个位置pos(M,N)。

（e）结合上步得到的位置pos(M,N)，选择辛格插值核长度N_s＝8，按式（29）计算格插值所需采样点位置q(M,N,n)，其中q(M,N,n)是大小为16384×2048×8的三维数组；

（f）结合式（29）利用辛格插值法，计算出经经斯托尔特插值后的二维复数组S_wave，由于二维数据是复数组，需要分别对S_wave的实部S_{wave_re}和虚部S_{wave_im}分别进行辛格插值法计算得出，具体操作按式（30）和式（31）进行。

步骤十：将步骤九得到的复数组S_wave按式（32）按列进行IFFT操作，得到方位时域距离频域复数组S_a-t；

步骤十一：将步骤十得到的复数组S_a-t同几何校正因子H₆相乘，得到经几何校正后的复数组S_geo，具体操作为：

（a）结合式（2）与式（22），计算大小为16384×1的几何校正因子H₆(M)，具体操作按式（33）进行；

（b）结合式（33），按式（34）计算经几何校正后的复数组S_geo；

步骤十二：将步骤十一得到的复数组S_geo按式（35）按行进行IFFT操作，得到最终的成像结果S_end。

经过上述步骤的成像处理，得到最终的场景成像结果如图4所示。表2给出了场景左上、中间、右下三个点目标的成像评估结果，图5、图6、图7分别给出了场景左上、中间、右下三个点目标的二维剖面图。

表二成像评估结果

根据表2评估结果以及图5～图7所示二维剖面图，可以得出：一方面，该成像方法在斜视角度为60度时仍然能够精确聚焦，说明本发明提出的方法大斜视条件下仍能精确聚焦；另一方面，该成像方法对于1m高分辨仍然能够精确聚焦，说明本发明提出的方法能够对目前较高分辨率实现精确聚焦。因此，本发明所提出的方法可以实现基于变重频技术的机载大斜视条带SAR成像精确成像，得到了高精度的成像结果。

Claims (1)

1.一种基于变重频技术的机载大斜视条带SAR成像方法，包括以下步骤：

步骤一：读入基于变重频技术的机载大斜视条带SAR二维原始回波仿真复数组S_start和成像参数；

S_start为N_a×N_r二维复数组，成像参数包括：方位向采样点数N_a，距离向采样点数N_r，信号采样率f_s，信号带宽Bw，脉冲宽度τ，调频斜率b，仿真中心时刻脉冲重复频率PRF_c，参考斜距R_ref，多普勒中心频率fd₀，多普勒调频率f_r0，卫星速度P_v，等效斜视角信号波长λ，信号载频f₀，信号传播速度c；

步骤二：读入变重频相关参数；

变重频相关参数包括一维方位向各采样时刻脉冲重复频率序列PRF_cha，一维方位向各采样时刻时间序列T_cha以及一维方位向各采样时刻瞬时斜距序列R_cha；其中，PRF_cha、T_cha、R_cha大小均为N_a×1；

步骤三：对二维原始回波仿真复数组S_start做距离位置补偿以及距离走动补偿处理，具体步骤为：

(a)构造一维序列M,N，其中M代表行，N代表列；

M＝[1,2,…,N^a]

(1)

N＝[1，2，…，N_r]

(b)计算按仿真中心时刻脉冲重复频率PRF_c进行均匀采样时的一维方位向各采样时刻序列T_uni(M)，并计算均匀采样时刻下的一维瞬时斜距序列R_uni(M)；

$T_{uni}\left(M\right)=\frac{M-N_a/2}{{\mathrm{PRF}}_c}---\left(2\right)$

(c)计算按PRF_cha采样时场景中心目标回波接收时刻与回波窗开启时刻的时延t_c(M)，并计算按PRF_c采样时场景中心目标回波接收时刻与回波窗开启时刻的时延t_u(M)；

其中，为向下取整函数；

(d)对二维回波仿真复数组S_start(M,N)做距离向傅里叶变换，即沿方位向进行快速傅里叶变换，得到二维复数组S_{loc_1}(M,N)；

S_{loc_1}(M,:)＝FFT(S_start(M,:))(6)

其中，S_start(M,:)表示S_start的第M行，S_{loc_1}(M,:)表示S_{loc_1}的第M行，FFT(·)表示对一维数组进行快速傅里叶变换；

(e)计算距离频域每列对应的距离频率数组f_τ(N)；

$f_{\mathrm{\&tau;}}\left(N\right)=\frac{N-N_r/2}{N_r}\&CenterDot;f_s---\left(7\right)$

(f)计算距离向位置补偿因子H₁(M,N)；

$H_1(M,N)=\exp \{-j\&CenterDot;2\mathrm{\&pi;}\&CenterDot;\&lsqb;(\frac{2R_{cha}\left(M\right)}{c}-\frac{2R_{uni}\left(M\right)}{c})-\left(t_c\right(M)-t_u(M\left)\right)\&rsqb;\&CenterDot;f_{\mathrm{\&tau;}}\left(N\right)\}---\left(8\right)$

其中j为单位复数；

(g)计算距离走动补偿因子H₂(M,N)；

H₂(M,N)＝exp{j·2π·f_τ(N)·[λfd₀T_cha(M)]}(9)

(h)将二维复数组S_{loc_1}(M,N)同距离向位置补偿因子H₁(M,N)与距离走动补偿因子H₂(M,N)相乘，得到二维复数组S_{loc_2}(M,N)；

S_{loc_2}(M,N)＝S_{loc_1}(M,N)·H₁(M,N)·H₂(M,N)(10)

(i)将二维复数组S_{loc_2}(M,N)做距离向傅里叶逆变换，即沿方位向进行快速傅里叶逆变换，得到调整距离向位置后的二维复数组S_loc(M,N)；

S_loc(M,:)＝IFFT(S_{loc_2}(M,:))(11)

其中，S_loc(M,:)表示S_start的第M行，S_{loc_2}(M,:)表示S_{loc_2}的第M行，IFFT(·)表示对一维数组进行快速傅里叶逆变换；

步骤四：方位向数据频谱搬移处理，具体可以分为以下几个步骤：

(a)构造方位向频谱搬移因子H₃(M)；

H₃(M)＝exp{j·2π·fd₀·T_cha(M)}(12)

(b)将二维复数组S_loc(M,N)同方位向频谱搬移因子H₃(M)相乘，得到方位向频谱搬移后的二维复数组S_base(M,N)；

S_base(M,N)＝S_loc(M,N)·H₃(M,N)(13)

步骤五：将步骤四得到的方位向频谱搬移后的二维复数组S_base做方位向拉格朗日插值处理，具体分为以下几个步骤：

(a)按列N计算均匀采样时刻序列T_uni(M)在不均匀采样时刻序列T_cha(M)中的位置p(M,N)；

具体为：以计算第一行第一列的位置p(1,1)为例，首先计算绝对差值Δt(n)＝|T_uni(1)-T_cha(n)|,n＝[1,2,…，N_a]，获取最小的绝对差值Δt_min和对应位置的n，若T_uni(1)＜T_cha(n)，p(1,1)＝n，若T_uni(1)≥T_cha(1)，p(1,1)＝n+1，以此类推，得到每一个位置p(M,N)；

(b)计算对二维均匀采样数组S_uni(M,N)进行拉格朗日插值所需的插值基函数L(M,N)；

$L(M,N)=\underset{\begin{array}{c}p=0\\ M\&NotEqual;p\end{array}}{\overset{N_l}{\&Pi;}}\frac{T_{uni}\left(M\right)-T_{cha}\&lsqb;p(M,N)-N_l/2+p\&rsqb;}{\begin{array}{c}{T}_{cha}\&lsqb;p(M,N)-N_l/2+M\&rsqb;-T_{cha}\&lsqb;p(M,N)-N_l/2+p\&rsqb;\\ N_l/2<pos\left(M,N\right)<N_a-N_l/2\end{array}}---\left(14\right)$

其中，N_l为拉格朗日插值阶数；

(c)结合拉格朗日插值基函数L(M,N)，分别对S_base(M,N)的实部S_{base_re}(M,N)和虚部S_{base_im}(M,N)分别进行拉格朗日插值法计算，得到方位向均匀采样的二维复数组S_uni(M,N)；

p(M,N)＜N_l/2,S_{uni_re}(M,N)＝S_{base_re}(M,N)

$N_l/2\&le;p(M,N)\&le;N_a-N_l/2,S_{uni\_re}(M,N)=\underset{k=0}{\overset{N_l}{\&Sigma;}}L\left(p\right(M,N)-N_l/2+k,N)\&CenterDot;S_{base\_re}\left(p\right(M,N)-N_l/2+k,N)---\left(15\right)$

p(M,N)＞N_a-N_l/2,S_{uni_re}(M,N)＝S_{base_re}(M,N_a)

p(M,N)＜N_l/2,S_{uni_im}(M,N)＝S_{base_im}(M,N)

$N_l/2\&le;p(M,N)\&le;N_a-N_l/2,S_{uni\_im}(M,N)=\underset{k=0}{\overset{N_l}{\&Sigma;}}L\left(p\right(M,N)-N_l/2+k,N)\&CenterDot;S_{base\_im}\left(p\right(M,N)-N_l/2+k,N)---\left(16\right)$

p(M,N)＞N_a-N_l/2,S_{uni_im}(M,N)＝S_{base_im}(M,N_a)

其中，S_{base_re}(M,N)是指二维复数组S_base第M行第N列数据的实部，S_{base_im}(M,N)是指二维复数组S_base第M行第N列数据的虚部，S_{uni_re}(M,N)是指二维数组S_uni第M行第N列数据的实部，S_{uni_im}(M,N)是指二维复数组S_uni第M行第N列数据的虚部；

步骤六：方位向数据频带恢复处理，具体可以分为以下几个步骤：

(a)构造方位向频带恢复因子H₄(M)；

H₄(M)＝exp{-j·2π·fd₀·T_uni(M)}(17)

(b)将二维复数组S_uni(M,N)同方位向频带恢复因子H₄(M)相乘，得到方位向频带恢复后的二维复数组S_new(M,N)；

S_new(M,N)＝S_base(M,N)·H₄(M)(18)

步骤七：将步骤六得到的复数组S_new做二维傅里叶变换，包括以下几个步骤：

(a)将步骤六得到的S_new(M,N)做距离向傅里叶变换，即沿方位向(按行)进行快速傅里叶变换(FFT)，得到方位时域距离频域复数组S_sp1(M,N)；

S_sp1(M,:)＝FFT(S_new(M,:))(19)

其中，S_sp1(M,:)表示S_sp1的第M行，S_new(M,:)表示S_new的第M行；

(b)将方位时域距离频域复数组S_sp1(M,N)做方位向傅里叶变换，即沿距离向(按列)进行快速傅里叶变换(FFT)，得到二维频域复数组S_sp(M,N)；

S_sp(:,N)＝FFT(S_sp1(:,N))(20)

其中，S_sp(:,N)表示S_sp的第N列，S_sp1(:,N)表示S_sp1的第N列；

步骤八：将步骤七得到的二维频域复数组S_sp(M,N)同一致压缩因子H₅(M,N)相乘，得到粗聚焦复数组S_bulk(M,N)，具体步骤如下：

(a)根据参考斜距R_ref计算最近斜距R_min；

$R_{\min }=R_{ref}-\frac{c}{2f_s}\&CenterDot;\frac{N_r}{2}---\left(21\right)$

(b)计算方位频域每行对应的方位频率一维数组f_a(M)；

$f_a\left(M\right)=\frac{M-N_a/2}{N_a}\&CenterDot;{\mathrm{PRF}}_c---\left(22\right)$

(c)计算二维一致压缩因子H₅(M,N)所需子因子ξ(M,N)；

(d)结合式(21)～式(23)计算大小为N_a×N_r的二维一致压缩因子H₅(M,N)；

(e)结合式(24)计算一致压缩后的粗聚焦数组S_bulk(M,N)；

S_bulk(M,N)＝S_sp(M,N)·H₅(M,N)(25)

步骤九：将步骤八得到的粗聚焦数组S_bulk(M,N)利用辛格插值法进行斯托尔特插值处理，得到精确聚焦二维复数组S_wave(M,N)，具体步骤如下：

(a)计算斯托尔特插值映射距离频率f′_τ(M,N)；

(b)遍历获取斯托尔特插值映射距离频率f′_τ(M,N)的最大值f′_τ,max与最小值f′_τ,min，并计算斯托尔特插值映射距离频率等分间隔Δf′_τ；

${\mathrm{\&Delta;f}}_{\mathrm{\&tau;}}^{\&prime;}=\frac{f_{\mathrm{\&tau;},max}^{\&prime;}-f_{\mathrm{\&tau;},min}^{\&prime;}}{N_r}---\left(27\right)$

(c)计算新距离频率均匀一维数组

$f_{\mathrm{\&tau;}}^u\left(N\right)=f_{\mathrm{\&tau;},min}^{\&prime;}+N\&CenterDot;{\mathrm{\&Delta;f}}_{\mathrm{\&tau;}}^{\&prime;}---\left(28\right)$

(d)计算二维频域复数组新距离频率均匀序列在每行对应的不均匀托尔特插值映射距离频率f′_τ(M,N)中的位置pos(M,N)。方法为按行进行下列操作：以计算第一行第一列的位置pos(1,1)为例，首先计算绝对差值获取最小的绝对差值Δk_min和对应位置的m，若pos(1,1)＝m-1，若pos(1,1)＝m，以此类推，得到每一个位置pos(M,N)；

(e)结合上步得到的位置pos(M,N)，计算辛格插值所需采样点位置q(M,N,n)；

$q(M,N,n)=\frac{f_{\mathrm{\&tau;}}^u\left(N\right)-f_{\mathrm{\&tau;}}^{\&prime;}(M,(pos\left(M,N\right)+n\left)\right)}{f_s/N_r},n=\&lsqb;-N_s/2,-N_s/2+1,...,N_s/2-1\&rsqb;---\left(29\right)$

其中，N_s是辛格插值核长度；

(f)结合式(29)利用辛格插值法，计算出经斯托尔特插值后的二维复数组S_wave(M,N)，由于S_wave(M,N)是复数组，需要分别对S_wave(M,N)的实部S_{wave_re}(M,N)和虚部S_{wave_im}(M,N)分别进行辛格插值法计算得出；

p(M,N)＜N_s/2,S_{wave_re}(M,N)＝S_{bulk_re}(M,N)

$N_s/2\&le;p(M,N)\&le;N_r-N_s/2,S_{wave\_re}(M,N)=\frac{\underset{n=-N_s/2}{\overset{N_s/2-1}{\&Sigma;}}{S}_{bulk\_re}(M,(pos\left(M,N\right)+n\left)\right)\&CenterDot;\sin c\left(q\right(M,N,n\left)\right)}{\underset{n=-N/2}{\overset{N/2-1}{\&Sigma;}}\sin c\left(q\right(M,N,n\left)\right)}---\left(30\right)$

p(M,N)＞N_r-N_s/2,S_{wave_re}(M,N)＝S_{bulk_re}(M,N_r)

p(M,N)＜N_s/2,S_{wave_im}(M,N)＝S_{bulk_im}(M,N)

$N_s/2\&le;p(M,N)\&le;N_r-N_s/2,S_{wave\_im}(M,N)=\frac{\underset{n=-N_s/2}{\overset{N_s/2-1}{\&Sigma;}}{S}_{bulk\_im}(M,(pos\left(M,N\right)+n\left)\right)\&CenterDot;\sin c\left(q\right(M,N,n\left)\right)}{\underset{n=-N/2}{\overset{N/2-1}{\&Sigma;}}\sin c\left(q\right(M,N,n\left)\right)}---\left(31\right)$

p(M,N)＞N_r-N_s/2,S_{wave_im}(M,N)＝S_{bulk_im}(M,N_r)

其中，sinc(·)是指插值函数S_{bulk_re}(M,N)是指二维数组S_bulk第M行第N列数据的实部，S_{bulk_im}(M,N)是指二维数组S_bulk第M行第N列数据的虚部，S_{wave_re}(M,N)是指二维数组S_wave第M行第N列数据的实部，S_{wave_im}(M,N)是指二维数组S_wave第M行第N列数据的虚部；

步骤十：将步骤九得到的复数组S_wave沿距离向(按列)进行快速傅里叶逆变换(IFFT)，得到方位时域距离频域复数组S_a-t；

S_a-t(:,N)＝IFFT(S_wave(:,N))(32)

其中，S_wave(:,N)表示S_wave的第N列，S_a-t(:,N)表示S_a-t的第N列；

步骤十一：将步骤十得到的方位时域距离频域复数组S_a-t同几何校正因子H₆相乘，得到经几何校正后的复数组S_geo，具体步骤如下：

(a)计算几何校正因子H₆(M)；

$H_6\left(M\right)=\exp \{-j2\mathrm{\&pi;}\&CenterDot;f_a\left(M\right)\&CenterDot;\frac{\mathrm{\&lambda;}\&CenterDot;{\mathrm{fd}}_0\&CenterDot;T_{uni}\left(M\right)}{c}\}---\left(33\right)$

(b)将方位时域距离频域复数组S_a-t同几何校正因子H₆(M)相乘；

S_geo(M,N)＝S_a-t(M,N)·H₆(M,N)(34)

步骤十二：将步骤十一得到的复数组S_geo沿方位向(按行)进行快速傅里叶逆变换(IFFT)，得到最终的成像结果S_end；

S_end(M,:)＝IFFT(S_geo(M,:))(35)

其中，S_geo(M,:)表示S_geo的第M行，S_end(M,:)表示S_end的第M行。