CN103142216A - 一种基于光声成像技术的多层介质声速计算的方法 - Google Patents
一种基于光声成像技术的多层介质声速计算的方法 Download PDFInfo
- Publication number
- CN103142216A CN103142216A CN201310113624XA CN201310113624A CN103142216A CN 103142216 A CN103142216 A CN 103142216A CN 201310113624X A CN201310113624X A CN 201310113624XA CN 201310113624 A CN201310113624 A CN 201310113624A CN 103142216 A CN103142216 A CN 103142216A
- Authority
- CN
- China
- Prior art keywords
- sound
- velocity
- iteration
- photoacoustic
- imaging technology
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000003384 imaging method Methods 0.000 title claims abstract description 31
- 238000005516 engineering process Methods 0.000 title claims abstract description 18
- 238000005316 response function Methods 0.000 claims abstract description 5
- 230000004044 response Effects 0.000 claims description 8
- 239000012141 concentrate Substances 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 2
- 238000001914 filtration Methods 0.000 abstract 3
- 230000000875 corresponding effect Effects 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 230000003287 optical effect Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 241001465754 Metazoa Species 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 2
- 230000002596 correlated effect Effects 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 102000001554 Hemoglobins Human genes 0.000 description 1
- 108010054147 Hemoglobins Proteins 0.000 description 1
- 241000699670 Mus sp. Species 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 210000004556 brain Anatomy 0.000 description 1
- 210000000481 breast Anatomy 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 210000005056 cell body Anatomy 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000010895 photoacoustic effect Methods 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 238000009774 resonance method Methods 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000002604 ultrasonography Methods 0.000 description 1
Images
Landscapes
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
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,最后计算I2(m,n)关于P[I(m,n)]的加权平均,即可得到图像中声源分布的信息,简称为声源集中强度。
本发明中,优选地,所述根据声源分布信息判定迭代完成与否的方法是用一个长度可变的数组存放声源集中强度值,每次产生出一个新的值后,将其添加至数组末端,若数组中的倒数第二个数是极大值,则停止迭代,输出最优的声速;否则,继续迭代。
附图说明
下面结合附图和具体实施方式对本发明做更进一步的具体说明,本发明的上述和/或其他方面的优点将会变得更加清楚。
图1是本发明方法中N型波的时域波形。
图2是本发明方法中N型波经***函数响应后的波形。
图3是本发明方法的一个具体实例所使用的组织结构图。
图4是本发明方法的流程图。
具体实施方式:
本发明公开了一种基于光声成像技术的多层介质声速计算的方法,包括以下步骤:
步骤一,根据传感器脉冲响应函数,对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据;
步骤二,设定组织中不同介质的初始声速及声速迭代范围;
步骤三,以逆滤波数据及所设定组织各层声速为基础,进行光声图像重建;
步骤四,从每次重建得到的光声图像中计算并提取图像中声源的分布;
步骤五,依据声源分布信息确定迭代是否完成,若完成则输出声速。
本发明中,步骤一,单个N型波经过传感器***函数响应后的时域信号可由公式(1)计算得到。
其中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(t)相对于pk(t)的延迟时间,pk(t)表示第k个传感器阵元所输出的信号。之后找到每路Rk(τ)中的极大值,并返回对应的τ值,记为这里的变量用下标k区分各传感器阵元,用上标表示同一传感器阵元中同一变量在不同时刻的值。
然后对于每路pk(t),认为为第一个N型波的响应到来的时刻,根据主瓣幅值定位其波形利用r(t)与f(t)波形关系的先验知识,还原出的波形,接着用pk(t)减去得对于认为为第二个N型波的响应到来的时刻,根据主瓣幅值定位其波形还原出的波形,接着用减去得依次进行下去即还原出一个传感器阵元所接收到的信号。
本发明中,步骤二,对可能的声速进行设定具体可按如下方式进行。
对于每层介质的声速Cmn,其中m表示层数,n表示该层声速相应迭代次数,图3给出了一个实际的组织模型。设置各层初始声速Cm0=1300m/s,随着该层的迭代次数增加,设定其声速为Cmn=(1300+n)m/s,其中0<n≤300。
本发明中,步骤三,根据步骤二中设定的一组声速值以及步骤一中还原出的接收数据,用延迟求和法进行图像重建,具体可使用公式(3)计算。
具体实现这个算法的关键问题是对于一个确定的找到其相关的这里,对于成像平面上的一个像素点,根据多层介质的几何尺寸以及假设的声速,计算出该点源声波传播到传感器阵列上孔径范围内的各阵元上的时间间隔以此为索引,返回对应传感器阵元中的数据,即为
本发明中,步骤四,提取每幅图像中声源的分布信息,具体可用公式(4)和(5)计算。
对每个像素点的灰度I(m,n),其中m,n分别表示像素点所在的行与列,令
再令
其中Γ是每幅图像中声源分布处单位像素点上的能量强度,表征声源能量的集中程度,σ是对灰度进行判断时所取的阈值。
本发明中,步骤五,依据Γ对声速迭代进行反馈控制,具体可用如下方法。
用长度可变的数组array来存放每次产生的Γ值,待产生新的值Γ后,将其添加至数组array的末端,之后寻找数组中的极大值点,若数组中倒数第二个元素对应于一个极大值,则返回值为0,停止声速的迭代,输出该极大值点所对应的各组织声速为最优解;否则,返回值为1,继续声速的迭代。
本发明的具体流程图如图4所示。
整个流程中,步骤一,光声数据通常是用传感器阵列接收人体组织在受到周期性强度调制的光的照射时所发出的超声波。可利用公式(1)、(2)以及所述操作对经过传感器输出的数据进行逆滤波。
整个流程中,步骤二,为确定不同介质中的声速,可利用所述操作对各层声速进行迭代。
整个流程中,步骤三,为了能够判断所取声速是否最优,需用公式(3)进行图像重建。
整个流程中,步骤四,步骤三中产生的图像需提取其中声源的分布信息,可利用公式(4)、(5)计算。
整个流程中,步骤五,可用所述操作在步骤四中提取出的声源分布信息的基础上,判断当前声速是否最优。
本发明提供了一种基于光声成像技术的多层介质声速计算的思路和方法,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。
Claims (6)
1.一种基于光声成像技术的多层介质声速计算的方法,其特征在于,包括以下步骤:
步骤一,根据传感器脉冲响应函数,对传感器接收到的原始数据进行逆滤波处理得到逆滤波数据;
步骤二,设定组织中不同介质的初始声速及声速迭代范围;
步骤三,以逆滤波数据及所设定组织各层声速为基础,进行光声图像重建;
步骤四,从每次重建得到的光声图像中计算并提取图像中声源的分布;
步骤五,依据声源分布信息确定迭代是否完成,若完成则输出声速。
2.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法,其特征在于,所述的逆滤波是在时域中进行的,从单个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,最后计算I2(m,n)关于P[I(m,n)]的加权平均并对系数归一化,即可得到图像中声源分布的信息,简称为声源集中强度。
6.根据权利要求1所述的一种基于光声成像技术的多层介质声速计算的方法,其特征在于,所述的根据声源分布信息判定迭代完成与否的方法是用一个长度可变的数据结构存放声源集中强度值,每次产生出一个新的值后,将其添加至 该数据结构末端,若数据结构中的倒数第二个数是极大值,则停止迭代,输出最优的声速;否则,继续迭代。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310113624.XA CN103142216B (zh) | 2013-04-03 | 2013-04-03 | 一种基于光声成像技术的多层介质声速计算的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310113624.XA CN103142216B (zh) | 2013-04-03 | 2013-04-03 | 一种基于光声成像技术的多层介质声速计算的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103142216A true CN103142216A (zh) | 2013-06-12 |
CN103142216B CN103142216B (zh) | 2014-11-12 |
Family
ID=48540768
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310113624.XA Expired - Fee Related CN103142216B (zh) | 2013-04-03 | 2013-04-03 | 一种基于光声成像技术的多层介质声速计算的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103142216B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105249993A (zh) * | 2015-11-16 | 2016-01-20 | 南京大学 | 一种通过光声成像选取最佳声速组优化超声成像的方法 |
CN111214213A (zh) * | 2020-02-13 | 2020-06-02 | 南京科技职业学院 | 一种适用于声速不均匀介质的光声断层成像方法 |
CN113777045A (zh) * | 2020-06-10 | 2021-12-10 | 复旦大学 | 基于单粒子多边定位追踪的超分辨功能性光声成像方法 |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3200902B2 (ja) * | 1991-12-24 | 2001-08-20 | 株式会社日立製作所 | 光音響信号検出方法及び装置 |
US6400450B1 (en) * | 2000-03-17 | 2002-06-04 | Fitel Usa Corp. | Method of qualifying a multimode optical fiber for bandwidth performance |
AU2002332365A1 (en) * | 2001-08-06 | 2003-02-24 | Vladimir Pavlovich Zharov | Optical method and device for spatially manipulating objects |
CN101214156A (zh) * | 2008-01-10 | 2008-07-09 | 复旦大学 | 声速不均匀介质热声成像的重建算法 |
CN101251413A (zh) * | 2008-04-17 | 2008-08-27 | 上海交通大学 | 采用边界元法重建循环平稳声源的方法 |
CN102306385A (zh) * | 2011-06-22 | 2012-01-04 | 复旦大学 | 任意扫描方式下光声成像的图像重建方法 |
CN102608036A (zh) * | 2012-03-20 | 2012-07-25 | 中北大学 | 基于声学透镜和传感器阵列的三维光声成像***及方法 |
US20120194646A1 (en) * | 2011-02-02 | 2012-08-02 | National Tsing Hua University | Method of Enhancing 3D Image Information Density |
-
2013
- 2013-04-03 CN CN201310113624.XA patent/CN103142216B/zh not_active Expired - Fee Related
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3200902B2 (ja) * | 1991-12-24 | 2001-08-20 | 株式会社日立製作所 | 光音響信号検出方法及び装置 |
US6400450B1 (en) * | 2000-03-17 | 2002-06-04 | Fitel Usa Corp. | Method of qualifying a multimode optical fiber for bandwidth performance |
AU2002332365A1 (en) * | 2001-08-06 | 2003-02-24 | Vladimir Pavlovich Zharov | Optical method and device for spatially manipulating objects |
CN101214156A (zh) * | 2008-01-10 | 2008-07-09 | 复旦大学 | 声速不均匀介质热声成像的重建算法 |
CN101251413A (zh) * | 2008-04-17 | 2008-08-27 | 上海交通大学 | 采用边界元法重建循环平稳声源的方法 |
US20120194646A1 (en) * | 2011-02-02 | 2012-08-02 | National Tsing Hua University | Method of Enhancing 3D Image Information Density |
CN102306385A (zh) * | 2011-06-22 | 2012-01-04 | 复旦大学 | 任意扫描方式下光声成像的图像重建方法 |
CN102608036A (zh) * | 2012-03-20 | 2012-07-25 | 中北大学 | 基于声学透镜和传感器阵列的三维光声成像***及方法 |
Non-Patent Citations (4)
Title |
---|
向良忠 等: "改进的同步迭代算法在光声血管成像中的应用", 《物理学报》 * |
吴丹 等: "光声成像中延迟求和方法和反投影重构方法的比较", 《无损检测》 * |
张弛 等: "声速不均匀介质的光声成形重建算法", 《光学学报》 * |
梁凤岗: "信号重构中的时域反滤波及其应用", 《振动、测试与诊断》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105249993A (zh) * | 2015-11-16 | 2016-01-20 | 南京大学 | 一种通过光声成像选取最佳声速组优化超声成像的方法 |
CN105249993B (zh) * | 2015-11-16 | 2018-01-02 | 南京大学 | 一种通过光声成像选取最佳声速组优化超声成像的方法 |
CN111214213A (zh) * | 2020-02-13 | 2020-06-02 | 南京科技职业学院 | 一种适用于声速不均匀介质的光声断层成像方法 |
CN111214213B (zh) * | 2020-02-13 | 2022-11-11 | 南京科技职业学院 | 一种适用于声速不均匀介质的光声断层成像方法 |
CN113777045A (zh) * | 2020-06-10 | 2021-12-10 | 复旦大学 | 基于单粒子多边定位追踪的超分辨功能性光声成像方法 |
CN113777045B (zh) * | 2020-06-10 | 2022-10-18 | 复旦大学 | 基于单粒子多边定位追踪的超分辨功能性光声成像方法 |
Also Published As
Publication number | Publication date |
---|---|
CN103142216B (zh) | 2014-11-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP7345208B2 (ja) | 学習非線形マッピングに基づく画像再構成方法 | |
CN102641137B (zh) | 使用幅度-相位调制超声波的粘弹性测量 | |
EP2501295B1 (en) | Image forming method | |
Liu et al. | Automatic mode extraction of ultrasonic guided waves using synchrosqueezed wavelet transform | |
Romano et al. | On the noninvasive determination of material parameters from a knowledge of elastic displacements theory and numerical simulation | |
CN102640014B (zh) | 图像生成装置和图像生成方法 | |
Besson et al. | Ultrafast ultrasound imaging as an inverse problem: Matrix-free sparse image reconstruction | |
EP2903530B1 (en) | Shear wave attenuation from k-space analysis system | |
CN104688224B (zh) | 一种应用于声学非均匀媒介磁声耦合成像重建方法 | |
CN105249993A (zh) | 一种通过光声成像选取最佳声速组优化超声成像的方法 | |
Brath et al. | Acoustic formulation of elastic guided wave propagation and scattering in curved tubular structures | |
CN101874744B (zh) | 用于长骨分析的超声导波参数测量方法 | |
CN109157215A (zh) | 一种基于***矩阵的磁感应磁声电导率图像重建方法 | |
KR102326149B1 (ko) | 모델-기반 이미지 재구성 방법 | |
Zhao et al. | Reconstruction of Lamb wave dispersion curves by sparse representation with continuity constraints | |
CN104013388B (zh) | 基于低频连续波的磁声耦合成像激励与检测方法及装置 | |
CN103142216B (zh) | 一种基于光声成像技术的多层介质声速计算的方法 | |
CN105395219A (zh) | 一种超声光声光声谱三模态成像***及方法 | |
CN114224387B (zh) | 基于超声多径信道特征参数感知的体脂率测量方法 | |
Nagatani et al. | Multichannel instantaneous frequency analysis of ultrasound propagating in cancellous bone | |
Almansouri et al. | Deep neural networks for non-linear model-based ultrasound reconstruction | |
KR20230145566A (ko) | 전파 반전을 이용한 반사 초음파 이미징 | |
JP5419727B2 (ja) | 画像形成方法及び音響波測定装置 | |
WO2019025510A1 (en) | METHOD AND DEVICE FOR CHARACTERIZING A WAVEGUIDE | |
CN104199013B (zh) | 一种能在有限水域内降低测试频率的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20141112 Termination date: 20150403 |
|
EXPY | Termination of patent right or utility model |