CC BY
59
Моделирование очагов динамических явлений на основе решения обратной задачи по геодезическим данным
Л.А. Назарова, Л.А. Назаров, М.П. Козлова1
Институт горного дела СО РАН, Новосибирск, 630091, Россия 1 Институт нефтегазовой геологии и геофизики СО РАН, Новосибирск, 630090, Россия
На основе решения обратной задачи для слоисто-неоднородной среды по геодезическим данным предложена методика определения параметров эквивалентного источника, описывающего аномальную зону (как прообраза готовящегося сейсмического события) в окрестности тектонического разлома.
Simulation of the foci of dynamic events on the basis of inverse problem solution using geodetic data
L.A. Nazarova, L.A. Nazarov, and M.P. Kozlova1
Institute of Mining SB RAS, Novosibirsk, 630091, Russia 1 Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia
An inverse problem is solved for the case of a layered heterogeneous medium using geodetic data. Based on the solution, we propose a technique for determining parameters of an equivalent source that describes an anomalous zone (as a prototype of the prepared seismic event) in the vicinity of a tectonic fault.
1. Введение
Быстрое развитие сети GPS, совершенствование программ обработки первичной информации и смещения пунктов земной поверхности (www.gpsworld.com, www.giscenter.ru) способствуют накоплению косвенных данных о геодинамических процессах в верхней части земной коры. В этой связи возникает задача интерпретации такой информации для определения характеристик деформационных процессов и их применения, в частности, для оценки состояния геосреды в зонах повышенного сейсмического риска.
Достигнутая к настоящему времени точность данных космической геодезии дает принципиальную возможность количественной оценки параметров очагов готовящихся сейсмических событий. В качестве моделей очагов рассматривают, как правило, зарождающийся или растущий разрыв [1] (по терминологии [2] — трещинное землетрясение), или аномальный участок нарушения сплошности (разлома) в породном массиве (верхней части земной коры) — тектоническое земле-
трясение [2]. Следуя концепции М.А. Садовского [3] о подавляющем преобладании числа землетрясений второго типа, а также приуроченности большинства сейсмических событий к разломным зонам [4], в настоящей статье рассмотрена модель очага тектонического типа.
Одним из методов определения напряженно-деформированного состояния среды, содержащей включение (занимающее некоторую область P), является построение фундаментального решения — функции Грина — с последующим интегрированием по P [5]. Идея работы — построить точечный источник, который «генерирует» в среде статическое поле (деформаций), аналогичное создаваемому аномальным участком нарушения. При этом известны только деформации на дневной поверхности («входная» информация), вычисленные, например, по данным GPS.
2. Постановка задачи
Рассмотрим прямоугольную область R = {0 < x < X, 0 < z < Z} в декартовой системе координат с нарушением
© Назарова ЛА., Назаров ЛА., Козлова М.П., 2GG8
сплошности D, наклоненным под углом в к вертикали (рис. 1). Всюду в R выполнены уравнения равновесия:
дахх даxz Л даxz даzz
—= 0, —= pg, (1)
дх dz dx dz
в R \ D среда упругая:
ахх = (Л + 2Ц)8хх + Л8zz .
а zz = Л8 хх + (Л + 2Ц )Е zz , (2)
ах2 = 2Ц8х2 ,
где а^- и 8jj = 0.5(ыг-,j + Ujj) — компоненты тензоров напряжений и деформаций; щ — смещения (i, j = x, z); Л и ц — параметры Ламе; p — плотность; g—ускорение свободного падения.
Материал нарушения описывается простыми уравнениями состояния [6]:
ап = K n Aun’ а t = Kt A«t> (3)
где аn и аt — нормальное и касательное напряжения на D; Kn и Kt — соответствующие жесткости; u п = = их cos в + uz sin в, ut =-ых sin в + uz cos в (рис. 1); знак А означает разницу смещений на берегах нарушения.
На dR зададим следующие граничные условия:
а xz (х ,0) = 0, а zz (х ,0) = 0, а^ (х, Z) = 0, uz (х, Z) = 0,
ых (0, z) = 0, аxz (0, z) = 0, (4)
ахх (х. z) = qpgz, аxz(X> z) = °, где q — коэффициент бокового распора. Заметим, что значения q < 1 соответствуют сбросовому тектоническому режиму, q > 1 — взбросовому [7].
Физические свойства области R, аналогом которой выбрано вертикальное сечение земной коры вкрест простирания основных разломных структур Байкальской рифтовой зоны, приведены в табл. 1 [8].
Параметры Ламе вычислялись по скоростям продольных VP и поперечных VS волн:
^ = pFs2, Л = pVp2 -2ц.
Жесткости Kn и Kt оценивались по формулам [9]:
E K =ц
—, Kt =—, h h где h — мощность нарушения (рис. 1).
K
(5)
c>xz = C>zz = 0 Нарушение сплошности X
х,
Участок j||ut F/ N О)
"зацепления" Q. Q-
zc берегов ^ f 7 II
-jmr Y N
О
о
1Ш м
i щ 1 М й
| Ц Вмещающая хс; Ц среда ю
м N N N : 0 Gxz = 0
3. Метод решения
Пусть в силу тех или иных причин (неоднородность контактирующих поверхностей, локальное повышение нормального напряжения) на нарушении D возник участок «зацепления» берегов P с повышенными деформационными и/или прочностными свойствами. Заметим, что такой подход к описанию очага соответствует модели подготовки землетрясения типа «stick-slip» [10]. Будем моделировать это явление увеличения жесткостей Kn и Kt в ^-1 раз (5 < 1).
В результате геомеханические поля в среде изменятся и, в частности, поле деформаций получит приращение А8j (х, z), горизонтальную составляющую которого у(х) = А8хх(х, 0) на свободной поверхности можно вычислить по данным GPS. Восстановить параметры аномальной зоны (местоположение, форму и размеры, физические свойства) по известной функции S достаточно проблематично. Поэтому применим обычный для сейсмологии подход [11]: построим эквивалентный точечный «источник», который бы создавал на поверхности z = 0 дополнительные горизонтальные деформации, близкие к у. Таким образом, необходимо подобрать тип источника, а также найти его параметры (рис. 1): амплитуду F и ориентацию (угол а) силы, координаты (хс, zс) точки ее приложения.
Задача (1)-(5) решалась методом конечных элементов с использованием программы МКЭЧК [12] при X= = 200 км, Z = 60 км, шаг дискретизации 0.25 км.
Были рассмотрены несколько типов сосредоточенных источников: одиночная сила, двойная сила без момента и с моментом. Численные эксперименты показали, что наиболее подходит последний тип, причем силы должны прикладываться на разных берегах нарушения (рис. 1), при этом момент равен Fh sin в. На рис. 2 для сравнения показано распределение приращения деформаций А8 хх на поверхности z = 0 при возникновении на нарушении аномальной зоны (штриховая линия) и при действии точечного источника типа «двойная сила с моментом» (сплошная линия). Можно видеть, что кривые имеют весьма схожую форму. Здесь и далее для краткости изложения значения параметров расчета приведены непосредственно на рисунках.
Таблица 1
Физические свойства расчетной области
Рис. 1. Постановка задачи и схема эквивалентного точечного источника
z, км VS, м/с /с м p, кг/м3
1 0-2 1000 2000 2500
2 2-6 3 500 5 800 2700
3 6-10 1900 3 500 2550
4 10-15 3 400 6100 2800
6 15-35 3 700 6400 2800
7 35-60 4 300 7 800 2800
40 "
20
-20 -
-40
-60
_ Участок нарушения
со "сцепленными" берегами
_ Заглубленный точечный
источник "двойная сила
- с моментом"
1\ DL = 0
А /' - д it \\ It /
\\ It \\ /| " '\/| \ 1 // /1 , // V 1
\ / - \ 1 1 1
V 1 1
/ = 3 км
zc = 10 км
1 1 1 1 1 = 0.25 i i i i
50
70
90 110
х, км
130
150
Рис. 2. Сравнение распределений приращений горизонтальных деформаций на дневной поверхности при выборе типа эквивалентного источника
Теперь введем целевую функцию
X 2
G(FX, Fz, Xc, zc) = J [Де XX (X, 0) - V (X )]2dx,
XI
где FX = F cos a, Fz = F sin a — компоненты силы; [ Xl5 x2] — участок свободной поверхности, где известны горизонтальные деформации. Точка минимума G до-
ф ф s|c ф
ставляет решение обратной задачи (FX , Fz , xc , zc ).
В [13] разработан эффективный алгоритм минимизации подобных функций, основанный на их линейнос-
* *
ти по FX и Fz. Компоненты силы FX и Fz определяются методом наименьших квадратов, число неизвестных аргументов сокращается до двух, и необходимо минимизировать функцию Ф(Xc, zc) = G(FX , Fz , Xc, zc). Численный анализ при варьировании параметров задачи в широком диапазоне (0.125 < £ < 0.5, 0.4 < q < 3, 1 < l < < 3 км, l — вертикальный размер аномальной зоны) показал, что Ф имеет единственный минимум (рис. 3), который может быть найден методом градиентного спуска [14].
<1 -1
' /I q = 0.4 zc = 10 км / = 3 км . 5 = 0.5 ^
1 1 1 ■! F Т 1 1 1 1
50
100
х, км
150
200
Рис. 3. Изолинии минимизируемого функционала Ф(Xc, zc)
Рис. 4. Примеры подбора параметров эквивалентного источника типа «двойная сила». Относительная ошибка аппроксимации входных данных: 5 = 82 (а); 3 % (б)
4. Обсуждение результатов
На рис. 4 в качестве иллюстрации приведены графики «входной» информации (функция у, штриховые линии) и дополнительных горизонтальных деформаций Де хх от точечного источника (серые линии) для двух «крайних» случаев: относительная точность 5 аппроксимации у изменяется от 82 до 3 %. Хотя следует отметить, что даже для невысоких 5 качественное поведение кривых идентично. Расчеты показали, что с увеличением I величина 5 уменьшается.
На рис. 5 для различных значений коэффициента бокового распора q приведено распределение амплитуды F и угла а в зависимости от величины £, которую можно ассоциировать с «контрастностью» аномальной зоны на фоне нарушения: чем меньше £, тем сильнее различаются их свойства. Как и следовало ожидать, увеличение £ ведет к уменьшению F, но практически не влияет (за исключением малых q) на ориентацию силы. В то же время, значение а быстро уменьшается от 55° при q = 0.4 до 5° при q = 3 (рис. 5, б).
Рис. 5. Зависимость параметров эквивалентного точечного источника от £ при различных коэффициентах бокового распора q
Рис. 6. Зависимость амплитуды эквивалентного точечного источника от коэффициента бокового распора q при различных значениях £: I = = 1 (а); 2 (б); 3 км (в)
На рис. 6 представлены зависимости F(q) для трех значений I при различных £. При фиксированных £ здесь наблюдается практически линейный рост амплитуды силы при увеличении размеров аномальной зоны и коэффициента бокового распора.
5. Заключение
Предложенная методика установления количественных значений параметров аномальных зон по данным измерений смещений на свободной поверхности позволяет оценить (привлекая дополнительные сведения о сейсмическом КПД [15]) магнитуду и определить фокальный механизм готовящегося сейсмического события.
Работа поддержана РФФИ (проект № 07-05-01020), Интеграционным проектом СО РАН и ДВО РАН № 87.
Литература
1. Костров Б.В. Механика очага тектонического землетрясения. -М.: Наука, 1975. - 234 с.
2. Шемякин Е.И., Курленя М.В., Кулаков Г.И. К вопросу о классификации горных ударов // ФТПРПИ. - 1986. - № 6. - С. 3-11.
3. Садовский М.А., Болховитинов Л.Г., Писаренко В.Ф. Деформирование среды и сейсмический процесс. - М.: Наука, 1987. - 101 с.
4. Шейдеггер А. Основы геодинамики. - М.: Недра, 1987. - 384 с.
5. Новацкий В. Теория упругости. - М.: Мир, 1975. - 872 с.
6. Barton N.R. Deformation phenomena in jointed rock // Geotechnique. -
1986. - V. 36. - No. 2. - Р. 147-167.
7. Назарова Л.А. Использование сейсмотектонических данных для оценки полей напряжений и деформаций земной коры // ФТПРПИ. - 1999. - № 1. - С. 28-36.
8. Дядъков П.Г., Назаров Л.А., Назарова Л.А. Детальное моделирование поля напряжений земной коры и динамической неустойчивости сейсмоактивных разломов при рифтогенезе // Геология и геофизика. - 1997. - № 12. - С. 2001-2010.
9. Юфин C.A. Механические процессы в природном массиве и их взаимодействие с подземными сооружениями / Автореф. дис. ... докт. техн. наук. - М.: МГИ, 1991. - 35 с.
10. Мячкин В.И. Процессы подготовки землетрясений. - М.: Наука, 1978. - 232 с.
11. Райс Дж. Механика очага землетрясения. - М.: Мир, 1982. -215 с.
12. Назарова Л.А. Напряженное состояние наклонно-слоистого массива горных пород вокруг выработки // ФТПРПИ. - 1985. - № 2. -С. 33-37.
13. Голъдин С.В., Назаров Л.А., Назарова Л.А., Козлова М.П. Оценка параметров очага готовящегося сейсмического события по данным о деформациях свободной поверхности // ФТПРПИ. - 2007. -№ 3. - С. 25-35.
14. Бахвалов Н.С. Численные методы. Т. 1. - М.: Наука, 1975. - 632 с.
15. Доброволъский И.П. Теория подготовки тектонического землетря-сения.-М.: Наука, 1991.-219 с.
Поступила в редакцию 29.10.2007 г.
