Scattering electromagnetic field and magnetotelluric field numerical simulation method based on coordinate scaling boundary condition, electronic equipment and computer readable storage medium

Abstract

The invention relates to a scattering electromagnetic field and magnetotelluric field numerical simulation method based on coordinate scaling boundary conditions. According to the method, the coordinate telescopic absorption boundary condition is added in finite difference numerical simulation solution of the magnetotelluric scattering secondary field to correct the control equation of the scattering secondary field, so that the influence of a truncation boundary can be effectively reduced, the solution space of model calculation is greatly compressed, and the efficiency of three-dimensional magnetotelluric forward modeling is improved; coordinate functions in different directions are set at different boundary layers through coordinate stretching and contracting absorbing boundary conditions, complex factors changing according to periods are adopted when longitudinal coordinate stretching and contracting factors of all the boundary layers are set, electromagnetic waves of different frequencies are rapidly absorbed, and the setting method further effectively compresses the size of a simulation grid area and improves the simulation accuracy. And the numerical simulation of the scattering secondary field and the numerical simulation of the magnetotelluric field can be quickly and accurately carried out.

Description

Technical Field

[0001] The present invention relates to the field of geophysical electromagnetic detection, and in particular to a numerical simulation method of a scattered electromagnetic field and a magnetotelluric field based on coordinate expansion and contraction boundary conditions, an electronic device, and a computer-readable storage medium. Background Art

[0002] Magnetotelluric (MT) is an important geophysical method that uses natural electromagnetic fields as field sources to study the electrical structure of underground media. It has been widely used in oil and gas, metal mineral exploration, and deep structural research. At present, with the development of computer technology and the needs of practical work, fast three-dimensional magnetotelluric forward and inversion simulation calculations have become a research hotspot. Among them, the staggered grid finite difference method has been widely used in three-dimensional magnetotelluric forward and inversion algorithms due to its advantages of fast calculation and strong adaptability. On the basis of the traditional finite difference method, the calculation speed and accuracy of three-dimensional electromagnetic field forward and inversion simulations can be effectively improved by improving the electromagnetic field distribution and introducing the electromagnetic field divergence correction method. In magnetotelluric numerical simulation, an important aspect is to set reasonable boundary conditions, please refer to Figure 1 The spatial range of the three-dimensional magnetotelluric solution is limited by six interfaces, the top is the air top boundary, the bottom is the deep bottom boundary, and the rest are four side boundaries. In order to meet these boundary conditions and avoid the influence of truncation boundaries, it is generally necessary to set a "sufficiently large" three-dimensional space in the three-dimensional electromagnetic field forward and inversion simulation to ensure that the normal derivative of the electromagnetic field component on the side boundary is 0, the abnormal electromagnetic field on the bottom boundary is 0, etc., which results in a large amount of physical storage space and a long calculation time.

[0003] To solve the above problems, most scholars have introduced the Perfectly Matched Layer (PML) absorbing boundary condition and the Anisotropic Perfectly Matched Layer (APML) absorbing boundary theory for electromagnetic field simulation calculations. Studies have shown that the APML absorbing boundary is more efficient and more stable than Berenger's PML absorbing boundary, and can be applied to electromagnetic field forward modeling of any frequency. However, although the existing PML absorbing boundary condition can effectively compress the simulation grid space and reduce the physical storage space, its application effect in three-dimensional magnetotelluric forward modeling is affected by many factors, such as the thickness of the boundary, the size of the boundary grid, and the difference in conductivity between the anomaly and the background field, especially the thickness of the boundary has a very large impact on the results; on the other hand, the addition of the PML boundary condition will cause the increase of the coefficient matrix condition number, which may lead to an increase in the number of iterative calculations and an increase in calculation time. Summary of the invention

[0004] Based on this, the purpose of the present invention is to provide a scattered electromagnetic field numerical simulation method based on coordinate telescopic boundary conditions, which has the advantages of effectively compressing the size of the magnetotelluric numerical simulation space and reducing the calculation time without affecting the calculation accuracy.

[0005] A numerical simulation method for scattered electromagnetic fields based on coordinate expansion and contraction boundary conditions comprises the following steps:

[0006] S1 decomposes the magnetotelluric field into the background primary field and the scattered secondary field;

[0007] S2 introduces coordinate stretching absorption boundary conditions into the scattered secondary field to correct the scattered secondary field. The corrected scattered secondary field control equation is:

[0008]

[0009] in,

[0010] Where: E p It represents the component of the earth's electric field in the background primary field, E s It represents the component of the earth electric field in the scattered secondary field, H s represents the component of the Earth's magnetic field in the scattered secondary field, i is an imaginary unit, ω is the angular frequency, μ is the magnetic permeability, σ is the electrical conductivity, represents the unit vector in the direction of the X axis, Represents the unit vector in the Y-axis direction, Represents the unit vector in the Z-axis direction, f x is the coordinate scaling factor in the X-axis direction, f y is the coordinate scaling factor in the Y-axis direction, f z is the coordinate scaling factor in the Z-axis direction;

[0011] S3 sets coordinate expansion factors in different directions for coordinate expansion absorbing boundary conditions and synthesizes coefficient matrices;

[0012] S4 loads a secondary field source according to the coefficient matrix and uses a finite difference numerical simulation method to solve the scattered secondary field.

[0013] Furthermore, the coordinate scaling factor f in the X-axis direction is x , the coordinate scaling factor f in the Y-axis direction y , the coordinate scaling factor f in the Z-axis direction z Respectively meet:

[0014] f x =f x (x);

[0015] f y =f y(y);

[0016] f z =f z (z);

[0017] In the formula, f x (x) is the coordinate function in the X-axis direction, f y (y) is the coordinate function in the Y-axis direction, f z (z) is the coordinate function in the Z-axis direction.

[0018] Furthermore, the coordinate expansion and contraction absorbing boundary condition is:

[0019] For the boundary layer perpendicular to the X-axis, its f x =a+iblog(T),f y =f y (y) = f z =f z (z);

[0020] For the boundary layer perpendicular to the Y axis, its f x =f x (x) = f z =f z (z), f y =a+iblog(T);

[0021] For the boundary layer perpendicular to the Z axis, its f x =f x (x) = f y =f y (y), f z =a+iblog(T).

[0022] Furthermore, the coordinate scaling factors in different directions in step S3 are set as follows:

[0023] For the boundary layer perpendicular to the X-axis, its f x =a+iblog(T),f y =f z =1;

[0024] For the boundary layer perpendicular to the Y axis, its f x =f z =1,f y =a+iblog(T);

[0025] For the boundary layer perpendicular to the Z axis, its f x =f y =1,f z =a+iblog(T).

[0026] Furthermore, the coefficient matrix satisfies:

[0027]

[0028] Where Λ represents the coefficient matrix.

[0029] Compared with the prior art, the numerical simulation method of scattered electromagnetic fields based on coordinate telescopic boundary conditions described in the present invention can effectively reduce the influence of truncation boundaries, greatly compress the solution space of model calculation, and improve the efficiency of three-dimensional magnetotelluric forward simulation by adding coordinate telescopic absorption boundary conditions in the finite difference numerical simulation solution of magnetotelluric scattered secondary fields to correct the control equation of the scattered secondary fields. The coordinate telescopic absorption boundary conditions are used to set coordinate functions in different directions in different boundary layers, and a complex factor that varies according to the period is used when setting the longitudinal coordinate telescopic factor of each boundary layer, so as to realize rapid absorption of electromagnetic waves of different frequencies. The setting method further effectively compresses the size of the simulation grid area, and can quickly and accurately perform numerical simulations of scattered secondary fields and magnetotelluric fields.

[0030] At the same time, the present invention also provides a method for numerical simulation of magnetotelluric fields based on coordinate expansion and contraction boundary conditions, comprising:

[0031] S5 solves the background primary field, and adds the background primary field to the scattered secondary field solved by any of the above-mentioned scattered electromagnetic field numerical simulation methods based on coordinate expansion and contraction boundary conditions to obtain a three-dimensional magnetotelluric field.

[0032] The beneficial effects of the numerical simulation method of the magnetotelluric field based on the coordinate expansion and contraction boundary condition described in the present invention are the same as those of the numerical simulation method of the scattered electromagnetic field based on the coordinate expansion and contraction boundary condition, which will not be described in detail here.

[0033] For better understanding and implementation, the present invention is described in detail below with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 This is a schematic diagram of the three-dimensional magnetotelluric numerical simulation grid;

[0035] Figure 2 for Figure 1 Schematic diagram of a grid unit in ;

[0036] Figure 3 is a flow chart of the present invention;

[0037] Figure 4 A schematic diagram of an absorption boundary of an example of the present invention;

[0038] Figure 5 A three-dimensional geoelectric model established based on the present invention;

[0039] Figure 6 To adopt Figure 5 The calculation results of magnetotelluric numerical simulation using a three-dimensional geoelectric model. DETAILED DESCRIPTION

[0040] The technical solution of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention.

[0041] When numerically simulating three-dimensional electromagnetic fields in a limited area, boundary conditions must be given at the truncated boundaries of the calculation area. The processing of boundary conditions is an important factor affecting the accuracy and efficiency of electromagnetic field calculations. In the secondary field method, when the six boundaries are far enough from the secondary field source, the field value of the magnetotelluric secondary field decays to zero at the boundary. Therefore, the boundary conditions of the secondary field method are relatively simple, but in order to meet the boundary conditions, the boundary needs to be set at "far enough" to avoid the influence of the truncated boundary, which requires a relatively large calculation area. In order to effectively reduce the solution space of three-dimensional magnetotelluric numerical simulation, the present application sets a coordinate telescopic boundary layer at the truncated boundary of the simulation grid area, so that when the secondary field propagates to the truncated boundary, it enters the coordinate telescopic boundary layer without reflection; this setting method effectively compresses the size of the simulation grid area, and can realize fast and accurate magnetotelluric numerical simulation calculations.

[0042] See also Figure 3 The numerical simulation method of scattered electromagnetic field and magnetotelluric field based on coordinate expansion and contraction boundary conditions of the present invention comprises:

[0043] S1 decomposes the magnetotelluric field into the background primary field and the scattered secondary field.

[0044] In the specific implementation, according to the superposition principle, the magnetotelluric field (E, H) is decomposed into the background primary field (E p ,H p ) and the scattered secondary field (E s ,H s ), specifically:

[0045]

[0046] In the formula, E p It represents the component of the earth's electric field in the background primary field, E s It represents the component of the earth's electric field in the scattered secondary field; H p represents the component of the Earth's magnetic field in the background primary field, H s Represents the component of the Earth's magnetic field in the scattered secondary field.

[0047] S2 introduces a coordinate stretching absorption boundary condition into the scattered secondary field to correct the scattered secondary field.

[0048] In specific implementation, according to the frequency range of magnetotelluric method, the displacement current is ignored and the time harmonic factor e is taken iωt, then the frequency domain Maxwell equations are:

[0049]

[0050] Where i is the imaginary unit, ω is the angular frequency, μ is the magnetic permeability, σ is the electrical conductivity, E is the electric field strength, and H is the magnetic field strength.

[0051] Substituting the background primary field and scattered secondary field decomposed in step S1 into the frequency domain Maxwell equations, the governing equation of the magnetotelluric scattered secondary field is:

[0052]

[0053] Introducing the coordinate expansion and contraction absorption boundary condition into the control equation of the secondary field of magnetotelluric scattering, the correction control method of the secondary field of magnetotelluric scattering is:

[0054]

[0055] in,

[0056] f x =f x (x), f y =f y (y), f z =f z (z)

[0057] In the formula, represents the unit vector in the direction of the X axis, Represents the unit vector in the Y-axis direction, Represents the unit vector in the Z-axis direction, f x is the coordinate scaling factor in the X-axis direction, f y is the coordinate scaling factor in the Y-axis direction, f z is the coordinate scaling factor in the Z-axis direction. Each coordinate scaling factor corresponds to a coordinate function in a different direction. x (x) is the coordinate function in the X-axis direction, f y (y) is the coordinate function in the Y-axis direction, f z (z) is the coordinate function in the Z-axis direction.

[0058] S3 sets coordinate scaling factors in different directions for the coordinate scaling absorbing boundary conditions to synthesize the coefficient matrix Λ.

[0059] In order to ensure that the magnetotelluric field enters the boundary layer without reflection and decays rapidly in the boundary layer, and to avoid the influence of the truncation boundary on the magnetotelluric numerical simulation calculation, experimental research shows that when the media on both sides of the boundary layer have the same conductivity and magnetic permeability and the transverse coordinate expansion factor f x 、fy When , the electromagnetic field can enter the boundary layer without reflection; further, in order to make the electromagnetic waves entering the absorbing boundary layer decay quickly, the longitudinal coordinate scaling factor f in the boundary layer is z Set to plural, specifically:

[0060] f z =a+iblog(T)

[0061] Where a and b are constants, and T is the period.

[0062] Scale the vertical coordinate by factor f z Set to a complex number adjusted according to the period change. When the period is small, the corresponding longitudinal coordinate scaling factor f z smaller; when the period is larger, the corresponding longitudinal coordinate expansion factor f z It becomes larger and can achieve the effect of quickly absorbing electromagnetic waves of different frequencies.

[0063] For example, see Figure 3 When the magnetotelluric field is incident on the boundary layer in the Z direction at any angle, the conductivity of the medium outside the boundary layer is σ1, the magnetic permeability is μ1, and the transverse coordinate expansion factor f x1 f y1 The conductivity of the medium in the boundary layer is σ2, the magnetic permeability is μ2, and the transverse coordinate expansion factor f x2 f y2 , so that the conductivity, magnetic permeability and transverse coordinate expansion factor on both sides of the boundary layer interface are the same, then σ1=σ2, μ1=μ2, f x1 =f x2 、f y1 =f y2 At the same time, the longitudinal coordinates in the boundary layer are scaled by the factor f z2 Set to complex and make f z1 =f z2 , specifically:

[0064] f z2 =a+iblog(T).

[0065] The coefficient matrix Λ satisfies:

[0066]

[0067] Specifically, in the three-dimensional electromagnetic numerical simulation, there are two boundary layers in the X-axis direction, the Y-axis direction, and the Z-axis direction. The values ​​of the coordinate scaling factors in the coefficient matrix for different directions are shown in Table 1:

[0068] Table 1

[0069]

[0070] S4 loads a secondary field source according to the coefficient matrix Λ, and uses a finite difference numerical simulation method to solve the scattered secondary field.

[0071]

[0072] In the formula, Represents the secondary field source.

[0073] S5 solves the background primary field, and adds the background primary field to the scattered secondary field solved in step S4 to obtain a three-dimensional magnetotelluric field.

[0074] In the specific implementation, the existing magnetotelluric layered medium analytical solution is used to accurately solve the background primary field.

[0075] Three-dimensional magnetotelluric numerical simulation test analysis

[0076] See also Figure 5 , Figure 5 A three-dimensional geoelectric model with three anomalies in a uniform half space with a resistivity of 100Ωm established based on the present invention:

[0077] ρ1 is a 10Ωm low-resistance anomaly with a size of 40km×5km×15km and a burial depth of 5km;

[0078] ρ2 is a 1Ωm low-resistance anomaly with dimensions of 15km×20km×5km and a burial depth of 15km;

[0079] ρ3 is a 10000Ωm high-resistance anomaly with dimensions of 15km×20km×30km and a burial depth of 20km.

[0080] See also Figure 6 , Figure 6 To adopt Figure 5 The calculation results of the three-dimensional geoelectric model shown in the figure are the numerical simulation results of magnetotelluric, among which Miensopust (finite difference method) is the reference result and SFM (quadratic field finite difference method) is the result of the present invention. Figure 6 It can be seen that the calculation results of the quadratic field method SFM are in good agreement with those of other classical algorithms, which further illustrates the reliability of the quadratic field algorithm SFM.

[0081] Based on the same inventive concept, the present application also provides an electronic device, which may be a terminal device such as a server, a desktop computing device or a mobile computing device (e.g., a laptop computing device, a handheld computing device, a tablet computer, a netbook, etc.). The device includes one or more processors and a memory, wherein the processor is used to execute a program to implement the scattered electromagnetic field simulation method based on coordinate telescopic boundary conditions and the magnetotelluric field numerical simulation method based on coordinate telescopic boundary conditions of the above-mentioned embodiment of the present invention; and the memory is used to store a computer program executable by the processor.

[0082] Based on the same inventive concept, the present application also provides a computer-readable storage medium, on which a computer program is stored, and the program is executed by a processor to generate instructions for the scattered electromagnetic field simulation method based on coordinate telescopic boundary conditions and the magnetotelluric field numerical simulation method based on coordinate telescopic boundary conditions in the above-mentioned embodiment.

[0083] The present application may take the form of a computer program product implemented on one or more storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing program code. Computer-usable storage media include permanent and non-permanent, removable and non-removable media, and information storage can be achieved by any method or technology. Information can be computer-readable instructions, data structures, modules of programs, or other data. Examples of computer storage media include but are not limited to: phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, read-only compact disk read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cassettes, tape disk storage or other magnetic storage devices or any other non-transmission medium that can be used to store information that can be accessed by a computing device.

[0084] The terms used in the embodiments of the present application are only for the purpose of describing specific embodiments, and are not intended to limit the embodiments of the present application. The singular forms of "a", "said" and "the" used in the embodiments of the present application and the claims are also intended to include plural forms, unless the context clearly indicates other meanings. It should also be understood that, unless otherwise specified, "multiple" and "several" refer to two or more; "and / or" refers to and includes any or all possible combinations of one or more associated listed items; "first", "second", "third", etc. are only used to distinguish, rather than to describe a specific order or sequence, and cannot be understood as indicating or implying relative importance. When the above description relates to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. In the description of the present application, for those of ordinary skill in the art, the specific meanings of the above terms in the present application can be understood according to the specific circumstances.

[0085] The above-mentioned embodiments only express several implementation methods of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the invention patent. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, and the present invention is also intended to include these modifications and modifications.

Claims

1. A numerical simulation method for scattered electromagnetic fields based on coordinate expansion and contraction boundary conditions, characterized in that: The following steps are involved: S1 decomposes the magnetotelluric field into the background primary field and the scattered secondary field; S2 introduces coordinate stretching absorption boundary conditions into the scattered secondary field to correct the scattered secondary field. The corrected scattered secondary field control equation is: in, Where: E p It represents the component of the earth's electric field in the background primary field, E s It represents the component of the earth electric field in the scattered secondary field, H s represents the component of the Earth's magnetic field in the scattered secondary field, i is an imaginary unit, ω is the angular frequency, μ is the magnetic permeability, σ is the electrical conductivity, represents the unit vector in the direction of the X axis, Represents the unit vector in the Y-axis direction, Represents the unit vector in the Z-axis direction, f x is the coordinate scaling factor in the X-axis direction, f y is the coordinate scaling factor in the Y-axis direction, f z is the coordinate scaling factor in the Z-axis direction; S3 sets coordinate expansion factors in different directions for coordinate expansion absorbing boundary conditions and synthesizes coefficient matrices; S4 loads a secondary field source according to the coefficient matrix and uses a finite difference numerical simulation method to solve the scattered secondary field.

2. The scattered electromagnetic field numerical simulation method based on coordinate expansion boundary conditions according to claim 1 is characterized in that: The coordinate scaling factor f in the X-axis direction x , the coordinate scaling factor f in the Y-axis direction y , the coordinate scaling factor f in the Z-axis direction z Respectively meet: f x =f x (x); f y =f y (y); f z =f z (z); In the formula, f x (x) is the coordinate function in the X-axis direction, f y (y) is the coordinate function in the Y-axis direction, f z (z) is the coordinate function in the Z-axis direction.

3. The method for numerical simulation of scattered electromagnetic fields based on coordinate expansion and contraction boundary conditions according to claim 2, characterized in that: The coordinate expansion and contraction absorbing boundary condition is: For the boundary layer perpendicular to the X-axis, its f x =a+iblog(T),f y =f y (y) = f z =f z (z); For the boundary layer perpendicular to the Y axis, its f x =f x (x) = f z =f z (z), f y =a+iblog(T); For the boundary layer perpendicular to the Z axis, its f x =f x (x) = f y =f y (y), f z =a+iblog(T).

4. The method for numerical simulation of scattered electromagnetic fields based on coordinate expansion and contraction boundary conditions according to claim 3 is characterized in that: The coordinate scaling factors in different directions in step S3 are set as follows: For the boundary layer perpendicular to the X-axis, its f x =a+iblog(T),f y =f z =1; For the boundary layer perpendicular to the Y axis, its f x =f z =1,f y =a+iblog(T); For the boundary layer perpendicular to the Z axis, its f x =f y =1,f z =a+iblog(T).

5. The scattered electromagnetic field numerical simulation method based on coordinate expansion boundary conditions according to any one of claims 1 to 4, characterized in that: The coefficient matrix satisfies: Where Λ represents the coefficient matrix.

6. The method for numerical simulation of scattered electromagnetic fields based on coordinate expansion and contraction boundary conditions according to claim 5, characterized in that: The quadratic field is solved as follows: In the formula, Represents the secondary field source.

7. A method for numerical simulation of magnetotelluric field based on coordinate expansion boundary conditions, characterized in that: include: S5 solves the background primary field, and adds the background primary field to the scattered secondary field solved by the scattered electromagnetic field numerical simulation method based on coordinate expansion and contraction boundary conditions described in any one of claims 1 to 6 to obtain a three-dimensional magnetotelluric field.

8. The method for numerical simulation of magnetotelluric field based on coordinate expansion boundary condition according to claim 7, characterized in that: The background primary field is solved using the magnetotelluric layered medium analytical solution.

9. An electronic device, comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the method for numerical simulation of scattered electromagnetic fields based on coordinate telescopic boundary conditions as described in any one of claims 1 to 6 is implemented, and / or the method for numerical simulation of magnetotelluric fields based on coordinate telescopic boundary conditions as described in claim 7 or 8 is implemented.

10. A computer-readable storage medium storing computer-executable instructions, characterized in that: The computer can execute the scattered electromagnetic field numerical simulation method based on coordinate telescopic boundary conditions as described in any one of claims 1 to 6, and / or implement instructions generated by the magnetotelluric field numerical simulation method based on coordinate telescopic boundary conditions as described in claim 7 or 8.

Images