CN107807392A - 一种自适应抗频散的分块时空双变逆时偏移方法 - Google Patents
一种自适应抗频散的分块时空双变逆时偏移方法 Download PDFInfo
- Publication number
- CN107807392A CN107807392A CN201710898194.5A CN201710898194A CN107807392A CN 107807392 A CN107807392 A CN 107807392A CN 201710898194 A CN201710898194 A CN 201710898194A CN 107807392 A CN107807392 A CN 107807392A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- time
- wave field
- double
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 230000005012 migration Effects 0.000 title claims abstract description 26
- 238000013508 migration Methods 0.000 title claims abstract description 26
- 230000003044 adaptive effect Effects 0.000 title claims abstract description 13
- 239000006185 dispersion Substances 0.000 title claims abstract description 13
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 32
- 230000008859 change Effects 0.000 claims abstract description 19
- 238000012545 processing Methods 0.000 claims abstract description 19
- 238000001514 detection method Methods 0.000 claims abstract description 12
- 238000001914 filtration Methods 0.000 claims abstract description 12
- 238000013316 zoning Methods 0.000 claims abstract description 11
- 238000003384 imaging method Methods 0.000 claims abstract description 9
- 230000008569 process Effects 0.000 claims abstract description 7
- 230000001419 dependent effect Effects 0.000 claims abstract description 4
- 238000005070 sampling Methods 0.000 claims description 17
- 238000009795 derivation Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 235000013399 edible fruits Nutrition 0.000 claims 1
- 238000010586 diagram Methods 0.000 description 8
- 238000004088 simulation Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 3
- 150000003839 salts Chemical class 0.000 description 3
- 230000005540 biological transmission Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 238000001615 p wave Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000010010 raising Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/20—Trace signal pre-filtering to select, remove or transform specific events or signal components, i.e. trace-in/trace-out
- G01V2210/27—Other pre-filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/324—Filtering
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
本发明涉及一种自适应抗频散的分块时空双变逆时偏移方法，包括以下步骤：1)对观测系统覆盖的计算区域做变网格划分后得到空间网格；2)采用变网格方法在对空间网格做可变处理的基础上，并考虑时间步长的变化；3)采用Lanczos滤波方法处理在变网格过渡区域中边点的波场值；4)对观测系统覆盖的计算区域利用分块时空双变算法，根据逆时偏移的流程进行震源波场的正向外推并存储震源波场；5)对该计算区域利用相同的分块时空双变算法，进行检波波场的反向外推并存储检波波场；6)应用互相关成像条件，对外推后的震源波场和检波波场进行相关成像，并对成像结果进行常规叠加和拉普拉斯滤波等叠后处理得到最终的逆时偏移结果。
Description
技术领域
本发明涉及一种逆时偏移方法，具体是关于一种自适应抗频散的分块时空双变逆时偏移方法。
背景技术
一般的有限差分地震逆时偏移方法基于笛卡尔坐标系中的规则网格，是处理非均匀介质的有效方法。但是，随着社会工业对油气资源的需求不断增长，地震勘探面临的地质条件日趋复杂，包括强纵横向变速区域、低速带、复杂构造区域及碳酸盐岩储层的小型孔缝洞等。当处理这些区域的地震资料时，为了保证计算精度和稳定性，网格间距必须取得很小，而这将导致波场外推计算储存量的增加和计算量的增大。从空间采样的角度考虑，最有效的提高模拟精度同时又降低计算机内存需求的方法，就是模拟的不同区域采用不同的网格步长，即变网格。
但上述可变网格方法仅仅是在空间网格上做可变处理，没有考虑时间上的变化。考虑稳定性条件的限制，时间步长由最小网格间距确定，因此在大网格区域采用小时间步长就会造成时间过采样，从而限制效率的进一步提升；更为重要的是，有研究表明：在大网格区域采用小时间步长，不仅不会提高精度，还可能会带来频散误差。前人的研究实现了基于交错网格的双变算法，并做了详细的误差分析，但他们都没有解决变网格方法在大采样时间下的稳定性问题，并且人为误差也较大。
由于实际地下常存在低速体、低降速带、碳酸盐岩深部构造等，所以采样记录时间较大，并且，局部变时间算法也增大了时间上的迭代次数，从而可能引起大时间采样下的不稳定。为了提高变网格算法的稳定性，将Lanczos滤波算子推广到时空双变算法中，不仅提高了稳定性，而且减弱了由于变网格引起的数值干扰。
前人对变网格算法做了深入研究，但他们的变网格区域都仅是一个区域。而实际地下介质异常复杂，常存在多个需网格加密的区域，并且这些区域可能相距较远。例如，可能同时存在近地表低降速带和地下深部低速体，此时如果仍然采用统一的变网格算法，则会造成对低速体之间高速区域的过采样，降低了模拟效率；并且，当两个低速体的速度有差异时，就需要不同的网格变化倍数，而采用统一变网格则需要根据速度更小的低速体来决定变网格倍数，这样就会造成对另一个低速体的过采样。
发明内容
针对上述问题，本发明的目的是提供一种兼顾计算精度和计算效率的自适应抗频散的分块时空双变逆时偏移方法。
为实现上述目的，本发明采取以下技术方案：一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，包括以下步骤：1)对用于采集地面地震的观测系统覆盖的计算区域做变网格划分后得到空间网格；2)采用变网格方法在对空间网格做可变处理的基础上，并考虑时间步长的变化；同时，利用波场滤波和三次样条平滑解决变网格方法在大采样时间下的稳定性问题；3)采用Lanczos滤波方法处理在变网格过渡区域中边点的波场值；4)对观测系统覆盖的计算区域利用分块时空双变算法，根据逆时偏移的流程进行震源波场的正向外推并存储震源波场；5)对观测系统覆盖的计算区域利用相同的分块时空双变算法，进行检波波场的反向外推并存储检波波场；6)应用互相关成像条件，对外推后的震源波场和检波波场进行相关成像，并对成像结果进行常规叠加和拉普拉斯滤波等叠后处理得到最终的逆时偏移结果。
在上述步骤2)中，采用式(1)中的各向同性非均匀介质的2D声波一阶速度—应力方程进行变网格下的波场计算：
其中，νx,νz分别表示横向和纵向的质点速度；p表示应力向量；ρ表示正演所用模型的密度；Vp表示正演所用模型的速度；t表示时间；x和z分别表示空间的横向和纵向坐标。
变网格方法的任意偶数阶精度差分近似式为：
其中，Dx表示x方向的空间导数；f(x,z)表示(x,z)坐标处的波场值；i表示序号；N表示计算阶数；Δ2i和Δ2i-1均表示x方向网格间隔的变化量；c2i和c2i-1均表示差分算子，由以下方程确定：
其中，n表示变时间倍数。
通过常规的有限差分算法计算得到大时间步长的波场值，然后利用大时间步长求出网格步长变化处的后一时刻的波场值，从而结合前一时刻的大时间步长的波场值，利用双线性插值公式(4)得到小时间步长各个时刻的波场值：
其中，f(i2)表示第i2个时刻的波场值，i2＝1,2,...,n-1；F0,F1分别表示大时间步长的前一时刻和后一时刻的波场值。
在上述步骤3)中，利用变网格过渡区域中边点周围若干个细网格点值计算该边点相应的粗网格点值：
其中，F(i3,j)表示粗网格点值；f(i3+m,j+l)表示细网格点值；i3表示横向网格点序号；j表示纵向网格点序号；m表示偏离第i3个网格点的横向网格数目；l表示偏离第j个网格点的纵向网格数目；k表示变网格倍数；ωml表示Lanczos滤波系数。
分块时空双变算法的具体实现步骤为：在每个时间步长，先更新加密区域以外的常规粗网格，再判断波场是否传递到第一精细区域：若传递到第一精细区域则利用时空双变原理进行精细处理；若未传递到，则采用粗网格更新，然后判断波场是否传递到第二精细区域：若传递到第二精细区域则利用时空双变原理进行精细处理；若未传递到，则采用粗网格更新；同理，依次类推对其他精细区域做相同处理。
本发明由于采取以上技术方案，其具有以下优点：1、本发明提出分块时空双变的思想，即可以同时存在多个变网格区域，且每个区域的变网格倍数可以自由选择，从而可以最大化的提高模拟效率。由于变网格区域之间是相互独立的，所以本发明方法易于推广到任意个时空双变区域。2、本发明处理在变网格过渡区域的波场传递问题时引入了Lanczos滤波方法，大大提高了变网格算法的稳定性，并将本发明方法拓展到时空双变的情况。3、本发明方法应用在处理含有多个变网格区域，得到了较好的效果，精度上可以达到和全局细网格一致，效率上相对全局细网格和统一变网格都有大幅度提升。
附图说明
图1是局部变时间示意图；
图2是复杂模型示意图；
图3(a)和(b)是低速体模型的速度场示意图；
图4(a)-(c)分别是全局细网格、统一变网格和分块变网格的炮记录的波形对比图；
图5是去掉直达波的波形对比图；
图6是全局细网格、统一变网格和分块变网格的计算耗时对比图；
图7(a)和(b)是含起伏分界面的低速体模型的速度场示意图；
图8(a)-(c)是分别采用全局粗网格、分块变网格和全局细网格的单炮记录示意图；
图9是不同网格的模拟数据波形对比图；
图10(a)-(c)是分块网格划分示意图；
图11是三维SEG/EAGE盐丘速度模型；
图12是三维SEG/EAGE窄方位角数据观测系统示意图；
图13(a)是三维显示方式的3D RTM偏移结果示意图；
图13(b)是二维显示方式的3D RTM偏移结果示意图。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。
本发明提供的自适应抗频散的分块时空双变逆时偏移方法，其包括以下步骤：
1)对用于采集地面地震的观测系统覆盖的计算区域做变网格划分后得到空间网格，网格尺寸可以根据偏移所用的速度场分布灵活确定。
2)采用变网格方法在对空间网格做可变处理的基础上，并考虑时间步长的变化；同时，利用波场滤波和三次样条平滑解决变网格方法在大采样时间下的稳定性问题。
采用式(1)中的各向同性非均匀介质的2D声波一阶速度—应力方程进行变网格下的波场计算：
其中，νx,νz分别表示横向和纵向的质点速度；p表示应力向量；ρ表示正演所用模型的密度；Vp表示正演所用模型的速度；t表示时间；x和z分别表示空间的横向和纵向坐标。
与传统交错网格不同，变网格方法的差分算子是空间变化的，网格间距△x不再为常数值，可记为Δix，它是网格步长的函数。变网格方法的任意偶数阶精度差分近似式为：
其中，Dx表示x方向的空间导数；f(x,z)表示(x,z)坐标处的波场值；i表示序号；N表示计算阶数；Δ2i和Δ2i-1均表示x方向网格间隔的变化量；c2i和c2i-1均表示差分算子，可由以下方程确定：
其中，n表示变时间倍数。
变时间算法的关键之处是小时间步长的更新计算，在每个精细时间层的边界部分都需要获得一个初始值。大时间步长的波场值易于通过常规的有限差分算法计算得到，然后利用大时间步长求出网格步长变化处的后一时刻的波场值，从而结合前一时刻的大时间步长的波场值，利用双线性插值公式(4)就可得到小时间步长各个时刻的波场值(如图1所示)：
其中，f(i2)表示第i2个时刻的波场值，i2＝1,2,...,n-1；F0,F1分别表示大时间步长的前一时刻和后一时刻的波场值。
3)变网格方法中需要粗细网格间的波场传递，此时若直接从细网格赋值给粗网格(直接传递法)，由于网格步长的突变，则可能产生较强的不稳定。具体原因是由于细网格允许的最小波长小于粗网格允许的最小波长，而一般情况下，细网格内传播的真实波长是介于两者之间的，当直接由细网格传递到粗网格时，因为真实波长小于粗网格允许的最小波长，从而可能会出现不稳定现象。因此，采用Lanczos滤波方法处理在变网格过渡区域中边点的波场值，即利用变网格过渡区域中边点周围若干个细网格点值计算该边点相应的粗网格点值：
其中，F(i3,j)表示粗网格点值；f(i3+m,j+l)表示细网格点值；i3表示横向网格点序号；j表示纵向网格点序号；m表示偏离第i3个网格点的横向网格数目；l表示偏离第j个网格点的纵向网格数目；k表示变网格倍数；ωml表示Lanczos滤波系数。这样，横向第i3个和纵向第j个网格点的波场值就可以通过横向第i3+m个和纵向第j+l个网格点的波场值计算实现。
4)对观测系统覆盖的计算区域利用分块时空双变算法，根据逆时偏移的流程进行震源波场的正向外推并存储震源波场。
当存在多个变网格区域时，每个变网格区域都是相互独立的、互不影响的，即一块区域需不需要加密处理与另一块区域没有关系。如图2所示，只有两个变网格区域A和B，则在每个时间步进，包括四种情况：A加密B不加密，A加密B加密，A不加密B不加密，A不加密B加密。考虑到A和B是两个独立事件，可以统一为：在每个时间步进，对A单独判断处理，对B单独判断处理，这样上述四种情况都包含了。
本发明提出分块时空双变算法，即可以同时存在多个变网格区域，且每个变网格区域的变网格倍数可以自由选择，从而可以最大化的提高模拟效率。分块时空双变算法的具体实现步骤为：在每个时间步长，先更新加密区域以外的常规粗网格，再判断波场是否传递到精细区域A：若传递到精细区域A则利用时空双变原理进行精细处理；若未传递到，则采用粗网格更新。然后判断波场是否传递到精细区域B：若传递到精细区域B则利用时空双变原理进行精细处理；若未传递到，则采用粗网格更新。同理，依次类推对其他精细区域(C、D...)做相同处理。
5)对观测系统覆盖的计算区域利用相同的分块时空双变算法，进行检波波场的反向外推并存储检波波场。
6)应用互相关成像条件，对外推后的震源波场和检波波场进行相关成像，并对成像结果进行常规叠加和拉普拉斯滤波等叠后处理得到最终的逆时偏移结果。
下面通过几个具体实施例来说明本发明的效果。
1、如图3-图6所示，是本发明应用到一个相距较远的低速体模型的试验。
如图3(a)所示，为内部含两个相距较远的低速体模型，上层P波速度为3000m/s，下层为3500m/s，左上角低速体的速度为1000m/s，右下角低速体速度为600m/s。如果使用常规算法，为了满足稳定性条件，网格间距和采样间隔都要取得很小，这样就造成了对其他区域的过采样，从而大大的降低了效率；而如果采用统一的变网格，由于两个低速体相距较远，就造成了对低速体之间区域的过采样问题，且都采用5倍变网格也会对左上角低速体过采样；最佳的处理方法是采用分块时空双变算法，即分别对两个低速体进行变网格处理，左上角的采用3倍变网格，右下角的采用5倍变网格(如图3(b))。为了对比分析，将上述三种方法(全局细网格、统一变网格、分块变网格)都进行了实现，其相应的单炮记录分别如图4(a)-(c)所示，从中可以看出三种方法的炮记录是相同的。为了更加细致的研究分析，选取几个单道波形比对图(如图图5所示，已去掉直达波)。从图5可以看出，三种方法得到的结果基本是相同的，仅有微小的差异，误差是可以忽略的。从而可以说明，分块时空双变算法是正确的、可靠的。
虽然三种方法的结果基本相同，但其存储量和计算量却有巨大差别，其计算耗时对比如图6所示，从中可以发现：统一变网格相对于全局细网格效率提高了4.9倍，分块细网格相对于统一变网格效率提高了12.5倍，分块变网格相对于全局细网格效率提高了61.2倍。因此采用分块变网格大大提高了效率，仅是统一变网格耗时的8％，全局细网格耗时的1.6％。
2)如图7-图9所示，是本发明应用到复杂近地表和深部构造模型的试验。
实际地质体中，近地表区域一般存在需要精细研究的起伏界面，并且深部也常含有低速体构造(如图7(a)所示)。若是采用统一变网格处理，由于二者相距较远，仍然会极大地浪费时间，因此分块变网格是一个较好的选择。
分别对近地表起伏界面和深部低速体做3倍双变处理，网格剖分如图7(b)所示。图8(a)-(c)为分别采用全局粗网格、分块变网格和全局细网格的单炮记录，可以看出，三者基本是相同的，只是全局粗网格得到的结果中含有较多的虚假绕射。同样取出多个单道记录进行对比分析(如图9所示)，可以发现，分块变网格和全局细网格吻合的非常好，都要比全局粗网格更加精确，这是由于粗网格对起伏界面的粗糙离散采样所致。
3)如图10所示，是本发明应用到更为复杂的模型的试验。
分块时空双变算法应用十分广泛，不仅仅局限于分布较远的多个地质目标体，在其他多种情况下也同样适用：(1)一个倾斜的大断裂或倾斜的切深较大的狭长裂缝(如图10(a)所示)。由于切深较大，精细研究时若采用统一变网格，则其加密如图10(b)所示，极大地降低了效率、增加了计算和存储。采用分块变网格，如图10(c)所示，将大断裂或大裂缝划分为5块，则加密区域仅是统一变网格的1/5，因此分块变网格的效率是统一变网格的5倍，即计算时间和内存相对于统一变网格节省了80％。并且，显然分块越多提高的效率也就越大。(2)模拟微小裂缝时。由于裂缝尺度较小，变网格倍数可以达到上千倍，此时两个裂缝即使相距很近，但上千倍变网格时仍然会极大地增加计算和存储量，而采用分块变网格，就可以仅仅针对每个裂缝的微小区域进行网格加密，不需要对两个裂缝之前的区域做加密处理，从而可以节省内存和计算耗时。
4)如图11所示，是本发明应用于三维声波逆时偏移的试验。
分块时空双变算法使得逆时偏移算法可以高效率地处理三维大数据，同时保证了高精度。如图11所示，为一个典型的三维复杂构造，模型中部存在一个巨大的盐丘，使得盐下的成像成为一个巨大的挑战。该模型在In-line(纵测线)，Cross-line(横测线)，深度方向采样分别为676，676，210，采样间隔均为20m。使用该模型的3C-NA(窄方位角)数据进行测试，其数据观测系统如图12所示。该数据大小约为6.25GB，共有4782炮，为1炮8线接收，接收线间距为40m，炮线间距为320m，时间采样为615，采样间隔为8ms。图13(a)和图13(b)为3DRTM成像结果，可以看到无论在In-line方向和Cross-line方向，RTM偏移都较好的恢复了该模型的正确构造形态。
综上所述，本发明提出了分块时空双变算法，相比于传统的双变算法极大地提高了效率、节省了内存和计算耗时。并且，在过渡区域波场传递时引入了Lanczos滤波算子，提高了变网格算法的稳定性、减弱了网格变化处的虚假反射。数值试验表明，本发明在处理含有多个变网格区域时，得到了较好的效果，精度上可以达到和全局细网格一致，效率上相对全局细网格和统一变网格都有大幅度提升。更进一步地，本发明还可以应用于倾斜大切深断裂和微小裂缝模拟中，从而对西部深层碳酸盐岩储层的高精度正演模拟及成像都有重要的现实意义。
上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。
Claims (6)
1.一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，包括以下步骤：
1)对用于采集地面地震的观测系统覆盖的计算区域做变网格划分后得到空间网格；
2)采用变网格方法在对空间网格做可变处理的基础上，并考虑时间步长的变化；同时，利用波场滤波和三次样条平滑解决变网格方法在大采样时间下的稳定性问题；
3)采用Lanczos滤波方法处理在变网格过渡区域中边点的波场值；
4)对观测系统覆盖的计算区域利用分块时空双变算法，根据逆时偏移的流程进行震源波场的正向外推并存储震源波场；
5)对观测系统覆盖的计算区域利用相同的分块时空双变算法，进行检波波场的反向外推并存储检波波场；
6)应用互相关成像条件，对外推后的震源波场和检波波场进行相关成像，并对成像结果进行常规叠加和拉普拉斯滤波等叠后处理得到最终的逆时偏移结果。
2.如权利要求1所述的一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，在上述步骤2)中，采用式(1)中的各向同性非均匀介质的2D声波一阶速度—应力方程进行变网格下的波场计算：
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>&rho;</mi>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>v</mi>
<mi>x</mi>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>p</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>x</mi>
</mrow>
</mfrac>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&rho;</mi>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>v</mi>
<mi>z</mi>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>p</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>z</mi>
</mrow>
</mfrac>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>p</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>t</mi>
</mrow>
</mfrac>
<mo>=</mo>
<msup>
<msub>
<mi>&rho;V</mi>
<mi>p</mi>
</msub>
<mn>2</mn>
</msup>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>v</mi>
<mi>x</mi>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>x</mi>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>v</mi>
<mi>z</mi>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>x</mi>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，νx,νz分别表示横向和纵向的质点速度；p表示应力向量；ρ表示正演所用模型的密度；Vp表示正演所用模型的速度；t表示时间；x和z分别表示空间的横向和纵向坐标。
3.如权利要求2所述的一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，变网格方法的任意偶数阶精度差分近似式为：
<mrow>
<msub>
<mi>D</mi>
<mi>x</mi>
</msub>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mo>&lsqb;</mo>
<msub>
<mi>c</mi>
<mrow>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>+</mo>
<msub>
<mi>&Delta;</mi>
<mrow>
<mn>2</mn>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>c</mi>
<mrow>
<mn>2</mn>
<mi>i</mi>
</mrow>
</msub>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>-</mo>
<msub>
<mi>&Delta;</mi>
<mrow>
<mn>2</mn>
<mi>i</mi>
</mrow>
</msub>
<mo>,</mo>
<mi>z</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，Dx表示x方向的空间导数；f(x,z)表示(x,z)坐标处的波场值；i表示序号；N表示计算阶数；Δ2i和Δ2i-1均表示x方向网格间隔的变化量；c2i和c2i-1均表示差分算子，由以下方程确定：
其中，n表示变时间倍数。
4.如权利要求3所述的一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，通过常规的有限差分算法计算得到大时间步长的波场值，然后利用大时间步长求出网格步长变化处的后一时刻的波场值，从而结合前一时刻的大时间步长的波场值，利用双线性插值公式(4)得到小时间步长各个时刻的波场值：
<mrow>
<mi>f</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>i</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>-</mo>
<msub>
<mi>i</mi>
<mn>2</mn>
</msub>
<mo>)</mo>
<msub>
<mi>F</mi>
<mn>0</mn>
</msub>
<mo>+</mo>
<msub>
<mi>i</mi>
<mn>2</mn>
</msub>
<msub>
<mi>F</mi>
<mn>1</mn>
</msub>
</mrow>
<mi>n</mi>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，f(i2)表示第i2个时刻的波场值，i2＝1,2,...,n-1；F0,F1分别表示大时间步长的前一时刻和后一时刻的波场值。
5.如权利要求4所述的一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，在上述步骤3)中，利用变网格过渡区域中边点周围若干个细网格点值计算该边点相应的粗网格点值：
<mrow>
<mi>F</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>i</mi>
<mn>3</mn>
</msub>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<mi>k</mi>
</mrow>
<mrow>
<mn>2</mn>
<mi>k</mi>
</mrow>
</munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>l</mi>
<mo>=</mo>
<mo>-</mo>
<mn>2</mn>
<mi>k</mi>
</mrow>
<mrow>
<mn>2</mn>
<mi>k</mi>
</mrow>
</munderover>
<msub>
<mi>&omega;</mi>
<mrow>
<mi>m</mi>
<mi>l</mi>
</mrow>
</msub>
<mi>f</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>i</mi>
<mn>3</mn>
</msub>
<mo>+</mo>
<mi>m</mi>
<mo>,</mo>
<mi>j</mi>
<mo>+</mo>
<mi>l</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
其中，F(i3,j)表示粗网格点值；f(i3+m,j+l)表示细网格点值；i3表示横向网格点序号；j表示纵向网格点序号；m表示偏离第i3个网格点的横向网格数目；l表示偏离第j个网格点的纵向网格数目；k表示变网格倍数；ωml表示Lanczos滤波系数。
6.如权利要求5所述的一种自适应抗频散的分块时空双变逆时偏移方法，其特征在于，分块时空双变算法的具体实现步骤为：在每个时间步长，先更新加密区域以外的常规粗网格，再判断波场是否传递到第一精细区域：若传递到第一精细区域则利用时空双变原理进行精细处理；若未传递到，则采用粗网格更新，然后判断波场是否传递到第二精细区域：若传递到第二精细区域则利用时空双变原理进行精细处理；若未传递到，则采用粗网格更新；同理，依次类推对其他精细区域做相同处理。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710898194.5A CN107807392A (zh) | 2017-09-28 | 2017-09-28 | 一种自适应抗频散的分块时空双变逆时偏移方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710898194.5A CN107807392A (zh) | 2017-09-28 | 2017-09-28 | 一种自适应抗频散的分块时空双变逆时偏移方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107807392A true CN107807392A (zh) | 2018-03-16 |
Family
ID=61591865
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710898194.5A Pending CN107807392A (zh) | 2017-09-28 | 2017-09-28 | 一种自适应抗频散的分块时空双变逆时偏移方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107807392A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112230277A (zh) * | 2020-09-30 | 2021-01-15 | 山东大学 | 基于柱坐标系的隧道地震波传播数值模拟方法及系统 |
| CN115081267A (zh) * | 2022-05-18 | 2022-09-20 | 内蒙古农业大学 | 基于优化的时空变步长有限差分地震波数值模拟方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102183790A (zh) * | 2011-02-12 | 2011-09-14 | 中国石油大学(华东) | 基于时空双变网格的弹性波正演模拟技术 |
| CN103926619A (zh) * | 2014-05-06 | 2014-07-16 | 王维红 | 一种三维vsp数据的逆时偏移方法 |
| CN104360381A (zh) * | 2014-10-20 | 2015-02-18 | 李闯 | 一种地震资料的偏移成像处理方法 |
| CN105388520A (zh) * | 2015-10-22 | 2016-03-09 | 中国石油化工股份有限公司 | 一种地震资料叠前逆时偏移成像方法 |
| CN105652320A (zh) * | 2015-12-30 | 2016-06-08 | 中国石油天然气集团公司 | 逆时偏移成像方法和装置 |
-
2017
- 2017-09-28 CN CN201710898194.5A patent/CN107807392A/zh active Pending
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102183790A (zh) * | 2011-02-12 | 2011-09-14 | 中国石油大学(华东) | 基于时空双变网格的弹性波正演模拟技术 |
| CN103926619A (zh) * | 2014-05-06 | 2014-07-16 | 王维红 | 一种三维vsp数据的逆时偏移方法 |
| CN104360381A (zh) * | 2014-10-20 | 2015-02-18 | 李闯 | 一种地震资料的偏移成像处理方法 |
| CN105388520A (zh) * | 2015-10-22 | 2016-03-09 | 中国石油化工股份有限公司 | 一种地震资料叠前逆时偏移成像方法 |
| CN105652320A (zh) * | 2015-12-30 | 2016-06-08 | 中国石油天然气集团公司 | 逆时偏移成像方法和装置 |
Non-Patent Citations (2)
| Title |
|---|
| 张慧 等: "基于双变网格算法的地震波正演模拟", 《地球物理学报》 * |
| 李振春 等: "一种稳定的高精度双变网格正演模拟与逆时偏移方法", 《石油物探》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112230277A (zh) * | 2020-09-30 | 2021-01-15 | 山东大学 | 基于柱坐标系的隧道地震波传播数值模拟方法及系统 |
| CN115081267A (zh) * | 2022-05-18 | 2022-09-20 | 内蒙古农业大学 | 基于优化的时空变步长有限差分地震波数值模拟方法 |
| CN115081267B (zh) * | 2022-05-18 | 2023-08-29 | 内蒙古农业大学 | 基于优化的时空变步长有限差分地震波数值模拟方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US5999488A (en) | Method and apparatus for migration by finite differences | |
| EP2869096B1 (en) | Systems and methods of multi-scale meshing for geologic time modeling | |
| CN102183790A (zh) | 基于时空双变网格的弹性波正演模拟技术 | |
| CN105319589B (zh) | 一种利用局部同相轴斜率的全自动立体层析反演方法 | |
| CN103926619B (zh) | 一种三维vsp数据的逆时偏移方法 | |
| CN108445538B (zh) | 基于反射地震资料建立深度域层q模型的方法和系统 | |
| CN105277978A (zh) | 一种确定近地表速度模型的方法及装置 | |
| CN107462924B (zh) | 一种不依赖于测井资料的绝对波阻抗反演方法 | |
| CN108508482A (zh) | 一种地下裂缝地震散射响应特征模拟方法 | |
| CN106353797A (zh) | 一种高精度地震正演模拟方法 | |
| CN107479092A (zh) | 一种基于方向导数的频率域高阶声波方程正演模拟方法 | |
| CN105911584B (zh) | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 | |
| CN106646645A (zh) | 一种新的重力正演加速方法 | |
| CN104237937A (zh) | 叠前地震反演方法及其系统 | |
| CN104975851B (zh) | 用于振幅随炮检距变化道集分析的油藏模型优化方法 | |
| CN110031898A (zh) | 数据优化方法及积分法叠前深度偏移方法 | |
| CN110389382A (zh) | 一种基于卷积神经网络的油气藏储层表征方法 | |
| CN110286410A (zh) | 基于绕射波能量的裂缝反演方法和装置 | |
| CN102830431B (zh) | 真地表射线追踪自适应插值方法 | |
| CN110261903B (zh) | 一种基于逆时能量聚焦的地下震源被动定位方法 | |
| CN107807392A (zh) | 一种自适应抗频散的分块时空双变逆时偏移方法 | |
| CN105182414B (zh) | 一种基于波动方程正演去除直达波的方法 | |
| CN106662665B (zh) | 用于更快速的交错网格处理的重新排序的插值和卷积 | |
| CN111948708B (zh) | 一种浸入边界起伏地表地震波场正演模拟方法 | |
| CN104597496B (zh) | 一种二维地震资料中速度的三维空间归位方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| CB02 | Change of applicant information | ||
| CB02 | Change of applicant information |
Address after: 100010 Beijing, Chaoyangmen, North Street, No. 25, No. Applicant after: China Offshore Oil Group Co., Ltd. Applicant after: CNOOC research institute limited liability company Address before: 100010 Beijing, Chaoyangmen, North Street, No. 25, No. Applicant before: China National Offshore Oil Corporation Applicant before: CNOOC Research Institute |
|
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20180316 |