CN103064111A - Micro seismic event recognition method based on morphological filtering - Google Patents

Abstract

The invention provides a micro seismic event recognition method based on morphological filtering. The method comprises the steps of conducting spectral analysis on micro seismic data; rectifying a demodulator used for receiving micro seismic signals; selecting structural elements, combining opening operation and closing operation to form a morphological noise filter, and conducting the morphological filtering on the micro seismic data going through spectral analysis; calculating the energy ratio of a long-time window and a short-time window according to the data going through the morphological filtering, and comparing the energy ratio and triggering threshold to identify whether the demodulator has micro seismic events; and conducting microseism pickup according to the micro seismic events. According to the method, the signal to noise ratio of the data can be improved through preprocessing and further filtering treatment; the micro seismic events can be detected quickly and accurately by utilization of the energy ratio method, the fluctuation type of the micro seismic events can be identified through polarization analysis, and first arrival can be picked up accurately by utilization of self-adaptation threshold and the long-time and short-time window energy ratio method.

Description

A kind of microearthquake event recognition method based on shape filtering

Technical field

The invention relates to the geophysical prospecting for oil technology, particularly about a kind of microearthquake event recognition method based on shape filtering.

Background technology

It is significant to oil-field development with the microearthquake method field produces to be carried out real-time dynamic monitoring.Because waterfrac treatment needs to obtain more fast the pressure break result, thereby fracturing effect is estimated, and then recovery scheme is in time adjusted, this just need to carry out real-time processing to waterfrac treatment microearthquake Monitoring Data and explain.Currently used microearthquake Monitoring Data disposal route is mainly take artificial cognition, first arrival inverting as the basis.Because the microearthquake Monitoring Data data recording time is long, data volume is large, adopt artificial cognition and the microearthquake validity event is carried out first arrival pick up and to expend very large manpower and time, process on the one hand a waterfrac treatment Monitoring Data and obtain the fracture distribution image and need the long time can not satisfy the demand of Real-Time Monitoring; This method can't be applicable to the oil and gas development monitoring on the other hand.In addition, the microearthquake signal is more weak, and signal to noise ratio (S/N ratio) is low, the identification difficulty of microearthquake event.

Summary of the invention

The invention provides a kind of microearthquake event recognition method based on shape filtering, with the automatic accurate first break picking of microearthquake.

To achieve these goals, the invention provides a kind of microearthquake event recognition method based on shape filtering, the method comprises: step 1: the microearthquake data are carried out spectrum analysis; Step 2: the wave detector that is used for receiving the microearthquake signal is proofreaied and correct; Step 3: the selecting structure element, will open computing and combine with closure operation and form the morphology noise filter, the described geological data through spectrum analysis is carried out shape filtering; Step 4: according to the energy Ratios that carries out data behind the shape filtering and calculate the length time window, described energy Ratios and activation threshold value are compared to identify described wave detector whether have the microearthquake event; Step 5: carry out microearthquake according to described microearthquake event and pick up.

Further, the wave detector that is used for receiving the microearthquake signal is proofreaied and correct, being comprised: the signal correction of the x component of wave detector on component direct P wave polarization direction, is corrected to the y component on the direction perpendicular to component direct P wave polarization direction.

Further, with the signal correction of the x component of wave detector on component direct P wave polarization direction, the y component is corrected on the direction perpendicular to component direct P wave polarization direction, comprise: obtain the polarization direction of the through compressional wave of first arrival that each wave detector receives and the angle theta of wave detector x component according to the energy of ripple, the microearthquake signal that then 2 horizontal components of each wave detector is received rotates θ corresponding to each wave detector.

Further, described selecting structure element comprises: the width of selecting structure element, height and shape.

Further, for multi-component and multilevel geophone, described step 4 also comprises: open a fenestella according to effective geological data of the microearthquake event of identifying at the wave detector of the rule that occurs between the phase neighboring trace outside described wave detector and detect.

Further, described step 5 comprises: utilize erosion algorithm to remove random noise; Choose sizeable semicircular structure element, after described microearthquake data take absolute value behind the shape filtering, carry out corrosion treatment; Choose the mean value of energy value of the signal after the corrosion treatment as the first arrival activation threshold value; Get forward window when sizeable at activation threshold value point place; Choose sizeable long short time-window, when described in the window take a sampled point as step-length, carry out the energy ratio between the long window and the short window computing; The point of determining the energy ratio maximum is the first arrival point.

The beneficial effect of the embodiment of the invention is, the present invention processes the signal to noise ratio (S/N ratio) that can improve data by pre-service and further filtering, utilize energy ratio function can detect rapidly and accurately the microearthquake event, and identify the fluctuation type of microearthquake event by polarization analysis, utilize automatically accurate first break picking of the adaptive threshold of design and energy ratio between the long window and the short window method.

Description of drawings

In order to be illustrated more clearly in the embodiment of the invention or technical scheme of the prior art, the below will do to introduce simply to the accompanying drawing of required use in embodiment or the description of the Prior Art, apparently, accompanying drawing in the following describes only is some embodiments of the present invention, for those of ordinary skills, under the prerequisite of not paying creative work, can also obtain according to these accompanying drawings other accompanying drawing.In the accompanying drawings:

Fig. 1 is that the embodiment of the invention flies microearthquake event recognition method process flow diagram;

Fig. 2 is the spectrum diagram of the microearthquake signal that receives of wave detector x, y, the z component of the embodiment of the invention;

Fig. 3 is one microearthquake synoptic diagram of the bandpass filtering front and back of the embodiment of the invention;

Fig. 4 is the three-component microearthquake signal record synoptic diagram that the embodiment of the invention extracts;

The variation synoptic diagram in three-component seismometer orientation when Fig. 5 is embodiment of the invention microearthquake observation;

Fig. 6 a is one section horizontal component synoptic diagram of perforation data;

Fig. 6 b is the hodograph synoptic diagram that draws according to Fig. 6 a;

Fig. 6 c is the curve synoptic diagram that draws according to energy criteria;

Fig. 7 a is the horizontal component synoptic diagram behind the rotation θ angle;

Fig. 7 b is the hodograph synoptic diagram of making according to horizontal component after the rotation;

Fig. 8 is the hodograph synoptic diagram that horizontal component location microearthquake the last period data are made;

Fig. 9 is the hodograph synoptic diagram of making after being transformed into horizontal component take component direct P wave polarization direction as the coordinate system of reference;

Figure 10 a is the corrosion signal schematic representation that utilizes the semicircular structure element;

Figure 10 b is the expansion signal schematic representation that utilizes the semicircular structure element;

Figure 11 is gray scale opening and closing synoptic diagram;

Figure 12 a is the opening and closing operation result schematic diagram;

Figure 12 b is make and break operation result synoptic diagram;

Figure 12 c is the average result synoptic diagram of opening and closing operation and make and break computing;

Figure 12 d is the comparing result synoptic diagram of opening and closing operation and make and break computing;

Figure 13 is one embodiment of the invention seismic wavelet model synoptic diagram;

Figure 14 is the signal schematic representation behind one embodiment of the invention different size structural element morphologic filtering;

Figure 15 is another embodiment of the present invention seismic wavelet model synoptic diagram;

Figure 16 is that the different semi-circular structure element morphology of radius is learned filtered signal schematic representation;

Figure 17 is another embodiment of the present invention seismic wavelet model synoptic diagram;

Figure 18 is sinusoidal structural element filter effect synoptic diagram;

Figure 19 is semicircular structure element filter effect synoptic diagram;

Figure 20 is flat-structure element filter effect synoptic diagram;

Figure 21 is triangular structure element filter effect synoptic diagram;

Figure 22 is pretreated one microearthquake schematic diagram data;

Figure 23 is to its closed through opening, closed schematic diagram data of opening after processing;

Figure 24 is the three-component energy Ratios curve synoptic diagram of certain road wave detector;

Figure 25 is the signal schematic representation after certain road original signal and the corrosion treatment;

Figure 26 is the rear signal energy value synoptic diagram of corrosion;

Figure 27 is signal energy figure after the local corrosion of amplifying;

Figure 28 is one section untreated original microearthquake schematic diagram data;

Figure 29 is through the schematic diagram data after the bandpass filtering treatment;

Figure 30 is the geological data synoptic diagram after corresponding diagram 4 conversion;

Figure 31 is one embodiment of the invention seismic wavelet model synoptic diagram;

Figure 32 to Figure 34 is that three kinds of wave filters are to the effect synoptic diagram of random noise compacting;

Figure 35 to Figure 38 is the pressing result figure of three kinds of wave filter paired pulses noises;

Figure 39 is pretreated one section microearthquake schematic diagram data;

Figure 40 is closed through opening, closed energy diagram of opening after processing;

Figure 41 is the energy Ratios curve synoptic diagram of the 12 road wave detectors that calculate of energy ratio function;

Figure 42 is the polarization coefficient curve synoptic diagram of 12 road wave detectors;

Figure 43 is energy ratio figure synoptic diagram;

Figure 44 is the theoretical model result schematic diagram of time difference method location;

Figure 45 is the ray tracing location algorithm theoretical model result schematic diagram based on polarization analysis.

Embodiment

For the purpose, technical scheme and the advantage that make the embodiment of the invention is clearer, below in conjunction with accompanying drawing the embodiment of the invention is described in further details.At this, illustrative examples of the present invention and explanation thereof are used for explanation the present invention, but not as a limitation of the invention.

In order to adapt to the real-time processing of microearthquake data, obtain fast and accurately the pressure break result, and then recovery scheme in time adjusted, present embodiment provides a kind of microearthquake event recognition method based on shape filtering, thereby provides accurately basic data for follow-up focus location.As shown in Figure 1, the method comprises:

Step S101: the microearthquake data are carried out spectrum analysis.

Step S102: the wave detector that is used for receiving the microearthquake signal is proofreaied and correct.

Step S101 and step S102 belong to microearthquake pretreatment step.The microearthquake data that field acquisition is arrived are the same with conventional earth's surface seismic section data, all will pass through first pre-service, namely the road editor, remove useless road, gain recovery, then just can carry out follow-up data and process.But because microearthquake signal self, its preconditioning technique is different from conventional disposal route, and the singularity of himself is arranged.

By step S101 the microearthquake data are carried out spectrum analysis, do bandpass filtering treatment, filter signal band noise in addition.

Fig. 2 is the frequency spectrum of three component microearthquake signals collecting in certain exploration area, and as can be seen from the figure, the microearthquake spectrum data signal that the present invention will process is probably at 10 ~ 600Hz; The microearthquake signal that is higher than 600Hz may be since ripple when in the stratum, propagating by the formation absorption loss.After knowing the spectral range of signal, can utilize bandpass filtering that the microearthquake data are processed.Fig. 3 is one microearthquake data before and after the bandpass filtering, and signal to noise ratio (S/N ratio) had obviously improved after contrast can be found out bandpass filtering treatment.

The wave detector of step S102 is proofreaied and correct and is comprised: comprising: the signal correction of the x component of wave detector on component direct P wave polarization direction, is corrected to the y component on the direction perpendicular to component direct P wave polarization direction.Detail is described below:

Be used for receiving the three-component seismometer of microearthquake signal, in its process of in well, transferring, although can make the Z perpendicular direction downward, but the orientation of two horizontal components that can not guarantee each wave detector is consistent, be actually at random, this microearthquake signal amplitude that has just caused each horizontal component on the different wave detectors to receive differs greatly, as shown in Figure 4, as can be seen from the figure, the different wave shape of horizontal component is very large in the original three-component microearthquake data, does not basically see any lineups.Therefore, before subsequent treatment, must be rotated the location to wave detector, make the horizontal component orientation of each wave detector consistent.

Concretism is: because source position and the speed of the microearthquake signal that perforating and fracturing excites are known, and energy comparison is large, signal to noise ratio (S/N ratio) is high, thereby first Seismic Direct Wave that just can utilize perforating and fracturing to excite, obtain the polarization direction of the through compressional wave of first arrival that each wave detector receives and the angle theta of wave detector x component according to the energy of ripple, the microearthquake signal that then 2 horizontal components of each wave detector is received rotates θ corresponding to each wave detector.

The wave detector horizontal component is carried out directional correction can only utilize the perforation data to ask for polarization angle.So-called perforation data refer to utilize perforating gun to penetrate the seismic event that sleeve-penetrating excites in the down-hole, and its energy comparison is strong, and first arrival is obvious.Utilizing the perforation data to carry out the horizontal component orientation is based upon on the prerequisite of such hypothesis: first component direct P ripple that passes to well seismometer from P ripple focus, its particle movement direction is consistent with direction of wave travel, all in the plane of being determined by focus and well, this component direct P wave polarization is linear, and its projection in surface level also is straight line.According to this hypothesis, just can be with the projection of component direct P wave polarization direction in surface level as a reference, measure the relative orientation of three axle wave detectors when observation horizontal component, and the signal of the horizontal component of observation is transformed into the consistent coordinate system that be projected as reference of component direct P wave polarization direction in surface level.Figure 5 shows that the orientation sketch of microearthquake observation three-component seismometer, the orientation of the horizontal component of wave detector is with degree of depth random variation, but the horizontal projection H of first P ripple _pRoughly constant.The horizontal component (x, y) of observation is transformed into H _pCoordinate projection system (x', y') for reference.Conversion formula is:

x'=xcosθ+y sinθ

y'=-xsinθ+ycosθ（1）

In the formula, θ is the angle of x and Hp, is called again component direct P wave polarization angle.

The energy criteria method is to ask for a kind of method of the angle of the polarization direction of the through compressional wave of first arrival and wave detector x component when making energy reach maximum.The time sequentially get a pair of sample value (x in the window _i, y _i), wherein i for from the time window starting point sequence number compiled of window terminal point then.If the expression formula of energy is:

$E\left(\mathrm{\θ}\right)=\underset{i}{\mathrm{\Σ}}{(x_i\cos \mathrm{\θ}+y_i\sin \mathrm{\θ})}^2---\left(2\right)$

When ENERGY E (θ) when getting maximal value corresponding angle be exactly the polarization angle of wave detector.

Fig. 6 a is one section horizontal component of perforation data, and Fig. 6 b is the hodograph that draws according to Fig. 6 a, and Fig. 6 c is the curve that draws according to energy criteria, maximal value is corresponding be exactly angle be exactly desired polarization angle.When energy is maximum as we can see from the figure in corresponding angle (about 74 degree) and the hodograph angle corresponding to the great direction of amplitude roughly consistent.

After obtaining polarization angle, according to formula (1) with seismologic record rotation of horizontal component θ angle, Fig. 7 a is the horizontal component behind the rotation θ angle, and Fig. 7 b is the hodograph of making according to horizontal component after the rotation, and the record of x' component just represents direct wave horizontal component Hp as seen from the figure.

Obtain the angle of 12 wave detectors according to the energy criteria method after, according to transformation for mula horizontal component is transformed into take component direct P wave polarization direction in the coordinate system of reference.Fig. 8 and Fig. 9 are respectively before the conversion and the hodograph synoptic diagram of making according to one section microearthquake data after the conversion, can find out that the horizontal component of conversion post-detection device has been corrected to same direction.

Step S103: the selecting structure element, will open computing and combine with closure operation and form the morphology noise filter, the described geological data through spectrum analysis is carried out shape filtering.

After the wave detector of the spectrum analysis of completing steps S102 and step S103 is proofreaied and correct, just can carry out the shape filtering operation.Below introduce in detail.

(1) mathematical morphology fundamental operation

The theoretical foundation of shape filtering is mathematical morphology.Mathematical morphology is that a kind of Nonlinear image processing take set theory and logical operation as the basis is theoretical, its basic thought is with the structural element with certain form target image to be carried out conversion, measure and extract wherein correspondingly-shaped, the simplified image data, keep the basic style characteristic of image, get rid of incoherent graphic structure, thereby reach the purpose of shape filtering.

The basic transformation of the mathematical morphology of bianry image has four kinds: i.e. burn into expansion, opening and closing, they can be combined into complicated image processing techniques.If A is the set of target image; B is the set of structural element, then:

Erosion operation is defined as: $\mathrm{A\ΘB}=\underset{b\∈B}{\mathrm{\∪}}A-b$

Dilation operation is defined as: $A\⊕B=\underset{b\∈B}{\mathrm{\∪}}A+b$

The unlatching operation definition is:

Closure operation is defined as: $A\·B=(A\⊕B)\mathrm{\ΘB}$

The unlatching computing is adopted to corrode afterwards first and is expanded, and plays a part separation, filtering, and the isolated part less than structural element all will be filtered, with Inhibitory signal peak value (positive pulse) noise; Closure operation adopts the post-etching that expands first, play fill a vacancy, the internal communication effect, with Inhibitory signal bottom-valley (negative pulse) noise, Figure 10 a is the corrosion signal schematic representation that utilizes the semicircular structure element, and Figure 10 b is the expansion signal schematic representation that utilizes the semicircular structure element.

The gray scale opening and closing also can simply be used geometric interpretation, and we discuss by Figure 11.Open f with b, i.e. f zero b can regard the lower edge that b is being pasted f as and be rolled into the other end from an end.Figure 11 (b) provides the several positions of b in unlatching, and figure c provides the result of open operation.Can find out from Figure 11 (c), all mountain peak its height and sharpnesses less than the diameter of b have all been weakened.In other words, when b pasting f lower when rolling, the position that does not contact with b among the f is all fallen with b and is contacted.Open operation commonly used is eliminated and is compared the less bright details of size with structural element in the reality, and keeps figure overall gray value and large bright area substantially unaffected.The erosion removal in concrete the 1st step little bright details and weakened simultaneously brightness of image, brightness of image that the expansion of second step has increased (the basic recovery) but again do not introduce the details of removing previously.With b closed f, i.e. f ● b, can regard the upper edge that b is being pasted f as and be rolled into the other end from an end.Figure 11 (d) provides the several positions of b in closure, and Figure 11 (e) has provided the result of closed procedure.Can find out that from Figure 11 (e) mountain peak does not have to change substantially, and all mountain valleys less than the diameter of b have obtained filling.In other words, when b pasting f upper when rolling, the mountain valley that does not contact with b among the f all is filled into b and contacts.Closed procedure commonly used is eliminated and is compared the less dark details of size with structural element in the reality, and keeps integral image gray scale and large dark areas substantially unaffected.Particularly, the expansion in the 1st step has been removed little dark details and has been strengthened simultaneously brightness of image, brightness of image that the corrosion in the 2nd step has weakened (the basic recovery) but again do not introduce the details of removing previously.

(2) mathematical morphology filter principle

Different from digital filtering with the shape filtering processing seismic data, the algorithm of shape filtering is not the algorithm of fixing, and decide along with the structural element of selecting, and the unlatching computing can filter the bur less than structural element, cuts off elongated overlap joint and plays the effect of separation.Closure operation can be filling than the little breach of structural element or hole, the interruption that overlap joint is short and play the connection effect.Because noise kind quasi-complexity in the geological data often needs the complex morphological filtering of many structural elements.

1. array mode is on the impact of filter effect

The anti-extendability of (image of doing opening operation always is positioned at the below of original image) and closed operation because the non scalable of opening operation (image of doing closed operation always is positioned at the top of original image), all there is the statistical bias phenomenon in two kinds of wave filters, this moment, open-close wave filter output amplitude was less than normal, and close-output amplitude of Kai wave filter is bigger than normal, in most situation, use separately the filter effect that they can not obtain, preferably adopt the average combined form of two kinds of wave filters.If input signal is

x(k)=s(k)+n(k) (k=1，...,N)

S in the formula (k) is original signal; N (k) is noise.Then filtering output signal y (k) is

y(k)=[OC(x(k))+CO(x(k))]/2

OC represents to carry out first opening operation in the formula, carries out closed operation again, and CO represents to carry out first closed operation and carries out opening operation again.

If original signal expression formula: y (t)=e ^50tSin2 π ft(is wherein: f=30, t=0～0.1 is spaced apart 0.0001), Figure 12 a, Figure 12 b, Figure 12 c and Figure 12 d provided respectively switching, make and break, the two on average and result schematic diagram of contrast.

2. the selection of structural element in the seismic processing

Filter effect is relevant with structural element, obtain desirable filter effect, need to be for problem choose reasonable structural element to be solved.The selection of structural element comprises the width (width of the field of definition of structural element) that determines structural element, highly (amplitude of structural element) and shape.The below will provide our result in the research of choice structure element.

The width of structural element is on the impact of filter effect:

The wavelet expression of test usefulness: y (t)=e ^-50t Sin 2 π ft

In the following formula: f=40, t=0～0.1 is spaced apart 0.0001.

Figure 13 is seismic wavelet model synoptic diagram, and the original signal of expression wavelet model is divided on the top of Figure 13, and the bottom divides the expression original signal to add random noise, s/n (useful signal energy/random noise energy)=57.403558.

Among Figure 14, it (d) is to be 2,5,10,40 filtered results with the flat-structure element width respectively that (a) and (b), (c) reach, and filtered s/n is followed successively by 108.144122,196.832655,260.819392,142.8286.

As seen from Figure 14, through behind the morphologic filtering, random noise is by filtering, and signal to noise ratio (S/N ratio) obviously increases, and the size of structural element is very large on the result impact.(c) filter effect is better than (b), (b) effect is better than (a), this is because their size of used structural element increases gradually, (c) structural element is of a size of 10 sampling points, (b) structural element is of a size of 5 sampling points, and structural element (a) is of a size of 2 sampling points.

But be not that the size of structural element is the bigger the better.The size of structural element is larger, and on the one hand, calculated amount can increase, and on the other hand, it may affect useful signal, sees (d).So, should be according to the situation of noise and useful signal the size (being number of samples) of suitable selecting structure element, size of structure element is too little, noise goes unclean, excessive, can damage useful signal again, so will be according to the size of the signal characteristic selecting structure element that will keep after the filtering.

For impulsive noise, be T if treat filtering impulsive noise breadth extreme, the sampling period is T _s, experimental result shows: the length M of structural element only need be a bit larger tham T/T _sGet final product.

The height of structural element is on the impact of filter effect:

After having considered the impact of width on filter effect of structural element, see again which type of impact is the amplitude of structural element have to morphologic filtering.

The wavelet expression of test usefulness: y=e ^-50tSin2 π ft

In the formula: f=40, t=0～0.1 is spaced apart 0.0001, and used structural element is the semi-circular structure element.

Figure 15 is seismic wavelet model synoptic diagram, and the original signal of expression wavelet model is divided on the top of Figure 15, and the bottom divides the expression original signal to add random noise, s/n (useful signal energy/random noise energy)=56.162128.

Among Figure 16, (a), (b), (c), (d) use respectively semi-circular structure element (0,0.0017,0.0020,0.0017,0), (0,0.0173,0.0200,0.0173,0), (0,0.1732,0.2000,0.1732,0), (0,1.7321,2.0000,1.7321,0) and filtered result, its radius is followed successively by 0.002,0.02, and 0.2,2, filtered s/n is followed successively by 267.392131,295.398691, and 216.562129,5.797541.As can be seen from Figure 16, impact is not very large to the amplitude of structural element on result, still, if the radius of structural element is large, so after the filtering, in beginning and ending place of signal, it is larger to distort, the waveform of signal after the s/n behind the integrated filter and the filtering, and filter effect (b) is quite a lot of.

Experimental result shows: the amplitude of structural element than signal value to be filtered to order of magnitude when young.

The shape of structural element is on the impact of filter effect:

Here we discuss four kinds of structural elements: linear, triangle, semicircle, four kinds of structural elements of sinusoidal are on the impact of result.

The wavelet expression of test usefulness: y=e ^-50tIn the sin 2 π ft formulas: f=30, t=0～0.1 is spaced apart 0.0001.

Figure 17 is seismic wavelet model synoptic diagram, and the original signal of expression wavelet model is divided on the top of Figure 17, and the bottom divides the expression original signal to add random noise, s/n (useful signal energy/random noise energy)=120.4241.

Figure 18 is with sinusoidal, semicircle, linear and four kinds of structural elements of triangle are suppressed random noise respectively to Figure 21.(a) expression sinusoidal structural element among Figure 18, (b) signal behind the expression sinusoidal structural element shape filtering; (a) expression semicircular structure element among Figure 19, (b) the filtered signal of expression semicircular structure element morphology; (a) expression flat-structure element among Figure 20, (b) signal behind the expression flat-structure element shape filtering.(a) expression triangular structure element among Figure 21, (b) the filtered signal of expression triangular structure element morphology.To the structural element of these several shapes, we have used two kinds of different sizes to carry out the filtering test, select the effective shape matching that carries out.Through the filtering of these several structural elements, noise obviously is pressed in the signal, and signal to noise ratio (S/N ratio) significantly improves, filtered signal to noise ratio (S/N ratio) difference: 1.0755e+003,1.0488e+003,661.4504,846.2939.

Can be found out that by experimental result in general, the shape of structural element will be decided according to the shape of the signal that will keep after the filtering, be used for the optional semicircle of structural element that Seismic signal filtering processes, sinusoidal etc. and the similar shape of waveform.

Owing to there are differences at structural form between signal and the noise, therefore can by choose suitable structural element, opening and closing be combined construction form noise filtering device microseismograms is processed.Figure 22 is pretreated one microearthquake data, and Figure 23 is to its closed through opening, closed data of opening after processing.Find out that by figure shape filtering has filtered the part noise, but do not change whole form, that is to say that the waveform of useful signal is substantially constant, the amplitude of useful signal, the feature of frequency are kept well.The method has higher actual application value.

Step S104: according to the energy Ratios that carries out data behind the shape filtering and calculate the length time window, described energy Ratios and activation threshold value are compared to identify described wave detector whether have the microearthquake event.

The energy spectrometer technology is most widely used a kind of automatic identification technology.Most of energy spectrometer is based on energy ratio between the long window and the short window method (STA/LTA) and realizes, he according to different components, different coordinates or not simultaneously the difference of the energy Ratios of window distinguish effective seismic signal.Because microearthquake event and ground unrest are different on energy feature, this is just for providing possibility in three-component seismogram identification microearthquake event.In for the situation of many components microearthquake data, can select suitable component to calculate according to the actual conditions of data.For one three-component record (x (t), y (t), z (t)), z (t) is vertical component, x (t), and y (t) is horizontal component, can be used for asking the combination of energy to have:

$E_{\mathrm{all}}\left(t\right)=\sqrt{z{\left(t\right)}^2+x{\left(t\right)}^2+y{\left(t\right)}^2}$

E _Z(t)=z(t)| （3）

$E_H\left(t\right)=\sqrt{x{\left(t\right)}^2+y{\left(t\right)}^2}$

These three kinds combinations are non-not to be three-component gross energy and, vertical component amplitude and horizontal component energy.

Window when in application process, select rolling, i.e. t at a time ₀, window when choosing with certain-length in its front and back, the mean value of window self-energy when asking for utilizes these the average energy value calculating energy ratios.

Energy ratio between the long window and the short window method (STA/LTA) fundamental function can be defined as:

R=STA(t)LTA(t) （4）

Wherein

$\mathrm{LTA}\left(t\right)=\underset{t=t_1}{\overset{t_0}{\mathrm{\&Sigma;}}}E\left(t\right)/|t_0-{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="DEST\_PATH\_HDA00002805414700021.png,DEST\_PATH\_HDA00002805414700151.png"\; state="\{\{state\}\}">t</figure-callout>}}_1|---\left(5\right)$

$\mathrm{STA}\left(t\right)=\underset{t=t_0}{\overset{t_2}{\mathrm{\&Sigma;}}}E\left(t\right)/|t_2-{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="DEST\_PATH\_HDA00002805414700061.png,DEST\_PATH\_HDA00002805414700071.png"\; state="\{\{state\}\}">t</figure-callout>}}_0|---\left(6\right)$

t ₁Window initial time during for length, t ₀Window stops constantly and the short time-window initial time during for length, t ₂For short time-window stops constantly.

Signal utilizes energy to portray the feature of microearthquake event in long short time-window, the variation tendency that window mean value (LTA) has been portrayed the signal background noise when then utilizing fundamental function to slide long short time-window calculating STA and LTA. length, short time-window mean value (STA) has been portrayed amplitude (energy) variation tendency of microearthquake signal, when signal arrives, STA changes soon than LTA, corresponding STA/LTA value has an obvious kick, when the ratio of short time-window mean value window mean value when long surpasses user-defined activation threshold value, just be judged to be the microearthquake event.Can be set as a fixed value about threshold value, in the situation of estimating microearthquake Monitoring Data signal to noise ratio (S/N ratio), usually can be made as 3/4 of signal to noise ratio (S/N ratio), under this thresholding, substantially can identify the microearthquake validity event that major part can be utilized.

Figure 24 is the three-component energy Ratios curve synoptic diagram of certain road wave detector.Because this automatic identifying method is based on the method for single-stage wave detector, and for the microseismograms of observing environment and low signal-to-noise ratio in the well of complexity, the signal that records on the single utilization one-level wave detector judges whether to exist the microearthquake useful signal mistake often to occur pick up situation.Thereby will utilize the aggregation of data on the multi-component and multilevel geophone to carry out the microearthquake automatic identification of effective time.If certain time detecting at the single-stage wave detector arrives the microearthquake signal, on other wave detectors, should also signal can occur so, therefore can open a fenestella at other wave detectors according to the occurrence law of effective microearthquake signal between the phase neighboring trace and detect.Just think signal if on multi-component and multilevel geophone, all detect signal, can avoid like this mistake that the very noisy interference causes on the single-stage wave detector to pick up.

In the situation that has default threshold value to be triggered, be merely able to judge existing effect event, but cause that the signal of triggering is that P ripple or S ripple just can't be judged.The validity event that identifies for judgement is P ripple or S ripple, comes these two kinds of ripples are distinguished with regard to the polarization direction that needs to analyze useful signal.

The microearthquake source location is unknown, but can determine according to certain prior imformation the spatial dimension of focus.If there is no medium anisotropy or the impact that produces shear wave splitting owing to the crack, polarizability between P ripple, SV ripple, this three of SH ripple is mutually perpendicular relation, and their polarization direction can be decomposed its maximal eigenvector of searching according to covariance matrix and be tried to achieve.

If the time window that N sampling point arranged, each sampling point is by three-dimensional coordinate X, Y, the definition of Z three-component.The time window [T1, T2] in, the mean value of sampling point coordinate is

$\stackrel{\&OverBar;}{z}=\frac{1}{N}\underset{i=n_1}{\overset{n_2}{\mathrm{\&Sigma;}}}{z}_i$ $\stackrel{\&OverBar;}{y}=\frac{1}{N}\underset{i=n_1}{\overset{n_2}{\mathrm{\&Sigma;}}}{y}_i$ $\stackrel{\&OverBar;}{x}=\frac{1}{N}\underset{i=n_1}{\overset{n_2}{\mathrm{\&Sigma;}}}{x}_i---\left(6\right)$

(n in the formula ₂-n ₁) Δ t=T ₂-T ₁, Δ t is sampling interval; N=n ₂-n ₁The+1st, so sampling number, the covariance matrix of sampled point is

$M=\left[\begin{array}{ccc}\mathrm{cov}[x,x]& \mathrm{cov}[y,x]& \mathrm{cov}[z,x]\\ \mathrm{cov}[x,y]& \mathrm{cov}[y,y]& \mathrm{cov}[z,y]\\ \mathrm{cov}[x,z]& \mathrm{cov}[y,z]& \mathrm{cov}[z,z]\end{array}\right]---\left(7\right)$

Wherein cov represents the covariance computing, and the computing of asking for the covariance of vector x (i) and y (i) is:

$\mathrm{cov}[x,y]=\frac{1}{N}\underset{i=n_1}{\overset{n_2}{\mathrm{\&Sigma;}}}(x_i-\stackrel{\&OverBar;}{x})\&CenterDot;(y_i-\stackrel{\&OverBar;}{y})---\left(8\right)$

Can obtain three eigenvalue λ by variance matrix ₁, λ ₂, λ ₃(λ ₁λ ₂λ ₃) and three proper vector Λ ₁, Λ ₂, Λ ₃The maximal eigenvector that eigenvalue of maximum is corresponding can be thought the direction of ripple vibration.

If there is no medium anisotropy or the impact that produces shear wave splitting owing to the crack, the polarizability between P ripple, SV ripple, this three of SH ripple is mutually perpendicular relation.P ripple direction of vibration is identical with the seismic ray direction of propagation, also just can think, vertical component is very little in the P ripple three-component record that the wave detector in the stratum that approaches with the focus distributed depth receives, and polarization direction and wave detector z component are near vertical.Identical with the vibrations plane of P ripple, but more be partial to the vertical direction of wave detector for the focus SV wave polarization direction of same position.P ripple if fractured layer position and the degree of depth of wave detector are staggered on the record, SV wave polarization can with opposite in the degree of depth situation.For the SH ripple, the plane of oscillation of its primary wave and plane, place, the direction of propagation and P wave polarization direction are orthogonal, and vertical component can not occur on the wave detector under stratified model.Before carrying out the validity event automatic Picking, according to the job design of hydraulic fracturing and the locus of monitor well fractured well, roughly estimate the scope of microearthquake focus, can estimate P ripple and SH wave polarization direction, in general these two the scope that overlaps or overlap can not occur very little.

Step S105: carry out microearthquake according to described microearthquake event and pick up.

Step 5 specifically comprises:

Step 501: utilize the erosion algorithm in the mathematical morphology can realize removing fast random noise, thus further first break picking.First the amplitude of microearthquake signal is taken absolute value allow its all become on the occasion of, make it meet the requirement of erosion algorithm.

Step 502: choose sizeable semicircular structure element, after described microearthquake data take absolute value behind the shape filtering, carry out corrosion treatment, the signal schematic representation after Figure 25 road original signal and the corrosion treatment.

Step 503: get the mean value of energy value of the signal after the corrosion treatment as the first arrival activation threshold value.Through a large number of experiments show that, average and more accurately to trigger primary wave as adaptive threshold, and can not trigger interference noise.Figure 26 corrodes rear signal energy value synoptic diagram, and the solid horizontal line of lowermost end is activation threshold value.

Step 504: get forward window when sizeable at activation threshold value point place, Figure 27 is signal energy figure after the local corrosion of amplifying, and for triggering moment point, the dotted line frame is the time window of getting to t0 constantly.

Step 505: choose sizeable long short time-window, when described in the window take a sampled point as step-length, carry out the energy ratio between the long window and the short window computing;

Step 506: the point of determining the energy ratio maximum is the first arrival point.Utilize the result of automatic Picking first arrival just can carry out the microquake sources location, mainly use the method for time difference method and ray casting co-located.

Below in conjunction with concrete example, microearthquake event recognition method of the present invention is described.

(1), bandpass filtering treatment

We have known the general frequency range of signal according to spectrum analysis, and we are set as 10-600Hz with frequency range when carrying out bandpass filtering.Figure 28 is one section untreated original microearthquake data, and Figure 29 is through the data after the bandpass filtering treatment.As seen from the figure, bandpass filtering has been removed the part random noise, and the signal ratio of data is improved significantly.

(2), the location of horizontal component

Because the orientation of two horizontal components of each wave detector is at random in the process that the three-component seismometer of reception microearthquake signal is transferred in well, before subsequent treatment, must be rotated the location to wave detector, make the horizontal component orientation of each wave detector consistent.Figure 30 is the geological data after corresponding diagram 4 conversion, because the wave detector horizontal component has been corrected on the direction, energy distribution is more even on wave detectors at different levels, and the microearthquake signal that obtains like this is just more continuous, simultaneously signal is enhanced.

(3) shape filtering is processed

In image processing and signal analysis, the method for smothing filtering has a lot, such as medium filtering and neighbour average filtering etc.Here we discuss, non-linear form filtering, and non-linear medium filtering, linear neighboring mean value is to the effect of random noise and impulsive noise compacting.

The wavelet expression of test usefulness: y=e ^-50tSin2 π ft

In the formula: f=30, t=0～0.1 is spaced apart 0.0001.

Figure 31 is seismic wavelet model synoptic diagram, and the original signal of expression wavelet model is divided on the top of Figure 31, and the bottom divides the expression original signal to add random noise, s/n (useful signal energy/random noise energy)=107.6534.

Figure 32 to Figure 34 is that three kinds of wave filters are to the effect synoptic diagram of random noise compacting.Signal to noise ratio (S/N ratio) after curve and processing can find out that three kinds of wave filters have some good effects, and random noise is pressed, and signal to noise ratio (S/N ratio) is improved, but on the whole, the morphological filter effect is better than median filter and neighbour average filtering device.

Figure 35 to Figure 38 is the pressing result figure of three kinds of wave filter paired pulses noises.Can find out that from filtered waveform and signal to noise ratio (S/N ratio) when removing impulsive noise, morphologic filtering method is better than other two kinds.Above experimental result shows that morphologic filtering can be used for the denoising of geological data fully.

It more than is the denoising effect that shape filtering is discussed with model, we illustrate the effect of shape filtering with real data now, Figure 39 is pretreated one section microearthquake schematic diagram data, can see obvious strong vibration noise is arranged on the first wave detector, this may be because wave detector does not push against on the borehole wall or backup power is inadequate, even to collapse appear in the borehole wall, so that wave detector and stratum or sleeve pipe coupling are bad, thereby produces stronger resonance.(Figure 40) can see after utilizing shape filtering to process, and strong energy resonance is obviously suppressed, thereby the explanation shape filtering can reach the purpose of denoising preferably, also can play the effect of strengthening weak signal.

The detection of (four) microearthquake event

After processing through filtering, utilize energy ratio function can identify fast and accurately the microearthquake signal, thereby be used for actual production.For the energy ratio function advantage is described, we carry out the identification of microearthquake event with energy ratio function and polarization coefficient method respectively.

Figure 41 is the energy Ratios curve of one section microearthquake data, Figure 42 is corresponding polarization coefficient curve, energy ratio function has maximum value at the signal place as seen from the figure, and P ripple, S wave energy enough separate, although maximum value on the polarization coefficient curve, also occurs, but larger error (the 2nd, 5,12 road) but appears on some seismic trace, and if two signal times poor when smaller, polarization coefficient can not separate identification fully.Arithmetic speed is convenient, energy ratio function also ratio polarization coefficient wants fast many, the microearthquake data of one section 10s of same operation, be Intel (R) Core (TM) Duo CPU E4500@2.20GHz2.20GHz at processor, in to save as on the computer of 1.00GB energy ratio function working time be 0.546000s, polarization coefficient working time is 1.822000s, and be 3.3 times of energy ratio function working time.Therefore energy ratio function can be identified the microearthquake event fast and accurately.In addition, utilize energy ratio function need to utilize the multiple tracks wave detector to carry out recognition detection.Can be seen if the microearthquake signal that by Figure 41 its can show at the multiple tracks wave detector (red frame place) rather than peak value occur together at certain, we just can avoid some to miss the possibility of picking up like this.And the identification of waveform can utilize polarization analyzing method, obtains the direction of the polarization of microseism, can judge more exactly the classification of waveform, is P ripple or S ripple.

(5) microearthquake first arrival automatic Picking

The precision of picking up in order to improve first arrival, need to carry out pre-service to the microearthquake data, thereby shape filtering can effectively be removed random noise and improve signal to noise ratio (S/N ratio), the method of tradition first break picking is to adopt fixed threshold, do not have the easy erroneous trigger of adaptability, therefore according to the characteristics of microearthquake data self, with the average of the Amplitude-squared value of pretreated data as adaptive threshold, can correctly trigger primary wave, improve and pick up precision.And window carries out the energy ratio between the long window and the short window computing when getting in part, trigger point place, has improved arithmetic speed.

Figure 43 is the first arrival in actual microearthquake data road and picks up the result, and as can be seen from the figure the energy ratio curve has an obvious peak value, and the time t1 that this peak value is corresponding is initial time.

(6) microearthquake focus location

Figure 44 is the theoretical model based on the time difference method location, " * " represents the wave detector position, and " o " represents actual source location, and the gained source location is calculated in "+" representative, calculate focus and coincide with actual source location, illustrate that the method picks up in the higher situation of precision the location in first arrival and stablize accurate.Figure 45 is the ray tracing location algorithm theoretical model based on polarization analysis, and " * " is the wave detector position, and " o " is the best convergence point of beam, is the source location, by the result as can be known the method in the obvious situation of P ripple first arrival, also can accurately locate.With a kind of method location the larger situation of error can appear separately, constrained each other if therefore two kinds of methods are joined together the location, can effectively reduce positioning error.

The beneficial effect of the embodiment of the invention is, the present invention processes the signal to noise ratio (S/N ratio) that can improve data by pre-service and further filtering, utilize energy ratio function can detect rapidly and accurately the microearthquake event, and identify the fluctuation type of microearthquake event by polarization analysis, utilize automatically accurate first break picking of the adaptive threshold of design and energy ratio between the long window and the short window method.

Above-described specific embodiment; purpose of the present invention, technical scheme and beneficial effect are further described; institute is understood that; the above only is specific embodiments of the invention; the protection domain that is not intended to limit the present invention; within the spirit and principles in the present invention all, any modification of making, be equal to replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (6)

1. microearthquake event recognition method based on shape filtering is characterized in that described method comprises:

Step 1: the microearthquake data are carried out spectrum analysis;

Step 2: the wave detector that is used for receiving the microearthquake signal is proofreaied and correct;

Step 3: the selecting structure element, will open computing and combine with closure operation and form the morphology noise filter, the described geological data through spectrum analysis is carried out shape filtering;

Step 4: according to the energy Ratios that carries out data behind the shape filtering and calculate the length time window, described energy Ratios and activation threshold value are compared to identify described wave detector whether have the microearthquake event;

Step 5: carry out microearthquake according to described microearthquake event and pick up.

2. method according to claim 1, it is characterized in that, the wave detector that is used for receiving the microearthquake signal is proofreaied and correct, being comprised: the signal correction of the x component of wave detector on component direct P wave polarization direction, is corrected to the y component on the direction perpendicular to component direct P wave polarization direction.

3. method according to claim 2, it is characterized in that, with the signal correction of the x component of wave detector on component direct P wave polarization direction, the y component is corrected on the direction perpendicular to component direct P wave polarization direction, comprise: obtain the polarization direction of the through compressional wave of first arrival that each wave detector receives and the angle theta of wave detector x component according to the energy of ripple, the microearthquake signal that then 2 horizontal components of each wave detector is received rotates θ corresponding to each wave detector.

4. method according to claim 1 is characterized in that, described selecting structure element comprises: the width of selecting structure element, height and shape.

5. method according to claim 1, it is characterized in that, for multi-component and multilevel geophone, described step 4 also comprises: open a fenestella according to effective geological data of the microearthquake event of identifying at the wave detector of the rule that occurs between the phase neighboring trace outside described wave detector and detect.

6. method according to claim 5 is characterized in that, described step 5 comprises:

Utilize erosion algorithm to remove random noise;

Choose sizeable semicircular structure element, after described microearthquake data take absolute value behind the shape filtering, carry out corrosion treatment;

Choose the mean value of energy value of the signal after the corrosion treatment as the first arrival activation threshold value;

Get forward window when sizeable at activation threshold value point place;

Choose sizeable long short time-window, when described in the window take a sampled point as step-length, carry out the energy ratio between the long window and the short window computing;

The point of determining the energy ratio maximum is the first arrival point.

Images