## Abstract

A method that estimates invisible cracks from the surface based on the surface deformation measured by digital image correlation (DIC) is being developed. An inverse problem is setup to estimate such invisible cracks from the surface deformation. Surface deformation, measured by the DIC method, contains noise. Inverse problems have ill-conditions. The regularization method applied in this study is an extension of the joint estimation maximum a posteriori (JE-MAP) method. The JE-MAP algorithm alternates between MAP method estimation and the grab-cut (GC) method to avoid ill-conditions. The physical constraints on displacement and the forces at the cracks and the crack perimeters (ligaments) are added to the MAP method. The displacement and load at the cracks and the ligaments have a cross-sparse relationship. The MAP method estimates the displacement or the load at the cracks and the ligaments. The estimated result varies greatly at the boundary between the cracks and the ligaments. This boundary is determined by the GC method based on the estimated result. This study amplified the changes at the boundary between the cracks and the ligaments in the estimated results. The amplified results were input into the GC method to improve the boundary-determination accuracy. The regularization method developed from the JE-MAP method was combined with DIC method to estimate the cracks in invisible locations. The method proposed in this study estimated cracks more accurately than L1-norm regularization in inverse problems where the observed data were strain distributions measured by the DIC method.

## 1 Introduction

Large electrical devices are regularly inspected for cracks to ensure their integrity. The inspection of cracks in invisible locations requires the disassembly of equipment, resulting in long shutdown periods. If cracks in invisible locations could be inspected without disassembling the equipment, the periods of equipment shutdown could be shortened.

Ultrasonic testing and X-ray inspection are widely used to inspect cracks in invisible locations. In ultrasonic testing, a probe is placed in contact with the surface to be inspected, ultrasonic waves are emitted and those reflected from a crack are measured to identify the crack's shape [1]. Reducing the size of the inspection system is difficult because the probe must generate and receive ultrasonic waves. In an X-ray inspection, an object is irradiated and a crack's shape is identified from the difference in the transmitted X-dose [2]. The inspection device is large and the number of objects to be inspected is limited. To inspect electrical equipment without disassembling it, the inspection device must be placed in a narrow space inside the equipment and be as small as possible to avoid damaging the equipment by making contact.

Camera-based visual inspection is a noncontact inspection method that can be easily miniaturized. In combination with the digital image correlation (DIC) method [3], which has attracted much attention in recent years, a camera can measure the deformation caused by a crack in an invisible area. A method was proposed to estimate cracks in invisible locations from the deformation measured by DIC [4]. In the above method, for each crack's shape, deformation data at visible locations are prepared from which a crack shape close to the measured deformation is selected. The estimation accuracy depends on the amount of data that must be prepared in advance, and it is difficult to prepare every possible crack shape.

The problem of estimating cracks in invisible locations from surface deformation encounters the inverse problem of estimating the displacement distribution of crack surfaces from surface deformation. The method for solving inverse problems with inadequacies is to predict the crack shape using topology optimization. The research using topology optimization proposes to leverage full-field response data obtained by DIC in a topology optimization framework to reconstruct the internal damage in members [5]. Furthermore, the above authors also propose the detection of interior anomalies of structural components, inferred from the discrepancy in constitutive properties such as elasticity modulus distribution of a three-dimensional heterogeneous/homogeneous sample, from limited full-field boundary measurements using three-dimensional-DIC [6]. The above authors then extend their research to demonstrate the feasibility and performance of the proposed method for a set of large-scale structural steel beams with and without buried defects using a full-field three-dimensional DIC sensor approach [7]. The above authors also propose a new strategy, to get a unique solution for a finite element model updating problem in detecting internal properties of a structure by topology optimization method, a novel strategy is proposed termed as octree partitioning algorithm [8]. The proposed method of using topology optimization is high computational costs because the finite element method is computed iteratively in the optimization.

Since the measured deformations contain noise, the inverse problem of estimating invisible cracks from the deformation is an ill-condition. In order to overcome this problem, Amaya et al. [9] proposed regularization by taking into account the physical relationship between the surface deformation gradient and the crack deformation and a physical constraint between the displacement and the force at a crack and its perimeter (ligament). In this study, the regularization method is the joint estimation maximum a posteriori (JE-MAP) method [10], which was developed to identify the location of body tissues from X-ray CT results. The JE-MAP algorithm alternates between MAP (method) estimation [11], which incorporates physical constraints on the a priori information, and the grab-cut (GC) method [12], an image segmentation algorithm, to avoid ill-conditions.

Therefore, in this study, a method to estimate cracks in invisible locations by combining a regularization method developed from the JE-MAP and a DIC method is being developed. The regularization method first adds to the MAP method the physical constraints of displacement and force at the cracks and ligaments. The physical constraint has a cross–sparse relationship between the displacement and the load at the crack and ligament, as described below. The grab-cut method is also improved to easily determine the boundary between cracks and ligaments based on the displacement or load predicted by the MAP method. Before being input to the grab-cut method, the displacements or loads are converted into an image that highlights the crack and ligament boundaries. We show the results of the estimation of cracks in invisible locations with a combination of the developed regularization method and a DIC method. First, an inverse problem is setup to estimate cracks. Next, the regularization method that extends the JE-MAP method is described. The specimen to be deformed is modeled using the finite element method to identify the relationship between the surface deformation and the crack for solving the inverse problem. The discretization method of the obtained results is described. In addition, the deformation distribution of specimens with and without cracks is measured by DIC method, and the observed data of the inverse problem are obtained from the deformation distributions. Finally, the cracks are estimated from the observed data using the existing and regularization methods, which are extensions of the JE-MAP method. The availability of this study is tested by comparing the estimation results with L1-norm regularization [13].

## 2 Regularization Method of the Inverse Problem

### 2.1 Inverse Problem Setting.

The inverse problem in this study is to estimate a crack propagating in the thickness direction from one side of a flat plate from the strain changes on the crack-free side. Figure 1 shows a schematic diagram of this study's inverse problem.

A virtual plane, which denotes the location of a crack to be estimated, is shown by a dashed flat line in Fig. 1. The virtual plane's width direction is the *X*-axis, its thickness direction is the *Y*-axis, and the length of the steel plate is the *Z*-axis. A uniform tensile load is applied in the *Z*-axis direction. The observed surface is where the strain changes are measured, shown in green in Fig. 1. The observed data of the inverse problem are the strain changes in the *Z*-direction of observed surface $\epsilon Z$. The unknown data of the inverse problem is the displacement *u* of the virtual plane in the *Z*-direction. The *Z*-direction strain of the virtual plane from the tensile loading changes with the presence or absence of a crack (Fig. 1). The Displacement *u* at the crack of the virtual plane is larger than at the ligament.

Here $\epsilon Z\u0303$ and **u** in bold are discretized vectors. $\beta $ is the discretized vector of the measurement errors. **H** is a constant matrix representing the relationship between $\epsilon Z$ and **u**. Constant matrix **H** is obtained by the finite element method or other methods.

### 2.2 Estimation of Displacement of Virtual Plane Using the Regularization Method

#### 2.2.1 Regularization Method Developed From JE-MAP Method.

The JE-MAP method proposed by our research group is a regularization method that alternately repeats the estimation of latent and physical variables to improve the estimation accuracy. The regularization method in this study is an extension of the JE-MAP method. First, the latent variables indicate the areas of cracks and ligaments on the virtual plane. Next, the physical constraints on the displacement and force at the cracks and ligaments are added to the MAP method. The physical constraint has a cross–sparse relationship between the displacement and load at the cracks and ligaments. Furthermore, the GC method is improved to predict the boundary between them.

Figure 2 shows a schematic diagram of the regularization method developed in this study. Figure 2(a) shows a schematic diagram of the strain changes distribution of the observed surface. Figure 2(b) shows a schematic diagram of the virtual plane. Figure 2(b-1) shows a schematic diagram of the displacement distribution and (b-2) shows a schematic diagram of the latent variable distribution. Constant matrix **H** is a matrix to calculate the strain changes of the observed surface from the displacement of the virtual plane. The physical variable is the displacement of the virtual plane. The virtual plane's likelihood distribution of the displacement was obtained from strain distribution $\epsilon Z\u0303$ of the observed surface. A key feature of the JE-MAP method is the introduction of latent variable **z** to rationally incorporate prior information. Latent variables are indicated by 0 and 1 for the region between the cracks and the ligaments on the virtual plane. The prior distribution of the virtual plane's displacement is obtained from distribution **z** of the latent variable and the physical constraints at the cracks and ligaments. Posterior distribution **u** of the displacement of the virtual plane is obtained from the likelihood and prior distributions using the MAP method. Distribution **z** of the latent variable is estimated by binarizing posterior distribution **u** of the virtual plane's displacement obtained by the MAP method using the GC method. The prior distribution is updated with the updated latent variable distribution **z**. Displacement distribution **u** of the virtual plane is updated by the MAP method from the updated prior and likelihood distributions. In this way, displacement distribution **u** and latent variable distribution **z** are repeatedly and alternately identified by the MAP and GC methods.

#### 2.2.2 Estimation by MAP Method and Likelihood Distribution.

**u**are the most frequent value of the posterior distribution. In other words, crack surface displacement $u*$ maximizes the posterior probability shown in Eq. (2). The posterior distribution is obtained from the likelihood and prior distributions. $L(u|\epsilon Z\u0303)$ is the likelihood distribution, and $p(u|z)$ is the prior distribution

**Σ**is their covariance matrix. The measurement error is assumed to be independent of the measurement position. $\mu w$ is 0. The variance of measurement error $\sigma w2$ is assumed to be constant. The covariance matrix becomes the following equation:

_{w}#### 2.2.3 Prior Distribution.

The following physical constraints on the displacement and force at the cracks and ligaments are reflected in the prior distribution of the MAP method.

*Z*-direction of the virtual plane is discontinuous in the crack region and continuous in the ligament region. The force in the

*Z*-direction of the virtual plane is zero in the crack region and nonzero in the ligament region. Therefore, the force in the virtual plane's

*Z*-direction is a sparse distribution. If the model is symmetric in the

*Z*-direction at the virtual plane, the displacement in the

*Z*-direction is also sparse at the virtual plane. The

*Z*-forces and displacements of the virtual plane have conflicting regions of zero at the cracks and ligaments. The conflicting sparse constraints are called cross sparsity. A schematic diagram of the cross-sparse constraint condition is shown in Fig. 3, which shows a model symmetrical in the

*Z*-direction in the virtual plane. The forces and displacements in the

*Z*-direction of the virtual plane are cross sparse at the cracks and ligaments (Section A). Cross-sparse constraints are reflected in the prior distribution of the MAP method by latent variable

*z*. The size of the vector of latent variables

*z*is identical to

**u**. The component of

*z*is 1 for the crack-opening region and 0 for the ligament region. Prior distribution $p(u|z)$ incorporates the cross sparsity of the displacements and the forces in the virtual plane by

**z**and is given by the following equation:

**f**is the discretized

*Z*-direction force vector of the virtual plane. Force vector

**f**has a relationship with displacement

**u**shown in the following equation:

**G**is a constant matrix representing the relationship between

**f**and

**u**. $f0$ is a force vector when there is no crack. $\mu u$ is the mean vector of the displacements on the virtual plane,

**Σ**is the covariance matrix of the displacements on the virtual plane, $\mu f$ is a mean vector of the forces on the virtual plane, and

_{u}**Σ**is a covariance matrix of the forces on the virtual plane. The displacements and loads of the virtual plane are independent of the virtual plane's discretized positions. The covariance matrices become the following equations:

_{f}**Σ**, $\mu f$, and

_{u}**Σ**

_{f}by latent variable

*z*. $\mu u$,

**Σ**, $\mu f$, and

_{u}**Σ**are defined by Eqs. (9) and (10). Here the model is assumed to be symmetric in the

_{f}*Z*-direction in the virtual plane

#### 2.2.4 Updating Latent Variable z.

Latent variable *z*, which represents the crack region as 1 and the ligament region as 0, is updated from the distribution of **u** or **f** on the virtual plane estimated by the MAP method. Latent variable **z** is estimated using the GC method, which divides the image into foreground and background elements. To estimate the variable with the GC method, the distribution of **u** or **f** on the virtual plane is converted to an image. In this paper, latent variable **z** is estimated from the distribution of *u* on the virtual plane. The GC method uses the opencv 3.4.1 algorithm [14].

To convert displacement **u** into an image, the displacement distribution is normalized by the maximum value of **u**. If the normalized displacement distribution is converted directly into an image and segmented, the rapid displacement change near the boundary between the crack and the ligament will be ignored and the crack area will be underestimated. To magnify the abrupt displacement change, the entire normalized displacement distribution is multiplied by a constant greater than 1. In a constant multiplied distribution, all values greater than 1 are assumed to be 1. Next, the distribution multiplied by the constant is converted to an image. The initial value of latent variable **z** is given as 0 or 1 across the entire virtual plane to avoid biasing the estimation results.

#### 2.2.5 Calculation Flow of Regularization Method Developed From JE-MAP Method.

The flow of the calculation that repeatedly updates the prior distribution is shown in Fig. 4. Latent variable **z** is updated by the GC method from the displacement distribution estimated by the MAP method. In this paper, the update of latent variable **z** is repeated a specified number of times. The displacement distribution $u*$, where the change in **z** is smaller than a predetermined criterion, is extracted from the results of repeated calculations. The strain distribution on the observed surface is obtained by multiplying each extracted displacement distribution $u*$ by a constant matrix **H**.

where Num denotes the number of discretized strain distributions in the observed surface.

The initial values of the likelihood and prior distributions, $\sigma w2$, *μ*_{fl}, $\sigma fl2$, *μ*_{uc}, and $\sigma uc2$, and the variances at $\mu f=0$ and $\mu u=0$, are given in advance. The values in this paper are given in Sec. 3.3.

## 3 Test Specimen for Crack Estimation

### 3.1 Shape of Test Specimen.

A specific example of an inverse problem treated in this section aims at estimating a crack in the thickness direction from one side of a flat plate from the strain changes on the crack-free side. The strain changes on the crack-free surface are measured by applying a tensile load to the plate. Figure 5 shows a schematic diagram of a flat plate. It is 50 mm wide and 24 mm thick. The estimated crack exists on the virtual plane, which is indicated by the red dashed line at the plate's center (Fig. 5). The virtual plane's width direction is the *X*-axis, its thickness direction is the *Y*-axis, and the steel plate's length direction is the *Z*-axis. The tensile load is applied uniformly in the *Z*-axis direction.

### 3.2 Calculation of Constant Matrices **H** and **G**.

Constant matrix **H** in Eq. (1) is used in the observed equation of the inverse problem. Constant matrix **G** in Eq. (6) is used to incorporate the cross sparsity into the prior distribution. In this paper, constant matrices **H** and **G** are obtained by the finite element method. Figure 6 shows a schematic diagram of the discretized observed surface and the virtual plane. The vectors of strain $\epsilon Z$, displacement **u**, and force **f** are also shown in Fig. 6. The finite element model is a 1/2 symmetric model of a flat plate where a virtual plane is a symmetry plane. In this paper, to measure the strain changes with the DIC method, a load of 1000 *μ* generated was applied to the observed surface of the flat plate in Fig. 5 without cracks. A load of 247 kN in the *Z*-direction was applied to the opposite side of the symmetry plane in the *Z*-direction. The discretized strains, displacements, and forces are the values at the nodes of the elements created in a grid by the finite element method. The virtual plane is divided into n elements in the *X*-direction and m elements in the *Y*-direction. The observed surface is divided into n elements in the *X*-direction and *p* elements in the *Z*-direction.

*Z*-direction strain of the observed surface. The variation in the

*Z*-direction strain is the difference in the strain between the conditions without/with cracks on the virtual plane. $(i,j)$ denotes the location of the nodes that are deemed to be cracked.

*i*is a value between 1 and

*n*, and

*j*is a value between 1 and

*m. u*(

*i*,

*j*) is a difference vector of the

*Z*-direction displacement of the virtual plane. The difference in the displacement in the

*Z*-direction is the difference in the displacement between the conditions without cracks and with them in the virtual plane.

*f*(

*i*,

*j*) is the force vector in the

*Z*-direction of the virtual plane under the crack condition. $\epsilon z(i,j)$,

*u*(

*i*,

*j*), and

*f*(

*i*,

*i*) are calculated for the number of nodes in the virtual plane.

*f*

_{0}is a force vector in the

*Z*-direction of the virtual plane under a crack-free condition. Constant matrices

**H**and

**G**are obtained from $\epsilon z(i,j)$,

*u*(

*i*,

*j*),

*f*(

*i*,

*i*), and

*f*

_{0}by Eqs. (12) and (13). Here $\beta $ in Eq. (1) is set to zero. On the virtual plane, the 50 mm wide are divided into 25 subdivisions and the 24 mm thick are divided into 12 subdivisions. On the observed surface, the 50 mm wide are divided into 25 subdivisions and the 24 mm thick are divided into 15 subdivisions. The material of the flat plate is S45C steel (C45 steel:ISO) with a proof stress of 490 MPa or higher. This S45C steel was made by Kobe Steel, LTD. in Japan. Young's modulus and Poisson's ratio of the plate are set to 206 GPa and 0.3. The finite element method solver is ansys 19.2 [15]

### 3.3 Parameters of Regularization Method Developed From Joint Estimation Maximum a Posteriori Method.

In this section, the parameters in Secs. 2.2.4 and 2.2.5 were determined. First, when updating the latent variables in Sec. 2.2.4, the constants that multiply the normalized displacement distribution are described. The constant to be multiplied denote is denoted *z*Const. *z*Const ranges from 3 to 6, which improves the crack-estimation accuracy of the numerical experiments using the finite element method.

*μ*

_{fl}, $\sigma fl2$,

*μ*

_{uc}, $\sigma uc2$, and the variance at $\mu u=0$, and the variance at $\mu f=0$. The variances at $\sigma fl2,\sigma uc2,\mu u=0$, and $\mu f=0$ are relative values to

*μ*

_{uc}and

*μ*

_{fl}, as shown in the following equations:

The value of *μ*_{fl} is the load applied to the flat plate divided by the number of divisions of the virtual plane. A range is set for *z*Const, $\sigma w2$, and *μ*_{uc} from which the optimum value is selected. Standard deviation $\sigma w2$ corresponds to the strain distribution's measurement error, whose strain distribution is in the range of 1 *μ*–100 *μ*. The initial value of *μ*_{uc} is the crack aperture in the range of 0.1 *μ* m–10 *μ* m. Bayesian optimization is used as the optimization method [16].

## 4 Strain of Observed Surface Measured by Digital Image Correlation Method

### 4.1 Measurement by Digital Image Correlation Method.

The strain changes due to the presence or absence of a crack is determined from a specimen's strains without/with a crack. Figure 7 shows the geometry of the specimens for which the strain was measured. Figure 7(a) shows a specimen without a crack; Fig. 7(b) shows a specimen with a crack. The geometry in Fig. 7 includes the part gripped by the testing machine to apply the load. The specimens are made of C45 quenched and tempered material. The specimen with a crack has a semi-elliptical slit that is 0.25 mm thick, 24 mm wide, and 9.8 mm depth. A semi-elliptical, longitudinal slit was made in the center of the specimen by electrical discharge machining. Its depth *c*-direction is in line with the specimen's thickness direction.

The strain-measuring surface of the specimen was given a random black-and-white pattern to measure the strain by the DIC method. The size of the black pattern is also random, with a maximum size of about 1 mm. The black-and-white random pattern was created using water-based acrylic paint. The surface where the strain is measured corresponds to the observed surface of the inverse problem. Figure 8 shows the situation due to the tensile load where the strain is measured by the DIC method. A 247 kN load was applied in the *Z*-direction (Fig. 8) by a servohydraulic fatigue-testing machine whose capacity is 500 kN. At a tensile load of 247 kN, the plastic zone around the semi-elliptical slit is small and satisfies the small scale yielding hypothesis. The images used in the DIC method were taken with a CCD camera (NICOND7200) that had a zoom lens with a focal length of 18–200 mm (AF-S DX NIKKOR18-200 mm f/3.5-5.6GEDVRII). The resolution of the images was set to 6000 × 4000 pixels. The images for the strain determination were taken at loads of 0 and 247 kN for all the tests. The strains were calculated from the correlation information between the 0 and 247 kN images and analyzed using vic-2d software for DIC from correlated solutions. The vic-2d was configured with a subset size of 151 pixels, steps of 25 pixels, and a filter of 51 pixels for strain determination.

Figure 9 shows the strain distribution in the *Z*-direction of the observed surface calculated by the DIC method at load of 247 kN. Figure 9(a) shows a specimen's strain distribution without a crack, and Fig. 9(b) shows it with a crack. The no-crack specimen has a strain distribution without any characteristic distribution. The strain distribution had a mean value of 1005 *μ* and a standard deviation of 26 *μ*. The specimens with cracks showed a decrease in strain near the strain distribution's center. The decrease in strain was caused by a crack on the opposite side of the observed surface.

### 4.2 Strain Changes as Observed Data for Inverse Analysis.

The observed data for the inverse problem is the strain changes obtained by the DIC method from images taken before and after a crack occurred in an electrical device. In this paper, the strain changes due to the crack are simulated by subtracting the strain with a crack from the strain of the specimen without a crack. The observed data, extracted from the strain distribution obtained in Sec. 4.1, denotes the amount of strain changes in the *Z*-direction with/without cracks. Figure 10 shows the difference between the strains without/with a crack in the *Z*-direction. The strain distribution in Fig. 10 is 0 in the *Z*-direction at the lower left and 0 in the *X*-direction at the left end of the specimen width. The observed data are calculated from the strain distribution in Fig. 10. The change in the strain distribution due to the crack is located in the entire *X*-direction around 32 mm of the *Z*-coordinate in Fig. 10. The observed data are the discrete values of the strain distribution shown in Fig. 10 divided into the element sizes shown in Sec. 3.2.

Figure 11 shows the method used to determine the location of the crack in the *Z*-direction based on the strain changes are Fig. 10. Figure 11(a) shows the method used to determine the *Z*-coordinate where the strain changes is maximum. The discretized strain change was divided into each *X*-coordinate. The maximum *Z*-coordinate is obtained for each divided strain change. Figure 11(b) shows the frequency distribution of the *Z*-coordinates with a maximum strain change. The virtual plane is located at the *Z*-coordinate with the highest frequency of the maximum strain change. The strain change's distribution is divided by the virtual plane into positive and negative directions of the *Z*-axis. The observed data are the amount of strain changes in the positive direction of the *Z*-axis. Figure 12 shows the strain distribution used as the observed data in the inverse analysis. The strain distribution is missing at the edge of the specimen in the *X*-direction. This missing distribution occurred because the strain could not be calculated due to the small number of patterns included in the subset of the DIC method. The range of the observed data in the *Z*-direction was identical to the observed surface shown in Sec. 3.2. The range of the observed data in the *X*-direction is the range measured by the DIC method. Constant matrix **H** is recalculated based on the range of the observed data and used in the inverse analysis in Sec. 5.

## 5 Estimated Results

### 5.1 Inverse Analysis by L1-Norm Regularization.

In this section, the inverse problem, in which the strain changes in Fig. 12 is the observed data, is solved by a generic method. The generic method is L1-norm regularization. Regularization parameter *λ* was $4.33\xd710\u22127$ by applying “the one standard error rule” [17] to the *k*-fold cross-validation results [18] with *k *=* *5. The generic method is shown in the Appendix. Figure 13 shows the displacement distribution in the *Z*-direction of the virtual plane analyzed by L1-norm regularization. In Fig. 13, the surface where the specimen crack opens is 0 in the *Y*-direction and the left end of the specimen width is 0 in the *X*-direction. The displacement in the *Z*-direction of the virtual plane is positive for the crack and zero for the no-crack part. Figure 14 shows the crack extracted from Fig. 13. The crack area is defined as the area of displacement greater than 10% of the maximum displacement of the virtual plane and is indicated by 1. Displacement other than that of the crack is indicated by 0. The zero points of x and y in Fig. 14 are identical to in Fig. 13.

Figure 15 shows the displacement distribution of the ground truth, which was obtained by applying the same load as in the test using the finite element method to a model with the same crack as the test specimen. The zero points of *X* and *Y* in Fig. 15 are identical to in Fig. 13.

The displacement distribution in Fig. 13 reflects the sparsity of the virtual plane since the displacement is estimated separately for the cracks and the ligaments. However, unlike Fig. 15, the estimated displacements are located away from the virtual plane's free edge. Furthermore, the maximum displacement in the *Z*-direction differs significantly from that in Fig. 15. The cracks extracted from the displacements shown in Fig. 14 also differ significantly from the cracks in Fig. 15. There are multiple cracks in the virtual plane when they are extracted from the displacement shown in Fig. 13. The generic method combines constant matrix **H** obtained by the finite element method and a solution method that introduces sparsity in the *Z*-direction displacement. It did not take into account the measurement error included in the strain changes of the observed data, resulting in a result different from the ground truth. The JE-MAP method, shown in the next section, is an inverse analysis method that takes into account the measurement error included in the strain changes of the observed data.

### 5.2 Inverse Analysis Using Regularization Method Developed From Joint Estimation Maximum a Posteriori Method.

In this section, the inverse problem with the strain changes in Fig. 12 as the observed data is solved by the JE-MAP method. GC used the opencv functions. To employ them in matlab, mexopencv [19] was used. The initial value of latent variable distribution *z* was set to 1 for the entire virtual plane. Table 1 shows the initial values of *z*Const, $\sigma w2$, and *μ*_{uc} used in the regularization method developed from the JE-MAP method. The initial values in Table 1 were obtained by Bayesian optimization, which is shown in the Appendix. The standard deviation of measurement error *σ _{w}* is close to the standard deviation of the strain distribution measured by DIC, $2.6\xd710\u22125$ (Sec. 4.1). The convergence criterion for latent variable

**z**was less than two elements of the latent variable whose values changed with the update. The MSE is calculated from the distribution of the displacements on the virtual plane where the latent variables converged. The estimated result is the displacement distribution of the virtual plane where the MSE was minimized.

σ_{w} | The initial value of μ_{uc} (mm) | zConst |
---|---|---|

$1.40\xd710\u22125$ | $7.27\xd710\u22124$ | 5.01 |

σ_{w} | The initial value of μ_{uc} (mm) | zConst |
---|---|---|

$1.40\xd710\u22125$ | $7.27\xd710\u22124$ | 5.01 |

Figure 16 shows the distribution of the displacement in the *Z*-direction of the virtual plane estimated by the regularization method developed from the JE-MAP method. Figure 17 shows the latent variables that were input to the prior distribution calculated in Fig. 16. The zero points of *x* and *y* in Figs. 16 and 17 are identical to those in Fig. 13. The maximum point of the displacement distribution in Fig. 16 is near the ground truth in Fig. 15, at the free edge of the virtual plane. The range of the latent variable crack in Fig. 17 is one in the virtual plane, similar to the ground truth in Fig. 15. The range of the cracks in Fig. 17 falls mostly within the range of the ground truth in Fig. 15. The regularization method developed from the JE-MAP method can estimate a crack's location.

The area of the largest crack estimated by the regularization method developed from the JE-MAP method is compared with the area estimated by the L1-norm regularization. The area was determined from the element number of cracked regions shown in Figs. 17 and 14. One element in the distribution figures in Figs. 13–17 has a width of 2 mm and a height of 2 mm. The area estimated by the regularization method developed from the JE-MAP method is 120 mm^{2}. The area estimated by the L1-norm regularization is 40 mm^{2}. The area estimated by Fig. 15 regularization is 224 mm^{2}. The area estimated by the L1-norm regularization is 18% of the ground truth. The area estimated by the regularization method developed from the JE-MAP method is 54% of the ground truth and closer to the ground truth than the L1-norm regularization. The results estimated by the regularization method developed from the JE-MAP method are closer to the ground truth in terms of the number of cracks, crack locations, and crack's area than the results estimated by the L1-norm regularization.

The estimated size of the crack is compared to the ground truth. The elements of the crack shown in Fig. 17 were compared to the elements with displacements in Fig. 15. The sizes estimated by the regularization method developed from the JE-MAP method are a maximum crack width of 22 mm and a maximum crack depth of 6 mm. The maximum width and depth could be estimated by 77% and 50%, respectively, for the ground truth in Fig. 15. Therefore, the regularization method developed from the JE-MAP method can estimate a crack's size and location from the observed data calculated by the DIC method.

We believe that the following factors contributed to the estimation of the regularization method developed from the JE-MAP method. First, the prior distribution of the MAP method considered the cross sparsity of displacement and force on the cracked surface by latent variables. Second, the likelihood distribution of the MAP method took into account the measurement error of the observed data. Next, the latent variables were automatically updated by the GC method from the estimation results of the MAP method. The displacement of the virtual plane was updated again with automatically updated latent variables. Finally, from the multiple displacements of virtual plane obtained by updating the latent variables, the displacements at which the latent variables converged are extracted. The estimation result is the displacement of the virtual plane that minimizes the MSE between the strain calculated from the extracted displacement of the virtual plane and constant matrix **H** and the measured strain.

The displacement distribution in Fig. 16 shows that the displacement in the crack-free area is about 20% of the maximum displacement in the crack. It is possible that the displacement occurred in the crack-free area due to the difference between the Gaussian distribution assumed in the likelihood distribution and the error distribution in the actual measurement. In the future, we will improve the displacement's estimation accuracy in the crack-free area.

## 6 Conclusions

This paper studied a method to estimate cracks in invisible locations by combining a regularization method developed from the JE-MAP and a DIC method. The results of this paper are as follows:

In the inverse problem in this study, the strain on the crack-free surface is set as the observed data and the displacement on the crack-containing plane as the unknown data.

The regularization method, which is a development of the JE-MAP method, introduces physical constraints on the displacements and the forces at the cracks and ligaments into the MAP method. The physical constraints are a cross-sparse relationship between the displacement and the load at the cracks and ligaments.

The regularization method, which is an extension of the JE-MAP method, improves the prediction of the crack and ligament boundaries by the grab-cut method. The displacement or load distribution is converted into an image that highlights the crack and ligament boundaries before being input into the grab-cut method.

The regularization method, which is an extension of the JE-MAP method, estimated the displacement on the crack-containing plane more accurately than L1-norm regularization in the inverse problem where the observed data were strain distributions measured by the DIC method.

The cracks obtained from the displacement on the crack-containing plane estimated by the regularization method, which is an extension of the JE-MAP method, are more accurate than those predicted by the L1-norm regularization.

A combination of the DIC method with the regularization method, which is an advanced version of the JE-MAP method, can estimate the size and location of cracks that are not visible from camera images.

In the future, the validity and reproducibility of the method proposed in this paper will be verified on observed data measured by DIC method for strains on surfaces of various geometries. In addition, the validity of the JE-MAP method will be verified for cracks of various geometries, e.g., multiple cracks.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

- DIC =
digital image correlation

**f**=discretized

*Z*-direction force vector of the virtual plane- $f0$ =
force vector when there is no crack

**G**=constant matrix representing the relationship between

**f**and**u**- GC =
grab-cut

**H**=constant matrix representing the relationship between $\epsilon Z$ and

**u**- MAP =
maximum a posteriori

- MSE =
mean-square-error

*u*=displacement of the virtual plane in the

*Z*-direction**u**=displacement vector of the virtual plane in the

*Z*-direction- $u*$ =
crack surface displacement maximizing the posterior probability shown in Eq. (2)

*z*Const =constants to be multiplied by the normalized displacement distribution when updating the latent variables

- $\beta $ =
discretized vector of the measurement errors

- $\epsilon Z$ =
strain changes in the

*Z*-direction of observed surface- $\epsilon Z\u0303$ =
vector of the measured strain changes in the

*Z*-direction of observed surface*λ*=regularization parameter of

*L*1-norm regularization- $\mu f$ =
mean vector of the forces on the virtual plane

- $\mu u$ =
mean vector of the displacements on the virtual plane

- $\mu w$ =
mean vector of the measurement errors

- $\sigma w2$ =
variance of the measurement error

- Σ
=_{f} covariance matrix of the forces on the virtual plane

- Σ
=_{u} covariance matrix of the displacements on the virtual plane

### L1-Norm Regularization and Determination of Regularization Parameter *λ* by Cross-Validation

**u**in the virtual plane is sparsely distributed since it is greater than zero at the crack and zero at the noncrack locations. L1-norm regularization [13], which is a method for estimating sparse unknowns, is a generic method in this paper. The solution is a displacement

**u**that satisfies the following equation:

*λ* denotes a regularization parameter. Regularization parameter *λ* was calculated by *k*-fold cross-validation with *k* = 5. Appendix Eq. (A1) was solved using the fast iterative shrinkage-thresholding algorithm method [20]. The method for determining regularization parameter *λ* is shown. The number of data partitions for cross-validation was set to 5. The observed data were partitioned by cvpartition, a data creation function for cross-validation [21] in matlab.

*λ* was searched for in the range of $1\xd710\u22128$ to $1\xd710\u22125$. One hundred *λ* were extracted from the search range at logarithmically evenly spaced. Cross-validation was performed for all the extracted *λ*. Appendix Fig. 18 shows the mean and standard deviation of MSE calculated by cross-validation. The vertical axis shows the MSE, and the horizontal axis shows *λ*. In Appendix Fig. 18, the mean and standard deviation of the MSE calculated by cross-validation for each value of *λ* are shown as circles and error bars. *λ* was determined to be $4.33\xd710\u22127$ by applying “the one standard error rule” to the results in Appendix Fig. 18.

### Determination of Parameters for Joint Estimation Maximum a Posteriori Method by Bayesian Optimization

*σ*,

_{w}*z*Const, and

*μ*

_{uc}of the JE-MAP method were determined by Bayesian optimization [22] of matlab function bayesopt. Appendix Eq. (B1) shows the objective function of the Bayesian optimization. $uJE\u2212MAP*$ in the objective function is the displacement vector in the

*Z*-direction of the virtual plane calculated by the JE-MAP method with input

*σ*,

_{w}*z*Const, and

*μ*

_{uc}

The displacement vector in the *Z*-direction of the virtual plane was calculated by the JE-MAP method, and *σ _{w}*,

*z*Const, and

*μ*

_{uc}were optimized in the following ranges. The following are the arguments of bayesopt that were changed from the default values.

Acquisition function: use “expected-improvement.” Enable modification of the acquisition function to escape a local objective function minimum. “Acquisition Function Name” = “expected-improvement-plus”.

Specify deterministic objective function: true (the objective function is specified deterministically).

Objective function evaluation limit: 60.

## References

**21**(6), pp.