多道的反射系数脉冲确定方法、存储介质及电子设备与流程

1.本发明属于地质勘探技术领域，尤其涉及一种多道的反射系数脉冲确定方法、存储介质及电子设备。


背景技术：

2.地震反射系数脉冲反演在推断地下构造和展现地震剖面中更多细节等方面扮演着重要的角色，也被广泛研究。反射系数脉冲反演的基本假设是地震记录可以看成是反射系数脉冲和地震子波的褶积，反射系数脉冲反演可以看成是高分辨率处理的一种类似于反褶积的方法。同时，由于地震反射系数脉冲与波阻抗密切相关，因此也可以被看成是波阻抗反演的基础。
3.近年来，随着油气藏的勘探目标从构造转向岩性油气藏，高分辨率地震资料的获取越来越受到关注。然而由于地震子波的固有带限特征，导致地震资料缺失高低频成分，传统的反褶积和谱白化方法，由于只能增强或恢复地震频带范围内的频谱信息，获得带限的反射系数脉冲，并不能从根本上拓宽地震资料的频带宽度，获得状的反射系数脉冲。难以准确描述地层边界和一些特殊地质体，如薄储层、透镜状、单砂体以及地层岩性尖灭位置。因此，需要提出新的能够获得宽频带的反射系数脉冲的反演方法。


技术实现要素：

4.本发明正是基于上述技术问题，提出了一种多道的反射系数脉冲确定方法、系统、存储介质及电子设备。
5.第一方面，本发明实施例提供了一种多道的反射系数脉冲确定方法，包括：
6.基于原始地震记录的反射系数脉冲的时间位置以及斜率，利用褶积模型，构建用于反射系数脉冲反演的多道正演模型；
7.基于所述多道正演模型，利用模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的反射系数脉冲的振幅幅值；
8.针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值；
9.基于满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得重构的反射系数脉冲序列。
10.可选地，所述多道正演模型为：
[0011][0012]
其中，s
(l+δl)
为第l+δl道地震道，k为非零反射系数脉冲的个数，w为地震子波，
为第l+δl道地震道的第k个反射系数脉冲的振幅幅值，为第l+δl道地震道的第k个反射系数脉冲的时间位置，l为目标地震道的标号，δl为目标地震道与相邻地震道的间隔，n
(l+δl)
为第l+δl道地震道的噪声，t为时间；
[0013]
其中，其中，floor为向下取整，为第l道地震道的第k个反射系数脉冲的斜率，且
[0014]
可选地，针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，包括：
[0015]
针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数计算获得的反射系数脉冲的时间位置以及振幅幅值对应的系统能量；其中，所述目标函数为：
[0016][0017]
其中，ja为系统能量，λ为阻尼系数，a
(l)
为第l道地震道的反射系数脉冲的振幅幅值向量，a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0018][0019]
其中，a
ik
为矩阵a中第i行第k列的元素；
[0020]
当所述系统能量达到预设阈值时，确定该次迭代过程中获得的反射系数脉冲的时间位置以及振幅幅值为满足预设条件的反射系数脉冲的时间位置和振幅幅值。
[0021]
可选地，所述基于所述多道正演模型，利用模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的反射系数脉冲的振幅幅值，包括：
[0022]
基于所述多道正演模型，分别利用具有不同的初始种子的模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的振幅幅值进行迭代；
[0023]
针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，包括：
[0024]
针对具有不同的初始种子的模拟退火反演算法每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，均利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，以获得每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值；
[0025]
基于满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得重构的反射系数脉冲序列，包括：
[0026]
分别基于每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得每个具有不同的初始种子的模拟退火反演算
法对应的反射系数脉冲序列；
[0027]
对获得的每个反射系数脉冲序列取平均值，得到平均反射系数脉冲序列，并将所述平均反射系数脉冲序列作为重构的反射系数脉冲序列。
[0028]
可选地，所述方法还包括：
[0029]
对所述重构的反射系数脉冲序列进行差异性滤波，获得地震子波；
[0030]
将所述重构的反射系数脉冲序列与所述地震子波进行褶积，得到分辨率比所述原始地震记录的分辨率更高的地震记录。
[0031]
可选地，对所述重构的反射系数脉冲序列进行差异性滤波，获得地震子波，包括：
[0032]
确定所述原始地震记录所在的工区的薄层厚度、纵波速度以及密度数据；
[0033]
根据所述薄层厚度，对所述纵波速度以及所述密度数据进行方波化，得到方波化的纵波速度以及密度数据；
[0034]
基于所述方波化的纵波速度以及密度数据，确定所述工区的测井反射系数脉冲；
[0035]
基于所述测井反射系数脉冲的频谱对所述重构的反射系数脉冲序列的频谱进行标定，并基于标定后的反射系数脉冲序列的频谱提取地震子波。
[0036]
可选地，所述基于所述测井反射系数脉冲的频谱对所述重构的反射系数脉冲序列的频谱进行标定，并基于标定后的反射系数脉冲序列的频谱提取地震子波，包括：
[0037]
利用滑动频率域窗口确定所述测井反射系数脉冲的频谱以及所述重构的反射系数脉冲序列的频谱在选定频段内的相关系数；
[0038]
根据所述相关系数在所述重构的反射系数脉冲序列的频谱的第一频率范围与第二频率范围的突变点上确定梯形滤波算子的起始与终止频率，并对确定初始与终止频率的梯形滤波算子进行归一化处理，获得归一化的梯形滤波算子；
[0039]
以归一化的梯形滤波算子与所述相关系数的乘积作为频率域滤波器的参数，对所述频率域滤波器进行反傅里叶变换，得到时间域的地震子波。
[0040]
第二方面，本发明实施例提供了一种多道的反射系数脉冲确定系统，包括：
[0041]
多道正演模型构建模块，用于基于原始地震记录的反射系数脉冲的时间位置以及斜率，利用褶积模型，构建用于反射系数脉冲反演的多道正演模型；
[0042]
迭代模块，用于基于所述多道正演模型，利用模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的反射系数脉冲的振幅幅值；
[0043]
确定模块，用于针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值；
[0044]
重构模块，用于基于满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得重构的反射系数脉冲序列。
[0045]
可选地，所述多道正演模型为：
[0046][0047]
其中，s
(l+δl)
为第l+δl道地震道，k为非零反射系数脉冲的个数，w为地震子波，
为第l+δl道地震道的第k个反射系数脉冲的振幅幅值，为第l+δl道地震道的第k个反射系数脉冲的时间位置，l为目标地震道的标号，δl为目标地震道与相邻地震道的间隔，n
(l+δl)
为第l+δl道地震道的噪声，t为时间；
[0048]
其中，其中，floor为向下取整，为第l道地震道的第k个反射系数脉冲的斜率，且
[0049]
可选地，所述确定模块具体用于：
[0050]
针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数计算获得的反射系数脉冲的时间位置以及振幅幅值对应的系统能量；其中，所述目标函数为：
[0051][0052]
其中，ja为系统能量，λ为阻尼系数，a
(l)
为第l道地震道的反射系数脉冲的振幅幅值向量，a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0053][0054]
其中，a
ik
为矩阵a中第i行第k列的元素；
[0055]
当所述系统能量达到预设阈值时，确定该次迭代过程中获得的反射系数脉冲的时间位置以及振幅幅值为满足预设条件的反射系数脉冲的时间位置和振幅幅值。
[0056]
可选地，所述迭代模块具体用于：
[0057]
基于所述多道正演模型，分别利用具有不同的初始种子的模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的振幅幅值进行迭代；
[0058]
所述确定模块具体用于：
[0059]
针对具有不同的初始种子的模拟退火反演算法每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，均利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，以获得每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值；
[0060]
所述重构模块具体用于：
[0061]
分别基于每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得每个具有不同的初始种子的模拟退火反演算法对应的反射系数脉冲序列；
[0062]
对获得的每个反射系数脉冲序列取平均值，得到平均反射系数脉冲序列，并将所述平均反射系数脉冲序列作为重构的反射系数脉冲序列。
[0063]
可选地，所述系统还包括：
[0064]
滤波模块，用于对所述重构的反射系数脉冲序列进行差异性滤波，获得地震子波；
[0065]
褶积模块，用于将所述重构的反射系数脉冲序列与所述地震子波进行褶积，得到分辨率比所述原始地震记录的分辨率更高的地震记录。
[0066]
可选地，所述滤波模块包括：
[0067]
确定单元，用于确定所述原始地震记录所在的工区的薄层厚度、纵波速度以及密度数据；
[0068]
方波化单元，用于根据所述薄层厚度，对所述纵波速度以及所述密度数据进行方波化，得到方波化的纵波速度以及密度数据；
[0069]
测井反射系数脉冲单元，用于基于所述方波化的纵波速度以及密度数据，确定所述工区的测井反射系数脉冲；
[0070]
标定单元，用于基于所述测井反射系数脉冲的频谱对所述平均反射系数脉冲的频谱进行标定，并从标定后的平均反射系数脉冲的频谱中提取地震子波。
[0071]
可选地，所述标定单元包括：
[0072]
滑动窗口单元，用于利用滑动频率域窗口确定所述测井反射系数脉冲的频谱以及所述重构的反射系数脉冲序列的频谱在选定频段内的相关系数；
[0073]
归一化单元，用于根据所述相关系数在所述重构的反射系数脉冲序列的频谱的第一频率范围与第二频率范围的突变点上确定梯形滤波算子的起始与终止频率，并对确定初始与终止频率的梯形滤波算子进行归一化处理，获得归一化的梯形滤波算子；
[0074]
提取单元，用于以归一化的梯形滤波算子与所述相关系数的乘积作为频率域滤波器的参数，对所述频率域滤波器进行反傅里叶变换，得到时间域的地震子波。
[0075]
第三方面，本发明实施例提供了一种存储介质，所述存储介质上存储有程序代码，所述程序代码被处理器执行时，实现如上述实施例中任一项所述的多道的反射系数脉冲确定方法。
[0076]
第四方面，本发明实施例提供了一种电子设备，所述电子设备包括存储器、处理器，所述存储器上存储有可在所述处理器上运行的程序代码，所述程序代码被所述处理器执行时，实现如上述实施例中任一项所述的多道的反射系数脉冲确定方法。
[0077]
在本发明实施例提供的一种多道的反射系数脉冲确定方法、系统、存储介质及电子设备，通过引入原始地震记录的反射系数脉冲的时间位置和斜率，构建了多道地震资料多道正演模型，并利用模拟退火反演算法同时搜索目标地震道的反射系数脉冲的时间位置和斜率。通过这样的非线性搜索，能够有效的提升了反演预设参数的选择鲁棒性，并且有效避免反演中野值的影响，提高了反射系数脉冲的反演精度和稳定性。
附图说明
[0078]
通过结合附图阅读下文示例性实施例的详细描述可更好地理解本公开的范围。其中所包括的附图是：
[0079]
图1示出了本发明实施例一提出的一种多道的反射系数脉冲确定方法的流程示意图；
[0080]
图2示出了模拟地震资料的反射系数脉冲反演示例图；
[0081]
图3示出了对加入野值的模拟地震资料进行反演的示意图；
[0082]
图4示出了本发明实施例二提出的一种多道的反射系数脉冲确定方法的流程示意图；
[0083]
图5示出了获得具有高分辨率的地震数据的示意图。
具体实施方式
[0084]
为使本发明的目的、技术方案和优点更加清楚，以下将结合附图及实施例来详细说明本发明的实施方法，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。
[0085]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
[0086]
随机稀疏反射系数脉冲反演在地层稀疏的假设前提下，将反射系数脉冲反演问题看作是两个问题的综合体，获得的状反射系数脉冲可以有效拓展地震资料的频谱。但是随机稀疏反射系数脉冲反演需要解决两个问题，第一个问题是一个非线性问题，即反射系数脉冲位置的确定问题；第二个问题是一个线性问题，即确定位置后的反射系数脉冲的幅值大小求取问题。这种反射系数脉冲反演思路是可以被看成是严格意义上的稀疏反射系数脉冲反演。这类反演方法存在两个问题。
[0087]
第一，这类反演的输入参数，包括阻尼因子和非零反射系数脉冲个数，通常根据经验由人工选定，而不合理的参数选择经常会引发过拟合问题。发现非零反射系数脉冲的个数设定会很大程度的影响反演的最终结果，不仅仅影响非零反射系数脉冲的位置，也会影响非零反射系数脉冲的振幅估计。因此需要谨慎的选取非零反射系数脉冲的个数预设值，并且在反演前进行大量实验。虽然通过在频率域引入有限更新率理论改善了选取预设非零反射系数脉冲的稳定性，但是预设非零反射系数脉冲个数当做是一个未知数，在估计反射系数脉冲的同时估计非零反射系数脉冲个数，然而使得反演维度增加，成本增大。
[0088]
第二，由于反演过程中，算法收敛取决于地震数据和反演反射系数脉冲的正演模拟之间的匹配程度，具有较大反射系数脉冲幅值的反射系数脉冲会对反演的走向影响巨大。因此，以幅值较大为特征的野值对这类反演的危害很大。
[0089]
针对上述技术问题，本公开通过引入原始地震记录的反射系数脉冲的时间位置和斜率，构建了多道地震资料的多道正演模型，并利用模拟退火反演算法同时搜索原始地震记录的反射系数脉冲的时间位置和斜率。通过这样的非线性搜索，能够有效的提升了反演预设参数的选择鲁棒性，并且有效避免反演中野值的影响，提高了反射系数脉冲的反演精度和稳定性。而且保留了反射系数脉冲估计的不确定性分析优势和识别薄层尖灭的能力，获得的高分辨率反射系数脉冲反映出了地震资料中的更多信息，有利于后续的解释工作。
[0090]
实施例一
[0091]
根据本发明的实施例，提供了一种多道的反射系数脉冲确定方法，图1示出了本发明实施例一提出的一种多道的反射系数脉冲确定方法的流程示意图，如图1所示，该多道的反射系数脉冲确定方法可以包括：步骤110至步骤130。
[0092]
在步骤110中，基于原始地震记录的反射系数脉冲的时间位置以及斜率，利用褶积模型，构建用于反射系数脉冲反演的多道正演模型。
[0093]
这里，地震记录可以描述为反射系数脉冲与地震子波的褶积。假设子波在某一时间和深度范围内是稳态的，那么褶积模型可以写成如下形式：
[0094]s(l)
(t)＝w(t)*r
(l)
(t)+n
(l)
(t)，l＝1,...,l，t＝1,...,n
ꢀꢀ
(1)
[0095][0096]
将计算式2代入计算式1中，可以得到：
[0097][0098]
其中，s
(l)
为目标地震记录，n为噪声，r为反射系数脉冲，为振幅，为时间位置，t为时间。
[0099]
则基于反射系数脉冲的时间位置以及斜率，结合褶积模型构建而成的多道正演模型可以表示为：
[0100][0101]
其中，s
(l+δl)
为第l+δl道地震道，k为非零反射系数脉冲的个数，w为地震子波，为第l+δl道地震道的第k个反射系数脉冲的振幅幅值，为第l+δl道地震道的第k个反射系数脉冲的时间位置，l为目标地震道的标号，δl为目标地震道与相邻地震道的间隔，n
(l+δl)
为第l+δl道地震道的噪声，t为时间；
[0102]
其中，其中，floor为向下取整，为第l道地震道的第k个反射系数脉冲的斜率，且
[0103]
值得说明的是，为了加入相邻道的约束，假设即相邻地震道的振幅大小变化小，并引入斜率其中：
[0104][0105]
其中，δt是时间采样间隔。则令
[0106][0107]
其中，floor(x)代表向下取整。
[0108]
用矩阵形式表示为：
[0109][0110]
其中，a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0111][0112]
其中，a
ik
为矩阵a中第i行第k列的元素。
[0113]
对于该多道正演模型，未知的参数有可以通过最小化地震资料与正演计算获得记录的残差进行最优化估计。
[0114]
在步骤120中，基于所述多道正演模型，利用模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，
利用最小二乘法确定所述时间位置以及斜率对应的反射系数脉冲的振幅幅值。
[0115]
这里，模拟退火反演算法源于统计热力物理学，通过模拟熔融状态下物体缓慢冷却达到结晶状态的物理过程。模拟退火反演算法的基本思想是：生成一系列参数向量模拟粒子的热运动，通过缓慢减小一个模拟温度的控制参数，使模拟的系统最终冷却结晶达到系统能量最小值的过程。模拟退火反演算法实质是利用了地球物理反演问题求解过程与熔化固体退火过程的相似性，开辟了地球物理反演的新途径，是非线性反演算法中一种比较常用的算法。
[0116]
值得说明的是，模拟退火反演算法的迭代过程实际上是通过设置一个初始温度、迭代次数进行迭代的。例如，第一次迭代，利用初始温度搜索到一个时间位置以及斜率，并基于初始温度生成一个新的温度；在第二次迭代中，基于新的温度搜索得到另一个时间位置以及斜率，直至满足预设迭代次数或者满足预设收敛条件，才停止迭代。
[0117]
而针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的反射系数脉冲的振幅幅值。
[0118]
值得说明的是，最小二乘法是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。通过利用最小二乘法对搜索到的反射系数脉冲的时间位置以及斜率进行计算，得到对应的振幅幅值。
[0119]
在步骤130中，针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值。
[0120]
这里，模拟退火反演算法每次迭代获得的反射系数脉冲的时间位置以及对应的振幅幅值均利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件，从而确定迭代获得的反射系数脉冲的时间位置以及对应的振幅幅值是否为可用的反射系数脉冲的时间位置以及对应的振幅幅值。
[0121]
在一个可选的实施方式中，步骤130中，针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，包括：
[0122]
针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数计算获得的反射系数脉冲的时间位置以及振幅幅值对应的系统能量；其中，所述目标函数为：
[0123][0124]
其中，ja为系统能量，λ为阻尼系数，a
(l)
为第l道地震道的反射系数脉冲的振幅幅值向量，a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0125][0126]
其中，a
ik
为矩阵a中第i行第k列的元素；
[0127]
当所述系统能量达到预设阈值时，确定该次迭代过程中获得的反射系数脉冲的时间位置以及振幅幅值为满足预设条件的反射系数脉冲的时间位置和振幅幅值。
[0128]
这里，考虑到一般地震资料中的噪音，在构建模拟退火反演算法的目标函数时，需要增加约束项，在本公开中，增加了对反演获得的反射系数脉冲的二范数约束，构建的模拟退火反演算法的目标函数为：
[0129][0130]
a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0131][0132]
其中，a
ik
为矩阵a中第i行第k列的元素。
[0133]
在求解目标函数进而获得反射系数脉冲之前，需要预先设定非零反射系数脉冲的个数，即k的值，一般大于人为预估的反射系数脉冲个数。然后通过模拟退火反演方法搜索原始地震记录的反射系数脉冲的时间位置以及斜率，然后利用最小二乘法获得目标地震道的反射系数脉冲的振幅大小。
[0134]
其中，通过使用模拟退火反演算法搜索时间位置以及斜率的同时，将ja最小化，振幅幅值则通过在每次模拟退火反演算法迭代过程中利用最小二乘法计算得到。
[0135]
值得说明的是，最小二乘法是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。
[0136]
在步骤140中，基于满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得重构的反射系数脉冲序列。
[0137]
这里，将反演得到的振幅幅值以及时间位置代入到反射系数脉冲表达式中，重建反射系数脉冲序列。
[0138]
其中，反射系数脉冲表达式为：
[0139][0140]
其中，r(t)为反射系数脉冲序列，k为非零反射系数脉冲的个数，τk为第k个反射系数脉冲的时间位置，ak为第k个反射系数脉冲的振幅幅值，δ(t)为函数。
[0141]
在本实施例中，通过引入原始地震记录的反射系数脉冲的时间位置和斜率，构建了多道地震资料多道正演模型，并利用模拟退火反演算法同时搜索原始地震记录的反射系数脉冲的时间位置和斜率。通过这样的非线性搜索，能够有效的提升了反演预设参数的选择鲁棒性，并且有效避免反演中野值的影响，提高了反射系数脉冲的反演精度和稳定性。
[0142]
下面，通过一个具体的示例对上述实施例的内容进行说明：
[0143]
图2示出了模拟地震资料的反射系数脉冲反演示例图，图2中(a)为反射系数脉冲模型，(b)为含噪声地震记录，(c)为结合多道地震资料的多道正演模型反演结果，(d)为单道反射系数脉冲反演结果。
[0144]
如图2所示，多道的反演结果图2(c)较单道的反演结果图2(d)能够识别较薄的界面信息，并且对噪声的抵抗能力也更强一些，图2(c)中由噪声引起的毛刺装的反射系数脉冲明显减少。
[0145]
图3示出了对加入野值的模拟地震资料进行反演的示意图，图3中(b)为加入野值的模拟地震资料，图3中(c)为利用本公开提出的方法进行反演的反演结果，图3中(a)为真实的反射系数脉冲。
[0146]
如图3所示，图3(c)展示了本公开提出的反演方法对野值的压制效果，由于很多稀疏反演方法力求用最小的反射系数脉冲个数达到更好的资料匹配，这种野值类噪声会对反演结果造成很大的影响，本公开提出的反演方法显然会降低这样的影响，获得了较为准确的反演结果。
[0147]
实施例二
[0148]
在上述实施例的基础上，本发明的实施例二还可以提供一种多道的反射系数脉冲确定方法。图4示出了本发明实施例二提出的一种多道的反射系数脉冲确定方法的流程示意图，如图4所示，该多道的反射系数脉冲确定方法可以包括：步骤210至步骤240。
[0149]
在步骤210中，基于原始地震记录的反射系数脉冲的时间位置以及斜率，利用褶积模型，构建用于反射系数脉冲反演的多道正演模型。
[0150]
这里，叠后地震资料可以描述为反射系数脉冲与地震子波的褶积。假设子波在某一时间和深度范围内是稳态的，那么褶积模型可以写成如下形式：
[0151]s(l)
(t)＝w(t)*r
(l)
(t)+n
(l)
(t)，l＝1,...,l，t＝1,...,n
ꢀꢀ
(1)
[0152][0153]
将计算式2代入计算式1中，可以得到：
[0154][0155]
其中，s
(l)
为目标地震记录，n为噪声，r为反射系数脉冲，为振幅，为时间位置，t为时间。
[0156]
则基于反射系数脉冲的时间位置以及斜率，结合褶积模型构建而成的多道正演模型可以表示为：
[0157][0158]
其中，s
(l+δl)
为第l+δl道地震道，k为非零反射系数脉冲的个数，w为地震子波，为第l+δl道地震道的第k个反射系数脉冲的振幅幅值，为第l+δl道地震道的第l个反射系数脉冲的时间位置，l为目标地震道的标号，δl为目标地震道与相邻地震道的间隔，n
(l+δl)
为第l+δl道地震道的噪声，t为时间；
[0159]
其中，其中，floor为向下取整，为第l道地震道的第k个反射系数脉冲的斜率，且
[0160]
值得说明的是，为了加入相邻道的约束，假设即相邻地震道的振幅大小变化小，并引入斜率其中：
[0161][0162]
其中，δt是时间采样间隔。则令
[0163][0164]
其中，floor(x)代表向下取整。
[0165]
用矩阵形式表示为：
[0166][0167]
其中，a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0168][0169]
其中，a
ik
为矩阵a中第i行第k列的元素。
[0170]
对于该多道正演模型，未知的参数有可以通过最小化地震资料与正演计算获得记录的残差进行最优化估计。
[0171]
在步骤220中，基于所述多道正演模型，分别利用具有不同的初始种子的模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的振幅幅值进行迭代。
[0172]
这里，模拟退火反演算法源于统计热力物理学，通过模拟熔融状态下物体缓慢冷却达到结晶状态的物理过程。模拟退火反演算法的基本思想是：生成一系列参数向量模拟粒子的热运动，通过缓慢减小一个模拟温度的控制参数，使模拟的系统最终冷却结晶达到系统能量最小值的过程。模拟退火反演算法实质是利用了地球物理反演问题求解过程与熔化固体退火过程的相似性，开辟了地球物理反演的新途径，是非线性反演算法中一种比较常用的算法。
[0173]
在步骤230中，针对具有不同的初始种子的模拟退火反演算法每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，均利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，以获得每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值。
[0174]
这里，关于利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值的计算过程在上述实施例中已详细说明，在此不再赘述。
[0175]
值得说明的是，在步骤230中实际上是通过调整模拟退火反演算法初始化时的初始种子，从而利用不同的初始种子进行多次反演，从而获得反射系数脉冲的多次的独立的估计结果，即多次实现反射系数脉冲反演结果。
[0176]
在步骤240中，分别基于每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得每个具有不同的初始种子的模拟退火反演算法对应的反射系数脉冲序列。
[0177]
在步骤250中，对获得的每个反射系数脉冲序列取平均值，得到平均反射系数脉冲
序列，并将所述平均反射系数脉冲序列作为重构的反射系数脉冲序列。
[0178]
这里，将反演得到的振幅幅值以及时间位置代入到反射系数脉冲表达式中，构建反射系数脉冲。
[0179]
其中，反射系数脉冲表达式为：
[0180][0181]
其中，r(t)为反射系数脉冲序列，k为非零反射系数脉冲的个数，τk为第k个反射系数脉冲的时间位置，ak为第k个反射系数脉冲的振幅幅值，δ(t)为函数。
[0182]
值得说明的是，多次反演得到的振幅幅值、时间位置以及斜率，可以构建出多个对应的反射系数脉冲。
[0183]
在本公开中，采用50次左右的独立实验后，将获得的反射系数脉冲的平均结果作为最终反射系数脉冲反演结果比较准确合理，这里平均反射系数脉冲结果即为获得的高分辨率反射系数脉冲反演结果。
[0184]
在一个可选的实施方式中，在步骤250之后，所述方法还包括：步骤260至步骤270。
[0185]
在步骤260，对所述重构的反射系数脉冲序列进行差异性滤波，获得地震子波。
[0186]
这里，由于平均反射系数脉冲的全频带信息并非准确，而只有低频、中频和中高频才是相对准确的。并且，获得的平均反射系数脉冲在低频、中频和中高频段的信噪比不同，因此需要对获得的平均反射系数脉冲进行差异性滤波，从而获得拓频后的地震子波。
[0187]
步骤270，将所述重构的反射系数脉冲序列与所述地震子波进行褶积，得到分辨率比所述原始地震记录的分辨率更高的地震记录。
[0188]
这里，通过利用多道的反射系数脉冲确定方法计算得到的平均反射系数脉冲以及滤波获得的宽频带的地震子波进行褶积，即可以获得具有高分辨率的地震数据。
[0189]
在一个可选的实施方式中，步骤260中，对所述重构的反射系数脉冲序列进行差异性滤波，获得地震子波，可以包括：
[0190]
步骤261，确定所述原始地震记录所在的工区的薄层厚度、纵波速度以及密度数据；
[0191]
步骤262，根据所述薄层厚度，对所述纵波速度以及所述密度数据进行方波化，得到方波化的纵波速度以及密度数据；
[0192]
步骤263，基于所述方波化的纵波速度以及密度数据，确定所述工区的测井反射系数脉冲；
[0193]
步骤264，基于所述测井反射系数脉冲的频谱对所述重构的反射系数脉冲序列的频谱进行标定，并基于标定后的反射系数脉冲序列的频谱提取地震子波。
[0194]
这里，可以根据研究区域地质资料和情况，确定薄层厚度，并选取研究区域内测井数据中纵波速度和密度数据。然后根据薄层厚度进行方波化，利用方波化后的纵波速度和密度数据计算得到测井反射系数脉冲，进而利用测井反射系数脉冲的频谱标定平均反射系数脉冲的频谱，从而基于标定后的平均反射系数脉冲的频谱中提取地震子波。
[0195]
在一个可选的实施方式中，步骤264中，基于所述测井反射系数脉冲的频谱对所述重构的反射系数脉冲序列的频谱进行标定，并基于标定后的反射系数脉冲序列的频谱提取
地震子波，包括：
[0196]
步骤2641，利用滑动频率域窗口确定所述测井反射系数脉冲的频谱以及所述重构的反射系数脉冲序列的频谱在选定频段内的相关系数；
[0197]
步骤2642，根据所述相关系数在所述重构的反射系数脉冲序列的频谱的第一频率范围与第二频率范围的突变点上确定梯形滤波算子的起始与终止频率，并对确定初始与终止频率的梯形滤波算子进行归一化处理，获得归一化的梯形滤波算子；
[0198]
步骤2643，以归一化的梯形滤波算子与所述相关系数的乘积作为频率域滤波器的参数，对所述频率域滤波器进行反傅里叶变换，得到时间域的地震子波。
[0199]
这里，归一化是一种简化计算的方式，即将有量纲的表达式，经过变换，化为无量纲的表达式，成为标量。即把数据变成(0,1)或者(1,1)之间的小数。主要是为了数据处理方便提出来的，把数据映射到0～1范围之内处理，更加便捷快速。把有量纲表达式变成无量纲表达式，便于不同单位或量级的指标能够进行比较和加权。
[0200]
频率域滤波是指先修改一副图像的傅里叶变换，然后计算其傅里叶反变换，从而得到处理结果，在此不再展开描述。
[0201]
值得说明的是，第一频率范围与第二频率范围的突变点指的是频谱中低频与高频的突变点。
[0202]
在本实施例中，通过引入原始地震记录的反射系数脉冲的时间位置和斜率，构建了多道地震资料多道正演模型，并利用模拟退火反演算法同时搜索原始地震记录的反射系数脉冲的时间位置和斜率。通过这样的非线性搜索，能够有效的提升了反演预设参数的选择鲁棒性，并且有效避免反演中野值的影响，提高了反射系数脉冲的反演精度和稳定性。以及通过多次不同的初始种子的退火模拟反演方法获得多个反射系数脉冲，并基于此求取平均反射系数脉冲，从而能够获得准确的反射系数脉冲。并且基于获得的平均反射系数脉冲提取宽频带的地震子波，从而将宽频带的地震子波与平均反射系数脉冲进行褶积，获得具有高分辨率的地震数据，从而实现提高地震资料的分辨率的目的。
[0203]
下面，通过一个具体的示例对上述实施方式进行说明：
[0204]
图5示出了获得具有高分辨率的地震数据的示意图，其中，图5中(a)图为地震资料，图5中(b)图为多次实现反射系数脉冲反演结果的示意图，图5中(c)图为利用平均反射系数脉冲获得的高分辨率地震剖面的示意图。
[0205]
如图5所示，该地区地下构造较为复杂，界面横向变化快，需要对地下构造和起伏进行精确的描述，利用的子波是通过井震标定获得的。地震记录和多次反演(50次)实现的平均结果，如图5(a)和图5(b)所示，反射系数脉冲反演后的地震剖面上更加清晰的展现了更多的地质构造和界面信息。图5(c)给出了利用宽频带的地震子波与平均反射系数脉冲获得的高分辨率地震记录，其中断层和构造的变化展现得更加清晰，该结果可以有效的改善地震解释的精度。证实了本公开提出的方法在实际地震资料中的应用潜力，对后续解释和储层的预测具有极高的指导意义。
[0206]
实施例三
[0207]
根据本发明的实施例，还提供了一种多道的反射系数脉冲确定系统，包括：
[0208]
多道正演模型构建模块，用于基于原始地震记录的反射系数脉冲的时间位置以及斜率，利用褶积模型，构建用于反射系数脉冲反演的多道正演模型；
[0209]
迭代模块，用于基于所述多道正演模型，利用模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的反射系数脉冲的振幅幅值；
[0210]
确定模块，用于针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数确定该反射系数脉冲的时间位置以及对应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值；
[0211]
重构模块，用于基于满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得重构的反射系数脉冲序列。
[0212]
可选地，所述多道正演模型为：
[0213][0214]
其中，s
(l+δl)
为第l+δl道地震道，k为非零反射系数脉冲的个数，w为地震子波，为第l+δl道地震道的第k个反射系数脉冲的振幅幅值，为第l+δl道地震道的第k个反射系数脉冲的时间位置，l为目标地震道的标号，δl为目标地震道与相邻地震道的间隔，n
(l+δl)
为第l+δl道地震道的噪声，t为时间；
[0215]
其中，其中，floor为向下取整，为第l道地震道的第k个反射系数脉冲的斜率，且
[0216]
可选地，所述确定模块具体用于：
[0217]
针对每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，利用预设的目标函数计算获得的反射系数脉冲的时间位置以及振幅幅值对应的系统能量；其中，所述目标函数为：
[0218][0219]
其中，ja为系统能量，λ为阻尼系数，a
(l)
为第l道地震道的反射系数脉冲的振幅幅值向量，a为n
×
k大小的矩阵，n为采样序列长度，该矩阵中的元素表示为：
[0220][0221]
其中，a
ik
为矩阵a中第i行第k列的元素；
[0222]
当所述系统能量达到预设阈值时，确定该次迭代过程中获得的反射系数脉冲的时间位置以及振幅幅值为满足预设条件的反射系数脉冲的时间位置和振幅幅值。
[0223]
可选地，所述迭代模块具体用于：
[0224]
基于所述多道正演模型，分别利用具有不同的初始种子的模拟退火反演算法进行迭代，搜索反射系数脉冲的时间位置以及斜率，针对每次搜索到的反射系数脉冲的时间位置以及斜率，利用最小二乘法确定所述时间位置以及斜率对应的振幅幅值进行迭代；
[0225]
所述确定模块具体用于：
[0226]
针对具有不同的初始种子的模拟退火反演算法每次获得的反射系数脉冲的时间位置以及对应的振幅幅值，均利用预设的目标函数确定该反射系数脉冲的时间位置以及对
应的振幅幅值是否为满足预设条件的反射系数脉冲的时间位置以及振幅幅值，以获得每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值；
[0227]
所述重构模块具体用于：
[0228]
分别基于每个具有不同的初始种子的模拟退火反演算法对应的满足预设条件的反射系数脉冲的时间位置以及振幅幅值，获得每个具有不同的初始种子的模拟退火反演算法对应的反射系数脉冲序列；
[0229]
对获得的每个反射系数脉冲序列取平均值，得到平均反射系数脉冲序列，并将所述平均反射系数脉冲序列作为重构的反射系数脉冲序列。
[0230]
可选地，所述系统还包括：
[0231]
滤波模块，用于对所述重构的反射系数脉冲序列进行差异性滤波，获得地震子波；
[0232]
褶积模块，用于将所述重构的反射系数脉冲序列与所述地震子波进行褶积，得到分辨率比所述原始地震记录的分辨率更高的地震记录。
[0233]
可选地，所述滤波模块包括：
[0234]
确定单元，用于确定所述原始地震记录所在的工区的薄层厚度、纵波速度以及密度数据；
[0235]
方波化单元，用于根据所述薄层厚度，对所述纵波速度以及所述密度数据进行方波化，得到方波化的纵波速度以及密度数据；
[0236]
测井反射系数脉冲单元，用于基于所述方波化的纵波速度以及密度数据，确定所述工区的测井反射系数脉冲；
[0237]
标定单元，用于基于所述测井反射系数脉冲的频谱对所述平均反射系数脉冲的频谱进行标定，并从标定后的平均反射系数脉冲的频谱中提取地震子波。
[0238]
可选地，所述标定单元包括：
[0239]
滑动窗口单元，用于利用滑动频率域窗口确定所述测井反射系数脉冲的频谱以及所述重构的反射系数脉冲序列的频谱在选定频段内的相关系数；
[0240]
归一化单元，用于根据所述相关系数在所述重构的反射系数脉冲序列的频谱的第一频率范围与第二频率范围的突变点上确定梯形滤波算子的起始与终止频率，并对确定初始与终止频率的梯形滤波算子进行归一化处理，获得归一化的梯形滤波算子；
[0241]
提取单元，用于以归一化的梯形滤波算子与所述相关系数的乘积作为频率域滤波器的参数，对所述频率域滤波器进行反傅里叶变换，得到时间域的地震子波。
[0242]
上述各功能模块的具体实施过程可参见上述实施例，在本实施例中不再重复赘述。
[0243]
实施例四
[0244]
根据本发明的实施例，还提供了一种存储介质，所述存储介质上存储有程序代码，所述程序代码被处理器执行时，实现如上述实施例任一项所述的多道的反射系数脉冲确定方法。
[0245]
实施例五
[0246]
根据本发明的实施例，还提供了一种电子设备，所述电子设备包括存储器、处理器，所述存储器上存储有可在所述处理器上运行的程序代码，所述程序代码被所述处理器
执行时，实现如上述实施例任一项所述的多道的反射系数脉冲确定方法。
[0247]
以上结合附图详细说明了本发明的技术方案，考虑到相关技术中，现有的反演方法如果出现不合理的参数选择经常会引发过拟合问题，不仅仅影响非零反射系数脉冲的位置，也会影响非零反射系数脉冲的振幅估计。本发明提供一种多道的反射系数脉冲确定方法、系统、存储介质及电子设备，通过引入原始地震记录的反射系数脉冲的时间位置和斜率，构建了多道地震资料多道正演模型，并利用模拟退火反演算法同时搜索原始地震记录的反射系数脉冲的时间位置和斜率。通过这样的非线性搜索，能够有效的提升了反演预设参数的选择鲁棒性，并且有效避免反演中野值的影响，提高了反射系数脉冲的反演精度和稳定性。
[0248]
在本技术所提供的几个实施例中，应该理解到，所揭露的装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。
[0249]
作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本发明实施例方案的目的。
[0250]
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以是两个或以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元形式实现。
[0251]
集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台电子设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read-only memory)、随机存取存储器(ram，random access memory)、磁碟或者光盘等各种可以存储程序代码的介质。
[0252]
虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的保护范围，仍须以所附的权利要求书所界定的范围为准。