CN103744072A - 基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法 - Google Patents

Abstract

本发明公开了一种基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，包括以下步骤：（1）将二维接收换能器阵列划分成若干个一级子阵，每个一级子阵作为一个基本单元，所有一级子阵组成二级子阵；结合二级子阵的波束方向图的旁瓣峰值和换能器的权重系数比，基于分布式并行子阵波束形成算法定义能量函数E(W,A)；（2）预先设定二级子阵的波束方向图的旁瓣峰值的目标阈值和换能器的权重系数比的阈值，利用模拟退火算法，进行二维接收换能器阵列的稀疏优化，得到需要开启的换能器数目的最小值。本发明的优化方法采用更少的换能器数目和更低的权重系数比，获得了相同的波束方向图性能，大幅度降低了系统的硬件复杂度和成本。

Description

基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法

技术领域

本发明涉及声纳阵列优化领域，具体涉及一种基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法。

背景技术

现有技术中的换能器阵列通常采用均匀分布或者周期分布的排列结构，再通过适当的窗函数对波束方向图进行加权处理，从而降低旁瓣峰值。

伴随海洋开发事业的迅猛发展，水下三维成像技术越来越追求高分辨率和实时操作，声纳换能器阵列的规模越来越庞大，换能器的数目越来越多，这必然会提高换能器阵列的设计成本、功耗和体积，因此，需要对换能器阵列的设计优化展开深入研究，稀疏阵列设计技术为解决该问题提供了一种有效的途径，并在雷达、声纳、无线通信和医学成像等领域获得越来越广泛的应用。

稀疏阵列是一种不等间距的阵列，通过在二维接收换能器的全阵列中去除掉一部分换能器，并对保留的换能器的位置和权重进行再次优化设计，从而降低换能器阵列的设计复杂度。

阵列稀疏通常会引起波束方向图旁瓣峰值的增大，如何在保证高稀疏率的同时有效控制旁瓣峰值，始终是稀疏阵列技术关注的热点。

模拟退火算法是一种启发式算法，由于其在模拟固体结晶的物理过程中，具备收敛速度快且具备全局最优解等特点，被广泛应用于稀疏阵列设计领域。

模拟退火算法在换能器稀疏阵列设计中，通过多次迭代寻找目标函数最小值的过程来模拟固体结晶，属于一种Markov链方法。在每次迭代过程中，如果新的配置参数引起目标函数减小，则接受当前配置并作为下一个Markov链的入口；相反，如果目标函数被增大，新配置参数要依据一定的概率分布被选择性地接受，概率大小取决于系统温度。伴随着迭代的不断进行，系统温度不断降低，直到冷却在某个最终状态，达到能量最小。

发明内容

本发明提供了一种基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，该优化方法与现有技术中的阵列稀疏算法相比较，采用更少的换能器数目和更低的权重系数比，获得了相同的波束方向图性能，大幅度降低了系统的硬件复杂度和成本。

一种基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，包括以下步骤：

（1）将二维接收换能器阵列划分成若干个一级子阵，每个一级子阵作为一个基本单元，所有一级子阵组成二级子阵；

为便于计算的进行，所述二维接收换能器阵列由N×N个间距为d的换能器构成，每个一级子阵由N_a×N_a个换能器构成，每个二级子阵由N_b×N_b个一级子阵构成，且N_a×N_b=N。

结合二级子阵的波束方向图的旁瓣峰值和换能器的权重系数比，基于分布式并行子阵波束形成算法定义能量函数E(W,A)；

以降低内存需求量和计算需求量为出发点，基于分布式并行子阵波束形成算法，并充分考虑波束方向图旁瓣峰值和换能器的权重系数比，进行能量函数的定义。

（2）预先设定二级子阵的波束方向图的旁瓣峰值目标阈值（一般为-22dB）和换能器的权重系数比的阈值（一般为3），利用模拟退火算法，进行二维接收换能器阵列的稀疏优化，得到需要开启的换能器数目的最小值。

模拟退火算法通过模拟固体退火结晶的原理，采用随机优化理论，相比现有技术，能够以更快的收敛速度获得更少的换能器的数目。

作为优选，所述能量函数E(W,A)的表达式如下：

$E(W,A)=l_1{\left(\underset{({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}(\frac{B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)}{B}-b_d)\right)}^2+l_2{A}^2+l_3{(R_0-R_d)}^2$

式中：B为分布式并行子阵波束形成算法中，二级子阵波束形成所得波束强度的最大值；

l₁，l₂，l₃分别代表每一项分量的权重系数；

是二级子阵波束形成后所得的最终波束强度；

b_d表示预期达到的旁瓣峰值；

A为二维接收换能器阵列中换能器的数目；

R_d表示预设的权重系数比值；

R₀表示当前权重系数矩阵计算所得的权重系数比值；

Ω表示单位矢量横纵坐标组合

取值的集合。

作为优选，所述Ω满足以下条件：

$\left\{\begin{array}{c}{\stackrel{\‾}{u}}_x^2+{\stackrel{\‾}{u}}_y^2>\mathrm{\λ}/D\\ \frac{B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)}{B}>b_d\end{array}\right.$

式中，D表示二维接收换能器阵列的孔径大小；

是二级子阵波束形成后所得的最终波束强度；

b_d表示预期达到的旁瓣峰值；

λ表示声纳回波信号的中心频率所对应的波长；

B为分布式并行子阵波束形成算法中，二级子阵波束形成所得波束强度的最大值。

作为优选，所述模拟退火算法的结束状态标准如下：

A(l-L+1)＝A(l),l＞L-1

其中，L为设定值，取值范围为[10,100]中的整数；

l为模拟退火算法的迭代次数；

A(l)为第l次迭代时，二维接收换能器阵列中换能器的数目。

作为优选，所述模拟退火算法中的初始温度T_start的取值≥1000。

本发明提出的系数阵列优化算法与其它算法相比较，采用更少的换能器数目和更低的权重系数比，获得了相同的波束方向图性能，采用模拟退火算法进行实时相控阵三维成像声纳系统的平面目标接收阵设计，能够有效减少换能器是数目，降低系统硬件复杂度和成本，具有重要的工程应用价值。

附图说明

图1为本发明目标阵列的示意图；

图2为本发明中子阵划分示意图；

图3为本发明中所采用的新方位角和仰视角的定义方式；

图4为本发明中二级子阵水平波束抽取示意图；

图5为本发明中的模拟退火算法的流程图。

具体实施方式

为了更为具体地描述本发明，下面结合附图和具体实施方式对本发明的方法做详细描述。

如图1所示，模拟退火算法以一个等间距均匀布阵的二维接收换能器阵列为优化目标，为简化数学表达形式又不失一般性，该二维接收换能器阵列采用方形布阵，由位于直角坐标系的xOy平面内的N×N个间距为d的换能器组成。

为了更加高效的应用模拟退火算法，图1中在方形的二维接收换能器阵列内划分出直径为D（D的大小通常与阵列边长相等）的圆形区域作为有效优化范围，圆形边界以外的换能器的权重系数均赋值为零，图1中空心圆表示权重系数为零的换能器，实心圆表示权重系数不为零的换能器。

将二维接收换能器阵列划分成若干个一级子阵，每个一级子阵作为一个基本单元，所有一级子阵组成二级子阵，如图2所示；

二级子阵为二维接收换能器阵列的一种虚拟表达方式，区别在于，二维接收换能器阵列以换能器为基本单元，而二级子阵以一级子阵为基本单元。

每个一级子阵由N_a×N_a个换能器组成，分别在水平Q_a个波束方向和垂直Q_a个波束方向进行波束形成，每个二级子阵由N_b×N_b个一级子阵组成，分别在水平Q_b个波束方向和垂直Q_b个波束方向进行波束形成。

每个一级子阵以二维接收换能器阵列的中心为坐标原点，那么第(m_a,n_a)号换能器的坐标r_m,n可由下式表示：

$r_{m,n}=(x_{m_a},y_{n_a},0)=((m_a-\frac{N_a+1}{2})d,(n_a-\frac{N_a+1}{2})d,0)$

式中，d为每行（或每列）中相邻两个换能器的间距；

m_a为一级子阵中各换能器的横坐标，取值为由1～N_a的自然数；

n_a为一级子阵中各换能器的纵坐标，取值为由1～N_a的自然数。

为减少相位偏移参数存储空间，DPS波束形成算法采用Palmese、Murino和Trucco提出的新方位角和仰视角定义（参见文献V.Murino and A.Trucco.Three-dimensional Image Generation and Processing in UnderwaterAcoustic Vision[J].Proceedings of IEEE.2000,88(12):1903～1948和M.Palmese and A.Trucco.Three-Dimensional Acoustic Imaging by Chirp ZetaTransform Digital Beamforming[J].IEEE Transactions on Instrumentationand Measurement.2009,58(7):2080～2086）。这种新方位角和仰视角定义方式被广泛应用于波束形成算法和稀疏阵列设计。

如图3所示，θ_a和θ_e分别表示波束方向的方位角和仰视角，此时波束方向的单位矢量u可表示如下：

$u=(u_x,u_y,u_z)=(\sin {\mathrm{\θ}}_{{\mathrm{ap}}_a},\sin {\mathrm{\θ}}_{{\mathrm{eq}}_a},\sqrt{{\cos }^2{\mathrm{\θ}}_{{\mathrm{ap}}_a}-{\sin }^2{\mathrm{\θ}}_{{\mathrm{eq}}_a}})$

式中：

表示一级子阵的任意一个波束方向的方位角，1≤p_a≤Q_a；

表示一级子阵的任意一个波束方向的仰视角，1≤q_a≤Q_a。

Q_a为一级子阵中的每行（或列）的波束数量。

如果使用v表示垂直于回波平面并指向声源的单位矢量，参考Vantrees在推导均匀矩形阵波束方向图时的描述（参见文献H.L.Vantrees.OptimumArray Processing.Part IV of Detection,Estimation,and Modulation Theory[M].New York:Wiley,2002），令v=(0,0,1)，则一级子阵波束形成的波束方向图|B1(W,u_x,u_y)|的表达式如下：

|B1(W,u_x,u_y)|

$=\left|\underset{m_a=1}{\overset{N_a}{\mathrm{\Σ}}}\underset{n_a=1}{\overset{N_a}{\mathrm{\Σ}}}{\mathrm{\ω}}_{(m_a,n_a)}{X}_{(m_a,n_a)}\left(k\right)\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}(r_{m,n}\·(v-u)))\right|$

$=\left|\underset{m_a=1}{\overset{N_a}{\mathrm{\Σ}}}\underset{n_a=1}{\overset{N_a}{\mathrm{\Σ}}}{\mathrm{\ω}}_{(m_a,n_a)}{X}_{(m_a,n_a)}\left(k\right)\exp (-j({\mathrm{\ψ}}_x(m_a-{\mathrm{\θ}}_{{\mathrm{ap}}_{\text{a}}})+{\mathrm{\ψ}}_y(n_a,{\mathrm{\θ}}_{{\mathrm{eq}}_a})))\right|$

$=\left|\underset{m_a=1}{\overset{N_a}{\mathrm{\Σ}}}\underset{n_a=1}{\overset{N_a}{\mathrm{\Σ}}}{\mathrm{\ω}}_{(m_a,n_a)}{X}_{(m_a,n_a)}\left(k\right)\exp (-j\frac{2\mathrm{\π}{f}_{0}d}{c}((m_a-\frac{N_a+1}{2})u_x+(n_a-\frac{N_a+1}{2})u_y))\right|$

式中：m_a为一级子阵中各换能器的横坐标，取值为由1～N_a的自然数；

n_a为一级子阵中各换能器的纵坐标，取值为由1～N_a的自然数；

表示换能器(m_a,n_a)的权重系数；

W表示换能器的权重系数矩阵；

λ表示声纳回波信号的中心频率所对应的波长；

f₀表示声纳回波信号的中心频率；

c代表声波在介质中的传播速度；

j表示虚部。

为每个换能器采集的原始数据进行DFT变换后的结果， $X_{(m_a,n_a)}\left(k\right)$ 的表达式如下：

$X_{(m_a,n_a)}\left(k\right)=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{x}_{(m_a,n_a)}\left(l\right)\exp (-j\frac{2\mathrm{\π}}{L}l\·k)$

式中，为换能器(m_a,n_a)采集的回波信号的原始数据；

(m_a,n_a)为一级子阵中各换能器的坐标；

k为DFT变换的线谱号；

L为DFT变换区间长度；

j表示虚部。

相位偏移参数表达式为：

${\mathrm{\ψ}}_x(m_a,{\mathrm{\θ}}_{{\mathrm{ap}}_a})=2\mathrm{\π}{f}_0\×\frac{(m_a-\frac{N_a+1}{2})d\sin {\mathrm{\θ}}_{{\mathrm{ap}}_a}}{c}$

${\mathrm{\ψ}}_y(n_a,{\mathrm{\θ}}_{{\mathrm{eq}}_a})=2\mathrm{\π}{f}_0\×\frac{(n_a-\frac{N_a+1}{2})d\sin {\mathrm{\θ}}_{{\mathrm{eq}}_a}}{c}$

由|B1(W,u_x,u_y)|的表达式可以看出，|B1(W,u_x,u_y)|满足以下条件：

$\left\{\begin{array}{c}B1(W,u_x,u_y)=B1(W,u_x,u_y)\\ B1(W,u_x,u_y)=B1(W,u_x,u_y)\end{array}\right.$

由这两个式子可以看出一级子阵的波束方向图的对称特性，因此波束方向的单位矢量u的横纵坐标取值范围可以限制如下：

u_x,u_y∈[0,1]。

基于每个换能器采集的原始数据进行DFT变换后的结果，分别针对一级子阵的每行换能器进行水平波束形成

然后利用水平波束形成结果，分别针对一级子阵中的每列换能器进行垂直波束形成 $\mathrm{bf}2({\mathrm{\θ}}_{{\mathrm{ap}}_a},{\mathrm{\θ}}_{{\mathrm{eq}}_a})$

对二级子阵进行波束形成时，首先对每个一级子阵的波束信号进行抽取，抽取分为水平抽取和垂直抽取，抽取方式分述如下：

3-1-1、在水平抽取过程中，从每个一级子阵水平方向的Q_a个波束面中选取与二级子阵欲成波束方向距离最近的一个平面。

图4为波束水平抽取示意图，如图4所示，实线表示一级子阵波束信号在水平方向的Q_a个波束平面，虚线表示二级子阵的波束方向。

以左侧第一个波束面为参考平面，二级子阵波束方向的水平标号为p_b，与参考平面的水平夹角为β；目标平面的水平标号为x，与参考平面水平夹角为α，水平抽取的过程是当β等于α时，求得x的数值，如下式所示：

$(p_b-1)\×\frac{{\mathrm{\θ}}_{\max }}{Q_b-1}=(x-1)\×\frac{{\mathrm{\θ}}_{\max }}{Q_a-1}$

θ_max表示二级子阵波束方向中，与参考平面夹角最大的波束方向与参考平面的夹角。

由此，可以得到目标平面的水平标号的表达式：

式中，

表示向下取整。

3-1-2、在垂直抽取过程中，从水平抽取所得的目标平面x中，选取与二级子阵欲成波束方向距离最近的一个波束，垂直抽取与水平抽取相类似，假定二级子阵波束方向的垂直标号为q_b，目标波束信号的垂直标号为y，则目标平面的垂直标号的表达式为：

经过波束抽取后，每个一级子阵的波束信号bf2(θ_ax,θ_ey)作为二级子阵的一个基元，参与二级子阵波束形成运算。

经过波束抽取后，来自于一级子阵的N_b×N_b个波束信号作为二级子阵的基元Beam(m_b,n_b)，重新编号为：

Beam(m_b,n_b)＝bf2(θ_ax,θ_ey)

二级子阵的波束方向与二维接收换能器阵列相同，同理可以推导出二级子阵波束形成的波束方向图

的表达式为：

$\left|B2\right({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y\left)\right|$

$=\left|\underset{m_b=1}{\overset{N_b}{\mathrm{\Σ}}}\underset{n_b=1}{\overset{N_b}{\mathrm{\Σ}}}\mathrm{Beam}(m_b,n_b)\exp (-j({\mathrm{\ψ}}_x(m_b,{\mathrm{\θ}}_{\mathrm{ap}})+{\mathrm{\ψ}}_y(n_b,{\mathrm{\θ}}_{\mathrm{eq}})))\right|$

$=\left|\underset{m_b=1}{\overset{N_b}{\mathrm{\Σ}}}\underset{n_b=1}{\overset{N_b}{\mathrm{\Σ}}}\mathrm{Beam}(m_b,n_b)\exp (-j\frac{2\mathrm{\π}{f}_0{N}_{a}d}{c}((m_b-\frac{N_b+1}{2}){\stackrel{\‾}{u}}_x+(n_b-\frac{N_b+1}{2}){\stackrel{\‾}{u}}_y)\right)|$

式中：m_b为二级子阵中各一级子阵的横坐标，取值为由1～N_b的自然数；

n_b为二级子阵中各一级子阵的纵坐标，取值为由1～N_b的自然数；

θ_ap表示二级子阵的任意一个波束方向的方位角，1≤p≤Q_b；

θ_eq表示二级子阵的任意一个波束方向的仰视角，1≤q≤Q_b。

Q_b为二级子阵中的每行（或列）的波束数量；

f₀表示声纳回波信号的中心频率；

c代表声波在介质中的传播速度；

j表示虚部。

二级子阵的相位偏移参数表达式为：

${\mathrm{\ψ}}_x(m_b,{\mathrm{\θ}}_{\mathrm{ap}})=2\mathrm{\π}{f}_0\×\frac{(m_b-\frac{N_b+1}{2})d\sin {\mathrm{\θ}}_{\mathrm{ap}}}{c}$

${\mathrm{\ψ}}_y(n_b,{\mathrm{\θ}}_{\mathrm{eq}})=2\mathrm{\π}{f}_0\×\frac{(n_b-\frac{N_b+1}{2})d\sin {\mathrm{\θ}}_{\mathrm{eq}}}{c}$

由

的表达式，可以得到

的对称性如下所示： $\left\{\begin{array}{c}B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)=B2(-{\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)\\ B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)=B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)\end{array}\right.$

基于分布式并行子阵波束形成算法的波束图表达式，并结合波束方向图的旁瓣峰值（Sidelobe Level Peak,SLP）和换能器的权重系数比（CurrentTaper Ratio,CTR），定义模拟退火算法的能量函数E(W,A)的表达式如下：

$E(W,A)=l_1{\left(\underset{({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}(\frac{B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)}{B}-b_d)\right)}^2+l_2{A}^2+l_3{(R_0-R_d)}^2$

式中：B为分布式并行子阵波束形成算法中，二级子阵波束形成所得波束强度的最大值，即 $B=\max \left(B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)\right);$

l₁，l₂，l₃分别代表每一项分量的权重系数，可以预先设定，通常l₁的取值范围为1000-100000，l₂和l₃的取值范围为0-1；

是二级子阵波束形成后所得的最终波束强度；

b_d表示预期达到的旁瓣峰值（SLP）；

A为二维接收换能器阵列中换能器的数目（即稀疏优化过程中开启的换能器的数目）；

R_d表示预设的权重系数比值（CTR）；

R₀表示当前权重系数矩阵计算所得的权重系数比值；

Ω表示单位矢量横纵坐标组合

取值的集合，该集合满足以下条件：

$\left\{\begin{array}{c}{\stackrel{\‾}{u}}_x^2+{\stackrel{\‾}{u}}_y^2>\mathrm{\λ}/D\\ \frac{B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)}{B}>b_d\end{array}\right.$

式中，D表示二维接收换能器阵列的孔径大小，也即图1中所示的圆形边界的直径值，该式表明

集合的取值范围排除波束方向图的主瓣区域。

当模拟退火算法的目标接收阵（最终优化得到的需要开启的换能器阵列）在Q×Q个波束方向进行波束形成时，根据u_x,u_y∈[0,1]以及

的取值范围的限定，可以选择Ω集合为一组离散坐标点组合，如下式所示：

${\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y\∈[\mathrm{0,1}/Q,\·\·\·,q/Q,\·\·\·,1],q\∈[0,Q]$

如图5所示，模拟退火算法过程主要包括以下几个步骤：

1）N×N二维接收换能器阵列参数初始化操作。首先给换能器的权重系数矩阵W初始化随机赋值为1或者0，初始化过程中，阵列圆形边界以外的换能器始终处于关闭状态，即赋值为0。由于当模拟退火算法Markov链的扰动配置引起能量函数增大时，新配置的参数要依据系统温度的玻尔兹曼概率分布被选择性地接受，因此，在初始化过程中，模拟退火算法的初始温度T_start取值应该足够大，以保证Markov链的第一个入口的扰动配置总是可以被接受。

2）在第l次迭代过程中，为了给系统引入微小的扰动配置，随机选择一个换能器，并将权重系数值赋值给中间变量ω_t；模拟退火算法采用随机排序的方式访问所有换能器，每个换能器在每次迭代过程中有且仅有一次被访问的机会，根据被访问的换能器的权重系数的不同，可以分别采取以下几种处理方式：

2-1、当被访问的换能器处于开启状态时（即权重系数不为0时），关闭该换能器并将其对应的权重系数赋值为0，此时权重系数矩阵W被改变，

2-1-1、如果重新计算所得的能量函数E(W,A)的数值相对于改变前减小，那么该换能器将保持关闭状态，并进入步骤3）；

2-1-2、如果重新计算所得的能量函数E(W,A)的数值相对于改变前增大，那么重新开启该换能器，并在权重系数矩阵W中引入以下扰动：

ω_(m,n)＝ω_t+unifrnd(-0.1,0.1)

式中，unifrnd(-0.1,0.1)是一个Matlab数学函数，表示任意选取-0.1和0.1之间的一个随机数。

2-2、当被访问的换能器处于关闭状态时（即权重系数等于0），该换能器将会按照以下的重生概率Pr(resurrection)被重新开启：

$\left\{\begin{array}{c}\mathrm{Pr}\left(\mathrm{resurrection}\right)=\\ l_4\·(1-\frac{4}{N\×N}(m-\frac{N+1}{2})(n-\frac{N+1}{2}))\end{array}\right.$

式中，l₄为重生概率的权重系数，取值为0-1；

m和n分别表示二维接收换能器阵列的换能器编号，1≤m≤N，1≤n≤N。

2-2-1、若换能器被重新开启，重新开启后的换能器会被随机赋予一个权重系数，并依此权重系数更新权重系数矩阵W。

2-2-2、若换能器没有被重新被开启，将会返回步骤2），继续处理下一个换能器。

2-3、完成步骤2-1和步骤2-2两步操作后，依据不同情况分别处理如下：

2-3-1、如果新的权重系数矩阵W_n减少了能量函数值，则接受该权重系数矩阵W_n，并将其作为下一个Markov链入口的权重系数矩阵W_l+1；

2-3-2、如果新的权重系数矩阵W_n增大了能量函数值，则根据系统温度T(l)的玻尔兹曼概率分布决定是否被接受为下一个Markov链入口的权重系数矩阵，玻尔兹曼概率分布Pr(W_l+1＝W_n)如下式所示：

$\mathrm{Pr}(W_{l+1}=W_n)=\left\{\begin{array}{cc}\exp \left(\frac{E_l-E_n}{\mathrm{bT}\left(l\right)}\right),& \mathrm{if}{E}_n>E_l\\ 1,& \mathrm{otherwise}\end{array}\right.$

式中，E_n表示新的能量函数值；

E_l表示第l次迭代后的能量函数值；

T(l)表示系统温度；

b为玻尔兹曼常数。

2-3-2-1、若新的权重系数矩阵W_n被接受，则依此权重系数更新权重系数矩阵W；

2-3-2-2、若新的权重系数矩阵W_n不被接受，则保持当前换能器的权重系数不变。

3）若所有换能器均被访问一次，则迭代次数l和系统温度T(l)将会按照以下数学表达式进行更新：

$T\left(l\right)=\left\{\begin{array}{cc}{T}_{\mathrm{start},}& l=1\\ 0.85\×T(l-1),& l>1\end{array}\right.$

否则返回步骤2），继续随机访问下一个换能器。

最终，模拟退火算法过程随着系统温度的降低，能量函数值E(W,A)取得最小值，系统达到结晶状态。

本发明方法定义模拟退火算法退火过程结束的标准如下式所示：

A(l-L+1)＝A(l),l＞L-1

由该式可知，在最后L次(设定值)迭代过程中，若目标阵列的换能器数量A不再减小，则退火过程结束，若没有达到该标准，则返回步骤2），继续进行退火处理，直至满足标准。

通常情况下，模拟退火算法迭代次数越多，能量函数的最终状态越稳定，因此在相控阵三维成像声纳系统的阵列稀疏化设计中，为了满足系统稳定性需求，L被选取为一个较大的数值。

本发明提出的系数阵列优化算法与其它算法相比较，采用更少的换能器数目和更低的权重系数比，获得了相同的波束方向图性能，采用模拟退火算法进行实时相控阵三维成像声纳系统的平面目标接收阵设计，能够有效减少换能器是数目，降低系统硬件复杂度和成本，具有重要的工程应用价值。

Claims (6)

1.一种基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，其特征在于，包括以下步骤：

（1）将二维接收换能器阵列划分成若干个一级子阵，每个一级子阵作为一个基本单元，所有一级子阵组成二级子阵；

结合二级子阵的波束方向图的旁瓣峰值和换能器的权重系数比，基于分布式并行子阵波束形成算法定义能量函数E(W,A)；

（2）预先设定二级子阵的波束方向图旁瓣峰值的目标阈值和换能器的权重系数比的阈值，利用模拟退火算法，进行二维接收换能器阵列的稀疏优化，得到需要开启的换能器数目的最小值。

2.如权利要求1所述的基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，其特征在于，所述能量函数E(W,A)的表达式如下：

$E(W,A)=l_1{\left(\underset{({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)\∈\mathrm{\Ω}}{\mathrm{\Σ}}(\frac{B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)}{B}-b_d)\right)}^2+l_2{A}^2+l_3{(R_0-R_d)}^2$

式中：B为分布式并行子阵波束形成算法中，二级子阵波束形成所得波束强度的最大值；

l₁，l₂，l₃分别代表每一项分量的权重系数；

是二级子阵波束形成后所得的最终波束强度；

b_d表示预期达到的旁瓣峰值；

A为二维接收换能器阵列中换能器的数目；

R_d表示预设的权重系数比值；

R₀表示当前权重系数矩阵计算所得的权重系数比值；

Ω表示单位矢量横纵坐标组合

取值的集合。

3.如权利要求2所述的基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，其特征在于，所述二维接收换能器阵列由N×N个间距为d的换能器构成，每个一级子阵由N_a×N_a个换能器构成，每个二级子阵由N_b×N_b个一级子阵构成，且N_a×N_b=N。

4.如权利要求3所述的基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，其特征在于，所述Ω满足以下条件：

$\left\{\begin{array}{c}{\stackrel{\‾}{u}}_x^2+{\stackrel{\‾}{u}}_y^2>\mathrm{\λ}/D\\ \frac{B2({\stackrel{\‾}{u}}_x,{\stackrel{\‾}{u}}_y)}{B}>b_d\end{array}\right.$

式中，D表示二维接收换能器阵列的孔径大小；

是二级子阵波束形成后所得的最终波束强度；

b_d表示预期达到的旁瓣峰值；

λ表示声纳回波信号的中心频率所对应的波长；

B为分布式并行子阵波束形成算法中，二级子阵波束形成所得波束强度的最大值。

5.如权利要求4所述的基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，其特征在于，所述模拟退火算法的结束状态标准如下：

A(l-L+1)＝A(l),l＞L-1

其中，L为设定值，取值范围为[10,100]中的整数；

l为模拟退火算法的迭代次数；

A(l)为第l次迭代时，二维接收换能器阵列中换能器的数目。

6.如权利要求5所述的基于模拟退火算法和分布式并行子阵波束形成算法的稀疏阵列优化方法，其特征在于，所述模拟退火算法中的初始温度T_start的取值≥1000。

Images