CN107153216A - 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 - Google Patents
确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 Download PDFInfo
- Publication number
- CN107153216A CN107153216A CN201710542740.1A CN201710542740A CN107153216A CN 107153216 A CN107153216 A CN 107153216A CN 201710542740 A CN201710542740 A CN 201710542740A CN 107153216 A CN107153216 A CN 107153216A
- Authority
- CN
- China
- Prior art keywords
- wave
- wave field
- field
- shear
- zero
- 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
- 239000013598 vector Substances 0.000 title claims abstract description 125
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000003860 storage Methods 0.000 title claims abstract description 25
- 238000000926 separation method Methods 0.000 claims abstract description 67
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 53
- 230000009466 transformation Effects 0.000 claims description 29
- 238000000354 decomposition reaction Methods 0.000 claims description 28
- 238000004590 computer program Methods 0.000 claims description 15
- 238000012545 processing Methods 0.000 claims description 11
- 238000005457 optimization Methods 0.000 claims description 5
- 238000004088 simulation Methods 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 3
- 239000002131 composite material Substances 0.000 claims description 2
- 230000006872 improvement Effects 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 12
- 238000003384 imaging method Methods 0.000 abstract description 8
- 230000006870 function Effects 0.000 description 91
- 238000010586 diagram Methods 0.000 description 17
- 238000004364 calculation method Methods 0.000 description 15
- 230000008878 coupling Effects 0.000 description 13
- 238000010168 coupling process Methods 0.000 description 13
- 238000005859 coupling reaction Methods 0.000 description 13
- 230000010287 polarization Effects 0.000 description 13
- 230000008859 change Effects 0.000 description 10
- 238000013459 approach Methods 0.000 description 8
- 230000000694 effects Effects 0.000 description 8
- 238000009826 distribution Methods 0.000 description 6
- 230000000644 propagated effect Effects 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 230000004044 response Effects 0.000 description 5
- 230000003287 optical effect Effects 0.000 description 4
- 230000000903 blocking effect Effects 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000005096 rolling process Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 239000003989 dielectric material Substances 0.000 description 1
- 238000007598 dipping method Methods 0.000 description 1
- 230000005684 electric field Effects 0.000 description 1
- 230000005672 electromagnetic field Effects 0.000 description 1
- 238000002955 isolation Methods 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000000844 transformation 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
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
本申请公开了一种确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质。该方法包括：构建地下炮点正传波场及检波点反传波场；判断传播正传波场及反传波场的介质类型；当介质为各向同性介质时，根据预设方程解耦算法，对正传波场和或反传波场进行纵横波场分离，当介质为非均匀各向异性时，根据预设变异函数IDW插值算法，对正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到多个方向波波场；从多个方向波波场中选择出多个方向波波场，并计算选择出的多个方向波波场的坡印廷矢量。本申请的方法能计算分解后的各方向波的坡印廷矢量，进而获取更加准确的角度域成像道集，提高成像质量。
Description
技术领域
本申请涉及地震勘探领域，尤其是涉及一种确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质。
背景技术
坡印廷(即Poynting)矢量是能流密度矢量。其中，坡印廷矢量用来表示单位时间内穿过垂直于此矢量方向的单位表面的电磁场能量。其矢量值为电场强度与磁场强度的叉乘，方位代表电磁波的传播方向。
Yoon和Marfurt(2003)将坡印廷矢量引入地震勘探领域，给出了地震波场中的坡印廷矢量的计算公式，并利用炮点波场的坡印廷矢量和检波点波场的坡印廷矢量输出波场的反射角，继而生成ADCIGs道集。其中，求取的ADCIGs道集可以对RTM算法进行约束。
Yoon和Marfurt(2006)提出使用坡印廷矢量法对互相关成像时的成像角度进行限制，本质上是消除特定角度的成像，对成像结果进行一个角度域的滤波，从而消除低频噪音。
现有技术中，坡印廷矢量代表整个波场中的能流密度的传播方向，求取的传播方向也是整个波场的传播方向。对于复杂介质或者各向异性介质中弹性波的相关计算，由于各种波型耦合，波前重叠，求取的整个波场的坡印廷矢量中的传播方向并不是正确的地震波的传播方向，因而影响ADCIGs道集的生成，从而影响使用ADCIGs道集进行成像的质量。
发明内容
本发明的目的在于提供了一种确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质，以解决上述的问题。
为了实现上述目的，本发明提供了一种确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质，该确定地震波场的坡印廷矢量的方法包括：
构建地下炮点正传波场以及检波点反传波场；
判断传播所述正传波场以及所述反传波场的介质类型；
当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；
对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
在一可选的实施方式中，所述对分离出的纵波波场和或横波波场在上下左右四个方向进行分解包括：
对分离出的纵波波场和横波波场进行上下方向的方向波分离，得到对应的上行方向波波场和下行方向波波场；
对上行方向波波场和下行方向波波场进行左右方向的方向波分离，得到对应的左上行、左下行方向波波场，和右上行、右下行方向波波场。
在一可选的实施方式中，所述对分离出的纵波波场进行上下方向的方向波分解，包括：
将分离出的纵波波场进行时间复数域拓展；
在Z方向对时间复数域拓展之后的纵波波场进行傅里叶变换，得到波数域纵波波场；
从所述波数域纵波波场中选择波数大于零的部分，波数小于零的部分置为零，进行Z方向上的反傅里叶变换，得到纵波波场的下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行Z方向上的反傅里叶变换，得到所述纵波波场的上行方向波波场。
在一可选的实施方式中，所述对上行方向波波场进行左右方向的方向波分解，包括：
将所述纵波波场的上行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述纵波波场的上行方向波波场进行傅里叶变换，得到波数域纵波波场；
从所述波数域纵波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到所述纵波波场的右上行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述纵波波场的左上行方向波波场。
在一可选的实施方式中，所述对下行方向波波场进行左右方向的方向波分解，包括：
将所述纵波波场的下行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述纵波波场的下行方向波波场进行傅里叶变换，得到波数域纵波波场；
从所述波数域纵波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向下的反傅里叶变换，得到所述纵波波场的右下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向下的反傅里叶变换，得到所述纵波波场的左下行方向波波场。
在一可选的实施方式中，所述对分离出的横波波场进行上下方向的方向波分解，包括：
将分离出的横波波场进行时间复数域拓展；
在Z方向对时间复数域拓展之后的横波波场进行傅里叶变换，得到波数域横波波场；
从所述波数域横波波场中选择波数大于零的部分，波数小于零的部分置为零，进行Z方向上的反傅里叶变换，得到横波波场的下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行Z方向上的反傅里叶变换，得到所述横波波场的上行方向波波场。
在一可选的实施方式中，所述对上行方向波波场进行左右方向的方向波分解，包括：
将所述横波波场的上行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述横波波场的上行方向波波场进行傅里叶变换，得到波数域横波波场；
从所述波数域横波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到所述横波波场的右上行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述横波波场的左上行方向波波场。
在一可选的实施方式中，所述对下行方向波波场进行左右方向的方向波分解，包括：
将所述横波波场的下行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述横波波场的下行方向波波场进行傅里叶变换，得到波数域横波波场；
从所述波数域横波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向下的反傅里叶变换，得到所述横波波场的右下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向下的反傅里叶变换，得到所述横波波场的左下行方向波波场。
在一可选的实施方式中，所述根据预设的变异函数IDW插值算法，对在所述目标地层的非均匀各向异性介质中传播的正传波场和或反传波场进行纵横波场分离，包括：
当在非均匀各向异性介质下时，利用傅里叶变换分别将各时刻的所述正传波场和或反传波场由空间域变换至波数域；
选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断优化；
在波数域利用所述分离算子对所述参考模型下的所述正传波场和或反传波场进行分离，得到所述参考模型下的分离后的所述纵横波波场；
将所述参考模型下的分离后的纵横波波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的纵横波波场，N为正整数；
在空间域，利用变异函数改进IDW插值算法并根据改进后的所述IDW插值算法计算N个参考模型的权重系数，对N个所述参考模型下的分离后的纵横波波场进行加权插值处理，得到最终分离后的纵波波场和横波波场。
在一可选的实施方式中，所述选取N个参考模型，包括：
遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到所述参考模型。
在一可选的实施方式中，所述参考模型的权重系数通过所述变异函数计算得到。
在一可选的实施方式中，所述变异函数的获取方法包括：
通过计算所述初始模型的各向异性参数分布的变异函数值，拟合得到所述初始模型的变异函数。
在一可选的实施方式中，所述矢量成像为弹性波逆时偏移成像。
本申请还提供了一种确定地震波场的坡印廷矢量的装置，包括：
构建单元，用于构建地下炮点正传波场以及检波点反传波场；
分离单元，用于判断传播所述正传波场以及所述反传波场的介质类型；当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
计算单元，从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
本申请还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
构建地下炮点正传波场以及检波点反传波场；
判断传播所述正传波场以及所述反传波场的介质类型；
当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；
对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
本申请的方法可以将耦合的弹性波进行纵横波的分离以及上下左右方向波的分解以得到多个方向波，然后在此基础上计算该多个方向波的坡印廷矢量；然后，利用其坡印廷矢量输出波场的反射角，继而生成ADCIGs道集，防止因不同的波型以及方向波的耦合和或叠加导致弹性波传播方向求取不准，影响角度道集的生成以及后期的速度分析处理。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法的流程图；
图2为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中采用变异函数IDW插值算法将获取到的所述弹性波分离为纵波和横波的流程图；
图3为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中采用变异函数IDW插值算法过程中如何得到自褶组合窗函数的流程图；
图4为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中将纵波分离为上行纵波和下行纵波的流程图；
图5为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中上行横波分离为左上行横波和右上行横波的流程图；
图6为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中下行横波分离为左下行横波和右下行横波的流程图；
图7为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中横波分离为上行横波和下行横波的流程图；
图8为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中上行横波分离为左上行横波和右上行横波的流程图；
图9为本申请的实施例提供的一种确定地震波场的坡印廷矢量的方法中下行横波分离为左下行横波和右下行横波的流程图；
图10(a)为本申请的实施例1中弹性波在x分量上的坡印廷矢量的示意图；
图10(b)为本申请的实施例1中分离出来的纵波在x分量上的坡印廷矢量的示意图；
图10(c)为本申请的实施例1中分离出来的横波在x分量上的坡印廷矢量的示意图；
图10(d)为本申请的实施例1中弹性波在z分量上的坡印廷矢量的示意图；
图10(e)为本申请的实施例1中分离出来的纵波在z分量上的坡印廷矢量的示意图；
图10(f)为本申请的实施例1中分离出来的横波在z分量上的坡印廷矢量的示意图；
图11(a)为本申请实施例2中的点脉冲波的坡印廷矢量；
图11(b)-11(c)为图11(a)中的点脉冲波在进行纵横波地分离后得到的坡印廷矢量；
图11(d)-11(g)为本申请实施例2中在将纵横波分离出的基础上再进行上下左右方向波的分离后计算得出的坡印廷矢量；
图12a为Marmousi2模型各向异性参数分布及选取的参考模型；
图12b为Marmousi2模型各向异性参数分布及选取的参考模型；
图12c为Marmousi2模型各向异性参数分布及选取的参考模型；
图12d为Marmousi2模型各向异性参数分布及选取的参考模型；
图13为本申请的实施例中变异函数IDW算法插值运算过程中选取的参考模型中的局部极值；
图14为本申请的实施例中变异函数IDW算法插值运算过程中选取的参考模型；
图15为本申请的实施例提供的坡印廷矢量计算装置的结构示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。
Wang和McMechan(2015)引入纵波应力，分别求取了各向同性介质中的纵波和横波的Poynting矢量，分别代表纵波和横波的传播方向。然而，本申请的发明人研究发现，其只是对纵横波进行解耦，以避免多波型的耦合对坡印廷矢量计算的影响，而没有考虑多个方向波重叠对对坡印廷矢量计算的影响。
此外，Wang等(2016)在求取坡印廷矢量矢量之前，分别进行纵横波的分离和上下方向波的分解，以避免多波前的混杂及多波型的耦合对坡印廷矢量计算的影响。然而，本申请的发明人研究发现，其纵横波分离方案是采用Zhang和McMechan(2010)的矢量分离算法，计算量巨大；且只是进行了上下两个方向的方向波分解，没有考虑左右方向波的分解；因此，对于存在垂直构造，其方案存在一定的局限，求取的坡印廷矢量在特定地质条件下可能存在不准确的情况。综合上述内容可知，不同的研究人员已经从不同的方面来试图解决坡印廷矢量计算的准确性，但各种方法的效果不佳或者无法在特定的条件下使用。
有鉴于此，为克服现有技术存在的缺陷，本申请的发明人提供了本申请的技术方案，具体请参见下文本申请实施例。
图1为本发明的实施例提供的确定地震波场的坡印廷矢量的方法的流程示意图，参照图1所示，可以包括以下步骤：
S101：构建地下炮点正传波场以及检波点反传波场；
S102：判断传播正传波场以及反传波场的介质类型；
S103：当介质为各向同性介质时，根据预设的方程解耦算法，对正传波场和或反传波场进行纵横波场分离；当介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对正传波场和或反传波场进行纵横波场分离；
S104：对分离出的纵波波场和横波波场在上下左右四个方向进行分解，得到所述纵波波场和/或所述横波波场的多个方向波波场；
S105：从所述多个方向波波场中选择出至少一方向波波场，并计算选择出的所述多个方向波波场的坡印廷矢量。
本申请的方法可以将弹性波进行纵横波的分离以及上下左右方向波分解，然后计算所有方向波波场的坡印廷矢量；然后，利用坡印廷矢量输出波场的反射角，继而生成对应的ADCIGs道集，此方法避免了波前重叠，波型耦合对波传播方向求取带来的影响，生成的ADCIGs道集更加准确，成像效果更好。
IDW(Inverse Distance Weighted)是一种常用而简便的空间插值方法。它以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样本点赋予的权重越大。例如，设平面上分布一系列离散点，已知其坐标和值为Xi，Yi,Zi(i＝1，2，…，n)通过距离加权值求z点值。
具体的，以Marmousi2模型为例，图12a到图12d以及图13为Marmousi2模型各向异性参数分布及选取的参考模型示意图。从上图可知，Marmousi2模型的各向异性参数的空间分布是不均匀的。在实际的计算过程中，各向同性介质不需要考虑两点之间的距离，因为对于各向同性介质，每个点都是一样的。但对于各向异性介质，每个点都是不同的。两点之间具有一定的距离，该距离代表的是一种相似性，假设越近的点相似性越大。因此，我们需要插值一个点，可以对它周围较近的点进行加权插值得到，越近的点权重越大。该方法是反距离插值(IDW)算法的一个应用，但该算法没有考虑空间变异性。因此，我们的方法还考虑空间变异性，并基于该变异性采用变异函数的IDW插值算法，进而能在纵横波分离时提高分离质量，从而提高成像质量。
其中，变异函数是两个点之间的距离函数，不仅描述了空间距离，也描述了两点的空间变异性，采用变异函数的插值计算可以更加地准确。变异函数的具体求解过程将会在下文阐述。
参照图2所示，S103步骤还包括以下几个子步骤：
S201：当在非均匀各向异性介质下时，利用傅里叶变换分别将各时刻的所述弹性矢量波场波场由空间域变换至波数域；
S202：选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断优化；
S203：在波数域利用所述分离算子对所述参考模型下的所述弹性波场进行分离，得到所述参考模型下的分离后的所述纵横波波场；
S204：将所述参考模型下的分离后的纵横波波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的纵横波波场，N为正整数；
S205：在空间域，利用变异函数插值理论，计算N个参考模型的权重系数，对N个所述参考模型下的分离后的纵横波波场进行加权插值处理，得到最终分离后的纵波波场和横波波场。
本发明实施例中，在波数域利用分离算子分离波场，在空间域运用基于变异函数的IDW算法插值运算得到最终分离后的波场，可以有效降低波场分离的计算复杂度。本发明实施例中，在波数域采用自褶积组合窗函数截断分离算子，再分离波场，可以提高波场分离的效率，可以保证截断得到的分离算子的准确度更高。
具体的，可以通过遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值R进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到上述参考模型(参看图14和图15所示)。
图3是本发明一实施例中获取自褶积组合窗函数的方法的流程示意图。如图3所示，可包括以下步骤：
S301：选择主瓣和旁瓣性能高于设定阈值的窗函数作为原始窗函数；
S302：对原始窗函数做L次自褶积运算得到自褶积后的窗函数，其中，L为正整数；
S303：对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积组合窗函数。
在上述实施例中，对选定的原始窗函数进行自褶积运算，自褶积运算使得原始窗函数的主瓣性能变差，旁瓣性能变好，然后根据当前的需要对原始窗函数和自褶积后的窗函数进行加权运算，从而得到一个满足要求的窗函数。利用该种自褶积组合窗函数的到的分离算子可以提高波场分离的精度和稳定性。
在上述步骤S302中，进行加权运算的加权系数可以具备特点：在确定需要主瓣性能优先的情况下，原始窗函数的加权系数高于自褶积后的窗函数的加权系数；在确定需要旁瓣性能优先的情况下，自褶积后的函数的加权系数高于原始窗函数的加权系数。以此在维持适当的主瓣宽度的前提下，可以尽可能增大旁瓣衰减。
为使本领域的技术人员更容易理解本发明的技术方案，下面将以一具体实施例说明变异函数IDW插值算法以及该步骤的子步骤。
对于非均匀各向异性而言，理论上每个采样点都需要根据相应的各向异性参数计算偏振向量，才能获得更好的分离结果，但是这样处理计算量太大，因此，可以从初始模型中选择出N个参考模型，分别在波数域进行矢量波场(弹性波)分离，然后反变换回空间域，根据事先计算好的权重系数，重构初始模型的分离结果。
因此，考虑到引入描述区域化变量的空间结构性变化和随机性变化的变异函数，代替两点之间的距离计算权重系数。具体的，可以按照以下公式计算初始模型各向异性参数分布的变异函数值：
其中，μ(h)表示变异函数值，Φ＝f[δ+2(ε-δ)sin2(α-θ)]sin 2(α-θ)，ε、δ表示TI介质各向异性的系数，θ表示TTI介质对称轴的倾角，表示插值参考点的位置增量，i＝1,2...n，其中，n表示插值参考点的个数。
在计算得到多个分散的变异函数值后，根据这多个变异函数值进行拟合，得到理论变异函数模型。
然后，按照以下公式计算各个参考模型的权重系数：
其中，wk表示权重系数，(εk,δk,θk)表示插值参考点。
上述表示位置增量，对于每一个都可以计算出一个以不同的的模为横坐标，为纵坐标，绘制点图，由绘制的点图可以发现：空间距离越小，变异性越小，越大，变异性越大，当增大到一个临界值R时，两个数值相互之间基本不会有太大的影响了，这也就表明在选取参考点，或者需要拟合的点的时候，应该在R范围内选取，否则选取的点没有太多的实际意义。例如：可以根据初始的参考模型计算得到多个分散的变异函数值，对得到的多个分散的变异函数值进行拟合得到理论变异函数模型，再根据拟合得到的变异函数模型确定变异函数中各个参数的参数值，遍历初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。
下面结合一个具体的实施例进行说明，然而值得注意的是，该具体实施例仅是为了更好地说明本发明，并不构成对本发明的不当限定。
基于偏振特性对矢量波场进行分离的过程中所存在的一个不可避免的问题是如何在保证计算效率的前提下，让波场更准确地投影到相应的偏振方向上。在波数域，称之为投影，而在空间域，称之为空间滤波。
在各向同性介质中，矢量波场的分离效果仅与微分算子的精度有关，精度越高，分离效果越好，同时计算量也增大。然而，在各向异性介质中，矢量波场的分离效果主要与各向异性参数和拟微分算子的精度这两个因素有关，考虑到在TI介质中，地震波动力学特征可以表示为各向同性部分与各向异性部分的叠加，可以通过Christoffel方程计算得到偏振向量，然后将其分解为波数向量与偏转向量之和。
其中，偏转向量是以波数为自变量的函数，具有空间分布特性，在波数域利用常规二项式窗函数和Gauss窗函数截断逼近拟微分算子，并在空间域利用IDW算法去插值偏振向量的各向异性部分，以降低波场分离的计算量，从而提高波场分离的效率与精度。
具体而言，影响基于偏振特性分离矢量波场的技术的因素主要有两个：其一是截断窗函数的幅频响应特性，窗函数的幅频响应的主瓣形状控制着过渡带的范围，也就是频谱覆盖范围，旁瓣的形状决定了差分算子逼近微分算子的偏差程度，主瓣和旁瓣性能直接影响到了逼近的精度；其二是插值算法的插值效果，对于一种插值算法，既要求其有较高的插值精度，又需要较少的计算量，IDW算法执行效率较高，但是仅考虑到两点之间的距离计算权重，因此插值精度并不是很高，因此需要需求一种在保证执行效率的情况下，又可以提高插值精度的插值算法。
在本例中，提出了一种矢量波场分离方法，主要是为了解决一下两个问题：一是截断窗函数的选择，二是插值算法的选择。
首先，因为窗函数幅度响应的主瓣和旁瓣性能直接影响着差分逼近拟微分算子的精度问题，如果想设计合适的窗函数，首先要明白主瓣和旁瓣的性能如何影响逼近精度，其次是要研究如何去控制主瓣和旁瓣。然后，是插值算法的选择，将空间插值算法的优化引入矢量波场分离技术，因为对于非均匀各向异性介质，如果每一个点都计算一次拟微分算子，分离矢量波场有更好的效果，但是，这将耗费巨大的计算量，假设拟微分算子的大小是n2，模型的大小是N2，那么对模型进行矢量波场分离的计算量是2n2N2，这远远大于m阶精度的有限差分法2mN2的计算量，n<m。因此需要利用一种插值算法，在空间域通过插值的形式重构初始模型的分离结果。
在本例中，所依据的基本原理是：在各向同性介质中，应用Helmholtz定理，分别对波场求取散度和旋度，以分离纵横波，有如下公式：
P＝▽·U,
S＝▽×U。
在波数域可以表示为：
P＝▽·U＝ikxUx+ikzUz；
S＝▽×U＝ikzUx-ikxUz。
由上述公式可以看出，在各向同性介质中，P波是矢量波场在波数方向的投影，S波是矢量波场在垂直波数方向的投影，P波的偏振向量为(kx,kz)，S波的偏振向量与其正交，因而为(kz,-kx)。
对于各向异性介质，通过各向异性介质相对应的Christoffel方程，也可以得到各向异性介质中纵横波的偏振向量，以二维TTI介质为例，TTI介质对应的Christoffel方程为：
其中：
其中，c11,c15,c33,c35,c55表示弹性系数张量，kx,kz表示归一化波数，Γ表示Christoffel矩阵，上述公式是典型的特征值和特征向量的问题，为使公式有非零解，就要使系数行列式为零，因此可以求得偏振向量P＝(Px,Pz)，这是kx,kz的函数，将iP反变换回空间域，便可以得到拟微分算子L，L也称作空间滤波算子。
现在引入表征TI介质(横向各向同性介质)各向异性的Thomsen系数VP0,VS0,ε,δγ,，其中，VP0,VS0,ε,δ与qP波和qSV波有关，VS0,γ与qSH波有关，因为qSH波是解耦的，因此矢量波场的分离仅需要VP0,VS0,ε,δ这四个参数。以二维TTI介质(具有倾斜对称轴的横向各向同性介质)为例，还需引入TI介质对称轴的倾角θ。在TI介质中，地震波动力学特征由两部分组成，第一个部分是各向同性的部分，第二个部分为各向异性的部分，可以将其近似表示为：
K≈Kiso+L(ε,δ)+Q(ε,δ,VS0) (4)
其中，K表示地震波的动力学特征，Kiso表示各向同性部分，ε＝δ＝0；L(ε,δ)+Q(ε,δ,VS0)表示各向异性部分，L(ε,δ)表示线性部分，Q(ε,δ,VS0)表示非线性部分。对于偏振向量，同样适用于上述近似公式，即，TI介质中的qP波和qSV波的偏振向量也由两部分组成：各向同性部分和各向异性部分。
对于各向同性部分，可以采用常规的有限差分算子的优化方法，即：使用窗函数截断方法或者最优化方法，这两种方法的本质都是类似的，都是希望在更高的波数范围达到一个较好的逼近精度。
对于各向异性部分，理论上每一个采样点都需要根据相应的各向异性参数计算偏振向量，以获得更好的分离效果，但这必然会带来海量的计算量，因此，考虑到采用插值的方式，即：从初始模型中选择N个参考模型，分别在波数域进行矢量波场分离，然后再反变换回空间域，再根据事先计算好的权重系数，重构初始模型的分离结果。既然是插值算法，插值的精度和插值算法的执行效率都是需要考虑的因素，IDW插值仅考虑到两点之间的距离计算权重，插值精度较差，且仅对均匀模型适用性较好，因此可以考虑引入描述区域化变量的空间结构性变化和随机性变化的变异函数，代替两点之间的距离，计算权重系数.
对上述内容具体阐述如下：
1、各向同性部分：
一个带限的连续信号f(x)在x＝0处的导数，可以被以一个均匀采样的信号fn表示为：
存在一个长度为N+1点的窗函数，N为偶数，去截断上述两个公式，并经过一些简单的处理，可以得到有限差分算子：
其中：
为了优化有限差分算子逼近微分算子的精度，可以选择Chebyshev窗函数去截断：
其中，r表示纹波率，代表旁瓣的衰减程度，N+1表示窗的长度，N为偶数。
不同窗函数幅频响应的主瓣和旁瓣性能影响着差分的逼近精度，具体影响方式包括：
1)主瓣大小和过渡带宽有关：主瓣窄，则过渡带窄，使用该窗函数截断逼近的有限差分算子的精度误差的谱覆盖范围大，可以用低阶算子达到高阶的精度，主瓣宽，则过渡带宽，使用该窗函数截断逼近的有限差分算子的精度误差的谱覆盖范围大。
2)旁瓣衰减和逼近精度稳定性的关系：旁瓣的衰减直接影响到了窗函数截断逼近的有限差分算子的精度误差的稳定性，旁瓣衰减越大，逼近精度误差波动越小，稳定性高，旁瓣衰减越小，逼近精度误差波动越大，稳定性低。
2、各向异性部分：
TI各向异性介质中，地震波的偏振向量的各向异性部分由线性部分和非线性部分构成，对于TTI介质，去掉非线性部分，也就是在弱各向异性的条件下，qP波的偏振角可以被近似表示为：
vp＝α+f[δ+2(ε-δ)sin2(α-θ)]sin 2(α-θ) (12)
其中，α表示相角，传播方向与Z轴的夹角，代表的是各向同性部分，θ表示TTI介质对称轴的倾角，
去除上述近似公式中的各向同性部分，只保留线性各向异性部分，可以得到：
Φ＝f[δ+2(ε-δ)sin2(α-θ)]sin 2(α-θ) (13)
其中，Φ和ε、δ成近似线性关系，和sin2(α-θ)成近似比例关系。
因此，可以设初始模型为m＝{ε，δ，θ}，在初始模型条件下偏振角的各向异性部分值为Φ，进一步的，可以条件选择N个参考模型mk，并计算得到每个参考模型对应的Φk，并通过以下公式插值得到Φ：
其中，wk表示权重系数。
IDW插值算法在计算权重时，仅考虑插值点与参考点的空间位置关系，也就是当不同参考点与插值点空间位置一致时，就有相同的权重系数，这种插值算法虽然执行效率较高，但是没有考虑参考点之间的差异。
在本例中，对权重系数的确定进行了改进，具体包括：假设Φ满足二阶平稳假设或者本征假设，使用变异函数μ(h)来描述Φ随着空间位置不同而变化的特性，按照以下公式计算得到变异函数值：
其中，表示位置增量，对于每一个都可以计算出一个以不同的的模为横坐标，为纵坐标，绘制点图，由绘制的点图可以发现：空间距离越小，变异性越小，越大，变异性越大，当增大到一个临界值R时，两个数值相互之间基本不会有太大的影响了，这也就说明在选取参考点，或者需要拟合的点的时候，应该在R范围内选取，否则选取的点没有太多的实际意义。这里描述得到的是零散的点，可以使用一些理论变异函数模型去拟合这些点，得到相关的参数。
引入变异函数R后可以得到权重系数wk：
因为μ(h)是距离的函数，因此将μ(h)引入IDW插值算法的权重计算中，可以获得更为可靠的插值结果。
对于插值参考点(εk,δk,θk)的选择，插值参考点的选择方法可包括：根据所述初始模型计算得到多个分散的变异函数值；对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；根据所述变异函数模型确定变异函数中的各个参数的参数值；遍历待数值模拟初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。
当然的，除了采用上述纵横波分离算法外，本发明的其它实施方式还提供其它类型的分离方法。例如，可以根据速度差异分离波场或者采用波动方程法等。因此，本申请并不对纵横波分离的算法进行限制，只要能将耦合纵横波分离的方法，都是符合本申请的要求。
参照图1所示，本申请的方法还包括：S104：经过纵横波分离之后的纵横波波场进行上下方向的方向波分解；
参照图4所示，经过纵横波分离之后的纵波波场进行上下方向的方向波分解，包括以下几个子步骤：
S1021：将分离出的纵波波场进行时间复数域拓展；
其中，复数域拓展策略是指将数据转化为复数的形式，进行各种计算的方法，该策略基于希尔伯特变换，可以有效地提高计算效率。由频谱分析可知，经过希尔伯特变换后，波场在频率域的振幅保持不变，相位倒转。对于经过复数域拓展的波场，可以发现在频率域，其振幅只有在正频率处有值，而且其相位保持不变。可以利用这个性质，将复数域拓展运用于全方位方向波分解，可以不必对时间进行傅里叶变换，只需要在空间方向进行傅里叶正反变换，有效较低计算量与存储量。
S1022：在Z方向对时间复数域拓展之后的纵波波场进行傅里叶变换，得到波数域纵波波场；
Kz表示纵横波场在Z方向上的波数。当地震波在地下介质传播时，定义向下为正，因此向下传播的方向波其Kz为正，向上传播的方向波，其Kz为负。
S1023：从所述波数域纵波波场中选择波数大于零的部分，波数小于零的部分置为零，进行Z方向上的反傅里叶变换，得到纵波波场的下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行Z方向上的反傅里叶变换，得到所述纵波波场的上行方向波波场。
参照图1所示，本申请的确定地震波场的坡印廷矢量的方法还包括步骤：
纵横波的上行方向波波场进行左右方向的方向波分解，分解为右上行方向波波场和左上行方向波波场；
在该步骤中，参照图5所示，所述纵波的上行方向波波场进行左右方向的方向波分解，分解为右上行方向波波场和左上行方向波波场，包括以下几个子步骤：
S1031：将分离出的所述纵波波场的上行方向波波场进行时间复数域拓展；
其中，Kx表示纵波纵横波场在X方向上的波数，定义向右为正，因此向右传播的方向波其Kx为正，向左传播的方向波，其Kx为负；
S1032：在X方向对时间复数域拓展之后的所述纵波波场的上行方向波波场进行傅里叶变换，得到波数域纵波波场；
S1033：从所述波数域纵波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到纵波波场的右上行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述纵波波场的左上行方向波波场。
参照图6所示，所述纵波的下行方向波波场进行左右方向的方向波分解，分解为右下行方向波波场和左下行方向波波场，包括以下几个子步骤：
S1041：将分离出的所述纵波波场的下行方向波波场进行时间复数域拓展；
S1042：在X方向对时间复数域拓展之后的所述纵波波场的下行方向波波场进行傅里叶变换，得到波数域纵波波场；
S1043：从所述波数域所述纵波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到纵波波场的右下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述纵波波场的左下行方向波波场。
参照图7所示，所述横波信号能在上下方向上被分离为上行横波信号和下行横波信号包括以下几个子步骤：
S1051：将分离出的横波波场进行时间复数域拓展；
S1052：在Z方向对时间复数域拓展之后的横波波场进行傅里叶变换，得到波数域横波波场；
S1053：从所述波数域横波波场中选择波数大于零的部分，波数小于零的部分置为零，进行Z方向上的反傅里叶变换，得到横波波场的下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行Z方向上的反傅里叶变换，得到所述横波波场的上行方向波波场。
该几个步骤中的复数域、Kz以及反傅里叶变换等已经在上文中进行解释，在此不再赘述。
在该步骤S104中，参照图8所示，所述上行横波信号能在左右方向上分离为左上行横波信号和右上行横波信号包括以下几个子步骤：
S1061：将分离出的所述横波波场的上行方向波波场进行时间复数域拓展；
S1062：在X方向对时间复数域拓展之后的所述横波波场的上行方向波波场进行傅里叶变换，得到波数域横波波场；
S1063：从所述波数域所述横波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到横波波场的右上行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述横波波场的左上行方向波波场。
该几个步骤中的复数域、Kx以及反傅里叶变换等已经在上文中进行解释，在此不再赘述。
参照图9所示，该步骤S104中，所述下行横波信号能在左右方向上分离为左下行横波信号和右下行横波信号包括以下几个子步骤：
S1071：将分离出的所述横波波场的下行方向波波场进行时间复数域拓展；
S1072：在X方向对时间复数域拓展之后的所述横波波场的下行方向波波场进行傅里叶变换，得到波数域横波波场；
S1073：从所述波数域所述横波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到横波波场的右下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述横波波场的左下行方向波波场。
该几个步骤中的复数域、Kx以及反傅里叶变换等已经在上文中进行解释，在此不再赘述。
参照图1所示，步骤S104为：对分离出的纵波波场和横波波场在上下左右四个方向进行分解，得到所述纵波波场和/或所述横波波场的多个方向波波场；
在该步骤中，耦合的弹性波波场根据纵横波偏振方向以及上下左右四个方向可以将耦合的弹性波场解耦并分为8个方向波(纵波4个方向波，横波4个方向波)。以炮点处的纵波为例，分离后的各波场记为
其中，S分别表示炮点处弹性波波场。上标p代表纵波。下标l,r代表左右，u,d代表上下。炮点处和检波点处的纵波按照上下左右四个方向分解后可以表示为：
同理，炮点处的横波按照上下左右四个方向分解后可以表示为：
以纵波为例，炮点处纵波的上下左右四个方向方向波可以通过以下公式进行表示：
其它位置以及其他类似的波场可以按照该公式的形式进行表示，在此不在赘述。
同时，为了降低计算量和存储量，可以将上述的数据拓展到复数域。对于拓展到复数域的数据，t＜0时值为零，因此可以不需要对数据进行时间Fourier变换，降低计算量，将公式(19)—公式(22)简化为：
在对单个耦合的弹性波场进行纵横波以及上下左右前后方向波的分解以及简化后，可以得到分离出来的纵横波场的8个方向波。利用分解后的各方向波计算坡印廷矢量，相对于现有技术中只对耦合弹性波场计算坡印廷矢量的实施方式相比，本申请的方法可以将单个的耦合弹性波波场解耦以及进行方向波分解，最终得到8个方向波波场，然后再对该8个方向波波场分别计算坡印廷矢量。
更进一步的，根据矢量叠加原理可知，波前重叠，波型耦合对波传播方向求取带来很大的影响，利用本申请的方法可以使求取的坡印廷矢量可以更好地代表波的传播方向，使求取的坡印廷矢量可以更好地代表波的传播方向，利用坡印廷矢量输出波场的反射角，使得角道集成像结果更加准确。例如，当将弹性波场分解为多个方向波以后，可采用以下公式(1)-(4)计算。
其中，P表示坡印廷矢量，下标x以及z表示方向波在空间位置的分量；v表示方向波的速度；τ表示应力张量。
实施例1：
参照图10(a)到图10(f)所示，分别为弹性波在纵横波分离以及在x分量以及z分量进行方向波分解前后的坡印廷矢量示意图。参照图10(a)到图10(c)所示，图10(a)为弹性波在x分量的坡印廷矢量示意图，图10(b)和图10(c)为将纵横波分离后的纵波以及横波在x方向分离后的坡印廷矢量示意图，可以看出坡印廷矢量得到很好地分离，其振幅和相位都得到了很好地保持。图10(d)到图10(f)将纵横波分离后的纵波以及横波在z方向分离后的坡印廷矢量示意图，与图10(d)到图(f)类似，在此不再赘述。
实施例2：
参照图11(a)到图11(g)所示，本实施例是采用点脉冲波响应去验证本申请的技术方案，可以看出图11(a)为该波的坡印廷矢量，图11(b)以及图(c)为将该波进行纵横波地分离后的坡印廷矢量，以及图11(d)到图11(g)为在将纵横波分离出的基础上再进行上下左右方向波的分离后计算得出的坡印廷矢量，可以看出采用本申请的方法可以将弹性波进行纵横波的分离以及上下左右方向波分分解，然后在此基础上计算所有方向波波场的坡印廷矢量；然后，利用坡印廷矢量输出波场的反射角，继而生成对应的ADCIGs道集。此方法避免了波前重叠，波型耦合对波传播方向求取带来的影响，生成的ADCIGs道集更加准确，成像效果更好。
参照图15所示，本申请还公开了一种确定地震波场的坡印廷矢量的计算装置，包括：构建单元401，用于构建地下炮点正传波场以及检波点反传波场；分离单元402，用于判断传播所述正传波场以及所述反传波场的介质类型；当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；计算单元403，从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
本申请还公开了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：构建地下炮点正传波场以及检波点反传波场；判断传播所述正传波场以及所述反传波场的介质类型；当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本申请时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(CPU)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(RAM)和/或非易失性内存等形式，如只读存储器(ROM)或闪存(flash RAM)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitory media)，如调制的数据波场和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。
本领域技术人员应明白，本申请的实施例可提供为方法、系统或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。
Claims (14)
1.一种确定地震波场的坡印廷矢量的方法，其特征在于，包括：
构建地下炮点正传波场以及检波点反传波场；
判断传播所述正传波场以及所述反传波场的介质类型；
当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；
对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
从所述多个方向波波场中选择出多个方向波波场，并计算选择出的所述多个方向波波场的坡印廷矢量。
2.如权利要求1所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述对分离出的纵波波场和或横波波场在上下左右四个方向进行分解包括：
对分离出的纵波波场和横波波场进行上下方向的方向波分离，得到对应的上行方向波波场和下行方向波波场；
对上行方向波波场和下行方向波波场进行左右方向的方向波分离，得到对应的左上行、左下行方向波波场，和右上行、右下行方向波波场。
3.如权利要求2所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述对分离出的纵波波场进行上下方向的方向波分解，包括：
将分离出的纵波波场进行时间复数域拓展；
在Z方向对时间复数域拓展之后的纵波波场进行傅里叶变换，得到波数域纵波波场；
从所述波数域纵波波场中选择波数大于零的部分，波数小于零的部分置为零，进行Z方向上的反傅里叶变换，得到纵波波场的下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行Z方向上的反傅里叶变换，得到所述纵波波场的上行方向波波场。
4.如权利要求3所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述对上行方向波波场进行左右方向的方向波分解，包括：
将所述纵波波场的上行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述纵波波场的上行方向波波场进行傅里叶变换，得到波数域纵波波场；
从所述波数域纵波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到所述纵波波场的右上行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述纵波波场的左上行方向波波场。
5.如权利要求3所述的确定地震波场的坡印廷矢量的方法，其特征在于，对下行方向波波场进行左右方向的方向波分解，包括：
将所述纵波波场的下行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述纵波波场的下行方向波波场进行傅里叶变换，得到波数域纵波波场；
从所述波数域纵波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向下的反傅里叶变换，得到所述纵波波场的右下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向下的反傅里叶变换，得到所述纵波波场的左下行方向波波场。
6.如权利要求2所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述对分离出的横波波场进行上下方向的方向波分解，包括：
将分离出的横波波场进行时间复数域拓展；
在Z方向对时间复数域拓展之后的横波波场进行傅里叶变换，得到波数域横波波场；
从所述波数域横波波场中选择波数大于零的部分，波数小于零的部分置为零，进行Z方向上的反傅里叶变换，得到横波波场的下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行Z方向上的反傅里叶变换，得到所述横波波场的上行方向波波场。
7.如权利要求6所述的确定地震波场的坡印廷矢量的方法，其特征在于，对上行方向波波场进行左右方向的方向波分解，包括：
将所述横波波场的上行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述横波波场的上行方向波波场进行傅里叶变换，得到波数域横波波场；
从所述波数域横波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向上的反傅里叶变换，得到所述横波波场的右上行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向上的反傅里叶变换，得到所述横波波场的左上行方向波波场。
8.如权利要求6所述的确定地震波场的坡印廷矢量的方法，其特征在于，对下行方向波波场进行左右方向的方向波分解，包括：
将所述横波波场的下行方向波波场进行时间复数域拓展；
在X方向对时间复数域拓展之后的所述横波波场的下行方向波波场进行傅里叶变换，得到波数域横波波场；
从所述波数域横波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行X方向下的反傅里叶变换，得到所述横波波场的右下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行X方向下的反傅里叶变换，得到所述横波波场的左下行方向波波场。
9.如权利要求1所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离，包括：
当在非均匀各向异性介质下时，利用傅里叶变换分别将各时刻的所述正传波场和或反传波场由空间域变换至波数域；
选取N个参考模型，根据所述参考模型计算分离算子，并在波数域利用自褶积组合窗函数对所述分离算子进行截断优化；
在波数域利用所述分离算子对所述参考模型下的所述正传波场和或反传波场进行分离，得到所述参考模型下的分离后的所述纵横波波场；
将所述参考模型下的分离后的纵横波波场由波数域变换至空间域，N个所述参考模型对应得到N个所述参考模型下的分离后的纵横波波场，N为正整数；
在空间域，利用变异函数改进IDW插值算法并根据改进后的所述IDW插值算法计算N个参考模型的权重系数，对N个所述参考模型下的分离后的纵横波波场进行加权插值处理，得到最终分离后的纵波波场和横波波场。
10.如权利要求9所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述选取N个参考模型，包括：
遍历待数值模拟初始模型的多个各向异性参数，并利用变异函数的临界值进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到所述参考模型。
11.如权利要求10所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述参考模型的权重系数通过所述变异函数计算得到。
12.如权利要求11所述的确定地震波场的坡印廷矢量的方法，其特征在于，所述变异函数的获取方法包括：
通过计算所述初始模型的各向异性参数分布的变异函数值，拟合得到所述初始模型的变异函数。
13.一种确定地震波场的坡印廷矢量的装置，其特征在于，包括：
构建单元，用于构建地下炮点正传波场以及检波点反传波场；
分离单元，用于判断传播所述正传波场以及所述反传波场的介质类型；当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
计算单元，从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
14.一种计算机存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现以下步骤：
构建地下炮点正传波场以及检波点反传波场；
判断传播所述正传波场以及所述反传波场的介质类型；
当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数IDW插值算法，对所述正传波场和或反传波场进行纵横波场分离；
对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
从所述多个方向波波场中选择出多个方向波波场，并计算选择出的方向波波场的坡印廷矢量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710542740.1A CN107153216B (zh) | 2017-07-05 | 2017-07-05 | 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710542740.1A CN107153216B (zh) | 2017-07-05 | 2017-07-05 | 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107153216A true CN107153216A (zh) | 2017-09-12 |
| CN107153216B CN107153216B (zh) | 2019-05-07 |
Family
ID=59796247
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710542740.1A Expired - Fee Related CN107153216B (zh) | 2017-07-05 | 2017-07-05 | 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107153216B (zh) |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107894613A (zh) * | 2017-10-26 | 2018-04-10 | 中国石油天然气集团公司 | 弹性波矢量成像方法、装置、存储介质及设备 |
| CN108196303A (zh) * | 2017-12-29 | 2018-06-22 | 中国石油天然气集团公司 | 弹性波波场分离方法、装置、存储介质及设备 |
| CN108345032A (zh) * | 2018-01-12 | 2018-07-31 | 中国科学技术大学 | 一种弱照明区域高信噪比偏移成像方法 |
| CN109212605A (zh) * | 2018-09-28 | 2019-01-15 | 中国科学院地质与地球物理研究所 | 拟微分算子储存方法及装置 |
| CN110596754A (zh) * | 2019-09-24 | 2019-12-20 | 中国矿业大学(北京) | 一种三维TTI介质qP波与qSV波波场模拟方法 |
| CN110749925A (zh) * | 2018-07-24 | 2020-02-04 | 中国石油化工股份有限公司 | 一种宽频逆时偏移成像处理方法 |
| CN111427984A (zh) * | 2020-03-24 | 2020-07-17 | 成都理工大学 | 一种区域地震概率空间分布生成方法 |
| CN112578455A (zh) * | 2019-09-30 | 2021-03-30 | 中国石油化工股份有限公司 | 一种空间波数混合域地震波场多方向分解方法及系统 |
| CN113031062A (zh) * | 2021-04-09 | 2021-06-25 | 中国海洋大学 | 一种基于波场分离的相关加权逆时偏移成像方法 |
| CN114444299A (zh) * | 2022-01-24 | 2022-05-06 | 中国科学院国家空间科学中心 | 一种基于距离加权多极展开方法的磁场重构方法 |
| CN114624766A (zh) * | 2022-05-16 | 2022-06-14 | 中国海洋大学 | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 |
| CN114924313A (zh) * | 2022-05-16 | 2022-08-19 | 中国海洋大学 | 基于行波分离的声波最小二乘逆时偏移梯度求取方法 |
| CN116381783A (zh) * | 2023-03-30 | 2023-07-04 | 东北石油大学 | 一种稳定高效的坡印廷矢量提取方法及系统 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130064431A1 (en) * | 2010-06-02 | 2013-03-14 | Graham A. Winbow | Efficient Computation of Wave Equation Migration Angle Gathers |
| CN104133241A (zh) * | 2014-07-31 | 2014-11-05 | 中国科学院地质与地球物理研究所 | 波场分离方法和装置 |
| CN105137486A (zh) * | 2015-09-01 | 2015-12-09 | 中国科学院地质与地球物理研究所 | 各向异性介质中弹性波逆时偏移成像方法及其装置 |
| CN106918840A (zh) * | 2017-05-05 | 2017-07-04 | 中国石油化工股份有限公司 | 基于角度波场分步提取的逆时偏移角道集成像方法 |
-
2017
- 2017-07-05 CN CN201710542740.1A patent/CN107153216B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130064431A1 (en) * | 2010-06-02 | 2013-03-14 | Graham A. Winbow | Efficient Computation of Wave Equation Migration Angle Gathers |
| CN104133241A (zh) * | 2014-07-31 | 2014-11-05 | 中国科学院地质与地球物理研究所 | 波场分离方法和装置 |
| CN105137486A (zh) * | 2015-09-01 | 2015-12-09 | 中国科学院地质与地球物理研究所 | 各向异性介质中弹性波逆时偏移成像方法及其装置 |
| CN106918840A (zh) * | 2017-05-05 | 2017-07-04 | 中国石油化工股份有限公司 | 基于角度波场分步提取的逆时偏移角道集成像方法 |
Cited By (20)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107894613B (zh) * | 2017-10-26 | 2019-07-26 | 中国石油天然气集团公司 | 弹性波矢量成像方法、装置、存储介质及设备 |
| CN107894613A (zh) * | 2017-10-26 | 2018-04-10 | 中国石油天然气集团公司 | 弹性波矢量成像方法、装置、存储介质及设备 |
| CN108196303A (zh) * | 2017-12-29 | 2018-06-22 | 中国石油天然气集团公司 | 弹性波波场分离方法、装置、存储介质及设备 |
| CN108196303B (zh) * | 2017-12-29 | 2019-10-01 | 中国石油天然气集团公司 | 弹性波波场分离方法、装置、存储介质及设备 |
| CN108345032A (zh) * | 2018-01-12 | 2018-07-31 | 中国科学技术大学 | 一种弱照明区域高信噪比偏移成像方法 |
| CN110749925A (zh) * | 2018-07-24 | 2020-02-04 | 中国石油化工股份有限公司 | 一种宽频逆时偏移成像处理方法 |
| CN109212605A (zh) * | 2018-09-28 | 2019-01-15 | 中国科学院地质与地球物理研究所 | 拟微分算子储存方法及装置 |
| CN110596754A (zh) * | 2019-09-24 | 2019-12-20 | 中国矿业大学(北京) | 一种三维TTI介质qP波与qSV波波场模拟方法 |
| CN112578455A (zh) * | 2019-09-30 | 2021-03-30 | 中国石油化工股份有限公司 | 一种空间波数混合域地震波场多方向分解方法及系统 |
| CN111427984B (zh) * | 2020-03-24 | 2022-04-01 | 成都理工大学 | 一种区域地震概率空间分布生成方法 |
| CN111427984A (zh) * | 2020-03-24 | 2020-07-17 | 成都理工大学 | 一种区域地震概率空间分布生成方法 |
| CN113031062A (zh) * | 2021-04-09 | 2021-06-25 | 中国海洋大学 | 一种基于波场分离的相关加权逆时偏移成像方法 |
| CN113031062B (zh) * | 2021-04-09 | 2022-01-28 | 中国海洋大学 | 一种基于波场分离的相关加权逆时偏移成像方法 |
| CN114444299A (zh) * | 2022-01-24 | 2022-05-06 | 中国科学院国家空间科学中心 | 一种基于距离加权多极展开方法的磁场重构方法 |
| CN114444299B (zh) * | 2022-01-24 | 2022-09-13 | 中国科学院国家空间科学中心 | 一种基于距离加权多极展开方法的磁场重构方法 |
| CN114624766A (zh) * | 2022-05-16 | 2022-06-14 | 中国海洋大学 | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 |
| CN114924313A (zh) * | 2022-05-16 | 2022-08-19 | 中国海洋大学 | 基于行波分离的声波最小二乘逆时偏移梯度求取方法 |
| CN114924313B (zh) * | 2022-05-16 | 2022-12-13 | 中国海洋大学 | 基于行波分离的声波最小二乘逆时偏移梯度求取方法 |
| CN116381783A (zh) * | 2023-03-30 | 2023-07-04 | 东北石油大学 | 一种稳定高效的坡印廷矢量提取方法及系统 |
| CN116381783B (zh) * | 2023-03-30 | 2023-10-03 | 东北石油大学 | 一种稳定高效的坡印廷矢量提取方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107153216B (zh) | 2019-05-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107153216B (zh) | 确定地震波场的坡印廷矢量的方法、装置以及计算机存储介质 | |
| CN107272058A (zh) | 成像方法、成像装置以及计算机存储介质 | |
| CN104133241B (zh) | 波场分离方法和装置 | |
| CN105137486B (zh) | 各向异性介质中弹性波逆时偏移成像方法及其装置 | |
| WO2015199800A1 (en) | Fast viscoacoustic and viscoelastic full-wavefield inversion | |
| CN105425289B (zh) | 确定低频波阻抗的方法和装置 | |
| CN107340540B (zh) | 弹性波场的方向波分解方法、装置以及计算机存储介质 | |
| CN107798156B (zh) | 一种频率域2.5维粘弹性波数值模拟方法及装置 | |
| KR20180067650A (ko) | 진폭 보존을 갖는 fwi 모델 도메인 각도 스택들 | |
| CN110954950B (zh) | 地下横波速度反演方法、装置、计算设备及存储介质 | |
| CN111766628A (zh) | 一种预条件的时间域弹性介质多参数全波形反演方法 | |
| CN107894613A (zh) | 弹性波矢量成像方法、装置、存储介质及设备 | |
| Tavakoli F et al. | Matrix-free anisotropic slope tomography: Theory and application | |
| CN110058303A (zh) | 声波各向异性逆时偏移混合方法 | |
| Vamaraju et al. | Enriched Galerkin finite element approximation for elastic wave propagation in fractured media | |
| CN110879412A (zh) | 地下横波速度反演方法、装置、计算设备及存储介质 | |
| CN113341455A (zh) | 一种粘滞各向异性介质地震波数值模拟方法、装置及设备 | |
| CN105182414B (zh) | 一种基于波动方程正演去除直达波的方法 | |
| de Jonge et al. | Debubbling seismic data using a generalized neural network | |
| Zuberi et al. | Mitigating nonlinearity in full waveform inversion using scaled-Sobolev pre-conditioning | |
| US11199641B2 (en) | Seismic modeling | |
| Sun et al. | Multiple attenuation using λ-f domain high-order and high-resolution Radon transform based on SL0 norm | |
| CN103901473A (zh) | 一种基于非高斯性最大化的双检信号上下行波场分离方法 | |
| CN114428343A (zh) | 基于归一化互相关的Marchenko成像方法及系统 | |
| CN105095555A (zh) | 一种速度场的无散平滑处理方法及装置 |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190507 |