本发明公开了一种基于单快拍MUSIC算法的距离向处理方法，它是基于运用单快拍MUSIC算法估计目标的位置，首先对原始回波信号进行去斜处理，然后根据单快拍MUSIC算法进行构造自相关矩阵、特征分解等过程；计算MUSIC谱并进行谱峰搜索，这样就能确定目标所在的距离单元格；最后使用最小二乘法计算目标所在距离单元格的值。本发明与传统的通过脉冲压缩方法进行距离向处理相比，本发明仅需要单个回波信号即可使用MUSIC算法进行距离向的处理，分辨率有较大的提高，能够实现距离向超分辨能力，本发明特别适用于目标较少，精度要求高的情况。
1.一种基于单快拍MUSIC算法的距离向处理方法，其特征是它包括以下几个步骤：步骤1、初始化SAR系统参数：
初始化SAR系统参数包括：雷达工作中心频率，记做fc；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做Br；雷达发射信号脉冲宽度，记做Tr；雷达发射信号的调频斜率，记做Kr；雷达接收系统的采样频率，记做fs；电磁波在空气中的传播速度，记做C；连续采样时间为TP，n为距离向时间序号，n＝1,2,…,Nr，n为自然数，Nr为距离向采样点总数；上述参数均为SAR系统标准参数，其中雷达中心频率fc，雷达载频波长λ，雷达发射基带信号的信号带宽Br，雷达发射信号脉冲宽度Tr，雷达发射信号调频斜率Kr，雷达接收系统的采样频率fs，波门延时为τ0，在SAR系统设计过程中已经确定；根据SAR成像系统方案和观测方案，距离向处理所需要的初始化成像系统参数均为已知；
在第t个距离向采样得到的原始回波数据，记为SaE(t)，t＝1,2,…,TP，t为自然数；在实际线阵SAR系统中，原始回波数据SaE(t)可由雷达系统数据接收机提供；在仿真线阵SAR成像过程中，原始回波数据SaE(t)根据雷达系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在雷达回波数据进行距离向处理之前，原始回波数据SaE(t)均已知；
初始化场景大小为[-X,X]，将原点设置场景中心，定义场景中有M个散射点，且目标都在距离向接收范围内，这M个散射点相对于场景中心位置为x1,x2,...,xM，这M个点的散射系数分别记为 通过公式 得到M个散射点所对应的回波时延，分别为其中，x为散射点相对于场景中心的距离，C为光速，原始回波数据SaE(t)通过采用标准合成孔径雷达原始回波仿真方法产生得到，即其中，A为回波信号的增益，这里将A设为1，根据采样频率fs对连续回波信号SaE(t)进行采样，得到离散回波信号SAE(n)，即其中，Nr为距离向采样点数；
其中，n＝1,2,…,Nr，n为自然数；上式中所使用的参数均为步骤1中初始化的参数；
通过标准的取共轭方法对参考信号href(n)进行处理，得到共轭参考信号通过将步骤2中的回波信号SAE(n)与共轭参考信号 相乘，实现去除回波信号相位中二次项，得到处理后的信号SE(n)，这个过程用公式表示为其中，n＝1,2,…,Nr，n为自然数， 为M个点的散射系数， 为M个散射点所对应的回波时延分，Kr为雷达发射信号的调频斜率，fs为雷达接收系统的采样频率，fc为雷达工作中心频率，SAE(n)为步骤2中得到的离散回波信号；
初始化一个Nr×Nr维零矩阵Ax，矩阵Ax的结构如下：矩阵Ax中的元素ax,y均为0，其中，x，y为自然数，x＝1,2,…,Nr，y＝1,2,…,Nr；
回波信号SE(n)表示为向量形式SE＝[SE(1),SE(2),…,SE(Nr)]，其中，SE(n)为步骤3中得到的n时刻的回波信号，n＝1,2,…,Nr；
选择向量SE＝[SE(1),SE(2),…,SE(Nr)]中第1个到第Nr-m+1个元素，这Nr-m+1个元素定义为向量SE(m,1)＝[SE(1),SE(2),…,SE(Nr-m+1)]，然后将向量SE(m,1)填入矩阵Ax中的第m行，第m列到第Nr列的元素位置，得到矩阵Ax(m,1)，其中，Ax是步骤4.1.1中定义的矩阵，矩阵Ax(m,1)的结果如下其中，SE(n)为步骤4.1.1中向量的元素SE＝[SE(1),SE(2),…,SE(Nr)]，n＝1,2,…,Nr-m+
1；ax,y为矩阵Ax中的元素，其中，x，y为自然数，x＝1,2,…,Nr，y＝1,2,…,Nr；
选择向量SE＝[SE(1),SE(2),…,SE(Nr)]中第2个到第m个元素，这m-1个元素定义为向量SE(m,2)＝[SE(2),SE(3),…,SE(m)]，然后将向量SE(m,2)中元素进行反向排序得到向量SE(m,-2)＝[SE(m),SE(m-1),…,SE(2)]，采用标准的取共轭的方法对向量SE(m,-2)进行处理得到共轭向量SE*(m,-2)＝[SE*(m),SE*(m-1),…,SE*(2)]，最后，将向量SE*(m,-2)填入矩阵Ax(m,1)的第m行，第1列到第m列的位置，得到矩阵Ax(m,2)，其中，Ax(m,1)是步骤4.1.2中的矩阵，矩阵Ax(m,2)的结果如下其中，SE(n)为向量SE＝[SE(1),SE(2),…,SE(Nr)]的元素，n＝1,2,…,Nr-m+1；SE*(n)为向量SE*(m,-2)＝[SE*(m),SE*(m-1),…,SE*(2)]中的元素，其中n＝2,3,…,m；ax,y为矩阵Ax中的元素，其中，x，y为自然数，x＝1,2,…,Nr，y＝1,2,…,Nr；
通过步骤4.1.2和步骤4.1.3完成对矩阵Ax第m行处理，得到矩阵Ax(m,2)，其中，Ax是步骤
4.1.1中定义的矩阵，Ax(m,2)是步骤4.1.3中定义的矩阵；
按照步骤4.1.2和步骤4.1.3的方法对矩阵Ax进行处理，详细过程为：Ax作为初始矩阵，首先对矩阵Ax的第一行进行处理，得到矩阵Ax(1,2)，此时只执行步骤4.1.2，然后对矩阵Ax(1,2)的第二行进行处理，即执行步骤4.1.2和步骤4.1.3，得到Ax(2,2)，依次进行下去，最后对矩阵 的第Nr行进行处理，得到 矩阵 的结构如下：将矩阵 命名为协方差矩阵Rx，采用标准的共轭转置方法得到共轭转置矩阵通过公式 得到自相关矩阵R；
根据步骤4.1.4中得到的自相关矩阵R，使用标准的特征分解方法对其进行处理，得到特征值集合D'和对应的特征向量集合V'，且特征值集合D'与特征向量集合V'是一一对应的；步骤4.2.2、提取特征向量：对特征值集合D'中的元素从大到小进行排序，得到特征值集合 其中， 为步骤4.2.1中得到的特征值集合D'中的元素，且 提取特征值集合D中最后Nr-K个元素，得到特征值集合 其中dn为特征值集合D，n为自然数，n＝K+1,K+2,...,Nr；
从特征向量集合V'中提取出Nr-K个与特征值集合D”对应的特征向量集合其中， 为特征向量集合V'中的元素，且
采用标准的施密特正交化方法对特征向量集合 进行处理，得到特征向量集合 其中un为使用标准的施密特正交化方法处理后的特征向量，n为自然数，n＝K+1,K+2,...,Nr；
根据标准的MUSIC算法和步骤4.2.2中得到的特征向量集合 通过公式 得到噪声子空间矩阵矩阵
连续采样时间为TP，采样点数设置为常数Nr，连续采样时间TP除以采样点数Nr得到采样间隔dτ；从τ0时刻开始到TP+τ0时刻结束，每过dτ时长记录一次当前时刻的值，将该时刻的值记为τn，得到时刻集合 其中，n＝1,2,...,Nr，τ0为波门延时，再将时刻集合 按照下标从小到大排序，得到离散时间向量τ
其中，τn为时间轴第n个采样点的时刻，n为自然数，n＝1,2,...,Nr，τ0为波门延时，即信号通过接收机系统自身的延时，波门延时τ0为已知量；
中τn为离散时间轴，n为自然数，n＝1,2,...,Nr；
其中，Kr为信号调频斜率，fS为雷达接收系统的采样频率，τk为离散时间，e为自然对数的底数，j表示虚部，-jnωm表示系数为n，频率为ωm的相位，其中n＝1,2,...,Nr，m＝1,
2,......,K；即，-jω1表示系数为1、频率为ω1的相位，-jω2表示系数为1、频率为ω2的相位，-jωK表示系数为1、频率为ωK的相位；-j(Nr-1)ω1表示系数为Nr-1、频率为ω1的相位，-j(Nr-1)ω2表示系数为Nr-1、频率为ω2的相位，-j(Nr-1)ωK表示系数为Nr-1、频率为ωK的相位；
对谱PMUSIC'(τ)中的元素PMUSIC'(τn)从大到小依次进行排序，其中n＝1,2,...,Nr，得到向量 PMUSIC'(τn')为PMUSIC'(τ)中的元素，其中n＝1,2,...,Nr，且 其中n＝1,2,...,Nr，根据目标个数K得到门限θ＝PMUSIC'(τM')，其中M为自然数，且K
根据步骤6.2中得到的门限θ和目标个数K，采用MATLAB中标准的谱峰搜索方法对PMUSIC(τ)1进行处理，得到PMUSIC(τ)1中最大的K个谱峰所对应的时延量τk；
然后通过公式 得到延时量τk所对应的目标所在单元格位置IDk，其中，τ0为波门延时，dτ为采样间隔；
通过公式 得到最小二乘估计矩阵中第n行，第k列元素的相位其中，n＝1,2,...,Nr，τk为延时量，k＝1,2,…,K，Kr为信号调频斜率，fS为雷达接收系统的采样频率；
其中， 为步骤7.1中得到的相位，n＝1,2,...,Nr，τk为延时量，k＝1,2,…,K，Kr为信号调频斜率，fS为雷达接收系统的采样频率，e为自然对数的底数，j表示虚部； 表示第n个采样点，延时为τk的相位，其中n＝1,2,...,Nr，m＝1,2,......,K；即，表示第1个采样点，延时为τ1的相位， 表示第1个采样点，延时为τ2的相位， 表示第1个采样点，延时为τK的相位； 表示第2个采样点，延时为τ1的相位， 表示第2个采样点，延时为τ2的相位， 表示第2个采样点，延时为τK的相位； 表示第Nr个采样点，延时为τ1的相位， 表示第Nr个采样点，延时为τ2的相位， 表示第Nr个采样点，延时为τK的相位；
采用标准的共轭转置方法得到最小二乘估计矩阵L1的最小二乘估计的共轭转置矩阵通过最小二乘估计矩阵L1左乘最小二乘估计的共轭转置矩阵 得到矩阵然后采用标准的矩阵求逆方法得到矩阵γ0的逆矩阵
通过公式 计算得到IDk处的值，其中，L1为最小二乘估计矩阵， 为最小二乘估计矩阵的共轭转置矩阵，SE为步骤3中去斜处理的回波数据向量SE(n)，IDk为步骤6中计算出的目标所在单元格的位置；
至此，已完成目标位置和幅度的估计，完成了距离向处理。
一种基于单快拍MUSIC算法的距离向处理方法
技术领域
[0001] 本技术发明属于雷达技术领域，它特别涉及了合成孔径雷达(SAR)成像技术领域。
背景技术
[0002] 合成孔径雷达(SAR)是一种高分辨率微波成像雷达，具有全天时和全天候工作的优点，已被广泛应用各个领域，如地形测绘、制导、环境遥感和资源勘探等。SAR应用的重要前提和信号处理的主要目标是通过成像算法获取高分辨、高精度的微波图像。其获得的高分辨微波图像已经被广泛应用于诸多领域，如生成数字高程图、观测火山活动和洪灾情况、监测陆地和海洋交通等。
[0003] 距离向处理是SAR成像过程中对信号进行的第一步处理方式。通过方位向处理，获得目标的一维图像。距离向处理传统方法有匹配滤波、去斜处理等方法，详见“皮亦鸣.杨建宇等.合成孔径雷达成像原理[M].电子科技大学出版社,2007”。
[0004] 阵列信号处理(ASP)作为信号处理的一个重要分支是在空间不同位置上以一定的规律布觉传感器，形成传感器阵列。阵列的配置由两部分内容组成：一是每个阵元的天线方向图；二是阵列的几何结构(即阵元的物理位置)，通常我们将其分成线性阵列、平面阵列和立体阵列三类。阵列信号处理被广泛应用于雷达、通信、声纳、测向、地震探测、射电天文和电子医疗工程等领域。
[0005] MUSIC算法是一种基于矩阵特征空间分解的方法。信号处理的观测空间分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值(噪声方差)对应的特征向量组成。通常使用MUSIC算法进行波达方向(DOA)估计。
[0006] 目前，用于SAR成像的MUSIC算法有部分研究，有的是需要先聚焦成像，再使用MUSIC算法提高分辨率，详见“张平,商建,杨汝良.一种有效的二维MUSIC超分辨SAR成像算法[J].系统仿真学报,2010,22(1):184-187”。有的是需要多个快拍，即多个距离向回波信号进行联合处理，且只计算出目标的位置，没有计算出信号的回波幅度，详见“张东浩.线阵三维SAR成像算法研究及仿真[D].电子科技大学,2010”。有的是计算出目标的位置及幅度，但是需要多个快拍才能完成计算，详见“成晨.三维合成孔径雷达超分辨成像方法研究[D].电子科技大学,2012”。
发明内容
[0007] 本发明提出了一种基于单快拍MUSIC算法的距离向处理方法，它主要运用单快拍MUSIC算法估计目标的位置，找到目标位置以后，再通过最小二乘法估计目标的幅度。通过这两步处理，就完成了SAR成像过程中的距离向处理。
[0008] 为了方便描述本发明的内容，首先作以下术语定义：
[0009] 定义1、合成孔径雷达(SAR)
[0010] 合成孔径雷达是将雷达是将雷达固定于载荷运动平台上，结合运动平台的运动以合成线性阵列以达到运动向的分辨率，再利用雷达波束向回波延时实现距离一维成像，从而实现观测目标二维成像的一种合成孔径雷达技术。
[0011] 定义2、标准的MUSIC算法
[0012] 标准的MUSIC算法是一种基于矩阵特征空间分解的方法。信号处理的观测空间可以分解为信号子空间和噪声子空间，显然这两个空间是正交的。信号子空间由阵列接收到的数据协方差矩阵中与信号对应的特征向量组成，噪声子空间则由协方差矩阵中所有最小特征值(噪声方差)对应的特征向量组成。通过MUSIC算法得到的频谱叫做MUSIC谱。详见“何子述.现代数字信号处理及其应用[M].清华大学出版社,2009”。
[0013] 定义3、标准合成孔径雷达原始回波仿真方法
[0014] 标准合成孔径雷达原始回波仿真方法是指给定雷达系统参数、平台轨迹参数以及观测场景参数等所需的参数条件下，基于合成孔径雷达成像原理仿真得到具有SAR回波信号特性的原始回波信号的方法，详细内容可参考文献“Lan G.Cumming Frank,H Wong.合成孔径雷达成像:算法与实现[M].电子工业出版社,2012”。
[0015] 将垂直于雷达平台运动的方向叫做距离向。对距离向进行离散化处理，离散化之后的距离向由多个小格组成，每个小格被称为距离单元。
[0016] 定义4、单快拍
[0017] 快拍数一般是指时域上的采样点数。单快拍是指全部阵元采样一次。协方差矩阵受到快拍数影响，如果快拍数少，则协方差矩阵估计不准确，那么噪声子空间不准，最终将导致MUSIC算法性能下降。
[0018] 定义5、标准的特征分解方法
[0019] 特征分解(Eigen decomposition)又被称为谱分解(Spectral decomposition)。它是将矩阵分解为由其特征值和特征向量表示的矩阵之积的方法。详见“矩阵理论”，黄廷祝编，高等教育出版社。
[0020] 定义6、标准的去斜处理方法
[0021] 去斜处理就是将输入的线性调频信号与一个调频率大小相同但符号相反的参考信号相乘，去掉2次相位，得到频率与信号输入位置成正比的单频信号。由于单频信号的频率与输入信号位置相关，可通过FFT变换实现对线性调频信号的压缩。
[0022] 定义7、最小二乘法
[0023] 最小二乘法是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。详见“Sorenson H W.Least-squares estimation:from Gauss to Kalman[J].IEEE spectrum,1970,7(7):63-68”。
[0024] 定义8、标准的施密特正交化方法
[0025] 施密特正交化(Schmidt orthogonalization)是求欧氏空间正交基的一种方法。在线性代数中，如果内积空间上的一组向量能够组成一个子空间，那么这一组向量就称为这个子空间的一个基。施密特正交化提供了一种方法，能够通过这一子空间上的一个基得出子空间的一个正交基，并可进一步求出对应的标准正交基。详见“线性代数”，同济大学数学系编，高等教育出版社。
[0026] 定义9、共轭转置
[0027] 共轭转置即对复数矩阵A进行转置并取共轭，记为AH，可通过标准的共轭转置方法计算得出，。详见“线性代数”，同济大学数学系编，高等教育出版社。
[0028] 定义10、标准的矩阵求逆方法
[0029] 假设矩阵A和矩阵B，如果AB＝E，其中E为单位矩阵，则称矩阵B是矩阵A的右逆矩阵，通常将矩阵B写为A-1，根据矩阵A，通过标准的矩阵求逆方法计算出矩阵A-1。详见“矩阵理论”，黄廷祝，钟守铭等，高等教育出版社。
[0030] 定义11、共轭复数
[0031] 共轭复数，两个实部相等，虚部互为相反数的复数互为共轭复数。详见“线性代数”，同济大学数学系编，高等教育出版社。
[0032] 定义12、标准的谱峰搜索方法
[0033] 谱峰搜索即找出信号频谱中谱峰的位置，可以将谱峰视为极大值点，在谱峰搜索的过程中，通常将设置一个门限值，小于门限值的频谱会被忽略。MATLAB中已有谱峰搜索函数，即’findpeaks’，可以直接调用，实现对MUSIC谱的谱峰搜索，本文中的谱峰搜索方法即基于MATLAB中的’findpeaks’函数。详见MATLAB中的帮助文档。
[0034] 定义13、标准合成孔径雷达原始回波仿真方法
[0035] 标准合成孔径雷达原始回波仿真方法是指给定雷达系统参数、平台轨迹参数以及观测场景参数等所需的参数条件下，基于合成孔径雷达成像原理仿真得到具有SAR回波信号特性的原始回波信号的方法，详细内容可参考文献“Lan G.Cumming Frank,H Wong.合成孔径雷达成像:算法与实现[M].电子工业出版社,2012”。
[0036] 本发明提供的一种基于单快拍MUSIC算法的距离向处理方法，它包括以下步骤：
[0037] 步骤1、初始化SAR系统参数：
[0038] 初始化SAR系统参数包括：雷达工作中心频率，记做fc；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做Br；雷达发射信号脉冲宽度，记做Tr；雷达发射信号的调频斜率，记做Kr；雷达接收系统的采样频率，记做fs；电磁波在空气中的传播速度，记做C；连续采样时间为TP，n为距离向时间序号，n＝1,2,…,Nr，n为自然数，Nr为距离向采样点总数；上述参数均为SAR系统标准参数，其中雷达中心频率fc，雷达载频波长λ，雷达发射基带信号的信号带宽Br，雷达发射信号脉冲宽度Tr，雷达发射信号调频斜率Kr，雷达接收系统的采样频率fs，波门延时为τ0，在SAR系统设计过程中已经确定；根据SAR成像系统方案和观测方案，距离向处理所需要的初始化成像系统参数均为已知；
[0039] 步骤2、初始化SAR回波信号：
[0040] 在第t个距离向采样得到的原始回波数据，记为SaE(t)，t＝1,2,…,TP，t为自然数；在实际线阵SAR系统中，原始回波数据SaE(t)可由雷达系统数据接收机提供；在仿真线阵SAR成像过程中，原始回波数据SaE(t)根据雷达系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在雷达回波数据进行距离向处理之前，原始回波数据SaE(t)均已知；
[0041] 初始化场景大小为[-X,X]，将原点设置场景中心，定义场景中有M个散射点，且目标都在距离向接收范围内，这M个散射点相对于场景中心位置为x1,x2,...,xM，这M个点的散射系数分别记为 通过公式 得到M个散射点所对应的回波时延，分别为其中，x为散射点相对于场景中心的距离，C为光速，原始回波数据SaE(t)通
过采用标准合成孔径雷达原始回波仿真方法产生得到，即
[0042]
[0043] 其中，A为回波信号的增益，这里将A设为1，根据采样频率fs对连续回波信号SaE(t)进行采样，得到离散回波信号SAE(n)，即
[0044]
[0045] 其中，Nr为距离向采样点数；
[0046] 采用标准的去斜处理方法，得到参考信号
[0047]
[0048] 其中，n＝1,2,…,Nr，n为自然数；上式中所使用的参数均为步骤1中初始化的参数；
[0049] 步骤3、对原始回波数据进行去斜处理：
[0050] 通过标准的取共轭方法对参考信号href(n)进行处理，得到共轭参考信号
[0051] 通过将步骤2中的回波信号SAE(n)与共轭参考信号 相乘，实现去除回波信号相位中二次项，得到处理后的信号SE(n)，这个过程用公式表示为
[0052]
[0053] 其中，n＝1,2,…,Nr，n为自然数， 为M个点的散射系数，为M个散射点所对应的回波时延分，Kr为雷达发射信号的调频斜率，fs为雷达
接收系统的采样频率，fc为雷达工作中心频率，SAE(n)为步骤2中得到的离散回波信号；
[0054] 步骤4、生成噪声子空间矩阵：
[0055] 步骤4.1、构造自相关矩阵：
[0056] 步骤4.1.1初始化矩阵及信号向量：
[0057] 初始化一个Nr×Nr维零矩阵Ax，矩阵Ax的结构如下：
[0058]
[0059] 矩阵Ax中的元素ax,y均为0，其中，x，y为自然数，x＝1,2,…,Nr，y＝1,2,…,Nr；
[0060] 回波信号SE(n)表示为向量形式SE＝[SE(1),SE(2),…,SE(Nr)]，其中，SE(n)为步骤3中得到的n时刻的回波信号，n＝1,2,…,Nr；
[0061] 步骤4.1.2、对矩阵Ax第m行处理的第1步：
[0062] 选择向量SE＝[SE(1),SE(2),…,SE(Nr)]中第1个到第Nr-m+1个元素，这Nr-m+1个元素定义为向量SE(m,1)＝[SE(1),SE(2),…,SE(Nr-m+1)]，然后将向量SE(m,1)填入矩阵Ax中的第m行，第m列到第Nr列的元素位置，得到矩阵Ax(m,1)，其中，Ax是步骤4.1.1中定义的矩阵，矩阵Ax(m,1)的结果如下
[0063]
[0064] 其中，SE(n)为步骤4.1.1中向量的元素SE＝[SE(1),SE(2),…,SE(Nr)]，n＝1,2,…,Nr-m+1；ax,y为矩阵Ax中的元素，其中，x，y为自然数，x＝1,2,…,Nr，y＝1,2,…,Nr；
[0065] 步骤4.1.3、对矩阵Ax第m行处理的第2步：
[0066] 选择向量SE＝[SE(1),SE(2),…,SE(Nr)]中第2个到第m个元素，这m-1个元素定义为向量SE(m,2)＝[SE(2),SE(3),…,SE(m)]，然后将向量SE(m,2)中元素进行反向排序得到向量SE(m,-2)＝[SE(m),SE(m-1),…,SE(2)]，采用标准的取共轭的方法对向量SE(m,-2)进行处理得到共轭向量SE*(m,-2)＝[SE*(m),SE*(m-1),…,SE*(2)]，最后，将向量SE*(m,-2)填入矩阵Ax(m,1)的第m行，第1列到第m列的位置，得到矩阵Ax(m,2)，其中，Ax(m,1)是步骤4.1.2中的矩阵，矩阵Ax(m,2)的结果如下
[0067]
[0068] 其中，SE(n)为向量SE＝[SE(1),SE(2),…,SE(Nr)]的元素，n＝1,2,…,Nr-m+1；SE*(n)为向量SE*(m,-2)＝[SE*(m),SE*(m-1),…,SE*(2)]中的元素，其中n＝2,3,…,m；ax,y为矩阵Ax中的元素，其中，x，y为自然数，x＝1,2,…,Nr，y＝1,2,…,Nr；
[0069] 步骤4.1.4、完成矩阵构造及计算自相关矩阵：
[0070] 通过步骤4.1.2和步骤4.1.3完成对矩阵Ax第m行处理，得到矩阵Ax(m,2)，其中，Ax是步骤4.1.1中定义的矩阵，Ax(m,2)是步骤4.1.3中定义的矩阵；
[0071] 按照步骤4.1.2和步骤4.1.3的方法对矩阵Ax进行处理，详细过程为：Ax作为初始矩阵，首先对矩阵Ax的第一行进行处理，得到矩阵Ax(1,2)，此时只执行步骤4.1.2，然后对矩阵Ax(1,2)的第二行进行处理，即执行步骤4.1.2和步骤4.1.3，得到Ax(2,2)，依次进行下去，最后对矩阵 的第Nr行进行处理，得到 矩阵 的结构如下：
[0072]
[0073] 将矩阵 命名为协方差矩阵Rx，采用标准的共轭转置方法得到共轭转置矩阵通过公式 得到自相关矩阵R；
[0074] 步骤4.2、生成噪声子空间矩阵：
[0075] 步骤4.2.1、对矩阵R进行特征分解：
[0076] 根据步骤4.1.4中得到的自相关矩阵R，使用标准的特征分解方法对其进行处理，得到特征值集合D'和对应的特征向量集合V'，且特征值集合D'与特征向量集合V'是一一对应的；
[0077] 步骤4.2.2、提取特征向量：
[0078] 对特征值集合D '中的元素从大到小进行排序，得到特征值集合其中， 为步骤4.2.1中得到的特征值集合D'中的元素，且
提取特征值集合D中最后Nr-K个元素，得到特征值集合
其中dn为特征值集合D，n为自然数，n＝K+1,K+2,...,Nr；
[0079] 从特征向量集合V'中提取出Nr-K个与特征值集合D”对应的特征向量集合其中， 为特征向量集合V '中的元素，且
与 一一对应；
[0080] 采用标准的施密特正交化方法对特征向量集合 进行处理，得到特征向量集合 其中un为使用标准的施密特正交化方法处理后
的特征向量，n为自然数，n＝K+1,K+2,...,Nr；
[0081] 步骤4.2.3、生成噪声子空间矩阵
[0082] 根据标准的MUSIC算法和步骤4.2.2中得到的特征向量集合通过公式 得到噪声子空间矩阵矩阵
[0083] 步骤5、计算MUSIC谱：
[0084] 连续采样时间为TP，采样点数设置为常数Nr，连续采样时间TP除以采样点数Nr得到采样间隔dτ；从τ0时刻开始到TP+τ0时刻结束，每过dτ时长记录一次当前时刻的值，将该时刻的值记为τn，得到时刻集合 其中，n＝1,2,...,Nr，τ0为波门延时，再将时刻集合 按照下标从小到大排序，得到离散时间向量τ
[0085] τ0≤τn≤TP+τ0
[0086] τn＝τ0+(n-1)·dτ
[0087] 其中，τn为时间轴第n个采样点的时刻，n为自然数，n＝1,2,...,Nr，τ0为波门延时，即信号通过接收机系统自身的延时，波门延时τ0为已知量；
[0088] 然后通过下面的公式得到MUSIC谱其中τn为离散时间轴，n为自然数，n＝1,2,...,Nr；
[0089]
[0090] 其中，
[0091]
[0092]
[0093] 其中，Kr为信号调频斜率，fS为雷达接收系统的采样频率，τk为离散时间，e为自然对数的底数，j表示虚部，-jnωm表示系数为n，频率为ωm的相位，其中n＝1,2,...,Nr，m＝1,2,......,K；即，-jω1表示系数为1、频率为ω1的相位，-jω2表示系数为1、频率为ω2的相位，-jωK表示系数为1、频率为ωK的相位；-j(Nr-1)ω1表示系数为Nr-1、频率为ω1的相位，-j(Nr-1)ω2表示系数为Nr-1、频率为ω2的相位，-j(Nr-1)ωK表示系数为Nr-1、频率为ωK的相位；
[0094] 步骤6、谱峰搜索：
[0095] 步骤6.1、归一化处理：
[0096] 使用标准的归一化方法对MUSIC谱进行归一化处理，得到谱
[0097] 步骤6.2、设置谱峰搜索门限：
[0098] 对谱PMUSIC'(τ)中的元素PMUSIC'(τn)从大到小依次进行排序，其中n＝1,2,...,Nr，得到向量 PMUSIC'(τn')为PMUSIC'(τ)中的元素，其中n＝1,2,...,Nr，且 其中n＝1,2,...,
Nr，根据目标个数K得到门限θ＝PMUSIC'(τM')，其中M为自然数，且K
[0099] 步骤6.3、计算目标位置：
[0100] 根据步骤6.2中得到的门限θ和目标个数K，采用MATLAB中标准的谱峰搜索方法对PMUSIC(τ)1进行处理，得到PMUSIC(τ)1中最大的K个谱峰所对应的时延量τk；
[0101] 然后通过公式 得到延时量τk所对应的目标所在单元格位置IDk，其中，τ0为波门延时，dτ为采样间隔；
[0102] 步骤7、估计目标所在距离单元的值：
[0103] 步骤7.1、计算相位：
[0104] 通过公式 得到最小二乘估计矩阵中第n行，第k列元素的相位 其中，n＝1,2,...,Nr，τk为延时量，k＝1,2,…,K，Kr为信号调
频斜率，fS为雷达接收系统的采样频率；
[0105] 步骤7.2、构造最小二乘估计矩阵：
[0106] 通过公式
[0107]
[0108] 得到最小二乘估计矩阵
[0109]
[0110] 其中， 为步骤7.1中得到的相位，n＝1,2,...,Nr，τk为延时量，k＝1,2,…,K，Kr为信号调频斜率，fS为雷达接收系统的采样频率，e为自然对数的底数，j表示虚部； 表示第n个采样点，延时为τk的相位，其中n＝1,2,...,Nr，m＝1,2,......,K；
即， 表示第1个采样点，延时为τ1的相位， 表示第1个采样点，延时为τ2的相位，表示第1个采样点，延时为τK的相位； 表示第2个采样点，延时为τ1的相位， 表示第2个采样点，延时为τ2的相位， 表示第2个采样点，延时为τK的相位； 表示第Nr个采样点，延时为τ1的相位， 表示第Nr个采样点，延时为τ2的相位， 表示第Nr个采样点，延时为τK的相位；
[0111] 步骤7.3、估计目标所在距离单元的值：
[0112] 采用标准的共轭转置方法得到最小二乘估计矩阵L1的最小二乘估计的共轭转置矩阵
[0113] 通过最小二乘估计矩阵L1左乘最小二乘估计的共轭转置矩阵 得到矩阵
[0114] 然后采用标准的矩阵求逆方法得到矩阵γ0的逆矩阵
[0115] 通过公式 得到IDk处的值，其中，L1为最小二乘估计矩阵， 为最小二乘估计矩阵的共轭转置矩阵，SE为步骤3中去斜处理的回波数据向量SE(n)，IDk为步骤6中计算出的目标所在单元格的位置；
[0116] 至此，已完成目标位置和幅度的估计，完成了距离向处理。
[0117] 本发明方法的主要思路是：首先，对原始回波信号进行去斜处理，然后，根据单快拍MUSIC算法进行构造自相关矩阵、特征分解等过程。然后，计算MUSIC谱并进行谱峰搜索。这样就能确定目标所在的距离单元格。最后，使用最小二乘法计算目标所在距离单元格的值。
[0118] 本发明的优点在于利用单快拍MUSIC算法进行目标所在距离单元格的估计，通过最小二乘法估计目标所在距离单元的值，本发明与通过脉冲压缩方法进行距离向处理相比，本发明仅需要单个回波信号即可使用MUSIC算法进行距离向的处理，分辨率有较大的提高，能够实现距离向超分辨能力，本发明适用于目标较少，精度要求高的情况。
附图说明
[0119] 图1为本发明流程图；
具体实施方式
[0120] 本发明主要采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2017a上验证正确。具体实施步骤如下：
[0121] 步骤1、初始化SAR系统参数：
[0122] 初始化SAR系统参数包括：雷达工作中心频率，记做fc＝10GHz；雷达载频波长，记做λ＝0.03m；雷达发射基带信号的信号带宽，记做Br＝100MHz；雷达发射信号脉冲宽度，记做Tr＝2.56μs；雷达发射信号的调频斜率，记做Kr＝3.91×1013Hz/s；雷达接收系统的采样频率，记做fs＝400MHz；电磁波在空气中的传播速度，记做C＝3×108m/s；连续采样时间为TP，n为距离向时间序号，n＝1,2,…,Nr，n为自然数，Nr＝1024为距离向采样点总数；上述参数均为SAR系统标准参数，其中雷达中心频率fc＝10GHz，雷达载频波长λ＝0.03m，雷达发射基带信号的信号带宽Br＝100MHz，雷达发射信号脉冲宽度Tr＝2.56μs，雷达发射信号调频斜率Kr＝3.91×1013Hz/s，雷达接收系统的采样频率fs＝400MHz，τ0为波门延时，在SAR系统设计过程中已经确定；根据SAR成像系统方案和观测方案，距离向处理所需要的初始化成像系统参数均为已知；
[0123] 步骤2、初始化SAR回波信号：
[0124] 在第t个距离向采样得到的原始回波数据，记为SaE(t)，t＝1,2,…,TP，t为自然数；在实际线阵SAR系统中，原始回波数据SaE(t)可由雷达系统数据接收机提供；在仿真线阵SAR成像过程中，原始回波数据SaE(t)根据雷达系统参数，采用标准合成孔径雷达原始回波仿真方法产生得到；在雷达回波数据进行距离向处理之前，原始回波数据SaE(t)均已知；
[0125] 初始化场景大小为[-19.2,19.2]，将原点设置场景中心，定义场景中有3散射点，且目标都在距离向接收范围内，这3个散射点相对场景中心位置为x1＝-4,x2＝0,x3＝1.5，这3个点的散射系数分别记为 通过公式 得到3个散射点所对应的回波时延，分别为 其中，x为散射点相对于场景中心的距离，C为光速，原始回波数据SaE(t)通过采用标准合成孔径雷达原始回波仿真方法产生得到，即
[0126]
[0127] 其中，A为回波信号的增益，这里将A设为1，根据采样频率fs＝400MHz对连续回波信号SaE(t)进行采样，得到离散回波信号SAE(n)，即
[0128]
[0129] 根据标准的去斜处理方法，得到参考信号
[0130]
[0131] 其中，n＝1,2,…,1024，n为自然数；上式中所使用的参数均为步骤1中初始化的参数；
[0132] 步骤3、对原始回波数据进行去斜处理：
[0133] 通过标准的取共轭方法对参考信号href(n)进行处理，得到共轭参考信号通过将步骤2中的回波信号SAE(n)与共轭参考信号 相乘，实现去除回波信号相位中二次项，得到处理后的信号SE(n)，这个过程用公式表示为
[0134]
[0135] 其中，n＝1,2,…,1024， 为3个点的散射系数， 为3个散射点所13
对应的回波时延分，Kr＝3.91×10 Hz/s为雷达发射信号的调频斜率，fs＝400MHz为雷达接收系统的采样频率，fc＝10GHz为雷达工作中心频率，SAE(n)为步骤2中得到的离散回波信号；
[0136] 步骤4、生成噪声子空间矩阵：
[0137] 步骤4.1、构造自相关矩阵：
[0138] 步骤4.1.1初始化矩阵及信号向量：
[0139] 初始化一个1024×1024维零矩阵Ax，矩阵Ax的结构如下：
[0140]
[0141] 矩阵Ax中的元素ax,y均为0，其中，x，y为自然数，x＝1,2,…,1024，y＝1,2,…,1024；
[0142] 回波信号SE(n)表示为向量形式SE＝[SE(1),SE(2),…,SE(1024)]，其中，SE(n)为步骤3中得到的n时刻的回波信号，n＝1,2,…,1024；
[0143] 步骤4.1.2、对矩阵Ax第m行处理的第1步：
[0144] 选择向量SE＝[SE(1),SE(2),…,SE(1024)]中第1个到第1025-m个元素，这1025-m个元素定义为向量SE(m,1)＝[SE(1),SE(2),…,SE(1025-m)]，然后将向量SE(m,1)填入矩阵Ax中的第m行，第m列到第1025列的元素位置，得到矩阵Ax(m,1)，其中，Ax是步骤4.1.1中定义的矩阵，矩阵Ax(m,1)的结果如下
[0145]
[0146] 其中，SE(n)为步骤4.1.1中向量的元素SE＝[SE(1),SE(2),…,SE(1024)]，n＝1,2,…,1025-m；ax,y为矩阵Ax中的元素，其中，x，y为自然数，x＝1,2,…,1024，y＝1,2,…,
1024；
[0147] 步骤4.1.3、对矩阵Ax第m行处理的第1步：
[0148] 选择向量SE＝[SE(1),SE(2),…,SE(1024)]中第2个到第m个元素，这m-1个元素定义为向量SE(m,2)＝[SE(2),SE(3),…,SE(m)]，然后将向量SE(m,2)中元素进行反向排序得到向量SE(m,-2)＝[SE(m),SE(m-1),…,SE(2)]，采用标准的取共轭的方法对向量SE(m,-2)进行处理得到共轭向量SE*(m,-2)＝[SE*(m),SE*(m-1),…,SE*(2)]，最后，将向量SE*(m,-2)填入矩阵Ax(m,1)的第m行，第1列到第m列的位置，得到矩阵Ax(m,2)，其中，Ax(m,1)是步骤4.1.2中的矩阵，矩阵Ax(m,2)的结果如下
[0149]
[0150] 其中，SE(n)为向量SE＝[SE(1),SE(2),…,SE(1024)]的元素，n＝1,2,…,1025-m；SE*(n)为向量SE*(m,-2)＝[SE*(m),SE*(m-1),…,SE*(2)]中的元素，其中n＝2,3,…,m；ax,y为矩阵Ax中的元素，其中，x，y为自然数，x＝1,2,…,1024，y＝1,2,…,1024；
[0151] 步骤4.1.4、完成矩阵构造及计算自相关矩阵：
[0152] 通过步骤4.1.2和步骤4.1.3完成对矩阵Ax第m行处理，得到矩阵Ax(m,2)，其中，Ax是(m,2)步骤4.1.1中定义的矩阵，Ax 是步骤4.1.3中定义的矩阵；
[0153] 按照步骤4.1.2和步骤4.1.3的方法对矩阵Ax进行处理，详细过程为：Ax作为初始矩阵，首先对矩阵Ax的第一行进行处理，得到矩阵Ax(1,2)，此时只执行步骤4.1.2，然后对矩阵Ax(1,2)的第二行进行处理，即执行步骤4.1.2和步骤4.1.3，得到Ax(2,2)，依次进行下去，最后对(1023,2) (1024,2) (1024,2)矩阵Ax 的第Nr行进行处理，得到Ax ，；矩阵Ax 的结构如下：
[0154]
[0155] 将矩阵Ax(1024,2)命名为协方差矩阵Rx，采用标准的共轭转置方法得到共轭转置矩阵 通过公式 得到自相关矩阵R；
[0156] 步骤4.2、生成噪声子空间矩阵：
[0157] 步骤4.2.1、对矩阵R进行特征分解：
[0158] 根据步骤4.1.4中得到的自相关矩阵R，使用标准的特征分解方法对其进行处理，得到特征值集合D'和对应的特征向量集合V'，且特征值集合D'与特征向量集合V'是一一对应的；
[0159] 步骤4.2.2、提取特征向量：
[0160] 对特征值集合D'中的元素从大到小进行排序，得到特征值集合D＝{d1,d2,......,d1024}，其中，d1,d2,......,d1024为步骤4.2.1中得到的特征值集合D'中的元素，且d1>d2>......>d1024；提取特征值集合D中最后1021个元素，得到特征值集合D”＝{d4,d5,......,d1024}，其中dn为特征值集合D，n为自然数，n＝4,5,...,1024；
[0161] 从特征向量集合V'中提取出1021个与特征值集合D”对应的特征向量集合V”＝{v4,v5,......,v1024}，其中，v4,v5,......,v1024为特征向量集合V'中的元素，且d4,d5,......,d1024与v4,v5,......,v1024一一对应；
[0162] 采用标准的施密特正交化方法对特征向量集合V”＝{v4,v5,......,v1024}进行处理，得到特征向量集合U＝u4,u5,…,u1024，其中un为使用标准的施密特正交化方法处理后的特征向量，n为自然数，n＝4,5,...,1024；
[0163] 步骤4.2.3、生成噪声子空间矩阵
[0164] 根据标准的MUSIC算法和步骤4.2.2中得到的特征向量集合U＝u4,u5,…,u1024，通过公式 得到噪声子空间矩阵矩阵
[0165] 步骤5、计算MUSIC谱：
[0166] 连续采样时间为TP，采样点数设置为1024，连续采样时间TP除以采样点数1024得到采样间隔dτ，从τ0时刻开始到Tp+τ0时刻结束，每过dτ时长记录一次当前时刻的值，将该时刻的值记为τn，得到时刻集合{τ1,τ2,......,τ1024}，其中τ0为波门延时，再将时刻集合{τ1,τ2,......,τ1024}按照下标从小到大排序，得到离散时间向量τ
[0167] τ＝[τ1 τ2 … τ1024]，τ0≤τn≤Tp+τ0
[0168] τn＝τ0+(n-1)·dτ
[0169] 其中，τn为时间轴第n个采样点的时刻，n为自然数，n＝1,2,...,1024，τ0为波门延时，即信号通过接收机系统自身的延时，波门延时τ0为已知量；
[0170] 然后通过下面的公式得到MUSIC谱PMUSIC(τ)＝[PMUSIC(τ1),PMUSIC(τ2),......,PMUSIC(τ1024)]，其中τn为离散时间轴，n为自然数，n＝1,2,...,1024；
[0171]
[0172] 其中，
[0173]
[0174]
[0175] 其中，Kr＝3.91×1013Hz/s为信号调频斜率，fs＝400MHz为雷达接收系统的采样频率，τk为离散时间，e为自然对数的底数，j表示虚部，-jnωm表示系数为n，频率为ωm的相位，其中n＝1,2,...,1024，m＝1,2,3；即，-jω1表示系数为1、频率为ω1的相位，-jω2表示系数为1、频率为ω2的相位，-jω3表示系数为1、频率为ω3的相位；-j1023ω1表示系数为1023、频率为ω1的相位，-j1023ω2表示系数为1023、频率为ω2的相位，-j1023ω3表示系数为1023、频率为ω3的相位；
[0176] 步骤6、谱峰搜索：
[0177] 步骤6.1、归一化处理：
[0178] 使用标准的归一化方法对MUSIC谱PMUSIC(τ)＝[PMUSIC(τ1),PMUSIC(τ2),......,PMUSIC(τ1024)]进行归一化处理，得到谱PMUSIC'(τ)＝[PMUSIC'(τ1),PMUSIC'(τ2),......,PMUSIC'(τ1024)]；
[0179] 步骤6.2、设置谱峰搜索门限：
[0180] 对谱PMUSIC'(τ)中的元素PMUSIC'(τn)从大到小依次进行排序，其中n＝1,2,...,1024，得到向量PMUSIC'(τ')＝[PMUSIC'(τ1'),PMUSIC'(τ2'),......,PMUSIC'(τ1024')]，PMUSIC'(τn')为PMUSIC'(τ)中的元素，其中n＝1,2,...,1024，且PMUSIC'(τ1')>PMUSIC'(τ2')>......>PMUSIC'(τ1024')，其中n＝1,2,...,1024，根据目标个数K＝3得到门限θ＝PMUSIC'(τM')，其中M为自然数，且3
[0181] 步骤6.3、计算目标位置：
[0182] 根据步骤6.2中得到的门限θ和目标个数K＝3，采用MATLAB中标准的谱峰搜索方法对PMUSIC(τ)1进行处理，得到PMUSIC(τ)1中最大的3个谱峰所对应的时延量τk，k＝1,2,3；
[0183] 然后通过公式 得到延时量τk所对应的目标所在单元格位置IDk，其中，τ0为波门延时，dτ为采样间隔；
[0184] 步骤7、估计目标所在距离单元的值：
[0185] 步骤7.1、计算相位：
[0186] 通过公式 得到最小二乘估计矩阵中第n行，第k列元素的相位 其中，Kr＝3.91×1013Hz/s为信号调频斜率，fs＝400MHz为雷
达接收系统的采样频率，τk为离散时间，k＝1,2,3；
[0187] 步骤7.2、构造最小二乘估计矩阵：
[0188] 通过公式
[0189]
[0190] 得到最小二乘估计矩阵
[0191]
[0192] 其中， 为步骤7.1中得到的相位，n＝1,2,...,1024，τk为延时量，k＝1,2,3，Kr＝3.91×1013Hz/s为信号调频斜率，fs＝400MHz为雷达接收系统的采样频率，e为自然对数的底数，j表示虚部； 表示第n个采样点，延时为τk的相位，其中n＝1,
2,...,1024，m＝1,2,3；即， 表示第1个采样点，延时为τ1的相位， 表示第1个采样点，延时为τ2的相位， 表示第1个采样点，延时为τ3的相位； 表示第2个采样点，延时为τ1的相位， 表示第2个采样点，延时为τ2的相位， 表示第2个采样点，延时为τ3的相位；
表示第1024个采样点，延时为τ1的相位， 表示第1024个采样点，延时为τ2的相位， 表示第1024个采样点，延时为τ3的相位；
[0193] 步骤7.3、估计目标所在距离单元的值：
[0194] 采用标准的共轭转置方法得到最小二乘估计矩阵L1的最小二乘估计的共轭转置矩阵
[0195] 通过最小二乘估计矩阵L1左乘最小二乘估计的共轭转置矩阵 得到矩阵
[0196] 然后采用标准的矩阵求逆方法得到矩阵γ0的逆矩阵
[0197] 通过公式 得到ID1，ID2，ID3处的值γ1,γ2,γ3，其中，L1为最小二乘估计矩阵， 为最小二乘估计矩阵的共轭转置矩阵，SE为步骤3中去斜处理的回波数据向量SE(n)，IDk为步骤6中计算出的目标所在单元格的位置；
[0198] 至此，已完成目标位置和幅度的估计，完成了距离向处理。
[0199] 估计出三个目标所在位置分别为-3.544米、0.019米和1.481米；归一化幅度分别为0dB，-24.9dB和-8.6dB，噪声水平为-67.5dB；至此，我们得到了目标的位置信息和幅度信息，整个方法结束。
[0200] 经过计算机仿真结果证明，本方法有效地估计出距离向上目标的位置信息以及幅度信息。
