Seismic underground fluid observation data interference identification method and system based on time-frequency response characteristics

Abstract

The application belongs to the technical field of data detection and processing, and specifically discloses a seismic underground fluid observation data interference identification method and system based on time-frequency response characteristics, which comprises the following steps: analyzing the time-frequency and space-time coupling relationship of multi-point observation data, and constructing a standardized overall measured response mode; establishing a forward physical simulation model composed of a parameterized geological structure network and a non-steady flow control equation; under an iterative optimization framework, the physical parameters of the simulation model are adjusted coordinately to drive the generated theoretical overall response mode to approximate the measured overall response mode; finally, the optimal solution is determined through a comprehensive scoring function taking into account the matching residual and the parameter adjustment amplitude. By introducing the cost function of the physical parameter adjustment as a regularization term, the solution space is constrained, so that the physical interpretability and convergence stability of the inversion result are improved while the matching accuracy is ensured.

Description

Technical Field

[0001] This invention belongs to the field of data detection and processing technology, and relates to a method and system for identifying interference in seismic underground fluid observation data based on time-frequency response characteristics. Background Technology

[0002] Subsurface fluid observation, as a crucial method for earthquake precursor monitoring, continuously monitors parameters such as well water level, water temperature, and hydrochemical composition to capture anomalous responses of subsurface fluid systems that may be caused by changes in the crustal stress field. In addition to information related to tectonic activity, the observation data typically contains responses caused by non-tectonic factors such as rainfall, air pressure changes, anthropogenic water injection, and solid tides. Identifying and separating these non-tectonic interferences from the observation data is a necessary prerequisite for extracting information related to tectonic activity.

[0003] In existing technical practices, methods for handling interference in underground fluid observation data mainly include threshold alarm methods, statistical correlation analysis methods, and conventional filtering methods. Threshold alarm methods determine anomalies by setting a fixed limit for data variation; an alarm is triggered when data fluctuations exceed this limit. Statistical correlation analysis calculates the temporal correlation coefficient between the observed data sequence and a known interference source sequence, such as rainfall or air pressure sequences, to determine the strength of the correlation. Conventional filtering methods, such as Fourier transforms or simple bandpass / bandstop filters, attempt to separate interference signals in specific frequency bands in the frequency domain.

[0004] However, the aforementioned existing technologies have significant technical shortcomings. Threshold-based methods rely solely on the magnitude of data changes, making it difficult to distinguish between changes of equal magnitude caused by different physical sources. Furthermore, fixed thresholds have limited adaptability in dynamically changing environments. Statistical correlation-based analysis methods depend on the knowledge of the interference source and the availability of synchronous observation data, and are primarily suitable for handling linear response relationships. Frequency-domain separation-based filtering methods assume that the target signal and interference signal are separable in the frequency domain; however, when the frequency bands of earthquake precursor signals and interference signals overlap, the filtering operation may distort the target signal while suppressing interference. Therefore, a more accurate technical solution for identifying and processing complex interference is still needed in this field. Summary of the Invention

[0005] To overcome the above-mentioned defects of the prior art and to achieve the above objectives, the present invention proposes the following technical solution: a method for identifying interference in seismic underground fluid observation data based on time-frequency response characteristics, comprising: S1, acquiring multi-source geological information and observation network topology of the target area, and constructing a parameterized geological structure network model with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges, wherein the parameters of the connecting edges have an initial parameter range.

[0006] S2. Monitor earthquake underground fluid observation data. When anomalies in the observation data are identified, combine the data from the external multi-source information database to generate at least one candidate physical interference scenario that includes the type, intensity, and spatiotemporal location of the interference.

[0007] S3. Using the parameters of the candidate physical disturbance scenarios as input, drive the parameterized geological structure network model to perform parallel simulation, and generate a theoretical overall response mode for each candidate physical disturbance scenario.

[0008] S4. Based on the measured data of all relevant observation points during the period of abnormal observation data, construct the overall measured response model.

[0009] S5. Match the measured overall response mode with all theoretical overall response modes, and through a feedback mechanism, with the goal of minimizing the difference between the measured overall response mode and the theoretical overall response mode, coordinately optimize the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model to obtain the matching judgment result.

[0010] S6. Based on the interference scenario finally determined in the matching judgment result, decouple the interference component corresponding to the scenario from the original observation data and output the purified fluid observation data.

[0011] The second aspect of the present invention provides a seismic underground fluid observation data interference identification system based on time-frequency response characteristics, comprising: a model building module for acquiring multi-source geological information and observation network topology of a target area, and constructing a parameterized geological structure network model with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges, wherein the parameters of the connecting edges have an initial parameter range.

[0012] The scene generation module is used to monitor earthquake underground fluid observation data. When anomalies in the observation data are identified, it combines data from an external multi-source information database to generate at least one candidate physical interference scene containing the interference type, intensity, and spatiotemporal location.

[0013] The simulation and deduction module is used to take the parameters of the candidate physical disturbance scenarios as input, drive the parameterized geological structure network model to perform parallel simulation, and generate a theoretical overall response mode for each candidate physical disturbance scenario.

[0014] The pattern building module is used to construct the overall measured response pattern based on the measured data of all relevant observation points during the period of observational data anomalies.

[0015] The matching and determination module is used to match the measured overall response mode with all theoretical overall response modes. Through a feedback mechanism, with the goal of minimizing the difference between the measured overall response mode and the theoretical overall response mode, the module performs collaborative optimization of the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model to obtain the matching and determination results.

[0016] The data decoupling module is used to decouple the interference components corresponding to the interference scenario determined in the matching judgment result from the original observation data and output the purified fluid observation data.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: (1) The present invention transforms the interference identification problem into a parameter inversion problem under the constraints of a physical model by constructing a parameterized geological structure network model and a set of groundwater unsteady flow control equations. Its output results are not only the classification of interference, but also a set of model parameters with clear physical meaning that can explain the observed phenomena. By introducing a cost function for adjusting the magnitude of physical parameters, the understanding space is further constrained, thereby reducing the ambiguity of the identification results and improving the stability of the results when dealing with complex or atypical interference events.

[0018] (2) This invention can reconstruct the theoretical response of the interference event at each observation point by reversing the optimal model parameters. Subtracting this theoretical response from the original observation data can separate the interference signal. The purified data obtained by this method has a higher signal-to-noise ratio than the original data, which provides a technical basis for subsequent analysis and extraction of target signals with small amplitudes that may be contained in the data.

[0019] (3) This invention stores each successfully identified interference scenario and its corresponding optimal model parameters into a case library. This case library can provide optimized initial parameters for subsequent processing of similar interference events, or serve as prior knowledge to constrain new inversion processes. This mechanism allows the system's model parameters to be continuously optimized as the number of application cases increases, thereby improving the computational efficiency and convergence stability for processing similar events. Attached Figure Description

[0020] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a schematic diagram of the implementation steps of the method of the present invention.

[0022] Figure 2This is a schematic diagram of the system module connections of the present invention. Detailed Implementation

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] Example 1

[0025] Please see Figure 1 As shown, the earthquake underground fluid observation data interference identification method based on time-frequency response characteristics proposed in this invention includes: S1, acquiring multi-source geological information and observation network topology of the target area, and constructing a parameterized geological structure network model with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges, wherein the parameters of the connecting edges have an initial parameter range.

[0026] In a preferred embodiment, a parameterized geological structure network model is constructed with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges. This includes: mapping observation wells and potential interference sources as network nodes with spatial attributes based on regional hydrogeological data, and establishing directed connecting edges representing groundwater flow channels based on aquifer distribution and fault structures, so as to construct the network topology of the geological structure network model. Each connecting edge is assigned physical parameters that characterize the hydrodynamic properties of the medium. These physical parameters include at least the permeability coefficient and the water storage coefficient. Based on regional stratigraphic lithology data or historical pumping test data, set initial parameter ranges for physical parameters; Based on a parameterized geological structure network model and combined with regional hydrological boundary conditions, a set of groundwater unsteady flow control equations is established to describe the dynamic response of the network model, thus forming a forward physical simulation model.

[0027] Specifically, this embodiment discloses a detailed engineering method for constructing a parameterized geological structure network model. This step aims to abstract the complex actual geological environment into a computable and optimizable mathematical graph model, providing a physical foundation for subsequent simulation of disturbance scenarios.

[0028] The first step is the digital mapping from physical geological entities to network nodes: The system first accesses multi-source geological information of the target area, which mainly includes regional geological structure maps, aquifer lithology distribution data, major fault zone distribution data, and borehole columnar sections of observation wells. The system utilizes Geographic Information System (GIS) technology to extract key spatial entities and abstract them into a set of network nodes, V. Nodes are strictly divided into two categories: the first category is "observation nodes," which directly correspond to existing seismic underground fluid observation wells, and their attributes include well depth and observed aquifer strata; the second category is "potential source nodes," which correspond to potential interference source locations or hydraulic boundaries, including major river sections, reservoirs, the center of known industrial pumping well groups, or the intersection of major water-conducting fault zones.

[0029] The second step involves constructing a logical topology based on hydrogeological patterns: analyzing the spatial location of each node and the attributes of its hydrogeological unit. If two nodes are located within the same connected aquifer, and the geological map shows no obstructing faults or impermeable dikes between them, a connecting edge E is established between the two nodes. This connecting edge physically represents the pressure conduction path or material transport channel of underground fluids. For areas with obvious regional groundwater flow directions, such as alluvial fans at the foot of mountains, directed edges are established; for still water environments in plains or reciprocating flow environments affected by tides, bidirectional edges are established. This constructs the basic topology diagram. This topology defines the possible paths for subsequent interference propagation.

[0030] The third step is to assign initial parameter ranges to the model's connecting edges: Based on the lithological descriptions in geological exploration data, including coarse sand, fractured limestone, and clay, the system consults standard hydrogeological parameter manuals to assign initial parameter ranges for core hydraulic parameters to each connecting edge. The main parameters include the equivalent permeability coefficient K and the water storage coefficient S. The parameter ranges for connecting edge (i,j) are defined as follows:

[0031]

[0032] in, This represents the edge connecting node i and node j in the connection model. The initial parameter range represents the equivalent permeability coefficient K of the connecting edge (i,j). The initial parameter range represents the water storage coefficient S of the connecting edge (i,j). and These are the minimum and maximum allowable values ​​for the equivalent permeability coefficient K, respectively. and These are the minimum and maximum allowable values ​​for the water storage coefficient S, respectively.

[0033] For example, when connecting two nodes in a "broken limestone" aquifer, the system determines the initial parameter range of its permeability coefficient based on prior knowledge. Set as m / s. This range is wider than specific single-point measurements to accommodate the uncertainty and heterogeneity of geological structures; however, it is also strictly limited by physical lithology to prevent parameter drift that violates physical principles during subsequent optimization, such as identifying dense clay layers as highly conductive channels. The resulting parameterized geological structure network model is mathematically represented as a weighted graph with parameter range constraints, which forms the physical framework for all subsequent simulations.

[0034] S2. Monitor earthquake underground fluid observation data. When anomalies in the observation data are identified, combine the data from the external multi-source information database to generate at least one candidate physical interference scenario that includes the type, intensity, and spatiotemporal location of the interference.

[0035] In a preferred embodiment, generating at least one candidate physical interference scenario that includes interference type, intensity, and spatiotemporal location includes: extracting preliminary morphological features of anomalies in observation data to form anomaly feature descriptors; Call a cause-effect knowledge template library that defines the mapping relationship between physical interference causes and data anomaly patterns, match the anomaly feature descriptors with the library, and obtain a list of potential physical interference causes; Using a list of potential physical interference causes as an index, external multi-source information databases are retrieved, and the matched external information events are instantiated into parameterized scene descriptions to form candidate physical interference scenarios. The data in the external multi-source information database includes at least two of the following: meteorological precipitation data, reservoir scheduling records, logs of manual water pumping activities, and information on surrounding earthquake sources.

[0036] Specifically, this embodiment details an engineering method for automatically and intelligently generating a series of physically plausible candidate physical interference scenarios when anomalies occur in observation data. The goal of this method is to transform a vague observation data anomaly into a set of explicit, parameterized physical event hypotheses that can be used for subsequent simulation and deduction, thereby providing initial input for interference identification.

[0037] The first step is to perform preliminary morphological classification of anomalies in the original observation data: The system maintains a dynamic baseline (e.g., a moving average) and its standard deviation for the time series data of each observation point based on a sliding time window (e.g., the past 24 hours). When the observation value at a certain moment deviates from this dynamic baseline by more than a preset threshold, such as 3 times the standard deviation, that point is determined to be the starting point of the anomaly, constituting an observation data anomaly.

[0038] The system then extracts the abnormal data segment from the anomaly initiation point to the point where the data regresses to the threshold, and calculates its preliminary morphological characteristics based on the data segment to form an anomaly feature descriptor. The calculation methods for these features include: Abnormal amplitude: Calculates the difference between the peak or valley value in the abnormal data segment and the dynamic baseline at the point of origin of the abnormality.

[0039] Anomaly duration: Calculates the time from the point of origin of the anomaly to the point where the data first stabilizes and returns to within the threshold.

[0040] Abnormal trend types: Classified by analyzing the characteristics of the first derivative (rate of change) of abnormal data segments. For example, a derivative segment that is consistently positive and has a large rate of change corresponds to a "step increase"; a sharp derivative pulse that is first positive and then negative corresponds to a "positive pulse"; and a derivative with periodic alternation of positive and negative values ​​corresponds to "periodic oscillation".

[0041] The second step is to establish a correlation between abstract data forms and specific physical causes: The system maintains a pre-defined cause-effect knowledge template library, which stores the mapping relationship between different physical disturbance causes and their typical data anomaly forms in a structured form. This template library is built based on historical data analysis and hydrogeological expert knowledge. For example, the template defines that "short-term heavy rainfall" usually corresponds to "gentle slope rise with a delay of several hours," while "nearby industrial pumping" corresponds to "periodic step-like descent related to working hours." To more accurately describe the mapping relationship, this embodiment provides the following example table of the data structure of the cause-effect knowledge template library:

[0042] The system will use the anomaly feature descriptor obtained in the previous step As input, a matching search is performed in the cause-effect knowledge template base to find all instances where the "effect" part matches... "Reason" entries with a similarity higher than a preset threshold (e.g., 0.7).

[0043] The third step is to instantiate the ambiguous "cause" entries into specific scenario descriptions: The system uses the matched "cause" entries as search indexes to automatically query the accessed external multi-source information databases. For example, if the matched cause is "regional rainfall," the system will query the meteorological database for rainfall records around the observation point in the period before the anomaly occurred. If a rainfall event that meets the criteria is found, the system will combine the actual parameters of the rainfall event, including rainfall amount, duration, and coverage, with the cause entry "regional rainfall" to generate a specific, parameterized scenario description. For example, "A rainfall event with an average intensity of 10 mm / h and a duration of 3 hours occurred in a range of 5-10 kilometers northeast of the observation point at a certain time." Each such specific scenario description is a candidate physical disturbance scenario. The system will repeat this process for all matched "cause" entries, eventually generating a set of candidate physical disturbance scenarios containing disturbance type, intensity, spatial location, and temporal parameters, and output it to the subsequent parallel simulation module.

[0044] Furthermore, this embodiment details the specific composition and acquisition methods of the external multi-source information upon which the system relies to generate candidate physical interference scenarios. The engineering purpose of this step is to provide diverse and reliable sources of background information for interference identification, ensuring that the system can consider as many potential physical causes as possible when making hypothesis reasoning.

[0045] The system uses standardized data interfaces to access and process at least two types of external multi-source information in real time or near real time. The selection of these information sources is based on the physical understanding of common influencing factors of underground fluid changes.

[0046] The first category is meteorological precipitation data. The system periodically requests data such as rainfall, snowfall, and air pressure within the coverage area of ​​the observation network from the regional meteorological data center. In engineering practice, this data is usually provided in a gridded format, with a spatial resolution of 1-5 kilometers and a temporal resolution of hours or minutes. After receiving this data, the system performs spatial interpolation to obtain the time series of meteorological parameters for each observation well location and potential recharge area. This data is a key basis for determining fluid changes caused by rainfall infiltration or air pressure effects.

[0047] The second category is data on water conservancy project activities. The system connects to the database of the regional water management department to obtain real-time water levels, gate openings, and scheduling plans for large reservoirs or rivers—in other words, reservoir scheduling records. For example, the system pays particular attention to the flood discharge or water storage operations of upstream reservoirs, which, through the hydraulic connection between surface water and groundwater, can affect the water levels in observation wells within hours to days. This data is crucial for identifying disturbances related to the activities of large surface water bodies.

[0048] The third category is data on human production activities, specifically logs of manual water pumping and injection activities. The system establishes connections with the data systems of surrounding industrial parks, mines, and farmland irrigation management areas to obtain log data such as pumping and injection start / stop times and flow rates from their production or irrigation wells. These activities typically cause rapid and significant changes in groundwater levels near the well points, making them common sources of strong interference. Obtaining these logs allows the system to accurately match the periodic fluctuations of observed data with the production activity cycle.

[0049] The fourth category is regional seismic geological activity information, namely, information on surrounding seismic sources. The system connects to the seismic network center to acquire real-time catalogs of earthquakes occurring in the surrounding area, typically within a range of hundreds of kilometers, including magnitude, epicenter location, focal depth, and time of occurrence. Seismic waves from far-field earthquakes, especially surface waves, can cause instantaneous changes in aquifer permeability, leading to coseismic responses or post-seismic effects in well water levels. This is an important type of non-artificial interference that needs to be distinguished from precursor anomalies.

[0050] The system unifies these heterogeneous external multi-source information by timestamp alignment and formatting, and stores them in an internal database that can be quickly retrieved, providing timely and accurate environmental background data support for subsequent scene generation modules.

[0051] S3. Using the parameters of the candidate physical disturbance scenarios as input, drive the parameterized geological structure network model to perform parallel simulation, and generate a theoretical overall response mode for each candidate physical disturbance scenario.

[0052] In a preferred embodiment, generating a theoretical overall response mode for each candidate physical disturbance scenario includes: converting the parameters of the candidate physical disturbance scenario into source and sink terms or boundary conditions of a parameterized geological structure network model. The theoretical time-series response data of all relevant observation points were calculated by numerically solving the governing equations of the parameterized geological structure network model. Joint time-frequency analysis and spatial correlation analysis are performed on the theoretical time-series response data to synthesize the overall theoretical response mode.

[0053] Specifically, this embodiment discloses detailed engineering steps for generating the corresponding theoretical overall response mode based on a given candidate physical disturbance scenario through forward simulation using a physical model. The goal of this method is to transform a text and parameters describing a physical event into a high-dimensional theoretical data volume that is completely consistent with the measured data structure and can be directly compared.

[0054] The first step involves transforming discrete scene parameters into continuous boundary conditions solvable by the physical model: The system receives a candidate physical disturbance scenario generated in the previous step, such as "pumping water at coordinates (X,Y) for a period of time T at a flow rate Q at a certain time." Simultaneously, the system invokes the constructed parameterized geological structure network model. Mathematically, this model is a system of partial differential equations describing groundwater flow, typically the Boussinesq equations for unsteady flow in unconfined aquifers. The system transforms the parameters in the candidate physical disturbance scenario, such as the flow rate Q and the duration t, into source and sink terms or second-type boundary conditions of the equation system, applying them to the corresponding nodes or regions in the network model. The governing equations are as follows:

[0055] in, It is the water head height, that is, the elevation of the groundwater surface, which can be approximated as the saturation thickness of the aquifer in non-confined aquifers; and These are the permeability coefficients in the x and y directions, characterizing the anisotropic ability of a medium to allow groundwater to pass through. Values ​​within a geologically reasonable range are considered, for example, when assuming a sandy stratum. exist ; The specific yield of an aquifer characterizes the ability of a non-confined aquifer to release or store water under gravity, such as sandy soil. Typical range is These are all inherent parameters of the parametric geological structure network model. Representing the source-sink intensity, i.e., the vertical recharge or discharge per unit area per unit time applied to the aquifer, it is derived from the candidate physical disturbance scenario. For example, for pumping activities at pumping point (X,Y) with a flow rate Q(t), this parameter is set by the following formula:

[0056] in, Represents the pumping flow rate over time, obtained from manual pumping activity logs or pump rated parameters in an external multi-source database. For example, the constant flow rate of a production well can be 0.02 m³ / s. (X,Y) represents the location coordinates of a point source or sink (such as a pumping well), obtained from well location coordinate data in a geographic information system. Representing Dirac A function is a mathematical tool used in physics and engineering to mathematically express a physical quantity that acts at a single point (such as the flow rate of a well in this scheme) in continuous spatial field equations. The negative sign indicates that pumping is the removal of water from the aquifer. This formula serves as a physical model based on the law of conservation of mass to quantitatively describe and predict the dynamic process of the evolution of groundwater head h in unconfined aquifers over time and space, driven by various source and sink terms W(x,y,t) (such as pumping and rainfall).

[0057] The second step involves solving the physical model to obtain the theoretical response time series for all observation points. After setting the initial conditions, boundary conditions, and the source and sink terms determined in the first step, the system uses numerical methods, such as the finite difference method or the finite element method, to solve the aforementioned partial differential equations. The solution time range covers the entire period of anomaly occurrence. The result is the head height value of each node in the parameterized geological structure network model at each time step. The system specifically extracts the model nodes corresponding to relevant observation points in the actual observation network and exports their head height value sequences throughout the simulation period, forming the theoretical time-series response data for each observation point. ,in Represents the observation point number. Represents time.

[0058] The third step is to combine multiple sets of one-dimensional theoretical time series into a unified theoretical overall response model that is equivalent to the measured data structure: the system processes each set of theoretical time series response data that has just been generated. Perform the exact same procedures as processing measured data to ensure fairness and consistency in subsequent comparisons. Specifically, this step includes: 1) Extracting theoretical time-frequency features: The system extracts theoretical time-series response data for each... A continuous wavelet transform (CWT) is performed to generate a two-dimensional time-frequency energy spectrum. Based on this spectrum, the system calculates a series of preset characteristic indices to form the theoretical time-frequency feature vector for each observation point. The vector includes at least: Main energy frequency: the frequency corresponding to the point of maximum energy density in the time frequency energy spectrum.

[0059] Energy concentration duration: the length of time during which energy exceeds a specific threshold.

[0060] Integrating energy in a specific frequency band within a physically defined frequency band, such as the solid tide response band or the rainfall response band.

[0061] 2) Calculate the theoretical spatiotemporal correlation: The system calculates the theoretical spatiotemporal correlation between any two observation points k and j. This relationship is quantified in the following ways: Theoretical time delay: calculated by measuring the theoretical time-series response data corresponding to observation points k and j. and Find the time offset corresponding to the peak value of the cross-correlation function, which is the theoretical signal propagation time delay.

[0062] Theoretical amplitude attenuation ratio: Calculate the time-frequency eigenvectors corresponding to observation points k and j. and The ratio of norms, for example, the Euclidean norm, i.e. This ratio reflects the theoretical attenuation of signal energy between two points.

[0063] 3) Constructing the theoretical overall response model: The system constructs the theoretical time-frequency feature vectors of all observation points. Theoretical spatiotemporal correlation with all observation point pairs The data is packaged together to form a structured data volume, which is the theoretical overall response mode corresponding to the candidate physical disturbance scenario i. .this The data structure is completely equivalent to the actual overall response pattern, and then it is sent to the matching and judgment module for comparison with the actual overall response pattern.

[0064] S4. Based on the measured data of all relevant observation points during the period of abnormal observation data, construct the overall measured response model.

[0065] In a preferred embodiment, constructing the measured overall response model includes: performing joint time-frequency transformation on the measured time-series data of all relevant observation points, and extracting the time-frequency feature vector of each observation point; Calculate the spatiotemporal correlation between time-frequency feature vectors at different observation points; By integrating all time-frequency feature vectors and their spatiotemporal correlations, a measured overall response mode is generated.

[0066] Specifically, in order to transform discrete, independent time-series observation data into a unified data structure that can comprehensively characterize the dynamic response of the system, i.e., the measured overall response pattern, this embodiment discloses a detailed construction method. This method upgrades the original one-dimensional time series to a structured data volume containing multiple dimensions of information such as frequency, time, and space, providing a high-dimensional, physically meaningful basis for subsequent matching with the simulation results of the physical model.

[0067] Step 1: Extract and integrate deep features from measured data of multiple observation points: The system retrieves all relevant observation points determined by the anomaly triggering step from the database. (k=1,2,...,N) Measured time series data during the abnormal period, such as water level data sequences. The system provides measured time-series data for each observation point. To perform joint time-frequency transformation, continuous wavelet transform (CWT) is typically used in engineering, and its formula is as follows:

[0068] in, These are the wavelet transform coefficients at the observation points, and the square of their moduli. This represents the energy density of the signal over time b and scale a (inversely proportional to frequency); It is an observation point The original measured signal; It is the mother wavelet function The complex conjugate wavelet basis function obtained after translation and scaling matches the local features of the signal at specific times and frequencies by shifting the mother wavelet by time (b) and scaling it by scale (a). Through this transformation, the original one-dimensional time-series data is converted into a two-dimensional time-frequency energy matrix. Based on this matrix, the system further calculates a series of preset feature indices to form the time-frequency feature vector for each observation point. This vector typically contains the principal energy frequency, the duration of energy concentration, and the integral energy within a specific physically meaningful frequency band. This formula decomposes non-stationary signals by performing inner product operations with a series of wavelet basis functions with different time-frequency localization characteristics, thereby revealing the dynamic laws governing the variation of signal energy with time and frequency.

[0069] The second step is to quantify the spatiotemporal coupling relationship between the responses of different observation points: the system does not view the characteristics of each point in isolation, but rather calculates the interrelationships between them. For any two observation points... and The system is based on their time-frequency feature vectors and Calculate spatiotemporal correlation This includes two key components: time delay. and response amplitude decay ratio .

[0070]

[0071]

[0072] in, The time delay required for a signal to propagate from observation point j to observation point k is represented by the time-frequency energy matrix of the two signals. and The value is obtained by two-dimensional cross-correlation calculation. For example, a value of 3600 seconds means that the signal response at point k is one hour later than that at point j. The mathematical operator represents the parameter value that maximizes the expression; in this case, it means finding the time shift that maximizes the correlation between the two time-frequency matrices. . The wavelet coefficient matrix represents the observation point j. The operations performed include: 1) taking the complex conjugate. ;2) Translate on the time axis . Represents a two-dimensional integral over the entire time-scale plane, used to compute the integral of two matrices over a given time-scale plane. Overall correlation. It is the ratio of the norms of the two time-frequency eigenvectors, reflecting the attenuation of signal energy during propagation. and They are vectors sum vector The Euclidean norm represents the overall intensity of the signal energy at observation points k and j.

[0073] The third step is to combine the extracted single-point features and inter-point relationships into a unified structured data volume: the system combines the time-frequency feature vectors of all relevant observation points. The spatiotemporal correlation between all observation point pairs , This data is organized into a comprehensive data structure. This structure is the measured overall response mode. Logically, it can be viewed as a weighted graph, where nodes are observation points and their attributes are time-frequency feature vectors. The edge's attribute is its spatiotemporal relationship. .this It fully describes the systematic response triggered by anomalous events throughout the entire observation network, including the distribution and evolution of energy in the three dimensions of time, frequency, and space.

[0074] S5. Match the measured overall response mode with all theoretical overall response modes, and through a feedback mechanism, with the goal of minimizing the difference between the measured overall response mode and the theoretical overall response mode, coordinately optimize the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model to obtain the matching judgment result.

[0075] In a preferred embodiment, the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model are coordinated and optimized, including: calculating the initial difference between the measured overall response mode and each theoretical overall response mode, and selecting the candidate physical disturbance scenario to be optimized based on the initial difference. The parameters of the candidate physical disturbance scenario to be optimized and the connection edge parameters of the parameterized geological structure network model associated with it are used as adjustable variables. Under the physical causal constraints determined by the initial parameter range and network topology, the adjustable variables are iteratively adjusted. The objective function of optimization is to minimize the difference between the newly generated theoretical overall response mode and the measured overall response mode.

[0076] Specifically, this embodiment discloses a detailed engineering method for identifying interference sources through collaborative optimization using a feedback mechanism. The core of this method is to transform static pattern matching into a dynamic parameter optimization process constrained by physical laws, thereby finding the most logically consistent explanation between uncertain geological models and incomplete observation data.

[0077] Step 1: Preliminary screening and ranking of multiple candidate physical interference scenarios: The system receives the measured overall response pattern generated in the previous steps. and multiple theoretical overall response modes The system calculates... With each Initial difference between This is used to quantify their similarity. The difference is calculated based on the Structural Similarity Index (SSIM) because it effectively assesses the combined similarity of two complex patterns in terms of structure, contrast, and brightness.

[0078] in, This represents the degree of difference between the i-th theoretical model and the measured model. Its value is in the range of [0,1]. The smaller the value, the more similar the two models are. It is the actual test mode. and the i-th theoretical model The structural similarity index between the two patterns has a value in the range of [0,1]. The larger the value, the more similar the two patterns are. It is a dimensionless parameter.

[0079] To apply SSIM to non-image pattern data P, the system performs the calculation as follows: Data matrixing: and The data is transformed into two numerical matrices of the same dimension. Specifically, the feature vectors of all nodes are... Arrange them in a predetermined order to form the first part of the matrix; [and] the relationships between all edges. (For example, delay and attenuation ratio) are also arranged in a predetermined order to form the second part of the matrix.

[0080] Block-based computation: SSIM is calculated using a multi-scale, window-based approach. Specifically, a fixed-size window is slid across the numerical matrix, and the similarity components of brightness, contrast, and structure are calculated within each window. These components are then weighted to obtain the SSIM value for that window.

[0081] Mean aggregation: final The value is the average of the SSIM values ​​for all windows.

[0082] After the calculation is completed, the system will process all... If the response is below a preset threshold, for example, a threshold of 0.5, the corresponding theoretical overall response mode, its associated candidate physical disturbance scenarios, and the parameterized geological structure network model used are packaged and marked as an optimization object for use in subsequent collaborative tuning steps.

[0083] The second step involves iterative bidirectional approximation for each selected optimization object. This process aims to minimize the difference between theory and experiment by fine-tuning the model parameters, simulating the process of experts repeatedly deducing and correcting based on physical laws. This "bidirectional approximation" is reflected in two aspects: first, searching for the optimal solution in the parameter space; and second, continuously verifying and constraining the search direction using physical laws during this process. For each optimization object selected in the first step, the following iterative optimization process is executed: 1) Define the optimization problem: Start an optimization solver, for example, the finite-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS) suitable for large-scale nonlinear optimization problems.

[0084] Adjustable variables: The candidate physical disturbance scene parameters (such as disturbance intensity and duration) within the optimization object, along with the connection edge parameters (such as permeability coefficient and specific yield) of the associated parameterized geological structure network model, are collectively defined as the variables that the optimization algorithm needs to adjust. .

[0085] Objective function: The objective function for optimization is defined as minimizing the newly generated theoretical global response mode after adjustment. Overall response mode as measured Difference between The method for calculating this difference degree D is the same as that in the first step. Exactly the same.

[0086] 2) Applying physical constraints: This iterative adjustment process is not an unconstrained mathematical fit, but is subject to strict physical constraints to ensure the realistic rationality of the adjustment results.

[0087] Parameter range constraint: All adjustments to the parameters of the connecting edges must be made within their preset initial physical parameter range. Within. For example, the initial parameter range of the equivalent permeability coefficient K for a certain edge is set to... Therefore, the value of K must not exceed this range during the optimization process. This constraint is achieved by directly inputting the variable limits into the L-BFGS solver.

[0088] Physical causality constraint: In each iteration, when the solver provides a new set of parameters And generate new theoretical models Then, the system will immediately perform a causality check. Specifically, the system will check... Time delay between observation points Does it still strictly conform to the topology of the parameterized geological structure network model it is based on, i.e., the response time of downstream observation points cannot be earlier than that of upstream points? If causal reversal occurs, a large penalty term is imposed on the objective function value of this iteration, thereby guiding the optimization algorithm to abandon this search direction.

[0089] 3) Determine convergence and output: When the objective function... The decrease is less than the preset convergence tolerance (e.g.) The iteration stops when the number of iterations reaches a preset maximum value (e.g., 200). The system outputs the final difference score obtained after the optimization converges. The complete set of optimized parameters corresponding to achieving this degree of difference .

[0090] This entire set of outputs will serve as input for subsequent matching and determination steps. Through this optimization process "guided" by physical laws, the present invention ensures that the final solution found is not only mathematically well-fitting to the observed data, but also physically reasonable and interpretable.

[0091] In a further preferred embodiment, the matching determination result includes: during the iterative adjustment process, recording the parameter modification amount and the degree of difference improvement corresponding to each round of adjustment to form an optimization process record; When the iterative adjustment meets the preset convergence condition, the comprehensive score of each candidate physical interference scenario is calculated based on the final difference, the quantitative evaluation of the rationality of parameter adjustment reflected in the optimization process record, and the quantitative evaluation of the complexity of the candidate physical interference scenario. The candidate physical interference scene with the highest comprehensive score is selected as the final interference scene in the matching judgment result.

[0092] Specifically, this embodiment details how, after the collaborative optimization process, a comprehensive evaluation is conducted to make a final judgment decision to determine the most likely interference scenario. This method aims to integrate evaluation indicators from multiple dimensions to avoid misjudgments caused by a single indicator, ensuring the scientific rigor and robustness of the final decision.

[0093] The first step is to quantify and record the behavior during the collaborative tuning process: During iterative adjustments, the system not only aims to reduce the degree of difference but also simultaneously records the process parameters related to each iteration. Specifically, after the nth iteration, the system records the parameter modification vector corresponding to that iteration. And the improvement in difference .

[0094] Parameter Modification Vector This vector contains the normalized changes of all adjusted model parameters (such as penetration coefficient and disturbance intensity) relative to their initial values. For the p-th parameter, its normalized change after the nth iteration is calculated as follows: ,in, It is the value of the p-th parameter after the nth iteration. This is the initial value of the p-th parameter. and These are the upper and lower limits of the preset physical reasonable range for the p-th parameter. In this way, the amount of modification to all parameters is unified into a dimensionless interval [0,1], which facilitates subsequent comprehensive evaluation.

[0095] Improvement in difference This value represents the effectiveness of the nth iteration in reducing the difference compared to the (n-1)th iteration. The calculation is as follows: ,in, and These are the differences after the nth and (n-1)th iterations, respectively.

[0096] The system will process data It is stored together with the new theoretical overall response pattern generated in each round to form a complete optimized trajectory record.

[0097] The second step, after the iteration converges, is to construct a multi-objective comprehensive scoring function to perform a final evaluation on each optimized candidate physical disturbance scenario: when the previous optimization step meets the preset convergence condition, such as the degree of improvement in difference after three consecutive iterations. Less than The optimization stops when the maximum number of iterations (e.g., 200) is reached. At this point, the system calculates the final comprehensive score for each optimization object i. :

[0098] in, It is the first The overall score of each optimization object. The weighting coefficients for the final difference degree, the rationality of parameter adjustment, and the complexity of candidate physical interference scenarios are set based on prior knowledge and expert experience, and their sum is 1. For example, . It is the first After optimization convergence, the degree of difference between the theoretical model and the measured overall response model of an optimization object is determined. The item directly quantifies the quality of the match; the higher the value, the better the match. The normalized norm of the cumulative parameter modification amount for the optimization object i during the entire tuning process is calculated as follows: First, obtain the cumulative parameter modification vector. ,in These are the optimized final parameters. These are the initial parameters before optimization; then, the vector is normalized to be unaffected by the number and dimensions of the parameters; finally, the Euclidean norm of the normalized vector is calculated. This term reflects the drastic degree of modification to the initial model assumptions required to achieve a match, therefore, The value represents the rationality of parameter adjustment; the higher the value, the smoother and more reasonable the adjustment. It is a quantitative evaluation of the complexity of the i-th candidate physical disturbance scenario itself. It is preset according to Occam's razor principle. For example, simpler or more common disturbance types such as single rainfall are set to score higher, such as 0.9; while complex or rare combined disturbance types such as rainfall superimposed on far-field structural stress changes are set to score lower, such as 0.4.

[0099] Step 3: Make a final decision based on the comprehensive score and output the results: The system's comprehensive score for all optimization objects. The candidate physical interference scenarios are sorted, and those with the highest scores are selected as the final identified interference scenarios behind the anomalies in the observed data. This determination, along with the final optimized model parameters and the corresponding final theoretical overall response mode, constitutes the matching determination result. This structured matching determination result will be passed to the subsequent data decoupling module as the basis for eliminating interference.

[0100] S6. Based on the interference scenario finally determined in the matching judgment result, decouple the interference component corresponding to the scenario from the original observation data and output the purified fluid observation data.

[0101] In a preferred embodiment, the interference components corresponding to the scene are decoupled from the original observation data, and the purified fluid observation data is output, including: parsing the theoretical interference signal waveform corresponding to the finally determined interference scene and for each observation point from the matching judgment result; Using the theoretical interference signal waveform as a reference signal, an adaptive filtering algorithm is used to remove components related to the morphology of the reference signal from the original observation data of the corresponding observation points, thus obtaining purified fluid observation data.

[0102] Specifically, this embodiment discloses a detailed engineering method for separating and removing corresponding interference components from the original observation data after obtaining the finally determined interference scenario. This method aims to decouple signals in an adaptive manner to ensure that, while removing known interference, it can retain as many weak unknown signals, such as earthquake precursors, as possible in the data.

[0103] Step 1: Reconstructing the Theoretical Interference Signal Based on the Matching Judgment Results: The system first determines the set of affected observation points that need signal removal. This set is a list of observation points whose original signals contain interference energy characteristics, identified through time-frequency analysis during the construction of the measured overall response model. Then, the system calls the matching judgment results, which include the finally determined interference scenario and the optimal model parameters after collaborative tuning. These optimal model parameters are configured in the parameterized geological structure network model, and a forward physical simulation is re-executed using the finally determined interference scenario as the input driving force (i.e., source-sink terms or boundary conditions). By solving the model's governing equations, the theoretical time-series response data caused by the interference scenario at each affected observation point in the aforementioned set is directly calculated. This data is then identified as the theoretical interference signal waveform specifically for that observation point, used for signal removal.

[0104] The second step involves adaptive decoupling of the original observation data to eliminate interference while preserving as many potential, unknown, and effective signals as possible: for each affected observation point... Instead of simply subtracting the theoretical interference signal waveform from the original observation data, the system employs an adaptive filtering algorithm. This algorithm uses the measured original data from the observation point as the input signal and the theoretical interference signal waveform as the reference signal.

[0105]

[0106] in, It is the final output of the system, specific to the observation point. The purified fluid observation data. It is an observation point Time series of raw measured fluid observation data during the anomaly. It is for the observation point The theoretical interference signal waveform time series. Key parameters. It is a dimensionless adaptive scaling factor. It is not a fixed constant, but is dynamically calculated and updated by an adaptive algorithm (such as the Least Mean Square algorithm, LMS) at each time step of the processed signal. Its setting is based on minimizing the output signal. The energy (i.e., mean square value) is used to achieve optimal interference cancellation, and its reasonable range can be set, for example, between 0.5 and 2.0, to cope with the uncertainty of the model. In the physical sense of signal processing, this is equivalent to... Remove reference signal The most morphologically relevant component. This adaptive mechanism can compensate for the unavoidable small amplitude differences between the theoretical model and the actual physical process, avoiding the problems of "over-subtraction" or "under-subtraction" caused by the incomplete accuracy of the model, thus protecting the weak real signals that may overlap with the morphological parts of the interfering signals. Finally, the purified fluid observation data from all relevant observation points. These data are combined to form a high-quality dataset that can be used for subsequent seismic analysis.

[0107] In a further preferred embodiment, after outputting the purified fluid observation data, a model self-optimization step is also included: the finally determined interference scenario, the collaboratively optimized parameterized geological structure network model connection edge parameters, and the measured overall response mode are encapsulated into a training sample and stored in the historical case library. Based on training samples accumulated in the historical case library, statistical corrections are made to the initial parameter range or the mapping relationship of cause-effect knowledge templates in the parameterized geological structure network model.

[0108] Specifically, this embodiment discloses a model self-optimization method that enables the entire interference identification system to have self-learning and evolution capabilities. Its engineering goal is to continuously and automatically calibrate and improve the prior knowledge model inside the system by utilizing the experience of each successful interference identification.

[0109] The first step is to solidify each successful interference identification process into a structured record that can be learned: After the system successfully outputs the purified fluid observation data, it does not discard the process information of this processing. Instead, it packages and encapsulates the detailed description of the finally determined interference scenario, the connection edge parameters of the parameterized geological structure network model finally determined after collaborative optimization, and the key information such as the matching degree and decoupling residuals that characterize the actual observation effect into a training sample with a timestamp and stores it in the system's historical case library.

[0110] The second step involves statistically correcting the initial parameter range of the parametric geological structure network model. This optimization process is automatically triggered periodically or after a certain number of training samples has been accumulated, such as exceeding 50 successful cases. This correction is applied to a specific, uncertain connection edge parameter in the model, such as the first... Equivalent permeability coefficient of the strip The system will retrieve all training samples involving this parameter from the historical case library, and extract the final value determined after collaborative tuning to form a parameter set. U represents the total number of "training samples" or "successful cases" retrieved from the historical case database that are related to the m-th connection edge parameter. The system performs statistical analysis on this parameter set and calculates its mean. and standard deviation Then, update the initial parameter range of the parameter accordingly.

[0111]

[0112] in, It is the first The new initial parameter range for each connected edge parameter; and These are the mean and standard deviation of its historical successful value set, respectively; This is a preset confidence coefficient, typically set to 2 or 3, representing a confidence interval of approximately 95% or 99%. Through this method, the model's prior knowledge gradually converges from a broad, fuzzy initial setting to a more precise range that better reflects the actual hydrogeological characteristics of the region, thereby improving the efficiency and accuracy of future interference identification.

[0113] The third step is to optimize the mapping relationship of the cause-effect knowledge template, enabling it to identify and learn new or atypical interference patterns. The system pays special attention to training samples in the historical case library that had low matching rates during the initial scene generation stage but were ultimately confirmed as correct interference sources through collaborative optimization. For such cases, the system considers the mapping relationship between the initial abnormal morphology of the observed data and the finally determined physical cause to be missing or inaccurate in the template library. Therefore, the system will supplement or update the cause-effect knowledge template with a new high-confidence mapping rule based on the combination of the case's initial morphological features and the finally determined interference scene, thereby enhancing the system's cognitive ability to recognize complex or unknown interference patterns and achieving self-expansion of knowledge.

[0114] Example 2

[0115] Please see Figure 2 As shown, based on Embodiment 1, the second aspect of the present invention provides an earthquake underground fluid observation data interference identification system based on time-frequency response characteristics, which includes: a model building module, a scene generation module, a simulation and deduction module, a pattern building module, a matching and determination module, and a data decoupling module.

[0116] The model building module is used to acquire multi-source geological information and observation network topology of the target area, and to construct a parameterized geological structure network model with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges, wherein the parameters of the connecting edges have an initial parameter range.

[0117] The scene generation module is used to monitor earthquake underground fluid observation data. When anomalies in the observation data are identified, it combines data from an external multi-source information database to generate at least one candidate physical interference scene containing the interference type, intensity, and spatiotemporal location.

[0118] The simulation and deduction module is used to take the parameters of the candidate physical disturbance scenarios as input, drive the parameterized geological structure network model to perform parallel simulation, and generate a theoretical overall response mode for each candidate physical disturbance scenario.

[0119] The pattern building module is used to construct the overall measured response pattern based on the measured data of all relevant observation points during the period of observational data anomalies.

[0120] The matching and determination module is used to match the measured overall response mode with all theoretical overall response modes. Through a feedback mechanism, with the goal of minimizing the difference between the measured overall response mode and the theoretical overall response mode, the module performs collaborative optimization of the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model to obtain the matching and determination results.

[0121] The data decoupling module is used to decouple the interference components corresponding to the interference scenario determined in the matching judgment result from the original observation data and output the purified fluid observation data.

[0122] The above content is merely an example and illustration of the concept of the present invention. Those skilled in the art can make various modifications or additions to the specific embodiments described, or use similar methods to replace them, as long as they do not deviate from the concept of the invention or exceed the scope defined by the present invention, and all such modifications and additions should fall within the protection scope of the present invention.

Claims

1. A method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics, characterized in that, include: S1. Obtain multi-source geological information and observation network topology of the target area, and construct a parameterized geological structure network model with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges, wherein the parameters of the connecting edges have an initial parameter range; S2. Monitor earthquake underground fluid observation data. When anomalies in the observation data are identified, combine the data from the external multi-source information database to generate at least one candidate physical interference scenario that includes the type, intensity, and spatiotemporal location of the interference. S3. Using the parameters of the candidate physical disturbance scenarios as input, drive the parameterized geological structure network model to perform parallel simulation and generate a theoretical overall response mode for each candidate physical disturbance scenario. Generate a theoretical overall response pattern for each candidate physical disturbance scenario, including: The parameters of the candidate physical disturbance scenarios are transformed into source and sink terms or boundary conditions of the parameterized geological structure network model. The theoretical time-series response data of all relevant observation points were calculated by numerically solving the governing equations of the parameterized geological structure network model. Joint time-frequency analysis and spatial correlation analysis are performed on theoretical time-series response data to synthesize the overall theoretical response mode; S4. Based on the measured data of all relevant observation points during the period of abnormal observation data, construct the overall measured response model; Constructing the measured overall response model, including: Perform joint time-frequency transformation on the measured time-series data of all relevant observation points, and extract the time-frequency feature vector of each observation point; Calculate the spatiotemporal correlation between time-frequency feature vectors at different observation points; By integrating all time-frequency feature vectors and their spatiotemporal correlations, a measured overall response mode is generated; S5. Match the measured overall response mode with all theoretical overall response modes, and through a feedback mechanism, with the goal of minimizing the difference between the measured overall response mode and the theoretical overall response mode, coordinate the optimization of the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model to obtain the matching judgment result. S6. Based on the interference scenario finally determined in the matching judgment result, decouple the interference component corresponding to the scenario from the original observation data and output the purified fluid observation data.

2. The method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics according to claim 1, characterized in that, A parameterized geological structure network model is constructed, with observation wells and potential disturbance sources as nodes and hydraulic connections as connecting edges, including: Based on regional hydrogeological data, observation wells and potential sources of disturbance are mapped as network nodes with spatial attributes, and directed connecting edges representing groundwater flow channels are established according to aquifer distribution and fault structures to construct the network topology of the geological structure network model. Each connecting edge is assigned physical parameters that characterize the hydrodynamic properties of the medium. These physical parameters include at least the permeability coefficient and the water storage coefficient. Based on regional stratigraphic lithology data or historical pumping test data, set initial parameter ranges for physical parameters; Based on a parameterized geological structure network model and combined with regional hydrological boundary conditions, a set of groundwater unsteady flow control equations is established to describe the dynamic response of the network model, thus forming a forward physical simulation model.

3. The method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics according to claim 1, characterized in that, Generate at least one candidate physical interference scenario that includes interference type, intensity, and spatiotemporal location, including: Extract preliminary morphological features of anomalies from observed data to form anomaly feature descriptors; Call a cause-effect knowledge template library that defines the mapping relationship between physical interference causes and data anomaly patterns, match the anomaly feature descriptors with the library, and obtain a list of potential physical interference causes; Using a list of potential physical interference causes as an index, external multi-source information databases are retrieved, and the matched external information events are instantiated into parameterized scene descriptions to form candidate physical interference scenarios. The data in the external multi-source information database includes at least two of the following: meteorological precipitation data, reservoir scheduling records, logs of manual water pumping activities, and information on surrounding earthquake sources.

4. The method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics according to claim 1, characterized in that, The parameters of the candidate physical disturbance scenarios and the connection edge parameters of the parameterized geological structure network model are jointly optimized, including: Calculate the initial difference between the measured overall response mode and each theoretical overall response mode, and select candidate physical disturbance scenarios to be optimized based on the initial difference. The parameters of the candidate physical disturbance scenario to be optimized and the connection edge parameters of the parameterized geological structure network model associated with it are used as adjustable variables. Under the physical causal constraints determined by the initial parameter range and network topology, the adjustable variables are iteratively adjusted. The objective function of optimization is to minimize the difference between the newly generated theoretical overall response mode and the measured overall response mode.

5. The method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics according to claim 4, characterized in that, The matching result includes: During the iterative adjustment process, the amount of parameter modification and the degree of improvement in difference corresponding to each round of adjustment are recorded to form an optimization process record; When the iterative adjustment meets the preset convergence condition, the comprehensive score of each candidate physical interference scenario is calculated based on the final difference, the quantitative evaluation of the rationality of parameter adjustment reflected in the optimization process record, and the quantitative evaluation of the complexity of the candidate physical interference scenario. The candidate physical interference scene with the highest comprehensive score is selected as the final interference scene in the matching judgment result.

6. The method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics according to claim 1, characterized in that, The interference components corresponding to this scene are decoupled from the original observation data, and the purified fluid observation data is output, including: The theoretical interference signal waveform corresponding to the finally determined interference scenario and for each observation point is extracted from the matching judgment results; Using the theoretical interference signal waveform as a reference signal, an adaptive filtering algorithm is used to remove components related to the morphology of the reference signal from the original observation data of the corresponding observation points, thus obtaining purified fluid observation data.

7. The method for identifying interference in seismic subsurface fluid observation data based on time-frequency response characteristics according to claim 1, characterized in that, After outputting the purified fluid observation data, a model self-optimization step is also included: The final determined interference scenario, the edge parameters of the parameterized geological structure network model after collaborative optimization, and the measured overall response mode are encapsulated into a training sample and stored in the historical case library. Based on training samples accumulated in the historical case library, statistical corrections are made to the initial parameter range or the mapping relationship of cause-effect knowledge templates in the parameterized geological structure network model.

8. A seismic subsurface fluid observation data interference identification system based on time-frequency response characteristics, characterized in that, include: The model building module is used to acquire multi-source geological information and observation network topology of the target area, and to construct a parameterized geological structure network model with observation wells and potential interference sources as nodes and hydraulic connections as connecting edges, wherein the parameters of the connecting edges have an initial parameter range. The scene generation module is used to monitor earthquake underground fluid observation data. When anomalies in the observation data are identified, it combines data from an external multi-source information database to generate at least one candidate physical interference scene containing the interference type, intensity, and spatiotemporal location. The simulation and deduction module is used to take the parameters of the candidate physical disturbance scenarios as input, drive the parameterized geological structure network model to perform parallel simulation, and generate a theoretical overall response mode for each candidate physical disturbance scenario. Generate a theoretical overall response pattern for each candidate physical disturbance scenario, including: The parameters of the candidate physical disturbance scenarios are transformed into source and sink terms or boundary conditions of the parameterized geological structure network model. The theoretical time-series response data of all relevant observation points were calculated by numerically solving the governing equations of the parameterized geological structure network model. Joint time-frequency analysis and spatial correlation analysis are performed on theoretical time-series response data to synthesize the overall theoretical response mode; The pattern construction module is used to construct the overall measured response pattern based on the measured data of all relevant observation points during the period of observation data anomaly. Constructing the measured overall response model, including: Perform joint time-frequency transformation on the measured time-series data of all relevant observation points, and extract the time-frequency feature vector of each observation point; Calculate the spatiotemporal correlation between time-frequency feature vectors at different observation points; By integrating all time-frequency feature vectors and their spatiotemporal correlations, a measured overall response mode is generated; The matching and determination module is used to match the measured overall response mode with all theoretical overall response modes. Through a feedback mechanism, with the goal of minimizing the difference between the measured overall response mode and the theoretical overall response mode, the module performs collaborative optimization of the parameters of the candidate physical disturbance scenario and the connection edge parameters of the parameterized geological structure network model to obtain the matching and determination results. The data decoupling module is used to decouple the interference components corresponding to the interference scenario determined in the matching judgment result from the original observation data and output the purified fluid observation data.

Images