Research on Seismic Vulnerability of High-Pier and Long-Span Bridges Based on Improved IMK Resilience Model

After the fragility curve is established, the probability of structural damage reaching each level of damage under the action of the ground motion can be determined according to the ground motion parameters, so as to calculate the direct and indirect loss caused by the structural damage and complete the earthquake damage prediction. This paper combines the improved IMK resilience model to study the seismic vulnerability of high-pier and long-span bridges. Moreover, this paper obtains the parameter calculation model based on the regression analysis of PEER’s 255 column specimen data. The improved IMK model needs to modify the elastic stiffness and strain hardening rate of the rotating spring to ensure the accuracy of the lateral stiffness of the component. The experimental research shows that the seismic vulnerability research model of high-pier and long-span bridges based on the improved IMK restoring force model has a certain analytical effect.


Introduction
The seismic design strength of existing ordinary small and medium-span bridges is generally based on the national unified seismic zoning. The seismic design parameters of the same structural type of bridges in the same zoning are the same, so the structural strength reserve when encountering an earthquake is the same [1]. However, during the actual earthquake, the seismic intensity of various places within the same seismic zone is not consistent or even significantly different, which leads to bridges with insufficient redundancy of individual structures that are prone to damage [2]. Therefore, it is necessary to use the probabilistic seismic hazard analysis method considering the temporal and spatial inhomogeneity of seismic activity to obtain the regional seismic hazard analysis results from the perspective of seismic subdivision. At the same time, it is necessary to combine dynamic time history analysis and structural vulnerability analysis to calculate the failure probability of the bridge under a given earthquake [3]. On this basis, the connectivity reliability of the traffic network after the earthquake is calcu-lated, the vulnerability of the traffic network under the influence of the earthquake is studied, and a calculation model for solving the optimal path is proposed. In order to further improve the calculation efficiency and reduce the amount of calculation storage, the Target Order algorithm is proposed and the effectiveness of the algorithm is verified in the actual road network. The data obtained from the research can provide references for future engineering construction planning in the region. At the same time, predicting the connectivity reliability of the transportation network after the earthquake will help the government to be effective and timely in making daily maintenance decisions and postearthquake disaster relief decisions. As a result, the loss caused by the earthquake is reduced as much as possible, and the overall earthquake prevention and disaster reduction capabilities of the region are improved.
Over the past few decades, bridge engineering, which is an important foundation for national development, has developed significantly. With the gradual advancement of the national infrastructure construction process, the number and span of newly built bridges have been greatly improved.
At the same time, relevant seismic codes for bridge engineering are gradually improving. However, bridge damage caused by frequent earthquakes is becoming more and more common. Once the bridge is damaged, the transportation network may be paralyzed; in the earthquake rescue and relief, the bridge acts as a lifeline, and the consequences of its damage will be more serious. The reasonable implementation of bridge seismic fortification needs to be established on the basis of reasonable bridge seismic damage prediction model. On the one hand, a reasonable bridge earthquake damage prediction model can make a more accurate assessment of local earthquake losses from a macro perspective, including preearthquake assessment to control possible losses within an acceptable range, and establish a reasonable earthquake early warning and monitoring mechanism. It can quickly formulate rescue plans after an earthquake to provide support for government decision-making; on the other hand, the bridge earthquake damage prediction model can also provide a more accurate reference for the seismic design of specific bridges. Although my country has carried out a lot of research on earthquake damage prediction models of highway bridges, the researches are usually only for single bridges, there are few studies on the seismic damage assessment of bridge networks, and there is still a lack of relatively unified highway bridge earthquake damage prediction models. This restricts the development of earthquake damage prediction systems for highway bridges to a certain extent. In contrast, some developed countries in the world (such as the United States) have established standardized earthquake damage prediction methods. The most typical representative is HAZUS, which was jointly promoted by FEMA and NIBS in 1997. In HAZUS, it is believed that the key to earthquake damage prediction of highway bridges is to classify bridges based on their seismic performance and define five failure states of bridges, including no damage, slight damage, medium damage, severe damage, and complete damage, and then be able to determine the vulnerability curve of each type of bridge in each failure state. The seismic vulnerability curve is a function of ground motion parameters as an independent variable (defined as a cumulative distribution function that obeys a lognormal distribution in HAZUS). This paper combines the improved IMK resilience model to study the seismic vulnerability of high-pier and long-span bridges and builds an intelligent model to provide a theoretical reference for the subsequent performance improvement of high-pier and long-span bridges.

Related Work
Literature [4] used expert investigation to analyze the seismic vulnerability of highway bridges and formed the seismic vulnerability curve of highway bridges. However, the limited number of samples limits the credibility of the research results. Literature [5] studied the contribution rate of various structural features to the seismic vulnerability of bridge structures and formed seismic vulnerability curves of railway bridges under the action of earthquakes of different intensities. The method proposed in [6] forms an earthquake vul-nerability curve. It is assumed that each seismic vulnerability curve is a standard log-normal cumulative distribution function with unknown location parameters and known proportional parameters, so as to consider both capacity and demand. It is assume that each seismic vulnerability curve is a standard log-normal cumulative distribution function with unknown location parameters and known constant scale parameters, so as to integrate the random uncertainties of both capacity and demand. Literature [7] formed the seismic vulnerability curve of bridge piers by using the bridge damage data observed in the Kobe earthquake. In the research, the log-normal distribution function of two parameters is used to characterize the seismic vulnerability curve, and the relevant parameters are estimated by the maximum probability method. Literature [8] formed the seismic vulnerability curve of the bridge structure based on the bridge damage data obtained in the Northridge earthquake using logistic regression analysis. The seismic vulnerability curve formed by the empirical analysis method is usually based on the actual seismic damage data and the corresponding statistical results of ground motion parameters and has high reliability. The damage curve is only suitable for situations similar to the data source. Therefore, the empirical vulnerability curve is difficult to promote and use, and it is necessary to establish a reliable theoretical analysis method for the seismic vulnerability curve of bridge structures.
When analyzing the seismic vulnerability of bridge structures, different researchers have adopted different analysis methods. Reference [9] combines the characteristics of the bridge structure to systematically study the seismic vulnerability of the pier column at the most vulnerable part of the regular beam bridge and proposes an analysis based on the numerical simulation to form the seismic vulnerability curve of the pier column method. Literature [10] formed seismic vulnerability curves of reinforced concrete piers based on traditional reliability methods and numerical simulation methods and compared the seismic vulnerability curves obtained by the two methods. Literature [11] uses the Monte Carlo sampling simulation method to form the seismic vulnerability curve of the bridge structure. The response of the structure is calculated by two different methods, namely, the nonlinear dynamic time history analysis method and the capacity spectrum method. The comparison results of the seismic vulnerability curves formed by the two methods show that the consistency of the seismic vulnerability curves of the bridge in the two states of severe damage and complete failure is not as good as that of the slightly damaged state. Literature [12] proposed a method of forming seismic vulnerability curve through numerical simulation. The bridge pier is idealized as a single-degreeof-freedom vibration system, the seismic records of the Hanshin earthquake and the Northridge earthquake are used as input to calculate the Park-Ang damage index of the pier, and then the Park-Ang damage index and ground motion strength measurement parameters are used to form the reinforced concrete pier Earthquake vulnerability curve. Literature [13] analyzed the seismic vulnerability of a four-span simply supported girder bridge, obtained the seismic response data of the structure through nonlinear dynamic 2 Journal of Sensors time history analysis, and used a logical model to determine that the structure surpassed the specified ground motion intensity measurement parameters and used a logic model to determine the conditional probability that the structure exceeded the specified ground motion intensity measurement parameters. In this method, whether a certain vector in the Bernoulli random variable is 1 or 0 is determined according to whether the bridge maintains a certain damage state, and the parameters in the model are determined by regression analysis. Literature [14] formed the seismic vulnerability curves of four typical bridges through a comprehensive analysis of bridges in the central and southeastern regions of the United States. In the literature [15], the nonlinear analysis model and artificial seismic records were used to first form the seismic vulnerability curve of each bridge component (pier or support), and then the first order reliability principle was used to obtain the seismic vulnerability of the entire bridge system and then the seismic vulnerability damage curve of the bridge is obtained by using the firstorder reliability principle. The research results show that among the four types of bridges, the most prone to damage are multispan simply supported beam bridges and multispan continuous steel beam bridges, and the least likely to be damaged are multispan continuous prestressed concrete beam bridges.

Improved Ibarra-Medina-Krawinkler Restoring Model
The hybrid simulation of seismic response of reinforced concrete frame structures and the analysis of their collapse resistance require a restoring force model that can effectively simulate the hysteresis characteristics of reinforced concrete members and accurate model parameters. The plastic dumpling model is a model often used in the nonlinear simulation of frame structures. It not only reflects the mechanical characteristics of the member but also is closely related to the material, restraint condition, and spatial layout of the member. Since decades, the plastic dumpling model has developed rapidly. The improved IMK model, as a plastic dumpling model with a trifold as the skeleton curve, introduces a degradation parameter β based on energy dissipation, which can consider a variety of degradation characteristics of the member under reciprocal loading. Compared with other plastic dumpling models, it can better simulate the hysteresis characteristics of reinforced concrete beams and columns effectively [16]. The main model skeleton parameters include elastic stiffness K e , yield moment M y , plastic corner θ cap,pl , postyield hardening stiffness M c /M y , and postpeak corner θ pc . An important feature of this model is the postpeak unloading section [17]. It can simulate the strain softening behavior associated with physical phenomena such as concrete crushing, steel buckling fracture, and bond damage. Figure 1 shows the skeleton curve of the modified IMK model. The improved IMK model includes three types of models, that is, the bilinear model, the peak pointing model, and the pinched model.
The bilinear model is based on the bilinear hysteresis rule with strain hardening, which retains the main part of the bilinear hysteresis rule and introduces the degradation section after the peak and the residual strength section. The bilinear model is shown in Figure 2.
After unloading from point 3, the model sets the strength at point 3 as the strength limit in the positive direction. When the load is unloaded from point 5 in the negative direction, the maximum value of the restoring force is the previously set strength limit value, that is, the vertical coordinate of point 6. The load is then loaded from point 6 to point 3 with О stiffness and then continues with the degraded segment stiffness.
The peak pointing model retains the basic hysteresis rule and adds the degraded section and the residual strength section after the peak, and the model is shown in Figure 3 [18].
Compared with the bilinear model, the peak pointing model is unloaded at the positive 2 points and then loaded from the point on the horizontal axis to the previously unloaded 2 points after a hysteresis half-loop, that is, 8 points and 2 points are the same point, without considering the strength degradation. When unloading, the reloading stiffness changes after each arrival at the horizontal axis, and its value is determined by the intersection point with the horizontal axis and the peak displacement point in the previous cycle together.
The pinch model is similar to the peak pointing model; however, its reloading phase consists of two parts, as shown in Figure 4. In the reloading process, it is first loaded to the breakpoint 8 and then loaded to the previous unloading point 2. The breakpoint is determined by the maximum permanent deformation and the maximum force in the direction of loading. The parameter k d is introduced in determining the transverse coordinate of the breakpoint, and the parameter k f is introduced in determining the longitudinal coordinate of the breakpoint. If the deformation δ per at the intersection of the reloading section with the transverse coordinate is greater than ð1 − k d Þ/δ per , and if strength 3 Journal of Sensors degradation is not considered, then the reloading path is loaded only for the straight line section from point 13 to the previous unloading point 10 [19].
The improved IMK model can depict four basic cyclic degradation modes of the member, that is, basic strength degradation, softening segment strength degradation, unloading stiffness degradation, and reloading stiffness degradation. The peak pointing model is used as an example to introduce the four degradation modes of the improved IMK model.
The cyclic degradation rate β is determined by the energy dissipation of the member under reciprocal loading. It assumes a reference value for the energy dissipation capacity of each member and does not consider the effect of the member loading history. The cyclic degradation rate β is calculated as follows [20]: Among them, β i denotes the degradation coefficient calculated in the cycle numbered i, E i denotes the reference value of the energy consumption of the component, E i = γ F y δ y , and γ is the coefficient calibrated from the experimental results. For each degradation mode is different, E i refers to the energy consumption of the component in the cycle numbered i, the second bias E i refers to the energy consumption of the component in the cycle numbered i − 1, c denotes the exponential term coefficient of the degradation rate, and c is taken from 1.0 to 2.0. If the loading cycle consists of a sequence of cycles of constant amplitude, then c = l indicates an almost constant degradation rate, while c = 2 indicates a slowing down of the degradation rate of the earlier cycles and an acceleration of the degradation rate of the later cycles [21].
Throughout the loading history, the value of β must be taken to satisfy 0 ≤ β i ≤ 1, that is    Journal of Sensors After each half turn of the cycle, the yield strength of the member will be degraded accordingly, and the degradation of the yield strength satisfies the expression as follows: Among them, F +/− i denotes the yield strength value after the cycle numbered i in positive and negative directions, and F +/− i=1 denotes the yield strength value before the cycle numbered i.
For each degradation parameter β s,i , which contains values in the positive direction and values in the negative direction, the basic intensity degradation rates in these two directions are calculated independently. That is, the value of F − i is updated after the end of each cycle in the positive direction, and the value of F + i is updated after the end of each cycle in the negative direction.
The degradation mode of the basic strength also includes the degradation of the hardened stiffness, and the form of the hardened stiffness degradation is as follows: Similar to the yield strength, the hardening stiffness is calculated independently in both positive and negative directions. Figure 5 illustrates the degradation pattern of the basic strength. After the loading path is loaded from 0 to 3, the degradation rate β s,i is calculated, at which time the negative yield strength of the member degrades from F − y to F − 1 and the negative stiffness of the member degrades from K s,0 to K s,1 − . After the loading path is loaded from 3 to 7, the value of the degradation rate β s,i is recalculated. At this time, the positive yield strength of the member is degraded from F + y to F + 1 , and the positive stiffness of the member is degraded from K s,0 to K s,1 + . Unlike the basic strength degradation, the softening segment strength degradation means that after the peak point, the slope of the softening segment remains constant. However, the intensity value of the intersection of the extension line and the vertical coordinate degrades with the number of cycles. The softening segment strength degradation equation is shown below: F +/− ref is the intensity value at the intersection of the extension of the softened segment with the vertical coordinate after the peak, and similarly, the positive and negative degradation rates β e,i are calculated independently. The softening segment strength degradation rate is calculated each time it intersects the horizontal axis, which may not affect the loading path in the early stages of the nonlinearity. Figure 6 shows the softening section strength degradation. This shows that the softening segment strength degradation rate is first calculated when loading to point 3, where the intersection of the softening segment extension with the vertical coordinate in the negative direction degrades from At loading up to 6, the softening segment strength degradation rate is recalculated, and the intersection of the positive softening segment extension with the longitudinal coordinate degrades from F + ref ,0 to F + ref ,1 . The unloading section contains negative unloading section and positive unloading section, the unloading stiffness degradation does not distinguish between positive and negative direction, and its stiffness degradation is calculated by the following formula: K u,i and K u,i−1 denote the unloading stiffness before and after the half-turn of the cycle numbered i, respectively.

Journal of Sensors
Among them, the zone is calculated by equation (1) and the appropriate cyclic degradation parameter γ k . Unlike the degradation rates of the remaining three cases, the unloading stiffness degradation rate is the only one that is calculated when the load is reversed in the nonlinear phase. In addition, the unloaded stiffness is updated in both positive and negative directions. Therefore, the degradation of the unloaded stiffness is twice as high as that of the other degradation types. Figure 7 illustrates the degradation of unloading stiffness. The unloading stiffness degrades from K e to K u,1 when loading to point 2 and then continues to degrade from K u,1 to K u,2 when loading to point 5.
The reloading stiffness degradation increases the magnitude of the target displacement. Depending on the direction of loading, the increased target displacement is defined by the maximum displacement of the last corresponding directional cycle. The reload stiffness degradation is only applicable to the peak pointing model and pinched model with the following expressions: Figure 8 shows an illustration of the reload stiffness degradation. The value of the reload stiffness degradation rate β a,i is updated each time it is unloaded to intersect with the horizontal axis. The cyclic half-turn target displacement δ +/− i,j−1 with the corresponding direction number i is calculated from the cyclic halfturn target displacement δ +/− i,j with that direction number i − 1. The improved IMK model requires the determination of seven parameters, and the seven parameters are initial stiffness K e , yield moment M y , plastic turning angle θ cap,pl , postyield hardening stiffness M c /M y , postpeak turning angle θ pc , normalized hysteretic energy dissipation parameter λ, and exponential term degradation parameter c.
The empirical prediction equations for the above parameters were established based on 255 column specimens in the PEER structural performance database, and the empirical prediction equations for the seven parameters are shown below.
At present, there is no good definition for the initial stiffness of steel composite members. Two forms are used to quantify the effective stiffness, that is, the cut-line stiffness value EI y past the yield point, and the cut-line stiffness EI suf 40 past the 40% yield value point, and the empirical prediction formulas for both forms are as follows: (1) The cut line stiffness past the yield point is calculated as follows: (3) The formula for calculating the secant stiffness at the point of over 40% yield value is as follows: Among them, P/A g f c ′ is the axial compression ratio of the column, and L s /H is the aspect ratio of the column.
The empirical prediction formula of yield bending moment, the prediction formula of bending strength, and the formula are calculated as follows: Among them, the yield curvature is the smaller value in equations (15) and (16).
(1) If the yield of the section is determined by the yield of the tension steel bar, the calculation formula of the yield curvature is as follows: (2) If the yield of the section is due to the significant nonlinear phenomenon of the compressed part of the concrete, that is, the strain of the compressed part of the concrete exceeds its limit ε e ≈ 1:8f c ′ /E c , then, the calculation formula of the yield curvature is calculated as follows: The depth k y of the compressed part when the section yields (normalized by d) is Among them, n = E s /E c . If the reason for the section yielding is formula (1), then, the calculation formulas of A and B are formula (18) and (19). If the cause of section yielding is equation (2), then, the calculation formulas for A and B are equations (20) and (21).
In the formula, ρ, ρ′, and ρ v , respectively, represent the reinforcement ratios of the longitudinal reinforcement of the compression part, the longitudinal reinforcement of the tension part, and the web reinforcement. δ ′ = d ′ /d, d ′ refers to the distance from the center of the steel bar in the compression part to the outer edge of the concrete in the compression part. b and d are the width and height of the section, respectively, and N refers to the axial load (compression is positive).
Studies have shown that the axial compression ratio and concrete strength are the key factors that determine the hardening stiffness. The hardening stiffness M/M is expressed as follows: The simplified constant equation for hardening stiffness after yielding is as follows: Among them, v is the axial compression ratio, f c ′ is the concrete axial compressive strength, and c units is the unit conversion variable. If the unit of f c ′ is MPa, then c units takes 1. If the unit of f c ′ is ksi, then c units takes 6.9.
The complete calculation formula of the plastic turning angle is as follows: The simplified calculation formula of the plastic turning angle is as follows:  (1 or 0), v represents the axial compression ratio of the member, s n is the buckling coefficient of the reinforcement, s n = ðs/d b Þ ð f y /100Þ 0:5 , d b is the diameter of the longitudinal reinforcement, f y is the yield strength of longitudinal bars, ρ sh is the area ratio of transverse steel bars in the column plastic dumpling zone, and ρ is the ratio of longitudinal bars.
Equations (24) and (25) are only applicable to the case where the cross-section is symmetrical. For the case of asymmetrical reinforcement, the analysis item of asymmetrical reinforcement is introduced into the formula, and the calculation formula of plastic rotation angle that can be used for asymmetrical reinforcement is obtained: Among them, v is the axial compression ratio, and ρ sh is the transverse steel area ratio. The improved IMK model contains four degradation modes, namely, basic strength degradation, softening section strength degradation, unloading stiffness degradation, and reloading stiffness degradation. Each degradation mode is defined by two parameters: the normalized capacity dissipation capacity (λ) and the exponential term of the cyclic degradation speed (c). In order to reduce the complexity, assuming that the two parameters of the four cases are the same, the cycle degradation parameters are reduced from 8 to 2. The calculation formula 29 is a simplified calculation formula of a, and the corresponding c is taken as 1.0.
The cyclic energy dissipation capacity a is incorporated into OpenSees in the form of an input parameter, and the calculation formula of ∧ is as follows: a is the parameter of normalized energy dissipation capacity, which is defined by the total energy of the component: In the above formula, when the model adopts the initial stiffness instead of EI y /EI g , the value of λ should be adjusted.
When the residual strength is 0, it sometimes causes errors in the OpenSees program. In order to avoid errors, the residual intensity can be taken as a nonzero positive number (0.01 is enough), so that errors in the program can be avoided in most cases.
For the flexural member that may deform to the inelastic range, a rotating spring of zero length is simulated at the end of the member, and the rotating spring at the end of the member (simulated by the zero-length unit) and the upper elastic unit are modeled in series. The overall stiffness of the member is a combination of the stiffness of the rotating spring and the stiffness of the elastic unit, and the expression is as follows: Among them, K mem is the stiffness of the member, K bc is the stiffness of the elastic unit, and K s is the stiffness of the rotating spring. Analyzing the above equation, we can see that when the stiffness of the elastic unit is much larger than the stiffness of the rotating spring, the stiffness of the member is approximately the stiffness of the rotating spring. When the stiffness of the elastic beam-column unit is much smaller than the stiffness of the rotating spring, the member stiffness is approximately the stiffness of the elastic beamcolumn unit. In Ibarra's analysis, n times the stiffness of the elastic beam-column unit is used to represent the stiffness of the rotating spring, that is, K s = nK bc . Bringing K a = nK bc into equations (33)-(35) are obtained. To make the stiffness of the member equal to the stiffness of the elastic unit, the cross-sectional moment of inertia of the elastic unit is the original cross-sectional moment of inertia of the elastic unit.
After correcting the stiffness of the rotating spring, the strain hardening ratio, which is the ratio of the hardened stiffness to the elastic stiffness, needs to be corrected as well. The angle of rotation Δθ mem of the member in the plastic phase consists of the angle of rotation of the rotating spring and the angle of rotation of the elastic unit, which is calculated as shown in equation (36). Equation (37)   Journal of Sensors α mem K mem of the member in the plastic phase, the expression (38) for the strain hardening coefficient of the rotating spring can be solved by combining it with equation (37).
Δθ mem = ΔM in 1 + nα s,s nα s,s K bc , ð37Þ Among them, Δθ s denotes the angle of rotation in the plastic phase of the rotating spring, Δθ bc denotes the angle of rotation in the plastic phase of the elastic unit, ΔM in is the increased bending moment in the plastic phase, K s,s and K bc denote the stiffness of the rotating spring and the elastic unit in the plastic phase, respectively, α s,s and α s,mem denote the strain hardening rate of the rotating spring and the member, respectively, and n is the ratio of the rotating spring stiffness K s to the elastic stiffness K bc set in the elastic phase.

Research on Seismic Vulnerability of High-
Pier and Long-Span Bridges Based on Improved IMK Resilience Model Figure 9 shows a high-pier long-span bridge model. Three types of ground motions are applied longitudinally to obtain the bridge response, and the maximum curvature ductility ratio of the control section is obtained. In this paper, PGA is used as the independent variable, and the curvature ductility ratio is the dependent variable, and the IDA curve that controls the section curvature ductility ratio is drawn, as shown in Figures 10 and 11.  It can be clearly seen from Figures 10 and 11 that with the continuous increase of PGA, the maximum curvature ductility ratio of the critical section of the bridge is also increasing, which shows that it is appropriate to use PGA to represent the index of ground motion intensity. Due to the differences in the characteristics of the three types of ground motions, the damage caused is quite different. Compared with the other two ground motions, the Landers wave has a greater seismic response, while the other two ground motions have little difference in response, which is related to the original intensity of the ground motions. Under the action of a single ground motion, with the continuous increase of PGA, the specific damage state of the key section of the bridge can be seen more intuitively from the damage state histogram.
From the above research, it can be seen that the seismic vulnerability research model of high-pier and long-span bridges based on the improved IMK restoring force model has a certain analytical effect.

Conclusion
Through the analysis and research on the seismic damage evolution and vulnerability of the bridge structure, it is possible to evaluate the seismic performance of the bridge, provide a basis for the seismic design of the bridge, provide a theoretical basis for the assessment of bridge damage after the earthquake, and provide a certain theoretical basis for post-disaster decision-makers to formulate emergency plans for earthquake resistance and disaster reduction. However, high-pier long-span continuous rigid-frame bridges have more complex seismic performance and seismic response compared to low-and medium-span bridges, such as traveling wave effects, pile-soil interaction, and various complex nonlinear problems are more obvious. This paper combines the improved IMK resilience model to study the seismic vulnerability of high-piers and long-span bridges and builds an intelligent model. The test results show that the seismic vulnerability research model of high-pier and long-span bridges based on the improved IMK restoring force model has a certain analytical effect.

Data Availability
The labeled dataset used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare no competing interests.   Journal of Sensors