استفاده از رویکرد‌های بنیادی زمین آمار به منظور شناسایی نواحی محتمل تغذیه‌کننده آبخوان دشت تهران-کرج

فصلنامه انسان و محیط زیست، شماره 47، زمستان 97

استفاده از رویکرد­های بنیادی زمین آمار به منظور شناسایی نواحی محتمل تغذیه­کننده آبخوان دشت تهران-کرج

بهزاد رایگانی^{[1] *}

behzad.rayegani@gmail.com

چکیده

آب زیرزمینی به خصوص در مناطق خشک و نیمه­خشک یکی از مهم­ترین منابع تامین­کننده آب شرب به حساب می­آید. برخلاف تصور عمومی در این مناطق، منابع آبی بیش­تر از نظر کیفیت با مشکل روبرو است تا کمیت، بنابراین پایش روند کیفیت آب و شناسایی منابع آلوده کننده آن، همواره به عنوان یکی از دغدغه­های اصلی پژوهش­گران این نواحی به شمار می­رود. در این پژوهش به منظور شناسایی نواحی آلوده­کننده آب زیرزمینی، اقدام به شناسایی نواحی محتمل تغذیه سفره آب آبخوان دشت تهران-کرج شده است. بدین منظور پس از استخراج آمار مربوط به سطح تراز آب زیرزمینی، از روش­های مختلف درون­یابی و زمین­آمار برای تولید تصویر سطح تراز استفاده شده است. دو نرم­افزار Gstat و تحلیل­گر زمین­آمار برای انجام این مطالعات بکار گرفته شده­اند، سپس عملکرد و توان هر کدام در تولید تصاویر سطح مورد ارزیابی قرار گرفته است. با استفاده از روش­های مختلف آماری، عملکرد مدل­های درون­یابی  مورد سنجش قرار گرفت. بر اساس نتایج این پژوهش نرم­افزار تحلیل­گر زمین­آمار انعطاف­ بهتری برای انجام تحلیل­های ویژه بر روی داده­ها نشان داد. بر اساس روش داده­های تعلیمی و آزمون، روش­های درون‌یابی بسیار به یکدیگر شبیه نشان دادند، اما خروجی تصاویر سطحی از منظور شبیه­سازی حرکت آب زیرزمینی و هیستوگرام تصاویر با یکدیگر بسیار متفاوت به نظر می­رسید. در نهایت بر اساس مقایسه خروجی مدل‌ها با واقعیت زمینی، مناسب­ترین روش، کریجینگ کلی نشان داد. بر همین اساس تصاویر سطح تراز در مقاطع مختلف زمانی بوسیله روش مذکور تولید و در نهایت با شبیه­سازی جهت جریان آب زیرزمینی محتمل­ترین نواحی تغذیه آبخوان شناسایی گردید.

کلمات کلیدی: زمین­آمار، درون‌یابی ، نرم­افزار Gstat، ماژول تحلیل­گر زمین­آمار، آب­زیرزمینی، جهت جریان آب زیرزمینی، ناحیه تغذیه­کننده آبخوان.

Applying fundamental geostatistical approaches to determine potential zones of groundwater recharge of Tehran-Karaj Plain Aquifer

Behzad Rayegani ^[2]*(Corresponding Author)

behzad.rayegani@gmail.com

Abstract

Groundwater, especially in arid and semi-arid is one of the main sources of drinking water. Contrary to public perception in these areas the water quality is more important than quantify. Therefore, monitoring of water quality and identification of pollution sources is one of the main concerns of researchers in this area. In this study, in order to identify sources of groundwater contamination, areas are likely aquifer recharge sources of Tehran-Karaj plain detected. So, after extraction of groundwater level data different interpolation and geostatistics methods are used to create surface images. The Gstat software and Geostatistical Analyst were used for this study then performance and ability of each one to produce surface images are evaluated. According to the results, the Geostatistical Analyst software has better flexibility to do the special analysis. Based on the training and test data interpolation methods were very similar, but the surface images regarding groundwater direction flow and histogram look very different. According to the results, universal kriging showed better performance. Accordingly, surface images of different time created by an appropriate method to simulate groundwater flow direction and eventually areas were identified which are more likely ground water recharge source of Tehran-Karaj aquifer.

Keywords: Geostatistics, Interpolation, Gstat Software, Geostatistical Analyst, Groundwater, Groundwater Direction Flow, Aquifer Recharge Area.

مقدمه

آب از منابع با ارزش طبیعی است که برای حیات بشر و سلامت زیست­بوم حیاتی می­باشد. از میان ذخایر آبی، آب زیرزمینی به‌خصوص در مناطق خشک و نیمه­خشک یکی از مهم‌ترین منابع تامین­کننده آب شرب است(1). بر خلاف تصور عمومی در این مناطق، منابع آبی بیش‌ترین از نظر کیفیت با مشکل روبرو می­باشند تا کمیت(2). بنابراین پایش روند کیفیت آب (3) و شناسایی منابع آلوده­کننده آن (1، 4) همواره به عنوان یکی از دغدغه­های اصلی پژوهش‌گران این نواحی به شمار می­رود. یکی از روشهای شناسایی منابع سطحی آلوده­کننده آب زیرزمینی، بررسی جهت جریان آب به منظور آشکارسازی واگرایی در آن است. به طور کلی وجود جریان‌های واگرا در منابع آب زیرزمینی یکی از نشانه­های تغذیه سفره آب می­باشد(5). اما یکی از مهم‌ترین چالش­ها در استفاده از این فرآیند، نیازمندی آن به یک تصویر پیوسته سطح تراز آب زیرزمینی است. برای تهیه این گونه تصاویر سطح تراز، از عمق و ارتفاع چاه­های مشاهداتی استفاده می­شود. به دلیل هزینه بسیار زیاد، معمولاً چاه­های مشاهداتی در سطح آبخوان از یکدیگر فاصله زیادی دارند و تهیه تصویر سطح تراز دقیق از آن‌ها یک چالش اساسی است. یکی از روش‌های حل این مشکل، استفاده از ابزارهای سامانه­های اطلاعات جغرافیایی می­باشد.

تکنیک­های درون‌یابی  مکانی[3] در سامانه­های اطلاعات جغرافیایی، ابزارهای بسیار قدرتمندی برای پیش­بینی مقادیر ارزشی تصاویر سطحی[4] به حساب می­آیند(6). در تعریف، درون‌یابی  مکانی فرآیندی است در آن ارزش ویژگی یک مکان خاص بوسیله ارزش‌ها مشخص برخی سایت‌های دیگر پیش­بینی می­شود(7). تمامی روش‌های درون‌یابی  از قانون توبلر[5] تبعیت می­کند، مطابق با این قانون تمامی مکان‌ها باهم در ارتباطتند اما مکان‌های نزدیک به هم ارتباط بیش‌تری نسبت به مکان‌های دور از هم دارند(7-9). بنابراین درفرآیند درون­یابی ارزش‌ها نامشخص، بر اساس ارزش‌ها مکان‌های مجاور یا نزدیک تعیین می­گردد. در تمامی روش‌های درون‌یابی  معنی " نزدیکی" در قانون توبلر تعریف می­شود. معمولاً روش‌های درون‌یابی  مکانی بر اساس فرضیه­هایشان به دو دسته کلی به چند صورت مختلف تقسیم می­شوند(7): کلی[6] یا محلی[7]: در روش محلی از میان داده­ها فقط یک بخش، آن هم به صورت فضای همسایگی برای محاسبه یک مجهول بکار می­رود؛ دقیق[8] یا تقریبی[9]: در روش دقیق، تصویر سطح ارزشی دقیقاً از روی ارزش‌ها معلوم عبور می­کند، در حالی‌که در روش‌های تقریبی، ارزش‌ها مکان‌های معلوم نیز دوباره تقریب زده می­شود؛ روش‌های قطعی (جبری[10]) یا روش‌های تصادفی[11] : در روش‌های قطعی فرض می­شود رفتار متغیر کاملاً مشخص است بنابراین خطاهای ممکن ارایه نمی­شود، در حالی‌که در روش‌های تصادفی به دلیل فرض بر تصادفی بودن متغیر، امکان ارایه برآوردهای احتمالی (احتمالات) وجود دارد. معمولاً در جداسازی روش‌های درون‌یابی  مکانی از شیوه تقسیم­بندی قطعی یا تصادفی استفاده می­شود(8،9). از جمله روش‌های قطعی می­توان به روش وزن­دهی بر اساس معکوس فاصله به اختصار[12]IDW که جزء روش‌های دقیق و محلی است؛ روش درون‌یابی  چندجمله­ای کلی یا به اختصار GPI[13] که از نوع روش‌های تقریبی و کلی است؛ درون‌یابی  بوسیله تابع شعاعی محور یا به اختصار RBF[14] که از روش‌های دقیق است(هم محلی و هم کلی را شامل می­شود)؛ و روش درون‌یابی  چندجمله­ای محلی یا LPI[15] که یک روش تقریبی و محلی است نام برد(8-11). معروفترین روش درون‌یابی  تصادفی، زمین­­آمار[16] می­باشد(12). استفاده از روش‌های زمین­آمار به خصوص روش‌های کریجینگ[17] که از جمله روش‌های محلی، تقریبی و تصادفی به شمار می­رود، در مطالعات علوم مختلف، سابقه چند ده ساله دارد(7، 12). دانیل کریج[18] پدر زمین­­آمار در معادن آفریقای جنوبی از یک روش آماری تجربی برای پیش­بینی درجه عیار سنگ از نمونه­های برگرفته از مناطق مجاور استفاده نمود که بعدها بوسیله مترون[19] حالت فرمولی به آن داده شد و به احترام کریج، نام این روش، کریجینگ لقب گرفت(13). در روش کریجینگ همانند IDW به منظور پیش­بینی ارزش مناطق اندازه­گیری نشده از وزن­دهی به ارزش مناطق مجاور استفاده می­شود، ولی برخلاف IDW وزن‌ها تنها بر اساس نزدیکی و فاصله نیستند، بلکه وزن‌ها از یک پلات فاصله-ارزشی خاص به نام شبه­واریوگرام[20] ناشی می­شوند(9) که در آن جهت نیز در نظر گرفته می­شود. شبه­واریوگرام از دو محور تشکیل شده است که در محور x آن فاصله جفت نمونه­ها قرار دارد و در محور y مقدار شبه­واریانس[21] یا متوسط گشتاور درجه دوم ارزش همان جفت نمونه را نشان می­دهد(شکل 1). پس از ترسیم شبه واریوگرام، داده­ها به یک­سری وقفه­های فاصله­ای[22] تقسیم می­شوند، متوسط شبه­واریانس برای هر وقفه بدست می­آید(بوسیله معادله 1) و در نهایت یک
معادله برای تعیین وزن‌ها به نقاط نماینده وقفه­های فاصله­ای و شبه واریانس، برازش می­یابد(شکل 1). به شکل نهایی منحنی بدست آمده واریوگرام گفته می­شود(شکل 1) که از آن برای تعیین وزن هر نمونه استفاده خواهد شد و دارای سه جزء اصلی است(شکل 1): Nugget، Sill و Range. دقت خروجی حاصل از روش درون‌یابی  کریجینگ به شدت به مدل برازش داده شده و مقدار عددی این سه جزء واریوگرام وابسته است(12،13). بنابراین یکی از مهم‌ترین چالشها در استفاده از روش زمین­آمار کریجینگ، تعریف صحیح مدل واریوگرام می­باشد.

به دلیل در نظر گرفتن جهت و فاصله در وزن­دهی، کاربرد زمین­آمار در مطالعه کمیت و کیفیت آب زیرزمینی کاملاً شناخته شده است(1-3، 11، 12، 14-19). با این وجود هنوز هم مدل واریوگرام و تعریف صحیح آن یکی از عمده­ترین دغدغه­های این پژوهش­ها باقی‌مانده است(12، 13). یکی از روش‌های حل این مساله، بکارگیری درست از قابلیت‌های نرم­افزارهای زمین­آمار موجود می­باشد.

معادله (1) شبه واریانس

شکل 1- شیوه نمایش جفت نمونه­ها در شبه­واریوگرام(محور x و y) و بخش­های مختلف یک واریوگرام(20)

نرم افزارهای زیادی برای انجام درون‌یابی ­های مکانی داده­های جفرافیایی به خصوص روش‌های زمین­آمار وجود دارد(9، 13، 21-24) و هر یک از این نرم­افزارها نقاط قوت و ضعف متفاوتی دارند(21). این وجود اساس استفاده از این نرم­افزارها تقریباً مشابه است و نکته مهم در کسب نتایج مطلوب از این نرم­افزارها، تسلط کافی داشتن بر مباحث نظری روش‌های زمین­آمار می­باشد(13).

از آنجا که هدف اصلی این پژوهش تهیه نقشه­های جهت جریان آب زیرزمینی به منظور شناسایی نواحی محتمل تغذیه سفره آب زیرزمینی است و همان‌گونه که اشاره شد، شیوه درونیایی مکانی به شدت تصویر سطح نهایی بدست آمده را متاثر می­سازد، بنابراین در این مطالعه به دقت بررسی خواهد شد، در مطالعات مشابه چه گام­هایی توصیه شده است و با در نظر داشتن این مراحل نتایج نهایی چگونه متاثر خواهند شد. در این بررسی از قابلیت­های دو نرم­افزار و دستورالعمل­های توصیه شده توسط آن‌ها استفاده شده است:

1. نرم افزار Gstat : یک نرم­افزار منبع باز که بر اساس بسیاری از مطالعات، توانایی آن در روش‌های مختلف کریجینگ به اثبات رسیده است(13، 21-24) و در حال حاضر نسخه گرافیکی آن در نرم افزار IDRISI وجود دارد(25، 26).
2. مدول Geostatistical Analyst در نرم­افزار ArcGIS که در بسیاری از مطالعات بکار گرفته شده است و از جمله قوی­ترین نرم­افزارهای زمین­آمار به شمار می­رود(9، 13، 21).

در این بررسی ضمن ارایه نتایج حاصل از این دو نرم­افزار، شیوه بکارگیری صحیح هر یک بر اساس موارد توصیه شده در منابع علمی نیز مورد اشاره قرار خواهد گرفت.

 مواد و روش­ها

 منطقه مورد مطالعه

آبخوان دشت تهران-کرج با مساحتی در حدود 2330 کیلومترمربع بین طول­های شرقی ̋ 627/46 ̒  44 ̊  50 تا ̋ 67/49 ̒  37 ̊  51 و عرض­های شمالی ̋ 545/15 ̒  17 ̊  35 تا  ̋ 8/05 ̒  55 ̊  35 در میان استانهای تهران و البرز قرار گرفته است(شکل 2). با توجه به قرارگیری دو کلان شهر تهران و کرج در محدوده آن، این سفره آب زیرزمینی یکی از مهم‌ترین آبخوان­های موجود در کشور می­باشد. ارتفاع محدوده­ی آبخوان بین 872 تا 1776 متر متغیر است و میزان ارتفاع از شمال آبخوان به سمت جنوب آن کاهش می­یابد. بدون در نظر گرفتن محدوده و حریم­ شهرها و روستاها و تنها بر اساس نوع استفاده از زمین، کاربری بیش از 51 درصد از وسعت آبخوان به کشاورزی و باغبانی اختصاص یافته است، بیش از 32 درصد از مساحت آن را مناطق مسکونی اشغال نموده­اند، در حدود 5/11 درصد از سطح آن را زمین­های رها شده و فاقد استفاده تشکیل می­دهد و در حدود 4 درصد نیز سهم مراتع و زمین­های جنگل­کاری می­باشد. در مجموع آبخوان تهران-کرج به طور تقریبی 9139646 نفر جمعیت شهرنشین را تحت پوشش قرار داده است (شکل 2).

شکل 2-  آبخوان دشت تهران-کرج و پراکنش شهرهای واقع در آن

 روش مطالعه

 استخراج اطلاعات تراز آب چاه­های مشاهده­ای

در درون آبخوان دشت تهران-کرج در حدود 168 چاه مشاهده­ای وجود دارد که برخی از آن‌ها به دلیل خشک شدن، اندازه­گیری نمی­شوند. آمار عمق این چاه­های مشاهده­ای در سال­های 1390 و 1389 از شرکت آب منطقه استان تهران دریافت شد. به منظور شناسایی داده­های پرت، توزیع فراوانی آمار هر یک این چاه­ها در طول­ سال ترسیم شد تا داده­های غلط شناسایی و حذف شوند. به منظور تهیه آمار تراز آب چاه­ها، یک­بار از میانگین اعداد استفاده شد تا اثر خطاهای ناخواسته از بین برود و یک­بار از آمار ماهیانه چاه­ها استفاده گردید. در نهایت ارتفاع هر چاه از مقدار عددی عمق آن کسر شد تا تراز چاه بدست آید. با توجه به تعدد ماه­ها و بررسی در دو سال مختلف، در ادامه فقط روش‌های بکار گرفته شده بر روی آمار میانگین و نتایج حاصل از آن اشاره خواهند شد. لازم به ذکر است به منظور شناسایی دقیق مناطق محتمل تغذیه سفره، مشابه با همین فرآیند­ها بر روی آمار ماهیانه نیز بکار گرفته شده است.

درون‌یابی  مکانی تراز چاه­ها به روش کریجینگ به کمک نرم­افزار Gstat

در فرآیندهای درون‌یابی  تراز چاه­ها بوسیله نرم­افزار Gstat، دقیقاً گام­های توصیه شده توسط راهنمای این نرم­افزار(22-27) و منابع اطلاعاتی مرتبط با زمین­آمار(7، 9، 13) طی شده است که در ادامه به طور خلاصه به آن‌ها خواهیم پرداخت.

بررسی شبه­واریوگرام و انتخاب فاصله و تعداد وقفه­های فاصله­ای بوسیله Spatial Dependence Modeler

در این ماژول در مرحله اول برای تمامی داده­های وکتوری تراز، مدل آزیموتی تصویر سطح بر اساس روش شبه واریوگرام ساخته شد تا بوسیله آن روند تغییرات در جهت­های مختلف مشخص گردد. بر اساس تصویر سطح مقدار فاصله و تعداد وقفه­ها به طور دستی تغییر نمود تا بیضی[23] آزیموتی تغییرات به وضوح مشخص گردد. بر این اساس در صورتی که در جهت­های مختلف، تغییر میزان وابستگی مکانی^[24] یا تغییرات آزیموتی ارزش پیکسل­های تصویر سطح، به شیوه یکسان نبود و الگوی خاصی در جهت­های مختلف دیده می­شد، از نظر الگوهای پیوستگی مکانی^[25] مدل ناهمسان^[26] بود. در این موارد برای محور طولی و عرضی بیضی آزیموتی تصویر سطح، واریوگرام ترسیم و هر یک از آن‌ها به صورت فایل جداگانه ذخیره می­گردید. در صورتی‌که در واریوگرام­های طولی و عرضی بیضی اختلاف اساسی در میزان Sill دیده شود، نشان می­دهد در داده­ها یک روند کلی[27] وجود دارد که باید عمل بی­روندسازی[28] بر روی داده­ها صورت پذیرد، بنابراین مدل، ناهمسان پهنه­ای^[29] دیده می­شود که برای مرحله دوم باید بر اساس آن عمل نمود(22، 25-27). در صورتی‌که میزان Sill در دو جهت عمود بر هم، تغییر زیادی نکند و فقط میزان Range تغییر نماید، مدل، ناهمسان دیده خواهد شد و نیازی به حذف روند در آن وجود ندارد، اما در مرحله دوم باید از تابع ناهمسان استفاده شود. به هر حال در هر دو مدل، نیاز است واریوگرام این دو جهت جداگانه ذخیره شود.

نکته بعدی در این ماژول تنظیم فاصله و تعداد وقفه­های فاصله­ای می­باشد. مطابق موارد پیشنهاد شده(13، 22، 25-27) با بررسی آماری به گونه­ای عمل شد که تعداد جفت­ها در وقفه­های ابتدایی که بسیار مهم و تاثیر گذار می­باشند کم‌تر از 30 نباشد، همچنین فاصله نهایی حاصل از ضرب تعداد وقفه در فاصله هر وقفه با توجه به فاصله چاه­های پیزومتری بیش‌تر از 30 کیلومتر نشود. هم‌زمان دقت شد بیضی شکل منظم خود را حفظ کند و در واریوگرامها نهایی حدوداً محل Range مشخص باشد.  همچنین برای اطمینان از دقت واریوگرام، بوسیله h-scatterplot در هر مرحله تنظیم فاصله و تعداد وقفه، بررسی می­گردید تا به‌خصوص در چند وقفه اول حالت خطی 45 درجه بیش‌تر وجود داشته باشد.

پس از تنظیم نهایی واریوگرامها ذخیره شدند تا در مرحله بعدی در برازش مدل از آن‌ها استفاده شود.

 برازش مدل نهایی بوسیله Model Fitting

بدون استثناء در تمامی ماه­ها و متوسط سالیانه مدل ناهمسان پهنه­ای وجود داشت، بنابراین مطابق با موارد توصیه شده بدین صورت عمل شد(22، 25، 26):

1. برای حذف اثر روندکلی و مشخص نمودن مدل ناهمسان پهنه­ای در ابتدا در زاویه بیش‌ترین تغییرات (محور عرضی بیضی مرحله قبل) در ساختار شماره 2[30] یک بیضی بی­نهایت با Range بسیار بزرگ و نرخ ناهمسانی[31] بسیار کم به واریوگرام محور طولی نسبت داده شد.
2. سپس بوسیله تنظیم واریوگرام­ها در ساختار شماره 3 و 4 و تنظیم Nugget در ساختار شماره 1 یک واریوگرام ترکیبی به شکلی به نقاط شبه­واریوگرام برازش داده شد که دقیقاً از چند نقطه اول عبور کند و به نقاط بعدی به صورت کم‌ترین مربعات خطا نزدیک باشد و ساختار Range در انتهای آن قابل مشاهده باشد. لازم به ذکر است این ساختارهای برای زاویه کم‌ترین تغییرات تنظیم گردید.
3. مطابق با راهنمای Gstat بهتر برازش به معنی عبور دقیق از تمامی نقاط یا داشتن کم‌ترین مربعات خطا در خط برازش داده شده نیست(22، 25، 26)، بلکه نیاز است بر اساس ماهیت نقاط و داشتن اطلاعات جانبی بهترین خط انتخاب شود. بنابراین با بررسی پراکنش چاه­های مشاهده­ای و نقاط شبه­واریوگرام فاصله مناسب برای Range انتخاب و بر همین اساس مدل برازش یافت. این مدل­ها برای مرحله بعدی ذخیره شدند.

ساخت تصاویر سطح تراز نهایی بوسیله Kriging and Simulation

در این پژوهش دو نوع روش زمین­آمار مورد بررسی قرار گرفته است. برای تمامی تخمین­ها و شبیه­سازی از کریجینگ معمولی[32] استفاده شده است. همچنین با توجه به ناهمسانی پهنه­ای در تمامی داده­ها روش کریجینگ کلی[33] نیز مورد بررسی قرار گرفته است. نکته مهم و تاثیرگذار در خروجی­ها، تعریف فضای همسایگی محلی[34] یا مشخص نمودن تعداد کمینه، بیشنه یا شعاع انتخاب نمونه­ها در فرآیند شبیه­سازی هر پیکسل می­باشد. با توجه به این­که در نرم­افزار Gstat تاثیر تعریف ناحیه همسایگی به صورت چشمی قابل ارزیابی نیست، برای این تعریف به دقت فاصله نقاط بررسی شد و بر آن اساس تعداد کمینه نمونه­برداری برابر 8 و تعداد بیشنه آن معادل 20 در نظر گرفته شد. همچنین به منظور بررسی دقت نقشه­خروجی بجزء فایل پیش­بینی، فایل واریانس نیز تولید گردید تا مشخص شود کدام نواحی پوشش مناسبی نداشته­اند و با توجه به این­که کریجینگ از جمله مدل‌های درون‌یابی  تقریبی است برای مشخص نمودن دقت شبیه­سازی از گزینه اعتبارسنجی از نوع تقاطع[35] نیز بکار گرفته شده است.

درون‌یابی  مکانی تراز چاه­ها به روش کریجینگ به کمک نرم­افزار Geostatistical Analyst

برای انجام این درون‌یابی ­ها از مدول ساحر زمین­آمار[36] جعبه­ابزار[37] تحلیل­گر زمین­آمار[38] استفاده شده است(9، 13). بر خلاف Gstat در اینجا بیش‌تر فرآیندها به طور خودکار یا نیمه خودکار قابل انجام است و همچنین برای انجام یک شبیه­سازی صحیح­تر می­توان بررسی­های متعددی به انجام رساند و بر همین اساس مدل نهایی را تنظیم نمود. در ادامه به طور خلاصه به این بررسی­ها خواهیم پرداخت.

بررسی توزیع فراوانی و نیازمندی داده­ها به تبدیلات عددی

یکی از شرایط برای آنکه روش‌های درون‌یابی  یک تصویر سطح با کم‌ترین میزان خطا تولید کنند، داشتن توزیع نرمال در داده­ها است(9، 13). بنابراین در این مطالعه از هر دو روش بررسی توزیع داده­ها یعنی هیستوگرام و Normal QQPlot استفاده گردید تا ضمن بررسی نرمال بودن داده­ها هم‌زمان بررسی شود تبدیل داده­ها به شکل لگاریتمی یاBox-Cox  می­تواند به نرمال­تر شدن داده­ها کمک نماید یا خیر. در صورتی‌که تبدیل اثر مثبتی بر نرمال شدن داده­ها نشان می­داد در ابتدا قبل از شبیه­سازی از آن برای تبدیل استفاده می­گردید.

بررسی وجود روند کلی در داده­ها و انتخاب شیوه مناسب بی­روندسازی

در صورتی‌که در داده­ها یک روند کلی دیده شود، نیاز است این روند از داده­ها حذف شود تا شبیه­سازی بر روی باقی‌مانده­هایی صورت گیرد که تغییرات محلی را بهتر نشان می­دهند(9، 13، 22، 26). برای حذف این روند در نرم­افزار Gstat  از روش ناهمسانی پهنه­ای استفاده می­شود(26) ولی در مدول تحلیل­گر زمین­آمار امکان شناسایی روندکلی و حذف آن قبل از شروع شبیه­سازی وجود دارد. بنابراین از مدول تحلیل روند[39] استفاده شد تا روندها شناسایی و در صورت نیاز حذف شون

کنکاش شبه­واریانس­ها به منظور بررسی همبستگی ارزشی-مکانی[40] در نمونه­ها

بر اساس قانون توبلر و فرض زمین­آمار هر چقدر فاصله یک جفت داده بیش‌تر شود، اختلاف ارزشی آن‌ها بیش‌تر می­شود(7-9، 13). شبه­واریوگرام اجازه می­دهد این رابطه یا به عبارت بهتر همبستگی مکانی-ارزشی بررسی شود. با بررسی نقاطی که در فاصله کم محور x (مولفه فاصله)، مقدار زیاد شبه­واریانس در محور y (متوسط گشتاور دوم اختلاف ارزشی) دارند یا بر عکس در فاصله زیاد، مقدار شبه واریانس بسیار کم است، می­توان جهت­هایی را شناسایی نمود که الگوی تغییرات ارزشی در آن‌ها به دلایل مختلف محیطی متفاوت می­باشد. بنابراین در اینجا با این بررسی، شناسایی مدل همسان یا ناهمسان شبیه­سازی امکان پذیر می­باشد که در مراحل بعدی فرآیند زمین­آمار از آن استفاده خواهد شد.

 بررسی مدل آزیموتی تصویر سطح به منظور استفاده از مدل ناهمسانی

همانند نرم­افزار Gstat در اینجا نیز امکان تشکیل مدل آزیموتی تصویر سطح وجود دارد، بنابراین به غیر از بررسی شبه­واریانس که در فاز قبل مورد آزمون قرار گرفت، روش تشکیل تصویر آزیموتی سطح به انتخاب نوع مدل، جهت مناسب و تعریف تعداد وقفه­ها و فاصله­ها کمک می­کند که در این پژوهش از این تکنیک نیز استفاده شده است.

 تعریف چشمی فضای همسایگی و بررسی اثر تغییر در پارامترهای آن

یکی از مهم‌ترین تفاوت­های تحلیل­گر زمین­آمار با نرم افزار Gstat، امکان تعریف فضای همسایگی محاسبه محلی به صورت چشمی و بررسی هم‌زمان تغییر اندازه آن در مدل خروجی می­باشد. بنابراین با دقت بهتری امکان تعریف این فضا وجود دارد. در این پژوهش از نوع بخش بندی چهار جهته، استفاده شد و در هر جهت کمینه معادل 2 و بیشینه معادل 5 در نظر گرفته شده است(تا تعداد نقاط با نرم افزار Gstat یکسان باشد).

 بررسی خطاهای تک­تک نقاط و خطای کل شبیه­سازی

در گام آخر پس از تنظیمات در ساحر زمین­آمار و تولید تصویر سطح نهایی، آمار خطای هر نوع و خطای کل در پنجره اعتبارسنجی تقاطع ذخیره شد تا در فرآیند ارزیابی صحت مورد استفاده قرار گیرد.

 دورنیابی مکانی تراز چاه­ها به کمک روش‌های قطعی

با توجه به هدف این پژوهش و اهمیت جهت جریان آب زیرزمینی در شناسایی درست نواحی محتمل تغذیه، بجز روش‌های زمین­آمار از روش‌های قطعی درون‌یابی  مرسوم شامل IDW، RBF و LPI نیز استفاده شد تا بر اساس بررسی­های عملکرد و اعتبارسنجی روش‌های درون‌یابی  مناسب­ترین آن‌ها جهت مدل‌سازی جهت جریان آب مورد استفاده قرار گیرد. بدین منظور از ماژول ساحر زمین­آمار در نرم­افزار ArcGIS استفاده شد.

 اعتبارسنجی روش‌های دورنیابی

بهترین روش اعتبارسنجی، مقایسه خروجی مدل‌سازی با اندازه­گیری­های نقاط جدید پیش­بینی­شده در محل می­باشد(9)، اما معمولاً چنین دسترسی­هایی امکان پذیر نیست(نیاز به حفر چاه مشاهده­ای دارد)، بنابراین متداول است داده­های مورد استفاده به دو دسته  تعلیمی[41] و آزمون[42] تقسیم شوند و مدل‌سازی ابتدا بوسیله داده­ها تعلیمی صورت پذیرد، سپس بر روی داده­های آزمون، اعمال شده و خروجی آن با واقعیت سنجیده می­شود(9). بنابراین به منظور بررسی عملکرد هر مدل درون‌یابی ، ابتدا بوسیله ماژول تولید زیرمجموعه[43] در ابزار تحلیل­گر زمین­آمار داده­ها به دو دسته تعلیمی و آزمون تفکیک شدند(به نسبت 70 به 30)، سپس مدل‌سازی­ها و عملیات­های درون‌یابی  بر اساس داده­تعلیمی انجام و با اعمال مدل بر روی داده آزمون، صحت مدل‌سازی مشخص می­گردید. بدین منظور برای بررسی دقت مدل‌سازی بر روی داده­های آزمون، در نرم­افزار ArcGIS از ابزار اعتبارسنجی GA Layer To Points استفاده و در نرم افزار IDRISI، ماژول رگرسیون یا REGRESS بکار گرفته شده است.

لازم به ذکر است در مورد روش‌های درون‌یابی  دقیق، تنها روش بررسی صحت، تفکیک داده­ها به دو دسته تعلیمی و آزمون است زیرا در این روشها ارزش پیش­بینی­ها دقیقاً با مقادیر ورودی یکسان می­باشند و امکان بررسی دقت مدل‌سازی وجود ندارد، اما در مورد روش‌های تقریب می­توان از آمار اعتبارسنجی به شکل تقاطع استفاده نمود.

همچنین بر آنکه از نظر عملی مشخص شود عملکرد هر روش به روش دیگر چقدر شباهت دارد، تصاویر سطح خروجی با یکدیگر از نظر آماری و خروجی مدل‌سازی جهت جریان آب مقایسه شدند و در نهایت با انتقال این خروجی­ها در نرم­افزار گوگل ارث، از نقطه­نظر کارشناسی هم عملکرد هر مدل درون‌یابی  مورد بررسی قرار گرفت.      

تعریف نواحی محتمل تغذیه سفره

بر اساس نتایج اعتبارسنجی برای هر داده مورد بررسی، بهترین روش درون‌یابی  مشخص و بر همان اساس یک مدل‌سازی درون‌یابی  مجدد با تمامی داده­ها انجام ­پذیرفت و تصویر سطح تراز مربوطه بدست ­آمد. در گام بعدی با استفاده از نرم­افزار ArcHydro Groundwater جهت جریان آب زیرزمینی بوسیله تصویر سطح مذکور مشخص و نواحی دارای جریان واگرا علامت زده می­شد. در نهایت بر اساس مقایسه نتایج در دو سال متمادی 1389 و 1390 محتمل­ترین مناطق تغذیه­کننده سفر آب مشخص گردید.

 نتایج

 مدل‌سازی توسط Gstat

با بررسی تصاویر سطح تراز در ماژول مدلساز وابستگی مکانی[44] مشخص گردید، تقریباً برای تمامی داده­ها مدل ناهمسان پهنه­ای باید بکار گرفته شود، زیرا در جهت­های مختلف مقدار Sill بسیار متفاوت بود(شکل 3).

شکل 3- تصویر آزیموتی سطح تراز (چپ) و شبه­واریوگرام (راست)در جهت­های مختلف داده متوسط سالیانه 89

در مورد روش کریجینگ معمولی، به منظور بکارگیری مدل ناهمسان پهنه­ای برای تمامی داده­ها از چهار ساختار برای ایجاد واریوگرام نهایی استفاده شده است. بر این اساس یکی از ساختارها به شکل یک بیضی بی­نهایت در جهت آزیموتی بیش‌ترین تغییرات با مقدار ناهمسانی بسیار ناچیز تعریف شده است(شکل 4 ساختار شماره 4).

در مورد روش کریجینگ کلی که برای مدل ناهمسان پهنه­ای توصیه شده است، واریوگرام بر اساس جهت بیش‌ترین وابستگی مکانی ایجاد و برای حذف روند کل از نمونه­ها، در ابتدا بوسیله ماژول روند، تصویر سطح روند کلی بدست آمد و در مدل نهایی به عنوان لایه رستری دوم بکار گرفته شد(25).

شکل 4- ساختارهای مختلف واریوگرام (شکل راست پایین) و واریوگرام ترکیبی (شکل راست بالا) داده متوسط سالیانه 89

مهم‌ترین مشکل استفاده از نرم­افزار Gstat مشخص نبودن اثر تعریف فضای همسایگی قبل از انجام شبیه­سازی می­باشد، بنابراین این فضا با بررسی چشمی پراکنش نقاط مشخص گردید(شکل 4). همچنین برای درک صحت شبیه­سازی از اعتبارسنجی به شکل تقاطع در مورد کل داده نیز استفاده شده است(شکل 4).

شکل 5- تعریف فضای همسایگی، مدل نهایی و اعتبارسنجی تقاطع در مورد داده متوسط سالیانه 89 و تصویر سطح نهایی

رابطه همبستگی بین پیش­بینی­ها و مشاهدات در مورد کل داده در زمان‌های مختلف بین 9494/0 تا 9497/0 متغیر بود. اما به منظور انتخاب مناسب­ترین شیوه درون‌یابی  تصویر سطح، از روش داده آزمون و تعلیمی استفاده شده است. در این حالت نیز همبستگی بین مشاهدات و شبیه­سازی­ها در روش کریجینگ بین 949/0 تا 953/0 متغیر نشان می­دهد.

مدل‌سازی توسط Geostatistical Analyst

با بررسی توزیع فراوانی داده­های تراز آب مشخص گردید این  توزیع در مورد تمامی تاریخها نرمال نیست و تبدیل داده نیز کمکی به نرمال­تر شدن آن‌ها نمی­کند(شکل 6). بنابراین تبدیل داده­ها صورت نپذیرفت.

شکل 6- بررسی هیستوگرام و  Normal QQPlot به منظور تشخیص نرمال بودن در مورد داده متوسط سالیانه 89

با بررسی روند همان‌گونه که در مورد نرم­افزار Gstat نیز بیان شد، مشخص گردید داده­ها داری روند است(شکل 7) و یک روند قوی خطی از بالا به پایین (محور y) و یک روند ضعیف­تر در محور x در مورد اغلب داده­ها وجود دارد. بنابرین برای حذف روند کلی، روند بالا به پایین که بسیار قوی به نظر می­رسید، در فرآیندهای بعدی مدل‌سازی حذف گردید(شکل 7).

شکل 7- بررسی روند کلی در مورد داده متوسط تراز سالیانه آب سال 89

در بررسی شبه واریوگرام نیز مشخص گردید با افزایش فاصله، متوسط گشتاور درجه دوم ارزشی داده های مربوط به غرب و شرق آبخوان افزایش نمی­یابد(شکل 8) و بیش‌ترین تغییرات مربوط روند بالا به پایین است.

شکل 8- بررسی شبه­واریوگرام در مورد داده متوسط تراز سالیانه آب سال 89

بنابراین با بررسی­های انجام شده بر روی داده مشخص گردید در مورد تمامی آن‌ها نیازی به تبدیل عددی نیست(کمکی به نرمال­تر شدن داده­ها نمی­کند)، روند بالا به پایین خطی باید از داده­ها حذف شود تا شبیه­سازی بر روی باقی‌مانده­ها صورت پذیرد، جهت در داده­ها بسیار مهم است، بنابراین مدل ناهمسان باید در نظر گرفته شود و جهت مدل‌سازی بر اساس بیضی تغییرات آزیموتی سطح و در جهت کم‌ترین تغییرات باید تعریف شود(شکل 9).

شکل 9-تعریف مدل ناهمسان و جهت بیضی در مورد داده متوسط تراز سالیانه آب سال 89

مهم‌ترین ویژگی ساحر زمین­آمار نسبت به Gstat امکان بررسی اثر تغییر فضای همسایگی در خروجی مدل‌سازی می­باشد(شکل 10). بر همین اساس با بررسی چشمی نقاط مختلف آبخوان، مشخص گردید انتخاب نوع بخش بندی به صورت چهارجهته و
داشتن کمینه 2 و بیشینه 5 در هر بخش آن، مناسب این پژوهش می­باشد. با این انتخاب امکان مقایسه خروجی ساحر زمین­آمار با Gstat نیز امکان پذیر می­باشد.

شکل 10- تعریف فضای همسایگی شبیه­سازی محلی در مورد داده متوسط تراز سالیانه آب سال 89

در مورد خروجی تحلیل­گر زمین­آمار رابطه همبستگی بین پیش­بینی­ها و مشاهدات در مورد کل داده در زمان‌های مختلف بین 9593/0 تا 968/0 متغیر بود. در روش داده آزمون و تعلیمی نیز همبستگی بین مشاهدات و شبیه­سازی­ها در روش کریجینگ بین 9494/0 تا 9768/0 متغیر نشان داد.

 مقایسه عملکرد دو نرم­افزار و روش‌های درون‌یابی

به منظور مقایسه عملکرد دو نرم­افزار و مقایسه عملکرد روش کریجینگ با سایر روش‌های درون‌یابی  در این پژوهش چندین مکانیزم مختلف بکار گرفته شده است. در روش اول بیش‌ترین همبستگی مشاهده شده در مورد داده­های آزمون مشاهده شده و پیش­بینی شده معیار قرار گرفته است(جدول 1). همان‌گونه که در جدول (1) مشاهده می­شود عملاً امکان ارزیابی عملکرد و گزینش بهترین روش درون‌یابی  یا نرم­افزار مطمئن­تر با این روش وجود ندارد(در واقع بیش‌ترین تفاوت همبستگی در حدود 02/0 می­باشد). مطابق با ضرایب همبستگی به نظر می­رسد تفاوت معنی­داری بین نرم­افزارها و روشها وجود ندارد. بنابراین با توجه به اهمیت تصویر سطح تراز در این پژوهش به منظور درک بهتر عملکرد روش‌های درون‌یابی  از مکانیزم­های دیگری نیز استفاده شده است.

جدول 1- بیش‌ترین همبستگی مشاهده شده در مورد داده­های آزمون

با توجه به برآورد تصویر سطحی در تمامی روشها، در ابتدا رابطه همبستگی میان تخمین­های هر روش در کل تصویر با روش‌های دیگر بدست آمد. هر چند مطابق با انتظارات بر اساس جدول (1)، اعداد ضرایب همبستگی بین روش­های مختلف، بسیار زیاد مشاهده گردید(جدول 2)، ولی تغییرات ابرپراکنش نقاط حول خط کم‌ترین مربعات بسیار شدید و محسوس نشان داده شده است(شکل 11) و هیستوگرام تصاویر خروجی هر مدل نیز این تغییرات آشکار عملکرد مدل‌ها را به خوبی نشان می­دهد(شکل 12).

جدول 2- مقایسه ضرایب همبستگی میان روشهای مختلف درون‌یابی  و زمین­آمار

شکل11- ابرپراکنش نقاط در رابطه رگرسیونی کریجینگ معمولی با RBF در مورد داده متوسط تراز سالیانه آب سال 89

شکل 12- مقایسه هیستوگرام تصاویر بدست آمده از روشهای مختلف در مورد داده متوسط تراز سالیانه آب سال 89

هر چند در روشهای مختلف درون‌یابی ، میانگین و انحراف معیار 31060 پیکسل واقع در آبخوان بسیار به یکدیگر نزدیک بود، ولی با استخراج این پیکسل­ها و بررسی آن‌ها در نرم­افزار SPSS مشخص گردید، بین خروجی روش کریجینگ معمولی با روشهای IDW، LPI، RBF و TIN در سطح یک درصد آزمون­های توکی، دانکن، LSD و دانت، اختلاف معنی­دار وجود دارد. هرچند بین روش کریجینگ معمولی با کریجینگ کلی بوسیله هیچ یک از آزمون­های مذکور اختلاف معنی­دار دیده نشد. بنابراین درست است که به لحاظ همبستگی خروجی مدل‌های مختلف درون‌یابی  بسیار شبیه به یکدیگر می­باشد، ولی عملکرد این روشها به لحاظ آماری کاملاً متفاوت از یکدیگر است. برای آنکه در عمل بتوان کارکرد روشهای درون‌یابی  را مورد ارزیابی قرار داد، در گام آخر برای تمامی تصاویر سطح بدست آمده از روشهای مختلف، جهت جریان آب زیرزمینی مشخص گردید.

شکل 13- مقایسه جهت جریان تصاویر بدست آمده از روشهای مختلف در مورد داده متوسط تراز سالیانه آب سال 89

بحث و نتیجه­گیری

با مقایسه شکل­های 11 تا 13 مشاهده می­گردد، با وجود دقت عملکرد بسیار مناسب تمامی روشهای درون‌یابی  (مطابق جدول 1)، تفاوت زیادی بین خروجی تصاویر سطحی حاصل از روشهای مختلف درون‌یابی  وجود دارد. با توجه به مشاهده این روند در تمامی داده­های مربوط به مقاطع زمانی مختلف، انتخاب روش دورنیابی مناسب بر اساس استدالال­های آماری بر اساس خطاهای محاسباتی عملاً امکان­پذیر نمی­باشد و در صورت اشتباه در انتخاب روش، خروجی می­تواند بسیار با واقعیت متفاوت باشد. بنابراین برای انتخاب روش مناسب درون‌یابی  در این تحقیق از روش تحلیل­کارشناسی بر اساس واقعیت موجود، استفاده شد. بر این اساس اطلاعات مربوط به شهرک­های صنعتی، نواحی مسکونی و زمین­های کشاورزی در کل ناحیه مطالعاتی جمع­آوری شد، مشخص گردید کدام ناحیه از روش چاه جذبی برای دفع پساب استفاده می­کند. در گام دیگر بر اساس تمامی روشهای درون‌یابی  و داده­های مقاطع زمانی مختلف و با استدلال واگرایی در جهت جریان، تمامی مناطق محتمل تغذیه آبخوان مشخص شدند. سپس با مقایسه تحلیل کارشناسی اطلاعات موجود و این خروجی­ها، بهترین روش انتخاب و نقشه نواحی محتمل تغذیه تهیه گردید(شکل 14). بر اساس این نتایج، روش زمین­آمار کریجینگ کلی عملکرد بهتری از روشهای درون‌یابی دیگر نشان ­داد.

بر اساس یافته­های این پژوهش هر چند برخی از مطالعات آبهای زیرزمینی بر برتری یک روش درون‌یابی  بر روش دیگر تاکید می­کنند(1، 6، 16، 27) ولی مطابق نتایج این تحقیق، این برتری در ظاهر وجود ندارد و تنها به کمک مقایسه خروجی با واقعیت می­توان با قطعیت بالاتری در زمینه برتری یک روش درون‌یابی  به جمع­بندی نهایی رسید.

همچنین نرم­افزار مورد استفاده و توان علمی پژوهشگر در بکارگیری قابلیت استفاده از آن(9، 13، 21، 22)، تاثیر زیادی بر تصمیم­گیری نهایی خواهد داشت. تعریف هوشمندانه واریوگرام بر اساس شرایط محیط و انتخاب فضای همسایگی مطابق با آن در خروجی تصویر سطحی بسیار اثرگذار است(13)، بنابراین بسیار مهم است پژوهشگر از واقعیت­های محیطی موجود، نهایت استفاده را برده باشد.

بر اساس خروجی­ها، نرم­افزار Gstat با وجود انعطاف­پذیری بالا، در مدل کریجینگ معمولی به دلیل از بین نبردن کامل روند، توان زیادی در مدل‌سازی به شکل ناهمسان پهنه­ای ندارد. هرچند می­توان در مورد این نوع از مدل­ها  از روش کریجینگ کلی استفاده نمود. در مقایسه، نرم­افزار تحلیل­گر زمین­آمار قابلیت‌های تصمیم­گیری بیش‌تری در اختیار پژوهشگر قرار می­دهد و در مورد مدل­های ناهمسان پهنه­ای، توان بی­روندسازی کلی را در مورد روش کریجینگ معمولی دارد. همچنین امکان مشاهده اثر تعریف فضای همسایگی در نرم­افزار Gstat به طور هم‌زمان برخلاف تحلیل­گر زمین­آمار وجود ندارد، بنابراین در این قبیل مطالعات استفاده از این تحلیل­گر توصیه می­شود.

بجزء روش­های زمین­آمار، سایر روش­های درون‌یابی  نیز عملکرد مناسبی نشان دادند، بنابراین باید تاکید شود، هیچ روش درون‌یابی  بر روش دیگر بدون در نظر گرفتن شرایط محیطی برتری خاصی ندارد(6، 8)، بنابراین توصیه می­شود کارکرد تمامی روش­ها در مورد مطالعه موردی مدنظر مقایسه و بر آن اساس انتخاب صورت پذیرد.

شکل 14- نقشه مناطق محتمل تغذیه آبخوان همراه با چاه­های آب مشاهده­ای که برای برنامه پایش توصیه شده است
 

تشکرو قدردانی

این پژوهش با تصویب و حمایت مالی حوزه معاونت محیط زیست انسانی سازمان حفاظت محیط زیست و همکاری اداره کل حفاظت محیط زیست استان البرز اجرا گردیده است که بدین وسیله سپاسگزاری می­شود.

منابع

1. Nwankwoala HO, Eludoyin OS, Obafemi AA, 2012. Groundwater Quality Assessment And Monitoring Using Geographic Information Sysytems (Gis) In Port Harcourt, Nigeria. Ethiopian Journal of Environmental Studies and Management,  vol. 5,  pp. 583-96.
2. Machiwal D, Jha MK, 2015. Identifying sources of groundwater contamination in a hard-rock aquifer system using multivariate statistical analyses and GIS-based geostatistical modeling techniques. Journal of Hydrology: Regional Studies,
3. Ahmed SS, Okour YH, Said EH, 2015. A Methodology Based on Advanced Modeling Techniques for Groundwater Monitoring and Management-Part A. International Journal For Research & Development In Technology,  vol. 3,  pp. 6-12.
4. Liu L, Yi L, Cheng X, 2014. Stochastic simulation of shallow aquifer heterogeneity and it’s using in contaminant transport modeling in Tianjin plains. Journal of Water Resources and Ocean Science,  vol. 3,  pp. 80-8.
5. Huang PM, Li Y, Sumner ME. Handbook of soil sciences : properties and processes. 2nd ed. Boca Raton, FL: CRC Press; 2011.
6. Gong G, Mattevada S, O’Bryant SE, 2014. Comparison of the accuracy of kriging and IDW interpolations in estimating groundwater arsenic concentrations in Texas. Environmental Research,  vol. 130,  pp. 59-69.
7. Karimi HA. Handbook of research on geoinformatics. Hershey: Information Science Reference; 2009. xxxiv, 481 p. p.
8. Di Piazza A, Conti FL, Noto LV, Viola F, La Loggia G, 2011. Comparative analysis of different techniques for spatial interpolation of rainfall data to create a serially complete monthly time series of precipitation for Sicily, Italy. International Journal of Applied Earth Observation and Geoinformation,  vol. 13,  pp. 396-408.
9. Johnston K, Ver Hoef JM, Krivoruchko K, Lucas N. Using ArcGIS geostatistical analyst: Esri Redlands; 2001.
10. Kebaili Bargaoui Z, Chebbi A, 2009. Comparison of two kriging interpolation methods applied to spatiotemporal rainfall. Journal of Hydrology,  vol. 365,  pp. 56-73.
11. Zhang K, Li H, Achari G, 2009. Fuzzy-stochastic characterization of site uncertainty and variability in groundwater flow and contaminant transport through a heterogeneous aquifer. Journal of Contaminant Hydrology,  vol. 106,  pp. 73-82.
12. Tapoglou E, Karatzas GP, Trichakis IC, Varouchakis EA, 2014. A spatio-temporal hybrid neural network-Kriging model for groundwater level simulation. Journal of Hydrology,  vol. 519,  pp. 3193-203.
13. Oliver MA, Webster R, 2014. A tutorial guide to geostatistics: Computing and modelling variograms and kriging. Catena,  vol. 113,  pp. 56-69.
14. Özdemir Y. Determination of surface water reseources changes by Multi Temporal Satellite Imagery and Investigations of groundwater by using of geoistatistical method.  30th EARSeL Symposium: Remote Sensing for Science, Education and Culture; Paris, France2010.
15. Theodossiou N, Latinopoulos P, 2006. Evaluation and optimisation of groundwater observation networks using the Kriging methodology. Environmental Modelling & Software,  vol. 21,  pp. 991-1000.
16. Hassan J, 2014. A Geostatistical approach for mapping groundwater quality (Case Study: Tehsil Sheikhupura). International Journal,
17. Candela L, Olea RA, Custodio E, 1988. Lognormal Kriging For The Assessment Of Reliability In Groundwater Quality Control Observation Networks. Journal of Hydrology,  vol. 103,  pp. 67-84.
18. Anderson F, 2014. Multivariate Geostatistical Model for Groundwater Constituents in Texas. International Journal of Geosciences,  vol. 5,  pp. 1609-17.
19. Ghadermazi J, Sayyad G, Mohammadi J, Moezzi A, Ahmadi F, Schulin R, 2011. Spatial Prediction of Nitrate Concentration in Drinking Water Using pH as Auxiliary Co-kriging Variable. Procedia Environmental Sciences,  vol. 3,  pp. 130-5.
20. Jensen JR. Introductory digital image processing : a remote sensing perspective. 3rd ed. Upper Saddle River, N.J.: Prentice Hall; 2005. xv, 526 p. p.
21. Negreiros J, Painho M, Costa A, Santos J, Lopes I, editors. Geostatistical Analysis: Software Flashpoint. Geocomputation Conference, Dublin, Ireland; 2007.
22. Pebesma EJ, Gstat user’s manual.
23. Pebesma EJ, Wesseling CG, 1998. Gstat: a program for geostatistical modelling, prediction and simulation. Computers & Geosciences,  vol. 24,  pp. 17-31.
24. Pebesma EJ, editor Gstat: multivariable geostatistics for S. Proceedings of DSC; 2003.
25. Eastman J. IDRISI Selva manual: Clark Labs; 2012. 324 p.
26. Eastman J. IDRISI Selva Tutorial: Clark Labs; 2012.
27. Jang C-S, Chen S-K, Kuo Y-M, 2013. Applying indicator-based geostatistical approaches to determine potential zones of groundwater recharge based on borehole data. Catena,  vol. 101,  pp. 178-87.