为了克服上述现有技术的不足，本发明的目的在于提出一种地震纵、横波波场分离与去噪的方法，其依赖地震信号之间的时距规律差异，通过两次正常时差校正处理的手段，分别把P-P波和P-S波校平，使之在横向上达到最佳相干性，目的是把有效信号转换到一种在横向上相干性更好的处理域中，然后分别两次通过奇异值分解，提取目标信号的奇异值重构信号，从而实现地震纵、横波波场分离与去噪，这样做的结果是目的性更强，直接针对感兴趣的地震纵、横波成分进行奇异值分解来实现波场分离与去噪的目的，同时避免了有效信号损失严重等缺陷性，具有地震纵、横波分离彻底等特点。
为了实现上述目的，本发明采用的技术和方案是：一种地震纵横波波场分离与去噪的方法，包括如下步骤：
第一步：将含有m道，每道有n个采样点的地震数据读取到二维数组F1中；
第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F1进行正常时差校正，得到L1。正常时差的确定方法为：
$\mathrm{\Delta t}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$
其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t0为零炮检距情况下的旅行时，vP为P波叠加速度；
第三步：对第一次正常时差校正后的地震数据L1进行奇异值分解：
$L_1=U_1{\Sigma }_1{V}_1^T$
其中，上角T表示转置，U1由L1L1T的特征值向量构成，V1由L1TL1的特征值向量构成，∑1由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：
${\Sigma }_1=\left[\begin{array}{cccc}{\sigma }_{\mathrm{1,1}}& & & 0\\ & {\sigma }_{\mathrm{1,2}}& & \\ & & O& \\ 0& & & {\sigma }_{1,m}\end{array}\right]$
其中σ1，1，σ1，2Λσ1，m为L1的奇异值；
第四步：对∑1进行处理，提取目标信号的奇异值∑2重构信号，得到L2；可将信号分为两类进行重构，具体如下：
(1)奇异值分解低通滤波，提取P-P波场：
$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\Sigma }}{\sigma }_{1,j}{u}_j{v}_j^T$
其中：上角T表示转置，LLp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p1为L1的秩，且1≤p1≤m，σ1，j为L1的第j个奇异值，uj为L1L1T的第j个特征向量，vj为L1TL1的第j个特征向量；
(2)奇异值分解高通滤波，分离出其它地震信号：
$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\Sigma }}{\sigma }_{1,j}{u}_j{v}_j^T$
其中：上角T表示转置，LHp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q1为L1的秩，且1≤q1≤m，m为地震总道数，σ1，j为L1的第j个奇异值，uj为L1L1T的第j个特征向量，vj为L1TL1的第j个特征向量；
p1和q1的选择取决于奇异值的相对大小，具体做法是通过奇异值σ1，j的下标j的函数曲线来确定；
根据实际滤波情况，L2为LLp1、LHP1之一；
第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的数据L2进行正常时差反校正后得到F2；对于其它波场，正常时差反校正后得到F3。正常时差的确定方法为：
$\mathrm{\Delta t}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$
其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t0为零炮检距情况下的旅行时，vP为P波叠加速度；
第六步：将正常时差反校正后的数据F2按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪；
第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F3进行正常时差校正，得到L3。正常时差的确定方法为：
$\Delta t_{\mathrm{PS}}=t_{\mathrm{PS}}-t_{0\mathrm{PS}}=\frac{1}{v_P}(\sqrt{x_P^2+z^2}-z)+\frac{1}{v_S}(\sqrt{x_S^2+z^2}-z)$
其中：ΔtPS为转换反射波正常时差，xP为震源到转换点的距离，xS为转换点到接收点的距离，tPS是偏移距为(xP+xS)的旅行时，t0PS为零炮检距情况下的旅行时，vP为P波叠加速度，vS为S波叠加速度，z为反射界面深度；
第八步：对第二次正常时差校正后的地震记录L3进行奇异值分解：
$L_3=U_3{\Sigma }_3{V}_3^T$
其中，其中U3由L3L3T的特征值向量构成，V3由L3TL3的特征值向量构成，∑3由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：
${\Sigma }_3=\left[\begin{array}{cccc}{\sigma }_{\mathrm{3,1}}& & & 0\\ & {\sigma }_{\mathrm{3,2}}& & \\ & & O& \\ 0& & & {\sigma }_{3,m}\end{array}\right]$
其中σ31，σ3，2Λσ3，m为L3的奇异值；
第九步：对∑3进行处理，提取目标信号的奇异值∑4重构信号，得到L4；可将信号分为两类进行重构，具体如下：
(1)奇异值分解低通滤波，提取P-S波场：
$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\Sigma }}{\sigma }_{3,j}{u}_j{v}_j^T$
其中：上角T表示转置，LLp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p2为L3的秩，且1≤p2≤m，σ3，j为L3的第j个奇异值，uj为L3L3T的第j个特征向量，vj为L3TL3的第j个特征向量；
(2)奇异值分解高通滤波，分离出其它噪声。
$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\Sigma }}{\sigma }_{3,j}{u}_j{v}_j^T$
其中：上角T表示转置，LHp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q2为L3的秩，且1≤q2≤m，m也地震总道数，σ3，j为L3的第j个奇异值，uj为L3L3T的第j个特征向量，vj为L3TL3的第j个特征向量；
p2和q2的选择取决于奇异值的相对大小，具体做法是通过奇异值σ3，j的下标j的函数曲线来确定；
根据实际滤波情况，L4为LLp2、LHP2之一；
第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的数据L4进行正常时差反校正后得到F4；对于其它噪声，正常时差反校正后得到F5。正常时差的确定方法为：
${\mathrm{\Delta t}}_{\mathrm{PS}}=t_{\mathrm{PS}}-t_{0\mathrm{PS}}=\frac{1}{v_P}(\sqrt{x_P^2+z^2}-z)+\frac{1}{v_S}(\sqrt{x_S^2+z^2}-z)$
其中：ΔtPS为转换反射波正常时差，xP为震源到转换点的距离，xS为转换点到接收点的距离，tPS是偏移距为(xP+xS)的旅行时，t0PS为零炮检距情况下的旅行时，vP为P波叠加速度，vS为S波叠加速度，z为反射界面深度；
第十一步：将正常时差反校正后的数据F4按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪；
第十二步：将正常时差反校正后的数据F5按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。
由于本发明采用两次正常时差校正与奇异值分解联合处理的手段对地震纵、横波场分离与去噪进行处理，这样做的结果是目的性更强，直接针对感兴趣的地震信号成分(地震纵波、横波或噪声)进行奇异值分解波场分离与去噪，同时避免了以往技术应用空间狭窄，有效信号损失严重等缺陷性，具有地震纵、横波场分离彻底和噪声剔除干净等优点。
附图说明
图1为本发明地震纵横波场分离与去噪处理结果图，其中(a)原始地震记录图，(b)依据纵波叠加速度和P-P波时距规律进行正常时差校正后的地震记录图，(c)第一次正常时差校正后的奇异值和提取重构P-P波的奇异值图，(d)提取的P-P波图，(e)依据纵波叠加速度和P P波时距规律进行正常时差反校正后的P-P波图，(f)剔除P-P波的地震波场图，(g)依据纵、横波的叠加速度和P-S波时距规律进行正常时差校正后的地震记录图，(h)第二次正常时差校正后的奇异值和提取重构P-S波的奇异值图，(i)提取的P-S波图，(j)依据纵、横波叠加速度和P-S波时距规律进行正常时差反校正后的P-S波图，(k)剔除的干扰噪声图，(l)依据纵、横波叠加速度和P-S波时距规律进行正常时差反校正后的噪声图；
图2为本发明地震纵波t0时间和对应的正常时差校正叠加速度参数图；
图3为本发明地震转换横波t0时间和对应的正常时差校正叠加速度参数图。
具体实施方法
下面结合附图对本发明进一步详细说明。
参见附图1、图2和图3，一种地震纵横波波场分离与去噪方法，包括如下步骤：
第一步：将含有m道，每道有n个采样点的地震数据读取到二维数组F1中；
第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F1进行正常时差校正，得到L1。正常时差的确定方法为：
$\mathrm{\Delta t}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$
其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t0为零炮检距情况下的旅行时，vP为P波叠加速度；
第三步：对第一次正常时差校正后的地震数据L1进行奇异值分解：
$L_1=U_1{\Sigma }_1{V}_1^T$
其中，上角T表示转置，U1由L1L1T的特征值向量构成，V1由L1TL1的特征值向量构成，∑1由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：
${\Sigma }_1=\left[\begin{array}{cccc}{\sigma }_{\mathrm{1,1}}& & & 0\\ & {\sigma }_{\mathrm{1,2}}& & \\ & & O& \\ 0& & & {\sigma }_{1,m}\end{array}\right]$
其中σ1，1，σ1，2Λσ1，m为L1的奇异值；
第四步：对∑1进行处理，提取目标信号的奇异值∑2重构信号，得到L2；可将信号分为两类进行重构，具体如下：
(1)奇异值分解低通滤波，提取P-P波场：
$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\Sigma }}{\sigma }_{1,j}{u}_j{v}_j^T$
其中：上角T表示转置，LLp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p1为L1的秩，且1≤p1≤m，σ1，j为L1的第j个奇异值，uj为L1L1T的第j个特征向量，vj为L1TL1的第j个特征向量；
(2)奇异值分解高通滤波，分离出其它地震信号：
$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\Sigma }}{\sigma }_{1,j}{u}_j{v}_j^T$
其中：上角T表示转置，LHp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q1为L1的秩，且1≤q1≤m，m为地震总道数，σ1，j为L1的第j个奇异值，uj为L1L1T的第j个特征向量，vj为L1TL1的第j个特征向量；
p1和q1的选择取决于奇异值的相对大小，具体做法是通过奇异值σ1j的下标j的函数曲线来确定；
根据实际滤波情况，L2为LLp1、LHP1之一；
第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的数据L2进行正常时差反校正后得到F2；对于其它波场，正常时差反校正后得到F3。正常时差的确定方法为：
$\mathrm{\Delta t}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$
其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t0为零炮检距情况下的旅行时，vP为P波叠加速度；
第六步：将正常时差反校正后的数据F2按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪；
第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F3进行正常时差校正，得到L3。正常时差的确定方法为：
${\mathrm{\Delta t}}_{\mathrm{PS}}=t_{\mathrm{PS}}-t_{0\mathrm{PS}}=\frac{1}{v_P}(\sqrt{x_P^2+z^2}-z)+\frac{1}{v_S}(\sqrt{x_S^2+z^2}-z)$
其中：ΔtPS为转换反射波正常时差，xP为震源到转换点的距离，xS为转换点到接收点的距离，tPS是偏移距为(xP+xS)的旅行时，t0PS为零炮检距情况下的旅行时，vP为P波叠加速度，vS为S波叠加速度，z为反射界面深度；
第八步：对第二次正常时差校正后的地震记录L3进行奇异值分解：
$L_3=U_3{\Sigma }_3{V}_3^T$
其中，其中U3由L3L3T的特征值向量构成，V3由L3TL3的特征值向量构成，∑3由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：
${\Sigma }_3=\left[\begin{array}{cccc}{\sigma }_{\mathrm{3,1}}& & & 0\\ & {\sigma }_{\mathrm{3,2}}& & \\ & & O& \\ 0& & & {\sigma }_{3,m}\end{array}\right]$
其中σ31，σ3，2Λσ3，m为L3的奇异值；
第九步：对∑3进行处理，提取目标信号的奇异值∑4重构信号，得到L4；可将信号分为两类进行重构，具体如下：
(1)奇异值分解低通滤波，提取P-S波场：
$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\Sigma }}{\sigma }_{3,j}{u}_j{v}_j^T$
其中：上角T表示转置，LLp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p2为L3的秩，且1≤p2≤m，σ3，j为L3的第j个奇异值，uj为L3L3T的第j个特征向量，vj为L3TL3的第j个特征向量；
(2)奇异值分解高通滤波，分离出其它噪声。
$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\Sigma }}{\sigma }_{3,j}{u}_j{v}_j^T$
其中：上角T表示转置，LHp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q2为L3的秩，且1≤q2≤m，m为地震总道数，σ3，j为L3的第j个奇异值，uj为L3L3T的第j个特征向量，vj为L3TL3的第j个特征向量；
p2和q2的选择取决于奇异值的相对大小，具体做法是通过奇异值σ3，j的下标j的函数曲线来确定；
根据实际滤波情况，L4为LLp2、LHP2之一；
第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的数据L4进行正常时差反校正后得到F4；对于其它噪声，正常时差反校正后得到F5。正常时差的确定方法为：
${\mathrm{\Delta t}}_{\mathrm{PS}}=t_{\mathrm{PS}}-t_{0\mathrm{PS}}=\frac{1}{v_P}(\sqrt{x_P^2+z^2}-z)+\frac{1}{v_S}(\sqrt{x_S^2+z^2}-z)$
其中：ΔtPS为转换反射波正常时差，xP为震源到转换点的距离，xS为转换点到接收点的距离，tPS是偏移距为(xP+xS)的旅行时，t0PS为零炮检距情况下的旅行时，vP为P波叠加速度，vS为S波叠加速度，z为反射界面深度；
第十一步：将正常时差反校正后的数据F4按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪；
第十二步：将正常时差反校正后的数据F5按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。
实施例：
以100道地震记录、采样点为2048的物理模型为例说明实施步骤：
第一步：将含有100道，每道有2048个采样点的地震数据读取到二维数组F1中。
第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F1进行正常时差校正，得到L1。正常时差的确定方法为：
$\mathrm{\Delta t}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$
其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t0为零炮检距情况下的旅行时，vP为P波叠加速度；
第三步：对第一次正常时差校正后的地震数据L1进行奇异值分解：
$L_1=U_1{\Sigma }_1{V}_1^T$
其中，上角T表示转置，U1由L1L1T的特征值向量构成，V1由L1TL1的特征值向量构成，∑1由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：
${\Sigma }_1=\left[\begin{array}{cccc}{\sigma }_{\mathrm{1,1}}& & & 0\\ & {\sigma }_{\mathrm{1,2}}& & \\ & & O& \\ 0& & & {\sigma }_{1,m}\end{array}\right]$
其中σ1，1，σ1，2Λσ1，m为L1的奇异值；
第四步：对奇异值系列∑1进行处理，提取目标信号的奇异值系列∑2重构信号，得到L2。可将信号分为两类进行重构，具体如下：
(1)奇异值分解低通滤波，取前10个奇异值重构信号，提取P-P波场：
$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\Sigma }}{\sigma }_{1,j}{u}_j{v}_j^T$
其中：上角T表示转置，LLp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p1为L1的秩，且1≤p1≤m，σ1，j为L1的第j个奇异值，uj为L1L1T的第j个特征向量，vj为L1TL1的第j个特征向量；
(2)奇异值分解高通滤波，取第11~100个奇异值重构信号，分离出其它地震信号：
$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\Sigma }}{\sigma }_{1,j}{u}_j{v}_j^T$
其中：上角T表示转置，LHp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q1为L1的秩，且1≤q1≤m，m为地震总道数，σ1，j为L1的第j个奇异值，uj为L1L1T的第j个特征向量，vj为L1TL1的第j个特征向量；
第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的信号进行正常时差反校正，回到原始的数据域中。对于P-P波，正常时差反校正后得到F2；对于其它波场，正常时差反校正后得到F3。
正常时差的确定方法为：
$\mathrm{\Delta t}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$
其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t0为零炮检距情况下的旅行时，vP为P波叠加速度；
第六步：将正常时差反校正后的数据F2按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪。
第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F3进行正常时差校正，得到L3。正常时差的确定方法为：
${\mathrm{\Delta t}}_{\mathrm{PS}}=t_{\mathrm{PS}}-t_{0\mathrm{PS}}=\frac{1}{v_P}(\sqrt{x_P^2+z^2}-z)+\frac{1}{v_S}(\sqrt{x_S^2+z^2}-z)$
其中：ΔtPS为转换反射波正常时差，xP为震源到转换点的距离，xS为转换点到接收点的距离，tPS是偏移距为(xP+xS)的旅行时，t0PS为零炮检距情况下的旅行时，vP为P波叠加速度，vS为S波叠加速度，z为反射界面深度；
第八步：对第二次正常时差校正后的地震记录L3进行奇异值分解：
$L_3=U_3{\Sigma }_3{V}_3^T$
其中，其中U3由L3L3T的特征值向量构成，V3由L3TL3的特征值向量构成，∑3由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：
${\Sigma }_3=\left[\begin{array}{cccc}{\sigma }_{\mathrm{3,1}}& & & 0\\ & {\sigma }_{\mathrm{3,2}}& & \\ & & O& \\ 0& & & {\sigma }_{3,m}\end{array}\right]$
其中σ31，σ3，2Λσ3，m为L3的奇异值；
第九步：对奇异值系列∑3进行处理，提取目标信号的奇异值系列∑4重构信号，得到L4。可将信号分为两类进行重构，具体如下：
(1)奇异值分解低通滤波，取前6个奇异值重构信号，提取P-S波场：
$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\Sigma }}{\sigma }_{3,j}{u}_j{v}_j^T$
其中：上角T表示转置，LLp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p2为L3的秩，且1≤p2≤m，σ3，j为L3的第j个奇异值，uj为L3L3T的第j个特征向量，vj为L3TL3的第j个特征向量；
(2)奇异值分解高通滤波，取第7~100个奇异值重构信号，分离出其它噪声：
$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\Sigma }}{\sigma }_{3,j}{u}_j{v}_j^T$
其中：上角T表示转置，LHp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q2为L3的秩，且1≤q2≤m，m为地震总道数，σ3，j为L3的第j个奇异值，uj为L3L3T的第j个特征向量，vj为L3TL3的第j个特征向量；
第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的信号进行正常时差反校正，回到原始的数据域中。对于P-S波，正常时差反校正后得到F4；对于其它噪声，正常时差反校正后得到F5。正常时差的确定方法为：
${\mathrm{\Delta t}}_{\mathrm{PS}}=t_{\mathrm{PS}}-t_{0\mathrm{PS}}=\frac{1}{v_P}(\sqrt{x_P^2+z^2}-z)+\frac{1}{v_S}(\sqrt{x_S^2+z^2}-z)$
其中：ΔtPS为转换反射波正常时差，xP为震源到转换点的距离，xS为转换点到接收点的距离，tPS是偏移距为(xP+xS)的旅行时，t0PS为零炮检距情况下的旅行时，vP为P波叠加速度，vS为S波叠加速度，z为反射界面深度；
第十一步：将正常时差反校正后的数据F4按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪。
第十二步：将正常时差反校正后的数据F5按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。