奇异值分解的定向地震波束形成方法
【专利摘要】本发明涉及一种奇异值分解的定向地震波束形成方法，采用小窗口奇异值分解方法提取低频、能量强的线性干扰，并从原始地震资料中滤除；然后根据勘探目标确定主波束方向，利用波束定向方法计算延时参数，获得定向地震记录；最后针对定向地震记录再利用奇异值分解方法重构奇异值，滤除定向地震记录中的随机背景噪声，获得高信噪比的定向地震记录。经试验，奇异值分解的定向地震波束形成方法能有效地压制地震记录中强线性干扰和背景随机噪声，并有效地增强了目标反射信号的能量值，实现了从降噪和增强信号两方面改善地震数据的信噪比，提高了地震勘探的数据质量。尤其是深部地震效果优于奇异值分解方法和定向地震波束形成方法。
【专利说明】
奇异值分解的定向地震波束形成方法
技术领域：
[0001] 本发明涉及一种定向地震波束合成方法，是针对深部地震数据采集以及人文活动 地区地震记录信噪比低的情况，为了压制强线性干扰和背景随机噪声并且增强有效信号能 量，提高地震记录信噪比，提出了一种奇异值分解的定向地震波束形成方法。
【背景技术】：
[0002] 波束定向思想最早源于相控雷达领域。由于波束定向方法能够有效加强目标体上 有效信号，很快被引进地震勘探领域。CN103984019A公开的《局部相关加权地震波束合成方 法》、CN104570121A公开了《定向地震波畸变信号消除方法》和CN104793243A公开了《基于N 次根叠加的定向地震波数据处理方法》采取三种不同方法解决了地震波束形成技术中主波 束方向外信号畸变的问题；CN103984007A公开了《定向地震波延时参数优化设计方法》和 CN104570097A公开了《基于离散微粒群算法的定向地震记录合成方法》采取两种方法解决 了如何优化选取延时参数，使不同地震记录中目标反射信号实现同相叠加的问题。上述及 当前已发表的相关专利和论文均致力于研究如何增强目标反射信号能量，对如何压制定向 地震记录中的强线性干扰和背景随机噪声没有进行针对性研究。
[0003] 奇异值分解方法被引进地震勘探领域，主要用于压制强线性干扰和背景随机噪 声。CN102338886A公开了《一种有效衰减三分量地震记录中面波的极化滤波方法》、 CN102193107A公开了《一种地震波场分离与去噪方法》和CN105319591A公开了《基于径向道 变换的SVD自适应面波压制方法》均实现了利用奇异值分解方法压制地震记录中的强线性 干扰;CN102854533A公开了《一种基于波场分离原理提高地震资料信噪比的去噪方法》，有 效地压制了地震记录中的背景随机噪声，CN104865603A公开了《一种针对倾斜层的SVD滤波 方法及装置》解决了在去噪处理过程中倾斜层经过滤波后波阻特征发生畸变的问题。上述 及当前已有相关奇异值分解文献对于奇异值分解压制噪声的研究均是在目标反射信号和 噪声能够彼此明显分离的基础上，并指出当地震记录的信噪比很低，目标反射信号和噪声 不能有效分离时，奇异值分解方法不再适用。
【发明内容】
：
[0004] 本发明的目的就在于针对深部地震数据采集以及人文活动地区地震资料信噪比 低时，结合奇异值分解方法压制噪声和定向地震波束形成方法增强目标反射信号能量的优 点，提出了奇异值分解的定向地震波束形成方法。
[0005] 本发明的主要思想是:通过奇异值分解的定向地震波束形成方法，实现对强线性 干扰和背景随机噪声的滤除以及目标反射信号能量的增强。本方法首先采用小窗口奇异值 分解方法提取低频、能量强的线性干扰，并从原始地震资料中滤除;然后根据勘探目标确定 主波束方向，利用波束定向方法计算延时参数，获得定向地震记录;最后针对定向地震记录 再利用奇异值分解方法重构奇异值，滤除定向地震记录中的随机背景噪声，获得高信噪比 的定向地震记录。本方法能够在去除强线性干扰和背景随机噪声的同时增强目标反射信号 的能量，实现从降噪和增强信号两方面改善地震记录信噪比。
[0006] 本发明是通过以下技术方案实现的：
[0007] 奇异值分解的定向地震波束形成方法，包括以下步骤：
[0008] a、在已有地震记录上选取待处理的地震记录Uo(t，x)，选取以Uo(t，x)为中心的连 续N个相邻炮点的地震记录集U= {U-(n-i)/2(t，x)，L，U-i(t，x)，Uo(t，x)，Ui(t，x)，L，U(n-ι)/2 (t，x)}，N为地震记录的个数(N彡3)，t为地震记录的时间集，x为地震记录的道集，Uo(t，x) 为中心炮点地震记录。设U g(t，x)为U中任一地震记录，则-(N-l)/2彡g<(N-l)/2,在Ug(t，x) 中根据去噪要求，选取待滤除的线性干扰(包括声波、面波、折射波、直达波等)方向上任意 两点A和B，记A的坐标为(ti X fs，XI)，B的坐标为(t2 X fs，X2)，ti和t2分别为A、B点坐标对应的 到时，xdPx2分别为A、B点坐标对应的道号，匕为采样率。按照公式(1)估算线性干扰的斜率
[0010] b、以估算的k为基准，分别在k前后等间隔选取m个值组成斜率集M，间隔大小1按照 公式(2)计算，m按照公式(3)计算：
[0013]其中，t'为线性干扰在地震记录上读取的时间宽度。选取的斜率集1=仏-1111，1^-il,L,k-l,k,k+l,L,k+il,L,k+ml}
[0014] c、当i = 0时，沿斜率k+i 1方向，在A点前后各选取n个数据点，组成数据体W，每个数 据点的横坐标5^按照公式(4)计算：
[0016] 其中，-η彡 j<n，定义数据体W= {Ug(y-n，xi-n)，L，Ug(y-i，xi-l)，U g(yo，xi)，Ug(yi，xi +1) ,L,Ug(yn,xi+n)}；
[0017] d、按照公式(5)计算数据体W相邻数据的振幅差值之和Δ Ei
[0019 ] e、分别令i = -m，L, -2, -1,1,2, L, m时，按照步骤c~d依次计算出斜率集Μ中剩余斜 率方向上数据体对应的振幅差值之和，组成集合S={AE-m，L，ΔΕ-i，L，ΔΕ-i，ΔΕο, AEhL， Λ Ei，L，Δ Em}，按照公式(6)计算出S中最小值Δ Emin，将Δ Emin对应的斜率方向k+il记为线性 干扰的方向：
[0020] AEmin=min{ ΔΕ-m，L，ΔΕ-i，L，ΔΕ-1，ΔΕο, AEi，L，AEi，L，ΔΕ」；(6)
[0021] f、以A点为中心，用与选取的斜率方向k+il平行且相距At个样点的两条直线作为
截取数据体X的上、下时窗边界， ，按照公式(7)截取包含线性干扰的小窗口数 据体X:
[0026] h、按照公式（10)对X '进行奇异值分解：
[0028]其中，uP是X'X'τ的本征值对应的本征向量组成的矩阵，v^X' TX'的本征值对应的 本征向量组成的矩阵，(或者Χ'ΤΧ'）本征值的非负平方根，即X'的奇异值。奇异值 按照递减的顺序排列，〇i>〇2>L>〇 p>L>〇r，r是矩阵X'的秩，1彡ρ彡r，Τ为转置符号； [0029] i、重构奇异值，将〇2,〇3丄，〇[)丄，(^均置为0，保留第一个奇异值 01，按照公式（11)恢 复出只包含线性干扰的数据体X"（不包含背景随机噪声）：
[0030] X"=〇lUlVlT (11)
[0031] 其中，U1为X'X'M^最大本征值对应的本征向量组成的矩阵，V1SX'TX'的最大本征 值对应的本征向量组成的矩阵；
[0032] j、将数据体X"按照公式（12)反校正回倾斜线性干扰，得到数据体Y，从地震记录Ug (t，x)中将Y减去：
[0023] g、按照公式(8)计算将X中倾斜线性干扰校正成水平同相轴的校正量At」，记校正 后的数据体为X'：
[0034] 1^、对地震记录1^(1：4)中所有线性干扰按步骤3~」逐一处理，将1]8(1：4)中所有线 性干扰去除，得到地震记录Rg(t，x);
[0035] 1、将地震记录集每炮地震记录均进行步骤a~k的处理，则U中每炮地震记录中线 性干扰均被去除，得到新地震记录集：
[0036] R= {R(n-1)/2(t，x)，L，R-i(t，x)，R〇(t，x)，Ri(t，x)，L，R(n-1)/2(t，x)};
[0037] 1]1、在目标反射信号中心点位置选取一点坐标(〖3\；^，13)，设震源点坐标为（1：〇\ 匕，^))，按照公式(13)计算地震主波束方向0_，按照公式(14)计算延时参数1:
[0040] 其中，D为检波器间距，d为炮间距，v为反射界面上覆层速度，通过测区资料获得；
[0041] n、延时后的地震记录集为:= χΚ?να+τ,χ)^，!^-1)/2'（七+0-1)以2 4)}，将1?'中邮包地震记录按照公式（15)叠加，合成 定向地震记录R"（t，x):
[0043] 其中，-(N-l)/2彡e 彡（N-l)/2;
[0044] 〇、设地震记录道数为H，利用相关法将R"（t，x)中每道信号与第一道信号做相关处 理，获得每道信号相对于第一道信号的延时差集合&扩={&1 1，&12，1^，&111，1^八七 4}，1<11<!1，按照公式（16)将矿^1)中目标反射信号校成水平同相轴信号，校平后的地 震记录记为r(t，x):
[0045] r(t，x) =R"（t_( Δ th-Δ ti)，xh) (16)
[0046] p、按照公式（17)将r(t，x)进行奇异值分解：
[0048] 其中，r '为r(t，x)的秩，Kq<r'，〇'q为r(t，x)r(t，x)T(或r(t，x) Tr(t，x))本征 值的非负平方根，（即r(t，x)的奇异值），u'q为r(t，x)r(t，χ)τ的本征值对应的本征向量组 成的矩阵，ν ' q为r (t，X) Tr (t，X)的本征值对应的本征向量组成的矩阵；
[0049] q、重构奇异值，将〇'3,〇'4丄，〇'/均置为0，保留前两个奇异值 〇'1，〇'2，按照公式 (18)恢复出只包含水平同相轴的数据体r'（t，x):
[0051] r、按照公式（19)将r'（t，x)中水平同相轴按照步骤η中集合At"反校回双曲同相 轴，得到新定向地震记录r"（t，x):
[0052] r"（t，x) =r，（t+( Δ t，h-Δ t，i)，xh) (19)
[0053] 其中，r"（t，x)即为最终定向地震记录。
[0054] 有益效果:经试验，本发明公开的奇异值分解的定向地震波束形成方法能有效地 压制地震记录中强线性干扰和背景随机噪声，并有效地增强了目标反射信号的能量值，实 现了从降噪和增强信号两方面改善地震数据的信噪比，提高了地震勘探的数据质量。针对 深部以及强人文地区地震资料信噪比低，强线性干扰发育，随机噪声干扰严重，目标反射信 号能量低的特征，奇异值分解的定向地震波束形成方法能有效地改善该类地区地震资料的 质量，效果优于奇异值分解方法和定向地震波束形成方法。
【附图说明】：
[0055] 图1两层水平层状介质模型
[0056] 图2奇异值分解的波束定向地震记录质量改善效果
[0057] a、原始地震记录；
[0058] b、基于奇异值分解的定向地震波束形成方法处理后地震记录
【具体实施方式】：
[0059]下面结合附图和实施例做进一步的详细说明：
[0060]在本实施例中以d = 10m，D = 10m，vi = 1500m/s，V2 = 2000m/s，fs = 1000为例进行深 部地震数据采集以及人文活动地区改善地震记录信噪比处理，但深部地震数据采集以及人 文活动地区应用奇异值分解的定向地震波束形成方法不受实施例中给出的参数的限制。 [00 61 ] a、在已有地震记录上选取待处理的地震记录Uo(t，x)，选取以Uo(t，x)为中心的连 续7个相邻炮点的地震记录集U= {U-3(t，x)，L，U-i(t，x)，Uo(t，x)，Ui(t，x)，L，U3(t，x)}，t为 地震记录的时间集，x为地震记录的道集，Uo (t，x)为中心炮点地震记录。设Ug (t，x)为U中任 一地震记录，贝>J-3彡g彡3，在Ug(t，x)中面波1方向上任意选取两点坐标A和B，记A为（269， 40)，B为(288，41)。按照公式(1)估算线性干扰的斜率
[0063] b、以估算的k为基准，分别在k前后等间隔选取m个值，间隔大小1按照公式（2)计 算，m按照公式(3)计算：
[0066] 其中，t' = 0· Is。经计算，选取的斜率集M= {k-501，L，k-il，L，k-l，k，k+l，L，k+il， L，k+501}，-50彡i彡50;
[0067] c、当i = 0时，沿斜率0 · 053+0 · 003i方向，在A点前后各选取50个数据点，组成数据 体W，每个数据点的横坐标η按照公式(4)计算：
[0069] 其中，_50彡 j彡50,数据体W= {Ug(y-5〇，χι_50)，L，Ug(y-i，xi-l)，Ug(yo，xi)，U g(yi，xi +1) ,L,Ug(y5〇,xi+50)}；
[0070] d、按照公式(5)计算数据体W相邻数据的振幅差值之和Δ Ei
[0072] e、分别令i = -50，L，-2，-1，1，2，L，50时，按照步骤c~d依次计算出斜率集Μ中剩余 斜率方向上数据体对应的振幅差值之和，组成集合S = { Δ E-5Q，L，Δ E-i，L，Δ E-i，Δ Εο，Δ Ei， L，Δ Ei，L，Δ E5Q}，按照公式(6)计算出S中最小值Δ Emin，将Δ Emin对应的斜率方向0.03记为 线性干扰的方向：
[0073] 〇.45=min{ ΔΕ-5q，L，ΔΕ-i，L，ΔΕ-1，ΔΕ。，AEi，L，AEi，L，ΔΕ5〇} (6)
[0074] f、以坐标点A为中心，用与选取的斜率方向0.03平行且相距100个样点的两条直线作为截取数据体X的上、下时窗边界，按照公式(7)截取包含线性干扰的小窗口数据体X:
[0076] g、按照公式(8)计算将X中倾斜线性干扰校正成水平同相轴的校正量△ 记校正后的数据体为X'：
[0079] h、按照公式（10)对X'进行奇异值分解：
[0081] i、重构奇异值，将〇2,〇3丄，〇[)丄，(^均置为0，保留第一个奇异值 01，按照公式（11)恢 复出只包含线性干扰的数据体X"（不包含背景随机噪声）：
[0082] X，，=〇lUlVlT (11)
[0083] 其中，m为X ' X 'τ的最大本征值对应的本征向量组成的矩阵，V1*X ' TX '的最
[0084] 大本征值对应的本征向量组成的矩阵；
[0085] j、将数据体X"按照公式（12)反校正回倾斜线性干扰，得到数据体Υ，从地震记录Ug (t，x)中将Y减去：
[0087] k、对地震记录Ug(t，x)中面波2和声波按步骤a~j逐一处理，将Ug(t，x)中所有线性 干扰去除，得到地震记录Rg(t，x);
[0088] 1、将地震记录集U中每炮地震记录均进行步骤a~k的处理，则U中每炮地震记录中 线性干扰均被去除，得到新地震记录集R= {R-3(t，x)，L，R-i(t，x)，R〇(t，x)，Ri(t，x)，L，R3 (t,x)}；
[0089] m、在目标反射信号中心点位置选取一点坐标（1020，60)，设震源点坐标为(0，0)， 按照公式(13)计算地震主波束方向0 max，根据Vl=1500m/s，按照公式(14)计算延时参数τ:
[0092] n、延时后的地震记录集为= X)，L，R3'（?+3τ，χ)}，将R'中7炮地震记录按照公式（15)叠加，合成定向地震记录R"（t，x):
[0094] 其中，-3彡e彡3;
[0095] 〇、设地震记录道数为101，利用相关法将R"（t，x)中每道信号与第一道信号做相关 处理，获得每道信号相对于第一道信号的延时差集合八扩={八11，八12，1^，八111，1^八 fioila^h^lOl,
[0096] 按照公式（16)将R"（t，x)中目标反射信号校成水平同相轴信号，校平后的地震记 录记为r(t，x):
[0097] r(t，x) =R"（t_( Δ th-Δ ti)，xh) (16)
[0098] p、按照公式（17)将r(t，x)进行奇异值分解：
[0100]其中，r ' = 101为r(t，x)的秩，Kq<101，〇 'q为r(t，x)r(t，x)T(或r(t，x)Tr(t， x))本征值的非负平方根，（即r(t，x)的奇异值），u'q为r(t，x)r(t，x)T的本征值对应的本征 向量组成的矩阵，v ' q为r (t，X) Tr (t，X)的本征值对应的本征向量组成的矩阵；
[0101] q、重构奇异值，将〇'3,〇'4，L，〇V均置为0,保留前两个奇异值o'U'2,按照公式 (18)恢复出只包含目标反射信号的数据体r'（t，x):
[0103] r、按照公式（19)将r'（t，x)中目标反射信号按照步骤〇中的集合At"反校回双曲 同相轴，得到新定向地震记录r"（t，x):
[0104] r"（t，x) =r '（t+( Δ t'h-Δ t'1)，xh) (19)
[0105] 其中，r"（t，x)即为最终定向地震记录。
【主权项】
1. 一种奇异值分解的定向地震波束形成方法，其特征在于，包括W下步骤： 曰、在已有地震记录上选取待处理的地震记录化(t，x)，选取W化(t，x)为中屯、的连续N个 相邻炮点的地震记录集 U=化-(N-l)/2(t，X)，L，U-l(t，X)，U〇(t，X)，Ul(t，X)，L，U(N-l)/2(t，X)}， N为地震记录的个数(N>3)，t为地震记录的时间集，X为地震记录的道集，Uo(t，x)为中屯、炮 点地震记录。设Ug(t，x)为U中任一地震记录，则-(N-l)/2《g《(N-l)/2,在Ug(t，x)中根据去 噪要求，选取待滤除的线性干扰(包括声波、面波、折射波、直达波等)方向上任意两点A和B， 记A的坐标为(tlXfs，Xl)，B的坐标为(t2Xfs，X2)，tl和t2分别为A、B点坐标对应的到时，Xl和 X2分别为A、B点坐标对应的道号，fs为采样率。按照公式(1)估算线性干扰的斜率b、 W估算的k为基准，分别在k前后等间隔选取m个值组成斜率集M，间隔大小1按照公式 (2)计算，m按照公式(3)计算：其中，t'为线性干扰在地震记录上读取的时间宽度。选取的斜率集M=化-ml, L，k-il, L,k-1 ,k,k+l ,L,k+il ,L,k+ml} c、 当i = 0时，沿斜率k+il方向，在A点前后各选取n个数据点，组成数据体W，每个数据点 的横坐标yd安照公式(4)计算：其中，-n《j《n，定义数据体W=化g(y-n，x广η)，L，Ug(y-i，x广 1)，Ug(yo，xi)，Ug(yi，xi+l)， l^，Ug(yn，xi+n)}; d、 按照公式巧)计算数据体W相邻数据的振幅差值之和Δ Ele、 分别令i = -m，L，-2，-1，1，2，L，m时，按照步骤c~d依次计算出斜率集Μ中剩余斜率方 向上数据体对应的振幅差值之和，组成集合S = { Δ E-m，L，Δ E-i，L，Δ Ε-1，Δ Εο，Δ El，L，Δ Ei， L，Δ Em}，按照公式(6 )计算出S中最小值Δ Emin，将Δ Emin对应的斜率方向k+ i 1记为线性干扰 的方向： AEmin=min{ AE-i,L ΔΕ-1, ΔΕο, ΔΕι,Ι^, AEi,L, AEm}; (6) f、 WA点为中屯、，用与选取的斜率方向k+il平行且相距At个样点的两条直线作为截取 数据体X的上、下时窗边界，令心>^，按照公式(7)截取包含线性干扰的小窗口数据体 X:g、 按照公式(8)计算将X中倾斜线性干扰校正成水平同相轴的校正量Δ tj，记校正后的 数据体为X' ：h、 按照公式(10)对X '进行奇异值分解：其中，Up是Χ'Χ'τ的本征值对应的本征向量组成的矩阵，Vp是Χ'τχ'的本征值对应的本征 向量组成的矩阵，曰Ρ是Χ'Χ'τ(或者χτχ'）本征值的非负平方根，旨帖'的奇异值。奇异值按照递 减的顺序排列，〇i>(52>L>〇p>L>〇r，r是矩阵X'的秩，l《p《r，T为转置符号； i、 重构奇异值，将〇2,〇3心〇。山〇朝置为0，保留第一个奇异值〇1，按照公式（11)恢复出 只包含线性干扰的数据体X"(不包含背景随机噪声）：其中，U1为Χ'Χ'τ的最大本征值对应的本征向量组成的矩阵，VI为χτχ'的最大本征值对 应的本征向量组成的矩阵； j、 将数据体X"按照公式(12)反校正回倾斜线性干扰，得到数据体Υ，从地震记录Ug(t，x) 中将Y减去：k、 对地震记录Ug(t，x)中所有线性干扰按步骤a~j逐一处理，将Ug(t，x)中所有线性干 扰去除，得到地震记录Rg(t，x); l、 将地震记录集每炮地震记录均进行步骤a~k的处理，则U中每炮地震记录中线性干 扰均被去除，得到新地震记录集： R= {R(N-i)/2(t，x)，LR-i(t，x)，R〇(t,x) ,Ri(t,x) ,LR(N-i)/2(t,x)}; m、 在目标反射信号中屯、点位置选取一点坐标（t3Xfs，X3)，设震源点坐标为（toXfs， xo)，按照公式(13)计算地震主波束方向0max，按照公式(14)计算延时参数τ:其中，D为检波器间距，d为炮间距，V为反射界面上覆层速度，通过测区资料获得； η、延时后的地震记录集为：R' = {R(N-i)/2'（t-(N-l)T/2，x)，L，R-i'（t-T，x)，R〇'（t，x), Rl'（t+τ，χ)，L，R(N-l)/2'（t+(N-l)τ/2，x)}，将R'中N炮地震记录按照公式（15)叠加，合成定向 地震记录R"(t，x):其中，-(N-l)/2《e《（N-l)/2; O、 设地震记录道数为H，利用相关法将R"(t，x)中每道信号与第一道信号做相关处理， 获得每道信号相对于第一道信号的延时差集合 At" = { At'i，At'2，L At'h，L，Δt'H}，l《h《H，按照公式（16)将R"(t，x)中目标反射 信号校成水平同相轴信号，校平后的地震记录记为r(t，x): ;r(t ,χ) =R" (t-( Δ th-A ti) ,xh) (16) P、 按照公式(17)将r (t，X)进行奇异值分解：其中，r'为r(t，x)的秩，，〇'q为；r(t，x);r(t，x)T(或；r(t，x)T；r(t，x))本征值的非 负平方根，（即r(t，x)的奇异值），u'q为Ht，x);r(t，x)T的本征值对应的本征向量组成的矩 阵，v'q为r(t，x)Tr(t，x)的本征值对应的本征向量组成的矩阵； q、 重构奇异值，将0'3,0'4，L，0V均置为0,保留前两个奇异值〇'1，〇'2,按照公式（18)恢 复出只包含水平同相轴的数据体r'（t，x):r、 按照公式（19)将r'（t，x)中水平同相轴按照步骤η中集合At"反校回双曲同相轴，得 到新定向地震记录r"(t,x): r" (t ,χ) =! '（t+( Δ t 'h- A t'1),祉）(19) 其中，r"(t，x)即为最终定向地震记录。
【文档编号】G01V1/36GK106094033SQ201610389602
【公开日】2016年11月9日
【申请日】2016年6月5日 公开号201610389602.X, CN 106094033 A, CN 106094033A, CN 201610389602, CN-A-106094033, CN106094033 A, CN106094033A, CN201610389602, CN201610389602.X
【发明人】姜弢, 邹艳艳, 徐学纯, 岳永高, 汪彦龙, 晁云峰, 贾海青
【申请人】吉林大学