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

2.1. 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⁻²), rotational angular velocity (rad s⁻¹), and rotational angles (rad) in the $x$, $y$, and $z$ 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 $ZYX$ coordinate system. The actual directly used vibration information was calculated using Eq. (1):

where $x\left(t\right)$, $y\left(t\right)$, $z\left(t\right)$, $rx\left(t\right)$, $ry\left(t\right)$, and $rz\left(t\right)$ are the actual used accelerations in the $x$, $y$, and $z$ axis and rotational angular velocities in the $x$, $y$, and $z$ axis in moment $t$, respectively. The above variables with a subscript 0 directly represent the corresponding vibration information output directly by the posture sensors. $\mathrm{\Delta }x$, $\mathrm{\Delta }y,$ and $\mathrm{\Delta }z$ are the rotational angles in $x$, $y$, and $z$ axis of the Tait-Bryant $ZYX$ coordinate system, respectively; $g$ is the local gravitational acceleration. The installation ensures $\mathrm{\Delta }z\equiv$ 0, whereas $\mathrm{\Delta }x$, $\mathrm{\Delta }y$, and g were chosen as the recorded maximum density value of the angle in the $x$ and $y$ axes, output directly by the posture sensors, and acceleration in the $z$ axis after the treatment in the right of Eq. (1), respectively.

Fig. 1Subjective scoring scale

Fig. 2Test vehicle and its layout: a) Baojun 730-type test vehicle; b) layout in the test vehicle

a)

b)

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 $z$, $x$, or $y$ 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 $x\left(t\right)$, $y\left(t\right)$, $z\left(t\right)$, $rx\left(t\right)$, $ry\left(t\right)$, and $rz\left(t\right)$ 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.

Fig. 3Vibration information in designed scenarios: a)-c) Crossing 1, 3, and 9 speed bumps at the speed of 40 km h−1, respectively; d), e) crossing 1 and 2 fiber-reinforced polymer (FRP) slab pairs at 40 km h−1, respectively; f)-h) crossing 5 FRP slab pairs with 40, 20, and 4 m spacing at 40 km h−1, respectively; i)-j) directly and indirectly crossing 1 manhole at 50 km h−1, respectively; k)-m) enduring 1, 2, and 3 sudden halts at the speed of 40 km h−1, respectively; n) enduring 1 heavy sway at 40 km h−1; o) enduring 3 consecutive heavy sways at 40 km h−1; p) crossing 1 FRP slab pair and then immediately enduring 1 sudden halt from 40 km h−1; q) crossing 2 FRP slab pairs with 4 m spacing and then immediately enduring 1 heavy sway at 40 km h−1; r), s) enduring 1 heavy sway and simultaneously crossing 1 and 2 FRP slab pairs at 40 km h−1, respectively. The unit for translational axes is m s-2, the unit for rotational axes is rad s-1

a)

b)

c)

d)

e)

f)

g)

h)

i)

j)

k)

l)

m)

n)

o)

p)

q)

r)

s)

Each scenario was experienced by 27-35 people (average 32.895, 61.280 % men). The age range was 20-50 years (average 30.910 years). The height range was 153-183 cm (average 169.714 cm). The weight range was 42-105 kg (average 65.920 kg). The testing sequences are 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).

Fig. 4Excitation configurations schematic diagram for scenarios containing at least one heavy sway

2.2. 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-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 $x$ and $y$ axis, $z$ 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 $k_r$ 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 $k_r$ 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 $a_{wj}\left(t\right)$ is the momentum in axis $j$ at moment $t$, which is frequency weighted acceleration in the translational axes or frequency weighted angular velocity in the rotational axes; $T$ is the total integration time; $\tau$ is the local integration time, its default value is 1 s; $t_0$ is any end moment; $k_j$ is the axial weight, where $k_x=k_y=k_z=$ 1, $k_{rx}=$ 0.63, $k_{ry}=$ 0.4, $k_{rz}=$ 0.2; $k_{aj}$ is the absolute balance weight in axis j, which is used to balance the weights between translational axes and rotational axes; we set $k_{ax}=k_{ay}=k_{az}=$ 1 for standardizing purposes. $k_{arx}=k_{ary}=k_{arz}=k_r$, $k_r$ is the absolute balance weight in the rotational axis, $k_r>$ 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. 5(a-e), (f-j), (k-o), (p-t), respectively. The values at the top of Fig. 5 are the $t$-values and $p$-values of a $t$-test between the below column and the adjacent right column. The columns for subjective scoring results are arranged so that t-values are always lower or equal to 0; this column sequence was also applied to the corresponding objective evaluation results. Fig. 5 only shows the major vibration axis information; the exceptions are group 4 scenarios due to the multiple major vibrational axes. We ensured the other translational axes only contained minor vibration information and are not demonstrated in Fig. 5. The objective evaluation results in the rotational axes are also not provided independently in Fig. 5 because both ISO 2631-1 and BS 6841 suggest that using vibration information in the translational axes is sufficient.

Fig. 5Boxplot 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

a)

b)

c)

d)

e)

f)

g)

h)

i)

j)

k)

l)

m)

n)

o)

p)

q)

r)

s)

t)

u)

v)

w)

x)

y)

z)

The t-test results are described as: (1) when $p\ge$ 0.05, there is no significant difference; (2) when 0.001 $\le p<$ 0.05, there is significant difference; (3) when $p<$ 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 $t$- values and $p$-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 $t$- value and $p$-value in the additional $t$-test between columns 1 and 6 in Fig. 5(e) are 4.947 and 5.64×10⁻⁶, 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 ($t=$ −5.058, $p=$ 6.356 × 10⁻⁶); (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.

2.3. 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 $E_1\left(t\right)$ using:

2. Find a consecutive interval that covers the largest $E_1\left(t\right)$ point in the residual timeline, can maximally integrate total $E_1\left(t\right)$, and is limited to a maximum length of local integration time τ. This procedure is shown in Eqs. (7) and (8):

where $E_{1,i}$ is the total translational axes energy for the $i$th local impact; $t_i$ is the end moment for the $i$th local impact, $t_0<$ 0; and $t_{i,max}$ is the moment when $E_1\left(t\right)$ in the residual timeline reaches its maximal value.

3. Remove the vibration information included in $E_{1,i}$ from the existing timeline.

4. If $\underset{t}{max}{\int }_{t-\tau }^t{E}_1\left(t\right)dt<d\times \underset{i}{max}\left(E_{1,i}\right)$, on the condition that $t\in \left[0,T\right]-\sum _{j=0}^{i-1}[t_j-\tau ,t_j]$ and decreasing factor $d\in \left(\mathrm{0,1}\right)$, then proceed to step 5. Otherwise, set $i=i+1$ and return to step 2.

5. If any $E_{1,i}$ is lower than $d\times \underset{i}{max}\left(E_{1,i}\right)$, delete this $E_{1,i}$.

6. Rename the left $E_{1,i}$ from its emerging sequence to ensure the $i$ sequence starts from 1. With each time step, add 1, then the slice corresponding to $E_{1,i}$ is defined as the $i$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 $t$- and $p$-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 $L_i$ is the window length for the $i$th local impact; $k_1$ and $k_2$ are the scale parameters, $k_1>$ 0 and $k_2>$ 0; $E_{2,i}$ is the total influence energy for the $i$th local impact, its integration interval is the same as $E_{1,i}$; $E_2\left(t\right)$ is the influence energy density at moment $t$; $k_{xi}$, $k_{yi}$, and $k_{zi}$ are the axial influence energy weight in the $x$, $y$, and $z$ axis, respectively. $k_{xi}>k_{zi}$, $k_{xi}>$ 0, $k_{yi}>$ 0, $k_{zi}=$ 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 $h\left(t\right)$ is the bathtub curve function; $\lambda$ is the failure rate parameter, $\lambda >$ 0; $\eta$ is the interval parameter, $\eta =$ 1 for standardizing purposes; and $\alpha$ and $\beta$ are the shape parameters. We set $\alpha =$ 0.2 and $\beta =$ 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:

where $w_i\left(t\right)$ is the vibration perception window for local impact $i$; $c$ is the intercept item, 0 $<c<$ 1 + 0.559$\lambda$. This range of $c$ can ensure the window height does not exceed 1; a higher $c$ indicates more vibration information beyond the reference local impact is integrated; $\sigma$ is the variance for the normal density curve; $t\text{'}$ is the moment when the bathtub curve as is expressed in Eq. (11) reaches its lowest point; and $R$ is the rectify coefficient, which is introduced for standardizing purposes and, combined with the range of $c$, 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 $V_{i,m,j}$ is the vibration information integrated by local impact $i$ from local impact $m$ on the $j$ axis, $m\ne i$; $n$ is the basic integration order, $n>$ 0.

Fig. 6Boxplot of VDV in the rotational axis for group 1: a) x-axis, b) y-axis, and c) z-axis

a)

b)

c)

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 $x$-axis than the other columns in group 1. Columns 6 and 7 in group 1 also contain much more vibration information in the rotational $y$-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 $z$-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 $k_z$ to 1 for standardizing purposes.

2.4. 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 $x$, $y$, $z$, rotational $x$, rotational $y$, and rotational $z$ axis are simultaneously put as the input variables. For COVV, the input variables are directly themselves.

The L₂₇(3¹³) 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 $k_x$ and $k_y$ are added to the parameters is because different and even the same standard has varied specification for these parameters. ISO 2631-1 sets $k_x=$1.4 and $k_y=$1.4 for the synthetic of translational RMS while $k_x=$1 and $k_y=$1 for PVTV; BS 6841 sets $k_x=$1 and $k_y=$1 for the synthetic of translational RMS; Mansfield and Maeda [31] suggest $k_x=$2.2 and $k_y=$2.4. The NMSEs for RMS, MTVV, VDV, PVTV and COVV in the parameter set in 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 $\tau =$1.2, $d=$0.4, $k_1=$0.5, $k_2=$1.8, $n=$1.5, $\lambda =$0.6, $\sigma =$0.02, $p=$0.6, $k_x=$0.8, $k_y=$1.4, $k_r=$14, $k_{xi}=$1.8, $k_{yi}=$1, the boxplot for the best parameter set COVV is shown in Fig. 7. Fig. 7 has the same serial number as Fig. 8, the $t$ 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. ${\mathrm{\Delta }}_1$ 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. $\left|{\mathrm{\Delta }}_1-{\mathrm{\Delta }}_2\right|$ 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 ($t=$–0.585, $p=$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 ($t=$–3.8498, $p=$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.

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 $x$ axis vibration information in the reference local impact, which means to increase the intercept item $c$.

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 $x$, $z$ 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 ($t=$–0.019, $p=$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 ($t=$–3.055, $p=$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 $z$ 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 $c$ in Eq. (12) should be replaced by $c\text{'}$ in Eq. (17). In order to further solve the third major problem in Fig. 7, we can lower the $k_y$, 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 $y$ axis is integrated:

where $c\text{'}$ 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 $c$ in Eq. (12) is replaced by $c\text{'}$ 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 $\tau =$ 0.8, $d=$0.2, $k_1=$0.6, $k_2=$1.6, $n=$1.6, $\lambda =$0.4, $\sigma =$0.01, $p=$0.6, $k_x=$0.8, $k_y=$0.6, $k_r=$15, $k_{xi}=$2.1, $k_{yi}=$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.

Fig. 7Boxplot 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. 8Boxplot of subjective scorings for all the designed scenarios: a) the results for scenario 1-10; b) the results for scenario 11-19

a)

b)

$t$ value and $p$ value between column 10 and 14 in Fig. 9 is –0.877, 0.385, respectively, it means increasing $c$ 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 $n$ 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 $z$ 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 $n_j$ is the integration order for axis $j$, $n_x\ne n_z$, $n_y\ne n_z$, $n_{rx}=n_{ry}=n_{rz}=n$.

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 $k_y$ 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 $n_{sj}$ is introduced, and Eq. (19) should be replaced by Eq. (21):

where $n_{sj}$ is the suppression order for axis $j$, $n_{sx}=n_{sy}=n_{sz}=$1, which means there is no suppression for the translational axes, $n_{srx}=n_{sry}=n_{srz}=n_{sr}$, $n_{sr}$ is the suppression order in the rotational axes, 0 $<n_{sr}<$1.

Fig. 9Boxplot 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-19

a)

b)

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₆₄(4²¹) 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(l). The end moments $t_1$, $t_2$ are 23.395 s, 63.095 s, respectively. The moments when the bathtub curve as is expressed in Eq. (11) reaches its lowest point $t\text{'}$ for the 1, 2th local impact are both 0.509. The translational axes energy $E_{\mathrm{1,1}}$, $E_{\mathrm{1,2}}$ are 13.616 m² s^-3, 12.640 m² s^-3, respectively. The influence energy $E_{\mathrm{2,1}}$, $E_{\mathrm{2,2}}$ are 30.339 m² s^-3, 28.202 m² s^-3, respectively. The window length $L_1$, $L_2$ are 215.707 s, 189.824 s, respectively. The rectify coefficient $R$ for the 1, 2th local impact are both 0.0018. As such, the $w_i\left(t\right)$ is fixed. For Eq. (18) in step 3, the ${\int }_{t_i-\tau }^{t_i}{a}_{wx}^{n_x}\left(t\right)dt\text{,}$${\int }_{t_i-\tau }^{t_i}{a}_{wy}^{n_y}\left(t\right)dt\text{,}$${\int }_{t_i-\tau }^{t_i}{a}_{wz}^{n_z}\left(t\right)dt\text{,}$${\int }_{t_i-\tau }^{t_i}{a}_{wrx}^{n_{rx}}\left(t\right)dt\text{,}$${\int }_{t_i-\tau }^{t_i}{a}_{wry}^{n_{ry}}\left(t\right)dt$, ${\int }_{t_i-\tau }^{t_i}{a}_{wrz}^{n_z}\left(t\right)dt$ can be calculated by using the non-weighted accelerations in the $x$, $y$, and $z$ axis and non-weighted rotational angular velocities in the $x$, $y$, and $z$ axis in moment $t$ as are calculated by Eq. (1) and using a Hilbert-Huang Transform (HHT) to calculate the frequency weighted ones [44]. When $i=$1, they are 2.195, 0.114, 0.294, 0.005, 0.058, 0.004, respectively. When $i=$2 they are 2.078, 0.130, 0.274, 0.008, 0.054, 0.001, respectively. With this calculation, the rectified intercept item $c\text{'}$ for the 1, 2th local impact can be easily defined according to Eq. (20) and they are 1.051, 1.053, respectively. In step 3 in Fig. 10, the vibration information integrated by local impact $i$ from local impact $m$ on the $j$ axis $V_{i,m,j}$ can be defined with the already calculated $w_i\left(t\right)$ and ${\int }_{t_i-\tau }^{t_i}{a}_{wj}^{n_j}\left(t\right)dt$. When $i=$1, $m=$2, $j=x$, $V_{i,m,j}=$1.943; when $i=$2, $m=$1, $j=x$, $V_{i,m,j}=$1.854; the results for other combinations of $i$, $m$, $j$ are not shown here for simplification purposes. Finally, according to Eq. (21), the maximum vibration information is integrated when $i=$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₆₄(4²¹) orthogonal test in Table 5.

Fig. 10Final flowchart for the calculation of COVV