The test instrument adopts the frozen soil dynamic and static triaxial tester FST-250 designed by Xi’an Lichuang Material Testing Technology Co., Ltd. (see Figure 1). The instrument adopts displacement control mode and is mainly composed of triaxial pressure system, XT5718ULT-E2000 cooling water circulating device, data acquisition system, cooling system of testing machine, etc. The temperature control range is -70°C∼+90°C, the temperature control accuracy is 0.5°C, the confining pressure range is 0～40MPa, the accuracy is 1%, and the maximum axial pressure is 250KN.

The microcomputer-controlled electro-hydraulic servo triaxial testing machine for frozen soil is mainly used for dynamic and static triaxial tests of conventional soil and frozen soil. The basic physical properties of the test soil are shown in Table 1. In order to reduce the heterogeneity and dispersion of natural silty clay, 1.25 mm sieve is used to screen out impurities such as plant fiber and gravel of natural silty clay. The distribution of soil with particle size less than 1.25 mm is shown in Table 2. Direct shear test is used to measure cohesion and internal friction angle of soil samples. The liquid limit and plastic limit in Table 1 are measured by LP-100D digital display soil liquid and plastic limit tester. The standard consolidation method is adopted in the consolidation test, and the compressibility factor and modulus are obtained. The maximum and minimum values are obtained from experiments, and the average values are obtained from calculation. According to the liquid plastic limit and particle size distribution, the test soil can be classified as silty clay.

The dry soil that has passed through 1 mm sieve is sprayed and mixed with soil according to the initial moisture content (15%, 17.5%, 20%, 22.5%, and 25%). Then the soil samples are sealed with plastic wrap for 24 h. The dry density of soil sample is controlled to 1.62 g/cm3. Soil samples with different water contents are divided into five layers and compacted into a three-part mold with an inner diameter of 61.8 mm and a height of 125 mm. After standing for a certain time, soil samples with different initial moisture contents are demolded. Then, the soil samples are packaged and labelled with plastic wrap. Finally, the soil samples are put into the high and low temperature alternating test chamber TEMP880.

The high and low temperature alternating test chamber TEMP880 designed by Shanghai-Huayi Equipment Co., Ltd. is used for the freeze-thaw cycle test. The test chamber is set at -35°C–35°C. The specific temperature time history curve in the high and low temperature alternating test chamber during the experiment is shown in Figure 2 below. The soil samples undergoe 4 h of freezing and 4 h of thawing, that is, a freeze-thaw cycle is performed every 8 h (Teng et al., 2022). Then, the soil samples are packaged with a thin latex sleeve, and place in a constant temperature test chamber with the same design temperature for 24 h after sticking the label. YH-10 aviation hydraulic oil in the triaxial pressure chamber is used to apply confining pressure and maintain temperature stability. Next, the soil samples with sufficient temperature control in the constant temperature test chamber are placed in the triaxial pressure chamber. When the temperature of the soil sample and hydraulic oil is coupled to the design temperature, the design confining pressure is applied for the test.

Because the time scale of soil freezing and thawing process in seasonally frozen soil area is considerably large, the reduced similarity ratio is adopted in the test. The similarity criteria based on the similarity theoretical model are as follows (Guo et al., 2010): F ο = a L 2 Γ (1) Where L represents the characteristic length, m. Γ represents the period of surface temperature change, h a represents the proportional coefficient.

The physical simulation characteristics of frozen silty clay are consistent with the test process and the actual process. The experiment lasted 8 h (4 h of freezing and 4 h of thawing), and the 1-year change process under the action of freezing and thawing cycles was simulated. The similarity calculation process is shown in Table 3:

Constant deformation loading is adopted during the test, and the shear rate is 1 mm·min-1. The test is terminated when the specimen strain reached 20% or when it failed. At -15°C, Triaxial tests are carried out on frozen silty clay nuder five freeze-thaw cycles, five initial moisture contents and three confining pressures. The change of soil sample stress is monitored in real time and recorded. After the triaxial test, the frozen soil sample is convex in the upper part and has expansion deformation in the radial direction. The frozen silty clay generally presents an oblique crack failure form (Xu et al., 2020; Niu et al., 2022; Zhang et al., 2022).

Long short-term memory neural network (LSTM) was created primarily to solve the problem of gradient disappearance and explosion. The biggest difference between LSTM and traditional neural network is the introduction of Forget gate, Input gate and Output gate. Another key is the change in the internal state, which is the core of each activated neuron and can be seen as a carrier for adding information or deleting information.

Step 1

Some information is selectively discarded by the forget gate, and LSTM receives the input of C t and h t − 1 . And output a number between 0 and 1 to the internal state S t . f t = σ W f C C t + W f h h t − 1 + b f (2)

Step 2

Some information in the internal state is determined to be stored. The input gate determines which values are updated, and then a candidate state S t is built by the input node g t . i t = σ W i C C t + W i h h t − 1 + b i   (3) g t = tanh W g C C t + W g h h t − 1 + b g   (4)

The original internal state is updated to S t by Eqs 2–4. S t = g t × i t + S t − 1 × f t   (5)

Step 3

The output information is determined by the output gate. o t = σ W o C C t + W o h h t − 1 + b o   (6)

Step 4

Based on the output gate and the new internal state, the new hidden layer state h t is output. h t = tan ⁡ h S t × o t   (7) where W and b represent the weight and bias, respectively, W f C is the weight of C t in the forget gate, and b f is the forget gate bias.

L1 and L2 regularization methods are often used in neural networks to prevent model overfitting. Some parameters in the loss function are limited by this method (Jiang et al., 2018; Huang et al., 2020). Neurons and fully connected layers are randomly deleted through Dropout technology to form new samples (Srivastava et al., 2014). The overfitting problem is effectively solved by this sample replacement effect. At the same time, the entire network is no longer very sensitive to the specific weights of neurons. Therefore, the generalization ability of the model has been greatly improved (Huang et al., 2020; Chang et al., 2022).

The steps of PSO-LSTM are as follows, and the flow chart is shown in Figure 3.

1) The dataset is divided into LSTM training set, validation set and test set.

According to the optimization parameters and constraints set above, the particle swarm optimization algorithm operation is performed. The fitness change curve is shown in Figure 4. It can be seen from Figure 4A that when the iteration exceeds 4 times, the fitness value quickly reaches the maximum. Therefore, the optimal solution is selected when the number of iterations (search steps) is 0–4, and the local optimal solution and the global optimal solution are output for each cycle. After 20 iterations, the global optimal solution can be obtained: [Number of neurons in the first layer, dropout, batch_size]=[61, 0.09, 95]. Therefore, PSO-LSTM structure optimization combination parameters are obtained and input.

It can be seen from the MSE attenuation curve of PSO-LSTM model in Figure 4B above that when the sample reaches 7.5 epochs, the MSE value tends to be stable. It shows that the validation set of the model has reached a good accuracy at this time. After determining the structural optimal parameters of PSO-LSTM, the weights and thresholds of each model are iteratively trained using the training set to obtain the final model. The test set is used to predict and compare each model (Chang et al., 2020; Huang et al., 2022).

75 groups of triaxial tests are conducted on frozen silty clay with initial moisture content(IMC) of 25%, 22.5%, 20%, 17.5%, 15%, freeze-thaw cycles (FTC) of 0, 3, 6, 9, 18 times and confining pressure(CP) of 0.1Mpa, 0.2Mpa, 0.3Mpa. It can be seen from Figure 5 that freeze-thaw cycles and confining pressure do not change the strain characteristics of frozen silty clay. With the increase of initial moisture content, there is hysteresis phenomenon in the failure strain of frozen silty clay.

It can be seen from Figure 5 that when the initial moisture content is 22.5% and 25%, the deviatoric stress shows three development stages with the increase of the axial strain: linear elastic growth stage, elastic-plastic development stage, and strain hardening stage. In the linear elastic growth stage, the deviatoric stress is dominated by cementation elements composed of pore ice and soil particles. Due to the strain is small, the pore ice does not crush and melt, and the structure does not have cracks. At this time, the overall deformation modulus is large, and the deviatoric stress increases linearly.

In the elastic-plastic stage, due to the large strain at this time, micro-cracks are generated in the pore ice. Accompanied by the pressure-melting phenomenon at the ice crack, part of the cementation element structure is destroyed and transformed into friction elements. Since the deformation modulus of friction element is smaller than that of cementation element, the overall deformation modulus decreases. The deviatoric stress increases slowly with the increase of axial strain. During the strain hardening stage, the number of cementation elements remains stable, and the friction element, as the main bearing structure, show linear growth trend with the increase of strain.

It can be seen from Figure 5 that when the initial moisture content is 15%, 17.5%, and 20%, the frozen silty clay shows strain softening properties. The pore ice is broken and local shear bands are formed. The number of cementing elements decreases rapidly, the number of friction elements increases, and the integrity of the structural block is destroyed. The soil samples shows strain softening state. And with the decrease of initial water content, the strain softening phenomenon is more obvious. And under the same working condition, different confining pressures do not change the evolution law of deviatoric stress.

The simplified frozen soil binary medium model based on damage parameters is shown in formula (8). The symbolic interpretation and derivation process of formula (8) are given (Zhang et al., 2022). σ 1 − σ 3 = 1 − E i − E E i − α E f E i ε 1 + E i − E E i − α E f ⋅ ε 1 a + b c ε 1 (8)

Since the research object in this paper is frozen silty clay, the local strain coefficient c is obtained from the stress-strain curve. The comparison between the fitted value of the formula and the actual test value is shown in Figure 6.

It can be seen from Figure 6 that the theoretical values of the constitutive model of frozen soil binary medium are very close to the measured values under different parameters. In addition, there is a high matching degree in the linear elastic stage, elastic-plastic stage and strain hardening stage of the stress strain curve. Therefore, it is further verified that the model can better simulate the curve relationship between stress and axial strain of silty clay in the low-temperature triaxial test. The process parameters of the binary medium model fitted by the frozen silty clay test are shown in Table 4.

In order to compare the calculation accuracy of the model and select the optimal model, the frozen silty clay binary medium model, random forest model, PSO-LSTM model proposed in this paper and the frozen soil triaxial test are compared and analyzed. The stress sharing curve of the cementation element and the friction element is made by Formula (1), and the specific comparison is shown in Figure7. The deviatoric stress-strain curve under IMC25-FTC18-CP0.2 is shown in Figure 7. It can be seen from Figure 7 and formula (1) that the deviatoric stress curve of frozen silty clay is formed by superposition of the stress sharing curve of cementation element and the stress sharing curve of friction element.

In the linear elastic stage, the evolution trend of deviatoric stress of frozen silty clay is more consistent with the sharing curve of cementation element. At this time, the deviatoric stress shared by the friction elements is less than 0.039 MPa. Therefore, the pore ice inside the cementation element is not crushed and melted in the linear elastic stage. The integrity of cementation element is good, and the cementation element mainly plays the role of bearing force under low temperature triaxial action.

In the elastic-plastic stage, the stress sharing curve of friction element shows a linear growth trend, and the stress sharing curve of cementation element begins to grow slowly. When the strain exceeds 2.7%, the difference between the stress sharing curve of cementation element and the total stress curve gradually increases. In the strain hardening stage, the stress sharing curve of cementation element tends to be stable. It shows that cementation element has yielded and stabilized at this time, and the stress sharing curve of the friction element still maintains linear growth trend. It shows that pore ice in the cemented element is crushed and melted under the action of low temperature triaxial. The local shear bands are formed, and the soil sample presents the failure form of oblique crack. In addition, it can be seen from Figure 7 that in the range of axial strain of 20%, the stress sharing curve of cementation element is always larger than stress sharing curve of the friction element. The peak value of the total stress is also the maximum value of the superposition of the stress sharing curve of the cementation element and the friction element.

It can be seen from Figure 7 that the calculated values of the frozen silty clay binary medium model and PSO-LSTM model are completely close to the test values. The constitutive relation of frozen silty clay can be better simulated by PSO-LSTM model. Because excellent matching degree can be achieved in linear elastic stage, elastic-plastic stage and strain hardening stage. The calculated values of the RF model and the test values show a high consistency in the strain hardening stage, but there is a large difference between the linear elastic stage and the elastic-plastic stage, which indicates that the model has a large error and still needs to be adjusted.

Based on the random forest model test set, the stress-strain relationship of frozen silty clay with different moisture content after the 18th freeze-thaw cycles is obtained In the constitutive relation curve obtained by selecting the random forest model, the working condition with lower R² is only shown in Figure 8. Therefore, the worst fitting condition of random forest model is selected as the test set of PSO-LSTM model. Finally, the universality and applicability of PSO-LSTM model can be well verified.

It can be seen from the results in Figure 8 that the strains corresponding to the maximum deviatoric stress errors of the RF model are 18.18%, 7.92%, 4.79%, and 15.04% under the initial moisture contents of 22.5, 20, 17.5, and 15%. The deviatoric stresses corresponding to the RF model are 6.16, 6.50, 5.76, and 2.07 MPa respectively, and the deviatoric stresses corresponding to the test are 5.80, 6.10, 5.45, and 2.41MPa respectively. The errors are 5.8%, 6.1%, 5.38%, and 16.4% respectively. Therefore, the trend of the random forest model curve is consistent with that of the test curve, but there is always a large deviation in the whole evolution process of the partial stress. The evaluation indexes of the initial training set and test set of the RF model are significantly different, and the training set is more excellent. It shows that RF model has over-fitting phenomenon in the training process, and the unsatisfactory test results are generated. Therefore, cross-validation is added to avoid overfitting. In the end, the maximum error of the model is still 16.4%.

The maximum error of PSO-LSTM model is only 1.56% under IMC20FTC18CP0.3 working condition. Compared with RF, the predicted value of the PSO-LSTM model constructed in this paper fits better with the actual measured value, and the relative error is smaller and more stable.

As can be seen from Figure 8, the evolution of deviatoric stress under each condition can be well simulated. The constitutive relationship predicted by the random forest model will have fluctuation segments, and it is prone to the problem of low strain sensitivity. The PSO-LSTM model can well simulate the evolution trend of deviatoric stress in each working condition.

Analysis of the data in Table 5 shows that the R² of the deviatoric stress prediction data of the frozen silty clay obtained by the RF model, the binary medium model and the PSO-LSTM model are all above 94%. Therefore, in the process of iterative training, the internal test laws of initial moisture content, freeze-thaw degradation effect and confining pressure are well excavated, and more reasonable constitutive prediction values are output. However, compared with the RF model and frozen soil binary medium model, the R² obtained by the PSO-LSTM model reaches 98.44%. At the same time, MSE and MAE are the lowest, which are 0.0078 and 0.01778 respectively, indicating that the model prediction effect is the best. The evaluation index of the initial moisture content of 25 is only shown in Table 5. Therefore, the evaluation index of the model test set under other conditions is given in the following Figure 9.

It can be seen from Table 6 that the accuracy of the optimal model under each working condition is above 90%. It can be seen from the comparison model that the accuracy of the PSO-LSTM model in each working condition is higher than that of the RF model. At the same time, RMSE and MAE errors are lower than RF model in all conditions, which is enough to show that the optimal model can well simulate the constitutive relation of frozen silty clay in all conditions.

The evaluation indicators of the RF and PSO-LSTM models on the validation set and the test set are relatively consistent, so the evaluation indicators of the validation set are omitted here. At the same time, the RF model and the PSO-LSTM model have not been fitted. In addition, the R² of PSO-LSTM for frozen silty clay constitutive prediction in the test set are all greater than 0.99. At the same time, RMSE, MSE and MAE are the smallest under each water content working condition, and they are the best in all networks. The better fitting accuracy and prediction performance of this model have been fully proved, and it is very accurate for the constitutive prediction of frozen silty clay. It also provides a new method for the prediction of frozen silty clay constitutive model.

As shown in Figure 10, the strength tends to be stable after 6, 9, 6, 9, and 9 freeze-thaw cycles at initial moisture content of 25, 22.5, 20, 17.5, and 15%, respectively. After 18 freeze-thaw cycles, the shear strength of frozen silty clay decreased by 2.66%, 11.85%, 18.83%, 16.79, and 29.02%, respectively. Under the same confining pressure, the shear strength of frozen silty clay decreases first and then tends to be stable with the increase of freeze-thaw cycles. With the increase of freeze-thaw cycles, the bearing structure of frozen silty clay changes from initial coarse particles and pore ice to more deteriorated fine particles and pore ice. That is, the positions of the pore ice in frozen silty clay are occupied by more fine particles. On the one hand, the thickness of pore ice is reduced, and thin pore ice is more likely to expand cracks under low temperature triaxial stress. On the other hand, the unfrozen water film on the particle surface shows lubrication effect. And with the increase of freeze-thaw cycles, the rolling friction and sliding friction between particles increase more, which makes the shear strength decrease more obviously. Under a certain confining pressure, the rolling friction and sliding friction between particles tend to be stable after 6 and 9 freeze-thaw cycles. So the connection, arrangement and stress history of the soil particles are changed by the freeze-thaw cycles. This deterioration effect decreases with the increase of freeze-thaw cycles and tends to be stable after a certain number of freeze-thaw cycles (6 or 9 times).

As shown in Figure 11, the deviatoric stress of frozen silty clay increases with the increase of confining pressure. This is very consistent with the actual project. And with the increase of the initial moisture content, the effect of confining pressure on the strength increase is no longer obvious. This is because in the third stage, the stress sharing curve of the friction element shows a linear growth trend. The stress sharing of the cemented element keeps the yield stable, and with the increase of the confining pressure, the strength of the cemented element for the yield stability is also greater. When the optimum moisture content is reached, the pore ice of cementation element tends to be saturated. At this time, with the increase of the initial moisture content, the expansion of the cracks generated by the thin pore ice is not obvious. Therefore, the influence of the confining pressure on the yield stability of the cemented element is gradually reduced.

As shown in Figure 12, the strength of frozen silty clay is the highest near the initial moisture content of 20%. With the increase of freeze-thaw times and the decrease of confining pressure, the strength of frozen silty clay decreases. Therefore, the most unfavorable working condition is IMC15-FTC18-CP0.1, and the highest strength working condition is IMC20-FTC18-CP0.3.

Compared with Figures 10, 11, it is obvious that the strength difference of frozen silty clay is the largest under different initial moisture content. That is, the initial moisture content has greater influence on the strength of frozen silty clay (Cao et al., 2022; Teng et al., 2022). This is because the effects of freezing and thawing and confining pressure on the connection and arrangement of soil particles and the stress history are not as significant as the initial moisture content. The bearing structures of soils with different initial moisture contents are quite different. As a result, the difference of pore ice thickness is large, and the difference of crack germination and expansion mode is increased. On the other hand, the difference of rolling friction and sliding friction between frozen silty clay particles under different initial moisture content is large, resulting in a large difference in shear strength.

The structural composition of unsaturated frozen soil (Teng et al., 2022) which influences factors of soil particle and pore ice force and soil strength change process is shown in Figure 13 below.

The structural model of frozen silty clay can be seen in Figure 14. According to the concept of the binary medium model (Zhang et al., 2022), cementation element is a stable ice-soil skeleton composed of clay particles and pore ice. The deformation modulus of cementation element is large, and it usually shows better integrity in case of small deformation. Under triaxial stress or low temperature, the pore ice inside the cemented element may be crushed and melted. Even soil particles can be broken. Local shear bands are gradually formed under deviatoric stress. At the same time, the cementation element transforms into the friction element.

There are a large number of cemented elements and friction elements in the frozen silty clay. And under the action of triaxial stress and low temperature, cementation element is gradually transformed into friction element. The strength characteristic difference of frozen soil is closely related to the quantity difference between cementation elements and friction elements.

The cementation element plays a major role in linear elastic stage, and the stress sharing curve of friction element can be ignored. In the elastic-plastic stage, the bearing capacity of the cementation element reaches the yield stability and continuously transforms into the friction element. During the strain hardening stage, the stress sharing curve of the cementation element remains stable, and the stress sharing curve of the friction element continues to increase linearly.

Compared with confining pressure and freeze-thaw cycles, the initial moisture content has the greatest influence on the evolution trend of deviatoric stress. Therefore, machine learning method is used to study the importance of features. It provides a new method for the study of the contribution weight of different factors to frozen silty clay constitutive curve. The contribution weight of deviatoric stress is shown in Figure 15.

The contribution weight of the four factors to the deviatoric stress can be seen from the importance map of random forest characteristics. According to the above experimental study, the trend of the stress-strain relationship curve will be changed by the different initial moisture content. The state of frozen silty clay before the failure will be determined by the initial moisture content, and the position of peak strain will also be changed. With the increase of initial moisture content, the failure stage changes from strain softening to strain stabilization, and then to strain hardening. At the same time, the failure stress will also lag. This also verifies the strong connection between the moisture water content and the stress evolution law.

In a certain range (IMC = 25, 22.5, 20, 17.5, 15%, FTC = 0, 3, 6, 9, 18 times, CP = 0.1, 0.2, 0.3 MPa), the influence factors of frozen silty clay stress from strong to weak are: initial moisture content > strain > confining pressure > freeze-thaw cycles. The contribution of confining pressure to the stress evolution law is small, which is because the freezing depth of the seasonal frozen area is less than 2 m. Therefore, the gradient and range of the confining pressure setting in this paper are small, so the contribution weight is small.

Although the binary constitutive model of frozen silty clay has obtained a good prediction effect. And the reference formula is provided, which is beneficial to the application in practical engineering. However, there are some shortcomings in the model. For example, the binary constitutive model of frozen silty clay in this paper is only suitable for the fitting of the constitutive relationship under the condition of IMC=25%, but it can not achieve a good fitting effect for the IMC=22.5%, 20%, 17.5%, and 15%. Therefore, the formula fitting model still needs to find the constitutive formula of similar curve law. The formula fitting software 1STOPT can also be used to perform a random search of constitutive relations. But some parameters in the formula cannot be explained mathematically and physically. Therefore, in-depth research and promotion are limited.

Based on the shortcoming that the above model is applicable to a single working condition, the RF model is established in this paper. The RF model can be applied to all initial moisture content conditions, but the model prediction accuracy is not satisfactory.

The PSO-LSTM model is established in this paper to solve the problems of the above models. There are many advantages, such as high prediction accuracy, wide application range, and no need to propose based on certain assumptions. The LSTM optimal combination input parameters are searched by PSO, and the parameter adjustment speed of the model for the data learning process of other soil types is greatly increased, which is conducive to popularization. Therefore, it is only necessary to carry out frozen soil triaxial tests of different soil types in various regions. PSO-LSTM model has extremely high accuracy in constitutive prediction, which can provide reference for frozen soil engineering application in cold regions.

Moreover, the deep learning algorithm has a strong fitting ability for multiple operating conditions and multiple factors. The data driven deep learning method will have an excellent fitting effect for more dimensional sample data. For other soil types, the triaxial test data are imported, and the high-precision constitutive relation curve can be obtained by modifying the training set.

In the Code for Design of Soil and Foundation of Building in Frozen Soil Region, the characteristic value fa of foundation bearing capacity is selected based on temperature and soil type. Only four types of frozen soil with ground bearing capacity characteristic values at -3°C∼-0.5°C are given in the specification. However, it is not suitable for saline frozen soil and frozen peat soil, and a large number of experiments need to be carried out to revise the specification. Therefore, the optimal model proposed in this paper has important generalization significance.

That is, a more comprehensive reference table of characteristic values of frozen soil bearing capacity is given by carrying out different working conditions tests on soils in various regions. The working conditions that have not been tested (IMC=16, 19, 21, 24%, FTC=12, 15 times, CP=0.05, 0.15, 0.25, 0.35Mpa) are predicted by the PSO-LSTM model to demonstrate the advantages of the model. As shown in Figure 16, the predicted strength law is consistent with the actual situation. It provides a new method for the constitutive study of frozen soil engineering in cold regions, and provides a favorable reference for the strength of remolded soil in seasonal frozen regions. Therefore, it is of great significance to carry out the model study in this paper.

Random forest and PSO-LSTM can provide a new method for frozen soil constitutive research. The constitutive model of frozen soil is developed based on mathematical skills, which can deeply learn the data of frozen soil triaxial test. Machine learning model has strong nonlinear mapping ability, and is not limited by frozen soil constitutive formula In the process of deep mining of multi-dimensional large sample data, there is usually a better fitting effect. Therefore, machine learning shows extremely high accuracy and wide applicability in frozen soil constitutive model prediction. It can be used to modify the characteristic value fa of foundation bearing capacity of the Code for Design of Soil and Foundation of Building in Frozen Soil Region.