CN103323831A - 一种基于czt变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法 - Google Patents

Abstract

本发明公开了一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，包括以下步骤：分别计算每个换能器接收到的声波中，与载波频率f_k对应的第k个DFT变换系数；利用所得变换系数，构建矩阵C；对应方位角方向，对矩阵C的每一列，进行一维离散卷积运算，得到P×N的矩阵MD；对应仰视角方向，对矩阵MD的每一行，进行一维离散卷积运算，得到P×Q的矩阵ME；利用矩阵ME，构建矩阵MF；计算矩阵MF中每个元素的模，得到三维摄像声纳波束。本发明在FFT运算过程中识别和去除不必要的运算操作，有效降低了三维摄像声纳波束形成方法所需的运算量。

Description

一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法

技术领域

本发明涉及相控阵三维摄像声纳技术领域，具体涉及一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法。 

背景技术

相控阵三维摄像声纳系统包含几千个接收换能器，这些换能器组成二维面阵，利用二维面阵接收到的声音信号以及波束形成算法生成的上万个波束信号，获得三维摄像声纳系统的图像。 

为了获得较高的距离分辨率和图像刷新率，相控阵三维摄像声纳系统必须具备强大的实时波束形成处理能力，这对水下信号处理机的实时运算性能、功耗和体积等都提出了非常高的要求，实时声纳波束形成算法所需要的巨大运算量是相控阵三维摄像声纳系统应用推广的重要瓶颈之一。 

基于Chirp Zeta变换（Chirp Zeta Transform，CZT）的波束形成算法是一种利用Chirp Zeta变换的频域波束形成算法，在具体的Chirp Zeta变换运算过程中，基于Chirp Zeta变换的波束形成算法能够使用快速傅里叶变换实现算法的加速，因而相对于传统的时域延时-相加波束形成算法和频域相移波束形成算法具有更高的运算效率。基于Chirp Zeta变换的频域波束形成算法同样能够被用于三维摄像声纳系统中的二维均匀矩形阵列中。 

相控阵三维摄像声纳系统中的M×N换能器阵列如图1所示，换能器阵列位于xoy平面，u为波束方向角，定义θ_a和θ_e分别为u对应方位的方位角和仰视角，其具体定义如图1所示。 

在图1所示的二维方位定义下，波束方向的单位向量u可以表示为 

$u=(\sin {\mathrm{\θ}}_a,\sin {\mathrm{\θ}}_e,\sqrt{{\cos }^2{\mathrm{\θ}}_a-{\sin }^2{\mathrm{\θ}}_e})$

以p和q(1≤p≤P；1≤q≤Q)分别表示方位角索引号和仰视角索引号；θ_ap 和θ_eq分别为第(p，q)号波束方向对应的方位角和仰视角。 

远场条件下，波束信号b(t，θ_ap，θ_eq)的频域表达式B(k，θ_ap，θ_eq)可以表示为： 

$B(k,{\mathrm{\θ}}_{\mathrm{ap}},{\mathrm{\θ}}_{\mathrm{eq}})=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{w}_{m,n}{S}_{m,n}\left(k\right)\exp (-j2\mathrm{\π}{f}_k(\frac{(m-1)d\sin {\mathrm{\θ}}_{\mathrm{ap}}}{c}+\frac{(n-1)d\sin {\mathrm{\θ}}_{\mathrm{eq}}}{c}))$

其中，k为DFT变换的频率索引号；w_m，n为第(m，n)个换能器的权重；S_m，n(k)为第(m，n)号换能器采集信号的长度为K的DFT系数；f_k 为第k个DFT谱线对应的信号频率，即f_k=k f_s/K（f_s为采样频率），k=0，1，...，K-1；d为换能器的水平间距（或垂直间距），c为声速。 

对于相控阵三维摄像声纳系统中的二维均匀矩形面阵，基于Chirp Zeta变换的波束形成算法表达式为 

$B(k,{\mathrm{\θ}}_{\mathrm{ap}},{\mathrm{\θ}}_{\mathrm{eq}})=$

$W_a^{\frac{-{(p-1)}^2}{2}}{W}_e^{\frac{-{(q-1)}^2}{2}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{w}_{m,n}{S}_{m,n}\left(k\right)A_a^{m-1}{A}_e^{n-1}\×W_a^{\frac{-{(m-1)}^2}{2}}{W}_e^{\frac{-{(n-1)}^2}{2}}{W}_a^{\frac{{(p-m)}^2}{2}}{W}_e^{\frac{{(q-n)}^2}{2}}$

其中， 

$\left\{\begin{array}{c}{A}_{\text{a}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ai}})\\ W_a=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_a\right)\\ \mathrm{\&Delta;}{s}_a=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{af}}-\sin {\mathrm{\&theta;}}_{\mathrm{ai}}}{P-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ai}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{af}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

$\left\{\begin{array}{c}{A}_{\text{e}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ei}})\\ W_e=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_e\right)\\ \mathrm{\&Delta;}{s}_e=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{ef}}-\sin {\mathrm{\&theta;}}_{\mathrm{ei}}}{Q-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ei}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{ef}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

这两式中，θ_ai和θ_af分别表示波束扫描范围的起始方位角和终止方位 角，θ_ei和θ_ef分别表示波束扫描范围的起始仰视角和终止仰视角，P为波束扫描范围内的方向角数目，Q为波束扫描范围内的仰视角数目。 

θ_ap和θ_eq的取值由下式决定： 

sinθ_ap＝sinθ_ai+Δs_a(p-1),p＝1,...,P 

sinθ_eq＝sinθ_ei+Δs_e(q-1),q＝1,...,Q 

定义两个M×N矩阵E(k)和矩阵D，矩阵E(k)中的每个元素为： 

$E(k,m,n)=w_{m,n}{S}_{m,n}\left(k\right)A_a^{m-1}{A}_e^{n-1}\&times;W_a^{-{(m-1)}^2/2}{W}_e^{-{(n-1)}^2/2}$

矩阵D中的每个元素为： 

$D(m,n)=W_a^{m^2/2}{W}_e^{n^2/2}$

B(k，θ_ap，θ_eq)可以利用矩阵E(k)和矩阵D的二维离散卷积运算实现，因此，能够借助于二维FFT算法降低运算量，相对其他常规的波束形成方法具有明显的优势。 

分裂基FFT算法是一种计算DFT变换的快速算法，相对于其他固定基FFT算法（Fixed radix FFT），分裂基FFT算法具有最小的运算量。 

对于一个长度为N_T的一维序列（N_T为2的r次幂），其按时间抽取基-2/4类型的分裂基FFT算法可以描述为： 

$X\left(2k\right)=\underset{n=0}{\overset{N_T/2-1}{\mathrm{\&Sigma;}}}[x\left(n\right)+x(n+N_T/2)]W_{N_T}^{2\mathrm{nk}}$

$X(4k+1)=\underset{n=0}{\overset{N_T/4-1}{\mathrm{\&Sigma;}}}\left\{\right[x\left(n\right)-x(n+N_T/2)]$

$\text{-j[x}(n+N_T/4)-x(n+3N_T/4)\left]\right\}{W}_{N_T}^n{W}_{N_T}^{4\mathrm{nk}}$

$X(4k+3)=\underset{n=0}{\overset{N_T/4-1}{\mathrm{\&Sigma;}}}\left\{\right[x\left(n\right)-x(n+N_T/2)]$

$\text{+j[x}(n+N_T/4)-x(n+3N_T/4)\left]\right\}{W}_{N_T}^{3n}{W}_{N_T}^{4\mathrm{nk}}$

其中， 

$W_{N_T}^{\mathrm{nk}}=\cos \left(\frac{2\mathrm{\&pi;}\&CenterDot;\mathrm{nk}}{N_T}\right)-j\&CenterDot;\sin \left(\frac{2\mathrm{\&pi;}\&CenterDot;\mathrm{nk}}{N_T}\right)$

因此，一个长度为N_T的DFT运算可以被分解为一个长度为N_T/2的DFT 运算和两个长度为N_T/4的DFT运算，这一分解过程逐级进行，完成整个DFT运算。 

此外，在实际应用中，为了进一步提高FFT运算的算法效率，可以根据实际所需的输入数据和输出数据对FFT运算过程进行优化，去除多余的数据路径和相应的运算操作，从而进一步提高FFT的运算效率，这种去除FFT过程中多余运算的技术就是剪除FFT技术。 

在三维摄像声纳系统中，基于二维Chirp Zeta变换的波束形成算法在利用二维FFT运算矩阵的二维卷积时，需要对矩阵进行补0扩展。例如，对于一个64×64的均匀矩形换能器阵列，当需要运算80×80个波束信号时，至少需要将64×64矩阵C(k)扩展成大小为(64＋80－1)×(64＋80－1)的矩阵。 

考虑到FFT具体实现时，通常要求FFT的长度为2的整数次幂，因此，应该选择将C(k)补0扩展为256×256的矩阵，然后再运算经过扩展后的矩阵C(k)的二维FFT，此时，FFT运算的输入数据中很大一部分数据为0；同时，由于真正需要的二维卷积输出数据只有80×80个，卷积运算所需的IFFT（或FFT）也只需要运算部分结果。 

如果能够根据实际所需的输入数据和输出数据对FFT进行优化，去除多余的数据路径和运算操作，可以进一步提高FFT的运算效率，从而减少实时波束形成算法所需的运算量。 

发明内容

本发明提供了一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，该方法依次对方位角方向和仰视角方向进行一维离散卷积运算，并且一维离散卷积运算通过一种剪除的分裂基FFT方法实现，在FFT运算过程中识别和去除不必要的运算操作，有效降低了三维摄像声纳波束形成方法所需的运算量。 

一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，包括以下步骤： 

（1）分别计算每个换能器接收到的声波中，与载波频率f_k对应的第k个DFT变换系数S_m,n(k)，得到M×N个S_m,n(k)； 

其中，m为换能器的水平索引号，m=1，2，...，M；n为换能器的垂直索引号，n=1，2，...，N；k=0，1，...，K-1，K为DFT变换的长度。 

三维摄像声纳系统中的换能器也称水听器，用于将接收到的声波信号转换为电信号，M和N分别为换能器阵列的列数和行数；M×N个换能器在水平和垂直方向上都以相同的间距均匀排列，且水平间距与垂直间距等。 

（2）利用M×N个S_m，n(k)，构建矩阵C，矩阵C中的每个元素为： 

$C(k,m,n)=w_{m,n}{S}_{m,n}\left(k\right)A_a^{m-1}{A}_e^{n-1}\&times;W_a^{-{(m-1)}^2/2}{W}_e^{-{(n-1)}^2/2}$

其中， 

$\left\{\begin{array}{c}{A}_{\text{a}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ai}})\\ W_a=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_a\right)\\ \mathrm{\&Delta;}{s}_a=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{af}}-\sin {\mathrm{\&theta;}}_{\mathrm{ai}}}{P-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ai}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{af}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

$\left\{\begin{array}{c}{A}_{\text{e}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ei}})\\ W_e=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_e\right)\\ \mathrm{\&Delta;}{s}_e=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{ef}}-\sin {\mathrm{\&theta;}}_{\mathrm{ei}}}{Q-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ei}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{ef}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

其中，w_m,n表示每个换能器的权重；j为虚数单位；f_k为载波频率；d为换能器的水平间距（由于换能器的水平间距与垂直间距相等，所以此处也可以为垂直间距）；c为水下声速；θ_ai表示声纳系统波束扫描范围的起始方位角；θ_af表示声纳系统波束扫描范围的终止方位角；P为波束扫描范围内的方位角数目；θ_ei表示声纳系统波束扫描范围的起始仰视角；θ_ef表示声纳系统波束扫描范围的终止仰视角；Q为波束扫描范围内的仰视角数目。 

（3）对应方位角方向，对矩阵C的每一列，进行一维离散卷积运算，保留一维离散卷积运算中的P个结果，得到P×N的矩阵MD。 

所述的一维离散卷积运算的公式如下： 

$F_1(p,n)=\underset{m-1}{\overset{M}{\mathrm{\&Sigma;}}}C(m,n)\&CenterDot;W_a^{\frac{{(p-m)}^2}{2}}$

其中，p的取值范围为1～P。 

为了提高一维离散卷积的运算效率，优选地，利用下述方法获得一维离散卷积运算的P个结果，具体包括以下步骤： 

3-1、利用下式定义序列h₁[l]，计算序列h₁[l]的离散傅里叶变换结果H[l']； 

$h_1\left[l\right]=W_a^{{(l-M)}^2/2},l=1,...,N_T$

其中，N_T在满足N_T≥(P+M-1)的同时，为2的最小整数次幂； 

3-2、针对矩阵C中的第n列，利用下式定义序列g_e[l]，利用分裂基FFT算法，计算序列g_e[l]的离散傅里叶变换结果G_e[l']； 

$g_e\left[l\right]=\left\{\begin{array}{cc}C(l,n),& l=1,...,M\\ 0,& l=M+1,...,N_T\end{array}\right.$

其中，C(l,n)为矩阵C中第l行第n列的元素； 

3-3、运算H[l']和G_e[l']的点积，并对点积做逆FFT变换，保留逆FFT变换中的第M～(P+M-1)个数据结果，即得到一维离散卷积运算中的P个结果。 

由于g_e[l]中仅有前M个元素为非0元素，因此，采用分裂基FFT的剪除算法，识别和去除不需要的运算操作，可以提高运算效率。 

（4）对应仰视角方向，对矩阵MD的每一行，进行一维离散卷积运算，保留一维离散卷积运算中的Q个结果，得到P×Q的矩阵ME。 

所述步骤（4）中的一维离散卷积运算的公式如下： 

$F_2(p,q)=\underset{n-1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_1(p,n)\&CenterDot;W_e^{\frac{{(q-n)}^2}{2}}$

其中，q的取值范围为1～Q。 

为了提高一维离散卷积的运算效率，优选地，利用下述方法获得一维离散卷积运算的Q个结果，具体包括以下步骤： 

4-1、利用下式定义序列h'₁[l]，计算序列h'₁[l]的离散傅里叶变换结果H'₁[l']； 

${h^{\&prime;}}_1\left[l\right]=W_e^{{(l-N)}^2/2},l=1,...,{N^{\&prime;}}_T$

其中，N’_T在满足N’_T≥(Q+N-1)的同时，为2的最小整数次幂； 

4-2、针对矩阵MD中的第p行，利用下式定义序列g'_e[l]，利用分裂基FFT算法，计算序列g'_e[l]的离散傅里叶变换结果G'_e[l']； 

${g^{\&prime;}}_e\left[l\right]=\left\{\begin{array}{cc}\mathrm{MD}(p,l),& l=1,...,N\\ 0,& l=N+1,...,{N^{\&prime;}}_T\end{array}\right.$

其中，MD(p,l)对应矩阵MD中第p行第l列的元素； 

4-3、运算H'₁[l']]和G'_e[l']的点积，并对点积做逆FFT变换，保留逆FFT变换中的第N～(Q+N-1)个数据结果，即得到一维离散卷积运算中的Q个结果。 

步骤3-2、步骤4-2中的FFT变换，以及步骤3-2、步骤4-2中对点积结果做的逆FFT变换（即IFFT），都通过分裂基FFT运算实现，且在运算过程中采用针对分裂基FFT的剪除算法，识别和去除不需要的运算操作，提高了运算效率。 

步骤3-1和步骤4-1都通过离线运算完成，事先存储在三维摄像声纳系统信号处理机的存储器内，步骤3-2、步骤3-3以及步骤4-2、步骤4-3在信号处理机中实时在线进行运算。 

所述FFT变换和逆FFT变换均通过分裂基FFT运算实现，且对于每个分裂基FFT运算，均根据具体的有效数据输入或输出特性，建立相应的剪除矩阵，所述剪除矩阵中的每个元素与分裂基FFT运算过程中的一个数据节点相对应，仅计算剪除矩阵中取值不为0的元素对应的数据节点。 

根据分裂基FFT运算的具体输入和输出分布特性，对分裂基FFT运算流程的所有数据节点，建立精确的剪除矩阵D_pr（该剪除矩阵可以人为设定，也可以依据一定的方法计算得到），然后，在分裂基FFT的实时运算中，利用D_pr识别和去除不需要的运算操作。 

在剪除矩阵D_pr中，每个元素与分裂基FFT运算流程的一个数据节点相对应；取值不为0（例如取值为1）的元素对应的数据节点为有效节点，运算流程中需要进行该数据节点的运算；取值为0的元素对应的数据节点 为无效节点，对应的运算可以直接略去；在有效节点的运算过程中，还利用D_pr中的相关元素，对该有效节点的运算过程进行合理的简化。 

（5）构建矩阵MF，矩阵MF中的每个元素为： 

$\mathrm{MF}(p,q)=\mathrm{ME}(p,q)\&times;W_a^{-{(p-1)}^2/2}{W}_e^{-{(q-1)}^2/2}$

其中，ME(p,q)对应矩阵ME中第p行第q列的元素。 

（6）计算矩阵MF中每个元素的模，得到声纳获取的各个方向上的波束强度，即得到三维摄像声纳波束。 

矩阵MF中的每个元素的模表明了一个波束方向上的波束强度，因此，可以直接用于相控阵三维摄像声纳系统的图像显示与处理。 

由每个元素对应的方位角和仰视角，以及该元素对应的波束强度，可以得到某一方向上的波束的强度信息，将所有元素对应的方向与波束强度进行合成，即可得到三维摄像声纳波束。 

本发明基于CZT变换和剪除分裂基快速傅里叶变换（Pruned Split-radix Fast Fourier Transform）的三维摄像声纳波束形成方法，依次对方位角方向和仰视角方向进行一维离散卷积运算，并且一维离散卷积运算通过一种剪除的分裂基FFT方法实现，在FFT运算过程中识别和去除不必要的运算操作，有效降低了三维摄像声纳波束形成方法所需的运算量。 

附图说明

图1为本发明中的M×N二维换能器阵列及方位定义示意图； 

图2为本发明基于CZT和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法的流程图； 

图3为本发明基于CZT和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法中的一维离散卷积运算过程示意图； 

图4为本发明中的分裂基FFT运算流程示意图，其中分裂基FFT长度N_T为16； 

图5为本发明基于CZT和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法中获得的二维波束方向图的结果示意图。 

具体实施方式

下面结合附图，对本发明基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法做详细描述。 

本实施例中设定M=N=48，换能器的水平间距和垂直间距均相等，需要运算128×128个方向的波束，仰视角θ_e和方位角θ_a的覆盖范围均为-25°～25°。 

相控阵三维摄像声纳系统的M×N换能器阵列如图1所示，图中θ_a为方位角，即波束与yoz平面的夹角；θ_e为仰视角，波束与xoz平面的夹角。。换能器阵列位于xoy平面，u为波束方向角，定义θ_a和θ_e分别为u对应方位的方位角和仰视角。 

如图2所示，一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，包括以下步骤： 

（1）分别计算每个换能器接收到的声波中，与载波频率f_k对应的第k个DFT变换系数S_m，n(k)，得到M×N个S_m，n(k)。载波频率f_k预先设定。 

其中，m为换能器的水平索引号，m=1，2，...，M；n为换能器的垂直索引号，n=1，2，...，N；k=0，1，...，K-1，K为DFT变换的长度。 

（2）利用M×N个S_m，n(k)，构建矩阵C，矩阵C中的每个元素为： 

$C(m,n)=w_{m,n}{S}_{m,n}\left(k\right)A_a^{m-1}{A}_e^{n-1}\&times;W_a^{-{(m-1)}^2/2}{W}_e^{-{(n-1)}^2/2}$

其中， 

$\left\{\begin{array}{c}{A}_{\text{a}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ai}})\\ W_a=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_a\right)\\ \mathrm{\&Delta;}{s}_a=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{af}}-\sin {\mathrm{\&theta;}}_{\mathrm{ai}}}{P-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ai}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{af}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

$\left\{\begin{array}{c}{A}_{\text{e}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ei}})\\ W_e=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_e\right)\\ \mathrm{\&Delta;}{s}_e=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{ef}}-\sin {\mathrm{\&theta;}}_{\mathrm{ei}}}{Q-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ei}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{ef}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

其中，w_m,n表示每个换能器的权重；j为虚数单位；f_k为载波频率；d为换能器的水平间距（或垂直间距）；c为水下声速；θ_ai表示声纳系统波束扫描范围的起始方位角；θ_af表示声纳系统波束扫描范围的终止方位角；P为波束扫描范围内的方位角数目；θ_ei表示声纳系统波束扫描范围的起始仰视角；θ_ef表示声纳系统波束扫描范围的终止仰视角；Q为波束扫描范围内的仰视角数目。 

（3）对应方位角方向，对矩阵C的每一列，进行一维离散卷积运算，保留一维离散卷积运算中的P个结果，得到P×N的矩阵MD； 

其中，一维离散卷积运算，如图3所示，包括以下步骤： 

3-1、利用下式定义序列h₁[l]，计算序列h₁[l]的离散傅里叶变换结果H[l']； 

$h_1\left[l\right]=W_a^{{(l-M)}^2/2},l=1,...,N_T$

其中，N_T在满足N_T≥(P+M-1)的同时，为2的最小整数次幂；例如，P+M-1=48+128-1，N_T≥175，由于2⁷=128，2⁸=256，则N_T为256。 

3-2、针对矩阵C中的第n列，利用下式定义序列g_e[l]，利用分裂基FFT方法，计算序列g_e[l]的离散傅里叶变换结果G_e[l']； 

$g_e\left[l\right]=\left\{\begin{array}{cc}C(l,n),& l=1,...,M\\ 0,& l=M+1,...,N_T\end{array}\right.$

其中，C(l,n)对应矩阵C中第l行第n列的元素； 

3-3、运算H[l']和G_e[l']的点积，并对点积做逆FFT变换，保留逆FFT变换y[l]中的第M～(P+M-1)个数据结果，即得到一维离散卷积运算中的P个结果。 

（4）对应仰视角方向，对矩阵MD的每一行，进行一维离散卷积运算，保留一维离散卷积运算中的Q个结果，得到P×Q的矩阵ME； 

其中，一维离散卷积运算，包括以下步骤： 

4-1、利用下式定义序列h'₁[l]，计算序列h'₁[l]的离散傅里叶变换结果H'₁[l']； 

${h^{\&prime;}}_1\left[l\right]=W_e^{{(l-N)}^2/2},l=1,...,{N^{\&prime;}}_T$

其中，N’_T在满足N’_T≥(Q+N-1)的同时，为2的最小整数次幂； 

4-2、针对矩阵MD中的第p行，利用下式定义序列g'_e[l]，利用分裂基FFT方法，计算序列g'_e[l]的离散傅里叶变换结果G'_e[l']； 

${g^{\&prime;}}_e\left[l\right]=\left\{\begin{array}{cc}\mathrm{MD}(p,l),& l=1,...,N\\ 0,& l=N+1,...,{N^{\&prime;}}_T\end{array}\right.$

其中，C(l,n)对应矩阵C中第l行第n列的元素； 

4-3、运算H'₁[l']]和G'_e[l']的点积，并对点积做逆FFT变换，保留逆FFT变换中的第N～(Q+N-1)个数据结果，即得到一维离散卷积运算中的Q个结果。 

步骤3-1和步骤4-1都通过离线运算完成，事先存储在三维摄像声纳系统信号处理机的存储器内，步骤3-2、步骤3-3以及步骤4-2、步骤4-3在信号处理机中实时在线进行运算。 

步骤3-2、步骤3-3以及步骤4-2、步骤4-3中，分裂基FFT的剪除方法根据分裂基FFT运算的具体输入和输出分布特性，对分裂基FFT运算流程的数据节点，建立精确的剪除矩阵D_pr，然后在分裂基FFT的实时运算中，利用D_pr识别和去除不需要的运算操作。 

图4为一个仅有部分输入有效的16点分裂基FFT的流程示意图，其输入序列中只有x(1)和x(2)为非0运算。图4中的W=cos(2π/16)-j×sin(2π/16)，此时，为了获得最终的FFT结果，只需要执行图4中实线所示的运算操作，而虚线对应的操作可以直接略去。与图4对应的剪除矩阵D_pr为 

$D_{\mathrm{pr}}={\left[\begin{array}{cccccccccccccccc}0& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0& 0& 0& 1& 1& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 1& 1& 0& 0& 1& 1& 0& 0& 1& 1& 0\\ 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1& 1\end{array}\right]}^T$

其中，上标“T”表示矩阵的转置。D_pr中每个元素与分裂基FFT运算流程的一个数据节点相对应；D_pr中取值为1的元素对应的数据节点为有效节点，运算流程中需要进行该节点的运算；D_pr中取值为0的元素对应的 数据节点为无效节点，对应的运算可以直接略去；在有效节点的运算过程中，还利用D_pr中的相关元素，对该有效节点的运算过程，进行合理的简化。因此，利用D_pr，可以有效的识别和去除分裂基FFT运算过程中的不必要操作。例如，在图4所示的例子中，由于D_pr(1，3)=0，因而针对数据节点DN₁的运算可以直接略去；由于D_pr(3，3)=1，D_pr(3，2)=1，D_pr(7，2)=0，因而针对数据节点DN₂的复数相加运算，可以简化为简单的赋值操作。 

图5给出了当入射方位为θ_a=θ_e=10°时，利用本发明中的波束形成方法获得的二维波束方向图。 

（5）构建矩阵MF，矩阵MF中的每个元素为： 

其中，ME(p,q)为矩阵ME中第p行第q列的元素； 

（6）计算矩阵MF中每个元素的模，得到声纳获取的各个方向上的波束强度，即得到三维摄像声纳波束。 

Claims (6)

1.一种基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，其特征在于，包括以下步骤：

（1）分别计算每个换能器接收到的声波中，与载波频率f_k对应的第k个DFT变换系数S_m，n(k)，得到M×N个S_m，n(k)；

其中，m为换能器的水平索引号，m=1，2，...，M；n为换能器的垂直索引号，n=1，2，...，N；k=0，1，...，K-1，K为DFT变换的长度；

（2）利用M×N个S_m，n(k)，构建矩阵C，矩阵C中的每个元素为：

$C(k,m,n)=w_{m,n}{S}_{m,n}\left(k\right)A_a^{m-1}{A}_e^{n-1}\&times;W_a^{-{(m-1)}^2/2}{W}_e^{-{(n-1)}^2/2}$

其中，

$\left\{\begin{array}{c}{A}_{\text{a}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ai}})\\ W_a=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_a\right)\\ \mathrm{\&Delta;}{s}_a=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{af}}-\sin {\mathrm{\&theta;}}_{\mathrm{ai}}}{P-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ai}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{af}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

$\left\{\begin{array}{c}{A}_{\text{e}}=\exp (-j2\mathrm{\&pi;}{f}_k\frac{d}{c}\sin {\mathrm{\&theta;}}_{\mathrm{ei}})\\ W_e=\exp \left(j2\mathrm{\&pi;}{f}_k\frac{d}{c}\mathrm{\&Delta;}{s}_e\right)\\ \mathrm{\&Delta;}{s}_e=\frac{\sin {\mathrm{\&theta;}}_{\mathrm{ef}}-\sin {\mathrm{\&theta;}}_{\mathrm{ei}}}{Q-1}\end{array}\right.,\begin{array}{c}\left|{\mathrm{\&theta;}}_{\mathrm{ei}}\right|<\frac{\mathrm{\&pi;}}{2}\\ \left|{\mathrm{\&theta;}}_{\mathrm{ef}}\right|<\frac{\mathrm{\&pi;}}{2}\end{array}$

其中，w_m,n表示每个换能器的权重；j为虚数单位；f_k为载波频率；d为换能器的水平间距；c为水下声速；θ_ai表示声纳系统波束扫描范围的起始方位角；θ_af表示声纳系统波束扫描范围的终止方位角；P为波束扫描范围内的方位角数目；θ_ei表示声纳系统波束扫描范围的起始仰视角；θ_ef表示声纳系统波束扫描范围的终止仰视角；Q为波束扫描范围内的仰视角数目；

（3）对矩阵C的每一列，进行一维离散卷积运算，保留一维离散卷积运算中的P个结果，得到P×N的矩阵MD；

（4）对矩阵MD的每一行，进行一维离散卷积运算，保留一维离散卷积运算中的Q个结果，得到P×Q的矩阵ME；

（5）构建矩阵MF，矩阵MF中的每个元素为：

$\mathrm{MF}(p,q)=\mathrm{ME}(p,q)\&times;W_a^{-{(p-1)}^2/2}{W}_e^{-{(q-1)}^2/2}$

其中，ME(p,q)为矩阵ME中第p行第q列的元素；

（6）计算矩阵MF中每个元素的模，得到声纳获取的各个方向上的波束强度，即得到三维摄像声纳波束。

2.如权利要求1所述的基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，其特征在于，所述步骤（3）中一维离散卷积运算的公式如下：

$F_1(p,n)=\underset{m-1}{\overset{M}{\mathrm{\&Sigma;}}}C(m,n)\&CenterDot;W_a^{\frac{{(p-m)}^2}{2}}$

其中，p的取值范围为1～P。

3.如权利要求1所述的基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，其特征在于，所述步骤（3）中的一维离散卷积运算，包括以下步骤：

3-1、利用下式定义序列h₁[l]，计算序列h₁[l]的离散傅里叶变换结果H[l']；

$h_1\left[l\right]=W_a^{{(l-M)}^2/2},l=1,...,N_T$

其中，N_T在满足N_T≥(P+M-1)的同时，为2的最小整数次幂；

3-2、针对矩阵C中的第n列，利用下式定义序列g_e[l]，利用分裂基FFT方法，计算序列g_e[l]的离散傅里叶变换结果G_e[l']；

$g_e\left[l\right]=\left\{\begin{array}{cc}C(l,n),& l=1,...,M\\ 0,& l=M+1,...,N_T\end{array}\right.$

其中，C(l,n)对应矩阵C中第l行第n列的元素；

3-3、运算H[l']和G_e[l']的点积，并对点积做逆FFT变换，保留逆FFT变换中的第M～(P+M-1)个数据结果，即得到一维离散卷积运算中的P个结果。

4.如权利要求1所述的基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，其特征在于，所述步骤（4）中的一维离散卷积运算的公式如下：

$F_2(p,q)=\underset{n-1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_1(p,n)\&CenterDot;W_e^{\frac{{(q-n)}^2}{2}}$

其中，q的取值范围为1～Q。

5.如权利要求1所述的基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，其特征在于，所述步骤（4）中的一维离散卷积运算，包括以下步骤：

4-1、利用下式定义序列h'₁[l]，计算序列h'₁[l]的离散傅里叶变换结果H'₁[l']；

${h^{\&prime;}}_1\left[l\right]=W_e^{{(l-N)}^2/2},l=1,...,{N^{\&prime;}}_T$

其中，N’_T在满足N’_T≥(Q+N-1)的同时，为2的最小整数次幂；

4-2、针对矩阵MD中的第p行，利用下式定义序列g'_e[l]，利用分裂基FFT方法，计算序列g'_e[l]的离散傅里叶变换结果G'_e[l']；

${g^{\&prime;}}_e\left[l\right]=\left\{\begin{array}{cc}\mathrm{MD}(p,l),& l=1,...,N\\ 0,& l=N+1,...,{N^{\&prime;}}_T\end{array}\right.$

其中，MD(p,l)对应矩阵MD中第p行第l列的元素；

4-3、运算H'₁[l']]和G'_e[l']的点积，并对点积做逆FFT变换，保留逆FFT变换中的第N～(Q+N-1)个数据结果，即得到一维离散卷积运算中的Q个结果。

6.如权利要求3或5所述的基于CZT变换和剪除分裂基快速傅里叶变换的三维摄像声纳波束形成方法，其特征在于，所述FFT变换和逆FFT变换均通过分裂基FFT运算实现，且对于每个分裂基FFT运算，均根据具体的有效数据输入或输出特性，建立相应的剪除矩阵，所述剪除矩阵中的每个元素与分裂基FFT运算过程中的一个数据节点相对应，仅计算剪除矩阵中取值不为0的元素对应的数据节点。

Images