CN107783078B - 一种波束-多普勒酉esprit多目标角度估计方法 - Google Patents
一种波束-多普勒酉esprit多目标角度估计方法 Download PDFInfo
- Publication number
- CN107783078B CN107783078B CN201710813630.4A CN201710813630A CN107783078B CN 107783078 B CN107783078 B CN 107783078B CN 201710813630 A CN201710813630 A CN 201710813630A CN 107783078 B CN107783078 B CN 107783078B
- Authority
- CN
- China
- Prior art keywords
- matrix
- domain
- time
- target
- doppler
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 41
- 239000011159 matrix material Substances 0.000 claims abstract description 143
- 238000009499 grossing Methods 0.000 claims abstract description 27
- 230000009466 transformation Effects 0.000 claims description 18
- 238000000354 decomposition reaction Methods 0.000 claims description 17
- 230000001427 coherent effect Effects 0.000 claims description 8
- 238000001514 detection method Methods 0.000 claims description 7
- 238000003491 array Methods 0.000 claims description 5
- 239000007983 Tris buffer Substances 0.000 claims description 3
- 238000011156 evaluation Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims 1
- 238000004422 calculation algorithm Methods 0.000 abstract description 39
- 238000004364 calculation method Methods 0.000 abstract description 5
- 238000005516 engineering process Methods 0.000 description 6
- 238000009825 accumulation Methods 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 230000002123 temporal effect Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 239000000243 solution Substances 0.000 description 3
- 238000007476 Maximum Likelihood Methods 0.000 description 2
- 238000000342 Monte Carlo simulation Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000001131 transforming effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明介绍了一种波一种波束‑多普勒酉ESPRIT多目标角度估计方法，主要为解决多目标的波达方向(DOA)估计问题。所述方法的一具体实施方式包括：首先对目标所在距离单元数据进行时域平滑，得到平滑后的L个时域快拍数据，然后采用中心共轭对称傅里叶变换矩阵将时域快拍数据变换至波束‑多普勒域，同时保留旋转不变结构，最后采用实值酉ESPRIT算法估计目标的DOA。该实施方式充分利用了信号的时域信息来改善空域参数估计性能，精度高且计算复杂度低。
Description
技术领域
本发明属于雷达技术领域，涉及基于时域平滑技术的波束-多普勒域二维酉ESPRIT算法。
背景技术
波束内多目标的分辨是地基和/或机载预警雷达对主波束内多个目标进行精细跟踪时需要解决的关键问题。多目标编队飞行时很容易落入同一距离单元或相邻距离单元，仅依靠脉冲压缩技术很难从距离上对这些目标进行分离。同时，由于预警雷达的相干积累时间通常较短，仅从多普勒维对径向速度接近的目标进行分离也很困难。本发明考虑从空间维通过高精度的角度测量来分辨主瓣内的多个目标。
经典测角方法主要包括基于和差波束的比幅测角和多波束比幅测角等方法。但是，当波束内存在多个目标时，传统测角方法将失效。实际上，可采用空域超分辨技术来实现波束内多目标的分辨。经典的空域超分辨算法主要包括ML、MUSIC、ESPRIT及相应的快速算法。这些超分辨算法只利用了一维空域信息，时域脉冲通常当做快拍来使用，当信号的波达方向夹角很小时，角度估计性能会变差，尤其是在低信噪比、少快拍数时，其性能下降非常严重。实际上对于预警雷达来说，目标信号的时域信息(例如各目标多普勒频率的差异)是可以利用的。在波达方向估计中，利用信号的时域信息，本质上是将空域一维参数估计问题转化为空时二维参数估计问题，一般情况下，信源间的空间距离会增大。因此，利用信号的时域结构信息，有可能改善空域参数的估计性能，特别是小夹角时的性能。
在空时二维数据结构下，该问题本质上是单快拍多目标超分辨DOA(direction ofarrival，高分辨波达方向)估计问题。可以采用最大似然方法来联合估计目标的DOA和多普勒频率，但该方法需要进行二维参数搜索，运算量很大。通过降维最大似然方法降低运算量，将DOA和多普勒频率联合估计问题解耦成两个一维顺序优化问题，但是该方法仍然需要一维参数搜索，并且容易收敛到局部最优解。实际上，基于子空间的空域超分辨算法也可以扩展应用于空时二维参数估计，例如二维MUSIC算法，但该算法也需要二维参数搜索。针对通信系统中多载频信号DOA估计问题，目前也有一些解决的方法，如Lemma AN等人提出基于ESPRIT的角度-频率联合估计(JAFE)方法，相对于其他子空间方法，该方法避免了二维参数搜索，运算量得以降低；而基于实数运算的酉-JAFE(U-JAFE)算法，进一步降低了运算量。但传统JAFE和U-JAFE算法需要在阵元-脉冲域估计信号的协方差矩阵并且进行特征分解，运算量仍然较大(特别当脉冲数和阵元数较多时)。近年来，稀疏信号恢复技术已被广泛应用于阵列信号处理领域。虽然该技术可扩展应用于空时二维参数估计，并且在低快拍数时具有较好的参数估计性能，但需要求解一个高维的优化问题，计算量也比较大。
发明内容
本发明的目的在于提高DOA估计精度同时兼顾运算量，在单快拍空时二维数据模型下，提出了一种基于时域平滑技术的波束-多普勒域二维酉ESPRIT算法。
实现本发明目的的技术方案是：首先通过时域平滑技术构造多个快拍数据。然采用中心共轭对称傅里叶变换矩阵将快拍数据变换至波束-多普勒域，同时保留旋转不变结构，最后采用实值酉ESPRIT算法估计目标的DOA。
本申请提出的一种波束-多普勒酉ESPRIT多目标角度估计方法，该方法包括：步骤1对目标所在距离单元数据进行时域平滑，得到平滑后的L个时域快拍数据，其中，第l个时域快拍数据可表示为xl，l＝1,...,L；步骤2利用目标检测阶段提供的角度和多普勒频率先验信息构造指向所述目标的空时二维波束域变换矩阵；步骤3利用所述空时二维波束域变换矩阵将平滑后的时域快拍数据xl变换到波束-多普勒域，得到空时二维波束域快拍数据yl，其中，l＝1,...,L；步骤4根据所述空时二维波束域快拍数据yl，构造实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，并对所述实值数据矩阵Y进行奇异值分解，得到P个大奇异值对应的信号子空间E和实值矩阵Φμ、Φν；步骤5对复值矩阵Φμ+jΦν进行特征分解，并估计信号的DOA和多普勒频率。
在一些实施例中，步骤1对目标所在距离单元数据进行时域平滑，得到平滑后的L个时域快拍数据，其中，第l个时域快拍数据可表示为xl，包括：假设时域子孔径的长度为M(我们将M称为时域平滑因子)、相干累计脉冲数为K，且M<K，那么通过前向时域平滑可获得L＝K-M+1个时域快拍数据，其中，第l(l＝1,...,L)个时域快拍数据可以表示为：
xl＝Jlx0
式中，x0为目标所在距离单元的空时快拍信号，
为第l个空时二维子孔径选择矩阵，其中，IN表示N×N维的单位矩阵，N表示接收端为由N个列子阵合成的等效均匀线阵，Jt,l＝[0M×(l-1),IM,0M×(K-l-M+1)]为M×K维的时域子孔径选择矩阵，表示选择第l个时域子孔径，t表示时域，0M×(l-1)表示M×(l-1)维的零元素矩阵，该矩阵中的元素全为零，IM表示M×M维的单位矩阵，0M×(K-l-M+1)表示M×(K-l-M+1)维的零元素矩阵，
表示直积。
在一些实施例中，步骤2利用目标检测阶段提供的角度和多普勒频率先验信息构造指向所述目标的空时二维波束域变换矩阵，包括如下步骤：
3a)假设波束域变换矩阵为Ws＝[ws,0,ws,1,...,ws,N-1]，则Ws的第n(n＝0,1,...,N-1)列为：
表示将第n个空域波束指向空间频率μ＝n(2π/N)，其中，符号T表示矩阵的转置，s表示空域。
3b)假设多普勒域变换矩阵为W t＝[wt,0,wt,1,...,wt,M-1]，则Wt的第m(m＝0,1,...,M-1)列为：
表示将第m个时域波束指向多普勒频率ν＝m(2π/M)。
3c)空时二维波束域变换矩阵可表示为：
在一些实施例中，步骤3利用所述空时二维波束域变换矩阵将平滑后的时域快拍数据xl变换到波束-多普勒域，得到空时二维波束域快拍数据yl，包括如下步骤：
4a)计算空时二维波束域MN×P维的导向矢量矩阵B：
式中，WH为所述空时二维波束域变换矩阵W的共轭转置，H为共轭转置符号，A＝[a1,a2,...,aP]为KN×P维的导向矢量矩阵，P为假设有P个相互统计独立的目标位于主波束且处于同一距离单元，
为第p(p＝1,...,P)个目标的空时二维导向矢量，其中，
为第p个目标的空域导向矢量，s表示空域，
为第p个目标的时域导向矢量，且μp＝2πdsin(θp)/λ，νp＝2πfpTr，其中，p＝1,...,P，d为列子阵间距，λ为波长，Tr为脉冲重复周期，θp和fp分别表示第p(p＝1,...,P)个目标的入射角和多普勒频率；
为MN×P维的子孔径导向矢量矩阵，
为第p个目标的空时二维子孔径导向矢量，其中，
为第p个目标的时域子孔径导向矢量。
式中，Bs＝[bs(μ1),bs(μ2),...,bs(μP)]为波束域导向矢量矩阵，其中，bs(μp)＝[bs,0(μp),bs,1(μp),...,bs,N-1(μp)]T(p＝1,...,P)，表示第p个目标的波束域导向矢量，bs(μp)的第n个元素即为第n个波束对应的输出，可表示为：
式中，Bt＝[bt(ν1),bt(ν2),...,bt(νP)]为多普勒域导向矢量矩阵，其中，bt(νp)＝[bt,0(νp),bt,1(νp),...,bt,M-1(νp)]T(p＝1,...,P)表示第p个目标的多普勒域导向矢量，bt(νp)第m个元素即为第m个波束对应的输出，可表示为：
4b)将经时域平滑后的第l个时域快拍数据xl变换至波束-多普勒域，可得空时二维波束域快拍数据yl：
式中，
为经时域平滑后第l个快拍对应的目标复幅度矢量，其中，diag{·}为向量对角矩阵化函数，s＝[s1,...,sP]T为P个目标经距离匹配滤波后的复幅度向量，sp为第p(p＝1,...,P)个目标的复幅度，n l＝Jln为第l个子孔径对应的噪声矢量，
为第l个子孔径对二维波束域的噪声矢量，n为KN×1维的零均值复高斯白噪声矢量，其协方差矩阵为σ2IKN，σ2表示白噪声功率，IKN表示KN×KN维的单位矩阵。
在一些实施例中，步骤4根据所述空时二维波束域快拍数据yl构造实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，并对所述实值数据矩阵Y进行奇异值分解，得到P个大奇异值对应的信号子空间E和实值矩阵Φμ、Φν，包括如下步骤：
5a)求取“旋转不变性”方程：
由上式可知，bs,n+1(μp)的分子与bs,n(μp)的分子互为相反数，bs,n(μp)相当于bs(μp)＝[bs,0(μp),bs,1(μp),...,bs,N-1(μp)]T式中第n项，bs,n+1(μp)相当于该式中第n+1项，由此可知bs(μp)的前N-1个相邻分量满足如下关系：
对上式进行三角函数变换可得：
由于第n＝0号波束与第n＝N-1号波束对应的空间频率分别为μp,0＝0和μp,N-1-2π＝(N-1)/(2π/N)-2π＝-2π/N，因此这两个波束在物理空间上是相邻的。bs(μp)的最后一个元素bs,N-1(μp)和第一个元素bs,0(μp)之间的关系为：
结合上述可以得到关于bs(μp)的N个“旋转不变性”方程，即：
式中，Γ1和Γ2为选择矩阵，其表达式如下：
5b)计算空时二维波束域的旋转不变性关系：
对于P个目标来说，波束空间导向矢量矩阵具有如下旋转不变性关系：
Γ1BsΩμ＝Γ2Bs
式中，Ωμ＝diag{[tan(μ1/2),tan(μ2/2),...,tan(μP/2)]T}为实值对角矩阵，其对角元素包含了P个目标的DOA信息；
同样的，根据波束域旋转不变性，计算得到多普勒域的旋转不变性关系：
Γ3BtΩν＝Γ4Bt
式中，Γ3、Γ4为选择矩阵，Γ3、Γ4与Γ1、Γ2的定义类似(需要将N替换为M)，Ων＝diag{[tan(ν1/2),tan(ν2/2),...,tan(νP/2)]T}为实值对角矩阵，其对角元素包含了P个目标的多普勒频率信息。
结合上述和Kronecker积的性质可得空时二维波束域的旋转不变性关系为：
其中，
和
为波束域选择矩阵，
和
为多普勒域选择矩阵。
5c)由波束-多普勒域数据(空时二维波束域快拍数据yl)可构造出一实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，其中，Re[]表示实部，Im[]表示虚部。
5d)对所述实值矩阵Y进行奇异值分解可得到实值信号子空间E，理论上该信号子空间可由波束-多普勒域导向矢量矩阵B的各列张成信号子空间，即：
E＝BT
式中，T为P×P维的实值非奇异矩阵。
5e)将B＝ET-1代入式
可得：
式中，Φμ＝T-1ΩμT，Φν＝T-1ΩνT，T-1是T的逆矩阵。
在一些实施例中，步骤5对复值矩阵Φμ+jΦν进行特征分解，并估计信号的DOA和多普勒频率，包括如下步骤：
6a)通过对矩阵Φμ+jΦν进行特征分解来实现空间频率μ和多普勒频率ν的估计和自动配对，即：
Φμ+jΦν＝T-1(Ωμ+jΩν)T
从特征值矩阵Ωμ+jΩν的实部和虚部可估计得到tan(μp/2)和tan(νp/2)(p＝1,...,P)的值，其中，p＝1,...,P。
根据估值得到的tan(μp/2)和tan(νp/2)(p＝1,...,P)，求得第p(p＝1,...,P)个目标空间频率μp和多普勒频率νp；
6b)第p个目标的入射角
和多普勒频率
可由下式估计：
式中，λp为Φμ+jΦν的第p个特征值，d为列子阵间距，λ为波长，Tr为脉冲重复周期。
本发明与现有技术相比具有以下优点：
1.本发明采用时域平滑技术获得多个空时二维快拍数据，然后利用预警雷达目标检测阶段提供的空时二维参数粗分辨信息(角度和多普勒频率先验信息)，形成少量指向目标潜在区域的空时二维波束。这样的处理方法使得参数估计可在低维空间进行，相对于现有的解决方法来说，显著地降低了运算复杂度。
2.在形成空时二维波束时，本发明采用了中心共轭对称傅里叶变换矩阵，可将数据变换到波束-多普勒域的同时保留旋转不变结构，然后采用实值酉ESPRIT算法估计目标的DOA估计方法。时域信息的DOA估计方法的引入可能会带来运算量的少量增加，但与此同时能够显著提高DOA估计精度。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显：
图1为根据本申请的一种波束-多普勒酉ESPRIT多目标角度估计方法的一个实施例的流程图；
图2是各种算法的DOA估计性能随信噪比变化的曲线；
图3是各种算法的运算复杂度对比(其中，K＝2N，K表示相干积累脉冲数，N表示接收端为由N个列子阵合成的等效均匀线阵)。
具体实施方式
下面结合附图和具体实施方案对本发明的技术方案作进一步详细地说明。可以理解的是，此处所描述的具体实施例仅仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
参考图1，示出了根据本申请的一种波束-多普勒酉ESPRIT多目标角度估计方法的一个实施例的流程图100。所述的波束-多普勒酉ESPRIT多目标角度估计方法，包括以下步骤：
步骤101，对目标所在距离单元数据进行时域平滑，得到平滑后的L个时域快拍数据，其中，第l个时域快拍数据可表示为xl，l＝1,...,L。
假设时域子孔径的长度为M(我们将M称为时域平滑因子)、相干累计脉冲数为K，且M<K，那么通过前向时域平滑可获得L＝K-M+1个时域快拍数据，其中，第l(l＝1,...,L)个时域快拍数据可以表示为：
xl＝Jlx0
式中，x0为目标所在距离单元的空时快拍信号，
为第l个空时二维子孔径选择矩阵，其中，IN表示N×N维的单位矩阵，N表示接收端为由N个列子阵合成的等效均匀线阵，Jt,l＝[0M×(l-1),IM,0M×(K-l-M+1)]为M×K维的时域子孔径选择矩阵，表示选择第l个时域子孔径，t表示时域，0M×(l-1)表示M×(l-1)维的零元素矩阵，该矩阵中的元素全为零，IM表示M×M维的单位矩阵，0M×(K-l-M+1)表示M×(K-l-M+1)维的零元素矩阵，
表示直积。
步骤102，利用目标检测阶段提供的角度和多普勒频率先验信息构造指向目标的空时二维波束域变换矩阵。
目标检测通常在PD(Pulse Doppler，经脉冲多普勒)处理后的距离-多普勒平面进行。常规PD处理的角度和多普勒频率分辨率都不能突破瑞利限，但能够提供目标大概位置的先验信息。因此，我们可以形成少量空时二维波束来覆盖目标所在角度和多普勒频带，这样目标的参数估计可以在低维的空时二维波束域内进行，运算复杂度将显著降低。这里我们假设目标位于第n个空域波束(对应Ws的第n列)和m个多普勒通道(对应Wt的第m列)。
该步骤可分解成如下几步：
3a)假设波束域变换矩阵为W s＝[ws,0,ws,1,...,ws,N-1]，则Ws的第n(n＝0,1,...,N-1)列为：
表示将第n个空域波束指向空间频率μ＝n(2π/N)，其中，符号T表示矩阵的转置，s表示空域。
3b)假设多普勒域变换矩阵为W t＝[wt,0,wt,1,...,wt,M-1]，则Wt的第m(m＝0,1,...,M-1)列为：
表示将第m个时域波束指向多普勒频率ν＝m(2π/M)。
3c)空时二维波束域变换矩阵可表示为：
步骤103，利用空时二维波束域变换矩阵将平滑后的时域快拍数据xl变换到波束-多普勒域，得到空时二维波束域快拍数据yl，其中，l＝1,...,L。
该步骤可分解成如下几步：
4a)计算空时二维波束域MN×P维的导向矢量矩阵B：
式中，WH为所述空时二维波束域变换矩阵W的共轭转置，H为共轭转置符号，A＝[a1,a2,...,aP]为KN×P维的导向矢量矩阵，P为假设有P个相互统计独立的目标位于主波束且处于同一距离单元，
为第p(p＝1,...,P)个目标的空时二维导向矢量，其中，
为第p个目标的空域导向矢量，
为第p个目标的时域导向矢量，且μp＝2πdsin(θp)/λ，νp＝2πfpTr，其中，p＝1,...,P，d为列子阵间距，λ为波长，Tr为脉冲重复周期，θp和fp分别表示第p(p＝1,...,P)个目标的入射角和多普勒频率；
为MN×P维的子孔径导向矢量矩阵，
为第p个目标的空时二维子孔径导向矢量，其中，
为第p个目标的时域子孔径导向矢量。
式中，Bs＝[bs(μ1),bs(μ2),...,bs(μP)]为波束域导向矢量矩阵，其中，bs(μp)＝[bs,0(μp),bs,1(μp),...,bs,N-1(μp)]T(p＝1,...,P)，表示第p个目标的波束域导向矢量，bs(μp)的第n个元素即为第n个波束对应的输出，可表示为：
式中，Bt＝[bt(ν1),bt(ν2),...,bt(νP)]为多普勒域导向矢量矩阵，其中，bt(νp)＝[bt,0(νp),bt,1(νp),...,bt,M-1(νp)]T(p＝1,...,P)表示第p个目标的多普勒域导向矢量，bt(νp)第m个元素即为第m个波束对应的输出，可表示为：
4b)将经时域平滑后的第l个时域快拍数据xl变换至波束-多普勒域，可得空时二维波束域快拍数据yl：
式中，
为经时域平滑后第l个快拍对应的目标复幅度矢量，其中，diag{·}为向量对角矩阵化函数，s＝[s1,...,sP]T为P个目标经距离匹配滤波后的复幅度向量，sp为第p(p＝1,...,P)个目标的复幅度，n l＝Jln为第l个子孔径对应的噪声矢量，
为第l个子孔径对二维波束域的噪声矢量，n为KN×1维的零均值复高斯白噪声矢量，其协方差矩阵为σ2IKN，σ2表示白噪声功率，IKN表示KN×KN维的单位矩阵。
步骤104，根据空时二维波束域快拍数据yl构造实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，并对实值数据矩阵Y进行奇异值分解，得到P个大奇异值对应的信号子空间E和实值矩阵Φμ、Φν。
该步骤可分解成如下几步：
5a)求取“旋转不变性”方程：
由上式可知，bs,n+1(μp)的分子与bs,n(μp)的分子互为相反数，bs,n(μp)相当于bs(μp)＝[bs,0(μp),bs,1(μp),...,bs,N-1(μp)]T式中第n项，bs,n+1(μp)相当于该式中第n+1项，由此可知bs(μp)的前N-1个相邻分量满足如下关系：
对上式进行三角函数变换可得：
由于第n＝0号波束与第n＝N-1号波束对应的空间频率分别为μp,0＝0和μp,N-1-2π＝(N-1)/(2π/N)-2π＝-2π/N，因此这两个波束在物理空间上是相邻的。bs(μp)的最后一个元素bs,N-1(μp)和第一个元素bs,0(μp)之间的关系为：
结合上述可以得到关于bs(μp)的N个“旋转不变性”方程，即：
式中，Γ1和Γ2为选择矩阵，其表达式如下：
5b)计算空时二维波束域的旋转不变性关系：
对于P个目标来说，波束空间导向矢量矩阵具有如下旋转不变性关系：
Γ1BsΩμ＝Γ2Bs
式中，Ωμ＝diag{[tan(μ1/2),tan(μ2/2),...,tan(μP/2)]T}为实值对角矩阵，其对角元素包含了P个目标的DOA信息。
同样的，根据波束域旋转不变性，计算得到多普勒域的旋转不变性关系：
Γ3BtΩν＝Γ4Bt
式中，Γ3、Γ4为选择矩阵，Γ3、Γ4与Γ1、Γ2的定义类似(需要将N替换为M)，Ων＝diag{[tan(ν1/2),tan(ν2/2),...,tan(νP/2)]T}为实值对角矩阵，其对角元素包含了P个目标的多普勒频率信息。
式中，Γ3和Γ4的表达式如下：
结合上述和Kronecker积的性质可得空时二维波束域的旋转不变性关系为：
其中，
和
为波束域选择矩阵，
和
为多普勒域选择矩阵。
5c)由波束-多普勒域数据(空时二维波束域快拍数据yl)可构造出一实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，其中，Re[]表示实部，Im[]表示虚部。
5d)对所述实值矩阵Y进行奇异值分解可得到实值信号子空间E，理论上该信号子空间可由波束-多普勒域导向矢量矩阵B的各列张成信号子空间，即：
E＝BT
式中，T为P×P维的实值非奇异矩阵；
5e)将B＝ET-1代入式
可得：
式中，Φμ＝T-1ΩμT，Φν＝T-1ΩνT，T-1是T的逆矩阵。
步骤105，对复值矩阵Φμ+jΦν进行特征分解，并估计信号的DOA和多普勒频率。
该步骤可分解成如下几步：
6a)通过对矩阵Φμ+jΦν进行特征分解来实现空间频率μ和多普勒频率ν的估计和自动配对，即：
Φμ+jΦν＝T-1(Ωμ+jΩν)T
从特征值矩阵Ωμ+jΩν的实部和虚部可估计得到tan(μp/2)和tan(νp/2)(p＝1,...,P)的值，其中，p＝1,...,P。
根据估值得到的tan(μp/2)和tan(νp/2)(p＝1,...,P)，求得第p(p＝1,...,P)个目标空间频率μp和多普勒频率νp。
6b)第p个目标的入射角
和多普勒频率
可由下式估计：
式中，λp为Φμ+jΦν的第p个特征值，d为列子阵间距，λ为波长，Tr为脉冲重复周期。
本发明的优点可通过仿真数据实验进一步说明。
1.仿真参数
在仿真实验中，我们考虑采用8个阵元的均匀线阵，阵元间距为d＝λ/2。脉冲重复频率Tr为1000Hz，相干积累脉冲数K为32个。假设两个等功率的窄带目标信号同时位于主波束内，这两个目标的角度和多普勒频率分别为25.9°，34.1°和463Hz，537Hz，其角度和多普勒间隔均为瑞利限的2/3。蒙特卡罗实验次数为Nm＝500，DOA估计的均方根误差作为衡量角度估计性能的参数，其定义为
其中，P表示相互统计独立目标的个数，θp为第p(p＝1,...,P)个目标的入射角，
表示第l次蒙特卡罗实验对第p个目标的DOA估计值。
2.仿真数据处理结果及分析
①与类ESPRIT算法的性能比较
选取B-UESPRIT和U-JAFE类ESPRIT算法与所提算法进行DOA估计性能对比。本发明所提算法简称为BD-UESPRIT(beam-Doppler unitary ESPRIT)，相应的时域平滑因子为0.4K，K为相干积累脉冲数，多普勒通道和波束数量分别为Kt＝3和Ks＝4。B-UESPRIT算法仅在空间波束域进行参数估计，时域脉冲当作快拍来用，相应的波束数量为4个。U-JAFE算法中的时域平滑因子与所提BD-UESPRIT算法相同。
图2给出了上述算法的DOA估计误差(RMSE of DOA)随信噪比(SNR)变化的关系以及CRB(Cramer-Rao bound，克拉美罗界)。从图2中可以看到，所提BD-UESPRIT算法的DOA估计精度要明显优于B-UESPRIT算法，因为所提算法充分利用了目标的时域信息。另外，所提算法在低信噪比时的估计精度要高于U-JAFE算法，在高信噪比时的性能相当。这是由于波束变换后目标能量被集中在少数几个波束内，信噪比得到了增强，因此所提算法的DOA估计精度要优于在阵元域估计参数的U-JAFE算法。
②与类ESPRIT算法的运算量比较
目标数量P＝3，波束和多普勒通道数量分别为Kt＝3和Ks＝4，时域平滑因子设为M＝0.4K。图3给出了U-JAFE，B-UESPRIT和BD-UESPRIT算法所消耗的实乘次数随系统自由度的变化关系，其中，K＝2N，K表示相干积累脉冲数，N表示接收端为由N个列子阵合成的等效均匀线阵。从图3中可以看出U-JAFE算法的运算量随着系统自由度的增加而急剧增加，而B-UESPRIT和BD-UESPRIT算法的运算量增加的速度很缓慢。这是由于B-UESPRIT和BD-UESPRIT算法的运算复杂度主要取决于空间波束和多普勒通道数，它们不随自由度的变化而变化。此外，由于本文算法利用了额外的时域自由度，因此所提BD-UESPRIT算法的运算复杂度要稍高于B-UESPRIT算法，但这是可以接受的，因为BD-UESPRIT算法的DOA估计精度得到了明显的提升。
2.2实验结论：
本发明提出了一种波束-多普勒域酉ESPRIT高精度DOA估计方法。该方法首先通过时域平滑获取多个快拍，然后利用中心共轭对称DFT矩阵将数据变换至波束-多普勒域，最后利用目标在波束-多普勒域的旋转不变性关系来求解目标的DOA和多普勒频率。仿真结果表明，相对于没有利用时域信息的DOA估计方法，本文方法在运算量少量增加的情况下能够显著提高DOA估计精度。
以上描述仅为本申请的较佳实施例以及对所运用技术原理的说明。本领域技术人员应当理解，本申请中所涉及的发明范围，并不限于上述技术特征的特定组合而成的技术方案，同时也应涵盖在不脱离所述发明构思的情况下，由上述技术特征或其等同特征进行任意组合而形成的其它技术方案。例如上述特征与本申请中公开的(但不限于)具有类似功能的技术特征进行互相替换而形成的技术方案。
Claims (6)
1.一种波束-多普勒酉ESPRIT多目标角度估计方法，其特征在于，所述方法包括如下步骤：
步骤1对目标所在距离单元数据进行时域平滑，得到平滑后的L个时域快拍数据，其中，第l个时域快拍数据可表示为xl，l＝1,...,L；
步骤2利用目标检测阶段提供的角度和多普勒频率先验信息构造指向所述目标的空时二维波束域变换矩阵；
步骤3利用所述空时二维波束域变换矩阵将平滑后的时域快拍数据xl变换到波束-多普勒域，得到空时二维波束域快拍数据yl，其中，l＝1,...,L；
步骤4根据所述空时二维波束域快拍数据yl，构造实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，并对所述实值数据矩阵Y进行奇异值分解，得到P个大奇异值对应的信号子空间E和实值矩阵Φμ、Φν；
步骤5对复值矩阵Φμ+jΦν进行特征分解，并估计信号的DOA和多普勒频率。
2.根据权利要求1所述的波束-多普勒酉ESPRIT多目标角度估计方法，其特征在于，步骤1所述的对目标所在距离单元数据进行时域平滑，得到平滑后的L个时域快拍数据，其中，第l个时域快拍数据可表示为xl，包括：
假设时域子孔径的长度为M、相干累计脉冲数为K，且M<K，那么通过前向时域平滑可获得L＝K-M+1个时域快拍数据，其中，第l(l＝1,...,L)个时域快拍数据可以表示为：
xl＝Jlx0
式中，x0为目标所在距离单元的空时快拍信号，
为第l个空时二维子孔径选择矩阵，其中，IN表示N×N维的单位矩阵，N表示接收端为由N个列子阵合成的等效均匀线阵，Jt,l＝[0M×(l-1),IM,0M×(K-l-M+1)]为M×K维的时域子孔径选择矩阵，表示选择第l个时域子孔径，t表示时域，0M×(l-1)表示M×(l-1)维的零元素矩阵，该矩阵中的元素全为零，IM表示M×M维的单位矩阵，0M×(K-l-M+1)表示M×(K-l-M+1)维的零元素矩阵，
表示直积。
3.根据权利要求2所述的波束-多普勒酉ESPRIT多目标角度估计方法，其特征在于，步骤2所述的利用目标检测阶段提供的角度和多普勒频率先验信息构造指向所述目标的空时二维波束域变换矩阵，包括如下步骤：
3a假设波束域变换矩阵为Ws＝[ws,0,ws,1,...,ws,N-1]，则Ws的第n(n＝0,1,...,N-1)列为：
表示将第n个空域波束指向空间频率μ＝n(2π/N)，其中，符号T表示矩阵的转置，s表示空域；
3b假设多普勒域变换矩阵为Wt＝[wt,0,wt,1,...,wt,M-1]，则Wt的第m(m＝0,1,...,M-1)列为：
表示将第m个时域波束指向多普勒频率ν＝m(2π/M)；
3c空时二维波束域变换矩阵可表示为：
4.根据权利要求3所述的波束-多普勒酉ESPRIT多目标角度估计方法，其特征在于，步骤3所述的利用所述空时二维波束域变换矩阵将平滑后的时域快拍数据xl变换到波束-多普勒域，得到空时二维波束域快拍数据yl，包括如下步骤：
4a计算空时二维波束域MN×P维的导向矢量矩阵B：
式中，WH为所述空时二维波束域变换矩阵W的共轭转置，H为共轭转置符号，A＝[a1,a2,...,aP]为KN×P维的导向矢量矩阵，P为假设有P个相互统计独立的目标位于主波束且处于同一距离单元，
为第p(p＝1,...,P)个目标的空时二维导向矢量，其中，
为第p(p＝1,...,P)个目标的空域导向矢量，
为第p(p＝1,...,P)个目标的时域导向矢量，且μp＝2πdsin(θp)/λ，νp＝2πfpTr，其中，p＝1,...,P，d为列子阵间距，λ为波长，Tr为脉冲重复周期，θp和fp分别表示第p(p＝1,...,P)个目标的入射角和多普勒频率；
为MN×P维的子孔径导向矢量矩阵，
为第p(p＝1,...,P)个目标的空时二维子孔径导向矢量，其中，
为第p(p＝1,...,P)个目标的时域子孔径导向矢量；
式中，Bs＝[bs(μ1),bs(μ2),...,bs(μP)]为波束域导向矢量矩阵，其中，bs(μp)＝[bs,0(μp),bs,1(μp),...,bs,N-1(μp)]T，表示第p(p＝1,...,P)个目标的波束域导向矢量，bs(μp)的第n个元素即为第n个波束对应的输出，可表示为：
式中，Bt＝[bt(ν1),bt(ν2),...,bt(νP)]为多普勒域导向矢量矩阵，其中，bt(νp)＝[bt,0(νp),bt,1(νp),...,bt,M-1(νp)]T表示第p(p＝1,...,P)个目标的多普勒域导向矢量，bt(νp)第m个元素即为第m个波束对应的输出，可表示为：
4b将经时域平滑后的第l个时域快拍数据xl变换至波束-多普勒域，可得空时二维波束域快拍数据yl：
式中，
为经时域平滑后第l个快拍对应的目标复幅度矢量，其中，diag{·}为向量对角矩阵化函数，s＝[s1,...,sP]T为P个目标经距离匹配滤波后的复幅度向量，sp(p＝1,...,P)为第p(p＝1,...,P)个目标的复幅度，nl＝Jln为第l个子孔径对应的噪声矢量，
为第l个子孔径对二维波束域的噪声矢量，n为KN×1维的零均值复高斯白噪声矢量，其协方差矩阵为σ2IKN，σ2表示白噪声功率，IKN表示KN×KN维的单位矩阵。
5.根据权利要求4所述的波束-多普勒酉ESPRIT多目标角度估计方法，其特征在于，步骤4所述的根据所述空时二维波束域快拍数据yl构造实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，并对所述实值数据矩阵Y进行奇异值分解，得到P个大奇异值对应的信号子空间E和实值矩阵Φμ、Φν，包括如下步骤：
5a求取“旋转不变性”方程：
由上式可知，bs,n+1(μp)的分子与bs,n(μp)的分子互为相反数，bs,n(μp)相当于bs(μp)＝[bs,0(μp),bs,1(μp),...,bs,N-1(μp)]T式中第n项，bs,n+1(μp)相当于该式中第n+1项，由此可知bs(μp)的前N-1个相邻分量满足如下关系：
对上式进行三角函数变换可得：
由于第n＝0号波束与第n＝N-1号波束对应的空间频率分别为μp,0＝0和μp,N-1-2π＝(N-1)/(2π/N)-2π＝-2π/N，因此这两个波束在物理空间上是相邻的，bs(μp)的最后一个元素bs,N-1(μp)和第一个元素bs,0(μp)之间的关系为：
结合上述可以得到关于bs(μp)的N个“旋转不变性”方程，即：
式中，Γ1和Γ2为选择矩阵，其表达式如下：
5b计算空时二维波束域的旋转不变性关系：
对于P个目标来说，波束空间导向矢量矩阵具有如下旋转不变性关系：
Γ1BsΩμ＝Γ2Bs
式中，Ωμ＝diag{[tan(μ1/2),tan(μ2/2),...,tan(μP/2)]T}为实值对角矩阵，其对角元素包含了P个目标的DOA信息；
同样的，根据波束域旋转不变性，计算得到多普勒域的旋转不变性关系：
Γ3BtΩν＝Γ4Bt
式中，Γ3、Γ4为选择矩阵，Γ3、Γ4与Γ1、Γ2的定义类似，即将N替换为M，Ων＝diag{[tan(ν1/2),tan(ν2/2),...,tan(νP/2)]T}为实值对角矩阵，其对角元素包含了P个目标的多普勒频率信息；
结合上述和Kronecker积的性质可得空时二维波束域的旋转不变性关系为：
其中，
和
为波束域选择矩阵，
和
为多普勒域选择矩阵；
5c由波束-多普勒域数据即空时二维波束域快拍数据yl，构造出一实值数据矩阵Y＝[Re[y1,…,yL],Im[y1,…,yL]]，其中，Re[]表示实部，Im[]表示虚部；
5d对所述实值矩阵Y进行奇异值分解可得到实值信号子空间E，理论上该信号子空间可由波束-多普勒域导向矢量矩阵B的各列张成信号子空间，即：
E＝BT
式中，T为P×P维的实值非奇异矩阵；
5e将B＝ET-1代入式
可得：
式中，Φμ＝T-1ΩμT，Φν＝T-1ΩνT，T-1是T的逆矩阵。
6.根据权利要求5所述的波束-多普勒酉ESPRIT多目标角度估计方法，其特征在于，步骤5所述的对复值矩阵Φμ+jΦν进行特征分解，并估计信号的DOA和多普勒频率，包括如下步骤：
6a通过对矩阵Φμ+jΦν进行特征分解来实现空间频率μ和多普勒频率ν的估计和自动配对，即：
Φμ+jΦν＝T-1(Ωμ+jΩν)T
从特征值矩阵Ωμ+jΩν的实部和虚部可估计得到tan(μp/2)和tan(νp/2)的值，其中，p＝1,...,P；
根据估值得到的tan(μp/2)和tan(νp/2)(p＝1,...,P)，求得第p(p＝1,...,P)个目标空间频率μp和多普勒频率νp；
6b第p个目标的入射角
和多普勒频率
可由下式估计：
式中，λp为Φμ+jΦν的第p个特征值，d为列子阵间距，λ为波长，Tr为脉冲重复周期。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710813630.4A CN107783078B (zh) | 2017-09-11 | 2017-09-11 | 一种波束-多普勒酉esprit多目标角度估计方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710813630.4A CN107783078B (zh) | 2017-09-11 | 2017-09-11 | 一种波束-多普勒酉esprit多目标角度估计方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107783078A CN107783078A (zh) | 2018-03-09 |
| CN107783078B true CN107783078B (zh) | 2020-09-22 |
Family
ID=61437880
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710813630.4A Expired - Fee Related CN107783078B (zh) | 2017-09-11 | 2017-09-11 | 一种波束-多普勒酉esprit多目标角度估计方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107783078B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109521393A (zh) * | 2018-11-05 | 2019-03-26 | 昆明理工大学 | 一种基于信号子空间旋转特性的波达方向估计算法 |
| WO2020097903A1 (zh) * | 2018-11-16 | 2020-05-22 | 华为技术有限公司 | 测角方法以及雷达设备 |
| CN109655804B (zh) * | 2019-01-24 | 2022-11-22 | 南京航空航天大学 | 一种基于奇异值分解的临近目标相对距离估计方法 |
| CN110161492B (zh) * | 2019-01-24 | 2020-12-08 | 北京机电工程研究所 | 舰船航向航速提取方法 |
| TWI693419B (zh) * | 2019-02-13 | 2020-05-11 | 國立交通大學 | 訊號處理方法 |
| CN109946664B (zh) * | 2019-03-06 | 2023-03-24 | 西安电子科技大学 | 一种主瓣干扰下的阵列雷达导引头单脉冲测角方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2001281326A (ja) * | 2000-03-29 | 2001-10-10 | Toyota Central Res & Dev Lab Inc | レーダ信号処理回路 |
| CN104865556A (zh) * | 2015-05-18 | 2015-08-26 | 哈尔滨工程大学 | 基于实域加权最小化l1范数方法的MIMO雷达系统DOA估计方法 |
-
2017
- 2017-09-11 CN CN201710813630.4A patent/CN107783078B/zh not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2001281326A (ja) * | 2000-03-29 | 2001-10-10 | Toyota Central Res & Dev Lab Inc | レーダ信号処理回路 |
| CN104865556A (zh) * | 2015-05-18 | 2015-08-26 | 哈尔滨工程大学 | 基于实域加权最小化l1范数方法的MIMO雷达系统DOA估计方法 |
Non-Patent Citations (2)
| Title |
|---|
| "Beamspace ESPRIT algorithm for bistatic MIMO radar";Guo Y D et al.;《Electronic Letters》;20121231;第47卷(第15期);全文 * |
| "非理想系统假设下阵列信号多信源多参数联合估计";吴湘霖;《中国优秀博硕士学位论文全文数据库（博士） 信息科技辑》;20070415;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107783078A (zh) | 2018-03-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107783078B (zh) | 一种波束-多普勒酉esprit多目标角度估计方法 | |
| CN103383452B (zh) | 分布式阵列目标到达角估计方法 | |
| CN110031794B (zh) | 一种基于差分共性阵重构的相干信源doa估计方法 | |
| CN105403856B (zh) | 基于嵌套式最小冗余阵列的波达方向估计方法 | |
| US20230152424A1 (en) | Method for eliminating one-bit signal harmonic false target and related component | |
| CN108872970B (zh) | 适用于一般等间距稀疏阵列单频信号波束形成的栅瓣判别方法 | |
| CN107505602A (zh) | 嵌套阵下基于dft的doa估计方法 | |
| CN103018730A (zh) | 分布式子阵波达方向估计方法 | |
| CN110346752B (zh) | 基于互质稀疏阵的无模糊测向方法 | |
| CN107092004A (zh) | 基于信号子空间旋转不变性的互质阵列波达方向估计方法 | |
| Fuchs et al. | Single-snapshot direction-of-arrival estimation of multiple targets using a multi-layer perceptron | |
| CN104730513A (zh) | 一种分级子阵聚焦mvdr波束形成方法 | |
| CN103605107B (zh) | 基于多基线分布式阵列的波达方向估计方法 | |
| CN104898119A (zh) | 一种基于相关函数的动目标参数估计方法 | |
| CN103364762B (zh) | 任意阵列流形的单基地mimo雷达波达方向估计方法 | |
| CN108398659B (zh) | 一种矩阵束与求根music结合的波达方向估计方法 | |
| Feng et al. | MIMO–monopulse target localisation for automotive radar | |
| CN109946672B (zh) | 基于被动孔径合成稀疏阵列的doa估计方法 | |
| CN115436896A (zh) | 快速的雷达单快拍music测角方法 | |
| CN109491009B (zh) | 一种光纤组合阵及基于光纤组合阵的栅瓣抑制方法 | |
| CN117420539A (zh) | 稀疏阵毫米波雷达频域波束降维快速联合超分辨估计方法 | |
| CN114609580A (zh) | 一种基于非圆信号的无孔互质阵列设计方法 | |
| CN112666558B (zh) | 一种适用于汽车fmcw雷达的低复杂度music测向方法及装置 | |
| Yang et al. | Motion-guided large aperture ULA to enhance DOA estimation: An inverse synthetic aperture perspective | |
| CN111458703A (zh) | 一种测量多目标横向速度的方法及系统 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200922 Termination date: 20210911 |