نشريه زمين شناسي مهندسي,جلد دوم، شماره ١بهار و تابستان ١٣٨٦

رهيافت مقايسهاي بر دو روش مدلسازي در پوسته البرز مرکزي

اميرهوشنگ رجايي، محمد مختاري:
پژوهشگاه بين المللي زلزله شناسي و مهندسي زلزله کيت پريستلي : دانشگاه کمبريج – انگلستان) تاريخ: دريافت ٢٥/٥/٨٥ پذيرش ٢٢/١٠/٨٦

چکيده
در ايـن مقاله نتايج مدل سازي ساختار سرعت امواج لرز هاي در پوسته با استفاده از روش
1268732149348

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

تحليل توابع گيرنده و روش تحليل برگردان ادغام ده توابع گيرنده با دادههاي پاشندگي سرعت
گروه امواج ريلي با نگرش مقايسه اي ارائه شده است. دادههاي مورد استفاده نگاشتهاي وقايع
دورلـرزه اي ثبـت شـده در يـک ايستگاه باند پهن نصب شده در منطقة البرز مرکزي را شامل
مـي شـود. مـ دل سرعت ب ه دست آمده از هر دو روش مدل سازي بيان گر تطابق و همخواني در
موقعيـت عمقـي ناپيوسـتگيهاسـت در حالـي که تفاوت هاي چشم گيري را در مقدار سرعت
لـرز هاي نشان مي دهد. همچنين تفاوت در نتايج مدل محاسبه اي حاصل از کاربرد ساختار هاي
مختلف سرعتي براي مدل اوليه در مدل سازي به روش تحليل برگردان توابع گيرنده يکتا نبودن
ايـن روش را مورد تأييد قرار مي دهد، در حالي که نتايج حاصل براي مقادير سرعت حاصل ازروش دوم داراي روند بسيار منطقي و با ثبات است.
مقدمه
روش تحليل برگردان توابع گيرنده١ به عنوان روش دقيقي كه قادر است ساختار سرعتي
پوسته را در يك بعد و زير محل هر ايستگاه لرزه نگاري مدل سازي كند شناخته شده است.

1. Receiver Function Inversion
1268732149348

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

توابع گيرنده نسبت به ناپيوستگيهاي سرعت امواج لرزهاي در پوسته حساسيت زياد ولينسبت به مقادير مطلق ساختار سرعت امواج لرز هاي حساسيت كم تري دارد. اين نكته بدينمعني است كه ساختار سرعت به دست آمده براي پوسته با روش تحليل توابع گيرنده به تنهايي مي تواند تا حد زيادي به مدل سرعت اوليه در آغاز مدل سازي وابسته باشد. در واقع نتايجحاصل از کاربرد اين تکنيک مي تواند تنوعي از ترکيب مدلهاي مختلف با مقادير سرعتي وضخامتي مختلف باشد. از طرف ديگر، اطلاعات سرعت فاز و سرعت گروه به دست آمده ازدادههاي امواج سطحي به مقادير متوسط سرعت مطلق امواج برشي حساسيت بيش تري داردتا ناپيوستگيهاي سرعت امواج لرز هاي در پوسته، لذا مدل سرعتي به دست آمده از اين روشبه صورت قابل توجهي مستقل از ساختار سرعتي مدل اوليه است. انتظار مي رود که يكمدل سازي سرعتي كه به طور رضايت بخشي از هر دو روش مدل سازي پردازش برگردان توابعگيرنده و پردازش برگردان سرعت فاز و يا گروه امواج سطحي را در برگيرد، هم مستقل ازمدل اوليه سرعتي و هم داراي قدرت تفكيك خوب در ناپيوستگيهاي لرز هاي واقع در اعماقمختلف پوسته باشد. از آن جا که در مناطقي با ساختار هاي پيچيده نظير رشته کوههاي البرزانتخاب روش مدل سازي مناسب باعث ارتقا قابل توجهي در نتايج بررسي هاي پوسته اي و لذاافزايش قطعيت در مقادير به دست آمده مي گردد، مقايسه بين دو روش مدل سازي بر اساسدادههاي منطقه مورد توجه قرار گرفت.
تحليل توابع گيرنده
يك تابع گيرنده از بخشي از يك لرزه نگاشت به دست مي آيد كه نشان دهندة پاسخ پوستةواقع در زير يك ايستگاه لرزه نگاري در مقابل عبور امواج P و يا S رويداد هاي دورلرزه اي است. تحليل توابع گيرنده در حوزة زمان براي مدل کردن فاز هاي تبديلي موج برشي وچندگانههاي٢ ثبت شده روي مؤلفه افقي شکل موجهايP است (Ammon, 1991). وقتي موجP از ميان پوسته و گوشتة بالايي عبور مي كند، با ناپيوستگيهاي سرعتي مواجه مي شودكه در اين ناپيوستگيها موجP مي تواند بازتابيده، انتقال يافته و يا به موجSV تبديل گردد.

2- Multiples
با فرض اينكه مسافت بين رخداد زمين لرزه(چشمه) و ايستگاه لرزه نگاري(گيرنده) بزرگ تر از ٣٠ درجه باشد(˚٣٠ ≥ ∆)، زاويه ورودي امواجP به سنگ كره در زير مكان ايستگاه بهوضعيت عمودي بسيار نزديك شده و لذا اين موج به صورت بسيار آشكار و قوي روي مؤلفةقائم لرزه نگار ثبت مي شود. از سوي ديگر امواج برشي تبديل شده به صورت بسيار واضح وقوي روي مؤلفههاي افقي لرزه نگار ثبت مي شود. در تحليل توابع گيرنده هر دو مؤلفه افقيلرزه نگاشت به گون هاي تبديل دوراني٣ مي يابد تا در امتداد موازي و امتداد عمود بر مسير موجاز گيرنده تا ايستگاه قرار گيرد. سيگنال ثبت شده توسط يك لرزه نگار مي تواند يك هم آميخت٤ در حوزة زمان و متشکل از سه عبارت اصلي به صورت زير ديده شود:
X i (t) = S(t)* I(t)* Ei (t) i=Z,R,T (۱)
1268732149348

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

كه در آنT , R , Z به ترتيب مؤلفههاي قائم، شعاعي و مماسي حركت،S(t) اثراتانتشار جبه و مجاور چشمه٥ ،I(t) پاسخ دستگاهي وE(t) پاسخ زمين درست در زير ايستگاه(ساختار سنگ كره در زير گيرنده) است. لانگستون (١٩٧٩) نشان داد كه به دليل زاويه تقريبﹰاقائم امواج ورودي، فاز هايSH رخداد هاي دورلرزه اي اثر بسيار ضعيفي روي مؤلفه قائملرزه نگاشت دارند. بنا بر اين:
X Z (t) ≈ S(t)*I(t) (۲)
يك تابع گيرنده كه در واقع پاسخ زمين در زيرگيرنده است مي تواند از طريق واهم آميخت٦ پاسخهاي دستگاهي و تابع چشمه كه به مؤلفة قائم مربوط مي شود با لرزه نگاشتهاي شعاعيو مماسي جداسازي شود. شكل١ فاز هايي كه به عنوان يك تعامل مابين پوسته تك لايه اي واقعدر بالاي يك نيم فضاي گوشته بالايي، به يك ايستگاه لرزه نگاري مي رسد را نشان مي دهد.
و Langston, 1977, ، Burdick and Langston, 1977:همچنين مراجعه کنيد به) ( Ammon, 1991

Rotation
Convolution
Near Source and Mantle Propagation Effects
Deconvolution
تكنيكهاي مختلفي تاکنون براي محاسبة توابع گيرنده معرفي شده است كه از جمله آنمي توان به فرايند هم ترازسازي در حوزة فركانس٧ (Clayton and wiggings, 1976)، رهيافت حوزه زمان٨(Sheehan etal,1995) و همچنين روش تكرار واهم اميخت٩ (Pablo Ligorria and Ammon, 1999) اشاره نمود. هر تكنيكي خود منفرداﹰ داراي
نكات قوت و ضعف است، ولي تكنيكهاي مختلف مزيتهاي خود را در وضعيتهايمتفاوت دارا هستند.
در اين تحقيق از روش تكرار واهم آميخت معرفي شده توسط پابلو ليگوريا و آمون(١٩٩٩) در تهيه توابع گيرنده استفاده شده است. براي هر يك از رخداد هاي دورلرزه اي انتخاب شدهابتدا يك پنجره زماني٩٠ ثانيه كه مركز آن درست روي رسيد موج مستقيمP (بازه زماني ٤٥- تا ٤٥+ ثانيه) منطبق مي شود انتخاب، و روي هرسه مؤلفة نگاشتهاي ثبت شده اعمال گرديد.
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

در اين تكنيك پردازش، يك شاخص درصدي از عدم انطباق بين مؤلفه افقي لرزه نگاشت وهم آميخت تابع گيرنده مربوط به آن با مؤلفه قائم لرزه نگاشت محاسبه مي شود. اگر اين شاخص عدم ان طباق متجاوز از١٠% باشد تابع گيرنده براي مرحله بعدي پردازش مورد استفاده قرارنخواهد گرفت.
-275081-267742شكل١. مثال مصنوعي از فاز هاي لرزهاي كه به عنوان يك تعامل مابين پوسته تك لايه اي واقع در بالاي يك نيم فضاي جبه بالايي، به ايستگاه لرزه نگاري مي رسد.
(aمجموعه اي از فاز هاي لرز ه اي که در مدل تک لايه اي در زير گيرنده توليد شده اند، b) توابع گيرنده شعاعي متناظر با مدل ساده شده تک لايه براي زمين،
c) مسير پرتو بين چشمه(کانون رخداد
دورلرز ه اي) و گيرنده(لرزه نگار)(Ammon, 1991).

Water level Stabilised frequancy Domain
Time Domain Approach
Iterative Deconvolution Method
اندازه گيريهاي پاشندگي امواج سطحي
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

براي درك ساختار متوسط پوسته و گوشته بالايي مي توان از اطلاعات به دست آمده ازاندازه گيريهاي پاشندگي امواج سطحي١٠ استفاده كرد. اندازه گيري سرعت گروه و سرعت فازامواج سطحي اطلاعاتي را در باره ساختار متوسط سرعتي از محيطي که موج از آن عبورمي کند شامل شده و به طور خاصي نسبت به ساختار سرعت موج برشي نيز حساس است. اندازه گيريهاي پاشندگي ابتدا با استفاده از روش اندازه گيري فاصله زماني مابين دو قله موج١١ توسط برون و همكاران(١٩٦٠) براي لرزه نگاشتهايي که به صورت آنالوگ ثبت شده و بهخوبي نيز پاشنده شده بودند مورد استفاده قرارگرفت. براي دادههاي ديجيتال، تحليل فركانس- زمان را لوشين و همكاران(١٩٩٩) معرفي كردند که همچنين توسط ساييل و عثمان شاهين(٢٠٠٠) نيز مورد استفاده قرار گرفت. همچنين روش دو ايستگاهي را كه گومبرگ ارائه كرد نيزمي تواند در استخراج اطلاعات پاشندگي سرعت فاز براي ساختار سرعت انتشار موج درمحدوده سنگ كره مابين دو ايستگاه زلزله نگاري مورد استفاده قرار گيرد(Priestley and Tilman,1999) .
وقتي مؤلفههاي فركانسي مختلف يك موج با سرعتهاي مختلف حركت مي كنند، قطارموج به يك شكل موج پالس گونه پهن تبديل مي شود(پراشيده مي شو د) به گونه اي كه مؤلفههايفركانسي مختلف به صورت مكاني از هم جدا مي شوند. به هر مؤلفه فركانسي يك سرعت فازوابسته است كه بيان گر سرعت موج در آن فركانس خاص است. وقتي شكل موجي كه ابتدابه صورت يك پالس انرژي است شروع به پاشندگي مي كند، مؤلفههاي فركانسي مختلف سيگنال شروع به تداخل با يكديگر كرده و به غير از زمانهاي مشخص و مخصوصي كهبه عنوان سرعت گروه موج تعريف مي شود، در بقية زمانها به طور كلي باعث از بينرفتن(ميرايي) انرژي مي شود؛ براي مثال، اگر فرض كنيم كه يك سيگنال لرز ه اي فقط انرژي دوفركانس ω1 و ω2 را شامل شود خواهيم داشت:
U(x,t) = Cos(ω1t − k1x) +Cos(ω2t − k2x) (۳)

Surface Wave Dispersion Measurements
Peak- trough
كه در آنU(x,t) جابه جايي مشاهده شده مربوط به سيگنال لرز ه اي،ω1= ω−δ ω ، k1= k− δ k ، ω2= ω+δ ω و k1= k+ δ k وω وk به ترتيب فركانس و عدد موج
متوسط هستند. سپس با جاي گذاري در رابطه خواهيم داشت:
U(x,t) = 2Cos(ω −tkx)Cos(δkx −δωt) (۴)
شكل موج منتج شده شامل يك سيگنال با يك فركانس متوسطω كه دامنه آن با فركانس
δω مدوله شده است. سرعت موج براي سيگنال با پريود كوتاه تر برابر است باc=ω / k يا سرعت فاز و براي پريود بزرگ تر برابر است با u=δω / δk يا سرعت گروه است
.(Lay and Wallace,1995)
تئوري پردازش برگردان
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

پردازش برگردان ١٢ توابع گيرنده و همچنين برگردان ادغام شده١٣ توابع گيرنده با دادههايپاشندگي سرعت گروه امواج ريلي براي پيدا كردن يك مدل سرعتي با استفاده از برنامةكامپيوتري زلزله شناسي تهيه شده هرمان و آمون(٢٠٠٣) به انجام رسيده است. پردازش برگردان تا يافتن يك مقدار حداقل براي تابع هدف زير ادامه مي يابد:
16459150
S =[(1− p)N r + pN s ]


که در آنS مقدار خطاي استاندارد مابين دادههاي پيش بيني شده ١٤ و مقادير مشاهد ه اي١٥، p يك مقدار وزن دهي، i Ori امين تابع گيرنده مشاهده شده،j Osj امين اندازه گيري پاشندگيامواج سطحي مشاهده شده، Pri وPrj مقادير پيش بيني شده متناظر براساس مدل سرعتيفعلي است.NT و NS نيز به ترتيب مقدار توابع گيرنده ومقادير پاشنده امواج سطحي و σri و σsj محدودههاي خطاي استاندارد ١٦به ترتيب برايiامين و jامين بخش از تابع گيرنده و
اطلاعات امواج سطحي است. براي يك مجموعه بزرگ داده با خطاي استانداردي كه به

Inversion
Inversion
Predicted Data
Obsereved Data
Standard Error Limits
درستي به آن اختصاص داده شده، حداقل مقدار مورد انتظارS=1.0 است . نسخه فعلي ازنرم افزار ب هكار رفته قادر به درنظر گرفتن اندازه گيريهاي خط اي استانداردσr براي توابعگيرنده ن مي باشد، لذا اين خطا معادل با محدودههاي±1SD (انحراف معيار ) مربوط به يكسري از توابع گيرنده برانبارش شده١٧ است.
براي اين كه دستورالعمل روش پردازش برگردان ادغام شده موفق باشد، ضروري است كهدادههايي كه براي دو روش به كار برده مي شود و مربوط به يك منطقه يك سان از سنگ كرهباشد. براي توابع گيرنده، دامنه تغييرات جانبي براي نمونهها به اين بستگي دارد كه چه مقدارفاز هاي تبديل شده فاز موج P مستقيم را دنبال مي كنند.
پردازش دادهها
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

دادههاي صحرايي مشتمل بر نگاشت وقايع دورلرزه اي مربوط به يک ايستگاه باند پهن سهمؤلفه اي نصب شده در منطقه البرز مرکزي به مختصات’١٣ ˚٥١ شرقي و ‘٠٤ ˚٣٦ شمالي است که طي سه مرحله مورد پردازش قرار گرفته است. در مرحله اول فرمت ركورد هاي ثبت شده با لرزه نگار از حالت فشرده خارج شده و رخدادهاي دورلرزه اي با بزرگاي بيش از ٥ ريشتر ازرکورد هاي ثبت شده استخراج شد. جدول١ مشخصات رخداد هاي انتخاب شده را نشانمي دهد. در مرحله دوم مؤلفههاي قائم و افقي توابع گيرنده توليد گرديد. در تهيه توابع گيرندهاز روش تكرار واهم آميخت استفاده شده است. سپس توابع گيرنده توليد شده که داراي كيفيتمطلوب بوده انتخاب و برانبارش شد. در مرحله سوم ابتدا مؤلفههاي شعاعي توابع گيرندهبرانبارش شده، با استفاده از روش تحليل برگردان١٨ پردازش شد و ساختار سرعتي پوسته دريك بعد و زير محل هر ايستگاه لرزه نگاري مدل سازي گرديد. در اين روش ابتدا يك مدلاوليه سرعت بر اساس اطلاعات مختلف براي ساختار پوسته زير ايستگاه مورد نظر فرضمي شود، سپس نگاشت تابع گيرندة مصنوعي متناظر با اين مدل اوليه سرعتي به دست آمده و بانگاشت تابع گيرنده واقعي مقايسه مي شود. اختلاف بين مقادير محاسبه شده و مقادير واقعيبراي اعمال تصحيح روي مدل اوليه سرعت به کار رفته و نگاشت تابع گيرندة مصنوعي مجدداﹰ

Stacked Receiver Functions
Inversion
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

براي مدل جديد سرعت به دست مي آيد، سپس اين نگاشت با نگاشت تابع گيرنده واقعيمقايسه و اختلاف آن مجدداﹰ براي تصحيح مدل سرعت مورد استفاده قرار مي گيرد. اين حلقهمحاسبات تا حصول بهترين برازش بين نگاشت مصنوعي و نگاشت واقعي تابع گيرنده تکرارمي شود. مدل سرعتي وابسته به نگاشت مصنوعي تابع گيرنده در آخرين مرحله از تکرارمحاسبات و يا به عبارت ديگر آخرين مدل سرعتي محاسبه شده به عنوان خروجي محاسباتدر نظر گرفته مي شود. در مدل سازي به روش ديگر دادهها به روش تحليل برگردان ادغام شدهتوابع گيرنده با دادههاي پاشندگي سرعت گروه امواج ريلي كه يکي از كامل ترين روشهاي دربررسي هاي پوسته اي با توجه به دادههاي موجود است مورد پردازش قرار گرفت.
جدول١ . مشخصات رخداد هاي دورلرزه اي انتخاب شده.
YEAR MO DA ORIG TIME Lot. Lon. DEP MAGNITUDE Scale
2003 5 21 184420.10 36.96 3.63 12 6.90 Ms
2003 5 21 185110.30 36.97 3.81 10 5.7 Mb
2003 5 26 092433.40 38.85 141.57 68 7.00 Mw
2003 5 26 192327.94 2.35 128.85 31 7.10 Ms
2003 5 26 231329.72 6.76 123.71 565 6.80 Mw
2003 5 26 232631.89 6.81 123.88 583 5.90 Mb
2003 6 10 084030.83 23.52 121.63 44 6.00 Mw
2003 6 15 192433.15 51.55 176.92 20 6.50 Mw
2003 6 16 220802.14 55.49 160.00 174 6.90 Mw
2003 6 23 121234.47 51.44 176.78 20 7.00 Ms
2003 6 28 152942.26 -3.33 146.15 10 6.3 Mw
2003 7 1 055225.92 4.53 122.51 653 6.00 Mw
2003 7 2 235226.28 42.32 144.84 23 5.80 Mb
2003 7 4 071644.72 76.37 23.28 10 5.70 Mb
2003 7 15 202750.53 -2.60 68.38 10 7.6 Ms
2003 7 19 212037.01 -8.68 111.23 56 5.90 Mb
2003 7 21 135358.49 -5.48 148.85 189 6.3 Mw
2003 7 23 163837.21 -15.59 -13.35 10 5.80 Mb
2003 7 24 102321.67 0.14 124.58 73 5.80 Mw
2003 7 25 093745.84 -1.53 149.69 24 6.4 Mw
2003 7 25 221329.97 38.42 141.00 6 6.10 Mw
YEAR MO DA ORIG TIME Lot. Lon. DEP MAGNITUDE Scale
2003 7 27 062531.95 47.15 139.25 470 6.80 Mw
2003 8 11 001909.28 1.14 128.15 10 6.00 Mw
2003 8 31 090132.12 10.54 146.35 57 5.80 Mb
2003 9 5 012301.96 5.32 95.90 124 5.90 Mw
2003 9 21 181613.41 19.92 95.67 10 6.90 Ms
2003 9 25 205713.39 41.73 143.64 33 5.7 Mb
2003 9 25 210800.03 41.77 143.59 33 7.40 Ms
2003 9 26 062657.15 42.16 144.67 33 5.90 Mb
2003 9 26 203822.10 41.99 144.58 33 6.00 Mw
2003 9 27 113325.08 50.04 87.81 16 7.50 Ms
2003 9 27 185246.98 50.09 87.76 10 6.60 Ms
2003 9 27 233606.90 44.67 150.35 33 5.90 Mb
2003 9 29 023653.14 42.45 144.38 25 6.50 Mw
2003 10 1 010325.24 50.21 87.72 10 7.10 Ms
2003 10 8 090655.34 42.65 144.57 32 6.70 Mw
2003 10 8 231517.49 42.21 144.69 33 5.90 Mb
2003 10 9 221913.85 13.76 119.94 33 6.2 Mw
2003 10 11 000849.14 41.92 144.36 33 5.9 Mb
2003 10 11 011131.18 43.97 148.21 51 6.2 Mb
2003 10 18 222713.25 0.44 126.10 33 6.4 Mw
2003 10 22 114530.84 -6.06 147.73 53 6.3 Mw
2003 10 25 124135.25 38.40 100.95 10 5.8 Mb
2003 10 25 124758.83 38.38 100.97 10 5.8 Mw
2003 10 28 214821.01 43.84 147.75 65 6.10 Mb
2003 10 31 010628.30 37.81 142.63 10 7.0 Mw
2003 11 1 131007.99 37.82 143.09 10 5.9 Mb
2003 11 9 192328.61 1.56 127.34 133 5.8 Mw
2003 11 11 184825.16 22.32 143.23 113 6.0 Mb
2003 11 12 002945.45 1.60 126.54 33 6.2 Ms
2003 11 12 082643.74 33.19 137.08 384 6.4 Mw
2003 11 14 184350.88 36.39 141.06 39 5.7 Mb
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

در اين پردازش، مدل ساده متشکل از يک نيم فضا با سرعت يك نواخت به عنوان مدل اوليةساختار سرعتي مورد استفاده قرار گرفته است. دادههاي پاشندگي سرعت گروه امواج ريليمورد استفاده در پردازش برگردان ادغام شده از رام و همكاران(٢٠٠٥) گرفته شده است.
تفاوت عمده الگوريتم محاسباتي در اين روش با روش تحليل توابع گيرنده در كنترل پروفيلمحاسبه اي براي ساختار سرعت پوسته براساس دادههاي پاشندگي است.

يافتهها
1268732149348Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017

Downloaded from jeg.khu.ac.ir at 11:39 IRST on Saturday October 28th 2017



قیمت: تومان

دسته بندی : زمین شناسی

دیدگاهتان را بنویسید