成像方法、成像装置以及计算机存储介质与流程
本申请涉及地震勘探领域，尤其是涉及一种成像方法、成像装置以及计算机存储介质。
背景技术：
弹性波地震勘探技术具有介质参数求取精度高、介质参数反演的确定性高、成像精度高、储层预测精度高和流体识别精度高等优势。由于弹性波方程包含更多信息，相较声波方程而言，其能够更加准确地描述地球介质的组成。因此，研究弹性波逆时叠前深度偏移技术对于改进地震成像质量的研究，以及解决一些地震勘探的新问题是至关重要的。
其中弹性波成像结果的耦合噪音，低频噪音和成像假象噪音的问题已经成为国内外地球物理勘探学界关注的热点与难点。
弹性波成像噪音主要来自于耦合噪音，低频噪音和成像假象噪音，如何压制这些噪音，多年来国内外学者也提出了多种技术手段实现噪音的去除，诸如解耦传播，上下方向波分解等技术方法，虽取得了一定的效果，但是弹性波成像质量改进不多。
zhang等(2007)推导出一种解耦的弹性波一阶速度应力方程，实现纵横波的分离。xiao和leaney(2010)引入了纵波应力，提出一种基于弹性波一阶速度应力方程的解耦方程组，该方程组在数学上与zhang等(2007)提出的解耦方程组是等价的，他们实现了纵横波解耦传播并对vsp数据进行成像。zhou等(2016)基于zhang等(2007)的工作，将解耦方程组由2维拓展到3维。fei等(2010)提出可以将上行和下行波场分离，然后对不同传播方向的波场分别应用成像条件，组合成像。
fei等(2015)提出成像假象的概念，他们认为liu(2011)的做法虽然能够较好地压制低频噪音，但是还是忽略了当速度变化剧烈时，后向散射造成的成像假象，即由炮点前向波场和检波点后向波场互相关成像造成的假象，称之为一次波假象(primaryfalseimage)，并由此提出了de-primaryrtm算法。王一博等(2016)在fei等(2015)工作的基础上，推导出了标量波方程上下左右4个象限方向波的分离。
以上方法在压制低频噪音，耦合噪音和成像假象噪音上有一定的效果，但解耦方程分离纵横波只能适用于各向同性介质，无法在各向异性介质中应用，上下方向波分解可以较好地压制由同方向传播的方向波相关产生的低频噪音，并且一定程度上压制由炮点前向波场和检波点后向波场互相关造成的假象噪音，但如果地下反射层倾角较大，并且目标位置分布复杂时，仍会产生较多的来自左右方向波的成像噪音，因此现存方法都存在一定的技术限制。
技术实现要素：
本发明的目的在于提供了一种成像方法、成像装置以及计算机存储介质，以解决上述的成像质量较差的问题。
为了实现上述目的，本发明提供了一种成像方法、成像装置以及计算机存储介质，该成像方法包括：
构建地下炮点正传波场以及检波点反传波场；
判断传播所述正传波场以及所述反传波场的介质类型；
当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数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插值算法，对所述正传波场和或反传波场进行纵横波场分离；
对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；
从所述多个方向波波场中选择至少两个方向波波场，并对选择出的方向波波场进行归一化互相关矢量成像。
该方法利用弹性波矢量波场的分离技术，即在各向同性介质中采用解耦方程法，在非均匀各向异性介质中采用基于变异函数的改进idw插值算法，以实现纵横波的分离，之后再应用傅里叶变换方法，实现纵横波全方位方向波的分解，在实现纵横波分离和全方位方向波分解的基础上，从多个方向波波场中选择至少两个方向波波场，并对选择出的方向波波场进行矢量成像，进而解决了现有技术中只能在标量波上进行上下左右方向波的分解并进行成像，而不能应用到矢量波(弹性波)上的问题。同时，该成像方法可以有效地压制弹性波逆时偏移过程中产生的耦合噪音，低频噪音以及成像假象噪音等，提高成像质量。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本申请的实施例提供的一种成像方法的流程图；
图2为本申请的实施例提供的一种成像方法中采用变异函数idw差值算法将获取到的所述弹性波分离为纵波和横波的流程图；
图3为本申请的实施例提供的一种成像方法中采用变异函数idw差值算法过程中如何得到自褶组合窗函数的流程图；
图4为本申请的实施例提供的一种成像方法中将纵波分离为上行纵波和下行纵波的流程图；
图5为本申请的实施例提供的一种成像方法中上行横波分离为左上行横波和右上行横波的流程图；
图6为本申请的实施例提供的一种成像方法中下行横波分离为左下行横波和右下行横波的流程图；
图7为本申请的实施例提供的一种成像方法中横波分离为上行横波和下行横波的流程图；
图8为本申请的实施例提供的一种成像方法中上行横波分离为左上行横波和右上行横波的流程图；
图9为本申请的实施例提供的一种成像方法中下行横波分离为左下行横波和右下行横波的流程图；
图10(a)为marmousi2模型弹性波逆时偏移局部成像结果；
图10(b)为marmousi2模型弹性波逆时偏移局部成像结果；
图11(a)为marmousi2模型的pp成像结果；
图11(b)为盐丘模型的pp成像结果；
图11(c)为盐丘模型的ps成像结果；
图12a为marmousi2模型各向异性参数分布及选取的参考模型；
图12b为marmousi2模型各向异性参数分布及选取的参考模型；
图12c为marmousi2模型各向异性参数分布及选取的参考模型；
图12d为marmousi2模型各向异性参数分布及选取的参考模型；
图13为本申请的实施例中变异函数idw算法插值运算过程中选取的参考模型中的局部极值；
图14为本申请的实施例中变异函数idw算法插值运算过程中选取的参考模型；
图15为本申请的实施例提供的弹性波逆时偏移成像装置的结构示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。
在逆时偏移成像领域，研究人员为了避免逆偏移成像过程中受到耦合噪音，低频噪音和成像假象噪音等的干扰做了许多努力。例如，fei等(2015)提出一种de-primaryrtm算法，基于该算法对弹性波进行上下方向波的分解，然后将分解后的波中的噪声去除。但其将上下行波的分离隐含在成像条件里，虽然可以实现上下行波的分离，但还是无法得到分离的波场。
wang等(2016)提出将纵横波分离与上下方向波分解相结合，可以有效压制耦合噪音，低频噪音和成像假象噪音，并显式分离了上下方向波；但在各向异性异性介质中使用zhang和mcmechan(2010)的算法，计算量过大。该方法不适用于三维情况(三维情况下要考虑该介质为各向异性异性)，因此其在地质勘探这种需要考虑三维情况下也应用受限。另外，wang等(2016)的工作没有考虑左右方向波的分解，只是进行了上下两个方向的方向波分解，对于存在垂直构造，其方案存在一定的局限。
王一博等(2016)在fei等(2015)工作的基础上，推导出了标量波方程上下左右4个象限方向波的分离，但该方法适用于标量波(声波)，不适用于矢量波场(弹性波)，因此也无法运用到弹性波逆时偏移成像上。
综合上述内容可知，不同的研究人员已经从不同的方面来试图解决成像过程中的噪声问题，但各种方法的效果不佳或者无法在特定的条件下使用。
为了解决成像过程中的噪声问题，本申请提供了一种成像方法。图1为本发明的实施例提供的成像方法的流程示意图，参照图1所示，可以包括以下步骤：
s101：构建地下炮点正传波场以及检波点反传波场；
s102：判断传播正传波场以及反传波场的介质类型；
s103：当介质为各向同性介质时，根据预设的方程解耦算法，对正传波场和或反传波场进行纵横波场分离；当介质为非均匀各向异性时，根据预设的变异函数idw插值算法，对正传波场和或反传波场进行纵横波场分离；
s104：对分离出的纵波波场和横波波场在上下左右四个方向进行分解，得到所述纵波波场和/或所述横波波场的多个方向波波场；
s105：从所述多个方向波波场中选择至少两个方向波波场，并对选择出的方向波波场进行矢量成像。
其中，在各向同性介质中应用解耦方程法以及在非均匀各向异性介质中应用基于变异函数的idw插值算法对正传/反传波场进行纵横波波场的分离。
本发明的实施例提供了一种成像方法。该方法利用弹性波矢量波场的分离技术，即在各向同性介质中采用解耦方程法，在非均匀各向异性介质中采用基于变异函数的改进idw插值算法，以实现纵横波的分离；之后应用复数域拓展方法，实现纵横波全方位方向波的分解，在实现纵横波分离和全方位方向波分解之后，从多个方向波波场中选择至少两个方向波波场，并对选择出的方向波波场进行矢量成像，进而解决了现有技术中只能在标量波上进行纵横波以及上下左右方向波的分离而不能运用于矢量波(弹性波)上的问题。同时，该成像方法可以有效地压制弹性波逆时偏移过程中产生的耦合噪音，低频噪音以及成像假象噪音等，提高成像质量。
具体来说，idw(inversedistanceweighted)是一种常用而简便的空间插值方法。它以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样本点赋予的权重越大。例如，设平面上分布一系列离散点，已知其坐标和值为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进行约束，选取在所述临界值范围内出现概率最大的参考点，根据所述参考点搜索得到上述参考模型(参看图13和图14所示)。
图3是本发明一实施例中获取自褶积组合窗函数的方法的流程示意图。如图3所示，可包括以下步骤：
s301：选择主瓣和旁瓣性能高于设定阈值的窗函数作为原始窗函数；
s302：对原始窗函数做l次自褶积运算得到自褶积后的窗函数，其中，l为正整数；
s303：对自褶积后的窗函数与原始窗函数进行加权运算，得到所述自褶积组合窗函数。
在上述实施例中，对选定的原始窗函数进行自褶积运算，自褶积运算使得原始窗函数的主瓣性能变差，旁瓣性能变好，然后根据当前的需要对原始窗函数和自褶积后的窗函数进行加权运算，从而得到一个满足要求的窗函数。利用该种自褶积组合窗函数的到的分离算子可以提高波场分离的精度和稳定性。
在上述步骤s302中，进行加权运算的加权系数可以具备特点：在确定需要主瓣性能优先的情况下，原始窗函数的加权系数高于自褶积后的窗函数的加权系数；在确定需要旁瓣性能优先的情况下，自褶积后的函数的加权系数高于原始窗函数的加权系数。以此在维持适当的主瓣宽度的前提下，可以尽可能增大旁瓣衰减。
为使本领域的技术人员更容易理解本发明的技术方案，下面将以一具体实施例说明变异函数idw差值算法以及该步骤的子步骤。
对于非均匀各向异性而言，理论上每个采样点都需要根据相应的各向异性参数计算偏振向量，才能获得更好的分离结果，但是这样处理计算量太大，因此，可以从初始模型中选择出n个参考模型，分别在波数域进行矢量波场(弹性波)分离，然后反变换回空间域，根据事先计算好的权重系数，重构初始模型的分离结果。
因此，考虑到引入描述区域化变量的空间结构性变化和随机性变化的变异函数，代替两点之间的距离计算权重系数。具体的，可以按照以下公式计算初始模型各向异性参数分布的变异函数值：
其中，μ(h)表示变异函数值，φ＝f[δ+2(ε-δ)sin2(α-θ)]sin2(α-θ)，
在计算得到多个分散的变异函数值后，根据这多个变异函数值进行拟合，得到理论变异函数模型。
然后，按照以下公式计算各个参考模型的权重系数：
其中，wk表示权重系数，(εk,δk,θk)表示插值参考点。
上述
下面结合一个具体的实施例进行说明，然而值得注意的是，该具体实施例仅是为了更好地说明本发明，并不构成对本发明的不当限定。
基于偏振特性对矢量波场进行分离的过程中所存在的一个不可避免的问题是如何在保证计算效率的前提下，让波场更准确地投影到相应的偏振方向上。在波数域，称之为投影，而在空间域，称之为空间滤波。
在各向同性介质中，矢量波场的分离效果仅与微分算子的精度有关，精度越高，分离效果越好，同时计算量也增大。然而，在各向异性介质中，矢量波场的分离效果主要与各向异性参数和拟微分算子的精度这两个因素有关，考虑到在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(α-θ)]sin2(α-θ)(12)
其中，α表示相角，传播方向与z轴的夹角，代表的是各向同性部分，θ表示tti介质对称轴的倾角，
去除上述近似公式中的各向同性部分，只保留线性各向异性部分，可以得到：
φ＝f[δ+2(ε-δ)sin2(α-θ)]sin2(α-θ)(13)
其中，φ和ε、δ成近似线性关系，和sin2(α-θ)成近似比例关系。
因此，可以设初始模型为m＝{ε，δ，θ}，在初始模型条件下偏振角的各向异性部分值为φ，进一步的，可以条件选择n个参考模型mk，并计算得到每个参考模型对应的φk，并通过以下公式插值得到φ：
其中，wk表示权重系数。
idw插值算法在计算权重时，仅考虑插值点与参考点的空间位置关系，也就是当不同参考点与插值点空间位置一致时，就有相同的权重系数，这种插值算法虽然执行效率较高，但是没有考虑参考点之间的差异。
在本例中，对权重系数的确定进行了改进，具体包括：假设φ满足二阶平稳假设或者本征假设，使用变异函数μ(h)来描述φ随着空间位置不同而变化的特性，按照以下公式计算得到变异函数值：
其中，
引入变异函数r后可以得到权重系数wk：
因为μ(h)是距离的函数，因此将μ(h)引入idw插值算法的权重计算中，可以获得更为可靠的插值结果。
对于插值参考点(εk,δk,θk)的选择，插值参考点的选择方法可包括：根据所述初始模型计算得到多个分散的变异函数值；对得到的多个分散的变异函数值进行拟合，得到理论变异函数模型；根据所述变异函数模型确定变异函数中的各个参数的参数值；遍历待数值模拟初始模型中不同数值的ε、δ、θ出现的概率，选取在拟合的变异函数临界值之内的，出现概率最大的点作为插值参考点。
当然的，除了采用上述纵横波分离算法外，本发明的其它实施方式还提供其它类型的分离方法。例如，可以根据速度差异分离波场或者采用波动方程法等。因此，本申请并不对纵横波分离的算法进行限制，只要能将耦合纵横波分离的方法，都是符合本申请的要求。
参照图1所示，本申请的方法的步骤s104：经过纵横波分离之后的纵横波波场进行上下方向的方向波分解中，参照图4所示，所述经过纵横波分离之后的纵波波场进行上下方向的方向波分解，包括以下几个子步骤：
s401：将分离出的纵波波场进行时间复数域拓展；
其中，复数域拓展策略是指将数据转化为复数的形式，进行各种计算的方法，该策略基于hilbert变换，可以有效地提高计算效率。由频谱分析可知，经过hilbert变换后，波场在频率域的振幅保持不变，相位倒转。对于经过复数域拓展的波场，可以发现在频率域，其振幅只有在正频率处有值，而且其相位保持不变。可以利用这个性质，将复数域拓展运用于全方位方向波分解，有效较低计算量。
s402：在z方向对时间复数域拓展之后的纵波波场进行傅里叶变换，得到波数域纵波波场；
其中，kz表示纵横波场在z方向上的波数。当地震波在地下介质传播时，定义向下为正，因此向下传播的方向波其kz为正，向上传播的方向波，其kz为负。
s403：从所述波数域纵波波场中选择波数大于零的部分，波数小于零的部分置为零，进行z方向上的反傅里叶变换，得到纵波波场的下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行z方向上的反傅里叶变换，得到所述纵波波场的上行方向波波场。
参照图1所示，本申请的成像方法还包括步骤：
纵横波的上行方向波波场进行左右方向的方向波分解，分解为右上行方向波波场和左上行方向波波场；
在该步骤中，参照图5所示，所述纵波的上行方向波波场进行左右方向的方向波分解，分解为右上行方向波波场和左上行方向波波场，包括以下几个子步骤：
s501：将分离出的所述纵波波场的上行方向波波场进行时间复数域拓展；
其中，kx表示纵波纵横波场在x方向上的波数，定义向右为正，因此向右传播的方向波其kx为正，向左传播的方向波，其kx为负；
s502：在x方向对时间复数域拓展之后的所述纵波波场的上行方向波波场进行傅里叶变换，得到波数域纵波波场；
s503：从所述波数域所述纵波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行x方向上的反傅里叶变换，得到纵波波场的右上行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行x方向上的反傅里叶变换，得到所述纵波波场的左上行方向波波场。
参照图6所示，所述纵波的下行方向波波场进行左右方向的方向波分解，分解为右下行方向波波场和左下行方向波波场，包括以下几个子步骤：
s601：将分离出的所述纵波波场的下行方向波波场进行时间复数域拓展；
s602：在x方向对时间复数域拓展之后的所述纵波波场的下行方向波波场进行傅里叶变换，得到波数域纵波波场；
s603：从所述波数域所述纵波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行x方向上的反傅里叶变换，得到纵波波场的右下行方向波波场；并从所述波数域纵波波场中选择波数小于零的部分，波数大于零的部分置为零，进行x方向上的反傅里叶变换，得到所述纵波波场的左下行方向波波场。
参照图7所示，所述横波信号能在上下方向上被分离为上行横波信号和下行横波信号包括以下几个子步骤：
s701：将分离出的横波波场进行时间复数域拓展；
s702：在z方向对时间复数域拓展之后的横波波场进行傅里叶变换，得到波数域横波波场；
s703：从所述波数域横波波场中选择波数大于零的部分，波数小于零的部分置为零，进行z方向上的反傅里叶变换，得到横波波场的下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行z方向上的反傅里叶变换，得到所述横波波场的上行方向波波场。
该几个步骤中的复数域、kz以及反傅里叶变换等已经在上文中进行了解释，在此不再赘述。
其中，在该步骤s104中，参照图8所示，所述上行横波信号能在左右方向上分离为左上行横波信号和右上行横波信号包括以下几个子步骤：
s801：将分离出的所述横波波场的上行方向波波场进行时间复数域拓展；
s802：在x方向对时间复数域拓展之后的所述横波波场的上行方向波波场进行傅里叶变换，得到波数域横波波场；
s803：从所述波数域所述横波波场的上行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行x方向上的反傅里叶变换，得到横波波场的右上行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行x方向上的反傅里叶变换，得到所述横波波场的左上行方向波波场。
该几个步骤中的复数域、kx以及反傅里叶变换等也已经在上文中进行了解释，在此不再赘述。
参照图9所示，该步骤s104中，所述下行横波信号能在左右方向上分离为左下行横波信号和右下行横波信号包括以下几个子步骤：
s901：将分离出的所述横波波场的下行方向波波场进行时间复数域拓展；
s1072：在x方向对时间复数域拓展之后的所述横波波场的下行方向波波场进行傅里叶变换，得到波数域横波波场；
s903：从所述波数域所述横波波场的下行方向波波场中选择波数大于零的部分，波数小于零的部分置为零，进行x方向上的反傅里叶变换，得到横波波场的右下行方向波波场；并从所述波数域横波波场中选择波数小于零的部分，波数大于零的部分置为零，进行x方向上的反傅里叶变换，得到所述横波波场的左下行方向波波场。
参照图1所示，步骤s105为：根据实际地质情况，选择特定组合的方向波波场，应用角度约束的矢量成像条件进行成像，得到地下成像结果；在该步骤中，弹性波根据纵横波的传播方向以及上下左右四个方向可以将弹性波分为8个方向波波场。以纵波为例，分离后的各波场记为
其中，s、r分别表示炮点处弹性波和检波点处弹性波。上标p代表纵波。下标l,r代表左右，u,d代表上下。炮点处和检波点处的纵波按照上下左右四个方向分解后可以表示为：
同理，炮点处和检波点处的横波按照上下左右四个方向分解后可以表示为：
以纵波为例，炮点处p波上下左右四个象限的方向波可以通过以下公式进行表示：
其它位置以及其他类似的波场可以按照该公式的形式进行表示，在此不在赘述。
同时，为了降低计算量和存储量，可以将上述的数据拓展到复数域。对于拓展到复数域的数据，ω＜0时值为零，因此可以不需要对数据进行时间fourier变换，降低计算量，将公式(19)—公式(22)简化为：
在对炮点正传以及检波点反传的弹性波场进行纵横波解耦以及上下左右方向的方向波分解以及简化后，可以得到分离出来的16个方向波波场。选择方向波组合进行互相关成像的结果要多于现有技术中只进行纵横波解耦以及上下方向波的分解得到8个方向波波场进行组合互相关成像的结果。
因此，采用本申请的方法可以具有较多的成像结果的选择。同时，当成像结果相对较多时，可以更容易地将成像结果中的噪声结果进行排除，也可以根据地质靶区的具体地质构造，排除一些特定的方向波组合。例如，ildld,ilulu,irdrd,iruru，这些相同下标的组合是产生成像低频噪声的主要来源；iluld,ilurd,iruld,irurd这些炮点波场上行波与检波点波场下行波，当存在强速度变化时，由于后向散射的存在，会导致一次波假象。其中，i表示成像结果。
对于水平构造，理论上，通过纵横波分离以及上下方位波分解就可以得到较好的成像结果，但是实际地质状况复杂，没有完全的水平或者垂直构造，仅仅进行上下方向波或者左右方向波是不够的，只有将纵横波分离与全方位方位波的分解结合起来，才可以压制特定的噪音，得到较好的成像结果。
当将存在噪声的成像结果排除后，可以将剩余的波形根据预设的方向波组合成像。例如，预设的方向波组合可以包括：炮点处的左下行纵波、检波点处的右上行纵波以及检波点处右上行横波。采用该组合可以得到如下的成像条件：
其中，m代表空间位置，θ代表adcigs中的离散角度，θp代表实际计算的角度，σ代表高斯函数的方差，
上述的方向波组合的成像效果相对较优，可以作为最终的成像结果。
当然的，除了上述的方向波组合外，还可以选择其他方向波组合进行逆时偏移成像。勘探人员可以根据实际的地形进行选择多组组合成像，再将各种方向波组合的成像结果进行比较，然后选择成像质量较高的作为结果。
虽然上文描述的过程流程包括以特定顺序出现的多个操作，但是，应当清楚了解，这些过程可以包括更多或更少的操作，这些操作可以顺序执行或并行执行。
实施例1：
参照图10(a)到图11(c)所示，为了验证本申请的成像方法的成像质量，实验人员利用该方法在marmousi2模型以及盐丘模型上成像。具体参看图10(a)和图10(b)所示，该图10(a)以及图10(b)为在marmousi2模型上pp成像以及ps成像的局部结果。其中，pp表示炮点处纵波和检波点处的纵波。ps表示炮点处纵波以及检波点处横波。
图11(a)到(c)分别为该方法在marmousi2模型上的pp成像结果；在盐丘模型上的pp成像结果以及在盐丘模型上的ps成像结果。其中pp表示炮点处纵波和检波点处的纵波。ps表示炮点处纵波以及检波点处横波。
通过上述的各图可以看出，将弹性波进行纵横波的分离的基础上又将其进行上下左右方向的分离，进而可以将弹性波分解为多个方向波波场，特定的方向波组合成像，将形成各种低频噪音以及假象噪音，将这些特定的方向波组合去除，选择合适的方向波组合进行互相关成像，使得采用本发明的成像方法的成像质量较优。
基于与图1所示的成像方法相同的发明构思，本申请还提供了一种弹性波成像装置。参照图15所示，该成像装置可以包括：构建单元401，用于构建地下炮点正传波场以及检波点反传波场；分离单元402，用于判断传播所述正传波场以及所述反传波场的介质类型；当所述介质为各向同性介质时，根据预设的方程解耦算法，对所述正传波场和或所述反传波场进行纵横波场分离，以及当所述介质为非均匀各向异性时，根据预设的变异函数idw插值算法，对所述正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和或横波波场在上下左右四个方向进行分解，得到所述纵波波场和或所述横波波场的多个方向波波场；成像单元403，用于从所述多个方向波波场中选择至少两个方向波波场，并对选择出的方向波波场进行归一化互相关矢量成像。
基于与图1所示的成像方法相同的发明构思，本申请还提供了一种计算机存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：获取目标地层的地下炮点正传波场以及检波点反传波场；根据预设的方程解耦算法，对在所述目标地层的各向同性介质中传播的正传波场和或反传波场进行纵横波场分离；并根据预设的变异函数idw插值算法，对在所述目标地层的非均匀各向异性介质中传播的正传波场和或反传波场进行纵横波场分离；对分离出的纵波波场和横波波场在上下左右四个方向进行分解，得到所述纵波波场和/或所述横波波场的多个方向波波场；从所述多个方向波波场中选择至少两个方向波波场，并对选择出的方向波波场进行矢量成像。
为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本申请时可以把各单元的功能在同一个或多个软件和/或硬件中实现。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
在一个典型的配置中，计算设备包括一个或多个处理器(cpu)、输入/输出接口、网络接口和内存。
内存可能包括计算机可读介质中的非永久性存储器，随机存取存储器(ram)和/或非易失性内存等形式，如只读存储器(rom)或闪存(flashram)。内存是计算机可读介质的示例。
计算机可读介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机的存储介质的例子包括，但不限于相变内存(pram)、静态随机存取存储器(sram)、动态随机存取存储器(dram)、其他类型的随机存取存储器(ram)、只读存储器(rom)、电可擦除可编程只读存储器(eeprom)、快闪记忆体或其他内存技术、只读光盘只读存储器(cd-rom)、数字多功能光盘(dvd)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读介质不包括暂存电脑可读媒体(transitorymedia)，如调制的数据波场和载波。
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。
本领域技术人员应明白，本申请的实施例可提供为方法、系统或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。
- 还没有人留言评论。精彩留言会获得点赞！