Identiﬁcation of Road Proﬁle Parameters from Vehicle Suspension Dynamics for Control of Damping

: Concept of symmetry covers physical link between road proﬁle form, vehicle dynamic characteristics, and speed conjunction. Symmetry frame between these items is asymmetric itself and has no direct expression, but it affects a vibration level on the vehicle and driving comfort. Usually, we can change only the vehicle’s speed to achieve desired vibrations level of the driver and passengers. Recently, vehicle dynamic characteristics can be changed depending on its damping system structure, but these solutions are limited by construction and control possibilities and evidently represented by symmetric dependency between road input and the resulting acceleration of the vehicle. The main limitation of this process is to have a reliable value of the existing road proﬁle that is mainly deﬁned by road category but unpredictable for each road distance. Functional road proﬁle calculations are provided in this article, where power spectral density (further-PSD) and waviness of the road play the main role in delineating road proﬁle parameters. Furthermore, the transfer function system was created using full car dynamic model analysis. Values on vehicle suspension’s effects on acceleration were obtained from vehicle speed and road roughness. Acceleration values and transfer function were used to calculate PSD value quickly and practically. This calculated result can be formed as a control value to the vehicle damping control. In addition, the provided methodology became useful to determine road quality for adjustment of vehicle suspension parameters and set safe driving characteristics, which became part of driver assistant systems or autonomous driving mode.


Introduction
Ride quality, vehicle handling, fuel consumption, and tire wear depend on road quality [1]. Road impact to vehicle and vehicle impact to the road degradation creates a symmetric entity with the reference frame represented by dynamic features of the vehicle. Bad road condition and potholes cause damage to vehicles, and regular monitoring of road quality is required to avoid the high costs of damage the vehicles [2]. Detailed information about road profile is key for transport authorities maintaining adequate ride quality for vehicles, vehicle designers and engineers developing suspension systems. Bad road profile requires a smart suspension system from passive to semi-active [3][4][5][6][7] and active control [8][9][10][11] in automotive technology. So far, a universal solution for vibrations control has not been achieved, but the possibility of controlling vehicle dynamic response allows for the adjustment of drivers' and passengers' felt vibrations [12]. Due to road irregularities' unpredictability, controlling the dynamic properties of the vehicle suspension becomes a complex task consisting of: (i) road irregularities' value reading process, (ii) generation of the control signals, and (iii) control of the regulation process [13]. Then, although vibrations control has its solutions, reading of the irregularities becomes a real challenge. Despite that, there are ways to read road profile. Different methods and equipment to measure the road surface characteristics are briefly discussed in [14][15][16]. They can be divided into three categories: system response-based estimation, direct measurement, and non-contact measurement [17,18]. A dozen various inexpensive response-based methods of road profile recognition are reviewed and presented in [19]. Some of them use accelerometers mounted in different places and transfer functions for profile reconstruction or roughness estimation [19][20][21][22]. There are papers where an accelerometer is used for control of smart suspension [6,23,24]. Research where road surface waviness used to control damping has not been published yet.
At the same time, we have found various roughness or ride quality indices used in literature, i.e., international roughness index (IRI), ride quality index (RQI), whole-body vibration (a wz ), ride number (RN) and others [14][15][16]25]. Some of them are associated with standards [15]. Standard ISO 8608:2016, which defines road surface profiles [26], provides two essential parameters: the degree of roughness G d (n 0 ) or G d (Ω 0 ) identified as power spectral density (further-PSD) of the vertical road profile displacement, and the waviness w. The main usefulness of road waviness value for the transportation community consists of [27]: (i) it modulates power distribution of the road elevation spectrum; (ii) short wavelengths (with w < 2) affect road surface, while long wavelengths (w > 2) affect the whole structure. On the other hand, they negatively influence vehicle handling and ride comfort. Waviness has a much higher impact on the total power than a degree of roughness G d alone, but vehicle velocity v and waviness w modulate PSD of the vertical vibrations as v w−1 . Also, the root mean square (further-RMS) value of particular vibration responses as well as IRI are proportional to v (w−1)/2 [28,29].
Use of acceleration value from wheel hub in our case is based on the possibility of collecting acceleration data practically by mounting the general purpose accelerometer (bandwidth not less than 1 Hz to 30 Hz, measurement range up to 50 m/s 2 ) on the lower end of the shock absorber. While the accelerometer is an inertial sensor, it requires no extra reference point to install it. Choice of waviness (w) as an essential road profile parameter better reflects real road surface and is possible to define with fewer efforts practically.
In contrast to the majority of existing publications, we are using "two-split" and "three-line" models to obtain three "waviness" values in the desired frequency range. Our research proves that this approach can be used to control damping parameters of car suspension. Moreover, it allows for realizing precise and smooth control of the damper damping coefficient, while [21] uses values of G d and w. Their approach reproduces the generated road profiles according to the assigned types very well, but using such a method to control shock absorbers will be insufficient since the damping coefficient cannot vary smoothly. The damping coefficient will be adjusted only when the profile value will exceed a certain range. Another group of scientists [30] used linear fitting in ten octaves. They showed good road classification for two-wheeled vehicles. However, this is not appropriate for controlling the damping coefficient due to the huge matrix of possible combinations.
As stated in [19], many studies use a quarter-vehicle model and passive suspension system, while transfer function and other methods mostly use the passive system. We plan to control our developed (semi-)active suspension using information obtained from a segment of passed road profile.
This paper provides a methodology to detect road profile characteristics from acceleration values, which can easily be detected for in-situ mounted sensors. It is possible to estimate the power spectral density of road profile from wheel hub acceleration using transfer functions. Obtained PSD allows for the adjustment of vehicle driving mode in the case of driver assistant systems or fully automated drive mode, set damping parameters for adjustable shock absorbers or provide driving hints for the driver. This data processing with derived transfer function is fast and requires lower processing power than methods presented in [6,30]. Therefore, it can be processed on board or even by an internal shock absorber processor. The methodology of data extraction is deeply related to the concept of symmetry due to our approach of implementing direct data transformation using polynomial; therefore, acceleration values can be processed backwards from vehicle Symmetry 2021, 13, 1149 3 of 11 body vibration trough a symmetric transformation procedure. Vehicle body and wheel hub vibrations are related with defined transfer functions. The provided methodology based on extraction of waviness values uses parameters from a full-car dynamic model, which already compared with experimental data [31].

Materials and Methods
The proposed methodology of road profile waviness is classified by obtained accelerations of a hub of each wheel, and requires some mathematical manipulation, which is provided below in the text. The road profile here is generated using the provided methodology, which ensures concordance of obtained profile to the real road surface of the corresponding category. A generation of longitudinal road profiles based on sums of various length sinusoidal waves with different amplitudes and phases has been applied by many authors [32][33][34][35][36][37][38]. The one-sided PSD spectrum was split into N equal parts according to frequency, and for each is generated a harmonic sample with central frequency, spectral amplitude. Random values of phase angles ϕ i were uniformly distributed in the range from 0 to 2π. To create artificial road profile samples were superimposed in the spatial domain according to Equations (1) and (2): where G d (Ω i ) is displacement PSD at angular spatial frequency Ω i , Z i is the amplitude of ith harmonic, Ω U , Ω L are upper and lower angular spatial frequencies in the PSD spectrum [rad/m], ∆Ω is the width of each frequency band, and Ω i = Ω L + (i − 1)∆Ω is the angular spatial frequency of ith harmonic or ith frequency band. In our case, the number of harmonic samples was N = 1000. For general on-road measurements, recommended lower and upper limits of angular spatial frequencies (=2πn) are 2π × 0.01 rad/m and 2π × 10 rad/m, respectively [26]. The required step size of road distance x, according to the data sampling theorem [39], is equal to ∆x ≤ 2π 2Ω U = 0.05 m. A smoother profile could be achieved by minimising steps. However, it only increases the required computational efforts. To get more various road profiles, instead of the linear fitting recommended by ISO standard [26] in Equation (2), we used two split fittings provided in [40]: The degree of roughness G d (Ω 0 ) is PSD value at the reference angular spatial frequency. The extreme spatial frequencies Ω 1 = 2π × 0.21 rad/m and Ω 2 = 2π × 1.22 rad/m produced a minimal error describing the Swedish road network [40]. A new random set of initial phases was generated at each 2π/Ω L distance to avoid a replication of the road profile. Additionally, we used the Welch type windowing function with the value of exponent equal to 10 to minimise the appearance of abrupt shifts at the connections between profile segments [41]. The road profiles were calculated through written MATLAB™ code. The profile for the right vehicle side (both front and rear wheel) was formed by reversing the order of the left side profile related to the vehicle. For comparison, identical sequences of random numbers were used to generate each profile.
Our developed dynamic model of the whole car [31] was used to simulate vehicle passing these profiles. The model included ten masses (car's mass, four wheels' masses, Symmetry 2021, 13, 1149 4 of 11 masses of four passengers and baggage box) and moments of inertia about the X and Y axes. Twelve generalised coordinates defined the behaviour of the modelled car: vertical axis coordinate-Z; rotation angles about X and Y axes-ϕ x and ϕ y ; z 1 , z 2 , z 3 , z 4 are the vertical coordinates of wheels; z 13 , z 14 , z 15 , z 16 are coordinates of the centres of mass of the driver and passengers, and z 18 is coordinate of the centre of baggage's mass. The model was build using Lagrange's equation of the second type. These equations were solved analytically, and after differentiating, we obtained the final system, consisting of the equations, presented in this form: where a, b, c, d, d* (with corresponding indexes) are coefficients of equations, derived from matrixes of stiffness, dissipation and inertia; q i is generalised coordinate, in respect to which differentiation was applied during the building of a system of equations; n is the number of generalised coordinates; and η 1 , η 2 , η 3 , η 4 are coordinates along which the system of the car is excited kinematically. Each overdot is a shorthand for a time derivative. The rear wheels are excited with the same road profile as the front wheels, with corresponding shifts by a distance equal to the distance between the front and rear wheel centres. The research was performed by simulating the car, Range Rover Evoque, with the driver and front passenger but without rear passengers or baggage (all parameters presented in supplementary materials Table S1). The second derivative of the ..   z 4 coordinates was taken as road generated acceleration values from road excitation for further analysis. Road PSD was estimated on the relationship between the road surface and vehicle acceleration via a transfer function: where H(.) and PSD road (.) or PSD acc (.) are the transfer function and PSD for a frequency Ω or n due to road profile or wheel accelerations, respectively. Waviness values were extracted from PSD spectra using the straight line for logarithmic values in the ranges specified for each waviness. To evaluate road surface PSD estimation quality, we used the normalised root mean squared error parameter (NRMSE). This is given by the RMS error between the estimated valueŷ and true value y, normalised to the RMS of road roughness true value, and expressed as a percentage [42]:

Results and Discussion
Initially, the samples of road profile with various w 1 , w 2 , w 3 waviness values and with the ISO road class B parameter G d (Ω 0 ) = 4 × 10 −6 m 3 , Ω 0 = 1 rad/m [26] have been generated. Simulated datasets were named using the naming rule as "prw 1 w 2 w 3 ", where w i can have the values of 1, 2 and 3. The length of each generated profile was 1000 m. Displacement PSD (G d ) spectra of profiles passing at speed v were obtained using fast Fourier transformation (FFT) with a 50% overlapping Hamming window of length N = 2 × f loor(1000[m]/v[m/s]/∆t[s]/2) points (with buffer size 1000 m) and time step ∆t = 0.01 s (sampling frequency 100 Hz). Implemented here, floor is the function rounding toward negative infinity. As is recommended by ISO standard [26], the smoothed PSD spectra of pr222 profile at car speed v = 20, 50, 70, 90 or 130 km/h are shown in Figure 1. From PSD spectra were estimated three cases of waviness values: (1) Using Equation (4)  (2) Using Equation (4) but after smoothing.
(3) After smoothing along a straight line using Equation (3), starting from octave with lower cut-off frequency 0.0557 cycles/m and centre frequency 0.0625 cycles/m up to octave with max value.
= 0.01 s (sampling frequency 100 Hz). Implemented here, floor is the function rounding toward negative infinity. As is recommended by ISO standard [26], the smoothed PSD spectra of pr222 profile at car speed v = 20, 50, 70, 90 or 130 km/h are shown in Figure 1. From PSD spectra were estimated three cases of waviness values: (1) Using Equation (4)   Waviness values estimated from pr222 profile at various speed using Hamming, Hanning, Bartlett or rectangular window in FFT are presented in Table 1. Increasing speed from 20 km/h to 130 km/h, the maximum frequency of one-sided PSD spectra decreases from 8.99 cycles/m down to 1.38 cycles/m. At the same time, powers of higher frequencies (see Figure 1) increase and this influences the decrease of waviness values w2, w3 (case 1 and 2) and w (case 3), increasing speed. Values of w3 without smoothing could be negative due to approximation of few points in the much-dispersed area (see inset of Figure 1, and columns of Table 1 with the header "1000 m") or not extracted at all such as when buffer size is 100 m. Values of w1 without smoothing are a bit higher than after smoothing. This could be explained by the fact that lower and higher frequencies are accounted for during the smoothing procedure. The smallest waviness deviations from theoretical values of ISO recommended ones are obtained using a rectangular window and smoothing procedure.  Waviness values estimated from pr222 profile at various speed using Hamming, Hanning, Bartlett or rectangular window in FFT are presented in Table 1. Increasing speed from 20 km/h to 130 km/h, the maximum frequency of one-sided PSD spectra decreases from 8.99 cycles/m down to 1.38 cycles/m. At the same time, powers of higher frequencies (see Figure 1) increase and this influences the decrease of waviness values w 2 , w 3 (case 1 and 2) and w (case 3), increasing speed. Values of w 3 without smoothing could be negative due to approximation of few points in the much-dispersed area (see inset of Figure 1, and columns of Table 1 with the header "1000 m") or not extracted at all such as when buffer size is 100 m. Values of w 1 without smoothing are a bit higher than after smoothing. This could be explained by the fact that lower and higher frequencies are accounted for during the smoothing procedure. The smallest waviness deviations from theoretical values of ISO recommended ones are obtained using a rectangular window and smoothing procedure.
Similar tendencies are observed, and mean waviness values in the distance of 1000 m are estimated with a buffer size of 100 m (see columns of Table 1 with the header "100 m"). However, a more detailed analysis of other profiles showed that the rectangular window gives the worst results only in some cases (see Tables S2-S4) but has the smallest error while others vary from case to case, which means the maximum difference between the values used to generate the profile and those extracted. It should be noted that a decrease of FFT buffer size increases the variation of waviness, but mean values remain close to the ones used in generation. Our previous investigation showed that the straight-line approximation of PSD of the road profile proposed by ISO standard doesn't contribute suitable information required to control damping parameters of car suspension [43]. In the light of previous observations, the third waviness index is not suitable for the control coefficient of damping of the suspension. Moreover, from the optimal damping coefficient calculations for various profiles in our previous research [31], we observed that the changing of the w 3 value influences only tens of the damping coefficient. Thus, we are using only values of w 1 and w 2 waviness and rectangular window in FFT calculation. Our aim is to get waviness value, not PSD; therefore, we are not using smoothing. Additionally, to minimise the calculation resources, it is possible to use, for example, MatLab function pwelch with a shorter vector of chosen frequencies.
Our next step is to obtain information about the passing road profile segment, and to evaluate PSD and extract waviness w 1 and w 2 . One should choose which is optimal from all road profile generation methods, requiring less computational resources and shorter passed distance, being more precise. This procedure demands a separate detailed study that is not our aim. We estimated PSD of road profile from wheel acceleration values using transfer functions obtained from the full-car model. Using Cramer's method and applying the principle of superposition to the system in the form of Equation (5), we can obtain the theoretical transfer functions for four accelerations for each wheel hub .. z i for each displacement input η j : where Laplace's complex variable s = √ −1Ωv = √ −12πnv. A more detailed form of H 11 is presented in the supplementary material Equation (S1). Calculated acceleration transfer functions H ii and cross-transfer functions of first point H 1j at 70 km/h are presented in Figure 2a. As can be expected, transfer functions are the same for front left and front right wheels or for the rear left and rear right wheels. Meantime cross-transfer function impacts are more than two orders of magnitude smaller. Consequently, it is enough to use only the H ii transfer function for further calculations. Comparison of the simulation transfer function H s 11 , evaluated using Equation (6) and PSD of input (displacement) and output (acceleration) of one wheel in our full-car model simulation, and theoretical transfer function varying car speed are shown in Figure 2b. One can observe good agreement between them in the frequency range interesting us. NRMSE of simulated and theoretical transfer functions at all velocities is less than 65% for all range values and decreases to 8% for a range from 0.21 to 1.22 cycles/m. Thus, using the transfer function at chosen frequency range in restoring road profile from wheel hub accelerations or without smoothing, one could expect precise results.
where Laplace's complex variable 1 1 2 s v nv       . A more detailed form of H11 is presented in the supplementary material Equation (S1). Calculated acceleration transfer functions Hii and cross-transfer functions of first point H1j at 70 km/h are presented in Figure 2a. As can be expected, transfer functions are the same for front left and front right wheels or for the rear left and rear right wheels. Meantime cross-transfer function impacts are more than two orders of magnitude smaller. Consequently, it is enough to use only the Hii transfer function for further calculations. Comparison of the simulation transfer function s 11 H , evaluated using Equation (6) and PSD of input (displacement) and output (acceleration) of one wheel in our full-car model simulation, and theoretical transfer function varying car speed are shown in Figure 2b. One can observe good agreement between them in the frequency range interesting us. NRMSE of simulated and theoretical transfer functions at all velocities is less than 65% for all range values and decreases to 8% for a range from 0.21 to 1.22 cycles/m. Thus, using the transfer function at chosen frequency range in restoring road profile from wheel hub accelerations or without smoothing, one could expect precise results. The transfer function method has some disadvantages. Usually, it is used with constant speed (a small-time span or buffer length can eliminate). In addition, few potholes on the road surface can significantly change estimated values for the segment. The provided methods, described in [44][45][46], can minimise road pothole impact to real road PSD values.
Any other method could be used for road PSD profile restoring from wheel hub acceleration values. For example, another possible way to identify the road profile, knowing the dynamic responses of the car (the chassis acceleration and deflection of the suspension), is a new fast and simple technique, called independent component analysis (ICA), based on the so-called inverse problem [47]. The advantages of the ICA are that the identification of road excitation is simpler and faster compared to the other estimation techniques. Also, it is robust to the sprung mass variation. Furthermore, the addition of noise The transfer function method has some disadvantages. Usually, it is used with constant speed (a small-time span or buffer length can eliminate). In addition, few potholes on the road surface can significantly change estimated values for the segment. The provided methods, described in [44][45][46], can minimise road pothole impact to real road PSD values.
Any other method could be used for road PSD profile restoring from wheel hub acceleration values. For example, another possible way to identify the road profile, knowing the dynamic responses of the car (the chassis acceleration and deflection of the suspension), is a new fast and simple technique, called independent component analysis (ICA), based on the so-called inverse problem [47]. The advantages of the ICA are that the identification of road excitation is simpler and faster compared to the other estimation techniques. Also, it is robust to the sprung mass variation. Furthermore, the addition of noise has a negligible effect on the estimation process. Therefore, there is no need to use a more complex full-car model to estimate the road profile [47].
Finally, we verified how waviness values extracted from PSD spectra of recovered road profile do meet values obtained from PSD spectra of generated road profile. Recovering road roughness data from generated PSD spectra expresses the concept of symmetry in the research idea and methodology. Minimal and maximal values of both w 1 , w 2 waviness and maximal NRMSE value of each profile at various speeds are presented in Table 2 (estimated using rectangular window without smoothing in FFT, buffer 100 m).

Conclusions
Provided analysis of vehicle suspension and wheel hub accelerations revealed an apparent possibility of extracting road characteristics from sensor data reliably. Sensing road parameters just with a simple technique opens new horizons in infrastructure observation using background data streamed from passing cars, detecting defect information, damage depth, and development in the time domain. This information will be useful for car producers as a method of vehicle operation monitoring in terms of safety and lifetime for the entire car body.
As the main achievement of this research, we should claim successful deciphering of accelerometer data into road category according to ISO classification and set damping characteristics of the vehicle suspension system. The created methodology with transfer function implementation in fast data processing provides a green light for in-situ mounted low computing power systems. Such installation inside shock-absorbers will transfer road characteristics to the central vehicle safety block and set shock absorber damping values to minimise vibration values of the vehicle body. These values, provided in Table 2, became extremely important in the case of driver assistant systems and autonomous driving systems.
Results of this research are still begging for continuation, because randomly distributed hills and holes on road surface requires a more detailed analysis of sensor signals; recent papers cover only standard defined road surface analysis.
Future research will evaluate non-standard road artefact evaluation and data fusion with a navigation system and issue the commands for vehicle suspension from obtained sensor data.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/sym13071149/s1, Figure S1: Dynamic model of the car body. Table S1: Range Rover Evoque simulation parameters. Table S2: Averaged waviness values of various profiles at 90 km/h speed and buffer size 100 m. Table S3: Waviness values of various profiles at 90 km/h speed and buffer size 1000 m. Table S4: Waviness values of various profiles using rectangular window and buffer size 1000 m. Equation (S1): Transfer function H 11 . Funding: This paper was supported by AutoDrive-"Advancing fail-aware, fail-safe, and failoperational electronic components, systems, and architectures for fully automated driving to make future mobility safer, affordable, and end-user acceptable" project, which has received funding from the European Union's Horizon 2020 research and innovation programme for 2014-2020, through the Electronic Component Systems for European Leadership (ECSEL) Joint Undertaking under grant agreement No. 737469.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.