A half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source

Abstract

The application relates to the field of geophysical numerical calculation, and particularly discloses a half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source. The half-space transient electromagnetic three-dimensional forward method based on the full-space initial field source is composed of six parts, i.e. an initial field source, a control equation, a three-dimensional discrete format, an iterative initial time and a time step, a model grid division mode, a ground-air boundary condition and a boundary condition of other boundaries. The ground-air boundary condition adopts an air layer setting method, thereby realizing forward calculation of a model with terrain and an aerial exploration model. The initial field source adopts a full-space initial field source, the allowed adjacent grid step ratio is increased from 1.2 of a traditional method to 2, and the allowed maximum and minimum grid step ratio is increased from 20 of the traditional method to be unlimited, so that the calculation amount is greatly reduced, the calculation efficiency is improved, and forward calculation of small-scale, large-scale and multi-scale models can be realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of geophysical numerical calculation, and particularly relates to a half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source. BACKGROUND

[0002] The transient electromagnetic method (TEM) is a time-domain artificial source electromagnetic detection method based on the principle of electromagnetic induction. It is a kind of geophysical exploration method for finding various geological targets by using a non-grounded loop or grounded wire to emit a primary field, and measuring the change of the induced secondary field with time generated by the underground medium after the primary field is turned off. The transient electromagnetic method is widely used in the fields of metal mineral exploration, groundwater and geothermal resource development, oil and gas resource exploration, mine hidden disaster body exploration, engineering and environmental geophysical exploration, and marine geophysical exploration, and the forward method is the basis of related theoretical research and application.

[0003] The existing half-space transient electromagnetic three-dimensional time-domain finite difference forward method is based on a half-space initial field source, and has the following disadvantages:

[0004] (1) Because it is based on a half-space initial field source, the upward continuation method is usually used for the ground-air boundary condition, which makes it difficult to handle terrain models and airborne exploration models.

[0005] (2) When the ratio of adjacent grid steps is greater than 1.2, and the ratio of the maximum and minimum grid steps is greater than 20, the iterative calculation result will diverge, resulting in incorrect calculation results. Therefore, it is usually required that the ratio of adjacent grid steps is less than 1.2, and the ratio of the maximum and minimum grid steps is limited to 20, which results in a large number of model grid points and large calculation amount, making it difficult to simulate small-scale and multi-scale models. SUMMARY

[0006] The present application aims to provide a half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source to solve the problems in the background art.

[0007] To achieve the above-mentioned purpose, the present application provides the following technical solutions:

[0008] A half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source, the method comprising: an initial field source, control equations and three-dimensional discrete format, iteration initial time and time step, model grid division method, boundary conditions of ground-air boundary and boundary conditions of other boundaries; the initial field source adopts a full-space initial field source.

[0009] As a further limitation of the technical scheme of the present application: the type of the initial field source is a full-space magnetic dipole source, or a full-space electric dipole source, or a superposition source of multiple full-space magnetic dipole sources, or a superposition source of multiple full-space electric dipole sources.

[0010] As a further limitation of the technical scheme of the present application: the resistivity parameter in the initial field source calculation formula adopts the earth resistivity value below the transmitting loop, or any resistivity value greater than the earth resistivity value.

[0011] As a further limitation of the technical scheme of the present application: the control equation and the three-dimensional discrete format adopt the time-domain finite difference method of the explicit time step method.

[0012] As a further limitation of the technical scheme of the present application: the resistivity parameter in the iteration initial time calculation formula and the time step calculation formula adopts the air layer resistivity value above the transmitting loop.

[0013] As a further limitation of the technical scheme of the present application: the model grid division method adopts a non-uniform grid.

[0014] As a further limitation of the technical scheme of the present application: the adjacent grid step length ratio in the non-uniform grid is not greater than 2, and the maximum and minimum grid step length ratio is not limited.

[0015] As a further limitation of the technical scheme of the present application: the boundary condition of the ground-air boundary adopts the method of setting an air layer.

[0016] Compared with the prior art, the method has the following advantages:

[0017] (1) Because the air-ground boundary processing method of setting an air layer is adopted, the forward calculation of the terrain model and the airborne exploration model can be realized.

[0018] (2) Because the full-space initial field source is adopted, the calculation error introduced by upward continuation calculation is avoided, the stability of the algorithm is greatly improved, the allowable adjacent grid step length ratio is increased from 1.2 of the traditional method to 2, and the maximum and minimum grid step length ratio is increased from 20 of the traditional method to unlimited. Therefore, a large model can be constructed with fewer grid numbers, thereby greatly reducing the calculation amount, improving the calculation efficiency, and realizing small-scale, large-scale and multi-scale model forward calculation. BRIEF DESCRIPTION OF DRAWINGS

[0019] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the embodiments or prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application.

[0020] Figure 1The application provides a half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source.

[0021] Figure 2 The application provides a method and a digital filtering method.

[0022] Figure 3 The application provides a relative error between a result of the method and a digital filtering solution.

[0023] Figure 4 The application provides a valley model diagram.

[0024] Figure 5 The application provides a hill model diagram.

[0025] Figure 6 The application provides a result of the method for the valley model and the hill model.

[0026] Figure 7 The application provides an airborne transient electromagnetic exploration model diagram.

[0027] Figure 8 The application provides a result of the method for the airborne transient electromagnetic exploration model.

[0028] Figure 9 The application provides a result of the method for a small-scale model, a large-scale model and a multi-scale model. DETAILED DESCRIPTION

[0029] Embodiment 1

[0030] In order to make the technical problems, technical solutions and beneficial effects of the application clearer, the application is further described in detail below in combination with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the application and do not limit the application.

[0031] Figure 1 The application provides a half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source, which comprises six parts of an initial field source, a control equation and a three-dimensional discrete format, an iterative initial time and a time step, a model grid division mode, a boundary condition of a ground-air boundary and a boundary condition of other boundaries.

[0032] The technical solution of the application is the same as that of the traditional method in the two aspects of the control equation and the three-dimensional discrete format and the boundary condition of other boundaries, and is different from that of the traditional method in the four aspects of the initial field source, the iterative initial time and the time step, the model grid division mode and the boundary condition of the ground-air boundary.

[0033] Specifically,

[0034] (1) Initial field source adopts full-space initial field source, that is, the uniform full-space electromagnetic field at an initial time after the transmitting current is turned off. The full-space initial field source is used for half-space forward.

[0035] (2) The type of initial field source is full-space magnetic dipole source, or full-space electric dipole source, or superposition source of multiple full-space magnetic dipole sources, or superposition source of multiple full-space electric dipole sources. The reference calculation formula of the full-space magnetic dipole source is as follows:

[0036]

[0037] In the formula, (x, y, z) is the position of the observation point, is the probability integral, ρ is the resistivity value of the uniform full-space medium.

[0038] (3) The resistivity parameter in the calculation formula of the initial field source adopts the earth resistivity value below the transmitting loop.

[0039] (4) The control equation and its three-dimensional finite difference discrete format are the same as the traditional method, and the time domain finite difference method with an explicit time step is adopted.

[0040] (5) The resistivity parameter in the calculation formula of the iteration initial time and the time step adopts the air layer resistivity value above the transmitting loop. The reference calculation formula of the iteration initial time is as follows:

[0041]

[0042] In the formula, μ0 is the vacuum permeability, σ1 is the air layer conductivity (the inverse of resistivity), and Δ1 is the grid step length at the loop.

[0043] The reference calculation formula of the time step is as follows:

[0044]

[0045] In the formula, α is a coefficient, Δ min is the minimum grid step length, μ0 is the vacuum permeability, σ1 is the air layer conductivity (the inverse of resistivity), and t is the sampling time of the last iteration step.

[0046] (6) The model grid division method adopts non-uniform grid, and the ratio of adjacent grid step lengths is not greater than 2, and the ratio of the maximum and minimum grid step lengths is not limited.

[0047] (7) The boundary condition of the ground-air boundary adopts the method of setting an air layer.

[0048] (8) The boundary condition of other boundaries is the same as the traditional method, and the Dirichlet, PML and other boundary conditions are adopted.

[0049] Based on the above method, the forward calculation of transient electromagnetic field can be completed through a certain calculation process. The calculation process includes setting the observation system and model parameters, calculating the initial field, iteratively calculating the secondary field, calculating the boundary condition, and outputting the calculation result when the number of iterations reaches the set maximum number of iterations. The present application realizes the forward calculation of the terrain model and the airborne exploration model by adopting the air-ground boundary processing method of setting the air layer. By adopting the full-space initial field source, the calculation error introduced by upward continuation calculation is avoided, the stability of the method is greatly improved, the allowed adjacent grid step ratio is increased from 1.2 of the traditional method to 2, and the allowed maximum and minimum grid step ratio is increased from 20 of the traditional method to unlimited. Therefore, the method of the present application can construct a large-scale model with fewer grid numbers, thereby greatly reducing the calculation amount, improving the calculation efficiency, and realizing the forward calculation of small-scale, large-scale and multi-scale models.

[0050] Example 2

[0051] In order to illustrate the calculation accuracy of the method of the present application, the following examples are given.

[0052] The model is a half-space model containing a low-resistance layer, the low-resistance layer has a depth of 50 m, a thickness of 30 m, a resistivity of 10 Ω·m, and a background resistivity of 100 Ω·m. The observation system is: the transmitting magnetic moment is 1 A·m 2 , the receiving point is located at the center of the transmitting loop, and the sampling time is 10 ms. The forward calculation is performed by using the method of the present application and the digital filtering method respectively.

[0053] Please refer to Figure 2 and Figure 3 , Figure 2 show the calculation results of the method of the present application and the digital filtering method, Figure 3 show the relative error of the calculation results of the method of the present application and the digital filtering solution. From Figure 2 and Figure 3 it can be seen that the calculation results of the method of the present application are consistent with the digital filtering solution, the relative error is within 3%, and the method of the present application has high calculation accuracy.

[0054] Example 3

[0055] In order to illustrate the applicability of the method of the present application in the forward calculation of the terrain model, the following examples are given.

[0056] The model is a valley model as shown in Figure 4 , the valley width is 260 m, the height is 70 m, and the ground resistivity is 100 Ω·m. The hill model is as shown in Figure 5 , the hill width is 260 m, the height is 70 m, and the ground resistivity is 100 Ω·m. The observation system is: the transmitting loop is located at the middle of the valley bottom and the hill top, the transmitting magnetic moment is 1 A·m 2The receiving point is located at the center of the transmitting loop, and the sampling time is 3ms.

[0057] Figure 6 The figure shows the calculation results of the method of this invention. In addition to the valley model and hill model, the calculation results of the flat land model are also given in the figure. Figure 6 It can be seen that in 10 -4 After milliseconds, the secondary field of the valley model is influenced by the highlands on both sides, making it stronger than that of the flatland model; the secondary field of the hill model is influenced by the air on both sides, making it weaker than that of the flatland model. The calculation results show that the method of this invention is applicable to forward modeling with terrain.

[0058] Example 4

[0059] To illustrate the applicability of the method of the present invention in airborne transient electromagnetic detection, the following calculation example is given.

[0060] The model is as follows: Figure 7 As shown, there are three layers: the first layer has a depth of 0-100m and a resistivity of 3Ω·m; the second layer has a depth of 100-400m and a resistivity of 20Ω·m; and the third layer has a depth of 400m and below and a resistivity of 3Ω·m. The observation system consists of a transmitter loop 30m above the ground and a transmitting magnetic moment of 1A·m. 2 The receiving point is located at the center of the transmitting loop, and the sampling time is 10ms. Figure 8 The calculation results demonstrate that the method of this invention is applicable to forward modeling of aeronautical models.

[0061] Example 5

[0062] To illustrate the applicability of the method of the present invention in forward modeling of small-scale, large-scale, and multi-scale models, the following examples are given.

[0063] The model comprises three three-dimensional low-resistivity blocks of different scales: Block 1 has a depth of 0.3 m and dimensions of x×y×z = 0.6 m × 0.6 m × 0.4 m; Block 2 has a depth of 51.1 m and dimensions of x×y×z = 102.2 m × 102.2 m × 51.2 m; Block 3 has a depth of 1638.3 m and dimensions of x×y×z = 3276.6 m × 3276.6 m × 1638.4 m. The background resistivity is 100 Ω·m, and the block resistivity is 1 Ω·m. The observation system uses a emitted magnetic moment of 1 A·m. 2 The receiving point is located at the center of the sending return line.

[0064] Figure 9 The figure simultaneously displays the anomalous responses of three low-resistivity blocks to illustrate the calculation results. In this example, the target volume scale spans five orders of magnitude, from 10 cm to 1 km. This demonstrates that the method of this invention is applicable to forward modeling of small-scale, large-scale, and multi-scale models.

[0065] The forward parameter setting of the embodiment is shown in Table 1, the model grid quantity is x×y×z=40×40×40, the model size is x×y×z=209715m×209715m×209715m, the minimum grid step is 0.1m, the adjacent grid step ratio is 2, the maximum grid step is 52429m, the maximum and minimum grid step ratio reaches 524290, far exceeding the limit of the traditional method 20. When the minimum grid step of the traditional method takes the same value (0.1m) as the method, the maximum grid step is limited to not more than 2m, and the grid quantity required to build the same size model exceeds 100000×100000×100000. If the grid quantity of the traditional method is calculated as 100000×100000×100000, the grid quantity is 1.5625×1010 times the grid quantity (40×40×40) of the method, that is, the calculation amount is 1.5625×1010 times the method. Therefore, the method can greatly reduce the calculation amount and improve the calculation efficiency.

[0066] Table 1 Forward parameter of multi-scale model

[0067]

[0068] The above only describes the preferred embodiments of the present application and is not intended to limit the present application. Any modification, equivalent replacement and improvement within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A half-space transient electromagnetic three-dimensional forward method based on a full-space initial field source, characterized in that, The method comprises: initial field source, control equation and three-dimensional discrete format, iteration initial time and time step, model grid division mode, boundary condition of ground-air boundary and boundary condition of other boundaries; the initial field source adopts a full-space initial field source, that is, a uniform full-space electromagnetic field at an initial time after the emission current is turned off, the full-space initial field source is used for half-space forward, and the following calculation formula is used when the full-space initial field source is a full-space magnetic dipole source: where (x, y, z) is the observation point position, is the probability integral, ρ is the uniform full-space medium resistivity value, the resistivity parameter in the initial field source calculation formula adopts the earth resistivity value below the transmitting loop, or adopts any resistivity value greater than the earth resistivity value; the control equation and the three-dimensional discrete format adopt the time-domain finite difference method of the explicit time step method; the iteration initial time calculation formula adopts: where μ0 is the vacuum magnetic permeability, σ1 is the air layer conductivity, and Δ1 is the grid step length at the loop; the time step calculation formula adopts: where α is the coefficient, Δ min is the minimum grid step length, μ0 is the vacuum magnetic permeability, σ1 is the air layer conductivity, and t is the sampling time of the previous step iteration; the model grid division method adopts a non-uniform grid, and the adjacent grid step length ratio is not greater than 2, and the maximum and minimum grid step length ratio is not limited; the boundary condition of the ground-air boundary adopts the method of setting an air layer.

2. The full-space initial field source based half-space transient electromagnetic 3D forward method of claim 1, wherein, The type of the initial field source is a full-space magnetic dipole source, or a full-space electric dipole source, or a superposition source of multiple full-space magnetic dipole sources, or a superposition source of multiple full-space electric dipole sources.

Images