本发明属于光声图像重建技术领域，公开了一种多声速自适应光声层析图像重建方法，包括步骤一，在超声层析成像的设备系统中用超声换能器发射超声波照射生物组织；步骤二，超声换能器接收生物组织发出的超声波信号；步骤三，采用弯曲射线模型的迭代重建算法对生物组织的声速分布进行重建；步骤四，将生物组织划分为若干个声速不同的区域；步骤五，在光声层析成像设备系统中用激光器发射脉冲激光照射组织；步骤六，超声换能器接收组织中光吸收点发出的光声信号；步骤七，计算光声信号传播到超声换能器的过程中所涉及的多声速；步骤八，多声速自适应应用到延时求和算法中，重建光声层析图像。本发明有效提高了光声成像重建图像的质量。(56)对比文件谷怀民;杨思华;向良忠.光声成像及其在生物医学中的应用.生物化学与生物物理进展.2006,(05),全文.孙正;贾艺璇.生物声学成像中声速不均匀性解决方法的研究进展.声学技术.2018,(05),全文.朱毅;李文超;张丽君;袁杰.拟合延时补偿的光声重建.光电工程.2015,(04),全文.
1.一种多声速自适应光声层析图像重建方法，其特征在于，包括以下步骤：步骤一，超声换能器发射超声波照射生物组织；
所述超声换能器既接收超声波也发射超声波，所述超声换能器阵列为环形，在超声换能器阵列中均匀选择多个超声换能器逐次发射超声波；
步骤三，采用弯曲射线模型的迭代重建算法对生物组织的声速分布进行重建；
针对步骤三中重建的声速分布，遍历每个声源点，与该点的声速相差在设定范围内的其他声源点归类为同一个区域，最后，遍历每个区域，以该区域内所有声源点的平均声速作为该区域的声速，进而，将一个生物组织划分为若干个声速不同的区域；
步骤六，超声换能器接收组织中光吸收点发出的光声信号；
步骤七，计算光声信号传播到超声换能器的过程中所涉及的多声速；
光声信号从光吸收点向组织表面传播的过程中经过多个声速不同的区域，最终被超声换能器接收，在步骤四划出的若干个不同声速区域的基础上，组织的声速区域受成像区域的限制而被离散成网格，设光声信号沿直线传播到超声换能器，传播向量记为 ，其中 表示光吸收点，表示超声换能器，同时 是传播路径的起点，是传播路径的终点，在每个区域里计算光声信号的传播路径时，选择与传播向量 最接近一致的那条路径，计算出每条传播路径以及与传播路径相对应的声速，最终得到整个传播路径上的多声速信息；
步骤八，多声速自适应应用到延时求和算法中，重建光声层析图像；
将多声速自适应应用到光声信号从吸收点到超声换能器的传播时间计算中，对每个吸收点来说，延时求和算法公式如下所示：，
其中， 表示基于延时求和算法重建的光声层析图像， 表示超声换能器 在时间 接收到的光声信号，表示光声信号从吸收点传播到超声换能器 所需的传播时间，是权重，多声速自适应应用到传播时间 的计算中，计算公式如下所示：，
其中 是所要求的传播时间， 为吸收点位置，为超声换能器 ，在声速 的区域上传播路径为 ， 是该区域的边界点，向量 最接近传播向量 ，即 和两向量夹角最小；同理，在声速为 的区域上传播路径为 ， 是该区域的边界点，同样地向量 最接近传播向量 ；直到传播路径的终点 所在的区域，该区域声速为 ，传播路径为 ；权重公式为：，
其中， 是吸收点到超声换能器的距离，是吸收点和超声换能器之间的夹角，在二维X‑Y平面内，以环形超声换能器阵列的圆心为坐标原点，假设吸收点 坐标为，超声换能器 坐标为 ，那么 计算公式为：，
2.根据权利要求1所述的多声速自适应光声层析图像重建方法，其特征在于，所述步骤二中超声换能器逐次接收步骤一中被测生物组织发出的超声波，超声波经过超声探头中的声波透镜到达超声换能器，超声换能器记录声压压强。
3.根据权利要求1所述的多声速自适应光声层析图像重建方法，其特征在于，所述步骤三中的弯曲射线模型的迭代重建算法如下：利用弯曲射线表征环形超声换能器阵列中的一个阵元发射超声信号，超声波在经过不同声速区域时发生折射，最终被超声换能器接收这一折射过程，公式如下：，
其中， 为超声信号到达超声换能器的第一到达时间向量；为被测生物组织的慢度向量； 为联系两者的系数矩阵，与超声换能器的空间分布以及组织的慢度分布 直接相关，声速重建的过程中，不断比较基于弯曲射线模型的理论声波到达时间和实验观测到达时间，直到两者偏差小于设定的阈值。
4.根据权利要求1所述的多声速自适应光声层析图像重建方法，其特征在于，所述步骤五中受照射的生物组织产生光声效应，满足方程：， (1)
其中 为生物组织热膨胀所产生的声压， 为比热容， 为等压体积热膨胀系数，是经验声速， 为由入射脉冲激光所激发的热源函数， ， 是生物组织的空间光吸收函数， 为光照射函数，激光打到吸收点和超声换能器开始接收信号为同时进行。
5.根据权利要求4所述的多声速自适应光声层析图像重建方法，其特征在于，脉冲激光照射被测生物组织下方或侧方的吸收点，超声换能器阵列所在平面与激光所在平面垂直。
6.根据权利要求4所述的多声速自适应光声层析图像重建方法，其特征在于，所述步骤六所采用的超声换能器与步骤二的超声换能器一致，所述超声换能器接收步骤五中光声效应产生的超声波经过超声探头中的声波透镜到达超声换能器，记录声压强度，设超声换能器数目为 ，记录第 个超声换能器由公式(1)中光声效应产生的声压信号为 , 。
一种多声速自适应光声层析图像重建方法
技术领域
[0001] 本发明属于光声图像重建技术领域，尤其涉及一种多声速自适应光声层析图像重建系统及方法。
背景技术
[0002] 光声层析成像是近年来新兴发展起来的一种无损医学成像方法，它结合了纯光学成像的高对比度特性和纯超声成像的高穿透深度特性，可以提供高分辨率和高对比度的活体组织成像，为研究生物组织的结构形态、生理特征、代谢功能、病例特征等提供了重要手段，在生物医学临床诊断以及在体组织结构和功能成像领域具有广泛的应用前景。光声层析成像的主要原理是光能到声能的转换，即光声效应。当一束脉冲激光照射到生物组织时，组织内部的光吸收体吸收激光能量，发生瞬时热弹性膨胀而激发出超声波，即光声信号。光声信号向组织表面传播，最终被超声换能器接收。根据接收到的光声信号的延时、幅度以及预设的系统声速，经计算机处理后重建出生物组织内部的光学特性分布，特别是光吸收系数分布。这些信息项可以被用于定量测量生物组织中的特定物质，比如血液里的氧合血红蛋白和脱氧血红蛋白，并在此基础上开展如乳腺癌的临床研究等。
[0003] 考虑到生物组织的复杂性，在光声层析成像时，为了简化问题，现有的方法通常假设光声信号在生物组织中的传播速度是常数，忽略了组织中的声速分布不均匀所带来的影响。研究表明，具有不同成分的生物组织通常具有不同的声学特性，其中软组织内部的声速差异可以达到10%。因此，光声层析图像重建算法中预设声速的不准确是造成重建图像质量不高的主要原因之一。
发明内容
[0004] 为解决上述技术问题，本发明的一种多声速自适应光声层析图像重建方法的具体技术方案如下：
[0005] 一种多声速自适应光声层析图像重建方法，包括以下步骤：
[0006] 步骤一，超声换能器发射超声波照射生物组织；
[0007] 步骤二，超声换能器接收生物组织发出的超声波信号；
[0008] 步骤三，采用弯曲射线模型的迭代重建算法对生物组织的声速分布进行重建；
[0009] 步骤四，将生物组织划分为若干个声速不同的区域；
[0010] 步骤五，激光器发射脉冲激光照射组织；
[0011] 步骤六，超声换能器接收组织中光吸收点发出的光声信号；
[0012] 步骤七，计算光声信号传播到超声换能器的过程中所涉及的多声速；
[0013] 步骤八，多声速自适应应用到延时求和算法中，重建光声层析图像。
[0014] 进一步的，步骤一中所述超声换能器既接收超声波也发射超声波，所述超声换能器阵列为环形，在超声换能器阵列中均匀选择多个超声换能器逐次发射超声波。
[0015] 进一步的，所述步骤二中超声换能器逐次接收步骤一中被测生物组织发出的超声波，超声波经过超声探头中的声波透镜到达超声换能器，超声换能器记录声压压强。
[0016] 进一步的，所述步骤三中的弯曲射线模型的迭代重建算法如下：
[0017] 利用弯曲射线表征环形超声换能器阵列中的一个阵元发射超声信号，超声波在经过不同声速区域时发生折射，最终被超声换能器接收这一折射过程，公式如下：
[0018] ，
[0019] 其中， 为超声信号到达超声换能器的第一到达时间向量；为被测生物组织的慢度向量；为联系两者的系数矩阵，与超声换能器的空间分布以及组织的慢度分布 直接相关，声速重建的过程中，不断比较基于弯曲射线模型的理论声波到达时间和实验观测到达时间，直到两者偏差小于设定的阈值。
[0020] 进一步的，所述步骤四包括如下具体步骤：
[0021] 针对步骤三中重建的声速分布，遍历每个声源点，与该点的声速相差在设定范围内的其他声源点归类为同一个区域，最后，遍历每个区域，以该区域内所有声源点的平均声速作为该区域的声速，进而，将一个生物组织划分为若干个声速不同的区域。
[0022] 进一步的，所述步骤五中受照射的生物组织产生光声效应，满足方程：
[0023] ， (1)
[0024] 其中 为生物组织热膨胀所产生的声压， 为比热容， 为等压体积热膨胀系数，是经验声速， 为由入射脉冲激光所激发的热源函数， ， 是生物组织的空间光吸收函数， 为光照射函数，激光打到吸收点和超声换能器开始接收信号为同时进行。
[0025] 进一步的，脉冲激光照射被测生物组织下方或侧方的吸收点，超声换能器阵列所在平面与激光所在平面垂直。
[0026] 进一步的，所述步骤六所采用的超声换能器与步骤二的超声换能器一致，所述超声换能器接收步骤五中光声效应产生的超声波经过超声探头中的声波透镜到达超声换能器，记录声压强度，设超声换能器数目为 ，记录第 个超声换能器由公式(1)中光声效应产生的声压信号为 , 。
[0027] 进一步的，所述步骤七中的光声信号从光吸收点向组织表面传播的过程中经过多个声速不同的区域，最终被超声换能器接收，在步骤四划出的若干个不同声速区域的基础上，组织的声速区域受成像区域的限制而被离散成网格，设光声信号沿直线传播到超声换能器，传播向量记为 ，其中 表示光吸收点，表示超声换能器，同时 是传播路径的起点，是传播路径的终点，在每个区域里计算光声信号的传播路径时，选择与传播向量 最接近一致的那条路径，计算出每条传播路径以及与传播路径相对应的声速，最终得到整个传播路径上的多声速信息。
[0028] 进一步的，所述步骤八将多声速自适应应用到光声信号从吸收点到超声换能器的传播时间计算中，对每个吸收点来说，延时求和算法公式如下所示：
[0029] ，
[0030] 其中， 表示基于延时求和算法重建的光声层析图像， 表示超声换能器在时间 接收到的光声信号，表示光声信号从吸收点传播到超声换能器 所需的传播时间，是权重，
[0031] 多声速自适应应用到传播时间 的计算中，计算公式如下所示：
[0032] ，
[0033] 其中 是所要求的传播时间， 为吸收点位置， 为超声换能器 ，在声速 的区域上传播路径为 ， 是该区域的边界点，向量 最接近传播向量 ，即和 两向量夹角最小；同理，在声速为 的区域上传播路径为 ， 是该区
域的边界点，同样地向量 最接近传播向量 ；直到传播路径的终点 所在的区域，该区域声速为 ，传播路径为 ；权重公式为：
[0034] ，
[0035] 其中， 是吸收点到超声换能器的距离，是吸收点和超声换能器之间的夹角，在二维X‑Y平面内，以环形超声换能器阵列的圆心为坐标原点，假设吸收点 坐标为，超声换能器 坐标为 ，那么 计算公式为：
[0036] ，
[0037] 计算公式为：
[0038] 。
[0039] 本发明的一种多声速自适应光声层析图像重建方法具有以下优点：
[0040] 本发明利用超声声速层析成像重建出生物组织的声速分布，划分出不同的声速区域，在此基础上计算出光声信号在传播过程中所涉及的多声速，把多声速自适应应用到光声层析图像重建中，避免了预设声速的不准确对光声成像造成的不利影响，从而有效提高光声成像重建图像的质量。
附图说明
[0041] 图1为本发明方法的流程图；
[0042] 图2为本发明方法中基于弯曲射线模型的迭代重建算法流程图；
[0043] 图3为本发明方法中计算光声信号在不同声速区域中的传播路径示意图；
[0044] 图4为使用k‑wave进行仿真验证的实验样品图；
[0045] 图5(a)为利用预设声速的光声层析图像重建效果图；
[0046] 图5(b)为利用本发明的多声速自适应光声层析图像重建效果图。
具体实施方式
[0047] 下面结合附图对本发明作进一步的说明。
[0048] 考虑到光声层析成像系统和超声层析成像系统具有类似的声信号接收和处理程序，因此可以以较低的软硬件成本实现两系统的融合。所以，为了克服现有光声层析成像技术中预设声速不准确的问题，本发明基于超声声速层析成像提供一种不依赖预设声速的、多声速能够自适应于不同生物组织的光声层析图像重建算法，该算法首先基于超声声速层析成像重建出生物组织的声速分布，划分出不同的声速区域，在此基础上计算出光声信号到超声换能器（超声传感器）的传播途径上所涉及的多声速，最后把多声速自适应应用到光声层析图像重建中，从而重建出质量更高的图像。
[0049] 如图1所示，本发明公开了一种多声速自适应光声层析图像重建方法，包括以下步骤：
[0050] 步骤一，在超声层析成像的设备系统中用超声换能器发射超声波照射生物组织；
[0051] 所述超声层析成像的设备系统可以完成光声层析图像的重建也可以完成超声层析图像的重建，该系统使用的超声换能器既可以接收超声波也可以发射超声波。系统中的超声换能器阵列为环形，超声换能器阵列围绕被测生物组织，它们之间填充重水。在超声换能器阵列中均匀选择若干个超声换能器逐次发射超声波。
[0052] 步骤二，超声换能器接收生物组织发出的超声波信号；
[0053] 超声换能器逐次接收步骤一中被测生物组织发出的超声波，超声波经过超声探头中的声波透镜到达超声换能器，声压压强被记录。所用的超声换能器数量越多，重建出图像的质量越高。
[0054] 步骤三，采用弯曲射线模型的迭代重建算法对生物组织的声速分布进行重建；
[0055] 算法原理如图2所示，环形超声换能器阵列中的一个阵元发射超声信号，超声波在经过不同声速区域时会发生折射，最终被超声换能器接收，这一折射过程可以利用弯曲射线表征，如下公式所示：
[0056] ，
[0057] 其中， 为超声信号到达超声换能器的第一到达时间向量；为被测生物组织的慢度（声速的倒数）向量； 为联系两者的系数矩阵，与超声换能器的空间分布以及组织的慢度分布 直接相关，声速重建的过程中，不断比较基于弯曲射线模型的理论声波到达时间和实验观测到达时间，直到两者偏差小于设定的阈值。
[0058] 步骤四，将生物组织划分为若干个声速不同的区域；
[0059] 针对步骤三中重建的声速分布，遍历每个声源点，与该点的声速相差在一定范围内（如<=10m/s）的其他声源点归类为同一个区域。最后，遍历每个区域，以该区域内所有声源点的平均声速作为该区域的声速。进而，将一个生物组织划分为若干个声速不同的区域。
[0060] 步骤五，在超声层析成像的设备系统中用激光器发射脉冲激光照射组织；
[0061] 超声层析成像的设备系统中，受照射的生物组织产生光声效应，满足方程：
[0062] ， (1)
[0063] 其中 为生物组织热膨胀所产生的声压， 为比热容， 为等压体积热膨胀系数，是经验声速， 为由入射脉冲激光所激发的热源函数， ， 是生物组织的空间光吸收函数， 为光照射函数，激光打到吸收点和超声换能器开始接收信号为同时进行。
[0064] 超声层析成像的设备系统中，激光器发出的激光脉冲的工作波长为1064nm，最大2
功率约为2J/cm。脉冲激光照射被测生物组织下方或侧方的吸收点。超声换能器围绕被测生物组织，它们之间填充重水。超声换能器阵列所在平面一般与激光所在平面垂直。
[0065] 步骤六，超声换能器接收组织中光吸收点发出的光声信号；
[0066] 步骤六所采用的超声换能器与步骤二的超声换能器一致，同样围绕被测生物组织，超声换能器和生物组织之间填充重水。超声换能器接收步骤五中光声效应产生的超声波，超声波经过超声探头中的声波透镜到达超声换能器，声压强度被记录。所用的超声换能器数量越多，重建出的图像质量越高。超声换能器的排列方式为环形。假设超声换能器数目为 ，记录第 个超声换能器由公式(1)中光声效应产生的声压信号为 , 。
[0067] 步骤七，计算光声信号传播到超声换能器的过程中所涉及的多声速；
[0068] 步骤七中的光声信号从光吸收点向组织表面传播的过程中经过多个声速不同的区域，最终被超声换能器接收。在步骤四划出的若干个不同声速区域的基础上，组织的声速区域受成像区域的限制而被离散成一定大小的网格，假设光声信号沿直线传播到超声换能器，传播向量记为 ，其中 表示光吸收点，表示超声换能器，同时 是传播路径的起点，是传播路径的终点，在每个区域里计算光声信号的传播路径时，选择与传播向量 最接近一致的那条路径。随着每条传播路径被计算出来，与传播路径相对应的声速也被计算出来，最终得到整个传播路径上的多声速信息。
[0069] 光声信号在每个声速区域里的传播路径需要单独计算，如图3所示，光声信号从光吸收点 传播到超声换能器 ,传播向量为 ，在声速区域2中光声信号的传播路径被计算为 ,这是因为向量 与传播向量 的夹角最小，即最接近传播方向。同理，求出在声速区域3中的传播路径为 ，在区域1中的传播路径为 。最终，光声信号总的传播路径为 ，对应的多声速为 。
[0070] 步骤八，多声速自适应应用到延时求和算法中，重建光声层析图像；
[0071] 多声速自适应应用到光声信号从吸收点（声源点）到超声换能器的传播时间计算中。对每个吸收点来说，延时求和算法公式如下所示：
[0072] ，
[0073] 其中， 表示基于延时求和算法重建的光声层析图像， 表示超声换能器在时间 接收到的光声信号，表示光声信号从吸收点传播到超声换能器 所需的传播时间，是权重。
[0074] 多声速自适应应用到传播时间 的计算中，计算公式如下所示：
[0075] ，
[0076] 其中 是所要求的传播时间， 为吸收点位置， 为超声换能器 ，在声速 的区域上传播路径为 ， 是该区域的边界点，向量 最接近传播向量 ，即和 两向量夹角最小；同理，在声速为 的区域上传播路径为 ， 是该区域的边界点，同样地向量 最接近传播向量 ；直到传播路径的终点 所在的区域，该区域声速为 ，传播路径为 。考虑超声换能器接收到的信号强度和吸收点与超声换能器之间的夹角成余弦关系，同时信号强度与距离成反比，所以本发明中权重公式设计为：
[0077] ，
[0078] 其中， 是吸收点到超声换能器的距离，是吸收点和超声换能器之间的夹角，在二维X‑Y平面内，以环形超声换能器阵列的圆心为坐标原点，假设吸收点 坐标为，超声换能器 坐标为 ，那么 计算公式为：
[0079] ，
[0080] 通过解析几何中两向量夹角计算公式可知， 计算公式为：
[0081] 。
[0082] 下面结合仿真实施用例，对本发明做出进一步的说明。
[0083] 仿真实施用例如下所示：
[0084] 采用k‑wave仿真工具生成如图4所示的实验样品，埋在样品中的拼成血管图样的三根头发（图中以1、2、3标注出）构成声源点，也是光吸收点，并且样品分为三个声速不同的区域，区域1声速为1400m/s，区域2声速为1510m/s，区域3声速为1590m/s，三根头发完全包含在区域3中。采用k‑wave仿真工具生成半环形超声换能器阵列，包含512个超声换能器，均匀分布在半环上。
[0085] 同样地，采用k‑wave仿真工具模拟超声发射，超声换能器阵列接收超声信号并送入上位机进行样品的声速分布重建和声速区域的划分。
[0086] 同样地，使用k‑wave仿真工具模拟脉冲激光器照射样品后头发发出光声信号并向外传播，超声换能器阵列接收光声信号并送入上位机进行多声速自适应光声层析图像的重建，重建的图像如图5(a)所示。
[0087] 为对本发明提出的方法和预设声速的方法进行光声成像效果的比较，将预设声速设置为1500m/s（上述样品中三个声速区域的平均声速），同样地代入到延时求和算法中进行光声层析图像的重建，重建的图像如图5(b)所示。
[0088] 由图5(a)和图5(b)可知，相比于预设声速的方法，本方法获得的光声层析图像伪影、畸变明显减少，图中三根头发清晰可见，图像质量得到了明显提高。
[0089] 以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
