本发明公开了一种基于指向性模型数值求解的远场窄带信号来波方向估计方法，该方法包括：第一步：获取M通道阵元数据x1(n),x2(n),…,xM(n)；第二步：根据给定的阵列参数自动生成等余弦间隔的预引导波束扫描角θk和总的预引导波束个数K；第三步：对已知频率f0的远场窄带信号进行频域波束形成，并求得该已知频率远场窄带信号的各预引导波束扫描角所对应的波束功率谱P(θk)；第四步：搜索波束功率谱P(θk)最大值所对应的预引导波束号kp；第五步：利用二分数值求解法求解信号来波方向相对偏差第六步，利用估计远场窄带信号的来波方向。该算法基于常规波束形成，实现简单，运算量小，且估计精度高，适用于实时的工程应用场合。
一种基于指向性模型数值求解的远场窄带信号来波方向估计
方法
技术领域
[0001] 本发明属于信号处理领域，尤其涉及一种基于指向性模型数值求解的远场窄带信号来波方向估计方法。
背景技术
[0002] 目标来波方向估计是雷达、声纳、声学、语音和无线通信等阵列信号处理应用中的一个基本问题，尤其在水声和电子侦察处理中扮演了极其重要的角色。此外，在现代移动通信系统中，来波方向估计在提高无线系统性方面有着重要应用。
[0003] 目前国内外学者提出了很多目标来波方向估计方法，如最小方差无失真响应、基于特征向量的来波方向估计算法(如MUSIC、Root‑MUSIC、RV‑Root‑MUSIC)。这些算法通常能达到很高的精度。然而，这些基于峰值搜索的算法需要很大的计算复杂度。此外还存在着估计性能低、复杂度高和相位模糊等问题。同时，这些算法对信号失配非常敏感。
[0004] 在波束形成方法中，常规波束形成是最常用的一种非自适应波束形成方法。常规波束形成使用一组固定的加权和时间延迟对来自阵列中传感器的信号进行组合叠加，这个叠加过程主要是使用关于传感器在空间中的位置和感兴趣的信号方向的信息。常规波束形成的优点是它对假定的和实际的信号波前之间的信号失配具有较强的鲁棒性，并且只需要少量样本即可稳健地估计目标来波方向，易于实现，因此在单目标场景中有着较好的性能。目前基于常规波束形成的来波方向估计方法主要有：峰值搜索算法和抛物线插值算法。
[0005] 谱峰搜索法利用常规波束形成结果使角谱最大化，并返回M个最大值的位置作为M个目标的来波方向估计。谱峰搜索法对于白噪声中单个远场、窄带和电源的来波方向估计而言，是已知的最佳地最大似然估计方法。谱峰搜索法的来波方向估计精度取决于预引导波束扫描角的扫描间隔。为了获得更高的估计精度，需要更多的预引导波束扫描角。然而，更多的预引导波束扫描角将带来更多的计算负担。此外，受方位估计瑞利极限的限制，仅通过增加预引导波束扫描角的个数是无法无限提高来波方向估计精度的。
[0006] 抛物线插值算法是基于对主瓣角谱的抛物线拟合，易于实现。对于大多数目标的来波方向估计，估计精度较高。然而，对于某些来波方向，角谱的主瓣与抛物线模型不匹配，在这种情况下，来波方向估计精度下降。
发明内容
[0007] 技术问题：本发明提出了一种利用给定均匀线阵的理论波束图的模型求解算法。该算法是基于二分法的数值解，采用二分法数值求解信号来波方向相对偏差从而精确地估计远场窄带信号的来波方向。由于采用的是理论波束模式模型，因此不存在模型失配。与抛物线插值算法相比，该算法对信号来波方向相对偏差具有更强的鲁棒性，并且具有更好的估计性能。
[0008] 技术方案：本发明的一种基于指向性模型数值求解远场窄带信号来波方向估计方法，包括以下步骤：
[0009] 第一步：获取待处理的均匀线阵列M通道阵元远场窄带信号采样数据序列x1(n)，x2(n),…,xM(n),n＝0,1,…,N‑1：从所述阵元基元接收N个采样点的实时采集数据作为待处理的数据序列x1(n),x2(n),…,xM(n),n＝0,1,…,N‑1；或从存储器中提取从某一时刻开始接收到的N个信号采样点的数据作为待处理的数据序列x1(n),x2(n),…,xM(n),n＝0,1,…,N‑1，所述的M为阵元通道数，M为大于等于2的正整数，N为检测到的远场窄带信号脉宽长度所对应的采样点个数，取值为2的整数次幂且要求N>2。
[0010] 第二步：根据给定的阵列参数自动生成等余弦间隔的预引导波束扫描角θk和总的预引导波束个数K，自动生成过程如下，
[0011] 首先应根据给定的阵列参数计算总的预引导波束个数K，即
[0012]
[0013] 其中M为均匀线列阵的阵元个数， 代表向上取整运算。
[0014] 然后，依据总的预引导波束个数K自动计算生成预引导波束扫描角θk，即[0015]
[0016] 其中arccos[]代表取反余弦运算，k为预引导波束号。
[0017] 第三步：对已知频率f0的远场窄带信号进行频域波束形成，并求得该已知频率远场窄带信号的各预引导波束扫描角θk所对应的波束功率谱P(θk)，计算过程如下：
[0018] 首先，计算待处理的均匀线列阵的M通道的远场窄带信号采样数据序列x1(n)，x2(n),…,xM(n),n＝0,1,…,N‑1的离散傅里叶变换，其中第m个通道的远场窄带信号采样数据序列xm(n)的离散傅里叶变换Xm(l)的计算方法为：
[0019]
[0020] 其中m为阵元通道号，m＝1,2,…,M，l为离散频率索引，j表示虚数单位，即[0021] 然后，依据各预引导波束扫描角θk对待处理的均匀线性阵列M通道的远场窄带信号采样数据序列的离散傅里叶变换X1(l),X2(l),…,XM(l),l＝0,1,…,N‑1进行移相处理，其中对第m个通道的远场窄带信号采样数据序列的离散傅里叶变换Xm(l)做与各预引导波束扫描角θk相对应的移相处理结果 的计算方法为：
[0022]
[0023] 其中m为阵元通道号，m＝1,2,…,M，l为离散频率索引，k为预引导波束号，k＝1,2,…,K，d为均匀线列阵的间距，f0为已知的信号频率，c为声音的传播速度。
[0024] 最后，计算已知频率f0所对应的远场窄带信号各预引导波束扫描角θk所对应的波束功率谱P(θk)，计算过程如下：
[0025]
[0026] 其中||代表取模值运算，M为均匀线列阵通道数，l0为信号已知频率f0所对应的离散频率索引，其值为：
[0027] l0＝round(Nf0/fs) 式(6)
[0028] 其中round()代表四舍五入运算，fs为设定的信号采样频率。
[0029] 第四步：搜索波束功率谱P(θk)最大值所对应的预引导波束号kp，搜索过程如下：
[0030] 根据下式搜索波束功率谱P(θk)最大值所对应的预引导波束号kp：
[0031]
[0032] 其中 表示在1≤k≤K范围内搜索波束功率谱P(θk)的最大值所对应的预引导波束号。
[0033] 第五步：利用二分数值求解法求解信号来波方向相对偏差 过程如下：
[0034] 首先，计算波束功率谱P(θk)主瓣的峰值和左右相邻的预引导波束所对应的波束功率谱值中较大的一个的比值
[0035]
[0036] 然后，进行参数初始化：设置二分数值求解法的搜索上限a0＝0，搜索下限b0＝0.5，初始化迭代次数q＝1，设置二分精度ε，则二分总迭代次数NB可以表示为：
[0037]
[0038] 这里 表示向上取整运算，ε为大于0且小于0.25的数，log2()表示以2为底的代数运算；
[0039] 最后，利用二分数值求解法对来波方向相对偏差 进行迭代求解，迭代求解过程如下：
[0040] (1)令第q次迭代得到的来波方向相对偏差值为
[0041] (2)计算理论波束功率谱的主瓣峰值和左右相邻的预引导波束所对应的波束功率谱值中较大的一个的比值，计算如下
[0042]
[0043] 其中,d为均匀线列阵的间距，f0为已知的信号频率，c为声音的传播速度。
[0044] (3)重置二分求解法的上下限，即：
[0045] 如果 则令a0＝a0， 否则令 b0＝b0；
[0046] (4)判断是否满足q＝NB，如果满足，转入(5)；否则，令q＝q+1，并返回(1)；
[0047] (5)计算求得来波方向相对偏差 即
[0048] 第六步：利用 估计远场窄带信号的来波方向 计算过程如下：
[0049] 如果1＜kp＜K，则
[0050]
[0051] 如果kp＝1，则
[0052]
[0053] 如果kp＝K，则
[0054]
[0055] 其中 为预引导波束扫描角等余弦间隔值。
[0056] 进一步的，本发明方法的第三步中，第m个通道的远场窄带信号采样数据序列xm(n)的离散傅里叶变换Xm(l)是采用快速傅里叶变换得到的。
[0057] 进一步的，本发明方法的第五步中，通过利用二分数值求解法对来波方向相对偏差 进行迭代求解，二分精度ε设置为0.00001。
[0058] 本发明方法首先获取M通道阵元数据，然后根据给定的阵列参数生成等余弦间隔的预引导波束扫描角和总的预引导波束个数，接着对已知频率的远场窄带信号进行频域波束形成，并求得各预引导波束扫描角所对应的波束功率谱，搜索最大值对应的预引导波束号，然后利用二分数值求解法求解信号来波方向的相对偏差，最后利用该偏差估计远场窄带信号的来波方向。
[0059] 有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：
[0060] 1.本发明的估计方法采用理论波束模式模型，因此不存在模型失配的问题。
[0061] 2.本发明采用二分法数值求解信号来波方向相对偏差从而精确地估计远场窄带信号的来波方向，与抛物线插值算法相比，该算法对相对来波方向偏差具有更强的鲁棒性，并且具有更好的估计性能。
[0062] 3.本发明仅利用角谱峰值附近的三个样本进行来波方向估计，因而具有与抛物线插值算法相当的计算量。
[0063] 4.本发明可以推广应用于任意几何阵列的来波方向估计。
附图说明
[0064] 图1所示为本发明的流程图。
[0065] 图2所示为M通道阵元接收远场窄带信号示意图。
[0066] 图3所示为实施例1中仿真M通道阵元的预引导波束扫描角所对应的波束功率谱。
[0067] 图4所示为实施例2中仿真M通道阵元的预引导波束扫描角所对应的波束功率谱。
具体实施方式
[0068] 下面结合实施例和说明书附图对本发明作进一步的说明。
[0069] 如图1所示，一种基于指向性模型数值求解远场窄带信号来波方向估计方法，包括以下步骤：
[0070] 第一步：获取待处理的均匀线阵列M通道阵元远场窄带信号采样数据序列x1(n)，x2(n),…,xM(n),n＝0,1,…,N‑1：从所述阵元基元接收N个采样点的实时采集数据作为待处理的数据序列x1(n),x2(n),…,xM(n),n＝0,1,…,N‑1；或从存储器中提取从某一时刻开始接收到的N个信号采样点的数据作为待处理的数据序列x1(n),x2(n),…,xM(n),n＝0,1,…,N‑1，所述的M为阵元通道数，M为大于等于2的正整数，N为检测到的远场窄带信号脉宽长度所对应的采样点个数，取值为2的整数次幂且要求N>2。
[0071] 第二步：根据给定的阵列参数自动生成等余弦间隔的预引导波束扫描角θk和总的预引导波束个数K，自动生成过程如下，
[0072] 首先应根据给定的阵列参数计算总的预引导波束个数K，即
[0073]
[0074] 其中M为均匀线列阵的阵元个数， 代表向上取整运算；
[0075] 然后，依据总的预引导波束个数K自动计算生成预引导波束扫描角θk，即[0076]
[0077] 其中arccos[]代表取反余弦运算，k为预引导波束序号。
[0078] 第三步：对已知频率f0的远场窄带信号进行频域波束形成，并求得该已知频率远场窄带信号的各预引导波束扫描角θk所对应的波束功率谱P(θk)，计算过程如下：
[0079] 首先，计算待处理的均匀线列阵的M通道的远场窄带信号采样数据序列x1(n)，x2(n),…,xM(n),n＝0,1,…,N‑1的离散傅里叶变换，其中第m个通道的远场窄带信号采样数据序列xm(n)的离散傅里叶变换Xm(l)的计算方法为：
[0080]
[0081] 其中m为阵元通道号，m＝1,2,…,M，l为离散频率索引，j表示虚数单位，即[0082] 然后，依据各预引导波束扫描角θk对待处理的均匀线性阵列M通道的远场窄带信号采样数据序列的离散傅里叶变换X1(l),X2(l),…,XM(l),l＝0,1,…,N‑1进行移相处理，其中对第m个通道的远场窄带信号采样数据序列的离散傅里叶变换Xm(l)做与各预引导波束扫描角θk相对应的移相处理结果 的计算方法为：
[0083]
[0084] 其中m为阵元通道号，m＝1,2,…,M，l为离散频率索引，k为预引导波束号，k＝1,2,…,K，d为均匀线列阵的间距，f0为已知的信号频率，c为声音的传播速度；
[0085] 最后，计算已知频率f0所对应的远场窄带信号各预引导波束扫描角θk所对应的波束功率谱P(θk)，计算过程如下：
[0086]
[0087] 其中||代表取模值运算，M为均匀线列阵通道数，l0为信号已知频率f0所对应的离散频率索引，其值为：
[0088] l0＝round(Nf0/fs) 式(6)
[0089] 其中round()代表四舍五入运算，fs为设定的信号采样频率。
[0090] 第四步：搜索波束功率谱P(θk)最大值所对应的预引导波束号kp，搜索过程如下：
[0091] 根据下式搜索波束功率谱P(θk)最大值所对应的预引导波束号kp：
[0092]
[0093] 其中 表示在1≤k≤K范围内搜索波束功率谱P(θk)的最大值所对应的预引导波束号。
[0094] 第五步：利用二分数值求解法求解信号来波方向相对偏差 过程如下：
[0095] 首先，计算波束功率谱P(θk)主瓣的峰值和左右相邻的预引导波束所对应的波束功率谱值中较大的一个的比值
[0096]
[0097] 然后，进行参数初始化：设置二分数值求解法的搜索上限a0＝0，搜索下限b0＝0.5，初始化迭代次数q＝1，设置二分精度ε，则二分总迭代次数NB可以表示为：
[0098]
[0099] 这里 表示向上取整运算，ε为大于0且小于0.25的数，log2()表示以2为底的代数运算；
[0100] 最后，利用二分数值求解法对来波方向相对偏差 进行迭代求解，迭代求解过程如下：
[0101] (1)令第q次迭代得到的来波方向相对偏差值为
[0102] (2)计算理论波束功率谱的主瓣峰值和左右相邻的预引导波束所对应的波束功率谱值中较大的一个的比值，计算如下
[0103]
[0104] 其中d为均匀线列阵的间距，f0为已知的信号频率，c为声音的传播速度。
[0105] (3)重置二分求解法的上下限，即：
[0106] 如果 则令a0＝a0， 否则令 b0＝b0；
[0107] (4)判断是否满足q＝NB，如果满足，转入(5)；否则，令q＝q+1，并返回(1)；
[0108] (5)计算求得来波方向相对偏差 即
[0109] 第六步：利用 估计远场窄带信号的来波方向 计算过程如下：
[0110] 如果1＜kp＜K，则
[0111]
[0112] 如果kp＝1，则
[0113]
[0114] 如果kp＝K，则
[0115]
[0116] 其中 为预引导波束扫描角等余弦间隔值。
[0117] 本发明方法首先获取M通道阵元数据，然后根据给定的阵列参数生成等余弦间隔的预引导波束扫描角和总的预引导波束个数，接着对已知频率的远场窄带信号进行频域波束形成，并求得各预引导波束扫描角所对应的波束功率谱，搜索最大值对应的预引导波束号，然后利用二分数值求解法求解信号来波方向的相对偏差，最后利用该偏差估计远场窄带信号的来波方向。
[0118] 本发明的实施例中，阵元接收远场窄带信号示意图如图2所示，仿真阵元接收远场窄带信号模型为：
[0119]
[0120] 其中A为均匀线性阵列M通道阵元接收的远场窄带信号的幅度， 为阵元接收远场窄带信号的初始相位，N为检测到的远场窄带信号脉宽长度所对应的采样点个数，取值为2的整数次幂且要求N>2，M为总的阵元数目，M为大于等于2的正整数，f0为远场窄带信号频率，fs为采样频率，d为相邻阵元之间的距离，c为远场窄带信号在介质中的传播速度，θT为信号来波方向，即待估计值，wm(n)为第m个阵元接收的相互独立的高斯白噪声，均值为0，方差2 2 2
为σ且大小由信噪比SNR决定：SNR＝10log10[A/(2σ)]。
[0121] 实施例1
[0122] 仿真信号参数分别设置为：信号幅度A＝1，初始相位 远场窄带信号脉宽长度为1s，接收信号长度为1.032s,采样点个数N＝4096，总的阵元数目M＝21，远场窄带信号频率f0＝1250Hz，采样频率fs＝4000Hz，相邻阵元的间距d＝0.6m，远场窄带信号在介质中的传播速度c＝1500m/s(声音在水中的传播速度)，信号来波方向θT＝50.7035°，SNR＝‑6dB。
[0123] 首先根据给定的阵列参数自动生成总的预引导波束个数K＝25，由此我们可以得到预引导波束扫描角θk，即[0 23.5565 33.5573 41.4096 48.1897 54.3147 60.0000 65.3757 70.5288 75.5225 80.4059 85.2198 90.0000 94.7802 99.5941 104.4775
109.4712 114.6243 120.0000 125.6853 131.8103 138.5904 146.4427 156.4435
180.0000]°。
[0124] 然后对已知频率f0的远场窄带信号进行频域波束形成，并求得该已知频率远场窄带信号的各预引导波束扫描角θk所对应的波束功率谱P(θk)，如图3所示。
[0125] 接着通过搜索波束功率谱P(θk)最大值，得到所对应的预引导波束号为kp＝5，其对应幅值为 相邻两分量幅值分别为由此根据式(8)可以计算主瓣的峰值和左右相邻的预引导波束所对应
的波束功率谱值中较大的一个的比值为
[0126]
[0127] 然后设置二分数值求解法的搜索上限a0＝0，搜索下限b0＝0.5，初始化迭代次数q＝1，设置对分精度ε＝0.00001，则根据式(9)可以求出二分总迭代次数为[0128]
[0129] 接着利用二分数值求解法对来波方向相对偏差 进行迭代求解，得到[0130] 由于 根据式(11)得到来波方向 为
[0131]
[0132] 估计来波方向相对误差为
[0133] 实施例2
[0134] 仿真信号参数分别设置为：信号幅度A＝1，初始相位 远场窄带信号脉宽长度为1s，接收信号长度为1.032s,采样点个数N＝4096，总的阵元数目M＝15，远场窄带信号频率f0＝1250Hz，采样频率fs＝4000Hz，相邻阵元的间距d＝0.6m，远场窄带信号在介质中的传播速度c＝1500m/s(声音在水中的传播速度)，信号来波方向θT＝115.0576°，SNR＝‑6dB。
[0135] 首先根据给定的阵列参数自动生成总的预引导波束个数K＝18，由此我们可以得到预引导波束扫描角θk，即[0 28.0725 40.1192 49.6798 58.0343 65.6843 72.8954 79.8358 86.6277 93.3723 100.1642 107.1046 114.3157 121.9657 130.3202 139.8808
151.9275 180.0000]°。
[0136] 然后对已知频率f0的远场窄带信号进行频域波束形成，并求得该已知频率远场窄带信号的各预引导波束扫描角θk所对应的波束功率谱P(θk)，结果如图4所示。
[0137] 接着通过搜索波束功率谱P(θk)最大值，得到所对应的预引导波束号为kp＝13，其对应幅值为 相邻两分量幅值分别为由此根据式(8)可以计算主瓣的峰值和左右相邻的预引导波束所对应
的波束功率谱值中较大的一个的比值为
[0138]
[0139] 然后设置二分数值求解法的搜索上限a0＝0，搜索下限b0＝0.5，初始化迭代次数q＝1，设置对分精度ε＝0.00001，则根据式(9)可以求出二分总迭代次数为[0140]
[0141] 接着利用二分数值求解法对来波方向相对偏差 进行迭代求解，得到由于 根据式(11)则来波方向 为
[0142]
[0143] 估计来波方向相对误差为
[0144] 从实施例1和例2的结果可以看出，本发明估计方法可以获得良好的估计精度，而且计算简单，计算量小，适用于高精度快速估计多阵元远场窄带信号来波方向的场合。
[0145] 上述实施例仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和等同替换，这些对本发明权利要求进行改进和等同替换后的技术方案，均落入本发明的保护范围。
