CN113221059B - 无需构造协方差矩阵的快速共轭梯度测向算法 - Google Patents

Abstract

本发明属于阵列信号处理技术领域，具体的说是一种实现了在不损失精度的情况下，显著降低了计算复杂度的无需构造协方差矩阵的快速共轭梯度测向算法，其特征在于，接收辐射源信号后，利用最小二乘法改进维纳霍夫方程，获得最小二乘维纳霍夫方程，最小二乘维纳霍夫方程为正态方程，使用正态方程构造仅与接收数据有关的参考信号，设置迭代初始值并进行共轭梯度迭代后，根据信源数估计出信号子空间并获得信号波达方向。

Description

无需构造协方差矩阵的快速共轭梯度测向算法

技术领域：

本发明属于阵列信号处理技术领域，具体的说是一种实现了在不损失精度的情况下，显著降低了计算复杂度的无需构造协方差矩阵的快速共轭梯度测向算法。

背景技术：

信号的波达方向估计是雷达、声呐、无线通信和无源定位等应用中经常遇到的重要研究课题。以多重信号分类和旋转不变子空间为代表的子空间类算法的提出，实现了传统空间谱估计向超分辨测角的飞跃，但MUSIC算法庞大的计算量和ESPRIT算法较低的估计精度阻碍了超分辨算法的工程化进度。针对此问题，促成了基于Krylov子空间的共轭梯度(Conjugate Gradient，CG)算法的诞生。

CG算法是最经典的基于Krylov子空间的波达方向(Direction of Arrival,DOA)估计方法之一。与传统的特征子空间类算法不同的是，CG算法利用Krylov子空间代替信号子空间，由梯度迭代实现了快速DOA估计的有效算法。由于CG算法在迭代初始值设置与迭代过程中都要计算阵列协方差矩阵(Array Covariance Matrix，ACM)，计算ACM的计算复杂度约为O(M²N)，在梯度迭代过程中若避免ACM的计算，则计算复杂度将会大大降低，在大快拍数时尤其明显。更重要的是，大多数现有的CG算法与MUSCI算法相比，在降低计算复杂度的同时还保持良好的估计性能。因此，相比传统方法，应该选择基于Krylov子空间的CG技术。

上述算法丰富了CG的理论内涵，但在实际工程应用中需要借助迭代估计出信号子空间。众所周知，迭代初始值的设置影响着算法的效率和准确性，不合适的处置可能导致迭代无法收敛或收敛到错误解。

发明内容：

本发明针对现有计算ACM引起CG算法高复杂度的问题，提出一种通过改进维纳霍夫方程将ACM利用接收数据矩阵代替来减少CG算法的计算复杂度，实现了在不损失精度的情况下，显著降低了计算复杂度的无需构造协方差矩阵的快速共轭梯度测向算法。

本发明可以通过以下措施达到：

一种无需构造协方差矩阵的快速共轭梯度测向算法，其特征在于，接收辐射源信号后，利用最小二乘法改进维纳霍夫方程，获得最小二乘维纳霍夫方程，最小二乘维纳霍夫方程为正态方程，使用正态方程构造仅与接收数据有关的参考信号，设置迭代初始值并进行共轭梯度迭代后，根据信源数估计出信号子空间并获得信号波达方向。

本发明所述利用最小二乘法改进维纳霍夫方程具体包括以下步骤：

步骤1-1：维纳霍夫方程为RW＝r_xd，化简维纳霍夫方程，R和r_xd分别表示为：

其中X为接收数据信号矩阵，有X＝[X(0),X(1),…,X(M-1)]，d为参考信号，有d＝[d(0),d(1),…,d(M-1)]；

将以上式子带入到维纳霍夫方程中，有

化简后得到X^HW＝d^H；

步骤1-2：定义最小二乘维纳霍夫方程，为r+X^HW＝d^H，其中，r为最小二乘残差矢量，最小二乘维纳霍夫方程为正态方程；

步骤1-3：为避免计算阵列协方差矩阵R，并且根据最小二乘维纳霍夫方程得到改进的过渡矢量，为v＝X^Hp；

步骤1-4：定义中间变量正态方程残差为q＝Xr，其中q只与阵列接收数据与参考信号有关。

本发明所述构造仅与接收数据有关的参考信号具体包括以下步骤：

步骤2-1：定义参考信号为

其中，a(θ)为搜索导向矢量，θ∈[-90°,90°]，d为参考信号；

步骤2-2：由于第k个阵元的接收信号数据矩阵表示为

则当θ＝θ_i，i∈{1，…，K}时，正态方程残差初始值q_cg,0表示为

可知，正态方程残差初始值q_cg,0为信号导向矢量的线性组合，则q_cg,0张成的K维的Krylov子空间为信号子空间的子集，表示为

span{Q_cg,K}≡span{A(θ)}，

其中

Q_cg,K＝{q_cg,0,q_cg,1,…,q_cg,K-1}；

步骤2-3：当θ≠θ_i，i∈{1，…，K}，正态方程残差初始值q_cg,0表示为：

此时正态方程残差矢量也为导向矢量的线性组合，是由K个信号矢量和搜索导向矢量组合的，即此时的正态方程残差矢量张成的K+1维Krylov子空间也为信号导向矢量张成的子空间的子集，表示为

span{Q_cg,K+1}≡span{A(θ),a(θ)}，

其中

Q_cg,K+1＝{q_cg,0,q_cg,1,…,q_cg,K-1,q_cg,K}。

本发明所述设置迭代初始值并进行共轭梯度迭代具体包括以下步骤：

步骤3-1：共轭梯度初始值设置中，初始权矢量为：w₀＝0，搜索向量初始值：p₁＝-r₁＝-d,其中d为参考信号；

步骤3-2：初始值设置中最小二乘残差矢量初始值为：r₀＝d^H-X^HW，搜索矢量初始值：p₁＝q₀，γ₀＝||q₀||²，正态方程残差矢量初始值为：q₀＝Xr₀；

步骤3-3：当迭代次数为k＝1,2,…,M，得到中间过渡变量：v_k＝X^Hp_k，确定步长：α_k＝γ_i-1/||v_k||²，得到权矢量：w_k＝w_k-1+α_kp_k，利用最小二乘残差：r_k＝r_k-1-α_kv_k，得到正态方程残差：q_k＝Xr_k，确定共轭方向向量：

本发明所述根据信源数估计出信号子空间并获得信号波达方向具体包括以下步骤：

步骤4-1：已知信源数为K，仅需要迭代K次就得到等价于信号子空间的Krylov子空间为：W_cg,K＝{w_cg,0,w_cg,1,…,w_cg,K-1}，MUSIC算法谱峰搜索函数为：

其中a(θ)为搜索导向矢量，θ∈(-90°,90°)，U_N为噪声子空间；

步骤4-2：信源数为先验已知的，则迭代次数选择为k＝1,2,…,K，新算法谱峰搜索函数为：

其中I_M为单位矩阵，W_cg,K为等价于信号子空间的Krylov子空间。

本发明接收辐射源信号表示如下：假设L个相互独立的阵元，以等间距d组成均匀线阵(ULA)，考虑空间中存在K个远场窄带信号入射到阵列。其中，假设K先验已知，d满足d≤λ2，以避免相位模糊，λ为窄带信号的波长，阵列接收辐射源信号为：

其中，A(θ)为L×K维的阵列流型矩阵，S(t)为K×1维的入射信号矢量，N(t)为L×1维的加性高斯白噪声矢量，α(θ)为A(θ)的列向量，可以表示为：

本发明实现了通过改进维纳霍夫方程，利用接收数据矩阵代替ACM，避免了设置迭代初始值与迭代过程中需要计算ACM而引起的高复杂度问题，为波达方向的工程化实现提供了技术支持。

附图说明：

附图1为本发明的流程图。

附图2为本发明和MUSIC算法、CG算法的谱峰搜索对比图，其中M＝10,SNR＝10dB，N＝100，K＝2，θ₁＝20°，θ₂＝40°。

附图3为RMSE随输入SNR的变化曲线图，其中M＝10，信噪比以5dB为步长，从-20dB变化到0dB,N＝200,K＝2,θ₁＝10°,θ₂＝30°。

附图4为本发明和MUSIC算法、CG算法、RV-CG算法的计算效率对比图，其中M从10到70，步长为10，SNR＝10dB,N＝200,K＝2,θ₁＝10°,θ₂＝30°。

具体实施方式：

本发明针对初始值与迭代过程都需要计算ACM引起算法高复杂度的问题，提出一种无需构造协方差矩阵DOA估计器，通过改进维纳霍夫方程将ACM利用接收数据矩阵代替来减少CG算法的计算复杂度，理论分析和试验结果表明，本发明在不损失精度的情况下，显著降低了计算复杂度，还可以有效避免迭代过程中的错误解，从而为CG的实际工程化提供理论参考，下面结合附图和实施例，进一步介绍本发明的技术方案。

如附图1所示，本发明是一种基于Krylov子空间的无需构造协方差矩阵的共轭梯度新方法，具体步骤如下：

第一步，利用天线阵列接收辐射源信号，所述第一步包括以下步骤：(1)假设L个相互独立的阵元，以等间距d组成均匀线阵(ULA)，考虑空间中存在K个远场窄带信号入射到阵列。其中，假设K先验已知，d满足d≤λ/2，以避免相位模糊，λ为窄带信号的波长。阵列接收辐射源信号为：

其中，A(θ)为L×K维的阵列流型矩阵，S(t)为K×1维的入射信号矢量，N(t)为L×1维的加性高斯白噪声矢量，α(θ)为A(θ)的列向量，可以表示为：

第二步，利用最小二乘法改进维纳霍夫方程，所述第二步包括以下步骤：

(1)维纳霍夫方程为RW＝r_xd，经典共轭梯度求解过程为：

首先设置初始权矢量为w₀＝0，搜索向量初始值为p₁＝-r₁＝-d，当k＝1,2…,M时，确定步长为:

过渡矢量为：v＝Rp_k，更新权矢量：w_k＝w_k-1+α_kp_k，得到梯度向量：r_k+1＝r_k+α_kRd_k，确定共轭方向向量：

搜索方向向量p_k是R共轭的，也就是说对所有的i≠j，都有

记矩阵P_k＝[p₁,p₂,…,p_k]是搜索方向向量组成的M×k维矩阵，则P_k是列满秩的。以上分析得到共轭梯度算法最多M步迭代就可以收敛到真实解，第k+1次梯度方向与P_k张成的k维空间正交：

经过k次迭代后，负梯度方向向量r₁，r₂，…，r_k相互正交。搜索方向p₁,p₂,…,p_k张成的空间、梯度方向张成的空间、矩阵R与初始梯度向量r₁形成的Krylov子空间相等：

由此得到信号子空间或噪声子空间，从而进行DOA估计。

(2)化简维纳霍夫方程，R和r_xd可以分别表示为：

其中X为接收数据信号矩阵，有X＝[X(0),X(1),…,X(M-1)]，d为参考信号，有d＝[d(0),d(1),…,d(M-1)]。

因此，将以上式子带入到维纳霍夫方程中，有

化简后得到X^HW＝d^H。

(3)此时的维纳霍夫方程中X并不满足厄米共轭矩阵结构，所以为了满足在Krylov子空间内对共轭梯度迭代的要求，定义最小二乘维纳霍夫方程，为r+X^HW＝d^H，其中，r为最小二乘残差矢量，最小二乘维纳霍夫方程为正态方程。

(4)在传统的共轭梯度算法中，过渡矢量定义为v＝Rp，

其中p为搜索矢量，为避免计算阵列协方差矩阵R，并且根据最小二乘维纳霍夫方程得到改进的过渡矢量，为v＝X^Hp。

(5)定义中间变量正态方程残差为q＝Xr，其中q只与阵列接收数据与参考信号有关。

第三步，构造仅与接收数据有关的参考信号，所述第三步包括以下步骤：

(1)定义参考信号为

其中，a(θ)为搜索导向矢量，θ∈[-90°,90°]，d为参考信号；

(2)由于第k个阵元的接收信号数据矩阵表示为

则当θ＝θ_i，i∈{1，…，K}时，正态方程残差初始值q_cg,0表示为

，从上式中可以看出正态方程残差初始值q_cg,0为信号导向矢量的线性组合，则q_cg,0张成的K维的Krylov子空间为信号子空间的子集，表示为

span{Q_cg,K}≡span{A(θ)}，其中Q_cg,K＝{q_cg,0,q_cg,1,…,q_cg,K-1}。

(3)当θ≠θ_i，i∈{1，…，K}，正态方程残差初始值q_cg,0表示为：

上式可以看出此时正态方程残差矢量也为导向矢量的线性组合，是由K个信号矢量和搜索导向矢量组合的，即此时的正态方程残差矢量张成的K+1维Krylov子空间也为信号导向矢量张成的子空间的子集，表示为

span{Q_cg,K+1}≡span{A(θ),a(θ)}，其中Q_cg,K+1＝{q_cg,0,q_cg,1,…,q_cg,K-1,q_cg,K}。

第四步，设置迭代初始值并进行共轭梯度迭代，所述第四步包括以下步骤：

(1)经典共轭梯度初始值设置中初始权矢量为：

w₀＝0，搜索向量初始值：p₁＝-r₁＝-d,其中d为参考信号。

(2)新算法中的初始值设置中最小二乘残差矢量初始值为：r₀＝d^H-X^HW，搜索矢量初始值：p₁＝q₀，γ₀＝||q₀||²，正态方程残差矢量初始值为：q₀＝Xr₀。

(3)当迭代次数为k＝1,2,…,M，得到中间过渡变量：v_k＝X^Hp_k，确定步长：α_k＝γ_i-1/||v_k||²，得到权矢量：w_k＝w_k-1+α_kp_k，利用最小二乘残差：r_k＝r_k-1-α_kv_k，得到正态方程残差：q_k＝Xr_k，确定共轭方向向量：

第五步，根据信源数估计出信号子空间并获得DOA，所述第五步包括以下步骤：

(1)已知信源数为K，仅需要迭代K次就得到等价于信号子空间的Krylov子空间为：W_cg,K＝{w_cg,0,w_cg,1,…,w_cg,K-1}。

(2)MUSIC算法谱峰搜索函数为：

其中a(θ)为搜索导向矢量，θ∈(-90°,90°)，U_N为噪声子空间。

(3)信源数为先验已知的，则迭代次数选择为k＝1,2,…,K，新算法谱峰搜索函数为：

其中I_M为单位矩阵，W_cg,K为等价于信号子空间的Krylov子空间。

实施例：

下面通过实施例仿真来说明本发明的性能：

仿真条件：假设采用10阵元的阵元间距为d＝λ2的ULA阵型，两个入射信号的方向为θ₁＝10°和θ₂＝30°。为了进一步评价本发明的性能，设置蒙特卡洛实验次数为500，将均方根误差(Root Mean Square Error,RMSE)作为评价指标。

仿真内容和结果：

仿真1，比较本发明、MUSIC算法和CG算法的谱峰搜索图，结果如图2所示。

图2中设置初始参考信号为d₀＝Xa(θ)/||Xa(θ)||，其中θ∈(-90°,90°)。从图中可以看出最小二乘CG算法的空间谱峰比经典CG算法和MUSIC算法的空间谱峰都要尖锐，具有良好的估计性能，并且WCC-CG算法的计算复杂度要低于经典CG算法，从而能够更好的应用到实际中。

仿真2，比较本发明与不同算法的RMSE随输入(Signal-to-NoiseRatio,SNR)的变化情况，其结果如图3所示。

从图3中可以看出RV-CG算法和WCC-CG算法的估计精度都几乎接近于MUSIC算法的估计精度，而RV-CG算法和WCC-CG算法为实值运算和不需要计算接收数据协方差矩阵与特征值分解，所以在精度几乎接近相同的情况下，RV-CG算法和WCC-CG算法的运算量相对而言是较小的，在-20dB到-7dB之间，WCC-CG算法的精度是要优于其它两种算法，在-5dB到0dB之间WCC-CG算法的性能较差一些，但是是在可以接受的范围内。经过对以上三种算法的对比，可以发现RV-CG算法和WCC-CG算法都具备一定的优势，进一步满足实际工程所需要的条件，并且新算法适用于低信噪比和小快拍数的情况。

仿真3，利用不同阵元数来比较本发明与不同算法的计算效率。通过在相同环境中运行MATLAB代码，由具有AMD Ryzen 5 3500U with Radeon Vega Mobile Gfx 2.10GHz和8.00GBRAM的PC端给出仿真结果。从CPU时间的角度对计算效率进行等效评估，结果如图4所示。

图4中是去掉谱峰搜索的算法计算时间，只记录参考信号的计算与梯度迭代得到信号子空间的过程。从图中可以看出CG算法及其改进CG算法的的计算速度是较快于MUSIC算法的，并且阵元数越大，算法的优势就越明显。其中WCC-CG算法的计算速度是最快的，其次为RV-CG算法和CG算法。

本发明实现了通过改进维纳霍夫方程，利用接收数据矩阵代替ACM，避免了设置迭代初始值与迭代过程中需要计算ACM而引起的高复杂度问题，为波达方向的工程化实现提供了技术支持。

Claims (4)

1.一种无需构造协方差矩阵的快速共轭梯度测向算法，其特征在于，接收辐射源信号后，利用最小二乘法改进维纳霍夫方程，获得最小二乘维纳霍夫方程，最小二乘维纳霍夫方程为正态方程，使用正态方程构造仅与接收数据有关的参考信号，设置迭代初始值并进行共轭梯度迭代后，根据信源数估计出信号子空间并获得信号波达方向；所述利用最小二乘法改进维纳霍夫方程具体包括以下步骤：

步骤1-1：维纳霍夫方程为RW＝r_xd，化简维纳霍夫方程，R和r_xd分别表示为：

其中X为接收数据信号矩阵，有X＝[X(0),X(1),…,X(M-1)]，d为参考信号，有d＝[d(0),d(1),…,d(M-1)]；

将以上式子带入到维纳霍夫方程中，有

化简后得到X^HW＝d^H；

步骤1-2：定义最小二乘维纳霍夫方程，为r+X^HW＝d^H，其中，r为最小二乘残差矢量，最小二乘维纳霍夫方程为正态方程；

步骤1-3：为避免计算阵列协方差矩阵R，并且根据最小二乘维纳霍夫方程得到改进的过渡矢量，为v＝X^Hp；

步骤1-4：定义中间变量正态方程残差为q＝Xr，其中q只与阵列接收数据与参考信号有关；

所述构造仅与接收数据有关的参考信号具体包括以下步骤：

步骤2-1：定义参考信号为

其中，a(θ)为搜索导向矢量，θ∈[-90°,90°]，d为参考信号；

步骤2-2：由于第k个阵元的接收信号数据矩阵表示为

则当θ＝θ_i，i∈{1，…，K}时，正态方程残差初始值q_cg,0表示为

可知，正态方程残差初始值q_cg,0为信号导向矢量的线性组合，则q_cg,0张成的K维的Krylov子空间为信号子空间的子集，表示为

span{Q_cg,K}≡span{A(θ)}，

其中

Q_cg,K＝{q_cg,0,q_cg,1,…,q_cg,K-1}；

步骤2-3：当θ≠θ_i，i∈{1，…，K}，正态方程残差初始值q_cg,0表示为：

此时正态方程残差矢量也为导向矢量的线性组合，是由K个信号矢量和搜索导向矢量组合的，即此时的正态方程残差矢量张成的K+1维Krylov子空间也为信号导向矢量张成的子空间的子集，表示为

span{Q_cg,K+1}≡span{A(θ),a(θ)}，

其中

Q_cg,K+1＝{q_cg,0,q_cg,1,…,q_cg,K-1,q_cg,K}。

2.根据权利要求1所述的一种无需构造协方差矩阵的快速共轭梯度测向算法，其特征在于，所述设置迭代初始值并进行共轭梯度迭代具体包括以下步骤：

步骤3-1：共轭梯度初始值设置中，初始权矢量为：w₀＝0，搜索向量初始值：p₁＝-r₁＝-d,其中d为参考信号；

步骤3-2：初始值设置中最小二乘残差矢量初始值为：r₀＝d^H-X^HW，搜索矢量初始值：p₁＝q₀，γ₀＝||q₀||²，正态方程残差矢量初始值为：q₀＝Xr₀；

步骤3-3：当迭代次数为k＝1,2,…,M，得到中间过渡变量：v_k＝X^Hp_k，确定步长：α_k＝γ_i-1/||v_k||²，得到权矢量：w_k＝w_k-1+α_kp_k，利用最小二乘残差：r_k＝r_k-1-α_kv_k，得到正态方程残差：q_k＝Xr_k，确定共轭方向向量：

3.根据权利要求1所述的一种无需构造协方差矩阵的快速共轭梯度测向算法，其特征在于，所述根据信源数估计出信号子空间并获得信号波达方向具体包括以下步骤：

步骤4-1：已知信源数为K，仅需要迭代K次就得到等价于信号子空间的Krylov子空间为：W_cg,K＝{w_cg,0,w_cg,1,…,w_cg,K-1}，MUSIC算法谱峰搜索函数为：

其中a(θ)为搜索导向矢量，θ∈(-90°,90°)，U_N为噪声子空间；

步骤4-2：信源数为先验已知的，则迭代次数选择为k＝1,2,…,K，新算法谱峰搜索函数为：

其中I_M为单位矩阵，W_cg,K为等价于信号子空间的Krylov子空间。

4.根据权利要求1所述的一种无需构造协方差矩阵的快速共轭梯度测向算法，其特征在于，接收辐射源信号表示如下：假设L个相互独立的阵元，以等间距d组成均匀线阵(ULA)，考虑空间中存在K个远场窄带信号入射到阵列， 其中，假设K先验已知，d满足d≤λ2，以避免相位模糊，λ为窄带信号的波长，阵列接收辐射源信号为：

其中，A(θ)为L×K维的阵列流型矩阵，S(t)为K×1维的入射信号矢量，N(t)为L×1维的加性高斯白噪声矢量，α(θ)为A(θ)的列向量，表示为：

Images