Acoustic Emission Analysis: Iterative and Non-iterative Localization Algorithms

Garam Lee, 06.09.2023 (Summer semester 2023)

Introduction

Acoustic Emission Analysis (AEA) is a passive, non-destructive testing technique used to monitor the health of structures or materials. AEA does not require a signal generator; instead, sensors detect signals emitted from the damage source when the damage was created (see Figure 1.1). This technique provides us insights into the origin, development, and location of damages. It finds applications in various fields, including materials science ^[12] , quality testing and monitoring ^[8], earthquake study ^[7] , and medicine (orthopedics) ^[10] .

One crucial aspect of AEA is an accurately determining the source location of acoustic emission. One method for determination involves utilizing the time that it takes for an acoustic emission signal to travel from the source to different sensors, as presented in Figure 1.2.

For a three-dimensional case, it can be mathematically represented as the following equation ^[6] :

$//<![CDATA[ \begin{array}{l}\sqrt{(x_i-x)^2+(y_i-y)^2+(z_i-z)^2}=v(t_i-t)\end{array} //]]>$     (1.1)

where:

$//<![CDATA[ \begin{array}{l}&x, y, z&&\text{the coordinates of the source}\\ &x_i, y_i, z_i&&\text{the coordinates of the i-th sensor}\\ &v&&\text{stress wave propagation velocity}\\ &t&&\text{the origin time of the event}\\ &t_i&&\text{the arrival time at the i-th sensor}\end{array} //]]>$

The function can be solved non-iteratively or iteratively. Their primary distinctions lie in the estimation approach and computational complexity.

Non-iterative Localization Algorithms

Non-iterative algorithms are easy to apply and provide rapid estimates, because they directly solve equation (1.1) and therefore do not need to successively calculate approximations. The non-iterative algorithms have limitations in dealing with wave velocity. Although the wave velocity is a significant parameter for the accuracy of localization ^[13]  and the average wave velocities are different depending on the region in the structure, non-iterative algorithms make an assumption that the wave velocities are same for all sensors.

Inglada's method

Inglada’s method assumes the use of a minimum number of sensors ^[3] . This method solves the equation (1.1) by the linearization. By squaring both sides of the equation (1.1) and subtracting the equation of the first sensor from the rest, which yield ^[3] :

$//<![CDATA[ \begin{array}{l}a_{i}x+b_{i}y+c_{i}z=e_{i}+f_{i}t\end{array} //]]>$                                        (2.1.1)

where:

$//<![CDATA[ \begin{array}{l}&a_{i}&&\text{$x_{i}-x_{1}$}\\ &b_{i}&&\text{$y_{i}-y_{1}$}\\&c_{i}&&\text{$z_{i}-z{1}$}\\&e_{i}&&\text{$0.5(R^{2}_{i}-R^{2}_{1}-v^{2}(t_{i}-t)^{2})$}\\&f_{i}&&\text{$v^{2}(t_{i}-t_{1})^{2}$}\\&R^{2}_{i}&&\text{$x^{2}_{i}+y^{2}_{i}+z^{2}_{i}$}\\&R^{2}_{1}&&\text{$x^{2}_{1}+y^{2}_{1}+z^{2}_{1}$}\end{array} //]]>$

The linearized equation (2.1) can be now solved for ,  and  in terms of t:

$//<![CDATA[ \begin{array}{l}\begin{bmatrix} x\\y\\z \end{bmatrix} = \begin{bmatrix} a_{2} & b_{2} & c_{2}\\ a_{3} & b_{3} & c_{3}\\ a_{4} & b_{4} & c_{4} \end{bmatrix}^{-1}\left (\begin{bmatrix} e_{2}\\e_{3}\\e_{4}\end{bmatrix}+\begin{bmatrix}f_{2}\\f_{3}\\f_{4}\end{bmatrix}t\right )\end{array} //]]>$               (2.1.2)

Substituting the equation (2.1.2) in the equation (1.1) results in a second order equation of t:

$//<![CDATA[ \begin{array}{l}At^{2}+Bt+C=0\end{array} //]]>$                                                          (2.1.3)

where:

$//<![CDATA[ \begin{array}{l}&A, B, C&&\text{constant}\end{array} //]]>$

The quadratic equation (2.1.3) yields two solutions, but only one of them is valid. To determine which one is valid, we need to either compare the calculated origin time of event with the observed arrival time or to use additional arrival time of an additional sensor.

The USBM method

The USBM method, also known as the United State Bureau of Mines method, uses the linearization to solve the equation (1.1). The method requires at least one more equation than the number of unknown variables and assumes that the velocity is known. It starts by letting the equation (1.1) as distance between the i-th sensor and the damage source ^{[3, 9]} :

$//<![CDATA[ \begin{array}{l}D_{i}=\sqrt{(x_{i}-x)^{2}+(y_{i}-y)^{2}+(z{i}-z)^{2}}\end{array} //]]>$               (2.2.1)

$//<![CDATA[ \begin{array}{l}D_{i}=v(t_{i}-t)\end{array} //]]>$                                                               (2.2.2)

Since there are one more unknown variable than the number of equations, the method subtracts the distance for  from the distance for  and arrange the resulted equation ^[3] :

$//<![CDATA[ \begin{array}{l}D_{i}-D_{1}=v(t_{i}-t_{1})\end{array} //]]>$                                                  (2.2.3)

$//<![CDATA[ \begin{array}{l}d_{i}+\frac{e_{i}}{d_{i}}-\frac{2(a_{i}x+b_{i}y+c_{i}z)}{d_{i}}=-D_{1}\end{array} //]]>$                                (2.2.4)

where:

$//<![CDATA[ \begin{array}{l}&a_{i}&&\text{$x_{1}-x_{i}$}\\ &b_{i}&&\text{$y_{1}-y_{i}$}\\ &c_{i}&&\text{$z_{1}-z{i}$}\\ &d_{i}&&\text{$v(t_{i}-t_{1}$}\\ &e_{i}&&\text{$x^{2}_{1}+y^{2}_{1}+z^{2}_{1}-(x^{2}_{i}+y^{2}_{i}+z^{2}_{i})$}\end{array} //]]>$

To linearize the non-linear equation (2.2.4), the equation (2.2.4) is subtracted from the one for the rest, which yields ^[3] :

$//<![CDATA[ \begin{array}{l}f_{i,1}x+f_{i,2}y+f_{i,3}z=h_{i}+g_{i}v^{2}\end{array} //]]>$                             (2.2.5)

where:

$//<![CDATA[ \begin{array}{l}&f_{i,1}&&\text{$2\left (\frac{a_{2}}{t_{2}}-\frac{a_{i}}{t_{i}} \right )$}\\ &f_{i,2}&&\text{$2\left (\frac{b_{2}}{t_{2}}-\frac{b_{i}}{t_{i}} \right )$}\\ &f_{i,3}&&\text{$2\left (\frac{c_{2}}{t_{2}}-\frac{c_{i}}{t_{i}} \right )$}\\ &h_{i}&&\text{$t_{2}-t_{i}$}\\ &g_{i}&&\text{$\frac{e_{2}}{t_{2}}-\frac{e_{i}}{t_{i}}$}\\ &i&&\text{$3,4,...,m$}\end{array} //]]>$

In matrix form, the linear system (2.2.5) can be written as:

$//<![CDATA[ \begin{array}{l}\textbf{$Ax=B$}\end{array} //]]>$                                                                          (2.2.6)

where:

$//<![CDATA[ \begin{array}{l}\textbf{A}= \begin{bmatrix} f_{3,1} & f_{3,2} & f_{3,3}\\ \vdots & \vdots & \vdots\\ f_{m,1} & f_{m,2} & f_{m,3} \end{bmatrix},\; \textbf{x}= \begin{bmatrix} x\\ y\\ z \end{bmatrix},\; \textbf{B}= \begin{bmatrix} h_{3}+g_{3}v^{2}\\ \vdots \\ h_{m}+g_{m}v^{2} \end{bmatrix}\end{array} //]]>$

Iterative Localization Algorithms

When a single velocity model is not suitable for localization of the AE source, iterative localization algorithms are used to overcome this limitation. These iterative localization algorithms evaluate the source location by initially applying a trial solution and then by progressively minimizing error. They offer flexibility and yield more accurate results ^[2] . However, their application is complex and demands greater computational resources. The outcomes are heavily reliant on factors such as the trial function, the number of iteration and the criteria for convergence.

The error  is caused from dissimilitude between the estimated solution and the real one and must be minimized. The estimated arrival time at the i-th sensor is determined by adding the time of event to the travel time from the origin to the sensor:

$//<![CDATA[ \begin{array}{l}e=t_{i}-t_{i, estimated}\end{array} //]]>$                                                        (3.1)

$//<![CDATA[ \begin{array}{l}t_{i, estimated}=t+T\end{array} //]]>$                                                        (3.2)

where:

$//<![CDATA[ \begin{array}{l}&t&&\text{time of event}\\ &T&&\text{travel time to i-th sensor}\\ &t_{i}&&\text{observed arrival time at i-th sensor}\\ &t_{i, estimated}&&\text{theoretically estimated arrival time at i-th sensor}\end{array} //]]>$

This article introduces two iterative methods: Geiger’s method (the derivative algorithm) and the Simplex method (the curve fitting technique). Their major difference lies in how the trial solution is updated.

Geiger's method

Geiger’s method ^[5]  assesses the error by utilizing the first-order Taylor approximation and numerically seeks an optimal solution through iterative updates to the coordinate for the next step. While this method converges rapidly, it can converge to a local minimum. And the number of iterations and the final solution are highly influenced by the initial trial solution, since the solution at each step is determined by the local neighborhood.

The observed arrival time can be determined with its estimated time and first-order derivatives ^{[4, 5]} :

$//<![CDATA[ \begin{array}{l}t_{i}=\frac{\sqrt{(x_{i}-x)^{2}+y_{i}-y)^{2}+(z_{i}-z)^{2}}}{v_{i}}+t\end{array} //]]>$                                              (3.1.1)

$//<![CDATA[ \begin{array}{l}t_{i}=t_{i, estimated}+\frac{\partial f_{i}}{\partial x}\Delta x+\frac{\partial f_{i}}{\partial y}\Delta y+\frac{\partial f_{i}}{\partial z}\Delta z+\frac{\partial f_{i}}{\partial t}\Delta t\end{array} //]]>$         (3.1.2)

where  is an trial solution and f(x, y, z, t) is the right side of equation (3.1.1). The error at each sensor, also known as channel residual, is defined:

$//<![CDATA[ \begin{array}{l}e_{i}(x,y,z,t)=t_{i}-t_{i, estimated}\end{array} //]]>$                                                 (3.1.3)

$//<![CDATA[ \begin{array}{l}e_{i}=\frac{\partial f_{i}}{\partial x}\Delta x+\frac{\partial f_{i}}{\partial y}\Delta y+\frac{\partial f_{i}}{\partial z}\Delta z+\frac{\partial f_{i}}{\partial t}\Delta t\end{array} //]]>$                              (3.1.4)

In matrix form, equation (3.1.4) can be written:

$//<![CDATA[ \begin{array}{l}\textbf{$E=F\cdot\Delta q$}\end{array} //]]>$                                                                               (3.1.5)

where:

$//<![CDATA[ \begin{array}{l}\textbf{E}=\begin{bmatrix} e_{1}\\ \vdots\\ e_{i} \end{bmatrix}, \; \textbf{F}= \begin{bmatrix} \frac{\partial f_{1}}{\partial x}} & \frac{\partial f_{1}}{\partial y}} & \frac{\partial f_{1}}{\partial z}} & \frac{\partial f_{1}}{\partial t}}\\ \vdots & \vdots & \vdots & \vdots\\ \frac{\partial f_{i}}{\partial x}} & \frac{\partial f_{i}}{\partial y}} & \frac{\partial f_{i}}{\partial z}} & \frac{\partial f_{i}}{\partial t}}\\ \end{bmatrix}, \; \textbf{$\Delta q$}= \begin{bmatrix} \Delta x\\\Delta y\\\Delta z\\\Delta t \end{bmatrix}\end{array} //]]>$

The unknown vector $//<![CDATA[ \begin{array}{l}\Delta{q}\end{array} //]]>$ is calculated: 

$//<![CDATA[ \begin{array}{l}\textbf{$\Delta{q}=F^{-1}E$}\end{array} //]]>$                                                                          (3.1.6)

If the determinant of the matrix F is equal to 0 or if fewer or more than four sensors are used, the least-squares methodcan be applied to solve the problem (3.1.5):

$//<![CDATA[ \begin{array}{l}\textbf{$F^{T}E=F^{T}F\Delta{q}$}\end{array} //]]>$    or   $//<![CDATA[ \begin{array}{l}\textbf{$\Delta{q}=(F^{T}F)^{-1}F^{T}E$}\end{array} //]]>$                   (3.1.7)

The solution of equation (3.1.7) is used to estimate the coordinate for the next iteration step and the process is repeated until the error is smaller than a predefined convergence criterion:

$//<![CDATA[ \begin{array}{l}f(q)=f(q_{i})+f^{'}(q_{i})(q_{i+1}-q_{i})\end{array} //]]>$                                       (3.1.8)

$//<![CDATA[ \begin{array}{l}0=f(q_{i})+f^{'}(q_{i})(q_{i+1}-q_{i})\end{array} //]]>$                                            (3.1.9)

$//<![CDATA[ \begin{array}{l}q_{i+1}=q_{i}-\frac{f(q_{i})}{f^{'}(q_{I})}}\end{array} //]]>$                                                               (3.1.10)

Simplex method

While Geiger’s method uses differentiation to find an optimal solution, the Simplex method approaches the problem geometrically. Compared to the Geiger’s method, this approach does not require complex numerical processes such as differentiation or matrices, and there is no risk of divergence. Another advantage is that an optimal solution can be obtained in fewer iterations ^[1]  and discrete problems can be solved ^[1] .

The Simplex method assumes that the arrival time at each sensor contains error and starts building an initial simplex figure with one more vertex than the dimension of the problem space ^[11] . For example, triangle for two-dimensional space and tetrahedron for three-dimensional one. Each vertex represents error calculated through the least squares method:

$//<![CDATA[ \begin{array}{l}total\,error=\sum e^{2}_{i}\end{array} //]]>$

The vertex with the maximum error is replaced by a new point, which is calculated according to the following mechanisms: reflection, expansion, contraction, and shrinkage (See Figure 3.3).

A reflected vertex is first created and if its error has a value within the range of the former minimum and maximum errors, the reflected vertex is accepted and replace the former maximum vertex. If the reflected vertex has error less than the former minimum, an expanded vertex is generated, and its error is calculated. If the error of the expanded vertex is less than the former minimum, the expanded vertex replaces the former maximum vertex. If the error is greater than the former minimum, the reflected vertex is accepted. If the reflected vertex has an error greater than the former maximum, then a contracted vertex is created. If the contracted vertex has an error less than the former maximum, the maximum vertex is replaced by the contracted one. If the error is greater than the maximum, then the figure is shrunk by replacing all vertices accept the minimum vertex. After a new vertex is adopted, the process is repeated until the size of shape is smaller than the predefined size.

Literature

1. Caceci, Marco S., Cacheris, William P., "Fitting curves to data." Byte 9.5 (1984): 340-362.
2. El Mountassir, Othmane, et al. "Quantification of the performance of iterative and non-iterative computational methods of locating partial discharges using RF measurement techniques." K. Gessel. Wiss. Goett.Electric Power Systems Research 143 (2017): 110-120.
3. Ge, Maochen. "Analysis of source location algorithms: Part I. Overview and non-iterative methods." Journal of Acoustic Emission 21 (2003).
4. Ge, Maochen. "Analysis of source location algorithms: Part II. Iterative methods." Journal of Acoustic Emission 21.1 (2003): 29-51.
5. Geiger, L., "Herdbestimmung bei Erdbeben aus den Ankunftszeiten." Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1910 (1910): 331-349.
6. Große, Christian U., "Einführung in die Zerstörungsfreie Prüfung im Ingenieurwesen." Lehrstuhl für Zerstörungsfreie Prüfung (2022): 182-204.
7. Grosse, Christian U., et al., eds. "Acoustic emission testing: Basics for research–applications in engineering". Springer Nature, (2021).
8. Lee, Dae-Eun, et al., "Precision manufacturing process monitoring with acoustic emission. " Springer London, (2006).
9. Leighton F., Blake W., "Rock noise source location techniques." Washington, DC: US Department of the Interior Information, Bureau of Mines (1970):
10. Li, Xiaodong, et al. "Micro/nanoscale mechanical and tribological characterization of SiC for orthopedic applications." Journal of Biomedical Materials Research Part B: Applied Biomaterials: An Official Journal of The Society for Biomaterials, The Japanese Society for Biomaterials, and The Australian Society for Biomaterials and the Korean Society for Biomaterials 72.2 (2005): 353-361.
11. Nelder, John A., and Mead, Roger. "A simplex method for function minimization." The computer journal 7.4 (1965): 308-313.
12. Wevers, Martine. "Listening to the sound of materials: acoustic emission for the analysis of material behaviour." Ndt & E International 30.2 (1997): 99-106.
13. Zhou, Z.; Rui, Y.; Zhou, J.; Dong, L.; Chen, L.; Cai, X.; Cheng, R. "A new closed-form solution for acoustic emission source location in the presence of outliers." Appl. Sci. (2018), 8, 949.