CC BY
17Учёные записки Крымского федерального университета имени В. И. Вернадского. География. Геология. Том 2 (68). № 4. 2016 г. С. 249-258.
УДК 530.344
РАСЧЕТ ВРЕМЕНИ ПРОБЕГА СЕЙСМИЧЕСКИХ ВОЛН ПРИ ЛИНЕЙНОМ ЗАКОНЕ ВОЗРАСТАНИЯ СКОРОСТИ И НАЛИЧИИ ГРАНИЦЫ РАЗРЫВА
СКОРОСТЕЙ Пивоварова Н. Б.
Государственное автономное учреждение Республики Крым «Крымский экспертный совет по оценке сейсмической опасности и прогнозу землетрясений» ГАУ КР «КЭС», г. Симферополь, Россия.
E-mail: pivnb@yandex.ru
Приведены формулы расчета времени пробега сейсмических волн от эпицентра до станции при условии линейного закона скоростей волн выше границы разрыва скоростей и постоянной скорости ниже разрыва. Рассмотрена задача нахождения координат очага землетрясения на основе времен регистрации сейсмических волн на нескольких станциях данного региона. В расчете используются одновременно времена прихода на станцию S- и Р-волн. Задача определения координат очага решается на основе минимизации суммарной невязки. Минимизация проводится с помощью соответствующей программы пакета Excel.
Ключевые слова: землетрясение, координаты очагов, поля скоростей. ВВЕДЕНИЕ
Спектр задач, связанных с распространением в толще Земли сейсмических сигналов, вызванных землетрясениями, может быть разделен на несколько последовательных задач: при заданном профиле скоростей и известных временах прихода сигнала на несколько станций определить координаты гипоцентра очага, по набору уточненных координат очагов уточнить параметры скоростного разреза. В данной работе ограничимся первой задачей нахождения координат очагов. Предполагаем, что скорость зависит только от вертикальной координаты z на масштабах, для которых влияние кривизны земной поверхности пренебрежимо мало, станции будем считать расположенными на поверхности Земли, а саму поверхность плоской.
1. ВОЗМОЖНЫЕ ТИПЫ ТРАЕКТОРИЙ
Используем стандартные обозначения: ось z направлена вниз, e - угол между вертикалью и направлением распространения сигнала, v(z) - скорость
sin e
распространения сигнала, p =--лучевой параметр.
v
Известно, что при плоском профиле скоростей p = const на траектории сигнала, поэтому:
d- = vz = vcose = ±v(z)^/l - p2v(z)2 ,
где знак «+» соответствует распространению сигнала вниз, а «-» - вверх.
249
Таким образом, вычисление горизонтального расстояния от эпицентра до станции (Ь) и времени прохождения сигнала (/) сводятся к вычислению интегралов:
Ь = | + ^ РУ({= |- + йт
л/1" Р М х)2 -1 - р М т)2
При этом интегрирование ведется как по нисходящей, так и по восходящей части траектории, при этом знак выбирается, разумеется, таким образом, чтобы значение интеграла в обоих случаях было положительным:
йх ру( т)
= tge = + у у '
~ 16 ^ ~ — I-
йт ф - р2у(т)2
Предполагая, что скорость сейсмической волны монотонно возрастает до некоторой глубины И&, на которой имеется зона скачкообразного возрастания скорости до величины \ъг, а ниже этой зоны скорость остается практически постоянной, можно представить себе несколько типов траекторий сигнала в зависимости от взаимного расположения гипоцентра и станции:
A) восходящая траектория с гипоцентром выше разрыва - сигнал все время распространяется вверх;
Б) нисходяще-восходящая траектория - гипоцентр расположен выше разрыва, сигнал распространяется вначале вниз, затем - вверх, постоянно оставаясь выше зоны разрыва;
B) «головная волна» - сигнал распространяется вниз до зоны разрыва, к которой подходит с таким значением угла е, что дальнейшее распространение происходит горизонтально по зоне разрыва, далее сигнал выходит из зоны вверх с тем же значением е;
Г) восходящая траектория с гипоцентром ниже зоны разрыва - сигнал распространяется вверх по прямой до зоны разрыва, где скачком меняет направление, далее распространяется вверх.
Вообще говоря, станция может принимать сигналы, распространяющиеся по траекториям разных типов, поэтому важным является вопрос - какие типы траекторий возможны при различных положениях гипоцентра и по какой траектории сигнал придет первым. В дальнейшем не будем рассматривать волны с эпицентром ниже зоны разрыва.
2. ФОРМУЛЫ РАСЧЕТА ВРЕМЕНИ ПРОХОЖДЕНИЯ СИГНАЛА
Для выбранных типов траекторий общие выражения для расстояния от эпицентра до станции и времени прохождения сигнала принимают следующий вид:
250
ру( т)йт _И йт
А) Ь ==\^Ъ t = |
0 л/1 - РЧт)2 0 у(т)^ - р2у(т)2 '
Б) Ь V ру( т)йт V ру(т)йт
Б) Ь = -1л „2,^2 + I
д/1 - рМт)2 И л/1 - рМт)2
Г йт Г йт
t = - --" + I
* у(т)^гр2у(т)2 И у(т)^1-р2у(т)2 '
где: [(у(И)]тах = —, то есть И тах(р - максимальная глубина траектории сигнала, р
при этом Итах < Иьг. Это означает, что вся траектория лежит выше зоны разрыва.
ИЬг / \ л ИЬг
В) Ь --Ь^ЙтЬ +
Ьг ру( т)йт + Ь + - ру( т)йт -0 л/1 - р2у(т)2 * - л/1 - р2у(т)2 ,
Ьг йт Ь, ^ йт
t = I -, "" + +
0 у(т- р2у(т)2 уЬг И у(т^1 - р2у(т)2
где: р = 1/уьг, ¿ьг - длина части траектории, проходящей вдоль разрыва. Естественным условием является Ььг>0.
Во всех вариантах решение заключается в определении лучевого параметра р из первого уравнения и подстановка его во второе, в случае же головной волны таким параметром является длина траектории вдоль разрыва Ььг. В общем случае эта задача допускает численное решение, в частных случаях возможно также аналитическое решение.
3. ФОРМУЛЫ РАСЧЕТА ВРЕМЕНИ ПРОБЕГА ПРИ ЛИНЕЙНОЙ МОДЕЛИ СКОРОСТИ
Рассмотрим простейшую модель профиля скорости, то есть линейную зависимость возрастания скорости с глубиной, при этом на глубине Н =35 км находится высокоскоростная граница. Будем рассматривать только очаги землетрясений, находящиеся на глубинах 2 < 35 км.
Пусть: у(т) = у0+ат. В этом случае интегралы вычисляются точно:
г йт -1 / 1 ^
I-, = = — агееп -
1 у{т)^1 - р 2у(т)2 « I ру( т) 251
I
J-
pv( z)dz
- Рыл2
= —Vl - p2v(z)2 .
y/l - p2v(z)2 ap Выражения для величин L и t приобретают тогда следующий вид:
aL = -,11/p2 - v02 -Ji/J2
A)
-2 - vl
at= arcch
f 1 ^
pv0
( 1 V
- arcch
pvH
aL = yj 1/p2 - v02 1/p
Б)
at = arcch
f 1 ^
pvo
22 2 - v2
VH f 1 ^,
+arcch
v pvh у 2Vv2 - v2(Hbr)
aL = aLjr + ^vl - vH - 2лhi - v2
aL,
at = —— + arcch
vbr + arcch vbr - 2 arcch vbr
V v(Hbr J
V v0 У V vH J
В)
Заметим, что в этих формулах Уи = у (И = Уо+аИ - скорость сигнала в эпицентре,у(Иьг) = Уо+аИьг - скорость сигнала на верхней границе зоны разрыва, а Уьг - на нижней ее границе, то есть предполагается Уьг>у(Иьг).
Первое уравнение для траекторий типа а) и б) разрешается следующим образом:
p2 =-
4a2 L2
(vH -v02-a2L2)2 + 4a2L2v2H'
подставив это выражение во второе уравнение, получаем после преобразований:
t = — arcch a
f v02 + vH + a2L2 ^
2v0vH
что можно выразить также как:
H2 + L2 =
1+a H
V vo У
2
402 Sh2 a,.
v
br
252
При этом траектории типа а) получаются при условии vH — v0 — а2L2 > 0 , а траектории типа б) - при обратном знаке неравенства. Максимальная глубина траектории типа б) дается выражением:
V2(И )-(v + aH Y - 1 fe — v0 —а2L)2 + 4а2LVH
V (H max )-(V0 +aH max ) - ~2 --7^2-
p 4а L
и из условия того, что траектория не опускается ниже слоя разрыва, получаем дополнительное условие для существования траектории типа б):
aL <Vv2(Hbr) — v0 +4v2(Hbr) — vH .
Условие существования «головной» волны в)
aL >V vl — Vo2 + V v2br — vH — ^ v2br — v 2( Hbr) .
Можно показать, что выражение в правой части последнего условия всегда меньше, чем выражение в правой части условия существования траектории б), таким образом, в некотором диапазоне расстояний существуют одновременно и прямая, и головная траектории.
Таким образом, в зависимости от горизонтального расстояния между эпицентром и станцией реализуются следующие типы траекторий:
• aL <ijvH — v2 - восходящая траектория;
• Jv2H — vo2 < aL <^¡vl—v¡■ +J v2br— vH — 2 л/ v2br — v2 (Hhi) - восходяще-
нисходящая траектория;
• Ы— v0 W vl— vH — ^ vl— v2 (Hbr) <aL v 2(Hbr) — v2 + V v2(Hbr) — v2b
- восходяще-нисходящая траектория и головная волна;
• аЬ <-у]у2(ИЬг) - у02 + у!у2 (ИЬг) - у2И - только головная волна.
Сравнение времен прохождения сигнала для диапазонов, допускающих несколько типов траекторий, а также их интенсивности не представляется возможным провести аналитически и требует численного расчета.
4. ОГРАНИЧЕНИЕ РАССМАТРИВАЕМОГО АЛГОРИТМА
Необходимо отметить ограничения описанного алгоритма. Оно связано с применимостью линейной модели скорости сейсмических волн. Модельные
253
расчеты показывают, что для достаточно глубоких землетрясений, когда траектория сигнала проходит практически через весь слой выше зоны разрыва, различие между линейной моделью и более точной, моделирующей зависимость скорости от глубины с помощью сплайна, практически отсутствуют, то есть линейная модель вполне приемлема. Для мелких землетрясений наблюдаются отличия на небольших расстояниях между эпицентром и принимающей станцией, когда траектория проходит в приповерхностном слое, имеющем другой характерный градиент скорости.
Приведем пример таблицы, которую следует заполнять при проведении расчета. При этом вверху расположена таблица для данных Р-волн, ниже аналогичная таблица для S-волн. Для каждой станции следует ввести в строки (t, с) экспериментальные времена пробега для Р- и S-волн, а также веса данных времен, оценивающие точность снятия времени сигнала на станции. Предлагается использовать для указания веса данного измерения цифры по шкале от 0 до 5. При этом коэффициент 5 соответствует фазам +iPg, -iPg, iPg, iSg, iPn, коэффициент 4 соответствуют фазам ePg, eSg, iSn, коэффициент 2 соответствует фазам e(Pg), e(Sg). Впрочем, специалисту по обработке сейсмических данных легче выставлять эти коэффициенты при различном качестве имеющихся данных.
Заметим, что если в таблице имеются столбцы с весом равным нулю, это означает отсутствие данных у этой конкретной станции, поэтому времена, которые стоят в этих столбцах в строке t, с, не участвуют в расчете координат очага. При заполнении таблицы следует особое внимание обращать на близкие станции. Как было указано выше формулы расчета времени пробега для прямых и головных волн различны. В данном расчете очаг находится в Керченско-Анапской зоне очагов, поэтому в таблице справа стоят станции «Керчь» (KERU) и «Анапа» (ANN), относительно близкие к очагу. Для них в таблице теоретическое время рассчитывается по формулам прямой волны, отличным от формул головной волны. Однако иногда достаточно трудно решить, какая волна на данной станции приходит первой: прямая или головная. Для этого в данном примере мы ввели станции АNN и KERU как в группу прямых, так и в группу головных волн. Теперь достаточно посмотреть, как проходит теоретический расчет для этих станций. Там, где станция ANN стоит среди станций, на которые приходят головные волны, в строке Lgoi стоит отрицательное число, что и подсказывает, что в данном случае на станции АNN регистрируется прямая волна. Это означает, что следует поставить вес 5 (или иной в соответствии с качеством пришедшего сигнала) для станции ANN в правой части таблицы, а в левой части выставить станции ANN вес равный 0. А вот для станции KERU ситуация оказалась противоположной, числа в колонке KERU положительные. Расчет для нее должен проводиться по формулам головной волны.
5. ПРИМЕР ТАБЛИЦЫ РАСЧЕТА КООРДИНАТ ОЧАГА
Приведем пример для реальной таблицы расчета координат очага. Верхняя таблица относится к данным Р-волны, нижняя к данным S-волны. Заметим, что в расчетной таблице справа стоят две близкие к очагу станции, в них теоретический
254
расчет времени идет по формулам прямой волны. Для далеких станций теоретический расчет времени идет по формулам головной волны. Пользователь вводит в расчетные таблицы только экспериментальные времена пробега t, c для каждой станции, которая зафиксировала данное землетрясение. В рассматриваемом примере станции «Ялта» (YAL), «Судак» (SUDU) и «Севастополь» (SEV) имеют большие экспериментальные времена пробега, для этих станции расчет теоретического времени происходит по формулам головной волны. Для станций «Керчь» и «Судак» такой определенности нет. Поэтому эти станции поставлены как в группу далеких, так и в группу близких станций. Просматриваем строку Lgoi и видим, что для станции KERU величина Д км (эпицентральное расстояние эпицентр
- станция) равно 10.5 км, а для станции ANN ситуация противоположная, величина Lgol отрицательная. Поэтому для станции ANN расчет теоретического времени пробега будет проведен по формулам прямой волны. Аналогично заполняется и Таблица 2, в которой стоят данные для S-волны. При заполнении данных о временах на станциях в строках t, c для Р- и S-волн программа автоматически вычисляет итоговую среднеквадратическую невязку для Р- и S-волн с учетом введенных коэффициентов качества снятия данных. Понятно, что данные тех столбцов, в которых w = 0 в расчетах не участвуют.
Поясним некоторые переменные. Д - эпицентральное расстояние эпицентр -станция, hbr = 35 км- глубина слоя с линейным изменением скорости по глубине, Lgol
- длина пути головной волны по границе раздела сред, Vo = 5.4 для Р- волны (V0 = 3.2 для S-волны) - скорости на глубине h = 0 для Р- и S- волн соответственно. V(hbr) = 6.6 (3.8) - скорость волн на верхней границе раздела сред, Vbr = 7.5 (4.4) -скорость движения головной волны по границе раздела сред, h - глубина очага, V(h)
- скорость в очаге.
Таблица 1.
Данные Р-волны
Станция KERU YAL SUDU ANN SEV KERU ANN
t, с 16.0 33.9 25.0 9.2 38.4 16.0 9.2
Х, км 195.4 12.3 79.2 263.5 -25.5 195.4 263.5
Y, км 101.2 9.7 54.3 53.7 16.1 101.2 53.7
Вес 4 1 4 0 1 0 5
Argl 0.17 0.16
Arg2 1.13 1.05
Д 83.7 214.4 150.3 48.0 251.85 83.7 48.0
Lgoll 8.2 138.8 75.75 -27.5 176.3 14.9 8.9
tteor 16.2 33.6 25.1 11.4 38.6 14.9 8.9
(tteor t) 0.04 0.08 0.01 5.01 0.05 1.28 0.07
255
Таблица 2.
Данные 5-волны
Станция KERU YAL SUDU ANN SEV KERU ANN
t, с 28.1 57.4 32.1 14.9 65.3 28.1 14.9
Х, км 195.4 12.3 79.2 263.5 -25.5 195.4 263.5
Y, км 101.2 9.7 54.3 53.7 16.1 101.2 53.7
Вес 4 1 4 0 1 0 5
Argl 0.29 0.28
Arg2 1.10 1.03
Д 83.7 214.4 150.3 48.0 251.85 83.7 48.0
LBoii 10.5 141.2 77.1 -25.2 178.6 25.4 15.2
tteor 27.8 57.5 25.1 19.7 66.0 25.4 15.2
(tteor-t) 0.08 0.01 0.01 23.0 0.52 7.51 0.10
При заполнении данных о временах на станциях в строках t, c для Р- и 5-волн программа автоматически вычисляет итоговую среднеквадратическую невязку для Р- и S-волн с учетом введенных коэффициентов качества снятия данных. Суммарная невязка находится в ячейке «Невязка».
Невязки для Р- и 5-волн складываются и находятся в ячейке «Невязка». Минимизация проводится программой Excel «Поиск решений», которая находится в разделе «Данные» программы Excel. При этом курсор должен стоять в ячейке, в которой находится суммарная невязка. Для осуществления процесса минимизации необходима начальная точка процесса. Если имеются координаты данного очага, вычисленные заранее. Можно поставить эти данные в колонку «Каталог» в качестве начальной точки. Если таких данных нет, то следует ставить в качестве начальной точки произвольную точку в районе данного очага, можно взять в качестве начальной точки координату ближайшей станции.
Расчет Каталог
44.32 44.33
33.12 33.00
26.59 7
0.34 0.73
-70.17 -79
-8.37 -8
Невязка | 0.34 |
Минимум ищется по параметрам очага: ф. X. к. Опыт работы с программой показал, что глубина является самым вариабельным параметром, поэтому рекомендуется проводить минимизацию невязки трижды: сначала по X и У. потом по глубине к и окончательно еще раз по всем трем параметром. Итоговая невязка
256
есть показатель того, насколько вычисленные координаты очага соответствуют исходным временам Р- и 5-волн.
ЗАКЛЮЧЕНИЕ
Для примера приведем результаты пересчета координат некоторых очагов по данной программе. В таблице 3 приводятся данные 7 землетрясений 2011 г. находящиеся в Керченско-Анапской зоне очагов. В первых трех колонках приведены данные о землетрясениях из каталога 2011 г. [1]. Далее приведены значения суммарной среднеквадратической невязки экспериментальных и теоретических времен. Эти землетрясения были переопределены с помощью программы, о которой рассказано выше.
Таблица 3.
Пример результатов переопределения очагов землетрясений
Дата Ф0° V h0, с 8, с Ф° Г h, км 8, с
12.03.11 44.70 36.81 20 0.3 44.71 36.77 20.9 0.2
17.03.11 43.39 36.13 31 0.4 43.41 36.02 30.0 0.3
10.07.11 45.63 32.95 12 0.3 45.71 33.12 12.2 0.1
24.07.11 45.65 33.19 20 1.6 46.17 34.14 35.0 1.2
26.08.11 43.89 35.23 25 0.2 43.78 35.26 35.0 0.1
15.09.11 45.19 36.85 11 0.4 45.18 36.84 11.0 0.4
25.10.11 44.60 36.96 10 0.25 44.59 36.96 21.7 0.1
Сравнение координат очагов каталога [1] и переопределенных очагов по программе, описанной выше, показывает небольшие отличия координат переопределенных очагов. Следует выделить землетрясение 24.07.11, у которого исходная невязка была заметно больше чем у остальных очагов. Пересчет координат по приведенной программе уменьшил невязку, хотя она и продолжает оставаться достаточно большой.
Список литературы
1. Свидлова В. А., Сыкчина З. Н., Козиненко Н. М., Антонюк Г. П., Сусин Д. А., Антонюк В. А., Курьянова И. В., Росляков А. В. Каталог и подробные данные о землетрясениях Крымско-Черноморского региона за 2011 г. // Сейсмологический бюллетень Украины за 2011 год. Севастополь: НПЦ «ЭКОСИ-Гидрофизика», 2011. С. 94-135.
THE CALCULATION OF SEISMIC WAVES IN A LINEAR LAW OF INCREASING SPEED AND AVAILABILITY OF BORDER BREAK SPEED
Pivovarova N. B.
State autonomous institution of Crimean Republic «Crimean Expert Council on Seismic Hazard Assessment and Earthquake Prediction» SAICR «CEC», Simferopol, Russian Federation E-mail: pivnb@yandex.ru
257
Range of problems related to propagation of seismic signals in the rock strata of the Earth, the signals are caused by earthquakes, may be separated to several consecutive tasks: to determine of earthquake hypocenter coordinates for a given velocity profile and the known times of arrival of a signal at several stations; to precise the parameters of cross section in seismic velocities using the set of adjusted focus coordinates. In this paper we restrict ourselves to the first task. It is assumed that the velocity depends only on vertical coordinate zone a scale for which the effect of the curvature of the earth's surface is negligible, the stations are considered as situated on the plane day surface. Station can receive signals propagating along the trajectories of different types, because the question is important, what types of trajectories are possible for different positions of the hypocenter and at what trajectory signal comes first.
Comparison of earthquake focus coordinates in catalogue and coordinates redefined following described program reveals small differences.
Assuming that the velocity grows monotonically to a depth of Hbr, which has a spasmodic speed increase to the value vbr and below this speed zone remains almost constant, one can imagine several types of trajectories of signal depending on the relative position of the epicenter and the station:
a) a rising trajectory with an epicenter above break, signal all the time covered up;
b) grinded with gradient-rising trajectory-epicentre located above the gap, the signal is spread first down, then up, constantly staying above the zone of gap;
c) «Headache» wave signal is propagated down to the gap to which suitable with that value of the angle (e) that further dissemination occurs horizontally on the zone, then break the signal goes from the zone up to the same value (e);
d) a rising trajectory with the epicenter zone below break signal is propagated up the line to break zone, where an abrupt changes of direction, then covered up.
In all cases the decision is to determine the radiation parameter p from the first equation and substituting it into the second, in the case of head waves this parameter is the length of trajectory along the gap Lbr.
Comparison of times passing signal for ranges allow several types of trajectories, as well as their intensity is not possible to carry out analytical and requires numerical calculation. In General, this task allows for numerical solution, in particular cases may also analytic solution.
Keywords: the earthquake, the coordinates of the centers, the velocity field, Poisson's ratio, Young's modulus
References
1. Svidlova V. A., Sykchina Z. N., Kozinenko N. M., Antonyuk G. P., Susin D. A., Antonyuk V. A., Kur'yanova I. V., Roslyakov A. V. Katalog i podrobnye dannye o zemletryaseniyah Krymsko-Chernomorskogo regiona za 2011 g. (Catalog and details of earthquakes Crimean Black Sea region for 2011) Sejsmologicheskij byulleten' Ukrainy za 2011 god. Sevastopol: NPC «EHKOSI-Gidrofizika», 2011, pp. 94-135. '
Поступила в редакцию 25.11.2016 г.
258
