Some problems of joint inversion of temporal gravity variations with data on land and seafloor displacements: a case study of the Tohoku-Oki earthquake of march 11, 2011

Full Text

ВВЕДЕНИЕ

Катастрофическое землетрясение Тохоку-Оки произошло 11 марта 2011 г. на участке зоны субдукции в районе северной части острова Хонсю и имело магнитуду M_W = 9.0$–$9.1. В области землетрясения много лет ведутся детальные наблюдения, развернуты плотные сейсмические и геодезические сети. Японская сеть GEONET насчитывает более 1200 станций GPS, включая 400 станций на севере о. Хонсю. Кроме того, имеется уникальная сеть донных транспондеров, впервые обеспечивших измерения вблизи глубоководного желоба, где произошли наибольшие смещения [Sato et al., 2011; Iinuma et al., 2012].

К настоящему времени построен целый ряд моделей поверхности сейсмического разрыва, определяющих ее геометрию и поле смещений (см. обзор [Tajima et al., 2013]). Эти модели базируются на различных комбинациях данных сейсмологии, ГНСС, гравиметрии, приливных станций, и поэтому полученные результаты различаются по положению области максимального смещения на поверхности разрыва и по его амплитуде. Модели, построенные на основе сейсмических волновых форм, с учетом или без учета данных GPS и волновых форм цунами, дают смещения на поверхности разрыва до 50 м (например, модель, опубликованная на сайте https://earthquake.usgs.gov/earthquakes). В работе [Pollitz et al., 2011] на основе только данных GPS максимальные смещения оценены в 33 м. Положение области максимальных смещений в различных моделях находится как ниже гипоцентра [Ammon et al., 2011; Ozawa et al., 2011; Simons et al., 2011], так и немного выше него [Pollitz et al., 2011; Shao et al., 2011] или же на большом удалении, в области океанического желоба [Lay et al., 2011; Iinuma et al., 2012]. В данной работе на основе анализа различных данных по землетрясению Тохоку-Оки проводится анализ разрешающей способности и согласованности информации, полученной различными методами.

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

В нашем исследовании мы использовали четыре группы данных. Это смещения на стационарных пунктах GPS, входящих в японскую систему GEONET; смещения отражающих площадок на земной поверхности и техногенных объектах, определенные по радарным снимкам спутников ENVISAT; смещения донных транспондеров, по результатам повторных локаций их положения с надводных кораблей [Sato et al., 2014]. Также были использованы гравитационные модели спутников ГРЕЙС, которые позволяют исследовать косейсмические и постсейсмические процессы в районах крупных землетрясений (например, [Михайлов и др., 2014]).

Косейсмические смещения на пунктах сети GEONET рассчитаны группой ARIA (JPL и Caltech) по данным японского агентства GSI. Данные системы GEONET весьма точны, среднеквадратическая погрешность горизонтальной компоненты смещений оценивается в 3 мм, вертикальной компоненты в 15 мм, но станции расположены на островах, далеко от очаговой зоны землетрясения. В работе использованы данные 59 станций, смещения на которых превосходили 50 см. Пункты определения смещений на дне океана расположены непосредственно над очаговой зоной. Горизонтальные смещения во время землетрясения здесь достигли 24 м [Sato et al., 2014]. Точность определения смещений оценивается в 0.5 м. Повторные измерения проводились с интервалом в несколько месяцев, поэтому они содержат и косейсмические, и постсейсмические смещения [Perfettini, Avouac, 2014].

Поле смещений в направлении на спутник (LOS $–$ от английского line-of-sight) для обширного района землетрясения Тохоку-Оки было получено нами на основе 5 парных интерферограмм, которые были рассчитаны для различных фреймов снимков, выполненных радарами с синтезированной апертурой (РСА) спутника ENVISAT с 347 нисходящего трека 19 февраля и 21 марта 2011 года. Счет выполнен в программном пакете ENVI SARscape, использована цифровая модель рельефа SRTM3 v4. Также как и данные сети GEONET, эти измерения характеризуют области, удаленные от очаговой зоны. Интервал спутниковой радарной съемки составляет около одного месяца, поэтому смещения содержат и постсейсмическую компоненту.

Гравитационные модели спутников ГРЕЙС равномерно покрывают океанические и континентальные области вокруг очаговой зоны. Косейсмические изменения гравитационного поля рассчитаны по временным рядам гравитационных моделей как разница двух линейных трендов до и после землетрясения. Для расчетов использованы гравитационные модели Центра космических исследований (CSR), университета штата Техас, США, которые содержат 96 сферических гармоник, что соответствует длине полуволны около 200 км. Разрешение этих моделей относительно невелико, но эти данные важны потому, что гравитационные аномалии в основном характеризуют изменение плотности пород литосферы и верхней мантии в результате изменения их напряженного состояния.

Проблема построения модели, основанной на всех имеющихся данных, состоит в том, что данные РСА-интерферометрии фиксируют смещения многочисленных отражающих объектов в направлении на спутник, менее многочисленные станции GPS измеряют три компоненты смещений $–$ на восток, север и по вертикали. Ежемесячные гравитационные модели слабо отражают смещения земной поверхности, но, в основном, характеризуют косейсмические изменения поля напряжений, а, следовательно, и плотности пород в огромных объемах мантии и литосферы. Естественным выходом из этой ситуации является согласование всех групп данных в рамках геодинамической модели косейсмических и постсейсмических процессов [Михайлов и др., 2007].

Для описания косейсмических процессов нами использована модель смещений на прямоугольной площадке в сферически расслоенной модели Земли. Решение задачи получено в работе [Pollitz, 1996] и реализовано в программном коде Static1 D, доступном на сайте Геологической службы США (USGS). Распределение параметров по радиусу Земли (упругие модули, плотность, вязкость) были заданы по модели PREM. Поле смещений и напряжений в литосфере и мантии определяется размерами и положением площадки и величиной смещений на ней по падению и по простиранию. При этом задача является линейной: поля напряжений и смещения дневной поверхности, связанные со смещениями в каждом направлении на каждой площадке, суммируются. Вначале нами была построена модель поверхности сейсмического разрыва, состоящая из двух плоскостей, каждая из которых была разделена на девять равных элементов. Для каждого элемента разбиения, для единичного смещения по падению и по простиранию, рассчитаны: смещения в пунктах GPS на восток и север (52 пункта); смещения в направлении на спутник в тех точках земной поверхности, где они определены по данным радарной интерферометрии (996 значений); скачок гравитационного поля по сетке 1×1° для области 130.5$–$149.5° в.д. и 31.5$–$50.5° с.ш. (400 точек); смещения на восток и север в точках расположения донных акустических транспондеров (5 пунктов). В результате получено 1510 уравнений с 36 неизвестными смещениями по падению и по простиранию на каждом элементе разбиения поверхности сейсмического разрыва.

Обратная задача по определению поля смещений на поверхности разрыва по данным о смещениях земной поверхности и геофизическим полям является некорректно поставленной по Адамару. Для решения таких задач применяются различные методы регуляризации [Тихонов, Арсенин, 1986]. При аппроксимации поверхности разрыва конечной системой прямоугольных элементов происходит замена непрерывного поля смещений сеточной функцией, что уже является регуляризацией задачи. Также решение определялось нами при условии, что угол подвижки rake на каждом элементе разбиения близок к значению, определенному по механизму очага. Так как тангенс угла подвижки равен отношению смещений по падению $D{S}_i$ к смещению по простиранию $S{S}_i,$ это условие делает минимизируемый функционал нелинейной функцией смещений по простиранию и падению. Чтобы сохранить линейность мы добавили в функционал невязки по всем группам данных следующее условие:

$$\begin{array}{c}{‖D{S}_i\cos \left(rake\right)-S{S}_i\text{sin}\left(rake\right)‖}_{L_2}+\\ +{‖D{S}_i-\overline{DS}‖}_{L_2}+{‖S{S}_i-\overline{SS}‖}_{L_2}\text{}=\text{min}.\end{array}$$

(1)

Здесь первый член требует, чтобы смещения происходили вдоль заданного направления с азимутом rake. Но этому условию удовлетворяют два вектора: один с углом rake, а второй с углом rake + 180°. Чтобы исключить не имеющие геодинамического смысла смещения в обратном направлении, в (1) введены еще два условия: смещения на каждом элементе близки к среднему смещению по падению $\overline{DS}$ и по простиранию $\overline{SS}$ по всей поверхности разрыва. Как показала практика вычислений, такое условие эффективно устраняет смещения в направлении rake + 180°.

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

Сразу после землетрясения на сайте USGS были опубликованы две модели, построенные Г. Хайесом (см. [Hayes et al., 2011]). Поверхность разрыва была аппроксимирована прямоугольником, который в первой модели был разбит на 2 элемента, а в более позднем варианте на 325 элементов (описание дано в работе [Михайлов и др., 2014]). Мы использовали эту модель, но разделили ее на 18 элементов (рис. 1) и вначале построили решение только по данным GPS (рис. 1а) и только по данным РСА-интерферометрии (рис. 1б).

 

Рис. 1. Решение только по данным GPS (а) и только по данным РСА-интерферометрии (б). Изолинии на рисунках в левой части (а) и (б) – смещения в направлении на спутник по данным РСА-интерферометрии (см. Цветные символы – расчетные значения). Синие и красные стрелки – смещения по данным сети GEONET и расчетные смещения, соответственно. Черные стрелки – смещения на элементах поверхности разрыва. Серый контур – береговая линия, синяя линия – положение океанического желоба. Красная звезда – эпицентр главного события. Справа на рисунках (а) и (б) показана косейсмическая гравитационная аномалия (цветовая шкала) и расчетная аномалия (зеленые изолинии) в микрогалах.

 

Решение только по данным ГНСС дает меньшие смещения в направлении на спутник, чем смещения, определенные по данным РСА-интерферометрии. Смещения в пунктах GPS совпадают полностью. Расчетная гравитационная аномалия имеет несколько меньшую амплитуду, чем реальные аномалии. Расчетная отрицательная аномалия сдвинута на северо-запад. В решении только по данным спутниковой интерферометрии расчетные смещения на пунктах GPS заметно больше измеренных, амплитуда гравитационной аномалии близка к реальным данным. Модель, построенная только по данным GPS, дает меньшие смещения, чем оценки по РСА-интерферометрии потому, что в данные интерферометрии частично вошли и постсейсмические смещения. Область максимальных смещений расположена вблизи эпицентра землетрясения.

 

Рис. 2. Решение по всем группам данных без учета донных станций (а) и с донными станциями (б) для модели поверхности разрыва, состоящей из 18 элементов. Красные и синие стрелки над восточной частью поверхности разрыва – расчетные и измеренные смещения для пяти донных станций. Коричневые изолинии на рис. (б) слева – смещения на поверхности разрыва в метрах. Максимальное смещение на поверхности разрыва на рис. (а) составляет 24 м, на рис. (б) – 44 м. Обозначения см. рис. 1.

 

На рис. 2а приведено решение в той же геометрии, но полученное путем подбора под данные GPS, РСА-интерферометрии и гравитационные аномалии. Ошибки подбора невелики и могут быть уменьшены, если разбить плоскости на бóльшее количество элементов (36 или 72). Однако включение в подбор смещений по донным станциям (рис. 2б) не дает приемлемого приближения, особенно велика ошибка подбора смещений донных станций.

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

Для того, чтобы существенно улучшить качество подбора, необходимо продлить модель поверхности разрыва до океанического желоба. В этом случае (рис. 3) максимальные смещения получаются в верхней части контакта плит и превосходят 55 м. При этом качество подбора по всем группам данных оказывается весьма высоким.

 

Рис. 3. Решение, полученное по всем группам данных, для модели поверхности разрыва, которая продлена до океанического желоба и содержит 72 элемента: (а) – смещения на элементах модели поверхности разрыва в метрах и интерполяция поля смещений. Максимальное смещение превосходит 55 м. Синяя линия показывает положение океанического желоба; (б) – расчетные и реальные смещения на пунктах GPS и донных станциях (стрелки) и на отражателях спутникового радарного сигнала (цветные символы и изолинии); (в) – расчетные и реальные гравитационные аномалии в микрогалах. Обозначения см. рис. 1 и рис. 2.

 

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ И ВЫВОДЫ

Использованные в работе данные ГНСС, донных станций, РСА-интерферометрии и спутников ГРЕЙС, в целом хорошо согласуются между собой. Однако данные GPS, РСА-интерферометрии о смещениях на островах и гравитационные модели ГРЕЙС не имеют высокого разрешения. Получаемая на их основе модель поверхности разрыва оказывается сильно сглаженной. При этом измеренные и расчетные смещения согласуются очень хорошо, что создает иллюзию того, что полученное решение близко к реальному процессу.

Для оценки устойчивости и разрешающей способности методов построения поверхности сейсмического разрыва в западной литературе обычно используется так называемый шахматный тест. В соседних элементах модели поверхности разрыва задаются противоположные смещения одинаковой амплитуды (если плюс покрасить белым, а минус черным, то получится шахматная доска) и от этих смещений рассчитываются теоретические данные $–$ поля смещений и/или гравитационные аномалии. Далее решается обратная задача и если в расчетных смещениях также чередуются положительные и отрицательные значения, то считается, что на реальных данных модель также даст хороший результат.

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

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

Финансирование работы

Работа выполнена при финансовой поддержке РФФИ (проект $№$ 18-05-00159).

Благодарности

Временные ряды смещений ARIA v0.3 получены нами от группы ARIA, работающей в JPL и Caltech (США). Эти смещения были рассчитаны на основании первичных данных сети GEONET RINEX, принадлежащей Geospatial Information Authority (GSI) Японии. Данные ГРЕЙС получены с сайта CSR, университета штата Техас, США. Снимки спутника ENVISAT предоставлены европейским космическим агентством. Авторы выражают благодарность всем перечисленным организациям за предоставленные данные, а также рецензентам, которые внимательно прочли рукопись и сделали ряд важных замечаний.

About the authors

V. O. Mikhailov

Schmidt Institute of Physics of the Earth, Russian Academy of Sciences; Faculty of Physics, Moscow State University

Author for correspondence.
Email: mikh@ifz.ru
Russian Federation, Moscow

E. P. Timoshkina

Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Email: mikh@ifz.ru
Russian Federation, Moscow

S. A. Khairetdinov

Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Email: mikh@ifz.ru
Russian Federation, Moscow

P. N. Dmitriev

Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Email: mikh@ifz.ru
Russian Federation, Moscow

I. M. Kartashov

Schmidt Institute of Physics of the Earth, Russian Academy of Sciences; Faculty of Physics, Moscow State University

Email: mikh@ifz.ru
Russian Federation, Moscow

V. B. Smirnov

Faculty of Physics, Moscow State University; Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Email: mikh@ifz.ru
Russian Federation, Moscow

References

Supplementary files