CN111611737A - Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium - Google Patents

Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium Download PDF

Info

Publication number
CN111611737A
CN111611737A CN202010426036.1A CN202010426036A CN111611737A CN 111611737 A CN111611737 A CN 111611737A CN 202010426036 A CN202010426036 A CN 202010426036A CN 111611737 A CN111611737 A CN 111611737A
Authority
CN
China
Prior art keywords
representing
sigma
electric field
conductivity
equation
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.)
Granted
Application number
CN202010426036.1A
Other languages
Chinese (zh)
Other versions
CN111611737B (en
Inventor
李健
柳建新
郭荣文
柳卓
刘嵘
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202010426036.1A priority Critical patent/CN111611737B/en
Publication of CN111611737A publication Critical patent/CN111611737A/en
Application granted granted Critical
Publication of CN111611737B publication Critical patent/CN111611737B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides a marine controllable source electromagnetic forward modeling method for a three-dimensional arbitrary anisotropic medium, which comprises the steps of reading a forward modeling, adding a gradient-divergence operator in a double-rotation equation to form a modified control equation, then obtaining a new discretization equation set, solving the linear equation set and the like, and can quickly and accurately obtain an electric field component Ea under the seax、EayAnd EazThe steps are simplified; the gradient-divergence operator is directly added into the dual-rotation equation, the modified dual-rotation equation can be simplified into a Laplace equation, the number of non-zero elements of each row of the coefficient matrix can be reduced, and the calculation time of an iterative solver can be obviously reduced; by adopting the forward modeling method, the maximum errors of the amplitude and the phase and the reference solution are small, and the precision is high; divergence correction is not required to be applied, and a control equation meets a divergence condition; when the frequency is very low, the algorithm of the invention has the advantagesThe potential is more pronounced.

Description

Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium
Technical Field
The invention relates to the technical field of geophysical, in particular to an ocean controllable source electromagnetic forward modeling method for a three-dimensional random anisotropic medium.
Background
The marine controllable source electromagnetic method becomes a new electromagnetic exploration technology and has great effect in the field of marine oil and gas exploration.
For conventional ocean electromagnetic field numerical simulation, the isotropic geoelectric model can well reflect the underground electric structure. However, for some complex depositional environments located on the seafloor (such as in thin interbed and fractured reservoirs), the study of conductivity anisotropy has important significance for the interpretation of marine controllable source electromagnetic data. In practical marine controlled source electromagnetic exploration problems, thin reservoirs are often characterized by conductivity anisotropy, which is characterized by vertical and dip anisotropy. The ocean controlled source electromagnetic forward modeling is a process of calculating corresponding geophysical response by an analytic or numerical method under the condition of given ocean medium distribution and excitation sources.
Through the geophysical forward modeling, the distribution rule of response under different geophysical models can be researched, so that the actual exploration work is guided. However, whatever the geophysical forward modeling method, the solution problem of the linear equation system is finally regressed, and there are two main methods for solving the linear equation system in the conventional method: the first method comprises the following steps: solving the sparse linear equation set by using a direct solution method, which is relatively stable, but needs a large amount of computer memory and calculation time; and the second method comprises the following steps: the krylov subspace iteration method is more attractive in solving a large sparse linear equation set. Many scholars have made intensive studies on three-dimensional electromagnetic forward modeling by using krylov subspace iteration, and have made great progress. However, in the marine environment, due to the influence of uneven conductivity distribution and anisotropy caused by submarine topography fluctuation and sedimentation factors, the geoelectric model is often very complex, and in addition, the calculation frequency of the marine controllable source electromagnetic method is relatively low, and the divergence condition is not satisfied, so that the discretely obtained linear equation set is difficult to converge.
Therefore, the forward modeling method which can accurately simulate the ocean complex anisotropic medium model and accelerate the convergence of the linear equation set has important research significance.
Disclosure of Invention
The invention aims to provide a controllable source electromagnetic forward modeling method which is simplified in steps and can quickly and accurately obtain an electrical structure under the sea, and the specific technical scheme is as follows:
a marine controllable source electromagnetic forward modeling method for a three-dimensional arbitrary anisotropic medium comprises the following steps:
reading a forward model to obtain initial parameters;
and step two, forming a double rotation degree equation expression 2 by combining the expression 1) of the Maxwell equation set:
Figure BDA0002498700020000021
Figure BDA0002498700020000022
wherein:
Figure BDA0002498700020000023
is Hamiltonian; e is an electric field; h is a magnetic field; j. the design is a squaresIs the current density vector of the external current; ρ is the free charge; is the dielectric permittivity; omega is angular frequency; σ is the conductivity tensor; μ is the magnetic permeability in vacuum; i denotes the imaginary part, i2＝-1；
Step three, firstly, adding a gradient-divergence operator to the double-rotation equation expression 2) in the step two to form a modified control equation expression 6); discretizing the modified control equation expression 6) to obtain a discretized control equation expression 7):
Figure BDA0002498700020000024
Figure BDA0002498700020000025
wherein: ebA primary field being an electric field E; eaThe secondary field of the electric field E represents the electric field vector value to be solved after discretization; lambda [ alpha ]1And λ2A control factor (like the regularization factor in regularized inversion); c represents a discretized rotation operator which is the mapping from the unit edge to the unit area;
Figure BDA0002498700020000026
representing the transposition of a discretized rotation operator C, representing the mapping from the cell area to the cell edge, G representing a gradient operator, representing a mapping from a node to an inner edge, Λ representing a diagonal matrix formed by a scaling factor, D representing a divergence operator, representing a mapping from an inner edge to an inner node, sigma representing a vector formed by conductivity elements, sigma being diag (sigma);
Figure BDA0002498700020000027
conductivity values on the edges Λ1Representing a diagonal matrix formed by the scaling factors; sigma1Vector representing the composition of the residual conductivity element, Σ1Δ σ represents a difference between the actual conductivity model and the electrical conductivity of the earth model;
step four, obtaining expressions of matrix coefficient A and constant b according to the discretization control equation expression 7), such as expression 8) and expression 9):
Figure BDA0002498700020000028
b＝diag(-iωμΔσ)Eb+GΛ1DΣ1Eb9)；
step five, solving the linear equation set Ae ═ b to obtain and output the electric field component Eax、EayAnd Eaz。
Preferably, in the above technical solution, the parameters in the first step are conductivity values of each grid unit, a frequency list, edge lengths in the x, y and z directions, and information of all grid nodes; the grid node information is a grid node number.
Preferably, in the above technical solution, the conductivity tensor σ in the second step is expressed by expression 4):
Figure BDA0002498700020000031
wherein: sigmaxxConductivity values representing the x-axis direction; sigmayyConductivity values representing the y-axis direction; sigmazzConductivity values in the z-axis direction; sigmaxyAnd σyxConductivity values representing the combined action of the x-axis and the y-axis are specifically: sigmaxyDenotes the current density formed in the y direction by applying an electric field in the x direction, and σyxRepresents the current density formed in the x direction by the electric field applied in the y direction; sigmaxzAnd σzxConductivity values representing the x-axis and z-axis co-action are specifically: sigmaxzDenotes the current density formed in the z direction by applying an electric field in the x direction, and σzxRepresents the current density formed in the x direction by the electric field applied in the z direction; sigmayzAnd σzyConductivity values representing the combined action of the y-axis and the z-axis are: sigmayzDenotes the current density formed in the z direction by applying an electric field in the y direction, and σzyIndicating the electric field applied in the z-direction and the resulting current density in the y-direction.
Preferably, in the above technical solution, in the step three: the geoelectricity model is a 1D geoelectricity model; ebIs a primary field obtained by calculation through a 1D geoelectric model, and a secondary field E is obtained by solvinga；λ1＝1/σ，λ2＝1/Δσ。
Preferably, in the above technical solution, the step five specifically includes the following steps:
step 5.1, setting an initial value ek＝e0Calculating the residual error rk＝b-AekWherein: e.g. of the type0Representation of assignment of iterative solutionAn initial value of the electric field of (a); e.g. of the typekRepresenting the electric field value of the kth iteration; r iskRepresenting the residual value of the kth iteration;
step 5.2, carrying out iterative loop, and outputting the electric field component Ea if the residual error meets the requirementx、EayAnd Eaz(ii) a Otherwise, k is taken as k +1, and the procedure returns to step 5.1.
Preferably, in the above technical solution, the iterative loop employs a bi-conjugate gradient iterative algorithm; residual value r of given double conjugate gradient method0(ii) a If the residual error reaches the requirement 10-8Or exceeds the maximum convergence time by 3000 times, and outputs the electric field component Eax、EayAnd Eaz。
By applying the technical scheme, the forward modeling method comprises the processes of reading a forward model, adding a gradient-divergence operator in a double-rotation equation to form a modified control equation, then obtaining a unique discretization equation set, solving a linear equation set and the like, and the electric field component Ea under the ocean can be quickly and accurately obtainedx、EayAnd EazThe steps are simplified; the gradient-divergence operator is directly added into the double-rotation equation, the modified double-rotation equation can be simplified into a Laplace equation, non-zero elements of each row can be reduced, and the calculation time of an iterative solver can be reduced; by adopting the forward modeling method, the maximum errors of the amplitude and the phase and the reference solution are small, and the precision is high; divergence correction is not required to be applied, and divergence conditions are automatically met; the algorithm of the present invention is more advantageous when the frequency is low (e.g. the maximum error of amplitude and phase is less than 10 at 1Hz, respectively-6(ii) a At 0.1Hz, the invention saves 63% of the running time and has outstanding effect compared with the more common CC-DC algorithm (the divergence correction algorithm is applied).
In addition to the objects, features and advantages described above, other objects, features and advantages of the present invention are also provided. The present invention will be described in further detail below with reference to the drawings.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this application, illustrate embodiments of the invention and, together with the description, serve to explain the invention and not to limit the invention. In the drawings:
FIG. 1 is a flow chart of a marine controllable source electromagnetic forward modeling method for a three-dimensional arbitrary anisotropic medium in example 1;
FIG. 2(a) is an initial global coordinate system;
FIG. 2(b) is a coordinate system after rotation about the z-axis based on FIG. 2 (a);
FIG. 2(c) is a coordinate system after rotation about the x-axis based on FIG. 2 (b);
FIG. 2(d) is a coordinate system after rotation about the z-axis based on FIG. 2 (c);
FIG. 3(a) is a sparse pattern diagram of a VTI anisotropic model coefficient matrix in the CC-DC algorithm;
FIG. 3(b) is a sparse pattern diagram of the coefficient matrix of the VTI anisotropic model in example 1;
FIG. 4 is a vertical sectional view of a deepwater complex VTI marine complex oil reservoir model in example 2;
FIG. 5 is a graph comparing the algorithm for applying divergence correction with the new algorithm set forth in example 2;
FIG. 6(a) is a graph comparing the calculation times of the algorithm applying the divergence correction and the algorithm proposed in example 2;
FIG. 6(b) is a graph comparing the number of iterations of the algorithm applying divergence correction and the algorithm set forth in example 2;
FIG. 7 is a comparison of the convergence process of these algorithms.
Detailed Description
Embodiments of the invention will be described in detail below with reference to the drawings, but the invention can be implemented in many different ways, which are defined and covered by the claims.
Example 1:
the forward calculation method of the embodiment is shown in fig. 1 in detail, and specifically includes the following steps:
reading the forward model, and obtaining the edge length of each grid unit along the x direction, the y direction and the z direction, the conductivity value of each grid unit and a frequency list. The reading forward model specifically comprises the following steps: generating a related three-dimensional complex geological model file (forward model) through a model generating function according to the requirements of the simulated underground geologic body; the generated geological model file (forward model) is read through a main program and stored in the form of data. The list of frequencies refers to the different frequencies entered, for example: 0.1Hz, 1Hz, 10Hz, are also read by the main program and stored in the form of data.
The forward model of the present embodiment is a smaller model of the marine VTI medium in order to understand the distribution of the non-zero elements and the change of the condition number of the model.
And step two, forming a double rotation degree equation expression 2 by combining the expression 1) of the Maxwell equation set:
Figure BDA0002498700020000051
Figure BDA0002498700020000052
wherein:
Figure BDA0002498700020000053
is Hamiltonian; e is an electric field; h is a magnetic field; j. the design is a squaresIs the external current density vector; ρ is the free charge; is the dielectric permittivity; omega is angular frequency; σ is the conductivity tensor; μ is the magnetic permeability in vacuum; i denotes the imaginary part, i2＝-1；
The conductivity tensor σ is expressed by expression 3):
Figure BDA0002498700020000054
wherein: sigmaxxConductivity values representing the x-axis direction; sigmayyConductivity values representing the y-axis direction; sigmazzConductivity values in the z-axis direction; sigmaxyAnd σyxConductivity values representing the combined action of the x-axis and the y-axis are specifically: sigmaxyDenotes the current density formed in the y direction by applying an electric field in the x direction, and σyxRepresents the current density formed in the x direction by the electric field applied in the y direction;σxzand σzxConductivity values representing the x-axis and z-axis co-action are specifically: sigmaxzDenotes the current density formed in the z direction by applying an electric field in the x direction, and σzxRepresents the current density formed in the x direction by the electric field applied in the z direction; sigmayzAnd σzyConductivity values representing the combined action of the y-axis and the z-axis are: sigmayzDenotes the current density formed in the z direction by applying an electric field in the y direction, and σzyIndicating the electric field applied in the z-direction and the resulting current density in the y-direction.
From the above, it can be seen that: the conductivity tensor is a symmetric, positive definite matrix, determined by 6 elements. Any one of the matrixes can obtain the conductivity tensor under the main reference system through three elementary transformations (matrix rotation), and only the conductivity diagonal element sigma of the matrix is used in the processxx、σyyAnd σzzAnd the other elements are zero, so that the conversion is also convenient for generating the forward model. FIG. 2 is a rotation parameter definition and rotation process from a main reference frame to a general measurement reference frame (i.e., Euler rotation diagram in the global coordinate system x-y-z, see FIGS. 2(a) - (d) for details), φ, θ and
Figure BDA0002498700020000066
the anisotropy run angle φ (rotation about the z-axis, where x becomes x ', y becomes y', z becomes z '), the tilt angle θ (rotation about the x-axis, where x' becomes x ", y 'becomes y", z' becomes z "), and the tilt angle θ, respectively, are in turn the corresponding rotation angles
Figure BDA0002498700020000067
(rotation about the z-axis, where x "becomes x '", y "becomes y '", and z "becomes z '"), the conductivity tensor σ is expressed using expression 4):
Figure BDA0002498700020000061
wherein: sigmaD＝diag(σxx,σyy,σzz) Diag is a diagonal function; rzAnd RxAre the sine and cosine of the angle between the principal axis and the global x-y-z coordinate system,
Figure BDA0002498700020000062
α represents variables;
for the tilt anisotropy model, expression 4) transforms to expression 5):
σ＝RTσDR 5)；
wherein: the elements of R are the sine and cosine of the angle between the principal axis and the global x-y-z coordinate system,
Figure BDA0002498700020000063
strike angle
Figure BDA0002498700020000064
And the inclination angle theta are both 0.
Step three, firstly, adding a gradient-divergence operator to the double-rotation equation expression 2) in the step two to form a modified control equation expression 6); discretizing the modified control equation expression 6), calculating an electric field value at the boundary of the forward model according to the forward model and the frequency list, and combining a field source item to obtain a discretized control equation expression 7):
Figure BDA0002498700020000065
Figure BDA0002498700020000071
wherein: ebA secondary field E calculated by a 1D geoelectric model as a primary field of the electric field EaSolving to obtain; eaThe secondary field of the electric field E represents the electric field vector value to be solved after discretization; lambda [ alpha ]1And λ2For the control factor (like the regularization factor in regularized inversion), λ1＝1/σ，λ 21/Δ σ, Δ σ represents the difference in conductivity between the actual conductivity model and the 1D geoelectric model; c represents a discretized rotation operator from unit edge to unitMapping of element areas;
Figure BDA0002498700020000072
representing the transposition of a discretized rotation operator C, representing the mapping from the cell area to the cell edge, G representing a gradient operator, representing a mapping from a node to an inner edge, Λ representing a diagonal matrix formed by a scaling factor, D representing a divergence operator, representing a mapping from an inner edge to an inner node, sigma representing a vector formed by conductivity elements, sigma being diag (sigma);
Figure BDA0002498700020000073
conductivity values on the edges Λ1Representing a diagonal matrix formed by the scaling factors; sigma1Vector representing the composition of the residual conductivity element, Σ1＝diag(Δσ)。
Step four, obtaining expressions of matrix coefficient A and constant b according to the discretization control equation expression 7), such as expression 8) and expression 9):
Figure BDA0002498700020000074
b＝diag(-iωμΔσ)Eb+GΛ1DΣ1Eb9)；
step five, solving the linear equation set Ae ═ b to obtain and output the electric field component Eax、EayAnd EazThe method specifically comprises the following steps:
step 5.1, setting an initial value ek＝e0Calculating the residual error rk＝b-AekWherein: e.g. of the type0An initial value representing an electric field imparted by the iterative solution; e.g. of the typekRepresenting the electric field value of the kth iteration; r iskRepresenting the residual value of the kth iteration;
step 5.2, adopting a double conjugate gradient iterative algorithm to carry out iterative loop, and giving a residual value r of the double conjugate gradient method0(ii) a If the residual error reaches the requirement 10-8Or exceeds the maximum convergence time by 3000 times, and outputs the electric field component Eax、EayAnd EazOutput electric field component Eax、EayAnd Eaz(ii) a Otherwise, k is taken as k +1, and the procedure returns to step 5.1.
By applying the technical scheme of the embodiment, the coefficient matrix A of the original equation is compared with the coefficient matrix of the modified control equation. For a VTI anisotropic model (6 × 6 × 8 grid), the sparse mode is shown in fig. 3, and fig. 3(a) is the result of applying divergence correction based on the dual rotation equation, i.e., CC-DC; FIG. 3(b) is a correction algorithm of the present invention, in which the gradient-divergence operator is directly added to our double rotation equation, i.e., CCGD. For the interleaved finite difference mesh, the double rotation operator has 13 non-zero elements in the standard second order mesh, i.e., there are at most 13 non-zero elements per row of our original equation. However, for the control equation after correction, in an air layer (where the conductivity is constant) or a non-abnormal body region, the corrected double rotation degree equation can be simplified into a laplace equation, and non-zero elements of each row can be reduced. The reduced number of matrix non-zero elements (from 6852 to 6077) may reduce the computation time of the iterative solver. More importantly, the condition number of the matrix is also greatly reduced. For the small model problem of FIGS. 3(a) - (b), the condition number drops from 1.851e +13 to 840.52. Thus, the Krylov solver can achieve convergence with fewer iterations.
Example 2:
the difference from example 1 is that the forward calculation uses a different model.
The forward modeling of the embodiment is a deep water complex VTI marine oil reservoir modeling, as shown in fig. 4. The seawater is isotropic, and the conductivity is 3.2 s.m-1The sea depth in the figure is 1000m, the height of the transmitting source from the sea bottom is 130m, the receiver is positioned on the surface of the sea bottom, the thickness of the sediment is 1.5km, the thickness of the oil reservoir is 100m, the model is discretized by using a 60 × 60 × 35 unit grid 58 × 58 × 27km3When the electric field source is located in the center portion of the model, the Dirichlet condition can be satisfied at the boundary. Using a discrete model of FD grid, the FD grid has a size of 100m in the x and y directions and a thickness of 80m in the z direction. The first layer of sediment conductivity is sigma considering the anisotropy of sediment and oil reservoirxx＝σyy＝1s·m-1,σzz＝0.5s·m-1Uniform half-space sediment conductivity of the seafloorxx＝σyy＝0.5s·m-1，σzz＝0.25s·m-1. The intermediate oil reservoir is a VTI model and the conductivity is sigmaxx＝σyy＝0.01s·m-1，σzz＝0.0025s·m-1。
To verify the correctness of the algorithm of the present invention, the forward simulation result of the method of the present invention (CCGD algorithm) was compared with the forward simulation result of the CC-DC algorithm with divergence correction applied, as shown in FIG. 5, the maximum errors of the amplitude and phase were each less than 10 at a frequency of 1Hz-6The high accuracy of the algorithm of the invention is proved.
Then, the method of the present invention (CCGD algorithm) in example 2 compares the calculation time and the number of iterations with the CC-DC algorithm, as shown in fig. 6, and in addition to the convergence process diagram of fig. 7, the three algorithms respectively represent a method directly based on the bispin equation (CC algorithm), a method of applying divergence correction to the bispin equation (CC-DC algorithm), and a method of the present invention of adding a gradient-divergence operator to the bispin equation (CCGD algorithm). Therefore, the following steps are carried out: (1) the ocean controllable source electromagnetic forward modeling method for the three-dimensional arbitrary anisotropic medium automatically meets the divergence condition; (2) the frequency commonly used in the ocean controllable source simulation is 0.1Hz-10Hz, the method is directly based on a double rotation equation (CC algorithm), and due to the existence of anisotropy, a coefficient matrix is more complex and is not converged in the frequency range, so that the CC algorithm is not shown; (3) furthermore, in this frequency range, the inventive method (CCGD algorithm) achieves convergence with the more common CC-DC algorithm, but the inventive method (CCGD algorithm) is faster than the CC-DC algorithm, as shown in FIG. 6 (b); (4) the method of the invention (CCGD algorithm) achieves a significant advantage as the frequency decreases, saving 63% of the run time at 0.1Hz compared to the more common CC-DC algorithm, as shown in FIG. 6(a), where 37s is used and 101s is used.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (6)

1. A marine controllable source electromagnetic forward modeling method for a three-dimensional arbitrary anisotropic medium is characterized by comprising the following steps:
reading a forward model to obtain initial parameters;
and step two, forming a double rotation degree equation expression 2 by combining the expression 1) of the Maxwell equation set:
Figure FDA0002498700010000011
▽×▽×E+iωμσE＝-iωμJs2)；
wherein ▽ is Hamiltonian, E is electric field, H is magnetic field, and JsIs the current density vector of the external current; ρ is the free charge; is the dielectric permittivity; omega is angular frequency; σ is the conductivity tensor; μ is the magnetic permeability in vacuum; i denotes the imaginary part, i2＝-1；
Step three, firstly, adding a gradient-divergence operator to the double-rotation equation expression 2) in the step two to form a modified control equation expression 6); discretizing the modified control equation expression 6) to obtain a discretized control equation expression 7):
▽×▽×Ea-λ1▽(▽·σEa)+iωμσEa＝-iωμΔσEb+λ2▽(▽·ΔσEb) 6)；
Figure FDA0002498700010000012
wherein: ebIs a primary field of an electric field E, EaIs a secondary field of the electric field E and represents the electric field vector value to be solved after discretization, Eb+Ea＝E；λ1And λ2Is a control factor; c represents the discretized rotationAn operator, which is the mapping from the unit edge to the unit area;
Figure FDA0002498700010000015
representing the transposition of a discretized rotation operator C, representing the mapping from a unit area to a unit edge, G representing a gradient operator, realizing the mapping from a node to an internal edge, Λ representing a diagonal matrix formed by a scaling factor, D representing a divergence operator, representing the mapping from the internal edge to the internal node, sigma representing a vector formed by conductivity elements, sigma being diag (sigma), and diag being a diagonal function;
Figure FDA0002498700010000013
conductivity values on the edges Λ1Representing a diagonal matrix formed by the scaling factors; sigma1Vector representing the composition of the residual conductivity element, Σ1Δ σ represents a difference between the actual conductivity model and the electrical conductivity of the earth model;
step four, obtaining expressions of matrix coefficient A and constant b according to the discretization control equation expression 7), such as expression 8) and expression 9):
Figure FDA0002498700010000014
b＝diag(-iωμΔσ)Eb+GΛ1DΣ1Eb9)；
step five, solving the linear equation set Ae ═ b to obtain the electric field component Eax、EayAnd Eaz。
2. The marine controlled-source electromagnetic forward modeling method for the three-dimensional arbitrary anisotropic medium according to claim 1, characterized in that the parameters in the first step are conductivity values of each grid cell, frequency lists, edge lengths in the x, y and z directions, and information of all grid nodes; the grid node information is a grid node number.
3. The marine controlled-source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic media according to claim 2, characterized in that the conductivity tensor σ in the second step is expressed by expression 3):
Figure FDA0002498700010000021
wherein: sigmaxxConductivity values representing the x-axis direction; sigmayyConductivity values representing the y-axis direction; sigmazzConductivity values in the z-axis direction; sigmaxyAnd σyxConductivity values representing the x-axis and y-axis co-action; sigmaxzAnd σzxConductivity values representing the x-axis and z-axis co-action; sigmayzAnd σzyThe values of the electrical conductivity are shown for the y-axis and z-axis co-action.
4. The marine controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic media according to claim 2, characterized in that in the third step: the geoelectric model is 1D geoelectric model, EbIs a primary field obtained by calculation through a 1D geoelectric model, and a secondary field E is obtained by solvinga；λ1＝1/σ，λ2＝1/Δσ。
5. The method for the controlled-source electromagnetic forward modeling of the three-dimensional arbitrary anisotropic medium sea according to claim 4, wherein the step five specifically comprises the following steps:
step 5.1, setting an initial value ek＝e0Calculating the residual error rk＝b-AekWherein: e.g. of the type0An initial value representing an electric field imparted by the iterative solution; e.g. of the typekRepresenting the electric field value of the kth iteration; r iskRepresenting the residual value of the kth iteration;
step 5.2, carrying out iterative loop, and outputting the electric field component Ea if the residual error meets the requirementx、EayAnd Eaz(ii) a Otherwise, k is taken as k +1, and the procedure returns to step 5.1.
6. Three-dimensional arbitrary anisotropy according to claim 5The marine controllable source electromagnetic forward modeling method of the sexual medium is characterized in that a double conjugate gradient iterative algorithm is adopted in an iterative loop; residual value r of given double conjugate gradient method0(ii) a If the residual error reaches the requirement 10-8Or exceeds the maximum convergence time by 3000 times, and outputs the electric field component Eax、EayAnd Eaz。
CN202010426036.1A 2020-05-19 2020-05-19 Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium Active CN111611737B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010426036.1A CN111611737B (en) 2020-05-19 2020-05-19 Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010426036.1A CN111611737B (en) 2020-05-19 2020-05-19 Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium

Publications (2)

Publication Number Publication Date
CN111611737A true CN111611737A (en) 2020-09-01
CN111611737B CN111611737B (en) 2022-05-20

Family

ID=72200706

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010426036.1A Active CN111611737B (en) 2020-05-19 2020-05-19 Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium

Country Status (1)

Country Link
CN (1) CN111611737B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113627027A (en) * 2021-08-17 2021-11-09 四川大学 Method and system for simulating electromagnetic field value of non-trivial anisotropic medium
CN114912310A (en) * 2022-04-11 2022-08-16 中南大学 Three-dimensional magnetotelluric numerical simulation method based on regularization correction equation
CN116842813A (en) * 2023-09-04 2023-10-03 中南大学 Three-dimensional geoelectromagnetic forward modeling method

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080071709A1 (en) * 2006-08-22 2008-03-20 Kjt Enterprises, Inc. Fast 3D inversion of electromagnetic survey data using a trained neural network in the forward modeling branch
CN106980736A (en) * 2017-04-11 2017-07-25 吉林大学 A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium
CN107845141A (en) * 2017-11-27 2018-03-27 山东大学 A kind of transient electromagnetic three-dimensional FDTD forward modeling multi-resolution meshes subdivision methods
CN108509693A (en) * 2018-03-13 2018-09-07 中南大学 Three-dimensional frequency domain controllable source method for numerical simulation
CN110058315A (en) * 2019-05-29 2019-07-26 中南大学 A kind of three dimensional anisotropic radio frequency magnetotelluric self-adapting finite element forward modeling method
CN110346835A (en) * 2019-07-22 2019-10-18 中国科学院地球化学研究所 Magnetotelluric forward modeling method, forward modeling system, storage medium and electronic equipment
CN110826283A (en) * 2019-11-15 2020-02-21 中南大学 Preprocessor and three-dimensional finite difference electromagnetic forward modeling calculation method based on preprocessor

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080071709A1 (en) * 2006-08-22 2008-03-20 Kjt Enterprises, Inc. Fast 3D inversion of electromagnetic survey data using a trained neural network in the forward modeling branch
CN106980736A (en) * 2017-04-11 2017-07-25 吉林大学 A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium
CN107845141A (en) * 2017-11-27 2018-03-27 山东大学 A kind of transient electromagnetic three-dimensional FDTD forward modeling multi-resolution meshes subdivision methods
CN108509693A (en) * 2018-03-13 2018-09-07 中南大学 Three-dimensional frequency domain controllable source method for numerical simulation
CN110058315A (en) * 2019-05-29 2019-07-26 中南大学 A kind of three dimensional anisotropic radio frequency magnetotelluric self-adapting finite element forward modeling method
CN110346835A (en) * 2019-07-22 2019-10-18 中国科学院地球化学研究所 Magnetotelluric forward modeling method, forward modeling system, storage medium and electronic equipment
CN110826283A (en) * 2019-11-15 2020-02-21 中南大学 Preprocessor and three-dimensional finite difference electromagnetic forward modeling calculation method based on preprocessor

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李健: "《一种频率域三维有限差分电磁正演的高效预处理器》", 《 2019年中国地球科学联合学术年会》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113627027A (en) * 2021-08-17 2021-11-09 四川大学 Method and system for simulating electromagnetic field value of non-trivial anisotropic medium
CN113627027B (en) * 2021-08-17 2023-04-11 四川大学 Method and system for simulating electromagnetic field value of non-trivial anisotropic medium
CN114912310A (en) * 2022-04-11 2022-08-16 中南大学 Three-dimensional magnetotelluric numerical simulation method based on regularization correction equation
CN114912310B (en) * 2022-04-11 2024-04-12 中南大学 Three-dimensional magnetotelluric numerical simulation method based on regularization correction equation
CN116842813A (en) * 2023-09-04 2023-10-03 中南大学 Three-dimensional geoelectromagnetic forward modeling method
CN116842813B (en) * 2023-09-04 2023-11-14 中南大学 Three-dimensional geoelectromagnetic forward modeling method

Also Published As

Publication number Publication date
CN111611737B (en) 2022-05-20

Similar Documents

Publication Publication Date Title
CN111611737B (en) Ocean controllable source electromagnetic forward modeling method for three-dimensional arbitrary anisotropic medium
da Silva et al. A finite element multifrontal method for 3D CSEM modeling in the frequency domain
Jahandari et al. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids
Soloveichik et al. Finite-element solution to multidimensional multisource electromagnetic problems in the frequency domain using non-conforming meshes
Kordy et al. 3-D magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on SMP computers–Part I: forward problem and parameter Jacobians
AU2011339017B2 (en) Constructing geologic models from geologic concepts
Key et al. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example
Jahandari et al. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids
AU2003272233B2 (en) Method of conditioning a random field to have directionally varying anisotropic continuity
CN106980736B (en) A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium
Davis et al. Efficient 3D inversion of magnetic data via octree-mesh discretization, space-filling curves, and wavelets
Jaysaval et al. Fast multimodel finite-difference controlled-source electromagnetic simulations based on a Schur complement approach
Persova et al. Improving the computational efficiency of solving multisource 3-D airborne electromagnetic problems in complex geological media
Peng et al. 3-D marine controlled-source electromagnetic modeling in electrically anisotropic formations using scattered scalar–vector potentials
Ouyang et al. Iterative magnetic forward modeling for high susceptibility based on integral equation and Gauss-fast Fourier transform
Li et al. A finite‐element time‐domain forward‐modelling algorithm for transient electromagnetics excited by grounded‐wire sources
Ye et al. 3-D adaptive finite-element modeling of marine controlled-source electromagnetics with seafloor topography based on secondary potentials
Rong‐Hua et al. 3‐D INVERSION OF FREQUENCY‐DOMAIN CSEM DATA BASED ON GAUSS‐NEWTON OPTIMIZATION
Mallesh et al. 3D gravity analysis in the spatial domain: Model simulation by multiple polygonal cross-sections coupled with exponential density contrast
Huang et al. 3D anisotropic modeling and identification for airborne EM systems based on the spectral-element method
Elías et al. Three-dimensional modelling of controlled source electro-magnetic surveys using non-conforming finite element methods
Wang et al. Efficient 2D modeling of magnetic anomalies using NUFFT in the Fourier domain
Zhu et al. 3D unstructured spectral element method for frequency-domain airborne EM forward modeling based on Coulomb gauge
Xie et al. Surface-geometry and trend modeling for integration of stratigraphic data in reservoir models
Jahandari et al. Forward modeling of direct-current resistivity data on unstructured grids using an adaptive mimetic finite-difference method

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