A Synergistic Method for Monitoring and Resource Utilization of Medium-Deep Groundwater

Abstract

This invention relates to the fields of geophysics and hydrogeochemistry, and discloses a collaborative method for monitoring and resource utilization of medium-deep groundwater. The method includes: acquiring physical field time-series data and hydrochemical isotope data from monitoring points; separating fault zone response signals using a bandpass filtering algorithm; constructing an augmented-dimensional coupled state vector based on the original physical parameter sequence and fault zone response signal sequence, and mapping it to a high-dimensional coupled state phase space to generate a physical field evolution trajectory; encoding hydrochemical and isotope data into geochemical fingerprint vectors and projecting them onto a mixing space to generate a chemical field evolution trajectory; calculating feature sequences from the two trajectories to generate physical trajectory stability features and chemical trajectory drift features; calculating the mutual information and phase synchronization coefficient of the two feature sequences to generate a phase space synchronization index; and using multi-source evidence fusion rules to determine the system state type. This method overcomes the problem of contradictory conclusions caused by independent diagnosis of physical and chemical fields, and has high diagnostic reliability.

Description

Technical Field

[0001] This invention relates to the fields of geophysics and hydrogeochemistry, and more specifically, to a method for the coordinated monitoring and resource utilization of medium-deep groundwater. Background Technology

[0002] In the long-term monitoring and resource development and utilization of medium-deep geothermal water, it is necessary to simultaneously acquire physical data (such as geothermal temperature, water pressure, and strain) and hydrochemical composition data of the coupled thermal-water-rock field to determine the operating status of the geothermal system. In complex areas such as the intersection of geothermal reservoirs and fault zones, the underground system exhibits highly nonlinear characteristics. The thermal conductivity of the fault zone exhibits a nonlinear response, causing weak precursor signals related to fault zone activity in the physical field evolution trajectory to be masked by normal fluctuation noise.

[0003] Current geothermal system condition diagnosis methods mainly fall into two categories: diagnosis based on physical field data and diagnosis based on chemical field data. Physical field diagnosis determines the system state by analyzing the changing trends of physical parameters such as geothermal temperature, water pressure, and strain. Chemical field diagnosis infers the dynamic changes in the water source mixing ratio by analyzing changes in water chemistry and isotopic composition. However, existing end-member mixing calculation methods can only provide the mixing ratio at the sampling time point and cannot track the continuous dynamic changes in chemical field characteristics.

[0004] The main drawback of existing technologies is that system state diagnosis based on physical and chemical field data can sometimes yield contradictory conclusions. This is because the two diagnostic methods operate independently, failing to consider the inherent correlation between the physical and chemical fields. When contradictory diagnostic results arise, it becomes impossible to determine the reliability of the results or whether the contradiction indicates a specific system anomaly. This deficiency in diagnostic methods limits the accuracy and reliability of geothermal system monitoring, becoming a key technical problem in assessing the state of hot water systems in complex areas. Summary of the Invention

[0005] This invention provides a collaborative method for monitoring and resource utilization of medium-deep groundwater, solving the technical problems of contradictory conclusions and insufficient diagnostic reliability caused by independent diagnosis of physical and chemical fields in related technologies.

[0006] This invention provides a method for the coordinated monitoring and resource utilization of medium-deep groundwater, comprising the following steps:

[0007] The physical field time series data and water chemical isotope data of the monitoring points are acquired. The low-frequency components related to the activity of the fault zone are separated from the physical field time series data using a bandpass filtering algorithm, and the fault zone response signal sequence and the original physical parameter sequence are generated.

[0008] Based on the original physical parameter sequence and the fault zone response signal sequence, an extended state space algorithm is used to construct an augmented coupled state vector, which is then mapped to a high-dimensional coupled state phase space to generate the physical field evolution trajectory.

[0009] The hydrochemical and isotopic data of each period are encoded into geochemical fingerprint vectors, projected onto a mixed space composed of end-member water source fingerprints, and the fingerprint drift trajectory curve is fitted using a trajectory smoothing algorithm to generate the chemical field evolution trajectory.

[0010] The local Lyapunov exponential sequence is calculated using a sliding window for the evolution trajectory of the physical field, and the motion direction and velocity characteristic sequences are calculated for the evolution trajectory of the chemical field, generating the stability characteristic sequence of the physical trajectory and the drift characteristic sequence of the chemical trajectory.

[0011] The physical trajectory stability feature sequence and the chemical trajectory drift feature sequence are time-aligned, the mutual information and phase synchronization coefficient of the two feature sequences at each time point are calculated, and the phase space synchronization index time series is generated by combining them.

[0012] Based on the range of phase space synchronicity index, physical trajectory stability feature sequence, and chemical trajectory drift feature sequence, the system state type is determined using multi-source evidence fusion rules, and collaborative diagnostic results and credibility assessment are output.

[0013] Furthermore, the frequency range of the bandpass filtering algorithm is determined based on the typical periodic characteristics of the fault zone activity, wherein the lower limit of the passband corresponds to the long-period activity characteristics of the fault zone, and the upper limit of the passband corresponds to the short-period response characteristics of the fault zone.

[0014] Furthermore, the dimension-enhanced coupled state vector is composed of the original physical parameter sequence and its multiple time-delay components and the fracture zone response signal sequence and its multiple time-delay components, wherein the time-delay parameter is determined based on the first zero-crossing point of the autocorrelation function.

[0015] Furthermore, when constructing the dimensionally enhanced coupled state vector, the extended state space algorithm assigns a higher weight coefficient to the fault zone response signal sequence than to the original physical parameter sequence.

[0016] Furthermore, the geochemical fingerprint vector includes the major anion-cation concentration ratio, stable isotope ratio, and trace element characteristic ratio.

[0017] Furthermore, the trajectory smoothing algorithm employs a combination of spline interpolation and Gaussian smoothing to generate smoothly transitioning trajectory segments between adjacent sampling points.

[0018] Furthermore, the phase space synchronization index is calculated as follows:

[0019] The mutual information between the two feature sequences is calculated based on the ratio of the joint probability distribution of the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence to their respective marginal probability distributions.

[0020] Simultaneously calculate the phase synchronization coefficient of the two feature sequences;

[0021] The phase space synchronization index is obtained by weighting and summing the normalized mutual information and the phase synchronization coefficient, with the sum of the weight coefficients being 1.

[0022] Furthermore, the system state types include synchronous stable state, synchronous changing state, physically dominated mutation state, chemically dominated drift state, and anomalous decoupling state;

[0023] Among them, when the phase space synchronization index is in the high value range and the local Lyapunov index and chemical trajectory drift rate are both in the low value range, it is determined to be a synchronous stable state.

[0024] When the phase space synchronicity index is in the high range and the local Lyapunov index and chemical trajectory drift rate increase synchronously, it is determined to be a synchronous change state.

[0025] When the phase space synchronicity index is in a low range and the local Lyapunov index is significantly increased while the chemical trajectory drift rate does not change significantly, it is judged to be a physically dominant mutation state.

[0026] When the phase space synchronicity index is in the low range and the chemical trajectory drift rate increases significantly while the local Lyapunov index does not change significantly, it is judged to be a chemically dominated drift state.

[0027] When the phase space synchronization index remains in a low range and the local Lyapunov index and chemical trajectory drift rate show an inverse relationship, it is determined to be an abnormal decoupling state.

[0028] Furthermore, it also includes the following steps:

[0029] A change point detection algorithm is applied to a local Lyapunov exponential sequence to identify potential precursor points of mutations and generate early warning markers for physical field mutations.

[0030] The rate of change of the projection components of the trajectory position to each endmember direction is calculated based on the chemical trajectory drift direction, and then converted into the rate of change of the contribution ratio of each endmember to generate a chemical field evolution trend marker.

[0031] By combining physical field mutation early warning markers and chemical field evolution trend markers with collaborative diagnostic results, mutation early warning levels are generated.

[0032] This invention provides a collaborative system for monitoring and resource utilization of medium-deep groundwater, used to execute the aforementioned method, including:

[0033] The signal separation module is used to acquire physical field time-series data and water chemical isotope data of the monitoring points, and uses a bandpass filtering algorithm to separate the fault zone response signal from the physical field time-series data, generating a fault zone response signal sequence and an original physical parameter sequence.

[0034] The physical field trajectory mapping module is used to construct an augmented-dimensional coupled state vector based on the original physical parameter sequence and the fault zone response signal sequence, and then map it to a high-dimensional coupled state phase space to generate the physical field evolution trajectory.

[0035] The chemical field trajectory mapping module is used to encode the hydrochemical and isotopic data of each period into geochemical fingerprint vectors, project them onto the mixed space composed of end-member water source fingerprints, and use the trajectory smoothing algorithm to fit the fingerprint drift trajectory curve to generate the chemical field evolution trajectory.

[0036] The feature extraction module is used to calculate the local Lyapunov exponential sequence for the physical field evolution trajectory, calculate the motion direction and velocity feature sequence for the chemical field evolution trajectory, and generate the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence.

[0037] The synchronization analysis module is used to time-align the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence, calculate the mutual information and phase synchronization coefficient, and generate a phase space synchronization index time series.

[0038] The collaborative diagnosis module is used to determine the system state type based on the phase space synchronization index, physical trajectory stability feature sequence, and chemical trajectory drift feature sequence, and to output collaborative diagnosis results and credibility assessment.

[0039] The beneficial effects of this invention are as follows:

[0040] The method for monitoring and resource utilization of medium-deep groundwater provided by this invention solves the problem of contradictory conclusions caused by independent diagnosis of physical and chemical fields by introducing a processing flow of fault zone response signal separation, two-phase spatial trajectory mapping and phase space synchronicity index fusion diagnosis, and achieves significant technical effects.

[0041] First, a bandpass filtering algorithm is used to separate the low-frequency components related to the activity of the fault zone from the physical field time series data, generate the fault zone response signal sequence, and construct an augmented coupled state vector together with the original physical parameter sequence. The weak precursor signals related to the activity of the fault zone in the physical field time series data can be effectively extracted from the background noise.

[0042] Secondly, the hydrochemical and isotopic data of each period are encoded into geochemical fingerprint vectors and projected onto a mixed spatial coordinate system composed of end-member water source geochemical fingerprints. The trajectory smoothing algorithm is used to fit the continuous chemical field evolution trajectory, and the dynamic change trend of the end-member mixing ratio can be continuously monitored.

[0043] Third, the mutual information and phase synchronization coefficients of the physical trajectory stability characteristic sequence and the chemical trajectory drift characteristic sequence are calculated, and the phase space synchronization index is generated in combination. The intrinsic correlation between the two independent diagnostic indicators can be quantitatively characterized, and the diagnostic results of the physical field and the chemical field can be mutually verified, thereby judging the reliability of the diagnostic results.

[0044] Therefore, through the synergistic effect of the above-mentioned technical means, the present invention provides dual evidence to support the diagnostic results, thus solving the problem of insufficient reliability in the status diagnosis of hot water systems in complex areas. Attached Figure Description

[0045] Figure 1 This is a flowchart of the method for coordinated monitoring and resource utilization of medium-deep groundwater in this invention;

[0046] Figure 2 This is a physical field time series data diagram of the present invention, showing the time series changes of ground temperature, water pressure and strain at the monitoring point within a 90-day period;

[0047] Figure 3 This is a graph showing the variation of geochemical fingerprint feature components, illustrating the main components in the geochemical fingerprint vector. and The changing trend;

[0048] Figure 4 This is a stacked area diagram of the changes in the water source contribution ratio of the present invention, showing the dynamic changes in the contribution ratio of three water sources (thermal storage water, shallow groundwater, and fault zone inflow water) during the monitoring period.

[0049] Figure 5 This is a line graph of the physical trajectory stability characteristics of the present invention, showing the change of the local Lyapunov exponent over time. This exponent is used to quantify the sensitivity of the physical field evolution trajectory to initial conditions.

[0050] Figure 6 This is a diagram showing the phase space synchronization index and diagnostic results of the present invention, illustrating the phase space synchronization index. The changes and threshold range of;

[0051] Figure 7 This is a heat map of the system status diagnosis process of the present invention. Detailed Implementation

[0052] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, some features described in the examples may be combined in other examples.

[0053] This embodiment provides a method for the coordinated monitoring and resource utilization of medium-deep groundwater, such as... Figure 1 As shown, it includes the following steps:

[0054] Step 1: Obtain physical field time series data and water chemical isotope data of the monitoring points, use bandpass filtering algorithm to separate the fault zone response signal, and generate the fault zone response signal sequence and the original physical parameter sequence.

[0055] In this step, time-series data of physical fields such as ground temperature, water pressure, and strain at monitoring points, as well as water chemistry and isotope analysis data from multiple sampling periods, are acquired. Since different physical quantities such as ground temperature, water pressure, and strain have different dimensions and numerical ranges, the time-series data of each physical field are first normalized to eliminate the influence of dimensional differences on subsequent calculations. A bandpass filter algorithm is then applied to the normalized physical field time-series data to separate the low-frequency components related to fault zone activity from the original data, outputting a fault zone response signal sequence. and the original physical parameter sequence .

[0056] It should be noted that the frequency range of the above bandpass filtering algorithm is determined based on the typical periodic characteristics of the fault zone activity, where the lower limit of the passband corresponds to the long-period activity characteristics of the fault zone, and the upper limit of the passband corresponds to the short-period response characteristics of the fault zone.

[0057] Step 2: Based on the original physical parameter sequence and the fault zone response signal sequence, construct an increased-dimensional coupled state vector using the extended state space algorithm, map it to a high-dimensional coupled state phase space, and generate the physical field evolution trajectory.

[0058] In this step, the original physical parameter sequence is... With the response signal sequence of the fault zone To achieve coupling, an extended state-space algorithm is used to construct an increased-dimensional coupled state vector. :

[0059]

[0060] in, The time delay parameter is determined based on the first zero crossing of the autocorrelation function. This represents the transpose operation of a vector.

[0061] Furthermore, the method for determining the number of delay terms in the dimension-enhanced coupled state vector is as follows: based on the embedding dimension criterion, the number of delay terms is determined... Set as ,in Let be the system dimension estimated by correlation dimension or other nonlinear dynamic methods. This setting ensures that the dimension-enhanced coupled state vector can fully represent the topological structure of the physical field evolution trajectory.

[0062] The dimension-enhanced coupled state vector is mapped to a high-dimensional coupled state phase space, and the state points at each time step are connected to form the physical field evolution trajectory. .

[0063] In this embodiment of the application, in order to enhance the sensitivity of the physical field evolution trajectory to the activity of the fault zone, the extended state space algorithm assigns a higher weight coefficient to the fault zone response signal sequence than the original physical parameter sequence when constructing the dimension-enhanced coupled state vector.

[0064] Furthermore, the weighting coefficients are assigned as follows: calculate the variance of the original physical parameter sequence and the fault zone response signal sequence, and set the weighting coefficient of the sequence with larger variance (usually the fault zone response signal sequence) to 1.5 to 2 times that of the original physical parameter sequence to enhance sensitivity.

[0065] Step 3: Encode the hydrochemical and isotopic data of each period into geochemical fingerprint vectors, project them onto a mixed spatial coordinate system composed of end-member water source geochemical fingerprints, and use a trajectory smoothing algorithm to fit the chemical field evolution trajectory.

[0066] In this step, the hydrochemical and isotopic data from each period are encoded into geochemical fingerprint vectors. ,in This represents the sampling period sequence number. Due to the significant differences in the numerical ranges of the ratios of various chemical indicators in the geochemical fingerprint vector, each ratio is normalized to eliminate the impact of these differences on the projection calculation. A hybrid spatial coordinate system is constructed based on the geochemical fingerprints of pre-calibrated end-source water sources (thermal reservoirs, shallow groundwater, fault zone inflow water, etc.). The normalized fingerprint vectors for each period are then... The projection is onto a hybrid spatial coordinate system composed of the geochemical fingerprints of the end-source water. A trajectory smoothing algorithm is applied to the projected points to fit the data, generating a continuous chemical field evolution trajectory. .

[0067] Furthermore, the method for projecting to the hybrid spatial coordinate system is as follows: The geochemical fingerprints of each end-source water source are used as basis vectors of the hybrid spatial coordinate system to construct a hyperpolyhedron with the geochemical fingerprints of the end-source water sources as vertices. For each period's fingerprint vector... Calculate the projected distance to each endmember direction, or use a hybrid proportional decomposition algorithm (such as non-negative least squares) to calculate. The mixing ratio on each endmember is used to determine the coordinates of the trajectory position in the mixed space coordinate system during that period.

[0068] It should be noted that the above geochemical fingerprint vector includes multidimensional chemical indicators such as the ratio of major anion and cation concentrations, the ratio of stable isotopes, and the characteristic ratio of trace elements.

[0069] Furthermore, the coding method for geochemical fingerprints is as follows: for sampling data from each period, the ratio of anion and cation concentrations is calculated (e.g., , (etc.), stable isotope ratios (such as...) , The ratios of trace element characteristics are combined and arranged in a fixed order to form a multidimensional vector. ,in The total number of chemical index ratios, each component ( () represents a ratio of a chemical index.

[0070] In this embodiment of the application, in order to overcome the problem of trajectory discontinuity caused by discrete sampling points, the trajectory smoothing algorithm adopts a processing method that combines spline interpolation and Gaussian smoothing to generate trajectory segments with smooth transitions between adjacent sampling points.

[0071] Step 4: Calculate the local Lyapunov exponent sequence for the physical field evolution trajectory using a sliding window, calculate the motion direction and velocity characteristic sequences for the chemical field evolution trajectory, and generate the physical trajectory stability characteristic sequence and the chemical trajectory drift characteristic sequence.

[0072] In this step, the evolution trajectory of the physical field is... The sliding window method is applied to calculate the local Lyapunov exponents within each window. This characterizes the stability of the physical field evolution during that period and generates a physical trajectory stability feature sequence. ,in This represents the time sequence corresponding to each sliding window.

[0073] Furthermore, the method for calculating the local Lyapunov exponent is as follows: select adjacent state points on the physical field evolution trajectory within the sliding window, calculate the rate of change of the distance between adjacent state points with time, take the logarithm of the rate of change and calculate the time average to obtain the local Lyapunov exponent at the corresponding moment of the window.

[0074] Furthermore, the size of the sliding window is set as: window width The number of monitoring sampling points is set to one-tenth to one-fifth to ensure that the window contains enough data points to obtain a stable estimate of the local Lyapunov exponent, while maintaining sensitivity to changes in the system state. The window slides point by point on the time axis, with each slide step being one sampling interval. .

[0075] Evolution trajectory of chemical field Calculate the motion direction vector of the trajectory point at each time step. and rate scalar :

[0076]

[0077] in, The norm (modulus) operation of a vector.

[0078] Generate chemical trajectory drift characteristic sequences:

[0079]

[0080] in This is a sequence of sampling times for the chemical field evolution trajectory.

[0081] Furthermore, the direction vector of motion The calculation method is as follows: for time... The trajectory point, calculate the distance between that point and the next moment. The position difference between trajectory points is normalized to obtain a motion direction vector with a magnitude of 1 pointing in the direction of trajectory motion. (Velocity scalar) By calculating the Euclidean distance between trajectory points at adjacent time points and dividing by the time interval... get.

[0082] It should be noted that the aforementioned local Lyapunov exponent is used to quantify the sensitivity of the phase space trajectory to small perturbations in initial conditions. Positive values ​​indicate that the system tends to be unstable, while negative values ​​indicate that the system tends to be stable.

[0083] Step 5: Time-align the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence, calculate the mutual information and phase synchronization coefficient of the two feature sequences at each time point, and generate a phase space synchronization index time series.

[0084] In this step, the physical trajectory stability feature sequence is analyzed. and chemical trajectory drift characteristic sequences Perform time alignment processing.

[0085] Furthermore, the time alignment method involves: [using the physical trajectory stability feature sequence...] At various times The corresponding local Lyapunov index With chemical trajectory drift characteristic sequence At various times corresponding rate Interpolation matching is performed, and the time axis is standardized so that the two sequences share a unified time coordinate. ( ).

[0086] For each aligned time point, calculate the mutual information of the two feature sequences. and phase synchronization coefficient .

[0087] The mutual information is calculated based on the ratio of the joint probability distribution of the two feature sequences to their respective marginal probability distributions, taking into account the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence.

[0088]

[0089] in, Represents the logarithmic function. For joint probability distribution, and These are marginal probability distributions, and the output is a sequence of mutual information. .

[0090] Furthermore, the joint probability distribution and marginal probability distribution are calculated as follows: the value ranges of the two feature sequences are divided into several discrete intervals of equal width, and the frequency of the value falling into each interval at each time point is counted. The probability is obtained by dividing the frequency by the total number of time points. The joint probability is the frequency of both sequence values ​​falling into the corresponding interval simultaneously divided by the total number of time points, and the marginal probability is the frequency of a single sequence value falling into the corresponding interval divided by the total number of time points.

[0091] Furthermore, the method for determining the number of discretization intervals is as follows: based on the Sturges formula or the equal bandwidth method, the value range of each feature sequence is uniformly divided into... There are several intervals, among which , This indicates the rounding up operation. Represents the logarithmic function with base 2. This represents the total number of sampling points in the feature sequence. This division method ensures statistical stability while avoiding over-refinement or over-coarseness.

[0092] Furthermore, in the formula and The definition of is: Physical trajectory stability feature sequence In the discretized interval index Chemical trajectory drift characteristic sequence rate component After discretization, the summation is performed on all possible interval index pairs. conduct.

[0093] To calculate the phase synchronization coefficient, the instantaneous phase is first extracted from the physical trajectory stability characteristic sequence and the chemical trajectory drift characteristic sequence, and then the phase difference between the two sequences at each time point is calculated. The phase synchronization coefficient is generated by measuring the stability of the phase difference. The phase synchronization coefficient is used to quantify the strength of the phase coupling between two characteristic sequences, and its value ranges from 0 to 1.

[0094] Furthermore, the instantaneous phase extraction method involves applying a Hilbert transform to the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence to obtain their respective analytical representations, from which the instantaneous phase is extracted. For analytical signals... ,in The original signal sequence to be analyzed. Represents the imaginary unit, and its instantaneous phase is ,in Represents the arctangent function. This indicates taking the imaginary part of a complex number. To take the real part of a complex number, for The Hilbert transform.

[0095] Furthermore, the stability measurement method for phase difference is as follows: calculate the standard deviation of phase difference within the sliding time window. The smaller the standard deviation, the more stable the phase difference and the closer the phase synchronization coefficient is to 1; the larger the standard deviation, the more unstable the phase difference and the closer the phase synchronization coefficient is to 0.

[0096] Furthermore, the sliding time window size used to calculate the standard deviation of the phase difference is: window width Set as This is done in accordance with the window used for calculating the local Lyapunov index, to ensure that the two sequences are measured for synchronicity on the same time scale.

[0097] Furthermore, the phase synchronization coefficient The calculation formula is: ,in Represents an exponential function. For time points The standard deviation of the phase difference within the corresponding sliding window, This is the normalization constant, typically taking the value of [value missing]. To ensure The value ranges from 0 to 1.

[0098] Since the mutual information value ranges to positive values ​​but its upper limit is not fixed, while the phase synchronization coefficient value ranges to a fixed range between 0 and 1, in order to perform weighted summation, it is necessary to adjust the mutual information sequence. Normalization is performed to adjust the value range to between 0 and 1. Then, a weighted summation is performed to generate the phase space synchronization index time series. :

[0099]

[0100] in, To normalize the mutual information content, This is the phase synchronization coefficient. and The weighting coefficients and , and The values ​​of are all between 0 and 1.

[0101] Furthermore, weighting coefficients and The method for determining it is as follows: calculate the normalized mutual information sequence. and phase synchronization coefficient sequence The variances of each sequence are weighted according to the inverse of the variance, that is, the sequence with smaller variance is assigned a larger weight, so as to improve the contribution of stable features to the phase space synchronization index.

[0102] Furthermore, weighting coefficients and It remains constant throughout the entire monitoring period, and its value range is constrained as follows: , and The specific calculation formula is as follows: , ,in and These are the variances of the normalized mutual information sequence and the phase synchronization coefficient sequence over the entire monitoring period, respectively.

[0103] Step 6: Using multi-source evidence fusion rules, determine the system state type based on the value range of the phase space synchronization index and the stability characteristics of the physical trajectory, and output collaborative diagnosis results and credibility assessment.

[0104] In this step, based on the phase space synchronization index The range of values ​​and the physical trajectory stability characteristic sequence Local Lyapunov index variation trend and chemical trajectory drift characteristic sequence The drift direction and rate are determined, and the system state is judged using multi-source evidence fusion rules.

[0105] Furthermore, the high-value and low-value ranges of the phase space synchronization index are defined as follows: the high-value range is... The low value range is

[0106] High value threshold Low threshold ,in and These represent the mean and standard deviation of the phase-space synchronicity index sequence over the entire monitoring period. The low range of the Lyapunov exponent in the physical trajectory stability characteristic sequence is defined as... , where the threshold The mean of the Lyapunov index sequence over the entire monitoring period. The low range of rates in the chemical trajectory drift characteristic sequence is defined as... , where the threshold This represents the mean of the rate sequence over the entire monitoring period.

[0107] Furthermore, the system state judgment rule is as follows:

[0108] when and and When this occurs, it is determined to be a synchronous stable state;

[0109] when and and When all values ​​exceed the mean in their respective historical sequences, they are considered to be in a state of synchronous change.

[0110] when and Exceeding one standard deviation of its mean Still below When this occurs, it is determined to be a physically dominant mutation state;

[0111] when and Exceeding one standard deviation of its mean Still below At that time, it was determined to be a chemically dominated drift state;

[0112] when The occurrence of consecutive occurrences exceeding the preset time threshold And during this period and When the changing trends show an inverse relationship, it is judged as an abnormal decoupling state.

[0113] Furthermore, a preset time threshold The method for determining this is as follows: set it as the monitoring sampling interval. Several times, usually taken That is, the continuous low value interval must contain at least 10 sampling periods. The criterion for determining the inverse relationship is: within this time period... The linear trend is monotonically increasing. The linear trend is either monotonically decreasing or vice versa.

[0114] The collaborative diagnostic results are output based on the judgment results. The reliability assessment method is as follows: the longer the duration and the larger the amplitude of the high value range of the phase space synchronicity index, the higher the reliability of the diagnostic results; conversely, the lower the duration and the larger the amplitude of the high value, the higher the reliability of the diagnostic results. When the value fluctuates near the threshold boundary or is short in duration, the reliability of the diagnostic results decreases accordingly.

[0115] Furthermore, the quantitative assessment of credibility is as follows: Let the total duration for which the phase space synchronicity index remains within the high-value range be... Calculate the time period average Then credibility Defined as:

[0116]

[0117] in This represents the duration of the entire monitoring period. The reliability value ranges from 0 to 1, with a value closer to 1 indicating higher reliability of the diagnostic result.

[0118] In addition to the steps described above, the following steps are also included:

[0119] Step 7: Apply the change point detection algorithm to the local Lyapunov index sequence to identify potential precursor points of mutations. At the same time, based on the chemical trajectory drift direction, convert the trajectory position into the rate of change of the contribution ratio of each endmember, and generate physical field mutation early warning markers and chemical field evolution trend markers.

[0120] In this step, the local Lyapunov exponential sequence in the physical trajectory stability feature sequence is analyzed. By applying a change point detection algorithm, potential precursor points of mutations in the sequence are identified, and physical field mutation early warning markers are generated.

[0121] Furthermore, the criteria for identifying potential precursor points of mutation are as follows: if the value of the local Lyapunov exponent sequence jumps abnormally at a certain moment, or if it shows a monotonically increasing trend over a period of time and the rate of change exceeds a preset threshold, then that moment or period of time is identified as a potential precursor point of mutation.

[0122] Furthermore, the criterion for identifying anomalous jumps is the change in the local Lyapunov exponent between adjacent time points. More than twice the standard deviation of the sequence, i.e. The preset threshold for the rate of change is set as follows: within a time window Within the region, the linear growth slope of the Lyapunov index exceeds... ,in For reference time length, it is usually taken as one-tenth of the monitoring cycle.

[0123] Meanwhile, based on chemical trajectory drift characteristic sequences The drift direction vector in The rate of change of the projection components of the trajectory position to each endmember direction is calculated, and then converted into the rate of change of the contribution ratio of each endmember to generate a chemical field evolution trend marker.

[0124] Furthermore, the conversion method for the rate of change of endmember contribution ratio is as follows: In the hybrid spatial coordinate system, the rate of change of distance between the current trajectory position and the position of each endmember is taken as the rate of change of the endmember contribution ratio. The greater the rate of decrease in distance, the faster the drift rate towards the endmember, and the higher the rate of increase in the contribution ratio of the endmember.

[0125] Physical field mutation warning markers and chemical field evolution trend markers are combined with the collaborative diagnostic results from step 6 as auxiliary information to generate mutation warning levels.

[0126] It is understood that data preprocessing methods known to those skilled in the art include data cleaning, data transformation, and data reduction. Data transformation includes type conversion and normalization and standardization. Although the dimensions and types of data were omitted in the description of the preceding embodiments, data preprocessing is a technical knowledge known to those skilled in the art and a prerequisite step in data processing. Therefore, the previously described well-known data preprocessing steps were not described independently.

[0127] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A method for the coordinated monitoring and resource utilization of medium-deep groundwater, characterized in that, Includes the following steps: The physical field time series data and water chemical isotope data of the monitoring points are acquired. The low-frequency components related to the activity of the fault zone are separated from the physical field time series data using a bandpass filtering algorithm, and the fault zone response signal sequence and the original physical parameter sequence are generated. Based on the original physical parameter sequence and the fault zone response signal sequence, an extended state space algorithm is used to construct an augmented coupled state vector, which is then mapped to a high-dimensional coupled state phase space to generate the physical field evolution trajectory. The hydrochemical and isotopic data of each period are encoded into geochemical fingerprint vectors, projected onto a mixed space composed of end-member water source fingerprints, and the fingerprint drift trajectory curve is fitted using a trajectory smoothing algorithm to generate the chemical field evolution trajectory. The local Lyapunov exponential sequence is calculated using a sliding window for the evolution trajectory of the physical field, and the motion direction and velocity characteristic sequences are calculated for the evolution trajectory of the chemical field, generating the stability characteristic sequence of the physical trajectory and the drift characteristic sequence of the chemical trajectory. The physical trajectory stability feature sequence and the chemical trajectory drift feature sequence are time-aligned. Based on the ratio of the joint probability distribution of the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence to their respective marginal probability distributions, the mutual information between the two feature sequences is calculated. Simultaneously calculate the phase synchronization coefficient of the two feature sequences; The normalized mutual information and the phase synchronization coefficient are weighted and summed, and the sum of the weight coefficients is 1, to generate the phase space synchronization index time series. Based on the range of values ​​of the phase space synchronization index, the physical trajectory stability feature sequence, and the chemical trajectory drift feature sequence, the system state type is determined by multi-source evidence fusion rules, and collaborative diagnostic results and credibility assessment are output. The high-value range of the phase space synchronization index is defined as follows: The low value range is defined as Among them, the high value threshold Low threshold , and These are the mean and standard deviation of the phase space synchronicity index sequence over the entire monitoring period, respectively. The low range of the Lyapunov exponent in the physical trajectory stability feature sequence is defined as follows: , where the threshold The mean of the Lyapunov index series over the entire monitoring period; The low range of rates in the chemical trajectory drift characteristic sequence is defined as... , where the threshold This represents the mean of the rate sequence over the entire monitoring period. The system state types include synchronous stable state, synchronous changing state, physically dominated mutation state, chemically dominated drift state, and anomalous decoupling state. When the phase space synchronization index is in the high range and the local Lyapunov index and chemical trajectory drift rate are both in the low range, it is determined to be a synchronous stable state. When the phase space synchronicity index is in the high range and the local Lyapunov index and chemical trajectory drift rate increase synchronously, it is determined to be a synchronous change state. When the phase space synchronicity index is in a low range and the local Lyapunov index is significantly increased while the chemical trajectory drift rate does not change significantly, it is judged to be a physically dominant mutation state. When the phase space synchronicity index is in the low range and the chemical trajectory drift rate increases significantly while the local Lyapunov index does not change significantly, it is judged to be a chemically dominated drift state. When the phase space synchronization index remains in a low range and the local Lyapunov index and chemical trajectory drift rate show an inverse relationship, it is determined to be an abnormal decoupling state.

2. The method according to claim 1, characterized in that, The frequency range of the bandpass filtering algorithm is determined based on the typical periodic characteristics of the fault zone activity, wherein the lower limit of the passband corresponds to the long-period activity characteristics of the fault zone, and the upper limit of the passband corresponds to the short-period response characteristics of the fault zone.

3. The method according to claim 1, characterized in that, The dimension-enhanced coupled state vector is composed of the original physical parameter sequence and its multiple time-delay components, and the fault zone response signal sequence and its multiple time-delay components. The time-delay parameters are determined based on the first zero-crossing point of the autocorrelation function.

4. The method according to claim 3, characterized in that, The extended state-space algorithm assigns a higher weighting coefficient to the fault zone response signal sequence than to the original physical parameter sequence when constructing the dimensionally increased coupled state vector.

5. The method according to claim 1, characterized in that, The geochemical fingerprint vector includes the major anion-cation concentration ratio, stable isotope ratio, and trace element characteristic ratio.

6. The method according to claim 1, characterized in that, The trajectory smoothing algorithm uses a combination of spline interpolation and Gaussian smoothing to generate trajectory segments with smooth transitions between adjacent sampling points.

7. The method according to claim 1, characterized in that, It also includes the following steps: A change point detection algorithm is applied to the local Lyapunov index sequence to identify potential precursor points of mutations and generate physical field mutation early warning markers. The rate of change of the projection components of the trajectory position to each endmember direction is calculated based on the chemical trajectory drift direction and converted into the rate of change of the contribution ratio of each endmember to generate a chemical field evolution trend marker. The physical field mutation early warning markers and chemical field evolution trend markers are combined with the collaborative diagnostic results to generate a mutation early warning level.

8. A collaborative system for monitoring and resource utilization of medium-deep groundwater, used to execute the method described in any one of claims 1-7, characterized in that, include: The signal separation module acquires physical field time-series data and hydrochemical isotope data from monitoring points, and uses a bandpass filtering algorithm to separate the fault zone response signal from the physical field time-series data, generating a fault zone response signal sequence and an original physical parameter sequence. The physical field trajectory mapping module constructs an augmented-dimensional coupled state vector based on the original physical parameter sequence and the fault zone response signal sequence using an extended state space algorithm, and maps it to a high-dimensional coupled state phase space to generate the physical field evolution trajectory. The chemical field trajectory mapping module encodes hydrochemical and isotope data from each period into geochemical fingerprint vectors, projects them onto a mixed space composed of end-member water source fingerprints, and uses a trajectory smoothing algorithm to fit the fingerprint drift trajectory curve to generate the chemical field evolution trajectory. The feature extraction module calculates the local Lyapunov exponent sequence for the physical field evolution trajectory and the motion direction and velocity feature sequences for the chemical field evolution trajectory, generating a physical trajectory stability feature sequence and a chemical trajectory drift feature sequence. The synchronization analysis module is used to time-align the physical trajectory stability feature sequence and the chemical trajectory drift feature sequence, calculate the mutual information and phase synchronization coefficient, and generate a phase space synchronization index time series. The collaborative diagnosis module is used to determine the system state type based on the phase space synchronization index, physical trajectory stability feature sequence, and chemical trajectory drift feature sequence, and to output collaborative diagnosis results and credibility assessment.

Images