Multi-field Coupled Inversion Analysis Method for Water-Rich Characteristics of Karst Disaster Sources

Abstract

This invention discloses a multi-field coupled inversion analysis method for water-rich characteristics of karst disaster sources. The method acquires tracer penetration curves and hydraulic head variation curves of the target area and extracts travel time features; constructs a unified parameter state vector including permeability coefficient, storage coefficient, and effective porosity, and establishes a forward model describing solute transport and groundwater flow processes; calculates the sensitivity matrix of the unified parameter state vector to observational data based on the forward model, and constructs a joint objective function including a spatial adaptive regularization term; iteratively solves the objective function using a prediction-correction alternating update strategy to obtain an optimized parameter distribution and delineate the location of disaster-causing structures. This invention solves the problems of unidirectional error accumulation and fragmented physical parameter coupling in traditional serial inversion by constructing a multi-physics synchronous coupling mechanism and adaptive constraint strategy, thus improving the identification accuracy and resolution of water-rich disaster-causing structures under complex geological conditions.

Description

Technical Field

[0001] This invention belongs to the field of underground engineering safety monitoring and disaster prevention and control, and in particular, it is a multi-field coupled inversion analysis method for water-rich characteristics of karst disaster sources. Background Technology

[0002] In underground engineering construction and deep resource development, accurately identifying water-rich and hazardous structures such as fault fracture zones and karst conduits can prevent water and mud inrush disasters. This relies on characterizing the hydrogeological parameters of the target area, especially the spatial distribution of parameters such as permeability coefficient, water storage coefficient, and effective porosity, which directly determine the seepage path of groundwater and the transport patterns of solutes. Therefore, obtaining high-resolution parameter field distributions helps to reveal the location and connectivity of hidden hazardous structures.

[0003] Currently, the mainstream techniques for integrated exploration combining tracer experiments and hydraulic tomography employ a sequential two-stage inversion strategy. The first stage uses travel time data from tracers to invert the slowness or velocity distribution of the formation through X-ray tomography (such as the SIRT algorithm), mapping it to a preliminary permeability field. The second stage uses this preliminary distribution as a priori model or initial condition, incorporating hydraulic head monitoring data for hydraulic tomography inversion to estimate hydrogeological parameters. This method, to some extent, utilizes the complementarity of the two types of data and is widely used in hydrogeological exploration.

[0004] However, the aforementioned serial inversion strategy suffers from problems in both theory and application, namely, the unidirectional propagation of errors and the disconnect between the coupling and physical processes. Therefore, further research and innovation are needed to address these issues in existing technologies. Summary of the Invention

[0005] The purpose of this invention is to provide a multi-field coupled inversion analysis method for the water-rich characteristics of karst disaster sources, in view of the problems existing in the prior art.

[0006] The technical solution, on the one hand, provides a multi-field coupled inversion analysis method for the water-rich characteristics of karst disaster sources, including:

[0007] Obtain the tracer penetration curve and water head change curve of the target area, extract features from the tracer penetration curve, and generate a travel time feature set;

[0008] A unified parameter state vector containing permeability coefficient, water storage coefficient and effective porosity is constructed to establish a forward model describing solute transport process and groundwater flow process;

[0009] Based on the travel time feature set, the head change curve, and the output of the forward model, the tracer observation residual and the head observation residual are calculated.

[0010] Based on the forward model, the sensitivity matrix of the unified parameter state vector to the travel time feature set and the head change curve is calculated, and a joint objective function containing a spatial adaptive regularization term is constructed based on this (sensitivity matrix).

[0011] The joint objective function is solved iteratively to obtain the optimized unified parameter state vector, and the location of the water-rich disaster-causing structure is delineated based on the optimized parameter distribution.

[0012] The joint objective function also includes the weighted sum of squares of the tracer observation residuals and the head observation residuals.

[0013] Beneficial effects: By constructing a multi-physics synchronous coupling mechanism and an adaptive constraint strategy, this invention solves the problems of unidirectional error accumulation and fragmented physical parameter coupling in traditional serial inversion, thereby improving the identification accuracy and resolution of water-rich disaster-causing structures under complex geological conditions. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0014] Figure 1 This is a flowchart of a multi-field coupled inversion analysis method for water-rich characteristics of karst disaster sources, provided in an embodiment of this application.

[0015] Figure 2 A flowchart illustrating the prediction-correction alternating update strategy provided in the embodiments of this application.

[0016] Figure 3 A flowchart illustrating the calculation of tracer travel time values ​​provided in an embodiment of this application.

[0017] Figure 4 This is a flowchart illustrating the location of water-rich disaster-causing structures based on optimized parameter distribution, provided as an embodiment of this application. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0019] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0020] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0021] Correspondingly, the model errors generated by the inversion, such as the path deviation in the high-permeability zone caused by the straight ray assumption, will solidify and propagate to later stages, and it will be difficult to use subsequent head data for retrospective correction, resulting in cumulative bias in the inversion results.

[0022] Furthermore, the slowness parameter in the solute transport process is physically dependent on the local permeability coefficient and hydraulic gradient. However, existing methods neglect dynamic physical coupling, which leads to distortion of parameter estimation in strongly heterogeneous media.

[0023] Based on this, existing methods mostly adopt globally uniform regularization constraints, which cannot adapt to the spatial sensitivity differences caused by uneven ray coverage density, and it is difficult to improve the resolution of key areas while ensuring stability.

[0024] To solve these problems, combined with Figures 1 to 4 The present invention will be specifically described through the following embodiments.

[0025] On the one hand, an exemplary scheme for multi-field coupled inversion analysis of water-rich characteristics of karst disaster sources is provided, describing the technical process from field data acquisition, preprocessing, modeling, joint inversion to delineation of disaster-causing structures. This solves the problem of overcoming the multiple solutions of single-physical-field inversion under complex geological conditions, enabling accurate location and identification of water-rich disaster-causing structures (such as karst caves, faults, and fracture zones).

[0026] Step 101: Obtain the tracer penetration curve and head change curve of the target area, extract features from the tracer penetration curve, and generate a travel time feature set.

[0027] In this embodiment, based on geological, geophysical, and drilling data of the engineering research area, the test site location is scientifically delineated, and a field tracing and pumping water test monitoring plan is formulated. For the tracing test, a multi-source monitoring system can be optionally set up, that is, tracers are released at multiple locations and monitored in multiple boreholes. Regarding the tracer release strategy, single-point release and single-point monitoring, single-point release and multi-point monitoring, or multi-source release and multi-point monitoring can be adopted to obtain solute transport information along multiple pathways.

[0028] To avoid interference between different tests, after the previous test is completed, the tracer concentration in the monitoring well must be monitored and confirmed to have recovered to near the initial background concentration before subsequent tests can be conducted. As an optional operating standard, a recovery is generally considered to be achieved when the monitored concentration change recovers to less than 5% of the initial peak concentration, and the next round of tests can be carried out. For water injection tests, the well depth of the filter pipe sections in each borehole should be set at the same or approximately the same height as much as possible to approximately satisfy the two-dimensional or quasi-three-dimensional flow field assumptions. In underground engineering such as tunnels, tunnels, and roadways, the monitoring position can also be flexibly arranged in a section or three-dimensional space in front of the working face by combining with advanced drilling construction.

[0029] Meanwhile, to ensure the accuracy of the inversion permeability coefficient and porosity, the collected head observation dataset should include unsteady flow pumping test data as much as possible.

[0030] After acquiring the raw monitoring data, the tracer breakthrough curve needs to be preprocessed and its features extracted. Since the raw data often contains environmental noise, the time-concentration curve is smoothed and denoised. Further, key temporal features are extracted from the denoised curve to generate a travel time feature set. The travel time feature set not only includes the peak arrival time of the tracer concentration but can also include time points corresponding to different cumulative mass recoveries, such as t0. _25% t _50% t _75% These parameters are used to reflect the speed of solute transport and dispersion characteristics in a medium. Similarly, for head change curves, water level fluctuation noise is removed by means of moving average, wavelet transform, or Kalman filtering, and then water level drawdown data at different times are extracted.

[0031] Step 102: Construct a unified parameter state vector containing permeability coefficient, water storage coefficient and effective porosity, and establish a forward model describing the solute transport process and groundwater flow process.

[0032] In other words, a unified parameter state vector containing permeability coefficient, water storage coefficient and effective porosity is constructed to establish a flow forward model describing the groundwater flow process and a solute transport forward model describing the solute transport process.

[0033] Accordingly, a unified parameter state vector was constructed to achieve simultaneous inversion of multiple physics fields. Specifically, three types of hydrogeological parameters—permeability coefficient, water storage coefficient, and effective porosity—were integrated into a set of unknowns to be solved. Establishing a forward model is the foundation of the inversion, and this model includes the following two coupled physical processes:

[0034] First, the flow model describing the groundwater flow process follows the unsteady flow control equations of saturated porous media.

[0035] Secondly, a solute transport model describing the migration and transformation of tracers in groundwater follows the convection-dispersion equation.

[0036] In a forward model, given the estimated values ​​of the parameters, the theoretical tracer penetration curve and head change curve can be simulated through numerical calculations, such as the finite difference method or the finite element method, i.e., the calculated values.

[0037] Step 103: Based on the travel time feature set, the head change curve, and the output of the forward model, calculate the tracer observation residual and the head observation residual.

[0038] Alternatively, based on the travel time feature set, the head change curve, and the travel time and head calculation values ​​output by the forward model, the tracer observation residuals and head observation residuals are calculated.

[0039] After obtaining the field-measured observations (i.e., travel time feature sets and head change curves) and the calculated values ​​output by the forward model, the differences between the two are calculated. Specifically, for tracer data, the difference between the observed and calculated values ​​for each monitoring channel at different characteristic time points is calculated to obtain the tracer observation residual; for head data, the difference between the observed and calculated drawdown for each monitoring well at different times is calculated to obtain the head observation residual. Furthermore, the deviation between the current parameter estimates and the actual geological environment is quantified, which is the direct driving force for subsequent optimization algorithms to update parameters.

[0040] Step 104: Calculate the sensitivity matrix of the unified parameter state vector to the travel time feature set and the head change curve based on the forward model, and construct a joint objective function containing a spatial adaptive regularization term based on this (sensitivity matrix).

[0041] The joint objective function also includes the weighted sum of squares of the tracer observation residuals and the head observation residuals.

[0042] The sensitivity matrix describes the extent to which a small change in each parameter to be inverted (such as the K value of a grid cell) will alter the observed data (such as the head or travel time at a monitoring point). This matrix is ​​typically calculated using the perturbation method or the adjoint state method.

[0043] In this step, a joint objective function is constructed. This function not only includes the weighted sum of squares of the tracer observation residuals and head observation residuals, but also introduces a spatially adaptive regularization term. Unlike traditional globally unified regularization, the spatially adaptive regularization term dynamically adjusts the constraint strength according to the sensitivity of parameters in different regions. Specifically, in regions with dense ray coverage and strong data constraints, the regularization weight is reduced to highlight the true information of the data; in regions with sparse ray coverage and weak data constraints, the regularization weight is increased to suppress spurious anomalies and numerical oscillations. This construction method balances the resolution and stability of the inversion.

[0044] Step 105: Iteratively solve the joint objective function to obtain the optimized unified parameter state vector, and delineate the location of the water-rich disaster-causing structure based on the optimized parameter distribution.

[0045] Building upon this, nonlinear optimization algorithms, such as the Gauss-Newton method or the Levenberg-Marquardt (LM) algorithm, are employed to minimize the joint objective function, which involves iterative updates. Starting from the initial parameter guesses, the unified parameter state vector is continuously corrected using the gradient descent direction until the objective function converges or reaches the preset number of iterations.

[0046] Once the optimized unified parameter state vector is obtained, the spatial distribution field of permeability coefficient, water storage coefficient, and effective porosity within the target area can be obtained. Based on the parameter distribution and combined with prior geological knowledge, the location of water-rich disaster-causing structures can be delineated. For example, areas with high permeability coefficients may indicate the presence of karst caves or water-conducting faults; while areas with low permeability coefficients may correspond to intact bedrock. Through comprehensive analysis of multiple parameters, hidden underground disaster sources can be identified more accurately, providing a scientific basis for safe construction.

[0047] On the other hand, an alternative implementation of an observation data preprocessing method is described, which employs a feature extraction algorithm to improve the quality of the joint inversion input data and address issues such as noise interference, baseline drift, and incomplete tracer curves in field experimental data. Accordingly, this method can be performed using the following steps:

[0048] Step 201: Extract features from the tracer penetration curve to generate a travel time feature set, specifically including:

[0049] Baseline correction was performed on the original tracer penetration curve using a multivariate Gaussian function model C(t) = ∑A _i exp(-(t-μ _i ) 2 / (2σ i 2 The net concentration curve was obtained by fitting the data.

[0050] In actual monitoring, background noise or instrument drift may cause tracer concentration readings to be non-zero. Therefore, it is necessary to subtract the background concentration and perform baseline correction. Furthermore, a multivariate Gaussian function model is used to fit the corrected curve to eliminate random fluctuations and extract physical trends.

[0051] Where C(t) represents the observed concentration, which is the superposition of multiple Gaussian components; A _i μ represents the amplitude of the i-th Gaussian component, reflecting the principal concentration of the solute along this path; _i The center time (mean) of the i-th Gaussian component reflects the average time of arrival; σ i The standard deviation of the i-th Gaussian component reflects the dispersion or bandwidth of the solute. By fitting these parameters using the nonlinear least squares method, a physically meaningful net concentration curve can be obtained, eliminating high-frequency noise.

[0052] Alternatively, the multivariate Gaussian function model can also be expressed as:

[0053] C(t) = ∑ i=1 n A _i exp(-(t-μ _i ) 2 / (2σ i 2 ));

[0054] Where C(t) is the fitted net concentration at time t; ∑ is the summation symbol; i is the index of the Gaussian component; n is the total number of Gaussian components; A _i Let be the magnitude of the i-th Gaussian component; exp be the exponential function; t be the time variable; μ _i σ is the center time (mean) of the i-th Gaussian component; i Let be the standard deviation of the i-th Gaussian component.

[0055] Furthermore, Monte Carlo analysis can be incorporated. Based on the standard deviation of measurement errors, random perturbations are added to the fitted curve to generate multiple sets of pseudo-observation data. By repeatedly extracting features from the pseudo-data, the standard deviation of the travel time for each feature can be statistically obtained, providing a statistical basis for subsequent calculation of data weights.

[0056] Step 202: Calculate the normalized cumulative mass recovery function F(t) based on the fitted net concentration curve. This function represents the proportion of the cumulative recovered tracer mass to the total recovered mass as of time t.

[0057] After obtaining the net concentration curve, the normalized cumulative mass recovery function is calculated through numerical integration. This function represents the ratio of the total mass of tracer flowing through the monitoring section before time t to the total mass flowing through the entire observation period. The specific calculation process can be expressed as F(t) = (∫0^t)^t. t C(Э)QdЭ) / (∫0 ∞ C(Э)QdЭ);

[0058] Where F(t) is the normalized cumulative mass recovery function at time t; ∫ is the integral sign; 0 is the lower limit of integration; t is the upper limit of integration (current time); ∞ is infinity; C(Э) is the net concentration of the integral variable at time Э; Q is the flow rate at the monitoring point; and dЭ is the integral differential term.

[0059] Under the assumption of a constant flow rate, this formula simplifies to the ratio of the time integrals of the concentration curve. F(t) is a function that monotonically increases from 0 to 1, shielding the influence of the absolute concentration amplitude and focusing more on the time distribution characteristics of solute transport, thus exhibiting high sensitivity for identifying dominant channels in porous media.

[0060] Step 203: Based on the preset mass recovery rate percentage p%, solve the equation F(t) _p% )=p / 100 determines the characteristic travel time t _p% The travel time features corresponding to different percentages are combined into a travel time feature set.

[0061] Where F is the normalized cumulative mass recovery function; t _p% This is the characteristic travel time corresponding to a mass recovery rate of p%; p is a preset mass recovery rate percentage value, such as 25, 50, 75; 100 is the denominator constant of the percentage.

[0062] To characterize the solute transport process, this step selects multiple percentage nodes to extract feature times. For example, p is set to 25, 50, 75, etc. Solving the equation involves finding the x-axis time corresponding to p / 100 on the cumulative mass curve. Here, t... _25% Representing early arrival time, it is highly sensitive to fast-moving channels with high permeability, such as karst conduits; t _50% When representing median travel, it reflects the average transport velocity of the medium; t _75% Representing late arrival times, these times better reflect the retention and back-diffusion effects within the matrix pores. Based on this, combining characteristic times into a travel-time feature set can provide multi-dimensional constraint information for inversion.

[0063] Step 204: When the acquired tracer penetration curve is an incomplete curve in which no concentration peak was detected, feature extraction also includes derivative peak time transformation processing, specifically:

[0064] Calculate the time derivative curve of a non-holonomic curve and identify the peak time t' of the derivative curve when it reaches its maximum value. _peak ;

[0065] Based on the analytical properties of the convection-diffusion equation, the transformation formula t is used. _peak =t' _peak ×φ estimates the time to peak concentration t _peak This will be added to the travel feature set.

[0066] Where φ is a transformation factor that depends on the Pecklet number. It is used to map the peak time of the derivative of the incomplete penetration curve to an equivalent travel time characteristic with the same caliber as the calculated travel time.

[0067] Alternatively, based on a pre-defined empirical conversion model or analytical approximation model, such as based on the analytical properties of the convection-diffusion equation, the concentration peak time can be estimated using the conversion formula and then added to the travel time feature set.

[0068] In field tests, due to time constraints or equipment malfunctions, a complete tracer curve may sometimes be unavailable, making it impossible to directly read the concentration peak time. To address this, this embodiment provides a remedial method based on derivative characteristics. Accordingly, the acquired incomplete curve (typically containing an ascending segment) is numerically differentiated to obtain the time derivative curve dC / dt. Since the moment of fastest concentration increase (inflection point) often precedes the peak time, the peak time t' of the derivative can usually be identified earlier. _peak Based on the analytical solution derivation of the one-dimensional convection-diffusion equation, it is known that there is a definite proportional relationship between the peak time of the derivative and the peak time of the concentration. The peak time of the concentration can be derived by using the conversion formula.

[0069] Optionally, the conversion formula is as follows:

[0070] t _peak =t' _peak ×φ;

[0071] Among them, t _peak The estimated peak concentration time; t' _peak φ represents the peak time of the identified derivative curve; φ is the conversion factor from the peak time of the derivative to the peak time of the concentration. The specific calculation of the conversion factor can be described as follows:

[0072] φ=sqrt(1+(1+4Pe -1 ) / 2);

[0073] Where φ is the conversion factor; sqrt is the square root symbol; 1, 2, and 4 are all constants; Pe is the Peckley number, representing the ratio of convection rate to diffusion rate; -1 This indicates the reciprocal operation.

[0074] In other words, Pe is the Peklay number, defined as Pe = vL / D _L Where v is the flow velocity, L is the transport distance, and D is the distance traveled. _L φ is the longitudinal dispersion coefficient. This factor reflects the relative strength of convection and dispersion: when convection dominates, i.e., Pe is large, φ approaches 1; when dispersion dominates, φ is greater than 1. This method allows the use of incomplete data, avoiding information waste.

[0075] On the other hand, the specific implementation process of constructing the joint objective function and allocating weights is described, especially the mathematical construction in the joint inversion, including the definition of the unified parameter state vector, the specific form of the objective function, and the adaptive weight allocation strategy based on data quality and sensitivity. This establishes the mathematical foundation for parameter optimization. In this embodiment, it specifically includes:

[0076] Step 301: The unified parameter state vector consists of the logarithmic values ​​of the permeability coefficient, the logarithmic value of the water storage coefficient, and the original value of the effective porosity of each grid cell in the target area.

[0077] This embodiment performs specific transformations on the parameters. Specifically, for the permeability coefficient and storage coefficient, since their values ​​often span multiple orders of magnitude and must be positive, a logarithmic transformation is used as the variable to be inverted. This not only ensures that the parameters after the inverse transformation are always positive, but also makes the parameter distribution closer to a normal distribution, which is beneficial to the convergence of the optimization algorithm. For the effective porosity, since its value range is usually between 0 and 1 and the variation is relatively small, its original value is kept as the variable to be inverted.

[0078] After discretizing the target region into M grid cells, the unified parameter state vector θ can be expressed as θ=[lnK _1 , ..., lnK _M lnS _1 , ..., lnS _M n _1 , ..., n _M ] T Its total dimension is 3M. ​​Where θ is the unified parameter state vector; ln is the natural logarithm function; K _1 , ..., K _M S represents the permeability coefficient of the first to Mth grid cells within the target area; _1 S _M n represents the water storage coefficient of the 1st to Mth grid cells within the target area. _1 , ..., n _M The effective porosity of the first to Mth grid cells within the target area; T This is the matrix transpose symbol.

[0079] Step 302: The joint objective function is constructed as the sum of the weighted sum of squared residuals of the tracer observation data, the weighted sum of squared residuals of the head observation data, and the spatial adaptive regularization term.

[0080] In some embodiments, the joint objective function Φ(θ) is specifically constructed according to the following formula:

[0081] Φ(θ)=∑ i=1 N [w _i_T (T _i_obs -T _i_cal (θ))] 2 +∑ j=1 P [w _j_H (H _j_obs -H _j_cal (θ))] 2 +R(θ);

[0082] Where θ is the unified parameter state vector; N is the total number of tracer travel-time data; T _i_obs T represents the i-th tracer travel time observation; _i_cal (θ) is the calculated value of the i-th tracer travel time based on the current parameter θ; w _i_T The weight of the i-th tracer travel time data; P is the total number of head observation data; H _j_obs H is the j-th head observation value; _j_cal (θ) represents the j-th head calculated based on the current parameter θ; w _j_H Let θ be the weight of the j-th head data; R(θ) is the spatial adaptive regularization term.

[0083] In other words, w _i_T w _j_H All can be determined based on the sensitivity matrix.

[0084] On the other hand, the joint objective function can also be described as:

[0085] Φ(θ)=∑ i=1 N [w i T (T i obs -T i cal (θ))] 2 +∑ j=1 P [w j H (H j obs -H j cal (θ))]2 +R(θ);

[0086] Where Φ(θ) is the joint objective function; ∑ is the summation symbol; i is the index of the tracer data, j is the index of the head data; N is the total number of tracer travel time data; w i T T represents the weight of the i-th tracer travel time data; i obs T represents the i-th tracer travel time observation; i cal (θ) represents the calculated travel time of the i-th tracer based on parameter θ; P is the total number of head observation data; w j H H represents the weight of the j-th head data point; j obs H is the j-th head observation value; j cal R(θ) is the calculated j-th head value based on parameter θ; R(θ) is the spatial adaptive regularization term.

[0087] This step presents the overall optimization objective for the joint inversion. The first term in the formula, ∑... i=1 N [w _i_T (T _i_obs -T _i_cal (θ))] 2 (or rather, ∑) i=1 N [w i T (T i obs -T i cal (θ))] 2 ), is the weighted sum of squared residuals of the tracer travel time data, used to constrain the combination of parameters related to solute transport, mainly permeability coefficient and effective porosity;

[0088] Second term ∑ j=1 P [w _j_H (H _j_obs -H _j_cal (θ))] 2 Or, in other words, ∑ j=1 P [w j H (H j obs -H j cal (θ))] 2, is the weighted sum of squared residuals of the head data, used to constrain the combination of parameters related to pressure transmission, mainly the permeability coefficient and storage coefficient. The two are fused additively to achieve complementarity of multiphysics information.

[0089] The third term, R(θ), is a regularization term used to handle the ill-posedness of the inversion and prevent overfitting.

[0090] Step 303: The data weights in the weighted residual sum of squares are determined based on the measurement error of the corresponding observation data and the Frobenius norm of the corresponding row vector in the sensitivity matrix.

[0091] In some embodiments, the weights of the tracer travel time data and the head data are specifically calculated using the sensitivity matrix and its corresponding measurement error, according to the following formula:

[0092] w _i_T =(||J _i_T || _F ) / (σ T 2 ×∑ k=1 N ||J _k_T || _F ), w _j_H =(||J _j_H || _F ) / (σ H 2 ×∑ k=1 P ||J _k_H || _F );

[0093] Among them, ||…|| _F J represents the Frobenius norm; _i_T σ is the row vector in the sensitivity matrix corresponding to the i-th tracer travel time data; T J represents the standard deviation of the measurement error for tracer travel time data. _j_H σ is the row vector corresponding to the j-th head data point in the sensitivity matrix; H Let be the standard deviation of the measurement error of the head data, and k be the summation index variable. We iterate through N tracer travel time data and P head observation data respectively.

[0094] In other words, the weights of the tracking travel time data can also be:

[0095] w _i_T =(||J _i_T || _F ) / (σ T 2 ×∑ k=1 N ||J_k_T || _F );

[0096] Among them, w _i_T The weight of the i-th tracer travel time data; ||…|| _F J is the Frobenius norm notation, representing the root of the sum of squares of the elements of a vector; _i_T σ is the row vector in the sensitivity matrix corresponding to the i-th tracer travel time data; T ∑ represents the standard deviation of the measurement error for the tracer travel time data; ∑ is the summation sign; k is the summation index variable; N is the total number of tracer travel time data; J _k_T is the row vector in the sensitivity matrix corresponding to the k-th tracer travel time data.

[0097] Alternatively, the weights of the head data can also be:

[0098] w _j_H =(||J _j_H || _F ) / (σ H 2 ×∑ k=1 P ||J _k_H || _F );

[0099] Among them, w _j_H The weight of the j-th head data; ||…|| _F J is the Frobenius norm notation; _j_H σ is the row vector corresponding to the j-th head data point in the sensitivity matrix; H J represents the standard deviation of the measurement error of the head data; P represents the total number of head data points; J represents the total number of head data points. _k_H This is the row vector corresponding to the k-th head data point in the sensitivity matrix.

[0100] This embodiment employs a weighted strategy that comprehensively considers data quality (measurement error) and information content (sensitivity) to allocate the contribution of different observational data in the inversion process. Specifically, the smaller the standard deviation of the measurement error of the tracer data and head data, the more accurate the data, and the greater its weight should be, reflecting reliance on high-quality data. _i_T || _F and ||J _j_H || _F Let Fi and Fj represent the Frobenius norms of the row vectors corresponding to the i-th travel time data and the j-th head data in the sensitivity matrix, respectively; that is, the roots of the sum of squares of the row vector elements. The larger the norm, the more sensitive the observed data is to changes in the model parameters, and the richer the constraint information it contains; therefore, the weights should be increased accordingly.

[0101] The denominator term serves a normalization function, making the weights of different types of data comparable in magnitude. For example, if a set of tracer data has a small measurement error and is highly sensitive to the penetration coefficient of key areas, this data will receive a larger weight, dominating the minimization of the objective function and forcing the model parameters to preferentially fit this data point. Conversely, if a data point has a large measurement error or is not very sensitive to parameters, its weight will be automatically reduced to decrease the interference of noise on the inversion results.

[0102] One example provides an optional implementation of a spatially adaptive regularization mechanism based on a sensitivity field. In traditional geophysical or hydrogeological inversions, regularization typically imposes a uniform smoothing constraint across the entire space, which can unnecessarily reduce the resolution of areas with dense ray coverage (information-rich areas), while ray-blank areas (information-poor areas) may generate false anomalies. In this embodiment, a dynamic weight matrix linked to the sensitivity distribution is constructed to achieve spatially differentiated control of the regularization intensity. Accordingly, it includes:

[0103] Step 401: The spatial adaptive regularization term is constructed based on the regularization weight matrix. The regularization weight matrix is ​​determined according to the comprehensive sensitivity field of each grid cell in the target area. The comprehensive sensitivity field consists of the comprehensive sensitivity index of each grid cell to the travel time feature set and the head change curve.

[0104] In other embodiments, the spatial adaptive regularization term R(θ) is constructed based on the first-order difference operator matrix and the iteratively updated regularization weight matrix, and the specific calculation formula is as follows:

[0105] R(θ)=θ T ×Λ _k ×L T ×L×θ;

[0106] Where θ is the unified parameter state vector; L is the first-order difference operator matrix, used to calculate the spatial gradient of the model parameters; T For transpose operation; Λ _k is the spatial adaptive regularization weight matrix in the k-th iteration step, used to weight the smoothing constraint strength at different spatial locations.

[0107] Or rather, Λ _k The spatial adaptive regularization weight matrix is ​​calculated dynamically based on the sensitivity matrix in the k-th iteration step and is used to weight the smoothing constraint strength at different spatial locations.

[0108] In this embodiment, the regularization term is used to suppress drastic spatial oscillations of the model parameters, ensuring the smoothness and stability of the inversion results. θ is the unified parameter state vector. L is the first-order difference operator matrix, a sparse matrix used to calculate the gradient of the model parameters in space, i.e., the difference in parameter values ​​between adjacent grid cells. T Let L be the transpose of L.

[0109] Furthermore, the formula actually defines a weighted seminorm, which measures the roughness of the parametric field in space. Where Λ _k This refers to the spatial adaptive regularization weight matrix in the k-th iteration. Unlike the traditional identity matrix, it is a diagonal matrix, with non-uniform values ​​on its diagonal, used to weight the smoothness constraint strength at different spatial locations. If the weight at a certain location is large, the smoothness of the parameters at that location and its neighborhood is subject to stronger constraints; conversely, if the weight is small, larger gradient changes in the parameters at that location are allowed, preserving high-frequency details.

[0110] Step 402, the comprehensive sensitivity index of the grid cell is defined as: the square root of the sum of the squares of the products of the partial derivatives of all travel time calculations and the corresponding travel time data weights for the grid cell, plus the sum of the squares of the products of the partial derivatives of all head calculations and the corresponding head data weights for the grid cell.

[0111] In other embodiments, the construction of the regularization weight matrix depends on the comprehensive sensitivity index of each grid cell within the target region; the comprehensive sensitivity index ξ of the m-th grid cell in the k-th iteration step. _k_m The specific calculation is performed according to the following formula:

[0112] ξ _k_m =sqrt(∑ i=1 N [( T _i_cal ) / ( θ _m )×w _i_T ] 2 +∑ j=1 P [( H _j_cal ) / ( θ _m )×w _j_H ] 2 );

[0113] Where, θ _m This represents the parameter component corresponding to the m-th grid cell in the unified parameter state vector; T _i_cal ) / ( θ_m () represents the partial derivative of the calculated value of the i-th travel time with respect to the parameter of the m-th grid cell; H _j_cal ) / ( θ _m ) represents the partial derivative of the j-th calculated head with respect to the m-th grid element parameter; w _i_T and w _j_H These are the corresponding weights.

[0114] To determine the spatial distribution of the regularization weights, it is necessary to quantify the ability of the observed data to constrain the parameters of each grid cell, i.e., to calculate the comprehensive sensitivity index. Wherein, θ _m This represents the parameter component corresponding to the m-th grid cell in the unified parameter state vector, such as the logarithmic value of the permeability coefficient of that cell. T _i_cal ) / ( θ _m The partial derivative of the calculated value of the i-th travel time with respect to the parameter of the m-th grid cell reflects the degree of influence of the change in the cell parameter on the travel time data; H _j_cal ) / ( θ _m ) represents the partial derivative of the j-th calculated head value with respect to the m-th grid element parameter, reflecting the impact on the head data. _i_T and w _j_H These are the corresponding data weights.

[0115] The formula uses the root mean square form to weighted couple the sensitivity information of the two physical fields (solute transport field and groundwater flow field). If a grid cell is located where multiple tracer rays pass through and is close to the head monitoring point, its corresponding partial derivative absolute value is large, and the calculated comprehensive sensitivity index will also be large, indicating that the cell is sufficiently constrained by data; conversely, if the cell is located in the monitoring blind zone, the comprehensive sensitivity index will approach zero.

[0116] Step 403: The regularization weight matrix is ​​a diagonal matrix, and the values ​​of its diagonal elements are negatively correlated with the comprehensive sensitivity index of the corresponding grid cell. This results in grid cells with lower comprehensive sensitivity indices having larger regularization weights, while grid cells with higher comprehensive sensitivity indices have smaller regularization weights.

[0117] In other embodiments, the regularization weight matrix is ​​a diagonal matrix, whose m-th diagonal element λ _mm_k The value is determined based on the normalized value of the comprehensive sensitivity index, and the specific calculation formula is as follows:

[0118] ξ _m_k '=(ξ _m_k -ξ _min_k ) / (ξ_max_k -ξ _min_k );

[0119] λ _mm_k =λ _0_k ×[1+γ×(1 / ξ _m_k ') β ];

[0120] Where, ξ _m_k ' represents the normalized relative sensitivity, also known as the normalized sensitivity; ξ _max_k and ξ _min_k λ represents the maximum and minimum values ​​of the comprehensive sensitivity index of all grid cells in the current iteration step, respectively. _0_k γ is the baseline regularization parameter for the current iteration step; γ is the regularization adjustment coefficient; β is the shape control parameter.

[0121] In other words, λ _0_k The baseline regularization parameter for the current iteration step is obtained by adaptive calculation based on the number of iterations.

[0122] In this step, a mapping rule from sensitivity to regularization weights is established. Accordingly, to eliminate the influence of numerical dimensions, the normalized relative sensitivity is calculated using the maximum and minimum values ​​of the comprehensive sensitivity index of all grid cells in the current iteration step, with a value ranging from 0 to 1.

[0123] Furthermore, regularized weights are constructed using an inverse power-law function. This formula is designed based on the adaptive logic of a weak model with sufficient data and a strong model with insufficient data. Specifically, when the normalized relative sensitivity is close to 1, i.e., in the high-sensitivity region, (1 / (ξ)) _m_k ') β When the regularization weight is close to 1, it approaches the baseline value. At this point, the inversion is mainly driven by the data residual term, which can characterize the boundary of geological anomalies. When the normalized relative sensitivity is close to 0, i.e., it is located in a low-sensitivity zone or a blind zone, (1 / (ξ)) _m_k ') β The sharp increase leads to a significant increase in regularization weights. At this point, the inversion is mainly driven by the smoothing constraint term, which forces the parameters in this region to transition smoothly and avoids artifacts.

[0124] According to one aspect of this application, a specific numerical calculation example is provided.

[0125] Suppose that in a certain iteration, the shape control parameter β is set to 1, the regularization adjustment coefficient γ is set to 1, and the baseline regularization parameter is set to 10. For mesh cell A located at the dense intersection of rays, its normalization sensitivity is calculated to be 0.8, then its regularization weight = 10 × [1 + 1 × (1 / 0.8)^1] = 22.5.

[0126] For grid cell B located at the edge of the ray coverage, its normalization sensitivity is only 0.05, so its regularization weight = 10 × [1 + 1 × (1 / 0.05)^1] = 210. It can be seen that cell B, which has weaker data constraints, is subject to a regularization intensity nearly 10 times that of cell A, suppressing the uncertainty of the blind zone.

[0127] Optionally, the shape control parameter β can also be set to 2 or greater to further widen the difference in regularization intensity between the high and low sensitivity regions.

[0128] Step 404: The baseline regularization parameter in the regularization weight matrix adopts an iterative adaptive decay strategy, and the value of the baseline regularization parameter gradually decreases as the number of iterations increases.

[0129] Alternatively, the baseline regularization parameter adopts an adaptive update strategy that decays with the number of iterations, and the specific calculation formula is as follows:

[0130] λ _0_k =λ _base ×[η _0 ×exp(-k / τ)+η _∞ ];

[0131] Where, λ _base The initial base regularization strength is set; k is the current iteration number; τ is the decay time constant; η _0 η is the initial decay factor; _∞ exp(...) is the asymptotic decay factor, and exp(...) is the exponential function.

[0132] In addition to spatial adaptation, this embodiment also introduces a temporal (iteration dimension) adaptation strategy. In the early stages of inversion, since the initial model may be far from the true solution and the sensitivity calculation based on the initial parameter field may be biased, a larger regularization strength is needed to stabilize the inversion path and prevent it from getting trapped in local minima.

[0133] As the number of iterations k increases, the parameter field gradually approximates the true distribution. At this point, the regularization constraint should be gradually relaxed to extract high-frequency information from the data. Wherein, λ _base η is the initial set of base regularization strength; τ is the decay time constant, controlling the rate of decay; _0 η is the initial decay factor; _∞ This is an asymptotic decay factor, ensuring that the regularization parameter does not drop to zero in the later stages of iteration, thus maintaining basic numerical stability. For example, setting λ... _base =100, η _0 =0.9, η _∞ =0.1, τ=5. In the first iteration, λ _0_1 ≈83.7; while in the 20th iteration, λ_0_20 ≈11.6.

[0134] It should be understood that the above dynamic decay strategy improves the convergence effect of complex structure inversion.

[0135] Optionally, the shape control parameter β ranges from 1 to 3, and can be set to 2. When β is small, such as β=1, the difference in regularization intensity between high-sensitivity and low-sensitivity areas is relatively gentle, suitable for scenarios with relatively uniform geological conditions; when β is large, such as β=3, the differentiation effect is greater, suitable for scenarios with obvious geological anomalies. The regularization adjustment coefficient γ ranges from 0.5 to 2, and can be set to 1. This parameter controls the amplification factor of the sensitivity feedback. The decay time constant τ is determined based on 1 / 4 to 1 / 3 of the expected total number of iterations. For example, when the expected total number of iterations is 20, τ can be set to 5 to 7.

[0136] Another example describes an alternative technical approach to resolving cyclic dependencies using an alternating prediction-correction update strategy, applicable to hydrogeological inversion processes. Specifically, this includes:

[0137] Step 501: The slowness parameter in the forward model depends on the hydraulic gradient field. The iterative solution process of the joint objective function adopts a prediction-correction alternating update strategy to solve the cyclic dependency between the slowness parameter and the permeability coefficient to be inverted.

[0138] In the theory of solute transport in porous media, the slowness *s* (the reciprocal of velocity) depends not only on the medium properties (effective porosity *n*, permeability coefficient *K*) but also on the external dynamic field (hydraulic gradient |▽H|). The specific physical relationship is *s* = *n* / (*K* × |▽H|). This also constitutes a logical deadlock: to invert *K*, the travel time residual needs to be calculated, and the travel time calculation depends on the slowness *s*; however, the calculation of the slowness *s* depends on the hydraulic gradient |▽H|; and the hydraulic gradient |▽H| itself is a function of *K*, with the flow field determined by the distribution of *K*. If |▽H| from the previous iteration is directly used or assumed to be constant in the iteration, a large nonlinear error will be introduced. Therefore, this embodiment designs a prediction-correction mechanism embedded in the inversion loop to decouple the flow field update from the parameter inversion.

[0139] Step 502, the prediction-correction alternating update strategy executes the following steps sequentially in each iteration:

[0140] Perform the prediction step, which is based on the penetration coefficient field K of the current iteration step. _k Solve the steady-state groundwater flow control equation ▽×(K) _k ▽H)=0, obtain the current head field H _k and hydraulic gradient field |▽H| _kIn other words, this method is based on the assumption of an equivalent porous medium model and is applicable to karst systems with well-developed fracture zones and diffuse flow-dominated environments.

[0141] Perform a correction step, which is based on the hydraulic gradient field, the permeability coefficient field of the current iteration step, and the effective porosity field n. _k Update the slowness value s of each grid cell. _m_k The specific calculation formula is as follows:

[0142] s _m_k =n _m_k / (K _m_k ×|▽H| _m_k );

[0143] Perform the inversion step, which involves recalculating the tracer travel time and its sensitivity matrix based on the updated slowness value, and performing optimization for the joint objective function to obtain the updated unified parameter state vector.

[0144] Where ▽ is the Hamiltonian operator; H is the head; |▽H| _m_k K represents the hydraulic gradient magnitude at the m-th grid cell in the k-th iteration. _m_k and n _m_k These are the permeability coefficient and effective porosity of the m-th grid cell in the k-th iteration, respectively.

[0145] In other words, the permeability field and effective porosity field of the current iteration step are extracted from the unified parameter state vector.

[0146] In the prediction phase, the permeability coefficient is temporarily frozen and treated as a known quantity, which is then substituted into the groundwater flow governing equation. For the steady-state flow case, the elliptic partial differential equation ▽×(K) is solved. _k ▽H)=0; for unsteady flow, the parabolic equation is solved by combining the time step. The head distribution of the entire field is obtained through numerical solution, such as the finite element method, and the hydraulic gradient modulus at the center of each grid cell is further calculated by numerical differentiation.

[0147] Furthermore, by utilizing the predicted updated hydraulic gradient field, combined with the current K... _k and n _k The slowness distribution of the entire field is recalculated according to the physical definition, ensuring that the slowness field is physically consistent with the current flow field and medium parameter field.

[0148] Based on this, using the corrected slowness field, a solute transport forward model (or ray tracing model) is run to calculate the tracer travel time and update the sensitivity matrix, since the sensitivity of travel time to K is directly related to slowness. The linearized objective function is constructed and minimized, and the parameter update amount is solved to obtain the parameters for the next iteration.

[0149] It should be understood that this step, through the PCI (Predict-Correct-Invert) loop, eliminates the computational bias caused by nonlinear coupling, ensuring the robust convergence of the inversion process.

[0150] Another example provides an alternative implementation of forward modeling based on curved ray tracing, suitable for highly heterogeneous media. In the presence of high-permeability structures such as karst conduits or faults, tracers tend to migrate along the path of least resistance, which deviates significantly from a straight line. Using curved ray tracing technology instead of the traditional straight-line approximation can improve the inversion accuracy in high-contrast media. Specifically, this approach includes:

[0151] Step 601, when calculating the tracer travel time in the inversion step, uses a curved ray tracing method based on Fermat's shortest time principle, specifically including:

[0152] Based on the current slowness distribution, the fast propagation method is used to solve the process function equation |▽T(x)|=s(x), and the shortest arrival time field of each point in the target area is calculated with the tracer source point as the starting point.

[0153] Starting from the monitoring point, backtrack along the negative gradient direction of the shortest arrival time field to determine the curved ray path connecting the tracer source point and the monitoring point;

[0154] The tracer travel time is calculated based on the sum of the products of the length of the curved ray path within each grid cell and the slowness value of the corresponding grid cell.

[0155] Where |▽T(x)| is the gradient magnitude of the arrival time field T(x) at spatial location x; s(x) is the slowness at spatial location x.

[0156] Accordingly, Fermat's principle states that the actual path of a wave (or solute front) propagating in a medium is the path with the shortest travel time. This step uses the functional equation |▽T(x)|=s(x) to describe this process, where T(x) represents the time it takes for the wavefront to arrive at spatial location x, and s(x) is the slow-motion field. To solve this equation, this embodiment may use the Fast Mover Method (FMM). FMM uses heap sorting to update grid nodes in ascending order of arrival time, obtaining the arrival times of the entire field in a single operation with a computational complexity of O(MlogM).

[0157] Where M corresponds to the total number of parameter grids in the inversion model, such as the number of discrete grid cells for the permeability coefficient; logM usually comes from the time overhead of binary search, balanced binary tree operation, or divide-and-conquer strategy; O(MlogM) corresponds to the representation of the algorithm's time complexity.

[0158] After obtaining the arrival time field, starting from the monitoring point, a streamline backtracking is performed along the direction of the fastest descent of the T(x) gradient, i.e., the direction of -▽T, until returning to the tracer source point. This backtracking trajectory is the curved ray path. Based on this, the travel time can be calculated by performing a line integral of the slowness along this path, which is represented in the discrete grid as the path length multiplied by the sum of the unit slowness. It should be understood that, compared to the linear model, this method can accurately capture the physical phenomenon of the tracer bypassing low-permeability zones and converging in high-permeability channels.

[0159] Step 602, the iterative solution process includes a collaborative update mechanism for the ray path and parameter field, specifically:

[0160] In each iteration or at a preset interval, the curved ray tracing method is re-executed based on the updated parameter distribution to obtain the updated curved ray path. Based on the intercept length of the updated curved ray path in each grid cell, the sensitivity matrix used to construct the joint objective function for the next iteration step is recalculated.

[0161] During the inversion process, as the parameter field is continuously updated, the slowness distribution of the medium also changes, inevitably leading to a change in the shortest propagation path. Therefore, a collaborative update mechanism must be adopted. This embodiment stipulates that, in each iteration of the LM algorithm, or every fixed number of iterations, to save computational resources, the FMM and backtracking algorithms must be rerun using the latest parameter field to refresh the geometric coordinates of the ray path. This mechanism ensures that the sensitivity matrix always reflects the current physical state, avoiding linearization errors introduced by a fixed path.

[0162] Step 603: Recalculate the sensitivity matrix used to construct the joint objective function for the next iteration step. Specifically, calculate the elements of the sensitivity matrix based on the curved ray path, using the following formula:

[0163] T _sr / (lnK _m )=-l _sr_m ×s _m ;

[0164] T _sr / n _m =l _sr_m ×1 / (K _m ×|▽H| _m );

[0165] in, T _sr / (lnK _m )and T _sr / n _m These are the sensitivity coefficients for the logarithmic permeability coefficient and effective porosity of the m-th grid cell during tracer travel, respectively; _sr_m s is the path length of the curved ray within the m-th grid cell; _m K represents the slowness of the m-th grid cell; _m and These are the permeability coefficient and hydraulic gradient modulus of the m-th grid cell, respectively.

[0166] This step presents a specific analytical formula for calculating the sensitivity coefficient based on the curved ray path, which forms the basis for constructing the Jacobian matrix. For the m-th mesh cell, its logarithmic penetration coefficient lnK... _m The change in the travel time T from the source point s to the monitoring point r _sr The effect is caused by the partial derivative ∂T _sr / ∂(lnK _m )describe.

[0167] According to the chain rule, this partial derivative is equal to the intercept length l of the curved ray within the element. _sr_m Multiply by the unit's slowness s _m The negative value of . If a cell is the path of a ray and has a large slowness (slow flow velocity), increasing the permeability coefficient of that cell (i.e. reducing the slowness) will shorten the travel time.

[0168] Similarly, effective porosity n _m The sensitivity coefficient is the path length l _sr_m ×1 / (K _m ×|▽H| _m It should be understood that l _sr_m The path length changes dynamically with each iteration, which reflects the characteristics of nonlinear tomography.

[0169] In some schemes, if the heterogeneity of the geological conditions is predicted to be weak, for example, the difference in permeability coefficients is within one order of magnitude, a straight ray approximation can be used to improve computational efficiency, i.e., assuming l _sr_m The length of the straight line segment connecting the source and the monitoring point within the unit is given. The above formula still applies, but the path length becomes a fixed value.

[0170] This embodiment describes an exemplary scheme for the comprehensive identification and engineering application of water-rich disaster-causing structures, particularly the interpretation of results and engineering applications in the later stages of inversion. This embodiment introduces a comprehensive disaster-causing index and multi-parameter joint identification rules to achieve accurate identification of underground hidden disaster source types.

[0171] Furthermore, the location of water-rich disaster-causing structures is delineated based on the optimized parameter distribution, specifically including:

[0172] Based on the optimized unified parameter state vector, the pressure conductivity distribution of the target area is calculated, where the pressure conductivity is the ratio of the permeability coefficient to the water storage coefficient. The comprehensive disaster index of each grid cell is calculated using the following formula:

[0173] I _haz (x)=α _1 ×K(x) / K _max +α _2 ×n(x) / n _max +α _3 ×D(x) / D _max ;

[0174] Where x represents the spatial location of the grid cell; K(x), n(x), and D(x) are the permeability coefficient, effective porosity, and pressure conductivity at that location, respectively; K _max n _max D _max These represent the maximum values ​​of the corresponding parameters within the target region; α _1 α _2 α _3 The weighting coefficients are preset and satisfy α. _1 +α _2 +α _3 =1;

[0175] Areas with a comprehensive disaster-causing index greater than a preset threshold are designated as water-rich disaster-causing structures, and the type of disaster-causing structure is identified based on the high and low combination characteristics of K(x), n(x), and D(x).

[0176] In other words, based on the optimized unified parameter state vector, the distribution of the pressure conductivity coefficient in the target area is calculated, where the pressure conductivity coefficient is the ratio of the permeability coefficient to the water storage coefficient.

[0177] The comprehensive disaster-causing index of each grid cell is calculated based on the permeability coefficient, effective porosity, and pressure conductivity coefficient.

[0178] Areas with a comprehensive disaster-causing index greater than a preset threshold are designated as water-rich disaster-causing structures. The type of disaster-causing structure is determined based on the combination of high and low permeability coefficient, effective porosity, and pressure conductivity coefficient.

[0179] In this step, the pressure permeability coefficient needs to be calculated. Physically, the pressure permeability coefficient is defined as D=K / S, reflecting the speed at which pressure waves propagate in a medium. Here, K and S correspond to the permeability coefficient and the water storage coefficient, respectively; in other words, they correspond to Kpermeability and Sstorage, respectively. _1 , ..., K _M S _1 S _M .

[0180] In the previous inversion, the distributions of permeability and storage coefficients were obtained, so the pressure conductivity field can be directly obtained through division. The pressure conductivity has a unique indicative role in the differentiation of fissure zones and karst caves.

[0181] Furthermore, this embodiment defines a comprehensive disaster-causing index I. _haz This index is a normalized dimensionless scalar that integrates fluid conductivity (characterized by permeability coefficient), water storage space size (characterized by effective porosity), and pressure transmission efficiency (characterized by pressure conductivity coefficient). It is normalized by dividing by the maximum value of each parameter across the entire field, eliminating differences in the dimensions and orders of magnitude of different physical quantities.

[0182] Furthermore, the weighting coefficient α _1 α _2 α _3 The selection of the permeability coefficient is usually determined based on engineering experience or the analytic hierarchy process (AHP). For example, in projects where the risk of sudden water inrush is the primary concern, the weight α of the permeability coefficient is... _1 It can be set to 0.5, while the weights of effective porosity and pressure conductivity are each set to 0.25.

[0183] Based on the calculated comprehensive disaster-causing index field, an anomaly identification threshold Z needs to be set. _th To delineate the danger zone. Optionally, the entire field I can be calculated first. _haz The mean and standard deviation of Z _th Set to mean + 2 standard deviations or mean + 3 standard deviations. For all I... _haz (x)>Z _th The continuous area is designated as a water-rich disaster-causing structure.

[0184] Furthermore, to guide specific engineering measures such as grouting, drainage, or bypass, this embodiment also provides a disaster-causing structural type identification rule based on the KnD (permeability coefficient - effective porosity - pressure conductivity) parameter combination characteristics. The specific identification logic is as follows:

[0185] If an abnormal area exhibits high K, high n, and high D characteristics, it is identified as a karst conduit or a large cave. This is because caves not only conduct water very quickly (high K and high D), but also have very high effective porosity (high n) due to the presence of cavities.

[0186] If it exhibits high K, medium n, and high D characteristics, it is identified as a fault fracture zone. Fault zones typically have high water conductivity and pressure transmission capacity, but their porosity is mainly composed of the gaps between fractured rocks, which is lower than that of cavities and sinkholes.

[0187] If the rock exhibits high K, low n, and medium D characteristics, it is identified as a zone with densely developed fractures. Although fracture flow has strong water conductivity, the overall porosity of the rock mass is low, and pressure transmission is limited by the connectivity of fractures.

[0188] In some scenarios, for a two-dimensional plane (including thickness), T represents the hydrodynamic conductivity, S represents the water storage capacity, and the pressure conductivity is defined as D. _h =T / S; while the hydraulic conductivity is the product of the permeability coefficient K and the effective thickness b of the aquifer, T=Kb; for a three-dimensional continuous medium, S _s The specific water storage capacity is represented by the pressure conductivity coefficient D. _h =K / S _s .

[0189] On the other hand, the equation can be further expressed as:

[0190] |▽T(x)|=s(x)=n(x) / (K(x)×|▽H(x)|);

[0191] In the formula, ▽ is the gradient symbol, |...| represents the magnitude, and ▽H(x) corresponds to the gradient vector of the head at position x.

[0192] On the other hand, to quantify the accuracy of the inversion results, the following technical indicators can be used:

[0193] Parameter reconstruction accuracy is calculated as the root mean square error (RMSE) and correlation coefficient (R) between the inverted parameter field and the true parameter field (if validation data is available). The smaller the RMSE and the closer R is to 1, the higher the inversion accuracy.

[0194] Goodness of fit is the normalized root mean square error (NRMSE) of the tracer travel-time residuals after the inversion is completed. _T Normalized root mean square error (NRMSE) of head residuals _H The possible convergence criterion is NRMSE. _T <0.1 and NRMSE _H <0.05.

[0195] Resolution evaluation, based on checkerboard testing or resolution matrix analysis, assesses the ability of inversion to identify geological anomalies at different scales, requiring the ability to distinguish anomalies that are 2 to 3 times larger than the grid cells.

[0196] According to another aspect of this application, the quality of the initial parameter guesses affects the convergence result during iterative inversion. Accordingly, if the initial values ​​are too far from the true solution, the LM algorithm may get trapped in local minima. To address this, an initialization strategy based on serial inversion is provided. Specifically, before initiating joint synchronous inversion, a simplified tomography is performed separately using tracer travel time data.

[0197] In this preprocessing stage, it is assumed that the ray is a straight line, the head data is ignored, and only the linear equation system ΔT=L is solved. *×Δs quickly yields a coarse distribution of slowness and permeability. This coarse distribution is then used as the initial value for the joint inversion. This strategy of coarse-to-fine and sequential-to-parallel approaches leverages the stability of linear inversion to provide a foundation for nonlinear inversion.

[0198] Or, in other words, in the formula, ΔT represents the deviation vector for tracking travel time, and L * Δs represents the ray length matrix, and Δs represents the deviation vector of the characteristic parameter to be determined.

[0199] According to another aspect of this application, for the curved ray tracing method, in some scenarios with relatively simple geological conditions, such as when the strata are approximately horizontal, the heterogeneity is weak, and the difference in permeability coefficient is less than one order of magnitude, a straight ray can be used as an optional implementation method.

[0200] In this scenario, the forward model assumes that the tracer propagates along a straight line connecting the source and monitoring points, with the path length being a fixed geometric constant that is no longer updated with parameter iterations. Although its accuracy is slightly lower than the curved ray method, this approach reduces computational load and is suitable for rapid preliminary investigations in engineering settings.

[0201] In some embodiments, the iterative solution process also includes convergence determination and exception handling mechanisms.

[0202] For convergence determination, the following two termination conditions can be set:

[0203] When the relative rate of decrease of the objective function in two consecutive iterations is less than a preset threshold (within the range of 10), -4 Up to 10 -3 When the inversion is complete, it is determined that the inversion has converged.

[0204] The iteration terminates when the number of iterations reaches the preset maximum number of iterations (ranging from 50 to 100).

[0205] For handling data anomalies, the following strategies are included:

[0206] For tracer concentration data, if the measured value at a certain moment exceeds the physically reasonable range, such as a negative value or more than twice the initial concentration, the data point will be marked as abnormal and removed or downweighted when calculating the residual.

[0207] For head data, if the rate of change of head at adjacent sampling times exceeds three times the standard deviation of the normal hydraulic response, an anomaly detection procedure is initiated, and if necessary, adjacent data are used for linear interpolation.

[0208] For the numerical stability of the sensitivity matrix, when the matrix condition number exceeds a preset threshold (which can be 10), 6 When the current iteration step is in the 0-5-fold range, the baseline regularization parameter is automatically increased to suppress solution oscillations caused by ill-conditioned conditions.

[0209] In other embodiments, the tracer travel time calculation is obtained in any of the following ways:

[0210] The penetration curve is obtained by numerical forward modeling based on the convection-diffusion equation and extracted according to the preset travel time feature rules;

[0211] When the Peklay number Pe is not less than the threshold Pe _0 When the convection-diffusion equation is converted into an equation function under the condition of convection dominance, the travel time is calculated.

[0212] When Pe≥Pe _0 When Pe is used, a functional equation is used for approximation; when Pe <Pe _0 Travel time was extracted using convection-diffusion forward modeling.

[0213] Among them, the equation of the process is the convection-dispersion equation in the case of Pe≥Pe _0 An approximate model of arrival time under convection-dominated conditions.

[0214] The conversion factor of the Pecklet number is used to map the peak time of the derivative of the incomplete penetration curve to an equivalent travel time characteristic with the same caliber as the travel time calculation.

[0215] In some other embodiments, parts of the method of the present invention may also be:

[0216] Optionally, the permeability coefficient field K(f{x}), water storage coefficient field S(f{x}), and effective porosity field n(f{x}) to be inverted can be combined into a unified parameter state vector θ=[lnK(x)]. _1 ), ..., lnK(x) _M ), lnS(x _1 ), ..., lnS(x _M ), n(x _1 ), ..., n(x) _M )] T

[0217] Where M is the number of grid cells, and the logarithm of the permeability coefficient and storage coefficient is taken to ensure non-negativity and improve numerical stability, K(x) _M Let be the permeability coefficient at the m-th grid cell, and lnK(x) _M S(x) is the natural logarithm of the permeability coefficient. _M Let ) be the water storage coefficient at the m-th grid cell, lnS(x) _M Similarly, n(x) _M ) represents the effective porosity at the m-th grid cell. This parameter itself has a value range of 0 to 1 and does not require additional logarithmic processing.

[0218] Furthermore, a tracer transport forward model is established, based on the convection-dispersion equation, to describe the spatiotemporal evolution of solute concentration, specifically:

[0219] n× C / t + ▽ × (vC) = ▽ × (D) _h ▽C);

[0220] Where C is the solute concentration, v is the actual flow rate, and D _h Where is the hydrodynamic dispersion coefficient, and n is the effective porosity of the porous medium. C / t represents the rate of change of solute concentration over time, and ∠ is the divergence operator. The travel time T is calculated based on numerical solutions. _cal =f _T (K, n, H).

[0221] Furthermore, a forward model of groundwater flow is established based on the unsteady flow governing equations, specifically expressed as follows:

[0222] S× H / t = ▽ × (K ▽ H) + W;

[0223] Where W represents the source and sink terms, or S represents the water storage coefficient of the porous medium. H / t is the rate of change of water head over time, H is the groundwater head, K is the permeability coefficient of the porous medium, ▽H is the water head gradient, and W is the source and sink term. Positive values ​​indicate the source of groundwater recharge (such as precipitation infiltration, river recharge, etc.), while negative values ​​indicate the discharge term of groundwater (such as artificial extraction, evaporation, etc.).

[0224] The calculated head response H is obtained based on numerical solution. _cal =f _H (K, S).

[0225] Optionally, the sensitivity matrix of the two types of observation data to each parameter can be calculated using either the perturbation method or the adjoint state method. That is, the sensitivity matrix of travel time T to the parameters is:

[0226] J _T =[ T _1 / lnK _1 , ..., T _1 / lnK _M , T _1 / n _1 , ..., T _1 / n _M ;

[0227] T _2 / lnK _1 , ..., T _2 / lnK _M , T _2 / n _1 , ..., T _2 / n _M ;

[0228] …;

[0229] T _N / lnK _1 , ..., T _N / lnK _M , T _N / n _1 , ..., T _N / n _M ];

[0230] Furthermore, the sensitivity matrix of the head H to the parameters can be expressed as:

[0231] J _H =[ H _1 / lnK _1 , ..., H _1 / lnK _M , H _1 / S _1 , ..., H _1 / S _M ;

[0232] H _2 / lnK _1 , ..., H_2 / lnK _M , H _2 / S _1 , ..., H _2 / S _M ;

[0233] …;

[0234] H _P / lnK _1 , ..., H _P / lnK _M , H _P / S _1 , ..., H _P / S _M ];

[0235] Where N is the number of travel time data points and P is the number of head observation data points.

[0236] Optionally, the joint objective function can be defined as the weighted sum of squares of the residuals of the two classes of data, i.e.:

[0237] Φ(θ)=∑ i=1 N [w _i_T (T _i_obs -T _i_cal (θ))] 2 +∑ j=1 P [w _j_H (H _j_obs -H _j_cal (θ))] 2 +λ×R(θ);

[0238] Where R(θ) is the regularization term and λ is the regularization parameter.

[0239] Furthermore, the weights of the travel time data are determined based on sensitivity and measurement error, namely:

[0240] w _i_T =(||J _i_T || _F ) / (σ T 2 ×∑ k=1 N ||J_k_T || _F );

[0241] Water head data weights, i.e.:

[0242] w _j_H =(||J _j_H || _F ) / (σ H 2 ×∑ k=1 P ||J _k_H || _F );

[0243] This weighting method gives greater weight to highly sensitive observation data while taking into account the impact of measurement accuracy.

[0244] Optionally, the L-Mt method can be used for iterative solution. Let the parameter vector of the k-th iteration be θ. _k The parameter update amount is:

[0245] Δθ _k =(J T ×W×J+μ _k ×I+λ×L T ×L) -1 ×J T ×W×δ _k ;

[0246] Where J is the Jacobian matrix, W is the weight matrix, and μ is a diagonal matrix. _k Let be the damping coefficient of the (LM) method in the k-th iteration, where k corresponds to the iteration number, I is the identity matrix, λ is the regularization coefficient (a global scalar), and L is the regularization matrix, typically a Laplace matrix or the identity matrix. -1 The matrix inversion operator, δ _k This is the residual vector for the k-th iteration, used to drive the iterative update of parameters.

[0247] Where the joint sensitivity matrix J = [(J _T J _H )] T ;

[0248] The weight matrix is ​​[w _1_T , ..., w _N_T ;w _j_H , ..., w _P_T ] T ; Corresponding to the data weights of the first N and first P times during the travel;

[0249] residual vector δ _k =[T _1_obs -T _1_cal (θ_k ), ..., T _N_obs -T _N_cal (θ _k ), H _1_obs -H _1_cal (θ _k ), ...,

[0250] (H _P_obs -H _P_cal (θ _k )];

[0251] Damping factor μ _k Adjust according to the adaptive strategy; that is, if the objective function decreases in this iteration, then μ... _k +1=μ _k / β; otherwise μ _k +1=μ _k ×β, where β>1 is an adjustment coefficient.

[0252] Perform the parameter update, i.e.:

[0253] Updated unified parameter state vector θ _k+1 =θ _k +Δθ _k ;

[0254] Optionally, iteration may stop when any of the following conditions are met:

[0255] The relative change in the objective function is less than the threshold;

[0256] The parameter change is less than the threshold.

[0257] Reaching the maximum number of iterations;

[0258] Output the permeability coefficient field, water storage coefficient field, and effective porosity field, and calculate the pressure conductivity coefficient.

[0259] Furthermore, based on the permeability field, the spatial anomaly of the logarithmic permeability coefficient is calculated, i.e.:

[0260] Z _K (x)=(lnK(x)-(lnK) ‾ ) / σ _lnK ;

[0261] Wherein, the spatial anomaly of the logarithmic permeability coefficient at spatial location x is given by lnK(x), and the natural logarithm of the permeability coefficient K(x) at spatial location x is given by (lnK). ‾ σ _lnK These are the mean and standard deviation of the logarithmic permeability coefficient, respectively.

[0262] Will satisfy Z _K (f{x})>Z _thThe area is designated as a hyperosmolarity zone, where Z _th This is the abnormal threshold, typically set to 2.0~3.0.

[0263] Furthermore, the comprehensive disaster-causing index is defined as:

[0264] I _haz (x)=α _1 ×(K(x)) / K _max +α _2 ×(n(x)) / n _max +α _3 ×(D(x)) / D _max ;

[0265] Where, α _1 +α _2 +α _3 =1 is the weighting coefficient, which is determined based on engineering experience or the analytic hierarchy process.

[0266] Based on the spatial distribution characteristics of I_{haz} and combined with the geological background, the disaster-causing structural type is identified. The specific determination logic is as follows:

[0267] High K, high n, and high D may indicate a cave or karst conduit.

[0268] High K, medium n, and high D may indicate a fault fracture zone;

[0269] High K, low n, and medium D may indicate a zone with densely developed fractures.

[0270] This application employs a multiphysics joint synchronous inversion strategy to construct a unified state vector containing permeability coefficient, storage coefficient, and effective porosity. Tracer and head data are simultaneously utilized through a joint objective function. This achieves data fusion, eliminates the unidirectional error accumulation inherent in traditional two-stage inversion, and leverages the complementary advantages of different data sets.

[0271] Furthermore, a prediction-correction alternating update strategy is introduced. Step-by-step decoupled computation helps achieve logical self-consistency of physical parameters during the inversion process. A curved ray tracing and collaborative update mechanism is employed. Utilizing Fermat's principle and the fast-progression method, the transport of solute along the path of least resistance is accurately simulated, improving the imaging accuracy of structures such as karst conduits.

[0272] Building upon this, a spatial adaptive regularization mechanism based on the sensitivity field is proposed. By dynamically adjusting the weights according to ray coverage and sensitivity strength, a balance is achieved between strong constraint in the blind zone and high resolution in the hot zone.

[0273] The optional embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solution of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A multi-field coupled inversion analysis method for water-rich characteristics of karst disaster sources, characterized in that, include: Obtain the tracer penetration curve and water head change curve of the target area, extract features from the tracer penetration curve, and generate a travel time feature set; A unified parameter state vector containing permeability coefficient, water storage coefficient and effective porosity is constructed to establish a forward model describing solute transport process and groundwater flow process; Based on the travel time feature set, the head change curve, and the output of the forward model, the tracer observation residual and the head observation residual are calculated. Based on the forward model, the sensitivity matrix of the unified parameter state vector to the travel time feature set and the head change curve is calculated, and a joint objective function including a spatial adaptive regularization term is constructed accordingly. The joint objective function is solved iteratively to obtain the optimized unified parameter state vector, and the location of the water-rich disaster-causing structure is delineated based on the optimized parameter distribution. The joint objective function also includes the weighted sum of squares of the tracer observation residuals and the head observation residuals; The unified parameter state vector consists of the logarithmic values ​​of the permeability coefficient, the logarithmic value of the water storage coefficient, and the original value of the effective porosity of each grid cell within the target area. The joint objective function Φ(θ) is specifically constructed according to the following formula: Φ(θ) = ∑ i=1 N [w _i_T (T _i_obs -T _i_cal (θ))] 2 +∑ j=1 P [w _j_H (H _j_obs -H _j_cal (θ))] 2 +R(θ); where θ is the unified parameter state vector; N is the total number of tracer travel-time data; T _i_obs T represents the i-th tracer travel time observation; _i_cal (θ) is the calculated value of the i-th tracer travel time based on the current parameter θ; w _i_T The weight of the i-th tracer travel time data; P is the total number of head observation data; H _j_obs H is the j-th head observation value; _j_cal (θ) represents the j-th head calculated based on the current parameter θ; w _j_H Let θ be the weight of the j-th head data point; R(θ) is the spatial adaptive regularization term; The weights of the tracer travel time data and the head data are calculated using the sensitivity matrix and its corresponding measurement error, according to the following formula: w _i_T =(||J _i_T || _F ) / (σ T 2 ×∑ k=1 N ||J _k_T || _F ), w _j_H =(||J _j_H || _F ) / (σ H 2 ×∑ k=1 P ||J _k_H || _F ); where, ||…|| _F J represents the Frobenius norm; _i_T σ is the row vector in the sensitivity matrix corresponding to the i-th tracer travel time data; T J represents the standard deviation of the measurement error for tracer travel time data. _j_H σ is the row vector corresponding to the j-th head data point in the sensitivity matrix; H Let J be the standard deviation of the measurement error of the head data, k be the summation index variable, and J be the standard deviation of the measurement error of the head data. _k_T J _k_H These are the row vectors corresponding to the k-th tracer travel time data and head data in the sensitivity matrix, respectively; The spatial adaptive regularization term R(θ) is constructed based on the first-order difference operator matrix and the iteratively updated regularization weight matrix. The specific calculation formula is: R(θ) = θ T ×Λ _k ×L T ×L×θ; where θ is the unified parameter state vector; L is the first-order difference operator matrix used to calculate the spatial gradient of the model parameters; T For transpose operation; Λ _k This is the spatial adaptive regularization weight matrix in the k-th iteration step, used to weight the smoothing constraint strength at different spatial locations; The construction of the regularization weight matrix depends on the comprehensive sensitivity index of each grid cell within the target region; the comprehensive sensitivity index ξ of the m-th grid cell in the k-th iteration step. _m_k Specifically, calculate according to the following formula: ξ _m_k =sqrt(∑ i=1 N [( T _i_cal ) / ( θ _m )×w _i_T ] 2 +∑ j=1 P [( H _j_cal ) / ( θ _m )×w _j_H ] 2 ); where θ _m This represents the parameter component corresponding to the m-th grid cell in the unified parameter state vector; T _i_cal ) / ( θ _m () represents the partial derivative of the calculated value of the i-th travel time with respect to the parameter of the m-th grid cell; H _j_cal ) / ( θ _m ) represents the partial derivative of the j-th calculated head with respect to the m-th grid element parameter; w _i_T and w _j_H These are the corresponding weights.

2. The method according to claim 1, characterized in that, The regularization weight matrix is ​​a diagonal matrix, and its m-th diagonal element λ _mm_k The value is determined based on the normalized value of the comprehensive sensitivity index, and the specific calculation formula is as follows: x _m_k '=(ξ _m_k -x _min_k ) / (ξ _max_k -x _min_k ); l _mm_k =λ _0_k ×[1+γ×(1 / ξ _m_k ') β ]; Where, ξ _m_k ' represents the normalized relative sensitivity; ξ _max_k and ξ _min_k λ represents the maximum and minimum values ​​of the comprehensive sensitivity index of all grid cells in the current iteration step, respectively. _0_k γ is the baseline regularization parameter for the current iteration step; γ is the regularization adjustment coefficient; β is the shape control parameter.

3. The method according to claim 1, characterized in that, The slowness parameter in the forward model depends on the hydraulic gradient field, and the iterative solution process of the joint objective function adopts a prediction-correction alternating update strategy.

4. The method according to claim 3, characterized in that, The prediction-correction alternating update strategy performs the following steps in each iteration: Based on the permeability coefficient field of the current iteration step, solve the steady-state groundwater flow control equation to obtain the current hydraulic head field and hydraulic gradient field; Based on the hydraulic gradient field, the permeability coefficient field of the current iteration step, and the effective porosity field, update the slowness value of each grid cell; Based on the updated slowness value, the tracer travel time and its sensitivity matrix are recalculated, and an optimization solution for the joint objective function is performed to obtain the updated unified parameter state vector.

5. The method according to claim 4, characterized in that, When calculating the tracer travel time, a curved ray tracing method based on Fermat's shortest time principle is used, specifically including: Based on the current slowness distribution, the fast-progression method is used to solve the process function equation, and the shortest arrival time field of each point in the target area is calculated with the tracer source point as the starting point. Starting from the monitoring point, backtrack along the negative gradient direction of the shortest arrival time field to determine the curved ray path connecting the tracer source point and the monitoring point; The tracer travel time is calculated based on the sum of the products of the length of the curved ray path within each grid cell and the slowness value of the corresponding grid cell.

Images