一种基于声波叠加的光声效应模拟与仿真方法与流程

1.本发明涉及一种光声效应模拟与仿真方法，属于光声效应、光声前向模型、仿真模型开发技术领域。


背景技术：

2.光声效应指的是物体在激光的激励下产生声信号的物理现象。基于光声效应，光声成像技术在近些年获得广泛的关注并取得飞速的发展。目前，光声成像技术已经在小动物成像、血氧饱和度定量、肿瘤良恶性诊断上取得了一定的成功。
3.光声成像的发展在一定程度上得益于光声仿真软件的开发与成熟。其中k
‑
wave是比较成熟的光声仿真软件，被广泛应用于与光声相关的各项工作中。k
‑
wave中最重要的功能在于其光声前向模型仿真，通过数值方法求解光声方程，从初始压力分布产生光声信号来模拟光声效应。但该数值方法的计算较为复杂，运算时间较长，特别是在采样率高、初始压力分布网格较大的情况下。k
‑
wave中运算速度慢的缺点大大缩小了该工具包的应用范围。


技术实现要素：

4.本发明要解决的技术问题是：光声仿真软件中用于求解光声方程的数值方法的计算较为复杂、运算时间较长。
5.为了解决上述技术问题，本发明的技术方案是提供了一种基于声波叠加的光声效应模拟与仿真方法，其特征在于，包括以下步骤：
6.步骤1、将光声仿真软件提供的初始压力分布网格分成n
×
m个像素点，共有j个处于不同位置的超声探头；
7.步骤2、设j＝1；
8.步骤3、设p＝1；
9.步骤4、通过对n
×
m个像素点的像素值的设置，使得初始压力分布网格中仅第p个像素点起作用；
10.步骤5、在设定好的采样率下，使用光声仿真软件获取第j个超声探头输出的超声信号，该超声信号即为标准信号s，并获得初始压力分布网格中第p个像素点的中心到第j个超声探头中心之间的距离d；
11.步骤6、根据第j个超声探头的位置，通过获得的标准信号s中第p个像素点到第j个超声探头的距离以及距离d计算第p个像素点产生光声信号的相位算子；
12.步骤7、根据标准信号s中第p个像素点的像素值和第p个像素点到第j个超声探头的距离计算第p个像素点产生光声信号的幅值算子；
13.步骤8、通过相位算子和幅值算子改变标准信号s的相位和幅值，获得第p个像素点在第j个超声探头产生的光声信号；
14.步骤9、p＝p+1，若p>n
×
m，则进入步骤10，否则返回步骤4；
15.步骤10、j＝j+1，若j>j，则进入步骤11，否则返回步骤3；
16.步骤11、获得j个处于不同位置的超声探头的整体的初始压力分布的光声信号。
17.优选地，步骤6中，设第j个超声探头的二维坐标为(x,y)，第p个像素点的二维坐标为(x
p
,y
p
)，则第p个像素点产生光声信号的相位算子的计算过程包括以下步骤：
18.计算标准信号s中第p个像素点到第j个超声探头的距离d
p
：
[0019][0020]
通过超声波在介质中的传播速度v和采样率f获得第p个像素点产生的光声信号的相位算子τ(d
p
,d)：
[0021][0022]
当相位算子τ(d
p
,d)>0时，进行相位右移；当相位算子τ(d
p
,d)<0时，进行相位左移。
[0023]
优选地，步骤7中，第p个像素点产生光声信号的幅值算子的计算过程包括以下步骤：
[0024]
使用负二次的幂指数函数来近似能量消散的过程，如下式(3)所示：
[0025][0026]
式(3)中，g(d
p
)表示能量消散带来的幅值衰减函数；k是可调整的超参数，在同样的扩散距离内，k越大则能量扩散越大；
[0027]
得到第p个像素点产生光声信号的幅值算子a(d
p
；p)，如下式(4)所示：
[0028][0029]
式(4)中，p表示标准信号s中第p个像素点的像素值。
[0030]
优选地，步骤8中，第p个像素点在第j个超声探头产生的光声信号表示为s
p
：
[0031]
s
p
＝a(d
p
；p)g〈s|τ(d
p
,d)>
ꢀꢀꢀꢀꢀꢀ
(5)
[0032]
式(5)中，g<s|τ(d
p
,d)>表示使用相位算子τ(d
p
,d)对标准信号s进行循环移位操作。
[0033]
本发明提出了将初始压力分布网格分成若干像素点，分别计算每个像素点单独产生的光声信号，接着根据波的叠加原理计算整体的光声信号。本发明提供的方法将求解复杂的光声方程转换为简单的循环叠加操作，大大降低了运行时间。在大部分情况下，本发明提供的方法运行时间均远小于k
‑
wave的计算时间。
附图说明
[0034]
图1为三个样品用k
‑
wave和我们的方法生成的光声信号及用延迟叠加算法重构的光声图像；
[0035]
图2为实验3的运行速度对比(单位：秒)。
具体实施方式
[0036]
下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人
员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。
[0037]
本发明提供了一种基于声波叠加的光声效应模拟与仿真方法，包括以下步骤：
[0038]
步骤1、将光声仿真软件提供的初始压力分布网格分成n
×
m个像素点，共有j个处于不同位置的超声探头；
[0039]
步骤2、设j＝1；
[0040]
步骤3、设p＝1；
[0041]
步骤4、通过对n
×
m个像素点的像素值的设置，使得初始压力分布网格中仅第p个像素点起作用；
[0042]
步骤5、在设定好的采样率下，使用光声仿真软件获取第j个超声探头输出的超声信号，该超声信号即为标准信号s，并获得初始压力分布网格中第p个像素点的中心到第j个超声探头中心之间的距离d；
[0043]
以k
‑
wave为例，在一个固定的网格中将第p个像素点的像素值设置成1，其他像素点的像素值设置成0，保证整个初始压力分布中只有第p个像素点起作用。在设定好的采样率下，在像素点的周围安置一个超声探头获取超声信号(标准信号s)；通过计算第p个像素点中心到第j个超声探头中心之间的距离获取距离信息d；
[0044]
步骤6、根据第j个超声探头的位置，通过获得的标准信号s中第p个像素点到第j个超声探头的距离以及距离d计算第p个像素点产生光声信号的相位算子；
[0045]
设第j个超声探头的二维坐标为(x,y)，第p个像素点的二维坐标为(x
p
,y
p
)，则第p个像素点产生光声信号的相位算子的计算过程包括以下步骤：
[0046]
计算标准信号s中第p个像素点到第j个超声探头的距离d
p
：
[0047][0048]
通过超声波在介质中的传播速度v和采样率f获得第p个像素点产生的光声信号的相位算子τ(d
p
,d)：
[0049][0050]
当相位算子τ(d
p
,d)>0时，进行相位右移；当相位算子τ(d
p
,d)<0时，进行相位左移；
[0051]
步骤7、根据标准信号s中第p个像素点的像素值和第p个像素点到第j个超声探头的距离计算第p个像素点产生光声信号的幅值算子；
[0052]
影响光声信号幅值有三个因素：(1)像素值本身的大小，像素值大小和信号幅值成正比关系。(2)信号的衰减，信号在媒介中往往以幂指数的形式衰减。(3)信号能量的扩散，平面波能量消散往往以幂指数形式扩散。因为本发明的仿真方法是基于理想环境的，所以本发明忽略信号在媒介中的衰减，只考虑像素值和能量的扩散，则有：
[0053]
第p个像素点产生光声信号的幅值算子的计算过程包括以下步骤：
[0054]
使用负二次的幂指数函数来近似能量消散的过程，如下式(3)所示：
[0055][0056]
式(3)中，g(d
p
)表示能量消散带来的幅值衰减函数；k是可调整的超参数，在同样的扩散距离内，k越大则能量扩散越大；
[0057]
得到第p个像素点产生光声信号的幅值算子a(d
p
；p)，如下式(4)所示：
[0058][0059]
式(4)中，p表示标准信号s中第p个像素点的像素值；
[0060]
步骤8、通过相位算子和幅值算子改变标准信号s的相位和幅值，获得第p个像素点在第j个超声探头产生的光声信号s
p
：
[0061]
s
p
＝a(d
p
；p)g<s|τ(d
p
,d)>
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0062]
式(5)中，g<s|τ(d
p
,d)>表示使用相位算子τ(d
p
,d)对标准信号s进行循环移位操作；
[0063]
步骤9、p＝p+1，若p>n
×
m，则进入步骤10，否则返回步骤4；
[0064]
步骤10、j＝j+1，若j>j，则进入步骤11，否则返回步骤3；
[0065]
步骤11、获得j个处于不同位置的超声探头的整体的初始压力分布的光声信号。
[0066]
为了说明本发明的有效性，使用三个图片作为光声的初始压力分布，分别是简单物体、简单血管和复杂血管。前两张图片是k
‑
wave中自带的示例图片。
[0067]
本发明首先展示k
‑
wave和我们方法生成的结果，如图1所示。
[0068]
图1的第一列展示了三个样品，第二列和第三列展示了k
‑
wave和本发明提供的方法产生的光声信号。第四列和第五列展示了生成的光声信号通过延迟叠加算法重构的光声图像。可以看出，本发明提供的方法几乎可以达到和k
‑
wave同样的效果。
[0069]
接着测试在不同条件下本发明提供的方法相对于k
‑
wave的速度优势。实验概括起来如表1所示：
[0070]
表1不同实验中的参数设置
[0071][0072]
实验1：改变初始压力分布图片的分辨率，从64到256，其运行速度在表2中展示。可以发现三个样品在不同分辨率条件下，本发明提供的方法都具有明显的优势。
[0073]
表2实验1的运行速度对比(秒)
[0074][0075]
实验2：改变采样率测试运行速度。采样率越大，运行时间越长，但是高采样率对k
‑
wave的影响远大于本发明提供的方法。本发明提供的方法仍具有优势，运行时间如表3所
示。
[0076]
表3实验2的运行速度对比(秒)
[0077][0078]
实验3：改变探头数目测试运行速度。探头数目并不会改变k
‑
wave的运行速度，但是由于需要重复步骤3
‑
9，所以探头数目会极大影响我们方法的运行时间。运行时间对比在图2中展示。对于样品1和样品2，即使探头数目高达512，本发明提供的方法仍具有优势。对于样品3，当探头数目超过384时，k
‑
wave比本发明的速度快。但是总体来说，在绝大部分情况下，本发明提供的方法仍有速度优势。
[0079]
实验4：改变探头到样品的距离测试运行速度。在k
‑
wave中，探头更远意味着需要一个更大的网格把探头包含在内，网格数变大会极大影响k
‑
wave的运行速度。相比而言，本发明提供的方法受探头距离影响不大。运行时间结果如表4所示。
[0080]
表4实验4的运行速度对比(秒)
[0081]