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

Abstract

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

Description

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

技术领域

本发明涉及一种基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法，属于光声断层成像技术。

背景技术

光声断层成像作为一种非侵入性的生物医学成像技术，近几年得到了飞速地发展。光声断层成像既具有声学方法对深层组织分辨率高的优点，又具有光学方法在功能成像、分子成像和成像对比度等方面的优势【1】。

然而，光声断层成像通常适用声学特性均匀的媒质，散射媒质的成像依然是光声成像技术面临的一大难题。在散射媒质中，接收到的光声信号中，除了包含来自于成像目标的直达波，还包含经由散射层多次散射后的散射波成分。由于散射体的无规律分布，散射波成分呈现随机性，难以提取图像信息，散射波会给光声图像带来斑点噪声、图像畸变。当散射波成分占据主导地位时，强散射波干扰将严重降低光声成像深度与成像质量、甚至使其不能够得到可辨识的图像。

为了解决这一问题，人们提出了各种适用于散射媒质的光声成像方案。例如，利用人们提出采用相关矩阵滤波的方法提升散射媒质中的光声成像质量【2】，然而该方法由于不能充分利用检测到的光声数据，存在成像区域窄、分辨率低等局限性。

参考文献

【1】L.V.Wang and S.Hu,“Photoacoustic tomography:in vivo imaging fromorganelles to organs,”Science 335(6075),1458–1462(2012)

【2】Wei Rui,Chao Tao,Xiaojun Liu,“Photoacoustic imaging in scatteringmedia by combining a correlation matrix filter with a time reversaloperator,”Optics Express 25(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个换能器单元的频域表达P_n(T,f)，该频域表达P_n(T,f)包含直达波信号P^D(T,f)与散射波信号P^S(T,f)两个部分；其中，T表示成像平面到换能器阵列声波传播的最短时间，也称为时间选通；Δt表示时间窗口的宽度，f表示频率；

所谓直达波是指从成像区域出发，未经任何散射体折射、散射，而直接传播到达换能器，并被检测到的信号；所谓散射波是指从从成像区域出发，经过散射体一次或多次的折射、散射后，而传播到达换能器，并被检测到的信号；

超声换能器阵列的信息向量记为P(T,f)＝[P₁(T,f),…,P_n(T,f),…,P_N(T,f)]，直达波信号记为散射波信号P^S(T,f)记为

第n个换能器单元的直达波信号满足：

其中，A₀＝Z^1/2，k＝2πf/c，c为媒质声速；Z＝cT表示成像平面到换能器阵列的垂直距离，也称为对应时间选通T的深度；X表示成像平面上任意点的x轴坐标值，该任意点表示为(X,Z)；

设成像平面上任意点(X,Z)和第n个换能器单元之间存在L条散射路径，则第n个换能器单元的散射波信号满足：

其中，A_l与s_l为第l条散射路径对应的幅值与相位参数；由于散射体为随机分布，因此A_l与s_l也体现出相应的随机性；

步骤32：利用信息向量P构建信息矩阵K：

信息矩阵K根据是否含有随机参数分为两个部分：等号右边第一项记作矩阵K^C，等号右边的剩余项记作矩阵K^R；

矩阵K^C不含有任何随机参数，即与散射体的分布位置无关，且矩阵元素间具有如下相关性：

其中，m＝1,2,…,N，1≤m±q≤N；表示矩阵K^C中位于反对角线且距离对角线距离为q的元素之间的相位差，表示矩阵K^C中位于(m–q,m+q)位置处的元素，表示矩阵K^C中位于(m,m)位置处的元素；

矩阵K^R中所有元素均含有随机参数，表现出随机性，因此矩阵元素之间不具有相关性。

具体的，所述步骤4中，信息矩阵K中散射波成分的滤除操作过程如下：

步骤41：不同于文献【1】中的方法，本发明专利本发明根据下式将信息矩阵K逆时针旋转45°，并对相应元素补零，以保证矩阵维度不变，得到两个对应矩阵A1＝{a1_u,v}、A2＝{a2_u,v}：

若N为奇数：

若N为偶数：

其中，k_m,n表示信息矩阵K中位于(m,n)位置处的元素，

矩阵A1的维度为N，矩阵A2的维度为N–1；将矩阵A1和矩阵A2统一用矩阵A表示，并用N_A表示矩阵A的维度；

步骤42：构造滤波矩阵对矩阵A进行滤波，滤波向量C的表达式为：

其中，为C的共轭转置；u＝1,2,…,N_A，

步骤43：对矩阵A进行滤波，将滤波后的矩阵A记为矩阵A^F：

A^F＝FA＝F(A^C+A^R)＝FA^C+FA^R

其中，矩阵A^C为矩阵A中不含随机参数的部分，矩阵A^R为矩阵A中含随机参数的部分；第一项FA^C＝A^C，矩阵A^C在滤波前后保持不变，而矩阵A^R在滤波用后被强有力地削弱，散射波部分被有效滤除；

步骤44：将矩阵A^F顺时针旋转45°，得到滤波后的信息矩阵K，将滤波后的信息矩阵K记为矩阵K^F：

若(m-n)是奇数，

若(m-n)是偶数，

其中，表示矩阵K^F中位于(m,n)位置处的元素，表示矩阵A1^F中位于((m-n)/2+(N+3)/4,(m+n)/2)位置处的元素，表示矩阵A2^F中位于((m-n-1)/2+(N+3)/4,(m+n-1)/2)位置处的元素，矩阵A1^F为滤波后的矩阵A1，矩阵A2^F为滤波后的矩阵A2；

矩阵K^F的大小仍为N×N，与信息矩阵K的一致，滤波前后没有信息维度的缺失。

具体的，所述步骤5中，根据滤波后的信息矩阵进行时间反转成像过程如下：

步骤51：对矩阵K^F作奇异值分解：

其中，U^F为N×N阶酉矩阵，它的行向量是矩阵K^F的左奇异向量；Λ^F为半正定N×N阶对角矩阵；为N×N阶酉矩阵，为V^F的共轭转置，V^F的行向量是矩阵K^F的右奇异向量；Λ^F对角线上的元素即为矩阵K^F的奇异值；

步骤52：计算成像平面上任意点(X,Z)与超声换能器阵列之间的直达波传输矩阵G，矩阵G中位于(m,n)位置处的元素g_m,n为：

其中，r_m,n为超声换能器阵列中的第n个换能器单元与第m个换能器单元之间的距离；

步骤53：最大的M个奇异值所对应的特征向量与直达波传输矩阵G应用时间反转算子，第i个奇异值在深度Z处的光声信号强度为：

其中，表示V^F的第i行奇异向量；将最大的M个奇异值所对应求出的光声信号强度叠加，获得深度Z处的光声图像 是I^F向量的第m个元素，它表示成像平面上与第m个换能器单元所对应位置的初步成像值。

具体的，所述步骤6中，修正参数生成与利用过程如下：

步骤61：在深度Z处，第m个换能器单元所对应的成像修正值为：

其中，为g_m,n的共轭；

步骤62：利用成像修正值对初步成像值进行修正，从而获得最终的成像值

即，最终的成像结果由滤波成像结果I^F与修正矩阵I^R对应元素相除得到。

有益效果：本发明提供的基于矩阵滤波和时间反转算子的光声断层成像的图像重构方法，具有如下优点：1、本发明利用超声换能器阵列接收到的信息矩阵，可以有效对声波信号进行滤波处理，减小散射媒质对成像效果的干扰，获得高质量的目标成像；2、本发明仅需一次激发与接收过程，操作简单高效，具有较高的实用性与广泛的适用性；3、本发明在滤波前后不会造成信息的丢失，能够对接收到的全体有效信息加以利用；4、本发明根据均匀媒质中的声传递函数，求得修正参数，降低由于有限视角超声信号探测导致的图像畸变、虚假对比度等图像缺陷。

附图说明

图1为基于矩阵滤波和时间反转算子的优化型光声断层成像的图像重构原理示意图；

图2为矩阵旋转、补零和滤波操作示意图；

图3为滤波前后特征空间矩阵对比图，3(a)为无散射层情况下接收到的信号矩阵K的实部图像，3(b)为在模型中加入散射层后信号矩阵K的实部图像，3(c)为采用文献【2】中滤波方法滤波后的矩阵K^F的实部图像，3(d)为采用本发明所提出的滤波方法滤波后的矩阵K^F的实部图像；

图4为传统滤波成像方法与优化型滤波成像方法成像结果对比图，4(a)为传统波束成形成像法的成像结果，4(b)为未经滤波情况下时间反转法的成像结果，4(c)为文献【2】中滤波成像方法的成像结果，4(d)为本发明所提出的优化型滤波成像方法的成像结果。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示为一种基于矩阵滤波和时间反转算子的光声断层成像的图像重构法的原理示意图，以对位于散射层后的小钢球成像为例，包括如下步骤。

步骤1：将40颗直径为0.8mm的小钢珠随机无规律地排布以作为随机散射层，置于超声换能器与目标区域之间。由于散射体的无规律分布，接收到的光声信号中，除了包含来自于成像目标的直达波，还包含经由散射层多次散射后的散射波成分。散射层用以产生干扰信号。将直径为0.8mm的四个小钢珠置于散射层后作为成像目标。使用脉冲激光器照射被成像目标，成像目标中的小钢珠吸收激光能量后，由于光声效应而向周围媒质辐射超声波(***如图1所示)。

步骤2：利用采样频率为20MHz的超声换能器阵列接收样品辐射出的超声信号，该超声换能器阵列具有N＝101个单元，相邻单元之间的间距为w＝0.5mm，因此第n个单元的横向坐标为x_n＝[n-(N+1)/2]w，此处N为奇数，当N为偶数时，本发明的方法仍然成立。

步骤3：对接收到的超声信号进行时间选通和短时傅里叶变换操作，获得其对应于频率f的信号分量，构造N×N的信息矩阵K。

步骤31：截取超声换能器接收到的处于时间范围内的信号，并对其作短时傅里叶变换，得到频域表达，记做P_n(T,f)，其中下标n表示对应第n个换能器单元。由N个换能器获得的频域信号构成信息向量P(T,f)＝[P₁(T,f),…,P_n(T,f),…,P_N(T,f)]。信息向量P(T,f)包含直达波信号P^D(T,f)与散射波信号P^S(T,f)两个部分(如图1所示)。直达波信号对应元素表达式满足：

其中，A₀＝Z^1/2，k＝2πf/c，c为媒质声速。相应的，散射波信号对应元素表达式满足：

其中，A_l与s_l为第l条散射路径对应的幅值与相位参数；由于散射体为随机分布，因此A_l与s_l也体现出相应的随机性。

步骤32：利用信息向量P构建信息矩阵K：

矩阵K根据是否含有随机参数分为两个部分。等号右边第一项记作K^C，该项不含有任何随机参数，即与散射体的分布位置无关，且矩阵元素间具有如下相关性：

与之相对，等号右边的剩余项中均含有随机参数，统记作K^R。K^R表现出随机性，元素之间不具有K^C所体现的相关性。

步骤4：对信息矩阵K进行45°逆时针旋转、补零、滤波、45°顺时针旋转还原等操作，滤除散射波成分，保留直达波成分。

步骤41：根据下式将信息矩阵逆时针旋转45°，得到对应矩阵A1、A2：

上式中m＝u+v–(N+1)/2，n＝v–u+(N+1)/2。矩阵A1的维度为N，A2的维度为N–1。为了叙述方便，下文统一用A表示A1或A2，并用N_A表示矩阵A的维度。

旋转、补零过程如图2所示，将原信号矩阵逆时针旋转45°。原矩阵的相干性由反对角线方向转移至列方向。然后，根据列数的奇偶性将旋转后的矩阵分为两个部分，浅绿色圆形部分对应奇数列，灰色菱形部分对应偶数列。最后，将两部分分离，然后各自补零(淡蓝色正方形部分)，从而得到矩阵A1、A2。

步骤42：构造滤波矩阵对A矩阵进行滤波。滤波向量C的表达式为：

其中u＝1,2,…,N_A，

步骤43：对矩阵A进行滤波，得到滤波矩阵A^F

A^F＝FA＝F(A^C+A^R)＝FA^C+FA^R

上式中第一项FA^C＝A^C，相关性项A^C在滤波前后保持不变，而随机项A^R在被滤波矩阵作用后被强有力地削弱，散射波部分被有效滤除。

步骤44：将矩阵A^F顺时针旋转45°，得到滤波后的信息矩阵K^F：

若(m-n)是奇数，

若(m-n)是偶数，

滤波后得到的信息矩阵K^F的大小仍为N×N，与原信号矩阵K的一致，滤波前后没有信息维度的缺失。

当选取时间T＝58μs(对应于位于[-5mm,87mm]的小钢珠)，中心频率f＝2.0MHz时，矩阵滤波前后的对比效果如图3所示。图3(a)为无散射层情况下，接收到的信号矩阵K的实部图像。此时矩阵信号仅含有直达波成分，图像展现出清晰的相关性。当我们在模型中加入散射层后，接收信号中的散射波成分占据主导地位，信号矩阵K的实部图像如图3(b)所示，此时图像表现为强烈的随机性。图3(c)为采用文献1中滤波方法滤波后的矩阵K^F的实部图像，图像中心区域的散射波干扰得到了一定程度的滤除，可以观察到相关性图案，然而缺失了四周区域的大量信息。图3(d)为采用本发明所提出的滤波方法滤波后的矩阵K^F的实部图像。可以观察到，一方面散射波成分的干扰得到了更有效的滤除，可以观察到清晰的相关性图案，结果与3(a)高度接近。另一方面，矩阵维度保持不变，没有出现信息的缺失，为之后的成像过程提供了更好的保障。

步骤5：对滤波后的矩阵应用时间反转法获得对应的传播区间的成像。

步骤51：对滤波矩阵K^F作奇异值分解：

步骤52：计算焦平面上的点与超声换能器阵列单元之间的直达波传输矩阵G，其元素表达式为：

式中r_m,n为换能器阵列的第n个单元与焦平面上第m个点之间的距离，m,n＝(1,2,…,N)。

步骤53：对最大的前几个奇异值所对应的特征向量与传输矩阵G应用时间反转算子，第i个奇异值在深度Z处的光声信号强度为：

将前M个非零奇异值所对应求出的光声信号强度叠加，获得对应深度的光声图像

步骤6：根据均匀媒质中的声传递函数，求得修正参数，降低由于有限视角超声信号探测导致的图像畸变、虚假对比度等图像缺陷。

步骤61：在指定的深度Z处，焦平面上第m个点所对应的成像修正值由下式获得：

步骤62：利用该修正值对初步成像值进行修正，从而获得最终的成像值最终的成像结果由滤波成像结果I^F与修正矩阵I^R对应元素相除得到：

步骤7：改变时间选通区域，重复步骤3～6，获得不同深度的光声成像，最终获得整个区域的光声图像。

图4为多种成像方法成像结果的对比图，图4(a)，4(b)分别为传统波束成形成像法以及未经滤波情况下时间反转法的成像结果。可以看出，这两种方法的成像质量较差，背景区域存在大量干扰图像与畸变，散射波成分对成像质量的影响十分明显。图4(c)为文献1中滤波成像方法的成像结果，光吸收体与背景区域的对比度明显提高，滤波过程对成像质量的改善效果得以体现。然而该方法仅对中心区域实现了成像，上下两侧区域的两个光吸收体未能得到成像。图4(d)为本发明所提出的优化型滤波成像方法的成像结果。由于有效信息的利用率得到了提升，光吸收体与背景区域的对比度进一步提高，同时获得了完整区域的成像结果，所有光吸收体均得到了有效成像。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (5)

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个换能器单元的频域表达P_n(T,f)，该频域表达P_n(T,f)包含直达波信号P^D(T,f)与散射波信号P^S(T,f)两个部分；其中，T表示成像平面到换能器阵列声波传播的最短时间，也称为时间选通；Δt表示时间窗口的宽度，f表示频率；

超声换能器阵列的信息向量记为P(T,f)＝[P₁(T,f),…,P_n(T,f),…,P_N(T,f)]，直达波信号记为散射波信号P^S(T,f)记为

第n个换能器单元的直达波信号满足：

其中，A₀＝Z^1/2，k＝2πf/c，c为媒质声速；Z＝cT表示成像平面到换能器阵列的垂直距离，也称为对应时间选通T的深度；X表示成像平面上任意点的x轴坐标值，该任意点表示为(X,Z)；

设成像平面上任意点(X,Z)和第n个换能器单元之间存在L条散射路径，则第n个换能器单元的散射波信号满足：

其中，A_l与s_l为第l条散射路径对应的幅值与相位参数；

步骤32：利用信息向量P构建信息矩阵K：

信息矩阵K根据是否含有随机参数分为两个部分：等号右边第一项记作矩阵K^C，等号右边的剩余项记作矩阵K^R；

矩阵K^C不含有任何随机参数，即与散射体的分布位置无关，且矩阵元素间具有如下相关性：

其中，m＝1,2,…,N，1≤m±q≤N；表示矩阵K^C中位于反对角线且距离对角线距离为q的元素之间的相位差，表示矩阵K^C中位于(m–q,m+q)位置处的元素，表示矩阵K^C中位于(m,m)位置处的元素；

矩阵K^R中所有元素均含有随机参数，表现出随机性，因此矩阵元素之间不具有相关性。

3.根权利要求1所述的基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，其特征在于：所述步骤4中，信息矩阵K中散射波成分的滤除操作过程如下：

步骤41：根据下式将信息矩阵K逆时针旋转45°，并对相应元素补零，以保证矩阵维度不变，得到两个对应矩阵A1＝{a1_u,v}、A2＝{a2_u,v}：

若N为奇数：

若N为偶数：

其中，k_m,n表示信息矩阵K中位于(m,n)位置处的元素，矩阵A1的维度为N，矩阵A2的维度为N–1；将矩阵A1和矩阵A2统一用矩阵A表示，并用N_A表示矩阵A的维度；

步骤42：构造滤波矩阵对矩阵A进行滤波，滤波向量C的表达式为：

其中，为C的共轭转置；u＝1,2,…,N_A，

步骤43：对矩阵A进行滤波，将滤波后的矩阵A记为矩阵A^F：

A^F＝FA＝F(A^C+A^R)＝FA^C+FA^R

其中，矩阵A^C为矩阵A中不含随机参数的部分，矩阵A^R为矩阵A中含随机参数的部分；第一项FA^C＝A^C；

步骤44：将矩阵A^F顺时针旋转45°，得到滤波后的信息矩阵K，将滤波后的信息矩阵K记为矩阵K^F：

若(m-n)是奇数，

若(m-n)是偶数，

其中，表示矩阵K^F中位于(m,n)位置处的元素，表示矩阵A1^F中位于((m-n)/2+(N+3)/4,(m+n)/2)位置处的元素，表示矩阵A2^F中位于((m-n-1)/2+(N+3)/4,(m+n-1)/2)位置处的元素，矩阵A1^F为滤波后的矩阵A1，矩阵A2^F为滤波后的矩阵A2；

矩阵K^F的大小仍为N×N，与信息矩阵K的一致，滤波前后没有信息维度的缺失。

4.根权利要求1所述的基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，其特征在于：所述步骤5中，根据滤波后的信息矩阵进行时间反转成像过程如下：

步骤51：对矩阵K^F作奇异值分解：

其中，U^F为N×N阶酉矩阵，它的行向量是矩阵K^F的左奇异向量；Λ^F为半正定N×N阶对角矩阵；为N×N阶酉矩阵，为V^F的共轭转置，V^F的行向量是矩阵K^F的右奇异向量；Λ^F对角线上的元素即为矩阵K^F的奇异值；

步骤52：计算成像平面上任意点(X,Z)与超声换能器阵列之间的直达波传输矩阵G，矩阵G中位于(m,n)位置处的元素g_m,n为：

其中，r_m,n为超声换能器阵列中的第n个换能器单元与第m个换能器单元之间的距离；

步骤53：最大的M个奇异值所对应的特征向量与直达波传输矩阵G应用时间反转算子，第i个奇异值在深度Z处的光声信号强度为：

其中，V_i ^F表示V^F的第i行奇异向量；将最大的M个奇异值所对应求出的光声信号强度叠加，获得深度Z处的光声图像 是I^F向量的第m个元素，它表示成像平面上与第m个换能器单元所对应位置的初步成像值。

5.根权利要求1所述的基于全矩阵滤波和时间反转算子的光声断层成像的图像重构方法，其特征在于：所述步骤6中，修正参数生成与利用过程如下：

步骤61：在深度Z处，第m个换能器单元所对应的成像修正值为：

其中，为g_m,n的共轭；

步骤62：利用成像修正值对初步成像值进行修正，从而获得最终的成像值

即，最终的成像结果由滤波成像结果I^F与修正矩阵I^R对应元素相除得到。