Rails play an important role in the maintenance and repair of railway lines. As a result of regularly inspecting the rail profile, as well as evaluating the state parameters of these rails, such as vertical wear and side wear [1–3], it is possible to gain a better understanding of not only the state of the rails but also how to grind the rails, as this is a crucial part of railway operations and maintenance [4,5]. The measurement of rail profiles using line-structured light vision is based on the principle of triangulation and features high speed, high precision, and noncontact. As a mainstream method of dynamic detection of rail profiles globally, it can detect the parameters of in-service rails, such as vertical wear and side wear [6–8]. In most cases, a line-structured light sensor is placed on each side of a rail in order to obtain the profile data of the left and right half-sections of the rail. The full-section profile of the rail is then produced by splicing these two half-sections together [9–11]. In this process, if the laser beams of the line-structured light sensors are not on the same plane on either side of the rail, the measurement profile will be distorted to a certain degree, resulting in errors in rail profiling. It is therefore necessary to accurately assess the laser plane attitude on both sides of the rail to ensure that the laser planes on both sides are coplanar in order to obtain high-precision full-section profile data. When the lasers on both sides are installed, they are incident on the calibration plate, forming a line of intersection between the two lasers. The calibration plate is marked with a scale line. In order to determine whether the laser planes are coplanar as required, it is necessary to visually observe the degree of coincidence between the intersection lines of the laser planes on both sides and the scale line of the calibration plate. Obviously, this method is limited to a qualitative analysis and cannot provide a quantitative or accurate assessment of the plane attitude of the lasers on both sides.

**Citation:** Wang, L.; Wang, H.; Han, Q.; Fang, Y.; Wang, S.; Wang, N.; Li, G.; Ren, S. A Laser Plane Attitude Evaluation Method for Rail Profile Measurement Sensors. *Sensors* **2023**, *23*, 4586. https://doi.org/10.3390/ s23104586

Academic Editor: Simone Borri

Received: 10 March 2023 Revised: 17 April 2023 Accepted: 21 April 2023 Published: 9 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Currently, little attention is being paid to the attitude of the laser planes on both sides when measuring rail profiles. In most existing studies, the focus is on adjusting the lasers in order to make them coplanar. As an example, Zhan et al. proposed a mechanism and method for the adjustment of laser planes. A mechanical adjustment mechanism could be used to translate and rotate the line-structured laser [12]. Through a high-precision manual translation stage, Chen et al. were able to scan the profiles of large workpieces [13]. In addition, some scholars have been able to align the laser planes on both sides by examining calibration techniques. In one study, Wu et al. were able to calibrate profile measurement components globally through the use of special mechanical parts, such as rotating arms [14]. In their study, Zhang et al. documented the positions of the laser profiler and markers by photographing them and converting the adjustment of coplanarity into the position adjustment of the laser profiler [15]. Ju et al. calibrated the laser planes by using the contour line calibration method [16]. Wang et al. also identified a method for correcting the error in rail profile measurement caused by non-coplanar lasers. Using projection transformation, a laser non-coplanarity correction model was proposed, as well as a reference coordinate system based on the longitudinal direction of rails. Half-section profile data were projected onto an auxiliary plane perpendicular to the longitudinal direction of the rail, and the projection profile was used to correct the measurement results [17]. Currently, few reports have been published on the use of the laser plane attitude evaluation method for both sides of the rail during the full-section measurement of the rail profile. In summary, existing methods cannot quantitatively and accurately evaluate the attitude of laser planes on both sides.

The author previously studied the distribution characteristics of laser non-coplanar error in rail profile measurement sensors and proposed a correction method for laser noncoplanar error and a calibration method for rail longitudinal parameters. However, during the installation process of rail profile measurement devices, the laser coplanar adjustment operation is mainly guided by observing the degree of collinearity of the laser line with the naked eye, which has a certain degree of blindness. In order to avoid this blindness and guide laser coplanar adjustment operations, this paper proposes a quantitative evaluation method for laser plane posture. When the measuring device is installed, three planar targets of different heights are used to obtain the attitude information of both laser planes in real time. Such information is used to establish a laser coplanarity evaluation criterion, which will be used to guide the alignment of the two line-structured lasers. As opposed to naked-eye evaluations with low precision, poor real-time performance, and subjectivity, the proposed method uses computer vision evaluations for laser coplanarity adjustment, with high precision, excellent real-time performance, and visualization, which can reduce calibration errors. In light of this, it is of great importance to improve the accuracy of the measurement of rail profiles in full sections.

#### **2. Basic Principle**

Figure 1 presents a schematic diagram of the full-section rail profile measurement based on line-structured light vision. Both sides of the rail are equipped with a linestructured light sensor, and the lasers of both sensors are positioned in the same plane in order to obtain the left and right half-sections of the rail. In order to obtain the full-section profile of the rail, the half-section profiles on both sides are spliced together according to the calibration parameters [18]. Using the scanning motion, the rail profile for the entire railway line can be measured in the full section. The degree of coplanarity of the laser planes on both sides of the rail is an important factor in determining the accuracy of a full-section rail profile measurement system based on line-structured light vision. Ideally, the laser planes on both sides should be coplanar in order to ensure the accuracy of rail profile measurements. As shown in Figure 2, the origin is the center of the top surface of the rail, the laser plane is the XOZ plane, and the horizontal direction is the X axis of the world coordinate system, O-XYZ. The non-coplanarity of the laser planes on both sides of the rail can therefore be expressed as the rotation of the laser plane around the X and

Z axes, respectively. In Figure 3, the lasers on both sides are not coplanar, in which the left laser plane is perpendicular to the longitudinal direction of the rail, but the right laser plane rotates around the Z axis, at an angle. In this case, the right laser plane does not remain perpendicular to the longitudinal direction of the rail, nor is it coplanar to the left laser plane. While the left camera still captures the profile data of the cross-section perpendicular to the rail longitudinal direction, the right camera captures the profile data not perpendicular to the rail longitudinal direction, referred to as the oblique-section profile data. In comparison with the cross-section profile of the rail, the oblique-section profile is stretched in a certain direction, and this stretching direction is directly related to the angle between the laser planes and the longitudinal direction of the rail. As a result of this stretching, the rail profile is distorted, causing deviations in the positioning of the feature points, leading to an increase in the measurement error of the rail wear. More generally, the measured profiles on both sides of the rail are not the cross-section profile of the rail when the laser planes are not coplanar or perpendicular to the longitudinal direction of the rail. In addition, the measurements will be distorted, resulting in greater errors in determining rail profile.

**Figure 1.** Schematic diagram of the full-section rail profile measurement based on line-structured light vision.

At present, the laser coplanarity of the rail profile measurement system is determined by observing the degree of overlap of two laser lines with the human eye. As shown in Figure 4, the surface of the aluminum alloy ruler is engraved with a long scale line. When installing the laser scan sensors on both sides, the laser planes on both sides are projected onto the same scale line of the aluminum alloy ruler, as shown in Figure 5. The intersection line of the laser planes on both sides and the scale line of the aluminum alloy ruler are observed with the naked eye to determine whether the laser planes on both sides meet the coplanar installation requirements based on the degree of overlap between the intersection line and the scale line. Obviously, this method evaluates the attitude of a two-dimensional plane through one-dimensional lines, which can only be qualitatively evaluated and cannot be quantitatively evaluated and has a certain degree of blindness. Due to the low level of visualization, it cannot effectively guide the coplanar installation operation of two laser scan sensors. The main purpose of this article is to propose a quantitative evaluation method

for laser plane attitude, which visualizes and quantifies the adjustment process of the laser plane, thereby guiding the coplanar installation operation of two laser scan sensors.

**Figure 2.** Schematic diagram of coplanar lasers on both sides: (**a**) coplanar lasers and (**b**) profile measurement results.

**Figure 3.** Schematic diagram of non-coplanar lasers on both sides: (**a**) non-coplanar lasers and (**b**) profile measurement results.

**Figure 4.** Aluminum alloy ruler.

**Figure 5.** Schematic diagram of the traditional laser plane evaluation method for rail profile measurement system.

#### **3. Materials and Methods**

As can be seen from the analysis above, the existing method is based on visual observation of whether the lasers on both sides are coplanar; however, there are many uncertainties involved in the evaluation process, which cannot guarantee the accuracy of calibration. Therefore, this study proposes a method for evaluating laser coplanarity based on fitting planes, which can serve as a guide for the adjustment of laser coplanarity on both sides. Figure 6 illustrates the visualized laser plane adjustment device for the full-section rail profile measurement system. It contains an additional target compared to the existing full-section rail profile measurement device. The special target consists of a convex calibration block and three planar calibration plates. There are three upper surfaces on the convex calibration block, and the upper surface in the center is higher than the upper surfaces on both sides. As shown in Figure 7, three calibration plates are placed on the three upper surfaces, which are numbered 1, 2, and 3, from left to right. Consequently, the corresponding target coordinate systems, *tcs*1, *tcs*2, and *tcs*3, are established, with the center of each calibration plate as the origin and the target plane as the XOY plane.

**Figure 6.** The visualized laser plane adjustment device for the full-section rail profile measurement system.

**Figure 7.** The coordinate system, which is composed of the convex calibration block and three planar targets.

An overview of the laser plane attitude evaluation for the full-section rail profile measurement system is presented in Figure 8. This system consists of a system calibration module, an image acquisition module, and a coplanarity evaluation module. The specific realization process for each module is described in more detail below.

The system calibration module is responsible for obtaining the internal and external parameters of the cameras. The calibration method previously described [19] is employed in this study to simultaneously collect images of the planar target calibration plates in different attitudes through the left and right cameras in order to obtain the cameras' internal parameters.

**Figure 8.** Flowchart of laser coplanarity evaluation of the full-section rail profile measurement system.

The convex calibration block is located in the common field of view of the left and right cameras, ensuring that the left laser plane intersects the target planes of *tcs*1 and *tcs*2, and the right laser plane intersects the target planes of *tcs*2 and *tcs*3. After the lasers on both sides are turned off, the left and right cameras are used to capture images of the convex calibration block. Due to occlusion, as shown in Figure 9, the left camera can capture the entire calibration plates No. 1 and No. 2, but only a portion of calibration plate No. 3, and the right camera can capture the entire calibration plates No. 2 and No. 3, but only a portion of calibration plate No. 1. The coordinate systems of the left and right cameras are expressed as *ccs*1 and *ccs*2, respectively, and *Rccs*1 *tcs*1 , *t ccs*1 *tcs*1 represent the rotation matrix and translation vector of the coordinate systems *ccs*1 and *tcs*1. Since the internal parameters of the cameras are known, based on the camera calibration method as previously described, the rotation matrix, *Rccs*1 *tcs*1 , *Rccs*1 *tcs*2 , and translation vector, *t ccs*1 *tcs*1 , *t ccs*1 *tcs*2 , of the coordinate systems from the left camera *ccs*1 to No. 1 and No. 2 targets, as well as the rotation matrix, *Rccs*2 *tcs*2 , *Rccs*2 *tcs*3 , and translation vector, *t ccs*2 *tcs*2 , *t ccs*2 *tcs*3 , of the coordinate systems from the right camera *ccs*2 to No. 2 and No. 3 targets are calculated based on the images of the convex calibration block.

**Figure 9.** The images of the convex calibration block captured by the left and right cameras.

The image acquisition module is used to acquire the real-time light stripe images of the convex calibration block during the process of adjusting the laser planes. During the calibration process, the positions of the convex calibration block, the planar calibration plates, and the cameras remain unchanged, while the lasers on both sides are turned on and adjusted as required. With a suitable exposure time, the cameras on both sides are used to collect the light strip images of the convex calibration block in real time. The light strip image sequence of the convex calibration block is denoted as follows:

$$I = \{I\_{\text{ij}}/i = 1, 2, j = 1, 2, 3 \dots n\} \tag{1}$$

where *i* = 1 is the light stripe image of the convex calibration block collected by the left camera, *i* = 2 is the light stripe image of the convex calibration block collected by the right camera, and *n* is the number of their respective light stripe images. The intersection line of the laser planes and the calibration plates form a light stripe, as shown in Figure 10. The light stripe image of the convex calibration block collected by the left camera shows the intersection lines between the left laser plane and the planar calibration plates No. 1 and No. 2, which are designated as *l*1 and *l*2, respectively. The light stripe image of the convex calibration block collected by the right camera shows the intersection lines between the right laser plane and the planar calibration plates No. 2 and No. 3, which are denoted as *r*2 and *r*3, respectively.

Based on the system calibration parameters and the convex calibration block light stripe image, the coplanarity evaluation module calculates the parameters of the left and right laser planes in order to determine whether the lasers on both sides are coplanar. The coplanarity evaluation module operates in six steps, according to the data processing flow.

**Figure 10.** The left and right light stripe images of the convex calibration block.

Step 1: The pixel coordinates of the light stripe centers are obtained by extracting the centers of the left and right light stripe images of the convex calibration block. In Figure 10, the light strip images of the convex calibration block collected by the cameras on both sides are shown, and they indicate the intersection lines *l1, l2*, *r*2, and *r*3. As follows, the light stripe centers are extracted using traditional algorithms (such as maximum value methods, grayscale center-of-gravity methods, Steger methods, template-matching methods, etc.), and the pixel coordinates of the left and right light stripe centers are obtained.

$$P\_{\bar{i}} = \left\{ (u\_{i\prime}v\_{\bar{i}})^T | 0 \le u \le width-1, 0 \le v \le he \text{light}-1, \bar{i} = 1, 2 \right\} \tag{2}$$

where the width is the image width, and the height is the image height; *P*1 is any point on the center of the left light stripe image; and *P*2 is any point on the center of the right light stripe image.

Step 2: The pixel coordinates of the left and right light stripe image centers of the convex calibration block are converted into the corresponding target coordinate system. Different calibration plates correspond to different external parameters. The coordinate transformation process is described below, using the light strip image of the convex calibration block captured by the left camera in Figure 10 as an example. First, as shown in Figure 10, the light stripe *l*1 of the calibration plate No. 1 and the light stripe *l*2 of the calibration plate No. 2 are located. Then, for the pixel coordinate of the light stripe center of the calibration plate No. 1, (*ui*, *vi*) *T*, the pixel coordinate of the center of the light stripe *l*1, (*ui*, *vi*) *T*, is transformed into the coordinate system of the No. 1 target according to the internal and the external parameters (*Rccs*1 *tcs*1 and *t ccs*1 *tcs*1 ) of the camera. Similarly, the pixel coordinate of the center of the light stripe *l*2, (*ui*, *vi*), is transformed into the coordinate system of the No. 2 target according to the internal and external parameters (*Rccs*1 *tcs*2 and *t ccs*1 *tcs*2 ) of the camera. A similar coordinate transformation is applied to the light stripe image of the convex calibration block captured by the right camera. As a result, we have the coordinate of the intersection line *l*1 in the coordinate system of target No. 1, *tcs*1 (*Pl*1 *tcs*1 = *xl*1 *tcs*1, *yl*1 *tcs*1, *zl*1 *tcs*1 *T* ); the coordinate of the intersection line *l*2 in the coordinate system of target No. 2, *tcs*2 (*Pl*2 *tcs*2 = *xl*2 *tcs*2, *yl*2 *tcs*2, *zl*2 *tcs*2 *T* ); the coordinate of the intersection line *r*2 in the coordinate system of target No. 2, *tcs*2 (*Pr*2 *tcs*2 = *xr*2 *tcs*2, *yr*2 *tcs*2, *zr*2 *tcs*2 *T* ); and the coordinate of the intersection line *r*3 in the coordinate system of target No. 3, *tcs*3 (*Pr*3 *tcs*3 = *xr*3 *tcs*3, *yr*3 *tcs*3, *zr*3 *tcs*3 *T* ).

Step 3: The coordinate system of target No. 2 is regarded as the world coordinate system, *wcs*, and the coordinates of the light stripe centers in the respective target coordinate systems are transformed into the world coordinate system. The transformation relationship between the coordinate systems *tcs*2, *tcs*1, and *tcs*3 is calculated. The rotation matrix, *Rtcs*2 *tcs*1, and translation vector, *t ccs*2 *tcs*1 , from the target coordinate system *tcs*2 to the target coordinate system *tcs*1 are calculated by Equation (3):

$$\begin{cases} \begin{array}{c} R\_{tcs1}^{tcs2} = \left(R\_{tcs2}^{ccs1}\right)^{-1} \cdot R\_{tcs1}^{ccs1} \\\\ t\_{tcs1}^{tcs2} = \left(R\_{tcs2}^{ccs1}\right)^{-1} \cdot \left(t\_{tcs1}^{ccs1} - t\_{tcs2}^{ccs1}\right) \end{array} \tag{3}$$

The rotation matrix and translation vector from the target coordinate system *tcs*2 to the target coordinate system *tcs*3 are calculated by Equation (4):

$$\begin{cases} \begin{array}{l} R\_{tcs3}^{tcs2} = \left(R\_{tcs2}^{ccs2}\right)^{-1} \cdot R\_{tcs3}^{ccs2} \\\\ t\_{tcs3}^{tcs2} = \left(R\_{tcs2}^{ccs2}\right)^{-1} \cdot \left(t\_{tcs3}^{ccs2} - t\_{tcs2}^{ccs2}\right) \end{array} \tag{4}$$

The world coordinates of the light stripe center on the intersection lines *l*2 and *r*2 are the same as the coordinates in the target coordinate system, that is, *Pl*2 *wcs* = *Pl*2 *tcs*2 and *Pr*2 *wcs* = *Pr*2 *tcs*2. The coordinates of the intersection lines *l*1 and *r*3 in the world coordinate system, *wcs*, namely *Pl*1 *wcs* = *xl*1 *wcs*, *yl*1 *wcs*, *zl*1 *wcs T* and *Pr*3 *wcs* = *xr*3 *wcs*, *yr*3 *wcs*, *zr*3 *wcsT* , are obtained by Equations (5) and (6), respectively. Similarly, we have the coordinate of the intersection line *l*2 in the world coordinate system (*Pl*2 *wcs* = *xl*2 *wcs*, *yl*2 *wcs*, *zl*2 *wcs T* ), the coordinate of the intersection line *r*2 in the world coordinate system (*Pr*2 *wcs* = *xr*2 *wcs*, *yr*2 *wcs*, *zr*2 *wcsT* ), and the coordinate of the intersection line *r*3 in the world coordinate system (*Pr*3 *wcs* = *xr*3 *wcs*, *yr*3 *wcs*, *zr*3 *wcsT* ).

$$
\begin{bmatrix} P\_{\text{ucc}}^{\text{II}} \\ 1 \end{bmatrix} = \begin{bmatrix} \mathbf{x}\_{\text{ucc}}^{\text{II}} \\ y\_{\text{ucc}}^{\text{II}} \\ z\_{\text{ucc}}^{\text{II}} \\ 1 \end{bmatrix} = \begin{bmatrix} R\_{\text{tcs1}}^{\text{tcs2}} & t\_{\text{tcs1}}^{\text{tcs2}} \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} P\_{\text{tcs1}}^{\text{I}} \\ 1 \end{bmatrix} = \begin{bmatrix} R\_{\text{tcs1}}^{\text{tcs2}} & t\_{\text{tcs1}}^{\text{tcs2}} \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x\_{\text{tcs1}}^{\text{II}} \\ y\_{\text{tcs1}}^{\text{II}} \\ z\_{\text{tcs1}}^{\text{II}} \\ 1 \end{bmatrix} \tag{5}
$$

$$
\begin{bmatrix} P\_{\text{wcs}}^{r3} \\ 1 \end{bmatrix} = \begin{bmatrix} x\_{\text{wcs}}^{r3} \\ y\_{\text{wcs}}^{r3} \\ z\_{\text{wcs}}^{r3} \\ 1 \end{bmatrix} = \begin{bmatrix} R\_{\text{tcs}3}^{\text{tcs2}} & t\_{\text{tcs3}}^{\text{tcs2}} \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} P\_{\text{tcs3}}^{\text{tcs2}} \\ 1 \end{bmatrix} = \begin{bmatrix} R\_{\text{tcs3}}^{\text{tcs2}} & t\_{\text{tcs3}}^{\text{tcs2}} \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x\_{\text{tcs3}}^{r3} \\ y\_{\text{tcs3}}^{r3} \\ z\_{\text{tcs3}}^{r3} \\ 1 \end{bmatrix} \tag{6}
$$

Step 4: In the world coordinate system, the two intersection lines (*l*1 and *l*2) of the left laser plane and the two intersection lines (*r*2 and *r*3) of the right laser plane are sequentially fitted to the plane, and the parameters of the left and right laser planes are obtained. The fitting process is described below, taking the intersection lines *l*1 and *l*2 as an example.

In *Pi wcs* = *xi wcs*, *yi wcs*, *zi wcsT* , *i* = 1, 2, 3 ... *k* is any point on the intersection lines *l*1 and *l*2 of the left laser plane and the calibration plates No. 1 and No. 2, where *k* = *m* + *n*, and *m* and *n* are the number of points on the intersection lines *l*1 and *l*2, respectively. The following matrix is thus constructed:

$$M = \begin{bmatrix} \mathbf{x}\_{\text{wcs}}^1 - \overline{\mathbf{x}}\_{\text{'}} \mathbf{x}\_{\text{wcs}}^1 - \overline{\mathbf{y}}\_{\text{'}} \mathbf{x}\_{\text{wcs}}^1 - \overline{\mathbf{z}} \\ \vdots \\ \mathbf{x}\_{\text{wcs}}^k - \overline{\mathbf{x}}\_{\text{'}} \mathbf{y}\_{\text{wcs}}^k - \overline{\mathbf{y}}\_{\text{'}} \mathbf{x}\_{\text{wcs}}^1 - \overline{\mathbf{z}} \end{bmatrix} \tag{7}$$

where *x* = 1 *k k* ∑ *i*=1 *xi wcs*, *y* = 1 *k k* ∑ *i*=1 *yi wcs*, *z* = 1 *k k* ∑ *i*=1 *zi wcs*. The point *y*(*x*, *y*, *z*) *T* is designated as the center of gravity of the plane. If *S* = *MT*·*M*, where *S* has three eigenvalues, then the eigenvector corresponding to the smallest eigenvalue is the normal of the fitting plane, *nl* . The left laser plane is constructed with (*x*, *y*, *z*) *T* as a point on the plane and the vector *nl* as the normal. Similarly, the intersection lines *r*2 and *r*3 of the right laser plane are fitted

to the plane to obtain the normal of the right laser plane, *nr* , and the right laser plane is constructed.

Step 5: The angle between the normals of the fitting planes on the left and right sides (*α*) and the distance between the planes (*d*) are calculated. The angle, *α*, can be calculated by Equation (8):

$$\alpha = \arccos\left(\frac{\boldsymbol{n}^l \cdot \boldsymbol{n}^r}{|\boldsymbol{n}^l| |\boldsymbol{n}^r|}\right) \tag{8}$$

In order to calculate the distance between the two planes, the corresponding plane coordinate systems *pcs*1 and *pcs*2 are established with the center of gravity of the fitting planes on the left and right sides as the origin and the normal direction as the Z axis. The rotation of the coordinate system around the Z axis will not affect the direction of the plane normal in the plane coordinate system. For simplicity, the rotation around the Z axis is set to 0 here. The rotation matrix and translation vector from the plane coordinate system *pcs*1 to the world coordinate system (*wcs*) are denoted as *Rpcs*1 *wcs* and *t pcs*1 *wcs* , and the rotation matrix and translation vector from the world coordinate system (*wcs*) to the plane coordinate system *pcs*2 are expressed as *Rwcs pcs*2 and *t wcs pcs*2. As a result, the corresponding homogeneous transformation matrices are constructed, as shown below.

$$H\_{\rm wcs}^{\rm pcs1} = \begin{bmatrix} R\_{\rm wcs}^{\rm pcs1} & t\_{\rm wcs}^{\rm pcs1} \\ 0 & 0 & 0 \end{bmatrix} \tag{9}$$

$$H\_{\rm pcs2}^{\rm wcs} = \begin{bmatrix} R\_{\rm pcs2}^{\rm wcs} & t\_{\rm pcs2}^{\rm wcs} \\ 0 & 0 & 0 \end{bmatrix} \tag{10}$$

Then the homogeneous transformation matrix of the coordinate system *pcs*1 and *pcs*2 can be expressed as follows:

$$H\_{pcs2}^{pcs1} = H\_{\text{wcs}}^{pcs1} \cdot H\_{pcs2}^{\text{uvcs}} \tag{11}$$

Therefore, for any point in the coordinate system *pcs*2, *Ppcs*2 = *xpcs*2, *ypcs*2, *zpcs*2 *T* , it can be transformed into the coordinate system *pcs*1 by Equation (12), where *Ppcs*1 = *xpcs*1, *ypcs*1, *zpcs*1 *T* is the corresponding coordinate of the point in the coordinate system *pcs*1.

$$
\begin{bmatrix} p\_{\rm pcs1} \\ 1 \end{bmatrix} = \begin{bmatrix} \chi\_{\rm pcs1} \\ \mathcal{Y}\_{\rm pcs1} \\ z\_{\rm pcs1} \\ 1 \end{bmatrix} = H\_{\rm pcs2}^{\rm pcs1} \cdot P\_{\rm pcs2} = H\_{\rm pcs2}^{\rm pcs1} \begin{bmatrix} \chi\_{\rm pcs2} \\ \mathcal{Y}\_{\rm pcs2} \\ z\_{\rm pcs2} \\ 1 \end{bmatrix} \tag{12}
$$

*N* points on the right laser plane are arbitrarily taken, as shown in Equation (13), and transformed into the coordinate system *pcs*1 through Equation (9). Then, we have Equation (14).

$$P\_{\rm pcs2}^{i} = \left(\mathbf{x}\_{\rm pcs2}^{i}, y\_{\rm pcs2}^{i}, z\_{\rm pcs2}^{i}\right)^{T}, i = 1, 2, 3 \dots N \tag{13}$$

$$P\_{p\text{cs1}}^{i} = \left(x\_{p\text{cs1}}^{i}, y\_{p\text{cs1}}^{i}, z\_{p\text{cs1}}^{i}\right)^{T}, i = 1, 2, 3 \dots N \tag{14}$$

Since the fitting plane on the left coincides with the *XOY* plane of the coordinate system *pcs*1, *zi pcs*1 , *i* = 1, 2, 3 ... *N* is the distance from these *N* points to plane 1, and the distance (*d*) from plane 2 to plane 1 can be expressed as follows:

$$d = \frac{1}{N} \sum\_{i=1}^{i=N} \left| z\_{p \text{cs1}}^i \right| \tag{15}$$

Step 6: Both the distance between the planes and the angle between the normals determine whether the lasers on both sides are coplanar. If the angle between the two planes is 0, it means that the two planes are parallel or coincident. If the distance from any point on one of the planes to the other plane is 0, the two planes coincide. To eliminate the influence of errors, two parameters are used to determine whether the lasers on both sides are coplanar, namely the angle (*α*) and the distance (*d*). As long as the angle and distance satisfy Equation (16), the coincidence degree of the left and right laser planes is high, and the laser planes on both sides of the rail are appropriate.

$$d \le T\_d \text{ \& } \alpha \le T\_d \tag{16}$$

where *Td*, *Tα* are the thresholds of the distance (*d*) and the angle (*α*), which can be determined in accordance with the accuracy requirements.

When the laser planes are adjusted, the left and right cameras collect real-time light stripe images of the convex calibration block and obtain the sequence of light stripes, as shown in Equation (2). Then, the coplanarity evaluation module processes the light strip image sequence of the convex calibration block in real time, calculates the laser plane parameters on both sides of the rail, draws the laser planes on both sides in real time in the display window, and displays the angle and distance between the normals. By utilizing Equation (16), the module determines if the laser planes on both sides of the rail are appropriate, allowing the full-section rail profile measurement system to visualize the alignment of both laser planes and effectively determine if the lasers are coplanar on both sides. The effect is shown in Figure 11.

**Figure 11.** Schematic diagram of the visualized adjustment of the laser planes on both sides of the rail.

#### **4. Results and Discussion**

#### *4.1. Experimental Design*

To confirm the advantages of the proposed laser coplanarity evaluation method, such as a strong real-time performance, high accuracy, and visualization, a laser coplanarity evaluation experiment was designed, and the experimental device shown in Figure 12 was constructed. Among them, the camera is the Ranger 3 high-speed camera produced by the German SICK company (Waldkirch, Germany), with a resolution of 2560 × 832 pixels, and the line laser is produced by the Canadian Osela company, with a wavelength of 660 nm and 450 nm, respectively. The laptop used for test is ThinkPad T480, Intel®, Core™, i7-8550U<EMAIL_ADDRESS>1.99 GHz, 32.0 GB memory, 64-bit operating system.

**Figure 12.** Experimental device for laser coplanarity evaluation: (**a**) calibration of the internal parameters of the cameras and (**b**) simplified convex calibration block and calibration plates.

First, the internal parameters of the left and right cameras are calibrated based on the planar target, and then the simplified convex calibration block is placed in the common field of view of the two cameras, ensuring that the left laser plane intersects the calibration plates No. 1 and No. 2, and the right laser plane intersects the calibration plates No. 2 and No. 3. Then, the lasers on both sides are turned off, and the images of the calibration plates are collected through the cameras on both sides. The result is shown in Figure 9. After that, the lasers on both sides are turned on, and the light strip images of the three calibration plates are collected in real time through the cameras on both sides. The result is shown in Figure 10. The parameters of the laser planes on both sides are obtained in real time from the cross-sectional laser images of the three calibration plates, and the two laser planes are drawn in real time in the program window. The lasers on both sides are adjusted according to the two laser planes displayed in real time in the program window until they satisfy the coplanarity standard. The following three groups of experiments were carried out:

The first group of experiments is a real-time verification experiment. On the basis of the abovementioned experimental device, the laser planes on both sides are continuously adjusted, and two cameras collect 500 target images each in real time. At the same time, the coplanar degree of the two laser planes is calculated by the proposed method, and the real-time performance of the proposed method is evaluated by the time of program brushing the new window.

The second group of experiments is the accuracy verification experiment. In this experiment, the accuracy of the proposed method is evaluated by the laser plane fitting error, laser plane distance measurement error, and laser plane angle measurement error. First, in the process of laser plane adjustment on both sides, 500 consecutive target images are collected by two cameras, corresponding to 500 × 2 = 1000 laser plane positions; therefore, the accuracy of the proposed method is evaluated with these 1000 laser plane fitting errors. As shown in Figure 13, we kept the camera and two targets of different heights stationary and placed the line laser on a precision displacement table and a rotating table. Then, the controlled line laser moved from 0 mm translational motion to 10 mm, with a step length of 1 mm. The light stripe images of two targets were collected at each position. With the position of 0 mm as the reference, the proposed method was used to calculate the distance between the laser plane at the other 10 different positions and the laser plane at the reference position. Finally, controlling the line laser, we rotated from 0 deg to 2 deg, with a step size of 0.200 deg. Similarly, two target light stripe images were collected at each position, and the angle between the remaining ten different laser planes and the reference position laser plane was calculated using the proposed method, using the 0 deg position as a reference. In order to compare the measurement accuracy of the rail profile after laser coplanar adjustment, a comparative experiment was designed using a standard worn rail as the measurement object. The standard worn rail are shown in Figure 14, with a vertical wear of 11.00 mm. When the laser planes of two laser scan sensors are not coplanar (in order to highlight the effect, the angle between the laser planes is about 3 deg), collect the full-section profile of the steel rail 20 times and calculate its vertical wear. Under the guidance of the proposed method, two laser scan sensors were adjusted to meet the installation requirements of coplanarity. The full profile of the steel rail was collected 20 times again, and its vertical wear was calculated.

**Figure 13.** Schematic diagram of accuracy verification experiment.

**Figure 14.** Standard worn rail with a vertical wear of 11.00 mm. (NO. JZG-MHZ02, horizontal wear 0 mm, vertical wear 11 mm).

The third group of experiments is the repeatability verification experiment. Keep the measuring device stationary, continuously collect 500 pairs of target images through two cameras, and also evaluate the repeatability of the proposed method with the fitting error of the laser plane.

#### *4.2. Experimental Result*

4.2.1. Real-Time Verification Experiment Results

Figure 15 exhibits the laser planes at three typical positions during the adjustment process. It can be seen that the program window displays the laser planes on both sides in real time and calculates the current coplanarity of the laser planes on both sides in real time, according to the standard. Therefore, as opposed to traditional methods that rely on visually observing laser beams, the proposed method allows for the visualization of the three-dimensional laser plane in real time, thus avoiding blindness caused by observing two-dimensional laser beams. Figure 16 shows the corresponding refresh window time of the two cameras. The average refresh time is 0.03 s, and the frame rate is about 33 frames/second, which can fully meet the real-time requirements.

**Figure 15.** The laser planes at three typical positions.

**Figure 16.** The time for the program to refresh the window.

#### 4.2.2. Accuracy Verification Experiment Results

The system calibration module calculates the transformation matrix between the camera coordinate system and the target coordinate system through a single planar target image, with known camera internal parameters. Figure 17 shows the reprojection error diagram. The average reprojection errors of the No. 1 calibration plate and the No. 2 calibration plate images obtained by the left camera are 0.037 pixel and 0.025 pixel, respectively. The average reprojection errors of the images of the No. 2 calibration plate and the No. 3 calibration plate obtained by the right camera are 0.026 pixel and 0.027 pixel, respectively, and the projection errors are within a reasonable range.

**Figure 17.** Reprojection errors of three calibration plates.

The laser plane fitting results and statistical values are given in Figure 18 and Table 1, respectively. It can be seen that, during the random change of the laser plane attitude, the average values of the laser plane fitting errors on both sides obtained by the proposed method are 0.062 mm and 0.054 mm, respectively, which are smaller than the laser plane calibration error of the rail profile measurement sensor.

**Table 1.** Statistical results of the laser plane fitting error (mm).

**Figure 18.** The results of laser plane fitting error.

The results of laser plane distance and angle measurement are shown in Table 2, where MV represents the measured value, AV represents the actual value, and ME represents the measurement error. It can be seen that when the laser plane is moving in a translational motion, the maximum value of the distance measurement error of the laser plane is 0.13 mm, and the average error is 0.09 mm. When the laser plane is rotating, the maximum value of the angle measurement error of the laser plane is 0.019 deg, and the average value is 0.009 deg. The measurement errors of plane distance and plane angle are within an acceptable range.

**Table 2.** Laser plane distance and angle measurement results.

The results of rail wear before and after coplanar adjustment are shown in Figure 19, and Table 3 provides the corresponding error statistics. It can be seen that, after the coplanarity adjustment, the average measurement error of the rail's vertical wear decreased from 0.279 mm to 0.037 mm. Therefore, the proposed method can guide the laser coplanar adjustment process and improve the accuracy of rail profile measurement.

**Figure 19.** Measurement results of rail wear before and after laser coplanar adjustment.

**Table 3.** Measurement results of rail wear before and after coplanar adjustment (mm).

Therefore, the proposed laser coplanar evaluation method has high accuracy. Based on the accurate acquisition of laser plane parameters on both sides, the degree of laser plane on both sides can be quantitatively evaluated, which is superior to the traditional qualitative evaluation method based on visual inspection.

#### 4.2.3. Repeatability Validation Test Results

The repeatability experimental results and statistical values are shown in Figure 20 and Table 4, respectively. It can be seen that when the laser plane is not adjusted, the average values of the fitting errors of the two laser planes obtained by the proposed method are 0.022 mm and 0.021 mm, respectively, and the standard deviations are 0.013 mm and 0.012 mm, respectively, indicating that the method has high repeatability.

**Figure 20.** The results of the repeatability verification experiment.

**Table 4.** Statistical results of the repeatability verification experiment (mm).

#### *4.3. Discussions*

Based on past experience, when the target area accounts for more than a quarter of the entire image, it can ensure that the target attitude evaluation has high accuracy. Conversely, when the target area accounts for a small proportion in the image, a certain degree of measurement accuracy will be lost. Due to the limitations of the on-site installation environment of the rail profile measurement sensor, larger targets cannot be used, and the target area accounts for significantly less than a quarter, which limits the accuracy of laser planar attitude assessment. Therefore, in response to the accuracy issue caused by the small proportion of target areas, the author intends to perform further research work around target optimization and attitude evaluation algorithm optimization in the next step.

Although using this method can ensure laser coplanar installation on both sides of the rail, the next question is how to ensure the laser plane perpendicular to the longitudinal direction of the rail. To address the issue, it is necessary to obtain the longitudinal direction vector of the rail in the camera coordinate system (or world coordinate system). We proposed a solution to calibrate the rail longitudinal direction in our previous research [8]. By combining the longitudinal vector parameters of the rail and the laser plane parameters, it is possible to further quantitatively evaluate whether the installation of laser scan sensors meets the requirements.

#### **5. Conclusions**

Considering the difficulty of evaluating laser plane attitude in the measurement of rail profiles, this paper presents a method for evaluating laser plane attitude using fitting planes, analyzes the impact of non-coplanar lasers on plane profile measurement results, elaborates on how the laser plane attitude calculation is calculated, and constructs the criterion for evaluating laser coplanarity. The experimental results indicate that the proposed method has the advantages of high accuracy, real-time performance, and excellent visualization. This method improves the laser coplanarity adjustment process from the original naked-eye evaluation with low precision, poor real-time performance and subjectivity to computer-based evaluation with high precision, and strong real-time performance and visualization. By doing so, it reduces the calibration error of the traditional calibration method and provides a theoretical basis for improving the accuracy and reliability of rail profile measurement systems. In on-site applications, factors such as changes in ambient light and target posture affect the on-site evaluation efficiency and calibration accuracy of the proposed method. Therefore, it is necessary to improve the evaluation algorithm for railway on-site application scenarios to further improve the efficiency and accuracy of on-site evaluation.

**Author Contributions:** Conceptualization, L.W. and S.W.; methodology, L.W., G.L. and S.R.; validation, N.W., Y.F. and L.W.; formal analysis, S.W.; investigation, S.W.; resources, H.W.; data curation, S.W.; writing—original draft preparation, L.W.; writing—review and editing, S.W.; visualization, L.W.; supervision, H.W.; project administration, Q.H.; funding acquisition, Q.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work described in this paper was partially supported by Science and technology research and development plan of China National Railway Group Co., Ltd., Grant NO. P2021G014; and Research and Development Plan of China Academy of Railway Sciences Co., Ltd., Grant NO. 2021YJ309.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request from the authors.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Communication* **A High-Precision Real-Time Pose Measurement Method for the Primary Lens of Large Aperture Space Telescope Based on Laser Ranging**

**Heng Shi \*, Junfeng Du, Lihua Wang, Jiang Bian, Guohan Gao, Dun Liu, Bin Fan and Hu Yang**

Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu 610209, China **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The aperture of space telescopes increases with their required resolution, and the transmission optical systems with long focal length and diffractive primary lens are becoming increasingly popular. In space, the changes in the pose of the primary lens relative to the rear lens group have a significant impact on the imaging performance of the telescope system. The measurement of the pose of the primary lens in real-time and with high-precision is one of the important techniques for a space telescope. In this paper, a high-precision real-time pose measurement method for the primary lens of a space telescope in orbit based on laser ranging is proposed, and a verification system is established. The pose change of the telescope's primary lens can be easily calculated through six high-precision laser distance changes. The measurement system can be installed freely, which solves the problems of complex system structure and low measurement accuracy in traditional pose measurement techniques. Analysis and experiments show that this method can accurately obtain the pose of the primary lens in real-time. The rotation error of the measurement system is 2 × 10−5 degrees (0.072 arcsecs), and the translation error is 0.2 μm. This study will provide a scientific basis for high-quality imaging of a space telescope.

**Keywords:** pose measurement; laser ranging; space telescope; high-precision

**1. Introduction**

A membrane optical telescope is a new transmission-type imaging space telescope based on special functional materials [1,2]. Theoretically, it overturns the traditional imaging systems, which rely on the primary lens with a curved surface. Technically, it breaks through the bottleneck of increasing the aperture, such as telescope weight, tolerance control of primary lens surface, and envelope size of primary lens [3,4]. Similar projects, such as the Membrane Optic Imager Real-Time Exploitation (MOIRE) program supported by US Defense Advanced Research Projects Agency (DARPA) [5–8], have been called the 21st-century space disruptive imaging technology.

In space, the relative pose (positions and orientations) of the optical elements, especially between the primary lens and the rear lens group, will change due to the changes in the thermal environment and the structure, which will seriously affect the imaging quality of the telescope [9–11]. Furthermore, it is found that the tilt and decentration of the primary lens will make the telescope imaging quality decrease obviously [12,13]. Therefore, the pose measurement of the telescope primary lens has become difficult to overcome urgently in membrane optical imaging technology [14,15].

Traditional pose measurement techniques usually adopt multiple sets of measuring equipment and sensors, which are complicated in structure and require high installation accuracy [16–18]. It is difficult to meet the pose measurement requirements of the primary lens of a large aperture transmission space telescope, and the calculation process is relatively cumbersome [19–21]. Due to excessive system devices and installation errors, more errors will be introduced [22,23].

**Citation:** Shi, H.; Du, J.; Wang, L.; Bian, J.; Gao, G.; Liu, D.; Fan, B.; Yang, H. A High-Precision Real-Time Pose Measurement Method for the Primary Lens of Large Aperture Space Telescope Based on Laser Ranging. *Sensors* **2023**, *23*, 4833. https://doi.org/10.3390/s23104833

Academic Editors: Michele Norgia, Qibo Feng, Jiakun Li and Qixin He

Received: 23 March 2023 Revised: 8 May 2023 Accepted: 11 May 2023 Published: 17 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This paper presents a high-precision real-time measurement method for the pose of the primary lens of a telescope based on laser ranging. First, the high-precision distance information can be obtained using the laser ranging method, and then the high-precision pose information of the primary lens can be obtained by calculation. Finally, the pose of the primary lens can be calculated by using six distances obtained by laser ranging, and the measurement system can be installed freely, which solves the problems of complex system structure and low measurement accuracy in traditional pose measurement techniques. The system has the advantages of a simple structure, convenient installation, and high measurement accuracy.

#### **2. Telescope System Design**

We take a long focal length diffraction imaging optical system as an example; the schematic diagram of the telescope optical system is shown in Figure 1. The diameter of the primary lens is 1 m, the primary lens is a plane, and there are annular microstructures prepared on the surface so that the light can converge on the relay mirror. The relay mirror is 300 mm in diameter, and the relay mirror and the rear lens group are fixed in a frame. The relay mirror is used to fold the optical path and shorten the total length of the system. The rear lens group is responsible for correcting aberrations and providing the final image. The primary lens is extended to the right place by trusses, and the distance between the primary lens and the relay mirror is about 3 m. In order to ensure the imaging quality of the telescope, the Modulation Transfer Function (MTF) @62.5 lp/mm should be up to 0.1. To meet the MTF requirement, the system needs to have good optical quality in terms of resolution, contrast, and uniformity. The MTF is a measure of the ability of the system to transfer contrast from the object to the image, and it depends on many factors, such as the aperture size, aberrations, and diffraction effects, which is a challenging requirement for a system with such a long focal length. Defocus tolerance is important for maintaining the focus of the system over time and under different environmental conditions. The focal depth is proportional to the square of the aperture diameter and the wavelength and inversely proportional to the focal length. Therefore, the defocus of the system should be less than 1/4 focal depth during the whole life cycle [24]. In summary, the long focal length diffraction imaging optical system requires careful design and manufacturing to achieve high MTF and low defocus tolerance. The use of annular microstructures, relay mirrors, and rear lens groups are effective ways to improve the optical performance of the system.

**Figure 1.** Schematic diagram of the telescope optical system.

As we all know, MTF is a key index in the evaluation criteria of an optical system. The index assigned by the telescope system to the MTF drop caused by the change of the pose of the primary lens is 0.98; that is, the MTF drop caused by the change of the pose of the primary lens is less than 2%. Optical analyses show that the decentration of the primary lens is most sensitive to the telescope optical system. Therefore, a tolerance of 1% is assigned to decentration for MTF drop, 0.5% for distance, and 0.5% for tilt. Through

Monte Carlo simulation, it can be known that in order to ensure that decentration caused MTF drop within 1%, decentration should be controlled within ±0.02 mm. Similarly, the distance should be controlled within ±0.1 mm, and the tilt should be controlled within ±0.02◦. According to Shannon's sampling theorem, their corresponding measurement accuracy requirements should be 1/5 of the control requirements. In order to achieve the requirement that MTF drops not more than 2%, the control requirements and measurement accuracy requirements of the pose of the primary lens are provided in Table 1.

**Table 1.** Influence of each degree of freedom of pose on imaging and measurement accuracy requirements.

#### **3. Principle of Measurement Method**

#### *3.1. Establishment of Optical Model of Measurement System*

The method is derived from the principle of the Stewart Gough platform, where the platform and base are connected by six legs at three points on the platform and six points on the base to form a hexapod geometry platform, as shown in Figure 2. The Stewart Gough platform is a mechanical system that uses actuators to drive changes in the length of legs inducing changes in the pose of the platform. For the pose measurement system, the platform legs are the laser optical paths, which constitute the laser optical paths between the primary lens and the relay mirror. In traditional Stewart Gough platforms, forward kinematics is used to map expected changes in platform pose to changes in leg length. For the pose measurement system of the primary lens of the telescope, the laser ranging lengths are known, and the pose of the primary lens must be determined. This is achieved by using an inverse kinematics technique, which generated a sensitivity matrix used to calculate primary lens pose changes from the laser ranging length changes.

**Figure 2.** Hexapod geometry Stewart Gough platform.

In this method, six laser collimators are fixed on the edge of the relay mirror, and three retroreflectors are fixed on the edge of the primary lens. Each retroreflector provides reflected light for the corresponding two collimators so that six laser distances can be obtained. According to the changes in the six measured distances and the initial coordinates of the collimators and retroreflectors, the pose of the primary lens can be calculated. The schematic diagram of the measurement principle is shown in Figure 3, where C1~C6 are the laser collimators and A1~A3 are the retroreflectors.

**Figure 3.** The schematic diagram of the pose measurement principle.

#### *3.2. Pose Calculation Algorithm*

The sensitivity matrix is the basis of the pose calculation technique. The sensitivity matrix defines the sensitivity of the pose changes of the primary lens relative to the relay mirror to changes in the laser length. The change in the laser length is defined as the sensitivity matrix multiplied by the changes in the pose

$$
\overrightarrow{\Delta l} = \mathcal{S} \cdot \overrightarrow{p} \tag{1}
$$

where *S* is the sensitivity matrix, → *p* is the pose vector which includes three translations and three rotations, and → Δ*l* is the laser length change. This is a simple linear algebraic equation, which is the basis of the Stewart Gough platform pose calculation.

The sensitivity matrix is formed using the coordinates data measured by a laser tracker. The laser path direction vectors must be established first. The coordinates of retroreflectors *ai*, and the coordinates of collimators *ci* are measured. The direction vectors *ei* of the laser are created

$$x\_i = \frac{c\_i - a\_i}{|c\_i - a\_i|}, (i = 1, \dots, 6) \tag{2}$$

The sensitivity matrix *S* is created by cross-multiplying the coordinates of the laser collimators *ci* with the direction vectors. The origin of coordinates is located in the center of the primary lens.

$$S = \begin{bmatrix} \mathfrak{c}\_1 \times \mathfrak{e}\_1 & \mathfrak{e}\_1 \\ \vdots & \vdots \\ \mathfrak{e}\_6 \times \mathfrak{e}\_6 & \mathfrak{e}\_6 \end{bmatrix} \tag{3}$$

The sensitivity matrix is a 6 × 6 matrix that defines the sensitivity of the system to the six degrees of freedom of rotation around and translation along the *x*-*, y*-, and *z*-axis of the primary lens.

After the sensitivity matrix is established, we measure six distances from the retroreflectors to the laser collimators *L*0*i* in the initial state. Then, the six distances from the retroreflectors to the laser collimators in real-time *Li* are measured. And the changes in distance are calculated

$$
\Delta l\_i = L\_i - L0\_{i\prime} \, (i = 1, \dots, 6) \tag{4}
$$

In forward kinematics, the sensitivity matrix is multiplied by the desired pose to determine the changes in path length. The inverse of the sensitivity matrix is needed to calculate the changes in the pose from the changes in laser path length. The changes in the pose of the primary lens relative to the relay mirror → *p* are calculated

$$
\overrightarrow{p} = \mathcal{S}^{-1} \overrightarrow{\Delta l} \tag{5}
$$

where *S*−1 is the inverse of the sensitivity matrix, → Δ*l* the changes of laser leg length, which can be obtained by two laser rangefinders.

#### **4. Experiments**

#### *4.1. Laser Rangefinders*

Laser rangefinders are the key equipment in experiments; the stability of their measured data largely determines the accuracy of experimental results. In addition, they have a long working distance range of up to 50 m, allowing for measurements to be taken over large distances. The laser rangefinders used in the experiment are two sets of IDS3005 laser rangefinders from ATTOCUBE Company (Munich, Germany), and the measurement principle is shown in Figure 4a. One rangefinder weighs less than 2 kg, which makes it a significant advantage for use in space telescopes. With a resolution of 1 pm and a repetition accuracy of 2 nm, one rangefinder has three laser channels and can measure three distances of data simultaneously. When two rangefinders are connected to a computer by a switch, six laser channels can be measured at the same time, which constitutes six laser paths of the measurement system.

**Figure 4.** Principle and accuracy test for the laser rangefinder. (**a**) Principle of the laser rangefinder, which is based on Fabry-Perot interference. (**b**) Accuracy test for the laser rangefinder, the laser ranging accuracy is ±0.5 μm.

IDS3005 laser rangefinder adopts the principle of Fabry-Perot interference [25], which has many advantages: The size of the laser collimators can be very small; the beam spacing can be arranged freely; directly measuring the metal surface; directly measuring the Angle by three beams; allowing a large target Angle deviation; little environmental impact; easy installation and debugging; used in vacuum, low temperature, radiation environment. These advantages make the measurement fairly simple. Each channel consists of a telecom fiber that is coupled to a collimator that transmits a beam of laser to the retroreflector. The backward reflection of the fiber tip is the reference arm of the rangefinder, and the measuring arm is the path through the collimator to the retroreflector and back again. The measured optical path difference (OPD) corresponds to the distance between the fiber tip and the retroreflector.

The measurement accuracy of the laser rangefinder was tested. The test was carried out in a standard laboratory environment, where the retroreflectors were fixed at about 3 m away from the collimators of the laser rangefinder, three channels of one laser rangefinder were measured in 5 min, and the measured distance data were recorded by a data recording software, as shown in Figure 4b. According to the measured data, it was determined that the laser rangefinder is accurate to within ±0.5 μm in such experimental conditions, with little variation between two separate rangefinders. This is consistent with the measurement accuracy data given by the manufacturer. Therefore, the error is acceptable and can meet the experimental requirements. Therefore, we can use these instruments for the pose measurement experiment.

Then, in order to determine the pose measurement accuracy under such laser ranging accuracy, the relationship between the laser ranging accuracy and the pose measurement accuracy was calculated by 10,000 Monte Carlo simulations, and the results are shown in Table 2. This experimental system is highly accurate and precise in terms of rotation around the x and y-axes, rotation around the *z*-axis, and translation along the x-, y-, and z-axes. For example, in the experimental system, when the laser ranging accuracy is ±0.5 μm, the accuracy of rotation around the x and y axes (tilt) Rx/Ry is 0.03", the accuracy of rotation around the z-axis (roll) Rz is 0.07", the accuracy of translation along the x and y axes (decentration) is 0.35 μm. The accuracy of translation along the *z*-axis (optical axis direction of the telescope) is 0.15 μm, respectively, indicating highly precise control of the optical axis direction of the telescope. Overall, this experimental system is highly precise and would be ideal for conducting experiments that require highly accurate measurements of position and orientation.

**Table 2.** Analysis accuracy (σ).

#### *4.2. Ground Test*

In order to verify the pose measurement method and its accuracy, an experimental system was built, as shown in Figure 5. For convenience, a frame was used instead of the primary lens of the telescope and its center was aligned with the center of the relay mirror. At the same time, the frame was fixed to a high-precision hexapod platform to provide an accurate reference pose. Three retroreflectors were attached to the edge of the frame. In addition, six laser collimators were fixed to the edge of the relay mirror and connected to the laser rangefinders by six fibers. Each retroreflector and collimator weigh only about 20 g, which makes the overall measurement system very light and flexible.

**Figure 5.** Test system. Three retroreflectors (red circles on the right of Figure 5) were attached to the edge of the frame, and six laser collimators (red circles on the left of Figure 5) were fixed to the edge of the relay mirror to verify the pose measurement method and its accuracy.

After the experimental equipment was installed, the frame was first adjusted to the initial position by the hexapod platform, and the six laser collimators were aligned with the corresponding retroreflectors. Six initial distances were recorded at the initial position, and then the hexapod platform was controlled to drive the frame; the change values of the six new distances were recorded, and the pose can be calculated by the algorithm described above. The pose obtained by this measurement method can be compared with the reference pose of the hexapod platform to verify the accuracy of this measurement system.

To determine the accuracy of rotation and translation of the pose measurement method, the hexapod platform was driven to rotate around the *x*-axis for step movement every 1 × 10−4 degree and move along the *x*-axis for step movement every 0.4 μm. The pose data were measured once a second in 5 min. Figure 6 shows the test results of 1 × 10−4 degrees and 0.4 μm step test of the measurement system, respectively. The error of rotation is 2 × 10−5 degrees (0.072 arcsec), and the error of translation is 0.2 μm. The spikes on the curve are caused by the fact that the hexapod platform has not been stabilized during the motion. The test results of rotation around the *y*/*z*-axis and translation along the *y*/*z*-axis are almost the same. The experimental results show that the accuracy of the measurement system for tilt reaches 1 × 10−5 degree level, and the measuring accuracy for decentration reaches 100 nanometers level, both of which are much better than the measurement requirements. The measurement and calculation were carried out on an ordinary laptop computer, and the single calculation time is 5~10 milliseconds, which achieves the purpose of high-precision real-time measurement.

**Figure 6.** Test results of 1 × 10−4 degrees and 0.4 μm step test of the measurement system. (**a**) The tilt error of the measurement system is 2 × 10−5 degrees, (**b**) The decentration error of the measurement system is 0.2 μm.

#### *4.3. Error Analysis*

During the experiment, there were some errors introduced by the environment and other factors. The analysis shows that the errors are about tens of nanometers. These error sources can be divided into three categories, namely, random errors, systematic errors and environmental errors.

The random errors include detector noise, optical noise, reading error, etc. For the measurement system, the typical value of these noises is about a few nanometers and can be further reduced by averaging the distance data. Systematic errors include some nonlinear periodic errors, which may be caused by optical crosstalk in the light source and electrical crosstalk between reference channels. These errors should be controlled by the design of the rangefinder system, which has little influence on the measurement results. Finally, the environmental errors are mainly caused by ambient temperature and humidity changes and air disturbance, which affect the stability of ranging data. Fortunately, we can reduce these errors caused by environmental factors through the environmental information compensation module of the laser rangefinder and finally control the errors caused by environmental factors within tens of nanometers.

Overall, the errors introduced by the environment and other factors can be mitigated through careful design of the rangefinder system and the use of compensation modules to reduce the impact of environmental factors. By minimizing these errors, the accuracy and reliability of the laser rangefinder can be greatly improved, making it an essential tool for a wide range of applications, from surveying and mapping to industrial positioning and robotics.

#### **5. Conclusions**

In this paper, a reasonable solution for pose measurement of a long focal length optical telescope in orbit is proposed. Considering the simplicity and accuracy, a measurement system is provided to realize high-precision real-time measurement. The step test results show that the rotation error of the measurement system is 2 × 10−5 degrees (0.072 arcsecs), the translation error is 0.2 μm, and the single calculation time is 5~10 milliseconds, achieving the purpose of high-precision real-time measurement. This method has the advantages of simple structure, flexible installation, lightweight, non-contact, high-precision and realtime, which lays a technical foundation for the pose measurement of the primary lens of the space telescope. This measurement method also can be applied to other optical systems and other related systems.

**Author Contributions:** Conceptualization, H.S.; methodology, H.S., J.D. and J.B.; formal analysis, H.S. and L.W.; investigation, H.S., G.G. and D.L.; original draft preparation, H.S.; review and editing, H.S., J.D., L.W., J.B., G.G., D.L., B.F. and H.Y.; funding acquisition, D.L., B.F. and H.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (61905254 and 62075220) and the National Key Research and Development Program of China (No. 2016YFB0500200).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

**Acknowledgments:** We thank the National Natural Science Foundation of China and the National Key Research and Development Program of China. The authors would like to thank the reviewers and the editors for their valuable comments and constructive suggestions that helped to improve the paper significantly.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Sensing Method for Wet Spraying Process of Tunnel Wall Based on the Laser LiDAR in Complex Environment**

**Degang Xu \*, Qing Song \*, Shiyu Fang and Yanrui Guo**

School of Automation, Central South University, Changsha 410083, China

**\*** Correspondence<EMAIL_ADDRESS>(D.X<EMAIL_ADDRESS>(Q.S.); Tel.: +86-0731-8887-9190 (D.X.); +86-187-1104-3216 (Q.S.)

**Abstract:** In tunnel lining construction, the traditional manual wet spraying operation is laborintensive and can be challenging to ensure consistent quality. To address this, this study proposes a LiDAR-based method for sensing the thickness of tunnel wet spray, which aims to improve efficiency and quality. The proposed method utilizes an adaptive point cloud standardization processing algorithm to address differing point cloud postures and missing data, and the segmented Lamé curve is employed to fit the tunnel design axis using the Gauss–Newton iteration method. This establishes a mathematical model of the tunnel section and enables the analysis and perception of the thickness of the tunnel to be wet sprayed through comparison with the actual inner contour line and the design line of the tunnel. Experimental results show that the proposed method is effective in sensing the thickness of tunnel wet spray, with important implications for promoting intelligent wet spraying operations, improving wet spraying quality, and reducing labor costs in tunnel lining construction.

**Keywords:** 3D LiDAR point cloud; intelligent detection methods; normalization; shotcrete

#### **1. Introduction**

With the rapid development of 3D laser scanning technology, the application of laser technology is rapidly expanding and offering decreased costs and increased accuracy. This technology has various applications, including road detection [1,2], object recognition [3,4], surface reconstruction [5,6] and tunnel detection [7,8].

In China, the total length of highways is reported to be 4,846,500 km, including 17,236.1 km of highway tunnels as of 2018 [9]. It is estimated that by 2030, the total number of tunnels in operation in China will reach 17,000, exceeding 30,000 km in length [10]. Therefore, the development of intelligent tunnels in construction and monitoring is becoming increasingly urgent. The use of 3D laser scanning technology in tunnel construction and monitoring can improve construction efficiency, ensure safety, and reduce labor costs.

Laser technology has become an indispensable tool in the intelligent development of tunnels [11,12]. The use of laser scanners for acquiring 3D data of excavation surfaces in tunnels was first proposed by Lemy et al. [13]. They determined the displacement of the excavation surface by comparing the point clouds obtained at different times. This study demonstrated the potential of LiDAR in data collection during tunnel excavation. In subsequent years, Fekete et al. [7,14] used 3D laser scanning for quality control in drill and blast tunnels. LiDAR scanning allowed for precise monitoring of excavation and support installation during construction. More recent research has further explored the potential of LiDAR in tunnel construction activities, such as rock mass classification [15], drill bit wear detection [16], and automatic surveying of tunnel sections [17]. However, research on tunnel shotcrete is still lacking. Shotcrete is commonly used in the construction of railway and highway tunnels [18]. Tunnel shotcrete spraying is a process of rapidly applying shotcrete to the rock or concrete surface to prevent tunnel collapse during excavation [19,20]. Currently, shotcrete spraying is a manual process, resulting in unstable construction quality

**Citation:** Xu, D.; Song, Q.; Fang, S.; Guo, Y. Sensing Method for Wet Spraying Process of Tunnel Wall Based on the Laser LiDAR in Complex Environment. *Sensors* **2023**, *23*, 5167. https://doi.org/10.3390/ s23115167

Academic Editor: Vittorio Passaro

Received: 13 April 2023 Revised: 15 May 2023 Accepted: 22 May 2023 Published: 29 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and slow construction progress. There is a need for further research on the intelligent construction of shotcrete in tunnels.

LiDAR-based 3D object detection is essential for automating tunnel shotcrete spraying, as it directly relates to understanding the tunnel environment. A previous study proposed a novel neural network structure based on LiDAR for detecting the area of tunnel shotcrete spraying [21], demonstrating the real-time monitoring of the tunnel profile and shotcrete area. Ranjbarnia [22] studied the effects of various parameters such as sprayed concrete thickness, soil geomechanical properties, tunnel depth, and fault plane dip angle using the 3D finite difference analysis algorithm and centrifuge physical model, but mainly focused on crossing faults in urban tunnels, which is not applicable for the construction analysis of tunnels in progress. Oreste [23] proposes a new calculation procedure based on the combined use of two calculation methods, convergence confinement and hyperstatic reaction methods, to analyze the factors of shotcrete and determine the trend of the lining safety factor over time. Unlike previous studies, our research focuses on detecting the thickness of shotcrete to enable intelligent detection for large spatial arch spraying processes in complex tunnel scenarios. To achieve this, we developed a LiDAR-based intelligent detection method. Our contribution lies in proposing a novel approach to address the specific challenge of shotcrete thickness detection.

In summary, the main contributions of this paper are as follows:

The structure of this paper is as follows: Section 2 presents the problem and describes the process of data acquisition. In Section 3, an innovative adaptive method for normalizing point cloud data is proposed. Section 4 discusses a fitting method used to obtain a mathematical model of the inner contour of a large-scene tunnel. The experimental results are presented in Section 5. Finally, the paper is concluded with suggestions for further research in Section 6.

#### **2. Material**

#### *2.1. Problem Description*

In modern tunneling operations, it is essential to acquire detailed information about the shotcrete thickness of lining and installed support structures to ensure safety, long-term stability, and quality control. The lining construction [24], including initial and secondary lining, is an important tunnel supporting structure. The initial lining refers to shotcrete after laying steel arch support [25], which is widely used as a support element in underground building construction.

Tunnel shotcrete spraying is a process of spraying concrete at a high velocity onto the surface of a rock or steel structure in order to prevent tunnel collapse during excavation [26]. Currently, shotcrete spraying is a manual process, which involves an operator manually operating the mechanical arm for spraying, as shown in Figure 1. This method is highly dependent on the operator's experience and skills, and the quality of the shotcrete largely depends on the operator. Furthermore, the manual operation exposes workers to a large amount of dust and heavy and dangerous work, which may cause health problems. There-

fore, an automatic tunnel lining construction system, as shown in Figure 2, is necessary to automate the tunnel construction process.

**Figure 1.** Current construction method and classification of wet shotcrete states. The red large circle indicates a magnification of the details in the small circle.

**Figure 2.** The solution of automatic tunnel lining construction system. The red box represents the components of "Intelligent Sensing for Large Spatial Arch Spraying Process", which include real-time point cloud acquisition, point cloud normalization, and intelligent perception for the steel arch. The arrows indicate the specific contents included in each component.

In the automatic tunnel lining construction system, the wet shotcrete process of the tunnel arch is divided into three states: the unsprayed-state, the spraying-state, and the sprayed-state [27,28]. The unsprayed-state is the state in which the shotcrete has not yet been sprayed, the spraying-state is the state in which the shotcrete is being sprayed, and the sprayed-state is the state in which the shotcrete has been applied, as shown in Figure 1. LiDAR is used to collect tunnel point cloud data in these three states, which is then processed using denoising, correction, clustering, and compensation techniques to extract the tunnel arch points.

Afterward, a wet shotcrete thickness description model is established to determine the depth of shotcrete needed for the tunnel arch surface to be sprayed continuously. With the prediction and motion planning of the mechanical arm, an automatic tunnel lining construction system is proposed. This paper focuses on the intelligent sensing for the large spatial arch spraying process in complex tunnel scenarios, which lays the foundation for prediction and motion planning of the system.

#### *2.2. Process of Data Acquirement*

The LiDAR used in this study is mounted on the mechanical arm of the concrete sprayer and integrates 16 laser/detector pairs in a compact housing for data collection [17]. In practice, the LiDAR is designed to be installed on the mechanical arm (as shown in Figure 3a). Since the shotcrete machine is generally oriented towards the working surface of the tunnel, during the wet spraying interval, when the mechanical arm stops moving, the LiDAR scans the tunnel to realize the measurement and perception of the environment. The original point cloud collected in reality is shown in Figure 3b, and there are mainly several issues with it:

(**a**) Data collection site.

(**b**) Data collected.

**Figure 3.** Data collection and analysis. (**a**) describes the data collection site. (**b**) displays a frame of tunnel point cloud data that includes various types of noise. The tunnel arch points are the target of what we want to extract.

To address the issues mentioned above with the collected tunnel point cloud data, an adaptive tunnel arch point extraction algorithm is proposed prior to sensing the shotcrete thickness of the lining.

#### **3. Method**

A point cloud is a vast collection of surface characteristic points of a target object that are directly obtained by LiDAR. However, due to the actual collection site conditions, noise interference is inevitable. Additionally, the presence of construction equipment, operators, rock waste, and other objects at the tunnel site may cause occlusion, which affects the quality of the point cloud to varying degrees. Moreover, the large volume of the tunnel point cloud and the existence of redundant information make it critical to extract the characteristics of the tunnel cross-section point cloud and reduce computational complexity. Therefore, standardizing the tunnel point cloud is crucial during the feature point extraction process to ensure accurate analysis and efficient computation.

#### *3.1. Data Denoising*

To filter out noise data caused by suspended concrete particles, we use a point cloud filtering algorithm based on the threshold neighborhood method proposed by Rusu and Cousins [29]. The algorithm selects neighboring points using a 3D Euclidean distance metric and terminates once a fixed number of neighbors *nnbrs* or all points within a bounding sphere of radius *rsph* have been found. The variables *nnbrs* and *rsph* control the size of the neighborhood selection. Subsequently, the mean *μ* and standard deviation *σ* of nearest neighbor distances are calculated, and points outside the *μ* ± *ϑ* · *σ* range are removed, where the parameter *ϑ* adjusts the sensitivity of the threshold. In our implementation, we have found the optimal value of *ϑ* to be 1, and *nnbrs* = 30.

Furthermore, based on experiments with multiple datasets, we have shown that the *μ* ± *σ* thresholds are effective in removing noise, where around 1% of the points are considered as noise.

#### *3.2. Adaptive Point Cloud Pose Normalization*

During the process of collecting radar data, the tunnel point cloud may be inclined to different degrees due to various installation reasons, as illustrated in the left schematic diagram in Figure 4. To enhance the stability and accuracy of tunnel wet spray state perception, it is crucial to transform the coordinates of the tunnel point cloud before the three-dimensional reconstruction of the tunnel. The pose normalization of the tunnel point cloud needs to be performed to adjust all the tunnel point clouds to the position shown on the right side of Figure 4. This process enables the alignment of the tunnel point cloud with the reference frame, facilitating a more precise three-dimensional reconstruction of the tunnel.

**Figure 4.** Schematic diagram of tunnel attitude standardization.

Considering the rigid body transformation characteristics of the tunnel, solving any part of the rotation transformation matrix can enable the complete rotation transformation of the tunnel point cloud. However, in the tunnel point cloud, the road point cloud exhibits significant planar characteristics, which can be utilized to correct the tunnel wall point cloud by solving the transformation matrix of the road point cloud. To extract the ground point cloud from the tunnel point cloud, we propose a planar extraction algorithm based on the M-estimator SAmple Consensus (MSAC) algorithm [30]. The MSAC algorithm optimizes the calculation method of the loss function of the RANSAC algorithm, addressing the RANSAC algorithm's sensitivity to the threshold *T* selection of interior points. The loss function value of the RANSAC algorithm is represented by *C*1. For the points inside the model, the loss function is 0, whereas for the points outside the model, a constant penalty is incurred. Thus, setting the threshold too high can result in poor estimation, while setting it too low can affect the robustness of the model.

$$\mathbb{C}\_1 = \sum\_i p\_1(e\_i^2) \tag{1}$$

where *e* is the error function, and *p*1 is the robust scale parameter, defined as

$$p\_1(\varepsilon^2) = \begin{cases} 0 & \varepsilon^2 < T^2 \\ \text{constant} & \varepsilon^2 \ge T^2 \end{cases} \tag{2}$$

To address this issue, Torr and Zisserman [31] proposed a new loss function *C*2 that can be minimized to obtain a more accurate model. The *C*2 loss function is defined as

$$C\_2 = \sum\_i p\_2(e\_i^2) \tag{3}$$

where *e* is the error function, and *p*2 is the robust standard error. The *p*2 function assigns a weight to each data point based on its distance to the model, with larger weights assigned to points that are closer to the model. This effectively reduces the influence of outliers in the data and improves the accuracy of the model estimation.

$$p\_2(e^2) = \begin{cases} e^2 & e^2 < T^2 \\ T^2 & e^2 \ge T^2 \end{cases} \tag{4}$$

To extract the plane of the point cloud data, the MSAC algorithm is utilized to obtain the plane equation, as shown in Equation (5). In this algorithm, a fixed penalty is given for out-of-model points, while for in-model points, the fitting effect of the model is considered to establish the most accurate model. Figure 5 shows the fitting result, which demonstrates that the extracted plane points correspond well with the road surface of the tunnel.

$$a\_P x + b\_P y + c\_P z + d\_P = 0 \tag{5}$$

where *ap*, *bp*, *cp*, and *dp* are the coefficients of the plane equation.

**Figure 5.** Fitting result of the MSAC algorithm for plane extraction. The pink plane represents the plane fitted by the algorithm.

In Figure 5, the point cloud is observed to be inclined at a certain angle with respect to the X, Y, and Z axes. However, if the normal vector of the road point cloud is found to be parallel to the Z-axis, it indicates that this part of the point cloud has already undergone the attitude standardization process.

In the cases where the normal vector of the road point cloud is not parallel to the Z-axis, a point cloud correction algorithm based on continuous projection is proposed in this section.

The proposed point cloud correction algorithm based on continuous projection involves a series of steps. Firstly, the point cloud is projected onto the YOZ, XOY, and XOZ planes sequentially. Next, the *α*-shape algorithm [32] is employed to determine the boundary points of the point cloud on each plane, which are then used to identify the center line of the point cloud in the plane. The rotation angle of the point cloud in the plane is determined by the declination angle from the coordinate axis, and the transformation matrix is obtained using the Rodrigue formula [33]. Finally, the point cloud is corrected according to the transformation matrix.

A common way to determine the boundary of a finite point set is through the *α*-shape algorithm [32]. For a finite point set, the algorithm forms a line segment between every two points and draws a circle with a diameter of the line segment. If one of the circles does not contain any other points except for the two points, then the two points are considered as two boundary points. The sum of these boundary points gives the boundary of the point cloud, as shown in Figure 6.

**Figure 6.** The alpha-shape algorithm criteria. The red points represent the point cloud data of the plane fitted by the MSAC algorithm in Figure 5, while the green points represent the boundary points detected by the *α*-shape algorithm. The points *pi* and *pj* on the left are the boundary points, whereas on the right, *pi* and *pj* are the internal points.

After obtaining the projected boundary points, a quadratic function is used to fit the boundary of the point cloud, as shown in Equation (6):

$$
\begin{bmatrix} y\_1 \\ y\_2 \\ \vdots \\ y\_n \end{bmatrix} = \begin{bmatrix} 1 & \boldsymbol{x}\_1 & \boldsymbol{x}\_1^2 \\ 1 & \boldsymbol{x}\_2 & \boldsymbol{x}\_2^2 \\ \vdots & \vdots & \vdots \\ 1 & \boldsymbol{x}\_n & \boldsymbol{x}\_n^2 \end{bmatrix} \begin{bmatrix} a\_q \\ b\_q \\ c\_q \end{bmatrix} + \begin{bmatrix} \boldsymbol{\varepsilon}\_1 \\ \boldsymbol{\varepsilon}\_2 \\ \vdots \\ \boldsymbol{\varepsilon}\_n \end{bmatrix} \tag{6}
$$

where *n* is the number of boundary points; *xi* and *yi* are the coordinates of the *i*-th boundary point; *aq*, *bq*, and *cq* are the coefficients of the quadratic function; and the *i* represent the residual error between the fitted curve and the original point cloud data.

To obtain the midline of the plane, we first use the uniform sampling method to obtain 50 sampling point sets on each boundary. For each sampling point, we record the normal line perpendicular to the boundary and the intersection point set with the opposite boundary. Then, we recalculate the intersection of the normal line at the intersection point on the boundary with the original boundary and record the intersection point set. Next, we take the midpoint of each line segment between the two intersection points and then calculate the midpoint of the two midpoints. The resulting set of midpoints is used to fit the midline using Equation (7). The entire process is illustrated in Figure 7.

$$y = k\_m x + b\_m \tag{7}$$

where *km* and *bm* are the parameters of a curve, and the RANSAC (RANdom SAmple Consensus) algorithm [34,35] is used to estimate the parameters of the model. Then, the inclination angle of the point cloud in the plane can be expressed by Equation (8).

$$\theta = \arctan(k\_m) \tag{8}$$

**Figure 7.** Determination of the midpoints on central axis. The green line represents one of the normal lines of the boundary. *M*1 is the midpoint of the line segment between the intersection point of the normal line with the opposite boundary and the boundary itself. The orange line represents the normal line of the opposite boundary. *M*2 is the midpoint of the line segment between the intersection point of the boundary's normal line with the original boundary and the boundary itself. *M* is the intersection point of the line segments between *M*1 and *M*2, which represents the point on the midline of the plane, and the direction vector is *v*.

Based on the Rodriguez rotation formula, the rotated vector *vrot* of any vector *v*0 in space, rotated by an angle *θ* around a given rotation axis*n*, can be expressed by Equation (9). This equation ensures the accuracy of the rotation operation.

$$\vec{v}\_{\rm rot} = M\_r(\vec{n}, \theta)\vec{v}\_0 = \vec{v}\_0 \cos \theta + (\vec{n} \times \vec{v}\_0) \sin \theta + (1 - \cos \theta)(\vec{n} \cdot \vec{v}\_0)\vec{n} \tag{9}$$

where *Mr* is the rotation matrix, which is defined as shown in Equation (10):

$$M\_{\varGamma}(\vec{u},\theta) = \begin{bmatrix} u^2(1-\cos\theta) + \cos\theta & uv(1-\cos\theta) - w\sin\theta & uv(1-\cos\theta) + v\sin\theta \\ uv(1-\cos\theta) + w\sin\theta & v^2(1-\cos\theta) + \cos\theta & uv(1-\cos\theta) - u\sin\theta \\ uv(1-\cos\theta) - v\sin\theta & vw(1-\cos\theta) + u\sin\theta & w^2(1-\cos\theta) + \cos\theta \end{bmatrix} \tag{10}$$

To investigate the impact of the continuous projection algorithm proposed in this section on point cloud correction, this paper conducts an analysis of the inclination angles and coordinate axes of 1000 frames of tunnel point clouds with varying inclination degrees, before and after correction. The results are presented in Figure 8: "·" represents the tilt angle before correction, and "\*" represents the tilt angle after correction. From the comparison of 1000 point clouds, it can be seen that there were varying degrees of tilt

before correction, with the maximum deviation angle exceeding 15° in all directions. After point cloud correction, the point cloud was corrected well in all directions, with a tilt angle not exceeding 2°. The local enlargement images of 400 frames to 1000 frames further demonstrate the effectiveness of the algorithm correction.

**Figure 8.** Comparison of deviation angle before and after correction. The x-axis represents different point cloud frames, while the y-axis represents the inclination angle.

#### *3.3. Adaptive Point Cloud Wall Normalization*

At the data collection site, construction equipment and workers are in continuous motion. Consequently, the point cloud collection process may result in some degree of obstruction of the tunnel wall, leading to data loss to varying extents. This can lead to subsequent issues in data processing and analysis. To address incomplete information, it is necessary to detect the area of the tunnel wall point cloud and interpolate missing data segments. This process is referred to as the wall standardization process of the tunnel point cloud in this paper.

In order to compensate for missing data in the tunnel wall, this paper presents an adaptive point cloud compensation algorithm based on an interpolation model. The algorithm consists of two main parts: automatic detection and automatic interpolation. During data collection, construction equipment and workers may obstruct the tunnel wall, leading to incomplete point cloud data. In order to address this issue, the proposed algorithm aims to automatically detect the missing areas of the tunnel wall point cloud and interpolate the missing data segments.

The adaptive point cloud compensation algorithm utilizes the first-order difference algorithm for the automatic detection of missing parts of the point cloud. Once the location of the missing point cloud is identified, the algorithm compensates for the missing segment through interpolation using the piecewise cubic Hermite interpolating polynomial (PCHIP) method, which is a type of piecewise polynomial interpolation that uses cubic Hermite polynomials to ensure the smoothness of the interpolated curve.

Assuming that *Pk*(*xk*, *yk*, *zk*) *T* and *Pk*+1(*xk*+1, *yk*+1, *zk*+1) *T* are two points of tunnel arch points after clustering, with a missing area between *Pk* and *Pk*+1, the algorithm checks if |*zk* − *zk*+1| ≥ *δ*, where

$$\delta = \frac{\sum\_{i=1}^{nums-1} |z\_{i+1} - z\_i|}{nums - 1} \tag{11}$$

Here, *nums* is the number of clusters of tunnel arch points, and *δ* is the threshold value for determining the missing area.

The cubic Hermite interpolation polynomial *H*3 is required to satisfy Equation (12):

$$\begin{cases} H\_3(\mathbf{x}\_k) = \varrho(\mathbf{x}\_k), H\_3(\mathbf{x}\_{k+1}) = \varrho(\mathbf{x}\_{k+1}) \\ H\_3(\mathbf{x}\_k) = \phi(\mathbf{x}\_k), H\_3(\mathbf{x}\_{k+1}) = \phi(\mathbf{x}\_{k+1}) \end{cases} \tag{12}$$

where *ϕ*(*x*) is the interpolation function, and *H*3(*x*) is the basis function of the piecewise cubic Hermite interpolation polynomial, which can be expressed as

$$H\_3(\mathbf{x}) = \psi\_k(\mathbf{x})\varphi(\mathbf{x}\_k) + \psi\_{k+1}(\mathbf{x})\varphi(\mathbf{x}\_{k+1}) + \phi\_k(\mathbf{x})\phi(\mathbf{x}\_k) + \phi\_{k+1}(\mathbf{x})\phi(\mathbf{x}\_{k+1}) \tag{13}$$

where *ψk*(*x*), *ψk*+1(*x*), *φk*(*x*), *φk*+1(*x*) is the cubic Hermitian interpolation basis function for nodes *xk* and *xk*+1, and they and their derivative must satisfy Equation (14).

$$\begin{bmatrix} \psi\_k(\mathbf{x}\_k)\_\prime & \psi\_k(\mathbf{x}\_{k+1})\_\prime & \dot{\psi}\_k(\mathbf{x}\_k)\_\prime & \dot{\psi}\_k(\mathbf{x}\_{k+1})\\ \psi\_{k+1}(\mathbf{x}\_k)\_\prime \,\_\prime \psi\_{k+1}(\mathbf{x}\_{k+1})\_\prime \,\_\prime \dot{\psi}\_{k+1}(\mathbf{x}\_k)\_\prime \,\_\prime \dot{\psi}\_{k+1}(\mathbf{x}\_{k+1})\\ \phi\_k(\mathbf{x}\_k)\_\prime \,\_\prime \quad \phi\_k(\mathbf{x}\_{k+1})\_\prime & \dot{\phi}\_k(\mathbf{x}\_k)\_\prime \quad \dot{\phi}\_k(\mathbf{x}\_{k+1})\\ \phi\_{k+1}(\mathbf{x}\_k)\_\prime \,\_\prime \phi\_{k+1}(\mathbf{x}\_{k+1})\_\prime \,\_\prime \dot{\phi}\_{k+1}(\mathbf{x}\_k)\_\prime \,\_\prime \dot{\phi}\_{k+1}(\mathbf{x}\_{k+1}) \end{bmatrix} = \mathbf{I}\_4\tag{14}$$

where **I**4 is a 4 × 4 identity matrix.

Under the constraints of Equation (14), *ψk*(*x*) and *φk*(*x*) can be constructed as follows:

$$\begin{cases} \psi\_k(\mathbf{x}) = \left(a\_\Psi \mathbf{x} + b\_\Psi\right) \left(\frac{\mathbf{x} - \mathbf{x}\_{k+1}}{\mathbf{x}\_k - \mathbf{x}\_{k+1}}\right)^2 \\\\ \phi\_k(\mathbf{x}) = a\_\phi(\mathbf{x} - \mathbf{x}\_k) \left(\frac{\mathbf{x} - \mathbf{x}\_{k+1}}{\mathbf{x}\_k - \mathbf{x}\_{k+1}}\right)^2 \end{cases} \tag{15}$$

The parameters *aψ*, *bψ*, and *aφ* can be obtained by using Equations (14) and (15) under the constraint of satisfying Equation (14). The expressions of *ψk*+1(*x*) and *φk*+1(*x*) can be constructed in a similar way. This results in an interpolation polynomial between the endpoints *Pk* and *Pk*+1, the interpolation results of which are shown in Figure 9. In Figure 9c,d, the red points represent the interpolated data points in the missing segments. It can be observed that the proposed algorithm not only preserves the general trend of the original data, but also achieves better data compensation accuracy with fewer wrongly interpolated data points, resulting in a superior data compensation effect.

**Figure 9.** Compensation of missing parts, the red dots represent the interpolation points of the algorithm. (**a**) shows an area of data loss, where the tunnel wall point cloud is missing. (**b**) shows the markers for the missing parts of the tunnel in the point cloud; specifically, the green asterisk (\*) point represents the starting point of the missing segment, while the blue asterisk (\*) point represents the end point of the missing segment. (**c**) shows the result of piecewise linear interpolation. (**d**) displays the interpolation result of our proposed algorithm.

#### **4. Model**

#### *4.1. Theoretical Basis*

4.1.1. The Gauss–Newton Iteration Method

The cross-section of a tunnel is typically designed as an ellipse [36]. There are two types of algorithms to obtain an ellipse equation: non-iterative and iterative algorithms. Examples of non-iterative algorithms include the Lagrange multiplier-based method proposed in [37]. Iterative algorithms include the Gauss–Newton algorithm-based method introduced in [38,39]. However, due to the complex tunnel environment during the construction phase, the collected point cloud data are not suitable for non-iterative algorithms. Therefore, in this study, Taylor series expansion was used to approximate the nonlinear regression model and improve the regression coefficients by multiple iterations and corrections until the minimum residual sum of squares was achieved. This method is referred to as the Gauss–Newton iteration method.

It is assumed that Equation (16) represents a nonlinear regression model of an elliptical cross-section of a tunnel:

$$\mathcal{G}\_i = f(\mathbf{x}\_i, r) + \varepsilon\_i \ (i = 1, 2, \dots, n) \tag{16}$$

where *r* = (*r*0,*r*1, ··· ,*rnr* ) *T* is an *nr* × 1 matrix of coefficients to be determined, and *εi* represents the error term, which follows a normal distribution. The total number of points to be fitted is denoted by *n*, and *xi* is the *x*-coordinate of the *i*-th point, while *y*ˆ*i* is the predicted value of *xi*.

To obtain an initial value of the regression coefficient *r*, let *g*(0) = *g*0(0), *g* (0) 1 , ··· , *g* (0) *nr T* . Taylor expansion is used at *g*(0) in Equation (16), and the second order and above partial derivative terms are omitted to obtain Equation (17). This approach replaces the nonlinear regression model with a series expansion, and the regression coefficients of the nonlinear regression model are then iteratively updated and corrected until the minimum residual sum of squares is obtained using the Gauss–Newton iteration method.

$$f(\mathbf{x}\_{i\prime}, r) = f(\mathbf{x}\_{i\prime} g^{(0)}) + \sum\_{j=0}^{n\_r - 1} \left[ \frac{\partial f(\mathbf{x}\_{i\prime} r)}{\partial r\_j} \right]\_{r = g^{(0)}} \left( r\_j - g\_j^{(0)} \right) \tag{17}$$

Equation (18) is obtained by combining Equations (16) and (17):

$$\mathbf{y}\_i^{(0)} \approx \sum\_{j=0}^{n\_r - 1} \mathbf{G}\_{ij}^{(0)} \mathbf{b}\_j^{(0)} + \varepsilon\_i, \quad (i = 1, 2, \cdots, n) \tag{18}$$

where

$$\begin{cases} y\_i^{(0)} = y\_i - f(x\_{i\prime}g^{(0)}) \\ G\_{ij}^{(0)} = \left[ \frac{\partial f(x\_{i\prime}r)}{\partial r\_j} \right]\_{r = g^{(0)}} \\ b\_j^{(0)} = \left( r\_j - g\_j^{(0)} \right) \end{cases} \tag{19}$$

Equation (18) can be written in a more simplified matrix form as Equation (20):

$$Y^{(0)} \approx G^{(0)} b^{(0)} + \varepsilon \tag{20}$$

where

$$\begin{cases} Y\_{n \times 1}^{(0)} = \begin{bmatrix} y\_1 - f(\mathbf{x}\_1, \mathbf{g}^{(0)}) \\ \vdots \\ y\_n - f(\mathbf{x}\_{n \times \mathbf{g}}, \mathbf{0}^{(0)}) \end{bmatrix} \\\\ \mathbf{G}\_{n \times n\_r}^{(0)} = \begin{bmatrix} \mathbf{G}\_{10}^{(0)} & \cdots & \mathbf{G}\_{1n\_r - 1}^{(0)} \\ \vdots & & \vdots \\ \mathbf{G}\_{n0}^{(0)} & \cdots & \mathbf{G}\_{nn\_r - 1}^{(0)} \end{bmatrix} \\\\ \mathbf{b}\_{n\_r \times 1}^{(0)} = \begin{bmatrix} \mathbf{b}\_0^{(0)} \\ \vdots \\ \mathbf{b}\_{n\_r - 1}^{(0)} \end{bmatrix} \end{cases} \tag{21}$$

Refine the correction coefficient *b*(0) using the least-squares method:

$$b^{(0)} = \left(G^{(0)T}G^{(0)}\right)^{-1}G^{(0)T}Y^{(0)}\tag{22}$$

The revised values for the regression coefficients *g*(1) can be obtained by using Equation (23):

$$\mathcal{g}^{(1)} = \mathcal{g}^{(0)} + b^{(0)} \tag{23}$$

To update the correction coefficient *b*(*s*) at the *s*-th iteration, we can use the leastsquares method, as shown in Equation (24). Then, the updated regression coefficients *g*(*s*+1) at the (*s* + 1)-th iteration can be obtained by adding *b*(*s*) to *g*(*s*). This iterative process is repeated until the SSR (sum of squares of residual) is below a certain tolerable error *K*, which is given by Equation (25). More specifically, the iterative process continues until *SSR*(*s*)−*SSR*(*s*−1) *SSR*(*s*) ≤ *K*, where *SSR*(*s*) and *SSR*(*s*−1) are the SSR values at the *s*-th and (*s* − 1)-th iterations, respectively:

$$\boldsymbol{b}^{(s)} = \left(\boldsymbol{G}^{(s)T}\boldsymbol{G}^{(s)}\right)^{-1}\boldsymbol{G}^{(s)T}\boldsymbol{Y}^{(s)}\tag{24}$$

$$SSR^{(s)} = \sum\_{i=0}^{n} \left[ y\_i - f\left(\mathbf{x}\_{i\prime} \mathbf{g}^{(s)}\right) \right]^2 \tag{25}$$

where *s* represents the number of iterations.

#### 4.1.2. Analysis of Common Fitting Models

In China, the majority of large-scale tunnels are constructed as arched structures. Consequently, when addressing the problem of fitting the inner contour of a tunnel, arched models, such as the circle, ellipse, and Lamé curve, are commonly employed.

(1) Circle: A circle is the most fundamental geometric shape. For any circular figure with center *Oc*(*xc*, *yc*) and radius *Rc*, the standard equation of the circle is given by

$$(\left(\mathbf{x} - \mathbf{x}\_c\right)^2 + \left(\mathbf{y} - \mathbf{y}\_c\right)^2 = \mathbf{R}\_c^2\tag{26}$$

Based on the findings reported in [40], circular structures are known to exhibit excellent pressure-bearing capacity. Therefore, tunnels excavated using shield machines commonly adopt circular structures.

(2) Ellipse: The actual tunnel environment is complex, and various factors such as geotechnical characteristics and surrounding rock mechanics need to be considered. Additionally, the deformation of the tunnel during use must be addressed. The elliptical structure can adjust its load capacity by changing the eccentricity and is commonly used in practical engineering. The equation for an ellipse is as follows:

$$\frac{\left(\left(x-\chi\_c\right)^2\right)^2}{a^2} + \frac{\left(y-y\_c\right)^2}{b^2} = 1\tag{27}$$

where *Oc*(*xc*, *yc*) is the center of the ellipse, *a* is the semimajor axis, and *b* is the semiminor axis.

The load capacity of elliptical structural tunnels is closely related to the flatness of the ellipse, which can be described mathematically by the eccentricity *e*. In practical engineering, the eccentricity of the ellipse can be adjusted to change its load capacity. When *e* is closer to 0, the load capacity of the ellipse is stronger. Conversely, when *e* is closer to 1, the flatter the ellipse is, and the weaker its load capacity. This relationship between eccentricity and load capacity is important to consider when designing tunnels with elliptical cross-sections.

(3) The Lamé curve: This is also known as the hyperellipse [41], which is an extension of the ellipse. It has been widely used in tunnel engineering due to its adjustable shape parameters and excellent structural performance. The equation of the Lamé curve is given by

$$\left|\frac{\mathbf{x}}{a}\right|^{\eta} + \left|\frac{\mathbf{y}}{b}\right|^{\eta} = 1\tag{28}$$

where *a* and *b* represent the major and minor axes of the Lamé curve, respectively, and *η* is the shape parameter that determines the shape of the curve.

By adjusting the values of *a* and *b*, symmetric closed curves such as rectangles, circles, and hyperellipses can be obtained. When *a* = 6.0 and *b* = 4.0, hyperellipse curves of different orders can be obtained, as shown in Figure 10.

**Figure 10.** Lamé curve.

From Figure 10, it can be observed that when 0 < *η* < 1, the curve is concave inward and takes the shape of a four-pointed star. When *η* = 1, the curve becomes a rhombus. For 1 < *η* < 2, the curve is convex, and the curvature increases as it approaches the vertices. When *η* = 2, the curve becomes an ellipse, which is a circle if *a* = *b*. For *η* > 2, the curve becomes a rectangle with rounded corners, and as *η* increases, it approaches a rectangle, which is also referred to as an ellipse in this case.

Based on the analysis of the circle, ellipse, and Lamé curve, it can be concluded that the Lamé curve is more suitable for fitting the inner contour of the tunnel section, depending on the specific shape of the tunnel. Therefore, this paper adopts the Lamé curve to fit the tunnel section.

#### *4.2. Extraction of Cross-Section Point Cloud*

The calculation of the central axis is a crucial step in obtaining the tunnel section, as it describes the direction of the tunnel, and each section is perpendicular to it. The central axis of the tunnel is typically calculated by projecting the point cloud of the inner wall of the tunnel. There are four commonly used methods to obtain the tunnel center line:

In this study, the fourth method was chosen due to its simplicity and effectiveness, as the boundary lines on both sides of the tunnel can be easily determined. The method involves determining the width of the tunnel by finding the boundary lines on both sides of the point cloud and then shifting each boundary line towards the center. The center line of the tunnel is finally obtained by taking the average value of the shifted boundary lines.

#### *4.3. Fitting of Cross-Section Point Cloud*

Based on the above analysis, it is apparent that the Lamé curve model is more suitable for fitting the inner contour of the tunnel in wet spraying due to its robustness. However, fitting the entire tunnel outline without distinction and estimating the error can increase the amount of calculation and produce a large fitting error due to the different wet spray conditions of each area. To account for these conditions, a method for fitting the inner contour of the tunnel using segmented Lamé curves is proposed in this paper.

Initially, the sequential sampling method is used to take *np* sampling points, and the Lamé curve is fitted to these points using the Gauss–Newton iteration method. The root mean square error (RMSE) of the fitting is then calculated, and if the threshold condition is met, the data are segmented, and sampling is performed from the left and right sides. The RMSE is recalculated until the threshold condition is exceeded, and the data segment fitted by the Lamé curve equation is obtained. Next, m sampling points are selected from the position where the data are interrupted, and this step is repeated until all points have been fitted.

$$RMSE = \sqrt{\frac{1}{n\_p} \sum\_{i=1}^{i=n\_p} (observed\_i - predicted\_i)^2} \tag{29}$$

For the inner contour curve of a specific section during the wet spraying process of the tunnel, the results of the fitting are presented in Figure 11.

**Figure 11.** Piecewise Lamé curve fitting results. The labels 1 to 5 represent the fitting results of five curves. The red points indicate the original tunnel data point cloud, while the other color points represent the fitting curves of different regions.

In Figure 11, the red point cloud represents the actual inner contour curve of a section at a certain position in the tunnel, while the other lines of different colors represent different fitted Lamé curve segments. These curve segments are numbered from left to right, and Table 1 shows the fitting parameters and RMSE of these five curves.

**Table 1.** Piecewise Lamé curve fitting parameters.

After analyzing the fitting parameters in Table 1, it can be observed that the Lamé curve coefficients vary for different sections. However, the root mean square error of each section is within the ideal range, indicating that the model fitting effect is satisfactory.

#### *4.4. Thickness Perception Model*

Figure 12 provides a front view of Figure 13, where the blue curve denotes the crosssection of the tunnel to be wet sprayed, while the red curve denotes the tunnel lining design line.

**Figure 12.** Front view of 3D wet shotcrete thickness description model. The blue curve denotes the cross-section of the tunnel to be wet sprayed, while the red curve denotes the tunnel lining design line.

**Figure 13.** The 3D wet shotcrete thickness description model. The blue curve *L* represents the intersection curve of the unsprayed tunnel surface point cloud and the section, while the red curve *L* represents the intersection curve of the designed surface and the section. Point *Pi* is a point in *L*, while point *P*- *i* is the corresponding mapping point of *Pi* in *L*- . The center of the tunnel at cross-section *y* = *yi* is denoted by *Oc*.

In Figure 13, the depth *di* to be wet sprayed at a point *Pi* on the tunnel section can be observed. The blue curve *L* represents the outline of the tunnel section, while the red curve *L*represents the inner outline of the tunnel lining design.

*Oc* denotes the central point of the section structure, while *P*- *i* is the intersection point between −−→*OcPi* and the inner contour line of the tunnel lining design. The thickness to be wet sprayed at point *P*- *i* is then calculated as *di* = 3 3*PiPi* - 3 3.

During the fitting of the tunnel contour using the segmental Lamé curve, each point *Pi* on the contour was mapped to a point *P*- *i* on the designed inner contour line of the tunnel. Therefore, *P*- *i* is a point on the Lamé curve and satisfies the equation of the Lamé curve for its corresponding segment. From the 3D wet shotcrete thickness description model, we have that the thickness to be wet sprayed at point *P*- *i* is given by *di* = *PiP*- *i* . Hence, *P*- *i* satisfies both the equation of the segmental Lamé curve and the equation of the distance between *Pi* and *P*- *i* . Therefore, we can write the two equations in a system of equations as follows:

$$\begin{cases} \left| \frac{\mathbf{x} - \mathbf{x}\_c}{a} \right|^\eta + \left| \frac{z - z\_c}{b} \right|^\eta = 1\\ \overrightarrow{O\_c P\_i} = \lambda \overrightarrow{O\_c P\_i} \end{cases} \tag{30}$$

where *λ* represents an arbitrary constant.

#### **5. Experiment**

The proposed algorithm has been implemented in a highway tunnel construction project in China, and the experimental results are presented in Figure 14.

**Figure 14.** Experimental results.

As previously described, there are three stages of shotcreting, each requiring a different depth of concrete to be sprayed onto the tunnel surface. These depths are shown in Figure 14 for each respective stage. For the tunnel area in its three states of undried spraying, wet spraying, and dried spraying, we sampled a typical area of 1 m × 3 m × 0.2 m from the experimental results for verification and analysis.

1.37% of the sampling points exceeded 4.5 cm, which is consistent with the on-site working conditions.

The consistency between the unsprayed depth of different states and the actual construction site indicates the reliability of the proposed algorithm.

To evaluate the accuracy of our model, the tunnel arch points extracted in the previous step were manually labeled, as shown in Figure 15.

**Figure 15.** Different types of points in a certain section of the tunnel. The TP points and FP points represent tunnel surface points that were labeled correctly and incorrectly, respectively. The FN points were the points that were falsely labeled as non-tunnel surface points.

In Figure 15, the variables *TP* (True-Positive) marked in blue and *FP* (False-Positive) marked in red indicate the number of points that were labeled correctly and incorrectly as tunnel surface points, respectively. The variable *FN* (False-Negative) marked in green represents the number of points that were falsely labeled as non-tunnel surface points. To evaluate the performance of the model, the precision, recall, and F-score criteria used by Yang et al. [42,43] are adopted.

$$\begin{cases} \textit{precision} = \frac{TP}{TP + FP'}\\ \textit{recall} = \frac{TP}{TP + FN'}\\ \textit{F}-score = 2 \cdot \frac{precision \cdot recall}{precision + recall} \end{cases} \tag{31}$$

Table 2 presents the performance evaluation of the proposed algorithm for tunnel surface extraction and compensation using the criteria mentioned above, including precision, recall, and F-score.

**Table 2.** Assessment of the depth of concrete required to be sprayed on the tunnel surface in six sample areas.

Furthermore, we compared the average precision, recall, and F-score rates of the results obtained by different methods and obtained intuitive comparison results, as shown in Table 3. Based on the experimental results, the proposed method demonstrated higher precision and recall rates than the control group.

**Table 3.** Comparison of precision and recall rate of each method.

The region-growing method [17] is applied to segment the rock surface from the tunnel point cloud. A curvature threshold is used with the region-growing algorithm to extract the rock surface of the tunnel. Additionally, the height threshold is used after DBSCAN to remove the miscellaneous points on the left and right walls of the tunnel. This method achieves an average precision, recall, and F-score of 80.7%, 79.5%, and 80.1%, respectively.

The elliptical cylinder model algorithm [36] differs from region growth algorithms in that it uses the central axis of the fit to divide the region into two parts. Subsequently, the elliptical fitting surface of the tunnel region is obtained through iteration to achieve the filtering of inner wall non-points. This method achieves an average precision, recall, and F-score of 83.7%, 82.9%, and 83.3%, respectively.

The continuous central axis is extracted by 2D projection [44], and then an interpolation algorithm based on quadric parametric surface fitting, using the BaySAC (Bayesian SAmpling Consensus) algorithm, is proposed to compute the cross-sectional point when it cannot be acquired directly from the tunnel points along the extraction direction of interest. This method achieves an average precision, recall, and F-score of 79.2%, 80.1%, and 79.6%, respectively.

By analyzing six sampled areas under different wet spraying conditions, our experimental results indicate that the algorithm achieved an average precision, recall, and F-score of 93.6%, 91.6%, and 92.6%, respectively. This demonstrates that our approach is better suited for analyzing complex wet spraying tunnel walls.

To verify the accuracy of the proposed method in determining the depth of concrete needed to be sprayed on the tunnel surface, a series of experiments was conducted in a real tunnel. The theoretical unsprayed depth *Di* obtained from engineering design was compared to the algorithmic results *Di* as shown in Figure 16.

The absolute average error *MAD* of the 1000 sampled areas is 0.989 cm, which demonstrates that the proposed algorithm meets the accuracy requirements of engineering design and confirms its reliability.

$$MAD = \frac{1}{n} \sum\_{i=1}^{n} |D\_i' - D\_i| \tag{32}$$

**Figure 16.** Comparison of the ideal unsprayed depth and the results obtained by our algorithm. The x-axis represents the sampling area, the y-axis represents the unsprayed depth, and the point marked with a star (✩) represents the ideal unsprayed depth. The point marked with an asterisk (\*) represents the result obtained by our algorithm, and the number next to the point represents the absolute error between the unsprayed depth obtained by the proposed algorithm and the ideal unsprayed depth.

#### **6. Conclusions**

In this paper, we proposed an algorithm for analyzing the area of interest in a tunnel point cloud. We used a continuous projection point cloud correction algorithm to process the attitude of the inclined point cloud during acquisition and developed an adaptive point cloud compensation algorithm to overcome data loss caused by occlusion. Due to the irregular cross-sections of the large tunnel scene and large space, we proposed fitting the segment Lamé curve to mathematically describe the inner contour line of the tunnel, instead of using a simple elliptical cylinder or cylinder model. We then compared the tunnel design line to analyze the thickness of the tunnel to be wet sprayed, allowing for accurate assessment of the tunnel construction. The proposed algorithm has been shown to effectively evaluate the depth of concrete required to be sprayed on the tunnel surface, with an absolute average error of 0.989cm, meeting the precise requirements of engineering design and demonstrating its reliability.

We plan to continue to conduct research on tunnel wet spraying processes, specifically analyzing different scenarios and studying the more complex construction process of long bend tunnels. Our aim is to establish a comprehensive system for monitoring the depth of wet spraying in tunnels, as well as to investigate the positioning issue of mobile LiDAR data to achieve higher real-time detection algorithms.

**Author Contributions:** D.X.: Conceptualization, Funding acquisition, Supervision, Writing—review & editing; Q.S.: Formal analysis, Methodology, Software, Validation, Writing—original draft, Writing review & editing; S.F.: Data curation, Investigation, Project administration; Y.G.: Resources. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China under grant number 61973320, the Key Research and Development Program of Hunan Province under grant number 2022GK2059, and Xinjiang Province under grant number<PHONE_NUMBER>.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

**Conflicts of Interest:** We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Dynamic Photoresponse of a DNTT Organic Phototransistor**

**Marcello Campajola 1,\*, Paolo Di Meo 1, Francesco Di Capua 1,2, Paolo Branchini 3 and Alberto Aloisio 1,2,4,5**

	- 4 CNR-SPIN, Via Campi Flegrei 34, 80078 Pozzuoli, Italy
	- 5 Task Force di Bioelettronica, University of Naples "Federico II", Via Cintia, 80126 Napoli, Italy
	- **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** The photosensitivity, responsivity, and signal-to-noise ratio of organic phototransistors depend on the timing characteristics of light pulses. However, in the literature, such figures of merit (FoM) are typically extracted in stationary conditions, very often from IV curves taken under constant light exposure. In this work, we studied the most relevant FoM of a DNTT-based organic phototransistor as a function of the timing parameters of light pulses, to assess the device suitability for real-time applications. The dynamic response to light pulse bursts at ~470 nm (close to the DNTT absorption peak) was characterized at different irradiances under various working conditions, such as pulse width and duty cycle. Several bias voltages were explored to allow for a trade-off to be made between operating points. Amplitude distortion in response to light pulse bursts was also addressed.

**Keywords:** organic field-effect transistors (OFET); organic phototransistor (OPT); dinaphtho-thienothiophene (DNTT); fast photoresponse; low-voltage operations OPT; random telegraph signals (RTS)

#### **1. Introduction**

Organic field-effect transistors (OFETs) are attracting the interest of the scientific community as a valid alternative to standard semiconductor-based devices. This is due to their several interesting properties, such as cost-effective fabrication over large areas, flexibility, and light weight. The recent progress in OFET fabrication techniques has led to the development of devices with high field-effect mobilities *μ*, low threshold voltages Vth, and high current on/off ratios Ion/off. OFETs can be successfully employed as building blocks in electronic circuits, as well as sensors in applications such as chemical and biological sensing, gas analysis, and pressure monitoring [1–4]. Several works also demonstrated how OFETs can be efficiently used as organic phototransistors (OPTs). The potential of such devices lies in the good photosensitive properties of organic semiconductors, combined with the transistor intrinsic amplification capability [5,6]. Fields of application of OPTs also include optical memory, light communication [7–9], and radiation detection and dosimetry in radiotherapy [10–12]. Several organic semiconductor materials have been investigated for OPT development, including TIPSpentacene and dinaphtho[2,3-b:2- ,3- -f]thieno[3,2-b- ]thiophene (DNTT). For those materials, field-effect mobility comparable to that of amorphous silicon, typically *μ*~1 cm2/Vs, low turn-on voltage, and on/off current ratio larger than 106 have been reported [13–15].

The origin of the photoresponse in OPTs is usually attributed to the increase in the minority carrier trap generation rate under light exposure [16,17]. For a p-type device (i.e., hole-transporting device), the accumulation of photogenerated electrons in trap states in the phototransistor channel causes a shift in the threshold voltage, resulting in multiple holes injected into the channel for a single photogenerated carrier pair. Such a mechanism is known as photoconductive gain. The transverse electric field induced by the gate plays the role of helping to segregate photogenerated electron–hole pairs, decreasing the probability of recombination. The photoconductive gain is recognized as a slow process:

**Citation:** Campajola, M.; Di Meo, P.; Di Capua, F.; Branchini, P.; Aloisio, A. Dynamic Photoresponse of a DNTT Organic Phototransistor. *Sensors* **2023**, *23*, 2386. https://doi.org/10.3390/ s23052386

Academic Editors: Qibo Feng, Jiakun Li and Qixin He

Received: 17 January 2023 Revised: 12 February 2023 Accepted: 17 February 2023 Published: 21 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

phototransistors that exhibit high photoconductive gain usually do not operate as fast as a photodiode with similar carrier transit time [16,18]. Hence, characterizing the dynamic photoresponse to transient light pulses is of paramount importance for devices of interest for real-time applications.

The organic thin-film charge photogeneration has been investigated in the range from picoseconds to seconds after excitation [19,20]. The optical response has been also investigated in depth in OPTs. Among others, we mention here studies on devices based on P3OT [21], PQT-12 [22], PBDFTDTBT [23], P3H3 [24], PDVT-8/PC61BM [25], DNTT [26–28], and TIPS-pentacene [15]. Most of these works studied the performance of OPTs under quasi-static light conditions at different irradiance, whereas a characterization of the device dynamic response to short light pulses at low irradiance has been barely investigated. Quasi-static characterization is useful for comparing the performance of different materials and devices; however, it is not sufficient for assessing their suitability for detection of short light pulses.

In this paper, we present a systematic characterization of a DNTT-based OPT, specifically developed for real-time radiation detection and dosimetry applications [12,29,30], where fast photoresponse to weak, fast light pulses is expected to be a key feature. Specifically, we report on the response to ~470 nm light pulses, under various timing conditions, at different irradiances and operating points. Several figures of merit (FoM), such as photosensitivity, responsivity, and signal-to-noise ratio, have been investigated. In all the explored regions, these quantities strongly depend on the timing characteristics of the light pulses. In this work, we have characterized them as a function of the pulse width, frequency, and duty cycle, studying the correlation with the bias voltages. We also have studied stress and distortion effects in the photoresponse when light pulse bursts are applied, to evaluate the exploitation of such devices in applications where random light pulse sequences are expected.

#### **2. Materials and Methods**

#### *2.1. Organic Phototransistor Layout*

The device under characterization is a thin-film transistor based on DNTT. The device was designed, engineered, and realized at the CNR-IMM laboratory (Rome, Italy). Details on the fabrication process are given in Refs. [29,30]. Here, we only recall some useful information. The OPT was fabricated in a bottom-gate/top-contact configuration on a 100 μm thick substrate of polyethylene-naphthalate.  Source and drain are made by 30 nm thick interdigitated finger electrodes of Au. The semiconductor is 50 nm thick. A dielectric layer, made of a 600 nm thick fluoropolymerbased material (CytopTM), separates the semiconductor from the gate. The whole structure is encapsulated within a 240 nm thick layer of Cytop.

OPT are arranged in a 2 × 2 matrix (see Figure 1a). Each OPT has a channel length L = 100 μm and a width W = 40 mm. The active areaA=L × W is 0.04 cm2.