CN101893719B - 一种地震纵横波波场分离与去噪方法 - Google Patents

Abstract

一种地震纵横波波场分离与去噪的方法，步骤为：将采样点的地震数据读取到二维数组中；将原始地震数据正常时差校正；对第一次正常时差校正后的地震数据奇异值分解；提取目标信号的奇异值重构信号；将重构后的数据进行正常时差反校正；将正常时差反校正后的数据按输入时的地震数据格式输出，实现地震P-P波分离与去噪；将分离P-P波后的地震记录正常时差校正；对第二次正常时差校正后的地震记录奇异值分解；提取目标信号的奇异值重构信号；重构后的数据正常时差反校正；将正常时差反校正后的数据按输入时的地震数据格式输出，实现地震P-S波场分离与去噪；将正常时差反校正后的数据按输入时的地震数据格式输出，具有地震纵、横波分离彻底的特点。

Description

一种地震纵横波波场分离与去噪方法

技术领域

本发明属于地震数据处理技术领域，特别是一种基于正常时差校正(Normal Move Out，简称NMO)与奇异值分解(Singular ValueDecomposition，简称SVD)联合实现地震纵、横波波场分离与去噪的方法和技术。

技术背景

地震波场非常复杂，单一类型地震波激发，当遇到弹性界面以后，波场会发生转换，形成转换波；另外，在野外地震数据采集时，外界各种干扰因素的影响和地震波自身的传播性质，均会带入许多规则和不规则的干扰信号，严重影响了地质信息的提取。地震纵、横波波场分离与去噪非常重要，它是提高地震资料信噪比与分辨率的前提和保障，同时是地震属性信息精细提取的基础。SVD滤波是从多道地震记录中优化提取相关信息，是利用奇异值作为正交基在信号空间正交分解的特征增强相干能量，压制噪声，即它是利用地震波的相干性差异来达到地震波场分离与去噪的方法，如果用较大的奇异值来重构数据就实现了对弱信号和随机噪音的消除；为了消除特定的同相轴，也可以舍去较大的奇异值来重构信号。现有方法几乎都只是用奇异值分解单一方法技术来对地震波场分离与去噪，尤其是不能分离地震纵波与横波，应用空间有限，而且有效信号容易受到破坏，从而妨碍了对地质信息的提取。

发明内容

为了克服上述现有技术的不足，本发明的目的在于提出一种地震纵、横波波场分离与去噪的方法，其依赖地震信号之间的时距规律差异，通过两次正常时差校正处理的手段，分别把P-P波和P-S波校平，使之在横向上达到最佳相干性，目的是把有效信号转换到一种在横向上相干性更好的处理域中，然后分别两次通过奇异值分解，提取目标信号的奇异值重构信号，从而实现地震纵、横波波场分离与去噪，这样做的结果是目的性更强，直接针对感兴趣的地震纵、横波成分进行奇异值分解来实现波场分离与去噪的目的，同时避免了有效信号损失严重等缺陷性，具有地震纵、横波分离彻底等特点。

为了实现上述目的，本发明采用的技术和方案是：一种地震纵横波波场分离与去噪的方法，包括如下步骤：

第一步：将含有m道，每道有n个采样点的地震数据读取到二维数组F₁中；

第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F₁进行正常时差校正，得到L₁。正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第三步：对第一次正常时差校正后的地震数据L₁进行奇异值分解：

$L_1=U_1{\mathrm{\Σ}}_1{V}_1^T$

其中，上角T表示转置，U₁由L₁L₁ ^T的特征值向量构成，V₁由L₁ ^TL₁的特征值向量构成，∑₁由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_1=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{1,1}}& & & 0\\ & {\mathrm{\σ}}_{\mathrm{1,2}}& & \\ & & O& \\ 0& & & {\mathrm{\σ}}_{1,m}\end{array}\right]$

其中σ_1，1，σ_1，2Λσ_1，m为L₁的奇异值；

第四步：对∑₁进行处理，提取目标信号的奇异值∑₂重构信号，得到L₂；可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，提取P-P波场：

$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₁为L₁的秩，且1≤p₁≤m，σ_1，j为L₁的第j个奇异值，u_j为L₁L₁ ^T的第j个特征向量，v_j为L₁ ^TL₁的第j个特征向量；

(2)奇异值分解高通滤波，分离出其它地震信号：

$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₁为L₁的秩，且1≤q₁≤m，m为地震总道数，σ_1，j为L₁的第j个奇异值，u_j为L₁L₁ ^T的第j个特征向量，v_j为L₁ ^TL₁的第j个特征向量；

p₁和q₁的选择取决于奇异值的相对大小，具体做法是通过奇异值σ_1，j的下标j的函数曲线来确定；

根据实际滤波情况，L₂为L_Lp1、L_HP1之一；

第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的数据L₂进行正常时差反校正后得到F₂；对于其它波场，正常时差反校正后得到F₃。正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第六步：将正常时差反校正后的数据F₂按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪；

第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F₃进行正常时差校正，得到L₃。正常时差的确定方法为：

$\mathrm{\Δ}{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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为(x_P+x_S)的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第八步：对第二次正常时差校正后的地震记录L₃进行奇异值分解：

$L_3=U_3{\mathrm{\Σ}}_3{V}_3^T$

其中，其中U₃由L₃L₃ ^T的特征值向量构成，V₃由L₃ ^TL₃的特征值向量构成，∑₃由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_3=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{3,1}}& & & 0\\ & {\mathrm{\σ}}_{\mathrm{3,2}}& & \\ & & O& \\ 0& & & {\mathrm{\σ}}_{3,m}\end{array}\right]$

其中σ₃₁，σ_3，2Λσ_3，m为L₃的奇异值；

第九步：对∑₃进行处理，提取目标信号的奇异值∑₄重构信号，得到L₄；可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，提取P-S波场：

$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₂为L₃的秩，且1≤p₂≤m，σ_3，j为L₃的第j个奇异值，u_j为L₃L₃ ^T的第j个特征向量，v_j为L₃ ^TL₃的第j个特征向量；

(2)奇异值分解高通滤波，分离出其它噪声。

$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₂为L₃的秩，且1≤q₂≤m，m也地震总道数，σ_3，j为L₃的第j个奇异值，u_j为L₃L₃ ^T的第j个特征向量，v_j为L₃ ^TL₃的第j个特征向量；

p₂和q₂的选择取决于奇异值的相对大小，具体做法是通过奇异值σ_3，j的下标j的函数曲线来确定；

根据实际滤波情况，L₄为L_Lp2、L_HP2之一；

第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的数据L₄进行正常时差反校正后得到F₄；对于其它噪声，正常时差反校正后得到F₅。正常时差的确定方法为：

${\mathrm{\Δ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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为(x_P+x_S)的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第十一步：将正常时差反校正后的数据F₄按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪；

第十二步：将正常时差反校正后的数据F₅按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。

由于本发明采用两次正常时差校正与奇异值分解联合处理的手段对地震纵、横波场分离与去噪进行处理，这样做的结果是目的性更强，直接针对感兴趣的地震信号成分(地震纵波、横波或噪声)进行奇异值分解波场分离与去噪，同时避免了以往技术应用空间狭窄，有效信号损失严重等缺陷性，具有地震纵、横波场分离彻底和噪声剔除干净等优点。

附图说明

图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为本发明地震纵波t₀时间和对应的正常时差校正叠加速度参数图；

图3为本发明地震转换横波t₀时间和对应的正常时差校正叠加速度参数图。

具体实施方法

下面结合附图对本发明进一步详细说明。

参见附图1、图2和图3，一种地震纵横波波场分离与去噪方法，包括如下步骤：

第一步：将含有m道，每道有n个采样点的地震数据读取到二维数组F₁中；

第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F₁进行正常时差校正，得到L₁。正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第三步：对第一次正常时差校正后的地震数据L₁进行奇异值分解：

$L_1=U_1{\mathrm{\Σ}}_1{V}_1^T$

其中，上角T表示转置，U₁由L₁L₁ ^T的特征值向量构成，V₁由L₁ ^TL₁的特征值向量构成，∑₁由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_1=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{1,1}}& & & 0\\ & {\mathrm{\σ}}_{\mathrm{1,2}}& & \\ & & O& \\ 0& & & {\mathrm{\σ}}_{1,m}\end{array}\right]$

其中σ_1，1，σ_1，2Λσ_1，m为L₁的奇异值；

第四步：对∑₁进行处理，提取目标信号的奇异值∑₂重构信号，得到L₂；可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，提取P-P波场：

$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₁为L₁的秩，且1≤p₁≤m，σ_1，j为L₁的第j个奇异值，u_j为L₁L₁ ^T的第j个特征向量，v_j为L₁ ^TL₁的第j个特征向量；

(2)奇异值分解高通滤波，分离出其它地震信号：

$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₁为L₁的秩，且1≤q₁≤m，m为地震总道数，σ_1，j为L₁的第j个奇异值，u_j为L₁L₁ ^T的第j个特征向量，v_j为L₁ ^TL₁的第j个特征向量；

p₁和q₁的选择取决于奇异值的相对大小，具体做法是通过奇异值σ_1j的下标j的函数曲线来确定；

根据实际滤波情况，L₂为L_Lp1、L_HP1之一；

第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的数据L₂进行正常时差反校正后得到F₂；对于其它波场，正常时差反校正后得到F₃。正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第六步：将正常时差反校正后的数据F₂按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪；

第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F₃进行正常时差校正，得到L₃。正常时差的确定方法为：

${\mathrm{\Δ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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为(x_P+x_S)的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第八步：对第二次正常时差校正后的地震记录L₃进行奇异值分解：

$L_3=U_3{\mathrm{\Σ}}_3{V}_3^T$

其中，其中U₃由L₃L₃ ^T的特征值向量构成，V₃由L₃ ^TL₃的特征值向量构成，∑₃由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_3=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{3,1}}& & & 0\\ & {\mathrm{\σ}}_{\mathrm{3,2}}& & \\ & & O& \\ 0& & & {\mathrm{\σ}}_{3,m}\end{array}\right]$

其中σ₃₁，σ_3，2Λσ_3，m为L₃的奇异值；

第九步：对∑₃进行处理，提取目标信号的奇异值∑₄重构信号，得到L₄；可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，提取P-S波场：

$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₂为L₃的秩，且1≤p₂≤m，σ_3，j为L₃的第j个奇异值，u_j为L₃L₃ ^T的第j个特征向量，v_j为L₃ ^TL₃的第j个特征向量；

(2)奇异值分解高通滤波，分离出其它噪声。

$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₂为L₃的秩，且1≤q₂≤m，m为地震总道数，σ_3，j为L₃的第j个奇异值，u_j为L₃L₃ ^T的第j个特征向量，v_j为L₃ ^TL₃的第j个特征向量；

p₂和q₂的选择取决于奇异值的相对大小，具体做法是通过奇异值σ_3，j的下标j的函数曲线来确定；

根据实际滤波情况，L₄为L_Lp2、L_HP2之一；

第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的数据L₄进行正常时差反校正后得到F₄；对于其它噪声，正常时差反校正后得到F₅。正常时差的确定方法为：

${\mathrm{\Δ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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为(x_P+x_S)的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第十一步：将正常时差反校正后的数据F₄按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪；

第十二步：将正常时差反校正后的数据F₅按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。

实施例：

以100道地震记录、采样点为2048的物理模型为例说明实施步骤：

第一步：将含有100道，每道有2048个采样点的地震数据读取到二维数组F₁中。

第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F₁进行正常时差校正，得到L₁。正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第三步：对第一次正常时差校正后的地震数据L₁进行奇异值分解：

$L_1=U_1{\mathrm{\Σ}}_1{V}_1^T$

其中，上角T表示转置，U₁由L₁L₁ ^T的特征值向量构成，V₁由L₁ ^TL₁的特征值向量构成，∑₁由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_1=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{1,1}}& & & 0\\ & {\mathrm{\σ}}_{\mathrm{1,2}}& & \\ & & O& \\ 0& & & {\mathrm{\σ}}_{1,m}\end{array}\right]$

其中σ_1，1，σ_1，2Λσ_1，m为L₁的奇异值；

第四步：对奇异值系列∑₁进行处理，提取目标信号的奇异值系列∑₂重构信号，得到L₂。可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，取前10个奇异值重构信号，提取P-P波场：

$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₁为L₁的秩，且1≤p₁≤m，σ_1，j为L₁的第j个奇异值，u_j为L₁L₁ ^T的第j个特征向量，v_j为L₁ ^TL₁的第j个特征向量；

(2)奇异值分解高通滤波，取第11~100个奇异值重构信号，分离出其它地震信号：

$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₁为L₁的秩，且1≤q₁≤m，m为地震总道数，σ_1，j为L₁的第j个奇异值，u_j为L₁L₁ ^T的第j个特征向量，v_j为L₁ ^TL₁的第j个特征向量；

第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的信号进行正常时差反校正，回到原始的数据域中。对于P-P波，正常时差反校正后得到F₂；对于其它波场，正常时差反校正后得到F₃。

正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第六步：将正常时差反校正后的数据F₂按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪。

第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F₃进行正常时差校正，得到L₃。正常时差的确定方法为：

${\mathrm{\Δ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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为(x_P+x_S)的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第八步：对第二次正常时差校正后的地震记录L₃进行奇异值分解：

$L_3=U_3{\mathrm{\Σ}}_3{V}_3^T$

其中，其中U₃由L₃L₃ ^T的特征值向量构成，V₃由L₃ ^TL₃的特征值向量构成，∑₃由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_3=\left[\begin{array}{cccc}{\mathrm{\σ}}_{\mathrm{3,1}}& & & 0\\ & {\mathrm{\σ}}_{\mathrm{3,2}}& & \\ & & O& \\ 0& & & {\mathrm{\σ}}_{3,m}\end{array}\right]$

其中σ₃₁，σ_3，2Λσ_3，m为L₃的奇异值；

第九步：对奇异值系列∑₃进行处理，提取目标信号的奇异值系列∑₄重构信号，得到L₄。可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，取前6个奇异值重构信号，提取P-S波场：

$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₂为L₃的秩，且1≤p₂≤m，σ_3，j为L₃的第j个奇异值，u_j为L₃L₃ ^T的第j个特征向量，v_j为L₃ ^TL₃的第j个特征向量；

(2)奇异值分解高通滤波，取第7~100个奇异值重构信号，分离出其它噪声：

$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₂为L₃的秩，且1≤q₂≤m，m为地震总道数，σ_3，j为L₃的第j个奇异值，u_j为L₃L₃ ^T的第j个特征向量，v_j为L₃ ^TL₃的第j个特征向量；

第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的信号进行正常时差反校正，回到原始的数据域中。对于P-S波，正常时差反校正后得到F₄；对于其它噪声，正常时差反校正后得到F₅。正常时差的确定方法为：

${\mathrm{\Δ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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为(x_P+x_S)的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第十一步：将正常时差反校正后的数据F₄按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪。

第十二步：将正常时差反校正后的数据F₅按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。

Claims (1)

1.一种地震纵、横波波场分离与去噪的方法，其特征在于，包括如下步骤：

第一步：将含有m道，每道有n个采样点的地震数据读取到二维数组F₁中；

第二步：依据纵波叠加速度和P-P反射波时距规律，将原始地震数据F₁进行正常时差校正，得到L₁，正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第三步：对第一次正常时差校正后的地震数据L₁进行奇异值分解：

$L_1=U_1{\mathrm{\Σ}}_1{V}_1^T$

其中，上角T表示转置，U₁由

的特征值向量构成，V₁由

的特征值向量构成，∑₁由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_1=\left[\begin{array}{cccccc}{\mathrm{\σ}}_{\mathrm{1,1}}& & & & & 0\\ & {\mathrm{\σ}}_{\mathrm{1,2}}& & & & \\ & & .& & & \\ & & & .& & \\ & & & & .& \\ 0& & & & & {\mathrm{\σ}}_{1,m}\end{array}\right]$

其中σ_1，1，σ_1，2…σ_1，m为L₁的奇异值；

第四步：对∑₁进行处理，提取目标信号的奇异值∑₂重构信号，得到L₂；可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，提取P-P波场：

$L_{\mathrm{Lp}1}=\underset{j=1}{\overset{p_1}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp1为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₁为L₁的秩，且1≤p₁≤m，σ_1，j为L₁的第j个奇异值，u_j为

的第j个特征向量，v_j为

的第j个特征向量；

(2)奇异值分解高通滤波，分离出其它地震信号：

$L_{\mathrm{Hp}1}=\underset{j=q_1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{1,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp1为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₁为L₁的秩，且1≤q₁≤m，m为地震总道数，σ_1，j为L₁的第j个奇异值，u_j为

的第j个特征向量，v_j为的第j个特征向量；

p₁和q₁的选择取决于奇异值的相对大小，具体做法是通过奇异值σ_1，j的下标j的函数曲线来确定；

根据实际滤波情况，L₂为L_Lp1、L_HP1之一；

第五步：依据纵波叠加速度和P-P反射波时距规律，将重构后的数据L₂进行正常时差反校正后得到F₂；对于其它波场，正常时差反校正后得到F₃，正常时差的确定方法为：

$\mathrm{\Δt}=t-t_0=\sqrt{t_0^2+\frac{x^2}{v_P^2}}-t_0$

其中：Δt为P-P反射波正常时差，x为偏移距，t是偏移距为x的旅行时，t₀为零炮检距情况下的旅行时，v_P为P波叠加速度；

第六步：将正常时差反校正后的数据F₂按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现地震P-P波分离与去噪；

第七步：依据纵、横波叠加速度和P-S反射波时距规律，将分离P-P波后的地震记录F₃进行正常时差校正，得到L₃，正常时差的确定方法为：

$\mathrm{\Δ}{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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为x_P+x_S的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第八步：对第二次正常时差校正后的地震记录L₃进行奇异值分解：

$L_3=U_3{\mathrm{\Σ}}_3{V}_3^T$

其中，其中U₃由的特征值向量构成，V₃由

的特征值向量构成，∑₃由奇异值构成，奇异值由大到小排列在矩阵的主对角线上：

${\mathrm{\Σ}}_3=\left[\begin{array}{cccccc}{\mathrm{\σ}}_{\mathrm{3,1}}& & & & & 0\\ & {\mathrm{\σ}}_{\mathrm{3,2}}& & & & \\ & & .& & & \\ & & & .& & \\ & & & & .& \\ 0& & & & & {\mathrm{\σ}}_{3,m}\end{array}\right]$

其中σ₃₁，σ_3，2…σ_3，m为L₃的奇异值；

第九步：对∑₃进行处理，提取目标信号的奇异值∑₄重构信号，得到L₄；可将信号分为两类进行重构，具体如下：

(1)奇异值分解低通滤波，提取P-S波场：

$L_{\mathrm{Lp}2}=\underset{j=1}{\overset{p_2}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Lp2为奇异值分解低通滤波后重构的信号，j为奇异值序号，p₂为L₃的秩，且1≤p₂≤m，σ_3，j为L₃的第j个奇异值，u_j为

的第j个特征向量，v_j为

的第j个特征向量；

(2)奇异值分解高通滤波，分离出其它噪声，

$L_{\mathrm{Hp}2}=\underset{j=q_2}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\σ}}_{3,j}{u}_j{v}_j^T$

其中：上角T表示转置，L_Hp2为奇异值分解高通滤波后重构的信号，j为奇异值序号，q₂为L₃的秩，且1≤q₂≤m，m为地震总道数，σ_3，j为L₃的第j个奇异值，u_j为的第j个特征向量，v_j为

的第j个特征向量；

p₂和q₂的选择取决于奇异值的相对大小，具体做法是通过奇异值σ_3，j的下标j的函数曲线来确定；

根据实际滤波情况，L₄为L_Lp2、L_HP2之一；

第十步：依据纵、横波叠加速度和P-S反射波时距规律，将重构后的数据L₄进行正常时差反校正后得到F₄；对于其它噪声，正常时差反校正后得到F₅，正常时差的确定方法为：

${\mathrm{\Δ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)$

其中：Δt_PS为转换反射波正常时差，x_P为震源到转换点的距离，x_S为转换点到接收点的距离，t_PS是偏移距为x_P+x_S的旅行时，t_0PS为零炮检距情况下的旅行时，v_P为P波叠加速度，v_S为S波叠加速度，z为反射界面深度；

第十一步：将正常时差反校正后的数据F₄按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震P-S波场分离与去噪；

第十二步：将正常时差反校正后的数据F5按输入时的地震数据格式输出，便完成了正常时差校正与奇异值分解联合实现了地震波场去噪处理。

Images