CN115186520A - Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium - Google Patents
Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium Download PDFInfo
- Publication number
- CN115186520A CN115186520A CN202211106686.3A CN202211106686A CN115186520A CN 115186520 A CN115186520 A CN 115186520A CN 202211106686 A CN202211106686 A CN 202211106686A CN 115186520 A CN115186520 A CN 115186520A
- Authority
- CN
- China
- Prior art keywords
- formula
- source
- dipole source
- layer
- vector
- 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
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
-
- 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
- G06F17/13—Differential equations
-
- 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/15—Correlation function computation including computation of convolution operations
-
- 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)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Computing Systems (AREA)
- Operations Research (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a time-frequency electromagnetic response simulation method of a vector dipole source in a horizontal layered geodetic medium, which is characterized in that a propagation coefficient matrix of a vector potential function caused by a universal electric dipole source and a universal magnetic dipole source in a rectangular coordinate system is established, and an electromagnetic field component of the earth surface is obtained through recursion solution; respectively calculating the electromagnetic fields of the horizontal dipole source and the vertical dipole source, and synthesizing to obtain the electromagnetic field response of any vector dipole source; the electric dipole source and the magnetic dipole source can be used for any direction in a horizontal layered ground; and calculating after obtaining the amplitude function of the dipole source in the uniform full-space potential function to obtain the time-frequency electromagnetic response of the dipole source in the horizontal layered earth medium. The invention can be applied to simulating electromagnetic disturbance of earthquake precursors in the laminar earth, an underground nuclear pollution source or an electromagnetic wave field radiated by an electric or magnetic vector dipole source of a coal mine roadway, and has good applicability and high application value.
Description
Technical Field
The invention provides a method for simulating and calculating an electromagnetic field generated by any vector electric dipole source and magnetic dipole source in a horizontal layered earth, in particular to a method for simulating time-frequency electromagnetic response of the vector dipole source in a horizontal layered earth medium, which simulates through a potential function propagation coefficient matrix and solves time domain and frequency domain electromagnetic field response and belongs to the technical field of earth electromagnetic application.
Background
The existing earthquake electromagnetic observation data show that the earthquake can generate electromagnetic radiation in the inoculation process to cause disturbance of an earth electromagnetic field in a larger range. The interpretation of the electromagnetic signals generated before the earthquake happens and the damage is caused has important significance for the earthquake prediction research. The source of the electromagnetic wave radiated before the earthquake rupture can be regarded as the superposition of an electric dipole source or a magnetic dipole source, and based on the superposition, the simulation of the electromagnetic response generated by a vector dipole source in the laminated earth is the basis for researching the electromagnetic image of the earthquake precursor. However, no algorithm for calculating the electromagnetic field quantity of the full space time domain and the frequency domain of the vector dipole source which is suitable for any direction exists in the work for calculating the electromagnetic field response of the dipole source in the laminated medium at present. Wait (1951) describes in The literature (Wait J R. The magnetic dipole over The magnetic dipole structured earth. Can J Phys, 195l,29: 577-592.) that The continuity of The horizontal electromagnetic field at The interface solves The propagation coefficient of The vector potential function caused by The magnetic source and realizes The calculation of The response of The magnetic dipole source in The horizontal layered medium. Xujianhua et al (1994) in literature (Xujianhua, zhudeby, huwenbao, lepeng. The coefficient recurrence relation of dipole field in multilayer medium. Oil geophysical prospecting, 1994, 29 (1): 69-74.) derived the layer system coefficient recurrence relation and the transfer matrix of the vector potential function caused by the magnetic dipole source in a similar way to Wait, and can be used to solve the electromagnetic field response generated by the presence of horizontal or vertical dipoles in the multilayer medium. The existing electromagnetic field response acquisition technology can only be used for a specific source or layer structure, has limitation on calculating the actual electromagnetic field in the geophysical application field, such as well logging and the electromagnetic field on the earth surface, and is difficult to solve the simulation calculation of the electromagnetic field response of seismic precursor electromagnetic radiation and a complex dipole source.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a simulation method of time-frequency electromagnetic response of a dipole source in a horizontal layered earth medium, which respectively calculates the electromagnetic fields of a horizontal dipole source and a vertical dipole source, then synthesizes and obtains the electromagnetic field response of any vector dipole source, can be applied to explain the electromagnetic disturbance phenomenon generated by the vector dipole source in seismic precursors in a layered earth, and further can be applied to the prediction and identification of underground nuclear pollution sources or electromagnetic radiation generated in coal mine tunnels in major earthquakes.
The invention improves the existing calculation method of the response of the magnetic dipole source in the horizontal layered medium and the method for solving the electromagnetic field response generated by the existence of horizontal or vertical dipoles in the multilayer medium by utilizing the layer system coefficient recurrence relation and the transfer matrix of the vector potential function caused by the electric and magnetic dipole sources. For different types of dipole sources, the time-frequency electromagnetic response of the dipole source in the horizontal laminar ground medium can be obtained only by obtaining the amplitude function of the dipole source in the uniform full-space potential function and then calculating.
The method of the invention comprises the following steps:
A. establishing a general expression of the electromagnetic response of a dipole source in a horizontal laminar earth medium:
a = represents a vector potential function of an electromagnetic field generated by an electric dipole source in any direction in a horizontal layerA z , Is a unit vector of the z-axis direction in a rectangular coordinate system; the vector potential function of the electromagnetic field generated by a magnetic dipole source is denoted as F =F z I.e. vector potential functions A and F are constrained to have only z-componentsA z AndF z 。
the components of the electromagnetic field excited by the electric dipole source in any direction in the horizontal layer are respectively expressed as formula 1 and formula 2:
In the formula, complex conductivity=σ+iωε,iRepresents the unit of an imaginary number and is,σas the electrical conductivity of the medium, the conductivity of the medium,εis the dielectric constant;E x 、E y 、E z respectively under a rectangular coordinate systemx、y、zAn axial electric field component;H x 、H y 、H z are respectively asx、y、zA magnetic field component of the shaft; whereinH z Is a TM polarized field;
the components of the electromagnetic field caused by the magnetic dipole source in any direction in the horizontal layer are expressed as formula 3 and formula 4, respectively:
Wherein the vertical component of the electric fieldE z The number of the carbon atoms is zero,E z is a TE polarization field;in order to be a resistance ratio of the resistor, ,irepresents the unit of an imaginary number and is,ω-a signal angular frequency;is medium magnetic permeability; in formulas 1 and 4kSatisfy the requirements ofIs the wave number of the medium;
B. calculating a vector potential function of an electromagnetic field, and respectively solving a solution of a primary field of a passive region potential function and a solution of a secondary field of an active layer potential function;
the method is considered in two cases, namely a solution of a vector potential function of a primary electromagnetic field of a passive region and a solution of a vector potential function of a secondary electromagnetic field of an active layer position;
B1. for the potential function of the passive region, a homogeneous differential equation can be solved to obtain:
In the formula, the upper wavy sign represents a fourier transform form of the corresponding potential function,A + andF + represents the corresponding uplink (+)zDirection) wave solution is performed,A - andF - then represents the corresponding downlink-zDirection) of the wave solution,k x andk y are respectivelyxAndythe number of spatial waves in a direction is,u 2 =k x 2 +k y 2 -k 2 ;
B2. solving a general solution of a potential function of the source layer:
for the source-containing layer, in addition to the complementary solution, a special solution of the non-homogeneous differential equation is added.
For any purposeNThe model of the layer of the earth,zin the axial direction; setting the center of the source as the origin of coordinates, the interface where the source is locatedz s =0, the interface is a virtual interface for computing the response of the source; is directly above the source-containing interfacezThe sequence number of the direction is marked as i from bottom to top 1 ,i 2 ,…,i L (ii) a Formation parameter corresponds to conductivityσ ii Dielectric constant ofε ii Thickness of the layerh ii (ii) a Is immediately negative below the source-containing interfacezThe sequence numbers of the layers in the directions are marked as j from top to bottom 1 ,j 2 ,…,j M The formation parameter corresponds to conductivityσ jj Dielectric constant ofε jj And thickness of each layerh jj (ii) a Is provided withA p AndF p for the amplitude of the source field, the characteristic solution of the non-homogeneous differential equation with source time is expressed as formula 7 and formula 8 according to the plane wave solution of the uniform full space:
In the form of a Fourier transform of a vector potential functionzComponent(s) ofAndproducing attenuation both above and below the source;Nin the course of the earth, the water is discharged,nis the layer number, the source-containing layer isn=i 1 Andn=j 1 ;
combining the special solutions of the sources of formulas 7 and 8 with the solution of the potential function of the passive region in B1, the hankel integral form of the potential function of the active layer is obtained, and is expressed as formulas 9 and 10:
In the formula (I), the compound is shown in the specification,A n + 、A n - 、F n + 、F n - in order to determine the coefficient to be determined,,,by coefficient factorβ n To distinguish between active and inactive regions,β n the values are as follows:
C. according to the potential function of the underground vertical electric dipole source, calculating to obtain the amplitude of the potential function:
will lie at the origin of coordinates with an electrode moment ofIdl The vertical electric dipole source density of (a) is expressed as:
In the formula (I), the compound is shown in the specification,is a unit vector in the z-axis direction; δ (d)x)、δ(y)、δ(z) Is a dirac function, is 1 only at the origin of coordinates (0, 0), otherwise is 0; the magnitude of the potential function for a vertical electric dipole source is expressed as:
The hankel integral form of the potential function of the vertical electric dipole source in any layer in the transform domain is obtained by the formula 9:
Propagation coefficientA n0 + AndA n0 - solving by a recursion relation; the recurrence relation includes:
C1. determining the top and bottom layer transfer matrices from the source and the medium parameters of each layer to obtain the propagation coefficients of the top and bottom layers, which are expressed as equations 14 and 15:
In the formula (I), the compound is shown in the specification,q2 x 2 transfer matrixQ4 coefficients, the coefficient value is determined by each layer of medium parameters;s i as a source vectorS i 2 coefficients of,s j As a source vectorS j 2 coefficients ofCoefficient value is determined by source;
And (3) obtaining a transfer matrix expression form of the source layer according to the recursion of the earth lamellar model:
C2. Determining a layer-by-layer propagation coefficient;
setting a virtual layer interface at the position of a source, and dividing the earth model into two conditions of an upper source condition and a lower source condition;
U i Is expressed as:
Wherein:
Ith L The layers being top and bottom interfaces and extending upwardly to an infinite extent, i.e. spaceA Li - If =0 has no downstream wave, then:
Under-source TM polarization modeU j Expressed in matrix form:
Similarly, TE polarizationThe transfer matrix form of the source layer of the mode is the same, whereinIs composed of,A p Is composed ofF p 。
C3. And (3) recursing layer by layer from the propagation coefficient matrixes of the top layer and the bottom layer to the interior of the model, obtaining the propagation coefficient of the source layer from the formulas 16-30, and calculating the propagation coefficient of each layerA n + AndA n - ,nis a layer number;
substituting the propagation coefficient into formula 13 to obtain potential functions of each layer;
substituting potential functions into equations 3-6 to obtain kernel functions of components of the electromagnetic field;
then, a digital filtering algorithm is applied to realize the Hankel transformation, and the frequency domain response of each field component is obtained;
D. establishing a potential function of a horizontal electric dipole source in the ground:
the current source density of a horizontal dipole source is expressed as:
In the formula (I), the compound is shown in the specification,is thatxUnit vector of direction of axis, δ: (x)、δ(y)、δ(z) Is a dirac function, is 1 only at the origin of coordinates (0, 0), otherwise is 0;
the potential function of a horizontal electric dipole source contained in any layer is expressed as:
In the formula (I), the compound is shown in the specification,the symbol is expressed asz>When the number 0 is a negative number,z<when 0, the plus number is taken;
the propagation coefficients are:
Using the propagation coefficient to obtain a potential function, then substituting the equations from 32 to 36 for the equations from 3 to 6 to obtain the kernel functions of each component of the electromagnetic field, and then applying a digital filtering algorithm to realize the Hankel transformation to obtain any layernFrequency domain response of each field quantity excited by the middle-level electric dipole source;
E. carrying out simulation calculation to obtain a vector electric dipole source field in any direction in the ground; the method comprises the following steps:
E1. decomposing a field source;
electric dipole momentP e =IdlAzimuth angleAnd the angle of inclinationThree parameters define any electric dipole in the ground, and the electric dipole moment is decomposed into horizontal componentsAnd a vertical componentP z Respectively is as follows:
For the layered medium, the field of the horizontal electric dipole source in any direction has rotational symmetry, and the field of the horizontal electric dipole source in any direction can be obtained through coordinate rotation; for the vertical electric dipole source, synthesizing through the calculation formula of the vertical dipole source response deduced in the step C;
for an azimuth angle ofThe coordinate parameters after rotation of the horizontal electric dipole source of (3) are as follows:
E2. Electromagnetic field synthesis;
an electric dipole moment ofP e Field vector of response of time-horizontal electric dipole source: (And) And field vector of vertical electric dipole source response (E) z And H z ) And (3) synthesizing to obtain the electromagnetic field components of any electric dipole source in the ground as follows:
E3. Performing inverse Fourier transform on the synthesized electromagnetic field to obtain an electromagnetic field component of a time domain;
therefore, the simulation calculation of the time-frequency electromagnetic response of the dipole source in the horizontal laminar earth medium is realized.
Compared with the prior art, the invention has the following beneficial technical effects:
the invention provides a calculation method of time domain and frequency domain electromagnetic field response of a vector dipole source in a horizontal layered earth, and derives a recursion formula and a solution algorithm suitable for an electric dipole source, a magnetic dipole source and a multilayer propagation coefficient; and a computing technology for realizing electromagnetic response of the line source and the surface source by utilizing multi-source combination. The invention can be applied to simulating the electromagnetic disturbance of earthquake precursors in the layered earth, the electromagnetic wave field radiated by an electric or magnetic vector dipole source of an underground nuclear pollution source or a coal mine tunnel, can be used for explaining the electromagnetic disturbance phenomenon generated by the vector dipole source in the earthquake precursors in the layered earth, and is further used for predicting and identifying the electromagnetic radiation generated in the underground nuclear pollution source or the coal mine tunnel in the earthquake, thereby having higher applicability and application value.
Drawings
FIG. 1 is a flow chart diagram of a time-frequency electromagnetic response simulation calculation method of a subsurface layered medium dipole source provided by the invention.
FIG. 2 is a schematic illustration of layers of a vertical dipole source in a horizontal laminar earth medium;
wherein the content of the first and second substances,AandFrepresenting the propagation coefficients of the vector potential functions induced by electric and magnetic dipole sources, the upper subscript + representing the electromagnetic waves going upwards and the upper subscript-representing the electromagnetic waves going downwards.
FIG. 3 is an exploded schematic view of a vector electric dipole source;
wherein the content of the first and second substances,andrespectively representing the azimuth and inclination of the vector dipole source,P e dipole source distance, decomposable into horizontal componentsAnd the perpendicular componentP z 。
FIG. 4 is a schematic view of the structure and parameters of a 4-layer earth model employed in an embodiment of the present invention;
wherein the content of the first and second substances,which is representative of the magnetic permeability of the medium,σthe conductivity of the medium is represented, and the lower corner marks the surface layer number.
FIG. 5 is an x-y (horizontal) profile of a medium-earth vector electric dipole source in response to 10Hz amplitude for each field magnitude,
wherein (a) - (c) are the x, y and z components of the electric field frequency domain response, respectively; (d) - (f) the x, y and z components of the magnetic field frequency domain response, respectively.
FIG. 6 is a graph of the x-z (vertical) profile of a medium-earth vector electric dipole source in response to 10Hz amplitude for each field magnitude,
wherein (a) - (c) are the x, y and z components of the electric field frequency domain response, respectively; (d) - (f) are the x, y and z components of the magnetic field frequency domain response, respectively.
FIG. 7 is an x-z (vertical) profile of the amplitude of each field quantity of 0.1s of a vector electric dipole source in the ground,
wherein (a) - (c) are the x, y and z components of the electric field response, respectively; (d) - (f) are the x, y and z components of the magnetic field response, respectively.
Detailed Description
The invention will be further described by way of examples, without in any way limiting the scope of the invention, with reference to the accompanying drawings.
The invention provides a simulation method of time-frequency electromagnetic response of a dipole source in a horizontal layered earth medium, which is suitable for an electric dipole source, a magnetic dipole source, a multilayer propagation coefficient recursion formula and a solving algorithm by deriving and realizing the simulation calculation of line source and surface source electromagnetic response by utilizing multi-source combination and has higher applicability and application value.
In specific implementation, a model of space atmosphere and stratum is set, the azimuth angle and inclination angle of the vector electric dipole source and the coordinate position of the source are set, and the position of the source is set in this case (x s ,y s ,z s ) And setting parameters of each layer. The invention relates to a method for simulating time-frequency electromagnetic response of a dipole source in a horizontal laminar earth medium, which comprises the following steps:
A. establishing a general expression for the electromagnetic response of a dipole source in a horizontal laminar earth medium:
because the interfaces of the layered geodetic media with different levels are respectively equal to the interfaces of the layered geodetic media with corresponding levelszPlane coincidence, so that partial differential equations of the surface electromagnetic field can be transformed into equations of equal-degreezOrdinary differential equation of plane. The general solution to the boundary problem is the sum of the idiomatic solution of the non-homogeneous differential equation and the complementary solution of the homogeneous equation. Let the vector potential function of the electromagnetic field generated by an electric dipole source be A =A z ,In a rectangular coordinate systemzThe unit vector in the direction of the axis,A z is the z component of A(ii) a And the vector potential function of the electromagnetic field generated by the magnetic dipole source is denoted as F =F z I.e. vector potential functions A and F are constrained to have only z-componentsA z AndF z . The components of the electromagnetic field excited by the electric dipole source are then written as:
In the formula, the complex conductivity=σ+iωε,σ-medium conductivity (S/m);εfor dielectric constant, takeε=ε 0 =8.854*10 -12 F/m。E x 、E y 、E z Respectively under a rectangular coordinate systemx、y、zThe component of the electric field in the axial direction,H x 、H y、 H z are respectively asx、 y、zMagnetic field component of the shaft due toH z The (vertical component of the magnetic field) is zero,H z referred to as the TM polarization field. Similarly, the components of the electromagnetic field induced by a magnetic type dipole source are respectively expressed as:
Wherein, due toE z (vertical component of electric field) is zero,E z known as the TE polarization field.In order to be a resistance ratio of the resistor,,irepresents an imaginary unit;ω-a signal angular frequency;is medium magnetic permeability; for general applications, take= 0 =4π*10 -7 H/m。
B. Calculating a vector potential function of an electromagnetic field of the electric dipole source and an electromagnetic field of the magnetic dipole source;
the method is divided into two cases, namely the solution of a vector potential function of a primary electromagnetic field of a passive region (homogeneous partial differential equation complementary solution) and the solution of a vector potential function of a secondary electromagnetic field of an active-layer-containing region (non-homogeneous partial differential equation specific solution);
B1. solving a homogeneous differential equation for the potential function of the passive region to obtain:
In the formula, the upper wavenumber represents the fourier transform form of the corresponding potential function,A + andF + represents the corresponding uplink (+)zDirection) wave solution is performed,A - andF - then represents the corresponding downlink (-zDirection) of the wave solution,k x andk y are respectivelyxAndythe number of spatial waves in a direction is,u 2 =k x 2 +k y 2 -k 2 . For passive regions of the laminar earth medium, i.e. except i in FIG. 2 1 Or j 1 The outer regions, which can be solved for in different layers, in combination with appropriate boundary conditions, can be derived the amplitude (constant) of each layer.
B2. Solving a general solution form of the potential function of the source containing layer:
for the source-containing layer, in addition to the complementary solution, a special solution of a non-homogeneous differential equation is added. Consider arbitrarilyNThe stratigraphic earth model is shown in figure 1,zaxially. Setting the center of the source as the origin of coordinates, the interface where the source is locatedz s =0, this interface is a virtual interface, mainly used to calculate the response of the source. Above the interface containing the source (positive)zDirection) are sequentially marked as i from bottom to top 1 ,i 2 ,…,i L The formation parameter corresponds to conductivityσ ii Dielectric constant ofε ii Thickness of the layerh ii (ii) a Below the source-containing interface (negative)zDirection) is sequentially marked as j from top to bottom 1 ,j 2 ,…,j M The formation parameter corresponds to conductivityσ jj Dielectric constant ofε jj And thickness of each layerh jj . Is provided withA p AndF p the amplitude of the source field is obtained by solving a uniform full-space plane wave, and a special solution of a non-homogeneous differential equation with source time is expressed as follows:
In the form of a Fourier transform of the vector potential functionzComponent(s) ofAndproducing attenuation both above and below the source, hereA p AndF p depending on the source,nrepresenting the horizon number. As shown in fig. 2NIn the stratigraphic earth model, the horizon containing source isn=i 1 Andn=j 1 . Combining the special solutions of the sources of formula 7 and formula 8 with the solution of the passive region potential function in B1, a hankel integral form of the potential function of the source-containing layer can be obtained, which is expressed as:
In the formula (I), the compound is shown in the specification,A n + 、A n - 、F n + 、F n - in order to determine the coefficient to be determined,,,here by coefficient factorsβ n Distinguishing between active and passive conditions, coefficientβ n The values are as follows:
C. calculating a potential function of a vertical electric dipole source in the earth:
a short current carrying wire, if its length is far less than the distance to the observation point, it can be equivalent to electric dipole source; if the electric dipole source is close to the observation point or the current-carrying conducting wire is long, the emission source of the current-carrying conducting wire can be regarded as the superposition of a plurality of dipole sources. At the origin of coordinates an electrode moment ofIdl Perpendicular electric dipole source density (A/m) 2 ) Can be expressed as:
In the formula (I), the compound is shown in the specification,is a unit vector in the z-axis direction, δ: (x)、δ(y)、δ(z) For the dirac function, it is 1 only at the origin of coordinates (0, 0), otherwise it is 0. The magnitude of the potential function for a vertical electric dipole source is expressed as:
The Hamkerr integral form of the potential function of the vertical electric dipole source contained in any layer in a transform domain is obtained by the formula 9:
Propagation coefficientA n + AndA n - the result is obtained from the recursion relation,β n see step B2. The recurrence relation includes:
C1. determining top layer and bottom layer transfer matrixes according to the source and each layer medium parameter to obtain the propagation coefficients of the top layer and the bottom layer:
In the formula (I), the compound is shown in the specification,q2 x 2 transfer matrixQ4 coefficients, the coefficient value is determined by each layer medium parameter;s i as a source vectorS i Two coefficients of,s j As a source vectorS j Two coefficients ofThe coefficient values are determined by the source.
And (3) obtaining a transfer matrix expression form of the source layer according to the ground laminar model recursion of the figure 2:
C2. Determining a layer-by-layer propagation coefficient recurrence relation, setting a virtual layer interface at the position of a source, and considering the earth model into two conditions of the upper part of the source and the lower part of the source.
The over-source TM polarization mode is represented in matrix form:
Ith L The layers being top and bottom interfaces and extending upwardly to an infinite extent, i.e. spaceA Li - If =0 has no downstream wave, then:
TM polarization mode under sourceU j Expressed in matrix form:
Similarly, the source layers of the TE polarization mode have the same transfer matrix form, and only need to be replacedIs composed of,A p Is composed ofF p 。
C3. And (3) carrying out successive recursion on the propagation coefficient matrixes of the top layer (formula 14) and the bottom layer (formula 15) to the interior of the model, obtaining the propagation coefficient of the source layer through formulas 16-30, and then calculating to obtain each layernPropagation coefficient ofA n + AndA n - ,nis a layer number; substituting the calculated propagation coefficient into formula 13 to obtain potential functions of each layer; finally, substituting the potential function into the formula 3 to the formula 6 to obtain the kernel function of each component of the electromagnetic field; then, a digital filtering algorithm is applied to realize the Hankel transformation, and the frequency domain response of each field component is obtained;
D. establishing a potential function of a horizontal electric dipole source in the ground:
the current source density of a horizontal dipole source can be expressed as:
In the formula (I), the compound is shown in the specification,is thatxUnit vector of direction of axis, δ: (x)、δ(y)、δ(z) For the dirac function, it is 1 only at the origin of coordinates (0, 0), otherwise it is 0. The potential function of a source containing a horizontal electric dipole in any layer can be expressed as:
In the formula (I), the compound is shown in the specification,the symbol represents whenz>When the number 0 is a negative number,z<when 0, the + number is taken. The propagation coefficients (levels) are:
Processing similar to the step C3, utilizing the propagation coefficient to obtain a potential function, then obtaining the kernel functions of all components of the electromagnetic field in the formulas 3 to 6, and finally applying a digital filtering algorithm to realize Hankel transformation to obtain any layernAnd responding to each field quantity frequency domain excited by the middle-level electric dipole source.
E. Simulating and calculating to obtain a vector electric dipole source field in any direction in the ground; the method comprises the following steps:
E1. field source decomposition
Electric dipole momentP e =IdlAzimuth angle of the magnetic fieldAnd the angle of inclinationThree parameters define an arbitrary electric dipole, as shown in FIG. 3, which can be decomposed into horizontal componentsAnd the perpendicular componentP z Respectively, respectivelyComprises the following steps:
The field of the horizontal electric dipole source with any orientation has rotation symmetry for the laminated medium, and the derived field can be usedxOn the basis of a calculation formula of the directional field, the field of the horizontal electric dipole source in any direction is directly obtained through coordinate rotation. For a vertical electric dipole source, the distribution of the field does not have rotational symmetry due to the ground boundary, and it needs to be synthesized by the calculation formula of the response of the vertical dipole source derived in step C.
For an azimuth angle ofThe coordinate parameters after rotation of the horizontal electric dipole source are as follows:
E2. Electromagnetic field synthesis
Respectively calculating the electric dipole moments asP e Field vector of response of time-horizontal electric dipole source: (And) And the field vector (E) of the response of the vertical electric dipole source z And H z ) Then, the field components of any electric dipole source in the ground can be synthesized as follows:
E3. And performing inverse Fourier transform on the synthesized electromagnetic field to obtain x, y and z directional components of the electromagnetic response of the dipole source time domain in the horizontal layered earth medium.
In specific implementation, the model of the atmosphere and the stratum in the space is set as shown in FIG. 4, and an azimuth angle is setAngle of inclination of =60 degreesA vector electric dipole source of =30 degrees is located in the ground (0,z s ) And the parameters of each layer of medium are respectively as follows: 1 = 2 = 3 = 0 ;σ 0 =10 –12 S/m;σ 1 =0.05 S/m、σ 2 =0.001 S/m、σ 3 =0.1 And (5) S/m. The coordinate settings are as shown, z-axis up, ground z =0; according to the existence of multiple underground interfaces, each layer is thickh 1 =2000 m,h 2 =5000m,. A vertical electric dipole source is located in layer 2,z s = -4000 m. Respectively obtaining the field components of the vertical electric dipole source and the horizontal electric dipole source by the step C and the step D, then decomposing the electric dipole moment and rotating the coordinate by the step E1, and then obtaining the field vector (C) responded by the horizontal electric dipole source by the step E2And) And the field vector (E) of the response of the vertical electric dipole source z And H z ) The synthesis results in an electromagnetic field response in the frequency domain. Inclination of unit electric dipole source arranged according to this example=30 degrees and azimuth angle=60 degrees each field amount of the vector source is synthesized using equations 37 to 41. The selection frequency of 10Hz is shown in fig. 4-5, for example. Where fig. 5 is taken from a near surface horizontal x-y section and fig. 6 is taken from an x-z longitudinal section (y =0 m). Finally, according to step E3, inverse fourier transform is performed on the synthesized electromagnetic field frequency domain response to obtain the magnitude of the electromagnetic field component in the time domain, which is shown in fig. 7 by taking the electromagnetic field response profile at the time of 0.1s as an example.
It is noted that the disclosed embodiments are intended to aid in further understanding of the invention, but those skilled in the art will appreciate that: various alternatives and modifications are possible without departing from the invention and scope of the appended claims. Therefore, the invention should not be limited to the embodiments disclosed, but the scope of the invention is defined by the appended claims.
Claims (3)
1. A time-frequency electromagnetic response simulation method of a vector dipole source in a horizontal layered geodetic medium is characterized in that a propagation coefficient matrix of a vector potential function caused by a universal electric dipole source and a universal magnetic dipole source in a rectangular coordinate system is established, and electromagnetic field components of the earth surface are obtained through recursion solution; respectively calculating the electromagnetic fields of the horizontal dipole source and the vertical dipole source, and synthesizing to obtain the electromagnetic field response of any vector dipole source; the electric dipole source and the magnetic dipole source can be used for any direction in a horizontal layered ground; calculating after obtaining the amplitude function of the dipole source in the uniform full-space potential function, and obtaining the time-frequency electromagnetic response of the dipole source in the horizontal layered earth medium; the method comprises the following steps:
A. establishing a general expression of the electromagnetic response of a dipole source in a horizontal laminar earth medium:
a = represents a vector potential function of an electromagnetic field generated by an electric dipole source in any direction in a horizontal layerA z ,Is a unit vector of the z-axis direction in a rectangular coordinate system; the vector potential function of the electromagnetic field generated by a magnetic dipole source is denoted F =F z I.e. vector potential functions A and F are constrained to have only z-componentsA z AndF z ;
the components of the electromagnetic field excited by the electric dipole source in any direction in the horizontal layer are respectively expressed as formula 1 and formula 2:
In the formula, complex conductivity=σ+iωε,iRepresents the unit of an imaginary number and is,σin order to be the conductivity of the medium,εis the dielectric constant;E x 、E y 、E z respectively under a rectangular coordinate systemx、y、zAn axial electric field component;H x 、H y 、H z are respectively asx、y、zA magnetic field component of the shaft; whereinH z Is a TM polarized field;
the components of the electromagnetic field caused by the magnetic dipole source in any direction in the horizontal layer are expressed as formula 3 and formula 4, respectively:
Wherein the vertical component of the electric fieldE z The number of the magnetic particles is zero,E z is a TE polarization field;in order to be a resistance ratio of the light, ,irepresents the unit of an imaginary number and is,ω-a signal angular frequency;is medium magnetic permeability; in formulae 1 and 4kSatisfy the requirements ofIs the wave number of the medium;
B. calculating a vector potential function of an electromagnetic field, and respectively solving a solution of a passive region potential function primary field and a solution of an active layer potential function secondary field; the method comprises the following steps:
B1. solving a homogeneous differential equation for the potential function of the passive region to obtain:
In the formula, the upper wavy sign represents a fourier transform form of the corresponding potential function,A + andF + represents the corresponding uplink namely +zThe wave solution of the direction is obtained,A - andF - represents the corresponding downlinkzThe wave solution of the direction is obtained,k x andk y are respectivelyxAndythe number of spatial waves in a direction is,u 2 =k x 2 +k y 2 -k 2 ;
B2. solving a general solution of the potential function of the active layer, comprising:
for any purposeNThe model of the layer of the earth,zin the axial direction; setting the center of the source as the origin of coordinates, the interface where the source is locatedz s =0, the interface being a virtual interface for computing the response of the source; is directly above the source-containing interfacezThe sequence number of the direction is marked as i from bottom to top 1 ,i 2 ,…,i L (ii) a Formation parameter corresponds to conductivityσ ii Dielectric constant ofε ii Thickness of sum layerh ii (ii) a Under the source-containing interface, i.e. negativezThe sequence number of the direction is marked as j from top to bottom 1 ,j 2 ,…,j M The formation parameter corresponds to conductivityσ jj Dielectric constant ofε jj And thickness of each layerh jj (ii) a Is provided withA p AndF p for the amplitude of the source field, the characteristic solution of the non-homogeneous differential equation with source time is expressed as formula 7 and formula 8 according to the plane wave solution of the uniform full space:
nIs composed ofNLayer number in the layer ground;
combining the special solutions of the sources of the formulas 7 and 8 with the solution of the potential function of the passive region in the step B1, obtaining the hankel integral form of the potential function of the active layer, which is expressed as the formulas 9 and 10:
In the formula (I), the compound is shown in the specification,A n + 、A n - 、F n + 、F n - in order to determine the coefficient to be determined,,,coefficient factorβ n For distinguishing active regions fromA passive region;
C. according to the potential function of the underground vertical electric dipole source, calculating to obtain the amplitude of the potential function:
the hankel integral form of the potential function of the vertical electric dipole source in any layer in the transform domain is obtained by the formula 9:
The magnitude of the potential function of a vertical electric dipole source is expressed as:
Wherein the electrode moment at the origin of coordinates isIdl ;Idl The vertical electric dipole source density of (a) is expressed as:
In the formula (I), the compound is shown in the specification,is a unit vector in the z-axis direction; δ (d)x)、δ(y)、δ(z) Is a dirac function;
propagation coefficientA n0 + AndA n0 - solving by a recursion relation; the method comprises the following steps:
C1. determining the top and bottom layer transfer matrices from the source and the medium parameters of each layer to obtain the propagation coefficients of the top and bottom layers, which are expressed as equations 14 and 15:
In the formula (I), the compound is shown in the specification,qfor transferring matricesQThe coefficients of (c);s i as a source vectorS i The coefficient of (a) is calculated,s j as a source vectorS j The coefficient of (a);
and (3) obtaining a transfer matrix expression form of the source layer according to the recursion of the earth lamellar model:
C2. Determining a layer-by-layer propagation coefficient;
setting a virtual layer interface at the position of a source, and dividing the earth model into a source upper part and a source lower part;
U i Is expressed as:
Wherein:
Ith (i) L The layer being the top-bottom interface and extending upwardly into an infinitely extending space, i.e.A Li - If =0 has no downstream wave, then:
Under-source TM polarization modeU j Expressed in matrix form:
The transfer matrix form of the source layer of TE polarization mode is: will be in the formulae from 25 to 30Instead, it is changed into,A p Instead, it is changed intoF p ;
C3. And (3) recursion is carried out layer by layer from the propagation coefficient matrixes of the top layer and the bottom layer to the interior of the model, the propagation coefficient of the source layer is obtained through formulas 16-30, and the propagation coefficient of each layer is obtained through calculationA n + AndA n - ;
substituting the propagation coefficient into the formula 13 to obtain the potential function of each layer;
substituting potential functions into equations 3-6 to obtain kernel functions of components of the electromagnetic field;
then, a digital filtering algorithm is applied to realize the Hankel transformation, and the frequency domain response of each field component is obtained;
D. establishing a potential function of a horizontal electric dipole source in the ground:
the current source density of a horizontal dipole source is expressed as:
In the formula (I), the compound is shown in the specification,is thatxA unit vector of direction of the axis; δ (d)x)、δ(y)、δ(z) Is a dirac function;
the potential function of a horizontal electric dipole source contained in any layer is expressed as:
In the formula (I), the compound is shown in the specification,the symbol represents whenz>When the number 0 is a negative number,z<when 0, the plus number is taken;
the propagation coefficient is expressed as:
Using the propagation coefficient to obtain a potential function, equation 32 to 36 are substituted for 3 to 6 to obtain the kernel functions of all components of the electromagnetic field, and then the Hankel transformation is realized by applying a digital filtering algorithm to obtain any layernFrequency domain response of each field quantity excited by the medium-level electric dipole source;
E. carrying out simulation calculation to obtain a vector electric dipole source field in any direction in the ground; the method comprises the following steps:
E1. decomposing a field source;
electric dipole momentP e =IdlAzimuth angleAnd the angle of inclinationDefining any electric dipole in earth, decomposing electric dipole moment into horizontal componentAnd a vertical componentP z Respectively as follows:
For an azimuth angle ofThe coordinate parameters after rotation of the horizontal electric dipole source are as follows:
E2. Electromagnetic field synthesis;
by an electric dipole moment ofP e Field vector of time-horizontal electric dipole source responseAndand the field vector E of the response of the vertical electric dipole source z And H z And (3) synthesizing to obtain the electromagnetic field components of any electric dipole source in the ground as follows:
E3. Performing inverse Fourier transform on the synthesized electromagnetic field to obtain an electromagnetic field component of a time domain;
therefore, the simulation calculation of the time-frequency electromagnetic response of the dipole source in the horizontal laminar earth medium is realized.
2. The method for simulating the time-frequency electromagnetic response of the vector dipole source in the horizontal laminar geodetic medium according to claim 1, wherein the step B2 represents the Fourier transform form of the vector potential function in the special solution of the equation of the time-inhomogeneous differential with sourcezComponent(s) ofAndproducing attenuation both above and below the source; the horizon containing the source isn=i 1 Andn=j 1 。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211106686.3A CN115186520B (en) | 2022-09-09 | 2022-09-09 | Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202211106686.3A CN115186520B (en) | 2022-09-09 | 2022-09-09 | Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115186520A true CN115186520A (en) | 2022-10-14 |
CN115186520B CN115186520B (en) | 2023-01-24 |
Family
ID=83524639
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202211106686.3A Active CN115186520B (en) | 2022-09-09 | 2022-09-09 | Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115186520B (en) |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140061507A1 (en) * | 2012-08-28 | 2014-03-06 | United States Government, As Represented By The Secretary Of The Navy | Broadband artificial dielectric with dynamically optical control |
CN110376655A (en) * | 2019-07-25 | 2019-10-25 | 中南大学 | The calculation method that dipole source electromagnetic field in any position responds under the conditions of layered geology |
CN112285435A (en) * | 2020-10-29 | 2021-01-29 | 中国舰船研究设计中心 | Equivalent simulation method of high-power magnetic field radiation source |
CN113076508A (en) * | 2021-04-02 | 2021-07-06 | 北京环境特性研究所 | Low-frequency near-field rapid calculation method based on vertical magnetic dipole in half space |
CN114239268A (en) * | 2021-12-16 | 2022-03-25 | 西北工业大学 | Method for acquiring cross-interface radiation field of underwater double-electric-dipole array based on Romberg |
-
2022
- 2022-09-09 CN CN202211106686.3A patent/CN115186520B/en active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20140061507A1 (en) * | 2012-08-28 | 2014-03-06 | United States Government, As Represented By The Secretary Of The Navy | Broadband artificial dielectric with dynamically optical control |
CN110376655A (en) * | 2019-07-25 | 2019-10-25 | 中南大学 | The calculation method that dipole source electromagnetic field in any position responds under the conditions of layered geology |
CN112285435A (en) * | 2020-10-29 | 2021-01-29 | 中国舰船研究设计中心 | Equivalent simulation method of high-power magnetic field radiation source |
CN113076508A (en) * | 2021-04-02 | 2021-07-06 | 北京环境特性研究所 | Low-frequency near-field rapid calculation method based on vertical magnetic dipole in half space |
CN114239268A (en) * | 2021-12-16 | 2022-03-25 | 西北工业大学 | Method for acquiring cross-interface radiation field of underwater double-electric-dipole array based on Romberg |
Non-Patent Citations (2)
Title |
---|
徐建华等: "多层介质中偶极子场的系数递推关系", 《石油地球物理勘探》 * |
翁爱华等: "基于电场不连续边界条件的层状介质电磁格林函数计算", 《吉林大学学报(地球科学版)》 * |
Also Published As
Publication number | Publication date |
---|---|
CN115186520B (en) | 2023-01-24 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Grayver et al. | 3D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation | |
CN110058317B (en) | Aviation transient electromagnetic data and aviation magnetotelluric data joint inversion method | |
Zhang et al. | Magnetotelluric strike rules | |
Revil et al. | Electrical properties of zeolitized volcaniclastic materials | |
CN110058315B (en) | Three-dimensional anisotropic radio frequency magnetotelluric adaptive finite element forward modeling method | |
Lajoie et al. | The electromagnetic response of a conductive inhomogeneity in a layered earth | |
Huang et al. | Bidimensional empirical mode decomposition (BEMD) for extraction of gravity anomalies associated with gold mineralization in the Tongshi gold field, Western Shandong Uplifted Block, Eastern China | |
Pedersen et al. | Groundwater exploration using combined controlled-source and radiomagnetotelluric techniques | |
Liu et al. | 2D sequential inversion of total magnitude and total magnetic anomaly data affected by remanent magnetization | |
Wilson et al. | Real-time 3D inversion of ultra-deep resistivity logging-while-drilling data | |
Saad | Understanding gravity gradients—a tutorial | |
CN111666534B (en) | Electrical random anisotropic electromagnetic field decoupling method | |
CN104992029B (en) | DISCRETE RANDOM MEDIUM modeling method in a kind of multiple dimensioned non-homogeneous lunar soil layer | |
Meju | Simple relative space–time scaling of electrical and electromagnetic depth sounding arrays: implications for electrical static shift removal and joint DC‐TEM data inversion with the most‐squares criterion | |
Yan et al. | Two‐dimensional magnetotelluric inversion using reflection seismic data as constraints and application in the COSC project | |
Papadopoulos et al. | Electrical resistivity tomography for the modelling of cultural deposits and geomophological landscapes at Neolithic sites: a case study from Southeastern Hungary | |
Liu et al. | Effects of electrical anisotropy on long-offset transient electromagnetic data | |
Zaher et al. | Integration of geophysical methods for groundwater exploration: A case study of El Sheikh Marzouq area, Farafra Oasis, Egypt | |
Wang et al. | Efficient finite-volume simulation of the LWD orthogonal azimuth electromagnetic response in a three-dimensional anisotropic formation using potentials on cylindrical meshes | |
CN115186520B (en) | Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium | |
Chen et al. | Induced polarization and magnetic responses of serpentinized ultramafic rocks from mid‐ocean ridges | |
Jol et al. | Stratigraphic imaging of the Navajo Sandstone using ground-penetrating radar | |
CN207673329U (en) | A kind of electric imaging logging image enhancement display processing system | |
Aghil et al. | Delineation of electrical resistivity structure using Magnetotellurics: a case study from Dholera coastal region, Gujarat, India | |
EP1070970B1 (en) | A method of three dimensional reconstructing a physical magnitude inside a borehole |
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 |