Parallel 3D forward modeling method for complex media using airborne natural source electromagnetic method

Abstract

The present invention belongs to the field of geophysical exploration, and specifically relates to a parallel three-dimensional forward modeling method for complex media of airborne natural source electromagnetic method. It includes: S1: constructing a three-dimensional geological geometric model of the study area; S2: generating an unstructured tetrahedral grid of the model; S3: setting anisotropic conductivity tensor and magnetic permeability of the unstructured tetrahedral grid of the model; S4: constructing the frequency domain electric field control equation; S5: using vector finite element to discretize the electric field control equation; S6: setting boundary conditions; S7: solving finite element linear equations; S8: calculating the three components of electric field and magnetic field at the measuring point and base station; S9: calculating the impedance, apparent resistivity, phase, and dipole data at the measuring point. The present invention provides a three-dimensional forward modeling algorithm for airborne natural field source electromagnetic method, and provides a basis for the theoretical research of airborne natural field source electromagnetic method and the inversion interpretation of actual observation data, so as to enable airborne natural field source electromagnetic method to better serve the application fields such as mineral exploration and engineering exploration.

Description

Technical Field

[0001] The invention belongs to the field of geophysical exploration, and in particular relates to a three-dimensional forward modeling method of an aerial natural field source electromagnetic method suitable for use under complex geological conditions. Background Art

[0002] Ground surveys are difficult to carry out in areas with complex geological environments, such as plateaus, mountains, glaciers, and swamps. Airborne natural field electromagnetic method is an emerging geophysical exploration method that uses natural orthogonal field sources to electromagnetically detect underground electrical structures. Airborne natural field electromagnetic method usually uses helicopters or unmanned aerial vehicle platforms equipped with magnetic sensors to observe different components of the magnetic field in the air, and at the same time, a small number of base stations are set up on the ground to observe electric or magnetic field components. Due to the use of aerial data collection, this method has the advantages of high efficiency and convenience. It is very suitable for complex geomorphic environments and large-scale surveys. It has great development potential in the fields of mineral exploration and engineering surveys. The three-dimensional forward modeling method of airborne natural field electromagnetic method is the basis of three-dimensional inversion, but the current research on three-dimensional forward modeling of airborne natural field electromagnetic method has the following shortcomings: (1) Airborne natural field electromagnetic method is different from traditional ground electromagnetic observation methods. The existing ground-based magnetotelluric three-dimensional forward modeling method cannot be directly used for airborne natural source electromagnetic method simulation; (2) The application area of ​​airborne natural source electromagnetic method usually has complex terrain and geological structure, but there is currently no three-dimensional forward modeling method of airborne natural field source electromagnetic method specifically for complex geological structures; (3) Anisotropic and high permeability media are widely present in rock strata. These two phenomena are not considered in the existing three-dimensional forward modeling method of airborne natural source electromagnetic, which affects the calculation accuracy of airborne natural source electromagnetic forward simulation.

[0003] Chinese patent application CN118884546A discloses an airborne electromagnetic forward modeling method and medium based on physical information neural network, including: determining a three-dimensional model of the airborne electromagnetic exploration area; delineating multiple airborne electromagnetic measuring points based on the three-dimensional model; constructing a physical information neural network; taking the three-dimensional coordinates of different airborne electromagnetic measuring points, the secondary electric field values ​​of the boundary measuring points, and the airborne electromagnetic system parameters as input, and training and optimizing the physical information neural network with the goal of minimizing the mean square error of the loss function; inputting the three-dimensional coordinates of the target airborne electromagnetic measuring 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 measuring point; and calculating the secondary magnetic field value of the target airborne electromagnetic measuring point using Faraday's law of electromagnetic induction. The method proposed in patent CN118884546A is used for controlled source airborne frequency domain electromagnetic method, and is not suitable for airborne natural source electromagnetic method.

[0004] Chinese patent application CN114880902A discloses a time-domain vector finite element forward modeling method and device for magnetic source transient electromagnetic method, the method includes: establishing a simulated geological model; applying an excitation source to the simulated geological model after meshing, obtaining the Maxwell equations of the simulated geological model to obtain the electric field control equation; using the weighted residual method to discretize the control equation and selecting the first-class Whitney basis function to establish the finite element equation; selecting the time difference format to time discretize the finite element equation; assembling the overall stiffness matrix and calling the PADISO solver to solve the large linear sparse equation group; determining the iterative time step according to the current shut-off characteristics of the excitation source, and then iterating with a variable time step to obtain each component of the electric field, and obtaining the electric field value of each node through the basis function; according to Faraday's law of electromagnetic induction, obtaining the transient electromagnetic observation value and apparent resistivity parameter. The method disclosed in CN114880902A is suitable for magnetic source transient electromagnetic method, but not for airborne natural source electromagnetic method. Summary of the invention

[0005] The purpose of the present invention is to overcome the deficiencies of the above-mentioned prior art and to provide a parallel three-dimensional forward modeling method for complex media of the aviation natural source electromagnetic method based on unstructured tetrahedral mesh modeling technology and vector finite element method, which makes up for the current situation that there is no method to uniformly solve the difficult problems of complex geological structure, conductivity anisotropy and magnetic permeability influence in aviation natural source electromagnetic forward modeling, thereby obtaining the real aviation natural field source electromagnetic response of complex earth media.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A parallel three-dimensional forward modeling method for complex media using an airborne natural source electromagnetic method comprises the following steps:

[0008] S1: Construct a 3D geological geometric model of the study area. Based on the CAD implicit modeling method, the 3D geological geometric model of the study area is constructed using terrain, stratum and geological body data.

[0009] S2: Generate unstructured tetrahedral mesh of the model; Use unstructured tetrahedral mesh to divide the geometric model, perform local encryption at the ground base station and the aerial measurement points, and use coarse mesh to divide the areas outside the core area, so as to generate an unstructured tetrahedral mesh of the geological model of the entire study area;

[0010] S3: Anisotropic conductivity tensor and permeability settings for unstructured tetrahedral meshes of the model; arbitrary anisotropic conductivity tensors and permeabilities are assigned to tetrahedral mesh cells of the model based on known geological, borehole and geophysical information.

[0011] Construct the frequency domain electric field control equation; after eliminating the magnetic field term from the frequency domain Maxwell equations, the following vector electric field control equation formula (1) can be obtained:

[0012] Formula (1): ,

[0013] Among them, 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;

[0014] S5: Use vector finite element to discretize the electric field control equation; use Galerkin finite element analysis to process the frequency domain electric field control equation formula (1) and obtain the weak form integral equation of the differential equation:

[0015] Formula (2):

[0016] Where Ω represents the simulated area, is the vector interpolation basis function, is the electric field to be solved in the study area.

[0017] 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 can be definition:

[0018] Formula (3): ,

[0019] in represents the e-th tetrahedral unit, and j represents the local number of the edge inside the tetrahedron.

[0020] 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:

[0021] Formula (4): ,

[0022] 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:

[0023] Formula (5): ,

[0024] Formula (6): ,

[0025] in, and denote the stiffness matrix and mass matrix respectively, is the vector 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 above unit integrals are calculated by analytical method.

[0026] S6: Set boundary conditions; Set Dirichlet boundary conditions. Apply excitation source to the top of the model and , the electric field at the side and bottom boundaries of the model can be obtained by one-dimensional forward modeling. Assume that the center coordinates of the edge at the boundary of the three-dimensional model are , the electric field at the edge center can be obtained by one-dimensional forward algorithm , project the one-dimensional analytical solution to the center of the edge of the three-dimensional model boundary:

[0027] Formula (7): ,

[0028] 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.

[0029] The boundary conditions are imposed by the "01" assignment method, and the sparse finite element linear equations to be solved are obtained:

[0030] Formula (8): ,

[0031] in, and Finite element coefficient matrices for two polarization modes The same, but the source direction on the right side of the equation different.

[0032] S7: Solve finite element linear equations; parallel calculation is used for different calculation frequencies, and the direct solver PARDIASO is used to solve the large finite element linear equations for each calculation frequency. First, the coefficient matrix Perform matrix decomposition and perform two back substitutions to obtain the solutions under two polarization modes. and .

[0033] S8: Calculate the three components of the electric field and magnetic field at the measuring point and the base station; solve the electric field component and magnetic field component of the measuring point. e When there are tetrahedral units, the electric field at the measuring point It can be expressed by a vector shape function as:

[0034] Formula (9):

[0035] According to Faraday's law, the magnetic field component at the measuring point It is expressed by the electric field and vector shape function:

[0036] Formula (10): .

[0037] S9: Calculate the impedance, apparent resistivity, phase, and dipole data at the measuring point; the response data of the airborne natural field source electromagnetic method includes the impedance tensor, apparent resistivity, phase, and dipole vector. The calculation expressions of the impedance tensor components of the airborne natural field source electromagnetic method are:

[0038] Formula (11): ,

[0039] ,

[0040] ,

[0041] ,

[0042] 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 two horizontal directions, and the subscripts "1" and "2" represent two polarization modes.

[0043] After calculating the impedance tensor, the apparent resistivity and phase can be obtained, namely:

[0044] Formula (12):

[0045] Formula (13):

[0046] The calculation expression of the tilt vector is:

[0047] Formula (14): , ,

[0048] 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, , , , , , These are all 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.

[0049] Preferably, the vector interpolation basis function described in step S4 is 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.

[0050] 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 finite element coefficient matrix calculation of each frequency, and the second layer is to use the internal OpenMP thread parallelism of the PARDIASO solver.

[0051] The present invention proposes a three-dimensional forward modeling method of airborne natural field source electromagnetic method applicable to complex geological structures, which overcomes the limitation that the existing airborne natural source electromagnetic method forward modeling technology cannot comprehensively consider the influence of complex geological structures, anisotropic media and magnetic permeability. Its innovations and advantages are as follows:

[0052] 1. The present invention can accurately simulate the influence of the dual parameters of anisotropic conductivity and magnetic permeability on the aerial natural source electromagnetic data, and is applicable to the anisotropic and magnetic medium effects that are widely present in real earth media.

[0053] 2. The present invention adopts three-dimensional geological modeling technology and a vector finite element method based on an unstructured tetrahedral grid, which can accurately simulate the aerial natural source electromagnetic response of real geological structures with complex shapes.

[0054] 3. The present invention adopts parallel technology to calculate observation data of different frequencies, and has a higher calculation speed.

[0055] The present invention provides a three-dimensional forward simulation algorithm for the aerial natural field source electromagnetic method, and provides a basis for the theoretical research of the aerial natural field source electromagnetic method and the inversion interpretation of actual observation data, so that the aerial natural field source electromagnetic method can better serve the application fields such as mineral exploration and engineering exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 Flow chart of the three-dimensional forward modeling method of airborne natural source electromagnetic method;

[0057] Figure 2 Unstructured grid and electrical parameters of layered models;

[0058] Figure 3 Comparison of the forward modeling results of the present invention's airborne natural source electromagnetic method and the analytical solution for the layered model (apparent resistivity-frequency), the ground base station position is (0, 0, 0) m, and the aerial measuring point position is (0, 0, -100) m;

[0059] Figure 4 Comparison of the forward modeling results of the present invention's airborne natural source electromagnetic method and the analytical solution for the layered model (phase-frequency), the ground base station position is (0, 0, 0) m, and the aerial measurement point position is (0, 0, -100) m;

[0060] Figure 5 The geometric structure of the model with terrain block anomaly. The black dots represent the locations of the aerial magnetic field measurement points, which are 200m above the ground.

[0061] Figure 6 Unstructured grid and electrical parameters with terrain block anomaly model;

[0062] Figure 7 Forward modeling results of the xy mode apparent resistivity of the airborne natural source electromagnetic method with a terrain block anomaly model. The ground base station is located at (0, -5000, 0) m and the Y=0m survey line.

[0063] Figure 8 Forward modeling results of apparent resistivity in yx mode of airborne natural source electromagnetic method with terrain block anomaly model, ground base station location (0, -5000, 0)m, Y=0m survey line.

[0064] Fig. 9 Forward modeling results of the xy mode phase rate of the airborne natural source electromagnetic method with a terrain block anomaly model, the ground base station position is (0, -5000, 0)m, and the Y=0m survey line.

[0065] Fig.10Phase forward modeling results of the yx mode of the airborne natural source electromagnetic method with a terrain block anomaly model, the ground base station position is (0, -5000, 0)m, and the Y=0m survey line. DETAILED DESCRIPTION

[0066] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0067] The structures, proportions, sizes, etc. illustrated in the drawings of this specification are only used to match the contents disclosed in the specification for people familiar with this technology to understand and read, and are not used to limit the limiting conditions for the implementation of the present invention, so they have no substantial technical significance. Any modification of the structure, change of the proportion relationship or adjustment of the size, without affecting the effects and purposes that can be achieved by the present invention, should still fall within the scope of the technical contents disclosed by the present invention. At the same time, the terms such as "upper", "lower", "left", "right", "middle" and "one" quoted in this specification are only for the convenience of description, and are not used to limit the scope of the implementation of the present invention. The change or adjustment of their relative relationship should also be regarded as the scope of the implementation of the present invention without substantially changing the technical contents.

[0068] See also Figure 1 This embodiment discloses a flowchart of a preferred three-dimensional forward modeling method of natural source electromagnetic method of aviation. In order to better illustrate the advantages and purposes of this embodiment, an example is given: Example 1

[0069] like Figure 1 As shown in the figure, the parallel three-dimensional forward modeling method of complex media by airborne natural source electromagnetic method includes the following steps:

[0070] S1: A geometric model was constructed based on the CAD implicit modeling method. The overall scope of the model was 5km×5km×6km, the air layer was 1km×5km×5km, the underground uniform half space was 5km×5km×5km, the first layer was 200m thick, and the second layer was 300m thick.

[0071] S2: Use unstructured tetrahedral mesh to divide the geometric model, perform local encryption at the ground base station coordinate (0,0,0)m and the aerial measurement point coordinate (0,0,-100)m, set the mesh size to 1m, and use coarse mesh to divide the area outside the core area, thus generating an unstructured tetrahedral mesh of the geological model of the entire study area. The number of tetrahedral units in the model is 99326, and the number of edges is 119118;

[0072] S3: Assign arbitrary anisotropic conductivity tensors and magnetic permeabilities to the model tetrahedral mesh cells based on known geological, borehole and geophysical information. The air layer resistivity of the model is set to 10 8ohmm, relative magnetic permeability is 1; the first layer main axis resistivity is (10 4 , 10 3 , 10 4 ) ohm, relative magnetic permeability is 1; the second layer main axis resistivity is (20,10,20)ohm, relative magnetic permeability is 2; the third layer resistivity is 100 ohm, relative magnetic permeability is 1.

[0073] S4: After eliminating the magnetic field term from the frequency domain Maxwell equations, the following vector electric field control equation can be obtained (1):

[0074] Formula (1): ,

[0075] Where, i= , ω represents the angular frequency, and the calculated frequency is 10 1 ~10 3 Take 11 equally logarithmically spaced frequency points, is the electric field, is the magnetic permeability, is the conductivity tensor, , , , , , , , , and are the conductivity components in each direction respectively.

[0076] S5: Galerkin finite element analysis is used to process the frequency domain electric field control equation (1) to obtain the weak form integral equation of the differential equation:

[0077] Formula (2): ,

[0078] Where Ω represents the simulated area, is the vector interpolation basis function, is the electric field to be solved in the study area.

[0079] Vector interpolation basis functions 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.

[0080] The electric field in the tetrahedral unit is discretized using the first-order linear vector basis function. The electric field at the center of each edge can be definition:

[0081] Formula (3): ,

[0082] in represents the e-th tetrahedral unit, and j represents the local number of the edge inside the tetrahedron.

[0083] 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:

[0084] Formula (4): ,

[0085] 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:

[0086] Formula (5): ,

[0087] Formula (6): ,

[0088] in, and denote the stiffness matrix and mass matrix respectively, is a vector 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 above unit integrals are calculated by analytical method.

[0089] S6: Set Dirichlet boundary conditions. Apply excitation source to the top of the model and , the electric field at the side and bottom boundaries of the model can be obtained by one-dimensional forward modeling. Assume that the center coordinates of the edge at the boundary of the three-dimensional model are , the electric field at the edge center can be obtained by one-dimensional forward algorithm , project the one-dimensional analytical solution to the center of the edge of the three-dimensional model boundary:

[0090] Formula (7): ,

[0091] 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.

[0092] The boundary conditions are imposed by the "01" assignment method, and the sparse finite element linear equations to be solved are obtained:

[0093] Formula (8): ,

[0094] in, and Finite element coefficient matrices for two polarization modes The same, but the source direction on the right side of the equation different.

[0095] 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 .

[0096] The parallel computing method includes two levels of parallelism. The first level is to assign a computing process to the finite element coefficient matrix calculation for each frequency, and the second level is to use the internal OpenMP thread parallelism of the PARDIASO solver.

[0097] S8: Solve the electric field component and magnetic field component of the measuring point. e When there are tetrahedral units, the electric field at the measuring point It can be expressed by a vector shape function as:

[0098] Formula (9):

[0099] According to Faraday's law, the magnetic field component at the measuring point It is expressed by the electric field and vector shape function:

[0100] Formula (10):

[0101] S9: The response data of the airborne natural field source electromagnetic method include impedance tensor, apparent resistivity, phase, and dipole vector. The calculation expressions of the impedance tensor components of the airborne natural field source electromagnetic method are:

[0102] Formula (11): ,

[0103] ,

[0104] ,

[0105] ;

[0106] 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 two horizontal directions, and the subscripts "1" and "2" represent two polarization modes.

[0107] After calculating the impedance tensor, the apparent resistivity and phase can be obtained, namely:

[0108] Formula (12): ,

[0109] Formula (13): ,

[0110] The calculation expression of the tilt vector is:

[0111] Formula (14): , ,

[0112] 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, , , , , , These are all 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.

[0113] Example 1 discloses a layered model such as Figure 2 As shown. 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 thickness is 200m, and the second layer thickness is 300m. The ground base station coordinates are (0,0,0)m, and the aerial measurement point coordinates are (0,0,-100)m. The air layer resistivity of the model is set to 10 8 ohmm, relative magnetic permeability is 1; the first layer main axis resistivity is (104 , 10 3 , 10 4 ) ohmm, relative magnetic permeability is 1; the second layer main axis resistivity is (20,10,20) ohm, relative magnetic permeability is 2; the third layer resistivity is 100 ohmm, relative magnetic permeability is 1. The calculation frequency is 10 1 ~10 3 Take 11 equally logarithmically spaced frequency points.

[0114] The mesh size at the model measurement points and base stations is set to 1m, the tetrahedral mesh size growth rate is 1.4, the number of tetrahedral elements in the model is 99326, and the number of edges is 119118. When using serial frequency calculation and 1 OpenMP thread, the time to solve each frequency finite element equation is 7.4 seconds.

[0115] Figure 3 , Figure 4 It shows that the forward modeling calculation results of the aviation natural source electromagnetic method of the present invention are consistent with the one-dimensional analytical solution, which proves the correctness of the three-dimensional forward modeling method of the aviation natural source electromagnetic method developed by the present invention, and can be used to simulate the influence of anisotropic conductivity and magnetic permeability. Example 2

[0116] like Figure 1 As shown in the figure, the parallel three-dimensional forward modeling method of complex media by airborne natural source electromagnetic method includes the following steps:

[0117] S1: A three-dimensional geological geometric model of the study area is constructed based on the CAD implicit modeling method. The overall range of the model is 50km×50km×50km, the core area range is 10km×10km×10km, the air layer thickness is about 25km, the underground uniform half-space depth is about 25km, and two abnormal blocks are embedded, both of which are 2km×2km×1km in size, and the center coordinates of the abnormal bodies are (2.5, 0, 0.5)km and (-2.5, 0, 0.5)km respectively;

[0118] S2: Use unstructured tetrahedral mesh to divide the geometric model, locally encrypt it at the ground base station and the aerial measuring point, and use coarse mesh to divide the area outside the core area, so as to generate an unstructured tetrahedral mesh of the geological model of the entire study area. The mesh size at the model measuring point and base station 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 units in the model is 336831, and the number of edges is 397325;

[0119] S3: Assign arbitrary anisotropic conductivity tensors and magnetic permeabilities to the model tetrahedral mesh elements. The air layer resistivity of the model is set to 108 ohmm, the underground background resistivity is 100 ohmm, the resistivities of the two anomalies are 10 ohmm and 1000 ohmm respectively, and the relative magnetic permeability of the whole model is 1.

[0120] S4: After eliminating the magnetic field term from the frequency domain Maxwell equations, the following vector electric field control equation formula (1) can be obtained:

[0121] Formula (1): ,

[0122] Where, i= , ω represents the angular frequency, and the calculated frequency is 10 1 ~10 3 Take 21 equally logarithmically spaced frequency points. is the electric field, is the magnetic permeability, is the conductivity tensor, , , , , , , , , and are the conductivity components in each direction respectively;

[0123] S5: Galerkin finite element analysis is used to process the frequency domain electric field control equation (1) to obtain the weak form integral equation of the differential equation:

[0124] Formula (2): ,

[0125] Where Ω represents the simulated area, is the vector interpolation basis function, is the electric field to be solved in the study area.

[0126] Vector interpolation basis functions 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.

[0127] The electric field in the tetrahedral unit is discretized using the first-order linear vector basis function. The electric field at the center of each edge can be definition:

[0128] Formula (3): ,

[0129] Among them, e represents the e-th tetrahedral unit, and j represents the local number of the edge inside the tetrahedron.

[0130] 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:

[0131] Formula (4): ,

[0132] 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:

[0133] Formula (5): ,

[0134] Formula (6): ,

[0135] in, and denote the stiffness matrix and mass matrix respectively, is a vector 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 above unit integrals are calculated by analytical method.

[0136] S6: Set Dirichlet boundary conditions. Apply excitation source to the top of the model and , the electric field at the side and bottom boundaries of the model can be obtained by one-dimensional forward modeling. Assume that the center coordinates of the edge at the boundary of the three-dimensional model are , the electric field at the edge center can be obtained by one-dimensional forward algorithm , project the one-dimensional analytical solution to the center of the edge of the three-dimensional model boundary:

[0137] Formula (7): ,

[0138] 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.

[0139] The boundary conditions are imposed by the "01" assignment method, and the sparse finite element linear equations to be solved are obtained:

[0140] Formula (8): ,

[0141] in, and The finite element coefficient matrix A of the two polarization modes is the same, but the source direction b on the right side of the equation is different.

[0142] 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 .

[0143] The parallel computing method 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. In this embodiment, 21 parallel processes are set and 2 OpenMP parallel threads are allocated to each process.

[0144] S8: Solve the electric field component and magnetic field component of the measuring point. e When there are tetrahedral units, the electric field at the measuring point It can be expressed by a vector shape function as:

[0145] Formula (9):

[0146] According to Faraday's law, the magnetic field component at the measuring point It is expressed by the electric field and vector shape function:

[0147] Formula (10):

[0148] S9: The response data of the airborne natural field source electromagnetic method include impedance tensor, apparent resistivity, phase, and dipole vector. The calculation expressions of the impedance tensor components of the airborne natural field source electromagnetic method are:

[0149] Formula (11): ,

[0150] ,

[0151] ,

[0152] ,

[0153] 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 two horizontal directions, and the subscripts "1" and "2" represent two polarization modes.

[0154] After calculating the impedance tensor, the apparent resistivity and phase can be obtained, namely:

[0155] Formula (12): ,

[0156] Formula (13):

[0157] The calculation expression of the tilt vector is:

[0158] Formula (14): , ,

[0159] 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, , , , , , These are all 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.

[0160] Example 2 discloses a geometric structure of a terrain block anomaly model such as Figure 5 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 underground uniform half-space depth is about 25km, and two anomaly blocks are embedded, both of which are 2km×2km×1km in size. The center coordinates of the anomaly blocks are (2.5, 0, 0.5)km and (-2.5, 0, 0.5)km respectively. The air layer resistivity of the model is set to 10 8ohmm, the underground background resistivity is 100 ohmm, the resistivity of the two anomalies is 10 ohmm and 1000 ohmm respectively, and the relative permeability of the whole model is 1. The calculation frequency is 10 1 ~10 3 Take 21 equally logarithmically spaced frequency points.

[0161] See also Figure 6 The 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.

[0162] Figure 6-Figure 10 It shows that the apparent resistivity and phase results calculated by the three-dimensional forward modeling of the vector finite element airborne natural source electromagnetic method adopted in this embodiment prove that the three-dimensional forward modeling method of the airborne natural source electromagnetic method developed by the present invention has the ability to simulate complex geological structures.

[0163] In the model of Example 2, when serial frequency calculation and 1 OpenMP thread are used, each frequency equation takes about 64.5 seconds to solve, and the time required to calculate all frequencies is 1360.3 seconds. Using the parallel computing strategy of the present invention, 21 parallel processes are set and 2 OpenMP parallel threads are allocated to each process. Each frequency equation takes about 38.6 seconds to solve, and the time required to calculate all frequencies is 52.6 seconds. The parallel computing rate is about 26 times that of serial computing, and the parallel acceleration will be greater as the number of computing frequencies increases. It shows that the three-dimensional forward simulation method of the natural source electromagnetic method of aviation adopted in this embodiment has high parallel efficiency.

[0164] Although the above describes the specific implementation mode of the present invention in conjunction with the accompanying drawings, it 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 on the basis of the technical solution of the present invention without creative work 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; The vector electric field control equation is: (1) Where i= , ω represents the angular frequency, E is the electric field, μ is the magnetic permeability, and σ is the conductivity tensor ,σ xx , σ xy , σ xz , σ yx , σ yy , σ yz , σ zx , σ zy , and σ zz are the conductivity components in each direction respectively; 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, apply excitation sources E1=(1, 0, 0) and E2=(0, 1, 0) to the top of the model, and obtain the electric fields at the side and bottom boundaries of the model through 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; The edge center formula of the three-dimensional model boundary is: (7) in, is the tangent unit vector of the edge of the 3D model boundary, E 3d is the electric field at the center of the edge, E 1d 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: YES 1,2 =b 1,2 (8) Among them, the finite element coefficient matrix A of the two polarization modes E1=(1, 0, 0) and E2=(0, 1, 0) is the same, but the source direction b on the right side of the equation is different; S7: Parallel computing is used for different calculation frequencies, and the direct solver PARDISO is used to solve the large finite element linear equations for each calculation frequency. The coefficient matrix A is first decomposed, and the solutions E1 and E2 in two polarization modes are obtained after two back substitutions. 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 weak form integral equation of the differential equation in step S5 is: （2） Where Ω represents the simulated area, N is the vector interpolation basis function, and E is the electric field to be solved in the study area.

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 vector interpolation basis function N is represented by a first-order linear basis function: Where 𝐿 is the node basis function, l is the length of the edge, and the subscripts "i1" and "i2" represent the starting and ending points of the edge.

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 electric field in the unstructured tetrahedron is discretized using a first-order linear vector basis function. The electric field E in the tetrahedron unit e The electric field E at the center of each edge j e definition: , Where e represents the e-th tetrahedral unit, and j represents the local number of the edge inside 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: [K-iωM σ ] E=b (4) Where b is the right-hand side of the equation after the boundary conditions are applied, K and M σ are stiffness matrix and mass matrix respectively, and the specific discretization form is: （5） （6） Among them, K and M σ and represent the stiffness matrix and mass matrix respectively, N is the vector interpolation basis function, Ω is the unit area to be solved, dv is the volume element, μ and σ represent the magnetic permeability and resistivity tensors of the unit respectively, the superscript "e" represents the e-th tetrahedral unit, and the subscripts "j" and "k" represent the edge numbers of the unit. The above unit integrals are calculated using the analytical method.

5. 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 layers of parallelism. The first layer is to allocate a computing process for the finite element coefficient matrix calculation of each frequency, and the second layer is to use the internal OpenMP thread parallelism of the PARDISO solver.

6. 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）。 7. 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） ， Among them, H x1 , H x2 , H y1 , H y2 is the horizontal component of the observed magnetic field observed in the air, E x1 , E x2 , E y1 , E y2 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.

8. The complex medium parallel three-dimensional forward simulation method of the airborne natural source electromagnetic method according to claim 7 is characterized by: 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 the air, H z1 , H z2 is the vertical component of the magnetic field in the air, H x1 , H x2 , H y1 , H y2 When the ground fixed base station observes the horizontal component of the magnetic field, H x1 , H x2 , H y1 , H y2 , H z1 , H z2 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