[发明专利]基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法

权利要求书

1.基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法，其特征在于：当短脉冲激光照射到待测的生物组织上时，组织吸收激光能量产生热，忽略热传导，此时光声成像的近似热方程为：

ρC p ∂ ∂ t T ( r , t ) = H ( r , t ) - - - ( 1 ) ]]>

其中r是三维空间位置矢量，t是时间，函数H(r,t)表示组织吸收的热量，函数T(r,t)表示组织升高的温度，ρ是组织的密度，Cp是比热；组织受热膨胀，从而形成一个初始声场，且初始声场将随着时间变化而变化；此时得到光声成像声场的运动方程和扩散方程：

ρ ∂ 2 ∂ t 2 u ( r , t ) = - ▿ p ( r , t ) - - - ( 2 ) ]]>

▿ · u ( r , t ) = - p ( r , t ) ρ · c 2 + β · T ( r , t ) - - - ( 3 ) ]]>

矢量函数u(r,t)表示声位移，函数p(r,t)表示声压，c是声速，β是等压膨胀系数，▽与▽·分别为梯度与散度的标识；联立式(1)-(3)，将函数H(r,t)写为光吸收分布函数A(r)和入射激光关于时间的分布函数I(t)的乘积，得：

▿ 2 p ( r , t ) - 1 c 2 ∂ 2 ∂ t 2 p ( r , t ) = - β C p A ( r ) ∂ ∂ t I ( t ) - - - ( 4 ) ]]>

式(4)是光声成像的基本方程；式(4)表示光声信号和组织的光吸收特性之间的关系，是描述光声成像的前向模型；

通常在实际成像过程中脉冲激光脉宽很短，理论上把光强函数假设为一个脉冲函数即I(t)＝δ(t)；用Green函数法对式(4)进行求解：

P ( r , i ) = β 4 πC p ∂ ∂ t ∫ dr ′ | r - r ′ | A ( r ′ ) δ ′ ( t - | r - r ′ | c ) - - - ( 5 ) ]]>

因为P0(r)＝ΓA(r′)，其中Γ是格留乃森系数因此得到初始光声信号与传感器接收到的光声信号的关系是：

P ( r , t ) = 1 4 πc 2 ∂ ∂ t [ 1 c t ∫ dr ′ P 0 ( r ′ ) δ ′ ( t - | 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重建问题转换为如下形式的最小二乘问题：

Ω = | | A X - b m e a s | | 2 2 + λ | | X | | 2 2 - - - ( 8 ) ]]>

在上式中，A为系统矩阵；bmeas是方程(7)中与b相对应的边界采集的超声信号；λ是正则化参数；||·||2表示L2范数；

系统矩阵A的计算，b的测量借助Matlab的K-Wave工具箱完成；

对A进行奇异值分解，方程(8)的解表示为:

X t i k h = ( V Σ T Σ V T + λ I ) - 1 V Σ T U T b m e a s = V Σ F Ub m e a s - - - ( 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正则化看作是指数滤波的近似。

2.根据权利要求1所述的基于Lanczos双对角化的快速指数滤波正则化光声成像重建方法，其特征在于：在对系统矩阵A进行奇异值分解时，考虑到光声成像系统矩阵比较大因此计算比较耗时；为此，先对系统矩阵进行Lanczos双对角化处理以对系统矩阵进行降维；对于M×N2维矩阵A，经过Lanczos双对角化处理后得到(k+1)*k维的下双对角矩阵B和两个标准正交矩阵Uk+1、Vk，其中：

Uk+1＝(u1,u2,...,uk+1)，Vk＝(v1,v2,...,vk)；且满足UTAV＝Bk；经过k步迭代后，原始问题的残差表示为：

rk＝AX-b＝AVkyk-β1u1＝Uk+1(Bkyk-β1e1) (11)

由于Uk+1是标准正交矩阵，因此原问题就等价于求解下面的最小化问题：

min y k | | B k y k - β 1 e 1 | | 2 2 - - - ( 12 ) ]]>

只要求出则原问题的解：

X*＝Vky* (13)。