A method for constructing an existing asphalt pavement damage model based on damage density

Abstract

The application discloses a kind of existing asphalt pavement damage model construction methods based on damage density, specifically related to the numerical simulation technical field of asphalt pavement structure damage assessment, for solving the problem that existing method is difficult to accurately, efficiently quantitatively evaluate and predict the internal damage state of asphalt pavement structure in service;Method by obtaining pavement material and geometric parameters and deflection data, a finite element numerical simulation model is established, the dynamic response under different damage states is simulated, and the phase space reconstruction is used to extract the attractor structure characteristics of dynamic system, and then the continuous evolution geometric constraint path of damage density in phase space is constructed, and finally the damage density of surface layer and base layer is calculated based on the path constraint combined with deflection data by optimization inversion algorithm, according to which the internal damage state of pavement is determined, realizes the quantitative evaluation of mechanism and data fusion of internal damage of asphalt pavement structure, and provides accurate basis for pavement maintenance decision.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for structural damage assessment of asphalt pavement, and in particular to a method for constructing an existing asphalt pavement damage model based on damage density. Background Technology

[0002] In the field of asphalt pavement structural damage assessment, existing technologies mainly rely on laboratory tests to obtain material parameters and establish theoretical models. Specifically, material performance tests are conducted in the laboratory through core sampling to obtain parameters such as dynamic modulus and resilient modulus, and based on these, mechanical response models of the pavement structure, such as finite element models, are constructed to simulate its performance.

[0003] However, due to significant differences between indoor test conditions and actual pavement service conditions in terms of specimen size, environmental complexity, and load simulation, models built based on such tests cannot accurately reflect the mechanical performance degradation process of existing pavement structures under real complex environments and load coupling. Furthermore, existing methods for obtaining pavement damage status, such as manual observation or image processing, are inefficient, while non-destructive testing technologies such as ground-penetrating radar can only identify existing internal defects and cannot quantitatively assess the damage status of the pavement structure, let alone predict its future performance evolution trend. Therefore, they cannot provide accurate damage degree criteria for pavement maintenance decisions. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a method for constructing a damage model of existing asphalt pavement based on damage density.

[0005] The technical solution of the present invention to solve the above-mentioned technical problems is as follows:

[0006] A method for constructing a damage model for existing asphalt pavements based on damage density, comprising:

[0007] S1. Obtain the material and geometric parameters of the target asphalt pavement structure, as well as the deflection data of the target asphalt pavement structure under its current service condition.

[0008] S2. Establish a finite element numerical simulation model of the target asphalt pavement structure based on material parameters and geometric parameters;

[0009] S3. The dynamic response of the target asphalt pavement structure under different preset damage states is simulated using a finite element numerical simulation model, and the dynamic system attractor structure features corresponding to each preset damage state are extracted using the phase space reconstruction method.

[0010] S4. Based on the attractor structure characteristics of the dynamic system, construct a geometric constraint path to characterize the continuous variation law of damage density in phase space;

[0011] S5. Based on geometric constraint paths combined with deflection data, the surface layer damage density and base layer damage density of the target asphalt pavement structure are calculated using an optimization inversion algorithm with geometric constraint paths as constraints.

[0012] S6. Based on the surface layer damage density and base layer damage density, determine the internal damage state of the target asphalt pavement structure.

[0013] Furthermore, S1 includes:

[0014] The dynamic modulus of the existing asphalt pavement surface layer and the resilient modulus of the base layer were obtained as material parameters through core sampling and indoor testing.

[0015] Ground-penetrating radar was used to detect the thickness of each structural layer of the existing asphalt pavement and the bonding state between each layer as geometric parameters.

[0016] The location and characteristic values ​​of the deflection basin of the existing asphalt pavement under its current service condition are obtained by using a falling weight deflectometer as deflection data.

[0017] Furthermore, S2 includes:

[0018] Input the material and geometric parameters into the multiphysics simulation platform;

[0019] In the multiphysics simulation platform, a finite element numerical simulation model of the existing asphalt pavement structure is established based on material parameters and geometric parameters. A viscoelastic damage constitutive model is applied to the asphalt surface layer in the finite element numerical simulation model, and a stress-related elastic damage constitutive model is applied to the base layer.

[0020] A falling weight deflection load was applied to the finite element numerical simulation model after constituting the model to establish a finite element numerical simulation model that can simulate the dynamic deflection time history response under different damage states.

[0021] Furthermore, S3 includes:

[0022] Using a finite element numerical simulation model, multiple preset damage states can be set by adjusting the damage density parameter in the finite element numerical simulation model.

[0023] Run the finite element numerical simulation model under each preset damage state to obtain the corresponding dynamic deflection time history response data.

[0024] For each preset damage state, the dynamic deflection time history response data is reconstructed using the delayed coordinate method;

[0025] In the reconstructed phase space, invariants used to characterize the attractor's geometry and dynamic properties are calculated to obtain the dynamic system attractor structural features corresponding to each preset damage state.

[0026] Furthermore, invariants used to characterize the attractor's geometry and dynamics are calculated, including: calculating the correlation dimension of the reconstructed phase space trajectory to characterize the geometry, and calculating the maximum Lyapunov exponent and recursive quantization analysis index to characterize the dynamics, thereby obtaining the attractor structural characteristics of the dynamic system.

[0027] Furthermore, S4 includes:

[0028] Collect the structural features of the dynamic system attractor corresponding to all preset damage states, and associate each preset damage state and its corresponding damage density parameter value with the structural features of the dynamic system attractor.

[0029] Using the attractor structural features of the dynamic system as coordinates, feature points corresponding to each preset damage state are marked in phase space;

[0030] By using curve fitting, feature points in phase space corresponding to different damage density parameter values ​​are fitted and connected to generate smooth and continuous paths.

[0031] The generated smooth continuous path is defined as a geometrically constrained path that characterizes the continuous variation of damage density in phase space.

[0032] Furthermore, using the structural features of the attractor of the dynamic system as coordinates, feature points corresponding to each preset damage state are marked in the phase space. This includes: constructing a multidimensional phase space using multiple structural features of the attractor of the dynamic system, such as the correlation dimension and Lyapunov index, as coordinate axes, and using a set of specific feature values ​​corresponding to each preset damage state as coordinates for positioning and marking in the multidimensional phase space.

[0033] Furthermore, S5 includes:

[0034] The previously obtained deflection basin location and deflection basin characteristic values ​​are used as deflection data input, and an objective function is constructed with surface layer damage density and base layer damage density as optimization variables.

[0035] The geometric constraint path is transformed into mathematical constraints on the optimization variables, where the geometric constraint path is achieved by defining the mapping relationship between damage density and the attractor structural characteristics of the dynamic system in phase space;

[0036] An artificial intelligence optimization algorithm is used to iteratively update the surface layer damage density and base layer damage density parameters within the solution space that satisfies mathematical constraints, until the objective function converges to a preset convergence threshold, and the final calculation results of surface layer damage density and base layer damage density are output.

[0037] Furthermore, the objective function is used to measure the difference between the inverted predicted deflection data and the actual deflection data.

[0038] Furthermore, S6 includes:

[0039] The surface layer damage density and base layer damage density at different locations, calculated in advance, will be used as discrete data points;

[0040] The discrete data points of surface layer damage density and base layer damage density are processed by spatial interpolation method to generate continuous damage density distribution fields covering the surface layer and base layer of the target asphalt pavement structure, respectively.

[0041] Based on the generated continuous damage density distribution field, the regions in the surface layer and base layer whose damage density exceeds the preset damage density threshold and their damage density values ​​are identified.

[0042] The area exceeding the preset damage density threshold and its combined damage density value are determined as the internal damage state of the target asphalt pavement structure.

[0043] The beneficial effects of this invention are:

[0044] 1. By constructing a system that integrates finite element simulation, dynamic system analysis, and constrained optimization inversion, a deep quantitative assessment of the internal damage state of existing asphalt pavement structures is achieved. Attractor structural features that characterize the nonlinear dynamic properties of the structure are extracted from the dynamic deflection response of the pavement, and a continuous evolution path of damage density in phase space is constructed based on this. This breaks through the limitations of traditional methods that rely solely on deflection amplitude or simple statistical features for inversion. It can extract deep state information closely related to the internal damage mechanism of the material from the complex dynamic response of the pavement system under load, thereby realizing independent, refined decoupling and quantitative inversion of the damage density of asphalt surface layer and base layer.

[0045] 2. By introducing geometric constraint paths based on the evolutionary laws of dynamic systems and using them as the core constraint conditions for optimizing the inversion algorithm, it is ensured that the inversion process not only relies on the matching of measured data and mechanical models, but also follows the physical continuity and regularity of damage development itself. This inversion paradigm that integrates mechanism and data effectively overcomes the problems of multiple solutions and large uncertainties in inversion results caused by model simplification or insufficient data in traditional methods. It significantly improves the physical rationality and reliability of damage assessment results. The resulting continuous damage density distribution field and the internal damage state determined accordingly provide direct, accurate and spatially visualized quantitative basis for pavement maintenance decisions. The entire process relies on numerical simulation and intelligent algorithms and belongs to the deep computational analysis and decision support of pavement structural performance data. Attached Figure Description

[0046] Figure 1 This is a flowchart of a method for constructing a damage model of existing asphalt pavement based on damage density according to the present invention.

[0047] Figure 2The flowchart for constructing a geometric constraint path to characterize the continuous variation of damage density in phase space is provided for this invention. Detailed Implementation

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

[0049] Example: Figure 1 This invention presents a method for constructing a damage model of existing asphalt pavement based on damage density, comprising:

[0050] S1. Obtain the material and geometric parameters of the target asphalt pavement structure, as well as the deflection data of the target asphalt pavement structure under its current service condition.

[0051] S2. Establish a finite element numerical simulation model of the target asphalt pavement structure based on material parameters and geometric parameters;

[0052] S3. The dynamic response of the target asphalt pavement structure under different preset damage states is simulated using a finite element numerical simulation model, and the dynamic system attractor structure features corresponding to each preset damage state are extracted using the phase space reconstruction method.

[0053] S4. Based on the attractor structure characteristics of the dynamic system, construct a geometric constraint path to characterize the continuous variation law of damage density in phase space;

[0054] S5. Based on geometric constraint paths combined with deflection data, the surface layer damage density and base layer damage density of the target asphalt pavement structure are calculated using an optimization inversion algorithm with geometric constraint paths as constraints.

[0055] S6. Based on the surface layer damage density and base layer damage density, determine the internal damage state of the target asphalt pavement structure.

[0056] S1. Obtain the material and geometric parameters of the target asphalt pavement structure, as well as the deflection data of the target asphalt pavement structure under its current service condition. Specifically, this is implemented as follows:

[0057] The dynamic modulus of the existing asphalt pavement surface layer and the resilient modulus of the base layer are obtained as material parameters through core drilling and laboratory testing. Specifically, at representative locations of the target asphalt pavement structure, such as 1 to 3 points per kilometer per lane, a complete cylindrical specimen with a diameter of 150 mm, containing both the asphalt surface layer and the base layer, is drilled using a core drill. This cylindrical specimen is then sent to the laboratory, where it is processed into standard beam-type specimens or cylindrical specimens that meet experimental requirements using a cutting and grinding machine. For the asphalt surface layer material, a four-point bending fatigue testing system is used to conduct dynamic mechanical property tests on the standard beam-type specimens under test conditions, such as 20 degrees Celsius, a loading frequency of 10 Hz, and 60 micro-strain levels. By measuring the stress and strain response of the specimen under different loading cycles, and calculating the dynamic modulus of the asphalt surface layer material under specific conditions based on the definition of dynamic modulus (the ratio of stress amplitude to strain amplitude), the unit is megapascal (MPa). For the base material, a uniaxial unconfined compressive strength test was conducted on the processed cylindrical specimens. During the test, the springback deformation of the specimens was measured, and the springback modulus was calculated using the formula for calculating the springback modulus (compressive strength divided by springback strain), with the unit being megapascals. These dynamic modulus and springback modulus values, directly measured through indoor tests, serve as the necessary surface layer and base material parameters for establishing the finite element numerical simulation model in subsequent steps.

[0058] Ground-penetrating radar (GPR) is used to obtain the thickness of each structural layer and the bonding state between layers of existing asphalt pavement as geometric parameters. In practice, a GPR antenna with a center frequency of 1.6 GHz is used, and multiple detection lines are arranged along the driving direction of the target asphalt pavement structure, such as the lane centerline and wheel tracks. The GPR antenna emits high-frequency electromagnetic pulses into the pavement structure and receives reflected wave signals from the interfaces of different material layers. The received reflected wave signals are preprocessed, including gain adjustment, background denoising, and filtering, and then time-depth conversion is performed. Time-depth conversion requires the relative permittivity of each structural layer material as input. This constant can be obtained by back-calculating from measured data on a calibrated road section with known thicknesses. Specifically, the propagation speed is obtained by dividing the known thickness by half the two-way travel time of the electromagnetic wave in that layer and then multiplying by the speed of light, thus calculating the permittivity. The processed radar image clearly shows the reflection interfaces of different layers. By reading the two-way travel time between adjacent interfaces and combining it with the electromagnetic wave propagation speed of the material in that layer, the thickness of each structural layer, such as the asphalt surface layer, base layer, and subbase layer, can be calculated, with the thickness measured in millimeters. Furthermore, the interlayer bonding state can be evaluated by analyzing the waveform characteristics, energy intensity, and continuity of reflected waves from the interlayer interface. For example, when interlayer bonding is good, the reflected wave signal is weak and continuous; when there are interlayer voids or poor bonding, a strong, discontinuous reflected signal will be generated at the interface. By identifying and interpreting these characteristics, areas with poor interlayer bonding can be qualitatively identified and located, thereby obtaining the geometric parameter of interlayer bonding state.

[0059] The location and characteristic values ​​of the deflection basin of an existing asphalt pavement under its current service condition are obtained using a falling weight deflectometer. The specific operation is as follows: The bearing plate of the falling weight deflectometer is placed on the surface of the target asphalt pavement structure at the point to be measured. A 120 kg falling weight is lifted and released by a hydraulic system, allowing it to fall freely from a set height and impact the bearing plate, thus applying a transient half-sine pulse load to the pavement. In the sensor array of the falling weight deflectometer, a high-precision displacement sensor located at the center of the bearing plate records the dynamic deflection time history curve of that point, i.e., the center deflection value. Simultaneously, a series of displacement sensors arranged at fixed intervals along the radial direction of the bearing plate, for example at distances of 0 mm, 200 mm, 300 mm, 450 mm, 600 mm, 900 mm, and 1200 mm from the center point, synchronously record the dynamic deflection time history curves at different distances. Each load measurement corresponds to a deflection basin location, which is precisely located and its coordinates are recorded using a GPS or surveying methods. Peak deflection values ​​are extracted from the time-history curves recorded by each sensor. The peak deflection value at the center point of the bearing plate is recorded as the central deflection. The peak deflection values ​​of each radial sensor are arranged by distance to form the deflection basin for that measuring point. The characteristic values ​​of the deflection basin are calculated from this series of peak deflection values, including the diameter and area of ​​the deflection basin. The diameter of the deflection basin is defined as the radial distance from the center of the bearing plate to where the deflection value decays to, for example, 50% of the central deflection value. The area of ​​the deflection basin is obtained by integrating the distribution curve of the peak deflection value along the radial direction using numerical integration methods, such as the trapezoidal integral method. These deflection basin location coordinates and the calculated characteristic values ​​of the deflection basin together constitute the deflection data characterizing the overall stiffness and structural integrity of the target asphalt pavement structure under the current load. All data obtained through core sampling, indoor tests, ground-penetrating radar detection, and falling-weighted deflectometer detection are archived according to a unified spatial reference coordinate system to ensure that material parameters, geometric parameters, and deflection data are accurately correlated in spatial location, providing complete and consistent input information for subsequent steps.

[0060] S2. Establish a finite element numerical simulation model of the target asphalt pavement structure based on material and geometric parameters. The specific implementation is as follows:

[0061] Material and geometric parameters are input into the multiphysics simulation platform. Specifically, the multiphysics simulation software is launched, and a new model file is created through its graphical user interface or via script commands. The dynamic modulus of the asphalt surface layer and the resilient modulus of the base layer, previously obtained through core sampling and laboratory tests, are input into the software's material library as material properties, with both dynamic and resilient modulus values ​​in MPa. Simultaneously, the thickness values ​​of each structural layer of the existing asphalt pavement, previously obtained through ground-penetrating radar detection, and the qualitative evaluation results regarding the interlayer bonding state, are used as input data for geometric modeling. The unit for each structural layer thickness is millimeters. The interlayer bonding state evaluation results can be used to guide the setting of contact interfaces in subsequent models; for example, areas evaluated as having good bonding are set as fully bound constraints, while areas evaluated as having poor bonding are set as frictional contacts that allow separation.

[0062] In a multiphysics simulation platform, a finite element numerical simulation model of an existing asphalt pavement structure is established based on material and geometric parameters. Specifically, according to the input thickness values ​​of each structural layer, the software's geometric modeling module creates three-dimensional solid geometries representing the asphalt surface layer, base course, subbase, and subgrade layer layer by layer by drawing rectangles and performing extrusion operations. These geometries are stacked vertically, and their thickness dimensions are strictly set according to the input thickness values. Subsequently, the created geometries are meshed using finite element methods. Hexahedral elements are selected as the mesh type, and the element size is set according to the requirements of computational accuracy and efficiency; for example, a smaller, finer mesh is used near the load-bearing area, while a larger, coarser mesh is used in areas further away. After meshing, the material properties stored in the material library are assigned to the corresponding geometric components; that is, the dynamic modulus property of the asphalt surface layer is assigned to the geometry representing the surface layer, and the resilient modulus property of the base course is assigned to the geometry representing the base course. For the material parameters of the subbase and subgrade, if they are not directly obtained through testing, they can be defined using engineering experience values.

[0063] A viscoelastic damage constitutive model is applied to the asphalt pavement in the finite element numerical simulation model. This constitutive model describes the time-dependent behavior and damage evolution of the asphalt mixture under load. In the software implementation, a viscoelastic material model is selected from the material model library, and damage criteria are superimposed. The viscoelastic behavior is described using a generalized Maxwell model, which consists of a series of parallel Maxwell elements, each containing an elastic spring and a viscous damper. Model parameters include the instantaneous elastic modulus, the long-term elastic modulus, and the relaxation time and modulus weighting coefficients for each Maxwell element. These parameters need to be calibrated through laboratory tests of the asphalt mixture. The specific method for obtaining these parameters is as follows: First, a uniaxial static creep test of the asphalt mixture is conducted. A constant compressive stress is applied to a standard specimen, and the complete curve of strain growth over time is continuously measured and recorded; this curve is called the creep compliance curve. Subsequently, based on the theoretical formula of the generalized Maxwell model, the macroscopic creep compliance is expressed as the sum of a series of negative exponential functions with the natural constant e as the base. Each exponential term corresponds to a Maxwell unit, with its time constant being the relaxation time and its coefficient being the modulus weighting coefficient. The experimentally measured creep compliance curve data is input into a parameter identification algorithm, such as using nonlinear least squares curve fitting. The algorithm aims to adjust the relaxation time spectrum and corresponding modulus weighting coefficients in the formula to minimize the sum of squared errors between the creep curve calculated by the theoretical model and the experimental curve. After iterative optimization, the algorithm outputs a set of optimal relaxation times and modulus weighting coefficients. The instantaneous and long-term elastic moduli can also be determined from the fitted curves. Damage evolution is described by a damage variable that varies from 0 to 1, where 0 represents no damage and 1 represents complete failure. The evolution law of the damage variable is calibrated through cyclic loading fatigue tests on asphalt mixtures. In strain-controlled fatigue tests, the curve of the specimen modulus decreasing with the number of loading cycles is recorded. The damage evolution equation is set as a function of the cumulative dissipated energy during the loading history. Specific coefficients in the equation are determined by fitting and comparing the modulus decay curve measured in experiments with the predicted curve of the theoretical damage model. The fitting process also uses an optimization algorithm to find the coefficient values ​​that best match the predictions with the experimental data. Finally, all the parameters obtained from the above calibration are input into the software material model, thus completing the application of the viscoelastic damage constitutive model to the asphalt pavement.

[0064] A stress-dependent elastic damage constitutive model is applied to the base layer in the finite element numerical simulation model. This constitutive model describes the change in modulus of the semi-rigid base layer material with stress state and the stiffness degradation caused by cracking. In software implementation, this constitutive relationship is defined through a user-defined material subroutine. First, the baseline elastic modulus of the material is defined, i.e., the modulus value of the material under no-damage and small deformation conditions. The key steps are establishing the functional relationship between the modulus and stress state, and determining the threshold and law of damage initiation and evolution. The modulus-stress correlation is determined through triaxial compression tests of base layer materials, such as cement-stabilized crushed stone, under different confining pressures. Through a set of tests, the stress-strain relationship of the material under different constant confining pressures and different deviatoric stress levels is obtained, and the trend of elastic modulus variation with mean stress or spherical stress is analyzed. Based on these discrete data points, function fitting methods are used, such as fitting to a linear function or a power function, to determine the specific mathematical relationship between the modulus enhancement coefficient and the stress invariant. The introduction of damage requires determining the stress threshold for material damage and the law of damage evolution. The initial damage threshold of a material is its tensile strength threshold, determined through direct tensile or splitting tests. A monotonically increasing tensile or radial compressive stress is applied to a standard specimen until fracture, and the maximum tensile stress the specimen can withstand before failure is recorded. This maximum value is defined as the tensile strength threshold of the material, a material constant with a definite physical meaning and unit, such as megapascals (MPa). In the custom constitutive model, this measured tensile strength threshold is preset as a fixed input parameter. The parameters in the damage evolution equation are calibrated by analyzing the stress-strain softening curve of the material after exceeding the tensile strength threshold. For example, in a uniaxial tensile test, the stress decrease curve after the specimen reaches the peak stress (the tensile strength threshold) is recorded. Assuming a damage evolution equation in an exponential decay form, the characteristic parameters in the equation are determined through curve fitting. Finally, the baseline elastic modulus, modulus-stress correlation function, tensile strength threshold, and damage evolution equation parameters are all incorporated into the logic of the user material subroutine. This subroutine calculates the current tensile stress value at each material integration point in real time and compares it with a pre-set tensile strength threshold parameter. When the current tensile stress value is less than the tensile strength threshold, the material is considered to be in a linear elastic or stress-dependent elastic stage. When the current tensile stress value reaches or exceeds the tensile strength threshold, damage is determined to have begun to occur, and the increment of the damage variable is calculated according to the damage evolution equation, thereby realizing the stress-dependent elastic damage constitutive behavior. Associating this subroutine with the geometric components representing the base layer completes the model assignment.

[0065] A falling weight deflection load is applied to the constitutive finite element numerical simulation model to establish a finite element numerical simulation model capable of simulating the dynamic deflection time history response under different damage states. Specifically, a time-dependent surface pressure load boundary condition is defined in the software. This load simulates the impact of a falling weight deflectometer, and its time function is defined based on the measured impact force time history curve of the falling weight deflectometer, typically simplified to a half-sine waveform. The peak load magnitude is calculated based on the falling weight mass, drop height, and bearing plate area using the momentum theorem and pressure formula. The load duration is directly read from the measured load time history curve. The load is applied to a circular area on the road surface, the diameter of which is the same as the diameter of the bearing plate. In the model setup, boundary conditions need to be defined to simulate the support of the actual roadbed, fully constraining all degrees of freedom on the bottom surface of the model and constraining the normal displacements of the four sides of the model to zero. The analysis type is set to transient dynamic analysis, using an explicit time integration algorithm. The time step is automatically calculated and set according to the stability condition formula based on the minimum size of the model mesh and the sound velocity in the material. The total analysis time is set to three to five times the load duration to ensure sufficient response decay. After completing all settings and performing the solution calculation, the finite element numerical simulation model is established. By modifying the internal variables representing the damage state in the material constitutive model in subsequent steps, such as directly setting initial values ​​of damage variables greater than zero in the material parameters of the asphalt surface layer and base layer to represent pre-existing damage, this model can be used to simulate different preset damage states and output dynamic deflection time history response data at each sensor location.

[0066] S3. The dynamic response of the target asphalt pavement structure under different preset damage states is simulated using a finite element numerical simulation model, and the dynamic system attractor structure features corresponding to each preset damage state are extracted using the phase space reconstruction method. The specific implementation is as follows:

[0067] Using a finite element numerical simulation model, multiple preset damage states are set by adjusting the damage density parameter in the finite element numerical simulation model. Specifically, an existing finite element numerical simulation model file is opened, which is equipped with a viscoelastic damage constitutive model and a stress-dependent elastic damage constitutive model, and subjected to a falling weight deflection load. The damage density parameter is concretized in the model as an internal state variable in the material constitutive model. For asphalt surface layer, this parameter is a damage variable defined in the viscoelastic damage constitutive model that continuously varies from 0 to 1; for base course, this parameter is a damage variable defined in the stress-dependent elastic damage constitutive model, also varying from 0 to 1. To systematically simulate the changes from a healthy state to different degrees of damage, a set of representative combinations of damage density parameter values ​​needs to be pre-planned. For example, a set of incremental numerical sequences starting from 0 is set for the surface layer damage variable and the base course damage variable, respectively. The damage variable sequence for the surface layer can be set to 0, 0.1, 0.2, 0.3, 0.4, and 0.5, while the damage variable sequence for the base layer can be set to 0, 0.05, 0.15, 0.25, 0.35, and 0.45. By combining the damage variable values ​​of the surface layer and the base layer, multiple different preset damage states can be defined. For example, state one is where the surface layer damage variable is equal to 0 and the base layer damage variable is equal to 0; state two is where the surface layer damage variable is equal to 0.1 and the base layer damage variable is equal to 0.05, and so on. In the preprocessing module of the finite element numerical simulation software, by modifying the initial value parameters of the damage variables in the corresponding material model definition, each set of planned damage variable values ​​is assigned to the asphalt surface layer and base layer materials in the model. Each set of assignments is saved as an independent model calculation case, thus completing the setting of multiple preset damage states.

[0068] The finite element numerical simulation model is run under each preset damage state to obtain the corresponding dynamic deflection time history response data. Specifically, in the solver module of the finite element numerical simulation software, independent calculation cases established for each preset damage state are loaded sequentially. For each case, it is submitted to the solver for transient dynamic calculation. The solver calculates the dynamic response of the entire pavement structure model over time under a set falling weight deflection load, based on explicit or implicit time integration algorithms. After the calculation is completed, displacement time history data at predefined sensor locations are extracted from the results file. These sensor locations are consistent with the actual deflection basin measuring points detected by the falling weight deflectometer, including, for example, the load center point and points at distances of 200 mm, 300 mm, 450 mm, 600 mm, 900 mm, and 1200 mm from the center point. The extracted displacement data is the dynamic deflection time history response data corresponding to the preset damage state. It is usually stored in the form of a time series, with the length and step size of the time series consistent with the solution settings. For example, the total time is 100 milliseconds, and the output time interval is 0.1 milliseconds, thus obtaining a time history curve containing 1000 data points. This process is repeated for each preset damage state, eventually obtaining multiple sets of dynamic deflection time history response data corresponding to all preset damage states.

[0069] For each preset damage state, the dynamic deflection time history response data is used to reconstruct the phase space using the delayed coordinate method. The purpose of phase space reconstruction is to reconstruct the multidimensional state space of the dynamic system from a single variable time series. The specific implementation steps are as follows: First, a representative single variable time series is selected from a set of dynamic deflection time history response data as the basis for reconstruction. Usually, the dynamic deflection time history of the load center point is selected because the response at this point is the most significant. This time series is denoted as a discrete sequence of length N. Next, two key parameters need to be determined: the delay time and the embedding dimension. The delay time can be determined using the autocorrelation function method. The autocorrelation function of the time series is calculated, which describes the linear correlation between the sequence and itself at different time delays. An autocorrelation function decline threshold is set for judgment. This threshold is defined as the initial autocorrelation function value multiplied by a coefficient, which is 1 minus 1 divided by the natural constant e. The calculation result is approximately 0.37. This means that when the calculated autocorrelation function value drops to about 37% of the initial value, the corresponding delay time is selected as a suitable delay time. In practical calculations, the calculated autocorrelation function values ​​corresponding to different delay times are sequentially compared with a pre-calculated autocorrelation function descent threshold. The delay time corresponding to the first autocorrelation function value less than or equal to this threshold is selected. The embedding dimension can be determined using the spurious nearest neighbor method. The basic idea is that when the embedding dimension is large enough, the geometric structure of the trajectory in phase space can be fully unfolded, and the neighboring points are the true neighbors of the dynamic system, rather than spurious neighbors generated by projection into a low-dimensional space. Therefore, a threshold for the proportion of spurious nearest neighbors needs to be set to determine whether the embedding dimension is sufficient. This threshold is an empirical value, usually set between 2% and 10%. To balance computational accuracy and stability, it is recommended to set the spurious nearest neighbor proportion threshold to 5%. In practical applications, it can be fine-tuned around this recommended value, and the final determination can be made by observing the stability of the calculation results of invariants such as the correlation dimension under different thresholds. Starting with an embedding dimension of 1, the number is gradually increased. For each embedding dimension, the proportion of false nearest neighbors (FNN) is calculated and compared to a set threshold. When the calculated FNN proportion first falls below this threshold, the current embedding dimension is considered sufficient. After determining the delay time and embedding dimension, reconstruction can proceed. For each point in the time series, it is used as the starting point. Points with subsequent delay times are selected, and so on, until the required number of embedding dimensions are obtained. The values ​​of these points are arranged sequentially to form a multidimensional vector, which represents a state point in the reconstructed phase space. This operation is performed on all possible starting points in the time series, generating a set of trajectory points in an m-dimensional phase space. This set of points constitutes the geometric representation of the dynamical system attractor in the reconstructed phase space.For each preset damage state, the dynamic deflection time history response data is independently processed using the above parameter calculation and reconstruction process to obtain the corresponding reconstructed phase space trajectory.

[0070] In the reconstructed phase space, invariants characterizing the attractor's geometry and dynamics are calculated to obtain the attractor structural features corresponding to each preset damage state. The calculation of invariants specifically includes calculating the correlation dimension, the maximum Lyapunov exponent, and recursive quantization analysis indices. The correlation dimension is calculated to characterize the geometric complexity and fractal characteristics of the attractor. The correlation dimension is calculated based on the trajectory point set in the reconstructed phase space. First, a small distance parameter is defined. Then, the proportion of point pairs in the phase space whose distance is less than this small distance parameter is counted to the total number of point pairs; this proportion is called the correlation integral. The small distance parameter takes a series of values ​​from small to large, and a series of corresponding correlation integral values ​​are calculated. According to fractal theory, the correlation integral and the small distance parameter have a linear relationship in a log-log coordinate system, and its slope is the correlation dimension. By plotting a series of data points on a log-log plot and fitting the slope of the straight line using linear regression, the correlation dimension of the attractor under the preset damage state can be obtained. The maximum Lyapunov exponent is calculated to characterize the dynamic system's sensitivity to initial conditions, i.e., its chaotic characteristics. The calculation of the maximum Lyapunov exponent is based on tracking the evolution of the phase space trajectory. An initial point is selected on the reconstructed phase space trajectory, and its nearest neighbor is found, with the initial distance between the two points calculated. Then, the trajectory is evolved forward for a period of time, and the new distance between corresponding points on the two evolved trajectories is calculated. The logarithm of the ratio of the initial distance to the new distance, divided by the evolution time, yields a local Lyapunov exponent estimate. This process is repeated with multiple different initial points on the trajectory, and the multiple local estimates are averaged to obtain the final maximum Lyapunov exponent estimate for the entire attractor. A positive maximum Lyapunov exponent usually indicates that the system exhibits chaotic characteristics. Recursive quantization analysis is performed to further characterize the recursive properties of the attractor. First, a recursive graph needs to be constructed. The recursive graph is a two-dimensional square matrix, with its rows and columns corresponding to time points on the phase space trajectory. A recursive distance threshold needs to be set during construction to determine whether two phase space points are close enough to be considered recursive. The recursive distance threshold is typically set as a fixed proportion of the average distance between all point pairs in the entire phase space trajectory point set, for example, 10% of the average distance. In specific calculations, the average distance of all point pairs is first calculated, and then this average is multiplied by 0.1 to obtain the specific value of the recursive distance threshold. When constructing the recursive graph, for any two time points, the Euclidean distance between their corresponding points in phase space is calculated, and this calculated distance is compared with the recursive distance threshold. If the calculated distance is less than or equal to the recursive distance threshold, a point is marked on the coordinates of the recursive graph, indicating that state recursion has occurred. Based on the generated recursive graph, multiple recursive quantitative analysis indicators can be calculated. Among them, the recursion rate refers to the proportion of all marked points in the recursive graph to the total number of points, reflecting the overall probability of system state recursion.Determinism refers to the proportion of points on segments forming continuous diagonal structures in the recursive graph out of all recursive points, reflecting the degree of determinism in the system's behavior. Laminar flow refers to the average length of these continuous diagonal segments, reflecting the average time the system maintains a certain state. By calculating these indicators, the recursive pattern of the attractor can be quantified from different perspectives. Finally, the correlation dimension, the maximum Lyapunov exponent, and a set of recursive quantification analysis indicators, including recursion rate, determinism, and laminar flow, calculated under a preset damage state are packaged together as the structural features of the dynamical system attractor corresponding to that state. Repeating the above complete calculation process for all preset damage states yields a complete set of dynamical system attractor structural feature data corresponding one-to-one with all preset damage states.

[0071] Figure 2 A flowchart of the present invention for constructing a geometric constraint path characterizing the continuous variation law of damage density in phase space is given. S4, Based on the attractor structure characteristics of the dynamic system, a geometric constraint path characterizing the continuous variation law of damage density in phase space is constructed, specifically as follows:

[0072] Collect the dynamic system attractor structural features corresponding to all preset damage states, and associate each preset damage state and its corresponding damage density parameter value with the dynamic system attractor structural features. In practice, a spreadsheet is created for systematic data management. From the calculation results of previous steps, extract the complete dynamic system attractor structural feature data corresponding to each preset damage state. This data includes the correlation dimension, maximum Lyapunov exponent, recursion rate, determinism, and laminarity for each state. Simultaneously, record the damage density parameter value corresponding to each preset damage state. This value consists of a set of specific numerical values, including surface layer damage variable values ​​characterizing the degree of asphalt surface layer damage and base layer damage variable values ​​characterizing the degree of base layer damage. In the data table, each preset damage state is treated as an independent data record row, explicitly containing the surface layer damage variable value, base layer damage variable value, and its strictly corresponding correlation dimension value, maximum Lyapunov exponent value, recursion rate value, determinism value, and laminarity value. Through this tabular organization, a one-to-one and traceable data association is established between the identity of each preset damage state and its damage density parameter and dynamic system attractor structural features.

[0073] Using the structural features of the attractor in the dynamic system as coordinates, feature points corresponding to each preset damage state are marked in the phase space. Specifically, the coordinate axes used to construct the phase space are first defined. Several key invariants from the structural features of the attractor in the dynamic system are selected as the basis for the coordinate axes; for example, the correlation dimension is chosen as the first coordinate axis, the maximum Lyapunov exponent as the second coordinate axis, and the recurrence rate as the third coordinate axis, thus constructing a three-dimensional phase space. The units of the coordinate axes are consistent with the physical dimensions of the feature values; for example, the correlation dimension is dimensionless, the maximum Lyapunov exponent is in units of seconds, and the recurrence rate is a percentage dimensionless number. In the phase space with the defined coordinate axis system, for each preset damage state recorded in the data table, its corresponding set of feature values—the correlation dimension value, the maximum Lyapunov exponent value, and the recurrence rate value—is taken as a coordinate data point. This coordinate data point is then plotted or calculated and located within the defined phase space, and an identifier is used to indicate the preset damage state number and its damage density parameter value corresponding to that point. Repeat this operation for all preset damage states, and finally mark multiple feature points in the phase space that are equal in number to all preset damage states.

[0074] Curve fitting is used to fit and connect feature points in phase space corresponding to different damage density parameter values, generating a smooth and continuous path. The core of this step is to find a mathematically continuous function or curve that best describes the distribution trend of feature points in phase space. The specific implementation process is as follows: First, a single-valued, continuously varying comprehensive damage index needs to be calculated for each preset damage state to characterize its overall damage degree and serve as a parameter variable for fitting. The calculation of the comprehensive damage index requires setting a fusion weight for the surface layer damage variable and the base layer damage variable. This weight is determined based on the relative contribution ratio of the surface layer and base layer to the overall bearing capacity in pavement structural mechanics. For example, by consulting pavement design specifications or finite element analysis based on typical structures, it can be determined that under standard load, the asphalt surface layer bears approximately 60% of the load stress, while the base layer bears approximately 40%. Therefore, the weight of the surface layer damage variable is set to 0.6, and the weight of the base layer damage variable is set to 0.4. The comprehensive damage index equals the surface layer damage variable value multiplied by 0.6 plus the base layer damage variable value multiplied by 0.4. After calculating the comprehensive damage index corresponding to each preset damage state, the index is used as the independent variable, and the values ​​of the feature points on each coordinate axis are used as the dependent variables to perform parametric curve fitting.

[0075] The specific method for determining the fusion weights of surface layer and base course damage variables is as follows: First, using the finite element numerical simulation model of the target asphalt pavement structure in a non-destructive state established in step S2, mechanical response calculations are performed under the same load conditions as those detected by a falling weight deflectometer. After the calculation, the absolute values ​​of the maximum principal stresses (i.e., tensile stresses) of all units in the entire asphalt surface layer area under load are extracted, and their volume-weighted average is calculated, denoted as the surface layer characteristic stress. Similarly, the absolute values ​​of the vertical compressive stresses of all units in the entire base course area under load are extracted, and their volume-weighted average is calculated, denoted as the base course characteristic stress. Then, the sum of the surface layer characteristic stress and the base course characteristic stress is calculated. Finally, the proportion of the surface layer characteristic stress to the total is calculated separately, serving as the fusion weight of the surface layer damage variable, and the proportion of the base course characteristic stress to the total is calculated, serving as the fusion weight of the base course damage variable. Through this method, the weights directly reflect the relative proportions of the load stress contributions actually borne by the surface layer and base course under specific loads, ensuring that the comprehensive damage index can reasonably characterize the overall damage state of the specific structure.

[0076] For each coordinate axis, a polynomial curve fitting method is used. For example, for the correlation dimension axis, a cubic polynomial function is fitted: correlation dimension value = C0 + C1 × ZS + C2 × ZS 2 +C3×ZS 3 Where C0 is a constant term, C1 is the coefficient of the first term, C2 is the coefficient of the second term, C3 is the coefficient of the third term, and ZS is the comprehensive damage index. The values ​​of coefficients C0, C1, C2, and C3 are determined using the least squares method. The implementation process of the least squares method is as follows: a target function is constructed, which is equal to the sum of the squares of the differences between the actual correlation dimension values ​​of all feature points and the theoretical values ​​calculated according to the polynomial function. The values ​​of C0 to C3 are adjusted through a mathematical optimization algorithm to minimize the value of the target function. An iteration termination condition is set for the fitting process, that is, the change in the target function value is less than a preset fitting convergence judgment threshold, for example, 1 multiplied by 10 to the power of -6. When the difference between the target function values ​​of two consecutive iterations is less than this threshold, the fitting is considered to be converged, and the final coefficient values ​​are output. The same fitting process is repeated independently for the maximum Lyapunov exponent axis and the recurrence rate axis. Finally, polynomial functions of the three coordinate axes with respect to the comprehensive damage index are obtained. The combination of these three functions defines a smooth and continuous parametric path in three-dimensional phase space.

[0077] The generated smooth, continuous path is defined as a geometrically constrained path characterizing the continuous variation of damage density in phase space. Specifically, the coefficients of the three polynomial functions obtained from the fitting process, along with the weighting rules for calculating the comprehensive damage index, are stored in a standard mathematical description file or function library. The mathematical definition of this path explicitly specifies how the theoretical coordinates corresponding to a specific comprehensive damage index in phase space should be calculated. When this constraint needs to be applied, for any given surface layer damage variable value and base layer damage variable value, the comprehensive damage index is first calculated according to the stored weighting rules. Then, this index is substituted into the three stored polynomial functions to calculate the theoretical coordinates on the correlation dimension axis, the maximum Lyapunov exponent axis, and the recurrence rate axis, respectively. This set of theoretical coordinates defines the expected location point in phase space corresponding to the damage state and conforming to the continuous evolution law. This computable theoretical trajectory, explicitly defined by mathematical functions, is formally defined as a geometrically constrained path. It provides a rigorous mathematical reference system for subsequent steps. The distance between any dynamic system feature point derived from measured data and its corresponding point on this path can be used as a quantitative measure to evaluate whether the inversion results conform to the physical laws of continuous damage evolution.

[0078] S5. Based on geometrically constrained paths combined with deflection data, an optimization inversion algorithm with geometrically constrained paths as constraints is used to calculate the surface layer damage density and base layer damage density of the target asphalt pavement structure. The specific implementation is as follows:

[0079] Using the previously acquired deflection basin locations and characteristic values ​​as input deflection data, an objective function is constructed with surface layer damage density and base layer damage density as optimization variables. This objective function measures the difference between the predicted deflection data and the actual deflection data. Specifically, the deflection data input includes a series of deflection basin location coordinates obtained through a falling weight deflectometer, and corresponding characteristic values ​​for each location, such as the central deflection value, deflection basin area, and deflection basin diameter. Surface layer damage density and base layer damage density are the two core optimization variables to be inverted, with values ​​continuously varying between 0 and 1. The first step in constructing the objective function is to establish an inversion prediction model. This model is a mapping function whose inputs are arbitrary surface layer and base layer damage density values, and whose output is the predicted deflection basin characteristic values. This prediction model is implemented by calling a pre-established, parameterized finite element numerical simulation model. Given a set of surface layer damage density values ​​and base layer damage density values, modify the initial state of the corresponding material damage variables in the finite element model, run transient dynamic simulation, calculate the dynamic deflection time history at each measuring point with the same measured deflection basin location, and extract the predicted deflection basin characteristic values, such as the predicted center deflection value, the predicted deflection basin area, and the predicted deflection basin diameter. The specific form of the objective function is constructed as a weighted sum of squared residuals. For the i-th deflection basin location, calculate the difference between the measured j-th deflection basin characteristic value and the corresponding characteristic value predicted by the model; this difference is called the residual; where i represents the index of the deflection basin location, and j represents the index of the deflection basin characteristic type. Square the residuals of all characteristic values ​​at all deflection basin locations, multiply by a corresponding weighting coefficient, and finally sum all the weighted residual squares to obtain the objective function value. The weighting coefficient is set to balance the differences in magnitude and importance of different characteristic values. The specific setting method is as follows: First, calculate the average value of each deflection basin characteristic value at all measured deflection basin locations. Then, the initial weights of each feature are set to the reciprocal of their average value to normalize the magnitude. Afterward, the relative importance of each feature is fine-tuned based on engineering experience. Finally, all weights are normalized so that their sum equals 1. The smaller the value of the objective function, the smaller the overall difference between the inverted predicted deflection data and the actual deflection data.

[0080] The geometric constraint path is transformed into mathematical constraints on the optimization variables. This is achieved by defining a mapping between damage density and the attractor structural characteristics of the dynamic system in phase space. The geometric constraint path, constructed and stored as a mathematical function in previous steps, defines a continuous trajectory that a legal damage state point should follow in a three-dimensional phase space comprised of the correlation dimension, the maximum Lyapunov exponent, and the recurrence rate. This trajectory is described by three polynomial functions concerning the comprehensive damage index. To transform this into mathematical constraints on the two optimization variables, surface layer damage density and base layer damage density, it is necessary to establish their relationship with their positions in phase space. First, for any set of candidate surface layer and base layer damage density values, the comprehensive damage index is calculated according to the same rules used when constructing the path. Then, this comprehensive damage index is substituted into the stored three polynomial functions to calculate the corresponding theoretical coordinates in phase space, namely the theoretical correlation dimension value, the theoretical maximum Lyapunov exponent value, and the theoretical recurrence rate value. Next, based on the candidate surface and base layer damage density values, a fast surrogate model is used to calculate the corresponding actual dynamic system attractor structure features. This surrogate model is a pre-trained multivariate nonlinear regression model. Its inputs are the surface and base layer damage density values, and its outputs are the predicted correlation dimension, maximum Lyapunov exponent, and recurrence rate. The data used to train this surrogate model comes from the dataset of all preset damage states from previous steps. This dataset contains the surface and base layer damage variable values ​​for each preset state, along with the corresponding correlation dimension, maximum Lyapunov exponent, and recurrence rate. The training process uses the least squares method to fit the coefficients of the regression model, minimizing the mean square error between the model's predicted feature values ​​and the true feature values ​​in the dataset. After obtaining the actual predicted feature values ​​corresponding to the candidate points, the Euclidean distance between them and the theoretical coordinate points in phase space is calculated. This distance reflects the degree to which the candidate damage state deviates from the preset continuous evolution law. Therefore, the mathematical constraint is defined as the Euclidean distance must be less than or equal to a preset maximum allowable deviation threshold. This maximum allowable deviation threshold is set based on the residual level during the path fitting stage. The specific setting method is as follows: First, for each preset damage state used when constructing the path, calculate the vertical distance from its actual feature point in phase space to the fitted path. Then, calculate the average of all these distances. Finally, set the maximum allowable deviation threshold to twice this average distance. This constraint is enforced as an inequality constraint in the optimization algorithm, ensuring that the searched solution not only matches the deflection data in terms of mechanical response but also has physical rationality in terms of the evolution law of the dynamic system.

[0081] An artificial intelligence optimization algorithm is employed to iteratively update the surface layer damage density and base layer damage density parameters within a solution space that satisfies mathematical constraints, until the objective function converges to a preset convergence threshold, outputting the final calculated results for the surface layer damage density and base layer damage density. Specifically, a particle swarm optimization algorithm is chosen as an example of an artificial intelligence optimization algorithm. During algorithm initialization, a certain number of candidate solutions are randomly generated within the domain of surface layer damage density and base layer damage density (range 0 to 1). Each candidate solution is called a particle, and all particles constitute a population. The position of each particle is represented by a pair of values: a surface layer damage density value and a base layer damage density value. Each particle also has a velocity vector to determine its direction and distance of movement in the next iteration. Algorithm parameters need to be preset, including population size, maximum number of iterations, inertia weight, individual learning factor, and social learning factor. For example, the population size is set to 30, and the maximum number of iterations is set to 100. The inertia weight controls the tendency of particles to maintain their previous velocity; the initial value can be set to 0.9 and linearly decreased to 0.4 with each iteration. The individual learning factor and social learning factor are typically both set to 2.0. In each iteration, the following steps are performed for each particle. First, check whether the surface damage density value and the base layer damage density value corresponding to the particle's current position satisfy the aforementioned mathematical constraints, i.e., whether the calculated feature distance is less than or equal to the maximum allowable deviation threshold. If not, a penalty is imposed on the particle, for example, by setting its objective function value to a very large number, such as 1 multiplied by 10 to the power of 10, thereby reducing its probability of being selected as an excellent example in the population. If the constraints are satisfied, the objective function calculation process is invoked, i.e., the predicted deflection data corresponding to the set of damage density values ​​is calculated using the finite element prediction model, and the objective function value is calculated by comparing it with the measured deflection data. Each particle records the position where its historical objective function value is the smallest, called the individual historical best position. The entire population also records the position where the objective function value is the smallest among all particles, called the global historical best position. At the end of each iteration, based on the individual historical best position and the global historical best position, the velocity and position of each particle are updated according to the velocity and position update formula of the particle swarm optimization algorithm. The velocity update formula is: the new velocity equals the inertia weight multiplied by the old velocity, plus the individual learning factor multiplied by a random number between 0 and 1 multiplied by the difference between the individual's historical best position and the current position, plus the social learning factor multiplied by a random number between 0 and 1 multiplied by the difference between the global historical best position and the current position. The position update formula is: the new position equals the old position plus the new velocity. After the update, it is ensured that the new position value remains within the domain of 0 to 1 for both surface layer damage density and base layer damage density. The iterative process continues, and the change in the objective function value corresponding to the global historical best position is monitored. A preset convergence threshold is set as the stopping criterion, which is based on the expected accuracy of the objective function and the computational cost. For example, when the objective function value is less than 1 × 10⁻⁶, the threshold is lowered. -4At that time, it was considered that the matching degree between the predicted deflection and the actual deflection was already high enough, so the preset convergence threshold was set to 1×10. -4 The convergence criterion is that the algorithm is considered converged when the absolute value of the change in the global historical optimal objective function value is less than this preset convergence threshold for 10 consecutive iterations. At this point, the surface layer damage density and base layer damage density values ​​corresponding to the global historical optimal position are output as the final inversion calculation result. This result represents the quantified surface layer and base layer damage densities of the target asphalt pavement structure that best match the theoretical deflection response to the measured data, under the premise of satisfying the dynamic system evolution geometric constraints.

[0082] S6. Based on the surface layer damage density and base layer damage density, determine the internal damage state of the target asphalt pavement structure. The specific implementation is as follows:

[0083] The surface layer damage density and base course damage density at different locations, calculated in advance, are used as discrete data points. Specifically, from the output of the previous optimization inversion calculation step, the surface layer damage density values ​​and base course damage density values ​​corresponding to all calculated deflection basin locations are extracted. Each deflection basin location is defined by its planar coordinates on the target pavement, such as using geodetic coordinates or local relative coordinates. Therefore, for the surface layer damage density, a set of data is obtained, where each data point contains a planar location coordinate and the corresponding surface layer damage density value at that location, which is a dimensionless number between 0 and 1. Similarly, for the base course damage density, another set of data is obtained, where each data point contains the same planar location coordinate and the corresponding base course damage density value at that location, which is also a dimensionless number between 0 and 1. These two sets of data together constitute a spatially discrete set of damage density data points bound to specific measurement points.

[0084] Spatial interpolation is used to process discrete data points of surface layer and base course damage densities to generate continuous damage density distribution fields covering both the surface layer and base course of the target asphalt pavement structure. The purpose of spatial interpolation is to estimate the values ​​at any location within the entire evaluation area based on known values ​​at finite discrete points, thus forming a spatially continuously varying field. First, the interpolation computation domain needs to be defined, i.e., the area of ​​the target asphalt pavement structure to be evaluated. This area is typically determined by the outer boundary of all deflection basin measuring points. Then, regular grid nodes are generated within this computation domain, with the grid size set according to the evaluation accuracy requirements, for example, a grid size of 0.5 m × 0.5 m. Next, the spatial interpolation algorithm is applied to the surface layer damage density data point set and the base course damage density data point set, respectively. Ordinary Kriging interpolation is selected as an example. Kriging interpolation not only considers the distance between the estimated point and known points but also reflects the spatial structure characteristics of the data through a variogram model. The first step in implementing interpolation is to calculate the experimental variogram. For the surface damage density dataset, the square of the difference between the distance and damage density value between all pairs of data points is calculated. These pairs are then grouped by distance intervals, and half the average of the squared differences within each distance group is calculated to obtain the experimental variogram value for that distance group. Then, a theoretical variogram model is selected to fit these experimental values; commonly used models include the spherical model, exponential model, and Gaussian model. The fitting process determines three key parameters of the theoretical variogram model: nugget value, sill value, and range. The nugget value represents the microscale variation or measurement error; the sill value represents the total spatial variation of the variable; and the range represents the maximum distance of spatial autocorrelation, beyond which the data no longer have spatial correlation. The fitting process typically uses the least squares method, adjusting the parameters of the theoretical model to minimize the sum of squared differences between the variogram values ​​calculated by the theoretical model and the experimental values. After determining the variogram model and parameters, for each point to be estimated on the grid, its interpolation weights are calculated using the Kriging equations. The Kriging equations consist of two parts: a spatial covariance matrix between all known points, calculated based on the variogram; and a spatial covariance vector between all known points and the current point to be estimated. Solving this linear system of equations yields weighting coefficients assigned to each known point, which satisfy the conditions of unbiasedness and optimality. Finally, the damage density value of each known point is multiplied by its corresponding weighting coefficient and summed to obtain the interpolated estimate of the point to be estimated. This calculation is performed on all nodes of the grid to obtain the damage density estimates for the surface layer and base layer at each grid node across the entire evaluation area. Spatially connecting these gridded data forms two continuous and spatially correlated damage density distribution fields covering the surface layer and base layer of the target pavement structure. These distribution fields are stored and visualized in the form of two-dimensional matrices or raster data.

[0085] Based on the generated continuous damage density distribution field, regions in the surface layer and base layer where the damage density exceeds a preset damage density threshold, along with their corresponding damage density values, are identified. The preset damage density threshold is a critical value used to determine whether the pavement structure has entered a damage state requiring attention or intervention. This threshold is divided into a surface layer damage density threshold and a base layer damage density threshold. These two thresholds are not fixed but are determined comprehensively based on the pavement structure's design standards, material properties, maintenance history, and management strategies. Specific setting methods can be based on historical experience data or structural reliability theory. For example, one method based on historical data is to collect a large amount of damage density inversion data from similar pavements when maintenance is required, perform statistical analysis on this data, calculate its mean and standard deviation, and then set the preset damage density threshold as the value obtained by adding one standard deviation to the mean, using this as the boundary requiring early warning. Another method based on mechanical analysis is to establish the limit state equation of the pavement structure, calculate the critical damage density corresponding to functional failure of the pavement structure under standard load using a finite element model, and use this critical value as the preset damage density threshold. For example, the damage density threshold for the surface layer can be calculated when the tensile strain at the bottom of the asphalt surface layer reaches the material fatigue limit, and the damage density threshold for the base layer can be calculated when the compressive strain at the top surface reaches the allowable value. After obtaining the specific threshold values, region identification begins. The identification process is completed by traversing every grid cell in both the surface layer and base layer continuous damage density distribution fields. For the surface layer distribution field, the estimated damage density value stored at each grid cell is checked and compared with the surface layer damage density threshold. If the estimated damage density value of a grid cell is greater than or equal to the surface layer damage density threshold, the grid cell is marked as a surface layer exceeding the threshold. The planar coordinates of all surface layer exceeding units and their specific estimated damage density values ​​are recorded. Similarly, for the base layer distribution field, the estimated damage density value of each grid cell is compared with the base layer damage density threshold to identify and record the location and value of base layer exceeding units. In addition, spatial clustering analysis can be performed on the identified out-of-standard units to merge spatially adjacent out-of-standard units into a continuous out-of-standard region, and calculate the geometric characteristics of each region, such as area, as well as the statistical characteristics of damage density within the region, such as maximum and average values.

[0086] The internal damage state of the target asphalt pavement structure is determined by combining the areas exceeding a preset damage density threshold with their combined damage density values. This "combination" refers to integrating, analyzing, and expressing the identified surface and base course out-of-range information to form a complete diagnostic report on the pavement's internal damage condition. The specific steps are as follows: First, the sets of identified surface and base course out-of-range areas are overlaid to examine their spatial relationship. Then, for each identified out-of-range area, regardless of whether it belongs to the surface or base course, key attribute information is extracted, including the area's location boundary, area, average damage density within the area, and maximum damage density within the area. Next, the severity of the damage can be graded based on the degree to which the damage density value exceeds the threshold. For example, multiple grading threshold ranges can be set: areas with damage densities between the corresponding threshold and the threshold plus 0.1 are defined as lightly damaged areas; areas between the threshold plus 0.1 and the threshold plus 0.3 are defined as moderately damaged areas; and areas exceeding the threshold plus 0.3 are defined as severely damaged areas. Each out-of-range area is then rated for severity based on its average damage density. Finally, all this information is integrated into a structured data file or visualization report. This final output comprehensively quantifies and locates the position, extent, depth, and severity of damage within the target asphalt pavement structure; it is defined as the internal damage state of the pavement structure. This state is a multi-dimensional, quantifiable assessment result that can directly provide precise data support and spatial guidance for subsequent maintenance and repair decisions.

[0087] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0088] It should be noted that this invention can be deployed on the device itself to realize embedded applications, or it can run on a PC or other terminal with a user interface, thereby meeting various hardware environments and usage requirements.

[0089] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. A computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions according to the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. Computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wireless or wired transmission; wired transmission methods include optical fiber, twisted pair, coaxial cable, etc.; wireless transmission includes infrared, microwave, etc. Computer-readable storage media can be any available medium that a computer can access or a data storage device such as a server or data center that contains one or more sets of available media. Available media can be magnetic media (e.g., floppy disks, hard disks, magnetic tapes), optical media (e.g., DVDs), or semiconductor media. Semiconductor media can be solid-state drives.

[0090] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and modules described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0091] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0092] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical modules; they may be located in one place or distributed across multiple network modules. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0093] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0094] If a function is implemented as a software module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0095] The above are merely specific embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0096] In conclusion, the above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for constructing a damage model for existing asphalt pavements based on damage density, characterized in that, include: S1. Obtain the material and geometric parameters of the target asphalt pavement structure, as well as the deflection data of the target asphalt pavement structure under its current service condition. S2. Establish a finite element numerical simulation model of the target asphalt pavement structure based on material parameters and geometric parameters; S3. The dynamic response of the target asphalt pavement structure under different preset damage states is simulated using a finite element numerical simulation model. The phase space reconstruction method is then used to extract the attractor structural features of the dynamic system corresponding to each preset damage state, including: Using a finite element numerical simulation model, multiple preset damage states can be set by adjusting the damage density parameter in the finite element numerical simulation model. Run the finite element numerical simulation model under each preset damage state to obtain the corresponding dynamic deflection time history response data. For each preset damage state, the dynamic deflection time history response data is reconstructed using the delayed coordinate method; In the reconstructed phase space, invariants used to characterize the attractor's geometry and dynamic properties are calculated to obtain the dynamic system attractor structural features corresponding to each preset damage state. S4. Based on the attractor structure characteristics of the dynamic system, construct a geometric constraint path to characterize the continuous variation of damage density in phase space, including: Collect the structural features of the dynamic system attractor corresponding to all preset damage states, and associate each preset damage state and its corresponding damage density parameter value with the structural features of the dynamic system attractor. Using the attractor structural features of the dynamic system as coordinates, feature points corresponding to each preset damage state are marked in phase space; By using curve fitting, feature points in phase space corresponding to different damage density parameter values ​​are fitted and connected to generate smooth and continuous paths. The generated smooth and continuous path is defined as a geometrically constrained path that characterizes the continuous variation of damage density in phase space; S5. Based on geometric constraint paths combined with deflection data, the surface layer damage density and base layer damage density of the target asphalt pavement structure are calculated using an optimization inversion algorithm with geometric constraint paths as constraints. S6. Based on the surface layer damage density and base layer damage density, determine the internal damage state of the target asphalt pavement structure.

2. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 1, characterized in that, S1 includes: The dynamic modulus of the existing asphalt pavement surface layer and the resilient modulus of the base layer were obtained as material parameters through core sampling and indoor testing. Ground-penetrating radar was used to detect the thickness of each structural layer of the existing asphalt pavement and the bonding state between each layer as geometric parameters. The location and characteristic values ​​of the deflection basin of the existing asphalt pavement under its current service condition are obtained by using a falling weight deflectometer as deflection data.

3. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 1, characterized in that, S2 include: Input the material and geometric parameters into the multiphysics simulation platform; In the multiphysics simulation platform, a finite element numerical simulation model of the existing asphalt pavement structure is established based on material parameters and geometric parameters. A viscoelastic damage constitutive model is applied to the asphalt surface layer in the finite element numerical simulation model, and a stress-related elastic damage constitutive model is applied to the base layer. A falling weight deflection load was applied to the finite element numerical simulation model after constituting the model to establish a finite element numerical simulation model that can simulate the dynamic deflection time history response under different damage states.

4. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 1, characterized in that, The calculation of invariants used to characterize the attractor's geometry and dynamics includes: calculating the correlation dimension of the reconstructed phase space trajectory to characterize the geometry, and calculating the maximum Lyapunov exponent and recursive quantization analysis index to characterize the dynamics, thereby obtaining the attractor structural characteristics of the dynamic system.

5. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 1, characterized in that, Using the structural features of the attractor of the dynamic system as coordinates, feature points corresponding to each preset damage state are marked in the phase space. This includes: constructing a multidimensional phase space using multiple structural features of the attractor of the dynamic system, such as the correlation dimension and Lyapunov index, as coordinate axes, and using a set of specific feature values ​​corresponding to each preset damage state as coordinates for positioning and marking in the multidimensional phase space.

6. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 1, characterized in that, S5 include: The previously obtained deflection basin location and deflection basin characteristic values ​​are used as deflection data input, and an objective function is constructed with surface layer damage density and base layer damage density as optimization variables. The geometric constraint path is transformed into mathematical constraints on the optimization variables, where the geometric constraint path is achieved by defining the mapping relationship between damage density and the attractor structural characteristics of the dynamic system in phase space; An artificial intelligence optimization algorithm is used to iteratively update the surface layer damage density and base layer damage density parameters within the solution space that satisfies mathematical constraints, until the objective function converges to a preset convergence threshold, and the final calculation results of surface layer damage density and base layer damage density are output.

7. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 6, characterized in that, The objective function is used to measure the difference between the inverted predicted deflection data and the actual deflection data.

8. The method for constructing a damage model of existing asphalt pavement based on damage density according to claim 1, characterized in that, S6 include: The surface layer damage density and base layer damage density at different locations, calculated in advance, will be used as discrete data points; The discrete data points of surface layer damage density and base layer damage density are processed by spatial interpolation method to generate continuous damage density distribution fields covering the surface layer and base layer of the target asphalt pavement structure, respectively. Based on the generated continuous damage density distribution field, the regions in the surface layer and base layer whose damage density exceeds the preset damage density threshold and their damage density values ​​are identified. The area exceeding the preset damage density threshold and its combined damage density value are determined as the internal damage state of the target asphalt pavement structure.

Images