A method for delineating fault structures in the sea-land transition zone based on aerial gravity data

Abstract

The present invention provides a method for delineating fault structures in a land-sea transition zone based on aerial gravity data, which relates to the fields of geophysics and geochemistry. The method fuses and splices aerial gravity data obtained from different sources to obtain regional aerial gravity field data covering a study area; uses the regional aerial gravity field data to obtain the depth of the Moho surface; uses notch analysis to perform boundary information enhancement processing on the regional aerial gravity field data to obtain a two-dimensional fault distribution; performs depth-weighted three-dimensional density structure inversion on the regional aerial gravity field data to obtain a three-dimensional density structure of an underground half space; performs three-dimensional boundary extraction on the three-dimensional density structure obtained by inverting the regional aerial gravity field data to obtain a three-dimensional fault structure distribution in the study area in the earth's crust; and comprehensively analyzes to obtain a scheme strategy for delineating and dividing faults in a land-sea transition zone based on aerial gravity data, thereby providing support for the division of fault structures in a land-sea transition zone, urban safety evaluation, and the like.

Description

Technical Field

[0001] The invention relates to the fields of geophysics and geochemistry, and in particular to a method for delineating fracture structures in a sea-land transition zone based on aerial gravity data. Background Art

[0002] The land-sea transition zone is a geographical unit where the ocean and land are compounded and intersected, and has the most superior resources and environmental conditions. The fault structure of the land-sea transition zone is often closely related to submarine landslides and earthquakes, and is a key parameter for the evaluation of land-sea geological disasters, and is of great significance to the stability evaluation of coastal cities.

[0003] At present, the fault structure in the sea is mainly identified and delineated by geophysical methods such as marine earthquakes and ship weights. At the same time, marine earthquake and ship weight measurements are mainly based on survey lines; the fault structure on land is often identified and delineated by geophysical methods such as ground gravity, earthquakes, and magnetotellurics, each with its own advantages and disadvantages. However, due to the influence of factors such as the depth of the intertidal zone and the sea-land transition zone, the geophysical and land geophysical measurements in the sea cannot carry out effective geophysical exploration in the "sea-land transition zone", which is the so-called "cannot come up from the sea and cannot go down to the land". The discontinuity of geophysical data in the sea-land transition zone restricts the research on basic geological issues such as the connection and extension of land faults in the sea.

[0004] In addition, satellite gravity can achieve full-area measurement and identify fault structures in the land-sea transition zone, but the resolution of satellite gravity is poor, and it can only preliminarily delineate some large fault structures, and cannot meet the needs of identifying fine fault structures. Therefore, a higher-resolution and more sophisticated geophysical exploration method is needed to meet the needs of delineating fault structures in the land-sea transition zone. Summary of the invention

[0005] In view of this, in order to solve the technical difficulties in identifying and delineating the fault structure in the land-sea transition zone, the purpose of the present invention is to propose a method for delineating the fault structure in the land-sea transition zone based on aerial gravity data. The method uses aerial gravity data as the main data and integrates multi-source information such as data splicing, Moho depth calculation, two-dimensional fault identification, three-dimensional density structure inversion and three-dimensional fault extraction. Especially in the land-sea transition zone in eastern China, where the fault structure is well developed and there are many large cities, the technology proposed in the present invention can economically and efficiently delineate and divide different types of faults in the land-sea transition zone, providing good technical support for land and sea geological disaster assessment, urban stability assessment, etc.

[0006] Based on the above purpose, in a first aspect, the present invention provides a method for delineating fault structures in a sea-land transition zone based on aerial gravity data, the method comprising the following steps:

[0007] The aerial gravity data obtained from different sources are fused and spliced ​​to obtain the regional aerial gravity field data covering the study area;

[0008] The depth of the Moho surface is obtained by fusion and splicing of regional airborne gravity field data;

[0009] The regional airborne gravity field data is enhanced by using the notch analysis to extract the boundary coefficient and obtain the two-dimensional fracture distribution.

[0010] Perform depth-weighted three-dimensional density structure inversion on regional aerial gravity field data to obtain the three-dimensional density structure of underground half-space;

[0011] The three-dimensional fracture structure distribution in the study area within the crust is obtained by extracting the three-dimensional boundary of the three-dimensional density structure obtained by inverting the regional airborne gravity field data.

[0012] By comprehensively analyzing the information of Moho surface depth, two-dimensional fault distribution, three-dimensional density structure and three-dimensional fault structure distribution, a strategy for delineating and dividing faults in the sea-land transition zone based on aerial gravity data was obtained.

[0013] As a further solution of the present invention, the aerial gravity data exceeding the preset measurement height among the gravity data obtained from different sources is selected as the reference data Δg A The height of the benchmark data is selected as the height of the final fused data. The benchmark data are unified in height, benchmark and scale by upward continuation, univariate linear regression model and interpolation, and then fused and spliced ​​to obtain the regional aerial gravity field data covering the study area.

[0014] As a further solution of the present invention, when the height is unified by upward extension, the calculation formula is:

[0015]

[0016] Where, Δg Bup is the gravity anomaly Δg B The aerodynamic gravity anomaly after the upward extension; h is the upward extension height, i.e. Δg A With Δg B The difference between the data acquisition heights; u and v are the frequencies in the x and y axis directions respectively; FT[ ] represents Fourier transform; x, y, z0 are the position coordinates of the measuring point.

[0017] As a further solution of the present invention, when the benchmark is corrected using univariate linear regression, the calculation formula is:

[0018] Δg Bup (x,y)=α+α1Δg A (x,y)+β

[0019] The parameters α and β of the univariate linear regression are calculated by the least squares method, and the aviation gravity anomaly Δg after the benchmark unification is calculated. B校 ; Among them, the aviation gravity anomaly Δg B校 The calculation formula is:

[0020] Δg B校 (x,y)=α+α1Δg Bup (x,y)+β

[0021] Then, by interpolation, Δg A (x,y) and Δg B校 The (x, y) data grids are unified, that is, after the scale is unified, the data are directly spliced ​​automatically through the Grid Knitting data splicing module in the geophysical data processing software Oasis montaj, and the hybrid method or stitching method is selected to obtain the regional airborne gravity field Δg covering the study area. 区 (x,y).

[0022] As a further solution of the present invention, when obtaining the Moho depth, a variable density model is used to calculate the Moho depth of the study area. The calculation formula of the Moho depth is:

[0023]

[0024] Where FT[ ] is the Fourier transform, h(x) is the Moho interface fluctuation value, z′ is the given average depth of the interface, Δρ0 is the density difference between the surface stratum and the Moho surface, and μ is the variable density factor.

[0025] As a further solution of the present invention, when obtaining the two-dimensional fracture distribution, the boundary information of the regional aerogravity field data is enhanced by using the two-dimensional notch analysis, and the statistical invariant M2 of the aerogravity surface is calculated. Then, the M2 surface is used as data input for enhancement processing, the boundary coefficient is extracted, the structural boundary is depicted, and the two-dimensional fracture distribution is divided. The calculation formula of the statistical invariant M2 is:

[0026] M2=m 20 +m 02

[0027]

[0028] Among them, m 20 、m 02 is the second-order spectral moment of the regional aerogravity field surface, n 20 、n 02 、n 11 are the three elements of the second-order spectral moment of the M2 surface, where is to find the divergence of the M2 data, It is used to take is the inverse sign of , and Mλ is the result parameter that describes the fracture and boundary.

[0029] As a further solution of the present invention, when performing depth-weighted three-dimensional density structure inversion on the regional aerial gravity field data to obtain the three-dimensional density structure of the underground half space, the regional aerial gravity field data Δg 区 (x, y) to perform three-dimensional physical structure inversion to obtain the density distribution information of the field source body. The objective function of the three-dimensional depth weighted inversion of the airborne gravity field is:

[0030]

[0031] Among them, d is the observed data; m is the model space; G is the sensitivity matrix; Fitting functionals to the data; is the model stabilization functional, W d and W m are the weight matrices for data and model, respectively;

[0032] In the process of depth-weighted 3D density structure inversion, the data weighting function W d Set as a constant diagonal matrix, Z is the depth weighted matrix, and the expression of the depth weighted matrix Z is:

[0033]

[0034] Where Z j represents the depth of the jth grid, β is a positive number related to the magnetic field decay rate, and is selected as 3 in the three-dimensional magnetic field inversion;

[0035] The conjugate gradient method was used to write a MATLAB language program. Under the condition of deep weighted regularization, the nonlinear conjugate gradient method was used to find the optimal solution of the objective function and obtain the three-dimensional density structure of the underground half space.

[0036] As a further solution of the present invention, when extracting the three-dimensional boundary of the three-dimensional density structure, the extraction is based on the three-dimensional second-order spectral moment and the statistical invariant method, and the anisotropy of the invariant intensity parameter is calculated through the three-dimensional statistical invariant intensity parameter. The covariance matrix of the three-dimensional second-order spectral moment of the three-dimensional density structure is:

[0037]

[0038] Among them, C is the covariance matrix, K x , K y , K z is the partial derivative of magnetic susceptibility with respect to the x, y, and z directions, i = 1, 2, ... M α ,j=1,2,…N α ,k=1,2,…l α ;

[0039] m 200 ,m 020 ,m 002 ,m 110 ,m 101 ,m 011 are the six elements of the three-dimensional second-order spectral moment;

[0040] L α 、M α 、N α is the size of the three-dimensional sliding window;

[0041] Perform eigenvalue decomposition on the covariance matrix C to obtain the diagonal matrix A, where the diagonal matrix A is:

[0042]

[0043] Use the elements a1, a2, and a3 in the diagonal matrix A to construct the local strength 3 M2, identification of boundary locations;

[0044] The local strength 3 The M2 data volume is used as input. The three elements b1, b2, and b3 of the diagonal matrix obtained by calculating the covariance matrix of the three-dimensional second-order spectral moment of the three-dimensional density structure and the diagonal matrix A are used to characterize the three-dimensional boundary coefficients of the fracture and identify the three-dimensional boundaries in the data volume. The three-dimensional boundary coefficients characterizing the fracture are:

[0045]

[0046] in is the local intensity 3 The divergence of the M2 data, It is used to take The inverse sign of .

[0047] In another aspect of the present invention, a computer device is provided, including a memory and a processor, wherein a computer program is stored in the memory, and when the computer program is executed by the processor, any one of the above methods for delineating fault structures in the sea-land transition zone based on aerial gravity data according to the present invention is executed.

[0048] In another aspect of the present invention, a computer-readable storage medium is provided, storing computer program instructions, which, when executed, implement any of the above methods for delineating fault structures in the sea-land transition zone based on aerial gravity data according to the present invention.

[0049] Compared with the prior art, the method for delineating the fault structure of the sea-land transition zone based on aerial gravity data proposed in the present invention realizes efficient and accurate fault structure identification and delineation by comprehensively processing aerial gravity data, multiple physical parameters and multi-source information, and has the following beneficial effects:

[0050] 1. Improved the accuracy of fault structure identification. The method of delineating fault structures in the land-sea transition zone based on aerial gravity data of the present invention can significantly improve the accuracy of fault structure identification. Due to the limitations of measurement range and environmental conditions, traditional geophysical methods often have the problem of data discontinuity in the application of land-sea transition zones. The present invention overcomes this limitation through the unified splicing of aerial gravity data and high-resolution inversion technology, and can obtain continuous and accurate gravity anomaly data, thereby improving the accuracy of fault structure identification.

[0051] 2. Enhanced boundary recognition capability. The present invention adopts a new three-dimensional fracture extraction method, which can extract the three-dimensional distribution morphology of three-dimensional fractures and enhance the recognition capability of fracture boundaries, significantly improve the clarity and accuracy of fracture boundaries, reduce the interference of high-frequency noise, and make fracture extraction more precise and reliable.

[0052] 3. Improved spatial resolution. Under the depth-weighted regularization condition, the present invention uses the nonlinear conjugate gradient method to invert the three-dimensional density structure, which can effectively improve the spatial resolution and obtain the underground three-dimensional density distribution map. Compared with the traditional two-dimensional fracture distribution extraction method, the three-dimensional inversion technology and the three-dimensional fracture extraction technology can more accurately reflect the spatial structure and distribution of fractures and provide more comprehensive geological information.

[0053] 4. Ability to make comprehensive use of multi-source information. The present invention comprehensively analyzes multi-source information such as Moho depth, two-dimensional fault distribution, three-dimensional density structure and three-dimensional fault distribution, and formulates a fault delineation and division strategy based on aerial gravity data. This comprehensive processing method not only improves the efficiency of data utilization, but also can more comprehensively and accurately reflect the characteristics and distribution of fault structures, providing strong technical support for the study of fault structures in the sea-continent transition zone.

[0054] 5. Wide application scope and high economic benefits. The method of the present invention is particularly suitable for the land-sea transition zone in eastern China where fault structures are well developed and large cities are numerous. In these areas, traditional geophysical methods are difficult to implement effectively. However, the present invention can economically and efficiently delineate and divide different types of fault structures through comprehensive processing of aerial gravity data, providing good technical support for land-sea geological disaster assessment, urban stability assessment, etc.

[0055] Moreover, compared with traditional geophysical methods, the method of the present invention has significant advantages in economy and efficiency. Aerial gravity data has a wide coverage and high measurement accuracy, and can obtain gravity anomaly data over a large area in a short time, reducing the cost and time of field measurement. By comprehensively processing multi-source information, the method of the present invention can efficiently and accurately delineate fault structures and provide reliable technical support.

[0056] In summary, the method of delineating the fault structure of the land-sea transition zone based on aerial gravity data of the present invention has the beneficial effects of improving recognition accuracy, enhancing boundary recognition ability, improving spatial resolution, comprehensively utilizing multi-source information, having a wide range of applications, and being economical and efficient, etc., which provides strong technical support and decision-making basis for the research and application of the fault structure of the land-sea transition zone.

[0057] These and other aspects of the present application will be more concise and understandable in the following description of the embodiments. It should be understood that the above general description and the following detailed description are only exemplary and explanatory and cannot limit the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention. For ordinary technicians in this field, other embodiments can be obtained based on these drawings without paying creative work.

[0059] In the figure:

[0060] Figure 1 The present invention is a flowchart of a method for delineating fault structures in a land-sea transition zone based on aerial gravity data according to an embodiment of the present invention. DETAILED DESCRIPTION

[0061] Below, the present application is further described in conjunction with the accompanying drawings and specific implementation methods. It should be noted that, under the premise of no conflict, the various embodiments or technical features described below can be arbitrarily combined to form a new embodiment.

[0062] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the embodiments of the present invention are further described in detail below in combination with specific embodiments and with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0063] At present, the fault structure in the sea is mainly identified and delineated by geophysical methods such as marine earthquakes and ship weights. At the same time, marine earthquake and ship weight measurements are mainly based on survey lines; the fault structure on land is often identified and delineated by geophysical methods such as ground gravity, earthquakes, and magnetotellurics, each with its own advantages and disadvantages. The discontinuity of geophysical data in the sea-land transition zone restricts the research on basic geological issues such as the connection and extension of land faults in the sea. The resolution of satellite gravity is poor, and only some large faults can be delineated.

[0064] In view of this, the present invention proposes a method for delineating the fault structure of the land-sea transition zone based on aerial gravity data. It is a comprehensive processing solution with aerial gravity data as the main data, integrating multi-source information such as data splicing, Moho surface fluctuation, two-dimensional fault identification, density structure and three-dimensional distribution of faults. Especially in the land-sea transition zone in eastern China, where the fault structure is well developed and there are many large cities, the technology proposed in the present invention can economically and efficiently delineate and divide different types of faults in the land-sea transition zone, providing good technical support for land and sea geological disaster assessment, urban stability assessment, etc.

[0065] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0066] The flowcharts shown in the accompanying drawings are only examples and do not necessarily include all the contents and operations / steps, nor must they be executed in the order described. For example, some operations / steps may also be decomposed, combined or partially merged, so the actual execution order may change according to actual conditions.

[0067] See also Figure 1 As shown, an embodiment of the present invention provides a method for delineating fault structures in a sea-land transition zone based on aerial gravity data, the method comprising the following steps:

[0068] Step S10, fusing and splicing the aerial gravity data obtained from different sources to obtain regional aerial gravity field data covering the study area;

[0069] Step S20, obtaining the Moho depth by using the regional aerial gravity field data obtained by fusion and splicing;

[0070] Step S30, using notch analysis to enhance the boundary information of the regional aerogravity field data, extract boundary coefficients, and obtain two-dimensional fracture distribution;

[0071] Step S40, performing depth-weighted three-dimensional density structure inversion on the regional aerial gravity field data to obtain the three-dimensional density structure of the underground half space;

[0072] Step S50, extracting the three-dimensional boundary of the three-dimensional density structure obtained by inverting the regional aerospace gravity field data, and obtaining the three-dimensional fracture structure distribution in the study area in the crust;

[0073] Step S60, by comprehensively analyzing the Moho depth, two-dimensional fault distribution, three-dimensional density structure and three-dimensional fault structure distribution, a strategy for delineating and dividing faults in the sea-land transition zone based on airborne gravity data is obtained.

[0074] The method for delineating the fault structure of the sea-land transition zone based on aerial gravity data in an embodiment of the present invention mainly includes obtaining regional aerial gravity anomalies covering the study area; obtaining the depth of the Moho surface; obtaining two-dimensional fault distribution; obtaining three-dimensional density structure of the underground half space; obtaining three-dimensional fault distribution in the crust, fault delineation and division strategy schemes, etc. The implementation steps include: 1) After unifying the height, benchmark and scale of the data, the regional aerogravity anomalies covering the study area are obtained by splicing; 2) The depth of the Moho surface is obtained by using the spliced ​​aerogravity field data; 3) The two-dimensional fault distribution is extracted using the notch analysis technology; 4) The depth-weighted three-dimensional density structure inversion of the regional aerogravity anomaly is performed, and the program is implemented using the MATLAB language. Under the condition of depth-weighted regularization, the nonlinear conjugate gradient method is used to cut the optimal solution of the objective function to obtain the underground three-dimensional density distribution; 5) For the three-dimensional density structure obtained by inverting the aerogravity data, a new three-dimensional fault extraction method is used to obtain the three-dimensional fault distribution in the crust. This method can extract faults from the three-dimensional data body, has the advantages of enhancing boundaries and reducing high-frequency noise, and significantly improves the vertical resolution; 6) On the basis of the four steps from step 2) to step 5), through the comprehensive utilization and analysis of the information such as the depth of the Moho surface, the two-dimensional fault structure distribution, the three-dimensional density structure and the three-dimensional fault distribution, the scheme strategy for the delineation and division of the faults in the sea-land transition zone based on the aerogravity data is obtained. The present invention is based on flight weight data and integrates a variety of physical parameters to provide support for the division of fault structures in the sea-land transition zone and urban safety evaluation.

[0075] Since aerial gravity data sometimes cannot completely cover the target study area, in order to better solve the fault structure problem of the land-sea transition zone, it is often necessary to merge and splice gravity data from different sources to form gravity field data covering the study area. Gravity data from different sources often have different scales, different collection heights and different benchmarks. In this embodiment, in step S10, the aerial gravity data with a higher measurement height among the gravity data obtained from different sources is selected as the benchmark data Δg AThe height of the benchmark data is selected as the height of the final fused data. The benchmark data are unified in height, benchmark and scale by upward continuation, univariate linear regression model and interpolation, and then fused and spliced ​​to obtain the regional aerial gravity field data covering the study area.

[0076] First, when the height is unified by upward extension, the calculation formula is:

[0077]

[0078] Where, Δg Bup is the gravity anomaly Δg B The aerodynamic gravity anomaly after the upward extension; h is the upward extension height, i.e. Δg A With Δg B The difference between the data acquisition heights; u and v are the frequencies in the x and y axis directions respectively; FT[ ] represents Fourier transform; x, y, z0 are the position coordinates of the measuring point.

[0079] Secondly, when the benchmark is corrected using univariate linear regression, the calculation formula is:

[0080] Δg Bup (x,y)=α+α1Δg A (x,y)+β

[0081] The parameters α and β of the univariate linear regression are calculated by the least squares method, and the aviation gravity anomaly Δg after the benchmark unification is calculated. B校 ; Among them, the aviation gravity anomaly Δg B校 The calculation formula is:

[0082] Δg B校 (x,y)=α+α1Δg Bup (x,y)+β

[0083] Then, by interpolation, Δg A (x,y) and Δg B校 The (x, y) data grids are unified, that is, after the scale is unified, the data are directly spliced ​​automatically through the Grid Knitting data splicing module in the geophysical data processing software Oasis montaj, and the hybrid method or stitching method is selected to obtain the regional airborne gravity field Δg covering the study area. 区 (x, y). The aerial gravity area data obtained in this step is the data basis of steps S20 to S50.

[0084] The Moho surface is the interface between the crust and the upper mantle, and is also an important density interface. The Moho surface uplift-depression junction zone is a deep transition zone between the regional orogenic belt and the stable block, and is often a regional plate junction zone or a deep fault development zone. The Bouguer gravity anomaly is closely related to the Moho surface. In this embodiment, in step S20, when obtaining the Moho surface depth, the variable density model is used to calculate the Moho surface depth of the study area. The calculation formula of the Moho surface depth is:

[0085]

[0086] Where FT[ ] is the Fourier transform, h(x) is the Moho interface fluctuation value, z′ is the given average depth of the interface, Δρ0 is the density difference between the surface stratum and the Moho surface, and μ is the variable density factor.

[0087] In this embodiment, in step S30, when obtaining the two-dimensional fracture distribution, the two-dimensional scratch analysis is used to perform boundary information enhancement processing on the regional aerogravity field data, and the statistical invariant M2 of the aerogravity surface is calculated. Then, the M2 surface is used as data input for enhancement processing, and the boundary coefficient is extracted to characterize the structural boundary and divide the two-dimensional fracture distribution. The calculation formula of the statistical invariant M2 is:

[0088] M2=m 20 +m 02

[0089]

[0090] Among them, m 20 、m 02 is the second-order spectral moment of the regional aerogravity field surface, n 20 、n 02 、n 11 are the three elements of the second-order spectral moment of the M2 surface, where is to find the divergence of the M2 data, It is used to take is the inverse sign of , and Mλ is the result parameter that describes the fracture and boundary.

[0091] In this embodiment, in step S40, when depth-weighted three-dimensional density structure inversion is performed on the regional aerial gravity field data to obtain the three-dimensional density structure of the underground half space, the regional aerial gravity field data Δg 区 (x, y) to perform three-dimensional physical structure inversion to obtain the density distribution information of the field source body. The objective function of the three-dimensional depth weighted inversion of the airborne gravity field is:

[0092]

[0093] Among them, d is the observed data; m is the model space; G is the sensitivity matrix; Fitting functionals to the data; is the model stabilization functional, W d and W m are the weight matrices for data and model, respectively;

[0094] In the actual depth-weighted 3D density structure inversion process, the data weighting function W d Set as a constant diagonal matrix, Z is the depth weighted matrix, and the expression of the depth weighted matrix Z is:

[0095]

[0096] Where Z j represents the depth of the jth grid, β is a positive number related to the magnetic field decay rate, and is selected as 3 in the three-dimensional magnetic field inversion;

[0097] Finally, the conjugate gradient method was used to write a MATLAB language program. Under the condition of deep weighted regularization, the nonlinear conjugate gradient method was used to find the optimal solution of the objective function and obtain the three-dimensional density structure of the underground half-space.

[0098] In this embodiment, in step S50, by performing 3D boundary extraction on the 3D density structure obtained by inverting the aerial weight data, the 3D fracture structure distribution in the study area can be obtained. The 3D fracture extraction method is a new method developed based on the 3D second-order spectral moment and its statistical invariant.

[0099] When extracting the three-dimensional boundary of the three-dimensional density structure, it is extracted based on the three-dimensional second-order spectral moment and statistical invariant method. The anisotropy of the invariant intensity parameter is calculated through the three-dimensional statistical invariant intensity parameter. The covariance matrix of the three-dimensional second-order spectral moment of the three-dimensional density structure is:

[0100]

[0101] Among them, C is the covariance matrix, K x , K y , K z is the partial derivative of magnetic susceptibility with respect to the x, y, and z directions, i = 1, 2, ... M α ,j=1,2,…N α ,k=1,2,…l α ;m 200 ,m 020 ,m 002 ,m 110 ,m 101 , m 011 are the six elements of the three-dimensional second-order spectral moment; L α 、M α 、N α is the size of the three-dimensional sliding window;

[0102] Perform eigenvalue decomposition on the covariance matrix C to obtain the diagonal matrix A, where the diagonal matrix A is:

[0103]

[0104] Use the elements a1, a2, and a3 in the diagonal matrix A to construct the local strength 3 M2, identification of boundary locations;

[0105] The local strength 3 The M2 data volume is used as input. The three elements b1, b2, and b3 of the diagonal matrix obtained by calculating the covariance matrix of the three-dimensional second-order spectral moment of the three-dimensional density structure and the diagonal matrix A are used to characterize the three-dimensional boundary coefficients of the fracture and identify the three-dimensional boundaries in the data volume. The three-dimensional boundary coefficients characterizing the fracture are:

[0106]

[0107] in is the local intensity 3 The divergence of the M2 data, It is used to take The new method of 3D fracture extraction based on 3D second-order spectral moment can not only enhance the weak boundary recognition ability and improve the lateral resolution, but also this local surface element statistical calculation method can effectively suppress high-frequency noise. At the same time, different from the traditional 2D boundary recognition method that can only obtain the plane boundary distribution, this method can obtain the 3D boundary distribution data body, which can not only characterize the plane distribution characteristics of the fracture at different depths, but also characterize the extension morphology of the fracture from shallow to deep.

[0108] In this embodiment, in step S60, the two-dimensional fracture and boundary distribution extracted in step S10 includes all the fractures and boundaries with gravity response in the deep and shallow parts of the study area. Generally, these lines with longer extensions and higher coefficient values ​​correspond to fractures; the steep edge bands of density on different depth slices of the three-dimensional density body obtained in step S20 correspond to fractures; the lines with longer extensions and higher coefficient values ​​on different depth slices of the three-dimensional fracture data body obtained in step S50 correspond to fractures.

[0109] Based on the Moho depth calculated in step S20, the two-dimensional fault and boundary distribution extracted in step S30, the three-dimensional density structure obtained by depth weighted inversion in step S40, and the three-dimensional fault distribution obtained in step S50, a delineation and division scheme for four types of faults in the sea-continent transition zone, namely, general faults, large faults, deep large faults (super-shell faults), and inferred faults, is established. The scheme is as follows:

[0110] (1) The two-dimensional fault and boundary distribution extracted in step S30 are used to delineate the plane distribution of all faults in the study area, including deep and shallow parts; (2) The depth of the Moho surface obtained by inversion is averaged and divided into three parts, which are respectively used as the depth of the upper crust, the depth of the middle crust and the depth of the lower crust; for example, the approximate depth of the Moho surface in the eastern sea-continent transition zone of my country is 30km, so by averaging this depth and dividing it into three parts, it can be roughly considered that the depth of the upper crust is 0-10km; the depth of the middle crust is 10-20km; the depth of the lower crust is 20-30km; (3) The three-dimensional density data body and the three-dimensional fault structure body are sliced ​​horizontally at the depths of the upper crust, the middle crust and the lower crust. For example, the inversion depth of the three-dimensional density body and the three-dimensional fault structure body in the eastern sea-continent transition zone of my country is 35km, and horizontal slices can be made every 5km, so that slices at the depths of 5km and 10km in the upper crust, 15km and 20km in the middle crust, 25km and 30km in the lower crust, and 35km in depth can be obtained.

[0111] If the fault identified in (1) is not present on the geological map of the study area, the fault is considered to be a hidden fault. If the fault identified in (1) is only reflected at a shallow depth in the upper crust of the three-dimensional density data body and the three-dimensional fault structure body, but not in the depth slices of the middle and lower crust, the fault is considered to be a general fault. For example, in the eastern sea-continent transition zone of my country, the fault identified in (1) is only displayed on slices of the upper crust of the three-dimensional density body and the three-dimensional fault structure body at a depth of 5 km or 10 km, so the fault is considered to be a general fault.

[0112] If the fault identified in (1) still responds on the depth slices of the middle crust and lower crust in the three-dimensional density data body and the three-dimensional fault structure body, then the fault is considered to be a large fault; for example, in the eastern sea-continent transition zone of my country, the fault identified in (1) is only displayed on the slices of the upper crust at a depth of 15km to 30km in the three-dimensional density body and the three-dimensional fault structure body, then the fault is considered to be a large fault.

[0113] If the fault identified in (1) is reflected in slices with depths of the upper, middle and lower crust, and also in slices with depths greater than the lower crust, then the fault is considered to be a deep fault (transcrustal fault). For example, in the eastern sea-continent transition zone of my country, the fault identified in (1) is only reflected in slices with depths of more than 30 km in the upper crust of the three-dimensional density body and the three-dimensional fault structure body, then the fault is considered to be a deep fault.

[0114] The method for delineating the fracture structure of the sea-land transition zone based on aerial gravity data of the present invention realizes efficient and accurate fracture structure identification and delineation by comprehensively processing aerial gravity data, multiple physical parameters and multi-source information, improves the fracture structure identification accuracy, enhances the boundary identification capability, improves the spatial resolution, and can comprehensively utilize multi-source information. The method for delineating the fracture structure of the sea-land transition zone based on aerial gravity data of the present invention can significantly improve the accuracy of fracture structure identification. Due to the limitations of measurement range and environmental conditions, traditional geophysical methods often have the problem of data discontinuity in the application of sea-land transition zones. The present invention overcomes this limitation through unified splicing of aerial gravity data, high-resolution inversion and three-dimensional fracture extraction technology, and can obtain continuous and accurate gravity anomaly data, thereby improving the accuracy of fracture structure identification.

[0115] The present invention adopts a new three-dimensional fracture extraction method, which can extract the three-dimensional distribution morphology of three-dimensional fractures and enhance the recognition ability of fracture boundaries, significantly improve the clarity and accuracy of fracture boundaries, reduce the interference of high-frequency noise, and make fracture extraction more precise and reliable. Under the condition of depth-weighted regularization, the present invention uses the nonlinear conjugate gradient method to invert the three-dimensional density structure, which can effectively improve the spatial resolution and obtain the underground three-dimensional density distribution map. Compared with the traditional two-dimensional fracture distribution extraction method, the three-dimensional inversion technology and the three-dimensional fracture extraction technology can more accurately reflect the spatial structure and distribution of fractures and provide more comprehensive geological information.

[0116] The present invention conducts a comprehensive analysis of multi-source information such as Moho depth, two-dimensional fault distribution, three-dimensional density structure and three-dimensional fault distribution, and formulates a fault delineation and division strategy based on aerial gravity data. This comprehensive processing method not only improves the efficiency of data utilization, but also can more comprehensively and accurately reflect the characteristics and distribution of fault structures, providing strong technical support for the study of fault structures in the sea-land transition zone. The method of the present invention is particularly suitable for the sea-land transition zone in eastern China where fault structures are well developed and there are many large cities. In these areas, traditional geophysical methods are difficult to implement effectively, and the present invention can economically and efficiently delineate and divide different types of fault structures through the comprehensive processing of aerial gravity data, providing good technical support for sea and land geological disaster assessment, urban stability assessment, etc.

[0117] Moreover, compared with traditional geophysical methods, the method of the present invention has significant advantages in economy and efficiency. Aerial gravity data has a wide coverage and high measurement accuracy, and can obtain gravity anomaly data over a large area in a short time, reducing the cost and time of field measurement. By comprehensively processing multi-source information, the method of the present invention can efficiently and accurately delineate fault structures and provide reliable technical support.

[0118] Therefore, the method of delineating the fault structure of the land-sea transition zone based on aerial gravity data in the present invention provides strong technical support and decision-making basis for the research and application of the fault structure of the land-sea transition zone by improving recognition accuracy, enhancing boundary recognition capability, improving spatial resolution, comprehensively utilizing multi-source information, having a wide range of applications, and being economical and efficient.

[0119] According to a second aspect of an embodiment of the present invention, a computer device is provided, including a memory and a processor, wherein a computer program is stored in the memory, and when the computer program is executed by the processor, the method of any one of the above embodiments is implemented.

[0120] The computer device includes a processor and a memory, and may also include an input system and an output system. The processor, memory, input system and output system may be connected via a bus or other means, and the input system may receive input digital or character information, and generate signal input related to the migration of the method for delineating the fracture structure of the sea-land transition zone based on airborne gravity data. The output system may include display devices such as a display screen.

[0121] As a non-volatile computer-readable storage medium, the memory can be used to store non-volatile software programs, non-volatile computer executable programs and modules, such as the program instructions / modules corresponding to the method for delineating the fracture structure of the sea-land transition zone based on aerial gravity data in the embodiment of the present application. The memory may include a program storage area and a data storage area, wherein the program storage area may store an operating system and an application required by at least one function; the data storage area may store data created by the use of the method for delineating the fracture structure of the sea-land transition zone based on aerial gravity data, etc. In addition, the memory may include a high-speed random access memory, and may also include a non-volatile memory, such as at least one disk storage device, a flash memory device, or other non-volatile solid-state storage device. In some embodiments, the memory may optionally include a memory remotely arranged relative to the processor, and these remote memories may be connected to the local module via a network. Examples of the above-mentioned network include, but are not limited to, the Internet, an intranet, a local area network, a mobile communication network, and combinations thereof.

[0122] The processor may be a central processing unit (CPU), a controller, a microcontroller, a microprocessor, or other data processing chip in some embodiments. The processor is generally used to control the overall operation of the computer device. In this embodiment, the processor is used to run the program code stored in the memory or process data. The processors of the multiple computer devices of the computer device of this embodiment execute various functional applications and data processing of the server by running the non-volatile software programs, instructions and modules stored in the memory, that is, the steps of the method for delineating the fracture structure of the sea-land transition zone based on aerial gravity data of the above method embodiment are implemented.

[0123] It should be understood that, to the extent that they do not conflict with each other, all the embodiments, features and advantages described above for the method for delineating fault structures in the land-sea transition zone based on aerial gravity data according to the present invention are also applicable to the method and storage medium for delineating fault structures in the land-sea transition zone based on aerial gravity data according to the present invention.

[0124] It will also be appreciated by those skilled in the art that various exemplary logic blocks, modules, circuits and algorithm steps described in conjunction with the disclosure herein can be implemented as electronic hardware, computer software or a combination of the two. In order to clearly illustrate this interchangeability of hardware and software, a general description has been given to the functions of various schematic components, blocks, modules, circuits and steps. Whether this function is implemented as software or hardware depends on specific applications and the design constraints imposed on the entire system. Those skilled in the art can implement the function in various ways for each specific application, but this implementation decision should not be interpreted as causing a departure from the disclosed scope of the embodiments of the present invention.

[0125] The above are exemplary embodiments disclosed in the present invention, but it should be noted that various changes and modifications may be made without departing from the scope disclosed in the embodiments of the present invention as defined in the claims. The functions, steps and / or actions of the method claims according to the disclosed embodiments described herein do not need to be performed in any particular order. In addition, although the elements disclosed in the embodiments of the present invention may be described or required in individual form, they may also be understood as multiple unless explicitly limited to the singular.

[0126] It should be understood that, as used herein, the singular form "a" or "an" is intended to include the plural form as well, unless the context clearly supports an exception. It should also be understood that, as used herein, "and / or" refers to any and all possible combinations of one or more of the items listed in association. The serial numbers of the embodiments disclosed in the above embodiments of the present invention are for description only and do not represent the advantages and disadvantages of the embodiments.

[0127] A person skilled in the art should understand that the discussion of any of the above embodiments is only exemplary and is not intended to imply that the scope of the disclosure of the embodiments of the present invention (including the claims) is limited to these examples; under the concept of the embodiments of the present invention, the technical features in the above embodiments or different embodiments can also be combined, and there are many other changes in different aspects of the embodiments of the present invention as above, which are not provided in detail for the sake of simplicity. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the embodiments of the present invention should be included in the protection scope of the embodiments of the present invention.

Claims

1. A method for delineating fault structures in the sea-land transition zone based on aerial gravity data, characterized in that: The following steps are involved: The aerial gravity data obtained from different sources are fused and spliced ​​to obtain the regional aerial gravity field data covering the study area; The depth of the Moho surface is obtained by fusion and splicing of regional airborne gravity field data; The regional airborne gravity field data is enhanced by using the notch analysis to extract the boundary coefficient and obtain the two-dimensional fracture distribution. Perform depth-weighted three-dimensional density structure inversion on regional aerial gravity field data to obtain the three-dimensional density structure of underground half-space; The three-dimensional fault structure distribution in the study area within the crust is obtained by extracting the three-dimensional boundary of the three-dimensional density structure obtained by inverting the regional airborne gravity field data; By comprehensively analyzing the Moho depth, two-dimensional fault distribution, three-dimensional density structure and three-dimensional fault structure distribution, a strategy for delineating and dividing faults in the sea-land transition zone based on airborne gravity data was obtained. Among them, the depth-weighted three-dimensional density structure inversion of the regional aerial gravity field data is performed to obtain the three-dimensional density structure of the underground half-space. 区 (x, y) to perform three-dimensional physical structure inversion to obtain the density distribution information of the field source body. The objective function of the three-dimensional depth weighted inversion of the airborne gravity field is: Among them, d is the observed data; m is the model space; G is the sensitivity matrix; Fitting functionals to the data; is the model stabilization functional, W d and W m are the weight matrices for data and model, respectively; In the process of depth-weighted 3D density structure inversion, the data weighting function W d Set as a constant diagonal matrix, Z is the depth weighted matrix, and the expression of the depth weighted matrix Z is: Where Z j represents the depth of the jth grid, β is a positive number related to the magnetic field decay rate, and is selected as 3 in the three-dimensional magnetic field inversion; The conjugate gradient method was used to write a MATLAB language program. Under the condition of deep weighted regularization, the nonlinear conjugate gradient method was used to find the optimal solution of the objective function and obtain the three-dimensional density structure of the underground half space.

2. The method for delineating fault structures in the sea-land transition zone based on aerial gravity data according to claim 1 is characterized in that: The aerial gravity data obtained from different sources that exceeds the preset measurement height is selected as the reference data Δg A The height of the benchmark data is selected as the height of the final fused data. The benchmark data are unified in height, benchmark and scale by upward continuation, univariate linear regression model and interpolation, and then fused and spliced ​​to obtain the regional aerial gravity field data covering the study area.

3. The method for delineating fault structures in the sea-land transition zone based on aerial gravity data according to claim 2 is characterized in that: When the height is unified by upward extension, the calculation formula is: Where, Δg Bup is the gravity anomaly Δg B The aerodynamic gravity anomaly after the upward extension; h is the upward extension height, i.e. Δg A With Δg B The difference between the data collection heights; u and v are the frequencies in the x and y axis directions respectively; FT[] represents Fourier transform; x, y, z0 are the position coordinates of the measuring point.

4. The method for delineating fault structures in the sea-land transition zone based on aerial gravity data according to claim 3 is characterized in that: When the benchmark is corrected using univariate linear regression, the calculation formula is: Δg Bup (x,y)=α+α1Δg A (x,y)+β The parameters α and β of the univariate linear regression are calculated by the least squares method, and the aviation gravity anomaly Δg after the benchmark unification is calculated. B校 ; Among them, the aviation gravity anomaly Δg B校 The calculation formula is: Δg B校 (x,y)=α+α1Δg Bup (x,y)+β Then, by interpolation, Δg A (x,y) and Δg B校 The (x, y) data grids are unified, that is, after the scale is unified, the data are directly spliced ​​automatically through the Grid Knitting data splicing module in the geophysical data processing software Oasis montaj, and the hybrid method or stitching method is selected to obtain the regional airborne gravity field Δg covering the study area. 区 (x,y).

5. The method for delineating fault structures in the sea-land transition zone based on aerial gravity data according to claim 4 is characterized in that: When obtaining the Moho depth, the variable density model is used to calculate the Moho depth of the study area. The calculation formula for the Moho depth is: Where FT[] is the Fourier transform, h(x) is the Moho interface fluctuation value, z′ is the given average depth of the interface, Δρ0 is the density difference between the surface stratum and the Moho surface, and μ is the variable density factor.

6. The method for delineating fault structures in the sea-land transition zone based on aerial gravity data according to claim 5 is characterized in that: When obtaining the two-dimensional fault distribution, the two-dimensional notch analysis is used to enhance the boundary information of the regional aerogravity field data, and the statistical invariant M2 of the aerogravity surface is calculated. The M2 surface is then used as data input for enhancement processing, the boundary coefficient is extracted, the structural boundary is depicted, and the two-dimensional distribution of the fault is divided. The calculation formula of the statistical invariant M2 is: M2=m 20 +m 02 Among them, m 20 、m 02 is the second-order spectral moment of the regional aerogravity field surface, n 20 、n 02 、n 11 are the three elements of the second-order spectral moment of the M2 surface, where is to find the divergence of the M2 data, It is used to take is the inverse sign of , and Mλ is the result parameter that describes the fracture and boundary.

7. The method for delineating fault structures in the sea-land transition zone based on aerial gravity data according to claim 6 is characterized in that: When extracting the three-dimensional boundary of the three-dimensional density structure, it is extracted based on the three-dimensional second-order spectral moment and statistical invariant method. The anisotropy of the invariant intensity parameter is calculated through the three-dimensional statistical invariant intensity parameter. The covariance matrix of the three-dimensional second-order spectral moment of the three-dimensional density structure is: Among them, C is the covariance matrix, K x , K y , K z is the partial derivative of magnetic susceptibility with respect to the x, y, and z directions, i = 1, 2, ... M α ,j=1,2,…N α ,k=1,2,…l α ;m 200 ,m 020 ,m 002 ,m 110 ,m 101 , m 011 are the six elements of the three-dimensional second-order spectral moment; L α 、M α 、N α is the size of the three-dimensional sliding window; Perform eigenvalue decomposition on the covariance matrix C to obtain the diagonal matrix A, where the diagonal matrix A is: Use the elements a1, a2, and a3 in the diagonal matrix A to construct the local strength 3 M2, identification of boundary locations; The local strength 3 The M2 data volume is used as input. The three elements b1, b2, and b3 of the diagonal matrix obtained by calculating the covariance matrix of the three-dimensional second-order spectral moment of the three-dimensional density structure and the diagonal matrix A are used to characterize the three-dimensional boundary coefficients of the fracture and identify the three-dimensional boundaries in the data volume. The three-dimensional boundary coefficients characterizing the fracture are: in is the local intensity 3 The divergence of the M2 data, It is used to take The inverse sign of .

8. A computer device, characterized in that: The method comprises a memory and a processor, wherein a computer program is stored in the memory, and when the computer program is executed by the processor, the method for delineating fault structures in a sea-land transition zone based on aerial gravity data as described in any one of claims 1 to 7 is executed.

9. A computer-readable storage medium storing computer program instructions, which, when executed, implement the method for delineating fault structures in the sea-land transition zone based on aerial gravity data as described in any one of claims 1 to 7.

Images