Next Article in Journal
Multi-Scale Remote Sensing Assessment of Ecological Environment Quality and Its Driving Factors in Watersheds: A Case Study of Huashan Creek Watershed in China
Next Article in Special Issue
A Deep Learning Gravity Inversion Method Based on a Self-Constrained Network and Its Application
Previous Article in Journal
Performance Analysis of Artificial Intelligence Approaches for LEMP Classification
Previous Article in Special Issue
Spreading of Localized Information across an Entire 3D Electrical Resistivity Volume via Constrained EMI Inversion Based on a Realistic Prior Distribution
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Efficient 3D Frequency Semi-Airborne Electromagnetic Modeling Based on Domain Decomposition

1
College of Geophysics, Chengdu University of Technology, Chengdu 610059, China
2
College of Geo-Exploration Sciences and Technology, Jilin University, Changchun 130026, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(24), 5636; https://doi.org/10.3390/rs15245636
Submission received: 23 October 2023 / Revised: 29 November 2023 / Accepted: 1 December 2023 / Published: 5 December 2023
(This article belongs to the Special Issue Multi-Data Applied to Near-Surface Geophysics)

Abstract

:
Landslides are common geological hazards that often result in significant casualties and economic losses. Due to their occurrence in complex terrain areas, conventional geophysical techniques face challenges in early detection and warning of landslides. Semi-airborne electromagnetic (SAEM) technology, utilizing unmanned aerial platforms for rapid unmanned remote sensing, can overcome the influence of complex terrain and serve as an effective approach for landslide detection and monitoring. In response to the low computational efficiency of conventional semi-airborne EM 3D forward modeling, this study proposes an efficient forward modeling method. To handle arbitrarily complex topography and geologic structures, the unstructured tetrahedron mesh is adopted to discretize the earth. Based on the vector finite element formula, the Dual–Primal Finite Element Tearing and Interconnecting (FETI-DP) method is further applied to enhance modeling efficiency. It obtains a reduced order subsystem and avoids directly solving huge overall linear equations by converting the entirety problem into the interface problem. We check our algorithm by comparing it with 1D semi-analytical solutions and the conventional finite element method. The numerical experiments reveal that the FETI-DP method utilities less memory and exhibits higher computation efficiency than the conventional methods. Additionally, we compare the electromagnetic responses with the topography model and flat earth model. The comparison results indicate that the effect of topography cannot be ignored. Further, we analyze the characteristic of electromagnetic responses when the thickness of the aquifer changes in a landslide area. We demonstrate the effectiveness of the 3D SAEM method for landslide detection and monitoring.

1. Introduction

Landslides are recognized as a global natural hazard that can result in the loss of human lives, damage to property, and destruction of infrastructures [1]. Landslide monitoring is one of the most significant measures to reduce landslide hazards [2]. Remote sensing and satellite interferometric synthetic aperture radar (InSAR) are the main monitoring methods, but they have difficulty identify potential landslides due to their inability to obtain underground characteristics [3]. The controlled-source electromagnetic (CSEM) method is an effective means of acquiring underground electrical information and has been successfully applied in engineering and environmental investigations [4,5,6,7]. By utilizing underground electrical information, the underground structure of the slope and groundwater characteristics can be obtained. The groundwater constructs the interaction with the geo-environment, and groundwater variation can destroy mechanical stability, which is a critical indicator for evaluating landslide stability [8,9].
The conventional ground CSEM methods suffer from low efficiency, especially in complex landslide geologic areas, due to the need to lay out transmitting and receiving equipment on the ground. In contrast, the SAEM method addresses this issue by installing a receiver on an unmanned aerial vehicle (UAV), which enables the efficient recording of the EM field. Additionally, a grounded source such as the transmitter is free from the restriction of the aircraft carrying capacity encountered by the airborne electromagnetic (AEM) method and ensures large transmission power [10]. Over the past few decades, the SAEM method has developed rapidly and has been widely used in engineering, hydrology, and disaster detection [11,12,13]. In actual SAEM exploration areas, the topography and subsurface structures can be arbitrarily complicated. To obtain the accurate solution, an effective 3D forward modeling method suitable for complex structures is crucial.
Over the past two decades, 3D EM forward modeling has made great progress. Finite difference (FD) [14,15,16], finite volume (FV) [17,18,19], and finite element (FE) [20,21,22] are the main categories of methods used for modeling 3D EM fields. For the SAEM method, Smirnova et al. [23] used the FD approach to study Kiruna Iron Ore fields. Wang et al. [24] utilized staggered grid finite difference method to develop SAEM software. He et al. [25] and Jing et al. [26] applied the staggered mimetic FV method to SAEM exploration. Based on the finite element method using a structural grid, Xue et al. [27] analyzed the detection ability of the SAEM method for underground cavities. Among these methods, structured hexahedron grids were utilized to divide the models, so it cannot accurately fit the topography and complex geologic structure [28]. The FV method can effectively resolve complex geological structures using orthogonal Delaunay–Voronoï unstructured grids, as its implementation is complicated and requires high mesh quality [29]. To be suitable for modeling complex terrain and geologic structures, the FE method, using unstructured grids, has been widely used in SAEM modeling [30,31].
The FE method with vector basic function was developed by Nédélec [32]. Due to the face that it is divergence-free in the element, the vector FE method has been widely used in CSEM and magnetotelluric (MT) with first-order edge elements [22,33,34]. Research has shown that there must be fine mesh to obtain high-precision solutions when using low-order elements [35]. To reduce the influence of mesh and ensure solution precision, the high-order elements for the FE method are introduced into CSEM and MT modeling [36,37,38,39]. However, for arbitrarily complex topography and geologic structure, highly dense grids are needed to accurately fit the complicated geometric shape [40,41]. The number of degrees of freedom (dof) in the calculation domain can reach millions or more. So, the FE method usually faces the problem of solving a huge sparse linear equation. The direct method is a mainstream approach in 3D EM modeling to solve the linear equations [42,43,44], while its computational complexity is the order O(N2) and consumes a huge amount of memory due to non-zero elements filling during factorization [45]. To reduce the scale of dof, some researchers have proposed h-adaptive refinement methods based on posterior error estimation [46,47]. Castillo-Reyes et al. developed the hp-refinement method by using the advantages of high-order elements, and they utilized the parallel strategy to improve computational efficiency [48]. On the other hand, Farhat et al. developed the FETI method, which combines the FE method and domain decomposition method (DDM) [49]. It divides the whole calculation domain into many nonoverlapping subdomains and introduces the Lagrange multipliers to couple all subdomains, which converts the entirety problem into the interface problem and avoids directly solving huge overall linear equations. Farhat et al. [50] forcedly imposed the continuity condition at the cross-points to improve algorithm stability, called the FETI-DP method. Li and Jin [51] developed the FETI-DP method with the vector FE method to simulate open-region electromagnetic propagation problems. Zhang and Jin [52] analyzed the parallel efficiency for the FETI-DP method. When the number of cores is more 400, its parallel efficiency is greater than 0.8. In geo-electromagnetic studies, Ren et al. [53] used the FETI-DP method to solve MT problems. Hui et al. [54] proposed a time-domain approach based on the FETI-DP method to 3D marine CSEM modeling. The numerical results showed that the FETI-DP method in the time domain has the same calculation accuracy as the conventional finite element method. Inspired by the successful application of the FETI-DP method, this paper developed SAEM forward modeling research based on the FETI-DP method, aiming to improve its computational efficiency.
In this paper, we adopted the FETI-DP method, which combines the FE and domain decomposition techniques, to solve frequency domain SAEM problems. We first introduced the technical details of the FETI-DP method. Then, to demonstrate the accuracy of our approach, we checked the results against 1D semi-analytical solutions. Furthermore, we compared the computational accuracy and efficiency of the algorithm used in this paper with the conventional FE method. Finally, considering the complex topography of landslide areas, we designed realistic models to analyze the influence of terrain on the SAEM responses and the identification ability of underground aquifers in landslide areas for the 3D frequency SAEM method.

2. Methods

2.1. Governing Equations

Taking quasi-static approximation, the frequency-domain Maxwell’s equations can be written as:
× E = i ω μ H ,
× H = σ E + J s ,
where E is the electric field, H is the magnetic field, ω is the angular frequency, μ is the magnetic permeability, σ is the conductivity, and J s is the imposed source current. Combining Equations (1) and (2), and eliminating H, we can obtain the curl–curl equation for electric field:
× × E + i ω μ σ E = i ω μ J s ,
To ensure a unique solution, we imposed the Dirichlet boundary condition on the outer surface of the whole modeling domain:
n × E = 0 on Ω ,
where n denotes the outer normal unit vector of the boundary Ω . Ω denotes the whole modeling domain.

2.2. Domain Decomposition Method

We divided the whole modeling domain into Ns nonoverlapping subdomains using METIS [55]. For each subdomain, there is approximately the same number of unknowns, so that they consume roughly the same computing resources to obtain high parallel efficiency for the FETI-DP method. Except for the outer boundary Ω , there is also an inner boundary among subdomains for the pth subdomain (shown in Figure 1). We utilize the unknown Newman boundary condition for the inner boundary:
n × 1 μ × ( × E p ) = Λ p on Γ p ,
where Γ p denotes the interface of pth subdomain, Λ p denotes an unknown value. We use the first order of vector basis functions [56] and apply Galerkin’s method and the boundary condition Equations (4) and (5), obtaining the following FE equation for pth subdomain
( A p + i ω B p ) E p = i ω S p λ p ,
where A p and B p are, respectively, stiffness and mass matrices for the pth subdomain, S p is the source term, and λ p is interface boundary term, called the Lagrange multiplier [50]. The explicit expressions are as follows:
A p = 1 μ Ω p   ( × N p ) ( × N p ) d V ,
B p = Ω p   N p σ N p d V ,
S p = Ω p   N p J s d V ,
λ p = Γ p N p Λ p d S .
For the convenience of representations, we rewrite Equation (6) as:
K p E p = f p λ p .
where K p = A p + i ω B p , f p = i ω S p . Due to the usage of unknown Newman boundary conditions, the Lagrange multiplier is unknown. So, in order to solve this equation, we need to first solve the Lagrange multiplier.
In each subdomain, the edges are divided into three categories: inner edges (in subdomain) denoted by I, interface edges (shared by two subdomains) denoted by f, and corner edges (shared by more than two subdomains) denoted by c. Figure 1 shows the diagram of edge classification.
The inner edges and interface edges are called the remaining edges, denoted by r. Then, Equation (11) can be written as:
[ K r r p K r c p ( K r c p ) T K c c p ] [ E r p E c p ] = [ f r p f c p ] [ λ r p λ c p ] .
To solve the Lagrange multiplier, we need to couple all subdomains, which converts the overall mesh problem into solving the interface Lagrange multiplier (shown in Figure 2).
We build the mapping between local variable and global variable by defining a signed Boolean matrix B r p and an unsigned Boolean matrix B c p , as follows:
λ r p = B r p λ f ,
E c p = B c p E c .
where λ f denotes the global interface Lagrange multiplier on the interface. Using B r p and B c p to rewrite Equation (12), we have:
K r r p E r p + K r c p B c p E c = f r p B r p λ f ,
( K r c p ) T E r p + K c c p B c p E c = f c p λ c p .
Combining Equations (15) and (16), and eliminating E r p , we can obtain:
( K c c p ( K r c p ) T ( K r r p ) - 1 K r c p ) K c p E c = f c p λ c p ( K r c p ) T ( K r r p ) - 1 f r p + ( K r c p ) T ( K r r p ) - 1 B r p λ f .
We take into account that p = 1 N s ( B c p ) T λ c p = 0 [51]. Multiplying both sides of Equation (17) by ( B c p ) T and assembling all subdomains, we can obtain:
K ˜ c c E c = p = 1 N s [ ( B c p ) T ( f c p ( K r c p ) T ( K r r p ) - 1 f r p + ( K r c p ) T ( K r r p ) - 1 B r p λ f ) ] ,
where
K ˜ c c = p = 1 N s [ ( B c p ) T ( K c c p ( K r c p ) T ( K r r p ) - 1 K r c p ) B c p ] .
According to the tangential continuity condition of the electric field, multiplying both sides of Equation (15) by ( B r p ) T and assembling all subdomains, we can obtain:
p = 1 N s [ ( B r p ) T ( D r r p ) 1 ( f r p B r p λ f D r c p B c p E c ) ] = 0 .
Eliminating E c from Equations (18) and (20), we can obtain:
F λ f = b ,
where
F = F r c K ˜ c c 1 F r c T + F r r ,
F r r = p = 1 N s [ ( B r p ) T ( K r r p ) 1 B r p ] ,
F r c = p = 1 N s [ ( B r p ) T ( K r r p ) 1 K r c p B c p ] ,
b = f ˜ r F r c K ˜ c c 1 f ˜ c ,
f ˜ r = p = 1 N s [ ( B r p ) T ( K r r p ) 1 f r p ] ,
f ˜ c = p = 1 N s [ ( B c p ) T ( f c p ( K r c p ) T ( K r r p ) - 1 f r p ) ] .
By solving Equation (21), we can obtain λ f . According to Equations (18) and (20), the electric field on all edges can be calculated. Figure 3 shows the flow diagram of the FETI-DP method.

2.3. Solutions

According to Equations (21)–(27), it is necessary to obtain K ˜ c c 1 and ( K r r p ) 1 when calculating the λ f , but calculating the inverse of the matrix is very difficult. In order to avoid this problem, we use adjoint forward modeling to calculate the product of the inverse of the matrix and a vector. K ˜ c c 1 and ( K r r p ) 1 are, respectively, the N c × N c and N r p × N r p matrices. N c is the number of global core edges, and N r p is the number of subdomain remaining edges. Due to its small size, the direct solver MUMPS [57] can be used.
In Equation (21), the coefficient matrix F is a dense Nf × Nf matrix, where Nf is the number of edges on global interfaces. If using a direct solver, it is necessary to construct the coefficient matrix F , which will consume a large amount of memory. And, in this process, we have to solve a large number of linear equations with coefficient matrix K r r p when calculating F r r and F r c , which is highly computationally complex. So, we adopt the generalized minimal residual (GMRES) [58] to iteratively solve the linear equations system. It can avoid constructing the coefficient matrix explicitly and we only need to calculate matrix F a few times and multiply a vector at each iteration. To accelerate the iterative convergence, we use P 1 as a preconditioner [51], which can be expressed as:
P 1 = p = 1 N s ( B r p ) T [ 0 0 0 Q f f p ] B r p ,
where:
Q f f p = K f f p ( K I f p ) T K I I p K I f p ,

3. Numerical Experiments

3.1. Accuracy Verification

To verify the accuracy of the algorithm for SAEM modeling, we designed a three layer model (as shown in Figure 4). In this model, the thickness of the first two layers is 100 m and 300 m. The resistivity of the air is 106 Ω m , while the resistivity of three underground layers is 100, 10, and 100 Ω m , respectively. We take the 1000 m long grounded electrical transmitter in our modeling, aligned in the x direction. The transmitting current is 1A. The receivers are aligned in the y direction and their x coordinates are the same as the coordinate of the transmitting source midpoint. The flight altitude of the receivers is 30 m above the surface. For the modeling, we divided the computational domain into 739,935 tetrahedral elements and then partitioned the mesh into 80 subdomains using METIS. The mesh and its partitions are shown in Figure 5. We calculated EM responses at 21 frequencies from 1 Hz to 10,000 Hz.
Figure 6 shows the comparison of the results between the FETI-DP and the 1D semi-analytical solutions for Bz response. It can be seen that our results match the 1D semi-analytical solutions well. The relative errors for amplitude are less than 3%. The maximum phase error is less than 0.8 degrees. The iterative convergence information for partial frequencies is given in Table 1.

3.2. Comparison with Finite Element Method

In the following, we compare the algorithm with the conventional FE method using the direct solver MUMPS. We designed a landslide model as show in Figure 7a,b. In this model, the resistivity of the air is set to 10−6 Ω m and the background’s resistivity is 100 Ω m . There is a low resistivity aquifer below the surface, whose resistivity is 10 Ω m . We take the 900 m long grounded electrical wire as the transmitting source in our modeling, aligned in the x direction. The receivers are aligned in the y direction and their x coordinates are x = 0 m. The flight altitude of the receivers is 30 m above the surface. We also calculated EM responses at 21 frequencies from 1 Hz to 10,000 Hz.
We used the same mesh for the FETI-DP and FE methods when calculating. Figure 8a,b shows the mesh and partitions. We divided the computational domain into 756,629 tetrahedral elements and 864,130 edges. Then, we partitioned the mesh into 80 subdomains. Figure 9 shows the comparison of the results between the FETI-DP method and FE method. The relative errors for amplitude are less than 0.25%. The maximum phase error is less than 0.08 degrees. The numerical results showed that the FETI-DP method, when used in frequency domain SAEM modeling, has the same calculation accuracy as the conventional FE method.
To compare the computational resource requirements of the algorithm proposed in this article and the conventional FE method, we refined the grid in Figure 8 using a predetermined threshold. We divided the computational domain into 2,025,999 tetrahedral elements and 2,357,140 edges. Then, we partitioned the mesh into 200 subdomains. We use a workstation with Intel(R) Xeon(R) CPU E5-2667 @ 3.20 GHz. The single-threaded calculation result for frequency at 1000 Hz is given in Table 2. From Table 2, it can be seen that the FETI-DP method has a faster computational speed than the conventional FE method when the number of grids increases. According to the characteristics of the FETI-DP method, the computational efficiency can be further improved when large-scale parallel computing is carried out.

3.3. Numerical Experiments on Landslide Model

In landslide areas, the topographic effect can be a serious problem. We designed a half space model with a topography as shown in Figure 10a. In this model, the resistivity of the air is set to 10−6 Ω m and the background’s resistivity is 100 Ω m .We take the 900 m long grounded electrical wire as the transmitting source in our modeling, aligned in the x direction and midpoint at (x = 150 m, y = 450 m). There are ten survey lines at a line spacing of 100 m from y = −600 m to 300 m. The receivers, at each survey line, are located from −500 m to 500 m with an interval of 50 m in the x direction. The flight altitude of the receivers is 30 m above the surface. The receivers cover topography surface as shown in Figure 10b. Taking the tenth survey line as an example, we study the influence of topography on SAEM responses. Figure 11 shows the responses of the landslide model and flat terrain model for frequencies of 10, 100, and 1000 Hz. From Figure 11, it can be seen that landslide terrain has a significant impact on electromagnetic response, especially at 1000 Hz, with a maximum relative anomaly of 70%, indicating that the impact of terrain cannot be ignored for SAEM in the frequency domain. And, it is difficult to summarize the rules about terrain effect for different frequencies, so it is necessary to simulate it with topography.
In order to analyze the identification ability of SAEM for landslide hazards, we add a three-layer aquifer into the half space model with topography, as shown in Figure 12, to simulate aquifer variation. The resistivity of the aquifers is 10 Ω m . The thickness of both layer 2 and layer 3 is 10 m. The resistivity of the air is set to 10−6 Ω m and the background resistivity is 100 Ω m . Figure 13 shows the information about sources and receivers. We take two 900 m long grounded electrical wires as the transmitting sources and compare their electromagnetic responses. The two sources are both aligned in the x direction. The midpoints of the two sources are, respectively, at (x = 150 m, y = 450 m) and (x = 150 m, y = 0 m). The transmission frequency is 1000 Hz. Figure 14 shows the electromagnetic responses for source 1 at y = −500 m survey line. Figure 14a show the amplitude of Bz for the model in Figure 12a, while Figure 14b is the relative percentage anomaly between the different layer model and half space model. From Figure 14, we can see that there is a significant difference in the amplitude of Bz between the aquifer model and half space model. As the thickness of the aquifer increases, the relative anomaly significantly rises, with a maximum increase of 10%. This result indicates that the SAEM method is effective for monitoring aquifer variations. The aquifer variation can destroy mechanical stability. So, the SAEM method has the potential to predict landslide hazards. Figure 15 shows the computed real part and imaginary part of the electric fields at y = −500 m. The electric field for the SAEM method is weaker in the conductive aquifer area, which demonstrates the fundamental basis for identifying a conductive body. Figure 16 shows the electromagnetic responses for source 2 at y = −500 m survey line. Figure 16a show the amplitude of Bz for the model in Figure 12a, while Figure 16b is the relative percentage anomaly between the different layer model and half space model. Comparing Figure 14 and Figure 16, we can see that source 2 has a larger relative percentage anomaly because it is closer to the anomaly body. When detecting the landslide hazards using the SAEM method, we can obtain a better result by setting the source on the slope.

4. Discussion

The advantage of the SAEM method is that the receivers are located on UAV, which enables the efficient recording of the EM field. Meanwhile, a grounded source can ensure that the large transmission power has a strong signal. To fit the arbitrarily complex topography and geologic structure, the FE method has to solve an equations system with a huge dof. The FETI-DP method combines the FE method and DD method. It divides the whole calculation domain into many nonoverlapping subdomains and introduces the Lagrange multipliers to couple all subdomains, which converts the entirety problem into the interface problem and avoids directly solving huge overall linear equations. When using the same mesh, the FETI-DP method and the FE method have almost consistent calculation accuracy. However, the FETI-DP method has high computational efficiency and low memory consumption, which means that the FETI-DP methos is a useful method for SAEM modeling.
In the landslide area, we can see from Figure 11 that the influence of topography for SAEM cannot be ignored, and it is difficult to summarize the rule to remove terrain effects. So, it is necessary to simulate with topography when calculating SAEM responses in landslide areas. The aquifer variation can destroy mechanical stability, which is a critical indicator for predicting landslide hazards. Utilizing a three-layer abnormal body model, we simulate the aquifer variation. From Figure 14, it can be seen that the Bz responses have a significant difference when the thickness of the aquifer changes by 10 m. This result indicates that the SAEM method is effective for monitoring aquifer variations and it has the potential to predict landslide hazards. By calculating the response of the source at different locations, we can see that it has a stronger single and relative percentage anomaly when the source is close to the aquifer. So, a better result can be obtained by setting the source on the slope when detecting landslide hazards using the SAEM method.
The induced polarization (IP) effect can be observed when the SAEM method is used with a grounded wire source. In the future, we will consider the IP effect when modeling.

5. Conclusions

In this paper, based on the unstructured finite element method and domain decomposition method, we have developed an effective 3D modeling scheme FETI-DP for the SAEM method in the frequency domain. We have verified our algorithm by comparing our results with 1D semi-analytical solutions. Further comparing our result with the conventional FE method, we have verified that the FETI-DP method has the same accuracy as the FE method when using the same mesh. However, the FETI-DP method has high computational efficiency and low memory consumption. Thus, the FETI-DP method can be a useful tool for SAEM modeling. The numerical results for a landslide model showed that EM responses are seriously affected by the topography, and it is difficult to summarize the rules, so it is necessary to simulate with topography. The final numerical experiment shows that the SAEM method can identify aquifer variation with the simulation of topography, which proves that it is a potential method for landslide monitoring. And, setting the source on the slope can lead to obtaining a better result.

Author Contributions

Conceptualization, Z.H., X.W. and C.Y.; methodology, Z.H. and C.Y.; software, Z.H.; formal analysis, Z.H., C.Y. and Y.L.; investigation, Z.H., X.W., C.Y. and Y.L.; writing—original draft preparation, Z.H.; writing—review and editing, Z.H., X.W., C.Y. and Y.L.; funding acquisition, X.W. and C.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This paper is financially supported by National Key R&D Program of China: (Project No. 2022YFC3003203) and National Natural Science Foundation of China: (Project No. 42304094).

Data Availability Statement

Data associated with this research are available and can be obtained by contacting the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lutgens, K.L.; Tarbuck, E.J. Foundations of Erath Sciences. Pearson International Edition; Pearson: London, UK, 2008. [Google Scholar]
  2. Xu, Q.; Zhu, X.; Li, W.; Dong, X.; Dai, K.; Jiang, Y.; Lu, H.; Guo, C. Technical progress of space-air-ground collaborative monitoring of landslide. Acta Geod. Cartogr. Sin. 2022, 51, 1416–1436. [Google Scholar]
  3. Xu, Q.; Lu, H.; Li, W.; Dong, X.; Guo, C. Types of potential landslide and corresponding identification technologies. Geomat. Inf. Sci. Wuhan Univ. 2022, 47, 377–387. [Google Scholar]
  4. Kafri, U.; Goldman, M. The use of the time domain electromagnetic method to delineate saline groundwater in granular and carbonate aquifers and to evaluate their porosity. J. Appl. Geophys. 2005, 57, 167–178. [Google Scholar] [CrossRef]
  5. Streich, R. Controlled-source electromagnetic approaches for hydrocarbon exploration and monitoring on land. Surv. Geophys. 2016, 37, 47–80. [Google Scholar] [CrossRef]
  6. Xue, G.; Zhang, L.; Hou, D.; Liu, H.; Ding, Y.; Wang, C.; Luo, X. Integrated geological and geophysical investigations for the discovery of deeply buried gold-polymetallic deposits in China. Geol. J. 2020, 55, 1771–1780. [Google Scholar] [CrossRef]
  7. Zhou, N.; Wei, X.; Zhang, S. Imaging of a shallow magma conduit system based on a high-power frequency-domain controlled-source electromagnetic survey. Geophysics 2023, 88, B47–B54. [Google Scholar] [CrossRef]
  8. Peng, D.; Xu, Q.; Zhang, X.; Xing, H.; Zhang, S.; Kang, K.; Qi, X.; Ju, Y.; Zhao, K. Hydrological response of loess slopes with reference to widespread landslide events in the Heifangtai terrace, NW China. J. Asian Earth Sci. 2019, 171, 259–276. [Google Scholar] [CrossRef]
  9. Huang, R.; Xu, Z.; Xu, M. Hazardous effects of underground water and extraordinary water flow-induced geohazards. Earth Environ. 2005, 33, 001–009. [Google Scholar]
  10. Lin, J.; Xue, G.; Li, X. Technological innovation of semi-airborne electromagnetic detection method. Chin. J. Geophys. 2021, 64, 2995–3004. (In Chinese) [Google Scholar]
  11. Maria, V.S.; Michael, B.; Christian, N.; Pritam, Y.; Wiebke, M.; Raphael, R.; Annika, S.; Stephan, C.; Maxim, Y.S.; DESMEX Working Group. A novel semiairborne frequency-domain controlled-source electromagnetic system: Three-dimensional inversion of semi-airborne data from the flight experiment over an ancient mining area near Schleiz, Germany. Geophysics 2019, 84, 281–292. [Google Scholar]
  12. Wu, X.; Xue, G.; Fang, G.; Li, X.; Ji, Y. The Development and applications of the semi-airborne electromagnetic system in China. IEEE Access 2019, 7, 104956–104966. [Google Scholar] [CrossRef]
  13. Kotowski, P.O.; Becken, M.; Thiede, A.; Schmidt, V.; Schmalzl, J.; Ueding, S.; Klingen, S. Evaluation of a semi-airborne electromagnetic survey based on a multicoper aircraft system. Geosciences 2022, 12, 26. [Google Scholar] [CrossRef]
  14. Newman, G.; Alumbaugh, D. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences. Geophys. Prospect. 1995, 43, 1021–1042. [Google Scholar] [CrossRef]
  15. Commer, M.; Newman, G. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources. Geophysics 2004, 69, 1192–1202. [Google Scholar] [CrossRef]
  16. Avdeev, D. Three-dimensional electromagnetic modelling and inversion from theory to application. Surv. Geophys. 2005, 26, 767–799. [Google Scholar] [CrossRef]
  17. Haber, E.; Ascher, U.M. Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients. SIAM J. Sci. Comput. 2001, 22, 1943–1961. [Google Scholar] [CrossRef]
  18. Haber, E.; Ruthotto, L. A multiscale finite volume method for Maxwell’s equations at low frequencies. Geophys. J. Int. 2014, 199, 1268–1277. [Google Scholar] [CrossRef]
  19. Caudillo-Mata1, L.A.; Haber, E.; Schwarzbach, C. An oversampling technique for the multiscale finite volume method to simulate electromagnetic responses in the frequency domain. Comput. Geosci. 2017, 21, 963–980. [Google Scholar] [CrossRef]
  20. Farquharson, C.; Miensopust, M. Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction. J. Appl. Geophys. 2011, 75, 699–710. [Google Scholar] [CrossRef]
  21. Cai, H.; Xiong, B.; Han, M.; Zhdanov, M. 3D controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method. Comput. Geosci. 2014, 73, 164–176. [Google Scholar] [CrossRef]
  22. Castillo-Reyes, O.; de la Puente, J.; Cela, J. PETGEM: A parallel code for 3D CSEM forward modeling using edge finite elements. Comput. Geosci. 2018, 119, 123–136. [Google Scholar] [CrossRef]
  23. Smirnova, M.; Juhojuntti, N.; Becken, M.; Smirnov, M.; the DESMEX WG. Exploring Kiruna Iron Ore fields with large-scale, semi-airborne, controlled-source electromagnetics. First Break 2020, 38, 35–40. [Google Scholar] [CrossRef]
  24. Wang, K.; Luo, W.; Cao, H.; Wang, X.; Lan, X.; Duang, C. 3-D inversion of UAV semi-airborne electromagnetic method in frequency domain. Chin. J. Geophys. 2021, 64, 1759–1773. (In Chinese) [Google Scholar]
  25. He, Y.; Xue, G.; Chen, W.; Tian, Z. Three-Dimensional Inversion of Semi-Airborne Transient Electromagnetic Data Based on a Particle Swarm Optimization-Gradient Descent Algorithm. Appl. Sci. 2021, 12, 3042. [Google Scholar] [CrossRef]
  26. Jing, X.; Li, X.; Cao, H.; Zhou, J.; Liu, W. The processing workflow of semi-airborne transient electromagnetic data using geological modeling and three-dimensional numerical simulation—A case history in Yunnan, China. J. Appl. Geophys. 2022, 203, 10467. [Google Scholar] [CrossRef]
  27. Xue, G.; Li, X.; Yu, S.; Chen, W.; Ji, Y. The Application of Ground-airborne TEM Systems for Underground Cavity Detection in China. J. Environ. Eng. Geophys. 2018, 23, 103–113. [Google Scholar] [CrossRef]
  28. Liu, Y.; Yin, C.; Qiu, C.; Hui, Z.; Ren, X.; Weng, A. 3-D inversion of transient EM data with topography using unstructured tetrahedral grids. Geophys. J. Int. 2019, 217, 301–318. [Google Scholar] [CrossRef]
  29. Lu, X.; Farquharson, C. 3D finite-volume time-domain modeling of geophysical electromagnetic data on unstructured grids using potentials. Geophysics 2020, 85, 221–240. [Google Scholar] [CrossRef]
  30. Ke, Z.; Liu, Y.; Su, Y.; Wang, L.; Zhang, B.; Ren, X.; Rong, Z.; Ma, X. Three-Dimensional Inversion of Multi-Component Semi-Airborne Electromagnetic Data in an Undulating Terrain for Mineral Exploration. Minerals 2023, 13, 230. [Google Scholar] [CrossRef]
  31. Nazari, S.; Rochlitz, R.; Günther, T. Optimizing Semi-Airborne Electromagnetic Survey Design for Mineral Exploration. Minerals 2023, 13, 796. [Google Scholar] [CrossRef]
  32. Nédélec, J. Mixed finite elements in R3. Numer. Math. 1980, 35, 315–341. [Google Scholar] [CrossRef]
  33. Zhang, Y.; Key, K. MARE3DEM: A three-dimensional CSEM inversion based on a parallel adaptive finite element method using unstructured meshes. In SEG Technical Program Expanded Abstracts 2016; SEG Library: Houston, TX, USA, 2016; pp. 1009–1013. [Google Scholar]
  34. Jahandari, H.; Farquharson, C. 3-D minimum-structure inversion of magnetotelluric data using the finite-element method and tetrahedral grids. Geophys. J. Int. 2017, 211, 1189–1205. [Google Scholar] [CrossRef]
  35. Han, X.; Ni, J.; Yin, C.; Zhang, B.; Huang, X.; Zhu, J.; Liu, Y.; Ren, X.; Su, Y. 3D airborne EM forward modeling based on finite-element method with goal-oriented adaptive octree mesh. Remote Sens. 2023, 15, 2816. [Google Scholar] [CrossRef]
  36. Schwarzbach, C.; Borner, R.; Spitzer, K. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics: A marine CSEM example. Geophys. J. Int. 2011, 187, 63–74. [Google Scholar] [CrossRef]
  37. Castillo-Reyes, O.; de la Puente, J.; García-Castillo, L.; Cela, J. Parallel 3-D marine controlled-source electromagnetic modelling using high-order tetrahedral Nédélec elements. Geophys. J. Int. 2019, 219, 39–65. [Google Scholar] [CrossRef]
  38. Castillo-Reyes, O.; Modesto, D.; Queralt, P.; Marcuello, A.; Ledo, J.; Amor-Martin, A.; García-Castillo, L. 3D magnetotelluric modeling using high-order tetrahedral Nédélec elements on massively parallel computing platforms. Comput. Geosci. 2022, 160, 105030. [Google Scholar] [CrossRef]
  39. Werthmüller, D.; Rochlitz, R.; Castillo-Reyes, O.; Heagy, L. Towards an open-source landscape for 3-D CSEM modelling. Geophys. J. Int. 2021, 227, 644–659. [Google Scholar] [CrossRef]
  40. Castillo-Reyes, O.; Queralt, P.; Marcuello, A.; Ledo, J. Land CSEM simulations and experimental test using metallic casing in a geothermal exploration context: Vallès basin (NE Spain) case study. IEEE Trans. Geosci. Remote Sens. 2021, 60, 4501813. [Google Scholar] [CrossRef]
  41. Castillo-Reyes, O.; de la Puente, J.; Cela, J. HPC geophysical electromagnetics: A synthetic VTI model with complex bathymetry. Energies 2022, 15, 1272. [Google Scholar] [CrossRef]
  42. Streich, R. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy. Geophysics 2009, 74, 95–105. [Google Scholar] [CrossRef]
  43. Li, J.; Farquharson, C.; Hu, X. 3D vector finite-element electromagnetic forward modeling for large loop sources using a total-field algorithm and unstructured tetrahedral grids. Geophysics 2017, 82, E1–E16. [Google Scholar] [CrossRef]
  44. Zhang, M.; Farquharson, C.; Lin, T. Comparison of nodal and edge basis functions for the forward modelling of three-dimensional frequency-domain wire source electromagnetic data using a potentials formulation. Geophys. Prospect. 2022, 70, 828–843. [Google Scholar] [CrossRef]
  45. Li, R. Robust High Performance Preconditioning Techniques for Solving General Sparse Linear Systems; University of Minnesota: Minneapolis, MN, USA, 2015. [Google Scholar]
  46. Ren, Z.; Kalscheuer, T.; Greenhalgh, S.; Maurer, H. A goal-oriented adaptive finite-element approach for plane wave 3D electromagnetic modeling. Geophys. J. Int. 2013, 187, 63–74. [Google Scholar]
  47. Yin, C.; Zhang, B.; Liu, Y.; Cai, J. A goal-oriented adaptive finite-element method for 3D scattered airborne electromagnetic method modeling. Geophysics 2016, 81, 337–346. [Google Scholar] [CrossRef]
  48. Castillo-Reyes, O.; Amor-Martin, A.; Botella, A.; Anquez, P.; García-Castillo, L. Tailored meshing for parallel 3D electromagnetic modeling using high-order edge elements. J. Comput. Sci. 2022, 63, 101813. [Google Scholar] [CrossRef]
  49. Farhat, C.; Roux, F.X. A method of finite element tearing and interconnecting and its parallel solution algorithm. Int. J. Numer. Methods Eng. 1991, 32, 1205–1227. [Google Scholar] [CrossRef]
  50. Farhat, C.; Lesoinne, M.; Letallec, P.; Daniel, R. FETI-DP: A dual–primal unified FETI method-part I: A faster alternative to the two-level FETI method. Int. J. Numer. Methods Eng. 2001, 50, 1523–1544. [Google Scholar] [CrossRef]
  51. Li, Y.; Jin, J. A vector dual-primal finite element tearing and interconnecting method for solving 3-D large-scale electromagnetic problems. IEEE Trans. Antennas Propag. 2006, 54, 3000–3009. [Google Scholar] [CrossRef]
  52. Zhang, K.; Jin, M. Parallel FETI-DP for efficient EM analysis of general objects and antenna arrays. In Proceedings of the 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), Turin, Italy, 7–11 September 2015; pp. 313–316. [Google Scholar]
  53. Ren, Z.; Kalscheuer, T.; Greenhalgh, S.; Maurer, H. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling. Geophysics 2014, 79, 255–268. [Google Scholar] [CrossRef]
  54. Hui, Z.; Yin, C.; Zhang, W.; Liu, Y.; Xiong, B. 3D MCSEM modeling based on domain-decomposition finite-element method. J. Appl. Geophys. 2022, 207, 104847. [Google Scholar] [CrossRef]
  55. Karypis, G.; Kumar, V. A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular Graphs. SIAM J. Sci. Comput. 1999, 20, 359–392. [Google Scholar] [CrossRef]
  56. Jin, J. The Finite Element Method in Electromagnetics; IEEE Press: Piscataway, NJ, USA, 2014. [Google Scholar]
  57. Amestoy, P.R.; Guermouche, A.; L’Excellent, J.Y.; Pralet, S. Hybrid scheduling for the parallel solution of linear systems. Parallel Comput. 2006, 32, 136–156. [Google Scholar] [CrossRef]
  58. Saad, Y. Iterative Methods for Sparse Linear System; SIAM: Philadelphia, PA, USA, 2003. [Google Scholar]
Figure 1. Edge classification. (a) Calculation domain divided into three subdomains p, j, k. (b) The mesh for calculation domain. The white edge in subdomain j represents the inner edge, denoted by I. The white edge on the interface of subdomain p and subdomain j represents the interface edge, denoted by f. The white edge on the intersection of subdomain p, subdomain j, and subdomain k represents the corner edge, denoted by c.
Figure 1. Edge classification. (a) Calculation domain divided into three subdomains p, j, k. (b) The mesh for calculation domain. The white edge in subdomain j represents the inner edge, denoted by I. The white edge on the interface of subdomain p and subdomain j represents the interface edge, denoted by f. The white edge on the intersection of subdomain p, subdomain j, and subdomain k represents the corner edge, denoted by c.
Remotesensing 15 05636 g001
Figure 2. Subdomain coupling. (a) Overall domain mesh divided into 27 subdomains; (b) the interface mesh of the 27 subdomains. The FETI-DP method only needs to calculate the equation for the interface mesh.
Figure 2. Subdomain coupling. (a) Overall domain mesh divided into 27 subdomains; (b) the interface mesh of the 27 subdomains. The FETI-DP method only needs to calculate the equation for the interface mesh.
Remotesensing 15 05636 g002
Figure 3. The flow diagram of the FETI-DP method.
Figure 3. The flow diagram of the FETI-DP method.
Remotesensing 15 05636 g003
Figure 4. A three layer model. The red line is the transmitter. The black points denote the receivers.
Figure 4. A three layer model. The red line is the transmitter. The black points denote the receivers.
Remotesensing 15 05636 g004
Figure 5. A three layer model. (a) Tetrahedral grids. The different colors denote different model properties; (b) partitions, the whole domain is divided into 80 subdomains. The different colors represent different subdomains.
Figure 5. A three layer model. (a) Tetrahedral grids. The different colors denote different model properties; (b) partitions, the whole domain is divided into 80 subdomains. The different colors represent different subdomains.
Remotesensing 15 05636 g005
Figure 6. Accuracy verification for the three layer model. (a) Comparison between our results and 1D semi-analytical solutions for amplitude; (b) relative errors for amplitude; (c,d) are for phase and relative error.
Figure 6. Accuracy verification for the three layer model. (a) Comparison between our results and 1D semi-analytical solutions for amplitude; (b) relative errors for amplitude; (c,d) are for phase and relative error.
Remotesensing 15 05636 g006
Figure 7. (a) A landslide model with an aquifer, the red line represents the transmitter, and the black points denote receivers; (b) a slice of the landslide model at y = −200 m.
Figure 7. (a) A landslide model with an aquifer, the red line represents the transmitter, and the black points denote receivers; (b) a slice of the landslide model at y = −200 m.
Remotesensing 15 05636 g007
Figure 8. A landslide model with an aquifer. (a) Tetrahedral grids of the model. The different colors denote different model properties; (b) partitions of the model, the whole domain is divided into 64 subdomains. The different colors represent different subdomains.
Figure 8. A landslide model with an aquifer. (a) Tetrahedral grids of the model. The different colors denote different model properties; (b) partitions of the model, the whole domain is divided into 64 subdomains. The different colors represent different subdomains.
Remotesensing 15 05636 g008
Figure 9. Comparison of FETI-DP with finite element method. (a) Comparison between our results and those of finite element for amplitude; (b) relative errors for amplitude; (c,d) are for phase and relative error.
Figure 9. Comparison of FETI-DP with finite element method. (a) Comparison between our results and those of finite element for amplitude; (b) relative errors for amplitude; (c,d) are for phase and relative error.
Remotesensing 15 05636 g009
Figure 10. (a) Topography model, the red line represents the transmitter, and the black points denote receivers; (b) topography surface that the receivers cover.
Figure 10. (a) Topography model, the red line represents the transmitter, and the black points denote receivers; (b) topography surface that the receivers cover.
Remotesensing 15 05636 g010
Figure 11. Comparison of SAEM responses for amplitude between flat terrain and topography. (ac) The responses of flat terrain model for frequencies of 10, 100, and 1000 Hz; (df) the responses of topography model for the same frequencies; (gi) relative percentage difference of responses for flat terrain model and topography model.
Figure 11. Comparison of SAEM responses for amplitude between flat terrain and topography. (ac) The responses of flat terrain model for frequencies of 10, 100, and 1000 Hz; (df) the responses of topography model for the same frequencies; (gi) relative percentage difference of responses for flat terrain model and topography model.
Remotesensing 15 05636 g011
Figure 12. (a) A three-layer aquifer model; (b) a slice of the aquifer model at y = −500 m.
Figure 12. (a) A three-layer aquifer model; (b) a slice of the aquifer model at y = −500 m.
Remotesensing 15 05636 g012
Figure 13. The topography surface, the red lines represent the transmitters, and the black points denote receivers.
Figure 13. The topography surface, the red lines represent the transmitters, and the black points denote receivers.
Remotesensing 15 05636 g013
Figure 14. Comparison of aquifer model and half space model for source 1 in Figure 13; (a) the amplitude of Bz at y = −500 m; (b) relative percentage anomaly.
Figure 14. Comparison of aquifer model and half space model for source 1 in Figure 13; (a) the amplitude of Bz at y = −500 m; (b) relative percentage anomaly.
Remotesensing 15 05636 g014
Figure 15. (a) Real and (b) imaginary parts of the electric fields at y = −500 m.
Figure 15. (a) Real and (b) imaginary parts of the electric fields at y = −500 m.
Remotesensing 15 05636 g015
Figure 16. Comparison of aquifer model and half space model for source 2 in Figure 13; (a) the amplitude of Bz at y = −500 m; (b) relative percentage anomaly.
Figure 16. Comparison of aquifer model and half space model for source 2 in Figure 13; (a) the amplitude of Bz at y = −500 m; (b) relative percentage anomaly.
Remotesensing 15 05636 g016
Table 1. Number of iterations for different frequencies.
Table 1. Number of iterations for different frequencies.
Frequencies (Hz)Number of IterationsRelative Residual Norm
1430.9 × 10−6
10420.7 × 10−6
100420.8 × 10−6
1000410.6 × 10−6
10,000400.8 × 10−6
Table 2. Comparison between the FETI-DP method and the conventional FE method using the same mesh.
Table 2. Comparison between the FETI-DP method and the conventional FE method using the same mesh.
MethodNumber of UnknownMemory (GB)Runtime (s)
FE864,13023.27596
FE2,357,140105.125506
FETI-DP864,13018.08508
FETI-DP2,357,14072.563012
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Hui, Z.; Wang, X.; Yin, C.; Liu, Y. Efficient 3D Frequency Semi-Airborne Electromagnetic Modeling Based on Domain Decomposition. Remote Sens. 2023, 15, 5636. https://doi.org/10.3390/rs15245636

AMA Style

Hui Z, Wang X, Yin C, Liu Y. Efficient 3D Frequency Semi-Airborne Electromagnetic Modeling Based on Domain Decomposition. Remote Sensing. 2023; 15(24):5636. https://doi.org/10.3390/rs15245636

Chicago/Turabian Style

Hui, Zhejian, Xuben Wang, Changchun Yin, and Yunhe Liu. 2023. "Efficient 3D Frequency Semi-Airborne Electromagnetic Modeling Based on Domain Decomposition" Remote Sensing 15, no. 24: 5636. https://doi.org/10.3390/rs15245636

APA Style

Hui, Z., Wang, X., Yin, C., & Liu, Y. (2023). Efficient 3D Frequency Semi-Airborne Electromagnetic Modeling Based on Domain Decomposition. Remote Sensing, 15(24), 5636. https://doi.org/10.3390/rs15245636

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop