WO2016138874A1 - 位场构造格架自动提取方法 - Google Patents

位场构造格架自动提取方法 Download PDF

Info

Publication number
WO2016138874A1
WO2016138874A1 PCT/CN2016/075626 CN2016075626W WO2016138874A1 WO 2016138874 A1 WO2016138874 A1 WO 2016138874A1 CN 2016075626 W CN2016075626 W CN 2016075626W WO 2016138874 A1 WO2016138874 A1 WO 2016138874A1
Authority
WO
WIPO (PCT)
Prior art keywords
magnetic
data
scale
anomaly
grid
Prior art date
Application number
PCT/CN2016/075626
Other languages
English (en)
French (fr)
Inventor
曹殿华
Original Assignee
中国地质科学院矿产资源研究所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from CN201510096981.9A external-priority patent/CN104658037B/zh
Priority claimed from CN201510303316.2A external-priority patent/CN104965232B/zh
Application filed by 中国地质科学院矿产资源研究所 filed Critical 中国地质科学院矿产资源研究所
Priority to RU2017133223A priority Critical patent/RU2664488C1/ru
Priority to CA2978500A priority patent/CA2978500C/en
Priority to US15/555,153 priority patent/US10884161B2/en
Priority to AU2016228027A priority patent/AU2016228027B2/en
Publication of WO2016138874A1 publication Critical patent/WO2016138874A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V11/00—Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00
    • G01V11/002—Details, e.g. power supply systems for logging instruments, transmitting or recording data, specially adapted for well logging, also if the prospecting method is irrelevant
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
    • G01V3/081—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices the magnetic field is produced by the objects or geological structures
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02—Details
    • G01V7/06—Analysis or interpretation of gravimetric records
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20036—Morphological image processing
    • G06T2207/20044—Skeletonization; Medial axis transform
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20064—Wavelet transform [DWT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20172—Image enhancement details
    • G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation

Definitions

  • the invention relates to an automatic extraction method for a gravity bit field and a magnetic field field structure frame, which is a technology for detecting geological structures based on gravity and magnetic anomaly data. More specifically, the present invention relates to the fields of wavelet analysis, image processing, geophysics, geology, mineral exploration, etc., and the method according to the present invention can be directly applied to the fields of mineral exploration and related geological surveys.
  • gravity and magnetic methods (hereinafter also referred to as gravity and magnetic) measurement methods are economical, rapid, and can cover many difficult-to-reach landscape areas, and play an increasingly important role in metal deposit exploration and evaluation.
  • the role especially the development of high-precision aeromagnetic measurement technology, makes the geological structure detection method based on magnetic anomaly data to control the formation of deposits, which is of great significance in all stages of metal deposit exploration selection to target location.
  • gravity and magnetic exploration methods are commonly applied to the direct detection of mineralized bodies with strong magnetism or high density, as well as geological structure interpretation and inversion with strong anomalies.
  • geological structure detection based on gravity field data and magnetic field data to control the formation of deposits is of great significance in the target location of metal deposits.
  • the automatic identification and extraction methods of gravity magnetic field structure information mainly include analytical signal method, Euler deconvolution method, and multi-scale edge detection method of bit field.
  • the common problem is that it is not sensitive to directional information and cannot be obtained. Complete and accurate abnormal boundary position.
  • the invention is entitled "Position Field Multi-Directional Multi-Scale Edge Detection Method", and the patent application No. 200810006676.6 discloses a multi-directional multi-scale edge detection method for a bit field, which is enhanced by a direction wavelet transform.
  • the directional information obtains the abnormal source boundary information in different directions, realizes the automatic extraction of the structural grid, and overcomes the analytical signal method, the Euler deconvolution method, the bit field multi-scale edge detection method, etc., which is insensitive to directional information.
  • the problem is a technical solution for the rapid inversion of the shallow crust three-dimensional structure based on the bit field data.
  • the calculated edge is not a single pixel point width, and the actual geographic range corresponding to the intersection of the edge and the different direction edges is large, resulting in low accuracy of the analysis result; (2) need to be based on different The edge of the scale is artificially vectorized, and the center line of the edge is taken as the abnormal boundary to generate the three-dimensional structure of the shallow crust, resulting in low work efficiency. (3) The three-dimensional structure of the shallow crust can only display structural information of different depths. Reflecting the structural belt and two Lateral lithological changes do not indicate structural deformation and activity intensity; (4) there is no clear definition of scale, that is, no specific geophysical properties are assigned to the scale. Therefore, it is desirable to provide a method capable of extracting a bit field construction grid with high precision based on acquired bit field data.
  • TMI Total Magnetic Intensity
  • a TMI anomaly can be obtained by performing a geomagnetic normal field correction on the TMI.
  • TMI anomalies have problems such as lateral offset, morphological deformation, and positive and negative values. It is often necessary to perform a reduction to pole process on the TMI anomaly data to eliminate the effects of such factors.
  • the morphological pole converts the observed TMI anomaly into a perpendicular magnetic anomaly in the case of perpendicular magnetization, which converts the observed TMI anomaly into an anomaly that can be measured at the north magnetic pole, thereby migrating the magnetic anomaly directly above the source region, facilitating the geological interpretation of the magnetic anomaly.
  • the object of the present invention is to provide an automatic extraction method for a gravity magnetic field structure grid, which can quickly obtain geological structure information for controlling the formation of a deposit, thereby realizing the target location of the metal deposit.
  • the invention further provides an automatic extraction method of a magnetic structure grid in a low latitude region, which can not only obtain linear structures (Lineaments) but also obtain a ring structure.
  • the present invention defines the configuration thus obtained as a structural framework.
  • the present invention uses a multi-directional multi-scale edge detection method to extract a gravity structure grid and a magnetic structure grid for the measured gravity bit field data and magnetic field data, and each of the obtained morphological skeleton algorithms will be obtained.
  • the scale structure grid is refined into a single pixel width, and the texture grids of different scales are superimposed with different color gradients to display different depth structure information, and a comprehensive structure grid map is generated; corresponding to each edge point of the grid structure based on different scales
  • the structural strength values of the gradient model extraction are superimposed and displayed on different gradient colors, and the density change and the magnetic change intensity information are highlighted to generate a comprehensive structural strength grid map.
  • the gravity and magnetic anomaly information of different depths in the study area can be obtained, Characterize the tectonic grid distribution information of different depth geological structures, the density variation of different depth tectonic frameworks and the magnetic variation intensity information, realize the identification and qualitative interpretation of the geological structure controlling the formation of the deposit, and determine the potential according to the prior knowledge of the research analysis area.
  • the type of deposit and the structural properties of the controlled deposits are screened for different types of tectonic grids to achieve the target location of the metal deposit.
  • the pre-processing includes performing a polarization calculation on the magnetic data to obtain a pseudo-magnetic anomaly or performing pseudo-gravity calculation to obtain a pseudo-gravity anomaly; and pre-processing the gravity data to obtain Bouguer gravity anomaly data.
  • the pre-processed gravity data or the two-dimensional directional wavelet transform of the magnetic data is connected to the edge of the gradient in the vertical direction of the gradient to form an edge.
  • the edges calculated in different directions of the same scale are obtained, and the obtained edge is obtained as the edge of the scale, so that the multi-directional edge detection of the potential field at each scale can be realized.
  • the obtained structural grids of various scales are superimposed to generate a comprehensive structural grid map reflecting different depth information.
  • the lateral offset of the different scale edges on the map reflects the occurrence information of the structural grid.
  • edges extracted after the height of the bit field are extended to different heights correspond to the structures of different depths, and the depth is half of the height after the upward extension (see the author is Jacobsen, BH, titled A case for upward continuation as a standard separation filter for potential- Field maps, journal name Geophysics, issue number v.52no.8, time 1987), can be used to characterize the texture grids at different depths.
  • the modulus of the gradient at each edge point on the calculated edge is taken as the intensity value at the edge point in the construction grid.
  • the intensity value of the edge points on the structural grid with intensity values reflects the lithological changes of the structural belt and both sides, reflecting the deformation and activity intensity of the structure.
  • the total horizontal derivative and the analytical signal of the Tilt derivative (TDR) of the TMI abnormal data are not affected by the magnetic dip angle, and the calculation result is independent of the magnitude of the magnetic dip value, and can be directly
  • the TMI anomaly data is used to calculate the total horizontal derivative of TDR and the analytical signal without performing polarization processing.
  • the calculation result of the analytical signal method will increase the magnetic anomaly range, loss of geological structure occurrence information and structural division information, and is not sensitive to the identification of geological structures.
  • the results of the total horizontal derivative of the current oblique derivative are expressed in the form of grid images or contour maps.
  • the traditional horizontal gradient-based edge detection method cannot obtain complete and accurate source bodies of magnetic anomalies because it does not consider the directional information of the data.
  • the invention can effectively identify and establish a magnetic tectonic framework in a low latitude region by performing multi-directional edge detection based on the total horizontal derivative of the oblique derivative at multiple scales on the pre-processed TMI anomaly data. Therefore, the automatic extraction method of the structural grid according to the present invention is particularly suitable for automatic extraction of the structural grid of magnetic measurement data in low latitude regions.
  • a method for automatically extracting a bit field structure grid comprising the steps of:
  • the preprocessed bit field data is extended upward by a plurality of predetermined heights to obtain a plurality of bit field data of corresponding scales;
  • the morphological skeletal algorithm is used to refine the calculated bit field edges of each scale into single pixel widths, and obtain a structural grid map of multiple corresponding scales.
  • the method further comprises overlaying the calculated structural grid maps of the plurality of respective scales to generate a comprehensive texture grid map.
  • the method further comprises using a mode of the gradient at each edge point on each scale structure grid map as the intensity value at the edge point in the scale structure grid map, and obtaining a plurality of structural scales of the corresponding scales Frame map.
  • the method further comprises overlaying the plurality of respective scales of structural strength grid maps to generate a composite structural strength grid map.
  • the bit field data is gravity bit field data or magnetic method bit field data
  • the preprocessing further comprises performing polarization calculation on the magnetic method data to obtain a polar magnetic abnormality or performing pseudo gravity calculation to obtain a pseudo Gravity anomaly; or preprocessing the gravity data to obtain a Bouguer gravity anomaly.
  • multi-directional edge detection is performed for the bit field data of each scale, including the following steps:
  • k(x, y, z) is the Green's function
  • the wavelet function in direction ⁇ is defined as:
  • the position (x, y), scale s and orientation ⁇ , f 0 (x, y ) in the direction of the two-dimensional wavelet transform W [f 0] (x, y, s, ⁇ ) and f z (x, y) of gradient Proportional, f 0 (x, y) in the direction of the two-dimensional wavelet transform W [f 0] (x, y, s, ⁇ ) may be f z (x, y) of the gradient To characterize.
  • the corresponding angular angle in the horizontal direction of the gradient is:
  • the edge point is the point where the modulo M[f z ](x, y, ⁇ ) has a local maximum along the radial direction Af z (x, y, ⁇ ),
  • the local maximum point of the gradient mode is connected along the vertical direction of the gradient to form an edge.
  • the edges are calculated in a plurality of different directions ⁇ , and the calculated edges are summed to obtain the edges of the height.
  • the edge field is extended to a plurality of predetermined heights and the extracted edges correspond to different depth configurations, and the obtained structural grid maps of each scale are superimposed to obtain a comprehensive structural grid map reflecting different depth information.
  • the height of the height is characterized by a gradation color to form the integrated texture grid map.
  • the magnitude of the intensity value is characterized by a gradation color to form the integrated structural strength grid map.
  • a method for automatically extracting a texture grid comprising the steps of:
  • TMI anomaly data or Bouguer gravity anomaly data Grid the obtained TMI anomaly data or Bouguer gravity anomaly data, and extend the meshed TMI anomaly data or Bouguer gravity anomaly data to a plurality of predetermined heights to obtain gridded TMI anomalies of multiple scales.
  • Data or Bouguer gravity anomaly data T h , h is the height after upward extension;
  • the TMI anomaly data of each scale or the tilt derivative TDR h of the Bouguer gravity anomaly data are calculated by using the meshed TMI anomaly data or the Bouguer gravity anomaly data T h of each scale respectively;
  • the horizontal gradient-based multi-directional edge detection is performed for each scale of the meshed TMI anomaly data or the oblique derivative of the Bouguer gravity anomaly data, and the magnetic or gravity anomaly source edge of each scale is obtained;
  • the morphological skeletal algorithm is used to refine the edge of the magnetic or gravity anomaly source of each scale to a single pixel width, and a structural grid map with multiple scales is obtained.
  • the method further comprises superimposing the calculated texture grid maps of the plurality of scales to generate a comprehensive structure grid map.
  • the method further comprises: extending the TMI anomaly data or the Bouguer gravity anomaly data to the plurality of predetermined heights and extracting the edges corresponding to the different depths, and superimposing the obtained structural grid maps of the respective depths to reflect different A comprehensive structural grid diagram of cutting depth information.
  • the horizontal gradient-based multi-directional edge detection is performed on the tilted derivative of the meshed TMI abnormal data or the Bouguer gravity abnormal data for each scale, respectively, including the following steps:
  • Tilt derivative TDR h in direction ⁇ and The directional derivatives are defined as:
  • the angle of the horizontal gradient is:
  • the edges are calculated in a plurality of different directions ⁇ , and the calculated edges are summed to obtain the edge of the magnetic or gravity abnormal source body of the corresponding scale.
  • the method further comprises characterizing the horizontal gradient of each edge point on each scale structure grid map to characterize the structural burial depth at the edge point of the scale structure grid map, to obtain a plurality of scale representation structures A structural grid map of buried depth.
  • the method further comprises overlaying the plurality of scales of the constructed burial depth of the constructed grid map to generate an integrated buried depth structure grid map.
  • the method further comprises calculating the three-dimensional analytical signal AS h based on the meshed TMI abnormality data or the Bouguer gravity abnormality data T h of each scale, respectively, and obtaining the AS h value of each edge point, thereby obtaining multiple scales.
  • a structural grid diagram that characterizes the magnetic or density at the edge points.
  • the method further comprises superimposing the plurality of scales of the structural lattice patterns characterizing the edge points of the magnetic points or the density to form an integrated magnetic strength structure grid diagram or a comprehensive density strength structure grid diagram.
  • the method further comprises performing a process of removing noise from the calculated oblique derivative before performing edge detection.
  • the method is suitable for automatic extraction of magnetic structure grids in low latitude regions; preferably, the method is applicable to magnetic measurement data of a region with a magnetic inclination angle of ⁇ 30°; further preferably, the method is applicable to magnetic Magnetic measurement data for an area with an inclination of ⁇ 20°.
  • the gravity and magnetic anomaly information at different depths of the study area can be obtained, and the control of deposit formation is realized.
  • the identification and qualitative interpretation of the geological structure according to the prior knowledge of the study area to determine the potential type of deposit and control the formation properties of the deposit, and screen the different types of tectonic grids to achieve the target location of the metal deposit.
  • the method for automatically extracting the structural grid by using the magnetic method measurement data in the low latitude region solves the problem that it is difficult to accurately obtain the structural information by using the magnetic method measurement data in the low latitude region in the prior art.
  • the structural grid diagram obtained according to the present invention visually characterizes the depth of construction, the depth of burial, the primary-secondary relationship, the delivery relationship, the magnetic strength, etc. for the geological solution compared to the existing grid image or contour map. Interpretation and prospecting have important information.
  • the automatic extraction method of the magnetic structure grid in the low latitude region proposed by the invention is also suitable for the automatic extraction of the magnetic structure grid in the middle and high latitude regions, and is also suitable for the automatic extraction of the gravity field structure grid.
  • the magnetic method is used to analyze and acquire the region of the structural grid, and the accuracy of automatically extracting the structural grid is improved, and the identification and qualitative interpretation of the geological structure for controlling the formation of the deposit can be realized, according to the research area.
  • the knowledge is used to determine the type of potential deposits and to control the structural properties of the deposits.
  • the different types of structural frameworks are screened to achieve the target location of the metal deposits.
  • FIG. 1 is a flow chart of a method for automatically extracting a bit field structure grid according to a first embodiment of the present invention
  • FIG. 2 is a block diagram of a single pixel width structure according to a first example of the present invention
  • Figure 3 is a diagram showing a comprehensive structure grid according to a first example of the present invention.
  • Figure 4 is a structural strength grid diagram in accordance with a first example of the present invention.
  • Figure 5 is a comprehensive structural strength grid diagram according to a first example of the present invention.
  • FIG. 6 is a flow chart of a method for automatically extracting a magnetic structure grid in a low latitude region according to a second embodiment of the present invention
  • FIG. 7 is a block diagram of a single pixel width structure according to a second example of the present invention.
  • Figure 8 is a diagram showing a comprehensive structure grid according to a second example of the present invention.
  • Figure 9 is a structural grid diagram reflecting the burial depth according to a second example of the present invention.
  • Figure 10 is a diagram showing a comprehensive buried depth structure grid according to a second example of the present invention.
  • Figure 11 is a diagram showing a grid structure reflecting magnetic strength according to a second example of the present invention.
  • Figure 12 is a diagram showing an integrated magnetic strength structure grid according to a second example of the present invention.
  • FIG. 1 is a flow chart of a method for automatically extracting a bit field structure grid according to a first embodiment of the present invention.
  • Figure 1 As shown, the method includes the following steps:
  • Step 101 Preprocess the measured gravity bit field data or magnetic field data.
  • the pseudo-magnetic anomaly is obtained by calculating the magnetic pole data to obtain a pseudo-magnetic anomaly or pseudo-gravity calculation.
  • the gravity data is preprocessed to obtain a Bouguer gravity anomaly.
  • step 102 the pre-processed gravity bit field data or the magnetic field data is processed by using a multi-directional edge multi-directional edge detection method at multiple scales to obtain edges of multiple scales.
  • the multi-directional multi-scale edge detection method of the bit field comprises extending the pre-processed gravity bit field data or the magnetic method bit field data to a plurality of predetermined heights to obtain a plurality of corresponding scale bit field data and a bit field for each scale respectively.
  • a method for calculating multi-directional edge detection for each scale field including the following steps:
  • k(x, y, z) is the Green's function.
  • the wavelet function in direction ⁇ can be defined as:
  • the two-dimensional directional wavelet transform of f 0 (x, y) can be represented by a gradient:
  • the position (x, y), scale s and orientation ⁇ , f 0 (x, y ) in the direction of the two-dimensional wavelet transform W [f 0] (x, y, s, ⁇ ) and f z (x, y) of gradient Proportional, f 0 (x, y) in the direction of the two-dimensional wavelet transform W [f 0] (x, y, s, ⁇ ) may be f z (x, y) of the gradient To characterize.
  • the corresponding angular angle in the horizontal direction of the gradient is:
  • the edge point is the point at which the modulo M[f z ](x, y, ⁇ ) has a local maximum along the radial direction Af z (x, y, ⁇ ).
  • the local maximum point of the gradient's mode is connected along the vertical direction of the gradient to form an edge.
  • the edges of the complete coverage of the two-dimensional plane are calculated, and the calculated edges are summed to obtain the edges of the scale.
  • step 103 the edge of each scale calculated by the morphological skeletal algorithm is refined into a single pixel width, and a grid map of each scale is obtained.
  • the actual geographical extent corresponding to the intersection of edges and edges in different directions is significantly smaller than that of the prior art, so that the resulting bit field structure grid is closer to the characteristics of the actual geological map identification structure.
  • the aspect information is clear, enhances readability, and on the other hand facilitates geological interpretation.
  • step 104 the calculated scale structures are stacked to generate a comprehensive structure grid map.
  • step 105 the edge intensity is reflected by the modulus of the gradient of each edge point on each scale grid.
  • the modulo M[f z ](x, y, ⁇ ) of the gradient at each edge point is constant, which is a fixed value, and the modulus of the gradient at the edge point represents the structural intensity value at the edge point.
  • Step 106 Stacking the plurality of scale structural strength grid maps to generate a comprehensive structural strength grid map.
  • the superimposed display of different gradient colors for the intensity values of the edges or depth edges highlights the magnetic and density variation intensity information.
  • the aeromagnetic data is subjected to block-shaped pole processing, and the processed data is spliced into a bit field grid file with a grid size of 500 meters.
  • the aeromagnetic data in this example is the aeromagnetic data collected multiple times in history. Data, measuring height in the range of 800-1200 meters.
  • the skeletal algorithm is used to refine the calculated edges of each scale to obtain the structural grid map of each scale.
  • Fig. 2 is a single-pixel width structure grid of a corresponding scale extracted by the skeletal algorithm for the upper extension edge of 5000 m. It can be seen that the bit field structure grid of the single pixel point width obtained by the automatic extraction according to the present invention is closer to the actual geological map recognition structure feature, and is convenient for geological interpretation. In addition, the information on the surface is clearer, which facilitates the overlay analysis of grids of different scales.
  • Figure 3 shows the scales of the scales obtained by extending the bit field upwards by 1000, 1500, 2000, 2500, 3000, 4000, 5000, 10000, 15000, 20000, 25000, and 30000 meters.
  • a comprehensive structural grid diagram formed by color gradual superposition. The respective heights after the extension are corresponding to the depth of the corresponding source region, and the comprehensive structure grid map can be used to characterize the structural grid information at different depths of the study area.
  • the gradient from gray to black represents the upward continuation height from low to high or the depth from shallow to deep.
  • the comprehensive structural grid diagram reflects the structural information at different depths.
  • Fig. 4 is the structural strength of the mode of the gradient at the edge point of the upper 5,000-meter elongate height grid, and the structural strength value represented by the mode which reflects the gradient from gray to black gradually increases. This figure reflects the change in structural strength at a corresponding scale.
  • Figure 5 is a comprehensive structural strength grid obtained by superimposing the edge strengths obtained by extending the potential fields by 1000, 1500, 2000, 2500, 3000, 4000, 5000, 10000, 15000, 20000, 25000, and 30000 meters.
  • Figure. It can be seen from the figure that at different scales, it is reflected that the main tectonic belts show a large structural strength, and the corresponding magnetic anomaly changes obviously, indicating that the main structural belts are all magnetic anomaly, and the corresponding depth is very Large, reflecting the structural belt controlling the deep magma-mineralization.
  • the structural grid map with depth and intensity information can clearly reflect the large-depth, long-length backbone structure in the region, and the shallower, relatively short-duration secondary structures, and the mutual The delivery relationship, and thus the structural grid map obtained by the method of the present invention, can help those skilled in the art to recognize the detected regional structural pattern.
  • Figures 3 and 5 Comparing Figures 3 and 5 with the surface mapping geological structure of the region, the spatial location and extent are very consistent, illustrating the accuracy and effectiveness of the automated extraction method in accordance with the present invention.
  • Figures 3 and 5 of the present invention add information on the three-dimensional extension, strength, structure, etc., which can help identify the concealed structural belt.
  • the method according to the present invention can quickly and accurately extract aeromagnetic data and gravity data.
  • a magnetic flow measurement data is taken as an example to specifically describe a flow chart of a method for automatically extracting a magnetic structure grid in a low latitude region according to a second embodiment of the present invention.
  • the method of the present invention is not limited to magnetic measurement data, but is also applicable to the automatic extraction of construction grids of gravity measurement data.
  • FIG. 6 is a flow chart of a method for automatically extracting a magnetic structure grid in a low latitude region according to the present invention. As shown in FIG. 6, the method includes the following steps:
  • step 601 the observation value obtained by the magnetic method is preprocessed, and the geomagnetic normal field (IGRF) is corrected to obtain the total magnetic field strength (TMI) abnormal data.
  • IGRF geomagnetic normal field
  • the TMI anomaly data is meshed.
  • the grid size is taken as 1/8 to 1/4 or the minimum dot pitch of the line spacing.
  • Step 602 Upward continuation a plurality of predetermined heights on the meshed TMI abnormal data T to obtain a plurality of meshed TMI abnormal data T h of corresponding scales, where h represents the height after the upward extension.
  • Step 603 Calculate the tilt derivative TDR h of the corresponding scale TMI abnormal data by using the meshed TMI abnormal data T h of each scale respectively:
  • VDR h and THDR h are the vertical first derivative and total horizontal derivative of the meshed TMI anomaly data T h , respectively:
  • Step 604 calculating multi-directional edge detection based on the horizontal gradient for each of the scaled derivatives TDR h of each scale.
  • Tilt derivative TDR h in direction ⁇ and The directional derivatives are defined as:
  • TDR_THDR The total horizontal derivative of the oblique derivative of the TMI anomaly data
  • the magnitude of the amplitude is independent of the magnitude of the magnetic dip.
  • edge points of the magnetic source body for the height h and the direction ⁇ are modulo Along the radial direction a point with a local maximum
  • the local maximum point of the gradient mode is connected along the vertical direction of the gradient to form an edge
  • the edges are calculated in a plurality of different directions ⁇ , and the calculated edges are summed to obtain the edge of the magnetic source body of the corresponding scale.
  • the subsequent multi-directional edge detection calculation is sensitive to noise.
  • a noise reduction process such as Gaussian filtering on the TDR h data having a large noise before calculating the multi-directional edge detection.
  • the calculation result is independent of the magnitude of the magnetic dip value.
  • the edge of the magnetic anomaly source obtained by the above steps is not affected by the magnetic dip angle, ie, is not affected by the latitude of the area to be analyzed. The influence of position can accurately characterize the structural lattice in low latitudes frame.
  • step 605 the calculated edge of each scale is separately processed into a single pixel width by using a morphological skeleton algorithm to obtain a structural grid map at multiple scales.
  • the actual geographical extent corresponding to the intersection of edges and edges in different directions is significantly smaller than that of the prior art, so that the structural grid is closer to the characteristics of the actual geological mapping structure.
  • the surface information is clear, the readability is enhanced, and the geological interpretation is facilitated on the other hand.
  • Step 606 The calculated structural grid maps of each scale are superposed to generate a comprehensive structural grid map.
  • the resulting structural grid maps of various scales are superimposed to generate a comprehensive structural grid map reflecting different depth information.
  • the lateral offset of the different scale edges on the graph reflects the occurrence information of the structural grid.
  • the meshed TMI anomaly data is extended to different heights and the extracted edges correspond to different depths.
  • the depth is half of the height after the upward extension (see the author is Jacobsen, BH, titled A case for upward continuation as a standard).
  • Separation filter for potential-field maps, journal name Geophysics, issue number v.52no.8, time 1987) can be used to characterize the texture of different cutting depths.
  • the obtained structural grid maps of the respective depths are superimposed to obtain a comprehensive structural grid map reflecting different cutting depth information.
  • the burial depth is constructed by using a model of the horizontal gradient of the TDR h at each edge point on the grid, and a structural grid map of the burial depth of the plurality of scales is obtained.
  • the mode of the horizontal gradient at each edge point The size does not change and is a fixed value.
  • the magnitude of TDR h is limited to between - ⁇ /2 and + ⁇ /2. Therefore, the mode of the horizontal gradient at each edge point
  • the size of the TIM has little to do with the magnitude of the TMI.
  • This value reflects the depth of the source region. The value is inversely proportional to the depth of the burial. The larger the value, the shallower the burial depth of the source region (see author Bruno Verduzco et al., titled New). Insights into magnetic derivatives for structural mapping, the journal name The Leading Edge, issue v.23no.2, time 2004).
  • the structure of the horizontal gradient of the oblique derivative TDR h reflects the thickness of the cover layer on the top of the structural belt, while the different depths corresponding to the upward extension of the plurality of predetermined heights reflect the depth of cut of the structural belt, which represents Construct the depth at which the belt extends downward.
  • the modulus of the gradient at the edge point Representing the relative burial depth of the edge at the edge point at the edge point, and establishing a structural grid map of the relative burial depth of the reaction structures of different scales or different depths.
  • the change in the magnitude of the modulus of the gradient of the same structural zone reflects the depth of burial of different parts of the structural zone.
  • Step 608 stacking the structural grid maps of the plurality of scales to characterize the buried depth to generate a comprehensive buried depth structure grid map.
  • modulus values of the TDR h horizontal gradients for each scale or depth edge point are superimposed with different gradient colors, highlighting the burial depth variation information of the structural grids with different cutting depth ranges.
  • Step 609 calculating a three-dimensional analysis signal AS h for the meshed TMI abnormal data T h of each height:
  • Calculated magnitude of the analytical signal AS h has a strong correlation with the magnitude TMI abnormalities, but magnetic tilt value irrespective of the size, can be used to indicate the position of the source region and the magnetic field strength of the magnetic anomalies.
  • step 610 the magnetic strength of the structure is indicated by the AS h value at each edge point on each scale grid.
  • the magnitude of the three-dimensional analytical signal AS h indicates the magnetic strength of the magnetic anomaly source, and the strong magnetic structure is usually closely related to mineralization.
  • Step 611 stacking the plurality of scale magnetic strength structure grid maps to generate a comprehensive magnetic strength structure grid map.
  • the aeromagnetic data used in this example has a measuring scale of 1:25000 and a measuring height of 70-120 meters.
  • the observations obtained from the aeromagnetic survey in the study area are preprocessed, and the geomagnetic normal field (IGRF) is corrected to obtain the total magnetic field strength (TMI) anomaly data.
  • the TMI abnormal data is meshed, and the grid size is taken as the value. 10 m.
  • the gridded TMI anomaly data is separately extended upward to obtain gridded TMI anomaly data T h of multiple scales, and the elevation heights are 100, 200, 300, 400, and 500 meters, respectively.
  • Multi-directional edge detection based on horizontal gradient is performed on the calculated TDR for each of the upper extension heights.
  • FIG. 7 is a single-pixel width magnetic structure grid with a scale of 300 M in the upper limit of the gridded TMI anomaly data, and the edge is refined by the bone algorithm. It can be seen that the structure frame of the single pixel dot width obtained by the automatic extraction according to the present invention is closer to the actual geological map identification structure feature, and is convenient for geological interpretation. In addition, the information on the surface is clearer, which facilitates the overlay analysis of the grids of different scales.
  • Fig. 8 is a comprehensive structure of color gradual superposition of structural grids of various scales obtained by extending the gridded TMI anomaly data up to 100, 200, 300, 400, and 500 meters upwards and then detecting the 64 directional edges.
  • the respective heights after the extension are corresponding to the depth of the corresponding source region, and the comprehensive structure grid map can be used to characterize the structural grid information at different depths of the study area.
  • the gradient from gray to black represents the upward extension height from low to high or the depth of cut from shallow to deep.
  • the comprehensive texture grid map reflects the construction information of different cutting depths.
  • Figure 9 is a structural burial depth characterized by a gradient of the gradient at the edge of the 300 m upper echelat, and the burial depth represented by the mode that reflects the horizontal gradient from gray to black gradually decreases. This figure reflects the variation of the burial depth of the structure at the corresponding scale.
  • Figure 10 is a comprehensive buried depth structure grid diagram obtained by superimposing the gridded TMI anomaly data on the elevation heights of 100, 200, 300, 400, and 500 meters, respectively. . It can be seen from the figure that at different scales, it reflects that the burial depth of the main structural belt varies greatly. Generally, the buried deep structure is thicker, and the corresponding TMI abnormal amplitude is lower, which is a concealed structure. The method contributes to the identification of hidden structures.
  • the three-dimensional analysis signal AS h is calculated for the meshed TMI abnormality data T h of each scale.
  • the calculated analytical signal AS h value is independent of the magnitude of the magnetic dip value and can be used to indicate the location of the magnetic anomaly source region and the magnetic field strength.
  • the AS h value at each edge point on each scale grid is used to indicate the structural magnetic strength.
  • Figure 11 is a structural grid diagram showing the magnetic strength using the AS h value at the edge point of the 300 m upper extension grid.
  • the gray-white to black gradient color indicates the gradual change of the structural band indicated by the analytical signal AS h value. Enhanced.
  • the magnetic scale structure maps of the plurality of scales are superposed to generate a comprehensive magnetic strength and weak structure grid diagram, as shown in FIG. 12, the AS h values of the scales or depth edges are gray-white to black gradient colors.
  • Overlay display highlights the magnetic change information of the texture grids in different depth ranges.
  • the magnitude of the three-dimensional analytical signal AS h indicates the magnetic strength of the anomalous source.
  • the strong magnetic structure is usually closely related to mineralization.
  • the structural grid map with the information of cutting depth and burial depth can clearly reflect the main structure with deep and long extension in the region, and the secondary structure with relatively shallow depth and relatively short extension.
  • the cover thickness of each configuration and the mutual delivery relationship, and thus the structural grid map obtained by the method of the present invention can help those skilled in the art to recognize the detected regional structural pattern.
  • the intersections along the magnetically strong, deeper structural belts, the different directions, and the transitional bends of the structural belts are important locations for the discovery of potential metal deposits.
  • the method according to the present invention can quickly and accurately extract and cut with magnetic measurement data.
  • Structural grid maps of depth, burial depth, magnetic strength, primary and secondary relationships, and delivery relationships help prospectors accurately and quickly identify potential metal deposits.

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Geology (AREA)
  • Environmental & Geological Engineering (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Electromagnetism (AREA)
  • Software Systems (AREA)
  • Geometry (AREA)
  • Computer Graphics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种位场构造格架自动提取方法。该方法包括以下步骤：对来自待研究区域的重力位场数据和/或磁法位场数据进行预处理；对经预处理的重力位场数据和/或磁法位场数据进行多个尺度下多方向的边缘检测，分别得到各尺度的边缘；采用形态学骨骼算法将计算得到的各尺度的边缘细化为单像素宽度，并且每个点具有深度和强度属性，得到综合构造格架图和综合构造强度格架图。可以实现对控制矿床形成的地质构造的识别和定性解释，根据研究区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。

Description

位场构造格架自动提取方法 技术领域
本发明涉及一种重力位场和磁位场构造格架自动提取方法，是一项基于重力和磁异常数据进行地质构造探测的技术。更具体地，本发明涉及小波分析、图像处理、地球物理、地质学、矿产勘查等领域，根据本发明的方法能直接应用于矿产勘探和相关地质调查领域等。
背景技术
随着航空物探技术的发展，重力和磁法(下文中也简称为重磁)测量方法具有经济、快速、可以覆盖很多难以到达的景观地区等优点，在金属矿床勘查评价中发挥越来越重要的作用，特别是高精度航磁测量技术的发展，使基于磁异常数据进行控制矿床形成的地质构造探测方法，在金属矿床勘探的选区到靶区定位的各个阶段都具有重要意义。目前在矿产勘查领域，重力和磁法勘探方法通常被应用于具有强磁性或高密度的矿化体的直接探测，以及具有强异常的地质构造解释和反演。大部分有色金属和稀有金属矿床不能直接形成明显重力和磁异常，但是控制矿床形成的地质构造可以产生异常差异。所以，基于重力位场数据和磁位场数据进行控制矿床形成的地质构造探测，在金属矿床靶区定位中具有重要意义。
目前重磁位场构造信息自动识别和提取方法主要包括解析信号法、欧拉反褶积法、位场多尺度边缘检测方法等，其共同存在的问题是对方向性的信息不敏感，不能得到完整、准确的异常边界位置。在本发明申请人的发明名称为“位场多方向多尺度边缘检测方法”，申请号为200810006676.6的专利申请中，公开了一种位场多方向多尺度边缘检测方法，通过方向小波变换增强了方向性信息，获得了不同方向的异常源边界信息，实现了构造格架的自动提取，克服了解析信号法、欧拉反褶积法、位场多尺度边缘检测方法等对方向性信息不敏感的问题，是一项基于位场数据进行浅部地壳三维结构快速反演的技术方案。但该方法仍存在以下问题：(1)计算得到的边缘不是单像素点宽度，边缘及不同方向边缘交点所对应的实际地理范围较大，导致分析结果精确度不高；(2)需要基于不同尺度边缘进行人工矢量化，取边缘中心连线作为异常边界，从而生成浅部地壳三维结构图，导致工作效率低；(3)该浅部地壳三维结构图只能显示不同深度的构造信息，不能反映构造带及两 侧岩性变化，不能指示构造的变形和活动强度；(4)对尺度没有明确定义，也就是没有对尺度这一特征赋予具体的地球物理属性。因此，需要提供一种能够基于已采集的位场数据高精度地提取位场构造格架的方法。
此外，受磁化方向影响，磁异常相对于重力异常更复杂。现代磁力仪通常记录的总磁场强度(Total Magnetic Intensity，缩写TMI)，相当于平行于地球主磁场方向的分量。通过对TMI进行地磁正常场校正可得到TMI异常。受斜磁化的影响，TMI异常存在侧向偏移、形态变形、正负值改变等问题。通常需要对TMI异常数据进行化极(Reduction to pole)处理来消除这样因素的影响。化极将观测的TMI异常转换成垂直磁化情况下的垂直磁异常，即将观测的TMI异常转换成可以在北磁极测量的异常，从而将磁异常迁移至源区正上方，便于磁异常的地质解释。
但是，由于受较小的磁倾角和噪声的影响，对低纬度(通常指磁倾角在±20°之间)地区的TMI异常数据进行化极难以得到可靠的垂直磁化磁异常数据。现有的基于化极后TMI异常数据进行磁构造信息自动识别和提取的方法，例如解析信号法、欧拉反褶积法、相位对称法(Phase symmetry)、位场多尺度边缘检测方法、位场多方向多尺度边缘检测方法等，均不适用于低纬度地区的磁构造信息自动识别。
因此，需要提供一种能够基于低纬度地区磁异常数据高精度地提取地质构造格架的方法。
发明内容
本发明目的在于针对现有技术的上述不足，提供一种重磁位场构造格架自动提取方法，能够快速得到控制矿床形成的地质构造信息，从而实现金属矿床靶区定位。本发明进一步提供一种低纬度地区磁构造格架自动提取方法，不仅可以得到线性构造(Lineaments)，还可以得到环形构造。本发明将由此得到的构造定义为构造格架(Structural framework)。
为实现这样的目的，本发明对测量得到的重力位场数据和磁位场数据采用多方向多尺度边缘检测方法提取重力构造格架和磁构造格架，采用形态学骨骼算法将得到的每一尺度构造格架细化处理成单像素宽度，将不同尺度的构造格架用不同颜色渐变叠加显示突出不同深度构造信息，生成综合构造格架图；将基于不同尺度构造格架各边缘点对应的梯度的模提取的构造强度值进行不同渐变颜色的叠加显示，突出密度变化和磁性变化强度信息，生成综合构造强度格架图。从而可以获得该研究区范围不同深度的重磁异常信息、 表征不同深度地质构造的构造格架分布信息、不同深度构造格架的密度变化和磁性变化强度信息，实现对控制矿床形成的地质构造的识别和定性解释，根据研究分析区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。
本发明的位场构造格架自动提取方法，包括如下的具体步骤：
1)对从研究区域重磁测量获得的重力数据或磁法数据进行预处理
该预处理包括对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常；对重力数据进行预处理得到布格重力异常数据。
2)对经预处理的重力位场数据或磁法位场数据进行多个尺度下位场多方向边缘检测，包括将位场数据上延多个预定高度后得到多个尺度的重力位场数据或磁法位场数据，并分别对所得到的每一尺度下的位场数据进行多方向边缘检测，得到各尺度下位场边缘。
针对上延多个预定高度后的每一尺度，选择不同的方向α进行边缘检测计算，可突出不同方向的边缘信息。为了能够达到完整覆盖，方向α取值为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n为大于或等于2的整数。针对每一方向，经预处理的重力数据或磁法数据的二维方向小波变换的模极大值点沿梯度的垂直方向连接得到的曲线构成边缘。获取同一尺度以不同方向计算得到的边缘，对所获取得到的边缘求并集作为该尺度的边缘，从而可实现各尺度下位场多方向的边缘检测。
3)采用形态学骨骼算法将计算得到的各尺度的边缘分别细化处理为单像素宽度，得到各尺度下的构造格架。
4)将计算得到的各尺度构造格架叠置生成综合构造格架图。
对得到的各尺度的构造格架进行叠置生成反映不同深度信息的综合构造格架图，图上不同尺度边缘的横向偏移反应了构造格架的产状信息。
将位场上延不同高度后提取的边缘对应于不同深度的构造，深度为向上延拓后高度的一半(见作者为Jacobsen,B.H.，题目为A case for upward continuation as a standard separation filter for potential-field maps，期刊名Geophysics，刊号v.52no.8，时间1987)，可以得到表征不同深度的构造格架图。
5)将计算得到的边缘上的每一边缘点处的梯度的模作为构造格架中该边缘点处强度值。赋有强度值的构造格架上边缘点强度值的高低反映了构造带及两侧岩性变化大小，反映了构造的变形和活动强度。建立反映不同深度信息的不同尺度的强度格架图，以及不同尺度叠加成综合构造强度格架图，揭 示了区域地质构造格局。
进一步，TMI异常数据的倾斜导数(Tilt derivative，缩写TDR)的总水平导数(Total horizontal derivatives)和解析信号(Analytic signal)均不受磁倾角影响，计算结果与磁倾角值大小无关，可直接对TMI异常数据进行TDR的总水平导数和解析信号的计算，而无需进行化极处理。但是，解析信号法的计算结果将增大磁异常范围，损失地质构造产状信息和构造分区信息，并且对地质构造的识别不敏感。然而，目前倾斜导数的总水平导数的结果均以网格图像或等值线图的形式表达，无法表征构造深度、主次关系、切割关系、磁性强弱等对于地质解释和找矿有重要意义的信息。传统的基于水平梯度的边缘检测方法因为没有考虑数据的方向性信息，而不能得到完整、准确的磁异常源体(Source bodies of magnetic anomalies)边界位置。本发明通过对经预处理得到的TMI异常数据进行多个尺度下基于其倾斜导数的总水平导数的多方向边缘检测，可有效识别和建立低纬度地区的磁构造格架。因此，根据本发明的构造格架自动提取方法特别适用于低纬度地区磁法测量数据的构造格架自动提取。
根据本发明的一个方面，提供一种位场构造格架自动提取方法，包括以下步骤：
对来自待研究区域的位场数据进行预处理；
将经预处理的位场数据向上延拓多个预定高度得到多个相应尺度的位场数据；
分别针对每一尺度的位场数据进行多方向边缘检测，得到多个相应尺度的位场边缘；
采用形态学骨骼算法将计算得到的各尺度的位场边缘分别细化为单像素宽度，得到多个相应尺度的构造格架图。
优选地，该方法进一步包括将计算得到的所述多个相应尺度的构造格架图叠置生成综合构造格架图。
优选地，该方法进一步包括将各尺度构造格架图上的每一边缘点处的梯度的模作为该尺度构造格架图中该边缘点处的强度值，得到多个相应尺度的构造强度格架图。
优选地，该方法进一步包括将所述多个相应尺度的构造强度格架图叠置生成综合构造强度格架图。
优选地，所述位场数据为重力位场数据或磁法位场数据，所述预处理进一步包括对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪 重力异常；或对重力数据进行预处理得到布格重力异常。
优选地，针对每一尺度的位场数据进行多方向边缘检测，包括以下步骤：
设尺度参数s＝z/z0，且z＞z0，z0代表测量高度，z代表向上延拓后的高度，定义高度为零的位置(x,y)处的重力异常或磁异常为f0(x,y)，
尺度s的平滑函数定义为：
Figure PCTCN2016075626-appb-000001
其中
Figure PCTCN2016075626-appb-000002
k(x,y,z)为格林函数，
在方向α的小波函数定义为：
Figure PCTCN2016075626-appb-000003
其中，D表示一阶导数；
在尺度s和位置(x,y)的情况下，重力异常或磁异常f0(x,y)在方向α的小波变换定义为：
Figure PCTCN2016075626-appb-000004
其中，*表示卷积运算，
由位场向上延拓公式已知
Figure PCTCN2016075626-appb-000005
fz(x,y)为f0(x,y)从高度零上延到高度z＝sz0处的重力异常或磁异常，通过将在测量高度z0测量得到的重力异常或磁异常
Figure PCTCN2016075626-appb-000006
向上延拓高度z-z0获得，
因此，
Figure PCTCN2016075626-appb-000007
进一步，在尺度s和位置(x,y)，重力异常或磁异常f0(x,y)在方向
Figure PCTCN2016075626-appb-000008
的小波变换定义为：
Figure PCTCN2016075626-appb-000009
则，f0(x,y)的二维方向小波变换用梯度表示为：
Figure PCTCN2016075626-appb-000010
Figure PCTCN2016075626-appb-000011
其中
Figure PCTCN2016075626-appb-000012
为二维梯度。
对于位置(x,y)、尺度s和方向α，f0(x,y)的二维方向小波变换W[f0](x,y,s,α)与fz(x,y)的梯度
Figure PCTCN2016075626-appb-000013
成正比，f0(x,y)的二维方向小波变换W[f0](x,y,s,α)可以用fz(x,y)的梯度
Figure PCTCN2016075626-appb-000014
来表征。
对于高度z,定义梯度
Figure PCTCN2016075626-appb-000015
的模为：
Figure PCTCN2016075626-appb-000016
该梯度相应的沿水平方向的辐角为：
Figure PCTCN2016075626-appb-000017
边缘点为模M[fz](x,y,α)沿辐角方向Afz(x,y,α)有局部极大值的点，
针对每一方向α，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘，
针对同一高度，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到该高度的边缘。
优选地，针对同一高度以多个不同的方向α计算边缘的步骤进一步包括，各方向α取值为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n为大于或等于2的整数，以完整覆盖二维平面。
优选地，位场向上延拓多个预定高度后提取的边缘对应于不同深度的构造，对得到的各尺度的构造格架图进行叠加得到反映不同深度信息的综合构造格架图。
优选地，以渐变颜色表征高度的高低，来形成所述综合构造格架图。
优选地，以渐变颜色表征强度值的大小，来形成所述综合构造强度格架图。
根据本发明的另一方面，提供一种构造格架自动提取方法，该方法包括以下步骤：
对来自待研究区域的磁法或重力测量数据进行预处理，得到总磁场强度TMI异常数据或布格重力异常数据；
将所得到的TMI异常数据或布格重力异常数据网格化，并将网格化的 TMI异常数据或布格重力异常数据向上延拓多个预定高度，得到多个尺度的网格化TMI异常数据或布格重力异常数据Th，h为向上延拓后的高度；
分别利用每一尺度的网格化TMI异常数据或布格重力异常数据Th计算各尺度的TMI异常数据或布格重力异常数据的倾斜导数TDRh；
分别针对每一尺度的网格化TMI异常数据或布格重力异常数据的倾斜导数，进行基于水平梯度的多方向边缘检测，得到各尺度的磁或重力异常源体边缘；
采用形态学骨骼算法将计算得到的各尺度的磁或重力异常源体边缘分别细化为单像素宽度，得到多个尺度的构造格架图。
优选地，该方法进一步包括将计算得到的所述多个尺度的构造格架图叠置生成综合构造格架图。
优选地，该方法进一步包括将TMI异常数据或布格重力异常数据向上延拓多个预定高度后提取的边缘对应于不同深度的构造，对得到的各深度的构造格架图进行叠加得到反映不同切割深度信息的综合构造格架图。
优选地，所述分别针对每一尺度的网格化TMI异常数据或布格重力异常数据的倾斜导数进行基于水平梯度的多方向边缘检测，包括以下步骤：
倾斜导数TDRh在方向α和
Figure PCTCN2016075626-appb-000018
的方向导数分别定义为：
Figure PCTCN2016075626-appb-000019
Figure PCTCN2016075626-appb-000020
其中，D表示一阶导数；
对于高度h和方向α，倾斜导数TDRh的水平梯度表示为：
Figure PCTCN2016075626-appb-000021
其中
Figure PCTCN2016075626-appb-000022
为水平梯度；
定义水平梯度
Figure PCTCN2016075626-appb-000023
的模为：
Figure PCTCN2016075626-appb-000024
该水平梯度的辐角为：
Figure PCTCN2016075626-appb-000025
则，对于高度h方向α的磁或重力异常源体边缘点为模
Figure PCTCN2016075626-appb-000026
沿辐角方向
Figure PCTCN2016075626-appb-000027
有局部极大值的点；
针对每一方向α，将倾斜导数TDRh的水平梯度的模的局部极大值点沿梯度的垂直方向连接，得到的曲线构成边缘；
针对同一高度h，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到相应尺度的磁或重力异常源体边缘，
其中，所述多个不同方向α取值为分别kπ/(2n-1)，k＝0,1,2,…,(2n-1)，n为大于或等于2的整数，以完整覆盖二维平面。
优选地，该方法进一步包括将各尺度构造格架图上的每一边缘点处的水平梯度的模表征该尺度构造格架图中该边缘点处的构造埋藏深度，得到多个尺度的表征构造埋藏深度的构造格架图。
优选地，该方法进一步包括将所述多个尺度的表征构造埋藏深度的构造格架图叠置生成综合埋藏深度构造格架图。
优选地，该方法进一步包括分别基于每一尺度的网格化TMI异常数据或布格重力异常数据Th计算三维解析信号ASh，得到各边缘点的ASh值，由此得到多个尺度的表征边缘点处磁性或密度强弱的构造格架图。
优选地，该方法进一步包括将所述多个尺度的表征边缘点磁性或密度强弱的构造格架图叠置生成综合磁性强弱构造格架图或综合密度强弱构造格架图。
优选地，该方法进一步包括，在进行边缘检测前对计算得到的倾斜导数进行去除噪声的处理。
优选地，该方法适用于低纬度地区的磁构造格架自动提取；优选地，该方法适用于磁倾角为±30°之间的区域的磁法测量数据；进一步优选地，该方法适用于磁倾角为±20°之间的区域的磁法测量数据。
根据本发明的方法，可以获得该研究区范围不同深度的重磁异常信息、表征不同深度地质构造的构造格架分布信息、不同深度构造格架的磁性和密度变化强度信息，实现对控制矿床形成的地质构造的识别和定性解释，根据研究区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。
本发明的利用低纬度地区磁法测量数据自动提取构造格架的方法，解决了现有技术中低纬度地区难以利用磁法测量数据准确获得构造信息的问题。根据本发明所得到的构造格架图相比于现有的网格图像或等值线图直观表征了构造切割深度、埋藏深度、主次关系、交割关系、磁性强弱等对于地质解 释和找矿有重要意义的信息。
本发明提出的低纬度地区磁构造格架自动提取方法同样适用于中高纬度地区的磁构造格架的自动提取，另外同样适用于重力位场构造格架的自动提取。
根据本发明的方法，延伸了利用磁法手段分析获取构造格架的区域，提高了自动提取构造格架的准确度，可以实现对控制矿床形成的地质构造的识别和定性解释，根据研究区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。
附图说明
从参照附图对实例性实施例的以下描述，本发明的进一步特征将变得明显。此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为根据本发明第一实施例的位场构造格架自动提取方法流程图；
图2根据本发明第一实例的单像素宽度构造格架图；
图3为根据本发明第一实例的综合构造格架图；
图4为根据本发明第一实例的构造强度格架图；
图5根据本发明第一实例的综合构造强度格架图；
图6为根据本发明第二实施例的低纬度地区磁构造格架自动提取方法流程图；
图7为根据本发明第二实例的单像素宽度构造格架图；
图8为根据本发明第二实例的综合构造格架图；
图9为根据本发明第二实例的反映埋藏深度的构造格架图；
图10根据本发明第二实例的综合埋藏深度构造格架图；
图11为根据本发明第二实例的反映磁性强弱构造格架图；
图12根据本发明第二实例的综合磁性强弱构造格架图。
具体实施方式
为了更好的理解本发明的技术方案，以下结合附图和具体实施例对本发明的实施方式作进一步的描述，但不作为对本发明的限定。
第一实施例
图1为本发明第一实施例的位场构造格架自动提取方法的流程图。如图1 所示，该方法包括如下步骤：
步骤101，对测量得到的重力位场数据或磁位场数据预处理。
对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常。
对重力数据进行预处理得到布格重力异常。
步骤102，采用多个尺度下位场多方向边缘检测方法对经预处理的重力位场数据或磁位场数据进行处理得到多个尺度的边缘。
位场多方向多尺度边缘检测方法包括将经预处理的重力位场数据或磁法位场数据向上延拓多个预定高度得到多个相应尺度的位场数据和分别针对每一尺度的位场进行多方向边缘检测得到各尺度位场边缘的步骤。
分别针对每一尺度的位场进行多方向边缘检测的计算方法，包括以下步骤：
定义地球表面点(x,y)处高程为零时的重力异常或磁异常为f0(x,y)。
设尺度参数s＝z/z0，且z＞z0，z0代表测量高度，z代表向上延后的高度。
尺度s的平滑函数定义为：
Figure PCTCN2016075626-appb-000028
其中
Figure PCTCN2016075626-appb-000029
k(x,y,z)为格林函数。
在方向α的小波函数可以定义为：
Figure PCTCN2016075626-appb-000030
其中，D表示一阶导数；
在尺度s和位置(x,y)的情况下，重力异常或磁异常f0(x,y)在方向α的小波变换定义为：
Figure PCTCN2016075626-appb-000031
其中，*表示卷积运算，
由位场向上延拓公式已知：
Figure PCTCN2016075626-appb-000032
fz(x,y)为f0(x,y)从高度零上延到高度z＝sz0处的重力异常或磁异常，通过 将在测量高度z0测量得到的重力异常或磁异常
Figure PCTCN2016075626-appb-000033
向上延拓高度z-z0获得，
因此，
Figure PCTCN2016075626-appb-000034
同样，在尺度s和位置(x,y)，重力异常或磁异常f0(x,y)在方向
Figure PCTCN2016075626-appb-000035
的小波变换定义为：
Figure PCTCN2016075626-appb-000036
f0(x,y)的二维方向小波变换可以用梯度来表示：
Figure PCTCN2016075626-appb-000037
Figure PCTCN2016075626-appb-000038
其中
Figure PCTCN2016075626-appb-000039
为二维梯度，
上式建立了任意高度z(z>z0)的重力异常或磁异常fz(x,y)的水平梯度与高度为零的重力异常或磁异常f0(x,y)的二维方向小波变换的联系。
对于位置(x,y)、尺度s和方向α，f0(x,y)的二维方向小波变换W[f0](x,y,s,α)与fz(x,y)的梯度
Figure PCTCN2016075626-appb-000040
成正比，f0(x,y)的二维方向小波变换W[f0](x,y,s,α)可以用fz(x,y)的梯度
Figure PCTCN2016075626-appb-000041
来表征。
对于向上延拓后的高度z,定义梯度
Figure PCTCN2016075626-appb-000042
的模为：
Figure PCTCN2016075626-appb-000043
该梯度相应的沿水平方向的辐角为：
Figure PCTCN2016075626-appb-000044
边缘点就是模M[fz](x,y,α)沿辐角方向Afz(x,y,α)有局部极大值的点。
针对每一向上延后的高度z，选择不同的方向α进行计算，可突出不同方向的边缘信息。为了能够达到二维平面完整覆盖，各方向α取值为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n为大于或等于2的整数。针对每一方向，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘。在同一上延后的高度下，实现完整覆盖二维平面的各方向α计算边缘，对计算得到的各边缘求并集得到该尺度的边缘。
对每一尺度的位场进行上述计算，可实现各尺度下位场多方向的边缘检 测。
步骤103，采用形态学骨骼算法对计算得到的各尺度的边缘进行细化处理成单像素宽度，得到各个尺度的格架图。
对上述计算得到的边缘图像采用Lam,L.,Seong-Whan Lee,and Ching Y.Suen,Thinning Methodologies-A Comprehensive Survey,IEEE Transactions on Pattern Analysis and Machine Intelligence,v.14,no.9,September 1992的细化格架算法将边缘细化成单像素点宽度。
通过对计算得到的边缘进行细化，边缘及不同方向边缘交点所对应的实际地理范围较现有技术明显缩小，使得到的位场构造格架更接近于实际地质填图识别构造的特点，一方面图面信息清晰，增强了可读性，另一方面也便于进行地质解释。
步骤104，将计算得到的各尺度构造格架叠置生成综合构造格架图。
位场上延不同高度后提取的各尺度的边缘叠置，得到待研究区域的位场综合构造格架图。
可将位场上延后的不同高度对应于不同源区深度，见作者为Jacobsen,B.H.，题目为A case for upward continuation as a standard separation filter for potential-field maps，Geophysics，v.52no.8，1987的文章，得到对应地球不同深度的构造。用渐变的颜色表示相应不同尺度或不同深度的格架，并将用不同的渐变颜色表征的各尺度构造格架叠加以突出显示不同深度的构造信息。通过如此叠加得到反映不同深度信息的综合构造格架图。
步骤105，用各尺度格架上各边缘点的梯度的模反映边缘强度。
针对不同的方向α，各边缘点处梯度的模M[fz](x,y,α)大小不变，为一固定值，取边缘点处梯度的模代表该边缘点处构造强度值，建立不同尺度或不同深度构造强度格架图。
步骤106，将所述多个尺度的构造强度格架图叠置生成综合构造强度格架图。
对各尺度或深度边缘的强度值进行不同渐变颜色的叠加显示，可突出显示磁性和密度变化强度信息。
第一实例
下面以采用云南西部航磁数据进行位场构造格架自动提取为例，对本发明的技术方案进行解释。
首先，对航磁数据进行分块化极处理，将处理后的数据拼接成位场网格文件，网格大小为500米。本实例中的航磁数据为历史多次分别采集的航磁 数据，测量高度在800-1200米范围内。
随后，将拼接的位场网格文件分别进行向上延拓处理得到高度上延后的位场数据，上延高度分别为1000、1500、2000、2500、3000、4000、5000、10000、15000、20000、25000、30000米，针对每一上延后高度进行32个方向的边缘检测，各方向α取值为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n＝5，得到各尺度的边缘。
随后，采用骨骼算法对计算得到的各尺度的边缘进行细化处理，得到各个尺度的构造格架图。
图2为上延5000米上延高度边缘经骨骼算法细化提取的相应一个尺度的单像素宽度构造格架。可以看出，根据本发明的自动提取得到的单像素点宽度的位场构造格架更接近于实际地质填图识别构造特点，便于进行地质解释。另外图面信息更清晰，便于不同尺度格架的叠加分析。
图3为将位场分别向上延拓1000、1500、2000、2500、3000、4000、5000、10000、15000、20000、25000、30000米上延高度后经32个方向边缘检测得到的各尺度的格架进行颜色渐变叠加形成的综合构造格架图。将上延后的各个高度对应于相应的源区深度，可将该综合构造格架图用于表征研究区域不同深度的构造格架信息。用从灰白色到黑色渐变色代表了从低到高的向上延拓高度或从浅入深的深度，综合构造格架图反映了不同深度的构造信息。
图4为上延5000米上延高度格架的边缘点处的梯度的模所反应的构造强度，从灰白到黑色反应了梯度的模所代表的构造强度值逐渐增大。该图反映了相应一个尺度下构造强度变化。
图5为将位场分别上延1000、1500、2000、2500、3000、4000、5000、10000、15000、20000、25000、30000米上延高度后得到的边缘强度进行叠加得到的综合构造强度格架图。从图中可以看到，在不同尺度上，反映出主要构造带均表现为构造强度较大，对应的磁异常变换明显，指示了主要构造带均为磁异常的突变带，同时对应的深度非常大，反映了构造带控制了深部岩浆-成矿作用。
另外，图3和5中，具有深度和强度信息的构造格架图能够清楚的反映出区域中深度大、延展长的主干构造，和深度较浅、延展相对短的次级构造，以及互相的交割关系，因此根据本发明的方法得到的构造格架图可以帮助本领域技术人员认识所探测的区域构造格局。
将图3和图5与该区域的地表填图地质构造对比，空间位置和范围非常吻合，说明了根据本发明的自动提取方法的准确性和有效性。此外，与区域 的地表填图地质构造对比，本发明的图3和图5增加了构造三维延伸、强度、结构等信息，可以帮助识别隐伏构造带。
沿大深度、高构造强度的构造带附近，不同方向构造的交汇部位以及构造带转折弯曲部位是发现潜在金属矿床的重要位置，根据本发明的方法利用航磁数据和重力数据可以快速准确提取得到具有深度和强度信息的构造格架图，帮助勘探人员准确快速地发现潜在的金属矿床。
第二实施例
下面以磁法测量数据为例，具体说明根据本发明第二实施例的低纬度地区磁构造格架自动提取方法的流程图。本领域技术人员可以理解，本发明的方法不仅限于磁法测量数据，同样适用于重力测量数据的构造格架自动提取。
图6为本发明的低纬度地区磁构造格架自动提取方法的流程图。如图6所示，该方法包括如下步骤：
步骤601，对磁法测量得到的观测值进行预处理，进行地磁正常场(IGRF)改正，得到总磁场强度(TMI)异常数据。
对TMI异常数据进行网格化，通常网格大小取值为测线间距的1/8到1/4或最小点距。
步骤602，对网格化TMI异常数据T向上延拓(upward continuation)多个预定高度得到多个相应尺度的网格化TMI异常数据Th，h代表向上延拓后的高度。
步骤603，分别利用每一尺度的网格化TMI异常数据Th计算相应尺度TMI异常数据的倾斜导数TDRh：
Figure PCTCN2016075626-appb-000045
其中，VDRh和THDRh分别为网格化TMI异常数据Th的垂向一阶导数和总水平导数:
Figure PCTCN2016075626-appb-000046
Figure PCTCN2016075626-appb-000047
步骤604，分别针对每一尺度的倾斜导数TDRh，计算基于水平梯度的多方向边缘检测。
倾斜导数TDRh在方向α和
Figure PCTCN2016075626-appb-000048
的方向导数分别定义为：
Figure PCTCN2016075626-appb-000049
Figure PCTCN2016075626-appb-000050
其中，D表示一阶导数；
对于高度h和方向α，倾斜导数TDRh的水平梯度表示为：
Figure PCTCN2016075626-appb-000051
其中
Figure PCTCN2016075626-appb-000052
为水平梯度；
定义水平梯度
Figure PCTCN2016075626-appb-000053
的模为：
Figure PCTCN2016075626-appb-000054
该梯度相应的辐角为：
Figure PCTCN2016075626-appb-000055
TMI异常数据的倾斜导数的总水平导数(缩写TDR_THDR)为：
Figure PCTCN2016075626-appb-000056
模
Figure PCTCN2016075626-appb-000057
的幅值大小与磁倾角的大小无关。
由此，对于高度h和方向α的磁源体边缘点为模
Figure PCTCN2016075626-appb-000058
沿辐角方向
Figure PCTCN2016075626-appb-000059
有局部极大值的点，
针对每一方向α，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘；
针对同一高度h，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到相应尺度的磁源体边缘。
为了能够达到二维平面完整覆盖，各方向α取值为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n为大于或等于2的整数，以完整覆盖二维平面。
由于对网格化TMI异常数据的TDR计算可能在计算结果中引入噪声，而其后的多方向边缘检测计算对噪声很敏感。为避免噪声对多方向边缘检测计算的影响，优选地，在计算多方向边缘检测前对于存在较大的噪声的TDRh数据进行例如高斯滤波的降噪处理。
因为TMI异常数据的倾斜导数的总水平导数不受磁倾角影响，计算结果与磁倾角值大小无关，采用上述步骤得到的磁异常源体边缘不受磁倾角的影响即不受待分析区域的纬度位置的影响，可以准确表征低纬度地区的构造格 架。
步骤605，采用形态学骨骼算法将计算得到的各尺度的边缘分别细化处理为单像素宽度，得到多个尺度下的构造格架图。
对上述计算得到的边缘图像采用Lam,L.,Seong-Whan Lee,and Ching Y.Suen,Thinning Methodologies-A Comprehensive Survey,IEEE Transactions on Pattern Analysis and Machine Intelligence,v.14,no.9,September 1992的细化格架算法将边缘细化成单像素点宽度。
通过对计算得到的边缘进行细化，边缘及不同方向边缘交点所对应的实际地理范围较现有技术明显缩小，使得到的构造格架更接近于实际地质填图识别构造的特点，一方面图面信息清晰，增强了可读性，另一方面也便于进行地质解释。
步骤606，将计算得到的各尺度的构造格架图叠置生成综合构造格架图。
对得到的各尺度的构造格架图进行叠置生成反映不同深度信息的综合构造格架图，图上不同尺度边缘的横向偏移反应了构造格架的产状信息。
将网格化的TMI异常数据上延不同高度后提取的边缘对应于不同深度的构造，深度为向上延拓后高度的一半(见作者为Jacobsen,B.H.，题目为A case for upward continuation as a standard separation filter for potential-field maps，期刊名Geophysics，刊号v.52no.8，时间1987)，可以得到表征不同切割深度的构造格架图。对得到的各深度的构造格架图进行叠加得到反映不同切割深度信息的综合构造格架图。
步骤607，用各尺度构造格架上各边缘点处的TDRh的水平梯度的模表征构造埋藏深度，得到多个尺度的表征构造埋藏深度的构造格架图。
针对不同的方向α，各边缘点处水平梯度的模
Figure PCTCN2016075626-appb-000060
大小不变，为一固定值。
由于三角函数arctan值的属性，不管VDRh和THDRh的幅值为多大，TDRh的幅值都限制在-π/2和+π/2之间。因此，各边缘点处水平梯度的模
Figure PCTCN2016075626-appb-000061
的大小与TMI异常幅值大小关系不大，该值反应了源区的埋藏深度，值的大小与埋藏深度成反比，值越大源区埋藏深度越浅(见作者Bruno Verduzco等，题目为New insights into magnetic derivatives for structural mapping，期刊名The Leading Edge,刊号v.23no.2，时间2004)。
倾斜导数TDRh的水平梯度的模反映的构造埋藏深度(Cover depth)代表了构造带上面覆盖层的厚度，而向上延拓多个预定高度对应的不同深度反映 了构造带的切割深度，代表了构造带向下延伸的深度。
取边缘点处梯度的模
Figure PCTCN2016075626-appb-000062
代表该边缘点处该尺度下的边缘相对埋藏深度，建立不同尺度或不同深度反应构造相对埋藏深度的构造格架图。对于同一个构造带的梯度的模值的大小变化，反应了该构造带不同部位的埋藏深度。
步骤608，将所述多个尺度的表征埋藏深度的构造格架图叠置生成综合埋藏深度构造格架图。
对各尺度或深度边缘点的TDRh水平梯度的模值采用不同渐变颜色的叠加显示，突出显示不同切割深度范围的构造格架的埋藏深度变化信息。
步骤609，对每一高度的网格化TMI异常数据Th计算三维解析信号ASh:
Figure PCTCN2016075626-appb-000063
计算得到的解析信号ASh幅值与TMI异常的幅值具有较强的相关性，但与磁倾角值大小无关，可以用来指示磁异常源区位置和磁场强度。
步骤610，用各尺度格架上各边缘点处ASh值指示构造的磁性强弱。
在不考虑磁异常源区埋藏深度的情况下，三维解析信号ASh的幅值大小指示了磁异常源的磁性强弱，磁性强的构造通常与成矿关系密切。
步骤611，将所述多个尺度的磁性强弱构造格架图叠置生成综合磁性强弱构造格架图。
对各尺度或深度边缘点的ASh值采用不同渐变颜色的叠加显示，突出显示不同深度范围的构造格架的磁性强弱变化信息。
第二实例
下面以采用位于低纬度地区的南美洲圭亚那地盾高精度航磁数据进行磁构造格架自动提取为例，对本发明的技术方案进行解释。
本实例中采用的航磁数据的测量比例尺为1:25000，测量高度在70-120米范围内。
首先，对研究区航磁测量得到的观测值进行预处理，进行地磁正常场(IGRF)改正，得到总磁场强度(TMI)异常数据，对TMI异常数据进行网格化，网格大小取值为10米。
随后，对网格化的TMI异常数据分别进行向上延拓处理得到多个尺度的网格化TMI异常数据Th，上延高度分别为100、200、300、400、500米。
对每一尺度的网格化TMI异常数据Th进行TDR计算；
分别针对每一个上延高度，对计算得到的TDR，进行基于水平梯度的多 方向边缘检测。针对每个TDR计算结果进行64个方向的边缘检测，各方向α取值分别为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n＝6，得到各尺度的边缘。
随后，采用骨骼算法对计算得到的各尺度的边缘进行细化处理，得到各个尺度的磁构造格架图。
图7为将网格化TMI异常数据上延300米上延高度，边缘经骨骼算法细化提取的相应一个尺度的单像素宽度磁构造格架。可以看出，根据本发明的自动提取得到的单像素点宽度的构造格架更接近于实际地质填图识别构造特点，便于进行地质解释。另外图面信息更清晰，便于不同尺度构造格架的叠加分析。
图8为将网格化TMI异常数据分别向上延拓100、200、300、400、500米上延高度后经64个方向边缘检测得到的各尺度的构造格架进行颜色渐变叠加形成的综合构造格架图。将上延后的各个高度对应于相应的源区深度，可将该综合构造格架图用于表征研究区域不同深度的构造格架信息。用从灰白色到黑色渐变色代表了从低到高的向上延拓高度或从浅入深的切割深度，综合构造格架图反映了不同切割深度的构造信息。
随后，计算各尺度格架上各边缘点处的TDRh的梯度的模以表征构造埋藏深度。图9为上延300米上延高度格架的边缘点处的梯度的模所表征的构造埋藏深度，从灰白到黑色渐变颜色反应了水平梯度的模所代表的构造埋藏深度逐渐减小。该图反映了相应一个尺度下构造埋藏深度变化。
随后，将所述多个尺度的表征相对埋藏深度的构造格架图叠置生成综合埋藏深度构造格架图。图10为将网格化的TMI异常数据分别上延100、200、300、400、500米上延高度后计算得到的表征相对埋藏深度构造格架图进行叠加得到的综合埋藏深度构造格架图。从图中可以看到，在不同尺度上，反映出主要构造带埋藏深度变化较大，通常埋藏深的构造即上覆地层较厚，对应的TMI异常幅值较低，为隐伏构造，因此该方法有助于隐伏构造的识别。
随后，对每一尺度的网格化TMI异常数据Th计算三维解析信号ASh。计算得到的解析信号ASh值与磁倾角值大小无关，可以用来指示磁异常源区位置和磁场强度。
该实例中，各尺度格架上各边缘点处ASh值用于指示构造磁性强弱。图11为上延300米上延高度格架的边缘点处采用ASh值指示磁性强弱的构造格架图，从灰白色到黑色渐变颜色表征了解析信号ASh值所指示的构造带磁性逐渐增强。
随后，将所述多个尺度的磁性强弱构造格架图叠置生成综合磁性强弱构 造格架图，如图12所示，各尺度或深度边缘的ASh值采用灰白色到黑色渐变颜色的叠加显示，突出显示不同深度范围的构造格架的磁性变化信息。
在不考虑源区埋藏深度的情况下，三维解析信号ASh的幅值大小指示了异常源的磁性强弱，磁性强的构造通常与成矿关系密切。
另外，同时观察图8和图10，可以看到具有切割深度和埋藏深度信息的构造格架图能够清楚的反映出区域中深度大延展长的主干构造、深度较浅延展相对短的次级构造、各构造的覆盖厚度以及互相的交割关系，因此根据本发明的方法得到的构造格架图可以帮助本领域技术人员认识所探测的区域构造格局。
沿磁性较强、深度较大的构造带附近、不同方向构造的交汇部位以及构造带转折弯曲部位是发现潜在金属矿床的重要位置，根据本发明的方法利用磁测数据可以快速准确提取得到具有切割深度、埋藏深度、磁性强弱、主次关系和交割关系等信息的构造格架图，帮助勘探人员准确快速地发现潜在的金属矿床。
以上借助优选实施例对本发明进行了详细说明，但是本发明不限于此。本技术领域技术人员可以根据本发明的原理进行各种修改。因此，凡按照本发明原理所作的修改，都应当理解为落入本发明的保护范围。

Claims (20)

  1. 一种位场构造格架自动提取方法，包括以下步骤：
    对来自待研究区域的位场数据进行预处理；
    将经预处理的位场数据向上延拓多个预定高度得到多个相应尺度的位场数据；
    分别针对每一尺度的位场数据进行多方向边缘检测，得到多个相应尺度的位场边缘；
    采用形态学骨骼算法将计算得到的各尺度的位场边缘分别细化为单像素宽度，得到多个相应尺度的构造格架图。
  2. 如权利要求1所述的位场构造格架自动提取方法，其特征在于，该方法进一步包括将计算得到的所述多个相应尺度的构造格架图叠置生成综合构造格架图。
  3. 如权利要求1所述的位场构造格架自动提取方法，其特征在于，该方法进一步包括：
    将各尺度构造格架图上的每一边缘点处的梯度的模作为该尺度构造格架图中该边缘点处的强度值，得到多个相应尺度的构造强度格架图。
  4. 如权利要求3所述的位场构造格架自动提取方法，其特征在于，该方法进一步包括将所述多个相应尺度的构造强度格架图叠置生成综合构造强度格架图。
  5. 如权利要求1所述的位场构造格架自动提取方法，其特征在于，所述位场数据为重力位场数据或磁法位场数据，所述预处理进一步包括
    对重力数据进行预处理得到布格重力异常；或
    对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常。
  6. 如权利要求1所述的位场构造格架自动提取方法，其特征在于，
    针对每一尺度的位场数据进行多方向边缘检测，包括以下步骤：
    设尺度s＝z/z0，且z＞z0，z0代表测量高度，z代表向上延拓后的高度，定义高度为零的位置(x,y)处的重力异常或磁异常为f0(x,y)，
    尺度s的平滑函数定义为：
    Figure PCTCN2016075626-appb-100001
    其中
    Figure PCTCN2016075626-appb-100002
    k(x,y,z)为格林函数，
    在方向α的小波函数定义为：
    Figure PCTCN2016075626-appb-100003
    其中，D表示一阶导数；
    在尺度s和位置(x,y)的情况下，重力异常或磁异常f0(x,y)在方向α的小波变换定义为：
    Figure PCTCN2016075626-appb-100004
    其中，*表示卷积运算，
    由位场向上延拓公式已知
    Figure PCTCN2016075626-appb-100005
    fz(x,y)为f0(x,y)从高度零上延到高度z＝sz0处的重力异常或磁异常，通过将在测量高度z0测量得到的重力异常或磁异常
    Figure PCTCN2016075626-appb-100006
    向上延拓高度z-z0获得，
    因此，Wα[f0](x,y,s)＝sDαfz(x,y)
    ＝(z/z0)Dαfz(x,y)；
    进一步，在尺度s和位置(x,y)，重力异常或磁异常f0(x,y)在方向
    Figure PCTCN2016075626-appb-100007
    的小波变换定义为：
    Figure PCTCN2016075626-appb-100008
    则，f0(x,y)的二维方向小波变换用梯度表示为：
    Figure PCTCN2016075626-appb-100009
    W[f0](x,y,s,α)＝(z/z0)▽fz(x,y,α)
    其中▽为二维梯度，
    对于高度z,定义梯度▽fz(x,y,α)的模为：
    Figure PCTCN2016075626-appb-100010
    该梯度相应的沿水平方向的辐角为：
    Figure PCTCN2016075626-appb-100011
    边缘点为模M[fz](x,y,α)沿辐角方向Afz(x,y,α)有局部极大值的点，
    针对每一方向α，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘，
    针对同一高度，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到相应尺度的位场边缘。
  7. 如权利要求6所述的位场构造格架自动提取方法，其特征在于，针对同一高度以多个不同的方向α计算边缘的步骤进一步包括，各方向α取值为kπ/(2n-1)，其中k＝0,1,2,…,(2n-1)，n为大于或等于2的整数，以完整覆盖二维平面。
  8. 如权利要求6所述的位场构造格架自动提取方法，其特征在于，位场向上延拓多个预定高度后提取的边缘对应于不同深度的构造，对得到的各尺度的构造格架图进行叠加得到反映不同深度信息的综合构造格架图。
  9. 如权利要求2或8所述的位场构造格架自动提取方法，其特征在于，以渐变颜色表征高度的高低，来形成所述综合构造格架图。
  10. 如权利要求4所述的位场构造格架自动提取方法，其特征在于，以渐变颜色表征强度值的大小，来形成所述综合构造强度格架图。
  11. 一种构造格架自动提取方法，该方法包括以下步骤：
    对来自待研究区域的磁法或重力测量数据进行预处理，得到总磁场强度(TMI)异常数据或布格重力异常数据；
    将所得到的TMI异常数据或布格重力异常数据网格化，并将网格化的TMI异常数据或布格重力异常数据向上延拓多个预定高度，得到多个尺度的网格化TMI异常数据或布格重力异常数据Th，h为向上延拓后的高度；
    分别利用每一尺度的网格化TMI异常数据或布格重力异常数据Th计算各尺度的TMI异常数据或布格重力异常数据的倾斜导数TDRh；
    分别针对每一尺度的网格化TMI异常数据或布格重力异常数据的倾斜导数，进行基于水平梯度的多方向边缘检测，得到各尺度的磁或重力异常源体边缘；
    采用形态学骨骼算法将计算得到的各尺度的磁或重力异常源体边缘分别细化为单像素宽度，得到多个尺度的构造格架图。
  12. 如权利要求11所述的构造格架自动提取方法，其特征在于，该方法进一步包括，将计算得到的所述多个尺度的构造格架图叠置生成综合构造格 架图。
  13. 如权利要求12所述的构造格架自动提取方法，其特征在于，该方法进一步包括，将网格化TMI异常数据或布格重力异常数据向上延拓多个预定高度后提取的边缘对应于不同深度的构造，对得到的各深度的构造格架图进行叠加得到反映不同切割深度信息的综合构造格架图。
  14. 如权利要求11所述的构造格架自动提取方法，其特征在于，所述分别针对每一尺度的网格化TMI异常数据或布格重力异常数据的倾斜导数进行基于水平梯度的多方向边缘检测，包括以下步骤：
    倾斜导数TDRh在方向α和
    Figure PCTCN2016075626-appb-100012
    的方向导数分别定义为：
    Figure PCTCN2016075626-appb-100013
    Figure PCTCN2016075626-appb-100014
    其中，D表示一阶导数；
    对于高度h和方向α，倾斜导数TDRh的水平梯度表示为：
    Figure PCTCN2016075626-appb-100015
    其中▽为水平梯度；
    定义水平梯度
    Figure PCTCN2016075626-appb-100016
    的模为：
    Figure PCTCN2016075626-appb-100017
    该水平梯度的辐角为：
    Figure PCTCN2016075626-appb-100018
    则，对于高度h方向α的磁或重力异常源体边缘点为模
    Figure PCTCN2016075626-appb-100019
    沿辐角方向
    Figure PCTCN2016075626-appb-100020
    有局部极大值的点；
    针对每一方向α，将倾斜导数TDRh的水平梯度的模的局部极大值点沿梯度的垂直方向连接，得到的曲线构成边缘；
    针对同一高度h，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到相应尺度的磁或重力异常源体边缘，
    其中，所述多个不同方向α取值为分别kπ/(2n-1)，k＝0,1,2,…,(2n-1)，n为大于或等于2的整数。
  15. 如权利要求14所述的构造格架自动提取方法，其特征在于，该方法进一步包括：
    分别将各尺度构造格架图上的每一边缘点处的水平梯度的模表征该尺度构造格架图中该边缘点处的构造埋藏深度，得到多个尺度的表征构造埋藏深度的构造格架图。
  16. 如权利要求15所述的构造格架自动提取方法，其特征在于，该方法进一步包括将所述多个尺度的表征构造埋藏深度的构造格架图叠置生成综合埋藏深度构造格架图。
  17. 如权利要求11所述的构造格架自动提取方法，其特征在于，该方法进一步包括，
    分别基于每一尺度的网格化TMI异常数据或布格重力异常数据Th计算三维解析信号ASh，得到各边缘点的ASh值，由此得到多个尺度的表征边缘点处磁性或密度强弱的构造格架图。
  18. 如权利要求17所述的构造格架自动提取方法，其特征在于，该方法进一步包括将所述多个尺度的表征边缘点磁性或密度强弱的构造格架图叠置生成综合磁性强弱构造格架图或综合密度强弱构造格架图。
  19. 如权利要求11所述的构造格架自动提取方法，其特征在于，该方法进一步包括，在进行边缘检测前对计算得到的倾斜导数进行去除噪声的处理。
  20. 如权利要求11所述的构造格架自动提取方法，其特征在于，该方法适用于低纬度地区的磁构造格架自动提取；优选地，该方法适用于磁倾角为±30°之间的区域的磁法测量数据；进一步优选地，该方法适用于磁倾角为±20°之间的区域的磁法测量数据。
PCT/CN2016/075626 2015-03-04 2016-03-04 位场构造格架自动提取方法 WO2016138874A1 (zh)

Priority Applications (4)

Application Number Priority Date Filing Date Title
RU2017133223A RU2664488C1 (ru) 2015-03-04 2016-03-04 Способ автоматического получения структурного строения из данных потенциального поля
CA2978500A CA2978500C (en) 2015-03-04 2016-03-04 Method for automatically extracting tectonic framework of potential field data
US15/555,153 US10884161B2 (en) 2015-03-04 2016-03-04 Method for automatically extracting structural framework from potential field data
AU2016228027A AU2016228027B2 (en) 2015-03-04 2016-03-04 Method for automatically extracting tectonic framework of potential field

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
CN201510096981.9 2015-03-04
CN201510096981.9A CN104658037B (zh) 2015-03-04 2015-03-04 一种位场构造格架自动提取方法
CN201510303316.2 2015-06-04
CN201510303316.2A CN104965232B (zh) 2015-06-04 2015-06-04 低纬度地区磁构造格架自动提取方法

Publications (1)

Publication Number Publication Date
WO2016138874A1 true WO2016138874A1 (zh) 2016-09-09

Family

ID=56848301

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2016/075626 WO2016138874A1 (zh) 2015-03-04 2016-03-04 位场构造格架自动提取方法

Country Status (5)

Country Link
US (1) US10884161B2 (zh)
AU (1) AU2016228027B2 (zh)
CA (1) CA2978500C (zh)
RU (1) RU2664488C1 (zh)
WO (1) WO2016138874A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117312898A (zh) * 2023-11-27 2023-12-29 山东省煤田地质规划勘察研究院 一种基于多重k均值聚类分析的找矿预测方法及系统

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108008459B (zh) * 2017-11-28 2019-11-01 北京中科地物能源技术有限公司 一种获得剩余重力异常的方法及装置
CN110135506B (zh) * 2019-05-20 2021-02-09 杭州电子科技大学 一种应用于web的七类皮肤肿瘤检测方法
CN110414060B (zh) * 2019-06-28 2023-01-03 中国地质大学（武汉） 一种基于四阶谱矩的位场边界识别方法
CN110703347B (zh) * 2019-10-24 2021-04-16 中国石油化工股份有限公司 基于构造背景的重力断裂影像识别方法
CN111325766B (zh) * 2020-02-20 2023-08-25 腾讯科技（深圳）有限公司 三维边缘检测方法、装置、存储介质和计算机设备
CN112328955B (zh) * 2020-10-16 2024-02-27 中国地质调查局沈阳地质调查中心 重磁数据的处理方法、存储介质及装置
CN113111500A (zh) * 2021-03-31 2021-07-13 中国地质大学（北京） 一种基于重磁物理场使用二维经验模态分解异常分析方法
CN113421194B (zh) * 2021-06-04 2022-07-15 贵州省地质矿产勘查开发局 一种根据布格重力异常图像提取隐伏断层的方法
CN113779013B (zh) * 2021-09-18 2024-01-30 核工业航测遥感中心 基于地质约束的小范围磁场数据补缺处理方法
US11815647B1 (en) * 2022-04-20 2023-11-14 Chinese Academy Of Geological Sciences Gravity inversion method and system based on meshfree method
CN115236755B (zh) * 2022-07-25 2023-10-03 中国自然资源航空物探遥感中心 基于张量特征值的航磁异常边界检测方法、装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7079135B2 (en) * 2002-03-07 2006-07-18 Samsung Electronics Co., Ltd. Method of wavelets-based multiresolution representation of three-dimensional image object
CN101256676A (zh) * 2008-01-31 2008-09-03 中国地质科学院矿产资源研究所 位场多方向多尺度边缘检测方法
CN102236108A (zh) * 2010-05-06 2011-11-09 中国石油天然气集团公司 一种磁性地表三维地形改正方法
CN102937725A (zh) * 2012-11-12 2013-02-20 中国科学院地质与地球物理研究所 一种基于过渡区与相叠合的位场异常边缘增强方法
CN104658037A (zh) * 2015-03-04 2015-05-27 中国地质科学院矿产资源研究所 一种位场构造格架自动提取方法
CN104965232A (zh) * 2015-06-04 2015-10-07 中国地质科学院矿产资源研究所 低纬度地区磁构造格架自动提取方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7477768B2 (en) * 1999-06-29 2009-01-13 The Research Foundation Of State University Of New York System and method for performing a three-dimensional virtual examination of objects, such as internal organs
RU2169384C1 (ru) * 1999-12-17 2001-06-20 Закрытое акционерное общество "Петербургская геофизическая компания" Способ поиска нефтегазовых месторождений
US7127104B2 (en) * 2004-07-07 2006-10-24 The Regents Of The University Of California Vectorized image segmentation via trixel agglomeration
US8041141B2 (en) * 2006-06-30 2011-10-18 The University Of Louisville Research Foundation, Inc. Method and software for shape representation with curve skeletons
CN101520518B (zh) 2008-02-25 2011-01-12 中国石油集团东方地球物理勘探有限责任公司 一种利用重磁电异常的组合特征识别火成岩岩性的方法
RU2401443C2 (ru) * 2008-03-17 2010-10-10 Общество с ограниченной ответственностью "Научно-производственное предприятие "Нейво" Способ обнаружения и отображения фигуры газонефтяной лог-трубки
US20110115787A1 (en) * 2008-04-11 2011-05-19 Terraspark Geosciences, Llc Visulation of geologic features using data representations thereof
NZ588826A (en) * 2008-05-22 2012-10-26 Exxonmobil Upstream Res Co Transforming seismic survey data to assess an area for hydrocarbon production potential
US8126247B2 (en) * 2009-05-19 2012-02-28 National Tsing Hua University Image preprocessing system for 3D image database construction
US9558588B2 (en) * 2012-06-26 2017-01-31 Schlumberger Technology Corporation Method for building a 3D model of a rock sample
CN103955007A (zh) 2014-05-20 2014-07-30 中国石油化工股份有限公司胜利油田分公司西部新区研究院 复杂山前构造带的综合建模方法及建立的地质结构模型
US10136869B2 (en) * 2016-03-25 2018-11-27 Perkinelmer Health Sciences, Inc. Systems and methods for characterizing a central axis of a bone from a 3D anatomical image

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7079135B2 (en) * 2002-03-07 2006-07-18 Samsung Electronics Co., Ltd. Method of wavelets-based multiresolution representation of three-dimensional image object
CN101256676A (zh) * 2008-01-31 2008-09-03 中国地质科学院矿产资源研究所 位场多方向多尺度边缘检测方法
CN102236108A (zh) * 2010-05-06 2011-11-09 中国石油天然气集团公司 一种磁性地表三维地形改正方法
CN102937725A (zh) * 2012-11-12 2013-02-20 中国科学院地质与地球物理研究所 一种基于过渡区与相叠合的位场异常边缘增强方法
CN104658037A (zh) * 2015-03-04 2015-05-27 中国地质科学院矿产资源研究所 一种位场构造格架自动提取方法
CN104965232A (zh) * 2015-06-04 2015-10-07 中国地质科学院矿产资源研究所 低纬度地区磁构造格架自动提取方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
CAO, DIANHUA: "Porphyry Copper Deposit Model and Exploration Technique in Zhongdian, Yunnan", CHINA DOCTORAL DISSERTATIONS FULL-TEXT DATABASE (BASIC SCIENCES, 15 February 2008 (2008-02-15), pages 53 - 54, ISSN: 1674-022X *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117312898A (zh) * 2023-11-27 2023-12-29 山东省煤田地质规划勘察研究院 一种基于多重k均值聚类分析的找矿预测方法及系统
CN117312898B (zh) * 2023-11-27 2024-03-15 山东省煤田地质规划勘察研究院 一种基于多重k均值聚类分析的找矿预测方法及系统

Also Published As

Publication number Publication date
CA2978500C (en) 2022-07-05
AU2016228027B2 (en) 2018-11-22
US20180052251A1 (en) 2018-02-22
CA2978500A1 (en) 2016-09-09
AU2016228027A1 (en) 2017-09-21
RU2664488C1 (ru) 2018-08-17
US10884161B2 (en) 2021-01-05

Similar Documents

Publication Publication Date Title
WO2016138874A1 (zh) 位场构造格架自动提取方法
Wang et al. A new edge recognition technology based on the normalized vertical derivative of the total horizontal derivative for potential field data
CN105510964B (zh) 复杂构造区低级序走滑断层的地震识别方法
EA018473B1 (ru) Способ представления геологических структур исследуемого района земли
Zahra et al. Application of high-pass filtering techniques on gravity and magnetic data of the eastern Qattara Depression area, Western Desert, Egypt
CN104965232B (zh) 低纬度地区磁构造格架自动提取方法
CN109425906B (zh) 一种磁异常探测矢量磁目标识别方法
Oliveira et al. EdgeDetectPFI: An algorithm for automatic edge detection in potential field anomaly images–application to dike-like magnetic structures
CN102879799A (zh) 多方位地震能量梯度差碳酸盐岩溶洞型储层识别方法
Zhang et al. NAV-Edge: Edge detection of potential-field sources using normalized anisotropy variance
CN109407161A (zh) 用于提取地球物理磁异常场边界的磁场刻痕分析方法
Khalifani et al. Generation of an efficient structural evidence layer for mineral exploration targeting
Beran et al. Detecting and classifying UXO
CN104658037B (zh) 一种位场构造格架自动提取方法
Pham A stable method for detecting the edges of potential field sources
Pregesbauer et al. An object oriented approach to automatic classification of archaeological features in magnetic prospection data
Prasad et al. Tectonic and structural elements of Southern Granulite Terrane, South India: Inferences from gravity and magnetic studies
Karbalaali et al. Identification of shallow geohazard channels based on discontinuity seismic attributes in the South Caspian Sea
Pei et al. Multi-scale edge detection method for potential field data based on two-dimensional variation mode decomposition and mathematical morphology
Argyriou A methodology for the rapid identification of neotectonic features using geographical information systems and remote sensing: A case study from Western Crete, Greece
Oyeniyi et al. An application of softsign function (SF) filter to low-latitude aeromagnetic data of Tafawa-Balewa Area, Northern Nigeria for geostructural mapping and tectonic analysis
Dhara et al. Automatic extraction and analysis of lineament features using ASTER and Sentinel 1 SAR data
Monsen et al. The automated interpretation of photorealistic outcrop models
Eshaghzadeh et al. Magnetic field interpretation using singular value decomposition method based on correlation coefficient of eigenimages
Ouadfeul et al. Multiscale analysis of geophysical signals using the 2D continuous wavelet transform

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16758489

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2978500

Country of ref document: CA

WWE Wipo information: entry into national phase

Ref document number: 15555153

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2016228027

Country of ref document: AU

Date of ref document: 20160304

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 2017133223

Country of ref document: RU

122 Ep: pct application non-entry in european phase

Ref document number: 16758489

Country of ref document: EP

Kind code of ref document: A1