مقایسه روش‌های درون‌یابی مکانی جهت پهنه‌بندی غلظت فلزات سنگین در خاک سطحی شهرستان آران و بیدگل

علوم و تکنولوژی محیط زیست، دورههجدهم، شماره یک، بهار 96

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

یونس سیفی^[1]

روح اله میرزایی^[2]*

 rmirzaei@kashanu.ac.ir

چکیده

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

روش بررسی: در مجموع 135 نمونه خاک سطحی از شهرستان آران و بیدگل جمع‌آوری شد و غلظت کادمیوم، روی و مس در نمونه‌ها تعیین گردید؛ سپس با استفاده از روش­های مختلف درون­یابی شامل کریجینگ معمولی، کوکریجینگ، وزن­دهی معکوس فاصله، چند جمله­ای موضعی و توابع پایه شعاعی، تغییرات مکانی غلظت عناصر در خاک سطحی شهرستان آران و بیدگل ارزیابی شد. ارزیابی متقابل و آماره­های MSE، MBE،RMSE ، NSE و PBIAS به منظور صحت­سنجی روش­ها استفاده شد.

یافته­­ها: طبق نتایج ارزیابی متقابل روش کریجینگ معمولی (مدل کروی) دارای بهترین کارایی برای برآورد غلظت کادمیوم و مس و روش کریجینگ معمولی (مدل نمایی) دارای بهترین کارایی برای برآورد غلظت روی در خاک این منطقه بود. همچنین با افزایش توان وزن­دهی در روش چند جمله­ای موضعی (LPI)، برای عناصر کادمیوم و مس، صحت درون­یابی افزایش و برای عنصر روی، صحت درون­یابی کاهش یافت؛ در حالی­که با افزایش توان وزن­دهی در روش وزن­دهی معکوس فاصله [3](IDW)، صحت درون­یابی برای عناصر کادمیوم و روی کاهش و برای عنصر مس افزایش یافت.

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

 واژه­های کلیدی: درون­یابی مکانی، زمین­آمار، آلودگی خاک، فلزات سنگین، آران و بیدگل.

Comparison of Spatial Interpolation Methods to Mapp Heavy Metals Concentrations in Surface Soil of Aran-O-Bidgol City

Younes seifi^[4]

 Rouhollah Mirzaei^[5]

rmirzaei@kashanu.ac.ir

Abstract

Background and Objective: Selecting a suitable interpolation method is a main factor of surface analysis and it is a difficult task in geostatistical analysis since different methods of interpolation can result in different surfaces and ultimately different results. Hence, the aim of this study is to evaluate the performance of different interpolation methods in determining the spatial concentration of Cd, Zn and Cu in the Surface soil of Aran-O-Bidgol City.

Method: A total of 135 surface soil (0-20 cm) samples were collected in the study area and the soil heavy metals concentrations were determined. The spatial distribution of heavy metals concentration in Surface Soil of Aran-O-Bidgol was evaluated using different interpolation methods including ordinary Kriging, CoKriging, Inverse Distance Weighting, Local Polynomials and Radial Basis Function. Cross validation and MSE, MBE, RMSE, NSE, PBIAS were applied to estimate their accuracy.

Findings: According to the results Ordinary Kriging (spherical model) had the best efficiency for estimating Cd and Cu concentrations and Ordinary Kriging (exponential model) had the best efficiency for estimating distribution pattern of Zn concentration in the soil of this region. The weight parameter has a significant influence on the accuracy of interpolation. The higher the order of local polynomial, the larger the accuracy of cross validation, whereas, the greater the weighting power of IDW for Cd and Zn, the greater error and larger accuracy of the interpolation will be for Cu.

Discussion and Conclusion: Although all of the considered interpolation methods had a high prediction accuracy of the mean content for soil heavy metals, but taking all cross validation indicators into consideration, Ordinary Kriging (OK) method shows a generally better performance than other methods.

Keywords: Spatial Interpolation, Geostatistics, Soil Pollution, Heavy Metals, Aran-O-Bidgos

مقدمه

الگوی مکانی متغیرهای محیطی به طور فزاینده­ای در مدیریت و علوم محیط­زیست مورد نیاز است. این اطلاعات معمولاً به آسانی در دسترس نیستند و دستیابی به آنها مخصوصأ در مناطق کوهستانی یا مناطق دریایی عمیق بسیار سخت و هزینه­بر است. توزیع مکانی متغیرهای محیطی اغلب از طریق منابع نقطه­ای اندازه­گیری می­شود. با این حال، مدیران محیط­زیستی برای تصمیم­گیری موثر و مطمئن، به داده­های مکانی پیوسته از سراسر ناحیه مورد نظر نیاز دارند. بنابراین روش­های درون­یابی مکانی، برای تخمین متغیرهای محیطی در مکان­های نمونه­برداری نشده به ابزاری ضروری تبدیل شده­اند. خصوصیات خاک از نقطه­ای به نقطه دیگر دارای تغییرات مکانی هستند و این تغییرات ممکن است به شکل ذاتی (تفاوت در مواد مادری) و یا غیر ذاتی (تفاوت در مدیریت منابع خاکی) روی دهد (1). دانستن خصوصیات خاک، به ویژه خصوصیاتی که از نظر تولید کشاورزی و مسائل اقتصادی مهم هستند، دارای اهمیت بسیار است (2). با توجه به انباشت مستمر فلزات­ سنگین در خاک، بررسی توزیع غلظت عناصر سنگین جهت پایش آلودگی خاک و حفظ کیفیت محیط‌زیست ضروری است. تغییر­پذیری خصوصیات خاک با این فرض که توزیع خصوصیات خاک در عرصه به صورت تصادفی است، اغلب توسط روش­های آماری کلاسیک بیان می­شود. در این روش­ها نتایج به دست آمده از اندازه­گیری نمونه­ها، مستقل از موقعیت فضایی آنها مورد بررسی قرار می­گیرد. بنابراین مقدار یک کمیت در یک نمونه هیچگونه اطلاعاتی درباره مقدار همان کمیت در نمونه­های دیگر به فواصل مختلف ارائه نمی­کند. در روش­های زمین­آمار می­توان بین مقادیر یک کمیت در جامعه نمونه­ها و فاصله و جهت قرار گرفتن نمونه­ها نسبت به یکدیگر ارتباط برقرار کرد. بنابراین در این روش ابتدا به بررسی وجود یا عدم وجود ساختار مکانی بین داده­ها پرداخته می­شود و در صورت وجود ساختار مکانی، تحلیل داده­ها انجام می­شود (3). انتخاب یک روش درون­یابی مناسب برای داده­ها آسان نیست و برای رسیدن به سطح پیوسته، هنوز روش جامعی وجود ندارد و روش­های درون­یابی مکانی در حال تکامل هستند و نمی­توان یک روش را برای هر منظور و هر منطقه تجویز کرد، چون عوامل مختلفی مانند: طراحی نمونه­برداری و توزیع فضایی نمونه­ها، ماهیت و کیفیت داده­ها، ارتباط بین متغیرهای اولیه و ثانویه، ارتباطات بین فاکتورهای مختلف و موارد دیگری از جمله انتخاب مدل مناسب تغییرنما و اندازه شبکه­ها روی صحت روش­ها تأثیرگذار می­باشند (4). به عبارت دیگر، هیچ پاسخ ساده­ای در رابطه با انتخاب یک روش درون­یابی مناسب وجود ندارد، چون " بهترین روش" فقط برای موقعیت­های خاص است (5). اکثر روش‌های درون‌یابی در مطالعات مربوط به پلات‌های آزمایشگاهی استفاده ‌شده و به ندرت مطالعات وابسته به مقیاس‌های بزرگ گزارش شده است. به هر حال برای پیش‌بینی روند مکانی آلاینده‌ها، رسیدن به نتایج توزیع مکانی مناسب به وسیله فقط یک روش سنتی از درون‌یابی آن هم به وسیله داده‌های جدا از هم مشکل است که در بررسی‌های اخیر روش‌های گوناگون درون‌یابی برای سنجش توزیع مکانی فلزات سنگین خاک استفاده می‌شود (6). در زمینه بررسی تغییرات مکانی خصوصیات خاک با استفاده از روش‌های درون‌یابی مکانی تحقیقات زیادی انجام شده است، اما نتایج به خوبی واضح و مشخص نیست. برخی از آن‌ها بیان داشتند که روش کریجینگ معمولی بهتر از سایر روش­ها عمل می‌کند (7، 8، 9، 10)، درحالی‌که برخی دیگر نشان دادند که کریجینگ معمولی بهتر از سایر روش‌ها نیست (11،12،13،14).

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

مواد و روش­ها

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

منطقه مورد بررسی، شهرستان آران و بیدگل واقع در شمالی­ترین ناحیه استان اصفهان و در همسایگی شهرستان­های کاشان، قم، نطنز و اردستان است (شکل1). این شهرستان با مساحت 6051 کیلومتر مربع در طول­های جغرافیایی 50 درجه و 15 دقیقه تا 52 درجه و 29 دقیقه طول شرقی و عرض­های جغرافیایی 33 درجه و 30 دقیقه تا 34 درجه و 27 دقیقه عرض شمالی و در حاشیه کویر مرکزی ایران قرار دارد. متوسط بارش سالیانه 100 میلی­متر و اقلیمی با تابستان­های گرم و خشک و زمستان­های سرد و خشک را دارا می­باشد. ارتفاع متوسط شهر از سطح دریا 912 متر می­باشد. این شهرستان با ضریب خشکی 4/5 جزء منطقه فراخشک محسوب می­گردد. عمده کاربری­های این منطقه دربرگیرنده­­­ کاربری صنعتی، مسکونی، اراضی بدون پوشش، زراعت آبی و باغات، جاده­ها و جنگل دست کاشت است. وجود صنایع گوناگون و کوره­های آجرپزی متعدد از یک سو و زمین­های کشاورزی از سوی دیگر، چهره­ای صنعتی- کشاورزی به شهر داده است. اراضی کشاورزی شهرستان آران و بیدگل بالغ بر 13500هکتار است که از این میزان حدود 2000 هکتار اختصاص به کشت محصولات دائمی و باغی داشته و حدود 11500هکتار زیر کشت انواع محصولات زراعی قرار دارد. از جمله مهم­ترین محصولات کشاورزی این شهرستان می­توان به گندم، جو، پنبه، پسته، صیفی، علوفه و ذرت اشاره نمود.

شکل 1-  موقعیت شهرستان آران و بیدگل در استان اصفهان و نقاط نمونه­برداری خاک در سطح شهرستان

Figure 1- The location of study area in Isfahan Province and the distribution of the soil samples

2. نمونه­برداری و آنالیز آزمایشگاهی

طی مطالعه میدانی 135 نمونه خاک سطحی از تمام سطح شهرستان آران و بیدگل در سلول­های هم­اندازه و از محل تقاطع شبکه­ها جمع­آوری شد. نمونه­برداری به روش سیستماتیک شبکه­ای صورت گرفت. پس از نقطه­یابی توسطGPS ، 5 زیرنمونه از هر گروه شبکه از عمق 20-0 سانتی­متری جمع­آوری شد، به این صورت که یک نمونه از نقطه تقاطع جمع­آوری شد و مابقی 4 نمونه از 4 قسمت به شعاع 50 متری از نقطه مرکزی نمونه­برداری برداشت شد. 5 زیرنمونه با هم مخلوط شدند تا نمونه ترکیبی خاک به دست آید. در شکل(1) موقعیت و تعداد نقاط نمونه‌گیری در منطقه مورد مطالعه آورده شده است. نمونه‌ها پس از خشک شدن در هوای آزاد (هوا خشک)، کوبیده شده و از الک 2 میلی‌متری عبور داده شدند. سپس قسمتی از هر نمونه خاک، با مش 100 الک شد. یک گرم از خاک الک شده با استفاده از ترکیب سه اسید HNO₃، HF و HClO₄ با نسبت 3:2:5 هضم شد و سپس محلول به دست آمده پس از عبور از کاغذ صافی واتمن شماره 42 با استفاده از آب دیونیزه به حجم 25 میلی­لیتر رسانده شد (15)؛ سپس غلظت عناصر کادمیم، روی و مس در نمونه‌های خاک با دستگاه جذب اتمی شعله (FAAS)  مدل Shimadzu, AA-670 مورد اندازه‌گیری قرار گرفت. خاک استاندارد SRM 2711 Montana II برای کنترل کیفیت و بررسی صحت اندازه­گیری مقدار فلزات سنگین کل در نمونه­ها استفاده شد. لازم به ذکر است در هنگام هضم به همراه هر گروه از نمونه­ها یک نمونه شاهد تهیه و همراه با سایر نمونه­ها مورد آنالیز قرار گرفت.

3. روش­های درون­یابی

3-1.کریجینگ

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

در این معادله z*(x) تخمین مقدار متغیر z در نقطه x و  وزن آماری اختصاص یافته به مقادیر z در نقطه x_iاست. شرط استفاده از این تخمین­گر آن است که متغیر z توزیع نرمال داشته باشد. در صورتی­که متغیر مورد نظر توزیع نرمال نداشته باشد، باید از کریجینگ غیرخطی استفاده نمود، یا با استفاده از روش­های تبدیل داده­ها، توزیع متغیر مورد نظر را تبدیل به نرمال نمود (16، 17). در این پژوهش از تابع نیم تغییرنما جهت نشان دادن تغییرات یک متغیر با در نظر گرفتن فاصله استفاده شد که معادله آن به صورت زیر است:

y(h) ²(2)                                                                            

که در آن y(h) مقدار نیم تغییرنما در فاصله h، Z(x_i+h) و Z(x_i) مقدار متغیر در نقاط (x_i+h) و x_i، N(h) تعداد جفت نمونه­های به­کار رفته به ازای هر فاصله h می­باشد. غالباً به 2y(h) واریوگرام و به y(h) سمی­واریوگرام اطلاق می­شود (18). در این معادله مقدار تغییرنما وابسته به فاصله بین مقادیر یک متغیر ناحیه­ای در دو نقطه است. اگر این مقدار وابسته به جهت نیز باشد، نیم تغییرنمای ناهمسانگرد و در غیر این صورت همسانگرد نامیده می­شود. در این پژوهش از مدل­های کروی[6]، نمایی[7]، گوسی[8]، خطی[9] و خطی دارای سقف[10] برای برازش نیم تغییرنما استفاده می­شود.

3-2. کوکریجینگ

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

(3)                                                                z*(x)   

 که در آن z*(xi) مقدار تخمین­ زده شده برای نقطه xi،  وزن مربوط به متغیر z،  وزن مربوط به متغیر کمکی y، z(xi) مقدار مشاهده شده متغیر اصلی و y(xk) مقدار مشاهده شده برای متغیر کمکی می­باشد.     

3-3. چند جمله­ای موضعی 

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

3-4. توابع پایه شعاعی    

از جمله روش­های درون­یابی است که در آن سطح تخمین از مقادیر مشاهده­ای عبور می­کند. از خصوصیات این روش که حالتی از شبکه عصبی مصنوعی می­باشد، این است که مقادیر بیشتر از حداکثر مشاهده­ای و یا کمتر از حداقل مشاهده­ای در سطح تخمین وجود دارد. در روش وزن­دهی معکوس فاصله نیز سطح تخمین از مقادیر مشاهده­ای عبور می­نماید، اما این روش هیچ­گاه مقادیر بیشتر از حداکثر و کمتر از حداقل داده­های مشاهده­ای را تخمین نمی­زند. روش تابع پایه شعاعی بر اساس پنج تابع مختلف درون­یابی داده­ها را انجام می­دهد که هسته اصلی توابع، مجموع مقادیر مربعات (h²+R²) است که به ترتیب، R عامل هموار کننده[11] و h ترسیم ناهمسانگرد[12] مرتبط با فاصله نقاط هستند (20). توابعی که در این تحقیق مورد استفاده قرار می­گیرد عبارت از: چند ربعی[13]، چند ربعی معکوس[14] و نواری کم ضخامت[15] هستند.

3-5. وزن­دهی معکوس فاصله

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

که در آن:  وزن نقطه i ام، Di فاصله نقطه i ام تا نقطه مجهول و α معادل توان وزن­دهی می­باشد.

4.2. اعتبارسنجی روش­های درون­یابی

به منظور مقایسه روش­های استفاده شده در این پژوهش و انتخاب مناسب­ترین روش درون­یابی، از روش اعتبارسنجی متقابل[16] استفاده شد. در این روش، در هر مرحله یک نقطه مشاهده­ای حذف شده و با استفاده از بقیه نقاط، آن نقطه برآورد می­شود. این کار برای همه نقاط مشاهده­ای تکرار می­شود، به طوری­که در آخر به تعداد نقاط مشاهده­ای، برآورد وجود خواهد داشت. همچنین با استفاده از مدل­های به­دست آمده، در تعدادی از نقاطی که اندازه­گیری وجود خواهد داشت، مقادیر تخمینی محاسبه می­شود. در نهایت با توجه به مقادیر مشاهده شده و برآورد شده، دقت هر روش با توجه به معیارهای آماری میانگین مجذور خطا (MSE)، میانگین اریب خطا (MBE)، خطای برآورد ریشه دوم میانگین مربعات (RMSE)، شاخص کارایی (NSE) Nash–Sutcliffe و درصد اریب خطا (PBIAS) محاسبه ­گردید. مقدار این معیارها با استفاده از روابط 5 تا 9 به دست آمد:  

در این معادلات z^*(x_i) مقدار برآورد شده در نقطه x_i، z (x_i) مقدار اندازه­گیری شده در نقطه x_i، o میانگین مقادیر اندازه­گیری شده و n تعداد نقاط می­باشد. در این تحقیق برای انجام تحلیل­های آماری از نرم­افزار SPSS نسخه 22 و برای انجام تحلیل­های زمین­آماری از نرم­­افزار­هایGS+  نسخه 1/5 و ArcGIS نسخه 1/10 استفاده شد.

نتایج و بحث

1.3. توصیف آماری داده­ها

قبل از انجام محاسبات زمین­آماری، نرمال بودن توزیع داده­ها توسط آزمون کولموگروف- اسمیرنوف بررسی و مشخص شد که غلظت عناصر کادمیم و مس در خاک از توزیع نرمال تبعیت می­کنند، ولی غلظت عنصر روی در خاک توزیع نرمال ندارد. داده­های این عنصر با روش لگاریتمی نرمال شدند. آمار توصیفی متغیرهای مورد مطالعه در جدول (1) آمده است. همان طور که نتایج نشان می­دهد مقدار کادمیم، روی و مس در منطقه مورد مطالعه به ترتیب از 30/0 تا 15/1 میلی­گرم در کیلو­گرم، 50/21 تا 75/165 میلی­گرم در کیلو­گرم و 45/1 تا 22/35 میلی­گرم در کیلوگرم در نوسان می­باشد. ضریب تغییرات غلظت کادمیم، روی و مس در منطقه به ترتیب: 25،1/45 و 9/43 درصد می­باشد، که نشان از تغییرپذیری متوسط غلظت کل فلزات مذکور در منطقه دارد. ضریب تغییرات کمتر از 10 درصد نشان از تغییرپذیری کم و ضریب تغیرات بیشتر از 90 درصد نشان از تغییرپذیری زیاد دارد (22).

جدول 1- خلاصه آماری غلظت کادمیم، روی و مس در خاک سطحی شهرستان آران و بیدگل (135n=)

Table 1- Descriptive statistical summary of the total concentrations of Cd, Cu and Zn in the topsoils of Aran-O-Bidgol (n=135)

2.3 برآورد غلظت عناصر با استفاده از روش­های درون­یابی

برای بررسی وجود یا عدم وجود همبستگی مکانی غلظت فلزات مذکور در منطقه مورد مطالعه طبق روش کریجینگ معمولی، پس از برازش چندین مدل نیم تغییرنما، بر پایه آماره مجذور مربعات خطا (RSS) و حداکثر ضریب همبستگی ((R²، برای فلزات کادمیم و مس، مدل کروی و برای فلز روی مدل نمایی به عنوان مناسب­ترین مدل برگزیده شد (شکل 2). پس از تهیه تغییرنمای سطحی متغیرها به منظور بررسی همسانگردی داده­ها، مشخص شد که متغیرها همسانگرد بوده و می­توان در محاسبات بعدی از تغییرنمای همه­جهته استفاده نمود. شکل­ (2) تغییرنمای همه­جهته متغیرهای کادمیم، روی و مس را نشان می­دهد. علاوه بر این، مشخصات مربوط به این تغییرنماها (مدل برازش شده، شعاع تأثیر، آستانه، اثر قطعه­ای و ...) در جدول (2) نشان داده شده است. از بررسی نسبت اثر قطعه­ای به آستانه می­توان نتیجه گرفت که هر سه متغیر مورد مطالعه دارای ساختار مکانی (وابستگی مکانی) قوی می­باشند، زیرا هرچه این نسبت کمتر باشد، نشان از آن دارد که ساختار مکانی بهتری برای متغیر وجود دارد. با توجه به نسبت اثر قطعه­ای به آستانه می­توان گفت: اگر این نسبت کمتر از 25/0 باشد، متغیر مورد نظر دارای ساختار مکانی قوی، بین 25/0 تا 75/0 ساختار مکانی متوسط و هنگامی که بزرگتر از 75/0 باشد، ساختار مکانی ضعیف است (23). این امر می­تواند نشان­دهنده غلبه بیشتر واریانس بخش ساختاری تغییرنما بر واریانس بخش تصادفی آن در مورد همه فلزات در تحقیق حاضر باشد که خود می­تواند دلیلی بر مناسب بودن الگوی نمونه­برداری، فواصل کوچک­تر نمونه­برداری و دقت نسبتاً مناسب تجزیه­های آزمایشگاهی باشد. در این برازش مقدار اثر قطعه­ای برای کادمیم 0001/0، روی 0001/0 و برای مس 001/0 به­دست آمد. در حالت ایده­ال باید مقدار اثر قطعه­ای صفر باشد، اما در واقعیت نیم­تغییرنماهای تجربی مقادیر بالاتر از صفر را نشان می­دهند، که این به دلیل ظهور جزء تصادفی متغیر می­باشد. اثر قطعه­ای به علت وجود تغییرات در فواصل نمونه­برداری و یا به دلیل خطای نمونه­برداری و اندازه­گیری بروز می­کند (3).

جدول 2- پارامترهاو نتایج ارزیابی متقابل بهترین مدل­های تغییرنمای برازش داده­شده به غلظت عناصر

Table 2- Cross-validation performance and the best-fitted semivariogram models and their parameters for soil heavy metals

شکل 2- تغییرنماهای تجربی غلظت عناصر کادمیوم، روی و مس در خاک سطحی

Figure 2- Experimental and fitted variogram models of the concentration of Cd, Cu and Zn in the surface soil

وجود همبستگی مناسب بین متغیرها لازمه روش کوکریجینگ است، به عبارت دیگر این روش در شرایطی می­تواند کارایی داشته باشد که بین متغیر مورد نظر و متغیر وابسته همبستگی معنی­داری وجود داشته باشد. در این تحقیق سه متغیر کادمیم، روی و مس به عنوان متغیر مستقل و سه متغیر pH، هدایت الکتریکی و درصد کربن آلی به عنوان متغیر وابسته در نظر گرفته شدند. نتایج محاسبه همبستگی خطی بین تمامی متغیرها به منظور تعیین بهترین متغیر کمکی در جدول (3) ارائه شده است. همان طور که نتایج نشان می­دهد، کربن آلی با کادمیم (در سطح 05/0)، روی (در سطح 01/0) و مس (در سطح 01/0) رابطه معنی­داری داشته است. بنابراین درصد کربن آلی به عنوان متغیر کمکی مناسب برای کادمیم، روی و مس انتخاب شد. پس از انجام درون­یابی­ها با استفاده از روش­های مورد­ نظر(کریجینگ معمولی، کوکریجینگ، چند جمله­ای موضعی، توابع پایه شعاعی و وزن­دهی معکوس فاصله)، آماره­های ارزیابی صحت روش­های درون­یابی­، برای متغیرها تعیین گردیدند. هرچه مقادیر میانگین مجذور خطا (MSE)، میانگین اریب خطا (MBE)، خطای برآورد ریشه دوم میانگین مربعات (RMSE) و درصد اریب خطا (PBIAS) به صفر نزدیک­تر باشند، موید آن است که صحت درون­یابی­ها بیشتر بوده است (13، 24، 25). همچنین هر چه مقدار شاخص کارایی (NSE) Nash–Sutcliffe به یک نزدیک­تر باشد، نشان­دهنده صحت بالای درون­یابی­ها می­باشد (24). مقایسه روش­های درون­یابی بر اساس این شاخص­ها در جدول (4) نشان داده شده­اند. در بین روش­های مختلف درون­یابی، روش کریجینگ معمولی (مدل کروی) برای عناصر کادمیم (0016/0RMSE=) و مس (039/0RMSE=) و روش کریجینگ معمولی (مدل نمایی) برای عنصر روی (978/0RMSE=) دارای کمترین میزان خطا و انحراف بود. از بین 14 روش مورد استفاده در این مطالعه، برای عناصر کادمیم، روی و مس، به ترتیب روش­های کوکریجینگ (464/2RMSE=)، چندجمله­ای موضعی (درجه4) (81/13RMSE=) و وزن­دهی معکوس فاصله (توان1) (422/9RMSE=) بالاترین میزان خطا را نشان دادند و نامناسب­ترین روش شناخته شدند (جدول4). نتایج به­دست آمده از جدول (4) نشان می­د­هد که با افزایش توان وزن­دهی در روش چند جمله­ای موضعی، برای عناصر کادمیوم و مس صحت درون­یابی افزایش و برای عنصر روی صحت درون­یابی کاهش می­یابد. در حالی­که با افزایش توان وزن­دهی در روش وزن­دهی معکوس فاصله، صحت درون­یابی برای عناصر کادمیوم و روی کاهش و برای عنصر مس افزایش می­یابد. همچنین طبق نتایج به­دست آمده از ارزیابی متقابل (جدول4) رتبه­بندی روش­ها (کریجینگ معمولی، کوکریجینگ، وزن­دهی معکوس فاصله، توابع پایه شعاعی و چند جمله­ای موضعی) بر اساس RMSE به قرار زیر است: برای عنصر کادمیوم پس از روش کریجینگ معمولی (مدل کروی)، به ترتیب روش­های توابع پایه شعاعی (مدل چندربعی)، وزن­دهی معکوس فاصله (توان1)، چند جمله­ای موضعی (توان4) و کوکریجینگ، برای عنصر روی پس از روش کریجینگ معمولی (مدل نمایی)، روش­های توابع پایه شعاعی (مدل چندربعی)، کوکریجینگ، وزن­دهی معکوس فاصله (توان1) و چند جمله­ای موضعی (توان4)، برای عنصر مس پس از روش کریجینگ معمولی (مدل کروی)، روش­های چند جمله­ای موضعی (توان4)، کوکریجینگ، توابع پایه شعاعی (مدل نواری کم ضخامت) و وزن­دهی معکوس فاصله (توان1) در رتبه­های بعدی قرار گرفتند.

روش­های کریجینگ در تخمین توزیع مکانی فلزات سنگین خاک، پتانسیل بالایی دارند و به عنوان روش مناسب برای درون­یابی و تهیه نقشه­های آلاینده­ها پیشنهاد شده­اند (26). این موضوع توسط مطالعات متعددی تاکنون تایید شده است از این
جمله در آزمایشی، موسوی و همکاران (27) برای پهنه­بندی غلظت فلزات آنتیموان، مس و کروم از روش کریجینگ معمولی (مدل نمایی) و برای فلز آرسنیک از روش کریجینگ معمولی (مدل کروی) استفاده کردند. برای ارزیابی تغییرات مکانی برخی از فلزات سنگین در خاک سطحی استان گلستان، میرزایی و همکاران (28) نیز از روش کریجینگ معمولی استفاده کردند، نتایج تحلیل زمین­آماری نشان داد که مدل نمایی، دایره­ای، نمایی و گوسی به ترتیب بهترین مدل برای تعیین تغییرپذیری مکانی کادمیم، آرسنیک، کروم و مس می­باشد. همچنین در مطالعات لادو[17] و همکاران (29)، دیانی و همکاران (30)، خداکرمی و همکاران (31) و محمودی و همکاران (32) نیز از روش کریجینگ معمولی استفاده شد و نتایج قابل قبولی به­دست آمد. به نظر می­رسد که روش مناسب زمین­آماری در برآورد یک متغیر، به نوع متغیر و عوامل منطقه­ای تأثیر­گذار بر آن بستگی دارد و نمی­توان روش منتخب در یک منطقه را به سایر مناطق تعمیم داد (21، 33).

جدول3- ضریب همبستگی خطی بین متغیرهای مطالعه شده

Table 3. Correlation coefficient matrix of variables in the studied soils

         ** همبستگی معنی­دار در سطح 01/0 و * همبستگی معنی­دار در سطح 05/0

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

شکل3- نقشه توزیع مکانی غلظت کادمیوم، روی و مس در خاک سطحی شهرستان آران و بیدگل

Figure 3- Spatial distributions of Cd, Cu and Zn concentrations (mg kg^-1) in surface soils of Aran-O-Bidgol

طبق نتایج به­دست آمده از این پژوهش، روش کوکریجینگ برای عناصر مورد مطالعه کارایی ضعیفی نشان داده است، به طوری­که از بین روش­های استفاده شده در این پژوهش، روش کوکریجینگ برای عناصر روی و مس در رتبه سوم و برای عنصر کادمیوم در رتبه آخر قرار گرفته است. از آنجاکه وجود همبستگی بالا بین متغیرها می­تواند باعث افزایش دقت روش کوکریجینگ شود؛ دلیل کارایی ضعیف روش کوکریجینگ در این مطالعه، می­تواند همبستگی ضعیف بین متغیرهای مورد مطالعه باشد (13، 34). به طور کلی، روش‌های درون‌یابی تعیین می‌کنند که چگونه اطلاعات نمونه‌های مجزا به نقشه‌های پیوسته تبدیل می‌شوند. دقت درون‌یابی بستگی به این دارد که روش درون‌یابی چگونه تغییرپذیری مکانی و همبستگی خصوصیات خاک را نشان می‌دهد (35). کریجینگ معمولی یک فیلتری است که جزئیات اطلاعات محلی را تعدیل می‌کند. در ناحیه‌ای با تغییرپذیری بالا که ارتباط مکانی ضعیف دارد، اثرات نرم‌سازی کریجینگ عادی خیلی نیرومند­تر است (36). صحت درون‌یابی مفهومی نسبی است، معیارها با اهداف درون‌یابی تغییر می‌کنند. دو هدف اصلی برای نقشه­سازی آلودگی خاک، تحلیل الگوی مکانی آلودگی خاک و شناسایی مناطق آلوده است. در هدف تحلیل الگوی مکانی آلودگی، نتایج پیش‌بینی روند مکانی کلی آلودگی خاک تا حدامکان باید دقیق باشد. شکی نیست که کریجینگ عادی دارای بیش‌ترین توانایی به منظور پیش‌بینی روند کلی آلودگی خاک است. به هر حال در هدف شناسایی مناطق آلوده، نیاز است تا تکنیک­های درون‌یابی خصوصیات محلی آلودگی خاک (به ویژه لکه‌های داغ و سرد محلی) را با دقت بیشتری پیش‌بینی کنند. مشخص است که تمام نتایج درون‌یابی دارای خطا هستند. شناسایی ناحیه‌ای به عنوان منطقه آلوده‌شده نباید انحصاراً متکی به نتایج درون‌یابی باشد. پیشنهاد می‌شود که زمینه طبیعی و فعالیت‌های انسانی قبل از تصمیم‌گیری در نظر گرفته شود. عدم قطعیت در ارزیابی آلودگی عمدتاً در ناحیه‌ای با تغییرپذیری محلی زیاد رخ می‌دهد (35). بنابراین نمونه‌برداری بیشتر در منطقه‌ای که دارای عدم قطعیت است توصیه می‌شود.

نتیجه­گیری

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

سپاسگزاری

این مطالعه بخشی از طرح پژوهشی به شماره 92011972 صندوق حمایت از پژوهشگران کشور (INFS) است. بدین وسیله از صندوق حمایت از پژوهشگران کشور سپاسگزاری می‌شود.

منابع

1. محمد­زمانی، س.، ایوبی، ش.ا و خرمالی، ف. بررسی تغییرات مکانی خصوصیات خاک و عملکرد گندم در بخشی از اراضی زراعی سرخنکلاته استان گلستان. مجله علوم و فنون کشاورزی و منابع طبیعی، 1386، شماره 40، صفحات 79 تا 92.
2. ایوبی، ش. ا.، محمد زمانی، س و خرمالی، ف. برآورد مقدار ازت کل خاک به کمک ماده آلی و با استفاده از روش­های کریجینگ، کوکریجینگ و کریجینگ- رگرسیون در بخشی از اراضی زراعی سرخنکلاته استان گلستان. مجله علوم کشاورزی و منابع طبیعی، 1386، جلد 14، صفحات 78 تا 87.
3. حسنی پاک، علی اصغر. زمین­آمار (ژئواستاتیستیک)، چاپ اول، انتشارات دانشگاه تهران، 1377، صفحه 314.
4. Mueller T.G., Pierce F.J., Schabenberger O., and Warncke D.D. 2001. Map quality for site-specific fertility management. Soil Sci. Soc. Am. J, vol. 65, pp. 1547–1558.
5. Isaaks E.H., and Srivastava R.M. 1989. Applied Geostatistics. Oxford University Press, New York.
6. Hengl T. 2007. A Practical Guide to Geostatistical Mapping of Environmental Variables. JRC Scientific and Technichal Reports. Office for Official Publication of the European Communities, Luxembourg.
7. کاظمی پشت مساری، ح.، طهماسبی سروستانی، ز.، کامکار، ب.، شتایی، ش و صادقی، س. ارزیابی روش­های زمین­آمار جهت تخمین و پهنه­بندی عناصر غذایی پرمصرف اولیه در برخی اراضی کشاورزی استان گلستان. نشریه دانش آب و خاک، 1391، جلد 22، شماره 1، صفحات 201- 218.
8. محمودی، ج.، زارعیان، ف.، جوادی، م. ر و خرسندی، ن. مقایسه کارایی چند روش زمین­آماری برای تخمین برخی ویژگی­های فیزیکی خاک. نشریه حفاظت منابع آب و خاک، 1391، سال اول، شماره 4، صفحات 67- 77.
9. موسوی فرد، س. م.، ممتاز، ح. ر و خداوردیلو، ح. کارایی برخی تخمین­گرهای زمین­آماری در میان­یابی و پهنه­بندی شماری از ویژگی­های کیفی خاک. نشریه حفاظت منابع آب و خاک، 1392، سال دوم، شماره 3، صفحات 57- 71.
10. Panagopoulos T., Jesus J., Antunes M.D.C., and Beltr J. 2006. Analysis of spatial interpolation for optimising management of a salinized field cultivated with lettuce. Europ J Agron, vol. 24, pp.1-10.
11. Chen T.B., Xie Y., Lei M., Yang J., Guo Q., Song B., and Zhou X. 2011. Spatial distribution of soil heavy metal pollution estimated by different interpolation method: Accuracy and uncertainty analysis. Chemosphere, vol. 82, pp. 468-476.
12. Jiachun S., Haizhen W., Jianming X., Jianjun W., Xingmei L., and Haiping Z. 2007. Spatial distribution of heavy metal in soil: A case study of changing, China. Environ Geol, vol. 52, pp.1-10.
13. نورزاده حداد، م.، مهدیان، م. ح و ملکوتی، م. ج. مقایسه کارایی برخی روش­های زمین­آماری به منظور بررسی پراکنش مکانی عناصر ریزمغذی در اراضی کشاورزی (مطالعه موردی: استان همدان). نشریه دانش آب و خاک، 1392، جلد 23، شماره 1، صفحات 71- 81. 
14. Robinson T.P., and metternicht G. 2006. Testing the performance of spatial interpolation techniques for mapping soil properties. computer and electronic in agriculter, vol. 50, pp. 97-108.
15. Li X.Y., Liu L.J., Wang Y.G., Luo G.P., Chen X., Yang X.L., GAO B., and He X.Y. 2012. Integrated assessment of heavy metal contamination in sediments from a coastal industrial basin. NE China, PLoS One.
16. محمدی، ج. پدومتری، آمار مکانی. نشر پلک، 1385، جلد دوم، صفحه 453.
17. حسنی پاک، علی اصغر. زمین­آمار (ژئواستاتیستیک)، چاپ دوم، انتشارات دانشگاه تهران، 1386، صفحه 314.
18. سمیعی، م.، کرمی، م و زارع، ل. ارزیابی روش­های زمین­آماری در پراکنش مکانی هدایت الکتریکی آب زیرزمینی در دشت فیروزآباد، اولین همایش ملی حفاظت و برنامه­ریزی محیط­زیست، 3 اسفندماه 1391.
19. Ruppert D.1996. Local polynomial regression and its applications in environmental statistics.
20. School of Operations Research and Industrial Engineering, Cornell University, New York.
21. Carlson R.E., and Foley T.A. 1991. The parameter R² in multiquadric interpolation. Computers and Mathematic Apply, vol. 21, pp. 29-42.
22. تقی­زاده مهرجردی، ر. ا.، زارعیان جهرمی، م.، محمودی، ش.، حیدری، ا و سرمدیان، ف. بررسی روش­های درون­یابی مکانی جهت تعیین تغییرات مکانی ویژگی­های کیفی آب­های زیرزمینی دشت رفسنجان. مجله علوم و مهندسی آبخیزداری ایران، 1387، سال 2، شماره 5، صفحات 63 تا 70.
23. Zhang X.Y., Sui Y.Y., Zhang X.D., Meng., and Herbert, K. S.J. 2007. Spatial variability of nutrient properties in black soil of northeast China. Pedosphere, vol. 17(1), pp.19–29.
24. Cambardella C.A., Moorman T.B., Novak J.M., Parkin T.B., Karlen D.L., Turco R.F., and Konopka  A.E. 1994. Field-scale variability of soil properties in Central Iowa soils. Soil Sci Soc Am J, vol. 58, pp.1501-1511.
25. Paul D.W., Peter F., Florian W., Shamita K., and Karl S. 2012. Comparison and evaluation of spatial interpolation schemes for daily rainfall in data scarce regions. Journal of Hydrology,vol. 464–465, pp.388–400.
26. Silvio J.G., Jonathan AL., Dennis W.H., Yann P., Jean C., and Jacques G. 2014. Mapping soil hydraulic conductivity and matric potential for water management of cranberry: Characterisation and spatial interpolation methods. biosystems engineering, vol. 128, pp.29-40.
27. Juang K.W., Lee D.Y., and Ellsworth T. R. 2001. Using rank-order geostatistics for special
28. interpolation of highly skewed data in heavy metal contaminated site. J. Environ. Qual, vol. 30, pp.894-903.
29. موسوی، ا.، سفیانیان، ع.، میرغفاری، ن و خداکرمی، ل. بررسی توزیع مکانی برخی فلزات سنگین در خاکهای استان همدان. مجله پژوهش­های خاک (علوم خاک و آب)، 1390، الف، جلد 25، شماره 4، صفحات 323 تا 336.
30. میرزایی، ر.، اسماعیلی ساری، ع.، قربانی، ه.، حافظی مقدس، ن.، همامی، م. ر و رضایی، ح. ر. پیش­بینی توزیع مکانی کادمیوم، آرسنیک، کروم و مس در خاک سطحی استان گلستان. مجله پژوهش­های محیط­زیست، 1392، سال 4، شماره 7، صفحات 35 تا 44.
31. Lado L.R., Hengl T., and Reuter H.I. 2008. Heavy metal in European soils: A geostatistical analysis of the foregs Geochemical database. Geoderma, vol. 148, pp.189-199.
32. دیانی، م.، محمدی، ج و نادری، م. تجزیه و تحلیل زمین­آماری غلظت سرب، روی و کادمیوم در خاک­های حومه سپاهان­شهر واقع در جنوب اصفهان. نشریه آب و خاک (علوم و صنایع کشاورزی)، 1388، جلد 23، شماره 4، صفحات 67 تا 76.
33. خداکرمی، ل.، سفیانیان، ع.، میرغفاری، ن.، افیونی، م و گلشاهی، ا. پهنه­بندی غلظت فلزات سنگین کروم، کبالت و نیکل در خاک­های سه زیر حوزه آبخیز استان همدان با استفاده از فناوری­های GIS و زمین­آمار. مجله علوم و فنون کشاورزی و منابع طبیعی، علوم آب و خاک، 1390، سال 15، شماره 58، صفحات 243 تا 254.
34. محمودی، ش.، محمدی، ج و نادری، م. توزیع آماری و مکانی برخی عناصر سنگین در کلاس­های اندازه­ای ذرات خاک سطحی جنوب شهر اصفهان. مجله پژوهش­های حفاظت آب و خاک، 1392، جلد 20، شماره 2، صفحات 1 تا 22.
35. دلبری، م.، خیاط خلقی، م و مهدیان، م. ح. ارزیابی روش­های زمین­آمار در برآورد هدایت هیدرولیکی خاک در مناطق شیب آب و پشت آب پایین دشت سیستان. مجله علوم کشاورزی ایران، 1383، جلد 35، شماره 1، صفحات 1 تا 13.
36. جعفری، م.، عسگری، ح. م.، معظمی، م.، بی نیاز، م و طهمورث، م. بررسی توزیع مکانی برخی از خصوصیات خاک با کاربرد روش­های زمین آماری. مجله پژوهش و سازندگی در زراعت و باغبانی، 1387، شماره 80، صفحات 177 تا 185.
37. Xie Y., Chen T., Lei M., Yang J., Guo Q., Song B., and Zhou X. 2011. Spatial distribution of soil heavy metal pollution estimated by different interpolation methods: Accuracy and uncertainty analysis. Chemosphere, vol. 82, pp.468–476.

38. Goovaerts P. 2000. Estimation or simulation of soil properties? An optimization problem with conflicting criteria. Geoderma, vol. 97, pp.165–186.