Model-Based Thermometry for Laser Ablation Procedure Using Kalman Filters and Sparse Temperature Measurements

— Objective: We implement a data assimilation Bayesian framework for the reconstruction of the spatiotemporal proﬁle of the tissue temperature during laser irradiation. The predictions of a physical model simulating the heat transfer in the tissue are associated with sparse temperature measurements, using an Unscented Kalman Filter. Methods: We compare a standard state-estimation ﬁltering procedure with a joint-estimation (state and parameters) approach: whereas in the state-estimation only the temperature is evaluated, in the joint-estimation the ﬁlter corrects also uncertain model parameters ( i.e. , the medium thermal diffusivity, and laser beam properties). We have tested the method on synthetic temperature data, and on the temperature measured on agar-gel phantom and porcine liver with ﬁber optic sensors. Results: The joint- estimation allows retrieving an accurate estimate of the temperature distribution with a maximal error < 1.5 ◦ C in both synthetic and liver 1D data, and < 2 ◦ C in phantom 2D data. Our approach allows also suggesting a strategy for optimizing the temperature estimation based on the positions of the sensors. Under the constraint of using only two sensors, optimal temperature estimation is obtained when one sensor is placed in proximity of the source, and the other one is non-symmetrical. Conclusion: The joint-estimation signiﬁcantly improves the predictive capability of the physical model. Signiﬁcance: This work opens new perspectives on the beneﬁt of data assimilation frameworks for laser therapy monitoring.


I. INTRODUCTION
L ASER Ablation (LA) has proven to be an important tool in the thermal-based treatment of focal neoplasms. This minimally invasive procedure has the advantage of being relatively low-risk, and has been proposed as an alternative therapy for the localized treatment of some tumors, like liver metastases and pancreatic adenocarcinoma, when the patient cannot undergo surgery or transplant [1]- [3]. LA represents also a ground-breaking therapy in other fields, such as pediatric neurosurgery [4]. Improvement of the procedure for inducing thermal damage within the tumor volume and safety margins, while spearing the surrounding healthy tissue and structures, is still a priority for research groups and companies working in this field. A critical concern in LA is the poor on-line control of the induced thermal effect, being the latter responsible for the effective tumor treatment. The treatment planning is usually based on manufacturer-initiated working algorithms in combination with operator experience. However, the shape and volume of the ablation zone are depending on physical and physiological tissue-specific parameters, such as blood perfusion and temperature-dependent properties [5]. For these reasons, the real-time monitoring of the temperature evolution could improve the therapy outcome and enable patient-specific treatment.
Two approaches for temperature monitoring are mostly used during LA procedures: contactless methods, relying on diagnostic imaging, and contact ones, based on the use of physical sensors. The most common contactless method is the magnetic resonance thermal-imaging (MRTI). Even though magnetic resonance (MR)-compatible fiber optics guide the laser beam inside the MR room, the interaction between the laser light and the tissue may cause the creation of cavitation bubbles. This phenomenon provokes artifacts on the thermometric image, which may impair the accurate temperature measurement of the target [6], [7]. Additionally, measurement errors such as drift of the magnetic field under temperature change and motion still represent a source of uncertainty in the temperature measurement [8]. Contact methods are accurate, highly resolved in space and do not require expensive equipment. In particular, fiber optic sensors with diverse sensing principle, such as fiber Bragg gratings (FBGs) technology and optical frequency domain reflectometry (OFDR) based on Rayleigh scattering, This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ allow for simultaneous multipoint temperature measurement with minimally invasive approach [9], [10]. Whereas FBGs allow for a typical spatial resolution ranging from 1 grating/cm to 1 grating/mm, OFDR-based Rayleigh scattering enables a distributed sensing approach for thermometry, characterized by a theoretical infinite number of sensing points for a given length [11]. It has been recently demonstrated that improvement of spatial resolution for distributed sensors in the sub-millimetric range determines an increase of noise in the signal, thus electing millimeter-resolved FBGs as the suitable choice for accurate thermometry in laser irradiated tissues [12]- [14]. FBGs are also biocompatible, immune from electromagnetic interferences and flexible; they can be easily embedded within needles to be inserted inside the organ under treatment [9]. However, the minimally-invasive clinical context imposes limitations on the ability to render a complete temperature map of the ablated zone. Indeed, typically not more than two or three needles can be inserted to monitor the procedure in real-time [3]. This restriction limits the spatial information that the clinician has on the treatment and the control on the extension of the thermal margins of the target.
Simulations based on the bioheat equation describe the heat transfer in tissues undergoing thermal therapies, but still suffer from the absence or inaccuracy of temperature-dependent and patient-specific parameters. However, a heat transfer model combined with the real-time monitoring of the tissue temperature in specific regions would allow a robust estimate of the procedure state [15], [16]. Following this approach, it would be possible to design also an accurate monitoring stage able to control the tissue temperature at the desired therapeutic levels. Kalman Filter theory is largely used in the biomedical field for estimating the state of the controlled quantity, and permits to manage noisy measurements and potential inaccuracies of the medical device, such as its placement inside the anatomical district [17]. In the scenario of the thermal treatments, Kalman Filter has been proposed to predict temperature during MR-guided thermal therapy delivery in the presence of noise or corrupted data [16], [18], [19]. Image-based thermometry can suffer from inaccuracy related to motion or artifact in the images, but can still benefit from a large number of information within the image field of view. When the procedure is not performed under image-based thermometry, alternative monitoring methods are needed. In this framework, the application of a data assimilation approach based on limited temperature information is of paramount importance.
In this work, we suggest reconstructing the temperature map through a data assimilation Bayesian framework. The predictions of a physical model simulating the heat-transfer in the tissue are associated with information retrieved from the temperature measurements using an Unscented Kalman Filter, to obtain a spatial and temporal augmentation of the temperature. Thermometry is performed with FBG sensors. We compare a standard, state-estimation filtering procedure, where only the temperature is evaluated, with a joint-estimation approach, where uncertain model parameters are also corrected through the filtering loop. We also suggest a strategy for choosing the positions of the sensors in order to optimize the temperature estimation. Our filtering methodology that combines theoretical predictions with temperature measurements is designed to tackle two important challenges which are imposed by the clinical context: at first, the complexity of the targeted biological system limits the scope of purely theoretical temperature predictions based on heat diffusion principles; secondly, the real-time temperature values measured by the FBGs -while being accurateare sparse.

II. MATERIALS AND METHODS
In this section, after some generalities regarding the physical model and the filtering scheme, we focus on three aspects: i) the comparison between simple (state only) and joint (state and parameters) estimation, ii) the tuning of the filter meta parameters, and iii) the criteria related to the placement of the temperature sensors. The experimental setup used to test our approach is also presented.

A. Numerical Model
We consider at first a one-dimensional homogeneous heat model to describe the essential physical and geometrical features of the phenomenon: where T is the temperature, t is the time, D is the tissue diffusivity, x is the spatial coordinate, β is the perfusion value and T b the blood temperature [18]. For laser-irradiated tissues, we consider a supplementary term T l , representing the additional temperature increase due to the laser source, as proposed by [20]. The 1D object is modeled as a grid composed of L-nodes at a distance of Δx (Fig. 1). The explicit discretized equation of the temperature at the i th node, where n represents the current time-step of duration Δt, becomes: We assume fixed (Dirichlet) Boundary Conditions at the grid extremities, corresponding to a constant ambient temperature.

B. Bayesian Filtering
The discretized (2) is an ideal representation of the physical process occurring in the biological tissue and may be affected by uncertainties and parameterization errors. The computed temperature is characterized by a Gaussian error ω i : At the same time, the estimation of the local temperature can also be affected by some random errors ν i , related to intrinsic sensors properties: T i n = T F BG i n + ν i . By combining the predicted temperature computed through the theoretical model with current information provided by the sensors, we can obtain an optimal estimation of the temperature value. Since our model will be validated in the ex vivo scenario, the term β has been set to zero [18].
Based on a general prediction-update scheme, Bayesian filters [21] allow estimating a random state X n by recursively updating its probability density function P(X n ). The random state is described by an uncertain process model X n = f n (X n−1 , ω n ), representing its behaviour in time, and noisy external observations Z n = g n (X n , ν n ), providing partial information on its current state. Given the previous estimation P (X n−1 | Z n−1 ) at discrete time n − 1, the first step consists in computing a prediction of the probability density function by propagation through the process model f n : P (X n |Z n−1 ). Then, the predicted estimate is updated using the current observation, to provide the posterior estimate: P (X n |Z n ). The final state of the system X n is thereafter computed according to some optimality criteria (e.g., Maximum A Posterior or expected value). Within such a Bayesian framework, we want to estimate the temperature of an object by combining an uncertain process model which describes the temperature distribution (3) with noisy external measurements provided by FBG sensors. In this work, we use an Unscented Kalman Filter [22] that applies a sampling based approach to handle a non-linear prediction and observation models. The original Kalman Filters, which were designed for linear systems, have been adapted to non-linear estimation. While Extended Kalman filters require explicit derivative of the model with respect to the parameter which cannot be obtained easily in our case, we opt for Unscented Kalman Filter. In a nutshell, the Unscented Kalman Filter samples the probability space in each step of the estimation and thus requires forward execution of the model (i.e., finite element simulation) for each sample.

C. State and Observations Vector
The state vector takes into account all the variables that are iteratively estimated by the filter. In the classic state-estimation approach, the state vector is composed of the temperature value T i of each node of the grid mesh X n ∈ R L×1 . In the state space, the state equation is presented in (4), shown at the bottom of the page, where ω n (T l , D) is the random error depending on the uncertainties on both the heating source and thermal diffusivity value, and α = D Δt Δx 2 , according to (2). To improve the estimation and the predictive power of the model, it is possible to consider an augmented state. In practice, the physical parameters representing the source of model uncertainties are estimated along with the temperature. In our case, we suppose to have a misknowledge on D, as well as on the parameters defining the additional temperature due to heating source T l : where x i is the distance between the laser and the i th node, σ is the standard deviation of the laser beam (described with a Gaussian distribution), and T l0 is function of the laser power and includes the parameters describing the optical behaviour of the medium [20]. Its purpose is to provide a general and simplified description of any additional changes in the target that are related to the laser-tissue interaction. The main potential sources of error depend on the accurate knowledge of the optical properties of the medium, which are supposed unknown in our approach. The augmented state vector can hence be formalised as: In the state space, the process model may be written as: Using FBG sensors, we retrieve the temperature value of M known locations of our object. The observations vector Z n ∈ R M ×1 is related to the state vector through a constant mapping between corresponding nodes: Z n = MX n + ν n . In practice, M is a diagonal matrix with non-null elements corresponding to the i th nodes of the grid mesh, where the temperature is measured with the sensor; ν n is the random Gaussian noise affecting the measurements. As above stated, both the process and the observation model are affected by some uncertainties and errors expressed respectively by ω n and ν n . Such quantities are taken into account within the prediction-correction loop [23], through the process noise covariance matrix Q n = E[ω n ω t n ] ∈ R N ×N , and the observations noise covariance matrix R n = E[ν n ν t n ] ∈ R M ×M . An appropriate initialization of such parameters is fundamental to an accurate estimate of the state variables.

D. State Vs. Joint Estimation
Classical Kalman Filters are designed to provide estimations of the state variables of the system. However, it is also possible to identify uncertain parameters of the model through the same prediction-correction scheme, via a joint-estimation algorithm. In this approach, the system state is augmented by the parameters of the model, and the filter simultaneously updates the system state together with the parameters. In this process, the corresponding diagonal values of the covariance matrix are initiated with the variances of the parameters.
In our system, three parameters of the (2) are considered uncertain (D, T l0 , 1/(σ 2 )). In the following, we compare the temperature estimation accuracy of state-estimation and jointestimation for all the explored scenarios. The computation time is directly related to the dimension S of the state vector, i.e., we need to perform S simulations per time step to estimate our system. Additionally, the simulation accuracy of (2) depends on the discretization of the domain (Δx). The following values were set: time resolution (Δt) of the simulations equal to 0.2 s, and Δx = 0.5 mm, considering that the ground truth data are provided with a time resolution of 1 s and a spatial resolution of 1 mm. To adapt the filter, the parameters are updated every 5Δt (thus, each 1 s). The code is implemented in Phyton; the characteristics of the system are: Intel(R) Core(TM) processor i7, 9th generation, CPU 2.60 GHz, RAM 16 GB.

E. Process and Measurement Noise
An important task in the Bayesian filtering framework is to determine the appropriate meta parameters, i.e., Q and R. A poor choice of these values may result in reducing the estimation accuracy of the dynamic process. Also, in the case of parameter identification, these values impact the convergence rate; thus, improper tuning results in slow convergence or alternatively overshooting. The value of R, which represents the noise of the measurement process, is relatively easy to determine. In our system, R reflects the accuracy of the temperature measurements, and it is evaluated at 0.1 • C. This value results from the accuracy of the optical spectrum interrogator (1 pm accuracy corresponding to 0.1 • C) [24]. Q represents the model noise, or the 'distance' between the physical phenomenon and its abstract mathematical representation. In our context, this 'distance' is related to effects not taken into account by the model, such The red axis corresponds to the laser beam axis, the blue and green vectors represent the location of two different observations. In this representation, the location P i of observation i is given by its polar coordinates (R i , θ i ). (b) Illustration of the estimation accuracy for three possible choices of P 1 and P 2 along the x axis. In this graph, the coordinates of each point correspond to {P 1 , P 2 } and its color corresponds to the temperature error obtained using sensors placed at P 1 and P 2 .
as potential heterogeneity of the medium properties, sensor location and boundary conditions. Having access to ground truth temperature values, we have adopted an empirical approach which consists of comparing the estimation accuracy for a large range of Q values and evaluating the sensitivity to Q as well as its optimal value, as discussed in section III-B. We use two metrics for computing the estimation accuracy. The mean and max errors represent respectively the average and the maximal, spatial and temporal difference between the estimation and the ground truth values of the temperature. The mean error is calculated as the difference between the temperature predicted with the specific method and the measured temperature, and it is averaged over the N nodes (i.e., sensors).

F. Observations Location
Besides measurement errors which are intrinsic to sensors properties, the quality of the data-driven correction is mostly related to the number and the placement of the observations inside the heated volume [25]. We propose a systematic analysis of the effect of the observation position on the filter estimation.
To perform this analysis, we represent the domain of interest Ω using a polar coordinate system (Fig. 2). We believe it is a natural choice as the diffusion process is an axial symmetric phenomenon, arising from the heat source generated by the laser. For clarity purpose, Fig. 2 illustrates the 2D case where the red axis corresponds to the laser beam, the blue and green vectors represent the location of two different observations. In this coordinate system, the location P i of observation i is represented as To make such a study computationally feasible, we need to discretize Ω in n different locations, and define a number of observations k. This leads to the following number of possible combinations For sake of simplicity, we choose to perform the analysis on the observation location for a one-dimensional scenario. In the 1D case, the location P i of an observation can be represented as P i = (R, θ ∈ {0, π}) and can be reduced to: As an example, the case of 1 or 2 sensing needles (which embed several FBGs) placed along a line passing through the laser beam center is considered. We discretize this line in 8 possible observation locations at R = −4, −3, −2, −1, 1, 2, 3, 4, and with the laser R = 0. In this scenario, which can be later generalized to a more realistic case in which the sensors do not pass through the laser beam center, we want to study the influence of the sensors placement onto the temperature estimation. Fig. 2(b) describes our findings, where (i, j) represent the location of the two FBG sensors (one called "position of observation 1, P 1 ," the other one "position of observation 2, P 2 "). The color corresponds to the relative temperature error obtained using FBG sensors placed at P 1 and P 2 . Based on this methodology, quantitative results will be discussed in section III-A and III-B.

G. Experimental Setup
We tested the proposed approach in two simplified settings which describe LA in biological tissues (Fig. 3). At first, we used a 3% agar gel phantom mimicking the thermal properties of the biological tissue ( Fig. 3(a)) [26]. The homogeneity of the phantom makes it suitable for a preliminary analysis of our data assimilation Bayesian framework. The second setup was based on fresh ex vivo porcine liver, obtained from a local butcher and stored in the fridge before the test (Fig. 3(b)). The ablation was performed with a diode laser (LuOcean Mini 4, Lumics, Berlin, Germany, 808 nm) and the laser light was guided through a quartz optical fiber connected to a collimator, positioned perpendicularly at 5 cm from the phantom and liver surfaces. Two cases have been reproduced: 1D case and 2D case. For the 1D case, a laser power of 3.5 W was irradiated on both the phantom and liver for 60 s. An array of 25 FBG sensors embedded in a single fiber and placed on the phantom surface monitored the temperature evolving during laser irradiation. Similarly, an array of 40 FBG sensors was arranged upon the liver surface. For the 2D case, a grid of 25 x 5 FBGs (5 fibers, each embedding 25 sensors, for a total of 125 sensors), placed at 2 mm distance, measured the phantom temperature during contactless irradiation (Fig. 3(c)). An optical spectrum interrogator (Micron Optics si255, Atlanta, USA, 1 pm accuracy) was used to interrogate the sensors and collect their optical output, which is function of the measured temperature. The initial samples temperature was measured by a K-type thermocouple and resulted to be 23±1 • C in all the settings. Details about the measurement systems can be found in previous works from the same group [12], [27]. For each setting, three experiments have been repeated in the same conditions (laser power and wavelength, time).

III. RESULTS AND DISCUSSION
In the following sections, we aim at demonstrating the efficiency of our physics-based Bayesian approach, evaluating filter performances on synthetic, phantom, and liver data. For all the above scenarios, we considered both simple state and state-parameters estimation. In general, using a Bayesian approach enables retrieving a more accurate temperature distribution compared to using a purely deterministic model.

A. Synthetic Data Results
We first considered a synthetic data set defined as a 1D grid of L-connected nodes, where the temperature of each node is described by (2). In general, having a known reference scenario allows outlining some initial hypothesis about model uncertainties and their impact on estimation accuracy. In particular, we have evaluated the sensitivity of the model to errors in physical parametrization. Model parameters can be known from literature within a range of physically coherent values: D ∈ [0.05, 0.3] mm 2 /s [28], T l0 ∈ [5.13, 15.0] • C, and 1/(σ 2 ) ∈ [0.11, 0.25] mm −2 . In light of that, we forecast several deterministic simulations, each characterized by a different initialization of physical parameters, whose values have been chosen within ranges above presented. Different temperature estimations result from different initializations of D, T l 0 , and 1/σ 2 (Fig. 4). These results highlight the sensitivity of the deterministic model to the three parameters. A correct knowledge of their initial value is paramount for an efficient estimation of the actual temperature.
To address this aspect, we compared the accuracy of a simple state-estimation with the accuracy of a joint-estimation (where data assimilation is performed along with temperature estimation), together with the estimation provided by the purely deterministic model. The ground truth scenario, which is the reference for comparing the quality of our estimation, is initialized with D real = 0.05 mm 2 /s, T l0 real = 5.13 • C, and 1/σ 2 = 0.11 mm −2 . Instead, incorrect parameters are initialized as: D = 0.175 mm 2 /s, T l0 = 10.6 • C, and 1/σ 2 = 0.18 mm −2 . External observations consist of the temperature values measured at two known locations on the grid (−2 mm and 2 mm around the center of the laser beam). The Bayesian approach enables achieving a more accurate temperature distribution compared to using a purely deterministic model (Fig. 5). In particular, whenever including the incorrect parameters within the state vector (purple line), we are able to accurately retrieve the ground truth distribution (blue line), as provided in Fig. 5(a). Fig. 5(b) shows the mean error respectively obtained by the joint-estimation, the state-estimation and the deterministic model with respect to the ground truth. The time evolution of the mean error further confirms that solely when performing parameter estimation, the estimation error converges through time (purple line). Although less accurate, performing simple state-estimation (green line), allows retrieving better results than the pure deterministic model (yellow line).
As above stated, having a better knowledge on physical parameters, allows retrieving an accurate estimate of the real temperature distribution with a maximal error <1.5 • C (purple line in Fig. 5(b)). In addition, performing the joint-estimation, where we include incorrect parameters within the state vector, allowed us to rapidly retrieve their real value. As shown in Fig. 6, starting from an erroneous initialization, the filter is able to estimate the real value of the parameters after a few iterations.
In conclusion, a Bayesian framework provides in general a better estimation than a purely deterministic model affected by errors and uncertainties. In particular, a joint-estimation, where the sources of error (i.e., incorrect parameters) are included within the state vector, allows significantly increasing the accuracy of the estimation.

B. Phantom Data Results
We carried out a second series of experiences on a agar phantom mimicking soft tissues (Fig. 3(a)). The FBG sensors are used to retrieve the actual temperature distribution along the surface, which is used as ground truth (Fig. 7).
The mean standard deviation on the central sensor (11 mm) is 2.1% along the whole trend in time (Fig. 7(a)) and is 2.5% considering the spatial profile at 60 s ( Fig. 7(b)), thus proving a good repeatability across the experiments [5]. For this reason, one of the tests has been randomly selected as the ground truth Fig. 6. Parameters Estimation. We parametrize the simulation using, for each parameter, the maximal value of the parameter range. During the assimilation process, the filter is initialized with the mean value of each parameter range. By including the model parameters in the state vector, it is possible to correct their value during the assimilation process. for the validation of our framework. Through our Bayesian framework, by combining the prediction model for temperature distribution and external observations obtained from two FBG sensors, we want to retrieve an optimal estimate of the phantom temperature. Generally speaking, the accuracy of filter estimation is not absolute but depends on both the filter parameters initialization and the quality of external observations. For that, we used this second data-set to empirically evaluate the tuning of the filter parameters and external observations characteristics.
Errors associated to estimated variables, which are due to model uncertainties, are taken into account by the filter through the model error covariance matrix (Q). Since a unique criteria to initialize such parameters does not exist, we wanted to test the sensitivity of our framework to different values of Q. In general, such quantity can be defined as a squared matrix with non-null elements on the diagonal Q = diag(Q T ) ∈ R N ×N . Fig. 8 shows the sensitivity analysis of our framework to different values of the parameter Q. The temperature error, Te, is calculated as the difference between the estimated temperature values and the ground truth, being this last the temperature measured by the FBGs. The estimated temperature values are obtained by updating the filter with two measurements, and the rest of the measurements are used as reference (ground truth) to evaluate the filter performance. In our case, choosing different values of Q T ∈ [0.01, 100], we saw that the model is not sensitive to small values of Q T (Fig. 8), and different initializations within a small range provide the same results (best results achieved for Q T = 1.0). Instead, estimation accuracy strongly decreases when significantly increasing the value of Q T . This means that our state vector is affected by a certain error due to uncertainties in model parametrization, which can be efficiently taken into account by initializing Q with small values. Overestimating the error associated to the estimated vector (i.e, Q T = 10, Q T = 100) implies a deterioration on filter performances.
External observations, providing information on the current state of our estimated variables, are provided by M sensors located at different emplacements of the tissue phantom. To complete this study, we explored the influence of the location of the observations on the temperature estimation. We recall here that, from a clinical standpoint, we are limited in the number of needles that can be placed in the organ to monitor the temperature. For this reason, in the following example, we limit the number of needles/observations to two.
Following the methodology described in Section II, Fig. 9 shows that the best temperature estimation is obtained when following two heuristic rules: first, one observation should be placed as close as possible to the heating source; secondly, the 2 observations should not be at the same distance from the laser beam center. The two needles can be placed on both sides of the laser source, but also on the same side, as long as they measure a temperature gradient. Considering only two observations (k = 2), i.e., using only two sensors, we empirically estimated the optimal locations which give the best temperature estimation in our specific scenarios. Since our material is homogeneous, the temperature diffusion is isotropic, following a radial pattern centered at the laser beam spatial origins. As a consequence, in this example we limit our observation locations to a line passing through the heat source. Fig. 9 reports our results.
We can see that the empirical rule derived from the synthetic data in section III also applies. Under the constraint of using only two observation points, optimal temperature estimations are obtained when the two sensors are placed at different distances from the heat source (Fig. 9). Our approach can estimate the tissue temperature, at any two dimensional locations within the phantom. When using the parameter estimation approach with synthetic data (Fig. 9), the minimum error corresponds to 1.6 • C, and it is obtained for the locations (2,1). The maximum error results to be 5.9 • C, when both locations are at −4. Considering the phantom, potential non-homogeneity of the material may result in important errors when using only one observation point. Adding a second observation, at a different distance from the heat source, reduces the error significantly. Indeed, the minimum error corresponds to 2.2 • C, and it is obtained when the observation are placed in non-symmetrical locations (i.e., 2,1). The maximum error results to be 21.4 • C, when both locations are at 2. These results prove the effect of the location of the observations on the quality of the estimation.
A more complex configuration based on 2D temperature measurement has been adopted to assess the performances of our Bayesian approach (Fig. 10). The temperature distribution measured after 55 s from the start of the irradiation is shown in Fig. 10(a). In this case, three observations, placed in (0,−3),(2,4),(0,1), have been used to update the filter. The mean error represented in Fig. 10(b) is calculated by averaging the node temperature error within the area heated by the laser. The joint-estimation ( Fig. 10(b), purple line) allowed retrieving the best results, with a mean error <2 • C during the laser heating. These results have been verified for arbitrary and different locations of observations within the laser beam region.

C. Liver Data Results
Lastly, we evaluated the accuracy of our approach on liver data ( Fig. 3(b)). The ground truth value of the temperature is obtained from the 40 FBG sensors used within the setup. Similarly to the experiments on the agar gel phantom, also for the ex vivo  a mean standard deviation <3% was registered, in agreement with previous observation [5]. In order to perform temperature estimation through our Bayesian approach, model parameters are initialized as D = 0.175 mm 2 /s, T l0 = 10.6 • C, and 1/σ 2 = 0.18 mm −2 , whereas external observations are obtained thanks to two FBG sensors positioned at emplacement (−3,2) on the virtual grid. Lastly, filter parameters are initialized with Q T = 1.0.
The achieved results are presented in Fig. 11. Experiments performed on liver further confirm the results obtained in the synthetic scenarios. A joint-estimation (purple line), in which incorrect parameters are included within the state vector, allowed retrieving the best results, with a maximum error lower than 1.3 • C, which is maintained over the treatment time. Conversely, the state-estimation and the model provide respectively an error of 5.1 • C and 23.5 • C, after 50 s of irradiation.
We have analysed the performance of the proposed Bayesian framework also with an increased number of observations, taking advantage of the quasi-distributed measurement properties of the FBGs (Fig. 12).
For instance, when using 5 ( Fig. 12(a)) and 10 ( Fig. 12(b)) observations, both the state-estimation and the parameters estimation have improved performance (the mean error is significantly reduced in comparison with results of Fig. 5(b)). The mean error obtained with the parameter-estimation is similar in both cases, and with a maximum value of 0.5 • C; on the other hand, the performance of the proposed framework with only state-estimation improves with the number of observations (the mean error decreases from 1 • C to 0.5 • C, green line).
The proposed approach sets the basis for the implementation of a patient-specific model for procedure dosimetry and temperature estimation based on the temperature information from FBGs. When compared to thermometry based on diagnostic imaging [29], [30], sensors provide limited information about the temperature profile, but they are more interesting from a clinical point of view for their usability and cost-effectiveness. To reconstruct the temperature map from FBGs-based measurement, we propose to use an Unscented Kalman Filter combined with a heat transfer model. In our data assimilation approach, the filter simultaneously predicts the temperature and estimates the model parameters, such as the thermal diffusivity of the medium and the properties of the laser beam. This approach simultaneously furnishes a time-varying temperature profile and improves the predictive capability of the model. Once the parameter-estimation process has converged (which only takes a few seconds) we can estimate the temperature at any point of the tissue with a maximum error <1.5 • C while using only two temperature measurements in the 1D case. This makes our approach clinically applicable. In the more complex scenario of 2D measurements, the temperature error provided by the parameter-estimation framework is <2 • C after 50 s of irradiation. This last finding suggests that when the complexity of the model and of the ground truth increases, the use of additional parameters to be estimated by the filter can further aid the correct temperature estimation.
The proposed framework is able to perform the parameters estimation in real-time, which is a fundamental ability for the future intraoperative monitoring of the laser treatment. Indeed, for the temperature acquisition of 90 s with a sampling time of 1 s in 1D configuration, the filter estimates the parameters in 21 s (we run the algorithm 20 times and 21 is the rounded average running time). For the 2D case (60 s of acquisition), the filter takes 18 s to perform the parameters estimation. These excellent performances in terms of computational cost will allow realizing our main long-term objective, which is the automatic control of the laser settings [14]. Here, we leverage the benefits of joint-estimation: the predictive capability of the model is progressively improved by learning system-specific parameters from current temperature measurements. The accurate parameterization of the model is a key factor for forecasting the impact of a change in the laser power on the tissue temperature distribution. While classical state-estimation exploits observations to compensate for model uncertainties, joint-estimation brings additional benefits. At first, refining the model parametrization improves the quality of the prediction step, which in turn results in a more accurate estimation by the prediction-correction loop. Secondly, a suitable predictive model allows performing better when real-time data is temporarily missing. Lastly, if the parameter values are constants over time, their calibration at the beginning of the process may result in a model which has acceptable performance without the necessity to rely on continuous measurements. We have chosen to address the most challenging scenario, in which the number of observations is strictly limited, in order to reduce the invasiveness of the monitoring. Thus, the results of our framework are extendable to other contact thermometry techniques, in which only a few single sensors can be used (e.g., thermocouples). The developed approach is valid for all the heating modalities, after considering the specific heat source and parameters of interest (e.g., electrical properties for RFA, acoustic properties for HIFU) [31], [32]. However, it is worth highlighting that using the temperature output measured by several FBGs embedded in a needle will improve the prediction capability of the filter, with a mean error of 0.5 • C (Fig. 12). Before being applied in a clinical context, our methodology needs to be further developed in three key areas: tissue heterogeneity, temperature range and dimension. Regarding tissue heterogeneity, our method will need to be extended to include: i) supplementary heat sources and related to blood flow and related parameters to be estimated by the filer, and ii) space and temperature-dependent parameter which describes heterogeneity in thermal diffusivity and its dependence with tissue temperature. This upgrade will account for the patient specific conditions, which include the tumor properties and the individual blood perfusion characteristics. It is worth noticing that a patient-specific system including all the material properties it is still a challenge in the computer medical simulation for intraoperative guidance [33]. For this reason, the development of the Bayesian framework for real-time parameterestimation enables to adjust the model parameters during the procedure, and to provide an estimation of effective parameters that improves the model predictions when the model is oversimplified. One of the limitations of this study is the restrained temperature range on which the developed Bayesian framework has been tested. The temperatures achieved in the experiments on phantom and liver used in this work (51 • C) include the optimal temperature for the hyperthermal treatment, i.e., 43 • C. Thus, future analysis will focus on expanding the temperature range up to ablative temperatures, e.g., 60 • C. Regarding the dimension of the model, the choice of a 1D representation is a simplification intended for a first analysis of the framework validity. The process of heating around the laser applicator tip is axial symmetric, thus the 1D model captures the essence of the physical phenomenon. The material heterogeneity and experimental conditions can affect the symmetry, so this simplification may be one of the causes of the large error between the 1D model and the ground truth. However, when moving from one to two dimensions, the approach implemented for the parameter-estimation remains the same, thus we proved that the results can be generalized to higher dimensions, but still in superficial laser irradiation. With all these premises, this work poses the preliminary basis towards the development of a more complex 3D model representing the final therapeutic application with the interstitial approach [24]. In this future step, two main aspects should be considered: the exploration of supervised training to refine the temperature and spatial accuracy, and the effect of the computation time for the successful use of the Bayesian framework. To keep the assimilation process as close to real-time as possible, we will investigate methods to both speed up the numerical simulations (using for example GPU accelerators [34]) and will explore reduced data assimilation variants such as ROUKF [35].

IV. CONCLUSION
We present a fast data assimilation Bayesian framework for the real-time reconstruction of the spatiotemporal profile of the tissue temperature, evolving under laser irradiation. An Unscented Kalman Filter associates the predictions of a heat transfer model with sparse real-time temperature measurements obtained with fiber optic-based thermometers. The filter simultaneously predicts the temperature and estimates the model parameters (medium thermal diffusivity and laser beam properties). This approach provides a time-varying temperature map and improves the predictive capability of the deterministic model. Results on different scenarios (synthetic data, phantom and liver) show that the joint estimation allows retrieving an accurate estimate of the temperature distribution, with a maximum error <1.5 • C in 1D setting and <2 • C in 2D case, which represent an optimal condition for the thermal therapies control. Our approach allows also solving a key aspect of the intraoperative monitoring, i.e., the optimal placement of the sensors in the target. We have demonstrated that good temperature estimations, when using only two sensors, are obtained when one sensor is placed in proximity of the source, and the other one is in nonsymmetrical position. Further optimization of the framework in terms of three-dimensionality of the domain, accounting for tissue heterogeneity and physiological heat sources, will make the methodology and findings of our work applicable in a clinical scenario, towards the implementation of an automatic control of the thermal therapy.