CN112668146A

CN112668146A - Improved field source position estimation method based on Euler deconvolution method practicability

Info

Publication number: CN112668146A
Application number: CN202011393146.9A
Authority: CN
Inventors: 陈青; 陈浩; 赖富强; 谭先锋; 汪佳蓓; 黄兆辉; 何加成; 黄小宇; 孙帅; 丁成艺
Original assignee: Chongqing University of Science and Technology
Current assignee: Chongqing University of Science and Technology
Priority date: 2020-12-03
Filing date: 2020-12-03
Publication date: 2021-04-16
Anticipated expiration: 2040-12-03
Also published as: CN112668146B

Abstract

Aiming at the Euler inversion calculation problem of horizontal positions and depths of geologic bodies with different burial depths, different thicknesses and different scales, the invention provides an extended inclination angle Euler deconvolution method (ETA _ Euler) independent of a construction index based on an Euler homogeneous equation and a normalized inclination angle method, so that on one hand, the divergence of a solution caused by selection of an empirical construction index is avoided, and on the other hand, the 'analytic singularity' existing in the calculation of each order derivative of the extended inclination angle (ETA) is avoided.

Description

Improved field source position estimation method based on Euler deconvolution method practicability

Technical Field

The invention relates to the field of geological exploration, in particular to a field source position estimation method based on improved practicability of an Euler deconvolution method.

Background

The Euler deconvolution is also called Euler inversion method, and is a position field inversion method capable of automatically estimating the position of a field source. The method is based on an Euler homogeneous equation, uses potential field abnormity and spatial derivatives, selects a construction index according to the shape of a field source, and automatically or semi-automatically defines the boundary and the depth of a geologic body through the solution of the Euler homogeneous equation. The method has the advantages that the prior information of the physical properties of the field source is not needed, the specific position of the geologic body can be quickly and effectively calculated by determining the construction index related to the property of the field source in advance, and the method is particularly suitable for analyzing and explaining large-area potential field data. In recent years, the Euler deconvolution method has attracted great attention due to its simplicity and rapidity. However, in practical application, due to the influence of the complexity of the geologic body and the superposition between field sources, the method has the problems of inversion result divergence, false solution, construction index determination under complex conditions and the like, and practical application of the method is limited.

(1) The limitation of source data. The conventional Euler deconvolution method is to directly use the original position field anomaly and the derivative thereof to carry out inversion calculation, and then to constrain the more divergent inversion result. The method has the advantages that: under the condition of less prior information, the boundary of the geologic body can be automatically or semi-automatically and quickly defined, and the depth estimation can be carried out on the latency source. The disadvantages are that: on one hand, in practical application, inversion calculation needs to be performed by using a construction index obtained empirically, which often results in inaccurate results or divergence of solutions; on the other hand, the original gravity anomaly field contains comprehensive information generated by a plurality of density non-uniform body field sources on the earth surface and underground, the interference of background field information exists, and deep and shallow source anomalies cannot be effectively distinguished, so that the inversion result is often dispersed, and the phenomenon of unconvergence of solution exists in the complex inversion of a plurality of geologic body models.

(2) And (4) constructing an index selection problem. The determination of the structural index is crucial in the solution process of the Euler deconvolution method, and needs to be selected according to the field source shape or the prior knowledge about the abnormal property, while the Euler result inverted by using the empirical structural index is often inaccurate or the solution diverges. Based on the derivation of the Euler homogeneous equation, the scholars propose that the solution of Euler deconvolution can be carried out by utilizing data sources such as a gravity high-order derivative, an analytic signal, a Tilt gradient, a normalized inclination angle and a Theta diagram. Has the advantages that: the Euler inversion calculation is carried out by utilizing the data source, and the convergence of the inversion result is improved to a certain extent due to the fact that shallow source abnormity is well highlighted and the influence of a background field is eliminated or weakened. Particularly, the Euler deconvolution based on Tilt gradient, normalized inclination angle and Theta diagram can deduce field source boundary and depth distribution without constructing indexes, thereby avoiding divergence and instability of Euler inversion solution caused by selection of the constructing indexes. The disadvantages are that: on one hand, the structural index still needs to be obtained by using experience based on the Euler deconvolution of the gravity higher-order derivative and the analytic signal, and the inversion result is still inaccurate or the solution is still diverged; on the other hand, based on the Tilt gradient, the normalized Tilt angle and the euler deconvolution of Theta diagram, although the problem of selecting the formation index is avoided, in the process of obtaining each order derivative, because the value is above the denominator of the horizontal total gradient (THDR), and the value may approach to 0 in the regional background field, the value may cause the "analytic singularity" of each order derivative.

(3) The problem of euler false solution. In practical application, due to the complexity of the geologic body and the superposition effect between field sources, the euler inversion result is easy to generate divergent and false solutions (pseudo solutions), so that the distribution of the field sources cannot be effectively identified. The main methods currently used to eliminate spurious solutions (pseudo solutions) are: estimating a depth standard deviation, a distance constraint evaluation criterion and a concentration constraint evaluation criterion, filtering by using a low-pass filter and the like. The method has the advantages that the distribution of the divergent solution can be eliminated to a certain degree; the disadvantages are that: firstly, the quality or vergence of the overall solution is difficult to be accurately evaluated, and secondly, under the influence of mutual superposition of field source anomalies, solutions of different depths cannot be effectively distinguished, so that a larger deviation exists in an inversion depth result.

Disclosure of Invention

The invention aims to solve the technical problem that in the prior art, the Euler result obtained by utilizing empirical structure index inversion is often inaccurate or the solution is diverged, and is influenced by the complexity of a geologic body and the superposition between field source bodies, so that the Euler inversion result is easy to generate divergence and false solution (pseudo solution), and the distribution of the field sources cannot be effectively identified.

The invention provides a field source position estimation method based on Euler deconvolution practicability improvement, which comprises the following steps:

s1, constructing a multi-scale gravity anomaly combination model, wherein the model comprises gravity anomalies of geologic bodies with different burial depths, different thicknesses and different scales;

s2, adding 2% of Gaussian noise into the gravity anomaly combined model, and performing Euler inversion;

s3, solving first and second derivatives of the combined model gravity anomaly in the x, y and z directions to obtain Vzx, Vzy, Vzz, Vzxx, Vzxy, Vzyy, Vzzx, Vzzy and Vzzz;

vzx is the first derivative of gravity anomaly in the x-direction, i.e.

Vzy is the first derivative of gravity anomaly in the y direction, i.e.

Vzz is the first derivative of gravity anomaly in the z-direction, i.e.

Vzxx is the first derivative of Vzx in the x-direction, i.e.

Vzxy is the first derivative of Vzx in the y-direction, i.e.

Vzyy is Vzy first derivative in the y direction, i.e.

Vzzx is the first derivative of Vzz in the x-direction, i.e.

Vzzy is the first derivative of Vzz in the y-direction, i.e.

Vzzz is the first derivative of Vzz in the z-direction, i.e.

S4, filtering each order derivative by upward continuation to smooth noise influence;

s5 applies the following equations, with the derivatives per order, per direction after the filtering process, to perform ETA and its derivative in the x direction kx, derivative in the y direction ky and derivative in the z direction kz:

the derivatives of ETA along the x-, y-, and z-directions were calculated using the following equations:

f represents total field strength, A represents analytic signal amplitude;

s6, the horizontal total gradient inclination angle TAHG is obtained by adopting the following formula:

THDR represents the horizontal overall gradient;

s7, taking the peak value of TAHG as a constraint, selecting the corresponding peak value of the derivative kx, ky and kz in each direction of ETA, and recording the peak value as kx _ peak, ky _ peak and kz _ peak;

s8, calculating by using kx _ peak, ky _ peak and kz _ peak according to the following formulas to obtain the horizontal position of the field source position and the burial depth parameters x0, y0 and z0 under the current sliding window; sliding the sub-window according to a preset step length, and repeatedly adopting the following formula to calculate until the whole area is covered;

s9, acquiring the boundary position distribution of the field source body and calculating the average burial depth.

The invention has the beneficial effects that:

(1) aiming at the Euler inversion calculation problem of horizontal positions and depths of geologic bodies with different burial depths, different thicknesses and different scales, the Euler inversion calculation method based on the Euler homogeneous equation and the normalized inclination angle method provides an extended inclination angle Euler deconvolution method (ETA _ Euler) independent of a construction index, on one hand, divergence of solutions caused by selection of empirical construction indexes is avoided, on the other hand, compared with other Euler deconvolution methods based on different data, the Euler inversion calculation method based on the structure index avoids 'analytic singularities' existing in the calculation of derivatives in all directions, can obtain more stable and continuous calculation results, and solves the divergence problem of inversion results. Meanwhile, the method does not depend on selection of an empirical construction index, has better practicability and is convenient to popularize and use in large-scale gravity magnetic field data processing.

(2) Because the accuracy of the horizontal position and the depth position of the solution estimated by the current main method for eliminating the false solution is not high, the method adopts a peak value constraint method based on the horizontal total gradient inclination angle (TAHG) to constrain the data points of Euler deconvolution, and improves the identification precision of geologic bodies with different burial depths, different thicknesses and different scales.

(3) Compared with an unconstrained Euler inversion result, the method adopts a horizontal total gradient inclination angle peak value constraint method, effectively constrains inversion data points, and enables the inversion result to be more convergent; meanwhile, the inversion result under the peak value constraint reduces the influence of the mutual superposition of the field source abnormity, so that the depth information of the deep field source is better reflected.

Drawings

FIG. 1 is a technical scheme of the present invention.

FIG. 2 is a schematic diagram of a combined model of adding 2% Gaussian noise

FIG. 3 is a comparison of Euler deconvolution processing results of a combined model with 2% Gaussian noise

Wherein, fig. 2(a) is a noisy gravity anomaly result of three cube models for testing the practical effect of the method of the present invention; the burial depths of the

models

1, 2 and 3 are 3km, 5km and 7km respectively. Fig. 2(b) shows the horizontal total gradient Tilt Angle (TAHG) results to define peak positions. FIG. 3(a) is the result of a conventional Euler inversion; (b) the vertical second derivative Euler inversion result is obtained; (c) the horizontal total gradient Europe inversion result is obtained; (d) analyzing the Euler inversion result of the signal; (e) is Tilt _ Euler inversion result; (f) is the TDX _ Euler inversion result; (g) the ETA _ Euler inversion result is obtained; (h) the results are the ETA _ Euler inversion under the TAHG peak constraint.

Detailed Description

The formula derivation process in the present invention is explained below.

(1) Basic principle of the extended Angle of inclination (ETA)

Tilt angle (Tilt deviation) is a method proposed for identifying field source boundaries of different burial depths, which uses the arctangent angle of the absolute value of the vertical gradient (VDR) of the total field strength f to the horizontal total gradient (THDR), defined as:

in formula (1):

wherein

And

the first derivatives of the total field strength f in the x, y and z directions, respectively. According to the nature of the arctan function, the tilt angle varies over a range of (- π/2, π/2), positive in the field source, negative on the periphery, and zero at the boundary. For a deep field source, the ratio of the vertical derivative and the horizontal derivative of the deep field source can still be large under the condition that the vertical derivative and the horizontal derivative of the deep field source are small, and therefore the inclination angle inversion result is less influenced by the buried depth. However, under regional background fields, the horizontal total gradient (THDR) at the tilt denominator may approach 0, resulting in tilt angle "analytic singularities".

The normalized slant angle method (TDX) is better able to identify the boundaries of a geologic volume than the slant angle method. TDX can be defined as:

the maxima of the TDX determine the spatial location of the field source boundary. TDX highlights the characteristic of gravity anomaly, reduces redundant information, and enables generated images to be visual, clear and easy to identify.

The derivatives of TDX along the x-, y-, and z-directions can be written as:

in the formula (4-6), the horizontal total gradient

Present in the denominator, may approach zero in the background field, which may lead to "resolving singularities". Therefore, the TDX equation is corrected by the method, and the expression can be written as follows:

the above equation is called extended tilt angle method (ETA). The derivatives of ETA along the x-, y-, and z-directions can be written as:

wherein

Solving for the derivatives in each direction of the improved tilt angleIn the process, only the analytic signal amplitude a is located at the denominator. However, even in the case where a tends to 0, there is no "analytic singularity" in each directional derivative formula (9-10) thereof, and the relationship is as follows:

thus, k can still be calculated_x ^’，k_y ^’，k_z ^’The value of (c).

(2) Basic principle of extended tilt angle Euler deconvolution (ETA _ Euler)

The basic expression of Euler deconvolution based on the Euler homogeneous equation is as follows:

wherein f is field source potential field anomaly; x, y and z are coordinates of the observation points; x is the number of₀、y₀、z₀Is a field source coordinate; n is a structural index (N ═ 1, 2, 3 … …), which is related to the geometric configuration of the field source, and is the decay rate of the field source abnormal intensity with the depth; b is called the region field or background field. Therefore, the method is based on Euler homogeneous equation, and uses the potential field anomaly and the spatial derivative thereof to determine the position and the depth of the anomalous field source by combining with the specific 'structure index' of the geologic body. However, because the field source types of the actual geological structure are complex and various, the accuracy and stability of the field source depth inversion solution are directly influenced by the correctness of the structure index value selection.

The derivatives of equation (12) along the x-, y-, and z-directions can be written as:

multiplying equation (15) by

(13) Multiplication by

Subtracting the former from the latter to obtain:

likewise, multiply equation (15) by

Equation (14) multiplied by

Then subtracting the latter from the former, we can get:

the formulas (15) are respectively multiplied by

And

then subtracting the latter from the former, we can get:

adding the equations (16-18) gives:

dividing equation (19) by

Obtaining:

combining the formulas (8-10), the above formula is known as follows:

we call the above equation ETA _ Euler deconvolution. The method has the advantages that: firstly, the method does not depend on the construction index N, so that the problem of divergence of a solution caused by improper selection of the construction index is effectively avoided; secondly, because the expression of the normalized inclination angle is improved, only the analytic signal amplitude A is positioned at the denominator, but even under the condition that A tends to be 0, the analytic singularities do not exist in all directional derivative formulas.

(3) Horizontal total gradient tilt angle peak value constraint method

In order to improve the convergence of the inversion result, the invention provides that the inversion data points are constrained by a horizontal total gradient inclination angle (TAHG) peak value, and the expression is as follows:

due to the arctangent characteristic, the TAHG transformation range is also (-pi/2, pi/2). The method not only effectively balances information from shallow and deep sources, but also maximizes at the field source boundary. Therefore, the data points of the euler deconvolution are constrained by the maximum range of TAHG, thereby improving the convergence of the inversion result.

The operation of the present invention will be described below by way of an example.

Example 1

As shown in fig. 1, the specific process of the present invention in the noise model test includes,

(1) constructing gravity anomaly combination models with different scales (depth, thickness and size);

(2) because the actual position field data contains certain noise, in order to test the applicability of the method, 2% of Gaussian noise is added into the gravity anomaly combination model, Euler inversion is carried out, and the stability of an inversion result is analyzed;

(3) calculating derivatives of each order and each direction: solving first and second derivatives of the combined model gravity anomaly in the x, y and z directions, namely obtaining Vzx, Vzy, Vzz, Vzxx, Vzxy, Vzyy, Vzzx, Vzzy and Vzzz;

(4) and (3) filtering treatment: because the gravity anomaly of the combined model contains noise, each order of derivative is interfered by the noise, particularly the second order derivative, and each order of derivative is subjected to filtering processing by upward continuation to smooth the influence of the noise;

(5) ETA and its directional derivative calculation: ETA and the derivatives (kx, ky, kz) in the directions of x, y and z are calculated by using the derivatives in each order and each direction after the filtering treatment, and the calculation is shown in formulas (7), (8), (9) and (10);

(6) horizontal total gradient tilt angle peak calculation: calculating a horizontal total gradient inclination angle (TAHG) by using the first-order derivatives in all directions after the filtering treatment, see formula (22); because TAHG enhances the information of deep weak anomaly, the information of deep and shallow anomaly is effectively highlighted, and the peak position of each geologic body boundary can be obtained by utilizing the maximum value of the contour line;

(7) calculation of the peak value of the derivative in each direction of ETA: selecting the peak values of kx, ky and kz corresponding to the peak values of TAHG as constraints, and recording the peak values as kx _ peak, ky _ peak and kz _ peak;

(8) ETA _ Euler calculation: substituting kx _ peak, ky _ peak and kz _ peak into a formula (21), setting a sliding window, and calculating the results of the position parameters x0, y0 and z0 of the field source body under the current sliding window, namely the horizontal position and the burial depth; sliding the sub-window according to a certain step length, and repeating ETA _ Euler calculation until the whole area is covered;

(9) and acquiring the boundary position distribution of the field source body and calculating the average burial depth.

FIG. 2(a) is a noisy gravity anomaly result for three cube models to test the utility of the method of the present invention; the burial depths of the

models

1, 2 and 3 are 3km, 5km and 7km respectively. Fig. 2(b) shows the horizontal total gradient Tilt Angle (TAHG) results to define peak positions. FIG. 3 is a graph of Euler inversion results processed by the present invention and their comparison with other Euler inversion results.

Fig. 3(g) is the result of inversion directly with ETA _ Euler, and the convergence of the solution is significantly better than the inversion based on the first 4 data, compared to the conventional Euler deconvolution result (fig. 3a), and the inversion based on the vertical second derivative europe (fig. 3b), the horizontal total gradient (fig. 3c), the analytic signal (fig. 3d), the Tilt _ Euler (fig. 3e), and the TDX _ Euler (fig. 3 f). However, the Tilt _ Euler and TDX _ Euler results are similar to the ETA _ Euler method results, and both reflect better convergence of the solution, but for the geological body No. 3 with the maximum buried depth (7km), the boundary contours of the two are not clear enough, while the boundary reflected by the ETA _ Euler method is more convergent and clear. Furthermore, the results of Tilt _ Euler and TDX _ Euler appear to be spurious solutions greater than 25km from the depth of the solution, compared to the depth results of the ETA _ Euler method, which are closer to the true values, but for geological volumes No. 1 (3km) and No. 2 (5km) with shallow burial depths, the depth results are not sufficiently well differentiated due to the influence of the anomalous superposition of the field sources. The inversion results of the ETA Euler method using the TAHG peak constraint are shown in fig. 3 (h). It can be seen that the depths of the geologic body with three different burial depths are clearly reflected, and the depth values are respectively 3.12 +/-0.47, 4.91 +/-0.29 and 6.99 +/-0.44 km, and the errors are respectively 4.07%, 1.80% and 0.10%.

TABLE 1 comparative analysis of depth results of various models calculated on the basis of different data

In summary, the model depth position calculated by the method is better in coincidence with a theoretical value, and compared with the traditional method, the obtained boundary position is obviously more convergent and more accurate; the accuracy of the identified depth position is higher, and particularly, the inversion result of the ETA _ Euler method utilizing the TAHG peak value constraint can clearly distinguish geologic body depths with different burial depths. The euler inversion depth results of different data bases are compared and shown in table 1.

The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A field source position estimation method based on improved Euler deconvolution practicability comprises the following steps: