Compaction degree measuring and calculating method based on multi-physics field joint inversion

Abstract

The invention discloses a compaction degree measuring and calculating method based on multi-physics field joint inversion. The method comprises the following steps: acquiring multi-physics field observation data and rock-soil compaction characteristic parameters of a to-be-detected area; establishing a mapping relationship between the moisture content and the observation data and between the dry density and the observation data by using a pre-established forward model; executing joint inversion iteration, and updating the state variable based on the observation residual error; in the iteration process, a physical feasible region is constructed by utilizing rock-soil compaction characteristics, and the physical feasible region is embedded into the inversion process as a constraint condition to constrain the value ranges of the dry density and the water content; and in response to the condition that the joint inversion iteration process meets a preset convergence condition, outputting the dry density after inversion convergence and calculating the compaction degree. According to the method, the problems that a traditional inversion solution violates the soil engineering compaction rule and cannot adapt to compaction under variable work conditions are solved, and the physical authenticity and precision of compaction quality detection are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geotechnical engineering quality detection, and particularly relates to a compaction degree calculation method based on multi-physical field joint inversion. BACKGROUND

[0002] Compaction degree is a core index for evaluating the quality of filling engineering such as embankment, roadbed and dam, and is directly related to the stability and durability of the engineering structure. With the improvement of engineering fine management demand, the multi-physical field joint inversion using ground penetrating radar (GPR) and seismic wave / surface wave technology has become an important technical means to realize wide coverage and non-destructive testing of soil. Through the fusion of complementary information of different physical fields, the multi-physical field joint inversion aims to solve the multi-solution problem of a single geophysical method, so as to more accurately invert the key parameters such as water content and dry density of soil.

[0003] The current joint inversion method is usually based on a data-driven framework, that is, the optimal solution is found by minimizing the residual between the observed data (such as radar electromagnetic wave travel time and surface wave dispersion curve) and the forward simulation data. The mainstream algorithm of the existing technology adopts the least square method or Bayesian inference, and uses an empirical formula (such as CRIM model or Topp formula) to establish the relationship between the physical parameters and the geophysical response. In the inversion process, only a simple upper and lower limit range constraint is usually set for the parameters, or a general mathematical regularization method (such as Tikhonov regularization) is used to smooth the distribution of the solution, mainly focusing on meeting the fitting accuracy of the observed data.

[0004] However, the existing multi-physical field joint inversion method still has technical defects in actual engineering application, mainly in the physical blindness of the inversion process and the static limitation of the constraint mechanism: specifically, the existing inversion algorithm lacks the cognition of the essential physical law of soil. As a special granular material, the dry density and water content of soil strictly follow the compaction characteristic (Proctor curve) law. The traditional method only pursues the minimization of data fitting residual, resulting in that the combination of dry density and water content obtained by inversion often falls outside the compaction curve, for example, the appearance of ultra-high dry density under extremely low water content. This mathematically optimal but physically unfeasible solution deviates from the real compaction state. In addition, the existing constraint model cannot adapt to the variable power compaction environment. In actual construction, the vibration frequency and speed of the road roller change in real time, resulting in that the maximum dry density limit of the soil dynamically fluctuates with the compaction energy. The existing technology still uses fixed and unchanged parameter constraints in the inversion process, ignoring the change of the energy limit of the soil state, resulting in the decrease of the inversion accuracy in over-compaction or under-compaction conditions. SUMMARY

[0005] The present application provides a compaction degree calculation method based on multi-physical field joint inversion, in order to solve the above problems existing in the prior art.

[0006] The technical scheme is a compaction degree calculation method based on multi-physical field joint inversion, comprising:

[0007] Obtaining multi-physical field observation data and geotechnical compaction characteristic parameters of a region to be detected;

[0008] Using a pre-constructed multi-physical field forward model to establish a mapping relationship between state variables and the multi-physical field observation data, the state variables at least including water content and dry density;

[0009] Based on the mapping relationship, performing a joint inversion iteration process, updating the state variables based on the residual between the multi-physical field observation data and the predicted output of the multi-physical field forward model; constructing a physical feasible region using the geotechnical compaction characteristic parameters, and embedding the physical feasible region as a constraint condition into the joint inversion iteration process to constrain the value range of the dry density and the water content;

[0010] In response to the joint inversion iteration process satisfying a preset convergence condition, outputting the dry density after inversion convergence, and using the dry density after inversion convergence and a preset reference maximum dry density to calculate the compaction degree.

[0011] According to an aspect of the present application, the multi-physical field observation data includes ground penetrating radar data and wave speed related data;

[0012] and simultaneously collected road roller real-time working condition data, the road roller real-time working condition data including vibration wheel equivalent normal force amplitude, rolling speed and vibration frequency.

[0013] According to an aspect of the present application, the multi-physical field forward model includes a dielectric constant forward model and a wave speed forward model;

[0014] The dielectric constant forward model is constructed based on a pre-configured complex refractive index model, and is used to define a functional relationship between the dielectric constant and the volume water content derived from the dry density and the water content;

[0015] The wave speed forward model is constructed based on the poroelastic theory or the empirically corrected Hardin formula, and is used to define a functional relationship between the wave speed and the saturation derived from the dry density and the water content, and the preset vertical effective stress.

[0016] According to an aspect of the present application, embedding the physical feasible region as a constraint condition into the joint inversion iteration process comprises:

[0017] Constructing a joint loss function including a data fitting term and a constraint penalty term;

[0018] Calculating the constraint violation degree of the dry density and the water content at the current iteration step with respect to the boundary of the physical feasible region;

[0019] A numerical value of the constraint penalty term is calculated based on the constraint violation degree, and the numerical value of the constraint penalty term is superimposed into the joint loss function, and the state variable is updated by minimizing the joint loss function.

[0020] According to an aspect of the present application, the constraint penalty term is constructed as a piecewise function with respect to the normalized constraint violation degree, and the piecewise function includes a soft constraint region, a transition region and a hard constraint region:

[0021] In the soft constraint region, the constraint penalty term is a quadratic function with respect to the normalized constraint violation degree, which is used to allow the adjustment of the inversion solution within the confidence range of the physical feasible region;

[0022] In the transition region, the constraint penalty term is a cubic polynomial function, which is used to maintain the continuity of the first derivative of the constraint penalty term and realize the smooth transition from soft constraint to hard constraint;

[0023] In the hard constraint region, the constraint penalty term is an exponential function, which is used to apply exponential growth penalty to the inversion solution exceeding the preset threshold.

[0024] According to an aspect of the present application, the physical feasible region is represented as a confidence region with statistical characteristics, and before performing the joint inversion iterative process, the confidence region is also constructed, specifically:

[0025] The pre-acquired indoor compaction test data is polynomially fitted to obtain a center line model of the compaction curve;

[0026] The fitting residual of the indoor compaction test data with respect to the center line model is analyzed to estimate the fitting standard uncertainty;

[0027] The pre-stored site correction coefficient reflecting the site construction variability is introduced, and the fitting standard uncertainty is expanded using the site correction coefficient to obtain a comprehensive standard uncertainty;

[0028] Based on the center line model and the comprehensive standard uncertainty, the upper and lower boundaries of the confidence region are determined.

[0029] According to an aspect of the present application, the physical feasible region is embedded into the joint inversion iterative process as a constraint condition, including:

[0030] A hard projection operator defined on the physical feasible region is constructed, which is used to calculate the minimum distance projection of a point in the state space to the physical feasible region;

[0031] In each update step of the joint inversion iterative process, the updated state variable is corrected using the hard projection operator;

[0032] When the updated state variable falls outside the physical feasible region, it is forced to be mapped to the boundary or inside of the physical feasible region to obtain a corrected state variable satisfying the physical constraint.

[0033] According to one aspect of the present application, the physical feasible region is an energy self-adaptive region with dynamically variable boundaries, and the method further comprises:

[0034] Based on the real-time working condition data of the road roller, a compaction energy index of the current spatial position is calculated;

[0035] The boundary shape and position of the physical feasible region are dynamically adjusted according to the compaction energy index, so that the dry density range covered by the physical feasible region moves with the increase of the compaction energy index.

[0036] According to one aspect of the present application, dynamically adjusting the boundary shape and position of the physical feasible region according to the compaction energy index comprises:

[0037] A logarithmic growth relationship between the maximum dry density and the compaction energy index, and a linear offset relationship between the optimal water content and the compaction energy index are established by using a pre-calibrated energy mapping model;

[0038] Based on the compaction energy index, the maximum dry density and the optimal water content in the rock-soil compaction characteristics parameters are updated in real time, and the center curve and the boundary bandwidth of the physical feasible region are reconstructed accordingly.

[0039] According to one aspect of the present application, the multi-physical field observation data contain time series observation data of the same detection area at different rolling times, and the method further comprises introducing a time series monotonicity constraint in the joint inversion iteration process, specifically:

[0040] The rolling time sequence corresponding to the time series observation data is identified;

[0041] Based on the rolling time sequence, a constraint condition is constructed to limit the dry density estimation value corresponding to the next rolling time to be not lower than the dry density estimation value corresponding to the previous rolling time, as the time series monotonicity constraint;

[0042] The time series monotonicity constraint is added as an additional penalty term to the objective function of the joint inversion iteration process.

[0043] Beneficial effects, the present application solves the problems that the traditional inversion solution violates the soil compaction rule and cannot adapt to variable working condition compaction, and improves the physical authenticity and precision of compaction quality detection. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 A step flowchart of a compaction degree measurement method based on multi-physical field joint inversion provided by an embodiment of the present application.

[0045] Figure 2 A step flowchart of embedding the physical feasible region as a constraint condition into the joint inversion iteration process provided by an embodiment of the present application.

[0046] Figure 3 A step flow chart for constructing a confidence domain provided by the embodiment of the present application is provided.

[0047] Figure 4 A step flow chart for dynamically adjusting the boundary shape and position of the physical feasible domain according to the compaction energy index provided by the embodiment of the present application is provided. DETAILED DESCRIPTION

[0048] In order for those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor should fall within the scope of protection of the present application.

[0049] It should be noted that the terms include and have as well as any variations thereof are intended to cover inclusive rather than exclusive inclusion, for example, a process, method, system, product or apparatus that comprises a list of steps or units is not necessarily limited to those steps or units that are clearly listed, but can include other steps or units that are not clearly listed or inherent to such processes, methods, products or apparatuses.

[0050] As shown in Figure 1 A compaction degree measurement method based on multi-physical field joint inversion, comprising the following steps:

[0051] Obtain multi-physical field observation data and geotechnical compaction characteristic parameters of the region to be detected.

[0052] In the embodiment, obtaining the multi-physical field observation data of the region to be detected specifically means collecting the response data of different physical fields inside the target engineering structure through field detection equipment. The multi-physical field observation data at least includes ground penetrating radar data and wave speed related data. The ground penetrating radar data can reflect the dielectric properties of the soil, and the wave speed related data can reflect the elastic modulus or shear stiffness properties of the soil. Obtaining the geotechnical compaction characteristic parameters specifically means obtaining the physical law parameters describing the specific relationship between the dry density and the water content of the soil. For example, the key parameters of the compaction curve, including the maximum dry density, the optimum water content and the shape coefficient of the curve, can be obtained through indoor standard Proctor compaction test. This provides a basic physical basis for subsequent construction of the physical feasible domain.

[0053] Using the pre-constructed multi-physical field forward model, a mapping relationship between the state variables and the multi-physical field observation data is established, and the state variables at least include the water content and the dry density.

[0054] Specifically, the state variables refer to unknown physical quantities that need to be solved in the inversion process. The water content w and the dry density p d As the core state variables, because these two parameters directly determine the engineering properties of the soil and are closely related to the compaction characteristics. The multi-physical field forward model refers to a mathematical model that can calculate the theoretical observation data according to the given state variables. Exemplarily, the multi-physical field forward model includes a dielectric constant forward model and a wave velocity forward model. The dielectric constant forward model is used to calculate the theoretical dielectric constant under the given water content and dry density, and then simulate the ground penetrating radar response; the wave velocity forward model is used to calculate the theoretical wave velocity under the given water content and dry density. By establishing a mapping relationship, the forward prediction from engineering parameters to geophysical responses is realized.

[0055] Based on the mapping relationship, a joint inversion iterative process is performed to update the state variables based on the residual between the multi-physical field observation data and the predicted output of the multi-physical field forward model.

[0056] In this embodiment, the joint inversion iterative process is a numerical optimization process that continuously adjusts the state variables to minimize the difference between the observation data and the predicted data. Specifically, based on the initial state variables, the multi-physical field forward model is used to calculate the predicted output; the residual between the predicted output and the actually obtained multi-physical field observation data is calculated; and the state variables are updated according to the size and gradient direction of the residual using an optimization algorithm. Common optimization algorithms include Gauss-Newton method, conjugate gradient method or L-BFGS-B algorithm, etc. Through repeated iteration, the predicted data gradually approaches the observation data, and the true state variables are inferred.

[0057] The geotechnical compaction characteristic parameters are used to construct a physically feasible region, and the physically feasible region is embedded as a constraint condition in the joint inversion iterative process to constrain the value range of the dry density and the water content.

[0058] Specifically, the physically feasible region refers to a reasonable value region in the water content-dry density state space that meets the geotechnical compaction physical law. Unlike traditional inversion which only relies on data fitting, this embodiment requires that the state variables in the inversion process must fall within or near this physically feasible region. Alternatively, the way to construct the physically feasible region can be based on a static compaction curve or a dynamic compaction curve that varies with compaction energy. The way to embed the physically feasible region as a constraint condition can be to add a penalty function term in the objective function to achieve soft constraint, or to use a projection operator to achieve hard constraint. The solution that can fit the observation data but violates the soil physics common sense is effectively eliminated, and the physical authenticity and reliability of the inversion result are improved.

[0059] In response to the joint inversion iterative process satisfying a preset convergence condition, the dry density after inversion convergence is output, and the compaction degree is calculated using the dry density after inversion convergence and a preset reference maximum dry density.

[0060] In the embodiment, the convergence condition includes that the relative change amount of the objective function value is less than a preset threshold value, or the update amount of the state variable is less than a preset threshold value. When the convergence condition is satisfied, the current state variable is considered as the optimal estimation value. Outputting the dry density after inversion convergence means that the dry density obtained through the final iteration is output as the detection result. Further, in order to evaluate the compaction quality, it is necessary to calculate the compaction degree. The reference maximum dry density is the denominator for calculating the compaction degree, and is usually taken from the fixed maximum dry density determined by the standard compaction test. The specific formula for calculating the compaction degree is the dry density after inversion convergence divided by the reference maximum dry density. It should be noted that whether the dynamically adjusted maximum dry density is used as a constraint in the inversion process, the reference maximum dry density should be uniformly used when finally evaluating the compaction degree, so that the evaluation standard conforms to the engineering specification and has comparability.

[0061] In the embodiment, the geotechnical compaction characteristics are parameterized as a physically feasible region and embedded in the inversion iteration, so that the geophysical data and the soil test rules are deeply integrated.

[0062] In one possible embodiment, the multi-physical field observation data includes ground penetrating radar data and wave speed related data; and real-time working condition data of the road roller is synchronously collected, and the real-time working condition data of the road roller includes an equivalent normal force amplitude of the vibration wheel, a rolling speed and a vibration frequency.

[0063] In the embodiment, the composition of the multi-source data is defined. The ground penetrating radar data specifically refers to the signals obtained by transmitting and receiving high-frequency electromagnetic waves by a ground penetrating radar antenna, such as radar profile, reflection wave travel time or amplitude information, etc. The wave speed related data specifically refers to physical quantities related to wave speed obtained by seismic wave or surface wave detection means, such as the phase velocity dispersion curve of Rayleigh wave. The real-time working condition data of the road roller is collected synchronously to provide a basis for subsequent analysis of compaction energy. The equivalent normal force amplitude F v (t) reflects the vertical force applied by the road roller on the soil; the rolling speed v(t) reflects the duration of energy action; and the vibration frequency f(t) reflects the loading frequency of energy. In addition, the number of rolling passes N, the width of the vibration wheel B and the thickness of the layer H can also be collected. Further, the multi-physical field observation data needs to be processed through spatial alignment and time synchronization, and mapped to the same detection grid.

[0064] In one specific embodiment, the collection parameters of the ground penetrating radar data include: a center frequency range of 400 MHz to 2 GHz, selected according to the detection depth requirement; a sampling point number of no less than 512 points / channel; and a scanning interval of 0.02 m to 0.10 m. The collection parameters of the wave velocity related data include: a frequency range of 5 Hz to 100 Hz for excitation during surface wave detection; a detector interval of 0.5 m to 2.0 m; and a sampling frequency of no less than 1000 Hz. The collection frequency of the roller compactor working condition data is no less than 10 Hz, so as to ensure that the transient response of the vibrating wheel can be captured.

[0065] In one exemplary embodiment, the multi-physical field forward modeling includes a dielectric constant forward model and a wave velocity forward model; the dielectric constant forward model is constructed based on a pre-configured complex refractive index model, and is used to define a functional relationship between the dielectric constant and the bulk water content derived from the dry density and the water content.

[0066] Specifically, the complex refractive index model (CRIM model) is used to establish the relationship between the dielectric constant and the three-phase (solid, liquid, and gas) components of the soil. The formula is as follows:

[0067] sqrt(ε) = (1 - n) * sqrt(ε s ) + θ * sqrt(ε w ) + (n - θ) * sqrt(ε a );

[0068] wherein ε represents the equivalent relative dielectric constant of the soil mixture; n represents the porosity, which can be calculated by the formula n = 1 - (ρ d / ρ s ), ρ d being the dry density and ρ s being the soil particle density; ε s represents the dielectric constant of the soil particle; θ represents the bulk water content, which can be calculated by the formula θ = w * (ρ d / ρ w ), w being the mass water content and ρ w being the density of water; ε w represents the dielectric constant of the pore water; and ε a represents the dielectric constant of air. The forward mapping of the state variables (w, ρ d ) to the observation variables (ε) is established.

[0069] The wave velocity forward model is constructed based on the poroelastic theory or the Hardin formula with empirical correction, and is used to define a functional relationship between the wave velocity and the dry density, the saturation derived from the water content, and the preset vertical effective stress.

[0070] Specifically, the Hardin formula is a commonly used empirical formula in geotechnical engineering to describe the small-strain shear modulus or shear wave velocity of soil. The Hardin formula modified by saturation is used to establish the relationship between shear wave velocity and state variables. The specific formula is as follows:

[0071] V s =A*(ρ d / ρ ref ) m *(σ' v / p a ) k *F(S r );

[0072] where V s represents the predicted shear wave velocity; A is an empirical coefficient related to soil quality; ρ d is the dry density; ρ ref is the reference dry density; m is the density index; σ' v is the vertical effective stress; p a is the atmospheric pressure; k is the stress index; S r is the saturation, which can be calculated by the formula S r =θ / n; F(S r ) is the saturation correction function, for example, F(S r )=1-λ*S r , λ is the weakening coefficient of water to the stiffness of soil skeleton. Quantitatively describes the physical mechanism that the increase of dry density increases the wave velocity, and the increase of water content (saturation) reduces the wave velocity. Alternatively, the vertical effective stress is calculated based on the preset thickness of the filling layer and the bulk density of the soil, or obtained by field measurement.

[0073] As Figure 2 shown, according to one aspect of the present application, the physically feasible region is embedded as a constraint condition in the joint inversion iterative process, including:

[0074] Construct a joint loss function including data fitting term and constraint penalty term.

[0075] In this embodiment, the total objective function of inversion is constructed. The joint loss function L total is composed of a data fitting term L data and a constraint penalty term, and can also include a prior information regularization term. The specific expression is as follows:

[0076] L total =L data +λ P *P(δ)+λ w *(w-w prior ) 2 / σ w 2;

[0077] wherein L data represents the fitting residual of multi-physical field observation data, usually in the form of weighted least squares, for example, L data =||y gpr -h gpr (x)|| 2 / σ g 2 +||y seis -h seis (x)|| 2 / σ s 2 ; y gpr is the measured observation data of ground penetrating radar (GPR); h gpr (x) is the forward model prediction value of ground penetrating radar; σ g is the standard deviation of ground penetrating radar observation data; y seis is the measured observation data of seismic wave (seis); h seis (x) is the forward model prediction value of seismic wave; σ s is the standard deviation of seismic wave observation data; λ P represents the weight coefficient of constraint penalty term; P(δ) represents the penalty function value calculated based on the constraint violation degree δ; λ w represents the weight of water content prior term; w prior represents the water content prior value; σ w represents the standard deviation of water content prior; w is the water content of rock-soil. By minimizing the joint loss function, the data fitting degree and the physical constraint condition can be considered at the same time.

[0078] As Figure 3 shown, in one possible implementation, the physically feasible region is represented as a confidence region with statistical characteristics, and before performing the joint inversion iteration process, the confidence region is also constructed, specifically:

[0079] Polynomial fitting is performed on the pre-acquired indoor compaction test data to obtain the center line model of the compaction curve;

[0080] The fitting residual of the indoor compaction test data relative to the center line model is analyzed to estimate the fitting standard uncertainty;

[0081] The pre-stored field correction coefficient reflecting the variability of field construction is introduced, and the fitting standard uncertainty is expanded by using the field correction coefficient to obtain the comprehensive standard uncertainty;

[0082] Based on the center line model and the comprehensive standard uncertainty, the upper and lower boundaries of the confidence region are established.

[0083] Exemplarily, based on n groups of indoor compaction test data (wi , p d，i ) is fitted with a quadratic polynomial center line by least square method:

[0084] p d，fit (w) = a * w 2 + b * w + c;

[0085] where p d，fit (w) is the fitted value of dry density corresponding to water content w; a, b, c are the fitting coefficients. The fitting residual e i = p d，i - p d，fit (w i ) is calculated, p d，i is the measured value of dry density of the geotechnical body measured by the i-th indoor compaction test; and the fitting standard uncertainty s fit is estimated:

[0086] s fit = sqrt(∑(e i 2 ) / (n-3));

[0087] where ∑ represents summation, n-3 is the degree of freedom, and n is the total number of indoor compaction tests. A field correction coefficient K is introduced to consider the difference between indoor tests and field environment, and the comprehensive standard uncertainty s total is calculated:

[0088] s total = sqrt(s fit 2 + (K * p dmax ) 2 );

[0089] where the typical value range of K is 0.02 to 0.05; p dmax is the maximum dry density. The upper and lower boundaries of the confidence interval are determined:

[0090] p d，upper (w) = p d，fit (w) + z α * s total ; p d，lower (w) = p d，fit (w) - z α * s total ;

[0091] where p d，upper (w) is the upper boundary value of the dry density confidence interval corresponding to water content w; p d，lower (w) is the lower boundary value of the dry density confidence interval corresponding to water content w; z α is the confidence level coefficient, for example, 1.96 corresponds to 95% confidence.

[0092] The dry density and water content of the current iteration step are calculated, and the constraint violation degree of the current iteration step is calculated with respect to the boundary of the physically feasible region.

[0093] Specifically, the degree of deviation of the current state variable from the physically feasible region is quantified. For the dry density ρ d and water content w of the current iteration, the normalized constraint violation degree δ is calculated with respect to the upper boundary ρ d，upper (w) of the confidence region. The specific calculation logic is as follows: if ρ d ≤ ρ d，upper (w), then δ = 0, indicating that the constraint is satisfied; if ρ d > ρ d，upper (w), then δ = (ρ d - ρ d，upper (w)) / σ total . Similarly, the violation degree with respect to the lower boundary can also be calculated. Here, the case of excessive dry density (exceeding the upper boundary) is taken as an example, because the virtual high dry density is a common non-physical phenomenon in inversion.

[0094] In a preferred embodiment, the constraint penalty term is constructed as a piecewise function with respect to the normalized constraint violation degree, including a soft constraint region, a transition region, and a hard constraint region:

[0095] In the soft constraint region, the constraint penalty term is a quadratic function with respect to the normalized constraint violation degree, which is used to allow adjustment of the inversion solution within the confidence range of the physically feasible region;

[0096] In other words, in the soft constraint region, the constraint penalty term is a quadratic function with respect to the normalized constraint violation degree, which is used to allow adjustment of the inversion solution within the confidence range of the physically feasible region within a preset tolerance range;

[0097] In the transition region, the constraint penalty term is a cubic polynomial function, which is used to maintain the continuity of the first derivative of the constraint penalty term and achieve a smooth transition from soft constraint to hard constraint;

[0098] In the hard constraint region, the constraint penalty term is an exponential function, which is used to impose an exponentially growing penalty on the inversion solution that exceeds the preset threshold.

[0099] In this embodiment, a specific penalty function form P(δ) is defined. In the soft constraint region (0 < δ ≤ δ soft ), a quadratic function is used:

[0100] P(δ) = 0.5*β*δ 2 ;

[0101] where β is the soft constraint coefficient, and δ soft is the soft constraint threshold, for example, 1.0. In the hard constraint region (δ > δ hard ), an exponential function is used:

[0102] P(δ) = M*exp(η*(δ-δ hard ))+C hard ;

[0103] where M, η are hard constraint parameters, C hard is a constant term of the hard constraint region exponential penalty function, δ hard is a hard constraint threshold, for example, take 2.5. In the transition region (δ soft <δ≤δ hard ), a cubic polynomial is adopted:

[0104] P(δ) = a3*δ 3 +a2*δ 2 +a1*δ+a0;

[0105] In order to ensure the continuity of the function and its first derivative at δ soft and δ hard , the coefficients a3, a2, a1, a0 need to be determined according to the boundary conditions. In particular, the cubic term coefficient γ (i.e. a3) needs to meet the smoothness condition. Through the three-section design, when the violation is small (within the confidence domain), the penalty is small and smooth, allowing data-driven fine-tuning; when the violation is too large (outside the confidence domain), the penalty increases sharply, forcing the solution to be pulled back to the physically reasonable interval.

[0106] The numerical value of the constraint penalty term is calculated based on the constraint violation degree, and the numerical value of the constraint penalty term is superimposed into the joint loss function, and the state variable is updated by minimizing the joint loss function.

[0107] Specifically, the calculated P(δ) is substituted into L total . The L-BFGS-B optimization algorithm is used to calculate the gradient of L total about the state variables w and p d , and the state variables are updated along the negative gradient direction. Since P(δ) is designed as a smooth function, the continuity of the gradient is guaranteed, which is conducive to the stable convergence of the optimization algorithm. The final inversion solution can not only fit the observation data well, but also meet the physical constraints of the soil compaction characteristics.

[0108] According to another aspect of the present application, embedding the physically feasible region as a constraint condition into the joint inversion iteration process can also be:

[0109] A hard projection operator defined on the physically feasible region is constructed, and the hard projection operator is used to calculate the minimum distance projection of a point in the state space to the physically feasible region.

[0110] In the present embodiment, the core tool for constraint execution is defined. The hard projection operator Π Ωis a geometric operator that can map any given state variable point (moisture content w, dry density p d ) to the nearest point within the physical feasible region Ω. It is equivalent to solving a constrained quadratic programming problem. Specifically, for the candidate state variable x update updated at the current iteration step, find the modified state variable x corrected such that the Euclidean distance or weighted Mahalanobis distance between them is minimized, and x corrected must lie inside or on the boundary of the physical feasible region Ω. The introduction of the hard projection operator enables the use of efficient algorithms such as the projected gradient descent method for the inversion process, ensuring convergence speed while strictly satisfying physical constraints.

[0111] In further embodiments, the physical feasible region is an energy-adaptive region with dynamically variable boundaries, and the method further comprises:

[0112] calculating a compaction energy index at the current spatial position based on real-time working condition data of the road roller.

[0113] Specifically, the compaction energy index E f is a physical quantity that reflects the mechanical work received by unit volume of soil. Based on the working condition data, the compaction energy index can be calculated by an integral formula. The specific calculation formula is as follows:

[0114] E f = (1 / (B*H))*integral(F v (t)*v(t)*dt);

[0115] where E f represents the equivalent compaction energy per unit volume; B represents the width of the vibrating wheel; H represents the current fill layer thickness; F v (t) represents the equivalent normal force amplitude applied by the vibrating wheel to the soil at time t, which is usually determined by the static line load and excitation force of the road roller; v(t) represents the rolling speed at time t; dt represents the time differential; integral represents the integral operation within the effective rolling time window. The discrete working condition parameters such as vibration frequency, speed, and excitation force of the road roller are integrated into a unified energy index, laying a foundation for subsequent establishment of the energy-density physical correlation.

[0116] According to the compaction energy index, the boundary shape and position of the physical feasible region are dynamically adjusted, so that the dry density range covered by the physical feasible region moves with the increase of the compaction energy index.

[0117] As Figure 4 shown, in a preferred implementation, dynamically adjusting the boundary shape and position of the physical feasible region according to the compaction energy index includes:

[0118] A logarithmic growth relationship between the maximum dry density and the compaction energy index, and a linear offset relationship between the optimum water content and the compaction energy index are established by using the pre-calibrated energy mapping model.

[0119] Based on the compaction energy index, the maximum dry density and the optimum water content in the soil compaction characteristics are updated in real time, and the center curve and the boundary bandwidth of the physical feasible region are reconstructed based on the updated maximum dry density and optimum water content.

[0120] In this embodiment, a quantitative mapping model from energy to compaction characteristics is established. According to the principle of soil mechanics, with the increase of compaction work (energy), the compaction curve of soil body will change in shape: the maximum dry density increases, and the corresponding optimum water content decreases (moves to the upper left). A specific empirical formula can be used to describe this rule. For the maximum dry density, a logarithmic growth model is used:

[0121] ρ dmax (E f ) = ρ dmax_std + k ρ *log(1 + E f / E ref );

[0122] where ρ dmax (E f ) represents the maximum dry density under the current energy; ρ dmax_std represents the maximum dry density under the standard compaction energy; k ρ is the density growth coefficient, reflecting the sensitivity of the soil to energy; log is the natural logarithm function; E ref is the reference energy value. For the optimum water content, a linear offset model is used:

[0123] w opt (E f ) = w opt_std -k w *(E f -E ref );

[0124] where w opt (E f ) represents the optimum water content under the current energy; w opt_std represents the optimum water content under the standard compaction energy; k w is the water content offset coefficient. Based on the real-time calculation of ρ dmax (E f ) and w opt (E f ), the current parabolic compaction center line is regenerated, and based on this center line, combined with the preset bandwidth Δ, the current dynamic physical feasible region Ω(E f). For example, when it is detected that the number of compaction passes of a certain area is increased, which leads to E f When the energy E

[0125] In each update step of the joint inversion iteration process, the updated state variable is corrected by using a hard projection operator; when the updated state variable falls outside the physically feasible region, the updated state variable is forced to be mapped to the boundary or inside of the physically feasible region, to obtain a corrected state variable that satisfies the physical constraint.

[0126] Specifically, the specific execution actions of the constraint are described. In the k-th step of the inversion iteration, an unconstrained intermediate solution x k+0.5 is obtained based on the gradient descent direction of the data residual. k+0.5 It is judged whether x f satisfies the current energy E k+1 corresponding to the feasible region condition g(x)≤0; wherein g(x) is a physical constraint function of the state variable x. If it is satisfied, x k+0.5 is directly set to x x If it is not satisfied, for example, the calculated dry density is much higher than the theoretical limit that can be reached by the current energy, a hard projection operator is called to solve the following optimization problem:

[0127] min k+0.5 ||x-x 2 s.t x∈Ω(E f );

[0128] Where x is the state variable vector to be solved. By solving the optimization problem, the point on the boundary of the feasible region closest to the intermediate solution is found as the new iteration solution x k+1 The prediction-correction mechanism ensures that the inversion trajectory is always limited within the physically reasonable channel, effectively suppressing the multi-solution.

[0129] Preferably, the method further comprises performing abnormality detection based on the hard projection operator, specifically:

[0130] When the updated state variable is corrected by using the hard projection operator, the Euclidean distance of the state variable before and after correction in the state space is calculated, denoted as the projection distance; the number of consecutive times that the projection distance exceeds a preset safety threshold in the joint inversion iteration process is counted; in response to the number of consecutive times exceeding a preset count threshold, a material abnormality or model mismatch alarm signal for the current detection area is generated.

[0131] In this embodiment, an intelligent diagnosis mechanism based on inversion process data is provided. The projection distance d projis an index to measure the degree of conflict between observation data and physical laws. The specific calculation formula is:

[0132] d proj = sqrt((w corrected -w update ) 2 +(p dcorrected -p dupdate ) 2 );

[0133] where w corrected and p dcorrected are the values after projection correction; w update and p dupdate are the values before correction. If d proj continues to be large, it indicates that the observation data is strongly inclined to produce a solution that violates the physical law, for example, showing extremely high dielectric constant and wave speed at low energy. It usually means that the preset physical model is invalid, such as sudden change of soil quality, high stone content, or sensor failure. By setting a continuous count threshold K, for example, K = 5, when the projection distance of 5 consecutive iterations all exceeds the safety threshold d th , the system automatically triggers an alarm to prompt the engineering personnel to review the area, rather than forcibly output the results that are corrected by the projection algorithm but may be distorted.

[0134] Further, in actual application, the following boundary conditions and abnormal situation processing also need to be considered: if the ground penetrating radar data or wave speed data is missing at a certain measurement point, spatial interpolation method can be used for completion, or the measurement point is skipped and only the inversion results of valid measurement points are output. If the observation data is obviously beyond the physical reasonable range, for example, the dielectric constant is less than 1 or the wave speed is negative, it should be excluded in the preprocessing stage and marked in the detection report. If the iteration reaches the preset maximum number, for example, 100 times, the current optimal estimation value should be output and a label with reduced confidence is added. For the material abnormal area identified by the projection distance anomaly detection, it is suggested to combine traditional methods such as core sampling for review and verification.

[0135] The embodiment takes advantage of the real-time feedback of energy data of modern intelligent road rollers, constructs a dynamically changing physical feasible region with compaction energy, and forces the inversion solution to meet this dynamic physical constraint through a hard projection operator, which is suitable for variable power conditions and large thickness of filling body compaction quality detection scenarios.

[0136] According to one aspect of the present application, the multi-physical field observation data includes time series observation data of the same detection area at different rolling times, and the method further comprises introducing a time series monotonicity constraint in the joint inversion iteration process, specifically:

[0137] identifying the rolling time sequence corresponding to the time series observation data;

[0138] Based on the compaction pass order, a constraint condition is constructed to limit the dry density estimation value corresponding to the next compaction pass not being lower than the dry density estimation value corresponding to the previous compaction pass, as a time sequence monotonicity constraint;

[0139] The time sequence monotonicity constraint is superimposed into the objective function of the joint inversion iterative process as an additional penalty term.

[0140] In the embodiment, a physically irreversible constraint is introduced for the multi-pass compaction scene. The compaction of the soil body is a densification process. Under normal construction, with the increase of the compaction pass t, the dry density p d (t) should be monotonically non-decreasing. In order to deal with measurement errors, a quadratic penalty function with a dead zone is preferably used to construct the monotonicity constraint term L mono . The specific formula is as follows:

[0141] L mono =∑(Q(Φ t ));

[0142] Wherein ∑ represents the summation of all adjacent passes t and t+1; p t is the quantity of violating monotonicity, defined as p d (t)-p d (t+1). The penalty function Q(p) is defined as: if p t ≤p tol , then Q(p t )=0; if p t >p tol , then Q(p t )=μ*(p t -p tol ) 2 ; p tol is the tolerance threshold (dead zone width) for tolerating small measurement fluctuations; and μ is the penalty coefficient. By minimizing L mono , the algorithm will punish those non-physical inversion paths in which the density of the next pass is significantly smaller than that of the previous pass, forcing the inversion result to maintain a reasonable growth trend in time sequence. In the solving strategy, the block coordinate descent method can be used to optimize the parameters of each pass in time sequence.

[0143] In an embodiment of the present application, the multi-physical field observation data covers a predetermined number of spatial distribution measurement points, and the method further comprises introducing a spatial continuity constraint in the joint inversion iterative process, specifically:

[0144] A measurement point adjacency graph is constructed based on the spatial coordinates of each measurement point, and a graph Laplacian matrix of the measurement point adjacency graph is calculated;

[0145] A graph regularization term is constructed using the graph Laplacian matrix, and the graph regularization term is used to constrain the difference between the dry density estimation values of adjacent measurement points;

[0146] The graph regularization term is added to the objective function of the joint inversion iterative process as an additional penalty term.

[0147] In this embodiment, a local smoothness constraint is introduced for the spatially distributed inversion results. In normal construction, the compaction characteristics of the soil gradually change with the spatial position and should not appear jagged fluctuations. The difference between adjacent nodes can be measured by constructing a graph regularization term R graph . Specifically, an adjacency graph G = (V, E) is constructed for all measurement points i, where V is the set of measurement points and E is the edge set. The adjacency weight W ij between measurement points i and j is defined, which is negatively correlated with the distance d ij between the measurement points, for example, calculated using a Gaussian kernel function. The graph Laplacian matrix L G = D - W is calculated, where D is the degree matrix. The graph regularization term is defined as:

[0148] R graph = p d T * L G * p d ;

[0149] where p d is the dry density vector of all measurement points; L G is the graph Laplacian matrix; and T T is the transpose. R graph is added to the objective function as an additional term, which serves to penalize large fluctuations in dry density between adjacent measurement points and improve the spatial continuity of the inversion results. To avoid over-smoothing, the weight coefficient λ graph of the graph regularization term can be decayed with the number of iterations, for example, λ graph_k = λ graph_0 *(0.95) k ; where λ graph_k is the weight coefficient of the graph regularization term for the kth iteration.

[0150] In some embodiments, the method further comprises introducing an identifiability constraint in the joint inversion iterative process to suppress the equivalence of the solutions of dry density and water content, specifically: based on the multi-physical field forward modeling, calculating the sensitivity matrix of the multi-physical field observation data to the dry density and water content; constructing an identifiability regularization term based on the coupling term of the sensitivity matrix; and adding the identifiability regularization term to the objective function of the joint inversion iterative process to penalize the inversion path that leads to similar physical field responses but large differences in state variables.

[0151] Specifically, for the water content w and dry density p dIn the joint inversion, sometimes changes in the two parameters can lead to opposite changes in the forward response (such as the dielectric constant), causing them to cancel each other out and resulting in uncertainty, i.e., an equivalent solution. An identifiability constraint term R can be constructed. id For example, the observation operator h is calculated with respect to the state variable x=[w, ρ]. d The sensitivity matrix (Jacobi matrix) J h Based on J h The column vector correlation is used to construct the penalty term:

[0152] R id =(J h T *J h ) wd ;

[0153] Among them (J) h T *J h ) wd Represent w and ρ in the Fisher information matrix d The cross terms (non-diagonal elements) are considered. A large cross term indicates a high correlation in the sensitivity directions of the two parameters, which can easily lead to confusion. This can be addressed by minimizing R... id It tends to seek inversion paths with better parameter decoupling to improve the uniqueness and reliability of the solution. It solves the parameter coupling problem in multiphysics inversion.

[0154] In one possible implementation, the compaction degree is calculated using the inverted converged dry density and a preset benchmark maximum dry density, including:

[0155] The maximum dry density determined by the pre-stored standard compaction test is used as the reference maximum dry density, which is a fixed value that does not change with compaction energy.

[0156] Calculate the ratio of the dry density after inversion convergence to the baseline maximum dry density, and output the ratio as the compaction degree, regardless of whether the physical feasible region is dynamically adjusted with compaction energy.

[0157] In this embodiment, although in order to obtain accurate dry density inversion values, the compressibility energy E may be used. f Dynamically changing ρ dmax (E f While serving as a process constraint, the compaction degree K, the final evaluation index, must be reverted to a unified benchmark. The benchmark's maximum dry density ρ dmax_std This refers to a fixed value determined in a standard compaction test or a heavy compaction test, which does not change with the instantaneous operating conditions of the roller on site. The calculation formula is:

[0158] K=(ρ d_inversion / ρ dmax_std )*100%;

[0159] Where ρ d_inversion The dry density estimate obtained from the inversion convergence; ρ dmax_std This refers to the standard maximum dry density determined in the laboratory. It avoids the problem of fluctuating evaluation benchmarks due to variations in roller operating conditions. For example, roller deceleration increases energy and theoretical maximum density, which in turn lowers the calculated compaction degree, making it unreasonable for engineering evaluation. The final output results are diverse, including a color-coded compaction degree map of the entire work area, a list of spatial coordinates of non-compliant areas, and a quality rating report for a specific station.

[0160] This embodiment solves the coordination problem between dynamic constraints in the inversion process and static standards for engineering evaluation, so that the output compaction index meets engineering specifications.

[0161] In another embodiment of this application, a detailed calculation process for the transition region coefficient in the three-segment penalty function is provided, specifically as follows:

[0162] Determine the segmented thresholds and basic parameters of the penalty function.

[0163] In this embodiment, a segmented threshold is set for the normalized constraint violation degree δ. Assume the soft constraint threshold δ... soft =1.0, corresponding to 1 standard deviation, hard constraint threshold δ hard =2.5, corresponding to 2.5 times the standard deviation. Set the quadratic term coefficient β=10.0 for the soft constraint region. Set the exponential term parameters M=100.0 and η=2.0 for the hard constraint region.

[0164] Construct the cubic polynomial function form of the transition region.

[0165] For example, the functional form of the transition region (1.0 < δ ≤ 2.5) is set as follows:

[0166] P trans (δ)=γ*(δ-δ soft ) 3 +k2*(δ-δ soft ) 2 +k1*(δ-δ soft )+k0;

[0167] To simplify the calculation, the origin of the coordinate system is shifted to δ. soft Here, we need to solve for the coefficients γ, k2, k1, and k0.

[0168] Solving for polynomial coefficients based on continuity conditions.

[0169] Specifically, the solution is obtained using boundary conditions of continuous function values ​​and continuous first derivative: Condition 1 (continuous function values ​​on the left boundary): P trans (δsoft )=P soft (δ soft )P soft (1.0) = 0.5 * 10.0 * 1.0 2 =5.0. Therefore, k0=5.0. Condition 2 (continuity of the left boundary derivative): P' trans (δ soft )=P' soft (δ soft )P' soft (δ) = β*δ. The derivative is 10.0 when δ = 1.0. P' trans (δ) Differentiate with respect to δ at δ soft The value at point P is k1. Therefore, k1 = 10.0. Condition 3 (continuity of the function value on the right boundary): P trans (δ hard )=P hard (δ hard )P hard (2.5) = M * exp(η * (2.5 - 2.5)) = M * 1 = 100.0. Substituting δ = 2.5, i.e., the increment Δδ = 1.5, we get equation A: γ * 1.5 3 +k2*1.5 2 +10.0*1.5+5.0=100.0; Condition 4 (continuity of right boundary derivative): P' trans (δ hard )=P' hard (δ hard )P' hard (δ)=M*η*exp(η*(δ-δ hard When δ=2.5, the derivative is 100.0*2.0*1=200.0. Substituting δ=2.5, we get equation B: 3*γ*1.5 2 +2*k²*1.5+10.0=200.0. Where P trans (δ soft P is the transition region penalty function. trans (δ) at the left boundary δ soft The function value at point P; soft (δ soft P is the penalty function for the soft constraint region. soft (δ) at δ soft The function value at point P' trans (δ soft P is the transition region penalty function. trans (δ) at δ soft The first derivative value at point P' soft (δ soft P is the penalty function for the soft constraint region. soft (δ) at δ softthe first derivative value of the penalty function P trans at the right boundary δ hard ;P trans the function value of the penalty function P hard at the right boundary δ hard ;P hard the function value of the penalty function P hard at the right boundary δ hard ;P trans the function value of the penalty function P hard at the right boundary δ trans ;P hard the first derivative value of the penalty function P hard at the right boundary δ hard ;P hard the first derivative value of the penalty function P hard at the right boundary δ

[0170] Solve simultaneously and determine the final function form.

[0171] In this embodiment, equations A and B are solved simultaneously to obtain γ and k2. From equation B, 6.75*γ + 3*k2 = 190.0; from equation A, 3.375*γ + 2.25*k2 = 80.0. Solving this linear equation set, the specific values of γ and k2 can be obtained. For example, if the calculation results in γ = 25.0, then this cubic term coefficient is the key parameter for achieving smooth transition. In the actual inversion code, this calculation process is automatically performed in the initialization stage, so that the gradient calculation in the inversion process does not jump, ensuring the stability of the L-BFGS-B algorithm.

[0172] This embodiment solves the mathematical implementation problem of ensuring that the penalty function is first-order derivative continuous (i.e., smooth transition) between the soft constraint region and the transition region, and between the transition region and the hard constraint region.

[0173] In an exemplary embodiment, assume that the initial conditions of a certain roadbed detection area are as follows: the equivalent dielectric constant ε obs = 12.5 measured by ground penetrating radar; the shear wave velocity V s，obs = 185 m / s obtained by surface wave detection; the reference maximum dry density ρ dmax_std = 1.85 g / cm 3 determined by indoor standard compaction test; the optimal water content w opt = 14.2%; the soil particle density ρ s = 2.68 g / cm 3 ; set the initial state variables: dry density ρ d 0 = 1.70 g / cm 3 , water content w 0= 15.0%. In the first iteration, the predicted dielectric constant is calculated using the CRIM model with p d 0 = 1.70 g / cm 3 and w 0 = 15.0%, the porosity n = 1 - 1.70 / 2.68 = 0.366 and the volumetric water content θ = 0.15 x 1.70 / 1.0 = 0.255 are calculated. Substituting into the CRIM model, the predicted dielectric constant ε = 11.8 is calculated. The predicted wave velocity is calculated using the Hardin formula, assuming a vertical effective stress σ' v = 30 kPa, the predicted shear wave velocity V s = 178 m / s is calculated. The residuals are calculated: dielectric constant residual = 12.5 - 11.8 = 0.7; wave velocity residual = 185 - 178 = 7 m / s. The state variable is updated based on gradient descent (step size a = 0.1), resulting in the unconstrained updated value p d 0.5 = 1.73 g / cm 3 , w 0.5 = 14.6%. According to the center line of the compaction curve p d，fit (14.6%) = 1.82 g / cm 3 , the upper bound p d，upper = 1.82 + 1.96 x 0.03 = 1.88 g / cm 3 . Since p d 0.5 = 1.73 < 1.88, the constraint is satisfied and no correction is needed. After the first iteration: p d 1 = 1.73 g / cm 3 , w 1 = 14.6%. After several iterations: the inversion result p d *= 1.78 g / cm 3 , w*= 14.3% is converged. The compaction degree K = 1.78 / 1.85 x 100% = 96.2% is calculated. The results show that the embodiment can accurately estimate the compaction degree while ensuring that the inversion solution falls within the physically feasible region of the compaction curve.

[0174] Overall, the compaction degree calculation method based on multi-physical field joint inversion includes: obtaining the multi-physical field observation data and the geotechnical compaction characteristic parameters of the region to be detected; a pre-constructed forward model is used to establish the mapping relationship between the water content and dry density and the observation data; joint inversion iteration is performed, and the state variable is updated based on the observation residual; in the iteration process, the geotechnical compaction characteristics are used to construct a physical feasible region, which is embedded as a constraint condition in the inversion process, and the state variable is forced to fall within the interval conforming to the physical law by using a penalty function or a projection operator; finally, the dry density after inversion convergence is output, and the compaction degree is calculated.

[0175] The application introduces a geotechnical compaction characteristic constraint mechanism, parameterizes the compaction curve measured by indoor test, constructs a physical feasible region of the state variable, and uses a penalty function or a geometric hard projection operator of soft and hard constraint switching to forcibly limit the inversion iteration process within the effective interval conforming to the soil rule. Pseudo-solutions that fit the observation data but violate physical common sense are effectively eliminated, and the physical authenticity of the detection result is improved. A dynamic inversion strategy based on energy self-adaptation is proposed, the compaction energy index of the road roller is calculated in real time, the boundary shape of the physical feasible region is dynamically adjusted, such as the maximum dry density increasing with the logarithm of energy. The deviation of the traditional static constraint under the over-pressing or under-pressing working condition is overcome, the intelligent adaptation of the detection algorithm to the complex construction process on site is realized, and the inversion accuracy under the variable energy environment is ensured.

[0176] The above describes the preferred embodiments of the application in detail, but the application is not limited to the specific details in the above embodiments, and various equivalent transformations of the technical solutions of the application can be made within the technical concept of the application, and these equivalent transformations all belong to the protection scope of the application.

Claims

1. A method for calculating compaction degree based on multiphysics joint inversion, characterized in that, include: Acquire multiphysics field observation data and soil compaction characteristic parameters of the area to be tested; Using a pre-built multiphysics forward model, a mapping relationship between state variables and multiphysics observation data is established. The state variables include at least water content and dry density. Based on the mapping relationship, a joint inversion iterative process is executed, and the state variables are updated based on the residual between the multiphysics observation data and the prediction output of the multiphysics forward model. The physical feasible region is constructed using the soil compaction characteristic parameters, and the physical feasible region is embedded as a constraint condition into the joint inversion iterative process to constrain the range of dry density and water content. In response to the convergence condition being met during the joint inversion iteration process, the dry density after inversion convergence is output, and the compaction degree is calculated using the dry density after inversion convergence and the preset maximum dry density.

2. The method according to claim 1, characterized in that, Multiphysics observation data includes ground-penetrating radar data and wave velocity correlation data; In addition, real-time operating data of the road roller is collected synchronously, including the amplitude of the equivalent normal force of the vibrating drum, the compaction speed, and the vibration frequency.

3. The method according to claim 1, characterized in that, Multiphysics forward modeling models include dielectric constant forward modeling and wave velocity forward modeling; The dielectric constant forward model is built based on a pre-configured complex refractive index model and is used to define the functional relationship between the dielectric constant and the volumetric water content derived from the dry density and water content. The wave velocity forward model is constructed based on the porosity elasticity theory or the empirically modified Hardin formula, and is used to define the functional relationship between wave velocity and dry density, saturation derived from water content, and preset vertical effective stress.

4. The method according to claim 1, characterized in that, Embedding the physical feasible region as a constraint in the joint inversion iterative process includes: Construct a joint loss function that includes a data fitting term and a constraint penalty term; Calculate the constraint violation of the dry density and moisture content of the current iteration step relative to the boundary of the physical feasible region; The constraint penalty term is calculated based on the constraint violation degree, and the value of the constraint penalty term is added to the joint loss function. The state variable is updated by minimizing the joint loss function.

5. The method according to claim 4, characterized in that, The constraint penalty term is constructed as a piecewise function of the normalized constraint violation rate, which includes a soft constraint region, a transition region, and a hard constraint region: In the soft-constraint region, the constraint penalty term is a quadratic function of the normalized constraint violation rate, which is used to allow the inverted solution to be adjusted within the confidence range of the physical feasible region; In the transition region, the constraint penalty term is a cubic polynomial function, which is used to maintain the continuity of the first derivative of the constraint penalty term and achieve a smooth transition from soft constraints to hard constraints. In the hard-constraint region, the constraint penalty term is an exponential function, used to impose an exponentially increasing penalty on inversion solutions that exceed a preset threshold.

6. The method according to claim 4, characterized in that, The physical feasible region is represented as a confidence region with statistical characteristics. Before performing the joint inversion iteration process, the confidence region is constructed, specifically as follows: Polynomial fitting was performed on the pre-acquired indoor compaction test data to obtain the centerline model of the compaction curve; Analyze the fitting residuals of indoor compaction test data relative to the centerline model and estimate the standard uncertainty of the fit; Introduce pre-stored field correction coefficients that reflect the variability of on-site construction, and use these field correction coefficients to expand the fitted standard uncertainty to obtain the comprehensive standard uncertainty; Based on the centerline model and the integrated standard uncertainty, the upper and lower boundaries of the confidence region are established.

7. The method according to claim 2, characterized in that, Embedding the physical feasible region as a constraint in the joint inversion iterative process includes: Construct a hard projection operator defined on the physical feasible region to compute the minimum distance projection from a point in the state space to the physical feasible region; In each update step of the joint inversion iteration process, the updated state variables are corrected using the hard projection operator; When the updated state variable falls outside the physical feasible region, it is forcibly mapped to the boundary or interior of the physical feasible region to obtain a corrected state variable that satisfies the physical constraints.

8. The method according to claim 7, characterized in that, The physically feasible region is an energy-adaptive region with dynamically variable boundaries, and the method further includes: Based on real-time operating data of the road roller, the compaction energy index at the current spatial location is calculated; The boundary shape and position of the physical feasible region are dynamically adjusted according to the compaction energy index, so that the dry density range covered by the physical feasible region shifts as the compaction energy index increases.

9. The method according to claim 8, characterized in that, The boundary shape and location of the physically feasible region are dynamically adjusted based on the compaction energy index, including: Using a pre-calibrated energy mapping model, a logarithmic growth relationship between maximum dry density and compaction energy index, and a linear offset relationship between optimum moisture content and compaction energy index are established. Based on the compaction energy index, the maximum dry density and optimum moisture content in the soil compaction characteristic parameters are updated in real time, and the center curve and boundary bandwidth of the physical feasible domain are reconstructed accordingly.

10. The method according to claim 1, characterized in that, Multiphysics observation data includes time-series observation data of the same detection area under different compaction passes. The method also includes introducing temporal monotonicity constraints during the joint inversion iteration process, specifically: Identify the order of compaction passes corresponding to time-series observation data; Based on the order of compaction passes, constraints are constructed to limit the estimated dry density value corresponding to the subsequent compaction pass from being lower than the estimated dry density value corresponding to the previous compaction pass, as a temporal monotonicity constraint. The temporal monotonicity constraint is added as an additional penalty term to the objective function of the joint inversion iterative process.

Images