Numerical simulation scheme of jointed rock masses using UAV photogrammetry and a disk-based discontinuous deformation analysis model

1.   Introduction

Unstable rock blocks and their combinations, which are cut using multiple sets of discontinuities, are generally evaluated as hazardous rock masses. Unstable rock masses characterized by a developed discontinuity network are widely distributed on steep slopes or cliffs in the southwestern mountainous areas of China ^[1]. However, owing to the development of China's economy, many large engineering projects are being increasingly conducted in such areas in accordance with government planning. The geological environments in these areas are extremely complex. Frequent landslide and rockfall disasters can result in significant societal and economic losses. If engineering construction is to be conducted in such areas, then preliminary investigations and surveys are necessary. Existing theoretical analyses and experimental results demonstrated the controlling impact of the discontinuity network on the mechanical behavior of rock masses ^[2]. Currently, based on the measurement of the discontinuity characteristics of rock masses, engineers can accurately assess their danger level and implement the necessary measures, such as early warnings or pre-reinforcements.

Regarding the theoretical analysis of jointed rock mass failure, researchers have primarily applied fracture mechanics and damage mechanics comprehensively to propose many strength criteria ^[3,4,5] and constitutive models ^[6,7,8] for assessing the mechanical properties of rock masses, thus promoting the development of rock mass analysis. Physical experiments are the most typical approach for investigating the mechanical properties of rocks. A series of geotechnical unit experiments, including uniaxial compression ^[9,10,11], direct shear ^[12,13], triaxial ^[14,15,16] and Brazilian ^[17,18] tests, were conducted using natural rock specimens or specimens with pre-set joints to investigate various aspects of the characteristics and mechanisms of jointed rock masses. Furthermore, computed tomography scanning ^[14,19,20], three-dimensional (3D) printing ^[17,21,22], acoustic emission ^[13,16,23], and other advanced technologies have been utilized in conjunction with the abovementioned tests to obtain more intuitive and effective test results. However, theoretical and experimental studies that consider the effects of rock discontinuities in terms of scale, irregularity, or complexity may differ significantly from actual discontinuity networks. It is still limited to use theoretical analysis or laboratory tests to accurately reflect the characteristics of jointed rock masses at the actual engineering scale. Therefore, it is of great significance to obtain the actual discontinuity characteristics of jointed rock masses from a more quantitative perspective and introduce them into the analysis.

As more abundant and accurate discontinuity characteristics of actual jointed rock masses must be obtained while ensuring the safety of the investigators, conventional geological exploration methods such as compass and tape measurements ^[24] are no longer adequate. In recent years, remote sensing technologies have developed rapidly and have gradually been applied to field surveys. Among them, flexible and low-cost unmanned aerial vehicle (UAV) digital photogrammetry can not only be used for surveying and monitoring topography and geomorphology ^{[25,26,27,28]}, but also for acquiring discontinuity characteristics related to jointed rock masses. Significant developments have been achieved in the application of UAVs equipped with compact digital cameras to measure basic information related to jointed rock masses, such as discontinuity orientation ^[29,30], discontinuity spacing ^[31] and discontinuity persistence ^[32]. This indicates that the technology presents significant potential and advantages for the acquisition of rock mass images in complex landforms ^[33]. In summary, combined with the current mature UAV photogrammetry, it is possible to quickly obtain high-resolution images of highly steep rock masses and extract the required discontinuity characteristics from a reconstructed model ^[34]. Actual and complete discontinuity networks are more conducive to the follow-up analysis of large rock masses.

As mentioned above, the typically used theoretical analysis methods and laboratory tests present limitations when analyzing engineering-scale jointed rock masses. By contrast, owing to significant enhancements in computing power, numerical simulation methods have been increasingly used in this field. Among the numerous numerical methods, the discontinuous deformation analysis (DDA) ^[35] is suitable for simulating discontinuous mechanical processes, such as friction, collision, and separation between blocks. In addition, it exhibits good computational stability and high accuracy. Moreover, the DDA has been developed and improved in recent years ^[36,37,38], allowing it to better simulate the continuous–discontinuous deformation process of engineering materials. Therefore, the improved DDA method can reveal the characteristics of jointed rock masses more effectively, thereby resulting in a more accurate simulation of the entire process of microcrack propagation and progressive failure in rock masses.

This study focuses on the Baima jointed rock masses located at K158 + 837 on the Chongqing–Huaihua Railway, which comprises several hazardous rock groups. To obtain high-quality rock slope models, we combine the UAV close-range photogrammetry technology with a 3D modeling method. In the workflow, an image of the rock slope was first generated using projection techniques, and the discontinuity characteristics were extracted from the image. Based on the actual size and the discontinuity network depicted in the image, a corresponding disk-based particle model was then established. Finally, the disk-based discontinuous deformation analysis (DDDA) method was used to conduct a numerical simulation analysis. The aforementioned workflow is intended to serve as a reference for acquiring and recognizing comprehensive image information of jointed rock masses and conducting numerical analyses at an engineering scale.

2.   Methodology

2.1.   Reconstruction of jointed rock masses

Photogrammetry enables the rapid acquisition of image data from actual rock masses, which can be used for modeling and subsequent analyses. Reconstructing a rock mass model from high-quality photographs with a certain degree of overlap mainly involves two key algorithms: structure from motion (SfM) and multi-view stereo (MVS) algorithms.

The SfM algorithm first extracts and matches the feature points on each photograph using a scale-invariant feature transformation algorithm, as shown in Figure 1. Subsequently, based on the feature matching data between photograph pairs, the camera orientation and scene structure information are automatically obtained through perspective projection matrix estimation, scene structure restoration, and iterative bundle adjustment procedures ^[39] to complete the sparse point cloud reconstruction of the target rock masses.

Figure 1.  Illustration of 3D spatial reconstruction based on SfM algorithm.

DownLoad: Full-Size Img PowerPoint

The MVS algorithm can process images from multiple viewpoints to identify the corresponding matching points from a photograph sequence to restore 3D information. Based on the sparse point cloud as well as the image and camera parameters, this algorithm can reconstruct more feature points from the image sequence and complete the dense point cloud reconstruction of the target rock masses.

2.2.   Modeling algorithm for jointed rock masses

The algorithm proposed herein for generating jointed rock mass models mainly involves particle filling and bonding in the model area, as well as the introduction of rock discontinuities and their cutting. The algorithm can establish a disk-based particle model with a discontinuity network to satisfy the modeling requirements of relevant numerical simulations.

Firstly, disks are used as the basic elements to generate particle clusters with arbitrary shapes. The particle-filling algorithm used in this study is based on the approach proposed by Cui and O'Sullivan ^[40], as shown in Figure 2(a). To enhance the algorithm performance, further improvements were realized based on the algorithm proposed by Zhang et al. ^[41]. Specifically, the model area was divided into triangular meshes and filled rapidly with non-overlapping disks. Subsequently, complex boundaries and convex or concave corners were processed to generate a disk-based particle model with a high density and coordination number. The improved disk-filling effect is illustrated in Figure 2(b).

Figure 2.  Schematic illustration of disk-based particle filling algorithm: (a) reference algorithm ^[40]; (b) infilling model based on improved algorithm.

DownLoad: Full-Size Img PowerPoint

Since the disks filled in the model area are discrete, no inherent bond strength exists between them, they must be bonded into clusters to simulate geotechnical materials. To achieve this, a similar method to the bonded particle model in the discrete element method (DEM) ^[42,43] was adopted in this study. Using two normal springs and one tangential spring, disks that are in contact with each other were bonded to achieve the required bond strength and the ability to withstand tensile forces as well as shear forces exceeding the maximum static friction. The bonding element formed between the two disks is shown in Figure 3. Using this bonding algorithm, all contacting disks can be bonded and assigned the necessary physical and mechanical parameters to accurately simulate geotechnical materials.

Figure 3.  The bonding element form.

DownLoad: Full-Size Img PowerPoint

To establish jointed rock mass models, the introduction and cutting of discontinuity network are necessary. In the two-dimensional (2D) simulation, the discontinuities are imported as polylines of certain width, including the coordinates of each point of the polyline and the width values of the beginning and end of each segment. Using a linear discontinuity as an example, the operating principle of cutting is briefly introduced. As shown in Figure 4(a), the width of discontinuity line MN is W, the length is L, and the inclination is θ, and the linear equation $Ax + By + C = 0$ can be obtained from the coordinates of M and N. The linear equations of the other boundary lines can be further obtained from the geometric relationship. Prior to cutting, the positional relationship between the disks and the discontinuity line is determined using the obtained linear equations. If the circle center of the disk is within the area encompassed by the discontinuity line width (for the linear discontinuity, the coverage area is rectangle M₁N₁N₂M₂ in Figure 4(a) plus the semicircle area at both ends), the bonding element associated with the disk will be weakened or deleted, thus weakening or eliminating the bond strength between the corresponding disks. This achieves the cutting effect required to simulate 2D discontinuities of a specified width. Figure 4(b) shows the cutting process of a linear discontinuity on a simple rectangular block.

Figure 4.  Schematic illustration of discontinuity cutting principle: (a) linear discontinuity; (b) cutting case.

DownLoad: Full-Size Img PowerPoint

2.3.   DDDA simulation method

Conventional discontinuous numerical methods typically employ polygons as the basic elements, which can result in complex contact problems, particularly for systems containing concave blocks and blocks with holes. However, using disks as basic elements can avoid these problems. Disks are bonded to approximate blocks of arbitrary shapes, rendering DDDA an effective method for simulating various problems ^[44,45]. Geotechnical disasters are typically characterized by the initiation, propagation, and penetration of fractures in materials, thus resulting in a large-scale displacement of the sliding mass. This process can be simulated by bonding and fracturing between particles in numerous particle aggregates using the DDA method ^[37]. In addition, researchers have calibrated the DDDA method using unit tests such as uniaxial compression and Brazilian splitting tests, demonstrating its feasibility and effectiveness in simulating rock failure ^[38]. In general, the DDDA method is highly practical for simulating the prevention and mitigation of geotechnical disasters. This method is based on the Newmark time integration scheme, and the control equation is established at each time step using the principle of minimum potential energy ^[38,44]. Subsequently, the control equation is solved using an implicit solution, and the open-close iteration is realized such that the basic unknown quantity of the disk can be obtained. In this process, the penalty method is employed to manage the interaction between the disks. The small displacements obtained at each time step are accumulated to calculate the final large displacement of the disk.

3.   Imaging and processing of jointed rock masses

3.1.   Engineering overview

This study focuses on the large jointed rock masses located in the upper section of the natural slope above the entrance of the Baima No.1 Tunnel on the Chongqing-Huaihua Railway. This section is characterized by a hard, medium-thick, steep, layered limestone slope with a height of approximately 141 m and an almost vertical gradient. The steep wall strike is NE60° and the rock masses exhibit a well-developed discontinuity network. Additionally, the vegetation on the slope is sparse. Figure 5 shows the location and characteristics of the study area.

Figure 5.  Engineering case overview: (a) UAV view; (b) close-up scenery of the jointed rock masses; (c) Baima No.1 Tunnel Entrance.

DownLoad: Full-Size Img PowerPoint

3.2.   UAV close-range photogrammetry

As shown in Figure 6, the operation flow of the UAV close-range photogrammetry for jointed rock masses involves the following steps: 1) preparation before photography; 2) capturing digital images of the jointed rock masses on-site; 3) using the SfM–MVS workflow to reconstruct dense point clouds from the images; and 4) generating the 3D scene model.

Figure 6.  Operation flow of UAV close-range photogrammetry.

DownLoad: Full-Size Img PowerPoint

The preparation before photography involves site surveys and route planning. On-site surveys are essential for ensuring the capture of high-quality images, and reasonable route planning can effectively achieve this goal. An automatic route photography scheme is generally preferred when a strong global positioning system (GPS) signal is available at the site. A vertical grid-type route survey scheme is typically used for geological masses on steep terrains. However, because the actual rock wall surface of jointed rock masses is not a vertical plane perpendicular to the horizontal plane and many large concave and convex rock masses are present, using only the vertical grid UAV survey scheme may result in large differences in the image resolution at different waypoints due to the varying distances between the UAV and wall during flight. Insufficient overlap of UAV images at the protruding and concave section of rock masses may lead to loopholes, stretching or blurring in the modeling process, which can affect the quality of the model and even cause modeling failure in certain cases. Setting an extremely small distance to ensure precision during the UAV flight may increase the risk of collisions with mountains. Therefore, a close-range photogrammetry scheme that involves maintaining the UAV at a relatively fixed distance from the actual surface of the rock masses during route planning is preferred. This approach ensures consistency in the image resolution across all sections of the rock wall as well as the modeling quality and UAV safety compared to the vertical grid-type UAV survey scheme.

Before obtaining digital images of rock masses using close-range photogrammetry, as shown in Figure 7(a), an approximate model of the entire area where the jointed rock masses are located must be established using a conventional route to facilitate the generation of the close-range route. Next, we denoted an open ground point accessible to the operator in the approximate model using GPS information as the take-off point, used the route planning tool to delimit the target rock masses for the UAV survey, and set the relevant parameters such as the camera focal length, minimum flight safety distance, number of waypoints, and overlap degree of UAV photographs. The overlap degree of the flight distance should ensure that the heading overlap is not less than 70% and the side overlap is not less than 60%. Based on the set parameters, the close-range route and camera orientation at each waypoint can be automatically calculated and generated. Through manual inspection, the inappropriate waypoints in the route can be adjusted, and then the final route can be derived. The specific UAV survey task was conducted by importing the close-range route into the UAV flight control software and identifying it as the flight reference route.

Figure 7.  3D modeling and detection of discontinuity network in jointed rock masses from projected image: (a) acquisition of rough model; (b) model obtained from the close-range photogrammetry; (c) identified and vectorized 2D discontinuity network.

DownLoad: Full-Size Img PowerPoint

After obtaining high-quality UAV photographs, the 3D scene geometry of the jointed rock masses was reconstructed using the SfM-MVS workflow. Their 3D point clouds and textured mesh models were obtained through photograph alignment, parameter calibration, point cloud encryption, mesh generation, texture generation, and other operational steps. The reconstructed 3D model of the jointed rock masses is shown in Figure 7(b).

3.3.   Projection of rock masses

The established high-precision 3D model can be used for various applications. To conduct the 2D numerical analysis, the reconstructed 3D model of jointed rock masses was used to generate an image of the rock slope using a projection operation. Orthogonal projection was performed to project the 3D model in a specified direction to ensure consistency in the rock masses size. The projection direction was selected to be approximately perpendicular to the rock surface, which was the main direction of the cameras along the route. Before projecting the model, it was first adjusted to an appropriate position and attitude through translation and rotation. Subsequently, the projection direction was defined. The projection direction is typically a 3D vector, which was assumed to be ${\boldsymbol{ D}} = ({x_d}, {y_d}, {z_d})$ in this study. The projection matrix was calculated according to the projection direction as follows:

Each point in the 3D model was represented as homogeneous coordinates, and the orthogonal projection matrix was multiplied by the homogeneous coordinates to obtain the projected coordinates:

Finally, the projected coordinates were converted from homogeneous coordinates to 2D coordinates, and the final coordinates of each point in the image were obtained as follows:

where p₁, p₂ and p₄ represent the elements in the first, second and fourth rows of matrix P, respectively.

The resulting image and vectorized discontinuity network are shown in Figure 7(c). Discontinuities in the 2D rock slope were orderly marked through the recognition and analysis of the image. A total of 153 discontinuities were identified and marked.

4.   Rock masses numerical simulation

4.1.   Modeling process

The projected jointed rock slope graph obtained through close-range photogrammetry contains the actual size information of the entire rock masses. The entire model area was initially filled with particles using the modeling algorithm for jointed rock masses. After performing triangulation based on the outer contours of the rock masses, the interior and boundaries of the model were filled with disks. The mesh was subdivided as much as possible while ensuring that the number of disks did not significantly affect the simulation calculation, thus ensuring the quality and compactness of the filled model.

The disk-based model after mesh division and filling is shown in Figure 8(a), (b), respectively, where the number of triangular meshes was 26,938, the number of disks was 41,979, the average radius of the disks was 0.542 m, and the average coordination number was 4. To simulate the rock slope, the disks were bonded into clusters, so that bond strength was reflected between the disks, and the physical and mechanical parameters of the bond element were assigned. This is illustrated in Figure 8(c), where the number of bonding elements was 81,665. Finally, the discontinuities were introduced into the model, which was then cut. Polylines were imported and corresponded one-to-one to the identification results of the discontinuity network. The width was set as the average value of the discontinuity width measured by randomly sampling multiple groups in the 3D reconstruction model. In this simulation, the specific input value for the discontinuity width was 0.4 m. The cutting algorithm was used to cut the model based on the imported discontinuity network, and the final modeling effect of the jointed rock masses is shown in Figure 8(d). By performing the steps above, we created a disk-based particle model that included the discontinuity network, which enabled us to realize preprocessing of the 2D jointed rock slope.

Figure 8.  Preprocessing of the jointed rock masses: (a) triangulation in the model area; (b) disks filling in the model; (c) bonding of the model; (d) formation of the discontinuity network.

DownLoad: Full-Size Img PowerPoint

4.2.   Simulation scheme

Numerical simulations were conducted on the Baima jointed rock masses model using DDDA to simulate the deformation and failure processes of the jointed rock masses under the action of self-weight. Three types of numerical simulation scenarios with different bond strengths were considered, as summarized in Table 1. In Scenario 1, the adopted strength parameter value was based on the relevant test values reported in the field exploration. To enable a comparative analysis, the strength parameter values were decreased and increased appropriately in Scenarios 2 and 3, respectively, based on the strength parameter in Scenario 1.

Table 1.  Physical and mechanical parameters of three numerical simulation scenarios.

Parameter	Scenario 1	Scenario 2	Scenario 3
Material density /(kg·m-3)	2750	2750	2750
Young's modulus /MPa	50,000	50,000	50,000
Shear modulus /MPa	20,000	20,000	20,000
Tensile strength /MPa	4.1	3.1	5.0
Internal friction /(°)	42.5	36.5	45.0
Cohesion /MPa	6.9	5.0	9.0

---

 | Show Table

DownLoad: CSV

4.3.   Simulation results analysis

The simulation results are presented in Figure 9. First, the initiation, propagation, and penetration of cracks in the jointed rock masses occurred gradually around the discontinuity network introduced via preprocessing in all three scenarios. This indicates that the discontinuity network formed during the initial cutting significantly affects the simulation results, highlighting the importance of identifying and extracting rock discontinuities in the early stages of rock mass analysis. Second, Figure 9 also shows a comparison of the simulation results for three scenarios. In Scenario 1, a few rock blocks near the bank collapsed and several relatively obvious crack propagations occurred at the edges of the jointed rock masses under the action of self-weight, as shown in Figure 9(a). If the bond strength in Scenario 1 is weakened, i.e., in the simulation of Scenario 2 (see Figure 9(b)), a large range of rock masses collapsed and slipped, and the area of significant crack growth increased significantly. In Scenario 3, the bond strength between the particles was improved compared with that in Scenario 1. As shown in Figure 9(c), the rock masses were stable, except for a small range of rock block disintegration.

Figure 9.  Numerical simulation results of the Baima jointed rock masses for three scenarios: (a) Scenario 1; (b) Scenario 2; (c) Scenario 3.

DownLoad: Full-Size Img PowerPoint

The development of microcracks varied during the simulation process. In Scenarios 1 and 3, almost all cracks caused by the failure of the bonding element were tensile cracks. After the crack on the far-shore side stopped developing, the crack on the near-shore side continued to expand, resulting in a more intensive crack distribution at the near-shore side. This indicates that, under the combined effect of the terrain where the rock masses were located and the discontinuity network, the rock masses near the bank were more susceptible to collapse due to multiple cracks reaching the critical condition. In Scenario 2, numerous cracks appeared in the upper section of the rock masses, and upper and lower penetrations formed rapidly around the discontinuity network. A wide range of rock masses collapsed, and a slip surface formed on the bank side. In contrast to the findings for Scenarios 1 and 3, a few shear cracks appeared in the lower-right area of the rock masses in Scenario 2. This indicates that shear cracks appeared later than tensile cracks in the simulation process, and that they primarily appeared near the sliding crack surface when the rock masses approached non-stability and failure.

The simulation results were analyzed using a post-processing program to extract the physical displacement and velocity data from different areas of the rock masses. The analysis results show that crack development occurred primarily in the lower area, and the damage to the near-shore and far-shore sides differed significantly. Therefore, three marked blocks (1, 2, and 3) were selected from the near-shore to the far-shore of the rock masses (Figure 10), and the displacement and velocity curves of these blocks were analyzed. Figure 11 shows the displacement changes of the marked blocks for the three scenarios, whereas Figure 12 shows the corresponding velocity changes. In all three scenarios, block 1 indicated a significantly earlier displacement and velocity increase compared to the other two marked blocks. For example, in Scenario 1, at the moment of 5.5 s, the displacement value of block 1 reached 72.55 m, which was much greater than the values of blocks 2 and 3 at this moment (0.018 and 0.015 m, respectively). This is because the small rock block corresponding to block 1 disintegrated and collapsed early or late in all the simulation processes, thus resulting in a parabolic displacement curve and a broken-line velocity curve with multiple slope changes corresponding to collisions with other blocks or terrain lines during the collapse. The falling times of block 1 in Scenarios 1, 2, and 3 were approximately 0.25, 0.1, and 2.3 s, respectively. A comparison of the displacement and velocity changes among the three marked blocks show that the failure of rock masses began at the near-shore side, and that the rock masses distant from the near-shore side were more stable.

Figure 10.  Location description of three marked blocks.

DownLoad: Full-Size Img PowerPoint

Figure 11.  Displacement curve of marked blocks in three scenarios: (a) Scenario 1; (b) Scenario 2; (c) Scenario 3.

DownLoad: Full-Size Img PowerPoint

Figure 12.  Velocity curve of marked blocks in three scenarios: (a) Scenario 1; (b) Scenario 2; (c) Scenario 3.

DownLoad: Full-Size Img PowerPoint

Figure 13 shows the displacement and velocity change curves of block 1 in the three simulation scenarios. The higher the bond strength assigned to the rock masses, the later was the collapse of block 1, resulting in a slower increase in the displacement and velocity. Therefore, under the same simulation time, block 1 in Scenario 2 reached a larger displacement faster than that in Scenario 3. This rule can be similarly derived from the change in velocity, i.e., the latest acceleration occurred in Scenario 3. However, the collision situation during the collapse of block 1 varied significantly among the three scenarios, thus resulting in different degrees of fluctuation in the velocity change curves. Specifically, in Scenario 1, block 1 indicated the most collisions and showed significant acceleration and deceleration from around the fourth second.

Figure 13.  Displacement and velocity change of marked block 1 in three scenarios: (a) displacement curve; (b) velocity curve.

DownLoad: Full-Size Img PowerPoint

The simulation results presented above can be used to further analyze the current state of the field reinforcement of the rock masses. The reduction in bond strength in the simulation can reflect the weakening of jointed rock masses due to weathering and erosion, whereas the improvement in strength can reflect the effect of reinforcement measures. Therefore, based on the comparison results, the necessary reinforcement measures must be implemented for the Baima large jointed rock masses to ensure traffic safety along the Chongqing–Huaihua Railway. This also corresponds to a series of in-situ pre-reinforcement works on the rock masses mentioned in relevant reports.

5.   Conclusions

In this study, a comprehensive workflow for obtaining and analyzing jointed rock masses at the engineering scale was proposed, and the following conclusions were inferred:

1) UAV close-range photogrammetry technology provides high-resolution and comprehensive image information of the Baima jointed rock masses, thus allowing for the identification of 153 discontinuities and the creation of a vector network containing actual size information.

2) Using the proposed jointed rock masses modeling algorithm, the disk-based particle model with a discontinuity network was preprocessed via the vectorized graph of rock masses. The order of magnitude of the disks reached 40,000.

3) DDDA was performed to numerically simulate and analyze the Baima jointed rock masses under multiple operating conditions; consequently, the instability and collapse characteristics of jointed rock masses and the development of microcracks in the rock were revealed. A comparative analysis highlighted the importance of the existing discontinuity networks and the significance of pre-reinforcement engineering.

Although this study successfully covered the entire process from information acquisition to numerical analysis, some issues still remain to be addressed. First, the recognition of discontinuity networks relied mainly on manual identification, which is inefficient and requires a certain level of geological expertise. Therefore, the recognition method should be further optimized by incorporating computer graphics. Additionally, to serve practical engineering applications more effectively, the proposed methods need to be applied to a wider range of high and steep rock masses, and a corresponding case database should be established.

Acknowledgments

This work was supported by the Fundamental Research Funds for the Central Universities of Central South University (No. 2022ZZTS0208), the Scientific and Technological Progress and Innovation Project of Hunan Provincial Department of Transportation (Nos. 202008; 202115), and the Chinese Scholarship Council (CSC, Nos. 202106370143). All financial supports are greatly appreciated.

Conflict of interest

The authors declare there is no conflict of interest.