CN113569447B - Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method - Google Patents
Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method Download PDFInfo
- Publication number
- CN113569447B CN113569447B CN202110761321.3A CN202110761321A CN113569447B CN 113569447 B CN113569447 B CN 113569447B CN 202110761321 A CN202110761321 A CN 202110761321A CN 113569447 B CN113569447 B CN 113569447B
- Authority
- CN
- China
- Prior art keywords
- matrix
- time domain
- equation
- transient electromagnetic
- area
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 115
- 230000001052 transient effect Effects 0.000 title claims abstract description 52
- 239000011159 matrix material Substances 0.000 claims abstract description 88
- 238000004088 simulation Methods 0.000 claims abstract description 34
- 230000002159 abnormal effect Effects 0.000 claims abstract description 28
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 19
- 230000008569 process Effects 0.000 claims abstract description 17
- 230000005684 electric field Effects 0.000 claims description 25
- 230000005672 electromagnetic field Effects 0.000 claims description 22
- 238000004364 calculation method Methods 0.000 claims description 13
- 230000000295 complement effect Effects 0.000 claims description 12
- 239000006185 dispersion Substances 0.000 claims description 11
- 238000005070 sampling Methods 0.000 claims description 8
- 230000035699 permeability Effects 0.000 claims description 7
- 230000004044 response Effects 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 4
- 238000005259 measurement Methods 0.000 claims description 4
- 230000006698 induction Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000008030 elimination Effects 0.000 claims 1
- 238000003379 elimination reaction Methods 0.000 claims 1
- 230000008901 benefit Effects 0.000 description 3
- 230000000903 blocking effect Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 238000006467 substitution reaction Methods 0.000 description 3
- 238000011835 investigation Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- APVPOHHVBBYQAV-UHFFFAOYSA-N n-(4-aminophenyl)sulfonyloctadecanamide Chemical compound CCCCCCCCCCCCCCCCCC(=O)NS(=O)(=O)C1=CC=C(N)C=C1 APVPOHHVBBYQAV-UHFFFAOYSA-N 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 239000013535 sea water Substances 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computing Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a transient electromagnetic three-dimensional fast forward modeling method based on a Shull compensation method, which is used for successfully introducing the Shull compensation method into transient electromagnetic three-dimensional forward modeling of a time domain finite volume method, and dividing a whole simulation area into a background area and an abnormal area, wherein the conductivity of the background area is not changed in the multi-forward modeling process, and only the conductivity of the abnormal area is updated. According to the regional division position, the forward equation can be decomposed into different parts by adopting the Shulting decomposition, and after the conductivity of the abnormal region is updated, the forward model updating only needs to solve the Shulting equation corresponding to the abnormal region, so that the efficiency of transient electromagnetic three-dimensional forward modeling is greatly improved. The invention is especially suitable for the situation that the conductivity distribution of the underground medium is well known and the unknown conductivity area only occupies a small part of the simulation area, and can reduce the scale of the coefficient matrix equation when the conductivity of the small area changes in the forward modeling process, thereby improving the simulation efficiency of the multi-model transient electromagnetic forward modeling.
Description
Technical Field
The invention relates to the technical field of geographic information, in particular to a transient electromagnetic three-dimensional fast forward method based on a Shuerbu method.
Background
Transient electromagnetic methods (Transient Electromagnetic Method, TEM), which may also be referred to as Time-domain electromagnetic methods (Time-Domain Electromagnetic Method, TDEM), are an important branch of geophysical electromagnetic exploration, which employ a grounded or return line source to emit a pulsed electromagnetic field into the subsurface and infer the distribution of subsurface conductivity by observing the pure secondary field induced by the subsurface medium. The method is widely applied to deep structure research of the earth, mineral resource exploration, hydrogeological investigation, environmental engineering investigation and the like due to the characteristics of high detection efficiency, low exploration cost, better resolution, flexible and changeable devices and the like.
Currently, transient electromagnetic three-dimensional forward modeling is mainly based on two strategies: 1) The frequency spectrum method, namely a frequency domain-to-time domain algorithm, firstly determines a frequency point distribution range according to an analog time range and a model subdivision condition, calculates frequency domain response of each frequency point, and then transfers to a time domain through inverse Fourier transform; 2) The direct time domain method is characterized in that the method starts from a time domain Maxwell equation set directly, the time domain and the space domain are discrete, and the transient electromagnetic field at each moment is solved step by step in an implicit or display iteration mode. The numerical methods can be classified into finite volume methods, finite difference methods, finite volume methods, finite element methods, and the like, and each of these numerical methods has its characteristics.
In transient electromagnetic three-dimensional forward modeling, no matter how numerical method is adopted, the dimension of each solved coefficient matrix equation corresponds to the overall subdivision of the model, and for a more complex three-dimensional problem, the scale of the coefficient matrix equation is larger, and the calculation time is longer. The transient electromagnetic inversion or time-shift transient electromagnetic forward modeling and other processes involve multiple forward modeling of models, and only the conductivity of a limited area between the models is updated and changed, but the conductivity of most areas is not changed. Limited by the original forward modeling method, the forward calculation of the updated model still requires the solution of the coefficient matrix equation corresponding to the whole segmentation, so that the calculation efficiency is low, and the actual application of transient electromagnetism is affected.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a transient electromagnetic three-dimensional fast forward modeling method based on a Shuerwho's complement method, which is introduced into multi-model transient electromagnetic three-dimensional simulation, aiming at reducing the scale of coefficient matrix equation when only small area conductivity changes in the forward modeling process, thereby improving the simulation efficiency of multi-model transient electromagnetic forward modeling.
The technical scheme adopted for solving the technical problems is as follows:
the invention provides a transient electromagnetic three-dimensional fast forward modeling method based on a Shultier method, which comprises the following steps:
step 1, reading in a model file and a data file of transient electromagnetic three-dimensional forward modeling in a simulation area;
step 2, obtaining a time domain Maxwell equation set of the model, and performing time domain implicit discrete on electromagnetic and magnetic fields of the time domain Maxwell equation set by adopting a backward Euler method;
step 3, performing space domain dispersion on the Maxwell equation set after the time domain implicit dispersion by adopting a simulated finite volume method based on a structured Yee grid to obtain a discrete form of a Helmholtz equation;
step 4, analyzing a Maxwell equation set with discrete time domain and space domain by a simulated finite volume method to obtain an integral coefficient matrix A at each iteration moment, and constructing a simulation equation of the whole simulation time region;
step 5, dividing the whole simulation area into a background area b and an abnormal area a according to the actual geological conditions and the existing information of the simulation area, classifying the subdivision grid, and dividing the whole coefficient matrix into a block matrixForm including abnormal region coefficient matrix A aa Background area coefficient matrix A bb And a cross-region coefficient matrix A ba 、A ab ；
Step 6, carrying out Shu 'S complement decomposition on the integral coefficient matrix of the primary model, storing the Shu' S complement matrix S, and solving SE at each iteration moment a ＝Y a And U bb E b ＝Y b -U ba E a Equation, obtaining electric field distribution E at sampling point, and storing matrix decomposition factor obtained by decomposing each iteration moment for forward modeling of subsequent updated modelSimulating;
step 7, obtaining transient electromagnetic responses of the model corresponding to each time channel at each measuring point through time domain interpolation and space domain interpolation in sequence;
step 8, performing simulated finite volume analysis on the updated model to obtain a coefficient matrix of the whole model subdivision region, and performing simulated finite volume analysis on a coefficient matrix A of the abnormal region aa Update is performed and S' =s-se:Sup>A is used aa +A′ aa Updating the Shu's complement matrix, while the coefficient matrix of the background area and the crossing area is kept unchanged;
step 9, adopting an equation direct solution to update a model to obtain a ShuerBu matrix equation SE a ＝Y a Solving to obtain electric field component E corresponding to the abnormal region a Then to equation U bb E b ＝Y b -U ba E a Solving to obtain electric field E of background area b ；
Step 10, returning to the step 7 to perform interpolation to obtain transient electromagnetic responses of the update model corresponding to each time channel at each measuring point;
and step 11, after the model is updated again, repeating the processes of step 8, step 9 and step 10 until the model is stopped being updated.
Further, in the step 1 of the present invention:
the model file and the data file comprise: model parameters, transceiver parameters, simulation time ranges, and related time step iteration parameters.
Further, in the step 2 of the present invention:
the time domain Maxwell equation set is expressed as:
wherein E and B are respectively the electric field intensity and the magnetic induction intensity, and E is expressed as V.m -1 B is T; mu and epsilon are respectively the permeability and the dielectric constant in vacuum, and mu unit is H.m -1 Sigma is the medium conductivity, and the unit is S/m; j (J) r (t) represents the current density of the external current source, and the unit is A.m -2 ；
Performing time domain dispersion on the time domain Maxwell equation set by adopting a backward Euler method, and considering the time step as δt, the equation at each discrete time point is expressed as follows:
wherein the superscript of each field value represents the iteration moment in the time domain iteration process.
Further, in the step 3 of the present invention:
performing space domain dispersion on the time domain Maxwell equation set by adopting a finite volume method to perform forward calculation;
by adopting a Yee grid structured subdivision mode, an electric field is defined at an edge, a magnetic field is defined at a face center, and a volume range definition mode is controlled by combining an electromagnetic field component, a limited volume discrete form corresponding to a limited volume discrete form can be obtained after limited volume discrete and analysis are carried out on a Maxwell equation set after time domain implicit discrete:
wherein ,Mfμ 、M eσ and Meε Respectively corresponding to the average matrix of the permeability, conductivity and permittivity of the medium,is a discrete source item; eliminating the magnetic field B at time i+1 i+1 To obtain the electric field E i+1 Discrete form of the time domain Helmholtz equation:
wherein, the electric field E is defined at the edge of the grid, and the magnetic field B is defined at the center of the grid cell.
Further, in the step 4 of the present invention:
simplifying the obtained time domain Helmholtz equation into:
A i+1 E i+1 ＝R i+1
wherein ,Ai+1 The coefficient matrix at the iteration moment of i+1 is a large sparse matrix, R i+1 Is the corresponding right-end item;
because the coefficient matrix A of each iteration moment is only relevant to the step length of the current moment and does not depend on the field value distribution, an equation capable of repeatedly utilizing the decomposition factors is adopted to directly solve the problem; meanwhile, for the comprehensive simulation efficiency and simulation precision, a self-adaptive step doubling algorithm is adopted to define the step; the method comprises the steps of firstly determining an initial time step delta t according to the earliest simulation time point, designating the number m of equal step sizes, attempting to double the step size to 2 delta t every m step sizes, judging whether the step size is doubled according to the precision requirement, and if not, continuing to iterate according to the original step size, and waiting for the next m step sizes.
Further, in the step 5 of the present invention:
dividing the whole model area into two parts, wherein one part is a background area, the conductivity of the area is a known condition, and the area is defined as b; secondly, an abnormal region, the conductivity of which is unknown and updated with model update, which is defined as a, and the time domain Helmholtz equation is further simplified into:
AE＝R
according to the subdivision condition of the areas a and b, the method is expressed as the following block matrix form:
wherein ,Abb And A is a aa The two block matrixes are respectively related to conductivity parameters of the background area b and the abnormal area a only and are influenced by grid subdivision; block matrix A ba And A is a ab The mutual influence between the two subdivision areas is represented, and is only controlled by mesh subdivision and area division, and is irrelevant to the medium conductivity parameter.
Further, in the step 6 of the present invention:
the method is characterized in that an equation direct solution is adopted for solving, LU decomposition is firstly needed to be carried out on a matrix A, and according to the situation of a partitioned matrix, the LU decomposition is expressed as follows:
wherein ,Lbb 、L ab and Laa in the form of a blocking matrix of a lower triangular matrix L, U bb 、U ba and Uaa In the form of a blocking matrix of an upper triangular matrix U; the following linear relationship between the matrices is obtained:
L bb U bb ＝A bb
L ab U bb ＝A ab
L bb U ba ＝A ba
and defines the following matrix:
the linear transformation is performed on the obtained matrix called as a Shu' er complement matrix:
the intermediate vector Y is introduced in the above decomposition process:
is converted into the following form:
and (3) disassembling to obtain:
SE a ＝Y a
U bb E b ＝Y b -U ba E a
solving SE a ＝Y a Obtaining an electric field component E corresponding to the abnormal region a a The method comprises the steps of carrying out a first treatment on the surface of the Solving U again bb E b ＝Y b -U ba E a Obtaining an electric field component E corresponding to the background region b b And the equation solving is completed.
Further, in the step 7 of the present invention:
firstly, carrying out time domain interpolation, carrying out logarithmic attenuation on a transient electromagnetic field along with the logarithmic time, carrying out time domain interpolation on the transient electromagnetic field by adopting double-logarithmic linear interpolation, wherein a time domain interpolation formula is as follows:
wherein ,tobs 、t i and ti+1 Respectively representing an observation time channel and two iteration moments adjacent to the observation actual measurement time, wherein f () represents a corresponding interpolation field value;
the spatial domain interpolation formula is:
wherein ,Bobj Interpolation of (a) first requires finding the phase with itAnd 8 adjacent z-direction magnetic field sampling points form 8 hexahedrons, and interpolation calculation is carried out according to the volume proportion of the 8 hexahedrons.
The invention has the beneficial effects that: the transient electromagnetic three-dimensional fast forward modeling method based on the Shulting method is particularly suitable for the situation that the conductivity distribution of an underground medium is well known and an unknown conductivity area only occupies a small part of a simulation area, such as time-shifting transient electromagnetic simulation, transient electromagnetic three-dimensional inversion and the like; the invention introduces a ShuerBuddha decomposition method into multi-model transient electromagnetic three-dimensional simulation, and aims to reduce the scale of a coefficient matrix equation when only small area conductivity changes in the forward modeling process, thereby improving the simulation efficiency of the multi-model transient electromagnetic forward modeling.
Drawings
The invention will be further described with reference to the accompanying drawings and examples, in which:
FIG. 1 is a schematic view of a spatially staggered grid of (a) an embodiment of the present invention, with an electric field E defined at the center of the edge and a magnetic field B defined at the center of the face. The conductivity tensor and the magnetic permeability are defined in the unit body center; (b) a control volume schematic of an x-direction electric field component; (c) a control volume schematic of the x-direction magnetic field component;
FIG. 2 is a flow chart of a time-step strategy and equation solving according to an embodiment of the present invention;
FIG. 3 is a vertical magnetic field B according to an embodiment of the present invention z A spatial domain interpolation schematic;
FIG. 4 is a schematic diagram of background and anomaly regions of an embodiment of the present invention; (a) The conductivity of the seawater or air is a known parameter, and the conductivity of the underground medium is an unknown parameter; (b) The conductivity information of the relevant surrounding rock is also known information besides the ocean or air conductivity is known;
FIG. 5 is a three-dimensional finite volume Ye interlaced mesh subdivision 2D schematic of an embodiment of the present invention; wherein the blue region is an abnormal region, and the white region is a background region;
FIG. 6 is a time-consuming comparison of a conventional direct solution and a Shuerbu solution according to an embodiment of the present invention;
fig. 7 is a flow chart of a method of an embodiment of the present invention.
Detailed Description
The present invention will be described in further detail with reference to the drawings and examples, in order to make the objects, technical solutions and advantages of the present invention more apparent. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the scope of the invention.
The following explanation of terms is referred to in the embodiments of the present invention:
transient electromagnetic method: the transient electromagnetic method can also be called a time domain electromagnetic method, the English of which is Transient Electromagnetic Method and is simply called TEM, is an important branch of the geophysical electromagnetic method, and generally has two forms of loop source emission and ground source emission. The method has the characteristics of high exploration efficiency, low cost, higher precision and the like;
decomposition of Shuerbu: english is Schur complete, the method is introduced in 1968 by Haynsworth, the method is an algebraic regional decomposition method, and the reduced order processing of solving a matrix equation set can be realized;
forward modeling: forward computing from model parameters to response data;
wherein the following capital bold letters represent vectors or matrices.
1. Time domain electromagnetic field control equation
The three-dimensional forward modeling of the transient electromagnetic method related by the invention needs to start from a time domain Maxwell equation set, and for a given space calculation region and time domain modeling range, the time domain Maxwell equation set is expressed as follows
Wherein E (V.m) -1 ) And B (T) is the electric field intensity and the magnetic induction intensity respectively; mu (H.m) -1 ) And epsilon are the permeability and dielectric in vacuum, respectivelyConstant, σ (S/m) is the medium conductivity; j (J) r (t)(A·m -2 ) Representing the current density of the applied current source.
Using the back-off euler method (Backward Euler Method) to perform time domain dispersion on equations (1) and (2), considering the time step δt, equations (1) and (2) at each discrete time point can be expressed as follows:
wherein the superscript of each field value indicates the iteration instant in the time domain iteration process.
The invention uses finite volume method (Finite Volume Method, FVM) to perform space domain discrete on the time domain Maxwell equation set to perform forward calculation. The method can accurately discrete a continuous differential equation while maintaining physical properties of a vector field in a continuous form. The method can adopt not only structured grid subdivision, but also unstructured subdivision, and has advantages compared with a limited difference method. The Yee grid structured subdivision mode shown in the figure 1 (a) (an electric field is defined at an edge, a magnetic field is defined at a face center) is adopted, and the limited volume discrete mode corresponding to the limited volume discrete mode can be obtained after the limited volume discrete mode and the analysis of the formulas (3) and (4) are combined with the electromagnetic field component control volume range definition mode shown in the figures 1 (b) and (c).
wherein Mfμ 、M eσ and Meε Respectively corresponding to the average matrix of the permeability, conductivity and permittivity of the medium,is a discrete source item. Substituting the expression (5) into the expression (6) eliminates the magnetic field B at the time of i+1 i+1 Can be obtained with respect to the electric field E i+1 Discrete form of the time domain Helmholtz equation:
2. initial field solution
In the three-dimensional time domain electromagnetic field simulation, the electromagnetic field at zero moment is used as an initial condition, and the electromagnetic field at each moment is obtained through iteration after time. The loop source transient electromagnetic method usually adopts a down step current excitation mode, and at the moment, stable direct current exists in the loop, and the current can be excited in the whole space to generate a stable static magnetic field, namely an initial field of the transient electromagnetic. Combining (5) and (6) and making δt approach infinity, the static magnetic field excited by the current in the loop can be obtained:
the magnetic field is required to be combined with the magnetic field under the condition of no dispersionAnd (5) carrying out stability solving. The staggered grid sampling mode adopted by the invention meets the conversion form of a rotation operator and a gradient operator, and the non-dispersion condition can be expressed as follows:
Grad T B 0 ＝0 (9)
both ends are multiplied byAnd adding with formula (8) to obtain:
solving the above equation can obtain the initial electromagnetic field distribution information required by the present invention.
3. Time step strategy
The time domain Helmholtz equation shown in the step (7) is iteratively solved after time, which is an essential step for solving the time domain electromagnetic field, and the simulation efficiency is directly influenced by reasonably selecting an equation solving method and a step length updating mode. Equation (7) can be reduced to:
A i+1 E i+1 ＝R i+1 (11)
wherein Ai+1 The coefficient matrix at the iteration moment of i+1 is a large sparse matrix, R i+1 Is the corresponding right-hand item. Because the coefficient matrix A of each iteration moment is only related to the step length of the current moment and does not depend on the field value distribution, the invention solves the formula (11) by adopting an equation direct solving method capable of repeatedly utilizing the decomposition factors. Meanwhile, in order to integrate simulation efficiency and simulation precision, the invention adopts a self-adaptive step doubling algorithm to define the step. The method comprises the steps of firstly determining an initial time step delta t according to the earliest simulation time point, designating the number m of equal step sizes, attempting to double the step size to 2 delta t every m step sizes, judging whether the step size is doubled according to the precision requirement, and if not, continuing to iterate according to the original step size, and waiting for the next m step sizes. The overall step strategy and equation solving scheme is shown in fig. 2.
4. Time domain electromagnetic field interpolation
By adopting the mode to carry out iterative solution on the equation (7) in the whole time simulation range, the electromagnetic field at each iteration moment sampling point can be obtained, the measuring point position and the measuring time channel in the field actual measurement are not overlapped with the electromagnetic field sampling point and the iteration moment in most cases, and at the moment, the space domain and the time domain interpolation are needed to obtain the needed electromagnetic field.
Considering the high efficiency of calculation, the invention firstly carries out time domain interpolation, and generally carries out logarithmic attenuation on the transient electromagnetic field along with the logarithm of time, so the invention adopts double-logarithmic linear interpolation to carry out time domain interpolation on the transient electromagnetic field, and a specific interpolation formula is as follows:
wherein tobs 、t i and ti+1 The observation time channel and two iteration moments adjacent to the observation actual measurement time are respectively represented, and f () represents the corresponding interpolation field value.
The spatial domain interpolation of the invention adopts a mode of weighted average of diagonal volumes, and the interpolation field obtained by the method is relatively stable and effective. As shown in FIG. 3, the magnetic field B in the z-direction z For example, B obj Firstly, 8 magnetic field sampling points in the z direction adjacent to the interpolation of the (b) are needed to be found and 8 hexahedrons are formed, and interpolation calculation is carried out according to the volume proportion of the 8 hexahedrons by adopting the following formula:
5. shuer correction modeling method
According to the invention, the Shultier compensation method is introduced into the transient electromagnetic three-dimensional forward modeling, so that the calculation efficiency of the multi-model forward modeling can be improved. In geophysical exploration, physical parameters of certain positions of a research area are known in advance, and as shown in fig. 4 (a), the conductivity of seawater or air is a known parameter, and the conductivity of an underground medium is an unknown parameter; as further shown in fig. 4 (b), in addition to the ocean or air conductivity being known, the conductivity information of the relevant surrounding rock is also known. In such cases, the known information may be used to improve forward efficiency using the schulp method.
The present invention divides the entire model area into two parts, one of which is defined as a background area (the area conductivity is a known condition) and the other of which is defined as an abnormal area (the area conductivity is unknown and updated with model update) and a, as shown in fig. 5. For ease of analysis, equation (11) is further simplified to:
AE＝R (14)
according to the subdivision of the areas a and b, the expression (14) can be expressed as a block matrix as follows:
wherein Abb And A is a aa The two block matrices are respectively related to conductivity parameters of the background area b and the abnormal area a only and are influenced by grid subdivision. Block matrix A ba And A is a ab The mutual influence between the two subdivision areas is represented, and is only controlled by mesh subdivision and area division, and is irrelevant to the medium conductivity parameter. Solving the formula (14) by adopting an equation direct solving method, firstly, carrying out LU decomposition on the matrix A, and according to the situation of the partitioned matrix shown in the formula (15), the LU decomposition can be expressed as follows:
wherein Lbb 、L ab and Laa in the form of a blocking matrix of a lower triangular matrix L, U bb 、U ba and Uaa In the form of a segmented matrix of an upper triangular matrix U. By the above equation, the following linear relationship between matrices can be obtained:
and defines the following matrix:
known as the schulk-complement matrix, the substitution of (16) and (18) into (15) and the linear transformation can be obtained:
the intermediate vector Y is introduced in the above decomposition process:
the equation (19) is converted into the following form:
the disassembly of formulas (20) and (21) can result in:
SE a ＝Y a (22)
U bb E b ＝Y b -U ba E a (23)
to solve the equation (15), the electric field component E corresponding to the abnormal region a can be obtained by solving the equation (22) a The method comprises the steps of carrying out a first treatment on the surface of the Solving (23) to obtain electric field component E corresponding to background region b b And the equation solving is completed. Notably the formula (23) U bb Since the solution (23) is a linear operation, the time and the computing resource consumption are small.
The formulas (15) to (23) are an integral process of the Shultier solving, and the integral process is only needed in the solving process of the first model. When the model is updated and the conductivity of the abnormal region a changes, the conductivity of the background region b does not change. At this time update A of the model bb 、L bb 、U bb 、A ab 、A ba 、L ab 、U ba The matrix does not change and the solving process of equation (21) does not change, i.e. the intermediate vector Y remains unchanged. To obtain the electric field component E 'corresponding to the abnormal region of the updated model' a The Shu' S complement matrix S of the updated model is updated and then solved (22). The update mode of the Shu's complement matrix of the new model is as follows:
S′＝S-A aa +A′ aa (24)
thus, the present inventionThe solution to the transient electromagnetic update model coefficient matrix (14) in the clear can be expressed as the following steps: (1) Obtaining a block matrix form of a coefficient matrix A through time domain finite volume analysis; (2) updating the Shultier matrix by adopting the method (24); (3) Solving equation (22) corresponding to the updated matrix S' to obtain the electric field E of the abnormal region a The method comprises the steps of carrying out a first treatment on the surface of the (3) Solving (23) by a back-substitution mode to obtain a background area electromagnetic field E b 。
The most time-consuming part of the above process is the solution of step (3), i.e. the sulf matrix (22), and the number of orders of the sulf matrix corresponds to the number of grid edges of the abnormal region, which is generally much smaller than the number of grids of the whole simulation region. Therefore, the reduction of the coefficient matrix equation can be completed in this way, so that the calculation efficiency is improved.
6. Inventive method test
In order to test the application effect of the method, a transient electromagnetic three-dimensional model is designed, and the solving efficiency of a single system equation is tested. The total mesh number of the model subdivision is 51 multiplied by 69 multiplied by 63, and three different mesh area sizes are set for the abnormal area a: (1) 15× 15×15; (2) 12× 12×12; (3) 12X 8. Table 1 compares the time of the sal's solution, step 1 in the table represents solving equation (21) for the next generation; step 2 represents solving the schulmi-patch matrix (22) and the back-substitution solution (23). For the current model, if a common solving mode is adopted, the total time required is 95.8s. Taking a 12 x 12 anomaly area grid as an example, the sul's complement solution of the first model requires 129.6s, and the solution of the update model only needs 12.6s, so that the time is greatly reduced compared with the time of a common solution mode.
The time-consuming difference between the ordinary solving mode and the Shuerbu solving mode is more visually compared by adopting the form of the histogram shown in fig. 6. It can be seen from the figure that the time consumption of solving the coefficient matrix for the first time by the Shuer solving mode is far less than that of the common solving mode, but the time consumption of updating the model is far less than that of the common solving mode, and the smaller the abnormal area is, the smaller the time consumption is, so that the computing efficiency can be effectively improved.
Table 1 schulb solving time statistics table
It is to be understood that the invention has been described in conjunction with the appended drawings, but is not limited to the specific embodiments described above, which are intended to be illustrative only and not limiting, and that modifications and variations may be resorted to by those skilled in the art, having the benefit of this description, and all such modifications and variations are intended to be within the scope of the invention as set forth in the appended claims.
Claims (7)
1. The transient electromagnetic three-dimensional fast forward modeling method based on the Shultier compensation method is characterized by comprising the following steps of:
step 1, reading in a model file and a data file of transient electromagnetic three-dimensional forward modeling in a simulation area;
step 2, obtaining a time domain Maxwell equation set of the model, and performing time domain implicit discrete on electromagnetic and magnetic fields of the time domain Maxwell equation set by adopting a backward Euler method;
step 3, performing space domain dispersion on the Maxwell equation set after the time domain implicit dispersion by adopting a simulated finite volume method based on a structured Yee grid to obtain a discrete form of a Helmholtz equation;
step 4, analyzing the Maxwell equation set with discrete time domain and space domain by using a simulated finite volume method to obtain an overall coefficient matrix at each iteration momentConstructing a simulation equation of the whole simulation time area;
step 5, dividing the whole simulation area into background areas according to the actual geological conditions and the existing information of the simulation areabAnd abnormal regionaClassifying the split grids, and dividing the whole coefficient matrix into block matricesForm, including abnormal region coefficient matrix->Background area coefficient matrix->And a cross region coefficient matrix +.>、/>；
Step 6, carrying out Shu's complement decomposition on the integral coefficient matrix of the primary model and storing the Shu's complement matrixSolving for +.>And->Equation to obtain electric field distribution at sampling pointEThe matrix decomposition factors obtained by decomposing at each iteration moment are stored and used for forward modeling of a subsequent updated model; in the step 6:
solving by adopting equation direct solution, firstly, the matrix needs to be solvedAPerforming LU decomposition, which is expressed as the following form according to the case of the block matrix:
wherein ,、/> and />Is a lower triangular matrixLIn the form of a block matrix>、/> and />Is an upper triangular matrixUIs divided into a block matrix form; the following linear relationship between the matrices is obtained:
and defines the following matrix:
the linear transformation is performed on the obtained matrix called as a Shu' er complement matrix:
introducing intermediate vectors in the above decomposition processY：
Is converted into the following form:
and (3) disassembling to obtain:
solving forObtaining an abnormal regionaCorresponding electric field component>The method comprises the steps of carrying out a first treatment on the surface of the Solving for->Obtaining background regionsbCorresponding electric field component>The equation solving is completed;
step 7, obtaining transient electromagnetic responses of the model corresponding to each time channel at each measuring point through time domain interpolation and space domain interpolation in sequence;
step 8, performing simulated finite volume analysis on the updated model to obtain a coefficient matrix of the whole model subdivision region, and performing simulated finite volume analysis on the coefficient matrix of the abnormal regionUpdate is performed and the formula +.>Updating the Shu's complement matrix, while the coefficient matrix of the background area and the crossing area is kept unchanged;
step 9, adopting equation direct solution to update the Shuerbu matrix equation of the modelSolving to obtain electric field component corresponding to the abnormal region>Then>Solving to obtain electric field of background area>；
Step 10, returning to the step 7 to perform interpolation to obtain transient electromagnetic responses of the update model corresponding to each time channel at each measuring point;
and step 11, after the model is updated again, repeating the processes of step 8, step 9 and step 10 until the model is stopped being updated.
2. The transient electromagnetic three-dimensional fast forward method based on the schulplement method according to claim 1, wherein in the step 1:
the model file and the data file comprise: model parameters, transceiver parameters, simulation time ranges, and related time step iteration parameters.
3. The transient electromagnetic three-dimensional fast forward method based on the schulplement method according to claim 1, wherein in the step 2:
the time domain Maxwell equation set is expressed as:
wherein , and />The electric field intensity and the magnetic induction intensity are respectively +.>The unit is->，/>The unit is T; /> and />Magnetic permeability and dielectric constant in vacuum, +.>The unit is->，/>The conductivity of the medium is S/m; />Represents the current density of the applied current source in +.>；
Performing time domain dispersion on the time domain Maxwell equation set by adopting a backward Euler method, and considering the time step asThe equation at each discrete point in time is expressed as follows:
wherein the superscript of each field value represents the iteration moment in the time domain iteration process.
4. A transient electromagnetic three-dimensional fast forward method based on a schulplement method according to claim 3, wherein in said step 3:
performing space domain dispersion on the time domain Maxwell equation set by adopting a finite volume method to perform forward calculation;
by adopting a Yee grid structured subdivision mode, an electric field is defined at an edge, a magnetic field is defined at a face center, and a volume range definition mode is controlled by combining an electromagnetic field component, a limited volume discrete form corresponding to a limited volume discrete form can be obtained after limited volume discrete and analysis are carried out on a Maxwell equation set after time domain implicit discrete:
wherein ,、/> and />Average matrix corresponding to permeability, conductivity and permittivity of medium, respectively,/->Is a source after discretizationAn item; elimination ofiMagnetic field +.1 time>Obtaining +.>Discrete form of the time domain Helmholtz equation:
wherein the electric fieldEDefined at the edges of the grid, magnetic fieldsBDefined at the grid cell centroid.
5. The transient electromagnetic three-dimensional fast forward method based on the schulp method according to claim 4, wherein in the step 4:
simplifying the obtained time domain Helmholtz equation into:
wherein ,is thatiThe coefficient matrix at +1 iteration moment is a large sparse matrix,>is the corresponding right-end item;
due to the coefficient matrix at each iteration instantThe method is only related to the step length at the current moment and does not depend on field value distribution, so that an equation capable of repeatedly utilizing a decomposition factor is adopted to directly solve; meanwhile, for the comprehensive simulation efficiency and simulation precision, a self-adaptive step doubling algorithm is adopted to define the step; the method is firstly based onThe earliest simulation time point does the initial time step +.>And specifies the number of equal stepsmEvery othermThe step size tries to double the step size to 2 +.>Judging whether the current step doubling is accepted or not according to the precision requirement, if not, continuing to iterate according to the original step, and waiting for the next stepmSub-step size.
6. The transient electromagnetic three-dimensional fast forward method based on the schulp method according to claim 5, wherein in step 5:
the whole model area is divided into two parts, one of which is the background area, and the conductivity of the area is a known condition and is defined asbThe method comprises the steps of carrying out a first treatment on the surface of the The second is an abnormal region, the conductivity of which is unknown and which is updated with model update, defined asaThe time domain Helmholtz equation is further simplified to:
according to the regionaAnd (3) withbIs expressed in the form of a block matrix as follows:
wherein ,and->The two block matrixes are respectively only connected with the background areabAnd abnormal regionaIs related to the conductivity parameters of the network and is influenced by mesh subdivision; block matrix->And->The mutual influence between the two subdivision areas is represented, and is only controlled by mesh subdivision and area division, and is irrelevant to the medium conductivity parameter.
7. The transient electromagnetic three-dimensional fast forward method based on the schulplement method according to claim 1, wherein in the step 7:
firstly, carrying out time domain interpolation, carrying out logarithmic attenuation on a transient electromagnetic field along with the logarithmic time, carrying out time domain interpolation on the transient electromagnetic field by adopting double-logarithmic linear interpolation, wherein a time domain interpolation formula is as follows:
wherein ,、/> and />Respectively represent two iteration moments adjacent to the observation time channel and the observation actual measurement time,representing the corresponding interpolated field value;
the spatial domain interpolation formula is:
wherein ,interpolation of (c) first requires finding 8 adjacent theretozThe directional magnetic field sampling points form 8 hexahedrons, and interpolation calculation is carried out according to the volume proportion of the 8 hexahedrons.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110761321.3A CN113569447B (en) | 2021-07-06 | 2021-07-06 | Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110761321.3A CN113569447B (en) | 2021-07-06 | 2021-07-06 | Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113569447A CN113569447A (en) | 2021-10-29 |
| CN113569447B true CN113569447B (en) | 2023-11-03 |
Family
ID=78163780
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110761321.3A Active CN113569447B (en) | 2021-07-06 | 2021-07-06 | Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113569447B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113779818B (en) * | 2021-11-15 | 2022-02-08 | 中南大学 | Three-dimensional geologic body electromagnetic field numerical simulation method, device, equipment and medium thereof |
| CN116151154B (en) * | 2023-04-18 | 2023-07-18 | 中南大学 | Soil groundwater pollutant migration simulation method and related equipment |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107121706A (en) * | 2017-05-08 | 2017-09-01 | 厦门大学 | Aviation transient electromagnetic electrical conductivity 3-d inversion method based on Bonn iterative method |
| CN107845141A (en) * | 2017-11-27 | 2018-03-27 | 山东大学 | A kind of transient electromagnetic three-dimensional FDTD forward modeling multi-resolution meshes subdivision methods |
| WO2021068527A1 (en) * | 2019-10-12 | 2021-04-15 | 中南大学 | Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey |
-
2021
- 2021-07-06 CN CN202110761321.3A patent/CN113569447B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107121706A (en) * | 2017-05-08 | 2017-09-01 | 厦门大学 | Aviation transient electromagnetic electrical conductivity 3-d inversion method based on Bonn iterative method |
| CN107845141A (en) * | 2017-11-27 | 2018-03-27 | 山东大学 | A kind of transient electromagnetic three-dimensional FDTD forward modeling multi-resolution meshes subdivision methods |
| WO2021068527A1 (en) * | 2019-10-12 | 2021-04-15 | 中南大学 | Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113569447A (en) | 2021-10-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2022227206A1 (en) | Fully convolutional neural network-based magnetotelluric inversion method | |
| CN113569447B (en) | Transient electromagnetic three-dimensional fast forward modeling method based on Shu' er compensation method | |
| CN112287534B (en) | NUFFT-based two-dimensional magnetic anomaly fast forward modeling method and device | |
| Iglesias et al. | Adaptive regularisation for ensemble Kalman inversion | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN113553748B (en) | Three-dimensional magnetotelluric forward modeling numerical simulation method | |
| CN110852025B (en) | Three-dimensional electromagnetic slow diffusion numerical simulation method based on super-convergence interpolation approximation | |
| US11782183B2 (en) | Magnetotelluric inversion method based on fully convolutional neural network | |
| Lu et al. | 3D finite-volume time-domain modeling of geophysical electromagnetic data on unstructured grids using potentials | |
| CN104360404A (en) | Magnetotelluric regularization inversion method based on different constraint conditions | |
| Achitouv et al. | Improving reconstruction of the baryon acoustic peak: The effect of local environment | |
| CN114896564B (en) | Transient electromagnetic two-dimensional Bayesian inversion method adopting adaptive Thiessen polygon parameterization | |
| CN117538945B (en) | Three-dimensional magnetotelluric multi-resolution inversion method, device, equipment and medium | |
| CN115755199A (en) | Practical unstructured grid three-dimensional electromagnetic inversion smoothing and regularization method | |
| CN111474594B (en) | Three-dimensional time domain aviation electromagnetic fast inversion method | |
| CN110119586B (en) | Axial conductivity anisotropy transient electromagnetic three-component three-dimensional FDTD forward modeling method | |
| CN116165722A (en) | Loop source transient electromagnetic three-dimensional rapid inversion method adopting Gaussian Newton method | |
| CN113325482B (en) | Time domain electromagnetic data inversion imaging method | |
| Boyle et al. | Methods for calculating the electrode position Jacobian for impedance imaging | |
| Sarakorn | 2-D magnetotelluric modeling using finite element method incorporating unstructured quadrilateral elements | |
| CN115563791B (en) | Magnetotelluric data inversion method based on compressed sensing reconstruction | |
| Song et al. | Memory-efficient method for wideband self-adjoint sensitivity analysis | |
| CN115113286B (en) | Aviation electromagnetic data fusion three-dimensional inversion method based on multi-component frequency domain | |
| Wang et al. | A constrained scheme for high precision downward continuation of potential field data | |
| CN116522713A (en) | Aviation electromagnetic data partition three-dimensional inversion method based on partition coordinate descent algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |