# Three-Dimensional Reconstruction of Tunnel Face Based on Multiple Images

Authors: Wenge Qiu; Liao Jian; Yunjian Cheng; Hengbin Bai
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The current geological sketch in tunnel engineering is mainly based on sketches of workers. However, geological sketch drawn by workers always offers fundamental data purely due to its drawing mode. A novel drawing method for geological sketch has been introduced using multiview photos in this process. The images of tunnel faces are taken from multiple angles, and every two pictures have overlaps. By measuring the distance between the camera and the tunnel face using a laser range finder, the photographic scale of each photo can be confirmed. SpeededUp Robust Features (SURF) is a good practice for detecting feature points, and the sparse point cloud is reconstructed from multiview photos by structure from motion (SFM). However, the sparse point cloud is not suitable for analysis for structural planes due to its sparsity. Therefore, patch-based multiview stereo (PMVS) is used to reconstruct dense point cloud from the sparse point cloud. After 3D reconstruction, the details of the tunnel face are recorded. The proposed technique was applied to multiview photos acquired in the Xiaosanxia railway tunnel and Fengjie tunnel in Chongqing, China. In order to record the geological conditions of the tunnel face quickly and accurately, Chengdu Tianyou Tunnelkey has developed a set of software and hardware integration system called CameraPad. Besides, CameraPad was used to collect the multiview photos of the tunnel face in the No. 1 Xinan railway tunnel in Jilin, China. By comparing with traditional and existing methods, the proposed method offers a more reductive model for geological conditions of the tunnel face.


---

## Body

## 1. Introduction

Tunnel engineering is concealed engineering, which is often affected by a variety of adverse geological conditions. As the geological conditions of the area that the tunnel traverses are becoming more and more complex, tunnel engineering has also become extremely challenging. And geological prospecting is becoming more and more critical. However, relying merely on the geological prospecting of the previous period is an utterly inadequate method in dealing with a difficult situation for the complex geological conditions. Therefore, the geological prospecting during the tunnel construction period will be more meaningful for tunnel engineering. The stratum lithology, geological structure, and structural plane occurrence in the tunnel construction period need to be recorded to provide raw data, which are used for the interpretation and analysis of geotectonic movement and the prediction of long-term stability of rock mass. At the same time, it also provides an essential geological foundation to optimize tunnel excavation and strengthening methods.Traditionally, the geological sketch of the tunnel face is often used to record the geological conditions during the tunnel construction period. The geological conditions of the tunnel face are often very complicated, and it is difficult to completely record the geological conditions in the traditional geological sketch.It is difficult to observe the tunnel face directly when the TBM and shield method are used in recent years, so many scholars have proposed many new techniques for collecting the geological conditions of the tunnel face. Yamamoto et al. developed the TBM Excavation Control System to acquire the accurate prediction of the geological conditions ahead of and surrounding the tunnel face in tunnel excavation by Tunnel Boring Machine (TBM). When using the TBM method, it is difficult to observe the tunnel face during the tunnelling process. Besides, it is difficult to obtain the geological conditions during the construction period, so such a method is proposed [1]. Li et al. used a digital optical borehole camera to directly observe the gaps in the rock mass through predrilled holes [2]. This method has been applied to the prediction of advanced geological conditions, which can make an excellent assessment of the geological conditions in front of the tunnel face.The equipment for the data collection of 3D reconstruction has also become lighter and more portable. Izadi et al. used a Kinect camera that can be held and moved freely to quickly perform 3D reconstruction of indoor scenes [3]. Track the position and posture of the sensor through the depth data using Kinect camera to reconstruct the 3D scene in time. This method has shown how to utilize low-cost equipment for rapid 3D reconstruction. At the same time, terrestrial laser scanning (TLS) technology is also used for the construction of the three-dimensional model of the tunnel face. Wang et al. proposed a new method for tunnel digital geological mapping and recording using 3D terrestrial laser scanning (TLS) [4]. The point cloud obtained from TLS was used to generate the most suitable structural plane on the tunnel face, and then the direction of the structural plane was calculated. The result has shown the digital geological sketch based on the unfolded image of the tunnel face, including the traces and directions of the fissures, classification of rock types, and water leakage. Cacciari and Futai proposed a method of using a terrestrial laser scanner (TLS) to generate discrete fracture networks (DFNs) of geological mapping [5]. The proposed method has effectively established the fractured rock masses of tunnel faces. Cacciari and Futai attempted to evaluate the current methods of geometrical discontinuity characterization by terrestrial laser scanning (TLS) in tunnel engineering [6]. By comparison, the mean trace lengths obtained from rectangular sampling windows are more suitable for estimating different geological surface orientations. Cheng et al. developed a general technique to automatically identify the different types of components using TLS [7]. Li et al. proposed deformation monitoring using TLS in Underground Station Cavity Clusters [8].Traditional geological prospecting methods also have incorporated many new technologies to predict geological conditions in front of the tunnel. McCrae and Cook proposed a geological mapping for tunnels at relatively shallow depths (<250 m) in strong discontinuous rocks, focusing on a means of locating possible instabilities and selecting support requirements on the basis of rock mass characterization. Therefore, characterizing the behaviour of discontinuities has become the major requirements in the geological mapping process [9]. Zarei et al. used geological mapping to evaluate high local groundwater inflow to a rock tunnel [10]. Due to the given hydrogeological assumptions and the simplification of heterogeneous media, the use of analytical and numerical tools to predict possible groundwater inflows from these characteristics often failed. Therefore, through detailed geological features, reliable groundwater inflow can be predicted. Xiong et al. proposed a regional scale 3D geological model for assessing tunnel engineering, based on regional field surveys of geological boundary and attitude [11]. Afterwards, according to the construction process, geological information was obtained from the tunnel face. This method dynamically adjusts the three-dimensional geological model according to the information obtained at different stages, providing a reference for tunnel engineering. Soldo et al. proposed the importance of site investigation of geological conditions for tunnel engineering. All the collected data must be interpreted inside a robust conceptual framework and primarily came from field surveys in geology, enhanced by new methods and technologies for investigating and modelling [12].As a result, the geological conditions of the tunnel face play an important role in the tunnel construction period. The tunnel design needs to be compatible with the actual geological conditions in front of the tunnel, achieving a safer and more economical construction.In order to preserve the geological data and present them in a three-dimensional visualization form, a method was proposed for the establishment of a three-dimensional geological model. A novel image preprocessing scheme can accurately and quickly take out the tunnel face area as the raw data for 3D modelling. Since it takes a long time to perform high-precision modelling of the tunnel face, the proposed image preprocessing scheme can effectively reduce the amount of calculation in the reconstruction process. The process of 3D reconstruction of the tunnel face is faster than before. At the same time, the photographing effect is poor due to the jitter in the handheld shooting. Besides, it is obviously problematic that the camera site is too close or too far to the tunnel face. The authors also designed a device named CameraPad to automatically take clear pictures of the tunnel face. The device can rotate the camera in three axes automatically to collect tunnel face photos after being placed in front of the tunnel as required. Due to the built-in gyroscope, the relative position of each photo can be obtained after rotating in three axes, and the camera exterior parameters can also be obtained. Compared with the traditional method of estimating exterior parameters through algorithms, it is more accurate.To sum up, this article can offer research objectives as follows: ① a novel method of rapid 3D reconstruction of the tunnel faces automatically; ② a new device to collect data for fast 3D reconstruction of the tunnel faces; ③ application of the proposed method.


## 2. Background and Related Work

In today’s tunnel engineering, geological conditions play an increasingly important role, affecting the days and cost of tunnel construction. Therefore, the record of geological conditions during construction has also been paid more and more attention. In order to record the geological conditions of the tunnel face quickly and accurately, Chengdu Tianyou Tunnelkey has developed a set of software and hardware integration system called CameraPad. As shown in Figure1, the main components of the CameraPad are divided into four parts: ① camera module; ② distance measuring module; ③ gyroscope module; and ④ pan/tilt/zoom (PTZ) module. Among them, the camera module is used to obtain multiview photos, and the distance measuring module composed of the laser range finder is used to measure the distance from the shooting position to the tunnel face to determine the photographing scale. Besides, the PTZ module is mainly used to measure and adjust attitude parameters, and the gyroscope is used to record the orientation of the lens while taking photos.Figure 1
Module composition and usage scenarios of CameraPad.
(a)(b)First of all, put CameraPad in a suitable position in front of the tunnel face. Besides, roughly align the lens to the centerline of the tunnel face. The device is turned on, and the laser range finder will be used as an aid to select a suitable position. Place two geological boards on the tunnel face in advance to convert the reconstructed model from the camera coordinate system to the geodetic coordinate system [7]. Take down the photosource equipped with the CameraPad and place it in a suitable position so that the tunnel face can be photographed clearly. After a one-button start, CameraPad will use the PTZ to adjust the posture to collect the image of the tunnel face.Then, CameraPad can be wirelessly connected and the data can be exported in the hotspot mode, and then the 3D reconstruction is performed.Finally, mark structural planes, joints, fissures, and lithology in the 3D reconstructed model.The algorithm flow of the software associated with CameraPad is shown in Figure2.Figure 2
The algorithm flow of the software associated with CameraPad.


## 3. Image Preprocessing

As shown in Figure3, multiple pictures from the tunnel face in Xiaosanxia railway tunnel in Chongqing were regarded as input in this process. Since multiple photos are used to reconstruct the 3D geological model of the tunnel face, the surrounding primary lining will be captured while taking pictures. In order to speed up the reconstruction process, the area of the captured tunnel face is used as the reconstructed image. Besides, the area of primary lining is taken as a non-reconstructed image area. Therefore, the 3D geological reconstruction of the tunnel face can speed up and reduce the amount of calculation. We extract the area of the tunnel face and use it as the basic data for 3D reconstruction, as shown in Figure 4(e).Figure 3
The images of the tunnel face obtained from the Xiaosanxia railway tunnel in Chongqing, China.
(a)(b)(c)(d)Figure 4
Extraction of tunnel face area from the image. (a) The original image of tunnel face. (b) The image after Gaussian filtering. (c) The image after region growth processing. (d) The area of tunnel face extracted. (e) The result after image preprocessing.
(a)(b)(c)(d)(e)First, the Gaussian filtering operation is performed on the picture of the tunnel face taken by the camera shown in Figure4(a). Bergholm proposed a Gaussian filter for noise reduction and edge detection [13]. As shown in Figure 5, in the 3 × 3 pixels’ area of the photograph, the pixel at the center point is the weighted average of all nine-pixel values, and the value of each pixel is calculated as follows:(1)Gx,y=12πσ2e−x2+y2/2σ2,where σ=6 is selected after multiple sets of data verification and x, y represent the position relative to the center of the area shown in Figure 5, with the origin at the center. The X-axis is horizontal to the right, and the Y-axis is vertically upward.Figure 5
Gaussian filter.Gaussian filter is a linear smoothing filter, which is widely used in the noise reduction process of image processing. In essence, Gaussian filtering is the process of weighted averaging the entire image. The value of each pixel is obtained by the weighted mean of itself and other pixel values in the neighborhood. The specific operation of Gaussian filtering is to scan each pixel in the image with a template, and the weighted average intensity of the pixels in the neighborhood determined by the template was used to replace the intensity of the center pixel of the template.As shown in Figure4(b), after the Gaussian filtering operation, the noise in the image is removed to facilitate subsequent image processing. The region growing algorithm is used to remove the areas that meet the preset threshold [14, 15], as shown in Figure 4(c). In this way, the area of primary lining in the picture can be roughly determined, as shown in Figure 4(d). Thus, the area of the tunnel face in the picture is further defined, as shown in Figure 4(e). We combined Gaussian filter with region growing algorithm to extract the tunnel face area in the images.


## 4. The Extraction and Matching of Feature Points

### 4.1. The Extraction of Feature Points

Feature points are used to describe the corresponding relationship that can be used to match two photos. Feature points should have apparent characteristics that are clearly different from nonfeature points and can be repeatedly detected in different photos. Feature points are essential for seeking the homonymy points in other images, and a lot of research studies have been done in the field of computer vision. In view of the inconvenience of taking photos at the construction site, the feature point extraction with scale invariance and rotation invariance is adopted.SURF (SpeededUp Robust Features) is an accelerated version based on the Scale Invariant Feature Transformation (SIFT) algorithm [16]. The SURF operator is much faster than the SIFT operator in detecting feature points of the image [17, 18]. In the case of multiple pictures in the same scene, the extracted feature points have better stability.Feature point detection based on Hessian matrix is the core of SURF:(2)detHx,σ=Lxxx,σLxyx,σLxyx,σLyyx,σ=Lxxx,σLyyx,σ−Lxyx,σ2.Among them,Lxxx,σ represents the convolution of the Gaussian second-order derivative ∂2gσ/∂x2 of the image at the position of the pixel x=x,y with a scale of σ; Lyyx,σ represents the convolution of the Gaussian second-order derivative ∂2gσ/∂y2 of the image at the position of the pixel x=x,y with the scale of σ; and Lxyx,σ represents the convolution of the Gaussian second-order derivative ∂2gσ/∂x∂y of the image at the position of the pixel x=x,y with the scale of σ.The SURF algorithm makes full use of the fast integration algorithm of the rectangular area to accelerate the convolution calculation of the Box Matrix and the source graph. Besides, the analysis of the Haar wavelet used in the feature vector extraction is also speeded up. Therefore, the SURF algorithm can be used in computer vision object recognition and 3D reconstruction.The initialization of Box Matrix is as follows:(1)
Given an image of sizeM × N, set a red rectangular box with a size of m × n, as shown in Figure 6.(2)
Create a new one-dimensional array of sizeM to store intermediate variables in the calculation process, represented by red pixels. The values of the red pixels are the sum of each column of pixels.(3)
Assuming that the coordinate (0, 0) on the top left corner of the red rectangular box is regarded as the center of the movement, slide the rectangular box from (0, 0) to (M-m, 0) pixel by pixel in the right. Whenever the red rectangular box is moved to a new position, calculate the sum of the pixels in the red rectangular box. After moving to the next row, push the red rectangular box to (0, 1) and slide the red rectangular box from (0, 1) to (M-m, 1) pixel by pixel in the right. Repeat this until the red rectangular box slides from (0, N-n) to the right pixel by pixel to (M-m, N-n).(4)
The calculation of the sum of pixels in the red rectangular box, first sum each column of pixels in the orange rectangular box, place the result in the red pixel and then sum the red pixels in the corresponding green rectangular box. Therefore, the result is the sum of the red rectangular box.Figure 6
The principle of Box Filter.The SURF (SpeededUp Robust Features) further simplifies the second derivative of the Gaussian function, using 9 × 9 Box Matrix instead of Gaussian second-order derivatives withσ=1.2, as shown in Figure 6. Happrox is denoted by Dxx, Dyy, and Dxy. The relative weights in Hessian matrix need to be further balanced with Lxy1.2FDxx9F/Lxx1.2FDxy9F=0.912≈0.9 for computational efficiency, where Dxx9F is the Frobenius norm of Dxx [17].After simplification, take the value of the Hessian matrix as follows:(3)detHapprox=DxxDyy−0.9Dxy2.The SURF feature point detection is also called Fast Hessian Detector. The Hessian matrix is the core of the entire SURF algorithm and is a second-order matrix composed of the second partial derivative of a real-valued function whose independent variable is a vector. The results of feature point matching are shown in Figure7.Figure 7
Feature point matching results. The blue lines in the figure are the matches that meet the preset threshold condition, and the red lines are the matches that do not meet the preset threshold condition.


### 4.2. Distribution and Matching of Feature Points in Main Direction

As shown in Figure8, the SURF algorithm counts the Harr wavelet features of the circular neighborhood of feature points. A 60° fan-shaped area is used in the circular neighborhood of the feature points, and the sums of the horizontal and vertical Harr wavelet features of all points in the 60° fan-shaped area are counted. And the fan-shaped area is rotated 60° around the feature point in turn, so as to traversal the six 60° fan-shaped regions of the circular neighborhood. Therefore, the Harr wavelet values in the six fan-shaped regions are summed to obtain six vector directions, and the vector direction with the largest modulus is the main direction of the feature point.Figure 8
Determination of the main direction of feature points.
(a)(b)(c)


### 4.3. The Generation of Feature Point Descriptors

After determining the main direction of the feature points, rotate the coordinate axis in the main direction. Select a 20σ × 20σ area surrounding a feature point and divide it into 4 × 4 small rectangular area, as shown in Figure 9. Then, calculate the 25 Haar wavelet values of size σ × σ in each small rectangular area. According to the principle of separation of positive and negative, calculate the sum of wavelet values in the x-direction and y-direction in the area: Σdx, Σdx, Σdx, Σdy. With eight gradient directions in each small rectangular area, a 4 × 4 × 8 = 128-dimensional real number vector is obtained [17]. After the feature points are extracted, the Euler distance of the two feature vectors can be calculated in the vector space of the feature points. This paper adopts Nearest Neighbor Distance Ratio (NNDR) for matching:(4)NNDR=DClosest PointDNext nearest point,where DClosest is the Euler distance between the feature point description vector and the closest feature point description vector and DNext nearest point is the Euler distance between the feature point description vector and the next nearest feature point description vector.Figure 9
The generation of feature point descriptors.
(a)(b)When NNDR is less than the preset threshold, the match is accepted.


## 4.1. The Extraction of Feature Points

Feature points are used to describe the corresponding relationship that can be used to match two photos. Feature points should have apparent characteristics that are clearly different from nonfeature points and can be repeatedly detected in different photos. Feature points are essential for seeking the homonymy points in other images, and a lot of research studies have been done in the field of computer vision. In view of the inconvenience of taking photos at the construction site, the feature point extraction with scale invariance and rotation invariance is adopted.SURF (SpeededUp Robust Features) is an accelerated version based on the Scale Invariant Feature Transformation (SIFT) algorithm [16]. The SURF operator is much faster than the SIFT operator in detecting feature points of the image [17, 18]. In the case of multiple pictures in the same scene, the extracted feature points have better stability.Feature point detection based on Hessian matrix is the core of SURF:(2)detHx,σ=Lxxx,σLxyx,σLxyx,σLyyx,σ=Lxxx,σLyyx,σ−Lxyx,σ2.Among them,Lxxx,σ represents the convolution of the Gaussian second-order derivative ∂2gσ/∂x2 of the image at the position of the pixel x=x,y with a scale of σ; Lyyx,σ represents the convolution of the Gaussian second-order derivative ∂2gσ/∂y2 of the image at the position of the pixel x=x,y with the scale of σ; and Lxyx,σ represents the convolution of the Gaussian second-order derivative ∂2gσ/∂x∂y of the image at the position of the pixel x=x,y with the scale of σ.The SURF algorithm makes full use of the fast integration algorithm of the rectangular area to accelerate the convolution calculation of the Box Matrix and the source graph. Besides, the analysis of the Haar wavelet used in the feature vector extraction is also speeded up. Therefore, the SURF algorithm can be used in computer vision object recognition and 3D reconstruction.The initialization of Box Matrix is as follows:(1)
Given an image of sizeM × N, set a red rectangular box with a size of m × n, as shown in Figure 6.(2)
Create a new one-dimensional array of sizeM to store intermediate variables in the calculation process, represented by red pixels. The values of the red pixels are the sum of each column of pixels.(3)
Assuming that the coordinate (0, 0) on the top left corner of the red rectangular box is regarded as the center of the movement, slide the rectangular box from (0, 0) to (M-m, 0) pixel by pixel in the right. Whenever the red rectangular box is moved to a new position, calculate the sum of the pixels in the red rectangular box. After moving to the next row, push the red rectangular box to (0, 1) and slide the red rectangular box from (0, 1) to (M-m, 1) pixel by pixel in the right. Repeat this until the red rectangular box slides from (0, N-n) to the right pixel by pixel to (M-m, N-n).(4)
The calculation of the sum of pixels in the red rectangular box, first sum each column of pixels in the orange rectangular box, place the result in the red pixel and then sum the red pixels in the corresponding green rectangular box. Therefore, the result is the sum of the red rectangular box.Figure 6
The principle of Box Filter.The SURF (SpeededUp Robust Features) further simplifies the second derivative of the Gaussian function, using 9 × 9 Box Matrix instead of Gaussian second-order derivatives withσ=1.2, as shown in Figure 6. Happrox is denoted by Dxx, Dyy, and Dxy. The relative weights in Hessian matrix need to be further balanced with Lxy1.2FDxx9F/Lxx1.2FDxy9F=0.912≈0.9 for computational efficiency, where Dxx9F is the Frobenius norm of Dxx [17].After simplification, take the value of the Hessian matrix as follows:(3)detHapprox=DxxDyy−0.9Dxy2.The SURF feature point detection is also called Fast Hessian Detector. The Hessian matrix is the core of the entire SURF algorithm and is a second-order matrix composed of the second partial derivative of a real-valued function whose independent variable is a vector. The results of feature point matching are shown in Figure7.Figure 7
Feature point matching results. The blue lines in the figure are the matches that meet the preset threshold condition, and the red lines are the matches that do not meet the preset threshold condition.


## 4.2. Distribution and Matching of Feature Points in Main Direction

As shown in Figure8, the SURF algorithm counts the Harr wavelet features of the circular neighborhood of feature points. A 60° fan-shaped area is used in the circular neighborhood of the feature points, and the sums of the horizontal and vertical Harr wavelet features of all points in the 60° fan-shaped area are counted. And the fan-shaped area is rotated 60° around the feature point in turn, so as to traversal the six 60° fan-shaped regions of the circular neighborhood. Therefore, the Harr wavelet values in the six fan-shaped regions are summed to obtain six vector directions, and the vector direction with the largest modulus is the main direction of the feature point.Figure 8
Determination of the main direction of feature points.
(a)(b)(c)


## 4.3. The Generation of Feature Point Descriptors

After determining the main direction of the feature points, rotate the coordinate axis in the main direction. Select a 20σ × 20σ area surrounding a feature point and divide it into 4 × 4 small rectangular area, as shown in Figure 9. Then, calculate the 25 Haar wavelet values of size σ × σ in each small rectangular area. According to the principle of separation of positive and negative, calculate the sum of wavelet values in the x-direction and y-direction in the area: Σdx, Σdx, Σdx, Σdy. With eight gradient directions in each small rectangular area, a 4 × 4 × 8 = 128-dimensional real number vector is obtained [17]. After the feature points are extracted, the Euler distance of the two feature vectors can be calculated in the vector space of the feature points. This paper adopts Nearest Neighbor Distance Ratio (NNDR) for matching:(4)NNDR=DClosest PointDNext nearest point,where DClosest is the Euler distance between the feature point description vector and the closest feature point description vector and DNext nearest point is the Euler distance between the feature point description vector and the next nearest feature point description vector.Figure 9
The generation of feature point descriptors.
(a)(b)When NNDR is less than the preset threshold, the match is accepted.


## 5. Three-Dimensional Reconstruction Based on Improved Structure-from-Motion (SFM) and Patch-Based Multiview Stereo (PMVS) Software

### 5.1. Improved Structure from Motion (SFM)

The structure-from-motion (SFM) algorithm is an algorithm for sparse point cloud based on collected disordered pictures, that is, to recover the three-dimensional structure of an object from motion (images taken at different times) [19]. In this process, it is necessary to estimate the exterior parameters involved with position parameters and attitude parameters when taking pictures and reconstruct the sparse point cloud with the intrinsic parameters of the camera.For the initial matching results of the feature vectors obtained above, a random sample consensus (RANSAC) is used to remove the invalid feature vectors [20]. Through accurate matching results, the homonymy points in multiple photos are also determined, as shown in Figure 10. Then, according to the principle of pinhole imaging [21], a series of projection equations are established. The difference values between the estimated projecting position and the actual projecting position of points are calculated through the projection equation set. These difference values are accumulated to construct the objective function [22]:(5)fKi,Di,Ri,ti,Xj=minKi,Di,Ri,ti,Xj∑i,jxji−x^ji2,where i represents the input order of photos; Ki represents the intrinsic parameter matrix of the camera that took the i photo; Di represents the radial distortion parameter; Ri|Ti represents the exterior parameter matrix of the camera taking the i photo (where Ri represents the rotation matrix and ti represents the translation matrix); and Xj=xjyjzj represents the spatial point.Figure 10
Schematic diagram of the three-dimensional reconstruction of SFM.The spatial information of the feature points is obtained by solving the above formula, and the intrinsic parameters are obtained by using calibration boards [23]. The exterior parameters, intrinsic parameters, and spatial information of the traditional SFM algorithm are all estimated by the algorithm, while the exterior parameters and intrinsic parameters of the proposed method are obtained externally. Only the spatial information is estimated by the algorithm, which will obviously greatly improve the accuracy and speed of 3D visualization. The extracted point cloud is shown in Figure 11 [24].Figure 11
The point cloud extracted after SFM.


### 5.2. Patch-Based Multiview Stereo (PMVS)

PMVS is a patch-based multiview stereo algorithm proposed by Furukawa in 2007 [25]. Due to its high accuracy and wide application scenarios, it is popularized on a large scale. The point cloud in Figure 11 is so sparse that the output result of the previous step is used as input for the dense reconstruction of the sparse point cloud.Each feature point that was successfully matched in the previous step is defined as a seed patch, which contains the following attributes:(6)cp⟵Triangulation fromfandf′,np⟵cpOIi⟶cpOIi⟶,Rp⟵Ii.Among them,cp represents the center of the patch p; np represents the normal vector of the patch p; Rp represents the reference picture of the patch p; f and f′ represent the detected feature points that match successfully; and OIi represents the camera optical center corresponding to the picture Ii.The photometric discrepancy functiongp is defined as(7)gp=1Vp/Rp∑I∈Vp/Rphp,I,Rp,where Vp denotes a set of images in which p is visible and hp,I,Rp is defined to be a pairwise photometric discrepancy function between picture I and the reference picture Rp of patch p.gp is used to select the best one from the set of matching feature points of patch p as the final patch. Next, the patch p is expanded, and the neighbor patch p′ is initialized with the information of patch p. When all the patch sets are expanded, the patches with weak grayscale consistency and geometric consistency are removed.Therefore, the feature pointf is obtained by improved structure from motion (SFM). Assuming that the corresponding patch is p, cp and np are obtained as variables through the above formula. The photometric discrepancy energy is evaluated many times in the energy minimization procedure performed by a conjugate gradient method to find the optimal patch as the final patch p. Given a patch p, we initialize p′ as a neighbor of p. Through a more complex optimization process to optimize cp′ and np′, the new patches are expanded. After the expansion, a filtering operation will be performed to remove the patches of weaker visibility consistency.Compared with Multiphoto Geometrically Constrained Matching (MPGC) and Maximization of the Posterior Marginal (MPM), the PMVS algorithm is more straightforward and faster in terms of the amount of calculation. Besides, the input parameters of PMVS are required less. The intrinsic parameter matrix and exterior parameter matrix of the image output by SFM can be used for dense reconstruction of the point cloud based on multiple images. Therefore, the method proposed in this paper uses PMVS for the dense reconstruction of the point cloud. The dense reconstruction of the point cloud based on PMVS is shown in Figure12. The final 3D reconstruction effect based on PMVS is shown in Figure 13.Figure 12
The point cloud extracted after dense reconstruction.Figure 13
The final result after 3D reconstruction.


## 5.1. Improved Structure from Motion (SFM)

The structure-from-motion (SFM) algorithm is an algorithm for sparse point cloud based on collected disordered pictures, that is, to recover the three-dimensional structure of an object from motion (images taken at different times) [19]. In this process, it is necessary to estimate the exterior parameters involved with position parameters and attitude parameters when taking pictures and reconstruct the sparse point cloud with the intrinsic parameters of the camera.For the initial matching results of the feature vectors obtained above, a random sample consensus (RANSAC) is used to remove the invalid feature vectors [20]. Through accurate matching results, the homonymy points in multiple photos are also determined, as shown in Figure 10. Then, according to the principle of pinhole imaging [21], a series of projection equations are established. The difference values between the estimated projecting position and the actual projecting position of points are calculated through the projection equation set. These difference values are accumulated to construct the objective function [22]:(5)fKi,Di,Ri,ti,Xj=minKi,Di,Ri,ti,Xj∑i,jxji−x^ji2,where i represents the input order of photos; Ki represents the intrinsic parameter matrix of the camera that took the i photo; Di represents the radial distortion parameter; Ri|Ti represents the exterior parameter matrix of the camera taking the i photo (where Ri represents the rotation matrix and ti represents the translation matrix); and Xj=xjyjzj represents the spatial point.Figure 10
Schematic diagram of the three-dimensional reconstruction of SFM.The spatial information of the feature points is obtained by solving the above formula, and the intrinsic parameters are obtained by using calibration boards [23]. The exterior parameters, intrinsic parameters, and spatial information of the traditional SFM algorithm are all estimated by the algorithm, while the exterior parameters and intrinsic parameters of the proposed method are obtained externally. Only the spatial information is estimated by the algorithm, which will obviously greatly improve the accuracy and speed of 3D visualization. The extracted point cloud is shown in Figure 11 [24].Figure 11
The point cloud extracted after SFM.


## 5.2. Patch-Based Multiview Stereo (PMVS)

PMVS is a patch-based multiview stereo algorithm proposed by Furukawa in 2007 [25]. Due to its high accuracy and wide application scenarios, it is popularized on a large scale. The point cloud in Figure 11 is so sparse that the output result of the previous step is used as input for the dense reconstruction of the sparse point cloud.Each feature point that was successfully matched in the previous step is defined as a seed patch, which contains the following attributes:(6)cp⟵Triangulation fromfandf′,np⟵cpOIi⟶cpOIi⟶,Rp⟵Ii.Among them,cp represents the center of the patch p; np represents the normal vector of the patch p; Rp represents the reference picture of the patch p; f and f′ represent the detected feature points that match successfully; and OIi represents the camera optical center corresponding to the picture Ii.The photometric discrepancy functiongp is defined as(7)gp=1Vp/Rp∑I∈Vp/Rphp,I,Rp,where Vp denotes a set of images in which p is visible and hp,I,Rp is defined to be a pairwise photometric discrepancy function between picture I and the reference picture Rp of patch p.gp is used to select the best one from the set of matching feature points of patch p as the final patch. Next, the patch p is expanded, and the neighbor patch p′ is initialized with the information of patch p. When all the patch sets are expanded, the patches with weak grayscale consistency and geometric consistency are removed.Therefore, the feature pointf is obtained by improved structure from motion (SFM). Assuming that the corresponding patch is p, cp and np are obtained as variables through the above formula. The photometric discrepancy energy is evaluated many times in the energy minimization procedure performed by a conjugate gradient method to find the optimal patch as the final patch p. Given a patch p, we initialize p′ as a neighbor of p. Through a more complex optimization process to optimize cp′ and np′, the new patches are expanded. After the expansion, a filtering operation will be performed to remove the patches of weaker visibility consistency.Compared with Multiphoto Geometrically Constrained Matching (MPGC) and Maximization of the Posterior Marginal (MPM), the PMVS algorithm is more straightforward and faster in terms of the amount of calculation. Besides, the input parameters of PMVS are required less. The intrinsic parameter matrix and exterior parameter matrix of the image output by SFM can be used for dense reconstruction of the point cloud based on multiple images. Therefore, the method proposed in this paper uses PMVS for the dense reconstruction of the point cloud. The dense reconstruction of the point cloud based on PMVS is shown in Figure12. The final 3D reconstruction effect based on PMVS is shown in Figure 13.Figure 12
The point cloud extracted after dense reconstruction.Figure 13
The final result after 3D reconstruction.


## 6. Application

The Xiaosanxia railway tunnel is located in Chongqing, China. In this tunnel, the CameraPad system developed by Chengdu Tianyou Tunnelkey was used to collect data. As shown in Figure13, the images of the tunnel face were mainly gathered to perform the three-dimensional geological reconstruction along the tunnel excavation line. As shown in Figure 13, the high-definition three-dimensional reconstruction model and the built-in gyroscope of CameraPad are very helpful for the identification of structural planes and the drawing of joints and fissures in the future. Not only that, but the 3D geological model in the geodetic coordinate system can also provide a basis for measuring the dip direction.The definition of the 3D geological model obtained by the proposed method can be adjusted by replacing the camera module of the CameraPad. The high-definition camera module can provide better results for identifying joint and fissures in the future. It can be seen that the high-definition camera module plays an essential role in recognition of the tunnel face with the crushed surrounding rock. Compared to 3D reconstruction from video flowing, CameraPad is also more suitable for use in tunnel engineering. The video flowing at 25 frames per second is used for 3D reconstruction [26], which is very convenient for some scenes that require online 3D reconstruction. Since online demand usually severely limits the available resolution of vidicon, this method also has the disadvantage of relatively low resolution.CameraPad developed by Chengdu Tianyou Tunnelkey was used to collect data in Fengjie railway tunnel in Chongqing, China. As shown in Figure14, the images of the tunnel face were mainly gathered to perform the 3D geological reconstruction along the tunnel excavation line.Figure 14
(a) The results of 3D reconstruction in Fengjie railway tunnel in Chongqing. (b) 3D geological model arranged according to the design of the rail line.
(a)(b)The proposed method placed these 3D geological models according to the tunnel design data and the mileage of the tunnel face, as shown in Figure14. In this way, the geological conditions during tunnel construction can be better represented and used as geological data storage. Compared with the traditional geological sketch, the proposed method not only is more intuitive but also can provide geological data for construction design more effectively. Through the geological conditions of the tunnel face, we can also make a prediction for the following geological conditions.CameraPad developed by Chengdu Tianyou Tunnelkey was used to collect data in No. 1 Xinan railway tunnel in Jilin, China. As shown in Figure15, the images of the tunnel face were mainly gathered to perform the 3D geological reconstruction along the tunnel excavation line.Figure 15
The results of 3D reconstruction in No. 1 Xinan railway tunnel in Jilin.The proposed method placed these 3D geological models according to the tunnel design data and the mileage of the tunnel face, as shown in Figure15.


## 7. Conclusions and Future Work

In this paper, the SURF algorithm was used to extract and match the feature points from the multiple photos of the tunnel face. Next, SFM is used to extract the sparse point cloud, and then the dense reconstruction of the point cloud is performed based on the patch-based multiview stereo algorithm. So far, the picture is pasted to the corresponding patch, and the three-dimensional reconstruction is completed. In order to speed up the process of 3D reconstruction, the proposed method preprocessed the multiple photos and extracted part of the tunnel face for reconstruction. The 3D geological model reconstructed in this way has the characteristics of full color and a high degree of precision. Compared with the traditional hand-drawn geological sketch, the method proposed in this paper is more efficient, and the geological conditions of the tunnel face can be entirely recorded by assigning values to the patches. Compared with the three-dimensional geological model obtained by the TLS method, the method proposed in this paper has the characteristics of full color, which makes it easier to identify the seepage and the stratum lithology. The proposed method can also be combined with BIM technology to build a BIM model during the construction period so that the record has a better manifestation.In this paper, the proposed method has the advantages of 3D visualization and high precision compared with traditional geological data preservation methods. A novel image preprocessing scheme can accurately and quickly take out the tunnel face area as the raw data for 3D modelling. Due to the built-in gyroscope, the relative position of each photo can be obtained after rotating in three axes, and the camera exterior parameters can also be obtained. Compared with the traditional method of estimating exterior parameters through algorithms, it is more accurate.There are still many problems to be solved. For example, the automatic identification of the attitude of rock formation and the structural plane of the three-dimensional geological model of the tunnel face are well down. The automatic identification of joints, fissures, and water leakage is also in future work. Of course, it is also hoped that in the future, automatic identification of rock formations and lithology will be achieved to complete the automatic 3D geological sketch.



---

*Source: 8837309-2021-04-19.xml*