CC BY
44
УДК 539.3
Оценка фокальных параметров предстоящего землетрясения по деформациям дневной поверхности
Л.А. Назарова, Л.А. Назаров, П.Г. Дядьков1, М.П. Козлова1,
В.П. Ильин2, Я.Л. Гурьева2
Институт горного дела СО РАН, Новосибирск, 630091, Россия 1 Институт нефтегазовой геологии и геофизики им. А.А. Трофимука СО РАН, Новосибирск, 630090, Россия 2 Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, 630090, Россия
Предложена модель очаговой зоны готовящегося сейсмического события в виде зоны на тектоническом нарушении, фрактальная размерность f которого больше, чем у остальной части разлома. Показано, что по деформациям свободной поверхности можно определить величину f и глубину этой зоны.
Ключевые слова: массив горных пород, землетрясение, очаг, нарушение сплошности, аномальная зона, фрактальная размерность, метод конечных элементов
Estimation of focal parameters of a forthcoming earthquake from day surface deformation
L.A. Nazarova, L.A. Nazarov, P.G. Dyadkov1, M.P. Kozlova1, V.P. Ilin2 and Ya.L. Gurieva2
Institute of Mining SB RAS, Novosibirsk, 630091, Russia 1 Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia 2 Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk, 630090, Russia
A model is proposed for the focal zone of a forthcoming seismic event in the form of a zone on a tectonic discontinuity whose fractal dimension f is larger than that of the rest part of the fault. It is shown that the fractal deformation and the depth of the zone can be estimated from deformation of a free surface.
Keywords: rock massif, earthquake, focus, discontinuity, anomalous zone, fractal dimension, finite element method
1. Введение
В настоящее время прогноз динамических событий различного масштабного уровня (от микроударов в породных массивах при проведении горных работ до землетрясений) основан на статистической обработке временных рядов сейсмических событий с использованием различных методик [1-4]. Результатами таких исследований являются, как правило, эмпирические соотношения между энергией (магнитудой) событий, их количеством, временем появления аномалий и размерами очага. Параметры соотношений могут меняться в зависимости от региона [5, 6]. Более детальную информацию о характеристиках готовящегося сейсмического события (время, местоположение очага и т.д.) с помощью
статистического анализа получить достаточно сложно, поэтому используется модельный подход, основанный на различных гипотезах о механизме подготовки землетрясений [1, 7].
Очаги подавляющего большинства землетрясений приурочены к разломным зонам литосферы [8, 9], поэтому гипотеза о том, что подготовка динамического события тектонического и «неотектонического» типа [10] происходит именно на нарушении сплошности, в достаточной мере обоснована. Одна из моделей, описывающих механизм формирования очаговой зоны, предполагает наличие на разломе аномальной зоны С, в окрестности которой под действием квазистатического изменения внешних усилий возникает область концент-
© Назарова Л.А., Назаров Л.А., Дядьков П.Г., Козлова М.П., Ильин В.П., Гурьева Я.Л., 2010
рации напряжений и при достижении ими предельного уровня происходит подвижка берегов с выделением части накопленной энергии в виде сейсмических волн [1113]. Дополнительные напряжения, обусловленные зоной С, можно зарегистрировать инструментально и попытаться оценить параметры С на основе этой информации.
Такой подход был реализован, например, в [14] для прямолинейного наклонного нарушения, когда зона С моделировалась участком с повышенной жесткостью, а «действие» С — эквивалентным источником типа «двойная сила с моментом». В настоящей работе в качестве аномальной зоны рассматривается участок разлома с иной фрактальной размерностью.
Фрактальная размерность линии L определяется формулой [15] f = - 1п И/ 1п I, где N— число сегментов длины I, «укладывающихся» на L; если L — прямая, то /' = 1. Фрактальная размерность оказалась удобным параметром, характеризующим одновременно и «изрезан-ность» линии, и степень дискретизации при описании геометрии последней. В [15] предложен метод генерации линий с заданной фрактальной размерностью, которые реализованы в [16, 17] с использованием различных законов случайного распределения угла наклона сегментов.
Цель настоящей работы — создать способ количественной оценки параметров зоны С по изменению поля смещений, замеренному на свободной поверхности.
2. Постановка задачи
В прямоугольной области G с размерами Хх У расположено нарушение сплошности Р (мощность А, угол падения в), участок которого на глубине h имеет фрактальную размерность /' > 1, остальная часть — прямолинейная (рис. 1). Ось у декартовой системы координат направим вертикально вниз, ось х — горизонтально. Рассмотрим плоское деформированное состояние, позволяющее описывать сбросовый и взбросовый тектонические режимы с помощью плоских моделей [18].
(4)
В G выполнены уравнения равновесия
3р, р +РЯ8г/ = 0, (1)
где <3^ — компоненты тензора напряжений (;', j = х, у); р — плотность; g — ускорение свободного падения; 8р — символ Кронекера; по повторяющимся индексам производится суммирование. Вмещающая среда (область G \ Р) — упругая, ее деформирование описывается законом Гука:
3ц = (Х8р + 2ц)Ер = 0, (2)
где Ер = 0.5(мг-, р + и л) — компоненты тензора деформаций, и{ — смещения; X, ц — параметры Ламе. Уравнения состояния Р примем в упрощенном виде [19]:
3 = КпР, т = К(R, (3)
3 и т — нормальное и касательное напряжения вдоль Р; Кп и К( — соответствующие жесткости; Р и R — конвергенция и проскальзывание берегов Р. На 8G сформулируем следующие условия: их (0, у) = 0, 3ху(0, у) = 0,
3хх (х> у) = Ч3У (у)> 3ху (х, у) = 0 3уу (х,0) = 0, 3ху (х, 0) = 0,
иу (х, У) = 0, 3 ху (х, У) = 0,
у — коэффициент бокового отпора; 3Г = рgy — литостатическое напряжение.
Примем характерное для многих сейсмоактивных регионов горизонтально-слоистое строение среды, прообразом которого выбрана окрестность гипоцентра Алтайского землетрясения 2003 г. Физические свойства пород [20] приведены в табл. 1.
Нормальная Кп и касательная К( жесткости вычислялись по предложенным в [21] соотношениям Кп = = Е/А, К( = ц/А, где Е — модуль Юнга вмещающих пород.
3. Тестирование численного алгоритма
Расчеты проводились методом конечных элементов с использованием кода 2МКЭЧК [22]. Принимая во внимание, что основная система линейных уравнений имела большую размерность (число неизвестных М ~ ~ 105-106), было выполнено тестирование алгоритма на задаче о напряженно-деформированном состоянии однородной упругой полуплоскости в поле гравитаци-
Таблица 1
Физические свойства расчетной области
Рис. 1. Схема расчетной области и граничные условия
у, км X, ГПа |Х, ГПа р, кг/м3
1 0-3 34 28 2 200
2 3-6 37 30 2 500
3 0 1 6 40 33 2 700
4 10-15 42 35 2 800
5 15-20 46 42 3 000
онных и тектонических сил [23], в которой поля смеще-
- 0 - „0 ний щ и напряжений а^ определяются аналитически:
u°( x y) =
u°(x, y) = 0.5(Ayy2 - Axx2),
(x y) = qav(y)> (5)
a<Xy(x y) = °
(x y) = av(y),
где Ax =pg(q-v)/E; Ay = pg(1 -qv)/E, v — коэффициент Пуассона. Решение (5) использовалось для формулировки краевых условий в тестовом примере:
- на свободной поверхности
аxy (x, 0) = аyy (x, 0) = 0, (6)
- на остальных участках границы G
ux be = u0, uy be = u0 • (7)
Тестовые расчеты проводились на базе Сибирского суперкомпьютерного центра в ИВМиМГ СО РАН на сервере Intel Itanium 2. Оказалось, что поле смещений, найденное из решения системы (1)-(3), (6) и (7), совпадает с (5), а напряжения незначительно отличаются в окрестности границ вследствие их осреднения по элементам в численном алгоритме. В табл. 2 приведено сравнение времен расчета, когда для решения системы линейных уравнений использовались метод Гаусса и явный итерационный метод неполной факторизации с применением модификации Айзенштадта и ускорением по методу сопряженных градиентов [24] (строки 1 и 2 соответственно). Сравнительное тестирование проводилось на сетках с увеличивающимися размерностями. Качество итерационного решения оценивалось по максимальной абсолютной разности элементов вектора решения, полученного прямым и итерационным методами. Для всех указанных расчетов это максимальное отклонение оказывалось не больше величины порядка 10-11. Тестирование показало, что на всех сетках итерационный алгоритм решения системы линейных уравнений работает быстрее, причем его преимущество увеличивается с ростом размерности системы; кроме того, начиная с некоторого M, для реализации прямого метода не хватает оперативной памяти. Кроме того, использование современных ресурсов Сибирского суперкомпью-терного центра представляется важным аспектом для получения численных решений системы линейных уравнений больших размерностей, т.к. ограничения на размер оперативной памяти и время получения решения
Таблица 2
Время расчета для различной размерности системы линейных уравнений
M 28 210 212 214 216 218 220
Время 1 0.005 0.059 0.81 12.2 190 4121 76 164
расчета, с 2 0.002 0.016 0.12 1.0 11 402 702
относятся к критическим моментам.
4. Параметрический анализ поля деформаций на свободной поверхности
Проанализируем поле деформаций на дневной поверхности при изменении параметров задачи. Из (2) и (4) следует, что при у = 0 v
е xy = 0 е xx = 1 _ v е yy,
поэтому остановимся только на горизонтальной компоненте тензора деформаций еxx, которую можно рассчитать по смещениям на свободной поверхности, полученным по данным космической геодезии.
Будем полагать, что возможные возмущения поля деформаций в рассматриваемой области вызваны наличием аномальной зоны C, поэтому анализировалась величина Деxx = е xx -е xxs, где последний член соответствует деформации в случае прямолинейного нарушения. Такой подход позволяет не только «ликвидировать» неопределенность, связанную с предположением прямолинейности участка F \ C, но и использовать данные GPS, по которым можно вычислить не абсолютные деформации на поверхности, а их приращения за некоторый промежуток времени.
Варьировались следующие параметры: глубина h, угол падения в и фрактальная размерность f аномальной зоны. Расчеты проводились при X = 30 км, Y = 20 км, s = 0.5 км, Д = 0.1 км, q = 0.4-2.0.
Замечание. Точность современных GPS-данных позволяет надежно определить деформации до 10-8 [25], поэтому диапазон варьируемых параметров выбирался так, чтобы абсолютная величина Деxx превышала указанную нижнюю границу.
На рис. 2 для различных значенийf приведен пример профилей C, сгенерированных при условии равномерного распределения угла отклонения сегмента на отрезке [-8, 8], где (как показано в [17])
8 = д/ 3sh[(2 - 2/ f )ln N ]/N.
Можно видеть, что возрастание f ведет к увеличению «изрезанности» линий.
На рис. 3 изображен фрагмент сетки конечных элементов в окрестности аномальной зоны с f = 1.01 (штриховые линии показывают берега прямолинейного разлома, для удобства восприятия показан только каждый четвертый узел).
. - — ____f=1.005
;=^СГ= uioT
f = гозх^-^—^
f=1.N
0.0 0.2 0.4 0.6 0.8 1.0
Рис. 2. Профили аномальной зоны при различных аномальных размерностях
Рис. 3. Фрагмент сетки конечных элементов в окрестности аномальной зоны
Распределение Аехх на свободной поверхности при различных в для h = 5 км, f = 1.1 и у = 1.5 показано на рис. 4: амплитуда приращения деформаций практически не зависит от угла падения нарушения, меняется только пространственное положение экстремумов.
На рис. 5 представлены графики Аехх (х, 0) для различных глубин h расположения аномальной зоны при /= 1.1, у = 2.0 и в = 60°. Как и ожидалось, амплитуда приращения горизонтальных деформаций затухает с ростом h примерно по квадратичному закону. Следует обратить внимание, что при прочих одинаковых условиях распределение Аехх практически не зависит от конфигурации С (которая генерируется случайным образом), а только от величины/. Это следствие принципа Сен-Венана, интерпретируемого с отличных от общепринятых позиций [26]. Аналогичный вывод можно сделать и из рис. 6, на котором показано распределение Ае хх на поверхности у = 0 для различных f при h = = 5 км, у = 0.8 и в = 70°. Кроме этого, здесь видно, что с ростом фрактальной размерности зоны С амплитуда приращений горизонтальных деформаций увеличивается.
В силу (4) величина Аехх пропорциональна коэффициенту бокового отпора у, количественную оценку которого можно осуществить по сейсмологическим данным [18].
Таким образом, одними из основных параметров, контролирующих изменение приращений деформаций на свободной поверхности, являются глубина h и фрак-
10 12 14 16 18 20
X, км
Рис. 4. Распределение приращений горизонтальных деформаций на свободной поверхности при различных углах падения нарушения
-3 -I-----------------------I--------------------.---------------------:---------------------:--------------------
10 12 14 16 18 20
X, км
Рис. 5. Распределение приращений горизонтальных деформаций на свободной поверхности при различной глубине расположения аномальной зоны
тальная размерность f аномальной зоны, поэтому введем функцию D(х, f, К) = Аехх (х, 0, f, К).
5. Решение обратной задачи
Рассмотрим целевую функцию
х2
ф(/, К) = |(Б(x, /, К -О0(x))2dx,
х1
где х1, х2 е [0, X]; D0 — входные данные на свободной поверхности. В расчетах в качестве последних использовались синтетические деформации (решение системы (1)-(4), полученное при некоторых значениях / = /0 и К = К0), на которые налагался мультипликативный шум:
^(х) = (1 + £)D(х, /0, и
где £ — случайная величина, равномерно распределенная на отрезке [-у, у].
На рис. 7 показаны линии уровня Ф (нормирован на максимальное значение) в двух диапазонах изменения фрактальной размерности: fе [1.03, 1.08] (/0 = 1.06, у = 0.1, рис. 7, а) и f е [1.003, 1.009] (/0 = 1.006, у = = 0.05, рис. 7, б). В обоих случаях К0 = 5 км, х1 = 12 км и х2 = 18 км. Можно видеть, что при / > /* (/* = 1.01) функция Ф имеет единственный минимум (который можно найти, например, методом градиентного спуска [27]) — обратная задача однозначно разрешима.
Иная ситуация имеет место при / < /* (рис. 7, б). Из-за внесенной даже небольшой погрешности (у = = 0.05) возникают «побочные» экстремумы, поэтому
10 12 14 16 18 20
X, км
Рис. 6. Распределение приращений горизонтальных деформаций на свободной поверхности при различной фрактальной размерности аномальной зоны
Рис. 7. Изолинии целевой функции Ф(к, f)
при использовании градиентных методов необходимо выполнять несколько циклов поиска с различными начальными приближениями, выбирая среди найденных локальных минимумов самый «глубокий».
Значение /* определено на основе численных экспериментов, оно имеет тенденцию к возрастанию при увеличении h и у.
6. Заключение
Предложенный подход позволяет при известных структуре и свойствах среды (которые можно установить с использованием геологических и геофизических методов), величине естественных горизонтальных напряжений (оцениваемых по сейсмологическим данным о фокальных механизмах сейсмических событий и их магнитуде) определить параметры расположенной на тектоническом нарушении аномальной зоны, ассоциируемой с очагом готовящегося горного удара или землетрясения.
Работа выполнена в рамках проектов РФФИ (№№ 09-05-00975, 10-05-01042) и интеграционного проекта СО РАН № 44.
Литература
1. СоболевГ.А. Основы прогноза землетрясений. - М.: Наука, 1993. -313 с.
2. Опарин В.Н., Тапсиев А.П., Востриков В.И. и др. О возможных причинах увеличения сейсмической активности шахтных полей рудников «Октябрьский» и «Таймырский» Норильского месторождения в 2003 г. Часть 1. Сейсмический режим // ФТПРПИ. -2004.- № 4. - С. 3-22.
3. Дядьков П.Г., Кузнецова Ю.М. Аномалии сейсмического режима перед сильными землетрясениями Алтая // Физ. мезомех. - 2008. -Т. 11. - № 1. - С. 19-25.
4. Физические основы прогнозирования разрушения горных пород при землетрясениях / Под ред. М.А. Садовского, Г.А. Соболева. -М.: Наука, 1987. - 127 с.
5. Завьялов АД. Среднесрочный прогноз землетрясений: основы, методика, реализация. - М.: Наука, 2006. - 255 с.
6. Востриков Г.А. Связь параметров графика повторяемости, сейсми-
ческого течения и очага землетрясения // Труды Геологического института РАН, Вып. 482. - М.: Изд-во Геологического института РАН, 1994. - 292 с.
7. Соболев Г.А., Пономарев А.В. Физика землетрясений и предвестники. - М.: Наука, 2003. - 270 с.
8. Буллен К. Введение в теоретическую сейсмологию. - М.: Мир, 1966. - 460 с.
9. Шейдеггер А. Основы геодинамики. - М.: Недра, 1987. - 384 с.
10. Шемякин Е.И, Курленя М.В., Кулаков Г.И. К вопросу о классификации горных ударов // ФТПРПИ. - 1986. - № 5. - С. 3-11.
11. Brace W.F., Byerlee J.D. Stick-slip as a mechanism for earthquakes // Science. - 1966. - V. 153. - P. 990-992.
12. Scholz C.H. The Mechanics of Earthquakes and Faulting. - Cambridge: Cambridge University Press, 1990. - 439 p.
13. Дядьков П.Г., Назаров Л.А., Назарова Л.А. Моделирование напряженного состояния земной коры в окрестности сейсмогенного разлома в центральной части Байкальского рифта // Геология и геофизика. - 1996. -№9.- С. 71-78.
14. Назарова Л.А., Назаров Л.А., Козлова М.П. Моделирование очагов динамических явлений на основе решения обратной задачи по геодезическим данным // Физ. мезомех. - 2008. - Т. 11. - № 1.-С. 51-54.
15. Mandelbrot B.B. The Fractal Geometry of Nature. - San Francisco: Freeman, 1983. - 460 p.
16. Seidel J.P., Haberfield C.M. Toward an understanding of joint roughness // Rock Mech. Rock Eng. - 1995. - V. 28. - No. 2. - P. 69-92.
17. Назаров Л.А., Назарова Л.А. О связи деформационных свойств нарушений сплошности в породном массиве и ее фрактальной размерности // ФТПРПИ. - 2008. - № 5. - С. 3-13.
18. Назарова Л.А. Использование сейсмотектонических данных для оценки полей напряжений и деформаций земной коры // ФТПРПИ. - 1999. - № 1. - С. 28-36.
19. Barton N.R. Deformation phenomena in jointed rock // Geotechnique. - 1986. - V. 36. - No. 2. - P. 147-168.
20. ЕмановА.А., Лескова Е.В. Особенности строения эпицентральной зоны Чуйского (Горный Алтай) землетрясения по данным метода сейсмической томографии с двойными разностями // Матер. II Межд. симп. «Активный геофизический мониторинг литосферы Земли». - Новосибирск: Изд-во СО РАН, 2005.
21. Юфин C.A. Механические процессы в природном массиве и их взаимодействие с подземными сооружениями: Автореф. дис.... докт. техн. наук. - М.: МГИ, 1991. - 35 с.
22. Назарова Л.А. Напряженное состояние наклонно-слоистого массива горных пород вокруг выработки // ФТПРПИ. - 1985. - № 2. -
С. 33-37.
23. Мусхелишвили Н.И. Некоторые основные задачи математической теории упругости. - М.: Наука, 1966. - 707 с.
24. Ильин В.П. Методы неполной факторизации для решения алгебраических систем. - М.: Физматлит, 1995. - 288 с.
25. http://www.gpsworld.com.
26. Новацкий В. Теория упругости. - М.: Мир, 1975. - 872 с.
27. Бахвалов Н.С. Численные методы. - М.: Наука, 1973. - Т.1. -
632 с.
Поступила в редакцию 15.04.2010 г.
Сведения об авторах
Назарова Лариса Алексеевна, д.ф.-м.н., гнс ИГД СО РАН, larisa@misd.nsc.ru Назаров Леонид Анатольевич, д.ф.-м.н., зав. лаб. ИГД СО РАН, naz@misd.nsc.ru Дядьков Петр Георгиевич, к.г.-м.н., доц., зав. лаб. ИНГГ СО РАН, dyadkovpg@ipgg.nsc.ru Козлова Марина Петровна, мнс ИНГГ СО РАН, kozlovamp@ipgg.nsc.ru Ильин Валерий Павлович, д.ф.-м.н., гнс ИВМиМГ СО РАН, ilin@lapasrv.sscc.ru Гурьева Яна Леонидовна, к.ф.-м.н., снс ИВМиМГ СО РАН, yana@lapasrv.sscc.ru
