CN102323583A - 一种超分辨线阵三维合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开一种基于MUSIC算法的线阵三维合成孔径雷达成像方法；首先采用脉冲压缩处理距离向数据，然后对每一个距离切面采用MUSIC算法聚焦场景信息，采用LS算法获得场景点的后向散射信息，最后对处理完的距离切面进行整合并显示出来，就得到了三维重建图。与传统的成像方法如BP算法、MF算法相比，本发明克服了瑞利限和天线长度有限的问题，在不提高天线的长度的同时，利用MUSIC方法来实现旁瓣的减弱，减少了旁瓣串扰对成像质量的影响，提高分辨率，实现了三维成像并提高了成像质量。

Description

一种超分辨线阵三维合成孔径雷达成像方法

技术领域

本发明属于雷达成像技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。

背景技术

合成孔径雷达(SAR)是一种高分辨率的微波成像系统，它依靠雷达和目标之间的相对运动来形成合成阵列来获得横向高分辨率，利用大带宽信号实现纵向高分辨率。三维SAR是基于常规SAR横向维和纵向维的基础之上外加切横向维，它也是依靠雷达和目标之间的相对运动来获得切横向分辨率。

三维合成孔径雷达(SAR)是一种新型合成孔径雷达技术。三维成像是三维成像合成孔径雷达(SAR)区别于其他遥感成像技术的特征。当前关于三维SAR技术的相关系统主要包括干涉SAR(InSAR)系统、曲线SAR(CSAR)系统和线阵SAR(LASAR)系统。与InSAR三维成像技术相比，由于干涉技术基于二维SAR图像，当二维SAR图像存在混叠时，干涉技术无法提取各散射点的高度信息。与CSAR技术相比，线阵SAR技术能有效解决平台运动模式的限制，避免了因平台轨迹精度低所带来的误差。

线阵三维SAR(LASAR)成像系统是一种新体制SAR系统，该系统工作在下视模式下，通过平台的运动和一线性阵列安装在载机的切航迹方向(与高度向-沿航迹方向平面垂直的方向)合成等效二维阵列，实现阵列方向的二维分辨。下视线阵SAR具备以对城区建筑物、街区和山区中的峡谷等地形变化剧烈区域成像的潜力。

由于线阵SAR切航迹向分辨率由线阵SAR切航迹向天线长度决定，通过增加天线长度的方法获得切航迹向高分辨率，大大增加了系统的硬件开销，而且，由于采样点个数有限，传统的匹配滤波器方法收到瑞利限的限制，分辨率很难提高。

多重信号分类(MUSIC)算法是1979年由R.O.Schmidt提出一种超分辨谱估计算法。MUSIC算法应用到线阵SAR有两点不足：(1)由于成像雷达的工作特点，雷达目标成像主要是以单帧雷达图像处理为主的。而MUSIC方法构造协方差矩阵需要回波的统计信息。(2)MUSIC谱是利用子空间的正交性构造的，是伪谱，因此不具有目标的后向散射特性。

发明内容

本发明的目的是为了克服现有的线阵SAR成像方法的切航迹向分辨率受到天线长度的限制，沿航迹向分辨率受到瑞利限的限制的缺点，提出了基于多重信号分类(MUSIC)算法的一种超分辨线阵三维SAR成像方法。该方法在距离向仍然采用脉冲压缩的方法，在切航迹向和沿航迹向采用二维MUSIC算法和最小二乘法获得场景点的位置信息和后向散射信息，降低旁瓣高度和主瓣宽度，从而达到了增加分辨率的目的。

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

定义1、线阵三维成像合成孔径雷达

线阵三维成像合成孔径雷达是将线性阵列天线固定于运动平台上，利用平台飞行速度和线阵长度以合成二维平面阵列，可以对测绘区域进行三维成像的一种合成孔径雷达系统。

定义2、线阵三维成像合成孔径雷达理论分辨率

线阵三维合成孔径雷达理论分辨率是指根据线阵三维合成孔径雷达系统参数，包括发射信号带宽，合成孔径长度以及线阵天线长度决定的三维合成孔径雷达所能达到的最大分辨率。

定义3、三维合成孔径雷达成像空间

三维合成孔径雷达成像空间是指运用成像算法由三维合成孔径雷达数据空间所得到的三维合成孔径雷达的图像空间。

定义4、三维合成孔径雷达数据空间

三维合成孔径雷达数据空间是指三维合成孔径雷达回波数据所构成的回波信号空间。

定义5、飞行孔径与慢时间

飞行孔径是指对于测绘场景中的一个散射点从收发波束共同照射到开始到发射波束或接收波束任意一个照射不到结束收发波束中心所走过的距离。

慢时间是指收发平台飞过一个飞行孔径所需要的时间，由于雷达以一定的重复周期T_r发射接收脉冲，在慢时刻n可以表示为一个离散化的时间变量t_s＝nT_r，n＝0，1，L，N，其中n是自然数。

定义6、延迟时间

延迟时间是发射机发射信号到接收机接收到信号这一段时间，记为τ，它由观测区域到该系统的距离R决定，τ＝(2R)/C，其中，C为光速。

定义7、合成孔径雷达标准距离压缩方法

合成孔径雷达标准距离压缩方法是指利用合成孔径雷达发射信号参数，采用匹配滤波技术对合成孔径雷达的距离向信号进行滤波的过程。详见文献“雷达成像技术”，保铮等编著，电子工业出版社出版。

定义8、去调频方法

去调频方法又称为STRETCH处理法、宽带压缩法或时频转换法，是针对线性调频信号提出的一种信号处理方法，在接收端设置一个参考函数，以场景中心为参考点，将参考函数与接收信号在时域共轭相乘，实现去调频。

定义8、多重信号分类算法(MUSIC)

MUSIC算法是一种基于矩阵特征空间分解的方法。从几何角度讲，信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值(噪声方差)对应的特征向量组成。MUSIC算法就是利用这两个互补空间之间的正交特性来估计空间信号的方位。噪声子空间的所有向量被用来构造谱，所有空间方位谱中的峰值位置对应信号的来波方位。MUSIC算法大大提高了测向分辨率，同时适应于任意形状的天线阵列，但是原型MUSIC算法要求来波信号是不相干的，MUSIC算法的详细流程参见文献“现代数字信号处理及其应用，何子述等著”。

定义9、最小二乘法

最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

定义10、空间平滑算法

空间平滑算法由Evans等人首先提出，利用各子阵间的互相关矩阵，构造协方差矩阵，是一种常用的处理相干源的预处理方法，具体步骤参见文献“AnImproved Spatial Smoothing Technique for Bearing Estimation in a Multipath Environment”。

定义11、滑动窗口方法

一般情况下，滑动窗口方法应用与二维图像(矩阵)空间滤波，可以根据实际需要自行选定窗口的大小，如5*5，7*7等，该方窗在整幅图像上滑动，从而计算出每个像素点滤波后的输出值。具体步骤参见文献“合成孔径雷达成像原理，皮亦鸣等著”。

定义12、合成孔径雷达发射机

合成孔径雷达发射机是指目前合成孔径雷达采用的向观测区域发射电磁信号的系统，主要包括信号发生器、混频器、放大器等模块。

定义13、合成孔径雷达接收机

合成孔径雷达接收机是指目前合成孔径雷达采用的接收观测区域回波的系统，主要包括混频器、放大器、模/数转换器、存储设备等。

定义14、三维合成孔径雷达成像空间

三维合成孔径雷达成像空间是指运用成像算法由三维合成孔径雷达数据空间所得到的三维合成孔径雷达的图像空间。

定义15、Matlab

Matlab是矩阵实验室(Matrix Laboratory)的简称，美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。具体用法详见文献“MATLAB 5手册”，EvaPart-Enander等编著，机械工业出版社出版。

本发明提供了基于多重信号分类(MUSIC)算法的一种超分辨线阵三维合成孔径雷达成像方法，它包括以下几个步骤，如附图2所示：

步骤1、初始化线阵三维成像合成孔径雷达成像系统参数：

初始化成像系统参数包括：平台速度矢量，记做

平台初始位置矢量，记做

雷达发射电磁波的波数，记做K_c，线阵天线各阵元相对平台中心的位置，记做

其中m为天线各阵元序号，为自然数，m＝0，1，...，M，M为线阵天线各阵元总数，雷达发射基带信号的信号带宽，记做B，雷达发射信号脉冲宽度，记做T_P，雷达接收波束持续宽度，记做T_o，雷达接收系统的采样频率，记做f_s，雷达系统的脉冲重复频率，记做PRF，雷达接收系统接收波门相对于发射信号发散波门的延迟，记做T_D，线性阵列天线长度，记做L，雷达的合成孔径长度，记做l；上述参数均为线阵三维成像合成孔径雷达系统的标准参数，其中，雷达发射电磁波的波数K_c，雷达发射基带信号的信号带宽B，雷达发射信号脉冲宽度T_P，雷达接收波门持续宽度T_o，雷达接收系统的采样频率f_s，线性阵列天线长度L，雷达的合成孔径长度l，雷达系统的脉冲重复频率PRF及接收系统接收波门相对于发射信号发散波门的延迟在线阵三维成像合成孔径雷达设计过程中已经确定；其中，平台速度矢量

及平台初始位置矢量

在线阵三维成像合成孔径雷达观测方案设计中已经确定；根据已有的线阵三维成像合成孔径雷达系统和线阵三维成像合成孔径雷达观测方案，线阵三维合成孔径雷达成像方法需要的初始化成像系统参数均为已知。

步骤2、线阵三维成像合成孔径雷达原始数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对合成孔径雷达距离向回波数据进行压缩，得到距离压缩后的线阵三维成像合成孔径雷达数据，记做

距离压缩后的线阵三维成像合成孔径雷达数据

是一个三维数组，包括距离维、沿航迹维、切航迹维；

步骤3、获取线阵三维成像合成孔径雷达图像空间的距离切片

将步骤2中得到的距离压缩后的线阵三维成像合成孔径雷达数据

中距离维上序号为1的二维矩阵提取出，得到第一个距离切片，记作第一个距离切片

是二维矩阵，第一个距离切片的行数为线阵阵元个数M，第一个距离切片

的列数为方位向采样点数N；

同理，将步骤2中得到的距离压缩后的线阵三维成像合成孔径雷达数据

中距离维上序号为2的二维矩阵提取出，得到第二个距离切片，记作

第二个距离切片

是二维矩阵，第二个距离切片

的行数为线阵阵元个数M，第二个距离切片

的列数为方位向采样点数N；

同理，将步骤2中得到的距离压缩后的线阵三维成像合成孔径雷达数据

中距离维上序号为i的二维矩阵提取出，得到第i个距离切片，记作第i个距离切片

是二维矩阵，第i个距离切片

的行数为线阵阵元个数M，第i个距离切片

的列数为方位向采样点数N，i是1～R之间的任意一个距离切片的序号，R为距离向采样个数，即距离切片数；

步骤4、获得第一个距离切片

的成像结果

对步骤3中获得第一个距离切片

进行成像处理，具体方法：

(1)采用去调频方法预处理第一个距离切片

采用传统的去调频方法补偿距离切片

相位中的二次项，具体方法是：

首先，第一个距离切片

中元素的表示形式为：

$s(n,m)=\underset{P}{\mathrm{\Σ}}{\mathrm{\σ}}_p\exp (-j\·2K_c\·R(n,m;p)){\mathrm{\χ}}^R(r_1-2R(n,m;p)),$

其中，p为第p个目标点，σ_p为第p个目标点的后向散射系数，j为虚数单位，K_c为雷达发射电磁波的波数，n为第n个慢时刻，为自然数，n的取值范围为1，2，L，N，m为第m个阵元，为自然数，m的取值范围为1，2，L，M，r₁为第一个距离切片

在距离向的值，χ^R为距离模糊函数，R(n，m；p)为目标点p在慢时刻n到第m阵元的距离，表示为： $R(n,m;p)=\sqrt{{({\mathrm{\ϵ}}_{\mathrm{xp}}-x_n)}^{\text{2}}\text{+}{({\mathrm{\ϵ}}_{\mathrm{yp}}-y_m)}^2+{(H-{\mathrm{\ϵ}}_{\mathrm{yp}})}^2},$ 其中(x_n，y_m，H)为慢时刻n时第m阵元的三维位置信息，第p个目标点的三维坐标信息为ξ＝(ξ_xp，ξ_yp，ξ_zp)^t，通过菲涅尔近似表示为， $R(n,m;p)=R_0(1+\frac{{({\mathrm{\ϵ}}_{\mathrm{xp}}-x_n)}^2+{y_m}^2-2y_m{\mathrm{\ϵ}}_{\mathrm{yp}}}{2R_0}),$ 其中R₀为

是R(n，m；p)在y-z平面的投影；

然后，将第一个距离切片

中的元素乘以乘子得到二维复正弦信号，记作

二维复正弦信号

中的元素表示为， $s(n,m)=\underset{P}{\mathrm{\Σ}}{\mathrm{\σ}}_p\exp (-j\·K_c\·(2\·x_n\frac{{\mathrm{\ϵ}}_{\mathrm{xp}}}{R_0}+2\·y_m\frac{{\mathrm{\ϵ}}_{\mathrm{yp}}}{R_0}));$

(2)采用空间滑动平均方法构造协方差矩阵

采用传统的空间滑动平均方法对二维复正弦信号

进行处理，具体方法是：

首先，设置二维滑动窗口的维数分别为P₁和P₂，其中，滑动窗口的行数P₁小于距离切片的行数M，滑动窗口的列数P₂小于距离切片的列数N，通过滑动窗口方法获得多个P₁×P₂子孔径二维矩阵，如附图3所示；

然后，对滑动窗口获得的每个子孔径二维矩阵P₁×P₂的二维回波数据进行一对一重排使子孔径二维矩阵P₁×P₂转化为P₁P₂×1的一维矢量，记为S_k，其中k为自然数，取值范围为k＝0，1，...，K，K为子孔径的个数：K＝(M-P₁+1)(N-P₂+1)，一对一重排的具体方法是：将矩阵中下标为(m，n)的元素的值作为矢量中下标为(mn，1)的元素值；

最后，构造P₁P₂×P₁P₂的协方差矩阵R_xx，构造公式为

为第k个子孔径的协方差矩阵：

()^H是共轭转置符号，J是转换矩阵，定义为 $J\left(\begin{array}{ccccc}0& 0& L& 0& 1\\ 0& 0& L& 1& 0\\ M& M& N& M& M\\ 0& 1& L& 0& 0\\ 1& 0& L& 0& 0\end{array}\right);$

(3)利用MUSIC算法对场景成像

首先，对协方差矩阵R_xx进行传统的特征值分解方法进行特征值分解，得到特征值λ₁，λ₂，L，λ_P1×P2，P₁×P₂-F个最小特征值对应的特征向量为u₁，u₂，L，

其中F为目标散射点个数；构造噪声子空间矩阵G，噪声子空间矩阵G的每一列为最小特征值对应的特征向量；

然后，在场景范围内遍历所有点位置坐标，计算MUSIC谱函数 $P_{\mathrm{MUSIC}}=\frac{a{({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}^{H}a({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}{a{({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}^{H}G{G}^{H}a({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}$ 的值，其中， $a{(\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)=\exp (-j\·K_c\·(2x_n\·\frac{{\mathrm{\ϵ}}_x}{R_0}+2y_m\·\frac{{\mathrm{\ϵ}}_y}{R_0}))$ 为行坐标x，纵坐标为y的场景采样点的位置(ε_x，ε_y)对应的导向矢量，谱峰位置对应场景点坐标就是散射点的坐标信息的估计值；

最后，将对应的场景采样点的谱函数的值赋值给坐标(x，y)，得到第一个距离切片的成像结果，记为

为一个二维矩阵，二维矩阵

行数为场景横向采样点数，二维矩阵

列数为场景纵向采样点数；

(4)采用LS算法获得场景点的后向散射信息

首先，根据阵元位置信息构造数据矩阵A，数据矩阵A的行数为线阵阵元个数M与慢时间的采样个数N的乘积，列数为总场景坐标个数，数据矩阵A中元素的值为单个场景采样点到单个阵元回波信号的相位信息；

然后，将第一个距离切片一对一重排，一对一重排的具体方法是：将矩阵中下标为(m，n)的元素的值作为矢量中下标为(mn，1)的元素值；得到观测向量y；

最后采用LS算法对观测向量y进行处理，得到后向散射系数矢量φ，处理公式为φ＝(A^HA)^-1A^Hy，φ中非0元素的值对应场景散射点的后向散射系数，获得场景点的后向散射信息σ_p，p＝1，L，F，p为自然数；

将场景点的后向散射信息σ_p赋值于第i个距离切片的距离切面的成像结果

中，i是1～R之间的任意一个距离切片的序号；具体方法是：

将第i个距离切片的成像结果

中第1个散射点对应的矩阵元素用σ₁替换，…，将第p个散射点对应的矩阵元素用σ_p替换，…，将第F个散射点对应的矩阵元素用σ_F替换；得到具有后向散射信息的距离切片的成像结果

步骤5、全场景成像

采用步骤4的方法处理第二个距离切片，得到第二个距离切片的成像结果，......，同理，采用步骤4的方法处理第i个距离切片，得到第i个距离切片的成像结果，i是1～R之间的任意一个距离切片的序号；最后，采用步骤4的方法处理第R个距离切片，得到第R个距离切片的成像结果，R为距离向采样个数；

建立R×N_x×N_y三维数据矩阵J，其中N_x×N_y为沿航迹向和切航迹向的场景采样个数，J的第一维为距离向，第二维为切航迹向，第三维为沿航迹向；

将距离切片的成像结果依次赋值给三维数据矩阵J中对应的两维矩阵，具体方法是：第一个距离切片中的元素H₁(x，y)赋值给三维数组J中的元素J(1，x，y)，......，同理，第i个距离切片中的元素H_i(x，y)赋值给三维数组J中的元素J(i，x，y)，i是1～R之间的任意一个距离切片的序号；第R个距离切片中的元素H_R(x，y)赋值给三维数组J中的元素J(R，x，y)，构造最终的成像结果矩阵；

经过上述步骤，即可得到基于MUSIC算法的超分辨率线阵三维成像合成孔径雷达图像；

本发明的创新点：在于针对现有技术中的三维线阵合成孔径雷达的成像过程中存在的旁瓣串扰现象导致系统成像质量下降，而在目标散射机理研究等应用中，对测量精度提出了很高的要求，本发明提供了一种新的基于MUSIC算法的线阵三维合成孔径雷达成像方法；该方法用MUSIC算法获得场景点的位置信息，用LS算法获得场景点的后向散射信息，降低旁瓣和主瓣宽度，提高分辨率，最终达到从而提高成像质量的目的。

本发明的基本原理：是用MUSIC算法代替原本的匹配滤波算法，首先采用脉冲压缩处理距离向数据，然后对每一个距离切面采用MUSIC算法聚焦场景信息，采用LS算法获得场景点的后向散射信息，从而达到主瓣宽度减小和旁瓣降低的目的，提高了成像质量，最后对处理完的距离切面进行整合并显示出来，就得到了三维重建图。

本发明优点：针对线阵SAR切航迹向和沿航迹向的分辨率受到线阵长度和瑞利限限制，利用MUSIC方法来实现旁瓣的减弱，减少了旁瓣串扰对成像质量的影响，从而提高了分辨率；此外，它还继承了线阵SAR成像系统的天线相位中心控制精度高和对复杂地形成像效果好的优点；本发明与传统的基于FFT原理的成像方法(BP算法，MF算法)相比，克服了瑞利限和天线长度有限的问题，在不提高天线的长度的同时，提高分辨率，实现了三维成像并提高了成像质量。

附图说明

图1为全激励线阵三维合成孔径雷达几何结构图

其中，为飞行器的速度矢量；X、Y、Z为三维坐标系，P_ω为一个场景点。

图2是本发明的结构流程图

图3是滑动平均方法原理图

图4是本发明具体实施方式成像仿真参数表

具体实施方式

本发明主要采用仿真实验的方式进行验证该系统模型的可行性，所有步骤、结论都在MATLAB7.0上验证正确。具体实施步骤如下：

步骤一：系统参数的设定

本具体实施方式所采用的系统参数详见图4。

步骤二、线阵三维成像合成孔径雷达原始数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对合成孔径雷达距离向回波数据

进行压缩，得到距离压缩后的线阵三维成像合成孔径雷达数据，记做

距离压缩后的线阵三维成像合成孔径雷达数据

是一个三维数组，包括距离维、沿航迹维、切航迹维；

步骤三、获取线阵三维成像合成孔径雷达图像空间的第一个距离切片步骤二中压缩后的回波数据

为一个三维数组，将

中距离维上序号为1的二维矩阵提取出，记作的行数为慢时间采样个数，m＝1，L，20，列数为线阵阵元个数，n＝1，L，20。

步骤四、获得第一个距离切片的成像结果

对步骤三中获得距离切片

进行成像处理，如附图2所示；

(1)利用去调频方法预处理距离切片

利用传统的去调频方法补偿距离切片

相位中的二次项，具体是：首先，距离切片

中元素的表示形式为：

$s(n,m)=\underset{P}{\mathrm{\Σ}}{\mathrm{\σ}}_p\exp (-j\·2K_c\·R(n,m;p)){\mathrm{\χ}}^R(r_1-2R(n,m;p)),$ 其中p为第p个目标点，目标点的三维坐标信息为ξ＝(ξ_xp，ξ_yp，ξ_zp)^t，σ_p为第p个目标点的后向散射系数，K_c为雷达发射电磁波的波数，n为第n个慢时刻，m为第m个阵元，r₁为第一个距离切片在距离向的值，R(n，m；p)为目标点p在慢时刻n到第m阵元的距离，表示为： $R(n,m;p)=\sqrt{{({\mathrm{\ϵ}}_{\mathrm{xp}}-x_n)}^{\text{2}}\text{+}{({\mathrm{\ϵ}}_{\mathrm{yp}}-y_m)}^2+{(H-{\mathrm{\ϵ}}_{\mathrm{yp}})}^2},$ 其中(x_n，y_m，H)为慢时刻n时第m阵元的三维位置信息，通过菲涅尔近似表示为， $R(n,m;p)=R_{0g}(1+\frac{{({\mathrm{\ϵ}}_{\mathrm{xp}}-x_n)}^2+{y_m}^2-2y_m{\mathrm{\ϵ}}_{\mathrm{yp}}}{2R_0}),$ 其中R₀为

是R(n，m；p)在y-z平面的投影，然后将

中的元素乘以乘子

转化为二维复正弦信号模型记作中的元素可以表示为， $s(n,m)=\underset{P}{\mathrm{\Σ}}{\mathrm{\σ}}_p\exp (-j\·K_c\·(2\·x_n\frac{{\mathrm{\ϵ}}_{\mathrm{xp}}}{R_0}+2\·y_m\frac{{\mathrm{\ϵ}}_{\mathrm{yp}}}{R_0}));$

(2)利用空间滑动平均方法构造协方差矩阵

采用空间滑动平均方法对二维复正弦信号

进行处理，具体方式是：首先，设置二维滑动窗口的维数分别为14和14，其中，滑动窗口的维数均小于距离切片的维数20，20，然后，对滑动窗口获得的每个子孔径14×14的二维回波数据进行一对一重排使其转化为156×1的一维矢量，即元素下标为(m，n)在矢量中表示为(mn，1)，最终，构造156×156的协方差矩阵R_xx；

(3)利用MUSIC算法对场景成像

首先，对协方差矩阵R_xx进行特征值分解，得到特征值λ₁，λ₂，L，λ₁₅₆，156-10个最小特征值对应的特征向量为u₁，u₂，L，u₁₄₆，其中10为目标散射点个数，构造噪声子空间矩阵G，G的每一列为最小特征值对应的特征向量，然后，在场景范围内扫描所有点位置坐标，计算MUSIC谱函数

的值，其中， $a{(\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)=\exp (-j\·K_c\·(2x_n\·\frac{{\mathrm{\ϵ}}_x}{R_0}+2y_m\·\frac{{\mathrm{\ϵ}}_y}{R_0}))$ 为场景采样点(ε_x，ε_y)对应的导向矢量，谱峰位置对应场景点坐标就是散射点的坐标信息的估计值，获得距离切面的成像结果，记为

(4)利用LS算法获得场景点的后向散射信息

首先，根据阵元位置信息构造数据矩阵A，A的行数为虚拟面阵的阵元个数，即线阵阵元个数M与慢时间的采样个数T的乘积，列数为总场景坐标个数，本发明中设定场景范围x轴范围为[-32m，32m]，y轴范围[-32m，32m]，采样点在x轴和y轴间隔都为为1m，这样场景范围内采样点个个数为65×65个，A中元素的值为单个场景点到单个阵元回波信号的相位信息；然后，将距离切片

一对一重排得到观测向量y，矢量化方法和步骤三中一致，最后利用LS算法，φ＝(A^HA)^-1A^Hy，φ中非0元素的值对应场景散射点的后向散射系数，获得场景点的后向散射信息σ_p，p＝1，L，10；

将场景点的后向散射信息赋值于MUSIC算法获得的成像结果

中最终得到距离切片的切片成像结果

步骤五、全场景成像

采用步骤四的方法处理第二个距离切片，同理，处理第i个距离切片，最后，处理第512个距离切片，512为距离向采样个数；建立512×65×65三维数据矩阵J其中65×65为沿航迹向和切航迹向的场景采样个数，J的第一维为距离向，第二维为切航迹向，第三维为沿航迹向，将距离切片的成像结果依次赋值给三维数组J中对应的两维矩阵，例如，第一个距离切片中的元素

赋值给三维数组J中的元素J(1，x，y)，同理，第i个距离切片中的元素赋值给三维数组J中的元素J(i，x，y)，第512个距离切片中的元素

赋值给三维数组J中的元素J(512，x，y)，构造最终的成像结果矩阵；

经过上述操作，即可得到基于MUSIC算法的超分辨率线阵三维成像合成孔径雷达图像。

通过本发明具体实施方式的仿真及测试，本发明所提出的基于MUSIC算法的线阵三维合成孔径雷达成像方法，与现有的三维合成孔径雷达成像图像相比，它实现降低旁瓣和主瓣宽度从而提高分辨率的目的。

Claims (1)

1.一种超分辨线阵三维合成孔径雷达成像方法，它包括以下几个步骤：

步骤1、初始化线阵三维成像合成孔径雷达成像系统参数：

初始化成像系统参数包括：平台速度矢量，记做

平台初始位置矢量，记做

雷达发射电磁波的波数，记做K_c，线阵天线各阵元相对平台中心的位置，记做

其中m为天线各阵元序号，为自然数，m＝0，1，...，M，M为线阵天线各阵元总数，雷达发射基带信号的信号带宽，记做B，雷达发射信号脉冲宽度，记做T_p，雷达接收波束持续宽度，记做T_o，雷达接收系统的采样频率，记做f_s，雷达系统的脉冲重复频率，记做PRF，雷达接收系统接收波门相对于发射信号发散波门的延迟，记做T_D，线性阵列天线长度，记做L，雷达的合成孔径长度，记做l；上述参数均为线阵三维成像合成孔径雷达系统的标准参数，其中，雷达发射电磁波的波数K_c，雷达发射基带信号的信号带宽B，雷达发射信号脉冲宽度T_P，雷达接收波门持续宽度T_o，雷达接收系统的采样频率f_s，线性阵列天线长度L，雷达的合成孔径长度l，雷达系统的脉冲重复频率PRF及接收系统接收波门相对于发射信号发散波门的延迟在线阵三维成像合成孔径雷达设计过程中已经确定；其中，平台速度矢量

及平台初始位置矢量在线阵三维成像合成孔径雷达观测方案设计中已经确定；根据已有的线阵三维成像合成孔径雷达系统和线阵三维成像合成孔径雷达观测方案，线阵三维合成孔径雷达成像方法需要的初始化成像系统参数均为已知；

步骤2、线阵三维成像合成孔径雷达原始数据进行距离压缩

采用传统的合成孔径雷达标准距离压缩方法对合成孔径雷达距离向回波数据

进行压缩，得到距离压缩后的线阵三维成像合成孔径雷达数据，记做

距离压缩后的线阵三维成像合成孔径雷达数据

是一个三维数组，包括距离维、沿航迹维、切航迹维；

步骤3、获取线阵三维成像合成孔径雷达图像空间的距离切片

将步骤2中得到的距离压缩后的线阵三维成像合成孔径雷达数据

中距离维上序号为1的二维矩阵提取出，得到第一个距离切片，记作

第一个距离切片

是二维矩阵，第一个距离切片的行数为线阵阵元个数M，第一个距离切片

的列数为方位向采样点数N；

同理，将步骤2中得到的距离压缩后的线阵三维成像合成孔径雷达数据

中距离维上序号为2的二维矩阵提取出，得到第二个距离切片，记作

第二个距离切片

是二维矩阵，第二个距离切片的行数为线阵阵元个数M，第二个距离切片

的列数为方位向采样点数N；

同理，将步骤2中得到的距离压缩后的线阵三维成像合成孔径雷达数据

中距离维上序号为i的二维矩阵提取出，得到第i个距离切片，记作

第i个距离切片是二维矩阵，第i个距离切片

的行数为线阵阵元个数M，第i个距离切片的列数为方位向采样点数N，i是1～R之间的任意一个距离切片的序号，R为距离向采样个数，即距离切片数；

步骤4、获得第一个距离切片

的成像结果

对步骤3中获得第一个距离切片

进行成像处理，具体方法：

(1)采用去调频方法预处理第一个距离切片

采用传统的去调频方法补偿距离切片

相位中的二次项，具体方法是：

首先，第一个距离切片

中元素的表示形式为：

$s(n,m)=\underset{P}{\mathrm{\Σ}}{\mathrm{\σ}}_p\exp (-j\·2K_c\·R(n,m;p)){\mathrm{\χ}}^R(r_1-2R(n,m;p)),$

其中，p为第p个目标点，σ_p为第p个目标点的后向散射系数，j为虚数单位，K_c为雷达发射电磁波的波数，n为第n个慢时刻，为自然数，n的取值范围为1，2，L，N，m为第m个阵元，为自然数，m的取值范围为1，2，L，M，r₁为第一个距离切片

在距离向的值，χ^R为距离模糊函数，R(n，m；p)为目标点p在慢时刻n到第m阵元的距离，表示为： $R(n,m;p)=\sqrt{{({\mathrm{\ϵ}}_{\mathrm{xp}}-x_n)}^{\text{2}}\text{+}{({\mathrm{\ϵ}}_{\mathrm{yp}}-y_m)}^2+{(H-{\mathrm{\ϵ}}_{\mathrm{yp}})}^2},$ 其中(x_n，y_m，H)为慢时刻n时第m阵元的三维位置信息，第p个目标点的三维坐标信息为ξ＝(ξ_xp，ξ_yp，ξ_zp)t，通过菲涅尔近似表示为， $R(n,m;p)=R_0(1+\frac{{({\mathrm{\ϵ}}_{\mathrm{xp}}-x_n)}^2+{y_m}^2-2y_m{\mathrm{\ϵ}}_{\mathrm{yp}}}{2R_0}),$ 其中R₀为

是R(n，m；p)在y-z平面的投影；

然后，将第一个距离切片

中的元素乘以乘子

得到二维复正弦信号，记作

二维复正弦信号

中的元素表示为， $s(n,m)=\underset{P}{\mathrm{\Σ}}{\mathrm{\σ}}_p\exp (-j\·K_c\·(2{\·x}_n\frac{{\mathrm{\ϵ}}_{\mathrm{xp}}}{R_0}+2\·y_m\frac{{\mathrm{\ϵ}}_{\mathrm{yp}}}{R_0}));$

(2)采用空间滑动平均方法构造协方差矩阵

采用传统的空间滑动平均方法对二维复正弦信号

进行处理，具体方法是：

首先，设置二维滑动窗口的维数分别为P₁和P₂，其中，滑动窗口的行数P₁小于距离切片的行数M，滑动窗口的列数P₂小于距离切片的列数N，通过滑动窗口方法获得多个P₁×P₂子孔径二维矩阵；

然后，对滑动窗口获得的每个子孔径二维矩阵P₁×P₂的二维回波数据进行一对一重排使子孔径二维矩阵P₁×P₂转化为P₁P₂×1的一维矢量，记为S_k，其中k为自然数，取值范围为k＝0，1，...，K，K为子孔径的个数：K＝(M-P₁+1)(N-P₂+1)，一对一重排的具体方法是：将矩阵中下标为(m，n)的元素的值作为矢量中下标为(mn，1)的元素值；

最后，构造P₁P₂×P₁P₂的协方差矩阵R_xx，构造公式为

为第k个子孔径的协方差矩阵：

()^H是共轭转置符号，J是转换矩阵，定义为 $J\left(\begin{array}{ccccc}0& 0& L& 0& 1\\ 0& 0& L& 1& 0\\ M& M& N& M& M\\ 0& 1& L& 0& 0\\ 1& 0& L& 0& 0\end{array}\right);$

(3)利用MUSIC算法对场景成像

首先，对协方差矩阵R_xx进行传统的特征值分解方法进行特征值分解，得到特征值λ₁，λ₂，L，λ_P1×P2，P₁×P₂-F个最小特征值对应的特征向量为u₁，u₂，L，

其中F为目标散射点个数；构造噪声子空间矩阵G，噪声子空间矩阵G的每一列为最小特征值对应的特征向量；

然后，在场景范围内遍历所有点位置坐标，计算MUSIC谱函数 $P_{\mathrm{MUSIC}}=\frac{a{({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}^{H}a({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}{a{({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}^{H}G{G}^{H}a({\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)}$ 的值，其中， $a{(\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y)=\exp (-j\·K_c\·(2x_n\·\frac{{\mathrm{\ϵ}}_x}{R_0}+2y_m\·\frac{{\mathrm{\ϵ}}_y}{R_0}))$ 为行坐标x，纵坐标为y的场景采样点的位置(ε_x，ε_y)对应的导向矢量，谱峰位置对应场景点坐标就是散射点的坐标信息的估计值；

最后，将对应的场景采样点的谱函数的值赋值给坐标(x，y)，得到第一个距离切片的成像结果，记为为一个二维矩阵，二维矩阵

行数为场景横向采样点数，二维矩阵

列数为场景纵向采样点数；

(4)采用LS算法获得场景点的后向散射信息

首先，根据阵元位置信息构造数据矩阵A，数据矩阵A的行数为线阵阵元个数M与慢时间的采样个数N的乘积，列数为总场景坐标个数，数据矩阵A中元素的值为单个场景采样点到单个阵元回波信号的相位信息；

然后，将第一个距离切片

一对一重排，一对一重排的具体方法是：将矩阵中下标为(m，n)的元素的值作为矢量中下标为(mn，1)的元素值；得到观测向量y；

最后采用LS算法对观测向量y进行处理，得到后向散射系数矢量φ，处理公式为φ＝(A^HA)^-1A^Hy，φ中非0元素的值对应场景散射点的后向散射系数，获得场景点的后向散射信息σ_p，p＝1，L，F，p为自然数；

将场景点的后向散射信息σ_p赋值于第i个距离切片的距离切面的成像结果

中，i是1～R之间的任意一个距离切片的序号；具体方法是：

将第i个距离切片的成像结果

中第1个散射点对应的矩阵元素用σ₁替换，…，将第p个散射点对应的矩阵元素用σ_p替换，…，将第F个散射点对应的矩阵元素用σ_F替换；得到具有后向散射信息的距离切片的成像结果

步骤5、全场景成像

采用步骤4的方法处理第二个距离切片，得到第二个距离切片的成像结果，......，同理，采用步骤4的方法处理第i个距离切片，得到第i个距离切片的成像结果，i是1～R之间的任意一个距离切片的序号；最后，采用步骤4的方法处理第R个距离切片，得到第R个距离切片的成像结果，R为距离向采样个数；

建立R×N_x×N_y三维数据矩阵J，其中N_x×N_y为沿航迹向和切航迹向的场景采样个数，J的第一维为距离向，第二维为切航迹向，第三维为沿航迹向；

将距离切片的成像结果依次赋值给三维数据矩阵J中对应的两维矩阵，具体方法是：第一个距离切片中的元素H₁(x，y)赋值给三维数组J中的元素J(1，x，y)，......，同理，第i个距离切片中的元素H_i(x，y)赋值给三维数组J中的元素J(i，x，y)，i是1～R之间的任意一个距离切片的序号；第R个距离切片中的元素H_R(x，y)赋值给三维数组J中的元素J(R，x，y)，构造最终的成像结果矩阵；经过上述步骤，即可得到基于MUSIC算法的超分辨率线阵三维成像合成孔径雷达图像。

Images