基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法

1.本发明涉及相控阵三维摄像声纳系统技术领域，具体涉及一种基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法。


背景技术：

2.相控阵三维摄像声纳系统使用一个窄带声脉冲透射整个水下场景，并利用一个二维均匀矩形换能器阵列接收场景中的背向散射回波信号产生阵列响应信号。这些阵列响应信号经由波束形成技术和实时图像处理技术处理可以得到高分辨率的水下三维图像。
3.为了提高图像的分辨率，同时保持一个较大的探测视角，二维均匀矩形换能器阵列通常具有较大的孔径，并且包含数千个阵元。此外，经由信号处理得到的波束图像中通常包含了上万个波束方向，这导致了巨大的计算负担。
4.传统的延时求和波束形成方法根据波束方向来调整阵元采样信号的时延并累加，能够精确的得到任意方向的波束结果，但其巨大的计算负担无法保证实时成像。采用频域直接波束形成方法能够将计算负担降低一个数量级，然而这仍然是难以承受的。因此，一些频域快速波束形成方法被提出，比如chirp-z变换(chirp zeta transform，czt)和非均匀傅里叶变换(nonuniform fast fourier transform，nufft)。这些快速波束形成方法以不同的方式利用了快速傅里叶变换的高效计算性能，其计算负担相比于频域直接波束形成方法降低了一到两个数量级。
5.然而，这些快速波束形成方法对于信号模型有着较为严格的要求，这在远场中是容易满足的。对于近场区域，通常采用菲涅尔近似模型以避免一些不满足要求的高阶参数，这导致了近场区域有效视野狭小。菲涅尔近似模型只能保证18
°
以内的波束形成精度，18
°
以外的波束结果存在较为明显的畸变。
6.专利文献cn109283536a公开了一种多波束测深声呐水体成像波束形成算法，在每次探测采样时间内，按照时间增益曲线对声波的传播损失进行补偿，经时间平均后得到当前探测水域背景噪声等级；对信号进行近场聚焦波束形成并根据当前背景噪声级，预估出当前快拍序号下的信源数目；对快拍数为1的信号向量进行协方差矩阵估计，通过对前后向平滑后的数据协方差矩阵进行矩阵重构，得到新的伪协方差矩阵；对伪协方差矩阵进行奇异值分解，利用常规波束形成输出结果配合阵列流型构造空间谱函数，得到多波束测深声呐水体成像结果。该技术波束成型方式的计算负担仍然大。


技术实现要素：

7.鉴于上述，本发明的目的是提供一种基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法该方法近场快速波束形成将完整的近场探测视野划分为多个子区域，分别构建每个子区域的信号模型，采用最小二乘加权和坐标轴旋转操作以保证模型精度，最后利用非均匀傅里叶变换(nufft)方法对每个子区域并行执行快速波束形成，对波束结果
重映射坐标后得到90
°
方向内高精度的波束图，其计算负担仍然保持在相同或更低的水平。
8.实施例提供的一种基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，包括以下步骤：
9.(1)将二维换能器阵列的观测范围的坐标系以预设的旋转角度进行旋转，并将旋转后的观测范围均匀划分为多个子区域，其中，旋转角度以顺时针方向为正方向，并为45
°
的奇数倍；
10.(2)简化每个子区域中换能器阵元的近场时延参数，利用最小二乘法根据简化的近场时延参数构建加权的近场时延参数，并对加权的近场时延参数进行求解以确定最优的加权系数；
11.(3)根据二维换能器阵列产生的阵列响应信号和带有加权系数的波束信号构建伴随非均匀傅里叶变换的算子，根据伴随非均匀傅里叶变换的算子并行执行非均匀傅里叶变换计算得到每个子区域内的波束；
12.(4)重映射每个子区域内的波束的坐标，得到高精度的波束图。
13.优选地，所述二维换能器阵列为包含n
×
n个换能器阵元的二维矩形阵列，该二维矩形阵列位于三维空间的xoy平面，二维矩形阵列的边分别与x和y坐标轴平行，第(n1,n2)个换能器阵元的坐标表示为(x
n1
,y
n2
)；
14.波束的聚焦距离r0为近场区域，所述观测范围由两个方向参数θa和θe表示，它们分别平行于x和y轴，并且取值范围为(-90
°
,90
°
)。
15.优选地，所述观测范围内具有p
×
q个波束需要被生成，对于第(p,q)个波束，其波束方向为(θ
ap
,θ
eq
)，波束信号b(r0,p,q)表示为：
[0016][0017]
其中，f
l
表示信号频率，r0表示波束在近场区域的聚焦距离，s
l
(n1,n2)表示第(n1,n2)个换能器阵元的频域采样数据，n表示每个方向的换能器数量，τ(r0,p,q,n1,n2)表示第(n1,n2)个换能器阵元对于波束方向(θ
ap
,θ
eq
)的时延，表示为：
[0018][0019]
其中，c表示声速，(x
n1
,y
n2
)表示第(n1,n2)个换能器阵元的坐标。
[0020]
优选地，步骤(1)中，将旋转后的观测范围均匀划分为3
×
3或4
×
4的子区域，具体操作如下：
[0021]
(1-1)将观测范围的坐标系以预设的旋转角度φ进行旋转后，建立的新方向参数θ
′a和θ
′e如下：
[0022][0023]
其中，限制方向参数θ
′a和θ
′e的取值范围为(-90
°
,90
°
)；
[0024]
(1-2)将方向参数θ
′a和θ
′e在取值范围内按照正弦域等分为三段或四段，以u和v表示每个子区域的索引，即：
[0025][0026]
或，
[0027][0028]
对于第(u,v)个子区域，其位置由sinθ
′a和sinθ
′e对应的取值范围确定。
[0029]
优选地，步骤(2)中，简化每个子区域中换能器阵元的近场时延参数，简化的近场时延参数τ
sp
(r0,p,q,n1,n2)表示为：
[0030][0031]
其中，r0表示波束在近场区域的聚焦距离，c表示声速，sinθ
′
ap
和sinθ
′
eq
是新坐标系中第(p,q)个波束方向对应的正弦域参数，且
[0032][0033][0034]
其中，sinθ
′
au
和sinθ
′
ev
是第(u,v)个子区域的中心，即对于3
×
3子区域有
[0035][0036]
或对于4
×
4子区域有
[0037][0038]
优选地，步骤(2)中，利用最小二乘法根据简化的近场时延参数构建加权的近场时延参数τ
ls
(r0,p,q,n1,n2)，表示为：
[0039][0040]
其中，g(n1,n2,u,v)＝[3ρ2(n1,n2,u,v)-v(n1,n2)]/2r0；
[0041]
根据最小二乘法对加权的近场时延参数τ
ls
(r0,p,q,n1,n2)求解得到加权系数k
1-k6的最优值以使子区域内时延误差最小。
[0042]
优选地，步骤(3)包括：
[0043]
(3-1)将p
×
q个波束放入到对应的子区域中，每个子区域内包含pu×qv
个波束，并用(pu,qv)索引，则对应波束的方向参数表示为：
[0044][0045][0046][0047]
(3-2)对于第(pu,qv)个波束，根据二维换能器阵列产生的响应信号和带有加权系数的波束信号构建伴随非均匀傅里叶变换的算子b(r0,pu,qv)，表示为：
[0048][0049]
其中，s(n1,n2)表示为响应信号，ξ(n1,n2,u,v)表示与(pu,qv)无关的加权系数，ωa(n1,n2,u,v)和ωe(n1,n2,u,v)表示非均匀节点坐标，且
[0050][0051][0052][0053]
其中，f
l
表示信号频率，r0表示波束在近场区域的聚焦距离，c表示声速；
[0054]
(3-2)预计算ξ(n1,n2,u,v)、ωa(n1,n2,u,v)和ωe(n1,n2,u,v)，以及非均匀傅里叶变换方法所需要的补偿系数和插值系数；
[0055]
(3-3)根据步骤(3-2)的计算结果，对伴随非均匀傅里叶变换的算子并行进行非均匀傅里叶变换计算得到每个子区域内的波束，即确定波束的方向参数和
[0056]
优选地，步骤(4)中，对于子区域内的第(pu,qv)个波束进行重映射，得到的重映射坐标表示为：
[0057][0058]
基于重映射坐标得到原坐标系下的波束图。
[0059]
与现有技术相比，本发明具有以下有益的技术效果：
[0060]
本发明提供的近场快速波束形成方法将二维换能器阵列的观测范围划分为多个子区域，利用坐标系旋转和最小二乘方法的加权保证了子区域内所有波束方向时延参数的精度；使用非均匀傅里叶变换在每个子区域内实现快速波束形成，以保持低的计算负担；所有子区域的快速波束形成可以并行执行，有利于简化计算结构。因此，本发明适用于相控阵三维摄像声纳系统的近场区域探测，提高了三维摄像声纳系统的成像质量和探测性能。
附图说明
[0061]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。
[0062]
图1为实施例提供的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法的流程示意图；
[0063]
图2为实施例提供的相控阵三维摄像声纳系统换能器阵列的近场信号模型示意图；
[0064]
图3为实施例提供的观测范围的坐标系旋转和子区域划分示意图；
[0065]
图4为实施例提供的采用本发明的近场快速波束形成方法得到的波束图；
[0066]
图5为实施例提供的基于传统菲涅尔近似模型的波束图；
[0067]
图6为实施例提供的不同阵列参数下传统菲涅尔近似模型和本发明的计算负担对比。
具体实施方式
[0068]
为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。
[0069]
为实现对近场区域的二维换能器阵列的快速成像，实施例提供了一种基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法。如图1所示，实施例提供的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，包括以下步骤：
[0070]
步骤1，将二维换能器阵列的观测范围的坐标系以预设的旋转角度进行旋转，并将旋转后的观测范围均匀划分为多个子区域。
[0071]
实施例中，假设一个具体的相控阵三维摄像声纳系统中的换能器阵列为一个50
×
50的二维均匀矩形阵列，阵元间距为λ/2，λ为声波的波长。如图2所示，阵列中心位于坐标原点，近场聚焦距离r0为1m，探测范围为z轴正半轴(φ＜90
°
)。在聚焦距离上设置4个方位为(0
°
,0
°
),(18
°
,0
°
),(-25
°
,-25
°
),(-34.5
°
,-55
°
)的散射目标。
[0072]
每个阵元用两个方向参数n1和n2索引，其坐标用(x
n1
,y
n2
)表示。波束方向由两个方向参数θa和θe表示，它们分别平行于x和y轴。假设在每个波束方向上需要计算200个波束以保证图像的高分辨率。
[0073]
步骤1-1，将观测范围的坐标系以预设的旋转角度进行旋转，并建立新的方向参数。
[0074]
实施例中，当将观测范围的坐标系顺时针旋转45
°
，建立新的方向参数θ
′a和θ
′e如下：
[0075][0076]
其中，限制方向参数θ
′a和θ
′e的取值范围为(-90
°
,90
°
)。
[0077]
原始的观测范围的坐标系以θa和θe的正弦为刻度，由于其取值范围为(-90
°
,90
°
)，θa和θe的正弦是单调的，不存在重叠。
[0078]
可以看到旋转前和旋转后的观测范围存在差异，由于观测范围φ＜90
°
区域是圆形的，处于两个坐标系中的重叠区域，因此坐标系旋转不会造成视野的缺失。在θa和θe方向设置200个波束与在θ
′a和θ
′e方向设置200个波束具有相同的波束密度。
[0079]
所有的快速波束形成方法都建立在矩形观测范围上以便于分离坐标参数，尽管少量的结果是无意义的，但这是必要的。
[0080]
步骤1-2，将方向参数θ
′a和θ
′e在取值范围内按照正弦域等分为四段，即将旋转后的观测区域均匀划分为3
×
3或4
×
4的子区域，以4
×
4子区域为例，旋转后的坐标系和子区域如图3所示。
[0081]
实施例中，以u和v表示每个子区域的索引，即：
[0082][0083]
并且每段中包含50个波束方向，一个子区域中具有50
×
50个波束方向。
[0084]
步骤2，简化每个子区域的近场时延参数，并利用最小二乘方法对简化参数加权。
[0085]
实施例中，简化每个子区域中换能器阵元的近场时延参数，利用最小二乘法根据简化的近场时延参数构建加权的近场时延参数，并对加权的近场时延参数进行求解以确定最优的加权系数。
[0086]
步骤2-1，对每个子区域中换能器阵元的近场时延参数进行简化。
[0087]
对于传统菲涅尔近似，其近场时延参数τ
fa
(r0,p,q,n1,n2)表示为：
[0088][0089]
它是足够简洁的，但只在范围φ＜18
°
内有效。本方法提出了更精确的时延参数，在新坐标系中，简化的近场时延参数τ
sp
(r0,p,q,n1,n2)表示为：
[0090][0091]
其中，r0表示波束在近场区域的聚焦距离，c表示声速，sinθ
′
ap
和sinθ
′
eq
是新坐标系中第(p,q)个波束方向对应的正弦域参数，且
[0092][0093][0094]
其中sinθ
′
au
和sinθ
′
ev
是第(u,v)个子区域的中心，即：
[0095][0096]
步骤2-2，利用最小二乘方法对简化参数加权并求解。
[0097]
简化的近场时延参数τ
sp
(r0,p,q,n1,n2)在子区域中心处具有最高精度，在子区域角落处的误差较大，为均衡子区域内所有方向的误差，构建加权的近似时延参数τ
ls
(r0,p,q,n1,n2)，表示为：
[0098][0099]
其中，g(n1,n2,u,v)＝[3ρ2(n1,n2,u,v)-v(n1,n2)]/2r0，根据最小二乘法对加权的近场时延参数τ
ls
(r0,p,q,n1,n2)可以得到k
1-k6的最优值使子区域内时延的总均方误差最小。
[0100]
本实施例中，部分子区域中的k
1-k6的最优值如表1所示。
[0101]
表1
[0102][0103]
步骤3，对二维换能器阵列产生的阵列响应信号执行nufft方法，并行计算每个子区域内的波束。
[0104]
实施例中，根据二维换能器阵列产生的响应信号和带有加权系数的波束信号构建伴随非均匀傅里叶变换的算子，根据伴随非均匀傅里叶变换的算子并行执行非均匀傅里叶变换计算得到每个子区域内的波束。
[0105]
具体地，步骤3包括：
[0106]
步骤3-1，将每个子区域内的波束信号构建伴随nufft算子的形式。
[0107]
实施例中，将p
×
q个波束放入到对应的子区域中，每个子区域内包含pu×qv
个波束，基于每个子区域内的波束索引(pu,qv)来表示波束方向，其中0≤pu≤49和0≤qv≤49。波束方向表示为和其中并且
[0108][0109]
对于第(pu,qv)个波束，根据二维换能器阵列产生的响应信号和带有加权系数的波束信号构建伴随非均匀傅里叶变换的算子b(r0,pu,qv)，表示为：
[0110][0111]
其中，s(n1,n2)表示为响应信号，ξ(n1,n2,u,v)是与(pu,qv)无关的加权系数，ωa(n1,n2,u,v)和ωe(n1,n2,u,v)是非均匀节点坐标，且
[0112][0113][0114][0115]
步骤3-2，预计算ξ(n1,n2,u,v)、ωa(n1,n2,u,v)和ωe(n1,n2,u,v)，以及nufft方法所需要的补偿系数和插值系数。
[0116]
实施例中，每个换能器阵元信号需要4
×
4个插值系数，每个波束结果需要计算一个补偿系数，这可以通过一些工具箱预先计算并提前储存，不需要实时计算。
[0117]
步骤3-3，根据步骤3-2的计算结果，对伴随nufft的算子并行进行nufft计算得到每个子区域内的波束，即确定波束的方向参数。
[0118]
实施例中，对每个子区域可以并行计算，以保持低计算负担。
[0119]
步骤4，重映射波束的坐标，得到高精度的波束图。
[0120]
实施例中，对于子区域内的第(pu,qv)个波束，针对以顺时针旋转45
°
，重映射坐标为：
[0121][0122]
基于重映射坐标得到原坐标系下的波束图，如图4所示。
[0123]
作为对比，图5给出了基于菲涅尔近似模型的波束图结果。
[0124]
表2给出了本实施例中，利用菲涅尔近似和本方法得到的四个散射体对应的波束方向的强度。
[0125]
表2
[0126][0127]
由表2可以看出，本发明提出方法对φ＜90
°
范围内的波束结果都保持了高波束强度，而菲涅尔近似模型仅在φ≤18
°
范围内的波束强度较高。
[0128]
对于不同的阵列参数，4
×
4子区域和3
×
3子区域的计算负担是不同的，假设一个n
×
n的二维均匀传感器阵列需要计算4n
×
4n个波束，所需的在线实操作数量如图6所示。可以看出4
×
4子区域具有和菲涅尔近似模型几乎相同的计算负担，而3
×
3子区域在一些情况下具有更低的计算量，相比于频域直接波束形成减少了一到两个数量级的计算负担。因此，本发明提出的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法在保持低计算负担的同时具有接近90
°
的视野范围，提高了三维摄像声纳系统在近场区域的成像质量和探测性能。
[0129]
以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

技术特征：
1.一种基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，包括以下步骤：(1)将二维换能器阵列的观测范围的坐标系以预设的旋转角度进行旋转，并将旋转后的观测范围均匀划分为多个子区域，其中，旋转角度以顺时针方向为正方向，并为45
°
的奇数倍；(2)简化每个子区域中换能器阵元的近场时延参数，利用最小二乘法根据简化的近场时延参数构建加权的近场时延参数，并对加权的近场时延参数进行求解以确定最优的加权系数；(3)根据二维换能器阵列产生的响应信号和带有加权系数的波束信号构建伴随非均匀傅里叶变换的算子，根据伴随非均匀傅里叶变换的算子并行执行非均匀傅里叶变换计算得到每个子区域内的波束；(4)重映射每个子区域内的波束的坐标，得到高精度的波束图。2.如权利要求1所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，所述二维换能器阵列为包含n
×
n个换能器阵元的二维矩形阵列，该二维矩形阵列位于三维空间的xoy平面，二维矩形阵列的边分别与x和y坐标轴平行，第(n1,n2)个换能器阵元的坐标表示为(x
n1
,y
n2
)；波束的聚焦距离r0为近场区域，所述观测范围由两个方向参数θ
a
和θ
e
表示，它们分别平行于x和y轴，并且取值范围为(-90
°
,90
°
)。3.如权利要求1所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，所述观测范围内具有p
×
q个波束需要被生成，对于第(p,q)个波束，其波束方向为(θ
ap
,θ
eq
)，波束信号b(r0,p,q)表示为：其中，f
l
表示信号频率，r0表示波束在近场区域的聚焦距离，s
l
(n1,n2)表示第(n1,n2)个换能器阵元的频域采样数据，n表示每个方向的换能器数量，τ(r0,p,q,n1,n2)表示第(n1,n2)个换能器阵元对于波束方向(θ
ap
,θ
eq
)的时延，表示为：其中，c表示声速，(x
n1
,y
n2
)表示第(n1,n2)个换能器阵元的坐标。4.如权利要求1所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，步骤(1)中，将旋转后的观测范围均匀划分为3
×
3或4
×
4的子区域，具体操作如下：(1-1)将观测范围的坐标系以预设的旋转角度φ进行旋转后，建立的新方向参数θ
′
a
和θ
′
e
如下：其中，限制方向参数θ
′
a
和θ
′
e
的取值范围为(-90
°
,90
°
)；(1-2)将方向参数θ
′
a
和θ
′
e
在取值范围内按照正弦域等分为三段或四段，以u和v表示每
个子区域的索引，即：或，对于第(u,v)个子区域，其位置由sinθ
′
a
和sinθ
′
e
对应的取值范围确定。5.如权利要求1所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，步骤(2)中，简化每个子区域中换能器阵元的近场时延参数，简化的近场时延参数τ
sp
(r0,p,q,n1,n2)表示为：其中，r0表示波束在近场区域的聚焦距离，c表示声速，sinθ
′
ap
和sinθ
′
eq
是新坐标系中第(p,q)个波束方向对应的正弦域参数，且第(p,q)个波束方向对应的正弦域参数，且其中，sinθ
′
au
和sinθ
′
ev
是第(u,v)个子区域的中心。6.如权利要求1或5所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，步骤(2)中，利用最小二乘法根据简化的近场时延参数构建加权的近场时延参数τ
ls
(r0,p,q,n1,n2)，表示为：其中，g(n1,n2,u,v)＝[3ρ2(n1,n2,u,v)-v(n1,n2)]/2r0；根据最小二乘法对加权的近场时延参数τ
ls
(r0,p,q,n1,n2)求解得到加权系数k
1-k6的最优值以使子区域内时延误差最小。7.如权利要求1或6所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，步骤(3)包括：(3-1)将p
×
q个波束放入到对应的子区域中，每个子区域内包含p
u
×
q
v
个波束，并用(p
u
,q
v
)索引，则对应波束的方向参数表示为：
(3-2)对于第(p
u
,q
v
)个波束，根据二维换能器阵列产生的响应信号和带有加权系数的波束信号构建伴随非均匀傅里叶变换的算子b(r0,p
u
,q
v
)，表示为：其中，s(n1,n2)表示为响应信号，ξ(n1,n2,u,v)表示与(p
u
,q
v
)无关的加权系数，ω
a
(n1,n2,u,v)和ω
e
(n1,n2,u,v)表示非均匀节点坐标，且坐标，且坐标，且其中，f
l
表示信号频率，r0表示波束在近场区域的聚焦距离，c表示声速；(3-2)预计算ξ(n1,n2,u,v)、ω
a
(n1,n2,u,v)和ω
e
(n1,n2,u,v)，以及非均匀傅里叶变换方法所需要的补偿系数和插值系数；(3-3)根据步骤(3-2)的计算结果，对伴随非均匀傅里叶变换的算子并行进行非均匀傅里叶变换计算得到每个子区域内的波束，即确定波束的方向参数和8.如权利要求1所述的基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，其特征在于，步骤(4)中，对于子区域内的第(p
u
,q
v
)个波束进行重映射，得到的重映射坐标表示为：基于重映射坐标得到原坐标系下的波束图。

技术总结
本发明公开了一种基于子区域模型和非均匀傅里叶变换的近场快速波束形成方法，包括以下步骤：将观测范围的坐标系顺时针旋转45


技术研发人员：蒋荣欣 王斐 刘雪松 辜博轩 周凡
受保护的技术使用者：浙江大学
技术研发日：2022.07.08
技术公布日：2022/11/1