本发明属于地震数值模拟技术领域，具体地说，涉及一种提高地震转换波数值模拟精度的方法。
背景技术：
传统的基于弹性波波动方程数值模拟得到的弹性波波场既含有纵波(p波)也含有横波(s波)，对p波或者s波作处理的之前必须先对两个作波场分离，这一过程很难确保p波和s波的波场不出现畸变，也很难确保两者能完全分离。近些年来以四分量接收为特点的海底节点技术在地震领域逐渐兴起，尤其是转换波的应用越来越受到地球物理学者的关注，转换波的对于寻找裂缝性油气藏、估计岩性等有非常重要的意义。
当前弹性波模拟只能得到混合的弹性波波场，包含p波和s波，传统的梯度旋度算子分离纵横波很难保证分离之后的波场不存在振幅、相位的畸变。
技术实现要素：
有鉴于此，本发明提供了一种提高地震转换波数值模拟精度的方法。
为了解决上述技术问题，本发明公开了一种提高地震转换波数值模拟精度的方法，包括以下步骤：
步骤1、设置原始的弹性波波动方程；
步骤2、在步骤1得到的原始的弹性波波动方程引入横波分量ss和sp，得到修正后的弹性波波动方程；
步骤3、通过数学公式推导证明方程等价之后求解步骤2中的修正后的弹性波动方程，得到ss波场和sp波场的表示式，通过程序实现ss波场和sp波场的表示式，即转换波波场和纵波波场。
进一步地，步骤1中的原始的弹性波波动方程如式(1)所示，式(1)包括式(1-1)、式(1-2)和式(1-3)，
其中，x、y、z表示空间的三个方向，u、v、w表示弹性波场在x、y、z三个方向上的分量，t表示传播时间；vp表示纵波传播速度，vs表示横波传播速度。
进一步地，步骤2中的引入横波分量ss和sp，得到修正后的弹性波波动方程具体为：
引入变量sp，ss
sp＝(up，vp，wp)(2)
ss＝(us，vs，ws)(3)
其中，up、vp、wp分别表示p波在x、y、z三个方向上的分量，us、vs、ws分别表示s波在x、y、z三个方向上的分量；
将式(2)、(3)带入式(1)，得到修正后的弹性波波动方程，如方程组(4)所示，方程式(4)包括式(4-1)-式(4-9)：
u＝up+us(4-1)
v＝vp+vs(4-2)
w＝wp+ws(4-3)
进一步地，步骤3中的通过数学公式推导证明方程等价具体为：
将方程(4)中的第(4-4)式、第(4-7)式代入第(4-1)式中，得到:
式(5)与方程组(1)中的第(1-1)式完全等价；
同理，通过将方程组(4)中的第(4-5)式、第(4-8)式代入第(4-2)式得到的方程组与方程组(1)中的第2式完全等价：将方程组(4)中的第(4-6)式、第(4-9)式代入第(4-2)式得到的方程组与方程组(1)中的第(1-1)式完全等价，即方程组(4)与方程(1)等价；
接着，证明引入的变量sp与ss能够分别代表纵波波场和横波波场：
因为sp旋度的二阶导数为0，推出sp旋度的一阶导数为常数，sp的旋度为线性函数或者为常数，根据波动方程的初始零解的条件得到，sp的旋度为0，即sp是无旋场，sp是纵波；同理，证明ss是横波。
进一步地，步骤3中的求解修正后的弹性波动方程，得到ss波场和sp波场具体为：通过有限差分方法求解方程式(4-1)-(4-9)，得到p波和s波的表示式，空间6阶：
其中，ui表示p波，uj表示s波，dh表示网格长度；通过程序实现p波和s波的表示式，得到转换波横波场和纵波波场。
与现有技术相比，本发明可以获得包括以下技术效果：
1)本发明获得完全分离的p、s波波场，波场干净，不存在纵横波相互干扰。
2)本发明分离之后的p、s波振幅、相位保持不变，有利于进一步的p、s传播特性研究。
3)本发明基于gpu加速，提升模拟效率。
当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1是本发明纵波速度模型图；
图2是本发明横波速度模型图；
图3是本发明0.5s时刻纵波波场快照图；
图4是本发明0.5s时刻转换波波场快照图；
图5是本发明p波单炮记录图；
图6是本发明s波单炮记录图。
具体实施方式
以下将配合实施例来详细说明本发明的实施方式，藉此对本发明如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。
本发明的思想是：本发明提出了等价的弹性波方程，该方程以p波和s波为变量，通过求解该波场方程得到完全分离的转换波波场，同时得到纵波波场。该方法在弹性波理论研究、转换波场特征分析以及转换波处理解释等方面能发挥重要的作用。
该转换波数值模拟的理论基础是通过变量的替换，在原先的弹性波波动方程中引入横波分量ss和sp,通过数学公式推导证明方程等价之后求解修正后的弹性波动方程，得到ss波场和sp波场，即转换波波场和纵波波场。具体数学推导过程如下：
步骤1、设置原始的弹性波波动方程：
其中，x、y、z表示空间的三个方向，u、v、w表示弹性波场在x、y、z三个方向上的分量，t表示传播时间；vp表示纵波传播速度，vs表示横波传播速度。
式(1)表示原始的弹性波波动方程，包括式(1-1)、式(1-2)和式(1-3)；
设置纵波速度模型，速度模型中的速度代表vp，如图1所示，该速度模型分为3层，从上向下，第一层为海水，速度为1500m/s，第二层的速度为2000m/s，第三层的速度为2500m/s；检波点在表面，炮点在表面，模拟实际地震数据采集模式；
设置横波速度模型，速度模型中的速度代表vs，如图2所示，每一层横波速度通过地质经验公式通过纵波速度转换得到，该速度模型分为3层，从上向下，第一层的速度为866m/s，第二层的速度为1155m/s，第三层的速度为1443m/s；检波点在表面，炮点在表面；
步骤2、引入横波分量ss和sp，得到修正后的弹性波波动方程(4)：
sp＝(up，vp，wp)(2)
ss＝(us，vs，ws)(3)
其中，up、vp、wp分别表示p波在x、y、z三个方向上的分量，us、vs、ws分别表示s波在x、y、z三个方向上的分量。
将式(2)、(3)带入式(1)，得到如下方程组(4)，方程式(4)包括式(4-1)-式(4-9)：
u＝up+us(4-1)
v＝vp+vs(4-2)
w＝wp+ws(4-3)
其中，up、vp、wp分别表示p波在x、y、z三个方向上的分量，us、vs、ws分别表示s波在x、y、z三个方向上的分量；x、y、z表示空间的三个方向，u、v、w表示弹性波场在x、y、z三个方向上的分量，t表示传播时间；vp表示纵波传播速度，vs表示横波传播速度。
步骤3、通过数学公式推导证明方程等价之后求解修正后的弹性波动方程，得到ss波场和sp波场，即转换波波场和纵波波场：
首先，证明方程(1)和方程式(4)等价。
将方程(4)中的第(4-4)式、第(4-7)式代入第(4-1)式中，得到:
其中，x、y、z表示空间的三个方向，u、v、w表示弹性波场在x、y、z三个方向上的分量，t表示传播时间。vp表示纵波传播速度，vs表示横波传播速度。up、us分别表示x分量上的纵波波场分离和横波波场分量；
容易发现，式(5)与方程组(1)中的第(1-1)式完全等价。
同理，通过将方程组(4)中的第(4-5)式、第(4-8)式代入第(4-2)式得到的方程组与方程组(1)中的第(1-2)式完全等价：将方程组(4)中的第(4-6)式、第(4-9)式代入第(4-2)式得到的方程组与方程组(1)中的第(1-1)式完全等价，即方程组(4)与方程(1)等价。
接着，证明引入的变量sp与ss能够分别代表纵波波场和横波波场：
其中，i，j，k，分别代表三维空间网格点x,y,z方向上的坐标，
因为sp旋度的二阶导数为0，可以推出sp旋度的一阶导数为常数，sp的旋度为线性函数或者为常数，根据波动方程的初始零解的条件可以得到，sp的旋度为0，即sp是无旋场，sp是纵波。同理，也可以证明ss是横波。
实施例1以二维情况为例，用有限差分法求解。
一种提高地震转换波数值模拟精度的方法，按照以下步骤进行：
步骤1、推导二维情况的完全分离的转换波数值模拟方程，该方程为：
u=up+us(4-1)
w＝wp+ws(4-3)
步骤2：通过有限差分方法求解该方程，得到p波和s波的表示式(空间6阶)：
其中，i,j表示参与计算的空间网格点的x、y方向上的坐标，与上下文中的i，j对应，a(n)表示差分系数；ui表示p波，uj表示s波，dh表示参与计算的相邻网格的距离。
步骤3、通过程序实现步骤2中的离散表达式，得到转换波横波场和纵波波场，具体为：
步骤3.1、对于二维各项同性介质中的弹性波波动方程
其中，x、y、z表示空间的三个方向，λ、μ表示拉梅常数，τzz、τxx、τxz表示应力的分量，vx、vz表示速度分量的分量。
将方程离散，将其展开，以
时间二阶近似得到
vx、vz表示速度分量的分量，δt表示采样率，m表示泰勒展开后取的阶数。
两式相减，得到2m阶的近似差分方程
vx表示速度分量的分量，δt表示采样率。
上述的n对等式中，对第n对等式的等号两边同时乘以系数cn(n)后相减，再将所有式求和，得到：
令上式右端仅保留
将上式代入(16)得
代入(15)得
同理可以得到所有的应力-位移方程的离散格式：
其中，i，j，k分别代表三维空间网格点x，y，z方向上的坐标，x、y、z表示空间的三个方向，λμ表示拉梅常数，τzz、τxx、τxz表示应力的分量，vx、vz表示速度分量的分量，δt表示采样率，n表示差分阶数，δx、δz表示横向、纵向的网格步长，cn(n)为系数。
至此，得到弹性波的模拟结果。
步骤3.2、通过公式(7)、(8)，从弹性波模拟得到的vx、vz分量提取纵横波，公式(7)、(8)的推导过程见上文。
其中，i,j表示某一网格点的坐标，a(n)表示差分系ui表示p波，uj表示s波，dh表示网格长度，也就是前面公式提到的δx、δz。
目前市场的转换波提取方法都是通过弹性波波场应用散度旋度算子分离得到纵横波波场，这种方法很难不使波场的振幅、相位信息畸变，对后期的数据研究不利，对此，本发明提出一种完全分离的纵横波数据模拟方法。相对比与传统的散度旋度算子分离纵横波，本发明方法得到的纵横波彼此完全分离，互不干扰，并且得到的纵横波波场保留原始的振幅、相位信息，对转换波特性研究、转换波数据处理、反演方法研究都有很大的益处。
图3是纵波的波场快照，图4是转换横波的波场快照图，由具体实施方案的步骤2、步骤3得到，通过两张图的对比，可以发现纵波和横波彼此完全分离，互不干扰，验证了本发明的有效性。
图5的纵波记录，图6是横波记录，由具体实施方案的步骤2、步骤3得到，图3、图4是通过波场快照对比，图5、图6通过最后得到的记录对比，p波和s波彼此完全分离，在领偏移距处不产生转换波、转换波在炮点两侧极性相反，这都与与实际现象相同，再一次证明了本发明的有效性。
由图3-图6可知，两种波场都是彼此完全独立的波场，纵波记录上不存在横波，横波记录上也不存在纵波。单炮记录一般是正演模拟的最后结果，从单炮上，纵横波也是完全独立的波场，不会存在传统的纵横波分离不干净导致纵波波场上有横波或者横波记录上有纵波的现象。
上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。