Abstract
Electromagnetic nondestructive evaluation of underground targets is of great significance for the safety of urban construction. Based on the accurate and efficient simulation of scattering, we can detect the underground targets successfully. As one of the most popular numerical methods in electromagnetics, surface integral equations solved by method of moments (MoM) are used to simulate the scattering from underground targets in this paper. The integral equation is discretized by RWG basis and Galerkin testing. Multilevel fast multipole algorithm (MLFMA) is used to decrease the computation complexity and memory cost. However, the octree used in MLFMA is not applied for rough surfaces and targets together; both the surface and target need to construct octree separately. Since the combination of MLFMA and ACA can build a more efficient method to compute scattering from underground targets, adaptive cross approximation (ACA) is used to compress the impedance matrix instead of MLFMA for the coupling action between the rough surface and target. That is to say that, when calculating the scattering of two targets, target self-interaction is suitable for MLFMA calculation and the coupling between targets is approximated by ACA. Numerical results demonstrate the accuracy and efficiency of our proposed method.
1. Introduction
Because the underground targets such as pipelines are important infrastructures in modern cities, the nondestruction evaluation has become an indispensable prerequisite for preventing the damage of these underground targets for construction. Obviously, it is necessary and meaningful to research on nondestructive evaluation technique based on scattering of electromagnetic wave for underground targets. At present, there are three kinds of numerical methods for electromagnetic simulation: finite element method (FEM) [1], finite difference time domain (FDTD) [2] based on differential equation, and method of moments (MoM) [3] based on integral equation. The calculation of the target scattering mainly uses the MoM of the integral equation, which has high computational precision. However, the computational complexity of the method of moments is for direct solution while the complexity for the iterative solution is , where is the number of unknowns. Accordingly, the high computational complexity posed a severe challenge to computer resources, such as storage and CPU. Thus, a wide range of fast algorithms is needed to be developed.
According to the principle of mathematics, the fast algorithm can be divided into three categories: the first one is based on fast Fourier transform (FFT) method [4, 5], the second is based on the fast multipole method (FMM) [6], and the last one is a class of pure algebraic methods based on low-rank matrix compression. Among them, the method of matrix compression realized by low-rank characteristic mainly uses the decomposition of matrix, which can save memory and accelerate the iteration of equations, such as singular-value decomposition (SVD) [7], QR decomposition [8], LU decomposition [9], and adaptive cross approximation (ACA) decomposition [10–13]. These pure algebraic matrix compression methods are applicable and maneuverable, which are widely used in the method of moments to accelerate the iteration of equations. In 2014, Li et al. proposed a method of compressing the impedance matrix of the method of moments (MoM) aiming at modeling of high-fidelity multiscale structures at low to the moderate frequencies, which can solve the modeling of multiscale structures efficiently [14]. Following that, a novel doubly hierarchical MOM for the analysis of large and multiscale structures is investigated by Li et al. that is effective for large structures with strongly nonuniform discretization [15].
The multilevel fast multipole algorithm (MLFMA) is one of the most attractive integral equation fast algorithms [16–18] of whose advantages are high precision and efficiency. Thus, it is widely used in electromagnetic scattering analysis of various complex targets. In the calculation of composite targets with metals and media, different media targets need to establish different octree trees. Due to the limitations of the Green function expansion, MLFMA becomes more difficult to calculate in this case. Therefore, other types of fast algorithms are further considered. Among the many matrix decomposition algorithms, ACA decomposition is widely used because of its unique advantages, such as only needing to know a part of the matrix information. Furthermore, ACA algorithm is also suitable for low-rank matrix compression. In [13], the use of ACA for multilevel tree structures has been illustrated to improve matrix iteration and storage. Because of the particular peculiarity of the ACA method, it can be used to solve this problem. The interaction of the two targets can be divided into the near field and far field based on the distance. That is, when the distance between the two targets is less than the target size, we assume that some of the basis functions interact as near field and the others interact as far field. Therefore, the octree can be used to divide the tree structure. According to the distance between the finest level groups of the two octets, we can judge that the two groups are the far field or near field. The near field is directly calculated using the method of moments, while the far field usually with low-rank features is using ACA acceleration method which would have a good compression effect.
2. Computation of Scattering of Targets under Rough Surfaces
2.1. Self-Action Accelerated by MLFMA
From the Maxwell’s equation, the surface integral equation (SIE) [19] can be established according to the equivalence principle and the boundary conditions of the electromagnetic field. For the perfect electrical conductor (PEC), the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) can be deduced by the boundary conditions of the electric field [3] and magnetic field [20], respectively. For the dielectric target, the well-known PMCHWT integral equation of the two-region medium surface can be established [21].
The surface integral equations can be discretized by method of moments, and the linear equations are built as follows: where is the impedance matrix, is the coefficient of RWG basis, and is the incident field.
Taking only one underground target for example, impedance matrix includes four parts: self-action of rough surface, self-action of the target, and coupling action between the surface and target, respectively, as written by 1. where diagonal parts and represent the self-action of the rough surface and the underground target, respectively, and the off-diagonal parts and denote the coupling actions between the rough surface and the target. Since each block of the matrix is completely independent, they can be calculated and stored separately. Each should be computed by the method of moments which makes the whole process inefficient with high complexity of .
In order to solve this problem, multilevel fast multipole method is used to accelerate the computation of impedance matrix elements and matrix vector product (MVP). However, when the distance between the surface and target is large, traditional MLFMA will generate many levels and lead to much computation time. Different from traditional MLFMA, our method uses more than one octree structure. A private octree structure is built for each surface or each target, which will lead to lesser levels. Accordingly, computation complexity of impedance matrix of self-action of the surface (target) will be , where represents the number of unknowns of the surface (target). And the memory cost for near field impedance matrix of self-action will be , either. Since off-diagonal matrices are computed by MoM, it is still time-consuming and a new technique is required to reduce the complexity.
2.2. Coupling Action Accelerated by ACA
For impedance matrices of coupling action between the surface and target, parts of them are near field and other parts are far field. Octrees can divide each off-diagonal impedance matrix into many submatrices. Each submatrix belongs to either the near field part or far field part. Far field submatrix has low-rank property and can be further compressed to save time and memory.
As a pure algebraic algorithm, adaptive cross approximation [10, 11] method can divide a low-rank matrix into two submatrices. As illustrated in Figure 1, it is obvious that the leading dimension of each submatrix equals the rank of the original matrix approximately. Compared with MLFMA, ACA is easier to be applied for compression of impedance matrix, because it does not need to decompose Green’s function. Besides, according to the low-rank property of the matrix, a matrix with a rank can be approximately decomposed as a product of matrices and . The formulation is listed as follows: where is the size of matrix and and are the sizes of the matrices and , respectively.
For off-diagonal matrices and shown in 3, they can be compressed by ACA method directly when the distance between the surface and target is large enough. When the distance is not large enough, some parts of the off-diagonal matrix are near field and other parts are far field. Obviously, only far field parts can be compressed by low-rank properties. Near field impedance elements are computed by MoM directly. To find the far field part of the off-diagonal matrix, more than one octree should be applied. One octree structure is used for the rough surface to divide the whole mesh into many groups. The same as the surface, each target requires a private octree structure, either. To model the action between a mesh group from the surface and a mesh group from a target, the distance between the two groups should be computed first. According to our experience, when the distance is larger than 0.4, where is the wavelength, ACA will be applied to compress the coupling action.
The algorithm of ACA allows to generate only few rows and columns of the matrix and approximates the rest of the matrix using only this information. According to the principle of ACA algorithm, only the elements of specific rows and columns should be computed. A large amount of time consumption for computing the entire matrix is saved, and the storage consumption is reduced from to . The smaller the rank of the matrix is, the better the matrix compression effect. The computational complexity reaches at low frequency and at medium and high frequencies [11].
2.3. Implementation Strategy
The implementation plan of our investigated method in this paper is depicted in Figure 2. There are two octrees, one is used for surface modeling and another one is used for a target, where the minimum length of each group is less than in the whole paper.
It should pay attention that there are many nonempty groups in the finest level of the octrees from each target and each nonempty group contains many edges and triangles. According to the tree structure of the target, the interaction between two nonempty groups in the finest level of the tree structure is compressed by ACA algorithm. In 2, the impedance matrices generated by the coupling interaction between the two objectives are the nondiagonal block matrices and , where is an interacting matrix from target 2 to target 1 and is an interacting matrix from target 1 to target 2, where the expression of is listed as follows:
Therein, element is the result of interaction from the first basis functions in target 2 to first basis functions in target 1 and is the result of interaction from the jth basis function in target 2 to the ith basis function in target 1. Suppose that the surfaces of the two targets are divided into and basis functions, respectively. Therefore, the matrix size is . The interaction between two nonempty groups in two targets will produce many submatrices , where equals the size of the nonempty group in target 2 multiplied by the size of the nonempty group in target 1. The impedance matrix generated by the interaction of the first nonempty group in target 2 to the first nonempty group in target 1 is denoted as , which can be decomposed into two matrices and by using the ACA decomposition method as shown in 5.
In the same way, we can decompose the submatrices generated by every two nonempty groups by ACA and the expression is expressed as follows:
Thus, the submatrix and Ne submatrix V can be obtained. Besides, since the number of basis functions included in each group is unequal, the size of the subimpedance matrix generated by the interaction of each two groups also would be greatly different. For example, if the basis function in the fifth group in target 1 is 200 and the number of basis functions in the 37th group in target 2 is 2, then the size of the matrix generated by the lease of the two nonempty groups is 200 × 2. The rank of the matrix is less than or equal to 2. According to the definition of ACA, the generated matrix is 200 × and the generated matrix size is × 2. Considering that ACA algorithm has better compression effects on lower-rank matrices, however, when the compressed matrix rows or columns are too small, the rank does not reach much less than the requirement of matrix dimensions. Generally speaking, when the row or column of the compressed matrix is less than 5, the traditional method of moments can be used to solve the calculation.
In the filling calculation of matrix, the ACA algorithm only needs to use the information of the partial basis functions in the nonempty group. That is, using the information of the partial rows and columns in the decomposed matrix, the original matrix can be approximated. According to the symmetry of the target mutual coupling, matrices and are symmetric and is the interaction of target 1 to target 2. Similarly, the compression of using the ACA algorithms also has the same principle and the same computational compression efficiency.
3. Numerical Results
In the simulation process, the test computer used has 64-bit operating systems, Windows 10 version. The parameter of the CPU is Intel(R) Core(TM) i5-6500U CPU @ 3.20 GHz 3.19 GHz. The double data used to support the findings of this study are available from the corresponding author upon request.
To verify the feasibility of the algorithm, the bistatic radar cross section (RCS) [22] of the metal pipe is simulated under the random rough surface. The resulting impedance matrices were iteratively solved using the generalized minimum residual (GMRES) solver, and the dimension of the Krylov subspace is set to 30. The tolerance of relative error is 10−3 [23].
3.1. Scattering of Two Spheres
First, place two spheres of PEC media in free space with the center of the sphere on the x-axis and symmetrically about the origin. The incident wave is a plane wave, the incident frequency is 300 MHz, and the incident angles are ° and °. The radiuses of two balls are set to 2 m and the center distance is 5 m. We can conclude that the nearest distance from the target surface is 1 m and the farthest distance is 9 m. In addition, we can obviously see that the surface distance of the target changes in a wide range and has small interactions when they are far away; thus, using ACA compression method of moments can achieve a good result. Conversely, they will have strong interactions when they are close with each other and this strong effect will lead to less low-rank matrix. In this circumstance, the grouping calculation has the ability to achieve different compression efficiencies with different rank matrices, even though the overall matrix will exhibit the characteristic of low rank. The number of levels for each MLFMA octree is 3. The box size is 0.2–0.4 with 153 edges per square wavelength. The tolerance of ACA is 10−3. Figure 3(a) illustrates the HH-polarized bistatic RCS.
(a)
(b)
We still calculate two spheres whose radii are 2 m; the center distance is 4.5 m. Other settings are the same as above. The HH-polarized bistatic RCS is shown in Figure 3(b). At this point, the nearest surface distance of the target is 0.5 m, equal to 1/2 of the wavelength, at which point the distance is very close.
Computational parameters of bistatic RCS of two spheres which have the same radius (2 m) and different centers of the distance are shown in Table 1. By using the traditional method of moments combined with MLFMA, our proposed new methods are MLFMA and ACA for MoM. Through the comparative numerical analysis demonstrated in Table 1, it can be obviously seen that the time of filling matrix and iteration time are significantly reduced and the calculation of the RCS curves almost coincides, which illustrated that the algorithm not only can improve the computational efficiency but also does not affect the calculation accuracy.
Then with the increase of the radius, the center distance of the sphere keeps three times the radius and the division scale remains unchanged. The two spheres with radii 0.5 m, 1 m, 2 m, 4 m, and 8 m were calculated. The numbers of levels for MLFMA are 1, 2, 3, 4, and 5, respectively. The MLFMA box size remained to be 0.2–0.4, and the tolerance of ACA remained to be 10−3. In this simulation process, the test computer used has 64-bit operating systems, Windows 10 version. The parameter of the CPU is Intel(R) Core(TM) i7-8700U CPU @ 3.20 GHz 3.19 GHz, and the RAM is 16.0 GB. We calculate the calculation time of the two spheres scattering, and the calculated memory consumption is shown in Figure 4. In the graph, the x-axis represents which is the number of unknowns.
(a)
(b)
In Figure 4, the comparison diagrams of matrix filling time, iteration time, and total time are listed. By analyzing the comparison diagram of the total time in the graph, we can see that the computational complexity is almost , while is the number of unknowns.
Finally, we calculate the two balls with unknown magnitude. Two spheres with an unknown value of 98304 are added. The radius of the sphere is set to 20 m and the center distance is 60 m. The number of levels for each MLFMA octree is 6, and the tolerance of ACA is 10−3. The incident frequency is 300 MHz. In this simulation process, the test computer used has 64-bit operating systems, Windows 10 version. The parameter of the CPU is Intel(R) Core(TM) i7-4770 U CPU @ 3.40 GHz 3.40 GHz. The HH-polarized bistatic RCS is shown in Figure 5. The filling time is 6798 s, the iteration time is 3662 s, and the total computation time is 10470 s.
3.2. Scattering of Underground Ball
In order to verify our finding that the MoM-MLFMA-ACA method has a better compression performance when calculating two PEC targets, the scattering of underground targets is tested as an example as follows.
Suppose that there is a ground, and the ground is not a smooth surface. Therefore, in the simulation process, it is necessary to simulate the surface of the rough program. Random rough surfaces are generated by spectral FFT method and Gauss spectrum [24]. Then the rough surface is set to 5 × 5 m, where a ball for the ideal conductor sphere center distance rough surface is 2 m underneath it with the radius of the ball being 0.5 m. Figure 6 illustrates a sketch map of the rough surface and metal ball, and also suppose that the incident wave is a plane wave and the incident angles are ° and °.
In order to compare the effects of different soils on scattering calculations, we give different soil parameters. For simplicity, the relative permeability of the soil is and sets of different relative dielectric constants are and . In addition, the RCS of different incident wave frequencies ( MHz, MHz) is calculated. When MHz and , the density of the upon rough surface is 949 edges per square wavelength and the number of levels for MLFMA is 3, while the below rough surface is 474 edges per square wavelength and the number of levels is 4. More detailed MLFMA data is shown in Table 2. The HH-polarized bistatic RCS is shown in Figure 7, and the calculation parameters of RCS are shown in Table 3, where we can see obviously that the unknown quantity is calculated to be 15280 and the iteration time is greatly compressed and the computation memory consumption is also decreased significantly. As the dielectric constant increases, the computation iteration time increases, while the time difference between the two methods remains large. The increase of the dielectric constant affects the matrix condition of the matrix, which increases the number of iteration steps of the matrix, as shown in Table 3, resulting in an increase in the iteration time. The memory compression performance decreases with the increase of dielectric constant. As the dielectric constant increases, the rank of most compression blocks in the impedance matrix increases, so the storage capacity of most matrix blocks increases and the overall storage capacity increases.
(a)
(b)
(c)
(d)
3.3. Scattering of Underground Pipelines
Consider a metal pipeline buried below the rough surface, with a radius of 0.5 m, a length of 4 m, and a coordinate of the center of the pipeline (, ). The rough surface is set to 6 × 6 m, the relative permittivity of soil epsilon is , relative permeability is , and ground is placed at position. The incident wave is a plane wave, the incident frequency is 100 MHz, and the incident angles are and . The rough surface density of the upper surface is 205 edges per square wavelength, and the number of levels for MLFMA is 3, while the lower surface is 103 edges per square wavelength and the number of levels is 4. The pipeline density is 202 edges per square wavelength, and the number of levels for MLFMA is 3. The bistatic RCS is calculated by simulation, as shown in Figure 8.
The depth of target burial is often uncertain, and different depths can achieve different compression efficiencies. For realizing the above calculation, different depths are set, while other conditions are unchanged. The calculated storage changes and time variations are shown in Table 4. As the distance increases, both the memory consumption and the fill matrix time are reduced in the new method. The new method is a combination of the ACA algorithm and the MoM-MLFMA method, in which ACA replaces the pure MoM used by the interaction block in the old method. Therefore, the storage of the nondiagonal block impedance matrix of the new method is greatly compressed. Moreover, the farther the target distance is, the more obvious the compression performance is.
4. Conclusions
In this paper, a hybrid method combined with both MLFMA and ACA is applied, which can simulate the electromagnetic scattering of the underground target efficiently. To verify the effectiveness of our investigated method, three kinds of actions are defined. The first kind is the self-action of the surface. Then the second kind is the self-action of each target. And the last is the coupling-action between the surface and targets. Self-action is compressed and accelerated by MLFMA while coupling-action is compressed and accelerated by ACA. According to the demonstrated numerical results of scattering from underground targets, our proposed method can reduce the impedance filling time and solution time, and at the meantime, the memory can also be saved greatly. In a word, the proposed MLFMA combined with ACA algorithm is an accurate and effective hybrid method for calculating the scattering of subsurface targets.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (61601185, 81460275) and the Natural Science Foundation of Jiangxi Province (20171ACB21040).