CN110174657B - 基于秩一降维模型和块矩阵恢复的波达方向估计方法 - Google Patents
基于秩一降维模型和块矩阵恢复的波达方向估计方法 Download PDFInfo
- Publication number
- CN110174657B CN110174657B CN201910516081.3A CN201910516081A CN110174657B CN 110174657 B CN110174657 B CN 110174657B CN 201910516081 A CN201910516081 A CN 201910516081A CN 110174657 B CN110174657 B CN 110174657B
- Authority
- CN
- China
- Prior art keywords
- matrix
- vector
- representing
- signal
- shift
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 305
- 238000011084 recovery Methods 0.000 title claims abstract description 73
- 238000000034 method Methods 0.000 title claims abstract description 43
- 239000013598 vector Substances 0.000 claims abstract description 212
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 29
- 238000005259 measurement Methods 0.000 claims abstract description 22
- 238000005457 optimization Methods 0.000 claims abstract description 11
- 230000008707 rearrangement Effects 0.000 claims description 5
- 241000543540 Guillardia theta Species 0.000 claims description 3
- 239000000654 additive Substances 0.000 claims description 3
- 230000000996 additive effect Effects 0.000 claims description 3
- 150000001875 compounds Chemical class 0.000 claims description 3
- 238000005192 partition Methods 0.000 claims 1
- 230000017105 transposition Effects 0.000 claims 1
- 238000004088 simulation Methods 0.000 description 9
- 238000000354 decomposition reaction Methods 0.000 description 3
- 238000000342 Monte Carlo simulation Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 101150040772 CALY gene Proteins 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
Landscapes
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于秩一降维模型和块矩阵恢复的波达方向估计方法,该方法的具体步骤为:建立雷达的接收信号模型,确定接收信号的测量矩阵;根据接收信号的测量矩阵,构造基于秩一降维的信号协方差矩阵模型;对基于秩一去噪的信号协方差矩阵模型进行稀疏重构,得到稀疏信号矢量;根据稀疏重构信号矢量,采用交替网格优化算法估计目标信源的波达方向。本发明能够重建无噪声信号协方差矩阵,避免了信号协方差矩阵和自由度的信息丢失,提高了对密集信源和多目标的DOA估计性能。
Description
技术领域
本发明属于雷达信号处理技术领域,尤其涉及基于秩一降维模型和块矩阵恢复的波达方向估计方法。
背景技术
众所周知,在非均匀噪声和有限快拍数的条件下,经典的波达方向(DOA)估计方法性能较差。
研究表明,用于波达方向估计的传统信号子空间方法,在具有足够快拍的高斯白噪声的条件下,对于波达方向估计的性能表现良好。利用稀疏重构算法的波达方向估计方法可以通过找到目标函数的稀疏解来实现信号的重建。顾名思义,稀疏重构算法计算量很大,对于快拍数的要求同样也很高。基于压缩感知理论的L1奇异值分解(L1-SVD)算法,该算法利用1范数约束和SVD可以在低信噪比环境下提高DOA估计性能。然而,上述稀疏重建算法计算复杂度大,需要的快拍数也相当多。基于具有反馈的零空间调整算法(NST+HT+FB)的交替网格优化(AGO)算法的提出很好的解决了快拍数和耗时问题。但是由于该方法会导致阵列孔径减小,因此造成DOA估计性能某种程度的损失。
基于子空间和稀疏重构的算法的优越性仅体现于均匀白噪声条件下,该方法在非均匀高斯白噪声条件下不能体现出其优势。可以通过相对于信号和噪声参数处理对数似然函数来迭代地确定输入信号的DOA估计。该方法包含的迭代过程和高度非线性优化问题导致耗时量大。
以迭代方式确定信号/噪声子空间的DOA估计方法的原理是通过求解对数似然函数或最小二乘问题来估计噪声协方差矩阵和信号/噪声子空间,然后可以从估计的信号/噪声子空间或噪声协方差矩阵完成DOA估计。但是该方法由于包含迭代过程导致耗时较大,而通过对阵列协方差矩阵进行矢量化并充分利用虚拟阵列结构能够改善耗时问题。但是该方法通过消除包括噪声方差的元素从而获得用于DOA估计的无噪声协方差稀疏向量,其中的移除操作可能导致信号协方差矩阵中包含的部分信息丢失。因此,现有传统的稀疏重构算法对于密集信源和多目标的DOA估计性能不是很理想。
发明内容
针对上述问题,本发明的目的是提出一种基于秩一降维模型和块矩阵恢复的波达方向估计方法。本发明基于秩一降维模型,得到降维的信号协方差矩阵,以消除非均匀噪声,再通过直接计算块矩阵,从降维的信号协方差矩阵中估计噪声协方差矩阵,从而重建无噪声信号协方差矩阵,避免了信号协方差矩阵和自由度的信息丢失;再采用稀疏重构和交替网格优化(AGO)算法,显著提高了对密集信源和多目标的DOA估计性能。
为了达到上述目的,本发明采用以下技术方案予以解决。
基于秩一降维模型和块矩阵恢复的波达方向估计方法,包括以下步骤:
步骤1,建立雷达的接收信号模型x(t),确定接收信号的测量矩阵X。
(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列,则t时刻接收到的目标信源的回波信号x(t),其表达式为:
x(t)=As(t)+n(t);
其中,A为导向矩阵,s(t)是信号波形矢量,n(t)=[n1(t),n2(t),…,nM(t)]T为零均值加性非均匀复数高斯白噪声向量,且n(t)~CN(0,Q),Q是与n(t)相关的噪声协方差矩阵。
(1.2)导向矩阵A的表达式为:
A=[a(θ1),a(θ2),a(θl),…,a(θL)];
其中,L是目标信源数,l=1,2,...,L;a(θl)是M×1的导向矢量,M代表雷达阵元数,θl代表第l个入射信号角度,θl∈Θ,Θ表示角度搜索范围,其表达式为:
a(θl)=[1,e-jα,…,e-j(M-1)α]T;
其中,α表示对应第l个DOA的阵列几何形状,α=2πdsin(θl)/λ,d代表阵元间距,λ代表波长;
s(t)是信号波形矢量,其表达式为:
s(t)=[s1(t),s2(t),…,sL(t)]T;
(1.3)根据接收信号模型x(t),得到接收信号的测量矩阵X:
其中,J是快拍数,tz表示第z个快拍时刻,x(tz)表示tz快拍时刻天线阵列所接收的M×1维回波信号数据,x1(tz)表示tz快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据,*表示共轭操作。
步骤2,根据接收信号的测量矩阵X,构造基于秩一降维模型的信号协方差矩阵,得到降维的信号协方差矩阵。
其中,d代表阵元间距,λ代表波长;
(3.3)将导向矩阵A表示为块矩阵形式:
其中,块矩阵A11、A21的维数分别为L×L,块矩阵A31的维数分别为(M-2L)×L。
其中,上标H表示矩阵的共轭转置。
步骤4,对无噪声信号协方差矩阵R0进行稀疏重构,得到稀疏重构信号矢量r。
(4.2)将信号矢量r'进行线性表示为:
r'=B(θ)P;
(4.3)根据实际回波信号在空域中的稀疏性,采用重构算法对信号矢量r'进行稀疏重构,得到稀疏重构信号矢量r:
其中,表示划分的一组字典角度,一般 是具有2(M-1)<<Nθ属性的过完备字典,且满足受限等距属性(RIP);表示字典角度为时的稀疏导向矢量,即第1列稀疏导向矢量,Nθ是字典数;表示稀疏信号功率矢量,其仅有L列非零系数,Nθ>>L;表示第1个稀疏信号功率。
步骤5,根据稀疏重构信号矢量r,采用交替网格优化算法,估计目标信源的波达方向。
与现有技术相比,本发明的有益效果为:
(1)本发明采用基于秩一降维模型获得降维的信号协方差矩阵,通过对该矩阵进行块矩阵恢复方法得到无噪声信号协方差矩阵,可以避免信号协方差矩阵和自由度的信息丢失。
(2)本发明通过对无噪声信号协方差矩阵进行矢量化稀疏重构,可以显著提高DOA估计性能和抗均匀噪声的鲁棒性。
附图说明
下面结合附图和具体实施例对本发明做进一步详细说明。
图1是本发明的实现流程图;
图2是本发明实施例中分别采用交替网格优化算法(AGO)、基于压缩感知理论的L1奇异值分解算法(L1-SVD)和He算法,在信噪比为0时的目标信源波达方向的估计结果图;
图3是本发明实施例中采用本发明方法在信噪比为0时的目标信源波达方向的估计结果图;
图4是本发明实施例中对网格优化算法(AGO)、基于压缩感知理论的L1奇异值分解算法(L1-SVD)、He算法和本发明方法,在不同信噪比条件下,目标信源波达方向的估计均方根误差与信噪比的关系图。
具体实施方式
下面结合附图对本发明的实施例及效果作进一步详细描述。
参照图1,本发明的实现步骤如下:
步骤1,建立雷达的接收信号模型x(t),确定接收信号的测量矩阵X。
(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列,则t时刻接收到的目标信源的回波信号x(t),其表达式为:
x(t)=As(t)+n(t)
其中,A为导向矩阵,s(t)是信号波形矢量,n(t)=[n1(t),n2(t),…,nM(t)]T为零均值加性非均匀复数高斯白噪声向量,且n(t)~CN(0,Q),Q是与n(t)相关的噪声协方差矩阵。
(1.2)导向矩阵A的表达式为:
A=[a(θ1),a(θ2),a(θl),…,a(θL)];
其中,L是目标信源数,l=1,2,...,L;a(θl)是M×1的导向矢量,M代表雷达阵元数,θl代表第l个入射信号角度,θl∈Θ,Θ表示角度搜索范围,其表达式为:
a(θl)=[1,e-jα,…,e-j(M-1)α]T;
其中,α=2πdsin(θl)/λ对应第l个DOA的阵列几何形状,d代表阵元间距,λ代表波长;
信号波形矢量s(t)的表达式为:
s(t)=[s1(t),s2(t),…,sL(t)]T;
(1.3)根据接收信号模型x(t),得到接收信号的测量矩阵X:
其中,J是快拍数,tz表示第z个快拍时刻,x(tz)表示tz快拍时刻天线阵列所接收的M×1维回波信号数据,x1(tz)表示tz快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据,*表示共轭操作。
(3.3)将导向矩阵A表示为块矩阵形式:
其中,块矩阵A11、A21的维数分别为L×L,块矩阵A31的维数分别为(M-2L)×L。
其中,上标H表示矩阵的共轭转置。
步骤4,对无噪声信号协方差矩阵R0进行稀疏重构,得到稀疏重构信号矢量r。
(4.2)将信号矢量r'进行线性表示为:
r'=B(θ)P;
其中,P表示信号功率矢量,P=[P1,P2,…,PL]T;B(θ)表示虚拟流形矩阵,B(θ)=[b(θ1),…,b(θl),…,b(θL)],第l列导向矢量b(θl)为
(4.3)根据实际回波信号在空域中的稀疏性,采用重构算法对信号矢量r'进行稀疏重构,得到稀疏重构信号矢量r:
步骤5,根据稀疏重构信号矢量r,采用交替网格优化算法估计目标信源的波达方向。
其具体实施步骤为:
(5.1)按照下式,计算导向矩阵A的右逆矩阵:
Φa=AH(A·AH)-1;
其中,Φa表示导向矩阵A的右逆矩阵,A表示导向矩阵,H表示共轭转置操作,上标-1表示求逆操作;
(5.2)按照下式,计算导向矩阵A的正交投影矩阵:
Φ=I-ΦaA;
其中,Φ表示导向矩阵A的正交投影矩阵,I表示单位矩阵,Φa表示导向矩阵A的右逆矩阵,A表示导向矩阵;
(5.3)按照下式,计算初始恢复矢量:
(5.4)估计信源数:
(5.4a)将动态信源数初始化为1;
(5.4b)按照下式,对第k次内循环中的恢复矢量进行排序操作:
其中,表示恢复矢量取模值后按降序重排的矢量,l表示外循环次数,k表示内循环次数,T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集,|·|表示取模值操作,sort(|·|,'descend')表示降序排列操作;
(5.4c)按照下式,计算第k+1次内循环中的恢复矢量:
(5.4d)按照下式,计算内循环相对误差值:
(5.4e)判断内循环相对误差值H2是否大于10-3,若是,则执行步骤(5.4b),否则,执行步骤(5.4f);
(5.4f)将第l次外循环中的动态信源数加1,用加1后的动态信源数作为下一次外循环中的动态信源数;
(5.4g)按照下式,计算失配相对误差:
其中,γl+1表示第l+1次外循环中的失配相对误差,l表示外循环次数,A表示导向矩阵,uk表示第k次内循环中的中间辅助矢量,X表示测量矩阵,||·||2表示取2范数操作;
(5.4h)按照下式,计算外循环相对误差值:
H1=γl+1-γl|;
其中,H1表示外循环相对误差值,γl+1和γl分别表示外循环次数为l+1和l时的失配相对误差,|·|表示取模值操作;
(5.4i)判断外循环相对误差值H1是否大于0.05,若是,则执行步骤(4.4b),否则,执行步骤(4.4j);
(5.4j)将外循环结束时的动态信源数值作为信源数的估计值;
(5.5)第一次估计目标到达角:
(5.5a)按照下式,寻找峰值位置矢量:
其中,pV表示恢复矢量取模值后按降序重排的峰值矢量,pI表示恢复矢量元素取模值降序重排后,恢复矢量原下标值经过重新排列得到的峰值位置矢量,findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.5b)将峰值位置矢量中的第一个元素值,放入信源位置矢量的第一个位置,将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置,直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值,得到最终信源位置矢量;
(5.6)按照下式,寻找峰值位置矢量:
其中,pV表示恢复矢量取模值后按降序重排的峰值矢量,pI表示恢复矢量元素取模值降序重排后,恢复矢量原下标值经过重新排列得到的峰值位置矢量,findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.7)按照下式,计算当前代价函数:
其中,F(M)表示当前代价函数,pV(s)表示峰值矢量pV中的第s个元素,s表示动态信源数,∑(·)表示求和,表示恢复矢量中的第i个元素, 表示不属于符号,表示峰值位置矢量pI的前个元素,表示信源数的估计值,|·|表示取模值操作;
(5.8)按照下式,计算当前峰值:
(5.9)计算左移代价函数:
(5.9a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围;
(5.9b)按照下式,计算左移导向矩阵:
A(L)=[a(θ1),a(θ2),…,a(θi),…]
其中,A(L)表示左移导向矩阵,M表示天线阵列中的阵元个数,θi∈Θ(L),∈表示属于符号,Θ(L)表示左移角度搜索范围,i表示搜索角度θi在左移角度搜索范围Θ(L)中的序号;
(5.9c)按照下式,计算左移导向矩阵的正交投影矩阵:
Φ(L)=I-(A(L))H(A(L)(A(L))H)-1A(L)
其中,Φ(L)表示左移导向矩阵A(L)的正交投影矩阵,I表示单位矩阵,A(L)表示左移导向矩阵,H表示共轭转置操作,-1表示求逆操作;
(5.9d)按照下式,计算左移恢复矢量:
(5.9e)按照下式,寻找左移峰值位置矢量:
其中,表示左移恢复矢量取模值后按降序重排的左移峰值矢量,表示左移恢复矢量元素取模值降序重排后,左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量,findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.9f)按照下式,计算左移代价函数:
其中,F(L)表示左移代价函数,表示左移峰值矢量中的第s个元素,s表示动态信源数,∑(·)表示求和,表示左移恢复矢量中的第i个元素, 表示不属于符号,表示左移峰值位置矢量的前个元素,表示信源数的估计值,|·|表示取模值操作;
(5.10)计算右移代价函数:
(5.10a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围;
(5.10b)按照下式,计算右移导向矩阵:
A(R)=[a(θ1),a(θ2),…,a(θi),…]
其中,A(R)为右移导向矩阵,M表示天线阵列中的阵元个数,θi∈Θ(R),∈表示属于符号,Θ(R)表示右移角度搜索范围,i表示搜索角度θi在右移角度搜索范围Θ(R)中的序号;
(5.10c)按照下式,计算右移导向矩阵的正交投影矩阵:
Φ(R)=I-(A(R))H(A(R)(A(R))H)-1A(R)
其中,Φ(R)表示右移导向矩阵A(R)的正交投影矩阵,I表示单位矩阵,A(R)表示右移导向矩阵,H表示共轭转置操作,-1表示求逆操作;
(5.10d)按照下式,计算右移恢复矢量:
(5.10e)按照下式,寻找右移峰值位置矢量:
其中,表示右移恢复矢量取模值后按降序重排的右移峰值矢量,表示右移峰值矢量元素取模值降序重排后,右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量,findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.10f)按照下式,计算右移代价函数:
其中,F(R)表示右移代价函数,表示右移峰值矢量中的第s个元素,s表示动态信源数,∑(·)表示求和,表示右移恢复矢量中的第i个元素, 表示不属于符号,表示右移峰值位置矢量的前个元素,表示信源数的估计值,|·|表示取模值操作;
(5.11)更新当前栅格参数:
(5.11a)判断左移代价函数是否同时大于当前代价函数和右移代价函数,若是,则执行步骤(5.11b),否则,执行步骤(5.11e);
(5.11b)将步骤(5.7)的当前代价函数的值更新为左移代价函数的值;
(5.11c)将恢复矢量的元素值更新为左移恢复矢量的元素值;
(5.11d)按照下式,更新目标到达角矢量:
(5.11e)判断右移代价函数是否同时大于当前代价函数和左移代价函数,若是,则执行步骤(5.11f),否则,执行步骤(5.12);
(5.11f)将当前代价函数的值更新为右移代价函数的值;
(5.11g)将恢复矢量的元素值更新为右移恢复矢量的元素值;
(5.11h)按照下式,更新目标到达角估计矢量:
(5.12)判断可调栅格步长△是否大于栅格优化精度ξθ,若是,则执行步骤(5.9),否则,执行步骤(5.13);
(5.13)按照下式,更新目标到达角估计矢量:
(5.14)判断动态信源数是否小于信源数,若是,则执行步骤(5.6),否则,执行步骤(5.15);
(5.15)获得第二次目标到达角估计矢量:
将步骤(5.13)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。
仿真实验
本发明的效果可通过以下仿真实验进一步说明。
(1)仿真参数:
本发明采用M=24个阵元,且阵元间距为半波长(d=λ2)的均匀线性阵列。采用在-90°至90°范围内具有1°间隔的一般网格。阈值τ和误差常数ξ分别为5M和5×10-4。假定非均匀噪声协方差矩阵为Q=diag{Qn,3*Qn,1.5*Qn},其中Qn=[2.0,10,2.5,5.0,0.5,1.5,3.0,5.0]。采用三个不相关的密集间隔信号源,其中DOA为{-5°,1°,3°}。
参数设置如表1:
表1***仿真参数
(2)仿真内容:
仿真1,在上述仿真参数下,分别采用AGO算法、L1-SVD算法和He算法在信噪比为0的条件下,对目标的回波信号进行波达方向估计,结果如图2所示。
由图2可以看出,当信号源密集放置时,AGO和He方法无法区分信号;此外,经典的L1-SVD方法性能显著恶化。
仿真2,在上述仿真参数下,采用本发明方法在信噪比为0的条件下,对回波信号进行波达方向估计,结果如图3所示。
由图3可以看出,本发明方法能够确定三个密集的不相关信号;比较图2和图3可以看出,本发明方法在旁瓣较低的情况下,由于获得了无噪声信号协方差矩阵,估计出了三个密集不相关信号的波达方向,说明本发明方法的波达方向估计性能优于AGO、L1-SVD和He方法。
在上述仿真参数下,分别采用AGO算法、L1-SVD算法、He算法和本发明方法在信噪比为分别为-6dB、-4dB、-2dB、0dB、2dB、4dB、6dB的条件下,对回波信号进行波达方向估计,并采用均方根误差RMSE对估计结果进行评价:
以上试验结果如图4所示,由图4可以看出,随着信噪比SNR的增大,四种方法的估计均方根误差RMSE均呈现下降趋势,且本发明方法的估计均方根误差在各信噪比时,较其他算法的估计均方根误差最小,说明本发明方法在不同信噪比的情况下,均具有较优的波达方向估计性能。
本领域普通技术人员可以理解:实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成,前述的程序可以存储于一计算机可读取存储介质中,该程序在执行时,执行包括上述方法实施例的步骤;而前述的存储介质包括:ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应以所述权利要求的保护范围为准。
Claims (8)
2.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法,其特征在于,步骤1的具体步骤为:
(1.1)设置信号接收雷达为包含M个阵元的均匀直线阵列,则t时刻接收到的目标信源的回波信号,即t时刻的接收信号x(t)的表达式为:
x(t)=As(t)+n(t);
其中,A为导向矩阵,s(t)是信号波形矢量,n(t)=[n1(t),n2(t),…,nM(t)]T为零均值加性非均匀复数高斯白噪声向量,且n(t)~CN(0,Q),Q是与n(t)相关的噪声协方差矩阵;
(1.2)根据接收信号模型x(t),得到接收信号的测量矩阵X:
其中,J是快拍数,tz表示第z个快拍时刻,x(tz)表示tz快拍时刻天线阵列所接收的M×1维回波信号数据,x1(tz)表示tz快拍时刻天线阵列中第一个天线阵元所接收到的回波信号数据,*表示共轭操作。
3.根据权利要求2所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法,其特征在于,所述导向矩阵A的表达式为:
A=[a(θ1),a(θ2),a(θl),…,a(θL)];
其中,L是目标信源数,l=1,2,...,L;M代表雷达阵元数,θl代表第l个入射信号角度,θl∈Θ,Θ表示角度搜索范围;a(θl)是M×1的导向矢量,其表达式为:
a(θl)=[1,e-jα,…,e-j(M-1)α]T;
其中,α表示对应第l个DOA的阵列几何形状,α=2πdsin(θl)/λ,d代表阵元间距,λ代表波长;
所述信号波形矢量s(t)的表达式为:
s(t)=[s1(t),s2(t),…,sL(t)]T;
其中,sL(t)表示第L个信号波形。
4.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法,其特征在于,步骤2按照以下步骤实施:
其中,d代表阵元间距,λ代表波长,θl代表第l个入射信号角度;
6.根据权利要求5所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法,其特征在于,所述通过块矩阵恢复出无噪声信号协方差矩阵R0,其具体步骤为:
(3.2)将导向矩阵A表示为块矩阵形式:
其中,块矩阵A11、A21的维数分别为L×L,块矩阵A31的维数为(M-2L)×L;
其中,上标H表示矩阵的共轭转置;
7.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法,其特征在于,步骤4按以下步骤实施:
(4.2)将信号矢量r′进行线性表示为:
r′=B(θ)P;
其中,P表示信号功率矢量,P=[P1,P2,…,PL]T;B(θ)表示虚拟流形矩阵,B(θ)=[b(θ1),…,b(θl),…,b(θL)],第l列导向矢量b(θl)为θl代表第l个入射信号角度,d代表阵元间距,λ代表波长;
(4.3)根据实际回波信号在空域中的稀疏性,采用重构算法对信号矢量r′进行稀疏重构,得到稀疏重构信号矢量r:
8.根据权利要求1所述的基于秩一降维模型和块矩阵恢复的波达方向估计方法,其特征在于,步骤5的按照以下步骤实施:
(5.1)按照下式,计算导向矩阵A的右逆矩阵:
Φa=AH(A·AH)-1,
其中,Φa表示导向矩阵A的右逆矩阵,A表示导向矩阵,H表示共轭转置操作,上标-1表示求逆操作;
(5.2)按照下式,计算导向矩阵A的正交投影矩阵:
Φ=I-ΦaA;
其中,Φ表示导向矩阵A的正交投影矩阵,I表示单位矩阵,Φa表示导向矩阵A的右逆矩阵,A表示导向矩阵;
(5.3)按照下式,计算初始恢复矢量:
(5.4)估计信源数:
(5.4a)将动态信源数初始化为1;
(5.4b)按照下式,对第k次内循环中的恢复矢量进行排序操作:
其中,表示恢复矢量取模值后按降序重排的矢量,l表示外循环次数,k表示内循环次数,T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集,|·|表示取模值操作,sort(|·|,′descend′)表示降序排列操作;
(5.4c)按照下式,计算第k+1次内循环中的恢复矢量:
(5.4d)按照下式,计算内循环相对误差值:
(5.4e)判断内循环相对误差值H2是否大于10-3,若是,则执行步骤(5.4b),否则,执行步骤(5.4f);
(5.4f)将第l次外循环中的动态信源数加1,用加1后的动态信源数作为下一次外循环中的动态信源数;
(5.4g)按照下式,计算失配相对误差:
其中,γl+1表示第l+1次外循环中的失配相对误差,l表示外循环次数,A表示导向矩阵,uk表示第k次内循环中的中间辅助矢量,X表示测量矩阵,||·||2表示取2范数操作;
(5.4h)按照下式,计算外循环相对误差值:
H1=|γl+1-γl|;
其中,H1表示外循环相对误差值,γl+1和γl分别表示外循环次数为l+1和l时的失配相对误差,|·|表示取模值操作;
(5.4i)判断外循环相对误差值H1是否大于0.05,若是,则执行步骤(4.4b),否则,执行步骤(4.4j);
(5.4j)将外循环结束时的动态信源数值作为信源数的估计值;
(5.5)第一次估计目标到达角:
(5.5a)按照下式,寻找峰值位置矢量:
其中,pV表示恢复矢量取模值后按降序重排的峰值矢量,pI表示恢复矢量元素取模值降序重排后,恢复矢量原下标值经过重新排列得到的峰值位置矢量,findpeaks(·,′descend′)表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.5b)将峰值位置矢量中的第一个元素值,放入信源位置矢量的第一个位置,将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置,直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值,得到最终信源位置矢量;
(5.6)按照下式,寻找峰值位置矢量:
其中,pV表示恢复矢量取模值后按降序重排的峰值矢量,pI表示恢复矢量元素取模值降序重排后,恢复矢量原下标值经过重新排列得到的峰值位置矢量,findpeaks(·,′descend′)表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.7)按照下式,计算当前代价函数:
其中,F(M)表示当前代价函数,pV(s)表示峰值矢量pV中的第s个元素,s表示动态信源数,∑(·)表示求和,表示恢复矢量中的第i个元素, 表示不属于符号,表示峰值位置矢量pI的前个元素,表示信源数的估计值,|·|表示取模值操作;
(5.8)按照下式,计算当前峰值:
(5.9)计算左移代价函数:
(5.9a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围;
(5.9b)按照下式,计算左移导向矩阵:
A(L)=[a(θ1),a(θ2),…,a(θi),…]
其中,A(L)表示左移导向矩阵,M表示天线阵列中的阵元个数,θi∈Θ(L),∈表示属于符号,Θ(L)表示左移角度搜索范围,i表示搜索角度θi在左移角度搜索范围Θ(L)中的序号;
(5.9c)按照下式,计算左移导向矩阵的正交投影矩阵:
Φ(L)=I-(A(L))H(A(L)(A(L))H)-1A(L)
其中,Φ(L)表示左移导向矩阵A(L)的正交投影矩阵,I表示单位矩阵,A(L)表示左移导向矩阵,H表示共轭转置操作,-1表示求逆操作;
(5.9d)按照下式,计算左移恢复矢量:
(5.9e)按照下式,寻找左移峰值位置矢量:
其中,表示左移恢复矢量取模值后按降序重排的左移峰值矢量,表示左移恢复矢量元素取模值降序重排后,左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量,findpeaks(·,′descend′)表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.9f)按照下式,计算左移代价函数:
其中,F(L)表示左移代价函数,表示左移峰值矢量中的第s个元素,s表示动态信源数,∑(·)表示求和,表示左移恢复矢量中的第i个元素, 表示不属于符号,表示左移峰值位置矢量的前个元素,表示信源数的估计值,|·|表示取模值操作;
(5.10)计算右移代价函数:
(5.10a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围;
(5.10b)按照下式,计算右移导向矩阵:
A(R)=[a(θ1),a(θ2),…,a(θi),…]
其中,A(R)为右移导向矩阵,M表示天线阵列中的阵元个数,θi∈Θ(R),∈表示属于符号,Θ(R)表示右移角度搜索范围,i表示搜索角度θi在右移角度搜索范围Θ(R)中的序号;
(5.10c)按照下式,计算右移导向矩阵的正交投影矩阵:
Φ(R)=I-(A(R))H(A(R)(A(R))H)-1A(R)
其中,Φ(R)表示右移导向矩阵A(R)的正交投影矩阵,I表示单位矩阵,A(R)表示右移导向矩阵,H表示共轭转置操作,-1表示求逆操作;
(5.10d)按照下式,计算右移恢复矢量:
(5.10e)按照下式,寻找右移峰值位置矢量:
其中,表示右移恢复矢量取模值后按降序重排的右移峰值矢量,表示右移峰值矢量元素取模值降序重排后,右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量,findpeaks(·,′descend′)表示寻找局部峰值并将局部峰值按降序排列,|·|表示取模值操作;
(5.10f)按照下式,计算右移代价函数:
其中,F(R)表示右移代价函数,表示右移峰值矢量中的第s个元素,s表示动态信源数,∑(·)表示求和,表示右移恢复矢量中的第i个元素, 表示不属于符号,表示右移峰值位置矢量的前个元素,表示信源数的估计值,|·|表示取模值操作;
(5.11)更新当前栅格参数:
(5.11a)判断左移代价函数是否同时大于当前代价函数和右移代价函数,若是,则执行步骤(5.11b),否则,执行步骤(5.11e);
(5.11b)将步骤(5.7)的当前代价函数的值更新为左移代价函数的值;
(5.11c)将恢复矢量的元素值更新为左移恢复矢量的元素值;
(5.11d)按照下式,更新目标到达角矢量:
(5.11e)判断右移代价函数是否同时大于当前代价函数和左移代价函数,若是,则执行步骤(5.11f),否则,执行步骤(5.12);
(5.11f)将当前代价函数的值更新为右移代价函数的值;
(5.11g)将恢复矢量的元素值更新为右移恢复矢量的元素值;
(5.11h)按照下式,更新目标到达角估计矢量:
(5.12)判断可调栅格步长Δ是否大于栅格优化精度ξθ,若是,则执行步骤(5.9),否则,执行步骤(5.13);
(5.13)按照下式,更新目标到达角估计矢量:
(5.14)判断动态信源数是否小于信源数,若是,则执行步骤(5.6),否则,执行步骤(5.15);
(5.15)获得第二次目标到达角估计矢量:
将步骤(5.13)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910516081.3A CN110174657B (zh) | 2019-06-14 | 2019-06-14 | 基于秩一降维模型和块矩阵恢复的波达方向估计方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910516081.3A CN110174657B (zh) | 2019-06-14 | 2019-06-14 | 基于秩一降维模型和块矩阵恢复的波达方向估计方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110174657A CN110174657A (zh) | 2019-08-27 |
CN110174657B true CN110174657B (zh) | 2022-12-27 |
Family
ID=67698455
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910516081.3A Active CN110174657B (zh) | 2019-06-14 | 2019-06-14 | 基于秩一降维模型和块矩阵恢复的波达方向估计方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110174657B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111814096B (zh) * | 2020-06-28 | 2023-10-20 | 海南大学 | 基于子空间拟合的加权块稀疏恢复的mimo雷达定位方法 |
CN113255579B (zh) * | 2021-06-18 | 2021-09-24 | 上海建工集团股份有限公司 | 一种施工监测异常采集数据自动识别与处理的方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103399291A (zh) * | 2013-07-22 | 2013-11-20 | 西安电子科技大学 | 基于快速稀疏恢复的超分辨波达方向估计方法 |
CN103954950A (zh) * | 2014-04-25 | 2014-07-30 | 西安电子科技大学 | 一种基于样本协方差矩阵稀疏性的波达方向估计方法 |
-
2019
- 2019-06-14 CN CN201910516081.3A patent/CN110174657B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103399291A (zh) * | 2013-07-22 | 2013-11-20 | 西安电子科技大学 | 基于快速稀疏恢复的超分辨波达方向估计方法 |
CN103954950A (zh) * | 2014-04-25 | 2014-07-30 | 西安电子科技大学 | 一种基于样本协方差矩阵稀疏性的波达方向估计方法 |
Non-Patent Citations (1)
Title |
---|
基于协方差矩阵降维稀疏表示的二维波达方向估计;李文杰等;《计算机应用》;20160810(第08期);第2197-2201页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110174657A (zh) | 2019-08-27 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110244272B (zh) | 基于秩一去噪模型的波达方向估计方法 | |
Zhou et al. | Direction-of-arrival estimation for coprime array via virtual array interpolation | |
CN110174658B (zh) | 基于秩一降维模型和矩阵补全的波达方向估计方法 | |
CN112363110B (zh) | 一种基于嵌套交叉偶极子阵列的无网格单比特doa估计方法 | |
CN109116293B (zh) | 一种基于离格稀疏贝叶斯的波达方向估计方法 | |
CN110109050B (zh) | 嵌套阵列下基于稀疏贝叶斯的未知互耦的doa估计方法 | |
CN111337893A (zh) | 一种基于实值稀疏贝叶斯学习的离格doa估计方法 | |
CN107544051A (zh) | 嵌套阵列基于k‑r子空间的波达方向估计方法 | |
Wu et al. | Direction-of-arrival estimation based on Toeplitz covariance matrix reconstruction | |
CN105335615B (zh) | 一种低复杂度的二维角度和极化参数联合估计方法 | |
Bohme | Source-parameter estimation by approximate maximum likelihood and nonlinear regression | |
CN104375133B (zh) | 一种空间二维doa的估算方法 | |
CN111239678A (zh) | 一种基于l型阵列的二维doa估计方法 | |
CN110174657B (zh) | 基于秩一降维模型和块矩阵恢复的波达方向估计方法 | |
KR101958337B1 (ko) | 신호의 도래각을 추정하는 방법 및 장치 | |
CN114720938A (zh) | 基于深度展开的大规模天线阵列单比特采样doa估计方法 | |
CN110941980B (zh) | 一种密集环境中基于压缩感知的多径时延估计方法及装置 | |
Feng et al. | An off-grid iterative reweighted approach to one-bit direction of arrival estimation | |
Sahnoun et al. | Tensor polyadic decomposition for antenna array processing | |
CN107544050A (zh) | 白噪声背景下一种构造自适应阈值估计信号源数目的方法 | |
CN112731280A (zh) | 互质阵列混合噪声环境下的esprit-doa估计方法 | |
Yang et al. | A correlation-aware sparse Bayesian perspective for DOA estimation with off-grid sources | |
CN107135026B (zh) | 未知互耦存在时基于矩阵重构的稳健波束形成方法 | |
CN115421119A (zh) | 一种基于收发翻转互质mimo雷达结构的doa估计方法 | |
CN114844544A (zh) | 一种基于低管秩张量分解的互质阵列波束成形方法、***及介质 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |