Soft Tissue Deformation Modeling in the Procedure of Needle Insertion: A Kriging-Based Method

The simulation and planning system (SPS) requires accurate and real-time feedback regarding the deformation of soft tissues during the needle insertion procedure. Traditional mechanical-based models such as the finite element method (FEM) are widely used to compute the deformations of soft tissue. However, it is difficult for the FEM or other methods to find a balance between an acceptable image fidelity and real-time deformation feedback due to their complex material properties, geometries and interaction mechanisms. In this paper, a Kriging-based method is applied to model the soft tissue deformation to strike a balance between the accuracy and efficiency of deformation feedback. Four combinations of regression and correlation functions are compared regarding their ability to predict the maximum deformations of ten characteristic markers at a fixed insertion depth. The results suggest that a first order regression function with Gaussian correlation functions can best fit the results of the ground truth. The functional response of the Kriging-based method is utilized to model the dynamic deformations of markers at a series of needle insertion depths. The feasibility of the method is verified by investigating the adaptation to step variations. Compared with the ground truth of the finite element (FE) results, the maximum residual is less than 0.92 mm in the Y direction and 0.31 mm in the X direction. The results suggest that the Kriging metamodel provides real-time deformation feedback for a target and an obstacle to a SPS.


Introduction
In needle-based minimally invasive surgeries, needle insertion errors always lead to a large probability of serious complications and increase the physical pain of patients [1][2][3][4]. Real-time feedback regarding soft tissue deformation during needle insertion, especially in the region around the target and obstacles, and the experience of the attending physicians affect the insertion accuracy [4]. The simulation and planning system (SPS) can provide pre-and intraoperative deformation information and help physicians to sharpen their skills [5], which can effectively improve the success rate of minimally invasive surgeries. Soft tissue visualization during the needle insertion procedure, which is the most important part of the SPS, is realized by a deformation model of the soft tissue [6].
The soft tissue deformation model has been widely studied in Refs. [7][8][9][10][11], of which continuum mechanical-based modeling methods such as the finite element method (FEM) are the most popular [11][12][13]. Dimaio et al. [13] established the soft tissue of a simple geometry using the FEM, in which the refined elements along the needle path were configured and the interaction force was experimentally loaded onto the node path. The material properties of the tissue were modeled as the linear elasticity. Misra et al. [14] reviewed and sorted out the tool-tissue interaction models based on the tissue material properties, which revealed that the hyperelastic  35:112 model had a relatively high accuracy. Kobayashi et al. [15] modeled the nonlinear and viscoelastic liver based on the finite element method. The complicated organ geometry and boundary conditions were taken into account in Refs. [16,17]. The complex boundary conditions and geometry of real organs, as well as the viscoelastic and heterogeneous material properties, increased the computational complexity. Furthermore, the needle insertion process includes material failure and a new internal surface formation; thus, the underlying physics is also an essential part. Some of the traditional FEM modeling methods, although able to provide fast calculations of tissue deformation, lead to unpleasing visual results during the needle insertion procedure due to their rough simplification of the needle puncture process. Takabi [18] presented a comprehensive survey mainly on the modeling of tissue cutting, focusing on the mechanical concepts and modeling techniques. Gokgol et al. [19] modeled the FEM cutting process by defining a fracture work, the nodal in the FEM model break when the cutting energy exceeded the fracture work. The thresholds can also be other failure parameters such as the maximum shear strain [20]. Oldfield and Dini used cohesive elements in the ABAQUS commercial software to describe the material failure [21]. However, all of these methods require the elements along the needle path to be sufficiently small, which may result in a large computational burden and instability in the calculations. Recently, the extended finite element method (X-FEM), typically utilized for crack propagation, was also applied to soft tissue cutting [22]. It reserved the initial mesh topology without very high-resolution grids in modeling the discontinuities, but an appropriate enrichment function was difficult to establish. To eliminate the mesh dependence in the traditional FEM, the meshless method was proposed for soft tissue deformation modeling, which is particularly suitable for simulating large deformation and cutting processes since it is not required to calculate equations on the grid [23,24]. However, the basis interpolation function in the meshless method is complex and difficult to construct.
A successful surgery simulation and planning system requires two contradictory basic processes: physical reality and real-time dynamic interaction. However, most existing needle insertion deformation models with high accuracy during needle insertion are simulated offline. It is difficult for the FEM or other methods to balance the acceptable image fidelity and real-time deformation feedback, which limits their applications in real clinical trainings [25,26]. In this paper, a Kriging metamodel is applied to model and predict the deformation of obstacles and targets inside soft tissue during the needle insertion procedure. The Kriging metamodel was addressed by the South African mining engineer Krige and developed by Matheron, a French mathematician [27]. It is a simplified computer-based simulation model that has an input/output function based on the surface response [28]. The innovation of this paper is to propose a needle insertion simulation model that can balance the computational burden and accuracy of needle insertion simulation, according to the soft tissue characteristics of needle insertion process. The advantages of this paper are twofold: First, the Kriging metamodel can reduce the large requirement for computing resources and offer real-time simulation, especially being able to overcome the challenges of computational techniques in complex simulation models [29]. Hence, depending on the accurate simulation dataset generated by other offline needle-tissue interaction models, the online Kriging model can be established rapidly without reducing the accuracy of the original model by much. Second, the computer-based simulation program is deterministic, that is, the same input corresponds to the same output, which cannot fully reflect the uncertainty of mechatronic systems [30,31]. In robotassisted needle insertion systems, complex factors such as the imaging equipment, physician's skill and patient's condition lead to uncertainty regarding tissue deformation. The Kriging metamodel with the random sampling method can reflect the uncertainties of robot-assisted needle-tissue interactions.
In our work, the input dataset is generated by Latin hypercube sampling (LHS), the parameters of which include the material properties of the tissues and needle, the geometrical properties of the needle and the solver parameters. Ten characteristic markers are selected to represent the targets and obstacles inside the tissue body. The output dataset is generated by running the needle-tissue coupling FE model presented in Ref. [32], which is used as the ground truth of the Kriging-based model. The validity and feasibility of the proposed Kriging-based model are analyzed, and the results suggest that the combination of the Kriging metamodel and the high-precision finite element model provides real-time deformation feedback for a target and an obstacle to the SPS.
The rest of this paper is organized as follows. Section 1 introduces the data preparation, the basic theory of the Kriging model and the maximum deformation modeling of needle insertion at a fixed depth. In Section 2, the functional-based Kriging model is applied to establish the deformation model at a series of depths, of which the adaptation to the step variation is also analyzed. Section 3 discusses the proposed Kriging model. Finally, conclusions and future work are given in Section 4.

Flowchart of Kriging Deformation Online Prediction
A flowchart of Kriging online deformation prediction is shown in Figure 1. In offline process, a physic-based model, such as finite element based needle insertion model is utilized to generate dataset. The dataset includes the input parameters and output deformations in different locations. The input parameters include material parameters of soft tissues, needle parameters and other parameters affecting the deformation results obviously. The dataset is provided and offline trained for Kriging meta model. Training process is for tuning Kriging parameters to best fit the deformation prediction model. Kriging model is in essence a computational model; thus, it can realize real-time prediction while maintain the accuracy of the physic-based model. In surgery simulation and planning, when input current parameters of tissue and needle, the Kriging model can online predict soft tissue deformation for surgery training or needle path planning.

Deformation Modeling of Needle Insertion at a Fixed Depth
A finite element model for the needle-tissue interaction was built by our group to simulate tissue deformation [32], in which a modified local constraints method was applied. In this paper, we aim to predict the motion of the target and obstacles inside the tissue in real time. For validation convenient, ten characteristic markers in the FE model are selected as the observation points that distributed in the whole computational domain to represent the deformation information of different parts, as shown in Figure 2. The observation points are chosen according to the distance from the needle body. In this paper, 10 characteristic markers are selected, which can be classi-    Table 1, wherein x 1 , x 2 , and x 11 are the material parameters of soft tissues, x 3 -x 7 are the needle parameters, x 8 -x 9 are the FEM solver parameters and x 10 is the friction coefficient between the needle and soft tissue. The boundary values are chosen from the empirical value and simulation results in Ref. [32].
The input variables are denoted as x = [x 1 , · · · , x 11 ] , which are employed to construct the design matrix of the Kriging-based model. Running the finite element model 20 times by using 20 groups of different combinations of the input variables x = [x 1 , · · · , x 11 ], the displacements at the 10 observed locations are collected as the output responses of the Kriging model, denoted by The soft tissue deformation occurs in both the x and y directions in the two-dimensional case; hence, the output response u xy is the resultant displacement in both directions and is written as where u x and u y are the displacements in the x and y directions, respectively. Since the resultant displacement increases as the needle insertion depth into the tissue increases, the maximum value of u xy is chosen as the output response of the Kriging-based prediction model.

Kriging Model Construction
The Kriging prediction model is an interpolation of known observation locations, of which the mean square error equals zero. The Kriging model is usually expressed as a combination of a polynomial and its deviation [27], written as where f (x) is the basis polynomial function, usually a regression function.
where F is the design matrix, written as: m is the number of groups of FE model, and in this application, m = 20. p is the number of basis functions, which is determined by the type of basis function. In our experiments, p = 12 for first order regression functions. The variable x i is obtained from the FE model by the Latin hypercube sampling method, and x i is normalized to the interval [0, 1] , as in where x i is the normalized data of x i and x L and x U are the minimum and maximum values of the variable x i , respectively.
The mean square error of the Kriging prediction model is formulated as [27] where is the maximum likelihood estimate of the variance.

Regression Function
The design matrix F m×p in the Kriging prediction model is where the basis function f (x) has several forms. Usually, constant and linear forms of the regression function are used to construct the basis function [33]. • Constant form, p = 1 , of which the design matrix F is a column vector: • Linear form, p = n + 1, of which the basis function is written as Eq. (8): where the component form of the estimate point x is expanded as x = [w 1 , w 2 , · · · , w n ] ∈ R n . In our experiment, both zero order (constant form, p = 1 ) and first order (linear form, p = 12 ) functions are adopted and compared.

Correlation Function
The Kriging model assumes that the correlation of the output is determined by the distance between the input variables. The correlation function of the input variables is written as the product of n one-dimensional correlation equations, as shown in Eq. (9): where d j = w j − x j and θ is the correlation coefficient. In the Kriging model, the most widely used correlation model is that shown in Eq. (10): where � · � denotes the Euclidean distance of d j . When p = 1 and p = 2 , the exponential and Gaussian correlation functions are as shown in Eq. (11) and Eq. (12): The correlation function decreases with the increase in the Euclidean distance d j , and a larger correlation coefficient θ leads to a rapid decline of the correlation function. Substituting the correlation function into Eq. (5), the mean square error is a function of the covariance σ 2 and the correlation coefficient θ . The optimal solution to the correlation coefficient θ * is converted into an unconstrained global optimization problem [34], as shown in Eq. (13): where |R| is the determinant value of the correlation matrix R.

Deformation Modeling of Needle Insertion at a Fixed Depth
The input and output data are prepared by running the finite element program 20 times, as in Ref. [32]. To establish a Kriging model, 19 sets of observation points [X in (1, :), · · · , X in (i − 1, :), X in (i + 1, :), · · · , X in (20, :)] are chosen as the input design sites. The colon in X in stands for the input observation points where there are 11 input variables in each observation point; hence, X in is a 19 × 11 matrix. The remaining set X in (i, :) is used as the test set. The zero order and first order regression functions and the Gaussian and exponential correlation functions are employed to construct the Kriging model for static soft tissue deformation. In this section, the maximum displacements of the marked points are regarded as the static soft tissue deformations. Since there are 10 markers in the Kriging model, the out-put of Kriging would be a 10 × 11 matrix. Given the input and output data of the model, the Kriging model would be built by two functions, namely the regression function and the correlation function to satisfy that the mean square error equals zero. Therefore, there are four kinds of Kriging model, namely, zero order Gaussian, zero order exponential, first order Gaussian and first order exponential, which are shown in Figure 3. In Figure 3, the experimental values represent the finite element simulation results. The abscissa indicates the number of observation locations, and the ordinate is the value of the maximum resultant displacement. The test set of each run is randomly selected. The upper-left test set is X in (5, :) , the upper-right test set is X in (6, :) , the lower-left test set is X in (7, :) , and the lower-right test set is X in (8, :) . From Figure 3, it is also suggested that the closer the tissue is to the needle body, the greater the deformation. To analyze the impacts of the regression function and correlation function on the model, the residual e of the simulation is defined as the difference between the finite element simulation value y and the Kriging prediction value ŷ , that is, e =ŷ − y . The maximum residual values of the zero order Gaussian, zero order exponential, first order Gaussian and first order exponential Kriging models are 7.5108 mm, 7.3327 mm, 3.2752 mm and 3.2994 mm, respectively. Furthermore, the mean residual of each output response is defined as in Eq.  Figure 4. In Figure 4, the mean residuals of the zero order Gaussian and exponential types are 1.5650 mm and 1.3992 mm, respectively. The mean maximum residuals of the first order Gaussian and exponential are 0.8941 mm and 0.9177 mm, respectively. From the prediction results, the first order regression function is better than the zero order regression function, and the Gaussian correlation function is slightly better than the exponential function from the enlarged part.

Sensitivity Analysis of the Parameters
The Kriging-based prediction model of soft tissue deformation takes into account the 11 input parameters in the finite element simulation process, and the effect of the parameters X = [x 1 , x 2 , · · · , x 11 ] on the tissue deformation Y = [y 1 , y 2 , · · · , y 10 ] is investigated with the first order Gaussian Kriging-based prediction model in this section.
In this section, the 'local' sensitivity for each variable is investigated, i.e., only one column is varied for each row selected, while the other columns are fixed. Since the 20 groups sample is randomly generated, there is no specific law for the order of rows. Hence, we construct the parameter test matrix by randomly selecting the four rows of data X(1, :) , X(7, :) , X(13, :) and X(18, :) . The varied variable is evenly distributed in the interval [0, 1] (the interval is divided into 100 equal parts), while the other input variables retain the same values. The sample size is 4×100×11, which are enough to reflect the sensitivity.
The sensitivity of parameters is reflected by the average sensitivity index (SI), which is the average change in the output response and input variable of the four selected groups of samples, that is, SI = 1 4 4 k=1 �u kj xy /�x ki (j = 1, · · · , 10, i = 1, · · · , 11) . When the output displacement increases with the variables, the SI is greater than zero. In the opposite case, the SI is less than zero. According to the range of the SI, the 11 input variables can be classified into three types, i.e., leading parameters, nonleading parameters and disturbed parameters, as shown in Figure 5, where the abscissa represents 11 variables and the ordinate represents the SI.
Since there are ten observation points inside the tissue, there are also ten sensitivity indices for each variable; thus, we choose a box plot to show the extremum (black line at both ends), the interquartile range (blue box) and the median (red line) of one group of the dataset. We classify the 11 variables into three categories based on the sensitivity index: The leading parameters include Young's modulus of the soft tissue x 1 , the soft-tissue Poisson ratio x 2 , the needle length x 6 and the coefficient of friction x 10 , of which the absolute value of the median of the SI is larger than 0.5  35:112 and the effect on different markers is consistent (indices are either all distributed above or below zero). For example, the SI of x 1 is always below zero, which means that the soft tissue deformation decreases as x 1 increases. These variables have an obvious impact on the deformation due to the large absolute value of the SI. By contrast, the deformation increases as x 2 increases, with x 2 having a relatively small impact on the deformation. The disturbed parameters include the angle of the needle tip x 3 , insertion angle x 4 , Young's modulus of the needle x 7 and area density x 11 . The effect of the disturbed parameters varies for different markers. For the disturbed parameters, the sensitivity index covers the yellow line, which means that the deformations of some markers increase with the parameters, while some markers decrease.
The nonleading parameters include the needle diameter radius x 5 and the Rayleigh damping coefficients α ( x 11 ) and β ( x 9 ), which means that the absolute value of the median sensitivity is less than 0.02. The nonleading parameters affect the output displacement slightly and can be regarded as constants.
Sensitivity analysis of the parameters can offer a thorough understanding of the impacts of the input parameters on soft tissue deformation. To ensure rigor of results, other rows of data are also tried in the same way, and the tendencies and slopes came to the same conclusion. The leading and disturbed input parameters play key roles in predicting tissue deformation, directly affecting the overall trend and accuracy of the prediction results. The nonleading parameters have little effect (SI<0.02) on the prediction results and can be removed from the input parameters.
The sensitivity analysis provides a basis for reducing the number of input parameters of the Kriging-based prediction model. For a more complex FE model with more input parameters, such as that for a nonlinear, viscoelastic tissue material, the sensitivity analysis of parameters is beneficial for selecting the optimal input parameters and reducing the computation time.

Kriging Model of the Functional Response
According to the continuum mechanics, the motion equation of a point inside the tissue O in the Euclidean space is expressed as Eq. (15): indices. When the Kriging metamodel is used to solve the functional response problem, the calculation of the correlation coefficient of the maximum likelihood estimation includes the inverse matrix and the determinant of the correlation matrix, which leads to the instability and time complexity of the computation [35]. The Kronecker product is employed to construct the correlation matrix and to reduce the computational burden of the Kriging method [36].
The input variables in the FE simulation experiment are x = [x 1 , x 2 , · · · , x n ] T , x i ∈ R n , i = 1, 2, · · · , m , and the output response is a function of time t . The output response y i = [y i 1 , y i 2 , · · · , y i r i ] T is recorded at each time where y(x, t) is the output response of the input variable x at time t. F(x, t) = [f 1 (x, t), · · · , f p (x, t)] T is a set of known polynomial basis functions, where usually f 1 (x, t) = 1 . β is an unknown basis function coefficient. z(x, t) is a zero-mean Gaussian random function, the covariance function of which is shown in Eq. (17): Assume that the correlation function R(x 1 − x 2 , t 1 − t 2 ) is the product of n one-dimensional correlation equations, as shown in Eq. (18): where R i (x i1 − x i2 ) is the correlation function of the ith set of variables and R T (t 1 − t 2 ) is the correlation function at time t.
The functional response is assumed to be the ( N × 1) -dimensional vector Y = [y T 1 , y T 2 , · · · , y T n ] T , and N = m i=1 r i . When the time step is r , then N = mr . The corresponding (N × n)-dimensional design matrix X = [X 1 , X 2 , · · · , X N ] is represented with the Kronecker notation as shown in Eq.  where ⊗ is the Kronecker product operator and 1 r i is the column vector with length r i = 1 . When the time steps are the same, X is an mr × n matrix. where F is written as where R X, t is an N × N correlation matrix, the elements of which are R(X i − X j , t * i − t * j ) . The correlation coefficient θ is optimized with Eq. (13). The optimization of the correlation coefficient requires a large number of solutions of R −1 X, t and R X, t so that the computational cost is extremely high.
In this paper, we consider the functional output on the regular grid, i.e., t 1 = · · · = t n = t , r 1 = · · · = r n = r . The correlation matrix R X, t of the functional response is written as Eq. (22): where R X is an (m × m) -dimensional correlation matrix whose elements are R(x 1 − x 2 ) and R t is an ( r × r) -dimensional correlation matrix whose elements are R T (t i − t j ) . The Gaussian correlation function is used to construct the correlation matrix, and the inverse of the correlation matrix is simplified to: Substituting this matrix into Eq. The computational complexity of the inversion matrix is reduced from O(n 3 m 3 ) to O(n 3 + m 3 ) with the Kronecker product operator. (19) x 11 x 12 · · · x 1n . . . . . . . . . . . .
Eq. (16) takes account of the coupling effect between the input variable x and time t . The basis function F (x, t) is identified by estimating the average value of the functional response.
In the case of r 1 = · · · = r n = r , we define e ·j = 1 m m i=1 (y ij − y i· ) and y i· = 1 r r l=1 y il . e = [e ·1 , · · · , e ·r ] T and y = [y 1· , · · · , y m· ] are fitted as shown in Eq. (24): The unknown polynomial basis functions k(x) and g(x) are constructed using the regression function. Furthermore, the basis function in Eq. (23) is written as

Deformation Modeling of Needle Insertion at a Series of Depth
One hundred sets of input variables are generated with LHS to construct the functional response Kriging model. After the normalized inverse operation, the finite element program runs 100 times to generate the deformations in the X and Y directions at 10 marked points for 31 time steps, as shown in Figure 2. For needle insertion at a constant speed, the insertion time is replaced with the depth d . Here, the 31 time steps are expressed as d = [1.0 mm, 2.0 mm, 5.0 mm, · · · , 89 mm] , and the increment of the insertion depth is 3 mm , except in the first step. Forty sets of data are randomly selected from 100 sets of records and used to construct the functional Kriging prediction model of soft tissue deformation, and an independent set of data is used to test the model. According to the former analysis, the polynomial basis functions are first order regression functions, and the correlation functions are Gaussian random functions. To evaluate the real-time performance of the Kriging module, the mean residual e i and the relative mean residual ε i at each time step (insertion depth) are defined for each observation location, written as Eq. (25): where r is the total number of puncturing steps, i.e., r = 31 , and i is the position number of the output response, where i = 1, · · · , 10.  The procedure of needle insertion is recorded from the time of contact between the soft tissue and the needle tip to the time at which a depth of 90 mm is reached inside the tissue. With the functional response Kriging model, the displacement in the X and Y directions can be predicted, and Figure 6 plots the displacements of N2 shown in Figure 2, which is close to the needle. In Figure 6, 'Kriging value' denotes the prediction results obtained using the functional response Kriging model, and 'FEM' denotes the results obtained with a finite element simulation. From Figure 6, the displacement fluctuations at the initial and final time steps are relatively large, while the fluctuations between these two time steps are smaller. The functional response Kriging model follows the law of tissue deformation during the needle insertion procedure. The model can also predict the displacements at other locations at different time steps very well. Figure 7 shows the displacements at N8, which is far away from the needle.
With Eq. (25), the mean residuals and relative mean residuals of the 10 locations are as listed in Table 2.
The fluctuation occurs at the initial time steps which is derived from the spring-back of the tissue surface as the needle punctured the tissue. And as the Kriging model is inherently a data-driven model which the interpolation function is applied to predict other unknown region, the prediction accuracy of the data which has large fluctuation would decrease. Overall, the Kriging model smooths the original finite element model, and the prediction accuracy of the initial stage is poor. The maximum residuals of the Kriging predictions and finite element calculations in the X and Y directions are 0.37 mm and 1.24 mm , respectively, and both occur in the fifth time step, at which the needle tip punctures the surface of the soft tissue. The average residuals in the X and Y directions occur within 0.05 mm and 0.12 mm , respectively, and the relative mean residual is at most 16%. The runtime is also compared. Both the Kriging model and FEM model run on the same computer with macOS Mojave 10.14.6, a 2.7 Hz Intel Core i5 and 8 GB of RAM, and the program is written in MATLAB 2019b. The average runtime of 5 runs is recorded, with the average time for the Kriging-based method being 0.0294 s and that of the FEM model being 4.6912 s. Note that neither of the two reported times includes a visualization of the soft tissue deformation process. The reported results indicate that the Kriging model runs approximately 160 times faster than the FEM model.

Adaptation of the Deformation Modeling to Step Variations
From the above analysis, the Kriging prediction model based on the observations at 31 time steps can reflect the soft tissue deformation at fixed (known) time steps. However, it is necessary to obtain the motion information of the sensitive area inside the soft tissue at any time step for both the path planning and the virtual training system. In this section, the adaptation to different time steps is studied. Additional insertion time steps are expressed as d 1 = [1.0 mm, 2.5 mm, 5.0 mm, · · · , 90 mm] , the insertion depth increment is 2.5 mm and the total depth achieved by the needle is 90 mm . The displacements of 10 observation locations at d 1 insertion time steps are recorded to compare with the Kriging model. Forty sets of data randomly selected from 100 sets of records are used to construct the real-time Kriging prediction model of soft tissue deformation, and an independent set of data is used to test the model. Figure 8 shows a comparison between the prediction and the computer-based experimental results at the N3 location. From Figure 8, the Kriging model with  a functional response can fit the computer-based experimental results very well and can reflect the fluctuation trends at the initial and end time steps. The maximum residual occurs at the 6th insertion time step, and its maximum residuals and relative mean residuals in the X direction are 0.22 mm and 21% , respectively. The maximum residuals and relative residuals in the Y direction are 0.55 mm and 8% , respectively. The maximum residuals and relative mean residuals of the remaining nine observation positions are listed in Table 3. From Table 3, the maximum residual in the X direction occurs at N1 and N7; its value is 0.31 mm . The maximum residual in the Y direction appears at N2; its value is 0.92 mm . The relative mean residual in the Y direction is less than 8% and has a relatively high accuracy. For the real insertion procedure, the displacements in the X direction are very small, sometimes approaching zero; hence, the relative mean residual is larger than 20% . However, the absolute mean residual of each location is less than 0.31 mm , which can meet the estimation accuracy.

Conclusions and Future Work
In this paper, a computer-based experimental analysis of the Kriging metamodel is presented to predict the deformations of soft tissues. The input data include the material properties of the tissues and needle, the geometrical properties of the needle, and the solver parameters, and they are sampled by the LHS method. The corresponding output dataset is generated by an accurate needle-tissue coupling FE model offline. Ten markers are used to represent the deformations at different positions inside the tissues. Unlike other mechanical-based simulation models, our model can consider both the accuracy and efficiency if high-precision data of the original dataset are well prepared. The results suggest that the proposed Krigingbased model with first order regression and Gaussian correlation functions can well reflect the mechanism of soft tissue deformation. The functional response Krigingbased model can provide feedback regarding deformation at a series of depths for the SPS, of which the time and space indices are both taken into account. The reported time indicates that the Kriging model runs approximately 160 times faster than the original FEM model. Moreover, we also verify that the model has excellent adaptation to step variations. Compared with the ground truth, the maximum residual is less than 0.92 mm , which can satisfy the requirements of the SPS.
The performance of the proposed Kriging-based model depends on the accuracy of the dataset generated by the mechanical-based simulation greatly. Therefore, future work will focus on combinations with a more accurate simulation model, considering complex tissue behavior ranging from hyperelasticity to viscoelasticity, and the needle deflection will be modeled at the same time. The sensitivity of additional parameters will be analyzed to improve the accuracy and efficiency of the Kriging-based model.