非线性射线走时层析成像方法、电子设备及介质与流程

1.本发明涉及地球物理勘探中的地震走时层析成像领域，更具体地，涉及一种非线性射线走时层析成像方法、电子设备及介质。


背景技术：

2.地震走时层析的目的是找到一个使理论合成走时与观测走时之差在某种意义下达到极小的“最佳”模型。目前基于线性反演思想的方法主要有两种实现方式。第一种是传统的求解层析方程组的方法。这种方法求解的方程组系数矩阵规模庞大，且非常稀疏，病态严重，需要采用正则化来保证反演的正确性与稳定性，因此反演精度严重依赖于方程组的求解方法。第二种是基于局部最优化思想，通过对目标函数的局部线性化近似获得速度的修改方向，然后通过迭代实现该非线性反演问题的线性化求解。但无论哪种实现方式，线性反演方法的主要缺点是对初始表层速度模型依赖性强，当初始表层速度模型选择不合适时，很容易陷入局部极值。尽管如此，这些方法在实际应用中仍可以取得不错的效果。
3.在第二种方法中，每一次迭代中速度修改方向的计算是关键。利用伴随状态法或散射积分法可以获得速度的修改方向。前者利用地震波的一次正传和数据残差的一次反传即可获得目标函数的梯度，无需存储大型矩阵，但其预条件在射线走时层析成像中不易实现。后者通过显式地计算核函数，并与走时残差向量相乘实现梯度的计算。但与传统的层析成像方法一样，散射积分法面临内存占用大的问题，特别是在利用hessian矩阵时。
4.hessian矩阵的利用在局部最优化反演中对提高反演精度与多参数解耦有非常重要的作用。核函数表达的hessian矩阵包含两项，利用准确的hessian矩阵构造的方向称为牛顿方向，只保留hessian第一项的矩阵称为近似hessian矩阵，相对应的方向为高斯-牛顿方向。实际应用中，无论是牛顿方向还是高斯-牛顿方向都涉及到hessian矩阵求逆的操作，因此往往难以准确构造。更通常的做法是利用近似hessian矩阵的对角阵来构造hessian矩阵的逆，所对应的方向是拟牛顿方向的一种。由于hessian矩阵的逆是直接作用在梯度上面，因此这种做法又称为梯度方向的预条件。预条件的实施与否对层析有很大的影响。
5.现有技术在波形反演研究中提出了一种改进的散射积分算法。该方法将大规模的核函数-向量乘运算表示为具有明确物理含义的向量-标量乘的累加运算，因此避免了大规模矩阵的存储，使得传统散射积分法波形反演更加实用化。同时，该方法可以方便地实现hessian矩阵主对角线元素的计算，进而可以方便实现梯度的预条件。但是，与传统的层析成像方法一样，散射积分法也面临内存占用大的问题，特别是在利用hessian矩阵时。
6.因此，有必要开发一种非线性射线走时层析成像方法、电子设备及介质。
7.公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。


技术实现要素：

8.本发明提出了一种非线性射线走时层析成像方法、电子设备及介质，其能够避免传统层析成像中的大规模矩阵的存储和预条件实现不方便的问题，进而实现带预条件的最速下降法地震走时层析成像。
9.第一方面，本公开实施例提供了一种非线性射线走时层析成像方法，包括：
10.步骤1：对原始地震数据进行去噪和反褶积的预处理，获得预处理地震数据；
11.步骤2：针对所述预处理地震数据进行初至拾取，获得初至波到达时；
12.步骤3：根据地下先验信息建立初始表层速度模型，设定反演参数；
13.步骤4：根据初至波到达时与所述初始表层速度模型通过散射积分法进行非线性射线走时层析成像反演，获得最终表层速度模型并成像。
14.优选地，所述步骤4包括：
15.步骤41：在所述初始表层速度模型上进行射线追踪，获取理论合成初至走时与射线路径，保存所述射线路径；
16.步骤42：判断所述理论合成初至走时与所述初至波到达时是否匹配，若是，则以所述初始表层速度模型作为所述最终表层速度模型，若否，则执行步骤43；
17.步骤43：逐个读取每一个炮检对对应的射线路径作为该炮检对对应的核函数；
18.步骤44：通过核函数-标量乘获得每一个炮检对所对应的梯度，进而计算整个观测系统的全局梯度，并计算所述全局梯度的预条件算子；
19.步骤45：计算预条件最速下降方向、方向更新步长与模型更新量；
20.步骤46：更新表层速度模型，以更新后的表层速度模型为初始表层速度模型，重复步骤41-46，直到获得所述最终表层速度模型。
21.优选地，所述步骤41包括：
22.步骤411：在所述初始表层速度模型网格剖分的网格边界上设定多个关键点；
23.步骤412：自激发点开始以矩形为波前，按惠更斯原理逐层向外扩展，直至扫描完成整个模型空间，再由外向里收缩，直至扫描到激发点，如此重复直至模型空间的最小旅行时场不变为止；
24.步骤413：检波点处的初至波到达时即为所述理论合成初至走时。
25.优选地，所述步骤44包括：
26.步骤441：获取每一个炮检对所对应的射线路径，进而获得所述射线路径与该炮检对对应的走时残差乘，得到该炮检对对应的梯度以及射线路径向量的元素平方向量；
27.步骤442：将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度；
28.步骤443：将所有的炮检对的核函数元素取平方并进行累加，计算全局梯度的预条件算子。
29.优选地，通过公式(1)将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度：
[0030][0031]
其中，k
t
δt为全局梯度，k
ij
为第i对炮检对对应的核函数在第j个模型网格内的核函数值，δti为旅行时差，k
ij
δti为第i对炮检对的梯度，k
t
为影响走时的核函数矩阵、δt为旅行时差。
[0032]
优选地，通过公式(2)计算全局梯度的预条件算子：
[0033][0034]
其中，h0为预条件算子。
[0035]
优选地，所述步骤45中，通过公式(3)计算预条件最速下降方向：
[0036][0037]
其中，p为预条件最速下降方向。
[0038]
优选地，所述步骤45中方向更新步长的计算方法为：
[0039]
根据先验信息获取所述初始表层速度模型与真实地下模型的速度差，将速度差除以预先设定的迭代次数，得到每次迭代可更新的最大速度值δv
max
，利用待定系数法求取所述方向更新步长t，使得max{|p
·
t|}＝δv
max
。
[0040]
作为本公开实施例的一种具体实现方式，
[0041]
第二方面，本公开实施例还提供了一种电子设备，该电子设备包括：
[0042]
存储器，存储有可执行指令；
[0043]
处理器，所述处理器运行所述存储器中的所述可执行指令，以实现所述的非线性射线走时层析成像方法。
[0044]
第三方面，本公开实施例还提供了一种计算机可读存储介质，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现所述的非线性射线走时层析成像方法。
[0045]
其有益效果在于：
[0046]
(1)无需存储核函数矩阵即可方便地实现目标函数梯度的计算，在任何时刻只需存储单个核函数—射线路径；
[0047]
(2)计算效率高，易于并行：将大规模的核函数-向量乘运算表示为具有明确物理含义的向量-标量乘的累加运算，无需利用svd或者lsqr方法解大规模矩阵方程组，计算时间大大减少；
[0048]
(3)方便预条件，使得计算效率进一步提高，反演精度也大大提高：该方法无需存储hessian矩阵即可方便实现预条件，预条件算子的应用可以加快收敛速度，并能够获得地
下深部介质的信息；
[0049]
(4)反演过程更稳定：优化的梯度导引算法和预条件算子的使用使得反演过程更稳定。
[0050]
本发明的方法和装置具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
[0051]
通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。
[0052]
图1示出了根据本发明的一个实施例的非线性射线走时层析成像方法的步骤的流程图。
[0053]
图2示出了根据本发明的实施例1的二维复杂起伏地表理论模型的示意图。
[0054]
图3示出了根据本发明的实施例1的初始表层速度模型的示意图。
[0055]
图4a和图4b分别示出了根据本发明的实施例1的水平位置9km和26km处的两个z分量炮记录的示意图。
[0056]
图5a和图5b分别示出了根据本发明的实施例1的本方法与sirt法层析经过40次迭代后的反演结果的示意图。
[0057]
图6a和图6b分别示出了根据本发明的实施例1的地下100m、200m处的速度曲线的示意图。
[0058]
图7示出了根据本发明的实施例2的初始表层速度模型的示意图。
[0059]
图8示出了根据本发明的实施例2的最终表层速度模型的示意图。
[0060]
图9示出了根据本发明的实施例2的层析静校正叠加剖面的示意图。
具体实施方式
[0061]
下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。
[0062]
本发明提供一种非线性射线走时层析成像方法，包括：
[0063]
步骤1：对原始地震数据进行去噪和反褶积的预处理，获得预处理地震数据；
[0064]
步骤2：针对预处理地震数据进行初至拾取，获得初至波到达时；
[0065]
步骤3：根据地下先验信息建立初始表层速度模型，设定反演参数；
[0066]
步骤4：根据初至波到达时与初始表层速度模型通过散射积分法进行非线性射线走时层析成像反演，获得最终表层速度模型并成像。
[0067]
在一个示例中，步骤4包括：
[0068]
步骤41：在初始表层速度模型上进行射线追踪，获取理论合成初至走时与射线路径，保存射线路径；
[0069]
步骤42：判断理论合成初至走时与初至波到达时是否匹配，若是，则以初始表层速度模型作为最终表层速度模型，若否，则执行步骤43；
[0070]
步骤43：逐个读取每一个炮检对对应的射线路径作为该炮检对对应的核函数；
[0071]
步骤44：通过核函数-标量乘获得每一个炮检对所对应的梯度，进而计算整个观测系统的全局梯度，并计算全局梯度的预条件算子；
[0072]
步骤45：计算预条件最速下降方向、方向更新步长与模型更新量；
[0073]
步骤46：更新表层速度模型，以更新后的表层速度模型为初始表层速度模型，重复步骤41-46，直到获得最终表层速度模型。
[0074]
在一个示例中，步骤41包括：
[0075]
步骤411：在初始表层速度模型网格剖分的网格边界上设定多个关键点；
[0076]
步骤412：自激发点开始以矩形为波前，按惠更斯原理逐层向外扩展，直至扫描完成整个模型空间，再由外向里收缩，直至扫描到激发点，如此重复直至模型空间的最小旅行时场不变为止；
[0077]
步骤413：检波点处的初至波到达时即为理论合成初至走时。
[0078]
在一个示例中，步骤44包括：
[0079]
步骤441：获取每一个炮检对所对应的射线路径，进而获得射线路径与该炮检对对应的走时残差乘，得到该炮检对对应的梯度以及射线路径向量的元素平方向量；
[0080]
步骤442：将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度；
[0081]
步骤443：将所有的炮检对的核函数元素取平方并进行累加，计算全局梯度的预条件算子。
[0082]
在一个示例中，通过公式(1)将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度：
[0083][0084]
其中，k
t
δt为全局梯度，k
ij
为第i对炮检对对应的核函数在第j个模型网格内的核函数值，δti为旅行时差，k
ij
δti为第i对炮检对的梯度，k
t
为影响走时的核函数矩阵、δt为旅行时差。
[0085]
在一个示例中，通过公式(2)计算全局梯度的预条件算子：
[0086][0087]
其中，h0为预条件算子。
[0088]
在一个示例中，步骤45中，通过公式(3)计算预条件最速下降方向：
[0089][0090]
其中，p为预条件最速下降方向。
[0091]
在一个示例中，步骤45中方向更新步长的计算方法为：
[0092]
根据先验信息获取初始表层速度模型与真实地下模型的速度差，将速度差除以预先设定的迭代次数，得到每次迭代可更新的最大速度值δv
max
，利用待定系数法求取方向更新步长t，使得max{|p
·
t|}＝δv
max
。
[0093]
具体地，步骤1：对原始地震数据进行去噪和反褶积的预处理，获得预处理地震数据；
[0094]
步骤2：针对预处理地震数据进行初至拾取，获得初至波到达时；
[0095]
步骤3：根据地下先验信息建立初始表层速度模型，设定反演参数；
[0096]
步骤4：根据初至波到达时与初始表层速度模型通过散射积分法进行非线性射线走时层析成像反演，获得最终表层速度模型并成像。
[0097]
其中，步骤4包括：
[0098]
步骤41：在初始表层速度模型上进行射线追踪，获取理论合成初至走时与射线路径，保存射线路径：
[0099]
步骤411：在初始表层速度模型网格剖分的网格边界上设定多个关键点；
[0100]
步骤412：自激发点开始以矩形为波前，按惠更斯原理逐层向外扩展，直至扫描完成整个模型空间，再由外向里收缩，直至扫描到激发点，如此重复直至模型空间的最小旅行时场不变为止，在扩展或收缩过程中，次级源来自于两个关键点之间，两个关键点之间的走时通过线性插值得到；
[0101]
步骤413：检波点处的初至波到达时即为理论合成初至走时；
[0102]
步骤42：判断理论合成初至走时与初至波到达时是否匹配，建立关于理论合成初至走时与初至波到达时的二范数形式的目标函数，判断是否满足预定的要求，若是，则以初始表层速度模型作为最终表层速度模型，若否，则执行步骤43；
[0103]
步骤43：逐个读取每一个炮检对对应的射线路径作为该炮检对对应的核函数；
[0104]
步骤44：通过核函数-标量乘获得每一个炮检对所对应的梯度，进而通过公式(1)计算整个观测系统的全局梯度，并通过公式(2)计算全局梯度的预条件算子：
[0105]
步骤441：获取每一个炮检对所对应的射线路径，进而获得射线路径与该炮检对对应的走时残差乘，得到该炮检对对应的梯度以及射线路径向量的元素平方向量；
[0106]
步骤442：将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度；
[0107]
步骤443：将所有的炮检对的核函数元素取平方并进行累加，计算全局梯度的预条件算子；
[0108]
步骤45：通过公式(3)计算预条件最速下降方向，根据先验信息获取初始表层速度模型与真实地下模型的速度差，将速度差除以预先设定的迭代次数，得到每次迭代可更新的最大速度值δv
max
，利用待定系数法求取方向更新步长t，使得max{|p
·
t|}＝δv
max
，并计算模型更新量；
[0109]
步骤46：更新表层速度模型，以更新后的表层速度模型为初始表层速度模型，重复步骤41-46，直到获得最终表层速度模型。
[0110]
本发明还提供一种电子设备，电子设备包括：存储器，存储有可执行指令；处理器，处理器运行存储器中的可执行指令，以实现上述的非线性射线走时层析成像方法。
[0111]
本发明还提供一种计算机可读存储介质，该计算机可读存储介质存储有计算机程
序，该计算机程序被处理器执行时实现上述的非线性射线走时层析成像方法。
[0112]
为便于理解本发明实施例的方案及其效果，以下给出四个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
[0113]
实施例1
[0114]
图1示出了根据本发明的一个实施例的非线性射线走时层析成像方法的步骤的流程图。
[0115]
如图1所示，该非线性射线走时层析成像方法包括：步骤1：对原始地震数据进行去噪和反褶积的预处理，获得预处理地震数据；步骤2：针对预处理地震数据进行初至拾取，获得初至波到达时；步骤3：根据地下先验信息建立初始表层速度模型，设定反演参数；步骤4：根据初至波到达时与初始表层速度模型通过散射积分法进行非线性射线走时层析成像反演，获得最终表层速度模型并成像。
[0116]
图2示出了根据本发明的实施例1的二维复杂起伏地表理论模型的示意图。
[0117]
图3示出了根据本发明的实施例1的初始表层速度模型的示意图。
[0118]
图4a和图4b分别示出了根据本发明的实施例1的水平位置9km和26km处的两个z分量炮记录的示意图。
[0119]
本实施例将以如图2所示的二维复杂起伏地表理论模型作为真实模型，初始表层速度模型如图3所示。该模型地表起伏落差达450m，近地表速度横向变化剧烈。模型被离散为4001
×
151个网格，空间离散间隔为10m
×
10m。利用弹性波方程模拟地表地震记录，并拾取z分量记录初至波走时作为观测数据。共模拟了619炮数据，第一炮位于5km处，炮点间隔40m。检波器均匀、对称地分布在炮点两端，水平间隔为20m。每炮201个检波器，最小炮检距为0m，最大炮检距为2.0km。水平位置9km和26km处的两个z分量炮记录如图4a、图4b所示。
[0120]
具体实施方式如下：
[0121]
(1)数据采集器采集地震波信号，并对原始地震数据进行去噪、反褶积等非破坏走时的预处理；
[0122]
(2)数据采集器将预处理后的地震波数据逐道输入到处理器中，在处理过的地震数据炮记录上进行初至拾取，以获得初至波到达时；
[0123]
(3)输入设备根据地下先验信息建立初始表层速度模型；
[0124]
(4)处理器在当前模型上进行射线追踪以获得理论合成初至走时及射线路径；
[0125]
(5)处理器判断理论合成走时与初至波到达时的匹配程度，如果“吻合”则保留当前模型并退出，并通过显示器显示模型，否则继续执行以下步骤；
[0126]
(6)处理器逐一读取步骤(4)算得的射线路径；
[0127]
(7)处理器通过核函数-标量乘获得每一个炮检对所对应的梯度，并累加形成整个观测系统所对应的梯度，同时计算核函数自身元素的平方并沿着所有炮检对累加求和获得梯度的预条件算子；
[0128]
(8)处理器计算预条件最速下降方向、方向更新步长、模型更新量，并更新模型重复上述步骤，获得最终表层速度模型并成像。
[0129]
图5a和图5b分别示出了根据本发明的实施例1的本方法与sirt法层析经过40次迭代后的反演结果的示意图。
[0130]
图6a和图6b分别示出了根据本发明的实施例1的地下100m、200m处的速度曲线的示意图。
[0131]
为了对比效果，散射积分法射线层析与传统sirt法层析经过40次迭代后的反演结果分别如图5a、5b所示。为了定量地对比不同方法的反演结果，抽取了地下100m、200m处的速度曲线，分别如图6a、图6b所示。由于模型左右两端并没有炮点覆盖，无法得到正确结果，因此仅显示了水平方向10km-30km的速度模型。对比不同深度的速度曲线可知，散射积分法的反演结果相比于传统层析成像方法有着较高的分辨率和精度。这可能是因为散射积分法稳定性较好，并且利用了预条件算子。在物理上，预条件算子能够补偿照明，从而反演更深部的信息并且提高反演的精确度。表1为两种方法的内存占用量与计算效率对比。
[0132]
表1
[0133] 内存占用量一次迭代反演所需时间传统层析233mb39.92s本方法19mb37.79s
[0134]
可以看出，散射积分法层析的内存占用量仅为传统层析的8.2％，计算效率也略高，体现了本方法的优势。
[0135]
实施例2
[0136]
图7示出了根据本发明的实施例2的初始表层速度模型的示意图。
[0137]
图8示出了根据本发明的实施例2的最终表层速度模型的示意图。
[0138]
图9示出了根据本发明的实施例2的层析静校正叠加剖面的示意图。
[0139]
本实施例将本发明提出的基于改进的散射积分法的非线性射线走时层析应用到我国西部四川盆地获取的实际资料中。模型共有2957
×
144个网格，网格间距15m
×
15m。共500炮，炮间距60m，第一炮在地表水平位置7185m处。每炮480道，均匀分布在炮点两端，道间距30m。该观测系统最大偏移距7185m，最小偏移距15m。具体实施流程与实施例1相似，其中初始表层速度模型如图7所示，显示器输出的最终表层速度模型如图8所示。层析静校正叠加剖面如图9所示，反射同相轴比较光滑和连续，也进一步验证了本发明可以得到合理的结果。
[0140]
实施例3
[0141]
本公开提供一种电子设备包括，该电子设备包括：存储器，存储有可执行指令；处理器，处理器运行存储器中的可执行指令，以实现上述非线性射线走时层析成像方法。
[0142]
根据本公开实施例的电子设备包括存储器和处理器。
[0143]
该存储器用于存储非暂时性计算机可读指令。具体地，存储器可以包括一个或多个计算机程序产品，该计算机程序产品可以包括各种形式的计算机可读存储介质，例如易失性存储器和/或非易失性存储器。该易失性存储器例如可以包括随机存取存储器(ram)和/或高速缓冲存储器(cache)等。该非易失性存储器例如可以包括只读存储器(rom)、硬盘、闪存等。
[0144]
该处理器可以是中央处理单元(cpu)或者具有数据处理能力和/或指令执行能力的其它形式的处理单元，并且可以控制电子设备中的其它组件以执行期望的功能。在本公开的一个实施例中，该处理器用于运行该存储器中存储的该计算机可读指令。
[0145]
本领域技术人员应能理解，为了解决如何获得良好用户体验效果的技术问题，本
实施例中也可以包括诸如通信总线、接口等公知的结构，这些公知的结构也应包含在本公开的保护范围之内。
[0146]
有关本实施例的详细说明可以参考前述各实施例中的相应说明，在此不再赘述。
[0147]
实施例4
[0148]
本公开实施例提供一种计算机可读存储介质，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现所述的非线性射线走时层析成像方法。
[0149]
根据本公开实施例的计算机可读存储介质，其上存储有非暂时性计算机可读指令。当该非暂时性计算机可读指令由处理器运行时，执行前述的本公开各实施例方法的全部或部分步骤。
[0150]
上述计算机可读存储介质包括但不限于：光存储介质(例如：cd－rom和dvd)、磁光存储介质(例如：mo)、磁存储介质(例如：磁带或移动硬盘)、具有内置的可重写非易失性存储器的媒体(例如：存储卡)和具有内置rom的媒体(例如：rom盒)。
[0151]
本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
[0152]
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

技术特征：
1.一种非线性射线走时层析成像方法，其特征在于，包括：步骤1：对原始地震数据进行去噪和反褶积的预处理，获得预处理地震数据；步骤2：针对所述预处理地震数据进行初至拾取，获得初至波到达时；步骤3：根据地下先验信息建立初始表层速度模型，设定反演参数；步骤4：根据初至波到达时与所述初始表层速度模型通过散射积分法进行非线性射线走时层析成像反演，获得最终表层速度模型并成像。2.根据权利要求1所述的非线性射线走时层析成像方法，其中，所述步骤4包括：步骤41：在所述初始表层速度模型上进行射线追踪，获取理论合成初至走时与射线路径，保存所述射线路径；步骤42：判断所述理论合成初至走时与所述初至波到达时是否匹配，若是，则以所述初始表层速度模型作为所述最终表层速度模型，若否，则执行步骤43；步骤43：逐个读取每一个炮检对对应的射线路径作为该炮检对对应的核函数；步骤44：通过核函数-标量乘获得每一个炮检对所对应的梯度，进而计算整个观测系统的全局梯度，并计算所述全局梯度的预条件算子；步骤45：计算预条件最速下降方向、方向更新步长与模型更新量；步骤46：更新表层速度模型，以更新后的表层速度模型为初始表层速度模型，重复步骤41-46，直到获得所述最终表层速度模型。3.根据权利要求2所述的非线性射线走时层析成像方法，其中，所述步骤41包括：步骤411：在所述初始表层速度模型网格剖分的网格边界上设定多个关键点；步骤412：自激发点开始以矩形为波前，按惠更斯原理逐层向外扩展，直至扫描完成整个模型空间，再由外向里收缩，直至扫描到激发点，如此重复直至模型空间的最小旅行时场不变为止；步骤413：检波点处的初至波到达时即为所述理论合成初至走时。4.根据权利要求3所述的非线性射线走时层析成像方法，其中，所述步骤44包括：步骤441：获取每一个炮检对所对应的射线路径，进而获得所述射线路径与该炮检对对应的走时残差乘，得到该炮检对对应的梯度以及射线路径向量的元素平方向量；步骤442：将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度；步骤443：将所有的炮检对的核函数元素取平方并进行累加，计算全局梯度的预条件算子。5.根据权利要求4所述的非线性射线走时层析成像方法，其中，通过公式(1)将所有炮检对的梯度累加求和，得到整个观测系统的全局梯度：其中，k
t
δt为全局梯度，k
ij
为第i对炮检对对应的核函数在第j个模型网格内的核函数值，δt
i
为旅行时差，k
ij
δt
i
为第i对炮检对的梯度，k
t
为影响走时的核函数矩阵、δt为旅行时差。
6.根据权利要求4所述的非线性射线走时层析成像方法，其中，通过公式(2)计算全局梯度的预条件算子：其中，h0为预条件算子。7.根据权利要求6所述的非线性射线走时层析成像方法，其中，所述步骤45中，通过公式(3)计算预条件最速下降方向：其中，p为预条件最速下降方向。8.根据权利要求2所述的非线性射线走时层析成像方法，其中，所述步骤45中方向更新步长的计算方法为：根据先验信息获取所述初始表层速度模型与真实地下模型的速度差，将速度差除以预先设定的迭代次数，得到每次迭代可更新的最大速度值δv
max
，利用待定系数法求取所述方向更新步长t，使得max{|p
·
t|}＝δv
max
。9.一种电子设备，其特征在于，所述电子设备包括：存储器，存储有可执行指令；处理器，所述处理器运行所述存储器中的所述可执行指令，以实现权利要求1-8中任一项所述的非线性射线走时层析成像方法。10.一种计算机可读存储介质，其特征在于，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现权利要求1-8中任一项所述的非线性射线走时层析成像方法。

技术总结
本申请公开了一种非线性射线走时层析成像方法、电子设备及介质。该方法可以包括：步骤1：对原始地震数据进行去噪和反褶积的预处理，获得预处理地震数据；步骤2：针对预处理地震数据进行初至拾取，获得初至波到达时；步骤3：根据地下先验信息建立初始表层速度模型，设定反演参数；步骤4：根据初至波到达时与初始表层速度模型通过散射积分法进行非线性射线走时层析成像反演，获得最终表层速度模型并成像。本发明避免了传统层析成像中的大规模矩阵的存储和预条件实现不方便的问题，实现带预条件的最速下降法地震走时层析成像。最速下降法地震走时层析成像。最速下降法地震走时层析成像。


技术研发人员：徐蔚亚 朱成宏 谢飞 魏哲枫 魏巍
受保护的技术使用者：中国石油化工股份有限公司石油勘探开发研究院
技术研发日：2022.03.25
技术公布日：2023/10/6