1.本发明涉及地球物理勘探技术领域，特别涉及一种基于点扩散函数的地震数据高分辨处理方法。
背景技术：2.高分辨和高保真地震反演成像，主要可以分为两大路线。第一条路线是基于bayes估计理论框架的线性反演。即所谓的最小二乘偏移成像，其为全波形反演的一个线性子问题。该方法的优点为对观测系统和数据规则程度的要求比较低，迭代求解过程中的偏移算子可以灵活地选择射线类或是波动类；该方法的缺点是计算效率较低，且对于实际复杂资料存在收敛性问题。保真成像的另一条发展路线为高频近似下的真振幅成像方法。
3.基于bayes估计理论框架的最小二乘偏移成像可以在数据域或成像域实现。一般地，数据域最小二乘偏移成像通常需要10次或者更多的迭代次数，其计算代价要高于常规的偏移成像一个数量级以上。另外，考虑初始模型的不准确、born近似正问题不能很好地模拟波场扰动、噪音不符合高斯分布、观测系统不规则、子波空变且未知等因素，先验信息的引入十分重要，如果没有有效的数据预条件或先验信息的引入，数据域最小二乘偏移成像的迭代收敛可能会很慢，甚至不收敛。对数据域最小二乘偏移成像的误差泛函在梯度等于零时的法方程进行求解，就是在成像域实现最小二乘偏移成像，其关键问题在于显式计算hessian矩阵。本质上，hessian矩阵的某一列对应着模型空间对应一点的点扩散函数，其取决于传播子波，观测系统以及背景速度，是一个空间有限展布的带限波包，具有局部性以及空变性。相较于数据域迭代求解的最小二乘偏移成像，成像域最小二乘偏移成像避免了偏移、反偏移以及步长的计算，总体计算量小、操作灵活。因此，成像域最小二乘偏移成像在实际应用中应该是更合理的选择。
4.然而，受震源子波未知、背景速度准确性和正演/偏移算子精度的影响，点扩散的计算结果不可能是精确的，为此需要构建一套修正不精确点扩散函数的成像域最小二乘偏移成像方法。
技术实现要素：5.本部分的目的在于概述本发明的实施例的一些方面以及简要介绍一些较佳实施例。在本部分以及本技术的说明书摘要和发明名称中可能会做些简化或省略以避免使本部分、说明书摘要和发明名称的目的模糊，而这种简化或省略不能用于限制本发明的范围。
6.鉴于上述存在的问题，提出了本发明。
7.因此，本发明解决的技术问题是：地震数据成像域最小二乘偏移成像中，点扩散函数的不精确所引起的成像假象问题。
8.为解决上述技术问题，本发明提供如下技术方案：一种基于点扩散函数的地震数据高分辨处理方法，包括：
9.输入地震记录、常规偏移结果和偏移速度，并基于所述地震记录的道头信息获取
观测系统；
10.基于震源子波、空间位置的旅行时和振幅解析计算初始点扩散函数；
11.基于所述点扩散函数对所述偏移结果进行反褶积处理，获取反褶积后的成像结果；
12.基于所述反褶积后的成像结果和所述常规偏移结果，更新点扩散函数，最终输出像域最小二乘偏移成像结果和更新后的点扩散函数；
13.作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述输入地震记录，包括
14.输入segy格式和su格式的地震记录。
15.作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述获取观测系统，包括：
16.获取观测系统的参数即炮点空间位置和检波点空间位置。
17.作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述震源子波，包括：
18.假设具有一般性意义的雷克子波作为震源子波参数代入至点扩散函数的解析公式中，用于计算初始点扩散函数。
19.作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述点扩散函数的解析公式，表示为：
[0020][0021]
其中，为空间点xi处的点扩散函数，r
ss
代表震源子波s(ω)的自相关； t表示旅行时，a表示振幅值，散射射线路径的振幅通过连接两条入射射线给出：a(xr；xj；xs)＝a(xj,xr)a(xj,xs)。
[0022]
作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述点扩散函数的解析公式，包括：
[0023]
对于hessian矩阵的计算量问题，引入wkbj近似下的渐进格林函数，得到hessian矩阵的近似表达式即点扩散函数的解析公式，表示为射线振幅和震源子波自相关的积分，这样只需要旅行时、振幅以及震源子波，就可以计算地下任一点的点扩散函数。
[0024]
作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述反褶积处理，包括：
[0025]
通过考虑反射系数的稀疏性和结构连续性，引入l1约束和全变分约束来构建求解反褶积处理对应的误差泛函jr(r)，进而实现反褶积处理；
[0026]
反褶积处理对应的误差泛函jr(r)表示为：
[0027][0028]
其中，*表示褶积算子，f(i)为第i次迭代的点扩散函数，r表示当前迭代中的反褶积，λ1、λ2和λ3表示超参，‖r‖
p
、‖r‖
tv
和δr分别表示对反褶积结果引入的lp约束、全变分约束和二阶导约束，p∈[0,2]。
[0029]
作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，
其中：所述常规偏移成像结果，包括：
[0030]
将常规偏移成像结果看作是点扩散函数与反射系数的褶积结果，对于一维情况，将褶积运算写为toeplitz矩阵乘积形式，对于二维或三维情况，将褶积运算与相关运算视作一对共轭算子，从而进行梯度更新；其中，相关运算具体为计算两个变量之间的相关系数。
[0031]
作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述更新点扩散函数，包括：
[0032]
引入点扩散函数的光滑约束μ1‖f‖q，构建误差泛函jf(f)来更新点扩散函数；
[0033]
此时的误差泛函jf(f)表示为：
[0034][0035]
其中，表示相关算子，r(i)为第i次迭代的反褶积结果，‖f‖q为对点扩散函数引入的lq约束，μ1为其对应的超参，q∈(1,2]。
[0036]
作为本发明所述的基于点扩散函数的地震数据高分辨处理方法的一种优选方案，其中：所述输出像域最小二乘偏移成像结果和更新后的点扩散函数，包括：
[0037]
判断点扩散函数与反褶积结果进行褶积与常规偏移成像结果误差是否小于给定的阈值；当点扩散函数与反褶积结果进行褶积与常规偏移成像结果误差小于给定阈值时输出像域最小二乘偏移成像结果和更新后的点扩散函数。
[0038]
本发明的有益效果：本发明从概率统计反演出发，将强非线性的参数估计反问题转化为一个凸性较好的像域最小二乘反演成像问题；对于hessian矩阵的计算量问题，引入wkbj近似下的渐进格林函数，得到hessian矩阵的近似表达式，实现对地下任一点的点扩散函数的计算；此外，本发明引入点扩散函数和反射系数的先验认识，构建点扩散函数和反射系数迭代更新的策略，得到更高分辨的成像结果。
附图说明
[0039]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。其中：
[0040]
图1为本发明一个实施例提供的一种基于点扩散函数的地震数据高分辨处理方法的整体流程图；
[0041]
图2为本发明一个实施例提供的反射系数剖面；
[0042]
图3为本发明一个实施例提供的常规偏移成像结果；
[0043]
图4为本发明一个实施例提供的初始点扩散函数的反褶积结果；
[0044]
图5为本发明一个实施例提供的高分辨成像结果；
[0045]
图6为本发明一个实施例提供的初始点扩散函数、真实点扩散函数和最后更新的点扩散函数测抽道对比图；
[0046]
图7为本发明一个实施例提供的真实点扩散函数的二维示意图；
[0047]
图8为本发明一个实施例提供的最后更新的点扩散函数的二维示意图；
[0048]
图9为本发明一个实施例提供的初始点扩散函数的二维示意图。
具体实施方式
[0049]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合说明书附图对本发明的具体实施方式做详细的说明，显然所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明的保护的范围。
[0050]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。
[0051]
其次，此处所称的“一个实施例”或“实施例”是指可包含于本发明至少一个实现方式中的特定特征、结构或特性。在本说明书中不同地方出现的“在一个实施例中”并非均指同一个实施例，也不是单独的或选择性的与其他实施例互相排斥的实施例。
[0052]
本发明结合示意图进行详细描述，在详述本发明实施例时，为便于说明，表示器件结构的剖面图会不依一般比例作局部放大，而且所述示意图只是示例，其在此不应限制本发明保护的范围。此外，在实际制作中应包含长度、宽度及深度的三维空间尺寸。
[0053]
同时在本发明的描述中，需要说明的是，术语中的“上、下、内和外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一、第二或第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
[0054]
本发明中除非另有明确的规定和限定，术语“安装、相连、连接”应做广义理解，例如：可以是固定连接、可拆卸连接或一体式连接；同样可以是机械连接、电连接或直接连接，也可以通过中间媒介间接相连，也可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0055]
实施例1
[0056]
参照图1，为本发明的一个实施例，提供了一种基于点扩散函数的地震数据高分辨处理方法，包括：
[0057]
s1：输入地震记录、常规偏移结果和偏移速度，并基于所述地震记录的道头信息获取观测系统；
[0058]
进一步的，输入segy格式和su格式的地震记录；
[0059]
应知道的是，segy格式和su格式的地震数据，每一道数据会有一个道头信息，记录了检波器和震源的相关信息。
[0060]
更进一步的，通过扫描地震数据的道头信息来获取观测系统参数；
[0061]
应知道的是，观测系统参数具体为炮点空间位置和检波点空间位置。
[0062]
s2：基于震源子波、空间位置的旅行时和振幅解析计算初始点扩散函数；
[0063]
进一步的，假设具有一般性意义的雷克子波作为震源子波参数代入至点扩散函数的解析公式中，用于计算初始点扩散函数；
[0064]
应说明的是，如果工区震源已知，也可以输入已知震源子波。
[0065]
更进一步的，对于hessian矩阵的计算量问题，引入wkbj近似下的渐进格林函数，得到hessian矩阵的近似表达式即点扩散函数的解析公式，表示为射线振幅和震源子波自相关的积分；
[0066]
点扩散函数解析公式表示为：
[0067][0068]
其中，为空间点xi处的点扩散函数，r
ss
代表震源子波s(ω)的自相关；t表示旅行时，a表示振幅值，散射射线路径的振幅通过连接两条入射射线给出：a(xr；xj；xs)＝a(xj,xr)a(xj,xs)。
[0069]
应说明的是，引入wkbj近似下的渐进格林函数，得到hessian矩阵的近似表达式后，只需要计算射线的传播时间和振幅以及震源子波即可计算地下任一点的点扩散函数；其中射线的传播时间和振幅可以通过射线追踪得到。
[0070]
s3：基于所述点扩散函数对所述偏移结果进行反褶积处理，获取反褶积后的成像结果；
[0071]
更进一步的，通过考虑反射系数的稀疏性和结构连续性，引入l1约束和全变分约束来构建求解反褶积处理对应的误差泛函jr(r)，进而实现反褶积处理；
[0072]
反褶积处理对应的误差泛函jr(r)表示为：
[0073][0074]
其中，*表示褶积算子，f(i)为第i次迭代的点扩散函数，r表示当前迭代中的反褶积，λ1、λ2和λ3表示超参，‖r‖
p
、‖r‖
tv
和δr分别表示对反褶积结果引入的lp约束、全变分约束和二阶导约束，p∈[0,2]。
[0075]
应知道的是，反褶积处理，就是把一个模糊的像，去除模糊，得到一个更清晰的像，反褶积后的成像结果，指的就是反褶积处理后的结果。
[0076]
应说明的是，反射系数的稀疏性，对应概率统计中长尾巴分布；地地质构造的连续性和光滑性，分别对应全变分正则化和二阶梯度正则化；基于以上两种先验认知并考虑实现后验概率密度最大化，引入l1约束和全变分约束来构建误差泛函jr(r)，进而实现反褶积处理。
[0077]
更进一步的，判断点扩散函数与反褶积结果进行褶积与常规偏移成像结果误差是否小于给定的阈值；若小于给定阈值则达到收敛状态，输出像域最小二乘偏移成像结果和更新后的点扩散函数，否则将进入更新点扩散函数步骤。
[0078]
s4：基于所述反褶积后的成像结果和所述常规偏移结果，更新点扩散函数，最终输出像域最小二乘偏移成像结果和更新后的点扩散函数。
[0079]
进一步的，将常规偏移成像结果看作是点扩散函数与反射系数的褶积结果，对于一维情况，将褶积运算写为toeplitz矩阵乘积形式，对于二维或三维情况，将褶积运算与相关运算视作一对共轭算子，从而进行梯度更新；
[0080]
应说明的是，相关运算是计算两个变量之间的相关系数，从而衡量相似程度。
[0081]
更进一步的，引入点扩散函数的光滑约束μ1‖f‖q，构建误差泛函jf(f)来更新点扩散函数；
[0082]
误差泛函jf(f)表示为：
[0083][0084]
其中，表示相关算子，r(i)为第i次迭代的反褶积结果，‖f‖q为对点扩散函数引入的lq约束，μ1为其对应的超参，q∈(1,2]。
[0085]
应说明的是，点扩散函数的光滑性，对应概率统计中类高斯分布；基于对点扩散函数的先验认知以及考虑后验概率密度最大化，从而引入点扩散函数的光滑约束μ1‖f‖q，来构建误差泛函jf(f)，进而实现对点扩散函数的更新。
[0086]
更进一步的，判断点扩散函数与反褶积结果进行褶积与常规偏移成像结果误差是否小于给定的阈值；若小于给定阈值则达到收敛状态，输出像域最小二乘偏移成像结果和更新后的点扩散函数，否则将进行反褶积处理步骤，直至收敛完成。
[0087]
应说明的是，迭代次数一版小于20次即可收敛。
[0088]
实施例2
[0089]
参照图2～9，为本发明的一个实施例，提供了一种基于点扩散函数的地震数据高分辨处理方法，为了验证本发明的有益效果，通过经济效益计算和仿真实验进行科学论证。
[0090]
首先，考虑溶洞、薄互层和倾斜地层构造，建立速度模型，计算对应的反射系数剖面，如图2所示，模型长乘宽为101乘101个样点，采样间隔5米；
[0091]
其次，使用点扩散函数进行二维褶积得到的结果作为常规偏移成像结果，如图3所示；
[0092]
接下来，采用如图9所示的二维高斯函数作为初始点扩散函数，基于初始点扩散函数和常规偏移结果进行反褶积处理，获得反褶积的成像结果，如图4所示；并基于反褶积的成像结果和常规偏移成像结果，进行点扩散函数的更新；
[0093]
最后，重复10次上述步骤后达到收敛状态，获得高分辨成像结果，如图 5所示。
[0094]
通过图6的初始点扩散函数、真实点扩散函数和最后更新的点扩散函数测抽道对比图以及图7～9的示意图可知本发明所述方法能够修正由于震源子波未知、背景速度和正演/偏移算子不精确引起的点扩散的计算结果的误差，避免错误点扩散函数所引起的错误构造成像，提高成像结果的分辨率。
[0095]
应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。
技术特征：1.一种基于点扩散函数的地震数据高分辨处理方法，其特征在于，包括：输入地震记录、常规偏移结果和偏移速度，并基于所述地震记录的道头信息获取观测系统；基于震源子波、空间位置的旅行时和振幅解析计算初始点扩散函数；基于所述点扩散函数对所述偏移结果进行反褶积处理，获取反褶积后的成像结果；基于所述反褶积后的成像结果和所述常规偏移结果，更新点扩散函数，最终输出像域最小二乘偏移成像结果和更新后的点扩散函数。2.如权利要求1所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述输入地震记录，包括输入segy格式和su格式的地震记录。3.如权利要求1或2所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述获取观测系统，包括：获取观测系统的参数即炮点空间位置和检波点空间位置。4.如权利要求3所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述震源子波，包括：假设具有一般性意义的雷克子波作为震源子波参数代入至点扩散函数的解析公式中，用于计算初始点扩散函数。5.如权利要求4所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述点扩散函数的解析公式，表示为：其中，为空间点x
i
处的点扩散函数，r
ss
代表震源子波s(ω)的自相关；t表示旅行时，a表示振幅值，散射射线路径的振幅通过连接两条入射射线给出：a(x
r
；x
j
；x
s
)＝a(x
j
，x
r
)a(x
j
，x
s
)。6.如权利要求5所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述点扩散函数的解析公式，包括：对于hessian矩阵的计算量问题，引入wkbj近似下的渐进格林函数，得到hessian矩阵的近似表达式即点扩散函数的解析公式，表示为射线振幅和震源子波自相关的积分，这样只需要旅行时、振幅以及震源子波，就可以计算地下任一点的点扩散函数。7.如权利要求6所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述反褶积处理，包括：通过考虑反射系数的稀疏性和结构连续性，引入l1约束和全变分约束来构建求解反褶积处理对应的误差泛函j
r
(r)，进而实现反褶积处理；反褶积处理对应的误差泛函j
r
(r)表示为：其中，*表示褶积算子，f
(i)
为第i次迭代的点扩散函数，r表示当前迭代中的反褶积，λ1、λ2和λ3表示超参，||r||
p
、||r||
tv
和δr分别表示对反褶积结果引入的lp约束、全变分约束和二阶导约束，p∈[0，2]。
8.如权利要求6或7所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述常规偏移成像结果，包括：将常规偏移成像结果看作是点扩散函数与反射系数的褶积结果，对于一维情况，将褶积运算写为toeplitz矩阵乘积形式，对于二维或三维情况，将褶积运算与相关运算视作一对共轭算子，从而进行梯度更新；其中，相关运算具体为计算两个变量之间的相关系数。9.如权利要求8所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述更新点扩散函数，包括：引入点扩散函数的光滑约束μ1||f||
q
，构建误差泛函j
f
(f)来更新点扩散函数；此时的误差泛函j
f
(f)表示为：其中，表示相关算子，r
(i)
为第i次迭代的反褶积结果，||f||
q
为对点扩散函数引入的lq约束，μ1为其对应的超参，q∈(1，2]。10.如权利要求9所述的基于点扩散函数的地震数据高分辨处理方法，其特征在于：所述输出像域最小二乘偏移成像结果和更新后的点扩散函数，包括：判断点扩散函数与反褶积结果进行褶积与常规偏移成像结果误差是否小于给定的阈值；当点扩散函数与反褶积结果进行褶积与常规偏移成像结果误差小于给定阈值时输出像域最小二乘偏移成像结果和更新后的点扩散函数。
技术总结本发明公开了一种基于点扩散函数的地震数据高分辨处理方法，包括：输入地震记录、常规偏移结果和偏移速度，基于地震记录的道头信息获取观测系统；基于震源子波、旅行时和振幅计算初始点扩散函数；基于点扩散函数对偏移结果做反褶积处理，获取反褶积后的成像结果；基于成像结果和常规偏移结果，更新点扩散函数，输出像域最小二乘偏移成像结果和更新后的点扩散函数；本发明将强非线性的参数估计反问题转化为一个凸性较好的像域最小二乘反演成像问题；引入WKBJ近似下的渐进格林函数解决Hessian矩阵的计算量问题，进而计算地下任一点的点扩散函数；引入点扩散函数和反射系数的先验认识，构建点扩散函数和反射系数迭代更新的策略，得到更高分辨的成像结果。得到更高分辨的成像结果。得到更高分辨的成像结果。
技术研发人员：王华忠 盛燊 冯波
受保护的技术使用者：同济大学
技术研发日：2022.06.24
技术公布日：2022/11/1
