CN103142216B - 一种基于光声成像技术的多层介质声速计算的方法 - Google Patents

Abstract

本发明公开了一种基于光声成像技术的多层介质声速计算的方法，包括以下步骤：根据传感器脉冲响应函数，对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据；设定组织中不同介质的初始声速及声速迭代范围；以逆滤波数据及所设定组织各层声速为基础，进行光声图像重建；从每次重建得到的光声图像中计算并提取图像中声源的分布，依据声源分布信息确定迭代是否完成，若完成则输出声速。本发明提供了一种简单无创的用于测量生物组织中声速的方法，计算简单，复杂度小，效果突出。

Description

一种基于光声成像技术的多层介质声速计算的方法

技术领域

本发明涉及一种基于光声成像技术的多层介质声速计算的方法，具体地说就是根据传感器脉冲响应函数，对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据，然后对各组织中可能的声速进行迭代，以逆滤波得数据和每组可能的声速为基础重建光声图像，对于每次得到的图像计算和提取其声源分布的信息，并以此判断该组声速是否为最优解，反馈控制声速的迭代，直到产生最优解为止。

背景技术

目前虽然声速的测量方法有很多，比如正弦连续波共振法、脉冲时延法、相位比较法，但它们往往需要精密的仪器设备，而且这些方法不适合生物组织的测量条件。同时生物组织中的介质并不唯一，甚至是多层的，而现有的方法大多只能测定出生物组织中的平均声速，因此对于生物组织(特别是活体组织)中各介质中声速的确定，尚无效果特别理想的方法。

基于光声效应的光声成像利用脉冲激光激发光声信号，并且通过超声探头检测光声信号，进而反演出组织内光吸收特性分布。光声成像利用生物组织内的光吸收特性分布来重构图像，因此具有光学成像的高对比度。此外，光声成像检测的是光声信号，所以对于深处组织成像，它又具有超声成像高空间分辨率的优点。由于以上原因，光声成像已经被广泛地应用于人体和动物组织成像中，例如小鼠脑部血红蛋白浓度的成像，动物关节成像，人体关节成像，人体乳腺成像。由于光声成像对生物组织无创，并且兼具了光学成像高对比度和超声成像高空间分辨率的优点，利用该技术对生物组织进行较为精准的成像，而后根据重建图像的反馈，可以有效地确定组织中多层介质的声速。

该方法在光声成像技术近乎完美的成像质量的基础上，利用计算机强大的计算能力确定组织中各介质的声速，实现了对生物组织的无创测量，且避免了精密仪器的使用及其所带来的繁复的操作。

发明内容

发明目的：本发明所要解决的技术问题是针对传统声速测量方法无法有效确定人体组织中各介质中声速的问题，提供一种基于光声成像技术的多层介质声速计算的方法。

为了解决上述技术问题，本发明公开了一种基于光声成像技术的多层介质声速计算的方法，包括以下步骤：

步骤一，根据传感器脉冲响应函数，对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据；

步骤二，设定组织中不同介质的初始声速及声速迭代范围；

步骤三，以逆滤波数据及所设定组织各层声速为基础，进行光声图像重建；

步骤四，从每次重建得到的光声图像中计算并提取图像中声源的分布；

步骤五，依据声源分布信息确定迭代是否完成，若完成则输出声速。

本发明中，优选地，所述逆滤波在时域中进行：先以单个N型波(N型波是指点源光声信号，得名于其近似于N形的时域波形)经过传感器系统函数响应后输出的时域波形为模板，对每个传感器阵元的输出信号进行互相关运算，以确定每个N型波到达传感器的时刻。之后对于各阵元的输出信号，从所确定的第一个N型波开始，确定其主瓣的幅值后，将其波形减去，接着用同样方法确定第二个N型波主瓣幅值，依次类推，以达到分隔出该路上单个N型波输出波形的目的。最后用单个N型波经过传感器系统函数响应后波形变化的先验知识，还原出每路信号在被传感器响应前的原始波形；

本发明中，优选地，所述根据反馈或初值设定声速的方法在初次迭代时设定所有组织声速为初值C，C取值在1300m/s至1600m/s之间；

本发明中，优选地，所述以声速为基础进行光声重建采用延迟求和法，用设定的组织中的介质声速以及介质的几何尺寸计算出声波传播时间，以此找出待重建的每个像素点在数据矩阵中所对应的一系列数据，并对这些数据进行加权求和；

本发明中，优选地，所述提取并计算重建图像中声源分布的方法是用一个变量P[I(m，n)]对每个像素点(m，n)进行处理，若其灰度值I(m，n)大于设定的阈值(可设为128)，则P[I(m，n)]置1，反之置0，最后计算I²(m，n)关于P[I(m，n)]的加权平均，即可得到图像中声源分布的信息，简称为声源集中强度。

本发明中，优选地，所述根据声源分布信息判定迭代完成与否的方法是用一个长度可变的数组存放声源集中强度值，每次产生出一个新的值后，将其添加至数组末端，若数组中的倒数第二个数是极大值，则停止迭代，输出最优的声速；否则，继续迭代。

附图说明

下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。

图1是本发明方法中N型波的时域波形。

图2是本发明方法中N型波经系统函数响应后的波形。

图3是本发明方法的一个具体实例所使用的组织结构图。

图4是本发明方法的流程图。

具体实施方式：

本发明公开了一种基于光声成像技术的多层介质声速计算的方法，包括以下步骤：

步骤一，根据传感器脉冲响应函数，对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据；

步骤二，设定组织中不同介质的初始声速及声速迭代范围；

步骤三，以逆滤波数据及所设定组织各层声速为基础，进行光声图像重建；

步骤四，从每次重建得到的光声图像中计算并提取图像中声源的分布；

步骤五，依据声源分布信息确定迭代是否完成，若完成则输出声速。

本发明中，步骤一，单个N型波经过传感器系统函数响应后的时域信号可由公式(1)计算得到。

$r\left(t\right)=\∫f\left(\mathrm{\τ}\right)\·h(t-\mathrm{\τ})\mathrm{d\τ}=\∫f\left(\mathrm{\τ}\right)(\frac{1}{2\mathrm{\π}}\∫H\left(\mathrm{\ω}\right)e^{\mathrm{j\ω}(t-\mathrm{\τ})}\mathrm{dt})\mathrm{d\τ}---\left(1\right)$

其中r(t)为单个N型波经过系统函数响应后的输出，f(t)为单个N型波的数学表达式，且本实施例中，取A＝40，B＝2，C＝1.5，D＝1，其时域波形如图1所示。系统函数ω为频率，BW为频带宽度，取3.5MHz。r(t)的波形如图2所示。

对比图1和图2，可以发现N型波经系统函数响应前后，其幅度峰值所在时刻不变，因此只要确定了某个N型波的输出响应的幅度峰值所在时刻，该N型波的幅度峰值时刻也随之确定。我们可以利用N型波的这个性质来定位单个N型波。同时由于传感器是一个线性系统，因此当输入N型波的幅值被放大一定倍数时，其经传感器系统函数的输出响应的幅值也被放大相同倍数。在定位了单个N型波的基础上，我们可以根据这个先验知识通过此N型波输出响应的峰值与上述r(t)的峰值之比来确定此放大倍数，进而还原出该N型波的全部波形。

接着以上述r(t)为模板，对每个传感器阵元输出的信号进行互相关运算，对第k个传感器阵元所接收信号进行互相关运算的结果可用公式(2)计算得到。

$R_k\left(\mathrm{\τ}\right)={\∫}_{-\∞}^{+\∞}r(t-\mathrm{\τ})p_k\left(t\right)\mathrm{dt}---\left(2\right)$

其中τ表示r(t)相对于p_k(t)的延迟时间，p_k(t)表示第k个传感器阵元所输出的信号。之后找到每路R_k(τ)中的极大值，并返回对应的τ值，记为这里的变量用下标k区分各传感器阵元，用上标表示同一传感器阵元中同一变量在不同时刻的值。

然后对于每路p_k(t)，认为为第一个N型波的响应到来的时刻，根据主瓣幅值定位其波形利用r(t)与f(t)波形关系的先验知识，还原出的波形，接着用p_k(t)减去得对于认为为第二个N型波的响应到来的时刻，根据主瓣幅值定位其波形还原出的波形，接着用减去得依次进行下去即还原出一个传感器阵元所接收到的信号。

本发明中，步骤二，对可能的声速进行设定具体可按如下方式进行。

对于每层介质的声速C_mn，其中m表示层数，n表示该层声速相应迭代次数，图3给出了一个实际的组织模型。设置各层初始声速C_m0＝1300m/s，随着该层的迭代次数增加，设定其声速为C_mn＝(1300+n)m/s，其中0＜n≤300。

本发明中，步骤三，根据步骤二中设定的一组声速值以及步骤一中还原出的接收数据，用延迟求和法进行图像重建，具体可使用公式(3)计算。

$A\left(\stackrel{\→}{r}\right)=\frac{\underset{k}{\mathrm{\Σ}}w(k,\stackrel{\→}{r})p_0(\stackrel{\→}{r},t+T(k,\stackrel{\→}{r}))}{\underset{k}{\mathrm{\Σ}}w(k,\stackrel{\→}{r})}---\left(3\right)$

其中是空间里的一个像素点，k是相关的阵元数，是权重因子，是延迟时间。

具体实现这个算法的关键问题是对于一个确定的找到其相关的这里，对于成像平面上的一个像素点，根据多层介质的几何尺寸以及假设的声速，计算出该点源声波传播到传感器阵列上孔径范围内的各阵元上的时间间隔以此为索引，返回对应传感器阵元中的数据，即为

本发明中，步骤四，提取每幅图像中声源的分布信息，具体可用公式(4)和(5)计算。

对每个像素点的灰度I(m，n)，其中m，n分别表示像素点所在的行与列，令

$P\left[I(m,n)\right]=\left\{\begin{array}{c}1,I(m,n)\&GreaterEqual;\mathrm{\&sigma;}\\ 0,I(m,n)<\mathrm{\&sigma;}\end{array}\right.---\left(4\right)$

再令

$\mathrm{\&Gamma;}=\frac{\mathrm{\&Sigma;P}\left[I(m,n)\right]I^2(m,n)}{\mathrm{\&Sigma;P}\left[I(m,n)\right]}---\left(5\right)$

其中Γ是每幅图像中声源分布处单位像素点上的能量强度，表征声源能量的集中程度，σ是对灰度进行判断时所取的阈值。

本发明中，步骤五，依据Γ对声速迭代进行反馈控制，具体可用如下方法。

用长度可变的数组array来存放每次产生的Γ值，待产生新的值Γ后，将其添加至数组array的末端，之后寻找数组中的极大值点，若数组中倒数第二个元素对应于一个极大值，则返回值为0，停止声速的迭代，输出该极大值点所对应的各组织声速为最优解；否则，返回值为1，继续声速的迭代。

本发明的具体流程图如图4所示。

整个流程中，步骤一，光声数据通常是用传感器阵列接收人体组织在受到周期性强度调制的光的照射时所发出的超声波。可利用公式(1)、(2)以及所述操作对经过传感器输出的数据进行逆滤波。

整个流程中，步骤二，为确定不同介质中的声速，可利用所述操作对各层声速进行迭代。

整个流程中，步骤三，为了能够判断所取声速是否最优，需用公式(3)进行图像重建。

整个流程中，步骤四，步骤三中产生的图像需提取其中声源的分布信息，可利用公式(4)、(5)计算。

整个流程中，步骤五，可用所述操作在步骤四中提取出的声源分布信息的基础上，判断当前声速是否最优。

本发明提供了一种基于光声成像技术的多层介质声速计算的思路和方法，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种基于光声成像技术的多层介质声速计算的方法，其特征在于，包括以下步骤：

步骤一，根据传感器脉冲响应函数，对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据；

步骤二，设定组织中不同介质的初始声速及声速迭代范围；

步骤三，以逆滤波数据及所设定组织各层声速为基础，进行光声图像重建；

步骤四，从每次重建得到的光声图像中计算并提取图像中声源的分布；

步骤五，依据声源分布信息确定迭代是否完成，若完成则输出声速。

2.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法，其特征在于，所述的逆滤波是在时域中进行的，从单个N型波经过传感器系统函数响应后的波形特点入手，通过互相关运算，依次确定各个N型波的到达时刻以及主瓣幅值，达到分隔出单个N型波输出波形的目的，还原出每个传感器的原始信号中各N型波之间的时延关系以及幅值大小，N型波是指点源光声信号，得名于其波形形状。

3.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法，其特征在于，初次迭代时设定所有组织声速为初值C，C取值在1300m/s至1600m/s之间。

4.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法，其特征在于，对于成像平面中的每个像素点，根据设定的组织声速以及各介质的几何尺寸计算出该像素点处声波由此传播到相关传感器阵元所需的时间，以此在逆滤波数据中找到与其对应的一系列声压数据，对这些数据加权求和即可重建出一个像素点。

5.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法，其特征在于，所述计算并提取重建图像中声源分布的方法是用一个变量P[I(m，n)]对每个像素点(m，n)进行处理，若其灰度值I(m，n)大于设定的阈值，则P[I(m，n)]置1，反之置0，最后计算I²(m，n)关于P[I(m，n)]的加权平均并对系数归一化，即可得到图像中声源分布的信息，简称为声源集中强度。

6.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法，其特征在于，所述的根据声源分布信息判定迭代完成与否的方法是用一个长度可变的数据结构存放声源集中强度值，每次产生出一个新的值后，将其添加至该数据结构末端，若数据结构中的倒数第二个数是极大值，则停止迭代，输出最优的声速；否则，继续迭代。