一种基于L型阵列的信源仰角和方位角估计方法技术领域
本发明涉及阵列信号处理领域，尤其是指一种信源仰角和方位角的估计方法。
背景技术
信源参数估计是阵列信号处理领域的主要研究内容，在无线通信、雷达、声呐诸多领域都有重要的应用价值。信源参数估计主要是对信源的仰角、方位角和距离估计。
目前存在的较高分辨率的算法都是基于子空间方法，子空间方法存在固有框架下的一些缺点。它们都需要较高的信噪比、较多的快拍数或者准确的信源先验信息。所以，这些劣势使得算法的实用性受到了很大的束缚。
近几年，稀疏重构的思想快速地卷入到阵列信号处理中，为上述出现的问题提供了很好的解决方法，稀疏重构算法具有高分辨率。较强的噪声鲁棒性和无需信源准确的先验信息，这为阵列信号处理领域信源参数估计又提供了新的估计方法。但是目前的稀疏重构方法都仅是局限于一维DOA(仰角或方位角)估计，在实际的应用中，当多个信源由同一个仰角(或者方位角)入射的时，仅仅估计信源的一个参数，有时候是很难将其区分开。所以对信源的仰角和方位角同时估计是理论的必然趋势，也具体实际应用的重要价值。
发明内容
本发明提供一种基于L型阵列的信源仰角和方位角估计方法，以解决目前的稀疏重构方法都仅是局限于一维DOA(仰角或方位角)估计的问题。
本发明采取的技术方案是包括下列步骤：
步骤一：有K个远场窄带入射信号，入射到2M个各向同性的传感器阵元组成的L型阵列，该阵列位于y轴与z轴的平面上，为了减少稀疏重构算法的计算复杂度，我们选取不同位置阵元组成的阵列，对其接收的数据进行处理；
步骤二：将位于z轴的阵元组成子阵一Z(t)，该子阵列仅仅包含仰角信息，这样分布式求解信源参数，降低了算法的求解维度，然后对其进行稀疏表示，利用非凸惩罚函数，通过非凸惩罚函数的加权约束方式可以更加有效地抑制噪声，提高参数的估计精度，通过优化工具包求解出仰角值
步骤三：将位于y轴的阵元组成子阵二Y(t)，接收的数据中是既包含仰角信息又包含方位角，将上一步求解出的仰角信息代入到数据模型中，对其得到数据模型稀疏表示，利用奇异值分解，降低算法的计算复杂度，奇异值分解后，利用优化工具包求解凸优化问题，求解出信源的方位角由谱峰位置可以完成仰角和方位角估计值的自动配对。
本发明一种实施方式是在步骤一中包括如下步骤：
(1)该接收数据中第k个信号的仰角为αk，方位角为βk，(1≤k≤K)，d为传感器之间的距离且等长，设位于坐标(0,md,nd)处的阵元为第(m,n)个阵元，(m,n)∈{(M-1,0),…(1,0),(0,0),(0,1),…,(0,M)}，m是位于y轴上第m个传感器，n是位于z轴上第n个传感器，则在某t个采样时刻，第(m,n)个阵元的接收信号为j是虚数单位，式中，sk(t)代表第k个入射信号，nm,n(t)为第(m,n)个传感器的噪声，γk和φk是仰角和方位角的函数，关系如下：γk＝-2πdsinαkcosβk/λ，φk＝-2πdcosαk/λ，其中，λ代表信号波长；
(2)下面将2M个传感器的接收数据写成向量形式：F(t)＝AS(t)+n(t)，S(t)＝[s1(t),…,sK(t)]T代表K个信号，n(t)是nm,n(t)的向量形式，方向矩阵A＝A(αk,βk)＝[a(α1,β1),…,a(αK,βK)]，其中，a(αk,βk)是方向矩阵A中的列向量，[·]T是矩阵的转置运算。
本发明一种实施方式是在步骤二中包括如下步骤：
(1)位于z轴的阵元组成子阵一Z(t)＝[f0,1(t),…,f0,M(t)]T＝A1S(t)+n1(t)，子阵一的方向矩阵A1表示为A1＝[a1(α1),…,a1(αK)]，其中k＝1,…,K；
(2)对子阵列一进行稀疏表示，其中，n1(t)＝[n0,1(t),…,n0,M(t)]T代表子阵一的接收噪声，Q1表示在仰角估计范围内划分的网格数，G(t)＝[g(t1),…,g(tT)]，g(t1)～g(tT)是具有相同稀疏结构的稀疏向量；
(3)利用非凸惩罚函数和DC分解将非凸函数分解为凸惩罚函数，非凸惩罚函数的表达式为，p(δ)＝ξlog(δ+ε)-ξlog(δ)，ξ是非凸惩函数表达式的正则化参数，可以去经验值，ε＝0.001，该非凸惩罚函数可以更加有效地抑制噪声，提高参数的估计精度，即得到：在这里，W表示由非凸惩罚函数迭代转化的加权系数，W＝ξ/(ξ+|gj|)，j＝1,…,Q1，通过凸优化工具包CVX迭代，求解出仰角值
本发明一种实施方式是在步骤三中包含如下步骤：
(1)位于y轴的阵元组成的阵列作为子阵二，其接受的数据为：Y(t)＝[f0,0(t),f1,0(t)…,fM-1,0(t),f0,1(t)]T＝A2S(t)+n2(t)，子阵二的方向矩阵A2表示为A2＝A2(αk,βk)＝[a2(α1,β1),…,a2(αK,βK)]，其中，已知接收的数据中是既包含仰角信息又包含方位角；
(2)将上一步求解出的仰角估计值代入到数据模型中，对其得到矩阵的输出稀疏表示，利用已经得到的仰角估计值辅助构造： A ~ 2 = [ A 2 ( α ^ 1 , θ ‾ ) , A 2 ( α ^ 2 , θ ‾ ) , ... , A 2 ( α ^ K , θ ‾ ) ] , ]]>其中， A 2 ( α ^ k , θ ‾ ) = [ a 2 ( α ^ k , θ 1 ) , a 2 ( α ^ k , θ 2 ) , ... , a 2 ( α ^ k , θ Q 2 ) ] , ]]>k＝1,…,K，Q2是在方位角估计范围划分的网格数。
(3)子阵二的稀疏表示为：H(t)是和G(t)一样具有相同稀疏结构的稀疏信号。对Y进行奇异值分解以分离信号子空间和噪声子空间：利用l1-SVD算法求解：这里，η代表正则化参数，最后通过凸优化工具包CVX求解上述式子得到方位角的估计值
本发明利用L型阵列，通过选取不同的阵元，组成两个新的子阵列，由两个子阵列的接收数据，可以看出其中一个子阵列仅仅包含仰角信息，另外一个既包含仰角信息又包含方位角信息。
本发明的优点：先利用子阵一列求解出仰角，利用仰角估计值去求解方位角。通过选取不同的子阵列，我们降低了稀疏重构在面对基于稀疏重构的多参数联合估计时的高维网格划分难的问题。利用两个子阵列的先后求解的过程，完成了参数的自动配对。特别要提及的是，在利用稀疏重构的数学工具求解的过程中，我们采用了非凸惩罚函数的形式，能更好地抑制噪声对信号的干扰，提高了估计的精度。
附图说明
图1是本发明L型阵列的示意图。
具体实施方式
步骤一：有K个远场窄带入射信号sk(t)，入射到由2M个各向同性的传感器组成的L型阵列，该阵列位于y轴与z轴的平面上，为了减少稀疏重构算法对仰角和方位角估计的计算复杂度，我们选取不同阵元组成的阵列，对其接收的数据进行处理，该接收数据中第k个信号的仰角为αk，方位角为βk，(1≤k≤K)，d为传感器之间的距离且等长，设位于坐标(0,md,nd)处的阵元为第(m,n)个阵元，(m,n)∈{(M-1,0),…(1,0),(0,0),(0,1),…,(0,M)}，m是位于y轴上第m个传感器，n是位于z轴上第n个传感器，则在某t个采样时刻，第(m,n)个阵元的接收信号为j是虚数单位，式中，sk(t)代表第k个入射信号，nm,n(t)为第(m,n)个传感器的噪声，γk和φk是仰角和方位角的函数，关系如下：γk＝-2πdsinαkcosβk/λ，φk＝-2πdcosαk/λ，其中，λ代表信号波长，下面将2M个传感器的接收数据写成向量形式：F(t)＝AS(t)+n(t)，S(t)＝[s1(t),…,sK(t)]T代表K个信号，n(t)是nm,n(t)的向量形式，方向矩阵A＝A(αk,βk)＝[a(α1,β1),…,a(αK,βK)]，其中，a(αk,βk)是方向矩阵A中的列向量，[·]T是矩阵的转置运算；
步骤二：将位于z轴的所有阵元组成阵列作为子阵一，其接收数据为Z(t)＝[f0,1(t),…,f0,M(t)]T＝A1S(t)+n1(t)，其中，n1(t)＝[n0,1(t),…,n0,M(t)]T代表子阵一的接收噪声，子阵一的方向矩阵A1表示为A1＝A1(αk)＝[a1(α1),…,a1(αK)]，其中，k＝1,…,K，该子阵列仅仅包含仰角信息，这样分布式求解信源参数，降低了算法的求解维度，然后对子阵列一进行稀疏表示，Q1表示在仰角估计范围内划分的网格数，G(t)＝[g(t1),…,g(tT)]，g(t1)～g(tT)是具有相同稀疏结构的稀疏向量，当入射信号的角度等于θi时，向量g(t1)、…、g(tT)的第i个分量值不等于零，而其它分量的值均为零；
本发明利用非凸惩罚函数和DC分解将非凸函数分解为凸惩罚函数，这样能保证全局最优解，该非凸惩罚函数的表达式为，p(δ)＝ξlog(δ+ε)-ξlog(δ)，ξ是非凸惩函数表达式的正则化参数，可以取经验值，ε＝0.001，该非凸惩罚函数可以更加有效地抑制噪声，提高参数的估计精度；
非凸惩罚函数被DC分解为两部分：
p(δ)＝p1(δ)-p2(δ)＝ξ|δ|-(ξ|δ|-ξlog(δ+ε)+ξlog(δ))
定义并带入到稀疏子阵一中，经过数学的迭代推导之后。即得到：
m i n | | Z - A 1 ( θ ~ ) G | | f 2 + W | | g ( l 2 ) | | 1 ]]>
在这里，W表示由非凸惩罚函数迭代转化的加权系数，W＝ξ/(ξ+|gj|)，j＝1,…,Q1，通过凸优化工具包CVX迭代，求解出仰角值
步骤三：将位于y轴的阵元组成的阵列作为子阵二，其接受的数据为：Y(t)＝[f0,0(t),f1,0(t)…,fM-1,0(t),f0,1(t)]T＝A2S(t)+n2(t)，n2(t)＝[n0,0(t),n1,0…,nM-1,0(t),n0,1(t)]T代表子阵二各个阵元上所叠加的噪声，子阵二的方向矩阵A2表示为A2＝A2(αk,βk)＝[a2(α1,β1),…,a2(αK,βK)]，其中，已知接收的数据中是既包含仰角信息又包含方位角，将上一步求解出的仰角估计值代入到数据模型中，对其得到矩阵的输出稀疏表示；
利用已经得到的仰角估计值辅助构造：
A ~ 2 = [ A 2 ( α ^ 1 , θ ‾ ) , A 2 ( α ^ 2 , θ ‾ ) , ... , A 2 ( α ^ K , θ ‾ ) ] , ]]>
其中，k＝1,…,K，Q2是在方位角估计范围划分的网格数；通过这样构造的稀疏表达式可以看出，相当于分别在已知的处，对方位角β划分网格，(从几何角度可以得到，是分别在每个得到的仰角估计值方向上，得到对方位角的一维搜索)。
则子阵二的稀疏表示为：
Y ( t ) = A ~ 2 H ( t ) + N 2 ( t ) ]]>
H(t)是和G(t)一样具有相同稀疏结构的稀疏信号；
对Y进行奇异值分解以分离信号子空间和噪声子空间：
Y＝UEVH
U和V分别是匹配Y的行和列的酉矩阵，E是对角矩阵。定义这里WK,0＝[IK,0]，其中，IK为K阶单位矩阵，0是K×(T-K)的零矩阵；有如下的表达式：
Y 11 = A ~ 2 H 11 + N 11 ]]>
定义表示矩阵H11的第(q,k)个元素，则代表矩阵H11的第q行的转置，即表示第q行向量的l2范数，定义该向量具有稀疏性，并保留原有的稀疏结构，其非零元素的位置代表了信号的入射角度；利用l1-SVD算法求解：
m i n | | Y 11 - A ~ 2 H 11 | | f 2 + η | | h ( l 2 ) | | 1 ]]>
这里，η代表正则化参数，最后通过凸优化工具包CVX求解上述式子得到方位角的估计值