Time domain differential pulse transient electromagnetic multi-resolution migration imaging method

Abstract

The invention discloses a time-domain differential pulse transient electromagnetic multi-resolution migration imaging method, which belongs to the technical field of electromagnetic detection, and comprises the following steps of: S1, exciting an electromagnetic field by adopting time-domain pseudo-random pulses with different pulse widths during a power supply period to obtain a multi-resolution transient electromagnetic diffusion field; s2, extracting a virtual wave field from the transient electromagnetic diffusion field by adopting preprocessing fine integral wave field inverse transformation; and S3, using a finite difference migration imaging technology to obtain imaging profiles of different pulse width time domain pulses, and obtaining high-resolution characteristics of different depth electrical interfaces through a superposition technology. According to the time domain differential pulse transient electromagnetic multi-resolution migration imaging method, high-resolution characteristics of an electrical interface from shallow to deep can be obtained, and the problem of multi-resolution detection is well solved.

Description

Technical Field

[0001] The invention relates to the technical field of electromagnetic detection, and in particular to a time domain differential pulse transient electromagnetic multi-resolution migration imaging method. Background Art

[0002] As an important geophysical exploration method, the electrical source transient electromagnetic method is widely used in metal mineral exploration, coalfield geological disaster investigation, hydrogeophysical exploration and tunnel advance prediction due to its advantages of high efficiency and deep detection. However, the transient electromagnetic field excited by traditional square waves has insufficient resolution and is difficult to achieve multi-resolution detection effects; in addition, the volume effect of the transient electromagnetic diffusion field is more obvious, and the migration imaging and inversion methods based on the diffusion equation have poor resolution for electrical interfaces.

[0003] Traditional field sources use square waves with fixed fundamental frequencies to excite transient electromagnetic fields, and their harmonic energy is mainly concentrated in the low-frequency band, resulting in limited detection resolution. In order to improve the detection capability of transient electromagnetic methods, some scholars have used the ON-Time emission waveform as a research entry point. According to relevant research, it is concluded that traditional square waves lack high-frequency harmonic components, which makes it difficult to fundamentally improve the detection resolution of time-domain electromagnetic methods. Relevant research on urban underground space imaging methods points out that by emitting time-domain pseudo-random pulses of different pulse widths during the ON-Time period for scanning detection, multi-resolution transient electromagnetic diffusion fields can be excited, thereby improving detection resolution.

[0004] Although some scholars have introduced multi-window scanning technology based on the pre-processing fine integration algorithm, which effectively improves the noise resistance of the virtual wave field, the event axis of the virtual wave field can usually only reflect the characteristics of the shallow strata well. As the underground structure becomes more and more complex, the event axis of the deep space virtual wave field cannot accurately represent the stratum morphology, and it is necessary to use the offset imaging method to spatially reposition it. Although subsequent scholars have proposed corresponding imaging methods, and all of them can obtain the electrical interface characteristics of the strata, it is difficult to achieve multi-resolution imaging from shallow to deep due to the limitation of the field source resolution. Summary of the invention

[0005] The purpose of the present invention is to provide a time-domain differential pulse transient electromagnetic multi-resolution migration imaging method, which can obtain high-resolution characteristics of electrical interfaces from shallow to deep, and better solve the multi-resolution detection problem.

[0006] To achieve the above object, the present invention provides a time domain differential pulse transient electromagnetic multi-resolution migration imaging method, comprising the following steps:

[0007] Step S1, during power supply, using time domain pseudo-random pulses with different pulse widths to excite the electromagnetic field to obtain a multi-resolution transient electromagnetic diffusion field;

[0008] Step S2, extracting a virtual wave field from the transient electromagnetic diffusion field by using a pre-processed fine integral wave field inverse transformation;

[0009] Step S3: using finite difference migration imaging technology to obtain imaging profiles of time-domain pulses with different pulse widths, and obtaining high-resolution features of electrical interfaces at different depths through superposition technology.

[0010] Preferably, step S1 specifically comprises:

[0011] The second-order backward difference algorithm BDF2 and the displacement inverse-Krylov subspace algorithm are used during power supply and after the current is turned off, respectively. For a given calculation area Ω, the Maxwell equations ignoring the displacement current are shown as follows:

[0012]

[0013] Where, e is the electric field intensity; b is the magnetic induction intensity; μ0 is the vacuum magnetic permeability; J s is the current density; σ is the conductivity of the formation; represents the differential operator; t represents the diffusion time;

[0014] The space is discretized using an unstructured tetrahedral grid, the magnetic induction intensity b is eliminated from formulas (1a) and (1b), and the Galerkin method is used to obtain the governing equation of the electric field, as shown below:

[0015]

[0016] Among them, N i and N j is the interpolation basis function within the unit; V represents the volume of the three-dimensional space; dV represents the volume integral; the derivative term in formula (2) is obtained by the second-order backward difference method:

[0017]

[0018] Among them, e n and represent the electric field intensity and current density at the nth moment respectively.

[0019] Preferably, let the matrix M σ =∫ V μ0σN i ·N j dV,M inner =∫ V μ0N i ·N j dV, the electric field strength during power supply is solved by the following formula:

[0020]

[0021] After the current is turned off, J in formula (2) s is zero, rewrite it as:

[0022]

[0023] The displacement inverse-Krylov subspace algorithm is used to solve formula (5), and the approximate expression of the solution is:

[0024]

[0025] Where m is the order of the displacement inverse-Krylov subspace; γ is the displacement parameter; I is the unit matrix; e0 is the electric field intensity after the current is turned off; the vector V m and H m are the basis vector matrix and the projection matrix respectively.

[0026] Preferably, step S2 specifically comprises:

[0027] The solution of the ill-conditioned linear equations is converted into the infinite integral of the matrix exponential function by using the precise integration method. At the same time, the multi-window scanning technology is introduced on the basis of the preprocessing precise integration method. The relationship between the transient electromagnetic diffusion field and the virtual wave field is shown as follows:

[0028]

[0029] Where, f and u represent diffusion field and virtual wave field respectively; τ represents virtual time; exp represents exponential function; r represents position vector;

[0030] The linear discretization of formula (7) yields:

[0031]

[0032] in, Formula (8) is converted into the following formula:

[0033] Ku=f (9);

[0034] Since the precise integration method requires the coefficient matrix K to be a symmetric positive definite matrix, formula (9) is converted into the following formula:

[0035] K T Ku=K T f (10);

[0036] In which, let the matrix A = K T K, multiply both sides of formula (10) by the preconditioning matrix Q, and let u = Py, then formula (10) is converted to the following formula:

[0037] QAPy=QK T f (11);

[0038] Where Q and P are both diagonal matrices; formula (11) can be simplified as follows:

[0039] My = b (12);

[0040] Where M = QAP, b = QK T f; Through the preprocessing fine integration method, the solution of formula (12) is written as the following iterative formula:

[0041]

[0042] Among them, y k represents the intermediate variable; ζ is the integration step; F(ζ) has the following expression:

[0043] F(ζ)=ζ[I+(-Mζ) / 2+(-Mζ) 2 / 6+(-Mζ) 3 / 24] (14);

[0044] After y is obtained by iterative formula (13), the virtual wave field is finally obtained by u=Py.

[0045] Preferably, step S3 includes the following sub-steps:

[0046] S301, finally obtaining a migration imaging result of the electrical interface by using a finite difference migration imaging technique;

[0047] S302, superimposing the obtained migration images with different resolutions to obtain a high-resolution imaging section.

[0048] Preferably, step S301 is specifically as follows: In a rectangular coordinate system, the control equation of the passive area is as follows:

[0049]

[0050] Since upgoing waves are used for migration imaging, formula (15) is converted into a two-dimensional wave equation:

[0051]

[0052] Where, let x'=x, z'=z, t'=t+z / v, then u(x,z,t)=u(x',z',t'), and formula (16) is converted to the following formula:

[0053]

[0054] Substituting formula (17) into formula (15), we obtain the wave equation in the floating coordinate system:

[0055]

[0056] Performing Fourier transform on x' and t' in formula (18) yields the first-order approximate equation:

[0057]

[0058] Among them, the velocity v is a constant, or a function of x and z v(x,z). Then, a symmetric implicit difference format is used to obtain the expression of the overall sampling. The imaging data is obtained through formal conversion, and these data are arranged in the horizontal and depth directions to finally obtain the offset imaging results of the electrical interface.

[0059] Therefore, the present invention adopts the above-mentioned time-domain differential pulse transient electromagnetic multi-resolution migration imaging method to obtain high-resolution characteristics of electrical interfaces from shallow to deep, which better solves the multi-resolution detection problem.

[0060] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 It is a flow chart of an embodiment of a time domain differential pulse transient electromagnetic multi-resolution migration imaging method of the present invention;

[0062] Figure 2 It is a schematic diagram of a symmetric implicit difference format of an embodiment of a time domain differential pulse transient electromagnetic multi-resolution migration imaging method of the present invention;

[0063] Figure 3 It is a model diagram of a hidden thrust fault in an embodiment of a time domain differential pulse transient electromagnetic multi-resolution migration imaging method of the present invention;

[0064] Figure 4 It is a time domain differential pulse transient electromagnetic multi-resolution migration imaging method embodiment of the present invention, which has different pulse widths and time domain pseudo-random pulses; Figure 4 (a) The pulse width is 0.4ms; Figure 4 (b) The pulse width is 2ms; Figure 4 (c) The pulse width is 5ms;

[0065] Figure 5 It is a virtual wave field of pseudo-random pulses with different pulse widths in the time domain of a time domain differential pulse transient electromagnetic multi-resolution migration imaging method embodiment of the present invention; Figure 5 (a) The pulse width is 0.4ms; Figure 5 (b) The pulse width is 2ms; Figure 5 (c) The pulse width is 5ms;

[0066] Figure 6 It is a time domain differential pulse transient electromagnetic multi-resolution migration imaging method embodiment of the present invention, which is the migration imaging of pseudo-random pulses with different pulse widths in the time domain; Figure 6 (a) The pulse width is 0.4ms; Figure 6 (b) The pulse width is 2ms; Figure 6 (c) The pulse width is 5ms; Figure 6 (d) in the figure is the image after superposition. DETAILED DESCRIPTION

[0067] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0068] Unless otherwise defined, technical or scientific terms used in the present invention shall have the common meanings understood by one having ordinary skills in the field to which the present invention belongs.

[0069] Embodiment 1

[0070] like Figure 1 As shown, the present invention provides a time domain differential pulse transient electromagnetic multi-resolution migration imaging method, comprising the following steps:

[0071] Step S1: For the three-dimensional forward simulation of pseudo-random pulse transient electromagnetic in the time domain, the present invention adopts the unstructured tetrahedron vector finite element method. In order to achieve efficient numerical simulation, pseudo-random pulses with different pulse widths in the time domain are used to excite the electromagnetic field during the power supply period to obtain a multi-resolution transient electromagnetic diffusion field; specifically:

[0072] The second-order backward difference algorithm BDF2 and the displacement inverse-Krylov subspace algorithm are used during power supply and after the current is turned off, respectively. For a given calculation area Ω, the Maxwell equations ignoring the displacement current are shown as follows:

[0073]

[0074] Where, e is the electric field intensity; b is the magnetic induction intensity; μ0 is the vacuum magnetic permeability; J s is the current density; σ is the conductivity of the formation; represents the differential operator; t represents the diffusion time;

[0075] The space is discretized using an unstructured tetrahedral grid, the magnetic induction intensity b is eliminated from formulas (1a) and (1b), and the Galerkin method is used to obtain the governing equation of the electric field, as shown below:

[0076]

[0077] Among them, N i and N jis the interpolation basis function within the unit; V represents the volume of the three-dimensional space; dV represents the volume integral; the derivative term in formula (2) is obtained by the second-order backward difference method:

[0078]

[0079] Among them, e n and represent the electric field intensity and current density at the nth moment respectively.

[0080] Let the matrix M σ =∫ V μ0σN i ·N j dV,M inner =∫ V μ0N i ·N j dV, the electric field strength during power supply is solved by the following formula:

[0081]

[0082] After the current is turned off, J in formula (2) s is zero, rewrite it as:

[0083]

[0084] The displacement inverse-Krylov subspace algorithm is used to solve formula (5), and the approximate expression of the solution is:

[0085]

[0086] Where m is the order of the displacement inverse-Krylov subspace; γ is the displacement parameter; I is the unit matrix; e0 is the electric field intensity after the current is turned off; the vector V m and H m are the basis vector matrix and the projection matrix respectively.

[0087] Step S2, using the pre-processed fine integral wave field inverse transformation to extract the virtual wave field from the transient electromagnetic diffusion field; specifically:

[0088] The virtual wave field can be interpreted from the transient electromagnetic diffusion field to achieve high-resolution migration imaging of the electrical interface of the formation. The solution of the ill-conditioned linear equations is converted into the infinite integral of the matrix exponential function by the fine integration method, which greatly improves the numerical accuracy of the virtual wave field. Considering that the field data contains background noise, the multi-window scanning technology is introduced on the basis of the pre-processing fine integration method, which can effectively improve the noise resistance of the virtual wave field inverse transformation method.

[0089] The relationship between the transient electromagnetic diffusion field and the virtual wave field is shown in the following formula:

[0090]

[0091] Where, f and u represent diffusion field and virtual wave field respectively; τ represents virtual time; exp represents exponential function; r represents position vector;

[0092] The linear discretization of formula (7) yields:

[0093]

[0094] in, Formula (8) is converted into the following formula:

[0095] Ku=f (9);

[0096] The linear equation group (9) is highly ill-conditioned. The present invention adopts a pre-processing precise integration method to solve the virtual wave field u. Since the precise integration method requires the coefficient matrix K to be a symmetric positive definite matrix, formula (9) is converted into the following formula:

[0097] K T Ku=K T f (10);

[0098] In which, let the matrix A = K T K, multiply both sides of formula (10) by the preconditioning matrix Q, and let u = Py, then formula (10) is converted to the following formula:

[0099] QAPy=QK T f (11);

[0100] Where Q and P are both diagonal matrices; formula (11) can be simplified as follows:

[0101] My = b (12);

[0102] Where M = QAP, b = QK T f; Through the preprocessing fine integration method, the solution of formula (12) is written as the following iterative formula:

[0103]

[0104] Among them, y k represents the intermediate variable; ζ is the integration step; F(ζ) has the following expression:

[0105] F(ζ)=ζ[I+(-Mζ) / 2+(-Mζ) 2 / 6+(-Mζ) 3 / 24] (14);

[0106] After obtaining y through the iterative formula (13), the virtual wave field is finally obtained from u = Py.

[0107] Generally, in the case of small formation dip angles and simple structures, the isophase axes of the virtual wave field can intuitively reflect the characteristics of geological structures. However, for complex geological structures, the isophase axes of the virtual wave field in the deep space cannot accurately represent the formation morphology, and migration imaging methods need to be used to perform spatial migration on it.

[0108] Step S3: Use the finite-difference migration imaging technology to obtain the imaging profiles of pulses with different pulse widths in the time domain, and obtain the high-resolution characteristics of the electrical interfaces at different depths through the stacking technology. It includes the following sub-steps: S301: Finally obtain the migration imaging result of the electrical interface through the finite-difference migration imaging technology; specifically: in the rectangular coordinate system, the control equation in the source-free region is shown as follows:

[0109]

[0110] Since the up-going wave is used for migration imaging, formula (15) is converted into a two-dimensional wave equation:

[0111]

[0112] Among them, for equation (16), perform a coordinate transformation, let x' = x, z' = z, t' = t + z / v, then u(x,z,t) = u(x',z',t'), and formula (16) is converted into the following formula:

[0113]

[0114] Substitute formula (17) into formula (15) to obtain the wave equation in the floating coordinate system:

[0115]

[0116] In the floating coordinate system, the term is avoided, which can not only reduce the calculation time, but also maintain the stability of the time difference. Perform Fourier transforms on x' and t' in formula (18), and the first-order approximation equation is obtained:

[0117]

[0118] Among them, the velocity v is a constant, or a function v(x,z) of x and z.

[0119] Assume that equation (19) exists in the space Ω(-∞ < x < ∞, 0 ≤ t ≤ T, 0 ≤ z < Z), and the following mixed conditions are given:

[0120]

[0121] in, is the wave field observed on the ground. By solving the above partial differential equation, we can obtain the wave field value u(x, Z, t) of the upgoing wave that has appeared at any point (x, z>0) in the underground space.

[0122] The present invention adopts a symmetric implicit difference format, such as Figure 2 As shown, the following symbols are then used to derive the formula:

[0123]

[0124] Where Δx and Δz are the spatial steps in the x and z directions respectively, and Δt is the time step. Figure 6 The center point in is differentiated, and the wave field of the center point is assumed to be And order Then equation (19) can be simplified as:

[0125]

[0126] for make Then we have:

[0127]

[0128] In summary, equation (19) can be written as follows:

[0129]

[0130] set up Then the discrete expression of equation (24) can be simplified to:

[0131]

[0132] In order to obtain the expression for the sampling at the whole point, let Then equation (25) can be transformed into:

[0133]

[0134] After rearranging formula (26), we can obtain:

[0135]

[0136] Where i = 1, 2, ... nx, nx is the number of nodes in the x direction; j = 1, 2, ... nt, nt is the number of time channels at each measuring point. Formula (27) is written in matrix form:

[0137]

[0138] in, is the wave field vector along the x direction at the depth layer z = nΔz and the time layer t = jΔt; and is the wave field vector of the corresponding layer. Matrix A, B and vector and The expression is as follows:

[0139]

[0140] The vector in equation (28) and In the first downward extrapolation process, it is the transient electromagnetic virtual wave field observed on the ground. When extrapolating backward, the wave field of the ground receiving point can be used as the source, and the calculation Time Knows is zero. Therefore, as long as the matrix A is invertible, equation (28) has a solution, namely:

[0141]

[0142] Since A is a matrix with dominant principal diagonal and α is always greater than zero, there is always |1+2α|>2|α|, that is, detA≠0, so the matrix A is reversible.

[0143] The wave field on the grid node of the underground space can be obtained by reverse extension to the underground through the above method. Then, the wave field value of the grid node with a search time equal to 2z / v is used as imaging data, and these data are arranged in the horizontal direction and depth direction to finally obtain the imaging result of the electrical interface.

[0144] S302, superimposing the obtained migration images with different resolutions to obtain a high-resolution imaging section.

[0145] The results are verified by numerical examples.

[0146] The imaging method is verified by the hidden thrust fault model. The model schematic diagram and the resistivity of the formation are shown in Figure 2. Figure 3 As shown. The geological model mainly includes two main geological structures, namely the shallow ductile rock folds and the deep concealed thrust faults. The present invention uses a grounding wire as the transmitting source with a length of 4km. The coordinates of the two ends of the transmitting source are (-2km, 0m, 0m) and (2km, 0m, 0m). The receiving survey line is parallel to the transmitting source and has a length of 3km. The coordinates of the two end points of the survey line are (-1.5km, 100m, 0m) and (1.5km, 100m, 0m) respectively. Figure 3There are three main low-resistance formations, namely 100Ω·m, 50Ω·m and 20Ω·m. Through a large number of experiments, time domain pseudo-random pulses with pulse widths of 0.4ms, 2ms and 5ms are selected as the transmission waveform, such as Figure 4 shown.

[0147] In order to get close to the field data and image the electrical interface, 5% random noise is first added to the diffuse field, and then the virtual wave field is extracted from the diffuse field using the preprocessing fine integration method, such as Figure 5 shown. Figure 5 (a) and Figure 5 (b) in the figure can clearly depict the interface characteristics of shallow ductile rock folds, but it is difficult to reflect the interface information of deep hidden thrust faults. As the pulse width increases to 5ms, Figure 5 (c) in the figure can show the interface information of deep strata. The event axis of the virtual wave field can only intuitively reflect the morphology of shallow strata, but cannot accurately depict the characteristics of deep strata. Therefore, it is necessary to use the migration imaging method to spatially locate it.

[0148] Figure 6 It is the migration imaging of pseudo-random pulses with different pulse widths in the time domain. It can be seen that the migration imaging results of pseudo-random pulses with different pulse widths in the time domain can effectively characterize the interface morphology of strata at different depths, such as Figure 6 (a)- Figure 6 By superimposing the imaging results of pseudo-random pulses in different time domains, a high-resolution imaging profile can be obtained, as shown in (c). Figure 6 As shown in (d), the red dashed line represents the electrical interface of the theoretical model. However, due to the volume effect and other factors, the offset imaging of shallow wrinkles is different from the theoretical model.

[0149] Therefore, the present invention adopts the above-mentioned time-domain differential pulse transient electromagnetic multi-resolution migration imaging method to obtain high-resolution characteristics of electrical interfaces from shallow to deep, which better solves the multi-resolution detection problem.

[0150] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A time domain differential pulse transient electromagnetic multi-resolution migration imaging method, characterized in that: The following steps are involved: Step S1, during power supply, using time domain pseudo-random pulses with different pulse widths to excite the electromagnetic field to obtain a multi-resolution transient electromagnetic diffusion field; Step S2, extracting a virtual wave field from the transient electromagnetic diffusion field by using a pre-processed fine integral wave field inverse transformation; Step S3: Based on step S2, the imaging profiles of time-domain pulses with different pulse widths are obtained by using finite difference migration imaging technology, and the high-resolution features of electrical interfaces at different depths are obtained by using superposition technology.

2. The time domain differential pulse transient electromagnetic multi-resolution migration imaging method according to claim 1 is characterized in that: Step S1 is specifically as follows: The second-order backward difference algorithm BDF2 and the displacement inverse-Krylov subspace algorithm are used during power supply and after the current is turned off, respectively. For a given calculation area Ω, the Maxwell equations ignoring the displacement current are shown as follows: Where, e is the electric field intensity; b is the magnetic induction intensity; μ0 is the vacuum magnetic permeability; J s is the current density; σ is the conductivity of the formation; represents the differential operator; t represents the diffusion time; Eliminate the magnetic induction intensity b from formulas (1a) and (1b) and use the Galerkin method to obtain the governing equation of the electric field, as shown below: Among them, N i and N j is the interpolation basis function within the unit; V represents the volume of the three-dimensional space; dV represents the volume integral; the derivative term in formula (2) is obtained by the second-order backward difference method: Among them, e n and represent the electric field intensity and current density at the nth moment respectively.

3. The time domain differential pulse transient electromagnetic multi-resolution migration imaging method according to claim 2 is characterized in that: Let the matrix M σ =∫ V μ0σN i ·N j dV,M inner =∫ V μ0N i ·N j dV, the electric field strength during power supply is solved by the following formula: After the current is turned off, J in formula (2) s is zero, rewrite it as: The displacement inverse-Krylov subspace algorithm is used to solve formula (5), and the approximate expression of the solution is: Where m is the order of the displacement inverse-Krylov subspace; γ is the displacement parameter; I is the unit matrix; e0 is the electric field intensity after the current is turned off; the vector V m and H m are the basis vector matrix and the projection matrix respectively.

4. The time domain differential pulse transient electromagnetic multi-resolution migration imaging method according to claim 3 is characterized in that: Step S2 is specifically as follows: The solution of the ill-conditioned linear equations is converted into the infinite integral of the matrix exponential function by using the precise integration method. At the same time, the multi-window scanning technology is introduced on the basis of the preprocessing precise integration method. The relationship between the transient electromagnetic diffusion field and the virtual wave field is shown as follows: Where, f and u represent diffusion field and virtual wave field respectively; τ represents virtual time; exp represents exponential function; r represents position vector; The linear discretization of formula (7) yields: in, Formula (8) is converted into the following formula: Ku=f (9); Since the precise integration method requires the coefficient matrix K to be a symmetric positive definite matrix, formula (9) is converted into the following formula: K T Ku=K T f (10); In which, let the matrix A = K T K, multiply both sides of formula (10) by the preconditioning matrix Q, and let u = Py, then formula (10) is converted to the following formula: QAPy=QK T f (11); Where Q and P are both diagonal matrices; formula (11) can be simplified as follows: My = b (12); Where M = QAP, b = QK T f; Through the preprocessing fine integration method, the solution of formula (12) is written as the following iterative formula: Among them, y k represents the intermediate variable; ζ is the integration step; F(ζ) has the following expression: F(ζ)=ζ[I+(-Mζ) / 2+(-Mζ) 2 / 6+(-Mζ) 3 / 24] (14); After y is obtained by iterative formula (13), the virtual wave field is finally obtained by u=Py.

5. The time domain differential pulse transient electromagnetic multi-resolution migration imaging method according to claim 4 is characterized in that: Step S3 includes the following sub-steps: S301, finally obtaining a migration imaging result of the electrical interface by using a finite difference migration imaging technique; S302, superimposing the obtained migration images with different resolutions to obtain a high-resolution imaging section.

6. The time domain differential pulse transient electromagnetic multi-resolution migration imaging method according to claim 5, characterized in that: Step S301 is specifically as follows: In a rectangular coordinate system, the control equation of the passive area is as follows: Since upgoing waves are used for migration imaging, formula (15) is converted into a two-dimensional wave equation: Where, let x'=x, z'=z, t'=t+z / v, then u(x,z,t)=u(x',z',t'), and formula (16) is converted to the following formula: Substituting formula (17) into formula (15), we obtain the wave equation in the floating coordinate system: Performing Fourier transform on x' and t' in formula (18) yields the first-order approximate equation: Among them, the velocity v is a constant, or a function of x and z v(x,z). Then, a symmetric implicit difference format is used to obtain the expression of the overall sampling. The imaging data is obtained through formal conversion, and these data are arranged in the horizontal and depth directions to finally obtain the offset imaging results of the electrical interface.

Images