A method for predicting the fracture energy of fiber asphalt mixture under small sample machine learning

Abstract

The application discloses a kind of machine learning under small sample of fiber asphalt mixture fracture energy prediction method, the method is carried out data enhancement by WGAN-GP conditional generation confrontation network to alleviate small sample predicament, based on fracture mechanics theory to key interface parameter Physical monotone constraint is applied to improve the physical consistency of model, and combined with SHAP explainability analysis and conformal prediction uncertainty quantification, build the fracture energy intelligent prediction framework with precision, explainability and reliability, provide technical support for the crack resistance of fiber asphalt mixture Performance evaluation and material design.

Description

Technical Field

[0001] This invention relates to the prediction of road material properties, and in particular to a machine learning method for predicting the fracture energy of fiber-reinforced asphalt mixtures under small sample conditions. Background Technology

[0002] As a typical multiphase composite material, fiber-reinforced asphalt mixtures exhibit significant complexity and nonlinearity in their fracture behavior, influenced by a combination of factors including fiber characteristics, asphalt properties, aggregate lithology, and gradation composition. Fracture energy, a core indicator of crack resistance, can be measured using methods such as disc tensile testing (DCT) and semicircular bending tensile testing (SCB). However, laboratory tests suffer from limitations such as complex operation, long cycles, and high costs, making it difficult to comprehensively cover various material combinations. With limited samples, traditional machine learning methods struggle to fully learn the distribution patterns and nonlinear mapping relationships of the data, leading to overfitting during model training and compromised generalization performance. Purely data-driven models often operate as "black boxes," lacking effective verification mechanisms to ensure their predictive behavior strictly conforms to the fundamental mechanical mechanisms of fiber-reinforced composites. The expected monotonic relationship between parameters with clear physical meaning, such as interfacial shear strength (IFSS) and interfacial bond strength (IFBS), and fracture energy is difficult to guarantee. Furthermore, existing technologies often focus on improving accuracy, neglecting model interpretability and uncertainty quantification, failing to meet the engineering requirements of transparency and safety margins in material reverse engineering. Summary of the Invention

[0003] The purpose of this invention is to address the aforementioned problems by proposing a method for predicting the fracture energy of fiber-reinforced asphalt mixtures using a combination of generative adversarial networks (GAN-GP) and physical constraint machine learning (PCM) under small sample conditions. Data augmentation is achieved through WGAN-GP to alleviate the small sample difficulty. Physical monotonic constraints are applied to key interface parameters based on fracture mechanics theory to improve the physical consistency of the model. Furthermore, SHAP interpretability analysis and conformal prediction uncertainty quantification are combined to construct an intelligent fracture energy prediction framework that combines accuracy, interpretability, and reliability, providing technical support for the crack resistance assessment and material design of fiber-reinforced asphalt mixtures.

[0004] This invention is achieved through the following technical solution: a machine learning method for predicting the fracture energy of fiber-reinforced asphalt mixtures under small sample conditions, comprising the following steps:

[0005] P1: The obtained samples related to the fracture energy of fiber-reinforced asphalt mixtures are divided into training and test sets using a stratified random partitioning method. The key features affecting the distribution of working conditions are used as stratification variables to divide the mixture into several subsets of working conditions. Samples are randomly drawn from each subset according to a preset ratio and assigned to the training and test sets respectively. After merging, the training set and independent test set are obtained. The test set is locked throughout the process and does not participate in subsequent data generation, model training, and hyperparameter optimization.

[0006] P2: Data augmentation of the training set is performed using Wasserstein Generative Adversarial Network and Gradient Penalty (WGAN-GP); the generator input is noise and a conditional vector (material type encoding), and the output is data containing input and output features; the discriminator input is the concatenation of real / generated data and the conditional vector, and the output is a discrimination score; a gradient penalty mechanism is used to satisfy the 1-Lipschitz continuity constraint; the consistency of the distribution between generated samples and real samples is verified by t-SNE visualization.

[0007] P3: Using Bayesian Optimized Histogram Gradient Boosting Regression Tree (BO-HGBR) as the evaluation benchmark, ablation experiments were conducted with different generation ratios; model training and testing were repeated multiple times at each ratio, and the mean and standard deviation of the evaluation index of the test set were statistically analyzed to determine the optimal generation ratio and the corresponding augmented training set size.

[0008] P4: Based on the optimal augmented training set, repeated randomized experiments were conducted to systematically compare multiple candidate models, including linear baseline models, ensemble baseline models, BP neural networks, and histogram gradient boosting regression trees; each model was optimized using a differentiated hyperparameter search strategy, and the optimal baseline model was determined based on prediction accuracy and stability.

[0009] P5: Based on Spearman correlation analysis and fracture mechanics theory, physical monotonic constraints are applied to key features. The constraint vector is then passed into the optimal baseline model to ensure that the model's predicted response to the specified features strictly maintains the preset monotonicity throughout the entire input space, thus constructing a physical-guided machine learning model.

[0010] P6: Leave-one-out cross-validation was used to test the unbiasedness of the physical guidance model. A linear regression model of measured and predicted values ​​was established to test the significant differences between the intercept and 0 and the slope and 1. The inductive conformal prediction method was used to quantify the model prediction uncertainty. The absolute prediction residuals on the calibration set were used as the inconsistency score. Two-sided prediction intervals for independent test set samples were constructed for the set nominal confidence level to evaluate the actual coverage and the average interval width.

[0011] P7: The SHAP method is used to perform global and local interpretability analysis and calculate the contribution of each feature to the prediction results. The cumulative local effect (ALE) method is used to analyze the marginal effect of key features. By calculating the average change of the feature in the predicted value within the local interval, the marginal relationship between the feature and the prediction target is revealed unbiasedly, and the physical consistency of the model behavior is verified.

[0012] Furthermore, the P1 process is implemented through the following steps:

[0013] Step 1: Data Collection

[0014] Fracture energy data of fiber-reinforced asphalt mixtures were obtained through multi-scale tests, with a total of 48 valid samples collected. Seven parameters were selected as input features, with fracture energy (G) as the primary parameter. f J / m 2 The output features are as follows; the specific parameter statistical analysis is shown in Table 1 below;

[0015] Table 1 Statistical Analysis of Modeling Parameters

[0016]

[0017] Step 2: Determining Hierarchical Variables and Dividing Working Condition Subsets

[0018] To prevent imbalance in the distribution of working conditions between the training and test sets, the maximum nominal aggregate size (NMAS) and the effective temperature range of asphalt (UTI) were selected as stratified variables. NMAS values ​​included two categories: 13.2 mm and 19.0 mm, and UTI values ​​included two categories: 86℃ and 98℃. The 48 samples were divided into four subsets based on the combination of these two variables: Subset A (NMAS = 13.2 mm, UTI = 86℃), Subset B (NMAS = 13.2 mm, UTI = 98℃), Subset C (NMAS = 19.0 mm, UTI = 86℃), and Subset D (NMAS = 19.0 mm, UTI = 98℃). Each subset had 12 samples, with a uniform distribution.

[0019] Step 3: Combine the real training set and the independent test set

[0020] Within each working condition subset, samples are randomly selected in a ratio of 80%:20% and assigned to the training set and the test set respectively; that is, 10 sets of samples are randomly selected from each subset to enter the training set, and the remaining 2 sets of samples are entered into the test set.

[0021] The training samples extracted from the four subsets are merged to obtain 38 training sets; the test samples extracted from the four subsets are merged to obtain 10 independent test sets; the test sets are locked throughout the process and do not participate in subsequent data generation, model training and hyperparameter optimization.

[0022] Step 4: Validation of the partition

[0023] Verify the consistency of the distributions of the training and test sets on continuous features; calculate IFSS, IFBS, The mean and standard deviation of features such as OAC and E are shown in Table 2 below. The results show that the relative deviation of the mean between the training set and the test set does not exceed 10%, and the relative deviation of the standard deviation is within 15%, indicating that the two sets of data have good homogeneity in feature distribution. The stratified random partitioning effectively ensures the consistency of the working condition distribution between the training set and the test set.

[0024] Table 2. Statistical Comparison of Continuous Features between Training and Test Sets

[0025]

[0026] Note: Relative bias = |Training set mean - Test set mean| / Training set mean × 100%.

[0027] Furthermore, the P2 process is implemented through the following steps:

[0028] Step 1: Training Data Preparation and Conditional Vector Construction

[0029] The 38 training sets obtained in the P1 process were used as the augmentation baseline for WGAN-GP; for each set of samples, 7 input features (IFSS, IFBS, UTI, ...) were extracted. (OAC, E, NMAS) and one output characteristic (fracture energy G) f The data forms an 8-dimensional original data vector. Simultaneously, a condition vector is constructed based on the material properties of each sample: fiber type (4 types: non-fiber, basalt fiber, glass fiber, polyester fiber) is converted into a 4-dimensional vector using unique thermal encoding, asphalt type (2 types: No. 70 base asphalt, SBS modified asphalt) is converted into a 2-dimensional vector, and gradation type (2 types: AC-13, AC-20) is converted into a 2-dimensional vector. The three are then concatenated to obtain an 8-dimensional condition vector.

[0030] Step 2: Constructing the Generator and Discriminator Networks

[0031] The generator and discriminator network structure of WGAN-GP is constructed, and the specific configuration is shown in Table 3 below;

[0032] Table 3. Network Structure of WGAN-GP Generator and Discriminator

[0033]

[0034] Note: 1 Noise dimension 100 + condition vector dimension 8 = 108; 2 Data dimensions: 7 input features + 1 output fracture energy G f = 8-dimensional; 3 Data dimension 8 + condition vector dimension 8 = 16;

[0035] The generator learns the conditional distribution of real data by combining random noise and conditional vectors. The discriminator is used to distinguish between real samples and generated samples. The larger its output value, the more likely the input is to come from the real distribution.

[0036] Step 3: Gradient penalty mechanism and loss function design

[0037] A gradient penalty mechanism is used to replace the traditional WGAN weight pruning to satisfy the 1-Lipschitz continuity constraint; during each discriminator update, the real sample x and the generated sample x are compared. Perform linear interpolation to obtain interpolated samples. ,in Uniformly distributed random numbers, Calculate the gradient norm at the interpolated samples and construct the gradient penalty term. , where the gradient penalty coefficient λ=10;

[0038] The discriminator loss function is defined as:

[0039] (1)

[0040] in: The discriminator loss function value. To generate sample distribution Expectations For the discriminator to generate samples The output score, To determine the distribution of real data Expectations Let λ be the output score of the discriminator for the real sample x, and λ be the gradient penalty coefficient. For the discriminator output of interpolated samples gradient, The L2 norm of the gradient;

[0041] The generator loss function is defined as:

[0042] (2)

[0043] in: The generator loss function value. To generate sample distribution Expectations For the discriminator to generate samples The output score;

[0044] Step 4: Adversarial Training and Parameter Update

[0045] The Adam optimizer was used, with the learning rates for both the generator and discriminator set to 0.0002, Adam-β1=0.5, and β2=0.9. An alternating training strategy was adopted, with the generator updated once after every 5 updates to the discriminator. The batch size was set to 16, and the total number of training rounds was 30,000.

[0046] Step 5: Generate sample quality verification

[0047] The t-SNE algorithm is used to reduce the dimensionality of real samples and generated samples to three-dimensional space for visualization. The distribution of 500% scale samples generated at Epoch=3000, 10000, 20000 and 30000 is compared with that of real samples. A 95% confidence ellipsoid is drawn to represent the core region of the sample distribution.

[0048] Furthermore, the P3 process is implemented through the following steps:

[0049] Step 1: Determine the evaluation benchmark model

[0050] Bayesian optimized histogram gradient boosting regression tree (BO-HGBR) was adopted as the unified evaluation benchmark for ablation experiments. This model has been fully validated in the field of asphalt mixture mechanical property prediction, with high prediction accuracy and stability. It can eliminate the interference of model selection on the optimization results of the generation amount and ensure the consistency of evaluation standards.

[0051] Step 2: Construct augmented training sets with different proportions

[0052] Based on the WGAN-GP model trained in the P2 process, data augmentation was performed on 38 sets of real training data according to a preset generation ratio, constructing 6 augmented training sets of different sizes:

[0053] 0% ratio: Only 38 real training sets were used.

[0054] 100% ratio: 38 real samples + 38 generated samples, totaling 76 samples.

[0055] 200% ratio: 38 real samples + 76 generated samples, totaling 114 samples.

[0056] 300% ratio: 38 real samples + 114 generated samples, totaling 152 samples.

[0057] 500% ratio: 38 real samples + 190 generated samples, totaling 228 samples.

[0058] 800% ratio: 38 real samples + 304 generated samples, totaling 342 samples;

[0059] Step 3: Bayesian hyperparameter optimization

[0060] For each scaled augmented training set, a Bayesian hyperparameter search is performed using the Optuna framework, with the optimization objective being the root mean square error (RMSE) on the validation set. The hyperparameter search space is defined as follows:

[0061] Number of trees (n_estimators): 50~300, integer type

[0062] Maximum depth (max_depth): 3~12, integer type

[0063] Learning rate: 0.01–0.3, log-uniform distribution

[0064] Subsample ratio: 0.6–1.0, uniformly distributed.

[0065] Feature sampling ratio (max_features): 0.6–1.0, uniformly distributed.

[0066] Minimum number of samples per leaf node (min_samples_leaf): 2~20, integer type

[0067] Each hyperparameter optimization is evaluated through 50 iterations, and the optimal parameter combination is selected for subsequent training.

[0068] Step 4: Design and execution of repeated experiments

[0069] To evaluate the statistical stability of the model's performance, repeated experiments were conducted for each generation ratio, as follows:

[0070] The optimal combination of hyperparameters determined by Bayesian optimization with this fixed ratio

[0071] The model training and testing process was repeated 30 times, with a different random seed each time, to ensure that the partitioning of the training and testing sets was random.

[0072] Calculate the prediction performance metrics, including the coefficient of determination (R²), for each iteration on the independent test set (10 groups). 2 ), Root Mean Square Error (RMSE) and Mean Absolute Error (MAE)

[0073] The mean and standard deviation of 30 repetitions were recorded as the final evaluation result of the model's generalization performance at this generation ratio.

[0074] Step 5: Determine the optimal generation ratio

[0075] For each generation ratio, 30 repeated experiments were conducted, with different random seeds used for model training and testing in each experiment. The coefficient of determination (R²) on the test set was then calculated. 2 The mean and standard deviation of the root mean square error (RMSE) and mean absolute error (MAE) are calculated. By comparing the prediction accuracy and stability of the model at different scales, the impact of the generated sample size on the model performance is comprehensively evaluated.

[0076] Furthermore, the P4 process is implemented through the following steps:

[0077] Step 1: Candidate Model Selection

[0078] Based on the optimal augmented training set determined by the P3 process, six representative models were selected for systematic comparison, including the linear baseline model Ridge Regression, the ensemble baseline model Random Forest (RF), the BP neural network model, and the gradient boosting tree model (HGBR).

[0079] The objective function of the ridge regression model is:

[0080] (3)

[0081] in: For the regression coefficient vector, Let the target variable (measured fracture energy value) be a vector. For the input feature matrix, The sum of squared residuals measures the accuracy of the model fit. This is the regularization coefficient, which controls the balance between model complexity and fitting accuracy. This is an L2 regularization term that penalizes excessively large coefficients to prevent overfitting.

[0082] The output of the random forest prediction is:

[0083] (4)

[0084] in: These are the predicted values ​​from the random forest. For the number of decision trees, Let t be the prediction result of the t-th decision tree for input x;

[0085] The formula for calculating the hidden layer output of a BP neural network is:

[0086] (5)

[0087] in: This represents the output of the j-th neuron in the hidden layer. ( ) represents the hidden layer activation function (using the log-sigmoid function). The connection weights are the connection weights from the i-th neuron in the input layer to the j-th neuron in the hidden layer. Let n be the input value of the i-th neuron in the input layer, and n be the number of neurons in the input layer (which is 7). The threshold value is the value of the j-th neuron in the hidden layer.

[0088] The formula for calculating the output of the BP neural network output layer is:

[0089] (6)

[0090] in: This is the output (predicted value) of the k-th neuron in the output layer. ( ) represents the activation function of the output layer (using a linear function). The connection weights from the j-th neuron in the hidden layer to the k-th neuron in the output layer are... is the output of the j-th neuron in the hidden layer, and m is the number of neurons in the hidden layer (8). The threshold value for the k-th neuron in the output layer;

[0091] The formula for calculating the overall error of a BP network is:

[0092] (7)

[0093] Where: E is the overall network error, and D is the desired output vector. This is the actual output vector. This represents the expected value (measured value) of the k-th output neuron. This is the actual (predicted) value of the kth output neuron.

[0094] The formula for the log-sigmoid activation function is:

[0095] (8)

[0096] in: Let x be the output value of the activation function, x be the input value of the activation function (the difference between the weighted sum of neurons and the threshold), and e be the natural constant.

[0097] Histogram gradient boosting regression trees are based on the forward distribution algorithm, and their additive model can be expressed as:

[0098] (9)

[0099] in: For the HGBR model, the input Predicted value The total number of decision trees, Let be the weight (learning rate) of the m-th tree. This represents the prediction result of the m-th regression tree for the input x;

[0100] Step 2: Data partitioning and experimental setup

[0101] Thirty repeated random partitioning experiments were conducted to eliminate the influence of a single data partition on the results. In each experiment, the optimal augmented training set was randomly divided into a training subset and a validation subset in a ratio of 80%:20%. The training subset was used for model training, and the validation subset was used for performance evaluation during the hyperparameter optimization process.

[0102] Step 3: Hyperparameter Optimization Strategy

[0103] Different hyperparameter search strategies were designed to suit the characteristics of different models. For the linear model Ridge and the ensemble models RF and HGBR, Bayesian optimization (BO) was used to efficiently find the optimal parameters within a preset space. For the BP neural network, genetic algorithm (GA), Bayesian optimization (BO), and particle swarm optimization (PSO) were introduced to perform global hyperparameter search to fully explore its nonlinear fitting potential. The hyperparameter search space settings for each model are shown in Table 4 below.

[0104] Table 4. Model Hyperparameter Search Space

[0105]

[0106] Step 4: Comparison of model training results and determination of the optimal baseline model

[0107] By comparing the mean values ​​and box plot distributions of the indicators from 30 repeated experiments for each model, the predictive performance and stability of the six models were evaluated, and the optimal baseline model was finally determined. The formulas for calculating each evaluation indicator are as follows:

[0108] (10)

[0109] in: The coefficient of determination measures how well the model explains the variance of the target variable; its value ranges from [0,1]. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample. The average of the measured values. ;

[0110] (11)

[0111] in: The mean squared error is calculated by squaring the error, which amplifies the weight of larger errors. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample;

[0112] (12)

[0113] in: It is the root mean square error, monotonically related to the mean square error (MSE), and retains the same dimensions as the original response variable (J / m²). 2 ), For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample;

[0114] (13)

[0115] in: The mean absolute error directly reflects the average magnitude of the prediction error. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample.

[0116] Furthermore, the P5 process is implemented through the following steps:

[0117] Step 1: Feature Correlation Analysis and Determination of Constraints

[0118] Based on the optimal augmented training set determined by the P3 process, the Spearman rank correlation coefficient is used to calculate and quantify the relationship between each input feature and the fracture energy (G). f The degree of monotonic correlation between the above parameters and the fracture energy; combined with the fracture mechanics theory of fiber-reinforced composites: interfacial shear strength (IFSS) reflects the load transfer capacity of the fiber-mortar interface; the higher the value, the more energy the fiber consumes during pull-out; interfacial bond strength (IFBS) characterizes the adhesion performance between the mortar and the aggregate surface; strengthening the interfacial bond can delay the propagation of cracks along the mortar-aggregate interface; fiber modulus (E) characterizes fiber stiffness; high-modulus fibers can more effectively transfer stress in the mixture and limit crack opening; thus, it is determined that the above parameters and fracture energy should satisfy a positive monotonic relationship, that is, the fracture energy does not decrease monotonically as the parameter value increases;

[0119] Step 2: Design of Physical Monotonic Constraint Vectors

[0120] Based on the analysis in step 1, positive monotonic constraints are applied to the three core interface parameters IFSS, IFBS, and E. For the remaining features (UTI, G, OAC, NMAS), due to nonlinear interaction effects, threshold effects, or theoretical optimal values, global monotonic constraints are not suitable. The monotonic constraint parameter monotonic of the BO-HGBR model (the optimal baseline model derived from P4) is set. cstFor vectors:

[0121]

[0122] The vector elements correspond to IFSS, IFBS, UTI, G, OAC, E, and NMAS respectively. 1 indicates a positive monotonic constraint (the larger the eigenvalue, the larger the predicted value), and 0 indicates no constraint.

[0123] Step 3: Physics-guided model construction

[0124] The aforementioned constraint vectors are input into the BO-HGBR model to construct a physically guided PG-BO-HGBR model. During model training, the HGBR algorithm checks whether the splitting direction violates the preset monotonic constraint each time a decision tree node splits: for features with positive monotonic constraints, if the predicted mean of the left child node is greater than the predicted mean of the right child node after the current node splits, the split is prohibited, and the algorithm automatically selects other candidate split points. Through this mechanism, it is ensured that the predicted response of each base learner and the final ensemble model to the specified feature strictly maintains the preset monotonic direction throughout the entire input space. The same Bayesian optimization strategy as the unconstrained BO-HGBR is used to optimize the hyperparameters of PG-BO-HGBR; the search space is the same as in step 3 of P3.

[0125] Step 4: Model training results and generalization ability verification

[0126] Under optimal hyperparameter configuration, the PG-BO-HGBR model is trained using the optimal augmented training set to obtain the final physically guided fracture energy prediction model; this model can be expressed as:

[0127] (14)

[0128] Where the function To integrate decision tree models;

[0129] The predicted and measured values ​​of PG-BO-HGBR on the training and test sets were compared with the results of the unconstrained BO-HGBR model to verify the performance of the PG-BO-HGBR model. In order to evaluate the generalization ability of the physically guided PG-BO-HGBR model in real-world conditions, 10 independent test sets of samples that did not participate in any training were used for verification.

[0130] Furthermore, the P6 process is implemented through the following steps:

[0131] Step 1: Unbiased Regression Test

[0132] The optimal augmented training set determined by the P3 process is used as the data basis for the unbiasedness test; leave-one-out cross-validation is used to obtain the cross-validation prediction value for each sample; specifically, one sample is reserved as the validation set each time, and the remaining samples are used for training. This process is repeated a certain number of times with the optimal augmented training set samples to obtain the cross-validation prediction value for all samples. ;

[0133] Paired data obtained from leave-one-out cross-validation Establish measured values Compared with the predicted value Linear regression model between:

[0134] (15)

[0135] The intercept 'a' reflects the translation bias of the model (if a ≠ 0, the predicted value has a constant offset), and the slope 'b' reflects the scaling bias of the model (if b ≠ 1, the predicted value is magnified or reduced as a whole). An ideal unbiased model should satisfy a ≈ 0 and b ≈ 1.

[0136] The parameters a and b were estimated using ordinary least squares, and their significance was tested: the intercept a was tested to see if it was significantly non-zero (p<0.05 indicates translation bias); the slope b was tested to see if it was significantly non-equal to 1 (by constructing the t-statistic t=(b-1) / SE(b), p<0.05 indicates scaling bias).

[0137] Step 2: Local weighted regression smoothing analysis

[0138] To further reveal potential nonlinear patterns or heteroscedasticity in the residuals, a locally weighted regression scatter smoothing method was used to fit the trend of residual variation with predicted values; the predicted values ​​were then used as the basis for further analysis. The x-axis represents the residual. Plot a scatter plot with the vertical axis as the ordinate and overlay a LOWESS smoothing curve. If the smoothing curve shows no systematic trend (such as a U-shape, inverted U-shape, or monotonic trend) near the zero line, it indicates that the model residuals are uniformly distributed and there is no obvious heteroscedasticity.

[0139] Step 3: Conformal prediction quantifies model uncertainty

[0140] The inductive conformal prediction method is used to quantify the model's prediction uncertainty. The optimal enhancement training set determined by the P3 process is divided into a training subset and an independent calibration set according to the principle of hierarchical random partitioning. The training subset is used for training the PG-BO-HGBR model, and the calibration set is used to construct the non-consistent score distribution. The hierarchical variables still use NMAS and UTI to ensure that the operating condition distribution of the calibration set and the training subset is consistent. The sample size ratio of the training subset and the calibration set is set to 70%:30% based on engineering experience.

[0141] The absolute prediction residuals of the model on the calibration set are used as the inconsistency score to characterize the out-of-sample prediction error distribution of the model; for the nominal confidence level of 1-α commonly used in engineering (90%), the residual quantiles are calculated using the rigorous theoretical formula of conformal prediction.

[0142] (16)

[0143] in, Given the sample size of the independent calibration set, the quantile values ​​of the calibration residuals are calculated based on this quantile, serving as the half-width baseline for the prediction interval; for the independent test set samples, two-sided prediction intervals are constructed for the point prediction values ​​of the two models respectively:

[0144] (17)

[0145] in Let i be the model prediction value for the i-th sample. To calibrate the theoretical quantile values ​​of the residuals;

[0146] Two core evaluation indicators are selected: actual coverage (measures interval reliability, the closer to the nominal confidence level, the better) and average interval width (measures interval accuracy, the smaller the better for the same coverage). The calculation formula is as follows:

[0147] (18)

[0148] in, The actual coverage rate measures the reliability of the prediction interval, with a value between 0 and 1. denoted as the number of samples in the independent test set, and denoted as i as the index number of the test sample. This represents the measured (true) fracture energy of the i-th test sample. Let be the predicted fracture energy value for the i-th test sample. The quantile values ​​calculated based on the calibration set residuals are used as the half-width reference for the prediction interval. This is an indicator function; it takes the value 1 when the condition within the parentheses is true, and 0 otherwise. Let be the two-sided prediction interval for the i-th test sample;

[0149] (19)

[0150] in, The average interval width measures the accuracy of the prediction interval; its unit is the same as the fracture energy (J / m). 2 ), denoted as the number of samples in the independent test set, and denoted as i as the index number of the test sample. This represents the measured (true) fracture energy of the i-th test sample. Let be the predicted fracture energy value for the i-th test sample. The quantile values ​​calculated based on the calibration set residuals are used as the half-width reference for the prediction interval. This is the upper bound of the prediction interval. The lower bound of the prediction interval is given; plot the prediction interval comparison charts of the PG-BO-HGBR and unconstrained U-BO-HGBR models on the independent test set, as well as the coverage curves at different nominal confidence levels.

[0151] Furthermore, the P7 process is implemented through the following steps:

[0152] Step 1: SHAP Global Feature Importance Analysis

[0153] Based on the trained PG-BO-HGBR model, the SHAP method is used to calculate the contribution of each feature of each sample to the fracture energy prediction. The SHAP value is derived from the Shapley value in cooperative game theory. For model f, the SHAP value of the j-th feature of sample x is... The calculation formula is:

[0154] (20)

[0155] in, Let be the SHAP value of the j-th feature, and N be the set of all features. Let S be the total number of features, and S be the subset of features that does not contain feature j. Let f(S) be the size of the subset, and f(S) represent the predicted value based on the feature subset S. Based on the SHAP values ​​of all samples, the mean absolute SHAP value of each feature is calculated as a global importance indicator.

[0156] (twenty one)

[0157] in, Let be the global importance score of the j-th feature, and n be the total number of samples. Let be the SHAP value of the j-th feature of the i-th sample; sort each feature from high to low importance score, draw the SHAP dependency graph and the bee colony graph, and use the feature value as the x-axis and the SHAP value as the y-axis to reveal the influence of features on the prediction results;

[0158] Step 2: Verification of physical consistency of cumulative local effects (ALE)

[0159] The cumulative local effect method is used to analyze the marginal effect of key features on fracture energy prediction; for feature x j The range of values ​​is divided into K equal-frequency intervals. The local effect is calculated in the k-th interval, and the accumulated effect is the ALE value.

[0160] (twenty two)

[0161] in, Features In taking values The cumulative local effect value at the location, To determine the number of interval divisions, Features The value falls on the A sample set within a certain interval The number of samples in this set. For the first The upper boundary value of each interval. For the first The lower boundary value of each interval. For the i-th sample, excluding features Other feature values, The model prediction function is used; by comparing the ALE curves of the PG-BO-HGBR and unconstrained BO-HGBR models on the three features IFSS, IFBS, and E, it is verified whether the ALE curve strictly maintains positive monotonicity throughout the entire feature value range.

[0162] Step 3: Identification of key influencing factors and analysis of their mechanisms of action

[0163] By integrating the results of SHAP and ALE analyses, the contribution ranking of each feature to fracture energy prediction was clarified. The nonlinear action modes of each key variable were identified through the SHAP dependency graph, including the critical threshold of UTI, the monotonically increasing law of IFSS and IFBS, and the negative influence mechanism of NMAS. The marginal response of the model within the value range of each feature was verified by combining ALE curves to see if it conforms to the physical mechanism of fiber-reinforced composite materials, thus providing a theoretical basis for optimizing the crack resistance of fiber-reinforced asphalt mixtures.

[0164] The beneficial effects of this invention are:

[0165] (1) Effectively solves the small sample problem and significantly improves prediction accuracy and stability. This invention introduces WGAN-GP for data augmentation, and through ablation experiments, determines the optimal generation ratio of 500%, expanding the original 38 training sets to 228 sets. The augmented test set R... 2 The mean increased from 0.9493 to 0.9559, and the standard deviation of 30 replicates decreased from 0.00254 to 1.11 × 10⁻⁶. -16 This effectively suppresses the randomness of data partitioning and overcomes the technical defect of overfitting caused by small samples.

[0166] (2) Embedding physical monotonic constraints to ensure that the predicted behavior conforms to the fracture mechanics mechanism. This invention applies positive monotonic constraints to the three core interface parameters IFSS, IFBS, and E to construct a physically guided PG-BO-HGBR model. A comparison of ALE curves shows that the unconstrained model has local non-physical fluctuations, while PG-BO-HGBR strictly maintains a positive monotonic relationship throughout the entire input space, improving the extrapolation reliability of the model under unknown working conditions.

[0167] (3) Achieving unbiased estimation and uncertainty quantification enhances the reliability of engineering applications. Leave-one-out test results show that there are no significant differences between the intercept and 0, and between the slope and 1 (p>0.05), indicating that the model basically achieves unbiased estimation. The conformal prediction achieves an actual coverage rate of 80.0% at a 90% confidence level, and the prediction error of high-value samples is reduced by 1.45 percentage points, meeting the safety margin requirements of pavement material design.

[0168] (4) Improve the transparency of model decision-making and provide a basis for reverse design of materials. SHAP analysis reveals that UTI is the dominant influencing factor of fracture energy, followed by IFBS and IFSS, verifying the core role of interface reinforcement effect in crack-resistant systems. Characteristic dependency relationship identifies the nonlinear action mode of each key variable, providing a quantitative basis for optimizing the crack resistance of fiber-reinforced asphalt mixtures. Attached Figure Description

[0169] Figure 1 This is a flowchart of the method of the present invention;

[0170] Figure 2 T-SNE distribution diagrams of generated samples and real samples at different stages of this invention;

[0171] Figure 3 Box plot of the performance index of the present invention in 30 repeated experiments;

[0172] Figure 4 This is a scatter plot comparing the predicted and measured values ​​of the PG-BO-HGBR model of this invention.

[0173] Figure 5 This is a comparison chart of the prediction and measured results on the independent validation set of this invention;

[0174] Figure 6 This is a comparison of the ALE effect curves of the U-BO-HGBR and PG-BO-HGBR of the present invention in terms of constraint characteristics;

[0175] Figure 7 This is a graph showing the variation of the LOOCV residual with the predicted value in this invention;

[0176] Figure 8 This is a comparison chart of the conformal prediction interval performance of the present invention;

[0177] Figure 9 This is a distribution diagram of the SHAP values ​​of the features of this invention. Detailed Implementation

[0178] The present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0179] This invention, through systematic multi-scale experiments, determined the interfacial properties and fracture energy of asphalt mixtures with different fibers, obtaining a total of 48 valid samples. The test materials covered three types of fibers (basalt fiber, glass fiber, and polyester fiber), two types of asphalt (70 base asphalt and SBS modified asphalt), three types of aggregates (limestone, basalt, and granite), and two gradations (AC-13 and AC-20). All tests followed current standard testing methods to ensure the accuracy and comparability of the data. The testing methods for each modeling parameter are as follows:

[0180] Fiber-asphalt interfacial shear strength (IFSS, MPa): determined by a single fiber pull-out test. A monofilament fiber is embedded in the asphalt binder at a predetermined embedment length and pulled out at a constant rate using a universal testing machine. The load-displacement curve is recorded, and the interfacial shear strength is calculated based on the peak load and the fiber embedment area. This index reflects the load transfer capacity of the fiber-binder interface; a higher value indicates a stronger interfacial bond.

[0181] Fiber-modified bitumen-aggregate interfacial bond strength (IFBS, MPa): tested using a pull-out test (BBS). Bitumen mortar is applied to the surface of an aggregate substrate. After curing, a pull-out head is used to bond the mortar. A vertical pull-out force is applied at a constant rate until failure, and the peak load is recorded. The interfacial bond strength is calculated based on the pull-out area. This index characterizes the adhesion performance between the mortar and the aggregate surface and is a key indicator of crack resistance.

[0182] Effective Temperature Range of Asphalt (UTI, °C): Determined by the asphalt performance grading (PG grading), UTI = High-temperature continuous grading temperature − Low-temperature continuous grading temperature. For example, the UTI of PG76-22 asphalt = 76 − (−22) = 98 °C. This index reflects the working temperature range of asphalt; a larger value indicates stronger temperature adaptability.

[0183] Complex shear modulus of asphalt mortar (G * (MPa): Measured using a dynamic shear rheometer (DSR) at 60℃ and 10 rad / s. Before testing, the asphalt mortar was heated to a fluid state and then poured into specimens of specified dimensions. This index reflects the viscoelastic stiffness of the mortar and affects the mixture's ability to resist deformation.

[0184] Optimal asphalt-aggregate ratio (OAC, %): Determined using the Marshall test method according to the "Technical Specification for Construction of Asphalt Pavement on Highways" (JTG F40-2004). For AC-13 and AC-20 gradations, the initial asphalt-aggregate ratio was determined, followed by Marshall compaction tests. The optimal asphalt-aggregate ratio was determined comprehensively based on indicators such as stability, flow value, void ratio, and asphalt saturation. This indicator represents the optimal amount of asphalt and affects the interfacial film thickness and filling state.

[0185] Fiber modulus (E, GPa): Calculated according to ASTM D2256 standard, this involves tensile testing of the fiber monofilament. The monofilament is clamped at both ends in a universal testing machine fixture and stretched at a constant rate until fracture. The tensile modulus is calculated based on the initial linear segment of the stress-strain curve. This index characterizes fiber stiffness and is an important parameter for fiber reinforcement.

[0186] Maximum nominal aggregate size (NMAS, mm): Determined based on gradation type. The maximum nominal aggregate size for AC-13 gradation is 13.2 mm, and for AC-20 gradation it is 19.0 mm. This indicator reflects the fineness of the aggregate and affects the internal structure of the mixture and the crack propagation path.

[0187] fracture energy (G) f J / m 2 The fracture energy was determined according to ASTM D7313 using a disc tensile test (DCT). Rotationally compacted specimens were cut into discs of specified dimensions and held at 25°C for at least 4 hours. Loading was then applied using a UTM-25 material testing machine in strain-controlled mode, with a compressive displacement at the opening (CMOD) rate set to 1 mm / min. The load-compressive displacement curve was recorded until the load decreased to 0.1 kN, at which point the test was terminated. The area envelope of the curve was calculated and divided by the product of the specimen thickness and the initial ligament length to obtain the fracture energy. This index comprehensively characterizes the crack resistance of the mixture; a higher value indicates stronger crack resistance.

[0188] The above seven input features cover interfacial properties, asphalt characteristics, mix proportions, fiber properties, and gradation, comprehensively describing the main influencing factors of fracture energy; the single output feature is the fracture energy of fiber-reinforced asphalt mixtures. These 48 sets of sample data provide a reliable foundation for subsequent data enhancement and modeling.

[0189] A machine learning method for predicting the fracture energy of fiber-reinforced asphalt mixtures with small sample sizes, such as... Figure 1 As shown, it includes the following steps:

[0190] P1: The obtained samples related to the fracture energy of fiber-reinforced asphalt mixtures are divided into training and test sets using a stratified random partitioning method. The key features affecting the distribution of working conditions are used as stratification variables to divide the mixture into several subsets of working conditions. Samples are randomly drawn from each subset according to a preset ratio and assigned to the training and test sets respectively. After merging, the training set and independent test set are obtained. The test set is locked throughout the process and does not participate in subsequent data generation, model training, and hyperparameter optimization.

[0191] The P1 process is implemented through the following steps:

[0192] Step 1: Data Collection

[0193] Fracture energy data of fiber-reinforced asphalt mixtures were obtained through multi-scale tests, with a total of 48 valid samples collected. Seven parameters were selected as input features, with fracture energy (G) as the primary parameter. f J / m 2 The output features are as follows; the specific parameter statistical analysis is shown in Table 1 below;

[0194] Table 1 Statistical Analysis of Modeling Parameters

[0195]

[0196] Step 2: Determining Hierarchical Variables and Dividing Working Condition Subsets

[0197] To prevent imbalance in the distribution of working conditions between the training and test sets, the maximum nominal aggregate size (NMAS) and the effective temperature range of asphalt (UTI) were selected as stratified variables. NMAS values ​​included two categories: 13.2 mm and 19.0 mm, and UTI values ​​included two categories: 86℃ and 98℃. The 48 samples were divided into four subsets based on the combination of these two variables: Subset A (NMAS = 13.2 mm, UTI = 86℃), Subset B (NMAS = 13.2 mm, UTI = 98℃), Subset C (NMAS = 19.0 mm, UTI = 86℃), and Subset D (NMAS = 19.0 mm, UTI = 98℃). Each subset had 12 samples, with a uniform distribution.

[0198] Step 3: Combine the real training set and the independent test set

[0199] Within each working condition subset, samples are randomly selected in a ratio of 80%:20% and assigned to the training set and the test set respectively; that is, 10 sets of samples are randomly selected from each subset to enter the training set, and the remaining 2 sets of samples are entered into the test set.

[0200] The training samples extracted from the four subsets are merged to obtain 38 training sets; the test samples extracted from the four subsets are merged to obtain 10 independent test sets; the test sets are locked throughout the process and do not participate in subsequent data generation, model training and hyperparameter optimization.

[0201] Step 4: Validation of the partition

[0202] Verify the consistency of the distributions of the training and test sets on continuous features; calculate IFSS, IFBS, The mean and standard deviation of features such as OAC and E are shown in Table 2 below. The results show that the relative deviation of the mean between the training set and the test set does not exceed 10%, and the relative deviation of the standard deviation is within 15%, indicating that the two sets of data have good homogeneity in feature distribution. The stratified random partitioning effectively ensures the consistency of the working condition distribution between the training set and the test set.

[0203] Table 2. Statistical Comparison of Continuous Features between Training and Test Sets

[0204]

[0205] Note: Relative bias = |Training set mean - Test set mean| / Training set mean × 100%.

[0206] P2: Data augmentation of the training set is performed using Wasserstein Generative Adversarial Network and Gradient Penalty (WGAN-GP); the generator input is noise and a conditional vector (material type encoding), and the output is data containing input and output features; the discriminator input is the concatenation of real / generated data and the conditional vector, and the output is a discrimination score; a gradient penalty mechanism is used to satisfy the 1-Lipschitz continuity constraint; the consistency of the distribution between generated samples and real samples is verified by t-SNE visualization.

[0207] The P2 process is implemented through the following steps:

[0208] Step 1: Training Data Preparation and Conditional Vector Construction

[0209] The 38 training sets obtained in the P1 process were used as the augmentation baseline for WGAN-GP; for each set of samples, 7 input features (IFSS, IFBS, UTI, ...) were extracted. (OAC, E, NMAS) and one output characteristic (fracture energy G) f The data forms an 8-dimensional original data vector. Simultaneously, a condition vector is constructed based on the material properties of each sample: fiber type (4 types: non-fiber, basalt fiber, glass fiber, polyester fiber) is converted into a 4-dimensional vector using unique thermal encoding, asphalt type (2 types: No. 70 base asphalt, SBS modified asphalt) is converted into a 2-dimensional vector, and gradation type (2 types: AC-13, AC-20) is converted into a 2-dimensional vector. The three are then concatenated to obtain an 8-dimensional condition vector.

[0210] Step 2: Constructing the Generator and Discriminator Networks

[0211] The generator and discriminator network structure of WGAN-GP is constructed, and the specific configuration is shown in Table 3 below;

[0212] Table 3. Network Structure of WGAN-GP Generator and Discriminator

[0213]

[0214] Note: 1 Noise dimension 100 + condition vector dimension 8 = 108; 2 Data dimensions: 7 input features + 1 output fracture energy G f = 8-dimensional; 3 Data dimension 8 + condition vector dimension 8 = 16;

[0215] The generator learns the conditional distribution of real data by combining random noise and conditional vectors. The discriminator is used to distinguish between real samples and generated samples. The larger its output value, the more likely the input is to come from the real distribution.

[0216] Step 3: Gradient penalty mechanism and loss function design

[0217] A gradient penalty mechanism is used to replace the traditional WGAN weight pruning to satisfy the 1-Lipschitz continuity constraint; during each discriminator update, the real sample x and the generated sample x are compared. Perform linear interpolation to obtain interpolated samples. ,in Uniformly distributed random numbers, Calculate the gradient norm at the interpolated samples and construct the gradient penalty term. , where the gradient penalty coefficient λ=10;

[0218] The discriminator loss function is defined as:

[0219] (1)

[0220] in: The discriminator loss function value. To generate sample distribution Expectations For the discriminator to generate samples The output score, To determine the distribution of real data Expectations Let λ be the output score of the discriminator for the real sample x, and λ be the gradient penalty coefficient. For the discriminator output of interpolated samples gradient, The L2 norm of the gradient;

[0221] The generator loss function is defined as:

[0222] (2)

[0223] in: The generator loss function value. To generate sample distribution Expectations For the discriminator to generate samples The output score;

[0224] Step 4: Adversarial Training and Parameter Update

[0225] The Adam optimizer was used, with the learning rates for both the generator and discriminator set to 0.0002, Adam-β1=0.5, and β2=0.9. An alternating training strategy was adopted, with the generator updated once after every 5 updates to the discriminator. The batch size was set to 16, and the total number of training rounds was 30,000.

[0226] Step 5: Generate sample quality verification

[0227] The t-SNE algorithm is used to reduce the dimensionality of real and generated samples to three-dimensional space for visualization. Samples generated at Epochs of 3000, 10000, 20000, and 30000 are compared with real samples to illustrate their distribution. A 95% confidence ellipsoid is plotted to represent the core region of the sample distribution. Figure 2 As shown.

[0228] P3: Using Bayesian Optimized Histogram Gradient Boosting Regression Tree (BO-HGBR) as the evaluation benchmark, ablation experiments were conducted with different generation ratios; model training and testing were repeated multiple times at each ratio, and the mean and standard deviation of the evaluation index of the test set were statistically analyzed to determine the optimal generation ratio and the corresponding augmented training set size.

[0229] The P3 process is implemented through the following steps:

[0230] Step 1: Determine the evaluation benchmark model

[0231] Bayesian optimized histogram gradient boosting regression tree (BO-HGBR) was adopted as the unified evaluation benchmark for ablation experiments. This model has been fully validated in the field of asphalt mixture mechanical property prediction, with high prediction accuracy and stability. It can eliminate the interference of model selection on the optimization results of the generation amount and ensure the consistency of evaluation standards.

[0232] Step 2: Construct augmented training sets with different proportions

[0233] Based on the WGAN-GP model trained in the P2 process, data augmentation was performed on 38 sets of real training data according to a preset generation ratio, constructing 6 augmented training sets of different sizes:

[0234] 0% ratio: Only 38 real training sets were used.

[0235] 100% ratio: 38 real samples + 38 generated samples, totaling 76 samples.

[0236] 200% ratio: 38 real samples + 76 generated samples, totaling 114 samples.

[0237] 300% ratio: 38 real samples + 114 generated samples, totaling 152 samples.

[0238] 500% ratio: 38 real samples + 190 generated samples, totaling 228 samples.

[0239] 800% ratio: 38 real samples + 304 generated samples, totaling 342 samples;

[0240] Step 3: Bayesian hyperparameter optimization

[0241] For each scaled augmented training set, a Bayesian hyperparameter search is performed using the Optuna framework, with the optimization objective being the root mean square error (RMSE) on the validation set. The hyperparameter search space is defined as follows:

[0242] Number of trees (n_estimators): 50~300, integer type

[0243] Maximum depth (max_depth): 3~12, integer type

[0244] Learning rate: 0.01–0.3, log-uniform distribution

[0245] Subsample ratio: 0.6–1.0, uniformly distributed.

[0246] Feature sampling ratio (max_features): 0.6–1.0, uniformly distributed.

[0247] Minimum number of samples per leaf node (min_samples_leaf): 2~20, integer type

[0248] Each hyperparameter optimization is evaluated through 50 iterations, and the optimal parameter combination is selected for subsequent training.

[0249] Step 4: Design and execution of repeated experiments

[0250] To evaluate the statistical stability of the model's performance, repeated experiments were conducted for each generation ratio, as follows:

[0251] The optimal combination of hyperparameters determined by Bayesian optimization with this fixed ratio

[0252] The model training and testing process was repeated 30 times, with a different random seed each time, to ensure that the partitioning of the training and testing sets was random.

[0253] Calculate the prediction performance metrics, including the coefficient of determination (R²), for each iteration on the independent test set (10 groups). 2 ), Root Mean Square Error (RMSE) and Mean Absolute Error (MAE)

[0254] The mean and standard deviation of 30 repetitions were recorded as the final evaluation result of the model's generalization performance at this generation ratio.

[0255] Step 5: Determine the optimal generation ratio

[0256] For each generation ratio, 30 repeated experiments were conducted, with different random seeds used for model training and testing in each experiment. The coefficient of determination (R²) on the test set was then calculated. 2 The mean and standard deviation of the root mean square error (RMSE) and mean absolute error (MAE) are calculated. By comparing the prediction accuracy and stability of the model at different scales, the impact of the generated sample size on the model performance is comprehensively evaluated.

[0257] P4: Based on the optimal augmented training set, repeated randomized experiments were conducted to systematically compare multiple candidate models, including linear baseline models, ensemble baseline models, BP neural networks, and histogram gradient boosting regression trees; each model was optimized using a differentiated hyperparameter search strategy, and the optimal baseline model was determined based on prediction accuracy and stability.

[0258] The P4 process is implemented through the following steps:

[0259] Step 1: Candidate Model Selection

[0260] Based on the optimal augmented training set determined by the P3 process, six representative models were selected for systematic comparison, including the linear baseline model Ridge Regression, the ensemble baseline model Random Forest (RF), the BP neural network model, and the gradient boosting tree model (HGBR).

[0261] The objective function of the ridge regression model is:

[0262] (3)

[0263] in: For the regression coefficient vector, Let the target variable (measured fracture energy value) be a vector. For the input feature matrix, The sum of squared residuals measures the accuracy of the model fit. This is the regularization coefficient, which controls the balance between model complexity and fitting accuracy. This is an L2 regularization term that penalizes excessively large coefficients to prevent overfitting.

[0264] The output of the random forest prediction is:

[0265] (4)

[0266] in: These are the predicted values ​​from the random forest. For the number of decision trees, Let t be the prediction result of the t-th decision tree for input x;

[0267] The formula for calculating the hidden layer output of a BP neural network is:

[0268] (5)

[0269] in: This represents the output of the j-th neuron in the hidden layer. ( ) represents the hidden layer activation function (using the log-sigmoid function). The connection weights are the connection weights from the i-th neuron in the input layer to the j-th neuron in the hidden layer. Let n be the input value of the i-th neuron in the input layer, and n be the number of neurons in the input layer (which is 7). The threshold value is the value of the j-th neuron in the hidden layer.

[0270] The formula for calculating the output of the BP neural network output layer is:

[0271] (6)

[0272] in: This is the output (predicted value) of the k-th neuron in the output layer. ( ) represents the activation function of the output layer (using a linear function). The connection weights from the j-th neuron in the hidden layer to the k-th neuron in the output layer are... is the output of the j-th neuron in the hidden layer, and m is the number of neurons in the hidden layer (8). The threshold value for the k-th neuron in the output layer;

[0273] The formula for calculating the overall error of a BP network is:

[0274] (7)

[0275] Where: E is the overall network error, and D is the desired output vector. This is the actual output vector. This represents the expected value (measured value) of the k-th output neuron. This is the actual (predicted) value of the kth output neuron.

[0276] The formula for the log-sigmoid activation function is:

[0277] (8)

[0278] in: Let x be the output value of the activation function, x be the input value of the activation function (the difference between the weighted sum of neurons and the threshold), and e be the natural constant.

[0279] Histogram gradient boosting regression trees are based on the forward distribution algorithm, and their additive model can be expressed as:

[0280] (9)

[0281] in: For the HGBR model, the input Predicted value The total number of decision trees, Let be the weight (learning rate) of the m-th tree. This represents the prediction result of the m-th regression tree for the input x;

[0282] Step 2: Data partitioning and experimental setup

[0283] Thirty repeated random partitioning experiments were conducted to eliminate the influence of a single data partition on the results. In each experiment, the optimal augmented training set was randomly divided into a training subset and a validation subset in a ratio of 80%:20%. The training subset was used for model training, and the validation subset was used for performance evaluation during the hyperparameter optimization process.

[0284] Step 3: Hyperparameter Optimization Strategy

[0285] Different hyperparameter search strategies were designed to suit the characteristics of different models. For the linear model Ridge and the ensemble models RF and HGBR, Bayesian optimization (BO) was used to efficiently find the optimal parameters within a preset space. For the BP neural network, genetic algorithm (GA), Bayesian optimization (BO), and particle swarm optimization (PSO) were introduced to perform global hyperparameter search to fully explore its nonlinear fitting potential. The hyperparameter search space settings for each model are shown in Table 4 below.

[0286] Table 4. Model Hyperparameter Search Space

[0287]

[0288] Step 4: Comparison of model training results and determination of the optimal baseline model

[0289] By comparing the mean values ​​of the indicators from 30 repeated experiments of each model with the box plot distribution (e.g.) Figure 3 As shown in the figure, the predictive performance and stability of the six models were evaluated, and the optimal baseline model was finally determined. The calculation formulas for each evaluation index are as follows:

[0290] (10)

[0291] in: The coefficient of determination measures how well the model explains the variance of the target variable; its value ranges from [0,1]. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample. The average of the measured values. ;

[0292] (11)

[0293] in: The mean squared error is calculated by squaring the error, which amplifies the weight of larger errors. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample;

[0294] (12)

[0295] in: It is the root mean square error, monotonically related to the mean square error (MSE), and retains the same dimensions as the original response variable (J / m²). 2 ), For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample;

[0296] (13)

[0297] in: The mean absolute error directly reflects the average magnitude of the prediction error. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample.

[0298] P5: Based on Spearman correlation analysis and fracture mechanics theory, physical monotonic constraints are applied to key features. The constraint vector is then passed into the optimal baseline model to ensure that the model's predicted response to the specified features strictly maintains the preset monotonicity throughout the entire input space, thus constructing a physical-guided machine learning model.

[0299] The P5 process is implemented through the following steps:

[0300] Step 1: Feature Correlation Analysis and Determination of Constraints

[0301] Based on the optimal augmented training set determined by the P3 process, the Spearman rank correlation coefficient is used to calculate and quantify the relationship between each input feature and the fracture energy (G). f The degree of monotonic correlation between the above parameters and the fracture energy; combined with the fracture mechanics theory of fiber-reinforced composites: interfacial shear strength (IFSS) reflects the load transfer capacity of the fiber-mortar interface; the higher the value, the more energy the fiber consumes during pull-out; interfacial bond strength (IFBS) characterizes the adhesion performance between the mortar and the aggregate surface; strengthening the interfacial bond can delay the propagation of cracks along the mortar-aggregate interface; fiber modulus (E) characterizes fiber stiffness; high-modulus fibers can more effectively transfer stress in the mixture and limit crack opening; thus, it is determined that the above parameters and fracture energy should satisfy a positive monotonic relationship, that is, the fracture energy does not decrease monotonically as the parameter value increases;

[0302] Step 2: Design of Physical Monotonic Constraint Vectors

[0303] Based on the analysis in step 1, positive monotonic constraints are applied to the three core interface parameters IFSS, IFBS, and E. For the remaining features (UTI, G, OAC, NMAS), due to nonlinear interaction effects, threshold effects, or theoretical optimal values, global monotonic constraints are not suitable. The monotonic constraint parameter monotonic of the BO-HGBR model (the optimal baseline model derived from P4) is set. cst For vectors:

[0304]

[0305] The vector elements correspond to IFSS, IFBS, UTI, G, OAC, E, and NMAS respectively. 1 indicates a positive monotonic constraint (the larger the eigenvalue, the larger the predicted value), and 0 indicates no constraint.

[0306] Step 3: Physics-guided model construction

[0307] The aforementioned constraint vectors are input into the BO-HGBR model to construct a physically guided PG-BO-HGBR model. During model training, the HGBR algorithm checks whether the splitting direction violates the preset monotonic constraint each time a decision tree node splits: for features with positive monotonic constraints, if the predicted mean of the left child node is greater than the predicted mean of the right child node after the current node splits, the split is prohibited, and the algorithm automatically selects other candidate split points. Through this mechanism, it is ensured that the predicted response of each base learner and the final ensemble model to the specified feature strictly maintains the preset monotonic direction throughout the entire input space. The same Bayesian optimization strategy as the unconstrained BO-HGBR is used to optimize the hyperparameters of PG-BO-HGBR; the search space is the same as in step 3 of P3.

[0308] Step 4: Model training results and generalization ability verification

[0309] Under optimal hyperparameter configuration, the PG-BO-HGBR model is trained using the optimal augmented training set to obtain the final physically guided fracture energy prediction model; this model can be expressed as:

[0310] (14)

[0311] Where the function To integrate decision tree models;

[0312] The performance of the PG-BO-HGBR model is verified by comparing the predicted and measured values ​​of PG-BO-HGBR on the training and test sets with the results of the unconstrained BO-HGBR model. Figure 4 As shown, and to evaluate the generalization ability of the physically guided PG-BO-HGBR model under real-world conditions, a test set of 10 independent samples that did not participate in any training was used for validation. Figure 5 As shown.

[0313] P6: Leave-one-out cross-validation was used to test the unbiasedness of the physical guidance model. A linear regression model of measured and predicted values ​​was established to test the significant differences between the intercept and 0 and the slope and 1. The inductive conformal prediction method was used to quantify the model prediction uncertainty. The absolute prediction residuals on the calibration set were used as the inconsistency score. Two-sided prediction intervals for independent test set samples were constructed for the set nominal confidence level to evaluate the actual coverage and the average interval width.

[0314] The P6 process is implemented through the following steps:

[0315] Step 1: Unbiased Regression Test

[0316] The optimal augmented training set determined by the P3 process is used as the data basis for the unbiasedness test; leave-one-out cross-validation is used to obtain the cross-validation prediction value for each sample; specifically, one sample is reserved as the validation set each time, and the remaining samples are used for training. This process is repeated a certain number of times with the optimal augmented training set samples to obtain the cross-validation prediction value for all samples. ;

[0317] Paired data obtained from leave-one-out cross-validation Establish measured values Compared with the predicted value Linear regression model between:

[0318] (15)

[0319] The intercept 'a' reflects the translation bias of the model (if a ≠ 0, the predicted value has a constant offset), and the slope 'b' reflects the scaling bias of the model (if b ≠ 1, the predicted value is magnified or reduced as a whole). An ideal unbiased model should satisfy a ≈ 0 and b ≈ 1.

[0320] The parameters a and b were estimated using ordinary least squares, and their significance was tested: the intercept a was tested to see if it was significantly non-zero (p<0.05 indicates translation bias); the slope b was tested to see if it was significantly non-equal to 1 (by constructing the t-statistic t=(b-1) / SE(b), p<0.05 indicates scaling bias).

[0321] Step 2: Local weighted regression smoothing analysis

[0322] To further reveal potential nonlinear patterns or heteroscedasticity in the residuals, a locally weighted regression scatter smoothing method was used to fit the trend of residual variation with predicted values; the predicted values ​​were then used as the basis for further analysis. The x-axis represents the residual. Plot a scatter plot with the vertical axis as the ordinate and overlay a LOWESS smoothing curve, as shown below. Figure 7 As shown; if the smooth curve has no systematic trend (such as U-shape, inverted U-shape or monotonic trend) near the zero line, it indicates that the model residuals are uniformly distributed and there is no obvious heteroscedasticity;

[0323] Step 3: Conformal prediction quantifies model uncertainty

[0324] The inductive conformal prediction method is used to quantify the model's prediction uncertainty. The optimal enhancement training set determined by the P3 process is divided into a training subset and an independent calibration set according to the principle of hierarchical random partitioning. The training subset is used for training the PG-BO-HGBR model, and the calibration set is used to construct the non-consistent score distribution. The hierarchical variables still use NMAS and UTI to ensure that the operating condition distribution of the calibration set and the training subset is consistent. The sample size ratio of the training subset and the calibration set is set to 70%:30% based on engineering experience.

[0325] The absolute prediction residuals of the model on the calibration set are used as the inconsistency score to characterize the out-of-sample prediction error distribution of the model; for the nominal confidence level of 1-α commonly used in engineering (90%), the residual quantiles are calculated using the rigorous theoretical formula of conformal prediction.

[0326] (16)

[0327] in, Given the sample size of the independent calibration set, the quantile values ​​of the calibration residuals are calculated based on this quantile, serving as the half-width baseline for the prediction interval; for the independent test set samples, two-sided prediction intervals are constructed for the point prediction values ​​of the two models respectively:

[0328] (17)

[0329] in Let i be the model prediction value for the i-th sample. To calibrate the theoretical quantile values ​​of the residuals;

[0330] Two core evaluation indicators are selected: actual coverage (measures interval reliability, the closer to the nominal confidence level, the better) and average interval width (measures interval accuracy, the smaller the better for the same coverage). The calculation formula is as follows:

[0331] (18)

[0332] in, The actual coverage rate measures the reliability of the prediction interval, with a value between 0 and 1. denoted as the number of samples in the independent test set, and denoted as i as the index number of the test sample. This represents the measured (true) fracture energy of the i-th test sample. Let be the predicted fracture energy value for the i-th test sample. The quantile values ​​calculated based on the calibration set residuals are used as the half-width reference for the prediction interval. This is an indicator function; it takes the value 1 when the condition within the parentheses is true, and 0 otherwise. Let be the two-sided prediction interval for the i-th test sample;

[0333] (19)

[0334] in, The average interval width measures the accuracy of the prediction interval; its unit is the same as the fracture energy (J / m). 2 ), denoted as the number of samples in the independent test set, and denoted as i as the index number of the test sample. This represents the measured (true) fracture energy of the i-th test sample. Let be the predicted fracture energy value for the i-th test sample. The quantile values ​​calculated based on the calibration set residuals are used as the half-width reference for the prediction interval. This is the upper bound of the prediction interval. The lower bound of the prediction interval is used; comparison plots of prediction intervals for the PG-BO-HGBR and unconstrained U-BO-HGBR models on independent test sets are drawn, along with coverage curves at different nominal confidence levels. The results are as follows: Figure 8 As shown.

[0335] P7: The SHAP method is used to perform global and local interpretability analysis and calculate the contribution of each feature to the prediction results. The cumulative local effect (ALE) method is used to analyze the marginal effect of key features. By calculating the average change of the feature in the predicted value within the local interval, the marginal relationship between the feature and the prediction target is revealed unbiasedly, and the physical consistency of the model behavior is verified.

[0336] The P7 process is implemented through the following steps:

[0337] Step 1: SHAP Global Feature Importance Analysis

[0338] Based on the trained PG-BO-HGBR model, the SHAP method is used to calculate the contribution of each feature of each sample to the fracture energy prediction. The SHAP value is derived from the Shapley value in cooperative game theory. For model f, the SHAP value of the j-th feature of sample x is... The calculation formula is:

[0339] (20)

[0340] in, Let be the SHAP value of the j-th feature, and N be the set of all features. Let S be the total number of features, and S be the subset of features that does not contain feature j. Let f(S) be the size of the subset, and f(S) represent the predicted value based on the feature subset S. Based on the SHAP values ​​of all samples, the mean absolute SHAP value of each feature is calculated as a global importance indicator.

[0341] (twenty one)

[0342] in, Let be the global importance score of the j-th feature, and n be the total number of samples. Let be the SHAP value of the j-th feature of the i-th sample; sort each feature by importance score from high to low, and draw the SHAP dependency graph and bee colony graph, as shown below. Figure 9 As shown, the influence of features on prediction results is revealed by plotting feature values ​​on the x-axis and SHAP values ​​on the y-axis.

[0343] Step 2: Verification of physical consistency of cumulative local effects (ALE)

[0344] The cumulative local effect method is used to analyze the marginal effect of key features on fracture energy prediction; for feature x j The range of values ​​is divided into K equal-frequency intervals. The local effect is calculated in the k-th interval, and the accumulated effect is the ALE value.

[0345] (twenty two)

[0346] in, Features In taking values The cumulative local effect value at the location, To determine the number of interval divisions, Features The value falls on the A sample set within a certain interval The number of samples in this set. For the first The upper boundary value of each interval. For the first The lower boundary value of each interval. For the i-th sample, excluding features Other feature values, The model prediction function is used; the ALE curve is verified by comparing the ALE curves of the PG-BO-HGBR and unconstrained BO-HGBR models on the three features of IFSS, IFBS, and E (e.g., ...). Figure 6 (As shown) Whether it strictly maintains positive monotonicity throughout the entire range of feature values;

[0347] Step 3: Identification of key influencing factors and analysis of their mechanisms of action

[0348] By integrating the results of SHAP and ALE analyses, the contribution ranking of each feature to fracture energy prediction was clarified. The nonlinear action modes of each key variable were identified through the SHAP dependency graph, including the critical threshold of UTI, the monotonically increasing law of IFSS and IFBS, and the negative influence mechanism of NMAS. The marginal response of the model within the value range of each feature was verified by combining ALE curves to see if it conforms to the physical mechanism of fiber-reinforced composite materials, thus providing a theoretical basis for optimizing the crack resistance of fiber-reinforced asphalt mixtures.

[0349] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A machine learning method for predicting the fracture energy of fiber-reinforced asphalt mixtures under small sample conditions, characterized in that: Includes the following steps: P1: The obtained samples related to the fracture energy of fiber-reinforced asphalt mixtures are divided into training set and test set according to the stratified random partitioning method; the key features affecting the distribution of working conditions are used as stratification variables to divide the working condition into several subsets. Each subset randomly selects samples according to a preset ratio and assigns them to the training set and test set respectively. After merging, the training set and independent test set are obtained. The test set is locked throughout the entire process and does not participate in subsequent data generation, model training, or hyperparameter optimization. P2: A Wasserstein generative adversarial network and gradient penalty are used to augment the training set. The generator takes noise and a conditional vector as input and outputs data containing both input and output features. The discriminator takes a concatenation of real / generated data and the conditional vector as input and outputs a discrimination score. A gradient penalty mechanism is used to satisfy the 1-Lipschitz continuity constraint. The consistency of the distribution between generated and real samples is verified by t-SNE visualization. P3: Using Bayesian optimized histogram gradient boosting regression tree as the evaluation benchmark, ablation experiments with different generation ratios were conducted; model training and testing were repeated multiple times at each ratio, and the mean and standard deviation of the evaluation index of the test set were statistically analyzed to determine the optimal generation ratio and the corresponding augmented training set size. P4: Based on the optimal augmented training set, repeated randomized experiments were conducted to systematically compare multiple candidate models, including linear baseline models, ensemble baseline models, BP neural networks, and histogram gradient boosting regression trees; each model was optimized using a differentiated hyperparameter search strategy, and the optimal baseline model was determined based on prediction accuracy and stability. P5: Based on Spearman correlation analysis and fracture mechanics theory, physical monotonic constraints are applied to key features. The constraint vector is then passed into the optimal baseline model to ensure that the model's predicted response to the specified features strictly maintains the preset monotonicity throughout the entire input space, thus constructing a physical-guided machine learning model. P6: Leave-one-out cross-validation was used to test the unbiasedness of the physical guidance model. A linear regression model of measured and predicted values ​​was established to test the significant differences between the intercept and 0 and the slope and 1. The inductive conformal prediction method was used to quantify the model prediction uncertainty. The absolute prediction residuals on the calibration set were used as the inconsistency score. Two-sided prediction intervals for independent test set samples were constructed for the set nominal confidence level to evaluate the actual coverage and the average interval width. P7: The SHAP method is used to perform global and local interpretability analysis, and the contribution of each feature to the prediction results is calculated. The cumulative local effect method is used to analyze the marginal effect of key features. By calculating the average change of the feature in the predicted value within the local interval, the marginal relationship between the feature and the prediction target is revealed unbiasedly, and the physical consistency of the model behavior is verified.

2. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions according to claim 1, characterized in that: The P1 process is implemented through the following steps: Step 1: Data Collection Fracture energy data of fiber-reinforced asphalt mixtures were obtained through multi-scale tests, with a total of 48 valid samples collected. Seven parameters were selected as input features, and fracture energy was used as the output feature. The specific parameter statistical analysis is shown in Table 1 below. Table 1 Statistical Analysis of Modeling Parameters ， Step 2: Determining Hierarchical Variables and Dividing Working Condition Subsets To prevent imbalance in the distribution of working conditions between the training and test sets, the maximum nominal aggregate size (NMAS) and the effective temperature range (UTI) of asphalt were selected as stratified variables. NMAS values ​​included two categories: 13.2 mm and 19.0 mm, and UTI values ​​included two categories: 86℃ and 98℃. The 48 samples were divided into four subsets based on the combination of these two variables: Subset A, NMAS = 13.2 mm, UTI = 86℃; Subset B, NMAS = 13.2 mm, UTI = 98℃; Subset C, NMAS = 19.0 mm, UTI = 86℃; Subset D, NMAS = 19.0 mm, UTI = 98℃. Each subset had 12 samples, with a uniform distribution. Step 3: Combine the real training set and the independent test set Within each working condition subset, samples are randomly selected in a ratio of 80%:20% and assigned to the training set and the test set respectively; that is, 10 sets of samples are randomly selected from each subset to enter the training set, and the remaining 2 sets of samples are entered into the test set. The training samples extracted from the four subsets were combined to obtain 38 training sets; the test samples extracted from the four subsets were combined to obtain 10 independent test sets. The test set is locked throughout the entire process and does not participate in subsequent data generation, model training, or hyperparameter optimization. Step 4: Validation of the partition Verify the consistency of the distribution of continuous features between the training set and the test set; Calculate IFSS, IFBS, The mean and standard deviation of features such as OAC and E are shown in Table 2 below. The results show that the relative deviation of the mean between the training set and the test set does not exceed 10%, and the relative deviation of the standard deviation is within 15%, indicating that the two sets of data have good homogeneity in feature distribution. The stratified random partitioning effectively ensures the consistency of the working condition distribution between the training set and the test set. Table 2. Statistical Comparison of Continuous Features between Training and Test Sets ， Note: Relative bias = |Training set mean - Test set mean| / Training set mean × 100%.

3. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions according to claim 2, characterized in that: The P2 process is implemented through the following steps: Step 1: Training Data Preparation and Conditional Vector Construction The 38 training sets obtained in the P1 process were used as the augmentation baseline for WGAN-GP; for each set of samples, seven input features were extracted: IFSS, IFBS, UTI, ... OAC, E, NMAS, and one output characteristic: fracture energy G f This forms an 8-dimensional original data vector. Simultaneously, a condition vector is constructed based on the material properties of each sample group: four fiber types (non-fiber, basalt fiber, glass fiber, and polyester fiber) are converted into 4-dimensional vectors using unique thermal encoding; two asphalt types (70 base asphalt and SBS modified asphalt) are converted into 2-dimensional vectors; and two gradation types (AC-13 and AC-20) are converted into 2-dimensional vectors. These three vectors are then concatenated to obtain an 8-dimensional condition vector. Step 2: Constructing the Generator and Discriminator Networks The generator and discriminator network structure of WGAN-GP is constructed, and the specific configuration is shown in Table 3 below; Table 3. Network Structure of WGAN-GP Generator and Discriminator ， Note: 1 Noise dimension 100 + condition vector dimension 8 = 108; 2 Data dimensions: 7 input features + 1 output fracture energy G f = 8-dimensional; 3 Data dimension 8 + condition vector dimension 8 = 16; The generator learns the conditional distribution of real data by combining random noise and conditional vectors. The discriminator is used to distinguish between real samples and generated samples. The larger its output value, the more likely the input is to come from the real distribution. Step 3: Gradient penalty mechanism and loss function design A gradient penalty mechanism is employed to satisfy the 1-Lipschitz continuity constraint; during each discriminator update, the real sample x and the generated sample x are compared. Perform linear interpolation to obtain interpolated samples. ,in Uniformly distributed random numbers, Calculate the gradient norm at the interpolated samples and construct the gradient penalty term. , where the gradient penalty coefficient λ=10; The discriminator loss function is defined as: （1） in: The discriminator loss function value. To generate sample distribution Expectations For the discriminator to generate samples The output score, To determine the distribution of real data Expectations Let λ be the output score of the discriminator for the real sample x, and λ be the gradient penalty coefficient. For the discriminator output of interpolated samples gradient, The L2 norm of the gradient; The generator loss function is defined as: （2） in: The generator loss function value. To generate sample distribution Expectations For the discriminator to generate samples The output score; Step 4: Adversarial Training and Parameter Update The Adam optimizer was used, with the learning rates for both the generator and discriminator set to 0.0002, Adam-β1=0.5, and β2=0.

9. An alternating training strategy was adopted, with the generator updated once after every 5 updates to the discriminator. The batch size was set to 16, and the total number of training rounds was 30,000. Step 5: Generate sample quality verification The t-SNE algorithm is used to reduce the dimensionality of real samples and generated samples to three-dimensional space for visualization. The distribution of 500% scale samples generated at Epoch=3000, 10000, 20000 and 30000 is compared with that of real samples. A 95% confidence ellipsoid is drawn to represent the core region of the sample distribution.

4. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions, as described in claim 3, is characterized in that: The P3 process is implemented through the following steps: Step 1: Determine the evaluation benchmark model Bayesian optimized histogram gradient boosting regression tree was adopted as the unified evaluation benchmark for ablation experiments. This model has been fully validated in the field of asphalt mixture mechanical property prediction, with high prediction accuracy and stability. It can eliminate the interference of model selection on the optimization results of the generation amount and ensure the consistency of evaluation standards. Step 2: Construct augmented training sets with different proportions Based on the WGAN-GP model trained in the P2 process, data augmentation was performed on 38 sets of real training data according to a preset generation ratio, constructing 6 augmented training sets of different sizes: 0% ratio: Only 38 real training sets were used. 100% ratio: 38 real samples + 38 generated samples, totaling 76 samples. 200% ratio: 38 real samples + 76 generated samples, totaling 114 samples. 300% ratio: 38 real samples + 114 generated samples, totaling 152 samples. 500% ratio: 38 real samples + 190 generated samples, totaling 228 samples. 800% ratio: 38 real samples + 304 generated samples, totaling 342 samples; Step 3: Bayesian hyperparameter optimization For each scaled augmented training set, a Bayesian hyperparameter search is performed using the Optuna framework, with the optimization objective being the root mean square error (RMSE) on the validation set. The hyperparameter search space is defined as follows: Number of trees: 50-300, integer type Maximum depth: 3~12, integer type Learning rate: 0.01–0.3, log-uniform distribution Sample size ratio: 0.6–1.0, uniformly distributed. Feature sampling ratio: 0.6–1.0, uniformly distributed. Minimum number of samples for leaf nodes: 2-20, integer type Each hyperparameter optimization is evaluated through 50 iterations, and the optimal parameter combination is selected for subsequent training. Step 4: Design and execution of repeated experiments To evaluate the statistical stability of the model's performance, repeated experiments were conducted for each generation ratio, as follows: The optimal combination of hyperparameters determined by Bayesian optimization with this fixed ratio The model training and testing process was repeated 30 times, with a different random seed each time, to ensure that the partitioning of the training and testing sets was random. Calculate the prediction performance metrics, including the coefficient of determination R, for each iteration on the independent test set. 2 Root mean square error (RMSE) and mean absolute error (MAE) The mean and standard deviation of 30 repetitions were recorded as the final evaluation result of the model's generalization performance at this generation ratio. Step 5: Determine the optimal generation ratio For each generation ratio, 30 repeated experiments were conducted, with different random seeds used for model training and testing in each experiment. The coefficient of determination R on the test set was then calculated. 2 The mean and standard deviation of the root mean square error (RMSE) and mean absolute error (MAE) are used to comprehensively evaluate the impact of the generated sample size on model performance by comparing the prediction accuracy and stability of the model at different scales.

5. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions according to claim 4, characterized in that: The P4 process is implemented through the following steps: Step 1: Candidate Model Selection Based on the optimal augmented training set determined by the P3 process, six representative models were selected for systematic comparison, including the linear baseline model Ridge Regression, the ensemble baseline model Random Forest, the BP neural network model, and the gradient boosting tree model. The objective function of the ridge regression model is: （3） in: For the regression coefficient vector, For the target variable vector, For the input feature matrix, The sum of squared residuals measures the accuracy of the model fit. This is the regularization coefficient, which controls the balance between model complexity and fitting accuracy. This is an L2 regularization term that penalizes excessively large coefficients to prevent overfitting. The output of the random forest prediction is: （4） in: These are the predicted values ​​from the random forest. For the number of decision trees, Let t be the prediction result of the t-th decision tree for input x; The formula for calculating the hidden layer output of a BP neural network is: （5） in: This represents the output of the j-th neuron in the hidden layer. ( ) represents the hidden layer activation function. The connection weights are the connection weights from the i-th neuron in the input layer to the j-th neuron in the hidden layer. Let be the input value of the i-th neuron in the input layer, and n be the number of neurons in the input layer, which is 7. The threshold value is the value of the j-th neuron in the hidden layer. The formula for calculating the output of the BP neural network output layer is: （6） in: The output of the k-th neuron in the output layer. ( ) represents the activation function of the output layer. The connection weights from the j-th neuron in the hidden layer to the k-th neuron in the output layer are... This represents the output of the j-th neuron in the hidden layer, where m is the number of neurons in the hidden layer, which is 8. The threshold value for the k-th neuron in the output layer; The formula for calculating the overall error of a BP network is: （7） Where: E is the overall network error, and D is the desired output vector. This is the actual output vector. Let be the expected value of the k-th output neuron. This represents the actual value of the k-th output neuron. The formula for the log-sigmoid activation function is: （8） in: Let x be the output value of the activation function, x be the input value of the activation function, and e be the natural constant. Histogram gradient boosting regression trees are based on the forward distribution algorithm, and their additive model can be expressed as: （9） in: For the HGBR model, the input Predicted value The total number of decision trees, Let the weight of the m-th tree be . This represents the prediction result of the m-th regression tree for the input x; Step 2: Data partitioning and experimental setup Thirty repeated random partitioning experiments were conducted to eliminate the influence of a single data partition on the results. In each experiment, the optimal augmented training set was randomly divided into a training subset and a validation subset in a ratio of 80%:20%. The training subset was used for model training, and the validation subset was used for performance evaluation during the hyperparameter optimization process. Step 3: Hyperparameter Optimization Strategy Different hyperparameter search strategies were designed to suit the characteristics of different models. For the linear model Ridge and the ensemble models RF and HGBR, Bayesian optimization was used to efficiently find the optimal parameters within the preset space. For the BP neural network, genetic algorithm, Bayesian optimization, and particle swarm optimization were introduced to perform global hyperparameter search to fully explore its nonlinear fitting potential. The hyperparameter search space settings for each model are shown in Table 4 below. Table 4. Model Hyperparameter Search Space ， Step 4: Comparison of model training results and determination of the optimal baseline model By comparing the mean values ​​and box plot distributions of the indicators from 30 repeated experiments for each model, the predictive performance and stability of the six models were evaluated, and the optimal baseline model was finally determined. The formulas for calculating each evaluation indicator are as follows: （10） in: The coefficient of determination measures how well the model explains the variance of the target variable; its value ranges from [0,1]. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample. The average of the measured values. ; （11） in: The mean squared error is calculated by squaring the error, which amplifies the weight of larger errors. n is the number of samples. Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample; （12） in: It is the root mean square error, monotonically related to the mean square error (MSE), and retains the same dimensions as the original response variable. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample; （13） in: The mean absolute error directly reflects the average magnitude of the prediction error. For the sample size, Let be the measured fracture energy of the i-th sample. Let be the predicted fracture energy value for the i-th sample.

6. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions according to claim 5, characterized in that: The P5 process is implemented through the following steps: Step 1: Feature Correlation Analysis and Determination of Constraints Based on the optimal enhanced training set determined by the P3 process, the Spearman rank correlation coefficient is used to calculate and quantify the relationship between each input feature and the fracture energy G. f The degree of monotonic correlation between them; Based on the fracture mechanics theory of fiber-reinforced composites: interfacial shear strength (IFSS) reflects the load transfer capacity of the fiber-mortar interface; the higher the value, the more energy the fiber consumes during pull-out. Interfacial bond strength (IFBS) characterizes the adhesion performance between the mortar and the aggregate surface; strengthening the interfacial bond can delay crack propagation along the mortar-aggregate interface. Fiber modulus (E) characterizes fiber stiffness; high-modulus fibers can more effectively transfer stress in the mixture, limiting crack opening. Therefore, it is determined that the above parameters and fracture energy should satisfy a positive monotonic relationship, that is, the fracture energy does not decrease monotonically as the parameter values ​​increase. Step 2: Design of Physical Monotonic Constraint Vectors Based on the analysis in step 1, positive monotonic constraints are applied to the three core interface parameters IFSS, IFBS, and E. For the remaining features UTI, G, OAC, and NMAS, due to nonlinear interaction effects, threshold effects, or theoretical optimal values, global monotonic constraints are not suitable. The monotonic constraint parameter monotonic is set for the BO-HGBR model. cst For vectors: ， The vector elements correspond to IFSS, IFBS, UTI, G, OAC, E, and NMAS respectively, with 1 indicating a positive monotonic constraint and 0 indicating no constraint. Step 3: Physics-guided model construction The aforementioned constraint vectors are input into the BO-HGBR model to construct a physically guided PG-BO-HGBR model. During model training, the HGBR algorithm checks whether the splitting direction violates the preset monotonic constraint each time a decision tree node splits: for features with positive monotonic constraints, if the predicted mean of the left child node is greater than the predicted mean of the right child node after the current node splits, the split is prohibited, and the algorithm automatically selects other candidate split points. Through this mechanism, it is ensured that the predicted response of each base learner and the final ensemble model to the specified feature strictly maintains the preset monotonic direction throughout the entire input space. The same Bayesian optimization strategy as the unconstrained BO-HGBR is used to optimize the hyperparameters of PG-BO-HGBR; the search space is the same as in step 3 of P3. Step 4: Model training results and generalization ability verification Under optimal hyperparameter configuration, the PG-BO-HGBR model is trained using the optimal augmented training set to obtain the final physically guided fracture energy prediction model; this model is expressed as: （14） Where the function To integrate decision tree models; The predicted and measured values ​​of PG-BO-HGBR on the training and test sets were compared with the results of the unconstrained BO-HGBR model to verify the performance of the PG-BO-HGBR model. In order to evaluate the generalization ability of the physically guided PG-BO-HGBR model in real-world conditions, 10 independent test sets of samples that did not participate in any training were used for verification.

7. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions, as described in claim 6, is characterized in that: The P6 process is implemented through the following steps: Step 1: Unbiased Regression Test The optimal augmented training set determined by the P3 process is used as the data basis for the unbiasedness test; leave-one-out cross-validation is used to obtain the cross-validation prediction value for each sample; specifically, one sample is reserved as the validation set each time, and the remaining samples are used for training. This process is repeated a certain number of times with the optimal augmented training set samples to obtain the cross-validation prediction value for all samples. ; Paired data obtained from leave-one-out cross-validation Establish measured values Compared with the predicted value Linear regression model between: （15） The intercept 'a' reflects the translational bias of the model, and the slope 'b' reflects the scaling bias of the model. An ideal unbiased model should satisfy a≈0 and b≈1. The parameters a and b are estimated using ordinary least squares, and their significance is tested: the intercept a is tested to see if it is significantly different from zero; the slope b is tested to see if it is significantly different from 1. Step 2: Local weighted regression smoothing analysis To further reveal potential nonlinear patterns or heteroscedasticity in the residuals, a locally weighted regression scatter smoothing method was used to fit the trend of residual variation with predicted values; the predicted values ​​were then used as the basis for further analysis. The x-axis represents the residual. Plot a scatter plot with the vertical axis as the ordinate and overlay a LOWESS smoothing curve. If the smoothing curve shows no systematic trend near the zero line, it indicates that the model residuals are uniformly distributed and there is no obvious heteroscedasticity. Step 3: Conformal prediction quantifies model uncertainty The inductive conformal prediction method is used to quantify the model prediction uncertainty; the optimal enhancement training set determined by the P3 process is divided into a training subset and an independent calibration set according to the hierarchical random partitioning principle; the training subset is used for PG-BO-HGBR model training, and the calibration set is used to construct the non-consistent score distribution. The stratified variables still use NMAS and UTI to ensure that the operating condition distribution of the calibration set and the training subset is consistent; the sample size ratio of the training subset to the calibration set is set to 70%:30% based on engineering experience. The absolute prediction residuals of the model on the calibration set are used as the inconsistency score to characterize the out-of-sample prediction error distribution of the model; for the nominal confidence level of 1-α commonly used in engineering (90%), the residual quantiles are calculated using the rigorous theoretical formula of conformal prediction. （16） in, Given the sample size of the independent calibration set, the quantile values ​​of the calibration residuals are calculated based on this quantile, serving as the half-width baseline for the prediction interval; for the independent test set samples, two-sided prediction intervals are constructed for the point prediction values ​​of the two models respectively: （17） in Let i be the model prediction value for the i-th sample. To calibrate the theoretical quantile values ​​of the residuals; Two core evaluation indicators are selected: actual coverage rate and average interval width. The calculation formula is as follows: （18） in, The actual coverage rate measures the reliability of the prediction interval, with a value between 0 and 1. denoted as the number of samples in the independent test set, and denoted as i as the index number of the test sample. Let be the measured fracture energy value of the i-th test sample. Let be the predicted fracture energy value for the i-th test sample. The quantile values ​​calculated based on the calibration set residuals are used as the half-width reference for the prediction interval. This is an indicator function; it takes the value 1 when the condition within the parentheses is true, and 0 otherwise. Let be the two-sided prediction interval for the i-th test sample; （19） in, This represents the average interval width, measuring the accuracy of the prediction interval; its unit is the same as the fracture energy. denoted as the number of samples in the independent test set, and denoted as i as the index number of the test sample. Let be the measured fracture energy value of the i-th test sample. Let be the predicted fracture energy value for the i-th test sample. The quantile values ​​calculated based on the calibration set residuals are used as the half-width reference for the prediction interval. This is the upper bound of the prediction interval. The lower bound of the prediction interval is given; plot the prediction interval comparison charts of the PG-BO-HGBR and unconstrained U-BO-HGBR models on the independent test set, as well as the coverage curves at different nominal confidence levels.

8. The method for predicting the fracture energy of fiber-reinforced asphalt mixtures using machine learning under small sample conditions according to claim 7, characterized in that: The P7 process is implemented through the following steps: Step 1: SHAP Global Feature Importance Analysis Based on the trained PG-BO-HGBR model, the SHAP method is used to calculate the contribution of each feature of each sample to the fracture energy prediction. The SHAP value is derived from the Shapley value in cooperative game theory. For model f, the SHAP value of the j-th feature of sample x is... The calculation formula is: （20） in, Let be the SHAP value of the j-th feature, and N be the set of all features. Let S be the total number of features, and S be the subset of features that does not contain feature j. Let f(S) be the size of the subset, and f(S) represent the predicted value based on the feature subset S. Based on the SHAP values ​​of all samples, the mean absolute SHAP value of each feature is calculated as a global importance indicator. （21） in, Let be the global importance score of the j-th feature, and n be the total number of samples. Let be the SHAP value of the j-th feature of the i-th sample; sort each feature from high to low importance score, draw the SHAP dependency graph and the bee colony graph, and use the feature value as the x-axis and the SHAP value as the y-axis to reveal the influence of features on the prediction results; Step 2: Verification of physical consistency of cumulative local effects The cumulative local effect method is used to analyze the marginal effect of key features on fracture energy prediction; for feature x j The range of values ​​is divided into K equal-frequency intervals. The local effect is calculated in the k-th interval, and the accumulated effect is the ALE value. （22） in, Features In taking values The cumulative local effect value at the location, To determine the number of interval divisions, Features The value falls on the A sample set within a certain interval The number of samples in this set. For the first The upper boundary value of each interval. For the first The lower boundary value of each interval. For the i-th sample, excluding features Other feature values, The model prediction function is used; by comparing the ALE curves of the PG-BO-HGBR and unconstrained BO-HGBR models on the three features IFSS, IFBS, and E, it is verified whether the ALE curve strictly maintains positive monotonicity throughout the entire feature value range. Step 3: Identification of key influencing factors and analysis of their mechanisms of action By integrating the results of SHAP and ALE analyses, the contribution ranking of each feature to fracture energy prediction was clarified. The nonlinear action modes of each key variable were identified through the SHAP dependency graph, including the critical threshold of UTI, the monotonically increasing law of IFSS and IFBS, and the negative influence mechanism of NMAS. The marginal response of the model within the value range of each feature was verified by combining ALE curves to see if it conforms to the physical mechanism of fiber-reinforced composite materials, thus providing a theoretical basis for optimizing the crack resistance of fiber-reinforced asphalt mixtures.

Images