«ТРЕХМЕРНЫЕ СКОРОСТНЫЕ МОДЕЛИ ЗЕМНОЙ КОРЫ ТЯНЬ-ШАНЯ НА ОСНОВЕ БИ-СПЛАЙН ПАРАМЕТРИЗАЦИИ И ТРИАНГУЛЯЦИИ ДЕЛОНЕ ...»

ИНСТИТУТ ДИНАМИКИ ГЕОСФЕР

РОССИЙСКОЙ АКАДЕМИИ НАУК

На правах рукописи

Усольцева Ольга Алексеевна

ТРЕХМЕРНЫЕ СКОРОСТНЫЕ МОДЕЛИ ЗЕМНОЙ КОРЫ

ТЯНЬ-ШАНЯ НА ОСНОВЕ БИ-СПЛАЙН

ПАРАМЕТРИЗАЦИИ И ТРИАНГУЛЯЦИИ ДЕЛОНЕ

Специальность 25.00.10 - геофизика, геофизические методы поисков полезных ископаемых

ДИССЕРТАЦИЯ

на соискание ученой степени кандидата физико-математических наук

Научный руководитель – доктор физико-математических наук Ирина Альфатовна Санина Москва - Оглавление.

СОКРАЩЕНИЯ, ТЕРМИНЫ И ОБОЗНАЧЕНИЯ

ВВЕДЕНИЕ

ГЛАВА 1 ИСТОРИЯ РАЗВИТИЯ СЕЙСМИЧЕСКОЙ ТОМОГРАФИИ И

СУЩЕСТВУЮЩИЕ ТОМОГРАФИЧЕСКИЕ МОДЕЛИ

1.1 Предыстория создания первых трехмерных моделей скоростей сейсмических волн

1.2 Некоторые существующие трехмерные модели скоростей Р и S волн.................. Выводы к главе

ГЛАВА 2 МЕТОДИКА ПОСТРОЕНИЯ ТРЕХМЕРНЫХ РЕГИОНАЛЬНЫХ

ТОМОГРАФИЧЕСКИХ МОДЕЛЕЙ ПО ДАННЫМ ВРЕМЕН ПРОБЕГА ОТ

ЛОКАЛЬНЫХ ЗЕМЛЕТРЯСЕНИЙ

2.1 Томография по временам пробега объемных сейсмических волн от землетрясений. Постановка задачи. Введение понятия параметризации функции.......... 2.2 Формулы для вычисления частных производных. Различные способы параметризации модели.

2.2.1 Блоковая параметризация.

2.2.2 Параметризация с использованием Би-сплайнов

2.2.3 Параметризация модели с помощью тетраэдров, внутри которых скорость меняется линейно непрерывно

2.2.4 Несколько слов о других способах параметризации модели

2.3 Используемые стандартные статистические характеристики для оценки качества полученных результатов

2.4 Алгоритмы построения сейсмических лучей

2.4.1 Построение лучей в сферически-симметричной среде (одномерный луч)... 2.4.2 Построение лучей в горизонтально-неоднородных средах

2.5 Методы обращения матриц. Вычисление матрицы разрешения. Решение обратной задачи

2.5.1 Нахождение решения уравнения (34) с помощью построения системы нормальных уравнений

2.5.2 Нахождение решения уравнения (34) с помощью сингулярного разложения матрицы B

2.5.3 Нахождение решения уравнения (34) методом LSQR

2.6 Методы локации землетрясений в одномерных скоростных моделях

2.7 Различные способы введения весовых коэффициентов

2.8 Основные этапы любого сейсмотомографического исследования

2.8.1 Сортировка данных

2.8.2 Построение оптимальной одномерной модели и перелокация событий в полученной одномерной модели. Вычисление станционных поправок

2.8.3 Выбор алгоритма и его тестирование

2.8.4 Построение трехмерной модели и перелокация событий, проверка правильности проведенных расчетов.

Выводы к главе

ГЛАВА 3 СРАВНИТЕЛЬНЫЙ И ОЦЕНОЧНЫЙ АНАЛИЗ

3.1 Сравнение различных численных алгоритмов для расчета времени пробега луча и вычисления траекторий сейсмических лучей

3.2 Сравнение различных способов параметризации модели

3.3 Сравнение различных способов обращения матриц

3.4 Оценка влияния процедуры пересчета лучей и перелокации событий после каждой итерации на значение RMSw.

3.5 Тестирование различных алгоритмов построения сейсмотомографических моделей

Выводы к главе.

ГЛАВА 4 ОБЩАЯ ГЕОФИЗИЧЕСКАЯ ХАРАКТЕРИСТИКА ТЯНЬ-ШАНЯ....... 4.1 Расположение исследуемого региона на топографической и тектонической картах

4.2 История развития региона.

4.3 Магнитные аномалии.

4.4 Аномалии силы тяжести, изостатическое состояние земной коры, плотность пород.

4.5 Тепловое поле

4.6 Активные разломы, скорости современных деформаций.

4.6.1 Скорости горизонтальных движений.

4.6.2 Скорости вертикальных движений.

4.7 Исследования Тянь-Шаня по сейсмологическим данным.

4.7.1 Общая информация.

4.7.2 Сейсмичность

4.7.3 Исследования методами глубинного сейсмического зондирования (ГСЗ). 4.7.4 Исследования земной коры всего Тянь-Шаня методом сейсмической томографии по данным объемных волн

4.7.5 Исследования земной коры Cеверного Тянь-Шаня методом сейсмической томографии по данным объемных волн

4.8 Детальное изучение зоны сочленения Чуйской впадины и Киргизского хребта.....

Выводы к главе

ГЛАВА 5 ОПИСАНИЕ ИСПОЛЬЗУЕМЫХ В РАБОТЕ ДАННЫХ И АНАЛИЗ

СЕЙСМИЧНОСТИ ТЯНЬ-ШАНЯ.

5.1 Общие сведения об используемых данных

5.2 Анализ сейсмичности и месторасположения станций на территории Тянь-Шаня по имеющимся данным

5.3 Отбор данных для построения одномерных и трехмерных томографических моделей

5.3.1 Отбор данных для территории Cеверного Тянь-Шаня (41.9-43.4 с.ш., 73.5в.д.)

5.3.2 Отбор данных для всей территории Тянь-Шаня

Выводы к главе

ГЛАВА 6 ПОЛУЧЕННЫЕ ОДНОМЕРНЫЕ И ТРЕХМЕРНЫЕ СКОРОСТНЫЕ

МОДЕЛИ ТЯНЬ-ШАНЯ

6.1 Скоростные модели Северного Тянь-Шаня

6.1.1 Построение одномерных моделей, расчет станционных поправок и перелокация событий в одномерной модели для территории Северного Тянь-Шаня.....

6.1.2 Трехмерные модели Р и S волн для территории Северного Тянь-Шаня, перелокация событий в трехмерной модели

6.1.3 Геологотектоническая интерпретация полученных результатов и сравнение с существующими моделями для территории Северного Тянь-Шаня

6.2 Исследование всей территории Тянь-Шаня с использованием локальных данных.

6.2.1 Построение одномерных моделей, расчет станционных поправок и перелокация событий в одномерной модели для территории всего Тянь-Шаня........ 6.2.2 Трехмерные модели Р и S волн для территории всего Тянь-Шаня, перелокация событий в трехмерной модели

6.2.3 Геологотектоническая интерпретация полученных результатов и сравнение с существующими моделями территории всего Тянь-Шаня.

Выводы к главе.

ВЫВОДЫ И ЗАКЛЮЧЕНИЕ

ПРИЛОЖЕНИЕ 1 СУЩЕСТВУЮЩЕЕ ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ ДЛЯ

РЕШЕНИЯ СЕЙСМОТОМОГРАФИЧЕСКИХ ЗАДАЧ

1. Программа Velest [21]

2. Программа Sphypit90 [43]

3 Программа Simulps [55]

4 Программа Fatomo [19].

ЛИТЕРАТУРА

Сокращения, термины и обозначения ТФР – Таласо-Ферганский разлом ГСЗ – глубинное сейсмическое зондирование МОСМ – минимальная одномерная скоростная модель DSW – derivative weighted sum (сумма взвешенных значений частных производных по всем лучам для данного параметра модели) RDE – resolution diagonal elements (диагональные элементы матрицы разрешения) Vp – стандартная ошибка для скоростей Р волн (измеряется в км/с) dVp/Vp – стандартная ошибка для возмущений скоростей Р волн (измеряется в %) Vp/Vs – стандартная ошибка для трехмерной модели распределения отношения Vp/Vs Р волна – продольная сейсмическая волна S волна – поперечная сейсмическая волна T=tреальн-tрассчетн – невязка времен пробега (разница между реальным и расчетным временем пробега луча);

- эпицентральное расстояние между двумя точками по дуге большого круга;

RMS – среднеквадратичная невязка по всем лучам RMSw – взвешенная среднеквадратичная невязка по всем лучам RMSсоб - среднеквадратичная невязка для данного события RMSwсоб – взвешенная среднеквадратичная невязка для данного события Введение Актуальность работы.

Несколько последних десятилетий сейсмическая томография является одним из наиболее распространенных и эффективных методов получения информации о скоростных свойствах пород внутри Земли. Существует большое количество различных сейсмотомографических алгоритмов.

Наличие этих алгоритмов с одной стороны существенно упрощает и убыстряет решение задачи восстановления трехмерной скоростной структуры в Земле. С другой стороны, не зная как устроены эти алгоритмы и какова их область применимости, трудно получить удовлетворительное решение. В этой связи особенно актуальным на сегодняшний день является сравнительное изучение физических основ различных алгоритмов, освоение этих алгоритмов и их усовершенствование для конкретного набора экспериментальных данных и особенностей геолого-тектонического строения изучаемого региона.

На данном этапе развития сейсмической томографии наиболее актуальным является исследование сложно-построенных регионов, например, такого горного массива, как Тянь-Шань. К настоящему моменту для территории Тянь-Шаня очевидно, что кора и мантия существенно неоднородны. Характер этих неоднородностей более сложный, с различной степенью контрастности, чем, например, в зонах субдукции и в местах расположения плюмов (нет ярко выраженной высокоскоростной области, связанной с погружающейся пластиной океанской литосферы, как в районах субдукции, или низкоскоростного канала, связанного с восходящей струей разогретой мантии, как над плюмами). Поэтому нужны более совершенные методы для восстановления этих неоднородностей. До сих пор при построении сейсмотомографических моделей территории Тянь-Шаня изучаемая область разбивалась на прямоугольные блоки с постоянной скоростью внутри или скоростная функция раскладывалась по полиномам Лежандра (гармонически-полиномиальный способ разложения). Оба эти способа имеют ряд существенных недостатков. При использовании гармонически-полиномиального способа часто сталкиваются с проблемой ложной экстраполяции искомой функции в слабо изученных районах. При разбивке исследуемой территории на блоки с постоянной скоростью не всегда удается правильно установить границы разноскоростных областей.

Сравнительно недавно в Тянь-Шаньском регионе наряду с аналоговой формой записи сейсмического сигнала, стало возможным производить запись сейсмического сигнала в цифровом виде. Первые цифровые станции на территории Тянь-Шаня появились в 1991 году. Количество цифровых станций на территории Тянь-Шаня с каждым годом увеличивается. В 1997гг. их было около 40. Точность определения времени вступления различных волн по цифровой волновой форме значительно выше, чем с использованием аналоговой сейсмограммы. Очевидно, что из двух скоростных моделей для одного и того же региона, та которая получена с использованием данных на цифровых станциях является более точной и детальной. Построение томографических моделей с использованием цифровых данных описано в работах [11] и [67]. В данной работе представлен анализ цифровых данных за более длительный период времени и при построении моделей используются эти данные в совокупности с наиболее совершенными сейсмотомографическими методиками.

Разработка и усовершенствование имеющихся сейсмотомографических алгоритмов для сложно-построенных сред, а также расчет трехмерных скоростных моделей с помощью этих алгоритмов и новых массивов данных для территории Тянь-Шаня актуальны при проведении сейсмологических исследований горного сооружения Тянь-Шаня. Новая информация, полученная на основе анализа рассчитанных вновь скоростных моделей совместно с новыми уточненными тектоническими картами [86; 109], в значительной степени дополняет существующие представления об особенностях геологического строения и динамических процессах, происходящих в этом регионе.

Основная цель.

Целью работы является изучение пространственных скоростных неоднородностей строения земной коры Тянь-Шаня на основе непрерывного способа параметризации модели по данным цифровых и аналоговых сейсмических станций.

Основные задачи

исследования.

- провести сравнительный анализ наиболее часто используемых различных сейсмотомографических методов, а также алгоритмов, разработанных автором на тестовых примерах;

экспериментального материала, относительное расположение источников и приемников, размеры территории и др.) разработать усовершенствованные алгоритмы томографической инверсии;

- для получения достоверной информации о скоростях сейсмических волн в зоне сочленения Киргизского хребта и Чуйской впадины построить трехмерные непрерывные скоростные модели верхнего этажа земной коры по Р и S волнам с использованием различных алгоритмов и различных наборов данных для этой территории;

- построить наиболее адекватную скоростную модель земной коры по Р волнам для территории всего Тянь-Шаня.

Научная новизна.

Предложены два новых алгоритма TomoTetraFD и TomoCubeFD для построения непрерывных и квазинепрерывных (функция скорости претерпевает разрыв на конкретных глубинах) моделей скоростей TomoTetraFD особенно эффективен в районах с существенно неравномерным расположением источников сейсмических волн. В отличие от многих существующих алгоритмов применение данных возможно не только для территорий локального масштаба (200*200 км), но и регионального (1000*1000 км).

сейсмотомографических программ: Sphypit90, Simulps14, Fatomo, TomoCubeFD, TomoTetraFD, в которых в качестве исходных данных сейсмических волн от локальных землетрясений. Три первых алгоритма (Sphypit90, Simulps14, Fatomo) активно применяются для построения трехмерных моделей скоростей объемных волн различных регионов. Два других (TomoCubeFD и TomoTetraFD) разработаны автором на базе существующего сейсмотомографического программного обеспечения. Даны рекомендации по поводу того, в каких случаях удобнее всего использовать тот или иной алгоритм.

зарегистрированных на аналоговых и цифровых станциях, впервые построена квазинепрерывная трехмерная скоростная модель Р волн всей территории Тянь-Шаня. Предложен и реализован на практике ряд тестов, подтверждающих устойчивость полученного результата. Проведена геологогеофизическая интерпретация полученного результата. Впервые проведена локация большого количества землетрясений и взрывов, произошедших на всей территории Тянь-Шаня методом сеточного поиска [34] в трехмерной скоростной модели.

Применен новый подход к построению трехмерных скоростных моделей P и S волн для верхней части коры под Северным Тянь-Шанем, а также для анализа отношения Vp/Vs в этом регионе. Этот подход позволяет выявить наиболее достоверные скоростные неоднородности, т.к. включает в себя построение целого множества трехмерных скоростных моделей с сейсмотомографических алгоритмов для Р волн и совместно для Р и S волн.

Применение такого подхода необходимо при наличии системы наблюдений, состоящей из малого количества станций. Очень важно также использование такого подхода, когда в распоряжении исследователя имеются только бюллетени сейсмических событий, а не сами волновые формы.

Защищаемые научные положения.

Использование Би-сплайн параметризации и триангуляции Делоне в томографических алгоритмах является эффективным при построении трехмерных скоростных моделей. Достаточно эффективным способом выявления действительно существующих трехмерных скоростных неоднородностей является использование жестких критериев по отбору данных и построение не одной, а большого количества скоростных моделей с использованием различных наборов данных и различных методов.

Обнаруженные скоростные неоднородности в верхней коре под ТяньШанем хорошо согласуются с геологической и тектонической картой данного региона.

Практическая значимость.

Полученные в настоящей работе трехмерные скоростные модели могут быть использованы для уточнения геолого-тектонического строения ТяньШаня, повышения точности локации эпицентров землетрясений и выявления характерных особенностей современного сейсмического режима исследуемой территории с последующей геофизической интерпретацией.

Установленные скоростные модели могут быть также использованы в интересах Международной системы сейсмического мониторинга (IMS) для контроля за соблюдением Договора о всеобъемлющем запрещении ядерных испытаний. Полученные скоростные разрезы могут быть использованы для расчета очагово-станционных сейсмических поправок SSSC для станций IMS Тянь-Шаньского региона. В работе [42] показано, как на основе высокоточной трехмерной скоростной модели для Индо-Пакистанского региона произведен расчет поправок SSSC и с их помощью проведена более уточненная локация событий по станциям IMS.

Рекомендации методической главы данной работы существенно упростят и ускорят процедуру выбора подходящего алгоритма для проведения локальных сейсмотомографических исследований коры и верхней мантии в других сейсмоактивных регионах земного шара. В свою очередь эти сейсмотомографические исследования позволят получить новую информацию о геолого-тектоническом строении и провести уточнение координат и времени в источнике землетрясений в этих регионах.

Глава 1 История развития сейсмической томографии и существующие томографические модели предшествующие построению первых трехмерных моделей скоростей Р и S сейсмических волн. Проводится краткий обзор работ, в которых построены сейсмотомографические модели, сильно повлиявшие на существующие сейчас представления о строении Земли и процессах, происходящих в ее недрах.

1.1 Предыстория создания первых трехмерных моделей скоростей сейсмических волн К концу шестидесятых годов прошлого века развитие наук о Земле вышло на ту стадию, когда стало возможным и необходимым построение трехмерных горизонтально неоднородных моделей скоростей Р и S сейсмических волн.

Построение трехмерных моделей было необходимым потому, что с использованием только одномерных скоростных моделей невозможно было объяснить некоторые процессы, которые, как предполагалось, происходят в Земле. К таким процессам относятся взаимодействие литосферных плит, конвекция в мантии, воздействие на окружающую среду восходящих горячих мантийных струй (плюмы).

Возможность построения таких моделей объясняется рядом причин.

Во-первых, был накоплен обширный экспериментальный материал, который позволил выявить основные черты внутреннего строения Земли в рамках одномерных скоростных моделей. Была создана довольно точная сферически симметричная модель внутреннего строения Земли и на ее основе рассчитаны таблицы времен прихода сейсмических волн в зависимости от эпицентрального расстояния – хорошо известные таблицы ДжеффрисаБуллена и Херрина, EASPEI-91. Во-вторых, была создана плотная мировая сеть сейсмологических наблюдений, большие сейсмические антенны (LASA в США в 1965г., NORSAR в Норвегии в 1972г.), развиты национальные системы сейсмологических наблюдений. В-третьих, была разработана методика построения одномерных скоростных моделей, а также методы построения сейсмических лучей. В-четвертых, развитие вычислительной техники позволило реализовать предложенные алгоритмы построения сейсмических лучей и расчета трехмерных моделей скоростного строения Земли по объемным волнам.

1.2 Некоторые существующие трехмерные модели скоростей Первые сейсмотомографические модели появились около 30 лет назад.

На сегодняшний день существует большое количество глобальных и региональных скоростных моделей Земли. Для построения скоростных моделей отдельных регионов используют данные времен пробега объемных (Р и S) и высокочастотных поверхностных волн. В глобальных моделях в качестве начальных данных часто используются: короткопериодные продольные объемные волны, прошедшие без отражений, отраженные от земного ядра, а также прошедшие через внешнее земное ядро;

длиннопериодные объемные волны с периодами от 45 до 200 сек;

поверхностные волны с периодами от 35 до 250 сек.; собственные колебания всей планеты (3-55 мин). Региональные скоростные модели строятся по данным времен пробега сейсмических волн как от локальных землетрясений или взрывов (т.е. источников, расположенных на исследуемой территории), так и по данным времен пробега от телесейсмических землетрясений или взрывов (удаленных на большие расстояния от региона).

Одна из первых работ, посвященных построению региональной скоростной модели по локальным данным, была проведена K.Aki и W.H.K.Lee в 1976 году [3]. Они построили трехмерную модель скоростей Р волн для верхних 15 км коры в районе разломов Сан Андреас и Калаверас в Калифорнии и произвели перелокацию источников в этом районе.

Впоследствии глубинные разломы в Калифорнии методом сейсмической томографии изучались в работах [7; 8; 57; 64]. После перелокации событий в трехмерной скоростной модели большинство гипоцентров землетрясений располагаются ближе к тектоническим разломным структурам. Границы высоко- или низкоскоростных аномалий часто совпадают с линиями разломов, которые были обнаружены с помощью геологических методов. В перечисленных работах для расчетов использовались данные от землетрясений, произошедших в верхней коре до глубин 10-15 км. В работе [46] использовались данные от землетрясений, произошедших на значительно больших глубинах до 250 км. С помощью этих данных С.Рокеру [46] удалось выявить низкоскоростную область, простирающуюся под Памиро-Гиндукушской зоной в Средней Азии до глубин 250 км.

По данным удаленных землетрясений в нашей стране первая трехмерная скоростная модель была построена А.С. Алексеевым и М.М.

Лаврентьевым на основе наблюдений времен пробега на сейсмологическом профиле Памир-Байкал [69]. Одной из основополагающих работ телесейсмической томографии, в которой подробно описывается методика построения скоростных моделей верхней мантии с использованием данных от удаленных землетрясений, также является работа [2].

По количеству работ в области глобальной сейсмической томографии особое место занимают работы американской школы Адама Дзивонского [5].

На сегодняшний день для большинства глобальных сейсмотомографических моделей в верхней мантии на глубинах 78-148 км с сейсмически активными районами связаны низкоскоростные аномалии, а с сейсмически спокойными регионами высокоскоростные аномалии.

Исследованию пространственного строения литосферы по данным поверхностных волн посвящены работы Т.Б. Яновской [122].

Помимо трехмерных моделей скоростей Р и S волн с помощью сейсмической томографии по локальным данным также строятся карты отношения Vp/Vs для различных глубин [64], которые несут в себе информацию о химическом составе, слагающих земные недра пород.

тектонических структур Европейским Геофизическим Союзом проводился ряд исследований с привлечением различных стран. На территории Германии, Дании, Швеции проводился глобальный эксперимент TOR (Transeuropaen Suture Zone) в 1996/97 годах, на территории Финляндии и России эксперимент Svekalapko в 1999 году. Основная цель этих экспериментов – более детальное изучение внутреннего строения в этих регионах за счет создания плотной сети переносных и стационарных сейсмических станций (при проведении эксперимента Svekapko станции располагались через каждые 50 км). Полученные данные были использованы для изучения горизонтальных неоднородностей не только методом сейсмической томографии, а также методом обменных волн и другими.

Наиболее яркими сейсмотомографическими результатами являются трехмерные скоростные модели для зон субдукции, разломных зон и областей плюмов. Ниже на Рис. 1а представлено вертикальное сечение трехмерной скоростной модели, полученной А.Горбатовым для зоны субдукции под Камчатским полуостровом [13]. Скоростные сейсмотомографические модели для зон субдукции под Новой Зеландией построены в работе [64], для Идзу-Бонинской зоны (Тихий океан) в [63], для района Памир-Гиндукуш в [46]. С помощью сейсмической томографии во всех работах (аналогично Рис. 1а) удается проследить, как субдуцирующая холодная высокоскоростная литосферная плита пересекает астеносферный слой и следует наклонно до поверхности нижней мантии (670 км). На Рис. 1а построена скоростная томографическая модель в районе субдукции до глубин 200 км. Трехмерные скоростные неоднородности над плюмами обсуждаются в работах [54] (Гавайский плюм) и [10] (Исландский плюм).

Рис. 1 Вертикальные сечения трехмерных скоростных моделей: а) под Камчатским полуостровом из [13], б) сейсмотомографическая модель Исландского плюма [10].

В данных работах удается проследить, как огромная капля перегретого низкоскоростного вещества поднимается из недр к коре (см. Рис. 1б).

В последнее время с помощью сейсмотомографического метода также (характеристика затухания сейсмических волн) [51; 27], анизотропных скоростных моделей [18] и трехмерных моделей температур вещества, слагающего земные недра.

Выводы к главе Описанные в этой главе сейсмотомографические модели показывают происходящих в Земле. Благодаря этим моделям для многих районов земного шара выявлены основные скоростные закономерности, как в глобальном, так и в региональном масштабе, связанные с асейсмичными и сейсмоактивными зонами. В последние годы сейсмическая томография развивается в направлении увеличения детальности исследований и повышения точности метода. Стало возможным построение трехмерных томографических моделей не только в районах с простым геологическим строением (платформы, щиты, континентальные впадины), но и восстановление горизонтальных скоростных неоднородностей в сложно-построенных районах с неочевидной геологической историей, например, таких как ТяньШань. До сих пор детально не изучена трехмерная скоростная структуры коры и верхней мантии под Тянь-Шанем и эта тема является предметом широкого обсуждения и научных дискуссий.

Глава 2 Методика построения трехмерных региональных томографических моделей по данным времен пробега от локальных землетрясений томографических моделей. Вторая глава посвящена теоретическим основам сейсмотомографического алгоритма. В данной главе с позиции автора подробно описана суть предложенного K.Aki и W.H.K.Lee в 1976 году подхода и дальнейшие модернизации этого подхода, проведенные большим количеством ученых. При описании каждой компоненты метода отмечается, проводилось ли автором только изучение этой компоненты, или же дополнительно тестирование и совершенствование.

2.1 Томография по временам пробега объемных сейсмических волн от землетрясений. Постановка задачи. Введение понятия параметризации функции.

Как уже упоминалось во введении, сейсмическая томография позволяет, используя данные времен пробега объемных сейсмических волн, рассчитывать трехмерные статистически обоснованные модели распределения скоростей сейсмических волн в коре и мантии.

Время пробега сейсмической волны определяется функционалом Ферма, в котором интегрирование выполняется вдоль экстремали (луча).

v(r) – скорость сейсмической волны в точке r, L – траектория сейсмического луча, также зависящая от скорости сейсмических волн в среде, координат источника. Координаты землетрясений и время, когда оно произошло, точно не известны, поэтому эти характеристики вносят ошибку в определение траектории луча. Координаты приемников обычно известны с высокой точностью. В более общем виде время пробега сейсмической волны T, распространяющейся от некоторого землетрясения с координатами xист yист zист и временем в источнике tист на станцию, можно представить как некую функцию от значений скоростей v(r) в различных точках пространства, координат землетрясений xист yист zист и времени в источнике tист.

Основная задача заключается в определении скорости v(r), координат землетрясений и времени в источнике по множеству измерений времен пробега T на поверхности для различных лучей.

Определить функцию v(r) в явном виде сложно, поэтому при решении сейсмотомографических задач вычисляют не саму функцию v(r), а только некоторые ее осредненные характеристики. Замена искомой модели v(r) упрощенной моделью, полученной из v(r) путем некоторого осреднения, называется параметризацией модели. Упрощенная скоростная модель, возникшая в процессе параметризации, представима в виде суммы конечного числа параметров, каждый из которых умножен на некоторую базисную функцию:

hk(r) – М базисных функций, - различные М коэффициентов при этих функциях. Базисные функции представляют собой набор стандартных функций, значения которых обычно изменяются от нуля до единицы. В сейсмической томографии часто осуществляется параметризация не функции скорости, а функции медленности (величины, обратной скорости) [44], или квадрата медленности [11]. Все дальнейшие выкладки этого пункта верны в не зависимости от того, какая из функций (скорость, медленность или квадрат медленности) раскладывается по набору базисных функций.

Учитывая (3) равенство (2) можно переписать следующим образом:

Следовательно, время пробега представляет собой функционал, зависящий от конечного числа параметров, описывающих искомую скоростную модель, от координат и времени в источнике.

В основу сейсмотомографического метода заложены принципы теории возмущений.

Предположим, что известна начальная приближенная скоростная модель v0(r), приближенные координаты события xист0 yист0 zист0 и время в источнике tист0. При параметризации начальной скоростной модели возможно базисным функциям. Время пробега сейсмической волны от землетрясений с координатами xист0 yист0 zист0 и временем в источнике tист0 в начальной скоростной модели обозначим T0.

Если искомая скоростная модель, координаты источников и время в источнике незначительно отличаются от нулевого приближения, то можно воспользоваться разложением Тейлора по малому параметру для функции многих переменных. Здесь необходимо отметить, что используются слагаемые только первого порядка малости и пренебрегается слагаемыми более высокого порядка малости. Получаем:

Разница между временем пробега T и T0 равна:

Основная сейсмотомографическая задача свелась к определению поправок к коэффициентам, соответствующим нулевой модели, поправок к координатам источника и поправки ко времени в очаге по разнице времен пробега T для различных лучей. Другими словами, необходимо решить систему линейных уравнений, количество уравнений в которой равно количеству сейсмических лучей, а неизвестными являются величины k, x, y, z и t.

Как было сказано выше, задачу нахождения поправок k, x, y, z и t, зная невязки T, возможно свести к линейной системе уравнений только в том случае, если при разложении функции T в (6) пренебречь членами с k2, x2, y2, z2, t2, k3, x3 и.т.д. Т.к. слагаемые, в которые искомые поправки входят во второй, третьей и более высоких степенях не учитывается, решение системы линейных уравнений типа (7) не является решением, поставленной задачи. Чтобы получить истинные значения k, x, y, z и t процедуру решения линейной системы уравнений типа (7) необходимо повторить итеративно несколько раз согласно формуле (8) до тех пор, пока значения искомых поправок на очередной итерации не будут пренебрежимо малы.

где i и i+1 – номера итераций. В этом случае искомые поправки k, x, y, z и t равны сумме поправок, полученных после каждой итерации (9).

Далее рассматриваются различные способы решения системы уравнений типа (7) только на первой итерации, т.к. методы, применяемые при решении систем на каждой итерации, одинаковы.

В матричном виде система уравнений типа (7) имеет вид N – количество лучей, М – количество коэффициентов при базисных функциях, К – количество землетрясений. Матрица А, элементами которой являются частные производные времени пробега по коэффициентам (производные Фреше), в дальнейшем будет называться матрицей из частных производных.

лучей i=1,…,N есть результат решения прямой задачи. Обратная задача уравнению (10), зная матрицы A, Hn и вектор T.

2.2 Формулы для вычисления частных производных.

Различные способы параметризации модели.

Обозначим функцию скорости - v(r), функцию медленности - s(r) и функцию квадрата медленности - s2(r). Выражение для времени пробега через функцию скорости приведено в (1), ниже приведены выражения для времени пробега через функцию медленности (11) и через функцию квадрата медленности (11):

При параметризации модели искомая функция раскладывается по некоторому набору базисных функций с определенными коэффициентами. В (12) приведен общий вид разложения функций s(r) и s2(r) по набору базисных функций.

Эти разложения аналогичны разложению (3) для функции v(r). Исходя из (3) и (12) с учетом (1) и (11) общий вид формул для вычисления производных Фреше для трех видов функций следующий:

Из формул (13) видно, что проще всего вычислить производные, если параметризовать функцию медленности. Производные, вычисляемые для функций s(r) и s2(r) всегда положительны, следовательно, матрица A в (10), которую они составляют, является положительно определенной. Если положительно определенной является матрица –A.

Параметризацию модели возможно проводить в классе кусочнопостоянных, в классе непрерывных, в классе один раз непрерывнодифференцируемых, дважды-непрерывно дифференцируемых (и.т.д. до бесконечности) функций. Ниже более подробно описаны те виды параметризаций, которые использовались автором в ходе работы.

2.2.1 Блоковая параметризация.

Блоковая параметризация применяется и в декартовой и в сферической системе координат и является наиболее распространенной, как в глобальной, так и в локальной томографии. Исследуемый объем разбивается на конечное число блоков. В декартовой системе координат блоки представляют собой параллелепипеды, в сферической системе координат сферические сегменты.

Считается, что искомая функция (v(r), s(r) или s2(r)) постоянна внутри каждого блока. Количество базисных функций равно количеству блоков.

Каждая базисная функция имеет вид:

Коэффициент, соответствующий k-ой базисной функции равен значению искомой функции (vk(r), sk(r) или s2k(r)) в k-ом блоке. При блоковой параметризации формулы (13) для вычисления производных Фреше принимают вид где lik – длина i-того луча в k-ом блоке. Знак интеграла исчезает, т.к.

подынтегральная функция постоянна. Блоковая параметризация используется для построения глобальных томографических моделей BDP98 (Boschi and Dziewonski), HWE97 (van der Hilst et al.) и других. Также блоковая параметризация применяется для построения региональных скоростных моделей в [3; 45; 44].

Одним из недостатков блоковой параметризации является то, что из-за наличия фиктивных вертикальных границ часто возникают проблемы при трассировки лучей и не всегда верно восстанавливаются границы разноскоростных блоков.

2.2.2 Параметризация с использованием Би-сплайнов Параметризация с использованием Би-сплайнов применяется только в декартовой системе координат, т.к. эта параметризация требует задания прямоугольной сетки. В исследуемом объеме задается сетка из точек. Каждая точка сетки - это точка пересечений трех прямых, первая из которых параллельна оси x, вторая параллельна оси y, третья параллельна оси z. Эти прямые отсекают на осях x,y и z соответственно некоторые отрезки. Длина отрезков, отсекаемых на оси x, может быть неодинакова, т.е. шаг сетки по оси x может быть неравномерным. По оси y и z также допустима неравномерность шага сетки. Для построения базисных функций в трехмерном пространстве сначала необходимо построить 3 одномерных функции, соответствующих каждой из осей. Трехмерная Би-сплайн базисная функция равна произведению трех одномерных базисных функций.

Теория построения одномерных Би-сплайнов изложена в [85]. П.

Фирбас в [112] считает, что Би-сплайновое описание среды является оптимальным при построении томографических моделей. Би-сплайн аппроксимация является аппроксимацией общего типа и поэтому не накладывает серьезных ограничений на форму возмущающей функции.

Одномерные и соответственно трехмерные Би-сплайны бывают 1, 2, 3 и.т.д.

порядка. Одномерный Би-сплайн 1-го порядка представляет собой ступенчатую функцию, т.е. такую функцию, которая в двух точках претерпевает разрыв 1-го рода (Рис. 2). Соответственно трехмерные Бисплайны первого порядка совпадают с базисными функциями, введенными при описании блоковой параметризации (14). Формулы для вычисления одномерных Би-сплайнов 2-го порядка (Би-сплайны 2-го порядка также называются линейными Би-сплайнами) и Би-сплайнов 4-го порядка (Бисплайны 4-го порядка также называются кубическими Би-сплайнами) для случая равномерной сетки приведены в [33]. На Рис. 2 представлен вид одномерных линейной и кубической базисной функции.

Рис. 2 Три различных типа одномерной базисной функции: а) Би-сплайн первого порядка (используется при блоковой параметризации), б) Би-сплайн второго порядка (используется при параметризации линейными Би-сплайнами), в) Би-сплайн четвертого порядка (при параметризации кубическими Би-сплайнами).

Ниже приводим формулы для вычисления одномерной базисной Бисплайн функции 4-го порядка Bix для случая неравномерной сетки. Точки сетки в одномерном пространстве обозначаются t1, t2, t3 и.т.д. Базисная функция Bi состоит из 4-х кубических полиномов. При сшивании этих полиномов приравниваются значения функций, значения первой и второй необходимость задания неравномерной сетки возникает очень часто.

Формульный вид сплайн функций 4-го порядка для неравномерной сетки изза его громоздкости в литературе найти трудно, поэтому эти формулы (17) были выведены автором самостоятельно.

Три первых и три последних Би-сплайн базисные функции, необходимые для описания функции на краях исследуемой территории, вычисляются по отличным от (17) формулам. Эти формулы имеют более сложный вид, они тоже были выведены автором, но здесь не приводятся, хотя они не менее важны. Без крайних модифицированных Би-сплайнов аппроксимация функции скорости в ограниченной области невозможна.

При Би-сплайн параметризации производные Фреше не упрощаются, как в случае блоковой параметризации, поэтому используют простейшие численные методы для вычисления интеграла (метод прямоугольника или метод трапеции).

Недостатком этого вида параметризации является появление выбросов при резких изменениях скорости [52]. Эти выбросы отчетливо видны на Рис.

3 на глубине 45 км.

Рис. 3Параметризация с помощью Би-сплайнов слоистой скоростной функции с резким изменением Параметризация линейными Би-сплайнами используется в сейсмотомографической программе Simulps (автор С.Thurber, D.EberhartPhillips). C использованием линейных Би-сплайнов построены скоростные модели в работах [7; 8; 30]. Параметризация кубическими Би-сплайнами применена при восстановлении скоростной структуры методами сейсмической томографии в работах [33; 52].

2.2.3 Параметризация модели с помощью тетраэдров, внутри которых скорость меняется линейно непрерывно Параметризация с помощью разбиения на тетраэдры применяется и в преимущество этого вида параметризации по сравнению с другими заключается в том, что точками сетки могут быть случайным образом выбранный набор точек. Весь исследуемый объем разбивается на тетраэдры, вершинами которых являются какие-либо четыре точки из имеющегося набора. Разбивку на тетраэдры возможно осуществить большим количеством способов. В данной работе используется разбивка на пирамиды, каждая из которых удовлетворяет условию Делоне. Это условие состоит в том, чтобы задать такую разбивку на тетраэдры, при которой внутрь сферы, описанной вокруг любого построенного тетраэдра, не попадало ни одной из заданных точек [105]. Считается, что разбиение на тетраэдры, удовлетворяющее условию Делоне, является наиболее равномерным из всех возможных.

Расчет базисных функций производится по формуле [49]:

где rk-положение точки k, r1, r2, r3 – положение трех других точек тетраэдра, содержащего точки r и rk. Параметризация модели с помощью тетраэдров применена в работах И.Ю.Кулакова [93; 92], а также в [64].

2.2.4 Несколько слов о других способах параметризации модели.

ортогональных многочленов Лежандра или Чебышева [73]. Такой способ параметризации применен в [108]. При построении глобальных моделей часто используется разложение функции в конечный ряд полностью использованием этого типа параметризации построены А.Морелли и А.Дзевонски [95].

Автором работы подробно изучены 4 различных типа параметризации, а также написано несколько программ для того, чтобы использовать тот или иной вид параметризации в процессе сейсмотомографической инверсии.

Краткая информация о степени участия автора в разработке того или иного способа параметризации представлена в Табл. 1.

Табл. 1 Степень участия автора в разработке четырех различных способов параметризации модели (И – изучен; П – запрограммирован автором; С – используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора); У – усовершенствован автором; Т – использовался при тестировании.) Способ параметризации модели Степень участия автора Кубическими Би-сплайнами (2.2.2) ИПСУТ 2.3 Используемые стандартные статистические характеристики для оценки качества полученных результатов Обычная и взвешенная среднеквадратичные невязки по всем лучам RMS и RMSw рассчитываются по формулам:

где Ti – разница между реальным и расчетным временем пробега (невязка) для i-того луча, N – количество лучей, Wi – итоговый весовой коэффициент, W1i – весовой коэффициент, зависящий от качества вступления волны, W2i весовой коэффициент, зависящий от величины невязки Ti, W3i - весовой коэффициент, зависящий от эпицентрального расстояния, W4i – весовой коэффициент, зависящий от качества регистрации на станции.

Обычная и взвешенная среднеквадратичные невязки для одного события RMSсоб и RMSwсоб вычисляются по формулам:

где Nk – количество лучей, зарегистрированных от k-го события, остальные обозначения такие же как и для (19).

Количество степеней свободы системы из N уравнений с M неизвестными равно:

Вариация данных вычисляется по формуле (22) и имеет размерность сек :

Вариация параметров модели вычисляется по формуле (23) и имеет размерность (км/сек)2:

где Vi – найденная i-ая поправка к первоначальному значению скорости, К – количество неизвестных скоростных параметров.

2.4 Алгоритмы построения сейсмических лучей 2.4.1 Построение лучей в сферически-симметричной среде (одномерный луч) Для однородной слоистой среды траекторией луча является ломаная линия, состоящая из отрезков прямых линий, соответствующих каждому слою. На границах слоев луч преломляется по закону Снеллиуса где - угол между лучом и радиусом, r- радиус Земли, с – скорость сейсмической волны в слое. Чтобы упростить расчеты используют плоскопараллельное преобразование для перехода из сферической системы координат в декартову [68]. Rz и vсфvпл переводятся по формулам:

Если задана непрерывная зависимость скорости от глубины, тогда для расчета траектории луча и времени пробега в сферической среде используют формулы где =r/c, p=rsin/c – лучевой параметр, r - радиус Земли, (p)-угол в радианах. Зависимости (26) выводятся из закона Снеллиуса в предположении бесконечно тонких слоев.

телесейсмической томографии при построении лучей в стандартных одномерных скоростных моделях – IASPEI91, PREM, Herrin model.

Например, в программе С.Рокера Sphypit90 (подпрограмма delcal2.f). Также построение лучей в горизонтально-неоднородных средах в локальной томографии иногда сводится к формулам (24) и (25). Например, метод, описанный в [56].

2.4.2 Построение лучей в горизонтально-неоднородных средах В неоднородных средах фронт (x,y,z) сейсмической волны в лучевом приближении в некоторый фиксированный момент времени t описывается уравнением эйконала:

2.4.2.1 Расчет лучевых трасс аналитически 1.Для некоторых функций v [72; 114] уравнение (27) решается аналитически. В [72] в аналитической форме приведена запись решения уравнения (27) для скоростных функций следующего вида:

Если скоростная функция представима в одном из видов (28), тогда фронт волны, распространяющейся от точечного источника с координатами представляет собой дугу окружности. В [114] представлено аналитическое решение уравнения (27) для функции v вида (29).

2.4.2.2 Расчет лучевых трасс численно. Методы изгиба траектории 1. Численный метод, предложенный С.H. Thurber и W.L.Ellsworth в 1980 году, сводится к построению лучей в одномерной слоистой среде. Суть метода изложена в [56], а также на русском языке в [90]. В трехмерной вычисляется своя специальная усредненная слоистая скоростная модель и затем рассчитывается луч в слоистой скоростной модели. Данный метод использования метода состоит в том, что его возможно использовать только в случае блоковой параметризации скоростной функции или функции горизонтально неоднородном скоростном поле траектория сейсмического луча, рассчитанная эти методом, является плоской, т.е. существует плоскость, которой принадлежат все точки рассчитанной траектории. С другой стороны скорость работы алгоритма очень высокая и даже если ПК Athlon 1,33ГГц и ОП 256 МБ занимает меньше минуты.

2. Метод окружностей (approximate ray tracing ART). Суть метода описана в работе [57]. Использовался в [108] для вычисления трехмерной скоростной модели коры Тянь-Шаня. Этот метод также запрограммирован в приближенной лучевой траектории. В работе не тестировался. Все точки результирующей траектории также как и в предыдущем методе принадлежат одной плоскости.

3. Метод изгиба траектории был предложен в работе [59], он широко распространен благодаря своей простоте и быстроте вычислений. Суть метода на русском языке также подробно описывается в [90]. Этот алгоритм используется при построении лучей в программе Simulps (подпрограмма minima.f). Также автором написана собственная программа для расчета сейсмотомографических моделей получено с использованием именно обсуждаемого подхода для вычисления траекторий сейсмических лучей, например в [7]. В отличие от методов 1 и 2, рассчитанная траектория сейсмического луча является трехмерной. Необходимым условием метода, описанного в [59] является существование первой производной функции скорости в любой точке пространства, что ограничивает множество скоростных полей, в которых можно использовать этот метод.

использоваться его две различные модификации: М3-1 и М3-2. В М3-1 в качестве нулевого приближения при построении лучей используется прямая линия. В М3-2 в качестве нулевого приближения используется траектория, рассчитанная с использованием метода, описанного в параграфе 2.4.2. (пункте 2).

В 1998 году K.Koketsu и S.Sekino [23] метод был усовершенствован.

Появилась возможность использовать данный алгоритм для построения лучей не в декартовых, а в сферических координатах, а также в моделях, где функция скорости в каких-то точках претерпевает разрыв, т.е. в этих точках первая производная функции скорости не определена. В этом обновленном варианте метод использовался в работе [14].

Общий недостаток, свойственный группе методов изгиба траектории, заключается в недостаточной чувствительности этих методов к резким изменениям скорости. Если в модели имеется скоростное включение, существенно изменяющее ход луча, при использовании метода изгиба существует вероятность того, что это резкое изменение скорости не отразится в траектории. Частичная демонстрация этого отрицательного свойства методов изгиба приведена в следующей главе.

2.4.2.3 Расчет лучевых трасс численно. Методы пристрелки (shooting) 4. Суть метода Перейры для расчета траекторий сейсмических лучей [39] состоит в численном решении системы дифференциальных уравнений первого порядка с нелинейными граничными условиями. Теория метода описана в [90]. Расчеты с использованием этого метода в работе не проводились. Существенным недостатком по сравнению с 1-3 и 5 является большая громоздкость вычислений.

5. Суть метода Рунге-Кутта описывается в [68]. В этом методе для системы дифференциальных уравнений вида (30) в которой (m1,m2,m3) –вектор медленности, (x1(s),x2(s),x3(s))- траектория луча, заданная параметрически, s - длина дуги луча, остальные обозначения такие же как и в (27), составляется некоторая конечно-разностная схема. Система дифференциальных уравнений (30) выводится из (27) [114]. Вообще метод Рунге-Кутта n-го порядка предполагает ошибку ограничения порядка hn+1, в работе поиск решения производился методом Рунге-Кутта второго порядка.

Сначала с помощью метода Рунге-Кутта рассчитывались траектории трех сейсмических лучей, вышедших под разными углами из источника. Зная точки выхода этих трех лучей на поверхность, можно вычислить частные производные функции изменения местоположения конечных точек луча на поверхности от углов выхода и найти поправки к углам выхода.

Существенное преимущество этого метода состоит в том, что его легко можно запрограммировать. В процессе работы этот метод использовался для проверки правильности работы других алгоритмов. Недостаток в том, что вычисления занимают существенно больше времени, чем при проведении расчетов с использованием всех остальных алгоритмов.

6. Параксиальный метод построения лучей был изучен и применялся в расчетах в ходе работы. Основы метода изложены в [619]. Программная реализация этого метода выполнена J.Viriux. F.Haslinger модифицировал параксиальный метод построения лучей также мог быть применен [17].

осуществлять параметризацию функции квадрата медленности с помощью набора кубических Би-сплайнов, т.к. для использования этого метода медленности в любой точке пространства. В процессе работы алгоритма источника и приходящего в какую-то точку на поверхности Земли. Затем корректируются углы выхода луча по отношению к отсчетным углам и рассчитываются поправки для всех точек траектории отсчетного луча. Далее производится проверка насколько точка выхода скорректированной траектории на поверхность ближе к приемнику, чем точка выхода, которая была до коррекции. Недостаток метода 6 в том, что в любой точке пространства необходимо знание не только первой, но и второй производной заданной функции скорости или квадрата медленности.

В методах 4 и 5 искомая траектория сейсмического луча состоит из составляется своя система дифференциальных уравнений, набор этих систем дифференциальных уравнений в 4 решается одновременно, а в последовательно. При использовании параксиального метода построения луча нет необходимости решать большое количество систем дифференциальных уравнений, т.к. выведены формулы, по которым, зная пространственную поправку к любой точке траектории отсчетного луча.

Благодаря этим формулам существенно сокращается время вычислений.

2.4.2.4 Расчет лучевых трасс численно. Методы конечных разностей 7. В вышеописанных методах 1-6 сначала производится построение сейсмических лучей от источника к приемнику, а затем расчет времени пробега вдоль этих лучей. В методах 7 и 8 используется другая тактика.

Сначала определяется время пробега от источника до любой точки в заданном объеме. Зная время пробега в каждой точке пространства, появляется возможность построить поверхность с одинаковым временем пробега (волновой фронт). Затем можно вычислить траекторию сейсмического луча, руководствуясь правилом, что сейсмический луч в каждый момент времени перпендикулярен волновому фронту.

Для того, чтобы вычислить время пробега от заданного источника до любой точки в пространстве, сначала весь исследуемый объем разбивается на достаточно мелкие блоки-кубики. Потом вычисляется значение медленности, характерное для каждого блока. Затем, представляя частные производные в (27) в виде конечных разностей, вычисляется время пробега в одной из вершин куба, зная времена пробега в других вершинах.

Для большинства представленных в данной работе скоростных моделей для расчета невязки и траектории лучей используется именно метод конечных разностей, а точнее программа time3d.c, написанная P.Podvin и I.Lecomte в 1991 году [41]. Согласно [64] также существует другая программная реализация метода конечных разностей, авторами которой являются Hole и Zelt (1995) (модифицирована Vidale (1990)).

Ошибка вычисления времени пробега, которая возникает в процессе применения метода конечных разностей, главным образом связана с размером блоков, на которые разбивается объем. Чем блоки мельче, тем ошибка меньше. С другой стороны при увеличении количества блоков требуется больше машинного времени для расчета поля времен пробега для данного источника и больше дискового пространства для хранения массивов.

8. Метод расчета сейсмического луча ненулевой толщины (fat ray).

Вычисляется множество F тех точек, которые при распространении сейсмической волны от источника, находятся внутри первой зоны Френеля.

Первой зоной Френеля называется область, в которой разность хода волн меньше чем половина длины волны. Согласно [19; 19] точка х принадлежит множеству F, если она удовлетворяет неравенству (31).

где Т – период сейсмической волны, tsr – время пробега вдоль сейсмического луча, tsx – время пробега от источника к точке x, trx – время пробега от приемника к точке x. Легко проверить какие точки принадлежат, а какие не принадлежат множеству F, если знать для всего исследуемого объема поле времен пробега, соответствующих данному источнику и поле времен пробега, соответствующих данному приемнику. Чтобы вычислить поля времен пробега, в данной работе применен метод [41], изложенной пунктом ранее.

Количество вычислений, производимых с использованием методов 7- значительно больше, чем с использованием методов 1-6. В этом заключается основная трудность работы с этими методами. Заметим дополнительно, что для построения «толстого» луча (метод 8) требуется приблизительно в два раза больше машинного времени, чем при построении луча методом 7. Т.о.

для проведения быстрых оценочных вычислений методы 7-8 не подходят.

Они нужны при проведении детального исследования высокой точности.

Большим плюсом методов 7 и 8 является также возможность использовать их в скоростных полях, которые имеют разрывы первого рода. Возможность вычислить достаточно быстро с помощью метода 7 поле времен пробега от данного источника для всех точек сетки очень важна при определении координат и времени в источнике, т.к. позволяет максимальным образом уменьшить значения RMSсоб или RMSwсоб (20). Ниже в Табл. 2 представлена степень участия автора в развитии восьми разных методов построения лучей.

Табл. 2 Степень участия автора в развитии восьми разных методов построения лучей (И – изучен; П – запрограммирован автором; С – используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора); Т – использовался при тестировании.) Название метода построения лучей Степень участия автора 2.5 Методы обращения матриц. Вычисление матрицы разрешения. Решение обратной задачи Матрицу (10) иначе можно переписать в виде (32). Смысл обозначений такой же, как и в (10).

Первые М строк матрицы А` представляют собой разреженную матрицу, остальные 4*K строк имеют вид блочно диагональной матрицы.

Матрица А` является переопределенной, т.е. N>M+4K. Т.о. точного решения системы (32) не существуют. Чаще всего решение системы (32) ищется по норме L2, т.е. находится такой вектор x из всех возможных, который удовлетворяет условию:

Из-за наличия ошибок в данных и ошибок, возникших при решении прямой задачи, удовлетворяющие условию (33) поправки к значениям скоростей могут во много раз превосходить начальные значения скоростей, а поправки к параметрам гипоцентров могут быть значительно больше реально допустимых. Поэтому были разработаны методы для учета имеющейся информации об ошибках в данных и сведений о скоростных свойствах в исследуемом регионе, а также методы для нахождения сглаженного решения, т.е. решения в котором все искомые поправки значительно меньше самих параметров, к которым они относятся.

Учет ошибок во входных данных становится возможен с введением ковариационной матрицы данных Сd, которая для (32) имеет размер (M+4K)(M+4K) и умножается слева на матрицу A` и вектор невязок T.

Процедура умножения на матрицу Сd иначе называют взвешиванием каждого уравнения в системе. Если матрица Сd является диагональной, то из правил матричного умножения следует, что произведение Сd на A` эквивалентно умножению всей строки матрицы A` на соответствующий диагональный элемент. Т.о. если известно, что точность данного измерения низкая, то все члены уравнения, соответствующего этому измерению домножаются на некоторое число и это уравнение, влияет меньше на общий ход решения, чем остальные. Более подробно различные способы взвешивания, применяемые в сейсмической томографии, описаны в разделе 2.7. В данной работе используется ковариационная матрица данных только диагонального вида.

Также существует ковариационная матрица параметров модели Сm, которая имеет размер NN. Эта матрица умножается справа на вектор данных. Если матрица Сm имеет диагональный вид, то умножение справа эквивалентно умножению целого столбца матрицы A` на соответствующий диагональный элемент. В работе используется только диагональная ковариационная матрица параметров модели.

модификаций системы (32) для получения более устойчивого и гладкого решения имеем систему (34):

Условие для нахождения решения z по норме L2 теперь с использованием обозначений из (34) выглядит следующим образом:

демфирующая матрица D1 является диагональной (36) и имеет размер (M+4K)(M+4K).

где 1,, М+4K – некоторые заданные числа, которые называются демферами, они могут быть равны или не равны между собой. Введение этой матрицы позволяет найти решение незначительно отличающееся от нуля.

Смысл введения второй демфирующей матрицы D2 состоит в том, чтобы возмущениями, соответствующими двум соседним блокам или двум соседним точкам сетки. Эта матрица является переопределенной и имеет следующий вид:

где µ – демфирующий коэффициент, 1 и k, 1 и l, j и l, k и l – номера соседних блоков или соседних точек сетки.

После введения демфирующих матриц имеет смысл вводить параметр, оценивающий насколько полученное решение близко к идеальному по норме L2, т.е. к решению, получаемому при выполнении условия (33). Такая оценка обеспечивается матрицей разрешения. Если решение уравнения (32) обозначить xист, а решение уравнения (34) xвыч, тогда для матрицы разрешения R будет верно равенство:

При нахождении решения (34) методом нормальных уравнений (см. 2.5.1) или с помощью сингулярного разложения (см.2.5.2) xвыч представимо в виде где F – некоторая матрица. С использованием (39) для матрицы R будет верно равенство:

Матрица разрешения R имеет размер (M+4K)(M+4K), значение каждого элемента матрицы R всегда лежит на отрезке [0,1]. Чтобы информацию, содержащуюся в матрице разрешения представить наглядно, строят карты диагональных элементов матрицы разрешения. Другой способ представления информации, скрытой в матрице R, – это расчет спрединговой функции [58;

33].

использованием (39) выражение для матрицы ковариации cov имеет вид где b – вектор ошибок, связанный с неточностями во входных данных.

Размер матрицы ковариации такой же, как и размер матрицы разрешения.

Ниже в 6.1.2 и в 6.2.2 построены карты стандартной ошибки полученного решения dVp/Vp, Vp/Vs. Под стандартной ошибкой в данном случае понимается квадратный корень диагональных элементов матрицы ковариации. При расчете стандартной ошибки точность определения времени вступления для всех лучей считается равной 0.1 сек. Т.о., в общем виде вычисляется по формуле:

Из формул (32) и (10) видно, что первые M столбцов матрицы A` образуют разреженную матрицу A, а последние 4*K столбца образуют блочно-диагональную матрицу, каждый блок которой представляет собой матрицу Hn.

Каждая матрица H n состоит из 4 столбцов и Nk строк (Nk – количество лучей от k-го землетрясения). Согласно [94] каждая матрица H n может быть сведена к верхней треугольной путем домножения H n на некоторую ортогональную матрицу Qn. Т.к. размер матрицы H n 4Nk, то, следовательно, после приведения ее к треугольному виду только 4 строки будут ненулевыми и Nk-4 будут нулевыми. Домножим первые N1 строк матрицы B на соответствующую матрицу Q1, следующие N2 строк на Q2 и.т.д. В результате получится некоторая матрица A, у которой в N-4*K строках элементы, соответствующие переменным hn, равны нулю. С помощью остальных 4*K строк матрицы A переменные hn определяются однозначно. Т.о. чтобы найти решение x, сначала необходимо обратить матрицу размером M(NK) и найти вектор, а затем последовательно, уже зная, вычислить для всех землетрясений (для всех n). Процедура нахождения ортогональных матриц Qn и расчет матрицы A далее в работе будет называться разделением переменных [38]. При использовании разделения переменных система (32) преобразуется в размер A M(N-4*K), размер r N-4*K. Ковариационные матрицы данных и параметров модели, а также демфирующие матрицы для уравнения (44) вводятся аналогичным образом, также как и для (32). Все приведенные ниже формулы соответствуют случаю, когда разделение переменных не производилось, т.е. решалась система уравнений (32). Аналогичные формулы легко могут быть выведены и для системы уравнений (44), т.е. случая, когда разделение переменных имело место. Метод разделения переменных [38] запрограммирован в программах Sphypit90 [43] (подпрограмма buildg2.f и u2.f), Simulps (автор С. Thurber, подпрограмма parsep.for). Этот метод применялся при расчете скоростных моделей в [44; 7] и в других работах.

2.5.1 Нахождение решения уравнения (34) с помощью построения системы нормальных уравнений Как было замечено выше, чтобы найти решение системы (32) по норме L2, необходимо найти такой вектор x, который удовлетворяет условию (33).

Дифференцируя условие (33) по xk получаем систему уравнений, которая называется системой нормальных уравнений.

В системе (45) количество неизвестных равно количеству переменных и равно M+4K. Эта система имеет единственное решение:

Матрица A`+ имеет M+4K строк и N столбцов. Эта матрица называется псевдообратной матрицей или матрицей Мура-Пенроуза.

Для модифицированного условия (35), если использовать обозначения, введенные в (34), система нормальных уравнений имеет вид:

Псевдообратная матрица A`+ с использованием обозначений (34) равна а вектор решения x с использованием обозначений (34) и (48) Из (38), (45), (34) и (48) следует, что матрица разрешения равна:

Из (48)и (49) видно, что решение x можно вычислить, только зная обратную матрицу к B T B + D T D. Для обращения этой матрицы возможно использовать различные методы. В рамках данной работы применялось разложение Холесского. Суть этого разложения в следующем.

Какой бы ни была матрица B (34) в любом случае матрица B T B + D T D является квадратной симметричной и неотрицательно определенной. Для любой квадратной симметричной и неотрицательно определенной матрицы P согласно [94; 87] существует верхняя треугольная матрица U, такая что верно равенство Если, например, представить матрицу B T B + D T D из (47) в виде (51):

то вектор z из (47) можно получить, решая две треугольные системы:

Матрицу разрешения (50) можно посчитать с помощью решения M+4K уравнений таким двухэтапным методом.

Если обозначить элементы матрицы U uij, а элементы матрицы P из (51) pij, тогда практически uij вычисляются по формуле из [94]:

В формуле (54) сумма считается равной нулю при i=1.

При использовании этого метода точное решение, удовлетворяющее условию (35), находится за конечное число шагов.

2.5.2 Нахождение решения уравнения (34) с помощью сингулярного разложения матрицы B В данной работе этот метод обращения матрицы использовался только для случая, когда где I – единичная матрица. Из курса матричного анализа известна следующая теорема [76]:

Для любой RT (R>T) матрицы S существует ортогональная RR матрица U, ортогональная TT матрица V и RT-диагональная матрица с диагональными элементами 12……0, такие что выполняется равенство (56).

Ортогональной называется такая матрица, у которой транспонированная и обратная матрицы равны между собой. Элементы матрицы 1, 2 и.т.д.

называются сингулярными числами. Если известно сингулярное разложение матрицы S и S – матрица полного ранга, найти обратную просто по формуле:

Т.о. представляя матрицу B в виде (55) и используя (56) уравнение (47) можно переписать в следующем виде:

Из (58) легко вывести формулу для вычисления z, а также с учетом (38) формулу для вычисления матрицы разрешения.

программный код R.C. Singleton, созданный в 1973 г (Simulps, подпрограмма fksvd.for). В указанной подпрограмме матрица B сначала преобразуется к верхней двудиагональной посредством последовательных умножений B справа на специально вычисленные матрицы. А затем производится сингулярное разложение полученной двудиагональной матрицы. [94] Также как и для предыдущего способа обращения в данном случае точное решение, удовлетворяющее условию (35) находится за конечное число шагов.

2.5.3 Нахождение решения уравнения (34) методом LSQR Этот способ решения системы (34) также используется только для случая, когда демфирующая матрица D представима в виде (55): Данный метод решения уравнения (34) является итерационным, т.е. решение ищется с помощью серии последовательных итераций. Итерации продолжаются до тех пор, пока норма L2 не будет меньше некоторого числа. Следовательно, решение, найденное с помощью этого метода является приближенным в отличие от двух предыдущих методов. В описываемом алгоритме в явном виде псевдообратная матрица A`+ не вычисляется, а сразу находится вектор решения. Т.о. невозможно вычислить матрицу разрешения R и матрицу ковариации cov, как это можно было сделать при использовании двух предыдущих методов. Этот недостаток описываемого метода хорошо известен и существуют работы [35], в которых предлагаются различные алгоритмы для приближенного вычисления матрицы разрешения. С другой стороны очень эффективным является использование метода LSQR для разреженных переопределенных матриц больших размеров (более неизвестных), т.к. скорость работы этого алгоритма при обращении таких матриц существенно выше, чем скорость работы двух предыдущих алгоритмов, а необходимый для хранения массивов объем памяти значительно меньше. Для получения решения этим методом используется программа, написанная C.C. Paige и M.A. Saunders в 1982 году. В этой программе согласно [53; 36] для нахождения решения на каждой итерации сначала необходимо вычислить некоторую систему базисных векторов.

Затем на p-ой итерации найти решение по норме L2 простейшей трехдиагональной системы из p+1 уравнений с p неизвестными, а затем провести еще одну процедуру умножения матрицы на вектор.

Табл. 3 Степень участия автора в развитии трех различных методов обращения матриц (И – изучен;

С – используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора); Т – использовался при тестировании.) Название метода обращения матрицы Степень участия автора

LSQR ИСТ

Теоретические основы решения обратных задач в геофизике вообще изложены в работе [121].

2.6 Методы локации землетрясений в одномерных скоростных моделях При определении координат вводится упрощающее предположение, что землетрясение – это точечный источник возбуждения поперечных и продольных сейсмических волн. Выше в 2.1 и далее в работе землетрясение рассматривается именно как точечный источник возбуждения сейсмических волн.

Задачи определения параметров гипоцентра и одномерной или трехмерной модели скоростей сейсмических волн в Земле, как видно из 2.1, являются взаимосвязанными и взаимозависимыми. Если неизвестны скорости сейсмических волн в исследуемом регионе, то нельзя с высокой точностью определить координаты землетрясения. Наоборот, если координаты землетрясений определены достаточно грубо, то траектории сейсмических лучей являются неверными и полученная скоростная модель абсолютно не соответствует реально существующей. Далее речь пойдет о методах, применяемых для локации событий в одномерной и трехмерной скоростных моделях. Более детально будут освещены методы, используемые автором в процессе работы по локации событий.

1. Самый простой и грубый способ определения эпицентрального расстояния до землетрясения по данным времен пробега волны Р и S, зарегистрированным на одной станции, изложен в [99]. Он заключается в вычислении эпицентрального расстояния по формуле:

где tP и tS – это время пробега волны Pи S соответственно. Если станция является трехкомпонентной, то азимут на событие возможно определить по отношению смещений продольной волны, зарегистрированной двумя перпендикулярными горизонтальными сейсмографами. Эпицентральное расстояние и азимут однозначно определяют точку, которая является гипоцентром данного землетрясения.

2. Большинство программ локации событий (Hypo71, Hypocenter и др.) основаны на использовании метода Гейгера, предложенного в 1910 году.

Например, в программе Hypoellipse сначала приближенные координаты событий и время в источнике определяются в однородной скоростной модели (скорость во всем пространстве предполагается равной V0) с помощью системы из К нелинейных уравнений вида:

где К – количество измерений одного события на разных станциях, x,y,z,t – неизвестные координаты события и время в источнике, xi,yi,zi,ti – координаты станции, регистрирующие i-ую волну и время вступления i-ой волны на станции (могут быть использованы различные типы волн). Эта система последовательного вычитания всех нижних уравнений из первого. Далее поправки к приближенно определенным координатам находятся с помощью системы уравнений (62).

где Nk – количество лучей от одного события на разные станции, Ti разница между реальным и расчетным временем пробега для i-го луча,, x, y, z – поправки к приближенно определенным координатам события и времени в источнике. Системы уравнений (61) и (62) чаще всего являются переопределенным, поэтому решение ищется по норме L2 (33) или какойнибудь другой норме. Если число измерений меньше 4-х, то неизвестные необходимо знать траектории сейсмических лучей в заданной слоистой скоростной модели. Эти производные связаны с углом выхода луча и скоростью в источнике следующим соотношением:

где i, i, i - углы между касательной к i-ому лучу в источнике и осью x, y, z соответственно, Vист – скорость сейсмической волны в источнике.

Программы, вычисляющие координаты и время в источнике по методу Гейгера, отличаются друг от друга введением дополнительных различных весовых матриц и весовых коэффициентов, различными способами обращения матриц. В ходе работы автор использовал программы Hypoellipse использованием метода Гейгера в трехмерной скоростной модели.

3. Также для локации событий применялся метод сеточного поиска. С развитием компьютерной техники, а также с появлением конечноразностных методов для вычисления времени пробега от источника до различных точек исследуемой области стало возможным не обращать матрицу H из (62), а искать абсолютный минимум обычным методом перебора. Метод перебора является более простым и надежным по сравнению с численным решением переопределенной системы линейных уравнений. Автором работы был запрограммирован алгоритм, предложенный в [34]. Локация событий производится как в одномерной, так и в трехмерной скоростной модели. Вычисленные функции RMSсоб или RMSwсоб после проведения локации этим методом действительно имеют меньшие значения, чем после проведения локации методом 1.

4. Также заслуживает внимания метод двойных разностей для определения локации событий, т.к. в некоторых случаях он является более эффективным, чем описанные выше методы. Он был предложен в 2000 году F. Waldhauser и W.Ellsworth [62]. Метод основан на решении линейной системы, каждое уравнение в которой имеет вид:

i – номер одного события, j – номер другого события, k – порядковый номер уравнения, drijk – двойная разность, остальные обозначения аналогичны (62).

Т.о. в предыдущих методах производился поиск таких значений x,y,z и t, которые уменьшают максимальным образом функцию RMSсоб или RMSwсоб, в этом методе предлагается минимизировать функцию всех возможных комбинаций двойных разностей. Этот метод позволяет производить более точную локацию событий, локализованных в одной области (кластере). Чаще всего наибольшее скопление землетрясений наблюдается в зонах разломов. В [62] методом двойных разностей произведена более точная локация событий, произошедших в районе Северного Хауварского (Hayward) разлома в землетрясения также образуют кластер. В работе [47] произведена перелокация событий методом двойных разностей для землетрясения, произошедшего 26 июля 2001 года на острове Скирос (Греция). При использовании метода двойных разностей контуры разломных зон или областей, в которых произошла разрядка накопленного за некоторый период геологического времени напряжения, удается выявить более точно. Заметим, что абсолютная точность локации всей этой области или группы событий попрежнему зависит от количества станций, регистрирующих его.

Преимущества этого метода перед более ранними в случаях применения в областях наибольшего скопления событий заключаются также в том, что 1. исчезает ошибка, связанная с неоднородной скоростной структурой непосредственно под станцией 2. локация не зависит от скоростных неоднородностей на пути от станции к паре событий, а зависит только от скоростных неоднородностей непосредственно вблизи перелоцируемых событий.

Авторами статьи [62] также была написана программа HypoDD для локации событий с использованием предложенного ими метода. В ходе данного исследования программа HypoDD была усовершенствована для случая локации событий в трехмерной среде с помощью построения лучей в скоростной модели с горизонтальными неоднородностями. Анализ сейсмичности территории Тянь-Шаня показал, что в этом районе существуют зоны, в которых плотность событий в несколько раз больше, чем во всем остальном пространстве. Для некоторых из этих зон в настоящей работе была произведена попытка локации событий методом двойных разностей. Ниже в Табл. 4 представлена степень участия автора в развитии трех различных методов локации событий.

Табл. 4 Степень участия автора в развитии трех различных методов локации событий (И – изучен; П – запрограммирован автором; С – используется при расчетах реальных скоростных моделей (может быть использован и не в программной реализации автора); У – усовершенствован автором; Т – использовался при тестировании.) Название метода локации событий Степень участия автора 2.7 Различные способы введения весовых коэффициентов Введение весовых коэффициентов является неотъемлемой частью любой томографической программы. Каждой станции может быть присвоен свой вес в зависимости от качества регистрации вступлений на этой станции, каждому событию может быть присвоен также свой весовой коэффициент в зависимости от качества вступлений сейсмических волн, пришедших только от этого события. Более распространено присвоение собственного веса каждому лучу в зависимости от того, насколько точно определено время определения времени вступления волны t. Наиболее часто используемым способом взвешивания является разделение всех лучей на весовые группы.

Например, если t0.07сек, то целочисленный вес этого луча wt=0, если 0.07~540 млн. лет) на поверхность [119]. Считается, что в эпоху каледонской складчатости (периоды ордовик, силур, ~500-410 млн. лет назад) произошло образование Северо-Тянь-Шаньской складчатой зоны, ограниченной с севера СевероТянь-Шаньским глубинным разломом, с юга линией Николаева (ТерскейКаратауским разломом), с запада Талассо-Ферганским (Рис. 18). Узкая область Срединного Тянь-Шаня, оконтуренная с разных сторон линией Николаева, Талассо-Ферганским и Атбаши-Иныльчекским разломами образовалась в эпоху Каледонской и Герцинской складчатости (~300- Млн. лет). Южнее до Гиссаро-Кокшаальского глубинного разлома расположена герцинская складчатая зона Кокшаал (~410-300 Млн. лет).

Рис. 18 Тектоническая карта из [100]: 1 – Таласо-Ферганский, 2 – система Северо-Тянь-Шаньских разломов, 3 – Гиссаро-Кокшаальский, ;4 – Атбаши-Иныльчекский, 5 – Линия Николаева.

К юго-западу от Талассо-Ферганского разлома основные системы хребтов Чаткало-Кураминский, Гиссаро-Алайский и Фергано-Атойнакский (Рис. 19) образовались частично в каледонско-герцинский, а в основном в герцинский этап.

горообразовательных процессов на территории Тянь-Шаня скорее всего не происходило. Исследование различных геофизических полей, которые кратко изложены в последующих параграфах, выявили ряд признаков (утончение коры, высокий тепловой поток, отрицательные изостатические рифтогенеза на юго-востоке и на севере Тянь-Шаня [119].

Рис. 19 На карту рельефа местности схематично нанесено местоположение некоторых впадин и горных хребтов Тянь-Шаня. Черными номерами обозначены впадины: 1 – Баткенская, 2 – Ферганская; 3 – Ноокатская; 4 – Чуйская; 5 – Нарынская; 6 – Суусамырская; 7 – Кочкорская; 8 – Джумгальская; 9 – Верхнее-Нарынская; 10 – Иссык-Кульская; 11 – Атбашинская; 12 – Алайская; 13 – Сарыджазская. Желтыми линиями и номерами обозначены хребты: 1 – Киргизский; 2 – Таласский Алатау; 3 – Ферганский; 4 – Кунгей Алатау; 5 – Заилийский Алатау; 6 – Джумгалтау;

7 – Кокшаал Тау; 8 – Алайский; 9 – Атбаши; 10 – Чаткальский; 11 – Атойнакский.

В работе [84], исследованы петрохимические и изотопно-геохимические данные Северного Тянь-Шаня, выдвинуто предположение, что 55 Млн. лет назад Тянь-Шань находился в условиях предрифтового режима под влиянием мантийного плюма, но дальнейшее его перерождение было приостановлено, начавшейся 45 Млн. лет назад коллизией Индостанской и Евразийской плит литосферы. Из работ [119; 100] видно, что под Нарынской впадиной происходит выгибание поверхности Мохо (кора утоньшается). С помощью исследований методом функции приемника [24] установлено, что в этой зоне коро-мантийный переход размытый, скорости поперечных волн пониженные.

Эти результаты говорят о возможном расположении в прошлом под впадиной горячей точки [78]. По мнению В.Н.Погребного и Т.М.Сабитовой [97] современная «динамика региона Высокой Азии определяется сочетанием процессов, связанных как со столкновением Индостанской и Евразийской литосферных плит, так и продолжающимся развитием плюма».

С эпохой альпийской складчатости связано появление следующих межгорных впадин справа от ТФР на севере: Чуйской, Суусамырской, Кочкорской, Джумгальской, на юге – Нарынской и Атбашинской, слева от ТФР – Ферганской и Алайской [119] (Рис. 19).

Далее в параграфах 4.4, 4.6, 4.7 также отмечается существенно разный характер аномалий ряда физических величин по обе стороны ТФР, что свидетельствует о различной глубинной природе геодинамических процессов в этих районах.

4.3 Магнитные аномалии.

Как и в большинстве районов мира, в Высокой Азии положительные магнитные аномалии T, вызванные магнитными неоднородностями в верхней мантии, наблюдаются в асейсмичных областях, т.е. на платформах:

Тибетской, Туранской, Индостанской [98], а отрицательные T – в высокосейсмичных зонах. Для исследования коры территории Тянь-Шаня имеет смысл рассматривать магнитные аномалии, связанные с магнитными неоднородностями на более мелких глубинах. Расчет магнитных аномалий T1, связанных с магнитными массами на глубинах 10-30 км показал, что отчетливо прослеживается разный характер аномалий T1 по обе стороны Северо-Тянь-Шаньского разлома. Севернее разлома аномалии положительные полосовидной формы, южнее наблюдаются отрицательные аномалии разнообразной формы. Слой пород, сохранивших магнитные свойства под впадинами имеет существенно большую мощность, чем под хребтами, поэтому положительные аномалии также соответствуют центральной и западной части оз. Иссык-Куль и центральной и западной частям Ферганской впадины.

4.4 Аномалии силы тяжести, изостатическое состояние земной коры, плотность пород.

Для всей территории Высокой Азии региональные минимумы силы тяжести характерны для горных систем: Тибета, Тянь-Шаня, Памира, Гиндукуша, в то время как с платформенными образованиями: Турнская и Таримская плита, Казахский и Индостанский щит, - связаны положительные региональные аномалии. Для Тянь-Шаня региональное гравитационное поле характеризуется общим понижением значений с север-северо-запада на югюго-восток. Это поле отражает поведение поверхности Мохоровичича и плотностные неоднородности верхней мантии. Для изучения плотностных неоднородностей в коре исследуются остаточные гравитационные поля второго и третьего порядка и их горизонтальные градиенты. Зоны высоких градиентов в основном имеют северо-восточное, северо-западное и субширотное направление. В.И.Кнауфом и К.Д.Помазковым было отмечено, что по этим же направлениям расположены зоны крупных разрывных нарушений. Зона высокого горизонтального градиента поля первого порядка на севере соответствует Северо-Тяньшаньскому разлому. Линия Николаева и Таласо-Ферганский разлом проявляются при исследовании зон высоких градиентов второго порядка.

Исследованию изостазии Тянь-Шаня посвящены работы [70; 83; 120] и многие другие. Артемьевым М.Е., Голландом В.Э. и Юдахиным Ф.Н.

установлено, что за счет вариаций мощности земной коры удается компенсировать 53-66% избыточной поверхностной нагрузки, оставшаяся часть компенсируется за счет изменения мощности разуплотненной мантии под корой до глубин ~80-90 км. Размер изостатических аномалий на ТяньШане ~ в 2-5 раз меньше, чем в горных системах, появившихся в период альпийской складчатости (например, на Кавказе.) На юго-востоке орогена интерес представляет узкая полоса отрицательных изостатических аномалий.

Эта полоса включает в себя Нарынскую, Верхнее-Нарынскую и Сарыджазскую впадины (Рис. 19). По гравитационным данным и данным времен пробега сейсмических волн под этими впадинами наблюдается существенное сокращение мощности земной коры и разуплотнение верхней мантии. В книге [119] также отмечается, что к юго-западу от ТалассоФерганского разлома в земной коре наблюдаются явные нарушения изостазии, в то время как к северо-востоку земная кора близка к состояние изостатического равновесия.

Для самого нижнего структурного этажа, сложенного архейпротерозойскими образованиями, среднее значение плотности около 2. г/см3 (>1650 Млн. лет). Средняя плотность пород рифей-вендского структурного этажа 2.73 г/см3 (540-1650 Млн. лет). Плотность пород, образованных в эпоху каледонской складчатости, 2.74 г/см3. Герцинскиий структурный этаж сложен породами (палеозойские породы, ~540-250 Млн.

лет) со средней плотностью 2.67 г/см3. Диапазон изменения плотности мезозойско-кайнозойских пород очень широк. Плотность варьируется от 1.5до 2.60-2.75 г/см3. В верхней части земной коры под Чуйской и ИссыкКульской впадинами плотность 2.8-2.85 г/см3 [119] (Рис. 19).

4.5 Тепловое поле Основные горные системы Альпийско-Гималайского и МонголоОхотского горных поясов выделяются повышенным значением теплового потока. Среднее значение теплового потока на Тянь-Шане на поверхности консолидированной коры, т.е. на нижней границе осадочного слоя около мВт/м2. Для сравнения тепловой поток на Памире на той же границе достигает значений 120 мВт/м2, на Кавказе 80-90 мВт/м2, в Копетдаке 70- мВт/м2, а на Казахской платформе 45 мВт/м2. Если исследовать значение теплового потока на глубине 1 км, то области повешенных значений 70- мВт/м2 находятся на севере и северо-западе исследуемого орогена, к югу тепловой поток уменьшается [118]. Выделяют две большие зоны аномально высокого теплового потока: Западно-Тянь-Шанскую (до 108 мВт/м2 на глубине 1 км) и Восточно-Тянь-Шанскую (до 134мВт/м2 на глубине 1 км), которые разделены зоной аномально низкого теплового потока – ЧуйскоНарынской (до 23 мВт/м2 на глубине 1 км) Западно Тянь-Шаньская зона состоит из Чаткала-Кураминской (70-108 мВт/м2 на глубине 1 км) аномалии на севере и Фергано-Алайской (70-75 мВт/м2 на глубине 1 км) на юге. Обе имеют северо-восточное простирание. Между этими аномалиями расположена область низкого теплового потока (до 32 мВт/м2 на глубине км) в Ферганской впадине. В Чуйско-Нарынской зоне выделяют ЮжноЧуйскую аномалию (23 мВт/м2) и Джумгал-Нарынскую (31 мВт/м2) северозападного простирания. В восточно-Тянь-Шаньской области самыми высокими значениями характеризуется Барскоонская аномалия (134 мВт/м2), занимающая центральную часть Иссык-Кульской впадины и обрамляющих их хребтов. Севернее также отмечают Кюнгейскую (70 мВт/м2) аномалию и Северо-Чуйскую (80 мВт/м2) [117].

Что касается связи зон глубинных разломов и аномалий теплового потока, то в [117] отмечается повышение теплового потока для глубинных разломов, испытывающих растяжение в новейшее время (Барскоонская аномалия). Важнейшая структурная линия Тянь-Шаня (линия Николаева) не находит отражения в тепловом поле.

Температура на поверхности Мохоровичича на территории Тянь-Шаня изменяется от 5000 до 16000С. Для центральной части орогена характерны температуры 8000-10000С. На северной и северо-западной окраинах расположены области аномально высоких температур 13000-16000С. В югозападной части исследуемой области находится область аномально низких температур на границе М 5000-6000С. На юге Восточно-Чуйской и востоке Иссык-Кульской впадины также расположены зоны аномально низких температур.

По геотермическим данным получены сведения о толщине литосферы в различных районах Тянь-Шаньского горного массива. К областям аномально низкой мощности литосферы относятся – Чаткало-Кураминская, простирающаяся вдоль северо–западной границы исследуемой территории, Северо-Чуйская, расположенная на северной окраине орогена, Кюнгейская, начинающаяся на северо-восточном побережье озера Иссык-Куль и уходящая на восток, Барскоонская, охватывающая западную часть озера Иссык-Куль и территории примыкающие к озеру с юга и Западно-Нарынская, локализованная в центральной части Тянь-Шаня. Выделены также три зоны аномально повышенной области теплового потока – Южно-Чуйская – примыкает к Северо-чуйской с юго-запада, Пржевальская – южнее Кюнгейской области, Фергано-Алайская – юго-западная часть изучаемого горного массива [117].

В заключении отметим, что в северных районах: Чуйском и ИссыкКульском проводилось гораздо больше измерений теплового потока, чем на юге, поэтому для этих областей информация является более детальной.

4.6 Активные разломы, скорости современных деформаций.

4.6.1 Скорости горизонтальных движений.

Современная средняя скорость укорочения Тянь-Шаньского горного массива вдоль направления север-юг по измерениям GPS составляет 20 мм/г [16] скорость сближения Индийской и Евроазиатской плиты около 50 мм/г.

Скорости горизонтальных и вертикальных движений, определенные по GPS, с помощью суммирования сейсмических моментов и геологическим методом, часто расходятся. Обычно геологический метод дает более низкие значения скоростей, чем два других. Ниже приведены скоростные характеристики разломов, полученные геологическим методом.

Карта активных разломов для изучаемой территории по [109] представлена на Рис. 20. Для одного из крупных разломов, ТаласоФерганского, имеющего протяженность около 400 км и северо-западное простирание, известно, что правосдвиговая компонента смещения в 10-30 раз превосходит взбросовую. По скоростям сдвиговых смещений разлом делится на 5 сегментов. Современные скорости сдвиговых перемещений и протяженности сегментов представлены в Табл. 7.

Табл. 7 Скорости сдвиговых перемещений Vсд для различных сегментов ТФР. Нумерация сегментов начинается с юга.

I II III IV V VI

Длина, км Vсд, мм/год Между I и V сегментом суммарное сдвиговое смещение за четвертичное время составило 10-12 км. Разломы Западно-Джунгарский и Джунгарский, которые не вошли в исследуемую территорию и расположены в северосеверо-восточном направлении от озера Иссык-Куль, тоже имеют северозападное простирание. Для этих разломов скорости современных правосдвиговых движений в 8-10 раз превосходят взбросовые и составляют мм/год для Западно-Джунгарского и 3-5 мм/год для Джунгарского (Рис. 20).

для Чаткальского (Пскемского) – скорость взброса 0,5 мм/год (Рис. 20).

К востоку от ТФР К.Е.Абдрахматовым выделена Нарын-Сонкульская разломная зона (Рис. 20). На южном фланге этой зоны ЗападноАкшийракский разлом, для которого скорость надвигания 2-3 мм/год, на севере Нарынский разлом, его скорость взброса 2 мм/год. Для СевероДжумгальского разлома имеющего северо-восточное направление скорость надвига до 3 мм/год. Для Южно-Сонкульского и Северо-Нарынского разломов взброс сочетается с превосходящим правым сдвигом (Рис. 20).

Скорости сдвига 1-3 мм/год и 1,5мм/год соответственно. Скорости взброса 0,1-0,5 и 0,15-0,2 мм/год.

На южном склоне Кунгей-Алатау (Рис. 20) для разлома СевероАксуйского скорость взброса 0.35-0.4 мм/год, скорость правого сдвига 0. мм/год. Для Заилийского и Алмаатинского взбросов (Рис. 20) средние скорости оценены величинами 0.5-0.8 и 0.25-0.4 мм/год 4.6.2 Скорости вертикальных движений.

В [119] отмечается, что помимо вертикальных движений, связанных с субмеридиальным сжатием, также существуют автономные вертикальные движения, которые по мнению Ф.Н.Юдахина вызваны разуплотнением мантии.

Характер вертикальных движений по обе стороны от ТалассоФерганского разлома за новейший период различны. К юго-западу от разлома общий размах движений 13-14 км (-6 км в Ферганской впадине и до +7 в Заалайском хребте), к северо-востоку 9-11 км (~-4 км в Чуйской и Иссык-Кульской впадинах и +4.5 в Киргизском хребте). Для Джунгарского разлома 8,4 мм/год.

4.7 Исследования Тянь-Шаня по сейсмологическим данным.

4.7.1 Общая информация.

Некоторые сведения о скоростях сейсмических волн для различных областей Тянь-Шаня возможно получить после анализа годографов для этих областей от землетрясений или взрывов. Первый годограф сейсмических волн для средней Азии был составлен Е.А.Розовой в 1936 г. Он состоял из трех ветвей, соответствующих определенным скоростным слоям в коре. В работе [101] этот годограф был скорректирован, выделены новые ветви, выявлены области на территории Киргизского Тянь-Шаня, где кора более неоднородна, а где менее. В книге [101] представлены горизонтальнослоистые скоростные модели для четырех различных профилей, полученные на основе построения системы встречных годографов. Эти профили пересекают всю территорию Киргизии в направлении юго-запад – северовосток и юго-восток – северо-запад. По годографам, построенным по данным площадных систем наблюдений, получена более детальная информация о скоростях верхних слоев коры до глубин 10-15км на севере Тянь-Шаня (хр.

Таласский Алатау, Киргизский хребет, Чуйская впадина), в центральной области (зона Токтогульского водохранилища и его окрестности) и в юговосточной части Ферганской впадины (Папанский р-н). Особенно ценна информация, полученная по данным от промышленных взрывов, произошедших в разные годы на территории Киргизского Тянь-Шаня.

Т.М.Сабитовой [101] проанализированы данные от Токтогульского взрыва (08.02.1975, ~730E, 420N), Актюзского (12.05.1982, севернее оз. Иссык-Куль), Торуайгырского (11.07.1974, севернее оз. Иссык-Куль), Сулюктинского (10.04.1977, горное обрамление Ферганской впадины), Наманганских (10.07.1974, 02.11.1974, 10.03.1975, ~71.50E, 410N, Узбекистан) и других. В работе выявлено, что северо-западная часть Киргизского Тянь-Шаня до глубин 10-15 км является более неоднородной, чем северо-северо-восточная.

На северо-западе в верхней коре обнаружены зоны высокоскоростных включений, а на север-северо-востоке возможно наличие зон низкоскоростных включений. Также в верхней части коры на юго-западе (юго-восточное обрамление Ферганской впадины) скорости сейсмических волн ниже, чем на севере и северо-западе. Автором предлагается четырехслойная модель коры, которая представлена в Табл. 8. Также в Табл.

8 представлена модель коры, предложенная в книге [119].

Табл. 8 Модель коры для территории Киргизского Тянь Шаня по [101] и [119]. Vp - скорость продольных сейсмичеких волн, Z - глубина залегания нижней границы.

Назв. слоя. «осадки» «гранит1» «гранит2» «базальт» Мантия Источник По разные стороны ТФР кровля «гранитного» слоя расположена на разных глубинах, к западу глубина залегания порядка 10-20 км, к востоку — 5-8 км.

По данным Т.М.Сабитовой и Е.М.Бутовской следует, что глубина залегания границы Конрада возрастает с севера на юг. На территориях, граничащих с Туранской плитой и Казахским щитом, она составляет 15-20 км, под горноскладчатыми сооружениями Тянь-Шаня возрастает до 35-40. Поверхность Мохоровичича в пределах Тянь-Шаньского блока в целом прогибается от 40км на периферии до 60-65 км внутри блока. К югу от важнейшей структурной линии Николаева мощность земной коры меньше, чем к северу.

В книге [101], а также в работах [96; 110] предполагается, что изменение мощности земной коры здесь происходит скачкообразно. В районе озера Сон-Куль разность залегания границы Мохо по разные стороны разлома возможно порядка 10-15 км, далее на северо-восток разность уменьшается до 10 км [101]. В районе Токтогульского водохранилища изменение мощности сейсмологическим данным не наблюдается.