CN102338884B - Elliptic window direction band-pass amplitude-preserved filtering data processing method in geophysical prospecting - Google Patents

Abstract

The invention discloses an elliptic window direction band-pass amplitude-preserved filtering data processing method in geophysical prospecting technology. The method comprises the following steps: gridding the gravity and magnetic data by an interpolation method; performing supplementary edging by use of the low-precision gravity and magnetic data, and performing fast Fourier forward transform; rotating the coordinate system so that the filter has directivity; calculating an elliptic window function so that the filter is an elliptic window in the coordinate system of a wave number domain; calculating an initial Hanning window filter and compensating for the initial Hanning window filter; in the wave number domain, multiplying the spectrum of the gravity and magnetic data by a band-pass filtering frequency response function to obtain the filtered wave number domain result; and performing Fourier inverse transform and anti-edging processing to obtain the gravity and magnetic data subjected to band-pass filtering. In the invention, the signal of the filter in a pass band is subjected to amplitude preservation but is totally inhibited in a stop band; the gravity and magnetic abnormal low-frequency background and high-frequency noise can be suppressed, and an intermediate-frequency useful signal is highlighted; and the method is especially suitable for highlighting the intermediate-frequency useful signal in the wave number domain which is an elliptic window.

Description

The logical width filtering data disposal route of protecting of fenestra ovalis direction band in physical prospecting

Technical field

The present invention relates to the geophysical exploration technology field, specifically the logical width filtering data disposal route of protecting of the fenestra ovalis direction band in geophysical survey.

Background technology

Geophysical survey comprises gravity, magnetic force, electrical method, seismic prospecting etc., gravity, magnetic force be develop the earliest, theoretical ripe, widely used two kinds of geophysical methods.Gravity, magnetic force method area coverage are wide, fast economical, are widely used in oil gas, Exploration of solid mineral resources, zone and the every field such as infrastructure research, environment and engineering investigation.Gravity, magnetic exploration comprise data sampling and processing and explain three parts.

Digital filtering is in the digital signal situation, signal to be done filtering to process, and in signal digital was processed, low pass, high pass, bandpass filtering were the most basic the most frequently used filtering, can design these wave filters in frequency field and time domain.

Frequency domain filtering is selected a certain scope on frequency domain, makes the signal in this scope pass through or remove, and can process signal by fast Fourier direct transform (FFT) conversion and FFT inverse transformation.Time-domain filtering is with the method for convolution filtering or recurrence, signal to be carried out digital processing according to needed signal on time domain.

Desirable bandpass filter is passed through all signal lossless in a certain frequency range, and allows the signal beyond this frequency range be completely blocked.Desirable bandpass filter has following two features: one, and bandwidth can be passed through wave filter at passband with interior signal undistortedly; Its two, the passband of wave filter and stopband are well-separated.

In fact, this ideal filter does not exist in practice, and they are a kind of model of problem analysis, and the wave filter in practical application is all that certain of ideal filter is approximate.Desirable frequency response is impossible realize, because its WAVENUMBER RESPONSE is suddenlyd change, therefore, the Gibbs phenomenon that wave filter produces can not be suppressed interference effectively, outstanding useful signal, but also may make the frequency spectrum of useful signal produce distortion.

In practical application, select a continuous function to come the analog bandpass filtering device can effectively eliminate Gibbs phenomenon.Generally, the zone of continuous function some frequency decay between passband and stopband, this part regional signal just has been subject to partly decay.The requirement of desirable Design of Bandpass is exactly that this continuous function should have steeper border, thereby makes attenuation region as much as possible little, so as this wave filter can be to greatest extent near the design of perfect passband.

According to the wave spectrum difference of local anomaly and regional anomaly, Dean (1958), Byerly (1965), Fraser etc. (1966) design digital filter according to the Spectral Characteristic of heavy magnetic anomaly.Darby and Davies (1967), Fuller (1967), and Zurflueh (1967) design two dimensional filter comes the gravity and magnetic data of gridding is extracted the residual gravity magnetic anomaly.

In wavenumber domain, for the radio-frequency component in throttle signal, existing fully smooth function, although can play very large suppression, but simultaneously also can suppress to a certain extent intermediate frequency and low-frequency component, this compacting to intermediate frequency and low-frequency component, be but we do not wish to occur.In order to reduce the compacting of this centering, low-frequency component, make the distortion of useful signal drop to as far as possible minimum, Hou Chongchu (1979,1981) (Hou Chongchu, a kind of frequency filtering method of suppressing interference, physical prospecting and geochemical exploration, 1979, NO.5; Hou Chongchu, FILTERING OF SMOOTH COMPENSATION, petroleum exploration, 1981, NO.2.) a kind of round and smooth low-pass filtering method of compensation that is called proposed, the method usually needs to suppress high frequency interference when explaining gravity and magnetic anomalies and studies.

Xiong Guangchu, Li Xuesheng (1979) (Xiong Guangchu, Li Xuesheng, the concept of directional spectrum and application thereof, the academic collection of thesis of international exchange geology---be that (six) Geophysical Ano Geochemical Exploration geomathematics is write in 26 International Geological Congresses, 1979) the fenestra ovalis anisotropic filter has been proposed, utilized it the long and short axle of semiaxis on two different directions on fenestra ovalis do not wait and the position in the oval center of circle be the center wave number of wave filter strengthen and two different directions of compacting on signal.The heavy magnetic anomaly that it can produce three degree bodies is effectively carried out filtering, particularly for moving towards different local field and regional fields, obvious separating effect is arranged.

Compensate round and smooth filtering and be a low-pass filter (function is limited) that proposes for the heavy magnetic high frequency interference of compacting, and more important in digital processing be the signal of outstanding a certain frequency range, that is need the bandpass filter that different filter windows are arranged.The subject matter of other wave filter except compensating round and smooth filtering has three, and the one, wave filter is not a continuous function, is also to be interrupted suddenly between passband and stopband, thus the Gibbs phenomenon of causing; The 2nd, being not suitable for giving prominence at wavenumber domain is the intermediate frequency useful signal of fenestra ovalis; The 3rd, passband and the attenuation region between stopband of wave filter are larger, passband and stopband do not have well-separated, and wave filter signal in passband is not protected width, and signal does not suppress fully in stopband, not only affect follow-up qualitative analysis, more seriously affect quantitative test.At present, the processing of geophysical information guarantor width is more and more important.

Based on fenestra ovalis anisotropic filter and compensation filter technology, the present invention proposes the logical width wave filter (component inventions of two technology) of protecting of fenestra ovalis direction band.By the compensation filter technology, passband and the attenuation region between stopband of wave filter are very little, are the logical width wave filters of protecting of a continuously differentiable fenestra ovalis direction band.Wave filter signal in passband is protected width, and signal suppresses fully in stopband.When the compensation number of times is zero, be common fenestra ovalis anisotropic filter, be also that common fenestra ovalis anisotropic filter is the logical special case of protecting the width wave filter of fenestra ovalis direction band.

Summary of the invention

The object of the invention be to provide a kind of in passband signal protect width, the logical width filtering data disposal route of protecting of the fenestra ovalis direction band in the physical prospecting that signal suppresses fully in stopband.

The present invention realizes as follows:

1) utilize gravity meter or magnetometer to gather gravity or magnetic data in the open air in the work area;

2) the heavy magnetic data of field acquisition is proofreaied and correct and corrected;

3) adopt interpolation method will weigh the magnetic Information Network and format, the heavy magnetic data of formation rule net f _a(x, y), wherein x, the y plane grid point coordinate of magnetic data of attaching most importance to respectively;

4) utilize the heavy magnetic data of low precision to replenish fringing; With the heavy magnetic data of rule mesh f _a(x, y) changes its line number into and counts and expand as 2 the heavy magnetic data of integer power f _b(x, y);

5) the heavy magnetic data f after utilizing fast Fourier direct transform FFT with fringing _b(x, y) transforms to wavenumber domain, obtains f _bThe frequency spectrum F of (x, y) _b(u, v), F _b(u, v)=FFT[f _b(x, y)], wherein u, v are respectively the wave number of x and y direction;

6) carry out the wavenumber domain rotation of coordinate with following formula, from wavenumber domain uov coordinate system rotation to new wavenumber domain u ' Ov ' coordinate system, thereby make wave filter have directivity:

$\left\{\begin{array}{c}{u}^{\′}=u\cos \mathrm{\α}-v\sin \mathrm{\α}\\ v^{\′}=u\sin \mathrm{\α}+v\cos \mathrm{\α}\end{array}\right.$

Wherein: u ', v ' are respectively the wave number of x ', y ' direction; α is angle or the rotation of coordinate angle of u axle and u ' axle;

7) calculate fenestra ovalis function R with following formula, making wave filter is fenestra ovalis at wavenumber domain u ' Ov ' coordinate:

$R={[\frac{{\left(\right|u^{\′}|-{u^{\′}}_0)}^2}{\mathrm{\Δ}{{u^{\′}}_0}^2}+\frac{{\left(\right|v^{\′}|-{v^{\′}}_0)}^2}{\mathrm{\Δ}{{v^{\′}}_0}^2}]}^{1/2}\·\mathrm{\π}$

Wherein: u ' ₀, v ' ₀Being respectively the center wave number, is the position coordinates of wavenumber domain elliptical center,

Δ u ' ₀, Δ v ' ₀Being respectively the passband half width, is half shaft length of wavenumber domain ellipse, and u ' ₀〉=0, v ' ₀〉=0, Δ u ' ₀＞0, Δ v ' ₀＞0;

8) calculate initial Hanning window wave filter with following formula:

9) with following formula, initial Hanning window wave filter is compensated:

Wherein: p is compensation number of times, positive integer;

Step 9) the compensation number of times p described in, 1≤p≤30.

10) at wavenumber domain, with frequency spectrum F _b(u, v) and bandpass filtering frequency response function

Multiply each other, the filtering that obtains filtered wavenumber domain is F as a result _c(u, v);

11) through fast Fourier inverse transformation IFFT, the heavy magnetic data f after acquisition fenestra ovalis bandpass filtering _c(x, y), f _c(x, y)=IFFT[F _c(u, v)];

12) anti-fringing is processed, from f _cObtain the gravity and magnetic data f after the bandpass filtering of former size of data in (x, y) _d(x, y).

Wave filter of the present invention signal in passband is protected width, signal suppresses fully in stopband, particularly owing to having adopted compensation technique, passband and the attenuation region between stopband of this wave filter are very little, can suppress heavy magnetic anomaly low frequency background and high frequency noise, outstanding intermediate frequency useful signal, being particularly suitable for giving prominence at wavenumber domain is the intermediate frequency useful signal of fenestra ovalis.

Signal in filter transmission band of the present invention has the width of guarantor feature, has directivity, because wave filter is a continuous function, thereby, there is not the Gibbs phenomenon due to the sudden change generation of wave filter.

The present invention not only can be used for the data of heavy magnetic data and process, and also can be used for the processing of the geophysical informations such as electrical method, earthquake.

Description of drawings

Fig. 1 is one dimension fenestra ovalis Spectral Characteristic (different compensation number of times, p=0 is common fenestra ovalis).(along with the increase of compensation number of times p, passband and the attenuation region between stopband of wave filter are more and more less, thereby wave filter has the width of guarantor feature, near theoretical wave filter spectral response feature)

Fig. 2 is two-dimensional elliptic window wave spectrum signature contrast (contrast (right figure) of common fenestra ovalis (left figure) and new method is seen by contrast, and passband and the attenuation region between stopband of new wave filter are very little, thereby wave filter has the width of guarantor feature).

Fig. 3 is the logical Spectral Characteristic of protecting the width wave filter of two-dimensional elliptic window direction band.

Fig. 4 is residual gravity anomaly figure after Bohai gulf basin work area filtering.

Fig. 5 is Bohai gulf basin work area 3D earthquake Tg structural map.

Embodiment

Describe the present invention in detail below in conjunction with accompanying drawing.

The present invention utilizes LCR type gravity meter to gather 1,/50,000 HIGH-PRECISION GRAVITY DATA in certain work area, according to following steps, gravity data is processed:

1) utilize gravity meter or magnetometer to gather gravity data in the open air in the work area;

2) gravity data of field acquisition is proofreaied and correct and corrected;

3) utilize Kriging interpolation method in Surfer software with the gravity data gridding, formation rule net gravity data f _a(x, y), wherein x, y are respectively the plane grid point coordinate of gravity data; At first choice menus order " Grid (grid) | Data (data) ", after selected data file in " Open (opening) " dialog box, click " Open (opening) ", " discrete data interpolation " dialog box appears; Selecting " Kriging (golden method of interpolation in gram) " to click " Ok " from " Gridding Method (gridding method) " group gets final product; The gridded data size is 130 row * 136 row;

4) utilize the gravity data of low precision to replenish fringing; With rule mesh gravity data f _a(x, y) changes its line number into and expands as 2 integer power gravity data f with counting _b(x, y); In order to reach best flange, utilize this to distinguish other gravity data and replenish fringing, with rule mesh gravity data f _a(x, y) changes its line number into and expands as 2 integer power gravity data f with counting _b(x, y); After fringing, the gridded data size is 256 row * 256 row;

5) the gravity data f after utilizing fast Fourier direct transform FFT with fringing _b(x, y) transforms to wavenumber domain, obtains f _bThe frequency spectrum F of (x, y) _b(u, v), F _b(u, v)=FFT[f _b(x, y)], wherein u, v are respectively the wave number of x and y direction;

6) carry out the wavenumber domain rotation of coordinate with following formula, from wavenumber domain uov coordinate system rotation to new wavenumber domain u ' Ov ' coordinate system, thereby make wave filter have directivity:

$\left\{\begin{array}{c}{u}^{\′}=u\cos \mathrm{\α}-v\sin \mathrm{\α}\\ v^{\′}=u\sin \mathrm{\α}+v\cos \mathrm{\α}\end{array}\right.$

Wherein: u ', v ' are respectively the wave number of x ', y ' direction; α is angle or the rotation of coordinate angle of u axle and u ' axle;

7) calculate fenestra ovalis function R with following formula, making wave filter is fenestra ovalis at wavenumber domain u ' Ov ' coordinate:

$R={[\frac{{\left(\right|u^{\′}|-{u^{\′}}_0)}^2}{\mathrm{\Δ}{{u^{\′}}_0}^2}+\frac{{\left(\right|v^{\′}|-{v^{\′}}_0)}^2}{\mathrm{\Δ}{{v^{\′}}_0}^2}]}^{1/2}\·\mathrm{\π}$

Wherein: u ' ₀, v ' ₀Being respectively the center wave number, is the position coordinates of wavenumber domain elliptical center,

Δ u ' ₀, Δ v ' ₀Being respectively the passband half width, is half shaft length of wavenumber domain ellipse, and u ' ₀〉=0, v ' ₀〉=0, Δ u ' ₀＞0, Δ v ' ₀＞0;

8) calculate initial Hanning window wave filter with following formula:

9) with following formula, initial Hanning window wave filter is compensated:

Wherein: p is compensation number of times, positive integer;

Step 9) the compensation number of times p described in, 1≤p≤30.

Fig. 1 is one dimension fenestra ovalis Spectral Characteristic, and in figure, p=0 is common one dimension fenestra ovalis wave spectrum curve.Along with the increase of compensation number of times p, passband and the attenuation region between stopband of wave filter are more and more less, near theoretical wave filter spectral response feature.

Fig. 2 is the contrast of two-dimensional elliptic window wave spectrum signature, and its left figure is common fenestra ovalis wave filter, and right figure is new wave filter, sees by contrast, and passband and the attenuation region between stopband of new method are very little.

Fig. 3 is the logical Spectral Characteristic of protecting the width wave filter of two-dimensional elliptic window direction band, as we can see from the figure, the filter field of wave filter is an ellipse, and wave filter has directivity, is a bandpass filter, attenuation region between passband and stopband is very little, thereby wave filter has the width of guarantor feature.

10) at wavenumber domain, with frequency spectrum F _b(u, v) and bandpass filtering frequency response function

Multiply each other, the filtering that obtains filtered wavenumber domain is F as a result _c(u, v);

11) through fast Fourier inverse transformation IFFT, the gravity data f after acquisition fenestra ovalis bandpass filtering _c(x, y), f _c(x, y)=IFFT[F _c(u, v)];

12) anti-fringing is processed, from f _cObtain the gravimetric data f after the bandpass filtering of former size of data in (x, y) _d(x, y).f _d(x, y) size of data is 130 row * 136 row;

Fig. 4 is the residual gravity anomaly figure that the Bohai gulf basin work area utilizes this method to obtain, and Fig. 5 is Bohai gulf basin work area 3-D seismics Tg structural map.As can be found from the comparison of Figures 4 and 5, that Fig. 4 sentences is ancient 3, the trend of residual gravity anomaly is the north-south to the west of white ancient 1 well, sentence to the east of ancient 3, white ancient 1 well the trend of residual gravity anomaly and be east northeast to, sentencing ancient 3 to there being a north-south tomography near white ancient 1 well.

Fig. 5 seismotectonics trend and gravity anomaly trend, high some position consistency confirm to have verified the achievement of gravity data sentencing ancient 3 to there being a north-south tomography near white ancient 1 well.

Above-mentioned discrete data gridding, the process that fringing is processed, anti-fringing is processed, this professional skill field analyst can skillfully complete.

Technology involved in the present invention also can be used for the data of other geophysical information and processes and the signal digital process field.

Claims (2)

1. the logical width filtering data disposal route of protecting of the fenestra ovalis direction band in a physical prospecting, realize as follows:

1) utilize gravity meter or magnetometer to gather heavy magnetic data in the open air in the work area;

2) the heavy magnetic data of field acquisition is proofreaied and correct and corrected;

3) adopt interpolation method will weigh the magnetic Information Network and format, the heavy magnetic data of formation rule net f _a(x, y), wherein x, the y plane grid point coordinate of magnetic data of attaching most importance to respectively;

4) utilize the heavy magnetic data of low precision to replenish fringing; With the heavy magnetic data of rule mesh f _a(x, y) changes its line number into and counts and expand as 2 the heavy magnetic data of integer power f _b(x, y);

5) the heavy magnetic data f after utilizing fast Fourier direct transform FFT with fringing _b(x, y) transforms to wavenumber domain, obtains f _bThe frequency spectrum F of (x, y) _b(u, v), F _b(u, v)=FFT ［ f _b(x, y) ］, wherein u, v are respectively the wave number of x and y direction;

6) carry out the wavenumber domain rotation of coordinate with following formula, from wavenumber domain uov coordinate system rotation to new wavenumber domain u'Ov' coordinate system, thereby make wave filter have directivity:

Wherein: u ', v ' are respectively the wave number of x ', y ' direction; α is angle or the rotation of coordinate angle of u axle and u ' axle;

7) calculate fenestra ovalis function R with following formula, making wave filter is fenestra ovalis at wavenumber domain u'Ov' coordinate:

Wherein: u ' ₀, v ' ₀Being respectively the center wave number, is the position coordinates of wavenumber domain elliptical center, Δ u ' ₀, Δ v ' ₀Being respectively the passband half width, is half shaft length of wavenumber domain ellipse, and u ' ₀〉=0, v ' ₀〉=0, Δ u ' ₀＞0, Δ v ' ₀＞0;

8) calculate initial Hanning window wave filter with following formula:

9) with following formula, initial Hanning window wave filter is compensated:

Wherein: p is compensation number of times, positive integer;

10) at wavenumber domain, with frequency spectrum F _b(u is v) with the bandpass filtering frequency response function

(u', v') multiplies each other, and the filtering that obtains filtered wavenumber domain is F as a result _c(u, v);

11) through fast Fourier inverse transformation IFFT, the heavy magnetic data f after acquisition fenestra ovalis bandpass filtering _c(x, y), f _c(x, y)=IFFT ［ F _c(u, v) ］;

12) anti-fringing is processed, from f _cObtain the gravity and magnetic data f after the bandpass filtering of former size of data in (x, y) _d(x, y).

2. method according to claim 1, the compensation number of times p described in step 9), 1≤p≤30.

Images