CN109856676A - 一种实现地震共反射面叠加参数优化的方法 - Google Patents
一种实现地震共反射面叠加参数优化的方法 Download PDFInfo
- Publication number
- CN109856676A CN109856676A CN201811451480.8A CN201811451480A CN109856676A CN 109856676 A CN109856676 A CN 109856676A CN 201811451480 A CN201811451480 A CN 201811451480A CN 109856676 A CN109856676 A CN 109856676A
- Authority
- CN
- China
- Prior art keywords
- indicate
- wave
- initial
- curvature radius
- radius
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Abstract
本发明提供一种实现地震共反射面叠加参数优化的方法，包括以下步骤：步骤S1，获取地震数据，基于地震数据获取初始反射点波前曲率半径、初始反射面波前曲率半径及倾角；步骤S2，通过所述初始反射点波前曲率半径、初始反射面波前曲率半径及倾角，基于目标函数，获取所述反射点波前曲率半径、初始反射面波前曲率半径及倾角的优化值。本发明提供一种实现地震共反射面叠加参数优化的方法，通过第一步先找出三个初始参数，搜索初始参数是单一逐个搜索，可使用简单的线性搜索；第二步是使用全局优化算法在初始参数附近搜索最优组合解，以实现CMAES算法解决CRS成像参数快速扫描问题。
Description
技术领域
本发明涉及地质勘测技术领域，更具体地，涉及一种实现地震共反射面叠加参数优化的方法。
背景技术
共反射面（common reflection surface，简称CRS）叠加是一种特殊的零偏移距的成像方法，具有不依赖于宏观速度模型和完全数据驱动的特点。由于CRS理论引入了反射面的概念,并将其限定在第一菲涅尔带范围内，从而能够在通过提高叠加次数来增强叠加剖面的信噪比的同时，又保持了一定的分辨率，所以CRS叠加方法更适应于低信噪比，低覆盖次数的地震资料。具体来说，共反射面元（CRS）叠加方法利用菲涅尔原理，扩大反射面元(对二维来说，扩大线元)增加覆盖次数，利用射线旁轴近似理论，将不同的共反射点道集校正到同一道集内，然后进行叠加，使反射信号的能量得以增强。
在二维情况下，需确定三个参数，相对确认CMP叠加速度来说稍微复杂，为了找出三个参数的最优组合，我们先找三参数初始值，再在初始值附件找出最优组合。初始值扫描使用线性扫描，以确定三参数的大概范围，之后再此基础之上需要使用全局优化算法来同时优化三参数。
二维条件下，依据旁轴射线旅行时射线理论（Schleicher等，1993），经典的CRS旅行时方程可以表达为（Jager等，2001）一个双曲型表达式，
(1)；
式中，为零炮检距旅行时，为零偏移距处的速度，为零偏移距对应的坐标，为炮检之间的中点坐标，为半偏移距，为零偏移距射线出射角，表示在IP波的曲率半径，表示在点N波的曲率半径。
传统的CMP叠加需要一个速度参数，而CRS叠加需要三个参数，此时参数搜索空间从一维上升到三维，搜索计算量成指数增长，计算量巨大。一维空间搜索使用线性搜索即可，但是三维空间需要更先进快速的优化搜索算法。目前使用较多的是模拟退火算法和基于梯度的算法，相对效率较低，容易陷入局部最小值，限制了CRS技术的推广应用。
发明内容
本发明提供一种实现地震共反射面叠加参数优化的方法，通过第一步先找出三个初始参数，搜索初始参数是单一逐个搜索，可使用简单的线性搜索；第二步是使用全局优化算法在初始参数附近搜索最优组合解，以实现CMAES算法解决CRS成像参数快速扫描问题。
根据本发明的一个方面，提供一种实现地震共反射面叠加参数优化的方法，包括以下步骤：
步骤S1，获取地震数据，基于地震数据获取初始反射点波前曲率半径、初始反射面波前曲率半径及初始倾角；
步骤S2，通过所述初始反射点波前曲率半径、初始反射面波前曲率半径及倾角，基于目标函数，获取所述波前曲率半径、反射面波前曲率半径及倾角的优化值。
在上述方案基础上优选，所述目标函数为：
；
且，；
；
；
其中，
表示待优化的反射点波前曲率半径；
表示待优化倾角；
表示待优化波前曲率半径的倒数，-1≤≤+1；
表示待优化的反射面波前曲率半径；
表示初始反射点波前曲率半径；
表示初始倾角；
表示初始反射面波前曲率半径的倒数，-1≤≤+1；
表示地震信号随时间变化的一致性函数；
A表示地震道振幅，是时间位置t的函数；
表示地震信号随时间变化的一致性组合函数；
M表示当前道集的道数；
K表示当前时间窗口内的地震信号样点数;
表示样点空间偏离初始参数模型的L2范数；
表示规则化的系数，0＜＜+1。
在上述方案基础上优选，所述位置时间t的计算表达式为：
；
其中：h表示偏移距；
表示待优化的反射点波前曲率半径；
表示待优化倾角；
表示垂直旅行双程时间；
表示地表速度。
本发明的一种实现地震共反射面叠加参数优化的方法，通过将目标函数设计成一致性函数与加权模型差的和的结构，使其在搜索参数时，在保证得到比较好的一致性条件下，又不偏离给定的初始模型参数（模型差异小），而且规则化参数，可用于调节偏离模型程度的大小。
本发明通过获取目标函数最小值时的反射点波前曲率半径、反射面波前曲率半径及倾角的优化值组合，使用CRS旅行时公式（1）就可以计算CRS旅行时间，就可以对CMP道集或超道集进行动较与叠加，生成最后的CRS叠加成果剖面，供地球物理学家和地质人员来做进一步的处理和解释。
附图说明
图1为本发明的某CDP点位置的三参数优化前与优化后的曲线图；
图2为本发明的某CDP位置动郊正时间随时间深度和偏移距变化图；
图3为本发明的单CDP叠加剖面图；
图4为本发明的CRS叠加剖面图；
图5为本发明的一种实现地震共反射面叠加参数优化的方法流程框图。
具体实施方式
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
请参阅图5所示，本发明提供了一种实现地震共反射面叠加参数优化的方法，包括以下步骤：
步骤S1，获取地震数据，基于地震数据计算获取初始反射点波前曲率半径、初始反射面波前曲率半径及初始倾角；
步骤S2，通过初始反射点波前曲率半径、初始反射面波前曲率半径及初始倾角，基于目标函数，获取所述反射点波前曲率半径、反射面波前曲率半径及倾角的优化值。
为了进一步的说明的本发明的技术方案，以下将详细的介绍本发明的步骤S1中，用于获取反射点波前曲率半径、反射面波前曲率半径及倾角的初始参数。
在常规地震数据处理中，我们需要人工拾取速度，而在CRS处理中需要的叠加速度就是常规叠加速度，但是测叠加速度可以依据CMP旅行时方程自动扫描实现。采用人工拾取的速度作为自动扫描的参照，在人工拾取速度的基础之上，扫描得到一致性函数最大值时对应的速度，这样可以保证得到的扫描速度的稳定性，同时又避免扫描到多次波的低速范围。
二维条件下，依据旁轴射线旅行时射线理论（Schleicher等，1993），经典的CRS旅行时方程可以表达为（Jager等，2001）一个双曲型表达式，
(1)；
其中，式中，为零炮检距旅行时，为零偏移距处的速度，为零偏移距对应的坐标，为炮检之间的中点坐标，为半偏移距，为零偏移距射线出射角，表示在点NIP波的曲率半径，表示在点N波的曲率半径。
由于CRS处理中叠加速度和直接相关，对于公式（1），假设，则该公式可变为：
（2）
对比叠加速度公式，
（3）
令，
（4）
变换后得到，
（5）
从上公式可以得出，获取和倾角，其他变量都是常数，可以计算。下一步获取扫描倾角。
在上一步叠加速度扫描之后，利用扫描的速度做叠加得到叠加剖面。注意倾角有正负，扫描前定义扫描方位如-40°到+40°。得到倾角后，结合第一步我们扫描得到的叠加速度，我们可以使用公式（5）计算。
同样，扫描N波的曲率半径，依然使用上一步得到的结果作为输入，即叠加剖面。使用简化后含的方程线性扫描N波的曲率半径。注意的范围是从负无穷到正无穷，可以通过线性扫描扫描其倒数而得到其初始值。
具体实施，是参数自动扫描速度，这里需要指出叠加速度在CRS处理中非常重要，但是无任何参数约束的速度模型自动扫描可能会存在很多问题，扫描速度前，提供参考速度，以参考速度为中心依据来设定扫描的范围，进行线性扫描，取一致性最好时对应的速度，再将得到了进行光滑，之后使用此速度做共中心点道集的动较和叠加，生成地震叠加剖面。生成剖面后，对剖面进行角度扫描，得到每个样点出的倾角值。倾角的扫描首先使用大角度间隔扫描，如2°左右，在确定大概角度范围后，使用小间隔进行二次扫描，如0.1°左右。在初始扫描得到速度和倾角后，依据公式（5）可以计算出NIP半径；假设公式（1）的h为0，即可扫描得到N波半径，取其倒数，极为N波曲率。
在获取三个初始值，倾角，N波的曲率半径，NIP波的曲率半径，接下来可以使用全局优化的算法，在初始值附近寻找最优参数组合。
其中，CMAES 是一种随机的与梯度无关的多维参数搜索的算法（Hansen，2011）。该算法并不基于梯度搜索，也没有导数计算过程，适合各种复杂搜索问题，如多局部最小问题，非线性非光滑问题，参数不可分离问题，不连续和含噪音问题等等。CMAES算法搜索是一个迭代的过程，使用随机采样方法不断的迭代更新随机变量，根据其特点可概括为三步：（1）依据样点均值和协方差矩阵对数据进行空间采样，（2）估算新的样点均值，（3）估算新的样点协方差。这三步依次不断的循环迭代更新，能够高效地搜索参数空间的多为多维解，而且收敛速度快。
这种特有的随机参数更新方式可用方程表示为，
（6）
（7）
（8）
~：表示随机变量分布；
：表示多维随机变量正态分布，均值为0，协方差矩阵；
：表示从第g+1代的第k个后代；
：表示第g代的均值；
：表示第g代总体方差或步长；
：表示第g代协方差矩阵；
：表示样点总数，通常大于或等于2；
：表示加权系数空间；
：表示空间样本；
：表示第g+1代空间实验协方差。
CRS叠加参数的搜索本质上就是一个含噪音，非线性，不连续，不可分离的参数优化问题，因而CMAES是一种适合于CRS参数搜索的非常好的优化方法。
使用CMAES搜索首先第一步是定义最小化目标函数，本发明的目标函数为：
（9）；
且， （10）；
（11）；
(12)；
其中，
表示待优化的反射点波前曲率半径；
表示待优化倾角；
表示待优化波前曲率半径的倒数，-1≤≤+1；
表示待优化的反射面波前曲率半径；
表示初始反射点波前曲率半径；
表示初始倾角；
表示初始反射面波前曲率半径的倒数，-1≤≤+1；
表示地震信号随时间变化的一致性函数；
A表示地震道振幅，是时间位置t的函数；
表示地震信号随时间变化的一致性组合函数；
M表示当前道集的道数；
K表示当前时间窗口内的地震信号样点数;
表示样点空间偏离初始参数模型的L2范数；
表示规则化的系数，0＜＜+1。
其中，本发明的位置时间t的计算表达式为：
；
其中：h表示偏移距；
表示待优化的反射点波前曲率半径；
表示待优化倾角；
表示垂直旅行双程时间；
表示地表速度。
以下将结合实施例具体说明如何采用CMAES基于目标函数来扫描优化三个初始参数。
在三个初始参数计算出来后，使用CMAES基于目标函数（9）来扫描优化后的三参数，这里以实例说明。图1是某CDP点位置的三参数优化前（蓝色）优化后（红色）的曲线图。相比较而言，红色曲线跳动更多，这个是由于CMAES是一种随机算法，为了寻找最优参数组合，搜索的参数在初始参数附近寻优达到更小的目标函数值。
在搜索到最优三个参数组合后，接下来是使用公式（1）计算得到曲面的旅行时间，图2显示的是某CDP位置动郊正时间图形，x轴是偏移距，图2中的时间参数是经过光滑后的，由于是随机的求解过程，计算的三参数含有随机成分，因而计算的时间需要做光滑。
得到了光滑后的动校正时间，即可对CDP道集数据进行动校正和叠加，图3显示的是单CDP叠加剖面图，图4是CRS叠加剖面图，属于本发明的最新成果，可以看出，图4的信号得到的增强，噪音得到了压制，图4成像效果得到了极大的改善。
最后，本申请的方法仅为较佳的实施方案，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (3)
1.一种实现地震共反射面叠加参数优化的方法，其特征在于，包括以下步骤：
步骤S1，获取地震数据，基于地震数据计算获得初始反射点波前曲率半径、初始反射面波前曲率半径及初始倾角；
步骤S2，通过所述初始反射点波前曲率半径、初始反射面波前曲率半径及初始倾角，基于目标函数，获取所述波前曲率半径、反射面波前曲率半径及倾角的优化值。
2.如权利要求1所述的一种实现地震共反射面叠加参数优化的方法，其特征在于，所述目标函数为：
；
且，；
；
；
其中，
表示待优化的反射点波前曲率半径；
表示待优化倾角；
表示待优化波前曲率半径的倒数，-1≤≤+1；
表示待优化的反射面波前曲率半径；
表示初始反射点波前曲率半径；
表示初始倾角；
表示初始反射面波前曲率半径的倒数，-1≤≤+1；
表示地震信号随时间变化的一致性函数；
A表示地震道振幅，是时间位置t的函数；
表示地震信号随时间变化的一致性组合函数；
M表示当前道集的道数；
K表示当前时间窗口内的地震信号样点数;
表示样点空间偏离初始参数模型的L2范数；
表示规则化的系数，0＜＜+1。
3.如权利要求2所述的一种实现地震共反射面叠加参数优化的方法，其特征在于，所述位置时间t的计算表达式为：
；
其中：h表示偏移距；
表示待优化的反射点波前曲率半径；
表示待优化倾角；
表示垂直旅行双程时间；
表示地表速度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811451480.8A CN109856676B (zh) | 2018-11-30 | 2018-11-30 | 一种实现地震共反射面叠加参数优化的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811451480.8A CN109856676B (zh) | 2018-11-30 | 2018-11-30 | 一种实现地震共反射面叠加参数优化的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109856676A true CN109856676A (zh) | 2019-06-07 |
| CN109856676B CN109856676B (zh) | 2021-03-12 |
Family
ID=66890441
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811451480.8A Active CN109856676B (zh) | 2018-11-30 | 2018-11-30 | 一种实现地震共反射面叠加参数优化的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109856676B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110954955A (zh) * | 2019-12-23 | 2020-04-03 | 西南石油大学 | 一种基于进化优化算法的地震随机反演方法 |
| CN112327364A (zh) * | 2020-11-02 | 2021-02-05 | 中国石油大学(华东) | 一种改进型一步法三维crs叠加方法 |
| CN113960668A (zh) * | 2021-10-21 | 2022-01-21 | 中国石油化工股份有限公司 | 基于叠前时间偏移的增强反射信息的方法及装置 |
| CN113960669A (zh) * | 2021-10-21 | 2022-01-21 | 中国石油化工股份有限公司 | 基于共成像点道集组合计算的反射信息增强方法及装置 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4878205A (en) * | 1987-07-23 | 1989-10-31 | Boris Gelchinsky | Method of analyzing the structure of a medium, particularly useful for seismic prospecting |
| CN105301648A (zh) * | 2014-07-31 | 2016-02-03 | 中国石油化工股份有限公司 | 一种获取共反射面元叠加参数的方法 |
| CN106483559A (zh) * | 2015-08-24 | 2017-03-08 | 中国石油化工股份有限公司 | 一种地下速度模型的构建方法 |
| KR20170040569A (ko) * | 2015-10-05 | 2017-04-13 | 주식회사 에프에스 | 지층구조의 해석력 향상 기능을 갖는 공통반사표면과 공통회절표면법을 이용한 탄성파 데이터의 복합처리방법 |
| CN107179543A (zh) * | 2016-03-11 | 2017-09-19 | 中国石油化工股份有限公司 | 对叠前数据进行规则化的方法和装置 |
| CN109557585A (zh) * | 2019-01-25 | 2019-04-02 | 中国石油大学(华东) | 一种基于差分进化算法的crs倾角分解方法 |
-
2018
- 2018-11-30 CN CN201811451480.8A patent/CN109856676B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4878205A (en) * | 1987-07-23 | 1989-10-31 | Boris Gelchinsky | Method of analyzing the structure of a medium, particularly useful for seismic prospecting |
| CN105301648A (zh) * | 2014-07-31 | 2016-02-03 | 中国石油化工股份有限公司 | 一种获取共反射面元叠加参数的方法 |
| CN106483559A (zh) * | 2015-08-24 | 2017-03-08 | 中国石油化工股份有限公司 | 一种地下速度模型的构建方法 |
| KR20170040569A (ko) * | 2015-10-05 | 2017-04-13 | 주식회사 에프에스 | 지층구조의 해석력 향상 기능을 갖는 공통반사표면과 공통회절표면법을 이용한 탄성파 데이터의 복합처리방법 |
| CN107179543A (zh) * | 2016-03-11 | 2017-09-19 | 中国石油化工股份有限公司 | 对叠前数据进行规则化的方法和装置 |
| CN109557585A (zh) * | 2019-01-25 | 2019-04-02 | 中国石油大学(华东) | 一种基于差分进化算法的crs倾角分解方法 |
Non-Patent Citations (2)
| Title |
|---|
| 杨锴等: "共反射面元叠加的应用实践", 《地球物理学报》 * |
| 黄亚飞 等: "求解全局优化问题的正交协方差矩阵自适应进化策略算法", 《计算机应用》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110954955A (zh) * | 2019-12-23 | 2020-04-03 | 西南石油大学 | 一种基于进化优化算法的地震随机反演方法 |
| CN112327364A (zh) * | 2020-11-02 | 2021-02-05 | 中国石油大学(华东) | 一种改进型一步法三维crs叠加方法 |
| CN112327364B (zh) * | 2020-11-02 | 2022-04-29 | 中国石油大学(华东) | 一种改进型一步法三维crs叠加方法 |
| CN113960668A (zh) * | 2021-10-21 | 2022-01-21 | 中国石油化工股份有限公司 | 基于叠前时间偏移的增强反射信息的方法及装置 |
| CN113960669A (zh) * | 2021-10-21 | 2022-01-21 | 中国石油化工股份有限公司 | 基于共成像点道集组合计算的反射信息增强方法及装置 |
| CN113960668B (zh) * | 2021-10-21 | 2024-04-16 | 中国石油化工股份有限公司 | 基于叠前时间偏移的增强反射信息的方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109856676B (zh) | 2021-03-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109856676A (zh) | 一种实现地震共反射面叠加参数优化的方法 | |
| US10330807B2 (en) | Device and method for displaying full azimuth angle domain image data | |
| Gizon et al. | Local helioseismology | |
| Wu et al. | Directional illumination analysis using beamlet decomposition and propagation | |
| Medvedev et al. | Relation of internal gravity wave anisotropy with neutral wind characteristics in the upper atmosphere | |
| CN106597540B (zh) | 高斯束偏移成像方法及装置 | |
| CN106094029A (zh) | 利用偏移距矢量片地震数据预测储层裂缝的方法 | |
| CN112883564B (zh) | 一种基于随机森林的水体温度预测方法及预测系统 | |
| CN105607122B (zh) | 一种基于全变分地震数据分解模型的地震纹理提取与增强方法 | |
| Abbad et al. | Automatic nonhyperbolic velocity analysis | |
| Xiang‐Bo et al. | Anisotropic Radon transform and its application to demultiple | |
| Liu et al. | Navigability analysis of local gravity map with projection pursuit-based selection method by using gravitation field algorithm | |
| CN106257309A (zh) | 叠后地震数据体处理方法及装置 | |
| Coimbra et al. | Offset-continuation stacking: Theory and proof of concept | |
| CN109581488A (zh) | 一种基于粒子群优化crs超道集地震波成像方法 | |
| CN109188516A (zh) | Radon域能量扫描叠加的微地震事件定位方法 | |
| CN109387872A (zh) | 表面多次波预测方法 | |
| CN112904407B (zh) | 复杂地形与干扰条件下的微动勘探方法 | |
| CN106569269B (zh) | 一种优选道加权速度谱计算方法 | |
| CN109655888A (zh) | 地震数据处理中光滑浮动基准面的定量选择方法及系统 | |
| Tang et al. | A stable and effective strategy for angle gather extraction from reverse time migration using traveltime gradient | |
| Wu et al. | Stacking common-image gathers via the flatten-illumination-similarity method | |
| Zhuang et al. | A computational method for wide-azimuth 3D dip-angle gathers using Gaussian beam migration | |
| Leung et al. | Liouville partial-differential-equation methods for computing 2D complex multivalued eikonals in attenuating media | |
| Zang et al. | Calibration method of petroleum underground prospects based on high precision gravity and magnetic exploration |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |