Abstract
Online soft tissue characterization is important for robotic-assisted minimally invasive surgery to achieve precise and stable robotic control with haptic feedback. This paper presents a new nonlinear recursive adaptive filtering methodology for online nonlinear soft tissue characterization. An adaptive unscented Kalman filter is developed based on the Hunt-Crossley model by windowing approximation to online estimate system and measurement noise covariances. To improve the accuracy of noise covariance estimations, a recursive formulation is subsequently developed for estimation of system and measurement noise covariances by introducing a weighting factor. This weighting factor is further modified to accommodate noise statistics of large variation which could be caused by rupture events and geometric discontinuities in robotic-assisted surgery. Simulations, experiments, and comparison analyses demonstrate that the proposed nonlinear recursive adaptive filtering methodology can characterize soft tissue parameters in the presence of system or measurement noise statistics in both small and large variations for robotic-assisted surgery. The proposed methodology can effectively estimate soft tissue parameters under system and measurement noises in both small and large variations, leading to improved filtering accuracy and robustness in comparison with UKF.
1. Introduction
Soft tissue properties are of great importance in robotic minimally invasive surgery to characterize the interaction between surgical tools and soft tissues for robotic control with force feedback. However, soft tissue properties are time-dependent and patient-specific. They are variable as per various tissue layers and regions, various organs, and various physiological conditions. In order to improve the robustness and stability of haptic control, it is necessary to characterize soft tissue properties dynamically [1].
A contact model describing the contact interaction between surgical tools and soft tissues is essential to characterization of soft tissue properties. The existing contact models can be generally classified into two categories. One is the continuum mechanics approach focusing on accurate characterization of the tool-tissue contact interaction. The typical examples in this category include the viscoelastic model [2] and finite element model (FEM) [3], where the mechanical behaviours of a soft tissue are accurately characterized based on continuum mechanics of elasticity. However, this approach is very complex and involves a large amount of computational load. Thus, it is only suitable for offline analysis, while unable to satisfy the real-time requirement of soft tissue characterization. The other is the analytical approach focusing on analytical establishment of the explicit relationship between force and displacement using a spring and damper. In spite of the less accuracy compared to the continuum mechanics approach, the analytical approach has the advantages of real-time performance, suitable for a robotic control purpose. Further, since the spring and damper represent the elasticity and viscosity of soft tissues, the approach still maintains certain accuracy for characterization of soft tissue behaviours. Thus, the existing studies are mainly dominated by the analytical approach, leading to various spring-damper models for online soft tissue characterization. The simple one is made up of a spring, unable to model the viscosity effect [4]. The Maxwell model consists of a serialized spring and damper [5]. The Kelvin-Voigt model consists of a parallel spring and damper [6]. The Kelvin-Boltzmann (KB) model improves the K-V model adding a serialized spring [7]. In general, the above models are a linear model, involving physical inconsistencies at the initial and end stages of contact with soft tissues [4]. Comparing to the above models, the H-C model is a nonlinear contact model and shows consistencies with physical contact. However, the use of such a nonlinear model requires a nonlinear estimation algorithm, which is a challenging research task. Accordingly, the use of the H-C model for accurate and real-time characterization of soft tissue properties has been very limited. A survey of the existing contact models can be found in [8, 9].
In order to achieve the real-time performance for soft tissue characterization, it also needs to develop an online estimation algorithm. So far, the existing studies are mainly dominated by a linear estimation algorithm such as the recursive least square (RLS) [10, 11] and Kalman filter [8]. However, a linear estimation algorithm is only suitable for linear contact models. In order to use a linear estimation algorithm with the nonlinear H-C model, it requires a linearization process. Since the linearization process of the nonlinear H-C model causes error for estimation, the resultant solution may be biased or even divergent [12]. Further, the Kalman filter also requires prior knowledge of system noise statistics. A survey on estimation and tuning of system noise statistics for Kalman filtering can be found in [13, 14].
The unscented Kalman filter (UKF) is an online nonlinear estimation algorithm, where the state mean and covariance are approximated in third-order accuracy using unscented transformation [15]. However, similar to the Kalman filter, the performance of UKF relies on accurate contact and measurement models. Currently, there have been few studies focusing on development of a UKF algorithm for characterization of nonlinear soft tissue properties. Xi et al. studied a reduced-order UKF based on a cubic Hermite FEM for characterization of myocardial stiffness [3]. However, this FEM is subject to the quasistatic assumption, which does not account for the contact dynamics. The reduced dimension of the state vector also deteriorates the estimation performance. The authors also studied a UKF based on the H-C model for soft tissue characterization [16], but without considering the UKF requirement of accurate system models. Just recently, the authors reported an improved UKF to handle the error of the H-C contact model using the scaling factor [17]. However, this method is based on the assumption that the measurement model is accurate, which does not correspond to the real-world situation with errors involved in both contact and measurement models. Further, it is also based on windowing approximation [18, 19] to use residuals within a small window of historical time points to calculate system noise characteristics at the current time point. In spite of the fast response to dynamic model noise, the estimation accuracy is limited due to the restricted data in the small time window. Recursion [20, 21] is a strategy to assess system characteristics at the present time point using all historical data, leading to higher accuracy than windowing approximation. Since all data are assembled via a recursive process, rather than stored in memory, this strategy also has a small computational load. However, compared to windowing approximation, the recursion strategy is slow in response to model noise due to the use of all available data. As windowing approximation and recursion complement with each other, it is necessary to study how to effectively combine both strategies together to fully leverage their individual advantages for improvement of the UKF performance.
This paper presents a new nonlinear recursive adaptive UKF to estimate the parameters of the H-C contact model. It combines windowing approximation with recursion to address the UKF limitation in requirement of accurate contact and measurement models. An adaptive UKF is developed based on windowing approximation to offer the classical UKF with the capability to online estimate the noise covariances of the H-C contact and measurement models. Subsequently, a recursive adaptive UKF is developed by incorporating a weighting factor in the adaptive filtering process to improve the estimation accuracy. This weighting factor is further modified to accommodate the noise statistics of large variation, which are occurring in robotic-assisted surgery due to rupture events and geometric discontinuities. Simulation and experimental results together with comparison analysis demonstrate the efficiency of the proposed method for online soft tissue characterization.
2. Hunt-Crossley Model
The Hunt-Crossley (H-C) model describes the dynamic mechanism of mechanical contact between soft tissues and surgical tools with a nonlinear spring and damper [22]: where , , , , , , and denote the interaction force, stiffness, damping coefficient, displacement, velocity, and power indices of displacement and velocity, respectively.
From (1), the system state equation is defined as where is the system state at time , is the system noise which is assumed as a white Gaussian noise with zero mean and covariance , is the system function, and is the time step.
The measurement function is defined as where denotes the measurement vector, is the measurement matrix, and is the measurement noise which is also assumed as a white Gaussian noise with zero mean and covariance .
3. Analysis of Conventional UKF
Let us briefly analyse the problems of the classical UKF. The conventional UKF procedure is as follows: (i)Initialization.
Calculate the initial estimated state and its associated covariance: where is the predefined initial state. (ii)Time update.
Select the sigma points of the estimated state at time (): where is the dimension of state vector , is the tuning parameter to control the spread of the sigma points around , and is the th column of the mean square root of .
Transform the sigma points with the state function, yielding a set of transformed samples:
Calculate predicted state and its associated covariance where (iii)Measurement update.
As the measurement model is linear, the measurement update is performed with the same equations as the Kalman filter.
Calculate predicted measurement and its associated covariance
Calculate the crosscovariance between predicted state and predicted measurement :
Determine the Kalman gain:
Update estimated state and its associated covariance : (iv)Go to step (i) for the next sample until all samples are processed
It can be seen from (8) that if the system noise covariance involves uncertainty, predicted state covariance will become inaccurate, leading crosscovariance given by (12) to become inaccurate. The inaccurate will further lead Kalman gain to become biased, thus degenerating the state estimate given by (15). In a similar way, if the measurement noise covariance involves uncertainty, predicted measurement covariance will become inaccurate, leading Kalman gain to be biased. Accordingly, the state estimate given by (14) will be corroded. From above, it is evident that without the accurate system noise covariances, the filtering solution of the conventional UKF will be biased or even divergent.
4. Recursive Adaptive Unscented Kalman Filter
In this paper, a new recursive adaptive UKF is developed by combining windowing approximation with recursion to the online estimate system and measurement noise covariances via covariance matching, which is to balance predicted state covariance and predicted measurement covariance with their theoretical values.
4.1. Estimation of System Noise Covariance
By introducing an adaptive scaling factor , we readily have the following relationship:
Substituting (16) into (8), the predicted state covariance can be further represented as
It can be seen from (17) that the adaptive scaling factor enables us to adjust the system noise covariance for online tuning the predicted state covariance.
In theory, the predicted state covariance is defined as
Since the real state is generally unknown, in (18) replacing the real state with the estimated state and applying windowing approximation to calculate the expectation, we readily have a suboptimal estimation for the predicted state covariance: where represents the number of the state vectors within the time window (, ) and is named the window size.
Define the innovation vector by
By (14), we readily have
Substituting (21) into (19) yields
Based on the above, the adaptive scaling factor can be attained by matching the predicted state covariance given by (17) with its theoretical value given by (22):
Thus, the adaptive scaling factor for the system noise covariance can be determined as where denotes the trace of a matrix.
It can be seen from (22) that is approximated using the information within the current window, and thus, it can account for dynamically changing conditions. However, due to the limited data available within the time window, the estimation accuracy is limited. In order to address this problem, the theoretical predicted state covariance is further expressed into a recursive form: where is called the weighting factor, which is commonly defined by [23, 24]: where is the number of the entire available time points.
It can be seen from (25) that is based on the entire available data rather than the limited data available within the current time window, leading to improved accuracy. is contributed by two portions. One is at the current time point, which is calculated by windowing approximation to use the data available within the current window to characterize dynamically changing conditions that occurred. The other is at the previous time point, which reflects the historical tendency. The contributions at both current and previous time points are combined via weighting factor . takes the full contribution at the current time point when , while taking the full contribution at the previous time point when . The smaller the value of is, the smaller the contribution at the current time point is and in turn the larger the contribution at the previous time point is.
4.2. Modified Weighting Factor
In robotic-assisted minimally invasive surgery, the variation of noise statistics is small when the dynamic interaction environment is relatively stable [25]. However, in case of abrupt changes such as geometrical discontinuities in soft tissue deformation, rupture events occurring in the penetration of different tissue layers, and unexpected contacts with rigid bones, a large variation in noise statistics will occur in the dynamic interaction environment [26].
The standard weighting factor defined by (26) cannot account for abrupt changes encountered in the interaction with soft tissues. Figure 1 plots the variation of the standard weighting factor in time series where an abrupt change in noise statistics has occurred at s. As shown in Figure 1(a), given by (26) starts from the maximum value of one at the beginning. Subsequently, it decreases drastically to a small value and gradually converges to zero in the time series. As remains at such a small value most of the time, the contribution at the present time point to is small while in turn the contribution at the previous time point is large. Due to the small contribution at the present time point most of the time, is unable to account for an abrupt change in system noise statistics. As shown in Figure 1(a), as still remains at a small value at time point s, it leads to the failure to account for the abrupt change of noise statistics occurring at this time point. Accordingly, given by (26) is only suitable for the case that system noise statistics are constant or have a small variation, making the proposed adaptive UKF unable to account for a large variation in system noise statistics.
(a) Standard weighting factor
(b) Modified weighting factor
In order to accommodate large variations in system noise statistics, the standard weighting factor defined by (26) is modified as where is a unit step function and is the index of the time point where an abrupt change occurs.
The initial value of parameter is zero. When an abrupt change occurs at the current time point, the parameter will be updated to , making become one again. Consequently, is calculated by taking the full contribution at the present time point to account for the abrupt change. Figure 1(b) shows the variation of the modified weighting factor where an abrupt change in system noise statistics occurred at s. Similar to , starts from the maximum value of one, where takes the full contribution at the present time point. Subsequently, the value drops drastically to and remains at a small value until an abrupt change occurs. With such a small value of , mainly follows the one at the previous time point to handle small variations in system noise covariance. When an abrupt change occurs at the time point s ( s), the parameter defined by (28) will be updated to , restoring back to the maximum value of one to account for the abrupt change that occurred. Subsequently, drastically drops to and remains at a small value until a new abrupt change occurs. This variation process of is repeated whenever an abrupt change occurs, allowing the adjustment of to account for an abrupt change by taking the full contribution at the present time point.
By substituting the weighting factor with the modified weighting factor, (25) becomes
Thus, the adaptive scaling factor is calculated by
4.3. Estimation of Measurement Noise Covariance
The measurement noise covariance can also be estimated in the similar way as the system noise covariance. By introducing an adaptive scaling factor , the measurement noise covariance can be represented as
Substituting (31) into (11), the predicted measurement covariance can be further expressed as
It can be seen from (32) that the adaptive scaling factor enables us to adjust the measurement noise covariance for online tuning the predicted measurement covariance.
Using windowing approximation [27], the theoretical predicted measurement covariance can be directly calculated as
Similar to (29), the theoretical predicted measurement covariance can be further expressed by a recursive form as where the modified weighting factor is also applied to handle both small and large variations in the measurement noise covariance.
By matching the modified predicted measurement covariance given by (32) and its theoretical value given by (34), we readily have
From (35), the adaptive scaling factor for estimation of measurement noise covariance is determined as
5. Performance Evaluation and Discussions
A prototype system has been implemented using the proposed methodology for online soft tissue characterization. This system takes measurements of contact force and displacement as input signals to online estimate the nonlinear parameters of the H-C model. Based on the estimated parameters of the H-C model, the contact forces with soft tissues are reconstructed from displacements and further compared with the input forces to calculate the estimation error. Simulations and experiments together with comparison analysis have been conducted to evaluate the performance of the proposed methodology in terms of (i) system noise covariance estimation, (ii) measurement noise covariance estimation, (iii) the effect of the weighting factor, and (iv) force reconstruction based on the estimated parameters of the H-C model.
5.1. Simulations and Analysis
Simulations were conducted to evaluate the performance of the proposed methodology in estimating system and measurement noise covariances together with associated force reconstructions. The time steps for the simulations were set to 10 ms. The initial state was defined as .
5.1.1. System Noise Covariance Estimation
For the performance evaluation in terms of system noise covariance estimation, it is supposed that the measurement noise statistics are exactly known. Without loss of generality, set and . Two cases are studied. One is that system noise covariance is constant or has a small variation, and the other is that system noise covariance involves abrupt changes.
(1) Constant or Small-Variation System Noise Covariance. Simulation trials were conducted under the same conditions by the adaptive UKF (AUKF) and recursive adaptive UKF (RAUKF) for the case that system noise covariance is constant or has a small variation. The measurement data were obtained from the H-C model by adding a white Gaussian noise with zero mean and the covariance of 0.036 N. The initial estimation value was N, the window size was , the real covariance of system noise was N, and . The simulation time was 90 s. As shown in Figure 2, the system noise covariance estimated by AUKF via windowing approximation involves large oscillations during the entire simulation time. This is because the limited data available within the small window restricts the estimation accuracy. In contrast, the system noise covariance estimated by RAUKF via recursion is close to the real value after 20 s. This is because the estimation via recursion is obtained based on the entire available data. The RMSE (Root Mean Square Error) of AUKF is 0.0125 N, while after 20 s, the RMSE of RAUKF is 0.0032 N. Thus, it is evident that the proposed RAUKF can significantly improve the AUKF accuracy for estimation of system noise covariance.
(a) AUKF via windowing approximation
(b) RAUKF via recursion
Figure 3 illustrates the estimation errors in terms of force reconstruction by both AUKF and RAUKF, based on the estimated system noise covariances as shown in Figure 2. It is clear that the estimation error of RAUKF is much smaller than that of AUKF. During the entire simulation time, since the system noise covariance estimated by AUKF involves large oscillations (see Figure 2(a)), the resultant estimation error also involves large oscillations. However, for RAUKF, the resultant estimation error becomes small and is further converged after the initial time period of 20 s, following the trend of the estimated system noise covariance. The maximum estimation error is 0.0394 N for RAUKF, but 0.1853 N for AUKF. The mean estimation error of AUKF is 0.0374 N, while that of RAUKF is 0.0053 N. The RMSE of RAUKF is 0.0070 N, while that of AUKF is almost nine times larger, which is 0.0475 N. Table 1 summarizes the estimation errors of both AUKF and RAUKF.
(2) System Noise Covariance Involving Abrupt Changes. In order to analyse the effect of the weighting factor on system noise covariance estimation, simulation trials were conducted under the same conditions by RAUKF via both standard and modified weighting factors for the case of system noise covariance involving abrupt changes. The input signals were obtained from the H-C model by adding a white Gaussian noise with zero mean. The noise covariance was set to 0.036 N within [0, 36 s] and 0.576 N within (36 s, 70 s], leading to an abrupt change from 0.036 N to 0.576 N at 36 s. The initial estimation value was N, the window size was set to , and the simulation time was 70 s.
Figure 4 shows the system noise covariances estimated by RAUKF via both standard and modified weighting factors. Although the estimation results with both weighting factors follow the reference value of 0.036 N within [0, 36 s], the estimation under the standard weighting factor does not follow the abrupt change from 0.036 N to 0.576 N at 36 s, leading to the mean error of 0.372 N after 36 s. In contrast, the estimation of system noise covariance under the modified weighting factor rapidly follows the abrupt change from 0.036 N to 0.576 N at 36 s, leading to the mean error of 0.0088 N after 36 s. This demonstrates that with the modified weighting factor, RAUKF can accommodate abrupt changes in system noise covariance.
(a) RAUKF with the standard weighting factor
(b) RAUKF with the modified weighting factor
Figure 5 shows the estimation errors in force reconstruction based on the estimated system noise covariances as shown in Figure 4. It can be seen that before 36 s, although the RAUKF estimation errors via both standard and modified weighting factors are small and converged to zero, the one via the modified weighting factor is even smaller and has a faster convergence speed. However, after 36 s, due to the biased estimate of the system noise covariance, the RAUKF estimation error via the standard weighting factor involves large oscillations, leading to the RMSE of 0.0749 N. In contrast, the RAUKF estimation error via the modified weighting factor is much smaller and converged to zero, leading to the RMSE of 0.0046 N. This demonstrates that the modified weighting factor enables RAUKF to account for abrupt changes in system noise covariance, leading to improved estimation accuracy. Table 2 summarizes the RAUKF estimation errors via both standard and modified weighting factors.
5.1.2. Measurement Noise Covariance Estimation
In order to evaluate the estimation of measurement noise covariance, it is assumed that the system noise statistics are exactly known. Without loss of generality, choose and . Similar to Section 5.1.1., two cases are studied. One is that the measurement noise covariance is constant or in a small variation, and the other is that the measurement noise covariance involves abrupt changes.
(1) Constant or Small-Variation Measurement Noise Covariance. Simulation trials were conducted under the condition that the measurement noise covariance is constant or has a small variation by both AUKF and RAUKF.
The input signals were obtained from the H-C model by adding a white Gaussian noise with zero mean and the covariance of 0.036 N. The initial estimation value was N. The window size was set to . The simulation time was 70 s. As shown in Figure 6, the estimate of measurement noise covariance by AUKF involves large oscillations during the entire simulation time, while the one by RAUKF is very close to the reference value after 15 s. The RMSE of RAUKF is 0.0051 N and is almost twice smaller than that of AUKF which is 0.0101 N. This demonstrates that the accuracy of RAUKF is significantly higher than that of AUKF for estimation of measurement noise covariance.
(a) AUKF via windowing approximation
(b) RAUKF via recursion
Figure 7 illustrates the estimation errors in force reconstruction by both AUKF and RAUKF based on the estimated measurement noise covariances as shown in Figure 6. It can be seen clearly that the estimation error of AUKF is much larger than that of RAUKF during the entire simulation time. In addition, after 15 s, the estimation error of RAUKF is close to zero, while that of AUKF is slightly smaller than the initial estimation error. This is because the estimated measurement noise covariance by windowing approximation is fluctuated during the entire simulation time (see Figure 6(a)), while that by RAUKF is converged to the reference value (see Figure 6(b)). The RMSE of RAUKF is 0.0286 N and is two-three times smaller than that of AUKF which is 0.0708 N. Further, the maximum and mean estimation errors of AUKF are about two times larger than those of RAUKF. The maximum estimation error is 0.2905 N for AUKF, while 0.145 N for RAUKF. The mean estimation error is 0.0551 N for AUKF, while 0.0102 N for RAUKF. Therefore, the estimation accuracy of RAUKF is significantly higher than that of AUKF. Table 3 compares the estimation errors of both AUKF and RAUKF.
(2) Measurement Noise Covariance Involving Abrupt Changes. In order to analyse the effect of the weighting factor on measurement noise covariance estimation, trials were conducted under the same conditions by RAUKF via both standard and modified weighting factors for the case of measurement noise covariance involving abrupt changes. The input signals were obtained from the H-C model by adding a white Gaussian noise with zero mean. The covariance of this measurement noise was set to 0.036 N within [0, 32 s] and 0.576 N within (32 s, 70s], leading to an abrupt change at 32 s. The initial estimation value was N. The window size was . The simulation time was 70 s. Figure 8 shows the system noise covariances estimated by RAUKF via both standard weighting factor and modified weighting factor. As shown in Figure 8, the estimations of measurement noise covariance under both weighting factors follow the reference value of 0.036 N within [0, 32 s]. However, the estimation under the standard weighting factor does not follow the abrupt change from 0.036 N to 0.576 N at 32 s, leading to the large mean error of 0.2662 N after 32 s. In contrast, the estimation under the modified weighting factor follows closely the abrupt change, leading to the mean error of 0.043 N after 32 s.
(a) RAUKF with the standard weighting factor
(b) RAUKF with the modified weighting factor
Figure 9 shows the estimation errors in force reconstruction based on the estimated measurement noise covariances as shown in Figure 8. Before the abrupt change at 32 s, both standard and modified weighting factors result in the similar estimation error. The mean error is 0.0025 N by the standard weighting factor and 0.0023 N by the modified weighting factor. However, after the abrupt change at 32 s, due to the biased estimate of the measurement noise covariance, the use of the standard weighting factor results in a large magnitude of oscillations in the estimation error curve, leading to the RMSE of 0.0907 N. In contrast, with the modified weighting factor, the magnitude of oscillations in the estimation error curve after 32 s is significantly decreased and gradually converged to zero, leading to the RMSE of 0.0444 N. Therefore, it is proved that with the modified weighting factor, RAUKF is able to accommodate abrupt changes in measurement noise covariance. Table 4 lists the estimation errors of RAUKF via both standard and modified weighting factors, where the mean error and maximum error via the standard weighting factor are 0.0603 N and 0.5493 N, while 0.0284 N and 0.2961 N via the modified weighting factor.
5.2. Experiments and Analysis
The performance evaluation is also conducted on two experimental cases in the presence of abrupt changes. One is robotic-assisted needle insertion in the presence of rapture events, where the measurement data was acquired from the literature. The other is the mechanical indentation in the presence of geometrical discontinuities to acquire the data of mechanical load and displacement on a phantom tissue.
5.2.1. Robotic-Assisted Needle Insertion
Trials were conducted to evaluate the performance of the proposed methodology for the case of robotic-assisted needle insertion. The experimental data of force and displacement on biological soft tissues during the process of needle insertion were obtained from the literature [28]. The input force profile involves four peaks at the displacements of around 20 mm, 29 mm, 32 mm, and 36 mm (see Figure 10). These peaks represent the rupture events encountered by the robotic needle during the insertion process. For comparison analysis, trials were conducted under the same conditions by UKF, AUKF, and RAUKF with the modified weighting factor. The initial state and noise covariance were set as , , and . The window size for both AUKF and RAUKF was set to .
(a) UKF
(b) AUKF
(c) RAUKF
Figure 10 shows the reconstructed forces by UKF, AUKF, and RAUKF, respectively. It can be seen that the UKF estimation results in large errors at the four peaks, where the estimation error at the first peak (at displacement of around 20 mm) is most significant, which is 0.738 N. An obvious deviation also remained in the displacement range from around 40 mm to 50 mm after the four peaks, leading to the mean error of 0.068 N and the RMSE of 0.15 N.
AUKF improves the estimation accuracy of UKF at the four peaks, leading to the mean error of 0.044 N and RMSE of 0.069 N. However, there are still pronounced errors at and after the four peaks. This is because the estimated noise covariance by AUKF involves large error as shown in the previous simulations. In contrast, the reconstructed forces by RAUKF are very close to the reference values at and after the four peaks, leading to the mean error of 0.020 N and RMSE of 0.039 N. Table 5 lists the mean and maximum errors as well as RMSEs of the three methods. It can be seen that although both AUKF and RAUKF outperform UKF, RAUKF has the smallest estimation errors (the mean error of 0.020 N, maximum error of 0.205 N, and RMSE of 0.039 N). The mean and maximum errors of RAUKF are more than three times smaller than those of UKF. The RMSE of RAUKF is almost four times smaller than that of UKF. The mean and maximum errors as well as RMSE of RAUKF are also much smaller than those of AUKF.
5.2.2. Mechanical Indentation with Geometrical Discontinuities
Mechanical indentation tests were also conducted to evaluate the performance of the proposed methodology for soft tissue characterization in the presence of geometrical discontinuities. The indentation tests were conducted on a phantom tissue sample using the DMA (Dynamic Mechanical Analyser, Seiko Instruments). As shown in Figure 11, the phantom tissue is made up of silicone rubbers (Ecoflex 0030), which have the similar characteristics as human soft tissues [29]. The phantom tissue is of cubic shape () to fit into the DMA machine. An indenter of 1 cm diameter is used to compress the phantom tissue vertically. In order to present geometrical discontinuities, the indenter is controlled to compress the phantom tissue with a given displacement and then held at the displacement for around 0.9 s. Subsequently, the indenter is controlled with the maximum speed (1,000,000 μm/min) to increase the displacement by 1 mm and then held at the new displacement for around 0.9 s, with the conduction of three times. The data of displacement and load were acquired from the DMA during the compression tests. As shown in Figure 12, the force profile has three abrupt changes (three steps) at around 0.9 s. 1.8 s, and 2.7 s, which represent geometrical discontinuities in robotic-assisted surgery. Trials were conducted based on measured forces and displacements under the same conditions by UKF, AUKF, and RAUKF with the modified weighting factor. The parameter settings for these methods are the same as that in the previous case of robotic-assisted needle insertion.
(a) The indenter and phantom tissue
(b) The appearance of DMA
(a) UKF
(b) AUKF
(c) RAUKF
Figure 12 shows the reconstructed forces by UKF, AUKF, and RAUKF from measured forces and displacements. It can be seen that the UKF estimation has large estimation errors at each abrupt change, leading to the maximum error of 0.6341 N, mean error of 0.0099 N, and RMSE of 0.0157 N, respectively. Comparing to the UKF estimation, the AUKF estimation follows the reference force curve, without involving a large estimation error. The mean error of AUKF is 0.0064 N, which is smaller than that of UKF. However, due to large error involved in noise covariance estimation, the estimation error of AUKF becomes large after around 2 s with a large peak error at 2.5 s where the input force is changed abruptly. However, due to the use of recursion, the estimation of RAUKF closely follows the reference force curve, especially at the three abrupt changes. Its RMSE is 0.0098 N, which is much smaller than those of UKF and AUKF. Table 6 compares the estimation errors of the three methods. It can be seen that in addition to the mean error, the maximum error and RMSE of RAUKF are also much smaller than those of UKF and AUKF.
6. Conclusions
This paper presents a new methodology for online nonlinear soft tissue characterization in robotic-assisted surgery. An adaptive UKF is established based on the nonlinear H-C model for online estimate soft tissue parameters without prior knowledge on noise statistics. Based on this, a recursive adaptive UKF with the modified weighting factor is further developed to improve the estimation accuracy and accommodate large variations in noise statistics. Simulation and experimental results as well as comparison analysis demonstrate that the proposed methodology can effectively estimate soft tissue parameters under system and measurement noises in both small and large variations, leading to improved accuracy and robustness in comparison with the conventional UKF.
Future work will focus on two aspects. One is to investigate the stability of the proposed methodology, including the effects of window size and initial state error on the filtering performance. The other is to incorporate the proposed methodology in robotic control and force feedback for robotic-assisted minimally invasive surgery. It is expected that new control strategies will be developed by integrating the proposed methodology into the control loop for precise and stable robotic control with force feedback in the presence of unknown system noises.
Nomenclature
| : | System state |
| : | Estimated state |
| : | Predicted state |
| : | th sigma point of estimated state |
| : | th sigma point transformed by system function |
| : | Real state |
| : | Predicted measurement |
| : | Measurement |
| : | System noise |
| : | System noise covariance |
| : | Measurement noise |
| : | Measurement noise covariance |
| : | Estimated state covariance |
| : | Predicted state covariance |
| : | Predicted state covariance in theory |
| : | Predicted measurement covariance |
| : | Predicted measurement covariance in theory |
| : | Crosscovariance of predicted state and measurement |
| : | Kalman gain |
| : | Adaptive scaling factor for system noise covariance estimation |
| : | Adaptive scaling factor for measurement noise covariance estimation |
| : | Dimension of state |
| : | Tuning parameter |
| : | Time step |
| : | Interaction force |
| : | Stiffness coefficient |
| : | Damping coefficient |
| : | Displacement |
| : | Velocity |
| : | Power of velocity |
| : | Power of displacement |
| : | Size of window |
| : | Number of time points |
| : | Weighting factor |
| : | Modified weighting factor |
| : | Time point when abrupt changes occur |
| : | Parameter for modified weighting factor |
| : | Innovation vector |
| : | Expectation |
| : | System function |
| : | Measurement matrix |
| : | Unit step function |
| N: | Newton. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
The work of this paper is supported by the Australian Research Council (ARC) Discovery Early Career Award (DECRA) (DE130100274).