基于球形粒子光脉冲激发效应的光声信号仿真方法、装置与流程

1.本发明涉及光声断层成像领域，尤其涉及一种基于球形粒子光脉冲激发效应的光声信号仿真方法、装置。


背景技术：

2.光声成像（photoacoustic imaging, 简称pai）是一种非入侵式和非电离式的新型生物医学成像方法。当脉冲激光照射到（热声成像则特指用无线电频率的脉冲激光进行照射）生物组织中时，组织的光吸收域将产生超声信号，我们称这种由光激发产生的超声信号为光声信号。生物组织产生的光声信号携带了组织的光吸收特征信息，通过探测光声信号能重建出组织中的光吸收分布图像。光声成像结合了纯光学组织成像中高选择特性和纯超声组织成像中深穿透特性的优点，可得到高分辨率和高对比度的组织图像，从原理上避开了光散射的影响，突破了高分辨率光学成像深度“软极限”，可实现深层活体内组织成像。
3.光声断层成像（photoacoustic computed tomography, 简称pact）是光声成像的一种成像模式，主要针对于较深层组织进行成像。是一种能够清晰获取物体二维断层图像或者三维立体图像的新兴成像技术。该种模式根据光声效应，利用非弹道激光使血流等特异性吸收光的物质产生声学信号，并通过超声换能器接收信号，反演重建出目标的结构和功能信息。在光声断层成像中，通常采用非聚焦大直径脉冲激光束实现组织表面的全场照明。光吸收体对入射光能的吸收导致被照射组织膨胀，随后快速产生宽带超声波。它们传播到组织表面，由机械扫描的非聚焦超声波接收器或接收器阵列检测到，经过放大处理后通过光声图像重建算法对光声信号的传播问题逆向求解。即通过超声换能器检测得到的原始光声信号，反推生物组织成像区域的光吸收分布，从而得到其光声图像。光声层析成像运用了生物组织结构固有的光学对比度，几乎适用于任何具有光吸收特性的生物组织成像。该技术结合了光学成像和超声成像的优点，能够实现分辨率和对比度都较高的生物组织成像。尤其是其非入侵式和非电离式的成像特点，在生物医学领域具有重要应用价值。
4.光声断层成像技术和设备的发展在一定程度上得益于光声信号仿真软件的研发和应用。目前，在多种光声信号仿真软件中，k-wave是开发较为完善、应用较为广泛的一款软件。k-wave软件是matlab和c++的开源声学工具箱，由bradley treeby和ben cox（伦敦大学学院）和jiri jaros（布尔诺理工大学）开发。该软件设计用于组织介质中的时域声学和超声仿真。仿真函数基于k空间伪光谱方法，核心是一个可以解释线性和非线性的波传播、非同质材料参数的任意分布和幂律声吸收的数学模型。最重要的功能在于其光声前向模型仿真，通过数值方法求解光声方程，从初始压力分布产生光声信号来模拟光声效应。但是，该软件采用的数值方法复杂，尤其在建立较大尺寸分析模型或者建立较高空间分辨率的模型时，运算时间会大大增加，甚至超出了计算机的运算能力，严重影响了仿真分析的效率，从而制约了该软件的应用场景。


技术实现要素：

5.针对现有技术不足，本发明的目的是提供一种基于球形粒子光脉冲激发效应的光声信号仿真方法，通过精简的数值计算方式大幅度提高运算效率，以解决k-wave软件在建立大尺寸或高空间分辨率分析模型时运算效率过低甚至因超过仿真平台能力不能计算的问题。
6.为达到上述目的，本发明采用的技术方案具体为：本发明实施例的第一方面提供了一种基于球形粒子光脉冲激发效应的光声信号仿真方法，所述方法包括以下步骤：s1，设定采样频率和采样时间，得到光声信号传播的时间矢量；根据待测体模所处的传输介质，设定光声信号的传输速度；s2，设置光声断层成像探测阵列，并计算光声断层成像探测阵列中各阵元的位置坐标；s3，根据待进行光声断层成像区域的重建分辨率，将待检测体模分解为若干球形粒子，并设定每个球形粒子的初始声压；s4，根据步骤s1得到的时间矢量、光声信号的传输速度和步骤s3设定的每个球形粒子的初始声压，计算球形粒子向四周辐射的光声信号的声压分布，并将球形粒子向四周辐射的光声信号的声压分布叠加转化为光声断层成像探测阵列中单个阵元接收到的各球形粒子构成的体模所产生的光声信号；s5，重复步骤s4，将光声断层成像探测阵列中每个阵元接收到的各球形粒子所产生的光声信号组成光声信号矩阵。
7.进一步地，所述步骤s1中设定采样频率和采样时间，得到光声信号传播的时间矢量包括：根据待测体模和需重建光声断层图像的尺寸、体模和探测阵元之间的最远距离以及需采样的深度在内的需求，设定采样频率和采样时间，从而得到光声信号传播的时间矢量。
8.进一步地，所述步骤s2中设置光声断层成像探测阵列包括：如果重建二维光声断层图像，则选择二维结构的光声断层成像探测阵列；如果重建三维光声断层图像，则选择三维结构的光声断层成像探测阵列。
9.进一步地，所述二维结构的光声断层成像探测阵列包括阵元均匀分布的环形阵列；所述三维结构的光声断层成像探测阵列包括碗状阵列。
10.进一步地，所述碗状阵列中各阵元通过斐波那契格点方法排布。
11.进一步地，所述步骤s3中根据待进行光声断层成像区域的重建分辨率，待检测体模分解为若干球形粒子，并设定每个球形粒子的初始声压包括：待进行光声断层成像区域的重建分辨率为，则将待测体模共均匀分解为m个半径均为的球形粒子,并设定各个粒子的初始声压,。
12.进一步地，所述步骤s4中计算光声断层成像探测阵列中单个阵元接收到的各球形粒子所产生的光声信号的公式如下：
其中，m为球形粒子的个数，为球形粒子的半径，为光声信号传播的时间矢量，为光声信号的传输速度，为阶跃函数，是阵元的位置坐标和球形粒子中心的距离，。
13.进一步地，所述步骤s5中将光声断层成像探测阵列中每个阵元接收到的各球形粒子所产生的光声信号组成光声信号矩阵，光声信号矩阵的公式如下：其中，为单个阵元接收到的各球形粒子所产生的光声信号，n为阵元个数。
14.本发明实施例的第二方面提供了一种基于球形粒子光脉冲激发效应的光声信号仿真装置，包括一个或多个处理器，用于实现上述的基于球形粒子光脉冲激发效应的光声信号仿真方法。
15.本发明实施例的第三方面提供了一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时，用于上述的基于球形粒子光脉冲激发效应的光声信号仿真方法。
16.与现有技术相比，本发明的有益效果为：本发明提供了一种基于球形粒子光脉冲激发效应的光声信号仿真方法，待检测体模分解为若干球形粒子，根据时间矢量、光声信号的传输速度和每个球形粒子的初始声压，计算光声断层成像探测阵列单个阵元接收到的光声信号，将每个阵元得到的光声信号进行拼接得到光声信号矩阵。本发明方法将待测体模受激光照射后产生的光声信号随时间的演化过程进行简化，直接考虑待测体模发出的所有光声信号经一定距离传播后在光声断层成像探测阵列单个阵元上的叠加效果。在不影响仿真信号的准确性的同时，有效降低了数值计算的复杂度，消除了光声信号仿真受到的来自光声断层成像区域的空间尺寸和分辨率的制约，可以快速得到待测体模发送到光声断层成像探测阵列各阵元的光声信号仿真数据，从而大大提高了仿真效率。
附图说明
17.为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
18.图1是本发明实施例中的光声信号数值计算仿真流程。
19.图2是本发明实施例中的光声断层成像探测阵列和体模示意图。
20.图3是本发明实施例中采用k-wave方法仿真光声信号的结果。
21.图4是本发明实施例中采用本发明方法仿真光声信号的结果。
22.图5是本发明实施例中采用k-wave方法仿真光声信号的xoy截面重建结果。
23.图6是本发明实施例中采用k-wave方法仿真光声信号的xoz截面重建结果。
24.图7是本发明实施例中采用k-wave方法仿真光声信号的yoz截面重建结果。
25.图8是本发明实施例中采用本发明方法仿真光声信号的xoy截面重建结果。
26.图9是本发明实施例中采用本发明方法仿真光声信号的xoz截面重建结果。
27.图10是本发明实施例中采用本发明方法仿真光声信号的yoz截面重建结果。
28.图11是本发明实施例提供的基于球形粒子光脉冲激发效应的光声信号仿真装置的结构框图。
具体实施方式
29.为使本技术的上述目的、特征和优点能够更加明显易懂，下面结合附图对本技术的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本技术。但是本技术能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本技术内涵的情况下做类似改进，因此本技术不受下面公开的具体实施的限制。
30.下面结合附图来具体描述本发明的优选实施例，其中，附图构成本技术一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。
31.本发明提出了一种基于球形粒子光脉冲激发效应的光声信号仿真方法，本发明方法的核心是将待检测体模分解为多个微小的球形粒子，根据光声成像相关理论，直接计算得出光声断层成像探测阵列各阵元接收到的上述粒子产生的光声信号。
32.如果要产生光声信号，通常由一个短脉冲的激光束照射吸收体，吸收体中的分子吸收光子后，当满足一定的条件时，吸收体分子的电子从低能级跃迁到高能级而处于激发态，而处于激发态的电子极不稳定，当电子从高能级向低能级跃迁时，会以光或热量的形式释放能量。在光声成像应用中通常会选择合适波长的激光作为激发源，使吸收的光子的能量转化为热能的效率最大，通常从光能转化为热能的效率可达到90%以上。释放的热量导致吸收体局部温度升高，温度升高后导致热膨胀进一步转换成压力升高，压力升高以超声波的形式传播，这各超声波就被称为光声信号。因此，光声信号的产生过程就是“光能”—“热能”—“机械能”的转化过程。
33.接下来，如图1所示，对光声断层成像探测阵列各阵元接收到的光声信号进行仿真计算主要包含如下步骤：s1，设定采样频率和采样时间，得到光声信号传播的时间矢量；根据待测体模所处的传输介质，设定光声信号的传输速度。
34.具体地，设定仿真所需的基本参数：根据待测体模所处的传输介质的特性，设定好在该种传输介质中光声信号的传输速度。并且根据数值计算的需求，如待测体模和需重建光声断层图像的尺寸，体模和探测阵元之间的最远距离以及需采样的深度等，设定采样频率和采样时间，从而得到光声信号传播的时间矢量。
35.s2，设置光声断层成像探测阵列，并计算光声断层成像探测阵列中各阵元的位置
坐标。
36.具体地，设计光声断层成像的探测阵列：根据仿真要求，如果重建二维光声断层图像，则选择二维结构的光声断层成像探测阵列，如阵元均匀分布的环形阵列；如果重建三维光声断层图像，则选择三维结构的光声断层成像探测阵列，如碗状阵列。此外，还可根据待测体模和需重建光声断层图像的尺寸，以及成像的分辨率等设置光声断层成像探测阵列，计算出光声断层成像阵列各阵元的位置坐标。
37.s3，根据待进行光声断层成像区域的重建分辨率，将待检测体模分解为若干球形粒子，并设定每个球形粒子的初始声压。
38.具体地，进行待测体模粒子分解的步骤包括：为了方便后续的数值计算，本发明依据波的叠加原理对待测体模受激光激发出的光声信号的物理描述进行简化。波的叠加原理是指，介质中同时存在几列波时，每列波能保持各自的传播规律而不互相干扰。在波的重叠区域里各点的振动的物理量等于各列波在该点引起的物理量的矢量和。由于要仿真的是各阵元接收到的待测体模传播光声信号的声压，该物理量为标量。因此，可以认为单个阵元处接收到的待测体模的光声信号的声压为构成待测体模的数个粒子传播到单个阵元处的光声信号的声压之和。如果待进行光声断层成像区域的重建分辨率为，则可将待测体模在空间坐标系中所占据的区域按照光声断层成像能够重建的最小尺度分解为m个半径均为的球形粒子,设定好各个粒子的初始声压,。
39.s4，根据步骤s1得到的时间矢量、光声信号的传输速度和步骤s3设定的每个球形粒子的初始声压，计算球形粒子向四周辐射的光声信号的声压分布，并将球形粒子向四周辐射的光声信号的声压分布叠加转化为光声断层成像探测阵列中单个阵元接收到的各球形粒子构成的体模所产生的光声信号。
40.本发明将球形粒子向四周辐射的光声信号的声压分布叠加转化为阵元接收到的各球形粒子所产生的光声信号，本发明方法将待测体模受激光照射后产生的光声信号随时间的演化过程进行精简，直接考虑待测体模发出的所有光声信号经一定距离传播后在光声断层成像探测阵列单个阵元上的叠加效果。
41.球形粒子向四周辐射的光声信号的声压分布的计算过程如下：当一个半径为的微小球形粒子被形式的光脉冲加热，球体内部就会产生初始声压。假设球体被无吸收的背景介质所包围，设表示球形粒子吸收体中心所在的位置，则球形粒子被形式的光脉冲照射激发并向粒子四周辐射的光声信号的声压分布可以用下列公示表示：其中，是阵元的位置坐标和球形粒子中心的距离，。是
光声信号的传播速度，是光声信号传播的时间。为阶跃函数，当时，；当时，。
42.将球形粒子向四周辐射的光声信号的声压分布叠加转化为光声断层成像探测阵列中单个阵元接收到的各球形粒子构成的体模所产生的光声信号的过程具体包括：叠加计算得到光声断层成像探测阵列第个阵元接收到的各球形粒子所产生的光声信号，公式如下：其中，m为球形粒子的个数，为球形粒子的半径，为光声信号传播的时间矢量，为光声信号的传输速度，为阶跃函数，是阵元的位置坐标和球形粒子中心的距离，。
43.s5，重复步骤s4，将光声断层成像探测阵列中每个阵元接收到的各球形粒子所产生的光声信号组成光声信号矩阵。
44.具体包括以下步骤：s501，设定阵元序号循环数。即先设定阵元序号循环数i=1。每重复步骤s4一次，则i=i=1。重复步骤s4，并判断循环次数，若i《n，n为阵元个数，则返回s4直至计算完所有阵元接收到的各球形粒子所产生的光声信号，否则进入s502。
45.s502，求解得到光声断层成像探测阵列光声信号矩阵：通过上述步骤，将光声断层成像探测阵列中每个阵元接收到的各球形粒子所产生的光声信号组成光声信号矩阵，最终得到光声断层成像探测阵列所有阵元接收到的光声信号矩阵：。
46.实施例1：本发明实例描述的是通过该方法仿真计算512元碗型阵列接收的来自阵列内部相互垂直的三根线段组成的体模的光声信号。受限于k-wave仿真对计算资源的要求，碗型阵列的半径设为20mm，包含阵列在内的区域设为40 mm *40 mm*40 mm，重建区域分辨率设为0.1mm。
47.使用本发明提供的方法按照如下步骤进行仿真：s1，设定仿真所需的基本参数。根据数值计算的需求，设定采样频率fs=50mhz和采
样时间t=50us，得到光声信号传播的时间矢量t=(0:1/fs:t)。并且，根据待测体模所处的传输介质的特性，设定光声信号的传输速度v=0.15cm/us。
48.s2，完成用于进行光声断层成像的探测阵列的设计。根据仿真要求，计算出光声断层成像阵列各阵元的位置坐标。该实施例中，如图2所示，仿真的是512元碗型超声探测阵列，要求各阵元均匀分布在底部具有规定半径圆孔的半球面上。该实施例中通过斐波那契格点方法构造出均匀分布在设定的碗状球面上的所有阵元。
49.s3，进行待测体模粒子分解。该实施例中仿真的体模为位于阵列内部分别平行于笛卡尔坐标系x、y、z坐标轴并且相互垂直的三个线段。假设待进行光声断层成像区域的重建分辨率为，则将三根线段构成的待测体模共均匀分解为m=603个半径均为的球形粒子,设定好各个粒子的初始声压,。本发明实施例1中，重建区域分辨率设为0.1mm，则球形粒子的半径为0.05mm。
50.s4，设定阵元序号循环数。为了通过循环方式获得各阵元的光声信号，先设定阵元序号循环数i=1。
51.s5，数值计算光声断层成像探测阵列单个阵元接收到的光声信号。通过数值计算得到光声断层成像探测阵列第个阵元接收到的各球形粒子所产生的光声信号：其中，。
52.s6，进行阵元信号计算循环判断。判断是否全部生成各阵元光声信号。首先i=i+1,若i《512，则返回s5，否则进入s7。
53.s7，求解得到光声断层成像探测阵列光声信号矩阵。通过上述步骤，最终可以得到光声断层成像探测阵列所有阵元接收到的光声信号矩阵：。
54.本发明首先对比上述方法和k-wave仿真的光声信号，如图3和图4所示。可以看出，本发明提供的方法可以得到和k-wave基本相同的仿真信号。两种方法产生的信号无论在信号波形还是幅度方面均十分接近。再对两种方法产生的仿真信号用同一种方法进行光声断层图像重建，如图5~图10所示，对两种方法产生的光声断层成像阵列接收的光声信号进行处理能够得到几乎相同的重建结果。因此，本发明提供的方法能够实现和k-wave同等质量的仿真计算效果。
55.进一步，对比两种方法对同一仿真条件的数值计算时间如表1所示：
表1：两种方法的数值计算时间比较表对同样规模的分析模型进行数值计算，本发明的仿真用时远远少于k-wave软件。可见本发明提供的方法的计算效率远高于k-wave软件。
56.进一步地，在进行上述仿真的计算机上继续开展数值计算，选择更大的仿真区域,或者选择更高的重建区域分辨率,则k-wave软件仿真会显示计算机内存不够的错误提示并且终止运算。而本发明提供的方法则可以继续正确进行。
57.综上所述，本发明所提出的基于球形粒子光脉冲激发效应的光声信号仿真方法，由于其仿真机理不受仿真区域尺寸和重建区域分辨率等因素的限制，能够在更大规模或更精细的仿真条件下具有较高的运算效率和运算能力，相较于k-wave软件等其他光声信号仿真方法具有明显优势。本发明提出的方法可以应用于较大规模复杂结构光声断层成像阵列的成像能力和成像性能的评估，可用于较大尺寸待测体模的光声断层图像重建结果分析，也可用于高分辨率光声断层图像重建情形下的数值分析。因此，本发明提出的方法在光声断层成像技术研究中具有较为广泛且实际的应用价值。
58.与前述基于球形粒子光脉冲激发效应的光声信号仿真方法的实施例相对应，本发明还提供了基于球形粒子光脉冲激发效应的光声信号装置的实施例。
59.参见图11，本发明实施例提供的一种基于球形粒子光脉冲激发效应的光声信号装置，包括一个或多个处理器，用于实现上述实施例中的基于球形粒子光脉冲激发效应的光声信号方法。
60.本发明基于球形粒子光脉冲激发效应的光声信号装置的实施例可以应用在任意具备数据处理能力的设备上，该任意具备数据处理能力的设备可以为诸如计算机等设备或装置。装置实施例可以通过软件实现，也可以通过硬件或者软硬件结合的方式实现。以软件实现为例，作为一个逻辑意义上的装置，是通过其所在任意具备数据处理能力的设备的处理器将非易失性存储器中对应的计算机程序指令读取到内存中运行形成的。从硬件层面而言，如图11所示，为本发明基于球形粒子光脉冲激发效应的光声信号装置所在任意具备数据处理能力的设备的一种硬件结构图，除了图11所示的处理器、内存、网络接口、以及非易失性存储器之外，实施例中装置所在的任意具备数据处理能力的设备通常根据该任意具备数据处理能力的设备的实际功能，还可以包括其他硬件，对此不再赘述。
61.上述装置中各个单元的功能和作用的实现过程具体详见上述方法中对应步骤的实现过程，在此不再赘述。
62.对于装置实施例而言，由于其基本对应于方法实施例，所以相关之处参见方法实施例的部分说明即可。以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本发明方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。
63.本发明实施例还提供一种计算机可读存储介质，其上存储有程序，该程序被处理器执行时，实现上述实施例中的基于球形粒子光脉冲激发效应的光声信号方法。
64.所述计算机可读存储介质可以是前述任一实施例所述的任意具备数据处理能力的设备的内部存储单元，例如硬盘或内存。所述计算机可读存储介质也可以是任意具备数据处理能力的设备，例如所述设备上配备的插接式硬盘、智能存储卡（smart media card，smc）、sd卡、闪存卡（flash card）等。进一步的，所述计算机可读存储介质还可以既包括任意具备数据处理能力的设备的内部存储单元也包括外部存储设备。所述计算机可读存储介质用于存储所述计算机程序以及所述任意具备数据处理能力的设备所需的其他程序和数据，还可以用于暂时地存储已经输出或者将要输出的数据。
65.本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。