Abstract
Most of the tunnel projects are related to the national economy and people’s livelihood, and their operational safety is of paramount importance. Tunnel safety accidents or hidden safety hazards often start from subtleties. Therefore, the identification of tunnel cracks is a key part of tunnel safety control. The development of computer vision technology has made it possible for the automatic detection of tunnel cracks. Aiming at the problem of low recognition accuracy of existing crack recognition algorithms, this paper uses an improved homomorphic filtering algorithm to dehaze and clear the collected images according to the characteristics of tunnel images and uses an adaptive median filter to denoise the grayscaled image. The extended difference of Gaussian function is used for edge extraction, and the morphological opening and closing operations are used to remove noise. The breakpoints of the binary image are connected after removing the noise to make the image more in line with the actual situation. Aiming at the identification of tunnel crack types, the block index is proposed and used to distinguish linear cracks and network cracks. Using the histogram-like method to distinguish linear cracks in different directions can well solve the mixed crack situation in an image. Compared with the traditional method, the recognition rate of the new algorithm is increased to 94.5%, which is much higher than the traditional crack recognition algorithm. The average processing time of an image is 5.2 s which is moderate, and the crack type discrimination accuracy is more than 92%. In general, the new algorithm has good prospects for theoretical promotion and high engineering application value.
1. Introduction
Tunnel engineering has complex forces and its structural safety is of great importance. Diseases such as cracks, water leakage, and lining deformation often occur in operating tunnels. If small cracks and deformations cannot be identified in time and treated scientifically, serious disasters or accidents may be induced. Therefore, the rapid and efficient identification and classification of tunnel cracks is essential to ensure the safety of tunnel operation.
The traditional manual inspection method not only consumes a lot of resources, but also due to subjective limitations in detection accuracy, the accuracy of crack judgment and the specific physical characteristics of the crack area (crack location, crack width, and length) cannot be precisely expressed. Combining computer vision and image processing technology to automatically detect tunnel lining cracks has become the main research direction.
The image processing methods of tunnel lining cracks can be divided into two categories: one is the deep learning method based on convolutional neural network and the other is the grayscale image processing method. The method of deep learning has a high accuracy for tunnel lining crack detection. However, it takes a lot of time to establish and train the model because of the need of disease annotation and a large number of images. The traditional grayscale image processing method is to determine the location of the cracks by the sudden change of the gray value. However, the surface of the tunnel lining has oil stains and water leakage, and the image is easily affected by noise, which will affect the accuracy of crack recognition. The shape of the fracture includes horizontal, vertical, oblique, and network. Currently, there is no method that can precisely distinguish cracks containing different directions. Therefore, it is also extremely important to develop relevant algorithms that can accurately classify fracture directions.
Many scholars have done a lot of research on the method of no-reference identification of tunnel cracks and have summarized many satisfactory results. Zhang [1] proposed a road crack detection method based on clustering-minimum spanning tree, using edge detection and image segmentation to jointly extract cracks in the image. Then, the extracted crack images are clustered, and the minimum spanning tree is used for detection. Xue and Yuan [2] put forward the problem of uneven illumination in tunnel images, and there will be differences in brightness and contrast between broken pixels and surrounding background pixels. They proposed a dynamic partition Gaussian crack detection algorithm based on the distribution of projection curves. The threshold and multiscale Gaussian factor are brought into the Gaussian detection model to complete the fracture extraction. Wang et al. [3] used the combination of Prewitt operator and Otsu algorithm for image segmentation, and the morphological algorithm is used to optimize the binary image. Using the different characteristic values of the fracture area and the impurity area to separate the fracture area has high engineering practical value. Wang et al. [4] proposed a new image preprocessing algorithm which combined global and local and a multistage filtering algorithm based on connected regions. The new algorithm has better processing results for crack images that contain a lot of complex noise, uneven illumination, and uneven contrast. Hoang and Nguyen [5] pointed out that if the fractures in an image are thin, the fine details in the image are very important for crack tracking. They constructed an automatic model for pavement crack detection and classification and used SVM (support vector machine), artificial neural network, and random forest to classify images. The results show that the accuracy of SVM is the highest. Based on the characteristics of coherence and directionality of fractures, Jia and Liu [6] designed an Arabic numeral algorithm for rough extraction of fractures. The iterative method is used to obtain the optimal threshold, and finally, the cracks are extracted finely through the morphologically connected domain mark. Fekri-Ershad [7] improved the LBP model to make it resistant to noise and suitable for multiresolution. The practical application of image surface defect detection shows that the proposed method has the advantages of high detection rate, low noise sensitivity, and low computational complexity.
The tunnel lining pictures acquired by the CCD camera generally have problems such as low contrast, poor lighting conditions, many types of noise, and the influence of tunnel surface texture. The existing tunnel lining image processing algorithms cannot take the above problems into consideration at the same time. Therefore, this paper proposes a computer vision-based inspection method for surface cracks in tunnel linings.
2. Image Processing Methods and Steps
The image filtering method is usually used for the preprocessing of the crack image of the tunnel lining surface, and different scholars use different methods to process the crack image of the tunnel lining [8–10]. However, the existing methods have reduced image details and blurred edges of fractures. In this paper, the improved homomorphic filtering is used to dehaze and deshadow the image to make the image clear. After the image is converted to grayscale, adaptive median filtering is used to denoise the image, and the XDoG (eXtended difference of Gaussians) to perform image edge extraction. Take morphology to open and close to remove small spots. Finally, connect the cracks that are still rough after morphological processing. The flowchart of tunnel lining crack extraction is shown in Figure 1.
2.1. Image Preprocessing
2.1.1. Image Dehaze
When the tunnel lining image is collected, due to environmental changes and the machine itself, the obtained image will be blurred than the original. To a certain extent, the details of the image will be lost, which will affect the further recognition of the image. Use defogging algorithm to make the overall image clearer, and the image details are more reflected. The most commonly used dehazing algorithm is histogram equalization. This method is relatively simple, but the results were not ideal. Sometimes grayscale changes will affect the amount of information in the image [11]. The overall performance of the retinex algorithm is good, and scholars have also summarized many satisfactory results through research. For example, scholars Tajeripour and Fekri-Ershad proposed a single-scale retinex algorithm [12], which has low computational complexity and is insensitive to noise, but it is not completely suitable for tunnel crack detection. This paper uses a method based on improved homomorphic filtering to dehaze and enhance the image. On the basis of the existing method, only one parameter needs to be changed to make the image clearer and easier for subsequent processing.
Homomorphic filtering is a processing method that combines image frequency filtering with grayscale transformation. The image illuminance/reflectance model is used as the basis of frequency domain processing, and the image quality is improved by compressing the image brightness range and enhancing the local contrast of the image which makes the image more in line with the nonlinear characteristics of the human eye for brightness [13]. The basic steps of homomorphic filtering are shown in Figure 2. The principle is in [14]. The image is composed of the product of the illumination component and the reflection component , as shown in the following equation:
Take the logarithm of both sides of (1), and then get
Do Fourier transform to (2) to get the following formula
Design filter transfer function to filter . Conventional transfer functions are divided into Gaussian, Butterworth, and exponential. Different filtering methods correspond to different filters. This paper uses Gaussian homomorphic filtering to process equation (3), and the optimized Gaussian transfer function is expressed as follows:
Among them, and are the parameters that control the amplitude range of the filter, the value range is , and the values are usually , according to experience. is the value of the sharpening parameter that controls the slope of the filter function, and the main function is to control the shape of the filter. is the cutoff frequency. This method needs to set the values of the four parameters of and and continue to experiment to obtain the optimal value, which leads to the selection of parameters that takes a lot of time and is inefficient. Aiming at the uneven illumination characteristics of the tunnel lining crack image, a parabolic filtering method that specifically amplifies the high-frequency signal is proposed to filter the image. Only one parameter is required to realize the filter control. The improved filter function form is shown in the following formula:
In the above formula, a is the maximum frequency gain coefficient value, generally taken as (1, 2]. After many tests in this article, the best value is selected as 1.5. d is the maximum frequency increase position, that is, the frequency at this position receives the maximum gain.
Finally, it is inversely transformed to the spatial domain and exponents are taken on both sides of the equation to obtain the relevant spatial filtering results. The improved homomorphic filtering method proposed in this paper is tested and compared with the algorithm based on dark channel prior [15] and the contrast limited adaptive histogram equalization (CLAHE) algorithm [16]. The experimental results are shown in Figure 3.
(a)
(b)
(c)
(d)
It can be obviously seen from the comparison in Figure 3 that the algorithm in this paper has advantages over the dark channel prior algorithm and the CLAHE algorithm. The improved homomorphic filtering algorithm optimizes the details of the image to a large extent, making the original blurry image clearer. The tunnel lining crack image processed by the dark channel prior algorithm is generally brighter, but it highlights useless details such as the concrete surface texture, which will affect the identification of cracks by the edge detection operator later. The CLAHE algorithm makes the overall image darker and reduces the image contrast. Although the sharpness of the image has increased, the overall darkness makes the crack characteristics inconspicuous and easy to mix with other areas. The difficulty of crack identification is increased, and the accuracy of crack identification is reduced.
2.1.2. Adaptive Median Filter Denoising
After defogging the input image, the contrast of the image is enhanced compared with the input image and the cracks appear more obvious. But there are still other noises, such as impulse noise caused by electromagnetic wave interference and noise caused by concrete surface texture. These noises will affect the subsequent steps, which in turn affects the identification of cracks. Therefore, it is necessary to filter the image after defogging. Zhang et al. [9] studied the application of adaptive median filter in tunnel cracks and proved that adaptive median filter has great advantages compared with traditional median filter and mean filter algorithm. It can not only reduce image noise but also protect the edges of the image without blurring the edges. Therefore, this paper uses an adaptive median filter algorithm to remove the noise in the image after grayscale defogging.
2.2. XDoG Extracts Crack Edges
The Laplacian of Gaussian (LoG) proposed by scholars Marr and Hildieth is a good edge detection operator, so it is widely used in the field of computer vision [17]. In the existing tunnel image defect detection methods, due to the particularity of the tunnel crack image, the Sobel operator and the Canny operator are more susceptible to noise than the LoG operator. It will cause more noise in the detection result, which is not conducive to subsequent processing. However, the LoG calculation process is relatively complicated, so it is very important to find a relatively simple calculation method. Research has shown that in many cases the results obtained by difference of Gaussian (DoG) filter and LoG are very similar, and relatively complex LoG filtering can be achieved by DoG approximation with extremely high accuracy [18]. Taking into account the relatively easy implementation of DoG and the advantages of less noise in the calculation results, this paper intends to choose DoG as the edge detector.
2.2.1. Gaussian Differential Filtering
The DoG filter is a change of the Gaussian filter and consists of the difference between two Gaussian functions. The specific expression form is shown in equation (6). σ is the standard deviation, which determines the smoothness of the Gaussian filter. The larger the σ, the better the smoothing effect.
The two-dimensional expression of the Gaussian function is given by
Since the image is a two-dimensional dataset, a two-dimensional Gaussian function expression must be used. The combination of x and y in equation (7) represents the corresponding value in the two-dimensional Gaussian function template. Substituting formula (7) into (6), we can get
2.2.2. Extended Gaussian Difference Function
Extended difference of Gaussian (XDoG) function [19] is an improved edge extraction algorithm based on DoG function. Compared with other tunnel image processing algorithms, it has unparalleled advantages. It can control the overall details of the image and control the image line thickness and overall brightness after edge extraction and other advantages, and the computational complexity is low. The basic calculation process of expanding the Gaussian difference function is to introduce the coefficient Υ on the basis of formula (6) to improve, as shown in the following formula:
Among them, σ is the standard deviation, which determines the smoothness of the Gaussian filter. The larger the value of σ, the thinner the edge of the image, σ2 = 1.5 × σ1. Mark the improved difference image as , and then analyze each pixel point. If the value is greater than or equal to ε, the value of this point is set as 1; otherwise, the function is added and subtracted, as shown in the following equation:
In addition to the σ variable that comes with the DoG function, XDoG newly introduces three parameters γ, Φ, and ε to control the parameters of the output image. γ controls the amount of detail in the entire output image. The larger the value of γ, the more details in the output image. In the image detection of tunnel lining cracks, the value of γ is usually [0.95, 1). When the γ value is small, the thinner cracks in the image are not displayed completely. In this paper, after many experiments, it is determined that when the γ value is 0.98, the output image will achieve the best effect. The specific image details are shown in Figure 4. The Φ value controls the brightness of the image. The smaller the Φ value, the darker the image. After specific experiments, it can be seen that the change in the Φ value has almost no effect on the image details. The value range can be as long as (200, 2000]. ε controls the shade of gray and has little effect on the recognition result whose range is [−0.01, 0.01].
(a)
(b)
(c)
(d)
2.3. Morphological Noise Removal
The fracture image of the tunnel lining processed by the XDoG can be clearly seen by the naked eye. However, due to the texture of the concrete surface and other influencing factors, there are a lot of noises around the cracks and the target boundary is relatively rough, and also there are noise spots of different sizes in the background area. The above noise will affect the subsequent crack recognition accuracy and crack labeling, so the method of morphology is used to remove the noise. First, perform the color inversion of the binary image on the image after the XDoG. Turn the background area to black and the cracks and noise area to white. Use Matlab’s built-in regionprops function to remove smaller noise points. Using the opening and closing operations in morphology for the image that has been denoised can better remove the remaining noise in the image. The specific implementation method is shown in the following formulas:
In the above formula, B is the collection of all structural elements, which can be decomposed into two parts and . A is the original image collection, X represents the image to be processed, represents the element in the set X, and XC represents the complement of the set X. If is an empty set, it is defined as a corrosion operation, which is calculated according to formula (11) to reduce or corrode the entire set. When is an empty set, it is defined as expansion. According to formula (12), the set is expanded or enlarged. The operating principle of the open operation is expressed as equation (13), that is, B performs corrosion calculation on A and then performs expansion calculation on the corrosion result. The closed operation principle is shown in equation (14), and B performs expansion calculation on A and then performs corrosion calculation on the expansion result. The result of multiple opening and closing operations is shown in Figure 5.
It can be seen from Figure 5 that after multiple morphological opening and closing operations, the noise near the crack and in the background has been basically removed. Compared with Figure 5, the overall width of the cracks is increased and clearer, but there are certain fractures, as shown in Figure 6(a). Using the breakpoint connection algorithm proposed by Wang et al. [20] to connect the cracks and fractures, the experimental results are shown in Figure 6(b). It can be seen from the figure that the fractured parts are connected, making the overall image more realistic.
(a)
(b)
3. Crack Direction Classification
Tunnel lining cracks can be roughly divided into two types: one is linear cracks, which can be divided into horizontal cracks, vertical cracks, and oblique cracks, and the other is network cracks. Different fractures have different effects on the tunnel lining force [21, 22], so it must be distinguished. The existing tunnel crack classification methods are mostly based on the projection method [23] and density [24] which are too idealistic. Peng et al. [25] give two premises. First, a picture can only have one crack target, and second, the smallest bounding rectangle of the crack has been found. However, although it has given a fracture treatment method, the applicable object of this method is pavement cracks. Pavement cracks are different from tunnel fractures to a certain extent, and its judgment method cannot be completely applied to the tunnel. This paper improves the existing algorithm and proposes two new indicators for judging the type of cracks.
3.1. Fracture Contour Extraction
Perform sliding window detection on images that have been binarized. The size of the detection window is 25 pixels × 25 pixels (each pixel corresponds to an actual distance of 1 mm), the image size is 1125 pixels × 750 pixels, and the red window itself occupies one pixel. If there are white pixels, the sliding detection window turns red for identification. If there are no white pixels, proceed to the next window inspection, and the resulting image after inspection is shown in Figure 7. The newly proposed method can not only save time and save the process of manual image annotation but also avoid the influence of color image noise on the accuracy of sliding detection. It can be clearly shown from the figure that each small red box contains white pixels, and the entire crack is completely surrounded by the red box, which is conducive to subsequent processing.
The region contour is extracted from the binary image that has been divided into subblocks. The extracted crack contour is formed by closing horizontal or vertical line segments connected end to end. The direction of the line segment is defined as horizontal and vertical. The subblock picture after region contour extraction is shown in Figure 8.
3.2. Identification of Crack Types
The extracted tunnel crack profile features can be used to classify cracks, such as the most commonly used area length, width, overall aspect ratio, and density features. Some scholars have proposed that it is not accurate enough to judge the types of cracks by relying on physical characteristics. Therefore, projection features and moment invariant features are introduced to distinguish the types of cracks. The projection feature is effective in identifying horizontal and vertical cracks, but it is not suitable for the discrimination of oblique cracks and network cracks. The calculation of the features of the moment invariants is complicated and easily affected by noise and image resolution power. Therefore, it is critical to develop a crack classification method with high recognition accuracy.
This paper introduces the block index to distinguish linear cracks and network cracks and introduces the histogram-like method to distinguish different linear fractures. In the crack image that has been contour extracted, there is a part of the blank area between the extracted contour and the image boundary, as shown in Figure 9, where the diagonal line is the contour image of the crack that has been removed. The overall image part outside the contour area is defined as Ω. The difference between linear cracks and network cracks lies in the number of subsets in the Ω set. Each unconnected irregular white block is automatically labeled starting from 1, and the set Ω contains elements Ω1, Ω2,… and so on. The maximum number of element subscripts represents the number of disconnected subblocks in the set. Through analyzing the block size of 100 pieces of linear crack pictures and network crack pictures, the results are shown in Table 1. It can be seen from the table that the Ω value of 4 is used as the dividing point between linear fractures and network fractures. If the Ω value is less than 4, it is considered as a linear fracture, and if the Ω value is greater than or equal to 4, it is classified as a network. crack.
For the classification of linear cracks, the method of histogram-like is introduced to identify the cracks. This method can identify oblique cracks well and can handle the situation of different types of cracks and two to three cracks of the same type in an image. The complete contours in the linear crack image that have been distinguished from the network cracks are labeled and represented by (the maximum i value in this paper is 3). Take the circumscribed rectangle of each contour, and establish a rectangular coordinate system at the lower left corner of each circumscribed rectangle. At this time, each small red box is in the rectangular coordinate system.
Suppose the length of the center of the small block from the x-axis (calculate the length from the vertical pixels from the x-axis to the center of the small block, and the obtained length is the pixel) is ai, then ai must be a multiple of 12, store all ai in set A. The distance between the center of the small block and the y-axis is bi. Since the small block is square, bi is also a multiple of 12, and all bi is stored in set B. Perform data analysis on the obtained sets A and B, and the results are shown in Table 2.
From the comparison of the data in the table, it can be seen that for horizontal fractures, the number of different elements in the dataset A is obviously less than that of the vertical and oblique cracks, while the values in the B set are almost different. For vertical cracks, it is the opposite of horizontal cracks, the value of B set is less, while the value of A set is almost different. The elements in the sets A and B of oblique cracks are different. Using the difference in the number of different elements in the set, a new parameter Δ is introduced to distinguish the fracture direction, where Δ = N/M. The number of different elements in set A is denoted as M, the number of different elements in set B is denoted as N, and Δ is the ratio of N to M. The difference of Δ value in different crack directions is shown in Table 3.
4. Algorithm Performance Analysis
This paper uses the GUI function of the visual development platform MatlabR2016a to process the collected images of the tunnel lining cracks offline. The computer parameters are Intel i7 10875H processor, 8G memory, 1 TB solid-state drive, 8 cores, and 16 threads.
The XDoG crack detection algorithm proposed in this paper is compared with another edge detection algorithm, namely, the gradient-improved Canny algorithm [26]. It is also compared with the traditional crack recognition algorithm based on SVM [27] and the crack detection algorithm based on maximum entropy segmentation [28] (with both black as the background and white as the detected tunnel cracks). Randomly select four images from the 600 image dataset to show the comparison effect with other algorithms. The original image is shown in Figure 10(a), which are horizontal cracks, vertical cracks, network cracks, and crack image with noise. The image based on the gradient-improved Canny algorithm is shown in Figure 10(b), the SVM-based crack recognition algorithm is shown in Figure 10(c), and the crack detection algorithm based on maximum entropy segmentation is shown in Figure 10(d). The crack detection algorithm proposed in this paper is shown in Figure 10(e).
(a)
(b)
(c)
(d)
(e)
Obvious contrast of the four algorithms can be seen from the above figure. Although the gradient-improved Canny algorithm can generally identify cracks, this method is most susceptible to the influence of concrete texture. The identification of cracks in certain tunnel linings is prone to generate a lot of complex noise, which has a great impact on subsequent treatment. The SVM-based crack recognition algorithm is better than the improved Canny algorithm and the maximum entropy segmentation result and can retain the structural characteristics of the crack as a whole. The linear noise filtering in the image is better, but the blocking noise filtering effect is relatively general. There are many discontinuities in the identified cracks, which are somewhat different from the actual situation. The crack detection algorithm based on maximum entropy segmentation has poor recognition results, but it has a good recognition effect for cracks with higher contrast. The background noise filtering effect is better than the Canny algorithm, but there are also cases where the block noise filtering is incomplete. Although some cracks can be identified, it will lead to missing cracks and corresponding noises in the background. The algorithm in this paper has been greatly optimized in terms of crack recognition accuracy and the number of background noise points and can basically reproduce the shape of tunnel cracks, but some impurities near the cracks may affect the recognition results.
In order to verify the accuracy of the algorithm, the algorithm proposed in this paper was used to compare the ordinary concrete crack image, tunnel crack image, and asphalt pavement crack image. The results are shown in Table 4. The theoretical algorithm in this paper shows that the crack recognition rate of ordinary concrete is as high as 99.0%. For tunnel cracks, due to the local leakage of water and other noise interference on the lining surface, the recognition rate is lower than that of ordinary concrete surface cracks, but the recognition rate can still reach 94.5%. For asphalt pavement cracks, because the asphalt has more textures, the recognition rate is relatively low, only 87.0%. After the inner surface of the tunnel is automatically photographed by a CCD camera, the complexity of the algorithm is a problem that must be considered due to the large number of pictures. Algorithm complexity can be expressed by algorithm time in a certain sense. The faster the processing speed of each picture, the higher the engineering application value. Table 5 shows the calculation time comparison of various algorithms in this paper. In terms of image processing time, the average processing time of the gradient-improved Canny algorithm is 1.2 s, but the processing effect is relatively worst. The algorithm based on SVM takes the longest time of 8.7 s, but the accuracy is relatively high while the maximum entropy segmentation requires the shortest time, and the average time is only 1.1 s. The mean time of this algorithm for processing an image is 5.2 s. Compared with the deep learning model based on the segNet architecture studied by this group, the recognition rate is only 4.2% lower, but the time spent is greatly reduced. The accuracy and performance of the algorithm can continue to be improved in the future.
The results of crack type recognition are shown in Table 6. Select 100 horizontal cracks, vertical cracks, oblique cracks, network cracks and mixed cracks each from the crack library to classify the crack directions. Because the characteristics of horizontal cracks and vertical cracks are relatively obvious, the recognition rate is the highest, reaching 98% and 96%, respectively. Because of the oblique cracks with fuzzy characteristics, the recognition accuracy is also low, only reaching 88%. Due to the obvious characteristics of network cracks, the accuracy rate can be as high as 94% under the identification of this paper. Mixed cracks are two or more cracks in a picture including horizontal, vertical, and oblique cracks, and the overall recognition accuracy can reach 90%.
5. Conclusion
(1)This paper proposes a new type of tunnel crack image processing algorithm. First, use improved homomorphic filtering to remove haze and shadow processing on the image and use grayscale image adaptive median filter preprocessing. The XDoG is used for edge extraction, and after morphological opening and closing operations are performed on the binary image, the fracture and breakpoint connection processing is performed.(2)Use a sliding window of 25 pixels × 25 pixels to detect the pixels of the processed image, and classify the types of tunnel cracks through the proposed new method of block index and histogram-like. The results prove that the newly proposed crack classification method can better distinguish between network cracks and linear cracks, and the recognition results are accurate. The accuracy of distinguishing the types of linear cracks is also high, and it can detect two or more cracks in one image.(3)Use MatlabR2016a to set the image processing GUI to process the image. The processing results show that the newly proposed crack detection algorithm can better identify most tunnel lining cracks, and the recognition rate can be as high as 94.5%, which is higher than the existing traditional crack identification methods. But there is still a certain gap with deep learning methods. The calculation time of each image is moderate, which has good engineering applicability.(4)In future research, for the accuracy of crack recognition, image semantic recognition can be carried out through methods related to deep learning to increase the recognition accuracy. For the classification of crack types, new physical parameters such as aggregation degree can be introduced on the basis of this article to conduct a more detailed analysis of image attributes, thereby improving the accuracy of crack classification.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Disclosure
Chunquan Dai and Kun Jiang are co-first authors.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
Chunquan Dai and Kun Jiang contributed equally to the article.
Acknowledgments
The authors are grateful for the financial support from the Humanities and Social Sciences Fund of the Ministry of Education (no. 20YJAZH022).