акустика / gan_ws_acoustical_imaging_techniques_and_applications_for_en

By substituting the total field into the Lippmann–Schwinger intergral [14], the final expression is given by

Hence, it is now possible to solve the Lippmann–Schwinger equation without omitting the scattered field, us (r ), which plays a significant role in determining the image’s resolution.

Upon arriving at equation (6.19), the next step is to create a reconstruction algorithm, which can solve the object function and generate an image. This can be a tedious process as the higher the resolution of the image, the more object functions, n(r ), will be needed. Hence, an efficient reconstruction algorithm is needed to solve the equation, which is represented in matrix form. GSPD would be used to incorporate the new model of scattering. One possible alternative is through the use of matrix formulation.

6.1.11Discrete Helmholtz Wave Equation

One other possible way of solving the Helmholtz wave equation is by first converting it into its discrete form. This discrete equation expresses each sample of the projection data as a summation of all scatterers within the object, that is,

When all the samples of each projection are taken into consideration, the above equation develops into a vector equation given by

and n (= 1, 2, . . .) is the total number of projection samples and m (= 1, 2, . . .) is the total number of scatterers.

The vector U , which contains the projection scattered field, can either be measured or simulated on a computer. For this research, the scattered field will be simulated, details of which will be covered in the next section.

Matrix A contains products of the Green’s function and the total field. Since the generated image will consist of the object in its surrounding medium, the values of r are not strictly coordinates of the scatterers; rather they are coordinates of both the scatterers and points located in the surrounding medium. The boundary of the object is taken to be a circle of radius R. This circle should be slightly bigger than the object to be examined so that all the scatterers found within this circle are considered to be the object itself. A circle is used since the interior view of the breast is circular in shape. The total field within this circle is taken to be the summation of both the scattered field (calculated previously through FA) and the incident field. The field outside the circle will only consist of the incident field.

It is necessary to solve vector N, which contains unknown values of the object function and the surrounding medium. When the number of samples and projections involved is small (e.g. 32 projections by 32 samples or less), the linear algebraic equation (6.21) can be solved conventionally by using the Gaussian elimination method or LU decomposition [20]. However, when the image resolution is expected to be 64 × 64 or above, such methods are not applicable since the size of matrix A (4096 rows by 4096 columns) is too large for storage in programs like Matlab R .

6.1.12Kaczmarz Algorithm

Many works involving equations with large matrices assume that a significant number of elements within the matrices are zero, and that sparse matrix techniques can be used. However, this is not the case here since all the elements of the Green’s function are nonzero.

One possible alternative to solve equations involving such large matrices with nonzero elements is through the use of a Kaczmarz algorithm [21]. In this algorithm, only one row of the matrix is operated at a time, thereby solving the problem of memory wastage.

This algorithm also guarantees convergence to a proper solution of the linear algebraic equation (6.21) and is able to satisfy the Helmholtz wave equation if proper discretization of all the functions is done. The Kaczmarz algorithm solves the linear equation U = A · N by considering each row of the vector equation to be representing a separate equation.

Therefore, a typical set of equations (e.g. 1 projection by 3 samples) is written as

Each of these equations actually represents a hyperplane in the three-dimensional space and the intersection of these planes represents the solution to these equations. The Kaczmarz

C

F

A

a11n1 + a12n2 = u1

n1

E

n2

B

D

Figure 6.8 Convergence of estimates in the Kaczmarz algorithm (Yu and Chan [19])

algorithm refines an initial guess by projecting its point onto another hyperplane in sequence. Every step of the iteration will bring the new point closer to the final solution.

If each row of matrix A is represented by ai (where i denotes the row number), each projection sample is represented by ui and each row of the object function is represented by nk (where k denotes the iteration number), then a better estimate of ui = ai · nk is given by

That the convergence of the Kaczmarz algorithm is guaranteed can be illustrated by the following example. For the case of hyperplanes in two-dimensional space (Figure 6.8), the initial guess n1 on line AB is first projected onto hyperplane CD: a11n1 + a12n2 = u1. This produces a new estimate (n2) on line CD. By projecting this new estimate onto the line EF, another new estimate (n3) is produced. Notice that each estimate is getting closer and closer to the solution of these equations (the intersection point of these two hyperplanes).

6.1.13Hounsfield Method

Although the Kaczmarz algorithm guarantees convergence, the speed of convergence is also a major concern here.

One factor that determines the speed of convergence is the interdependency of the row equations. If the equations generate hyperplanes that are perpendicular to each other, the correct answer can actually be calculated by one iteration. If they are parallel to each other, the speed of convergence will be very slow and more iterations are needed. In this case, the best

obtained by grouping together those equations with hyperplanes that are separated by the largest angle. However, to find a pair of equations that is least interdependent, the computational work involved is tremendous since each hyperplane will have to be measured against all other hyperplanes.

This is where the useful relationship between the scatterers and sampling points (as expressed in the GSPD) can be exploited. By rearranging the equations according to the steps shown below, a suboptimum arrangement with minimum work can be obtained. Focusing on a scatterer at a time, two rows of equations containing the scatterer’s highest and lowest probabilities in relation to a sampling point will be grouped together to form the first two rows of an imaginary matrix. By rearranging the equations this way, the chance of having two adjacent hyperplanes is minimum. This would imply that the angle between these two hyperplanes is also large and thus an improvement in the rate of convergence can be obtained. If the above procedure is not possible due to the lack of relationships, the equations can also be suboptimized by permutating the sequence of equations randomly.

Once the object function is known, the amplitude, real or imaginary values of the object function are then multiplied by a colour or grey scale factor to obtain the pixel values of the image. Based on these pixel values, a tomographic image is generated.

6.1.15Fractal Algorithm Using Frequency Domain Interpolation

Although the aforementioned method using the space domain interpolation can successfully reconstruct a tomographic image, this algorithm is too time-consuming to be used in practical medical imaging equipment. Hence, we propose that the frequency domain interpolation be implemented with the FA of the scattered field. The following section shall provide the readers with background of the algorithm used.

6.1.16Derivation of Fractal Algorithm’s Final Equation Using Frequency Domain Interpolation

The conventional diffraction tomography treatment uses the Born approximation which produces a poor approximation. In this work, instead of substituting Born’s approximation into the Lippmann–Schwinger integral equation, we substitute the total field into the Lippmann–Schwinger integral equation.

Recall that equation (6.18) is given by

60

50

40

30

20

10

Figure 6.10 Image obtained by fractal approximation (Yu and Chan [19])

in Figure 6.10, does improve the resolution compared to that using the Born approximation, which ignores the scattered field.

C1×1 . . . . . . . . . C1×64

. . . . . . . . . . . . . . . . . . . . .

Fractal Appro = . . . . . . . . . . . . . . . . . . . . .

C64×1 . . . . . . . . . C64×64

=

Equations (6.30) and (6.31) are the root mean square values of the pixels of the images formed by Born approximation and FA. As shown by the numerical values, the FA produces a smaller value with respect to the original image. This clearly shows that the FA is a better technique.

Generally, K0 is taken as 1 or 2 in the Bessel Function of Third Order (Hankel function) [24]. If K0 is changed to 300, Figure 6.12 shows that the image obtained is distorted, and the main ellipses within the Phantom have been ignored during the calculation of the Bessel function. Figure 6.13 is the mesh plot of Figure 6.12. Hence, a value of 1 or 2 is suggested for K0.

60

50

40

30

20

10

Figure 6.11 Image obtained by Born approximation (Yu and Chan [19])

6.1.18Comparison between Born and Fractal Approximations

Generally, images from both approximations suffer from severe artefacts for two primary reasons. First, the number of iterations is insignificant in bringing the initial guess of the object function closer to the final result. Second, the image resolution is based on only 64 × 64 projection samples. This is considered very low in resolution for medical or image processing standards. Artefacts could therefore be reduced by having more rounds of iterations and a higher resolution [3].

60

50

40

30

20

10

Figure 6.12 Image obtained by fractal approximation of Hankel function K0 = 300, Pcolour plot (Yu and Chan [19])

Figure 6.13 Image obtained by fractal approximation of Hankel function K0 = 300, mesh plot (Yu and Chan [19])

Although the limited rounds of iterations are not able to bring out the progressive convergence of the Kaczmarz algorithm, the differences between the Born approximation and FA are still noticeable in the images. An obvious difference is the presence of the boundary of the main ellipse in the images based on a FA. This marks a clear improvement from the Born approximation since the boundary of the object is established at such an early stage of iteration. Another improvement is the reduction of the magnitude of artefacts within the object. This is shown by the intensity of the colour, which reduces when it is within the object’s boundary. The colour intensity of images based on the Born approximation, however, remains the same [3].

6.2Nonclassical Nonlinear Acoustical Imaging

6.2.1Introduction

There is an area of nonlinear acoustical imaging known as nonclassical nonlinear acoustical imaging, which involves nonclassical contact acoustic nonlinearity (CAN) spectra. A traditional view of nonlinear ultrasonics is concerned with the classical idea of elastic wave distortion due to material nonlinearly: the waveform distortion caused by a variation in the local velocity accumulates with propagation distance and provides the progressive transition of a harmonic wave into a sawtooth or N-type wave. As a result, the spectrum acquires higher harmonics of the fundamental frequency which deliver information on the material. In the free-from-defects media, the material nonlinearity is quite low and normally only a few harmonics are observable, so that classical nonlinear nondestructive testing (NDT) is basically ‘second harmonic’ NDT.

Nonclassical nonlinear acoustical imaging is useful in the area of material characterization and NDT. In these two applications, two types of acoustic nonlinearity are of most importance, namely mesoscopic (hysteretic) nonlinearity and CAN. The origin of the mesoscopic nonlinearity is in the bond system of inclusions like cracks, grain contacts, dislocations and so on in heterogeneous materials that results in a strongly nonlinear and hysteretic stress–strain relation. The NDT applications are mainly based on the slow dynamic behaviour of the hysteretic materials: a strong mechanical impact is followed by a slow recovery. The latter is an indication of the presence of defects and is a robust methodology for a pass/nonpass evaluation in NDT [19].

The CAN [25] generally manifests in a wide class of damaged materials and is caused by the mechanical constraint between the fragments of planar defects that make their vibrations extremely nonlinear. The cracked defects therefore efficiently generate multiple ultraharmonics and support multiwave interactions. Another contribution to nonlinear vibration spectrum comes from the resonance properties of planar defects and brings the nonlinear resonance scenario with ultra-subharmonics (USB) spectra into elastic wave-defect interactions.

6.2.2Mechanisms of Harmonic Generation via CAN

Clapping Mechanism

We consider a prestressed crack with static stress σ o driven with longitudinal acoustic traction σ2, which is sufficiently strong to provide clapping of the crack interface. The clapping nonlinearity derives from the asymmetrical dynamics of the contact stiffness, which due to clapping, is higher in a compression phase than that for tensile stress when the crack is assumed to be supported only by an edge-stress.

The above behaviour of a clapping interface can be described approximately by the following piecewise stress (σ )–strain (ε) relation [26]:

where H (ε) is the Heaviside unit step function; C = [C − (dσ /dε)ε>0], and C is the contact material (linear) stiffness.

The bimodular prestressed contact driven by a harmonic acoustic strain ε(t ) = ε0 cos νt is similar to a ‘mechanical diode’ and results in a pulse-type modulation of its stiffness C(t ) (Figure 6.14). It also provides an unconventional nonlinear waveform distortion: a halfperiod rectified output (Figure 6.14) instead of the sawtooth-like pulses in classical materials. Since C(t ) is a pulse-type periodic function of the driving frequency ν (Figure 6.14, right), the nonlinear part of the spectrum induced in the damaged area σ NL (t ) = C (t ) · ε (t ) contains a number of its higher harmonics nν (both odd and even orders) whose amplitudes are modulated by the sinc-envelope function [26]:

An = C τ ε0 [sin c ((n + 1) τ ) − 2 cos (π τ ) sin c (n τ ) + sin c ((n − 1) τ )]

(6.33)

where τ = τ /T (τ = (T /π ) arccos (εo/ε0 )) is the normalized modulation pulse length.