CN111812708B - 地震波成像方法及装置 - Google Patents
地震波成像方法及装置 Download PDFInfo
- Publication number
- CN111812708B CN111812708B CN201910287978.3A CN201910287978A CN111812708B CN 111812708 B CN111812708 B CN 111812708B CN 201910287978 A CN201910287978 A CN 201910287978A CN 111812708 B CN111812708 B CN 111812708B
- Authority
- CN
- China
- Prior art keywords
- stratum
- seismic wave
- wave imaging
- seismic
- model
- 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.)
- Active
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 220
- 238000005457 optimization Methods 0.000 claims abstract description 88
- 238000000034 method Methods 0.000 claims abstract description 24
- 230000015572 biosynthetic process Effects 0.000 claims description 38
- 239000002245 particle Substances 0.000 claims description 22
- 238000004422 calculation algorithm Methods 0.000 claims description 21
- 238000004590 computer program Methods 0.000 claims description 13
- 238000001514 detection method Methods 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 12
- 230000006870 function Effects 0.000 claims description 8
- 238000003860 storage Methods 0.000 claims description 7
- 238000010586 diagram Methods 0.000 description 14
- 238000004364 calculation method Methods 0.000 description 12
- 230000004044 response Effects 0.000 description 6
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 150000003839 salts Chemical class 0.000 description 5
- 230000010363 phase shift Effects 0.000 description 4
- 238000007796 conventional method Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000005012 migration Effects 0.000 description 3
- 238000013508 migration Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 241000764238 Isis Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种地震波成像方法及装置，该方法包括：根据地层速度模型，建立地震波成像参数优化模型；对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；根据多个地层的成像数据，获得地震波成像结果。本发明可以进行地震波成像，成像精度高。
Description
技术领域
本发明涉及一种地震波成像技术，尤其涉及一种地震波成像方法方法及装置。
背景技术
地震波成像技术是地震数据处理的重要环节，对于油气资源勘探至关重要。单程波偏移算法作为地震偏移成像的一种重要分支，是一种将声波方程分解为上行波方程和下行波方程分别进行波场延拓计算的一种方法。现有技术包括在假设速度只在深度方向变化的情况下的单程波相移偏移算法。单程波相移偏移算法在介质存在横向速度变化的情况下，应用效果并不理想。为解决该问题，现有技术提出了相移加插值偏移算法，该算法在波场深度延拓计算中选择多个参考速度分别进行相移计算然后利用数学工具进行插值计算。另外，现有技术还提出了裂步傅里叶偏移算法，该算法的实现是在频率-空间和频率-空间域的交替进行的，通过在频率-空间域中一个相位校正来实现在复杂速度变化介质中的延拓计算。为了适应强横向速度变化介质下的构造准确成像，现有技术提出了傅里叶有限差分偏移算法。相对于裂步傅里叶偏移算法，傅里叶有限差分偏移算法的成像角度得到了大大的增加，对于复杂构造的成像质量也提升不少。现在技术根据扰动理论推导了广义屏传播算子，在同性介质中对比了各阶广义屏算子与理论波场的误差并将高阶广义屏传播算子推广到了TI介质中，建立高阶广义屏比较完整的理论基础。
在相位屏传播算子研究的基础上，为了准确描述地震波场在复杂介质中的传播，发展了单程波高阶广义屏传播算子。常规的高阶广义屏传播算子包括由Taylor级数展开确定的系数、地层的速度，在利用广义屏传播算子进行偏移成像的时候，一般选择速度模型的最小速度为参考速度。当速度的横向变化比较大时，参考速度与最大速度的差异比较大，图1为利用广义屏传播算子计算的相位信息与准确的相位值。从图1中可见，当参考速度只有模型速度的50％时，随着传播角度的增加，传统广义屏传播算子的计算误差也随之增大；即使由Taylor级数展开确定的系数采用高阶展开式(20阶)，传统广义屏传播算子的计算误差也比较大，即精度低，而且采用高阶展开式也会增加计算成本，不利于实际应用。
为了克服地震波成像时常规单程波传播算子在横向速度变化较大介质中大角度成像的问题，现有技术提出了基于粒子群优化算法的高阶广义屏偏移成像。该技术方案没有采用常规的Taylor级数展开式计算高阶广义屏传播算子中的系数，而是采用粒子群优化算法在全局最优的思想下计算高阶广义屏传播算子中的系数。但该技术虽然对常规的广义屏传播算法进行了一定程度的改进，使得计算精度得到了一定程度的提高，但是在计算高阶广义屏传播算子中的系数系数的时候并没有考虑实际地层的速度变化情况，无法根据实际速度变化情况给出最优解，成像精度还是较低。
综上所述，现有技术在进行地震波成像时的精度低。
发明内容
本发明实施例提出一种地震波成像方法，用以进行地震波成像，成像精度高，该方法包括：
根据地层速度模型，建立地震波成像参数优化模型；
对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；
根据多个地层的成像数据，获得地震波成像结果。
本发明实施例提出一种地震波成像装置，用以进行地震波成像，成像精度高，该装置包括：
优化模型获得模块，用于根据地层速度模型，建立地震波成像参数优化模型；
地震波成像数据获得模块，用于对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；
地震波成像结果获得模块，用于根据多个地层的成像数据，获得地震波成像结果。
本发明实施例还提出了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述地震波成像方法。
本发明实施例还提出了一种计算机可读存储介质，所述计算机可读存储介质存储有执行上述地震波成像方法的计算机程序。
在本发明实施例中，根据地层速度模型，建立地震波成像参数优化模型；对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；根据多个地层的成像数据，获得地震波成像结果。本发明实施例中，获得了考虑地层速度的地震波成像参数优化模型，通过将每一地层的速度输入至该地震波成像参数优化模型，求解该地震波成像参数优化模型，可获得优化后的地震波成像参数，根据优化后的地震波成像参数，可获得最终的地震波成像结果，由于上述过程考虑了实际的地层速度，且对每一地层的地震波成像参数进行了优化，因此最后获得的地震波成像结果的精度高。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为利用广义屏传播算子计算的相位信息与准确的相位值；
图2为本发明实施例中地震波成像方法的流程图；
图3为本发明实施例中地震波成像方法的详细流程图；
图4为本发明实施例中一个横向速度变化较大的地层速度模型；
图5为优化后的广义屏传播算子对应的脉冲响应；
图6为常规广义屏传播算子对应的脉冲响应；
图7为盐丘速度模型的示意图；
图8为采用本发明提出的方法获得的地震波成像结果；
图9为采用常规方法得到的地震波成像方法；
图10为本发明实施例中地震波成像装置的示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
图2为本发明实施例中地震波成像方法的流程图，如图2所示，该方法包括：
步骤201，根据地层速度模型，建立地震波成像参数优化模型；
步骤202，对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；
步骤203，根据多个地层的成像数据，获得地震波成像结果。
本发明实施例中，获得了考虑地层速度的地震波成像参数优化模型，通过将每一地层的速度输入至该地震波成像参数优化模型，求解该地震波成像参数优化模型，可获得优化后的地震波成像参数，根据优化后的地震波成像参数，可获得最终的地震波成像结果，由于上述过程考虑了实际的地层速度，且对每一地层的地震波成像参数进行了优化，因此最后获得的地震波成像结果的精度高。
在一实施例中，在根据地层速度模型，建立地震波成像参数优化模型之前，还包括：
根据地层速度模型，获得地震波的广义屏传播算子模型；
根据地层速度模型，建立地震波成像参数优化模型，包括：
根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型。
在一实施例中，地震波成像参数为广义屏传播算子；
求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数，包括：
求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值；
根据该地层的广义屏传播算子的系数的最优值，获得该地层的优化后的广义屏传播算子。
在一实施例中，根据地层速度模型，建立地震波成像参数优化模型包括：
采用如下公式，根据地层速度模型，获得地震波的广义屏传播算子模型：
其中，Pi GSP为第i层的广义屏传播算子；
vi＝[v1,i,v2.,i,…,vM,i]为第i层的地层速度，M为第i层地层点数；
vr,i为地层速度模型中第i层的参考地层速度，为vi的最小值或平均值；
ω为角频率；
kx,ky分别为x方向和y方向的波数；
aj,i为第i层的广义屏传播算子的第j个系数；n为系数的总数；
将公式(1)进行归一化，可以写成：
可知公式(2)是以
为未知数xi求取的泰勒级数高阶展开的近似式。
其中，
的泰勒级数近似式如下：
分析可知广义屏传播算子；的计算精度很大程度上取决于
近似式的计算精度。
因此，采用如下公式，根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型：
其中，Ii(a1,a2，…,aj，…,an)为地层的地震波成像参数优化模型的第i层的目标函数；
xMIN,i和xMAX,i分别为以第i层的地层速度计算得到的
的最小值和最大值。
当公式(4)的目标函数的值最小时，可以获得a1,i,a2,i，…,aj,i，…,an,i中的最优值，以求取更好的逼近效果。
在一实施例中，采用粒子群算法，求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值。
具体实施时，粒子群算法的具体公式如下：
Vj(t+1)＝ωVj(t)+c1 r1(Pj(t)-aj(t))+c2 r2(Pg(t)-aj(t)) (7)
aj(t+1)＝aj(t)+Vj(t+1) (8)
Vj(t+1)＝Vmax,Vj(t+1)＞Vmax (9)
Vj(t+1)＝Vmin,Vj(t+1)＜Vmin (10)
其中，aj(t)为粒子j在搜索次数t时的广义屏算子的系数；
Vj(t)为粒子j在搜索次数t时的搜索增量；
ω为惯性权值；
r1和r2为学习因子；
Vmax和Vmin为每个粒子限制搜索的最大速度和最小速度(并非地层速度)；
以上惯性权值、学习因子、每个粒子限制搜索的最大速度和最小速度均为粒子群算法的初始化参数；公式(7)和公式(8)分别为速度更新公式和位置更新公式。
采用粒子群算法，求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值的具体过程如下：
首先，初始化粒子群算法中的初始化参数；
然后，迭代更新粒子群的速度和位置，直至达到终止条件；
以公式(9)和公式(10)为约束条件，根据公式(7)和公式(8)更新粒子群的速度和位置。所述终止条件包括包括设定迭代次数或粒子j在搜索次数t时的广义屏算子的系数或搜索增量或最大迭代次数达到了设定值。
当然，可以理解的是，以上粒子群算法仅为举例，还可以采用其他可以获得该地层的广义屏传播算子的系数的最优值的方法，例如遗传算法等智能优化算法，相关变化例均应落入本发明的保护范围。
在一实施例中，根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据，包括：
获得该地层的检波点波场数据和震源波场数据；
根据该地层的优化后的广义屏传播算子，对该地层的检波点波场数据和震源波场数据进行延拓处理；
根据该地层的延拓处理后的检波点波场数据和震源波场数据，进行该地层的地震波成像，获得该地层的地震波成像数据。
图3为本发明实施例中地震波成像方法的详细流程图，如图3所示，该方法包括：
步骤301，根据地层速度模型，建立地震波成像参数优化模型；
步骤302，将第i层的地层速度输入至地震波成像参数优化模型；
步骤303，求解第i层的地层对应的地震波成像参数优化模型，获得第i层的广义屏传播算子的系数的最优值；
步骤304，根据第i层的广义屏传播算子的系数的最优值，获得第i层的优化后的广义屏传播算子；
步骤305，获得第i层的检波点波场数据和震源波场数据；
步骤306，根据第i层的优化后的广义屏传播算子，对第i层的检波点波场数据和震源波场数据进行延拓处理；
步骤307，根据第i层的延拓处理后的检波点波场数据和震源波场数据，进行第i层的地震波成像，获得第i层的地震波成像数据；
步骤308，判断第i层是否为最大深度的地层，若第i层为最大深度的地层，执行步骤310，否则执行步骤309；
步骤309，令i＝i+1，执行步骤302；
步骤310，根据所有地层的成像数据，获得地震波成像结果。
下面给出一具体实施例，说明采用本发明提出的方法计算的优化后的广义屏传播算子的优点。
图4为本发明实施例中一个横向速度变化较大的地层速度模型，如图4所示，该地层速度模型的网格间距dx＝5.0m，dz＝5.0m，时间采用间隔dt＝0.001s，时间采样点为1000个，用于计算的理论子波主频f0＝30Hz，计算的频率范围为1Hz-80Hz参考地层速度取地层速度的最小值。
首先根据该地层速度模型建立地震波成像参数优化模型，获得公式(4)-公式(6)表示的地震波成像参数优化模型。
然后，对每一地层，将该地层的地层速度输入至公式(4)-公式(6)表示的地震波成像参数优化模型，a1,i,a2,i，…,aj,i，…,an,i中的n＝4；采用公式(7)-公式(10)表示的粒子群算法，求解地震波成像参数优化模型，其中，公式(7)-公式(10)表示的粒子群算法中，学习因子r1和r2均为1.5，惯性权重ω＝0.7，最大迭代次数为200次，最后可得到粒子j在搜索次数t时的广义屏算子的系数的最优值，根据广义屏算子的系数的最优质可最终计算得到优化后的广义屏传播算子。图5为优化后的广义屏传播算子对应的脉冲响应，图6为常规广义屏传播算子对应的脉冲响应，如图5和图6所示，在强横向速度变化的速度模型(即横向速度变化较大的速度模型)中，采用本发明实施例计算的优化后的广义屏传播算子的脉冲响应与准确的脉冲响应在大角度区域吻合的更好，体现了本发明实施例相对于常规广义屏传播算子在大角度成像上的效果更好。
表1为本实施例中计算的广义屏算子的系数的最优值(第1-10层)，表2为常规广义屏传播算子的系数的值(第1-10层)，两个表的数值是不一样的，从而可以体现出图5和图6的差异。
表1本实施例中计算的广义屏算子的系数的最优值(第1-10层)
| 地层层数 | a<sub>1</sub> | a<sub>2</sub> | a<sub>3</sub> | a<sub>4</sub> |
| 第1层 | 0.5605 | 0.5799 | -0.2185 | -0.8707 |
| 第2层 | -0.8203 | 0.7807 | -0.1247 | -0.8026 |
| 第3层 | 0.9444 | 0.7709 | -0.7480 | -0.6223 |
| 第4层 | 0.4817 | 0.2718 | -0.0057 | -0.8778 |
| 第5层 | 0.4810 | -0.4943 | 0.3761 | -0.8109 |
| 第6层 | 0.5011 | -0.3265 | 0.5538 | -0.9889 |
| 第7层 | 0.0578 | 0.6784 | 0.0025 | -0.9811 |
| 第8层 | 0.8357 | 0.3077 | -0.0270 | -0.9488 |
| 第9层 | 0.5000 | 0.3160 | -0.7800 | -0.2912 |
| 第10层 | -0.5631 | 0.5363 | 0.2360 | -0.9852 |
表2常规广义屏传播算子的系数的值(第1-10层)
| 地层层数 | a<sub>1</sub> | a<sub>2</sub> | a<sub>3</sub> | a<sub>4</sub> |
| 第1层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第2层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第3层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第4层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第5层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第6层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第7层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第8层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第9层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第10层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
下面给出另一实施例，说明本发明提出的地震波成像方法的具体应用。
本实施例利用盐丘速度模型来说明本发明提出的地震波成像方法在强横向速度变化介质模型中的成像差异，图7为盐丘速度模型的示意图。盐丘速度模型的特点是在模型的中间存在一个高速体，其速度为4500m/s，该高速体的速度为周围围岩速度的两倍，因此对高速体下方构造的成像是本实施例的重点。在本实施例中速度模型的网格间距dx＝20.0m，dz＝4.0m，时间采用间隔dt＝0.001s，时间采样点为3000个，计算的频率范围为4Hz-40Hz，参考地层速度取地层速度的最小值。
首先根据该地层速度模型建立地震波成像参数优化模型，获得公式(4)-公式(6)表示的地震波成像参数优化模型。
然后，对每一地层，将该地层的地层速度输入至公式(4)-公式(6)表示的地震波成像参数优化模型，采用公式(7)-公式(10)表示的粒子群算法，求解地震波成像参数优化模型，其中，公式(7)-公式(10)表示的粒子群算法中，学习因子c1＝1.5,c2＝1.5，惯性权重ω＝0.7，最大迭代次数为200次。
最后可得到粒子j在搜索次数t时的广义屏算子的系数的最优值。
表3为本实施例中计算的广义屏算子的系数的最优值(第501-510层)，表4为常规广义屏传播算子的系数的值(第501-510层)。
表3本实施例中计算的广义屏算子的系数的最优值(第501-510层)
| 地层层数 | a<sub>1</sub> | a<sub>2</sub> | a<sub>3</sub> | a<sub>4</sub> |
| 第501层 | -0.8636 | 0.8403 | 0.0645 | -0.9073 |
| 第502层 | 0.9253 | -0.3008 | 0.4885 | -0.9830 |
| 第503层 | -0.4187 | -0.3418 | 0.8850 | -0.9984 |
| 第504层 | 0.7386 | 0.4873 | -0.1035 | -0.9026 |
| 第505层 | 0.8749 | 0.0550 | -0.1010 | -0.7118 |
| 第506层 | -0.6565 | 0.6414 | -0.0841 | -0.6995 |
| 第507层 | 0.6562 | 0.2156 | 0.1757 | -0.9629 |
| 第508层 | -0.5149 | 0.9233 | 0.0078 | -0.9944 |
| 第509层 | 0.7862 | 0.0621 | 0.1105 | -0.8643 |
| 第510层 | 0.5296 | 0.9497 | -0.4227 | -0.8652 |
表4常规广义屏传播算子的系数的值(第501-510层)
| 地层层数 | a<sub>1</sub> | a<sub>2</sub> | a<sub>3</sub> | a<sub>4</sub> |
| 第501层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第502层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第503层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第504层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第505层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第506层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第507层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第508层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第509层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
| 第510层 | -0.5000 | -0.1250 | -0.0625 | -0.0391 |
根据广义屏算子的系数的最优质可最终计算得到优化后的广义屏传播算子；
获得该地层的检波点波场数据和震源波场数据；
根据该地层的优化后的广义屏传播算子，对该地层的检波点波场数据和震源波场数据进行延拓处理；
根据该地层的延拓处理后的检波点波场数据和震源波场数据，进行该地层的地震波成像，获得该地层的地震波成像数据。
图8为采用本发明提出的方法获得的地震波成像结果，图9为采用常规方法得到的地震波成像方法，图8和图9中箭头处为高速体下方的一个断层，可以看到：与常规方法相比，本发明提出的方法获得的高速体下方的断层断面的清晰度更高，即成像精度更好；图8和图9中的虚线框为高速体下方的水平底界面，可以看到，本发明提出的方法可减小盐丘这类高速介质下方成像同向轴扭曲的现象，使得成像效果更稳定，成像效果更好。
在本发明实施例提出的地震波成像方法中，根据地层速度模型，建立地震波成像参数优化模型；对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；根据多个地层的成像数据，获得地震波成像结果。本发明实施例中，获得了考虑地层速度的地震波成像参数优化模型，通过将每一地层的速度输入至该地震波成像参数优化模型，求解该地震波成像参数优化模型，可获得优化后的地震波成像参数，根据优化后的地震波成像参数，可获得最终的地震波成像结果，由于上述过程考虑了实际的地层速度，且对每一地层的地震波成像参数进行了优化，因此最后获得的地震波成像结果的精度高。
基于同一发明构思，本发明实施例还提供了一种地震波成像装置，如下面的实施例所述。由于这些解决问题的原理与地震波成像方法相似，因此装置的实施可以参见方法的实施，重复之处不在赘述。
图10为本发明实施例中地震波成像装置的示意图，如图10所示，该装置包括：
优化模型获得模块1001，用于根据地层速度模型，建立地震波成像参数优化模型；
地震波成像数据获得模块1002，用于对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；
地震波成像结果获得模块1003，用于根据多个地层的成像数据，获得地震波成像结果。
在一实施例中，该装置还包括广义屏传播算子模型获得模块1004，用于：根据地层速度模型，获得地震波的广义屏传播算子模型；
优化模型获得模块1001具体用于：根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型。
综上所述，在本发明实施例提出的地震波成像装置中，优化模型获得模块，用于根据地层速度模型，建立地震波成像参数优化模型；地震波成像数据获得模块，用于对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；地震波成像结果获得模块，用于根据多个地层的成像数据，获得地震波成像结果。本发明实施例中，获得了考虑地层速度的地震波成像参数优化模型，通过将每一地层的速度输入至该地震波成像参数优化模型，求解该地震波成像参数优化模型，可获得优化后的地震波成像参数，根据优化后的地震波成像参数，可获得最终的地震波成像结果，由于上述过程考虑了实际的地层速度，且对每一地层的地震波成像参数进行了优化，因此最后获得的地震波成像结果的精度高。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (6)
1.一种地震波成像方法，其特征在于，包括：
根据地层速度模型，建立地震波成像参数优化模型；
对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；
根据多个地层的成像数据，获得地震波成像结果；
在根据地层速度模型，建立地震波成像参数优化模型之前，还包括：根据地层速度模型，获得地震波的广义屏传播算子模型；
根据地层速度模型，建立地震波成像参数优化模型，包括：根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型；
根据地层速度模型，建立地震波成像参数优化模型包括：
采用如下公式，根据地层速度模型，获得地震波的广义屏传播算子模型：
其中，Pi GSP为第i层的广义屏传播算子；
vi＝[v1,i,v2.,i,…,vM,i]为第i层的地层速度，M为第i层地层点数；
vr,i为地层速度模型中第i层的参考地层速度，为vi的最小值或平均值；
ω为角频率；
kx,ky分别为x方向和y方向的波数；
aj,i为第i层的广义屏传播算子的第j个系数；n为系数的总数；
采用如下公式，根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型：
其中，Ii(a1,i,a2,i，…,aj,i，…,an,i)为地震波成像参数优化模型的第i层的目标函数；
xMIN,i和xMAX,i分别为以第i层的地层速度计算得到的
的最小值和最大值；
地震波成像参数为广义屏传播算子；
求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数，包括：求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值；根据该地层的广义屏传播算子的系数的最优值，获得该地层的优化后的广义屏传播算子。
2.如权利要求1所述的地震波成像方法，其特征在于，求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值，包括：
采用粒子群算法，求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值。
3.如权利要求1所述的地震波成像方法，其特征在于，根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据，包括：
获得该地层的检波点波场数据和震源波场数据；
根据该地层的优化后的广义屏传播算子，对该地层的检波点波场数据和震源波场数据进行延拓处理；
根据该地层的延拓处理后的检波点波场数据和震源波场数据，进行该地层的地震波成像，获得该地层的地震波成像数据。
4.一种地震波成像装置，其特征在于，包括：
优化模型获得模块，用于根据地层速度模型，建立地震波成像参数优化模型；
地震波成像数据获得模块，用于对每一地层，将该地层的地层速度输入至地震波成像参数优化模型，求解地震波成像参数优化模型，获得该地层的优化后的地震波成像参数；根据该地层的优化后的地震波成像参数，获得该地层的地震波成像数据；
地震波成像结果获得模块，用于根据多个地层的成像数据，获得地震波成像结果；
所述装置还包括广义屏传播算子模型获得模块，用于：根据地层速度模型，获得地震波的广义屏传播算子模型；
优化模型获得模块具体用于：根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型；
采用如下公式，根据地层速度模型，获得地震波的广义屏传播算子模型：
其中，Pi GSP为第i层的广义屏传播算子；
vi＝[v1,i,v2.,i,…,vM,i]为第i层的地层速度，M为第i层地层点数；
vr,i为地层速度模型中第i层的参考地层速度，为vi的最小值或平均值；
ω为角频率；
kx,ky分别为x方向和y方向的波数；
aj,i为第i层的广义屏传播算子的第j个系数；n为系数的总数；
采用如下公式，根据地层速度模型和地震波的广义屏传播算子模型，建立地震波成像参数优化模型：
其中，Ii(a1,i,a2,i，…,aj,i，…,an,i)为地震波成像参数优化模型的第i层的目标函数；
xMIN,i和xMAX,i分别为以第i层的地层速度计算得到的
的最小值和最大值；
地震波成像参数为广义屏传播算子；
地震波成像数据获得模块具体用于：求解地震波成像参数优化模型，获得该地层的广义屏传播算子的系数的最优值；根据该地层的广义屏传播算子的系数的最优值，获得该地层的优化后的广义屏传播算子。
5.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至3任一所述方法。
6.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有执行权利要求1至3任一所述方法的计算机程序。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910287978.3A CN111812708B (zh) | 2019-04-11 | 2019-04-11 | 地震波成像方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910287978.3A CN111812708B (zh) | 2019-04-11 | 2019-04-11 | 地震波成像方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111812708A CN111812708A (zh) | 2020-10-23 |
| CN111812708B true CN111812708B (zh) | 2022-08-30 |
Family
ID=72844298
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910287978.3A Active CN111812708B (zh) | 2019-04-11 | 2019-04-11 | 地震波成像方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111812708B (zh) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105911587A (zh) * | 2016-04-22 | 2016-08-31 | 中国地质大学（北京） | 一种利用单程波算子的双程波叠前深度偏移方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102313902B (zh) * | 2011-06-01 | 2013-06-19 | 西南石油大学 | 一种基于切比雪夫展开的广义屏叠前深度偏移方法 |
| KR101549388B1 (ko) * | 2014-10-17 | 2015-09-02 | 한국지질자원연구원 | 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법 |
-
2019
- 2019-04-11 CN CN201910287978.3A patent/CN111812708B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105911587A (zh) * | 2016-04-22 | 2016-08-31 | 中国地质大学（北京） | 一种利用单程波算子的双程波叠前深度偏移方法 |
Non-Patent Citations (3)
| Title |
|---|
| High-order generalized screen propagator migration based on particle swarm optimization;He Run et al.;《Applied Geophysics》;20170331;第14卷(第01期);第65-66、69页 * |
| Modeling and imaging with the scalar generalized-screen algorithms in isotropic media;Jerome H. Le Rousseau et al.;《GEOPHYSICS》;20011031;第66卷(第5期);第1551-1568页 * |
| 稳定的保幅高阶广义屏地震偏移成像方法研究;刘定进等;《地球物理学报》;20120731;第55卷(第7期);第2402-2411页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111812708A (zh) | 2020-10-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111639429A (zh) | 基于切比雪夫多项式谱的水下声场数值模拟方法、系统及介质 | |
| CN112698390B (zh) | 叠前地震反演方法及装置 | |
| CN111766628A (zh) | 一种预条件的时间域弹性介质多参数全波形反演方法 | |
| CN116774292B (zh) | 一种地震波走时确定方法、系统、电子设备及存储介质 | |
| CN111597753B (zh) | 数据深度变化特征自适应的二维电阻率反演方法及系统 | |
| CN109633749A (zh) | 基于散射积分法的非线性菲涅尔体地震走时层析成像方法 | |
| CN112231974A (zh) | 基于深度学习的tbm破岩震源地震波场特征恢复方法及系统 | |
| CN108828659B (zh) | 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置 | |
| CN109670195B (zh) | 融合单井敏感性局域化的EnKF油藏辅助历史拟合方法 | |
| CN111812708B (zh) | 地震波成像方法及装置 | |
| bin Waheed et al. | A holistic approach to computing first-arrival traveltimes using neural networks | |
| CN109490954B (zh) | 波场正演模拟方法及装置 | |
| CN113866827B (zh) | 一种解释性速度建模地震成像方法、系统、介质和设备 | |
| CN116577821A (zh) | 基于神经网络的时移地震预测方法和系统 | |
| CN105698799A (zh) | 一种提高捷联惯导系统姿态精度的预处理最优fir滤波器 | |
| CN112731520A (zh) | 一种基于结构张量扩散滤波的全波形反演方法及系统 | |
| CN111624653B (zh) | 地震速度自动拾取方法及装置 | |
| CN113219534B (zh) | 一种叠前深度偏移速度质控方法、装置、介质及电子设备 | |
| CN113219533B (zh) | 一种叠前时间偏移速度建模方法、装置、介质及电子设备 | |
| US20220171083A1 (en) | Systems and methods to correct seismic data for undesired effects caused by near surface velocity variations by correcting for long-wavelength statics | |
| CN114002735B (zh) | 弯曲射线叠前时间偏移速度求取方法及装置 | |
| CN117607951A (zh) | 一种适用于三维隧道模型的多尺度全波形反演方法及系统 | |
| CN112014875B (zh) | 叠前地震反演方法及装置 | |
| CN111596346B (zh) | 弹性波速度反演方法和装置 | |
| Plotnitskii et al. | DL-fused elastic FWI: Application to marine streamer data |
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 |