CN107656237B

CN107656237B - 一种多信源频率与doa联合检测的方法及其装置

Info

Publication number: CN107656237B
Application number: CN201710657595.1A
Authority: CN
Inventors: 黄翔东; 刘明卓; 李长滨
Original assignee: Tianjin University
Current assignee: Tianjin University
Priority date: 2017-08-03
Filing date: 2017-08-03
Publication date: 2020-12-01
Anticipated expiration: 2037-08-03
Also published as: CN107656237A

Abstract

本发明公开了一种多信源频率与DOA联合检测的方法及其装置，包括：对L路信号样本做M_f点DFT，利用Tsui频谱校正器对DFT结果进行频率和相位校正，得到校正后的D组频率、相位、幅值的参数组，利用参数组构造出D个方向矢量；将D个方向矢量按照最小距离进行分组匹配，得出D个信源参数匹配信息；根据D个信源参数匹配信息构造频率余数组，并将频率余数组带入闭式鲁棒中国余数定理模型进行重构，得出频率估计值；根据D个信源参数匹配信息得出L‑1个相位差，并构造相位余数组，将相位余数组与重构模值组带入闭式鲁棒CRT重构得到中间参数，计算出DOA估计值。本发明实现了多对目标的高精度频率与DOA联合估计。

Description

一种多信源频率与DOA联合检测的方法及其装置

技术领域

本发明涉及信号分析处理技术领域，尤其涉及在观测多信源数入射到阵列天线的情况下，基于一种新型互素稀疏阵列，对其进行欠采样处理，并进行频谱校正与中国余数定理重构，从而实现对多目标频率与DOA(到达角)的高精度联合测量。

背景技术

多目标信源的频率与DOA联合估计是阵列信号处理领域的研究热点，在雷达^[1]、无线通信^[2]、电子战等领域有着广泛应用。目前现有的空间谱估计算法多是基于子空间分解谱估计算法(如MUSIC(多重信号分类算法)、ESPRIT(旋转不变信号参数估计技术)算法等)，然而这类算法主要存在两个方面的不足。一是需要满足奈奎斯特采样速率，即阵元传感器的采样速率需要大于等于最大入射信号频率的2倍。这也同时导致了所采用的均匀线性阵列的阵元间距不能大于信号波长的一半。因此，在高频段下，阵元采样器的功率难以满足实际需求；二是，必须将阵元密集排列，导致阵元间产生较大耦合，从而对信源参数的估计精度产生较大影响。

针对多目标信号的频率与DOA联合估计问题，目前研究者提出了一些新的研究方法：文献[3]提出了一种基于均匀线阵，利用长短基线概念设定多组采样速率，对接收信号多次进行ESPRIT算法分解转换，实现到达角和频率联合估计的算法。文献[4]提出了一种树形结构(FSF(频率-空间-频率)-MUSIC)，即先后3次利用MUSIC算法对接收到的协方差矩阵进行转换，再利用子空间特性实现频率和DOA联合估计。但是这些算法只适用于低频段信号(<10MHz)。文献[5]利用均匀线阵将2-18GHz接收信号划分为多个1GHz的子频段，再对每路子信号以250MHz进行时间欠采样，利用PRO(专业)-ESPRIT生成的信息向量实现对信源的DOA与频率估计。但是该算法各子频段存在一定时延，且欠采样率不高。

近年来稀疏阵列吸引了越来越多的关注，这是因为稀疏阵列打破了原有传统均匀线阵的阵元布置限制，并利用差分思想扩展了阵列的有效孔径，提高了阵列自由度。如：较早提出的最小空洞阵列^[6]和最小冗余阵列^[7]均利用N个阵元将自由度提高至N²，但是这两稀疏阵列没有闭合的推导形式，因此未能进一步发展和应用。文献[8]提出的嵌套阵列虽然具有闭合推导形式并能实现对目标DOA估计，但该阵列只是实现了局部稀疏化，其中一部分阵元仍旧按照经典的均匀线阵紧密布置，因而阵元间耦合依旧比较严重。文献[9]基于互素感知理论提出了互素阵列，虽然在一定程度上改善了阵列稀疏度，但仍然存在局部阵元耦合问题。需要指出的是，上述该些稀疏阵列无法实现频率与DOA的联合估计。

发明内容

本发明提供了一种多信源频率与DOA联合检测的方法及其装置，本发明通过空时域欠采样，并结合谱校正^[10]与闭式鲁棒中国余数定理^[11]，实现了多对目标的高精度频率与DOA联合估计，节约了硬件成本，提高了数据利用率与估计精度，详见下文描述：

一种多信源频率与DOA联合检测的方法，所述方法包括以下步骤：

对L路信号样本做M_f点DFT，利用Tsui频谱校正器对DFT结果进行频率和相位校正，得到校正后的D组频率、相位、幅值的参数组，利用参数组构造出D个方向矢量；

将D个方向矢量按照最小距离进行分组匹配，得出D个信源参数匹配信息；

根据D个信源参数匹配信息构造频率余数组，并将频率余数组带入闭式鲁棒中国余数定理模型进行重构，得出频率估计值；

根据D个信源参数匹配信息得出L-1个相位差，并构造相位余数组，将相位余数组与重构模值组带入闭式鲁棒CRT重构得到中间参数，计算出DOA估计值。

所述L路信号具体为：

设置包含L个阵元的稀疏线性天线阵列，每个阵元位置布置两个ADC采样器；

各阵元的两个ADC采样器分别以f_s1,f_s2两种采样速率对入射信号进行并行欠采样，每个阵元采集的快拍数为M_f。

所述方法还包括：多信源参数分组匹配。

其中，所述方法多信源参数分组匹配具体为：

对于同一信源，构造得到的导向矢量距离非常近，即：

其中，D为检测目标数目；c_i为匹配分组；确定两路序列中D个谱峰对应情况，并将同一信源参数匹配分组。

其中，所述根据D个信源参数匹配信息构造频率余数组，并将频率余数组带入闭式鲁棒中国余数定理模型进行重构，得出频率估计值具体为：

根据多信源参数分组匹配、以及单信源频率估计获取多个信源频率估计。

所述根据D个信源参数匹配信息得出L-1个相位差，并构造相位余数组，将相位余数组与重构模值组带入闭式鲁棒CRT重构得到中间参数，计算出DOA估计值具体为：

根据多信源参数分组匹配、以及单信源DOA估计获取多个信源DOA估计。

所述单信源频率估计具体为：

对在各阵元上经过欠采样得到的L路采样样本作M_f点DFT得到X_l(k)，对其进行谱峰搜索，得到谱峰位置，进而获取左谱峰修正参数和右谱峰修正参数；

根据左谱峰修正参数和右谱峰修正参数得到左修正序列和右修正序列，并计算差值，得出校正后的谱峰位置，得出频偏估计值；

根据频偏估计值，得出相位、频率、幅值的校正估计值。

所述单信源DOA估计具体为：

将CRT重构得到的

和已得到的波长估计值

带入下式，即可得到到达角估计值：

其中，

为相位重构量；

为中间量。

一种多信源频率与DOA联合检测的装置，所述检测装置包括：稀疏阵列天线、ADC

采样器、DSP、输出驱动及显示装置，

当远场窄带信号入射到稀疏阵列体天线上，各阵元上的两个ADC采样器以不同步采样速率进行并行欠采样，将得到的数据输入到DSP器件中；

DSP器件进行处理，最终得到多个入射目标的频率与DOA估计值，并在输出驱动及显示装置上显示结果。

本发明提出的基于互素稀疏阵列的多目标频率和DOA联合测量方法，若将其应用于实际工程领域将具有以下功能优势：

1)降低硬件成本

本发明所提出的互素稀疏阵列，不需要按照传统均匀线阵将传感器按照半波长紧密排列，只需通过少数阵元即可实现多目标频率与DOA估计。同时，本发明打破了传统奈奎斯特采样速率限制，不需要进行下变频即可实现对高频段信号的采样，进一步降低了硬件要求。

2)提高数据利用率

本发明只需要利用各阵元上两个ADC(数模转换器)同时进行欠采样得到的数据，按照所提算法对其进行处理，即可得出频率与DOA估计结果，避免了多次采样与长时间采样。

3)测量精度高

由于本发明采用了频谱校正与闭式鲁棒中国余数定理，具有很高的估计精度。

附图说明

图1为单信源频率与DOA联合估计流程图；

图2为互素稀疏阵列天线布置示意图；

图3为频谱校正流程图；

图4为两种采样速率下参数对应示意图；

图5为频率检测成功概率和SNR(信噪比)关系曲线；

图6为DOA检测成功概率和SNR关系曲线；

图7为频率测量的均方误差关系曲线；

图8为DOA测量的均方误差关系曲线；

图9为本发明的硬件系统结构图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

实施例1

101：设置包含L个阵元的稀疏线性天线阵列，每个阵元位置布置两个ADC采样器；

其中，从某一特定时刻开始，各阵元的两个ADC采样器分别以f_s1,f_s2两种采样速率对入射信号进行并行欠采样，每个阵元采集的快拍数为M_f。

具体实现时，采样速率根据实际应用中的需要进行设定，本发明实施例对此不做限制。

102：对得到的L路信号样本做M_f点DFT，利用Tsui频谱校正器^[12]对DFT结果进行频率和相位校正，得到校正后的D组频率、相位、幅值的参数组，利用参数组构造出D个方向矢量；

其中，该校正器为本领域技术人员所公知，本发明实施例对此不做赘述。

103：将步骤102中得到的D个方向矢量按照最小距离进行分组匹配，得出D个信源参数匹配信息；

104：根据D个信源参数匹配信息构造频率余数组，并将频率余数组带入闭式鲁棒CRT(中国余数定理)模型进行重构，得出频率估计值；

105：根据D个信源参数匹配信息得出L-1个相位差，并构造相位余数组，将相位余数组与重构模值组带入闭式鲁棒CRT重构得到中间参数，计算出DOA估计值。综上所述，本发明实施例在空间欠采样情况下，借助互素稀疏阵列实现频率和DOA联合估计；在多阵元并行单次时间欠采样情况下实现频率和DOA联合估计；借助频谱校正和中国余数定理实现多目标源的频率和DOA联合估计。

实施例2

下面结合具体的计算公式、实例对实施例1中的方案进行进一步地介绍，详见下文描述：

201：阵元结构与信号模型；

如图1所示，布置阵元个数为L的互素稀疏阵列，且在每个阵元位置布置两个具有不同采样速率的ADC采样器(采样速率分别为f_s1＝M_fη₁,f_s2＝M_fη₂，其中η₁和η₂是互素整数，M_f为正整数)。

假设一载波频率为f₀的远场窄带信号以角度θ₀入射到天线阵列上，则第l个阵元与第l+1个阵元所接收到的信号相位差为：

其中，λ₀＝c/f₀为信号波长，d_l表示阵元间距，c为光速。

由于本方法提出的阵列结构稀疏度较高，阵元间距d_l>>λ₀/2，因而相位差会产生以2π为周期的模糊度，即：

其中，n_l是未知整数，称之为折叠整数(即对相位差除以2π会产生整数和余数，其中折叠整数可以看作是所得到的整数)。

将式(1)(2)合并可得：

接下来引进一个中间量：

其中，M_θ是一个正整数，{Γ₁,...,Γ_L}为互素整数组。阵元间距可以按下式确定：

考虑在实际工程中的应用情况，可以将入射的远场窄带信号表示为如下形式：

其中，s₀(t)为复指数信号；A₀为幅值。

由于噪声必可避免，因而第l个阵元所接受到的信号可以表示为：

其中，x_l(t)为接收信号；

为入射信号相位；ξ_l(t)为加性高斯白噪声。

202：单信源频率估计；

一、频谱校正流程如图2所示，现概括如下：

1)对在各阵元上经过欠采样得到的L路采样样本作M_f点DFT得到X_l(k)，对其进行谱峰搜索，得到谱峰位置k_c；

其中，X_l(k)＝DFT(x_l(t)),k＝0,...,M_f-1，令k_L＝k_c-0.5,k_R＝k_c+0.5；k_L为左谱峰修正参数；k_R为右谱峰修正参数；X_l(k)为做傅里叶变换得到的信号序列。

2)根据步骤1)得出的k_L，k_R计算

其中，M为采样长度；n为采样点数；X_L为左修正序列；X_R为右修正序列；x(n)为接收到的信号序列。

3)通过判定步骤2)得出的X_L和X_R大小关系，计算中间量β(即X_L和X_R取差值)，并得出校正后的谱峰位置

确定得出频偏估计值

其中[·]为四舍五入运算；

4)根据步骤3)得出的频偏估计值

得出相位、频率、幅值的校正估计值。

上述相位、频率、幅值的校正估计值的计算为本领域技术人员所公知，本发明实施例对此不做赘述。

二、单信源频率估计

如图3所示，在任意第l个阵元上，利用两个ADC对入射信号(式(7)所表示信号)分别以采样速率f_s1和f_s2进行欠采样，得到两个长度为M_f的序列：

其中，x_l,1(n)、x_l,2(n)为两种采样速率采样后得到的两路信号序列；ξ₁(n)为采样器的加性高斯白噪声；ξ₂(n)为另一采样器的加性高斯白噪声。

由于欠采样的缘故，上式中的f₀/f_s1和f₀/f_s2必然都大于1，因此：

其中，

为第一路信号(一采样器)产生的频率余数；

为第二路信号(另一采样器)产生的频率余数；

正整数。

则式(8)可以进一步表示为：

令

将其带入式(9)可得

其中，

为重构模值组；η_m为互素整数；

为重构余数；

为式(9)中的频率余数。

式(11)完全符合闭式鲁棒中国余数定理模型，其中

为模值，

和

为余数，n₁,n₂为折叠整数。

现将频率估计过程概括如下：

1)在任意第l个阵元上，利用两个ADC，分别以f_s1和f_s2两种采样速率对入射信号进行欠采样；

2)将步骤1)得到的采样序列按照图3所示的谱校正流程进行频谱校正与参数估计，分别得到3个参数的估计集合

3)利用步骤2)得到的估计参数构造频率余数

将各项参数带入闭式鲁棒中国余数定理模型进行重构，对得到L个频率粗略估计值做统计平均，得到频率估计值

203：单信源DOA估计；

如图3所示，在阵列上第l个阵元和第l+1个阵元所接收到的相位差可以表示为：

由于本发明实施例中阵元间距远大于半波长，因此会导致相位差中包含以2π的整周模糊，即

其中，n_l表示折叠整数，故

(mod表示模除操作)

将式(5)带入到式(14)中，有

其中，

为式(4)计算得到的中间量；Γ_l为互素数组中第l个素数。

令

显然

Q_θ＝n_lM_l+r_l,l＝1,...,L-1 (17)

将CRT重构得到的

和已得到的波长估计值

带入下式，即可得到到达角估计值：

其中，多信源频率估计的过程与步骤202中单信源频率估计、以及多信源DOA估计的与步骤203中单信源DOA估计的主要区别在于多信源估计中需要解决参数分组匹配问题，其他过程的原理是相同的。

步骤202和步骤203中详细论述了单信源频率估计、以及单信源DOA的估计原理，单信源频率估计与DOA估计是基础，在理解了单信源频率与DOA估计的基础上，进一步理解步骤204中的多信源参数是如何分组匹配的，即可明白本发明实施例提供的多信源估计原理。

204：多信源参数分组匹配。

在多信号入射情况下，对接收到的序列做DFT后得到的谱线中包含多个谱峰，由于每个阵元上布置了两个具有不同速率(f_s1,f_s2)的ADC采样器，同一信源在DFT谱X_l,1(k),X_l,2(k)中的谱峰位置是完全不同的，如图4所示。因此，需要将相同信源的参数信息进行分组匹配。

对于同一信号，其幅值与相位参数是成对的(如式(7)所示)，利用上述特性可以识别同一信源在两路ADC采样序列上的对应谱峰位置。同时为了克服非零频偏导致的DFT谱线泄露，本发明实施例利用加窗DFT(hannning窗)和频谱校正对DFT谱线进行处理。

将通过频谱校正器后得出的参数对

构造导向矢量，

为幅值估计值；

为相位估计值，令：

ADC1得到的DFT序列X_l,1(k)中可得出一个导向矢量

同理ADC2得到的DFT序列X_l,2(k)中可得出另一个导向矢量

其中，

为第一路信号序列中第i个估计幅值；

为第一路信号序列中第i个估计相位；

为第二路信号序列中第p个估计幅值；

为第二路信号序列中第p个估计相位。

对于同一信源，其构造得到的导向矢量距离会非常近，即：

其中，D为检测目标数目；c_i为匹配分组。

最终确定两路序列中D个谱峰对应情况，并将同一信源参数匹配分组。接下来参照单信源估计即可得到多目标的频率与DOA估计值。

综上所述，本发明实施例在空间欠采样情况下，借助互素稀疏阵列实现频率和DOA联合估计；在多阵元并行单次时间欠采样情况下实现频率和DOA联合估计；借助频谱校正和中国余数定理实现多目标源的频率和DOA联合估计。

实施例3

按照图1所示布置阵元数目L＝3的互素稀疏阵列天线，设置单位波长λ＝0.15，M_θ＝2，对应{Γ₁,Γ₂}＝{5,6}，根据式(5)得出阵元间距为d₁＝0.9,d₂＝0.75。设置{η₁,η₂}＝{5641,5647}，M_f＝512，对应的ADC1余ADC2的采样速率分别为：f_s1＝2888192Hz,f_s2＝2891264Hz。

假设4个远场窄带信号入射到阵列上，其参数设置如下表所示：

首先对本发明的多目标估计能力进行检验，令SNR＝25dB，估计结果如下表所示：

实验结果表明本发明能够准确估计4个信源的频率与DOA，并具有较高的估计精度。

下面分别用检测概率和均方根误差(root mean square error,RMSE)对本发明的参数估计性能进行考察。假设各阵元接收特性相同且相互独立，噪声为零均值的高斯白噪声。

设置SNR∈[5:25]dB，在每个SNR情况下进行1000次Monte Carlo实验。以第一个信源为检验目标，当频率估计值满足

时则认为测量成功,否则视为失败.同理,当DOA估计值

满足

时视为检测成功。图5和图6分别给出了频率检测概率和DOA检测概率随SNR变化的曲线。如图5所示，当SNR＞22dB，频率估计成功概率达到100％。如图6所示，当SNR＞24dB，DOA检测成功率达到了100％。

设置SNR∈[20:40]dB，在每个SNR情况下进行1000次Monte Carlo实验。图7和图8分别给出了本发明实现频率与DOA估计的均方根误差随SNR变化的曲线。如图7所示，当SNR＞26dB，频率估计的均方根误差达到10Hz；如图8所示，当SNR＞26dB，DOA估计的均方根误差下降到0.1。

实施例4

本发明实施例是与实施例1和2中的方法部分对应的装置部分，参见图9，该装置具体为：

如图9所示，本发明实施例所提出的互素稀疏阵列的频率与DOA联合测量装置包括：稀疏阵列天线、ADC采样器、DSP、输出驱动及显示装置。

测量过程如下：当远场窄带信号入射到稀疏阵列体天线上，各阵元上的两个ADC采样器以不同步采样速率对其进行并行欠采样，将得到的数据输入到DSP器件中，按照实施例1和2中所提出的算法进行处理，最终得到多个入射目标的频率与DOA估计值，并在显示装置上显示结果。

本发明实施例中，DSP作为本装置的核心元件，其内部程序存储器所储存的算法是整个测量系统的精度、复杂度和稳定性的主要因素，在整个测量过程中承担以下功能：

1、运行核心算法，对各阵元ADC采样得到的信号序列进行频谱校正，并调用闭式鲁棒中国余数定理，从而实现频率与DOA的检测；

2、将测量结果输入到驱动和显示模块。

参考文献

[1]L.Xu,J.Li,and P.Stoica,Target detection and parameter estimationfor MIMO radar systems[J].IEEE Transactions on Aerospace and ElectronicSystems,2008,44(3):927–939.

[2]Y.Li,N.SESHADRI,S.ARIYAVISITAKUL,Channel estimation for OFDMsystems with transmitter diversity in mobile wireless channels[J].IEEEJournal on Selected areas in communications,1999,17(3):461–471.

[3]LEMMA A N,VAN DER VEEN A-J,DEPRETTERE E F.Analysis of joint angle-frequency estimation using ESPRIT[J].IEEE Transactions on Signal Processing,2003,51(5):1264-83.

[4]LIN J-D,FANG W-H,WANG Y-Y,et al.FSF MUSIC for joint DOA andfrequency estimation and its performance analysis[J].IEEE Transactions onSignal Processing,2006,54(12):4529-42.

[5]ZOLTOWSKI M D,MATHEWS C P.Real-time frequency and 2-D angleestimation with sub-Nyquist spatio-temporal sampling[J].IEEE Transactions onSignal Processing,1994,42(10):2781-94.

[6]TAYLOR H,GOLOMB S W.Rulers part I[J].Univ.Southern Calif.,LosAngeles,CSI Tech.Rep,1985(85-05):01.

[7]MOFFET A.Minimum-redundancy linear arrays[J].IEEE Transactions onantennas and propagation,1968,16(2):172-175.

[8]PAL P,VAIDYANATHAN P P.Nested arrays:A novel approach to arrayprocessing with enhanced degrees offreedom[J].IEEE Transactions on SignalProcessing,2010,58(8):4167-4181.

[9]VAIDYANATHAN P P,PAL P.Sparse sensing with co-pprime samplers andarrays[J].IEEE Transactions on Signal Processing,2011,59(2):573-86.

[10]ZHANG F,GENG Z,YUAN W.The algorithm ofinterpolating windowed FFTfor harmonic analysis of electric power system[J].IEEE Transactions on powerdelivery,2001,16(2):160-164.

[11]LI X,LIANG H,XIA X-G.A robust Chinese remainder theorem with itsapplications in frequency estimation from undersampled waveforms[J].IEEETransactions on Signal Processing,2009,57(11):4314-22.

[12]Tsui J B.Digital techniques for wideband receivers[M].Raleigh:SciTech Publishing,2004:449-461.

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。