A full-zone apparent resistivity imaging method based on sensitivity of ground-to-air frequency domain electromagnetic method

Abstract

The application provides a ground-air frequency domain electromagnetic full-area apparent resistivity imaging method based on sensitivity, which comprises uniform half-space discretization, calculation of underground electric field response value based on magnetic vector potential function, calculation of air magnetic field response based on tensor Green function, calculation of coupling sensitivity based on air vertical magnetic field response, calculation of full-area apparent resistivity based on quasi-static approximation formula of vertical magnetic induction intensity of uniform ground space, and calculation of full-area apparent resistivity combined with each system parameter. The application solves the problem of volume effect influence on the existing ground-air frequency domain data interpretation effect, and through the calculation of the coupling sensitivity of the underground grid, the calculation of the full-area apparent resistivity based on the coupling sensitivity can be realized, the coupling relationship between the measured magnetic field at the receiving point and the underground grid is considered, the volume effect is corrected, the apparent resistivity imaging precision is improved, and then the ground-air frequency domain data interpretation effect is improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of geophysical exploration, and relates to a full-area apparent resistivity imaging method, in particular to a full-area apparent resistivity imaging method based on sensitivity in the ground-air frequency domain electromagnetic system. BACKGROUND

[0002] The ground-air frequency domain electromagnetic system lays a transmitting source on the ground, and an unmanned flight platform carries a receiving device to collect the magnetic field component in the air, so as to realize the detection of underground electrical structure. Due to the separated device form, the ground-air frequency domain electromagnetic system can realize rapid non-contact continuous signal collection in complex terrain areas. The ground-air frequency domain electromagnetic system has large transmitting power, large exploration depth and high safety, and has gradually become a research hotspot in electromagnetic detection in recent years and is widely applied to production practice.

[0003] In order to realize the detection of underground electrical structure, the collected magnetic field data need to be interpreted to obtain the underground apparent resistivity distribution map. The apparent resistivity imaging has the characteristics of simple and efficient processing flow and is widely applied to the interpretation of ground-air frequency domain electromagnetic data. However, the ground-air frequency domain detection system separates the transmitting device and the receiving device, and the collected magnetic field data contain the influence of volume effect. The traditional apparent resistivity imaging based on apparent depth does not correct the influence of volume effect, so the interpretation accuracy is limited, which leads to inaccurate calculation of the ground-air frequency domain apparent resistivity. SUMMARY

[0004] In order to solve the above technical problems and improve the apparent resistivity imaging accuracy of the ground-air frequency domain electromagnetic detection system, the application provides a full-area apparent resistivity imaging method based on sensitivity in the ground-air frequency domain electromagnetic system, which comprises the following steps:

[0005] Step one, the uniform half space is gridded, a space rectangular coordinate system is established with the transmitting source as the origin, and grid division is performed, wherein the length of the region in the x, y and z directions is obtained by analyzing the curve of the magnetic field response with the length of the underground divided region, and the length of the grid in the x, y and z directions is obtained by calculating the magnetic field response curve under different grid lengths and taking the grid length when the magnetic field response tends to be stable as the length of the underground grid division;

[0006] Step two, the electric field response value of each grid underground is calculated based on the magnetic vector potential function, and the calculation formula is as follows:

[0007]

[0008]

[0009] In the formula, A x , A zrespectively, are the x and z components of the geomagnetic vector potential function; I1 is the electric dipole moment, J0(λρ) is the zeroth order Bessel function, λρ is the argument of the Bessel function, ρ = (x 2 +y 2 ) 1 / 2 , (x, y, z) is the coordinate of the observation point, λ is the Hankel integral variable; k0 and k1 are the wave numbers of the air and the underground medium, respectively, k0 = ω(ε0μ0) 1 / 2 , k1 = (iωμ0σ1) 1 / 2 ; μ0 is the vacuum permeability, μ0 = 4π × 10 -7 H / m;

[0010] Applying the Foster identity:

[0011]

[0012] where I0 and K0 are the first and second kind modified Bessel functions, respectively, r = (ρ 2 +z 2 ) 1 / 2 ;

[0013] Applying the Sommerfeld identity:

[0014]

[0015]

[0016]

[0017]

[0018] In the above formula, N and P represent the Foster integral formula and the Sommerfeld integral formula, respectively;

[0019] The electric field expression of each grid in the underground is obtained from the relationship between the magnetic vector potential and the electromagnetic field:

[0020]

[0021] Step three, based on the tensor Green function, the magnetic field response of each grid at the observation point (receiving point) in the air is calculated:

[0022] Substitute the electric field response of each grid in the underground obtained from formula (7) into the tensor Green function:

[0023] H s (r r ) = ∫∫∫ V' G H (r r ​r)J(r)dv' (8)

[0024] wherein

[0025]

[0026] The current density J(r) is calculated from E calculated by formula (7), wherein r r and r represent the position vectors of the aerial observation point (receiving point) and the underground grid point, respectively, to obtain the magnetic field response of each grid point at the aerial observation point (receiving point);

[0027] Step four, the coupling sensitivity of each grid with respect to the aerial observation point is calculated based on the magnetic field response of each grid at the aerial observation point (receiving point):

[0028]

[0029] wherein H z (r r , f) is the secondary magnetic field z component generated by the underground grid point r at the aerial receiving point r r when the frequency is f, i, j, and k are the grid numbers in the x direction, the y direction, and the z direction, respectively.

[0030] Step five, the apparent resistivity of the whole area is calculated by bisection method iteration from the observed magnetic field response at the aerial observation point combined with the quasi-static approximation formula of the spatial vertical magnetic induction intensity of the uniform earth:

[0031]

[0032] wherein J1(λρ) is the first-order Bessel function, and (x, y, z) is the aerial receiving point coordinate.

[0033] Step six, the apparent resistivity of the whole area based on the coupling sensitivity weighting of each grid underground is calculated from the apparent resistivity of the whole area calculated in step five combined with the coupling sensitivity distribution under the same system parameters in step four, and the calculation formula is: (xi,yj,zk)

[0034]

[0035] wherein D n is the nth observation point, f Q is the Qth frequency point, ρ (Dn,fQ) is the apparent resistivity of the whole area calculated by iterative calculation of the magnetic field response at the observation point D n when the frequency is f Q is the apparent resistivity of the whole area calculated by iterative calculation of the magnetic field response at the observation point D n when the frequency is f Q under the system parameters of (D i ​​y j ,z k ) is the sensitivity coefficient at the receiving point, (x i ,y j ,z k ) is the underground grid point coordinate.

[0036] Further, the region x direction edge length is 4000-4500m, the y direction edge length is 2500-3000m outside the transmitting-receiving distance, and the z direction edge length is 2500-3000m in step one.

[0037] The present application has the following advantages:

[0038] The present application provides a full-area apparent resistivity imaging method based on sensitivity of ground-to-air frequency domain electromagnetic, which calculates the sensitivity of different underground regions to the magnetic field at the receiving point to represent the contribution of the magnetic field at the receiving point; when calculating the underground apparent resistivity from the magnetic field, the sensitivity is used for weighted summation of the apparent resistivity, which solves the problem of volume effect affecting the interpretation effect in the existing ground-to-air frequency domain data interpretation, and through the calculation of the coupling sensitivity of the underground grid, the full-area apparent resistivity calculation based on the coupling sensitivity can be realized, the coupling relationship between the measured magnetic field at the receiving point and the underground grid is considered, the volume effect is corrected, the apparent resistivity imaging precision is improved, and the ground-to-air frequency domain data interpretation effect is improved. BRIEF DESCRIPTION OF DRAWINGS

[0039] Figure 1 It is a method flowchart of the present application;

[0040] Figure 2 It is a grid diagram of the uniform half-space of the present application;

[0041] Figure 3 It is a curve diagram of the change of the magnetic field response with the edge length of the divided region of the present application;

[0042] Figure 4 It is a curve diagram of the change of the magnetic field response with the edge length of the discretization grid of the present application. DETAILED DESCRIPTION

[0043] The present application provides a full-area apparent resistivity imaging method based on sensitivity of ground-to-air frequency domain electromagnetic, which comprises the following steps:

[0044] Step one, grid the uniform half-space, establish a space rectangular coordinate system for the underground region with the transmitting source as the origin, and perform grid division, wherein through analyzing the curve of the change of the magnetic field response with the edge length of the underground divided region, the edge length of the region with the stable magnetic field response is taken as the x, y and z direction edge length of the underground divided region; through calculating the magnetic field response curve under different grid edge lengths, the grid edge length when the magnetic field response tends to be stable is taken as the underground grid division edge length;

[0045] The area x direction side length is 4000m, the y direction side length is 3000m outside the distance, and the z direction side length is 3000m. The grid length is 50m.

[0046] Step two, based on the magnetic vector potential function, the electric field response value of each grid in the underground is calculated, and the calculation formula is as follows:

[0047]

[0048]

[0049] In the formula, A x and A z are the x component and z component of the magnetic vector potential function in the ground respectively; Il is the electric dipole moment, J0(λρ) is the zero order Bessel function, λρ is the Bessel function argument, ρ=(x 2 +y 2 ) 1 / 2 , (x, y, z) is the coordinate of the observation point in the underground, and λ is the Hankel integral variable; k0 and k1 are the wave numbers of the air and the underground medium respectively, k0=ω(ε0μ0) 1 / 2 , k1=(iωμ0σ1) 1 / 2 ; μ0 is the vacuum magnetic permeability, μ0=4π×10 -7 H / m;

[0050] The Foster identity is applied:

[0051]

[0052] In the formula, I0 and K0 are the first type and the second type modified Bessel functions respectively, r=(ρ 2 +z 2 ) 1 / 2 ;

[0053] The Sommerfeld identity is applied:

[0054]

[0055] , we get:

[0056]

[0057]

[0058] In the above formula, N and P represent the Foster integral formula and the Sommerfeld integral formula respectively;

[0059] The electric field expression of each grid in the subsurface is obtained from the relationship between the magnetic vector potential and the electromagnetic field:

[0060]

[0061] Step three, the magnetic field response of each grid at the observation point (receiving point) in the air is calculated based on the tensor Green function:

[0062] The electric field response of each grid in the subsurface obtained from equation (7) is substituted into the tensor Green function:

[0063] H s (r r )=∫∫∫ V' G H (r r ,r)J(r)dv' (8)

[0064] In the formula:

[0065]

[0066] The current density J(r) is calculated from the E calculated by equation (7), and r r and r respectively represent the position vectors of the observation point (receiving point) in the air and the grid point in the subsurface, and the magnetic field response of each grid point at the observation point (receiving point) in the air is obtained;

[0067] Step four, the coupling sensitivity of each grid with respect to the observation point in the air is calculated based on the magnetic field response of each grid at the observation point (receiving point) in the air:

[0068]

[0069] In the formula, H z (r r ,r,f) is the secondary magnetic field z component produced by the subsurface grid point r at the air receiving point r r when the frequency is f, i, j, and k are the grid numbers in the x direction, y direction, and z direction, respectively.

[0070] Step five, the observed magnetic field response at the observation point in the air is combined with the quasi-static approximation formula of the spatial vertical magnetic induction intensity of the uniform earth, and the apparent resistivity of the whole area is calculated by bisection method iteration:

[0071]

[0072] In the formula, J1(λρ) is the first order Bessel function, and (x, y, z) is the coordinate of the receiving point in the air.

[0073] Step six, the full-zone apparent resistivity calculated from step five, combined with the coupling sensitivity distribution under the same system parameters obtained from step four, calculates the full-zone apparent resistivity weighted by the coupling sensitivity at each grid in the subsurface ρ (xi,yj,zk) , the calculation formula is:

[0074]

[0075] In the formula, D n is the nth observation point, f Q is the Qth frequency point, ρ (Dn,fQ) is the full-zone apparent resistivity obtained by iterative calculation of the magnetic field response at the observation point D n , the frequency is f Q , is the sensitivity coefficient at the grid (x i , y j , z k ) under the (D n , f Q ) system parameters, and (x i , y j , z k ) is the subsurface grid point coordinates.

Claims

1. A sensitivity-based full-zone apparent resistivity imaging method for ground-based frequency-domain electromagnetic, characterized in that: The method comprises the following steps: Step one, uniform half-space is gridded, a space rectangular coordinate system is established with the emission source as the origin, and the grid is divided, wherein the length of the x, y and z directions of the underground divided region is obtained by analyzing the curve of the magnetic field response with the length of the underground divided region, and taking the length of the region where the magnetic field response tends to be stable; the length of the grid is obtained by calculating the magnetic field response curve under different grid lengths, and taking the length of the grid where the magnetic field response tends to be stable as the length of the underground grid division; Step two, the electric field response value of each grid in the underground is calculated based on the magnetic vector potential function: where A x and A z are the x and z components of the geomagnetic vector potential, respectively; I1is the electric dipole moment, J0(λρ) is the zeroth order Bessel function, λρis the argument of the Bessel function, ρ = (x 2 +y 2 ) 1 / 2 , (x,y,z) is the coordinate of the observation point in the subsurface, λ is the Hankel integral variable; k0and k1are the wave numbers in the air and subsurface media, respectively, k0= ω(ε0μ0) 1 / 2 , k1= (iωμ0σ1) 1 / 2 ; μ0is the vacuum permeability, μ0= 4π x 10 -7 H / m; u j = (λ 2 +k j 2 ) 1 / 2 , j = 0, 1. The Foster identity is applied: where I0and K0are the first and second kind modified Bessel functions, respectively, r = (p 2 +z 2 ) 1 / 2 ; The Sommerfeld identity is applied: The above formula is obtained: In the above formula, N and P represent the Foster integral formula and the Sommerfeld integral formula respectively; The electric field expression of each grid in the underground is obtained from the relationship between the magnetic vector potential and the electromagnetic field: Step three, the magnetic field response of each grid at the observation point in the air is calculated based on the tensor Green function: The electric field response of each grid in the underground obtained from formula (7) is substituted into the tensor Green function: H s (r r )=∫∫∫ V' G H (r r ,r)J(r)dv' (8) In the formula, The current density J(r) is calculated from E calculated from equation (7) where r r where r and r represent the position vectors of the airborne observation point and the subsurface grid point, respectively, to obtain the magnetic field response at the airborne observation point due to each grid point. Step four, the coupling sensitivity of each grid with respect to the observation point in the air is calculated from the magnetic field response of each grid at the observation point in the air: where H z (r r is the secondary magnetic field z-component produced by the underground grid point r at the aerial receiving point r r at frequency f, and i, j, k are the grid indices in x, y, z directions, respectively. Step five, the apparent resistivity of the whole region is calculated by the observed magnetic field response at the observation point in the air, combined with the quasi-static approximation formula of the vertical magnetic induction intensity of the uniform earth space, through the bisection method iterative calculation: In the formula, J1(λρ) is the first order Bessel function, and (x, y, z) is the coordinate of the receiving point in the air; Step six, the full-zone apparent resistivity calculated from step five, combined with the coupling sensitivity distribution under the same system parameters obtained in step four, calculates the full-zone apparent resistivity ρ weighted based on the coupling sensitivity at each grid in the underground (xi,yj,zk) The calculation formula is: In the formula, D n For the nth measurement point, f Q For the Q-th frequency point, ρ (Dn,fQ) For observation point D n At this location, the frequency is f Q The apparent resistivity of the entire region obtained by iterative calculation of the magnetic field response. For in (D) n ,f Q Under system parameters, the mesh (x) i ,y j ,z k Sensitivity coefficient at (x) i ,y j ,z k () represents the coordinates of the underground grid points.

2. The full-zone apparent resistivity imaging method based on sensitivity of ground-to-air frequency domain electromagnetic according to claim 1, characterized in that: In step one, the length of the x direction of the region is 4000-4500m, the length of the y direction is 2500-3000m outside the transmitting-receiving distance, the length of the z direction is 2500-3000m, and the length of the divided grid is 30-50m.

Images