Probabilistic estimates of hypocenters from the data of Kamchatka seismic network stations

Full Text

ВВЕДЕНИЕ

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

Данная работа направлена на разработку метода, позволяющего наиболее достоверно оценить качество определения координат гипоцентров и времени в источнике. Неопределенность в определении гипоцентра может существенно увеличиваться для землетрясений, зарегистрированных малым числом станций или в случае, когда большинство станций находится в стороне от эпицентра землетрясения, не обеспечивать приемлемый азимутальный охват эпицентра. Первая ситуация типична для слабых вулканических землетрясений и для тектонических землетрясений средней силы на краях зоны ответственности Камчатской сети станций [Чебров и др., 2018], а вторая $–$ для поверхностных землетрясений в зоне субдукции. Таким образом, проблема аккуратной оценки качества определения гипоцентров актуальна для большинства землетрясений, регистрируемых в сейсмически- и вулканически активных районах Дальнего Востока России.

Практически все методы локации гипоцентров основаны на минимизации невязок или разностей между расчетными и фактическими временами прихода сейсмических волн на станции регистрации. В качестве целевой функции часто используется сумма квадратов невязок. В связи с нелинейной зависимостью времен прихода от расстояния между эпицентром и станцией, эта задача оптимизации принципиально нелинейна, и, соответственно, возможны ситуации, когда нелинейность приводит к существенной неоднозначности решения, вплоть до существования более одного минимума целевой функции. Эффективные методы локации землетрясений появились в 60-х годах 20-го века благодаря применению компьютеров. Как обсуждается в работе [Husen, Hardebeck, 2010], на первых этапах в связи с ограниченной мощностью компьютеров задача оптимизации решалась итеративными методами. Хорошо известными примерами такого подхода являются программа HYPO71 [Lee, Lahr, 1975] или метод, долгое время используемый на Камчатке [Гусев, 1979].

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

Поэтому в последние годы все более популярными становятся статистические методы локации землетрясений [Sambridge, Kennet, 1986; Lomax et al., 2000; Minson, Lee, 2014]. Такие методы не ограничиваются поиском экстремума целевой функции, а оценивают ее распределение во всем пространстве возможных положений гипоцентра и времен в источнике. Таким образом, выявляется область достоверных решений. Для этого часто используется комбинация имеющихся моделей времен пробега волн и соответствующих ошибок в качестве априорной информации для расчета распределения апостериорной вероятности положения конкретных гипоцентров по конкретным данным о временах вступлений P- и S-волн на конкретных станциях с помощью теоремы Байеса. Таким образом, вместо одной точки наиболее вероятного значения гипоцентра будут получаться облака точек, для каждой из которых рассчитана вероятность того, что именно в ней находится реальный гипоцентр. В этой работе демонстрируется использование теоремы Байеса для локации землетрясений в Камчатской зоне субдукции и в районе активных вулканов Камчатки на основе данных Камчатского филиала ФИЦ «Единая геофизическая служба РАН». После описания теоретических основ метода и практической реализации алгоритма, мы обсуждаем, как неоптимальное распределение станций может приводить к неоднозначности решений разных типов землетрясений.

ИСХОДНЫЕ ДАННЫЕ

Информационной базой исследований являются результаты обработки землетрясений за 2010$–$2018 гг. по данным системы детальных сейсмологических наблюдений Камчатского филиала (КФ) ФИЦ ЕГС РАН [Чебров и др., 2013]. Для хранения и доступа к информации о землетрясениях в КФ созданы информационные ресурсы: http://www.emsd.ru/ts/all.php [Чебров и др., 2010] и Единая информационная система сейсмологических данных (ЕИС СД), http://www.emsd.ru/sdis [Чеброва и др., 2015]. В работе нами использованы времена вступлений P- и S-волн, а также каталог землетрясений за рассматриваемый период, составленный по расчетам программой DIMAS [Дрознин, Дрознина, 2010] и включающий определение гипоцентров более 150 000 тектонических и вулканических событий. Расчеты теоретических времен пробега волн проведены на основе регионального годографа P- и S-волн из пакета программ [Мельников, 1990], используемого в лаборатории исследования сейсмической и вулканической активности КФ ФИЦ ЕГС РАН для определения параметров тектонических землетрясений и, в связи с отсутствием уточненных локальных скоростных моделей, вулканических землетрясений. Исключение составляют землетрясения Северной и Авачинской групп вулканов, для которых используются локальные годографы [Сенюков, 2006].

МЕТОД РАСЧЕТОВ

Основной целью данной работы является представление расчетных гипоцентров тектонических и вулканических землетрясений Камчатки не в традиционном виде точечных оценок, а в виде трехмерного распределения вероятности положения гипоцентра. В качестве распределения вероятности мы используем апостериорное четырехмерное распределение вероятности положения гипоцентра и времени в очаге, рассчитываемое по формуле Байеса на основе априорных представлений о возможном решении, фактических данных о временах вступлений (в нашем случае волн P и S), физической модели скоростей сейсмических волн (годографа) и модели распределения невязок фактического и расчетного времен прохождения волн [Tarantola, Valette, 1982]. В общем виде такое решение представляется [Minson, Lee, 2014]:

 

$$p\left(\phi ,\lambda ,h,t\text{|}d\right)=F\cdot p\left(d\text{|}\phi ,\lambda ,h,t\right)p\left(\phi ,\lambda ,h,t\right),$$

(1)

где: d $–$ вектор данных о всех зарегистрированных временах вступлений сейсмических волн P и S на всех станциях: $p\left(\phi ,\lambda ,h,t\text{|}d\right)$ $–$ искомая апостериорная вероятность оценок положения гипоцентра и времени в очаге; $p\left(\phi ,\lambda ,h,t\right)$ $–$ априорное распределение возможного решения; $p\left(d\text{|}\phi ,\lambda ,h,t\right)$ $–$ вероятность реализации значений времен вступлений при данном положении гипоцентра и времени в очаге и при определенной модели ошибок, называемая правдоподобием; F $–$ нормировочный коэффициент, равный величине, обратной интегралу от $p\left(d\text{|}\phi ,\lambda ,h,t\right)p\left(\phi ,\lambda ,h,t\right)$ по всем возможным значениям $\left(\phi ,\lambda ,h,t\right).$ Отметим, что теорема Байеса позволяет пересчитать функцию правдоподобия данных для конкретного гипоцентра и времени в очаге в вероятность положения гипоцентра и времени в очаге при конкретных данных о временах вступления.

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

Функция правдоподобия $p\left(d\text{|}\phi ,\lambda ,h,t\right)$ обычно предполагается [Tarantola, Valette, 1982; Minson, Lee, 2014] имеющей вид многомерного нормального распределения:

 

$$p\left(d\text{|}\phi ,\lambda ,h,t\right)=\frac{1}{{\left(2\pi \right)}^{N_d}\sqrt{C_{\chi }}}{e}^{-\frac{1}{2}{\left(d-\widehat{d}\right)}^T{C}_{\chi }{}^{-1}\left(d-\widehat{d}\right)},$$

(2)

где $N_d$ $–$ размерность вектора данных d (число зарегистрированных вступлений); $\widehat{d}$ $–$ вектор расчетных времен вступления волн, зависящий от параметров $\left(\theta ,\phi ,h,t\right),$ и, соответственно, $d-\widehat{d}$ $–$ вектор невязок; $C_{\chi }$ $–$ матрица ковариации невязок; символ Т обозначает транспонирование вектора. Если предположить независимость невязок и общую величину ошибки, то матрица ковариации имеет диагональный вид со значениями, равными дисперсии невязок. Решение, соответствующее максимуму правдоподобия (2) в этом случае эквивалентно решению методом наименьших квадратов минимизацией среднеквадратичной невязки вступлений волн [Tarantola, Valette, 1982]. В разных вариантах вероятностного подхода для локации гипоцентров также обычно используется правдоподобие в форме (2). При этом по-разному решается проблема определения матрицы ковариации $C_{\chi }.$ В работе [Tarantola, Valette, 1982] постулируется корреляция времен вступления на разных станциях, зависящая от расстояния между станциями, но теоретическая ошибка определения времени пробега волны для разных станций предполагается равной и независящей от гипоцентрального расстояния. Эта ошибка принята равной 0.2 с. В работе [Minson, Lee, 2014], посвященной переопределению гипоцентров по историческим данным с использованием более совершенного годографа, предполагается, что ошибки расчета времени пробега волны пренебрежимо малы по сравнению с точностью определения времен вступления. В работе [Lomax et al., 2000] отмечается сложность реалистичного определения матрицы ковариации. Вместе с тем, очевидно, что размеры «эллипса ошибок» определения гипоцентра кардинальным образом зависят, как минимум, от величины ошибок (диагональные элементы матрицы ковариации $C_{\chi }).$ Поэтому для получения адекватных представлений об ошибке определения гипоцентра необходимы реалистичные определения матрицы ковариации $C_{\chi }.$

В нашей работе мы предполагаем зависимость диагональных элементов матрицы от расстояния между гипоцентром и соответствующей станцией. Мы предполагаем также, что недиагональные члены масштабируются значениями соответствующих диагональных членов, но допускается также зависимость от расстояний до гипоцентра и/или между станциями соответствующих коэффициентов корреляции. Построение матрицы ковариации несколько осложняется тем фактом, что для части станций могут отсутствовать времена вступлений либо S-, либо P-волн. Ниже приведен общий вид вектора невязок и матрицы ковариации $C_{\chi }$ с учетом сделанных предположений. Нижний индекс i соответствует какой-либо станции, верхний индекс $–$ волне P или S. В начале следуют $N^{PS}$ пар строк и столбцов матрицы и компонент вектора невязок, соответствующих станциям, для которых имеются вступления и P-, и S-волн, затем $N^P$ элементов, соответствующие станциям, для которых имеются вступления только волны P и, наконец, $N^S$ элементов, соответствующих станциям, для которых имеются вступления только волны S. Общее число вступлений волн составляет $2N^{PS}\text{}+N^P\text{}+N^S\text{},$ а число станций $N^{PS}\text{}+N^P\text{}+N^S.$ Вектор невязок составлен из разностей фактического и расчетного времен вступления волн. Расчетное время определено как $t+\delta t_P$ для волны P или $t+\delta t_S$ для волны S, где времена $\delta t_P$ и $\delta t_S$ пробега волн P и S от точки $\left(\phi , \lambda ,h\right)$ до станции рассчитываются по годографу.

 

$\begin{array}{c}\overrightarrow{ d}-\overrightarrow{\widehat{d}}=\left\{t_P^{1^{}}-\delta t_P^1,t_S^1-\delta t_S^1,t_P^2-\delta t_P^2,t_S^2-\delta t_S^2,\dots ,\right.\\ t_P^{N^{SP}}\text{}\text{}-\delta t_P^{N^{SP}}\text{},t_S^{N^{SP}}-\delta t_S^{N^{SP}}\text{},t_P^{N^{SP}+1}\text{}-\delta t_P^{N^{SP}+1},\dots ,\\ t_P^{N^{SP}+N^P}\text{}-\delta t_P^{N^{SP}+N^P}\text{},t_S^{N^{SP}+N^P+1}\text{}-\delta t_P^{N^{SP}+N^P+1},\dots ,\\ \left.t_S^{N^{SP}+N^P+N^S}\text{}-\delta t_P^{N^{SP}+N^P+N^S}\right\}\end{array}$ (3)

$$

(4)

 

Здесь мы предполагаем, что значения диагональных членов матрицы зависят только от расстояния ${\Delta }_i$ от точки $\left(\phi ,\lambda ,h\right)$ до i-ой станции. Для удобства расчетов здесь в качестве меры расстояния мы используем расчетное время $\delta t_i$ пробега P-волны; ${\rho }_{SP}$ $–$ коэффициент корреляции невязок разных типов волн. Здесь мы также пока предполагаем, что этот коэффициент может зависеть от расстояния. Ковариационные члены матрицы, соответствующие каждым двум разным станциям i, j, масштабируются квадратным корнем из соответствующих диагональных членов. Соответствующие коэффициенты $r_{PP},$ $r_{SS},$ $r_{PS},$ имеющие смысл коэффициентов корреляции, предполагаются зависимыми как от расстояния между этими станциями ${\delta }_{i,j},$ так и от среднего расстояния ${\Delta }_{i,j}=1/2({\Delta }_i+{\Delta }_j).$ Исследование зависимости параметров ${\sigma }_P, {\sigma }_S, {\rho }_{PS}$ от гипоцентрального расстояния и параметров $r_{PP}, r_{SS}, r_{PS}$ от среднего гипоцентрально расстояния и расстояния между станциями, а также построение модели невязок составляет предмет следующего раздела.

В том случае, если можно пренебречь членами матрицы ковариации, соответствующими разным станциям, и учитывать только корреляцию невязок P- и S-волн на одной и той же станции, матрица $C_{\chi }$ в уравнении (2) приобретает кусочно-диагональный вид, что позволяет получить в явном виде обратную матрицу, и тогда уравнение (1) может быть представлено в линейном виде:

 

$$\begin{array}{c}p\left(\phi ,\lambda ,h,t\right)=F_2{\displaystyle \prod }_{i=1}^{N^{PS}}\frac{1}{2\pi {\sigma }_p{\sigma }_s\sqrt{\left(1-{\rho }_{PS}{}^2\right)}}\times \\ \times e^{-\frac{1}{1-{\rho }_{PS}{}^2}\left\{\frac{\delta p_i{}^2}{2{\sigma }_p{}^2}-\rho \frac{\delta p_i\delta s_i}{{\sigma }_p{\sigma }_s}+\frac{\delta s_i{}^2}{2{\sigma }_s{}^2}\right\}}\times \\ \times {\displaystyle \prod }_{j=1}^{N^P}\frac{1}{\sqrt{2\pi }{\sigma }_p}{e}^{-\frac{\delta p_j{}^2}{2{\sigma }_p{}^2}}{\displaystyle \prod }_{k=1}^{N^S}\frac{1}{\sqrt{2\pi }{\sigma }_s}{e}^{-\frac{\delta s_k{}^2}{2{\sigma }_s{}^2}}.\end{array}$$

(5)

Здесь F2, аналогично (1), $–$ нормировочный коэффициент. Величины ${\sigma }_p,$ ${\sigma }_s$ и ${\rho }_{PS}$ могут зависеть от взаимного расположения данной станции и точки с координатами $\left(\phi ,\lambda ,h,t\right).$

В качестве априорного распределения при использовании теоремы Байеса наиболее часто используется так называемый «неинформативный приор» $–$ равномерное распределение параметра в заданных пределах. Мы знаем только, в каких границах должен находиться параметр, но не знаем, какие его значения более вероятны. В задаче локации гипоцентров обычно известно, в каких пределах должен находиться гипоцентр. В данной работе задача решается итеративно в два этапа, так как оперативной памяти компьютера хватает для ограниченного размера области при заданной разрешающей способности решения задачи. Иначе говоря, для больших областей возможно нахождение решений лишь на редкой пространственно-временной сетке. Мы используем идею подхода, получившего в англоязычной литературе обозначение “importance sampling method” [Lomax et al., 2000].

Для ускорения процесса счета на каждом шаге итерации используется пространственно-временная сетка размерностью 41 × 41 × 41 × 41. На первом шаге первого этапа границы равномерного априорного распределения по пространству задаются в виде параллелепипеда, в горизонтальной проекции описывающего границы зоны ответственности Камчатского филиала ФИЦ ЕГС РАН [Левина, 2008; Чебров, 2018] с некоторым запасом, а по глубине от $–$5 до 700 км. По времени априорные пределы интервала значений задаются как минимально и максимально возможные для заданной точки с учетом теоретического времени пробега. На каждом последующем шаге первого этапа происходит уменьшение района в 2 раза с центром в точке максимума $p\left(\phi ,\lambda ,h,t\right)$ до достижения заданной точности (величины шага). На втором этапе происходит постепенное увеличение района путем смещения сетки с сохранением шага по всем параметрам и ведется подсчет апостериорной вероятности $p\left(\phi ,\lambda ,h,t\right).$ Процесс останавливается, когда вероятность положения гипоцентра в границах текущей области сетки превышает 99%. Такая процедура обеспечивает, чтобы внутри области интеграл от ненормированной апостериорной вероятности как минимум в 100 раз превышал аналогичный интеграл вне области, если бы расчет можно было производить с той же разрешающей способностью в изначальных границах.

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

МОДЕЛЬ РАСПРЕДЕЛЕНИЯ НЕВЯЗОК P- И S- ВОЛН

Ошибки невязок и матрица ковариации $C_{\chi }$ в целом имеют две аддитивные компоненты: «теоретическую» $C_T$ и «фактическую» $C_t,$ $C_{\chi }=C_T+C_t$ [Tarantola, 1982; Гусев, 1979]. Теоретическая компонента обусловлена неточностью годографа и неоднородностью среды. Фактическая компонента связана с неточностью данных. В современных условиях, когда координаты сейсмических станций известны с высокой точностью, время на сейсмограммах калибруется по GPS, фактическая компонента может быть обусловлена только неточным определением времен вступления. На практике разделить компоненты очень сложно, поэтому в данной работе мы пошли по пути исследования эмпирических закономерностей, используя большой фактический материал. Из Камчатского каталога выбраны землетрясения с энергетическим классом ${\text{K}}^{ф}$ ≥ 7 (класс по S-волне, определенный по номограмме С.А. Федотова [Федотов, 1972]), для которых в станционной базе данных имеется не менее 5 вступлений P и не менее 5 вступлений S. Составлена база данных «оптимальных» невязок $–$ разности расчетного (с учетом гипоцентра и времени в очаге по данным каталога) и фактического времен вступления сейсмических волн ${\delta }^p\text{},{\delta }^s.$ База данных включила более 106 значений невязок.

Задача исследования невязок в данном разделе $–$ определить зависимость параметров ${\sigma }_P, {\sigma }_S, {\rho }_{PS}$ от гипоцентрального расстояния (расстояния от предполагаемого гипоцентра с координатами $\left(\phi ,\lambda ,h\right)$ до соответствующей станции) и параметров $r_{PP}, r_{SS}, r_{PS}$ от среднего гипоцентрального расстояния и от расстояния между соответствующими станциями.

Диагональные элементы матрицы ковариации $C_{\chi }$ могут зависеть от взаимного расположения очага землетрясения и станции [Гусев, 1979]. Удобной мерой такого взаимного расположения в процедуре локации гипоцентра является расчетное время $\delta t_p$ прохождения волны P, поскольку эта величина непосредственно используется в расчетах. Как оказалось (рис. 1), разброс невязок возрастает с увеличением $\delta t_p.$

Невязки P- и S-волн на одной и той же станции нельзя считать независимыми. При расчетах гипоцентра этот эффект может оказаться весьма существенным, особенно в тех случаях, когда на части станций зарегистрировано лишь одно из времен вступления волн P или S. В любом случае, если корреляция невязок не учитывается, то искусственно завышается точность определения гипоцентра. В принципе, коррелированность невязок может меняться от станции к станции и в зависимости от местоположения очага землетрясения, при этом вклад могут давать как «теоретическая», так и «фактическая» составляющие ковариации. Здесь мы предполагаем, что эти детали несущественны по сравнению с усредненной моделью корреляции невязок P- и S-волн. На рис. 2 показана зависимость невязок для волны S в зависимости от невязки волны P для той же станции и того же землетрясения, построенная по всей базе данных невязок, описанной выше. Коэффициент корреляции ${\rho }_{PS}$ составляет 0.55.

 

Рис. 1. Зависимость разности расчетного и фактического времен вступления волн P (а) и S (б) от расчетного времени волны P по всем станциям и всем землетрясениям с классом 2010–2018 гг.

 

Чтобы исследовать зависимость от величины $\delta t_p$ значений ошибок ${\sigma }_p,$ ${\sigma }_s$ и коэффициента корреляции ${\rho }_{PS},$ мы упорядочили базу данных по возрастанию $t_g^p$ и затем в скользящих окнах, включающих 1000 пар невязок P и S на одной и той же станции и от одного и того же землетрясения, вычисляли среднеквадратичные отклонения и выборочный коэффициент корреляции, а также среднее значение $\delta t_p$ для группы данных. Соответствующие зависимости представлены в двойном логарифмическом масштабе на рис. 3.

 

Рис. 2. Значения разности расчетного и фактического времен вступления волны S от соответствующей разности волны P по всем станциям и всем землетрясениям с классом 2010–2018 гг.

 

Рис. 3. Зависимость ошибок (а), (б) и коэффициента корреляции невязок (в) от расчетного времени прохождения волны P от гипоцентра до станции. Прямые линии отмечают линейную аппроксимацию зависимостей: ${\mathit{\sigma }}_{\mathbf{P}}\left(\mathbf{\delta }{\mathbf{t}}_{\mathbf{p}}\right)\mathbf{=}\mathbf{0}\mathbf{.14}\mathbf{ }\mathit{\delta }{\mathit{t}}_{\mathbf{P}}^{\mathbf{0}\mathbf{.42}}$ и ${\mathit{\sigma }}_{\mathbf{S}}\left(\mathbf{\delta }{\mathbf{t}}_{\mathbf{p}}\right)\mathbf{=}\mathbf{0}\mathbf{.16}\mathbf{ }\mathit{\delta }{\mathit{t}}_{\mathbf{P}}^{\mathbf{0}\mathbf{.53}}$.

 

Как видно из рисунка, зависимости ошибок ${\sigma }_p$ и ${\sigma }_s$ от $\delta t_p$ хорошо аппроксимируются в двойном логарифмическом масштабе прямыми линиями при $\delta t_p\ge 10 c\text{.}$ Отклонение от этих зависимостей при малых значениях ${\sigma }_p$ и ${\sigma }_s$ в сторону больших значений можно интерпретировать как относительное возрастание вклада «фактической» компоненты в матрицу ковариации: по-видимому, ошибка определения времен вступления волн P и S составляет, соответственно, около 0.3 и 0.5 c. Коэффициент корреляции ${\rho }_{PS}$ невязок P и S, как оказалось (рис. 3в), слабо зависит от величины $t_g^p.$ Поэтому в качестве окончательной модели зависимости ошибок ${\sigma }_p$ и ${\sigma }_s$ от величины $\delta t_p$ и модели величины ${\rho }_{PS}$ мы приняли соотношения:

${\sigma }_p=\max \left\{\mathrm{0.3,} 0.14 \delta t_p^{0.42}\right\}$

$${\sigma }_s=\text{max}\left\{\mathrm{0.5,}0.16 \delta t_p^{0.53}\right\}$$

(6)

${\rho }_{PS}=\mathrm{0.55.}$

При анализе коэффициентов, соответствующих парам различных станций следует учитывать две разные ситуации. Корреляция невязок на близких станциях должна быть малой, если расстояние от обеих станций до гипоцентра невелико. При больших гипоцентральных расстояниях, наоборот, корреляция очень близких станций должна быть близкой к единице. При этом, однако, «фактическая» составляющая ошибок может снижать эту корреляцию из-за относительно существенной ошибки определения времен вступления. В данной работе мы моделируем первую ситуацию условием $(\delta \ge \Delta ),$ где $\delta$ $–$ расстояние между парой станций, Δ $–$ среднее расстояние от гипоцентра до каждой из двух станций. В этом случае мы принимаем значения коэффициентов $r_{pp}, r_{ss}$ и $r_{ps}$ равными 0. Для второй ситуации $(\delta <\Delta ),$ в скользящем окне по величине $\delta$ шириной 0.1 град, мы оцениваем выборочную ковариацию для невязок, нормированных на соответствующие ошибки, определяемые соотношением (6). Благодаря нормировке коэффициенты $r_{pp,} r_{ss}$ и $r_{ps}$ безразмерны и имеют смысл коэффициентов корреляции. При этом нормировка учитывает увеличение ошибок и, соответственно, всех ковариационных членов с увеличением гипоцентрального расстояния. Результаты представлены на рис. 4.

Окончательно в качестве модели мы приняли следующие соотношения:

 

$$\begin{array}{c}{r}_{pp}\left(\delta |\right)=r_{ss}\left(\delta \right)=\left\{\begin{array}{c}0.55 e^{-\delta /0.15}, \delta <\Delta \\ \mathrm{0,} \delta \ge \Delta \end{array}\right.,\\ r_{ps}\left(\delta \right)=\left\{\begin{array}{c}0.3 e^{-\delta /0.15}, \delta <\Delta \\ \mathrm{0,} \delta \ge \Delta \end{array}\right..\end{array}$$

(7)

 

Рис. 4. Зависимость оценок коэффициентов и от расстояния между станциями (град) при условии Кругами отмечены оценки параметра квадратами – и ромбами Сплошная линия отмечает экспоненциальную аппроксимацию зависимостей: штриховая линия – зависимости

ПРИМЕРЫ АПОСТЕРИОРНЫХ РАСПРЕДЕЛЕНИЙ ПОЛОЖЕНИЯ ГИПОЦЕНТРОВ

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

На рис. 5 показано решение для Охотоморского землетрясения 24.05.2013, M_W = 8.3 [Сильные…, 2014]. Большое число станций, зарегистрировавших землетрясение и хорошее окружение станциями эпицентра землетрясения, позволило получить весьма точную оценку, несмотря на большую глубину очага землетрясения. Линейный размер 95% доверительной области составляет около 10 км в горизонтальной плоскости и около 30 км по глубине. Аналогичное по точности решение получено и для неглубокого Ближне-Алеутского землетрясения 17.07.2017, M_W = 7.8 (рис. 6). Линейный размер 95% доверительной области для эпицентра составляет около 5 км и для глубины очага около 30 км. Следует отметить, что данные оценки характеризуют предельную возможную точность для определения гипоцентров землетрясений в этих районах при существующей сети сейсмических станций.

 

Рис. 5. Апостериорные распределения вероятностных оценок положения гипоцентра и времени в очаге Охотоморского землетрясения 24.05.2013, M_W = 8.3. На карте слева треугольниками показаны станции, данные которых использованы для решения, область решения отмечена перекрестием линий и линиями уровня (область решения очень мала в масштабе карты). На диаграмме справа показаны линии уровня апостериорной вероятности для трех проекций, а также линии уровня маргинального распределения для времени в очаге. Линии уровня от 0 до 95% показаны оттенками серого цвета.

 

Рис. 6. Апостериорные распределения вероятностных оценок положения гипоцентра и времени в очаге Ближне-Алеутского землетрясения 17.07.2017, M_W = 7.8. Обозначения как на рис. 5.

 

При небольшом числе станций и плохом окружении эпицентра станциями любая точечная оценка гипоцентра может оказаться неоднозначной из-за протяженности области оптимального решения. На рис. 7 показана доверительная область решения, определенная по разработанной методике для одного из афтершоков Ильпырского землетрясения 13.03.2013, M_W = 5.8 [Сильные…, 2014]. Линейная форма области максимума решения объясняет обнаруженный ранее артефакт линейного расположения облака афтершоков этого землетрясения [Салтыков, Дрознина, 2014]. Такая форма доверительной области решения связана с неудачным расположением станций относительно эпицентра. Точечные решения для разных афтершоков оказываются случайным образом разбросанными вдоль этой области.

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

 

Рис. 7. Апостериорные распределения вероятностных оценок положения гипоцентра и времени в очаге афтершока Ильпырского землетрясения 13.03.2013, M_W = 5.8. Данный афтершок магнитуды ML = 4.0 произошел через 1 ч 19 мин после основного толчка. Обозначения как на рис. 5.

 

Рис. 8. Апостериорные распределения вероятностных оценок положения гипоцентра и времени в очаге сейсмического события в районе вулкана Удина 24.01.2018, ML = 1.5. Обозначения как на рис. 5.

 

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

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

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

Полученные доверительные области расположения гипоцентра для двух сильнейших землетрясений вблизи Камчатки за последние 10 лет (Охотоморского, M_W = 8.3 2013 г. и Ближне-Алеутского, M_W = 7.8 2017 г.) демонстрируют наилучшую возможную точность решений при существующей конфигурации сейсмической сети. Для более точного определения гипоцентра, очевидно, необходимо развитие сети станций, включая установку донных станций. Развитие сети требуется также для исследования вулканической сейсмичности.

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

В настоящей работе использовался региональный годограф P- и S-волн [Мельников, 1990], соответствующей одномерной скоростной модели среды. Такой же подход используется в настоящее время при рутинной обработке в Камчатском филиале. В «одномерном» годографе времена пробега волн зависят только от глубины источника и эпицентрального расстояния, но не от географического положения источника и станций. Использование одномерного годографа приводит к систематическим ошибкам в теоретических временах вступлений волн. Такие систематические ошибки могут быть связаны с неоднородностями внутри Земли. В этом случае они несут полезную информацию, которая может использоваться в томографии. Важно отметить, что в рамках предложенного здесь метода определения гипоцентров одномерный характер годографа не имеет принципиального значения. Вместо него можно использовать теоретические времена вступлений, вычисленные на основе трехмерных скоростных моделей. При таком подходе, расчетные времена δt_p и δt_s в уравнении (3) должны быть взяты не из годографа, а из таблицы, рассчитанной по трехмерной модели на сетке для всех возможных положений источников на всех станциях. При современных мощности и объемах памяти компьютеров такой подход в принципе возможен. Хотя эта проблема выходит за рамки данной статьи, хотелось бы отметить, что в ближайшем будущем может иметь смысл опробовать для улучшения точности гипоцентров современные трехмерные томографические модели, полученные на основе объемных сейсмических волн для Камчатского региона в целом [Gorbatov et al., 1999; Gontovaya et al., 2007; Koulakov et al., 2011a] и их более детальные версии для вулканических районов [Gontavaya et al., 2004; Koulakov et al., 2011b; 2017]. Сейсмические скорости в наиболее близких к поверхности слоях могут быть определены при помощи поверхностно-волновой томографии [Гордеев и др., 2009; Koulakov et al., 2014; Яновская, 2015; Droznina et al., 2017].

ВЫВОДЫ

В данной работе для обработки данных землетрясений Камчатки и прилегающих районов, зарегистрированных сетью станций детальных сейсмологических наблюдений КФ ФИЦ ЕГС РАН, впервые используется статистический метод определения положения гипоцентра и времени в очаге, в котором решением является распределение вероятности этих величин.

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

Для фактических расчетов предложена двухуровневая итерационная схема дискретизации пространства-времени землетрясений, существенно ускоряющая расчет, что в будущем позволит использовать метод при рутинной обработке в КФ ФИЦ ЕГС РАН. Метод применим и в других районах при условии построения адекватной для региона модели невязок.

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

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

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

Исследование выполнено при поддержке гранта Министерства образования и науки $№$14.W03.31.0033 в рамках Госзадания ФИЦ ЕГС РАН по теме НИР $№$ АААА-А16-116070550057-7 «Проведение непрерывного сейсмологического, геофизического и геодинамического мониторинга на глобальном, федеральном и региональном уровнях: разработка и внедрение новых технологий обработки и системного анализа больших объемов данных».

About the authors

D. V. Droznin

Kamchatka Branch, Geophysical Survey, Russian Academy of Sciences

Author for correspondence.
Email: shebalin@mitp.ru
Russian Federation, Petropavlovsk-Kamchatskii

S. Ya. Droznina

Kamchatka Branch, Geophysical Survey, Russian Academy of Sciences

Email: shebalin@mitp.ru
Russian Federation, Petropavlovsk-Kamchatskii

S. L. Senyukov

Kamchatka Branch, Geophysical Survey, Russian Academy of Sciences

Email: shebalin@mitp.ru
Russian Federation, Petropavlovsk-Kamchatskii

D. V. Chebrov

Kamchatka Branch, Geophysical Survey, Russian Academy of Sciences

Email: shebalin@mitp.ru
Russian Federation, Petropavlovsk-Kamchatskii

N. M. Shapiro

Institut de Physique du Globe de Paris; Schmidt Institute of Physics of the Earth, Russian Academy of Sciences

Email: shebalin@mitp.ru
Russian Federation, Paris; Moscow

P. N. Shebalin

Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences

Email: shebalin@mitp.ru
Russian Federation, Moscow

References

Supplementary files