0 引　言

随着目标主动声呐特征在水下探测识别领域的广泛应用，人们对水下目标散射传递函数仿真速度的要求越来越高。目前，水下目标散射传递函数仿真方法主要有有限元/边界元^{[1 − 2]}、有限差分^[3]、T矩阵^[4]以及板块元等方法^{[5 − 6]}等。其中，板块元方法是一种适用于中高频段的近似方法，它将积分运算转化为代数运算，计算时间正比于板块数量一次方，而其他方法计算时间正比于板块节点数量的2～3次方。为了进一步提高计算速度，范军等^[7]提出了基于可视化图形声学计算方法（GRACO）的板块元加速算法，借助计算机图像学光照渲染算法解决面元之间遮挡快速判断问题。李威等^[8]进一步将GRACO和弹跳射线法相结合，修正了凹表面板块之间多次散射问题。上述研究中的GRACO方法均是在GPU的图形渲染功能的基础上对目标进行遮挡判断和散射声场的计算，计算上限仍然受限制于计算机CPU的主频，GPU（Graphic Processing Unit）在并行计算中的核心优势没有完全发挥出来。针对多个频点计算问题，贾波等^[9]将水下目标近场散射声场分解为随频率变化的快速起伏项和缓变项的乘积，再利用切比雪夫插值计算各频点缓变项，将散射传递函数计算速度提高了3～5倍。

为进一步提高声学计算中的计算效率，基于图像处理器（GPU）的并行计算方法逐渐体现出了其计算优势。余立志等^[10]利用GPU的通用计算技术，实现了可控响应功率算法在声源定位计算中的应用，相比基于CPU的计算方法，效率提高了约20倍。吕云飞等^[11]实现了基于GPU的FIR分数时延波束形成算法，相较于CPU计算有80多倍的加速。受此启发，本文在对水下目标散射声场计算理论和并行计算可行性分析的基础上，提出了一种面向散射传递函数多频点计算的GPU处理方法，为目标散射传递函数多频点快速计算提供了新思路。

1 理论方法

水下目标散射声场满足Helmholtz积分公式：

$ {p_s}(\overrightarrow r ) = \frac{1}{{4{\text π} }}\int_s {{p_s}(\overrightarrow {{r_s}} )\frac{{\partial G(\overrightarrow r ,\overrightarrow {{r_s}} )}}{{\partial n}} - jw\rho {v_s}G(\overrightarrow r ,\overrightarrow {{r_s}} ){\rm d}S} 。$

(1)

式中：$ {p_s} $、$ {v_s} $分别为表面散射场声压和表面散射场法向振速；$ G(\overrightarrow r ,\overrightarrow {{r_s}} ) $为格林函数；$ w $、$ \rho $分别为角频率和密度；$ {e^{jwt}} $为时间因子。中高频情况下目标表面可以近似满足刚性边界：

$ {v_s}(\overrightarrow {{r_s}} ) = - \frac{1}{{jw\rho }}\frac{{\partial {p_i}(\overrightarrow {{r_s}} )}}{{\partial n}} 。$

(2)

式中：$ {p_i} $为入射声场声压。求解目标表面声场声压$ {p_s} $是计算远场散射声场的前提。

根据Helmholtz表面积分公式：

$\begin{split} {p_s}(\overrightarrow r ) =& \frac{1}{2{\text π} }\mathop {\lim }\limits_{r \to {r_s}} \int_S {{p_s}(\overrightarrow {{r_s}} )\frac{\partial G(\overrightarrow r ,\overrightarrow {{r_s}} )}{\partial n}} -\\ & jw\rho {v_s}(\overrightarrow {r_s} )G(\overrightarrow r ,\overrightarrow {r_s} ) {\text d}S 。\end{split}$

(3)

利用式（3）能够求解表面散射声场的声压。随着分析频率逐渐升高，目标表面的离散网格也大量增加，计算效率也随之降低。

为了快速求解表面散射场声压，范军等^[6]提出并发展了以物理声学近似为理论基础的板块元方法。板块元方法的核心思想是为把目标表面划分为若干个板块（尺度小于声波波长三分之一），根据Kirchhoff近似方法计算每个位于“亮区”板块散射声场，求和得到目标散射传递函数。

将目标表面划分为板块元，在收-发合置情况下，散射波势函数可表示为：

$ \phi = - \frac{A}{{2{\text π} }}\int\limits_S {{e^{i2{k_0}r}}V(\alpha )} \left( {\frac{{i{k_0}r - 1}}{{{r^3}}}\cos \alpha } \right){\rm d}S 。$

(4)

式中，A为入射波的幅值；$ \alpha $为目标表面声波的入射角；r是表面上一点到源点的距离；$ V(\alpha ) $为目标表面反射系数；$ {k_0} $代表波数。

在远场条件下（$ {k_0}r \gg 1 $），每个板块元散射声场可以表示为：

$ {I_m} = V(\alpha ) \cdot w\sum\limits_{n - 1}^K { \displaystyle\frac{{{e^{2i{k_0}({x_n}u + {y_n}u)}}({p_{n - 1}} - {p_n})}}{{(2{k_0}u + 2{k_0}{p_{n - 1}})(2{k_0}u + 2{k_0}{p_n}v)}}}。$

(5)

式中：$ {p_n} = \displaystyle \frac{{{y_{n + 1}} - {y_n}}}{{{x_{n + 1}} - {x_n}}} $ ，$ {p_0} =\displaystyle \frac{{{y_1} - {y_3}}}{{{x_1} - {x_3}}} $ ；K为板块元顶点数。本文计算中取为三角形板块（见图1），$ \overrightarrow {{r_0}} = u\overrightarrow i + v\overrightarrow j + w\overrightarrow k $为板块参考点到接收点的单位矢量，$ ({x_n},{y_n}) $为多边形顶点的坐标。

总散射声场可以表示为：

$ I = \sum\limits_{i = 1}^G {\sum\limits_{j = 1}^M {\left[ {V({\alpha _{ij}}')w' \sum\limits_{n = 1}^3 { \displaystyle\frac{{{e^{2i{k_0}({x_n}'{u_{ij}}' + {y_n}'{v_{ij}}')}}({p_{n - 1}} - {p_n})}}{{({u_{ij}}' + {p}_{n - 1}{v_{ij}}')({u_{ij}}' + {p_n}{v_{ij}}')}}} } \right]} }，$

(6)

式中：$ {x_n}' $、$ {y_n}' $为第$ (i,j) $个板块变换到二维平面的坐标。

根据式（4）第m个板块的近场散射声场可表示为：

$ {\phi _{s,m}} = - \frac{1}{{2{\text π} }}\frac{{ik{R_1} - 1}}{{R_1^3}}{e^{ik{R_1}}} \times w \sum\limits_{m = 1}^3 {\frac{{{e^{ - i({x_m}u' + {y_m}v')}}({p_{m - 1}} - {p_m})}}{{(u' + {p_{m - 1}}v')(u' + {p_m}v')}}} ，$

(7)

式中，$ {R_1} $为板块等效中心到发射点的距离。

2 面向多频点散射传递函数计算的GPU处理方法

受限于中央处理器（CPU）的处理速度，在板块数与频点数乘积数量级较大的情况下往往需要大量的计算时间。取板块数为M，频点数为N，时间t(s)，图2显示了不同排水量水下航行器板块数随频率的变化规律，图3展示了计算时间随板块数与频点数乘积的变化规律。

可以看出，10吨级中小型航行器网格可以到10¹⁸@30kHz，单一频段散射声场传递函数计算时间需要10¹²s@GHz。受制造工艺以及散热问题困扰，CPU处理速度提升空间越来越小，为此提高散射传递函数计算速度需要探索新的解决方案。

2.1 CUDA平台

GPU作为显示芯片常用于处理多边形3D渲染、视频像素级转换等方面，在硬件结构上存在大量的晶体管来为并行线程提供计算支持。相较于CPU，GPU更适合处理大量的重复计算。在此基础上，NVIDIA公司研发了并行计算的计算统一设备架构（Compute Unified Device Architecture，CUDA），主要用于NVIDA开发的GPU异构系统。此外，CUDA架构还提供了丰富的库和工具，用于简化并行编程的开发过程。在CUDA的编程模型中，程序员将计算任务划分为多个线程块（thread block），每个线程块包含多个线程（thread），并发执行不同的任务，利用GPU的并行计算能力加速计算过程。

2.2 基于切比雪夫多项式的频间插值技术

由式（3）每个板块散射声场可以分为随频率变化的快速起伏项和缓变项，即

$ {\phi _{s,m}} = A\phi _{s,m}^{(1)}\phi _{s,m}^{(2)} ，$

(8)

其中，

$ A = - \frac{1}{{2{\text π} }}\frac{{ikR - 1}}{{R_1^3}}，$

(9)

$ \phi _{s,m}^{(1)} = {e^{i2k{R_1}}} ，$

(10)

$ \phi _{s,m}^{(2)} = w\sum\limits_{m = 1}^3 {\frac{{{e^{ - i({x_m}u' + {y_m}v')}}({p_{m - 1}} - {p_m})}}{{(u' + {p_{m - 1}}v')(u' + {p_m}v')}}} 。$

(11)

当$ k{R_1} \gg 1 $时，相位因子较小，频率不是引起振荡的主要原因。在$ \phi _{s,m}^{(1)} $中的$ 2k{R_1} $项的值一般较大，因此$ \phi _{s,m}^{(1)} $随频率变化的起伏较大，是散射声场中随频率变化的快速起伏项。$ \phi _{s,m}^{(2)} $中的相位项为$ 2k({x_m}u + {y_m}u) $，由于板块元方法的计算特点所致板块的尺寸划分较小，在高频情况下，$ 2k({x_m}u + {y_m}u) $的值也相对较小，从而$ \phi _{s,m}^{(2)} $随频率变化的起伏较小，是散射声场中随频率变化的缓变函数。

根据切比雪夫逼近理论，对于给定定义域$ k \in [{k_m},{k_n}] $的目标函数$ f(k) $，其函数值可由逼近公式表示为：

$ f\left(k\right)={\sum }_{l=1}^{n}{c}_{l}{T}_{l-1}\tilde{\left(k\right)}-\frac{{c}_{1}}{2}。$

(12)

式中：$ {T_l}\tilde{(k)} $为$ l $阶第1类切比雪夫多项式。

$ \tilde{k}=\frac{2k-({k}_{m}+{k}_{n})}{{k}_{m}-{k}_{n}}，$

(13)

$ {c}_{l}=\frac{2}{n}{\sum }_{i=1}^{n}f\left({k}_{i}\right)\mathrm{cos}\left[(l-1)\frac{i-0.5}{n}{\text π} \right] 。$

(14)

式中：$ {k_i} $为$ [{k_m},{k_n}] $中的切比雪夫节点。

散射声场缓变项$ \phi _{s,m}^{(2)} $的计算，只需计算切比雪夫节点处的值，其他频率点的值可由给定阶数的切比雪夫逼近得到。

2.3 CUDA并行化方法

基于GPU适用于大规模计算的特点，利用CUDA编程模型结合切比雪夫逼近的方法，对水下目标在不同频率下的目标强度进行加速计算。

基于CUDA的计算特点提出2种GPU部署方法进行比较，其一为利用CUBLAS库的矩阵运算方法：以板块为单位建立数据结构，在GPU端计算所有板块，根据式（14）得到每个板块的$ l $阶切比雪夫插值系数$ c_{l} $。在设备端分配每一个板块一个线程进行计算，最终分别得到散射声场在宽带范围内的$ m \times n $阶缓变项矩阵和$ n \times m $阶快速变化项矩阵。

图4中的快速项和缓变项矩阵分布如图5所示，每个元素分别对应不同频点下不同板块的散射散射声场。利用式（8）和CUDA中的矩阵计算库CUBLAS，可以将图4中的矩阵快速相乘计算，计算得到的矩阵取第1列即为所取频带范围内的散射声场。具体而言，可以按照以下步骤进行计算：

步骤1　将图5中的矩阵分别转换为CUDA中的矩阵形式，可以使用CUDA提供的矩阵数据类型定义复数矩阵。

步骤2　调用CUBLAS提供的矩阵相乘函数，使用cublasZgemm函数进行复数矩阵相乘，该函数在CUBLAS V2版本中可以被使用，能够实现高效的矩阵计算。

步骤3　计算得到的矩阵包含了所有频率范围内的散射声场，取结果矩阵的第1列即可得到所取频带范围内的散射声场，即

$ {A_{i1}} = \sum\nolimits_{k = 1}^n {\phi _{ik}^{(1)}\phi _{ki}^{(2)}} 。$

(15)

由于第1种部署方法只需要矩阵相乘后的第1列数值，因此提出第2种根据计算得到的快速项和缓变项矩阵进行向量相乘的部署方法。具体步骤如下：

步骤1　在图5的2个矩阵中，分别取矩阵A、B的第1列和第1行，利用CUBLAS库的向量运算函数，相乘后得到单一频点的散射声场传递函数，如图6所示

步骤2　创建核函数，分配每个线程计算不同频点的散射声场传递函数，最终得到整个频带内的传递函数。

如图7所示，对2种不同的部署方法进行比较后可得，矩阵计算方法相较于后一种向量部署方法的计算速度有80%的提升，因此在后续仿真中，采用第1种部署方式来比较GPU的优化性能。

3 实验仿真及结果 3.1 仿真参数及设备参数

本文使用的CPU的硬件配置为AMD4800H，主频为2.90 GHz；GPU的硬件配置为NVIDIA公司的GeeForce RTX 2060，显存容量为6 GB。仿真对象分别采用排水量为100～10000 kg的航行器。仿真时采用收发合置的方式，入射方向向量和潜航器水平轴心的夹角取45°。保持模型比例不变，通过划分不同的表面板块数量来对比CPU和GPU配置下的计算时间，得到并行算法的加速比。

3.2 仿真结果

对10～30 kHz范围内对不同体积的潜航器模型进行网格划分，得到利用板块元方法计算目标强度时的板块数。图8对比了步长400 Hz情况下GPU、CPU计算不同排水量的潜航器模型散射声场传递函数的运行时间。

表1展示了在相同频点数不同板块数量情况下，GPU-CUDA框架相对于CPU计算的加速比。根据表中的数据可以看出，采用CUDA技术后目标的散射声场传递函数计算速度得到了显著提升，且提升幅度随着网格增加而增大。在当前的硬件配置下，当板块数与频点数乘积达到10⁷时，计算加速比可超过80。

4 结　语

本文介绍了一种基于GPU（图形处理器）的CUDA技术，用于实现宽带目标回波的快速计算方法。通过和传统的基于CPU的计算方法相比较，本文验证了CUDA技术在宽带回波计算中的加速效果。算例仿真表明，随着板块数量和频点数量的增加，目标散射声场传递函数的计算时间也会快速增加。本文针对多频点的板块元计算在CUDA平台提出了基于矩阵和向量的2种部署方法，对在计算量较大的情况下，矩阵部署方法的加速效果更明显。仿真结果表明，基于GPU的并行计算方法可以有效提高多频点条件下目标强度的计算速度，从而缩短宽带回波仿真的计算时间而且矩阵部署方法的加速效果较向量部署方法更加明显。本文的研究仅限于单GPU条件下的性能效果测试，为了进一步提高并行算法的加速效果，下一步可以开展基于多GPU的并行加速方法研究，以拓展CUDA技术的应用范围，实现更高效、更快速的水下目标回波计算。