CN113946955B - 基于融合中心反馈信息的多目标贝叶斯波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法，包括：根据融合中心的反馈信息，确定多个目标波达方向的先验分布信息；根据阵列接收信号和所述先验分布信息，构造最小均方误差估计量；对所述最小均方误差估计量进行离散近似，以得到近似积分；基于蒙特卡洛法对先验分布信息对应的随机数和所述近似积分进行计算，以确定每个目标波达方向的估计值。本发明能够在低信噪比、少快拍数情况下，提高多目标波达方向估计性能。

Description

基于融合中心反馈信息的多目标贝叶斯波达方向估计方法

技术领域

本发明属于信号处理领域，具体涉及一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法。

背景技术

波达方向估计在雷达、声纳和无线通信中，有很广泛的应用。

目前，在传统的波达方向估计方法中常规波束形成法的分辨率通常会受到瑞利限的限制，为了克服该问题，提出了子空间类方法，如MUSIC算法。尽管子空间类法能够突破瑞利限的限制，能够在极大的提高目标的分辨率的同时保证估计性能，但是该类方法在进行波达方向估计的过程中需要计算阵列接收信号的相关矩阵。为了准确计算相关矩阵，需要大量快拍数，因此，计算过程受到快拍数的影响极大，尤其是当快拍数较少或者接收信号相干时，阵列接收信号的相关矩阵会发生秩亏，从而严重影响子空间类算法的估计性能，甚至无法完成估计。此外，子空间类算法的估计性能受信噪比影响极大，当信噪比降低时，估计性能也会急剧下降。

进一步的，传统的波达方向估计方法在建模时，均假定目标的波达方向是一个确定的值，需要通过阵列接收信号进行估计，导致波达方向估计方法受到极大限制。

另外，针对多目标的波达方向估计问题，现有技术中通常采用两种估计方法：一种是最大后验估计法，另一种是传统的最小均方误差估计法。其中，最大后验估计法会遇到多维最大值求解问题，其算法复杂度会随着目标数量的增加而呈指数增加；传统的最小均方误差估计法中包含多重积分，其解析求解过程复杂。因此，现有技术中针对多目标的波达方向估计方法存在耗时过长、估计准确度低等问题。

发明内容

为了解决现有技术中存在的上述问题，本发明提供了一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法。本发明要解决的技术问题通过以下技术方案实现：

一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法，应用于传感器网络，所述传感器网络中包括融合中心和均匀线列阵，所述方法包括：步骤1：根据融合中心的反馈信息，确定多个目标波达方向的先验分布信息；步骤2：根据阵列接收信号和所述先验分布信息，构造最小均方误差估计量；步骤3：对所述最小均方误差估计量进行离散近似，以得到近似积分；步骤4：基于蒙特卡洛法对所述先验分布信息对应的随机数和所述近似积分进行计算，以确定每个目标波达方向的估计值。

在本发明的一个实施例中，所述步骤1包括：步骤1-1：获取融合中心的反馈信息，其中，所述反馈信息中包括K个目标对应的位置预测值及其协方差矩阵Q_ui，所述x、y表示直角坐标系中的x坐标轴、y坐标轴；步骤1-2：假设均匀线列阵与x坐标轴平行，根据每个目标对应的所述位置预测值/>计算目标的波达方向的预测值θ_i ^o，以及根据每个目标对应的所述协方差矩阵Q_ui计算目标的波达方向的预测方差/>表示为：

其中，b_i为参数，表示为：

s_x和s_y分别表示阵列几何中心的位置坐标；步骤1-3：假设目标波达方向满足高斯分布，根据所述预测值θ_i ^o和所述预测方差确定目标波达方向的先验分布信息，表示为：

其中，θ表示波达方向，θ＝[θ₁,θ₂，…，θ_K]^T， diag(·)表示由向量形成的对角阵。

在本发明的一个实施例中，所述步骤2包括：步骤2-1：设所述均匀线列阵中包括N个阵元，将阵列接收信号x(t)表示为：

x(t)＝A(θ)s(t)+n(t)，

其中，t＝1,2,…,M，M表示阵列快拍数；s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T，表示t时刻的第i个来波信号，a_ti和φ_ti表示t时刻的第i个来波信号的幅度和相位，j表示虚数单位；n(t)～CN(0,Q_n)表示均值为零且协方差矩阵为Q_n的复值圆对称高斯白噪声；A(θ)表示阵列流形向量，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]，其中，a(θ_i)表示为：

a(θ_i)的第k个元素代表第i个信号在第k个阵元的时间延迟所带来的相位变化，d表示阵元间距，λ表示接收信号的波长；步骤2-2：令变量且将阵列接收信号x(t)改写为：/>基于上式，获取目标波达方向的似然函数，表示为：

其中，步骤2-3：基于贝叶斯公式，对所述先验分布信息进行计算，以得到后验分布信息/>表示为：

则得到最小均方误差估计量，表示为：

本发明的有益效果：

本发明能够根据融合中心的反馈信息进行多目标贝叶斯波达方向估计，具体的，利用先验分布信息和阵列接收信号构造最小均方误差估计量，能够将被估计量视为需估计出其现实的随机变量，在估计时引入被估计量的先验分布信息，并且得到后验概率密度函数，最后得到波达方向的最小均方误差统计量，并通过蒙特卡洛方法进行数值求解。

本发明相比于最大后验估计方法，能够避免多维搜索；相比于传统最小均方误差估计方法，能够基于蒙特卡洛方法降低求解复杂度，从而提高目标波达方向估计性能。另外，本发明还能够在低信噪比、少快拍数情况下，进一步提高多目标波达方向估计性能。

以下将结合附图及实施例对本发明做进一步详细说明。

附图说明

图1是本发明实施例提供的一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法流程示意图；

图2是本发明实施例提供的另一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法流程示意图；

图3是本发明实施例提供的在信噪比为15dB时的仿真实验估计结果示意图；

图4是本发明实施例提供的在信噪比为-15dB时的仿真实验估计结果示意图；

图5是本发明实施例提供的在信噪比为-15dB时的10000次估计结果误差统计直方图；

图6是本发明实施例提供的在不同信噪比情况下的单目标波达方向估计性能结果示意图。

具体实施方式

下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。

实施例

请参见图1，图1是本发明实施例提供的一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法示意图，应用于传感器网络，所述传感器网络中包括融合中心和均匀线列阵，所述方法包括：

步骤1：根据融合中心的反馈信息，确定多个目标波达方向的先验分布信息。

可选的，所述步骤1包括：

步骤1-1：获取融合中心的反馈信息，其中，所述反馈信息中包括K个目标对应的位置预测值及其协方差矩阵Q_ui，所述x、y表示直角坐标系中的x坐标轴、y坐标轴。

步骤1-2：假设均匀线列阵与x坐标轴平行，根据每个目标对应的所述位置预测值计算目标的波达方向的预测值θ_io，以及根据每个目标对应的所述协方差矩阵Q_ui计算目标的波达方向的预测方差/>表示为：

其中，b_i为参数，表示为：

s_x和s_y分别表示阵列几何中心的位置坐标。

步骤1-3：假设目标波达方向满足高斯分布，根据所述预测值θ_i ^o和所述预测方差确定目标波达方向的先验分布信息，表示为：

其中，θ表示波达方向，θ＝[θ₁,θ₂，…，θ_K]^T， diag(·)表示由向量形成的对角阵。

步骤2：根据阵列接收信号和所述先验分布信息，构造最小均方误差估计量。

可选的，所述步骤2包括：

步骤2-1：设所述均匀线列阵中包括N个阵元，将阵列接收信号x(t)表示为：

x(t)＝A(θ)s(t)+n(t) (5)

其中，t＝1,2,…,M，M表示阵列快拍数；s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T，表示t时刻的第i个来波信号，a_ti和φ_ti表示t时刻的第i个来波信号的幅度和相位，j表示虚数单位；n(t)～CN(0,Q_n)表示均值为零且协方差矩阵为Q_n的复值圆对称高斯白噪声；

A(θ)表示阵列流形向量，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]，其中，a(θ_i)表示为：

a(θ_i)的第k个元素代表第i个信号在第k个阵元的时间延迟所带来的相位变化，d表示阵元间距，λ表示接收信号的波长。

可选的，一般取d＝λ/2。

步骤2-2：令变量 且/>将阵列接收信号x(t)改写为：

基于上式，获取目标波达方向的似然函数，表示为：

其中，

步骤2-3：基于贝叶斯公式，对所述先验分布信息进行计算，以得到后验分布信息表示为：

结合θ的先验分布p(θ)，利用贝叶斯公式可以计算得到其后验分布

进一步，则得到最小均方误差估计量，表示为：

可以看出，该估计量是关于未知量的非线性函数的复杂积分，很难进行闭式解析求解，因此，本发明通过以下步骤3将对该积分进行离散近似后，通过蒙特卡洛方法给出数值解。

步骤3：对所述最小均方误差估计量进行离散近似，以得到近似积分。

可选的，所述步骤3包括：

步骤3-1：将所述最小均方误差估计量中的分子看作变量关于θ的先验分布取期望，记为/>以及将所述最小均方误差估计量中的分母看作变量关于θ的先验分布取期望，记为/>

则将所述最小均方误差估计量改写为：

步骤3-2：对所述改写后的最小均方误差估计量进行离散近似求平均，以得到近似积分，表示为：

其中，N_m表示离散求和的点数。

通常，点数越大，近似越精确。

步骤4：根据蒙特卡洛法对所述先验分布信息对应的随机数和所述近似积分进行计算，以确定每个目标波达方向的估计值。

可选的，所述步骤4包括：

步骤4-1：令：

解得：

步骤4-2：将步骤4-1解得的带入/>再将该/>代入所述近似积分中，化简可得：

其中，是仅关于θ^(k)的函数，θ^(k)表示样本。

由于函数也包含了自变量/>因此不能直接计算，需要将/>表示为θ^(k)的函数。

步骤4-3：根据所述先验分布信息的概率密度函数，得到N_m个样本。

步骤4-4：将所述N_m个样本代入步骤4-2中的公式，计算得到目标波达方向的估计值

综上，本发明的有益效果：

本发明能够根据融合中心的反馈信息进行多目标贝叶斯波达方向估计方法，具体的，利用先验分布信息和阵列接收信号构造最小均方误差估计量，能够将被估计量视为需估计出其现实的随机变量，在估计时引入被估计量的先验分布信息，并且得到后验概率密度函数，最后得到波达方向的最小均方误差统计量，并通过蒙特卡洛方法进行数值求解。

本发明相比于最大后验估计方法能够避免多维搜索，相比于传统最小均方误差估计方法，能够基于蒙特卡洛方法降低降低求解复杂度。另外，本发明能够在低信噪比、少快拍数情况下，提高多目标波达方向估计性能。

进一步地，基于仿真实验验证本发明的有益效果：

实验条件：

设置均匀线列阵中有10个阵元，取阵元间距d为半波长值，即假设一共有两个目标且信号频率均为f＝1000Hz，采样频率fs＝8192Hz，目标所处坐标位置通过式(1)换算可得两个信号的波达方向分别为/>度，/>度，假设通过式(2)计算两个波达方向的预测方差分别为σ_θ1＝2度，σ_θ2＝1度。

实验过程：

参见图3，图3中快拍数为30，信噪比设置为15dB。

将阵列接收信号x(t)按照式(7)构造然后根据式(4)的先验分布信息产生N_m＝800个样本θ^(k)，k＝1,2,…,800。然后将构造的/>和样本θ^(k)，k＝1,2,…,800带入公式(14)进行计算/>那么/>就是第一个目标的波达方向，/>就是第二个目标的波达方向。

从图3可以看出MUSIC方法和本发明的提出的MMSE(最小均方误差，Minimum MeanSquared Error)估计方法都可以准确的估计出两个目标方向。图4重复了图3的仿真实验，只是信噪比设置为-15dB，远低于图3的信噪比。可以看到此时图4中的MUSIC算法仅能估计一个目标，因此已经失效，而本发明的方法依然可以估计出两个目标位置。验证了本发明的方法在低信噪比时的优异性能。

图5分别给出已有MUSIC算法和本发明方法的10000次仿真实验估计误差结果的统计直方图。可以看到在同样条件下，总体而言本发明的MMSE方法的估计误差更集中在0度附近，意味着估计结果更集中在真实角度附近。而MUSIC算法的估计误差结果尽管在0度附近有所集中，但是其估计误差分布广泛，比如在-60°到60°范围外依然还存在，这意味着MUSIC方法的估计结果的可靠性不高。特别是针对目标二，MUSIC算法在0度附近的次数少于1000次，而本发明的方法误差在0度附近的次数约为1100次。因此可以直观看出来，本发明的MMSE方法的估计精度和估计可靠度优于MUSIC算法。

图6对比MUSIC算法，给出了MMSE方法的单目标波达方向估计性能结果，性能指标由估计的均方误差(RMSE)衡量，指标定义如下：

式中，θ_l表示第l次实验的估计位置，L＝500表示总实验次数。信噪比从-30dB变为20dB。图6中快拍数设置为10，两种方法在信噪比较低(-5dB附近)时，MUSIC算法的估计结果已经远离BCRB(贝叶斯克拉美罗下界，Bayesian Cramer-Rao Bound)，而本发明提出的MMSE方法在该信噪比依然可以达到BCRB。验证了本方法在低信噪比条件下的优越性能。

需要说明的是，如在仿真实验中，随机数根据仿真软件产生标准的正态分布(期望为0，方差为1)，然后变换到给定期望和方差的正态分布。所述仿真软件如matlab。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (2)

1.一种基于融合中心反馈信息的多目标贝叶斯波达方向估计方法，其特征在于，应用于传感器网络，所述传感器网络中包括融合中心和均匀线列阵，所述方法包括：

步骤1：根据融合中心的反馈信息，确定多个目标波达方向的先验分布信息；

步骤2：根据阵列接收信号和所述先验分布信息，构造最小均方误差估计量；

步骤3：对所述最小均方误差估计量进行离散近似，以得到近似积分；

步骤4：根据蒙特卡洛法对所述先验分布信息对应的随机数和所述近似积分进行计算，以确定每个目标波达方向的估计值；

所述步骤1包括：

步骤1-1：获取融合中心的反馈信息，其中，所述反馈信息中包括K个目标对应的位置预测值i＝1,2,…,K及其协方差矩阵Q_ui，所述x、y表示直角坐标系中的x坐标轴、y坐标轴；

步骤1-2：假设均匀线列阵与x坐标轴平行，根据每个目标对应的所述位置预测值i＝1,2,…,K计算目标的波达方向的预测值θ_i ^o，以及根据每个目标对应的所述协方差矩阵Q_ui计算目标的波达方向的预测方差/>表示为：

其中，b_i为参数，表示为：

s_x和s_y分别表示阵列几何中心的位置坐标；

步骤1-3：假设目标波达方向满足高斯分布，根据所述预测值θ_i ^o和所述预测方差确定目标波达方向的先验分布信息，表示为：

其中，θ表示波达方向，θ＝[θ₁,θ₂，…，θ_K]^T， diag(·)表示由向量形成的对角阵；

所述步骤2包括：

步骤2-1：设所述均匀线列阵中包括N个阵元，将阵列接收信号x(t)表示为：

x(t)＝A(θ)s(t)+n(t)，

其中，t＝1,2,…,M，M表示阵列快拍数；s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T，表示t时刻的第i个来波信号，a_ti和φ_ti表示t时刻的第i个来波信号的幅度和相位，j表示虚数单位；n(t)～CN(0,Q_n)表示均值为零且协方差矩阵为Q_n的复值圆对称高斯白噪声；

A(θ)表示阵列流形向量，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]，其中，a(θ_i)表示为：

a(θ_i)的第k个元素代表第i个信号在第k个阵元的时间延迟所带来的相位变化，d表示阵元间距，λ表示接收信号的波长；

步骤2-2：令变量 且/>将阵列接收信号x(t)改写为：

基于上式，获取目标波达方向的似然函数，表示为：

其中，

步骤2-3：基于贝叶斯公式，对所述先验分布信息进行计算，以得到后验分布信息表示为：

则得到最小均方误差估计量，表示为：

所述步骤3包括：

步骤3-1：将所述最小均方误差估计量中的分子看作变量关于θ的先验分布取期望，记为/>以及将所述最小均方误差估计量中的分母看作变量/>关于θ的先验分布取期望，记为/>则将所述最小均方误差估计量改写为：

步骤3-2：对所述改写后的最小均方误差估计量进行离散近似求平均，以得到近似积分，表示为：

其中，N_m表示离散求和的点数。

2.根据权利要求1所述的方法，其特征在于，所述步骤4包括：

步骤4-1：令：

解得：

步骤4-2：将步骤4-1解得的带入/>再将该/>代入所述近似积分中，化简可得：

其中，是仅关于θ^(k)的函数，θ^(k)表示样本；

步骤4-3：根据所述先验分布信息的概率密度函数，得到N_m个样本；

步骤4-4：将所述N_m个样本代入步骤4-2中的公式，计算得到目标波达方向的估计值