CN106094033A - The orientation seismic beam forming method of singular value decomposition - Google Patents

Abstract

The present invention relates to the orientation seismic beam forming method of a kind of singular value decomposition, use the linear disturbance that wicket singular value decomposition method extracts low frequency, energy is strong, and filter from original seismic data；Then determine main beam direction according to exploration targets, utilize beam-forming method computation delay parameter, it is thus achieved that orientation earthquake record；Finally for orientation earthquake record recycling singular value decomposition method reconstruct singular value, filter the random background noise in orientation earthquake record, it is thus achieved that the orientation earthquake record of high s/n ratio.Through test, the orientation seismic beam forming method of singular value decomposition can suppress strong linear disturbance and background random noise in earthquake record effectively, and effectively enhance the energy value of target echo, achieve the signal to noise ratio improving geological data in terms of noise reduction and enhancing signal two, improve the quality of data of seismic prospecting.Especially deep earthquake effect is better than singular value decomposition method and orientation seismic beam forming method.

Description

The orientation seismic beam forming method of singular value decomposition

Technical field:

The present invention relates to a kind of orientation seismic beam synthetic method, be for deep earthquake data acquisition and humanity activities The situation that Earthquakes record signal to noise ratio is low, in order to suppress strong linear disturbance and background random noise and strengthen useful signal energy Amount, improves S/N ratio of seismic records, it is proposed that the orientation seismic beam forming method of a kind of singular value decomposition.

Background technology:

Beam-forming thought comes from phased array radar field the earliest.Owing to beam-forming method can effectively be strengthened on objective body Useful signal, is quickly introduced into field of seismic exploration.Disclosed " the local correlation weighting seismic beam synthesis side of CN103984019A Method ", CN104570121A discloses " directionally seismic wave distorted signal removing method " and CN104793243A discloses " based on N The directionally seismic wave data processing method of secondary root superposition " take three kinds of distinct methods to solve main ripple in seismic beam formation technology The problem of Shu Fangxiang external signal distortion；CN103984007A disclose " directionally seismic wave delay parameter Optimization Design " and CN104570097A discloses " orientation earthquake record synthetic method based on Discrete Particle Swarm algorithm " and takes two kinds of methods to solve How optimum option delay parameter, makes the problem that in different earthquake record, target echo realizes in-phase stacking.On address The current Patents delivered and paper are all devoted to study how to strengthen target echo energy, to how to suppress orientation Strong linear disturbance and background random noise in earthquake record do not carry out specific aim research.

Singular value decomposition method is introduced into field of seismic exploration, is mainly used in suppressing strong linear disturbance and background is made an uproar at random Sound.CN102338886A disclose " polarized filtering method of face ripple in a kind of effective attenuation three-component coefficient ", CN102193107A discloses " a kind of seismic wave field separates and denoising method " and CN105319591A discloses " based on radial direction road The SVD self adaptation surface wave pressing method of conversion " all achieve utilize in singular value decomposition method compacting earthquake record strong linear Interference；CN102854533A discloses " a kind of denoising method improving seismic data signal to noise ratio based on wave field separation principle ", has Having suppressed the background random noise in earthquake record to effect, CN104865603A discloses " a kind of SVD filtering for dipping bed Method and device " solve the dipping bed problem that rear wave resistance feature is distorted after filtering during denoising.Above-mentioned And be currently correlated with singular value decomposition document for the research of singular value decomposition compacting noise be all in target echo and On the basis of noise can separate considerably from one another, and point out when the signal to noise ratio of earthquake record is the lowest, target echo and noise When can not efficiently separate, singular value decomposition method is the most applicable.

Summary of the invention:

The purpose of the present invention is that for deep earthquake data acquisition and humanity activities Earthquakes data signal to noise ratio Time low, strengthen the excellent of target echo energy in conjunction with singular value decomposition method compacting noise and orientation seismic beam forming method Point, it is proposed that the orientation seismic beam forming method of singular value decomposition.

Main idea is that: by the orientation seismic beam forming method of singular value decomposition, it is achieved to strong linear Interference filters and the enhancing of target echo energy with background random noise.This method is initially with wicket singular value The linear disturbance that decomposition method extracts low frequency, energy is strong, and filter from original seismic data；Then determine according to exploration targets Main beam direction, utilizes beam-forming method computation delay parameter, it is thus achieved that orientation earthquake record；Finally for orientation earthquake record Recycling singular value decomposition method reconstruct singular value, filters the random background noise in orientation earthquake record, it is thus achieved that high s/n ratio Orientation earthquake record.This method can strengthen target echo while removing strong linear disturbance and background random noise Energy, it is achieved in terms of noise reduction and enhancing signal two, improve S/N ratio of seismic records.

The present invention is achieved by the following technical solutions:

The orientation seismic beam forming method of singular value decomposition, comprises the following steps:

A, on existing earthquake record, choose pending earthquake record U₀(t x), chooses with U₀(t, x) centered by company The earthquake record set U={U of continuous N number of adjacent shot point_-(N-1)/2(t,x),L,U_-1(t,x),U₀(t,x),U₁(t,x),L,U_(N-1)/2 (t, x) }, N is the number (N >=3) of earthquake record, and t is the time collection of earthquake record, and x is the road collection of earthquake record, U₀(t,x) Centered by shot point earthquake record.If U_g(t, x) is arbitrary earthquake record in U, then-(N-1)/2≤g≤(N-1)/2, at U_g(t,x) Middle according to denoising requirement, choose on linear disturbance (including sound wave, face ripple, refracted wave, direct wave etc.) direction to be filtered any 2 A and B, the coordinate of note A is (t₁×f_s,x₁), the coordinate of B is (t₂×f_s,x₂), t₁And t₂It is respectively A, B point coordinates corresponding Then, x₁And x₂It is respectively the Taoist monastic name that A, B point coordinates is corresponding, f_sFor sample rate.Slope according to formula (1) estimation linear disturbance

$k=\frac{x_2-x_1}{f_s\×(t_2-t_1)}---\left(1\right)$

B, on the basis of the k of estimation, before and after k, choose m value respectively at equal intervals form slope collection M, gap size l according to Formula (2) calculates, and m calculates according to formula (3):

$l=\frac{x_2-x_1}{f_s\×(t_2-t_1)-1}-\frac{x_2-x_1}{f_s\×(t_2-t_1)}---\left(2\right)$

$m=\frac{t^{\′}\×f_s}{2}---\left(3\right)$

Wherein, t' is the time width that linear disturbance reads on earthquake record.The slope collection M={k-ml chosen, L, k- Il, L, k-l, k, k+l, L, k+il, L, k+ml} ,-m≤i≤m；

C, as i=0, along slope k+il direction, before and after A point, respectively choose n data point, form data volume W, every number The abscissa y at strong point_jCalculate according to formula (4):

$y_j=\frac{j}{k+il}+t_1\×f_s---\left(4\right)$

Wherein ,-n≤j≤n, define data volume W={U_g(y_-n,x₁-n),L,U_g(y_-1,x₁-1),U_g(y₀,x₁),U_g(y₁,x₁ +1),L,U_g(y_n,x₁+n)}；

D, according to formula (5) calculate data volume W adjacent data amplitude difference sum Δ E_i

${\mathrm{\ΔE}}_i=\underset{j=-n}{\overset{n-1}{\Σ}}{U_g}^2(y_{j+1},x_1-(j+1\left)\right)-{U_g}^2(y_j,x_1-j);---\left(5\right)$

E, make i=-m, L respectively ,-2 ,-1, when 1,2, L, m, calculate in slope collection M residue successively tiltedly according to step c～d The amplitude difference sum that on rate direction, data volume is corresponding, composition set S={ Δ E_-m,L,ΔE_-i,L,ΔE_-1,ΔE₀,ΔE₁,L, ΔE_i,L,ΔE_m, calculate minima Δ E in S according to formula (6)_min, by Δ E_minCorresponding slope direction k+il is designated as linearly The direction of interference:

ΔE_min=min{ Δ E_-m,L,ΔE_-i,L,ΔE_-1,ΔE₀,ΔE₁,L,ΔE_i,L,ΔE_m}； (6)

F, centered by A point, with parallel with slope direction k+il chosen and at a distance of two straight line conducts of t sampling point of Δ Window border during data intercept body X upper and lower, orderThe wicket number comprising linear disturbance is intercepted according to formula (7) According to body X:

$X=\left[\begin{array}{ccccc}{U}_g(y_{-n}-\mathrm{\Δ}t,x_1-n)& L& U_g(y_0-\mathrm{\Δ}t,x_1)& L& U_g(y_n-\mathrm{\Δ}t,x_1+n)\\ M& M& M& M& M\\ U_g(y_{-n},x_1-n)& L& U_g(y_0,x_1)& L& U_g(y_n,x_1+n)\\ M& M& M& M& M\\ U_g(y_{-n}+\mathrm{\Δ}t,x_1-n)& L& U_g(y_0+\mathrm{\Δ}t,x_1)& L& U_g(y_n+\mathrm{\Δ}t,x_1+n)\end{array}\right]---\left(7\right)$

G, according to formula (8) calculate inclined linear interference correction in X is become correction amount delta t of horizontal lineups_j, note correction After data volume be X':

Δt_j=y_j-y₀ (8)

$X^{\′}=\left[\begin{array}{ccccccc}{U}_g(y_{-n}+{\mathrm{\Δt}}_{-n}-\mathrm{\Δ}t,x_1-n)& L& U_g(y_{-n}+{\mathrm{\Δt}}_{-1}-\mathrm{\Δ}t,x_1-1)& U_g(y_0-\mathrm{\Δ}t,x_1)& U_g(y_1+{\mathrm{\Δt}}_1-\mathrm{\Δ}t,x_1+1)& L& U_g(y_n+{\mathrm{\Δt}}_n-\mathrm{\Δ}t,x_1+n)\\ M& & & M& M& & M\\ U_g(y_{-n}+{\mathrm{\Δt}}_{-n},x_1-n)& L& U_g(y_{-n}+{\mathrm{\Δt}}_{-n},x_1-1)& U_g(y_0,x_1)& U_g(y_1+{\mathrm{\Δt}}_1,x_1+1)& L& U_g(y_n+{\mathrm{\Δt}}_n,x_1+n)\\ M& & & M& M& & M\\ U_g(y_{-n}+{\mathrm{\Δt}}_{-n}+\mathrm{\Δ}t,x_1-n)& L& U_g(y_{-n}+{\mathrm{\Δt}}_{-1}+\mathrm{\Δ}t,x_1-1)& U_g(y_0+\mathrm{\Δ}t,x_1)& U_g(y_1+{\mathrm{\Δt}}_1+\mathrm{\Δ}t,x_1+1)& L& U_g(y_n+{\mathrm{\Δt}}_n+\mathrm{\Δ}t,x_1+n)\end{array}\right]---\left(9\right)$

H, according to formula (10) X' carried out singular value decomposition:

$X^{\′}=\underset{p=1}{\overset{r}{\Σ}}{\mathrm{\σ}}_p{u}_p{v_p}^T---\left(10\right)$

Wherein, u_pIt is X'X'^ΤCharacteristic value corresponding eigenvector composition matrix, v_pIt is X'^ΤThe characteristic value of X' is corresponding The matrix of eigenvector composition, σ_pIt is X'X'^Τ(or X'^ΤX') the non-negative square root of characteristic value, the i.e. singular value of X'.Singular value According to the order arrangement successively decreased, σ₁＞ σ₂＞ L ＞ σ_p＞ L ＞ σ_r, r is the order of matrix X', and 1≤p≤r, Τ are transposition symbol；

I, reconstruct singular value, by σ₂,σ₃,L,σ_p,L,σ_rAll it is set to 0, retains first singular value σ₁, extensive according to formula (11) Appear again and only comprise the data volume X of linear disturbance " (not comprising background random noise):

X "=σ₁u₁v₁ ^Τ (11)

Wherein, u₁For X'X'^ΤDominant eigenvalue corresponding eigenvector composition matrix, v₁For X'^ΤThe maximum intrinsic of X' The matrix of the eigenvector composition that value is corresponding；

J, by data volume X " according to formula (12) counter correct back to inclined linear interference, obtain data volume Y, from earthquake record U_g (t, x) in Y is deducted:

$Y=\left[\begin{array}{ccccccc}{U_g}^{\′}(y_{-n}-\mathrm{\Δ}t,x_1-n)& L& {U_g}^{\′}(y_{-n}-\mathrm{\Δ}t,x_1-1)& {U_g}^{\′}(y_0-\mathrm{\Δ}t,x_1)& {U_g}^{\′}(y_1-\mathrm{\Δ}t,x_1+1)& L& {U_g}^{\′}(y_n-\mathrm{\Δ}t,x_1+n)\\ M& & & M& M& & M\\ {U_g}^{\′}(y_{-n},x_1-n)& L& {U_g}^{\′}(y_{-n},x_1-1)& {U_g}^{\′}(y_0,x_1)& {U_g}^{\′}(y_1,x_1+1)& L& {U_g}^{\′}(y_n,x_1+n)\\ M& & & M& M& & M\\ {U_g}^{\′}(y_{-n}+\mathrm{\Δ}t,x_1-n)& L& {U_g}^{\′}(y_{-n}+\mathrm{\Δ}t,x_1-1)& {U_g}^{\′}(y_0+\mathrm{\Δ}t,x_1)& {U_g}^{\′}(y_1+\mathrm{\Δ}t,x_1+1)& L& {U_g}^{\′}(y_n+\mathrm{\Δ}t,x_1+n)\end{array}\right]---\left(12\right)$

K, to earthquake record U_g(t, x) in all linear disturbance process one by one by step a～j, by U_g(t, x) middle institute is wired Property interference remove, obtain earthquake record R_g(t,x)；

L, every for earthquake record set big gun earthquake record is all carried out the process of step a～k, then every big gun earthquake record center line in U Property interference be all removed, obtain new earthquake record set:

R={R_(N-1)/2(t,x),L,R_-1(t,x),R₀(t,x),R₁(t,x),L,R_(N-1)/2(t,x)}；

M, choose a point coordinates (t in target echo center position₃×f_s,x₃), if focal point coordinate is (t₀× f_s,x₀), calculate earthquake main beam direction θ according to formula (13)_max, according to formula (14) computation delay parameter τ:

${\mathrm{\θ}}_{\max }=\arctan \left(\frac{(x_3-x_0)\×D}{(t_3-t_0)\×v/2}\right)---\left(13\right)$

$\mathrm{\τ}=\frac{d{\mathrm{sin\θ}}_{max}}{v}---\left(14\right)$

Wherein, D is detector interval, and d is shot interval, and v is reflecting interface overlying strata speed, obtains by surveying district's data；

Earthquake record set after n, time delay is: R'={R_(N-1)/2'(t-(N-1)τ/2,x),L,R_-1'(t-τ,x),R₀'(t, x),R₁'(t+τ,x),L,R_(N-1)/2' (t+ (N-1) τ/2, x) }, by N big gun earthquake record in R' according to formula (15) superposition, synthesis Orientation earthquake record R " (t, x):

$R^{\′\′}(t,x)=\underset{e=-(N-1)/2}{\overset{(N-1)/2}{\Σ}}{R_e}^{\′}(t+e\mathrm{\τ},x)---\left(15\right)$

Wherein ,-(N-1)/2≤e≤(N-1)/2；

O, set Seismic Traces number as H, utilize method of correlation by R " (t, x) in per pass signal do relevant place to first signal Reason, it is thus achieved that per pass signal is relative to the delay inequality set Δ t of first signal "={ Δ t'₁,Δt'₂,L,Δt'_h,L,Δt '_H, 1≤h≤H, according to formula (16) by R " (t, x) in target echo school become horizontal lineups signal, the ground after smoothing Shake record be designated as r (t, x):

R (t, x)=R " (t-(Δ t_h-Δt₁),x_h) (16)

P, according to formula (17) by r (t, x) carries out singular value decomposition:

$r(t,x)=\underset{q=1}{\overset{r^{\′}}{\Σ}}{{\mathrm{\σ}}^{\′}}_q{u^{\′}}_q{v}^{\′}{{}_q}^T---\left(17\right)$

Wherein, r' is r (t, order x), 1≤q≤r', σ '_qFor r (t, x) r (t, x)^Τ(or r (t, x)^ΤR (t, x)) intrinsic The non-negative square root of value, (i.e. r (t, singular value x)), u'_qFor r (t, x) r (t, x)^ΤEigenvector group corresponding to characteristic value The matrix become, v'_qFor r (t, x)^ΤR (t, the matrix of the eigenvector composition that characteristic value x) is corresponding；

Q, reconstruct singular value, by σ '₃,σ'₄,L,σ'_r' be all set to 0, retain the first two singular value σ '₁,σ'₂, according to formula (18) the data volume r'(t only comprising horizontal lineups is recovered, x):

$r^{\′}(t,x)=\underset{q=1}{\overset{2}{\Σ}}{{\mathrm{\σ}}^{\′}}_q{u^{\′}}_q{v}^{\′}{{}_q}^T---\left(18\right)$

R, according to formula (19) by r'(t, x) in horizontal lineups according to step n is gathered Δ t " anti-school returns hyperbolic homophase Axle, is newly oriented earthquake record r " (t, x):

R " (t, x)=r'(t+ (Δ t'_h-Δt'₁),x_h) (19)

Wherein, r " (t x) is and is ultimately oriented earthquake record.

Beneficial effect: through test, the orientation seismic beam forming method of singular value decomposition disclosed by the invention can be effectively Strong linear disturbance and background random noise in compacting earthquake record, and effectively enhance the energy value of target echo, real Show the signal to noise ratio improving geological data in terms of noise reduction and enhancing signal two, improve the quality of data of seismic prospecting.For Deep and strong man's literary composition Earthquakes data signal to noise ratio are low, and strong linear disturbance is grown, and random noise disturbance is serious, and target reflection is believed Number feature that energy is low, the orientation seismic beam forming method of singular value decomposition can improve such Earthquakes data effectively Quality, effect is better than singular value decomposition method and orientation seismic beam forming method.

Accompanying drawing illustrates:

Fig. 1 two-layer layer-cake model

The beam-forming earthquake recording quality of Fig. 2 singular value decomposition improves effect

A, original seismic data；

Earthquake record after the process of b, orientation seismic beam forming method based on singular value decomposition

Detailed description of the invention:

It is described in further detail with embodiment below in conjunction with the accompanying drawings:

In the present embodiment with d=10m, D=10m, v₁=1500m/s, v₂=2000m/s, f_sCarry out deep as a example by=1000 Portion's earthquake data acquisition and humanity activities area are improved S/N ratio of seismic records and are processed, but deep earthquake data acquisition and people The orientation seismic beam forming method of literary composition movable area application singular value decomposition is not limited by the parameter be given in embodiment.

A, on existing earthquake record, choose pending earthquake record U₀(t x), chooses with U₀(t, x) centered by company The earthquake record set U={U of continuous 7 adjacent shot points_-3(t,x),L,U_-1(t,x),U₀(t,x),U₁(t,x),L,U₃(t, x) }, t is The time collection of earthquake record, x is the road collection of earthquake record, U₀(t, x) centered by shot point earthquake record.If U_g(t is x) to appoint in U One earthquake record, then-3≤g≤3, at U_g(t, x) in arbitrarily choose on face ripple 1 direction two point coordinates A and B, note A be (269, 40), B is (288,41).Slope according to formula (1) estimation linear disturbance

$k=\frac{41-40}{288-269}=\frac{1}{19}\≈0.053---\left(1\right)$

B, by estimation k on the basis of, before and after k, choose m value respectively at equal intervals, gap size l according to formula (2) meter Calculating, m calculates according to formula (3):

$l=\frac{1}{18}-\frac{1}{19}\≈0.003---\left(2\right)$

$m=\frac{0.1\×1000}{2}=50---\left(3\right)$

Wherein, t'=0.1s.It is computed, slope collection M={k-50l, L, k-il, L, k-l, k, k+l, L, the k+il chosen, L, k+50l} ,-50≤i≤50；

C, as i=0, along slope 0.053+0.003i direction, before and after A point, respectively choose 50 data points, form data Body W, the abscissa y of each data point_jCalculate according to formula (4):

$y_j=\frac{j}{0.053+0.003i}+269---\left(4\right)$

Wherein ,-50≤j≤50, data volume W={U_g(y_-50,x₁-50),L,U_g(y_-1,x₁-1),U_g(y₀,x₁),U_g(y₁,x₁ +1),L,U_g(y₅₀,x₁+50)}；

D, according to formula (5) calculate data volume W adjacent data amplitude difference sum Δ E_i

${\mathrm{\ΔE}}_i=\underset{j=-50}{\overset{49}{\Σ}}{U_g}^2(y_{j+1},x_1-(j+1\left)\right)-{U_g}^2(y_j,x_1-j)---\left(5\right)$

E, make i=-50, L respectively ,-2 ,-1,1,2, L, when 50, calculate successively in slope collection M according to step c～d and remain The amplitude difference sum that in slope direction, data volume is corresponding, composition set S={ Δ E_-50,L,ΔE_-i,L,ΔE_-1,ΔE₀,ΔE₁, L,ΔE_i,L,ΔE₅₀, calculate minima Δ E in S according to formula (6)_min, by Δ E_minCorresponding slope direction 0.03 is designated as The direction of linear disturbance:

0.45=min{ Δ E_-50,L,ΔE_-i,L,ΔE_-1,ΔE₀,ΔE₁,L,ΔE_i,L,ΔE₅₀} (6)

F, centered by coordinate points A, with parallel with the slope direction 0.03 chosen and at a distance of two straight lines of 100 sampling points As window border during data intercept body X upper and lower, intercept the wicket data volume X comprising linear disturbance according to formula (7):

$X=\left[\begin{array}{ccccc}{U}_g(y_{-50}-100,x_1-50)& L& U_g(y_0-100,x_1)& L& U_g(y_{50}-100,x_1+50)\\ M& M& M& M& M\\ U_g(y_{-50},x_1-50)& L& U_g(y_0,x_1)& L& U_g(y_{50},x_1+50)\\ M& M& M& M& M\\ U_g(y_{-50}+100,x_1-50)& L& U_g(y_0+100,x_1)& L& U_g(y_{50}+100,x_1+50)\end{array}\right]---\left(7\right)$

G, according to formula (8) calculate inclined linear interference correction in X is become correction amount delta t of horizontal lineups_j, note correction After data volume be X':

Δt_j=y_j-y₀ (8)

$X^{\′}=\left[\begin{array}{ccccccc}{U}_g(y_{-50}+{\mathrm{\Δt}}_{-50}-\mathrm{\Δ}t,x_1-50)& L& U_g(y_{-1}+{\mathrm{\Δt}}_{-1}-\mathrm{\Δ}t,x_1-1)& U_g(y_0-\mathrm{\Δ}t,x_1)& U_g(y_1+{\mathrm{\Δt}}_1-\mathrm{\Δ}t,x_1+1)& L& U_g(y_{50}+{\mathrm{\Δt}}_{50}-\mathrm{\Δ}t,x_1+50)\\ M& & & M& M& & M\\ U_g(y_{-50}+{\mathrm{\Δt}}_{-50},x_1-50)& L& U_g(y_{-1}+{\mathrm{\Δt}}_{-1},x_1-1)& U_g(y_0,x_1)& U_g(y_1+{\mathrm{\Δt}}_1,x_1+1)& L& U_g(y_{50}+{\mathrm{\Δt}}_{50},x_1+50)\\ M& & & M& M& & M\\ U_g(y_{-50}+{\mathrm{\Δt}}_{-50}+\mathrm{\Δ}t,x_1-50)& L& U_g(y_{-1}+{\mathrm{\Δt}}_{-1}+\mathrm{\Δ}t,x_1-1)& U_g(y_0+\mathrm{\Δ}t,x_1)& U_g(y_1+{\mathrm{\Δt}}_1+\mathrm{\Δ}t,x_1+1)& L& U_g(y_{50}+{\mathrm{\Δt}}_{50}+\mathrm{\Δ}t,x_1+50)\end{array}\right]---\left(9\right)$

H, according to formula (10) X' carried out singular value decomposition:

$X^{\′}=\underset{p=1}{\overset{101}{\Σ}}{\mathrm{\σ}}_p{u}_p{v_p}^T---\left(10\right)$

I, reconstruct singular value, by σ₂,σ₃,L,σ_p,L,σ_rAll it is set to 0, retains first singular value σ₁, extensive according to formula (11) Appear again and only comprise the data volume X of linear disturbance " (not comprising background random noise):

X "=σ₁u₁v₁ ^Τ (11)

Wherein, u₁For X'X'^ΤDominant eigenvalue corresponding eigenvector composition matrix, v₁For X'^ΤX' is

The matrix of the eigenvector composition that big characteristic value is corresponding；

J, by data volume X " according to formula (12) counter correct back to inclined linear interference, obtain data volume Y, from earthquake record U_g (t, x) in Y is deducted:

$Y=\left[\begin{array}{ccccccc}{U_g}^{\′}(y_{-50}-\mathrm{\Δ}t,x_1-50)& L& {U_g}^{\′}(y_{-1}-\mathrm{\Δ}t,x_1-1)& {U_g}^{\′}(y_0-\mathrm{\Δ}t,x_1)& {U_g}^{\′}(y_1-\mathrm{\Δ}t,x_1+1)& L& {U_g}^{\′}(y_{50}-\mathrm{\Δ}t,x_1+50)\\ M& & & M& M& & M\\ {U_g}^{\′}(y_{-50},x_1-50)& L& {U_g}^{\′}(y_{-1},x_1-1)& {U_g}^{\′}(y_0,x_1)& {U_g}^{\′}(y_1,x_1+1)& L& {U_g}^{\′}(y_n,x_1+50)\\ M& & & M& M& & M\\ {U_g}^{\′}(y_{-50}+\mathrm{\Δ}t,x_1-50)& L& {U_g}^{\′}(y_{-1}+\mathrm{\Δ}t,x_1-1)& {U_g}^{\′}(y_0+\mathrm{\Δ}t,x_1)& {U_g}^{\′}(y_1+\mathrm{\Δ}t,x_1+1)& L& {U_g}^{\′}(y_{50}+\mathrm{\Δ}t,x_1+50)\end{array}\right]---\left(12\right)$

K, to earthquake record U_g(t, x) in face ripple 2 and sound wave process one by one by step a～j, by U_g(t, x) middle institute is linear Interference is removed, and obtains earthquake record R_g(t,x)；

L, big gun earthquake record every in earthquake record set U is all carried out the process of step a～k, then in U in every big gun earthquake record Linear disturbance is all removed, and obtains new earthquake record set R={R_-3(t,x),L,R_-1(t,x),R₀(t,x),R₁(t,x),L,R₃ (t,x)}；

M, choose a point coordinates (1020,60) in target echo center position, if focal point coordinate is (0,0), Earthquake main beam direction θ is calculated according to formula (13)_max, according to v₁=1500m/s, according to formula (14) computation delay parameter τ:

Earthquake record set after n, time delay is R'={R_-3'(t-3τ,x),L,R_-1'(t-τ,x),R₀'(t,x),R₁'(t+τ, x),L,R₃' (t+3 τ, x) }, by 7 big gun earthquake records in R' according to formula (15) superposition, synthesis orientation earthquake record R " (t, x):

$R^{\′\′}(t,x)=\underset{e=-3}{\overset{3}{\Σ}}{R_e}^{\′}(t+e\mathrm{\τ},x)---\left(15\right)$

Wherein ,-3≤e≤3；

O, set Seismic Traces number as 101, utilize method of correlation by R " (t, x) in per pass signal do relevant to first signal Process, it is thus achieved that per pass signal is relative to the delay inequality set Δ t of first signal "={ Δ t'₁,Δt'₂,L,Δt'_h,L,Δ t'₁₀₁, 1≤h≤101,

According to formula (16) by R " (t, x) in target echo school become horizontal lineups signal, after smoothing earthquake note Record be designated as r (t, x):

R (t, x)=R " (t-(Δ t_h-Δt₁),x_h) (16)

P, according to formula (17) by r (t, x) carries out singular value decomposition:

$r(t,x)=\underset{q=1}{\overset{101}{\Σ}}{{\mathrm{\σ}}^{\′}}_q{u^{\′}}_q{v}^{\′}{{}_q}^T---\left(17\right)$

Wherein, r'=101 is r (t, order x), 1≤q≤101, σ '_qFor r (t, x) r (t, x)^Τ(or r (t, x)^Τr(t, X)) the non-negative square root of characteristic value, (i.e. r (t, singular value x)), u'_qFor r (t, x) r (t, x)^ΤIntrinsic corresponding to characteristic value The matrix of vector composition, v'_qFor r (t, x)^ΤR (t, the matrix of the eigenvector composition that characteristic value x) is corresponding；

Q, reconstruct singular value, by σ '₃,σ'₄,L,σ'_r'All be set to 0, retain the first two singular value σ '₁,σ'₂, according to formula (18) the data volume r'(t only comprising target echo is recovered, x):

$r^{\′}(t,x)=\underset{q=1}{\overset{2}{\Σ}}{{\mathrm{\σ}}^{\′}}_q{u^{\′}}_q{v}^{\′}{{}_q}^T---\left(18\right)$

R, according to formula (19) by r'(t, x) in target echo according to the set Δ t in step o " anti-school returns hyperbolic Lineups, are newly oriented earthquake record r " (t, x):

R " (t, x)=r'(t+ (Δ t'_h-Δt'₁),x_h) (19)

Wherein, r " (t x) is and is ultimately oriented earthquake record.

Claims (1)

1. the orientation seismic beam forming method of a singular value decomposition, it is characterised in that comprise the following steps:

A, on existing earthquake record, choose pending earthquake record U₀(t x), chooses with U₀(t, x) centered by N continuous The earthquake record set U={U of adjacent shot point_-(N-1)/2(t,x),L,U_-1(t,x),U₀(t,x),U₁(t,x),L,U_(N-1)/2(t, x) }, N is the number (N >=3) of earthquake record, and t is the time collection of earthquake record, and x is the road collection of earthquake record, U₀(t, x) centered by big gun Point earthquake record.If U_g(t, x) is arbitrary earthquake record in U, then-(N-1)/2≤g≤(N-1)/2, at U_g(t, x) middle according to going Make an uproar requirement, choose any two points A and B on linear disturbance (including sound wave, face ripple, refracted wave, direct wave etc.) direction to be filtered, The coordinate of note A is (t₁×f_s,x₁), the coordinate of B is (t₂×f_s,x₂), t₁And t₂It is respectively A, B point coordinates corresponding then, x₁With x₂It is respectively the Taoist monastic name that A, B point coordinates is corresponding, f_sFor sample rate.Slope according to formula (1) estimation linear disturbance

$k=\frac{x_2-x_1}{f_s\×(t_2-t_1)}---\left(1\right)$

B, on the basis of the k of estimation, before and after k, choose m value respectively at equal intervals form slope collection M, gap size l is according to formula (2) calculating, m calculates according to formula (3):

$l=\frac{x_2-x_1}{f_s\×(t_2-t_1)-1}-\frac{x_2-x_1}{f_s\×(t_2-t_1)}---\left(2\right)$

$m=\frac{t^{\′}\×f_s}{2}---\left(3\right)$

Wherein, t' is the time width that linear disturbance reads on earthquake record.Slope collection M={k-ml, L, the k-il chosen, L, k-l, k, k+l, L, k+il, L, k+ml} ,-m≤i≤m；

C, as i=0, along slope k+il direction, before and after A point, respectively choose n data point, form data volume W, each data point Abscissa y_jCalculate according to formula (4):

$y_j=\frac{j}{k+il}+t_1\×f_s---\left(4\right)$

Wherein ,-n≤j≤n, define data volume W={U_g(y_-n,x₁-n),L,U_g(y_-1,x₁-1),U_g(y₀,x₁),U_g(y₁,x₁+1), L,U_g(y_n,x₁+n)}；

D, according to formula (5) calculate data volume W adjacent data amplitude difference sum Δ E_i

${\mathrm{\ΔE}}_i=\underset{j=-n}{\overset{n-1}{\mathrm{\Σ}}}{U_g}^2(y_{j+1},x_1-(j+1\left)\right)-{U_g}^2(y_j,x_1-j);---\left(5\right)$

E, make i=-m, L respectively ,-2 ,-1, when 1,2, L, m, calculate residue slope side in slope collection M successively according to step c～d The amplitude difference sum that upwards data volume is corresponding, composition set S={ Δ E_-m,L,ΔE_-i,L,ΔE_-1,ΔE₀,ΔE₁,L,ΔE_i, L,ΔE_m, calculate minima Δ E in S according to formula (6)_min, by Δ E_minCorresponding slope direction k+il is designated as linear disturbance Direction:

ΔE_min=min{ Δ E_-m,L,ΔE_-i,L,ΔE_-1,ΔE₀,ΔE₁,L,ΔE_i,L,ΔE_m}； (6)

F, centered by A point, with parallel with slope direction k+il chosen and at a distance of two straight lines of t sampling point of Δ as intercepting Window border during data volume X upper and lower, orderThe wicket data volume comprising linear disturbance is intercepted according to formula (7) X:

$X=\left[\begin{array}{ccccc}{U}_g(y_{-n}-\mathrm{\Δ}t,x_1-n)& L& U_g(y_0-\mathrm{\Δ}t,x_1)& L& U_g(y_n-\mathrm{\Δ}t,x_1+n)\\ M& M& M& M& M\\ U_g(y_{-n},x_1-n)& L& U_g(y_0,x_1)& L& U_g(y_n,x_1+n)\\ M& M& M& M& M\\ U_g(y_{-n}+\mathrm{\Δ}t,x_1-n)& L& U_g(y_0+\mathrm{\Δ}t,x_1)& L& U_g(y_n+\mathrm{\Δ}t,x_1+n)\end{array}\right]---\left(7\right)$

G, according to formula (8) calculate inclined linear interference correction in X is become correction amount delta t of horizontal lineups_j, after note correction Data volume is X':

Δt_j=y_j-y₀ (8)

$X^{\′}=\left[\begin{array}{ccccccc}{U}_g\left(y_{-n}+{\mathrm{\Δt}}_{-n}-\mathrm{\Δ}t,x_1-n\right)& L& U_g\left(y_{-1}+{\mathrm{\Δt}}_{-1}-\mathrm{\Δ}t,x_1-1\right)& U_g\left(y_0-\mathrm{\Δ}t,x_1\right)& U_g\left(y_1+{\mathrm{\Δt}}_1-\mathrm{\Δ}t,x_1+1\right)& L& U_g\left(y_n+{\mathrm{\Δt}}_n-\mathrm{\Δ}t,x_1+n\right)\\ M& & & M& M& & M\\ U_g\left(y_{-n}+{\mathrm{\Δt}}_{-n},x_1-n\right)& L& U_g\left(y_{-1}+{\mathrm{\Δt}}_{-1},x_1-1\right)& U_g\left(y_0,x_1\right)& U_g\left(y_1+{\mathrm{\Δt}}_1,x_1+1\right)& L& U_g\left(y_n+{\mathrm{\Δt}}_n,x_1+n\right)\\ M& & & M& M& & M\\ U_g\left(y_{-n}+{\mathrm{\Δt}}_{-n}+\mathrm{\Δ}t,x_1-n\right)& L& U_g\left(y_{-1}+{\mathrm{\Δt}}_{-1}+\mathrm{\Δ}t,x_1-1\right)& U_g\left(y_0+\mathrm{\Δ}t,x_1\right)& U_g\left(y_1+{\mathrm{\Δt}}_1+\mathrm{\Δ}t,x_1+1\right)& L& U_g\left(y_n+{\mathrm{\Δt}}_n+\mathrm{\Δ}t,x_1+n\right)\end{array}\right]---\left(9\right)$

H, according to formula (10) X' carried out singular value decomposition:

$X^{\′}=\underset{p=1}{\overset{r}{\Σ}}{\mathrm{\σ}}_p{u}_p{v_p}^T---\left(10\right)$

Wherein, u_pIt is X'X^'TCharacteristic value corresponding eigenvector composition matrix, v_pIt is X^'TThe intrinsic that the characteristic value of X' is corresponding The matrix of vector composition, σ_pIt is X'X^'T(or X^'TX') the non-negative square root of characteristic value, the i.e. singular value of X'.Singular value is according to passing The order arrangement subtracted, σ₁＞ σ₂＞ L ＞ σ_p＞ L ＞ σ_r, r is the order of matrix X', and 1≤p≤r, T are transposition symbol；

I, reconstruct singular value, by σ₂,σ₃,L,σ_p,L,σ_rAll it is set to 0, retains first singular value σ₁, recover according to formula (11) Only comprise the data volume X of linear disturbance " (not comprising background random noise):

X "=σ₁u₁v₁ ^T (11)

Wherein, u₁For X'X^'TDominant eigenvalue corresponding eigenvector composition matrix, v₁For X^'TThe dominant eigenvalue pair of X' The matrix of the eigenvector composition answered；

J, by data volume X " according to formula (12) counter correct back to inclined linear interference, obtain data volume Y, from earthquake record U_g(t,x) Middle Y is deducted:

$Y=\left[\begin{array}{ccccccc}{U_g}^{\′}(y_{-n}-\mathrm{\Δ}t,x_1-n)& L& {U_g}^{\′}(y_{-1}-\mathrm{\Δ}t,x_1-1)& {U_g}^{\′}(y_0-\mathrm{\Δ}t,x_1)& {U_g}^{\′}(y_1-\mathrm{\Δ}t,x_1+1)& L& {U_g}^{\′}(y_n-\mathrm{\Δ}t,x_1+n)\\ M& & & M& M& & M\\ {U_g}^{\′}(y_{-n},x_1-n)& L& {U_g}^{\′}(y_{-1},x_1-1)& {U_g}^{\′}(y_0,x_1)& {U_g}^{\′}(y_1,x_1+1)& L& {U_g}^{\′}(y_n,x_1+n)\\ M& & & M& M& & M\\ {U_g}^{\′}(y_{-n}+\mathrm{\Δ}t,x_1-n)& L& {U_g}^{\′}(y_{-1}+\mathrm{\Δ}t,x_1-1)& {U_g}^{\′}(y_0+\mathrm{\Δ}t,x_1)& {U_g}^{\′}(y_1+\mathrm{\Δ}t,x_1+1)& L& {U_g}^{\′}(y_n+\mathrm{\Δ}t,x_1+n)\end{array}\right]---\left(12\right)$

K, to earthquake record U_g(t, x) in all linear disturbance process one by one by step a～j, by U_g(t, x) middle institute is linear dry Disturb removal, obtain earthquake record R_g(t,x)；

L, every for earthquake record set big gun earthquake record is all carried out the process of step a～k, then U linearly does in every big gun earthquake record Disturb and be all removed, obtain new earthquake record set:

R={R_(N-1)/2(t,x),L,R_-1(t,x),R₀(t,x),R₁(t,x),L,R_(N-1)/2(t,x)}；

M, choose a point coordinates (t in target echo center position₃×f_s,x₃), if focal point coordinate is (t₀×f_s, x₀), calculate earthquake main beam direction θ according to formula (13)_max, according to formula (14) computation delay parameter τ:

${\mathrm{\θ}}_{max}=\arctan \left(\frac{(x_3-x_0)\×D}{(t_3-t_0)\×v/2}\right)---\left(13\right)$

$\mathrm{\τ}=\frac{d{\mathrm{sin\θ}}_{max}}{v}---\left(14\right)$

Wherein, D is detector interval, and d is shot interval, and v is reflecting interface overlying strata speed, obtains by surveying district's data；

Earthquake record set after n, time delay is: R'={R_(N-1)/2'(t-(N-1)τ/2,x),L,R_-1'(t-τ,x),R₀'(t,x), R₁'(t+τ,x),L,R_(N-1)/2' (t+ (N-1) τ/2, x) }, by N big gun earthquake record in R' according to formula (15) superposition, synthesis orientation Earthquake record R " (t, x):

$R^{\′\′}(t,x)=\underset{e=-(N-1)/2}{\overset{(N-1)/2}{\Σ}}{R_e}^{\′}(t+e\mathrm{\τ},x)---\left(15\right)$

Wherein ,-(N-1)/2≤e≤(N-1)/2；

O, set Seismic Traces number as H, utilize method of correlation by R " (t, x) in per pass signal do relevant treatment with first signal, Obtain the delay inequality set relative to first signal of the per pass signal

Δ t "={ Δ t'₁,Δt'₂,L,Δt'_h,L,Δt'_H, 1≤h≤H, according to formula (16) by R " (t, x) in target reflection Signal school become horizontal lineups signal, the earthquake record after smoothing be designated as r (t, x):

R (t, x)=R " (t-(Δ t_h-Δt₁),x_h) (16)

P, according to formula (17) by r (t, x) carries out singular value decomposition:

$r(t,x)=\underset{q=1}{\overset{r^{\′}}{\Σ}}{{\mathrm{\σ}}^{\′}}_q{u^{\′}}_q{v}^{\′}{{}_q}^T---\left(17\right)$

Wherein, r' is r (t, order x), 1≤q≤r', σ '_qFor r (t, x) r (t, x)^T(or r (t, x)^TR (t, x)) characteristic value non- Negative square root, (i.e. r (t, singular value x)), u'_qFor r (t, x) r (t, x)^TCharacteristic value corresponding eigenvector composition square Battle array, v'_qFor r (t, x)^TR (t, the matrix of the eigenvector composition that characteristic value x) is corresponding；

Q, reconstruct singular value, by σ '₃,σ'₄,L,σ'_r'All be set to 0, retain the first two singular value σ '₁,σ'₂, extensive according to formula (18) Appear again and only comprise the data volume r'(t of horizontal lineups, x):

$r^{\′}(t,x)=\underset{q=1}{\overset{2}{\Σ}}{{\mathrm{\σ}}^{\′}}_q{u^{\′}}_q{v}^{\′}{{}_q}^T---\left(18\right)$

R, according to formula (19) by r'(t, x) in horizontal lineups according in step n gather Δ t " hyperbolic lineups are returned in anti-school, To new orientation earthquake record r " (t, x):

R " (t, x)=r'(t+ (Δ t'_h-Δt'₁),x_h) (19)

" (t x) is and is ultimately oriented earthquake record for wherein, r.