تخمین کریجینگ مقاطع به‌منظور مدلسازی و ارزیابی ذخیره معدنی (مورد مطالعه: کانسار سرب و روی عمارت)

نوع مقاله : مقاله پژوهشی

نویسندگان

گروه مهندسی معدن، دانشکده مهندسی علوم زمین، دانشگاه صنعتی اراک، اراک، ایران

چکیده

مقدمه
مدلسازی توزیع فضایی عیار و تخمین ذخیره کانسار، مهم­ترین مسئله و هدف اصلی عملیات اکتشاف است. این نوع مدلسازی با استفاده از روش­های مختلفی انجام می­شود که به میزان، نوع و روش کارهای اکتشافی انجام­ شده و اطلاعات اکتشافی در دسترس، بستگی دارد. در شرایطی که ماده معدنی یک کانسار بصورت قائم و یا با شیب زیاد باشد؛ وضعیت ژئومورفولوژی منطقه به­گونه­ای باشد که منطقه کوهستانی بوده، دارای ارتفاعات بلند و توپوگرافی خشن باشد؛ استقرار دکل و دستگاه عملیات حفاری و در نتیجه انجام حفاری اکتشافی مشکل یا غیرممکن باشد و نیز میزان ارتفاع و حجم روباره (باطله سطحی) پوشاننده روی سطح ماده معدنی خیلی زیاد باشد، عملیات اکتشاف کانسارها بوسیله کارهای اکتشافی زیرزمینی همانند تونل­های اکتشافی صورت می­گیرد. در واقع در چنین شرایطی در مورد اغلب مواد معدنی، در راستای قائم پیوستگی وجود دارد و تغییرات، بیشتر در راستای سطوح افقی است. در پژوهش حاضر مدلسازی بلوکی و تخمین ذخیره زمین­آماری کانسار سرب و روی عمارت واقع در استان مرکزی با روش لاگ­کریجینگ با رویکردی متفاوت و تکنیکی ظریف صورت گرفته است. برای دست­یابی به هدف ابتدا پردازش، مدلسازی و تخمین زمین­آماری داده­های عیارسنجی مجموع سرب و روی حاصل از حفر تونل­های اکتشافی در ترازهای ارتفاعی مختلف به­طور جداگانه و بصورت دوبعدی صورت گرفته، سپس به­کمک نتایج حاصل از این مرحله، فرآیند تخمین برای کل کانسار بصورت سه­بعدی انجام شده است. چنین رویکردی برای تخمین زمین­آماری عیار و ذخیره یک کانسار که با کارهای اکتشافی زیرزمینی مورد اکتشاف قرار گرفته است، قبلاً در هیچ پژوهشی گزارش نشده است.
مواد و روش­ها
کانسار مورد مطالعه و داده­های اکتشافی در دسترس
کانسار سرب و روی عمارت در فاصله حدود 45 کیلومتری جنوب­غرب شهرستان اراک، بین طول­های جغرافیایی ¢30 °49 تا ¢45 °49 شرقی و عرض­های جغرافیایی ¢45 °33 تا ¢0 °34 شمالی در منطقه­ای با ارتفاع ۲۱۸۰ متر از سطح آب دریا قرار دارد. توپوگرافی منطقه عمارت بسیار ناهموار می­باشد؛ چینه­شناسی یکنواخت، چین­خوردگی شدید، فقدان سنگ­های آذرین و استراتی­باندبودن از ویژگی­های زمین­شناسی منطقه است که جهت چین­خوردگی­ها مطابق با روند چین­خوردگی زاگرس می­باشد.
کانسار سرب و روی عمارت بر روی پهنه تکتونیکی سنندج- سیرجان و کمربند فلززایی سرب و روی ملایر- اصفهان واقع شده است. واحدهای زمین­شناسی منطقه عموماً واحدهای رسوبی کربناته و شیل و مارن متعلق به دوره ژوراسیک تا کرتاسه همراه با چین­خوردگی و گسلش است که از دامنه­های شیب­دار تشکیل شده­اند. به­طورکلی کانسار سرب و روی عمارت یک ناودیس شکنجی با امتداد شمال­غرب - جنوب­شرق است که طول آن 5/1 کیلومتر و عرض آن از 250 متر تا 850 متر متغیر است. در کانسار عمارت، ماده معدنی سرب و روی درون یک لایه سیلیسی آهکی در مرز یک نوار سنگ­آهک توده­ای خاکستری تیره­رنگ در کمرپایین و یک لایه شیلی کرتاسه در کمربالا قرار دارد. این تشکیلات متعلق به کرتاسه تحتانی یا میانی است. در کانسار سرب و روی عمارت فعالیت­های زیادی شامل حفر تونل­های اکتشافی و استخراجی افقی و دستک­ها با طول مجموع 11000 متر در شش تراز ارتفاعی با اختلاف ارتفاع حدود 10 متر صورت گرفته و داده­های عیارسنجی عناصر سرب و روی 1238 نمونه برداشت­شده از تونل­ها در دسترس است.
روش کریجینگ مقاطع
از نقطه­نظر ابعادی، در مورد داده­های برداشت­شده دوبعدی همانند داده­های عیارسنجی نمونه­های سطحی، عملیات پردازش و تخمین زمین­آماری با روش کریجینگ، بصورت دوبعدی انجام می­شود. مواردی هم وجود دارند که در آنها اگرچه برداشت داده­ها بصورت سه­بعدی است، اما بهتر است تخمین زمین­آماری ابتدا در راستای یک سری سطوح دوبعدی انجام شود، سپس نتایج تخمین حاصل از سطوح مختلف، با یکدیگر تلفیق شده و مدل سه­بعدی نهایی تولید گردد. به این نوع از تخمین کریجینگ، کریجینگ مقطعی گفته می­شود که می­تواند در راستای یک سری مقاطع افقی یا قائم صورت گیرد؛ البته استفاده از آن در مورد مقاطع افقی منطقی­تر و مطلوب­تر است. از آنجایی ­که در این­گونه موارد در فاصله بین دو سطح متوالی و مجاور هم، عملاً هیچ نمونه و در نتیجه هیچ­گونه داده­ای وجود ندارد و تمامی داده­های مورد پردازش و تخمین، مربوط به نمونه­های برداشت­شده از سطوح دوبعدی هستند، بنابراین میزان ارتباط فاصله­ای و فضایی داده­های هر سطح با یکدیگر بیشتر است و میزان پیوستگی کانسار در هر سطح، بهتر نمایان می­شود.
روش لاگ­کریجینگ
اگر توزیع آماری داده‌های مورد استفاده از نوع نرمال نباشد، نمی‌توان روش­های تخمین کریجینگ خطی را به­کار گرفت، زیرا در این حالت اثر تناسب واریانس با میانگین وجود خواهد داشت و ناهمسانگردی دروغین در واریوگرام­های امتدادی پدیدار می­شود. در این صورت بهتر است که داده‌ها با یک روش تبدیل مناسب همانند تبدیل لگاریتمی نرمال شوند تا بتوان روش­های خطی را برای تخمین به ­کار برد. بعد، عملیات تخمین با روش کریجینگ معمولی (یا ساده) بر روی لگاریتم داده­ها اعمال می­شود و سپس مقادیر تخمین­زده­شده با یک تبدیل معکوس، به مقادیر واقعی تبدیل می­گردد. به این روش کریجینگ غیرخطی، لاگ­کریجینگ گفته می­شود.
نتایج و بحث
بر اساس پردازش آماری داده­های عیارسنجی مجموع سرب و روی برای ترازهای ارتفاعی 2032، 2024، 1998، 1988، 1978 و 1968-1964 متر به­طور جداگانه، توزیع داده­ها از نوع لاگ­نرمال بود که با تبدیل­های لگاریتمی دو و سه­پارامتری به حالت نرمال تبدیل شدند. همچنین رسم نمودار انحراف­معیار در مقابل میانگین داده­های عیارسنجی مجموع سرب و روی ترازهای ارتفاعی مختلف کانسار، اثر تناسب واریانس با میانگین را نشان داد. برای تجزیه و تحلیل ساختار فضایی منطقه، واریوگرام­های امتدادی افقی برای هر تراز ارتفاعی به­صورت مجزا با آزیموت­هایی در بازه 25 تا 135 درجه ترسیم شد و دو واریوگرام افقی بهتر در راستاهای عمود بر هم انتخاب شدند. با توجه به لاگ­نرمال بودن توزیع داده­های عیارسنجی، به­منظور اجتناب از پیدایش ناهمسانگردی دروغین، واریوگرافی برای مقادیر داده­های عیارسنجی مجموع سرب و روی تبدیل­یافته انجام شد. تمامی مدل­های تئوری واریوگرام منطبق بر واریوگرام­های تجربی، از نوع کروی بوده و واریوگرام­های امتدادی در جهات مختلف دارای سقف یکسان ولی شعاع تاثیر متفاوت هستند؛ بنابراین منطقه مورد مطالعه دارای ناهمسانگردی از نوع هندسی می­باشد. به ­کمک واریوگرافی ترازهای ارتفاعی مختلف و تعیین شعاع­های بیضی جستجو، تخمین عیار هر تراز ارتفاعی با روش زمین­آماری لاگ­کریجینگ مقاطع برای بلوک­های 10×10 متر انجام شد و نقشه هم­عیار حاصل از فرایند تخمین ترسیم گردید. مطابق نقشه­های 
هم­عیار در اغلب ترازهای ارتفاعی کانسار سرب و روی عمارت، عیار ماده معدنی در نیمه شرقی از نیمه غربی بالاتر است. سپس با رویکرد اعتبارسنجی متقابل به روش جک­نایف و ترسیم نقشه­های خطای تخمین کریجینگ هر تراز ارتفاعی، اعتبار فرآیند تخمین اثبات شد.
 
در این روش اعتبارسنجی بر اساس مدل واریوگرام، هر بار یکی از داده­های ورودی (معلوم) با استفاده از نمونه‌های همسایگی اطراف آن نمونه به روش کریجینگ تخمین زده می‌شود و مقادیر تخمینی با مقادیر واقعی مقایسه می‌شوند. به ­عبارت دیگر هر مقـدار معلومی با فـرض ایـنکه مقدار آن مجهول است، تخمین زده می­شود. در این راستا، میزان ضریب تعیین رگرسیون بین مقادیر عیار واقعی با تخمینی در تمامی موارد به جز تراز ارتفاعی 2024، مقداری بیش از 5/0 می­باشد که نشان­دهنده همبستگی خوب بین داده­هاست. بنابراین تخمین صورت گرفته از درجه اعتبار مطلوبی برخوردار است. برای ایجاد مدل بلوکی کانسار، از یک مدل شبکه­ای اولیه به ­ابعاد 265×530×1110 متر با سلول­های 10×10×10 متری و داده­های عیارسنجی تخمینی ترازهای ارتفاعی مختلف به ­کمک الگوریتم عکس فاصله وزن­دار پیشرفته استفاده شده است. ابعاد بلوک­ها براساس اندازه و ابعاد و مساحت سطح گسترش ماده معدنی در هر تراز ارتفاعی و نیز فاصله بین ترازهای ارتفاعی کانسار، 10×10×10 متر انتخاب شد. در الگوریتم عکس فاصله وزن­دار پیشرفته امکان وزن­دهی فاصله با توان متفاوت در جهات مختلف وجود دارد. در این مورد به داده­های در راستاهای افقی (تونل­های اکتشافی در هر تراز ارتفاعی) به­ دلیل تغییرپذیری بیشتر، وزن دو و در راستای قائم (فاصله بین ترازهای ارتفاعی مختلف) وزن یک نسبت داده شد. درنهایت رده­های مختلف ذخایر قطعی، احتمالی و ممکن کانسار به­ازای هفت عیارحد 4، 5/4، 5، 5/5، 6، 5/6 و 7 درصد تعیین شد. براساس نمودار عیار- تناژ کل کانسار برای رده­های مختلف ذخایر قطعی، احتمالی و ممکن، نمودارهای ذخیره قطعی و متوسط عیار نظیر ذخیره قطعی کانسار به ­ازای عیارحدهای مختلف، دارای تغییرات کم هستند و با فاصله اندک از یکدیگر قرار دارند. تغییرات اندک و حالت تقریباً افقی این دو نمودار به دلیل بالا بودن عیار مجموع سرب و روی کانسار در بخش ذخیره قطعی است که با استفاده از فعالیت­های اکتشافی بیشتری اکتشاف شده و خطای تخمین آن نیز کمتر است.
نتیجه­‌گیری
در این پژوهش از رویکردی نوین برای تخمین زمین­آماری عیار و میزان ذخیره کانسار سرب و روی عمارت استفاده شد. برای دست­یابی به هدف، براساس داده­کاوی ترازهای ارتفاعی مختلف کانسار، ابتدا تخمین دوبعدی زمین­آماری با روش لاگ­کریجینگ مقاطع برای هر تراز ارتفاعی به­طور جداگانه صورت گرفت و سپس مدل بلوکی سه­بعدی کانسار ایجاد شد. در واقع در این پژوهش داده­های اکتشافی موجود در سطوح افقی در ترازهای ارتفاعی مختلف کانسار سرب و روی عمارت، از طریق واریوگرافی و تخمین زمین­آماری دوبعدی، به کمک بلوک­بندی و ساخت مدل بلوکی، به فضای سه­بعدی تعمیم داده شد. نتایج پژوهش حاضر نشان می­دهند که در برخی از موارد برحسب شرایط کانسار مورد مطالعه و نوع فعالیت­های اکتشافی انجام­شده، می­توان عملیات مدلسازی و تخمین ذخیره کانسار را با استفاده از روش­های نوین و دقیق زمین­آماری به­طور ساده­تر و البته با دقت مطلوب انجام داد. این فرآیند، با رویکرد کاهش بعد فضای تخمین و سپس تبدیل فضای تخمین از دو به سه­بعد انجام می­گیرد. نتایج این پژوهش برای تمام کاربران علوم زمین شامل زمین­شناسان، مهندسین اکتشاف و استخراج معدن مفید خواهد بود.
 

کلیدواژه‌ها

موضوعات


عنوان مقاله [English]

Estimation of cross-sectional kriging for modeling and evaluation of ore reserve (Case study: Emarat Pb-Zn deposit)

نویسندگان [English]

  • Reza Ahmadi
  • Mohammad Saleh Ahmadi
Department of Mining Engineering, College of Earth Sciences Engineering, Arak University of Technology, Arak, Iran
چکیده [English]

Introduction
Modeling the spatial distribution of grade and reserve estimation are the most important issues and the main goal of exploration operation. This kind of modeling depending on the amount, type and method of carried out exploration works and available exploration information, is performed using various methods. In the following situations, underground exploratory works such as exploration tunnels are carried out: 1- vertical or steep slope mineral deposits, 2- mountainous geomorphology with high altitudes and rough topography, 3- difficult or impossible situation of drilling, and 4- very high thickness or volume of overburden on the mineral deposit. In fact, in such conditions, for the most of mineral deposits there is continuity in the vertical direction whereas the changes are mostly along the horizontal surfaces. In the present research, block modeling and geostatistical reserve estimation of the Emarat Pb-Zn deposit located in the Markazi province have been carried out using log-kriging method with a different approach and a delicate technique. To achieve the goal, firstly, processing, modeling and geostatistical estimation of total Pb-Zn assay data obtained from exploratory tunnels drilled at the various elevation levels were performed as two-dimensional, separately. Then, using the gained results, three-dimensional estimation process for whole of the deposit was done. Such an approach for geostatistical estimation of a deposit grade and reserve discovered by underground exploratory works, has not been reported in any research before.
 
Materials and Methods
Study deposit and available exploration data
The Emarat Pb-Zn deposit is located about 45 km southwest of Arak city, between longitudes 49°30' to 49°45' East and latitudes 33°45' to 34°0' North, in an area with the altitude of 2,180 m above sea level. The topography of the Emarat region is very uneven. Uniform stratigraphy, severe folding, absence of igneous rocks and strati-banding are the geological characteristics of the region, with the strike of folding conforming to the trend of the Zagros folding. The Emarat Pb-Zn deposit is located on the Sanandaj-Sirjan tectonic zone and the Malayer-Isfahan lead and zinc metallogenic belt. The geological units of the region are generally carbonate, shale, and marl sedimentary units from the Jurassic to Cretaceous periods, with folding and faulting, which are composed of steep slopes. In general, the Emarat Pb-Zn deposit is a synclinorium trending northwest-southeast, with a length of 1.5 km and a width ranging from 250 to 850 m. In the Emarat deposit, lead and zinc minerals are located within a calcareous siliceous layer interface of a dark gray massive limestone band in the footwall and a Cretaceous shale layer in the hanging-wall. This formation belongs to the lower or middle Cretaceous.
In the Emarat Pb-Zn deposit, many activities containing drilling of horizontal exploratory- exploitational tunnels and crosscuts with a total length of 11,000 m in six elevation levels with a height difference of about 10 m have been carried out, whereas 1,238 Pb-Zn assay data of samples taken from the tunnels are available.
 
Cross-sectional kriging method
From the dimensional viewpoint, in the case of two-dimensional acquisitioned data, such as assay data of surface samples, processing operation and geostatistical estimation with Kriging method is done as two-dimensional. There are also cases in which, although the data acquisition is three-dimensional, it is better to perform the geostatistical estimation first on two-dimensional surfaces, then the estimation results obtained from different levels are combined with each other and the final three-dimensional model is produced. This type of kriging estimation is called cross-sectional kriging, which can be done along a series of horizontal or vertical sections. Of course, its use in horizontal sections is more logical and desirable. Since in such cases there is practically no sample and consequently no data in the gap between two consecutive adjacent surfaces, and all processed and estimated data is related to samples taken from two-dimensional surfaces, therefore, the distance and spatial relationship of the data at each level is higher. So that the continuity of the deposit is better revealed at each level.
Log-Kriging method
If the statistical distribution of the used data is not normal, linear kriging estimation methods cannot be employed, because in this case, there will be a correlation between variance and mean and pseudo-anisotropy appears in the strike variograms. In such conditions, it is better to normalize the data with a suitable transformation method such as logarithmic so that linear methods can be used for estimation. Next, the estimation operation is applied on the logarithm of data with the ordinary (or simple) kriging method, and then the estimated values are converted to real values with an inverse transformation. This non-linear kriging method is called log-kriging.
 
Results and Discussion
Based on the statistical processing of total Pb-Zn assay data for the elevation levels of 2032, 2024, 1998, 1988, 1978 and 1964-1968 m separately, data distribution was known lognormal, normalized with two and three parameters logarithmic transformation. Also, plotting standard deviation against the total Pb-Zn assay data of the various elevation levels revealed dependence of the variance with mean. To analyze the spatial structure of the region, horizontal strike variograms with azimuths ranging from 25 to 135 degrees were drawn for each elevation level separately, as well as two better horizontal variograms in perpendicular directions were selected. Due to the lognormal distribution of the assay data, to avoid the appearance of pseudo-anisotropy, the variography was performed for the transformed total lead and zinc assay data. All theoretical variogram models fitted over empirical variograms are spherical type, and strike variograms in various directions have the same sill but different ranges. Therefore, the studied area has geometric anisotropy. Using variography of the different elevation levels and determination of ellipse search radii, grade of each elevation level was estimated with geostatistical cross-sectional log-kriging method for 10 by 10 m blocks as well as isograde map from estimation process was drawn. According to the isograde maps at most elevation levels of the Emarat Pb-Zn deposit, the mineral deposit grade in the eastern half is higher than the western part. Afterward, the estimation process was validated through cross-validation approach with the well-known jackknife method and kriging estimation error maps of each elevation level. In such validation method based on the variogram model, each time one of the input data (known) is estimated through the Kriging method using neighboring samples around that sample, while the estimated values are compared with the actual ones. In other words, any known data is estimated by assuming that its value is unknown. In this regard, the determination coefficient of regression between actual and estimated grade values in all cases except elevation level of 2024 is more than 0.5, indicating good correlation between the data. Thus, the estimation has a good validity.
To create block model of the deposit, a prototype grid model with size of 1110×530×265 m, cell size of 10×10×10 m and estimated assay data of different elevation levels through the advanced inverse distance weighted (AIDW) algorithm, was used. The size of blocks was set to 10*10*10 m based on the size and extension area of mineral deposit in each elevation level, as well as distance of the elevation levels. In AIDW algorithm, it is possible to weight the distance with different powers in various directions. In this case, weights of two and one were assigned to the data in the horizontal (exploration tunnels at each elevation level) and vertical directions (distance between different elevation levels) respectively, due to greater variety in the horizontal direction. At the end, different categories of proven, probable and prospected reserves were determined for seven cut-off grades of 4, 4.5, 5, 5.5, 6, 6.5 and 7 percent.
According to the tonnage-grade diagram of the total deposit for different categories of proven, probable and prospected reserves, the proven reserve and related average grade diagrams reveal little changes for various cut off grades locating a short distance from each other. The low changes and almost horizontal state of both two graphs are due to the high total grade of lead and zinc in the proven reserve category, discovered using more exploration activities with less estimation error.
 
Conclusion
In this research, a new approach was used for geostatistical estimation of grade and ore reserve, in the Emarat Pb-Zn deposit. To achieve the goal, based on data mining of different elevation levels of the deposit, first a 2D geostatistical estimation was carried out using the cross-sectional log-kriging method for each elevation level, separately. Afterward, 3D block model of the deposit was created. In fact, in this study, the exploration data available on horizontal surfaces at the different elevation levels of the Emarat Pb-Zn deposit were generalized to 3D space through variography and 2D geostatistical estimation, by creating block model. The results of the research show that in some cases according to the conditions of the studied deposit and type of carried out exploration activities, modeling and estimation of ore reserve is possible more simply with optimal accuracy through the relatively new and accurate geostatistical methods. This scenario was performed by the approach of dimension reduction of estimation space and then converting the estimation space from two to three dimensions. The results of this research will be useful for all earth sciences users, including geologists, mining exploration and exploitation engineers.
 

کلیدواژه‌ها [English]

  • Cross-sectional log-kriging
  • Geostatistical estimation
  • Elevation level
  • Tunnel
  • Emarat Pb-Zn deposit
Ahmadi, R. and Ehsan-nejad, J., 2022. Employing non-linear geostatistical estimation methods for grade modeling and ore reserve estimation of Yazd, Aliabad copper deposit based on selecting optimal size of blocks, Earth and Statistics, v. 2(1), p. 1-12.
Ahmadi, R. and Sadat Koodehi, S.M., 2018. Classification and reserve estimation of Robat Arregije Pb-Zn deposit, Khomein Township, Markazi Province, using geostatistical methods. New Findings in Applied Geology, v. 12(24), p. 39-53 (In Persian). 
Ahmadi, R., 2010. Application of statistical patterns for ore reserve estimation emphasis to Ali-abad, Yazd copper mine. Arak University of Technology, Arak, Report 1, 102 p (In Persian). 
Ahmadi, R., 2011. Comparison of the results of linear and non-linear geostatistical methods for modeling and evaluation of Saveh North-Narbaghi copper ore reserve. Quarterly Iranian Journal of Geology, v. 14(56), p. 43-59 (In Persian). 
Ataeepour, M., 2019. Principles of 2D ore-body modelling. Amirkabir University of Technology (Tehran Poly technique), Tehran, 326 p (In Persian). 
Dimitrakopoulos, R., 2020. Ore reserve estimation and strategic mine planning: stochastic models and optimizations with case studies, Springer-Verlag New York Inc., 325 p.
Ehya, F., Lotfi, M. and Rasa, I., 2010. Emarat carbonate-hosted Zn–Pb deposit, Markazi Province, Iran: A geological, mineralogical and isotopic (S, Pb) study, Journal of Asian Earth Sciences, v. 37, p. 186-194.
Erickson, A.J., 1992. Geological interpretation, modeling and representation. In: H. Hartman (Editor), SME Mining Engineering Handbook. SME-AIME, New York, p. 333-343.
Faraji, K., 2009. Exploitation plan of Emarat lead and zinc deposit: A report by Shahin. Industrial and Mining Company, Iran, 30 p.
Hassani-Pak, A.A. and Sharafodin, M., 2001. Exploration data analysis. Tehran University Press, 987 p (In Persian). 
Hassani-pak, A.A., 1998. Geostatistics. Tehran University Press, 314 p. (In Persian).
Jones, O., Aspandiar, M.F., Dugdale, A. and Smith, B., 2019. The business of mining: mineral deposits, exploration and ore-reserve estimation. London, CRC press, 194 p.
Jun, S. and Dongsheng, T., 2012. The Jackknife and bootstrap (Springer series in statistics), Springer, 534 p.
Karimpour, M.H. and Saadat, S., 2004. Applied Economic Geology. Mashhad University Press, Mashhad (In Persian). 
Madani, H., 1995. Basics of Geostatistics. Amirkabir University of Technology- Tafresh branch, Tafresh, 659 p (In Persian).
Melakpour, H., 2009. Report of exploration operation of the Emarat mine, Shahin Industrial and Mining Company (In Persian). 
Rajabi, A., Rastad, E. and Canet, C., 2012. Metallogeny of Cretaceous carbonate-hosted Zn–Pb deposits of Iran: geotectonic setting and data integration for future mineral exploration, International Geology Review, v. 54(14), p. 1649-1672.
Rastad, E., 1981. Geological, mineralogical, and ore facies investigations on the Lower Cretaceous stratabound Zn-Pb (Ba-Cu-) deposits of the Irankuh Mountain range, Esfahan, West Central Iran, Ruprecht-Karls-Universität Heidelberg, 334 p.
Rendu, J.M., 1981. An introduction to geostatistical methods of mineral evaluation, South African Institute of Mining and Metallurgy monograph series, Johannesburg, 84 p.
Rossi, M.E. and Deutsch, C.V., 2014. Mineral resource estimation, 1st Edition, Springer Dordrecht Heidelbrg New York London, 332 p.
Sadat Koodehi, S.M., 2017. Reserve estimation of Khomein-Robat Pb-Zn deposit using geometrical and geostatistical methods. MSc. thesis of mining engineering, Arak University of Technology, 128 p (In Persian). 
Yates, S.R. and Warrick, A.W., 1987. Estimating soil water content using Co-Kriging, Soil Science Society of America Journal, v. 51, p. 23-30.