基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法技术

本发明专利技术公开了一种基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，首先，根据超声换能器阵列接收到的光声信号构建信息矩阵，信息矩阵中既包括含有直达波的信息，又包含有造成干扰的散射波的信息；然后，根据直达波部分所特有的逆对角线性质，通过矩阵45°旋转、矩阵滤波、矩阵45°逆向旋转还原等操作，可以滤除信息矩阵中的散射波成分；最后，利用时间反转算子，可以从滤波后的信息矩阵中获得光吸收体所发射声波强度的空间分布图像。本发明专利技术提出的方法能够有效地滤除散射波的干扰，同时完整地保留了有效探测信息，从而有效地提高散射媒质内及散射媒质后光声图像重构准确性，具有较高的易用性和广泛的适用性。

【技术实现步骤摘要】
基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法
本专利技术涉及一种基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法，属于光声断层成像技术。
技术介绍
光声断层成像作为一种非侵入性的生物医学成像技术，近几年得到了飞速地发展。光声断层成像既具有声学方法对深层组织分辨率高的优点，又具有光学方法在功能成像、分子成像和成像对比度等方面的优势【1】。然而，光声断层成像通常适用声学特性均匀的媒质，散射媒质的成像依然是光声成像技术面临的一大难题。在散射媒质中，接收到的光声信号中，除了包含来自于成像目标的直达波，还包含经由散射层多次散射后的散射波成分。由于散射体的无规律分布，散射波成分呈现随机性，难以提取图像信息，散射波会给光声图像带来斑点噪声、图像畸变。当散射波成分占据主导地位时，强散射波干扰将严重降低光声成像深度与成像质量、甚至使其不能够得到可辨识的图像。为了解决这一问题，人们提出了各种适用于散射媒质的光声成像方案。例如，利用人们提出采用相关矩阵滤波的方法提升散射媒质中的光声成像质量【2】，然而该方法由于不能充分利用检测到的光声数据，存在成像区域窄、分辨率低等局限性。参考文献【1】L.V.WangandS.Hu,“Photoacoustictomography:invivoimagingfromorganellestoorgans,”Science335(6075),1458–1462(2012)【2】WeiRui,ChaoTao,XiaojunLiu,“Photoacousticimaginginscatteringmediabycombiningacorrelationmatrixfilterwithatimereversaloperator,”OpticsExpress25(19),22840.
技术实现思路
专利技术目的：为了提升声散射媒质中光声断层成像的图像质量，本专利技术提供一种基于矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法，以有效地提升声散射媒质中的光声成像的图像质量。具体说是根据直达声波所对应信息矩阵的逆对角线特性，通过矩阵45°逆时针旋转、补零、矩阵滤波、矩阵45°顺时针旋转还原等操作，滤除散射波的干扰信号，然后进行时间反转成像，从而减少散射声波的干扰、提升散射媒质中的光声成像的图像质量。技术方案：为实现上述目的，本专利技术采用的技术方案为：基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，包括如下步骤：步骤1：使用脉冲激光照射被成像样品，样品中的光吸收体吸收脉冲激光的能量后，由于光声效应而向周围媒质辐射超声波；步骤2：利用超声换能器阵列接收样品辐射出的超声波信号，所述超声换能器阵列具有N个沿x轴方向直线排布的换能器单元，相邻两个换能器单元的中心间距为w，第n个换能器单元的横向坐标为：步骤3：对接收到的超声波信号进行包括时间选通、短时傅里叶变换在内的操作，构造大小为N×N的信息矩阵K；步骤4：对信息矩阵K顺序进行包括45°逆时针旋转、补零、滤波、45°顺时针旋转还原在内的操作，滤除散射波成分，保留直达波成分；步骤5：对滤波后的信息矩阵进行奇异值分解，对最大的M个奇异值所对应的特征向量应用时间反转算子，获得对应深度的光声图像；步骤6：根据均匀媒质中的声传递函数，求得修正参数，降低由于有限视角超声信号探测导致的包括图像畸变、虚假对比度在内的图像缺陷；步骤7：改变时间选通区域，重复步骤3～6，获得不同深度的光声成像，最终获得整个区域的光声图像。具体的，所述步骤3中，信息矩阵K的构造过程如下：步骤31：截取超声换能器接收到的处于时间窗口内的信号，并对其作傅里叶变换，得到第n个换能器单元的频域表达Pn(T,f)，该频域表达Pn(T,f)包含直达波信号PD(T,f)与散射波信号PS(T,f)两个部分；其中，T表示成像平面到换能器阵列声波传播的最短时间，也称为时间选通；Δt表示时间窗口的宽度，f表示频率；所谓直达波是指从成像区域出发，未经任何散射体折射、散射，而直接传播到达换能器，并被检测到的信号；所谓散射波是指从从成像区域出发，经过散射体一次或多次的折射、散射后，而传播到达换能器，并被检测到的信号；超声换能器阵列的信息向量记为P(T,f)＝[P1(T,f),…,Pn(T,f),…,PN(T,f)]，直达波信号记为散射波信号PS(T,f)记为第n个换能器单元的直达波信号满足：其中，A0＝Z1/2，k＝2πf/c，c为媒质声速；Z＝cT表示成像平面到换能器阵列的垂直距离，也称为对应时间选通T的深度；X表示成像平面上任意点的x轴坐标值，该任意点表示为(X,Z)；设成像平面上任意点(X,Z)和第n个换能器单元之间存在L条散射路径，则第n个换能器单元的散射波信号满足：其中，Al与sl为第l条散射路径对应的幅值与相位参数；由于散射体为随机分布，因此Al与sl也体现出相应的随机性；步骤32：利用信息向量P构建信息矩阵K：信息矩阵K根据是否含有随机参数分为两个部分：等号右边第一项记作矩阵KC，等号右边的剩余项记作矩阵KR；矩阵KC不含有任何随机参数，即与散射体的分布位置无关，且矩阵元素间具有如下相关性：其中，m＝1,2,…,N，1≤m±q≤N；表示矩阵KC中位于反对角线且距离对角线距离为q的元素之间的相位差，表示矩阵KC中位于(m–q,m+q)位置处的元素，表示矩阵KC中位于(m,m)位置处的元素；矩阵KR中所有元素均含有随机参数，表现出随机性，因此矩阵元素之间不具有相关性。具体的，所述步骤4中，信息矩阵K中散射波成分的滤除操作过程如下：步骤41：不同于文献【1】中的方法，本专利技术专利本专利技术根据下式将信息矩阵K逆时针旋转45°，并对相应元素补零，以保证矩阵维度不变，得到两个对应矩阵A1＝{a1u,v}、A2＝{a2u,v}：若N为奇数：若N为偶数：其中，km,n表示信息矩阵K中位于(m,n)位置处的元素，矩阵A1的维度为N，矩阵A2的维度为N–1；将矩阵A1和矩阵A2统一用矩阵A表示，并用NA表示矩阵A的维度；步骤42：构造滤波矩阵对矩阵A进行滤波，滤波向量C的表达式为：其中，为C的共轭转置；u＝1,2,…,NA，步骤43：对矩阵A进行滤波，将滤波后的矩阵A记为矩阵AF：AF＝FA＝F(AC+AR)＝FAC+FAR其中，矩阵AC为矩阵A中不含随机参数的部分，矩阵AR为矩阵A中含随机参数的部分；第一项FAC＝AC，矩阵AC在滤波前后保持不变，而矩阵AR在滤波用后被强有力地削弱，散射波部分被有效滤除；步骤44：将矩阵AF顺时针旋转45°，得到滤波后的信息矩阵K，将滤波后的信息矩阵K记为矩阵KF：若(m-n)是奇数，若(m-n)是偶数，其中，表示矩阵KF中位于(m,n)位置处的元素，表示矩阵A1F中位于((m-n)/2+(N+3)/4,(m+n)/2)位置处的元素，表示矩阵A2F中位于((m-n-1)/2+(N+3)/4,(m+n-1)/2)位置处的元素，矩阵A1F为滤波后的矩阵A1，矩阵A2F为滤波后的矩阵A2；矩阵KF的大小仍为N×N，与信息矩阵K的一致，滤波前后没有信息维度的缺失。具体的，所述步骤5中，根据滤波后的信息矩阵进行时间反转成像过程如下：步骤51：对矩阵KF作奇异值分解：本文档来自技高网...

【技术保护点】
1.一种基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，其特征在于：包括如下步骤：步骤1：使用脉冲激光照射被成像样品，样品中的光吸收体吸收脉冲激光的能量后，由于光声效应而向周围媒质辐射超声波；步骤2：利用超声换能器阵列接收样品辐射出的超声波信号，所述超声换能器阵列具有N个沿x轴方向直线排布的换能器单元，相邻两个换能器单元的中心间距为w，第n个换能器单元的横向坐标为：

【技术特征摘要】
1.一种基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，其特征在于：包括如下步骤：步骤1：使用脉冲激光照射被成像样品，样品中的光吸收体吸收脉冲激光的能量后，由于光声效应而向周围媒质辐射超声波；步骤2：利用超声换能器阵列接收样品辐射出的超声波信号，所述超声换能器阵列具有N个沿x轴方向直线排布的换能器单元，相邻两个换能器单元的中心间距为w，第n个换能器单元的横向坐标为：步骤3：对接收到的超声波信号进行包括时间选通、短时傅里叶变换在内的操作，构造大小为N×N的信息矩阵K；步骤4：对信息矩阵K顺序进行包括45°逆时针旋转、补零、滤波、45°顺时针旋转还原在内的操作，滤除散射波成分，保留直达波成分；步骤5：对滤波后的信息矩阵进行奇异值分解，对最大的M个奇异值所对应的特征向量应用时间反转算子，获得对应深度的光声图像；步骤6：根据均匀媒质中的声传递函数，求得修正参数，降低由于有限视角超声信号探测导致的包括图像畸变、虚假对比度在内的图像缺陷；步骤7：改变时间选通区域，重复步骤3～6，获得不同深度的光声成像，最终获得整个区域的光声图像。2.根权利要求1所述的基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，其特征在于：所述步骤3中，信息矩阵K的构造过程如下：步骤31：截取超声换能器接收到的处于时间窗口内的信号，并对其作傅里叶变换，得到第n个换能器单元的频域表达Pn(T,f)，该频域表达Pn(T,f)包含直达波信号PD(T,f)与散射波信号PS(T,f)两个部分；其中，T表示成像平面到换能器阵列声波传播的最短时间，也称为时间选通；Δt表示时间窗口的宽度，f表示频率；超声换能器阵列的信息向量记为P(T,f)＝[P1(T,f),…,Pn(T,f),…,PN(T,f)]，直达波信号记为散射波信号PS(T,f)记为第n个换能器单元的直达波信号满足：其中，A0＝Z1/2，k＝2πf/c，c为媒质声速；Z＝cT表示成像平面到换能器阵列的垂直距离，也称为对应时间选通T的深度；X表示成像平面上任意点的x轴坐标值，该任意点表示为(X,Z)；设成像平面上任意点(X,Z)和第n个换能器单元之间存在L条散射路径，则第n个换能器单元的散射波信号满足：其中，Al与sl为第l条散射路径对应的幅值与相位参数；步骤32：利用信息向量P构建信息矩阵K：信息矩阵K根据是否含有随机参数分为两个部分：等号右边第一项记作矩阵KC，等号右边的剩余项记作矩阵KR；矩阵KC不含有任何随机参数，即与散射体的分布位置无关，且矩阵元素间具有如下相关性：其中，m＝1,2,…,N，1≤m±q≤N；表示矩阵KC中位于反对角线且距离对角线距离为q的元素之间的相位差，表示矩阵KC中位于(m–q,m+q)位置处的元素，表示矩阵KC中位于(m,m)位置处的元素；矩阵KR中所有元素均含有随机参数，表现出随机性，因此矩阵元素之间不具有相关性。3.根权利要求1所述的基于全矩阵滤波和时间反转算子的光声断层成像的...

【专利技术属性】
技术研发人员：陶超，芮纬，刘晓峻，
申请(专利权)人：南京大学，
类型：发明
国别省市：江苏,32