Evaluation of passenger comfort with road field test multi-axis vibration

Abstract. Using objective vibration evaluation to produce results highly consistent with real road ride comfort is challenging. The deficiencies in traditional evaluation indices, adopting an average operator, maximum operator, or cumulative operator as the main vibration information integration logic, are reported here through 19 designed road field tests in which major vibration information distribution covers all axes and vibration information is distributed in spacetime in various patterns. A new evaluation index which adopted a combination of maximum and cumulative operator is proposed to overcome these deficiencies and an interactive mechanism which standardized the process of selecting vibration information distributed among axes and spacetime is devised between the localized major vibrations. The results show that the proposed road ride comfort evaluation index is more consistent and accurate than the evaluation indices proposed by ISO 2631-1 and can be used more generally.


Introduction
People in commuter or industrial vehicles are exposed to various types of whole-body vibrations. These vibrations can lead to feelings of discomfort, which are not merely unpleasant, but are also related to health risks [1]. Whole-body vibrations can lead to harmful effects in the cardiovascular, musculoskeletal, and central nervous systems, as well as other areas of the human body [2,3]. Therefore, evaluating these vibrations is crucial, and vibration reduction methods can be developed based on that, such as an optimization design method for a body mounting system developed by Liu et al. [4] and an optimization scheme for the powertrain mount system of an electric vehicle proposed by Xin et al. [5]. The success of those vibration evaluation methodologies rely heavily upon their accuracy and applicability; however, a lack of sufficient theoretical and practical support has led to inaccurate and even contradictory evaluation results compared to true values. This is a problem because recent investigations in the automobile industry revealed that customers are expressing greater needs in terms of improved road ride comfort [6], and a prerequisite to successfully meeting their needs is an appropriate evaluation.
Current road ride comfort evaluation methods were established based on laboratory test results, partly because the test conditions are likely to be strictly controlled, so the results are more comparable. Evidence shows that field tests do not produce independent major vibration information for the longitudinal, lateral, vertical, roll, pitch, and yaw axes, complicating analyses [7,8]. As such, few evaluation methods have been established based on field test results, and this is increasingly becoming a problem. Laboratory tests mostly adopt consecutive sinusoidal excitation as the input, whereas most excitations in actual field tests are non-sinusoidal and non-consecutive -they are random or quasi-random [9,10], so the suitability of directly applying laboratory-based methods to actual field test scenarios is in doubt. Mansfield et al. [11] found that the correlation strength between several vibration evaluation methods and normalized subjective discomfort scores is strongly influenced when excitation type changes from non-random to JOURNAL OF VIBROENGINEERING random or quasi-random. Dickey et al. [12] showed that different combinations of major vibration axes can lead to different levels of discomfort, which is not captured by current methodologies. How actual field test results differ from laboratory based theoretical predictions needs to be investigated in more depth.
Specifically, the road ride comfort evaluation indices can be divided into subjective and objective indices. The subjective evaluation indices require a scoring scale and corresponding semantic levels; a group of people assign their subjective scores after being exposed to vibrations in actual road field tests, then the scores are processed to form a unique value that represents road ride comfort. Subjective indices are generally accurate and are regarded as the true value of road ride comfort because they provide a direct output of human experience. When employing subjective methods, a minimal number of people must be involved to ensure statistical stability. The American Society of Testing Materials (ASTM) specifies a rating panel size of 15 members based on a 0.3 mean panel rating (MPR) maximum error with a normal distribution [13]. Nakamura found that the evaluation results by the American Association of State Highways Officials (AASHO) subjective scoring method deviate from true values by more than 0.4 units in less than 5 % possibility when the number of participants reaches 20 [14]. People who participate in the test must also be representative of the population in terms of age, sex, and other factors so the results will be more representative. These requirements mean subjective methods are time consuming and expensive, and are therefore unable to be scaled up for large-scale use. Subjective methods are a black box process; the quantitative relationship between objective vibration information and subjective scores is unclear, so they cannot directly provide further guidance for health risk prevention.
Objective evaluation indices are more widely used because they involve a white box process, and the data collection and processing are performed more automatically, which lowers the cost. Objective evaluation indices typically integrate vibration information in the entire timeline or frequency zone to represent real road ride comfort; the integration methods vary and lead to different results. First, an average operator is adopted as the main logic to integrate vibration information. The international roughness index (IRI) is defined as the cumulative displacement of the suspension system when a standard vehicle travels 1 km at a speed of 80 km h −1 [15]; it is transformed to create ride quality indexes (RQI), such as China RQI [16] and Minnesota RQI [17]. IRI has its own deficiencies: (1) it is based on displacement, whereas ride comfort is more related to acceleration or other physical quantities [18][19][20]; and (2) the human body is sensitive to wavelengths between 0.5 and 2.5 m [21]. IRI is sensitive to wavelengths between 2.4 and 15 m [22]; therefore, it cannot appropriately reflect human response to vibration frequencies. IRI normalizes the integration interval to 1 km irrespective of the real interval, so its integration mechanism uses an average operator. Root mean square (RMS) [23] values of frequency weighted acceleration and average absorbed power (AAP) [24] also adopt this average logic. Griffin [25] suggested that using an average operator in the time zone is oversimplistic because it can result in excessive estimation for short durations or apparently unreasonably low estimation for long durations; periods containing different magnitudes, transients, and shocks in varying vibration exposure must be defined. Standardizing the total integration time can solve this problem, but it is almost impossible to achieve [10,26]. Secondly, a maximum operator is adopted as the main logic to integrate vibration information. The maximum transient vibration value (MTVV) is a typical example; it only integrates the largest vibration information within one second [23]. However, the vibration information beyond the maximum one second part can also significantly affect human discomfort, so MTVV is unsuitable [2,11]. Vibration greatness (VG), proposed by Miwa [27], integrates vibration information beyond the maximum part, but the integration process is purely in the frequency domain, which is inevitably accompanied by a loss of information on location and geometric characteristics for transients that are essential for evaluating road ride comfort [28]. Thirdly, a cumulative operator is adopted as the main logic to integrate vibration information. Indices using this logic only increase when new vibration information is introduced and thus do not suffer from large shifts due to total integration time alone. Vibration dose value 3 (VDV) is a typical example and is currently considered a good index because its high integration order of four is deemed to appropriately consider shocks, which is based on frequency weighted acceleration capable of representing human response to frequencies [11,29]. Griffin [26] stated that there is currently no evidence to suggest a measure that is more accurate than VDV. Total absorbed power [30] and total daily acceleration dose [31] also adopt a cumulative operator as the main vibration information integration logic. Total absorbed power does not need frequency weightings because it is energy-based. The total daily acceleration dose only considers vibration in the vertical axis and can only provide guidance for spinal risk prevention. Indices using a cumulative operator to integrate vibration information in the time domain are at risk of over integration. For example, adults may suffer considerable vibrations during their lifetime when riding in vehicles. Although this means that the calculation results for these indices are high, the passenger may not currently feel uncomfortable. The International Standard Organization (ISO) specifies that the VDV and total daily acceleration dose should be limited in a whole day, which suggests that the longer the vibration exposure, the more inaccurate the two indices. Fourth, a combination of different operators is adopted as the main logic to integrate vibration information. Weighted longitudinal profile (WLP) uses a combination of variance and range to evaluate road ride comfort [28]; however, it cannot completely avoid the deficiencies of both operators. One variance can correspond to infinite combinations of wavelengths and amplitudes, which can lead to diverse levels of road ride comfort [32]. This is also the case even if a peak value is fixed. As such, determining how vibration information should be integrated to better represent real road ride comfort is a problem to be overcome.
In this study, we aimed to establish a more general objective road ride comfort evaluation index. When applied to various road field test scenarios, the proposed index is consistent and accurate. This process is based on field test results, which are deemed to more accurately represent actual situations. The designed scenarios cover major vibration information in all axes because the current trend is to use multi-axis vibration information [33,34]. The main findings are that localized major vibrations and their interactions, rather than all vibration information, should be integrated to appropriately represent road ride comfort, and the interactive forms vary according to localized characteristics.

Field test configurations
A subjective scoring method was developed in this research. Considering the accuracy required in the analysis and the operational ease, the scoring scale was designed as shown in Fig. 1. Only 4 semantic levels are involved, and the minimal interval is 0.1. The objective vibration data were collected by JY-901-type or HWT-901B-type posture sensors (Witmotion, Shenzhen, China). They collect the same vibration information; the difference is that the data transmission line is fixed to the latter, thus lowering the risk of connection failure. The posture sensors collect translational axis acceleration (m s −2 ), rotational angular velocity (rad s −1 ), and rotational angles (rad) in the , , and directions at a frequency of 200 Hz. The sensors were attached to the central surface of the seat of a Baojun-730-type test vehicle (SACI-GM, Shanghai, China) using Fixate gel pads (Frosted Concepts, Bristol, Australia). The test vehicle and its layout are shown in Fig. 2. The vibration information collected by the posture sensor in the right front seat also represents the information in the right rear seat. The sensors were inertial type and set onto an inclined surface, and adopted a Tait-Bryan coordinate system. The actual directly used vibration information was calculated using Eq. (1): where ( ), ( ), ( ), ( ), ( ), and ( ) are the actual used accelerations in the , , and axis and rotational angular velocities in the , , and axis in moment , respectively. The above variables with a subscript 0 directly represent the corresponding vibration information output directly by the posture sensors. Δ , Δ , and Δ are the rotational angles in , , and axis of the Tait-Bryant coordinate system, respectively; is the local gravitational acceleration. The installation ensures Δ ≡ 0, whereas Δ , Δ , and g were chosen as the recorded maximum density value of the angle in the and axes, output directly by the posture sensors, and acceleration in the axis after the treatment in the right of Eq. (1), respectively. The test scenarios were divided into 4 groups. The scenarios in the first, second, and third group were deemed as producing major translational vibration information only in the , , or axis, respectively. Scenarios in the fourth group were deemed as producing major vibration information in more than one translational axis. Vibration information productions in the rotational axes were not especially designed, they occur with vibration information in the translational axes. Fig. 3 shows ( ) , ( ) , ( ) , ( ) , ( ) , and ( ) in the designed scenarios. The first group is shown in Fig. 3(a-j). The speed bumps in scenarios represented by Fig. 3(a-c) were trapezium-shaped, 110 mm in the upper base, 300 mm in the lower base, 45 mm high, and 3300 mm wide. Each FRP slab in scenarios represented by Fig. 3(d-h) was 300 mm in length, 55 mm high, and 1000 mm wide. Two FRP slabs formed a pair that was crossed by the test vehicle simultaneously. Fig. 3(i, j) show the scenarios of crossing an abnormally raised manhole, only wheels on one side crossed the manhole, this side was called the directly crossing side, the other side was called the indirectly crossing side. The second group is shown in Fig. 3(k-m), the third group is shown in Fig. 3(n, o), and the fourth group is shown in Fig. 3(p-s). The schemes in Fig. 3 (n, o, q-s) are shown in Fig. 4. In these designed scenarios, our major considerations were to include as many excitation types and major vibrational axes and vibration information distribution characteristics as possible. ISO 2631-1 [23] or BS (British Standard) 6841 [35] adopt the crest factor, defined as peak value divided by RMS, as a judgement index for choosing a specific ride comfort evaluation index, suggesting that the vibration information distribution characteristics essentially affect the performances of the indexes. Therefore, the spacing between excitations and the number of excitations is considered. Note that the vibration information shown in Fig. 3 is only a representative; it might not be reproduced precisely because the movements were conducted manually by the drivers. Vibration information on only one side is provided because symmetry characteristics were designed; the exceptions are the scenarios represented by Fig. 3(i, j), where they are highly asymmetric. randomly assigned and the testers were not told what excitation would they experience to avoid the Hawthorne effect [36]. All subjects provided their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of Tongji University (2019tjdx294).

Analysis of the road field test results
RMS, MTVV, VDV, and pseudo vibration total value (PVTV) were chosen as the reference road ride comfort evaluation indexes for the road field test. The calculation equations for RMS, MTVV, and VDV are shown in Eqs. (2)(3)(4), respectively [23]. The calculation equation for PVTV is shown in Eq. (5). The frequency weightings for these indexes are provided by ISO 2631-1, in which and axis, axis, and rotational axes respectively adopt a distinct frequency weightings. We chose these indexes because RMS, MTVV, VDV, and vibration total value (VTV) are currently widely used. RMS, MTVV, and VDV adopt an average operator, a maximum operator, and a cumulative operator as the main logic to integrate vibration information, respectively. VTV is a further development for RMS, which includes vibration information in the rotational axes, whereas RMS, MTVV, and VDV do not. These indexes are contradictory on this point but are included in the same ride comfort evaluation standard ISO 2631-1, this is another reason for why we put these indexes together for comparisons.
PVTV is different from VTV on the point that in the rotational axes it is based on frequency weighted rotational angular velocity rather than frequency weighted rotational angular acceleration. Using PVTV rather than VTV is because the posture sensors directly collect angular velocities, calculating angular accelerations directly from discrete angular velocities would lead to excessive large results. The alternatives for absolute balance weight in the rotational axis is chosen as 1 to 10 with an interval of 1, these alternatives do not allow the rotational axis contribution to exceed 50 %, as is suggested by ISO 2631-1. The optimal for PVTV is 1, and this parameter set PVTV is chosen for comparisons. The optimize model is chosen as recursive partitioning and regression trees (RPART) [37], the input is PVTV and testers, the output is subjective scorings. Testers were added to the input variables because different people respond largely to the same kind of vibration. Adding this variable can help to largely reduce the fitting errors [38]. For the RPART model, the minimum number of observations in a node for a split to be attempted was set to 20, any split that does not decrease the overall lack of fit by a factor of complexity parameter of 0.015 was not attempted, the maximum depth of any node of the final tree was set to 10 (the root node was counted as depth 0), and the number of cross-validations was set to 10. The normalized mean square error (NMSE) was chosen as the measurement of fitting errors, which is the main measurement of the performance of the objective evaluation indexes. The reference model for NMSE is the mean value model: where ( ) is the momentum in axis at moment , which is frequency weighted acceleration in the translational axes or frequency weighted angular velocity in the rotational axes; is the total integration time; is the local integration time, its default value is 1 s; is any end moment; is the axial weight, where is the absolute balance weight in axis j, which is used to balance the weights between translational axes and rotational axes; we set = = = 1 for standardizing purposes. = = = , is the absolute balance weight in the rotational axis, > 0. ISO only suggests the calculation for RMS, MTVV, and VDV in the translational axes, which we expanded to the rotational axes to enable comparison in this research.
The subjective scoring results and objective evaluation results for groups 1, 2, 3, and 4 are shown in Fig       S e r i a l 1 2 x) y) z) Fig. 5. Boxplot of subjective scoring results and objective evaluation results in the designed scenarios: a), f), k), p) the subjective scoring results for groups 1, 2, 3, and 4, respectively; b)-e) RMSz, MTVVz, VDVz, and PVTV for group 1, respectively; g)-j) RMSx, MTVVx, VDVx, and PVTV for group 2, respectively; l)-o) RMSy, MTVVy, VDVy, and PVTV for group 3, respectively; q)-s) RMSx, RMSy, RMSz for group 4, respectively; t)-v) MTVVx, MTVVy, MTVVz for group 4, respectively; w)-y) VDVx, VDVy, VDVz for group 4, respectively; and z) PVTV for group 4 The t-test results are described as: (1) when ≥ 0.05, there is no significant difference; (2) when 0.001 ≤ < 0.05, there is significant difference; (3) when < 0.001, there is highly significant difference. Wasserstein and Lazar [39], on behalf of the American Statistical Association (ASA), stated that traditional thresholds for p-values are empirical and should not be strictly followed. In this case, when subjective scoring results and objective evaluation results are highly contradictory, we must comment; if we do not provide any comment, we use p-values as an auxiliary judging measure. The column sequence number and -values and -values in Fig. 5 are marked red when this column is suspected to have a major problem. The values in group 4, except for those in Fig. 5(z), are all marked green because the corresponding scenarios involve more than one major vibrational axis; therefore, using vibration information in only one axis to judge the total discomfort would be inappropriate. RMS, MTVV, VDV, and PVTV performed poorly in some of the designed scenarios. For group 1, when compared to Fig. 5(a): (1) Fig. 5(b) shows RMS highly significantly over-integrated the vibration information in the scenarios represented by columns 2 and 3, and highly significantly under-integrated vibration information in scenarios represented by columns 7 and 10; (2) Fig. 5(c) shows MTVV produced evaluation results consistent with the subjective evaluation results except for highly significantly under-integrating the vibration information in the scenarios represented by column 7; (3) Fig. 5(d) shows VDV highly significantly over-integrated the vibration results in the scenarios represented by column 2 and 3, and highly significantly under-integrated vibration information in scenarios represented by columns 7, 8, and 10; (4) the -value and -value in the additional -test between columns 1 and 6 in Fig. 5(e) are 4.947 and 5.64×10 −6 , respectively. Therefore, Fig. 5(e) shows PVTV highly significantly over-integrated the vibration information in scenarios represented by columns 1-3, and highly significantly under-integrated the vibration information in the scenarios represented by columns 7 and 10.
For group 2, the following contradictory evaluation results were obtained: (1) the subjective scorings in column 3 are significantly higher than in column 2, whereas the RMS x in column 3 is significantly lower than the RMS x in column 2; (2) the subjective scorings in column 3 are highly significantly higher than the subjective scorings in column 1, whereas the MTVV x in column 3 is highly significantly lower than the MTVV x in column 1 ( = −5.058, = 6.356 × 10 −6 ); (3) the subjective scorings in column 3 are significantly higher than subjective scorings in column 2, whereas VDV x in column 3 is highly significantly lower than VDV x in column 2; (4) the subjective scorings in column 3 are significantly higher than subjective scorings in column 2, whereas there is no significant difference between the PVTV in columns 2 and 3.
For group 3, when compared to Fig. 5(k): (1) MTVV highly significantly under-integrated the vibration information in the scenario represented by column 2; and (2) VDV highly significantly over-integrated the vibration information in the scenario represented by column 1.
For group 4, the RMS y in column 4 is highly significantly higher than that in column 3, whereas we found no significant difference in MTVV y or VDV y between columns 3 and 4. The mechanism into these discrepancies must be interpreted, and a new road ride comfort evaluation index capable of overcoming these discrepancies needs to be developed to achieve better consistency and accuracy.

Establishing a new ride comfort evaluation index
Columns 1-3 in Fig. 5(a-e) provide the vibration information introduced by only one speed bump integrated to represent real discomfort; otherwise, over-integration occurs. So, the localized vibration information introduced by each speed bump was extracted first and then we decided which information was going to be integrated. Localized characteristics are essential for road ride comfort evaluation because vibration information distribution in an actual road riding process is often highly uneven [40,41]. A local impact concept is proposed in the following steps to capture the major vibration information in localized segments: 1. Calculate translational axes energy density ( ) using: 2. Find a consecutive interval that covers the largest ( ) point in the residual timeline, can maximally integrate total ( ), and is limited to a maximum length of local integration time τ. This procedure is shown in Eqs. (7) and (8): where , is the total translational axes energy for the th local impact; is the end moment for the th local impact, < 0; and , is the moment when ( ) in the residual timeline reaches its maximal value.
3. Remove the vibration information included in , from the existing timeline.
6. Rename the left , from its emerging sequence to ensure the sequence starts from 1. With each time step, add 1, then the slice corresponding to , is defined as the th local impact.
Steps 4 and 5 above were designed according to the Beibo law: when someone experiences a strong excitation, then other excitations become negligible [42].
We found no significant difference between columns 4 and 5 in Fig. 5(a), which means the number of FRP slab pairs in these scenarios do not affect the total discomfort and is accordance with the results in columns 1-3 in Fig. 5(a, c). However, when the time spacing between FRP slab pairs is sharply reduced, as shown in Fig. 3(e, f), the total discomfort highly significantly increased, as shown in columns 5 and 9 in Fig. 5(a). So, localized major vibrations do not always individually represent total discomfort; they can interact with each other and add to total discomfort when their time spacings are short enough. In this case, a vibration perception window concept is proposed to quantify the conditional interactions between localized major vibrations, which assumes localized vibrations can interact with other vibrations within a limited time window. Columns 1 and 4 in Fig. 5(d) show that the vibration intensity produced by the FRP slab pair was higher than that produced by the speed bump. Fig. 3(c, f) show the time spacings between the FRP slab pairs were larger than those between the speed bumps, yet the FRP slab pairs were able to interact with each other and highly significantly add to total discomfort, whereas the speed bumps were not, as shown in columns 1 and 2 in Fig. 5(a), which suggests the FRP slab pair was able to interact with a wider range of vibration information than the speed bump. In this case, the window length was designed to increase with increasing localized vibration intensity, as shown in Eq. (9). Fig. 3(e, l) each contain the same two designed types of excitations the time spacings between the excitations in the two scenarios were designed to be the same. The vibration intensity represented by VDV in the main vibration axis was also not significantly different, as shown in column 5 in Fig. 5(d) and column 2 in Fig. 5(i). The -and -values were 0.235, 0.815, respectively. This led to different responses in the subjective scorings. The number of FRP slab pairs in Fig. 3(e) did not significantly add to total discomfort, as shown in columns 4 and 5 in Fig. 5(a), whereas the number of sudden stops in Fig. 3(l) highly significantly added to total discomfort, as shown in columns 1 and 2 in Fig. 5(f), which suggested different influence energy weights among translational axes, as expressed by: where is the window length for the th local impact; and are the scale parameters, > 0 and > 0; , is the total influence energy for the th local impact, its integration interval is the same as , ; ( ) is the influence energy density at moment ; , , and are the axial influence energy weight in the , , and axis, respectively. > , > 0, > 0, = 1 for standardizing purposes. We found no significant difference between columns 8 and 9 in Fig. 5(a), which indicated the interaction intensity did not sharply decrease with the increase in time spacing within a certain range, so we introduced the bathtub curve proposed by Yu [43], as shown in Eq. (11), to characterize the window shape: where ℎ( ) is the bathtub curve function; is the failure rate parameter, > 0; is the interval parameter, = 1 for standardizing purposes; and and are the shape parameters. We set = 0.2 and = 0.1 for symmetry purposes. Columns 8-10 in Fig. 5(a-e) show that only MTVV as able to show that the high concentration of vibration information alone can highly significantly add to total discomfort; therefore, a normal distribution density curve was added to the center of the vibration perception window to assign more weight to the vibration information close to the integration center. This is modeled by: , 0 , where ( ) is the vibration perception window for local impact ; is the intercept item, 0 < < 1 + 0.559 . This range of can ensure the window height does not exceed 1; a higher indicates more vibration information beyond the reference local impact is integrated; is the variance for the normal density curve; ′ is the moment when the bathtub curve as is expressed in Eq. (11) reaches its lowest point; and is the rectify coefficient, which is introduced for standardizing purposes and, combined with the range of , can ensure the window height is always equal to 1. The quantitative interactions between local impacts is shown in Eq. (15). A local impact integrates weighted vibration information in other local impacts through its own vibration perception window. The weights are the values of the vibration perception window that correspond exactly to the central location of other local impacts: where , , is the vibration information integrated by local impact from local impact on the axis, ≠ ; is the basic integration order, > 0. a) b) c) Fig. 6. Boxplot of VDV in the rotational axis for group 1: a) -axis, b) -axis, and c) -axis Columns 4-6 in Fig. 5(a, c, d) show the values in column 6 are highly significantly underestimated. Columns 7 and 8 in Fig. 5(a, b, d, e) show the values in column 7 are highly significantly underestimated. The underestimation in columns 6 and 7 in group 1 occurs under the condition that the z axis vibration information in group 1 was found to be dominant in the translational axes. Therefore, the rotational axis VDV in group 1 was plotted (Fig. 6) to verify if the rotational vibration information in columns 6 and 7 of group 1 was significantly different from that in the other columns. Fig. 6 shows that columns 6 and 7 in group 1 contain much higher vibration information in the rotational -axis than the other columns in group 1. Columns 6 and 7 in group 1 also contain much more vibration information in the rotational -axis compared with the scenario that included only one FRP slab pair (column 4 in group 1), and all the columns in group 1 contain minor vibration information in the rotational -axis. This result showed that the vibration information in the rotational axis is needed to appropriately evaluate actual ride comfort. As such, a characteristic overall vibration value (COVV) was established by Eq. (16), which involves multi-axis vibration information, adopts the Beibo law in that it only uses vibration information in the local impacts, and selects only one reference local impact to avoid repeated integration: where we set to 1 for standardizing purposes.

Verification and revision of the COVV model
RPART model is used to build the connection between objective vibration information and subjective scorings. The parameter settings are the same as that for the optimization of PVTV. The testers are also added to the input variables. The response variable is the subjective scorings. RMS, MTVV, and VDV do not provide a synthetic value for all axes, so for each of these indexes, their values in the , , , rotational , rotational , and rotational axis are simultaneously put as the input variables. For COVV, the input variables are directly themselves.
The L 27 (3 13 ) orthogonal test table is used to optimize the parameters in the COVV model. The parameters to be optimized and their alternative values are shown in Table 1, the initial alternative values are chosen empirically. The reason why and are added to the parameters is because different and even the same standard has varied specification for these parameters.  Table 1 are shown in Table 2. The best parameter set for COVV in Table 2 is the 22th parameter set, the NMSE reaches the lowest value to 0.340, the parameter set is = 1.  Fig. 7. Fig. 7 has the same serial number as Fig. 8, the values between the below and the neighboring right column are always lower or equal to 0 in Fig. 8. Fig. 8 includes all the designed scenarios. Δ in Fig. 8 is defined as the difference between the mean values in the below column and the neighboring right column divided by the mean value in the neighboring right column, the subscription 1 means the values are concerned with the subjective scorings, and if the subscription is 2, it means the values are concerned with the objective evaluation indexes. |Δ − Δ | is another auxiliary judging measure for the performances of the objective evaluation indexes. It describes the discrepancy between subjective scoring results and objective evaluation results. It's ideal value is 0. There are obviously three major problems in Fig. 7: (1) there is no significant difference between column 10 and 14 ( = -0.585, = 0.562), which is contradictory to the results in Fig. 8 where the values in column 10 is significantly lower than the values in column 14; (2) the values in column 14 should be highly significantly higher than the values in column 11, as shown in Fig. 8, whereas the evaluation results is totally on the contrary ( = -3.8498, = 3.388×10 -4 ); (3) too much vibration information is integrated in scenarios represented by column 13, 15, 17, 19, especially when compared to column 18 which also includes the heavy sway excitation.

JOURNAL OF VIBROENGINEERING
To solve the first and the second major problems in Fig. 9 at the same time, we can assign more weight to the non reference local impacts when there is more axis vibration information in the reference local impact, which means to increase the intercept item .  (27) 0.418 COVV(*) means the *th parameter set in the L27(3 13 ) orthogonal test in Table 1 In this way, relatively more vibration information is integrated in column 14 rather than column 11 in Fig. 7 because the major vibrational axis in scenarios represented by column 14, 11 in Fig. 7 are , axis, respectively. Relatively more vibration information is integrated in column 14 rather than column 10 in Fig. 7 because the excitation numbers in scenario represented by column 14 in Fig. 7 is higher, as shown in Fig. 4(l), (m). To partly solve the third major problem in Fig. 9, we can assign more weight to the non reference local impacts when there is more y axis vibration information in the reference local impact. In this way, relatively more vibration information is integrated in column 18 rather than column 13, 15, 17, 19 in Fig. 7 because there are more major vibration durations beyond the reference local impact in scenario represented by column 18 in Fig. 7, as shown in Fig. 4(n), (o), (q), (r), (s). More consistency will also be achieved between column 13 and 15 in Fig. 7, there is no significant difference between the two columns ( = -0.019, = 0.985), which is contradictory to the results in Fig. 8 where the values in column 13 is significantly lower than the values in column 15 ( = -3.055, = 0.004), relatively more vibration information would be integrated to column 15 because the corresponding scenario includes one more FRP slab pair excitation. Note that we cannot assign more weight to the non reference local impacts when there is more axis vibration information in the reference local impact, because it would lead to highly contradictory results with Fig. 8 when compare column 11 and 12 in Fig. 7 to column 8-10, 14 in Fig. 7, much relatively more vibration information would be integrated to column 11 and 12 in Fig. 7. According to the above analysis, the in Eq. (12) should be replaced by ′ in Eq. (17). In order to further solve the third major problem in Fig. 7, we can lower the , because the translational axes contribution ratios in column 13, 15, 17, 19 in Fig. 7 are 0.761, 0.747, 0.743, 0.754, respectively, suggesting too much vibration information in the axis is integrated: where ′ is the rectified intercept item.
Combining the analysis between Fig. 7 and 10 and range analysis results for Table 1, the alternatives in Table 1 are changed into the alternatives in Table 3. Note that the range analysis is only an auxiliary judging measurement, because it alone can result in local optimum results. When the in Eq. (12) is replaced by ′ in Eq. (17), the NMSEs for COVV in parameter set shown in Table 3 are shown in Table 4. The best parameter set for COVV in Table 3 is the 4th parameter set, the NMSE is 0.354, the parameter set is = 0.8, = 0.2, = 0.6, = 1.6, = 1.6, = 0.4, = 0.01, = 0.6, = 0.8, = 0.6, = 15, = 2.1, = 1.9, the boxplot for the best parameter set COVV is shown in Fig. 9, the column sequence in this boxplot is in accordance with that in Fig. 8.   a) b) Fig. 7. Boxplot of COVV in the 22th parameter set in Table 2 for all the designed scenarios: a) the results for scenario 1-10; b) the results for scenario 11-19 a) b) Fig. 8. Boxplot of subjective scorings for all the designed scenarios: a) the results for scenario 1-10; b) the results for scenario [11][12][13][14][15][16][17][18][19] value and value between column 10 and 14 in Fig. 9 is -0.877, 0.385, respectively, it means increasing for the sudden halt excitation is not enough, more duration rather than peak value should be stressed for this excitation type, therefore the integration order should be lowered. This is the same for the heavy sway excitation type, because the comparison between column 17 and 18 in Fig. 9 show that the ride comfort in the scenario containing 3 consecutive heavy sways is highly underestimated on condition that the vibration perception window already gives enough weight to the vibration information beyond the reference integration location, more duration rather than peak value should be stressed, the peak value comparisons for scenarios including at least a heavy sway are shown in Fig. 5(m) and (u). Note that there is no evidence in Fig. 9 to show that the current integration order for the axis has a major problem. In this case, the results suggest different integration order for different translational axes, so Eq. (15)-(17) should be rectified into Eq. (18)- (20), respectively: where is the integration order for axis , ≠ , ≠ , = = = .  Table 3 The rotational axis contribution ratios in column 13, 15, 17, 19 in Fig. 9 are 0.532, 0.491, 0.473, 0.478, respectively, they are too high and are contradictory to the suggestion in ISO 2631-1, BS 6841, and many other ride comfort evaluation indexes in that the major contribution should come from the translational axes. If the integrated vibration information in the rotational axes should remains the same, then Fig. 9 suggest a even lower to satisfy consistency, then the contribution proportion in the rotational axis would be even higher. The rotational axis contribution ratios in column 7, 8 in Fig. 9 are 0.408, 0.236, respectively, and as is discussed before, these two columns especially need to integrate vibration information in the rotational axes. So the above results suggest vibration information integration in the rotational axis should be suppressed only when their values are high. In this case, a suppression order is introduced, and Eq. (19) should be replaced by Eq. (21): where is the suppression order for axis , = = = 1, which means there is no suppression for the translational axes, = = = , is the suppression order in the rotational axes, 0 < < 1. a) b) Fig. 9. Boxplot of COVV in the 4th parameter set in Table 4 for all the designed scenarios: a) the results for scenario 1-10; b) the results for scenario [11][12][13][14][15][16][17][18][19] The final flowchart for the calculation of COVV is shown in Fig. 10, which is derived from the calculation processes as shown in section 2.3 and 2.4. The red arrows display the path of the main flow, the black arrows display the path of the subflow, the red frame display the final step and final output.
The scenario as shown in Fig. 3(l) is taken as an example for the calculation of COVV under the flowchart given by Fig. 10. The 25th parameter set in the L 64 (4 21 ) orthogonal test in Table 5 is adopted because it is the best parameter set for the final rectified COVV model. When step 1 and 2 in Fig. 10 is finished, the locations and the corresponding interaction windows of the calculated a total of 2 local impacts are shown in Fig. 12 , the maximum vibration information is integrated when = 1, as is visualized by the red mark in Fig. 12(l), the COVV is calculated as 5.171 with the aforementioned calculation results and the 25th parameter set in the L 64 (4 21 ) orthogonal test in Table 5.
1. Find out all the local impacts following the steps described in section 2.3 2. Calculate the interaction window w i (t) for local impact i using Equation (12) 2.1 In Equation (12), the end moment for the i th local impact t i is defined by the calculation process of local impacts, the moment when the bathtub curve reaches its lowest point t' is calculated by Equation (13), the interactive window length for the i th local impact L i is calculated by Equation (9,10), the bathtub curve function h(t) is calculated by Equation (11), the rectify coefficient R is calculated by Equation (14). 3. Calculate the vibration information integrated by local impact i from local impact m on the j axis using Equation (18) 2.2 In Equation (12), the intercept item c need to be replaced by the rectified intercept item c' as is calculated by Equation (20). 4. For each local impact i, integrate its vibration information both in the translational and rotational axes and additionally integrate the vibration information it gained from its interaction with other local impacts. When the integrated amount reaches the maximum, output as COVV. This whole step is expressed by Equation (21).

Results for the final rectified COVV model
Given the revision in section 2.4, the comparing results between Fig. 8 and 9, the range analysis in Table 4, the final spreadsheet come to Table 5 and the results are shown in Table 6. Note that the orthogonal test changes from L 27 (3 13 ) to L 64 (4 21 ) because more parameters are involved. The best parameter set for COVV in Table 6 is the 25th parameter set, the NMSE is 0.325, the parameter set is  Table 6 is shown in Fig. 11, the column sequence in the boxplot is the same as Fig. 8. The range of |Δ − Δ | between the results in Fig. 8 and 11 is [0.004, 0.448], on average is 0.118. The visualizations for the interactions between the local impacts in all the designed scenarios when parameter set for COVV is the 25th in Table 6 are shown in Fig. 12. The width between the two vertical lines in the central area of a vibration perception window in Fig. 12 is the local integration time . a) b) Fig. 11. Boxplot of COVV in the 25th parameter set in Table 6  COVV 1 + 0 + 0 + 0 + 0 9 + 0 + 0 + 0 + 0 3 + 0 + 0 + 0 + 0 0 + 1 + 0 + 0 + 0 0 + 0 + 0 + 1 + 0 0 + 2 + 0 + 0 + 0 0 + 0 + 1 + 0 + 0 + in d ir e c t 0 + 0 + 1 + 0 + 0 + d ir e c t 0 + 1 + 0 + 1 + 0 + is o la ti o n 0 + 0 + 0 + 2 + 0  Fig. 12. Visualizations for the interactions between the local impacts in all the designed scenarios when parameter set for COVV is the 25th in Table 6: a)-s) correspond to Fig. 4(a-s), respectively   Table 5 4. Discussion   Fig. 11 is generally in consistent with Fig. 8, which indicate the 25th parameter set COVV in Table 6 has good performances. The NMSE for the 25th parameter set COVV in Table 6 is 0.325 and is 20.538 % lower than that of RMS which is the best evaluation index given by ISO 2631-1 for the designed scenarios. PVTV performs worst among the reference indexes, its NMSE is 68.615 % higher than that of the 25th parameter set COVV in Table 6 and 33.985 % higher than that of RMS. VDV, MTVV respectively give only1.467 %, 7.579 % higher NMSE than RMS, and according to their varied performances, they cannot completely replace or be replaced by RMS.
The range and average |Δ − Δ | between Fig. 8 and 11 is also generally acceptable. |Δ − Δ | between Fig. 8 and 11 reaches its highest level of 0.448 between column 2 and 3, but it is hard to adjust because if the is lowered to ensure less vibration information integration in column 2, then vibration information integration in column 12 would be not enough, as shown in Fig. 12(c) and (f). The |Δ − Δ | between column 2 and 3 in Fig. 5(b), (d) are 0.485, 0.506, respectively, which indicate RMS and VDV performs even worse when compared to the 25th parameter set COVV in Table 6 in scenarios represented by these two columns. But 0.448 is still too high and suggest the method to define be further optimized in future research. Nevertheless, this is still able to show that there is no significant difference between column 1 and 2 in Fig. 11, the |Δ − Δ | is only 0.104, which is not the case for RMS, VDV, and PVTV, the |Δ − Δ | between column 1 and 2 in Fig. 5(b), (d), (e) are 0.356, 0.346, 0.254, respectively, they are too high and show that both RMS, VDV, and PVTV can give highly significantly contradictory results between the scenarios represented by Fig. 4(a) and (c). In Fig. 11, the values in column 16 is highly significantly higher than the values in column 11 and 12, which shows the 25th parameter set COVV in Table 6 gives consistent evaluation results with the results in Fig. 8, whereas this is not the case for RMS, VDV, and PVTV, as shown in column 8-10 in Fig. 5(b), (d), (e), the values between these columns are not significantly different or not highly significantly different. In this case, adding a normal distribution density curve to characterize the vibration perception window has gained its effect, as shown in Fig. 12(h). Note that the normal distribution density curve can add to the difference between column 16 and column 11 or 12 in Fig. 11 is also due to its function to avoid excessive vibration information in scenarios represented by Fig. 12(f) and (g). In Fig. 11, the values in column 7 is not significantly different to the values in column 6, which shows the 25 th parameter set COVV in Table 6 gives consistent evaluation results with the results in Fig. 8, the |Δ − Δ | is only 0.052, whereas this is not the case for MTVV and VDV, as shown in column 5, 6 in Fig. 5(c), (d), the values between the two columns are highly significantly different, the |Δ − Δ | reaches a high level of 0.648 and 0.959, respectively. In this case, involving rotational axis vibration information integration is obviously necessary. The 25th parameter set COVV in Table 6 also suggest the way to integrate rotational axis vibration information as shown in Eq. (21) can perform far better than the way provided by PVTV for all the designed scenarios, PVTV is a simple expansion for RMS and suffers from the difficulty in standardizing the total integration time , however, this is not the case for COVV. None of RMS, MTVV, VDV, PVTV is able to give highly consistent evaluation results with subjective scoring results in group 3, as shown in Fig. 5(f-j). RMS, VDV, and PVTV all used all the vibration information in the entire timeline, but they all fail to express that the real ride comfort in scenario represented by column 3 in Fig. 5(f) is significantly higher then the real ride comfort in scenario represented by column 2 in Fig. 5(f). Fig. 5(i) shows the integration order 4 for VDV is too high, it overly emphasized the peak values in column 2 or highly significantly underestimated the durations in column 3, the peak values comparisons are shown in Fig. 5(h), the durations comparisons are shown in Fig. 3(k-m). Fig. 5(g) shows the integration order 2 for RMS is still too high, because the values in column 2 are significantly higher than the values in column 3, this is on condition that the total integration time in scenarios represented by Fig. 3(k-m) are designed to be the same. Fig. 5(j) shows the values in column 2 is not significantly different to the values in column 3 and is contradictory to the results in Fig. 5(f) where the values in column 2 is significantly lower to the values in column 3, which also suggest the integration order for PVTV should be lowered. Fig. 5(h) shows MTVV neglected the vibration information beyond the maximum part of a local integration time , which highly significantly underestimated the durations in scenario represented by column 3. The values between column 5 and 10 in Fig. 11 are highly significantly different ( = -25.161, < 2.2×10 −16 ), the values between column 10 and 14 in Fig. 11 are significantly different ( = -2.712, = 0.010), which show the 25th parameter set COVV in Table 6 is able to give highly consistent evaluation results with subjective scoring results in group 3. For the 25th parameter set COVV in Table 6, the is 0.82 and is lower than 2 or 4, which is able to emphasize more of the durations rather than the peak values of the sudden halt excitation type and thus gives more consistent evaluation results accordingly; the vibration perception windows covers the 3 sudden halt excitations and assigns comparable weight to the sudden halts beyond the selected reference sudden halt, as shown in Fig. 12(m), which shows the rectified interception item ′ has gained its effect. COVV is contradictory to MTVV in that it also integrate vibration information beyond a local integration time , which helps to integrate necessary vibration information and achieves better consistency.

Conclusions
The devised COVV model achieves greater consistency and accuracy than RMS, MTVV, VDV, and PVTV in evaluation of road ride comfort. Its interactive mechanism is able to express the facts: (1) highly concentration of vibration information alone can add to total discomfort; (2) the interactive strength between localized major vibrations does not change sharply in a range neither too close nor too far away, which means the total discomfort remain stable in this range; (3) not all vibrations in the entire timeline need to be integrated to express real road ride comfort, only the segment where a major vibration together with its interactions with other major vibrations reaches the maximum should be used; (4) the interactive form vary according to the characteristics of a major vibration, a major vibration can interact with a wider range of other major vibrations when its intensity is higher, a major vibration can have more interactive strength with other major vibrations if it contains more vibration information in the or axis; (5) current integration order 2 for RMS, MTVV, PVTV and 4 for VDV are too high; (6) integration order should vary among different translational axes, rotational axis integration order should be suppressed when rotational axis vibration information is high. The single road ride comfort evaluation index COVV can be apply to various actual scenarios where vibration information distribution characteristics vary among axes or spacetime.

Supplementary materials
All the materials, data, and protocols associated with this research are available online at http://dx.doi.org/10.17632/w5zbjpdrsw.2.