具体实施方式

以下描述用于揭露本发明以使本领域技术人员能够实现本发明。以下描述中的优选实施例只作为举例，本领域技术人员可以想到其他显而易见的变型。

逆时偏移成像是利用反演得到的速度模型，对采集数据做逆时偏移成像，取得最终影像。

逆时偏移(reverse time migration)最早于1978年由法国人hemon提出，由Baysal、Whitemore、McMechan等人应用于纵波资料的叠后偏移，获得好效果。

逆时偏移成像的方法是通过双程波波动方程在时间域上对人工给予的震源子波正向传播和接收到的地震资料进行反向传播，结合成像条件实现偏移(Claerbout，1971)。由此可见，一方面由于逆时偏移需同时计算震源正向传播和地震资料的反向传播过程中各个时刻的波场，因其各自时间的延拓方向不同，而致使实际应用中必须存储其中一个方向的传播过程，这就是何以需要甚大额外存储空间的原因所在。如今的地震资料采集已经发展到了三维甚至四维，其数据量更是十分庞大，逆时偏移在推向实际应用时尚需面对此类业内皆知的困难。另一方面，相较于单程波方程的波场延拓而言，逆时偏移运用的是双程波波动方程进行波场延拓，因其避免了上下行波的分离处理，因而成为最准确的成像算法，且不受倾角的限制，并能实现回转波和多次波成像。正基于此，地球物理学家从来没有放弃过对逆时偏移的追求。研究的焦点也集中于提高逆时偏移的计算效率和减少偏移中的存储I/O问题。成像条件的制定是地震偏移成像算法的关键之一，它直接影响成像的效果和计算成本。成像条件由Claerbout提出，他同时提出两种成像条件：相关成像条件和褶积成像条件，又有学者称这两种成像条件分别为乘法成像条件和除法成像条件。

波动方程的高阶有限差分近似：

三维声波方程的表达如下：

其中，t代表时间向量坐标，x,y,z代表三维空间向量坐标，P表示位移函数，V表示速度函数。分别代表利用时间二阶中心有限差分近似利用空间高阶中心有限差分近似可得：

其中，

再利用N阶中心差分近似得：

其中：

分别用Δx,Δy,Δz表示差分网格的间距、f_max表示子波频率的最大值、V_min表示速度模型的最小值，那么差分格式的频散条件可表示为：

其中，n表示一个波长用n个样点来采样，h表示最大网格间距。

那么对于该算法而言，其稳定性条件可以表示为：

褶积成像条件(除法成像条件)的表达式可以表示为：

其中，x,y,z代表三维空间向量的坐标，U(x,z,t)为上行波波场，D(x,z,t)为下行波波场,dt为延拓时间步长，I(x,z)表示成像结果，tmax表示最大旅行时。其物理含义为反射波能量除以入射波能量可以求出反射系数，但是这种成像条件的实际应用中是有稳定性问题存在的(有可能除零)。

所以，通常采用的是相关成像条件(乘法成像条件)，其表达式可以表示为：

双程波逆时偏移中以外推观测波场P_R(x,z,t)取代上行波，以震源外推波场P_S(x,z,t)取代下行波，因此公式(6)式变为：

在公式(7)式中的被积函数P_R(x,z,t)P_S(x,z,t)表示t时刻对整个波场做一次成像运算，积分说明像空间I(x,z)中的像是各时间步所成的像的叠加。所以，互相关成像条件充分利用了成像信息，在增强成像信号的同时也有效压制了成像噪音。逆时偏移互相关成像条件的成像步骤是：

首先要完成一次波动方程正演计算，并保存每一时间步的波场信息。接着，逆时外推记录波场，存每一时间步的波场信息。最后，分别读取保存的同一时刻的震源波场和记录波场作成像运算。即P_R(x,z,t)P_S(x,z,t)，再累加入成像空间。

逆时偏移成像的条件需要使用在同一时刻的震源波场(经过激发和正演模拟波的传播)和记录波场(逆时反传播回地下的波场)，因为一个是正传波场，另一个是反传波场，如果想同时得到相同时刻的两个波场，则必须存储其中一个波场的传播过程，即每一时刻的波场分布，需要消耗甚大的存储资源，这个要求在实际操作中是难以满足的。目前存在一些解决策略，主要包括：

1)策略：最容易想到的方法就是震源和记录同时从T_n(表示n*dt时刻)延拓到T₀，这样就不需要额外的存储空间，并且可以达到最小计算量O(2N)的量级，然而T_n时刻的震源波场不是已知数据，所以不能直接实现；

2)策略：存储T₀到T_n的波场，然后反传记录波场，当传播的某个时刻时读取那个时刻的激发点波场。此方法也可以达到O(2N)计算量的要求，但是需要巨大的磁盘空间来存储波场，即使磁盘空间可以满足需求但是因此产生的I/O时间消耗十分庞大；

3)策略：只存储少数几个时刻的震源波场，当记录波场反向延拓的过程中，利用这些存储的震源波场插值来近似当前时刻的波场然后成像，但是这种处理方式对波场的插值首先是不准确的，并且计算量也随之增加；

4)策略：首先把震源波场正向延拓到T_n和T_n-1两个时刻，把这两个波场作为初始条件，与记录波场同时反向延拓，并且随延拓过程中应用成像条件，这样就不需要额外的存储空间。这种策略只可以应用在密度为常数的介质中，并且激发点边界条件是Dirichlet边界条件。但是这种波场的延拓假设是在无限大空间内传播的，而通常的计算方法是加入了人工边界的，包括自由边界条件、吸收边界条件、指数衰减边界条件和完美匹配层边界条件等。其中除自由边界条件外，其他三种边界条件都会破坏波场的完整性，使得波场的延拓不可逆。为此只能在适当的时刻在适当的空间位置上插入波场能量以补充波场的完整性，这样就需要存储每个时刻的位于人工边界处的波场值，当计算规模增大后，存储量也会快速上升，尽管小于2)策略中需要的存储量，但是也是通常不能满足实际的要求。

5)策略：首先正向延拓波场，并且记录多组检查点(P_n，P_n-1)作为中间时刻的波场初始条件，然后反向延拓检波点波场到某时刻，用最近的检查点来正向传播得到该时刻的震源波场，然后应用成像条件，当应用PML(完全匹配，Perfectly Matched Layer)边界条件的时候波场还需要额外的存储空间，并且需要大量的重复计算过程，最大的重复计算可能会达到O(mN)的量级，其中m为检查点的个数。Griewank利用这个思路提出优化的检查点的选取方法，可以使得重复计算的次数减少到O(Nlog^m)的量级尽管如此，逆时偏移的消耗时间和存储的矛盾也并未得到完全解决，但是该方法是目前最经济的计算策略。

在上述5种逆时偏移的解决策略中可以看出计算量级与存储规模是一对矛盾，要是想得到最小的计算量就会达到最大的存储量，反之亦然。然而在策略四中不需要额外的存储空间，而其他策略都要涉及到存储规模与计算量级的平衡问题。在策略四中，人工边界是主要的影响因素，然而自由边界条件又会引入很多的成像假象，其他的边界条件又不复合该策略的要求只能舍弃。

在本发明中，主要针对上述的4)策略。2009年Robert提出随机边界模型，其思想是消除人工边界自由边界条件反射波的相干性，使边界反射不能成像。具体实现过程为：将有限空间外扩一定的距离，然后在外扩的空间内填充随机速度，从而形成随机边界速度模型，当波传到随机速度区域时波前面将被随机化，使波场变成随机噪音反传回真实速度区域，破坏了边界反射的相干性使得边界反射不能成像。

构造随机边界函数如下：

其中，x,y,z代表空间向量坐标，为边界点的随机速度函数，V(x,y)为边界点的原始速度函数，r是随机数，d为速度点与内层边界的空间距离。

可以运用处理单元，例如CPU/GPU作为偏移计算核心，采用随机速度边界模型，将激发点波场传播至T_n时刻，然后利用T_n时刻的波场作为初始条件，同时反传激发点波场和检波点波场，并且同时利用成像条件进行成像。进而避免了额外的存储空间，虽然激发点波场的重复计算一次，但利用GPU作为波场延拓的数值计算核心其耗时与大规模的磁盘I/O相比还是非常经济的。

根据上述内容，请参阅图1，本发明的医学影像逆时偏移成像方法包括以下步骤：

S100)获取模型参数，根据模型参数，确定时间延拓步长以及偏移孔径范围，在本实施例中，模型参数包括偏移参数和单炮数据，且可以是先根据偏移参数确定延拓步长，再根据单炮数据确定偏移孔径范围，另外，单炮数据引用了地震勘探中shot的概念，在这里可以对应一次波场发生装置例如超声波换能器发射后，所有接受波场发生装置记录的数据为一个shot，也就是一个单炮数据。

S200)在孔径范围内添加随机速度边界，生成孔径范围内的随机速度边界模型。

S300)确定最大接收时间，根据最大接收时间，获取反向传播最后两个时刻的正演震源波场，即T_n和T_n-1两个时刻的正演震源波场，最大接收时间是通过时间延拓步长确定，具体的，请参阅图2，包括以下步骤：

S301)生成震源子波。

S302)根据震源子波和随机速度边界模型，利用有限差分的方法正演震源波场到最大接收时间(GPU加速)。

S303)获取反向传播最后两个时刻的正演震源波场(GPU加速)。

S400)自最大接收时间，将单炮数据反向传播回地下，具体的，请参阅图3，包括以下步骤：

S401)对单炮数据规则化，即得到规则化的单炮数据。

S402)根据孔径范围，生成孔径范围内的真实速度模型。

S403)根据规则化的单炮数据和真实速度模型，自最大接收时间，将自地表反向传播接收的单炮数据再反向传播回地下(GPU加速)。

S500)根据逆时偏移成像的条件(GPU加速)，进行初步的医学成像。

S600)存储初步的医学成像结果。

S700)去除初步的医学成像结果的成像噪音。

S800)恢复成像结果的空间坐标位置信息，得到最终的医学成像结果。

进一步的，请参阅图4，本发明提供的医学影像逆时偏移成像装置包括：

参数获取模块，用于获取模型参数，根据模型参数，确定时间延拓步长以及偏移孔径范围，在本实施例中，模型参数包括偏移参数和单炮数据，且可以是先根据偏移参数确定延拓步长，再根据单炮数据确定偏移孔径范围，另外，单炮数据引用了地震勘探中shot的概念，在这里可以对应一次波场发生装置例如超声波换能器发射后，所有接受波场发生装置记录的数据为一个shot，也就是一个单炮数据。

模型生成模块，用于在孔径范围内添加随机速度边界，生成孔径范围内的随机速度边界模型。

波场获取模块，用于确定最大接收时间，根据最大接收时间，获取反向传播最后两个时刻的正演震源波场，即T_n和T_n-1两个时刻的正演震源波场，最大接收时间是通过时间延拓步长确定，具体的，请参阅图5，包括以下步骤：

子波生成单元，用于生成震源子波。

正演单元，用于根据震源子波和随机速度边界模型，利用有限差分的方法正演震源波场到最大接收时间(GPU加速)。

获取单元，用于获取反向传播最后两个时刻的正演震源波场(GPU加速)。

反向传播模块，用于自最大接收时间，将单炮数据反向传播回地下，具体的，请参阅图6，包括以下步骤：

规则化单元，用于对单炮数据规则化，即得到规则化的单炮数据。

模型生成单元，用于根据孔径范围，生成孔径范围内的真实速度模型。

反向传播单元，用于根据规则化的单炮数据和真实速度模型，自最大接收时间，将自地表反向传播接收的单炮数据再反向传播回地下(GPU加速)。

初步成像模块，用于根据逆时偏移成像的条件(GPU加速)，进行初步的医学成像。

存储模块，用于存储初步的医学成像结果。

去噪模块，用于去除初步的医学成像结果的成像噪音。

成像确定模块，用于恢复成像结果的空间坐标位置信息，得到最终的医学成像结果。

本发明的有益效果为：

本发明可以实现低成本高分辨率的医学成像，可以以超声成像的硬件成本，实现核磁共振的成像精度(精度可达0.5毫米)，甚至更高精度。相比CT，没有辐射。比起传统超声成像，可以利用所有反射、折射、绕射、多次等复杂波场，成像精度更高，对于特殊部位(如头颅)传统超声波无法精确成像的情况下，依然可以提供高精度准成像结果。使用到本发明后，采集效率高，一次采集一般需要2-5秒，对人体运动不是很敏感。

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明的范围内。