Methods for Differentiating and Correlating Seismic Ground Motion Parameters of Upper and Lower Reservoirs in Large Pumped Storage Power Stations

Abstract

This invention relates to the field of seismic wave analysis technology, specifically a method for distinguishing and correlating ground motion parameters of the upper and lower reservoirs of a large pumped storage power station. The method includes the following steps: acquiring three-dimensional wave velocity and coordinate data of the upper and lower reservoirs through remote sensing mapping and drilling tests; modeling and simulating the wave field to output acceleration response; extracting phase abrupt change points, calculating density functions to generate a comparison table; statistically analyzing the offset and density distribution of abrupt change points; calculating normalized offset rates to form a discriminant dataset; extracting peak acceleration and response spectrum parameters to construct a feature set; training a random forest regression model to establish a mapping relationship; and introducing topographic propagation coefficients for calibration to generate a correlation table. In this invention, by jointly modeling topography and wave velocity to reconstruct the influence of structural differences on seismic wave paths, segmenting acceleration sequences to analyze abrupt change characteristics, and extracting peak values ​​and spectral points to quantify response differences, the method improves the analytical capability and matching degree of ground motion characteristics.

Description

Technical Field

[0001] This invention relates to the field of seismic wave analysis technology, and in particular to a method for distinguishing and correlating ground motion parameters of the upper and lower reservoirs of a large pumped storage power station. Background Technology

[0002] The field of seismic wave analysis technology involves the acquisition, identification, decomposition, and parameter extraction of seismic waves generated during seismic activity to reveal the propagation characteristics of seismic waves, the response of geological structures, and the reaction patterns of engineering facilities to seismic waves. Core aspects of this technology include the determination and analysis of ground motion parameters, the inversion of seismic wave propagation paths, the identification of seismic wave spectral characteristics, the simulation and calculation of seismic wave fields, and the dynamic analysis of structural responses. This field is widely used in engineering seismic safety assessments, seismic fortification research of major infrastructure, geological structural exploration, and earthquake disaster prediction. It is highly interdisciplinary, involving theoretical support and practical methods from multiple disciplines such as geophysics, geotechnical engineering, and structural dynamics.

[0003] The traditional method for distinguishing and correlating seismic motion parameters between the upper and lower reservoirs of large-scale pumped storage power stations involves separating and identifying the seismic response characteristics of the upper and lower reservoirs based on seismic monitoring data during seismic wave analysis, and further establishing the correlation between their parameters. This technique primarily addresses the issue of significant differences in the amplitude, frequency, and duration of seismic motions during propagation when pumped storage power stations are constructed in mountainous or hilly terrain, due to the substantial differences in terrain conditions between the upper and lower reservoirs. Traditional methods use single-point seismic monitoring data combined with empirical wave selection criteria and spectral matching analysis. By comparing and analyzing the acceleration, velocity, and displacement time-history waveform characteristics in existing seismic records from the upper and lower reservoirs, and combining the response spectrum calculation results, the differences in their seismic motion parameters are determined, and empirical or statistical correlations are established accordingly.

[0004] Existing technologies rely on empirical rules to identify features in seismic ground motion data from upper and lower reservoirs. However, due to limitations in the number of monitoring points and terrain complexity, they cannot effectively distinguish and represent seismic waves when there are spatial variations in their propagation paths. In particular, in scenarios where there are abrupt changes in geological structure or significant changes in wave velocity gradients between upper and lower reservoirs, time-history data are prone to misalignment of key time periods and waveform distortion. The spectral content lacks a stable reference basis, resulting in non-universal applicability of the established parameter correlations. This further affects the accuracy of ground motion feature extraction and the effectiveness of response difference determination, creating a significant risk of error in the selection of seismic fortification design parameters. Summary of the Invention

[0005] To address the limitations of existing technologies that rely on empirical rules for feature identification of upper and lower reservoir ground motion data, such as the inability to effectively distinguish and represent data when seismic wave propagation paths exhibit spatial variability due to the limited number of monitoring points and terrain complexity, especially in scenarios with abrupt geological changes or significant variations in wave velocity gradients between the upper and lower reservoirs, time-history data is prone to misalignment of key time periods and waveform distortion. Furthermore, the lack of a stable reference basis for spectral content results in non-universal parameter correlations, further impacting the accuracy of ground motion feature extraction and the effectiveness of response difference determination. This poses a significant error risk in the selection of seismic fortification design parameters. Therefore, this invention provides a method for distinguishing and correlating ground motion parameters between the upper and lower reservoirs of a large pumped storage power station. The technical solution is as follows:

[0006] On the one hand, a method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of large pumped storage power stations is provided. This method includes:

[0007] S1: Obtain three-dimensional coordinate data of the upper and lower reservoirs through remote sensing mapping, obtain wave velocity structure data through drilling tests, construct a three-dimensional geological model by geological body layering based on the geological profile and lithological characteristics of the engineering area, mesh the model based on the three-dimensional wave finite element method, carry out wave field simulation, and output acceleration response sequence.

[0008] S2: Divide the acceleration response sequence into time windows, extract phase change points and statistically analyze the time series, use the phase change point density function to calculate the unit time distribution, and generate a phase change point time series comparison table;

[0009] S3: Based on the phase mutation point time series comparison table, calculate the time series offset of the upper and lower database mutation points within the same time window, filter the concentrated distribution period of mutation points based on the density threshold, and combine the normalized offset rate of the response duration of the upper and lower databases to calculate their joint index and obtain the upper and lower database cointegration discrimination dataset.

[0010] S4: Based on the upper and lower library cointegration discrimination dataset, mark the time window segments that meet the thresholds of time offset, time period intersection ratio, and duration offset rate, extract the acceleration peak value and key period point values ​​of the response spectrum to obtain the differential parameter feature set.

[0011] As a further aspect of the present invention, the phase change point density function quantifies the distribution characteristics of phase change points per unit time, revealing the drastic changes in structural response and their temporal patterns.

[0012] The three-dimensional geological model includes geological interface elevation data, elastic modulus parameters, and grid division parameters. The acceleration response sequence includes time-domain waveforms, Fourier spectral components, and energy attenuation coefficients. The phase abrupt change point time series comparison table includes phase abrupt change frequency, time interval standard deviation, and cumulative frequency histogram. The upper and lower database cointegration discrimination dataset includes the mean time series offset, time period overlap rate, and normalized offset rate threshold. The differential parameter feature set includes peak acceleration spectrum values, feature period point set, and differential index matrix.

[0013] As a further aspect of the present invention, the specific steps of S1 include:

[0014] S101: Acquire remote sensing images and topographic mapping data, combine lidar scanning and photogrammetry results to extract three-dimensional coordinate points, perform spatial registration and digitization processing, and generate topographic spatial node coordinate values.

[0015] S102: Drilling points are set according to the coordinate values ​​of the topographic spatial nodes, layer information and wave velocity data are collected and the profile is reconstructed. Based on the geological profile and lithological characteristics, the three-dimensional geological body is layered and modeled. Wave velocity mapping is completed by combining the three-dimensional wave finite element modeling method to obtain the node wave velocity distribution value.

[0016] The nodal wave velocity distribution value refers to the numerical distribution of the velocity of the medium at the node in the three-dimensional model in response to the propagation of seismic waves, reflecting the response capability of the geological medium at the node to wave propagation;

[0017] S103: Construct a three-dimensional wave velocity structure model based on the nodal wave velocity distribution values, set boundary conditions and source parameters, perform wave propagation simulation, and obtain the acceleration response sequence.

[0018] As a further aspect of the present invention, the specific steps of S2 include:

[0019] S201: Based on the acceleration response sequence, a fixed time interval is set to divide the time into multiple independent time windows. The response value changes within the window are scanned sequentially, phase change points are identified and marked with time positions, change time data are recorded, and a phase change time list is generated.

[0020] S202: Call the phase mutation time list, divide adjacent time intervals into segments by unit time, count the number of mutations in multiple segments and calculate the density value, and generate a mutation point distribution density sequence.

[0021] S203: Based on the mutation point distribution density sequence, compare the density value similar segments, screen synchronous mutation time points and record the sequence number, and generate a phase mutation point time series comparison table;

[0022] The phase abrupt change point refers to the instantaneous point at which the phase continuity is broken in the wave response sequence.

[0023] As a further aspect of the present invention, the specific steps of S3 include:

[0024] S301: Based on the phase mutation point time series comparison table, mutation point time series matching is performed according to the same time window conditions. The time offset of the upper and lower library mutation points in each time window is extracted as the mutation point offset. The differentiated phase windows are processed in segments and multiple offset sequences are calculated to obtain the upper and lower library mutation time offset sequences.

[0025] S302: Call the mutation time offset sequence of the upper and lower libraries, count the frequency of mutation points within the time period, filter the concentrated distribution segments with frequencies exceeding the benchmark according to the mutation point density threshold, remove low-density segments, obtain the concentrated distribution interval of mutation points in the upper and lower libraries, and generate the concentrated distribution time period interval value.

[0026] The recommended range for the mutation point density threshold is 2 to 5 mutation points / second, which is selected based on real-time engineering monitoring and statistical experience.

[0027] S303: Based on the values ​​of the concentrated distribution time period interval, and taking into account the normalized offset rate of the mutation points of the upper and lower databases within the concentrated distribution interval, generate a cointegration discrimination dataset for the upper and lower databases according to the weighting method.

[0028] The normalized offset rate refers to the ratio of the upper and lower database response time offset to the reference response duration, and is dimensionless.

[0029] As a further aspect of the present invention, the specific steps of S4 include:

[0030] S401: Based on the cointegration discrimination dataset of the upper and lower libraries, calculate the time tag offset, filter out samples that exceed the time offset threshold, judge according to the intersection ratio of time periods and the duration offset rate, and generate the offset time window filtering interval.

[0031] The recommended range for the time offset threshold is 0.1s to 0.5s, which is selected based on real-time engineering monitoring statistical experience.

[0032] S402: Call the sample acceleration sequence in the offset time window screening interval, extract the continuous change segment based on the fixed window, identify the peak position, and perform moving average calculation on the acceleration fluctuation value after Z-score normalization, calculate the peak amplitude and the periodic fluctuation change rate parameter between adjacent peaks in the sequence, and obtain the peak time series parameter quantity.

[0033] S403: Based on the peak time in the peak time series parameters, obtain the corresponding time response spectrum sequence, calculate the mean square error and offset rate of the response values ​​within the main period segment, perform matching processing in combination with the sequence number, and integrate to generate a differentiated parameter feature set.

[0034] As a further aspect of the present invention, the periodic fluctuation rate of change parameter is calculated using the following formula:

[0035] ;

[0036] in, The parameter representing the rate of change of periodic fluctuations, Represents the total number of acceleration peaks. Representing the The acceleration amplitude of each peak is a dimensionless parameter after Z-score normalization. Representing the The acceleration amplitude of each peak is a dimensionless parameter after Z-score normalization. Representing the The time labels corresponding to each peak are dimensionless parameters after Z-score normalization. Representing the The time labels corresponding to each peak are dimensionless parameters after Z-score normalization. Representative of the first The peaks represent the average acceleration within the central window, and are dimensionless parameters after Z-score normalization. Representative of the first The peaks represent the average acceleration within the central window, and are dimensionless parameters after Z-score normalization. To avoid the dimensionless minimum constant term with a denominator of zero.

[0037] As a further aspect of the present invention, the method includes step S5:

[0038] S5: Based on the differential parameter feature set, construct a training dataset with the response difference between the upper and lower libraries as input, use a random forest regression model to learn the nonlinear mapping relationship between the mutation point density and the response difference, form a cointegration discrimination model, and introduce the terrain propagation path difference coefficient for model calibration to generate an association mapping table.

[0039] The cointegration discrimination model includes decision tree depth parameters, feature criticality ranking, and terrain attenuation factor. The association mapping table includes propagation path correction coefficients, difference weight matrix, and model calibration parameters.

[0040] As a further aspect of the present invention, the specific steps of S5 include:

[0041] S501: Based on the differential parameter feature set, extract the response frequency difference value and the earthquake energy ratio, construct a differential dataset and normalize it, train a random forest model, and generate a predicted value for the density of mutation points.

[0042] S502: Call the predicted value of mutation point density, match the terrain grid path according to the multidimensional response index residual, calculate the propagation path length difference and complete the adjustment, and generate the path length difference correction value.

[0043] S503: Based on the path length difference correction value, select the density offset and response difference interval, perform nonlinear surface fitting to construct a mapping relationship, integrate the fitting results and match the initial parameter feature set to generate an association mapping table.

[0044] As a further aspect of the present invention, the propagation path length difference is calculated using the following formula:

[0045] ;

[0046] in, This represents the difference in propagation path length, in meters (m). This represents the number of nodes involved in the calculation within the terrain grid path. Representative node The residual value of the multidimensional response index at the location, in m / s 2 , Representative node The ratio of earthquake energy at a location Representative node The difference in terrain elevation at each location, in meters. Representative node The corresponding predicted mutation point density, in points / km 2 , The mean of the predicted mutation point density within the path, points / km 2 , The standard deviation of the predicted density of mutation points within the path, in points / km. 2 .

[0047] The beneficial effects of the technical solutions provided by the embodiments of the present invention include at least the following:

[0048] By jointly modeling topographic data and wave velocity structure, a three-dimensional geological scene with spatial resolution is generated, which can restore the impact of the differences in the actual topographic structure between the upper and lower reservoirs on the propagation path of seismic waves. Combined with the refined segmentation of acceleration response time series and the assessment of the density distribution of phase change points, the ability to capture the change law of response characteristics is enhanced. By comparing and analyzing the time offset and duration normalization deviation of the response change points between the upper and lower reservoirs and setting a discrimination threshold to screen key time windows, the ability to quantitatively extract non-uniform propagation characteristics is enhanced. The extracted parameter peaks and spectral points cover key response sections, which can fully express the differences in the response characteristics of ground motion under different tectonic backgrounds, and significantly improve the analytical ability of ground motion response law and the matching degree of characteristic parameters. Attached Figure Description

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

[0050] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0051] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0052] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning.

[0053] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0054] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0055] Please see Figure 1 This invention provides a method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station. The processing flow of this method may include the following steps:

[0056] S1: Obtain three-dimensional coordinate data of the upper and lower reservoirs through remote sensing mapping, obtain wave velocity structure data through drilling tests, construct a three-dimensional geological model by geological body layering based on the geological profile and lithological characteristics of the engineering area, mesh the model based on the three-dimensional wave finite element method, carry out wave field simulation, and output acceleration response sequence.

[0057] S2: Divide the acceleration response sequence into time windows, extract phase change points and statistically analyze the time series, use the phase change point density function to calculate the unit time distribution, and generate a phase change point time series comparison table;

[0058] S3: Based on the phase mutation point time series comparison table, calculate the time series offset of mutation points in the upper and lower libraries within the same time window, filter the concentrated distribution period of mutation points based on the density threshold, and combine the normalized offset rate of the response duration of the upper and lower libraries to calculate their joint index and obtain the upper and lower library cointegration discrimination dataset.

[0059] S4: Based on the cointegration discrimination dataset of the upper and lower libraries, mark the time windows that meet the thresholds of time offset, time period intersection ratio, and duration offset rate, extract the acceleration peak value and key period point values ​​of the response spectrum, and obtain the differential parameter feature set;

[0060] S5: Based on the differential parameter feature set, a training dataset with the response difference between the upper and lower databases as input is constructed. A random forest regression model is used to learn the nonlinear mapping relationship between the mutation point density and the response difference, forming a cointegration discrimination model. The terrain propagation path difference coefficient is introduced for model calibration, and an association mapping table is generated.

[0061] The three-dimensional geological model includes geological interface elevation data, elastic modulus parameters, and grid division parameters. The acceleration response sequence includes time-domain waveforms, Fourier spectral components, and energy attenuation coefficients. The phase abrupt change point time-series comparison table includes phase abrupt change frequency, time interval standard deviation, and cumulative frequency histogram. The upper and lower database cointegration discrimination datasets include the mean time-series offset, time period overlap rate, and normalized offset rate threshold. The differential parameter feature set includes peak acceleration spectrum values, feature periodic point set, and differential index matrix. The cointegration discrimination model includes decision tree depth parameters, feature importance ranking, and terrain attenuation factor. The association mapping table includes propagation path correction coefficients, differential weight matrix, and model calibration parameters.

[0062] Specifically, the steps of S1 are as follows:

[0063] S101: Acquire remote sensing images and topographic mapping data, combine lidar scanning and photogrammetry results to extract three-dimensional coordinate points, perform spatial registration and digitization processing, and generate topographic spatial node coordinate values.

[0064] In the process of acquiring remote sensing imagery and topographic mapping data, it is necessary to first obtain high-resolution remote sensing images covering the upper and lower reservoir areas of the power station, and then use a GNSS real-time dynamic differential positioning system to collect control point coordinates. The image data is then geo-registered to ensure a clear correspondence between the image and the actual ground in the geographic coordinate system. Next, multiple scans are performed using a UAV equipped with a LiDAR system to acquire dense point cloud data. Combined with the deployment of ground control points, multi-view photogrammetry is used to acquire real-scene images. High-density three-dimensional coordinates are obtained through structured light reconstruction and image matching techniques. The point set and coordinate point set undergo error correction processing, including removing redundant edge data of the collected points and correcting elevation abrupt changes, to ensure that the spatial accuracy meets the terrain accuracy control within a range of ±0.05m. Then, representative nodes are extracted from key locations such as the upper reservoir dam crest, lower reservoir dam foundation, and the inlet and outlet. The three-dimensional coordinates of each node are expressed in latitude and longitude (WGS84) and geodetic elevation, and gridded at a spatial resolution of 0.5m × 0.5m to generate a digital terrain model (DTM). For example, for the upper reservoir area of ​​a power station, a 10km... 2 The area contains approximately 12 million terrain points. After spatial registration, 3,000 feature nodes were selected as analysis targets and labeled accordingly. , , The values ​​and specific location information are shown in Table 1.

[0065] Table 1: Three-dimensional coordinates of feature nodes in the upper database

[0066]

[0067] As shown in Table 1, the extracted nodes are all distributed in the area near the top of the upper reservoir dam. Some nodes are used as a basis for subsequent drilling point layout. After the coordinate data is processed, the points are imported into the GIS system for spatial relationship verification. The relative elevation difference and slope relationship between multiple points are analyzed. For example, if the elevation difference between nodes N001 and N005 is determined to be 625.32-620.13=5.19m and the horizontal distance is 57.24m, then the slope is... This method filters areas with a slope greater than 0.05 as potentially unstable slope areas, further providing reference information for drilling point placement. Spatial registration error verification is performed using control point residuals. If the residual between the measured control point value and the corresponding point in the image exceeds 0.1m, recalibration is required. The judgment rule is: assuming the measured point elevation is... The elevation of the image points is ,when When m is reached, outliers are considered to exist and need to be removed; after the spatial node coordinate values ​​are generated, the export format is uniformly adopted as GeoTIFF and LAS point cloud dual format to be compatible with the direct reading of subsequent finite element mesh modeling tools.

[0068] S102: Drilling points are set according to the coordinate values ​​of the topographic spatial nodes, layer information and wave velocity data are collected and the profile is reconstructed. Based on the geological profile and lithological characteristics, the three-dimensional geological body is layered and modeled. Wave velocity mapping is completed by combining the three-dimensional wave finite element modeling method to obtain the node wave velocity distribution value.

[0069] In setting drilling points based on the coordinates of topographic spatial nodes, it is necessary to combine the key nodes marked in the aforementioned digital terrain model and prioritize areas with terrain slopes between 0.05 and 0.25 and high terrain abrupt changes at 0.5km intervals. During the point setting process, first, the data is imported into the GIS system to generate a point layer. Then, representative areas are selected based on the slope stability analysis results. For example, in the upper reservoir dam crest edge area, areas with elevation gradients exceeding 3.5m / 50m are selected to set points. Subsequently, drilling operations are performed for each point, with a drilling depth of 25m. Soil and rock samples are collected every 2m, recording the lithology, porosity, and moisture content of the corresponding layers. Wave velocity is also collected using a wave velocity testing device. The acquisition method adopts a vertical source excitation and receiver synchronous monitoring approach. Five P-wave and S-wave propagation times are collected at each depth point, and the average value is taken as the wave velocity parameter for that point. For example, at a drilling depth of 12m, the recorded P-wave propagation time is 1.12ms and the distance is 2.5m; therefore, the P-wave velocity is:

[0070] ;

[0071] The wave velocity at depth points is calculated in this manner to form a one-dimensional wave velocity profile data structure. Then, interpolation is performed on a spatial profile point set formed by multiple borehole points. A trilinear interpolation method is used to construct a continuous wave velocity field based on the wave velocity values ​​of adjacent points in space, ensuring the model's complete continuity in space. Simultaneously, wave velocity data mapping is achieved by matching the wave velocity values ​​of spatial nodes with their corresponding coordinates. To avoid the impact of boundary disturbances on data validity, interpolation filtering conditions are set: when the standard deviation of the wave velocity of any node in its three neighborhoods exceeds 300 m / s, the data for that point is discarded and re-interpolated. After the wave velocity data is processed, it is imported into the three-dimensional wave propagation modeling module, and the wave velocity of each spatial node is... Import the data mesh and construct a three-dimensional wave finite element unit. In the mesh generation, the side length of each element is set to 5m to ensure matching with the feature preservation under the original 0.5m spatial resolution. At the same time, the medium parameters within each element are set to be consistent, representing the longitudinal wave velocity. transverse wave velocity and density After constructing the model, boundary absorption conditions are set. For example, a free field boundary is used to set the vibration propagation direction to the negative Z-axis, the source point is set at an elevation of 635m in the middle of the upper reservoir dam crest, the excitation frequency is 20Hz, the duration is set to 2s, the source form is Ricker wavelet, and the simulation step size is 0.001s. Finally, the simulation is executed to complete the mapping of wave velocity values ​​at each node and output them in a three-dimensional vector format for easy use in subsequent seismic response calculations.

[0072] S103: Construct a three-dimensional wave velocity structure model based on the nodal wave velocity distribution values, set boundary conditions and source parameters, perform wave propagation simulation, and obtain the acceleration response sequence;

[0073] In constructing a three-dimensional wave velocity structure model based on node wave velocity distribution values, the first step is to standardize the data structure of the generated three-dimensional coordinates and wave velocity data, with each node recording... The information for a six-tuple, with the number of nodes controlled between 10,000 and 30,000, is organized hierarchically using a spatial octree structure for easier subsequent loading and computation. For example, the information for a certain node is... This represents a longitudinal wave velocity of 2360 m / s, a transverse wave velocity of 1320 m / s, and a medium density of 2.35 g / cm³. 3 After completing the model construction, boundary conditions need to be set. Fully absorbing boundary conditions are used at the upper and lower boundaries to avoid interference from reflected waves in the simulation results. Free boundary conditions are set at the side boundaries to reflect the real-time free response behavior of the slope. Simultaneously, source parameters are set, selecting the source location as nodes (118.924013, 29.342741, 623.55). The source excitation mode is a Ricker wave with a center frequency of 20Hz, an initial amplitude of 1.0mm, and a vertically downward vibration direction. Response values ​​are recorded every 0.01s in the model. During the simulation, the explicit finite difference method is used to solve the wave propagation process. The velocity and acceleration responses at key nodes are recorded at each time step. The response intensity is calculated based on the total energy received by each node. The energy calculation method is as follows: the acceleration value is integrated over time to obtain the velocity, which is then multiplied by the mass to obtain the kinetic energy. The kinetic energy value is then combined with the number of nodes to output a spatial energy distribution map. For example, the acceleration of a node at t=0.25s is 0.35m / s². 2 Given a mass of 2.35 kg (converted from density and mesh volume), the kinetic energy is:

[0074] ;

[0075] If the velocity obtained after integration is 0.22 m / s, then the kinetic energy is... This process is repeated to obtain the acceleration response sequence. The response results of each node are exported as a time series file in CSV format, with each column representing the time, acceleration X, Y, and Z components. This results in a complete seismic response database for the analysis and differentiation of ground motion parameters in the upper and lower reservoir areas.

[0076] Specifically, the steps of S2 are as follows:

[0077] S201: Based on the acceleration response sequence, a fixed time interval is set to divide the time into multiple independent time windows. The response value changes within the window are scanned sequentially, phase change points are identified and marked with time positions, change time data are recorded, and a list of phase change times is generated.

[0078] Based on the acceleration response sequence, the complete response data is divided into multiple independent time windows with fixed time intervals, each window length set to 0.5 seconds. This is suitable for the real-time response frequency variation range of the upper and lower reservoirs of a pumped storage power station. Before windowing processing, the raw acceleration data is first subjected to baseline correction and filtering to remove zero drift and high-frequency noise before being stored in the processing matrix. The acceleration data within each time window is scanned sequentially, and an instantaneous phase sequence is constructed using complex analytic signals. The instantaneous phase value at each moment is obtained using Hilbert transform. The absolute value of the phase difference between consecutive time points is calculated, and a phase change identification threshold is set. That is, if the phase difference between any two consecutive time points is greater than If a point is identified as a phase abrupt change, this threshold is determined based on the frequency changes in the seismic excitation response, as shown in Table 2. The time locations corresponding to these abrupt changes are further recorded, and the set of time locations is stored in an array. If the sampling frequency within a certain time window is 100Hz, then a 0.5-second window will contain 50 data points. By judging the phase abrupt change at each point, a sequence of abrupt change time points can be formed. For example, if an abrupt change point occurs in a time series... If the mutation point is not found, it is added to the mutation time list. During the mutation point marking process, the phase difference operation is achieved by calculating the phase difference pairwise from consecutive data points and comparing it with a threshold. For example, when the th mutation point is found to be a mutation point, the phase difference is calculated pairwise from consecutive data points and compared with a threshold. The instantaneous phase of each data point is , No. The data points are Then execute If the judgment operation is true, it is marked as a mutation point, and... corresponding time point The data is recorded in a list, resulting in a final list of phase transition times summarized from multi-window scans.

[0079] Table 2: Statistical Table of Phase Abrupt Change Thresholds (Unit: radians)

[0080]

[0081] As shown in Table 2, the statistics of the maximum phase difference values ​​of the upper and lower reservoir acceleration response signals under typical seismic conditions show that the phase difference at most abrupt change points exceeds [a certain value]. Therefore, it is appropriate to set the mutation judgment threshold as follows: Under this setting, the accuracy of mutation identification can reach over 93%. This result indicates that the identification of phase mutation points depends on point-by-point detection of instantaneous phase changes in the acceleration sequence, and the setting of the phase difference threshold has a crucial impact on the determination of mutation points, thus affecting subsequent mutation density calculations and phase point comparisons.

[0082] S202: Call the phase mutation time list, divide adjacent time intervals into segments by unit time, count the number of mutations in multiple segments and calculate the density value, and generate a mutation point distribution density sequence.

[0083] Call the obtained list of phase transition times The mutation points are sorted according to the time series to establish a unified time axis. On this time axis, the unit time interval is set to 1 second. The entire time interval covered by the acceleration response sequence is divided. Assuming the total duration of the current response record is 60 seconds, it is divided into 60 time intervals, each numbered... Then iterate through each time point in the mutation time list. Determine the time period number in which it occurs, i.e., if Seconds are counted in the first... In counting mutations over time periods, this process is repeated until the entire list of mutation times has been traversed, resulting in a one-dimensional integer array of length 60. Each element Indicates the first The number of mutations recorded within a second is further normalized into a density sequence. Each density value is obtained through the formula Calculation, where A second is the length of a unit of time interval, therefore the density value here corresponds to the number of mutations. For example, if a certain time interval... If the number of internal mutations is 4, then the corresponding density is: After calculating the density values, the density sequences are compared by difference to analyze the density trend changes over multiple consecutive time periods. In this process, a high density interval (greater than 3.0) and a low density interval (less than or equal to 1.0) are defined. The judgment criteria are based on previous test statistics and earthquake frequency characteristics. If the density value is greater than 3.0 in three consecutive time periods, the segment is marked as a dense abrupt change segment; if the density is less than 1.0 in three consecutive time periods, it is marked as a sparse abrupt change segment. For example, if the density value is less than 3.0 in three consecutive time periods... If the density value is high, it is considered to be a period of concentrated synchronous strong response events. This period can be used as a basis for subsequent synchronous mutation analysis. The density value and the corresponding time period number are stored in a two-dimensional array in sequence to form a mutation density distribution sequence table. This table will serve as the basic data for screening synchronous mutation points in the next stage.

[0084] S203: Based on the distribution density sequence of mutation points, compare segments with similar density values, screen synchronous mutation time points and record sequence numbers to generate a phase mutation point time series comparison table.

[0085] A phase abrupt change point refers to the instantaneous point at which the phase continuity is broken in a wave response sequence;

[0086] Based on the mutation point distribution density sequence obtained in the previous step The density values ​​of different time periods are compared to identify time segments with highly similar density values. Cross-analysis is then performed on their abrupt change points. First, the density values ​​are iterated through by time period number. For each time period, the density value is compared with its two adjacent time periods, and the absolute value of the difference is calculated. and Set a similarity threshold If both are less than Then determine the current time period. These three time periods, which have highly similar densities to their adjacent time periods, are added to the comparison sequence as a density similarity group. Further comparisons are made with the corresponding set of mutation time points within this density similarity group, recording the specific time value of each mutation point and its original sampling sequence number. The system determines whether there are overlapping mutation points between two time periods. Let the two time periods be... and The mutation time sets are respectively and If there exists any ,in If a synchronous mutation point exists between two time periods, it is considered to be a synchronization point, and a synchronous mutation time point pair is constructed. The time period number of each time period is stored in the phase change point time series comparison table, thus completing the work of screening synchronous change time points and collecting sequence numbers. For example, when and For periods of similar density, each includes a mutation time point. and Within the set error range If there is no synchronization point between the two, it will not be recorded. If in a certain sample and If a pair is identified as a synchronous mutation pair, then the pair at that time point... , The results are stored in the results table, ultimately generating a complete set of phase abrupt change comparison results.

[0087] Specifically, the steps of S3 are as follows:

[0088] S301: Based on the phase mutation point time series comparison table, mutation point time series matching is performed according to the same time window conditions. The time offset of the upper and lower library mutation points in each time window is extracted as the mutation point offset. The differential phase window is processed in segments and multiple offset sequences are calculated to obtain the upper and lower library mutation time offset sequences.

[0089] Based on the phase abrupt change point time series comparison table, the abrupt change point pairs are grouped and organized according to their corresponding time period numbers, and then the same time window condition is set as follows. The time information of mutation points in the upper and lower libraries within each group is extracted and denoted as follows: and For each pair of mutation points, the time offset is calculated separately, i.e., obtained through interpolation. To improve processing efficiency, mutation point pairs are pre-indexed according to time sequence. Let the times of the first pair of mutation points be... , The offset is then calculated as Record this value at the beginning of the offset sequence, and repeat the above interpolation operation sequentially for each pair of mutation points to form a complete offset sequence. The response data, totaling 60 seconds, was then divided into 60 time segments at 1-second intervals, each segment denoted as _____. to Traverse the time points in the offset sequence, record the time interval to which they belong, forming a set of offset values ​​per second. For example, if the offset time point is at 23.88 seconds, it belongs to... The offset values ​​are stored in the 23rd data segment. Local set extraction is performed on each data segment to establish multiple time period offset segment subsets. The offset data within each segment is grouped by iterative traversal to finally obtain 60 sets of time segment offset values. Combined with each set of offset values, the time offset distribution sequence of the mutation points of the upper and lower databases is constructed.

[0090] S302: Call the mutation time offset sequence of upper and lower databases, count the frequency of mutation points within the time period, filter the concentrated distribution segments with frequencies exceeding the baseline based on the mutation point density threshold, remove low-density segments, obtain the concentrated distribution interval of mutation points in upper and lower databases, and generate the concentrated distribution time period interval value.

[0091] The recommended range for the mutation point density threshold is 2 to 5 mutation points / second, which is selected based on real-time engineering monitoring and statistical experience.

[0092] The generated time offset distribution sequence is called to perform statistical operations on the frequency of abrupt changes within a time period. Specifically, this involves counting the number of abrupt change pairs in each offset set. A threshold range for abrupt change density is set to 2 to 5 abrupt changes per second. This threshold is based on statistical analysis of earthquake monitoring data from 12 pumped storage power stations in China over the past three years. In 95% of valid earthquake events, the abrupt change density in both the upper and lower reservoirs was distributed between 2.5 and 4.8 abrupt changes per second. Therefore, this threshold range is determined by rounding down. The frequency of abrupt changes in each segment is then analyzed. If the number of abrupt change pairs in a given time period is greater than or equal to 2, it is considered a concentrated segment; otherwise, it is considered a low-density segment. This process completes the filtering of abrupt change density within a time period. Further, the set of numbers for the concentrated time periods is extracted. For example, if there are 3 abrupt change pairs in the 12th second segment, the corresponding offset value is... The frequency of this segment is 3, which meets the density threshold condition. This segment is recorded as a concentrated distribution segment. This process is repeated for 60 time periods to perform density screening, ultimately forming a set of concentrated distribution time period interval values. The specific time period numbers and the frequency of their mutation points are shown in the table below:

[0093] Table 3: Statistics on Concentrated Distribution Time Periods

[0094]

[0095] As shown in Table 3, the time periods that meet the mutation point density threshold condition are marked as concentrated distribution segments, and only such time periods are retained for cointegration analysis in subsequent analyses.

[0096] S303: Based on the values ​​of the concentrated distribution period interval, and taking into account the normalized offset rate of the mutation points of the upper and lower databases within the concentrated distribution interval, generate the upper and lower database cointegration discrimination dataset according to the weighting method.

[0097] Normalized offset rate refers to the ratio of the upper and lower database response time offset to the reference response duration, and is dimensionless.

[0098] Based on the obtained concentrated distribution time interval values For each segment, the offset value of the recorded mutation point Normalization is performed using the following formula:

[0099] ;

[0100] in, Indicates the normalized offset rate. For a given mutation point pair, the absolute offset time is given. For reference response duration, this example uses... This makes the time offset dimensionless. For example, if a certain abrupt change point corresponds to an offset of... The normalized offset rate is After normalizing the offset values ​​within each concentrated time segment, a normalized offset rate sequence is generated for each corresponding segment. Let the set of offset values ​​in the 12th second segment be... The normalized offsets are then:

[0101] ;

[0102] Record the normalization results, and then generate the upper and lower database cointegration discrimination dataset using a weighted average method. Let the normalized offset rate sequence of each segment be denoted as . The weighting method uses an equal-weighted summation, meaning the mean of each segment is:

[0103] ;

[0104] Taking the 12-second segment as an example, its average normalized offset is:

[0105] ;

[0106] The same operation is performed on the concentrated distribution segments, ultimately generating a set of cointegration offset evaluation values ​​for the corresponding concentrated segments. This result can serve as the basic input dataset for cointegration analysis of mutation points in the upper and lower libraries during the synchronous response period. The result indicates that time periods with smaller normalized offsets (e.g., below 0.001) suggest smaller differences in response times between the upper and lower libraries, providing a reference for subsequent screening of synchronous response events.

[0107] Specifically, the steps of S4 are as follows:

[0108] S401: Based on the cointegration discrimination dataset of upper and lower libraries, calculate the time label offset, filter out samples that exceed the time offset threshold, judge according to the intersection ratio of time periods and the duration offset rate, and generate the offset time window filtering interval.

[0109] The recommended range for the time offset threshold is 0.1s to 0.5s, which is selected based on experience in real-time engineering monitoring statistics.

[0110] Based on the upper and lower reservoir cointegration discrimination dataset, the acceleration sequences collected by the sensors in the upper and lower reservoirs need to be time-aligned first. This process begins by extracting the time label of each data sample and calculating its time offset between the two sequences. For example, if the upper reservoir sample time is 12.53 seconds and the lower reservoir sample time is 12.68 seconds, the offset is 0.15 seconds. This value is compared with a set time offset threshold of 0.3 seconds. Based on the statistical results of the seismic response time difference between the upper and lower reservoirs in the original monitoring data of a pumped storage power station, 95% of the sample offsets are distributed in the range of 0.1 seconds to 0.5 seconds. The median is taken as the recommended threshold. Since 0.15 seconds is less than the threshold, the sample is retained; otherwise, it is discarded. The discarding operation specifically involves iterating through the sample index list, performing conditional judgments, filtering out data points that do not meet the offset threshold, and recording their index numbers. Subsequently, in the dataset that meets the threshold condition, data points with the same index number from both sequences are extracted. The data segment is cross-referenced based on the percentage of intersection between time periods and the duration offset rate. The percentage of intersection is calculated by dividing the length of the intersection between two time periods by the minimum length of the time period. For example, if the time period for the upper database sample is 12.53-13.23 seconds and the lower database sample is 12.68-13.08 seconds, then the intersection is 12.68-13.08 seconds, with a length of 0.4 seconds. The minimum length of the time period is 0.55 seconds, and the percentage of intersection is 0.4 / 0.55≈72.7%. Meanwhile, the duration offset rate is calculated by dividing the difference between the lengths of the two time periods by their average value. The offset rate is calculated as (0.7-0.4) / ((0.7+0.4) / 2)=0.3 / 0.55≈54.5%. If the percentage of intersection is higher than 50% and the offset rate is lower than 60%, then the sample enters the offset time window filtering interval. Finally, a structured data list containing information such as sequence number, start and end time, and offset is formed as the filtering interval for the next step.

[0111] S402: Call the offset time window to filter the sample acceleration sequence in the interval, extract the continuously changing segment based on the fixed window, identify the peak position, and perform the moving average calculation on the acceleration fluctuation value after Z-score normalization. Calculate the peak amplitude and the periodic fluctuation change rate parameter between adjacent peaks in the sequence to obtain the peak time series parameter quantity.

[0112] The system uses an offset time window to filter samples within a specified interval. The acceleration sequences are segmented using a fixed time window, for example, with a window width of 0.2 seconds and a sliding step of 0.05 seconds. Data within the window is extracted sequentially from the starting point. Within each window, it checks for continuously varying segments. If a segment exhibits amplitude variations exceeding 1.5 times the standard deviation threshold, it is marked as a continuously varying segment. Local extrema are also detected as peak positions. The acceleration sequences are then processed using Z-score normalization to convert them into a standard normal distribution. This involves subtracting the mean from each data value and dividing by the standard deviation to obtain the original acceleration amplitude. Calculate the mean Standard deviation After standardization, it becomes: Original time series Calculate the mean Standard deviation The standardized time stamp is: Window acceleration mean Calculate the mean Standard deviation After standardization, it becomes: Next, the time interval between adjacent peaks and the difference in peak amplitude are identified, and the periodic fluctuation rate parameter is calculated using the formula:

[0113] ;

[0114] in, The parameter representing the rate of change of periodic fluctuations, Represents the total number of acceleration peaks. Representing the The acceleration amplitude of each peak is a dimensionless parameter after Z-score normalization. Representing the The acceleration amplitude of each peak is a dimensionless parameter after Z-score normalization. Representing the The time labels corresponding to each peak are dimensionless parameters after Z-score normalization. Representing the The time labels corresponding to each peak are dimensionless parameters after Z-score normalization. Representative of the first The peaks represent the average acceleration within the central window, and are dimensionless parameters after Z-score normalization. Representative of the first The peaks represent the average acceleration within the central window, and are dimensionless parameters after Z-score normalization. To avoid the dimensionless minimum constant term with a denominator of zero.

[0115] Taking peak 1 and peak 2 as examples, the standardized parameters are: , Substitute the smallest positive number into the calculation:

[0116] ;

[0117] The calculation result is the periodic fluctuation rate parameter between the two peaks. The average of multiple peaks can be used to characterize the overall fluctuation periodicity.

[0118] S403: Based on the peak time in the peak time series parameters, obtain the response spectrum sequence at the corresponding time, calculate the mean square error and offset rate of the response values ​​within the main period segment, perform matching processing in combination with the sequence number, and integrate to generate a set of differentiated parameter features;

[0119] Based on the peak times in the peak time series parameters, the original response spectrum sequence at the corresponding time is retrieved from the structural monitoring system. In practice, this is done using standardized time tags. Used as an index to locate real-time points in the original time-domain data. For example, when , , real-time Obtain the response spectrum data within a 0.1-second window before and after this moment, extract the acceleration response sequence of the main period segment (0.5-5Hz), assuming this segment contains 6 frequency points, and the normalized acceleration value is... Calculate its mean square error:

[0120] ;

[0121] Simultaneously calculate the offset rate and take the normalized time label of the current peak. Compared to the previous peak tag The difference ratio, for example, peak 2 ( ) relative to peak 1 ( The offset rate is:

[0122] ;

[0123] Finally, the peak numbers (P1, P2), root mean square error (0.86), and offset rate (88.9%) were integrated into structured entries to generate a complete set of differential parameter features. This process ensures that the calculations are based on standardized parameters, eliminating the influence of dimensional differences on the analysis results.

[0124] Specifically, the steps in S5 are as follows:

[0125] S501: Based on the differential parameter feature set, extract the response frequency difference value and the earthquake energy ratio, construct the differential dataset and normalize it, train the random forest model, and generate the predicted value of the mutation point density.

[0126] Extract the response frequency difference value and the earthquake energy ratio from the differential parameter feature set, and extract the frequency response corresponding to each wave peak from the differential parameter feature set. With earthquake energy An array, assuming peak numbers are 1 to 4, with the following associated data:

[0127] ;

[0128] Calculate the frequency difference value using an array of frequency differences between adjacent peaks:

[0129] ;

[0130] The earthquake energy ratio is calculated as follows:

[0131] ,

[0132] Construct a difference dataset in array form:

[0133] ;

[0134] The differential dataset was then normalized using min-max normalization. , The normalized value is calculated as follows:

[0135] ;

[0136] ;

[0137] Finally, the normalized array is... Input the random forest model for training and prediction, and generate an array of predicted mutation point densities, for example, the predicted values ​​are... .

[0138] S502: Call the predicted value of mutation point density, match the terrain grid path according to the residual of multidimensional response index, calculate the difference in propagation path length and complete the adjustment, and generate the path length difference correction value.

[0139] With the predicted value array As a node density value Call the multidimensional response residual corresponding to each path node Earthquake energy ratio Obtain the terrain height difference of the nodes. Example data is as follows:

[0140] ;

[0141] and the mean of predicted values ​​within the path Standard deviation ;

[0142] The propagation path length difference is calculated using the following formula:

[0143] ;

[0144] in, This represents the difference in propagation path length, in meters (m). This represents the number of nodes involved in the calculation within the terrain grid path. Representative node The residual value of the multidimensional response index at the location, in m / s 2 , Representative node The ratio of earthquake energy at a location Representative node The difference in terrain elevation at each location, in meters. Representative node The corresponding predicted mutation point density, in points / km 2 , The mean of the predicted mutation point density within the path, points / km 2 , The standard deviation of the predicted density of mutation points within the path is represented by points / km².

[0145] Calculate item by item:

[0146] ;

[0147] ;

[0148] ;

[0149] Average and take the absolute value:

[0150] ;

[0151] The difference in propagation path length is approximately 0.0089m, and this difference is then passed on to the next step.

[0152] S503: Based on the path length difference correction value, select the density offset and response difference interval, perform nonlinear surface fitting to construct a mapping relationship, integrate the fitting results and match the initial parameter feature set to generate an association mapping table;

[0153] Path length difference correction value and node prediction density offset array and the range of node response differences Perform nonlinear surface fitting and construct respectively based on Density offset PY For the three-dimensional scatter plot of the independent variable, a quadratic polynomial is used for fitting, and the fitting form is:

[0154] ;

[0155] Substitute the node data into the fit and construct a system of equations:

[0156] ;

[0157] ;

[0158] ;

[0159] Solve the above system of three quadratic equations, assuming the fitting coefficients are... Store it in the association mapping table and match it with the initial parameter feature set (including peak number, , The predicted density and other parameters are merged to form the final associated mapping item, which is used in the subsequent path correction process.

[0160] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for distinguishing and correlating seismic ground motion parameters of upper and lower reservoirs of large-scale pumped storage power stations, characterized in that, Includes the following steps: S1: Obtain three-dimensional coordinate data of the upper and lower reservoirs through remote sensing mapping, obtain wave velocity structure data through drilling tests, construct a three-dimensional geological model by geological body layering based on the geological profile and lithological characteristics of the engineering area, mesh the model based on the three-dimensional wave finite element method, carry out wave field simulation, and output acceleration response sequence. S2: Divide the acceleration response sequence into time windows, extract phase change points and statistically analyze the time series, use the phase change point density function to calculate the unit time distribution, and generate a phase change point time series comparison table; S3: Based on the phase mutation point time series comparison table, calculate the time series offset of the upper and lower database mutation points within the same time window, filter the concentrated distribution period of mutation points based on the density threshold, and combine the normalized offset rate of the response duration of the upper and lower databases to calculate their joint index and obtain the upper and lower database cointegration discrimination dataset. S4: Based on the upper and lower library cointegration discrimination dataset, mark the time windows that meet the thresholds of time offset, time period intersection ratio, and duration offset rate, extract the acceleration peak value and key period point values ​​of the response spectrum to obtain the differential parameter feature set; The method includes step S5: S5: Based on the differential parameter feature set, construct a training dataset with the response difference between the upper and lower libraries as input, use a random forest regression model to learn the nonlinear mapping relationship between the mutation point density and the response difference, form a cointegration discrimination model, and introduce the terrain propagation path difference coefficient for model calibration to generate an association mapping table. The cointegration discrimination model includes decision tree depth parameters, feature key ranking and terrain attenuation factor, and the association mapping table includes propagation path correction coefficients, difference weight matrix and model calibration parameters; The specific steps of S5 include: S501: Based on the differential parameter feature set, extract the response frequency difference value and the earthquake energy ratio, construct a differential dataset and normalize it, train a random forest model, and generate a predicted value for the density of mutation points. S502: Call the predicted value of mutation point density, match the terrain grid path according to the multidimensional response index residual, calculate the propagation path length difference and complete the adjustment, and generate the path length difference correction value. S503: Based on the path length difference correction value, select the density offset and response difference interval, perform nonlinear surface fitting to construct a mapping relationship, integrate the fitting results and match the initial parameter feature set to generate an association mapping table.

2. The method for distinguishing and correlating the ground motion parameters of upper and lower reservoirs of large pumped storage power stations according to claim 1, characterized in that, The phase change point density function quantifies the distribution characteristics of phase change points per unit time, revealing the dramatic changes in structural response and their temporal patterns. The three-dimensional geological model includes geological interface elevation data, elastic modulus parameters, and grid division parameters. The acceleration response sequence includes time-domain waveforms, Fourier spectral components, and energy attenuation coefficients. The phase change point time series comparison table includes phase change frequency, time interval standard deviation, and cumulative frequency histogram. The upper and lower database cointegration discrimination dataset includes the mean time series offset, time period overlap rate, and normalized offset rate threshold. The differential parameter feature set includes peak acceleration spectrum values, feature period point set, and differential index matrix.

3. The method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station according to claim 1, characterized in that, The specific steps of S1 include: S101: Acquire remote sensing images and topographic mapping data, combine lidar scanning and photogrammetry results to extract three-dimensional coordinate points, perform spatial registration and digitization processing, and generate topographic spatial node coordinate values. S102: Drilling points are set according to the coordinate values ​​of the topographic spatial nodes, layer information and wave velocity data are collected and the profile is reconstructed. Based on the geological profile and lithological characteristics, the three-dimensional geological body is layered and modeled. Wave velocity mapping is completed by combining the three-dimensional wave finite element modeling method to obtain the node wave velocity distribution value. S103: Construct a three-dimensional wave velocity structure model based on the nodal wave velocity distribution values, set boundary conditions and source parameters, perform wave propagation simulation, and obtain the acceleration response sequence.

4. The method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station according to claim 3, characterized in that, The specific steps of S2 include: S201: Based on the acceleration response sequence, a fixed time interval is set to divide the time into multiple independent time windows. The response value changes within the window are scanned sequentially, phase change points are identified and marked with time positions, change time data are recorded, and a phase change time list is generated. S202: Call the phase mutation time list, divide adjacent time intervals into segments by unit time, count the number of mutations in multiple segments and calculate the density value, and generate a mutation point distribution density sequence. S203: Based on the mutation point distribution density sequence, compare the density value similar segments, screen synchronous mutation time points and record the sequence number, and generate a phase mutation point time series comparison table; The phase abrupt change point refers to the instantaneous point at which the phase continuity is broken in the wave response sequence.

5. The method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station according to claim 4, characterized in that, The specific steps of S3 include: S301: Based on the phase mutation point time series comparison table, mutation point time series matching is performed according to the same time window conditions. The time offset of the upper and lower library mutation points in each time window is extracted as the mutation point offset. The differentiated phase windows are processed in segments and multiple offset sequences are calculated to obtain the upper and lower library mutation time offset sequences. S302: Call the mutation time offset sequence of the upper and lower libraries, count the frequency of mutation points within the time period, filter the concentrated distribution segments with frequencies exceeding the benchmark according to the mutation point density threshold, remove low-density segments, obtain the concentrated distribution interval of mutation points in the upper and lower libraries, and generate the concentrated distribution time period interval value. S303: Based on the values ​​of the concentrated distribution time period interval, and taking into account the normalized offset rate of the mutation points of the upper and lower databases within the concentrated distribution interval, generate a cointegration discrimination dataset for the upper and lower databases according to the weighting method. The normalized offset rate refers to the ratio of the upper and lower database response time offset to the reference response duration, and is dimensionless.

6. The method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station according to claim 5, characterized in that, The specific steps of S4 include: S401: Based on the cointegration discrimination dataset of the upper and lower libraries, calculate the time tag offset, filter out samples that exceed the time offset threshold, judge according to the intersection ratio of time periods and the duration offset rate, and generate the offset time window filtering interval. S402: Call the sample acceleration sequence in the offset time window screening interval, extract the continuous change segment based on the fixed window, identify the peak position, and perform moving average calculation on the acceleration fluctuation value after Z-score normalization, calculate the peak amplitude and the periodic fluctuation change rate parameter between adjacent peaks in the sequence, and obtain the peak time series parameter quantity. S403: Based on the peak time in the peak time series parameters, obtain the corresponding time response spectrum sequence, calculate the mean square error and offset rate of the response values ​​within the main period segment, perform matching processing in combination with the sequence number, and integrate to generate a differentiated parameter feature set.

7. The method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station according to claim 6, characterized in that, The periodic fluctuation rate of change parameter is calculated using the following formula: ； in, The parameter representing the rate of change of periodic fluctuations, Represents the total number of acceleration peaks. Representing the The acceleration amplitude of each peak is a dimensionless parameter after Z-score normalization. Representing the The acceleration amplitude of each peak is a dimensionless parameter after Z-score normalization. Representing the The time labels corresponding to each peak are dimensionless parameters after Z-score normalization. Representing the The time labels corresponding to each peak are dimensionless parameters after Z-score normalization. Representative of the first The peaks represent the average acceleration within the central window, and are dimensionless parameters after Z-score normalization. Representative of the first The peaks represent the average acceleration within the central window, and are dimensionless parameters after Z-score normalization. To avoid dimensionless constant terms with a denominator of zero.

8. The method for distinguishing and correlating the seismic motion parameters of the upper and lower reservoirs of a large pumped storage power station according to claim 1, characterized in that, The propagation path length difference is calculated using the following formula: ； in, This represents the difference in propagation path length, in meters (m). This represents the number of nodes involved in the calculation within the terrain grid path. Representative node The residual value of the multidimensional response index at the location, in m / s 2 , Representative node The ratio of earthquake energy at a location Representative node The difference in terrain elevation at each location, in meters. Representative node The corresponding predicted mutation point density, in points / km 2 , The mean of the predicted mutation point density within the path, points / km 2 , The standard deviation of the predicted density of mutation points within the path, in points / km. 2 .

Images