CN115186520B - 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 PDF

Info

Publication number
CN115186520B
CN115186520B CN202211106686.3A CN202211106686A CN115186520B CN 115186520 B CN115186520 B CN 115186520B CN 202211106686 A CN202211106686 A CN 202211106686A CN 115186520 B CN115186520 B CN 115186520B
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.)
Active
Application number
CN202211106686.3A
Other languages
Chinese (zh)
Other versions
CN115186520A (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.)
Peking University
Original Assignee
Peking 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 Peking University filed Critical Peking University
Priority to CN202211106686.3A priority Critical patent/CN115186520B/en
Publication of CN115186520A publication Critical patent/CN115186520A/en
Application granted granted Critical
Publication of CN115186520B publication Critical patent/CN115186520B/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
    • 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 comprises the steps of establishing 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, and recursively solving to obtain an electromagnetic field component of the earth surface; 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 precursor in the laminar earth, underground nuclear pollution source or electromagnetic wave field radiated by electric or magnetic vector dipole source of coal mine tunnel, and has good applicability and high application value.

Description

Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium
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 a vector dipole source in a horizontal layered earth medium, which simulates through a potential function propagation coefficient matrix and solves the electromagnetic field response of a time domain and a frequency domain 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 adaptive to any direction exists in the current work for calculating the electromagnetic field response of the dipole source in the laminated medium. 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 method for simulating time-frequency electromagnetic response of a dipole source in a horizontal layered earth medium, which respectively calculates 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 explaining the electromagnetic disturbance phenomenon generated by the vector dipole source in seismic precursors in a layered earth, and further can be applied to large earthquake prediction and identification of electromagnetic radiation generated in underground nuclear pollution sources or coal mine tunnels.
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 layered earth 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 for the electromagnetic response of a dipole source in a horizontal laminar earth medium:
let A =denote the vector potential function of the electromagnetic field generated by the electric dipole source in any direction in the horizontal laminar groundA z
Figure DEST_PATH_IMAGE001
，
Figure 252649DEST_PATH_IMAGE001
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
Figure 802580DEST_PATH_IMAGE001
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:
Figure 385745DEST_PATH_IMAGE002
formula 1
Figure 679323DEST_PATH_IMAGE003
Formula 2
In the formula, the complex conductivity
Figure 267431DEST_PATH_IMAGE004
=σ+iωε，iRepresents the unit of an imaginary number and is,σin order to be the conductivity of the medium,εis a dielectric constant；E x 、E y 、E z Respectively under rectangular coordinate systemx、y、zAn electric field component in the axial direction;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 type dipole source in any direction in the horizontal layer ground are expressed as formula 3 and formula 4, respectively:
Figure 988262DEST_PATH_IMAGE005
formula 3
Figure 419243DEST_PATH_IMAGE006
Formula 4
Wherein the vertical component of the electric fieldE z The number of the magnetic particles is zero,E z is a TE polarization field;
Figure 391879DEST_PATH_IMAGE007
in order to be a resistance ratio of the light,
Figure 693547DEST_PATH_IMAGE008
，irepresents the unit of an imaginary number,ω-a signal angular frequency;
Figure 726225DEST_PATH_IMAGE009
is medium magnetic permeability; in formulas 1 and 4kSatisfy the requirements of
Figure 644502DEST_PATH_IMAGE010
Is the wave number of the medium;
B. calculating a vector potential function of the electromagnetic field, and respectively solving a solution of a primary field of the passive region potential function and a solution of a secondary field of the active layer potential function;
the method is considered in two cases, namely a solution of a passive region primary electromagnetic field vector potential function and a solution of an active region secondary electromagnetic field vector potential function;
B1. for the potential function of the passive region, a homogeneous differential equation can be solved to obtain:
Figure 14304DEST_PATH_IMAGE011
formula 5
Figure 75539DEST_PATH_IMAGE012
Formula 6
In the formula, the upper wavy sign represents a fourier transform form of the corresponding potential function,A + andF + represents the corresponding uplink (+)zDirection) of the wave solution,A - andF - then represents the corresponding downlink-zDirection) wave solution is performed,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 a 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 numbers of the directions are marked as i from bottom to top in sequence 1 ，i 2 ，…，i L (ii) a Formation parameter corresponds to conductivityσ ii Dielectric constant ofε ii Thickness of sum layerh ii (ii) a Is immediately negative below the source-containing interfacezThe 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 to the amplitude of the source field, according to a uniform sumThe spatial plane wave solution, the solution of the source-containing time-inhomogeneous differential equation, is expressed by the following equations 7 and 8:
Figure 138173DEST_PATH_IMAGE013
formula 7
Figure 153533DEST_PATH_IMAGE014
Formula 8
In the form of a Fourier transform of a vector potential functionzComponent(s) of
Figure 123763DEST_PATH_IMAGE016
And
Figure 540969DEST_PATH_IMAGE017
producing attenuation both above and below the source;Nin the course of the earth, the water is discharged,nis a layer number, the source-containing layer isn=i 1 Andn=j 1 ；
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 B1 to obtain the Hankel integral form of the potential function of the active layer, which is expressed as the formulas 9 and 10:
Figure 774504DEST_PATH_IMAGE018
formula 9
Figure 277161DEST_PATH_IMAGE019
Formula 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,
Figure 254344DEST_PATH_IMAGE020
，
Figure 296030DEST_PATH_IMAGE021
，
Figure 231625DEST_PATH_IMAGE022
by coefficient factorβ n To distinguish between active and inactive regions,β n the values are as follows:
Figure 221578DEST_PATH_IMAGE023
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
Figure 736873DEST_PATH_IMAGE024
The vertical electric dipole source density of (a) is expressed as:
Figure 987726DEST_PATH_IMAGE025
formula 11
In the formula (I), the compound is shown in the specification,
Figure 704009DEST_PATH_IMAGE026
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:
Figure 40313DEST_PATH_IMAGE027
formula 12
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:
Figure 359298DEST_PATH_IMAGE028
formula 13
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:
Figure 340024DEST_PATH_IMAGE029
formula 14
Figure 351842DEST_PATH_IMAGE030
Formula 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 medium parameter;s i as a source vectorS i 2 coefficients of
Figure 909863DEST_PATH_IMAGE031
，s j As a source vectorS j 2 coefficients of
Figure 406441DEST_PATH_IMAGE032
Coefficient values are determined by the source;
and (3) obtaining a transfer matrix expression form of the source layer according to the recursion of the earth lamellar model:
Figure 631886DEST_PATH_IMAGE033
formula 16
Figure 549026DEST_PATH_IMAGE034
Formula 17
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 above the source and below the source;
binding 16 to 17 over-source TM polarization mode
Figure 867DEST_PATH_IMAGE035
The propagation coefficient is written as:
Figure 927235DEST_PATH_IMAGE036
formula 18
U i Is expressed as:
Figure 616973DEST_PATH_IMAGE037
formula 19
Wherein:
Figure 236174DEST_PATH_IMAGE038
formula 20
Figure 768786DEST_PATH_IMAGE039
Formula 21
Figure 843052DEST_PATH_IMAGE040
Formula 22
Figure 777510DEST_PATH_IMAGE041
Formula 23
Ith 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:
Figure 302033DEST_PATH_IMAGE042
formula 24
TM polarization mode under sourceU j Expressed in matrix form:
Figure 695843DEST_PATH_IMAGE043
formula 25
Figure 964013DEST_PATH_IMAGE044
Formula 26
Figure 752977DEST_PATH_IMAGE045
Formula 27
Figure 182822DEST_PATH_IMAGE046
Formula 28
Figure 565393DEST_PATH_IMAGE047
Formula 29
Figure 637254DEST_PATH_IMAGE048
Formula 30
Similarly, the source layers of the TE polarization mode have the same transfer matrix form, in which
Figure 280725DEST_PATH_IMAGE049
Is composed of
Figure 22416DEST_PATH_IMAGE050
，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:
Figure 751337DEST_PATH_IMAGE051
formula 31
In the formula (I), the compound is shown in the specification,
Figure 361310DEST_PATH_IMAGE052
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:
Figure 265812DEST_PATH_IMAGE053
formula 32
Figure 37459DEST_PATH_IMAGE054
Formula 33
In the formula (I), the compound is shown in the specification,
Figure 519256DEST_PATH_IMAGE055
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:
Figure 306821DEST_PATH_IMAGE056
formula 34
Figure 924884DEST_PATH_IMAGE057
Formula 35
Figure 867433DEST_PATH_IMAGE058
Formula 36
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 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 angle
Figure 711892DEST_PATH_IMAGE059
And the angle of inclination
Figure 663667DEST_PATH_IMAGE060
Three parameters define any electric dipole in the ground, and the electric dipole moment is decomposed into horizontal components
Figure 136237DEST_PATH_IMAGE061
And a vertical componentP z Respectively is as follows:
Figure 656211DEST_PATH_IMAGE062
formula 37
Figure 847021DEST_PATH_IMAGE063
Formula 38
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 of
Figure 336908DEST_PATH_IMAGE059
The coordinate parameters after rotation of the horizontal electric dipole source are as follows:
Figure 929564DEST_PATH_IMAGE064
formula 39
E2. Electromagnetic field synthesis;
an electric dipole moment ofP e Field vector of time-horizontal electric dipole source response: (
Figure 620439DEST_PATH_IMAGE065
And
Figure 32966DEST_PATH_IMAGE066
) 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:
Figure 983603DEST_PATH_IMAGE067
formula 40
Figure 430765DEST_PATH_IMAGE068
Formula 41
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 the calculation technology of the electromagnetic response of the line source and the surface source is realized by utilizing multi-source combination. The invention can be applied to simulating the electromagnetic disturbance of the earthquake precursor in the layered earth, the electromagnetic wave field radiated by the electric or magnetic vector dipole source of the underground nuclear pollution source or the coal mine tunnel, can be used for explaining the electromagnetic disturbance phenomenon generated by the vector dipole source in the earthquake precursor in the layered earth, is further used for predicting and identifying the electromagnetic radiation generated in the underground nuclear pollution source or the coal mine tunnel in the great earthquake, and has 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,
Figure 417175DEST_PATH_IMAGE059
and
Figure 723523DEST_PATH_IMAGE060
respectively representing the azimuth and inclination of the vector dipole source,P e dipole source pitch, resolvable into horizontal components
Figure 555212DEST_PATH_IMAGE061
And a vertical 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,
Figure 997826DEST_PATH_IMAGE069
which is representative of the 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 vector electric dipole source in the ground 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 an x-z (vertical) profile of a vector electric dipole source in the ground 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. 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 the 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 the example (x s ，y s ，z s ) And sets 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 generic 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 the electric dipole source be A =A z
Figure 155138DEST_PATH_IMAGE070
，
Figure 683203DEST_PATH_IMAGE071
In a rectangular coordinate systemzThe unit vector in the direction of the axis,A z is the z component of A
Figure 318583DEST_PATH_IMAGE072
(ii) a And the vector potential function of the electromagnetic field generated by the magnetic dipole source is denoted as F =F z
Figure 5917DEST_PATH_IMAGE026
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:
Figure 708031DEST_PATH_IMAGE073
formula 1
Figure 848025DEST_PATH_IMAGE074
Formula 2
In the formula, complex conductivity
Figure 896884DEST_PATH_IMAGE075
=σ+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 electric field component 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 known as the TM polarization field. Similarly, the components of the electromagnetic field caused by a magnetic dipole source are respectively represented as:
Figure 438724DEST_PATH_IMAGE076
formula 3
Figure 813204DEST_PATH_IMAGE077
Formula 4
Wherein, due toE z (vertical component of electric field) is zero,E z known as the TE polarization field.
Figure 174916DEST_PATH_IMAGE078
In order to be a resistance ratio of the resistor,
Figure 152099DEST_PATH_IMAGE079
，irepresents an imaginary unit;ω-a signal angular frequency;
Figure 689391DEST_PATH_IMAGE080
is medium magnetic permeability; for general applications, take
Figure 93827DEST_PATH_IMAGE069
=
Figure 677255DEST_PATH_IMAGE069
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 the vector potential function of the primary electromagnetic field of the passive region (complementary solution of homogeneous partial differential equation) and the solution of the vector potential function of the secondary electromagnetic field of the active-containing layer (special solution of non-homogeneous partial differential equation);
B1. solving a homogeneous differential equation for the potential function of the passive region to obtain:
Figure 458129DEST_PATH_IMAGE011
formula 5
Figure 348463DEST_PATH_IMAGE012
Formula 6
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 laminar earth medium, i.e. except i in FIG. 2 1 Or j 1 The outer regions can be solved in different inner regions, and the amplitude (constant) of each layer can be deduced by combining with proper boundary conditions.
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 source-containing interface (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 And a layerThickness ofh 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:
Figure 923801DEST_PATH_IMAGE013
formula 7
Figure 994525DEST_PATH_IMAGE014
Formula 8
In the form of a Fourier transform of the vector potential functionzComponent(s) of
Figure 188877DEST_PATH_IMAGE081
And
Figure 559815DEST_PATH_IMAGE082
producing 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 source-containing horizon 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:
Figure 306054DEST_PATH_IMAGE018
formula 9
Figure 5020DEST_PATH_IMAGE019
Formula (II)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,
Figure 127697DEST_PATH_IMAGE020
，
Figure 353142DEST_PATH_IMAGE021
，
Figure 145649DEST_PATH_IMAGE022
here by coefficient factorsβ n Distinguishing between active and passive conditions, coefficientβ n The values are:
Figure 190965DEST_PATH_IMAGE083
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 an 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 the coordinates, the electrode moment isIdl
Figure 851753DEST_PATH_IMAGE084
Perpendicular electric dipole source density (A/m) 2 ) Can be expressed as:
Figure 571186DEST_PATH_IMAGE085
formula 11
In the formula (I), the compound is shown in the specification,
Figure 659227DEST_PATH_IMAGE024
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. Potential of vertical electric dipole sourceThe function magnitude is represented as:
Figure 191840DEST_PATH_IMAGE086
formula 12
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:
Figure 797265DEST_PATH_IMAGE087
formula 13
Propagation coefficientA n + AndA n - the result is obtained from the recursion relation,β n see step B2. The recurrence relation includes:
C1. determining a top layer transfer matrix and a bottom layer transfer matrix according to the source and each layer medium parameter to obtain the propagation coefficients of the top layer and the bottom layer:
Figure 997302DEST_PATH_IMAGE088
formula 14
Figure 256245DEST_PATH_IMAGE089
Formula 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 medium parameter;s i as a source vectorS i Two coefficients of
Figure 151520DEST_PATH_IMAGE090
，s j As a source vectorS j Two coefficients of
Figure 685269DEST_PATH_IMAGE091
The 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:
Figure 474234DEST_PATH_IMAGE033
formula 16
Figure 779444DEST_PATH_IMAGE034
Formula 17
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 above-source TM polarization mode is represented in matrix form:
Figure 286649DEST_PATH_IMAGE036
formula 18
Figure 358510DEST_PATH_IMAGE037
Formula 19
Figure 381742DEST_PATH_IMAGE092
Formula 20
Figure 513646DEST_PATH_IMAGE039
Formula 21
Figure 242567DEST_PATH_IMAGE093
Formula 22
Figure 993486DEST_PATH_IMAGE041
Formula 23
Ith L The layer being the top-bottom interface and extending upwardly into an infinitely extending space, i.e.A Li - =0 no downgoing wave, then:
Figure 491463DEST_PATH_IMAGE042
formula 24
Under-source TM polarization modeU j Expressed in matrix form:
Figure 528689DEST_PATH_IMAGE043
formula 25
Figure 885852DEST_PATH_IMAGE044
Formula 26
Figure 33937DEST_PATH_IMAGE094
Formula 27
Figure 652000DEST_PATH_IMAGE046
Formula 28
Figure 735494DEST_PATH_IMAGE095
Formula 29
Figure 704587DEST_PATH_IMAGE096
Formula 30
Similarly, the source layers of the TE polarization mode have the same transfer matrix form, and only need to be replaced
Figure 656362DEST_PATH_IMAGE097
Is composed of
Figure 502833DEST_PATH_IMAGE098
，A p Is composed ofF p 。
C3. And (3) recursion is carried out layer by layer from the propagation coefficient matrixes of the top layer (formula 14) and the bottom layer (formula 15) to the interior of the model, the propagation coefficient of the source layer is obtained through formulas 16 to 30, and then each layer is obtained through calculationnPropagation 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:
Figure 881862DEST_PATH_IMAGE099
formula 31
In the formula (I), the compound is shown in the specification,
Figure 72672DEST_PATH_IMAGE052
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:
Figure 703505DEST_PATH_IMAGE100
formula 32
Figure 296160DEST_PATH_IMAGE101
Formula 33
In the formula (I), the compound is shown in the specification,
Figure 846090DEST_PATH_IMAGE102
the symbol represents whenz>When the number 0 is a negative number,z<when 0, the + number is taken. The propagation coefficients (levels) are:
Figure 399562DEST_PATH_IMAGE103
formula 34
Figure 693140DEST_PATH_IMAGE104
Formula 35
Figure 140302DEST_PATH_IMAGE105
Formula 36
Processing similar to the step C3, calculating a potential function by using the propagation coefficient, calculating each component kernel function of the electromagnetic field in the steps of 3 to 6, and finally realizing Hankel transformation by applying a digital filtering algorithm to obtain any layernAnd the frequency domain response of each field quantity 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
Figure 861133DEST_PATH_IMAGE059
And the angle of inclination
Figure 636323DEST_PATH_IMAGE060
Three parameters define an arbitrary electric dipole, as shown in FIG. 3, which can be decomposed into horizontal components
Figure 733592DEST_PATH_IMAGE061
And a vertical componentP z Respectively is as follows:
Figure 300839DEST_PATH_IMAGE106
formula 37
Figure 566473DEST_PATH_IMAGE063
Formula 38
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 field distribution does not have rotational symmetry due to the ground boundary, and it needs to be synthesized by the calculation formula of the vertical dipole source response derived in step C.
For an azimuth angle of
Figure 953592DEST_PATH_IMAGE059
The coordinate parameters after rotation of the horizontal electric dipole source are as follows:
Figure 854552DEST_PATH_IMAGE064
formula 39
E2. Electromagnetic field synthesis
Respectively calculating the electric dipole moments asP e Field vector of response of time-horizontal electric dipole source: (
Figure 151672DEST_PATH_IMAGE065
And
Figure 214306DEST_PATH_IMAGE066
) And field vector of vertical electric dipole source response (E) z And H z ) Then, the field component of any electric dipole source in the earth can be obtained by synthesis as follows:
Figure 88721DEST_PATH_IMAGE067
formula 40
Figure 403159DEST_PATH_IMAGE068
Formula 41
E3. And performing inverse Fourier transform on the synthesized electromagnetic field to obtain x, y and z directional components of electromagnetic response of a dipole source time domain in the horizontal laminar earth medium.
In specific implementation, the model of the space atmosphere and the stratum is set as shown in FIG. 4, and an azimuth angle is set
Figure 679420DEST_PATH_IMAGE107
Angle of inclination of =60 degrees
Figure 912955DEST_PATH_IMAGE108
A 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:
Figure 415612DEST_PATH_IMAGE009
1 =
Figure 392795DEST_PATH_IMAGE009
2 =
Figure 523562DEST_PATH_IMAGE009
3 =
Figure 567479DEST_PATH_IMAGE009
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; depending on the presence of multiple subsurface layers each having a thicknessh 1 =2000 m，h 2 =5000m，
Figure 416486DEST_PATH_IMAGE109
. A vertical electric dipole source is located in layer 2,z s and = 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, decomposing the electric dipole moment and rotating the coordinate by the step E1, and obtaining the horizontal electric dipole source by the step E2Field vector of response: (
Figure 197361DEST_PATH_IMAGE065
And
Figure 182634DEST_PATH_IMAGE066
) And field vector of vertical electric dipole source response (E) 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
Figure 898917DEST_PATH_IMAGE060
=30 degrees and azimuth angle
Figure 235221DEST_PATH_IMAGE059
=60 degrees each field amount of the vector source is synthesized using equations 37 to 41. The selection frequency of 10Hz is shown in FIGS. 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 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; 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 laminar geodetic medium; the method comprises the following steps:
A. establishing a general expression of the dipole source electromagnetic response in the horizontal laminar earth medium:
a = is recorded as a vector potential function of an electromagnetic field generated by an electric dipole source in any direction in a horizontal laminar earth mediumA z
Figure 633535DEST_PATH_IMAGE002
，
Figure 979197DEST_PATH_IMAGE003
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
Figure 381360DEST_PATH_IMAGE004
I.e. defining vector potential functions A and F 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 laminar earth medium are respectively expressed as formula 1 and formula 2:
Figure DEST_PATH_IMAGE005
formula 1
Figure 112555DEST_PATH_IMAGE006
Formula 2
In the formula, complex conductivity
Figure DEST_PATH_IMAGE007
=σ+iωε，iRepresents the unit of an imaginary number,σ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 laminar earth medium are expressed as formula 3 and formula 4, respectively:
Figure 119301DEST_PATH_IMAGE008
formula 3
Figure DEST_PATH_IMAGE009
Formula 4
Wherein the vertical component of the electric fieldE z The number of the carbon atoms is zero,E z is a TE polarization field;
Figure 686680DEST_PATH_IMAGE010
in order to be a resistance ratio of the resistor,
Figure DEST_PATH_IMAGE011
，irepresents the unit of an imaginary number,ωis the signal angular frequency;
Figure 486008DEST_PATH_IMAGE012
is medium magnetic permeability; in formulae 1 and 4kSatisfy the requirement of
Figure 743814DEST_PATH_IMAGE013
Is the wave number of the medium;
B. calculating a vector potential function of the 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:
Figure DEST_PATH_IMAGE014
formula 5
Figure 986708DEST_PATH_IMAGE015
Formula 6
In the formula, the upper wavy sign represents a fourier transform form of the corresponding potential function,A + andF + represents the corresponding uplink, i.e. positivezThe wave solution of the direction is determined,A - andF - representing a corresponding downstream, i.e. negativezThe 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 is a virtual interface for computing the response of the source; is directly above the source-containing interfacezThe sequence numbers of the directions are marked as i from bottom to top in sequence 1 ，i 2 ，…，i L (ii) a Formation parameter corresponds to conductivityσ ii Dielectric constant ofε ii Thickness of sum layerh ii (ii) a Is immediately negative below the source-containing interfacezThe 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:
Figure DEST_PATH_IMAGE016
formula 7
Figure 556230DEST_PATH_IMAGE017
Formula 8
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:
Figure DEST_PATH_IMAGE018
formula 9
Figure 907052DEST_PATH_IMAGE019
Formula 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,
Figure DEST_PATH_IMAGE020
，
Figure 878419DEST_PATH_IMAGE021
，
Figure DEST_PATH_IMAGE022
coefficient factor ofβ n For distinguishing active and inactive regions;
C. according to the potential function of a vertical electric dipole source in the ground, calculating to obtain the amplitude of the potential function:
the hankel integral form of the potential function containing the vertical electric dipole source in any layer in the transform domain is obtained by the formula 9:
Figure 416848DEST_PATH_IMAGE023
formula 13
The magnitude of the potential function of a vertical electric dipole source is expressed as:
Figure DEST_PATH_IMAGE024
formula 12
Wherein the electrode moment at the origin of coordinates isIdl
Figure DEST_PATH_IMAGE026
；Idl
Figure 817873DEST_PATH_IMAGE026
The vertical electric dipole source density of (a) is expressed as:
Figure 365529DEST_PATH_IMAGE027
formula 11
In the formula (I), the compound is shown in the specification,
Figure DEST_PATH_IMAGE028
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:
Figure 942135DEST_PATH_IMAGE029
formula 14
Figure DEST_PATH_IMAGE030
Formula 15
In the formula (I), the compound is shown in the specification,qfor transferring matricesQThe coefficient of (a);s i as a source vectorS i The coefficient of (a) is determined,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:
Figure 776099DEST_PATH_IMAGE031
formula 16
Figure DEST_PATH_IMAGE032
Formula 17
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;
over-source TM polarization mode
Figure 805366DEST_PATH_IMAGE033
The propagation coefficient is expressed as:
Figure DEST_PATH_IMAGE034
formula 18
U i Is expressed as:
Figure 891134DEST_PATH_IMAGE035
formula 19
Wherein:
Figure 837093DEST_PATH_IMAGE036
formula 20
Figure 514062DEST_PATH_IMAGE037
Formula 21
Figure DEST_PATH_IMAGE038
Formula 22
Figure 785554DEST_PATH_IMAGE039
Formula 23
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:
Figure DEST_PATH_IMAGE040
formula 24
Under-source TM polarization modeU j Expressed in matrix form:
Figure 534067DEST_PATH_IMAGE041
formula 25
Figure 741057DEST_PATH_IMAGE042
Formula 26
Figure DEST_PATH_IMAGE043
Formula 27
Figure 933135DEST_PATH_IMAGE044
Formula 28
Figure 327208DEST_PATH_IMAGE045
Formula 29
Figure DEST_PATH_IMAGE046
Formula 30
The transfer matrix form of the source layer of TE polarization mode is: will be represented by the formulae 25 to 30
Figure 879412DEST_PATH_IMAGE047
Instead, it is changed into
Figure DEST_PATH_IMAGE048
，A p Instead, it is changed intoF 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 formulas 16-30, and calculating the propagation coefficient of each layerA n + AndA n - ；
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 the horizontal dipole source is expressed as:
Figure 19537DEST_PATH_IMAGE049
formula 31
In the formula (I), the compound is shown in the specification,
Figure DEST_PATH_IMAGE050
is thatxA unit vector of direction of the axis; δ (d)x)、δ(y)、δ(z) Is a dirac function;
the potential function of the horizontal dipole source contained in any layer is expressed as:
Figure 631784DEST_PATH_IMAGE051
formula 32
Figure DEST_PATH_IMAGE052
Formula 33
In the formula (I), the compound is shown in the specification,
Figure 247573DEST_PATH_IMAGE053
the symbol is expressed asz>When the number 0 is a negative number,z<when 0, the plus number is taken;
the propagation coefficient is expressed as:
Figure DEST_PATH_IMAGE054
formula 34
Figure 354201DEST_PATH_IMAGE055
Formula 35
Figure DEST_PATH_IMAGE056
Formula 36
Using propagation coefficients to obtain potential functions, substituting formulas 3 to 6 with formulas 32 to 36 to obtain kernel functions of each component of the electromagnetic field, and then applying a digital filtering algorithm to realize Hankel transformation to obtain any layernFrequency domain response of each field quantity excited by the middle horizontal 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 angle
Figure 863679DEST_PATH_IMAGE057
And the angle of inclination
Figure DEST_PATH_IMAGE058
Defining any electric dipole in earth, decomposing electric dipole moment into horizontal component
Figure 129051DEST_PATH_IMAGE059
And a vertical componentP z Respectively is as follows:
Figure DEST_PATH_IMAGE060
formula 37
Figure 232136DEST_PATH_IMAGE061
Formula 38
For an azimuth angle of
Figure 657301DEST_PATH_IMAGE057
The rotated coordinate parameters of the horizontal electric dipole source are as follows:
Figure DEST_PATH_IMAGE062
formula 39
E2. Electromagnetic field synthesis;
by an electric dipole moment ofP e Field vector of time-horizontal electric dipole source response
Figure 162232DEST_PATH_IMAGE063
And
Figure DEST_PATH_IMAGE064
and the field vector E of the vertical electric dipole source response z And H z And (3) synthesizing to obtain the electromagnetic field component of any electric dipole source in the ground as follows:
Figure 335856DEST_PATH_IMAGE065
formula 40
Figure DEST_PATH_IMAGE066
Formula 41
E3. Performing inverse Fourier transform on the synthesized electromagnetic field to obtain an electromagnetic field component of a time domain;
therefore, analog calculation of 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 a vector dipole source in a horizontal laminar geodetic medium as claimed in claim 1, wherein the step B2 represents the Fourier transform form of a vector potential function in the special solution of a time-inhomogeneous differential equation with a sourcezComponent(s) of
Figure 50871DEST_PATH_IMAGE067
And
Figure DEST_PATH_IMAGE068
producing attenuation both above and below the source; the horizon containing the source isn=i 1 Andn=j 1 。
3. the method for simulating the time-frequency electromagnetic response of the vector dipole source in the horizontally-laminated geodetic medium according to claim 1, wherein in step B2,β n the values are as follows:
Figure 499301DEST_PATH_IMAGE069
。
CN202211106686.3A 2022-09-09 2022-09-09 Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium Active CN115186520B (en)

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 CN115186520A (en) 2022-10-14
CN115186520B true 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 (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239268A (en) * 2021-12-16 2022-03-25 西北工业大学 Method for acquiring cross-interface radiation field of underwater double-electric-dipole array based on Romberg

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8784704B2 (en) * 2012-08-28 2014-07-22 The United States Of America 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
CN112285435B (en) * 2020-10-29 2022-04-22 中国舰船研究设计中心 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

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239268A (en) * 2021-12-16 2022-03-25 西北工业大学 Method for acquiring cross-interface radiation field of underwater double-electric-dipole array based on Romberg

Also Published As

Publication number Publication date
CN115186520A (en) 2022-10-14

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
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
CN110058315B (en) Three-dimensional anisotropic radio frequency magnetotelluric adaptive finite element forward modeling method
Teasdale Methods for understanding poorly exposed terranes: the interpretive geology and tectonothermal evolution of the western Gawler Craton
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
Papadopoulos et al. Electrical resistivity tomography for the modelling of cultural deposits and geomophological landscapes at Neolithic sites: a case study from Southeastern Hungary
Zhang et al. Regional gravity survey and application in oil and gas exploration in China
Lee et al. Three-dimensional facies architecture and three-dimensional calcite concretion distributions in a tide-influenced delta front, Wall Creek Member, Frontier Formation, Wyoming
Ringrose et al. Permeability estimation functions based on forward modeling of sedimentary heterogeneity
Tu et al. Joint focusing inversion of marine controlled-source electromagnetic and full tensor gravity gradiometry data
Yang et al. Mantle flow in the vicinity of the eastern edge of the Pacific‐Yakutat slab: Constraints from shear wave splitting analyses
CN115186520B (en) Time-frequency electromagnetic response simulation method of vector dipole source in horizontal layered geodetic medium
Ramdani et al. How in-place volumes of subsurface reservoir models are impacted by using 3D high-resolution outcrop analogue data. A case study using depositional architectural heterogeneity of stromatoporoid/coral buildups of the Hanifa Fm, Saudi Arabia
Jol et al. Stratigraphic imaging of the Navajo Sandstone using ground-penetrating radar
Aghil et al. Delineation of electrical resistivity structure using Magnetotellurics: a case study from Dholera coastal region, Gujarat, India
Chen et al. Induced polarization and magnetic responses of serpentinized ultramafic rocks from mid‐ocean ridges
Malik et al. Ground‐Penetrating Radar Investigations along Hajipur Fault: Himalayan Frontal Thrust—Attempt to Identify Near Subsurface Displacement, NW Himalaya, India
Thiel et al. Full 3D reservoir mapping using deep directional resistivity measurements
Hoversten et al. Subsalt imaging via seaborne electromagnetics
CN113267830A (en) Two-dimensional gravity gradient and seismic data joint inversion method based on non-structural grid
Thiel et al. Fast 2D inversion for enhanced real-time reservoir mapping while drilling

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