买专利卖专利找龙图腾，真高效！ 查专利查商标用IPTOP,全免费！专利年费监控用IP管家,真方便！

申请/专利权人：北京工业大学

摘要：本发明公开了基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法，属于医学图像处理领域。当短脉冲激光照射到待测的生物组织上时，组织吸收激光能量产生热，忽略热传导，构建光声成像的近似热方程；基于MATLAB的开源工具箱来求解光声波方程。在成像域的边界处的光声信号波被一组换能器收集。收集这些信号的过程表示为时变因果系统。使用k‑Wave工具箱来构建该系统的系统矩阵，脉冲响应IR针对完整成像域被逐个像素记录。本方法能够保证成像质量的同时提高成像速度不但可以对光声信号进行准确重建，同时有极高的计算效率。

主权项：基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法，其特征在于：当短脉冲激光照射到待测的生物组织上时，组织吸收激光能量产生热，忽略热传导，此时光声成像的近似热方程为：CptTr,t=Hr,t---1]]其中r是三维空间位置矢量，t是时间，函数Hr,t表示组织吸收的热量，函数Tr,t表示组织升高的温度，ρ是组织的密度，Cp是比热；组织受热膨胀，从而形成一个初始声场，且初始声场将随着时间变化而变化；此时得到光声成像声场的运动方程和扩散方程：2t2ur,t=-pr,t---2]]ur,t=-pr,tc2+Tr,t---3]]矢量函数ur,t表示声位移，函数pr,t表示声压，c是声速，β是等压膨胀系数，▽与▽·分别为梯度与散度的标识；联立式1‑3，将函数Hr,t写为光吸收分布函数Ar和入射激光关于时间的分布函数It的乘积，得：2pr,t-1c22t2pr,t=-CpArtIt---4]]式4是光声成像的基本方程；式4表示光声信号和组织的光吸收特性之间的关系，是描述光声成像的前向模型；通常在实际成像过程中脉冲激光脉宽很短，理论上把光强函数假设为一个脉冲函数即It＝δt；用Green函数法对式4进行求解：Pr,i=4Cptdr|r-r|Art-|r-r|c---5]]因为P0r＝ΓAr′，其中Γ是格留乃森系数因此得到初始光声信号与传感器接收到的光声信号的关系是：Pr,t=14c2t1ctdrP0rt-|r-r|c---6]]使用k‑Wave：基于MATLAB的开源工具箱来求解光声波方程；在成像域的边界处的光声信号波被一组换能器收集；收集这些信号的过程表示为时变因果系统；使用k‑Wave工具箱来构建该系统的系统矩阵，脉冲响应IR针对完整成像域被逐个像素记录，对于尺寸为N*N的成像区域按堆栈顺序矢量化为N2*1，N为成像区域的边长；用x表示，也就是在逆问题中重建的值；将时变光声数据堆叠在M*1维向量中，用b表示测量向量；系统矩阵A每列表示图像的相应像素的脉冲响应IR；光声成像的前向模型归纳为：Ax＝b     7考虑到系统矩阵A具有严重的病态性，同时由于测量数据的不足，导致式7是一个病态的欠定方程；为减弱方程的病态性，PAT重建问题使用正则化方法求解；基于正则化思想，将初始光声信号的求解即PAT重建问题转换为如下形式的最小二乘问题：=||AX-bmeas||22+||X||22---8]]在上式中，A为系统矩阵；bmeas是方程7中与b相对应的边界采集的超声信号；λ是正则化参数；||·||2表示L2范数；系统矩阵A的计算，b的测量借助Matlab的K‑Wave工具箱完成；对A进行奇异值分解，方程8的解表示为:Xtikh=VTVT+I-1VTUTbmeas=VFUbmeas---9]]其中U＝u1,…,un，V＝v1,…,vn是左右奇异正交矩阵，且满足UTU＝VTV＝In，∑＝diagσ1,…,σn，σ1≥…≥σn＞0为A的奇异值，φi是正则化滤波因子其形式是实质是在最小二乘解的基础上增加了约束高频分量的滤波因子；但正则化参数对重建的光声图像有着重要的影响；指数滤波因子与Tikhonov正则化奇异值分解的方法类似，不同的是滤波因子是：φi＝1‑exp‑σi2λ；通过比较这两种滤波因子发现，当σ2＜＜λ时，Tikhonov正则化看作是指数滤波的近似。

全文数据：基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法技术领域[0001]本发明属于医学图像处理领域，涉及一种基于Lanczos双对角化的指数滤波正则化快速光声成像重建方法。背景技术[0002]医学影像作为一种辅助的医疗手段，可以帮助医生检测、判定、认识和研究疾病。传统的医学影像技术有计算机断层成像（ComputedTomography，CT、超生成像UltrasoundImaging，USI、核磁共振成像（MagneticResonanceImaging，MRI、扩散光层析成像diffuseopticaltomography，D0T这些成像技术的发展极大地推动了生命科学的进步，但他们也有一定的局限性。[0003]光声成像（以下简称PAT作为一种新兴的医学成像技术，是一种基于生物组织内部光学吸收差异，以超声为媒质的无损生物光子成像方法。PAT结合了纯光学成像高对比度和纯超声成像高穿透深度特性的优点，以超声探测器探测光声信号代替光学成像中的光子检测，从原理上避免了光学散射的影响，可以在保持生物样本自然状态下得到高对比度、高分辨率的组织图像。由于不同生理状态的生物组织对光的吸收不同，光声图像同时反映了组织代谢的差异和病变特征。近年来光声成像技术被广泛应用于肿瘤检测、血管成像、脑部结构和肾功能成像等领域。[0004]由于在实际测量过程中边界数据有限以及测量的光声信号中不可避免的混有噪声，导致基于模型的PAT逆问题在数学上是一个高度病态的问题，外界微小的测量扰动都会给重建结果带来很大的变化。因此，降低PAT的病态性，如何快速、准确地重建光声信号是PAT研究的重点及热点。[0005]为了降低PAT重建问题的病态性，在重建时可以使用正则化的求解方法将光声重建问题转变成一个非线性的最优化问题。正则化方法最早在上世纪60年代由苏联科学院院士Tikhonov提出，用来求解形如Fx=y方程的不适定问题。正则化的核心思想是在求解问题过程中引入先验信息，从而得到对原问题有意义的解。目前有L2范数、总变分（totalvariation，TV范数和Lp范数的正则化项形式。最常用的方法是基于Tikhonov的正则化重建方法，但其准确性依赖于正则化参数的选取。近年来又提出了指数滤波方法的正则化算法，不仅可以有效地去除重建的光声图像中的高频噪声可以看作一个低通滤波器），而且还具有对正则化参数选择偏向较小的优点。在数据有限的情况下提高了图像质量，但是由于系统矩阵维数较大，导致计算时间过长，并不能实现实时性的需求。因此，可以通过先降维再进行计算的方法进行重建。发明内容[0006]本发明的目的在于提出了一种基于Ianczos双对角化的指数滤波的正则化重建方法，以保证成像质量的同时提高成像速度。[0007]为实现上述目的，本发明采用的技术方案为一种基于lanczos双对角化的指数滤波的正则化重建方法，当短脉冲激光照射到待测的生物组织上时，组织吸收激光能量产生热，忽略热传导，此时光声成像的近似热方程为：[0009]其中r是三维空间位置矢量，t是时间，函数Hr,t表示组织吸收的热量，函数Tr,t表示组织升高的温度，P是组织的密度，CP是比热。组织受热膨胀，从而形成一个初始声场，且初始声场将随着时间变化而变化。此时得到光声成像声场的运动方程和扩散方程：[0012]矢量函数ur，t表示声位移，函数pr，t表示声压，c是声速，β是等压膨胀系数，▽与▽·分别为梯度与散度的标识。联立式（1-3，将函数Hr，t写为光吸收分布函数Ar和入射激光关于时间的分布函数It的乘积，得：[0014]式4是光声成像的基本方程。式4表示光声信号和组织的光吸收特性之间的关系，是描述光声成像的前向模型。[0015]通常在实际成像过程中脉冲激光脉宽很短，理论上把光强函数假设为一个脉冲函数即It=δt;用Green函数法对式4进行求解：[0017]因为PQr=ΓΑ〆），其中Γ是格留乃森系娄，因此得到初始光声信号与传感器接收到的光声信号的关系是：[0019]使用k-WaVe:基于MATLAB的开源工具箱来求解光声波方程。在成像域的边界处的光声信号波被一组换能器收集。收集这些信号的过程表示为时变因果系统。使用k-Wave工具箱来构建该系统的系统矩阵，脉冲响应IR针对完整成像域被逐个像素记录，对于尺寸为N*N的成像区域按堆栈顺序矢量化为N2*1，N为成像区域的边长。用X表示，也就是在逆问题中重建的值。将时变光声数据堆叠在M*1维向量中，用b表示测量向量。系统矩阵A每列表示图像的相应像素的脉冲响应IR。[0020]光声成像的前向模型归纳为：[0021]Ax=b7[0022]考虑到系统矩阵A具有严重的病态性，同时由于测量数据的不足，导致式⑺是一个病态的欠定方程。为减弱方程的病态性，PAT重建问题使用正则化方法求解。[0023]基于正则化思想，将初始光声信号的求解即PAT重建问题转换为如下形式的最小二乘问题：[0025]在上式中，A为系统矩阵;bmeas是方程⑺中与b相对应的边界采集的超声信号;λ是正则化参数；11·Ih表示L2范数。[0026]系统矩阵A的计算，b的测量借助Matlab的K-Wave工具箱完成。[0027]对A进行奇异值分解，方程⑻的解表示为：[0029]其中U=ui，···，un，V=VI，···，vn是左右奇异正交矩阵，且满足=CliagO1^On,O1多…彡ση〇为A的奇异值：，Φί是正则化滤波因子其形式是。实质是在最小二乘解的基础上增加了约束高频分量的滤波因子。但正则化参数对重建的光声图像有着重要的影响。指数滤波因子与Tikhonov正则化奇异值分解的方法类似，不同的是滤波因子是，通过比较这两种滤波因子发现，当σ2λ时，Tikhonov正则化看作是指数滤波的近似。[0030]在对系统矩阵A进行奇异值分解时，考虑到光声成像系统矩阵比较大因此计算比较耗时。为此，先对系统矩阵进行Lanczos双对角化处理以对系统矩阵进行降维。对于MXN2维矩阵A，经过Lanczos双对角化处理后得到k+1*k维的下双对角矩阵B和两个标准正交矩阵Uk+i、Vk，其中：[0032]且满足UtAV=Bk。经过k步迭代后，原始问题的残差表示为：[0035]由于Uk+1是标准正交矩阵，因此原问题就等价于求解下面的最小化问题：[0037]只要求出.V;则原问题的解：附图说明[0039]图1为仿体初始的光声信号分布，图中面积较大的黑色方形区域为成像区域，黑色区域中面积较小的两个白色区域为光声信号的初始分布初始压力分布为IkPa;[0040]图2为通过本方法重建的光声信号分布结果。具体实施方式[0041]下面根据具体实施示例与附图对本发明进行说明。[0042]首先，通过MatIab的工具箱κ-wave建立成像区域。本实验使用方形区域的成像区域，初始光声信号分布如图1所示。在仿真实验中，假定介质具有均匀的超声特性，并且在传播过程中没有声音的吸收和散射。声速为1500ms软组织的近似值）。试验中共设置了40个超声波检测器，将这些检测器等距离的放置在一个半径为22毫米的圆上，成像域的大小设置为101X101像素每个像素大小为0.Imm。在实验中，光声信号采集的时间间隔是50ns，米集时间是25ys。[0043]使用基于Lanczos双对角化的快速指数滤波正则化重建方法进行光声信号的重建，在实验中设置迭代次数K=25，正则化参数λ为6*10_3，实验证明在取该正则化系数值时能得到较好的重建效果。[0044]本发明使用基于Lanczos双对角化的快速指数滤波正则化重建方法。迭代次数K=25,通过本方法计算，可以得出如图2所示的光声信号重建结果，重建时间为23s。实验结果表明，本方法不但可以对光声信号进行准确重建，同时有极高的计算效率。

权利要求：1.基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法，其特征在于：当短脉冲激光照射到待测的生物组织上时，组织吸收激光能量产生热，忽略热传导，此时光声成像的近似热方程为：其中r是三维空间位置矢量，t是时间，函数Hr，t表示组织吸收的热量，函数Tr，t表示组织升高的温度，P是组织的密度，〇是比热；组织受热膨胀，从而形成一个初始声场，且初始声场将随着时间变化而变化;此时得到光声成像声场的运动方程和扩散方程：矢量函数ur，t表示声位移，函数pr，t表示声压，c是声速，β是等压膨胀系数，▽与▽•分别为梯度与散度的标识;联立式（1-3，将函数Hr，t写为光吸收分布函数Ar和入射激光关于时间的分布函数It的乘积，得：式4是光声成像的基本方程;式4表示光声信号和组织的光吸收特性之间的关系，是描述光声成像的前向模型；通常在实际成像过程中脉冲激光脉宽很短，理论上把光强函数假设为一个脉冲函数即I⑴=δt;用Green函数法对式⑷进行求解：因为P〇r=rAr〇，其中Γ是格留乃森系数，因此得到初始光声信号与传感器接收到的光声信号的关系是：使用k-Wave:基于MATLAB的开源工具箱来求解光声波方程;在成像域的边界处的光声信号波被一组换能器收集;收集这些信号的过程表示为时变因果系统;使用k-Wave工具箱来构建该系统的系统矩阵，脉冲响应IR针对完整成像域被逐个像素记录，对于尺寸为N*N的成像区域按堆栈顺序矢量化为N2*1，N为成像区域的边长;用X表示，也就是在逆问题中重建的值;将时变光声数据堆叠在M*1维向量中，用b表示测量向量;系统矩阵A每列表示图像的相应像素的脉冲响应IR;光声成像的前向模型归纳为：Ax=b7考虑到系统矩阵A具有严重的病态性，同时由于测量数据的不足，导致式7是一个病态的欠定方程;为减弱方程的病态性，PAT重建问题使用正则化方法求解；基于正则化思想，将初始光声信号的求解即PAT重建问题转换为如下形式的最小二乘问题：在上式中，A为系统矩阵;bmeas是方程7中与b相对应的边界采集的超声信号;λ是正则化参数；II·I|2表示L2范数；系统矩阵A的计算，b的测量借助Matlab的K-Wave工具箱完成；对A进行奇异值分解，方程⑻的解表示为：其中U=Ul，…，Un，V=VI，…，Vn是左右奇异正交矩阵，且满足UtU=VtV=Ιη，Σ=为A的奇异值，，Φi是正则化滤波因子其形式是·，实质是在最小二乘解的基础上增加了约束高频分量的滤波因子;但正则化参数对重建的光声图像有着重要的影响;指数滤波因子与Tikhonov正则化奇异值分解的方法类似，不同的是滤波因子是：=;通过比较这两种滤波因子发现，当〇2λ时，Tikhonov正则化看作是指数滤波的近似。2.根据权利要求1所述的基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法，其特征在于:在对系统矩阵A进行奇异值分解时，考虑到光声成像系统矩阵比较大因此计算比较耗时；为此，先对系统矩阵进行Lanczos双对角化处理以对系统矩阵进行降维；对于MXN2维矩阵A，经过Lanczos双对角化处理后得到k+1*k维的下双对角矩阵B和两个标准正交矩阵Uk+1、Vk，其中：且满足UtAV=Bk;经过k步迭代后，原始问题的残差表示为：由于认+1是标准正交矩阵，因此原问题就等价于求解下面的最小化问题：只要求出¥则原问题的解：

百度查询： 北京工业大学 基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法