CN110244352B - Variable density-based crust thickness gravity inversion method - Google Patents
Variable density-based crust thickness gravity inversion method Download PDFInfo
- Publication number
- CN110244352B CN110244352B CN201910502882.4A CN201910502882A CN110244352B CN 110244352 B CN110244352 B CN 110244352B CN 201910502882 A CN201910502882 A CN 201910502882A CN 110244352 B CN110244352 B CN 110244352B
- Authority
- CN
- China
- Prior art keywords
- data
- density
- thickness
- crust
- gravity anomaly
- Prior art date
- Legal status (The legal status 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 status listed.)
- Expired - Fee Related
Links
- 230000005484 gravity Effects 0.000 title claims abstract description 143
- 238000000034 method Methods 0.000 title claims abstract description 34
- 230000008859 change Effects 0.000 claims abstract description 60
- 230000008021 deposition Effects 0.000 claims description 39
- 238000001914 filtration Methods 0.000 claims description 15
- 230000009466 transformation Effects 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000001228 spectrum Methods 0.000 claims description 5
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 5
- 230000005856 abnormality Effects 0.000 claims description 4
- 238000005553 drilling Methods 0.000 claims description 4
- 230000000007 visual effect Effects 0.000 claims description 4
- 239000004744 fabric Substances 0.000 claims description 3
- 239000013049 sediment Substances 0.000 description 7
- 230000006870 function Effects 0.000 description 6
- 239000000463 material Substances 0.000 description 5
- 241001417093 Moridae Species 0.000 description 2
- 206010027977 Morose Diseases 0.000 description 2
- 238000004097 X-ray Buerger Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 241000227135 Atrina Species 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000009189 diving Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000007620 mathematical function Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a variable density-based CRUST thickness gravity inversion method, which comprises the steps of obtaining CRUST model CRUST1.0 CRUST thickness and density data and upper mantle density data, wherein the global space resolution is 1 degree multiplied by 1 degree; establishing a shell mantle density difference which changes along with the horizontal position according to the density and the thickness of the upper, middle and lower shells of the shell model CRUST1.0 and the density data of the upper mantle; acquiring gravity anomaly caused by fluctuation change of the Mohuo surface from the Buge gravity anomaly; and (3) inverting the depth of the Mohuo surface by the gravity anomaly caused by the fluctuation change of the Mohuo surface and obtaining the thickness of the earth crust. Compared with the prior art, the method and the device have the advantages that the disclosed data are used as the basis, the practical shell mantle density model is established, the method and the device are easy to technically realize, the inversion result is more in line with the actual geological characteristics, and the accuracy is higher.
Description
Technical Field
The invention relates to the field of geoscience, in particular to a crust thickness gravity inversion method based on variable density.
Background
The change of the earth gravity field is related to the uneven distribution of the density of the materials from the vicinity of the earth surface to the deep part of the earth, and the density distribution of the materials is the reflection of the structure and the change of the earth system, so that the observation of the earth gravity field and the study of the gravity change (namely gravity anomaly) caused by the uneven distribution of the density of the underground materials have important scientific significance for understanding the structure of the earth layer, the dynamic process of the earth, the resources and the environmental change. In addition, the gravity observation is simple, the data acquisition is convenient, and relatively speaking, the gravity exploration has the advantages of economy, large exploration depth, rapid acquisition of information on the area and the like, so that the gravity exploration is widely applied.
The gravity anomaly inversion has important values for deep structure research, resource and energy investigation, earthquake prediction research and the like. According to the calculated parameters. In contrast, the inversion of gravity anomaly is divided into inversion of physical properties and inversion of physical properties. The density distribution of the field source body is determined according to the density distribution; the latter is based on given density, the characteristic point coordinates of the underground space subdivision unit are inverted to obtain the fluctuation form of the upper and lower interfaces of the geologic body. In recent years, inversion of the depth of the mojojoba surface by using gravity anomaly to further obtain the thickness of the earth crust has become a main research content of gravity inversion.
When the depth of the Mohuo surface and the thickness of the crust are inverted by utilizing gravity anomaly, the density difference between the crust and the mantle is a key parameter. Currently, most scholars perform inversion using a constant density mode. Liu Zuihui et al (1983) when inversion of depth of Mohuo surface in south sea area is carried out by using gravity anomaly, average density of crusta is 2.67g/cm3The material in the upper mantle is 3.27g/cm3The density difference of the instant mantle is-0.6 g/cm3. Xudeqiong and Jiangjiazhen (1989) taking crusta with average density of 2.73g/cm3(Shell mantle density difference is-0.54 g/cm3) The depth of the surface of the Mohuo in the North of the south sea was inverted, Tsai et al (2005), and the average densities of the crust and the upper mantle, respectively, at the time of inverting the surface of the Mohuo in the North of the south sea were 2.84g/cm3And 3.28g/cm3The density difference of the shell mantle is-0.44 g/cm3The density values are also adopted to respectively invert the depths of the mohuo surfaces of the south sea and the mohuo surfaces of the east sea, the south sea and the adjacent areas thereof, and the density difference of the shell taking mantle of the atrina hutch and the like (2008) is-0.30 g/cm3The depths of the mohoheds in the northeast of the south sea were inverted. Subsequently, when many scholars invert the depth of the Mohuo surface, the density difference of most shell-taking mantle is-0.50 g/cm3(Liu Jian Hua, 1993; Wang et al, 2017; Wu Bian Ma et al, 2017; Unnikrishnan et al, 2018) which differ only in that the average density of the crust and the average density of the upper mantle differ, e.g. the average density of the crust and the upper mantle taken from Liu Jian Hua (1993) are 2.82g/cm3And 3.32g/cm3The average densities of the crust and upper mantle, taken from Unnikrishnan et al (2018), were 2.80g/cm, respectively3And 3.30g/cm3. In addition to the above values, Trungent (2004) gave a mantle density difference of-0.40 g/cm when inverting the depth of Mohuo surface in south China sea3Korea and so on (2014) also choose the density difference to be-0.40 g/cm3The depths of the Mohuo surfaces of Taiwan and adjacent regions are inverted. Liand Wang (2016) has shell density difference of-0.32 g/cm3The mohol face depth was inverted for the east and southeast regions of asia.
Although the above constant density values solve the inversion problem of the mohuo surface depth and the crust thickness of the research area to a certain extent, when the area to be inverted is large, the difference of the crust density cannot be accurately depicted by only using one constant value. To this end, some scholars use mathematical functions to model the variation of density with depth and refer to the mojohn plane depth inversion. For example, Feng Juan et al (2014) invert the surface depth of the muhols in North China by using an index density-depth function, while Chenyan et al (2017) invert the surface depth of the muhols in the North Suzhou basin by using an index density-depth function; zhang En et al (2015) inverted the depth of the Mohuo surface in Chuandian by using a parabolic density-depth function. In addition to these 2 density-depth functions, other density-depth functions (such as polynomial) can also be used for the mojohn plane depth inversion, but no relevant research has been found yet. In addition, when the depth of the surface of the moore in the sea area is inverted, in order to improve the inversion accuracy, the thermal disturbance gravity anomaly can be corrected before inversion (Liand Wang, 2016; Wu Zhao et al, 2017; Liu shou et al, 2018) to obtain more accurate surface gravity anomaly of the moore.
The above variable density function takes into account the variation of the shell mantle density difference along the depth, which is more accurate than the constant density. However, the density changes only in the vertical direction, and when the structure of the research area is complex (such as the diving zone, the central ridge and other areas), the change of the density in the horizontal direction cannot be ignored, so that when the structure of the crust in the middle of the south China sea basin is explained by using gravity, such as Yangguanyu (2011), the section explanation is performed by using the combination of the density changes in the horizontal direction. And Liujun and the like (2009) utilize data such as earthquake profiles, sonar buoys and the like as control points, and adopt transverse density change to perform gravity profile fitting to obtain a crust structure in the north of the south China sea. Besides sea areas, the gravity profile fitting by using the transverse density change is also a main research means in the study of the land crust structure. ZyongYong et al (2014) establish crust density structures in the middle and lower reaches of Yangtze river and adjacent areas; the deep crustal structure of the eastern part of Ordos, Zhongqinling, Sichuan was studied by the Wangqian et al (2015); zhang et al (2018) studied the deep crustal structure of the great igneous province of the mount of the Emei mountain by density fitting of the Lijiang-Guiyang gravity profile.
The density difference changing along the section better reflects the transverse change characteristic of the shell mantle density difference, is more consistent with the actual geological characteristic, but is only a two-dimensional result and cannot present the regional characteristics of the variation of the depth of the muhuo surface and the thickness of the crust. In addition, the process of profile gravity inversion is essentially to continuously adjust a model to perform forward fitting, and the key point is to establish an initial model by using more accurate prior information (such as a deep reflection seismic profile and the like), which is difficult to realize in three-dimensional inversion, and to establish a variable shell mantle density difference model by other methods so as to improve the accuracy of inversion results of regional mohuo surface depth and crust thickness, and the method has important practical significance for deep structure research, regional structure research and the like.
Disclosure of Invention
The invention aims to solve the defects in the prior art, and provides a crust thickness gravity inversion method based on variable density.
In order to achieve the purpose, the invention is implemented according to the following technical scheme:
a variable density-based crust thickness gravity inversion method comprises the following steps:
s1, acquiring crustal model CRUST1.0 crustal thickness and density data and upper mantle density data with the global coverage spatial resolution of 1 degree multiplied by 1 degree;
s2, establishing a shell mantle density difference which changes along with the horizontal position according to the density and thickness of the upper, middle and lower earth shells and the density data of the upper mantle of the earth shell model CRUST 1.0;
s3, acquiring gravity anomaly caused by fluctuation change of the Mohuo surface due to the Buge gravity anomaly;
and S4, inverting the depth of the Mohuo surface by the gravity anomaly caused by the fluctuation change of the Mohuo surface and obtaining the thickness of the earth crust.
Further, the specific step of S1 is:
s11, opening a global CRUST model CRUST1.0 database with the spatial resolution of 1 degree multiplied by 1 degree, downloading a CRUST1.0 data file CRUST1.0.tar.gz in a Download selection, wherein the file CRUST1.0.tar.gz comprises codes written in FORTRAN language and is used for extracting CRUST thickness and density data and upper mantle density data;
s12, running and downloading a getCN1xyz.f file in the file by using a Microsoft Visual Studio platform, reading a thickness file of 8 density layers (a water layer, an ice layer, an upper deposition layer, a middle deposition layer, a lower deposition layer, an upper crust, a middle crust and a lower crust) covering the world and density data files of the 8 density layers and an upper mantle from the data file, wherein the thickness file covers the 8 density layers in latitude and longitude ranges of 179.5 degrees W-179.5 degrees E and 89.5 degrees S-89.5 degrees N;
and S13, extracting thickness and density data of each density layer in the research area range from the data file of S12 according to the specific latitude and longitude range of the research area, and arranging the thickness and density data of each density layer into two files of the thickness and density of each density layer.
Further, the specific step of S2 is:
s21, carrying out projection transformation on the thickness and density data of each density layer arranged in the S13, and projecting the geographic coordinates into plane rectangular coordinates;
s22, calculating the shell mantle density difference of each point at different horizontal positions in the research area by using the following formula:
where i denotes different horizontal position points, Δ ρiIs the difference in shell density at the ith point,
the density values of the upper crust, the middle crust, the lower crust and the upper mantle at the ith point are respectively,
the thicknesses of the upper crust, the middle crust and the lower crust at the ith point are respectively;
s23 for Δ ρiData gridding is carried out to obtain and researchAnd (3) carrying out filtration on the shell mantle density difference grid data of the area grid gravity anomaly data grids with the same size to obtain the filtered shell mantle density difference data.
Further, the specific step of S3 is:
s31, acquiring sedimentary deposit thickness data of a research area according to shallow earthquake, well drilling or comprehensive geophysical inversion results, and acquiring sea area sedimentary deposit thickness data and land area sedimentary deposit thickness data from sedimentary deposit thickness data in a CRUST1.0 data file through the NOAA of the national ocean and atmosphere administration; the thickness data resolution of the sea area sedimentary deposit is 5 'multiplied by 5', and the sea area sedimentary deposit contains global sea area data;
s32, performing projection transformation on the settled layer thickness data, projecting the geographic coordinates into plane rectangular coordinates, and then gridding the settled layer thickness data, wherein the grid size is consistent with the grid size of the grid-laying gravity anomaly data;
s33, calculating gravity anomaly caused by thickness change of the deposition layer, and obtaining bottom interface depth data of the deposition layer by using terrain data and deposition layer thickness data during calculation; the topographic data is subjected to regional topographic elevation data ETOPO1 obtained through NOAA (national oceanic and atmospheric administration), the resolution of the topographic data is 1 'x 1', the topographic data is converted into plane rectangular coordinates, the data is gridded, and the size of the grid is consistent with that of a grid of the grid gravity anomaly data; and then calculating the gravity anomaly caused by the thickness change of the deposition layer by using a frequency domain density interface gravity anomaly forward formula:
wherein G is a universal gravitation constant;
representing a two-dimensional wavenumber field;
when the gravity abnormality caused by the thickness change of the sedimentary deposit is calculated, the density difference between the sedimentary deposit and the crust
Establishing according to CRUST1.0 data; then, calculating the frequency spectrum of a gravity field caused by the thickness of the deposit layer by utilizing the frequency domain density interface gravity anomaly forward modeling, and obtaining a gravity anomaly value through Fourier inversion;
s34, subtracting gravity anomaly caused by thickness change of the deposition layer from the Bragg gravity anomaly to obtain gravity anomaly after influence of the deposition layer is eliminated, filtering the gravity anomaly to eliminate gravity anomaly of local density inhomogeneity in the crust, and then carrying out background adjustment on the filtered gravity anomaly by utilizing gravity anomaly corresponding to a known depth point of the Mohuo surface obtained by deep reflection seismic profile data to obtain gravity anomaly caused by fluctuation change of the Mohuo surface.
Further, the specific step of S4 is:
s41, determining the average depth of the mohowl surface of the research area according to the data of the deep reflection seismic section, the seabed seismic section OBS and the CRUST 1.0;
s42, taking gravity anomaly caused by fluctuation change of the Mohol surface, filtered shell density difference and average Mohol surface depth as basic data, and obtaining the depth of the Mohol surface by utilizing a frequency domain density interface inversion method Parker-Oldenburg;
s43, filtering the terrain data of the research area, eliminating local terrain fluctuation change to objectively reflect the change characteristics of the thickness of the crust, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research area.
In addition, as a further improved technical solution of the present invention, the specific step of S4 may further be:
s41, determining the average depth of the mohowl surface of the research area according to the data of the deep reflection seismic section, the seabed seismic section OBS and the CRUST 1.0;
s42, after a shell mantle density difference of transverse change is established, inverting the depth of a Mohollo surface by using a density interface inversion method;
s421, calculating the average shell mantle density difference delta rho according to the shell mantle density difference of the whole research area0Calculating the deviation rho of the shell mantle density differencei＝Δρi-Δρ0；
S422, according to the Moire surface depth data and the average Moire surface depth of the CRUST1.0 database and the deviation of the shell mantle density difference, positively calculating the gravity anomaly caused by the uneven shell mantle density by using a formula (1), and eliminating the gravity anomaly caused by the fluctuation change of the Moire surface;
s423 to eliminate gravity anomaly and average shell screen density difference delta rho after shell screen density is not uniform0Obtaining the depth of a Mohol surface by using a frequency domain density interface inversion method Parker-Oldenburg as basic data;
s43, filtering the terrain data of the research area, eliminating local terrain fluctuation change to objectively reflect the change characteristics of the thickness of the crust, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research area.
Compared with the prior art, the method establishes the shell mantle density difference which changes transversely, is used for inversion of the depth of the Mohuo surface and the thickness of the earth crust, can accurately present the features of the fluctuation of the Mohuo surface and the thickness of the earth crust, and adopts the constant shell mantle density difference which is not consistent with the actual geological features and can not objectively reflect the thickness change of the earth crust in the whole area only along the two-dimensional shell mantle density difference of the section when the research area range is larger or the area with more complex construction.
Drawings
FIG. 1 is a flow chart of the variable density-based crust thickness gravity inversion method provided by the invention.
Fig. 2 is a topographical map of the ridges and neighborhoods of the southeast ocean of the implementation example indian ocean.
Fig. 3 is a graph of the bragg gravity anomaly of the spine and neighborhood of the eastern south-east ocean of the example of implementation.
FIG. 4 is a contour plot of shell density differences for the ridges and adjacent regions of the southeast ocean of the example of implementation.
Fig. 5 is a contour map of the thickness of the sedimentary layers of the ridges and neighborhoods in the southeast ocean of the implementation example indian ocean.
Fig. 6 is a contour plot of density differences of sedimentary layers and crust in the ridges and neighborhoods of the eastern south-east ocean of india, an example implementation.
Fig. 7 is a contour plot of gravity resulting from the thickness of the sedimentary layers in the ridges and neighborhoods of the example of implementation of the indian ocean southeast.
FIG. 8 is a Mohuo surface gravity anomaly contour plot of ridges and neighborhoods in the southeast ocean of the implementation example.
FIG. 9 is a Mohuo surface depth contour plot resulting from implementing an example Indian southeast ocean mid-ridge and neighborhood inversion.
FIG. 10 is a contour plot of crust thickness from an example implementation of the inversion of the ridges and neighborhoods in the southeast ocean, India.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. The specific embodiments described herein are merely illustrative of the invention and do not limit the invention.
As shown in fig. 1, the method for inversion of crust thickness gravity based on variable density provided by the invention comprises the following steps:
s1, acquiring crustal model CRUST1.0 crustal thickness and density data and upper mantle density data with the spatial resolution of 1 degree multiplied by 1 degree covering the whole world:
s11, opening a global CRUST model CRUST1.0 database (website: igppweb. ucsd. edu/. about. gai/CRUST 1.html), downloading a CRUST1.0 data file CRUST1.0.tar.gz in Download selection, wherein the file CRUST1.0.tar.gz comprises codes written in FORTRAN language for extracting CRUST thickness and density data and upper mantle density data;
s12, running and downloading a getCN1xyz.f file in the file by using a Microsoft Visual Studio platform, reading a thickness file of 8 density layers (a water layer, an ice layer, an upper deposition layer, a middle deposition layer, a lower deposition layer, an upper crust, a middle crust and a lower crust) covering the world and density data files of the 8 density layers and an upper mantle from the data file, wherein the thickness file covers the 8 density layers in latitude and longitude ranges of 179.5 degrees W-179.5 degrees E and 89.5 degrees S-89.5 degrees N;
and S13, extracting thickness and density data of each density layer in the research area range from the data file of S12 according to the specific latitude and longitude range of the research area, and arranging the thickness and density data of each density layer into two files of the thickness and density of each density layer.
S2, establishing shell density difference along with horizontal position change according to the density and thickness of the upper, middle and lower earth CRUST of the earth CRUST model CRUST1.0 and the density data of the upper earth CRUST:
s21, carrying out projection transformation on the thickness and density data of each density layer arranged in the S13, and projecting the geographic coordinates into plane rectangular coordinates;
s22, calculating the shell mantle density difference of each point at different horizontal positions in the research area by using the following formula:
where i denotes different horizontal position points, Δ ρiIs the difference in shell density at the ith point,
the density values of the upper crust, the middle crust, the lower crust and the upper mantle at the ith point are respectively,
the thicknesses of the upper crust, the middle crust and the lower crust at the ith point are respectively;
s23 for Δ ρiCarrying out data gridding to obtain shell mantle density difference grid data with the same size as a grid of cloth grid gravity abnormal data in a research area, and then filtering the shell mantle density difference data to obtain a filtered shellMantle density difference data.
S3, acquiring gravity anomaly caused by Mohuo surface fluctuation change due to the Buge gravity anomaly:
s31, acquiring sedimentary stratum thickness data of a research area according to shallow seismic, well drilling or comprehensive geophysical inversion results, and if the data are insufficient, acquiring sea sedimentary stratum thickness data and acquiring land sedimentary stratum thickness data according to sedimentary stratum thickness data in a CRUST1.0 data file by the NOAA (website address: www.ngdc.noaa.gov/mgg/sedthick /); the thickness data resolution of the sea area sedimentary deposit is 5 'multiplied by 5', and the sea area sedimentary deposit contains global sea area data;
s32, performing projection transformation on the settled layer thickness data, projecting the geographic coordinates into plane rectangular coordinates, and then gridding the settled layer thickness data, wherein the grid size is consistent with the grid size of the grid-laying gravity anomaly data;
s33, calculating gravity anomaly caused by thickness change of the deposition layer, and obtaining bottom interface depth data of the deposition layer by using terrain data and deposition layer thickness data during calculation; the topographic data obtains topographic elevation data ETOPO1 (website: https:// www.ngdc.noaa.gov/mgg/global. html) of a research area through NOAA (national oceanic and atmospheric administration), the resolution of the topographic data is 1 '× 1', the topographic data is converted into plane rectangular coordinates, gridding is carried out on the data, and the size of the grid is consistent with that of a grid of grid-laying gravity anomaly data; and then calculating the gravity anomaly caused by the thickness change of the deposition layer by using a frequency domain density interface gravity anomaly forward formula:
wherein G is a universal gravitation constant;
representing a two-dimensional wavenumber field;
when the gravity abnormality caused by the thickness change of the sedimentary deposit is calculated, the density difference between the sedimentary deposit and the crust
Establishing according to CRUST1.0 data; then, calculating the frequency spectrum of a gravity field caused by the thickness of the deposit layer by utilizing the frequency domain density interface gravity anomaly forward modeling, and obtaining a gravity anomaly value through Fourier inversion;
s34, subtracting gravity anomaly caused by thickness change of the deposition layer from the Bragg gravity anomaly to obtain gravity anomaly after influence of the deposition layer is eliminated, filtering the gravity anomaly to eliminate gravity anomaly of local density inhomogeneity in the crust, and then carrying out background adjustment on the filtered gravity anomaly by utilizing gravity anomaly corresponding to a known depth point of the Mohuo surface obtained by deep reflection seismic profile data to obtain gravity anomaly caused by fluctuation change of the Mohuo surface.
S4, inverting the depth of the Mohol surface and obtaining the thickness of the crust by the gravity anomaly caused by the fluctuation change of the Mohol surface:
s41, determining the average depth of the mohowl surface of the research area according to the data of the deep reflection seismic section, the seabed seismic section OBS and the CRUST 1.0;
s42, taking gravity anomaly caused by fluctuation change of the Mohol surface, filtered shell density difference and average Mohol surface depth as basic data, and obtaining the depth of the Mohol surface by utilizing a frequency domain density interface inversion method Parker-Oldenburg;
s43, filtering the terrain data of the research area, eliminating local terrain fluctuation change to objectively reflect the change characteristics of the thickness of the crust, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research area.
S42 in this step may also be replaced with the following scheme, specifically, after setting up the difference in shell density in lateral variation, inverting the mojolt face depth using a density interface inversion method:
s421, calculating the average shell density difference according to the shell density difference of the whole research areaΔρ0Calculating the deviation rho of the shell mantle density differencei＝Δρi-Δρ0；
S422, according to the Moire surface depth data and the average Moire surface depth of the CRUST1.0 database and the deviation of the shell mantle density difference, the gravity anomaly caused by the uneven shell mantle density is calculated by using the following forward modeling, and the gravity anomaly caused by the fluctuation change of the Moire surface is eliminated:
wherein G is a universal gravitation constant;
representing a two-dimensional wavenumber field;
when the gravity abnormality caused by the thickness change of the sedimentary deposit is calculated, the density difference between the sedimentary deposit and the crust
Establishing according to CRUST1.0 data; then, calculating the frequency spectrum of a gravity field caused by the thickness of the deposit layer by utilizing the frequency domain density interface gravity anomaly forward modeling, and obtaining a gravity anomaly value through Fourier inversion;
s423 to eliminate gravity anomaly and average shell screen density difference delta rho after shell screen density is not uniform0Obtaining the depth of a Mohol surface by using a frequency domain density interface inversion method Parker-Oldenburg as basic data;
s43, filtering the terrain data of the research area, eliminating local terrain fluctuation change to objectively reflect the change characteristics of the thickness of the crust, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research area.
In order to make those skilled in the art better understand the technical solution of the present application, the technical solution of the present invention is explained below by taking inversion of the depths of the mohuo surface and the thickness of the crust in the central and adjacent regions of the southeast ocean, indian as an example (the research range is 75 ° E-150 ° E, 30 ° S-65 ° S).
Firstly, downloading a CRUST1.0 data file CRUST1.0.tar.gz through a global CRUST model CRUST1.0 database, and running a getCN1xyz.f file in the downloaded file by using a Microsoft Visual Studio platform so as to read the thickness of each density layer (water layer, ice layer, upper deposition layer, middle deposition layer, lower deposition layer, upper CRUST, middle CRUST and lower CRUST) covering the whole world and the density data of the 8 density layers and the upper mantle. And respectively extracting the thickness and density data of each density layer in the research range of the hinde and the adjacent region of the hinde in India south-east China (75-150 degrees E, 30-65 degrees S), sorting, and then performing projection transformation to obtain the thickness and density data of each density layer in the hinde and the adjacent region of the hinde in India south-east China in a plane rectangular coordinate system.
The terrain data and the Bragg gravity anomaly data of the ridges and the adjacent regions in the southeast ocean of India are downloaded and come from the global gravity model WGM2012 (the website: http:// bgi. omp. obs-mip.fr/data-products/Grids-and-models/WGM2012) in the example. Fig. 2 is a topographical view of the ocean center and neighborhood of the southeast ocean of the example indian ocean, where the vast majority of the water depth in the ocean exceeds 3500m and the ocean floor topography at the ocean center is a pronounced uplift caused by the expansion of the ocean hulls and the upwelling of the vallecular material. Fig. 3 is a graph of the buerger gravity anomaly in the ridges and adjacent regions of the southeast ocean, indian, where the overall trend of the buerger gravity anomaly reflects the fluctuation of the mohol surface, and thus the negative value of the gravity anomaly is relatively low in the land area (e.g., australia in the northwest corner of the figure), and the negative value of the gravity anomaly is relatively high in the ocean area (e.g., large area in the middle of the figure), which reflects that the depth of the mohol surface is relatively large in the land area, and the depth of the mohol surface is relatively shallow in the ocean area.
And then calculating the shell mantle density difference of the horizontal position points of the ridge and the adjacent areas in the southeast ocean of India according to a weighted average density calculation formula, wherein the formula is as follows:
gridding the calculation result to obtain shell mantle density difference grid data with the same size as the grid of the cloth grid gravity anomaly data of the research area, filtering the shell mantle density difference data to eliminate local density nonuniformity for eliminating the influence of local density change, and finally obtaining the shell mantle density difference of the south eastern ocean of India and the adjacent areas as shown in figure 4. The density difference of the shell mantle is larger in the land area and is-0.5 g/cm3On the left and right, the difference in the density of the mantle in the ocean area is small, mostly less than-0.4 g/cm3This is due to the relatively complete structure of the crust in the land area, while the upper crust is generally missing in the ocean area.
Since the middle and adjacent areas of the southeast ocean of India lack information of shallow earthquake, well drilling and the like, and the sediment thickness data of the research area can not be obtained, the sea area sediment thickness data is obtained through the National Ocean and Atmosphere Administration (NOAA) (website: www.ngdc.noaa.gov/mgg/settick /), and the sediment thickness data of the research area is supplemented by the sediment thickness data in the CRUST1.0 data file. The thickness of the sedimentary layer of the ridges and the adjacent areas in the southeast ocean of India is obtained by data stitching and gridding as shown in FIG. 5. The area with the largest sedimentary deposit thickness is located in the offshore area of Australia, the maximum thickness exceeds 7km, the sedimentary deposit thickness in other areas is smaller, and the ocean area is generally less than 100 m.
The density difference between the deposition layer and the substrate obtained by the method for calculating the same shell mantle density difference is shown in fig. 6, the depth data of the bottom interface of the deposition layer is obtained by using the data of the terrain and the thickness of the deposition layer, the terrain and the bottom interface of the deposition layer are respectively used as an upper interface and a lower interface, and the gravity anomaly caused by the thickness change of the deposition layer is calculated by using a frequency domain density interface gravity anomaly forward formula, wherein the formula is as follows:
wherein G is a universal gravitation constant;
representing a two-dimensional wavenumber field;
the frequency spectrum of the gravity field caused by the thickness of the deposition layer is calculated by using the formula, and then the gravity anomaly caused by the thickness change of the deposition layer is obtained through Fourier inversion, as shown in figure 7. Since the sediment layer in the ocean area is very thin, the gravity anomaly caused by the sediment layer is less than-10 mGal in a large part of the area, and the gravity anomaly caused by the sediment layer at the position with the largest thickness can reach-160 mGal at the maximum.
Subtracting the gravity anomaly caused by the thickness change of the deposition layer from the Bragg gravity anomaly to obtain the gravity anomaly after eliminating the influence of the deposition layer, filtering the anomaly to eliminate the gravity anomaly of the local density inhomogeneity in the CRUST, and then performing background adjustment on the filtered gravity anomaly by using the gravity anomaly corresponding to the known depth point of a part of the Moros surface in the CRUST1.0 to obtain the gravity anomaly caused by the fluctuation change of the Moros surface as shown in FIG. 8. The gravity anomaly is basically consistent with the overall region of the Bragg gravity anomaly, is different only in a local region and is mainly expressed in a region near the ocean ridge.
The average depth of the mohol surface of the research area is determined according to CRUST1.0 data and the like, and the depth of the mohol surface is obtained by taking the gravity anomaly, the shell mantle density difference and the average depth of the mohol surface caused by fluctuation change of the mohol surface as basic data and utilizing a frequency domain density interface inversion method (namely a Parker-Oldenburg method) for inversion, and is shown in figure 9. According to an inversion result, the depths of the Mohuo surfaces of the ridges and the adjacent regions in the southeast ocean of India ocean of the implementation example are 10-38 km, the region with the largest depth is positioned in Australia, the depth of the Mohuo surfaces generally exceeds 30km, the depths of the Mohuo surfaces of most regions in the India ocean are close to 10km, the depths of the Mohuo surfaces of individual small regions are relatively large and close to 20km, and the depth can be reflected by micro-plots in the big ocean.
Filtering the terrain data of the ridges and the adjacent regions in the southeast ocean of India and the south of east China to eliminate local terrain fluctuation, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research region as shown in figure 10. According to the calculation result, the thickness of the crust in the ridge and the neighborhood in the southeast ocean of India ocean is 6-38 km, the thickness of the crust in Australia generally exceeds 32km, the thickness of the crust in most regions in the Atlantic ocean is 7km, and the thickness of the crust of individual small blocks exceeds 18km, which is probably reflected by micro-blocks in the ocean. In addition, the thickness of the crust at the ocean ridge is 8-9 km, which is slightly larger than the thickness of the crust at the side of the ocean ridge, and the thickness is presumed to be caused by slow expansion of the ocean ridge and accumulation of the substances of the mantle.
The embodiment further verifies that the density difference of the upper part and the lower part of the interface (namely the Mohol surface) of the crust and the mantle, which changes along with the position, is established according to the density data of different layer positions of the crust, and the depth of the Mohol surface is inverted by utilizing the gravity density interface inversion technology, so that the thickness of the crust is obtained. Therefore, the technology provided by the invention has the characteristics of simplicity, easiness in operation and higher accuracy. Compared with the existing inversion method adopting the constant density, the inversion method has higher accuracy of the inversion result, better accords with the geological characteristics of the crust, and can provide key technical support for deep structure research and regional research.
The technical solution of the present invention is not limited to the limitations of the above specific embodiments, and all technical modifications made according to the technical solution of the present invention fall within the protection scope of the present invention.
Claims (5)
1. The crust thickness gravity inversion method based on variable density is characterized by comprising the following steps:
s1, acquiring crustal model CRUST1.0 crustal thickness and density data and upper mantle density data with the global coverage spatial resolution of 1 degree multiplied by 1 degree;
s2, establishing a shell mantle density difference which changes along with the horizontal position according to the density and thickness of the upper, middle and lower earth shells and the density data of the upper mantle of the earth shell model CRUST 1.0;
s3, acquiring gravity anomaly caused by fluctuation change of the Mohuo surface due to the Buge gravity anomaly;
s4, inverting the depth of the Mohuo surface and obtaining the thickness of the earth crust by the gravity anomaly caused by the fluctuation change of the Mohuo surface,
the specific steps of S1 are as follows:
s11, opening a global CRUST model CRUST1.0 database with the spatial resolution of 1 degree multiplied by 1 degree, downloading a CRUST1.0 data file CRUST1.0.tar.gz in a Download selection, wherein the file CRUST1.0.tar.gz comprises codes written in FORTRAN language and is used for extracting CRUST thickness and density data and upper mantle density data;
s12, running and downloading a getCN1xyz.f file in the file by using a Microsoft Visual Studio platform, and reading thickness files of 8 density layers covering the world in latitude and longitude ranges of 179.5-degree E, 89.5-degree S-89.5-degree N and density data files of the 8 density layers and the upper mantle from the data files;
s13, extracting thickness and density data of each density layer in the research area range from the data file of S12 according to the specific latitude and longitude range of the research area, and arranging the thickness and density data of each density layer into two files,
the specific steps of S2 are as follows:
s21, carrying out projection transformation on the thickness and density data of each density layer arranged in the S13, and projecting the geographic coordinates into plane rectangular coordinates;
s22, calculating the shell mantle density difference of each point at different horizontal positions in the research area by using the following formula:
where i denotes different horizontal position points, Δ ρiIs the difference in shell density at the ith point,
the density values of the upper crust, the middle crust, the lower crust and the upper mantle at the ith point are respectively,
the thicknesses of the upper crust, the middle crust and the lower crust at the ith point are respectively;
s23 for Δ ρiAnd carrying out data gridding to obtain shell mantle density difference grid data with the same size as the cloth grid gravity anomaly data grid in the research area, and then filtering the shell mantle density difference data to obtain the filtered shell mantle density difference data.
2. The gravity inversion method for crust thickness based on variable density according to claim 1, wherein the specific steps of S3 are as follows:
s31, acquiring sedimentary deposit thickness data of a research area according to shallow earthquake, well drilling or comprehensive geophysical inversion results, and acquiring sea area sedimentary deposit thickness data and land area sedimentary deposit thickness data from sedimentary deposit thickness data in a CRUST1.0 data file through the NOAA of the national ocean and atmosphere administration; the thickness data resolution of the sea area sedimentary deposit is 5 'multiplied by 5', and the sea area sedimentary deposit contains global sea area data;
s32, performing projection transformation on the settled layer thickness data, projecting the geographic coordinates into plane rectangular coordinates, and then gridding the settled layer thickness data, wherein the grid size is consistent with the grid size of the grid-laying gravity anomaly data;
s33, calculating gravity anomaly caused by thickness change of the deposition layer, and obtaining bottom interface depth data of the deposition layer by using terrain data and deposition layer thickness data during calculation; the topographic data is subjected to regional topographic elevation data ETOPO1 obtained through NOAA (national oceanic and atmospheric administration), the resolution of the topographic data is 1 'x 1', the topographic data is converted into plane rectangular coordinates, the data is gridded, and the size of the grid is consistent with that of a grid of the grid gravity anomaly data; and then calculating the gravity anomaly caused by the thickness change of the deposition layer by using a frequency domain density interface gravity anomaly forward formula:
wherein G is a universal gravitation constant;
representing a two-dimensional wavenumber field;
when the gravity abnormality caused by the thickness change of the sedimentary deposit is calculated, the density difference between the sedimentary deposit and the crust
Establishing according to CRUST1.0 data; then, calculating the frequency spectrum of a gravity field caused by the thickness of the deposit layer by utilizing the frequency domain density interface gravity anomaly forward modeling, and obtaining a gravity anomaly value through Fourier inversion;
s34, subtracting gravity anomaly caused by thickness change of the deposition layer from the Bragg gravity anomaly to obtain gravity anomaly after influence of the deposition layer is eliminated, filtering the gravity anomaly to eliminate gravity anomaly of local density inhomogeneity in the crust, and then carrying out background adjustment on the filtered gravity anomaly by utilizing gravity anomaly corresponding to a known depth point of the Mohuo surface obtained by deep reflection seismic profile data to obtain gravity anomaly caused by fluctuation change of the Mohuo surface.
3. The gravity inversion method for crust thickness based on variable density according to claim 1, wherein the specific steps of S4 are as follows:
s41, determining the average depth of the mohowl surface of the research area according to the data of the deep reflection seismic section, the seabed seismic section OBS or CRUST 1.0;
s42, taking gravity anomaly caused by fluctuation change of the Mohol surface, filtered shell density difference and average Mohol surface depth as basic data, and obtaining the depth of the Mohol surface by utilizing a frequency domain density interface inversion method Parker-Oldenburg;
s43, filtering the terrain data of the research area, eliminating local terrain fluctuation change to objectively reflect the change characteristics of the thickness of the crust, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research area.
4. The gravity inversion method for crust thickness based on variable density according to claim 1, wherein the specific steps of S4 are as follows:
s41, determining the average depth of the mohowl surface of the research area according to the data of the deep reflection seismic section, the seabed seismic section OBS or CRUST 1.0;
s42, after a shell mantle density difference of transverse change is established, inverting the depth of a Mohollo surface by using a density interface inversion method;
s421, calculating the average shell mantle density difference delta rho according to the shell mantle density difference of the whole research area0Calculating the deviation rho of the shell mantle density differencei＝Δρi-Δρ0；
S422, according to the Moholt surface depth data and the average Moholt surface depth of the CRUST1.0 database and the deviation of shell mantle density difference, calculating the gravity anomaly caused by uneven shell mantle density by using a frequency domain density interface gravity anomaly forward formula, and eliminating the gravity anomaly caused by fluctuation change of the Moholt surface;
s423 to eliminate gravity anomaly and average shell screen density difference delta rho after shell screen density is not uniform0Obtaining the depth of a Mohol surface by using a frequency domain density interface inversion method Parker-Oldenburg as basic data;
s43, filtering the terrain data of the research area, eliminating local terrain fluctuation change to objectively reflect the change characteristics of the thickness of the crust, and then subtracting smooth terrain fluctuation from the depth of the Mohuo surface to obtain the thickness of the crust of the research area.
5. The variable density-based crust thickness gravity inversion method according to claim 1, wherein: the 8 density layers in the S12 are water layer, ice layer, upper deposition layer, middle deposition layer, lower deposition layer, upper crust, middle crust and lower crust.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910502882.4A CN110244352B (en) | 2019-06-11 | 2019-06-11 | Variable density-based crust thickness gravity inversion method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910502882.4A CN110244352B (en) | 2019-06-11 | 2019-06-11 | Variable density-based crust thickness gravity inversion method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110244352A CN110244352A (en) | 2019-09-17 |
| CN110244352B true CN110244352B (en) | 2020-12-25 |
Family
ID=67886681
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910502882.4A Expired - Fee Related CN110244352B (en) | 2019-06-11 | 2019-06-11 | Variable density-based crust thickness gravity inversion method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110244352B (en) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112666607B (en) * | 2019-10-16 | 2024-06-25 | 中国石油天然气集团有限公司 | Method and device for inverting thickness distribution of yellow soil layer by gravity |
| CN111045098B (en) * | 2019-12-27 | 2022-06-24 | 核工业北京地质研究院 | Method for picking up underground deep structure information |
| CN111337993A (en) * | 2020-03-30 | 2020-06-26 | 中国地质科学院 | Variable density and variable depth constraint-based gravity density interface inversion method |
| CN111999778A (en) * | 2020-07-13 | 2020-11-27 | 国家海洋信息中心 | Antarctic continental Mohuo surface depth inversion method based on satellite gravity gradient data |
| CN113552644B (en) * | 2021-07-05 | 2022-11-11 | 中国科学院地质与地球物理研究所 | Density determination method and device and electronic equipment |
| CN113985490B (en) * | 2021-09-22 | 2023-05-05 | 中国人民解放军战略支援部队信息工程大学 | Method and device for carrying out surface gravity simulation by utilizing terrain and crust density data |
| CN115373024B (en) * | 2022-08-09 | 2023-04-18 | 中国科学院南海海洋研究所 | Method and device for inverting passive land edge crustal structure based on stratum recording settlement |
| CN117741795A (en) * | 2023-12-21 | 2024-03-22 | 中国自然资源航空物探遥感中心 | Method for jointly inverting rock ring density by using multi-source data |
| CN117555025B (en) * | 2024-01-11 | 2024-04-02 | 应急管理部国家自然灾害防治研究院 | Multi-layer crust structure inversion method based on gravity data |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2497576B (en) * | 2011-12-15 | 2017-06-14 | Statoil Petroleum As | ASEP+D Method:identifying anomalous areas of the earth's lower crust |
| CN104459795A (en) * | 2014-12-08 | 2015-03-25 | 中国科学院南海海洋研究所 | Depth-varying-to-density earth crust extension coefficient thermal calibration gravity anomaly retrieval method |
| CN104835201B (en) * | 2015-05-08 | 2017-07-04 | 华东师范大学 | A kind of method simulated on digital earth software platform and show Global Crustal Structure |
| US10474767B2 (en) * | 2016-01-26 | 2019-11-12 | Saudi Arabian Oil Company | Gravity modeling a rifted continental margin |
| CN106886047A (en) * | 2017-02-28 | 2017-06-23 | 中国地质大学（北京） | A kind of method of receiver function and gravity Inversion CRUSTAL THICKNESS and ripple ratio |
| CN108919338A (en) * | 2018-05-28 | 2018-11-30 | 中国地震局地震预测研究所 | Method based on earth's surface gravity and GNSS observation data prediction reservoir-induced earthquake |
-
2019
- 2019-06-11 CN CN201910502882.4A patent/CN110244352B/en not_active Expired - Fee Related
Non-Patent Citations (3)
| Title |
|---|
| "MohoIso: A MATLAB program to determine crustal thickness by an isostatic and a global gravitational model";Mohammad Bagherbandi;《Computers & Geosciences》;20121231(第44期);第177-183页 * |
| "The development and evaluation of the Earth Gravitational Model 2008（EGM2008）";Pavlis N K等;《J. Geophys. Res》;20121231;第117卷;第B04406-1-38页 * |
| "Update on CRUST1.0-A1-degree Global Model of Earth’s Crust;Laske G等;《Geophys. Res. Abstracts》;20131231;第15卷;第2013-2658页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110244352A (en) | 2019-09-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110244352B (en) | Variable density-based crust thickness gravity inversion method | |
| Rovere et al. | The analysis of Last Interglacial (MIS 5e) relative sea-level indicators: Reconstructing sea-level in a warmer world | |
| Zhu et al. | Land subsidence due to groundwater withdrawal in the northern Beijing plain, China | |
| DiBiase et al. | Deltaic deposits at Aeolis Dorsa: Sedimentary evidence for a standing body of water on the northern plains of Mars | |
| Adams et al. | Digital characterization of thrombolite-stromatolite reef distribution in a carbonate ramp system (terminal Proterozoic, Nama Group, Namibia) | |
| CN105425315B (en) | A kind of small scale depositional phase ancient landform inversion method in weak structural deformation area | |
| Han et al. | Pervasive low-velocity layer atop the 410-km discontinuity beneath the northwest Pacific subduction zone: Implications for rheology and geodynamics | |
| Bosman et al. | The first ultra-high resolution Digital Terrain Model of the shallow-water sector around Lipari Island (Aeolian Islands, Italy) | |
| CN105868482B (en) | A kind of deposition phase spends the projectional technique and device of palaeohigh a little | |
| Lee et al. | Stratigraphy of late Quaternary deposits using high resolution seismic profile in the southeastern Yellow Sea | |
| Cella et al. | High-resolution geophysical 3D imaging for archaeology by magnetic and EM data: the case of the iron age settlement of Torre Galli, Southern Italy | |
| CN104268848A (en) | Ocean internal wave velocity monitoring method | |
| Dávila et al. | The evolution of the high‐elevated depocenters of the northern Sierras Pampeanas (ca. 28˚ SL), Argentine broken foreland, South‐Central Andes: the Pipanaco Basin | |
| Khakim et al. | Detection of localized surface uplift by differential SAR interferometry at the hangingstone oil sand field, Alberta, Canada | |
| Cuffaro et al. | The Ventotene Volcanic Ridge: a newly explored complex in the central Tyrrhenian Sea (Italy) | |
| Callens et al. | Temporally stable surface mass balance asymmetry across an ice rise derived from radar internal reflection horizons through inverse modeling | |
| Pappa et al. | Gravity, magnetics and geothermal heat flow of the Antarctic lithospheric crust and mantle | |
| Lodolo et al. | Gravity map of the Isla Grande de Tierra del Fuego, and morphology of Lago Fagnano | |
| Dahm et al. | Combining geophysical data sets to study the dynamics of shallow evaporites in urban environments: application to Hamburg, Germany | |
| Liu et al. | Holocene coastal morphologies and shoreline reconstruction for the southwestern coast of the Bohai Sea, China | |
| Ishitsuka et al. | Ground uplift related to permeability enhancement following the 2011 Tohoku earthquake in the Kanto Plain, Japan | |
| Urlaub et al. | Gravity crustal models and heat flow measurements for the Eurasia Basin, Arctic Ocean | |
| Akawwi | Geomorphology using geographic information system and globel mapper | |
| Nunes et al. | Reserve evaluation of a fault-conditioned aquifer: the Barreiras Aquifer in the coastal region of NE Brazil | |
| Michels et al. | Geophysical expression of the Leka Ophiolite, Norway modeled from integrated gravity, magnetic and petrophysical data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201225 Termination date: 20210611 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |