Concealed porphyry delineation based on nonlinear three-dimensional densitydifference inversion: An example in the Beiya mine area, Western Yunnan, China

Jian Yang, Sheng-xian Liang, Qiao Wang, Wei Zhang, Jing Guo, Guo-zhong Liao

Chengdu Center, China Geological Survey, Ministry of Natural Resources, Chengdu 610081, China

Keywords:

Three-dimensional density inversion

Concealed porphyry

Gold deposit

Mineral resources exploration engineering

Beiya mine area

Yunnan Province

China

A B S T R A C T

Intermediate acid-complex rock masses with low-density characteristics are the most important prospecting sign in the Beiya area, of western Yunnan province, and provide a physical basis for good gravity exploration.It is usually difficult to obtaining solutions in connection with actual geological situations due to the ambiguity of the conventional gravity-processing results and lack of deep constraints.Thus, the three-dimensional (3D) inversion technology is considered as the main channel for reducing the number of solutions and improving the vertical resolution at the current stage.The current study starts from a model test and performs nonlinear 3D density-difference inversion called “model likelihood exploration”, which performs 3D inversion imaging and inversion of the known model while considering the topographic effects.The inversion results are highly consistent with those of the known models.Simultaneously, we consider the Beiya gold mine in Yunnan as an example.The nonlinear 3D densitydifference inversion technology, which is restricted by geological information, is explored to obtain the 3D density body structure below 5 km in the mine area, and the 3D structure of the deep and concealed rock masses are obtained using the density constraints of the intermediate-acid-complex rock masses.The results are well consistent with the surface geological masses and drilling-controlled deep geological masses.The model test and examples both show that the 3D density-difference nonlinear inversion technology can reduce inversion ambiguity, improve resolution, optimize the inversion results, and realize“transparency” in deeply concealed rock masses in ore-concentrated areas,which is useful in guiding the deep ore prospecting.

1.Introduction

The Beiya gold-polymetallic deposit situated in the western margin of the Yangtze Platform in western Yunnan Province, China, is in the southeastern margin of the Tibetan Plateau (Yan JY et al., 2018).It is an important and typical alkali-rich porphyry-hydrothermal- and contact metasomatictype polymetallic deposit in the alkali-rich porphyry polymetallic metallogenic zone along three rivers, namely,Jinsha River, Nu River, and Lancang River, in the southwestern areas.By the end of 2016, the proven reserves in the gold mine reached more than 160 t, with a total iron-ore reserve that exceeds 3×107t.It is already a super large gold ore and a large-sized iron ore with many associated metals such as silver, lead, zinc, and copper (Yan JG et al., 2003; He ZH et al., 2013).In recent years, following in-depth geological and geophysical research and exploration,mineralization has occurred in Xiangshuihe, Huangcaopo, and Ma’anshan (Sun Y et al., 2010; Niu HB et al., 2015; Sun L et al., 2018), which shows that great prospecting potential exists beyond the deep and known deposits.The known deposits or mineralization sites in the Beiya Mine area are closely related to the Himalayan alkali-rich porphyry (quartz syenite porphyry, ivernite, etc.).The ore-controlled rock masses,tectonics, and strata are integrated as one, and from the center of the rock masses to the areas beyond consists of porphyry copper-gold minerals (body), banded skarn iron-copper-gold ore masses, and vein-shaped and vesicular-filled penetrationtype iron-gold ore body.Therefore, It is essential to understand the distribution characteristics of the deep-rock mass, strata, and structural distribution for deep-metallogenic prognosis.

Alkali-rich porphyry rock masses are low-density masses compared with peripheric the Beiya Formation carbonate rocks and Mount Emei basalt.Using massive gravity data, we aim to search for local low-gravity abnormity which resulted from alkali-rich porphyry rock masses to realize direct prospecting.According to previous studies (Sun WK et al.,1997; Wang JY, 1998; Fedi M et al., 1999; Cao LM et al.,2011; Deng Z et al., 2012), gravity data are an added component to the gravity field generated by an underground field source.When an underground object is complicated, the two-dimensional (2D) data-processing method based on fieldsource separation and different height extensions is poor in terms of both fineness and resolution, thus lacking profound information.In recent years, gravity and magnetic inversion has developed from 2D to three-dimensional (3D) with the progress in geophysical inversion theory and computer technology.The 3D inversion is more suitable for actual distribution status of underground geologic bodies than the 2D inversion; thus, the 3D inversion is better than the 2D section inversion (Nabighian MN et al., 2005; Hirt C and Kuhn M, 2014).Williams (2008), used a theoretical model,and proved the effectiveness and feasibility of the gravity and magnetic 3D inversion with geological restraint, and attained good results in the exploration of the nickel deposit in Perseverance, Western Australia.

The 3D inversion of gravity data is generally divided into two categories, namely, physical and form inversions.The latter often result in the single and multiple partial gravity and magnetic abnormalities and inversion interpreters need to have adequate experience and waste a lot of time and efforts to separate them..Some other inversion algorithms such as the fast Fourier transform (Young Hong Shin.et al., 2006; Zhang WH et al., 2019) and apparent density imaging methods (Hao TY et al., 1998; Xu SZ et al., 2007) are fast in terms of inversion speed, but the inversion results are not sufficiently accurate for complexly distributed density bodies.Meanwhile,the physical inversion suffers from two weaknesses, namely,low calculation efficiency (Hearst RB and Morris WA, 2001;Yao CL et al., 2007; Tai ZH et al., 2011) and distortion in the equations (Li YG and Oldenburg DW, 1996, 1998; Wang YS et al., 2018).The current work introduces the nonlinear 3D gravity density-difference inversion algorithm named “model likelihood exploration,” which performs a 3D density model inversion test using constraints.We employ the 1:250000 gravity data in the Beiya Mine area to carry out the 3D density inversion.The model test and example results show that the nonlinear 3D density-difference inversion can reduce the inversion ambiguity, improve the resolution, optimize the inversion results, identify the 3D patterns and stratum 3D distribution of the rock masses in the Beiya mine, realize“transparency ” of the target, and provide information and basis for deep-mining prospecting.

2.Geological and geophysical characteristics of the study area

2.1.Geological background of the study area

Beiya mine is located to the east of the joints between three Class I tectonic units, namely, Dege-Zhongdian,Yangtze, and Lanping-Simao landmasses (Fig.1) between the Jinsha River-Honghe, Binchuan-Chenghai, and Lijiang-Muli fractures.The mining area with continental tensioned-rift environment is located in a depression and fold zone of the Lijiang Platform edge west of the Yangtze quasi-platform,which is classified as a Dali-Lijiang centrally mineralized area in the middle and southern sections of the Yangtze metallogenic zone.The outcropping stratum in the study area includes the Permian upper Mountain Emei formation (P2β),lower Triassic Qingtianbao formation (T1q), Triassic middle Beiya formation (T2b), and Quaternary.The Triassic Middle Beiya and Qingtianbao formations entirely contact each other,which mainly consist of dolomitic-sand clastic limestone,dolomite, iron-sand clastic limestone, argillaceous limestone,and vermicular limestone, among others, as the main ore layer of the mine.This mine has an abundant ore-deposit level,which indicates the distribution of the ore body to some extent.The lower Triassic Qingtianbao formation outcropping east of the study area and distributed in a south-north direction, is integrated with the Mount Emei basalt and it has a yellowish green-dark gray color and contains green thin-tomedium layered and hornfelsed graywacke.The upper Permian Mount Emei basalt outcrops in the eastern part of the study area and is composed of dark green basalt, which is very much weathered and fractured on the surface.The outcropping magmatic rock in the Beiya area is dominated by hypabyssal-intrusive alkali-rich porphyry.The Beiya syncline is the main fold in the mine area, located at the southern warping end of the Songgui synclinorum, which is a secondorder structural units of the Heqing-Songgui synclinorum.The axis direction is north-northeast, which is a wide and gentle short-axis syncline.The fault mainly comprises two formations (one is almost in the north-south direction and the other is almost in the east-west direction), which are distributed on both wings of the Beiya syncline.The almost north-south fault formation is an internally controlled mine fault, and the almost east-west fault group is an ore-breaking fault.

Fig.1.Geological map of the Beiya area (modified from Regional Geology of Yunnan Province, 2012).

2.2.Physical properties of the rock and ores

To study the different stratum- and lithology-density characteristics in the Beiya area, we collected 305 rock samples and cores, which included limestone, dolomite,sandstone, ballast, quartz syenite porphyry, and other.The density parameters were determined.Table 1 lists the statistical results.

Table 1.Statistics of the physical parameters of rocks (ores).

The stratum density in the area gradually increases from new to old strata.The Middle Triassic Beiya formation carbonate is widely distributed in the area with a comprehensive density of 2.64×103kg/m3.The Quaternary clay-rock density is 1.92 × 103kg/m3.The Lower Triassic Qingtianbao formation features a high density of siltstone at 2.73×103kg/m3, which is significantly higher than the previously measured siltstone.The quartz syenite porphyry has the lowest density of 2.24×103kg/m3except for the Quaternary.The basalt density of Mt Emei is the highest at 2.87×103kg/m3.In summary, we compared the alkali-rich porphyry rock masses with the Beiya formation carbonate rocks and Emei basalt and determined that it had a lowdensity body, which was slightly more than the surface Quaternary density value (Fig.2).Therefore, the low-gravity anomaly in this area was mainly caused by the geological structure associated with the rock-mass distribution and magma activity, which could describe 3D distribution characteristics.

Fig.2.The comparision map of density of alkali-rich porphyry with carbonate rock and Emei basalt

2.3.Characteristics of the gravity field in the Beiya Mine

Fig.3a shows that the Bouguer gravity anomaly in Beiya Mine is negative (from -252 to -238 × 10-5m/s2), and the measured area reaches up to 42 km2.The figure shows that the gravity anomaly in the northern part of the area is low, and that in the southern part is high.The anomaly is distributed in an arc-gradient zone from north to south.The Wandongshan Mine section is in the middle and northern low-anomaly zones, and the Hongnitang and Jingouba Mine sections are located in the anomaly-gradient zone.

Fig.3.Gravity anomaly of Beiya mine.a-Bouguer anomaly; b-residual anomaly; c- regional anomaly

The Bouguer gravity anomaly reflects a macro geological phenomenon and is used to understand the distribution information of regional-, local-, deep-, and shallow-density bodies(Nowell D.A.G, 1999; Kou XY et al., 2006; Liang XT et al., 2016; Jiang XD et al., 2019).The Bouguer gravity anomaly is subjected to anomaly separation, and the gravityanomaly separation should be performed according to different geological conditions using different separation methods such as tendency analysis, matched separation, and other methods for comparison(Li Y et al., 2018; Li QL et al.,2019).The regional field and residual anomaly obtained by a moving-average field window of 1.3 km × 1.3 km is consistent with the spatial position and pattern alignment of an outcropped geological body, and we adopt this method to calculate the remaining regional anomalies, as shown in Fig.3b and Fig.3c.According to the analysis of the geological data and physical parameters, the gravity anomaly is known to vary because of the underlying rock-mass scale and basalt.Thus, the height and weight anomaly is caused by the basal ballast, whereas the middle regional low and heavy anomalies are caused by the deep underlying rock masses.According to the residual-gravity-anomaly map, local low anomalies formed by the low-density rock masses are observed in the Wandongshan and Hongnitang mine sections.However, the local anomaly is quite different from that of the actual geological conditions at the basal of the eastern basalt in the mine area because the data-processing method based on the field-source separation and extensions of various heights are poor in terms of fineness and resolution and lack constraints.Thus, the results are inconsistent with those in actual geological conditions.

3.Method and principle of nonlinear 3D densitydifference inversion

The 3D inversion of gravity data is generally divided into two categories: physical and form inversions (Qi G et al.,2012, 2014; Yan YJ et al., 2014; Yan L et al., 2018).The latter is time consuming and requires much effort and inversion interpreters need to have adequate experience,which result in the difficulty in separating the single and multiple partial gravity and magnetic abnormalities.Some other inversion algorithms such as the fast Fourier transform(Young Hong Shin.et al., 2006; Anderson ED et al., 2014)and apparent density-imaging methods (Xu SZ et al., 2007)have fast inversion speed.However, the inversion results are not sufficiently accurate for complexly distributed density bodies.The physical inversion suffers from two weaknesses,i.e., low calculation efficiency (Boszczuk P et al., 2011; Yao CL et al., 2007; Boulanger O et al., 2001) and distortion in the equations (Li Y and Oldenburg DW, 1996, 1998).In the present work, we introduce a nonlinear 3D density-difference inversion algorithm named “model likelihood exploration”,which was proposed by Spanish scholar Antonio GC (2002).The gravity data are imaged through 3D inversion, and the geological effects are considered in the inversion.

The nonlinear 3D gravity density-difference inversion“model likelihood exploration ” algorithm divides the underground into fixed units with several dimensions (Fig.4).The density difference in the element is a fixed value, and what distinguishes the unit division from the traditional division is that the size of the element vertically varies and increases with the increase in the horizontal depth.It shows a“pyramid ” profile on a 2D section map, which greatly decreases the number of units and enormously improves the 3D inversion solution speed.Second, we know that the amplitude of the surface gravity and magnetic data attenuates with the increase in the field-source depth.Thus, the “model likelihood exploration” algorithm uses a random search to avoid this effect while avoiding large matrix storage and solving linear equations.

Fig.4.Three-dimensional gravity forward modeling cell sketch map

We assume the following: (1) there aren(i=1,2,···N)surface observation points, (2) the model space is discretized into a fixed-size cell element, (3) the density value of each element is equal toσj(j=1,2,···M), and (4)gravity response Δ giat surface observation pointPi(xi,yi,zi)can be expressed as

whereaijis the kernel function of the gravity response of thejth element at theilocation on the surface, which is irrelevant to the specific-density value of the element.The analytical formula is expressed as follows (Antonio GC et al.,2000, 2002; Chen ZX, 2012; Chen H et al., 2015):

The inversion is based on the surface-observed gravity response to calculate the field-source distribution.Typically,the number of surface observation data is less than that of the models, i.e.,N＜M.In other words, the inversion problem is underdetermined.In addition to considering data covariance matrixQD, covariance matrix elements(ei(i=1,2,···n))are considered, which are standard deviations of the gravity observations.

We also consider remaining solution model setmin addition to the abovementioned model gridding and data errors.In fact, we can consider some prior density difference(positive, negative, and zero density differences) to fill a certain cell element.After the filling, the model anomaly further increases to build an entire model anomaly.

The specific algorithm is described as follows.Starting with a certain (random) element, a possible density difference must be explored so that the corresponding surface gravity response is ideally fitted with the observations at certain proportion factorf1.Then, we start searching the second element to continue, wherefn＞fn-1···＞f1.Whenfn=1, the algorithm stops, and the inversion results are exported.The mathematical expression of the algorithm is expressed as follows:

We assume that the initial model isand

Then, if we search for stepk+1 and the (k+1)th unit, the density of the firstkth unit iswhere Δρjis the varying value of the density, which can be positive,negative, or zero.Then,

In the equation,Riis the regional gravity anomaly, which can be obtained using the first-order polynomial, andfis an unknown scaling factor.Their specific methods can be found in the literature (Antonio GC, 2002).To reduce the inherent non-uniqueness and instability of the inversion, an objective function of the inversion can be expressed using the Tikhonov regularization principle, which takes into account the datafitting and model-constraint terms, as follows:

wherem=(Δρ1,Δρ2,···Δρj)T, λ is a regularization factor,andQDandQmare the data and model covariance matrices,respectively.

Thus, at stepk+1, the algorithm exploration corresponds to the likelihood of positive Δ ρjand negative Δ ρjof the element, which determines thatEk+1is a minimum.Newis calculated.The process is repeated untilf=1 where the algorithm is terminated.

4.Theoretical-model test

Through theoretical-model experiments, we can understand the application effects of the adopted densitytomography inversion method(Meng XH et al., 2012; Meng ZH et al., 2018).We assume that the underground half-space contains three residual-density bodies.The specific model parameters are listed in Table 2.First, the surface Bouguer anomaly of the forward computation model abnormally responds, and the “model likelihood exploration” algorithm is utilized for inversion of the solution model for comparison with the actual model.The sampling density of the forward computation data is 20×50 m.Figs.5 and Fig.6 show a set of theoretical models and the inversion results of the 3D density.We can see that the slices at various depths of the inversion result are comparable with those of the real model, which basically reflect the density contours of the real model.The solution of the likelihood exploration algorithm of the inversion model is reliable and can be used in practical application.

Table 2.Theoretical model parameters for three-dimensional inversion.

5.Three-dimensional density inversion examples in Beiya mine

Considering the effective resolution of the 1:10000 gravity data and the inversion work objectives in the Beiya fine sectioning area, the adopted inversion model is sectioned in the north, east, and vertical directions.The dimensions of the sectioned blocks are 100 m×100 m×100 m.The data area is 46.15 km2, and the depth is 3 km below the ground.The inversion physical constraints and parameters are listed in Table 3.

Table 3.Three-dimensional inversion parameter table of the Beiya fine sectioning area.

In the inversion process, the surface property model is used, and the surface geological information is used as the empirical data, which can be easily obtained.In this experiment, the surface geological map of the mine area is obtained, and the density parameters of all target geological bodies are set (Fig.7).The inversion process is forcefully restrained.

Fig.7.Geophysical model of the surface geology.

5.1.Three-dimensional density inversion results

According to the physical characteristics of the mining area, the rock masses feature a low density.Therefore, the low-gravity anomalies generated by the rock mass can be extracted according to the density constraint.The 3D inversion density model reveals the 3D distribution characteristics of the rock mass (Fig.8).Fig.9 shows the 3D effect obtained from the mine 3D density inversion results.Fig.10 shows the extracted surface gravity anomaly after the 3D density inversion.Compared with the surface geology in the mine area (Fig.1), the high-gravity anomaly is generated by the ballast rocks east of the mine area.The low-gravity anomaly generated by the deep rock masses in the Wandongshan and Hongnitang Mine sections correspond well.In addition, from the point of view of the extracted surface gravity anomaly, we assume that the Hongnitang and Wandongshan rock masses are integrated into one at this depth.

Fig.8.Three-dimensional density structure visualization

Fig.9.Surface gravity anomaly of the 3D inversion

Fig.10.Three-dimensional density structure of the low-density body.a-Overall structure; b-local structure.

To verify our inversion results and confirm if the inferred deep concealed rock masses are accurate, we use the 3D density inversion results to consider the line55 exploration profile of the Hongnitang mine section in the Beiya mine area(Fig.11a, b).Compared with the actual geological profile(Fig.12), we find that the entire density structure is characterized by “low west and high east”.The east highgravity assumption is caused by the ballast, and the western“funnel” low-gravity anomaly is considered as caused by the concealed Hongnitang and Dashadi porphyry.However,because of the medium- and low-density characteristic interruption of the Triassic Qingtianbao formation sandstone,some differences are observed between the low-gravityanomaly and actual buried porphyry forms.Clearly, by using the aforementioned nonlinear 3D gravity density-difference inversion “model likelihood exploration” algorithm, the 3D density structure is obtained, which can accurately indicate the position, size, and deep extension of the buried porphyry.Further, clues for direct and indirect prospecting in the area can be provided.

Fig.11.Three-dimensional density inversion slice of line 55 in the Hongnitang section.a-3D density structure model; b-inversion profile of line 55.

Fig.12.Line55 prospecting line profile map in the Hongnitang section.

5.2.Three-dimensional distribution characteristics of stratum

From the 3D density model slices (Fig.13a), we can see that the Beiya Mine density is distributed in an arc-gradient distribution.The mine area is abnormally low in the north,whereas the southern anomaly is relatively high.It is distributed in an abnormal arc-gradient distribution from north to south.The eastern high-density body is caused by the basal basalt, the middle low-density anomalies are caused by surface Quaternary and deep porphyry rock masses, and the western density is gradually reduced from deep to shallow stratum.From the 3D density structure, the mine structure is extended in a nearly south-north direction, whereas the Beiya syncline is obvious.The syncline kernel is located in the Wandongshan and Hongnitang Mine sections.Fig.13b shows the 3D distribution characteristics of the high-density body in the Beiya Mine volume (＞1.73 g/cm3), which mainly reflect the ups and downs of the deep basal ballast formation.We can see that the basal basalt formations are uplifted in a north-south direction in the eastern Bijiashan area with large deep extensions.The alignment is consistent with the formation tectonics of the Beiya Mine.These formations are sporadically distributed in blocks in the western and central parts, and the surface only outcrops at the western boundary of the mine area, showing the tectonic action in the mine area,such as folds and extrusions, which have seriously destroyed the basal ballast formation.

Fig.13.Three-dimensional density inversion map in the Beiya mine.a-Inversion slice profile; b-Three-dimensional density structure of highdensity body

5.3.Three-dimensional distribution characteristics of porphyry

The low gravity anomalies generated by the porphyry are extracted according to the low-density characteristics of the porphyry and surface constraints.The 3D inversion density model reveals the 3D distribution characteristics of porphyry(Fig.14).We map clearly shows that porphyry exists in Wandongshan, Hongnitang, and Dashadi, which corresponds very well to the porphyry controlled by the exploration line.In addition, from the 3D morphology of the porphyry, we infer that the Hongnitang porphyry and Wandongshan porphyry are deeply connected through Beiya.Further exploration and verification are needed in the later stage.

Fig.14.Three-dimensional distribution of porphyry in Beiya Location.

In addition, to study the relationship between the porphyry and mineralization in the Beiya Mine area, the exploration line profiles of Wandongshan and Hongnitang Mine sections are cut using the 3D gravity inversion, and the known ore bodies are projected onto the 3D gravity inversion slice map(Fig.15).We find that both iron and gold deposits in Wandongshan and iron deposit in Hongnitang are closely related to the porphyry.After a comprehensive study, we consider that minerals are not necessarily present where porphyry exists in the Beiya Mine area.However, porphyry may be present where minerals exist.Porphyry is closely related to mineralization.Therefore, the main purpose of prospecting in the Beiya Mine area and its periphery is to search for porphyry, which also demonstrates a large response to the 3D inversion of gravity.

Fig.15.Three-dimensiona density inversion imaging (ore body).

6.Conclusions

According to the inversion theoretical test for nonlinear 3D density difference and the gravity 3D inversion results in Beiya Mine, we can draw the following conclusions.

(i) The nonlinear inversion results of the 3D density differences using the surface geological information are obviously better than those performed by the 2D dataprocessing methods based on the field-source separation and extensions at different heights, which greatly reduce the datafitting difference and improve the correlation coefficient of the inversion model under a real situation.In the inversion process, the surface and drilling empirical information is called surface because it can provide restraints in the depth direction.It can greatly improve the local anomalyinterpretation precision and vertical resolution, which effectively reduces ambiguity.

(ii) According to the study on the relationship between the fine physical properties and lithology, the nonlinear gravity 3D density-difference inversion study of the Beiya Mine area is used as an example, which is supported by additional surface empirical information.Some results are obtained.Most of the quartz syenite porphyry, ivernite, and rocks with low-density characteristics can be clearly observed, and a corresponding relationship exists compared with the known surface outcropping and drilling control.The extensive basal ballast uplifting is consistent with the distribution pattern and alignment of the high-density body.

(iii) We consider that minerals do not necessarily exist where porphyry is present in the Beiya Mine area.However,porphyry exists where minerals are present.Porphyry is closely related to mineralization.Therefore, the main purpose of prospecting in the Beiya Mine area and its periphery is to search for porphyry, which also demonstrates a great response to the 3D inversion of gravity.

(iv) The physical difference of some rock units in this test is small, such as the Qingtianbao formation sandstone and Beiya formation carbonate rocks, which cannot be identified from the inversion results.To solve the current problems of inadequate gravity-based inversion, we need to improve the inversion precision, reduce the inversion ambiguity,strengthen the system measurement, and study the physical properties to learn the physical combination of different rock masses and the corresponding physical properties of the strata.

Acknowledgement

The authors would like to thank the China Geological Survey (DD20190033) and National Natural Science Foundation (41804144) for the financial support, Yunnan Gold and Mineral Group Co., Ltd.for providing the original geological information, and the reviewers for providing valuable comments.