A robust high-resolution geological boundary identification method based on potential field data

Abstract

The present invention relates to a robust high-resolution geological body boundary recognition method based on potential field data, belonging to the field of geophysical exploration technology. The method comprises the following steps: denoising the potential field data and using the denoised data as input for edge detection. f 1. Calculation f 1, we get f 2. Calculation f 2 Harris filter field, we get f 3. Yes f 3 Perform two-dimensional Gaussian filtering to eliminate possible singular points and obtain f 4. Yes f 4 Further ILTHG filtering is performed to obtain f 5. Extraction f The maximum point position set of 5 is used as the final geological body boundary recognition result. The invention has the advantages of strong anti-interference ability, can produce high-resolution geological body edge detection results, and is not easy to produce false boundaries. To a certain extent, it solves the trade-off problem between noise resistance and resolution faced by the existing edge detection method, and the geological body boundary recognition effect is better than the traditional edge detection method.

Description

Technical Field

[0001] The invention relates to a robust high-resolution geological body boundary recognition method based on potential field data, belonging to the technical field of geophysical exploration. Background Art

[0002] Potential field edge detection technology is a long-standing and widely used method for identifying geological boundaries. This technology can highlight information that is not easily perceived in the original potential field data and plays an important role in structural interpretation at various scales, such as ore body detection, fault identification, delineation of different geological bodies, and division of regional structural frameworks. Based on past research, potential field edge detection can be mainly divided into three categories: (1) derivative-based, (2) phase or ratio-based, and (3) statistical / sliding window-based. Among them, horizontal derivatives and vertical derivatives (VDR method) are the basis of derivative methods. For example, the earlier proposed total horizontal derivative (THDR method) is composed of horizontal derivatives in different directions. The tilt angle (TA method) is a typical phase or ratio-based method. Recently, the ILTHG method and the sliding window-based NHF method have also been proposed.

[0003] The study found that hybrid edge detection filters that use both sequential combinations and ratios have higher horizontal or vertical resolution than traditional single edge detectors based on derivatives, ratios, or sliding windows, and are less likely to produce artifacts. However, hybrid filters with high-order derivatives, such as the ILTHG method, are often more susceptible to noise interference. In this case, the effect is not as good as the NHF filter based on statistical time windows. In general, although the potential field edge detection technology is constantly developing, most of the current filters have the following limitations: (1) It is difficult to balance the response to field sources of different depths, sizes, and intensities; (2) Low-order edge detectors have limited resolution, and their resistance to noise decreases rapidly with the increase of the derivative order, that is, they face a compromise between resolution and noise resistance. Summary of the invention

[0004] In view of the above-mentioned defects in the prior art, the present invention proposes a new robust high-resolution geological body boundary recognition method based on potential field data, which can improve the anti-interference ability, produce high-resolution edge detection results, and is not prone to produce false boundaries.

[0005] The present invention is implemented by adopting the following technical solution: The present invention provides a robust high-resolution geological body boundary recognition method based on potential field data, comprising the following steps:

[0006] Step 1: Determine whether the signal-to-noise ratio of the potential field data meets the requirements. If the signal-to-noise ratio is too low, perform upward extension or other denoising processing on it, and use the denoised potential field data as the input data for edge detection. f 1;

[0007] Step 2: Calculate input data f The vertical first derivative of 1 is:

[0008] S1: The calculation formula is as follows:

[0009] (1)

[0010] This step is generally performed in the frequency-wavenumber domain, and then converted to the spatial domain through inverse Fourier transform to obtain f 2;

[0011] Step 3: f 2 Perform normalized Harris filtering:

[0012] S2: In any 3×3 window, calculate f 2 Along x , y The derivative matrix M of the direction is used to extract f The local structural characteristics of 2 are as follows:

[0013] (2)

[0014] in, n Corresponding to small displacement f 2 the number of neighbors along rows, columns, and diagonals, and They are f 2 in x Direction and y Directional derivative of the square;

[0015] S3: Calculate the response function of Harris filter. The specific formula is:

[0016] (3)

[0017] Where Det(M) represents the determinant of M, represents the trace of M, is a parameter for adjusting sharpness or smoothness. Using the upper envelope of NHF to normalize itself, we get f 3;

[0018] Step 4: f 3 Perform two-dimensional Gaussian smoothing to further eliminate potential singular points. f 3 If there is no singular point, this step can be ignored;

[0019] S4: Two-dimensional Gaussian smoothing is specifically expressed as

[0020] (4)

[0021] Among them, Gaussian() represents the two-dimensional Gaussian distribution function, so we get f 4;

[0022] Step 5: f 4. Perform ILTHG filtering, as follows:

[0023] S5: Calculate separately f 4 in x direction, y Direction and z Directional derivative of direction , and ;

[0024] S6: , and Substitute into the ILTHG filter formula:

[0025] (5)

[0026] Among them, exp represents the exponential function with the natural constant e as the base, is the power index, which is generally between 1 and 10;

[0027] Step 6: Calculate formula (5) to obtain f 5. Extraction f The maximum point position collection of 5 is taken as the final geological body boundary recognition result.

[0028] The beneficial effects of the present invention are as follows: by adopting the robust high-resolution geological body boundary recognition method based on potential field data described in the present invention, by integrating the advantages of multiple edge detection filters, a high-resolution geological body boundary recognition result can be obtained, and false boundaries are not easily generated. It is still robust under a certain noise intensity, and the problem that the existing method cannot take into account both high resolution and strong noise resistance. The noise resistance and boundary resolution are better than those of traditional methods, and the boundary recognition results are reliable. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 is a flow chart of the present invention;

[0030] Figure 2 It is a 3D theoretical model diagram;

[0031] Figure 3 Gravity anomaly maps generated for the model;

[0032] Figure 4 (a) is the result of geological boundary recognition by VDR method in the absence of noise;

[0033] Figure 4 (b) is the result of geological boundary identification using the THDR method without noise;

[0034] Figure 4 (c) is the result of TA method geological boundary identification in the noise-free case;

[0035] Figure 4 (d) is the result of geological boundary identification using the NHF method without noise;

[0036] Figure 4 (e) is the result of geological boundary identification by ILTHG method in the noise-free case;

[0037] Figure 4 (f) is a diagram showing the geological boundary identification result of the method of the present invention in the absence of noise;

[0038] Figure 5 (a) is the random noise image to be added;

[0039] Figure 5 (b) is the gravity anomaly map after adding random noise;

[0040] Figure 6 (a) is the result of VDR method geological boundary recognition after adding random noise;

[0041] Figure 6 (b) is the result of geological boundary identification using the THDR method after adding random noise;

[0042] Figure 6 (c) is the result of TA method geological boundary identification after adding random noise;

[0043] Figure 6 (d) is the result of geological boundary identification using the NHF method after adding random noise;

[0044] Figure 6 (e) is the result of geological boundary identification using the ILTHG method after adding random noise;

[0045] Figure 6 (f) is a geological boundary identification result diagram of the method of the present invention after adding random noise; DETAILED DESCRIPTION

[0046] In order to make the purpose and technical solution of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings. Figure 1 As shown, the following steps are included:

[0047] Step 1: Determine whether the signal-to-noise ratio of the potential field data meets the requirements. If the signal-to-noise ratio is too low, perform upward extension or other denoising processing on it, and use the denoised potential field data as the input data for edge detection. f 1;

[0048] Step 2: Calculate input data f The vertical first derivative of 1 is:

[0049] S1: The calculation formula is as follows:

[0050] (1)

[0051] This step is generally performed in the frequency-wavenumber domain, and then converted to the spatial domain through inverse Fourier transform to obtain f 2;

[0052] Step 3: f 2 Perform normalized Harris filtering:

[0053] S2: In any 3×3 window, calculate f 2 Along x , y The derivative matrix M of the direction is used to extract f The local structural characteristics of 2 are as follows:

[0054] (2)

[0055] in, n Corresponding to small displacement f 2 the number of neighbors along rows, columns, and diagonals, and They are f 2 in x Direction and y Directional derivative of the square;

[0056] S3: Calculate the response function of Harris filter. The specific formula is:

[0057] (3)

[0058] Where Det(M) represents the determinant of M, represents the trace of M, is a parameter for adjusting sharpness or smoothness. Using the upper envelope of NHF to normalize itself, we get f 3;

[0059] Step 4: f 3 Perform two-dimensional Gaussian smoothing to further eliminate potential singular points. f 3 If there is no singular point, this step can be ignored;

[0060] S4: Two-dimensional Gaussian smoothing is specifically expressed as

[0061] (4)

[0062] Among them, Gaussian() represents the two-dimensional Gaussian distribution function, so we get f 4;

[0063] Step 5: f 4. Perform ILTHG filtering, as follows:

[0064] S5: Calculate separately f 4 in x direction, y Direction and z Directional derivative of direction , and ;

[0065] S6: , and Substitute into the ILTHG filter formula:

[0066] (5)

[0067] Among them, exp represents the exponential function with the natural constant e as the base, is the power index, which is generally between 1 and 10;

[0068] Step 6: Calculate formula (5) to obtain f 5. Extraction f The maximum point position collection of 5 is taken as the final geological body boundary recognition result.

[0069] The model test of the present invention is explained and illustrated below in conjunction with specific implementation methods.

[0070] Embodiment 1:

[0071] In order to illustrate the implementation idea and process of this method and prove the effectiveness of the method, a three-dimensional model is established for testing and compared with the results of traditional edge detection methods.

[0072] S1: Establish Figure 2 The theoretical three-dimensional model shown. The model consists of three prisms with positive density and three prisms with negative density. Prism ① is 15 km long, 50 km wide, 5 km high, and has a density difference of -0.3 g / cm 3 ; Prism ② is 15 km long, 50 km wide, 5 km high, and has a density difference of 0.35 g / cm3 ; Prism ③ is 20 km long, 20 km wide, 5 km high, and has a density difference of 0.4 g / cm 3 ; Prism ④ is 80 km long, 20 km wide, 5 km high, and has a density difference of 0.3 g / cm 3 ; Prism ⑤ is 60 km long, 60 km wide, 5 km high, and has a density difference of -0.5 g / cm 3 ; Prism ⑥ is 20 km long, 20 km wide, 5 km high, and has a density difference of -0.4 g / cm 3 .

[0073] S2: Within the scope of this model, gravity observation points are arranged along the horizontal surface (z=0 km) with a spacing of 2 km between the observation points.

[0074] S3: Calculate the gravity anomaly produced by this model f 1. Such as Figure 3 shown.

[0075] S4: Calculate input data f 1, we get f 2.

[0076] S5: Yes f 2 Perform normalized Harris filtering and get f 3.

[0077] S6: Yes f 3 Perform two-dimensional Gaussian smoothing and ILTHG filtering to obtain f 4 and f 5.

[0078] S7: Extraction f The maximum point position collection of 5 is taken as the final geological body boundary recognition result.

[0079] In order to illustrate the effect of the method of the present invention, the geological body boundary recognition results of the present invention are compared with the results of existing methods (including VDR method, THDR method, TA method, NHF method and ILTHG method). Figure 4 (a~e) are the geological boundary identification results of VDR method, THDR method, TA method, NHF method and ILTHG method, respectively. Figure 4(f) is the geological boundary identification result of the present invention. It can be clearly seen from the figure that the VDR method corresponds well to the real block boundary position near the zero value point. However, not all zero value points indicate the real boundary position. Therefore, VDR has ambiguity and false boundaries. The maximum value of the THDR method indicates the position of the boundary, but has a lower resolution and is less sensitive to deep source small targets (prisms ③). The performance of the TA method is roughly equivalent to that of the VDR method, and it inherits the shortcomings of the VDR method, that is, there are obvious excess zero-crossing points, which leads to the generation of false boundaries. The NHF method has a stronger resolution than the THDR method, but it also cannot detect small blocks in the deep. In contrast, the ILTHG method and the method of the present invention are better able to balance the amplitude intensity of field sources of different sizes and depths, and all six block boundaries are detected. Compared with the ILTHG method, the edge line of the method of the present invention is finer and the block boundary performance is sharper, indicating a higher resolution.

[0080] Embodiment 2:

[0081] In order to further illustrate the noise resistance and robustness of this method, the gravity anomaly generated by the model in Example 1 is added as follows Figure 5 The random noise shown in (a) was tested and compared with the results of existing methods (including VDR method, THDR method, TA method, NHF method and ILTHG method). The gravity anomaly after adding random noise is shown in Figure 5 The gravity anomaly after adding random noise is extended upward by 4 km and used as the input data for edge detection. f 1. Afterwards, the specific steps of the second embodiment are the same as S4 to S7 of the first embodiment. Figure 6 (a~e) are the geological boundary recognition results of VDR method, THDR method, TA method, NHF method and ILTHG method after adding random noise. Figure 6 (f) is the geological boundary recognition result after adding random noise in the present invention.

[0082] It can be seen that the results of the VDR method and the THDR method are relatively fuzzy, and there is a boundary misalignment at the position where prism ③ and prism ⑤ meet. The results of the TA method are similar to those of the above two methods, but it seems to be more affected by noise. The NHF method is less affected by noise, but its resolution is still limited. The boundary resolution of the ILTHG method is relatively high, but there are false trivial spots, and the large noise level causes the instability of the ILTHG method. Compared with the ILTHG method, the method of the present invention greatly reduces the possibility of detecting discrete false points as boundaries. Compared with the NHF method, the method of the present invention significantly improves the horizontal resolution of edge detection, and the positioning of the boundary is clearer and more reliable. It can be seen that the method of the present invention is superior to existing edge detection methods.

[0083] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. All equivalent modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the patent scope of the present invention.

Claims

1. A robust high-resolution geological body boundary recognition method based on potential field data, characterized in that: The steps include: Step 1: Determine whether the signal-to-noise ratio of the potential field data meets the requirements. If the signal-to-noise ratio is too low, perform upward extension of the gravity anomaly or other denoising processing to suppress random noise, and use the denoised potential field data as the input data f1 for edge detection; Step 2: Calculate the vertical first-order derivative of the input data f1: S1: The calculation formula is as follows: This step is performed in the frequency-wavenumber domain and then converted to the spatial domain through inverse Fourier transform to obtain f2; Step 3: Normalize f2 using Harris filtering: S2: In any 3×3 window, the derivative matrix M of f2 along the x and y directions is calculated to extract the local structural features of f2. The specific formula is: where n corresponds to the number of neighbors of f2 along rows, columns, and diagonals within a small displacement, and are the derivatives of f2 in the x and y directions respectively; S3: Calculate the response function of Harris filter. The specific formula is: NHF=Det(M)+μTr 2 (M) (3) Where Det(M) represents the determinant of M, Tr(M) represents the trace of M, and μ is a parameter for adjusting sharpness or smoothness. The upper envelope of NHF is normalized to itself to obtain f3. Step 4: Perform two-dimensional Gaussian smoothing on f3 to further eliminate potential singular points; if f3 no longer has singular points, directly use f3 as f4, ignore this step, and jump to step 5; S4: Two-dimensional Gaussian smoothing is specifically expressed as f4=Gaussian(f3) (4) Wherein, Gaussian() represents a two-dimensional Gaussian distribution function, thereby obtaining f4; Step 5: Perform ILTHG filtering on f4, as follows: S5: Calculate the directional derivatives of f4 in the x, y and z directions respectively and S6: and Substitute into the ILTHG filter formula: Among them, exp represents the exponential function with the natural constant e as the base, and α is the power exponent, which is between 1 and 10; Step 6: Calculate formula (5) to obtain f5, and extract the maximum value point position set of f5 as the final geological body boundary recognition result.

Images