基于加权二阶稀疏贝叶斯的相干信源离栅DOA估计方法与流程

基于加权二阶稀疏贝叶斯的相干信源离栅doa估计方法
技术领域
1.本发明属于雷达信号处理领域，具体涉及一种基于加权二阶稀疏贝叶斯的相干信源离栅doa估计方法。


背景技术：

2.波达方向(direction of arrival,doa)估计是阵列信号处理的一个基本问题，也是雷达、声呐等许多领域的重要任务之一。在doa估计技术的发展过程中，人们先后提出了许多的算法，包括传统的doa估计算法(capon算法、前向预测算法、最大熵算法、最小模算法等)、基于子空间的算法和最大似然估计等。其中，子空间类算法的估计精度完全依赖于信号子空间和噪声子空间能否被正确表示，当快拍数减少、信噪比降低以及目标信号之间相关性等情况下，该类方法估计性能会大大降低。
3.随着压缩感知(cs)理论的兴起，针对空域的稀疏性，研究者们将稀疏表示理论应用到doa估计中，并提出许多算法。其中，最具代表性的算法是l1
‑
svd算法。然而上述算法是假设信号源目标角度正好在观测矩阵上的前提下，当存在偏差，即off
‑
grid问题时估计性能就会降低。
4.为了减少网格偏移所带来的误差，学者们对off
‑
grid场景进行了大量的研究，近年来陆续提出了针对off
‑
grid问题的doa估计方法。其中，具有代表性的是文献中yang等人基于贝叶斯思想提出的off
‑
grid稀疏贝叶斯推论(ogsbi)算法，该算法实现了粗网格条件下的高精度doa估计。zhang等人在ogsbi算法基础上，结合music算法，构建信源的先验权矢量，提出了加权稀疏贝叶斯(ogwsbi)算法。较ogsbi算法相比，ogwsbi算法提高了迭代时的收敛速度与估计精度。然而ogwsbi算法在对相干信源进行估计时，会由于music算法的局限性，在空间谱中造成doa估计的漏报，进而造成存在信源的地方不能得到有效加权，降低参数的估计精度。因此，需要构造新的加权系数，为迭代提供一个正确的先验信息。另外，由于导引矢量泰勒展开中，高阶项的缺失会造成参数估计精度的降低，从而需要加入高阶项，进一步提高估计精度。


技术实现要素：

5.为了解决导引矢量泰勒展开中高阶项的缺失及加权稀疏贝叶斯算法运用在相干信源时对doa估算精度不高的问题，本发明提供了一种基于加权二阶稀疏贝叶斯的相干信源离栅doa估计方法。该算法在ogwsbi算法基础上，首先为了降低导引矢量一阶泰勒展开而带来的模型误差，将导引矢量扩展至二阶泰勒展开；然后利用酉变换既对接收信号进行解相干处理，恢复协方差矩阵的秩，又将估计模型由复数域转换为实数域，提高算法的估计精度以及运算效率。
6.本发明为实现上述目的，采用如下技术方案：
7.一种基于加权二阶稀疏贝叶斯的相干信源离栅doa估计方法，所述doa估计方法用于以下模型：接收信号为其稀疏模型表达式为y
c
＝φ(β)s n,其中为信
号源，β为角度偏离量，表示噪声，m为阵元数，n为网格划分个数，t为快拍数，k为信号源个数，φ为由导引矢量构成的过完备字典；
8.所述doa估计方法包括如下步骤：
9.步骤s1：对接收信号y
c
酉变换实数域转换及信号解相干处理，从而得到酉变换后的接收信号y；
10.步骤s2：对酉变换后的接收信号y奇异值分解，得到y
sv
及s
sv
，结合music算法构造加权矢量w；
11.步骤s3：对由二阶导引矢量构建的过完备字典的参数矩阵酉变换并结合加权稀疏贝叶斯推论求解超参数α0、α及偏移量β，并输出doa估计值
12.进一步地，步骤s3包含以下步骤：
13.步骤s31：利用二阶泰勒展开的真实信源到达角度导引矢量构建过完备字典：
14.其中，a
c
为观测矩阵：为观测矩阵：是分布在[0,π]上的空间角度离散采样值，n为网格划分个数；二阶泰勒展开的导引矢量为的空间角度离散采样值，n为网格划分个数；二阶泰勒展开的导引矢量为为与真实信号角度θ
k
最接近的网格值，为a(θ)在处的一阶导数，为a(θ)在处的二阶导数，由组成的矩阵为成的矩阵为组成的矩阵
[0015]
步骤s32：对过完备字典中的a
c
、b
c
及c
c
进行酉变换得到酉变换后的过完备字典
[0016]
步骤s33：用λ＝diag(α)、p(α
n
)＝γ(α|1,w
n
ρ)，ρ为用户选择性参数及信号源s
sv
的后验概率分布的后验概率分布得到信号的均值以及方差μ(t)＝α0∑φ
h
y
sv
(t),t＝1,....k，∑＝(α0φ
h
φ λ
‑1)
‑1，利用酉变换后的过完备字典φ更新信号的均值μ(t)以及方差∑，并基于信号的均值μ(t)以及方差∑更新超参数α0、α及偏移量β；判断是否达到最大迭代次数i
max
或者||α
i 1
‑
α
i
||2/||α
i2
＜τ，τ为容限参数，若两者均不满足，利用更新后偏移量β更新过完备字典φ，并利用更新后的过完备字典φ继续更新信号的均值μ(t)以及方差∑从而更新超参数α0、α及偏移量β；若满足其中一个或者两个条件，则结束迭代，并输出doa估计值
[0017]
进一步地，更新后的超参数α0、α分别为：
[0018][0019]
[0020]
其中，μ＝[μ(1),...μ(k)]＝α0∑φ
h
y
sv
，ρ＝ρ/k，b,c
→
0。
[0021]
进一步地，更新后的偏移量β为：
[0022][0023]
其中，
[0024][0025][0026]
进一步地，步骤s1包括如下步骤：
[0027]
步骤s11：对接收信号y
c
构建增广矩阵形成中心埃尔米特矩阵：
[0028]
其中，y
c*
为y
c
的复共轭，
[0029]
步骤s12：对中心埃尔米特矩阵y
aug
酉变换处理：其中，y表示y
aug
经过酉变换后，从复数域变换为实数域的矩阵，u为酉矩阵；
[0030]
当m为偶数时，u
m
可表示为：
[0031]
其中，i为m/2
×
m/2的单位矩阵，j为m/2
×
m/2的交换矩阵，其副对角线为1，其他元素为0；
[0032]
m为奇数时，u
m
可表示为：
[0033][0034]
其中，i和j的维数为(m
‑
1)/2；
[0035]
同理，u
2t
为2t维酉矩阵。
[0036]
进一步地，步骤s32中对过完备字典中的a
c
、b
c
及c
c
进行酉变换的a、b及c分别如下所示：
[0037][0038][0039][0040]
进一步地，步骤s2包含以下步骤：
[0041]
对酉变换后的接收信号y进行奇异值分解，即y＝u
s
∑
s
v
st
u
e
∑
e
v
et
,得到信号子空间u
s
、噪声子空间u
e
；令；令为实数化后的信号源，得到y
sv
、s
sv
；
[0042]
则倒music空间谱公式为:记权值w＝[w1,
…
,w
n
]
t
,同时对权值进行归一化得到加权矢量w：w＝w/min(w)。
[0043]
有益效果：利用酉变换对接收信号进行解相干处理，恢复协方差矩阵的秩；同时为了降低导引矢量一阶泰勒展开而带来的模型误差，将导引矢量扩展至二阶泰勒展开，并对导引矢量构成的过完备字典进行酉变换，从而实现接收信号稀疏模型y
c
＝φ(β)s n的实数化，最后结合加权稀疏贝叶斯推论求解偏移量，本发明的doa算法具有估计精度以及运算效率高的特点。
附图说明
[0044]
图1为本发明一实施例中基于加权二阶稀疏贝叶斯的相干信源离栅doa估计方法(iogwsbi)的流程图；
[0045]
图2为导引矢量一阶泰勒展开模型条件下，ogwsbi算法在相干以及非相干信源下角度估计剖面图；
[0046]
图3为酉变换处理前后music算法相干信源角度估计剖面图；
[0047]
图4为ogwsbi算法与iogwsbi的相干信源角度估计剖面图；
[0048]
图5为ogwsbi算法与iogwsbi算法doa估计均方根误差与信噪比关系图；
[0049]
图6为ogwsbi算法与iogwsbi算法doa估计均方根误差与快拍数关系图；
[0050]
图7为ogwsbi算法与iogwsbi算法运行时间与信噪比关系图。
具体实施方式
[0051]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0052]
如图1所示，一种基于加权二阶稀疏贝叶斯的相干信源离栅doa估计方法，所述doa估计方法用于以下模型：假设k个窄带远场信号以角度[θ1,θ2,...θ
k
]
t
入射到由m个全向接收阵元组成的均匀线阵上。入射信号的波长为λ，相邻阵元间距为d＝λ/2。是分布在[0,π]上的空间角度离散采样值，n为网格划分个数。当出现off
‑
grid情况时，由于入射信号的角度θ
k
没有完全落在观测矩阵的网格上，因此与观测矩阵上的角度存在偏移量β。接收信号为其稀疏模型表达式为y
c
＝φ(β)s n,其中为信号源，表示噪声，m为阵元数，n为网格划分个数，t为快拍数，k为信号源个数，φ为由导引矢量构成的过完备字典；现有技术中一般导引矢量采用一阶泰勒级数近似；
[0053]
所述doa估计方法包括如下步骤：
[0054]
步骤s1：对接收信号y
c
酉变换实数域转换及信号解相干处理，从而得到酉变换后的接收信号y；
[0055]
在步骤s1中接收信号y
c
进行酉变换时此时的偏移量β取值为0，具体酉变换过程如下：
[0056]
步骤s11：对接收信号y
c
构建增广矩阵形成中心埃尔米特矩阵：
[0057]
其中，y
c*
为y
c
的复共轭，
[0058]
步骤s12：对中心埃尔米特矩阵y
aug
酉变换处理：其中，y表示y
aug
经过酉变换后，从复数域变换为实数域的矩阵，u为酉矩阵；
[0059]
当m为偶数时，u
m
可表示为：其中，i为m/2
×
m/2的单位矩阵，j为m/2
×
m/2的交换矩阵，其副对角线为1，其他元素为0；
[0060]
m为奇数时，u
m
可表示为：
[0061][0062]
其中，i和j的维数为(m
‑
1)/2,i为m/2
×
m/2的单位矩阵，j为m/2
×
m/2的交换矩阵，其副对角线为1，其他元素为0。
[0063]
同理，u
2t
为2t维酉矩阵。
[0064]
由于酉变换还可以恢复接收信号协方差矩阵的秩，使其与信号源个数k相等，即rank(r)＝k，进而实现解相干性能。
[0065]
步骤s2：对酉变换后的接收信号y奇异值分解，得到y
sv
及s
sv
，结合music算法构造加权矢量w；；
[0066]
具体步骤如下：对酉变换后的接收信号y进行奇异值分解，即y＝u
s
∑
s
v
st
u
e
∑
e
v
et
,得到信号子空间u
s
、噪声子空间u
e
；令；令；令为实数化后的信号源，得到y
sv
、s
sv
；
[0067]
则倒music空间谱公式为:记权值w＝[w1,
…
,w
n
]
t
,同时对权值进行归一化得到加权矢量w：w＝w/min(w)。
[0068]
步骤s3：对由二阶导引矢量构建的过完备字典的参数矩阵酉变换并结合加权稀疏贝叶斯推论求解超参数α0、α及偏移量β，并输出doa估计值步骤s3包含以下步骤：
[0069]
步骤s31：利用二阶泰勒展开的真实信源到达角度导引矢量构建过完备字典：
[0070]
其中，a
c
为观测矩阵：为观测矩阵：是分布在[0,π]上的空间角度离散采样值，n为网格划分个数；二阶泰勒展开的导引矢量为的空间角度离散采样值，n为网格划分个数；二阶泰勒展开的导引矢量为为与真实信号角度θ
k
最接近的网格值，为a(θ)在处的一阶导数，为a(θ)在处的二阶导数，由组成
的矩阵为的矩阵为组成的矩阵
[0071]
步骤s32：对过完备字典中的a
c
、b
c
及c
c
进行酉变换得到酉变换后的过完备字典
[0072]
对于均匀线阵，导引矢量满足中心埃尔米特特性，步骤s32中对过完备字典中的a
c
、b
c
及c
c
进行酉变换的a、b及c分别如下所示：
[0073][0074][0075][0076]
通过对接收信号以及观测矩阵进行酉变换，使得整个估计模型由复数域转换为实数域，从而有效降低计算量。
[0077]
步骤s33：用λ＝diag(α)、p(α
n
)＝γ(α|1,w
n
ρ)，ρ为用户选择性参数及信号源s
sv
的后验概率分布的后验概率分布得到信号的均值以及方差μ(t)＝α0∑φ
h
y
sv
(t),t＝1,....k，∑＝(α0φ
h
φ λ
‑1)
‑1，利用酉变换后的过完备字典φ更新信号的均值μ(t)以及方差∑，并基于信号的均值μ(t)以及方差∑更新超参数α0、α及偏移量β；判断是否达到最大迭代次数i
max
或者||α
i 1
‑
α
i
||2/||α
i
||2＜τ，τ为容限参数，若两者均不满足，利用更新后偏移量β更新酉变换后的过完备字典φ，并利用更新后的过完备字典φ继续更新信号的均值μ(t)以及方差∑从而更新超参数α0、α及偏移量β；若满足其中一个或者两个条件，则结束迭代，并输出doa估计值
[0078]
在步骤s33中，经过步骤s2获得加权矢量w后需要将其融合到稀疏贝叶斯框架中。根据稀疏贝叶斯算法对于信号分布，假设其满足：
[0079]
其中，s
sv
(t)表示s
sv
的第t列，对角矩阵λ＝diag(α)，α＝[α1,...,α
n
]
t
，α
n
控制信号源的稀疏性，即只有存在信源的方向上α
n
为非零数。将权矢量加入到其概率密度假设中，假设其满足gamma分布：p(a
n
)＝γ(a|1,w
n
ρ)。其中，γ(
·
)为gamma函数，ρ为用户选择性参数，w
n
为第n个权值系数。由信号源s
sv
的后验概率分布得到信号的均值以及协方差满足：
[0080]
μ(t)＝α0∑φ
h
y
sv
(t),t＝1,....k
ꢀꢀꢀ
(2)
[0081]
∑＝(α0φ
h
φ λ
‑1)
‑1ꢀꢀꢀꢀꢀ
(3)
[0082]
更新后的超参数α0、α分别为：
[0083]
[0084][0085]
公式(4)和公式(5)中，μ＝[μ(1),...μ(k)]＝α0∑φ
h
y
sv
，，ρ＝ρ/k，b,c
→
0。
[0086]
求解角度偏离量β的方法为，通过em算法，使e{logp(y
sv
|s
sv
,α0,β)p(β)}最大，使下式最小：
[0087][0088]
则偏移量β的更新公式通过计算得到：
[0089][0090]
其中，其中，
[0091][0092]
为了进一步简化偏移量运算，令为变换矩阵，是第i行为1，其它行为0的n
×
1列向量。β
‑
k
表示偏移量β中除第k个元素以外的部分，令对β
new
进行变换得到：
[0093][0094]
其中，d1为式β
t
p1β中与β
k
无关的部分。类似的：
[0095]
v
1t
β＝(v
1t
)
k
β
k
d2
[0096][0097]
[0098][0099]
其中，d2、d3、d4、d5都是与β
k
无关的部分，则偏移量β
k
为下式的实值根：
[0100][0101]
另外，由于偏离量β与信源具有相同的稀疏结构，因此，p和v可以降维至k或者k
×
k维。公式(7)中f
′
(β
k
)所得根为一次迭代下的偏移量，而在每次迭代过程中，需要计算均值、协方差矩阵以及更新超参数α0,α,β，且终止迭代条件为达到最大迭代次数i
max
或者||α
i 1
‑
α
i
||2/||α
i
||2＜τ。
[0102]
超参数α0、α及偏移量β的更新过程如下：
[0103]
(1)输入初始值：β＝0，ρ＝0.01，c＝d＝1
×
10
‑4,,
[0104]
(2)根据公式(1)得到酉变换后的过完备字典φ，根据公式(2)及(3)得到信号源的均值μ(t)以及协方差∑；
[0105]
(3)利用均值μ(t)以及协方差∑结合公式(4)及公式(5)得到更新后的超参数α0、α，利用公式(7)更新偏离量β；
[0106]
(4)判断是否达到最大迭代次数i
max
或者||α
i 1
‑
α
i
||2/||α
i
||2＜τ，若两者均不满足，利用更新后偏移量β根据公式(1)更新酉变换后的过完备字典φ，即转至步骤(1)；若满足其中一个或者两个条件，则结束迭代，并输出doa估计值足其中一个或者两个条件，则结束迭代，并输出doa估计值
[0107]
其中，表示第j个信号源的网格值，β
′
j
为第j个信号源的网格偏离量。
[0108]
接下来以具体实施例进行说明，假设k＝2个窄带远场信号以角度[θ1,θ2,...θ
k
]
t
入射到由m＝8个全向接收阵元组成的均匀线阵上。入射信号的波长为λ，相邻阵元间距为d＝λ/2，是分布在[0,π]上的空间角度离散采样值，n为网格划分个数，设置n＝91，迭代次数i
max
＝2000，容限参数τ＝10
‑3，具体参数如表1所示：
[0109]
表1系统仿真参数
[0110]
参数名称参数数值阵元数(m)8阵元间隔(d)λ(波长)/2角度范围0
°
～180
°
角度间隔(r)2
°
信源个数(k)2
[0111]
图2是ogwsbi算法在相干以及非相干信源的角度剖面图，从图2中可以得到，当信号源是非相干信源时，ogwsbi算法可以精准估计出信号来向；当信号源是相干信源时，估计性能迅速降低。
[0112]
图3是music算法在酉变换前后估计相干信源的角度剖面图。从图3中可以看出，在经过酉变换处理后，协方差的秩恢复为rank(r)＝k，信号子空间和噪声子空间得以正确表示，music算法能够高精度估计出相干信源，从而可以提供正确的加权矢量。
[0113]
本次发明基于以上导引矢量二阶泰勒展开离栅模型以及估计模型的酉变换对ogwsbi算法加以改进，图4所示为ogwsbi与本文所提出算法iogwsbi的在snr＝10db下一次蒙特卡诺实验空间谱对比图，表2为其相对应的结果。
[0114]
从图4和表2可以发现，iogwsbi的估计性能要明显优于ogwsbi算法。这表明，在信号源相干时，iogwsbi算法由于酉变换去相关处理，仍能保持优越的doa估计性能。
[0115]
表2 doa结算结果
[0116] θ1θ2true signal64.8
°
86.5
°
ogwsbi65.06
°
87
°
iogwsbi64.88
°
86.55
°
[0117]
图5是在快拍数l＝100下，doa估计均方根误差随信噪比变化的曲线。图6是在snr＝10db下，doa估计均方根误差随快拍数变化的曲线。由图5以及图6可以看出，两个算法随着信噪比以及快拍数的增加，均方根误差逐渐降低。但由于iogwsbi算法进行了去相关处理，使得加权矢量得以正确表示。因此，在同等信噪比或者快拍数的条件下，所提算法要优于ogwsbi算法，且估计性能提高约50％。
[0118]
iogwsbi算法由酉变换、奇异值分解和加权稀疏重构3部分组成。酉变换计算量由于其只涉及加减运算而忽略不计。奇异值分解的计算复杂度为o(max(m2t,mt2))。可以发现，不管是酉变换还是奇异值分解计算都只需要计算一次，而本算法是迭代类算法，算法复杂度取决于加权稀疏重构部分。ogwsbi算法单次迭代的计算复杂度为o(max(mn2)),且在相干信源下，由于错误的加权信息而导致迭代次数的增加，进一步增加了计算复杂度。iogwsbi由于采用酉变换，一次实数域的运算仅为之前复数域的1/4，且迭代次数由于正确的加权信息而大大减少。另外，iogwsbi算法由于采用导引矢量二阶泰勒展开模型，运算复杂度较一阶模型有一定的提高，但整体运算速度仍优于ogwsbi算法。
[0119]
图7是在快拍数l＝100下，算法运行时间随信噪比变化的曲线。从图7中可以看出，当信噪比在0～20db变化时，随着信噪比的提高，两个算法的运算时间不断降低，而iogwsbi算法运行效率较ogwsbi提高约60％。因此，可以得出，在相干信源条件下，iogwsbi算法是一个高效、高精度的离栅doa估计算法，更加适用于对精度要求更高的场景。
[0120]
本发明在ogwsbi算法基础上，通过对估计模型采取酉变换处理，不仅能对相干信源进行去相干处理，解决了ogwsbi算法在相干信源下估计性能不足的问题，还可以使算法从复数域转换为实数域进行运算，显著提高了算法运行效率。另外，为了进一步提高估计精度，还在导引矢量一阶泰勒展开的基础上进行二阶泰勒展开，从而减少了信号的拟合误差。本发明所提算法在精度以及效率上进行了有效的改进，拓宽了加权稀疏贝叶斯算法的实际应用场景。
[0121]
以上所述仅为本发明中的具体实施方式，并不局限于此，任何在本发明精神和原则之内的，所做的变换或替换，均包含在本发明的保护范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。