Solutions of inverse problems with potential application for breast tumour detection using microwave measurements

This paper presents one-dimensional and two-dimensional microwave inverse computing methods to detect an internal object using measurements based on a signal applied from the surface of the host material. The modelling of our application system has been aimed towards the in vivo detection of a breast tumour, in particular, and to enable the calculation of the tumour size and its distance from the surface of the breast. However, our approach is also applicable for more general foreign object identification. Complex backscattered electromagnetic waves characterise the relations of the internal properties of the host material. Forward and backscattered signals are used to calculate the impedance and reflection coefficients as a function of the applied microwave frequency. In the study of onedimensional modelling, we discuss the approach to identifying a foreign object hidden inside the host material and we present a method for computing the distance to the object from the surface of the host. Subsequently, a cylindrical coordinate system is used for two-dimensional modelling. A method to compute the size of the object (up to one millimetre in radius) is discussed. Computation of unknown electrical and non-electrical parameters using front-end microwave application is challenging but it is feasible.


Introduction
Breast cancer is a non-skin malignancy in women and is the most prevalent cause of female cancer mortality [1]. In vivo methods of early-stage breast tumour detection make a significant contribution towards the likelihood of successful treatment. Microwave technology can be used as a non-invasive method to detect breast tumours in early stages. Finding accurate and stable solutions to inverse problems is the most important task in object detection using microwave measurements. The microwave signal we use in this application is harmless, because it has a low power and therefore does not damage the normal cells. Screening mammography is the most effective method used at present for breast cancer detection but it suffers from a number of drawbacks such as: high false-positive and falsenegative rates, a possible risk factor, discomfort for the examinee and their difficulty in tolerating breast compression [2,3]. In our approach, the microwave signal applied from the surface of the breast skin requires a minimal compression of the breast for accurate measurements.
Our earlier research for internal property measurement of dairy products and food samples have shown that microwave imaging is feasible using a dielectric permittivity profile obtained from a suitable measuring system [4,5]. Further, ex vivo measurements taken using the Keam Holdem VE2 analyser [6] have shown that a tumour has a significant difference in complex dielectric permittivity to that of healthy breast tissue. As the malignant tumours have increased protein hydration, they have a significant contrast in dielectric properties with normal breast tissue [7,8]. Results of the investigation performed using microwave radar technology [9] also illustrate the opportunities for active microwave sensing in the breast. The contrast between the normal and malignant breast tissue can provide a significant difference in the backscattered microwave energy [10] but, for a time-domain approach using the ultra wideband radar technique, it can be challenging to obtain robust signal processing methods and appropriate hardware is needed. In the study of confocal microwave imaging for breast imaging [7], the breast is modelled with planar and cylindrical configurations and methods are developed to detect and localise tumours in three dimensions. In this approach, a finite difference time domain method is used to compute the backscattered data. Another development uses optical tomography to obtain medical imaging for diagnosis of possible breast tumours [11]. In this approach, the solutions to the time-dependent diffusion equation are found using data produced from only one source but with many detectors. However, these methods too may have practical limitations when used for breast tumour detection.
We propose non-invasive methods that can be used to identify a very small tumour inside the breast. Our approach analyses the behaviour of the microwave signal and computes the distance from the surface of the host material. The type of signal used in this application is a uniform plane wave which penetrates through the non-homogeneous internal structure. In practice, there are losses due to finite conductivity and lossy dielectric but these are usually very small and can be neglected [12,13]. The behaviour of the signal with different material properties have led to general equations which can be obtained from the well-known theory of electromagnetic wave propagation [13,14]. Those equations contain information on the electrical and magnetic properties of the internal structure and can be used to develop algorithms to compute the unknown parameters of the internal object. The front-end microwave measurement provides us with the required information which is needed to identify and then to compute these parameters. We discuss a simple and cost effective approach for detecting a tumour in the first instance using microwave measurements. With further development, this method has a potential for breast screening. Earlier, we provided in Ref. [15] evidence that the method works in the two-dimensional situation, and this is extended further here. The derivation is given here for completeness.
The X-ray technique is very common as a diagnostic tool but it needs the signal to penetrate fully through the breast and the quality of the received signal is dependent upon breast compression. In the method, we propose, breast compression is not required because the information about the tumour, if any, is found using the difference between forward and reflected signals measured at the same antenna position. We use a continuous microwave signal rather than a discrete signal as used for the conventional X-ray approach.
The basic model for the microwave measurement system is shown in figure 1. The measurement system provides a microwave signal to the antenna system and this signal is then excited into the host material. The backscattered signal from the internal structure of the host is received by the same antenna system and is returned to the measurement system for analysis. The measured data is then processed using the reconstruction algorithms. In both transmitted and received signals, amplitude and phase changes are expected and will be measured a number of times at different frequencies. There are two main parts to our study. One is to solve the inverse problem using the reflection coefficient measurements based on the multi-layered plane wave reflection. The other approach is to consider the internal object as a 'wave scatterer' and solve the two-dimensional inverse problem in cylindrical coordinates to compute the unknowns. We have begun with simple canonical geometries in order to illustrate the general approach. The one-dimensional study is straightforward but helps us to understand the practical difficulties with accuracy when working with plane wave measurements. It is important because the computed results of the forward and inverse algorithms provide an insight into the subsequent two-dimensional and three-dimensional cases. The two-dimensional study lays the groundwork for the practical computation of the inverse method using a simple microwave measurement system. Although the forward problem is simple, the inverse problem in a two-dimensional application may have practical complications due to the complexity of multiple scattering [16], multi-path and diffraction effects in non-homogeneous internal structure of the host material.
The frequency of the microwave signal must have a constant value during the measurement time but it may be changed to another value for subsequent measurements. The reflection coefficient [17] at the front surface of the model for any given profile is given by where f represents the frequency of the microwave signal, Z in ( f) is the complex electrical impedance into the surface of the breast skin and Z 0 is the complex electrical impedance of the measurement system. In order to develop the model, we initially develop the forward and reflected wave theory using the formulation appropriate to the geometry (one-and two-dimensional situations) considered. We then, with field equations, use the theoretical results to find the dimensions of an object from the effect of the scattered reflected wave by the familiar inverse problem technique. This is explained in following sections.
In the first case considered, the internal structure of the breast is modelled as thin layers (described in section 2). A basic analytical system is constructed for forward and inverse computations to find front-end impedance and the thickness of the layers. In section 3, a more realistic approach for microwave scattering is developed using a cylindrical coordinate system. An inverse algorithm for calculating an unknown tumour's size and its location is presented and the error and stability is discussed. For a general outline of inverse techniques see Ref. [18].
There are more sophisticated models for inverse estimation using finite bases, that is, pixels representing the material properties, which may be applicable in microwave tomography [19,20]. Our approach is to find the size and the location of the object using microwave measurements. As the object is relatively small in size compared to the wave length of the signal, a fair approximation is for the object to be a circular in shape. We seek a simple and cost effective method for a dielectric reconstruction of an internal object using the frequency domain with phase and amplitude measurements.
In this study, we first use the numerically generated data to test our inverse algorithm. Then, the process is extended to error and stability analysis. By doing this, we investigate the possible range of guess values which can be used when calculating the unknowns from measurement results.
Our computer-programmed algorithms have been verified with laboratory microwave experiments. In this experimental application, a conducting cylinder was placed at different distances in front of an antenna. The amplitude and phase of the reflection coefficient (the ratio of incoming and outgoing microwave signals in complex form) were measured using a network analyser. The amplitude and phase of the measured reflection coefficients with respect to the position of the cylinder have been plotted in figures 2 and 3, respectively. Also, the calculated reflection coefficients using our forward equation have been plotted in the same graph for easy comparison. Overall, there is a good agreement between calculated model predictions and experimental data (figures 2 and 3). More details of the experimental measurements and the calculation of unknowns using the measured data are explained in Ref. [21].

One-dimensional study
In the following sections, we analyse the time and space dependence of the microwave signal within the application model. The plane wave reflection model is shown in figure 4. The internal structure of the host material is represented using a number of regions, each of which is homogeneous. In particular, we assume that the electrical properties are constant over each of these regions. For the one-dimensional model, the regions are a number of thin rectangular layers which have been cascaded to form the host material. The layers inside the model are specified with individual material properties. These are permittivity 1, permeability m and conductivity s and they characterise the media with electric flux, magnetic flux and the electric current, respectively. When the microwave signal is applied from the front, it penetrates through the layers and, if the properties of any two layers differ from each other, it reflects back from the boundary between them. Similarly, looking from the electrical view point, it can be observed that each of these layers must have individual impedances in the presence of uniform plane waves. In order to find the reflection from the surface of the host, it is necessary to perform a series of impedance transformations at the layer boundaries. The impedance transformation towards the front end can be seen as a belt with n cascaded strips, as shown in figure 4.

Impedance transformation
The front-end impedances of the layers are indicated as Z (·) (looking from the front) and the first and last layer impedances are taken as Z in and Z nþ1 , respectively. We consider the host internal structure to be lossless (s ¼ 0) for the electromagnetic waves and therefore the wave propagation depends only upon the complex value of the propagation constant [22]. For simplicity, in the reminder of this paper the magnetic permeability m is assumed to be unity, but this is not a restriction for this application. The recursive equation to find the electrical impedance [14] at the front of the nth layer of the model, that has a width of d n , is where h n ( f) is the intrinsic or characteristic impedance of the nth layer given by h 0 is the intrinsic impedance in a vacuum and is approximately equal to 377 Ohm and j ¼ ffiffiffiffiffiffi ffi 21 p . The 1 n ( f) is the relative permittivity of the nth layer and b n ( f) represents the phase constant of the nth layer given by where c is the velocity of light given as approximately 3 £ 10 8 m/s and f is the frequency of the microwave signal. Given the characteristic impedance and the propagation constant of any layer and the load impedance of the succeeding layer, the front-end impedance of that layer can be calculated using equation (2). Starting with the known impedance of the last layer (the deepest layer of the model), the layer impedances can be computed from layer to layer up to the surface of the front layer. Finally, this result is used to compute the reflection coefficient using equation (1).
We consider three layers having thicknesses of d n21 , d n and d nþ1 . The front-end impedance of the (n 2 1)th layer at frequency f i can be found by substituting Z n in the equation that is obtained for Z n21 using equation (2). Then at frequency f i , where Z i , b i and h i represent the front-end impedance, phase constant and characteristic impedance, respectively (at the frequency f i ). Equation (5) calculates the front-end impedance without knowing the front-end impedance Z n of the middle layer. Similarly this procedure can be continued up to the first layer of our model to calculate the front-end impedance, Z in ( f). The final equation of this process would be a large and complicated equation with a number of unknowns. We can obtain i equations for i different frequencies, and these equations can be used to find the i unknowns.
If, instead, the reflection coefficients are found from measurement, we can find the frontend impedance of the host material and calculate alternative unknown quantities, for example, the distances (as explained below). Again using the equations (2) and (5), we may obtain two different equations for two frequencies (i ¼ 1 and i ¼ 2) to find two unknowns in the Z n and Z n21 layers.

Layer thickness calculation
Suppose that we know the front-end impedance, Z in ( f), by practical measurement of the host material using the microwave antenna system. Then, the next task is to find the distance to the scattering object from the surface of the host. The total distance is the sum of the individual widths of each layer. This is calculated using inverse equations derived from equation (5).
If we take the three layers, the distance to the (n 2 1)th layer from the (n þ 1)th layer can be obtained from where b n and b n21 can be calculated using equation (4). By following the above procedure, we can find a similar equation for even more than three layers. In equation (6), we have two unknowns, d n21 and d n . For simplicity, we assume Z nþ1 ¼ 0. Then using equation (6), we obtain a general equation as where f i is the microwave frequency applied at each measurement, and i ¼ 1,2,. . ., m.
With more than one unknown, it is not possible to obtain a direct solution using a single equation, therefore we use m equations produced from m trials each applying a different frequency into the system. In fact, for equation (7) to have a solution Z n21 has to be purely imaginary, i.e. a complex number with real part equal to zero. In order to then find the unknowns, we use an algorithm based on Newton's iterative method [23].
Let our unknowns be the thicknesses of the layers. Taking m ¼ n frequencies then, the set of equations are: We form the n £ n Jacobian matrix of the above system of equations where J k;l ¼ ›F k =›d l and k, l ¼ 1,2,. . .,n.
Using the initial guess for X ð0Þ ¼ ðd ð0Þ 1 ; d ð0Þ 2 ; . . .d ð0Þ n Þ T , we carry out a multidimensional Newton's method [24] to search for the solution to F(X) ¼ 0. The solution of the above system often needs several iterations, the number of iterations dependent mainly upon the number of unknowns and the value of the initial guess.  (2), the front-end impedance of 10 layers was calculated recursively. Starting with the last layer, Z nþ1 , the calculation is carried out from layer to layer up to the first layer, Z 1 ¼ Z in . The simulation results are shown in figure 5. The plot (i) is for ten layers each having electrical properties: m ¼ 1, s ¼ 0 and 1 ¼ 10.
We assume the medium is lossless. In all the layers, the permittivity 1 was equated to 10 which is assumed to be similar to that of the normal breast tissue [2,3]. However, these values of the electrical properties vary with the frequency in use at the time of measurements. The plot (ii) is similar except that the last layer's permittivity is set to be similar to that of a breast tumour (1 10 ¼ 50; There is a significant difference in front-end impedance when the permittivity of the last layer is 50 rather than 10 as in the other layers. With our detection system, there should be a high probability of identifying an internal object with significantly different material properties to its surroundings.

Distance calculation.
Finding the layer thicknesses is important because their arithmetic addition will give the distance from the surface to the layer boundary of interest. For a test example, we used our algorithm to compute the thicknesses of two layers (n ¼ 2) within our model. Using equation (5), we first calculate the front-end impedances of the (n 2 1)th layer for two different frequencies f 1 ¼ 2 GHz and f 2 ¼ 2.2 GHz. The thickness of the second and first layers are taken to be 0.004 and 0.002 m, respectively. As the f 1 and f 2 are close in value, we can set 1 i ð f 1 Þ ¼ 1 i ð f 2 Þ ¼ 10 for i ¼ 1, 2. From the values of Z 1,1 and Z 2,1 at the two different frequencies, we may estimate the values d 1 and d 2 using two equations of the form of equation (7). That is where The solutions for d 1 and d 2 have been computed using Newton's method in the mathematical software package MATLAB and are plotted in figure 6 (Newton's method is discussed further in the two-dimensional study). The two graphs show that the approximations to d 1 and d 2 in F 1 and F 2 rapidly approach the exact values of

Two-dimensional study
Here, we consider the internal object to be a conducting circular cylinder with a radius in millimetres. The microwave signal application system is the same as figure 1, but the model of the host material is different from the one-dimensional case. The two-dimensional model with cylindrical coordinates is shown in figure 7. Here, we analyse the propagation of the microwave signal and the scattering effect at the cylinder boundary (the circle with radius a in figure 7). We only consider the boundary condition at the cylinder in the model because both forward and reflected signals are measured simultaneously at the antenna front-end. The antenna is connected to the measuring system (figure 1) where the reflection coefficient is calculated using the forward and backward signals at the antenna A 1 . The scattered wave for another skin position (e.g. at A 2 and A 3 ) can be found by the same procedure.
We consider this to be a single scattering problem and construct the electromagnetic wave equations for the forward and scattered fields using the solutions to the scalar Helmholtz equation [25]. Our discussion of the two-dimensional model is restricted to the conducting cylinder. The extension to the case of a non-conducting cylinder having dielectric with parameters m and 1 similar to those of a breast tumour can be developed as in [26 -28]. In this latter case, the boundary conditions of the normal and scattering fields on the cylinder are based on the wave impedance that we discussed in the one-dimensional model.

The forward problem in the two-dimensional study
A plane wave incident upon the host material can be expressed in terms of cylindrical waves [25]. The incident wave at the host material is z-polarised and travelling in the x direction as shown in figure 5. The forward incident wave of frequency f is where p is the radial distance from the centre of the cylinder, f is the angle with respect to the x direction and k is the wave number of the medium given by As the wave is finite at the origin and periodic in f, of period 2p, equation (12) becomes where J n is the Bessel function [29] of the first kind. Due to scattering at the cylinder boundaries, for the outward-travelling waves, the scattered field is where H ð2Þ n is the Hankel function [25,29] of second kind. The total field is the sum of the incident and scattered field, that is, Since we have zero electrical conductivity in the host material and zero field intensity inside the conducting cylinder, then we have zero electric field components on the outer surface of the cylinder by continuity of the field [13,25]. Considering the boundary condition on the cylinder E z ¼ 0, the total field at the point o given by equation (16) can be re-written using equations (14) and (15).
Equation (17) The computational cost of the inverse process depends on the number of arithmetic operations required to calculate (18) accurately. It can be reduced by combining terms with positive and negative values of n. The new equation is Equation (19) is the forward equation we use to find the unknowns for the two-dimensional model. Using a similar approach to that used in the one-dimensional case, E z ( f) can be measured and then equation (19) may be numerically inverted in order to determine a and d.

The inverse problem in the two-dimensional study
In the forward problem a and d are known at a single location of the cylinder and the amplitude values of the scattering field at each frequency have to be computed as explained elsewhere [15,28]. In the inverse problem, the field quantity at the receiver is known by measurement and a and d are the unknowns. In the practical situation, we can measure E z ( f) for two frequencies f 1 ¼ 2.0 GHz and f 2 ¼ 2.2 GHz. Then, we compute the unknowns, a and d, using the inverse method. Similarly to the one-dimensional case, the general equation has two constituent equations of the form where we now use the subscripts 1 and 2 to indicate two different frequencies. In full, where c 1 is the right hand side of the forward equation when the frequency is f 1 , and c 2 is the right-hand side of the forward equation when the frequency is f 2 . Here, E 1 and E 2 are the two field components to be obtained from microwave measurements. In this analytical study, we have used reasonable values for a and d (a ¼ 0.002 m and d ¼ 0.04 m) to calculate these two values. In this study, we used the electrical parameters s ¼ 0, m ¼ 1 and 1 ¼ 10 within the host material (outside the cylinder), and E 0 is considered to be equal to unity. As there are Bessel and Hankel functions inside the summation of the equation (19), it is necessary to handle this equation carefully to obtain the correct answer for E z ( f). Calculations of the field vectors and the subsequent verification using the experimental measurements have been discussed in Ref. [21]. Solution of the inverse problem can be used to find the size and the location of the internal object and the routine is summarised as follows. We start from an initial approximation to our unknowns, X (0) ¼ (a 0 , d 0 ) T . Subsequent improvements to this approximation X (N) ¼ (a N , d N ) T , N ¼ 1,2,. . ., are obtained by the following steps: (1) Compute the field component based on the incident and the scattered field. The circular boundary of the cylinder is the cause of the potential scatter and the angle f determines the receiver location for the measuring system. The wave number k is frequency dependent and there exists a field component for each frequency. (2) Form the difference of the field vectors by subtracting the calculated field components from the measured electric fields. (3) Construct the n £ n Jacobian matrix J l,m , l, m ¼ 1, 2,. . ., n, which is necessary in Newton's method to find the minimum difference in (2) above. (4) Obtain the correction vector and form the implicit function for the vector X (N) to update the computed values of unknowns in the vector X (N21) [23,24]. (5) Repeat the steps 1 -4 until the vector X (N) satisfies some suitable stopping criterion.
Following the above procedure, the two unknowns were computed. The result is shown in figure 8.   DE 2 converge towards zero as a and d approach 0.002 and 0.04 m, respectively (approximately after 12 iterations). When we determine E z ( f), we use equation (19), truncating the series at successively larger values of n until we obtain a steady value for E z ( f). The value of n will depend on the values of frequency, distance d, radius a and the electrical parameters of the medium. Care must be taken when testing for convergence of equation (19) since the solution oscillates with respect to n [28]. The accuracy of the initial guess can also be a limitation.
We carried out a number of simulations to test our inverse algorithm for convergence. First, using the forward method, the values of E z ( f 1 ) and E z ( f 2 ) were calculated for the values of a and d equal to 0.002 and 0.04 m respectively. Those results were used in our inverse algorithm to test for the convergence using a range of initial (guess) values of a and d. The selected range starts from a ¼ 0.001 m and d ¼ 0.03 m and changed by 0.0002 and 0.002 up to a ¼ 0.003 m and d ¼ 0.05 m respectively. Every pair of initial values was tested separately, and the number of iterations, N, required for convergence is counted.
The results are displayed in figure 10. The x and y axes represent the initial values of a and d, respectively. Every grid point corresponds to a pair of a and d initial values. The number of iterations (N) required for convergence is shown at the grid point. When the initial values are far away from the true values we require a large number of iterations. Furthermore, there exists a range beyond which we cannot expect any accuracy in the convergence. In our test, a and d could vary by up to^50 and^22.5% from their actual values, respectively, and still converge. In our example, we found that a ¼ 0.0010:0.0030 and d ¼ 0.030:0.048 is a safe range for convergence within a reasonable number of iterations. We suggest that in general 20 iterations are needed in order to determine whether the process is within a safe range before restarting the iteration process with alternative starting values for a and d. It is important to display the result with double precision in order to identify the exact solution.
We carried out an error analysis to study the robustness of the inverse computation method with respect to measurement errors in the forward system. We added errors into the E z ( f 1 ) and E z ( f 2 ) simulated values to find the corresponding errors in a and d values. Those values are tabulated in table 1. The percentage error in a is quite large. We should note that our original value of a is small compared to d (d is 20 times larger than a). However in general, d is less sensitive to measurement errors than a.

Discussion and conclusions
The approach used in this study has demonstrated the feasibility of calculating the unknown size and position of a tumour inside the breast. Our inverse method is relatively simple when compared to other approaches, but has the properties of fast convergence capability and computational simplicity.
In the one-dimensional model, a layer of high permittivity equivalent to that of a breast tumour is shown to provide a significant difference in front-end impedance. Similarly, the algorithm developed for three layers indicates that the determination of the distance from the breast skin surface to the cancer is computationally feasible using in vivo microwave measurements.
As we have previously explained in Ref. [15], the two-dimensional study of the microwave inverse computing method has proved capable of estimating the unknowns. In the twodimensional model, we have taken the object boundary to be conducting and so the reflection coefficient at this boundary is equal to unity. For a non-conducting cylinder, the reflection coefficient would change relative to the electrical properties on the two sides of the boundary. For this, the equation (19) is modified to include the effect of the reflection at the nonconducting cylinder boundary. Using the wave impedances of the two regions, the reflection coefficient can be calculated at the boundary. Parameters for the properties of both regions appear inside the new equation [28,30].
Inverse problems are notorious for being ill-posed yielding (in this case) non-unique solutions, see Hadamard [31]. This is Hadamard's original paper describing the concept of well-and ill-posedness. This will be reflected by the Jacobian matrix in the linearisation process described earlier having a high condition number and being nearly singular. This can occur here and the system of equations can have non-unique values for a and d. We took care within the computational scheme to ensure that the outputs gave unique answers for these unknowns, which were both feasible (being both real and positive) and of the appropriate sizes. This worked at the practical level. However, further work is needed to finally characterise the degree to which this problem is indeed ill-posed.
Reflected waves may be received at the antenna due to other internal scattering mechanisms. It is clear that the amplitude of the received reflected signal at the antenna due to the tumour is significantly higher than the majority of these [6]. Therefore, the measured data with and without a tumour should be easily distinguished despite these other effects. However, the reflection from the chest wall may have to be considered separately, although this effect can be included in the forward equation [21,28].
In practice, some form of calibration could be performed to reduce the influence of measurement error. This could consist, for example, of normalising the measurement data with measurements taken where it is known that there are no scattering objects present. Once the Table 1. Result of the error analysis (as in Ref. [15]).

Measurement error (%)
Error in a (%) Error in d (%) calibration procedure is completed then analysis of the measured data has a potential first to identify the presence of the tumour and, if present, then to find its size and location inside the breast. This approach could be extended to use more than two frequencies in an over-determined system. In this study, we tested and studied the computational stability of the inverse algorithm with the minimum requirement for different frequencies. This has been shown to be adequate; however, one could use more than two frequencies, and subsequently more than two equations, to find the two unknowns, and this might be more accurate. There will be practical limits on the number of different frequencies which can be used as appropriate frequencies to provide a high contrast of the measured data with and without the object. This frequency selection procedure is in Ref. [28]. The best range of microwave frequencies depends upon a number of factors in the type of application [4 -6]. We used 2.0 and 2.2 GHz frequencies for this study.
Future work will incorporate three-dimensional modelling with a similar approach. This will provide a means to compute the size and position of the breast tumour accurately. The algorithms developed at this stage will be further validated in the laboratory.