• Ei tuloksia

3.2 Bayesian approximation error approach for modelling of errors

3.2.1 Uncertainties in the speed of sound

In PAT, the speed of sound in the imaged object is commonly considered as a known nuisance parameter. However, the accurate speed of soundc, that can be constant or spatially varying, is typically not available in practical applications. Therefore, some nominal value ec is used in the inverse problem, which can cause modelling errors. However, these errors can be compensated by expressing the observation model utilising the Bayesian approximation error modelling [182] in the form

pt = K(ec)p0+e+e

= K(ec)p0+n, (3.17)

where e = K(c)p0−K(ec)p0 is the approximation error describing the modelling error, or discrepancy, between the exact forward model with the accurate speed of soundcand the approximate forward model with an inexact speed of soundecand n=e+eis the total error.

When the total errorn is approximated to be mutually independent of the un-known p0, the observation model (3.17) leads to a likelihood density

π(pt|p0) =πn(pt−K(ec)p0), (3.18) whereπnis the probability density of the total errorn. Further, the posterior density can be written as

π(p0|pt)∝π(p0)πn(pt−K(ec)p0). (3.19)

Let us assume that noisee is Gaussian distributed e ∼ N(ηee), and approx-imate the modelling error e as a Gaussian e ∼ N(ηee). Thus, the total error n is a Gaussian distribution n ∼ N(ηnn), where ηn = ηe+ηe and Γn = Γee

leading to a Gaussian approximation for the likelihood. With this likelihood, a Gaussian prior, and a linear observation model, the posterior density is a Gaussian N(ηp0|ptp0|pt)given by the mean and covariance

ηp0|pt = Γp0|ptK(ec)TΓ−1n (ptηn) +Γ−1p0ηp0

(3.20)

Γp0|pt = K(ec)TΓ−1n K(ec) +Γ−1p

0

−1

. (3.21)

The realisation of the approximation erroreis unknown since its value depends on the actual unknown p0and inaccurately known nuisance parameterc. However, samples for the modelling error can be simulated utilising ’teaching’ distributions of unknowns, which are probability densities including the prior information about the unknowns, as follows. First, a set of samples {p(l)0 ,l = 1, ...,Ns}and {c(l),l = 1, ...,Ns}are drawn from the teaching distribution of the initial pressure and speed of sound, respectively. Then, samples of the approximation error are computed using

e(l)=K(c(l))p(l)0 −K(ec)p0(l). (3.22) From these samples, the meanηeand covarianceΓeof the approximation error can be estimated as

ηe= 1 Ns

Ns l=1

e(l) (3.23)

Γe = 1 Ns−1

Ns

l=1

(e(l)ηe)(e(l)ηe)T. (3.24)

4 Summary of the results

This chapter summarizes the main results of the publications constituting this thesis.

First the Bayesian approach is validated with numerical simulations in Section 4.1.

These results are based on publicationsIandII. Secondly, results of the experimen-tal validation of the approach are presented in Section 4.2. These results are based on publications IIand III. Finally, in Section 4.3, results of numerical simulations on modelling of uncertainties in the speed of sound using the approximation error modelling are presented. These results are based on publicationIV.

4.1 EVALUATING THE BAYESIAN APPROACH TO PAT WITH SIM-ULATIONS

The feasibility of the Bayesian approach to PAT was first evaluated with numerical simulations. In publicationI, simulations were done in 2D whereas in publicationII simulations were conducted in 3D. In the 2D simulations, the approach was inves-tigated utilising two different distributions of a Gaussian prior in various imaging situations including limited view measurement setups. In addition, different sensor properties such as a size and bandwidth were considered. In these simulations, the entire posterior density was determined and inspected. In addition to the poste-rior standard deviation, the reliability of the estimates was assessed by computing marginal densities and credible intervals. In the 3D simulations, the performance of the approach was demonstrated in different sensor geometries, and point estimates for the image reconstruction and uncertainty quantification were computed. In all simulations, the quantitative accuracy of the posterior mean estimates was evalu-ated by computing the relative errors of the estimates with respect to the true initial pressure distribution using

Ep0 =100%· kp0−p0,MAPk

kp0k , (4.1)

where p0is the simulated initial pressure distribution and p0,MAP is the estimated value.

4.1.1 2D

In the 2D simulations, a square domain with a side length of 10 mm was used.

The medium was assumed to be non-attenuating with a constant speed of sound 1500 m/s. The simulated initial pressure distribution is shown in Figure 4.1. Sen-sors were arranged around the domain (full view), on two adjacent sides of the domain (two side) or on one side of the domain (one side). Ideal point-like sen-sors (ideal sensen-sors), sensen-sors of width 0.5 mm (finite sized sensen-sors), and sensen-sors with limited bandwidth (bandlimited sensors) were investigated. In the case of the ban-dlimited sensors, the frequency response of the sensors was modeled as a bandpass filter with−6 dB bandwidth of 8 MHz and central frequency of 6 MHz. The simu-lated measurement data was generated by solving the wave equation (2.4) using the

Figure 4.1: The simulated (true) initial pressure distribution given in arbitrary units. The white square and circle indicate the locations where the marginal densi-ties are plotted and the dashed line indicates the location where the credible interval is plotted.

k-space time domain method and adding a 1% Gaussian noise to it. In the data sim-ulation, the domain was discretised into 300×300 square pixels (pixel side length of 33µm) and 283 temporal samples (sampling frequency 20 MHz)were used.

In the inverse problem, the mean and covariance of the posterior density were computed using (3.5) and (3.6), respectively. In addition, credible intervals were computed on a diagonal cross section of the domain (a dashed line in Figure 4.1) us-ing (3.8). Furthermore, marginal densities were calculated in two locations (a square and circle in Figure 4.1) using (3.9). In the computations, the domain was discretised into 120×120 square pixels (pixel side length of 83µm). The Ornstein-Uhlenbeck and white noise priors were used as the prior model for the initial pressure. For both priors, the mean was set as the expected mean value of the initial pressure ηp0 = 5 and the standard deviation was set asσp0 = 2.5, which means that 99.7%

of the initial pressure values were expected to be normally distributed within the range[−2.5, 12.5]. For the Ornstein-Uhlenbeck prior, the characteristic length scale l=1.25 mm was used. In addition, the accurate noise statistics were assumed.

The simulations show that the sensor geometry affects the initial pressure esti-mates and their uncertainties. In the full view sensor geometry, the inverse problem is only mildly ill-posed. Thus, the estimates of the initial pressure distribution are qualitatively similar with true one (first column of Figure 4.2). In addition, the es-timates are quantitatively accurate and the relative errors are smallest (Table 4.1).

Furthermore, the standard deviations are small indicating small uncertainty of the estimates (first column of Figure 4.2). In limited view sensor geometries, the inverse problem becomes more ill-posed. Thus, the estimates of the initial pressure dis-tribution suffer from artefacts and distortions (second and third column of Figure 4.2). The distortion of the estimates increases at larger distances from the sensors.

In addition, the quantitative accuracy of values in these distorted areas is reduced.

This results into higher relative errors (Table 4.1). However, the uncertainty of the estimates also increases, especially in the distorted areas, as the number of detection edges decreases. In all cases, the true value of the initial pressure is within the error limits obtained in the Bayesian approach (Figures 4.3 and 4.4). This indicates that the Bayesian approach can provide reliable uncertainty estimates.

The simulations also indicate that the prior has an influence on the solution of the inverse problem (Figure 4.2). In the full view sensor geometry case, the

Figure 4.2: The posterior mean (top block) and standard deviation (bottom block) obtained using the Ornstein-Uhlenbeck prior (first row) and the white noise prior (second row) in the case of the ideal sensors. The columns from left to right represent the full view (first column), two side (second column) and one side (third column) sensor geometries. The red dots in the first row images indicate the locations of the sensors.

Table 4.1: The relative errors Ep0(%) of the estimated posterior mean obtained using the full view, two side and one side sensor geometries in 2D. The ideal sen-sors (ideal), finite sized sensen-sors (finite sized), and bandlimited sensen-sors (bandlim-ited) were considered using the Ornstein-Uhlenbeck (OU) and/or white noise (WN) prior.

Ideal Finite sized Bandlimited

OU WN OU OU

Full view 13 13 23 14

Two side 15 16 33 19

One side 34 36 56 51

Figure 4.3: The true initial pressure distribution (black solid line) and the posterior mean (red dashed line) with 99.7% credible intervals (purple area) on a diagonal cross-section shown in Figure 4.1 obtained using the Ornstein-Uhlenbeck prior (first row) and the white noise prior (second row) in the case of the ideal sensors. The columns from left to right represent the full view (first column), two side (second column) and one side (third column) sensor geometries.

posterior mean estimates are almost identical and uncertainty estimates differ only slightly. However, the significance of the prior increases in the case of the limited view sensor geometries. In addition, the simulations with different noise levels in publicationIindicate that proper prior information is important in the case of noisy measurements.

In addition, the simulations show that sensor properties such as a size and band-width influence the estimates of the initial pressure and their uncertainties (Figure 4.5). In the case of the finite sized sensors, the quality of the estimates is reduced and the uncertainty is increased. This can be seen very clearly in the results ob-tained using the limited view sensor geometries (second and third column of Figure 4.5). The reduced quality cause significant increase in the relative errors (Table 4.1).

Also, the use of the bandlimited sensors reduces the accuracy of the solution, but the effect is not that significant as for the finite sized sensors.

Figure 4.4:The marginal probability densities of the posterior densities at locations denoted by a square (first row) and circle (second row) in Figure 4.1 in the case of the ideal sensors. The columns presents results obtained using the full view (first column), two side (second column), and the one side (third column) sensor geometry. Shown in the graphs are the true initial pressure (vertical black line) together with the marginal densities of the posterior density obtained using the Ornstein-Uhlenbeck prior (blue solid line) and white noise prior (red dotted line).

Figure 4.5: The posterior mean (top block) and standard deviation (bottom block) obtained using the finite sized sensors (first row) and bandlimited sensors (second row). The columns from left to right represent the full view (first column), two side (second column) and one side (third column) sensor geometries. The red dots in the images indicate the locations of the sensors.

Figure 4.6: The simulated (true) initial pressure distribution given in arbitrary units. The left image shows the contour surface that indicates the areas where the parameter has value 1 or more. The three images on the right represent maximum intensity projections along axis directionsx, y, andz. The red square indicates the location where the marginal densities are plotted.

4.1.2 3D

In the 3D simulations, the simulation domain was a cube with a side length of 10 mm and a discretisation of 304×304×304 cubic voxels (voxel side length 33µm). The speed of sound was set to 1500 m/s and attenuation was not considered. The true simulated initial pressure distribution is illustrated in Figure 4.6. Three sensor ge-ometries of same type as in the 2D simulations but extended to 3D were considered.

In the data simulation, 849 time samples (sampling frequency 60 MHz) were simu-lated using thek-space time domain method. Further, the simulated pressure signals were corrupted by a 1% Gaussian noise.

In the inverse problem, the mean of the posterior was computed by solving the system of equations (3.13). In addition, the covariance of the posterior density and marginal densities were calculated in one location (a square in Figure 4.6) using (3.16) and (3.9), respectively. For the computations, the domain was discretised into 204×204×204 cubic voxels (voxel side length 49µm). The prior model was chosen to be based on the Ornstein-Uhlenbeck prior. The mean of the prior was set to the value of the background (ηp0 = 0), the standard deviation was set toσp0 = 2 and the characteristic length scale was set to l = 0.49 mm. This means that 99.7% of the initial pressure values are expected to be in the range between −6 and 6. In addition, the noise statistics were assumed to be known accurately. For comparison, a time reversal solution was computed.

The simulations show that the Bayesian approach can be used to provide esti-mates of the initial pressure distribution also in 3D. In addition, the results of the 3D simulations follow the trend of the 2D results. That is, the most accurate esti-mates of the initial pressure are obtained using the full view sensor geometry and the accuracy of the estimates reduces in the limited view geometries (Figures 4.7 and 4.8 and Table 4.2). In addition, the uncertainty of the estimates increases as the number of detection surfaces decreases (Figure 4.9). The reduction of the accuracy is more significant in 3D. This is due to increasing ill-posedness of the inverse problem compared to 2D case in a limited view geometry. Comparing the estimates obtained using the Bayesian approach and time reversal, it can be seen that the main features of the estimates are the same (Figures 4.7 and 4.8). However, the Bayesian approach

Figure 4.7: The estimate of the initial pressure distribution obtained using the Bayesian approach in the full view (first row), two side (second row) and one side (third row) sensor geometries. The left image shows the contour surface that indi-cates the areas where the parameter has value 1 or more. In these images, the gray area indicate the locations of the sensors. The three images on the right represent maximum intensity projections along axis directionsx,y, andz.

Figure 4.8: The estimate of the initial pressure distribution obtained using the time reversal in the full view (first row), two side (second row) and one side (third row) sensor geometries. The left image shows the contour surface that indicates the areas where the parameter has value 1 or more. In these images, the gray area indicate the locations of the sensors. The three images on the right represent maximum intensity projections along axis directionsx,y, andz.

Figure 4.9:The marginal probability densities of the posterior densities at location denoted by a square in Figure 4.6. Shown in the graph are the true initial pressure (vertical black line) together with the marginal densities of the posterior density obtained using the full view (blue solid line), two side (red dotted line) and one side sensor geometries (black dashed line).

Table 4.2: The relative errorsEp0(%)of the estimated initial pressure distribution obtained using the full view, two side and one side sensor geometries in 3D. The estimates were computed using the Bayesian approach and time reversal.

Bayesian Time reversal

Full view 4 4

Two side 25 62

One side 48 80

seems to tolerate limited view artefacts better than the time reversal. The relative errors support these findings (Table 4.2).

4.2 EXPERIMENTAL VALIDATION OF THE BAYESIAN APPROACH