CN103472449B - 一种基于mimo图像域的bp宽带合成方法 - Google Patents

Abstract

本发明公开了一种基于MIMO图像域BP宽带合成方法，其特点是首先利用MIMO体制下多个发射机同时发射中心频率步进的LFM信号，对接收到的各个窄带回波信号经BP成像得到各个低分辨率的子图像后，再对多个低分辨的子图像进行图像域的相参叠加以获得高分辨距离向的SAR图像。这种方法与现有的数据域的宽带合成方法相比，大大简化了宽带合成过程中的通道不一致性误差补偿步骤，在图像域进行宽带合成避免了数据域合成带来的子带重叠误差校正繁琐的问题，在降低运算量的同时实现了SAR距离向的高分辨率成像。本发明可适用于系统硬件条件(发射端瞬时发射带宽和接收机采样率)受限条件下的SAR距离向高分辨率成像的技术领域。

Description

一种基于MIMO图像域的BP宽带合成方法

技术领域：

本发明属于合成孔径雷达(Synthetic Aperture Radar，SAR)高分辨成像技术领域，它特别涉及到了系统硬件条件(发射端瞬时发射带宽和接收机采样率)受限条件下的SAR距离向高分辨率成像的技术领域。

背景技术：

合成孔径雷达是一种具有高分辨率的微波成像雷达，具有全天时和全天候工作的优点，并具有一定的植被和地面穿透能力，因而被广泛应用于军事及民用领域。在军事侦察、地质普查、地形测绘和制图、灾情预报等应用时，为了获取观测区域更多的信息，需对观测区域高分辨成像，实现高分辨率一直是SAR发展的重要目标。SAR利用载机的运动来模拟大孔径天线实现方位向高分辨；通过发射大带宽的线性调频信号(Linear Frequency Modulation，LFM)并利用脉冲压缩(Pulse Compression)技术获得距离向高分辨率。传统雷达系统由于受硬件复杂度和成本制约难以发射很大带宽的信号以获得距离向高分辨率，带宽合成技术就是利用在发射端发射多个中心频率步进的LFM信号，并在接收端将各个接收机接收的回波信号经过相参叠加形成等效的大带宽LFM信号，从而实现距离向高分辨率成像的技术。详见“邓云凯、陈倩、祁海明、郑慧芳、刘亚东.一种基于频域子带合成的多发多收高分辨率SAR成像算法[J].电子与信息学报.2011，33(5)：1082-1087”。

后向投影(Back Projection，BP)算法首先将合成孔径雷达原始数据沿距离向进行距离压缩(脉冲压缩)，然后通过选择不同慢时间观测空间中任意像素点在距离压缩后SAR数据空间中的数据，补偿方位向多普勒相位，并进行相干累加，最终获得各像素点散射系数的成像算法。由于BP算法可以有效的补偿运动误差，针对雷达平台抖动和非匀速直线运动带来的运动误差可以精确补偿以实现合成孔径雷达的精确成像，非常适用于宽带合成技术中对高分辨率的要求。详见“师君.双基地SAR与线阵SAR原理及成像技术研究[D].电子科技大学博士论文.2009”。

图像域宽带合成技术就是将各个接收机接收的中心频率步进的窄带LFM回波信号各自进行BP成像获得低分辨率的子图像并在图像域对各子带图像进行相参叠加以获得等效宽带信号的高分辨率距离向的SAR图像的技术。

多输入多输出(Multiple Input Multiple Output，MIMO)技术最早是由Marconi于1908年提出的用多天线来抑制信道衰落的。该技术可以解决传统宽带合成系统中单发单收SISO(Single-Input Single-Output)系统顺次发射步进频信号导致脉冲重复频率(Pulse RepetitionFrequency，PRF)过低引起测绘带受限的问题。

发明内容：

为了提高SAR图像距离向分辨率同时降低雷达系统对发射机发射信号带宽和接收机采样率的要求，本发明提出了一种基于MIMO图像域BP宽带合成方法，其特点是首先利用MIMO体制下多个发射机同时发射中心频率步进的LFM信号，对接收到的各个窄带回波信号经BP成像得到各个低分辨率的子图像后，再对多个低分辨的子图像进行图像域的相参叠加以获得高分辨距离向的SAR图像。这种方法与现有的数据域的宽带合成方法相比，大大简化了宽带合成过程中的通道不一致性误差补偿步骤，避免了数据域合成带来的子带重叠误差校正繁琐的问题，在降低运算量的同时实现了SAR距离向的高分辨率成像。

为了方便描述本发明的内容，首先作以下术语定义：

定义1、距离向

在雷达系统中，将雷达测距的范围划分成若干小的区域并将其编号，每个编号代表一个距离向。

定义2、方位向

将雷达扫描空间均匀划分为若干等分，每一等分为一个方位向。

定义3、快时间和慢时间

将距离向和方位向的时间分别定义为快时间和慢时间，“快”和“慢”是相对而言的。一般而言距离向的时间变化的快而方位向的时间变化的较慢。

定义4、斜距平面

本发明中的斜距平面是指由雷达速度方向矢量和斜距方向矢量构成的空间平面。

定义5、脉冲压缩

脉冲压缩是一种现代雷达信号处理技术，简单来说就是雷达发射宽脉冲，然后再接收端“压缩”为窄脉冲，从而改善雷达的两种性能：作用距离和距离分辨率。详见“皮亦鸣，杨建宇，付毓生，杨晓波.合成孔径雷达成像原理.第一版.电子科技大学出版社.2007.3”。

定义6、升采样

本发明中的升采样(Up sampling)技术是指将信号变换到频域后在原始信号的每两个采样点之间插L-1个零，然后再变换到时域的过程，其中L为升采样倍数。

定义7、脉冲重复频率

本发明中的脉冲重复频率(Pulse Repetition Frequency，PRF)是指雷达1秒钟内发射脉冲的个数，即脉冲重复的频率，在这里也即方位向的采样率。

定义8、步进频率

本发明中的步进频率(Step Frequency)是指发射脉冲信号的相邻窄带信号中心频率的间隔。

定义9、BP算法

后向投影(Back Projection，BP)算法是一种首先将合成孔径雷达原始数据沿距离向进行距离压缩，然后从距离压缩后的数据中选择空间中各散射点回波，并进行相干累加以恢复各散射点的散射系数的成像算法。详见“师君.双基地SAR与线阵SAR原理及成像技术研究[D].电子科技大学博士论文.2009”。

定义10、快速傅里叶变换

计算离散傅里叶变换的一种快速算法，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。详见“程乾生.数字信号处理.北京大学出版社，北京，2003”。

定义11、线性调频信号

线性调频(Linear Frequency Modulation，LFM)信号又称为Chirp信号，它是研究最早而且应用最广泛的一种脉冲压缩信号。采用线性调频脉冲压缩技术的雷达可以同时获得远的作用距离和高的距离分辨力。线性调频信号的主要优点是所用的匹配滤波器对回波的多普勒频移不敏感，即使回波信号有较大的多普勒频移，仍能用同一个匹配滤波器完成脉冲压缩，这将大大简化信号处理系统。详见“皮亦鸣，杨建宇，付毓生，杨晓波.合成孔径雷达成像原理.第一版.电子科技大学出版社.2007.3”。

定义12、MIMO系统

多输入多输出(Multiple Input Multiple Output，MIMO)技术最早是由Marconi于1908年提出用多天线来抑制信道衰落的。在本发明中，MIMO系统就是在发射端有多个发射天线(通道)用于发射雷达信号，同样的在接收端有多个接收天线(通道)用于接收雷达回波的系统。

本发明提供了一种基于MIMO图像域的BP宽带合成的方法，该方法的步骤如下：

步骤1、用于基于MIMO图像域的BP宽带合成方法相关参数的初始化

初始化的参数均为已知，且初始化的参数如下：所有的坐标都是以笛卡尔直角坐标形式给出；光速为C；雷达发射线性调频信号，脉冲重复频率为PRF，脉冲重复周期为PRI，MIMO阵列天线个数为N，也即每个脉冲重复周期内雷达发射N个中心频率步进的窄带(子带)线性调频信号，这N个天线同时也应用于接收各自信号反射的回波；雷达发射的各个子带信号的带宽为B_pn；雷达发射的各个子带信号的时宽为T_pn；发射信号的调频斜率记为γ，发射脉冲各子带信号的步进频率为Δf；经宽带合成后的宽带信号的中心频率为f₀；各子带信号的中心频率记为f_c(n)，f_c(n)＝f₀+nΔf，n＝0，1，…，N-1；发射平台中心Pt₀；场景中心P_center；雷达参考斜距R₀；雷达在距离向上窄带信号的采样频率F_pn；雷达载机平台速度为V；合成孔径长度为L_sar；合成孔径时间为T_sar，一个合成孔径时间内距离向采样点数为N_r(便于FFT计算，在此取值为2的整数次幂)，方位向采样点数为N_a(便于FFT计算，在此取值为2的整数次幂)；快时间为一个长度为N_r的行矢量， $t=\frac{[-\frac{N_r}{2},-\frac{N_r}{2}+1,...,\frac{N_r}{2}-2,\frac{N_r}{2}-1]}{F_{\mathrm{pn}}};$ 慢时间ts为一个长度为N_a的行矢量， $\mathrm{ts}=\frac{[-\frac{N_a}{2},-\frac{N_a}{2}+1,...,\frac{N_a}{2}-2,\frac{N_a}{2}-1]}{\mathrm{PRF}};$ BP成像场景距离向(斜距平面)和方位向点数分别为M和K；BP成像场景距离向(斜距平面)和方位向网格间距为dx和dy；目标位置(x₀，y₀)。

步骤2、获取雷达回波数据

第n个子带信号回波数据为一个N_a×N_r的矩阵，记为S⁽ⁿ⁾，其表现形式如下：

其中第n个子带信号的第jj个快时间、第ii个慢时间的雷达回波数据记为 $S_{\mathrm{ii},\mathrm{jj}}^{\left(n\right)}=\mathrm{rect}\left(\frac{(t\left(\mathrm{jj}\right)-\mathrm{\τ}\left(\mathrm{ii}\right))}{T_{\mathrm{pn}}}\right)\·\exp \left(\mathrm{j\π\γ}{(t\left(\mathrm{jj}\right)-\mathrm{\τ}\left(\mathrm{ii}\right))}^2\right)\·\exp (-j2\mathrm{\π}{f}_c\left(n\right)\·\mathrm{\τ}\left(\mathrm{ii}\right)),$ (n＝1，…，N；ii＝1，2，…，N_a；jj＝1，2，…，N_r)，N、N_a、N_r分别为子带数、方位采样点数、距离向采样点数，其中τ(ii)为第ii个方位时刻雷达天线到场景目标点的双程距离延迟时间，ts(ii)为ts的第ii个元素，exp(·)为以自然底数e为底的指数函数，j为虚数单位(下同)，rect(t/T_pn)为(t、T_pn分别为距离向快时间、子带时宽)矩形窗函数，满足

步骤3、对雷达回波数据矩阵每一行(距离向)进行脉冲压缩

取出步骤2中所获得的雷达回波数据构造距离向参考函数H_r，利用传统的脉冲压缩方法对S⁽ⁿ⁾的每一行(距离向)进行脉冲压缩，脉冲压缩后的矩阵记为PS⁽ⁿ⁾，PS⁽ⁿ⁾＝IFFT(FFT((S⁽ⁿ⁾)′)·FFT(H_r))，其中(S⁽ⁿ⁾)′表示S⁽ⁿ⁾的转置，PS⁽ⁿ⁾的第ii行第jj列记为其中n＝1，…，N；ii＝1，2，…，N_a；jj＝1，2，…，N_r。其中H_r为一个N_a×N_r的矩阵，其表现形式如下：

H_r＝ones(N_a，1)×rect(t/T_pn)·exp(jπγt²)

其中ones(N_a，1)产生N_a×1的全1矩阵，FFT(·)、IFFT(·)分别为快速傅里叶变换和快速傅里叶逆变换。

步骤4、对每个距离向的所有回波数据矩阵进行升采样L倍

对步骤3中脉冲压缩后的N个子带、N_a个方位向数据矩阵统一做如下处理：

取出步骤3中处理后的第n个子带信号的第ii个方位向的脉压后的数据n＝1，…，N；ii＝1，2，…，N_a，表示PS⁽ⁿ⁾矩阵的第ii行(第ii个方位向)数据。

首先将第n个子带信号的第ii个方位向的脉压后的数据变换到频域得到数据记为 然后将数据中间插零，得到数据矩阵记为 ${\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},:}^{\left(n\right)}}}^{\′}=[\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},1}^{\left(n\right)}},\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},2}^{\left(n\right)}},\·\·\·,\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},\frac{N_r}{2}}^{\left(n\right)}},\mathrm{SO},\stackrel{\‾}{P{S}_{\mathrm{ii},\frac{N_r}{2}+1}^{\left(n\right)}},\stackrel{\‾}{P{S}_{\mathrm{ii},\frac{N_r}{2}+2}^{\left(n\right)}},\·\·\·,\stackrel{\‾}{P{S}_{\mathrm{ii},N_r}^{\left(n\right)}}],$ 其中n=1，…，N；ii=1，2，…，N_a，其中S0为1行(L-1)·N_r的零矩阵；

然后将数据矩阵记为变换到时域 ${\mathrm{RS}}_{\mathrm{ii},:}^{\left(n\right)}=\mathrm{IFFT}\left({\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},:}^{\left(n\right)}}}^{\′}\right),(n=1,...,N;\mathrm{ii}=\mathrm{1,2},...,N_a).$

步骤5、对每个方位向数据矩阵进行相关叠加进行BP成像

对步骤4中升采样后的N个子带、N_a个方位向数据矩阵统一做如下处理：

步骤5.1取出步骤4中第n个子带、第ii个升采样后的方位向数据矩阵（其中n=1，…，N；ii=1，2，…，N_a)，对第n个子带、第ii个升采样后的方位向数据矩阵做距离向FFT变换到距离一多普勒域，即其中是对的转置操作；然后找到第ii个方位时刻雷达平台的位置P_t=V·ts与成像场景中的点T(mm,kk)(mm=1,2，…，M；kk=1,2，…，K，其中T(mm,kk)=[(mm-M/2)·dx，(kk-K/2)·dy])的距离历史为记为R(ii)， $R\left(\mathrm{ii}\right)=\sqrt{{(P_t(2,\mathrm{ii})-T\left(2\right))}^2+{(T\left(1\right)+R_r)}^2}$ (T(1)表示场景中第(mm,kk）个点的横坐标，T(2)表示场景中第(mm，kk)个点的纵坐标)，其中R_r为斜距平面内的场景参考斜距，并计算第ii个方位向用于BP成像搜索的距离向编号，记为ID， $\mathrm{ID}=\mathrm{round}((\frac{2R\left(\mathrm{ii}\right)-{2R}_0}{c}*F_{\mathrm{pn}}+\frac{N_r}{2})*L+0.5),$ 其中round(·)为取四舍五入的函数；最后对第ii个方位向、ID对应的距离向数据补偿方位向多普勒相位得到的图像矩阵记为I⁽ⁿ⁾，I⁽ⁿ⁾的第mm行第kk列记为 (其中λ(n)=C/f(n)为第n个子带信号的波长)，

步骤5.2按照步骤5.1对场景每个点(mm，kk)(mm=1,2，…，M；kk=1,2，…，K，M和K分别为BP成像场景距离向(斜距平面)和方位向点数)，补偿N_a个方位向的数据(ii=1,2，…，N_a；n=1，…，N)的多普勒相位，并做相参叠加就可得到低分辨距离向的子带图I⁽ⁿ⁾，n=1，…，N，I⁽ⁿ⁾为一个M×K的矩阵。

步骤6、对每个距离向低分辨的子图像(图像域)作相位补偿

对步骤5中处理后的N个距离向低分辨的子图像矩阵统一做如下处理：

取出步骤5中第n个距离向低分辨的子图像矩阵I⁽ⁿ⁾，n=1，…，N，对第n个距离向低分辨的子图像矩阵I⁽ⁿ⁾做相位补偿得到补偿后的子图像记为 其中Φ为相位补偿函数，其表达式为Φ=exp(-j4πf_c(n)·r/C)，其中r=C·t/2为斜距平面对应的距离轴。

步骤7、对每个子图像(图像域)作频移操作

构造第n个子图像对应的频移参考函数S_ref(n)=exp(j4π(n-(N+1)/2)·Δf·r/C)，其中n=1，…，N，将该参考函数与子图像矩阵相乘进行频移操作得到频移后的数据记为 $\stackrel{\‾}{I_c^{\left(n\right)}}=I_c^{\left(n\right)}\·S_{\mathrm{ref}}\left(n\right),$ 其中n=1，…，N。

步骤8、对频移后的子图像作FFT操作变换到波数域

对步骤7中处理后的N个子图像数据矩阵统一做如下处理：

取出步骤7中第n个子数据矩阵n=1，…，N，利用快速傅里叶变换方法对的每一行进行快速傅里叶变换，得到处理后的M行K列的矩阵记为 ${\stackrel{\‾}{I_c^{\left(n\right)}}}^{\′}=\mathrm{FFT}\left({\left(\stackrel{\‾}{I_c^{\left(n\right)}}\right)}^{\′}\right),n=1,...,N,$ 其中是对的转置矩阵。

步骤9、对波数域各子带图像进行相参叠加

取出步骤8中得到的N个子图像矩阵进行相干叠加得到波数域的等效宽带信号记为I_full， $I_{\mathrm{full}}={\mathrm{\Σ}}_{n=1}^N{\stackrel{\‾}{I_c^{\left(n\right)}}}^{\′},n=1,...,N.$

步骤10、将合成的波数域的宽带图像变换到空域

对步骤9中合成的波数域的M个距离向行向量统一做如下处理：

取出步骤9中第mm个距离向的行向量I_full(mm，：)，其中mm=1，2，…，M，I_full(mm，：)表示I_full的第mm行的所有数据，将其进行快速傅里叶逆变换到空域即可得到高分辨距离向的图像记为Image，Image(mm,：)=IFFT(I_full(mm，：))，mm=1，2，…，M，则Image即为最终高分辨距离向图像矩阵。

本发明的创新点在于将宽带合成技术应用到图像域，本发明提供一种基于MIMO图像域BP宽带合成的方法。该方法首先利用传统BP成像算法对MIMO阵列天线接收的各子带回波信号进行成像得到距离向低分辨率的子带图像，然后在图像域对各个低分辨率的子带图像进行频移叠加得到高分辨率距离向的图像。这种方法大大简化了宽带合成过程中相位补偿操作。

本发明的优点：本发明利用MIMO阵列天线在获得宽测绘带的同时不会造成方位模糊；应用BP算法进行子带成像大大简化了运动误差补偿步骤；在图像域进行宽带合成避免了数据域合成带来的子带重叠误差校正繁琐的问题。

附图说明

图1为MIM0阵列成像几何关系图

直角坐标系中X，Y，Z轴分别代表距离向、方位向、高度向；P(X₀，Y₀)为地面散射点(目标)对应的位置；PRI为脉冲重复周期；ii代表第ii个慢时间时刻；V为雷达平台的运动速度；f₁、f₂、f₃分别表示3个天线发射信号的中心频率。

图2为传统BP算法的流程图

PRF为脉冲重复频率；APC(Antenna Phase Center,APC)轨迹即为天线相位中心轨迹；(u，v)表示成像场景对应坐标为(u，v)的点；为第i个方位时刻天线相位中心位置与成像场景中的点(u，v)的距离，为场景中的点(u，v)对应的位置。

图3为系统流程图

图中T_x/R_x表示收/发天线；B_p1、B_p2、B_p3分别表示3个子带信号的带宽；Irmge1、Image2、Image3表示3幅对应的子图像。

具体实施方式

本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2010b上验证正确。具体实施步骤如下：

步骤1、用于基于MIMO图像域的BP宽带合成方法相关参数的初始化

初始化的参数均为已知，且初始化的参数如下：所有的坐标都是以笛卡尔直角坐标形式给出；光速为C=3×10⁸m/s；雷达发射线性调频信号，脉冲重复频率记为PRF,PRF=500Hz，脉冲重复周期记为PRI，MIMO天线个数为N=3，也即每个脉冲重复周期内雷达发射3个中心频率步进的窄带(子带)线性调频信号，这3个天线同时也应用于接收各自信号反射的回波；雷达发射的各个子带信号的带宽为B_pn=300MHz；雷达发射的各个子带信号的时宽为T_pn=1μs；发射信号的调频斜率记为γ，发射脉冲各子带信号的步进频率为Δf=300MHz；经宽带合成后的宽带信号的中心频率为f₀=10GHz；各子带信号的中心频率记为f_c(n)，f_c(n)=f₀+nΔf=(10+n×0.3)×10⁹，n=0，1,2；发射平台中心Pt₀=[0,0,4000]；场景中心P_cener=[3000,0,0]；雷达参考斜距R₀=5000；雷达在距离向上窄带信号的的采样频率F_pn=390MHz；雷达载机平台速度为V=[0，100,0]；合成孔径长度为L_sar=400m；合成孔径时间记为 一个合成孔径时间内距离向采样点数为N_r=512，方位向采样点数为N_a=512；快时间为慢时间为BP成像场景距离向(斜距平面)和方位向点数分别为M=30和K=30；BP成像场景距离向(斜距平面)和方位向网格间距为dx=0.5m和dy=0.5m；目标位置(5000，0)。

步骤2、获取雷达回波数据

第n个子带信号回波数据为一个512×512的矩阵，记为S⁽ⁿ⁾，其表现形式如下：

其中第n个子带信号的第jj个快时间、第ii个慢时间的雷达回波数据记为， $\begin{array}{c}{s}_{\mathrm{ii},\mathrm{jj}}^{\left(n\right)}=\mathrm{rect}((\frac{\mathrm{jj}}{3.9}\×{10}^{-8}-\frac{\mathrm{ii}}{500})\×{10}^6)\·\exp (\mathrm{j\π}3\×{10}^2\×{(\frac{\mathrm{jj}}{3.9}\×{10}^{-8}-\frac{\mathrm{ii}}{500})}^2)\×\\ \exp (-j2\mathrm{\π}(10+n\×0.3)\×{10}^9\×\mathrm{\τ}\left(\mathrm{ii}\right)),(n=\mathrm{1,2,3};\mathrm{ii}=\mathrm{1,2},\·\·\·,512;\mathrm{jj}=\\ \mathrm{1,2},\·\·\·,512)\end{array},,$ 其中τ(ii)为第ii个方位时刻雷达天线到场景目标点的双程距离延迟时间， $\mathrm{\τ}\left(\mathrm{ii}\right)=2\sqrt{{5000}^2+10000\×{\left(\frac{\mathrm{ii}}{500}\right)}^2}/3{\×10}^8),$ exp(·)为以自然底数e为底的指数函数，j为虚数单位(下同)，为平方根操作(下同)，rect(t/T_pn)(其中t、T_pn分别为距离向快时间、子带时宽)为矩形窗函数，矩形窗函数，满足

步骤3、对雷达回波数据矩阵每一行(距离向)进行脉冲压缩

取出步骤2所得到的所有雷达回波数据构造距离向参考函数H_r，利用脉冲压缩方法对S⁽ⁿ⁾的每一行(距离向)进行脉冲压缩得到脉冲压缩后的数据记为PS⁽ⁿ⁾，PS⁽ⁿ⁾＝IFFT(FFT((S⁽ⁿ⁾)′)·FFT(H_r))，其中(S⁽ⁿ⁾)′表示S⁽ⁿ⁾的转置，PS⁽ⁿ⁾的第ii行第jj列的数据记为其中n＝1，2，3；ii＝1，2，…，512；jj＝1，2，…，512。其中H_r为一个512×512的矩阵，其表现形式如下：

$\begin{array}{c}{H}_r=\mathrm{ones}\left(\mathrm{512,1}\right)\×\mathrm{rect}\left(\frac{[-256,-255,\·\·\·\mathrm{254,255}]}{390}\right)\×\exp (\mathrm{j\π}\×3\\ \×{10}^{-2}\×{\left(\frac{[-256,-255,\·\·\·\mathrm{254,255}]}{3.9}\right)}^2)\end{array}$

其中ones(512，1)产生512×1的全1矩阵，FFT(·)、IFFT(·)分别为快速傅里叶变换和快速傅里叶逆变换。

步骤4、对每个距离向的所有回波数据矩阵进行升采样8倍

对步骤3中脉冲压缩后的N个子带、N_a个方位向数据矩阵统一做如下处理：

取出步骤3中处理后的第n个子带信号的第ii个方位向的脉压后的数据n＝1，2，3；ii＝1，2，…，512，表示PS⁽ⁿ⁾矩阵的第ii行(第ii个方位向)数据(下同)。首先将这行数据变换到频域得到的数据记为 然后将该行数据中间插零后得到的数据记为 ${\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},1}^{\left(n\right)}}}^{\′}=[\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},1}^{\left(n\right)}},\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},2}^{\left(n\right)}},\·\·\·\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},256}^{\left(n\right)}},S0,\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},257}^{\left(n\right)}},\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},258}^{\left(n\right)}},\·\·\·\·,\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},511}^{\left(n\right)}},\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},512}^{\left(n\right)}}],$ 其中n＝1，2，3；ii＝1，2，…，512，其中S0为1行3584列的零矩阵；最后将其变换到时域得到的矩阵记为 ${\mathrm{RS}}_{\mathrm{ii},1}^{\left(n\right)}=\mathrm{IFFT}\left({\stackrel{\‾}{{\mathrm{PS}}_{\mathrm{ii},1}^{\left(n\right)}}}^{\′}\right),(n=\mathrm{1,2,3};\mathrm{ii}=\mathrm{1,2},...,512).$

步骤5、对每个方位向数据矩阵进行相关叠加进行BP成像

对步骤4中升采样后的3个子带、512个方位向数据矩阵统一做如下处理：

首先取出步骤4中第n个子带、第ii个升采样后的方位向数据矩阵(其中n＝1，2，3；ii＝1，2，…，512)对其做距离向FFT变换到距离-多普勒域，得到的数据记为 其中是对的转置操作；然后找到第ii个方位时刻雷达平台的位置 $P_t=\left[\mathrm{0,100,0}\right]\×\frac{[-256,-255,...\mathrm{254,255}]}{500}$ 与成像场景中的点T(mm，kk)(mm＝1，2，…，30；kk＝1，2，…，30，其中T(mm，kk)＝[(mm-15)×0.5，(kk-15)×0.5])                             的距离历史记为R(ii)， $R\left(\mathrm{ii}\right)=\sqrt{{(100\×\frac{\mathrm{ii}}{500}-(\mathrm{kk}-15)\×0.5)}^2+{((\mathrm{mm}-15)\×0.5+5000)}^2}$ 并计算第ii个方位向用于BP成像搜索距离向编号记为ID，

$\begin{array}{c}\mathrm{ID}=\\ \mathrm{round}(((\sqrt{{(100\×\frac{\mathrm{ii}}{500}-(\mathrm{kk}-15)\×0.5)}^2+{((\mathrm{mm}-15)\×0.5+5000)}^2}-\\ 5000)\×2.6+256)\×8+0.5)\end{array},$ 其中round(·)为取四舍五入的函数；最后对第ii个方位向、ID对应的距离向数据补偿方位向多普勒相位后得到的矩阵记为I⁽ⁿ⁾，I⁽ⁿ⁾的第mm行第kk列记为则

$\begin{array}{c}{I}_{\mathrm{mm},\mathrm{kk}}^{\left(n\right)}=\stackrel{\‾}{{\mathrm{RS}}_{11,\mathrm{ID}}^{\left(n\right)}}\×\\ \exp (j4\mathrm{\π}\sqrt{{(100\×\frac{\mathrm{ii}}{500}-(\mathrm{kk}-15)\×0.5)}^2+{((\mathrm{mm}-15)\×0.5+5000)}^2}/0.3/\\ (10+n\×0.3))\end{array},$ 按上述步骤对场景中任一点(mm，kk)(mm＝1，2，…，30；kk＝1，2，…，30)补偿512个方位向的数据的多普勒相位，并做相参叠加就可得到低分辨距离向的子带图I⁽ⁿ⁾，n＝1，2，3，I⁽ⁿ⁾为一个M×K的矩阵。

步骤6、对每个距离向低分辨的子图像(图像域)作相位补偿

对步骤5中处理后的3个距离向低分辨的子图像矩阵统一做如下处理：

取出步骤5中第n个距离向低分辨的子图像矩阵I⁽ⁿ⁾，n＝1，2，3，对其做相位补偿得到补偿后的子图像得到的数据记为 其中Φ为相位补偿函数，其表达式为 $\mathrm{\Φ}=\exp (-j4\mathrm{\π}\×(10+n\×0.3)\×10\×\frac{[-256,-255,...\mathrm{254,255}]}{7.8}),n=\mathrm{1,2,3}.$

步骤7、对每个子图像(图像域)作频移操作

构造第n个子图像对应的频移参考函数 $S_{\mathrm{ref}}\left(n\right)=\exp (j4\mathrm{\π}(n-2)\×\frac{[-256,-255,...\mathrm{254,255}]}{2.6}),$ 其中n＝1，2，3，将该参考函数与子图像矩阵相乘进行频移操作得到的数据记为 $\stackrel{\‾}{I_c^{\left(n\right)}}=I_c^{\left(n\right)}\·S_{\mathrm{ref}}\left(n\right).$

步骤8、对频移后的子图像作FFT操作变换到波数域

对步骤7中处理后的3个子图像数据矩阵统一做如下处理：

取出步骤7中第n个子数据矩阵利用快速傅里叶变换方法对的每一行进行快速傅里叶变换，得到处理后的M行K列的矩阵记为 其中是对的转置矩阵。

步骤9、对波数域各子带图像进行相参叠加

取出步骤8中得到的3个子图像矩阵进行相干叠加得到波数域的等效宽带信号数据记为I_full， $I_{\mathrm{full}}={\mathrm{\Σ}}_{n=1}^N{\stackrel{\‾}{I_c^{\left(n\right)}}}^{\′},n=\mathrm{1,2,3}.$

步骤10、将合成的波数域的宽带图像变换到空域

对步骤9中合成的波数域的30个距离向行向量统一做如下处理：

取出步骤9中第mm个距离向的行向量I_full(mm，：)，mm＝1，2，…，30，其中I_full(mm，：)为I_full的第mm行，将其进行快速傅里叶逆变换到空域即可得到高分辨距离向的图像数据记为Image，Image(mm，：)＝IFFT(I_full(mm，：))，Image(mm，：)为Image的第mm行，mm＝1，2，…，30，则Image即为最终高分辨距离向图像矩阵。

通过本发明的具体实施可以看出，本发明通过利用MIMO阵列天线各自同时发射中心频率步进的LFM信号并同时接收的模式，解决了用于传统宽带合成SISO系统顺序发射步进频信号带来的不可调和的方位向高分辨率和大测绘带的矛盾；运用BP对各子带信号进行成像，大大简化了由于平台抖动和非直线运动带来的运动误差补偿步骤；在图像域进行宽带合成避免了数据域合成带来的子带重叠误差校正繁琐的问题。另外由于对每个通道的数据都做了相同的处理，可以结合GPU并行化处理，提高系统的实时性。

Claims (1)

1.一种基于MIMO图像域的BP宽带合成的方法，其特征是它包括以下步骤：

步骤1、用于基于MIMO图像域的BP宽带合成方法相关参数的初始化

初始化的参数均为已知，且初始化的参数如下：所有的坐标都是以笛卡尔直角坐标形式给出；光速为C；雷达发射线性调频信号，脉冲重复频率为PRF，脉冲重复周期为PRI，MIMO阵列天线个数为N，也即每个脉冲重复周期内雷达发射N个中心频率步进的窄带即子带线性调频信号，这N个天线同时也应用于接收各自信号反射的回波；雷达发射的各个子带信号的带宽为B_pn；雷达发射的各个子带信号的时宽为T_pn；发射信号的调频斜率记为γ，发射脉冲各子带信号的步进频率为Δf；经宽带合成后的宽带信号的中心频率为f₀；各子带信号的中心频率记为f_c(n)，f_c(n)＝f₀+nΔf，n＝0，1，…，N-1；发射平台中心Pt₀；场景中心P_center；雷达参考斜距R₀；雷达在距离向上窄带信号的采样频率F_pn；雷达载机平台速度为V；合成孔径长度为L_sar；合成孔径时间为T_sar，一个合成孔径时间内距离向采样点数为N_r，便于FFT计算，在此取值为2的整数次幂，方位向采样点数为N_a，便于FFT计算，采样点数为取值为2的整数次幂；快时间为一个长度为N_r的行矢量， $t=\frac{[-\frac{N_r}{2},-\frac{N_r}{2}+1,\·\·\·,\frac{N_r}{2}-2,\frac{N_r}{2}-1]}{F_{\mathrm{pn}}};$ 慢时间ts为一个长度为N_a的行矢量， $\mathrm{ts}=\frac{[-\frac{N_a}{2},-\frac{N_a}{2}+1,\·\·\·,\frac{N_a}{2}-2,\frac{N_a}{2}-1]}{\mathrm{PRF}};$ BP成像场景距离向即斜距平面和方位向点数分别为M和K；BP成像场景距离向即斜距平面和方位向网格间距为dx和dy；目标位置(x₀，y₀)；

步骤2、获取雷达回波数据

第n个子带信号回波数据为一个N_a×N_r的矩阵，记为S⁽ⁿ⁾,其表现形式如下：

其中第n个子带信号的第jj个快时间、第ii个慢时间的雷达回波数据记为 $S_{\mathrm{ii},\mathrm{jj}}^{\left(n\right)}=\mathrm{rect}\left(\frac{(t\left(\mathrm{jj}\right)-\mathrm{\τ}\left(\mathrm{ii}\right))}{T_{\mathrm{pn}}}\right)\·\exp \left(\mathrm{j\π\γ}{(t\left(\mathrm{jj}\right)-\mathrm{\τ}\left(\mathrm{ii}\right))}^2\right)\·\exp (-j2\mathrm{\π}{f}_c\left(n\right)\·\mathrm{\τ}\left(\mathrm{ii}\right)).$

n＝1，…，N；ii＝1，2，…，N_a；jj＝1，2，…，N_r，N、N_a、N_r分别为子带数、方位采样点数、距离向采样点数，其中τ(ii)为第ii个方位时刻雷达天线到场景目标点的双程距离延迟时间，ts(ii)为ts的第ii个元素，exp(·)为以自然底数e为底的指数函数，j为虚数单位，rect(t/T_pn)为矩形窗函数，满足

t、T_pn分别为距离向快时间、子带时宽；

步骤3、对雷达回波数据矩阵每一行即距离向进行脉冲压缩

取出步骤2中所获得的雷达回波数据构造距离向参考函数H_r，利用传统的脉冲压缩方法对S⁽ⁿ⁾的每一行即距离向进行脉冲压缩，脉冲压缩后的矩阵记为PS⁽ⁿ⁾，PS⁽ⁿ⁾＝IFFT(FFT((S⁽ⁿ⁾)′)·FFT(H_r))，其中(S⁽ⁿ⁾)′表示S⁽ⁿ⁾的转置，PS⁽ⁿ⁾的第ii行第jj列记为其中n＝1，…，N；ii＝1，2，…，N_a；jj＝1，2，…，N_r，其中H_r为一个N_a×N_r的矩阵，其表现形式如下:

H_r＝ones(N_a，1)×rect(t/T_pn)·exp(jπγt²)

其中ones(N_a，1)产生N_a×1的全1矩阵，FFT(·)、IFFT(·)分别为快速傅里叶变换和快速傅里叶逆变换；

步骤4、对每个距离向的所有回波数据矩阵进行升采样L倍

对步骤3中脉冲压缩后的N个子带、N_a个方位向数据矩阵统一做如下处理：

取出步骤3中处理后的第n个子带信号的第ii个方位向的脉压后的数据n＝1，…，N；ii＝1，2，…，N_a，表示PS⁽ⁿ⁾矩阵的第ii行即第ii个方位向数据；首先将第n个子带信号的第ii个方位向的脉压后的数据变换到频域得到数据记为 然后将数据中间插零，得到数据矩阵记为 其中

n＝1，…，N；ii＝1，2，…，N_a，其中S0为1行(L-1)·N_r的零矩阵；然后将数据矩阵记为变换到时域n＝1，…，N；ii＝1，2，…，N_a；

步骤5、对每个方位向数据矩阵进行相关叠加进行BP成像

对步骤4中升采样后的N个子带、N_a个方位向数据矩阵统一做如下处理：

步骤5.1取出步骤4中第n个子带、第ii个升采样后的方位向数据矩阵其中n＝1，…，N；ii＝1，2，…，N_a,对第n个子带、第ii个升采样后的方位向数据矩阵做距离向FFT变换到距离-多普勒域，即其中是对的转置操作；然后找到第ii个方位时刻雷达平台的位置P_t＝V·ts与成像场景中的点T(mm，kk)＝[(mm-M/2)·dx，(kk-K/2)·dy]，mm＝1，2，…，M，kk＝1，2，…，K距离历史为记为R(ii)，T(1)表示场景中第(mm，kk)个点的横坐标，T(2)表示场景中第(mm，kk)个点的纵坐标,其中R_r为斜距平面内的场景参考斜距，并计算第ii个方位向用于BP成像搜索的距离向编号，记为ID，其中round(·)为取四舍五入的函数；最后对第ii个方位向、ID对应的距离向数据补偿方位向多普勒相位得到的图像矩阵记为I⁽ⁿ⁾，I⁽ⁿ⁾的第mm行第kk列记为 $I_{\mathrm{mm},\mathrm{kk}}^{\left(n\right)}=\stackrel{\‾}{{\mathrm{RS}}_{\mathrm{ii},\mathrm{ID}}^{\left(n\right)}}\·\exp (j4\mathrm{\πR}\left(\mathrm{ii}\right)/\mathrm{\λ}\left(n\right)),$ 其中λ(n)＝C/f(n)为第n个子带信号的波长，

步骤5.2按照步骤5.1对场景每个点(mm，kk)，mm＝1，2，…，M；kk＝1，2，…，K，M和K分别为BP成像场景距离向即斜距平面和方位向点数，补偿N_a个方位向的数据ii＝1，2，…，N_a；n＝1，…，N的多普勒相位，并做相参叠加就可得到低分辨距离向的子带图I⁽ⁿ⁾，n＝1，…，N，I⁽ⁿ⁾为一个M×K的矩阵；

步骤6、对每个距离向低分辨的子图像即图像域作相位补偿

对步骤5中处理后的N个距离向低分辨的子图像矩阵统一做如下处理：

取出步骤5中第n个距离向低分辨的子图像矩阵I⁽ⁿ⁾，n＝1，…，N，对第n个距离向低分辨的子图像矩阵I⁽ⁿ⁾做相位补偿得到补偿后的子图像记为 其中Φ为相位补偿函数，其表达式为Φ＝exp(-j4πf_c(n)·r/C),其中r＝C·t/2为斜距平面对应的距离轴；

步骤7、对每个子图像即图像域作频移操作

构造第n个子图像对应的频移参考函数S_ref(n)＝exp(j4π(n-(N+1)/2)·Δf·r/C),其中n＝1，…，N，将该参考函数与子图像矩阵相乘进行频移操作得到频移后的数据记为 $\stackrel{\‾}{I_c^{\left(n\right)}}=I_c^{\left(n\right)}\·S_{\mathrm{ref}}\left(n\right),$ 其中n＝1，…，N；

步骤8、对频移后的子图像作FFT操作变换到波数域

对步骤7中处理后的N个子图像数据矩阵统一做如下处理：

取出步骤7中第n个子数据矩阵n＝1，…，N，利用快速傅里叶变换方法对的每一行进行快速傅里叶变换，得到处理后的M行K列的矩阵记为

n＝1，…，N,其中是对的转置矩阵；

步骤9、对波数域各子带图像进行相参叠加

取出步骤8中得到的N个子图像矩阵进行相干叠加得到波数域的等效宽带信号记为I_full， $I_{\mathrm{full}}={\mathrm{\Σ}}_{n=1}^N{\stackrel{\‾}{I_c^{\left(n\right)}}}^{\′},$ n＝1，…，N；

步骤10、将合成的波数域的宽带图像变换到空域

对步骤9中合成的波数域的M个距离向行向量统一做如下处理：

取出步骤9中第mm个距离向的行向量I_full(mm，：)，其中mm＝1，2，…，M，I_full(mm，：)表示I_full的第mm行的所有数据，将其进行快速傅里叶逆变换到空域即可得到高分辨距离向的图像记为Image，Image(mm，：)＝IFFT(I_full(mm，：))，mm＝1，2，…，M，则Image即为最终高分辨距离向图像矩阵。