Ground stress azimuth earthquake prediction method based on transverse wave velocity change rate

Abstract

The invention relates to the technical field of earthquake prediction, in particular to a crustal stress azimuth earthquake prediction method based on a transverse wave speed change rate, which comprises the following steps: calculating a speed value of a transverse wave and a propagation distance, extracting a difference value to generate a change rate set, calculating an offset rate, screening a continuous overrun section, identifying a disturbance path section, and performing interpolation restoration to generate a purification sequence. Extracting forward and reverse propagation trends and counting direction identifiers, and generating a crustal stress main stress direction map; according to the method, a change rate set is generated through a sliding window segmentation speed sequence, data continuity is enhanced, local features are reserved, a sudden change area is screened and accurately positioned in combination with a continuous offset rate, multi-path lithology difference identification and disturbance repair are carried out, non-structural interference is eliminated, and forward and reverse path cross validation is carried out to screen a direction consistent section, so that the discrimination robustness is improved; the velocity direction is dynamically associated with the structure trend, synchronous analysis of the stress field and fault activity is realized, multi-stage processing is combined with dynamic threshold and direction verification, and the derivation precision and timeliness are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake prediction, and particularly to an earthquake prediction method based on the azimuth of in-situ stress using the change rate of shear wave velocity. Background Art

[0002] The technical field of earthquake prediction includes a system of methods for predicting the time, location, intensity and other characteristics of earthquakes, which is a key research direction at the intersection of seismology and geology. The core contents of this technical field include the monitoring and identification of earthquake precursor phenomena, the statistical analysis of seismic activity, the monitoring of changes in physical parameters of underground media, and the evolutionary analysis of tectonic stress fields. This field usually involves the analysis of seismic wave propagation characteristics, the inversion derivation of stress field distribution, and the dynamic monitoring of underground fault activities, and widely combines information such as geophysical observation data, microseismic activity records, and the change trend of elastic parameters of underground media to construct a multi-factor correlation model for judging the possible earthquake occurrence area and potential risk level.

[0003] Among them, the earthquake prediction method based on the azimuth of in-situ stress using the change rate of shear wave velocity refers to a method that uses the change rate of the underground shear wave propagation velocity over time or space to infer the direction of the principal in-situ stress within the underground rock mass, and uses this as an important basis for judging seismic activity. The technical matters mainly targeted by this patent theme include the identification of the corresponding relationship between the shear wave velocity change characteristics and the in-situ stress field, the extraction method of the principal stress direction, and the analysis of its evolutionary trend in spatial distribution. Specifically, by obtaining regional seismic wave data, comparing the change rate of the shear wave velocity in different propagation paths and time periods, and based on the functional relationship between stress and wave velocity in elastic mechanics theory, the principal stress direction is calculated and derived, thus forming a set of in-situ stress azimuth discrimination processes for earthquake activity prediction.

[0004] The prior art relies on a static inversion model of the relationship between seismic wave propagation characteristics and stress fields, failing to combine dynamic segmentation of data with a sliding time window, resulting in the masking of local stress change characteristics by the overall trend and affecting the timeliness of in-situ stress direction identification. When screening velocity anomaly sections based on a fixed threshold, the verification of the direction consistency of the continuous sliding window offset rate is ignored, making it vulnerable to random noise interference and causing misjudgment, thereby reducing the positioning accuracy of candidate sections. Existing methods lack a mechanism for comparing multi-path lithology parameter differences and repairing disturbed path segments, and are unable to distinguish data fluctuations caused by structural and non-structural factors, leading to systematic biases in the stress field inversion results. Traditional directional analysis is limited to the discrimination of the trend of a single propagation path and does not introduce cross-verification of the forward and reverse path direction identifiers, making it difficult to exclude the influence of the multi-path wave velocity superposition effect and reducing the credibility of the principal stress direction discrimination. The lack of association between the static structural strike and the dynamic velocity direction results in a decoupling risk between the stress field evolution trend and the actual activity state of the fault. For example, in a strike-slip fault area, the offset amount of the principal stress direction may be misjudged, affecting the accuracy of seismic risk level assessment. The prior art relies on the analysis of a single geophysical parameter and fails to construct a multi-factor dynamic association model, restricting the refined analysis of the spatial distribution characteristics of the in-situ stress field in complex tectonic environments. Summary of the Invention

[0005] An object of the present invention is to solve the deficiencies existing in the prior art and to propose a method for seismic prediction of in-situ stress azimuth based on the shear wave velocity change rate.

[0006] To achieve the above object, the present invention adopts the following technical solutions: A method for seismic prediction of in-situ stress azimuth based on the shear wave velocity change rate, comprising the following steps:

[0007] S1: In the monitoring of active faults, calculate the wave velocity value based on the shear wave arrival time recorded by the observation station and the propagation path distance, divide the velocity sequence using a sliding time window with a fixed length, extract the velocity difference between adjacent sliding window segments, and generate a set of shear wave velocity change rate sequences;

[0008] S2: Based on the set of shear wave velocity change rate sequences, extract continuous sliding window segments, calculate the velocity range within the window segment, calculate the offset rate using the median of the range as a reference value, and count the window segments where the offset rate exceeds the upper limit of the range three times continuously and in the same direction, and output candidate sections with velocity mutations;

[0009] S3: Based on the candidate sections with velocity mutations, extract the lithology parameters of the paths within the multi-path intersection area of the shear wave, calculate the lithology difference between adjacent observation points, compare it with the disturbance threshold value, identify the disturbed path segments, and use the velocity difference between adjacent segments for linear interpolation repair to generate a purified shear wave change rate sequence segment;

[0010] S4: Based on the purified shear wave rate of change sequence segment, extract the velocity change trends of the forward and reverse propagation paths, count the direction identifiers of the velocity differences within three consecutive sliding windows, screen the path segments where the forward and reverse directions are consistent, and output the shear wave direction consistency segment identifier.

[0011] As a further solution of the present invention, the shear wave velocity rate of change sequence set includes the fixed-length sliding window division of the velocity sequence, the extraction of the velocity differences between adjacent sliding window segments, and the generation of the shear wave velocity rate of change classification. The velocity mutation candidate segment includes the construction of the maximum and minimum velocity differences in the range interval, the determination of the median reference value, the calculation of the sliding window offset rate, and the statistics of the offset rates of three consecutive segments. The purified shear wave rate of change sequence segment includes the extraction of lithology parameters, the comparison of the lithology difference perturbation threshold, the identification of the perturbed path segment, and the linear interpolation repair of the velocity differences between adjacent segments. The shear wave direction consistency segment identifier includes the extraction of the velocity trends of the forward and reverse propagation paths, the statistics of the direction identifiers of the velocity differences, and the screening of the path segments with consistent directions.

[0012] As a further solution of the present invention, the specific steps of S1 are as follows:

[0013] S101: Based on the shear wave arrival time and propagation path distance data recorded by the observation stations in the active fault monitoring area, calculate the shear wave propagation velocity values corresponding to the propagation paths, and integrate the calculation results of all paths to generate a shear wave propagation velocity set;

[0014] S102: Call the shear wave propagation velocity set, slide and divide continuous window segments along the time axis according to the preset window length, perform mean processing on the velocity values within each window segment, and generate a sliding window velocity sequence set;

[0015] S103: Extract the mean data of adjacent window segments in the sliding window velocity sequence set, calculate the velocity differences between window segments, group and classify the differences according to the propagation path number, and generate a shear wave velocity rate of change sequence set.

[0016] As a further solution of the present invention, the specific steps of S2 are as follows:

[0017] S201: Based on the shear wave velocity rate of change sequence set, extract consecutive sliding window segments with a fixed length, summarize the velocity rate of change differences within the segments, sort them according to the difference size, locate the median position and extract the corresponding value to generate a velocity difference median reference value;

[0018] S202: According to the velocity difference median reference value, calculate the offset rate of the velocity rate of change differences of the sliding window segments to form a complete offset rate sequence, screen the segments greater than the upper limit of the difference set and extract the offset direction to generate a super-threshold offset direction sequence;

[0019] S203: Based on the sequence of super-threshold offset directions, perform consistency judgment on three consecutive offset directions, extract the number intervals of the sliding window segments that meet the conditions and integrate them into a set of number ranges, and generate candidate sections of velocity mutation.

[0020] As a further solution of the present invention, the specific steps of S3 are as follows:

[0021] S301: Based on the candidate sections of velocity mutation, extract the velocity data of the paths within the shear wave multipath intersection area and the formation structure information of the corresponding measuring points, classify and combine the lithology density, porosity and elastic modulus of the path segments, and generate the set value of the lithology parameters of the path segments;

[0022] S302: According to the set value of the lithology parameters of the path segments, compare the difference in lithology parameters between adjacent observation points with the set disturbance threshold value, identify the numbers of the path segments that exceed the threshold, and generate a difference sequence of disturbed path segments after integration;

[0023] S303: Call the difference sequence of the disturbed path segments, repair the disturbed path segments with reference to the velocity difference between adjacent undisturbed path segments, and update the shear wave change rate of the observation points in the path segments to obtain a purified shear wave change rate sequence segment.

[0024] As a further solution of the present invention, the specific calculation formula for comparing the difference in lithology parameters between adjacent observation points with the set disturbance threshold value is:

[0025] ;

[0026] Wherein, represents the dynamic difference value of the lithology parameters between the i-th and (i - 1)-th observation points, represents the standardized value of the k-th type of lithology parameter at the i-th observation point, represents the standardized value of the k-th type of lithology parameter at the (i - 1)-th observation point, represents the weight coefficient of the j-th type of environmental factor, represents the normalized influence factor of the j-th type of environmental factor at the i-th observation point, represents the spatial distance between the i-th and (i - 1)-th observation points, represents the path segment curvature correction coefficient, represents the distance attenuation exponent.

[0027] As a further solution of the present invention, the specific steps of S4 are as follows:

[0028] S401: Based on the purified shear wave change rate sequence segment, extract the velocity change direction of the shear wave propagation path in the time window, identify the change of the path direction within adjacent windows, and generate a velocity change direction sequence;

[0029] S402: Call the speed change direction sequence, extract the path direction sequence in a three-segment sliding window manner, identify whether the directions within each group of windows are consistent, record the path segment information that meets the conditions, and obtain the label set of consecutive sliding window direction consistent intervals;

[0030] S403: According to the label set of consecutive sliding window direction consistent intervals, filter the shear wave propagation paths with consistent directions, identify the path segments with constant directions, and generate shear wave direction consistent section identifiers.

[0031] As a further solution of the present invention, the specific calculation formula for filtering the shear wave propagation paths with consistent directions is:

[0032] ;

[0033] where, represents the mean value of all shear wave propagation direction angles within the i-th sliding window, n represents the total number of windows in the label set of consecutive sliding window direction consistent intervals, represents the balance factor based on the dynamic distribution of the shear wave propagation direction angle amplitude, and i represents the index number of the sliding window.

[0034] As a further solution of the present invention, the method further includes:

[0035] S5: Based on the shear wave direction consistent section identifier, extract the included angle between the section speed change direction and the structural strike, construct a path unit vector through angle correction, and integrate the direction vectors of all paths to form a map of the in-situ stress principal stress direction;

[0036] The map of the in-situ stress principal stress direction includes the calculation of the included angle between the speed change direction and the structural strike, the correction and construction of the path unit vector, and the integration of the direction vectors.

[0037] As a further solution of the present invention, the specific steps of S5 are as follows:

[0038] S501: Based on the shear wave direction consistent section identifier, extract the speed change directions of the shear wave paths within the section, combine the strike azimuth information of the structural belt, analyze the included angle relationship between the path direction and the structural strike, and generate the included angle value between the path and the structural strike;

[0039] S502: According to the included angle value between the path and the structural strike, perform direction correction on the path direction vector, unify the numerical structure of the direction vector, and integrate all path results to obtain the path unit vector set data;

[0040] S503: Based on the path unit vector set data, extract the spatial distribution trend of the direction vectors, calculate the angle mean value after classifying and counting by direction, and construct a map of the in-situ stress principal stress direction angle.

[0041] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0042] In the present invention, by generating a change rate set through sliding time window segmentation of the velocity sequence, the data continuity is enhanced and local features are retained. Combining the continuous offset rate screening to accurately locate the mutation area, identifying the lithology differences of multiple paths and repairing the disturbances, excluding non-structural interferences, cross-verifying the forward and reverse paths to screen the sections with consistent directions, improving the discrimination robustness, dynamically associating the velocity direction with the structural strike, realizing the synchronous analysis of the stress field and fault activities, and combining multi-stage processing with dynamic threshold and direction verification to improve the derivation accuracy and timeliness. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 It is a schematic diagram of the step flow of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0044] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the drawings and 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.

[0045] In the description of the present invention, it should be understood that the orientation or positional relationship indicated by the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be construed as limiting the present invention. In addition, in the description of the present invention, "a plurality of" means two or more, unless otherwise specifically defined.

[0046] Please refer to Figure 1 , the in-situ stress azimuth seismic prediction method based on the shear wave velocity change rate, comprising the following steps:

[0047] S1: In the active fault monitoring area, obtain the shear wave arrival time and propagation path distance of the observation stations, calculate the shear wave velocity value, divide the velocity sequence by a sliding time window with a fixed length, extract the velocity difference between adjacent sliding window segments, and classify and generate a shear wave velocity change rate sequence set;

[0048] S2: Based on the shear wave velocity change rate sequence set, extract continuous sliding window segments, calculate the maximum and minimum velocity differences within each group of segments, construct a range interval and take the median as a reference value, calculate the offset rate between the sliding window and the median, and count the sliding window segments with the offset rate exceeding the upper limit of the range and the same direction in three consecutive segments, and output the candidate sections of velocity mutation;

[0049] S3: Based on the velocity mutation candidate section, extract the lithology parameters of the paths within the shear wave multipath intersection area, calculate the lithology difference between adjacent observation points, compare it with the perturbation threshold value, identify the perturbed path segments, and use the velocity difference between adjacent segments for linear interpolation repair to generate a purified shear wave rate of change sequence section;

[0050] S4: Based on the purified shear wave rate of change sequence section, extract the velocity change trends of the forward and backward propagation paths, count the velocity difference direction identifiers within three consecutive sliding windows, screen the path segments where the forward and backward directions are consistent, and output the shear wave directionality consistent section identifier;

[0051] S5: Based on the shear wave directionality consistent section identifier, extract the angle between the velocity change direction of the section and the tectonic strike, construct the path unit vector through angle correction, and integrate the direction vectors of all paths to form the map of the principal stress direction of in-situ stress.

[0052] The shear wave rate of change sequence set includes the fixed-length sliding window division of the velocity sequence, the extraction of the velocity difference between adjacent sliding window segments, and the classification generation of the shear wave rate of change. The velocity mutation candidate section includes the construction of the maximum and minimum velocity differences in the range interval, the determination of the median reference value, the calculation of the sliding window offset rate, and the statistics of the offset rates for three consecutive segments. The purified shear wave rate of change sequence section includes the extraction of lithology parameters, the comparison of the lithology difference perturbation threshold, the identification of perturbed path segments, and the linear interpolation repair of the velocity difference between adjacent segments. The shear wave directionality consistent section identifier includes the extraction of the velocity trends of the forward and backward propagation paths, the statistics of the velocity difference direction identifiers, and the screening of path segments with consistent directions. The map of the principal stress direction of in-situ stress includes the calculation of the angle between the velocity change direction and the tectonic strike, the correction and construction of the path unit vector, and the integration of the direction vectors.

[0053] The specific steps of S1 are as follows:

[0054] S101: Based on the shear wave arrival time and propagation path distance data recorded by the observation stations within the active fault monitoring area, calculate the shear wave propagation velocity value corresponding to the propagation path, and integrate the calculation results of all paths to generate a shear wave propagation velocity set;

[0055] In the active fault monitoring area, after the observation stations are deployed, it is necessary to collect and mark their specific spatial coordinates. For example, the station numbered A is located at a certain longitude and latitude, and the station numbered B is located in a similar longitude and latitude area, forming a station network to cover the target monitoring range. In an actual earthquake event, the arrival time of the shear wave corresponding to a specific station is extracted by the seismic wave recording equipment. For example, in a certain event, the shear wave is recorded to arrive at Station A at 13.45 seconds. At the same time, based on the earthquake source location information and the geographic coordinate system, the straight-line propagation distance from the epicenter to the station can be calculated. For example, this distance is 10.2 kilometers. Based on this, the shear wave propagation speed of the corresponding path can be expressed as the propagation path distance divided by the shear wave arrival time, obtaining a value representing the shear wave propagation speed of the path. Continue to extract data with other earthquake events and other station combinations. For example, the measured distance from the epicenter to Station B is 15.3 kilometers, and the arrival time is 18.6 seconds. After processing, the corresponding speed of this path is obtained. The arrival times and propagation path distances collected in each event are uniformly converted into speed values according to the above method, and after summarization, a shear wave propagation speed set is formed. This set contains the speed values corresponding to multiple paths and multiple event combinations, constituting the basic sample set of the shear wave propagation speed in this monitoring area.

[0056] S102: Call the shear wave propagation speed set, slide along the time axis according to the preset window length to divide continuous window segments, and perform mean processing on the speed values within each window segment to generate a sliding window speed sequence set;

[0057] After calling the shear wave propagation speed set formed above, it is necessary to set a sliding window segment with a fixed length on the time axis and perform window division processing by setting the sliding step size. For example, if the length of each window segment is set to 60 seconds and the sliding step size is set to 30 seconds, then multiple overlapping window segments are generated by advancing along the time axis. Each window segment contains several speed records. For example, a certain window segment contains 10 speed values, and these 10 values need to be averaged to obtain the average propagation speed within this time period. According to this method, all sliding window segments are processed one by one to generate a complete window average speed sequence. During the operation, for the situation where there are too few speed data in some window segments, a minimum sample number limit is set. For example, when the number of records is less than 3, the data of this segment is discarded to ensure the representativeness of the mean value. The window division strategy can adopt the sliding ratio setting method. For example, the sliding step size is half of the window length to ensure coverage continuity. After the above steps of processing, a complete sliding window speed sequence set is formed, providing a stable data sequence basis for subsequent dynamic change analysis.

[0058] S103: Extract the mean data of adjacent window segments in the sliding window speed sequence set, calculate the speed difference between window segments, and group and classify the differences according to the propagation path number to generate a shear wave speed change rate sequence set;

[0059] After extracting the means of two adjacent window segments in the sliding window speed sequence, calculate the speed difference between them. For example, if the mean of the previous window is 0.803 and the mean of the next window is 0.790, the speed difference between the two segments is negative. Take the absolute difference as the analysis basis. Classify and count these differences according to the propagation path number. For example, if there are 10 sets of speed change records under Path-001, summarize them into the classification of this path. Subsequently, perform continuity processing on the change sequence under this path. The three-point moving average method can be used to process each difference in the sequence to eliminate the influence of short-term abnormal fluctuations and improve data stability. When dividing whether it is a significant change, set the threshold boundary of the significant change. For example, in historical statistics, the speed change rates of a large number of paths are concentrated below 0.018, and 0.020 can be set as the boundary point to distinguish whether it is significant. If the change value of any path within a certain period exceeds this boundary, it is regarded as a strong fluctuation phenomenon of this path during this period. Finally, summarize the speed differences between each segment in all paths to form a set of speed change rate sequences organized by path number, which serves as the analysis basis for the dynamic characteristics of shear wave propagation.

[0060] The specific steps of S2 are as follows:

[0061] S201: Based on the shear wave speed change rate sequence set, extract continuous sliding window segments of a fixed length, summarize the speed change rate differences within the segments, sort them according to the difference size, locate the median position, and extract the corresponding value to generate the median reference value of the speed difference.

[0062] Extract equal-length sliding window segments from the continuous data sequence of the shear wave speed change rate. First, the window length needs to be determined. For example, select every five observation points to form a segment and slide with a step size of one observation point to continuously extract these sliding window segments from the entire speed change rate sequence. Within each segment, the recorded shear wave speed change rate values are arranged in chronological order, and the difference between the maximum and minimum values is calculated one by one to obtain the speed change rate difference of this segment. For example, if the data in a certain segment is 0.03, 0.02, 0.04, 0.01, 0.05 in sequence, the difference is the maximum value 0.05 minus the minimum value 0.01, which is equal to 0.04. After collecting the difference sets of all sliding window segments, sort them from smallest to largest, and select the value at the middle position as the reference. If the total number of data is even, take the average of the two middle numbers. For example, if the summary result of a certain time is 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, the median position is between the third and fourth items, and the reference value is 0.035. This value is used as a reference standard for subsequent comparison and evaluation.

[0063] S202: According to the median reference value of the speed difference, calculate the offset rate of the speed change rate difference of the sliding window segment to form a complete offset rate sequence, screen the segments greater than the upper limit of the difference set, and extract the offset direction to generate a sequence of offset directions exceeding the threshold.

[0064] After the reference value is determined, traverse each sliding window segment again to calculate the offset rate of the difference between its velocity change rate and the reference value. The offset rate is defined as the growth or reduction ratio of the current segment difference relative to the reference value. If the difference of a certain segment is 0.05 and the reference value is 0.035, its offset rate is expressed as the ratio of the deviation of this segment to the reference value, that is, this difference is about 42% higher than the reference value. After calculating all the sliding window segments in sequence, a complete offset rate sequence is formed. To identify significant abnormal changes, a critical value for judgment needs to be set, and this value can be set based on the distribution characteristics of the difference set. For example, a combination of the upper quartile and the interquartile range is used. The segments exceeding this upper limit are regarded as abnormal, and their offset directions are recorded. If the difference is higher than the reference value, the direction is marked as positive, and if it is lower, it is negative. Finally, the offset directions of all segments exceeding the threshold are summarized in chronological order to obtain a direction sequence for subsequent analysis.

[0065] S203: Based on the sequence of offset directions beyond the threshold, judge the consistency of the offset directions of three consecutive segments, extract the number intervals of the sliding window segments that meet the conditions and integrate them into a set of number ranges to generate candidate sections of velocity mutation;

[0066] In the obtained sequence of offset directions, check segment by segment whether there is a situation where the directions are continuously consistent. If the direction marks of more than three consecutive segments are all positive or all negative, it can be determined that there is a stable consistent trend in this sequence. At this time, extract the numbers of the sliding window segments corresponding to these continuously consistent directions and merge them into a number range. For example, if the direction sequence is positive, positive, positive, negative, negative, positive, positive, then numbers 1 to 3 and 6 to 7 can be regarded as two segments with continuous offset directions respectively. When judging, it is necessary to ensure that the length of the direction mark reaches at least the set minimum number of segments to avoid judgment deviation caused by accidental errors. After the number interval is mapped to the actual time, candidate regions of velocity mutation are generated for subsequent spatial positioning or time series research. This process is based on the numbers of the sliding window segments, traverses and integrates all sequence intervals that meet the consistency conditions in sequence to form the final set of candidate sections of velocity mutation.

[0067] The specific steps of S3 are as follows:

[0068] S301: Based on the candidate sections of velocity mutation, extract the velocity data of the paths within the shear wave multipath intersection region and the formation structure information of the corresponding measuring points, classify and combine the lithology density, porosity and elastic modulus of the path segments to generate a set value of lithology parameters for the path segments;

[0069] After identifying the path area with a significant mutation in the recognition speed, first, according to the spatial variation of the shear wave propagation speed in the seismic data, the part of the path where the speed change rate exceeds the preset critical value is screened. For example, when the shear wave speed change rate is higher than 15%, it can be regarded as the speed mutation section. Subsequently, the data of each observation point covered by the shear wave propagation path are extracted from these sections, and the corresponding formation structure information of these points is obtained from the existing geological profiles or drilling data. Such structure information usually includes the rock layer name, rock layer thickness, burial depth, and other physical properties. Next, combined with the typical lithology parameters recorded in the geological database, the lithology category of the formation where each observation point is located is classified, and the corresponding basic parameters such as density, porosity, and elastic modulus are extracted according to the classification results. For example, a path may contain two formations, a clay layer and a sand layer, which are classified into the corresponding categories through lithology comparison, and their physical parameters are obtained by looking up the table respectively. If the lithology of a certain observation point is silty clay, its typical density can be obtained as 1.85 grams per cubic centimeter, porosity as 42%, and elastic modulus as 150 MPa; while the sand layer may have a density of 2.10 grams per cubic centimeter, porosity of 35%, and elastic modulus of 320 MPa. In the process of constructing the lithology parameter set of the path segment, it is also necessary to supplement the missing or abnormal observation point parameters, usually by using the method of linear interpolation, which is derived by averaging the parameters of adjacent measurement points to ensure the integrity of the parameter set and the logical consistency between data. The finally formed parameter set will be used as the basis for subsequent difference identification and perturbation analysis.

[0070] S302: According to the lithology parameter set value of the path segment, compare the difference in lithology parameters between adjacent observation points with the set perturbation threshold value, identify the path segment numbers exceeding the threshold, and generate a perturbation path segment difference sequence after integration;

[0071] The specific calculation formula for comparing the difference in lithology parameters between adjacent observation points with the set perturbation threshold value is:

[0072] ;

[0073] where represents the dynamic difference value of lithology parameters between the i-th and i - 1-th observation points, represents the standardized value of the k-th type of lithology parameter of the i-th observation point, represents the standardized value of the k-th type of lithology parameter of the i - 1-th observation point, represents the weight coefficient of the j-th type of environmental factor, represents the normalized influence factor of the j-th type of environmental factor of the i-th observation point, represents the spatial distance between the i-th and i - 1-th observation points, represents the path segment curvature correction coefficient, represents the distance attenuation exponent;

[0074] Normalized value of lithology parameter: and are the normalized values of the lithology parameter of the k-th type at the i-th and (i - 1)-th observation points, obtained by normalizing the actual monitoring data.

[0075] Example: The original value of the resistivity of sandstone at the i-th observation point in a certain area is 150 Ω·m. After normalization (formula , where Ω·m, Ω·m), we get ; The resistivity of sandstone at the (i - 1)-th observation point is 120 Ω·m, and after normalization, we get .

[0076] Weight coefficient of environmental factor: is the weight coefficient of the j-th type of environmental factor, determined based on the evaluation of geological experts and the results of principal component analysis.

[0077] Example: The environmental factors in the area include rock humidity (j = 1), slope (j = 2), and fracture density (j = 3). The weights are set as , , respectively, based on the fact that the influence of humidity on lithology stability accounts for 40%, and the slope and fractures each account for 30%.

[0078] Normalized influence factor: represents the normalized result of the measured value of the j-th type of environmental factor at the i-th observation point.

[0079] Example: The humidity sensor measures the humidity at the i-th point as 25%, and after normalization ( ); The slope is measured as 15°, and after normalization ( ); The fracture density is 8 pieces / m², and after normalization ( ).

[0080] Spatial distance: is the spatial distance between the i-th and (i - 1)-th observation points, calculated by the Euclidean distance from the GPS coordinates.

[0081] Example: The coordinate differences between two points are Δx = 30 m and Δy = 40 m, then m.

[0082] Curvature correction coefficient: , obtained by fitting based on historical path curvature data (when the curvature increases by 0.1 m -1 , increases by 0.05, in the interval [0.5, 1.5]); Distance attenuation exponent: , in the interval [0.3, 0.7].

[0083] Formula calculation and derivation:

[0084] Calculate the absolute value of the difference in lithology parameters:

[0085] ;

[0086] Calculate the weighted sum of environmental factors:

[0087] ;

[0088] Calculate the numerator part:

[0089] ;

[0090] Calculate the denominator part:

[0091] ;

[0092] Calculate the dynamic difference value:

[0093] ;

[0094] Parameter summary: represents the dynamic difference value of lithology parameters, and is the standardized lithology parameter, is the weight of the environmental factor, is the measured value of the normalized environmental factor, is the spatial distance, is the curvature correction coefficient, is the distance attenuation exponent;

[0095] Analysis of the example results: It shows that the difference in lithology parameters between the i-th and (i - 1)-th observation points is in a low disturbance state after environmental and spatial corrections. If the disturbance threshold is set to 0.05, then this path segment is not marked as a disturbed segment and needs to be excluded when integrated into the sequence.

[0096] S303: Call the difference sequence of the disturbed path segment, repair the disturbed path segment with the velocity difference of the adjacent undisturbed path segment as a reference, and update the shear wave change rate of the observation points in the path segment to obtain a purified shear wave change rate sequence segment;

[0097] After obtaining the difference sequence of the disturbed path segments, it is necessary to first determine the actual position of each disturbed segment and find its adjacent undisturbed path segments to obtain the reference value. For the reference path segment, the change in the shear wave velocity between two adjacent observation points can be calculated as the velocity change reference for repairing the velocity information of the disturbed path segment. During the repair process, the velocity of each observation point within the disturbed segment needs to be estimated according to the velocity change trend of the reference segment. The estimation method can be to start from an adjacent undisturbed point and derive the velocity value based on the change rate and the distance to the disturbed point. For example, when the velocities of two points on the left reference path segment are 850 m / s and 870 m / s respectively, and the distance between them is 20 m, the derived change rate is 1 m / s per meter increase. Then, if the disturbed point is 15 m away from the left reference point, the velocity of this point can be estimated to be 875 m / s. And so on, the velocities of all observation points within the entire disturbed segment are repaired to form a new velocity sequence. The change in the repaired velocity can be analyzed by taking differences again to ensure that the overall trend is smooth and there are no sudden changes. This process is applicable to data cleaning and noise rejection scenarios, especially in path segments where seismic data has discontinuities, strong interference, or local abnormal signals, and has a prominent application effect. The purified shear wave change rate sequence will be used in the subsequent modeling and analysis stages, ensuring high data quality and consistency.

[0098] The specific steps of S4 are as follows:

[0099] S401: Based on the purified shear wave change rate sequence segment, extract the velocity change direction of the shear wave propagation path in the time window, identify the change in the path direction within adjacent windows, and generate a velocity change direction sequence;

[0100] To obtain the change rate information generated when the seismic shear wave propagates in the sensor array, it is necessary to extract a specific frequency band from the original seismic waveform, such as the fluctuation data with a frequency range between 5 and 20 Hz. The original waveform is first processed by denoising and filtering to eliminate environmental interference and measurement errors. Then, the propagation delay change rate is calculated based on the waveform arrival time difference between multiple sensors and the distance between the sensors. The processed change rate sequence needs to be unified in dimension and made comparable among multiple paths through standardization means. To explore the velocity change trend in the path direction, a time window of a certain length can be set, and the change rate sequence within the window is aggregated and its change direction in the spatial path is analyzed. The change in the propagation velocity of each path during this time period can be transformed into a direction vector, which represents the velocity change trend of the path. When the time window slides forward to an adjacent time period, the change in the included angle of the direction vectors of the same path between adjacent windows needs to be compared. If the included angle is greater than a specific angle threshold, such as 30 degrees, it is judged that the direction has changed, otherwise it is considered that the direction remains the same. The direction change situations of all paths in consecutive time periods are recorded one by one, and finally a complete velocity direction change sequence is formed, providing a basis for subsequent analysis.

[0101] S402: Call the speed change direction sequence, extract the path direction sequence using a three-segment sliding window method, identify whether the directions within each group of windows are consistent, record the information of the path segments that meet the conditions, and obtain the label set of consecutive sliding window direction-consistent intervals;

[0102] After importing the time series of path direction changes, it is divided into three adjacent but continuously sliding time periods, which can be set to 10 seconds each, with a sliding interval of 5 seconds. For each path, its direction features are extracted within each period of time, and the path directions within these three periods are compared. If the angular difference between the path directions of any two of the three segments does not exceed a set range, such as 20 degrees, then record this path segment as having a consistent direction. During the successive forward movement of the sliding window, all path segments that meet the direction-consistent conditions are recorded. Each record contains content such as the path number, direction features, and time interval. Then, according to the continuity on the time axis, overlapping or adjacent records are merged to form multiple sets of time intervals with consistent directions. The corresponding path segment identifiers are summarized within each set to obtain the label set of time segments with consistent path directions, providing support for path screening.

[0103] S403: According to the label set of consecutive sliding window direction-consistent intervals, screen the shear wave propagation paths with consistent directions, identify the path segments with constant directions, and generate the identifier of the shear wave direction-consistent section;

[0104] The specific calculation formula for screening the shear wave propagation paths with consistent directions is:

[0105] ;

[0106] where, represents the mean value of all shear wave propagation direction angles (unit: radian) within the i-th sliding window, n represents the total number of windows in the label set of consecutive sliding window direction-consistent intervals, represents the balance factor based on the dynamic distribution of the shear wave propagation direction angle amplitude (defined as ), and i represents the index number of the sliding window (from 1 to n);

[0107] Parameter definition and data source:

[0108] : Record the mean value of the shear wave propagation direction angle (unit: radian) within the sliding window through seismic monitoring equipment. The monitoring equipment records the waveform at a sampling frequency of 1000 times per second and calculates the arithmetic mean of all shear wave direction angles within the window. For example, the first window measures radians, the second window measures radians, and so on.

[0109] n: The total number of windows in the set of intervals with consistent continuous sliding window directions. It is calculated based on the length of the original seismic signal (e.g., 120 seconds in duration) and the sliding window length (e.g., 10 seconds per window, 5 - second step). .

[0110] : The balance factor is defined as . It is calculated through monitoring data , for example radians, and substituting it gives , because , finally taking .

[0111] Calculation process of the example:

[0112] Taking n = 5 windows as an example, assume The measured values are {0.52, 0.48, 0.50, 0.49, 0.51} radians:

[0113] Calculate radians;

[0114] Calculate ;

[0115] Calculate the numerator ;

[0116] Calculate the first term ;

[0117] Calculate the second term ;

[0118] Combine the results ;

[0119] Basis for parameter setting and numerical rationality:

[0120] The threshold of (about 31.4 radians) is based on the statistical empirical value of the angular fluctuation range of the seismic shear - wave direction (in most scenarios radians), ensuring a reasonable dynamic adjustment range.

[0121] The sliding window length of 10 seconds and the step of 5 seconds refer to the seismic signal processing standard (such as the method recommended by USGS), balancing the time resolution and calculation efficiency.

[0122] Relevance between results and steps:

[0123] In the example, H = 0.7947. When H≥0.75, it is determined that the direction consistency index meets the standard, and the corresponding shear - wave direction - consistency section identifier is generated. This result indicates that the fluctuation of the shear - wave propagation direction within the current window sequence is small, meeting the requirements of direction constancy.

[0124] The specific steps of S5 are as follows:

[0125] S501: Based on the identification of the section with consistent shear wave directionality, extract the velocity change direction of the shear wave paths within the section. Combine the strike azimuth information of the structural belt to analyze the angular relationship between the path direction and the structural strike, and generate the angle value between the path and the structural strike;

[0126] After obtaining the section with consistent shear wave directionality, it is necessary to conduct a directionality consistency analysis on the shear wave signal data obtained along the survey line. Usually, the shear wave waveforms received by multiple seismic measurement points are used. By comparing the main direction change trends of the waveforms between adjacent measurement points, sections with good continuity and consistent direction are identified. This type of analysis can be completed with the help of seismic data processing platforms such as SeisSpace or Petrel. Within the identified section, multiple shear wave propagation paths are selected. By using the picking results of the first arrival time and the geometric length of the paths, calculate the velocities of the shear wave on each path. Then, based on the change direction of these velocity values between different paths, judge the propagation trend direction of the shear wave within the section. If the velocity of a certain path is relatively larger, its propagation direction is more likely to be preferentially selected by the shear wave. For example, assume that the shear wave propagation velocity between measurement points A and B is 2900 m / s, and the velocity between A and C is 3100 m / s. Then it is considered that the shear wave tends to propagate along the direction from A to C. Next, introduce the strike information of the structural belt. Based on the azimuth angle of the known structural belt, compare the path direction with the structural strike and measure the angular relationship between them. Calculate the angle value through the spatial relationship between the path direction coordinates and the structural strike direction. For example, if the path direction is northeast and the structural strike is due east, the angle between the two is 45 degrees. In this way, generate the angle value between the path and the structural strike for subsequent direction adjustment.

[0127] S502: According to the angle value between the path and the structural strike, correct the direction of the path direction vector, unify the numerical structure of the direction vector, and integrate all the path results to obtain the path unit vector set data;

[0128] The included angle data between the path and the structural strike can be used to adjust the directions of each path, so that the direction information of all paths has the same direction expression structure in the unified coordinate system. First, all the included angles between the paths and the structural strike need to be divided into multiple intervals. For example, it is set that 0 to 30 degrees is an interval, 30 to 60 degrees is another interval, and so on. Then, select the path direction with the smallest included angle in each interval as the standard direction of this interval. Next, all the path directions in this interval need to be adjusted to be consistent with the standard direction. The adjustment method can be achieved through vector projection or vector rotation to make its direction consistent with the standard vector. For example, in a direction interval, if the selected standard direction is 45 degrees east of north, other direction vectors need to be adjusted to the same direction. Subsequently, all the corrected direction vectors are standardized, that is, converted into unit vectors, indicating that the path only has the direction attribute and no distance influence, forming a complete set of unit path vectors. Each vector in the set has a unified length of 1 and only records the direction, laying a foundation for subsequent statistical analysis.

[0129] S503: Based on the data of the set of path unit vectors, extract the spatial distribution trend of the direction vectors, calculate the angle mean value after classifying and counting by direction, and construct the angle atlas of the principal stress direction of the in-situ stress;

[0130] After the set of unit vectors is established, statistical analysis needs to be carried out on the direction information in the set. First, extract the direction angle represented by each unit vector. The direction angle can be expressed as the angle value rotated clockwise from the due north direction. All the direction angles are uniformly divided between 0 and 180 degrees. For example, every 10 degrees is an angle interval, and count the number of path vectors in each interval to form a direction distribution density table. Further, extract the angle intervals where some directions gather more as the key analysis objects. In each direction interval, sum up the direction angles of all the vectors in this interval, and combine the quality information of the path corresponding to this vector, such as the signal-to-noise ratio or the path length, and set it as the direction weight. Average each group of direction angles according to the weight. For example, in the interval of 30 to 40 degrees, if there are three path vectors, whose direction angles are 32 degrees, 34 degrees, and 38 degrees respectively, and the corresponding signal-to-noise ratios are 10, 12, and 8 respectively, then the weighted average value is 34.4 degrees. This value represents the dominant direction angle of the paths in this interval. After sorting out the dominant angles of all intervals, generate the angle atlas of the principal stress direction of the in-situ stress in the region, forming a direction layer with spatial distribution characteristics.

[0131] The above are only the preferred embodiments of the present invention, and do not limit the present invention in other forms. Any person skilled in the relevant art may use the technical content disclosed above to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, as long as it does not depart from the technical solution content of the present invention, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for predicting earthquakes based on the azimuth of ground stress based on the shear wave velocity change rate, characterized in that: The following steps are involved: S1: In active fault monitoring, the wave velocity value is calculated based on the arrival time of the shear wave recorded by the observation station and the distance of the propagation path. The velocity sequence is divided by a fixed-length sliding window, and the velocity difference between adjacent sliding window segments is extracted to generate a wave velocity change rate sequence set; S2: Based on the shear wave velocity change rate sequence set, extract continuous sliding window segments, calculate the velocity range in the window segment, calculate the offset rate with the median of the range as a reference value, count the window segments where the offset rate exceeds the upper limit of the range three times in a row and in the same direction, and output the velocity mutation candidate segment; S3: Based on the velocity mutation candidate section, the lithology parameters of the path in the shear wave multi-path intersection area are extracted, the lithology difference between adjacent observation points is calculated, and the lithology difference is compared with the disturbance threshold value to identify the disturbance path segment, and the velocity difference between adjacent segments is used for linear interpolation repair to generate a purified shear wave change rate sequence segment; S4: Based on the purified shear wave change rate sequence segment, the velocity change trend of the forward and reverse propagation paths is extracted, the velocity difference direction identifiers in three consecutive sliding windows are counted, the path segments with consistent forward and reverse directions are screened, and the shear wave directionality consistent segment identifier is output.

2. The method for predicting earthquake azimuth based on shear wave velocity change rate according to claim 1, characterized in that: The shear wave velocity change rate sequence set includes fixed-length sliding window division of velocity sequence, extraction of velocity differences between adjacent sliding window segments, and classification generation of shear wave velocity change rates. The velocity mutation candidate segment includes construction of maximum and minimum velocity differences between extreme ranges, determination of median reference values, calculation of sliding window offset rates, and statistics of offset rates for three consecutive segments. The purified shear wave velocity change rate sequence segment includes extraction of lithology parameters, comparison of lithology difference disturbance thresholds, identification of disturbance path segments, and linear interpolation repair of velocity differences between adjacent segments. The shear wave directionality consistent segment identification includes extraction of velocity trends of forward and reverse propagation paths, statistics of velocity difference direction identification, and screening of direction consistent path segments.

3. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 1, characterized in that: The specific steps of S1 are: S101: based on the shear wave arrival time and propagation path distance data recorded by the observation stations in the active fault monitoring area, the shear wave propagation velocity value corresponding to the propagation path is calculated, and the calculation results of all paths are integrated to generate a shear wave propagation velocity set; S102: calling the shear wave propagation velocity set, sliding along the time axis to divide continuous window segments according to a preset window length, performing mean processing on the velocity values ​​in each window segment, and generating a sliding window velocity sequence set; S103: extracting mean value data of adjacent window segments in the sliding window velocity sequence set, calculating velocity differences between window segments, grouping and classifying the differences according to propagation path numbers, and generating a shear wave velocity change rate sequence set.

4. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 3, characterized in that: The specific steps of S2 are: S201: extracting a fixed-length continuous sliding window segment based on the shear wave velocity change rate sequence set, summing up the velocity change rate differences in the segment, locating the median position after sorting by the difference values, and extracting the corresponding values, thereby generating a velocity difference median reference value; S202: performing an offset calculation on the speed change rate difference of the sliding window segment according to the speed difference median reference value to form a complete offset rate sequence, screening segments greater than the upper limit of the difference set and extracting offset directions to generate an over-threshold offset direction sequence; S203: Based on the above-threshold offset direction sequence, consistency judgment is performed on three consecutive offset directions, and sliding window segment number intervals that meet the conditions are extracted and integrated into a number range set to generate a speed mutation candidate segment.

5. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 4, characterized in that: The specific steps of S3 are: S301: based on the velocity mutation candidate section, extract the velocity data of the path in the shear wave multi-path intersection area and the stratigraphic structure information of the corresponding measuring points, classify and combine the lithology density, porosity and elastic modulus of the path section, and generate the path section lithology parameter set value; S302: comparing the lithology parameter difference between adjacent observation points with a set disturbance threshold value according to the path segment lithology parameter set value, identifying the path segment numbers exceeding the threshold value, and integrating them to generate a disturbance path segment difference sequence; S303: calling the disturbance path segment difference sequence, repairing the disturbance path segment with reference to the velocity difference of the adjacent unperturbed path segment, updating the shear wave change rate of the observation point in the path segment, and obtaining a purified shear wave change rate sequence segment.

6. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 5, characterized in that: The specific calculation formula for comparing the lithology parameter difference between adjacent observation points with the set disturbance threshold value is: ； in, represents the dynamic difference value of lithological parameters between the i-th and i-1-th observation points, represents the standardized value of the kth type of lithological parameter at the ith observation point, represents the standardized value of the kth type of lithological parameter at the i-1th observation point, represents the weight coefficient of the jth environmental factor, represents the normalized impact factor of the jth type of environmental factor at the ith observation point, represents the spatial distance between the i-th and i-1 observation points, represents the path segment curvature correction coefficient, Represents the distance decay exponent.

7. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 5, characterized in that: The specific steps of S4 are: S401: extracting the velocity change direction of the shear wave propagation path in the time window based on the purified shear wave change rate sequence segment, identifying the change of the path direction in the adjacent window, and generating a velocity change direction sequence; S402: calling the speed change direction sequence, extracting the path direction sequence by using a three-segment sliding window method, identifying whether the directions in each group of windows are consistent, recording the path segment information that meets the conditions, and obtaining a set of continuous sliding window direction consistent interval labels; S403: According to the continuous sliding window direction consistent interval label set, the shear wave propagation paths with consistent directions are screened, the path segments with constant directions are identified, and the shear wave directionality consistent segment identifiers are generated.

8. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 7, characterized in that: The specific calculation formula for the shear wave propagation path with the screening direction being consistent is: ； in, represents the mean of all shear wave propagation direction angles in the i-th sliding window, n represents the total number of windows in the set of consecutive sliding window direction consistent interval labels, represents the balance factor based on the dynamic distribution of the angular amplitude of the shear wave propagation direction, and i represents the index number of the sliding window.

9. The method for predicting earthquakes based on the in-situ stress azimuth based on the shear wave velocity change rate according to claim 1, characterized in that: The method further comprises: S5: Based on the segment identification with consistent shear wave directivity, extract the angle between the segment velocity change direction and the structural strike, construct the path unit vector through angle correction, integrate the direction vectors of all paths, and form a principal stress direction map of ground stress; The in-situ stress principal stress direction map includes calculation of the angle between the velocity change direction and the structural strike, construction of path unit vector correction, and integration of direction vectors.

10. The method for predicting earthquakes based on the in-situ stress orientation based on the shear wave velocity change rate according to claim 9, characterized in that: The specific steps of S5 are: S501: based on the segment identifier of the consistent shear wave directivity, extract the velocity change direction of the segment shear wave path, combine the strike azimuth information of the structural belt, analyze the angle relationship between the path direction and the structural strike, and generate the angle value between the path and the structural strike; S502: According to the angle between the path and the structural strike, the direction vector of the path is corrected, the numerical structure of the direction vector is unified, and all path results are integrated to obtain path unit vector set data; S503: Based on the path unit vector set data, extract the spatial distribution trend of the direction vector, classify and count the directions, and calculate the angle mean to construct a geostress principal stress direction angle map.

Images