(19)中华人民共和国国家知识产权局
(12)发明专利申请
(10)申请公布号 CN 106037795 A(43)申请公布日 2016.10.26
(21)申请号 201610333423.4(22)申请日 2016.05.17
(71)申请人 南京大学
地址 210093 江苏省南京市鼓楼区汉口路
22号南京大学1019信箱(72)发明人 袁杰 朱昀浩 孙辉 尤琦
余双春 王学鼎 陶超 刘晓峻 (51)Int.Cl.
A61B 8/00(2006.01)
权利要求书2页 说明书6页 附图2页
CN 106037795 A(54)发明名称
一种优化全波反演法重建声速图像的方法(57)摘要
本发明公开了一种优化全波反演法重建声速图像的方法，包括以下步骤：利用超声数据采集设备，得到传感器所在位置的实测声压信号；在计算机上完全模拟发射、接收过程作为正向函数，该函数输入变量为声速图矩阵，输出变量为传感器对应位置所记录的声压信号与的实测声压信号的差的范数；进入迭代过程：给出一个声速图像作为初始图像，使其作为自变量带入到上述正向函数中，得到输出范数；计算出正向函数在当前声速图位置的一阶梯度图；通过拟牛顿法解决最优化问题，更新原有声速图，得到新的声速图；将新的声速图作为初始图像重复迭代过程，直到得到满足条件退出迭代过程，并得到当前的声速图像。
CN 106037795 A
权 利 要 求 书
1/2页
1.一种优化全波反演法重建声速图像的方法，其特征在于，包括如下步骤：步骤一，利用超声数据采集设备，依次完成M个信号源单独发射、所有的M个传感器接收的过程，直到所有信号源都发射过一次，得到传感器所在位置的实测声压信号gm(r，t)，其中r为传感器坐标向量参数，t为时间参数，m表示发射次序，m＝1，2，...，M，M同时也表示总共的发射次数；
步骤二，根据步骤一所用设备在计算机上完全模拟步骤一中的发射、接收过程作为正向函数，使之输入变量为声速图矩阵c，输出变量为传感器对应位置所记录的声压信号gm(r，t)与步骤一中的实测声压信号gm(r，t)的差的范数L；
步骤三，开始迭代过程，首次迭代中给出声速图像c0作为初始图像，也即c＝c0，作为步骤二中的正向函数的输入量，得到输出变量L；
步骤四，计算出正向函数在声速图c位置的一阶梯度图J；步骤五，通过拟牛顿法解决最优化问题，利用步骤四中的梯度图J更新声速图c，得到新的声速图；
步骤六，将新的声速图作为初始图像重复迭代过程，重复步骤三到六，直到得到满足条件退出迭代过程；
步骤七，得到满足步骤六的条件时对应的声速图像c即为优化的全波反演法重建的声速图像。
2.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤一中发射接收过程遵循以下原则：M个传感器和M个信号源成对地位于同一位置，每次只有一个信号源发射，所有传感器接收，传感器共M个，那么发射次数也应为M次。
3.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤二中的正向函数模拟了计算区域、传播介质、发射设备、接收设备以及声波传播过程。
4.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤二中的正向函数，也即模拟过程中，接收到的声压信号与波源、传播介质满足方程：
其中r为坐标向量参数，t为时间变量参数，pm(r，t)为声压，m表示发射的顺序次数，m＝1，2，...，M，M表示总共的发射次数，sm(r，t)为超声信号，c(r)为声速图矩阵，c(r)有n＝N×N×Nz个格点，时间变量参数t在计算中需要离散化表示；在传感器点记录声压信息的采样时间间隔以dt表示，T表示时间总长度。
5.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤二中范数可以表示为：
其中gm，gm分表表示步骤一得到的实测声压和在接下来的步骤三得到的模拟声压；与pm(r，t)区别是，gm表示在测量位置的声压信号，也就是r为传感器位置时的pm(r，t)。
6.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤三中得到L，求解声速图的问题就转化为在接下来的步骤中解决如下表达式表述的问题：
2
CN 106037795 A
权 利 要 求 书
2/2页
该式表达了一个非线性的最优化问题，即为最后要求得的矩阵。
7.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤三中除了得到输出值L，函数产生的声压信息pm(r，t)也在此步骤中得到并保留，并被用于接下来步骤中的运算。
8.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤四中c指的是当前声速图，如果为首次迭代，那么c＝c0。
9.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤四中的梯度的计算可以用伴随状态法计算。
10.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步
BFGS法均可以在此应用，步长根据线性搜索法确定，具体骤五的拟牛顿法两种算法DFP法、
而言线性搜索法包括了进退法和黄金分割法。
11.根据权利要求1所述的一种优化全波反演法重建声速图像的方法，其特征在于，步骤六中阈值的设定是相对的，在具体实验分析中而定的，除此以外，还可以设定固定的迭代次数，利用迭代的过程不断更新图像的方法重建声速图像。
3
CN 106037795 A
说 明 书
一种优化全波反演法重建声速图像的方法
1/6页
技术领域
[0001]本发明属于医学超声成像及处理领域，特别是一种针对声速图像重建的优化方法。
背景技术
[0002]经过几十年的研究，全波反演重建技术在医学领域的应用也日趋深入。例如，在广泛应用于乳腺癌检测的超声检测层析成像的重建过程中，全波反演的重建方法得到的图像具有极具优化的空间分辨率，这一点也强过传统的Ray-Based重建技术，Ray-Based将声波传播假设为具有光的折射特性，所以不能重建出全波形所具有的特征。然而由于全波反演方法本身需要相当长的计算时间和巨大的计算机存储资源，这项技术没有被广泛的应用于医学临床应用。因此，需要一种优化的方法来提高全波反演重建声速图像方法的效率。发明内容
[0003]发明目的：本发明所要解决的技术问题是针对现有全波反演重建声速图像速度慢，效率低的不足，提供了一种从算法的角度提高重建速度和效率，并保证图像质量的优化方法。
[0004]为了解决上述技术问题，本发明公开了一种利用拟牛顿法降低运算量并提高效率的全波反演重建声速图像方法，包括如下步骤：[0005]步骤一，依次完成M个信号源单独发射、所有的M个传感器利用超声数据采集设备，接收的过程，直到所有信号源都发射过一次，得到传感器所在位置的实测声压信号gm(r，t)，其中r为传感器坐标向量参数，t为时间参数，m表示发射次序，m＝1，2，...，M，M同时也表示总共的发射次数；[0006]步骤二，根据步骤一所用设备在计算机上完全模拟步骤一中的发射、接收过程作为正向函数，使之输入变量为声速图矩阵c，输出变量为传感器对应位置所记录的声压信号gm(r，t)与步骤一中的实测声压信号gm(r，t)的差的范数L；[0007]步骤三，开始迭代过程，首次迭代中给出声速图像c0作为初始图像，也即c＝c0，作为步骤二中的正向函数的输入量，得到输出变量L；[0008]步骤四，计算出正向函数在声速图c位置的一阶梯度图J；[0009]步骤五，通过拟牛顿法解决最优化问题，利用步骤四中的梯度图J更新声速图c，得到新的声速图；[0010]步骤六，将新的声速图作为初始图像重复迭代过程，重复步骤三到六，直到得到满足条件退出迭代过程；[0011]步骤七，得到满足步骤六的条件时对应的声速图像c即为优化的全波反演法重建的声速图像。
[0012]本发明中，优选地，步骤一中所用信号源用于超声波的发射，传感器用做声压信号也即gm的记录、保存，其中下标m表示发射次序所有的信号源和传感器成对地位于同一位
4
CN 106037795 A
说 明 书
2/6页
置。
本发明中，优选地，步骤一中发射接收过程遵循以下原则：每次只有一个信号源发
射，其余传感器(包括其本身)接收，传感器共M个，那么发射次数也应为M次。[0014]本发明中，优选地，步骤一中选用环形传感器围绕目标，或双排线型传感器置于目标两端，共有M个传感器依次排列，其位置坐标可由r对应取值唯一表示。[0015]本发明中，优选地，步骤二中的正向函数模拟了计算区域、传播介质、发射设备、接收设备以及声波传播过程。[0016]具体而言，正向函数具体包括定义了计算区域、传播介质、发射设备、接收设备。计算区域的定义是传播介质的定义的基础，二者可以表示为大小一致的矩阵。定义计算区域大小n＝N×N×Nz，表示将传播介质划分为n个区域，我们称其每个区域为一个格点，其中N和Nz分别表示长或宽和厚度。传播介质在其基础上具体化了每个格点的参数，包括但不限于声速、密度、衰减系数等信息。发射设备、接收设备的两种功能由信号源和传感器来完成，信号源和传感器位置相同并位于计算区域中，计算区域的边界经过处理后，声波会在边界被完全吸收。
[0017]本发明中，优选地，步骤二中的声速图为一个矩阵，也即一个数字化的传播介质模型，声速只是传播介质的一部分信息，其余信息待定为一般参数，一般是目标周围介质的对
该声速图矩阵每一点的值为对应位置的声速值。应参数，传播介质中应包含目标信息，
[0018]本发明中，优选地，步骤二中的正向函数为模拟过程，正向函数输入变量为声速图矩阵，该函数可以通过给定的声速图直接得到模拟声压信息，进而可得范数L。[0019]本发明中，优选地，步骤二中的正向函数，也即模拟的声波传播过程中，接收的声压信号与波源、传播介质满足方程：
[0020][0021][0013]
其中r为坐标向量参数，t为时间变量参数，pm(r，t)为声压，m表示发射的顺序次数，m＝1，2，...，M，M表示总共的发射次数，sm(r，t)为超声信号，c(r)为声速图矩阵，c(r)有n＝N×N×Nz个格点。
由于奈奎斯特定律，上述的超声信号sm(r，t)的频率不能超过
其
[0022]
中c0在表示声速(m/s)，dx代表每一格点换算成实际的边长(m)。
[0023]本发明中，优选地，步骤二中的正向函数，时间变量参数t在计算中需要离散化表示。在传感器点记录声压信息的采样时间间隔以dt(s)表示，T(s)表示时间总长度。[0024]本发明中，优选地，步骤二中的正向函数，也即模拟过程中，r可以为任意坐标，不限于传感器的位置。故而，在模拟过程中，不仅仅可以得到传感器位置的声压信息，包括可以得到整个计算区域所有位置的声压信息。[0025]本发明中，优选地，步骤二中范数可以表示为：
[0026][0027]
其中gm，gm分表表示步骤一得到的实测声压和在接下来的步骤三得到的模拟声压。与pm(r，t)区别是，gm表示在测量位置的声压信号，也就是r等于传感器位置坐标时的pm(r，
5
CN 106037795 A
说 明 书
3/6页
t)。
[0028]
本发明中，优选地，步骤三中得到L，求解声速图的问题就转化为在接下来的步骤
中解决如下表达式表述的问题：
[0029][0030][0031]
该式表达了一个非线性的最优化问题，即为最后要求得的矩阵。
本发明中，优选地，步骤三中除了得到输出值L，函数产生的声压信息pm(r，t)也在
此步骤中得到并保留，并被用于接下来步骤中的运算。[0032]本发明中，优选地，步骤四中c指的是当前声速图，如果为首次迭代，那么c＝c0，否则应按上次迭代结果带入。[0033]本发明中，优选地，步骤四中的梯度的计算可以用伴随状态法(adjoint state method)计算并说明如下。首先，L计算式表述如下：
[0034][0035]
其中Ωm(N×N×Nz)表示整个包括目标在内的计算区域，T表示时间总长度。那么根据伴随状态法和L的表述式，L的梯度计算式表述如下：
[0036][0037][0038][0039]
其中qm(r，t)是下面波动方程的解。
其中τt)可以表述为τt)＝gm(r，T-t)-gm(r，T-t)。那么计算中L关于c的梯m(r，m(r，度用J表达如下：
[0040][0041]
其中n表示对应的声速图像中每一个离散化格点，同时qm，pm对应qm(r，t)和pm(r，
t)中坐标参数r为对应格点位置时的值。[0042]本发明中，优选地，步骤五的拟牛顿法两种算法DFP法、BFGS法均可以在此应用，在本问题中具体表述如下：
[0043]拟牛顿法用近似的Hesse矩阵Bk代替真实的Hesse矩阵。那么Bk在每次迭代中的更新公式如下，在DFP法中记
[0044][0045][0046][0047][0048]
其中，c(k)表示第k次迭代中的声速图，J(k)表示第k次迭代中的梯度图。
在BFGS法中，Bk的计算公式有所不同并区别如下：
那么迭代公式，也即最终得到的c(k+1)应为：
6
CN 106037795 A[0049][0050]
说 明 书
4/6页
c(k+1)＝c(k)+akpk其中
ak为步长根据线性搜索法确定。具体而言线性搜索法
包括了进退法和黄金分割法。
[0051]本发明中，优选地，步骤六中利用迭代的过程不断更新图像的方法重建声速图像。[0052]本发明中，优选地，步骤六中阈值的设定是相对的，在具体实验分析中而定的，除此以外，还可以设定固定的迭代次数。[0053]本发明中，针对全波反演法现有的低效问题，进行了一定的优化，由于拟牛顿法充分利用一阶梯度而得到近似的二阶梯度矩阵大大提高了效率，又在GPU的应用下，计算时间可以从一周左右缩短至几个小时。
附图说明
[0054]下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的有关叙述和其他方面的优点将会变得更加清楚。[0055]图1为本发明流程图。
[0056]图2为本发明系统截面示意图。
具体实施方式
[0057]如图1所示，本发明公开了一种优化全波反演法重建声速图像的方法，包括以下步骤：
[0058]步骤一，依次完成M个信号源单独发射、所有的M个传感器利用超声数据采集设备，接收的过程，直到所有信号源都发射过一次，得到传感器所在位置的实测声压信号gm(r，t)，其中r为传感器坐标向量参数，t为时间参数，m表示发射次序，m＝1，2，...，M，M同时也表示总共的发射次数；[0059]步骤二，根据步骤一所用设备在计算机上完全模拟步骤一中的发射、接收过程作为正向函数，使之输入变量为声速图矩阵c，输出变量为传感器对应位置所记录的声压信号gm(r，t)与步骤一中的实测声压信号gm(r，t)的差的范数L；[0060]步骤三，开始迭代过程，首次迭代中给出声速图像c0作为初始图像，也即c＝c0，作为步骤二中的正向函数的输入量，得到输出变量L；[0061]步骤四，计算出正向函数在声速图c位置的一阶梯度图J；[0062]步骤五，通过拟牛顿法解决最优化问题，利用步骤四中的梯度图J更新声速图c，得到新的声速图；[0063]步骤六，将新的声速图作为初始图像重复迭代过程，重复步骤三到六，直到得到满足条件退出迭代过程；[00]步骤七，得到满足步骤六的条件时对应的声速图像c即为优化的全波反演法重建的声速图像。
[0065]如图2所示，实验中传感器共M个在计算区域内环形排列。目标中心尽量置于计算区域的中心位置。图2所示意的结构既是实验设备结构，也是重建过程中的正向函数的计算区域及其相关信息的示意图。[0066]本实例中，步骤一中选用了既可用于超声波的发射，也可用做声压信息(也即指得
7
CN 106037795 A
说 明 书
5/6页
到声压信号gm)的记录、保存传感器。发射接收过程遵循以下原则：每次只有一个传感器发射，其余传感器(包括其本身)接收，接收传感器为M个，发射次数也为M次。[0067]本实例中，步骤一中选用环形传感器围绕目标，其位置坐标可由r唯一表示。[0068]本实例中，步骤二中的初始声速图选为一个矩阵，由普通的Ray-based方法得到。该声速图矩阵每一点的值为对应位置的声速值。[0069]在本实例的首次迭代中，正向函数输入变量为初始声速图矩阵c0。[0070]本实例中，由于奈奎斯特定律，超声信号sm(r，t)的频率不能超过
其中c0在表示声速(m/s)，dx代表每一格点换算成实际的边长(m)。故在本实例
中，发射超声信号设定为频率f＝0.8Mhz的正弦函数与高斯核函数的乘积，c0＝1510m/s并满足上述要求。格点长度大小设为波长的三分之一。具体实际的计算区域设定为N＝256，Nz＝5，N对应0.128m实际长度。那么dx＝0.5mm[0071]本实例中，范数表示为：
[0072]
其中gm，gm分表表示步骤一得到的实测声压和在接下来的步骤三得到的模拟声压。[0074]本实例中，步骤四中的梯度的计算可以用伴随状态法(adjoint state method)计算，首次迭代中c＝c0，否则根据迭代确定c的取值，L计算式表述如下：
[0075][0076]
[0073]
其中Ωm(N×N×Nz)表示整个包括目标在内的计算区域，T表示时间总长度。T设定
根据伴随状态法和L的表述
为足够超声波传播的实践长度。dt的选取设定为式，L的梯度计算式表述如下：
[0077][0078][0079][0080]
其中qm(r，t)是下面波动方程的解。
其中τt)可以表述为τt)＝gm(r，T-t)-gm(r，T-t)。那么计算中L关于c的近m(r，m(r，似梯度用J表达如下：
[0081]
其中n表示对应的声速图像中每一个离散化格点，同时qm，pm对应qm(r，t)和pm(r，
t)中坐标参数r为对应格点位置时的值。[0083]本实例中，步骤五的BFGS拟牛顿法在本问题中具体表述如下：[0084]在BFGS法中，Bk的计算公式如下：
[0082]
8
CN 106037795 A[0085][0086][0087][0088]
说 明 书
6/6页
那么迭代公式，也即最终得到的c(k+1)应为：c(k+1)＝c(k)+akpk其中
ak为步长根据线性搜索法确定。具体而言线性搜索法包括了进
退法和黄金分割法两种方法共同应用，进退法求出区间范围[a，b]，黄金分割法在区间[a，
b]之间得出合适的步长ak。[00]本实例中，步骤六中利用迭代的过程不断更新图像的方法重建声速图像。[0090]本实例中，步骤六中阈值的设定为1×10-9，每次迭代后保存结果。迭代过程中可以随时停止并查看结果，并且可以继续上次的结果运行。[0091]本实例流程图参照图1。
[0092]本发明提出了一种优化全波反演法重建声速图像的方法，应当指出，步骤一种涉及的实验设备型号形式不对本专利构成；传感器位置可以是环形也可以是其他形式，具体使用何种分布形式不对本专利构成；模拟过程中的计算区域设定的大小等非关键参数，不对本专利构成。应当指出，对于本技术领域的普通人员来说，在不脱离发明原理的前提下还可以做出若干改进和润饰，这些也应视为本发明的保护范围。另外，本实例中未明确的各组成部分均可用现有技术加以实现。
9
CN 106037795 A
说 明 书 附 图
1/2页
图1
10
CN 106037795 A
说 明 书 附 图
2/2页
图2
11
