RU2463627C1 - Method of monitoring region with network of seismic stations - Google Patents

Abstract

FIELD: physics.

SUBSTANCE: network of seismic stations which use pressure sensors as primary velocity sensors are placed in an earthquake-prone region. The ratios of velocities of longitudinal and transverse waves in the earth's crust are determined at different azimuth angles of the route of each of the network of stations. Discrete samples of measurements are obtained for each of the routes in the same time intervals. The Fourier spectra of the obtained samples are determined with determination of weighted average frequencies. The moment of change in weighted average frequencies is identified with the occurrence of an anomalous zone of an imminent earthquake on the measured route. Bearings of the anomalous zone are taken as points of intersection of the radius vectors of velocities of several stations on the measured routes. Signal energy of the calculated spectra is determined and the behaviour of its variation is monitored. The size of the anomalous zone is calculated. The time and magnitude of the expected seismic impact are predicted.

EFFECT: longer time interval for preventive forecast, high accuracy of determining forecast parameters of a seismic impact.

6 dwg

Description

The invention relates to the field of seismology and may find application in national systems for observing and processing data from geophysical measurements for predicting earthquakes.

A reliable forecast of earthquakes is possible when establishing their root cause. One of the root causes of earthquakes is the pumping of the earth's crust by the additional energy of an upward flow of gases (H ₂ , He, CH ₄ ). The fact of the emanation of gases from the earth's crust on the eve of a seismic shock [see, for example, Patents RU No. 2204852, 2003; No. 2275659, 2006; No. 2302020, 2007]. Pumping the earth's crust with additional energy leads to a buildup of the earthquake source, accompanied by the propagation of ultra-low lithospheric waves from the hypocenter [see Scientific discovery No. 365. “The phenomenon of buildup of the earthquake source before a seismic impact”, Russian Academy of Natural Sciences, 2008].

The famous "Method for predicting earthquakes." Patent RU No. 2170446, 2001, by measuring the ultralow lithospheric waves of the buildup of the earthquake source by space means of the Navstar (GPS) system - an analogue.

In the analogous method, receiving stations of the space navigation system are placed in a seismic region, spaced on an extended measuring base, they perform continuous high-precision measurement of the coordinates of the receiving points of the receiving stations, register the moment of occurrence of periodic deviations Δx _i , Δy _i , Δz _{i of the} coordinates of the points and track changes in these deviations in time, the hypothetical center of the focus is calculated as the point of intersection of the radius vectors in space, the length and the guides of the cosines, which are determined from the relations

calculate the impact time, counted from the moment of the appearance of periodic deviations of coordinates:

and magnitude

where T is the coordinate deviation period, hour, dek r is the natural logarithm of the ratio of the coordinate deviation amplitudes of two adjacent periods; d, I - regression coefficients.

The disadvantages of the analogue method are:

- the presence of a hidden section of insensitivity (dead zone), limited by the mean square error of the GPS means, reducing the forecast interval;

- inaccuracy of the regression dependences of calculating the time of impact and magnitude due to the dependence of the period (T) of the coordinate deviation from time.

The closest analogue to the claimed method is the "Method of short-term earthquake prediction." Patent RU No. 2181205, G01V, 9/00, 2002

In the closest analogue method, the seismic background is recorded in the form of a continuous sequence of discrete samples of the signal amplitude Δ (t) at two points spaced along the coordinates, the Fourier spectrum of the sequence of measurement samples with the report size in each sample N≥2F _max / σ is calculated, autocorrelation functions are calculated In (τ) the signals of the samples, the correlation interval τ is determined, the onset of the change in the parameter τ is recorded, and when it is continuously monitored, the delay time Δτ of the phase change of this sign between the two points is calculated traverse direction cosine wave arrival hearth ultralow

determine the hypothetical center of the focus as the point of intersection on the beam of the radius vectors of points with the cosine of the angle at the apex:

(час) и время удара t_x=2,3Т₀, где F_mах - максимальная частота спектра сейсмического фона, Гц; σ - среднеквадратическая ошибка вычисления спектра Фурье по дискретной выборке измерений; α - длина базы между двумя пунктами, м; B₁(0), В₂(0) -значения автокорреляционных функций в нуле для каждого пункта; V -скорость сейсмических волн в земной коре, м/с.the period T _{0 of the} parameter τ is calculated and magnitude is predicted by its magnitude

(hour) and impact time t _x = 2,3Т ₀ , where F _max is the maximum frequency of the spectrum of the seismic background, Hz; σ is the standard error of the calculation of the Fourier spectrum from a discrete sample of measurements; α is the length of the base between two points, m; B ₁ (0), B ₂ (0) -values of autocorrelation functions at zero for each item; V is the speed of seismic waves in the earth's crust, m / s.

The disadvantages of the closest analogue include:

- the autocorrelation interval of the signal τ is determined primarily by the width of the spectrum of the signal of the seismic background, therefore, the measurement of the parameter Δτ is technically difficult to implement, since it practically does not depend on the size of the base;

- the inaccuracy of the regression dependences of determining the magnitude and time of impact, in particular, is known from the Gutenberg – Richter relations that the longer the time T ₀ , the magnitude of the expected impact should be greater, moreover, the estimated period T _{0 is} not constant, but a function of time T ₀ (t) the buildup of the earthquake source.

The problem solved by the claimed method is to extend the time interval of the proactive forecast, increase the accuracy and reliability of determining the predicted parameters of seismic impact by existing ground-based means by extracting hidden information from the ratio of the velocities of longitudinal and transverse waves along the selected propagation paths.

The technical result is achieved by the fact that the method for monitoring a region by a network of seismic stations includes route measurements of the ratio of the propagation velocities of longitudinal V _p and transverse V _s waves in the earth's crust (V _p / V _s ) under different azimuths of the tracks of each network station using primary velocity sensors pressure sensors placed in mutually orthogonal measurement planes V _p and V _s , obtaining discrete measurement samples (V _p / V _s ) for each of the tracks at the same time intervals [V _p / V _s (t _i )], calculation Fourier spectra Paul scientific samples with the determination of the weighted average frequencies of the samples, the identification of the moment of change in the average frequencies with the occurrence of the anomalous zone of the prepared earthquake on the measured path, direction finding of the anomalous zone as the point of intersection of the velocity radius vectors of several stations on the measured paths, the determination of the signal energy of the calculated spectra as

- постоянная времени затухания сигнала на измеряемой трассе;and tracking the dynamics of its change in the form of an exponent E (t) = E ₀ * exp (-t / T), calculating the size of the anomalous zone R, km ≂ 1 / E _ust predicting the magnitude of the expected seismic impact М≂ln [R, km] and impact time, measured from the moment the weighted average frequency of the Fourier spectrum changes: t _y ~ 4.7T, where

- the time constant of the signal attenuation on the measured path;

Δt is the time interval between the calculated definitions of the signal energy E ₁ and E ₂ ;

E ₀ is the signal energy until the anomaly occurs;

E _set ≅E ₀ * е ^-4.7 is the steady-state value of the energy of the signal passing through the anomaly on the eve of the seismic shock.

The invention is illustrated by drawings, where

figure 1 - network diagram of the seismic stations of the Kamchatka geophysical test site, which worked out the claimed method and direction finding anomalies as the point of intersection of the radius of the vectors of the seismic stations;

figure 2 - implementation of the functions Vp / Vs (t) for individual tracks in the same measurement intervals;

figure 3 - Fourier spectra of samples of measurements of the monitored functions;

figure 4 - dynamics of changes in the energy of the signals on the measured paths;

5 is a functional diagram of a device that implements the method;

6 is a summary graph of the monitoring results.

The technical essence of the invention is as follows. Reliable prediction can be provided by those methods that are based on measuring the precursors of the root cause of earthquakes. One of the root causes of earthquakes is the degassing of the earth's crust and its pumping with additional energy of an upward flow of high pressure gases. Saturation of the earth's crust with a gas component changes its viscoelastic characteristics. The dispersion of rock density leads to a dispersion of the propagation velocities of lithospheric waves, i.e. a change in the relationship between the propagation velocities of longitudinal V _p and transverse V _s waves.

The propagation velocity of longitudinal waves in the medium

transverse

where E is the Young's modulus of elasticity, G is the shear modulus, ρ is the density of the medium [see Physics course G.A. Zisman, O.M. Todes, vol. I, pp. 301-302, publishing house Science. M., 1964]. As a function of the sensitivity of changes in the viscoelastic characteristics of the medium, the function V _p / V _{s is} selected, which varies in time and in space from the azimuth of the measurement path. Implementations of the function V _p / V _s before and after the seismic impacts on the measured paths are illustrated by graphs of figure 2.

It is known [see, for example, R. Duda, P. Hart. “Pattern Recognition and Scene Analysis”, trans. from English M.: Mir, 1976, pp. 319-331], that the maximum information about the ongoing processes is contained in the image of the object. The image of the object of the earth's crust is the signal spectrum of the measured function.

To obtain a spectrum image, mathematical procedures are used for the spectral analysis of recorded functions by expanding them into a trigonometric series that provides the highest display accuracy [see, for example, N. S. Piskunov. Differential and integral calculus for VTU3, textbook, vol. II, 5th edition. M .: Nauka, 1964, Fourier Series, pp. 180-182, 218-221].

рассчитывается из соотношения:By definition, the Fourier spectrum F (jω) of a function

calculated from the ratio:

t ₂ -t ₁ is the interval of sample changes. In accordance with the Kotelnikov theorem, for an adequate representation of the sample of measurements by the frequency spectrum F (jω), it is necessary that the sampling interval Δt

time function samples satisfied the condition

(F _max is the maximum frequency of the spectrum of the seismic signal), and the sample size was at least 600 samples.

When this condition is met, the smallest changes in the spectrum of the selected function will be detected.

Spectrum calculation is carried out by fast Fourier transform (FFT) algorithms according to standard programs included in the package of specialized PC software such as MATHCAD, ER MAPPER, (see, for example, specialized software MATHCAD 6.0 PLUS, 2nd stereotyped edition. M .: Information and Publishing House "Filin", 1997, p.441).

(фиг.2) до и после состоявшихся сейсмических ударов иллюстрируются графиками фиг.3.Calculated Fourier spectra of samples of functions

(figure 2) before and after the seismic impacts are illustrated by graphs of figure 3.

увеличивается, что отражается в спектре как увеличение его средневзвешенной частоты.The occurrence on the measured path of the gas-saturated volume of the earth's crust of the zone of the prepared earthquake changes the parameters of the recorded signal: the width of the spectrum, the proportions in the ratio of the amplitudes of the spectral lines, the weighted average frequency of the spectrum, and the amount of attenuation. With the beginning of the transition to seismic shock of the process, the rate of change of function

increases, which is reflected in the spectrum as an increase in its weighted average frequency.

The moment of change in the weighted average frequency of the spectrum is identified with the occurrence of an anomalous zone of the prepared earthquake on the measured path and the beginning of the seismic process.

After detecting anomalies on the measurement path, it is direction-finding. It is known (see, for example, Korn G., Korn T. "A Handbook of Mathematics for Scientists and Engineers" translated from English. M.: Science, 1971, section "Analytical Geometry", pp. 73-74) that the position of the radius of the vector in space is completely determined by its cosine guides.

From analytical geometry it is known that the cosine of a vector is equal to the ratio of its projection onto a given axis to the length of the vector. The length of the total velocity vector ( _for V) is found as the square root of the sum of its projections:

where V _x , V _y are the longitudinal wave velocities for a given seismic station along the x, y coordinates.

For each seismic station, the direction is measured.

The coordinates of the anomaly are found as the point of intersection of the radius of the velocity vectors from several seismic stations, as illustrated in Fig. 1.

After detecting the anomaly and its direction finding on the measured paths, the parameters of the expected seismic impact are forecasted. Why track the change in signal energy.

In accordance with the Rayleigh law, the energy of the signal spectrum is calculated as:

The change in signal attenuation on the measured path over time is determined as the difference | ΔE | = E (t ₂ ) - E (t ₁ ) where t ₁ , t ₂ is the time of the measurement samples for which the Fourier spectrum was calculated.

A well-known (see "Short-term forecast of catastrophic earthquakes using radiophysical ground-space methods. Conference reports, OIFZ named after O. Schmidt. M .: RAS, 1998, p. 9) an approximate dependence of the size of the zone of the prepared earthquake R [km] from the magnitude (M) of the expected seismic impact.

R [km] ≈ exp (M)

It is also known that the physical quantity itself and its rate of change are connected by a differential equation of the first degree, the general solution of which is exponential dependence (see N. S. Piskunov. "Differential and integral calculus for VTU3-s." M.: Nauka, 1976 ., p. 458). The exponential dependence has the property that, from its discrete values, the whole function can be restored.

The exponential dependence of the attenuation of the signal energy on the measured paths with increasing size (R) of the anomalous zone is illustrated in the graph of figure 4

The magnitude of the signal attenuation is judged on the size of the anomalous zone. The calculation of the size of the anomalous zone [R, km] to assess the expected magnitude (M) of seismic shock is given below in the example implementation of the method.

An example implementation of the method

The implementation of the method was carried out at the Kamchatka Geophysical Test Site RAS. The layout of the landfill is illustrated in figure 1. The range includes 6 seismic stations: Russkaya (RUS), Petropavlovsk (PET), Shipunsky (SPN), Cape Kozlova (MKZ), Krutoberegovo (KVT), Avacha (AVH).

Seismic monitoring stations record the arrival of longitudinal and transverse waves from each earthquake that occurred in the focal zone located 150 km from the coast of the Kamchatka Peninsula.

The focal zone of Kamchatka is a junction of continental and oceanic plates, the last of which carries out movement under the continental plate. To the north, in the monitoring zone of the KVT stations, there is a junction of these plates with a more northern one, i.e. triple articulation occurs, the seismic activity of which is very high. The focal zone of earthquakes is from the north of Kamchatka to the south, right up to Indonesia. In the monitoring of the environment by the stations RUS, PET, KVT, MKZ, SPN, AVH, data obtained only from local earthquakes were used.

Functional diagram of a measurement system that implements the method is illustrated in Fig.5.

The measurement system (figure 5) contains six seismic stations located in the Kamchatka region (1 ... 6). Each of the seismic stations, for the azimuths of the measurement paths in the focal zone, contains a group of two three-component seismic receivers, (7, 8) the sensitivity axes of the same name components of which are mutually parallel. The seismic receivers are located on the measuring base (9), which is a multiple of the average wavelength of the measured spectrum. The bases are perpendicular to the measured paths. The seismic receivers are installed buried coaxially with each other.

и сигнал пеленга аномалии (V_px/g_лV_p). Результирующие сигналы от канального коммутатора 13 оцифровывают в АЦП 14, 15 записывают в буферное ЗУ 16 для последующего ввода данных в компьютер 17, в стандартном наборе элементов: процессор 18, оперативное ЗУ 19, винчестер 20, дисплей 21, принтер 22, клавиатура 23. Компьютеры 17 объединены в локальную вычислительную сеть и обеспечена синхронизация их работы во времени.Through the built-in amplifiers 10, 11, 12, the signals from the output of the geophones are fed to the corresponding inputs of the channel switch 13. The channel switch implements the function of the relations of the sums of signals

and a bearing signal of the anomaly (V _px / g _l V _p ). The resulting signals from the channel switch 13 are digitized in the ADC 14, 15 are recorded in the buffer memory 16 for subsequent data input into the computer 17, in a standard set of elements: processor 18, operational memory 19, hard drive 20, display 21, printer 22, keyboard 23. Computers 17 are integrated into a local area network and their work is synchronized over time.

All elements of the device are existing technical developments. The device used elements of vibration measuring equipment from Bruel & Kjairi ENDEVCO, Denmark, the following models: built-in amplifiers 10, 11, 12 models 2626, 2628; channel switch, ADC, buffer memory (elements 13, 14, 15, 16) - multifunction block 3560-L.

- скорость колебательного движения, u - смещение частиц от положения равновесия в процессе колебательного движения. Поскольку давление во фронте сейсмической волны пропорционально скорости, для ее измерения на сейсмостанциях используют датчики давления. Известен класс устройств мембранного типа, позволяющих измерять инфразвуковые колебания земной коры (см., например, «Кондуктометрический датчик колебаний». Патент RU №2055352, G.01.N, 27/02, 1996 г.). Датчик имеет высокую чувствительность и линейность выходной характеристики в области инфрачастот 0,1…100 Гц, быстродействие 100 ms. Промышленная разработка кондуктометрического датчика представляет собой трехкомпонентный сейсмоприемник (см., например, «Способ определения гипоцентра нефтегазового месторождения». Патент RU №2150719, G01V 1/100, 2000 г.). Он может использоваться в любой пространственной ориентации, т.е. осуществлять измерения во взаимно ортогональных плоскостях.A seismic wave is characterized by an excess pressure P [n / m ² ] = ρc | ν |, where ρ is the density of the medium, c is the speed of the acoustic waves,

is the speed of oscillatory motion, u is the displacement of particles from the equilibrium position in the process of oscillatory motion. Since the pressure at the front of the seismic wave is proportional to the speed, pressure sensors are used to measure it at seismic stations. A known class of membrane-type devices that can measure infrasonic vibrations of the earth's crust (see, for example, "Conductometric vibration sensor". Patent RU No. 2055352, G.01.N, 27/02, 1996). The sensor has high sensitivity and linearity of the output characteristic in the field of infra-frequencies 0.1 ... 100 Hz, speed 100 ms. The industrial development of the conductivity sensor is a three-component seismic receiver (see, for example, “Method for determining the hypocenter of an oil and gas field.” Patent RU No. 2150719, G01V 1/100, 2000). It can be used in any spatial orientation, i.e. take measurements in mutually orthogonal planes.

The analysis of seismic data was carried out for the period from January 1, 2001 to December 31, 2004. The time intervals during which the changes in the state of the environment were considered were ~ 16, 8, 4, 2 days before the earthquake and, respectively, 2, 4, 8, 16 days after. Signal processing at each station was carried out using registered seismic waves taking into account the direction of their arrival, i.e. different azimuths were chosen (the azimuth can have a value from 0 to 360) and the decision depended on where the earthquake wave came from. For the KVT station, two azimuth intervals were chosen: the first from 0 to 130, the second from 50 to 110; an independent study was conducted for each azimuth interval.

For MKZ stations - from 70 to 130 and from 160 to 220. SPN - 70-130 and 110-190. PET-90-150, 100-130 and 180-230. RUS-60-110, 90-200 and 180-230.

The implementation of functions signals (V _p / V _s) 2 ^x stations (kW, MKZ) illustrated in Figure 2 graph. Samples of discrete samples of signals for stations were carried out in time intervals of 16, 8, 4 days before and after the earthquakes. Fourier spectra of these samples are illustrated in figure 3. When the earth's crust is saturated with a gas component, the medium becomes dispersive. The occurrence of such an anomaly on the measured paths leads to a change in the amplitude relations between the harmonics of the spectrum, as illustrated in Fig.3. The weighted average frequency of the calculated spectra, respectively, was F _cp1 (16 days) ≂ 5 Hz, F _cp2 (8 days) ≂ Hz, F _cp3 (4 days) ≂ 8 Hz. The signal energy of the measurement samples calculated by the Rayleigh ratio was respectively equal to: E ₁ (16 days) ≂0.306 mW, E ₂ (8 days) ≂ 0.0407 mW, E ₃ (4 days) ≂ 0.0203 mW. Since the measured signal represents the ratio of two functions (V _p / V _s ), the measured value of the signal energy (E) does not depend on the absolute values of the values of V _p , V _s , i.e. is a robust parameter. The signal attenuation time constant on the measured path, for the exponential process:

Estimated steady-state value of the signal energy on the eve of a seismic shock

An earthquake precursor sign is an increase in the size of the zone of the prepared earthquake (anomalous zone). As a first approximation, we can assume that the signal energy on the propagation path decreases inversely with E _min ~ 1 / R _{max to the} radius R of the anomalous zone, as illustrated in Fig. 1. Where does the prepared earthquake zone come from?

The expected magnitude of the seismic impact is ln [R, km] = M, ln [147 km] ≈ 5.1 points.

A summary graph of the calculated values of earthquakes for all six stations of the network in the interval of observation times 2001 - 2003 (in the observed azimuths of the stations) is illustrated in Fig.6. The convergence of the predicted results (with an accuracy of up to the ^2nd sign) corresponds to the posterior data of the earthquakes. The effectiveness of the method is characterized by high reliability, the scale of the territories covered by monitoring and a large interval of time proactive forecast.

Claims (1)

A method for monitoring a region by a network of seismic stations, including route measurements of the ratio of the propagation velocities of longitudinal V _p and transverse V _s waves in the earth's crust (V _p / V _s ) under different azimuths of the tracks of each network station using pressure sensors placed as primary velocity sensors in mutually orthogonal measurement planes V _p and V _s , obtaining discrete measurement samples (V _p / V _s ) for each trace at the same time intervals [V _p / V _s (t _i )], calculation of the Fourier spectra of the obtained samples with the definition of medium Addressing frequency samples, the identification of the change of average frequency with the occurrence of the anomalous zone being prepared earthquake measured track direction finding abnormal area as the intersection of radius vectors point velocities by multiple stations on the measured slopes, the determination signal energy spectra calculated as the

and tracking the dynamics of its change in the form of an exponential E (t) = E ₀ * exp (-t / T), calculating the size of the anomalous zone R, km ≂ 1 / E _mouth , forecasting the magnitude of the expected seismic impact M ≂ ln [R, km] and the time of the impact, measured from the moment the weighted average frequency of the Fourier spectrum changes: t _y ~ 4.7T, where

- the time constant of the signal attenuation on the measured path;
Δt is the time interval between the calculated definitions of the energy of the signal E ₁ and
E ₂ ;
E ₀ is the signal energy until the anomaly occurs;
E _set ≅E ₀ * е ^-4.7 is the steady-state value of the energy of the signal passing through the anomaly on the eve of the seismic shock.

Images