آموزش مدل RUSLE در GIS برای برآورد فرسایش خاک (فاکتورهای R، K، LS، C، P با داده ایران)

ایران از نظر فرسایش خاک در وضعیت نگران‌کننده‌ای قرار دارد و برآوردهای مختلف، متوسط فرسایش در بسیاری از حوضه‌های زاگرس و البرز را چند برابر حد قابل تحمل (حدود ۱ تا ۵ تن در هکتار در سال) نشان می‌دهند. رسوب حاصل، مخازن سدها را پر می‌کند و خاک حاصلخیز را از بین می‌برد. مدل RUSLE به دلیل سادگی و نیاز به داده‌های در دسترس، پرکاربردترین مدل برآورد فرسایش در مطالعات آبخیزداری و پایان‌نامه‌های ایران است. با این حال، محتوای فارسی موجود یا فقط فرمول کلی را می‌دهد یا فروش کد است. در این راهنما فرمول هر فاکتور، منبع داده رایگان آن برای ایران، دستور دقیق Raster Calculator و یک مثال عددی کامل را می‌آوریم. این مطلب مکمل صفحه پژوهشی برآورد خطر فرسایش خاک در سایت است و روی اجرای عملی تمرکز دارد.

مدل RUSLE چیست و چه محدودیت‌هایی دارد؟

معادله جهانی فرسایش خاک (USLE) را Wischmeier و Smith در دهه ۱۹۶۰ و ۱۹۷۰ بر اساس هزاران کرت آزمایشی در آمریکا ارائه کردند و نسخه بازنگری‌شده آن (RUSLE) در سال ۱۹۹۷ منتشر شد. مدل فرسایش ورقه‌ای و شیاری را برآورد می‌کند، نه فرسایش آبراهه‌ای، خندقی یا توده‌ای.

آماده‌سازی: هم‌ترازی همه لایه‌ها

پیش از ساخت فاکتورها، یک «رستر مرجع» تعریف کنید: DEM حوضه در سیستم UTM (مثلاً EPSG:32639) با سلول ۳۰ متر. همه فاکتورها باید همین اندازه سلول، همین Extent و همین Snap Raster را داشته باشند؛ در غیر این صورت ضرب لایه‌ها جابه‌جایی پیکسلی ایجاد می‌کند. در ArcGIS Pro از Environments گزینه‌های Cell Size، Snap Raster و Mask را روی DEM و مرز حوضه تنظیم کنید. توضیح کامل در هم‌ترازی و بازنمونه‌گیری رستر. مرز حوضه را هم از استخراج حوضه آبریز از DEM بگیرید.

فاکتور R: فرسایندگی باران

تعریف اصلی R به شدت بارش ۳۰ دقیقه‌ای (EI30) نیاز دارد که در بیشتر ایستگاه‌های ایران موجود نیست. بنابراین از روابط تجربی بر اساس بارش سالانه یا ماهانه استفاده می‌شود. رابطه Renard و Freimund (۱۹۹۴) برای بارش سالانه P (میلی‌متر):

R = 0.0483 × P^1.61                       (P < 850 mm)
R = 587.8 − 1.219 × P + 0.004105 × P²     (P ≥ 850 mm)

برای مثال، ایستگاهی با بارش سالانه ۴۰۰ میلی‌متر: R = 0.0483 × 400^1.61 ≈ 747؛ و برای ۲۵۰ میلی‌متر حدود ۳۵۰. گام‌ها در GIS:

1. میانگین بارش سالانه دوره بلندمدت (حداقل ۲۰ سال) را برای ایستگاه‌های داخل و اطراف حوضه آماده کنید.
2. R را برای هر ایستگاه در جدول محاسبه کنید.
3. با IDW یا کریجینگ رستر R را بسازید (راهنما در مقایسه IDW و کریجینگ). در حوضه‌های کوهستانی، ابتدا بارش را با گرادیان ارتفاعی روی DEM بسازید و بعد R را از رستر بارش حساب کنید.

اگر ایستگاه کم دارید، داده ماهواره‌ای CHIRPS (حدود ۵ کیلومتر) جایگزین مناسبی است، اما در کوهستان کم‌برآورد دارد و باید با ایستگاه‌ها تصحیح شود. شاخص فورنیه اصلاح‌شده (MFI) بر اساس بارش ماهانه هم در مطالعات ایران رایج است. برای محاسبه بارش متوسط حوضه از ایستگاه‌ها، پلیگون تیسن را ببینید.

فاکتور K: فرسایش‌پذیری خاک

K به بافت، ماده آلی، ساختمان و نفوذپذیری خاک بستگی دارد. وقتی داده آزمایشگاهی کامل ندارید، رابطه Williams (مدل EPIC) فقط با درصد شن (SAN)، سیلت (SIL)، رس (CLA) و کربن آلی (C) کار می‌کند:

fcsand = 0.2 + 0.3 × exp(−0.256 × SAN × (1 − SIL/100))
fcl-si = (SIL / (CLA + SIL))^0.3
forgc  = 1 − 0.25 × C / (C + exp(3.72 − 2.95 × C))
fhisand = 1 − 0.7 × SN1 / (SN1 + exp(−5.51 + 22.9 × SN1)),   SN1 = 1 − SAN/100
K(US) = fcsand × fcl-si × forgc × fhisand
K(SI) = K(US) × 0.1317

1. لایه‌های درصد شن، سیلت، رس و کربن آلی عمق ۰ تا ۳۰ سانتی‌متر را از SoilGrids (۲۵۰ متر) دانلود کنید. دقت کنید SoilGrids مقادیر را با ضریب ذخیره کرده است (مثلاً رس بر حسب گرم در کیلوگرم، یعنی تقسیم بر ۱۰ برای درصد).
2. هر چهار لایه را به رستر مرجع ۳۰ متری بازنمونه‌گیری کنید.
3. فرمول‌ها را مرحله‌به‌مرحله در Raster Calculator اجرا کنید.

مقادیر معمول K در SI بین ۰٫۰۱ (خاک شنی یا رسی سنگین) تا ۰٫۰۶ (خاک سیلتی) است. اگر نقشه‌های خاک مؤسسه تحقیقات خاک و آب یا نمونه‌برداری میدانی دارید، آن‌ها بر SoilGrids ارجح‌اند.

فاکتور LS: محاسبه از DEM

فاکتور LS معمولاً بیشترین تأثیر را بر الگوی مکانی فرسایش دارد و حساس‌ترین بخش محاسبه است. رابطه Moore و Burch (۱۹۸۶) با جریان تجمعی کار می‌کند:

LS = (As / 22.13)^m × (sin β / 0.0896)^n      m = 0.4 تا 0.6 ، n = 1.3
As = FlowAccumulation × CellSize

1. DEM را با Fill اصلاح کنید و Flow Direction و Flow Accumulation را بسازید (مراحل در آموزش استخراج حوضه).
2. شیب را به درجه بسازید؛ راهنما در ایجاد نقشه شیب در GIS.
3. در Raster Calculator ArcGIS Pro:

Power(("FlowAcc" * 30) / 22.13, 0.4) * Power(Sin("Slope_deg" * 0.01745) / 0.0896, 1.3)

مثال: پیکسلی با جریان تجمعی ۳ سلول (As = ۹۰ متر) و شیب ۱۲ درجه: (90/22.13)^0.4 × (sin12° / 0.0896)^1.3 ≈ 1.75 × 2.99 ≈ 5.24.

در QGIS، الگوریتم SAGA «LS Factor» یا «LS-Factor, Field Based» همین محاسبه را با روش‌های مختلف (Moore، Desmet و Govers) انجام می‌دهد. شناخت دقیق DEM را در مدل رقومی ارتفاع چیست؟ و روش‌های تولید DEM بخوانید.

فاکتور C: پوشش گیاهی و مدیریت

دو روش رایج دارید:

روش ۱: از NDVI

C = exp(−α × NDVI / (β − NDVI))      α = 2 ، β = 1   (Van der Knijff و همکاران)

برای NDVI برابر ۰٫۴۵: C = exp(−2 × 0.45 / 0.55) ≈ 0.195. NDVI را از میانگین فصل رشد چند سال (نه یک تصویر) بسازید؛ روش در محاسبه شاخص NDVI. مقادیر منفی NDVI (آب) را قبل از محاسبه ماسک کنید و C را بین ۰ و ۱ محدود کنید.

روش ۲: جدول کاربری اراضی

نقشه کاربری را طبق تهیه نقشه کاربری اراضی بسازید و با Reclassify یا Join جدول، آن را به C تبدیل کنید. ترکیب دو روش (C جدول برای کشاورزی و C از NDVI برای مرتع و جنگل) معمولاً واقع‌گرایانه‌ترین نتیجه را می‌دهد.

فاکتور P: عملیات حفاظتی

در بسیاری از حوضه‌های ایران به دلیل نبود نقشه عملیات، P برابر ۱ فرض می‌شود که حالت محافظه‌کارانه است. اگر چنین فرضی می‌کنید، آن را صریحاً در گزارش بنویسید. تراس‌های باغی و دیوارچینی‌های سنتی را می‌توان از تصاویر با تفکیک بالا شناسایی کرد.

ترکیب فاکتورها و طبقه‌بندی فرسایش

"R" * "K" * "LS" * "C" * "P"

خروجی بر حسب تن در هکتار در سال است. برای پیکسل مثال ما (R = ۷۴۷، K = ۰٫۰۳۲، LS = ۵٫۲۴، C = ۰٫۱۹۵، P = ۱): A = 747 × 0.032 × 5.24 × 0.195 × 1 ≈ 24.4 t/ha/yr که در کلاس «خیلی زیاد» قرار می‌گیرد.

برای اجرای تکرارپذیر کل زنجیره، همه مراحل را در ModelBuilder بسازید و میانگین فرسایش هر زیرحوضه را با Zonal Statistics as Table استخراج کنید.

مثال حوضه: از فرسایش ناخالص تا رسوب

فرض کنید برای حوضه‌ای ۱۲۵۰ کیلومتر مربعی میانگین فرسایش ۱۴٫۶ تن در هکتار در سال به دست آمده است:

• فرسایش ناخالص کل: 125,000 ha × 14.6 = 1,825,000 t/yr.
• نسبت تحویل رسوب با رابطه Vanoni: SDR = 0.4724 × Area(km²)^−0.125 = 0.4724 × 1250^−0.125 ≈ 0.19.
• رسوب ویژه برآوردی در خروجی: 1,825,000 × 0.19 ≈ 347,000 t/yr.
• اگر ایستگاه رسوب‌سنجی خروجی حوضه بار معلق سالانه را حدود ۳۱۰ هزار تن نشان دهد، نسبت برآورد به مشاهده ۱٫۱۲ است که برای مدلی در این مقیاس قابل قبول است. اختلاف بیش از ۲ برابر نشانه مشکل در LS یا C است.

اولویت‌بندی زیرحوضه‌ها برای عملیات آبخیزداری

خروجی ارزشمند برای مدیران منابع طبیعی، نقشه پیکسلی نیست؛ فهرست اولویت‌دار زیرحوضه‌هاست. با Zonal Statistics as Table میانگین فرسایش هر زیرحوضه را بگیرید و با مساحت آن ترکیب کنید:

زیرحوضه SB1 با فقط ۱۷ درصد مساحت، بیش از یک‌چهارم فرسایش را تولید می‌کند و بازگشت سرمایه عملیات حفاظتی در آن بیشتر است. درون همین زیرحوضه، پیکسل‌های کلاس «خیلی زیاد» را با لایه کاربری روی هم بگذارید: اگر بیشترشان دیم روی شیب بالای ۱۵ درصد است، تبدیل دیم کم‌بازده به باغ یا مرتع پیشنهاد اصلی می‌شود. این تحلیل سناریو را می‌توانید با تغییر لایه C (مثلاً کاهش C دیم از ۰٫۴ به ۰٫۱) و اجرای دوباره مدل، به صورت عددی نشان دهید.

RUSLE در گوگل ارث انجین

برای حوضه‌های بزرگ یا مقایسه چندساله، کل زنجیره را می‌توان در GEE اجرا کرد: CHIRPS برای R، SoilGrids (در کاتالوگ جامعه کاربران) برای K، Copernicus DEM یا SRTM برای LS و سری زمانی سنتینل ۲ یا لندست برای C. محدودیت اصلی، محاسبه Flow Accumulation است که در GEE ابزار بومی ندارد؛ معمولاً از محصول آماده MERIT Hydro یا HydroSHEDS برای جریان تجمعی استفاده می‌شود. اگر با محیط GEE آشنا نیستید، آشنایی با گوگل ارث انجین را ببینید. برای پایان‌نامه، اجرای ترکیبی رایج است: C و R در GEE و LS و ترکیب نهایی در ArcGIS Pro.

حساسیت‌سنجی: کدام فاکتور نتیجه را تعیین می‌کند؟

چون مدل ضربی است، خطای ۲۰ درصدی در هر فاکتور مستقیماً ۲۰ درصد خطا در فرسایش می‌دهد؛ اما منبع عدم قطعیت فاکتورها یکسان نیست. یک آزمون ساده: هر فاکتور را جداگانه ۲۰ درصد افزایش دهید، یا روش محاسبه آن را عوض کنید و تغییر میانگین فرسایش حوضه را ثبت کنید.

نتیجه: LS و انتخاب‌های مربوط به DEM بیشترین اثر را دارند؛ پس آن‌ها را با دقت بیشتری مستند کنید.

چک‌لیست گزارش مدل RUSLE در پایان‌نامه

• منبع، تفکیک و دوره آماری داده هر فاکتور.
• رابطه دقیق هر فاکتور با مرجع علمی و واحد SI.
• نقشه هر پنج فاکتور با راهنمای عددی، در کنار نقشه نهایی.
• تنظیمات LS: روش، توان‌ها و حداکثر طول دامنه.
• جدول مساحت هر کلاس فرسایش و میانگین زیرحوضه‌ها.
• مقایسه با داده رسوب‌سنجی یا مطالعات پیشین و حساسیت‌سنجی.

نقشه‌ها را طبق آموزش طراحی لی‌اوت نقشه در یک صفحه شش‌قسمتی (پنج فاکتور و نتیجه) آماده کنید؛ این قالب در داوری‌ها بسیار مورد پسند است. اگر می‌خواهید فرسایش را در کنار سایر معیارها برای مکان‌یابی عملیات آبخیزداری وزن‌دهی کنید، آموزش AHP در GIS را ببینید.

خطاهای رایج در مدل RUSLE و راه حل

منبع اروپایی خوب برای روش‌ها و داده‌های فاکتورها، صفحه RUSLE2015 مرکز داده خاک اروپا (ESDAC) است.

پرسش‌های متداول درباره مدل RUSLE

مطالب مرتبط