Abstract
Track geometry is an important parameter reflecting the safety of bridge railway, which is of great significance for the structure control of railway track. The measurement of track geometry is an indispensable part of structural health monitoring for bridge railway. However, existing methods for track static track geometry measurement are inadequate in terms of accuracy and efficiency. In order to improve the measurement method, this paper proposes a rapid and precise measurement method based on IMU/odometer aided by multiple total stations, which has the advantages of high efficiency, high accuracy, and is less affected by the environment. The proposed method uses a track measurement trolley (TMT) integrated with IMU and an odometer to collect railway track data on the bridge, and multiple total stations are used to synchronously measure the three-dimensional coordinates of the TMT. After coordinate system registration and time synchronisation, the data are then fused using a filtering algorithm to finally obtain the static geometry of the railway track. The proposed method was validated on the NanSha Port Railway Bridge across the Xijiang River in Foshan, China, and the results showed that the proposed method can make a single on-orbit measurement of a 1000 m length of large-span bridge railway track within half an hour, with an accuracy on the order of millimetres. The proposed method has obvious advantages in terms of efficiency and accuracy and is therefore an effective new measurement method with practical significance.
1. Introduction
Bridges are an important component of railway track. In recent years, with the popularity of railways, more and more large-span railway bridges have emerged, and the demand for structural health monitoring of the railway track located on bridges is also growing rapidly. With the increase in the length of railway bridges, the difficulty of ensuring the stability and safety of the bridges also increases, which in turn increases the difficulty of maintaining railway track safety. The static geometry of railway track refers to the shape and three-dimensional coordinates of the track, which can be used to analyse the safety index parameters of the track such as superelevation, alignment, irregularity, and twist. Track geometry is a very important type of data in track safety operation and maintenance, which has a significant impact on the travelling speed, smoothness, operation safety, and service life of the wheel-rail system [1, 2]. Therefore, the measurement of static geometry of railway track on a large-span bridge is important for the structural health monitoring of railway tracks, which is an urgent problem to be solved.
Measurement of the static geometry of the railway track on large-span bridges faces several challenges: (1) Measurement speed problem. The shape of railway track on a large-span bridge can change on the order of centimetres due to disturbances caused by temperature, wind, and other external factors. It has been shown that the maximum displacement of the bridge reaches 226 mm due to temperature differences between day and night [3]. So, the measurements need to be fast enough, so as to avoid data invalidation caused by significant shape change of the railway track itself. (2) Intervisibility problem. The slope of railway bridge is generally between 2‰ and 4‰ and presents an arch-shaped profile, resulting in the fact that the maximum height difference of the bridge is more than 1 metre. For example, for a 1100 m long arch-shaped bridge with a slope of 3.5‰, the maximum height difference is 1.925 m. So, the intervisibility from one end of a large-span bridge to the other will likely not be guaranteed, and the measurement range of the optical instrumentation may not cover the entire shape of the bridge’s railroad track. (3) Accuracy problem. Under the precondition of fast enough measurement speed, it is also necessary to control the measurement accuracy on the millimetre scale. In the middle of the bridge far away from the ground, the traditional measurement method of track geometry selects the CPIII (track control network of China high speed railway) control point on the bridge as the data. However, the location of the CPIII control points will change with the change of the bridge shape, which leads to a decrease in the measurement accuracy of the track geometry. In conclusion, the measurement method needs to balance both accuracy and efficiency, and there is an urgent need for a rapid and precise measurement method to solve the above difficulties in measuring the static geometry of railway track on large-span bridge structures.
Currently, there are multiple methods for measuring the geometry of bridge railway track, but they are not entirely suitable for large-span bridge railway track. One method is to use camera-based visual sensors for long-term observation and calculate the bridge track geometry based on photos [4–6], which has advantages such as simple installation, no need for special targets, and flexible measurement locations. However, this method requires good lighting conditions and has a relatively short effective measurement distance, only suitable for small-span bridges.
The track inspection vehicles are a mature method for measuring the geometry of the track and are divided into large inspection trains for dynamic track geometry [7–10] and portable inspection trolleys for track static geometry [11–15]. The former has the advantages of high speed, high efficiency, and comprehensive inspection items, and the maximum operating speed reaches 380 km/h. The latter carries out absolute positioning by total station (TS), which is light, flexible, and convenient, with positioning accuracy up to 1-2 mm, and can reflect the track geometry more truly.
However, there are still some inadequacies of the rail inspection vehicle for the railway track of large-span bridges. Large inspection trains are usually used for the dynamic measurement of railway track and cannot be the measurement means for the track static geometry of large-span bridge railway. In addition, the large inspection train is expensive with small number, and the inspection time needs to be planned in advance, which makes it difficult to work whenever the track needs to be inspected. Rail inspection trolley has the following inadequacies. (1) The rail inspection trolley needs to reset up every about 70 m. Each station needs to observe 6∼8 CPIII control points, and the trolley has to be stationary at each rail sleeper to measure 10 seconds. Therefore, the track inspection trolley working speed is only about 150 m/h [1], which means that only about 300 m of track geometry can be measured in 2 h skylight time. (2) The absolute position accuracy of trolley is highly dependent on the CPIII control points. However, the CPIII points will move with the bridge shape change, and the relative point position changes up to 6 mm [16], which is difficult to guarantee the positioning accuracy of trolley. (3) The frequent resetting up leads to “overlap error” [17–19], which affects the accuracy of the final results.
In response to the shortcomings of the track inspection trolley based on total station positioning, GNSS-based positioning trolley has been proposed [20–22]. It has the advantage of not requiring frequent resetting up, thus improving the measurement efficiency and having a relative accuracy of better than 1 mm. However, the absolute accuracy of this method is limited by the GNSS positioning accuracy, which is generally in the centimetre scale, and the GNSS signal is easily interfered or obscured in complex scenes such as bridges, tunnels, and valleys [23, 24], so GNSS-based positioning trolley is not fully applicable to large-span bridge railway tracks.
From the above analysis, it is clear that the current method is difficult to fully adapt to the rapid and precise measurement of the static geometry of railway tracks on large-span bridges, and a measurement method with high efficiency, high accuracy, and low influence by the environment is needed.
In this paper, we propose an IMU/odometer-based method aided by multiple total stations for rapid and precise static geometry measurement of large-span bridge railway track. Specifically, we use a track measurement trolley (TMT) integrated with high-accuracy IMU and an odometer to collect track shape data along the railway, while using multiple total stations to synchronously observe the TMT to obtain its 3D coordinate data. Then, we register the coordinate systems and synchronise the multisensor data. After fusion, we can obtain the static geometry of the railway track.
The advantages of the method are the following. First, the TMT can collect track geometry data with high accuracy and frequency. Second, multiple total stations observe the TMT in sections to achieve intervisibility everywhere, and the absolute positioning of the trolley can be observed throughout the whole process. Third, after coordinate system registration and time synchronisation, the track static geometry is obtained by multisensor data, which realises inertial navigation positioning aided by multiple total stations and odometer, thus effectively suppressing the divergence error of IMU.
On-site test shows that the IMU/odometer-based measurement method aided by multiple total stations proposed in this paper has good feasibility and is able to make a single on-orbit measurement of a 1000 m length of large-span bridge railway track within half an hour, with an accuracy on the order of millimetres. This method has obvious advantages in terms of accuracy and efficiency and is an effective new type of measurement means, which is of good practical significance for the structural health monitoring and safety maintenance of bridge railway tracks.
2. Materials and Methods
The flowchart of IMU/odometer-based static geometry measurement method aided by multiple total stations is shown in Figure 1. In Figure 1, the dashed rectangles represent the major parts of the proposed method, the solid rectangles represent the data obtained from measurement or calculation, the rounded rectangles represent the instrument hardware for collecting the data, and the notched rectangles represent the data processing algorithm.
The overall process is divided into two major parts: data acquisition and data calculation. In the data acquisition process, the track measurement trolley (TMT) collects data related to the track alignment, where the inertial measurement unit (IMU) collects the acceleration and angular velocity of the track measurement trolley, the odometer collects the speed, and the total stations observe the 3D coordinates of the trolley. In the data calculation process, we first register the coordinate systems of multiple total stations, and then synchronise the time of multiple total stations and IMU/odometers. We fuse the data from multiple sensors with Kalman fusing filter and Rauch–Tung–Striebel (RTS) filter. In the end, we obtain the track static geometry.
2.1. Data Acquisition
2.1.1. Hardware System
The hardware used in the method of this paper mainly includes track measurement trolley and total stations. We use TS60 total station with tracking measurement function, which has 2″ setup accuracy with electronic level. The continuous tracking accuracy of TS60 is 0.5″ in terms of angle and 3 mm + 1.5 ppm in terms of distance. For a 1000 m long railway bridge, a single total station tracks TMT halfway, whose point observation accuracy is at least 3 mm + 1.5 ppm 0.5 km = 3.75 mm. It meets the required accuracy of the method. For multiple total stations that have been set up and levelled, the maximum vertical axis difference between them is 4 seconds. If two total stations are placed at the starting point and the ending point of a large-span railway bridge, the maximum total station lap error caused by the vertical axis angle is 0.5 km sin (4″) = 9.7 mm in the elevation direction. Therefore, it is necessary to perform coordinate system registration for multiple total stations to unify the coordinate systems of each total station.
The TMT, shown in Figure 2(a), is equipped with a combined IMU/odometer measurement system (white cylindrical machine in the figure) developed by Shenzhen University. The IMU has an accuracy of 0.01°/h, and the odometer encoder is mounted on the wheels of the TMT and connected to the system. A column is applied to place a prism to obtain the 3D coordinates of the trolley itself.
(a)
(b)
2.1.2. On-Site Data Acquisition
Bridge expansion joints are located at both ends of the bridge and divide the bridge into three parts, named the bridge main body between the expansion joints and the two approach bridges outside the expansion joints. The bridge expansion joint has the function of buffering and absorbing the displacement of the bridge main body, so it is considered that the vibration of the bridge main body caused by wind, water flow, trolley movement, etc., has negligible effect on the approach bridges.
Before collecting data, the TMT and total station were placed on top of the bridge deck, and the position relationship diagram is shown in Figure 2(b). The total stations are placed on the approach bridges near the expansion joints and set up through the CPIII control points on the approach bridges, thus avoiding the CPIII location errors caused by the shaking of the bridge main body. The TMT is placed on top of the track in the bridge main body, and the prism of the trolley needs to ensure the intervisibility. While the trolley is moving along the track and collecting track shape data, one total station observes the prism to obtain the 3D coordinates of the trolley.
The data acquisition process is shown in Figure 3, where Mirror rotating location 1/2 is the position where the trolley is exactly incompatible with total station TS 01/02 during the pushing process. The initial position of the TMT is on the railway track of bridge main body near the TS01, with the prism facing TS01. We turn on the data acquisition instrument and push the trolley along the track toward TS02 on the opposite bridge end. At the Mirror rotating location 1, we change the direction of the prism so that it faces TS02 and continue to push in the original direction. When the trolley reaches the expansion joint near TS02, i.e., the return point, we keep the direction of the prism unchanged, turn the direction of the trolley, and push it toward TS01. When the trolley reaches Mirror rotating location 2, we change the direction of the prism again so that it faces TS01 and continue to push the trolley until it returns to the initial position.
From Figure 3, it can be seen that the overall motion route of the TMT is a round trip of track alignment, performing a repeat measurement of track alignment. The two rotating mirror positions divide the motion route of the trolley into three parts. The process from the initial position Expansion joint 1 to Mirror rotating location 1 is named Run01, the process from Mirror rotating location 1 to Mirror rotating location 2 is named Run02, and the process from Mirror rotating location 2 to the initial position Expansion joint 1 is named Run03. Obviously, in Run01 and Run03 parts, TS01 observes the 3D coordinates of the prism of the trolley, and in Run02 part, TS02 observes.
2.2. Data Calculation
2.2.1. Coordinate System Registration
The 3D coordinate data of the TMT come from several different total stations. Although the total station has been set up and levelled before the measurement, the vertical axes of the total stations are not completely parallel to each other, and there are errors in the CPIII points used for setting up the station, thus the coordinate systems of multiple total stations are not unified. Therefore, the coordinate systems of multiple total stations need to be registered, so that the respective coordinate systems can be unified into the same coordinate system. Coordinate system registration for total stations is similar to point cloud registration, whose principle is to minimise the square sum of distances to the corresponding points between two point clouds by rotation and translation [25–27].
The bridge railway track is generally a straight line in the horizontal plane, while it shows an undulating shape in the elevation direction. Therefore, we establish the coordinate system by taking the track as the X-axis, the mileage as the X-coordinate, the elevation direction as the Z-axis, and the direction perpendicular to the track as the Y-axis. The shape of a series of 3D points acquired by the total station is a curve extending along the track. The mileage of the curve is corrected and normalised so that the mileage of points measured by multiple total stations can be unified. Therefore, we believe that multiple total station data do not need to be adjusted in the mileage, i.e., X coordinate, between multiple total stations, so we take X axis as the reference and adjust the error in Y and Z-axes directions to complete the coordinate system registration for total stations.
The focus of point cloud registration is to determine the relationship of corresponding points between two point clouds. Since the rotation angle is very small, the mileage change of the points on the point cloud curve after rotation is also very small. So, we take the part of the total station point cloud curve with overlapping mileage as the area to be aligned, and the points with the same mileage as the corresponding points. Specifically, the mileage range of overlapping area is , the point cloud of the reference curve located in the overlap region is , the point cloud of the target curve located in the overlap region is , and are a point in the point cloud and , respectively. and share the same mileage coordinates , is elevation, and the subscript denotes the number of pairs of corresponding points in the overlapping region.
After the above preprocessing of mileage correction, mileage normalisation, and corresponding point relationship determination, the coordinate system registration for the total stations can be performed. We take the mileage-elevation plane coordinate system, i.e., X-Z plane, as an example, the coordinate system registration consists of two steps. The first step is to rotate, as shown in Figure 4, the reference total station point cloud data are selected, and the target point cloud is rotated counter-clockwise by an angle around its own total station centre point , so that the two point cloud curves with the same mileage are parallelled. The variance of the elevation residuals (i.e., the distance in the point cloud registration) of the corresponding points is minimised. The formulas are as follows:where and are the mileage coordinate and elevation coordinate of the point in the rotated target curve , respectively. and are the elevation residuals between the corresponding points after rotation and their average values. is the variance of the elevation residuals. The value is obtained by solving the above equation, and then, the target point cloud curve is rotated according to the value.
The second step is translation, as shown in Figure 5. The target point cloud that has been rotated in the previous step is translated up or down along the elevation direction with the translation distance . The purpose is to minimise the square sum of the elevation residuals of the corresponding points in the overlapping part. The formulas are as follows:
After the above rotation and translation, the target point cloud curve and the reference point cloud curve are in the best matching state, as shown in Figure 6. The two curves in the overlapping area coincide and are both in the coordinate system of the reference point cloud curve, thus completing the coordinate system registration for the total stations in the X-Z plane.
The coordinate system registration for total stations in the mile and direction perpendicular to the track plane, i.e., the X-Y plane, is similar to that of the X-Z plane. However, the difference is that no rotation is required. It is enough to translate the target point cloud curve in the Y-axis direction.
2.2.2. Time Ssynchronisation
Due to the long distance apart, it is difficult to connect the total stations to the IMU/odometer on the TMT by wired cable or by wireless WLAN, and there is no hardware-level time synchronisation solution [28, 29]. We adopt synchronisation algorithm to unify the time systems of the total stations and IMU/odometer.
Generally, the IMU/odometer system on the TMT will be turned on during the whole measurement time, with continuous and complete time stamp. While each total station only observes a part of the route, its time stamp is discontinuous, only accounting for a part of the whole measurement time. Therefore, the IMU/odometer system is chosen as the reference base to adjust the time stamp of each total station for synchronisation.
The time synchronisation method is a three-step process. First, the velocity of the prism motion is calculated based on the 3D coordinates and time stamps of the points observed by the total stations. Then, we adjust the time stamp of the prism velocity profile, slide it over the IMU/odometer system time stamp, as shown in the left panel of Figure 7, and calculate the Pearson correlation coefficient between the prism velocity and the odometer velocity in the same time interval. We obtain the function between the time adjustment amount and the correlation coefficient . Finally, the best time adjustment amount corresponding to the maximum correlation coefficient is taken, as shown in the right panel of Figure 7. The Pearson correlation coefficient is calculated using the following formula:where and are the velocity calculated by the total station and the odometer velocity, respectively, and are the standard deviations of them. The total station time stamp is adjusted by an overall amount , thus completing the time synchronisation between this total station and the IMU/odometer system.
2.2.3. Multisource Data Fusion
The data from IMU, odometer, and total station are fused by Kalman filtering [2, 28, 30–32] to obtain the static geometry of the track, i.e., the three position coordinates of the track and the track shape. The system state model, error model, and observation model of Kalman filter are expressed as follows:where and are the system state vector and its covariance matrix at the previous moment, respectively. and are the system state vector and its covariance matrix at the present moment, respectively. and are the state transition matrix and measurement matrix, respectively. is the measurement vector at the present moment. and are the process noise vector and observation noise vector, both following normal distributions with mean 0 and covariance matrices and , respectively.
According to the actual situation, the system state matrix and the observation matrix are designed as shown in following formulas. contains all kinds of errors of the IMU/odometer system and is the position error of the prism.where are attitude error vector, velocity error vector, position error vector, residual gyroscope bias, residual accelerometer bias, and scale coefficient error of odometer, respectively. are the predicted and measured control point coordinates of prism, respectively, is the predicted position of IMU, is estimated attitude matrix, is position vector of the prism with respect to the IMU, is measurement value of the prism from total station, and is the measurement noise of total station.
In order to obtain the system state and covariance at this moment, we calculate the Kalman gain matrix from the system state matrix and its covariance matrix at the previous moment. After continuous forward iteration, the system state of the whole track route can be calculated.
The data are backward filtered using RTS algorithm [33–36]. The principle of RTS algorithm is similar to Kalman filtering, but the computation is in the opposite direction. The RTS algorithm computes the gain matrix from the state and covariance , at the next moment and computes , at the present moment. The formulas are as follows:
3. Results and Discussion
The experimental object of the method in this paper is the NanSha Port Railway Bridge across the Xijiang River in Foshan, China, as shown in Figure 8, which is located at the junction of Shunde District, Foshan City, and Pengjiang District, Jiangmen City, Guangdong Province. The bridge belongs to the double-tower, double-rope surface steel box hybrid girder cable-stayed bridge, with a total length of 1117.5 m, as shown in Figure 9. The two main towers of the bridge are 208 m and 200 m high, respectively. The bridge is the world’s largest cable-stayed bridge with the largest structural span in a two-line freight railway. Our team, Shenzhen University’s dynamic and precise engineering surveying team, measured the railway track geometry of this bridge from October 24, 2022, to October 26, 2022, with a daily measurement skylight time of 21:00∼23:00.
3.1. Results and Analysis of Coordinate System Registration
The data obtained from two total stations at both ends of the bridge are shown in Figure 10. Since the railway track is a straight line running east to west, we set the east direction, north direction, and elevation direction as X, Y, and Z direction, respectively. From Figure 10, it can be calculated that the total length of track alignment is 1112 m, and the overlapping part of Run01 and Run02 has a mileage length of 191 m.
Before and after registration, the residuals and their trend-fitting lines of the overlapping parts of Run01 and Run02 are shown in Figures 11 and 12, in elevation direction and in northward direction, respectively. To show the effect of the registration, we perform a statistical analysis for the residuals of the overlapping parts of Run01 and Run02, as shown in Tables 1 and 2. We choose these statistical parameters for comparison: maximum absolute value, average, standard deviation of the residuals, and the slope of trend-fitting straight line of the residuals, where maximum absolute value, average, and standard deviation are used to describe how far apart the two curves are after registration, and the slope is used to describe how parallel the two curves are.
As can be seen from the table, in the elevation direction and in the north direction (perpendicular to the track direction), there are significant decreases in the maximum absolute value, average, and standard deviation after registration, and the slope of trend-fitting straight line tends to be closer to 0. It can be seen that the coordinate system registration for total stations not only improves the problem of excessive residuals between corresponding points but also corrects the non-horizontal trend of the residuals to a horizontal trend floating up and down around 0.
Combining with the mileage length, the relative accuracy of the overlapping part of Run01 and Run02 in the elevation direction is 6.1 mm/191 m = 1/31311.5 = 31.9 ppm (parts per million). The relative accuracy in the north direction is 7.7 mm/191 m = 1/24805.2 = 40.3 ppm. Therefore, after the coordinate system registration for the total stations, the point clouds of the two total stations are in the same coordinate system and have high registration accuracy.
3.2. Results and Analysis of Time Ssynchronisation
Figures 13 and 14 show the results before and after the time synchronisation of TS01 and TS02, respectively. It can be found that the vibration amplitude of the velocity calculated by the total station is larger than that of the IMU/odometer system. This is due to the fact that the sampling rate of the total station is about 10 Hz, lower than the IMU/odometer. Nevertheless, it can be seen from the figure that both IMU/odometer and total station have the same trend of velocity change.
(a)
(b)
(a)
(b)
After time synchronisation, the best correlation coefficient between TS01 velocity and IMU/odometer velocity is 0.94, and the best correlation coefficient between TS02 velocity and IMU/odometer velocity is 0.92. Therefore, it has good effect that the time stamps of the two total stations have been adjusted to best match with the IMU/odometer system.
3.3. Results and Analysis of Data Fusion
The result of the data fusion calculation of IMU, odometer and total stations, i.e., the track static geometry is shown in Figure 15, from which the track shape and track 3D coordinates can be obtained.
The total mileage of the round-trip measurement is 2224.4 m, and the total time spent is 3215.5 seconds, equal to approximately 54 min. It means that a single on-orbit measurement time is 27 minutes, and the speed is 2471.5 m/h, resulting in a significant increase in measurement efficiency. On the other hand, the sampling rate of the IMU was set to 500 Hz, and a total of more than 1.12 million sample points were acquired. This implies an increase in the data acquisition efficiency, which leads to a more coherent and detailed description of the track geometry.
We analyse the accuracy of the repeated part of Run02 itself, i.e., the part of Run02 from Mirror rotating location 1 to Expansion joint 2. The repeat part is divided into two parts, the former and the latter, using Expansion joint 2 as the dividing point. The elevation and north residuals between the corresponding points with the same mileage of the two parts are calculated and analysed, as shown in Figure 16.
(a)
(b)
The mileage of the repeat part of Run02 itself is 428.5 m. To show the effect of the data fusion, we perform a statistical analysis for the residuals of the repeat part of Run02 itself, as shown in Tables 3 and 4. We choose these statistical parameters for comparison: maximum absolute value, average, standard deviation, and the relative accuracy of the residuals. As can be seen from the table, there is a significant decrease in all these statistical parameters. It is concluded that the fused track geometry data has a high accuracy.
4. Conclusions and Prospects
The track static geometry measurement is important for the structural health monitoring of railway track on the large-span bridge. Existing methods have deficiencies and cannot fully meet the needs of fast and precise measurement. In this paper, we propose a rapid and precise measurement method based on IMU/odometer aided by multiple total stations for large-span bridge railway track static geometry. After coordinate system registration, time synchronisation, filter fusion, the track geometry is obtained.
The advantages of this method include that it can make a single on-orbit measurement of a 1000 m length of large-span bridge railway track within half an hour, with an accuracy on the order of millimetres. Meanwhile, the method does not rely on GNSS signals or unstable CPIII control points located in the bridge main body. Consequently, compared with other track geometry measurement methods based on total station, level and camera-based visual sensor, the method proposed in this paper has better measurement performance and dependability.
The method in this paper has been implemented on the NanSha Port Railway Bridge across the Xijiang River in Foshan, successfully obtaining track shape, track 3D coordinates and other information. The results show that the measurement efficiency can reach more than 2400 m/h. The maximum error of the track geometry is −4.5 mm in the elevation direction and −3.9 mm in the direction perpendicular to the track. The relative measurement accuracy of the track geometry is 2.6 ppm in the elevation direction and 2.1 ppm in the direction perpendicular to the track. On-site test shows that the accuracy, efficiency, sampling rate, and reliability of this method have obvious advantages compared with the existing methods and is of good practical significance for structural health monitoring and safety maintenance of bridge railway track.
In the future, the method in this paper should be developed in the direction of higher speed under the precondition of millimetre scale accuracy. Higher measurement speed can further reduce the disturbances caused by temperature, wind, and other factors on the shape of the bridge and improve the measurement validity.
Data Availability
The nature of the data is the fact that our team worked with the company that owns the bridge railway to obtain the results. The data that support the findings of this study can be accessed from the first author, Shiwang LV, upon reasonable request. The data are not owned by our team alone, but is related to the company that owns the bridge railway. So there are restrictions on data access.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (U1934215) and Young Scientist Fund of National Natural Science Foundation of China (41901412).