Aviation natural source electromagnetic method complex medium parallel three-dimensional forward modeling method

Abstract

The invention belongs to the field of geophysical exploration, and particularly relates to an aviation natural source electromagnetic method complex medium parallel three-dimensional forward modeling method. Comprising the following steps: S1, constructing a three-dimensional geological geometric model of a research area; s2, generating a model unstructured tetrahedral mesh; s3, setting anisotropic conductivity tensor and magnetic conductivity of the model unstructured tetrahedral mesh; s4, constructing a frequency domain electric field control equation; s5, adopting a vector finite element discrete electric field control equation; s6, setting boundary conditions; s7, solving a finite element linear equation; s8, calculating three components of an electric field and a magnetic field at the measuring point and the base station; and S9, calculating impedance, apparent resistivity, phase and tipper data at the measuring point. According to the method, a three-dimensional forward modeling algorithm is provided for the aviation natural field source electromagnetic method, a basis is provided for research of the aviation natural field source electromagnetic method theory and inversion interpretation of actual observation data, and then the aviation natural field source electromagnetic method can better serve the application fields of mineral exploration, engineering exploration and the like.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration, specifically relating to a three-dimensional forward modeling method for electromagnetic methods of airborne natural field sources applicable to complex geological conditions. Background Technology

[0002] In areas with complex geological environments such as plateaus, mountains, glaciers, and swamps, ground exploration is difficult. Airborne natural field electromagnetic methods, which utilize natural orthogonal field sources to electromagnetically detect underground electrical structures, are an emerging geophysical exploration technique. Airborne natural field electromagnetic methods typically employ helicopters or unmanned aerial vehicles (UAVs) carrying magnetic sensors to observe different components of the magnetic field from the air, while simultaneously deploying a small number of ground stations to observe electric or magnetic field components. Due to its use of airborne data acquisition, this method offers advantages in efficiency and convenience, making it highly suitable for complex geological environments and large-scale exploration, and it has enormous development potential in fields such as mineral exploration and engineering surveys. The three-dimensional forward modeling method of the airborne natural field electromagnetic method is the basis of three-dimensional inversion. However, the current research on the three-dimensional forward modeling method of the airborne natural field electromagnetic method has the following shortcomings: (1) The airborne natural field electromagnetic method is different from the traditional ground electromagnetic observation method. The existing ground magnetotelluric three-dimensional forward modeling simulation method cannot be directly used for the simulation of the airborne natural source electromagnetic method; (2) The application area of ​​the airborne natural source electromagnetic method usually has complex terrain and geological structure, but there is currently no airborne natural source electromagnetic method three-dimensional forward modeling simulation method specifically for complex geological structures; (3) Anisotropic and high magnetic permeability media are widely present in rock strata. The existing airborne natural source electromagnetic three-dimensional forward modeling method does not consider these two phenomena, thus affecting the calculation accuracy of the airborne natural source electromagnetic forward modeling simulation.

[0003] Chinese patent application CN118884546A discloses an airborne electromagnetic forward modeling method and medium based on a physical information neural network, comprising: determining a three-dimensional model of the airborne electromagnetic exploration area; delineating multiple airborne electromagnetic measurement points based on the three-dimensional model; constructing a physical information neural network; training and optimizing the physical information neural network with the goal of minimizing the mean square error of the loss function, using the three-dimensional coordinates of different airborne electromagnetic measurement points, the secondary electric field values ​​of boundary measurement points, and airborne electromagnetic system parameters as inputs; inputting the three-dimensional coordinates of the target airborne electromagnetic measurement point into the trained and optimized physical information neural network for frequency domain airborne electromagnetic three-dimensional forward modeling to obtain the secondary electric field value of the target airborne electromagnetic measurement point; and calculating the secondary magnetic field value of the target airborne electromagnetic measurement point using Faraday's law of electromagnetic induction. The method proposed in patent CN118884546A is applicable to controlled-source airborne frequency domain electromagnetic methods but not to airborne natural-source electromagnetic methods.

[0004] Chinese patent application CN114880902A discloses a time-domain vector finite element forward modeling method and apparatus for transient electromagnetic methods using magnetic sources. The method includes: establishing a simulated geological model; applying an excitation source to the meshed simulated geological model to obtain the Maxwell's equations and the electric field control equations; discretizing the control equations using the weighted residual method and selecting the first type of Whitney basis functions to establish the finite element equations; using a time-difference scheme to discretize the finite element equations over time; assembling the overall stiffness matrix and using the PADISO solver to solve a large-scale linear sparse equation set; determining the iteration time step based on the current-off characteristics of the excitation source, then iterating with a variable time step to obtain the electric field components, and obtaining the electric field values ​​at each node through the basis functions; and obtaining the transient electromagnetic observations and apparent resistivity parameters according to Faraday's law of electromagnetic induction. The method disclosed in CN114880902A is applicable to transient electromagnetic methods using magnetic sources, but not to electromagnetic methods using airborne natural sources. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and provide a parallel three-dimensional forward modeling method for complex media based on unstructured tetrahedral mesh modeling technology and vector finite element method using airborne natural source electromagnetic method. This method makes up for the lack of a unified solution to the problem of the influence of complex geological structure, conductivity anisotropy and magnetic permeability in airborne natural source electromagnetic forward modeling, thereby obtaining the real airborne natural field source electromagnetic response of complex Earth media.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: A parallel three-dimensional forward modeling method for complex media using electromagnetic methods from natural sources in airborne applications, comprising the following steps: S1: Construct a three-dimensional geological geometric model of the study area; Based on the implicit modeling method of CAD, construct a three-dimensional geological geometric model of the study area using topographic, stratigraphic and geological body data; S2: Generate an unstructured tetrahedral mesh for the model; use an unstructured tetrahedral mesh to partition the geometric model, locally refine the mesh at ground base stations and aerial measurement points, and use a coarse mesh outside the core area to generate an unstructured tetrahedral mesh for the geological model of the entire study area. S3: Setting the anisotropic electrical conductivity tensor and magnetic permeability of the unstructured tetrahedral mesh of the model; Based on known geological, borehole and geophysical information, assign arbitrary anisotropic electrical conductivity tensors and magnetic permeability to the tetrahedral mesh elements of the model.

[0007] Construct the frequency domain electric field control equations; after eliminating the magnetic field term from the Maxwell's equations in the frequency domain, the following vector electric field control equation formula (1) can be obtained: Formula (1): , Where, i= ω represents angular frequency. For electric field, Permeability, For conductivity tensor, , , , , , , , , and These are the conductivity components in each direction; S5: The electric field control equation is discretized using vector finite element method; the frequency domain electric field control equation (1) is processed using Galerkin finite element analysis to obtain the weak integral form of the differential equation: Formula (2):

[0008] Where Ω represents the simulated region, For vector interpolation basis functions, Let represent the electric field to be solved within the study area.

[0009] The electric field within an unstructured tetrahedron is discretized using first-order linear vector basis functions. It can be determined by the electric field at the center of each edge. definition: Formula (3): , in Let represent the e-th tetrahedral element, and j represent the local number of the edge inside the tetrahedron.

[0010] Using the first vector Green's theorem, the discrete form of equation (2) is obtained in each tetrahedral element. A sparse set of finite element coefficient equations is assembled using the unstructured tetrahedral mesh of the model, electrical conductivity, and magnetic permeability parameters: Formula (4): , in, These are the right-hand terms of the equation after applying boundary conditions. and These are the stiffness matrix and the mass matrix, respectively, and their specific discrete forms are as follows: Formula (5): , Formula (6): , in, and Let these represent the stiffness matrix and the mass matrix, respectively. For vector basis functions, For the cell region to be solved, It is a volume infinitesimal element. and Let represent the permeability and resistivity tensors of the element, respectively. The superscript "e" indicates the e-th tetrahedral element, and the subscripts "j" and "k" indicate the edge number of the element. The integral of the above elements is calculated using analytical methods.

[0011] S6: Set boundary conditions; set Dirichlet boundary conditions. Apply an excitation source to the top of the model. and The electric fields at the lateral and bottom boundaries of the model can be obtained through one-dimensional forward modeling. Let the center coordinates of the edges at the boundary of the three-dimensional model be... The electric field at the center of the edge can be obtained through a one-dimensional forward modeling algorithm. Project the one-dimensional analytical solution onto the center of the edge of the three-dimensional model boundary: Formula (7): , in, It is the tangential unit vector of the boundary edge of the 3D model. It is the electric field at the center of the edge. It is the electric field component obtained by one-dimensional forward modeling.

[0012] By applying boundary conditions using the "01" assignment method, a sparse system of linear finite element equations to be solved is obtained: Formula (8): , in, and Finite element coefficient matrices for two polarization modes The same, but the source direction on the right-hand side of the equation is... different.

[0013] S7: Solve the finite element linear equations; parallel computation is used for different computation frequencies, and the direct solver PARDIASO is used to solve the large-scale finite element linear equation system for each computation frequency. First, the coefficient matrix is... Matrix decomposition and two back-substitutions yield solutions for the two polarization modes. and .

[0014] S8: Calculate the three components of the electric and magnetic fields at the measuring point and the base station; solve for the electric and magnetic field components at the measuring point, where the measuring point is at the [missing information]. e When the electric field at the measuring point is within a tetrahedral element... It can be represented by vector shape functions as: Formula (9):

[0015] According to Faraday's law, the magnetic field component at the measuring point... It can be expressed by the electric field and vector shape function as follows: Formula (10): .

[0016] S9: Calculate the impedance, apparent resistivity, phase, and tilt vector at the measurement point; the response data for the airborne natural field source electromagnetic method includes the impedance tensor, apparent resistivity, phase, and tilt vector. The calculation expressions for each component of the impedance tensor in the airborne natural field source electromagnetic method are as follows: Formula (11): , , , , in, , , , This refers to the horizontal component of the observed magnetic field from aerial observations. , , , For the observation of the horizontal component of the electric field from a ground-based fixed base station, the subscripts “x” and “y” represent two horizontal directions, and the subscripts “1” and “2” represent two polarization modes.

[0017] After calculating the impedance tensor, the apparent resistivity and phase can be obtained, i.e.: Formula (12):

[0018] Formula (13):

[0019] The expression for calculating the tilt vector is: Formula (14): , , When only the vertical component of the magnetic field is observed from the air... , This represents the vertical component of the airborne magnetic field. , , , For ground-based fixed base stations, the horizontal component of the magnetic field is observed; when observing the three components of the magnetic field from the air... , , , , , All are magnetic field components observed from the air. Subscripts "x", "y", and "z" represent the directions of the three components, and subscripts "1" and "2" represent the two polarization modes.

[0020] Preferably, the vector interpolation basis function described in step S4 Represented using first-order linear basis functions: Where L is the node basis function, l is the length of the edge, and the subscript " "and" "Represents the start and end points of the edge."

[0021] Preferably, the parallel computing method described in step S6 includes two layers of parallelism. The first layer is to allocate a computing process for the calculation of the finite element coefficient matrix for each frequency, and the second layer is to use the OpenMP thread inside the PARDIASO solver for parallelism.

[0022] This invention proposes a three-dimensional forward modeling method using airborne natural source electromagnetic methods applicable to complex geological structures. It overcomes the limitations of existing airborne natural source electromagnetic forward modeling techniques that cannot comprehensively consider the influence of complex geological structures, anisotropic media, and magnetic permeability. Its innovations and advantages are as follows: 1. This invention can accurately simulate the influence of anisotropic electrical conductivity and magnetic permeability on electromagnetic data of airborne natural sources, and is applicable to the influence of anisotropic and magnetic media that are widely present in real Earth media.

[0023] 2. This invention employs three-dimensional geological modeling technology and a vector finite element method based on unstructured tetrahedral meshes to accurately simulate the electromagnetic response of aerospace natural sources with complex geological structures.

[0024] 3. This invention uses parallel technology to calculate observation data at different frequencies, resulting in high calculation speed.

[0025] This invention provides a three-dimensional forward modeling algorithm for the airborne natural field source electromagnetic method, laying the foundation for the theoretical research of the airborne natural field source electromagnetic method and the inversion interpretation of actual observation data. This will enable the airborne natural field source electromagnetic method to better serve applications such as mineral exploration and engineering exploration. Attached Figure Description

[0026] Figure 1 Flowchart of the three-dimensional forward modeling method for electromagnetic methods of airborne natural sources; Figure 2 Unstructured mesh and electrical parameters of layered models; Figure 3 Comparison of forward modeling results and analytical solutions (apparent resistivity-frequency) of the present invention using the electromagnetic method of airborne natural sources for layered models, with ground base station location (0, 0, 0) m and airborne measuring point location (0, 0, -100) m; Figure 4Comparison of forward modeling results and analytical solutions (phase-frequency) of the layered model using the electromagnetic method of the airborne natural source of the present invention; ground base station location (0, 0, 0) m; airborne measuring point location (0, 0, -100) m. Figure 5 The geometric structure of the topographic block anomaly model, with black dots representing the locations of airborne magnetic field measurement points, which are 200m above the ground; Figure 6 Unstructured mesh and electrical parameters of a model with topographic block anomalies; Figure 7 Apparent resistivity forward modeling results of the airborne natural source electromagnetic method xy mode with topographic block anomaly model, ground base station location (0, -5000, 0) m, Y=0 m survey line.

[0027] Figure 8 Apparent resistivity forward modeling results of the electromagnetic method yx model with topographic block anomaly model, ground base station location (0, -5000, 0) m, Y=0 m survey line.

[0028] Figure 9 Phase rate forward modeling results of the XY mode of the electromagnetic method for airborne natural source with topographic block anomaly model, ground base station location (0, -5000, 0) m, Y=0 m survey line.

[0029] Figure 10 Phase forward modeling results of the airborne natural source electromagnetic method yx model with topographic block anomaly model, ground base station location (0, -5000, 0) m, Y=0m survey line. Detailed Implementation

[0030] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0031] The structures, proportions, and sizes illustrated in the accompanying drawings are merely for illustrative purposes and to aid those skilled in the art in understanding and reading the invention. They are not intended to limit the scope of the invention and therefore have no substantial technical significance. Any modifications to the structure, changes in proportions, or adjustments to size, provided they do not affect the effectiveness or purpose of the invention, should still fall within the scope of the technical content disclosed herein. Furthermore, the terms "upper," "lower," "left," "right," "middle," and "one" used in this specification are merely for clarity and not intended to limit the scope of the invention. Changes or adjustments to their relative relationships, without substantially altering the technical content, should also be considered within the scope of the invention's implementation.

[0032] See Figure 1This embodiment discloses a flowchart of a preferred three-dimensional forward modeling method for electromagnetic methods using airborne natural sources. To better illustrate the advantages and objectives of this embodiment, examples are provided: Example 1

[0033] like Figure 1 As shown, the parallel three-dimensional forward modeling method for complex media using the electromagnetic method of airborne natural sources includes the following steps: S1: Construct a geometric model based on the implicit modeling method of CAD. The overall range of the model is 5km×5km×6km, the air layer is 1km×5km×5km, the underground uniform half space is 5km×5km×5km, the first layer is 200m thick, and the second layer is 300m thick. S2: The geometric model is meshed using an unstructured tetrahedral mesh. Local refinement is performed at the ground base station coordinates (0,0,0) m and the aerial measuring point coordinates (0,0,-100) m, with the mesh size set to 1 m. A coarse mesh is used for meshing outside the core area, thereby generating an unstructured tetrahedral mesh for the geological model of the entire study area. The number of tetrahedral elements in the model is 99326, and the number of edges is 119118. S3: Based on known geological, borehole, and geophysical information, assign arbitrary anisotropic electrical conductivity tensors and magnetic permeability to the tetrahedral mesh elements of the model. The air resistivity of the model is set to 10. 8 Oh mm, relative permeability is 1; the resistivity of the first layer principal axis is (10 4 , 10 3 , 10 4 The first layer has a resistivity of (20, 10, 20) ohm and a relative permeability of 1; the second layer has a resistivity of (20, 10, 20) ohm and a relative permeability of 2; the third layer has a resistivity of 100 ohm and a relative permeability of 1.

[0034] S4: After eliminating the magnetic field term from Maxwell's equations in the frequency domain, the following vector electric field control equation formula (1) can be obtained. Formula (1): , Where, i= ω represents angular frequency, and the calculated frequency is 10. 1 ~10 3 Take 11 frequency points with equal logarithmic intervals. For electric field, Permeability, For conductivity tensor, , , , , , , , , and These represent the conductivity components in each direction.

[0035] S5: Using Galerkin finite element analysis, the frequency domain electric field control equation (1) is processed to obtain the weak integral form of the differential equation: Formula (2): , Where Ω represents the simulated region, For vector interpolation basis functions, Let represent the electric field to be solved within the study area.

[0036] Vector interpolation basis function Represented using first-order linear basis functions: Where L is the node basis function, l is the length of the edge, and the subscript " "and" "Represents the start and end points of the edge."

[0037] The electric field within a tetrahedral element is discretized using first-order linear vector basis functions. It can be determined by the electric field at the center of each edge. definition: Formula (3): , in Let represent the e-th tetrahedral element, and j represent the local number of the edge inside the tetrahedron.

[0038] Using the first vector Green's theorem, the discrete form of equation (2) is obtained in each tetrahedral element. A sparse set of finite element coefficient equations is assembled using the unstructured tetrahedral mesh of the model, electrical conductivity, and magnetic permeability parameters: Formula (4): , in, These are the right-hand terms of the equation after applying boundary conditions. and These are the stiffness matrix and the mass matrix, respectively, and their specific discrete forms are as follows: Formula (5): , Formula (6): , in, and Let these represent the stiffness matrix and the mass matrix, respectively. For vector basis functions, For the cell region to be solved, It is a volume infinitesimal element. and Let represent the permeability and resistivity tensors of the element, respectively. The superscript "e" indicates the e-th tetrahedral element, and the subscripts "j" and "k" indicate the edge number of the element. The integral of the above elements is calculated using analytical methods.

[0039] S6: Set Dirichlet boundary conditions. Apply an excitation source at the top of the model. and The electric fields at the lateral and bottom boundaries of the model can be obtained through one-dimensional forward modeling. Let the center coordinates of the edges at the boundary of the three-dimensional model be... The electric field at the center of the edge can be obtained through a one-dimensional forward modeling algorithm. Project the one-dimensional analytical solution onto the center of the edge of the three-dimensional model boundary: Formula (7): , in, It is the tangential unit vector of the boundary edge of the 3D model. It is the electric field at the center of the edge. It is the electric field component obtained by one-dimensional forward modeling.

[0040] By applying boundary conditions using the "01" assignment method, a sparse system of linear finite element equations to be solved is obtained: Formula (8): , in, and Finite element coefficient matrices for two polarization modes The same, but the source direction on the right-hand side of the equation is... different.

[0041] S7: Parallel computation is used for different computation frequencies. The direct solver PARDIASO is used to solve the large-scale linear finite element equation system for each computation frequency. First, the coefficient matrix is... Matrix decomposition and two back-substitutions yield solutions for the two polarization modes. and .

[0042] The parallel computing approach includes two layers of parallelism. The first layer allocates a computing process for the calculation of the finite element coefficient matrix for each frequency, and the second layer uses the OpenMP threads inside the PARDIASO solver for parallelism.

[0043] S8: Solve for the electric and magnetic field components at the measuring point, where the measuring point is at the [missing information]. e When the electric field at the measuring point is within a tetrahedral element... It can be represented by vector shape functions as: Formula (9):

[0044] According to Faraday's law, the magnetic field component at the measuring point... It can be expressed by the electric field and vector shape function as follows: Formula (10):

[0045] S9: The response data for the electromagnetic method using airborne natural field sources includes impedance tensor, apparent resistivity, phase, and tilt vector. The calculation expressions for each component of the impedance tensor in the electromagnetic method using airborne natural field sources are as follows: Formula (11): , , , ; in, , , , This refers to the horizontal component of the observed magnetic field from aerial observations. , , , For the observation of the horizontal component of the electric field from a ground-based fixed base station, the subscripts “x” and “y” represent two horizontal directions, and the subscripts “1” and “2” represent two polarization modes.

[0046] After calculating the impedance tensor, the apparent resistivity and phase can be obtained, i.e.: Formula (12): , Formula (13): , The expression for calculating the tilt vector is: Formula (14): , , When only the vertical component of the magnetic field is observed from the air... , This represents the vertical component of the airborne magnetic field. , , , For ground-based fixed base stations, the horizontal component of the magnetic field is observed; when observing the three components of the magnetic field from the air... , , , , , All are magnetic field components observed from the air. Subscripts "x", "y", and "z" represent the directions of the three components, and subscripts "1" and "2" represent the two polarization modes.

[0047] Example 1 discloses a layered model, such as Figure 2As shown. The overall model range is 5km × 5km × 6km, with an air layer of 1km × 5km × 5km, and a uniform underground half-space of 5km × 5km × 5km. The first layer is 200m thick, and the second layer is 300m thick. The ground base station coordinates are (0,0,0)m, and the aerial measuring point coordinates are (0,0,-100)m. The resistivity of the air layer in the model is set to 10. 8 Oh mm, relative permeability is 1; the resistivity of the first layer principal axis is (10 4 , 10 3 , 10 4 The first layer has a resistivity of (20, 10, 20) ohms and a relative permeability of 1; the second layer has a resistivity of (20, 10, 20) ohms and a relative permeability of 2; the third layer has a resistivity of 100 ohms and a relative permeability of 1. The calculation frequency is 10⁻⁶. 1 ~10 3 Eleven frequency points with equal logarithmic intervals were selected.

[0048] The mesh size at the model measurement points and base station was set to 1m, the tetrahedral mesh size growth rate was 1.4, the number of tetrahedral elements in the model was 99326, and the number of edges was 119118. Using serial frequency calculation and one OpenMP thread, the time to solve the finite element equation for each frequency was 7.4 seconds.

[0049] Figure 3 , Figure 4 The results of the forward modeling simulation of the electromagnetic method of airborne natural sources in this invention are consistent with the one-dimensional analytical solution, proving the correctness of the three-dimensional forward modeling simulation method of the electromagnetic method of airborne natural sources developed in this invention. It can be used to simulate the effects of anisotropic conductivity and permeability. Example 2

[0050] like Figure 1 As shown, the parallel three-dimensional forward modeling method for complex media using the electromagnetic method of airborne natural sources includes the following steps: S1: A three-dimensional geological geometric model of the study area was constructed based on the implicit modeling method of CAD. The overall range of the model is 50km×50km×50km, the core area is 10km×10km×10km, the air layer thickness is about 25km, the depth of the uniform underground half space is about 25km, and two anomalous blocks are embedded, each with a size of 2km×2km×1km, and the coordinates of the center of the anomalous blocks are (2.5,0,0.5)km and (-2.5,0,0.5)km, respectively; S2: The geometric model is meshed using an unstructured tetrahedral mesh. Local mesh refinement is applied at ground base stations and aerial measurement points, while a coarse mesh is used outside the core area, thus generating an unstructured tetrahedral mesh for the entire geological model of the study area. The mesh size at model measurement points and base stations is set to 50m, the maximum mesh size in the core area is set to 1000m, the tetrahedral mesh growth rate is set to 1.4, the number of tetrahedral elements in the model is 336,831, and the number of edges is 397,325. S3: Assign arbitrary anisotropic conductivity tensors and magnetic permeability to the tetrahedral mesh elements of the model. Set the air layer resistivity of the model to 10. 8 The background resistivity is 100 ohmm, the resistivity of the two anomalies is 10 ohmm and 1000 ohmm respectively, and the relative permeability of the entire model is 1.

[0051] S4: After eliminating the magnetic field term from Maxwell's equations in the frequency domain, we can obtain the following vector electric field control equation (1): Formula (1): , Where, i= ω represents angular frequency, and the calculated frequency is 10. 1 ~10 3 Take 21 frequency points with equal logarithmic intervals. For electric field, Permeability, For conductivity tensor, , , , , , , , , and These are the conductivity components in each direction; S5: Using Galerkin finite element analysis, the frequency domain electric field control equation (1) is processed to obtain the weak integral form of the differential equation: Formula (2): , Where Ω represents the simulated region, For vector interpolation basis functions, Let represent the electric field to be solved within the study area.

[0052] Vector interpolation basis function Represented using first-order linear basis functions: Where L is the node basis function, l is the length of the edge, and the subscript " "and" "Represents the start and end points of the edge."

[0053] The electric field within a tetrahedral element is discretized using first-order linear vector basis functions. It can be determined by the electric field at the center of each edge. definition: Formula (3): , Where e represents the e-th tetrahedral element, and j represents the local number of the edge inside the tetrahedron.

[0054] Using the first vector Green's theorem, the discrete form of equation (2) is obtained in each tetrahedral element. A sparse set of finite element coefficient equations is assembled using the unstructured tetrahedral mesh of the model, electrical conductivity, and magnetic permeability parameters: Formula (4): , in, These are the right-hand terms of the equation after applying boundary conditions. and These are the stiffness matrix and the mass matrix, respectively, and their specific discrete forms are as follows: Formula (5): , Formula (6): , in, and Let these represent the stiffness matrix and the mass matrix, respectively. For vector basis functions, For the cell region to be solved, It is a volume infinitesimal element. and Let represent the permeability and resistivity tensors of the element, respectively. The superscript "e" indicates the e-th tetrahedral element, and the subscripts "j" and "k" indicate the edge number of the element. The integral of the above elements is calculated using analytical methods.

[0055] S6: Set Dirichlet boundary conditions. Apply an excitation source at the top of the model. and The electric fields at the lateral and bottom boundaries of the model can be obtained through one-dimensional forward modeling. Let the center coordinates of the edges at the boundary of the three-dimensional model be... The electric field at the center of the edge can be obtained through a one-dimensional forward modeling algorithm. Project the one-dimensional analytical solution onto the center of the edge of the three-dimensional model boundary: Formula (7): , in, It is the tangential unit vector of the boundary edge of the 3D model. It is the electric field at the center of the edge. It is the electric field component obtained by one-dimensional forward modeling.

[0056] By applying boundary conditions using the "01" assignment method, a sparse system of linear finite element equations to be solved is obtained: Formula (8): , in, and The finite element coefficient matrix A is the same for both polarization modes, but the source direction b on the right-hand side of the equation is different.

[0057] S7: Parallel computation is used for different computation frequencies. The direct solver PARDIASO is used to solve the large-scale linear finite element equation system for each computation frequency. First, the coefficient matrix is... Matrix decomposition and two back-substitutions yield solutions for the two polarization modes. and .

[0058] The parallel computing approach includes two layers of parallelism. The first layer allocates a computation process for each frequency's finite element coefficient matrix calculation. The second layer utilizes OpenMP threads within the PARDIASO solver for parallel processing. This embodiment sets up 21 parallel processes, with each process allocated two OpenMP parallel threads.

[0059] S8: Solve for the electric and magnetic field components at the measuring point, where the measuring point is at the [missing information]. e When the electric field at the measuring point is within a tetrahedral element... It can be represented by vector shape functions as: Formula (9):

[0060] According to Faraday's law, the magnetic field component at the measuring point... It can be expressed by the electric field and vector shape function as follows: Formula (10):

[0061] S9: The response data for the electromagnetic method using airborne natural field sources includes impedance tensor, apparent resistivity, phase, and tilt vector. The calculation expressions for each component of the impedance tensor in the electromagnetic method using airborne natural field sources are as follows: Formula (11): , , , , in, , , , This refers to the horizontal component of the observed magnetic field from aerial observations. , , , For the observation of the horizontal component of the electric field from a ground-based fixed base station, the subscripts “x” and “y” represent two horizontal directions, and the subscripts “1” and “2” represent two polarization modes.

[0062] After calculating the impedance tensor, the apparent resistivity and phase can be obtained, i.e.: Formula (12): , Formula (13):

[0063] The expression for calculating the tilt vector is: Formula (14): , , When only the vertical component of the magnetic field is observed from the air... , This represents the vertical component of the airborne magnetic field. , , , For ground-based fixed base stations, the horizontal component of the magnetic field is observed; when observing the three components of the magnetic field from the air... , , , , , All are magnetic field components observed from the air. Subscripts "x", "y", and "z" represent the directions of the three components, and subscripts "1" and "2" represent the two polarization modes.

[0064] Example 2 discloses a geometric structure of a model with terrain block anomalies, as follows: Figure 5 As shown. The overall model range is 50km × 50km × 50km, the core area is 10km × 10km × 10km, the air layer thickness is approximately 25km, the uniform underground half-space depth is approximately 25km, and two anomalous blocks are embedded, each measuring 2km × 2km × 1km, with the center coordinates of the anomalous blocks being (2.5, 0, 0.5)km and (-2.5, 0, 0.5)km, respectively. The resistivity of the air layer in the model is set to 10. 8 The background resistivity is 100 ohms, the resistivity of the two anomalies is 10 ohms and 1000 ohms respectively, and the relative permeability of the entire model is 1. The calculation frequency is 10 ohms. 1 ~10 3 Take 21 frequency points with equal logarithmic intervals.

[0065] See Figure 6The grid size at the model measurement points and base stations is set to 50m, the maximum grid size in the core area is set to 1000m, the tetrahedral grid growth rate is set to 1.4, the number of tetrahedral elements in the model is 336831, and the number of edges is 397325.

[0066] Figures 6-10 The apparent resistivity and phase results obtained from the three-dimensional forward modeling of the airborne natural source electromagnetic method using the vector finite element method in this embodiment demonstrate that the three-dimensional forward modeling method of the airborne natural source electromagnetic method developed in this invention has the ability to simulate complex geological structures.

[0067] In Example 2, using serial frequency calculation and one OpenMP thread, the solution time for each frequency equation is approximately 64.5 seconds, and the total time required to calculate all frequencies is 1360.3 seconds. Using the parallel computing strategy of this invention, with 21 parallel processes and two OpenMP parallel threads allocated to each process, the solution time for each frequency equation is approximately 38.6 seconds, and the total time required to calculate all frequencies is 52.6 seconds. The parallel computing speed is approximately 26 times that of serial calculation, and the parallel acceleration increases further with the number of frequencies being calculated. This demonstrates that the airborne natural source electromagnetic three-dimensional forward modeling simulation method used in this embodiment has high parallel efficiency.

[0068] While the specific embodiments of the present invention have been described above in conjunction with the accompanying drawings, this is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art without creative effort based on the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A parallel three-dimensional forward modeling method for complex media using airborne natural source electromagnetic method, characterized by: The following steps are involved: S1: Based on the CAD implicit modeling method, a three-dimensional geological geometric model of the study area is constructed using terrain, stratum and geological body data; S2: Use unstructured tetrahedral meshes to divide the geometric model, locally encrypt it at the ground base station and the aerial measurement points, and use coarse meshes outside the core area to generate an unstructured tetrahedral mesh of the geological model of the entire study area; S3: Assign arbitrary anisotropic conductivity tensors and magnetic permeabilities to the unstructured tetrahedral mesh of the model based on known geological, borehole, and geophysical information; S4: After eliminating the magnetic field term from the frequency domain Maxwell equations, we get the vector electric field governing equations; S5: The vector electric field governing equations are processed using Galerkin finite element analysis to obtain the weak form integral equation of the differential equation; S6: Set Dirichlet boundary conditions and apply excitation source to the top of the model and , the electric fields at the lateral and bottom boundaries of the model are obtained by one-dimensional forward modeling; Assume that the center coordinates of the edge at the boundary of the 3D model are , the electric field at the edge center is obtained by one-dimensional forward algorithm , project the one-dimensional analytical solution to the edge center of the three-dimensional model boundary; S7: Parallel computing is used for different calculation frequencies. The direct solver PARDIASO is used to solve the large linear equations of finite element at each calculation frequency. The coefficient matrix Perform matrix decomposition and perform two back substitutions to obtain the solutions under two polarization modes and ; S8: Solve the electric field component and magnetic field component of the measuring point; S9: Solve the response data of the airborne natural field source electromagnetic method, including impedance tensor, apparent resistivity, phase, and dipole vector.

2. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 1 is characterized in that: The vector electric field control equation of step S4 is: (1), Where, i= , ω represents the angular frequency, is the electric field, is the magnetic permeability, is the conductivity tensor, , , , , , , , , and are the conductivity components in each direction respectively.

3. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 2 is characterized by: The weak form integral equation of the differential equation in step S5 is: (2), where Ω represents the simulated area, is the vector interpolation basis function, is the electric field to be solved in the study area.

4. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 3 is characterized by: The vector interpolation basis function It is expressed using a first-order linear basis function: Where L is the node basis function, l is the length of the edge, and the subscript " "and" ” represents the starting point and end point of the edge.

5. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 4 is characterized in that: The electric field in the unstructured tetrahedron is discretized using a first-order linear vector basis function. The electric field at the center of each edge definition: ,in represents the e-th tetrahedral unit, j represents the local number of the edge within the tetrahedron; The discretized form of equation (2) is obtained in each tetrahedral element using the first vector Green's theorem, and the model unstructured tetrahedral mesh, conductivity and permeability parameters are assembled into a sparse finite element coefficient equation set: （4）， in, is the right side of the equation after the boundary conditions are applied, and are stiffness matrix and mass matrix respectively, and the specific discretization form is: （5） （6）， in, and denote the stiffness matrix and mass matrix respectively, is the vector interpolation basis function, is the unit area to be solved, is the volume element, and They represent the magnetic permeability and resistivity tensors of the unit respectively. The superscript "e" represents the e-th tetrahedral unit. The subscripts "j" and "k" represent the edge numbers of the unit. The analytical method is used to calculate the above unit integrals.

6. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 1 is characterized by: The edge center formula of the three-dimensional model boundary in step S6 is: (7), in, is the tangent unit vector of the boundary edge of the 3D model, is the electric field at the center of the edge, is the electric field component obtained by one-dimensional forward modeling; The "01" assignment method is used to impose boundary conditions, and the sparse finite element linear equations to be solved are obtained: （8）， in, and Finite element coefficient matrices for two polarization modes The same, but the source direction on the right side of the equation different.

7. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 1 is characterized by: The parallel computing method described in step S7 includes two levels of parallelism. The first level is to allocate a computing process for the finite element coefficient matrix calculation of each frequency, and the second level is to use the internal OpenMP thread parallelism of the PARDIASO solver.

8. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 1 is characterized by: In step S8, the measuring point is e When there are tetrahedral units, the electric field at the measuring point It is represented by a vector shape function: （9）， According to Faraday's law, the magnetic field component at the measuring point It is expressed by the electric field and vector shape function: （10）。 9. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 1 is characterized by: In step S9, the calculation expressions of the impedance tensor components of the aerial natural field source electromagnetic method are as follows: ， ， ， （11）， ， in, , , , is the horizontal component of the observed magnetic field observed in the air, , , , The horizontal component of the electric field is observed by a ground fixed base station. The subscripts "x" and "y" represent the two horizontal directions, and the subscripts "1" and "2" represent the two polarization modes.

10. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 9, characterized in that: In step S9, after calculating the impedance tensor, the apparent resistivity and phase can be obtained, namely: （12） （13） The calculation expression of the tilt vector is: , （14）， When only the vertical component of the magnetic field is observed in space, , is the vertical component of the magnetic field in the air, , , , The horizontal component of the magnetic field is observed by a fixed base station on the ground; when the three components of the magnetic field are observed in the air, , , , , , All are magnetic field components observed in the air; the subscripts "x", "y", and "z" represent the directions of the three components, and the subscripts "1" and "2" represent the two polarization modes.

Images