弹性波矢量成像方法、装置、存储介质及设备与流程

本发明涉及地震勘探技术领域，尤其涉及一种弹性波矢量成像方法、装置、存储介质及设备。

背景技术：

近年来计算机得到了飞速提高，这使得叠前逆时偏移成为解决复杂构造和强速度对比地区成像的有力工具。地震波在地下弹性介质传播中既包含纵波，也包含横波，若仅考虑纵波，则是默认地下介质为液态，不符合实际介质情况。近年来，随着采集技术的进步，多分量资料的记录得到了一定发展。使用较多的一种传统处理方法是，先在地表对多分量资料进行分离，分离为纵波、横波资料，然后再使用标量方程进行纵波和横波的偏移及成像。严格来讲，因为纵、横波地震资料中存在sp、ps转换波，所以直接用标量声波方程对多分量资料进行延拓处理并不能正确地反映地下介质的弹性、矢量信息。另一种传统处理方法是直接利用笛卡尔分量分别进行互相关成像，然而，因为纵横波波长不同且存在转换波效应，所以存在的串扰噪声以及非相干同相轴的叠加会得到很不理想的叠加成像剖面。这就使得有必要使用弹性波方程对多分量资料进行延拓，并在传播过程中进行波场分离，然后分别对纯p和s波场互相关成像，得到物理意义更明确的成像结果。

在地球物理学中，广泛存在地震各向异性。这种性质往往由地球内部的薄互层、定向排列的裂缝和裂隙、地层沉积压实等因素引起。地震资料处理中有必要考虑到各向异性的因素。将弹性逆时偏移推广到各向异性中，波场分离技术是一个关键环节。目前这方面的研究正如火如荼。在传统各向同性介质中，根据亥姆霍茨理论，波场纵、横波可以通过散度、旋度的求取得以分离。但是由于在各向异性介质中波传播的方向发生了偏离，其振动方向并不再是纵波与其相同、横波与其垂直的关系，所以使用亥姆霍兹分解并不能很好地将各向异性介质中的耦合纵横波分离。

对于各向异性介质，传统做法是通过在波数域求解christoffel方程，然后将其变换回空间域，得到拟散度、拟旋度算子，进而实现波场分离；在非均匀介质中，通过将局部极化方向向量定义为非稳态空间滤波器，得到拟散度、拟旋度算子。或者是做一些近似来减少计算量。改进的方法也大都是基于散度和旋度的思路，以上方法的缺陷是均需在规则网格的空间域完成而且涉及多次求导破坏波场信息，同时不能很好地应用于非规则网格的地震数值模拟算法，例如谱元法、格子法。

技术实现要素：

本发明提供了一种弹性波矢量成像方法、装置、存储介质及设备，以提高成像效率和准确度。

本发明实施例提供一种弹性波矢量成像方法，包括：利用成像区域中各点各向异性介质的设定方向的p波速度、设定方向的s波速度以及各向异性参数，建立群角与极化角之间的转换关系；应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，得到震源正传波场和炮记录反传波场，以及相应的坡印廷矢量；按设定时刻间距，利用坡印廷矢量和群角与极化角之间的转换关系，对震源正传波场和炮记录反传波场进行波场分解；利用多分量地震数据中所有炮对应的波场分解结果进行各向异性介质中弹性波矢量成像。

本发明实施例还提供一种弹性波矢量成像装置，包括：转换关系建立单元，用于：利用成像区域中各点各向异性介质的设定方向的p波速度、设定方向的s波速度以及各向异性参数，建立群角与极化角之间的转换关系；波场延拓单元，用于：应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，得到震源正传波场和炮记录反传波场，以及相应的坡印廷矢量；波场分解单元，用于：按设定时刻间距，利用坡印廷矢量和群角与极化角之间的转换关系，对震源正传波场和炮记录反传波场进行波场分解；弹性波矢量成像单元，用于：利用多分量地震数据中所有炮对应的波场分解结果进行各向异性介质中弹性波矢量成像。

本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现所述方法的步骤。

本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现所述方法的步骤。

本发明实施例的弹性波矢量成像方法、装置、存储介质及设备，不对接收的多分量资料进行波场分解，而是应用弹性波方程直接进行多分量资料延拓深度，并利用建立的群角与极化角之间的转换关系进行波场分解和成像，以此能够获得更准确的纵波及转换波界面反射信息。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：

图1是本发明实施例的弹性波矢量成像方法的流程示意图；

图2是本发明一实施例中建立群角与极化角之间的转换关系的方法流程示意图；

图3是本发明一实施例中计算得到群角与极化角之间的转换关系的方法流程示意图；

图4是本发明一实施例中进行多分量地震数据的波场深度延拓的方法流程示意图；

图5是本发明一实施例中得到多炮可用于逆时偏移的炮记录的方法流程示意图；

图6是本发明一实施例中进行波场分解的方法流程示意图；

图7是本发明一实施例中根据群角值和群角与极化角之间的转换关系获取极化角角度的方法流程示意图；

图8是本发明一实施例中进行各向异性介质中弹性波矢量成像的方法流程示意图；

图9是本发明一实施例的弹性波矢量成像方法流程示意图；

图10是本发明一实施例中成像区域内沿对称轴方向的p波及s波速度和各向异性参数的构造示意图；

图11是本发明一实施例中正传波场快照示意图；

图12是本发明一实施例中反传波场快照示意图；

图13是图12的反传波场的p波和s波分解结果示意图；

图14是图10地质构造的弹性波矢量成像的偏移成像剖面示意图；

图15是本发明一实施例中角道集示意图；

图16是本发明一实施例的弹性波矢量成像装置的结构示意图；

图17是本发明一实施例中转换关系建立单元的结构示意图；

图18是本发明一实施例中群角与极化角转换关系生成模块的结构示意图；

图19是本发明一实施例中波场延拓单元的结构示意图；

图20是本发明一实施例中波场分解单元的结构示意图；

图21是本发明一实施例中极化角角度生成模块的结构示意图；

图22是本发明一实施例中弹性波矢量成像单元的结构示意图；

图23是本发明一实施例的计算机设备的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

图1是本发明实施例的弹性波矢量成像方法的流程示意图。如图1所示，本发明实施例的弹性波矢量成像方法，可包括：

步骤s110：利用成像区域中各点各向异性介质的设定方向的p波速度、设定方向的s波速度以及各向异性参数，建立群角与极化角之间的转换关系；

步骤s120：应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，得到震源正传波场和炮记录反传波场，以及相应的坡印廷矢量；

步骤s130：按设定时刻间距，利用坡印廷矢量和群角与极化角之间的转换关系，对震源正传波场和炮记录反传波场进行波场分解；

步骤s140：利用多分量地震数据中所有炮对应的波场分解结果进行各向异性介质中弹性波矢量成像。

在上述步骤s110中，该设定方向可以是多种不同方向，例如可以是对称横方向和沿对称轴方向。实施例中，可以利用各向异性速度建模方法，计算得到成像区域中各点各向异性介质的设定方向的p波速度、设定方向的s波速度以及各向异性参数。针对不同的地震波，例如p波和sv波，可以分别建立其各自的群角与极化角之间的转换关系。极化角可指质点振动极化角。

在上述步骤s120中，多分量地震数据可以是二分量地震数据或三分量地震数据。该步骤中，可采用逆时偏移方法进行波场深度延拓，如此一来，可精确考虑复杂速度、构造及各向异性，实现复杂地质条件下的准确成像。

通过上述步骤s130和步骤s140，可针对地下成像区域中的各点、各时刻、各种波场、各种波进行波场分解和成像。

本发明实施例的方法可以是各向异性介质中弹性波矢量成像方法，可应用于地震勘探中多分量反射地震资料处理，可针对二维、二分量采集地震资料的偏移成像方法。

本实施例中，不是对接收的多分量资料进行波场分解，而是应用弹性波方程直接进行多分量资料延拓深度，并利用建立的群角与极化角之间的转换关系进行波场分解和成像，以此能够获得更准确的纵波及转换波界面反射信息。应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，能够精确考虑复杂速度、构造及各向异性，实现复杂地质条件下的准确成像。

图2是本发明一实施例中建立群角与极化角之间的转换关系的方法流程示意图。如图2所示，在上述步骤s110中，利用成像区域中各点各向异性介质的设定方向的p波速度、设定方向的s波速度以及各向异性参数，建立群角与极化角之间的转换关系的方法，可包括：

步骤s111：对christoff方程求特征值，并结合设定方向的p波速度、设定方向的s波速度以及各向异性参数，分别针对p波和sv波，计算得到以相角为自变量的相速度，其中，设定方向包括对称横方向和沿对称轴方向；

步骤s112：分别针对p波和sv波，将相速度和相速度对相角的偏导数代入群角与相角之间的关系式，计算得到群角与相角之间的转换关系；

步骤s113：分别针对p波和sv波，在christoff方程中求解极化角，并结合各向异性参数、相速度及群角与相角之间的转换关系，计算得到群角与极化角之间的转换关系。

在上述步骤s111中，根据christoff方程可以求得两个特征值，可分别对应p波的相速度和sv波的相速度。将设定方向的p波速度、设定方向的s波速度以及各向异性参数代入所得特征值中可得到以相角为自变量的相速度。

在上述步骤s112中，将p波的相速度和p波的相速度对p波的相角的偏导数代入p波对应的群角与相角之间的关系式，可计算得到p波对应的群角与相角之间的转换关系。将sv波的相速度和sv波的相速度对sv波的相角的偏导数代入sv波对应的群角与相角之间的关系式，可计算得到sv波对应的群角与相角之间的转换关系。

在上述步骤s113中，可以利用christoff方程分别求解p波的极化方向向量对应的极化角度和sv波的极化方向向量对应的极化角度，所得极化角度可以是以相角为自变量。将群角与相角之间的转换关系代入所得极化角度，并将各向异性参数、相速度代入所得极化角度，可以得到群角与极化角之间的转换关系。

图3是本发明一实施例中计算得到群角与极化角之间的转换关系的方法流程示意图。如图3所示，在上述步骤s113中，分别针对p波和sv波，在christoff方程中求解极化角，并结合各向异性参数、相速度及群角与相角之间的转换关系，计算得到群角与极化角之间的转换关系的方法，可包括：

步骤s1131：分别针对p波和sv波，在设定角度范围内按设定角度间隔取群角角度值；

步骤s1132：分别针对p波和sv波，在christoff方程中求解极化角，并结合各向异性参数、相速度、群角角度值及群角与相角之间的转换关系，计算得到群角和极化角之间的数值对应表，作为群角与极化角之间的转换关系。

实施例中，群角和极化角之间的数值对应表包含各向异性参数范围内的群角和极化角，每一组各向异性参数均对应一组群角和极化角的对应关系。

实施例中，设定角度范围可以是0～90度。设定角度间隔可以例如是1度。

因为波场每更新一次，坡应廷矢量就重新算一次，那么对应的群角也就更新一次。如果用解析的方法去求群角和极化角的关系，在每一个时间步更新时，就需要求一次。所以，在本实施例中，提前将各向异性参数范围内的群角和极化角求取，得到的群角和极化角之间的数值对应表，可以基于查表的方法进行波场分解，如此一来，在处理复杂各向异性介质的地震数据时基本不会增加计算量，从而能够提高计算效率。

实施例中，波场延拓过程可以基于查群角和极化角之间的数值对应表和插值进行。

实施例中，在上述步骤s111中，分别针对p波和sv波，对christoff方程求特征值，并结合设定方向的p波速度、设定方向的s波速度以及各向异性参数，计算得到以相角为自变量的相速度的方法，可包括：分别针对p波和sv波，将相速度和相速度对相角的偏导数代入群角与相角之间的关系式，计算得到以群角为自变量且以相角为因变量的关系式；分别针对p波和sv波，将以群角为自变量且以相角为因变量的关系式作为群角与相角之间的转换关系；或者，分别针对p波和sv波，在设定角度范围内按设定角度间隔取群角角度值，并将群角角度值代入以群角为自变量且以相角为因变量的关系式，计算得到群角和相角之间的数值对应表，作为群角与相角之间的转换关系。

图4是本发明一实施例中进行多分量地震数据的波场深度延拓的方法流程示意图。如图4所示，在上述步骤s120中，应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，得到震源正传波场和炮记录反传波场，以及相应的坡印廷矢量的方法，可包括：

步骤s121：对多分量地震数据进行预处理，得到多炮可用于逆时偏移的炮记录；

步骤s122：由每炮的炮记录的频宽生成宽频子波，并利用宽频子波进行数值模拟，生成成像区域的边界记录；

步骤s123：针对每炮的边界记录和炮记录，利用弹性波方程数值解法进行弹性波场反传播，并在不同时刻将边界记录赋值给成像区域的边界位置，得到不同时刻的震源正传波场及相应的坡印廷矢量，在不同时刻将炮记录赋值给检波点位置，得到不同时刻的炮记录反传波场及相应的坡印廷矢量。

实施例中，上述步骤s122的具体实时方式可以是：确定炮的频谱，将其转换到时间域得到一个类似雷克子波的子波，并利用该子波进行正演模拟生成边界记录。

实施例中，可以采用宽时间和宽空间记录生成的边界记录，用于波场延拓。以此可以降低弹性逆时偏移的带来内存的消耗和数据的读取带来计算量的巨大消耗的问题。

图5是本发明一实施例中得到多炮可用于逆时偏移的炮记录的方法流程示意图。如图5所示，在上述步骤s121中，对多分量地震数据进行预处理，得到多炮可用于逆时偏移的炮记录的方法，可包括：

步骤s1211：用测线在地表记录人工震源激发的多分量地震数据，其中，多分量地震数据为二分量地震数据；

步骤s1212：从二分量地震数据中抽取包含垂直速度分量和水平速度分量的炮记录；

步骤s1213：对抽取的炮记录衰减面波和直达波，得到可用于逆时偏移的炮记录。

本实施例中，通过上述步骤s1213可以得到仅包含反射波的炮记录，用于逆时偏移。

实施例中，上述步骤s1212中，从二分量地震数据中抽取包含垂直速度分量和水平速度分量的炮记录后，可按炮点位置排序并存储炮记录。

图6是本发明一实施例中进行波场分解的方法流程示意图。如图6所示，在上述步骤s130中，按设定时刻间距，利用坡印廷矢量和群角与极化角之间的转换关系，对震源正传波场和炮记录反传波场进行波场分解的方法，可包括：

步骤s131：分别针对p波和sv波，利用设定时刻间距的坡印廷矢量计算得到相应的群角值；

步骤s132：分别针对p波和sv波，根据群角值和群角与极化角之间的转换关系，获取成像区域中各点的极化角角度；

步骤s133：分别针对p波和sv波，将极化角角度代入波场分解公式，计算得到伪p波和伪s波，作为波场分解结果。

实施例中，在上述步骤s133中，波场分解公式可以是：

uqp＝uxcosυ+uzsinυ

uqs＝-uxsinυ+uzcosυ

其中，uqp和uqs分别表示伪p波和伪s波，ux表示p波或sv波的x轴分量，uz表示p波或sv波的z轴分量，υ表示p波的极化角或sv波极化角。

图7是本发明一实施例中根据群角值和群角与极化角之间的转换关系获取极化角角度的方法流程示意图。如图7所示，在上述步骤s132中，分别针对p波和sv波，根据群角值和群角与极化角之间的转换关系，获取极化角角度的方法，可包括：

步骤s1321：分别针对p波和sv波，利用群角值查找群角与极化角之间的转换关系的数值对应表，得到相应的极化角值；

步骤s1322：分别针对p波和sv波，利用得到的极化角值进行插值，得到成像区域中各点的极化角角度。

本实施例中，通过依据群角值查表得到极化角值，并进而进行插值得到成像区域中所有点的极化角角度，可以显著降低计算量。

图8是本发明一实施例中进行各向异性介质中弹性波矢量成像的方法流程示意图。如图8所示，在上述步骤s140中，利用多分量地震数据中所有炮对应的波场分解结果进行各向异性介质中弹性波矢量成像的方法，可包括：

步骤s141：在设定相角范围内，利用波场分解结果得到多分量地震数据中单炮的成像剖面、偏移角道集及照明强度；

步骤s142：将多分量地震数据中所有炮的成像剖面、偏移角道集及照明强度分别累加，得到累加的成像剖面、累加的偏移角道集及累加的照明强度；

步骤s143：用累加的照明强度除以累加的成像剖面得到各向异性介质中弹性波矢量成像的成像剖面，并由累加的偏移角道集得到各向异性介质中弹性波矢量成像的偏移角道集。

在上述步骤s141中，偏移角道集可以包括pp波的偏移角道集和ps波的偏移角道集。该设定相角范围例如可以是[-π/4,π/4]。通过在设定相角范围内获取单炮的成像剖面，可以达到压制低频噪音的效果。通过在设定相角范围内获取照明强度，可以防止照明补偿过量。

实施例中，可以利用包括类似步骤s111和步骤s112的方法得到群角与相角之间的转换关系，然后结合类似步骤s1131所取的群角角度值得到群角与相角的数据对应表，进而可以在设定相角范围内取相角角度，用于获取单炮的成像剖面、偏移角道集及照明强度。

具体实施例中，多分量地震数据可以是二分量反射地震资料，二分量可以包括垂直和水平方向。设定方向可包括对称横方向和沿对称轴方向。偏移角道集可以包括pp波的偏移角道集和ps波的偏移角道集。图9是本发明一实施例的弹性波矢量成像方法流程示意图。如图9所示，该具体实施例中，弹性波矢量成像方法，可包括：

步骤s210：用测线在地表记录人工震源激发的二分量反射地震资料，抽取包含垂直和水平速度分量的炮记录，按炮点位置排序存放到计算机中；

步骤s220：对记录的炮记录衰减面波和直达波；利用各向异性速度建模方法，求得成像区域中各点各向异性介质的对称横方向、沿对称轴方向的p波和s波速度和各向异性参数；

步骤s230：利用成像区域中各点各向异性介质的对称横方向、沿对称轴方向的p波和s波速度和各向异性参数，建立相速度和群速度角度以及质点振动极化角的角度转换表；

步骤s240：对炮记录依次求解，由炮记录的频宽生成宽频子波，数值模拟生成边界记录，生成子波褶积炮记录；

步骤s250：针对每炮的边界记录和子波褶积炮记录(针对步骤s240得到的每炮的边界记录和步骤s220得到的衰减后的炮记录)，利用弹性波方程数值解法进行弹性波场反传播，在不同时刻将边界记录赋值给步骤s240指定的边界位置，可得到不同时刻的震源正传波场，在不同时刻将子波褶积炮记录赋值给检波点位置，得到不同时刻的炮记录反传波场；

步骤s260：针对每炮的各个时刻的震源正传波场和炮记录反传波场，按指定的时刻间距，进行波场分解和成像，给出单炮的纵波和转换波成像剖面、偏移角道集和照明强度；

步骤s270：将全部炮的纵波和转换波成像剖面、偏移角道集和照明强度进行累加，用累加的照明强度除以累加的纵波和转换波成像剖面，得到各向异性介质中弹性波矢量成像的成像剖面和偏移角道集。

成像剖面和偏移角道集可用于指示地下构造的形态、断裂部位和反射特征，用于确定地下生、储油构造和识别油气储层。

具体实施例中，数值对应表的具体实现过程可如下：

在各向异性介质中，群速度的方向和相速度方向不相同。在二维情况下，群角和相角之间的关系可以表示为：

其中，ψ表示群速度的角度，θ表示相速度的角度，v是波场传播相速度。

由公式(1)求出的是相角和群角的对应关系。为了便于后续求解群角和极化角之间的对应关系，实施例中，可根据由公式(1)求出的相角和群角的对应关系得到群角和相角的对应关系。发明人通过分析发现公式(1)在给定范围内是一个单调函数，所以在实施例中，可以用数值的方法求出公式(1)的反函数。

实施例中，数值求解反函数的代码可如下：

其中，ang_p[]中存储的是相角按每一度的间隔对应的群角，求取完反函数以后ang[]中存储的是群角按每一度的间隔对应的相角。

在求解反函数的过程中，还需要用到相速度v以及其导数在后续内容中会给出求解相速度v和相速度对相角偏倒数的具体实施例。实施例中，根据求得的相速度v和相速度对相角偏倒数，可以在每一组各向异性参数下得到群角(例如[0-90]度)和相角(例如[0-90])的对应关系：

θp＝f(ψp),θsv＝f(ψsv)(2)

其中，θp和θsv分别表示p波和sv波的相角，ψp和ψsv分别表示p波和sv波的群角。

可以将公式(2)所示的关系暂记作表一。根据表一可以进一步得到数据对应表一。实施例中，在成像过程得到群角后可以通过查表得到对应的相角产生共成像点偏移道集。

实施例中，在数据对应表一中可只存0-90度群角/相角，90-180度的群角/相角可以变换到0-90度的群角/相角区间求解。

实施例中，可以通过求解christoff方程求解各向异性参数对应的相速度。实施例中，可以通过求解christoff方程求解2dvti介质对应的相速度，tti介质对应的相速度可以通过角度旋转到vti介质的数值表得到。求解christoff方程如下：

其中，c11、c55、c13、c33表示弹性张量，ρ表示介质密度，kx和kz分别表示x方向和z方向的波数，u1和u3分别表示极化方向向量在x方向和z方向的分量，ω表示角速度，v表示波速。

考虑到且对于后续求解无影响，实施例中，为了便于化简，可将kx和kz记kx＝cosθ，kz＝sinθ，那么公式(3)可变为：

令公式(4)的系数行列式为零：

得到对应的两个特征值，分别对应p波对应的相速度和sv波对应的相速度：

其中，公式(6)中，v表示p波对应的相速度或sv波对应的相速度，θ表示p波对应的相角或sv波对应的相角。

公式(6)中的参数有如下关系：

其中，vp0和vsv0分别表示p波和s波垂直各向同性面的相速度，ε表示p波各向异性参数的度量，δ表示纵横波速度比的相关参数。可以利用公式(7)对公式(6)化简后，得到如下的相速度公式：

其中，vp(θ)表示以相角θ为自变量的p波相速度，vsv(θ)表示以相角θ为自变量的sv波相速度。根据相速度函数表达式(公式(8)和(9))，可以直接求导得到相速度对相角的偏导数：

通过上述具体实施方式可以得到以相角为自变量的p波相速度vp(θ)、p波相速度对相角的偏导数以相角为自变量的sv波相速度vsv(θ)及sv波相速度对相角的偏导数

实施例中，将求得的相应地代入ψp＝f(θp)和ψsv＝f(θsv)(由公式(1)得到)。可利用数值方法求得其函数与反函数θp＝f^-1(ψp)和θsv＝f^-1(ψsv)的数值表，其中，θp和θsv分别表示p波对应的相角和sv波对应的相角。

实施例中，按1度的间隔从0到90度取群角的角度值。根据群角的角度值查数值表θp＝f(ψp)和θsv＝f(ψsv)，可以得到相速度方向(相角)θp和θsv。将查表得到的相角θ(包括θp和θsv)代入christoff方程。

通过求解p波和sv波极化方向向量(u1,u3)来得到对应的p波极化角度υp和sv波极化角度υsv。极化角度的求解如下：

其中，υ表示p波极化角度或sv波极化角度。

其他实施例中，可以直接用相速度角度(p波对应的相角θp和sv波对应的相角θsv)作为中间变量，不求解极化角(p波对应的极化角υp和sv波对应的极化角υsv)与群角(p波对应的群角ψp和sv波对应的群角ψsv)的对应关系的显式，而是直接得到群角和极化角的对应关系，以此可以避免在计算过程中求解方程。即直接构造一个新的数值表，记作表二，可以直接用来查询：

υp＝f'(ψp),υsv＝f'(ψsv)(13)

上述各公式中，符号f、f、f'为函数符号，可以表示不同的函数关系。

通过上述具体实施过程可以得到群角和极化角的对应关系，具体地，可以是群角和极化角的数值对应表，或者可以是以群角为自变量且以极化角为因变量的函数关系。

根据上述对应关系，结合各向异性介质的对称横方向、沿对称轴方向的p波和s波速度和各向异性参数ε、δ，可以运行程序生成相应间隔的数值表，此数值表是一个四维表，按对群角和极化角关系影响大小排序为δ，ε，纵横波速度比群角ψ。经过我们研究分析纵横波速度比取成群角取为sinψ，插值时候误差小，因此需要选取的点就少，相应的数值表存储就小。

实施例中，可以在生成相应的数值表以后，对地震数据进行预处理，得到可用于偏移的炮记录。实施例中，可以确定炮的频谱，将其转换到时间域得到一个类似雷克子波的子波，并利用该子波进行正演模拟生成边界记录。实施例中，可以根据炮的总数，按例如每十炮为一组，分为若干组，并行读取炮记录的子波频率，生成正演模拟的震源子波。

为了降低弹性逆时偏移的带来内存的消耗和数据的读取带来计算量的巨大消耗的问题。对边界记录可以采用宽时间和宽空间记录，具体实现方式例如可以是：

若检波点波场的采用间隔为1ms，实施例中，可以在满足奈奎斯特采样定理的前提下，将炮点的正传波场按例如4ms的采样间隔来存储，反传的时候将检波点波场进行重采样的采样间隔例如为4ms，从而实现宽时间存储，如此一来，相当于把震源波场的正传存储量削减了4倍，因此可以大大降低存储量和读取量。实施例中，对于边界空间的存储我们可以按类似的方法实现宽空间存储，例如，若正常的正演模拟空间间隔为5m，那么可以在边界记录的时候选取10m空间间隔的边界点记录正传波场，这样可以进一步降低存储量。实施例中，在实现坡应廷矢量的求取、查表、波场分离、成像以及角道集的抽取时，可以选择将非规则网格上的点投影到宽空间的规则网格上来实现宽空间存储，以此可以减小相应的存储量和计算量。

实施例中，在弹性逆时偏移波场延拓过程中，每个空间点的坡印廷矢量可以通过如下公式计算得到：

其中，p为坡印廷矢量，τ为应力张量，υ为质点速度矢量，vx、vy、vz分别表示在x、y、z方向上的质点速度分量。

实施例中，为提高算法的稳定性，可从局部最小二乘准则出发，得到新的能流密度单位向量：

其中，表示单位向量，∑ω表示在以成像点为中心的邻域ω中相加求和。实施例中，利用群角的定义可以得到：

ψ＝arctan(nx/nz)(16)

其中，ψ表示群角，nx和nz分别表示单位向量在x方向和z方向上的分量。

实施例中，可根据公式(16)所得到的群角查表υp＝f'(ψp)υsv＝f'(ψsv)，再插值得到极化角，并通过下式完成各向异性波场分离：

uqp＝uxcosυ+uzsinυ，uqs＝-uxsinυ+uzcosυ(17)

其中，uqp和uqs分别表示伪p波和伪s波，ux表示p波或sv波的x轴分量，uz表示p波或sv波的z轴分量，υ表示p波的极化角或sv波极化角。

实施例中，一切与角度有关的运算都以单位向量间运算完成。因为求取角度需要反三角函数，该运算占用较多时钟周期，且只能取到0到π，另外的0到-π需借助其他的判断准则。单位向量的运算则不存在这些问题，相对方便和快捷，而且用单位向量来限制成像角度很方便。

实施例中，根据坡印廷矢量的方向传播特性，可以在成像条件上加衰减因子，从而可得到低频噪音的压制效果。实施例中，设定的角度阈值可为|α|＜π/4，以此可相当于主动进行低频噪音切除。因为前文中用到向量的运算，没必要直接求取反三角函数，而是直接用该角度的余弦值进行判断，即用cos2|α|＞0进行判断。与拉普拉斯算子相比，本实施例的方法更有优势，能够避免求导损坏波场信息。实施例中，类似地，可以在设定角度范围，例如|α|＜π/4范围内照明补偿，以此可以防止补偿过量。

实施例中，可以抽取各向异性介质中弹性逆时偏移共成像点角道集，以更好获得岩性信息。非转换波道集抽取相对简单，转换波道集抽取相对复杂。在二维情况下，可以使用一种简单方求取转换波的张角。实施例中，在求取群角时可以先将坡印廷矢量单位化，然后通过单位向量之间的运算得到群角，求出群角以后通过查表获得对应相角，以此可以避免因直接求取反三角函数的计算开销相对较大的问题。

实施例中，进行角道集的抽取。在完成波场分离步骤以后，在成像条件中，可以约定等式右端第一项为入射波场，第二项为反射波场。则pp波、ps波的角道集抽取方式分别为：

其中，ipp表示pp波的角道集，ips表示ps波的角道集，x表示x方向位移变量，θ表示相角，t表示时间变量，psrc(x,t)表示入射p波，prec(x,t)表示反射p波，srec(x,t)表示反射s波。mpp表示pp波的角度范围，mps表示ps波的角度范围，表示入射波单位矢量，表示反射波单位矢量。

实施例中，对公式(18)和(19)，关于入射波和界面法线的夹角求取以及转换波的极性校正，可以为：

其中，符号“.”表示向量间的点乘。

实施例中，对于转换波的极性校正，可以使用方法如下：

其中，符号“×”表示向量间的叉乘。在共成像点角道集抽取过程中，因为低频噪音存在于α的大角度，所以可以不对所有角度进行输出。实施例中，可以设定角度阈值为α＜π/4，以此能够主动进行低频噪音切除。实施例中，为节省计算内存及存储空间，根据炮点检波点格林函数互易定律，在进行完转换波极性校正后，输出pp、ps的单边道集即可满足要求。对于最终偏移，可以由上式(18)和(19)直接在共成像点角度域叠加得到：

其中，和分别表示单炮pp和ps波的叠加剖面。

图10是本发明一实施例中成像区域内沿对称轴方向的p波及s波速度和各向异性参数的构造示意图。如图10所示，图中数值分别指示所在地层的p波速度vp、s波速度vs、各向异性参数ε和各向异性参数δ。地层位置l1处，δ＝0.1322；地层位置l2处，ε＝0.2124，δ＝0.1966，vp＝1791.85m/s，vs＝2687.78m/s；地层位置l3处，ε＝0.2124，δ＝0.1966；地层位置l4处，δ＝0.1922；地层位置l5处，vp＝3444.68m/s，vs＝2296.46m/s；地层位置l6处，ε＝0.2686，δ＝0.1966；地层位置l7处，vp＝4958.49m/s，vs＝3305.66m/s，ε＝0.1562；地层位置l8处，ε＝0.2686；地层位置l9处，ε＝0.2124；地层位置l10和l11处，vp＝3444.68m/s，vs＝2296.46m/s；地层位置l12处，vp＝4201.59m/s，vs＝2801.06m/s。图11是本发明一实施例中正传波场快照示意图，图11显示炮点坐标xs＝5000m的炮记录在1.2s时刻正传波场的波场快照，其中，a)部分是垂直分量的波场快照，b)部分是水平分量的波场快照。图12是本发明一实施例中反传波场快照示意图。图12显示炮点坐标xs＝5000m的炮记录在1.2s时刻反传波场的波场快照，其中a)部分是垂直分量的波场快照，b)部分是水平分量的波场快照。图13是图12的反传波场的p波和s波分解结果示意图，其中，a)部分是分解得到的qp波的波场快照，b)部分是分解得到的qs波的波场快照。图14是图10地质构造的弹性波矢量成像的偏移成像剖面示意图，其中，a)部分是纵波成像剖面，b)部分是转换波成像剖面。图15是本发明一实施例中角道集示意图，图15显示水平位置标xs＝5000m处的偏移角道集，其中，a)部分是纵波成像的角道集，b)部分是转换波成像的角道集。根据图10～图15可以看出，根据本发明实施例的方法可以有效进行各向异性介质中弹性波矢量成像。

本发明基于角度域的各向异性波场分离方法进行初步尝试。其基本思想为：根据模型的参数范围，提前建立群角和极化角的数值表，在波场延拓过程中由坡印廷矢量方向求取能流密度方向，即群速度方向；求取每个空间点的群角，并确定该空间点位置的各向异性参数，通过数值表来查找指定参数间隔的提前计算好的极化角，由相邻点来插值该点极化角，完成在角度域的各向异性波场分解。因tti介质可视作vti介质经旋转得到，可只建立vti介质的数值表，然后，tti介质中的极化方向，由vti介质经角度旋转后，查数值表完成。

本发明实施例的各向异性介质中弹性波矢量成像的方法，可应用于提高各向异性波场分离计算效率和减小存储量，利用群角和极化角的高效数值表在角度域完成波场分离，能够解决传统拟导数算子的存储量大，多次计算傅里叶变换消耗计算时间等诸多问题，同时避免了直接由波动方程偏移进行速度建模所需要的多次求导计算，保持了波场相位和振幅信息。该方法在节省计算量的同时可查表获得相角来限制成像角度和照明补偿角度改善弹性逆时偏移的成像效果，便于后续的解释和处理工作，能够拓展和提升弹性波矢量成像方法的应用范围和效果。能够同时解决传统方法计算效率低和存储量大的问题。利用数值表中群角和相角的关系，可以直接限制成像角度和照明补偿角度，相当于拉普拉斯滤波。能够同时避免求导，保持波场信息，对提高成像的信噪比，改善成像效果有很大的益处。能确定和寻找地下生、储油构造，为油田的实际生产提供储量估计与确定井位。对地下复杂地区的油气、矿产资源勘探有重要应用价值。

基于与图1所示的弹性波矢量成像方法相同的发明构思，本申请实施例还提供了一种弹性波矢量成像装置，如下面实施例所述。由于该弹性波矢量成像装置解决问题的原理与弹性波矢量成像方法相似，因此该弹性波矢量成像装置的实施可以参见弹性波矢量成像方法的实施，重复之处不再赘述。

图16是本发明一实施例的弹性波矢量成像装置的结构示意图。如图16所示，本发明实施例的弹性波矢量成像装置，可包括：转换关系建立单元310、波场延拓单元320、波场分解单元330及弹性波矢量成像单元340，上述各单元顺序连接。

转换关系建立单元310，用于：利用成像区域中各点各向异性介质的设定方向的p波速度、设定方向的s波速度以及各向异性参数，建立群角与极化角之间的转换关系；

波场延拓单元320，用于：应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，得到震源正传波场和炮记录反传波场，以及相应的坡印廷矢量；

波场分解单元330，用于：按设定时刻间距，利用坡印廷矢量和群角与极化角之间的转换关系，对震源正传波场和炮记录反传波场进行波场分解；

弹性波矢量成像单元340，用于：利用多分量地震数据中所有炮对应的波场分解结果进行各向异性介质中弹性波矢量成像。

图17是本发明一实施例中转换关系建立单元的结构示意图。如图17所示，转换关系建立单元310，可包括：相速度获取模块311、群角与相角转换关系生成模块312及群角与极化角转换关系生成模块313，上述各模块顺序连接。

相速度获取模块311，用于：对christoff方程求特征值，并结合设定方向的p波速度、设定方向的s波速度以及各向异性参数，分别针对p波和sv波，计算得到以相角为自变量的相速度，其中，设定方向包括对称横方向和沿对称轴方向；

群角与相角转换关系生成模块312，用于：分别针对p波和sv波，将相速度和相速度对相角的偏导数代入群角与相角之间的关系式，计算得到群角与相角之间的转换关系；

群角与极化角转换关系生成模块313，用于：分别针对p波和sv波，在christoff方程中求解极化角，并结合各向异性参数、相速度及群角与相角之间的转换关系，计算得到群角与极化角之间的转换关系。

图18是本发明一实施例中群角与极化角转换关系生成模块的结构示意图。如图18所示，群角与极化角转换关系生成模块313，包括：群角间隔取值模块3131和群角和极化角数值对应表生成模块3132，二者相互连接。

群角间隔取值模块3131，用于：分别针对p波和sv波，在设定角度范围内按设定角度间隔取群角角度值；

群角和极化角数值对应表生成模块3132，用于：分别针对p波和sv波，在christoff方程中求解极化角，并结合各向异性参数、相速度、群角角度值及群角与相角之间的转换关系，计算得到群角和极化角之间的数值对应表，作为群角与极化角之间的转换关系。

图19是本发明一实施例中波场延拓单元的结构示意图。如图19所示，波场延拓单元320，可包括：炮记录获取模块321、边界记录获取模块322及波场延拓模块323，上述各模块顺序连接。

炮记录获取模块321，用于：对多分量地震数据进行预处理，得到多炮可用于逆时偏移的炮记录；

边界记录获取模块322，用于：由每炮的炮记录的频宽生成宽频子波，并利用宽频子波进行数值模拟，生成成像区域的边界记录；

波场延拓模块323，用于：针对每炮的边界记录和炮记录，利用弹性波方程数值解法进行弹性波场反传播，并在不同时刻将边界记录赋值给成像区域的边界位置，得到不同时刻的震源正传波场及相应的坡印廷矢量，在不同时刻将炮记录赋值给检波点位置，得到不同时刻的炮记录反传波场及相应的坡印廷矢量。

图20是本发明一实施例中波场分解单元的结构示意图。如图20所示，波场分解单元330，可包括：群角值生成模块331、极化角角度生成模块332及波场分解模块333，上述各模块顺序连接。

群角值生成模块331，用于：分别针对p波和sv波，利用设定时刻间距的坡印廷矢量计算得到相应的群角值；

极化角角度生成模块332，用于：分别针对p波和sv波，根据群角值和群角与极化角之间的转换关系，获取成像区域中各点的极化角角度；

波场分解模块333，用于：分别针对p波和sv波，将极化角角度代入波场分解公式，计算得到伪p波和伪s波，作为波场分解结果。

图21是本发明一实施例中极化角角度生成模块的结构示意图。如图21所示，极化角角度生成模块332，可包括：极化角值查找模块3321和极化角值插值模块3322，二者相互连接。

极化角值查找模块3321，用于：分别针对p波和sv波，利用群角值查找群角与极化角之间的转换关系的数值对应表，得到相应的极化角值；

极化角值插值模块3322，用于：分别针对p波和sv波，利用得到的极化角值进行插值，得到成像区域中各点的极化角角度。

图22是本发明一实施例中弹性波矢量成像单元的结构示意图。如图22所示，弹性波矢量成像单元340，可包括：单炮成像结果生成模块341、成像结果累加模块342及最终成像结果生成模块343。

单炮成像结果生成模块341，用于：在设定相角范围内，利用波场分解结果得到多分量地震数据中单炮的成像剖面、偏移角道集及照明强度；

成像结果累加模块342，用于：将多分量地震数据中所有炮的成像剖面、偏移角道集及照明强度分别累加，得到累加的成像剖面、累加的偏移角道集及累加的照明强度；

最终成像结果生成模块343，用于：用累加的照明强度除以累加的成像剖面得到各向异性介质中弹性波矢量成像的成像剖面，并由累加的偏移角道集得到各向异性介质中弹性波矢量成像的偏移角道集。

本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述各实施例所述方法的步骤。

本发明实施例还提供一种计算机设备，如图23所示，该计算机设备400可包括存储器410、处理器420及存储在存储器410上并可在处理420器上运行的计算机程序，所述处理器420执行所述程序时实现上述各实施例所述方法的步骤。

综上所述，本发明实施例的弹性波矢量成像方法、装置、存储介质及设备，不对接收的多分量资料进行波场分解，而是应用弹性波方程直接进行多分量资料延拓深度，并利用建立的群角与极化角之间的转换关系进行波场分解和成像，以此能够获得更准确的纵波及转换波界面反射信息。应用弹性波方程在成像区域进行多分量地震数据的波场深度延拓，能够精确考虑复杂速度、构造及各向异性，实现复杂地质条件下的准确成像。实施例中，基于查表的方法进行波场分解，如此一来，在处理复杂各向异性介质的地震数据时基本不会增加计算量，从而能够提高计算效率。实施例中，通过在设定相角范围内获取单炮的成像剖面，可以达到压制低频噪音的效果。通过在设定相角范围内获取照明强度，可以防止照明补偿过量。

在本说明书的描述中，参考术语“一个实施例”、“一个具体实施例”、“一些实施例”、“例如”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。各实施例中涉及的步骤顺序用于示意性说明本发明的实施，其中的步骤顺序不作限定，可根据需要作适当调整。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。