A three-dimensional controlled source electromagnetic forward modeling method and system

Abstract

The present invention discloses a three-dimensional controlled source electromagnetic forward modeling method and system, the method comprising: meshing a calculation area including a measurement area; then performing three-dimensional controlled source electromagnetic forward modeling based on the electromagnetic field excited by the field source to obtain the corresponding electric field intensity in the measurement area; wherein, the forward modeling process uses the Lanczos algorithm combined with a single optimized repeated pole to construct a rational Krylov subspace; then uses the rational Krylov subspace to construct a rational Lanczos approximation of the electric field intensity to be determined, and finally solves the electric field intensity based on the rational Lanczos approximation. The present invention further sets 4 to 6 layers of uniaxial anisotropic PML layers, and fixes the lowest frequency UPML parameters at the same time, ensuring that the mass matrix and the curl matrix of the forward modeling equation do not change with the frequency. The entire forward modeling process greatly reduces the memory consumption of the forward modeling solution, and can realize three-dimensional multi-frequency CSEM fast forward modeling solution on an ordinary notebook.

Description

Technical Field

[0001] The invention belongs to a frequency domain three-dimensional controlled source electromagnetic forward modeling method in the field of geophysical exploration, and specifically relates to a three-dimensional controlled source electromagnetic forward modeling method and system. Background Art

[0002] With the needs of national economic construction and the advancement of science and technology, electromagnetic exploration has received attention and development, and has played a great role in solving resource and environmental problems, determining deep structures, searching for deep ore bodies, and solving hydrological engineering problems.

[0003] At present, most of the one-dimensional or two-dimensional forward and inversion methods are used to interpret the controlled source field data. Due to the three-dimensional nature of the actual geological structure, one-dimensional and two-dimensional inversion are prone to interpretation deviations. The three-dimensional inversion interpretation depends on the calculation accuracy and speed of the three-dimensional forward modeling. Due to the complexity of the three-dimensional geological body, there is no analytical expression for the three-dimensional forward modeling problem, so it is necessary to use numerical discrete methods such as the finite element method to solve it.

[0004] In the forward modeling of controlled source electromagnetics, the number of frequencies determines the number of times the equations are solved. Existing three-dimensional forward modeling technology of multi-frequency controlled source electromagnetics often uses OMP and MPI-based solvers for parallel solving. When combined with efficient and stable direct solvers such as PARDISO or MUMPS for equation solving, high-speed forward calculation is obtained by consuming a large amount of memory. It is suitable for high-performance servers and computer clusters, but ordinary laptops are difficult to meet the computing requirements. Summary of the invention

[0005] The purpose of the present invention is to provide a three-dimensional controlled source electromagnetic forward modeling method and system for the problem that the above-mentioned existing high-speed forward modeling calculation consumes a large amount of memory, is suitable for high-performance servers and computer clusters, but is difficult to be applied to ordinary laptops. Among them, the three-dimensional controlled source electromagnetic forward modeling method constructs a Krylov subspace, and then projects a large-scale coefficient matrix into the subspace to obtain a reduced-order matrix of smaller dimensions. The inverse problem of the reduced-order matrix can be quickly solved, and then a rational approximation of the multi-frequency forward modeling problem can be quickly obtained. Through the above-mentioned conversion, the computational efficiency of the forward modeling process can be improved and memory consumption can be reduced.

[0006] In one aspect, the present invention provides a three-dimensional controlled source electromagnetic forward modeling method, which comprises the following steps:

[0007] Meshing the computational area including the measurement area;

[0008] Then, based on the electromagnetic field excited by the field source, a three-dimensional controlled source electromagnetic forward modeling is performed to obtain the corresponding electric field strength in the measurement area;

[0009] Among them, the forward modeling process is to use the Lanczos algorithm combined with a single optimized repeated pole to construct a rational Krylov subspace; then use the rational Krylov subspace to construct a rational Lanczos approximation of the electric field strength to be calculated, and finally solve the electric field strength based on the rational Lanczos approximation, and the rational Lanczos approximation is expressed as:

[0010] E≈||X|| M V m+1 g ω (A m+1 )e1

[0011] In the formula, T is the symbol of matrix transposition, E is the electric field strength to be determined, X = M -1 b, M is the mass matrix, b is the source term on the right side of the finite element linear equations of the controlled source electromagnetics, V m+1 =[v1,..,v m+1 ] represents a set of orthogonal basis vectors of the Krylov subspace, v1, v m+1 is the orthogonal basis vector V m+1 The first and m+1th orthogonal basis in , m is the dimension of the Krylov subspace, g ω (A m+1 ) is the transfer function, A m+1 A uses an orthogonal basis V m+1 The orthogonal projection matrix on the Krylov subspace is: A=M -1 C, C is the curl matrix.

[0012] In the rational Lanczos approximation constructed by the present invention, A m+1 A uses an orthogonal basis V m+1 In the Krylov subspace Q m+1 The orthographic projection matrix on A. m+1 The dimension of is much smaller than the dimension of A, and can be directly inverted using solvers such as PARDISO or MUMPS, thereby greatly reducing the amount of computation required to solve the equation. Therefore, the present invention constructs a Krylov subspace, and then projects a large-scale coefficient matrix onto the subspace to obtain a reduced-order matrix of smaller dimension. The inverse problem of the reduced-order matrix can be quickly solved, thereby quickly obtaining a rational approximation of the multi-frequency forward problem. Among them, especially when constructing a rational Krylov subspace in combination with a single optimized repeated pole, the constructed rational Krylov subspace is a Krylov subspace shared by multiple frequencies, and the construction of the subspace only requires solving an equation once, which can greatly save computing memory and improve computing efficiency.

[0013] Further optionally, the single optimal repeated imaginary pole in the constructed rational Krylov subspace is set as:

[0014]

[0015] In the formula, ξ j ,ξ0 is the jth pole and the 0th pole;ω min and ω max They are the minimum frequency value and the maximum frequency value within the calculation frequency range of the multi-frequency controllable source.

[0016] It should be understood that when the frequency range is relatively wide, the frequency may also be segmented, and a rational Krylov subspace is constructed using a single optimized repeated virtual pole for each frequency segment.

[0017] Further optionally, the constructed Krylov subspace is expressed as:

[0018]

[0019] In the formula, represents a rational Krylov subspace, m+1 is the dimension of the subspace, span{} represents the expansion space, ξ0 represents the 0th pole, ξ j represents the jth extreme point.

[0020] Further optionally, before meshing the calculation area including the measurement area, the method further includes:

[0021] Arranging multiple uniaxial anisotropic PML layers outside the measurement area to form a UPML medium, thereby forming a calculation area;

[0022] Among them, the magnetic permeability of the UPML medium Dielectric constant and conductivity It is expressed as:

[0023]

[0024] Where [Λ] is the principal axis anisotropy tensor matrix, μ0 is the vacuum magnetic permeability, ε0 ​​is the vacuum dielectric constant, σ is the dielectric conductivity, and s y ,s x ,s z All are UPML metrics;

[0025] When the electromagnetic wave is incident on the UPML medium perpendicular to the z-axis, s y =s x =1, When an electromagnetic wave is incident on a UPML medium perpendicular to the x-axis, s y =s z=1, When an electromagnetic wave is incident on the UPML medium perpendicular to the y-axis, s x =s z =1, k w , σ w , α w is the characteristic parameter that characterizes the absorption performance of the UPML layer, i is the imaginary unit, and ω is the angular frequency.

[0026] The present invention takes into account that in the conventional forward modeling method of the electromagnetic field, the background field method only needs to solve the abnormal electromagnetic field of each measuring point, and the initial field is loaded by one-dimensional forward modeling. Although it has the characteristics of high precision, it is not conducive to considering the posture of the actual field source and the calculation of complex terrain models. The total field method takes into account the loading method of the field source during the calculation process, which is closer to the layout of the actual field source, but the Dirichlet boundary condition is often used in numerical calculations. Although the boundary condition is easier to load, it only needs to expand the grid, but the grid boundary requires 5 to 10 times the skin depth to ensure that the electromagnetic wave is fully attenuated. Due to the loading of the Dirichlet boundary condition, the simulation area reaches tens of kilometers or even hundreds of kilometers, while the actual measurement area of ​​interest is often only a few kilometers, so this boundary causes too many grids at the boundary, resulting in too many unknowns and increased difficulty in solving equations.

[0027] In view of the above problems existing in the prior art, the present invention further utilizes the good absorption characteristics of uniaxial anisotropic perfectly matched layers (UPML), which can simulate unbounded regions through smaller boundaries and is very easy to implement in finite element programming. The present invention attempts to use UPML media to truncate the boundaries of controlled source electromagnetic forward modeling calculations, thereby reducing the scale of the equations and speeding up the solution of the equations.

[0028] Further optionally, the UPML medium is composed of 4-6 uniaxial anisotropic PML layers.

[0029] Further optionally, for a multi-frequency controllable source, the parameters of the UPML medium are fixed to the UPML parameters of the lowest frequency, and the mass matrix M and the curl matrix C at different frequencies are the mass matrix M and the curl matrix C at the lowest frequency.

[0030] In a second aspect, the present invention provides a system based on the above three-dimensional controlled source electromagnetic forward modeling method, comprising:

[0031] A meshing module, used to mesh the calculation area including the measurement area;

[0032] The forward modeling module is used to perform three-dimensional controllable source electromagnetic forward modeling based on the electromagnetic field excited by the field source to obtain the corresponding electric field strength in the measurement area;

[0033] Among them, the forward modeling process is to use the Lanczos algorithm combined with a single optimized repeated imaginary pole to construct a rational Krylov subspace; then use the rational Krylov subspace to construct a rational Lanczos approximation of the electric field strength to be calculated, and finally solve the electric field strength based on the rational Lanczos approximation, and the rational Lanczos approximation is expressed as:

[0034] E≈||X|| M V m+1 g ω (A m+1 )e1

[0035] In the formula, T is the symbol of matrix transposition, E is the electric field strength to be determined, X = M -1 b, M is the mass matrix, b is the right-hand source term, V m+1 =[v1,..,v m+1 ] represents a set of orthogonal basis vectors, v1, v m+1 is the orthogonal basis vector V m+1 The first and m+1th orthogonal basis in g ω (A m+1 ) is the transfer function, A m+1 A uses an orthogonal basis V m+1 In the Krylov subspace Q m+1 The orthogonal projection matrix on , A = M -1 C, C is the curl matrix.

[0036] In a third aspect, the present invention provides an electronic terminal, comprising:

[0037] one or more processors;

[0038] a memory storing one or more computer programs;

[0039] The processor calls the computer program to implement:

[0040] A three-dimensional controlled source electromagnetic forward modeling method.

[0041] In a fourth aspect, the present invention provides a readable storage medium storing a computer program, wherein the computer program is called by a processor to implement:

[0042] A three-dimensional controlled source electromagnetic forward modeling method.

[0043] Beneficial Effects

[0044] The invention provides a three-dimensional controlled source electromagnetic forward modeling method, which is based on a model reduction algorithm of a rational Krylov subspace. By selecting the optimal pole, a multi-frequency shared Krylov subspace is constructed, and then a large-scale coefficient matrix is ​​projected onto the subspace to obtain a reduced-order matrix of smaller dimension. The inverse problem of the reduced-order matrix can be quickly solved, and then a rational approximation of the multi-frequency forward modeling problem can be quickly obtained. When a single optimal pole is selected, the construction of the subspace only requires solving a first-order equation, which can greatly save computing memory and improve computing efficiency. In addition, the memory consumption is very small, and it is suitable for ordinary laptop computers, and can provide technical support for the three-dimensional rapid inversion of controlled source electromagnetics.

[0045] In a further preferred embodiment of the present invention, a UPML absorption boundary is set. Due to the good absorption performance of the UPML medium to electromagnetic waves, the entire calculation area is greatly reduced, which can improve the speed of solving the equation and further reduce the computer's memory consumption. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 Schematic diagram of the reduced calculation area.

[0047] Figure 2(a), Figure 2(b) and Figure 2(c) are the calculation results of the horizontal component of the electric field in a uniform half space.

[0048] Figure 3 It is a universal low-resistance anomaly model.

[0049] Figures 4(a), 4(b) and 4(c) are the comparison results between the method provided by the present invention and the conventional finite element calculation results, the frequency is 3 Hz. Figure 4(a) is the real part comparison, Figure 4(b) is the imaginary part comparison result, and Figure 4(c) is the relative error curve. DETAILED DESCRIPTION

[0050] The present invention provides a three-dimensional controlled source electromagnetic forward modeling method, the purpose of which is to improve the forward modeling efficiency and reduce memory consumption. To this end, a model reduction algorithm based on a rational Krylov subspace selects the optimal pole to construct a multi-frequency shared Krylov subspace, and then projects a large-scale coefficient matrix into the subspace to obtain a reduced-order matrix of smaller dimension. The inverse problem of the reduced-order matrix can be quickly solved, and then a rational approximation of the multi-frequency forward modeling problem is quickly obtained. Under the above core mechanism, the present invention further sets a UPML absorption boundary. Due to the good absorption performance of the UPML medium for electromagnetic waves, the entire calculation area is greatly reduced, which further improves the speed of solving the equation and reduces the computer's memory consumption. It should be understood that the forward modeling method of this embodiment will be described in combination with the above two optimization points.

[0051] About UPML layer:

[0052] In controlled source electromagnetic forward modeling, in order to simulate the actual unbounded area and ensure that the electromagnetic wave is fully attenuated, the Dirichlet boundary condition is often used. Although this method is easy to load, it only needs to expand the grid, but the grid boundary requires 5 to 10 times the skin depth, and the simulation area reaches tens of kilometers or even hundreds of kilometers. However, the actual measurement area of ​​interest is often only a few kilometers. Therefore, this boundary causes too many grids at the boundary, resulting in too many unknowns and difficulty in solving equations. In order to reduce the boundary, the present invention introduces UPML medium to truncate the calculation boundary.

[0053] In this embodiment, for the observation area / measurement area, 4 to 6 uniaxial anisotropic PML layers (UPML) are set outside the measurement network and field source positions calculated in the forward modeling. The truncation boundary should be a certain distance away from the measurement network and the field source and cannot be truncated directly. Therefore, this embodiment designs the number of PML layers (generally 4 to 6 layers are sufficient) and PML parameters according to the calculated frequency. In order to be applied to all frequencies, this embodiment fixes the parameters of the lowest frequency. In order to ensure the absorption effect of the PML layer on all frequencies, the frequency width should not be too large (the lowest and highest frequencies differ by 1 to 3 orders of magnitude). The following will provide a detailed description of the UPML medium parameters.

[0054] The magnetic permeability of the UPML introduced in the present invention is Dielectric constant and conductivity Its value affects the composition of the mass matrix M and the curl matrix C in the following formula (3), which is expressed as:

[0055]

[0056]

[0057] Where [Λ] is the principal axis anisotropy tensor matrix, μ0 is the vacuum magnetic permeability, ε0 ​​is the vacuum dielectric constant, σ is the dielectric conductivity, and s y ,s x ,s z All are UPML metrics;

[0058] When electromagnetic waves propagate in a three-dimensional medium, there are three different directions of UMPL / non-UPML interfaces, namely, interfaces perpendicular to the x-axis, y-axis, and z-axis. When an electromagnetic wave is incident on a UPML medium perpendicular to the z-axis, s y =s x =1, When an electromagnetic wave is incident on a UPML medium perpendicular to the x-axis, s y =s z =1, When an electromagnetic wave is incident on the UPML medium perpendicular to the y-axis, s x =s z=1, i is the imaginary unit, ω is the angular frequency. κ w , σ w , α w is a characteristic parameter to characterize the absorption performance of the UPML layer. As the UPML medium moves away from the UMPL / non-UPML region, κ w , σ w Gradually increases from 0 to the corresponding maximum value, α w It gradually decreases from the maximum value to 0, and its maximum value needs to be determined by numerical experiments. It should be understood that using the coefficients of the above formula, when an electromagnetic wave in any direction is incident on the UPML / non-UPML with the xy plane / xz plane / yz plane as the interface, it can enter the UPML layer without reflection.

[0059] After the UPML layer is introduced, the calculation area is meshed. In order to solve the electromagnetic field excited by the artificial field source, the finite element method is a commonly used numerical simulation method in the controlled source electromagnetic method (CSEM) forward modeling. For example, when using hexahedral meshes for spatial discretization, the air and the ground containing the anomaly are divided into hexahedral elements, and then the finite element linear equations of the controlled source electromagnetic are obtained as follows:

[0060] (C+iωM)·E=-iωb (3)

[0061] Where b is the source term on the right, M is the mass matrix, C is the curl matrix, and E is the electric field strength to be determined. i is the imaginary unit, and ω is the angular frequency. It should be understood that when the UPML layer is introduced, the permeability of the UPML layer is Dielectric constant and conductivity The specific structures of the mass matrix M and the curl matrix C are common knowledge in the art, so they are not described in detail.

[0062] As can be seen from the above formula, due to the introduction of UPML medium, the mass matrix M and the curl matrix C of different frequencies will change with the frequency. In order to construct the Krylov subspace of multiple frequencies, it is necessary to fix C and M. Then, the present invention ensures that the mass matrix M and the curl matrix C of the forward equation do not change with the frequency by fixing the UPML parameters of the lowest frequency. In this embodiment, the lowest frequency is set to a fixed frequency of ω0, and then equation (3) becomes:

[0063] (C(ω0)+iωM(ω0))·E=-iωb (4)

[0064] Right now

[0065] (C+iωM)·E=-iωb (5)

[0066] About rational Krylov subspace: The present invention uses Lanczos algorithm combined with a single optimized repeated imaginary pole to construct a rational Krylov subspace.

[0067] Among them, the rational Krylov subspace is expressed as:

[0068]

[0069] In the formula, represents the rational Krylov subspace, m+1 is the dimension of the subspace. J represents the jth vector of the Krylov subspace. span{} represents the expansion space. ξ0 represents the 0th pole, ξ j represents the jth pole. From equation (5), we can see that ξ j The number of determines the number of decompositions of the matrix equation. In order to reduce the decomposition of the matrix equation, the present invention uses a single optimized repeated imaginary pole (ω min and ω max are the minimum and maximum values ​​in the frequency range respectively), then only one matrix decomposition is needed to construct the subspace. Then the original coefficient matrix is ​​projected onto the subspace, which greatly reduces the dimension of the coefficient matrix.

[0070] In order to simplify the construction of the subspace, let A = M in equation (5) -1 C, X = M -1 b. Construct transfer function g ω (A) = -iω(A+iωI) -1 , I is the identity matrix, since the transfer function g ω (A) is about the inverse matrix of a large sparse matrix, which is quite difficult to solve directly. Therefore, the Krylov subspace method can be used to construct an efficient approximation r for such matrices. m (A), the electric field strength E is expressed as:

[0071] E≈r m (A)·X (7)

[0072] In the formula, q m (z)=(z-ξ1)(z-ξ2)...(z-ξ m ),ξ i (i=1,...,m) are m poles, p m is a non-zero polynomial of order m. Among them, solving r m (A) is equivalent to solving A in the m+1 dimensional Krylov subspace Since the coefficient matrix M is not a unit matrix, if the traditional Lanczos orthogonal algorithm is directly used to solve the orthogonal matrix, the calculation will become complicated. In order to make the calculation process relatively simple, an M-orthogonal Lanczos algorithm, that is, M-orthogonal operation, can be used to obtain the standard column orthogonal vector V m+1 The algorithm needs to define the corresponding induced norm and inner product. For any vectors x and y, we have:

[0073]

[0074] Among them, y H is the conjugate transpose of y.

[0075] Use the M-orthogonal Lanczos algorithm to construct the orthogonal basis V m+1 The process can be described by the recursive expression given in equation (9). First, let the initial column vector v1 = X / ||X|| M , and the subsequent m column vectors v j+1 , (j=1,...,m) will be generated in each recursive calculation step, so that a set of orthogonal basis vectors V of the rational Krylov subspace can be constructed m+1 =[v1,..,v m+1 ]. The recursive expression is as follows:

[0076]

[0077] Among them, ξ j is the extreme point, h i,j =<(A-ξ j I) -1 v j , v i > M ,h i,j > 0, so that ||v j+1 || M =1. And the standard orthogonal matrix constructed by this method satisfies: Among them I m+1 is the unit matrix. Then the rational Lanczos approximation (7) of equation (3) can be expressed as:

[0078] E≈||X|| M V m+1 g ω (A m+1 )e1. (10)

[0079] In the formula, The projection matrix is ​​expressed as: A m+1 A uses an orthogonal basis V m+1Orthogonal projection matrix on Krylov subspace. A m+1 The dimension of X is much smaller than that of A, and can be directly inverted using solvers such as PARDISO or MUMPS, thereby greatly reducing the amount of computation required to solve the equation. The calculation of X is essentially the solution of an equation, namely, MX = b. Although it involves the inversion of the M matrix, since the M matrix is ​​a positive definite symmetric sparse matrix, the calculation using the Pardiso solver often only takes a few seconds. V m+1 is a rectangular matrix, which is gradually expanded during the construction of the subspace and directly stored in the memory for easy access. m+1 Then, the Pardiso solver is used to solve g ω (A m+1 ). The remaining operations are all vector and matrix product operations, which can be quickly solved using interfaces such as self-written programs, compiler built-in functions, and MKL library functions.

[0080] Based on the above principle, in this embodiment, a UPML layer is set for the measurement area to form a calculation area, a grid is divided for the calculation area, and then a three-dimensional controlled source electromagnetic forward modeling is performed for the electromagnetic field excited by the artificial field source to obtain the corresponding electric field strength in the measurement area. The forward modeling process is to solve the electric field strength through formula (10).

[0081] It should be understood that in order to improve the computing speed and reduce memory consumption, this embodiment calculates a single optimal pole according to the frequency range, assembles a sparse matrix of the finite element linear equation system, and uses the Lanczos orthogonalization algorithm to construct a Krylov subspace (10 to 30 dimensions). Then calculate the projection of the coefficient matrix A on the subspace to obtain the reduced-order matrix A m+1 , and finally solve formula (10) to get the electric field strength.

[0082] In some feasible embodiments, the present invention further provides a system based on the above three-dimensional controlled source electromagnetic forward modeling method, which includes a meshing module and a forward modeling module. The meshing module is used to mesh the calculation area including the measurement area, and the forward modeling module is used to perform three-dimensional controlled source electromagnetic forward modeling based on the electromagnetic field excited by the field source to obtain the corresponding electric field strength in the measurement area.

[0083] In some implementation processes, before meshing, the meshing module further arranges multiple layers of uniaxial anisotropic PML layers outside the measurement area to form a UPML medium, thereby forming a calculation area.

[0084] The forward modeling module is used to deduce the electric field intensity based on the formula (10) constructed based on the constructed Krylov subspace. The specific formula deduction process refers to the content of the aforementioned method.

[0085] Among them, please refer to the corresponding process of the aforementioned method for the specific implementation process of each unit module. It should be understood that the specific implementation process of the above unit modules refers to the method content, and the present invention will not be specifically repeated here, and the division of the above functional module units is only a division of logical functions. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. At the same time, the above integrated unit can be implemented in the form of hardware or in the form of software functional units.

[0086] In some implementations, the present invention provides an electronic terminal comprising: one or more processors and a memory storing one or more computer programs, wherein the computer programs in the memory are called by the processor to implement: steps of a three-dimensional controlled source electromagnetic forward modeling method.

[0087] The specific implementation process is as follows: meshing the calculation area including the measurement area, and performing three-dimensional controlled source electromagnetic forward modeling based on the electromagnetic field excited by the field source to obtain the corresponding electric field strength in the measurement area.

[0088] In some implementation processes, before meshing, the meshing module further arranges multiple layers of uniaxial anisotropic PML layers outside the measurement area to form a UPML medium, thereby forming a calculation area.

[0089] Based on the electromagnetic field excited by the field source, the three-dimensional controlled source electromagnetic forward modeling is performed to obtain the corresponding electric field strength in the measurement area. The electric field strength is deduced based on the formula (10) constructed based on the constructed Krylov subspace. The specific formula deduction process refers to the content of the above method.

[0090] The terminal also includes: a communication interface for communicating with external equipment and performing data interactive transmission.

[0091] The memory may include a high-speed RAM memory, and may also include a non-volatile memory, such as at least one disk memory.

[0092] If the memory, processor and communication interface are implemented independently, the memory, processor and communication interface can be connected to each other through a bus and communicate with each other. The bus can be an industrial standard architecture bus, a peripheral device interconnection bus or an extended industrial standard architecture bus, etc. The bus can be divided into an address bus, a data bus, a control bus, etc.

[0093] Optionally, in a specific implementation, if the memory, processor and communication interface are integrated on a chip, the memory, processor and communication interface can communicate with each other through an internal interface.

[0094] For the specific implementation process of each step, please refer to the description of the above method.

[0095] It should be understood that in the embodiments of the present invention, the processor may be a central processing unit (CPU), and the processor may also be other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc. The memory may include a read-only memory and a random access memory, and provide instructions and data to the processor. A portion of the memory may also include a non-volatile random access memory. For example, the memory may also store information about the type of device.

[0096] In some implementations, a readable storage medium is provided, which stores a computer program, and the computer program is called by a processor to implement: steps of a three-dimensional controlled source electromagnetic forward modeling method.

[0097] The specific implementation process is as follows: meshing the calculation area including the measurement area, and performing three-dimensional controlled source electromagnetic forward modeling based on the electromagnetic field excited by the field source to obtain the corresponding electric field strength in the measurement area.

[0098] In some implementation processes, before meshing, the meshing module further arranges multiple layers of uniaxial anisotropic PML layers outside the measurement area to form a UPML medium, thereby forming a calculation area.

[0099] Based on the electromagnetic field excited by the field source, the three-dimensional controlled source electromagnetic forward modeling is performed to obtain the corresponding electric field strength in the measurement area. The electric field strength is deduced based on the formula (10) constructed based on the constructed Krylov subspace. The specific formula deduction process refers to the content of the above method.

[0100] For the specific implementation process of each step, please refer to the description of the above method.

[0101] The readable storage medium is a computer-readable storage medium, which may be an internal storage unit of the controller described in any of the foregoing embodiments, such as a hard disk or memory of the controller. The readable storage medium may also be an external storage device of the controller, such as a plug-in hard disk, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), etc. equipped on the controller. Further, the readable storage medium may also include both an internal storage unit of the controller and an external storage device. The readable storage medium is used to store the computer program and other programs and data required by the controller. The readable storage medium may also be used to temporarily store data that has been output or is to be output.

[0102] Based on this understanding, the technical solution of the present invention, in essence or in other words, the part that contributes to the prior art, or all or part of the technical solution can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned readable storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, etc., various media that can store program codes.

[0103] Effect verification:

[0104] like Figure 1 As shown, the present invention forms a UPML medium by setting Figure 1 The area shown is reduced, the calculation area including the measurement area and the field source is reduced, and the external expansion area is greatly reduced.

[0105] As shown in FIGS. 2(a), 2(b) and 2(c), the observation area is -1 km ≤ x ≤ 1 km, 4 km ≤ y ≤ 6 km, the resistivity of the homogeneous half-space is 100 Ω·m, and the 100-m field source is arranged along the x-axis at the origin of coordinates. After loading the PML, the grid boundary is reduced from -50 km < x, y, z < 50 km to -6 km ≤ x ≤ 6 km, -3 km ≤ y ≤ 10 km, -5 km ≤ z ≤ 10 km. The number of grid unknowns is reduced from 553834 to 270828. Forward modeling calculations are performed at 21 frequencies. The calculation time of the method provided by the present invention is only 50 s, and the memory consumption is 3.2 G. The calculation time of the conventional finite element method is 2427 s, and the memory consumption is 8.5 G. The computer is a personal PC (6 cores, 16 G). FIG. 2(a) is at 1 Hz, FIG. 2(b) is at 2 Hz, and FIG. 2(c) is at 50 Hz, showing the relative error distributions between the method provided by the present invention, the conventional finite element method and the analytical solution. It shows that the method provided by the present invention can meet the error requirements.

[0106] As Figure 3 shown, the source length is 100 m, the center coordinates of the low-resistivity anomaly are (1000 m, 0 m, 250 m), the resistivity is 5 Ω·m, and the background resistivity is 50 Ω·m. After loading the PML medium, the grid boundary is reduced to -3 km ≤ x ≤ 4 km, -4 km ≤ y ≤ 4 km, -5 km ≤ z ≤ 10 km. The number of grid unknowns is reduced from 662823 to 352067. Forward modeling calculations are also performed at 21 frequency points. The calculation time of the conventional finite element method is 2867 s, and the memory consumption is 10.5 G; the calculation time of the method provided by the present invention is only 54 s, and the memory consumption is 4.2 G.

[0107] As shown in FIGS. 4(a), 4(b) and 4(c), the comparison results between the method provided by the present invention and the calculation results of the conventional finite element method at a frequency of 3 Hz. FIG. 4(a) shows the comparison of the real part, FIG. 4(b) shows the comparison of the imaginary part, and FIG. 4(c) shows the relative error curve.

[0108] It should be emphasized that the examples described in the present invention are illustrative rather than restrictive. Therefore, the present invention is not limited to the examples described in the specific embodiments. Any other embodiments obtained by those skilled in the art according to the technical solutions of the present invention, without departing from the spirit and scope of the present invention, whether modified or replaced, also belong to the protection scope of the present invention.

Claims

1. A three-dimensional controlled source electromagnetic forward modeling method, Features: The following steps are involved: Meshing the computational area including the measurement area; Then, based on the electromagnetic field excited by the field source, a three-dimensional controlled source electromagnetic forward modeling is performed to obtain the corresponding electric field strength in the measurement area; Among them, the forward modeling process is to use the Lanczos algorithm combined with a single optimized repeated pole to construct a rational Krylov subspace; then use the rational Krylov subspace to construct a rational Lanczos approximation of the electric field strength to be calculated, and finally solve the electric field strength based on the rational Lanczos approximation, which is expressed as: E≈||X|| M V m+1 g ω (A m+1 )e 1 In the formula, T is the symbol of matrix transposition, E is the electric field strength to be determined, X = M -1 b, M is the mass matrix, b is the source term on the right side of the finite element linear equations of the controlled source electromagnetics, V m+1 =[v 1 ,..,v m+1 ] represents a set of orthogonal basis vectors of the Krylov subspace, v 1 , v m+1 is the orthogonal basis vector V m+1 The first and m+1th orthogonal basis in, m+1 is the dimension of the Krylov subspace, g ω (A m+1 ) is the transfer function, A m+1 A uses an orthogonal basis V m+1 The orthogonal projection matrix on the Krylov subspace is: A=M -1 C, C is the curl matrix.

2. The method according to claim 1, Features: The single optimal repeated imaginary pole in the constructed rational Krylov subspace is set as: In the formula, ξ j ,ξ 0 is the jth pole and the 0th pole; ω min and ω max They are the minimum frequency value and the maximum frequency value within the calculation frequency range of the multi-frequency controllable source.

3. The method according to claim 1, Features: The constructed Krylov subspace is expressed as: In the formula, represents the rational Krylov subspace, m+1 is the dimension of the subspace, span{} represents the expansion space, ξ 0 represents the 0th pole, ξ j represents the jth extreme point.

4. The method according to claim 1, Features: Before meshing the calculation area including the measurement area, the method further includes: Arranging multiple uniaxial anisotropic PML layers outside the measurement area to form a UPML medium, thereby forming a calculation area; Among them, the magnetic permeability of the UPML medium Dielectric constant and conductivity It is expressed as: Where [Λ] is the principal axis anisotropy tensor matrix, μ 0 is the vacuum permeability, ε 0 is the vacuum dielectric constant, σ is the dielectric conductivity, s y ,s x ,s z All are UPML metrics; When the electromagnetic wave is incident on the UPML medium perpendicular to the z-axis, s y =s x =1, When an electromagnetic wave is incident on a UPML medium perpendicular to the x-axis, s y =s z =1, When an electromagnetic wave is incident on the UPML medium perpendicular to the y-axis, s x =s z =1, k w , σ w , α w is the characteristic parameter that characterizes the absorption performance of the UPML layer, i is the imaginary unit, and ω is the angular frequency.

5. The method according to claim 4, Features: The UPML medium is composed of 4-6 uniaxial anisotropic PML layers.

6. The method according to claim 4, Features: For a multi-frequency controllable source, the parameters of the UPML medium are fixed to the UPML parameters of the lowest frequency, and the mass matrix M and the curl matrix C at different frequencies are the mass matrix M and the curl matrix C at the lowest frequency.

7. A system based on the method according to any one of claims 1 to 6, Features: include: A meshing module, used to mesh the calculation area including the measurement area; The forward modeling module is used to perform three-dimensional controllable source electromagnetic forward modeling based on the electromagnetic field excited by the field source to obtain the corresponding electric field strength in the measurement area; Among them, the forward modeling process is to use the Lanczos algorithm combined with a single optimized repeated imaginary pole to construct a rational Krylov subspace; then use the rational Krylov subspace to construct a rational Lanczos approximation of the electric field strength to be calculated, and finally solve the electric field strength based on the rational Lanczos approximation, and the rational Lanczos approximation is expressed as: E≈||X|| M V m+1 g ω (A m+1 )e 1 In the formula, T is the symbol of matrix transposition, E is the electric field strength to be determined, X = M -1 b, M is the mass matrix, b is the right-hand source term, V m+1 =[v 1 ,..,v m+1 ] represents a set of orthogonal basis vectors, v 1 , v m+1 is the orthogonal basis vector V m+1 The first and m+1th orthogonal basis in g ω (A m+1 ) is the transfer function, A m+1 A uses an orthogonal basis V m+1 In the Krylov subspace Q m+1 The orthogonal projection matrix on , A = M -1 C, C is the curl matrix.

8. An electronic terminal, Features: include: one or more processors; a memory storing one or more computer programs; The processor calls the computer program to implement: The steps of the method according to any one of claims 1 to 6.

9. A readable storage medium, Features: A computer program is stored, which is called by a processor to implement: The steps of the method according to any one of claims 1 to 6.

Images