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

摘要

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

主权项

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]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)做距离向傅里叶变换，即沿方位向进行快速傅里叶变换，得到二维复数组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_{\tau }\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\pi ·[(\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_{\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_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,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)；$\begin{array}{c}L(M,N)=\underset{\underset{M\ne p}{p=0}}{\overset{N_l}{\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}{\Sigma }}L(p(M,N)-N_l/2+k,N)·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}{\Sigma }}L(p(M,N)-N_l/2+k,N)·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}·\frac{N_r}{2}---\left(21\right)$（b）计算方位频域每行对应的方位频率一维数组f_a(M)；$f_a\left(M\right)=\frac{M-N_a/2}{N_a}·{\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}}_{\tau }^{\prime }=\frac{f_{\tau ,\max }^{\prime }-f_{\tau ,\min }^{\prime }}{N_r}---\left(27\right)$（c）计算新距离频率均匀一维数组$f_{\tau }^u\left(N\right)=f_{\tau ,\min }^{\prime }+N·{\mathrm{\Delta f}}_{\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_{\tau }^u\left(N\right)-f_{\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}{\Sigma }}{S}_{\mathrm{bulk}\_\mathrm{re}}(M,(\mathrm{pos}(M,N)+n))·\sin c\left(q(M,N,n)\right)}{\underset{n=-N/2}{\overset{N/2-1}{\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}{\Sigma }}{S}_{\mathrm{bulk}\_\mathrm{im}}(M,(\mathrm{pos}(M,N)+n))·\sin c\left(q(M,N,n)\right)}{\underset{n=-N/2}{\overset{N/2-1}{\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\pi ·f_a\left(M\right)·\frac{\lambda ·{\mathrm{fd}}_0·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行。