Multi-Lane Trafﬁc Load Clustering Model for Long-Span Bridge Based on Parameter Correlation

: Trafﬁc loads are the primary external loads on bridges during their service life. However, an accurate analysis of the long-term effect of the operating trafﬁc load is difﬁcult because of the diversity of trafﬁc ﬂow in terms of vehicle type and intensity. This study established a trafﬁc load simulation method for long-span bridges based on high authenticity trafﬁc monitoring data, and an improved k -means clustering algorithm and Correlated variables Sampling based on Sobol sequence and Copula function (CSSC) sampling method. The monitoring trafﬁc data collected through a weigh-in-motion (WIM) system was processed to generate a multi-lane stochastic trafﬁc ﬂow. The dynamic response of a prototype suspension bridge under a trafﬁc load was analyzed. The results show that the trafﬁc load can be divided into clusters with identical distribution characteristics using a clustering algorithm. Combined with CSSC sampling, the generated trafﬁc ﬂow can effectively represent daily trafﬁc and vehicle characteristics, which improves the accuracy of the assessment of the loads long-term effect. The dynamic response of the bridge to different trafﬁc ﬂows varied signiﬁcantly. The maximum and minimum vertical displacement of the main girder was 0.404 m and 0.27 m, respectively. The maximum and minimum bending stresses of the short suspender were 50.676 MPa and 28.206 MPa, respectively. The maximum equivalent bending stress and axial stress were 16.068 MPa and 10.542 MPa, respectively, whereas the minimum values were 9.429 MPa and 8.679 MPa, respectively. These differences directly inﬂuence the short and long-term evaluation of bridge components. For an accurate evaluation of the bridge operation performance, the trafﬁc ﬂow density must be considered.


Introduction
With the rapid development of highway transportation, long-span bridges have been constructed across mountains, valleys, and rivers, owing to their mechanical characteristics and excellent spanning performance. The major external loads a bridge must support when it is in use are wind and traffic loads [1]. Long-span bridges have large loading areas and complex load composition, and the reciprocating action of traffic loads may cause a vibration response. The service life of some vibration-sensitive components was significantly reduced. The longitudinal vibration displacement of structures caused by an external load may lead to fatigue of the expansion joints and other ancillary components. The fluctuation in the stress of the suspender may accelerate fatigue degradation. It has been demonstrated that the component fatigue life is closely related to the stress amplitude and cycle numbers of structural vibrations under traffic flow [2,3]. As the intensity of the traffic load increases, the degradation of bridge components has become a crucial factor affecting bridge safety [4]. Thus, an accurate simulation of traffic load is required for long-term evaluation based on the bridge's dynamic response.
The analysis under the design load could not accurately reflect the actual operating condition because of the diversity of the traffic flow in terms of vehicle type and intensity. Traffic load has an impact on the dynamic characteristics and vibration response of bridges. The development of weigh-in-motion (WIM) technology has improved the analysis and evaluation of bridge under traffic loads [5]. For long-span bridges, the analysis using the traffic load, which is closer to the actual service condition, should consider the impact of time, region, and other factors. Several studies have been conducted on traffic load simulation methods. Crespo developed a continuous random traffic flow model for bridges, simulated and generated the relevant traffic flow based on the statistical characteristic of measured traffic data, and extrapolated the long-term effect from the short-term results [6]. Chen pointed out that vehicle load excitation has a spatial correlation that affects the response spectrum obtained from spatially uncorrelated white noise excitation. The Poisson distribution can be used to simulate the vehicle flow load [7]. Chen simplified the random traffic load into a moving load column, and analyzed the impact of loading length and vehicle type on the load response of long-span bridges based on cellular automata theory [8]. Jabri considered the randomness of driver selection in the distribution of vehicle headways, ensuring nonnegative traffic flow density, and further improved the cellular automata-based random traffic flow simulation method [9]. After extensive development, the bridge random traffic load was used to simulate the vehicle mass parameters, traffic load characteristics, and driving behavior.
Traffic flow simulation methods primarily use macro and micro traffic flow simulations. Macro traffic flow focuses on the simulation of traffic load characteristics, such as traffic flow parameters, vehicle type, and vehicle weight, which have been proven to follow a random distribution [10]. Micro-traffic flow focuses on microscopic driving behavior simulations, such as lane change and acceleration. Macro traffic flow is used in bridge dynamic analysis when evaluating the long-term load response of long-span bridges. Wang generated random traffic flow data using the Monte Carlo method (MC) based on the vehicle type, vehicle weight, and axial distance and used it for fatigue analysis of bridges [11]. Zong fitted the distribution of vehicle load parameters based on WIM data and calculated the load effect distribution characteristics of the mixed-state vehicle load model [12]. Han adopted the Monte Carlo method, considering vehicle load parameters, to establish a vehicle bridge coupling analysis framework under random vehicle flow [13]. Usually, a bridge in service can be evaluated by adopting measured traffic flow loading or sampling simulation loading. The measured traffic flow data are authentic but cannot fully reflect the statistical characteristics of long-term traffic flow. Parameter simulation integrating the measured data can better reflect the long-term characteristics of traffic flow [14], wherein the most important aspect of random traffic flow simulation is to reproduce the characteristics of the actual traffic load. Existing research has established diverse traffic load simulation methods. However, the integration effects of long-term traffic flow are constrained because the density of traffic flow varies over the course of days and years, and traffic flow with different intensities is typically divided according to the specific limits in the processing of traffic data. Parameters such as the axle load and wheelbase of vehicles have potential correlation laws because of the connecting shaft. Because economic activities differ significantly across China [15], the load effect of high proportion of multi-axle loaded vehicles in heavy industrial zones cannot be ignored. The generated random traffic flow load effect may deviate from the actual situation if these factors are neglected. A precise traffic flow simulation that considers the difference in daily traffic intensity and parameter correlation should be developed to assess accurately the response and degradation of bridge components.
In this study, a traffic load simulation method for long-span bridges based on highauthenticity WIM data was proposed. The k-means clustering algorithm is improved by introducing Critical Importance Through Intercritical Correlation (CRITIC) weight and sample weight and is used to process the WIM traffic data for the acquisition and simulation of traffic flow intensity characteristics. Subsequently, the correlated variables sampling method based on Sobol sequence and Copula theory (CSSC) was established to conduct the sampling simulation of the parameter correlation of vehicle samples. Finally, a refined stochastic traffic flow was used to analyze the long-span bridge response under the traffic flow of different clusters. This study can provide a reference for the evaluation of long-span suspension bridges under traffic loads.

Method
The traffic load involves the vehicle mass characteristics, vehicle flow speed, traffic volume, and other parameters, particularly for long-span bridges with long loading intervals. Based on the requirement of evaluating long-span bridges under traffic flow load, a traffic flow simulation method was developed in this study by fusing a clustering algorithm with the improved k-means method. The measured traffic data were used to extract the characteristic parameters that are sensitive to the load effect of long-span bridges. The clustering algorithm is used to intelligently divide the traffic flow and assure the applicability of the classification of complex and high-dimensional traffic flow features.

Improved k-Means Clustering Algorithm
The basic concept of the clustering algorithm is to divide the data into multiple heaps, each of which has the same clustering center. k-means clustering is arguably the most popular clustering algorithm, and its principle is to initialize the k cluster centers. Then, the class of samples under each cluster is summed in accordance with the computational distance between the sample and cluster center, as the Euclidean distance shown in Equation (1).
where x it and µ kt denote the coordinates of two points in space. For each sample i, the cluster belongs to what calculate by Equation (2).
where x (i) and µ k denote the coordinates of the sample point i and cluster center k. The cluster center point is updated by the ratio of the sample features to the number of samples, as shown in Equation (3).
where |C l | denotes the number of elements in cluster l, 1 < l < k, and X i denotes the i th element in cluster l. The entire clustering process iterates until the distortion function converges, that is, the sum of squares of the Euclidean distances from all samples to the center of their clusters satisfies the requirements. The minimum value of objective function is set as the convergence condition, that is, the minimum sum of the squared errors, as shown in Equation (4).
where x i is the i th sample point, and µ k is the central point of k th cluster.
In multi-dimensional sample clustering, there are differences in the characteristics of each dimension. For samples with high dimensions, there may be potential correlation laws between certain dimensions [16,17]. Using the traffic flow sample as an example, the traffic volume, traffic flow speed, and traffic flow density conform to classical traffic flow, and there is an apparent correlation between them. Discrete sample points, such as special traffic flow conditions caused by accidents and other emergencies, also interfere with clustering results. The traditional k-means algorithm ignores sample point characteristics and dimensional correlation, resulting in the unsatisfactory clustering effect. Therefore, to improve k-means algorithm, dimensional importance through CRITIC weight and sample weight are introduced.
To eliminate the influence of the order of magnitude and dimension in the data on clustering, the data are first standardized such that each feature parameter has the same scale. Average normalization is used, as expressed in Equation (5) [18]: where Z ij is i th normalized value of the j th sample vector, and x ij is i th sample of the j th vector. The CRITIC weight method determines the weight of indicators through the variability and correlation of sample indicators, and uses correlation coefficients to express the correlation between indicators [19]. For the indicator to be evaluated, the stronger the correlation with other indicators, the less the conflict between the indicator and other indicators, and the more the same information is reflected, such that the weight assigned to this indicator should be reduced. The weight of each dimension is determined using Equation (6).
where z ij is the normalized sample, r ij is the correlation coefficient between evaluation indexes i and j, and m is the dimension of the sample.
To weaken the influence of sample discreteness on the clustering effect, relatively small weights are assigned to discrete sample points. The specific process calculates the distance between each sample point and the surrounding points to reflect the weight of the sample in the entire sample set. The weight of each sample point was calculated using Equation (7).
where δ is an assumed determination coefficient, which can be determined from reference [20]. x j is the sample point and x i is the sample point around x j . Then, the dimension and sample point weight are introduced when calculating the distance between the sample points, as shown in Equation (8).
where µ i is the center of the cluster, w m is the weight of m th dimension, p j is the weight of the sample point j.
The objective minimum square error of the k-means algorithm improved by the information entropy weight and sample weight is modified as follows: where µ i is the center of the cluster, w m is the weight of m th dimension, p j is the weight of the sample point j.
The detailed process of cluster center establishment is: (1) For N sample points, calculate the total distance between each point and other points, and select the point with the largest result as the first cluster center µ 1 ; (2) Calculate the distance from other samples to µ 1 and select the point farthest from the cluster center µ 1 as the second cluster center µ 2 ; (3) Select the center point of the third initial cluster at the point with the largest distance from the first two points, and repeat this step until k initial cluster centers are selected.

Correlated Variables Sampling Based on Sobol Sequence and Copula Function
Multi-axle vehicles typically adopt coupling design of axles, of which the distribution of wheelbase and axle load may be relevant. Typically, the sampling simulation is conducted separately based on the statistical characteristics of each axle, which causes a deviation between the generated sample and the actual vehicle model. In addition, the load effect was not sufficiently precise. Therefore, a simulation method of parameter correlation was established to simulate the traffic load accurately.
Various types of correlation coefficients can reflect the correlation between the parameters [21]. The Pearson and Spearman correlation coefficients are used to determine the correlation of the random variables. The scope of application of the two methods is different. The Pearson correlation coefficient was calculated using Equation (10).
where x i and y i are the elements of variables X and Y, n is the number of variables. The Pearson correlation coefficient ρ can accurately reflect the linear correlation between the random variables. Assuming that ρ L is the Pearson correlation coefficient matrix of the normal random vector L with mean µ L and standard deviation σ L , according to the theory of probability statistics, L can be obtained by linear transformation from the standard normal distribution vector L S as shown in Equation (11). If the standard normal distribution vector L S which has a correlation coefficient matrix ρ L can be obtained, the vector L can also be obtained.
where L S is standard normal distribution random vector with correlation coefficient matrix ρ L . The correlation coefficient matrix ρ L is a symmetric positive definite matrix with diagonal element 1. A lower triangular matrix C and its transposition matrix can be obtained by Cholesky decomposition, and L S can be obtained by multiplying a group of mutually independent standard normal distribution vector Z S and matrix C, as shown in Equation (12): . . .
where C is a lower triangular matrix obtained by the Cholesky decomposition of ρ L , Z S is an independent standard normal distribution vector. Based on the above principles, a sampling simulation of the relevant data can be performed. Each axle load or wheelbase of the vehicle has its own distribution characteristics; therefore, a Copula function is introduced to connect the edge distribution of the axle load [22,23]. According to Sklar theory, if the joint probability distribution function F whose edge distribution satisfies F 1 (u 1 ), F 2 (u 2 ), · · · , F M (u M ), there must be a Copula function that satisfies Equation (13). where C is Copula function, F M is the edge distribution of F. The Copula function C is expressed as Equation (14) C(u 1 , u 2 , · · · , u M ) = F F −1 The Copula function usually includes five functions with different characteristics. Based on the existing research [24], the characteristics of symmetric distribution can be considered by normal Copula function. Combined with the correlation coefficient of the edge distribution, the joint distribution function is expressed by Equation (15).
where R is the Pearson correlation coefficient matrix of the edge distributions, Φ R is the normal Copula function whose correlation coefficient matrix is R, and Φ −1 is the inverse function of the standard normal distribution.
Based on probability integral transformation, the data of a random variable with any given continuous distribution can be converted into a random variable with a standard uniform distribution (see Equation (16)). If X 1 and X 2 are random variables, U 1 and U 2 are their cumulative probability distribution functions, and thus, U 1 and U 2 are subject to a uniform distribution.
Therefore, if the mean, variance, and the Pearson correlation coefficient matrix are known, any set of related normal distribution random variables can be obtained by converting a set of unrelated standard normal distribution random variables. Subsequently, the random number is converted to a random number with uniform distribution on [0, 1], and the correlation coefficient matrix remains unchanged. Finally, through the inverse function of the cumulative distribution function, it can be transformed into a random number sample satisfying the marginal distribution of the original random variable.
The use of the Pearson correlation coefficient has the necessary condition that random variables conform to a continuous normal distribution. However, it is difficult to idealize the actual samples to follow a normal distribution; therefore, the Spearman rank correlation coefficient should be introduced to describe the correlation of samples. The Spearman rank correlation coefficient was calculated using Equation (17).
where d i is the rank difference between variable X and Y, n is the number of variables. Because the Spearman rank correlation coefficient always exists and differs from the Pearson correlation coefficient, when applied to the joint normal distribution random vector (X,Y), the Pearson correlation coefficient can be transformed using Equation (18): Because the random traffic flow simulation process should be sampled based on the characteristics of the parameter distribution, the crucial sampling process is the generation of random numbers. Modern computer-based random number generation algorithms are quasi-random and restricted to a specific cycle. Therefore, when number of cycles is exceeded, it produces repetition and aggregation, resulting in invalid samples. Sobol sequence sampling is a sequence with a uniform distribution than a pseudo-random sequence in a given space, which is characterized by favorable stability [25]. A Sobol sequence with a low difference is a set of sequences formed based on direction number v i . Each direction number is expressed as a decimal value.
where m i is positive and odd numbers less than 2 i . The generation of v i is based on polynomial Equation (20) f = x d + a 1 x n−1 + · · · + a n−1 x + a n (20) where a 1 , a 2 , · · · , a n−1 are the given coefficients.
where indicates rounding down, ⊕ is binary bitwise XOR. The equivalent recursive formula for m is Sobol sequence sampling can be expressed as The parameters n, m i , a can be obtained from the generation matrix C = (d, n, a, m i ), where d represents the dimension of C. Each dimension of the Sobol sequence has a different generation matrix C. Existing research has found that the generation matrix C has the highest dimension of 21,201 [26]. Figure 1 displays the distribution of the two-dimensional sampling results of Monte Carlo sampling and Sobol sequence sampling. The results illustrate that samples of the Sobol sequence are more evenly distributed in the sample space than those of Monte Carlo sampling. sequence with a low difference is a set of sequences formed based on direction number vi. Each direction number is expressed as a decimal value.
where i m is positive and odd numbers less than 2 i .
The generation of vi is based on polynomial Equation (20) are the given coefficients.
where     indicates rounding down, ⊕ is binary bitwise XOR.
The equivalent recursive formula for m is is the binary representation of i. The parameters n, mi, a can be obtained from the generation matrix C = (d, n, a, mi), where d represents the dimension of C. Each dimension of the Sobol sequence has a different generation matrix C. Existing research has found that the generation matrix C has the highest dimension of 21,201 [26]. Therefore, in this study, Sobol sequence sampling is used to generate samples subject to uniform distribution in [0, 1] space and then transformed to a standard normal distribution vector through its inverse function using Equation (24).
is the inverse function of standard normal distribution. Through the function of Pearson correlation coefficients, R, samples subjected to the multivariate standard normal distribution function can be obtained. According to Equation (16), they can be converted into random numbers that obey a uniform distribution on [0, 1], without a change in R. The samples obeying the edge distribution of the Therefore, in this study, Sobol sequence sampling is used to generate samples X i subject to uniform distribution in [0, 1] space and then transformed to a standard normal distribution vector through its inverse function using Equation (24). where Φ −1 is the inverse function of standard normal distribution. Through the function Φ R of Pearson correlation coefficients, R, samples subjected to the multivariate standard normal distribution function can be obtained. According to Equation (16), they can be converted into random numbers that obey a uniform distribution on [0, 1], without a change in R. The samples obeying the edge distribution of the original random variables are then obtained using the inverse function of the cumulative distribution. The detailed steps are shown in Figure 2. original random variables are then obtained using the inverse function of the cumulative distribution. The detailed steps are shown in Figure 2.

Multi-Lane Traffic Load Clustering Model
Vehicle parameters, such as vehicle weight, axle weight, and axle base, and their parameter distribution conform to certain statistical regularities [27]. Based on the established clustering and sampling method, random traffic flow can be simulated according to the macro distribution characteristics obtained by the traffic monitoring system. The sampling process for traffic flow is shown in Figure 3. The traffic parameters were clustered using the improved k-means method to achieve traffic volume, traffic speed, and vehicle type proportion. Then, the axle distance and axle weight of each vehicle are generated by the CSSC method, where the proportion of the hourly traffic volume and vehicle type in each lane obtains the number of each vehicle type. Finally, each vehicle is combined with randomly generated parameter samples based on the statistical characteristics, and the sequence of vehicle samples is randomly sorted. Accurate analysis and evaluation of the traffic flow load effect can be achieved based on the proportion of each cluster in the operation period.

Bridge Information
A long-span suspension bridge was selected as the research subject to analyze the dynamic response under simulated traffic flow. The bridge was a single-span steel-concrete composite girder suspension bridge. The section layout is shown in Figure 4. The span of the main cable was arranged as (250 + 838 + 215) m, and the sagittal span ratio of the midspan main cable was 1/10. The central transverse spacing of the two main cables in midspan is 26.0 m. Two flexible central buckles are set near both sides of the middle span of each main cable to form a cable-beam connection. The stiffening girder adopted a steel-concrete composite structure wherein the steel beam was attached to the concrete deck using shear nails, as shown in Figure 5. The full width of the stiffening girder is 33

Multi-Lane Traffic Load Clustering Model
Vehicle parameters, such as vehicle weight, axle weight, and axle base, and their parameter distribution conform to certain statistical regularities [27]. Based on the established clustering and sampling method, random traffic flow can be simulated according to the macro distribution characteristics obtained by the traffic monitoring system. The sampling process for traffic flow is shown in Figure 3. The traffic parameters were clustered using the improved k-means method to achieve traffic volume, traffic speed, and vehicle type proportion. Then, the axle distance and axle weight of each vehicle are generated by the CSSC method, where the proportion of the hourly traffic volume and vehicle type in each lane obtains the number of each vehicle type. Finally, each vehicle is combined with randomly generated parameter samples based on the statistical characteristics, and the sequence of vehicle samples is randomly sorted. Accurate analysis and evaluation of the traffic flow load effect can be achieved based on the proportion of each cluster in the operation period. original random variables are then obtained using the inverse function of the cumulative distribution. The detailed steps are shown in Figure 2.

Multi-Lane Traffic Load Clustering Model
Vehicle parameters, such as vehicle weight, axle weight, and axle base, and their parameter distribution conform to certain statistical regularities [27]. Based on the established clustering and sampling method, random traffic flow can be simulated according to the macro distribution characteristics obtained by the traffic monitoring system. The sampling process for traffic flow is shown in Figure 3. The traffic parameters were clustered using the improved k-means method to achieve traffic volume, traffic speed, and vehicle type proportion. Then, the axle distance and axle weight of each vehicle are generated by the CSSC method, where the proportion of the hourly traffic volume and vehicle type in each lane obtains the number of each vehicle type. Finally, each vehicle is combined with randomly generated parameter samples based on the statistical characteristics, and the sequence of vehicle samples is randomly sorted. Accurate analysis and evaluation of the traffic flow load effect can be achieved based on the proportion of each cluster in the operation period.

Bridge Information
A long-span suspension bridge was selected as the research subject to analyze the dynamic response under simulated traffic flow. The bridge was a single-span steel-concrete composite girder suspension bridge. The section layout is shown in Figure 4. The span of the main cable was arranged as (250 + 838 + 215) m, and the sagittal span ratio of the midspan main cable was 1/10. The central transverse spacing of the two main cables in midspan is 26.0 m. Two flexible central buckles are set near both sides of the middle span of each main cable to form a cable-beam connection. The stiffening girder adopted a steel-concrete composite structure wherein the steel beam was attached to the concrete deck using shear nails, as shown in Figure 5. The full width of the stiffening girder is 33

Bridge Information
A long-span suspension bridge was selected as the research subject to analyze the dynamic response under simulated traffic flow. The bridge was a single-span steel-concrete composite girder suspension bridge. The section layout is shown in Figure 4. The span of the main cable was arranged as (250 + 838 + 215) m, and the sagittal span ratio of the midspan main cable was 1/10. The central transverse spacing of the two main cables in midspan is 26.0 m. Two flexible central buckles are set near both sides of the middle span of each main cable to form a cable-beam connection. The stiffening girder adopted   A three-girder bridge model was established using ANSYS 18.0 software to si the structural characteristics. The BEAM4 element was used to simulate the main st small stringer, and the pylon. The LINK10 element was used to simulate the cabl ponents. The girder, pylon, and cable components have 1836 BEAM4 elements, 82 B elements, and 279 LINK10 elements, respectively. Only the mass of bridge deck pav was considered as it contributed minimally to the stiffness of the stiffening gird stress stiffening of the suspenders was conducted in accordance with the measured

Traffic Monitoring Data
Traffic flow on bridges is a complicated and changeable load with different densities, volumes, and speeds. These parameters directly influence the respo bridges. In this study, traffic flow data collected using the WIM system was used lyze the bridge's dynamic response under vehicle load. The monitoring data includ hicle data for 31 days. Vehicles were divided into five types according to the num axles, that is, T-1 to T-6. Figures 6 and 7 show the hourly traffic volume and its distri in each lane, respectively. Lane 3 is the passing lane occupied by two-axle ve whereas multi-axle vehicles mostly drive in lane 1; the hourly traffic volume in ea exhibits considerable differences for different traffic speeds. Researchers typical sider traffic density as an indicator to simulate based on long-term overall statist sults [28,29].     A three-girder bridge model was established using ANSYS 18.0 software to si the structural characteristics. The BEAM4 element was used to simulate the main s small stringer, and the pylon. The LINK10 element was used to simulate the cab ponents. The girder, pylon, and cable components have 1836 BEAM4 elements, 82 B elements, and 279 LINK10 elements, respectively. Only the mass of bridge deck pa was considered as it contributed minimally to the stiffness of the stiffening gird stress stiffening of the suspenders was conducted in accordance with the measure

Traffic Monitoring Data
Traffic flow on bridges is a complicated and changeable load with differen densities, volumes, and speeds. These parameters directly influence the respo bridges. In this study, traffic flow data collected using the WIM system was used lyze the bridge's dynamic response under vehicle load. The monitoring data inclu hicle data for 31 days. Vehicles were divided into five types according to the num axles, that is, T-1 to T-6. Figures 6 and 7 show the hourly traffic volume and its distr in each lane, respectively. Lane 3 is the passing lane occupied by two-axle v whereas multi-axle vehicles mostly drive in lane 1; the hourly traffic volume in ea exhibits considerable differences for different traffic speeds. Researchers typical sider traffic density as an indicator to simulate based on long-term overall statist sults [28,29].  A three-girder bridge model was established using ANSYS 18.0 software to simulate the structural characteristics. The BEAM4 element was used to simulate the main stringer, small stringer, and the pylon. The LINK10 element was used to simulate the cable components. The girder, pylon, and cable components have 1836 BEAM4 elements, 82 BEAM4 elements, and 279 LINK10 elements, respectively. Only the mass of bridge deck pavement was considered as it contributed minimally to the stiffness of the stiffening girder. The stress stiffening of the suspenders was conducted in accordance with the measured force.

Traffic Monitoring Data
Traffic flow on bridges is a complicated and changeable load with different traffic densities, volumes, and speeds. These parameters directly influence the response of bridges. In this study, traffic flow data collected using the WIM system was used to analyze the bridge's dynamic response under vehicle load. The monitoring data included vehicle data for 31 days. Vehicles were divided into five types according to the number of axles, that is, T-1 to T-6. Figures 6 and 7 show the hourly traffic volume and its distribution in each lane, respectively. Lane 3 is the passing lane occupied by two-axle vehicles, whereas multi-axle vehicles mostly drive in lane 1; the hourly traffic volume in each lane exhibits considerable differences for different traffic speeds. Researchers typically consider traffic density as an indicator to simulate based on long-term overall statistical results [28,29].   Figure 9 shows the scatter distribution of axle weight and distan cles. It is observed that the axle weight and axle distance of the tw noticeable correlation, whereas the correlation between the front axle    Figure 9 shows the scatter distribution of axle weight and distan cles. It is observed that the axle weight and axle distance of the tw noticeable correlation, whereas the correlation between the front axle    Figure 9 shows the scatter distribution of axle weight and distance for six-axle vehicles. It is observed that the axle weight and axle distance of the two rear axles have a noticeable correlation, whereas the correlation between the front axle and the rear axle is weak.   Figure 9 shows the scatter distribution of axle weight and distance for six-axle vehicles. It is observed that the axle weight and axle distance of the two rear axles have a noticeable correlation, whereas the correlation between the front axle and the rear axle is weak.

Simulated Traffic Flow
Based on the clustering algorithm, six parameters in the traffic flow data of the three lanes, including traffic volume Q and traffic flow speed V, were selected as cluster analysis parameters. Table 1 lists the correlation coefficient for each parameter, the closer it is to 1, the stronger the correlation is. The weight of each dimension is calculated using Equation (6). The quality of the clustering results was measured using the Davies-Bouldin index [30], which tended to be stable when k reached 5. Figure 10 shows the distribution of the clustering results. Figure 11 displays the distribution of the traffic volume, traffic speed, and vehicle proportion in cluster 1. Although the proportion of vehicle types is not clustered, the distribution of vehicle types affects traffic volume and speed. The distribution of traffic parameters in the cluster showed a unimodal distribution after clustering, indicating the same distribution characteristics.

Simulated Traffic Flow
Based on the clustering algorithm, six parameters in the traffic flow data of the three lanes, including traffic volume Q and traffic flow speed V, were selected as cluster analysis parameters. Table 1 lists the correlation coefficient for each parameter, the closer it is to 1, the stronger the correlation is. The weight of each dimension is calculated using Equation (6). The quality of the clustering results was measured using the Davies-Bouldin index [30], which tended to be stable when k reached 5. Figure 10 shows the distribution of the clustering results. Figure 11 displays the distribution of the traffic volume, traffic speed, and vehicle proportion in cluster 1. Although the proportion of vehicle types is not clustered, the distribution of vehicle types affects traffic volume and speed. The distribution of traffic parameters in the cluster showed a unimodal distribution after clustering, indicating the same distribution characteristics.

Simulated Traffic Flow
Based on the clustering algorithm, six parameters in the traffic flow data of the three lanes, including traffic volume Q and traffic flow speed V, were selected as cluster analysis parameters. Table 1 lists the correlation coefficient for each parameter, the closer it is to 1, the stronger the correlation is. The weight of each dimension is calculated using Equation (6). The quality of the clustering results was measured using the Davies-Bouldin index [30], which tended to be stable when k reached 5. Figure 10 shows the distribution of the clustering results. Figure 11 displays the distribution of the traffic volume, traffic speed, and vehicle proportion in cluster 1. Although the proportion of vehicle types is not clustered, the distribution of vehicle types affects traffic volume and speed. The distribution of traffic parameters in the cluster showed a unimodal distribution after clustering, indicating the same distribution characteristics.   Following the clustering analysis, distribution fitting was performed for the characteristics of the vehicle flow parameters in each cluster ( Table 2). Table 2 lists the characteristic parameters of each cluster, which obey the normal, log-normal, and Weibull distributions. Significant differences were observed in the values of the parameters of each cluster, proving the effectiveness of clustering analysis in dealing with the distribution characteristics of traffic flow. The traffic volume of Cluster 5 was the largest, and that of Cluster 2 was the smallest. The parameter distributions of the clusters differ from each other, especially in terms of traffic volume. It is determined by the composition of vehicle types, wherein the proportion of multi-axle vehicles play a decisive role. In Lane 3, which is the passing lane, the share of two-axle vehicles accounts for more than 99%, and the difference in the proportion of vehicle types is neglected. According to the distribution parameters of each cluster, the multi-lane traffic flow can be simulated and used to analyze the bridge's response.  Following the clustering analysis, distribution fitting was performed for the characteristics of the vehicle flow parameters in each cluster ( Table 2). Table 2 lists the characteristic parameters of each cluster, which obey the normal, log-normal, and Weibull distributions. Significant differences were observed in the values of the parameters of each cluster, proving the effectiveness of clustering analysis in dealing with the distribution characteristics of traffic flow. The traffic volume of Cluster 5 was the largest, and that of Cluster 2 was the smallest. The parameter distributions of the clusters differ from each other, especially in terms of traffic volume. It is determined by the composition of vehicle types, wherein the proportion of multi-axle vehicles play a decisive role. In Lane 3, which is the passing lane, the share of two-axle vehicles accounts for more than 99%, and the difference in the proportion of vehicle types is neglected. According to the distribution parameters of each cluster, the multi-lane traffic flow can be simulated and used to analyze the bridge's response.  A clustering algorithm is used to separate the traffic composition into different clusters based on the state parameters to consider the characteristics of the traffic flow. The samples in each cluster have consistent parameter characteristics. Based on the distribution characteristics of traffic flow obtained by cluster analysis, the total number of vehicles, speed of traffic flow, and proportion of vehicle types in each cluster are sampled by the Monte Carlo method, and the distance between vehicles can be obtained by classical traffic flow theory, as shown in Equation (25).
where Q denotes the traffic volume, K denotes the traffic density, and V denotes the traffic speed. Figure 12 shows 3000 axle weight and distance samples of six-axle vehicles generated by the CSSC method, which is similar to Figure 9. This demonstrates that the CSSC sampling method can adequately represent the correlation of the original samples. Q = K·V (25) where Q denotes the traffic volume, K denotes the traffic density, and V denotes the traffic speed. Figure 12 shows 3000 axle weight and distance samples of six-axle vehicles generated by the CSSC method, which is similar to Figure 9. This demonstrates that the CSSC sampling method can adequately represent the correlation of the original samples.

Response under Traffic Load
To analyze the vehicle-bridge interaction, the traffic loads of different clusters were generated based on the proposed method. The process is as follows: (1) collect traffic load characteristics using the WIM system and then cluster them; (2) fit the characteristics of traffic flow parameters based on cluster analysis and generate random traffic flow samples using CSSC sampling; and (3)

Response under Traffic Load
To analyze the vehicle-bridge interaction, the traffic loads of different clusters were generated based on the proposed method. The process is as follows: (1) collect traffic load characteristics using the WIM system and then cluster them; (2) fit the characteristics of traffic flow parameters based on cluster analysis and generate random traffic flow samples using CSSC sampling; and (3) feed traffic flow into the traffic flow bridge coupling analysis system to determine the dynamic response of the bridge. Based on the proposed traffic simulation method, the traffic volumes of different clusters were generated according to their hourly characteristics, and the axle load and wheelbase of vehicles were simulated accordingly. The generated traffic volumes for each cluster were 211, 705, 535, 1071, and 1240. The opposite lane generates a traffic flow to achieve two-way loading. The bridge response under each cluster of traffic flow can be determined by the established vehicle-bridge coupling analysis system, wherein the principle and vehicle parameter dynamic model parameters can be observed in the published literature [31]. Road surface roughness is described by a power spectral density function, that can be generated through Fourier inversion [32]. Figures 13 and 14 show the response of the girder bending moment and vertical displacement at the midspan. To demonstrate the response difference under different cluster traffic flows more intuitively, responses under Clusters 1, 2, and 5 were selected from low to high traffic volumes. The differences in the responses under different clusters can be observed; the girder bending moment under Cluster 2 is the greatest and reaches 28,937.6 kN · m. The maximum vertical displacement of the main beam reached 0.404 m, in Cluster 3. Figure 15 shows the longitudinal displacement response of the girder end. As the amplitude of the longitudinal reciprocating motion of the main beam is small, the difference in the time-history curve is not apparent, and the cumulative displacement is shown in Figure 15b. The growth rate of the cumulative displacement in Cluster 5 was significantly higher than that of the other clusters. This difference may directly influence the fatigue life of bridge components, such as expansion joints. Figures 16 and 17 show the time history response of axial and bending stress variation of a short suspender at midspan. As the bending stress cannot be extracted by the Link 10 element, it was calculated according to the Wyatt formula and experimental results [33,34]. The stress has differences such that peak values appear in different clusters. The maximum bending stress reached 50.676 MPa in Cluster 5, whereas the minimum value of 28.206 MPa appears in Cluster 4. The maximum axial stress reached 67.115 MPa in Cluster 1, which has the smallest traffic volume. It is primarily due to the large dead weight of long-span suspension bridges and the small spacing of suspenders. The axial force on a suspender is primarily affected by the weight of a single vehicle. The bending stress is affected by the traffic load effect.
to high traffic volumes. The differences in the responses under different clusters can be observed; the girder bending moment under Cluster 2 is the greatest and reaches 28,937.6 kN m  . The maximum vertical displacement of the main beam reached 0.404 m, in Cluster 3. Figure 15 shows the longitudinal displacement response of the girder end. As the amplitude of the longitudinal reciprocating motion of the main beam is small, the difference in the time-history curve is not apparent, and the cumulative displacement is shown in Figure 15b. The growth rate of the cumulative displacement in Cluster 5 was significantly higher than that of the other clusters. This difference may directly influence the fatigue life of bridge components, such as expansion joints.
(a) (b)  Figure 15 shows the longitudinal displacement response of the girder end. As the amplitude of the longitudinal reciprocating motion of the main beam is small, the difference in the time-history curve is not apparent, and the cumulative displacement is shown in Figure 15b. The growth rate of the cumulative displacement in Cluster 5 was significantly higher than that of the other clusters. This difference may directly influence the fatigue life of bridge components, such as expansion joints.
(a) (b)    Figures 16 and 17 show the time history response of axial and bending stress variation of a short suspender at midspan. As the bending stress cannot be extracted by the Link 10 element, it was calculated according to the Wyatt formula and experimental results [33,34]. The stress has differences such that peak values appear in different clusters. The maximum bending stress reached 50.676 MPa in Cluster 5, whereas the minimum value of 28.206 MPa appears in Cluster 4. The maximum axial stress reached 67.115 MPa in Cluster 1, which has the smallest traffic volume. It is primarily due to the large dead weight of long-span suspension bridges and the small spacing of suspenders. The axial sults [33,34]. The stress has differences such that peak values appear in different clusters. The maximum bending stress reached 50.676 MPa in Cluster 5, whereas the minimum value of 28.206 MPa appears in Cluster 4. The maximum axial stress reached 67.115 MPa in Cluster 1, which has the smallest traffic volume. It is primarily due to the large dead weight of long-span suspension bridges and the small spacing of suspenders. The axial force on a suspender is primarily affected by the weight of a single vehicle. The bending stress is affected by the traffic load effect.
(a) (b) These results show an extreme response under different traffic flows. To explain the difference in long-term effects under different traffic flows, the variable-amplitude stress was transformed into an equivalent stress. Figure 18 shows the suspender stress under Cluster 1, processed using the rain-flow counting method. When the stress ranges and number of cycles is obtained, the equivalent stress range can be calculated based on the Paris criterion, as in Equation (26).  These results show an extreme response under different traffic flows. To explain the difference in long-term effects under different traffic flows, the variable-amplitude stress was transformed into an equivalent stress. Figure 18 shows the suspender stress under Cluster 1, processed using the rain-flow counting method. When the stress ranges and number of cycles is obtained, the equivalent stress range can be calculated based on the Paris criterion, as in Equation (26).
where ∆σ j and n ij · are the stress range and the corresponding number of cycles, m is the coefficient of Paris law, and has the value 3.
The equivalent stresses of suspender bending and axial stresses are listed in Table 3. The maximum equivalent bending stress and axial stress were 16.068 MPa and 10.542 MPa, whereas the minimum values were 9.429 MPa and 8.679 MPa, respectively. The maximum equivalent stress appears in Cluster 3, and is not consistent with the peak values. Because the stress response is closely related to the fatigue life of suspenders, the response under different traffic flows must be considered in the fatigue analysis of suspenders because of the noticeable difference. The established traffic flow simulation method can consider traffic flow effects at different intensities, which is beneficial for achieving an accurate evaluation in the long-term evaluation of bridge component performance.
where j s  and ij n  are the stress range and the corresponding number of cycles, m is the coefficient of Paris law, and has the value 3. The equivalent stresses of suspender bending and axial stresses are listed in Table 3. The maximum equivalent bending stress and axial stress were 16.068 MPa and 10.542 MPa, whereas the minimum values were 9.429 MPa and 8.679 MPa, respectively. The maximum equivalent stress appears in Cluster 3, and is not consistent with the peak values. Because the stress response is closely related to the fatigue life of suspenders, the

Conclusions
In this study, a random traffic load simulation method for long-span bridges is established based on monitoring traffic data. To analyze the load level composition during service, an improved k-means clustering algorithm was established to process traffic characteristics. Consequently, a CSSC sampling method that considers parameter correlation was established based on Sobol sequence sampling. Finally, the generated stochastic traffic flow was used to analyze the dynamic response of a long-span suspension bridge under the traffic flow of different clusters. The following conclusions were drawn: 1. An improved k-means clustering algorithm was constructed by introducing the CRITIC weight and sample weight to reduce the correlation between sample dimensions and the interference of discrete points. Based on the clustering algorithm, the traffic flow can be separated into clusters with different traffic volumes and speeds. The key parameters of traffic flow and the proportion of vehicle types in each cluster subject to specific distribution characteristics, which are optimized from the multi-peak distribution to the single-peak distribution.
2. The CSSC sampling method was established based on the Sobol sequence and Copula joint function theory. Based on the correlation coefficient matrix of the original sample, the simulation of the vehicle axle load and axle distance considering parameter correlation can be achieved by CSSC. Combined with the clustering results, a refined stochastic traffic flow simulation method is established and used to reconstruct the traffic load. The loading traffic flow can effectively avoid the simple division of traffic intensity and consider the daily traffic characteristics.
3. There are evident differences in the bridge responses under different clusters of traffic flow. Typical indicators, such as the girder bending moment, girder displacement, and suspender stress, show noticeable differences in extreme values under different clusters. The equivalent stresses exhibit clear differences, and the maximum equivalent stress appears in different clusters from that of the peak values. These differences directly influenced the long-term evaluation of these components. The traffic flow density must be considered to achieve an accurate evaluation in the long-term evaluation of bridge component performance.
The multi-lane traffic load model has the capability to model long-term traffic variation trends and characteristic parameters. However, if the bridge is extremely long, aerodynamic forces should be considered. The interaction between traffic loading and aerodynamic loading could further compromise the safety of a bridge; further research and relevant investigations are required in the analysis and evaluation of bridge responses under external loads.  Data Availability Statement: Some or all data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.