基于CUDA的波数域三维超声全矩阵成像方法

基于cuda的波数域三维超声全矩阵成像方法
技术领域
1.本发明属于工业超声无损检测技术领域，涉及一种基于cuda的波数域三维超声全矩阵成像的实现方法。


背景技术：

2.作为重要的无损检测技术，超声检测从一开始的无损探伤、无损检测的定性分析，已经转向定量无损评价的方向发展，这也对超声检测技术的定量检测特性提出了更高的要求。随着现代工业结构的复杂化以及工业设备运行条件的极端化，对工件与设备质量的要求也越来越高，同时也要求超声检测技术能够发现和检测出更加细微的缺陷，并且对检出缺陷的位置、形态和类型都有了更高的要求。在这种情况下，传统的常规超声已经无法满足要求。相控阵换能器是由多组小型的晶片阵元按照一定的几何方式排列而成，通常有线阵、面阵以及环形阵列等结构。相控阵超声检测技术凭借其波束的柔性合成与自适应控制能力，在检测效率、检测范围、声束可达性和灵敏度方面具有明显的技术优势，也使其在近二十年来得到了快速发展和广泛应用。
3.在相控阵超声无损检测中，全矩阵数据采集通过相控阵换能器依次激励单个换能器阵元发射超声波，然后用所有阵元接收回波信号并存储，充分利用了相控阵超声换能器中全部阵元发射接收的信号数据，包含了换能器阵元从不同角度入射到被测工件内的声波信息，拥有更加全面的被测工件内部信息，对于缺陷位置、尺寸和形态的描述更全面、更准确。尤其是当缺陷位置不可知时，相控阵超声全矩阵成像无需反复调整声束偏转角度与聚焦位置，只需按照原有的数据采集方式即可获得未知缺陷的全面信息。这一特性大大提高了相控阵超声成像的检测精度。并且，全矩阵成像技术不需要预先计算延时聚焦法则，也不需要根据不同的缺陷类型和缺陷位置调整波束合成方案，只需按照全矩阵采集方式对所有被测工件进行检测与信号收集，然后通过上位机软件进行全区域聚焦成像，即可一次性完成一定区域内所有缺陷的检测。
4.相比于传统的二维超声成像检测，三维成像以其对于空间缺陷更为全面准确的表征，更为广泛的成像区域获得越来越多的关注。二维图像往往只能获取试件的单个截面信息，不能获得试件缺陷整体的评估，而在石油管道、流水线检测和螺纹紧固件等检测方面，三维成像逐渐取代二维，成为不可或缺的无损检测手段。三维图像可以有效定位缺陷在试件中的空间位置，能够渲染出缺陷的形状和特征，便于定量无损评估试件的破坏部位和损坏程度。
5.全矩阵三维成像具备诸多便利，但是其更大的数据量和计算消耗给后处理能力提出了更高的要求，严重影响了成像的效率和精度。相较于传统的时域算法，波数域算法借助于快速傅里叶变换较低的计算复杂度，大幅提高全矩阵三维成像的计算效率，而采用较为精确的八线性插值，可以有效弥补插值所带来的相位误差，提高感兴趣区域内缺陷成像的精度。为了在保证成像精度的基础上进一步提高计算效率，满足实时在线检测的要求，充分利用快速傅里叶变换并行计算和灵活的移植能力，波数域算法可以进一步利用具备cuda平
台的gpu来加速运算过程，并且借助于vtk模块实现gpu的二维渲染，达成原始数据到图像显示与键鼠交互的完整流程。


技术实现要素：

6.本发明的目的在于提供一种基于cuda的波数域三维超声全矩阵成像方法，采用本发明的方法可以提高三维全矩阵成像的效率和精度，渲染二维最大密度投影，实现实时在线三维超声无损检测的目标。根据所成三维图像的大小和范围设置五维快速傅里叶变换的基本参数，根据面阵超声换能器的中心频率和带宽设置波数域的截断滤波，基于weyl’s identity获得原始数据坐标系和图像坐标系之间的波数转换关系，通过八线性插值将波数域原始数据映射到波数域图像，采用三维快速反傅里叶变换得到所需的三维图像。通过vtk模块渲染三维图像的最大密度投影，并通过交互窗口显示，三维图像的观察视角和距离可以实时改变，数据可以实时刷新，有利于实时三维超声无损定量评价。
7.本发明采用的技术方案是：
8.本发明首先提供了一种基于cuda的波数域三维超声全矩阵成像方法，其包括如下步骤：
9.步骤一：设置面阵超声换能器，以直接接触的方式与各向同性均匀介质的试件耦合，使用全矩阵模式采集原始三维全矩阵数据，确定所用试件的声速、面阵超声换能器的参数以及所成三维图像的大小范围，以面阵超声换能器左下角的阵元的中心为坐标原点，以面阵超声换能器的两边为x轴和y轴，以垂直换能器平面的声波传播方向为z轴，建立原始数据坐标系；
10.步骤二：确定五维傅里叶变换每个维度的计算长度，基于weyl’s identity编写cuda核函数，计算原始数据坐标系和图像坐标系之间的波数转换关系；
11.步骤三：拷贝cpu内存中的原始三维全矩阵数据至gpu显存中，于cuda核函数中将三维全矩阵数据重排列为五维全矩阵数据，数据索引按维度组成；使用行-列算法分解五维傅里叶变换，按照矩阵转置的坐标关系改变数据在存储空间中的索引，将五维傅里叶变换转换为计算多次二维和三维傅里叶变换；使用快速五维傅里叶变换，将原始数据转换到波数域，并且对频率维度做矩形窗截断；
12.步骤四：对于波数域五维矩阵中发射维度的每组三维数据，按照步骤二中原始数据坐标系和图像坐标系之间的波数转换关系为插值索引，选取空间维度上相邻的若干个点为参照，使用插值算法分别计算波数域五维矩阵在三维图像下的stolt映射；
13.步骤五：对于步骤四中计算得出的波数域stolt映射，提取其中每组三维数据分别对齐叠加，获得三维波数域图像矩阵，做三维快速反傅里叶变换，将图像矩阵由波数域重新变换到空间域，生成三维超声全矩阵图像，截断感兴趣区域以外的部分；
14.步骤六：对于步骤五所生成三维超声全矩阵图像，将其从gpu显存中拷贝到cpu内存中，压入vtk模块的数据处理流水线，使用最大密度投影渲染为二维图像并显示在窗口，并且能够实时更新数据，交互改变视角位置和焦点远近。
15.作为本发明的优选实施方案，步骤一中所述的试件为均匀各向同性介质，其声速设为c，使用的面阵超声换能器的阵元数分布为其中离散下标w1,w2分别代表面阵超声换能器在x轴和y轴方向上的阵元数量，而x轴和y轴的阵元中心间距pitch分别为
试件及成像区域摆放在z轴正半轴，成像区域的像素点数在x，y和z轴方向分别设为n
x
,ny,nz，其像素点的间距分别为d
x
,dy,dz；成像区域的原点与原始数据坐标系原点重合；获得的原始三维全矩阵超声数据记为e
raw
(t0,v,u)。
16.作为本发明的优选实施方案，所述的步骤二具体为：
17.由深度、宽度和像素间隔求出五维傅里叶变换每个维度的所需要的点数值，分别为n
t
,其中离散下标t代表时间，v1,v2代表接收阵元在x和y轴方向的分量，u1,u2代表发射阵元在x和y轴方向的分量；使用全矩阵采集进行成像时，有
18.通过weyl’s identity和包含格林函数的前向传播模型，所述的weyl’s identity在三维空间中的具体表达式为：
[0019][0020]
其中，g(ω,x,y,z)为三维空间中平面波分解的格林函数，k为空间波数，k
x
,ky分别为空间波数k在x和y轴上的分量，ω为角频率，x,y,z分别为x，y和z轴方向上的空间位移；根据成像系统简化的前向传播模型，发射-接收信号对，也即接收数据即为e(t,v1,v2,u1,u2)，
[0021]
其关于时间的频谱e(ω,v1,v2,u1,u2)可以用如下关系式计算：
[0022][0023]
其中，函数g为格林函数，f(x,y,z)为前向传播模型中的缺陷散射点分布，也即所需要获得的三维图像。将weyl’s identity代入前向传播模型公式，以v1,v2,u1,u2为变量改写关系式：
[0024][0025][0026]
其中，k
v1
,分别为接收阵元在x和y轴方向的波数分量以及发射阵元在x和y轴方向的波数分量，f(k
x
,ky,kz)为三维图像f(x,y,z)的波数域形式。因此，e(t,v1,v2,u1,u2)的五维傅里叶变换形式可以表示为如下关系式：
[0027][0028]
其中又有如下替换关系式：
[0029][0030]
[0031][0032]
将上式改写为用k
x
,ky,kz来表示k,为了与原始数据的空间波数做区分，这里以kg,k
v1g
,k
v2g
代替：
[0033]kv1g
＝k
x-k
u1g
,
[0034]kv2g
＝k
y-k
v1g
,
[0035][0036]
其中k
x
,ky,kz分别为空间波数k在x、y、z轴上的分量，其中离散下标x,y,z分别代表x，y和z轴，k
v1
,分别为接收阵元在x和y轴方向的波数分量以及发射阵元在x和y轴方向的波数分量。
[0037]
对于给定的面阵探头参数和初始化成像区域，其空间波数和各个方向的波数分量已知，如下关系式表示：
[0038][0039][0040][0041][0042][0043][0044][0045][0046][0047]fs
为接收全矩阵信号的采样频率。根据步骤三做傅里叶变换之后根据探头带宽做矩形截断，其截断最小最大频率分别设为f
min
,f
max
，相对应的角频率范围ω∈[2πf
min
,2πf
max
]。因此截断之后的角频率为：
[0048][0049]
根据步骤三和步骤四，为了充分利用插值映射的原始数据，三维图像在z轴方向的空间波数kz设为如下关系式：
[0050][0051]
编写cuda核函数计算波数在原始数据坐标系k,和图像坐标系k
x
,ky,kz的转换关系，作为波数域成像计算的参照索引。
[0052]
作为本发明的优选实施方案，所述的步骤三具体为：三维原始数据e
raw
(t0,v,u)转换为五维数据e(t,v1,v2,u1,u2)，其中索引t0代表原始数据中的时间，索引v,u分别代表全矩阵采集的接收和发射阵元。表示为离散矩阵形式为e
raw
[t0][v][u]和e[t][v1][v2][u1][u2]，其数据索引的关系式如下：
[0053]
t＝t0[0054][0055][0056][0057][0058]
其中，函数floor代表向下取整的函数；使用行-列算法分解五维傅里叶变换，转换为分别执行多次三维和二维傅里叶变换，将傅里叶变换的结果及中间步骤记为en(t,v1,v2,u1,u2)其计算流程如下：
[0059]
(4)次二维点快速傅里叶变换：
[0060][0061]
(5)转置数据矩阵
[0062][0063]
(6)次三维n
t
,点快速傅里叶变换：
[0064][0065]
截断数据，保留ω∈[2πf
min
,2πf
max
]的部分，设其点数为
[0066]
作为本发明的优选实施方案，步骤四中所述的插值算法，也即stolt反变换，需要锁定步骤二中前向传播模型公式可以改写为：
[0067][0068]
由于缺陷成像检测更关注相对幅值，忽略正比系数。代入步骤二中预先计算的波数参量，由于空间波数在原始数据和图像数据的转换非线性，并且锁定k
u1
,k
u2
之后仅剩ω,k
v1
,k
v2
三个参量，本发明采用八线性插值为插值算法，通过相应cuda核函数完成计算。
[0069]
作为本发明的优选实施方案，所述的步骤五中波数域中沿所有发射阵元叠加，并做三维快速反傅里叶变换，其表达式如下所示：
[0070][0071][0072]
其中，f(z,x,y)即为所计算的三维超声图像，截断感兴趣区域之外的部分，得到感兴趣区域之内的图像i(z,x,y)。
[0073]
作为本发明的优选实施方案，步骤六中所述的压入vtk模块的数据处理流水线，使用最大密度投影渲染为二维图像，具体包括如下步骤：
[0074]
a：初始化vtk模块参数；
[0075]
b：三维图像数据做滤波操作；
[0076]
c：使用最大密度投影将三维图像数据渲染为二维图像；
[0077]
d：设置vtk交互窗口，在窗口中实时显示和刷新二维图像。
[0078]
本发明的有益效果主要表现在：
[0079]
(1)本发明利用快速傅里叶变换将常规的时域延时叠加转换为波数域处理，大幅降低以往时域算法的计算复杂度，从而提高三维超声成像的效率和精度，实现实时高帧率三维超声全矩阵成像。
[0080]
(2)本发明采用了高效简洁的波数域计算公式，能够整体计算图像中所有像素点的值，公式稳定性高，方法简便易于理解，避免常规时域算法分别延时叠加计算每个像素点的低效率，从而从原理上提高三维超声全矩阵成像的速度。
[0081]
(3)本发明采用了波数域的频谱处理技术，降低了滤波等提升图像质量的手段的计算时间，从而有效降低计算时间，以较小的计算代价获取较高的成像效果。
[0082]
(4)本发明采用了自主编写和修改的插值方法，解决难以在软件中实现复杂插值方法的问题，从而有效解决三维波数域超声成像的软件实现问题，并且其计算速度较快，同样降低了计算消耗。
[0083]
(5)本发明利用cuda的gpu并行处理能力进一步提高三维成像的效率，由于所提出的利用快速傅里叶变换转换为波数域处理、频谱处理技术和插值方法均在本发明中以cuda平台的方法实现，主要计算在gpu中进行，解决常规cpu方法中串行计算的低效率，大幅节省成像时间。
[0084]
(6)本发明的算法具备极强的便捷性和可移植能力，采用具备cuda平台的计算设备，较之普通fpga(现场可编程门阵列)避免烧写代码和代码更新升级周期相对较长的问题，便于纠错和个人定制，使得能够在复杂灵活的现场环境实现三维实时超声成像。
[0085]
(7)本发明的算法不仅仅提供三维超声图像，更能够在获得三维超声图像之后，进一步渲染出二维映射图像以便显示在屏幕，同时具备实时更新图像源和相机视角的能力，渲染图窗可以实时交互，这避免了常规方法只计算三维图像无法给用户提供直观实时的成像检测能力，大大有利于用户在现场及后处理中的读图与判别能力。
附图说明
[0086]
图1是本发明方法整体流程图。
[0087]
图2是三维检测模型的安装位置示意图。
[0088]
图3是使用vtk模块实时渲染显示二维图像流程图。
[0089]
图4是使用三维实时超声成像及实时渲染二维图像结果图。
具体实施方式
[0090]
以下结合附图对本发明作进一步说明。
[0091]
本发明的实施例如下：
[0092]
如图1所示，本发明利用一种基于cuda的波数域三维超声全矩阵成像方法的操作流程可以分为以下几个步骤：
[0093]
步骤一：设置面阵相控阵探头，以直接接触的方式与各向同性均匀介质的试件耦合，使用全矩阵模式采集原始数据。各部件及相关参数如图2所示，超声换能器为8
×
8的正方形面阵探头，每个阵元的中心间距在x和y轴方向的分量分别为d
x
＝1.5mm,dy＝1.5mm，每个阵元的形状为正方形，阵元与阵元间的缝隙在x轴和y轴方向分别为s＝0.1mm。试件选取为金属铝，其声速为c＝6350m/s，三维成像区域为面阵探头下方，其原点与面阵探头的原点重合，深度为hz，宽度在x和y轴方向上分别为h
x
和hy，每个方向的像素间隔分别为δ
x
,δy,δz，由此计算出每个方向的像素数n
x
＝ceil(h
x
/δ
x
),ny＝ceil(hy/δy),nz＝ceil(hz/δz)。
[0094]
步骤二：由深度、宽度和像素间隔求出五维傅里叶变换每个维度的所需要的点数值，分别为n
t
,通过weyl’s identity和包含格林函数的前向传播模型，推导出波数域中原始坐标系和图像坐标系的坐标轴波数转换关系：
[0095][0096][0097][0098]
具体实施需要改变变量的组成形式，将上式改写为用k
x
,ky,kz来表示k,为了与原始数据的空间波数做区分，这里以kg,k
v1g
,k
v2g
代替，其中下标的g,v1g,v2g代表在成像坐标系下的空间波数，x轴方向的波数和y轴方向的波数：
[0099]kv1g
＝k
x-k
u1g
,
[0100]kv2g
＝k
y-k
v1g
,
[0101][0102]
由于超声信号的带限特性，根据探头带宽做矩形截断，其截断最小最大频率分别设为f
min
,f
max
，相对应的角频率范围ω∈[2πf
min
,2πf
max
]。因此截断之后的角频率为：
[0103][0104]
步骤三：将原始全矩阵三维数据变换为五维矩阵，通过多次二维和三维傅里叶变换代替直接五维傅里叶变换，将数据转换至波数域。
[0105]
所述的三维原始数据e
raw
(t0,v,u)转换为五维数据e(t,v1,v2,u1,u2)，表示为离散矩阵形式为e
raw
[t0][v][u]和e[t][v1][v2][u1][u2]，其数据索引的关系式如下：
[0106]
t＝t0[0107][0108][0109][0110][0111]
使用行-列算法分解五维傅里叶变换，转换为分别执行多次三维和二维傅里叶变换，将傅里叶变换的结果及中间步骤记为en(t,v1,v2,u1,u2)其计算流程如下：
[0112]
(7)次二维点快速傅里叶变换：
[0113][0114]
(8)转置数据矩阵
[0115][0116]
(9)次三维n
t
,点快速傅里叶变换：
[0117][0118]
截断数据，保留ω∈[2πf
min
,2πf
max
]的部分，设其点数为
[0119]
步骤四：计算得到的波数域数据通过kg,k
v1g
,k
v2g
的转换关系，对于每一次发射所对应的三维数据矩阵，分别进行八线性插值计算，转换至对应的图像坐标系，x、y和z轴方向的插值索引关系式如下：
[0120]
(1)所有数据点依次并行读取，首先获取插值所需的立方体在三个方向x，y和z轴上坐标最小的点：
[0121][0122][0123]
indz＝floor((k
z-min(k))/(k(1)-k(0)))
[0124]
(2)由基准点分别计算立方体剩余7个点，截断超过的坐标索引为边界值。
[0125]
(3)根据锁定的分别代入，计算插值所用每个坐标轴方向的公共乘积系数borda,bordb,bordc,a,b,c：
[0126]
borda＝(1-ind
x
)/n
x
[0127]
bordb＝(1-indy)/ny[0128][0129]
a＝k
v1g-ind
x
[0130]
b＝k
v2g-indy[0131]
c＝k
g-indz[0132]
(4)锁定由公共系数和波数域数据计算八线性插值结果：
[0133][0134]
步骤五：用步骤四得到的每组发射对应的插值后的图像坐标系下的波数域数据，对齐叠加，和为一组三维数据矩阵，该矩阵做三维快速反傅里叶变换即得到三维超声全矩阵图像，根据每个方向的感兴趣区域大小n
x_roi
,n
y_roi
,n
z_roi
截断图像。
[0135]
步骤五中矩阵做三维快速反傅里叶变换，其表达式如下所示：
[0136][0137][0138]
其中，f(z,x,y)即为所计算的三维超声图像，截断感兴趣区域之外的部分，得到i
(z,x,y)。
[0139]
步骤六：将数据由gpu显存中拷贝至内存中，使用vtk模块的最大密度投影渲染二维图像，在数据映射过程中缩小网格间距，提升图像质量，设置相机视角和图像窗口，增加窗口交互功能。渲染过程中为了提高渲染效果，使用10倍升采样因子。本实施例使用vtk模块实时渲染显示二维图像的流程如图3所示，具体包括：
[0140]
a：初始化模块参数，设置vtk数据原点和网格间距，设置vtk体积映射的内存空间和采样间距，设置vtk体积的属性，包括颜色和不透明度等。
[0141]
b：三维图像数据中间处理，进行升降采样和滤波等操作，得到处理后的三维图像。
[0142]
c：处理后的三维图像数据使用最大密度投影渲染为二维图像数据，需要将体积映射和属性添加至vtk体绘制变量。并且设置vtk渲染变量及相机视角，将三维图像数据导入渲染变量，获得二维图像数据。
[0143]
d：设置vtk交互窗口变量，导入二维图像数据，根据原始全矩阵数据的变化，实时显示和刷新窗口图像，相机视角同样可以实时改变。
[0144]
设计一个斜向排列具有5个孔洞型缺陷的方形金属试件，5个孔洞大小一致，在倾斜方向上等距排列。以n
t
＝1024,为例,选择三维成像区域的像素数分别为n
x
＝64,ny＝64,nz＝512，五维傅里叶变换之后频率维度截断最小最大频率设为f
min
＝3mhz,f
max
＝7mhz。选取缺陷所处的大致范围为感兴趣区域，其大小为n
x_roi
＝64,n
y_roi
＝64,n
z_roi
＝300，通过波数域的计算处理，使用vtk渲染，成像结果如图4所示，可以完整有层次地利用像素点的强度在图像中展现设计的5个孔洞的缺陷，并且试件的方形底面也能完整的计算出来，图像质量较高，没有条纹感。
[0145]
由于快速傅里叶变换的低算法复杂度和gpu的并行计算能力，可以实现大范围实时高帧率三维超声全矩阵成像、渲染和显示。
[0146]
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明的保护范围应以所附权利要求为准。