RU2254599C1 - Method of predicting parameters of earthquake - Google Patents

Abstract

FIELD: GEOPHYSICS.

SUBSTANCE: method can be used for monitoring of environment. According to the method images of seismic center can be received from two polarization mutually orthogonal receiving channels. Synthesized image array if formed from pixel-to-pixel relation of signal amplitudes in polarization channels. Area of relief of synthesized array is calculated in relation to geometrical area of photograph and changes in the relation are traced by means of set of photographs made. Parameters of earthquake are predicted from numerical characteristics of relation to be traced.

EFFECT: improved reliability of prediction.

5 dwg

Description

The invention relates to seismology and may find application in national space observation systems for remote monitoring of natural environments.

In recent years, several new, previously inaccessible for observation, signs have been revealed - precursors of future earthquakes based on the use of remote sensing methods of the Earth's surface from space. Among them are such signs as: the buildup of the earthquake source on the eve of the impact, monitored by the Navstar space system (Patent RU No. 2.170.446, 2001 - analogue), the formation of a vertical electrostatic field above the epicentral region of the source, with a voltage of several kV / m ( Patent RU No. 2.205.432, 2003 - analogue). The well-known "Method for detecting earthquake foci", Patent RU No. 2.209.452, 2003 - analogue.

In the analogue method, an image of the underlying surface in the violet-blue portion of the visible range is obtained in the form of a digital matrix of the luminance signal I (x, y) from spatial coordinates, the contours in the image are isolated by spatial differentiation methods, the functions of the fractal dimension of the image signal inside the selected contours and the gradient field of directions of lineaments, form a pattern of the directions of lineaments inside the contours and along the geometry of the contour, the pattern of the azimuths of the directions of lineaments, fra The total dimension of the image inside the contour is judged on the belonging of the identified anomaly to the space of the earthquake source.

A disadvantage of the known solution is the impossibility of accurately predicting the parameters of the impending impact on the operations of the analogue method.

The closest analogue to the claimed technical solution is the "Method for predicting earthquakes", Patent RU No. 2.213.359, 2003, G.01.V, 9/00, 2003

методами пространственного дифференцирования выделяют контуры на синтезированном изображении, проводят фрактальный и линиаментный анализ сигналов фрагментов изображения внутри выделенных контуров и формируют градиентное поле направлений линиаментов, идентифицируют очаг землетрясения на изображении по геометрии контура, фрактальной размерности и узору рисунка поля линиаментов, вычисляют средневзвешенную сумму азимутов линиаментов и отслеживают изменение этой функции по серии последовательных изображений идентифицированного очага, фиксируют момент появления периодических отклонений функции средневзвешенной суммы азимутов линиаментов, рассчитывают период (Т) и коэффициент (k) наклона касательных отслеживаемой функции, прогнозируют Т время удара

отсчитывамое от начала периодических отклонений функции и магнитуду удара из соотношения lgt_y[cут]=0,54M-3,37, где u₀ - предельная величина сдвиговой деформации, при которой происходит разрыв земной коры.In the closest analogue method, images of the underlying surface A ₁ (x, y) are obtained; And ₂ (x, y) along two reception channels mutually orthogonal in polarization, form a synthesized matrix by pixel-by-pixel combination of the resulting images with the signal amplitude at each point

using spatial differentiation methods, isolate the contours in the synthesized image, conduct fractal and lineament analysis of the signal fragments of the image inside the selected paths and form a gradient field of lineament directions, identify the earthquake in the image by the geometry of the contour, fractal dimension and pattern of the line field pattern, calculate the weighted average sum of line azimuths track the change in this function by a series of consecutive images identified about Aha, fixed time of occurrence of the periodic deviations of the average weighted sum function liniamentov azimuth is calculated period (T) and the coefficient (k) of inclination of tangent tracking functions predict impact time T

reckoned from the beginning of periodic deviations of the function and impact magnitude from the relation lgt _y [s] = 0.54M-3.37, where u ₀ is the ultimate value of the shear strain at which the crust breaks.

The disadvantages of the method of the closest analogue are:

- low sensitivity associated with a limited (unit grad.) range of change of the traced sign - a precursor;

- mathematical inaccuracy, since the envelope of the amplitude of the buildup increases not exponentially, but exponentially.

The problem solved by the claimed method is to increase the sensitivity of the method, statistical stability and reliability of forecasting by calculating the surface area of the surface of the synthesized matrix and tracking the dynamics of changes in this area over a series of obtained images of the focus.

вычисляют постоянную экспоненты Т наблюдаемого процесса, рассчитывают ожидаемое время удара: (t_y)≅4,6T и его магнитуду (М) из соотношения lgt_y[сут]=0,54М-3,37, где:The problem is solved in that in a method for predicting earthquake parameters, including obtaining focal images from two survey channels mutually orthogonal in polarization, forming a synthesized image matrix from pixel-by-pixel amplitudes of polarization channel signals, processing the synthesized image matrix, and additionally, shooting by means of the orbital complex is carried out in mode with equal time intervals Δt between adjacent turns, calculate the mathematical expectation of image signals polarization channels, the synthesized matrix is formed from pixel-by-pixel ratios of signal amplitudes with a higher mathematical expectation to a smaller one, sequentially, from the beginning of the synthesized matrix, the image is divided into fragments of four adjacent pixels, each quadrangle is divided into two adjacent triangles with a diagonal, and their areas are calculated using the Heron formula finding the area S _p of the relief surface of the image as the sum of the mosaic areas of triangles all the fragments of the time-series images hearth tracking dynamics functions - D = S _p / S _0, the roughness of the relief, depending on

calculate the constant of the exponent T of the observed process, calculate the expected impact time: (t _y ) ≅4.6T and its magnitude (M) from the relation logt _y [days] = 0.54M-3.37, where:

t ₂ -t ₁ = Δt, time interval between two consecutive surveys, [hour];

D ₁ , D ₂ are the calculated functions of the roughness of the relief at moments t ₁ t ₂ ;

D ₀ - the limit of the roughness of the relief, at which there is a "rupture" of the focus;

S ₀ is the area of the synthesized image matrix from | m × n | pixels.

The invention is illustrated by drawings, where:

figure 1 - a change in the albedo (A) of the surface and the polarizability (φ) of the light flux from the intensity of the electrostatic field in the atmosphere above the source;

figure 2 - visualization of the synthesized image matrix (printout with a PC);

figure 3 - function changes in the roughness of the relief D = S _p / S ₀ at various T;

4 is a sequence of calculating the surface area of the surface of the synthesized matrix as the sum of the areas of the mosaic of triangles;

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

The technical essence of the proposed method is as follows. On the eve of the impact, due to the asymmetry of the application of compression forces and elastic recoil relative to the fault, a moment arises leading to the buildup of the earthquake source. The characteristics of the transient oscillatory process: the period of oscillations and the rate of increase of the buildup amplitude (see, for example, “Method for predicting earthquakes”, Patent RU No. 2.204.852, 2003) contain information about the parameters of the impending impact. The larger the vibrational mass of the source and the stored potential energy, the longer the oscillation period and magnitude of the upcoming earthquake. The buildup of the focus is accompanied by the release of ionized gases from the earth’s crust, the formation of an uncompensated electric charge in the atmosphere and the appearance of a vertical electrostatic field of several kV / m over the epicentral region [see, for example, “Short-term forecast of catastrophic earthquakes using radiophysical ground-space methods”, Reports conferences, OFIF them. O. Yu. Schmidt, RAS, M., 1998, p. 27, Fig. 1].

The presence of a vertical electrostatic field above the source surface causes the polarization of water vapor molecules contained in the atmosphere. Based on the classical principle of the interaction of an electromagnetic field with matter [see, for example, Refraction of Light. Physical Encyclopedic Dictionary, edited by A.M. Prokhorova, Izov. Encyclopedia, Moscow, 1983, p. 168] electrons and atoms of a substance undergo a forced oscillation under the influence of a light wave. The refractive index of a substance (atmosphere) for an incident light flux depends both on the concentration of the secondary emitters in the substance and on the ratio of the wavelengths of the incident light flux and the intrinsic radiation of the secondary vibrators.

The greater the magnitude of the electrostatic field above the surface of the source, the higher the density of charges induced by the polarization of water molecules, the greater the dielectric constant ε and the refractive index of the medium.

As a result of the secondary reradiation of the supplying light flux by dipole-oriented water molecules, a decrease in the surface albedo above the earthquake source is observed. Moreover, if the wavelength of the incident light flux approaches the length of the water molecules, which takes place in the blue-violet part of the visible range, then the refractive index increases further.

Resonant re-emission of light by dipole-orienting water molecules leads to polarization of the waves of the reflected light flux. This means that the albedo of the surface of the earthquake source in magnitude and direction in space is proportional to the magnitude of the electric field and the direction of the electric lines of force of this field above the source. The change in the albedo of the source surface (A%) and the polarization of the light flux (φ °) from the magnitude of the vertical electrostatic field in the atmosphere is illustrated in FIG. 1. When receiving the reflected light flux through two channels mutually orthogonal in polarization, the difference in signal polarization is converted into the difference in pixel amplitudes of identical image sections in these channels. To emphasize the contrast, pixel-by-pixel combination of the obtained images is carried out in the reception channels orthogonal in polarization. The mathematical expectation of the amplitudes of the image signal in each channel is preliminarily calculated. The synthesized image matrix is formed from pixel-by-pixel ratios of signal amplitudes with a higher mathematical expectation to a smaller one.

Natural sunlight is not polarized. Therefore, the ratio of the amplitudes of the pixels of the images of the unpolarized light flux is close to unity, and the surface topography of the resulting image matrix is almost smooth. Figure 2 shows that on the periphery of the focus the image fury is small, and the tone is almost dark. Sufficient roughness (roughness) of the synthesized image, on the contrary, indicates the polarizability of the light flux reflected from the surface of the focus. The tone of the image in figure 2 in the center of the focus is lighter. The degree of indentation of the amplitude of the pixels of the synthesized image matrix of the source contains information about the parameters of the expected impact. The higher the intensity of the vertical electric field above the source, the higher the polarizability of the reflected light flux and the greater the roughness of the relief (or turbulence of the resulting signal).

It is proposed to use the surface area of the surface of the synthesized image matrix as a statistically stable sign - a harbinger of an impending earthquake.

Since the relief area in each case depends on the area of the observed focus (image matrix sizes of | m × n | pixels), a relative value is introduced that characterizes only the "roughness" of the relief: D = S _p / S ₀ , where S _p is the calculated area topography, S ₀ is the area of the rectangle equal to the product of the number of rows (m) by the number of columns (n) and the area of one image pixel. The interval of change of function D is wide enough, which determines the high sensitivity of the method to the entered indicator.

The buildup of the earthquake source (and with it the intense release of ionized gases into the atmosphere) continues until the dynamic pressure of the vibrational mass exceeds the tensile strength of the earth's crust, at which its rupture occurs, cleavage, accompanied by an impact. Since the "tearing up" of the earthquake source occurs at the same values of the ultimate strength of the earth's crust, there is a limit (D ₀ ) to which the roughness function of the image relief associated with this process tends. The dynamics of the change in the function of the relief area contains information on the magnitude and time of the impact. The relationship between process dynamics (speed) and the function itself is represented by a differential equation of the first degree. It is known [see, for example, N. S. Piskunov, “Differential and Integral Calculus” for Technical Schools, textbook, 5th ed., M., Nauka, 1964, pp. 443-451], that a general solution of the differential equation The first degree is the exhibitor. The initial conditions for the numerical solution of the differential equation are determined by two consecutive in time, images of the tracked focus. From the properties of the exhibitor follows, see figure 3):

where t ₂ -t ₁ = Δt, time interval

and between two consecutive shootings, an hour;

D ₁ D ₂ - calculated functions of the roughness of the relief at moments t ₁ t ₂ ;

T is the exponent constant (initial conditions).

The value of D _{0 is} calculated a posteriori by a series of earthquakes.

The set of solutions of the original differential equation for various values of the constant exponent is illustrated in Fig.3. The exponent constant characterizes the transitional oscillatory process from the moment the start of the buildup of the source to the impact. The larger the vibrational mass of the source, the greater the process time constant T. There is a Gutenberg-Richter formula relating the lifetime of the precursor sign (t) to the magnitude of the expected earthquake M [see, for example, “Short-term forecast of catastrophic earthquakes using radiophysical, ground-space methods”, Conference reports, UIFZ im. O. Yu. Schmidt, RAS, M., 1988, p. 10].

For seismic waves, this dependence has the form:

lgt [s] = 0.54M-3.37.

The lifetime of the precursor sign is determined from the properties of the exponent. With a probability of 0.99, the exponent reaches its limit value at t = 4.6T. Thus, the impact time, measured from the moment of violation of the "smoothness" of the relief, will be the value of t _y = 4,6T.

The constant T itself is found from the above dependence by calculating the values of the relief roughness D ₁ and D ₂ at the shooting times t ₁ , t ₂ .

The entire relief area is represented by a mosaic of triangles. The area of each triangle is calculated by the Heron formula:

- полупериметр сторон. В свою очередь длину каждой из сторон треугольника находят по теореме Пифагора при известных разрешениях на пиксел (Δ_x, Δ_у) и шаге квантования амплитуды сигнала (А) синтезированной матрицы:Where:

- half perimeter of the parties. In turn, the length of each side of the triangle is found by the Pythagorean theorem with known pixel resolutions (Δ _x , Δ _y ) and the quantization step of the signal amplitude (A) of the synthesized matrix:

The sequence of calculating the relief area is illustrated in figure 4. The image matrix, sequentially, from the beginning, is divided into fragments of four adjacent pixels. The main diagonal 1-4 (from top to bottom) divides the quadrangle into two adjacent triangles.

The area of each of the triangles is calculated, and then the result is summed over the entire mosaic of triangles. The calculation is carried out by software on a PC. The calculation program is shown in the example implementation of the method.

An example implementation of the method.

The inventive method can be implemented according to the scheme of figure 5. The functional diagram of the device of FIG. 5 contains an orbital station 1, of the ISS type, with digital coaxial digital cameras orthogonal to the polarization 2,3, type BYP-30 (Japan), mounted on its board, receiving reflected light flux from the underlying surface.

The images received in the reception channels are recorded on the Niva type 5 on-board tape recorder. The survey of the planned seismic regions is carried out according to the programs or one-time commands laid down in the onboard control complex 6 transmitted from the control center (MCC) 7.

The captured (by target designation from the MCC) image frames of the surface of the monitored foci are transmitted via the operational communication channel (8) to the ground-based relay station 9, where they are recorded on an Arktur video recorder 10. The compiled information files (along with service signs, time taken, coordinates) are transferred to the Geophysical Center of the Ministry of Emergencies 11, where a long-term archive of 12 all captured earthquake sources is created, based on streamers, such as FT-120. The processing of the obtained images of the foci, according to the operations of the proposed method, is carried out on a PC 13 in a standard set of elements: a processor 14, an operational memory 15, a hard drive 16, a display 17, a printer 18, a keyboard 19. According to the results of image processing, highlighting and tracking the dynamics of the sign of the precursor create a database of reference data 20, which is displayed on the server 21 of the Internet.

Obtained in the polarization channels of the survey 2, 3 images are injected into a PC 13, from which the synthesized matrix is formed. The calculation of the synthesized matrix by the well-known algorithm of the pixel-by-pixel ratio of amplitudes is represented by the standard mathematical operation included in the specialized software package MATH CAD 7.0 PLUS [see, for example, MATH CAD, the second stereotyped publication, Informizdat, house. Filin, M., 1997, Vectorization of matrix elements, p. 211]. The original synthesized matrix (figure 2) contained 512 rows and 480 columns obtained in the video channels with a resolution of Δx = Δy = 120 m per pixel. Moreover, the surface area of the outbreak on the ground is

S ₀ = [512 × 120 × 480 × 120] ≈356 km ² .

The amplitude of the pixels was quantized in a standard scale of 0 ... 256 levels.

Since the calculation uses the relative size of the relief surface S / S ₀ , the area of the triangle is calculated in units of the product of the pixel resolution by the image depth, i.e. a quantization step equal to 1 m. Calculation of the relief area is carried out by the software method, according to a specially developed program.

The text of the program for calculating the relief area.

Program SqrGeron;

Const maxx = 1000;

Var S: real;

Ss: string;

St1, st2: array [1 ... maxx] of byte;

F: text;

i, j kolx: integer;

dx, dy, dh: real;

procedure SqrSqr (var S: real; x1, x2, y2: real;

h11, h12, h21, h22: byte);

var dS, dS1, dS2: real;

function Geron (x1, y1, h1, x2, y2, h2, x3, y3, h3: real): real;

var 112, 123, 131, p: real;

begin

112: = sqrt (sqr (x1-x2) + sqr (y1-y2) + sqr (h1-h2));

123: = sqrt (sqr (x2-x3) + sqr (y2-y3) + sqr (h2-h3));

131: = sqrt (sqr (x3-x1) + sqr (y3-y1) + sqr (h3-h1));

p: = (112 + 123 + 131) / 2;

Geron: = sqrt (p * (p-112) * (p-123) * (p-131));

end;

begin

dS1: = 0;

if (h11 <> 255) and (h12 <> 255) and (h21 <> 255) then

dS1: = dS1 + Geron (x1, y1, h11 * dh, x2, y1, h12 * dh, x1, y2, h21 * dh);

if (h22 <> 255) and (h12 <> 255) and (h21 <> 255) then

dS1: = dS1 + Geron (x2, y2, h22 * dh, x2, y1, h12 * dh, x1, y2, h21 * dh);

dS2: = 0;

if (h11 <> 255) and (h12 <> 255) and (h22 <> 255) then

dS2: = dS2 + Geron (x1, y1, h11 * dh, x2, y1, h12 * dh, x2, y2, h22 * dh);

if (h22 <> 255) and (h11 <> 255) and (h21 <> 255) then

dS2: = dS2 + Geron (x2, y2, h22 * dh, x1, y1, h11 * dh, x1, y2, h21 * dh);

if dS1> 0 then dS: = dS1;

if (dS2> 0) and (dS2 <dS) then dS: = dS2;

S: = S + dS

end;

begin {main}

S: = 0;

repeat

writeln (Enter file name, please ');

readln (ss);

assign (f, ss);

{SI-} reset (f); {SI +}

i: = IOresult;

until i = 0;

writeln ('Enter dx, dy, dh');

readln (dx, dy, dh);

kolx: = l;

if not eof (f) then

while not eoln (f) or (kolx> maxx) do begin

read (f, st1 [kolx]); kolx: = kolx + 1;

end;

if eoln (f) then begin;

readln (f);

j: = 1;

while not eof (f) do begin

for i: = 1 to kolx do read (f, st2 [i]);

for i: = 1 to kolx-1 do

SqrSqr (S, 0, dx, 0, dy, st1 [i], st1 [i + 1], st2 [i], st2 [i + 1]);

readln (f);

St1: = st2;

j: = j + 1;

end;

end else writeln ('Very short array');

close (f);

assign (f, 'result.txt');

rewrite (f);

writeln (f, 'file name -', ss);

writeln (f, dx = ', dx: 10: 5,' dy: 10: 5, 'dh =', dh: 10: 5);

writeln (f, 'Strok -'j: 5,' Stolb - ', kolx);

writeln (f, 'Area of region is', S: 10: 2);

writeln ('file name -', ss);

writeln ('dx =', dx: 10: 5, 'dy =', dy: 10: 5, 'dh =', dh: 10: 5);

writeln ('Strok -', j: 5, 'Stolb -', kolx);

writeln ('Area of region is', S: 10: 2);

close (f)

end. {main}

в полученных снимках составило 2,2. Откуда постоянная экспоненты T=3/ln2,2=4,25 часов.The survey was carried out in the visibility range of the ISS from a ground station at the beginning of the zone (at the first turn) and at the end of the zone (at the third turn). The time interval between two consecutive surveys was: 1.5 hours per revolution × 2 = 3 hours. Attitude

in the resulting images was 2.2. From where does the exponent constant T = 3 / ln2.2 = 4.25 hours.

The expected impact time t _y = 4.25 × 4.6 = 19.5 = 0.82 days. The expected impact magnitude is lg0.82 = 0.54M-3.37. M = 6.1 points. A large range of changes in the precursor attribute in the form of the surface area of the synthesized matrix (about 12) and the representativeness of the measurement sample [(volume m × n) samples] provide high sensitivity and statistical stability of the forecasting method considered.

Claims (1)

A method for predicting earthquake parameters, including obtaining focal images from two survey channels mutually orthogonal in polarization, generating a synthesized image matrix from pixel-by-pixel amplitudes of the polarization channel signals, processing the synthesized image matrix, characterized in that the shooting by means of the orbital complex is carried out in a mode with equal time intervals Δt between adjacent turns, calculate the mathematical expectation of the image signals of the polarization channels, syn ezirovannuyu matrix form of pixel-wise relationships amplitude signal with great expectation to smaller sequentially from the beginning of the synthesized matrix divide images into fragments of four adjacent pixel, divide each rectangle a diagonal of two adjacent triangles and Heron's formula calculating their area, find relief area (S _p ) image surfaces as the sum of the areas of the mosaic of triangles of all fragments, according to a series of successive time-domain images of the focus, the dynamics of change is monitored function D = S _p / S ₀ surface roughness, depending

calculate the exponent constant T of the observed process, calculate the expected impact time t _at ≅4.6T and its magnitude (M) from the relation logt _y [days] = 0.54M-3.37, where t ₂ -t ₁ = Δt is the time interval between two consecutive shoots, h; D _1, D ₂ are the calculated functions of the roughness of the relief at moments t ₁ , t ₂ ; D ₀ is the roughness limit of the relief area at which the “vaporization” of the focus occurs; S ₀ is the area of the synthesized image matrix from | m × n | pixels.

Images