基于L 型延迟阵列调制宽带转换器的信号载频和二维到达角联合估计*

姜思仪 付宁 乔立岩 彭喜元

(哈尔滨工业大学电子与信息工程学院, 哈尔滨 150001)

1 引 言

阵列信号处理已被广泛应用于雷达[1]、声纳[2,3]和无线通信[4,5]等领域, 通过按一定样式在空间中分布的传感器阵列接收空间电磁信号, 并进行一系列处理得到空时频各域的参数.在实际应用中, 频率 和 二 维 到 达 角(two-dimensional direction of arrival, 2D-DOA)是对电磁波进行识别的重要特征, 对载频和2D-DOA 的联合估计在阵列信号处理领域备受关注.

随着信号频率范围越来越宽, 基于传统奈奎斯特定理的信号采集方法给采样、存储和处理设备带来了巨大的压力.近来新兴的压缩感知(compressed sensing, CS)理论[6-9]能够在少量样本的情况下准确地恢复稀疏信号, 该理论大大降低了采样率, 具有很高的应用价值.目前已有不少学者将CS 理论应用到模拟信号的采样中[10-15].其中, 由Mishali 和Eldar[14,15]提出的MWC 系统是一种专门针对多频带模拟信号采样的模拟信息转换系统.

为了降低阵列中多个通道同时采样产生的庞大数据量, 近来许多学者提出了基于欠采样的载频和到达角联合估计方法, 黄翔东等[16]利用互素稀疏阵列, 结合闭式中国余数定理实现了欠采样下频率和DOA 的联合估计.沈志博等[17]提出一种CS奇异值分解算法, 但配对步骤较为繁琐.Liu 等[18,19]提出了一种基于多陪集结构的联合估计方法, 分别基于子空间分解和正则分解/平行因子分析(CANDECOMP/PARAFAC, CP)[20]方法, 但接收阵列结构较为复杂.接着, 他们又进行了扩展[21,22], 利用稀疏阵列减少了所需的阵元个数, 但这种结构对ADC 仍有比较高的要求.Stein 等[23,24]将MWC技术与阵列信号处理相结合, 提出了L 型阵列MWC 结构来接收信号, 并提出了一种基于旋转不变子空间(estimating signal parameter via rotational invariance techniques, ESPRIT)方法的联合载波DOA 恢复算法.Cui 等[25]提出了一种基于均匀线性MWC 阵列由两阶段估计和参数配对构成的估计方法, 算法复杂且精度有限.Chen 等[26]提出了一种基于多个MWC 结构的DOA 和载频联合估计方法, 鲁棒性较好, 但该结构所需的通道数较多, 硬件成本高.然而在实际中电磁波存在于三维空间中, 使用2D-DOA 来描述目标的方向更为准确, 但目前对基于欠采样的载频和2D-DOA联合估计研究较少.陈玉龙等[27]提出了一种基于L 型阵列的2D-DOA 估计的欠采样算法, 但所需的天线阵元个数较多.Esmaeil 等[28]提出了一种基于双L 型阵列结构的网格CS 估计方法, 但是参数估计精度受网格的大小限制.他们还提出了一种贝叶斯压缩感知算法[28], 方法不依赖网格但步骤复杂且计算量较大.

从以上分析可以发现, 现有的方法大多都只考虑一维DOA 的估计, 需要假设目标同时位于一个平面上, 然而在实际应用中目标一般分散在三维空间中, 因此需要2D-DOA 估计.而基于欠采样的载频与2D-DOA 联合估计方法研究较少, 现有的方法中存在结构冗余, 精度受网格限制和算法复杂等问题.因此在压缩采样条件下对信号的2D-DOA和频率参数联合估计还有待研究.

L 型阵列结构有设计简单, 阵列方向图主瓣窄, 估计性能较好等优点, 因此本文提出了一种基于MWC 技术的L 型延迟阵列结构, 在L 型阵列结构的基础上增加延迟通道来增加可估计到达角的维度, 利用延迟通道与未延迟通道的采样值之间的已知的相位差可直接估计载频, 进而计算2DDOA 参数, 无需额外的载频和2D-DOA 的参数配对操作, 避免了配对所引入的误差和复杂度的提升.并结合L 型延迟阵列结构的特点构造相关矩阵和三线性模型, 提出了两种联合载波频率和2DDOA 估计方法.其中一种基于ESPRIT 算法, 另一种基于CP 分解技术.基于ESPRIT 的方法有计算量小, 无需谱峰搜索的优点, 适用于实时性要求较高的应用场合; 基于CP 分解技术的方法利用三维数据的结构特性, 有鲁棒性较好的优点, 适用于信噪比较低的应用场合, 本文主要基于这两种方法进行研究.本文第2 节介绍了信号模型; 第3 节首先介绍了L 型延迟阵列接收结构, 接着给出了接收信号模型; 第4 节首先对所提的两种方法的原理和步骤进行了详细的阐述.接着给出了参数估计条件, 最后对复杂度进行了详细的分析.第5 节给出了本文所提算法的验证结果以验证该算法的正确性和有效性.

2 问题描述

考虑M个远场目标发出的窄带信号si(t)ej2πfi,其中i=1,2,··· ,M.信号的频率范围为F=[-fNyq/2,fNyq/2] , 其中fNyq为奈奎斯特频率.si(t)为互不相关的基带信号, 带宽不超过B且远小于调制载频fi.假设信源发出的信号以方位角θi ∈(-90°,90°) 和俯仰角φi ∈(0°,90°) 的入射方向被天线阵列接收, 为了避免空间模糊, 假设载频、方位角与俯仰角满足[24]:

本文的目标是设计一种基于MWC 技术的阵列结构, 使其能够从欠采样样本中完全估计上述信号的载频, 俯仰角和方位角参

数.

3 L 型延迟阵列与接收信号模型

3.1 L 型延迟阵列结构

提出图1 所示的L 型延迟天线阵列结构, 利用延迟通道的采样值可以先估计载频, 进而计算2D-DOA 参数, 避免了额外的参数配对操作.该结构由两个分别沿x轴和y轴正方向均匀分布的N个线性天线阵列{x1,x2,··· ,xN}和{y1,y2,··· ,yN}组成, 两个轴在原点共用同一个天线阵元.相邻阵元之间的间距d满足d≤c/fNyq, 其中c为光速.

图1 L 型阵列MWC 结构图Fig.1.L shaped array MWC.

每个天线阵元后依次连接混频器、低通滤波器和采样模块, 类似MWC 通道结构.在x轴的每个阵元后增加一个延迟通道, 在混频之前对接收信号进行固定延时τ, 如图2 所示.

图2 x 轴延迟通道结构图Fig.2.x-axis delay channel.

3.2 接收信号模型

定义xn(t)和yn(t) 分别为x轴和y轴的第n个天线阵元接收的信号,(t) 为xn(t) 经过延迟τ后的信号, 由于目标信号为窄带信号, 有si(t+τ)≈si(t)[24], 则接收信号可以表示为

图3 接收信号频谱示意图Fig.3.Spectrum of received signal.

阵元接收信号在混频器与周期为Tp=1/fp的± 1 伪随机序列p(t) 相乘后, 再经过截止频率为fp/2 的低通滤波器, 只留下Fp=[-fp/2,fp/2] 内的频谱.滤波后的信号经过傅里叶变换可得:

滤波后的信号以fs=fp的速率进行低速采样得到观测值, 将其傅里叶变换后写为矩阵形式可得:

其 中X(f) 为N×1 的矩阵, 第n个元素为Xn(ej2πfTs).W(f) 是M×1 的 矩 阵, 第i个 元 素为矩阵为阵列流型矩阵, 其中

定义x轴通道的观测值为x[k]=[x1[k],··· ,xN[k]]T, 对(4)式取离散时间傅里叶逆变换, 在时域有:

类似地, 对于y轴和x轴延迟通道的采样值y[k]=[y1[k], y2[k]··· , yN[k]]T和有:

本文的目的是从欠奈奎斯特采样值x[k],y[k]和x′[k] 中 估计信号的载频俯仰角和方位角参数.

4 载频与二维DOA 联合估计方法

本节给出了基于图1 接收结构所提的两种载频与二维DOA 联合估计方法, 一种方法基于ESPRIT 技术, 另一种方法基于CP 分解.对两种方法的原理和步骤进行了详细的阐述.接着给出了参数估计条件, 最后对复杂度进行了详细的分析.

4.1 基于ESPRIT 的估计方法

首先给出基于ESPRIT 的参数联合估计方法,将采样值分为前N—1 项和后N—1 项组成的两个子矩阵, 并构造如下的相关矩阵:

其中Ax1和Ay1分别为Ax和Ay的前N—1 行,

构造如(9)式的矩阵R,

对R进行奇异值分解, 得到前M个非零奇异值对应的左奇异向量U=[U1;U2;U3;U4;U5;U6].U是一个 6 (N -1)×M的矩阵,Ui为 (N -1)×M的矩阵,i=1,2,··· ,6.通过推导可得矩阵Ui和对角阵Φx,Φy和Ψ之间的关系, 定义矩阵V1,V2和V3,

对矩阵 (V1+V2+V3) 进行特征值分解, 这样可以保证矩阵Φx,Φy和Ψ中元素顺序对应, 可得到特征向量矩阵, 进而通过(13)式—(15)式计算可以得到

目标信号的载频fi, 方位角θi和俯仰角φi就可以计算,

其中 ∠ (·) 为求复数的相位角, 综上所述, 现将基于ESPRIT 的估计方法的步骤总结为如下.

输入各通道Q快拍采样值x[k] ,y[k]和x′[k] ,其中k=1,··· ,Q;

输出载频俯仰角和方位角

步骤1根据(8)式和(9)式计算协方差矩阵R;

步骤2对R进行奇异值分解, 得到左奇异向量U, 根据(10)式—(12)式计算矩阵V1,V2和V3;

步骤3对矩阵 (V1+V2+V3) 进行特征值分解, 根据(13)式—(15)式计算得到

步骤4根据(16)式—(18)式计算载频俯仰角和方位角

4.2 基于CP 分解的估计方法

CP 分解[17]可以将一个高维的张量分解成多个秩为1 的张量的和, 从而降低待处理参数的维度.为了使(9)式中的矩阵R满足CP 分解的模型,定义一个大小为 6×M的矩阵R=[R1R2···RM] , 其第i列向量Ri为

由于信号满足互不相关, 因此自相关矩阵Rw中除对角线元素外均为0.定义一个三阶张量χ(N-1)×(N-1)×6, 其正向切片Xk=Rk, k=1,2,...,6,Rk可以表示为

通过观察可以发现(20)式符合典型的CP 分解模型, 阵列流型矩阵Ax1和Ay1均为范德蒙德矩阵, 当N≥M+1 时,Ax1和Ay1均 为满k-秩 等于M.由于信源满足(1)式, 由(19)式的定义可以看出,矩阵R一定有任意两列独立, 即kR≥2.因此有

满足CP 分解的唯一性条件[17], 在不考虑列置换模糊和尺度模糊的条件下, 因子矩阵Ax1,Ay1和R可以从张量χ中唯一恢复.通过固定除一个因子矩阵外的所有矩阵的方式, 将CP 转化为以下的线性最小二乘的问题:

其中X(n)为张量X的n-模式矩阵化, 因子矩阵和分别为Ax1,Ay1和R的估计值.利用交替最小二乘法进行循环求解, 可解得三个因子矩阵.对中的元素计算可得(23)式—(25)式, 为了方便将其记为

目标信号的载频fi, 方位角θi和俯仰角φi参数就可以计算:

综上所述, 现将基于CP 分解的估计方法的步骤总结为如下.

输入各通道Q快拍采样值x[k],y[k]和x′[k] ,其中k=1,··· ,Q;

输出载频俯仰角和方位角

步骤1根据(8)式和(20)式计算三线性模型χ;

步骤2利用交替最小二乘法对三线性模型χ进行求解, 可求得因子矩阵

步骤3利用中元素根据(23)式—(25)式计算得和

步骤4根据(26)式—(28)式计算载频俯仰角和方位角

4.3 参数估计条件

由于L 型延迟阵列MWC 结构是基于MWC理论[14,15]展开的, 因此估计条件有一定的一致性,结合文献[14]中MWC 理论的重构条件给出如下的参数成功估计的条件.

对于满足(1)式的M个互不相关的窄带远场信号使用图1 所示的共有2N—1 个阵元的L 型延迟阵列MWC 结构接收, 当满足以下条件:

1)阵元间距d≤c/fNyq, 阵元个数N≥M+1 ,通道延迟τ≤1/fNyq;

2)fs≥fp≥B;

4.4 复杂度分析

接下来分析基于ESPRIT 的参数估计算法的复杂度.计算矩阵R大约需要N2×Q次复乘, 其中Q为每通道快拍数,N为一个轴上的阵元个数.对R进行奇异值分解大约需要 6 (N -1)3次复乘,对矩阵 (V1+V2+V3) 进 行特征值分解, 需要M3次复乘,M为信源个数.因此, 该算法的时间复杂度大约为O(N3+N2Q+M3) , 所需的空间复杂度约O(N2+M2).

同样, 对基于CP 分解的参数估计算法的复杂度进行分析.交替最小二乘算法的每次循环中都要对3 个矩阵分别进行求解, 其中计算R⊙Ay1大约需 要 (N -1)2×M次 复 乘, 计 算R⊙Ay1的 伪 逆则需要2(N-1)×M2+M3次复乘,计算X(1)与的乘积则需要(N-1)3×M次复乘.假设交替最小二乘算法循环I次, 则CP 分解的时间复杂度大约为O(IM ·N3+INM2+N2Q) , 其大小和迭代循环次数有直接关系, 迭代的次数由设定的阈值和最大循环次数有关.在计算过程中需要存储矩阵X(k)和因子矩阵等, 所需的空间复杂度约为O(N3).可以看出, CP 分解方法的单次循环的计算复杂度与ESPRIT 方法相当, 但由于CP分解方法是求解优化问题, 需要通过不断重复迭代直到达到收敛条件, 而ESPRIT 方法通过一次特征值分解即可求得解析解.因此无论是计算的时间复杂度还是空间复杂度, 基于ESPRIT 的方法均小于CP 分解方法.

表1 不同方法的复杂度对比Table 1.Complexity comparison of different methods.

文献[28]中的基于CS-OMP 的2D-DOA 与载频联合估计方法, 利用双L 型阵列MWC 阵列结构进行信号接收, 用OMP 算法进行参数估计,其所需的时间复杂度大约为O(NM2P2+N3+N2(M+Q)), 其中P为分割的网格个数.空间复杂度大约为O(NP2+N2+NM) , 可见CS-OMP算法复杂度很大程度上由网格的大小决定, 网格越小, 参数估计精度虽然会变高, 但时间复杂度和空间复杂度会急速增长.将3 种方法的复杂度总结如表1 所示.图4 为M= 3,Q= 160, 迭代次数I=50 时3 种算法的时间复杂度, 结果表明ESPRIT方法的复杂度最小, 但是从第5 节仿真实验可以看出CP 分解方法在鲁棒性方面有更好的表现.

图4 不同方法在M = 3, Q = 160, 迭代次数I = 50 时的复杂度对比图Fig.4.Multiplications comparison vs.N with M = 3, Q =160, and I = 50.

5 仿真实验

在本节中, 对本文提出的基于L 型延迟阵列MWC 结构的ESPRIT 算法、CP 分解算法以及CS-OMP 方法[28]进行对比分析, 探究3 种估计方法在不同阵元数、快拍数以及不同信噪比情况下的参数估计表现.实验中使用的信号均根据2 节中定义的窄带信号模型由仿真生成.定义

分别作为载频和DOA 参数估计准确度的评价指标[24], 其中和φˆi分别为载频、方位角和俯仰角的估计值.

设置信号个数M=3 , 奈奎斯特频率fNyq=10 GHz, 带宽B= 150 MHz.伪随机序列p(t) 每周期65 个点, 周期频率fp= 1.1B= 154 MHz.调制载频fi、方位角θi和俯仰角φi均在定义范围内随机选取.设置每个轴的阵元数N=6 , 间距d=0.03 m.低通滤波器的截止频率为fp/2=77 MHz,每通道采样率为fs=fp.

首先对传感器阵元个数对估计效果的影响进行探究.设置信噪比SNR=10 dB, 快拍数Q= 100,每个轴的天线个数从临界最小个数N=M+1=4到N=10 递增, 步进值为1.仿真结果如图5 所示,可以看出3 种方法随着天线个数的增加, 参数估计效果逐渐变好, 这主要因为传感器数量的增加提高了系统对噪声的鲁棒性, 并使其能够处理更多的源信号.同时可以发现在信噪比为10 dB 的情况下,CP 分解方法的估计效果要比ESPRIT 算法好, 因为CP 分解方法可以充分利用多维信息, 直接处理接收数据形成的张量模型, 具有更好的性能.结合复杂度分析的结论可知, 相比另外两种方法,ESPRIT 方法的计算复杂度更低, 运算量较小.而基于双L 型阵列的CS-OMP 方法估计精度受网格大小限制, 因此效果不如提出的两种方法.

图5 不同阵元个数下参数估计效果 (a) 载频估计效果; (b) 方位角估计效果; (c) 俯仰角估计效果Fig.5.Performance of estimated parameters under different N: (a) Performance of carrier frequency; (b) performance of azimuth angle; (c) performance of elevation angle.

图6 不同快拍数Q 下参数估计效果 (a) 载频估计效果; (b) 方位角估计效果; (c) 俯仰角估计效果Fig.6.Performance of estimated parameters under different Q : (a) Performance of carrier frequency; (b) performance of azimuth angle; (c) performance of elevation angle.

图7 不同信噪比下参数估计效果 (a) 载频估计效果; (b) 方位角估计效果; (c) 俯仰角估计效果Fig.7.Performance of estimated parameters under different SNR: (a) Performance of carrier frequency; (b) performance of azimuth angle; (c) performance of elevation angle.

接下来探究快拍个数Q的大小对参数估计精确度的影响.设置信噪比为10 dB, 阵元个数N= 10,Q的值从40 到110, 步进值为10.参数估计结果如图6 所示, 可以看出随着快拍数Q的增大, 各方法的信号参数估计效果变好.所提的阵列结构的总阵元个数为2N— 1=19 个, 文献[28]中的双L 型阵列的总阵元个数为3N— 1=29 个.由于所提方法使用阵元个数较少, 因此在快拍数较少的时候参数估计的表现不如CS-OMP 方法, 但是在快拍数大于80 时, 所提的两种方法的载频和方位角参数估计效果均比CS-OMP 方法好, 这是由于CS-OMP方法的参数估计精度受网格大小的限制.

最后探究系统信噪比大小对各个方法估计精度的影响.设置信噪比从—20 dB 到30 dB 递增, 步进值为5 dB, 快拍数Q= 100, 阵元数N= 10.实验结果如图7 所示, 随着信噪比的增大, 3 个方法的参数估计效果均越来越好.通过对比可以发现,基于ESPRIT 的参数估计方法在信噪比小于20 dB时, 鲁棒性不如CP 分解方法, 这是因为ESPRIT方法需要先计算通道采样值间的协方差函数, 而CP 分解方法可以直接处理多维的信息, 具有更好的性能.由于CS-OMP 方法阵元个数比所提的阵列结构多, 因此在信噪比较低时效果较好.但在信噪大于5 dB 时, CS-OMP 方法的精度受网格限制, 参数估计的误差较大.同时CS-OMP 算法只能重构θi ∈(0°,90°) 内的角度, 提出的两种方法均可以恢复在 (-90°,90°) 之内的角度, 范围较大, 更加适合于实际应用.

从仿真结果和复杂度分析结果可以看出, 提出的基于ESPRIT 的方法相对计算量更小, 实时性更高, 但估计效果比CP 分解方法稍差.提出的CP 分解方法的相对鲁棒性更好, 但计算量大约是ESPRIT 方法的IM倍.CS-OMP 算法在信噪比较低时对2 D-DOA 的估计效果较好, 但信噪比较高时精度受限于网格大小, 计算复杂度随网格的增加会呈指数增长.同时提出的方法可以估计的方位角范围为 (-90°,90°) , 而CS-OMP 算法可估计范围仅为 ( 0°,90°) , 因此提出的算法更加具有实用性.

6 结 论

本文提出了一种L 型延迟阵列MWC 的接收结构, 通过增加延迟结构直接估计信号载频, 进而估计2D-DOA 参数, 无需额外的参数配对, 避免了配对引入的误差和复杂度的提升.接着在该阵列结构特点的基础上分别提出了两种参数估计方法, 一种基于ESPRIT 方法, 另一种基于CP 分解技术,针对两种方法进行了复杂度分析和仿真, 探究了阵元间距、快拍数以及信噪比对方法估计精度的影响, 并对提出的两种方法以及文献[28]中的基于网格的CS-OMP 方法进行对比分析, 从实验结果中可以总结出所提方法均可从欠奈奎斯特样本中估计信号的载频和二维DOA 参数.基于ESPRIT方法的计算量较小, 但不足之处是低信噪比条件下表现欠佳, 适用于主动雷达等要求实时性的场合;而基于CP 分解技术的方法鲁棒性更好, 适用于低信噪比的场合, 例如被动雷达、遥感等被动感知的目标探测领域, 但所需计算量较大, 实时处理困难.本文研究了空间电磁波信号在欠采样条件下的载频和二维到达角估计, 如果在已知信号的某些特征的情况下, 还可以进一步改进来减小算法复杂度或提升抗噪性等.未来在减小算法复杂度方面可以结合信号已知信息来合理的选取初始值, 或者利用梯度下降等方法加快收敛速度, 减少循环次数; 在增强参数估计精度方面可以利用先验信息设计合理的去噪算法; 也可以利用互质阵列、嵌套阵列等稀疏阵列来减少接收阵列的阵元个数, 降低硬件复杂度和成本.