CC BY
70
УДК 550.34.01
Т. А. Тушко
АЛГОРИТМ РАСЧЕТА ГОДОГРАФА ДЛЯ СРЕДЫ, СОДЕРЖАЩЕЙ СЛОЙ С ЛИНЕЙНОЙ ЗАВИСИМОСТЬЮ СКОРОСТИ ОТ ГЛУБИНЫ
Рассматривается алгоритм расчета годографа рефрагированной и головной волн в заданной модели среды. Дано подробное обоснование используемых аналитических решений, вытекающих из общей постановки прямой кинематической задачи сейсмики в лучевом приближении. Предполагается использование алгоритма в задачах гипоцентрии.
Ключевые слова: годограф, лучевое приближение, лучевой параметр, рефрагированная волна, головная волна, гипоцентр.
Во многих случаях трехмерное скоростное строение земной коры удобно аппроксимировать конечным набором одномерных моделей, поставленных в соответствие различным участкам исследуемой среды. Модель подобного типа была построена для описания неоднородностей земной коры на территории Тывы и юга Красноярского края [1], характеризующихся высокой сейсмоактивностью, с целью решения задач, связанных с обработкой наблюдений за местными землетрясениями. Анализ сейсмической изученности рассматриваемой части Алтае-Саянской складчатой области позволил выделить участки (блоки) земной коры с примерно однородным строением среды, и каждому блоку поставить в соответствие несложную скоростную модель, предполагающую линейный закон возрастания скорости распространения упругих волн с глубиной в пределах всей земной коры или значительной ее части. Скорость в верхах мантии предполагается постоянной, а граница Мохо между корой и мантией - горизонтальной. Параметры такой модели являются результатом осреднения скоростных разрезов, полученных различными методами исследований.
Применение математической модели среды к задачам гипоцентрии предполагает наличие соответствующего алгоритма решения прямой кинематической задачи. Данная работа посвящена обоснованию алгоритма расчета годографа рефрагированных и головных сейсмических волн в условиях заданной модели.
Рассматриваемый алгоритм строится на основе аналитических решений, описывающих волновой процесс в лучевом приближении для плоской Земли. В работах [2; 3] дано уравнение годографа для среды с линейным возрастанием скорости с глубиной, при условии, что источник расположен на дневной поверхности. Расчетные формулы, применимые к сейсмическому лучу от источника, расположенного на глубине, не приводятся. Однако именно этот случай представляет интерес для нахождения координат гипоцентра землетрясения, особенно его глубины. В работе подробно рассмотрены вопросы расчета лучевого параметра и времени движения рефрагирован-ной волны от заглубленного источника до заданной точки полупространства. Особое внимание уделено вопросу идентификации лучевой траектории в зависимости от эпицентрального расстояния и скоростных
характеристик среды. Описанный алгоритм решения прямой задачи принципиально отличается от численных аналогов, основанных на варьировании угла выхода луча из источника либо на идее динамического программирования.
Расчет времени движения рефрагированной волны между двумя точками среды. Пусть задана декартова система координат. Плоскость (х, y) совпадает с дневной поверхностью, ось (0, z) направлена вглубь Земли. Задано положение источника M0( х0, y0, z0) и принимающей волну станции
M ( х, y, z), расположенной на поверхности Земли (z = 0) на эпицентральном расстоянии
Л = \/(х-х0)2 +(y-y0)2 от источника.
Скорость распространения упругих волн в среде зависит только от глубины z и описывается следующей функцией:
v( z) = v0 + a-z, (1)
где v0 и а - параметры модели.
Требуется найти время распространения упругой волны из точки М0 в точку M . В геометрическом (лучевом) приближении оно рассчитывается как время движения волны вдоль некоторой траектории -луча, соединяющего точки среды M0 и M. В каждой своей точке луч перпендикулярен к фронту волны и согласно принципу Ферма является траекторией, обеспечивающей наименьшее время пробега. Если скорость в среде не зависит от х и y, то луч лежит в вертикальной плоскости, проходящей через точки
M0 и M.
Время движения волны в среде и скоростная функция в ней связаны интегрально [2; 3]:
M J
Т = jf-, (2)
M0 v(z)
где ds - элемент дуги луча.
Из геометрических соображений ясно, что ds = dz /cose(z), где e(z) - угол между касательной к лучу и вертикалью в точке луча на глубине z. Кроме того, отношение sin e(z)/ v(z) = p называется луче-
вым параметром p и согласно закону Снеллиуса является константой вдоль всей лучевой траектории. Выразив отсюда cose(z) = -у/1 -p2 • v2(z) и подставив в выражение (2) получим
dz
Т = í
v(z>y/l - p2 • v2(z)
(З)
1 pv
Т=-í
rv V
d %
pv0
(4)
Полученный интеграл является табличным [4]. Воспользовавшись таблицей и вернувшись к исходным параметрам интегрирования, получим
Т = -- ln
а
1 + ^ 1 - p2 • (v0 + а • z)2 p • (Vo + а^ z)
(5)
где z0 и г - пределы интегрирования по глубине.
Прежде чем вычислить интеграл (5), рассмотрим с помощью рис. 1 возможные траектории луча с тем, чтобы правильно задать пределы интегрирования.
Т( Mol, M ) = -ln
(
1 + 71 - p2 • v2(z)
p •v(z)
zo =0
= Aln
z0 v(zгу1 - p
С целью вычислить интеграл воспользуемся (1) и перейдем к интегрированию по параметру dv = а • dz . Далее введем еще одну замену. Учитывая, что p • v = sin e, и обозначив p • v = E, можно записать
dv = dE¡p, что позволит перейти к интегрированию по угловому параметру. Выражение (3) в этом случае приобретает вид
і+>/т
- p2 • vo2 p • V(z)
p • vo
1 + yj 1 - p2 • V2(z)
Учитывая, что в точке наибольшего погружения 2 угол между лучом и нормалью составляет 90°, выразим скорость как у(2) = 1/р и, подставив, получим
Т( Mol, M ) = —ln
2 2 p • v0
p • vo
(б)
Легко показать, что полученное выражение совпадает с известным уравнением годографа рефрагированной волны, приведенным в [2]:
2
Т (A) = — Arsh
а
( A-а^ v2• vo У
где Д - эпицентральное расстояние между точками М01 и М.
Перейдем к рассмотрению другой траектории луча -(М03,М). Этот луч содержит лишь восходящую часть дуги, на которой лежат точки М03 и М. Время движения волны по этой траектории также рассчитываем исходя из выражения (5):
Т( M 0З, M ) = -1п
і+ут
2 2 - p • V
Л
p • V
^ (1 + 41 - p2 • v0 )(vo + а zo )
, -(^ + 41 -p2 •( + аzo)2)
(7)
= — ln
а
Исходя из соображений симметрии, время движения по дуге (М02, М) найдем как разность
Т(М01,М) - Т(М03,М):
Рис. 1. Лучевые траектории (М01, М02, М03) - возможные
положения источника; М - приемник; 2 - точка наибольшего погружения луча
Известно, что в одномерной среде с линейным законом возрастания скорости с глубиной лучи являются дугами окружности. Поэтому луч, соединяющий точки поверхности (на рис. 1 это М01 и М), будет симметричен относительно точки своего наибольшего погружения в среду 2 . Время движения по такому лучу, согласно формуле Герглотца-Вихерта [3], равно удвоенному времени движения от поверхности до точки 2. Применим это к (5):
Т(Mo2M) = а ln
(1 + 41 - p2 •v2 )х
1 x(l + 41 - p2 •( + а zo )2 У
p2 • vo •( + аzo)
(8)
На основе выражений (6)-(8) можно строить алгоритм решения прямой задачи. Однако для этого прежде всего требуется найти значение лучевого параметра р.
Значение лучевого параметра в источнике. Напомним, что в принятой модели среды луч, соединяющий точки М 0 (источник) и М (приемник), является дугой окружности радиуса К с центром в некоторой точке С, расположенной на высоте к = у0/а над дневной поверхностью.
o
z=0
Рис. 2. К оценке лучевого параметра
Значение лучевого параметра p = sin e(z)/v(z) определяется углом выхода луча из источника, обеспечивающего его попадание в точку M. Воспользуемся геометрическими построениями в плоскости луча.
Рассмотрим два возможных варианта положения источника (рис. 2): источник M01 расположен на поверхности и совпадает с началом координат (0,0) (слева); источник M02 имеет глубину z0 > 0 (справа). Ось (0, z) направлена вглубь, ось (0, D) проходит через эпицентр (проекция источника на дневную поверхность) и точку M, в которой находится приемник. Эпицентральное расстояние А считается известным. Точка zp есть проекция точки наибольшего погружения луча z на ось (0, D).
Пусть глубина источника z0 = 0. Угол e выхода луча из точки M01 связан соотношением
sin e = sin(90°-P) = cosp = h/R,
где h = v0/a, R = 4(А / 2)2 + (v0 / а)2. Вычислим чение лучевого параметра в точке M01:
зна-
p = • -1 = ■
1
«-•R vo а^(Д /2)2 + (v0 /а)2 2
(9)
i
а2 • Д2 + 4v2
Если глубина источника z0 > 0, то угол выхода луча из точки M02 (рис. 2, справа) можно задать выражением sin e = sin (90° - Р) = cos р = (h + z0) / R.
Чтобы найти радиус дуги R, воспользуемся тем, что эпицентральное расстояние Д нам известно. Представим его в виде суммы: Д = D1 + D2, где D1 -расстояние, пройденное лучом по горизонтали до точки максимального погружения, D2 - после нее. Радиус можно выразить двояко - как расстояние от точки C до точки M02 и от точки C до точки M:
R = 4h2 + D22 и R =4(h + z0)2 + D2.
Выразив Д как (Д - В2) и приравняв выражения для К найдем Б2, а затем и радиус дуги:
r = , h2 + +Д2 +2h ^z0)2
4Д2
Учитывая, что h = v0/а, найдем p:
h + zn
p=
cos p _ v(z0) R • (v0 +а • z0) 2
(10)
Í
4vq + (Д • а + z0 • Д 1 • (2v0 + а • z0 ))2
Случай, когда источник находится ближе к приемнику, чем точка наибольшего погружения дуги луча (траектория (М03, М) на рис.1), является симметричным по отношению к рассмотренному. Поэтому формула (10) будет справедлива и для этого случая.
Лучевая траектория и глубина погружения луча. Зная координаты источника и приемника, можно определить форму лучевой траектории, соединяющей их, что позволит правильно применить формулы (6)-(10) при построении алгоритма решения прямой задачи. Форма (или тип) траектории определяется положением точки максимального погружения относительно отрезка дуги, соединяющего источник с приемником. Если глубина источника 20 > 0, точка может как принадлежать отрезку, так и находиться вне этого отрезка. На рис. 1 это траектории (М02, М) и (М03,М). Используем критерий распознавания типа траектории, основанный на сравнении заданной глубины источника с некоторым расчетным ее значением, являющимся функцией известных нам данных (20, Д, у0, а). Проиллюстрируем этот подход (рис. 3).
Три источника расположены на оси (0,2) на разной глубине (20, 21 и 22). Все они находятся на эпи-центральном расстоянии Д от приемника, расположенного в точке М. Источникам соответствуют лучи (дуги) с индексами 0, 1, 2 соответственно.
Дуга 1 пересекает ось (0, г) в точке г = г1 под прямым углом, так как является частью окружности с центром С1, расположенным
радиуса R = 4Д + к2
на оси (0, г). Точка максимального погружения дуги здесь совпадает с положением источника г = г1. Поэтому луч, вышедший строго горизонтально из источника, попадет в точку М.
С2 c0
2 1 1 ° Г - -ї~ M , h
Z0 D
Z1
""''í
Z 1
z1 = R - к =
4v0 + (а*Д)2 - v0
а
(11)
Подставив значение лучевого параметра, найденное по формуле (9), если глубина источника г0 = 0, и по формуле (10), если г0 > 0, получим соответственно
д/(д-А)2 + 4уд~ - 2Ур
2а
Рис. 3. Лучевые траектории, проведенные из точки М
Если источник поместить выше точки г1, например в точку г0, то лучевая траектория, соединяющая источник и приемник, будет содержать точку наибольшего погружения дуги 0. Если же источник расположить на оси (0, г) ниже точки г1 (в точке г2 ), то лучевая траектория не будет содержать точку наибольшего погружения дуги 2. Наша задача: определить тип траектории для заданных источника и приемника.
Если рассматривать величину г1 как некоторое переключательное значение глубины источника, то можно построить критерий идентификации типа лучевой траектории. Для заданных параметров среды и эпицентрального расстояния вычислим значение г1:
д/4у2 + (Д • а + z0 • Д 1 • (2v0 + а • z0 ))2 - 2v0 z2 = ------------------------------------------.
2 2а
Сравнение рассчитанной величины Z с глубиной подошвы линейного слоя позволит внести необходимые ограничения или поправки в алгоритм решения прямой задачи.
Расчет головной волны от границы Мохо. Головными называют волны, часть траектории которых состоит из скольжения вдоль границы раздела двух сред. Это возникает, когда синус угла падения волны на границу равен отношению скоростей в средах. Лучевой параметр для траектории такого типа равен p = 1/ v2, где v2 - скорость распространения волн под границей Мохо.
Представим лучевую траекторию головной волны как сумму трех отрезков (рис. 4), время движения вдоль которых рассчитывается отдельно. Введем обозначения: Т1 - время погружения луча из источника М0 до границы Мохо; Т2 - время восхождения луча от границы до приемника, расположенного в точке М; Т3 - время скольжения луча вдоль границы.
Введем параметр zz = z1 - z0, где z0 - глубина рассматриваемого источника.
Тогда случаям zz < 0, zz > 0, zz = z1 будут соответствовать свои типы траекторий, соединяющих источник с приемником, и, следовательно, свои наборы аналитических решений для расчета лучевого параметра и времени пробега луча из числа полученных выше.
Рассмотрим одно важное замечание. При расчете времени движения волны, рефрагированной в слое с линейной зависимостью скорости от глубины, предполагалось, что луч не погружается ниже границы слоя. Это касается лучевых траекторий, для которых параметр zz > 0. Глубину наибольшего погружения луча z можно найти из условия, что луч в этой точке направлен под прямым углом к вертикали:
sin e( z) 1 1
Рис. 4. Головная волна
Граница предполагается горизонтальной и залегает на глубине гм.
7] находим на основе решения (5):
1 (1 + 41 -Р2 • (у0 +“• 2о)2 )• (у0 +“• 2м)
7] =—• 1п-----------------/—. (12)
(v0 +а*Zo)-(і + 41- Р2 •(vo + а*zm )2
T2 - в соответствии с (7):
Í
T — ln
а
1 1 I1 + 41 -Р2 *v0 )*(v0 +а*ZM)
I •I1 +41 - Р2 *(v0 +а*ZM )2 )
. (13)
Р = ■
v( Z) v( Z) v0 +а • Z
Чтобы найти Т3, необходимо знать, какое расстояние Ь прошел луч вдоль границы. Из рис. 4 видно, что Ь = Д- Д - Д2, где Д - эпицентральное расстояние. Пусть йх - расстояние, пройденное падаю-
щим на границу лучом по горизонтали и соответствующее элементу дуги луча ds. Из геометрических соображений ясно, что dx = ds • sin e, где e - угол луча с вертикалью. В то же время dx = dz¡cos e, где dz - отрезок пути по глубине. Выразим путь луча по горизонтали:
ZM i Zm
( M qM]
,) =í
dz
cose(z)
• sin e
= í
dz • p • v(z)
Vі - P2 • v2(z)
. (14)
Введя подстановки, аналогичные тем, что при вычислении интеграла (3), получим
X,
( Mq,Mi)
i zM
ы-
v%_
Vm?
)
^•а V /
M
где £ = р-V (2).
Подставляя в качестве пределов интегрирования глубину залегания границы 2 = 2м и глубину источника 2 = 20, получим путь по горизонтали для нисходящей части луча:
D = V2•(>/1 -P2 • (Vq +а• zq)2 -- Vі-P2 • (v0 + аzM)2).
(15)
Интегрируя (14) от 2 = 2м до 2 = 0, найдем путь для восходящей части луча:
V2-(V1 -р2 -VI2 ^>/1—P2^(VГ+a-ZM)I). (16)
Общее время пробега луча Т найдем исходя из решений (12)-(16) как сумму
Т = Т + Т2 + Т3 = Т + Т2 + (Д - -Д - —2 ) / V2 ,
при условии, что (Д- Д - Д) >0.
Таким образом, исходя из общей постановки прямой кинематической задачи сейсмики в лучевом приближении, дано обоснование необходимых аналитических решений для разработки и программной реализации алгоритма расчета годографа в принятой модели среды. Особое внимание уделено случаю, когда источник сейсмических волн расположен на глубине.
Предложен критерий, позволяющий однозначно идентифицировать форму лучевой траектории, соединяющей источник с приемником, при расчете годографа рефрагированной волны от заглубленного источника.
Рассмотрен подход к расчету годографа головной волны от подошвы линейного слоя.
Представленный в работе алгоритм может быть также использован для нахождения численного решения обратных задач, требующих многократного решения прямой задачи. Алгоритм основан на аналитических решениях, обладающих устойчивостью и не дающих большой погрешности в вычислениях.
Библиографические ссылки
1. Тушко Т. А., Осеев В. Г., Пилимонкин Н. С. Построение скоростной модели среды для решения задач гипоцентрии // Соврем. методы обработки и интерпретации сейсмологических данных : материалы V Междунар. сейсмологич. школы. Обнинск, 2010. С. 209-214.
2. Саваренский Е. Ф., Кирнос Д. П. Элементы сейсмологии и сейсмометрии. М. : Гос. изд-во тех.-теор. литературы, 1955.
3. Аки К., Ричардс П. Количественная сейсмология: теория и методы : в 2 т. М. : Мир, 1983.
4. Двайт Г. Б. Таблицы интегралов и другие математические формулы. М. : Наука, 1978.
z
о
T. A. Tushko
ALGORITHM OF CALCULATION OF THE HODOGRAPH CONTAINING A LAYER WITH A LINEAR DEPENDENCE OF VELOCITY ON DEPTH
An algorithm of calculation of hodograph of refracted and head waves in an assumed environment model is considered. A detailed argumentation of the applied analytical solutions, derived from the universal arrangement of the direct amplitude time problem in ray approximation, is presented. The algorithm is offered to be used in hypocenter location problems.
Keywords: hodograph, ray approximation, ray parameter, refracted wave, head wave, hypocenter.
© Тушко Т. А., 2011
