مقدمه : 
باتوجه به نقش مواد معدنی در تامین نیازهای اولیه صنایع مختلف، جستجوی آنها با روشهای کارآمد ضروری بنظر می رسد. امروزه صرفا با بکارگیری اطلاعات زمین شناسی سطحی نمی توان منابع معدنی را جستجو نمود. روشهای ژئوفیزیکی از جمله روشهای مناسبی می باشند که درپی جوئی منابع مذکور کارایی بالایی دارند. تقریبا در تمامی مراحل عملیات اکتشاف از این روشها که ارزان، قابل اعتماد و در بسیاری موارد باعث کاهش ریسک های بزرگ سرمایه گذاری می گردند، بهره گرفته می شود.
در این علم و علوم وابسته، برداشتها به صورت سیستماتیک و تحت شرایط کنترل شده ای انجام می شوند و نتایج آنها به صورت مقادیر عددی تحت عنوان داده های مشاهده ای ارائه می شوند. برای هر نوع نتیجه گیری از این داده ها ، باید رابطه بین توزیع خواص فیزیکی سیستم مورد مطالعه (زمین) و پاسخ ژئوفیزیکی قابل مشاهده معلوم گردد. سیستم معادلاتی که این رابطه را نشان می دهد، تئوری مستقیم (Forward Theory) نامیده می شود. استنباط خواص فیزیکی سیستم مورد مطالعه از روی داده های مشاهده ای ی روند معکوس است وابزار کمی مورد نیاز برای رسیدن به آن، تئوری معکوس (Inverse Theorry) خوانده می شود. از آنجا که نتایج خام در اندازه گیری های ژئوالکتریکی نمی توانند هیچگونه اطلاعات کمی درباره عمق حقیقی ، ابعاد و یا توزیع مقاومت های ویژه حقیقی الکتریکی در زیرزمین بیان کنند ، لذا مدل سازی داده ها ، اجتناب ناپذیر است .
در این مقاله سعی برآنست تا موقعیت توده های زیرسطحی با استفاده از روش مدل سازی معکوس پارامتری تعیین گردد و نتایج به دست آمده با نتایج حاصل از مدل سازی های هموار دو بعدی که قبلاً انجام شده است (۵) مقایسه و سپس محل حفاری های پیشنهادی ارائه شود .
مدل سازی معکوس پارامتری
بدلیل سرعت پایین مدل سازی پیشرو و برای کاهش ریسک عملیات اکتشاف ، مدل سازی های معکوس ، جایگاه ویژه ای در امور اکتشافی و بخصوص مسائل ژئوفیزیکی یافته اند . در مدل سازی معکوس دادده های مشاهده ای به عنوان ورودی به سیستم داده می شوند و سپس با انجام محاسبات مربوط به تئوری معکوس سازی بر روی آنها ، پارامترهای مدل به عنوان خروجی دریافت می شوند . چون در نظر است در مطالعه حاضر از مدل سازی معکوس پارامتری (Parametric ) برای تعبیر و تفسیر داده ها استفاده گردد چگونگی آن توضیح داده می شود .
معکوس سازی داده های الکتریکی و الکترومغناطیسی یک مسأله غیر خطی است . از آنجا که در غالب موارد داده های برداشت شده کم دقت ( حاوی نویز ) غیر هماهنگ و گاهاً متناقضند ، در نتیجه تعداد بی شماری مدل مساعد برای این داده ها وجود دارد . در نتیجه پاسخ مسائل معکوس سازی یکتا نمی باشد . هدف از معکوس سازی ، مشخص کردن مدلی است که مشاهدات ما را به بهترین وجه تشریح نماید و محدودیت های وارد شده از طرف فیزیک مسأله را پاسخگو باشد . تاکنون روش های پارامتری متعددی برای حل این گونه مسائل ارائه و ابداع شده است [۷,۸,۹,۱۰,۱۳,۱۲,۱۵,۱۷,۱۸] .
روش قدرتمند کمترین مربعات برای حل مسائل غیر خطی بصورت تکراری استفاده می شود . اما در بیشترین این تکرارها ، برای همگرایی قابل قبول مدل ، نیاز به یک حدس اولیه (Initial guess ) مناسب است . البته هیچ تضمینی هم وجود ندارد که همه قسمت های طراحی شده به سمت مدل واقعی همگرا شوند . در مسائل غیر خطی داده های تجربی (d) به پارامترهای مدل (m) وابسته اند ؛ و یک تابع غیر خطی (f) (از مدل پیشرو ) موجود است ؛ لذا :
(۱) d=f(m)+e
که در آن e بردار خطا می باشد .
هدف تعیین مدلی از زمین است که پاسخی مشابه داده های اندازه گیری شده داشته باشد . این مساله را می توان با کمینه کردن یک مساله بهینه سازی ، فرمول سازی نمود ، که در آن محاسبه خطای تکراری های مختلف نیز مفید است . مثلاً اگر تعداد خطاهای ( ) داده ها ، n باشد و خطاها دارای توزیع گوسی باشند ، یک ماتریس قطری وزنی به نام W همانند رابطه زیر وجود دارد :
(۲) 
از این ماتریس برای مقایسه داده های مشاهده شده ، به منظور ممانعت از کاربرد تخمین های ضیعیف استفاده می شود ، تا اینکه داده ها یک روند استاندارد را طی کنند . حال باید مقدار خطا در رابطه (۱) کمینه شود . از آنجا که e بصورت برداری تعریف شده ، باید اندازه آن مشخص شود . پس :
(۳) SSq=||e||=eTe=(Wd-Wf(m))T(Wd-Wf)m))
برای مسائل خطی از این نوع ، روش استاندارد کمترین مربعات می تواند به صورت تکراری ، در آراستن یک مدل حدسی اولیه بکار رود . فرض این است که مدل خطی و حدس اولیه m0 است . بسط تیلور مرتبه اول برای آن نوشته می شود ؛ داریم :
(۴) 
یا
(۵) f(m)=f(m0)+Ax
 ماتریس مشتقات جزئی مربوط به پارامترهای مدل و X=(m-m0) شامل بردار پارامترهای تصحیح شده به مقادیر معلوم است .
همچنین رابطه (۳) را می توان به صورت رابطه زیر نوشت :
(۶) SSq=(Wy-WAX)T(Wy-WAX)
بردار y=(d-f(m0)) تفاوت میان جوابهای تجربی و مقادیر محاسبه شده برای m0 است .
در روش رگرسیون ریج یا کمترین مربعات تعدیل شده ، طول پرشها را با تصحیح X و بوسیله قراردادن یک شرط در رابطه فوق می توان کنترل نمود . حال باید تابع زیر کمینه شود :
(۷) =(Wy-WAX)T(Wy-WAX)+β(XTX-L2) φ
(L2) محدود کننده (Constrained ) مانع تغییر زیاد پارامترها و β فاکتور تعدیل (Damping factor ) خوانده می شود .
کمینه سازی وقتی رخ می دهد که مشتقات φ نسبت به پارامترهای قابل تصحیح برابر صفر گردد .
این عمل رابطه زیر را نتیجه می دهد .
(۸) ((WA)TWA+ βI)X=(WA)TWY
که در آن I ماتریس واحد است . و در نهایت پارامترهای تصحیح شده از رابطه زیر بدست می آید :
(۹) X=((WA)TWA+ Βi)-1(WA)TWY
از آنجان که رگرسیون ریج ، روشی برای برآورد ضرایب رگرسیون داده های غیر متعامد است که در آن ضرایب برآورد شده در مقایسه با برآوردهای کمترین مربعات با روشهای معمول ، میانگین مربعات خطای کوچکتری دارند ؛ لذا مقادیر برآورد شده با این روش ، در مقایسه با روشها معمولی ، به مقادیر حقیقی ضرایب رگرسیون نزدیکترند . از طرفی برآوردهای کمترین مربعات ریچ پایدار است و تحت تاثیر تغییرات جزئی داده ها قرار نمی گیرد (۱۰)
میزان برتری ضرایب رگرسیون برآورد شده با روش ریج بر ضرایب رگرسیون برآورد شده با روش های معمولی به مقادیر حقیقی این ضرایب در مدل بستگی دارد . چون مقادیر حقیقی مجهول می باشند ؛ استفاده از این روش در مواقعی که بین داده های هم راستایی چندگانه بالایی وجود دارد ، توصیه می شود . فرمول تکرار مورد استفاده در روش رگرسیون ریج عبارتست از :
(۱۰) mk+1=mk+((WA)TWA+ βI)-1(WA)TWY
Mk مدل بدست آمده در تکرار kام و Y,A در mk قابل ارزیابی هستند . عملیات تکرار تا زمانی ادامه می یابد که مقصود نهایی حاصل شود . در این تحقیق برای مدل سازی پارامترهای دو بعدی داده ها از نرم افزار [۱۱] RESIXIP2DI استفاده شده است .
مثال عملی
به منظر دستیابی به ذخایر مس ، با توجه به شواهد زمین شناسی و رخنمونهای دیده شده از ماده معدنی در ترانشه حفر شده در یکی از اندیس های معدنی بر دسیر کرمان ( پاینده ) ، سعی شد در مناطقی که آنومالی های اولیه گسترده تری دارند ، عملیات ژئوفیزیکی انجام شود تا نتایج معتبرتری حاصل شود .
سنگهای موجود در ناحیه را اغلب سنگهای آتشفشانی – رسوبی ائوسن تشکیل می دهند . همچنین نهشته های کربناتی و یا آواری نئوژن به صورت سنگهای آتشفشانی در منطقه مورد مطالعه رخنمون دارند . حرکات زمین ساختی در پایان میوسن میانی نهشته های موجود در منطقه را تحت تاثیر قرارداده و سبب چین خوردگی آنها با روند کلی شمال غرب – جنوب شرق گشته اند (۱) .
ظاهراً در پایان این چین خوردگی ها یا کمی پس از آن فعالیت های ماگمایی سبب شده توده های دیوریتی از راه درز و شکاف های حاصل از گسلش در سنگهای پوست نفوذ کند و متعاقب آن محلول های گرمابی نیز از طریق درزه و شکافهای حاصل از فعالیتهای تکتونیکی به درون سنگهای دربرگیرنده نفوذ کرده و از راه متاسوماتیسم با سنگهای دربرگیرنده ، کانی سازی مس را به صورت پراکنده و رگه ای سبب شوند (۴)
نتایج مطالعات بیانگر آنست که در این ناحیه کانی زایی های مس اغلب در زون آلتره و یا زون های گسله قرار دارند و در بعضی از موارد و تنها بر حسب تصادف در ولکانیک های دگرسان شده خارج این زون دیده می شوند . کانه های متداول زون اکسیده ، مالاکیت و آزوریت هستند اما کانه های کوولیت و کالکوپیریت در محل تداخل گسلها با سنگ آهک ظهور پیدا می کنند .
با توجه به شواهد فوق الذکر بنظر می رسد بهترین روش مطالعه ، روش اندازه گیری پلاریزاسیون القایی (IP ) و در درجه بعد اندازه گیری مقاومت ویژه الکتریکی (RS ) باشد طراحی برداشت و بر روی آنومالی های اولیه صورت گرفت و کار برداشت ژئوفیزیکی در ناحیه تقریباً در جهت عمود بر ساختارهای زمین شناسی عمومی منطقه و لایه های مشاهده شده در ترانشه حفر شده با سه پروفیل به طول تقریبی ۵۰۰ متر طراحی و داده های مقاومت ویژه و پلاریزاسیون القایی بطور همزمان با استفاده از دستگاه ABEM SAS1000 ادامه یافت (۴)
بحث
برای بدست آوردن تصویر درست از ساختارهای زیرسطحی در امتداد پروفیل های اندازه گیری شده لازم است کلیه داده های خام پس از اعمال تصحیحات لازم برای کاهش اثرات توپوگرافی و حذف برخی داده های ناجور مدل سازی شوند ، تا گسترش جانبی و عمقی زون های کانی سازی شده با دقت بهتری تعیین گردد . لذا ابتدا با استفاده از نرم افزار RES2DINV (14) داد ها تصحیح توپو گرافی شده سپس با ملاحظه طبیعت داده های خام و اعمال تصحیحات لازم ، مدل سازی معکوس هموار دو بعدی بصورت همزمان و ترکیبی برای داده های مقاومت ویژه و پلاریزاسیون القایی در هر یک از پروفیل ها انجام شد و نتایج مدل سازی مورد تفسیر قرار گرفت .
براساس نتایج مربوط به مدل سازی هموار دو بعدی برای یکی از پروفیل ها (PAYANDEG00 ) یک نقطه حفاری بعمق ۹۰ متر بصورت عمود بر سطح و در فاصله ۲۷۰ متری از مبدا پروفیل ( شکل ۱ ) پیشنهاد شده (۵)
اما از آنجا که در معکوس سازی هموار ، هاله ای از توده آنومالی که بتدریج به زمینه تقلیل می یابد نشان داده می شود ، تعیین محل دقیق کانی ساز و مرز افتراق توده های مختلف مشکل و گاهاً مدل سازی پارامتری در این است که در این روش با مدل سازی توام (JOINT ) داده های IP و RS مدلی ساخته می شود که پاسخ آن در تطابق خوبی با داده های اندازه گیری شده باشد .
در حالی که در مدل سازی هموار برای هر یک از سری داده های IP و RS مدلی جداگانه ارائه می گردد که ممکن است نتایج آنها با هم در تطابق خوبی نباشند (شکل ۱)
نتایج مدل سازی معکوس پارامترهای روی داده های پروفیل مذکور (payandeh00 ) در شکل ۲ و جدول ۱ نشان داده شده است .
همان طوری که از شکل و جدول فوق پیداست دو ناحیه کانی سازی شده ( توده ۱ و توده ۲ ) با مقاومت ویژه پایین در مدل وجود دارد که با توجه به شواهد زمین شناسی و اطلاعات حاصل از ترانشه های اکتشافی احتمالاً توده های حاوی کانی سازی مس می باشند .
اگر چه وزن کانی سازی شده اول دارای گسترش جانبی بیشتر و اندازه بزرگتری می باشد و در عمق کمتری نسبت به زون ۲ قرار گرفته است ولی در عوض توده دومی از قابلیت پلاریزاسیون ( بارپذیری ) بیشتری برخوردار می باشد که ویژگی مهمی برای کانی سازی مس است . بنا به دلایل فوق پیشنهاد می شود در مرحله اول یک حفاری در فاصله ۳۲۰ متری از مبدا پروفیل به عمق ۶۰ متر عمود بر سطح انجام شود تا حضور کانی سازی توده ۲ مورد بررسی قرار گیرد. علاوه برآن با توجه به عمق کم توده ۱ و گسترش جانبی زیاد آن پیشنهاد می گردد تا ابتدا ضمن حفر ی تا دو ترانشه و گرفتن نمونه های لازم آثار کانی سازی مورد مطالعه قرار گیرد و در صورت دریافت نتایج مثبت، یک گمانه در فاصله ۱۶۰ متری مبدا پروفیل حفر و تا عمق ۴۵ متری ادامه پیدا کند تا ضخامت این زون کانی سازی با دقت بیشتری تعیین گردد.
نتیجه گیری و پیشنهادات:
با توجه به تئوری ارائه شده و مثال عملی مطروحه، لازم است برای تعیین و تفسیر داده های پلاریزاسیون القایی و مقاومت ویژه ، مدل سازی معکوس روی آنها انجام شود؛ چرا که تنها دراین صورت می توان موقعیت توده مورد نظر را بخوبی تعیین و محل های حفاری را با اطمینان بیشتر مشخص نمود. بررسی ها نشان می دهد که روش ها ی معکوس سازی هموار نسبت به روش های پارامتری به نویز حساسیت کمتری دارد. در صورتی که اگر از رگرسیون ریج در مدل سازی پارامتری استفاده شود؛ پایداری این روش در حضور نویز افزون می گردد و با یک حدس اولیه ضعیف هم به سوی مدل واقعی همگرا می شود.
مقایسه مدل های بدست آمده با مدل سازی هموار و مدل بدست آمده با مدل سازی پارامتری نشان می دهد که اولا توده اصلی کانی سازی (با باپذیری بالا) بطور مایل (در شکل ۱) در حد فاصل ۲۰۰تا۳۲۰ متری از مبدا پروفیل و در بازه عمقی ۴۰تا۷۰ متری از سطح قرار دارد در حالی که آثار همین توده کانی سازی در شکل ۲ (توده ۲) در یک گستره ۱۰۰ متری (از ۲۷۰ تا ۳۷۰ متری پروفیل و با عمق متوسط ۴۵ متری وجود دارد در حالی که آثار همین توده در شکل ۱ در یک گستره بیشتر (از ۸۰ تا ۳۳۰ متری ) و با حداکثر عمق ۵۵ متر مشاهده می شود.
نتایج بررسی های صورت گرفته موید این است که مدل سازی های معکوس هموار اگرچه در تطابق خوبی با نتایج مدل سازی پارامتری می باشند ولی حد وحدود مناطق کانی سازی شده به دست آمده با روش مدل سازی پارامتری دقیق تر بوده و از این نظر برای پیشنهادهای مناسب تر حفاری می باشند. در نهایت آنجا که با استفاده از روش پلاریزاسیون القایی نمی توان نوع کانی سازی فلزی یا غیرفلزی را تشخیص داد، پیشنهاد می شود در اکتشاف چنین کانی سازی هایی در صورت امکان از روش پلاریزاسیون القایی طیفی (SIP) بمنظور تفکیک کانی سازی فلزی از غیرفلزی استفاده شود[۲و۳و۵].
منابع :
مرادزاده،علی- عرب امیری، علیرضا
[۱] سازمان زمین شناسی و اکتشافات معدنی ایران،۱۹۹۵ . نقشه زمین شناسی ۱۰۰۰۰۰: ۱ چهارگوشه چهار گنبد
[۲] عرب امیری، ع، ۱۳۷۹. ارزیابی کانیسازی روی و سرب در غرب معدن تپه سرخ از مجتمع معادن ایرانکوه به کمک روشهای زمین شناسی، ژئوشیمیایی و ژئوفیزیکی،پایان نامه کارشناسی ارشد، دانشگاه صنعتی اصفهان
[۳] عرب امیری،ع، ۱۳۷۹. کاربرد آنالیز طیفی داده های پلاریزاسیون القایی در اکتشاف کانسازی فلزی، سمینار کارشناسی ارشد، دانشگاه صنعتی اصفهان
[۴] مرادزاده ،ع،و عرب امیری،ع.۱۳۸۲. طرح پژوهشی مطالعات ژئوفیزیکی مقاومت ویژه و پلاریزاسیون القایی منطقه معدنی بردسیر، دانشگاه صنعتی شاهرود
[۵] مرادزاده ،ع. و عرب امیری،ع.۱۳۸۲. تعبیر وتفسیر و مدل سازی دو بعدی داده های مقاومت ویژه و پلاریزاسیون القایی اندیس معدنی پاینده، بیست و دومین گردهمایی علوم زمین، سازمان زمین شناسی و اکتشافات معدنی کشور.
[۶] Beogard-Yugoslavia, 1973. Exploration for ore deposits in Kerman region. Geological survey of Iran. [7] Broyden, C.G., 1965. A class of methods for solving nonlinear simultaneous equations. Mathematics of computation, 577-593.
[8] Constable, S.C., Parker, R.L. and Constable, C.G., 1987. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52, 289-300.
[9] Inman, J.R., Ryu, J. and Ward, S.H., 1973. Resistivity inversion. Geophysics, 38, 1088-1108. [ 10] Inman, J.R., 1975. Resistivity inversion with ridge regression. Geophyslcs, 40, 789-797.
[ 1 1] Interpex, Ltd., 1999. Users manual of RESlX IP2DI, 2-D resistivity and induced polarization data interpretation software. Interpex Ltd., USA.
[12] Loke, M.H. and Barker, R.D., 1996. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method. Geophysical prospecting, 131-152.
[ 13] Loke, M.H., 1999. Electrical imaging surveys for environmental: A practical guide to 2-D and 3-D surveys and engineering studies.
[ 14J Loke, M.H., 2001. Manual of RES2DINV, rapid 2-D resistivity and IP inversion using the least square method, GEOTOMO software.
[ l5] Marquardt, D.W., 1963. An algorithm for least -squares estimation of nonlinear parameters. J. Soc. Indust. Appl. Math., 11.
[16] Meju, M.A., 1992. An effective ridge regression procedure for resistivity data inversion. Computer and geosciences, 18, 99-118.
[17] Smith, N.C. and Vozoff; K., 1984. Two-dimensional DC resistivity inversion for dipole-dipoles data. lEEE transactions on geosciences remote sensing, GE-22.
[18] Vozoff, K. and Jupp, D.L.B., 1975. Joint inversion of geophysical data. Geophys.J.Roy.Astro. Soc., 977991