线性走时插值射线追踪算法的改进

卢江波　方志

摘要：在LTI（Linear Traveltime Interpolation）射线追踪算法基础上提出的扩张收缩扫描算法能正确追踪直达波、绕射波和回波的射线路径，但其存在计算效率低、收敛速度慢的问题.采用交叉扫描方式对扩张收缩扫描算法进行改进，并由此提出了基于交叉扫描方式的扩张收缩扫描改进算法.理论分析及数值模拟结果表明：改进算法在保留了原扩张收缩扫描算法所有优点的同时，具有更高的计算效率；当模型网格尺寸划分较细时，改进算法在计算效率上的优势更为显著.

关键词：射线追踪；LTI算法；交叉扫描；扩张收缩扫描法；改进算法

中图分类号：TU317 文献标识码：A

射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用.目前常用射线追踪方法主要有有限差分解程函方程法[1-2]、最短路径法[3-4]以及LTI（Linear Travel-time Interpolation）射线追踪算法[5-6]等.实验表明[5]，LTI算法在走时计算以及射线路径追踪上比其它方法（如有限差分解程函方程法）更为快速、精确.

LTI射线追踪算法由Asakawa等人提出[5]，该算法以走时线性变化为前提，分两步进行.第一步，向前计算最小走时：先将模型划分为若干规则的单元，并将各单元边界划分为若干节段，然后，根据走时线性变化的假定以及走时最小原理（Fermat原理），求出从发射点到接收点的最小走时；第二步，根据得到的最小走时，反向追踪射线路径.但是该算法在向前计算走时时，没有考虑射线的逆向传播，不能追踪回波[9]，进而影响射线追踪的精度.针对这一问题，不少学者提出了改进算法，张东等人提出了循环计算LTI改进算法[7]，王浩全提出了交叉扫描LTI改进算法[8]，黄靓等人提出了扩张收缩扫描改进算法[9-10].在考虑射线的逆向传播时，文献[7-8]采用从发射点所在列向外逐列逆向扫描的方式，文献[9-10]则采用从边界列（行）向发射点所在列（行）进行逆向扫描的方式.在考虑射线的逆向传播上，文献[9-10]提出的改进算法更为合理，但是该算法存在计算效率低、收敛速度慢的问题.因为在LTI算法中，接收点要得到其最小走时，需其理论射线路径与模型单元边界所有相交节点都已得到其最小走时，这些相交节点的最小走时可能来自列扫描，也可能来自行扫描.根据本文作者对文献[5，9-10]算法步骤的理解以及对文献[9-10]中数值算例的研究，发现文献[9-10]不仅增加了从边界进行逆向传播的射线，即文献[9-10]中的收缩扫描过程，而且将文献[5]中既考虑邻列也考虑来自邻行入射射线的逐列扫描方式改为逐列（行）扫描只考虑邻列（行）入射射线的方式，这点可由文献[9]算例模型1的计算结果得出：若文献[9]的逐列（行）扫描过程既考虑邻列又考虑邻行的入射射线，则模型1经一次扩张扫描即得接收点的最小走时，此时相对误差应为0.247‰，而不应为按行列分开扫描方式经一次扩张扫描得到的计算结果3.3‰注：事实上，模型1接收点在按行列分开扫描的方式下，经扩张收缩再扩张后能够得到其最小走时，文献[9]给出的是一次扩张扫描结果，这应该是文献[9]作者的疏忽..文献[9]采用行列分开扫描的方式使得文献[9-10]在能正确追踪回波的同时，接收点为获取其最小走时也进行了较多的无效扫描，降低了算法的计算效率.由以上分析可知，若采用行列交叉扫描的方式对扩张收缩扫描算法进行改进，将有效提高扩张收缩扫描算法的计算效率.

同时，由于绕射波以及回波的存在，文献[5]算法步骤3对每列都进行水平边界节点最小走时搜索，其意义并不明确.此外，由于扩张收缩扫描算法相比文献[5]增加了收缩扫描过程，在逐列扫描过程中能够考虑上行或下行首波的最小走时，所以文献[5]中的逐行扫描可以省略.最后，由于存在收缩扫描过程，在计算竖向边界各节点最小走时时，若按文献[5]算法步骤5的计算方法，则存在着重复无效扫描.因此，也有必要对文献[5]的扩张扫描过程进行简化和改进.

基于此，本文在前人研究的基础上提出了基于交叉扫描方式的扩张收缩扫描新算法，该算法以扩张收缩扫描算法为基础，结合文献[5]在逐列扫描过程中进行交叉扫描的思想，改进扩张收缩扫描算法中逐列扫描的具体计算方法，提高了算法的计算效率，减少了迭代次数；同时，保留了原扩张收缩扫描算法的所有优点.

2）计算发射点S所在单元边界上各节点的走时，并记录次级源.假定此模型在各单元的边界上均只划分两个节段（图5），则各单元均有8个节点.根据单元节点以及发射点在整体坐标系中的坐标，可求得发射点S所在单元边界上各节点的走时，并将S点记为各节点的次级源.

以I点的走时计算为例（仅考虑通过下边界GE到达I节点的射线），易知，满足I点最小走时要求的射线可能来自GE中的任一节段.此时，根据第1节给出的计算公式和计算方法分别计算出射线通过GF节段和FE节段时I点的最小走时，取两个最小走时的较小值作为I点的最小走时，并记录相应的次级源.以同样步骤求出EJLG单元其它节点的走时，完成该单元的计算.这种通过单元下边界节点走时计算其它节点走时的过程，称为向上扫描.其它通过上边界、左边界、右边界的情况分别称为向下扫描、向右扫描以及向左扫描.

算例3的计算结果表明：改进算法保留了扩张收缩扫描算法能正确处理射线逆向传播的优点，并且具有更快的收敛速度.

5结语

理论分析以及数值算例表明，本文提出的基于交叉扫描方式的扩张收缩扫描算法，不仅具有原扩张收缩扫描算法的所有优点，而且在不增加单次扩张扫描计算量的前提下，通过改变算法的扫描顺序，提高了算法的计算效率，加快了算法的收敛速度.特别当模型网格划分数较多时，改进算法在计算效率方面的优势更为显著.

参考文献

[1]VIDALE J. Finitedifference calculationoftravel times[J].Bulletin of the Seismological Society of America， 1988， 78（6）： 2062-2076.

[2]QIN F， LUO Y， OLSEN K B， et al. Finitedifferencesolution of the eikonal equation alongexpandingwavefronts[J].Geophysics， 1992，57（3）：478-487.

[3]MOSERT J. Shortest path calculation of seismic rays[J]. Geophysics， 1991， 56（1）： 59-67.

[4]刘洪，孟凡林，李幼铭.计算最小走时和射线路径的界面网全局方法[J]. 地球物理学报，1995，38（6）：823-832.

[5]ASAKAWA E， KAWANAKA T. Seismic ray tracing using linear traveltimeinterpolation[J]. Geophysical Prospecting， 1993， 41（1）： 99-111.

[6]赵改善，郝守玲，杨尔皓，等. 基于旅行时线性插值的地震射线追踪算法[J]. 石油物探， 1998， 37（2）： 14-24.

[7]张东，谢宝莲，杨艳，等. 一种改进的线性走时插值射线追踪算法[J]. 地球物理学报，2009， 52（1）： 200-205.

[8]WANG Haoquan. An improved method of linear traveltime interpolation ray tracing algorithm[J].Acta Physica PolonicaSeries A， 2010， 118（4）： 521.

[9]黄靓，黄政宇. 线性插值射线追踪的改进方法[J]. 湘潭大学学报：自然科学版， 2002， 24（4）： 105-108.

[10]黄靓. 混凝土超声波层析成像的理论方法和试验研究[D]. 长沙： 湖南大学土木工程学院， 2008：33-35.

摘要：在LTI（Linear Traveltime Interpolation）射线追踪算法基础上提出的扩张收缩扫描算法能正确追踪直达波、绕射波和回波的射线路径，但其存在计算效率低、收敛速度慢的问题.采用交叉扫描方式对扩张收缩扫描算法进行改进，并由此提出了基于交叉扫描方式的扩张收缩扫描改进算法.理论分析及数值模拟结果表明：改进算法在保留了原扩张收缩扫描算法所有优点的同时，具有更高的计算效率；当模型网格尺寸划分较细时，改进算法在计算效率上的优势更为显著.

关键词：射线追踪；LTI算法；交叉扫描；扩张收缩扫描法；改进算法

中图分类号：TU317 文献标识码：A

射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用.目前常用射线追踪方法主要有有限差分解程函方程法[1-2]、最短路径法[3-4]以及LTI（Linear Travel-time Interpolation）射线追踪算法[5-6]等.实验表明[5]，LTI算法在走时计算以及射线路径追踪上比其它方法（如有限差分解程函方程法）更为快速、精确.

LTI射线追踪算法由Asakawa等人提出[5]，该算法以走时线性变化为前提，分两步进行.第一步，向前计算最小走时：先将模型划分为若干规则的单元，并将各单元边界划分为若干节段，然后，根据走时线性变化的假定以及走时最小原理（Fermat原理），求出从发射点到接收点的最小走时；第二步，根据得到的最小走时，反向追踪射线路径.但是该算法在向前计算走时时，没有考虑射线的逆向传播，不能追踪回波[9]，进而影响射线追踪的精度.针对这一问题，不少学者提出了改进算法，张东等人提出了循环计算LTI改进算法[7]，王浩全提出了交叉扫描LTI改进算法[8]，黄靓等人提出了扩张收缩扫描改进算法[9-10].在考虑射线的逆向传播时，文献[7-8]采用从发射点所在列向外逐列逆向扫描的方式，文献[9-10]则采用从边界列（行）向发射点所在列（行）进行逆向扫描的方式.在考虑射线的逆向传播上，文献[9-10]提出的改进算法更为合理，但是该算法存在计算效率低、收敛速度慢的问题.因为在LTI算法中，接收点要得到其最小走时，需其理论射线路径与模型单元边界所有相交节点都已得到其最小走时，这些相交节点的最小走时可能来自列扫描，也可能来自行扫描.根据本文作者对文献[5，9-10]算法步骤的理解以及对文献[9-10]中数值算例的研究，发现文献[9-10]不仅增加了从边界进行逆向传播的射线，即文献[9-10]中的收缩扫描过程，而且将文献[5]中既考虑邻列也考虑来自邻行入射射线的逐列扫描方式改为逐列（行）扫描只考虑邻列（行）入射射线的方式，这点可由文献[9]算例模型1的计算结果得出：若文献[9]的逐列（行）扫描过程既考虑邻列又考虑邻行的入射射线，则模型1经一次扩张扫描即得接收点的最小走时，此时相对误差应为0.247‰，而不应为按行列分开扫描方式经一次扩张扫描得到的计算结果3.3‰注：事实上，模型1接收点在按行列分开扫描的方式下，经扩张收缩再扩张后能够得到其最小走时，文献[9]给出的是一次扩张扫描结果，这应该是文献[9]作者的疏忽..文献[9]采用行列分开扫描的方式使得文献[9-10]在能正确追踪回波的同时，接收点为获取其最小走时也进行了较多的无效扫描，降低了算法的计算效率.由以上分析可知，若采用行列交叉扫描的方式对扩张收缩扫描算法进行改进，将有效提高扩张收缩扫描算法的计算效率.

同时，由于绕射波以及回波的存在，文献[5]算法步骤3对每列都进行水平边界节点最小走时搜索，其意义并不明确.此外，由于扩张收缩扫描算法相比文献[5]增加了收缩扫描过程，在逐列扫描过程中能够考虑上行或下行首波的最小走时，所以文献[5]中的逐行扫描可以省略.最后，由于存在收缩扫描过程，在计算竖向边界各节点最小走时时，若按文献[5]算法步骤5的计算方法，则存在着重复无效扫描.因此，也有必要对文献[5]的扩张扫描过程进行简化和改进.

基于此，本文在前人研究的基础上提出了基于交叉扫描方式的扩张收缩扫描新算法，该算法以扩张收缩扫描算法为基础，结合文献[5]在逐列扫描过程中进行交叉扫描的思想，改进扩张收缩扫描算法中逐列扫描的具体计算方法，提高了算法的计算效率，减少了迭代次数；同时，保留了原扩张收缩扫描算法的所有优点.

2）计算发射点S所在单元边界上各节点的走时，并记录次级源.假定此模型在各单元的边界上均只划分两个节段（图5），则各单元均有8个节点.根据单元节点以及发射点在整体坐标系中的坐标，可求得发射点S所在单元边界上各节点的走时，并将S点记为各节点的次级源.

以I点的走时计算为例（仅考虑通过下边界GE到达I节点的射线），易知，满足I点最小走时要求的射线可能来自GE中的任一节段.此时，根据第1节给出的计算公式和计算方法分别计算出射线通过GF节段和FE节段时I点的最小走时，取两个最小走时的较小值作为I点的最小走时，并记录相应的次级源.以同样步骤求出EJLG单元其它节点的走时，完成该单元的计算.这种通过单元下边界节点走时计算其它节点走时的过程，称为向上扫描.其它通过上边界、左边界、右边界的情况分别称为向下扫描、向右扫描以及向左扫描.

算例3的计算结果表明：改进算法保留了扩张收缩扫描算法能正确处理射线逆向传播的优点，并且具有更快的收敛速度.

5结语

理论分析以及数值算例表明，本文提出的基于交叉扫描方式的扩张收缩扫描算法，不仅具有原扩张收缩扫描算法的所有优点，而且在不增加单次扩张扫描计算量的前提下，通过改变算法的扫描顺序，提高了算法的计算效率，加快了算法的收敛速度.特别当模型网格划分数较多时，改进算法在计算效率方面的优势更为显著.

参考文献

[1]VIDALE J. Finitedifference calculationoftravel times[J].Bulletin of the Seismological Society of America， 1988， 78（6）： 2062-2076.

[2]QIN F， LUO Y， OLSEN K B， et al. Finitedifferencesolution of the eikonal equation alongexpandingwavefronts[J].Geophysics， 1992，57（3）：478-487.

[3]MOSERT J. Shortest path calculation of seismic rays[J]. Geophysics， 1991， 56（1）： 59-67.

[4]刘洪，孟凡林，李幼铭.计算最小走时和射线路径的界面网全局方法[J]. 地球物理学报，1995，38（6）：823-832.

[5]ASAKAWA E， KAWANAKA T. Seismic ray tracing using linear traveltimeinterpolation[J]. Geophysical Prospecting， 1993， 41（1）： 99-111.

[6]赵改善，郝守玲，杨尔皓，等. 基于旅行时线性插值的地震射线追踪算法[J]. 石油物探， 1998， 37（2）： 14-24.

[7]张东，谢宝莲，杨艳，等. 一种改进的线性走时插值射线追踪算法[J]. 地球物理学报，2009， 52（1）： 200-205.

[8]WANG Haoquan. An improved method of linear traveltime interpolation ray tracing algorithm[J].Acta Physica PolonicaSeries A， 2010， 118（4）： 521.

[9]黄靓，黄政宇. 线性插值射线追踪的改进方法[J]. 湘潭大学学报：自然科学版， 2002， 24（4）： 105-108.

[10]黄靓. 混凝土超声波层析成像的理论方法和试验研究[D]. 长沙： 湖南大学土木工程学院， 2008：33-35.

摘要：在LTI（Linear Traveltime Interpolation）射线追踪算法基础上提出的扩张收缩扫描算法能正确追踪直达波、绕射波和回波的射线路径，但其存在计算效率低、收敛速度慢的问题.采用交叉扫描方式对扩张收缩扫描算法进行改进，并由此提出了基于交叉扫描方式的扩张收缩扫描改进算法.理论分析及数值模拟结果表明：改进算法在保留了原扩张收缩扫描算法所有优点的同时，具有更高的计算效率；当模型网格尺寸划分较细时，改进算法在计算效率上的优势更为显著.

关键词：射线追踪；LTI算法；交叉扫描；扩张收缩扫描法；改进算法

中图分类号：TU317 文献标识码：A

射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用.目前常用射线追踪方法主要有有限差分解程函方程法[1-2]、最短路径法[3-4]以及LTI（Linear Travel-time Interpolation）射线追踪算法[5-6]等.实验表明[5]，LTI算法在走时计算以及射线路径追踪上比其它方法（如有限差分解程函方程法）更为快速、精确.

LTI射线追踪算法由Asakawa等人提出[5]，该算法以走时线性变化为前提，分两步进行.第一步，向前计算最小走时：先将模型划分为若干规则的单元，并将各单元边界划分为若干节段，然后，根据走时线性变化的假定以及走时最小原理（Fermat原理），求出从发射点到接收点的最小走时；第二步，根据得到的最小走时，反向追踪射线路径.但是该算法在向前计算走时时，没有考虑射线的逆向传播，不能追踪回波[9]，进而影响射线追踪的精度.针对这一问题，不少学者提出了改进算法，张东等人提出了循环计算LTI改进算法[7]，王浩全提出了交叉扫描LTI改进算法[8]，黄靓等人提出了扩张收缩扫描改进算法[9-10].在考虑射线的逆向传播时，文献[7-8]采用从发射点所在列向外逐列逆向扫描的方式，文献[9-10]则采用从边界列（行）向发射点所在列（行）进行逆向扫描的方式.在考虑射线的逆向传播上，文献[9-10]提出的改进算法更为合理，但是该算法存在计算效率低、收敛速度慢的问题.因为在LTI算法中，接收点要得到其最小走时，需其理论射线路径与模型单元边界所有相交节点都已得到其最小走时，这些相交节点的最小走时可能来自列扫描，也可能来自行扫描.根据本文作者对文献[5，9-10]算法步骤的理解以及对文献[9-10]中数值算例的研究，发现文献[9-10]不仅增加了从边界进行逆向传播的射线，即文献[9-10]中的收缩扫描过程，而且将文献[5]中既考虑邻列也考虑来自邻行入射射线的逐列扫描方式改为逐列（行）扫描只考虑邻列（行）入射射线的方式，这点可由文献[9]算例模型1的计算结果得出：若文献[9]的逐列（行）扫描过程既考虑邻列又考虑邻行的入射射线，则模型1经一次扩张扫描即得接收点的最小走时，此时相对误差应为0.247‰，而不应为按行列分开扫描方式经一次扩张扫描得到的计算结果3.3‰注：事实上，模型1接收点在按行列分开扫描的方式下，经扩张收缩再扩张后能够得到其最小走时，文献[9]给出的是一次扩张扫描结果，这应该是文献[9]作者的疏忽..文献[9]采用行列分开扫描的方式使得文献[9-10]在能正确追踪回波的同时，接收点为获取其最小走时也进行了较多的无效扫描，降低了算法的计算效率.由以上分析可知，若采用行列交叉扫描的方式对扩张收缩扫描算法进行改进，将有效提高扩张收缩扫描算法的计算效率.

同时，由于绕射波以及回波的存在，文献[5]算法步骤3对每列都进行水平边界节点最小走时搜索，其意义并不明确.此外，由于扩张收缩扫描算法相比文献[5]增加了收缩扫描过程，在逐列扫描过程中能够考虑上行或下行首波的最小走时，所以文献[5]中的逐行扫描可以省略.最后，由于存在收缩扫描过程，在计算竖向边界各节点最小走时时，若按文献[5]算法步骤5的计算方法，则存在着重复无效扫描.因此，也有必要对文献[5]的扩张扫描过程进行简化和改进.

基于此，本文在前人研究的基础上提出了基于交叉扫描方式的扩张收缩扫描新算法，该算法以扩张收缩扫描算法为基础，结合文献[5]在逐列扫描过程中进行交叉扫描的思想，改进扩张收缩扫描算法中逐列扫描的具体计算方法，提高了算法的计算效率，减少了迭代次数；同时，保留了原扩张收缩扫描算法的所有优点.

2）计算发射点S所在单元边界上各节点的走时，并记录次级源.假定此模型在各单元的边界上均只划分两个节段（图5），则各单元均有8个节点.根据单元节点以及发射点在整体坐标系中的坐标，可求得发射点S所在单元边界上各节点的走时，并将S点记为各节点的次级源.

以I点的走时计算为例（仅考虑通过下边界GE到达I节点的射线），易知，满足I点最小走时要求的射线可能来自GE中的任一节段.此时，根据第1节给出的计算公式和计算方法分别计算出射线通过GF节段和FE节段时I点的最小走时，取两个最小走时的较小值作为I点的最小走时，并记录相应的次级源.以同样步骤求出EJLG单元其它节点的走时，完成该单元的计算.这种通过单元下边界节点走时计算其它节点走时的过程，称为向上扫描.其它通过上边界、左边界、右边界的情况分别称为向下扫描、向右扫描以及向左扫描.

算例3的计算结果表明：改进算法保留了扩张收缩扫描算法能正确处理射线逆向传播的优点，并且具有更快的收敛速度.

5结语

理论分析以及数值算例表明，本文提出的基于交叉扫描方式的扩张收缩扫描算法，不仅具有原扩张收缩扫描算法的所有优点，而且在不增加单次扩张扫描计算量的前提下，通过改变算法的扫描顺序，提高了算法的计算效率，加快了算法的收敛速度.特别当模型网格划分数较多时，改进算法在计算效率方面的优势更为显著.

参考文献

[1]VIDALE J. Finitedifference calculationoftravel times[J].Bulletin of the Seismological Society of America， 1988， 78（6）： 2062-2076.

[2]QIN F， LUO Y， OLSEN K B， et al. Finitedifferencesolution of the eikonal equation alongexpandingwavefronts[J].Geophysics， 1992，57（3）：478-487.

[3]MOSERT J. Shortest path calculation of seismic rays[J]. Geophysics， 1991， 56（1）： 59-67.

[4]刘洪，孟凡林，李幼铭.计算最小走时和射线路径的界面网全局方法[J]. 地球物理学报，1995，38（6）：823-832.

[5]ASAKAWA E， KAWANAKA T. Seismic ray tracing using linear traveltimeinterpolation[J]. Geophysical Prospecting， 1993， 41（1）： 99-111.

[6]赵改善，郝守玲，杨尔皓，等. 基于旅行时线性插值的地震射线追踪算法[J]. 石油物探， 1998， 37（2）： 14-24.

[7]张东，谢宝莲，杨艳，等. 一种改进的线性走时插值射线追踪算法[J]. 地球物理学报，2009， 52（1）： 200-205.

[8]WANG Haoquan. An improved method of linear traveltime interpolation ray tracing algorithm[J].Acta Physica PolonicaSeries A， 2010， 118（4）： 521.

[9]黄靓，黄政宇. 线性插值射线追踪的改进方法[J]. 湘潭大学学报：自然科学版， 2002， 24（4）： 105-108.

[10]黄靓. 混凝土超声波层析成像的理论方法和试验研究[D]. 长沙： 湖南大学土木工程学院， 2008：33-35.