Stress orientation prediction method based on shear wave velocity change rate

Abstract

The present application relates to the technical field of earthquake prediction, in particular to a method for predicting the direction of ground stress based on the rate of change of shear wave velocity, comprising the following steps: calculating the velocity value of shear wave and propagation distance, extracting the difference value to generate a change rate set, calculating the offset rate and screening the continuous over-limit section, identifying the disturbance path section for interpolation repair to generate a purified sequence, extracting the forward and backward propagation trend and counting the direction identifier, and generating a ground stress principal stress direction atlas; the present application generates a change rate set by dividing the velocity sequence with a sliding window, enhances the data continuity and retains the local features, combines the continuous offset rate to screen and accurately locate the mutation area, identifies the multi-path lithology difference and disturbance repair, excludes non-structural interference, cross-verification of the forward and backward paths screens the direction consistent section, improves the discrimination robustness, dynamically correlates the velocity direction and the structural trend, realizes the synchronous analysis of the stress field and fault activity, combines the dynamic threshold and direction verification in multiple stages, and improves the derivation accuracy and timeliness.

Description

Technical Field

[0001] This invention relates to the field of earthquake prediction technology, and in particular to a method for predicting the azimuth of ground stress based on the rate of change of shear wave velocity. Background Technology

[0002] Earthquake prediction technology encompasses a system of methods for predicting the time, location, and intensity of earthquakes, and is a key research direction at the intersection of seismology and geology. The core content of this technology includes the monitoring and identification of earthquake precursors, statistical analysis of seismic activity, monitoring of changes in subsurface physical parameters, and evolutionary analysis of tectonic stress fields. This field typically involves the analysis of seismic wave propagation characteristics, the inversion and derivation of stress field distribution, and the dynamic monitoring of subsurface fault activity. It also extensively integrates geophysical observation data, microseismic activity records, and information on the changing trends of elastic parameters in the subsurface medium to construct multi-factor correlation models to determine the possible earthquake occurrence areas and potential risk levels.

[0003] The seismic prediction method based on the rate of change of shear wave velocity refers to the method of inferring the direction of principal stress within the underground rock mass by utilizing the rate of change of the propagation velocity of underground shear waves with time or space, and using this as an important basis for judging seismic activity. The main technical issues addressed in this patent include the identification of the correspondence between the characteristics of shear wave velocity variation and the geostress field, the extraction method of the principal stress direction, and the analysis of its spatial distribution evolution trend. Specifically, by acquiring regional seismic wave data, comparing the rate of change of shear wave velocity along different propagation paths and time periods, and calculating and deriving the principal stress direction based on the functional relationship between stress and wave velocity in elasticity theory, a set of geostress orientation discrimination procedures for seismic activity prediction is formed.

[0004] Existing technologies rely on static inversion models that relate seismic wave propagation characteristics to stress fields, failing to incorporate dynamically segmented data using sliding time windows. This leads to local stress variation characteristics being masked by the overall trend, affecting the timeliness of stress direction identification. When screening velocity anomaly segments based on fixed thresholds, the lack of directional consistency verification of continuous sliding window offset rates makes them susceptible to random noise interference, resulting in misjudgments and reduced accuracy in candidate segment location. Existing methods lack multi-path lithological parameter difference comparison and disturbance path segment repair mechanisms, failing to distinguish between data fluctuations caused by tectonic and non-tectonic factors, leading to systematic biases in stress field inversion results. Traditional directional analysis is limited to single propagation path trend discrimination, failing to introduce cross-validation of forward and reverse path direction markers, making it difficult to exclude the influence of multi-path wave velocity superposition effects and reducing the reliability of principal stress direction discrimination. The lack of correlation between static tectonic strike and dynamic velocity direction leads to a risk of decoupling between stress field evolution trends and actual fault activity. For example, in strike-slip fault areas, the principal stress direction offset may be misjudged, affecting the accuracy of seismic risk level assessment. Existing technologies rely on analysis of single geophysical parameters and fail to construct multi-factor dynamic correlation models, which limits the refined analysis of the spatial distribution characteristics of geostress fields in complex tectonic environments. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity, comprising the following steps:

[0007] S1: In active fault monitoring, wave velocity values ​​are calculated based on the arrival time of shear waves and the distance of the propagation path recorded by the observation stations. A fixed-length sliding time window is used to divide the velocity sequence, and the velocity difference between adjacent sliding window segments is extracted to generate a set of wave velocity change rate sequences.

[0008] S2: Based on the shear wave velocity change rate sequence set, extract continuous sliding window segments, calculate the velocity range within the window segment, calculate the offset rate with the median of the range as a reference value, count the window segments in which the offset rate exceeds the upper limit of the range three times in a row and the direction is consistent, and output the candidate segments of velocity change abruptness.

[0009] S3: Based on the candidate segments of velocity mutation, extract the lithological parameters of the paths within the cross region of the shear wave multipath, calculate the lithological difference between adjacent observation points, compare it with the disturbance threshold value, identify the disturbance path segment, use the velocity difference between adjacent segments for linear interpolation repair, and generate the purified shear wave rate of change sequence segment.

[0010] S4: Based on the purified shear wave rate of change sequence segment, extract the velocity change trend of the forward and reverse propagation paths, statistically analyze the velocity difference direction indicators within three consecutive sliding window segments, filter out path segments with consistent forward and reverse directions, and output the shear wave directionality consistent segment indicators.

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

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

[0013] S101: Based on the arrival time and propagation path distance data of shear waves recorded by observation stations in 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 set of shear wave propagation velocities.

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

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

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

[0017] S201: Based on the set of shear wave velocity change rate sequences, extract continuous sliding window segments of fixed length, summarize the velocity change rate differences within the segments, sort them by the size of the differences, locate the median position, extract the corresponding values, and generate the median reference value of the velocity difference.

[0018] S202: Based on 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, filter out segments that are greater than the upper limit of the difference set and extract the offset direction to generate an over-threshold offset direction sequence.

[0019] S203: Based on the above-threshold offset direction sequence, perform consistency judgment on three consecutive offset directions, extract the sliding window segment number intervals that meet the conditions and integrate them into a number range set to generate candidate segments for velocity mutation.

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

[0021] S301: Based on the candidate velocity change zone, extract the velocity data of the path and the stratigraphic structure information of the corresponding measuring point in the cross region of the shear wave multipath, classify and combine the lithological density, porosity and elastic modulus of the path segment, and generate the lithological parameter set value of the path segment.

[0022] S302: Based on the set of lithological parameters of the path segment, compare the difference in lithological parameters between adjacent observation points with the set disturbance threshold value, identify the path segment number that exceeds the threshold, and generate a disturbance path segment difference sequence after integration;

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

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

[0025] ;

[0026] in, This represents the dynamic difference in lithological parameters between the i-th and i-1th observation points. This represents the standardized value of the k-th type of lithological parameter at the i-th observation point. This represents the standardized value of the k-th type of lithological parameter at the (i-1)-th observation point. The weight coefficient representing the j-th type of environmental factor. The normalized influence factor representing the j-th type of environmental factor at the i-th observation point. This represents the spatial distance between the i-th and i-1th observation points. Represents the curvature correction factor for the path segment. This represents the distance decay index.

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

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

[0029] S402: Call the velocity 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 path segment information that meets the conditions, and obtain the label set of consecutive sliding window direction consistent intervals;

[0030] S403: Based on the set of labels for the continuous sliding window direction-consistent intervals, filter out transverse wave propagation paths with consistent directions, identify path segments with constant directions, and generate transverse wave direction-consistent interval identifiers.

[0031] As a further aspect of the present invention, the specific calculation formula for the transverse wave propagation path with consistent screening direction is as follows:

[0032] ;

[0033] in, represents the mean angle of all transverse wave propagation directions within the i-th sliding window, and n represents the total number of windows in the set of consecutive sliding window intervals with consistent directions. Represents the balance factor based on the dynamic distribution of the angular amplitude of the transverse wave propagation direction, and i represents the index number of the sliding window.

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

[0035] S5: Based on the segment identifier with consistent transverse wave directionality, 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 the geostress principal stress direction map.

[0036] The geostress principal stress direction map includes the calculation of the angle between the velocity change direction and the structural trend, the construction of path unit vector correction, and the integration of direction vectors.

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

[0038] S501: Based on the shear wave directionality consistent section identifier, extract the velocity change direction of the shear wave path within the section, combine it with the orientation information of the structural zone, analyze the angle relationship between the path direction and the structural orientation, and generate the angle value between the path and the structural orientation.

[0039] S502: Based on the angle between the path and the construction direction, the path direction vector is corrected, the numerical structure of the direction vector is unified, and all path results are integrated 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, classify and statistically analyze them according to direction, calculate the mean angle, and construct the geostress principal stress direction angle map.

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

[0042] In this invention, a set of rate of change is generated by segmenting the velocity sequence through a sliding time window, which enhances data continuity and preserves local features. Combined with continuous offset rate screening, abrupt change zones are accurately located. Multi-path lithological difference identification and disturbance repair are performed to eliminate non-tectonic interference. Cross-validation of forward and reverse paths is used to screen sections with consistent directions, which improves the robustness of the judgment. The velocity direction is dynamically correlated with the structural strike, which enables synchronous analysis of stress field and fault activity. Multi-stage processing combined with dynamic thresholds and direction verification improves the accuracy and timeliness of the derivation. Attached Figure Description

[0043] Figure 1 This is a schematic diagram of the steps of the present invention. Detailed Implementation

[0044] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0045] In the description of this invention, it should be understood that the terms "length," "width," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicating orientation or positional relationships, are based on the orientation or positional relationships shown in the accompanying drawings and are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, in the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0046] Please see Figure 1 A seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity includes the following steps:

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

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

[0049] S3: Based on the candidate segments of velocity mutation, extract the lithological parameters of the paths in the cross region of the shear wave multipath, calculate the lithological difference between adjacent observation points, compare it with the disturbance threshold value, identify the disturbance path segment, use the velocity difference between adjacent segments for linear interpolation repair, and generate the purified shear wave rate of change sequence segment.

[0050] S4: Based on the purified shear wave rate of change sequence segment, extract the velocity change trend of the forward and reverse propagation paths, statistically analyze the velocity difference direction indicators within three consecutive sliding window segments, filter the path segments with consistent forward and reverse directions, and output the shear wave directionality consistent segment indicators.

[0051] S5: Based on the segment identifier with consistent shear wave directionality, 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 the geostress principal stress direction map.

[0052] The shear wave velocity change rate sequence set includes fixed-length sliding window division of the velocity sequence, extraction of velocity difference between adjacent sliding window segments, and generation of shear wave velocity change rate classification. The velocity abrupt change candidate segment includes the construction of the maximum and minimum velocity difference in the range interval, determination of the median reference value, calculation of sliding window offset rate, and statistical analysis of offset rates for three consecutive segments. The purified shear wave velocity change rate sequence segment includes lithological parameter extraction, comparison of lithological difference disturbance thresholds, identification of disturbance path segments, and linear interpolation repair of velocity difference between adjacent segments. The shear wave directionality consistency segment identification includes extraction of velocity trends in forward and reverse propagation paths, statistical analysis of velocity difference direction identification, and screening of directionally consistent path segments. The geostress 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.

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

[0054] S101: Based on the arrival time and propagation path distance data of shear waves recorded by observation stations in 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 set of shear wave propagation velocities.

[0055] In active fault monitoring areas, after the observation stations are deployed, their specific spatial coordinates need to be collected and labeled. For example, station A can be located at a certain longitude and latitude, and station B can be located in a similar longitude and latitude region, forming a station network covering the target monitoring area. In actual earthquake events, the arrival time of the shear wave corresponding to a specific station is extracted using seismic wave recording equipment. For example, in a certain event, station A recorded the shear wave arrival at 13.45 seconds. At the same time, based on the source location information and 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 velocity 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 velocity of the path. Further data extraction is performed by combining other earthquake events with other stations. For example, station B's measured distance from the epicenter is 15.3 kilometers, and the arrival time is 18.6 seconds. After processing, the velocity corresponding to this path is obtained. The arrival time and propagation path distance collected from each event are uniformly converted into velocity values ​​according to the above method. After being summarized, a set of shear wave propagation velocities is formed. This set contains velocity values ​​corresponding to multiple paths and multiple event combinations, which constitute the basic sample set of shear wave propagation velocities in the monitoring area.

[0056] S102: Call the transverse wave propagation velocity set, slide along the time axis to divide the continuous window segment according to the preset window length, perform average processing on the velocity value in each window segment, and generate a sliding window velocity sequence set;

[0057] After calling the aforementioned set of transverse wave propagation velocities, a fixed-length sliding window segment needs to be set on the time axis. The window segmentation is performed by setting a sliding step size; for example, if each window segment is 60 seconds long and the sliding step size is 30 seconds, multiple overlapping window segments are generated by advancing forward along the time axis. Each window segment contains several velocity records. If a window segment contains 10 velocity values, these 10 values ​​need to be averaged to obtain the average propagation velocity within that time period. This method is used to process all sliding window segments one by one, generating a complete window average velocity sequence. During the operation, for some window segments containing too little velocity data, a minimum sample size limit is set; if there are fewer than 3 records, the data in that segment is discarded to ensure the representativeness of the mean. The window segmentation strategy can adopt a sliding ratio setting method, for example, a sliding step size of half the window length to ensure coverage continuity. After the above steps, a complete set of sliding window velocity sequences is formed, providing a stable data sequence foundation for subsequent dynamic change analysis.

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

[0059] After extracting the mean of two adjacent window segments in the sliding window velocity sequence, the velocity difference between them is calculated. For example, if the mean of the previous window is 0.803 and the mean of the next window is 0.790, the velocity difference between the two segments is negative, and the absolute difference is used as the basis for analysis. These differences are categorized and statistically analyzed according to the propagation path number. For example, if there are 10 sets of velocity change records under the Path-001 path, they are summarized into that path category. Subsequently, the change sequence under this path is processed for continuity. A three-point moving average can be used to process each difference in the sequence to eliminate the influence of short-term abnormal fluctuations and improve data stability. When classifying whether a change is significant, a threshold limit for significant change is set. For example, if the velocity change rate of a large number of paths in historical statistics is concentrated below 0.018, 0.020 can be set as the boundary point for distinguishing whether it is significant. Any path whose change value exceeds this limit within a certain period of time is considered to have a strong fluctuation phenomenon in that path during that period. Finally, the velocity differences between each segment in all paths are summarized to form a set of velocity change rate sequences organized by path number, which serves as the basis for analyzing the dynamic characteristics of transverse wave propagation.

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

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

[0062] To extract equal-length sliding window segments from a continuous data sequence of shear wave velocity change rates, the window length must first be determined. For example, five observation points can be selected to form a segment, and the sliding window can be operated with a step size of one observation point. These sliding window segments are continuously extracted from the entire velocity change rate sequence. Within each segment, the recorded shear wave velocity change rate values ​​are arranged in chronological order, and the difference between the maximum and minimum values ​​is calculated for each segment to obtain the velocity change rate difference. For example, if the data in a segment are 0.03, 0.02, 0.04, 0.01, and 0.05, the difference is the maximum value of 0.05 minus the minimum value of 0.01, which equals 0.04. After collecting and sorting the differences from all sliding window segments, from smallest to largest, the median value is selected as the baseline. If the total number of data points is even, the average of the two median values ​​is taken. For example, if the summary results are 0.01, 0.02, 0.03, 0.04, 0.05, and 0.06, the median is between the third and fourth values, and the baseline value is 0.035. This value serves as a reference standard for subsequent comparisons and evaluations.

[0063] S202: Based on the median reference value of the velocity difference, the offset rate of the velocity change rate difference of the sliding window segment is calculated to form a complete offset rate sequence. Segments that are greater than the upper limit of the difference set are selected and the offset direction is extracted to generate an over-threshold offset direction sequence.

[0064] After determining the baseline value, each sliding window segment is iterated again, and the offset rate is calculated based on the difference between its rate of change and the baseline value. The offset rate is defined as the proportion of the increase or decrease of the current segment's difference relative to the baseline value. If a segment's difference is 0.05 and the baseline value is 0.035, its offset rate is expressed as the proportion of the segment's deviation to the baseline value, meaning the difference is approximately 42% higher than the baseline value. After calculating all sliding window segments sequentially, a complete offset rate sequence is formed. To identify significant abnormal changes, a judgment threshold needs to be set. This value can be set based on the distribution characteristics of the difference set, such as using a combination of the higher quartile and the interquartile range. Segments exceeding this upper limit are considered abnormal, and their offset direction is recorded. If the difference is higher than the baseline value, the direction is marked as positive; otherwise, it is marked as 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 over-threshold offset direction sequence, the consistency of three consecutive offset directions is judged, the sliding window segment numbering intervals that meet the conditions are extracted and integrated into a numbering range set to generate candidate segments for velocity mutation.

[0066] In the acquired offset direction sequence, each segment is checked for continuous and consistent directions. If three or more consecutive segments have positive or negative direction markers, the sequence can be considered to have a stable and consistent trend. At this point, the sliding window segment numbers corresponding to these continuous and consistent directions are extracted and merged 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 considered as two segments with continuous offset directions. During the judgment, it is necessary to ensure that the length of the direction markers reaches at least the set minimum number of segments to avoid judgment bias caused by occasional errors. After mapping the number range to actual time, candidate velocity mutation regions are generated, facilitating subsequent spatial localization or temporal series studies. This process, based on the sliding window segment numbers, sequentially traverses and integrates all sequence intervals that meet the consistency condition to form the final set of candidate velocity mutation regions.

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

[0068] S301: Based on the candidate segments of velocity abrupt change, extract the velocity data of the path and the stratigraphic structure information of the corresponding measuring points in the cross region of the shear wave multipath, classify and combine the lithological density, porosity and elastic modulus of the path segments, and generate the set of lithological parameter values ​​of the path segments.

[0069] After identifying path regions with significant velocity abrupt changes, the process first filters out sections where the velocity change rate exceeds a preset threshold based on the spatial variation of shear wave propagation velocity in the seismic data. For example, a shear wave velocity change rate higher than 15% can be considered a velocity abrupt change segment. Subsequently, data from each observation point covered by the shear wave propagation path is extracted from these segments, and the corresponding stratigraphic structure information is obtained from existing geological profiles or drilling data. This structural information typically includes stratum name, stratum thickness, burial depth, and other physical properties. Next, using typical lithological parameters recorded in the geological database, the lithology of the strata at each observation point is classified, and basic parameters such as density, porosity, and elastic modulus are extracted based on the classification results. For example, a path may contain clay and sand strata; these are classified into their respective categories through lithological comparison, and their physical parameters are obtained from tables. If the lithology of an observation point is silty clay, its typical density is 1.85 g / cm³, porosity is 42%, and elastic modulus is 150 MPa; while sand layers may have a density of 2.10 g / cm³, porosity of 35%, and elastic modulus of 320 MPa. In constructing the lithological parameter set for the path segment, parameters for missing or anomalous observation points need to be supplemented. This is typically done using linear interpolation, derived by averaging parameters from adjacent observation points, ensuring the integrity of the parameter set and the logical consistency between data. The final parameter set will serve as the basis for subsequent difference identification and disturbance analysis.

[0070] S302: Based on the set values ​​of lithological parameters of the path segment, compare the difference of lithological parameters between adjacent observation points with the set disturbance threshold value, identify the path segment number that exceeds the threshold, and generate a disturbance path segment difference sequence after integration;

[0071] The specific calculation formula for comparing the difference in lithological parameters between adjacent observation points with the set disturbance threshold is as follows:

[0072] ;

[0073] in, This represents the dynamic difference in lithological parameters between the i-th and i-1th observation points. This represents the standardized value of the k-th type of lithological parameter at the i-th observation point. This represents the standardized value of the k-th type of lithological parameter at the (i-1)-th observation point. The weight coefficient representing the j-th type of environmental factor. The normalized influence factor representing the j-th type of environmental factor at the i-th observation point. This represents the spatial distance between the i-th and i-1th observation points. Represents the curvature correction factor for the path segment. Represents the distance decay index;

[0074] Standardized values ​​of lithological parameters: and The standardized values ​​of the k-th type of lithological parameters at the i-th and i-1 observation points are obtained by normalizing the actual monitoring data.

[0075] Example: The original resistivity of sandstone at the i-th observation point in a certain area is 150 Ω·m. After standardization (formula...) ,in Ω˙m, Ω˙m) The resistivity of the sandstone at the (i-1)th observation point is 120 Ω·m, which, after standardization, yields... .

[0076] Environmental factor weighting coefficients: The weight coefficients for the j-th type of environmental factors are determined based on the assessment of geological experts and the results of principal component analysis.

[0077] Example: Regional environmental factors include rock moisture (j=1), slope (j=2), and fracture density (j=3), with weights set as follows: , , The basis for this is that humidity accounts for 40% of the impact on lithological stability, while slope and fissure each account for 30%.

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

[0079] Example: A humidity sensor measures the humidity at point i to be 25%. Normalize ( The slope measurement was 15°, normalized ( The fracture density is 8 fractures / m², normalized ( ).

[0080] Spatial distance: Let be the spatial distance between the i-th and i-1-th observation points, calculated using GPS coordinates in Euclidean form.

[0081] Example: If the coordinate difference between two points is Δx = 30m and Δy = 40m, then m.

[0082] Curvature correction factor: Based on historical path curvature data fitting, (for every 0.1m increase in curvature)... -1 , Increase by 0.05, in the range [0.5, 1.5]); Distance decay exponent: The range is [0.3, 0.7].

[0083] Formula calculation derivation:

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

[0085] ;

[0086] Calculate the weighted sum of environmental factors:

[0087] ;

[0088] Calculate the numerator:

[0089] ;

[0090] Calculate the denominator:

[0091] ;

[0092] Calculate the dynamic difference value:

[0093] ;

[0094] Parameter summary: Represents the dynamic differences in lithological parameters. and To standardize lithological parameters, The weights of environmental factors, These are normalized measured values ​​of environmental factors. For spatial distance, This is the curvature correction factor. The distance decay exponent;

[0095] Analysis of the results of the example: This indicates that the difference in lithological parameters between observation points i and i-1 is in a low-disturbance state after environmental and spatial correction. If the disturbance threshold is set to 0.05, this path segment is not marked as a disturbance segment and needs to be excluded when integrating it into the sequence.

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

[0097] After obtaining the difference sequence of the disturbed path segments, the actual location of each disturbed segment must first be determined, and its adjacent undisturbed path segments must be found to obtain a reference value. For the reference path segment, the change in its shear wave velocity between two adjacent observation points can be calculated as a velocity change benchmark to repair 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 based on the velocity change trend of the reference segment. The estimation method can be to derive the velocity value from the adjacent undisturbed point according to its rate of change and the distance between the point and the disturbed point. For example, if the velocities of two points in the left reference path segment are 850 m / s and 870 m / s respectively, and the distance between them is 20 meters, and the derived rate of change is 1 m / s per meter, then if the distance between the disturbed point and the left reference point is 15 meters, the velocity of that point can be estimated to be 875 m / s. Similarly, the velocity repair is performed on all observation points within the entire disturbed segment to form a new velocity sequence. The velocity changes after repair can be analyzed again using difference analysis to ensure that the overall trend is smooth and without abrupt changes. This process is suitable for data cleaning and noise removal scenarios, and its effectiveness is particularly pronounced in seismic data segments with discontinuities, strong interference, or localized anomalous signals. The cleaned shear wave rate of change sequence will be used in 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 rate of change sequence segment, extract the velocity change direction of the shear wave propagation path in the time window, identify the change of the path direction in adjacent windows, and generate a velocity change direction sequence.

[0100] To obtain information on the rate of change of seismic shear waves propagating in a sensor array, specific frequency bands, such as fluctuation data in the range of 5 to 20 Hz, need to be extracted from the original seismic waveform. The original waveform is first processed by denoising and filtering to eliminate environmental interference and measurement errors. Next, the propagation delay rate of change is calculated based on the waveform-time difference and distance between multiple sensors. The processed rate of change sequence needs to be standardized to ensure comparability across multiple paths. To investigate the velocity change trend along the path, a time window of a certain length can be set, and the rate of change sequences within the window can be aggregated and their direction of change along the spatial path analyzed. The propagation velocity change of each path within this time period can be converted into a direction vector, representing the velocity change trend of the path. When the time window slides forward to an adjacent time period, the angle between the direction vectors of the same path in adjacent windows needs to be compared. If the angle is greater than a specific angle threshold, such as 30 degrees, it is determined that the direction has changed; otherwise, the direction is considered to remain consistent. The direction changes of all paths within continuous time periods are recorded one by one, ultimately forming a complete velocity direction change sequence, providing a basis for subsequent analysis.

[0101] S402: Call the velocity change direction sequence, extract the path direction sequence using a three-segment sliding window method, identify whether the directions within each window are consistent, record the path segment information that meets 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, each 10 seconds long with a 5-second sliding interval. For each path, its directional features are extracted within each time period, and the path directions in these three time periods are compared. If the angle difference between any two path directions in the three periods does not exceed a set range (e.g., 20 degrees), the path segment is recorded as having the same direction. As the sliding window moves forward successively, all path segments that meet the directional consistency condition are recorded. Each record includes the path number, directional features, and time interval. Then, overlapping or adjacent records are merged based on continuity on the time axis to form multiple sets of time intervals with consistent direction. The corresponding path segment identifiers are summarized within each set to obtain a set of time segment labels for path direction consistency, providing support for path filtering.

[0103] S403: Based on the set of labels for consecutive sliding window directions with consistent orientation, filter out transverse wave propagation paths with consistent orientation, identify path segments with constant orientation, and generate transverse wave directionality consistent segment identifiers.

[0104] The specific calculation formula for selecting transverse wave propagation paths with consistent direction is as follows:

[0105] ;

[0106] in, represents the mean (in radians) of the angles of all transverse wave propagation directions within the i-th sliding window, and n represents the total number of windows in the set of consecutive sliding window intervals with consistent directions. The balance factor (defined as) represents the dynamic distribution of the amplitude of the angle based on the propagation direction of the transverse wave. ), i represents the index number of the sliding window (from 1 to n);

[0107] Parameter definition and data source:

[0108] The average angle (in radians) of the propagation direction of shear waves within a sliding window is recorded using seismic monitoring equipment. The monitoring equipment records waveforms at a sampling frequency of 1000 times per second and calculates the arithmetic mean of the angles of all shear wave directions within the window. For example, the first window measures... radians, measured in the second window Radius, and so on.

[0109] n: The total number of windows in the continuous sliding window interval label set with consistent direction. Calculated based on the original seismic signal length (e.g., duration 120 seconds) and the sliding window length (e.g., 10 seconds per window, step size 5 seconds). .

[0110] The balance factor is defined as follows: Calculated based on monitoring data ,For example radians, substituting them gives ,because Ultimately .

[0111] Example calculation process:

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

[0113] calculate radian;

[0114] calculate ;

[0115] Calculate the numerator ;

[0116] Calculate the first term ;

[0117] Calculate the second term ;

[0118] Merge results ;

[0119] Parameter setting basis and numerical rationality:

[0120] threshold (Approximately 31.4 radians) Based on statistical experience regarding the range of angular fluctuations in the direction of seismic shear waves (in most scenarios) (radians), ensure The dynamic adjustment range is reasonable.

[0121] The sliding window length of 10 seconds and the step size of 5 seconds are based on seismic signal processing standards (such as the USGS recommended method) to balance time resolution and computational efficiency.

[0122] Relationship between results and steps:

[0123] In the example, H=0.7947. When H≥0.75, the directional consistency index is considered met, and a corresponding transverse wave directional consistency segment identifier is generated. This result indicates that the transverse wave propagation direction fluctuation within the current window sequence is small, meeting the requirement of directional constancy.

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

[0125] S501: Based on the identification of the shear wave directionality consistent section, extract the velocity change direction of the shear wave path within the section, combine the orientation information of the tectonic zone, analyze the angle relationship between the path direction and the tectonic orientation, and generate the angle value between the path and the tectonic orientation.

[0126] After identifying sections with consistent shear wave directionality, it is necessary to perform directional consistency analysis on the shear wave signal data acquired along the seismic line. This typically involves using shear wave waveforms received from multiple seismic points. By comparing the main directional changes of waveforms between adjacent points, sections with good continuity and consistent directionality are identified. This type of analysis can be performed using seismic data processing platforms such as SeisSpace or Petrel. Within the identified sections, multiple shear wave propagation paths are selected. Using the first arrival times and path geometric lengths, the velocity of the shear wave on each path is calculated. Then, based on the direction of these velocity values ​​changing between different paths, the propagation trend of the shear wave within the section is determined. If the velocity on a certain path... A larger angle indicates a greater likelihood that the transverse wave will preferentially choose its propagation direction. For example, if the transverse wave propagation speed between measuring points A and B is 2900 m / s and between A and C is 3100 m / s, then the transverse wave is considered to tend to propagate along the direction from A to C. Next, the orientation information of the tectonic zone is introduced. Based on the known azimuth of the tectonic zone, the path direction is compared with the tectonic orientation, and the angle between them is measured. The angle value is calculated through the spatial relationship between the path direction coordinates and the tectonic orientation. For example, if the path direction is northeast and the tectonic orientation is due east, then the angle between the two is 45 degrees. The angle value between the path and the tectonic orientation is generated in this way for subsequent direction adjustment.

[0127] S502: Based on the angle between the path and the construction direction, the path direction vector is corrected, the numerical structure of the direction vector is unified, and all path results are integrated to obtain the path unit vector set data.

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

[0129] S503: Based on the path unit vector set data, extract the spatial distribution trend of the direction vectors, classify and statistically analyze them according to the direction, calculate the mean angle, and construct the geostress principal stress direction angle map.

[0130] After establishing the unit vector set, the directional information in the set needs to be statistically analyzed. First, the directional angle represented by each unit vector is extracted. The directional angle can be represented as the angle value rotating clockwise from due north. All directional angles are uniformly divided between 0 and 180 degrees, for example, each 10-degree interval. The number of path vectors in each interval is counted to form a directional distribution density table. Further, angle intervals with a high concentration of certain directions are extracted as key analysis objects. Within each directional interval, the directional angles of all vectors in that interval are summarized, and combined with the quality information of the path corresponding to that vector, such as signal-to-noise ratio or path length, are set as directional weights. Each group of directional angles is averaged according to the weights. For example, in the 30 to 40 degree interval, if there are three path vectors with directional angles of 32 degrees, 34 degrees, and 38 degrees, and corresponding signal-to-noise ratios of 10, 12, and 8, the weighted average is 34.4 degrees. This value represents the dominant directional angle of the path in that interval. After sorting out the dominant angles of all intervals, an angle map of the principal stress direction in the region is generated, forming a directional layer with spatial distribution characteristics.

[0131] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for earthquake prediction based on the rate of change of shear wave velocity, characterized by, Includes the following steps: S1: In active fault monitoring, wave velocity values ​​are calculated based on the arrival time and propagation path distance of shear waves recorded by observation stations. A fixed-length sliding time window is used to divide the velocity sequence, and the velocity difference between adjacent sliding time window segments is extracted to generate a set of shear wave velocity change rate sequences. S2: Based on the shear wave velocity change rate sequence set, extract continuous sliding window segments, calculate the velocity range within the window segment, calculate the offset rate with the median of the range as a reference value, count the window segments in which the offset rate exceeds the threshold three times in a row and the offset direction is consistent, and output the candidate segments of velocity abrupt change. S3: Based on the candidate velocity mutation zone, extract the lithological parameters of the path within the cross region of the shear wave multipath, calculate the lithological difference between adjacent observation points, compare it with the disturbance threshold value, identify the disturbed path segment, use the velocity difference of adjacent undisturbed path segments for linear interpolation repair, and generate the purified shear wave rate of change sequence segment. S4: Based on the purified shear wave rate of change sequence segment, extract the velocity change trend of the forward and reverse propagation paths, statistically analyze the velocity difference direction indicators within three consecutive sliding windows, filter the path segments with consistent forward and reverse directions, and output the shear wave directionality consistent section indicators. The specific steps are as follows: S401: Based on the purified transverse wave rate of change sequence segment, extract the velocity change direction of the transverse wave propagation path in the time window, identify the change of the path direction in adjacent windows, and generate a velocity change direction sequence. S402: Call the velocity 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 path segment information that meets the conditions, and obtain the set of labels for consecutive sliding window directions that are consistent. S403: Based on the set of labels for the continuous sliding window with consistent direction, filter the transverse wave propagation paths with consistent direction, identify the path segments with constant direction, and generate transverse wave directionality consistent segment identifiers. The method also include, S5: Based on the shear wave directionality consistent section identifier, extract the angle between the section velocity change direction and the structural strike, construct the path unit vector through angle correction, integrate the direction vectors of all paths, and form the geostress principal stress direction angle map. The specific steps are as follows: S501: Based on the identification of the shear wave directionality consistent section, extract the velocity change direction of the shear wave path of the section, combine it with the orientation information of the structural zone, analyze the angle relationship between the path direction and the structural orientation, and generate the angle value between the path and the structural orientation. S502: Based on the angle between the path and the construction direction, the path direction vector is corrected, the numerical structure of the direction vector is unified, and all path results are integrated to obtain the path unit vector set data. S503: Based on the path unit vector set data, extract the spatial distribution trend of the direction vectors, classify and statistically analyze them according to direction, calculate the mean angle, and construct the geostress principal stress direction angle map.

2. The seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Based on the arrival time and propagation path distance data of shear waves recorded by observation stations in 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 set of shear wave propagation velocities. S102: Call the shear wave propagation velocity set, slide along the time axis to divide it into continuous window segments according to the preset window length, perform average processing on the velocity values ​​in each window segment, and generate a sliding window velocity sequence set; S103: Extract the mean data of adjacent window segments from the sliding window velocity sequence set, calculate the velocity difference between window segments, group and classify the differences according to the propagation path number, and generate a transverse wave velocity change rate sequence set.

3. The seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity according to claim 2, characterized in that, The specific steps of S2 are as follows: S201: Based on the set of shear wave velocity change rate sequences, extract continuous sliding window segments of fixed length, summarize the velocity change rate differences within the segments, sort them by the size of the differences, locate the median position, extract the corresponding values, and generate a median reference value for the velocity difference. S202: Based on the median reference value of the velocity difference, calculate the offset rate of the velocity change rate difference of the sliding window segment to form a complete offset rate sequence, filter out segments that are greater than the threshold of the difference set and extract the offset direction to generate a sequence of offset directions exceeding the threshold. S203: Based on the above-threshold offset direction sequence, perform consistency judgment on three consecutive offset directions, extract the sliding window segment number intervals that meet the conditions and integrate them into a number range set to generate candidate segments for velocity mutation.

4. The seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity according to claim 3, characterized in that, The specific steps of S3 are as follows: S301: Based on the candidate velocity change zone, extract the velocity data of the path and the stratigraphic structure information of the corresponding measuring point in the cross region of the shear wave multipath, classify and combine the lithological density, porosity and elastic modulus of the path segment, and generate the lithological parameter set value of the path segment. S302: Based on the set of lithological parameters of the path segment, compare the difference in lithological parameters between adjacent observation points with the set disturbance threshold value, identify the path segment number that exceeds the threshold, and generate a disturbance path segment difference sequence after integration; S303: Call the disturbance path segment difference sequence, use the velocity difference between adjacent undisturbed path segments as a reference to repair the disturbance path segment, update the shear wave change rate of the observation points in the path segment, and obtain the purified shear wave change rate sequence segment.

5. The seismic prediction method for geostress azimuth based on the rate of change of shear wave velocity according to claim 4, characterized in that, The specific calculation formula for comparing the difference in lithological parameters between adjacent observation points with the set disturbance threshold value is as follows: ； in, This represents the dynamic difference in lithological parameters between the i-th and i-1th observation points. This represents the standardized value of the k-th type of lithological parameter at the i-th observation point. This represents the standardized value of the k-th type of lithological parameter at the (i-1)-th observation point. The weight coefficient representing the j-th type of environmental factor. The normalized influence factor representing the j-th type of environmental factor at the i-th observation point. This represents the spatial distance between the i-th and i-1th observation points. Represents the curvature correction factor for the path segment. This represents the distance decay index.

Images