CC BY
11
УДК 004.942
А. Г. Марчук
Институт вычислительной математики и математической геофизики СО РАН пр. Акад. Лаврентьева, 6, Новосибирск, 630090, Россия
Новосибирский государственный университет ул. Пирогова, 2, Новосибирск, 630090, Россия
E-mail: mag@omzg.sscc.ru
МИНИМИЗАЦИЯ ПОГРЕШНОСТЕЙ ПРИ ЧИСЛЕННЫХ РАСЧЕТАХ ВОЛНОВЫХ ЛУЧЕЙ И ФРОНТОВ ЦУНАМИ *
При расчетах кинематики волновых лучей и фронтов в неоднородных средах с переменной скоростью распространения возмущений наиболее эффективными являются методы, базирующиеся на принципе Гюйгенса. В данном подходе траектории лучей восстанавливаются после численного расчета времен пробега волны из источника во все остальные точки расчетной области. Однако тот факт, что построенные таким методом волновые лучи проходят только через узлы расчетной сетки, приводит к погрешностям, которые не зависят от шага расчетной сетки. Эти погрешности можно минимизировать путем описываемой в работе процедуры уточнения траектории луча и времен распространения цунами между соседними узлами сетки. Тестовые расчеты иллюстрируют заметное улучшение точности вычислений времен пробега цунами.
Ключевые слова: кинематика волны цунами, время прихода волны, узлы расчетной сетки, волновой луч, оптимизация траектории луча.
Исследование волновых лучей и фронтов цунами численными методами позволяет решать некоторые задачи оперативного прогноза цунами. В частности, задачу оценки времен прихода волны в заданные точки акватории и побережья. Для численного моделирования цунами используется прямоугольная расчетная сетка, в узлах которой задается глубина океана. Идея поиска оптимального пути в области, где скорость проводимости задана в узлах прямоугольной сетки, принадлежит, вероятно, Danzig [1], который предлагал этот метод для решения задач, совершенно не связанных с геофизикой. Позднее такой подход применялся для нахождения фронтов сейсмических волн. В нашей стране этот подход для решения задач, связанных с волнами цунами, вероятно, впервые применил Карев [2]. Существует несколько способов нахождения волновых лучей, являющихся, по сути, экстремалями (т. е. трассами, по которым возмущение приходит из одной точки области в другую за кратчайшее время) [3]. Некоторые методы неэффективны при решении краевой задачи для волнового луча или при расчетах в областях со сложной топологией (острова, проливы и т. п.). Другие не очень приспособлены к работе с сеточной цифровой батиметрией.
Наиболее универсальными в этом плане являются методы, основанные на принципе Гюйгенса. Они приспособлены к работе именно на прямоугольных сетках, а также без проблем позволяют строить волновой луч между двумя заданными точками расчетной области. Этот
* Работа поддержана контрактом Министерства образования и науки Российской Федерации № 14.В37.21.0643, грантом РФФИ № 12-07-00406 и интеграционным проектом СО РАН № 117.
Марчук А. Г. Минимизация погрешностей при численных расчетах волновых лучей и фронтов цунами // Вестн. Новосиб. гос. ун-та. Серия: Информационные технологии. 2013. Т. 11, вып. 3. С. 27-36.
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2013. Том 11, выпуск 3 © А. Г. Марчук, 2013
принцип, по сути, состоит в том, что все точки акватории, до которых волновое возмущение пришло к данному моменту времени, становятся в свою очередь источниками волнового возмущения и излучают волновую энергию во все стороны. Таким образом, алгоритм расчета строится на переборе всех прилегающих к волновому фронту расчетных узлов сетки (в которые возмущение еще не пришло) и нахождении времен прихода волны в эти узлы путем минимизации суммы времен пробега волны от очага до каждой из соседних точек (куда волна к этому моменту уже пришла) и от этой соседней точки до рассматриваемого узла сетки. Поясним это рисунком. Пусть расчет кинематики цунами проводится в области с прямоугольной расчетной сеткой. Это означает, что в узлах этой сетки известны значения глубин и требуется найти времена добегания волн цунами от заданного источника (одного или нескольких узлов сетки) до всех остальных узлов, до которых волна в принципе может дойти. Изобразим схематически на рис. 1 участок расчетной области. Здесь черными квадратиками обозначены те узлы расчетной сетки, в которые возмущение от источника к данному моменту времени уже пришло, и нам известны времена прихода волны в эти узлы. Нам требуется найти время добегания волны от источника до точки А (см. рис. 1). По отношению к точке А соседними узлами, где время прихода волны уже известно, будут точки В, С, В, Е, ¥, О и Н.
в ^ А
1 1 1 1 ^л ■ / 1 1 \ 1
. е. / 0 1 1 V | | н
Рис. 1. Схема работы алгоритма, основанного на принципе Гюйгенса
Пусть времена добегания до них равны соответственно ТВ, Тс, Тв, Те, Т¥, То, и Тн. Скорость распространения волны цунами зависит только от глубины и определяется формулой Лагранжа
V=4ёН, (1)
где я - ускорение силы тяжести, Н - глубина. Считается, что между соседними узлами сетки изменение глубины происходит по линейному закону. Найдем время пробега волны цунами между двумя точками, расстояние между которыми равно Ь, а глубина меняется от величины Нх до величины Н2. Введем вспомогательную величину а = (Н2 -Нх) / Ь (угол наклона дна между этими точками). Тогда время пробега запишется в виде
ь
Т = /
(И
1
/- /—11
^ • (н +1 • & а) ^ • & а 01
I +
I+н
tgа
1/2 (1
(
Н 2 -Л/ Н1
_Н_
\0
Н 2 - Н1
л/я^а фяа _= 2 Ь
Следовательно, время движения волны цунами между соседними узлами расчетной сетки равно расстоянию между ними, деленному на среднее арифметическое скоростей цунами в
этих узлах. Таким образом, для того, чтобы найти время добегания от источника до точки А (см. рис. 1), требуется найти минимум из семи времен Т, (/ = 1, ..., 7):
Т = Тв +
2-Ах
Т = Т +
12
2 -V (Ах)2 + (Ау)2
4ёнС+4§На
Т = Т +
2-Ау
Т = Т +
Т = Т + 16
2 -V (Ах)2 + (Ау)2
2 (Ах)2 + (2Ау)2
Т = Т +
15 1И ^
Т7 = Тр +
2 -у] (Ах)2 + (2Ау)2
4ёйН+у1ша '
2 -у] (2Ах)2 + (Ау)2
(2)
4ШР+4§На
где Ах , Ау - шаги сетки в горизонтальном и вертикальном направлениях, а НА, Нв, НС, Н0, НЕ, Нр, Н0, НН - значения глубины в соответствующих точках. Минимальное из значений Т, (} = 1, ..., 7) даст нам время пробега волны цунами от очага до точки А. Таким способом можно поочередно найти времена прихода волны цунами во все узлы расчетной сетки.
Скажем несколько слов о понятии соседних точек. Соседние точки по отношению к рассматриваемой определяются так называемым шаблоном. Наиболее широко используется шаблон (например, [2]), состоящий из шестнадцати точек, окружающих точку А расчетной сетки, в которую требуется определить время прихода волны. В таком шаблоне соседними к точке А с сеточными координатами (/, /) являются шестнадцать узлов расчетной сетки, индексы которых отличаются от сеточных координат точки А на единицу и некоторые из узлов, одна из координат которых различается на двойку (рис. 2). Таким образом, времена пробега из каждого из них в точку А вычисляются по шестнадцати лучам (отрезкам). Поэтому этот шаблон иногда называют шестнадцатилучевым (шестнадцатиточечным).
] •
И*
]-2-
Рис. 2. Шестнадцатиточечный шаблон для вычисления времени пробега цунами на прямоугольной расчетной сетке
В случае если в какой-либо из 16 точек шаблона время добегания цунами уже известно, то можно найти время прихода волны и в точку А (центр шаблона). Затем, многократно перебирая все точки области, постепенно вычисляем времена прихода цунами из источника (один или несколько узлов расчетной сетки, где в начальный момент задается время, равное нулю) во все остальные точки расчетной области.
В принципе для нахождения времени прихода в какой-либо узел сетки можно использовать полученное непосредственно перед рассмотрением этого очередного узла значение в одной из соседних точек. Однако это может привести к ошибкам (иногда значительным). Если применять описанный алгоритм «прямолинейно», т. е. в случае, если хотя бы в одной из шестнадцати точек шаблона становится известным время добегания, то, используя формулу (2), вычисляется время добегания в рассматриваемый узел расчетной сетки (в центральную точку шаблона). Тогда результаты расчета во всей области будут сильно зависеть от направления перебора уз-
лов расчетной сетки. В частности, индекс 7 можно менять как от минимального значения до максимального, так и наоборот. То же относится к индексу / (положение точек расчетной сетки по вертикали).
Рассмотрим конкретный пример. В прямоугольной расчетной области ЛБСО дно имеет вид наклонной плоскости, при этом глубина линейно возрастает от стороны ЛО к стороне БС. Пусть в центре области в точке М(70, /0) расположен точечный источник цунами (рис. 3). Пусть направлениями перебора расчетных узлов являются направления слева направо и постепенное движение вниз. Это означает, что индекс 7 меняется от единицы на левой границе (ЛБ) до значения 7тах (на стороне СО), а индекс / меняется от единицы на верхней границе (БС) до /тах на нижней (ЛО). В начале процесса расчета имеется единственная точка расчетной сетки с известным временем добегания в нее - точка с координатами (70, /0), в которой значение времени добегания равно нулю. После начала перебора всех расчетных узлов области во всех точках, где индекс/ меньше /0 - 2 (т. е. в верхней части области), никаких значений времен добегания вычислено не будет, так как в окрестностях всех этих узлов (см. шаблон на рис. 2) не оказалось ни одной, где время добегания было бы известно. При переборе точек (/0 - 2)-й строки в окрестности точки (70 - 1, /0 - 2) окажется источник цунами - М с координатами (70, /0), т. е. появится возможность вычислить время добегания в эту точку. Согласно формуле (2) имеем
2 (V (Ах)2 + (2Ау)2
Т
0-1, л-2
I I Т1—1 —I— I #гм
1, Л -2
>/Н 70 +4 ёН 0 -1,.
Далее становится возможным определить времена прихода цунами во все точки этой строки, которые расположены правее точки (7 0 - 1, /0 - 2), так как очевидно, что каждая из этик точек будет иметь в своей окрестности как минимум одну точку с известными временами прихода туда волны. Продолжая эти рассуждения, можно установить, что после первого прохода (перебора) области в указанной ранее последовательности получится картина, изображенная на рис. 3.
Рис. 3. Состояние расчетной области после первого прохода (перебора) расчетных узлов. Серым цветом закрашена та часть области, где найдены времена добегания
Серым цветом закрашена та часть расчетной области, где времена добегания определены. Следовательно, в оставшейся части области времена еще не вычислены. Рассмотрим точку Л, расположенную у правой границы области и лежащую в той же строке, что и источник цунами М. Согласно алгоритму расчетное время добегания туда волны примерно равно времени движения возмущения вдоль прямой МЛ:
I
Т ■■
где I - расстояние между точками М и Л, НЛ - глубина в точке N. Однако, как показано в [4], самой быстрой траекторией движения волны из точки М в точку N является циклоида, которая схематически изображена на рис. 3. Следовательно, практически во всех расчетных узлах (/ + 2)-й строки, расположенных правее источника, времена добегания будут найдены неверно!
А так как алгоритм не предполагает пересмотра найденных значений времен добегания, то по окончании расчета во всей области (после многократных проходов) мы будем иметь ошибочную карту изохрон цунами.
При использовании других направлений просмотра точек ошибочные значения времен добегания могут появиться в других частях расчетной области. Для того чтобы исключить подобное, необходимо многократно перебирать расчетные узлы и во время каждого такого прохода вычислять времена добегания лишь в точках, непосредственно прилегающих к точкам, в которые к этому моменту волна уже пришла (т. е. в которых определено время прихода). Это можно реализовать следующими действиями. Перед началом расчетов выбирается некоторый условный шаг по времени At, величина которого приближенно определяется выражением
. [Ax/JgHZ
At ~ min < _,
[aj/^H:
где Hmax - максимальное значение глубины в рассматриваемой области, а Ax и Ay - расстояния между узлами расчетной сетки в горизонтальном и вертикальном направлениях. Во время первого прохода (просмотра точек расчетной области) значения времен добегания в новых точках присваиваются только в том случае, когда вычисленное там по описанному алгоритму значение времени добегания будет меньше, либо равно At. На втором проходе эта временная граница становится равной 2At, а на просмотре расчетной области с номером n времена, присваиваемые в новых точках, будут ограничиваться величиной n ■ At. Таким образом, для безошибочного определения времен добегания во все расчетные точки необходимы многократные переборы точек расчетной области, при этом во время некоторых проходов ни одной новой точки с вычисленным в ней временем прихода к уже имеющимся может не добавиться.
Такой алгоритм сопряжен с многократным перебором точек расчетной области, что влечет за собой достаточно большие вычислительные ресурсы (при расчетах в больших областях). Например, расчет в сеточной области размером 3 000 на 2 000 узлов требует до 10 минут на персональном компьютере с процессором Intel Core 2 Duo 2500 МГц. А при расчете трансокеанских цунами на достаточно мелкой сетке количество расчетных узлов может вырасти в несколько раз. Поэтому возникает стремление как-то оптимизировать процесс расчета. Оптимизацию можно проводить по двум направлениям: сокращать количество перебираемых узлов во время каждого прохода области и уменьшать число арифметических действий при вычислении времени пробега цунами между точками шаблона. Объем вычислений, приходящийся на каждую точку, можно легко сократить следующим образом. В расчетных формулах (2) нигде не требуются значения собственно глубины. Везде используется только скорость распространения волны в данной точке, вычисляемая из значения глубины по формуле Лагранжа (1). Поэтому можно заранее перед началом расчета вычислить значения скоростей волны во всех точках расчетной области и ввести в оперативную память компьютера массив значений этих скоростей. Что касается сокращения количества перебираемых точек во время каждого прохода, то здесь можно все время держать в поле зрения текущие размеры области, куда распространилась волна, и вести перебор узлов только в прямоугольнике, описанном вокруг этой области. Например, если в какой-то момент индексы всех узлов, куда дошла волна, содержатся в интервалах
i . < i < i ;
min max'
i ■ < i < i
J min J J max'
тогда при следующем просмотре области не обязательно перебирать все точки расчетной области, а можно ограничиться перебором узлов, где индексы i, j ограничены интервалами
(i ■ - 2) < i < (i + 2);
min max
(jmin - 2) < j < (jmax + 2).
Очевидно, что такой прием может дать значительное убыстрение расчетов только на начальной стадии процесса, когда область возмущения занимает небольшую часть всей расчетной области. После того, как область возмущения приблизится ко всем границам расчетной
области, данный прием потеряет смысл. Однако при расчетах кинематики цунами в обширных областях такой метод оптимизации может значительно ускорить дело.
Следует заметить, что при расчете кинематики волнового фронта с помощью метода, основанного на принципе Гюйгенса, сами волновые лучи остаются не найденными. Для того чтобы построить волновой луч, соединяющий две заданные точки расчетной области, предлагается следующий алгоритм. В одну из заданных точек помещаем точечный источник волн и по методу, использующему принцип Гюйгенса, находим времена добегания волны из этого источника до всех возможных точек расчетной области. В процессе расчета мы формируем массив, в который записывается информация о том, из каких конкретно соседних узловых точек волна пришла в каждую из расчетных точек. После окончания расчета становится возможным восстановить траекторию наискорейшего движения волны из точки источника до любого из других расчетных узлов области. При этом восстановленный луч будет состоять из таких же отрезков, как и расчетный шаблон (см. рис. 2). Это значит, что при использовании шестнадцатиточечного шаблона имеется лишь 16 вариаций направления сегментов волнового луча. Проиллюстрируем процедуру восстановления траектории волнового луча на примере.
Ю 11 12 13 14 15 16
(11,11) (13,12)
(10,10)
(14,12)
ИСТО' шик
Рис. 4. Схема восстановления волнового луча между точками (10, 10) и (16, 11)
На рис. 4 приведен небольшой фрагмент расчетной области. Пусть источник находится в точке с сеточными координатами (10, 10). Требуется найти траекторию волнового луча, соединяющего источник с точкой (16, 11). Пусть при расчете времен добегания цунами от источника до остальных точек расчетной области выяснилось, что в точку (16, 11) наибыстрейшим образом волна пришла из точки (14, 12). Следовательно, в точку (14, 12) волна прибыла из точки (13, 12). И так далее вплоть до источника (см. рис. 4). Эти данные в процессе расчета времен прихода волны сохраняются в двух вспомогательных двумерных числовых массивах. Если проследить эту цепочку от точки приемника (16, 11) и до источника (10, 10), то таким образом восстанавливается траектория волнового луча (серая ломаная линия на рис. 4). В качестве примера на рис. 5 представлен волновой луч, который восстановлен описанным методом при расчете времен добегания волны цунами от точечного источника у входа в Авачинскую бухту до трех других точек у побережья Камчатки. Однако построенные таким методом волновые лучи не являются гладкими линиями, а представляют собой ломаные, составленные из отрезков, пространственное направление которых имеет только 16 возможных вариаций (соответствующих лучам на рис. 2). Здесь следует заметить, что при расчетах в реальных акваториях сетки привязаны к географическим координатам и поэтому не являются в строгом смысле прямоугольными. Ячейки имеют вид трапеций, где при удалении от экватора длина шага вдоль параллелей уменьшается пропорционально косинусу широты. Вблизи экватора длина географического градуса вдоль меридианов и параллелей примерно одинакова.
Рис. 5. Волновые лучи у побережья Камчатки
Эти траектории можно считать первым приближением действительных волновых лучей. Для большей точности предлагается проводить процедуру уточнения траектории. Суть ее состоит в том, что по окончании расчета каждую из точек, через которые проходит волновой луч (первое приближение), перемещаем вдоль линий сетки с целью нахождения нового положения точки, при котором время пробега вдоль измененной ломаной будет минимальным. И так проделываем последовательно с каждой из точек первого приближения волнового луча. В результате траектория луча становится более гладкой, так как она уже не обязательно проходит через узлы исходной расчетной сетки.
V
40 -35 -
30 -25 -
20 -15 -
10 - ;
5 -:
т
о
Рис. 6. Сравнение формы волнового луча до (слева) и после оптимизации
В качестве примера на рис. 6 приводится волновой луч над параболическим дном до процедуры оптимизации (слева) и после нее (справа). Здесь глубина возрастает пропорционально квадрату расстояния от левой границы расчетной области. Форма «оптимизированного» луча хорошо соответствует точному решению, которое имеет вид отрезка окружности [4]. Похожий алгоритм в свое время был предложен Брэддоком [5; 6]. Там также при помощи вариаций траектории волнового луча искался путь, дававший минимальное время пробега вдоль него. Однако за начальное приближение обычно бралась прямая линия (геодезическая) между источником и приемником. В этом случае результатом оптимизации мог стать волновой луч, дающий локальный минимум времени пробега, в то время как оптимальный путь может проходить совсем в другом месте.
Аналогичную процедуру уточнения можно проводить в процессе расчета времен добега-ния. Ее схема изображена на рис. 7, являющегося фрагментом рис. 2.
Рис. 7. Схема процедуры оптимизации при расчете времени прихода волны в точку А
Пусть требуется найти время прихода волны в точку А, используя при этом значение времени в точке С. В ходе расчета выяснено, что в точку С волна быстрее всего дошла из расчетной точки В. Далее, варьируя положение точки С1, находим такое ее положение, которое дает минимум времени пробега волны вдоль ломаной ВС1А. Именно это значение приплюсовывается к значению времени в точке В для нахождения времени добегания волны в точку А. И так для каждой точки области. Таким образом, уже в ходе численных расчетов времен прихода волны в узлы расчетной сетки происходит корректировка (уточнение) этих времен.
Рис. 8. Результат расчета изохрон цунами от точечного Рис. 9. Визуализация изохрон цунами над ровным источника над ровным дном без оптимизации дном с применением оптимизации
Наиболее жестким тестом для методов расчета кинематики волновых фронтов является расчет изохрон волны цунами от круглого (точечного) источника в области с постоянной глубиной. Использование обычного алгоритма с шестнадцатиточечным шаблоном приводит к тому, что вместо окружностей получаются многоугольники (рис. 8). Описанная процедура уточнения в значительной степени позволяет исправить картину (рис. 9). При этом разница между оптимизированным и неоптимизированным массивом времен прихода может превышать в некоторых точках одну минуту на каждый час распространения цунами.
Может сложиться впечатление, что волновой луч (оптимальную траекторию) между двумя заданными точками расчетной области можно легко найти чисто вариационным методом, считая первым приближением луча прямую линию, соединяющую эти точки. К сожалению, не все так просто. Во-первых, потребуется большое количество «проходов» вдоль траектории с корректировкой местоположения точек волнового луча. Во-вторых, что более существенно, таким
способом удастся отыскать локальную экстремаль, которая не обязательно является наибыстрейшим путем распространения цунами. Пример такого рода уже приводился. В заключение можно заметить, что предложенный метод может хорошо работать и в трехмерном случае, например, при решении задач сейсмики.
Рис. 10. Изохроны волны цунами от очага 29.09.2009 и ее записи глубоководными регистраторами системы DARTS
На рис. 10 представлены изохроны волны цунами (положение волнового фронта в определенные моменты времени) через каждые 300 с при ее распространении над донным рельефом из базы данных ОЕБСО. Жирные изолинии обозначают время, кратное 1 500 с. Видно, что волна от точечного источника, расположенного в точке эпицентра землетрясения, доходит до трех глубоководных регистраторов за 3 879 с до БЛЯТ51425, за 4 026 с до БЛЯТ51426 и за 7 613 с до станции БЛЯТ54401. Реально зафиксированные времена вступления цунами этими регистраторами на 150-200 секунд меньше полученных в ходе численного расчета. Это объясняется тем, что в численном эксперименте время прихода вычислялось от точечного источника, а учитывая то, что при таких глубинах, как в области эпицентра землетрясения, скорость волн составляет 250-300 м/с, при размерах очага порядка 70-100 км эта разница во временах добегания (вычисленных и реальных) цунами практически исчезнет.
Список литературы
1. Dantzig G. On the Shortest Route Through a Network // Management Science. 1960. Vol. 6. No 2. P. 187-190.
2. Карев В. Ю., Судаков А. Н., Чубаров Л. Б. Методы расчета карт изохрон // Методы интерпретации экспериментальных данных и численные расчеты. Красноярск: КГУ, 1987. С. 75-81.
3. Марчук А. Г. Численные методы расчета кинематики волн цунами // Математические проблемы геофизики: Численные исследования геофизических задач: Сб. науч. тр. / Под ред. А. С. Алексеева; Вычислительный центр СО АН СССР. Новосибирск, 1988. С. 69-90.
4. Marchuk A. G., Moskalensky E. D. Some Analytical Solutions for Tsunami Wave Rays and Front // Bulletin of the Novosibirsk Computing Center. Series: Mathematical Modeling in Geophysics, 2010. No. 13. P. 113-126.
5. Braddock R. D., Doilibi P., Voss G. Wave Charting for Tsunamis // Proceedings Internetional Tsunami Symposium / Ed. by T. S. Murty and W. J. Rapatz. Sidney, British Columbia, Canada, 1985. P.264-274.
6. Braddok R. D. Tsunami Propagation over Large Distances. Tsunamis in the Pacific Ocean / Ed. by W. M. Adams. Honolulu, Hawaii: East West Center Press, 1970. P. 285-303.
Материал поступил в редколлегию 26.05.2013
A. G. Marchuk
MINIMIZING COMPUTATIONAL ERRORS OF TSUNAMI WAVE-RAY AND TRAVEL TIME
Accuracy of computational methods for tsunami kinematics is very important in obtaining reliable results. Methods based on Huygens principle are the most effective for computations on the rectangular computational grids. The present study proposes modification in the methodology of determining tsunami travel-times and of wave-ray paths. An approximate ray trace path can be developed from a source origin point to any other point on a computational grid by solving directly the problem - and thus obtain the tsunami travel times. The initial ray approximation can be optimized with the use of an algorithm that calculates all potential variations and applies corrections to travel-time values.
Keywords: tsunami wave kinematics, wave travel time, computational grid, wave-ray, ray path optimization.
References
1. Dantzig G. On the Shortest Route Through a Network // Management Science. 1960. Vol. 6. No. 2. P. 187-190.
2. Karev V. U., Sudakov A. N., Chubarov L. B. Metody rascheta kart izohron // Metody interpretacii eksperimentalnyh dannyh i chislennye raschety. Krasnoyarsk: KGU, 1987. S. 75-81.
3. Marchuk A. G. Chislennye metody rascheta kinematiki voln cunami // Matematicheskie problemy geofiziki: Chislennye issledovaniya geofizicheskih zadach: Sb. nauch. tr. / Pod red. A. S. Alekseeva; Vychislitelnyy centr SO AN SSSR. Novosibirsk, 1988. S. 69-90.
4. Marchuk A. G., Moskalensky E. D. Some Analytical Solutions for Tsunami Wave Rays and Front // Bulletin of the Novosibirsk Computing Center. Series: Mathematical Modeling in Geophysics, 2010. No. 13. P. 113-126.
5. Braddock R. D., Doilibi P., Voss G. Wave Charting for Tsunamis // Proceedings Internetional Tsunami Symposium / Ed. by T. S. Murty and W. J. Rapatz. Sidney, British Columbia, Canada, 1985. P. 264-274.
6. Braddok R. D. Tsunami Propagation over Large Distances. Tsunamis in the Pacific Ocean / Ed. by W. M. Adams. Honolulu, Hawaii: East West Center Press, 1970. P. 285-303.
