PML-coupled three-dimensional magnetotelluric anisotropy parallel finite element forward modeling method

Abstract

The invention relates to a PML-coupled three-dimensional magnetotelluric anisotropy parallel finite element forward modeling method. The method comprises the following steps: according to a continuity condition of an electromagnetic field and a dispersion relation of a plane wave in TE and TM modes, determining a non-reflection condition which needs to be met by a PML; determining a parameter form of the PML meeting stable absorption performance by utilizing a relationship between an incident field and a transmission field in the lossy medium; constructing an improved wide-frequency-band PML parameter tensor according to a non-reflection condition, a parameter form meeting stable absorption performance and preset characteristic parameters; and realizing finite element forward modeling under a preset multi-level parallel architecture by using the improved wide-frequency-band PML parameter tensor to obtain a finite element forward modeling result. The method is accurate, efficient and wide in application range.

Description

Technical Field

[0001] The present application relates to the technical field of three-dimensional electromagnetic numerical simulation, and in particular to a three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML. Background Technique

[0002] In resource exploration, the electromagnetic differences of media are mainly utilized to explore the electrical structure underground or under the sea by measuring the variation laws of natural or artificial electromagnetic fields. Electromagnetic forward simulation provides technical support for inversion. Fast, efficient, and accurate forward simulation can provide reference for the interpretation of actual data, help understand geological / submarine structures (such as identifying hydrocarbon-bearing properties based on the huge resistivity difference between hydrocarbon-bearing reservoirs and the surrounding water-saturated formations), and reduce the risks and costs of onshore exploration and offshore drilling.

[0003] At present, 2D electromagnetic forward and inverse modeling are relatively mature, while 3D forward and inverse modeling are not yet perfect, and the response laws of complex geological bodies still need to be further understood. Therefore, high-efficiency and high-precision 3D forward modeling has become a research hotspot and difficulty. In 3D magnetotelluric forward modeling, the methods widely used at present mainly include the integral equation method, the finite difference method, and the finite element method, etc. The integral equation method only discretizes the abnormal body and is generally used for the calculation and analysis of the electromagnetic field of simple models. The finite difference method uses differences to approximate differentials and has certain difficulties in dealing with arbitrarily undulating terrain. The finite element method has the advantages of flexible mesh generation and the ability to accurately model complex boundaries, and can be used for models with complex physical property distributions or geometric features, and has become one of the mainstream methods for magnetotelluric forward modeling. The finite element method can be divided into nodal finite elements and vector finite elements. Vector finite elements can automatically satisfy the continuity of the tangential electric or magnetic field and the condition that the divergence is zero, thus avoiding the problem of spurious solutions that may occur in nodal finite elements. When using the finite element method for electromagnetic numerical simulation, whether using the total field method, the secondary field method or the potential method, it is necessary to truncate the open-domain space to reduce the computational requirements. To avoid the influence brought by truncation, it is usually necessary to extend the outer boundary by a sufficient distance (generally several skin depths) and then use Dirichlet, Neumann or mixed boundary conditions, which will greatly increase the computational resources and computational time consumed. Therefore, on the premise of ensuring the computational accuracy, improving the boundary conditions of magnetotelluric forward modeling and compressing the numerical simulation solution space are very helpful for improving the forward modeling efficiency. The Perfectly Matched Layer (PML) is a free-space simulation method proposed by Berenger in 1994 based on the field splitting theory. It is a specially designed virtual medium that can absorb electromagnetic waves without reflection. Berenger's PML is not perfect. The control equation used in the PML is a non-Maxwell equation, and the absorption effect on evanescent waves is not ideal enough. Chew and Weedon systematically analyzed Berenger's PML from the perspective of complex "coordinate stretching", thus providing a theoretical basis. Sacks et al. and Gedney proposed the Uniaxial Anisotropic Perfectly Matched Layer (UPML) based on Maxwell's equations, which does not require field splitting and can better absorb evanescent waves. Kuzuoglu and Mittra introduced a complex frequency shift stretching operator and proposed the Complex Frequency Shifted PML (CFS-PML), which can effectively absorb grazing waves and low-frequency waves.Roden and Gedney proposed the Convolutional Perfectly Matched Layer (CPML) by calculating the convolution terms appearing in the time-domain PML equations through recurrence formulas, which improved the computational efficiency. After continuous improvement, PML has been widely used in recent years. Some scholars have significantly improved the performance of the absorption boundary using PML based on the complex frequency-shifted stretching tensor. Some scholars have applied CFS-PML to the finite element time-domain simulation of ground penetrating radar (GPR) in Debye dispersive media, improving the simulation accuracy. Some scholars have used the staggered-grid finite-difference method with coupled PML boundary conditions to improve the computational accuracy and efficiency of three-dimensional magnetotelluric forward modeling. Some scholars have applied CFS-PML to high-order finite elements to solve three-dimensional GPR numerical simulation problems. Some scholars have applied PML to the model reduction algorithm, reducing the dimension of the Krylov subspace and improving the convergence effect. Some scholars have proposed GMI-PML, avoiding the problem of improper time-step synchronization and achieving higher absorption performance. UPML is derived from the high-frequency wave equation and has strong absorption ability in the low-frequency diffusion field, which will cause large numerical reflections. CFS-PML is also applicable in the low-frequency diffusion field, but the absorption effect is affected by the conductivity and frequency of the model, and the PML parameters need to be adjusted according to the model. To address these problems, some scholars have provided a wide-band PML metric, improving the truncation effect on the diffusion field and further expanding the application scenarios of PML.

[0004] However, the above technologies or methods have two drawbacks, namely, the PML is limited to a single conductivity parameter variation scenario and the calculation efficiency is relatively low. The current PML is limited to a single conductivity parameter variation scenario and does not consider the permeability parameter and more complex anisotropic scenarios. However, in practical scenarios, the magnetotelluric response is usually jointly affected by the conductivity and permeability of the medium. Many studies have shown that in high magnetic regions, the influence of permeability on the magnetotelluric response cannot be ignored. As the magnetic susceptibility increases, the electromagnetic forward response is significantly affected by it. The influence of the magnetic susceptibility parameter on the TM mode is stronger than that on the TE mode, and high-resistance anomalies are more susceptible to the influence of magnetic susceptibility than low-resistance anomalies. Some scholars have pointed out that electromagnetic data contains inherent information on the depth distribution of magnetic susceptibility. Others have used three-dimensional models of conductivity and permeability anisotropy to compare the influence of different magnetic susceptibilities on apparent resistivity. The traditional three-dimensional electromagnetic forward calculation has low efficiency. When performing electromagnetic forward modeling of anisotropic media based on vector finite elements, numerical simulation methods all need to truncate the infinite space to reduce the calculation requirements. In electromagnetic forward modeling, the inner boundary conditions are automatically satisfied, while for the outer boundary, to avoid the influence of truncation, boundary conditions are usually set after extending a sufficient distance. According to the type of boundary conditions, they can be divided into the first-kind boundary condition, the second-kind boundary condition, and the third-kind boundary condition. The first-kind boundary condition is also called the Dirichlet boundary condition, which directly imposes known values on the boundary. The second-kind boundary condition is also called the Neumann boundary condition, which gives the normal derivative value of the physical quantity on the boundary. The third-kind boundary condition is called the Robin boundary condition, and what is given on the boundary is a combination of the physical quantity and its normal derivative value, so it is also called the mixed boundary condition. Even for the third-kind boundary condition with relatively high calculation efficiency, the boundary still needs to be set far away from the anomaly, and the calculation scale is large. In addition, the problem scale of three-dimensional electromagnetic forward modeling is large and the calculation time is long. Using the direct solution method can obtain higher solution accuracy, but this will further exacerbate the occupation of system resources such as computer memory. The common method is to use multi-frequency point parallel solution to accelerate the solution speed, but the acceleration potential in aspects such as in-frequency parallelism and distributed storage has not been fully explored. Summary of the Invention

[0005] Based on this, it is necessary to provide a three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML that is accurate, efficient, and has a wide application range for the above technical problems.

[0006] A three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method coupled with PML, the method comprising: According to the continuity conditions of the electromagnetic field and the dispersion relations of plane waves in TE and TM modes, the non-reflection conditions that the PML needs to satisfy are determined; the parameter form for the PML to satisfy stable absorption performance is determined by using the relationship between the incident field and the transmitted field in a lossy medium; An improved wide-band PML parameter tensor is constructed based on the non-reflection conditions, the parameter form for satisfying stable absorption performance, and the preset characteristic parameters; The finite element forward modeling is implemented using the improved wide-band PML parameter tensor under the preset multi-level parallel architecture to obtain the finite element forward modeling results.

[0007] For the above three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method with coupled PML, in this application, the non-reflection conditions are determined based on the electromagnetic field continuity conditions and the TE and TM mode dispersion relations, and the parameter form is determined by combining the relationship between the incident field and the transmitted field in a lossy medium, so that the improved PML can consider the changes of both conductivity and permeability parameters simultaneously. In high magnetic regions, the permeability has a significant impact on the magnetotelluric response. The improved PML model is no longer limited to a single conductivity scenario and is applicable to more complex anisotropic models. At the same time, the parameter form determined in this application for satisfying stable absorption performance enables the improved wide-band PML to maintain stable absorption performance, overcomes the defects of traditional PML at low frequencies, and improves the overall absorption performance. The finite element forward modeling is implemented using the improved PML model under the multi-level parallel architecture, fully exploiting the acceleration potential such as in-band parallelism and distributed storage. Compared with traditional methods, the number of grids, calculation time, and memory occupancy are significantly reduced. By comprehensively constructing the improved PML parameter tensor based on the non-reflection conditions, absorption performance improvement, and preset characteristic parameters, its performance is more stable, and the accuracy and reliability of the calculation results can be ensured under different parameter changes and complex models. Description of the Drawings

[0008] Figure 1 It is a schematic flow chart of a three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method with coupled PML in an embodiment; Figure 2 It is a schematic diagram of plane wave incidence in an embodiment; Figure 3 It is a schematic diagram of the division of each region of the PML in a three-dimensional model in an embodiment; Figure 4 It is a schematic diagram of multi-level parallelism in another embodiment; Figure 5 It is a process diagram of the vector finite element solution for the coupled PML boundary conditions in an embodiment; Figure 6 It is a geometric schematic diagram of an experimental model in an embodiment; Figure 7A schematic diagram of apparent resistivity, phase and error calculated by model 1 using the improved PML and Dirichlet boundary conditions of the present application in one embodiment; Figure 8 A schematic diagram of apparent resistivity, phase and error calculated by using the PML and Dirichlet boundary conditions of the present application for model 2 in an embodiment; Figure 9 is an XOY cross-sectional diagram of apparent resistivity and phase in an XY mode and a YX mode in one embodiment; Figure 10 A schematic diagram of the parallel acceleration results of PML and the traditional method within a single frequency point in one embodiment; Figure 11 Schematic diagram of parallel acceleration of PML and traditional methods at 8 frequency points in one embodiment. DETAILED DESCRIPTION

[0009] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0010] In one embodiment, Figure 1 As shown, a three-dimensional magnetotelluric anisotropy parallel finite element forward modeling method coupled with PML is provided. The method can be used for scenarios with dual parameter changes of conductivity and permeability and three-dimensional electromagnetic forward modeling containing arbitrary anisotropic anomalies, and can significantly improve the truncation effect, including the following steps: Step 102, according to the continuity condition of the electromagnetic field and the dispersion relation of the plane wave in the TE and TM modes, determine the non-reflection condition that the PML needs to meet; and use the relationship between the incident field and the transmitted field in the lossy medium to determine the parameter form of the PML that satisfies the stable absorption performance.

[0011] like Figure 2 The figure shows the schematic diagram of a plane wave incident on the PML from the simulation area. To simplify the derivation process, in the diffuse field, the displacement current is ignored and the time harmonic factor is taken as , and assuming that PML is a uniaxial anisotropic medium, the anisotropic parameter and about Axis rotational symmetry, in tensor form (1) The incident surface is , the interface is , The area is the simulation area. The area of ​​the simulation is the PML. The conductivity of the simulation area is , the magnetic permeability is , the electrical conductivity and magnetic permeability of the PML are and . is the incident angle with respect to the normal direction, represents the incident wave vector, represents the reflected wave vector, represents the transmitted wave vector. For the convenience of derivation, the plane wave is represented in the form of the wave vector and the position vector, that is . In this way, the double curl equation can be represented by using the phase matching principle and the operator substitution relationship. To make and have non-zero solutions, it is necessary to satisfy that the coefficient matrix is not full rank. According to this condition, the dispersion relations of and modes are obtained. Combining the continuity conditions of the tangential components of the electric field and the magnetic field at the interface, the non-reflection conditions that the PML needs to satisfy are comprehensively solved. The purpose of using the PML is to truncate the calculation region as efficiently as possible. Therefore, on the basis of the first step, it is also necessary to further determine the form of the PML parameters so that it not only satisfies the non-reflection conditions but also has good absorption effects. Through the parameter relations determined by the non-reflection conditions and the relations between the electric field and the magnetic field in the lossy medium, the parameter form that satisfies the stable absorption performance is further solved as (2) This form is applicable to both high-frequency wave fields and low-frequency diffusion fields and can be used in a wide frequency band. Finally, the three-dimensional parameter matching matrix of the wide-band PML that meets the conditions can be obtained as (3) Among them, , , are the characteristic parameters of the PML. The finite element method discretizes the simulation space, which will cause the discontinuity on both sides of the interface between the simulation region and the PML, resulting in spurious wave reflection. Therefore, in the direction perpendicular to the interface, the absorption of the PML should increase gradually from zero. In this application , , parameters adopt a gradual exponential spatial distribution. When far from the interface, and values gradually increase, and value gradually decreases. Usually, the total thickness of the PML is limited and the number of layers is not too many. When the grid cell size of the PML is uniform, the spatial distributions of the characteristic parameters , , are as follows (4) Among them, is the permittivity in vacuum, is the thickness of the PML grid cell, represents the distance from the center of each PML cell to the interface between the simulation region and the PML, is the total thickness of the PML, is a constant. At the outer boundary of the PML, the perfect electric conductor boundary condition is adopted to terminate the entire calculation region. In Equation (4), , , are respectively , , 's maximum values. In order to control the absorption ability of the PML and balance the reflection error of the PEC boundary and the numerical discretization error, the optimal values of , , are denoted as , , . When the outer boundary of the PML is truncated by the PEC condition, the parameter is a function of the PML depth . Define the reflection coefficient of the PML as (5) Among them, is a gradual distribution, which can be written as a function with the coefficient , that is , where is the thickness of each layer of the PML. Write the integral in Equation (5) in discrete form (6) Among them is the total number of layers of the PML on each side, is the number of each layer of the PML. Equation (6) shows that the reflection coefficient of the PML is independent of the thickness and is related to the number of layers of the PML. This shows that the absorption performance of the PML is not sensitive to the PML thickness. As long as the number of PML layers remains unchanged, the thickness of the PML can be set as needed, so it can be applied to different scenarios. This application simultaneously considers the conductivity and permeability parameters, and the obtained PML can be used in scenarios with double-parameter changes. Since the non-reflection condition and attenuation condition of the PML are theoretically only related to the medium of the adjacent simulation region, it can be used for truncating the simulation region with arbitrary anisotropic anomalies inside.

[0012] Step 104: Construct an improved wide-band PML parameter tensor based on the non-reflection condition, the parameter form satisfying stable absorption performance, and the preset characteristic parameter spatial distribution.

[0013] The improved wide-band PML parameter tensor is no longer limited to a single conductivity scenario and is applicable to more complex anisotropic models. At the same time, the parameter form determined in this application that satisfies stable absorption performance enables the improved PML to maintain stable absorption performance in the low-frequency diffusion field, overcomes the defects of traditional PML at low frequencies, and improves the overall absorption performance.

[0014] Step 106: Use the improved wide-band PML parameter tensor to perform finite element forward modeling under the preset multi-level parallel architecture to obtain the finite element forward modeling result.

[0015] In a three-dimensional model, as Figure 3 shown, the PML region has a total of 6 plane regions, 12 edge regions, and 8 corner regions, and the parameter forms of each PML region are shown in Table 1.

[0016] Table 1

[0017] When calculating frequency points (as Figure 4 shown), the program runs within the total communication domain, using nodes and processes. Since the calculations for each frequency point do not interfere with each other, considering the sufficient parallelism within the frequency points, this application divides the total communication domain into sub-communication domains according to the number of frequency points. Each sub-communication domain contains frequency points to be calculated, and each frequency point is assigned processes. When calculating each frequency point, multiple processes run simultaneously to accelerate the matrix assembly, right-hand side construction, and matrix solution processes; the matrix is distributed and stored on different nodes responsible for calculating that frequency point to reduce the memory requirements of a single node. The petsc distributed toolkit is used for distributed processing in the matrix assembly and right-value construction parts, and the superlu-dist solver combined with petsc is used for distributed parallel solution in the matrix solution part.

[0018] The process of using the improved PML model to perform finite element forward modeling under the preset multi-level parallel architecture is as Figure 5 shown. When truncating the boundary in an open-domain problem, applying a perfectly matched layer to absorb the outgoing traveling wave will have a better effect. Therefore, here the calculation of the total electric field is decomposed into the independent calculations of the primary field and the secondary field . The relationship is (7) In the total field, the Maxwell equations and the double curl equation of the electric field are respectively (8) (9) where is the relative magnetic permeability of the simulation region, is the conductivity of the simulation region, is the permittivity of the simulation region, is the displacement current, is the identity matrix. In the primary field, the Maxwell equations (governing equations) in the simulation region and the PML region are (10) where is the background conductivity, is the background relative magnetic permeability. The double curl equation of the electric field can be expressed as (11) From equations (7), (9) and (11), the double curl equation of the secondary electric field can be obtained as (12) In the simulation region, is equal to the identity matrix; in the PML region, takes the form shown in Table 1. Integrating the entire region using the Galerkin weighted residual method, and using vector identities and the divergence theorem, finally we can obtain (13) where is the variation of , is the element coefficient matrix, . The corresponding global matrix equation is (14) where is the total system coefficient matrix, . After solving for the electric field according to equation (14), the magnetic field can be further obtained.

[0019] The field sources used in the magnetotelluric method of this application are and two orthogonal sources. Source is , , ; source is , , . Source The calculated electric and magnetic fields are , , , ; source The calculated electric and magnetic fields are , , , . Then the tensor impedance calculation formula is (15) From equation (15), the apparent resistivity and phase in different modes can be obtained as (16) For the vector finite element with coupled PML boundary conditions, only the PML constraint needs to be added to the control equation, and the PML parameter matching matrix is multiplied when establishing the element stiffness matrix. The other parts are the same as the ordinary vector finite element method.

[0020] For the above three-dimensional magnetotelluric anisotropic parallel finite element forward modeling method with coupled PML, in this application, the reflectionless condition is determined based on the electromagnetic field continuity condition and the TE and TM mode dispersion relations, and the parameter form is determined by combining the relationship between the incident field and the transmitted field in the lossy medium, so that the improved PML can consider the changes of both conductivity and permeability parameters simultaneously. In the high magnetic region, the permeability has a significant impact on the magnetotelluric response. The improved PML model is no longer limited to a single conductivity scenario and is applicable to more complex anisotropic models. At the same time, the parameter form determined in this application to satisfy the stable absorption performance enables the improved PML to maintain stable absorption performance in the low-frequency diffusion field, overcomes the defects of traditional PML at low frequencies, and improves the overall absorption performance. The finite element forward modeling is realized using the improved wide-band PML model under the multi-level parallel architecture, fully exploiting the acceleration potential such as in-band parallelism and distributed storage. Compared with the traditional method, the number of grids, calculation time, and memory occupancy are significantly reduced. By comprehensively constructing the improved PML based on the reflectionless condition, absorption performance improvement, and preset characteristic parameters, its performance is more stable, and the accuracy and reliability of the calculation results can be ensured under different parameter changes and complex models.

[0021] In one embodiment, according to the electromagnetic field continuity condition and the dispersion relations of plane waves in TE and TM modes, the reflectionless condition that the PML needs to satisfy is determined, including: The plane wave is represented in the form of a wave vector and a position vector, that is , using the phase matching principle and the operator The substitution relationship represents the double curl equation. To make and If there is a non-zero solution, it is necessary to satisfy that the coefficient matrix is rank-deficient. Based on the condition that the coefficient matrix is rank-deficient, the dispersion relations in the TE and TM modes are obtained. Then, combined with the continuity conditions of the tangential components of the electric and magnetic fields at the interface, the non-reflection conditions that the PML needs to satisfy are comprehensively solved.

[0022] In one embodiment, the non-reflection conditions that the PML needs to satisfy are

[0023] where represents the anisotropic parameter tensor, represents the parameter, represents , and directions.

[0024] In one embodiment, the parameter form for the PML to satisfy stable absorption performance in the low-frequency diffusion field is determined by using the relationship between the incident field and the transmitted field in a lossy medium, including: The parameter form for the PML to satisfy stable absorption performance is determined by using the relationship between the incident field and the transmitted field in a lossy medium as

[0025] In the low-frequency diffusion field, it is simplified to

[0026] where , , represent different characteristic parameters of the PML, is the conductivity, is the permeability, represents the angular frequency.

[0027] In one embodiment, the preset characteristic parameters are:

[0028] where is the permittivity in vacuum, is the thickness of the PML grid cell, represents the distance from the center of each PML cell to the interface between the simulation region and the PML, is the total thickness of the PML, is a constant, , , are respectively , , the maximum values of, represents the permeability.

[0029] In one of the embodiments, the pre-set multi-level parallel architecture includes: When calculating frequency points, within the total communication domain, nodes and processes are used to divide the total communication domain into sub-communication domains according to the number of frequency points. Each sub-communication domain contains frequency points to be calculated, and each frequency point is assigned processes. When calculating each frequency point, multiple processes run simultaneously to accelerate the matrix assembly, right-hand side term construction, and matrix solution processes; the matrix is distributed and stored on different nodes responsible for calculating this frequency point. In the matrix assembly and right-hand side term construction parts, the petsc distributed toolkit is used for distributed processing, and in the matrix solution part, the superlu-dist solver is combined with petsc for distributed parallel solution.

[0030] In a specific embodiment, the air-ground model 1, the air-sea-ground model 2, and the different magnetic permeability model group 3 are used. All numerical experiments are completed on a supercomputer. As Figure 6 shown, the geometric size of the simulation area of the model is 2000 m × 2000 m × 2000 m, the background resistivity is 100 , and the background relative magnetic permeability is 1. The anomaly is located at the center of the simulation area, with a top burial depth of 600 m and a size of 800 m × 400 m × 400 m. The simulation area and the anomaly are evenly divided into hexahedral units of 100 m × 100 m × 100 m. When using PML, 6 layers of PML are set on the outside of the simulation area, with each layer having a thickness of 0.05 m, and the outermost layer is truncated using the perfect electric conductor condition. The parameter settings are all , , ; when using the traditional grid extension method, the simulation area is extended by about 5 skin depths on the outside (at this extension distance, the secondary field at the outer boundary decays to less than 1%), and the size growth rate of the extended grid is 1.25. How to find the optimal parameters of PML is prior art and will not be elaborated here. Figure 6 is a schematic diagram of the model geometry. The blue wireframe is the calculation area, the purple wireframe is the simulation area, the area between them is the PML or extension area, and the middle green part is the anomaly.

[0031] Model 1 is a simplified geological model. The top 4 layers of the background medium in the simulation area are air layers, with a resistivity of , and a relative magnetic permeability of 1. The anomaly is anisotropic, with the principal axis resistivities being , , , and the Euler rotation angles being , , ; The principal axis magnetic susceptibilities (principal axis magnetic permeability = principal axis magnetic susceptibility + 1) are 0.3, 0.2, and 0.1 respectively, and the Euler rotation angles are , , respectively. The test frequency is 0.125 Hz. Figure 7 shows the apparent resistivity, phase, and error calculated by using the improved PML of this application and the Dirichlet boundary condition respectively. The test frequency is 0.125 Hz. ResXY-PML represents the calculation result of the PML of this application, ResXY-Diri represents the calculation result of the Dirichlet boundary condition (reference value), RE represents the relative error (‰), represents the absolute error ( °), the first row shows the calculation result of the ground at Y = 50 m along the X-direction survey line, and the second to fourth rows are schematic diagrams of the XOY cross-section Model 2 is a simplified layered ocean model. Different from Model 1, there are two seawater layers under the air layer of the background medium in the simulation area, and the resistivity is , and the relative magnetic permeability is 1. The anomaly is anisotropic, and the principal axis resistivities are , , respectively, and the Euler rotation angles are , , respectively; the principal axis magnetic susceptibilities are 0.2, 0.1, and 0.3 respectively, and the Euler rotation angles are , , respectively. The test frequency is 0.1 Hz. Figure 8 shows the apparent resistivity, phase, and error calculated by using the PML of this application and the Dirichlet boundary condition respectively. The test frequency is 0.1 Hz. ResXY-PML represents the calculation result of the PML of this application, ResXY-Diri represents the calculation result of the Dirichlet boundary condition (reference value), RE represents the relative error (‰), represents the absolute error ( °), the first row shows the calculation result of the seabed at Y = 50 m along the X-direction survey line, and the second to fourth rows are schematic diagrams of the XOY cross-section The model group 3 is a set of simplified geological models used to test the effects caused by different magnetic susceptibilities, in order to illustrate the significance of the PML of the present application that considers both conductivity and permeability parameters compared to the PML that only considered conductivity parameters in the past. They are only different from model 1 in the magnetic susceptibility of the anomaly body, and the rest of the settings are the same. Moreover, the same solution method with the PML of the present application as the boundary condition is adopted. Here, according to the actual situation of the natural world, four settings of negative magnetic susceptibility, zero magnetic susceptibility, low magnetic susceptibility, and high magnetic susceptibility are considered for comparison. The main axis magnetic susceptibilities of the corresponding anomaly bodies are: -0.1, -0.2, -0.3; 0, 0, 0; 0.3, 0.2, 0.1; 0.7, 0.6, 0.5. Figure 9 It shows the XOY cross-sectional diagrams of the apparent resistivity and phase in their XY mode and YX mode. Figure 9 Model group 4, XOY cross-sectional diagrams of the apparent resistivity and phase under different permeabilities. The test frequency is 0.125 Hz. ResXY represents the apparent resistivity in the XY mode, PhaseXY represents the phase in the XY mode, and ms = -0.1 -0.2 -0.3 means the main axis magnetic susceptibilities of the anomaly body are -0.1, -0.2, and -0.3 respectively. In order to fully test the acceleration effects brought by multi-level parallelism and the PML boundary condition, two groups of tests of single frequency point and 8 frequency points are carried out here. The calculation model is as Figure 6 shown. In the traditional method, the number of unknowns is 888822, and the test frequency is 0.125 Hz. Figure 10 shows the parallel acceleration results of the PML and the traditional method within a single frequency point. The test frequency is 0.125 Hz. The blue solid line with circles represents the parallel acceleration ratio of the PML method of the present application, the black solid line with triangles represents the parallel acceleration ratio of the traditional method, and the red dashed line with squares represents the total acceleration ratio of the parallel PML method compared to the serial traditional method. Figure 11 shows the parallel acceleration of the PML and the traditional method at 8 frequency points. The test frequency is 0.125 Hz. The blue solid line with circles represents the parallel acceleration ratio of the PML method of the present application, the black solid line with triangles represents the parallel acceleration ratio of the traditional method, and the red dashed line with squares represents the total acceleration ratio of the parallel PML method compared to the serial traditional method. Table 2 shows the performance comparison between the PML of the present application and the traditional method. Table 3 shows the comparison of parallel calculation times at a single frequency point. Table 4 shows the comparison of parallel calculation times at 8 frequency points.

[0032] Table 2

[0033] Table 3

[0034] Table 4

[0035] At Figure 7 and8 Among them, (a) and (b) respectively show the apparent resistivity in the XY mode and the YX mode, and (c) and (d) respectively show the phase in the XY mode and the YX mode. The first row of the figure shows the apparent resistivity and phase curves of the survey line at the ground and at the seabed where Y = 50 m. The second and third rows show the XOY cross-sectional views of the ground apparent resistivity and phase. The fourth row shows the relative error of the apparent resistivity and the absolute error of the phase. The results show that the calculation results of the two are in very good agreement. In Figure 7, the maximum relative error of the apparent resistivity in the two modes is less than 0.2‰, and the absolute error of the phase ; in Figure 8, the maximum relative error of the apparent resistivity in the two modes is less than 0.03‰, and the absolute error of the phase . This indicates that the PML of the present application is applicable to the forward modeling of homogeneous media and anisotropic anomalies, and the forward modeling of layered media and anisotropic anomalies. And only 6 layers, a total of 0.3 m of PML absorption is required to achieve almost the same calculation results as the first-type boundary condition using the large-distance continuation method, with very high calculation efficiency and calculation accuracy.

[0036] In Figure 9 , (a) to (d) respectively show the calculation results from negative magnetic susceptibility to high magnetic susceptibility. The first to fourth rows respectively show the apparent resistivity and phase in the XY mode and the YX mode. It can be observed from the figure that as the magnetic susceptibility increases, the color of the low-resistance dark area gradually fades, while the high-resistance yellow area becomes brighter. The phase cross-sectional view has the same changing trend. This indicates that it is very necessary to improve the PML that only considers the conductivity parameter in the past and propose a PML that considers the changes of both conductivity and magnetic permeability parameters.

[0037] The parallel test results of the two methods at a single frequency point and 8 frequency points are as Figure 10 and Figure 11 shown. The PML speedup curve represents the ratio of the single-process calculation time to the parallel calculation time of the PML method of the present application. The traditional speedup curve represents the ratio of the single-process calculation time to the parallel calculation time of the traditional method. The total speedup curve represents the ratio of the single-process calculation time of the traditional method to the parallel calculation time of the PML method of the present application, which is the acceleration effect brought by the combined application of the hierarchical parallelization design and the PML boundary condition of the present application. Figure 10The results show that the PML speedup reaches its maximum value of approximately 5.46 at 32 processes; the traditional speedup is approximately 6.73 at 128 processes; the total speedup reaches its maximum value of approximately 85.15 at 32 processes. When the number of processes is from 1 to 4, the acceleration effects brought by the hierarchical parallel design for the two methods are almost the same. When the number of processes is from 4 to 32, the speedups of both methods increase rapidly with the increase in the number of processes. However, the PML speedup curve of this application is above the traditional method's speedup curve. This may be because, among the additional system overheads brought by the parallelization design, the advantage of the PML method of this application in saving computing resources is transformed into additional parallel acceleration effects. When the number of processes is from 32 to 168, the speedup of the traditional method still has a certain increase, but the improvement effect brought by increasing the number of processes becomes weaker and weaker. While the PML method's speedup of this application no longer increases but begins to slightly decrease. At this time, the computing time of the PML method of this application has decreased from 106.2 seconds for a single process to 19.44 seconds for 32 processes. This is because the additional system overheads brought by dividing the calculation into more processes are excessive, resulting in a slight decrease in the speedup.

[0038] For the total speedup, it can be divided into three parts, namely, the boundary condition acceleration brought by the PML of this application, the parallel acceleration brought by the hierarchical parallelization design, and the overhead acceleration brought by the smaller parallel overhead of the PML method. When the number of processes is from 1 to 4, the "slope" of the total speedup curve is almost the same as that of the other two curves, indicating that the parallel acceleration effects are almost the same at this time, and the boundary condition acceleration causes the total speedup curve to be above the other two curves; when the number of processes is from 4 to 32, in addition to the boundary condition acceleration, the total speedup and the PML speedup have almost the same parallel acceleration and overhead acceleration. However, further observing the three curves, it can be found that the overhead acceleration effect of the total speedup and the PML speedup is the best when the number of processes is from 4 to 8, and the overhead acceleration effect gradually decreases when the number of processes is from 8 to 32; when the number of processes is from 32 to 128, for the same reason, the total speedup curve and the PML speedup curve have the same trend.

[0039] Figure 11The results show that: The PML speedup ratio basically reaches its maximum value at 512 processes, approximately 37.07; the traditional speedup ratio is approximately 33.62 at 512 processes; the total speedup ratio basically reaches its maximum value at 512 processes, approximately 649.63. When the number of processes is from 1 to 32, as the number of processes increases, the PML speedup ratio increases rapidly, while the traditional speedup ratio increases slowly. This may be because the traditional method requires more system resources, and with the current number of processes, compared with the parallel overhead, the acceleration potential cannot be released. When the number of processes is from 32 to 256, all three speedup ratio curves rise rapidly as the number of processes increases. When the number of processes is from 256 to 512, the PML speedup ratio and the total speedup ratio curves gradually slow down, which may be due to the same reason as the speedup ratio bottleneck occurring during single-frequency point parallel testing. In the figure, the PML speedup ratio curve of this application is always above the traditional method speedup ratio curve. This may be because as the number of calculation frequency points increases, more calculation resources are required, resulting in a more obvious acceleration effect of the overhead. It can be seen from Table 2 that in the 3 test models, the grid number ratio between the two exceeds 8 times, the time ratio exceeds 13 times, and the average memory occupancy ratio exceeds 12 times. This indicates that when the same truncation effect is achieved, compared with the traditional boundary conditions, the PML of this application has obvious advantages in terms of grid number, calculation time, and memory occupancy. It can be seen from Tables 3 and 4 that the hierarchical parallel design of this application has obvious acceleration effects for both single-frequency points and multi-frequency points. Compared with the traditional boundary conditions, the maximum speedup ratio of the parallelized PML of this application is approximately 85.24 at single-frequency points and approximately 649.63 at 8 frequency points. Further observation reveals that the PML of this application can achieve almost the maximum acceleration effect with a relatively small number of processes (for example, the acceleration effects at 32 processes and 128 processes are almost the same at single-frequency points, and the acceleration effects at 128 processes and 512 processes are almost the same at 8 frequency points), while the traditional boundary conditions require more processes to achieve the maximum acceleration effect (for example, more than 128 processes are required at single-frequency points, and more than 512 processes are required at 8 frequency points). This further demonstrates the advantage of the multi-level parallelized PML method of this application in saving system resources.

[0040] Previous PMLs were all derived based on conductivity parameters and did not consider the influence of permeability. Here, the present application proposes a three-dimensional magnetotelluric conductivity and permeability anisotropy parallel finite element forward modeling method with coupled PML boundary conditions. Different from the past, the present application considers the permeability parameter, proposes a PML applicable to the variation of both conductivity and permeability parameters, improves it for the low-frequency band, and applies it to the anisotropic model. In addition, in order to further improve the three-dimensional forward modeling efficiency, a multi-level parallel design is also carried out. Through numerical experiments, the following conclusions are obtained: 1. The PML method proposed in the present application is correct, efficient, and has high precision, and is applicable to homogeneous models and layered models with conductivity and permeability anisotropic anomalies; 2. Compared with the traditional boundary conditions, the number of PML grids in the present application is reduced by more than 85%, and the calculation time and memory occupation are reduced by more than 90%, improving the efficiency of three-dimensional electromagnetic forward modeling; 3. Compared with the traditional method with 888822 unknowns, through multi-level parallelization design, the acceleration ratio of the PML in the present application is about 85.24 when using 32 processes to solve a single frequency point, and about 649.63 when using 512 processes to solve 8 frequency points, and the maximum acceleration effect can be exerted with fewer process numbers; 4. The parameters of the PML in the present application have a wide range of use, and the same set of optimal parameters has good truncation effects in models with different sizes, different mesh divisions, and different test frequencies in the present application; 5. The change of magnetic susceptibility will affect the apparent resistivity and phase. As the magnetic susceptibility increases, the apparent resistivity and phase will also be affected and become larger, and even show completely different characteristics.

[0041] Compared with traditional boundary conditions, the Perfectly Matched Layer (PML) is a more efficient and accurate truncation method. However, the current PML is limited to the scenario of single conductivity parameter variation and cannot be applied to the scenarios of dual-parameter variation of conductivity and permeability and complex anisotropic models. Therefore, this application simultaneously considers the dual parameters of conductivity and permeability as well as anisotropy, proposes a PML applicable to complex anisotropic models, and further combines MPI to design a multi-level parallel scheme, realizing the 3D magnetotelluric conductivity and permeability anisotropic parallel vector finite element forward modeling with coupled PML boundary conditions. By comparing with the results of predecessors, it is verified that the PML boundary conditions proposed in this application have the advantages of high efficiency, high precision, and stable performance. The numerical experimental results of several models in this application show that compared with traditional boundary conditions, the number of PML grids in this application is reduced by more than 85%, and the calculation time and memory occupation are reduced by more than 90%; compared with the traditional method with 888822 unknowns, by applying the PML and multi-level parallelization design in this application, the problem of low calculation efficiency of 3D electromagnetic forward modeling is alleviated. When using 32 processes to solve a single frequency point, the speedup ratio is about 85.24, and when using 512 processes to solve 8 frequency points, the speedup ratio is about 649.63, and the maximum acceleration effect can be achieved with fewer process numbers. The PML proposed in this application has a wider scope of application and better effects, and has a broader application prospect.

[0042] It should be understood that although Figure 1 the steps in the flowchart of Figure 1 are shown in sequence according to the indication of the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless there is a clear indication in this application, the execution of these steps has no strict order limit, and these steps can be executed in other orders. Moreover,

[0043] The technical features of the above embodiments can be combined arbitrarily. For the sake of concise description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope recorded in this specification.

[0044] The above-described embodiments merely represent several implementation manners of the present application. The description thereof is relatively specific and detailed, but it should not be construed as a limitation on the scope of the invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several modifications and improvements can still be made, and these all fall within the protection scope of the present application. Therefore, the protection scope of the present application shall be subject to the appended claims.

Claims

1. A three-dimensional magnetotelluric anisotropy parallel finite element forward modeling method coupled with PML, characterized in that: The method comprises: According to the continuity conditions of the electromagnetic field and the dispersion relations of the plane wave in the TE and TM modes, the non-reflection conditions that the PML needs to meet are determined; the relationship between the incident field and the transmitted field in the lossy medium is used to determine the parameter form of the PML to meet stable absorption performance; Determine an improved broadband PML parameter tensor according to the non-reflection condition, the parameter form satisfying the stable absorption performance, and the preset characteristic parameter spatial distribution; The improved wide-band PML parameter tensor is used to implement finite element forward modeling under a pre-set multi-level parallel architecture to obtain a finite element forward modeling result.

2. The method according to claim 1, characterized in that According to the continuity conditions of the electromagnetic field and the dispersion relations of the plane wave in the TE and TM modes, the non-reflection conditions that the PML needs to meet are determined, including: The representation of a plane wave takes the form of a wave vector and a position vector, i.e. , using the phase matching principle and operator Substituting the relation to express the double curl equation, we need to make and If there is a non-zero solution, it is necessary to satisfy the condition that the coefficient matrix is ​​not rank-full. Based on the condition that the coefficient matrix is ​​not rank-full, the dispersion relation in TE and TM modes is obtained, and then combined with the continuity condition of the tangential components of the electric and magnetic fields at the interface, the no-reflection condition that PML needs to satisfy is comprehensively solved.

3. The method according to claim 2, characterized in that The reflection-free condition that needs to be satisfied in solving the PML is in, represents the anisotropy parameter tensor, Indicates the parameters, express , and direction.

4. The method according to claim 3, characterized in that The relationship between the incident field and the transmitted field in the lossy medium is used to determine the parameter form of the PML that satisfies stable absorption performance, including: Using the relationship between the incident field and the transmitted field in the lossy medium, the parameter form of the PML that satisfies the stable absorption performance is determined as follows: in, , , Represents different characteristic parameters of PML, is the conductivity, is the magnetic permeability, Represents the angular frequency.

5. The method according to claim 4, characterized in that The pre-set characteristic parameters are: in, is the dielectric constant in vacuum, is the thickness of the PML grid cell, represents the distance from the center of each PML unit to the simulation area-PML interface, is the total thickness of the PML, is a constant, , , They are , , The maximum value of Represents magnetic permeability.

6. The method according to claim 1, characterized in that The pre-set multi-level parallel architecture includes: In calculation When the frequency point is within the total communication domain, use nodes and process, the total communication domain is divided into sub-communication domains, each of which contains frequency points to be calculated, each frequency point is assigned When calculating each frequency point, multiple processes run simultaneously to accelerate the matrix assembly, right-hand term construction, and matrix solution process; the matrix is ​​distributed and stored on different nodes responsible for calculating the frequency point. The petsc distributed toolkit is used for distributed processing in the matrix assembly and right-hand term construction parts, and the superlu-dist solver combined with petsc is used for distributed parallel solution in the matrix solution part.

Images