*2.2. Samples*

Cube concrete samples of 150 × 150 × 150 mm<sup>3</sup> were prepared. They matured in water for 28 days and after another 60 days, smaller samples were water cut (Figure 3) to a size of 120 × 120 × 30 mm3. The cut sample surfaces of 120 × 120 mm<sup>2</sup> were further processed. Sandblasting was chosen from many available surface treatment methods. This method is a popular process for concrete processing in technical applications and it transforms the surface randomly in an isotropic way. Sandblasting was performed using a standardized nozzle, pressure and aggregate grain size in order to repeat the experiment in the future. During the sandblasting with sand grains with d = 0.10–2.00 mm and pressure 0.6 <sup>N</sup>/mm2, an attempt was made to maintain a uniform regime, i.e., constant distance of the nozzle from the surface and constant speed of its movement to obtain even surfaces, as presented in Figure 4.

**Figure 3.** A 60 × 60 mm<sup>2</sup> fragment of the water cut concrete sample surface.

**Figure 4.** Stereoscopic image of the sandblasted surface (together, the left and right images create a 3D effect while keeping eyes as "looking into the distance").

The surface area ratios of the concrete components of the aggregate and cement matrix were estimated on the sand processed surfaces with 2D image analysis. The aggregate and the cement matrix ratios were assessed, as presented in Figure 5 and Table 2, in areas A, B, C and D. The ratio of the aggregate surface varied from 41.81% to 49.48%.

**Figure 5.** The 2D model of the concrete surface presented in Figure 3 (green—acement matrix, white—aggregates). (**A**–**D**): four segments of the sample.


**Table 2.** Sample presented in Figure 5—aggregate and cement matrix areas and ratios.

Afterward, five surfaces of the five samples were used in the research. One of the surfaces, P0, was not processed (mold touching left) and the others (P20, P21, P22, P23) were sandblasted after water cutting.

### *2.3. Background of the Proposed Method*

Based on samples of limited dimensions, the presented method detects surface parameters and extends the results to a larger concrete area. The developed approach includes a combination of the following methods:


In this paper, determination of the surface parameters, i.e., point cloud coordinates and profile height (*x,y,z*), was obtained using close-range photogrammetry. The entire procedure of the assessment of the concrete shear strength parameters is presented in Figure 6.

**Figure 6.** Block diagram of the shear strength parameter assessment.

### *2.4. Data Acquisition Employing Photogrammetry*

Terrestrial close-range photogrammetry (CRP) is widely used in civil engineering. Professional photogrammetric cameras equipped with high-quality lenses take metric photos with known orientations in the adopted reference system. These photos include elements of internal orientation, which allows researchers to read the coordinates of the photographed objects. Based on these coordinates, it is possible to control the shape of structural elements and their location in space. Nowadays, another approach is emerging in photogrammetry which involves the use of the so-called non-metric cameras with cheaper lenses that are characterized by greater distortions (e.g., radial, tangential). Geometric distortions are minimized in the camera calibration process based on the analysis of many photos presenting special calibration charts (e.g., black and white chessboard), as presented in Figure 7a. The lack of elements of internal and mutual orientation is replaced by a large number of photos taken, which overlap to a large extent (covering almost the same area as the object). Automatic detection of the tie points on adjacent photos allows researchers to estimate the camera position and orientation in space. It allows them to build a 3D model of the measured object. To obtain the correct scale of the model, as well as to optimize the estimated camera position and provide appropriate georeferencing, the control points are established on the object. The known coordinates of the control points (usually determined using an electronic total station) also provide the possibility of checking the accuracy of the created 3D model of the object. The obtained model's accuracy depends on many of the following factors: size of the object, shooting distance, camera matrix resolution and lens parameters, number of control points, their arrangemen<sup>t</sup> and coordinate accuracy, as well as the lighting, colour and texture of the photographed object.

**Figure 7.** Stages of close-range photogrammetry (CRP) model preparation: (**a**) black and white chessboard for camera calibration; (**b**) a sample on the stand with markers (six coded control points); (**c**) estimated locations of the camera (blue rectangles); (**d**) dense point cloud obtained for the concrete sample.

In the literature, there are many examples of photogrammetry applications in construction and civil engineering. The usefulness of CRP in crack detection and the accurate 3D modelling of medieval bridge elements as well as for exterior material characterization were analysed in [25]. An example of a CRP application for the deformation analysis of a concrete beam during a loading experiment is presented in [26]. Precise estimation of the vertical deflections and horizontal displacements of the concrete beam during the load test in laboratory conditions showed that the photogrammetry technique was capable of monitoring both static and dynamic deformations [27]. Nowadays, photogrammetry is increasingly used for measurement of the vibrating structures. The most current trends in this technique are point tracking, digital image correlation and targetless approaches, as well as the integration of photogrammetry with other measurement techniques used in structural dynamics (e.g., interferometry) [28]. Both the point tracking technique and the development of a 3D model of the object based on a set of photos require an e ffective algorithm for recognizing the control point markers on the pictures. The aforementioned desire to integrate various measurement techniques usually involves combining photogrammetry with laser scanning. A coupled 3D laser scanning and digital image correlation system for the geometry acquisition of a railway tunnel was proposed and tested in [29]. When integrating measurement methods, attention should be paid to the accuracy achieved by the individual techniques and, in particular, this accuracy should be checked in real conditions prevailing at the construction site [30].

In the presented research, the photos of the concrete surface samples (120 × 120 mm2) were taken with a Nikon D800 non-metric camera (Nikon CEE GmbH, Warsaw, Poland) with a 50 mm single-focal-length lens with a resolution of 7360 × 4912 pixels. To calculate the pre-calibration parameters (in Agisoft Lens software, Version 0.4.2 beta 64 bit (build 2399), Agisoft LLC, St. Petersburg, Russia [31]), 10 photos of a black and white chessboard were taken (Figure 7a). Next, for each concrete sample, several dozen photos were taken using a special flat stand and six markers (control points) (Figure 7b). The coordinates of the control points in the local coordinate system were determined using the least square method (LSM) through the adjustment of 15 linear measurements. The horizontal position accuracy of the control points after adjustment was no greater than 1 mm. For each concrete sample, the tie points on the photos were found with the highest accuracy requirement and all photos were aligned in the Agisoft PhotoScan Professional software (Version 1.2.4 64 bit (build 2399), Agisoft LLC, St. Petersburg, Russia) [32]. Examples of the estimated positions of the camera when pictures of the concrete sample were taken are presented in Figure 7c as blue rectangles. The obtained tie points were filtered using the values of reprojection error, reconstruction uncertainty and projection accuracy. Based on control points, the optimization of camera alignment, as well as the scaling and georeferencing of the samples, were performed. The last step was to generate a dense point cloud for each concrete sample, an example of which is presented in Figure 7d.

The CSV format files were generated from the results of the photogrammetric measurement of the concrete samples. Each file contained a collection of points representing the sample surface (3D model) supplemented with colour information components in the form of records (*<sup>x</sup>*,*y*,*z*,R,G,B). The collected *z* coordinates of the points were normalized in reference to an average value of *z* or in reference to a surface trend if necessary.

### *2.5. Geostatistical Models of the Sample Surfaces*

To fit the theoretical semivariograms to empirical data from 3D surface models, geostatistical methods were used. The semivariograms measured the spatial autocorrelation of the acquired sample points. The measured height, i.e., the "*z*" values, were used to build an empirical semivariogram with the LSM and the Gauss–Newton algorithm as a nonlinear fitting method. The general form of the semivariograms according to Equation (2) was used [33], where for each *h*, half of the mean value of the squared di fference (*z*(*d*) − *z*(*d* + *h*)) is defined as semivariance with a squared length unit.

$$\gamma(h) = \frac{1}{2N} \sum\_{i=1}^{N} \left( z(d) - z(d+h) \right)^2 \quad \{d, d+h\} \in \mathcal{P} \tag{2}$$

where *z(d)* are normalized "z" coordinates at location *d, h* is the lag distance, i.e., distance that separates the two analysed locations and *N*—number of tested pairs {*d, d* + *h*} from space *P*.

For any two locations in a small distance on a surface, a small value of a measure (*z*(*d*) − *z*(*d* + *h*))<sup>2</sup> is expected. With the increasing lag distance *h* between points, the similarity of the measure decreased. This natural behaviour could be described by various theoretical relationships, as presented in Table 3 and Figure 8. The nugge<sup>t</sup> model is used for discontinuous characteristics or areas of limited local extend. Spherical, exponential, Gaussian and linear-plateau models are important among the commonly used semivariograms in engineering applications.

**Model Name Semivariogram Equation No.** nugge<sup>t</sup> γ(*h*) = 0 ⇒ *h* = 0 *s* ⇒ *h* > 0 (3) linear with sill γ(*h*) = *sh r* ⇒ *h* ≤ *r s* ⇒ *h* > *r* (4) spherical γ(*h*) = ⎧ ⎪⎪⎨ ⎪⎪⎩ *s* #1.5 *h r* − κ - *h r* 3 \$ ⇒ *h* ≤ *r s* ⇒ *h* > *r* (5) exponential γ(*h*) = s - 1 − e −h r (6) logarithmic γ(*h*) = 0 ⇒ *h* = 0 *s* log(*h* + *r*) ⇒ *h* > 0 (7)

**Table 3.** Basic semivariogram models.

*r*—range [L], typically constant value limiting the zone of mutual correlation of points in the model, *s*—sill is the model constant [L2], *h*—lag distance, variable in the model function [L], κ—model constant typical as 0.5 [·].

**Figure 8.** *Cont.*

**Figure 8.** Theoretical semivariogram models (graphs): (**a**) nugget; (**b**) expotential; (**c**) linear; (**d**) logarithmic; (**e**) spherical.

Notation and semivariance fittings are presented in Figure 9, where the empirical points and fitted model are presented. When a pair of points appears in an interval, the lag distance between them increases linearly. This concept introduces a statistical measure of "z" variability based on a distance between points which is di fferent from subjective deterministic measures. Moreover, some specific features of semivariograms are commonly used to describe natural observations:


**Figure 9.** Semivariogam fitting and notation.

#### *2.6. Geostatistical Model Fitting to Empirical Semivariogram—Surface Semivariogram Model (SSM)*

Semivariogram modelling governs a step between spatial description and spatial prediction. Geostatistical analysis provides many semivariogram functions to model the empirical semivariogram, as presented in Table 3 and Figures 8 and 9.

The selection of a theoretical model and its fitting procedure is crucial to ge<sup>t</sup> a satisfactory prediction of unsampled points. The maximum likelihood or least square regression (including the weighted version) can be used to fit the experimental semivariance. A small and comparable

sum of the deviations indicates comparable performance of the models. Therefore, the selection of the semivariogram model is a prerequisite for better performance. It is possible to have several models to choose from. The selection of a satisfactory model requires the balancing of goodness of fit and model complexity. Taken together, a likelihood ratio approach leading to a chi-square test was used in this study.

Semivariograms were prepared with two methods for the CRP normalized surface models:


At this stage, theoretical semivariograms were fitted to the empirical semivariograms with the use of LSM with uniform weights, as presented in Figure 10.

**Figure 10.** Algorithm No. 1 in R language. Data collection and semivariogram fitting to a sample surface (# starts description line).

Algorithm 1 in the R language [22] contains the main part of the fitting code for the spherical model, based on the *Gstat* package [34,35] and SP spatial pack library [36]. The fitted model according to Equation (5) is fully described with two parameters: *r* (range) and *s* (sill), for the zero value of the nugge<sup>t</sup> effect. Both parameters were subsequently used in the process of the generation of a random field described in the next section.

### *2.7. Generation of the Random Field*

As mentioned before, semivariograms were used to perform simulations of both the profile and the surface by using the random field theory for the Gaussian process, correlated with the semivariogram. For the field generation, the sequential simulation algorithm was used and coded in the R language [22], as presented in Figure 11. It was assumed that the modelled phenomenon could be described employing the Gaussian random field. The use of the Gaussian random field in Euclidean space assumes that a random process is stationary and isotropic. The symmetrical non-negative correlation function and zero mean value were also assumed. Moreover, the points of the generated samples were evenly distributed on the surfaces and on the profiles.

The random field generator presented in Figure 11 was also based on the *Gstat* package [34,35]. The sequential algorithm operates well on fields of large dimensions (i.e., more than 10<sup>7</sup> points). To approximate the conditional distribution at a location, this effective method uses only data and simulated values from a local neighborhood. The process is described by parameter *nmax*. The larger the *nmax,* the better approximation is with the time consumption increasing exponentially. The selection of the nearest *nmax* data or previously simulated points is performed in *Gstat* with a bucket quadrate neighbourhood search algorithm [37].

**Figure 11.** Algorithm No. 2 in R language. Random field generator. (# starts description line).

In the sequential Gaussian simulation, a random path through the locations was assumed because with a regular path, the local approximation of conditional distribution (using only the *nmax* nearest) may have caused spurious correlations. *Gstat* reuses expensive results (neighborhood selection and solution to the kriging equations) for each of the subsequent simulations when multiple realisations of the same field are requested. There is a considerable speed gain in the simulating of multi-fields in a single call compared to several hundred calls, each for simulating a single field. The random number generator used was the native random number generator of the R environment. Besides this, the reproduction of the sampling results (i.e., fixing randomness) was assured by using the random number seed with the *set.seed()* function.

As a result of the random field generation, a cloud of points with the coordinates (*x,y,z*) was obtained. Examples of the various generated random profiles based on the spherical model, that fits well with relatively even surfaces, are presented in Figure 12.

**Figure 12.** Examples of virtual profiles based on the spherical model with parameters: (**a**) *s* = 0.01 mm2, *r* = 1 mm, (**b**) *s* = 0.1 mm2, *r* = 5 mm, (**c**) *s* = 0.01 mm2, *r* = 5 mm, (**d**) *s* = 0.1 mm2, *r* = 1 mm. All profiles are supplied with histograms of *z*. (Horizontal axis - index is a vector of distance in 100 μm).

### *2.8. Surface Roughness Parameters*

The random field containing the "*z*" coordinates obtained from the generated virtual surface could be used to calculate any traditional roughness parameters such as average roughness (*Ra*), mean peak height (*Rpm*), mean valley depth (*Rvm*) and many others. However, to explore the available experimental results, [5,10] the mean valley depth (*Rvm*) was used. *Rvm* is defined as the average of the maximum valley depth from each sampling length and is given by:

*Rvm* = 1 *n n i*=1 *vi* , (8)

where *vi* is the maximum valley depth at each sampling length *i* and *n* is the number of equal sampling lengths.

Then, a relationship between the roughness parameter (e.g., *Rvm*) and semivariogram parameter sill *s* (Figure 9) was obtained with the regression analysis algorithm presented in Figure 13. Evaluation of the roughness parameters is not a necessary step when full test data (i.e. (*<sup>x</sup>*,*y*,*<sup>z</sup>*) coordinates) are available because semivariogram parameters can be directly correlated with the strength test results.

**Figure 13.** Algorithm No. 3 in R language. Nonlinear least squares (NLS) estimates of the s parameters. NLS use as a default the Gauss–Newton algorithm.

### *2.9. Concrete Interface Strength*

The final step in the presented method was the transformation of the roughness parameters and, more importantly, the semivariogram parameters into contact bond strength. Factors *c* and μ are required (Equation (1)) if [3] is used. According to the results presented in [5,10], relationships between the mean valley depth (*Rvm*) and the factors are expressed by the following:

$$\mathcal{L} = 1.062 \, R\_{\text{v}\text{m}}^{-0.145} \text{\textdegree} \tag{9}$$

$$
\mu = 1.366 \, R\_{\text{vm}}^{-0.041} \,\text{s} \tag{10}
$$

However, relationships between the semivariogram parameter *s* (sill) and the *c*, μ factors can be achieved with the method presented in Figure 13.

### **3. Results and Discussion**

The five concrete surface samples (P0, P20, P21, P22, P23) described in Section 2.2 were examined with the use of CRP. Based on the control points (Figure 7b), the optimization of the camera alignment, as well as the scaling and georeferencing of the models, were performed. The final number of tie

points, as well as the obtained model accuracy, are summarized in Table 4. The total error calculated based on six control points did not exceed 0.32 mm and the total error of model scaling calculated based on six scale bars did not exceed 0.16 mm for all five concrete samples. The last step of CRP was to generate a CSV file with the dense point cloud for each concrete sample, as presented in Figure 14. The files contained sets of coordinates and color records (x,y,z,R,G,B). Discrete projection of the surface with a resolution of 100 points per 1 mm (2540 dpi) in the sample plane results in a single file size of approximately 40 MB.

**Table 4.** The number of tie points and assessment of obtained accuracy for the CRP models of the concrete surface samples.


**Figure 14.** The perspective view of the segments of point clouds obtained from CRP as models of concrete surfaces for samples P0, P20, P21, P22 and P23 (the scale bars are given in meters).

The resulting three-dimensional surface models of the concrete samples were slightly jagged at the edges. For this reason, the concrete sample models were cut symmetrically on each side to a final size of 110 mm × 110 mm and further calculations were performed on the cut samples with the densities presented in Table 4.

For each sample, theoretical semivariograms were fitted to the empirical data with the use of LSM with uniform weights (Algorithm 1 in Figure 10). The "z" coordinates were normalized in reference to the average value and the final semivariogram type was selected based on the shape of the empirical semivariogram, as presented in Figure 15 and Table 5. This selection included the theoretical semivariogram type and method of point selection (each pair from surface or combination along parallel lines).

**Figure 15.** Empirical and theoretical semivariograms for samples: (**a**) P0, (**b**) P20, (**c**) P21, (**d**) P22, (**e**) P23.


**Table 5.** Semivariogam fitting results.

1 cop-x—combination of points along parallel lines, along the x-direction; 2 cop-x45—combination of points along parallel lines, along 45 deg to the x-direction; 3 cop-y—combination of points along parallel lines, along the y-direction.

It is noteworthy that the surface processing type did not change the semivariogram model because the models for samples P0 and P20 were the same and the semivariogram model changed among sandblasted surfaces. This suggests that factors other than the observed parameters can contribute to the semivariogram type.

The theoretical semivariograms presented in Table 5 were subsequently used to generate virtual surfaces, as presented in Figures 16 and 17. Furthermore, from these results, the regression models of the mean valley depth *Rvm*, cohesion factor *c* and friction factor μ were approximated as a power function of sill *s*:

$$R\_{\rm vm} = B\_1 \cdot \mathbf{s}^{A\_1} + \mathbf{C}\_1 \boldsymbol{\omega} = B\_2 \cdot \mathbf{s}^{A\_2} + \mathbf{C}\_2 \boldsymbol{\mu} = B\_3 \cdot \mathbf{s}^{A\_3} + \mathbf{C}\_3. \tag{11}$$

where (A, B, C) is a set of fitting vectors.

**Figure 16.** Single random field realisations of the original surfaces P0, P20, P21, P22, P23 (scale values are *z* coordinates in mm).

**Figure 17.** A 3D image of the random field realisation of the original surface P20.

In Table 6, the regression parameters can be found, and in Figure 18, the fitting results of factors *c* and μ can be found.


**Table 6.** Parameters of regression models *Rvm*, *c* and μ.

**Figure 18.** Results of fitting cohesion—*c* (**a**) and friction—μ (**b**) factors to the sill—*s* parameter.

The estimation of factors *c* and μ as functions of *s* changes Equation (1) to:

$$\upsilon\_{\rm Rdi} = c(s) \cdot f\_{\rm cld} + \mu(s) \cdot \sigma\_{\rm n} + \rho \cdot f\_{\rm yd} \cdot (\mu(s) \cdot \sin a + \cos a) < 0.5 \,\upsilon \, f\_{\rm cd} \,\tag{12}$$

where *vRdi* is the shear strength at the concrete-to-concrete interface; *c(s)* and μ*(s)* are the factors that depend on the interface geostatistical parameter *s* model; ρ is the reinforcement ratio; σ*n* is the external normal stress acting on the interface; *v* is a strength reduction factor; and *fcd* is the concrete compressive strength.

If the reinforcement is absent, Equation (12) leads to:

$$
\omega\_{Rdi} = \mathcal{c}(s) \cdot f\_{\text{ctd}} + \mu(s) \cdot \sigma\_n \tag{13}
$$

or

$$\frac{\upsilon\_{Rdi}}{f\_{ctd}} = c(s) + \mu(s) \cdot \frac{\sigma\_{ll}}{f\_{ctd}}.\tag{14}$$

Equation (14)'s plots are presented in Figure 19 for the range of normal stress σ*n* ∈ [0; 0.5] *fctd*. This type of plot could be included in any formal design document.

**Figure 19.** Plots of the relative interface shear strength as a function of the geostatistical sill—s parameter. (**a**) half-logaritmic scale of sill—*s*; (**b**) linear scale of sill—*<sup>s</sup>*.

In the absence of normal stress <sup>σ</sup>*n,* the value of *vRdi*/*fctd* gives the cohesion factor *c* value which, according to [3], varies from 0.025 to 0.4 (and even 0.5 for indented surface). Since the surface of the sample P0 was not processed, it could be regarded as smooth or even very smooth according to [3] *c* = 0.025 − 0.20. The minimum value of *c* in Figure 19 is greater than 0.2 which indicates that the values of *c* in [3] might be too conservative. However, where the factor μ is concerned, it does not change very much with surface texture, according to [3] (0.5–0.7 and 0.9 for indented surface). Considering that <sup>σ</sup>*n*/*fctd* = 1, the plots in Figure 19 give *c* + μ values, which for the smooth or very smooth surface is 0.525–0.8. However, in Figure 19 the minimum value of *c* + μ = 1.1 which again might be too conservative.
