基于交替方向惩罚的频率一致宽带波束形成方法及系统

1.本发明属于波束形成技术领域，尤其涉及基于交替方向惩罚的频率一致宽带波束形成方法及系统。


背景技术：

2.宽带波束形成是声学领域的重要研究内容，包括频率不变波束形成(fib)。频率不变波束形成技术可在接受目标信号范围内形成一定宽度的波束，以抑制信号失真，在实际中可用于水声中波场的采样。
3.为满足实际需要，超指向波束形成技术被开发出来。当线性阵列的阵元间距趋于0时，其方向性因子(df)在端射方向达到最大值，此时紧凑的孔径易于嵌入到其他器件中，同时提供高的方向增益。然而，超指向波束形成器在低频时仍然对白噪声和阵列不匹配很敏感。为解决该问题，许多设计波束形成器的优化问题被开发出来，优化问题常用的求解方法之一是内点法，但内点法往往会耗费大量的计算时间。
4.近年来，交替方向乘子法(admm)在求解优化问题方面也取得了一定的成就，可以实现高效求解非凸非线性优化问题，但存在不收敛问题，且运算过程也有待简化。


技术实现要素：

5.本发明的目的是提供基于交替方向惩罚的频率一致宽带波束形成方法及系统。
6.本发明的技术方案如下：
7.基于交替方向惩罚的频率一致宽带波束形成方法，包括：
8.步骤一，以hh(ωn)γ(ωn)h(ωn)最小为目标函数，结合无失真约束、主瓣约束、副瓣约束以及白噪声增益约束构建优化问题；其中，h(ωn)表示角频率ωn对应的权值向量，γ(ωn)表示伪相干矩阵；
9.步骤二，利用交替方向惩罚法求解优化问题获得最优权值向量，该步骤进一步包括：
10.2.1引入似实值助变量和对偶辅助向量重新表达优化问题；
11.2.2引入各线性约束的缩放对偶变量构造优化问题的增广拉格朗日函数；
12.2.3对权值向量、似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量进行交叉迭代求解；
13.2.4重复子步骤2.3直至达到收敛条件，此时的权值向量即最优的权值向量；
14.步骤三，将最优的权值向量与导向矢量相乘结果作为频率一致宽带波束形成。
15.进一步的，步骤一所构建的优化问题如下：
[0016][0017]
其中：
[0018]
h(ωn)表示角频率ωn对应的权值向量，γ(ωn)表示伪相干矩阵；
[0019]hh
(ωn)d(ωn,cosθd)＝1为无失真约束，d(ωn,cosθd)为角频率ωn对应的导向矢量，θd为期望信号的入射方向角；
[0020]hh
(ωn)h(ωn)≤γ(ωn)为白噪声增益约束，γ(ωn)为角频率ωn处白噪声增益的上限；
[0021]
分别为主瓣约束和旁瓣约束；
[0022]
分别为主瓣区域和副瓣区域的波束图向量，dm(ωn)、ds(ωn)分别为角频率ωn下覆盖主瓣区域和副瓣区域的导向矢量所组成的导向矩阵，ε1(ωn)和ε2(ωn)分别为主瓣和副瓣对应的正误差参数；
[0023]
以及子步骤2.1中重新表达的优化问题如下：
[0024][0025]
其中：s1(ωn)和s2(ωn)为引入的似实值助变量，z(ωn)为引入的对偶辅助向量。
[0026]
进一步的，子步骤2.2中构造的增广拉格朗日函数如下：
[0027][0028]
其中：ρ1、ρ2、ρ3、ρ4均为自定义的惩罚因子，也即迭代步长，ρ1＞0，ρ2＞0，ρ3＞0，ρ4》0；μ(ωn)为拉格朗日乘子；
[0029]
μ(ωn)为线性约束hh(ωn)d(ωn,cosθd)＝1对应的缩放对偶变量；
[0030]
v(ωn)为线性约束z(ωn)＝h(ωn)对应的缩放对偶变量；
[0031]
w(ωn)为线性约束对应的缩放对偶变量；
[0032]
r(ωn)为线性约束对应的缩放对偶变量。
[0033]
基于上述构造的增广拉格朗日函数，子步骤2.3进一步包括：
[0034]
2.3a初始化权值向量、似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量；
[0035]
2.3b基于似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，通过求解无约束优化问题，对权值向量进行更新；
[0036]
2.3c基于权值向量、似实值助变量、惩罚因子、缩放对偶变量的当前值，仅在与对偶辅助向量相关的约束条件下求解优化问题，对偶辅助向量进行更新；
[0037]
2.3d基于权值向量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，仅在与似实值助变量相关的约束条件下求解优化问题，对似实值助变量进行更新；
[0038]
2.3e基于权值向量、似实值助变量、对偶辅助向量、惩罚因子的当前值，更新缩放对偶变量；
[0039]
2.4f基于权值向量、、似实值助变量、对偶辅助向量、缩放对偶变量的当前值，更新惩罚因子。
[0040]
进一步的，子步骤2.3b中，更新后的权值向量
[0041]
其中，h
(k+1)
(ωn)为更新后的权值向量；ξ1和ξ2如下：
[0042][0043]
[0044]
进一步的，子步骤2.3c中，所述与对偶辅助向量相关约束条件为其中，z(ωn)为对偶辅助向量，γ(ωn)为角频率ωn处白噪声增益的上限；
[0045]
更新后的对偶辅助变量z
k+1
(ωn)如下：
[0046][0047]
其中：
[0048]
进一步的，子步骤2.3d中，先基于权值向量、第二似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，仅在与第一似实值助变量相关的约束条件下求解优化问题，对第一似实值助变量进行更新；
[0049]
再基于权值向量、第一似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，仅在与第二似实值助变量相关的约束条件下求解优化问题，对第二似实值助变量进行更新。
[0050]
进一步的，与第一似实值助变量相关的约束条件为更新后的第一似实值助变量其中：
[0051][0052]
与第二似实值助变量相关的约束条件为更新后的第二似实值助变量其中：
[0053][0054]
进一步的，子步骤2.4中收敛条件为：
[0055]
迭代次数达到预设的最大迭代次数或同时满足以及其中，δ1和δ2均为公差参数阈值，
均为各线性约束条带的当前迭代原始残差，为当前迭代的对偶残差。
[0056]
本发明提供的基于交替方向惩罚的频率一致宽带波束形成系统，包括：
[0057]
第一模块，用来以hh(ωn)γ(ωn)h(ωn)最小为目标函数，结合无失真约束、主瓣约束、副瓣约束以及白噪声增益约束构建优化问题；其中，h(ωn)表示角频率ωn对应的权值向量，γ(ωn)表示伪相干矩阵；
[0058]
第二模块，用来利用交替方向惩罚法求解优化问题获得最优权值向量；
[0059]
所述第二模块进一步包括子模块：
[0060]
第一子模块，用来引入似实值助变量和对偶辅助向量重新表达优化问题；
[0061]
第二子模块，用来引入各线性约束的缩放对偶变量构造优化问题的增广拉格朗日函数；
[0062]
第三子模块，用来对权值向量、似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量进行交叉迭代求解；
[0063]
第四子模块，用来判断重复第三子模块直至达到收敛条件，此时的权值向量即最优的权值向量；
[0064]
第三模块，用来将最优的权值向量与导向矢量相乘结果作为频率一致宽带波束形成。
[0065]
和现有技术相比，本发明具有如下优点和有益效果：
[0066]
本发明设计了频率一致波束形成器，可以确保近似频率一致以及频率相关的白噪声增益值。目前优化问题求解方法普遍存在效率低且非收敛的问题，加之在不同频率上搜索惩罚参数也存在一定难度，本发明对交替方向法进行改进，提出了一种交替方向惩罚法来求解优化问题。该方法利用残差和惩罚因子间的关系来减少迭代次数，可以更高效的求解优化问题，尤其是在求解具大量阵元优化问题时，求解效率更是得到极大提高。
附图说明
[0067]
图1为本发明方法流程示意图；
[0068]
图2为本发明方法在仿真实验中不同惩罚系数下得到的目标函数df与迭代次数的关系；
[0069]
图3为本发明方法在不同惩罚系数下得到的对偶残差eh(ωn)与迭代次数的关系；
[0070]
图4为本发明方法在不同惩罚系数下得到的原始残差e
r1
(ωn)与迭代次数的关系；
[0071]
图5为本发明方法在不同惩罚系数下得到的原始残差e
r3
(ωn)与迭代次数的关系；
[0072]
图6为本发明方法在不同惩罚系数下得到的原始残差e
r4
(ωn)与迭代次数的关系；
[0073]
图7为本发明方法在仿真条件m＝8，f＝200hz，下的波束图；
[0074]
图8为本发明方法在m＝8，δ＝1cm仿真条件下的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0075]
图9为非相干法在m＝8，δ＝1cm仿真条件下的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0076]
图10为mn法在m＝8，δ＝1cm仿真条件下的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0077]
图11为ice法在m＝8，δ＝1cm仿真条件下的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0078]
图12为本发明方法、非相干法、mn法、ice法的df曲线；
[0079]
图13为本发明方法、非相干法、mn法、ice法的wng曲线；
[0080]
图14为本发明方法在m＝8，δ＝1cm仿真条件下α为0.2时的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0081]
图15为本发明方法在m＝8，δ＝1cm仿真条件下α为0.45时的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0082]
图16为本发明方法在m＝8，δ＝1cm仿真条件下α为0.6时的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0083]
图17为本发明方法在m＝8，δ＝1cm仿真条件下α为0.8时的宽带波束图和在f＝2500hz下的极坐标波束图，其中，图(a)为宽带波束度，图(b)为极坐标波束图；
[0084]
图18为α取0.2、0.45、0.6、0.8时本发明方法的df；
[0085]
图19为α取0.2、0.45、0.6、0.8时本发明方法的wng曲线。
具体实施方式
[0086]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图对本发明具体实施方式进行详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明。此外，下面所描述的具体实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0087]
本文中所涉及字符参数含义如下：
[0088]
小写的粗体符号表示向量，大写的粗体符号表示矩阵，(大写)。为虚数单位；(
·
)
t
、(
·
)h、(
·
)
*
分别表示转置算子、共轭转置算子、复共轭算子；vec(
·
)表示向量化运算符；im表示大小m
×
m的单位矩阵；||
·
||2表示二范数；|
·
|表示模值；o(
·
)表示算法复杂度的朗道符号(landau符号)。
[0089]
为便于理解，下面将对本发明所基于的技术原理及所采用技术方案分别详细说明。
[0090]
一、信号模型构建
[0091]
考虑基于均匀间隔的线性阵列的宽带波束形成结构，远场源信号(平面波)以声速c＝340m/s入射到由全方位传感器组成的阵列上。两个连续传感器间的距离等于δ，远场源
信号到阵列的入射方位角采用θ表示，导向矢量d(ω,cosθ)由式(1)给出：
[0092][0093]
式(1)中：
[0094]
ω表示角频率，ω＝2πf，f表示时间频率，f＞0，则波长λ＝c/f；
[0095]
τ0＝δ/c是入射方位角θ＝0
°
处两个连续传感器间的时间差；
[0096]
e表示自然对数，m表示阵列中传感器数量，即阵列中阵元数量。
[0097]
阵列的波束图的定义见式(2)，波束图用来描述阵列对来自θ方向信号的响应。
[0098][0099]
式(2)中；
[0100]
h(ω)表示角频率ω对应的复数权值向量，h1(ω)、h2(ω)、
…hm
(ω)为第m个传感器的权值，m取1、2、
…
m。
[0101]
假设期望信号从入射方位角θd传播，则观测信号向量y(ω)为：
[0102][0103]
式(3)中：ym(ω)表示第m个传感器的观测信号；d(ω,cosθd)表示期望信号源所在方向的导向矢量；x(ω)表示期望信号；vn(ω)表示加性噪声信号矢量。
[0104]
通常波束形成是通过对每个传感器信号应用时间滤波器，并对滤波后的信号求和来完成。在频域中，相当于向每个传感器的输出添加一个复数权重，并对整个阵元直径求和，即：
[0105][0106]
式(4)中，z(ω)表示波束形成器输出信号；x
fd
(ω)表示经过滤波器的所需信号，v
rn
(ω)表示残余噪声。
[0107]
本发明的第一个目标是设计超定向波束形成器，其波束图频率一致，同时尽可能接近期望的功率波束图；第二个目标是找到一种快速收敛方法来求解。下面将分别对优化问题构建和求解方法进行介绍。
[0108]
二、优化问题构建
[0109]
给定参考波束图即期望的功率波束图，使用均方误差约束来获得所需的频率一致波束图，可写成：
[0110]
[0111]
式(5)中，ω和θ分别表示感兴趣的角频率范围和入射角度范围；ε
ls
为正公差参数，表示合成波束图和参考波束图间的误差。
[0112]
虽然式(5)中合成波束图考虑到所有的入射角度和角频率，但不能考虑不同角频率下获得频率一致波束的能力差异，主瓣区域的频率不变性比副瓣区域的频率不变性更重要，这意味着必须优先考虑以确保主瓣的频率不变性。因此，统一离散角频率和入射角度，引入i个方向和j个角频率单元，将主瓣和副瓣分别划分为l和s个角度。式(5)重写为：
[0113][0114]
式(6)中，上标m表示主瓣，上标s表示副瓣，θm表示主瓣区域，θs表示副瓣区域，c是表示微分dωdθ的正常数；
[0115]
表示主瓣区域波束图向量，θ
l
表示第l个主瓣角度，l＝1,2,...,l，表示第l个主瓣角度的波束图；
[0116]dm
(ωn)表示大小m
×
l的导向矩阵，即角频率ωn下覆盖主瓣区域的导向矢量所组成的导向矩阵，其表示见式(7)：
[0117]dm
(ωn)＝[d(ωn,cosθ1),d(ωn,cosθ2),...,d(ωn,cosθ
l
)]
ꢀꢀꢀꢀꢀꢀꢀ
(7)
[0118]
式(7)中，d(ωn,cosθ
l
)表示角频率ωn的信号源在入射方位角θ
l
的导向矢量。则为副瓣区域波束图向量，ds(ωn)为角频率ωn下覆盖副瓣区域的导向矢量所组成的导向矩阵。
[0119]
关于约束条件，无失真约束是必须的，无失真约束条件如下：
[0120][0121]
式(8)中，ωn表示第n个频率，h(ωn)则为角频率ωn对应的复数权值向量，d(ωn,cosθd)为角频率ωn对应的导向矢量。
[0122]
无失真约束意味：任何到达的信号将沿着导向矢量d(ωj,cosθd)在不失真情况下通过波束形成器。
[0123]
将式(6)的第一个求和项分解为多个联合约束，得到式(9)所示的主瓣约束，对有：
[0124][0125]
类似地，副瓣约束见式(10)，对有：
[0126][0127]
式(9)～(10)中，ε1(ωn)和ε2(ωn)分别为主瓣和副瓣对应的正误差参数。
[0128]
由于频率不变性决定了各种性能指标，且主瓣比副瓣的影响更大，根据式(5)和(6)，建立以下关系：
[0129][0130]
评估阵列对某些缺陷(如传感器噪声)灵敏度的最方便方法是噪声增益(wng)，白噪声增益的表达式为：
[0131][0132]
白噪声增益也是一种在存在不相关白噪声的情况下阵列增益的度量的指标。假设信号不失真通过滤波器，考虑白噪声增益可以等效于白噪声增益的分母，即可将白噪声增益的值限定在某个阈值之上，对有：
[0133][0134]
它保证了阵列的鲁棒性，γ(ωn)表示在角频率ωn处白噪声增益的最大值上限。
[0135]
为更好地抑制声场中的噪声，需要优化组合多个麦克风，并考虑方向增益。使用的一个常用量是方向性因子，对球形漫反射噪声，方向性因子定义为：
[0136][0137]
式(14)中：γ(ωn)为大小m
×
m的伪相干矩阵，对应于球各向同性噪声；γ(ωn)中第(i，l)个元素
[0138]
本发明拟找到一种波束形成器，其有着理想的df且波束图是是fi的，同时可灵活调节波束图的频率一致性，结合约束条件构建如下优化问题：
[0139]
[0140]
式(15)中：表示目标函数，以hh(ωn)γ(ωn)h(ωn)最小为目标；其他为约束条件。
[0141]
对于上述优化问题，一般处理方法采用内点法，但优化问题中所有约束相互耦合，因为都为向量h(ωn)的函数，采用内点法的计算复杂度为o(m
3.5
(l+s))，这会耗费大量的计算时间，因此需要采用一种快速收敛方法来解决该问题。
[0142]
三、基于交替方向惩罚的频率一致波束形成方法
[0143]
引入似实值助变量和引入对偶辅助向量z(ωn)＝h(ωn)，重构式(15)的优化问题，对有：
[0144][0145]
利用目标函数和线性约束以缩放形式来构造优化问题(16)的增广拉格朗日函数：
[0146][0147]
式(17)中：ρ1、ρ2、ρ3、ρ4均为自定义的惩罚因子，也即迭代步长，ρ1＞0，ρ2＞0，ρ3＞0，ρ4》0；μ(ωn)为拉格朗日乘子；
[0148]
μ(ωn)为线性约束hh(ωn)d(ωn,cosθd)＝1对应的缩放对偶变量；
[0149]
v(ωn)为线性约束z(ωn)＝h(ωn)对应的缩放对偶变量；
[0150]
w(ωn)为线性约束对应的缩放对偶变量；
[0151]
r(ωn)为线性约束对应的缩放对偶变量。
[0152]
更新在每次迭代中进行，而不是在双重上升法中直到内部循环收敛的多次更新。更重要的是，本发明交替方向惩罚在admm框架上完成，不仅限于以上常见的优点。交替方向惩罚法基于原始残差更新惩罚因子和缩放对偶变量，进一步提高收敛能力，这种省时的实现策略可以确保交替方向惩罚法有更高概率找到比admm和内点法更优的可行点。
[0153]
本发明利用交替方向惩罚法(admp)求解优化问题的过程如下：
[0154]
以d(ωn,cosθd)、dm(ωn)、ds(ωn)、γ(ωn)、ε1(ωn)、ε2(ωn)为输入，h(ωn)为输出，根据h(ωn)可形成波束图。求解过程如下：
[0155]
s100：初始化：令迭代次数k为0，最大迭代次数设为k
max
；初始化参数其中，a＝1,2,3,4；参数一般初始化为0或经验值。
[0156]
s200：通过求解式(18)所示的无约束优化问题，来更新h(ωn)：
[0157][0158]
zk(ωn),μk(ωn),vk(ωn),wk(ωn),rk(ωn)为当前参量值，忽略与h(ωn)无关的项，则优化问题简化为：
[0159][0160]
每次迭代中都有h(ωn)闭合解，h(ωn)更新公式如下：
[0161][0162]
其中：
[0163][0164][0165]
步骤300：由于最小化问题中不能解决含有四次项的式子，因此引进对偶辅助变量z(ωn)，使用h
k+1
(ωn),μk(ωn),vk(ωn),wk(ωn),rk(ωn)解决以下优化问题，来更新z(ωn)：
[0166][0167]
将式(17)代入式(23)并移除无关项，式(23)改写为：
[0168][0169]
采用式(24)更新对偶辅助变量z(ωn)：
[0170][0171]
其中：
[0172]
步骤400：使用h
k+1
(ωn),z
k+1
(ωn),s
2k
(ωn),μk(ωn),vk(ωn),wk(ωn),rk(ωn)解决式(25)所示优化问题并更新s1(ωn)：
[0173][0174]
将式(17)插入(25)并移除无关项，则优化问题被改写为：
[0175][0176]
s1(ωn)的更新公式如下：
[0177][0178]
其中，
[0179]
同理可得s2(ωn)的更新公式如下：
[0180][0181]
其中，
[0182]
步骤500：采用式(29)～(30)更新缩放对偶变量μ(ωn)、v(ωn)、w(ωn)、r(ωn)：
[0183][0184][0185][0186][0187]
式(29)～(32)中涉及变量定义如下：
[0188][0189]
其中，η2、η3、η4为足够大的正数。
[0190]
步骤600：采用式(34)更新惩罚因子：
[0191][0192]
其中，其中，
[0193]
参数δ
a1
、δ
a2
的值均应接近1，该值与收敛速度及相应的残值有关，而惩罚因子的更新与原始残差有关。具体而言，如果不随迭代次数减少，则使用更大的惩罚参数
使接近零，从而以高概率找到行点；否则，保持不变。
[0194]
s700：重复迭代步骤s200～s600，直至k达到k
max
；或直至同时满足以及其中，δ1和δ2均为公差参数阈值，例如可设δ1＝10-4
，δ2＝10-3
，具体可根据实际需求调整δ1和δ2值。
[0195]
本步骤中对有：
[0196][0197]
式(24)中：均为第(k+1)次迭代时的原始残差，表示第(k+1)次迭代时的对偶残差。
[0198]
上标k和k+1分别表示该变量为迭代更新前后的值。
[0199]
分析本发明所采用算法的计算复杂度，每次迭代中对h(ωn)、ξ1、ξ2的更新，计算复杂度分别为o(m2l)、o(m3)和o(ml)，步骤1的计算复杂度共计o(ml+m2l+m3)；步骤2的计算仅涉及简单的向量运算，其复杂度为o(m)；步骤3是传统admm算法的一部分，与步骤1～2相同，每次迭代的复杂度为o(ml+ms)；步骤4与步骤5相似，每次迭代的复杂度为o(2m+ml+ms)。因此，本发明算法每次迭代的计算复杂度共计o(m3+m2l+3ml+4ms+5m)。
[0200]
四、仿真实验
[0201]
本仿真实验在均匀线阵比较本发明方法、非相干法、最小范数法(简记为mn法)、不相关拓展法(简记为ice法)的性能，非相干法是以hh(ωn)h(ωn)最小为目标函数，以为约束条件来构建优化问题，并采用内点法求解。mn法是基于n阶差分麦克风的最小范数设计，其将(差分麦克风)dma的阶数与传感器数量不关联化，并在固定阶数下增加阵元数。ice法与非相干法类似，唯一区别在于ice法中误差参数ε1(ωn)、ε2(ωn)和阵元间距δ相对增加。
[0202]
将这四种方法用于在f
low
＝200hz和f
high
＝4480hz间的频率范围内，设计频率一致的宽带波束图。假设用于相应时域接收信号的窗口分析的持续时间t＝25毫秒，频率分辨率
f＝1/t＝40hz。因此，频点数量最后，考虑为均匀线阵，使用8个全向麦克风，阵元间距
[0203]
通过推导三阶锐心形差分波束形成器，将理论差分波束图的零点设置为三个，并基于最大化df设计一个三阶锐心形波束图，阶数n＝3，获得其理论波束图
[0204][0205]
式(25)中，a
n,n
是确定n阶差分麦克风不同方向图的实归一化系数值，n＝0,1,...,n；an是长度向量，an＝[a
n,0
,a
n,1
,...a
n,n
]
t
。
[0206]
该三阶锐心形波束图可表示为b
3,hd
(θ)，b
3,hd
(θ)为期望的波束图：
[0207]b3,hd
(θ)＝-0.14-0.57cosθ+0.57cos2θ+1.15cos3θ
ꢀꢀꢀꢀꢀꢀꢀꢀ
(26)
[0208]
建立矩阵dm(ωn)和ds(ωn)，令离散角δθ＝1
°
，ρ
l
＝θ
l
＝60
°
，ρs＝θs＝120
°
，即主瓣区域是[-60
°
,-60
°
]，旁瓣区域是[-180
°
,-61
°
]∪[61
°
,180
°
]∪[61
°
,180
°
]。
[0209]
对于本发明方法，设置参数为：ε1(ωn)＝0.1，ε2(ωn)＝0.2，α(ωn)＝0.4，δr1＝δr2＝δr3＝δr4＝10-6
，δ1＝10-6
，δ2＝10-3
，δ
11
＝0.9985，δ
21
＝0.9995，δ
31
＝0.9995，δ
41
＝0.9995，δ
12
＝1.025，δ
22
＝1.005，δ
32
＝1.005，δ
42
＝1.005。对于非相干法，δra、δ1、δ2取值同本发明，但δ
11
＝0.95，δ
12
＝1.05。
[0210]
求解时目标函数在一定范围收敛。为测试迭代次数对收敛性的影响，令采用本发明方法进行求解，仿真条件为：f＝200hz，m＝8，δ＝1cm。不同收敛范围下的收敛结果见图2，原始残差e
r1
(ωn)与迭代次数的关系见图3，原始残差e
r3
(ωn)与迭代次数的关系见图4，原始残差e
r4
(ωn)与迭代次数的关系见图5，对偶残差eh(ωn)与迭代次数的关系见图6，仿真波束图见图7。
[0211]
图3～6给出了残差随迭代次数的变化规律，从图中可以看出，残差岁迭代次数而降低，且四个残差最终都趋于0，得到了收敛。更准确来说，在一定迭代次数后，原始残差变为0，而对偶残差接近0。通过比较对收敛速度的影响，从图2～6可以看出，少量改变对收敛速度有显著影响，而值的改变对收敛速度只有少量影响，在迭代过程中比对目标函数收敛速度影响更大，于是令但
[0212]
图2给出了七种下的目标值df，从图2可以看出，随着迭代次数增加，均可收敛到11.7db～11.9db。因为本发明方法所构建的优化问题是在整个宽带上，处理问题的效率很重要，因此仅考虑目标函数的值是不够的，更希望能在目标函数值和收敛速度上折中。目标函数值与惩罚因子没有直观的关系,这也是由于惩罚系数也随着迭代次数变化，这点不同于传统的admm。见图7，图中虚线图形代表本发明方法所求波束图，实线图像代表期望波束图；从图7可以观察到，在特定的惩罚参数
下设计的波束图，低频时非常接近期望波束图，这意味着本发明可以很好地求解所构建的优化问题。
[0213]
见图8～11，分别为上述四种方法在m＝8，δ＝1cm仿真条件下的宽带波束图和在f＝2500hz下的极坐标波束图。其中图8(b)、图9(b)、图10(b)、图11(b)展示了特定频率f＝2500hz下极坐标波束图(见虚线图形)与期望波束图(见实线图形)(即，三阶锐心形波束图)的比较。本发明方法、非相干法、mn法的阵元线阵布局相同，ice法的误差参数和阵元间距在非相关方法的基础上发生变化。从图中可以看出，mn法所得波束图与期望波束图最相似，而ice法所得波束图具有更高旁瓣，本发明利用ice法可找到一个有效的白噪声输出功率上限，但会牺牲波束图频率一致，以便能在特定范围内获得效果较好的频率一致波束图。对于有着更好白噪声增益性能的非相关法，其所得波束图的波瓣与期望波束图波瓣几乎相似，但代价是旁瓣和主瓣间的误差相对较大。
[0214]
见图8(a)、图9(a)、图10(a)、图11(a)，本仿真实验中，期望波束图在200hz和2500hz间的频率范围内具有三个零点，具体而言，在200hz和1000hz间，波束图与期望波束图案几乎相似，第二个和第三个零点的角度逐渐增大，第三个零点从154
°
增加到180
°
直至消失。总体而言，随着频率增加，零点的数量大约从f＝2500hz开始减少到两个，从深到浅变化。对于本发明方法和非相干法，所得宽带波束图几乎频率一致，且在主瓣区域误差ε1(ωn)较小，在旁瓣区域误差ε2(ωn)较大。mn法实现了最完美的频率一致波束图，波束方向图在所有频率的主瓣和旁瓣范围内都是一致的；对于ice法，对于大部分频率都有两个零点，第二个零点在3000hz和4480hz不明显，频率一致性较差。
[0215]
图12～13为四种方法在m＝8，δ＝1cm仿真条件下的df和wng曲线，通过与期望波束图的df和wng曲线比较，本发明方法和不相干法具有与期望波束图最相似的结果。
[0216]
本发明方法具有比不相干法更高效的运算速度，同时在64gb内存的amd ryzen 9 3900x 12core上运行本发明方法和不相干法，在不包括参数搜索过程的情况下，本发明方法大约需要15秒运行时间，而不相干法则运行了将近90秒，几乎是本发明方法的六倍。虽然mn法可获得频率一致的相对最优df，这意味着它的波束图频率一致性完美，但会遭受低频白噪声的放大，对于较少的线性阵列阵元，df几乎恒定，最高可达8000hz；但随着m的增加，df从大约f＝2000hz开始不断减小；但在高频下，本发明方法可获得更高的wng。
[0217]
图14～图17展示了在折中参数α为0.2、0.45、0.6、0.8时本发明方法所获得的宽带波束图和f＝2500hz下的极坐标波束图，在极坐标波束图中，实现图形为期望波束图。图18～图19展示了α为0.2、0.45、0.6、0.8时本发明方法的df和wng曲线。从图14～19可以看出，本发明fi超指向性波束形成器在wng和df间取得了很好的折中，特别是α为0.45时。随着α增加，以牺牲wng为代价，df得到改善，并且波束图零点更清晰。
[0218]
本发明所提出的频率不变波束形成器，所形成波束图主瓣区域几乎频率不变，同时还可以实现所需的wng值，并具有一定的灵活性以最大化df。与不相关法相比，本发明在显著降低计算复杂度情况下还实现了良好的性能，并且优化方法可用于所有方向和感兴趣频率。
[0219]
上述实施例所述是用以具体说明本发明，文中虽通过特定的术语进行说明，但不能以此限定本发明的保护范围，熟悉此技术领域的人士可在了解本发明的精神与原则后对其进行变更或修改而达到等效目的，而此等效变更和修改，皆应涵盖于权利要求范围所界
定范畴内。

技术特征：
1.基于交替方向惩罚的频率一致宽带波束形成方法，其特征是，包括：步骤一，以h
h
(ω
n
)γ(ω
n
)h(ω
n
)最小为目标函数，结合无失真约束、主瓣约束、副瓣约束以及白噪声增益约束构建优化问题；其中，h(ω
n
)表示角频率ω
n
对应的权值向量，γ(ω
n
)表示伪相干矩阵；步骤二，利用交替方向惩罚法求解优化问题获得最优权值向量，该步骤进一步包括：2.1引入似实值助变量和对偶辅助向量重新表达优化问题；2.2引入各线性约束的缩放对偶变量构造优化问题的增广拉格朗日函数；2.3对权值向量、似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量进行交叉迭代求解；2.4重复子步骤2.3直至达到收敛条件，此时的权值向量即最优的权值向量；步骤三，将最优的权值向量与导向矢量相乘结果作为频率一致宽带波束形成。2.如权利要求1所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：步骤一所构建的优化问题如下：其中：h(ω
n
)表示角频率ω
n
对应的权值向量，γ(ω
n
)表示伪相干矩阵；h
h
(ω
n
)d(ω
n
，cosθ
d
)＝1为无失真约束，d(ω
n
，cosθ
d
)为角频率ω
n
对应的导向矢量，θ
d
为期望信号的入射方向角；h
h
(ω
n
)h(ω
n
)≤γ(ω
n
)为白噪声增益约束，γ(ω
n
)为角频率ω
n
处白噪声增益的上限；分别为主瓣约束和旁瓣约束；分别为主瓣区域和副瓣区域的波束图向量，d
m
(ω
n
)、d
s
(ω
n
)分别为角频率ω
n
下覆盖主瓣区域和副瓣区域的导向矢量所组成的导向矩阵，ε1(ω
n
)和ε2(ω
n
)分别为主瓣和副瓣对应的正误差参数；以及子步骤2.1中重新表达的优化问题如下：
其中：s1(ω
n
)和s2(ω
n
)为引入的似实值助变量，z(ω
n
)为引入的对偶辅助向量。3.如权利要求1所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：子步骤2.2中构造的增广拉格朗日函数如下：其中：ρ1、ρ2、ρ3、ρ4均为自定义的惩罚因子，也即迭代步长，ρ1＞0，ρ2＞0，ρ3＞0，ρ4＞0；λ(ω
n
)为拉格朗日乘子；μ(ω
n
)为线性约束h
h
(ω
n
)d(ω
n
，cosθ
d
)＝1对应的缩放对偶变量；v(ω
n
)为线性约束z(ω
n
)＝h(ω
n
)对应的缩放对偶变量；w(ω
n
)为线性约束对应的缩放对偶变量；r(ω
n
)为线性约束对应的缩放对偶变量。基于上述构造的增广拉格朗日函数，子步骤2.3进一步包括：2.3a初始化权值向量、似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量；2.3b基于似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，通过求解无约束优化问题，对权值向量进行更新；2.3c基于权值向量、似实值助变量、惩罚因子、缩放对偶变量的当前值，仅在与对偶辅助向量相关的约束条件下求解优化问题，对偶辅助向量进行更新；2.3d基于权值向量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，仅在与似实值助变量相关的约束条件下求解优化问题，对似实值助变量进行更新；
2.3e基于权值向量、似实值助变量、对偶辅助向量、惩罚因子的当前值，更新缩放对偶变量；2.3f基于权值向量、似实值助变量、对偶辅助向量、缩放对偶变量的当前值，更新惩罚因子。4.如权利要求1所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：子步骤2.3b中，更新后的权值向量其中，h
(k+1)
(ω
n
)为更新后的权值向量；ξ1和ξ2如下：如下：5.如权利要求3所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：子步骤2.3c中，所述与对偶辅助向量相关约束条件为其中，z(ω
n
)为对偶辅助向量，γ(ω
n
)为角频率ω
n
处白噪声增益的上限；更新后的对偶辅助变量z
k+1
(ω
n
)如下：其中：6.如权利要求3所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：子步骤2.3d中，先基于权值向量、第二似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，仅在与第一似实值助变量相关的约束条件下求解优化问题，对第一似实值助变量进行更新；再基于权值向量、第一似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量的当前值，仅在与第二似实值助变量相关的约束条件下求解优化问题，对第二似实值助变量进行更新。7.如权利要求6所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：所述与第一似实值助变量相关的约束条件为更新后的第一似实值助变量其中：
所述与第二似实值助变量相关的约束条件为更新后的第二似实值助变量其中：8.如权利要求3所述的基于交替方向惩罚的频率一致宽带波束形成方法，其特征是：子步骤2.4中收敛条件为：迭代次数达到预设的最大迭代次数或同时满足以及其中，δ1和δ2均为公差参数阈值，均为公差参数阈值，均为各线性约束条带的当前迭代原始残差，为当前迭代的对偶残差。9.基于交替方向惩罚的频率一致宽带波束形成系统，其特征是，包括：第一模块，用来以h
h
(ω
n
)γ(ω
n
)h(ω
n
)最小为目标函数，结合无失真约束、主瓣约束、副瓣约束以及白噪声增益约束构建优化问题；其中，h(ω
n
)表示角频率ω
n
对应的权值向量，γ(ω
n
)表示伪相干矩阵；第二模块，用来利用交替方向惩罚法求解优化问题获得最优权值向量；所述第二模块进一步包括子模块：第一子模块，用来引入似实值助变量和对偶辅助向量重新表达优化问题；第二子模块，用来引入各线性约束的缩放对偶变量构造优化问题的增广拉格朗日函数；第三子模块，用来对权值向量、似实值助变量、对偶辅助向量、惩罚因子、缩放对偶变量进行交叉迭代求解；第四子模块，用来判断重复第三子模块直至达到收敛条件，此时的权值向量即最优的权值向量；第三模块，用来将最优的权值向量与导向矢量相乘结果作为频率一致宽带波束形成。

技术总结
本发明公开了基于交替方向惩罚的频率一致宽带波束形成方法及系统，该方法包括：步骤一，以h


技术研发人员：巩朋成 张峻嘉 徐开彦 李利荣 吴方祥 戴晗 李婕 贺章擎
受保护的技术使用者：湖北工业大学
技术研发日：2022.07.12
技术公布日：2022/11/1