**1. Introduction**

The installation of structural health monitoring systems (SHMSs) has emerged as an effective approach with which to gauge bridge responses, detect potential damage, and evaluate the overall condition of bridges [1,2]. To date, numerous SHMSs comprising a multitude of sensors have been set up on a range of modern bridges, notably the Qingma bridge and the Hong Kong–Zhuhai–Macao bridge [1]. Such systems accumulate extensive data during the operational phase of these structures. The data amassed are analysed using various data-driven assessment algorithms, including machine learning, deep learning, and other intelligence approaches, which are proposed based on the dataset recorded via the SHMS [3]. These data-driven algorithms aid in the extraction of crucial features' parameters to evaluate the health condition of bridges. Despite the SHMS's capability in detecting anomalies in loading and responses, offering real-time information for timely safety assessments after disasters and extreme events, its application is largely limited to long-span or keynote bridges. This limitation is due to the high costs associated with its instalment and maintenance. Additionally, anomalies within the data are inevitable due to potential sensor faults, transmission failures, and the aging of sensor components [4–6]. These challenges necessitate manual inspection as the primary method for monitoring the health status of small- to medium-span bridges.

Moving force identification (MFI) plays a crucial role in the health state monitoring of bridges and can provide valuable references for bridges' management and maintenance during their service period [7–9]. Given recent reports of bridge collapse accidents due to overload, it is crucial for bridge managers to have timely knowledge of moving forces as vehicles cross bridges. While an installed bridge weigh-in-motion device can directly measure moving forces, these devices need to be installed during the bridge construction stage, and their functionality can be affected by the pavement status [10,11]. Moreover, the installation and maintenance of such devices are generally expensive, making it impractical to deploy these devices on every bridge, especially considering the large quantity of smalland medium-span bridges.

To this end, there has been growing interest in the MFI technique over the past few decades. This technique identifies moving forces using bridge responses captured by installing sensors such as strain gauges, displacement sensors, and accelerometers [12–21]. Law and Chan proposed a time domain method to identify moving loads by combining the use of bending moments and acceleration responses [12], while Zhu and Law developed a time domain method for a multi-span continuous bridge using strain and displacement measurements [13]. Liu et al. proposed a load identification method utilising a new global kernel function matrix to reduce ill-posedness and improve the identification accuracy under the displacement response [14,15]. The literature [16] also shows that the identification of the moving load can be determined with the bending moment and accelerations. A notable contribution was made by Yu and Chan, who proposed a frequency–time domain method for identifying moving vehicle axle loads, relying on the bending moment responses [17]. Lage et al. also proposed a method using a concept of response transmissibility based on measuring the displacement. Their methods include two steps: initially identifying the number of forces and their respective locations, followed by reconstructing the load vector [18]. Aucejo and De proposed a space–force multiplicative regularisation method to avoid the preliminary definition of any regularisation parameter [19]. He et al. developed a frequency domain method for load bounds' identification in uncertain structures [20]. These studies showcase the major identification techniques categorized into two primary methods: time domain [12–16] and frequency domain methods [17–21].

For the time domain method, the motion equation is initially transformed into a second-order differential equation in modal coordinates, corresponding to the structural response and load. Subsequently, the differential equation is decoupled via the convolution integral within the time domain [12–16]. By contrast, the frequency domain method starts by transforming the time domain motion equation into its frequency domain counterpart using Fourier transformation. The subsequent stage involves solving the equation in the frequency domain. The dynamic spectrum is determined based on the relationship between the transfer function matrix and the system response spectrum [17–21].

However, owing to unknown initial conditions and state variables, the MFI presents a typical ill-posed problem, indicating that the inverse identification problems may not have unique solutions [22]. The discrete MFI equations are ill-conditioned, and the identified results are sensitive to external environmental noise. To overcome the ill-posed problem, effective techniques, including regularisation techniques, have been adopted [23]. For instance, Chen and Chan proposed various methods, such as the truncated generalized singular value decomposition method (TGSVD), piecewise polynomial truncated singular value decomposition (PP-TSVD) method, modified preconditioned conjugate gradient (M-PCG) method, and preconditioned least square QR factorization method, to solve the ill-posed problems in identifying the moving force from the response of the bridge deck [24–27].

In further research, Chen et al. compared these four methods to evaluate their overall performance through numerical simulations and laboratory verification [28]. Recently, Chen proposed a modified truncated singular value decomposition (MTSVD) method for moving force identification, aiming to overcome the ill-posed problems, and a comparative study was conducted with its conventional counterparts: the SVD and TSVD methods [29]. Pan and his co-authors proposed MFI methods to address the ill-posed problem [30–33]. In 2018, Pan et al. proposed a hybrid moving force identification method based on weighted L1-norm regularisation and redundant concatenated dictionary [30]. The dictionary, including trigonometric functions and rectangular functions, serves to match the principal features of unknown moving forces, while the weighted L1-norm regularisation method forms the MFI equation. Later, Pan et al. introduced a matrix regularisation-based method for solving large-scale inverse problems for force identification. Moreover, they proposed a constrained sparse regularisation-based method, taking into account unknown moving forces and initial conditions [31,32]. Recently, Pan et al. devised an equivalent load-based method for identifying the gross weight of a vehicle moving on a beam-like bridge [33]. Zhou et al. proposed an integral time domain method, effectively eliminating errors generated in the discrete unit impulse response function [34]. Zhang et al. proposed a fresh MFI method based on learning dictionary with double sparsity to a fixed dictionary for expressing moving forces with sufficient sparsity. The sparse K-singular value decomposition (K-SVD) algorithm was employed to realize the learning process [35]. He et al. developed a time domain MFI method that considered the uncertainty in the finite element model and verified the effectiveness and superiority of the method using numerical and experimental examples [36].

The aforementioned literature establishes that the proposed MFI methods can be transformed into a problem of solving linear equations, applicable to both time domain and frequency domain methods [24–27,30–33]. The Newmark-β method, valued for its stability, has been widely applied to forward structural dynamic problems [37]. In recent years, this method has been applied to the solution of inverse problems. For instance, Liu et al. firstly proposed a Newmark-based method to identify dynamic loads, demonstrating smaller cumulative errors and a superior recognition accuracy compared to the traditional state-space method [38]. Jiang et al. developed an inverse Newmark-based algorithm for estimating pavement roughness and proposed a dynamic load identification method for continuous systems based on the Newmark-β method [39,40]. Pourzeynali et al. introduced a moving load identification method based on the Newmark-β method and generalized Tikhonov regularisation method and verified the feasibility by numerical and experimental studies [41].

A review of the existing studies reveals that the majority of sensors are typically arranged in critical sections of the bridge to capture the bridge dynamic responses [24–27,30–33,37–41]. However, traditional identification methods exhibit lower efficiency in the actual bridge process, making it challenging to swiftly identify moving loads for numerous small- and medium- span bridges [42–44]. To address this issue, the vehicle scanning method (VSM) proposed by Yang has recently gained widespread attention [45–47]. The method necessitates the placement of sensors in a movable detection vehicle which traverses the bridge under testing, thus capturing vibration data to identify modal parameters indicative of the bridge's health status, such as the frequency, damping, and mode [46]. Otherwise, there are no related studies identifying the moving contact force using VSM. In this context, this study aims to propose a new MFI method inspired by the core idea of the VSM, with sensors installed on the moving vehicle rather than the bridge. Section 2 of this paper introduces the theoretical background of the method. The use of Tikhonov regularisation to address the ill-posed problem is discussed in Section 3. The finite element method for vehicle–bridge interaction analysis and a detailed identification procedure for the proposed method are outlined in Sections 4 and 5. Section 6 describes the employment of a single degree-of-freedom (Dof) vehicle crossing a simply supported bridge to verify the accuracy of the algorithm. Section 7 presents a parametric analysis of several typical external factors, while Section 8 concludes the study.

#### **2. Theoretical Background of the Problem**

The vehicle model is simplified to a single DOF moving mass with a support spring to better clarify the moving load force identification method proposed in this study. The bridge is set according to simple support boundary conditions. The mechanical model of the vehicle moving over the bridge is shown in Figure 1. The vehicle runs at a constant speed of *v*. As depicted in Figure 1, *EI* represents bridge bending stiffness, *mb* represents unit length weight, *L* represents bridge span, *v*(*x*, *t*) represents bridge vertical displacement, and *uc* represents contact displacement at the current time step; for the vehicle model, *mv* represents vehicle mass, *cv* represents the vehicle damping coefficient, and *kv* represents the vehicle stiffness coefficient.

**Figure 1.** Mechanical model of the vehicle–bridge system.

Based on Dalembert's principle, the vehicle dynamics equation can be expressed as:

$$m\_v \ddot{y}\_v(t) + c\_v \dot{y}\_v(t) + k\_v y\_v(t) = P(t) \tag{1}$$

where *<sup>P</sup>*(*t*) is the contact force between the vehicle and bridge; and .. *yv*(t), . *yv*(*t*), and *yv*(t) are the acceleration, velocity, and displacement of the vehicle, respectively.

Based on the Newmark algorithm, the interrelationship between the dynamic responses of moment *ti* and moment *ti*+<sup>1</sup> can be established as:

$$
\dot{y}\_v^{i+1} = \dot{y}\_v^i + (1 - \gamma) \Delta t \ddot{y}\_v^i + \gamma \Delta t \dddot{y}\_v^{i+1} \tag{2}
$$

$$y\_v^{i+1} = y\_v^i + \Delta t \dot{y}\_v^i + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{y}\_v^i + \beta \Delta t^2 \ddot{y}\_v^{i+1} \tag{3}$$

where *γ* = 0.5, *β* = 0.25; and the superscript "*i*" and "*i* + 1" denote the vehicle response at the *ti* and *ti*+<sup>1</sup> moments, respectively. For simplicity, the same, simplified notations will be adopted in the following equations.

By rewriting Equation (3), the acceleration at the *ti*+<sup>1</sup> moment can be obtained as:

$$\ddot{y}\_v^{i+1} = \frac{1}{\beta \Delta t^2} \left( y\_v^{i+1} - y\_v^i \right) - \frac{1}{\beta \Delta t} \dot{y}\_v^i - \left( \frac{1}{2\beta} - 1 \right) \ddot{y}\_v^i \tag{4}$$

Substituting Equation (4) into Equation (2) yields the velocity response at the *ti*+<sup>1</sup> moment, illustrated as:

$$\dot{y}\_v^{i+1} = \frac{\gamma}{\beta \Delta t} \left( y\_v^{i+1} - y\_v^i \right) + \left( 1 - \frac{\gamma}{\beta} \right) \dot{y}\_v^i + \left( 1 - \frac{\gamma}{2\beta} \right) \Delta t \ddot{y}\_v^i \tag{5}$$

At the *ti*+<sup>1</sup> moment, the vehicle's motion equation can be expressed as:

$$m\_v \ddot{y}\_v^{i+1} + c\_v \dot{y}\_v^{i+1} + k\_v y\_v^{i+1} = P(t\_{i+1}) \tag{6}$$

Subsequently, by substituting Equations (4) and (5) into Equation (6) and merging the similar items, Equation (6) can be simplified as:

$$
\hat{\mathbb{K}}y\_v^{i+1} = \hat{P}\_{i+1} \tag{7}
$$

where *K*ˆ and *P*ˆ *<sup>i</sup>*+<sup>1</sup> represent the equivalent stiffness and external load, respectively, as illustrated:

$$
\mathcal{K} = k\_v + \frac{1}{\beta \Delta t^2} m\_v + \frac{\gamma}{\beta \Delta t} c\_v \tag{8}
$$

$$\begin{split} \hat{P}\_{i+1} &= P(t\_{i+1}) + \left[ \frac{1}{\beta \Delta t} m\_{\upsilon} + \frac{\gamma}{\beta \Delta t} c\_{\upsilon} \right] \dot{y}\_{\upsilon}^{i} + \left[ \frac{1}{\beta \Delta t} m\_{\upsilon} + \left( \frac{\gamma}{\beta} - 1 \right) c\_{\upsilon} \right] \dot{\bar{y}}\_{\upsilon}^{i} \\ &+ \left[ \left( \frac{1}{2\beta} - 1 \right) m\_{\upsilon} + \frac{\Delta t}{2} \left( \frac{\gamma}{\beta} - 1 \right) c\_{\upsilon} \right] \dot{\bar{y}}\_{\upsilon}^{i} \end{split} \tag{9}$$

For Equation (7), the displacement response *yi*+<sup>1</sup> *<sup>v</sup>* at the end of the time step can be obtained as:

$$y\_v^{i+1} = A\_0 P(t\_{i+1}) + A\_d y\_v^i + A\_v \dot{y}\_v^i + A\_a \ddot{y}\_v^i \tag{10}$$

with

$$A\_0 = \mathcal{K}^{-1},$$
 
$$(11a)$$

$$A\_d = \mathcal{K}^{-1} \left[ \frac{1}{\beta \Delta t^2} m\_v + \frac{\gamma}{\beta \Delta t} c\_v \right] \tag{11b}$$

$$A\_{\upsilon} = \mathcal{R}^{-1} \left[ \frac{1}{\beta \Delta t} m\_{\upsilon} + \left( \frac{\gamma}{\beta} - 1 \right) c\_{\upsilon} \right],\tag{11c}$$

$$A\_d = \mathcal{K}^{-1} \left[ \left( \frac{1}{2\beta} - 1 \right) m\_\upsilon + \frac{\Delta t}{2} \left( \frac{\gamma}{\beta} - 1 \right) c\_\upsilon \right] \tag{11d}$$

After obtaining the displacement response *<sup>y</sup>i*+<sup>1</sup> *<sup>v</sup>* , the velocity response . *y i*+1 *<sup>v</sup>* at this moment can easily be calculated using Equation (5). Additionally, the constant coefficients matrices in Equation (5) are converted into unit ones in terms of *I* = *K*ˆ−1*K*ˆ. After combining the same terms, Equation (5) can be unified into an expression similar to Equation (10), expressed as: ..

$$\ddot{\boldsymbol{y}}\_{\upsilon}^{i+1} = B\_0 \boldsymbol{P}(t\_{i+1}) + B\_d \boldsymbol{y}\_{\upsilon}^i + B\_{\upsilon} \dot{\boldsymbol{y}}\_{\upsilon}^i + B\_d \ddot{\boldsymbol{y}}\_{\upsilon}^i \tag{12}$$

with

$$B\_0 = \frac{\gamma}{\beta \Delta t} \hat{K}^{-1},\tag{13a}$$

$$B\_d = -\frac{\gamma}{\beta \Delta t} \hat{\mathbf{K}}^{-1} \mathbf{K}\_\prime \tag{13b}$$

$$B\_{\upsilon} = \frac{\gamma}{\beta \Delta t} \hat{K}^{-1} \left[ \left( \frac{\beta \Delta t}{\gamma} - \Delta t \right) k\_{\upsilon} + \frac{1}{\gamma \Delta t} m\_{\upsilon} \right],\tag{13c}$$

$$B\_{a} = \frac{\gamma}{\beta \Delta t} \hat{\mathcal{K}}^{-1} \left[ \left( \frac{\beta \Delta t^{2}}{\gamma} - \frac{\Delta t^{2}}{2} \right) k\_{\text{v}} + \left( \frac{1}{\gamma} - 1 \right) m\_{\text{v}} \right] \tag{13d}$$

Similarly, the acceleration response .. *y i*+1 *<sup>v</sup>* can be simplified to the following expression by substituting Equation (10) into Equation (4), expressed as:

$$\dot{\mathbf{y}}\_{v}^{i+1} = \mathbb{C}\_{0}\mathbf{P}(t\_{i+1}) + \mathbb{C}\_{d}\mathbf{y}\_{v}^{i} + \mathbb{C}\_{v}\dot{\mathbf{y}}\_{v}^{i} + \mathbb{C}\_{d}\ddot{\mathbf{y}}\_{v}^{i} \tag{14}$$

with

$$\mathcal{L}\_0 = \frac{1}{\beta \Delta t^2} \hat{K}^{-1},\tag{15a}$$

$$\mathbf{C}\_d = -\frac{1}{\beta \Delta t^2} \hat{\mathbf{K}}^{-1} \mathbf{K}\_\prime \tag{15b}$$

$$\mathcal{L}\_{\upsilon} = -\frac{1}{\beta \Delta t^2} \hat{K}^{-1} (c\_{\upsilon} + \Delta t k\_{\upsilon}),\tag{15c}$$

$$\mathcal{C}\_{\rm a} = \frac{1}{\beta \Delta t^2} \hat{\mathcal{K}}^{-1} \left[ (\gamma - 1) \Delta t c\_v - \beta \Delta t^2 \left( \frac{1}{2\beta} - 1 \right) k\_v \right] \tag{15d}$$

Combing Equations (10), (12) and (14), one yields:

$$
\begin{bmatrix} \dot{y}\_v^{i+1} \\ \dot{y}\_v^{i+1} \\ \dot{y}\_v^{i+1} \end{bmatrix} = \begin{bmatrix} A\_0 \\ B\_0 \\ C\_0 \end{bmatrix} P(t\_{i+1}) + \begin{bmatrix} A\_d & A\_v & A\_a \\ B\_d & B\_v & B\_a \\ C\_d & C\_v & C\_a \end{bmatrix} \begin{bmatrix} \dot{y}\_v^i \\ \dot{y}\_v^i \\ \dot{y}\_v^i \end{bmatrix} \tag{16}
$$

Equation (16) clearly indicates that a recursion formula for the vehicle's response between the *ti* and *ti*+<sup>1</sup> moments is established. On this basis, the displacement, velocity, and acceleration responses of the vehicle at the *ti* moment are written using the above recursion formula as:

$$
\begin{bmatrix} \dot{y}\_{\upsilon}^{i} \\ \dot{y}\_{\upsilon}^{i} \\ \dot{y}\_{\upsilon}^{i} \end{bmatrix} = \sum\_{j=0}^{i-1} \begin{bmatrix} A\_{d} & A\_{\upsilon} & A\_{d} \\ B\_{d} & B\_{\upsilon} & B\_{d} \\ \mathbb{C}\_{d} & \mathbb{C}\_{\upsilon} & \mathbb{C}\_{d} \end{bmatrix}^{j} \begin{bmatrix} A\_{0} \\ B\_{0} \\ \mathbb{C}\_{0} \end{bmatrix} P\_{i} + \begin{bmatrix} A\_{d} & A\_{\upsilon} & A\_{d} \\ B\_{d} & B\_{\upsilon} & B\_{d} \\ \mathbb{C}\_{d} & \mathbb{C}\_{\upsilon} & \mathbb{C}\_{d} \end{bmatrix}^{i} \begin{bmatrix} \dot{y}\_{\upsilon}^{0} \\ \dot{y}\_{\upsilon}^{0} \\ \dot{y}\_{\upsilon}^{0} \\ \dot{y}\_{\upsilon}^{0} \end{bmatrix} \tag{17}
$$

where *y*<sup>0</sup> *<sup>v</sup>*, . *y* 0 *<sup>v</sup>*, and .. *y* 0 *<sup>v</sup>* are the initial displacement, velocity, and acceleration responses, respectively.

In the context of this study, the acceleration transducers are installed on the test vehicle. The vehicle's acceleration response is recorded first, followed by the application of the proposed algorithm to retrieve the moving loads through the recorded data. Supposing the output response of the test vehicle is represented by vector **<sup>y</sup>**(*t*) <sup>∈</sup> **<sup>R</sup>***m*×1, and *<sup>m</sup>* is the number of measured responses,

$$Y = R\_d \ddot{y}\_v + R\_v \dot{y}\_v + R\_d y = \begin{bmatrix} R\_d \\ & R\_v \\ & & R\_d \end{bmatrix} \begin{bmatrix} y\_v \\ \dot{y}\_v \\ \dot{y}\_v \end{bmatrix} = R \begin{bmatrix} y\_v \\ \dot{y}\_v \\ \dot{y}\_v \end{bmatrix} \tag{18}$$

where **R** =*diag*! *Rd Rv Ra* " . Given the installation of a single acceleration transducer on the test vehicle, *Ra* = 1, *Rd* = *Rv* = 0.

Substituting Equation (17) into Equation (18), one yields the following equation:

$$Y(t\_i) = \sum\_{j=0}^{i-1} \mathbf{R} \begin{bmatrix} A\_d & A\_v & A\_d \\ B\_d & B\_v & B\_d \\ \mathbb{C}\_d & \mathbb{C}\_v & \mathbb{C}\_d \end{bmatrix}^j \begin{bmatrix} A\_0 \\ B\_0 \\ \mathbb{C}\_0 \end{bmatrix} P\_i + \mathbf{R} \begin{bmatrix} A\_d & A\_v & A\_d \\ B\_d & B\_v & B\_d \\ \mathbb{C}\_d & \mathbb{C}\_v & \mathbb{C}\_u \end{bmatrix}^i \begin{bmatrix} \mathcal{Y}\_v^0 \\ \mathcal{Y}\_v^0 \\ \mathcal{Y}\_v^0 \end{bmatrix} \tag{19}$$

Assuming the null initial conditions, i.e., *y*<sup>0</sup> *<sup>v</sup>* <sup>=</sup> . *y* 0 *<sup>v</sup>* <sup>=</sup> .. *y* 0 *<sup>v</sup>* = 0, a new matrix **H***<sup>k</sup>* is defined as:

$$\mathbf{H}\_k = \mathbf{R} \begin{bmatrix} A\_d & A\_v & A\_a \\ B\_d & B\_v & B\_a \\ C\_d & C\_v & C\_a \end{bmatrix}^k \begin{bmatrix} A\_0 \\ B\_0 \\ C\_0 \end{bmatrix} \tag{20}$$

Using the matrix **H***k*, Equation (19) can be written in matrix convolution form within the time interval [*t*<sup>1</sup> *tn*], expressed as:

$$\mathbf{Y} = \mathbf{H}\_L \mathbf{P} \tag{21}$$

where

$$\mathbf{Y} = \begin{Bmatrix} \mathbf{y}(t\_1) \\ \mathbf{y}(t\_2) \\ \vdots \\ \mathbf{y}(t\_n) \end{Bmatrix},\tag{22a}$$

$$\mathbf{H}\_{L} = \begin{bmatrix} \mathbf{H}\_{0} & 0 & \dots & 0 \\ \mathbf{H}\_{1} & \mathbf{H}\_{0} & \dots & 0 \\ \dots & \dots & \dots & \dots \\ \mathbf{H}\_{n-1} & \mathbf{H}\_{n-2} & \dots & \mathbf{H}\_{0} \end{bmatrix} \tag{22b}$$

$$\mathbf{P} = \begin{cases} \mathbf{P}(t\_1) \\ \mathbf{P}(t\_2) \\ \dots \\ \mathbf{P}(t\_n) \end{cases} \tag{22c}$$

Equation (21) can be directly solved to acquire the load vector **P** as:

$$\mathbf{P} = \left(\mathbf{H}\_L^T \mathbf{H}\_L\right)^{-1} \mathbf{H}\_L^T \mathbf{Y} \tag{23}$$

Given that the linear correlation between the columns of matrix **H***<sup>L</sup>* is large, the determinant of **H***<sup>T</sup> <sup>L</sup>***H***<sup>L</sup>* is close to zero, implying that the direct solution for Equation (23) becomes an unsettled problem. In other words, the load vector cannot be calculated using Equation (23). To overcome this obstacle, the Tikhonov regularisation approach is utilised in the subsequent section.

#### **3. Tikhonov Regularisation**

In practice, the equation used for identifying a moving load often manifests illposedness due to the coefficient matrix. The ill-posedness problem is exhibited as a loss of rank and the presence of pathology; thus, the paramount task lies in solving these issues for the force identification. When data collected at a high sampling frequency from measurement points exceed the number of moving loads, the problem of rank loss no longer exists. Consequently, the critical aspect for identifying the moving load centres on solving the pathological problem. To this end, this study employs Tikhonov regularisation to solve the instability encountered during the process of moving load identification.

As depicted in Equation (23), the solution for identifying the moving load is **P** =- **H***<sup>T</sup> <sup>L</sup>***H***<sup>L</sup>* −<sup>1</sup> **H***<sup>T</sup> <sup>L</sup>***Y**. Nevertheless, the inversion process for **<sup>H</sup>***<sup>T</sup> <sup>L</sup>***H***<sup>L</sup>* could become unstable due to the significant impacts of the numerous coefficient conditions and external noise in the measured data. As a result of such instabilities, the solution of the load vector **P** deviates greatly from the actual vehicle axle load. In addressing the pathological problem, the Tikhonov regularisation method applies a weak smoothness constraint and selects the approximate solutions of the equation. The actual process for Tikhonov regularisation thus involves solving the minimum value of the following equation:

$$
\hat{\mathbf{y}} \text{ min} = \left\{ \| \mathbf{H}\_L \mathbf{P} - \mathbf{Y} \|\_2 + \lambda^2 \| \mathbf{P} \|\_2 \right\} (\lambda \ge 0) \tag{24}
$$

where · <sup>2</sup> signifies the two-dimensional norm of the vector, and *λ* signifies a regularisation parameter. From Equation (24), it can be observed that the Tikhonov regularisation method actually serves to find the minimum of the residual two-dimensional norm **H***L***P** − **Y** 2 and the two-dimensional norm of the solution.

The L-curve principle is employed to obtain the optimal solution and to identify the optimal regular parameter *λ*. As a raph-based parameter selection method, the L-curve principle involves plotting two quantities, the solution norm log **P***λ* <sup>2</sup> (y axis) and the residual norm log **H***L***P***<sup>λ</sup>* − **Y** <sup>2</sup> (x axis), on the logarithmic coordinate scale to find the equilibrium point, illustrated in Figure 2. At the corner point of the L-curve, both the solution norm and the residual norm are maintained in a balanced dynamic equilibrium, indicating that the point represents the optimal regularisation parameter value.

**Figure 2.** L-curve.

#### **4. Finite Element Method for Vehicle–Bridge Interaction Analysis**

A vehicle–bridge interaction (VBI) element with road surface roughness r was used, illustrated in Figure 3. The test vehicle is modelled as a sprung mass *mv* supported by a spring of stiffness *kv* and a dashpot of damping coefficient *cv*. The coordinate *xc* represents the position of the contact point in the element local coordinate. The equation of motion for the VBI element can be expressed as:

 *mv* 0 **0** [*mb*] \* .. *y* + *v*.. *ub* , - + *cv* <sup>−</sup>*cv*{*N*}*<sup>T</sup> c* <sup>−</sup>*cv*{*N*}*<sup>c</sup>* [*cb*] <sup>+</sup> *cv*{*N*}*c*{*N*}*<sup>T</sup> c* .\* . *y* + *v*. *ub* , - + *kv* <sup>−</sup>*kv*{*N*}*<sup>T</sup> <sup>c</sup>* − *vcv*{*N* }*T c* −*kv*{*N*}*<sup>c</sup>* [*kb*] + *vcv*{*N*}*c*{*N* }*T <sup>c</sup>* <sup>+</sup> *kv*{*N*}*c*{*N*}*<sup>T</sup> c* .\* *yv* {*ub*} - = \* *cvvr <sup>c</sup>* + *kvrc* −(*cvvr <sup>c</sup>* + *kvrc* + *mvg*){*N*}*<sup>c</sup>* - (25)

where *yv* and {*ub*} represent the displacement of the vehicle and that of the bridge element, respectively, and a dot represents the derivative with respect to time *t*; [*mb*], [*cb*], and [*kb*] denote the mass, damping, and stiffness matrices of the bridge element, respectively; *rc* represents the road surface roughness at the contact point *xc*; and {*N*}*<sup>c</sup>* represents the cubic Hermitian function evaluated at the contact point *xc*, whose detailed expression is based on the literature.

**Figure 3.** VBI element.

The global equation of motion for the entire bridge system is formulated by assembling all the bridge elements (with and without the vehicle) and imposing the boundary conditions, i.e., zero vertical displacement and moments, at the two ends. Subsequently, the Newmark-*β* method, with *β* = 0.25 and *γ* = 0.5, is adopted for the response to the coupled dynamic equation.

#### **5. Identifying Procedures for the Proposed Method**

The proposed identification method primarily includes the following three steps:

**Step 1**: The motion equation of the vehicle–bridge coupling system under a moving test vehicle is established. The vehicle's displacement, velocity, and acceleration responses are obtained via numerical calculation.

**Step 2**: *Ai*, *Bi*, and *Ci* (*i* = 0, *d*, *v*, and *a*) are calculated based on the properties of the test vehicle. According to the sensor type and its layout position on the vehicle, the matrix is determined, and the vector **Y** and matrix **H***<sup>k</sup>* are calculated using the vehicle response obtained in Step 1.

**Step 3**: Substituting the vector **Y** and matrix **H***<sup>k</sup>* obtained in Step 2 into Equation (23) and using the Tikhonov regularisation method to solve the ill-conditioned problem, the moving force as the test vehicle traverses on the bridge can finally be identified.

#### **6. Numerical Verifications**

In order to authenticate the reliability of the proposed method, a numerical verification is conducted within this section. This is achieved by employing the identification procedure delineated in Section 5. The contact response, obtained through numerical calculation, is taken as the *real* value, whereas the contact response identified in reverse is taken as the *predicted* one. The accuracy of the proposed method is verified by comparing these two values. The properties of the vehicle and bridge are listed in Table 1. The vehicle traverses the bridge at a constant speed of 10 m/s. A random road roughness sample with Grade A is generated based on the ISO 8606 specification, with the simulated road roughness sample illustrated in Figure 4. It is assumed that the vehicle maintains a zero initial state, with a step size of 0.001 s adopted for the numerical calculation. Employing the vehicle and bridge parameters as listed in Table 1, the vehicle acceleration response is demonstrated in Figure 5.


**Table 1.** Properties of the vehicle and bridge.

**Figure 4.** Random road roughness sample.

**Figure 5.** Acceleration response of the test vehicle.

To evaluate the predicted accuracy of the proposed method, the root mean square error (RMSE) between the predicted value and the real one is used herein, calculated as:

$$\text{RMSE}(y, \hat{y}) = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (y\_i - \hat{y}\_i)^2} \tag{26}$$

where *yi* and *y*ˆ*<sup>i</sup>* denote the predicted and real contact forces at time step *i*, respectively, and *n* stands for the total number of time steps. As derived from the definition of the RMSE, the value is smaller, while the predicted accuracy is higher for the proposed method.

Figure 6 presents a comparison of the moving load identification results between the real and predicted values. It can be observed from the figure that the results for the predicated moving load align well with the real ones, despite the influence of road roughness. The RMSE value is 0.1379, being close to null for the case considered, implying that the proposed method has a high accuracy in identifying the moving contact force. Given the aforementioned comparison results, the reliability of the proposed method is verified.

**Figure 6.** Comparison results of the contact force for the vehicle–bridge system.

Furthermore, the conventional time domain method (CTM) is utilised to further validate the feasibility of the proposed method. Differing from the method presented here, the CTM involves installing sensors on the bridge to record the bridge acceleration responses as the vehicle passes over it. Two different sensor deployments are considered for the CTM, referred to as Case 1 and Case 2. The specific deployments for both cases are listed in Table 2.



Figure 7 depicts the comparison results of the moving contact force for various identified methods. It can be observed that the contact force identified using the proposed method coincides closely with Case 2 using the CTM and the real values, whereas a distinct disparity is observed with the prediction result of Case 1. In comparison with the CTM, the proposed method utilises fewer sensors to obtain the moving contact force with high precision. However, a greater number of sensors must be adopted to ensure the identification accuracy of the CTM. Consequently, the comparison results further validate the accuracy and advantages of the proposed method.

**Figure 7.** Comparison results of the contact force for various identified methods.

#### **7. Parametric Analyses of Various Influencing Factors**

The accuracy of the proposed method was verified in the preceding section. This section further studies the robustness of the method, performing a sensitivity analysis of various external factors, i.e., road roughness, running speed, and noise level, that may affect the recognition accuracy.

#### *7.1. Effect of Road Roughness*

Using the simulation method outlined in the ISO 8606 specification, random road roughness samples for Grades B and C were generated, as illustrated in Figure 8. It becomes evident that pavement classes B and C possess greater amplitudes and higher frequencies in comparison to Class A. Figure 9 showcases a comparison of the moving force identification results at various pavement roughness levels. It can be observed that the vehicle–bridge contact force gradually rises with the increasing road surface grade. Notwithstanding the heightened roughness of classes B and C, the moving load identification employing the proposed method remains consistent with the real one. The RMSE values for Grades B and C stand at 0.1380 and 0.1384, respectively. Comparatively, the identified accuracy for these cases is slightly reduced compared to the result for Grade A, as depicted in Section 6, yet continues to present high precision. Consequently, the proposed method exhibits minimal susceptibility to road roughness in identifying the moving load. Furthermore,

higher pavement grades lead to a more significant acceleration response, favouring the identification of the moving contact force.

**Figure 8.** Road roughness samples for various class levels.

**Figure 9.** Comparison results of contact force for the vehicle–bridge system at various class levels: (**a**) Class B. (**b**) Class C.
