基于酉重构子空间的一维DOA估计

金彦亮，闾儒坤，汪小勇，郑国莘

1.上海大学通信与信息工程学院，上海200444

2.上海大学上海先进通信与数据科学研究院，上海200444

3.卡斯柯信号有限公司，上海200070

4.上海轨道交通无人驾驶列控系统工程技术研究中心，上海200434

通信感知一体化是一种基于软硬件资源共享，协同实现通信和感知功能的新型信号处理技术[1]，多输入多输出（multiple input multiple output,MIMO）技术在毫米波通信和雷达上的运用给通信感知一体化系统带来了全新的突破[2]。无线通信系统的频段向毫米波等高频段的发展，这将与传统感知系统的频段产生重叠。因此，在同一频段实现通信与感知，避免频谱之间的干扰，提升频谱利用率，成为通信感知一体化系统的发展目标[3]。

通信感知一体化场景下的感知任务通常包括对目标进行测距、测角、测速等。信号波达方向（direction of arrival,DOA）估计是通信感知一体化系统目标参数估计中的重要组成部分，准确的DOA 估计对建立大规模MIMO 通信感知系统的信道模型是至关重要的[4]，同时也有助于实现感知功能中的高精度定位[5-6]。常用的角度参数估计方法以基于信号子空间与噪声子空间正交性的多重信号分类（multiple signal classification,MUSIC）算法[7]和基于信号子空间旋转不变性的旋转不变子空间[8]（estimating signal parameters via rotational invariance techniques,ESPRIT）算法为代表，其中MUSIC 算法突破了瑞利限在空间分辨率上的限制[9]。然而仅利用噪声子空间信息的传统MUSIC 算法在低信噪比、小节拍和低阵元数目等恶劣情况下有较大的估计偏差。为了充分利用子空间信息，提高系统估计精度，降低算法复杂度和运行时间，文献[10] 采用一种应用于MUSIC 算法的酉变换，降低了计算复杂度。文献[11] 提出双平行线阵（double parallel linear array,DPLA）的DOA 算法，避开了协方差矩阵特征值分解，降低了计算复杂度，然而该算法只通过线性变换得到噪声子空间，在低阵元数目条件下，表现出较大的估计偏差。文献[12] 利用新特征值矩阵对应的特征向量得到噪声子空间，虽然其充分利用了噪声数据，但是忽视了信号数据本身的分辨能力。文献[13]通过把加权伪噪声子空间投影引入了传统MUSCI 算法，确保子空间投影算法高分辨能力的同时提高了算法对小节拍和强弱信号源情况下的估计稳定性。文献[14] 利用信号子空间的抗噪性，给出了信号子空间投影加权DOA 算法，该算法较于MUSIC 算法具有更好的稳健性和空间分辨率。文献[15] 提出了子空间降维MUSIC（subspace dimension reduced MUSIC,SDR-MUSIC），通过降低噪声子空间的维数来减少运行时间。文献[16] 提出了基于量子行为粒子群优化（quantum-behaved particle swarm optimization,QPSO）的信号子空间重构增强检测算法，对低信噪比和小节拍数有很强的适应性。

以上工作虽然在一定程度上改善了算法的性能，但都限于改变信号或者噪声子空间的投影比例，没有利用信号源之间的特征信息。为了充分利用信号特征向量信息，本文提出了一种基于酉重构子空间的改进算法。该算法根据文献[17] 中提出的酉变换技术将阵元接收数据实数化，在后续的信号处理过程中将复数域的数据转化到实数域中，从而显著降低计算复杂度。协方差矩阵特征分解得到的信号子空间是由与信号源数目相同的较大特征值所对应的特征向量构成，通过逐次剔除较大信号特征向量得到重构子空间，避免了正交分量较大的信号源对其余信号的影响，同时结合校正矩阵得到新的空间谱函数，最后与信号子空间投影（signal subspace projection,SSP）算法[18]联合估计。

仿真结果表明，本文算法相较于MUSIC 算法和SSP 算法在低信噪比、低阵元数目和小节拍数等恶劣条件下，空间谱峰更尖锐。在相同的仿真条件下，统计性能达到检测误差不超过1◦时，本文算法所需信噪比和节拍数较MUSIC 算法分别减少7 dB 和80 次。在入射方位间隔小且信号源相干的情况下，本文算法能够在信源方位间隔仅为5.5◦的情况下达到100% 的分辨率。

1 算法理论依据

1.1 阵列信号处理模型

均匀线性阵列（uniform linear array,ULA）模型如图1 所示，假设空间有M个天线阵元组成ULA，阵元间距d相等且不大于信源半波长λ，即d≤λ/2，其中，λ=2πc/f，f为信源中心频率，c为光速。

现以阵元1 为参考阵元，假设有D个不相关的远场窄带信号源入射到ULA 上，信号方向为θi(i=1,2,···,D)。假设接收到的加性噪声为平稳、零均值的高斯白噪声，方差为且信号与噪声相互独立，则第k个阵元在第n次节拍输出的信号矢量形式xk(n) 为

式中：k=1,2,···,M；n=1,2,···,L。ak(θi) 为阵元k对第i个信号的方向矢量，由此可得窄带远场信号的数学模型为

式中：M为阵元数目；X(n) 为M×1 维阵列输出矢量；S(n) 为入射信号的D×1 维信号矢量；A ∈CM×D为M×D维阵列流型矩阵；C 为复数域空间；a(θi) 为M×1 维方向矢量；N(n) 为M×1 维噪声向量。

1.2 传统MUSIC 算法原理

对阵列接收数据X做自相关得到协方差矩阵RX

式中：(·)H是矩阵的共轭转置。由于信号与噪声相互独立，数据协方差矩阵可以分解成信号和噪声两部分[19]

式中：ARSAH是信号部分；RS是信号协方差矩阵；RN是噪声协方差矩阵。由于接收到的加性噪声为平稳、零均值的高斯空间白噪声，方差为，则

对协方差矩阵RX进行特征分解得

式中：U为特征向量矩阵；对角阵Σ=diag[λ1,λ2,···,λD,···,λM] 由特征值组成且特征值满足λ1≥λ2≥··· ≥λD≥λD+1=λD+2=···=λM=；M为阵元数目。由于阵列流型矩阵A是范德蒙矩阵，当D个互不相关的远场窄带信号从不同方向入射到ULA，矩阵A的列向量线性无关，并且信号协方差矩阵RS为非奇异阵，则

由于信号协方差矩阵RS为正定矩阵，因此矩阵ARSAH的特征值为正值，即共有D个正特征值。因为>0，所以矩阵RX的M个特征值均为正值。因此，较大的D个特征值是信号特征值，其余是M-D重的噪声特征值。设M-D重的特征值对应的噪声特征向量为ei(i=D+1,···,M)，则

将式(5) 右乘ei得

又λmin=，所以

由于阵列流型矩阵A是范德蒙矩阵，信号协方差矩阵RS是正定矩阵，因此

即信号协方差矩阵RS的噪声特征向量ei(i=D+1,···,M)与阵列流型矩阵A的列向量相互正交。设M-D个噪声特征值对应特征向量构成噪声子空间为UN=[vD+1,vD+2,···,vM]，由D个信号特征值对应特征向量构成信号子空间为US=[v1,v2,···,vD]，则协方差矩阵RX改写为

由文献[20] 可知US与UN分别表示与信号子空间和噪声子空间对应的正交基向量。并且由式(11) 可知，矩阵A中的方向矢量也与噪声子空间正交即

则传统MUSCI 算法的谱估计公式为

对式(14) 进行谱峰搜索，在理想情况下，空间谱图像会在对应的入射角方位上形成D个尖锐谱峰。

2 基于酉变换的DOA 估计算法

式(3) 得到的阵列信号协方差矩阵是一个复数矩阵，对复数协方差矩阵特征分解的运算相比于实数运算，在节拍数和阵元数目较多的情况下，运算量呈指数增长。本文利用酉变换将复矩阵及其特征向量转换为实矩阵，在不牺牲精度的情况下降低计算成本[21]，并保留原始入射信号的来波方向信息。定义一个复矩阵Υ ∈CM×D且满足以下等式

则称Υ是中心Hermitian 矩阵。式中：Jp是p阶的反对角线上为1，其他元素均为0 的交换方阵；(·)∗代表复共轭运算。

定义酉矩阵Q，表达式为

则中心Hermitian 矩阵Υ可以通过酉变换转化为实数矩阵，具体表达式为

在本节中，先对协方差矩阵RX进行空间平滑[22]操作，使之成为中心Hermitian 矩阵。已知协方差相关矩阵RX，那么中心Hermitian 矩阵Rcentre表示如下

当该阵列由偶数M个天线阵元组成时，对应的酉变换复数阵Q为

式中：I为阶单位矩阵；J为阶交换矩阵。利用式(17) 将Rcentre实数化得到矩阵Rreal表示如下

式中：(·)T代表转置。对实数矩阵Rreal进行奇异值分解得

式中：Ureal=[w1,w2,···,wM] 为特征向量空间矩阵；Σreal为特征值构成的对角阵，ΣNreal=diag[λ1-real,λ2-real,···,λM-real] 且λ1-real≥λ2-real≥···≥λM-real，且酉变换后方向矢量表示为areal(θ)=QTa(θ),i=(1,2,···,D)，则对应的谱估计表达式为

本节通过构造酉变换矩阵，对协方差矩阵进行变换，在实数域进行特征值分解。同时，环境中的相干信号会导致信号协方差矩阵秩的亏损，进而导致方向矢量与噪声子空间不完全正交，严重影响了MUSIC 算法的分辨力。式(18) 通过空间平滑，恢复信号协方差矩阵的秩[23]，提高了本文算法在处理相干信号时的鲁棒性。

3 基于重构子空间的DOA 估计算法

传统MUSIC 算法利用噪声子空间与方向矢量正交特性来估计信号入射角度。实际接收数据经过数字接收机采样得到，通常由最大似然估计值ˆR代替协方差矩阵即

式中：L为采样节拍数。在低信噪比、低阵元数目和小节拍数等恶劣情况下，采样协方差矩阵与理想协方差矩阵偏差会增大，MUSIC 算法对多目标的检测能力下降。尤其是当相邻信号以较小的角度间隔入射时，空间谱峰会出现互相干扰，甚至混叠在一起的情况，严重降低了MUSIC 算法的分辨率。因此本文基于酉变换得到的实数协方差矩阵式(20)，结合SSP 算法[18]，改变方向矢量在信号子空间上的投影分量的大小，然后通过依次剔除较大信号特征向量，并使用对应的特征值校正矩阵重构子空间，削弱相邻目标之间的混叠现象，提高算法准确率和分辨率。

3.1 信号子空间投影

由文献[18] 可知，在低信噪比、低阵元数目等条件下，通过采样协方差矩阵特征分解得到的信号子空间受噪声干扰的影响小，对小节拍数引起的误差较小。本节利用酉变换后的信号特征值的倒数对角矩阵校正信号向量子空间，弱化较大的信号特征向量投影分量后与式(22)联立得到新的空间谱函数表达式为

式中：ΣSreal=diag[λ1-real,λ2-real,···,λD-real]，USreal=[w1,w2,···,wD]。

从式(24) 可以看出，该算法是遍历方向矢量在酉变换后的信号子空间投影与原空间谱的叠加，有效利用了信号子空间的信息，提高了算法在非理想条件下的分辨率。

3.2 重构子空间算法

重构子空间算法是在酉变换后的实数域特征向量的基础上，根据特征值大小，依次剔除较大信号特征值对应的特征向量，削弱较大信号特征分量对其余信号的干扰，得到重构子空间，并根据特征值大小对重构子空间进行校正，得到一组新的空间估计谱，最后与SSP 算法联合估计，算法的基本步骤如下。

根据式(21) 对实数采样协方差矩阵特征分解得到特征向量空间矩阵和对应的特征值对角阵

依次对式(25) 和(26) 剔除最大信号特征向量和信号特征值，第一次递减最大信号特征值λ1-real以及对应的特征向量w1得到重构子空间以及校正矩阵

第二次递减最大信号特征值λ2-real对应的特征向量w2得到子空间以及校正矩阵

以此类推，第D-1 次递减最大信号特征值λD-1-real对应的特征向量wD-1得到子空间以及校正矩阵

将重构子空间重构谱函数与式(22) 相结合，得到新的空间估计谱公式为

最后与式(24) SSP 算法进行联合估计得到一种新的空间谱估计算法，其谱函数为

由式(31) 可以看出，该算法增加了方向矢量对较小信号特征向量子空间的扫描次数，在削弱较大信号特征分量的同时重构子空间，并利用校正矩阵加强较小的信号特征分量，有效地避免大信号特征分量对邻近信号的干扰。

3.3 计算复杂度分析

与MUSIC 算法相比，本文算法增加了重构子空间矩阵Φ和校正矩阵Z的计算，并且在式(31) 中构造空间谱函数时复用了重构子空间。但是酉变换技术的运用能够将复数域运算转化为同等维数矩阵下的实数域运算，使得本文算法在性能高于其他算法的同时，计算量不会大幅升高。

算法计算复杂度一般通过矩阵间的复乘次数来衡量[24]。假设M为阵元数目，L为抽样节拍数，D为信号源数目，J为谱峰搜索步长。通过分析可知，MUSIC 及SSP 算法对阵列输出矢量进行协方差的计算量与阵列天线规模和抽样节拍数有关，其所需的计算复杂度为O(M2L)。直接通过对协方差矩阵RX进行特征向量分解得到信号和噪声子空间所需的计算复杂度为O(M3)。MUSIC 算法谱峰搜索阶段[25]所需的计算量为O((180/J)(M+1)(M-D))，SSP算法谱峰搜索阶段所需的计算量为O{(180/J)[(M+1)(M+MD2-D)]}。本文算法通过式(18) 和(20) 构造实数协方差矩阵，增加的计算量[26]为O(M3+(2M)2)。入射信号源个数D对本文算法峰搜索阶段计算量的影响较大，为O{(180/J)[(M+1)(M+MD2-D)+M(M+1)(M-i)2]}。

另外酉变换技术的运用，使得本文算法在特征值分解和谱峰搜索阶段所需的计算量与复数情况相比，计算量可以降低到原来的25%，并且本文算法在达到较高的估计精度时，只需要较低节拍数和阵元数目，因此在对阵列输出矢量进行协方差运算的计算量更少。

表1 3 种算法计算复杂度对比Table 1 Comparison for computing complexity of three DOA estimation algorithms

4 仿真实验与结果

阵元数目、节拍数和信噪比是影响DOA 算法分辨率的重要因素，为了进一步验证所提算法在低信噪比、低阵元数目和小节拍数情况下，对检测入射间隔小的信号的有效性。本文在非相干信号的环境下，对所提算法进行多次Monte Carlo 实验，并且在相同实验条件下与MUSIC 算法和SSP 算法进行多方面对比。最后，在环境中存在相干信号的情况下，本文对3种算法进行实验仿真，以验证所提算法在处理相干信号情况时的有效性。

4.1 算法空间谱图像对比

通过仿真实验对比了3 种算法的分辨能力。假设入射信号为3 个不相干且波长为λ的远场窄带信号，分别以-3.0◦、2.0◦和50.0◦的方向入射到阵元数目为8、阵元间距为λ/2 的均匀线阵上。信号与噪声相互独立，且噪声为加性高斯白噪声，信噪比为0 dB。节拍数为50。图2(a) 为3 种算法的空间谱函数单次实验对比图。由图2(a) 可以看出，在低信噪比、低阵元数目和小节拍数条件下，MUSIC 算法和SSP 算法可以检测出入射角为50◦的来波信号，却无法分辨出两个相邻入射信号，两种算法在该场景下失效。而本文算法在-3◦和2◦处精确检测出两个目标，并且有较为尖锐的谱峰，初步验证了本文算法在信号入射间隔小情况下的分辨性能优于其他两种算法。

图2 算法空间谱及分辨概率变化关系图Figure 2 Comparison of spatial spectrum and probability of resolution between algorithms

图2(b) 展示了分辨概率与方位角度间隔的关系曲线，进一步说明了本文算法对入射间隔较小的信号具有更好的分辨概率。假设2 个远场窄带非相干信源，入射方向分别为θ1和θ2，其中θ1是随机角度，θ2=θ1+∆且∆从0◦∼10.0◦以0.5◦为间隔增长，且对每一个间隔进行500 次Monte Carlo 实验。对任一次Monte Carlo 试验，如果估计出的两信号角度和满足则称两信号正确分辨[27]，分辨概率指的是正确分辨次数占实验总数的百分比。由图2(b) 可以看出，3 种算法的分辨成功率随信源方位间隔的增大而提高，且本文算法的分辨成功率始终高于MUSIC 算法和SSP 算法，本文算法能够在信源方位间隔仅为4◦的情况下达到100% 的分辨概率。

下面将本文算法、MUSIC 算法和SSP 算法仿真统计结果随信噪比、节拍数、阵元数目变化进行多方面量化对比。在以下仿真中，均采用ULA 阵列，噪声为加性高斯白噪声。为体现本文算法在低信噪比、低阵元数目和小节拍数等非理想条件下对入射间隔小的信号分辨率的优势，本文采用均方根误差（root mean squared error,RMSE）来衡量3 种算法在检测入射角度相近的两个不相干信号的估计误差。基于Monte Carlo 实验的D个入射信号角度估计值的均方根误差定义为

式中：ρ为Monte Carlo 实验的次数；为第ρ次实验中第D个信号入射角度的估计值；θD为该入射信号的实际角度。

4.2 均方根误差随信噪比变化对比

假设信号源为两个不相干且波长为λ的远场窄带信号，分别以10◦和14◦的方向入射到阵元数目为10、阵元间距为λ/2 的均匀线阵上，节拍数为40。设置信噪比从-5 dB 按1 dB的步长增加到15 dB，对3 种算法进行5 000 次独立仿真实验并比较RMSE，如图3 所示。

图3 均方根误差随信噪比变化曲线Figure 3 RMSE contrast curves as the function of SNR

由图3 可以看出，随着信噪比的增大，3 种算法的RMSE 均随之减小，其中本文算法的RMSE 下降速度快，误差最小。要达到检测误差不超过1◦，本文算法仅需信噪比5 dB，比SSP 算法减少5 dB，比传统MUSIC 算法减少7 dB。这说明本文算法在低信噪比的情况下有更好的估计精度。

4.3 均方根误差随节拍数、阵元数目变化对比

入射角度不变，在信噪比为5 dB、阵元数目为10 的条件下，设置节拍数从10 按步长为5 增加到100，对3 种算法进行5 000 次独立仿真实验，并比较RMSE，如图4 所示。

图4 均方根误差随节拍数变化曲线Figure 4 RMSE contrast curves as the function of snapshots

在信噪比5 dB 条件下，随着节拍数的增大，3 种算法的RMSE 均随之减小。由图4 可以看出，在系统RMSE 均为1◦时，本文算法节拍数比SSP 算法降低80 次，大幅度减少了系统估计所需要的采样次数。随着节拍数的增加，本文算法与SSP 算法的性能差距逐渐减小，但始终比另外两种算法RMSE 要小。

在信噪比为5 dB、节拍数为40 的条件下，设置阵元数目从6 按步长为2 增加到20，进行5 000 次独立仿真实验，并比较RMSE，如图5 所示。

图5 均方根误差随阵元数目变化曲线Figure 5 RMSE contrast curves as the function of array number

由图5 可以看出，在阵元数目变化的情况下，本文算法的RMSE 始终小于其他两种算法，特别是当阵元数目较少时，本文算法的RMSE 明显小于传统MUSIC 算法和SSP 算法。

4.4 相干信号空间谱对比

本实验通过仿真对比本文算法和另外两种算法在处理相干信号方面的性能。假设3 个远场窄带信号，分别从-10◦、20◦和50◦的方向入射到阵元数目为8、阵元间距为λ/2 的均匀线阵上。其中，20◦和50◦的入射信号是相干信号，信噪比为5 dB，节拍数为50。图6(a) 为3 种算法处理相干信号时空间谱函数单次实验对比图，从图6(a) 中可以看出，MUSIC 算法和SSP 算法谱函数只在-10◦方向上产生了一个尖锐的谱峰，无法对相干信号有效估计。而本文算法在结合空间平滑和酉变换技术之后，能够得到正确的谱峰，估计精度高，初步验证了本文算法适用于相干入射信号的角度估计。

图6 相干信号算法空间谱及分辨概率变化关系图Figure 6 Comparison of coherent signals spatial spectrum and probability of resolution between algorithms

假设仿真实验中为一阵元数目为8、阵元间距为λ/2 的ULA，有2 个远场窄带相干信号分别以波达角度θ1和θ2入射到此阵列，其中θ1是随机角度，θ2=θ1+∆且∆从0◦以0.5◦为步长间隔增长至10.0◦。每一方位间隔均进行500 次Monte Carlo 实验，对每一算法进行仿真并统计在不同方位间隔下的对相干信号的分辨概率，结果如图6(b) 所示，可以看出，在信号相干时，本文算法的分辨成功率随信源方位间隔的增大而提高，其他两种算法却无法对相干信号进行有效的DOA 估计。

5 结语

传统MUSIC 算法在低信噪比、低阵元数目和小节拍数等非理想环境下对入射间隔较小的信号检测能力严重下降，难以适应实际工程环境。SSP 算法虽然减小了较大信号特征向量的投影分量，但是无法在邻近信号处形成尖锐的谱峰。针对以上问题，本文通过特征值校正后的重构子空间构造新的空间谱函数，减少空间谱峰混叠的现象，并与SSP 算法联合估计，提高算法的分辨率。在阵列输出矢量数据处理阶段，本文采用空间平滑和酉变换技术，降低了算法的计算量，并且针对性地解决了MUSIC 算法和SSP 算法在相干信号环境下失效的问题。

目标方位参数估计算法的性能与通信感知一体化系统的时效性、精确性等紧密相关，直接影响通信感知一体化服务质量。并且通信感知场景中信道环境复杂，大规模阵列天线的运用也使得传统算法的计算量呈指数级增加，本文方法能够有效避免因信号相干带来的DOA估计性能下降的影响，具有较高的目标方位估计精度和较低的计算复杂度。通信感知一体化场景下目标感知参数估计有助于节省通信系统资源开销，提高系统通信感知效率，因此未来的工作集中在将DOA 估计算法的应用到轨道交通[28]、5G 车联网[29]、无人机通信与传感[30]等领域中。此外，低复杂度参数估计、实时高精度参数估计、多感知节点联合参数估计等都是通信感知一体化参数估计技术的未来研究重点。