时−空欠采样下多个信号的频率和DOA 联合估计

2019-06-24郜丽鹏李佳林

应用科技 2019年3期

郜丽鹏，李佳林

哈尔滨工程大学信息与通信工程学院，黑龙江哈尔滨150001

随着电子战技术的日益发展，电磁环境也变得越发复杂，针对多信号的多参数估计是雷达、声呐、电子战等领域的研究热点，能够快速、准确地侦察检测多个雷达信号的频率、方向信息，对于电子对抗、雷达干扰技术而言至关重要[1−2]。在经典的空间谱估计算法MUSIC和ESPRIT的基础上提出的改进算法可以完成对多个信号的频率和DOA联合估计及配对[3−4]，但是这些算法的运算量都较大，且均需要满足信号在时域和空域上的理论前提。时域方面，香农采样定理表明，奈奎斯特采样率是被采样信号的最大频率的2倍，比如为了适应电子战2~18GHz的工作频段，需要的采样频率会很高，必须相应地增加模数转换器（analog-to-digital converter，ADC）的硬件成本和整体结构的功耗，甚至在某些特定情况下是不可实现的。同样，在空域方面若要适应电子战的工作频段，为满足相邻阵元的间距不得大于信号波长的二分之一的采样定理，需要密集的安装阵元天线，这样会带来相邻阵元间的互耦和硬件安装困难等问题，进而影响测向精度。因此，亟需寻求一种在时−空域欠采样条件下，能够对多信号多参数（频率、DOA）联合估计的新方法。

文献[5−6]提出了一种基于闭式中国余数定理（CRT）的频率、DOA联合估计方法，建立了基于闭式CRT的频率和DOA估计的模型，利用不同速率ADC采样样本（离散傅里叶变换，Discrete Fourier Transform）变换后搜索的谱峰位置构建余数完成频率估计，通过谱峰位置求得信号到达每个阵元的初相位及其相位差构建余数完成DOA估计，实现了高精度的联合估计。但该算法是针对单目标信号的参数估计，限制了其应用场合。文献[7]在闭式CRT的频率和DOA估计模型的基础上，提出了针对多目标信号的参数估计，利用信号的幅度不一致的信息，在多个信号经不同速率ADC采样样本DFT变换后得到的多个谱峰中，寻找属于同一信号的一组谱峰以构建频率余数和相位余数完成多个信号的频率和DOA联合估计，但该方法在有2个及2个以上信号的幅度一致的条件下会失效。文献[8]提出了一种能够重构多个整数的中国余数定理，并且将其应用在窄带信号的频率估计上，但其未考虑当信噪比较低时或样本数较少时谱峰位置估计不准确的问题，导致其测频精度的降低。

本文建立了一种基于重构多个整数的CRT算法和AM估计器的多信号的频率、DOA参数联合估计模型，利用多整数估计的CRT算法对稀疏线阵欠采样后的少量样本做DFT、谱峰校正，余数构建等，完成多个信号的频率估计，并根据估计完成的频率值和已经构建的余数信息在多组相位差中确定每个信号对应的一组相位差，完成频率和相位差的参数配对，进一步使用闭式CRT算法鲁棒地估计出多个信号的DOA值。本文算法只需对各个阵元接收到的多个信号做一次并行欠采样，并且谱校正后得到的谱峰信息可同时用于频率估计和DOA估计，满足目前复杂电磁环境下需要快速参数估计的需求。仿真结果证明了本文算法的高精度和低信噪比下的顽健性。

1 重构多个整数的顽健CRT算法

模数 {m1=ΓM1,m2=ΓM2,···,ms=ΓMs}由 s个整数构成，其中 M1,M2,···,Ms互素。假设 n个整数X1,X2,···,Xn均匀分布在( L,G]内，其中G

1.1 余数误差情况下的顽健重构算法

首先考虑只有一个带重构整数的情况，先前许多 CRT算法[9−11]的理论前提为|r˜il−ril|=然而在实际应用中，有可能会发生的情况，这会导致进一步使得CRT算法失效。文献[12]给出了这种问题的解决方法，其主要思想是将每个模数中的误差 ∆il调整为一致。首先根据 ml=ΓMl，可以得到：

式中 β为整数，因此

可以发现，只要重构出 Q就可以求得待重构的整数 X。接下来将重点放在如何鲁棒、顽健地重构整数 Q 上。分别用 a和 b表示集合最小和最大的元素。接下来分以下几种情况进行讨论：

1.2 重构多个整数的算法

对于多个整数的情况，为了方便讨论，假设待重构整数的集合 {Xi}(i=1,2,···,n)和模数的集合{ml|l=1,2,···,s}都是以升序排列的，会得到 n×s个余数 {ri,l}，如文献[13]所提到的，2个及2个以上的整数重构的难点是整数和余之间的对应关系是未知的，即由每个模数得到的 n个余数与 n个整数的对应关系是未知的，所以无法使用CRT算法进行重构。

对于这种问题，首先概述主要的算法思想。由于最大的模值 ms>d+2δ，所以可以在区间[X1− δ,Xn+δ]找到 X ，并且通过余数集 r˜is和 ⟨X⟩ms正确地估计 {Xi−X}的值。因∑此，选择估计待重构整数集的平均值 ¯，即X同样地，在整数集与余数集对应关系未知的前提下，可以得到每个整数模除以同一模值的余数的和并通过 r˜和上述的CRT算法顽健的重l构整数集的和为了顽健地重构需要假设 M并且属于同一模值的余数误差的总和需要满足

接下来就是如何通过 X˜¯ 和 r˜is顽健地重构各个整数 Xi。由于考虑到若有3个整数A、 B 和 m，如果 |A−B|

引理对于任意的 f∈R ，在区间 [f,f+ms)内至多有 n个元素 κp,κp+1,···,κp+n−1，其中p∈{1,2,···,h−n+1}。

关于引理的证明见文献[8]。

2 欠采样的频率估计方法

首先给出系统模型，如图1所示，由 s(s>2)个阵元天线组成非均匀线阵，假设有 m个同时到达的远场信号，不失一般性，将其表示为单频复指数形式：

式中： Ai和 fi分别为各个信号的幅度和频率；φi为信号的初相位； i=1,2,···,m。假设各个信号的频率fi都不相同，由于噪声和阵列阵元的间距，第l个阵元接收到的信号为：

式中： φil为第i个信号到达第l个阵元的初相位；ζl(t)为高斯白噪声； l=1,2,···,s。从任一时刻开始，s个阵元分别以 fs1∼ fss的采样速率对同时到达的 m个信号进行欠采样，并且得到了 s路样本，假设每路的样本数为 N ，则第l个阵元得到的样本序列为：

图1 稀疏阵列模型

对 s路欠采样样本序列做DFT，并进行谱峰搜索，得到谱峰位置 kil，则各个信号的频率 fi可表示为：

式中： nil为模糊倍数； δil为频偏小数，且 |δil|⩽0.5，也是影响测频精度的重要因素，δil可由AM估计器进行校正，进而得到更精确的余数；i=1,2,···,M;l=1,2,···,s。可将式（3）转变为中国余数定理模型

式中：待测的信号频率 fi对应CRT中的被除数 Xi；采样频率 fsl对应CRT中的模数 ml=ΓMl，Γ为最大公约数， Ml为互素的整数；i =1,2,···,M ；l =1,2,···,s。而重构所需的余数为

对于 m个信号的情况，在 s路样本序列DFT变换后，都会搜索得到 m个谱峰，而在 m个信号的幅度有2个及2个以上是相同的情况下，不能在每路 m个谱峰中寻找对应同一信号的 s个谱峰，即不能构造出 m组余数 {ri1,ri2,···ris}，进而分别进行 m次闭式CRT算法重构出 m个信号的频率。

考虑到第1章提出的重构多个整数的CRT算法，其前提条件为 ms>d+2δ，即模数组中的最大模值需要大于多个整数中最大值与最小值的差值，结合本节的频率估计模型，即最大的欠采样速率 fss需要大于信号带宽（ fm−f1），所以该CRT算法适用于窄带信号中的多个频率估计，算法流程如下：

1）对 s路欠采样样本做 N点的DFT变换，谱峰搜索后得到谱峰位置 kil， i=1,2,···,m ， l=1,2,···,s，即每路样本得到 m个谱峰位置。

2）校正每路样本的谱峰位置对应的频偏小数，对于第l 路欠采样样本，将 m频偏小数的值进行初始化=0， i=1,2,···,m，按式（5）进行计算得到迭代值

式中： q=±0.5； Re(·)表示取其实部。将其迭代2次以上的收敛值做为频偏小数的估计值，i=1,2,···,m ，l=1,2,···,s。

3）将谱峰位置 kil、频偏小数代入式（4）得到校正后带有误差频率余数估计值 r˜il， i=1,2,···,m，l=1,2,···,s。

4）计算对应同一模数 ml的余数估计值的和

3 欠采样的DOA估计方法

3.1 基于CRT算法的DOA估计模型

阵列结构如图1所示，由相位干涉仪原理可知，由信号到达阵元的波程差会产生相位差，在任一时刻，阵元 j 和阵元 j+1的接收到的信号 Si的相位差为：

式中：i=1,2,···,m；j=1,2,···,s−1；λi为信号Si的波长； θi为信号Si的入射角度；dj为阵元间距，当则会产生相位模糊问题，即：

式中：nij为模糊倍数； ∆φij为真实情况下测得的相位差。由式（6）、（7）可得：

式中：C为非负常数；Mθ为任意正整数； Γθ1∼ Γθ(s−1)为互素的整数。令待重构的整数Q别定义模为相位差构造的对应CRT中的余数为

所以，只要我们得到每个信号在阵元间的相位差，就可以通过CRT算法重构出Qθi，在频率值和相位差余数配对后，通过式（9）求得DOA的估计值，其中 λ0的值可由第2章中测得频率结果得到。

3.2 相位差余数与信号频率的配对

考虑到DFT变换后谱峰位置对应的相位为信号Si到达阵元l的初相位，记为 φpil。由第2章分析可知，对于m个信号同时到达s个阵元的情况，可以通过s路样本序列得到m×s个谱峰位置。同样我们可以得到m×s个初相位，即每个阵元可以得到m个初相位。只有在相邻阵元的初相位中，选择属于同一信号频率的2个初相位，才能得到正确的相位差，进而得到正确的DOA估计值。

不同整数Xi模除同一模值m得到的余数为ri，定义P0为余数ra和rb不相等的概率，其中i=1,

表1中s为模值的个数。由式（10）可知，当模值m相对较大时，余数ra和rb不相等的概率极大，特别是对应第2章中的欠采样采样频率，一般都在100MHz以上，此时P0几乎为1。所以，可以认为，对于任一阵元的样本序列sl(n)，通过DFT变换后的谱峰位置得到的频率余数值

式中i=1,2,···,m。

表1 s个阵元采集m个信号得到的余数及对应初相位

如表1所示，每列的余数都是不相同的，同时，我们由谱峰位置得也可以到对应的初相位 φpil。

在频率估计后，m个信号的频率估计值f˜i和s个阵元的采样频率fsl均为已知的，故可以通过f˜imodfsl，l=1,2,· · ·,s，得到信号Si的s个余数ril，并依次用余数ri1∼ris，在表中每列找到r˜i1∼r˜is对应的的位置，进而得到信号Si到达s个阵元的s个初相位 φpi1∼ φpis，并求得信号Si的s−1个相位差 ∆φij，j=1,2,···s−1。同理，依次可以得到m个信号的m×(s−1)个相位差。所以此时第1章的重构多个整数的C RT算法适用，而考虑到该算法的理论前即要求通过谱峰位置得到所有信号到达同一阵元的初相位的误差的和小于DOA估计模型中的这个条件对于相位测量是比较苛刻的。但与频率估计模型不同的一点是，对于DOA估计来说，经过上述方法的参数配对，属于同一信号带有DOA信息的待重构整数Qθi的

s个相位差余数个待重构整数与余数集的对应关系是已知的。所是已知的，也就是每以我们可以分别使用m次传统的闭式CRT算法求得m个整数Qθi的值，这样相较于直接重构多个整数的算法，我们就把相位测量误差的限度降低了m倍。故DOA估计算法流程如下：

1）对于第l路样本，根据谱峰位置对应的幅值求得个信号到达阵元l的初相位 φpil,i=1,2,···,m，并由频率估计过程中得到的频偏小数通过式（11）对信号初相位进行校正：

式中 i=1,2,···,m。

式中 j=1,2,···,l−1。

3）通过相位差 {∆φij,j=1,2,···,l−1}构成的相位差余数 {rθij,j=1,2,···,l−1}和闭式CRT算法求得信号 Si的 DOA估计值θ˜i，分别对 m组相位差余数使用CRT算法即可求得 m个信号的DOA估计值。

4 仿真实验

设有3个远场信号同时到达阵列，阵列由4个阵元组成，假设信号的频率在1 500∼2150MHz，根据本文提出的欠采样下的频率和DOA联合估计方法，最大模值需大于多个整数中最大值与最小值的差的条件。阵元中最大的采样频率 fs4应设置为 650MHz，选取对应CRT算法参数中的最大公约数 Γ=50×106，则一组互素模数中最大模数M4=13。为满足时域欠采样的条件，选取3个小于 M4的模数构成互素的模数组[M1,M2,M3,M4]=[3,5,11,13]，所以对应的4个阵元ADC的采样频率分别为 fs1=150MHz、 fs2=250MHz、 fs3=550MHz、fs4=650MHz。显然阵元ADC的采样速率远小于信号的频率，属[于时域欠]采样。非负常数 C取为0.1，互质参数 Γθ1,Γθ2,Γθ3=[3,5,7]，即阵元间距为[d1,d2,d3]=[0.35m,0,21m,0.15m]，小于入射信号的半波长，属于空域欠采样。

本节实验采用检测概率 pd和均方根误差（root mean square error,RMSE）对算法的顽健性、抗噪性以及参数测量精度进行考量。假设每个阵列的接收特性相同并且是相互独立的，欠采样样本序列的长度 N=2000，噪声是均值为零的高斯白噪声，在每个信噪比（SNR）条件下进行1000次Monte Carlo实验。

4.1 仿真实验1

在时域欠采样的条件下，单独对多信号的频率估计的顽健性、抗噪声性能进行考察。在该仿真中，固定3个信号的入射角均为38°，设信号的频率 fi,i=1,2,3，在15002~2150MHz内随机选取。

图2 不同信噪比下的频率检测概率

4.2 仿真实验2

接下来验证多个信号的频率、DOA联合估计算法的精度性能。该仿真中，固定3个信号频率和入射角度 Si[fi,θi]分别为且同时到达阵列，基于AM估计器的频率估计的RMSE仿真曲线如图3所示。从图3可以看出，本文提出的基于多整数估计CRT的频率估计算法在各个信噪比下均有着较高的精度，当信噪比大于10dB后，3个信号的频率估计误差均小于1kHz，即小于信号频率的

图3 不同信噪比下3个信号频率估计的均方根误差

且由于信号频率估计的精确度较高，我们通过频率估计值和ADC采样率可以得到准确的余数，并通过频率估计中谱峰位置得到的余数正确配对同一信号到达4个阵元的初相位。图4为在信号和相位参数配对后，基于闭式CRT算法分别对3个信号的DOA估计的RMSE曲线。可以看出，在各个信噪比条件下均表现出较高的DOA估计精度，并且在(0,90°]测角范围内，信号入射角度越大，DOA估计误差越大，与理论DOA估计误差一致。

图4 不同信噪比下3个信号DOA估计的均方根误差

4.3 仿真实验3

考察频率和DOA联合估计时，DOA估计的顽健性和抗噪声性能。在该仿真中，固定3个信号的频率分别为1.573、1.900、2.016GHz，信号角度在(0,90°]内随机选取。若DOA估计值则认为检测成功；否则为检测失败。在DOA估计过程中，采用重构多个整数的CRT和独立3次使用闭式CRT2种算法做比较，得到的检测概率随信号比SNR变化的仿真曲线如图5所示。正如3.2节的讨论，重构多个整数的CRT算法在DOA估计时对允许的相位测量误差的限度较低，其在低信噪比条件下的性能不佳，而在参数配对后单独3次使用闭式CRT算法在信噪比大于−9dB后，检测概率均能达到90%以上，可以看出其较好的顽健性。

图5 不同信噪比下DOA的检测概率

4.4 仿真实验4

最后对比基于MUSIC的空时二维谱频率和DOA联合估计算法的性能和算法耗时。与仿真2相同，固定3个信号频率和入射角度 Si[fi,θi]分别为33°]，且同时到达阵列。对于空时二维谱算法，阵列阵元间距保持不变，4个阵元的采样速率均为650MHz，满足时域和空域为欠采样的条件。角度在(0,90°]的范围内以0.1°为步长，频率在(1500MHz,2150MHz]以 0.1MHz为步长搜索得到二维谱如图6所示。

图6 基于MUSIC的空时二维谱

可以看出由于时域和空域的欠采样，空间二维谱会出现伪峰，无法对3个入射信号进行正确的频率和DOA估计。且基于MUSIC的空时二维谱的估计算法的仿真时间为71.985s，并且会随着搜索步长的减小而增大。而本文提出的时−空欠采样下的频率和DOA联合估计算法的仿真时间为0.037s，极大地减少了算法运算时间，适合快速时变的工作场景。