本发明属于一种成像重建方法，特别涉及一种光声断层成像重建方法。
背景技术：
光声成像是一种新兴的成像技术，为临床成像提供了潜在的视角和机遇。光声断层成像是光声成像的一个重要分支，其优势在于成像深度较深，可达数厘米，符合临床成像的基本需求。目前已有基于环形阵列声学换能器(全角度探测)的光声断层成像用于乳腺的检查，然而环形声学探测器受其固有形状结构限制，难以用于与更多部位的光声成像和临床检查。基于手持式(有限角度探测)的声学换能器，其优势在于灵活便携，同时便于与目前临床上常用的超声成像结合，是光声断层成像的发展趋势。光声断层成像区别于其他光声成像技术的最显著特点在于光声断层成像需要将采集到的射频数据进行重建才能获得光声图像，目前常用的重建方法包括时间反演法、滤波反投影法，这两种方法对于全角度探测的光声断层成像具有较好效果；然而对于手持式的声学换能器，由于只采集了有限角度的信号，上述两种方法难以实现精确的重建成像。另一方面，二维矩阵声学换能器是一种新兴的手持式声学换能器，可以直接进行三维成像，目前还没有可以对二维矩阵声学换能器进行重建的方法。因此，目前亟需一种新型的光声重建算法，实现基于二维矩阵声学换能器的光声有限角度探测下的精确重建成像。
技术实现要素：
本发明的目的在于针对上述已有技术的不足，提供一种基于频域波数域的光声断层成像重建方法，该方法适用于二维矩阵声学换能器，可以在频率域和波数域内进行光声信号的三维重建，同时降低了有限角度下信号采集不全的影响，重建计算量较小，此外由于不是基于对时域信号的空间投影，大大减少了投影伪迹的影响。
为了实现上述目的，本发现的基于频域波数域的光声断层成像重建方法，包含以下步骤：
(1)、获取二维矩阵声学换能器采集到的射频数据r(x，y，z＝0，t)，以及二维矩阵声学换能器的系统参数和待测物体的声速分布v(z′)。其中二维矩阵声学换能器是基于线阵的超声换能器，射频数据r(x，y，z＝0，t)是实时采集的或事先保存好的，x，y是沿二维矩阵声学换能器的位置坐标，且x，y是相互垂直的，z是垂直二维矩阵声学换能器方向的位置坐标，t是表示时间顺序的坐标，二维矩阵声学换能器的系统参数为二维矩阵声学换能器的采样频率fs、二维矩阵声学换能器的阵元数量n、二维矩阵声学换能器的阵元间距pt。
(2)、对采集到的射频数据r(x，y，z＝0，t)进行时间增益补偿，得到补偿射频数据s(x，y，z＝0，t)。其中时间增益补偿的具体方法为将射频数据r(x，y，z＝0，t)的每一列上的数值乘以加权函数y(t)，y(t)＝at 1。即s(x，y，z＝0，t)＝r(x，y，z＝0，t)×(at 1)，1＜a＜5；由于a大于1，随着时间的增大，射频数据会被相应地方法，以抵消声学信号在传播过程中的损失。
(3)、对补偿射频数据s(x，y，z＝0，t)按x，y和t进行三维傅里叶变换，得到频域数据φ(kx，ky，z＝0，f)，再对φ((kx，ky，z＝0，f)进行插值，得到空间波数域数据φ′(kx，ky，kz，t＝0)，其中kx是沿x方向的波数，ky是沿y方向的波数，kz是沿z方向的波数，f是频率。由波动方程可知，z＝0平面(二维矩阵声学换能器所处位置)的频域数据φ(kx，ky，z＝0，f)与探测平面内的初始位置(t＝0)的空间波数域数据φ′(kx，ky，kz，t＝0)相差一个相位因子
(4)、对空间频域数据φ′(kx，ky，kz，t＝0)进行三维傅里叶逆变换，得到空间幅度数据s′(x，y，z，t＝0)。
(5)、对空间幅度数据s′(x，y，z，t＝0)进行解调取包络处理，再进行对比度调节，最终获得光声重建图像s*，其中解调取包络的方法为正交解调法或希尔伯特解调法，对比度调节方法，可以是对数压缩法，也可以是将解包络后得到的数据做n次幂，其中0＜n＜2。
(6)、对声速分布v(z′)进行优化，将优化得到的v2(z′)数据代入步骤(3)至步骤(5)中，得到最终组织对比度ctr最高的光声重建图像s2*；
其中声速分布v(z′)进行优化的具体步骤为：
(1)、设置优化次数计数器t＝0，设置最大优化次数t，t的取值范围为3-10，随机生成n个声速分布v1(z′)作为初始样本集合k1，n的取值范围为300-1000，其中每个声速分布中声速的范围随机在1350-1650m/s以内；
(2)、对上述步骤中的每个声速分布v1(z′)对光声图像进行重建，重建步骤为前文所述的光声断层成像重建方法中步骤(1)至步骤(5)所述的重建方法，得到n个光声重建图像si*，计算每个光声重建图像si*的ctr，ctr的具体计算方式为光声重建图像si*中目标物体与背景的灰度值之比；
(3)、对n个计算得到的ctr按数值大小进行降序排序，选择前30％较大ctr对应的声速分布构成样本集合ka，选择前30％-60％较大ctr对应的声速分布构成样本集合kb，ka和kb的样本数量均为0.3n；
(4)、从ka中随机选择一个声速分布样本va(z′)，从kb中随机选择一个声速分布样本vb(z′)，将两个声速分布样本的均值作为融合后的声速分布vc(z′)；
(5)、重复步骤(4)，得到0.2n个融合后的声速分布vc(z′)，将步骤(4)ka和kb中未被选择的声速分布与融合后的声速分布vc(z′)一起构成样本集合k2，k2的样本数量均为0.4n；
(6)、以k2作为初始样本集合重复步骤(2)-(6)，不断更新声速分布的样本集合；
(7)、当t＝t时，完成对声速分布v(z′)的优化，选择最终样本集合中ctr最大的样本作为最终的声速分布v2(z′)。
本发明与现有技术相比，具有如下有益效果：
对于基于手持式的二维矩阵声学换能器的光声断层重建，本发明提出的重建方法具有重建图像分辨率高、信噪比高、伪影少，重建计算速度快的优点。
附图说明
图1是本发明的流程图；
图2是二维矩阵声学换能器采集示意图；
图3是基于仿真的点声源射频数据重建结果；
图4是基于仿真的血管射频数据重建结果的二维截面；
图5是基于真实探测的血管模型射频数据重建结果的二维截面。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。
如图1所示，是本发明一种基于频域波数域的光声断层成像重建方法的流程图。
步骤(1)，获取二维矩阵声学换能器采集到的射频数据r(x，y，z＝0，t)，以及二维矩阵声学换能器的系统参数和待测物体的声速分布v(z′)。其中二维矩阵声学换能器是基于线阵的超声换能器，射频数据r(x，y，z＝0，t)是实时采集的或事先保存好的；如图2所示，是二维矩阵声学换能器采集示意图，x，y是沿二维矩阵声学换能器的位置坐标，且x，y是相互垂直的，z是垂直二维矩阵声学换能器方向的位置坐标，t是表示时间顺序的坐标；二维矩阵声学换能器的系统参数为二维矩阵声学换能器的采样频率fs、二维矩阵声学换能器的阵元数量n、二维矩阵声学换能器的阵元间距pt、声速分布v2(z′)。重建图像的宽度为pt×n，重建图像的高度由射频数据的长度、二维矩阵声学换能器的采样频率fs以及声速共同决定。
步骤(2)，对采集到的射频数据r(x，y，z＝0，t)进行时间增益补偿，得到补偿射频数据s(x，y，z＝0，t)。其中时间增益补偿的具体方法为将射频数据r(x，y，z＝0，t)的每一列上的数值乘以加权函数y(t)，y(t)＝2t 1。即s(x，y，z＝0，t)＝r(x，y，z＝0，t)×(2t 1)；随着时间的增大，射频数据会被相应地方法，以抵消声学信号在传播过程中的损失。
步骤(3)，再对补偿射频数据s(x，y，z＝0，t)按x和t进行二维傅里叶变换，得到频域数据φ(kx，ky，z＝0，f)，再对φ((kx，ky，z＝0，f)进行插值，得到空间波数域数据φ′(kx，ky，kz，t＝0)，其中kx是沿x方向的波数，ky是沿y方向的波数，kz是沿z方向的波数，f是频率。
由波动方程可知，z＝0平面(二维矩阵声学换能器所处位置)的频域数据φ(kx，ky，z＝0，f)与探测平面内的初始位置(t＝0)的空间波数域数据φ′(kx，ky，kz，t＝0)相差一个相位因子
此处需要强调的是频域数据φ(kx，ky，z＝0，f)与原始采集到的射频数据互为傅里叶变换对，而空间波数域数据φ′(kx，ky，kz，t＝0)与待重建的图像互为傅里叶变换对。因此建立了已知数据和待重建图像之间的联系。作为优选的，插值方法为加权反演插值：
φ′(kx，ky，kz，t＝0)＝∫φ(kx，ky，z＝0，f)m(kx，ky，kz，f)df，其中：
步骤(4)，再对空间频域数据φ′(kx，ky，kz，t＝0)进行二维傅里叶逆变换，得到空间幅度数据s′(x，y，z，t＝0)。
步骤(5)，采用希尔伯特解调法对空间幅度数据s′(x，y，z，t＝0)进行解调取包络处理，再进行对数压缩调节对比度，获得光声重建图像s*。
步骤(6)，对s*计算组织对比度ctr，ctr是图像中组织的图像亮度与图像中背景位置的图像亮度的比值，对声速分布v(z′)进行优化，将优化得到的v2(z′)数据代入步骤(3)至步骤(5)中，得到最终组织对比度ctr最高的光声重建图像s2*；
其中声速分布v(z′)进行优化的具体步骤为：
(1)、设置优化次数计数器t＝0，设置最大优化次数t，t的取值范围为3-10，随机生成n个声速分布v1(z′)作为初始样本集合k1，n的取值范围为300-1000，其中每个声速分布中声速的范围随机在1350-1650m/s以内；
(2)、对上述步骤中的每个声速分布v1(z′)对光声图像进行重建，重建步骤为前文所述的光声断层成像重建方法中步骤(1)至步骤(5)所述的重建方法，得到n个光声重建图像s1*，计算每个光声重建图像s1*的ctr，ctr的具体计算方式为光声重建图像s1*中目标物体与背景的灰度值之比；
(3)、对n个计算得到的ctr按数值大小进行降序排序，选择前30％较大ctr对应的声速分布构成样本集合ka，选择前30％-60％较大ctr对应的声速分布构成样本集合kb，ka和kb的样本数量均为0.3n；
(4)、从ka中随机选择一个声速分布样本va(z′)，从kb中随机选择一个声速分布样本vb(z′)，将两个声速分布样本的均值作为融合后的声速分布vc(z′)；
(5)、重复步骤(4)，得到0.2n个融合后的声速分布vc(z′)，将步骤(4)ka和kb中未被选择的声速分布与融合后的声速分布vc(z′)一起构成样本集合k2，k2的样本数量均为0.4n；
(6)、以k2作为初始样本集合重复步骤(2)-(6)，不断更新声速分布的样本集合；
(7)、当t＝t时，完成对声速分布v(z′)的优化，选择最终样本集合中ctr最大的样本作为最终的声速分布v2(z′)。
如图3所示，是基于仿真的点声源射频数据重建结果，原始图像是通过k-wave仿真软件生成的，其中采样频率fs＝30mhz、二维矩阵声学换能器的阵元数量n＝192、二维矩阵声学换能器的阵元间距pt＝0.2mm、声速v＝1540m/s。利用二维矩阵声学换能器采集到射频数据r(x，y，z＝0，t)，分别进行时间反演重建和本发明方法的重建；结果表明，本发明方法的重建结果的横向分辨率显著优于时间反演重建方法。
图4所示，是基于仿真的血管射频数据重建结果的二维截面，原始数据是通过k-wave仿真软件生成的，其中采样频率fs＝30mhz、二维矩阵声学换能器的阵元数量n＝192、二维矩阵声学换能器的阵元间距pt＝0.2mm、声速v＝z 1500m/s。利用二维矩阵声学换能器采集到射频数据r(x，y，z＝0，t)，分别进行滤波反投影重建和本发明方法的重建；结果表明，本发明方法的重建图像伪影远远小于滤波反投影重建，同时信噪比显著优于滤波反投影重建方法。
图5是基于真实探测的血管模型射频数据重建结果的二维截面，原始数据是利用血管模型实际采集得到的，其中采样频率fs＝50mhz、二维矩阵声学换能器的阵元数量n＝192、二维矩阵声学换能器的阵元间距pt＝0.25mm、声速v＝z 1500m/s。利用二维矩阵声学换能器采集到射频数据r(x，y，z＝0，t)，分别进行滤波反投影重建和本发明方法的重建；结果表明，本发明方法的重建图像伪影远远小于滤波反投影重建，同时具有较好的轴向分辨率。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。尽管本发明已进行了一定程度的描述，明显地，在不脱离本发明的精神和范围的条件下，可进行各个条件的适当变化。可以理解，本发明不限于所述实施方案，而归于权利要求的范围，其包括所述每个因素的等同替换。
1.一种基于频域波数域的光声断层成像重建方法，其特征在于，由以下步骤组成：
(1)、获取二维矩阵声学换能器采集到的射频数据r(x，y，z＝0，t)、二维矩阵声学换能器的系统参数和待测物体的声速分布v(z′)，其中的二维矩阵声学换能器采集到的射频数据r(x，y，z＝0，t)为实时采集或者提前采集后保存，x和y是沿二维矩阵声学换能器的位置坐标，x和y为相互垂直，z是垂直二维矩阵声学换能器方向的位置坐标，t是表示时间顺序的坐标；
(2)、对采集到的射频数据r(x，y，z＝0，t)进行时间增益补偿，得到补偿射频数据s(x，y，z＝0，t)，其中时间增益补偿方法为将射频数据r(x，y，z＝0，t)的每一列上的数值乘以加权函数y(t)，y(t)＝at 1，其中a的取值范围是1-5；
(3)、将补偿射频数据s(x，y，z＝0，t)沿着x方向，y方向和t方向进行三维傅里叶变换，得到φ(kx，ky，z＝0，f)，再对φ(kx，ky，z＝0，f)进行插值，得到空间波数域数据φ’(kx，ky，kz，t＝0)，其中空间波数域数φ’(kx，ky，kz，t＝0)＝∫φ(kx，ky，z＝0，f)m(kx，ky，kz，f)df，
(4)、对空间频域数据φ′(kx，ky，kz，t＝0)进行三维傅里叶逆变换，得到空间幅度数据s′(x，y，z，t＝0)；
(5)、对空间幅度数据s′(x，y，z，t＝0)进行解调取包络处理，再进行对比度调节，获得光声重建图像s*，其中解调取包络的方法为正交解调法或希尔伯特解调法，对比度调节方法为对数压缩法或者是将解包络后得到的数据做n次幂，其中0＜n＜2；
(6)、对声速分布v(z′)进行优化，将优化得到的v2(z′)数据代入步骤(3)至步骤(5)中，得到最终组织对比度ctr最高的光声重建图像s2*；
2.如权利要求1所述的一种基于频域波数域的光声断层成像重建方法，其特征在于上述步骤(6)中所述对声速分布v(z′)进行优化的具体步骤为：
(1)、设置优化次数计数器t＝0，设置最大优化次数t，t的取值范围为3-10，随机生成n个声速分布v1(z′)作为初始样本集合k1，n的取值范围为300-1000，其中每个声速分布中声速的范围随机在1350-1650m/s以内；
(2)、按照上述步骤的每个声速分布v1(z′)对光声图像进行重建，重建步骤为权利要求1中步骤(1)至步骤(5)所述的重建方法，得到n个光声重建图像s1*，计算每个光声重建图像s1*的ctr，ctr的具体计算方式为光声重建图像s1*中目标物体与背景的灰度值之比；
(3)、对n个计算得到的ctr按数值大小进行降序排序，选择前30％较大ctr对应的声速分布构成样本集合ka，选择前30％-60％较大ctr对应的声速分布构成样本集合kb，ka和kb的样本数量均为0.3n；
(4)、从ka中随机选择一个声速分布样本va(z′)，从kb中随机选择一个声速分布样本vb(z′)，将两个声速分布样本的均值作为融合后的声速分布vc(z′)；
(5)、重复步骤(4)，得到0.2n个融合后的声速分布vc(z′)，将步骤(4)中未被选择ka和kb的声速分布与融合后的声速分布vc(z′)一起构成样本集合k2，k2的样本数量均为0.4n；
(6)、以k2作为初始样本集合重复步骤(2)-(6)，不断更新声速分布的样本集合；
(7)、当t＝t时，完成对声速分布v(z′)的优化，选择最终样本集合中ctr最大的样本作为最终的声速分布v2(z′)。
技术总结