A Stratum Structure Inversion Algorithm Based on Microseismic Detection

Abstract

The invention relates to a stratum structure inversion algorithm based on micro-motion detection, which solves the technical problems of slow calculation speed, low inversion result accuracy and insufficient ability to describe the details of stratum characteristics. The invention adopts step 1 of data acquisition; step 2 of processing the acquired original micro-motion signal data, including noise removal, filtering and normalization operations; step 3 of extracting features that can effectively characterize the relevant information of stratum structure from the pre-processed data; step 4 of constructing electromagnetic characteristic parameter vectors according to the physical properties of each layer of the stratum such as conductivity, magnetic permeability and dielectric constant, clarifying the relationship between the electromagnetic characteristic parameter vector and the stratum structure and the value law, and establishing a forward model; step 5 of performing inversion calculation by executing at least two distributed parallel inversion calculation models; step 6 of performing correlation calculation on the inversion results, and determining the optimal inversion result as the output technical solution, which better solves the problem and can be used in micro-motion detection.

Description

Technical Field

[0001] The present invention relates to the field of micro-motion detection, and in particular to a stratum structure inversion algorithm based on micro-motion detection. Background Art

[0002] In many fields such as geological exploration, engineering construction and underground resource development, it is crucial to accurately obtain information on the stratigraphic structure. Although traditional stratigraphic structure detection methods, such as drilling and seismic exploration, can provide relevant information to a certain extent, they each have obvious limitations.

[0003] The drilling method can directly obtain underground core samples, so as to accurately know the rock characteristics of the formation at the drilled location. However, drilling is expensive and inefficient, and can only obtain local information about the drill hole location, making it difficult to fully reflect the formation structure over a large area. The detection of deep formations faces many technical difficulties and high cost challenges.

[0004] Seismic exploration is one of the most widely used geophysical exploration methods. It artificially excites seismic waves and uses the differences in the propagation characteristics of seismic waves in different stratigraphic media to infer the stratigraphic structure. However, seismic exploration requires professional source equipment for seismic wave excitation. The equipment is complex and requires high operation requirements. At the same time, in complex geological environments, such as urban areas or places with a large number of artificial interference bodies underground, seismic wave signals are easily interfered, which affects the accuracy of the interpretation results.

[0005] With the continuous development of science and technology, micromotion detection technology has gradually attracted attention. Micromotion refers to weak ground vibrations caused by nature or human activities, which contain rich stratigraphic information. However, existing micromotion detection technologies often lack efficient and accurate algorithm support when analyzing and processing micromotion signals to invert stratigraphic structures. Conventional algorithms may have problems such as slow calculation speed, low inversion accuracy, and insufficient ability to depict details of stratigraphic characteristics when dealing with complex stratigraphic structures and large amounts of micromotion data. They cannot meet the urgent needs of actual engineering and scientific research for high-precision and rapid acquisition of stratigraphic structure information.

[0006] In summary, there is an urgent need for an innovative stratigraphic structure inversion algorithm based on micro-seismic detection, which can overcome the shortcomings of existing detection methods and related algorithms, obtain stratigraphic structure information more efficiently and accurately, and provide strong support for the development of related fields. Summary of the invention

[0007] The technical problem to be solved by the present invention is that the calculation speed is slow, the inversion result precision is not high, and the ability to describe the details of the formation characteristics is insufficient in the prior art; a new formation structure inversion algorithm for micro-seismic detection is provided, which has the characteristics of fast calculation speed, high inversion result precision, and strong ability to describe the details of the formation characteristics;

[0008] In order to solve the above technical problems, the technical solutions adopted are as follows:

[0009] A stratum structure inversion algorithm based on micro-seismic detection, characterized in that the algorithm includes:

[0010] Step 1, data acquisition, includes rationally arranging a sensor array containing multiple micro-motion sensors in the target detection area, which can fully cover and obtain micro-motion signals reflecting the formation characteristics;

[0011] Step 2, data preprocessing: processing the collected raw micro-motion signal data, including noise removal, filtering, and normalization operations;

[0012] Step 3, feature extraction: Extract features that can effectively characterize information related to the stratigraphic structure from the preprocessed data, including mean features Variance characteristics Spectral characteristics; Spectral characteristics include spectrum peak frequency f peak，i , spectrum bandwidth B i ;

[0013] Step 4, establishing a forward model;

[0014] Step 5, performing inversion calculation: performing inversion calculation by executing at least two distributed parallel inversion calculation models;

[0015] Step 6, performing correlation calculation on the inversion results obtained in step 5, and determining the most accurate and optimal inversion result as output;

[0016] Step 4 includes:

[0017] Step 4.1 Define the relationship between electromagnetic characteristic parameters and stratum structure; let the electromagnetic characteristic parameter vector related to stratum structure be θ em =[θ em1 ,θ em2 ,…θ emD ] T ;

[0018] Where D is the dimension of the electromagnetic characteristic parameter;

[0019] Step 4.2: transform and map the mean, variance, spectrum features, etc. extracted in step 2;

[0020] Suppose the micro-motion detection feature input vector after transformation and mapping is Where K is the transformed feature dimension, and its elements are obtained by processing the features obtained in the preprocessing and feature extraction steps;

[0021] Among them, the mean feature conversion is the mean feature extracted for the i-th sensor Converted into parameters related to electromagnetic forward modeling Among them, a1 and b1 are predetermined coefficients;

[0022] Variance feature transformation is the variance feature extracted for the i-th sensor Convert to parameter Among them, a2 and b2 are predetermined coefficients used to adjust the variance characteristics to suit electromagnetic forward modeling;

[0023] The spectrum feature conversion is for the spectrum peak frequency f peak，i , converted to parameter F em3 =a3f peak，i +b3; for spectrum bandwidth B i , converted to parameter F em4 =a4B i +b4; where a3, b3, a4, and b4 are also predetermined coefficients used to incorporate the spectrum characteristics into the input parameters of the electromagnetic forward model;

[0024] Step 4.3, with F em As the input feature vector, θ em As electromagnetic property parameter vector, fast forward modeling algorithm is executed, including:

[0025] (1) Spatial discretization, defining the problem space: setting the three-dimensional space range, including geometric dimensions such as length, width, and depth, denoted as V strata ;

[0026] Determine the distribution of electromagnetic physical quantities including electric field intensity E and magnetic field intensity H that need to be solved in the three-dimensional space of the formation;

[0027] Generate unstructured grid: According to the geometric shape of the stratum structure and the complexity of the electromagnetic characteristics, use a suitable unstructured grid generation algorithm to generate the spatial region V strata Discretize; divide the continuous three-dimensional space of the stratum into irregular unit structures, and the discretized grid is Ω strata , the node set in the grid is recorded as N strata , the unit set is recorded as E strata ; For each node i∈N strata , record its spatial coordinates (x i ,y i , z i );

[0028] Establish the finite element linear equations: Based on Maxwell's equations, on the unstructured grid Ω strataThe vector finite element method is applied to discretize the above. By using appropriate vector basis functions to expand the electromagnetic field and substituting them into the variational form of Maxwell's equations, the finite element linear equations for the unknown electromagnetic coefficients are obtained:

[0029] Ax = b;

[0030] Among them, A is a large sparse complex coefficient matrix, whose elements are related to the electromagnetic characteristic parameters of the grid unit and the integral calculation of the basis function, and the correlation is predefined; x is an unknown vector containing the electromagnetic coefficients to be solved at each node; b is a known right-hand side vector, which is related to the source term and boundary conditions;

[0031] (2) Single-pole model reduction; select the optimal pole, in the frequency range [ω min ,ω max ], analyze the characteristics of the coefficient matrix A of the equation group with frequency, and select a single optimal pole p through a predefined optimization criterion. opt ; Using the selected optimal pole p opt , transform the original large sparse complex coefficient equation system Ax = b; let the transformed equation system be: A real x real =b real ;

[0032] Among them, A real is the transformed real coefficient matrix, which has the same dimension as the original matrix A, but the coefficients are real numbers; x real is the corresponding real unknown vector; b real is the transformed real right-hand side vector;

[0033] (3) Lanczos algorithm constructs an orthogonal basis, including initializing the Lanczos algorithm: using the transformed real coefficient matrix A real As input, set the initial vector υ1 of the Lanczos algorithm; usually a random vector that satisfies the normalization condition ||υ1||=1 can be selected;

[0034] Gradually iterate to construct an orthogonal basis: for k = 1, 2, ..., m, where m is a pre-set subspace dimension, much smaller than the matrix A real The order of ; perform the following iterative steps:

[0035] Calculate the vector w k =A real υ k ;

[0036] For vector w k Perform orthogonal processing to obtain the orthogonal vector υ k+1 ;

[0037] The specific orthogonalization process is as follows: β k =||w k ||,

[0038] Among them, α k and β k are the coefficients generated during the orthogonalization process;

[0039] After m iterations, a set of orthogonal vectors {υ1,υ2,…,υ m+1}, which form an orthogonal basis of the rational Krylov subspace;

[0040] For each orthogonal vector v i (i=1,2,…,m+1), calculate its norm ||υ i ||, then the vector υ i Divide by its norm to get the normalized orthogonal vector

[0041] (4) Calculate the orthogonal projection, including constructing the projection matrix: Use the normalized orthogonal basis vectors obtained Construct a projection matrix P; the column vectors of the projection matrix P are composed of orthogonal basis vectors, that is;

[0042] Compute the orthogonal projection of the initial matrix:

[0043] Calculate the initial matrix A real Orthogonal projection A onto Krylov subspace proj ; According to the definition of projection matrix and the calculation formula of orthogonal projection, we have:

[0044] A proj =P T A real P;

[0045] Similarly, for the right-hand vector b real Also project and get b proj =PTb real ;

[0046] (5) Fast forward solution, including the use of projection matrix solution:

[0047] Using the projection matrix A proj To replace the initial matrix A real Solve it; get the projected system of equations:

[0048] A proj x proj =b proj ;

[0049] Since the projection matrix A proj The dimension is much smaller than the initial matrix A real So the speed of solving this system of equations will be greatly improved;

[0050] Solve the projected system of equations: Use the linear equations solving algorithm to solve the projected system of equations A proj x proj =b proj , and get the solution vector x proj ;

[0051] Restore the original solution: Through the inverse projection operation, the solution vector x proj Restore the solution vector x to the original space; specifically, according to the properties of the projection matrix P, we have:

[0052] x=Px proj The solution to the electromagnetic characteristics of the stratigraphic structure related to microtremor detection will serve as the basis for subsequent inversion calculations.

[0053] Working principle of the present invention: The present invention realizes formation structure inversion based on micro-motion detection, and mainly includes the following key steps and principles: First, data acquisition is carried out, and signals reflecting formation characteristics are obtained by arranging micro-motion sensors. Then, data preprocessing is performed to remove noise, filter and normalize the data to improve data quality. The focus is on step 4 to establish a forward model: on the one hand, an electromagnetic characteristic parameter vector is constructed based on the physical properties of each layer of the formation, such as conductivity, magnetic permeability, and dielectric constant, to clarify its relationship with the formation structure and the law of value selection. On the other hand, the characteristics of the micro-motion detection data after preprocessing and feature extraction are mapped to the input parameters of the forward algorithm through a specific transformation. Subsequently, the forward algorithm is executed, and through operations such as spatial discretization, single-pole model reduction, Lanczos algorithm to construct orthogonal basis and calculate orthogonal projection, etc., a rapid solution of the electromagnetic characteristics of the formation structure related to micro-motion detection is achieved, providing a basis for subsequent inversion calculations, thereby accurately inverting the formation structure.

[0054] In the above scheme, further, the data collection in step 1 specifically includes:

[0055] Arrange multiple micro-motion sensors in the target detection area to form a sensor array to ensure that it can fully cover and obtain micro-motion signals that sufficiently reflect the characteristics of the formation; the number of arranged micro-motion sensors is N, and the data collected by each sensor at different times constitutes a time series; the data collected by the i-th sensor at the j-th time is x ij , where i = 1, 2, ..., N, j = 1, 2, ..., M, where M represents the length of the time series of data collected by each sensor;

[0056] Furthermore, the data preprocessing includes noise removal, filtering, and normalization operations;

[0057] The noise removal includes setting the original micro-motion signal data as x ij , the concave-convex coefficient obtained after the concave-convex transformation is W ij (s, k), where s represents the scale parameter and k represents the translation parameter;

[0058] The threshold λ can usually be determined based on experience or some statistical method, such as using the general threshold formula: Where σ is an estimate of the standard deviation of the noise, NM is the total number of data points, and the number of sensors N is multiplied by the length of each sensor data sequence M;

[0059] Threshold the bump coefficient:

[0060]

[0061] Then the processed concave and convex coefficients are transformed by inverse concave and convex transformation. Reconstruct the denoised signal

[0062] First, the convex-concave transformation is used in the noise removal process. By setting a reasonable threshold to process the convex-concave coefficient, various types of noise in the original micro-motion signal can be effectively filtered out, the essential characteristics of the signal can be restored, and the data purity can be improved. The filtering operation uses the Butterworth filter to remove specific frequency interference and highlight the effective micro-motion frequency band. Finally, normalization processing is performed to map the data to a specific interval, so that the data collected by different sensors are comparable, laying a good foundation for the subsequent accurate feature extraction and inversion calculation, and greatly improving the accuracy and reliability of the inversion results as a whole.

[0063] Furthermore, the mean feature in step 3 is: for the normalized data sequence y after preprocessing of the i-th sensor i1 ,y i2 , …, y iM , its mean characteristic The calculation formula is: in Represents the mean feature of the i-th sensor data;

[0064] Variance feature: Normalized data y after preprocessing for the i-th sensor i1 ,y i2 , …, y iM sequence, whose variance characteristic The calculation formula is: in Represents the variance feature of the i-th sensor data;

[0065] Spectrum characteristics:

[0066] The normalized data sequence y after preprocessing of the i-th sensor i1 ,y i2 , …, y iM Perform discrete Fourier transform DFT to obtain its spectrum representation Y i (k), where k = 0, 1, ..., M-1 represents the frequency point, and the discrete Fourier transform formula is:

[0067] Then we can get the spectrum Y i (k) Extract the spectrum peak frequency and spectrum bandwidth.

[0068] Furthermore, the step 5 comprises:

[0069] Step 5.1, initialize inversion parameters and population:

[0070] Define the parameter vector to be optimized: Let the parameter vector to be optimized (i.e. the parameters related to the formation structure, which is the same as the parameter vector θ in the forward model) em Correspondingly) is X=[x1,x2,…,x D ] T , where D is the dimension of the parameter to be optimized;

[0071] Determine the population size: Set the population size involved in the inversion calculation to m, that is, there are m individuals. The m individuals can be understood as different combinations of initial guess values ​​participating in the inversion process;

[0072] Initialize individual position and velocity: For each individual i (i = 1, 2, ..., m), at the initial time t = 0, define its position vector as X i (0) = [x i1 (0), x i2 (0),…,x iD (0)] T ; If the concept of speed is involved in the subsequent update process (similar to some algorithms based on dynamic update), the speed vector is initialized to V at the same time i (0) = [v i1 (0), v i2 (0),…,v iD (0)] T , the initial velocity can be set to zero vector or random vector according to the specific situation;

[0073] Initialize individual optimal position and global optimal position: For each individual i, initialize its individual optimal position vector (i.e. the optimal position found by the individual during the search process) as P i (0) = [p i1 (0), p i2 (0),…,piD (0)] T , which can be set as the initial position vector X i (0); At the same time, initialize the global optimal position vector of the entire population to P g (0) = [p g1 (0), p g2 (0),…,p gD (0)] T , which can be set as the initial position vector of an individual, such as P g (0) = X1(0);

[0074] Step 5.2, iterative update process;

[0075] Step 5.3, after the maximum number of iterations T max After the iteration update, the final inversion result is the global optimal position vector P in the population. g (T max ), which represents the best estimate of the parameters related to the formation structure obtained by inversion calculation;

[0076] Further, step 5.2 includes:

[0077] Set the number of iterations and related parameters: Determine the maximum number of iterations for the inversion calculation as T max , and set other related parameters, such as inertia weight as w, individual learning factor as c1, and social learning factor as c2;

[0078] Iteratively update individual positions and related parameters: starting from t = 0, perform T max Iteration update;

[0079] Generate a new solution by updating the velocity vector: For each individual i, at the tth iteration, according to the current position vector X i (t) Generate new solutions;

[0080] First, calculate the new velocity vector V according to the velocity update formula i (t+1), such as: ij (t+1)=wυ ij (t)+c1r1(p ij (t)-x ij (t))+c2r2(p gj (t)-x ij (t)); where r1 and r2 are random numbers in the interval [0, 1];

[0081] Then calculate the new position vector X according to the position update formula i (t+1), such as: x ij (t+1)=xij (t)+v ij (t+1);

[0082] Calculate the fitness function value: For each individual i, according to the generated new solution X i (t+1), calculate its fitness function value; the fitness function uses the error between the forward model simulation result and the actual collected data as the measurement standard; let the actual collected preprocessed data be d obs , the data obtained by forward model simulation is d sim , then the fitness function f(X) can be defined as:

[0083] Where N is the number of sensors and M is the length of each sensor data sequence;

[0084] Update individual optimal positions and global optimal positions, including:

[0085] Update the optimal position of the individual: For each individual i, if the fitness function value f(X i (t+1)) is better than the fitness function value f(P i (t)), then update the individual optimal position vector to P i (t+1)=X i (t+1);

[0086] Update the global optimal position: among all individuals, if the fitness function value f(X i (t+1)) is better than the fitness function value f(P g (t)), then update the global optimal position vector to P g (t+1)=X i (t+1).

[0087] Further, step 5.2 includes:

[0088] Set the number of iterations and related parameters: Determine the maximum number of iterations for the inversion calculation as T max , and set the initial temperature to T0 and the cooling coefficient to α (0<α<1);

[0089] Iteratively update individual positions and related parameters: starting from t = 0, perform T max Iteration update;

[0090] Generate new solutions (by adding random perturbations): For each individual i, at the tth iteration, in the current solution X i (t), a new solution X′ is generated by adding a random perturbation i (t), such as x′j (t) = x j (t)+Δx j (t)

[0091] Where Δx j (t) is a random perturbation on the jth parameter, generated by a random method such as normal distribution, such as Δx j (t) = ε·N(0,1), where ε is the disturbance amplitude and N(0,1) represents a normal distribution with a mean of 0 and a standard deviation of 1;

[0092] Calculate the fitness function value: For each individual i, according to the generated new solution X′ i (t), calculate its fitness function value; the fitness function uses the error between the forward model simulation result and the actual collected data as the measurement standard; let the actual collected preprocessed data be d obs , the data obtained by forward model simulation is d sim , then the fitness function can be defined as f(X):

[0093]

[0094] Where N is the number of sensors and M is the length of each sensor data sequence;

[0095] Update individual optimal positions and global optimal positions. Including:

[0096] Update the optimal position of the individual: For each individual i, if the fitness function value f(X′ i (t)) is better than the fitness function value f(P i (t)), then update the individual optimal position vector to P i (t+1)=X′ i (t);

[0097] Update the global optimal position: among all individuals, if the fitness function value f(X′ i (t)) is better than the fitness function value f(P g (t)), then update the global optimal position vector to P g (t+1)=X′ i (t).

[0098] Furthermore, the step 6 comprises:

[0099] Assume that the inversion results obtained through different runs are, R2 = [r 21 , r 22 ,…,r 2D ] T , ..., Rn =[r n1 , r n2 ,…,r nD ] T R1=[r 11 , r 12 ,…,r 1D ] T , where the dimension D is consistent with the dimension of the parameter to be optimized;

[0100] Correlation calculation: For two vectors R i and R j , the calculation formula of its correlation coefficient r is:

[0101] in, is the vector R i The mean of is the vector R j The mean of

[0102] If there are multiple inversion results, let the inversion result set be {R1, R2, …, R n}, first calculate the Pearson correlation coefficient matrix C between the two pairs = [c ij ] n×n , where c ij For R i and R j The Pearson correlation coefficient.

[0103] Then the sum of the correlations between each inversion result and all other inversion results is calculated, that is,

[0104] Finally, the inversion result with the largest sum of correlations is selected as the optimal inversion result output.

[0105] Beneficial effects of the present invention: The stratigraphic structure inversion algorithm based on micro-motion detection proposed in the present invention has many significant beneficial effects compared with traditional stratigraphic detection and related algorithms. In terms of data processing, through the data preprocessing step, the concave-convex transformation is used to remove noise, Butterworth filter filtering and normalization processing are used to effectively improve the data quality. Noise removal can accurately filter out the interference components in the original micro-motion signal and restore the true signal characteristics; filtering further highlights the effective frequency band information; normalization makes the data of different sensors comparable, laying a solid foundation for subsequent precise analysis. The feature extraction link can comprehensively mine multiple features such as mean, variance, spectrum, etc. from the preprocessed data. These rich and targeted features can more carefully characterize the stratigraphic structure characteristics, making the inversion results more in-depth and accurate.

[0106] When establishing the forward model, through operations such as spatial discretization, single-pole model reduction, and Lanczos algorithm to construct orthogonal basis, not only can the electromagnetic relationship between micro-vibration signals and stratum structure be accurately simulated, but also the calculation efficiency is greatly improved, and complex forward calculations can be completed in a relatively short time, providing strong support for rapid inversion. The inversion calculation part uses multiple optimization algorithm ideas for parallel processing, combined with the correlation analysis of the inversion results to select the optimal result, further improving the accuracy and reliability of the inversion, and can more accurately obtain stratum structure information, thereby providing more effective decision-making basis in many fields such as geological exploration and engineering construction, reducing costs and improving work efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0107] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0108] Figure 1 , Schematic diagram of the stratigraphic structure inversion algorithm based on microseismic detection. DETAILED DESCRIPTION

[0109] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with embodiments; it should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention;

[0110] Example 1

[0111] This embodiment provides a stratum structure inversion algorithm based on micro-seismic detection, such as Figure 1 , the algorithm includes:

[0112] Step 1, data acquisition: multiple micro-motion sensors are reasonably arranged in the target detection area to form a sensor array to ensure that it can fully cover and obtain micro-motion signals that sufficiently reflect the characteristics of the formation; the number of arranged micro-motion sensors is N, and the data collected by each sensor at different times constitutes a time series; the data collected by the i-th sensor at the j-th time is x ij , where i = 1, 2, ..., N, j = 1, 2, ..., M, where M represents the length of the time series of data collected by each sensor;

[0113] Step 2, data preprocessing: processing the collected raw micro-motion signal data, including noise removal, filtering, and normalization operations;

[0114] Step 3, feature extraction: extracting features that can effectively characterize information related to the stratigraphic structure from the preprocessed data, including mean features, variance features, and spectrum features; spectrum features include spectrum peak frequency and spectrum bandwidth;

[0115] Step 4: Establish the forward model:

[0116] Step 4.1 defines the relationship between electromagnetic characteristic parameters and stratum structure;

[0117] Assume that the electromagnetic characteristic parameter vector related to the formation structure is θ em =[θ em1 ,θ em2 ,…,θ emD ] T , where D is the dimension of the electromagnetic characteristic parameters;

[0118] Assume that the stratum consists of L layers. For the lth layer (l = 1, 2, ..., L), its conductivity σ l , magnetic permeability μ l and dielectric constant ∈ l The value range and change rule of are as follows:

[0119] Conductivity

[0120] where σ l0 is the base conductivity of the first layer (which can be determined based on factors such as the common rock type in this layer), α l and β l It is a coefficient related to the formation characteristics (which can be obtained through experiments or fitting of existing geological data) and is the depth coordinate relative to a reference plane (such as the surface);

[0121] Magnetic Permeability

[0122] where μ l0 is the base permeability of the first layer (similarly determined by rock type, etc.), γ l and δ l is a coefficient related to formation characteristics, z l is the depth coordinate; γ l

[0123] Dielectric constant

[0124] where ∈ l0 is the base dielectric constant of layer l, and l is a coefficient related to formation characteristics, z l is the depth coordinate; then the electromagnetic characteristic parameter vector θ related to the formation structure em The elements in can be calculated according to the above formula, based on the number of layers and the parameters of each layer;

[0125] Step 4.2 Determine the correspondence between micro-motion detection features and electromagnetic characteristic parameters

[0126] Match the features of the micro-motion detection data after preprocessing and feature extraction with the input of the electromagnetic fast forward algorithm;

[0127] For example, the extracted mean, variance, spectrum characteristics, etc. are transformed and mapped in a certain way (the specific transformation method can be determined according to the actual situation and physical principles) as part of the relevant input parameters in the electromagnetic fast forward algorithm, similar to source terms, boundary conditions, etc.;

[0128] Suppose the micro-motion detection feature input vector after transformation and mapping is Where K is the transformed feature dimension, and its elements are obtained by performing specific processing on the features obtained in the preprocessing and feature extraction steps. The specific processing includes:

[0129] Mean feature conversion: the mean feature extracted for the i-th sensor Convert it into the parameter F related to electromagnetic forward modeling em1 (Assuming it is the first element in the corresponding relationship, the order can be adjusted according to the specific situation), through the following formula:

[0130] Where a1 and b1 are coefficients determined based on experimental data or theoretical analysis, and their purpose is to adjust the mean characteristics of the micro-motion signal to a range and form suitable for electromagnetic forward modeling input;

[0131] Variance feature conversion: Variance feature extracted for the i-th sensor Convert to parameter F em2 , the formula is:

[0132] Similarly, a2 and b2 are determined coefficients used to adjust the variance characteristics to suit electromagnetic forward modeling;

[0133] Spectral feature conversion:

[0134] For the spectrum peak frequency f peak，i , converted to parameter F em3 , the formula is: F em3 =a3f peak，i +b3

[0135] For the spectrum bandwidth B i , converted to parameter F em4 , the formula is: F em4 =a4B i +b4

[0136] Here, a3, b3, a4, and b4 are also coefficients determined according to actual conditions, so that the spectrum characteristics can be reasonably integrated into the input parameters of the electromagnetic forward model;

[0137] Through the above conversion and mapping methods, the characteristics of micro-motion detection data are converted into input parameters suitable for electromagnetic fast forward algorithm;

[0138] Step 4.3, with F em As the input feature vector, θ em As an electromagnetic characteristic parameter vector, the specific steps are as follows:

[0139] (1) Spatial discretization (using unstructured grids)

[0140] Define the problem space:

[0141] Considering the stratum space area involved in the micro-motion detection, set its three-dimensional space range, including its geometric dimensions including length, width, and depth, denoted as V strata

[0142] Determine the distribution of electromagnetic physical quantities that need to be solved, such as electric field intensity E, magnetic field intensity H, etc. in the three-dimensional space of the formation. The distribution of these electromagnetic physical quantities is closely related to the formation structure and micro-motion detection signals;

[0143] Generate an unstructured grid:

[0144] According to the geometric shape of the stratum structure and the complexity of the electromagnetic characteristics, a suitable unstructured grid generation algorithm (such as Delaunay triangulation algorithm, etc.) is used to generate the spatial region V strata Discretize

[0145] The continuous three-dimensional space of the stratum is divided into a series of irregular units (such as tetrahedral units, etc.), and the discretized grid is denoted as Ω strata , the node set in the grid is recorded as N strata , the unit set is recorded as E strata ; For each node i∈N strata , record its spatial coordinates (x i ,y i , z i );

[0146] Establish the finite element linear equation system:

[0147] Based on Maxwell's equations, on an unstructured grid Ω strata The vector finite element method is applied to discretize the above; by using appropriate vector basis functions to expand the electromagnetic field (such as first-order edge element basis functions, etc.), and substituting it into the variational form of Maxwell's equations, after a series of derivations (the detailed derivation process is omitted here, and you can refer to relevant finite element theory books), the finite element linear equations for the unknown electromagnetic quantity coefficients are obtained:

[0148] Ax=b

[0149] Among them, A is a large sparse complex coefficient matrix, whose elements are related to the electromagnetic characteristic parameters of the grid unit (such as conductivity, magnetic permeability, etc.) and the integral calculation of the basis function; x is an unknown vector, which contains the electromagnetic quantity coefficients to be solved at each node; b is a known right-hand side vector, which is related to the source term (such as the equivalent source term of the micro-motion detection signal source, etc.) and the boundary conditions;

[0150] (2) Single-pole model reduction

[0151] Select the optimal point:

[0152] In the frequency range [ω min ,ω max ], analyze the characteristics of the coefficient matrix A of the equation group changing with frequency, and select a single optimal pole p through some optimization criteria (such as minimizing the error criterion, maximizing the convergence speed criterion, etc., which can be determined according to the actual problem). opt ; The selection of this pole should make the subsequent order reduction process able to simplify the solution of the equation system to the greatest extent while ensuring accuracy;

[0153] Transform the system of equations into:

[0154] Using the selected optimization point p opt , transform the original large sparse complex coefficient equation group Ax = b; specifically, by introducing a transformation matrix related to the poles (whose specific form is determined according to the pole model), the original equation group is transformed into a new real number equation group; let the transformed equation group be: A real x real =b real

[0155] Among them, A real is the transformed real coefficient matrix, which has the same dimension as the original matrix A, but the coefficients are real numbers; x real is the corresponding real unknown vector; b real is the transformed real right-hand side vector;

[0156] (3) Lanczos algorithm to construct orthogonal basis

[0157] Initialize the Lanczos algorithm:

[0158] The transformed real coefficient matrix A real As input, set the initial vector υ1 of the Lanczos algorithm; usually you can choose a random vector (the random vector must satisfy the normalization condition ||υ1||=1) or choose an initial vector with a specific meaning (such as a vector related to the source term, etc.) according to the physical characteristics of the problem;

[0159] Construct an orthogonal basis iteratively:

[0160] The iterative process of Lanczos algorithm is performed, for k = 1, 2, ..., m (where m is a pre-set subspace dimension, much smaller than the matrix A real The order of ), iterate according to the following steps:

[0161] Calculate the vector w k =A real υ k ;

[0162] Orthogonalize the vector wk to obtain the orthogonal vector υ k+1 ; The specific orthogonalization process is as follows:

[0163] β k =||w k ||,

[0164] Among them, α k and β k are the coefficients generated during the orthogonalization process;

[0165] After m iterations, a set of orthogonal vectors {υ1,υ2,…,v m+1}, which form an orthogonal basis of the rational Krylov subspace;

[0166] In order to further improve the effectiveness and stability of the orthogonal basis, we can normalize the obtained orthogonal basis; for each orthogonal vector υ i (i=1,2,…,m+1), calculate its norm ||υ i ||, then the vector υ i Divide by its norm to get the normalized orthogonal vector In this way, the orthogonal basis after normalization It can play a better role in subsequent calculations and avoid problems such as numerical instability caused by large differences in vector lengths;

[0167] (4) Calculate orthogonal projection

[0168] Construct the projection matrix:

[0169] Using the normalized orthogonal basis vectors Construct a projection matrix P; the column vectors of the projection matrix P are composed of orthogonal basis vectors, that is;

[0170] Compute the orthogonal projection of the initial matrix:

[0171] Calculate the initial matrix A real Orthogonal projection A onto Krylov subspace proj; According to the definition of projection matrix and the calculation formula of orthogonal projection, we have:

[0172] A proj =P T A real P

[0173] Similarly, for the right-hand vector b real Also project and get b proj =PTb real ;

[0174] (5) Fast forward modeling

[0175] Solve using the projection matrix:

[0176] Now, we use the projection matrix A proj To replace the initial matrix A real Solve it; get the projected system of equations:

[0177] A proj x proj =b proj

[0178] Since the projection matrix A proj The dimension is much smaller than the initial matrix A real So the speed of solving this system of equations will be greatly improved;

[0179] Solve the projected system of equations:

[0180] Use a suitable linear equations solving algorithm (such as Gaussian elimination method in direct method, conjugate gradient method in iterative method, etc., which can be selected according to specific circumstances) to solve the projected equations A poj x proj =b proj , and get the solution vector x proj ;

[0181] Recover the original solution:

[0182] Finally, through the back projection operation, the solution vector x proj Restore the solution vector x to the original space; specifically, according to the properties of the projection matrix P, we have:

[0183] x=Px proj

[0184] In this way, the electromagnetic characteristics of the stratigraphic structure related to the micro-seismic detection are solved, which will serve as part of the basis for the subsequent inversion calculation;

[0185] Step 5, perform inversion calculation: design multiple distributed parallel inversion calculation algorithms by combining different optimization algorithm ideas to find the formation structure parameters that best match the actual collected data;

[0186] Step 6, performing correlation calculation on multiple inversion results obtained in step 5 (for example, by running the step 5 algorithm with different initial settings multiple times) to determine the most accurate and optimal inversion result as output.

[0187] Preferably, step 2 comprises:

[0188] 1. Remove noise

[0189] Wavelet transform is used for denoising; wavelet transform can decompose the signal into wavelet coefficients of different scales. By setting a suitable threshold, the wavelet coefficients smaller than the threshold are set to zero, and then wavelet reconstruction is performed to obtain the denoised signal.

[0190] Assume the original micro-motion signal data is x ij , the wavelet coefficients obtained after wavelet transform are W ij (s, k), where s represents the scale parameter and k represents the translation parameter;

[0191] The threshold λ can usually be determined based on experience or some statistical method, such as using the general threshold formula: Where σ is the estimated value of the noise standard deviation (which can be obtained by analyzing the wavelet coefficients of the original data), NM is the total number of data points (the number of sensors N multiplied by the length of each sensor data sequence M);

[0192] Threshold the wavelet coefficients:

[0193]

[0194] Then the processed wavelet coefficients are transformed into Reconstruct the denoised signal

[0195] 2. Filtering

[0196] By filtering, we can further remove the interference signals in a specific frequency range and highlight the frequency band of the micro-motion signal we are interested in. The Butterworth filter is used for filtering.

[0197] Assume that the discrete time series of the original micro-motion signal data is x ij [n], where n = 0, 1, ..., M-1 represents discrete time points.

[0198] For a low-pass Butterworth filter, the transfer function is:

[0199] where z is a discrete complex variable, N b is the order of the filter, a k is the filter coefficient, which can be calculated through the design formula based on the required cutoff frequency and other parameters.

[0200] Filter the original data to get the filtered signal y ij [n], implemented through discrete convolution: y ij [n] = x ij [n]*h[n]

[0201] where h[n] is the unit impulse response of the filter, which can be obtained by z-inverting the transfer function H(z), that is, converting the discrete complex variable into a form that can be used for calculations.

[0202] 3. Normalization

[0203] Mapping data to a specific interval facilitates subsequent data processing and analysis, such as making data collected by different sensors comparable. The common minimum-maximum normalization method is used.

[0204] Suppose the micro-motion signal data after noise removal and filtering is x ij (Here, in order to simplify the symbols, we still use x ij Indicates that it should actually be the result after the previous step).

[0205] Calculate the minimum value in the data and maximum Represent the minimum and maximum values ​​among all sensors and all times, respectively.

[0206] The minimum-maximum normalization formula is used to normalize the data to the interval [0, 1]:

[0207] Among them, y ij is the normalized data.

[0208] Step 3 includes:

[0209] Mean feature: Normalized data sequence y after preprocessing for the i-th sensor i1 ,y i2 , …, y iM , its mean characteristic The calculation formula is: in Represents the mean feature of the i-th sensor data;

[0210] Variance feature: Normalized data y after preprocessing for the i-th sensor i1 ,y i2 , …, y iMsequence, whose variance characteristic The calculation formula is: in Represents the variance feature of the i-th sensor data;

[0211] Spectrum characteristics:

[0212] The normalized data sequence y after preprocessing of the i-th sensor i1 ,y i2 , …, y iM Perform discrete Fourier transform DFT to obtain its spectrum representation Y i (k), where k = 0, 1, ..., M-1 represents the frequency point, and the discrete Fourier transform formula is:

[0213] Then we can get the spectrum Y i (k) Extract the spectrum peak frequency and spectrum bandwidth.

[0214] Spectrum peak frequency: Among them, the spectrum Y i The frequency point k corresponding to the maximum value in (k) max , f s is the sampling frequency,

[0215] Spectral bandwidth: Spectral bandwidth is defined as the frequency range that contains a certain percentage (e.g. 95%) of the total energy. First, calculate the total energy of the spectrum Then, the spectrum energy is accumulated from the low frequency end until the accumulated energy reaches 95% of the total energy. At this time, the corresponding frequency point is f low，i ; Similarly, the spectrum energy is accumulated from the high frequency end until the accumulated energy reaches 95% of the total energy. At this time, the corresponding frequency point is f high，i Then the spectrum bandwidth B i For: B i =f high，i -f low，i .

[0216] Preferably, step 5 comprises:

[0217] Step 5.1, initialize inversion parameters and population:

[0218] Define the parameter vector to be optimized: Let the parameter vector to be optimized (i.e. the parameters related to the formation structure, which is the same as the parameter vector θ in the forward model) em Correspondingly) is X=[x1,x2,…,x D ] T , where D is the dimension of the parameter to be optimized;

[0219] Determine the population size: Set the population size involved in the inversion calculation to m, that is, there are m individuals. The m individuals can be understood as different combinations of initial guess values ​​participating in the inversion process;

[0220] Initialize individual position and velocity: For each individual i (i = 1, 2, ..., m), at the initial time t = 0, define its position vector as X i (0) = [x i1 (0), x i2 (0),…,x iD (0)] T ; If the concept of speed is involved in the subsequent update process (similar to some algorithms based on dynamic update), the speed vector is initialized to V at the same time i (0) = [υ i1 (0),υ i2 (0),…,υ iD (0)] T , the initial velocity can be set to zero vector or random vector according to the specific situation;

[0221] Initialize individual optimal position and global optimal position: For each individual i, initialize its individual optimal position vector (i.e. the optimal position found by the individual during the search process) as P i (0) = [p i1 (0), p i2 (0),…,p iD (0)] T , which can be set as the initial position vector X i (0); At the same time, initialize the global optimal position vector of the entire population to P g (0) = [p g1 (0), p g2 (0),…,p gD (0)], which can be set as the initial position vector of an individual, such as P g (0) = X1(0);

[0222] Step 5.2, iterative update process;

[0223] Step 5.3, after the maximum number of iterations T max After the iteration update, the final inversion result is the global optimal position vector P in the population. g (T max ), which represents the best estimate of the relevant parameters of the formation structure obtained by inversion calculation.

[0224] Preferably, step 5.2 comprises:

[0225] Set the number of iterations and related parameters: Determine the maximum number of iterations for the inversion calculation as T max, and set other related parameters, such as inertia weight as w, individual learning factor as c1, and social learning factor as c2;

[0226] Iteratively update individual positions and related parameters: starting from t = 0, perform T max Iteration update;

[0227] Generate a new solution by updating the velocity vector: For each individual i, at the tth iteration, according to the current position vector X i (t) Generate new solutions;

[0228] First, calculate the new velocity vector V according to the velocity update formula i (t+1), such as: ij (t+1)=wυ ij (t)+c1r1(p ij (t)-x ij (t))+c2r2(p gj (t)-x ij (t)); where r1 and r2 are random numbers in the interval [0, 1];

[0229] Then calculate the new position vector X according to the position update formula i (t+1), such as: x ij (t+1)=x ij (t)+υ ij (t+1);

[0230] Calculate the fitness function value: For each individual i, according to the generated new solution X i (t+1), calculate its fitness function value; the fitness function uses the error between the forward model simulation result and the actual collected data as the measurement standard; let the actual collected preprocessed data be d obs , the data obtained by forward model simulation is d sim , then the fitness function f(X) can be defined as:

[0231] Where N is the number of sensors and M is the length of each sensor data sequence;

[0232] Update individual optimal positions and global optimal positions, including:

[0233] Update the optimal position of the individual: For each individual i, if the fitness function value f(X i (t+1)) is better than the fitness function value f(P i (t)), then update the individual optimal position vector to P i (t+1)=X i (t+1);

[0234] Update the global optimal position: among all individuals, if the fitness function value f(X i (t+1)) is better than the fitness function value f(P g (t)), then update the global optimal position vector to P g (t+1)=X i (t+1).

[0235] Preferably, step 5.2 may also include:

[0236] Set the number of iterations and related parameters: Determine the maximum number of iterations for the inversion calculation as T max , and set the initial temperature to T0 and the cooling coefficient to α (0<α<1);

[0237] Iteratively update individual positions and related parameters: starting from t = 0, perform T max Iteration update;

[0238] Generate new solutions (by adding random perturbations): For each individual i, at the tth iteration, in the current solution X i (t), a new solution X′ is generated by adding a random perturbation i (t), such as x′ j (t) = x j (t)+Δx j (t)

[0239] Where Δx j (t) is a random perturbation on the jth parameter, generated by a random method such as normal distribution, such as Δx j (t) = ε·N(0,1), where ε is the disturbance amplitude and N(0,1) represents a normal distribution with a mean of 0 and a standard deviation of 1;

[0240] Calculate the fitness function value: For each individual i, according to the generated new solution X′ i (t), calculate its fitness function value; the fitness function uses the error between the forward model simulation result and the actual collected data as the measurement standard; let the actual collected preprocessed data be d obs , the data obtained by forward model simulation is d sim , then the fitness function can be defined as f(X):

[0241]

[0242] Where N is the number of sensors and M is the length of each sensor data sequence;

[0243] Update individual optimal positions and global optimal positions. Including:

[0244] Update the optimal position of the individual: For each individual i, if the fitness function value f(X′ i (t)) is better than the fitness function value f(P i (t)), then update the individual optimal position vector to P i (t+1)=X′ i (t);

[0245] Update the global optimal position: among all individuals, if the fitness function value f(X′ i (t)) is better than the fitness function value f(P g (t)), then update the global optimal position vector to P g (t+1)=X′ i (t).

[0246] Preferably, step 6 comprises:

[0247] Assume that the inversion results obtained through different runs are, R2 = [r 21 , r 22 , ..., r 2D ] T , ..., R n =[r n1 , r n2 , ..., r nD ] T R1=[r 11 , r 12 , ..., r 1D ] T , where the dimension D is consistent with the dimension of the parameter to be optimized;

[0248] Correlation calculation: For two vectors R i and R j , the calculation formula of its correlation coefficient r is:

[0249] in, is the vector R i The mean of is the vector R j The mean of

[0250] If there are multiple inversion results, let the inversion result set be {R1, R2, …, R n}, first calculate the Pearson correlation coefficient matrix C between the two pairs = [c ij ] n×n , where c ij For Ri and R j The Pearson correlation coefficient.

[0251] Then the sum of the correlations between each inversion result and all other inversion results is calculated, that is,

[0252] Finally, the inversion result with the largest sum of correlations is selected as the optimal inversion result output.

[0253] Although the above describes the illustrative specific embodiments of the present invention so that those skilled in the art can understand the present invention, the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, all inventions and creations using the concepts of the present invention are protected.

Claims

1. A stratum structure inversion algorithm based on micro-seismic detection, characterized in that: The algorithm includes: Step 1, data acquisition, includes rationally arranging a sensor array containing multiple micro-motion sensors in the target detection area, which can fully cover and obtain micro-motion signals reflecting the formation characteristics; Step 2, data preprocessing: processing the collected raw micro-motion signal data, including noise removal, filtering, and normalization operations; Step 3, feature extraction: Extract features that can effectively characterize information related to the stratigraphic structure from the preprocessed data, including mean features Variance characteristics Spectral characteristics; Spectral characteristics include spectrum peak frequency f peak，i , spectrum bandwidth B i ; Step 4, establishing a forward model; Step 5, performing inversion calculation: performing inversion calculation by executing at least two distributed parallel inversion calculation models; Step 6, performing correlation calculation on the inversion results obtained in step 5, and determining the most accurate and optimal inversion result as output; Step 4 includes: Step 4.1 Define the relationship between electromagnetic characteristic parameters and stratum structure. Suppose the electromagnetic characteristic parameter vector related to stratum structure is θ em =[θ em1 ,θ em2 ,…,θ emD ] T ; Where D is the dimension of the electromagnetic characteristic parameter; Step 4.2: transform and map the mean, variance, spectrum features, etc. extracted in step 2; Suppose the micro-motion detection feature input vector after transformation and mapping is Where K is the transformed feature dimension, and its elements are obtained by processing the features obtained in the preprocessing and feature extraction steps; Among them, the mean feature conversion is the mean feature extracted for the i-th sensor Converted into parameters related to electromagnetic forward modeling Among them, a1 and b1 are predetermined coefficients; Variance feature transformation is the variance feature extracted for the i-th sensor Convert to parameter Among them, a2 and b2 are predetermined coefficients used to adjust the variance characteristics to suit electromagnetic forward modeling; The spectrum feature conversion is for the spectrum peak frequency f peak，i , converted to parameter F em3 =a3f peak，i +b3; for spectrum bandwidth B i , converted to parameter F em4 =a4B i +b4; where a3, b3, a4, and b4 are also predetermined coefficients used to incorporate the spectrum characteristics into the input parameters of the electromagnetic forward model; Step 4.3, with F em As the input feature vector, θ em As electromagnetic property parameter vector, fast forward modeling algorithm is executed, including: (1) Spatial discretization, defining the problem space: setting the three-dimensional space range, including geometric dimensions such as length, width, and depth, denoted as V strata ; Determine the distribution of electromagnetic physical quantities including electric field intensity E and magnetic field intensity H that need to be solved in the three-dimensional space of the formation; Generate unstructured grid: According to the geometric shape of the stratum structure and the complexity of the electromagnetic characteristics, use a suitable unstructured grid generation algorithm to generate the spatial region V strata Discretize; divide the continuous three-dimensional space of the stratum into irregular unit structures, and the discretized grid is Ω strata , the node set in the grid is recorded as N strata , the unit set is recorded as E strata ; For each node i∈N strata , record its spatial coordinates (x i ,y i , z i ); Establish the finite element linear equations: Based on Maxwell's equations, on the unstructured grid Ω strata The vector finite element method is applied to discretize the above. By using appropriate vector basis functions to expand the electromagnetic field and substituting them into the variational form of Maxwell's equations, the finite element linear equations for the unknown electromagnetic coefficients are obtained: Ax = b; Among them, A is a large sparse complex coefficient matrix, whose elements are related to the electromagnetic characteristic parameters of the grid unit and the integral calculation of the basis function, and the correlation is predefined; x is an unknown vector containing the electromagnetic coefficients to be solved at each node; b is a known right-hand side vector, which is related to the source term and boundary conditions; (2) Single-pole model reduction; select the optimal pole, in the frequency range [ω min ,ω max ], analyze the characteristics of the coefficient matrix A of the equation group with frequency, and select a single optimal pole p through a predefined optimization criterion. opt ; Using the selected optimal pole p opt , transform the original large sparse complex coefficient equation system Ax = b; let the transformed equation system be: A real x real =b real ; Among them, A real is the transformed real coefficient matrix, which has the same dimension as the original matrix A, but the coefficients are real numbers; x real is the corresponding real unknown vector; b real is the transformed real right-hand side vector; (3) Lanczos algorithm constructs an orthogonal basis, including initializing the Lanczos algorithm: using the transformed real coefficient matrix A real As input, set the initial vector v1 of the Lanczos algorithm; usually a random vector that satisfies the normalization condition ‖v1‖=1 can be selected; Gradually iterate to construct an orthogonal basis: for k = 1, 2, ..., m, where m is a pre-set subspace dimension, much smaller than the matrix A real The order of ; perform the following iterative steps: Calculate the vector w k =A real v k ; For vector w k Perform orthogonal processing to obtain the orthogonal vector v k+1 ; The specific orthogonalization process is as follows: β k =‖w k ‖, Among them, α k and β k are the coefficients generated during the orthogonalization process; After m iterations, we get a set of orthogonal vectors {v1, v2, …, v m+1 }, which form an orthogonal basis of the rational Krylov subspace; For each orthogonal vector v i (i=1,2,…,m+1), calculate its norm ||v i ||, then the vector v i Divide by its norm to get the normalized orthogonal vector (4) Calculate the orthogonal projection, including constructing the projection matrix: Use the normalized orthogonal basis vectors obtained Construct a projection matrix P; the column vectors of the projection matrix P are composed of orthogonal basis vectors, that is; Compute the orthogonal projection of the initial matrix: Calculate the initial matrix A real Orthogonal projection A onto Krylov subspace proj ; According to the definition of projection matrix and the calculation formula of orthogonal projection, we have: A proj =P T A real P; Similarly, for the right-hand vector b real Also project and get b proj =P T b real ; (5) Fast forward solution, including the use of projection matrix solution: Using the projection matrix A proj To replace the initial matrix A real Solve it; get the projected system of equations: A proj x proj =b proj ; Since the projection matrix A proj The dimension is much smaller than the initial matrix A real So the speed of solving this system of equations will be greatly improved; Solve the projected system of equations: Use the linear equations solving algorithm to solve the projected system of equations A proj x proj =b proj , and get the solution vector x proj ; Restore the original solution: Through the inverse projection operation, the solution vector x proj Restore the solution vector x to the original space; specifically, according to the properties of the projection matrix P, we have: x=Px proj The solution to the electromagnetic characteristics of the stratigraphic structure related to microtremor detection will serve as the basis for subsequent inversion calculations.

2. The stratum structure inversion algorithm based on micro-seismic detection according to claim 1 is characterized in that: The data collection in step 1 specifically includes: Arrange multiple micro-motion sensors in the target detection area to form a sensor array to ensure that it can fully cover and obtain micro-motion signals that sufficiently reflect the characteristics of the formation; the number of arranged micro-motion sensors is N, and the data collected by each sensor at different times constitutes a time series; the data collected by the i-th sensor at the j-th time is x ij , where i = 1, 2, ..., N, j = 1, 2, ..., M, where M represents the length of the time series of data collected by each sensor.

3. The stratum structure inversion algorithm based on micro-seismic detection according to claim 1 is characterized in that: The data preprocessing includes noise removal, filtering, and normalization operations; The noise removal includes setting the original micro-motion signal data as x ij , the concave-convex coefficient obtained after the concave-convex transformation is W ij (s, k), where s represents the scale parameter and k represents the translation parameter; The threshold λ can usually be determined based on experience or some statistical method, such as using the general threshold formula: Where σ is an estimate of the standard deviation of the noise, NM is the total number of data points, and the number of sensors N is multiplied by the length of each sensor data sequence M; Threshold the bump coefficient: Then the processed concave and convex coefficients are transformed by inverse concave and convex transformation. Reconstruct the denoised signal 4. The stratum structure inversion algorithm based on micro-seismic detection according to claim 1 is characterized in that: In step 3 Mean feature: Normalized data sequence y after preprocessing for the i-th sensor i1 ,y i2 , …, y iM , its mean characteristic The calculation formula is: in Represents the mean feature of the i-th sensor data; Variance feature: Normalized data y after preprocessing for the i-th sensor i1 ,y i2 ,…y iM sequence, whose variance characteristic The calculation formula is: in Represents the variance feature of the i-th sensor data; Spectrum characteristics: The normalized data sequence y after preprocessing of the i-th sensor i1 ,y i2 ,…y iM Perform discrete Fourier transform DFT to obtain its spectrum representation Y i (k), where k = 0, 1, ..., M-1 represents the frequency point, and the discrete Fourier transform formula is: Then we can get the spectrum Y i (k) Extract the spectrum peak frequency and spectrum bandwidth.

5. The stratum structure inversion algorithm based on micro-seismic detection according to claim 1 is characterized in that: The step 5 comprises: Step 5.1, initialize inversion parameters and population: Define the parameter vector to be optimized: Let the parameter vector to be optimized (i.e., the parameters related to the formation structure, and the parameter vector θ in the forward model) be em Correspondingly) is X=[x1,x2,…,x D ] T , where D is the dimension of the parameter to be optimized; Determine the population size: Set the population size involved in the inversion calculation to m, that is, there are m individuals. The m individuals can be understood as different combinations of initial guess values ​​participating in the inversion process; Initialize individual position and velocity: For each individual i (i = 1, 2, ..., m), at the initial time t = 0, define its position vector as X i (0) = [x i1 (0), x i2 (0),…,x iD (0)] T ; If the concept of speed is involved in the subsequent update process (similar to some algorithms based on dynamic update), the speed vector is initialized to V at the same time i (0) = [v i1 (0), v i2 (0),…,v iD (0)] T , the initial velocity can be set to zero vector or random vector according to the specific situation; Initialize individual optimal position and global optimal position: For each individual i, initialize its individual optimal position vector (i.e. the optimal position found by the individual during the search process) as P i (0) = [p i1 (0), p i2 (0),…,p iD (0)] T , which can be set as the initial position vector X i (0); At the same time, initialize the global optimal position vector of the entire population to P g (0) = [p g1 (0), p g2 (0),…,p gD (0)] T , which can be set as the initial position vector of an individual, such as P g (0) = X1(0); Step 5.2, iterative update process; Step 5.3, after the maximum number of iterations T max After the iteration update, the final inversion result is the global optimal position vector P in the population. g (T max ), which represents the best estimate of the relevant parameters of the formation structure obtained by inversion calculation.

6. The stratum structure inversion algorithm based on micro-seismic detection according to claim 5 is characterized in that: Step 5.2 includes: Set the number of iterations and related parameters: Determine the maximum number of iterations for the inversion calculation as T max , and set other related parameters, such as inertia weight as ω, individual learning factor as c1, and social learning factor as c2; Iteratively update individual positions and related parameters: starting from t = 0, perform T max Iteration update; Generate a new solution by updating the velocity vector: For each individual i, at the tth iteration, according to the current position vector X i (t) Generate new solutions; First, calculate the new velocity vector V according to the velocity update formula i (t+1), such as: v ij (t+1)=ωv ij (t)+c1r1(p ij (t)-x ij (t))+c2r2(p gj (t)-x ij (t)); where r1 and r2 are random numbers in the interval [0, 1]; Then calculate the new position vector X according to the position update formula i (t+1), such as: x ij (t+1)=x ij (t)+v ij (t+1); Calculate the fitness function value: For each individual i, according to the generated new solution X i (t+1), calculate its fitness function value; the fitness function uses the error between the forward model simulation result and the actual collected data as the measurement standard; let the actual collected preprocessed data be d obs , the data obtained by forward model simulation is d sim , then the fitness function f(X) can be defined as: Where N is the number of sensors and M is the length of each sensor data sequence; Update individual optimal positions and global optimal positions, including: Update the optimal position of the individual: For each individual i, if the fitness function value f(X i (t+1)) is better than the fitness function value f(P i (t)), then update the individual optimal position vector to P i (t+1)=X i (t+1); Update the global optimal position: among all individuals, if the fitness function value f(X i (t+1)) is better than the fitness function value f(P g (t)), then update the global optimal position vector to P g (t+1)=X i (t+1).

7. The stratum structure inversion algorithm based on micro-seismic detection according to claim 5 is characterized in that: Step 5.2 includes: Set the number of iterations and related parameters: Determine the maximum number of iterations for the inversion calculation as T max , and set the initial temperature to T0 and the cooling coefficient to α (0<α<1); Iteratively update individual positions and related parameters: starting from t = 0, perform T max Iteration update; Generate new solutions (by adding random perturbations): For each individual i, at the tth iteration, in the current solution X i (t), a new solution X′ is generated by adding a random perturbation i (t), such as x′ j (t) = x j (t)+Δx j (t) Where Δx j (t) is a random perturbation on the jth parameter, generated by a random method such as normal distribution, such as Δx j (t) = ε·N(0,1), where ε is the disturbance amplitude and N(0,1) represents a normal distribution with a mean of 0 and a standard deviation of 1; Calculate the fitness function value: For each individual i, according to the generated new solution X′ i (t), calculate its fitness function value; the fitness function uses the error between the forward model simulation result and the actual collected data as the measurement standard; let the actual collected preprocessed data be d obs , the data obtained by forward model simulation is d sim , then the fitness function can be defined as f(X): Where N is the number of sensors and M is the length of each sensor data sequence; Update individual optimal positions and global optimal positions. Including: Update the optimal position of the individual: For each individual i, if the fitness function value f(X′ i (t)) is better than the fitness function value f(P i (t)), then update the individual optimal position vector to P i (t+1)=X′ i (t); Update the global optimal position: among all individuals, if the fitness function value f(X′ i (t)) is better than the fitness function value fP of the current global optimal position g (t)), then update the global optimal position vector to P g (t+1)=X′ i (t).

8. The stratum structure inversion algorithm based on micro-seismic detection according to claim 5 is characterized in that: The step 6 comprises: Assume that the inversion results obtained through different runs are, R2 = [r 21 , r 22 ,…,r 2D ] T , ..., R n =[r n1 , r n2 ,…,r nD ] T R1=[r 11 , r 12 ,…,r 1D ] T , where the dimension D is consistent with the dimension of the parameter to be optimized; Correlation calculation: For two vectors R i and R j , the calculation formula of its correlation coefficient r is: in, is the vector R i The mean of is the vector R j The mean of If there are multiple inversion results, let the inversion result set be {R1, R2, …, R n }, first calculate the Pearson correlation coefficient matrix C between the two pairs = [c ij ] n×n , where c ij For R i and R j The Pearson correlation coefficient. Then the sum of the correlations between each inversion result and all other inversion results is calculated, that is, Finally, the inversion result with the largest sum of correlations is selected as the optimal inversion result output.

Images