本发明提供一种基于波动方程的医学超声CT多参数图像重建方法，包括：利用超声采集装备发射、接收超声，获取真实三维空间中接收阵元的观测声压信号；利用多参数黏声波动方程仿真模拟上述过程，获取对应于真实三维空间中接收阵元的合成声压信号；利用反卷积目标函数和二范数目标函数比较观测声压信号和合成声压信号的差异，并求解一阶梯度；利用一阶梯度求取新的参数矩阵；将新的参数矩阵输入到多参数黏声波动方程进行迭代运算，获取下一个新的参数矩阵得到符合条件的多参数重建图像。本发明采用FWI算法，引入多参数黏声波动方程和反卷积目标函数，利用全波形信息精准重构多种人体组织参数的图像，获得优于传统射线方法的高分辨率医学超声图像。
1.一种基于波动方程的医学超声CT多参数图像重建方法，其特征在于，所述方法包括：利用超声CT数据采集装备发射、接收超声，获取真实三维空间数据；所述数据是所述真实三维空间中接收阵元所在位置的观测声压信号；
利用多参数黏声波动方程仿真模拟所述数据采集装备发射、接收超声的过程，获取对应所述真实三维空间中接收阵元所在位置的合成声压信号；所述多参数黏声波动方程的输入变量为原始声速图、原始声阻抗图和原始衰减品质因子图所组成的参数矩阵；所述多参数黏声波动方程，即为正演函数，给定发射阵元，即位置xs的信息，输入声速图、声阻抗图、声衰减图的参数矩阵，通过数值仿真方法，包括有限差分、有限元方法，离散地求解在整个计算空间x和时间t上的声压信息p(x，t)，然后提取出接收阵元xr处合成声压信号dcal(xr，t；
xs)＝R(x，xr)p(x，t；xs)，最大记录时间T应记录相距最远的发射和接收阵元组合产生的穿过目标人体组织的透射波信号；
利用反卷积目标函数和二范数目标函数求解所述观测声压信号和所述合成声压信号的差异；所述反卷积目标函数，即为反演函数，使得全波形反演(Full waveform inversion，FWI)算法的结果更接近全局最小值，求解实测数据与仿真数据之间的时间域维纳滤波器，通过对该滤波器加权设计目标函数，有效突出两种数据间基于波动方程的渡越时间差异，弱化正演振幅不准确对反演的影响；
将所述反卷积目标函数和所述二范数目标函数通过伴随状态法计算得到关于所述参数矩阵的一阶梯度；基于所述一阶梯度和所述参数矩阵计算得到新的参数矩阵；将所述新的参数矩阵转换为新声速图、新声阻抗图和新衰减品质因子图；
将所述新的参数矩阵输入到所述多参数黏声波动方程，结合所述反卷积目标函数、所述一阶梯度进行迭代运算，获取下一个新的参数矩阵；所述迭代运算停止条件为所述合成声压信号与所述观测声压信号达到一致。
2.根据权利要求1所述的多参数图像重建方法，其特征在于，所述超声CT数据采集装备具备多层环形阵列，单次激发后多层阵列同时采集数据，获取具有纵向孔径的三维容积数据；
所述数据采集装备为含有多个阵元的观测阵列换能器；所述多个阵元的每一个阵元能发射、接收超声信号；
所述观测阵列换能器选用环形探头、对穿线阵探头、对穿面阵探头和对穿凹阵探头。
3.根据权利要求2所述的多参数图像重建方法，其特征在于，所述多个阵元的每一个阵元所在位置的观测声压信号为dobs(xr，t；xs)，其中，xs为发射阵元的坐标位置，xr为接收阵元的坐标位置，t表示接收时间。
4.根据权利要求1所述的多参数图像重建方法，其特征在于，所述多参数黏声波动方程，其三维方程公式如下：其中，描述声场的质子速度为vi，应力为σij，合成声压p＝∑kσkk/3；描述介质物理性质的参数密度为ρ，声速为Vp，声阻抗Ip＝ρVp描述，衰减系数由品质因子Q通过滞弹系数Yl≈‑1ylQ 引入，该滞弹系数采用了L个在参考频率ωl下的记忆变量 如不考虑衰减，则最后两个方程及与Yl相关的系数都会被舍去，退化回传统的各向同性声波方程，其二维形式通过舍去与z相关项退化形成，变量vi，σij，ρ，Vp，Yl，均为定义在计算成像空间x和时间t上的图像。
5.根据权利要求1所述的多参数图像重建方法，其特征在于，所述反卷积目标函数表示为：其中，A(τ)为权重函数，W表示维纳滤波器，
在反演时，权函数A(τ)的选择至关重要，权函数决定算法的鲁棒性，|τ|/τmax、τ和高斯函数作为选择；多尺度反演的特点通过设置权重函数中的τ值上限来实现，因为该上限表征目标函数考虑的观测数据和仿真数据的最大渡越时间差；该上限应结合分频多尺度反演选取，在相应频段数据主频倒数的3‑5倍之间，随着频率逐渐升高，τ上限逐渐减小，再利用传统二范数目标函数Cl2迭代提高分辨率，
6.根据权利要求1所述的多参数图像重建方法，其特征在于，所述反卷积目标函数关于所述参数矩阵所述一阶梯度 通过伴随状态法计算，其中，伴随态波场算子μ为方程 的解；改变反卷积目标函数不影响
正传波场w的计算；但是，伴随态波场μ是通过反传 获取而偏导数 会随
反卷积目标函数改变而变化；在二范数目标函数情况下，所述偏导数为：
7.根据权利要求1所述的多参数图像重建方法，其特征在于，所述迭代视为求解求优化问题的过程，在第0次迭代输入初始的所述参数矩阵m0，即基于先验知识的对声参数模型的估计；
在检测乳腺软组织的情况下，假设起始所述参数矩阵m0等效于均匀的水介质。
8.根据权利要求1所述的多参数图像重建方法，其特征在于，所述迭代更新采用梯度下降法，沿着下降方向d迭代更新所述参数矩阵的值，mk+1＝mk+αkdk，
其中，α利用线搜索方法获取，而d 在不同优化策略，包括最速下降法、共轭梯度法、L‑BFGS法、牛顿法中有不同的计算方法，最速下降法利用目标函数关于模型的梯度的负方向作为下降方向，即 共轭梯度、L‑BFGS方法需要利用前置迭代中的梯度信息改善当前梯度获取下降方向，牛顿法需估计海瑟矩阵，即目标函数关于梯度的二阶导数，来约束当前梯度获取下降方向。
9.根据权利要求1所述的多参数图像重建方法，其特征在于，所述迭代运算停止条件选用所述反卷积目标函数值下降的阈值或者迭代的次数。
10.根据权利要求1所述的多参数图像重建方法，其特征在于，求解所述多参数黏声波动方程、所述反卷积目标函数的代码利用CPU进行验证，然后移植到GPU上提高效率以应对时间敏感的超声CT医疗成像临床诊断任务需求。
一种基于波动方程的医学超声CT多参数图像重建方法
技术领域
[0001] 本申请涉及医学超声成像技术领域，尤其涉及一种基于波动方程的医学超声CT多参数图像重建方法。
背景技术
[0002] 传统的超声CT图像重建一般分为基于回波的回声强度图像重建和基于透射波的声速、衰减系数重建。回声强度图像分辨率较高，但仅可在组织明显病变、已形成具有声属性差异的病灶时辅助医生诊断。透射波方法可重建声速、衰减系数等功能参数，对于监测病变组织早期的功能参数变化具有重要意义，但是分辨率较低，极端情况下难以区分背景组织与病变组织的特征差异。因此，结合多类型的声波信号，包括但不限于回波和透射波，重建高分辨率的功能参数，对于病变的早期成像和诊断具有重要意义。
[0003] 多类型信号的超声CT图像重建有基于射线理论和波动理论的两种方法。基于射线理论的多波型重建一般采用对前述两种图像进行融合的方法，受限于两者分辨率的差距，融合过程需要适当的处理，否则可能引入错误的病灶信息。基于波动方程的重建可以直接利用全部的波形数据重建声速、声阻抗、衰减系数等多种参数图像，并具有较高分辨率，辅助区分病灶效果较好。全波形反演(Full waveform inversion，FWI)的成像方法，结合前沿计算设备其计算量在可接受范围内，恰当地设计处理流程也可保证其方法的鲁棒性，目前来看是一种稳定、实用、高分辨率的多特征图像重建方法。
[0004] 由于基于FWI的多参数图像重建方法需要引入复杂的多参数波动方程(如声速、声阻抗、衰减等)，因此现在国内外基于FWI的超声CT图像重建方法主要关注声速单参数重建。此外，FWl由于是一种局部最优化方法，存在局部最小值的问题，影响结果准确性。因此，需要一种基于多参数波动方程的精准FWI超声CT图像重建方法，来获取高分辨率的多参数图像辅助病灶判别和诊断。
发明内容
[0005] 本申请的目的在于解决基于FWI获取高分辨率的多参数超声CT精准图像的问题。
[0006] 为实现上述目的，本申请提供了一种基于波动方程的医学超声CT多参数图像重建方法，所述包括以下步骤：
[0007] 利用超声CT数据采集装备发射、接收超声，获取真实三维空间数据；所述数据是所述真实三维空间中接收阵元所在位置的观测声压信号；
[0008] 利用多参数黏声波动方程仿真模拟所述采集设备发射、接收超声的过程，获取对应所述真实三维空间中接收阵元所在位置的合成声压信号；所述多参数黏声波动方程的输入变量为原始声速图、原始声阻抗图和原始衰减品质因子图所组成的参数矩阵；
[0009] 利用反卷积目标函数和二范数目标函数求解所述观测声压信号和所述合成声压信号的差异；
[0010] 将所述反卷积目标函数和所述二范数目标函数通过伴随状态法计算得到关于所述参数矩阵的一阶梯度；基于所述一阶梯度和所述参数矩阵计算得到新的参数矩阵；将所述新的参数矩阵转换为新声速图、新声阻抗图和新衰减品质因子图；
[0011] 将所述新的参数矩阵输入到所述多参数黏声波动方程，结合所述反卷积目标函数、所述一阶梯度进行迭代运算，获取下一个新的参数矩阵；所述迭代运算停止条件为所述合成声压信号与所述观测声压信号达到一致。
[0012] 优选的，所述超声CT数据采集装备具备多层环形阵列，单次激发后多层阵列同时采集数据，获取具有纵向孔径的三维容积数据；
[0013] 所述采集数据装备为含有多个阵元的观测阵列换能器；所述多个阵元的每一个阵元能发射、接收超声信号；
[0014] 可选的，所述观测阵列换能器选用环形探头、对穿线阵探头、对穿面阵探头和对穿凹阵探头。
[0015] 优选的，所述多个阵元的每一个阵元所在位置的观测声压信号为dobs(xr，t；xs)，其中，xs为发射阵元的坐标位置，xr为接收阵元的坐标位置，t表示接收时间。
[0016] 优选的，所述多参数黏声波动方程，即为正演函数，其三维方程公式如下：
[0017]
[0018] 其中，描述声场的质子速度为vi，应力为σij，合成声压p＝∑kσkk/3；描述介质物理性质的参数密度为ρ，声速为Vp，声阻抗Ip＝ρVp描述，衰减系数由品质因子Q通过滞弹系数引入，该滞弹系数采用了L个在参考频率 下的记忆变量 如不考虑衰减，则最后两个方程及与 相关的系数都会被舍去，退化回传统的各向同性声波方程，其二维形式可通过舍去与z相关项退化形成，变量vi，σij，ρ，Vp， 均为定义在计算成像空间x和时间t上的图像；
[0019] 给定发射阵元(位置xs)的信息，输入声速图、声阻抗图、声衰减图的参数矩阵，可以通过数值仿真方法(如有限差分、有限元等方法)离散地求解在整个计算空间x和时间t上的声压信息p(x，t)，然后提取出接收阵元xr处合成声压信号dcal(xr，t；xs)＝R(x，xr)p(x，t；xs)，最大记录时间T应可记录相距最远的发射和接收阵元组合产生的穿过目标人体组织的透射波信号。
[0020] 优选的，所述反卷积目标函数，即为反演函数，使得FWI算法的结果更接近全局最小值，求解实测数据与仿真数据之间的时间域维纳滤波器，通过对该滤波器加权设计目标函数，可有效突出两种数据间基于波动方程的渡越时间差异，弱化正演振幅不准确对反演的影响，所述反卷积目标函数可表示为：
[0021]
[0022] 其中，A(τ)为权重函数，w表示维纳滤波器，
[0023]
[0024] 在反演时，权函数A(τ)的选择至关重要，权函数可以决定算法的鲁棒性，|τ|/τmax、2
τ和高斯函数均可作为选择。多尺度反演的特点可以通过设置权重函数种中的τ值上限来实现，因为该上限表征目标函数可考虑的观测数据和仿真数据的最大渡越时间差。该上限应结合分频多尺度反演选取，在相应频段数据主频倒数的3‑5倍之间，由于频率逐渐升高因此也实现了τ上限逐渐减小的效果，再利用传统二范数目标函数Cl2迭代提高分辨率，[0025]
[0026] 所述传统二范数目标函数在稳定性、准确性和高效性等方面还有待改进：在数据中缺少低频信息、初始模型离真实模型较远的情况下，优化算法难以有效收敛到全局最小值获取精准的超声CT图像重建。
[0027] 优选的，所述反卷积目标函数关于所述参数矩阵所述一阶梯度 可通过伴随状态法计算，
[0028]
[0029] 其中，伴随态波场算子μ为方程 的解。改变反卷积目标函数不影响正传波场w的计算。但是，伴随态波场μ是通过反传 获取，而该偏导数会随反卷积目标函数改变而变化。在二范数目标函数情况下，该偏导数为：
[0030]
[0031] 在反卷积条件下，该偏导数为：
[0032]
[0033] 将 通过波动方程反向传播即可计算伴随态波场算子μ。
[0034] 优选的，所述迭代可视为求解求优化问题的过程，在第0次迭代输入初始的所述参数矩阵m0，即基于先验知识的对声参数模型的估计；
[0035] 可选的，在检测乳腺等软组织的情况下，可以假设起始所述参数矩阵m0等效于均匀的水介质。
[0036] 优选的，所述迭代更新一般采用梯度下降法，沿着下降方向d迭代更新所述参数矩阵的值，
[0037] mk+1＝mk+αkdk，
[0038] 其中，αk一般利用线搜索方法获取，而dk在不同优化策略(如：最速下降法、共轭梯度法、L‑BFGS法、牛顿法)中有不同的计算方法，最速下降法一般利用目标函数关于模型的梯度的负方向作为下降方向，即 共轭梯度、L‑BFGS等方法需要利用前置迭代中的梯度信息改善当前梯度获取下降方向，牛顿法需估计海瑟矩阵，即目标函数关于梯度的二阶导数，来约束当前梯度获取下降方向。
[0039] 优选的，所述迭代运算停止条件可以选用所述反卷积目标函数值下降的阈值或者迭代的次数。
[0040] 优选的，求解所述多参数黏声波动方程、所述反卷积目标函数的代码可利用CPU进行验证，然后移植到GPU上提高效率以应对时间敏感的超声CT医疗成像临床诊断任务需求。
[0041] 本申请提供一种基于波动方程的医学超声CT高分辨率图像重建方法，针对精准多参数图像重建的问题引入黏声波动方程(含声速、声阻抗、衰减参数)和反卷积目标函数，提供了一种从算法的角度利用全波形信息精准重构多种人体组织参数图像的方法，并保证显著优于传统射线方法的高分辨率图像质量，为病灶诊断提供影像学和功能参数辅助诊断依据，包括但不限于应用在筛查质密乳腺内微小肿瘤这类传统标准钼靶法难以早期诊断的场景。
附图说明
[0042] 为了更简单说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0043] 图1为本申请实施例中基于波动方程的医学超声CT多参数图像重建方法的流程示意图；
[0044] 图2为本申请一个实施例中基于FWI的超声CT图像重建流程；
[0045] 图3为本申请一个实施例中人体组织仿体和基于FWI的超声CT多参数图像。
具体实施方式
[0046] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例，本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围[0047] 图1中，本发明提出的一种基于波动方程的医学超声CT多参数图像重建方法的步骤如下：
[0048] 步骤S110：利用超声CT数据采集装备发射、接收超声，获取真实三维空间数据；数据是真实三维空间中接收阵元所在位置的观测声压信号；
[0049] 步骤S120：利用多参数黏声波动方程仿真模拟采集设备发射、接收超声的过程，获取对应真实三维空间中接收阵元所在位置的合成声压信号；多参数黏声波动方程的输入变量为原始声速图、原始声阻抗图和原始衰减品质因子图所组成的参数矩阵；
[0050] 步骤S130：利用反卷积目标函数和二范数目标函数求解所述观测声压信号和所述合成声压信号的差异；
[0051] 步骤S140：将所述反卷积目标函数和所述二范数目标函数通过伴随状态法计算得到关于所述参数矩阵的一阶梯度；基于所述一阶梯度和所述参数矩阵计算得到新的参数矩阵；将所述新的参数矩阵转换为新声速图、新声阻抗图和新衰减品质因子图；
[0052] 步骤S150：将新的参数矩阵输入到多参数黏声波动方程，结合反卷积目标函数、一阶梯度进行迭代运算，获取下一个新的参数矩阵；迭代运算停止条件为合成声压信号与观测声压信号达到一致。
[0053] 在一个实施例中，超声CT数据采集装备具备多层环形阵列，单次激发后多层阵列同时采集数据，获取具有纵向孔径的三维容积数据；
[0054] 采集数据装备为含有多个阵元的观测阵列换能器；多个阵元的每一个阵元能发射、接收超声信号；
[0055] 可选的，观测阵列换能器选用环形探头、对穿线阵探头、对穿面阵探头和对穿凹阵探头。
[0056] 在一个实施例中，多个阵元的每一个阵元所在位置的观测声压信号为dobs(xr，t；xs)，其中，xs为发射阵元的坐标位置，xr为接收阵元的坐标位置，t表示接收时间。
[0057] 在一个实施例中，多参数黏声波动方程，即为正演函数，其三维方程公式如下：
[0058]
[0059] 其中，描述声场的质子速度为vi，应力为σij，合成声压p＝∑kσkk/3；描述介质物理性质的参数密度为ρ，声速为Vp，声阻抗Ip＝ρVp描述，衰减系数由品质因子Q通过滞弹系数引入，该滞弹系数采用了L个在参考频率 下的记忆变量 如不考虑衰减，则最后两个方程及与 相关的系数都会被舍去，退化回传统的各向同性声波方程，其二维形式可通过舍去与z相关项退化形成，变量vi，σij，ρ，Vp， 均为定义在计算成像空间x和时间t上的图像；
[0060] 给定发射阵元(位置xs)的信息，输入声速图、声阻抗图、声衰减图的参数矩阵，可以通过数值仿真方法(如有限差分、有限元等方法)离散地求解在整个计算空间x和时间t上的声压信息p(x，t)，然后提取出接收阵元xr处合成声压信号dcal(xr，t；xs)＝R(x，xr)p(x，t；xs)，最大记录时间T应可记录相距最远的发射和接收阵元组合产生的穿过目标人体组织的透射波信号。
[0061] 在一个实施例中，反卷积目标函数，即为反演函数，使得FWI算法的结果更接近全局最小值，求解实测数据与仿真数据之间的时间域维纳滤波器，通过对该滤波器加权设计目标函数，可有效突出两种数据间基于波动方程的渡越时间差异，弱化正演振幅不准确对反演的影响，反卷积目标函数可表示为：
[0062]
[0063] 其中，A(τ)为权重函数，w表示维纳滤波器，
[0064]
[0065] 在反演时，权函数A(τ)的选择至关重要，权函数可以决定算法的鲁棒性，|τ|/τmax、2
τ和高斯函数均可作为选择。多尺度反演的特点可以通过设置权重函数种中的τ值上限来实现，因为该上限表征目标函数可考虑的观测数据和仿真数据的最大渡越时间差。该上限应结合分频多尺度反演选取，在相应频段数据主频倒数的3‑5倍之间，由于频率逐渐升高因此也实现了τ上限逐渐减小的效果，再利用传统二范数目标函数Cl2迭代提高分辨率，[0066]
[0067] 传统二范数目标函数在稳定性、准确性和高效性等方面还有待改进：在数据中缺少低频信息、初始模型离真实模型较远的情况下，优化算法难以有效收敛到全局最小值获取精准的超声CT图像重建。
[0068] 在一个实施例中，反卷积目标函数关于所述参数矩阵一阶梯度 可通过伴随状态法计算，
[0069]
[0070] 其中，伴随态波场算子μ为方程 的解。改变反卷积目标函数不影响正传波场w的计算。但是，伴随态波场μ是通过反传 获取，而该偏导数会随反卷积目标函数改变而变化。在二范数目标函数情况下，该偏导数为：
[0071]
[0072] 在反卷积条件下，该偏导数为：
[0073]
[0074] 将 通过波动方程反向传播即可计算伴随态波场算子μ。
[0075] 在一个实施例中，迭代可视为求解求优化问题的过程，在第0次迭代输入初始的参数矩阵m0，即基于先验知识的对声参数模型的估计；
[0076] 可选的，在检测乳腺等软组织的情况下，可以假设起始参数矩阵m0等效于均匀的水介质。
[0077] 在一个实施例中，迭代更新一般采用梯度下降法，沿着下降方向d迭代更新所述参数矩阵的值，
[0078] mk+1＝mk+αkdk，
[0079] 其中，αk一般利用线搜索方法获取，而dk在不同优化策略(如：最速下降法、共轭梯度法、L‑BFGS法、牛顿法)中有不同的计算方法，最速下降法一般利用目标函数关于模型的梯度的负方向作为下降方向，即 共轭梯度、L‑BFGS等方法需要利用前置迭代中的梯度信息改善当前梯度获取下降方向，牛顿法需估计海瑟矩阵，即目标函数关于梯度的二阶导数，来约束当前梯度获取下降方向。
[0080] 在一个实施例中，迭代运算停止条件可以选用反卷积目标函数值下降的阈值或者迭代的次数。
[0081] 在一个实施例中，求解多参数黏声波动方程、反卷积目标函数的代码可利用CPU进行验证，然后移植到GPU上提高效率以应对时间敏感的超声CT医疗成像临床诊断任务需求。
[0082] 图2中，以单个环形阵列利用全矩阵采集探测乳腺仿体为例，图3(a)为本实施例中人体组织仿体俯视图，在假定水浸观测的条件下，具体说明本申请的实施方式。选用观测系统为单个环形阵列，通过水耦合的方式探测前述乳腺仿体声学参数，仿体含有1.5‑2.0mm微型缺陷，系统发射脉冲中心频率此处为1.0MHz，具体流程如下：
[0083] 步骤1：使用环形阵列包围仿体后，利用全矩阵的采集方式获取Ns×Ne×Nt大小的射频观测数据dobs，进行预处理。其中，数据采集装备为含Ne个阵元的阵列换能器，Ns为发射次数，Nt为时间采样点个数，且Ns≤Ne。
[0084] 步骤2：选定基于声速密度参数的二维声波方程，采用二阶时间差分、四阶空间交错网格差分的偏微分方程解法，结合卷积完美匹配层消除计算边界回声，通过库朗数确定时间和空间网格尺寸，并重采样观测数据以保证相同的时间采样率。然后利用水的介质属0 3
性对模型进行初始化，初始模型中声速 为均匀的1490m/s，密度ρ为1000g/cm ，即最后基于此初始模型估算激发阵元处的子波时间函数s(t)，然后在激发阵
元处注入子波函数生成正传波场w，此后可在接收位置处提取仿真数据dcal。此处，可记录可压缩后的正传波场或其边界，边界可用来反向重构正传波场。
[0085] 步骤3：选定目标函数，计算目标函数关于仿真数据的偏导数 此偏导数也称为残差。将此残差通过反向传播波动方程的方法计算伴随波场μ，且在对应的时间步上，将反传波场μ与记录或重构的正传波场 做零延迟的互相关，并与系数矩阵作用。反传完成后即可获得目标函数关于模型参数的梯度
[0086] 步骤4：利用步骤3获取的梯度通过优化算法(如：最速下降法、共轭梯度法、LBFGS、牛顿法等)计算下降方向，并与预条件算子矩阵作用，再利用线搜索方法获取更新步长。此k+1 k k时，在k次迭代时，可利用处理过的下降方向更新现有参数模型，即m ＝m +αkd 。之后回到步骤3迭代循环，直到达到预定的终止条件。
[0087] 步骤5：前述步骤2‑4为内循环，实际操作中，可通过对观测数据滤波形成不同中心频率的观测数据。然后，从低频开始逐渐升高频率，在每个频段可以进行一次步骤2‑4的处理，下一个频段采用上一个频段的输出模型作为初始模型。不同频段也可采用不同目标函数。
[0088] 步骤6：达到中止条件，停止迭代。
[0089] 步骤7：输出多参数图像，图3(b)为本施例中基于FWI的超声CT声阻抗图像，图3(c)为本施例中基于FWI的超声CT声速图像。
[0090] 本发明提出了一种基于波动方程的超声CT成像，具体来说，即是基于FWI的超声CT多参数图像重建系统和方法，该方法通过利用多参数黏声波动方程的数值解，可以充分利用全部波形的信息(反射波、透射波、绕射波、折射波等)重建多种声参数(声速、声阻抗、衰减系数等)的高分辨率图像。多参数波动方程的数值解法可在CPU上做开发，部署在GPU上用于实际的探测任务。构建创新的多尺度反卷积目标函数并提出相应目标函数在实际重建中的具体使用策略。通过着重比较仿真数据和观测数据的动力学特征的思路提高了模型重建的鲁棒性。
[0091] 需要说明的是，在此提供的方法不与任何特定计算机、虚拟装置或者其它设备固有相关。各种通用装置也可以与基于在此的示教一起使用。根据上面的描述，构造这类装置所要求的结构是显而易见的。此外，本发明也不针对任何特定的编程语言。应当明白，可以利用各种编程语言实现在此描述的本发明内容，并且上面对特定语言、系统功能模块的调用所做的描述仅仅是为了披露发明的最佳实施方式。
[0092] 在此处所提供的说明书中，说明了大量的具体细节。然而，能够理解，本发明的实施例可以在没有这些具体细节的情况下完成实现。在一些示例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。
[0093] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若对本发明的这些修改和变型属于本发明权利要去及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。
