Application of multi-attribute matching technology based on geological models for sedimentary facies:A case study of the 3rd member in the Lower Jurassic Badaowan Formation,Hongshanzui area,Junggar Basin,China

Zong-Quan Yao ,Fan Yang ,Dlqiati Jianatayi ,Wi Wang ,Yang Gao ,L Sun ,Wn-Fng Wang ,Bo-Yu Mu ,Pin-Bai Pan

a School of Geology and Mining Engineering,Xinjiang University,Urumqi,Xinjiang 830047,China

b Key Laboratory of Central Asian Orogenic Belts and Continental Dynamics,Xinjiang University,Urumqi,Xinjiang 830047,China

c Research Institute of Peleoleum Exploration and Development,PetroChina,Beijing,100083,China

d National 305 Project Office,Urumqi,Xinjiang 830000,China

e Beijing Sunshine Geo-Tech Co.Ltd.,Beijing,100080,China

fCNOOC Research Institute Co.,Ltd.,Beijing,100028,China

Keywords:Clustering analysis Mapping the sedimentary facies Normal distribution Multi-attribute matching Hongshanzui area

ABSTRACT Regarding high drilling costs,an effort should be made to substantially reduce the drilling operation.To achieve this goal,exploration and development stages should be carried out precisely with maximum information acquired from the reservoir.The use of multi-attribute matching technology to predict sedimentary system has always been a very important but challenging task.To resolve the challenges,we utilized a quantitative analysis method of seismic attributes based on geological models involving high resolution 3D seismic data for sedimentary facies.We developed a workflow that includes core data,seismic attribute analysis,and well logging to highlight the benefit of understanding the facies distribution in the 3rd Member of the Lower Jurassic Badaowan Formation,Hongshanzui area,Junggar Basin,China.1) Data preprocessing.2) Cluster analysis.3) RMS attribute based on a normal distribution constrains facies boundary.4) Mapping the sedimentary facies by using MRA (multiple regression analysis)prediction model combined with the lithofacies assemblages and logging facies assemblages.The confident level presented in this research is 0.745,which suggests that the methods and data-mining techniques are practical and efficient,and also be used to map facies in other similar geological settings.?2021 The Authors.Publishing services by Elsevier B.V.on behalf of KeAi Communications Co.Ltd.This

1.Introduction

With the increasing difficulty of oil and gas exploration,the identification of effective sand bodies becomes a key.Information obtained from cores samples and well logs are localized and represents properties from a small portion of the reservoir.And such samples are generally very limited in availability due to the expense of drilling those wells.Seismic attribute technology can extract information from seismic database that is otherwise hidden in the data.Additionally,it can enhance the use and value of geophysics for predicting,characterizing,and interpreting lithological features,which is now a prime focus of the petroleum industry (Sidney,1997).Rummerfield (1954) applied seismic attributes to oil and gas exploration for the first time and accurately predicted the relevant faults in the fracture region.During the 1970s and 1980s,the seismic attributes most used in petroleum exploration were amplitude-based instantaneous attributes (Justice et al.,1985;Balch,2012).However,in the 1990s,seismic attribute technology dramatically advanced in every aspect of hydrocarbon exploration and development (Chopra and Marfurt,2005).Its application ranged from computations of single-trace instantaneous event attributes to more complex multitrace windowed seismic event attribute extractions to the generation of seismic attribute volumes(Alao et al.,2014;Bueno et al.,2014;Herrera Volcan et al.,2015;Torabi et al.,2016).

Although the application of conventional seismic have achieved some results.It is difficult to establish the corresponding relationships between attributes and reservoirs,as the application of seismic attributes is often uncertain with multiple interpretations.Moreover,the workload of multiple attributes is affected by subjective factors (Xu,2010).Therefore,conventional seismic interpretation can no longer meet the current needs of oil and gas exploration.To this end,scholars put forward seismic attribute matching technology.It is no longer based on experience or statistical methods,whereas the quantitative analysis technology of multi-attribute matching is used to determine the plane distribution of sand body,which is an essential step for the construction of reliable geologic models.

Multi-attribute matching technology is used to fuse single standardized spatial attributes in intensity,hue,saturation (IHS)coordinate system by a mathematical algorithm.Compared with the 256 colors of the original single-attribute display,the pixels of the fusion display are approximately 16 million dpi.At present,multi-attribute matching technology can be divided into RGB attribute fusion,cluster analysis attribute matching,multiple linear regression attribute matching,well attribute matching,and seismic attribute matching based on fuzzy logic (Li et al.,2015).There are generally two types of methods for multi-attribute technology:empirical method and mathematical method(Zhu et al.,2003;Gan et al.,2018).Uwe Strecker (2004) applied this technology to carbonate reservoirs.Wei et al.(2013) used seismic multi-attribute matching technology to analyze sedimentary and reflection characteristics,and determined the distribution of sedimentary facies.Xie et al.(2014) proposed a seismic attribute matching method based on contourlet transformation and applied it to a meandering channel.Chen and Cheng (2014) believed that a single seismic attribute could not fully reflect the geological problem,so they proposed the multi-attribute matching method of RGB-HIS to greatly improve the accuracy of sandstone target identification.Yang et al.(2020) predicted the thickness of the N873 sandbody based on the seismic attribute matching.

Above all,multi-attribute matching technologies improve the image resolution for determining thickness and macroscopic sedimentary facies.And these methods mainly based on geophysical principles,and lack prediction guided by geological models(And and Rowell,2012).Furthermore,study area contains coarsegrained reservoirs with conglomerates,fast changing lithofacies,strong reservoir heterogeneity (Yao et al.,2017),unclear rock physics,complex relationships between seismic reflection intensity and the rock (Pan et al.,2017),and difficult to define the fan body boundaries (Xiao et al.,2016).Additionally,the spatial configuration of the relationship between attributes and sedimentary microfacies cannot be accurately predicted (Chen et al.,2009).Thus,these difficult makes promoting the application of seismic multi-attribute matching technologies and multiple regression analysis (MRA) for mapping facies.

The aim of this study is to eliminate the bias of single-attribute reflections in interpreting facies distribution,and the mutual interference of related information among similar attributes,as well as in the case when core sample data are sparse or limited.This study focuses on the 3rd Member of the Lower Jurassic Badaowan Formation,Hongshanzui area based on MRA,RMS attribute based on normal distribution.Then,the result is matched with selected attributes to improve the reliability of the predicted reservoir distribution.A novel characterization method is proposed for mapping guidance and interpretation of sedimentary facies with combined well and seismic data.

2.Geological setting

The Junggar Basin is in the Xinjiang Uygur Autonomous Region.The boundary of the basin is a large thrust fault zone,with overlying strata,pitch-out,and denudation,forming many types of traps (Tao et al.,2006).Influenced by the late Carboniferous and early Permian tectonic movements,the boundary of the Junggar Basin was restricted by the surrounding mountains in the late Paleozoic,and the strata from the middle Carboniferous to the Quaternary were deposited(Novikov,2013).As one of the large oilbearing basins in northwest China,the Junggar Basin contains many hydrocarbon fields of the Permian,Triassic,and Jurassic age.

The Hongshanzui oil field is located in the northwest Junggar Basin approximately 30 km from the southern part of Karamay City.Vegetation is scarce,and the terrain is flat,with an average altitude of 280 m,and the area is dominated by Gobi dunes (Xing,2008).During the Hercynian and Indochinese periods,affected by compression from the northwest to the southeast,the fault block produced multiple stages of thrust faults,and caused the Badaowan Formation unconformity to overlie the underlying strata (Fig.1).The study area of the Badaowan Formation is divided into four sections from bottom to top:the 5th member,4th member,3rd member and 1st member,with the 2nd member being undeveloped.Specifically,the 3rd member is the target section.This single member is approximately 8-20m thick.The 5th member and 3rd member consist of gravelly sandstone,the 4th member is sandstone,and the 1st member is composed of coal layers.

3.Data and methods

3.1.Data

The study area is approximately 60 km2.More than 100 wells have been drilled for oil and gas exploration,and including 11 wells were cored.However,the wells are mainly distributed around the study area,with few wells in the center.Qualitative petrographic descriptions are available,and the logs presented in the dataset are gamma ray (GR),density (RHOB),resistivity (ILD),and sonic (DT).

The study area is covered by a 3D seismic survey.One great advantage of 3D seismic over 2D is the ability to extract 3D multitrace seismic attributes along defined events using different spatial pattern configurations.The method adopted in this research is to invert the seismic data to rock properties,and then interpret these rock properties in terms of reservoir properties.There are twenty kinds of seismic attributes in this research,but they are not good for calculation.Let a be the INST FREQ SLOPE,b be the REF STR SLOPE,…,u be the ISO CJROON (Table 1).

Table 1 Correspongding relationship between the seismic attributes and number.

3.2.Methods

A seismic geomorphology workflow was used for attribute analysis based on geological models for mapping facies(Fig.2).The method can be divided into five steps:1)Data preprocessing.2)Rmode cluster analysis.3) Use the max-min normalization method to normalize the RMS attribute for attribute extraction.4) MRA,ifR2is greater than 0.7,stop(if not,repeat the attribute optimization step untilR2greater than 0.7).5) Match the RMS threshold limit value and the gravelly sandstone content again to obtain the matched gravelly sandstone content contour map.Such analysis is primarily used for oil prospecting in new regions where no well information is available (Saggaf et al.,2000).

Fig.1.Study area location and structure of Hongshanzui area.Study area is indicated by red rectangle.The northwestern region of Hongshanzui area is Zaire Mountain,the eastern and western regions are Zhongguai uplift,the southeastern region is Shawan sag and Chepaizi uplift,which is a complex and large fault block.

3.2.1.Data preprocessing

1) Core correction

Because core analysis and logging data are collected independent of each other,the lithology at corresponding depths is not recorded at the same depth.Before establishing a logging interpretation model,the core should be returned to its original position(Yao et al.,2017).The DEN curve and apparent density,the AC curve and porosity were matched to each other.The results show that the porosity and AC curve have a high degree of anastomosis (Fig.3).The AC curve and porosity were selected for core correction.

2) Logging crossplot

The study area developed multi-component conglomerates which contain many radioactive materials.Thus,the relationship between the GR curve and the volcanic rock is not well established.In this paper,RT and SP curves were selected to construct the RT-SP crossplot(Xu et al.,2012),and the DEN curve is used as a reference to assist with lithological discrimination(Fig.4).

The crossplot is used to discriminate the lithology of the study area.A conglomerate is indicated by RT ＞11 Ω·m and SP ＜-50 mV.When RT is 6-11 Ω·m and SP is-60～-35 mV,sandstone is present.When RT is 5-7 Ω·m and SP is -30～-15 mV,siltstone is present.When RT is 4-6 Ω·m and SP is-15～-5 mV,mud is present.Coal is indicated by RT ＞15 Ω·m and DEN＜1.9 g/cm3(Table 2).

Table 2 Lithological discrimination criteria of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.

Based on these criteria,the lithological discrimination,without well cuttings was carried out using the crossplot.The method paves the way for mapping sedimentary facies.

3) Synthetic seismograms

Synthetic seismograms were generated for wells to link logs(in the depth domain) to seismic data (in the time domain) and to observe the seismic characteristics of gravelly sandstone within the area.It is the basis of reservoir prediction and attribute analysis.

Fig.2.The workflow of multi-attribute matching technology.

Well to seismic ties were performed by establishing a correlation between the seismic and synthetic seismograms with adjusting T-D functions through stretch and squeeze (Leszczynski and Steinberger,2014).50 wells were selected to develop the synthetic seismograms in this paper.

Fig.5 shows the result of the Hongshan 3 well,the wave groups correspond well,the relationship between the strong and weak amplitudes is consistent,and the purpose of synthetic seismograms is realized.

3.2.2.Multiple regression analysis (MRA)

We can seek to derive a set of constantsa0,a1,a2,…,amfor the following linear combination ofxj(j=1,2,…,m).

where the constantsa0,a1,a2,…,amare deduced using regression criteria and calculated by the successive regression analysis of MRA(Morrison,2005).

4.Lithofacies and logging facies assemblages

Lithofacies assemblage is a regular assemblage of lithofaices types in section,which represents the vertical assemblage characteristics of sedimentary environment formed under different sedimentary environment in this area.Yao(2021)has been divided into thirteen lithofacies types based on lithological characteristic analysis,including grain size,sedimentary structure,and support forms.At the same time,a regular combination of six kinds of lithofacies assemblages are established,named FA-1 to 6,which reflect the comprehensive characteristics of the sedimentary environment.

Well logging data contain abundant subsurface rock information with good continuity,and high longitudinal resolution,and the data are easy to collect,thus,these data are commonly used in geology.In this study,the traditional method used to determine the number of logging facies by lithology description is improved,and the clustering method is adopted.Based on the lithology description of 11 core wells in the study area,the logging is calibrated,and eight types of logging facies assemblages are identified,and characterized by the amplitude,shape,and top-bottom contact relationship (Fig.6).

5.Discussion

Though,there are many methods for mapping sedimentary facies.One of the most commonly used sedimentary faices mapping methods is single factor analysis and multi-factor integrated mapping (Feng,2004).Lin et al.(2017) proposed quantitative sedimentary faices mapping using the thickness content combination of main lithologies.With the complication of geological conditions and the increase of difficulty in exploration,it is no longer possible to work out map facies based on a single geological parameter or seismic data.

Fig.3.Core correction in the study area.The DEN logging and apparent density,AC and porosity were matched.

Fig.4.Lithology crossplot of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.

Identification of multi-attributes through interpreting data from different parts of the reservoir and classifying them based on seismic response properties played a key role in this research.The seismic attributes were evaluated by comparing the attributes to each other,and subsequently,the useful attribute was selected,and the attributes that were not useful were discarded.It is important to consider the following points in this regard.

5.1.Cluster analysis

Cluster analysis has been widely used in natural and social sciences since the 1970s(Everitt B.2011).Q-mode and R-mode cluster analysis techniques are the most popular.Because this research is mainly aimed at geological parameter dimension reduction,rather than for sample reduction,R-mode cluster analysis is selected instead of Q-mode cluster analysis(Shi.2010).

The common seismic attribute extraction methods are:interlayer extraction and time window extraction along the layer.The layer thickness of the study area is thin,and lateral variation is not uniform.To avoid interference from the adjacent attribute,21 attributes were extracted along the selected horizon by using interlayer extraction.The maximum correlation of a single attribute is 0.0181,and the minimum is 0.0021(Fig.7).Both are unable to meet the prediction requirement,thus,it is necessary to perform cluster analysis on the 21 attributes,which in order to decrease the number of factors and reduce the amount of data(Shi,2010).

The cluster analysis process includes the following:

(1) Max-min normalization

Due to the difference between calculation and extraction methods of seismic attributes,the dimensional and order of magnitude variation ranges are large.Therefore,it is necessary to max-min normalization seismic attributes before clustering analysis,to eliminate systematic errors.A high level of normalization is adopted to handle the large magnitude range and variable seismic attribute parameters,and the seismic attribute values of different dimensions were normalized over the interval [0,1] 1using the following formula:

whereX2is the max-min normalized seismic attribute value,X1is actual seismic attribute value,andXmax，Xminare the maximum and minimum values of the seismic attribute,respectively.

(2) Clustering statistics

Clustering analysis is the process of finding a set of similar elements in a data set,called a“cluster”.Objects in the same cluster have a high degree of similarity,while those in different clusters are quite different (Wang et al.,2012).The number of classes (NOC) is an important factor that influences the classification process.A low NOC value possibly obscures important information.The recalibration distance of clustering reflects the correlation between lithology and seismic attributes.The closer the calibration distance is,the higher the correlation between these two attributes will be.To avoid interference between information caused by attribute clustering and inaccurate lithology prediction results,it is necessary to select one of several attributes with a high correlation coefficient.Searching for and identifying seismic attribute anomalies are often time consuming.

Fig.5.Synthetic seismograms of Hongshan well 3.

The Hungarian clustering algorithm was used to analyze the 21 extracted attributes,and the pedigree diagram of cluster analysis was obtained (Fig.8).According to principle that the distance between each center of gravity must be far,21 attributes are divided into four classes,as classification with more than five classes did not add more meaningful information.

To effectively identify similar potential reservoirs in neighboring non-drilled areas,it is necessary to optimize the properties that are significant to the sandbody in this area,namely,attribute optimization.There are many ways to optimize attributes.Cooke et al.(1999) proposed the empirical method.Zhao and Chen (2005)suggested the dual attribute optimization method of K-L transformation.The idea of the research involves all attributes in the calculation and gradually eliminates the attributes with poor correlation with the target through various optimization algorithms,to reduce attribute dimensionality,the disadvantage of methods are low efficiency and highly subjective factors.To ensure optimal dimensions and comprehensive information from the attribute combination,the Bucong algorithm was adopted in this study(Chen et al.,1997).The specific steps are as follows:

1) First,the attributes that were divided into four classes are ranked from largest to smallest in terms of correlation,namely class I＞class II＞class III＞class IV.

2) Second,in the attribute list classes (I1,I2,…,Ii),the correlation and standard error of a single attribute is calculated separately.According to the error analysis,the attribute with the smallest error is selected as single optimal attribute of class I,as I1.

3) Third,I1and class II attribute lists (II1,II2,…,IIi) are used to march each attribute,and select the property combination with minimum error as the second property of the optimal combination,and the other optimal combination of the two attributes is called the subprime attribute,i.e.,II1.

4) Fourth,to match each attribute of (I1,II1) and class III attribute

lists(III1,III2,…,IIIi),property combination with minimum error is selected as the third property of the optimal combination,and the third optimal combination of the three attributes is called the subprime attribute,i.e.,III1,which is the optimal combination of three attributes (I1,II1,III1).

5) Finally,the above steps were repeated until the matching correlation coefficient of the gravelly sandstone content was greater than 0.7.Then,the calculation was stopped,and the obtained attributes were the optimal attribute combination.

Through the above method,when the confidence levelR2=0.745 (Table 3),the calculation was stopped and the optimal attribute set was obtained (Table 4).It can represent variation in seismic signals related to the specific facies or lithology.

5.2.MRA prediction model

Linear analysis cannot adequately handle the inherent nonlinearity of the seismic data.Moreover,methods that rely on multivariate statistics are inflexible and require a large amount of statistical data (Social,2004).The predicted value can be obtained by establishing a mathematical model between the variable and the independent variable,and substituting the value of the independent variable into the regression model under the condition of conforming criteria (Muller and Stewart,2006).

The following step is the original data package obtained after correlation analysis by SPSS software(Table 4).Finally,eight kind of optimized seismic attribute values that are considered important reflections for the reservoir were input in the model,and the formula obtained using multiple regression,which corresponds to Eq.(2),is:

Fig.6.Logging facies assemblages of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.

From multiple regression,we have obtained matching the gravelly sandstone content(Table 5).The correlation between true value and matching value is 0.745(Fig.9).These predictions imply that the MRA is reliable.

Table 3 Correlation analysis of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.

Table 4 Correlation coefficient of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area(B value).

Table 5 Validation of multiple regression prediction model of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.

5.3.RMS attribute based on normal distribution

The RMS attribute is sensitive to the lithology and lithofacies,so it is important to interpret some special lithologic bodies(Kong and Bian,2001).According to the crossplot analysis,the correlation coefficient between the RMS attribute and the thickness of gravelly sandstone is 0.7802(Fig.10a),and the gravelly sandstone content is 0.8056 (Fig.10b).It is indicated that the RMS attribute is a quantitative positive correlation with not only the thickness of gravelly sandstone,but also the content of gravelly sandstone.

Fig.7.Scatter plots of correlation analysis between different seismic attributes and the matching content of gravelly sandstone.a.average instant phase attribute;b.total amplitute attribute.

Fig.8.Hierarchical cluster analysis of seismic attributes of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.Seismic attributes corresponding to serial numbers in the pedigree diagram are shown in Table 1.

Therefore,the RMS attribute is used to determine the boundary of sandbodies.How is the relationship between the RMS attribute and the sandbodies boundary built? This study adopts one of theprobabilities distributions types for random variables:normal distribution (Kong and Bian,2001).

The RMS attribute is normally distributed to obtain a normal distribution curve (Fig.11).The boundary of the sandbodies is confirmed.The section characterized with high amplitude in some locations,indicates the presence of sandbodies.Compared with RMS attribute without constrained(Fig.12a),the distribution range of sandbodies are constrained by the RMS attributes which is mormally distributed(Fig.12b).

In the RMS attribute diagram,the red and yellow colors represent the effective sandbodies,and the blue and green colors represent mud or few sandbodies.

6.Mapping the sedimentary facies

There is no doubt that seismic attributes provide unique interpretations regarding the reservoir parameters,and sandbodies distribution is reflected by the gravelly sandstone contour map.However,it cannot describe the internal structure of facies.Thus,it is necessary to map plane sedimentary facies based on geological models.

With the migration of the deposition center,the gravelly sandstone content is generally around 0.7 (Fig.13),and the facies zone evolves from two lobes to one lobe(Fig.14).The distribution of the fan delta plain consists of FA-1 with medium-thick gravelly sandstone and the logging facies assemblage is bell to box shaped and thick,which is in the area of h32042-Hong92-h15021-0718.The interdistributary channel consists of FA-4,with finger shaped response in the logging.

As the plain expands outward,the sandbodies become thin and disappear along the southern area,where the gravelly sandstone content decrease,about 0.5 (Fig.14),and the fan delta front develops.Among them,the subaqueous distributary channel inner zone is made up of FA-3,and FA-2 in the outer zone.The logging facies are usually boxed to bell shaped,and the amplitude is lower than that of the fan delta plain.The subaqueous interdistributary channel is similar to the interdistributary channel,but the color is darker and the layers are thinner,whereas in Hong43,CH7209,and0096.The mouth bar is composed of FA-5,and the log facies is funnel shaped and thick(Fig.14).

Fig.9.Correlation analysis between well locations and matching gravelly sandstone.

The pro-fan delta with gravelly sandstone content of less than 0.3,and is defined as black,silty sandstone with FA-6,and the logging facies assemblage is finger shaped,which is located around well HS3,00451,and h3429.

7.Conclusion

This paper has used multiple seismic attributes matching technology based on geological models quite successfully for mapping sedimentary facies,and the following basic conclusions are obtained:

1) Using the R-mode cluster analysis,including max-min normalization and clustering statistics,seven optimal attribute combinations were obtained.The MRA prediction model by the formula:y=-1.066-0.388r-0.221n+0.345l+0.037m-1.092q-1.553s+1.621p+4.374t was uniquely determined.

Fig.10.Correlation between RMS and the gravelly sandstone.a.Thickness of gravelly sandstone;b.Gravelly sandstone content.

Fig.11.The upper and lower limit after RMS limits the threshold value.Let the area of the whole closed normal distribution curve as 1,matching gravelly sandstone content is 0.46,namely the proportion of the area delineated by the red line in the figure,the corresponding RMS threshold value can be obtained.The lower limit a is 3247.341,the upper limit b is 13741.631,and the boundary of sand body is determined.

Fig.12.The change of boundary for the RMS attribute.a.without constrained by the normal distribution;b.Constrained by the normal distribution.

Fig.13.Matching contour map of gravelly sandstone.

2) Through RMS attribute based on normal distribution,the threshold value of the upper and lower limits for the RMS were obtained.The lower limit a is 3247.341,the upper limit b is 13741.631,and the boundary of sand body is determined.So,the RMS attributes distribution could more truly reflect the response of the sand body distribution.

3) The seismic geomorphology workflow used in this research successfully integrated seismic interpretation attribute analysis,geostatistical simulation of well log facies,and lithofacies assemblages to build a geologic scenario for mapping sedimentary facies.If gravelly sand bodies are analyzed in a more comprehensive way,then the exploration risks can be reduced.Multiattribute matching technology have been proven in this research to provide valuable information between good and poor reservoirs.

Fig.14.Sedimentary facies of the 3rd Member in the Lower Jurassic Badaowan Formation,Hongshanzui area.Notes:B-bottom,T-top,Cc-Coarsen-grained conglomerate,Mc-Middle-grained conglomerate,Fc-Fine-grained conglomerate,Cs-Coarsen-grained sandstone,Ms-Middle-grained sandstone.

Acknowledgements

This study was supported by the National Natural Science Foundation of China (41902109),Tianshan Youth Program(2020Q064),National Major Projects (2017ZX05001004),Tianshan Innovation Team Program (2020D14023).It also received aid from the research team at the China University of Geosciences,Beijing.