Compaction degree measurement method based on multi-physical field joint inversion

Abstract

The application discloses a compaction degree measuring and calculating method based on multi-physical field joint inversion, comprising the following steps: obtaining multi-physical field observation data and geotechnical compaction characteristic parameters of a region to be detected; using a pre-constructed forward model to establish a mapping relationship between water content, dry density and observation data; performing joint inversion iteration, updating state variables based on observation residuals; in the iteration process, using geotechnological compaction characteristics to construct a physically feasible region, and embedding it as a constraint condition into the inversion process to constrain the value range of dry density and water content; in response to the joint inversion iteration process meeting a preset convergence condition, outputting the dry density after inversion convergence and calculating the compaction degree. The application solves the problems that the traditional inversion solution violates the soil compaction rule and cannot adapt to variable power compaction, and improves the physical authenticity and accuracy of compaction quality detection.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering quality testing technology, and in particular to a method for calculating compaction degree based on multi-physics field joint inversion. Background Technology

[0002] Compaction degree is a core indicator for evaluating the quality of embankment, roadbed, and dam construction projects, directly affecting the stability and durability of the engineering structure. With the increasing demand for refined engineering management, multi-physics field joint inversion using ground-penetrating radar (GPR) and seismic / surface wave technologies has become an important technical means to achieve wide-coverage, non-destructive testing of soil. This technology aims to solve the problem of multiple solutions in single geophysical methods by integrating complementary information from different physical fields, thereby more accurately inverting key parameters such as soil moisture content and dry density.

[0003] Current joint inversion methods are typically based on a data-driven framework, seeking the optimal solution by minimizing the residuals between observed data (such as radar electromagnetic wave travel time and surface wave dispersion curves) and forward modeling data. Most mainstream algorithms employ least squares or Bayesian inference, using empirical formulas (such as the CRIM model or Topp formula) to establish the relationship between physical parameters and geophysical responses. During the inversion process, only simple upper and lower bounds are typically set for the parameters, or common mathematical regularization techniques (such as Tikhonov regularization) are used to smooth the distribution of solutions, primarily focusing on satisfying the fitting accuracy of the observed data.

[0004] However, existing multiphysics joint inversion methods still have technical shortcomings in practical engineering applications, mainly reflected in the physical blindness of the inversion process and the static limitations of the constraint mechanism: Specifically, existing inversion algorithms lack an understanding of the essential physical laws of soil. As a special granular material, soil's dry density and moisture content strictly follow the compaction characteristics (Proctor curve). Traditional methods only pursue minimizing the data fitting residuals, resulting in the inverted dry density and moisture content combinations often falling outside the compaction curve, such as extremely high dry density at extremely low moisture content. This mathematically optimal but physically infeasible solution deviates from the true compaction state. In addition, existing constraint models cannot adapt to the varying work conditions of compaction in the field. In actual construction, the vibration frequency and speed of the road roller change in real time, causing the maximum dry density limit of the soil to fluctuate dynamically with the compaction energy. Existing technologies still use fixed parameters to constrain the inversion process, ignoring the change of energy on the soil's state limit, leading to a decrease in inversion accuracy under over- or under-compaction conditions. Summary of the Invention

[0005] The purpose of this invention is to provide a method for calculating compaction degree based on multi-physics field joint inversion, in order to solve the above-mentioned problems existing in the prior art.

[0006] Technical solution: A method for calculating compaction degree based on multi-physics field joint inversion, comprising:

[0007] Acquire multiphysics field observation data and soil compaction characteristic parameters of the area to be tested;

[0008] 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.

[0009] 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.

[0010] 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.

[0011] According to one aspect of this application, multiphysics observation data includes ground-penetrating radar data and wave velocity correlation data;

[0012] 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.

[0013] According to one aspect of this application, the multiphysics forward model includes a dielectric constant forward model and a wave velocity forward model;

[0014] 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.

[0015] 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.

[0016] According to one aspect of this application, the physical feasible region is embedded as a constraint in the joint inversion iterative process, including:

[0017] Construct a joint loss function that includes a data fitting term and a constraint penalty term;

[0018] 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;

[0019] 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.

[0020] According to one aspect of this application, the constraint penalty term is constructed as a piecewise function of normalized constraint violation, the piecewise function comprising 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 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;

[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 achieve a smooth transition from soft constraints to hard constraints.

[0023] 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.

[0024] According to one aspect of this application, the physically feasible region is represented as a confidence region with statistical characteristics, and the process further includes constructing the confidence region before performing the joint inversion iterative process, specifically:

[0025] Polynomial fitting was performed on the pre-acquired indoor compaction test data to obtain the centerline model of the compaction curve;

[0026] Analyze the fitting residuals of indoor compaction test data relative to the centerline model and estimate the standard uncertainty of the fit;

[0027] 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;

[0028] Based on the centerline model and the integrated standard uncertainty, the upper and lower boundaries of the confidence region are established.

[0029] According to one aspect of this application, the physical feasible region is embedded as a constraint in the joint inversion iterative process, including:

[0030] 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;

[0031] In each update step of the joint inversion iteration process, the updated state variables are corrected using the hard projection operator;

[0032] 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.

[0033] According to one aspect of this application, the physically feasible region is an energy-adaptive region with dynamically variable boundaries, and the method further includes:

[0034] Based on real-time operating data of the road roller, the compaction energy index at the current spatial location 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 shifts as the compaction energy index increases.

[0036] According to one aspect of this application, dynamically adjusting the boundary shape and location of the physically feasible domain based on the compaction energy index includes:

[0037] 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.

[0038] 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.

[0039] According to one aspect of this application, the multiphysics observation data includes time-series observation data of the same detection area under different compaction passes, and the method further includes introducing a time-series monotonicity constraint during the joint inversion iteration process, specifically:

[0040] Identify the order of compaction passes corresponding to time-series observation data;

[0041] 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.

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

[0043] Beneficial effects: This invention solves the problems of traditional inversion solutions violating the laws of soil compaction and being unable to adapt to compaction under varying work conditions, thereby improving the physical authenticity and accuracy of compaction quality detection. Attached Figure Description

[0044] Figure 1 A flowchart illustrating the steps of a method for calculating compaction degree based on multiphysics field joint inversion provided in this application embodiment.

[0045] Figure 2 A flowchart illustrating the steps of embedding the physical feasible region as a constraint in the joint inversion iterative process, as provided in this application embodiment.

[0046] Figure 3 A flowchart illustrating the steps for constructing a confidence region provided in this application embodiment.

[0047] Figure 4 A flowchart illustrating the steps of dynamically adjusting the boundary shape and position of the physically feasible domain based on the compaction energy index, as provided in this application embodiment. Detailed Implementation

[0048] 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.

[0049] It should be noted that the terms include and have, and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0050] like Figure 1 As shown, a method for calculating compaction degree based on multiphysics joint inversion includes the following steps:

[0051] Acquire multiphysics field observation data and soil compaction characteristic parameters of the area to be tested.

[0052] In this embodiment, acquiring multiphysics field observation data of the area to be detected specifically refers to collecting response data of different physical fields within the target engineering structure using on-site testing equipment. This multiphysics field observation data includes at least ground-penetrating radar (GPR) data and wave velocity correlation data. GPR data reflects the dielectric properties of the soil, while wave velocity correlation data reflects the elastic modulus or shear stiffness of the soil. Acquiring soil compaction characteristic parameters specifically refers to obtaining physical parameters describing the specific relationship between the dry density and moisture content of the soil. For example, key parameters of the compaction curve, including maximum dry density, optimum moisture content, and curve shape coefficients, can be obtained through indoor standard Proctor compaction tests. This provides a fundamental physical basis for subsequently constructing the physically feasible region.

[0053] 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.

[0054] Specifically, state variables refer to the unknown physical quantities that need to be solved during the inversion process. We choose moisture content w and dry density ρ. d These two parameters are core state variables because they directly determine the engineering properties of the soil and are closely related to its compaction characteristics. A multiphysics forward model is a mathematical model that can calculate theoretical observation data based on given state variables. For example, multiphysics forward models include dielectric constant forward models and wave velocity forward models. The dielectric constant forward model is used to calculate the theoretical dielectric constant at a given water content and dry density, thereby simulating the ground-penetrating radar response; the wave velocity forward model is used to calculate the theoretical wave velocity at a given water content and dry density. By establishing a mapping relationship, a forward prediction from engineering parameters to geophysical response is achieved.

[0055] Based on the mapping relationship, a joint inversion iterative process is performed, and the state variables are updated based on the residuals between the multiphysics observation data and the predicted output of the multiphysics 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 observed and predicted data. Specifically, based on the initial state variables, the predicted output is calculated using a multiphysics forward model; the residual between the predicted output and the actual acquired multiphysics observation data is calculated; and the state variables are updated using an optimization algorithm based on the magnitude of the residual and the gradient direction. Commonly used optimization algorithms include the Gauss-Newton method, the conjugate gradient method, or the L-BFGS-B algorithm. Through repeated iterations, the predicted data gradually approximates the observed data, thus inferring the true state variables.

[0057] A physical feasible region is constructed using soil compaction characteristic parameters, and this physical feasible region is embedded as a constraint condition in the joint inversion iterative process to constrain the range of dry density and moisture content.

[0058] Specifically, the physically feasible region refers to the reasonable range of values ​​within the moisture content-dry density state space that conforms to the physical laws of soil compaction. Unlike traditional inversion methods that rely solely on data fitting, this embodiment mandates that the state variables during the inversion process must fall within or near this physically feasible region. Optionally, the physically feasible region can be constructed based on a static compaction curve or a dynamic compaction curve that varies with compaction energy. Embedding the physically feasible region as a constraint can achieve soft constraints by adding a penalty function term to the objective function, or hard constraints by using a projection operator. This effectively eliminates solutions that, while fitting the observed data, violate the fundamental principles of geotechnical physics, thus improving the physical accuracy and reliability of the inversion results.

[0059] 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.

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

[0061] This embodiment achieves a deep integration of geophysical data and geotechnical test patterns by parameterizing the compaction characteristics of soil and rock into a physically feasible domain and embedding it into inversion iteration.

[0062] In one possible embodiment, the multiphysics observation data includes ground-penetrating radar data and wave velocity correlation data; and synchronously acquired real-time operating data of the road roller, including the amplitude of the equivalent normal force of the vibrating drum, the compaction speed, and the vibration frequency.

[0063] In this embodiment, the composition of multi-source data is defined. Ground-penetrating radar (GPR) data specifically refers to signals obtained by transmitting and receiving high-frequency electromagnetic waves using a GPR antenna, such as radar profiles, reflected wave travel time, or amplitude information. Wave velocity-related data specifically refers to physical quantities related to wave velocity obtained through seismic wave or surface wave detection methods, such as the phase velocity dispersion curve of Rayleigh waves. Synchronously acquiring real-time operating data of the road roller provides a basis for subsequent analysis of compaction energy. The equivalent normal force amplitude F of the vibratory drum... v The multiphysics observation data (v(t)) reflects the vertical force applied to the soil by the roller; the compaction speed v(t) reflects the duration of the energy application; and the vibration frequency f(t) reflects the frequency of energy loading. In addition, parameters such as the number of compaction passes N, the width of the vibratory roller B, and the pavement thickness H can also be collected. Furthermore, the multiphysics observation data needs to be spatially aligned and temporally synchronized to be uniformly mapped to the same detection grid.

[0064] In one specific implementation, the acquisition parameters for ground-penetrating radar data include: a center frequency range of 400MHz to 2GHz, selected according to the required detection depth; a sampling point count of no less than 512 points / channel; and a scanning spacing of 0.02m to 0.10m. The acquisition parameters for wave velocity-related data include: for surface wave detection, an excitation frequency range of 5Hz to 100Hz; a detector spacing of 0.5m to 2.0m; and a sampling frequency of no less than 1000Hz. The acquisition frequency for road roller operating condition data is no less than 10Hz to ensure the capture of the transient response of the vibrating wheel.

[0065] In an exemplary embodiment, the multiphysics forward model 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 the functional relationship between the dielectric constant and the volumetric water content derived from the dry density and water content.

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

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

[0068] Where ε represents the equivalent relative permittivity of the soil mixture; n represents the porosity, which can be expressed by the formula n=1-(ρ d / ρ s ) Calculate ρ d For dry density, ρ s ε is the density of soil particles; s θ represents the dielectric constant of soil particles; θ represents the volumetric water content, which can be expressed by the formula θ=w*(ρ d / ρ w ) Calculate, where w is the mass moisture content, ρ w ε is the density of water; w ε represents the dielectric constant of pore water. a Let represent the dielectric constant of air. State variables (w, ρ) are established. d The positive mapping from the observed variable (ε) to the observed variable (ε).

[0069] 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.

[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. This paper uses a saturation-corrected Hardin formula 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 properties; ρ d ρ is the dry density; ref The reference dry density is m; the density exponent is σ'. v For vertical effective stress; p a Atmospheric pressure; k is the stress exponent; S r Saturation can be expressed by the formula S. r =θ / n calculation; F(S) r ) is the saturation correction function, for example, F(S) r )=1-λ*S r λ is the weakening coefficient of water on the stiffness of the soil skeleton. It quantitatively describes the physical mechanism by which increased dry density leads to increased wave velocity, and increased water content (increased saturation) leads to decreased wave velocity. Optionally, the vertical effective stress is calculated based on a preset fill layer thickness and soil unit weight, or obtained from field measurements.

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

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

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

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

[0077] Where L data The fitting residuals of multiphysics observation data are typically expressed using weighted least squares, such as Li. data =||y gpr -h gpr (x)|| 2 / σ g 2 +||y seis -h seis (x)|| 2 / σ s 2 ;y gpr This refers to measured observation data from ground-penetrating radar (GPR); h gpr (x) represents the forward model prediction of ground-penetrating radar; σ g y represents the standard deviation of ground-penetrating radar observation data. seis This refers to measured observation data of seismic waves (seis); h seis (x) represents the forward model prediction of the seismic wave; σ s λ is the standard deviation of the seismic wave observation data. P λ represents the weight coefficient of the constraint penalty term; P(δ) represents the penalty function value calculated based on the constraint violation degree δ; w Indicates the weight of the water content prior test term; w prior Indicates the initial a priori value of water content; σ w denoted by , where represents the standard deviation of the prior a priori water content; w is the water content of the soil and rock mass. By minimizing this joint loss function, both the goodness of data fit and physical constraints can be considered simultaneously.

[0078] like Figure 3 As shown, in one possible implementation, the physical feasible region is represented as a confidence region with statistical characteristics. Before performing the joint inversion iterative process, the confidence region is constructed, specifically as follows:

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

[0080] Analyze the fitting residuals of indoor compaction test data relative to the centerline model and estimate the standard uncertainty of the fit;

[0081] 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;

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

[0083] For example, based on n sets of indoor compaction test data (wi , ρ d，i The centerline of the quadratic polynomial is fitted using the least squares method:

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

[0085] Where ρ d，fit (w) represents the fitted dry density value corresponding to moisture content w; a, b, and c are the fitting coefficients. Calculate the fitting residual e. i =ρ d，i -ρ d，fit (w i ), ρ d，i The measured dry density of the soil and rock mass obtained from the indoor compaction test of the i-th group; and the standard uncertainty of the fit σ is estimated. fit :

[0086] σ 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 factor κ is introduced to account for the differences between indoor and field environments, and the overall standard uncertainty σ is calculated. total :

[0088] σ total =sqrt(σ fit 2 +(κ*ρ dmax ) 2 );

[0089] The typical value range of κ is 0.02 to 0.05; ρ dmax To determine the maximum dry density, define the upper and lower boundaries of the confidence region:

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

[0091] Where ρ d，upper (w) represents the upper boundary value of the confidence region for dry density corresponding to moisture content w; ρ d，lower (w) represents the lower boundary value of the dry density confidence region corresponding to moisture content w; z α This is the confidence level coefficient; for example, 1.96 corresponds to a 95% confidence level.

[0092] 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.

[0093] Specifically, it quantifies the degree to which the current state variable deviates from the physically feasible region. For the dry density ρ of the current iteration... d Given the water content w, calculate its relative to the upper boundary ρ of the confidence region. d，upper The normalization constraint violation degree δ of (w). 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 relative to the lower boundary can also be calculated. This is illustrated using the example of excessive dry density (exceeding the upper boundary), as artificially 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 of normalized constraint violation, the piecewise function 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 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;

[0096] Alternatively, in the soft-constraint region, the constraint penalty term is a quadratic function of the normalized constraint violation degree, which is used to allow the inverted solution to be adjusted within a preset tolerance range within the confidence range of the physical feasible region;

[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 constraints to hard constraints.

[0098] 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.

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

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

[0101] Where β is the soft constraint coefficient, δ soft This is a soft constraint threshold, for example, 1.0. In the hard constraint region (δ>δ),... hard Using an exponential function:

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

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

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

[0105] To ensure that the function and its first derivative are within δ soft and δ hard For a continuous solution, the coefficients a3, a2, a1, and a0 need to be determined based on the boundary conditions. In particular, the cubic coefficient γ (i.e., a3) must satisfy the smoothness condition. Through a three-stage design, when the violation is small (within the confidence region), the penalty is small and smooth, allowing for data-driven fine-tuning; when the violation is too large (outside the confidence region), the penalty increases sharply, forcibly pulling the solution back to the physically reasonable range.

[0106] 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.

[0107] Specifically, substituting the calculated P(δ) into L total In the middle, L is calculated using optimization algorithms such as L-BFGS-B. total Regarding state variables w and ρ d The gradient is calculated, 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 beneficial to the stable convergence of the optimization algorithm. The final inversion solution can both fit the observation data well and satisfy the physical constraints of the soil compaction characteristics.

[0108] According to another aspect of this application, embedding the physically feasible region as a constraint in the joint inversion iterative process can also result in:

[0109] Construct a hard projection operator defined on the physical feasible region. The hard projection operator is used to compute the minimum distance projection from a point in the state space to the physical feasible region.

[0110] In this embodiment, the core tool for constraint enforcement is defined: the hard projection operator Π. ΩIt is a geometric operator that can convert any given state variable point (moisture content w, dry density ρ) into a geometric variable. d This maps to the nearest point within the physically feasible region Ω. This is equivalent to solving a constrained quadratic programming problem. Specifically, for the candidate state variable x updated in the current iteration step... update Find the corrected state variable x corrected This minimizes the Euclidean distance or weighted Mahalanobis distance between the two, and x corrected It must lie inside or on the boundary of the physically feasible region Ω. The introduction of the hard projection operator allows the inversion process to employ efficient algorithms such as the projection gradient descent method, ensuring convergence speed while strictly satisfying physical constraints.

[0111] In a further embodiment, the physically feasible region is an energy-adaptive region with dynamically variable boundaries, and the method further includes:

[0112] Based on real-time operating data of the road roller, the compaction energy index at the current spatial location is calculated.

[0113] Specifically, the compaction energy index E f It is a physical quantity that reflects the mechanical work performed per unit volume of soil. Based on working condition data, the compaction energy index can be calculated using 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 B represents the equivalent compaction energy per unit volume; H represents the width of the vibratory roller; F represents the current thickness of the fill layer; v (t) represents the amplitude of the equivalent normal force exerted by the vibratory roller on the soil at time t, which is usually determined by the static linear load and excitation force of the roller; v(t) represents the compaction speed at time t; dt represents the time derivative; integral represents the integral operation within the effective compaction time window. By integrating the discrete working condition parameters of the roller, such as vibration frequency, speed, and excitation force, into a unified energy index, a foundation is laid for establishing the physical correlation between energy and density.

[0116] 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.

[0117] like Figure 4 As shown, in a preferred implementation, dynamically adjusting the boundary shape and position of the physically feasible region based on the compaction energy index includes:

[0118] 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.

[0119] 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 based on the updated maximum dry density and optimum moisture content.

[0120] In this embodiment, a quantitative mapping model from energy to compaction characteristics is established. According to soil mechanics principles, as compaction work (energy) increases, the compaction curve of the soil changes shape: the maximum dry density increases, and the corresponding optimum moisture content decreases (shifts to the upper left). This pattern can be described using specific empirical formulas. For the maximum dry density, a logarithmic growth model is adopted:

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

[0122] Where ρ dmax (E f ) represents the maximum dry density at the current energy; ρ dmax_std k represents the maximum dry density at standard compaction energy. ρ is the density growth coefficient, reflecting the soil's sensitivity to energy; log is the natural logarithm function; E ref The reference energy value is used. For the optimum moisture content, a linear migration model is employed:

[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 k represents the optimum moisture content at the standard compaction energy. w This is the moisture content offset coefficient. Based on real-time calculated ρ dmax (E f ) and w opt (E f The current parabolic impact centerline is regenerated, and based on this centerline, combined with the preset bandwidth Δ, the current dynamic physical feasible region Ω(E) is constructed. fFor example, when an increase in the number of compaction passes by a road roller in a certain area is detected, leading to E... f When the dry density is significantly increased, the algorithm will automatically raise the upper limit of the dry density constraint in that region, allowing the inversion to obtain higher dry density results and avoiding the peak-shaving error that may be caused by traditional static constraints.

[0125] 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 variables fall outside the physical feasible region, the updated state variables are forcibly mapped to the boundary or interior of the physical feasible region to obtain corrected state variables that satisfy the physical constraints.

[0126] Specifically, it describes the concrete actions taken to enforce the constraints. In the k-th step of the inversion iteration, an unconstrained intermediate solution x is obtained solely based on the gradient descent direction of the data residuals. k+0.5 Determine x k+0.5 Does the current energy E meet the requirements? f The corresponding feasible region condition is g(x) ≤ 0; where g(x) is the physical constraint function on the state variable x. If satisfied, then directly let x... k+1 =x k+0.5 If the conditions are not met, for example, if the calculated dry density is much higher than the theoretical limit achievable with the current energy, then the hard projection operator is used to solve the following optimization problem:

[0127] min x ||xx k+0.5 || 2 st x∈Ω(E f );

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

[0129] Preferably, the method further includes anomaly detection based on the execution process of the hard projection operator, specifically:

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

[0131] In this embodiment, an intelligent diagnostic mechanism based on inversion process data is provided. Projection distance d projIt is an indicator that measures the degree of conflict between observational data and physical laws. The specific calculation formula is:

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

[0133] Where w corrected and ρ dcorrected This is the value after projection correction; w update and ρ dupdate This is the value before correction. If d proj A consistently large value indicates that the observed data strongly tends to produce solutions that violate physical laws, such as exhibiting extremely high dielectric constants and wave velocities at low energies. This usually signifies the failure of the pre-set physical model, such as a sudden change in soil properties, excessively high rock content, or sensor malfunction. By setting a continuous counting threshold K, for example, K=5, when the projected distance exceeds a safe threshold d for 5 consecutive iterations... th When this happens, the system will automatically trigger an alarm, prompting engineers to review the area, instead of forcibly outputting results that have been corrected by the projection algorithm but may be distorted.

[0134] Furthermore, in practical applications, the following boundary conditions and anomaly handling should also be considered: If ground-penetrating radar data or wave velocity data is missing at a certain measurement point, spatial interpolation methods can be used to complete it, or the measurement point can be skipped and only the inversion results of valid measurement points can be output. If the observed data significantly exceeds the physically reasonable range, such as a dielectric constant less than 1 or a negative wave velocity, it should be removed during the preprocessing stage and marked in the detection report. If the convergence condition is not met after reaching the preset maximum number of iterations, such as 100 iterations, the current optimal estimate should be output with a marker indicating reduced confidence. For material anomaly areas identified by projection distance anomaly detection, it is recommended to combine traditional methods such as core drilling for verification.

[0135] This embodiment utilizes the feature of modern intelligent road rollers that can provide real-time feedback of energy data to construct a physically feasible domain that dynamically changes with compaction energy. It also forces the inversion solution to satisfy this dynamic physical constraint through a hard projection operator, making it suitable for scenarios involving compaction quality testing of fill bodies with varying work conditions and large thicknesses.

[0136] According to one aspect of this application, the multiphysics observation data includes time-series observation data of the same detection area under different compaction passes, and the method further includes introducing a time-series monotonicity constraint during the joint inversion iteration process, specifically:

[0137] Identify the order of compaction passes corresponding to time-series observation data;

[0138] 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.

[0139] The temporal monotonicity constraint is added as an additional penalty term to the objective function of the joint inversion iterative process.

[0140] In this embodiment, a physically irreversible constraint is introduced for multi-pass compaction scenarios. Soil compaction is a densification process; under normal construction conditions, as the number of compaction passes (t) increases, the dry density (ρ) increases. d (t) should be monotonically non-decreasing. To handle measurement errors, it is preferable to use a quadratic penalty function with a dead zone to construct the monotonicity constraint term L. mono The specific formula is as follows:

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

[0142] Where ∑ represents the summation over all adjacent iterations t and t+1; Φ t A quantity that violates monotonicity is defined as ρ. d (t)-ρ d (t+1). The penalty function Q(Φ) is defined as: if Φ t ≤Φ tol Then Q(Φ t If )=0; if Φ t >Φ tol Then Q(Φ t )=μ*(Φ t -Φ tol ) 2 ;Φ tol L is the tolerance threshold (dead zone width), used to tolerate small measurement fluctuations; μ is the penalty coefficient. This is achieved by minimizing L... mono The algorithm penalizes non-physical inversion paths whose density is significantly lower in subsequent passes than in previous passes, forcing the inversion results to maintain a reasonable growth trend over time. In terms of solution strategy, a block coordinate descent method can be used to optimize the parameters of each pass sequentially over time.

[0143] In one embodiment of this application, multiphysics observation data covers a predetermined number of spatially distributed measurement points. The method further includes introducing spatial continuity constraints during the joint inversion iteration process, specifically:

[0144] Construct an adjacency graph of each measuring point based on its spatial coordinates, and calculate the graph Laplace matrix of the adjacency graph.

[0145] A graph regularization term is constructed using the graph Laplacian matrix. This graph regularization term is used to constrain the difference in dry density estimates between adjacent measurement points.

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

[0147] In this embodiment, a local smoothing constraint is introduced for the spatial distribution inversion results. During normal construction, the compaction characteristics of the soil vary gradually with spatial location and should not exhibit sawtooth-like fluctuations. This can be achieved by constructing a graph regularization term R. graph This is used to measure the differences between adjacent nodes. Specifically, an adjacency graph G=(V,E) of all measured points i is constructed, where V is the set of measured points and E is the set of edges. The adjacency weight W between measured points i and j is defined. ij Adjacency weight and measurement point distance d ij Negative correlation, for example, can be calculated using a Gaussian kernel function. Calculate the graph Laplacian matrix L. G =DW, where D is the degree matrix. Define the graph regularization term:

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

[0149] Where ρ d L is the dry density vector for all measuring points; G The graph is a Laplace matrix; T This is a transpose. R is... graph The graph regularization term is added to the objective function as an additional term to penalize large abrupt changes in dry density between adjacent measurement points, thereby improving the spatial continuity of the inversion results. To avoid over-smoothing, the weight coefficient λ of the graph regularization term... graph It can decay with the number of iterations, for example, λ. graph_k =λ graph_0 *(0.95) k ;where λ graph_k The weight coefficients of the graph regularization term in the k-th iteration are given.

[0150] In some embodiments, the method further includes introducing an identifiability constraint during the joint inversion iteration process to suppress the equivalence of solutions for dry density and moisture content. Specifically, this involves: calculating the sensitivity matrix of multiphysics observation data to dry density and moisture content based on a multiphysics forward model; constructing an identifiability regularization term based on the coupling term of the sensitivity matrix; and superimposing the identifiability regularization term into the objective function of the joint inversion iteration process to penalize inversion paths that result in similar physical field responses but significant differences in state variables.

[0151] Specifically, regarding moisture content w and dry density ρ 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 at point P; trans (δ hard P is the transition region penalty function. trans (δ) on the right boundary δ hard The function value at point P; hard (δ hard P is the penalty function for the hard constraint region. hard (δ) at δ hard The function value at point P' trans (δ hard P is the transition region penalty function. trans (δ) at δ hard The first derivative value at point P' hard (δ hard P is the penalty function for the hard constraint region. hard (δ) at δ hard The value of the first derivative at that point.

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

[0171] In this embodiment, equations A and B are solved simultaneously to obtain γ and k2. Rearranging equation B: 6.75*γ + 3*k2 = 190.0; Rearranging equation A: 3.375*γ + 2.25*k2 = 80.0. Solving this system of linear equations yields the specific values ​​of γ and k2. For example, if γ = 25.0, then the coefficient of this cubic term is the key parameter for achieving a smooth transition. In the actual inversion code, this calculation process is automatically executed during the initialization phase, ensuring that gradient calculations during the inversion process do not jump, thus guaranteeing the stability of the L-BFGS-B algorithm.

[0172] This embodiment solves the mathematical problem of ensuring that the penalty function maintains the continuity of the first derivative (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 the initial conditions of a roadbed detection area are as follows: the equivalent dielectric constant ε measured by ground penetrating radar is... obs =12.5; Shear wave velocity V obtained from surface wave detection s，obs =185m / s; the benchmark maximum dry density ρ determined by the indoor standard compaction test. dmax_std =1.85g / cm 3 ;Optimal moisture content w opt =14.2%; soil particle density ρ s =2.68g / cm 3 ; Set initial state variables: dry density ρ d 0 =1.70g / cm 3 Moisture content w 0=15.0%. In the first iteration, the predicted dielectric constant was calculated using the CRIM model, derived from ρ. d 0 =1.70g / cm 3 and w 0 =15.0%, the calculated porosity n = 1 - 1.70 / 2.68 = 0.366, and the volumetric water content θ = 0.15 × 1.70 / 1.0 = 0.255. Substituting into the CRIM model, the predicted dielectric constant ε = 11.8 is obtained. The predicted wave velocity is calculated using the Hardin formula, assuming a vertical effective stress σ'. v =30kPa, substituting into the formula, we obtain the predicted shear wave velocity V. s =178m / s. Calculate the residuals: dielectric constant residual = 12.5 - 11.8 = 0.7; wave velocity residual = 185 - 178 = 7m / s. Update the state variables based on gradient descent (step size α = 0.1) to obtain the unconstrained update value ρ. d 0.5 =1.73g / cm 3 w 0.5 =14.6%. Based on the centerline ρ of the compaction curve. d，fit (14.6%) = 1.82 g / cm³ 3 upper boundary ρ d，upper =1.82 + 1.96 × 0.03 = 1.88 g / cm³ 3 Due to ρ d 0.5 =1.73 < 1.88, the constraint is satisfied, no correction is needed. After the first iteration: ρ d 1 =1.73g / cm 3 w 1 =14.6%. After several iterations: convergence yields the inversion result ρ. d *=1.78g / cm 3 w*=14.3%. The calculated compaction degree K=1.78 / 1.85×100%=96.2%. The results show that this embodiment can achieve accurate estimation of compaction degree while ensuring that the inversion solution falls within the physical feasible region of the compaction curve.

[0174] In summary, the compaction degree calculation method based on multi-physics joint inversion includes: acquiring multi-physics observation data and soil compaction characteristic parameters of the area to be tested; establishing the mapping relationship between water content and dry density and observation data using a pre-constructed forward model; performing joint inversion iteration and updating state variables based on observation residuals; constructing a physically feasible region using soil compaction characteristics during the iteration process and embedding it as a constraint condition into the inversion process, forcing the state variables to fall within the range that conforms to physical laws through a penalty function or projection operator; and finally outputting the dry density after inversion convergence and calculating the compaction degree.

[0175] This invention introduces a constraint mechanism based on the compaction characteristics of soil and rock. By parameterizing the compaction curves obtained from indoor tests, a physically feasible region of state variables is constructed. A penalty function that switches between soft and hard constraints, or a geometric hard projection operator, is used to force the inversion iteration process to be confined within an effective range that conforms to geotechnical laws. This effectively eliminates spurious solutions that fit observational data but violate physical principles, improving the physical accuracy of the test results. A dynamic inversion strategy based on energy adaptation is proposed. By calculating the compaction energy index of the road roller in real time, the boundary shape of the physically feasible region is dynamically adjusted, such as the maximum dry density increasing logarithmically with energy. This overcomes the bias of traditional static constraints under over- or under-compaction conditions, enabling the detection algorithm to intelligently adapt to complex on-site construction processes and ensuring inversion accuracy under variable energy environments.

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

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. The process of embedding the physical feasible region as a constraint into the joint inversion iterative process includes: constructing a joint loss function containing data fitting terms and constraint penalty terms; calculating the constraint violation degree of the dry density and moisture content of the current iteration step relative to the boundary of the physical feasible region; calculating the value of the constraint penalty term based on the constraint violation degree, and superimposing the value of the constraint penalty term into the joint loss function; and updating the state variables by minimizing the joint loss function. The constraint penalty term is constructed as a piecewise function of the normalized constraint violation degree. The piecewise function 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 degree, 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 realize a smooth transition from soft constraints to hard constraints; In the hard constraint region, the constraint penalty term is an exponential function, which is used to impose an exponentially increasing penalty on the inverted solution that exceeds a preset threshold.

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, 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.

5. 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.

6. The method according to claim 5, 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.

7. The method according to claim 6, 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.

8. 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