CN103698810A - 混合网最小走时射线追踪层析成像方法 - Google Patents

Abstract

本发明公开了一种混合网最小走时射线追踪层析成像方法，它是将二维复杂外部几何边界及内部复杂空间结构的成像区域，根据其几何结构分解成规则矩形子域和不规则子域，对规则矩形子域采用矩形网格剖分，对不规则区域则采用三角网格剖分；采用混合网波行面扩展法计算声波初至走时，追踪出任意接收点至激发点处的射线路径；将射线路径信息形成矩阵方程；再采用外推序列加速代数重构技术求解矩阵方程得到模型的速度分布，最终形成规则的网格成像数据。本发明有效解决单种网格剖分存在的问题以及正演模拟时存储量大、计算时间长、正演模拟难度增大、效率低、层析反演时方程性态差、求解困难问题。

Description

混合网最小走时射线追踪层析成像方法

技术领域

本发明涉及地学层析成像技术，主要应用于岩体工程地质声波探测、孔间地质构造探测及各种人工结构体缺陷检测中的计算机数据图像处理。 

背景技术

目前，地学层析成像技术以射线代数重建法为主，称为射线层析成像方法。 

具体方法是，在地层中打两孔，一孔设有信号激发装置，另一孔设有信号接收装置（单道或多道），形成扇形观测系统，测量探测区域介质速度。通过改变激发点和接收排列的位置，组成密集交叉的射线网络，然后根据射线的疏密程度及成像精度划分规则的成像单元，运用射线追踪理论，采用反演计算方法形成被测区域的声波波速图像。根据射线理论，初至旅行时t与探测区域介质速度v(x,z)有下列关系 

$t=\underset{R\left(v\right)}{\∫}\frac{1}{v(x,z)}\mathrm{ds}---\left(1\right)$

式(1)中，R(v)为发射点到接收点间的路径，为速度v(x,z)的函数；ds为沿射线路径R(v)的距离增量。将射线穿过的区域离散成网格模型，则可建立如下反演控制方程： 

AX=B   (2) 

式(2)即为层析成像成像方程，其中，A是M×N阶矩阵，M为观测射线条数，N是单元个数，A的元素a_ij是第i次观测中传播路径与第j个网格有关的微分系数，i=1，2，…，M，j=1，2，…，N；X是N维列向量，其元素x_j是第j个网格中的慢度（速度v_j的倒数）；B是M维列向量，其元素t_i是第i次观测得到的初至声波走时。 

式(1)中的初至旅行时与速度分布的关系是非线性的，层析反演需要对(1)式进行线性化处理。迭代一般从某一初始速度模型v₀开始，建立观测旅行时与理论旅行时之差δt满足的线性方程，解出速度的扰动量δv，修正速度模型v＝v₀+δv，得到新的速度模型。矩阵A是通过射线追踪算法追踪激发点与接收点之间的射线路径而形成的，成像效果的优劣取决于射线追踪方法。 

由于工程地质探测现场条件的复杂性，所探测研究对象的外部几何形状大部分情况下是不规则的形状；内部地质异常体在空间结构上更加复杂，有的呈不规则柱状体（如溶洞、陷落柱等），有的呈条带状（如断层、裂缝等）、有的呈层状（如软弱夹层等）。这就对模型参数化提出更高的要求，以满足复杂结构体的探测。 

目前模型参数化一般都采用单一的矩形网或三角网，称为单一网格射线追踪层析成像法。 

在该方法中，矩形网格参数化时，网格大小和形状通常是相同的，这样可使得网格间检索和数据点定位简单；但这样的处理方式也有很多缺点： 

（1）模型剖分的灵活性差，矩形网格大小和形状通常是规则的，这样就不能根据地质构造的复杂程度进行区别对待，为迁就对复杂构造区域的采样精度，矩形网格的尺度需要很小，而网格剖分的均一性要求在构造简单区域必须以同样的尺度划分，增大了剖分网格的数目。 

（2）对速度界面的描述精度差，速度模型和界面模型不一致。界面与网格边界不重合，在用速度值的相对大小表征界面的位置时只能用锯齿状的分布来逼近平滑变化的界面。为了保证逼近的精度，则网格的尺寸必须很小，因而导致网格的数目很多。 

在该方法中，相对于矩形网格参数化，三角网格参数化具有如下优点： 

（1）模型剖分的灵活性强，剖分方式非常灵活，网格大小和形状可根据探测 区域外部几何形状及地下不同区域构造的复杂程度而灵活设置，在构造简单区域采用大网格剖分，在构造复杂区域采用小网格剖分，总体上网格数目少； 

（2）对速度界面的描述精度高，速度模型和界面模型具有一致性。三角网格的弯曲边界和界面的延伸方向完全一致，界面位于三角网格间的分界面上，对速度界面的描述精度高。 

但在该方法中，三角网格参数化也存在如下缺点： 

（1）正演模拟时存储量大、计算时间长。模型剖分的灵活性差导致剖分网格的数目很多，与正演模拟结合时，由于在Frechet微商计算中要计算射线在每个网格内的出射点和射线长度，由于网格数目众多，因此在总体上存储量大、计算耗时长。 

（2）用于层析反演时存储量大、计算时间长、方程性态差、求解困难；网格参数化与层析反演相结合时，更多的网格意味着更多的未知参数，因此需要更多的射线，则待求解方程组的维数大大增加。此外更多的未知参数意味着参数矢量中欠定的和位于零空间的分量数目更多，性态大大变坏，层析反演的难度增加，通常需要加入正则化约束改变方程组的性态。而正则化约束的加入又增加了存储量和计算量。 

发明内容

为了克服地学层析成像技术在模型参数化时的单一网格射线追踪层析成像法存在的不足，本发明提供一种混合网射线追踪层析成像法。 

为达到上述目的，本发明采取的技术方案是： 

一种混合网射线追踪层析成像法，是采用计算机程序实现的，其特征在于，它是将二维复杂外部几何边界及内部复杂空间结构的成像区域，根据其几何结构分解成规则矩形子域和不规则子域，对规则矩形子域采用矩形网格剖分，而对不规则区域则采用三角网格剖分，从而将整个区域网格剖分成包括矩形网和三角网的混合网；采用混合网波行面扩展法计算声波初至走时，追踪出任意接收点至激发点处的射线路径，射线路径包含了从接收点至源点的坐标、所在网格单元等信 息；将射线路径信息形成矩阵方程，则区域层析成像则归结为求解大型病态线性方程组；再采用外推序列加速代数重构技术求解矩阵方程，得到模型的速度分布，即最终形成规则的网格成像数据。 

下面对上述混合网射线追踪层析成像法进一步描述，详细步骤如下： 

第1步、数据采集，即拾取各震源至接收换能器的初至走时t_m； 

第2步：建立成像区域几何结构； 

根据边界几何特征、先验信息等建立成像区域的背景网，包括初始网格密度的设置等，为混合网剖分打下基础，同时建立成像区域的相三维对坐标系统； 

第3步：建立成像区域初至走时信息； 

初至走时信息包括：每一震源点坐标、对应每一震源点接收排列的三维坐标和走时t_m； 

第4步：将第2步成像区域的相对坐标系统进行旋转变换和切面投影，以有利于二维成像操作； 

第5步：建立混合型区域网格； 

首先对二维复杂外部几何边界的地学层析成像区域、以及内部复杂空间结构的地学层析成像区域的几何结构进行分解，分解成规则矩形子域和不规则子域，对规则矩形子域采用矩形网格剖分，对不规则区域则采用三角网格剖分，从而使整个成像区域网格成为矩形网和三角网混合组成的区域网格； 

按背景网及当前速度模型设置的网格密度函数值进行剖分，形成网格单元节点坐标数据文件； 

第6步：建立网格各节点的拓扑关系，为射线追踪打下基础； 

上述的网格各节点包括：源节点、接收节点、网格单元顶点节点、网格边插入节点； 

第7步：初始速度模型输入； 

当已知先验速度，例如已知声速测井资料等信息时，应首先输入先验速度信 息建立初始速度模型；否则，按均匀速度信息建立初始速度模型； 

第8步：对每个源点进行射线追踪，得到相应源点下各节点的次级源及次级源所在的网格单元；射线追踪步骤如下： 

第8.1步：波行面扩展 

从震源或一个相对最小走时波前点出发，确定这一节点的当前最小走时单元，如三角网子域时三角形单元为当前最小走时单元，如四边形子域时四边形所在单元为当前最小走时单元； 

先计算当前最小走时单元节点的最小旅行走时，再计算该单元相邻单元节点旅行走时；将当前单元和相邻单元节点的旅行走时相比较，如果旅行走时无变化时，对当前这一最小走时单元作出“休眠”标志，再寻找下一个最小走时单元作为当前单元，以此类推，遍历所有单元； 

在上述遍历所有单元过程中，把遍历单元作为一个波行平面，当该波行平面遇回转波时，先前已经“休眠”的单元将通过单元之间的传递而“激活”，再次成为当前最小走时单元；当所有波行面的节点走时都不再变化、且波行平面中所有的单元均处于“休眠”状态时，波行面扩展完毕； 

第8.2步：利用计算出的各节点的旅行时和方向信息确定射线路径，从接收点开始，根据各节点走时及次级源信息，进行向源检索，拾取从接收点到源点的射线路径； 

第9步：设置成像反演迭代的终止条件，终止条件可以是迭代次数，也可以是预先设置走时残差； 

第10步：接收点的向源检索； 

对每个源点对应的接收点进行向源检索，得到射线路径，射线路径包含了从接收点至源点的坐标和所在网格单元信息，所在网格单元信息包括网格单元的编号、次级源节点坐标和射线长度，据此形成Jacobi矩阵A； 

由于Jacobi矩阵A为一稀疏矩阵，采用压缩存储，只需记录射线所经历的网 格单元及在此网格单元的长度即可； 

判断成像反演迭代终止条件，若满足终止条件，则转下一步；否则转第8步； 

若想根据射线密度、速度分布信息重新对成像区域离散化时，则转第5步；如此反复进行，直至满足设置的终止条件； 

第11步：反演求解 

根据第10步形成的Jacobi矩阵A，计算理论走时t_c和走时残差，走时残差δt＝t_m-t_c；建立反演方程Aδv＝δt；采用外推加速代数重建反演算法解方程求取δv，修改模型v＝v₀+δv；终止条件可以是迭代次数，也可以是预先设置走时残差； 

第12步：成像结果输出 

对于规则网格成像结果输出时，只需输出网格各单元中心坐标及单元的速度值即可； 

对于不规则网格成像结果输出时，采用“采样”输出的办法；即在不规则网格单元的横向和纵向上等间距的设置多个采样点，输出采样点的坐标及单元的速度值，最终形成规则网格成像数据。 

本发明的积极效果是，采用混合网射线层析成像方法可有效解决单种网格剖分存在的问题：对不规则区域来说，有效地解决了矩形网剖分对复杂区域网格参数化灵活性差，速度界面描述精度低等问题；对规则区域来说，采用矩形网剖分可大大减少网格数目，从而可有效地解决正演模拟时存储量大、计算时间长、正演模拟难度增大、效率低、层析反演时方程性态差、求解困难等问题。 

附图说明

图1是以三角网射线追踪为例的波行平面从(a)到(h)的扩展过程图； 

图2是实施例的三角网节点次级源检索图。 

具体实施方式

下面参照附图，对本发明的混合网射线追踪层析成像法进一步描述。 

第1步、数据采集，即拾取各震源至接收换能器的初至走时t_m； 

第2步：建立成像区域几何结构； 

根据边界几何特征、先验信息等建立成像区域的背景网，包括初始网格密度的设置等，为混合网剖分打下基础，同时建立成像区域的相对三维坐标系统； 

第3步：建立成像区域初至走时信息； 

初至走时信息包括：每一震源点坐标、对应每一震源点接收排列的坐标、和走时t_m； 

第4步：将第2步成像区域的相对三维坐标系统进行旋转变换和切面投影，成为二维坐标系统； 

第5步：建立混合型区域网格； 

首先对二维复杂外部几何边界的地学层析成像区域、以及内部复杂空间结构的地学层析成像区域的几何结构进行分解，分解成规则矩形子域和不规则子域，对规则矩形子域采用矩形网格剖分，对不规则区域则采用三角网格剖分，从而使整个成像区域网格成为矩形网和三角网组成的区域； 

按背景网及当前速度模型设置的网格密度函数值进行剖分，形成网格单元节点坐标数据文件； 

第6步：建立混合网格各节点的拓扑关系，为射线追踪打下基础； 

上述的混合网格各节点包括：源节点、接收节点、网格单元顶点节点、网格边插入节点； 

第7步：初始速度模型输入； 

当已知先验速度，例如已知测井资料获得的声速信息时，应首先输入先验速度信息，建立初始速度模型；否则，按均匀速度信息建立初始速度模型； 

第8步：对混合网每个源点进行射线追踪，得到相应源点下各节点的次级源及次级源所在的网格单元，混合网射线追踪的详细步骤将根据图1图2在以后描述； 

第9步：设置成像反演迭代的终止条件，终止条件可以是迭代次数，也可以是预先设置走时残差； 

第10步：接收点的向源检索； 

对混合网每个源点对应的接收点进行向源检索，得到射线路径；射线路径包括从接收点至源点的坐标、该源点所在网格单元的信息（网格单元的编号、次级源节点坐标和射线长度），据此形成Jacobi矩阵A； 

由于Jacobi矩阵A为一稀疏矩阵，采用压缩存储方式，只需记录射线所经历的网格单元编号及射线在此网格单元中的长度即可； 

判断成像反演迭代终止条件，若满足终止条件，则转下一步；否则转第8步； 

若想根据射线密度、速度分布信息重新对成像区域离散化时，则转第5步；如此反复进行，直至满足设置的终止条件； 

第11步：反演求解： 

根据第10步形成的Jacobi矩阵A，计算理论走时t_c和走时残差，走时残差δt＝t_m-t_c，建立反演方程Aδv＝δt；采用外推加速代数重建反演算法解方程求取δv，修改模型v＝v₀+δv；终止条件可以是迭代次数，也可以是预先设置走时残差； 

上述外推加速代数重建反演算法解方程求取δv的详细步骤如下： 

1、建立线性方程组 

Ax＝b   (3) 

其中，A为已知方阵A＝(a_ij)_n×n，b是一组给定的n维向量b＝(b₁,b₂,…,b_n)^T，未知量x＝(x₁,x₂,…,x_n)^T； 

方程(3)中的第i个方程可以写成： 

(a_i,x)＝b_i(i＝1,2,…,n)   (4) 

其中，a_i是A的第i行向量的转置，即a_i＝(a_i1,a_i2,…,a_in)^T；b_i是向量b的第i个 分量，(·，·)表示向量内积。(4)式的每一个方程都代表n维空间中的一个超平面，而a_i是该超平面的法线向量； 

假定已知x^(k)，则对一个固定的i，有式 

x^(k+1)＝x^(k)+βa_i   (5) 

(5)式表明x^(k+1)是从x^(k)出发在第i个超平面的法线方向上选取的，令x^(k+1)应该满足： 

(a_i,x^(k+1))＝b_i   (6) 

将(5)式中的x^(k+1)代入(6)式得： 

$\mathrm{\β}=\frac{b_i-(a_i,x^{\left(k\right)})}{(a_i,a_i)}=\frac{b_i-(a_i,x^{\left(k\right)})}{{\left|\right|a_i\left|\right|}^2}----\left(7\right)$

从而有： 

$x^{(k+1)}=x^{\left(k\right)}+\frac{b_i(a_i,x^{\left(k\right)})}{{\left|\right|a_i\left|\right|}^2}{a}_i---\left(8\right)$

这时，可以把ART迭代方法写成如下形式： 

其中，λ_k叫做松驰因子（0＜λ_k＜2）； 

2、解方程组： 

（1）由迭代初值x⁽⁰⁾开始，按格式（10）迭代l步，得到x⁽⁰⁾，x⁽¹⁾，…，x^(l)，l（l≧2）为选定的某一正整数； 

（2）利用Aiten外推公式得到加速序列： 

${\stackrel{\‾}{x}}_j^{\left(l\right)}=x_j^{\left(l\right)}-\frac{{(x_j^{\left(l\right)}-x_j^{(l-1)})}^2}{x_j^{\left(l\right)}-{2x}_j^{(l-1)}+x_j^{(j-2)}}(j=\mathrm{1,2},\·\·\·,n)----(11)$

其中 $x^{\left(l\right)}={(x_1^{\left(l\right)},x_2^{\left(l\right)},...,x_n^{\left(l\right)})}^T,$ ， ${\stackrel{\‾}{x}}^{\left(l\right)}={({\stackrel{\‾}{x}}_1^{\left(l\right)},{\stackrel{\‾}{x}}_2^{\left(l\right)},...,{\stackrel{\‾}{x}}_n^{\left(l\right)})}^T$

（3）以

为迭代初值，反复上述（1）、（2）两步，直到满足收敛条件为止； 

第11步：成像结果输出； 

对于规则网成像结果的输出，直接输出各单元中心坐标及单元的速度值即可； 

对于不规则网成像结果的输出，采用“采样”输出的办法；“采样”输出时，在不规则网格单元的横向和纵向上等间距的设置采样点，这样一个网格单元可有多个采样点，最终形成规则的网格成像数据。 

下面根据附图1、图2，以三角网射线追踪为例，对第8步中所述的射线路径追踪计算方法详细描述如下： 

第一步：波行面扩展： 

利用波行平面，将波行平面的扩展视为一个递归过程，图1展示了波行面从(a)到(h)的扩展过程，源点在左上角，粗线三角单元表示当前三角单元； 

（1）设源点的走时为零，除源点外的其他节点的走时为无穷大。左上角为1、2号三角形的公共节点，只需将1号三角形纳入波行面三角单元存储容器，2号三角形单元会自行遍历，并使1号三角单元处于“运算”状态，作出存在于波行面的标志。 

（2）如图1(a)，在波行面三角单元存储容器中查找处于“运算”状态的三角单元中的最小走时节点，此时波行面三角单元存储容器只有1号三角单元，置1号三角单元作为当前三角单元，并将1号三角单元及其相邻三角单元2～6纳入波行面三角单元存储容器，标记为“运算”状态，同时做出存在于波行面的标志。 

（3）利用“波前邻域点最小走时及次级源位置计算”方法，计算、比较当前的1号最小走时三角单元各节点旅行时(从1号三角单元最小走时节点开始)；再计算、比较最小走时三角单元的2～6相邻三角单元各节点的旅行时，当相邻的2～6三角单元中的任一节点走时有变化时，将这一相邻三角单元作出“运算”标志；再一次计算、比较当前的1号最小走时三角单元节点旅行时，当1号最小走时三角单元各节点走时无变化时，对1号三角单元作出“运算暂停”标志。 

（4）如图1(b)，与第（2）步类似，在波行面三角单元存储容器中查找处于“运算”状态的三角单元中的最小走时节点，这时波行面三角单元存储容器中有1、2、3、4、5、6号三角单元，1号为“运算暂停”状态，2～6三角单位处于“运算”状态，则必找到左上角节点，这个节点位于2号三角单元。置2号三角单元作为当前最小走时三角单元，并将2号三角单元及其相邻三角单元1、3～7纳入波行面三角单元存储容器，并做出存在于波行面的标志，实际上1～6已存在于波行面三角单元存储容器，2号三角形只扩展7号三角形单元。 

（5）与第（3）步类似，计算、比较2号三角单元节点旅行时(从2号三角单元最小走时节点开始)；再计算、比较最小走时三角单元的1、3～7相邻三角单元节点的旅行时，当1、3～7三角单元中的任一节点走时有变化时，将这一三角单元作出“运算”标志；再一次计算、比较当前的2号最小走时三角单元节点旅行时，当2号三角单元各节点走时无变化时，对2号三角单元作出“运算暂停”标志。 

图1(c)～图1(h)的3～8号三角单元的扩展与1、2号三角单元扩展类似，即：3号三角单元的相邻三角单元为1、2、4、5、6、8、9、12、13，实际扩展了8、9、12、13单元；4号三角单元的相邻三角单元为1、2、3、5、7、8、10、11，实际扩展了10、11单元；5号三角单元的相邻三角单元为1、2、3、4、6、7、8、9、10、11、12、13，没有扩展；6号三角单元的相邻三角单元为1、3、5、8、9、12，13、14，实际扩展了14单元；7号三角单元的相邻三角单元为2、4、5、8、10、11，15，实际扩展了15单元；8号三角单元的相邻三角单元为3、4、5、6、7、9，10，11、12、13、16、17、18，实际扩展了16、17、18单元。 

以此方式遍历所有三角单元，当所有波行面的节点走时都不再变化，且所有波行面中的三角单元均处于“运算暂停”状态，波行面扩展完毕。 

与一些将波前刻画成一条曲线不同，本算法在于：将所遍历的三角单元作为一个波行平面，以波行面代替波前面；一方面通过最小走时三角单元的与其相邻三角单元间的“运算”传递扩展波行面，另一方面当遇回转波时，也通过三角单元的相邻三角单元之间的传递将先前已经“运算暂停”状态的三角单元激活为“运算” 状态，而再次成为最小走时三角单元，无需另行“收缩”；这样就有效的解决了波的扩展和回传的问题。 

第二步：波传播路径的向源拾取； 

任意点到源点的射线路径是通过节点最小走时和次级源信息向源追索完成的；图2表示节点E的次级源F的检索；根据E的位置，有以下几类节点的次级源检索： 

（1）当节点E位于三角单元内部时，如图2(a)，通过双曲线近似和最小走时搜索从三角单元三边上找到次级源，算出该节点走时； 

（2）当节点E是三角顶点或三角边的插入点时（图2(b)），可直接拾取次级源； 

（3）当节点E位于三角边插入点A、C之间时，检索E的次级源F要依据A的次级源B与C的次级源D。这时主要会遇到如图2(c)、(d)两种情况：(c)情形是A的次级源B与C的次级源D在一个三角单元内；(d)情形就比较复杂，由于C点是三角顶点，C的次级源D可能存在于其它三角单元内。对图2(c)、(d)两种情况，也通过双曲线近似和最小走时搜索从三角单元相应边上寻找次级源，并算出该节点走时。 

Claims (1)

1.一种混合网最小走时射线追踪层析成像方法，是采用计算机程序实现的，其特征在于，它是将二维复杂外部几何边界及内部复杂空间结构的成像区域，根据其几何结构分解成规则矩形子域和不规则子域，对规则矩形子域采用矩形网格剖分，而对不规则区域则采用三角网格剖分，从而将整个区域网格剖分成包括矩形网和三角网的混合网；采用混合网波行面扩展法计算声波初至走时，追踪出任意接收点至激发点处的射线路径，射线路径包含了从接收点至源点的坐标、所在网格单元信息；将射线路径信息形成矩阵方程，则区域层析成像则归结为求解大型病态线性方程组；再采用外推序列加速代数重构技术求解矩阵方程，得到模型的速度分布，即最终形成规则的网格成像数据；详细步骤如下：

第1步、数据采集，即拾取各震源至接收换能器的初至走时t_m；

第2步：建立成像区域几何结构；

根据边界几何特征、先验信息等建立成像区域的背景网，包括初始网格密度的设置，为混合网剖分打下基础，同时建立成像区域的相三维对坐标系统；

第3步：建立成像区域初至走时信息；

初至走时信息包括：每一震源点坐标、对应每一震源点接收排列的三维坐标和走时t_m；

第4步：将第2步成像区域的相对坐标系统进行旋转变换和切面投影，以利于二维成像操作；

第5步：建立混合型区域网格；

首先对二维复杂外部几何边界的地学层析成像区域，以及内部复杂空间结构的地学层析成像区域的几何结构进行分解，分解成规则矩形子域和不规则子域，对规则矩形子域采用矩形网格剖分，对不规则区域则采用三角网格剖分，从而使整个成像区域网格成为矩形网和三角网混合组成的区域网格；

按背景网及当前速度模型设置的网格密度函数值进行剖分，形成网格单元节点坐标数据文件；

第6步：建立网格各节点的拓扑关系，为射线追踪打下基础；

上述的网格各节点包括：源节点、接收节点、网格单元顶点节点、网格边插入节点；

第7步：初始速度模型输入；

当已知先验速度，应首先输入先验速度信息建立初始速度模型；否则，按均匀速度信息建立初始速度模型；

第8步：对每个源点进行射线追踪，得到相应源点下各节点的次级源及次级源所在的网格单元；射线追踪步骤如下：

第8.1步：波行面扩展

从震源或一个相对最小走时波前点出发，确定这一节点的当前最小走时单元；

先计算当前最小走时单元节点的最小旅行走时，再计算该单元相邻单元节点旅行走时；将当前单元和相邻单元节点的旅行走时相比较，如果旅行走时无变化时，对当前这一最小走时单元作出休眠标志，再寻找下一个最小走时单元作为当前单元，以此类推，遍历所有单元；

在上述遍历所有单元过程中，把遍历单元作为一个波行平面，当该波行平面遇回转波时，先前已经休眠的单元将通过单元之间的传递而激活，再次成为当前最小走时单元；当所有波行面的节点走时都不再变化、且波行平面中所有的单元均处于休眠状态时，波行面扩展完毕；

第8.2步：利用计算出的各节点的旅行时和方向信息确定射线路径，从接收点开始，根据各节点走时及次级源信息，进行向源检索，拾取从接收点到源点的射线路径；

第9步：设置成像反演迭代的终止条件，终止条件可以是迭代次数，也可以是预先设置走时残差；

第10步：接收点的向源检索；

对每个源点对应的接收点进行向源检索，得到射线路径，射线路径包含了从接收点至源点的坐标和所在网格单元信息，所在网格单元信息包括网格单元的编号、次级源节点坐标和射线长度，据此形成Jacobi矩阵A；

由于Jacobi矩阵A为一稀疏矩阵，采用压缩存储，只需记录射线所经历的网格单元及在此网格单元的长度即可；

判断成像反演迭代终止条件，若满足终止条件，则转下一步；否则转第8步；

若想根据射线密度、速度分布信息重新对成像区域离散化时，则转第5步；如此反复进行，直至满足设置的终止条件；

第11步：反演求解

根据第10步形成的Jacobi矩阵A，计算理论走时t_c和走时残差，走时残差δt＝t_m-t_c；建立反演方程Aδ_v＝δt；采用外推加速代数重建反演算法解方程求取δv，修改模型v＝v₀+δv；终止条件是迭代次数或预先设置的走时残差；

第12步：成像结果输出

对于规则网格成像结果输出时，只需输出网格各单元中心坐标及单元的速度值即可；

对于不规则网格成像结果输出时，采用采样输出的办法；即在不规则网格单元的横向和纵向上等间距的设置多个采样点，输出采样点的坐标及单元的速度值，最终形成规则网格成像数据。

Images