本发明涉及互质阵中一种基于联合实值子空间的波达角估计方法，属于阵列信号处理。
背景技术：
波达角(Direction of arrival,DOA)估计是阵列信号处理的一个重要研究方向，其在雷达和移动通信等多个领域有着广泛应用。目前很多经典的超分辨率DOA估计方法已经被提出，如多重信号分类(Multiple signal classification,MUSIC)方法，借助旋转不变技术的参数估计(Estimation of signal parameters via rotational invariance technique,ESPRIT)方法，传播算子方法(Propagator method,PM)以及旋转不变性PM方法等等。然而为了避免估计模糊问题，这些方法全是采用阵元间距不大于半波长的紧凑型阵列，有限的DOF导致了有限的DOA估计分辨率，且存在天线间互耦的干扰。
作为一种非均匀稀疏阵列，嵌套阵在虚拟阵列域能获得远大于实际物理阵元数的自由度(degrees of freedom,DOF)，从而增加系统容量和空间分辨率，然而嵌套阵中依然存在一个紧凑的子阵，较近的天线间距依然会受到天线间互耦的影响，使DOA估计性能降低。另一种当下比较热门的非均匀稀疏阵列是互质阵，其由两个均匀线性子阵列构成，彼此的阵元数和阵元间距存在互质性。互质阵同样能在虚拟域获得较大的DOF，而且相对于嵌套阵，互质阵的阵元分布更加稀疏，因此对互耦问题有较强的鲁棒性。为了解决较大的阵元间距带来的DOA估计模糊问题，文献“DOA estimation with combined MUSIC for coprime array”中的联合MUSIC方法基于互质阵子阵间的互质性，从两个子阵MUSIC谱峰的交集确定了无模糊的DOA估计，但其需要两次全局谱峰搜索，复杂度高。为了降低复杂度，文献“Partial spectral search-based DOA estimation method for co-prime linear arrays”中的局部搜索MUSIC方法利用了子阵的均匀性，将角度域估计问题转变到相位域，从而仅需局部搜索即可恢复出所有的MUSIC谱峰，然后再根据互质性解模糊。然而无论是联合MUSIC还是局部搜索MUSIC，均是对子阵进行分别处理，因此在最终求交集解模糊的时候，多个DOA估计(包括真实解和模糊解)之间会存在干扰，尤其在低信噪比时容易出错。为了避免这个问题，文献“Improved DOA estimation algorithm for co-prime linear arrays using root-MUSIC algorithm”将改进求根MUSIC方法(Li J,Jiang D.Joint elevation and azimuth angles estimation for L-shaped array)应用到了互质阵中，首先通过全体阵列的信号子空间，构造两个子阵的关系矩阵，然后基于全体阵列的噪声子空间利用求根方式求解DOA估计，该方法不但实现了两个子阵输出结果的自动配对，还进一步提升了DOA估计性能。然而，全体阵列的特征分解以及基于全体噪声子空间的估计函数均涉及到较高的复杂度。
技术实现要素：
本发明所要解决的技术问题是提供一种互质阵中基于联合实值子空间的DOA估计方法。首先，我们将互质阵的阵列流形进行了改进，使两个子阵是中心对称的，这样通过酉变换即可使协方差变为实值的，从而大大降低复杂度。此外，我们无需全体阵列的协方差，仅需两个实值的子协方差即可求取两个子阵间的关系矩阵，避免额外配对，进一步降低了复杂度。最后基于两个子阵的噪声子空间构造了联合子空间的DOA估计代价函数。相比于改进求根MUSIC方法，所提算法复杂度大大降低，但DOA估计性能却能几乎保持不变。
本发明为解决上述技术问题采用以下技术方案：
本发明提供一种互质阵中一种基于联合实值子空间的波达角估计方法，所述互质阵的第一和第二子阵的结构满足中心对称，该方法的具体步骤为：
步骤1，构造第一和第二子阵之间的互协方差以及第二子阵的自协方差；
步骤2，通过酉变换对步骤1中的第一和第二子阵之间的互协方差以及第二子阵的自协方差进行实值变换；
步骤3，对步骤2中第一和第二子阵之间的互协方差的实值变换结果进行奇异值分解，得到对应第一子阵的噪声子空间和对应第二子阵的噪声子空间；
步骤4，根据步骤2中第一和第二子阵之间的互协方差以及第二子阵的自协方差的实值变换结果，构造关系矩阵；
步骤5，根据步骤3中对应第一子阵的噪声子空间和对应第二子阵的噪声子空间以及步骤4中的关系矩阵，构造联合实值子空间的角度估计函数；
步骤6，利用求根方式、解的均匀分布性以及关系矩阵，分别得到对应第二子阵的一组DOA估计解以及对应第一子阵的一组DOA估计解；
步骤7，利用第一和第二子阵的互质性，将步骤6中得到的两组DOA估计解中最靠近的两个解的平均作为最终DOA估计。
作为本发明的进一步技术方案，步骤1中，第一和第二子阵之间的互协方差为第二子阵的自协方差为其中，A1＝[a1(α1),...,a1(αK)]、A2＝[a2(β1),...,a2(βK)]分别为第一和第二子阵的方向矩阵，K为空间存在的远场信号的不同方向数，αk＝Nπsinθk，βk＝Mπsinθk，M、N分别为第一子阵、第二子阵的阵元数，M与N为互质的整数，θk为第k个信号的DOA，k＝1,2,…,K，为对角矩阵，表示第k个信号的能量，σ2为噪声能量，IN为N×N的单位矩阵。
作为本发明的进一步技术方案，步骤2中，第一和第二子阵之间的互协方差的实值变换结果为Rr＝H1RsH2H，第二子阵的自协方差的实值变换结果为R2r＝H2RsH2H，其中，分别为A1、A2的酉变换结果，QM和QN分别表示M阶和N阶的酉矩阵。
作为本发明的进一步技术方案，偶数和奇数维度下的酉矩阵分别定义为：和其中，Ip表示p×p的单位矩阵，Πp表示p×p的反向单位矩阵，p为自然数。
作为本发明的进一步技术方案，步骤3具体为：
Rr的奇异值分解为Rr＝UΛVT，其中，U、Λ和V分别为左奇异向量、奇异值和右奇异向量；
对应第一子阵的噪声子空间Un由对应最小M-K个奇异值的左奇异向量构成；
对应第二子阵的噪声子空间Vn由对应最小N-K个奇异值的右奇异向量构成。
作为本发明的进一步技术方案，步骤4中关系矩阵G＝RrR2r+＝H1RsH2H(H2RsH2H)+＝H1H2+，其中，(·)+表示矩阵广义逆。
作为本发明的进一步技术方案，步骤5中联合实值子空间的角度估计函数为
作为本发明的进一步技术方案，步骤6具体为：
利用求根方式对联合实值子空间的角度估计函数进行求解，具体为：
由是一个范德蒙矢量，令z＝ejβ，联合实值子空间的角度估计函数改写为a2(z)＝[1,z,…,zN-1]T，是两个噪声子空间的联合；
利用求根方式对进行求解，在求根结果中选取单位圆内且最靠近单位圆的K个根zk,k＝1,…,K；
由βk＝Mπsinθk，从第二子阵得到一个DOA估计其中，angle(·)表示取相位；
根据所有M个解以2/M为间隔的均匀分布性，得到对应第二子阵的一组DOA估计解其中，u为随着m的变化而变化的整数并保证位于区间[-1,1]之内；
由以及h1(αk)＝Gh2(βk)，得到
由且得到通过a1(αk)相邻元素之间的相位差即可得到对应第一子阵的一个DOA估计为其中a1,m(αk)和a1,m+1(αk)分别表示a1(αk)的第m个和第m+1个元素；
根据所有N个解以2/N为间隔的均匀分布性，得到对应第一子阵的一组DOA估计解其中，v为随着n的变化而变化的整数并保证位于区间[-1,1]之内。
作为本发明的进一步技术方案，第一和第二子阵之间的互协方差以及第二子阵的自协方差通过有限快拍数进行估计：
其中，T表示快拍数，x1(t)、x2(t)分别为第一子阵、第二子阵的输出。
本发明采用以上技术方案与现有技术相比，具有以下技术效果：
1.对互质阵列结构的改进，使其为中心对称的；
2.基于对称的阵列结构，我们将子协方差均转为实值，从而降低了运算复杂度；
3.无需全体阵列的特征分解和噪声空间，我们仅通过实值的子协方差即可得到关系矩阵，避免额外配对，进一步降低了复杂度；
4.最后我们基于两个子阵的实值噪声子空间，构造了DOA估计函数，获得了来自两个子阵的DOA估计信息，并基于互质性进行了解模糊；
5.相比于当下先进的算法，我们提出的算法能大大降低复杂度，但DOA估计性能却几乎保持不变。
附图说明
图1是基本互质阵和改进互质阵阵列结构示意图；
图2是本发明的方法流程示意图；
图3是改进求根MUSIC方法与本发明所提方法的算法复杂度对比示意图；
图4是本发明所提方法在前100次实验下的DOA估计结果；
图5是本发明所提方法与改进求根MUSIC算法的估计性能对比示意图。
具体实施方式
下面结合附图对本发明的技术方案做进一步的详细说明：
本发明提出了互质阵中一种基于联合实值子空间的波达角(Direction of arrival,DOA)估计方法。我们首先对基本互质阵进行了阵列结构的改进，使其为中心对称的，这样阵列协方差即可通过酉变换变成实值的，从而大大降低运算复杂度。之后通过两个子阵构造两个子协方差，并进一步得到子阵间的关系矩阵。基于该关系矩阵，构造联合子空间波达角估计函数。最后从两个子阵的输出的交集中取得无模糊的DOA估计。由于关系矩阵的存在，两个子阵的输出是自动配对的，并且我们联合了两个子空间，所以估计性能得到保证。相比当下先进的算法，我们提出的方法因为仅需要实值特征分解和低维度的子协方差，所以复杂度大大降低，但在DOA估计性能方面却能几乎保持不变。因此本发明在雷达、无线通信中的定位系统中可大大降低复杂度，节省硬件成本，提高响应速度。
如图1所示的本互质阵和改进互质阵阵列结构，图1上半部分给出了基本互质阵的阵列结构，其由两个稀疏的均匀线性阵列(Uniform linear array,ULA)构成，子阵1有M个阵元，阵元间距为Nd，而子阵2有N个阵元，其阵元间距为Md，d为单位间距，一般设为半波长，M与N为互质的整数。该阵列结构不满足中心对称，所以无法进行实值变换。所以我们在图1下半部分给出了改进的互质阵列结构，使两个子阵此时关于Z轴是中心对称的(这里需要说明的是，图1给出的仅为一个示例，在本发明方法的实际应用中只需要保证两个子阵是中心对称即可)。假设此时空间存在来自K个不同方向的远场信号，则两个子阵的输出为
x1(t)＝A1s(t)+n1(t) (1)
x2(t)＝A2s(t)+n2(t) (2)
其中s(t)＝[s1(t),...,sK(t)]T表示K个信源基带信号；n1(t)和n2(t)表示各自阵列上的加性高斯白噪声；A1＝[a1(α1),...,a1(αK)]和A2＝[a2(β1),...,a2(βK)]分别表示子阵1和子阵2的方向矩阵，其列称为导向矢量，表示为
其中αk＝Nπsinθk,βk＝Mπsinθk,θk表示第k个信号的DOA。
实值变换
首先根据(1)-(2),构造两个子阵之间的互协方差为
其中为对角矩阵，包含K个信号能量(Gu J F,Wei P.Joint SVD of Two Cross-Correlation Matrices to Achieve Automatic Pairing in 2-D Angle Estimation Problems)。因为此时的两个子阵是中心对称的，我们可以使用酉变换进行实值转换，从而降低复杂度。奇数和偶数维度下的酉矩阵定义如下：
以子阵1的导向矢量为例，其酉变换之后变为实值导向矢量(假设M为奇数)
对应的实值方向矩阵为同理，子阵2的方向矩阵酉变换之后变为所以(5)式中的互协方差变为
实际中为了保证互协方差的实值性，我们可以对上式取实部，即
基于(8)式中的互协方差，其奇异值分解(Singular value decomposition,SVD)为
Rr＝UΛVT (9)
其中U,Λ和V分别表示左奇异向量，奇异值和右奇异向量。对应最小(M-K)个奇异值的左奇异向量构成子阵1的噪声子空间Un，而对应最小(N-K)个奇异值的右奇异向量构成了对应子阵2的噪声子空间Vn。所以两个子阵的角度信息可从以下代价函数获得
UnHh1(α)＝0 (10)
VnHh2(β)＝0 (11)
但是需要注意的是，此时α和β是分别获得，那么这就会带来两个问题：1.α和β是乱序的，此时需要额外的配对步骤，从而增加复杂度；2.α和β分别是基于各自子阵获得的，相当于只使用了部分阵列孔径，因此估计性能有待提升。所以在进行角度估计之前，我们需要得到两个子阵之间的关系，从而对两个子阵进行联合使用。
关系矩阵构造
构造子阵2的自协方差为
其中σ2为噪声能量，同样对其进行酉变换，并进行噪声消除
其中σ2可从的N-K个较小特征值的平均获得([16]Yunhe C.Joint estimation of angle and Doppler frequency for bistatic MIMO radar)，这里需要注意酉矩阵满足正交性，即所以酉变换不会对噪声能量产生影响。根据(8)式和(13)式，可以构造关系矩阵为
G＝RrR2r+＝H1RsH2H(H2RsH2H)+＝H1H2+ (14)
那么G矩阵给出了两个子阵方向矩阵之间的关系，即
H1＝GH2 (15)
需要注意的是，G矩阵的构造，是基于(8)式和(13)式中的子协方差，相比全体阵列的协方差，其构造复杂度大大降低，此外，这两个子阵还是实值的，复杂度进一步降低。
基于联合子空间的波达角估计
基于G矩阵,将h1(α)＝Gh2(β)代入(10)式,那么基于联合实值子空间的角度估计函数可构造为
(16)式可用求根方式进行求解，因为其中是一个范德蒙矢量。令z＝ejβ,那么其中a2(z)＝[1,z,…,zN-1]T.那么(16)式可以写为
其中是两个噪声子空间的联合。根据(17)式的多项式求根结果，我们选取单位圆内且最靠近单位圆的K个根zk,k＝1,…,K。
因为而βk＝Mπsinθk，所以从子阵2(即β)得到的DOA估计为
从上式可以看到，因为子阵2较大的阵元间距(M>1)，zk的相位值会超出[-π,π]，所以存在相位模糊，我们根据zk从(18)式求出的只是其中一个模糊解，还存在(M-1)个解这M个解应满足所以他们之间的关系([11]SUN F,LAN P and GAO B.Partial spectral search-based DOA estimation method for co-prime linear arrays)为
其中u是个整数，其随着m的变化而变化，并保证位于区间[-1,1]之内。所以所有M个解之间存在均匀分布的关系，且间距为2/M,根据这个关系以及(18)式中的一个解，我们就可以恢复出所有M个解。
现在我们从子阵2中获得了一组解接下来，我们需要考虑来自子阵1的DOA估计。因为以及又因为导向矢量之间满足关系h1(αk)＝Gh2(βk)，所以得到
而且则
其中的因为是标量，可以通过归一化消除。通过a1(αk)相邻元素之间的相位差即可得到对应子阵1的角度估计为
同理，上式为一模糊解，根据所有N个解以2/N为间隔的均匀分布性，以及(22)式获得的其中一个解，同样可以获得一组解
现在我们获得了关于DOA估计的两组解。这两组解包含许多错误解，但也包含有真实的DOA估计。根据ZHOU C,SHI Z,GU Y,et al.DECOM:DOA estimation with combined MUSIC for coprime array，由于互质阵中的互质性，无模糊的DOA可以从两组解的共同解，即交集获得。
但实际中，(5)式和(12)式中的协方差只能通过有限快拍数进行估计，即
其中T表示快拍数。由于残余噪声的影响，最终两组解中往往不存在完全重合的解，所以我们可以从两组解中最接近的两个解的平均作为最终无模糊DOA的估计。因为(18)和(22)获得的解均是基于zk,k＝1,...,K，所以是自动配对的，即这两组解均是对应第k个信源。这样在求解最接近的两个解的时候，不会受到来自其余DOA的影响。
本发明互质阵中一种基于联合实值子空间的波达角(Direction of arrival,DOA)估计方法的步骤，如图2所示：
1.根据(23)-(24)式，构造两个子协方差Rc和R2；
2.基于(8)式和(13)式，对两个子协方差进行实值变换得到Rr和R2r，并根据(9)式对Rr进行奇异值分解，得到对应子阵1的噪声子空间Un和对应子阵2的噪声子空间Vn；
3.由G＝RrR2r+得到关系矩阵G；
4.根据Un，Vn以及关系矩阵G得到(16)式中的联合子空间角度估计函数f(β)；
5.利用求根方式得到对应子阵2的一个DOA估计，并基于(19)式利用解的均匀分布性，得到一组解
6.利用关系矩阵，同理得到对应子阵1的一组DOA估计解最后利用互质性从两组解中最靠近的两个解的平均作为最终DOA估计。
我们所提算法无需额外配对和谱峰搜索，其计算复杂度主要集中在子协方差的构造和关系矩阵的求取以及最后的求根。需要注意的是算法中的子协方差均是实值的，所以计算复杂度只有相应复数矩阵运算的1/4。最终需要总复乘次数约O(MNT+N2T+0.25*(N3+2MN2)+N)次。而改进求根MUSIC方法需要复乘次数约O((M+N)2T+(M+N)3+MNK+2N2K+K3+N)。
在图3中我们给出了典型参数配置下两个算法复杂度的对比，其中M＝4,N＝3,K＝2,T＝100。从图3可以发现，我们提出的算法复杂度只有改进求根MUSIC算法的40％。
实验仿真
实验中，我们采用所提的改进互质阵阵列结构，如图1中所示，其中M＝4,N＝3,以及K＝2个信源，其DOA分别为θ1＝20°和θ2＝35°。各协方差矩阵分别通过T＝100个快拍进行估计。为了衡量DOA估计性能，定义求根均方误差(root mean square error,RMSE)为
其中是第l次蒙特卡洛实验时DOAθk的估计结果，我们共计进行了L＝500次实验。
图4给出了所提算法在前100次实验下的DOA估计结果，可以发现所提算法始终能有效地、准确地估计出两个信源的DOA。
图5则给出了所提算法与改进求根MUSIC算法的估计性能对比，纵轴为RMSE，横轴为SNR，可以发现二者几乎重合，说明所提算法拥有和改进求根MUSIC算法几乎一样的性能，但所提算法的复杂度却大大降低。
本文提出了互质阵中一种基于联合实值子空间的DOA估计算法。我们所做的关键工作可总结如下：
1.对互质阵列结构的改进，使其为中心对称的；
2.基于对称的阵列结构，我们将子协方差均转为实值，从而降低了运算复杂度；
3.无需全体阵列的特征分解和噪声空间，我们仅通过实值的子协方差即可得到关系矩阵；
4.最后我们基于两个子阵的实值噪声子空间，构造了DOA估计函数，获得了来自两个子阵的DOA估计信息，并基于互质性进行了解模糊；
5.相比于当下先进的算法，我们提出的算法能大大降低复杂度，但DOA估计性能却几乎保持不变。
以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。