Вулканология и сейсмология, 2020, № 1, стр. 330-48

ВВЕДЕНИЕ

Прогноз таких природных катастроф, как извержения вулканов и разрушительные землетрясения, является фундаментальной научной проблемой, с древнейших времен привлекающей внимание исследователей. Однако несмотря на пристальное внимание, эта проблема по-прежнему далека от разрешения. С современным состоянием проблемы можно ознакомиться в работе [Encyclopedia …, 2016].

Данная работа представляет собой развитие методов саморазвивающихся процессов и картирования сейсмической активности по плотности сейсмического потока [Тихонов, 2006, 2009]. Автор этих строк, являясь первоначальным разработчиком обоих вышеупомянутых методов, ни в коем случае не претендует на окончательное решение при их помощи проблемы прогноза природных катастроф, преследуя более скромные цели – оценить прогнозные возможности уравнения динамики саморазвивающихся природных процессов (уравнения ДСПП, см. ниже) и, в случае положительных перспектив, довести вышеупомянутые методы до стадии пригодности к практическому использованию. Подобная формулировка задачи предполагает постановку двух групп вопросов. Первая (оценочная) группа определяет правомерность самой постановки задачи прогноза при помощи указанных методов (насколько адекватно использование метода саморазвивающихся процессов при прогнозе потока сейсмической энергии? как соотносятся с прогнозируемостью сейсмического потока входящие в его состав сильные землетрясения?). Вторая группа вопросов соответствует изучению возможности практического использования прогнозируемости сейсмического потока (возможно ли заблаговременное выделение зон активизации, предшествующей сильным землетрясениям? возможен ли количественный прогноз усиления сейсмической активности в этих зонах? и др.). Естественно, что постановка практических вопросов возможна лишь при условии правомерности подобного прогностического подхода, т.е. после получения ответов на вопросы первой группы. Поэтому первоначально предлагается абстрагироваться от прогноза сильных землетрясений и сконцентрировать внимание на прогнозируемости потока энергии землетрясений E, затем рассмотреть взаимосвязь прогнозируемости этого параметра с сильными землетрясениями, и, при наличии взаимосвязи, обсудить перспективы ее использования в прогнозах сильных землетрясений и афтершокового затухания сейсмичности.

Для объективной оценки как прогнозных возможностей, так и перспектив их практического использования подобный анализ предполагается выполнить по нескольким регионам на основе (по возможности) различных источников сейсмических данных. В частности, подобное исследование выполнено для сейсмичности Камчатки по данным регионального каталога [Малышев, 2019], для потока сейсмической энергии северо-западного обрамления Тихого океана [Малышев, Малышева, 2018] и Северной Америки по данным каталога Геологической службы США, а также для сейсмичности Японии по данным каталога Японского метеорологического агентства.

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

МЕТОДЫ ИССЛЕДОВАНИЯ

В качестве математической модели в работе используется уравнение ДСПП [Малышев, 1991, 2005]

где x – любая неубывающая количественная характеристика процесса, x' и x'' – ее производные по времени, k – коэффициент пропорциональности, а показатели степени λ и α определяют нелинейность процесса соответственно в окрестностях стационарного состояния (x' ≈ $x_{0}^{'}$) и на значительном от него удалении (x' $ \gg $ $x_{0}^{'}$). Случай x' $ \gg $ $x_{0}^{'}$ интересен для прогноза потенциально катастрофических процессов. Поэтому в качестве аппроксимационной модели имеет значение уравнение:

Логический смысл уравнения (1) сводится к сделанному в терминах теории подобия предположению [Малышев, 2000], что в случае саморазвивающихся процессов “силы”, возникающие при отклонении системы от стационарного состояния, формируются за счет развития самой системы и пропорциональны разности “энергии движения” системы в текущем и стационарном состояниях: F_x = a|E_x – E₀|^γ или m_xx'' = a|m_x(x')²/2 – ${{{{m}_{x}}\left( {x{\kern 1pt} '} \right)_{0}^{{\text{2}}}} \mathord{\left/ {\vphantom {{{{m}_{x}}\left( {x{\kern 1pt} '} \right)_{0}^{{\text{2}}}} {\text{2}}}} \right. \kern-0em} {\text{2}}}$|^γ. Здесь m_x, F_x и E_x – соответственно “мера инертности”, “сила” и “энергия движения” системы по параметру x. Несложные преобразования последнего выражения приводят к нелинейному дифференциальному уравнению второго порядка (1) при λ = 2, α = 2γ и k = $am_{x}^{{{\gamma }--{\text{1}}}}.$ В связи с вышесказанным на основе уравнения (1) представляется возможным предложить универсальную методику прогноза количественных характеристик саморазвивающихся природных процессов.

Уравнение (2) соответствует уравнению, предложенному Б. Войтом [Voight, 1988] для описания динамики нарастания вулканической активности в преддверии кульминации извержения. Уравнение Б. Войта используется в методе прогноза разрушений FFM (Forecasting Failure Method), однако в ряде современных работ [Bell et al., 2011, 2013, 2016] утверждается, что данный метод необъективен и неточен даже для ретроспективного анализа. Это утверждение свидетельствует, что исследователи, пытаясь использовать уравнение Б. Войта, не смогли добиться устойчивого моделирования даже на стадии аппроксимации. Причина неудачи в данном случае – проблема с выбором критерия оптимизации [Малышев, 2016]: метод наименьших квадратов, применяемый исследователями, неэффективен.

В отличие от уравнения (1), решения которого в общем случае не представимы в явном виде и требуют численного интегрирования, уравнение (2) легко решается аналитически, причем его решения представляют собой либо собственно линейную зависимость (k = 0), либо сводятся к линейным зависимостям при логарифмировании.

Поиск наилучшей аппроксимации осуществляется посредством оптимизации по максимуму коэффициента упорядоченности K_reg = [n(x_n – x₁)(t_n – – t₁)/Σ(Δx_iΔt_i)]^0.5. Здесь n – число точек на аппроксимируемом участке фактических данных; (x_n – x₁) и (t_n – t₁) – диапазоны изменения фактических данных на аппроксимируемом участке соответственно по параметру x и времени t (выполняют функции нормирования обеих координат на диапазон изменений от 0 до 1), Δx_i и Δt_i – отклонения каждой точки фактических данных от расчетной кривой соответственно по оси абсцисс и по оси ординат. На практике удобен десятичный логарифм коэффициента упорядоченности – уровень упорядоченности L_reg = lg(K_reg).

Под прогнозируемостью здесь и далее понимается нахождение фактических данных “будущего” в полосе допустимых ошибок относительно расчетной кривой в ее экстраполяционной части. Для оценки прогнозируемости используется среднее отклонение σ фактических точек от расчетной кривой по нормали в координатах, нормированных на диапазон от 0 до 1. Затем аппроксимация экстраполируется в “будущее” до тех пор, пока нормальное расстояние каждой последующей (прогнозируемой) фактической точки до расчетной кривой находится в полосе допустимых ошибок (±3σ).

Количественная оценка дальности прогноза определяется через величину прогнозной дистанции Δ = {[(t_p – t_n)/(t_n – t₁)]² + [(x_p – x_n)/(x_n – x₁)]²}^0.5, где x_p и t_p – значения параметра и времени предельного прогнозируемого события, x_n и t_n – соответствующие значения для “текущего” события и x₁ и t₁ – для начального события в опорной (для аппроксимации и последующего прогноза) последовательности. Проекции прогнозной дистанции Δ на оси координат характеризуют дальность прогноза (прогнозные дистанции) по времени Δ_t и параметру Δ_x. Для оценки качества прогноза используется относительная прогнозная дистанция Δ_rel = Δ/σ или ее десятичный логарифм – уровень прогнозируемости L_p = lg(Δ_rel).

Для дифференцированной оценки прогнозной статистики по активизации и затуханию, а также для определения важных для прогноза значений показателя степени нелинейности α в уравнениях (1) и (2) используется коэффициент прогнозной нелинейности K_pn = Δ_rel × ${\text{lg}}\left| {{{x_{p}^{'}} \mathord{\left/ {\vphantom {{x_{p}^{'}} {x_{n}^{'}}}} \right. \kern-0em} {x_{n}^{'}}}} \right|,$ где $x_{p}^{'}$ = x'((t_n + t_p)/2) – прогнозируемая на середину прогнозного интервала скорость изменения параметра, а $x_{n}^{'}$ = x'(t_n) – ее текущее значение. Более подробно методика оценки прогнозируемости изложена в работе [Малышев, 2016].

При пространственном анализе сейсмических данных оценка прогнозируемости осуществляется по фиксированным сферическим гипоцентральным выборкам с радиусами 1.5, 3, 7.5, 15, 30, 60 и 150 км. Выборки распределены по широте, долготе и глубине с шагом смещения, в 1.5 раза меньшим радиуса выборки (т.е. соответственно 1, 2, 5, 10, 20, 40 и 100 км), что обеспечивает пространственное перекрытие выборок и исключает пропуск данных для прогностических оценок.

Ради определенности под сильными землетрясениями понимаются такие землетрясения, которые в кумулятивном распределении по энергии превышают порог 99.9% от общего числа землетрясений. В связи с различиями в регистрируемости сейсмического потока в разных районах земного шара данное определение применяется по фиксированным сферическим гипоцентральным выборкам с радиусом 100 км, в которых зарегистрировано более 1000 землетрясений. В том случае, если в выборке зарегистрировано от 100 до 1000 землетрясений, под сильным землетрясением понимается самое сильное землетрясение выборки. В выборках со слабой (от 10 до 100) и очень слабой (менее 10 землетрясений) регистрируемости сейсмического потока землетрясений сильные землетрясения в составе выборок не выделяются.

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

В качестве исходных данных в работе используется Всемирный каталог Геологической службы США по состоянию на 1 января 2017 г.¹¹ Прогнозируемость потока сейсмической энергии Южной Европы и Средиземноморья изучалась в пределах координат 30°–55° по широте при долготе –15°…+45°. В рассматриваемых пределах на глубинах от –2.0 до 634.0 км расположены гипоцентры 132710 землетрясений с магнитудой M = = 0–7.8 при ее среднем значении 2.9 и модальном – 2.7. В распределении толчков по магнитуде при M = 0 имеется пик из 1056 событий на фоне их полного отсутствия для M < 0 и единичных событий с магнитудой M = 0.1–0.3 (9 событий). Представляется вероятным, что эти толчки представлены землетрясениями, для которых магнитуда не установлена и ее нулевое значение является артефактом. Тем не менее, реальные толчки с M = 0 обычны для каталога Геологической службы США и хорошо соответствуют распределению землетрясений по магнитуде для других регионов, например, Северной Америки. Поэтому какая бы то ни было фильтрация подобных толчков в используемом алгоритме не предусмотрена, и они обрабатываются как есть. Предполагается, что полученные оценки прогнозируемости имеют дополнительный потенциал к улучшению за счет уточнения данных для толчков с M = 0.

В соответствии с приведенным выше определением в числе землетрясений Южной Европы и Средиземноморья выделяются 217 сильных толчков с M = 3.4–7.8. В качестве параметра x рассматривается сумма энергии землетрясений E. При этом энергия одиночного землетрясения в Дж оценивается согласно имеющейся зависимости между его магнитудой M и энергетическим классом K [Kanamori, 1977]: K = 1.5M + 4.8.

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

Сводная статистика ретропрогнозных оценок приведена в табл. 1. В большинстве прогнозных определений (72.2%) прогнозная дистанция Δ более чем в 3 раза превышает величину средней ошибки σ, т.е. подобные прогнозные определения рассматриваются как значимые. Максимальные уровни прогнозируемости L_p соответствуют превышению прогнозной дистанции Δ над средним отклонением σ более чем на 6 десятичных порядков, причем максимальная прогнозируемость потока сейсмической энергии отмечается на средних гипоцентральных радиусах выборок (7.5–15 км). Средние уровни прогнозируемости соответствуют превышению прогнозной дистанции Δ над средним отклонением σ на 1.2 порядка. Средневзвешенные значения прогнозных дистанций по времени Δ_t составляют десятки и сотни дней для активизации сейсмичности и несколько лет для ее затухания.

На рис. 1 отображены значения прогнозной нелинейности²² для всех определений, у которых |K_pn| > 1, а также отображены экстремумы прогнозной нелинейности, связанные с сильными землетрясениями. Как можно видеть, экстремумы прогнозной нелинейности, свойственные сильным землетрясениям, как правило, соответствуют аналогичным экстремумам K_pn для сейсмического потока в целом.

Статистические данные прогнозно значимого показателя степени α демонстрируют (см. табл. 1, рис. 2), что нелинейность как активизации сейсмического процесса, так и его затухания определяется классом гиперболических функций (1 < α < 2). При этом форма зависимостей близка к равнобокой гиперболе (α = 1.5) для зависимостей активизации и смещена к значениям α = 1.8 для затухания. Тем не менее достаточно большой разброс значений показателя α не позволяет ограничиться равнобокой гиперболой в качестве математической модели для ретропрогнозных оценок даже в случае зависимостей активизации.

В табл. 2 приведена общая статистика ретропрогнозов сильных землетрясений Южной Европы и Средиземноморья по потоку сейсмической энергии. Форшоковая прогнозируемость фиксируется для 84 из 217 сильных землетрясений региона. Она характеризуется ~1.6 тысячами ретропрогнозных определений, для которых сильное землетрясение оказывается в полосе ошибок (±3σ) экстраполяционной части прогнозной зависимости (рис. 3). Форшоковая прогнозируемость сильных землетрясений по потоку сейсмической энергии в количестве прогнозируемых землетрясений начинает проявляться на малых (1.5–7.5 км) радиусах гипоцентральных выборок, быстро возрастает на средних (15 и 30 км) радиусах, затем более плавно увеличивается на радиусе 60 км и несколько снижается при радиусе 150 км. Более высокие уровни прогнозируемости и прогнозной нелинейности отмечаются на средних гипоцентральных радиусах (7.5 и 15 км). Средняя прогнозная дистанция по времени в форшоковой прогнозируемости составляет десятки и сотни дней, что соответствует региональным показателям для активизации сейсмичности в целом.

Статистика ретропрогнозов затухания сейсмичности после сильных землетрясений приведена в табл. 3. Афтершоковая прогнозируемость по выделяющейся сейсмической энергии прослеживается для 200 из 217 сильных землетрясений и характеризуется ~72 тысячами ретропрогнозных определений, что почти на два порядка больше по сравнению с форшоковой прогнозируемостью. Афтершоковая прогнозируемость затухания сейсмичности после сильных толчков появляется на малых радиусах (1.5 и 3 км), быстро нарастает на средних (от 7.5 до 30 км) радиусах, достигает максимума при радиусе 60 км, а затем несколько снижается при 150 км. Наиболее высокие максимальные уровни прогнозируемости и прогнозной нелинейности отмечаются для выборок с гипоцентральным радиусом 7.5 км, тогда как наиболее высокие средние уровни – на радиусах 3 км. По мере увеличения радиусов прослеживается тенденция к уменьшению средних уровней прогнозируемости и прогнозной нелинейности. Средняя прогнозная дистанция по времени в афтершоковой прогнозируемости на два порядка выше аналогичного параметра форшоковой прогнозируемости: она увеличивается с полугода на радиусе 1.5 км до 6–8 лет на радиусах 7.5–60 км, а затем уменьшается до 3.5 лет на радиусе 150 км.

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

Таблицы 4 и 5 содержат сведения по статистике соответственно форшоковой и афтершоковой прогнозируемости некоторых сильных землетрясений Южной Европы и Средиземноморья по потоку сейсмической энергии. Таблицы 6 и 7 содержат характеристики ретропрогнозных зависимостей с экстремальной прогнозной нелинейностью соответствующих землетрясений из табл. 4 и 5. Рисунок 4 иллюстрирует пространственное распределение сильных землетрясений Южной Европы и Средиземноморья и их форшоковую и афтершоковую прогнозируемость.

Высокая прогнозируемость сейсмического потока (сохранение аппроксимационных тенденций при их экстраполяции в будущее, см. табл. 1) свидетельствует, что уравнение (2) при использовании описанной методики адекватно моделирует динамику сейсмичности по всем рассмотренным потоковым характеристикам. Это делает возможным применение уравнения (2) для дифференцирования потоковых характеристик сейсмичности, т.е. для определения скоростей и ускорений сейсмического потока. В свою очередь анализ пространственного распределения этих производных позволяет определить локализацию имеющихся тенденций к изменению сейсмичности. Поскольку сильные землетрясения приурочены к экстремумам прогнозируемой нелинейности сейсмического потока (см. рис. 1), то пространственные локализации с экстремальными тенденциями к его изменчивости представляют интерес для прогноза как самих сильных землетрясений, так и афтершокового затухания сейсмичности, т.е. именно в этих локализациях имеют смысл экстраполяции существующих тенденций в будущее. По сути мы получаем комбинацию методов саморазвивающихся процессов и картирования сейсмической активности по плотности сейсмического потока [Тихонов, 2006, 2009].

Потенциальную эффективность изложенной методики отчасти подтверждают ранее полученные данные [Малышев, 2014]: на картах ускорений сейсмической активности, построенных при сканировании с шагом 10 км (гипоцентральный радиус выборок 15 км) зона предстоящего Восточно-Японского землетрясения (11.03.2011 г., M = 9.0) начинает отчетливо выделяться за несколько дней до толчка. Зона подготовки имеет размеры примерно 100 × 80 × 40 км, что существенно превышает шаг сканирования. Это также согласуется с результатами данной работы (см. табл. 2), свидетельствующими о том, что прогнозируемость сильных землетрясений проявляется на средних (от 7.5 до 30 км) и больших (60 и 150 км) радиусах гипоцентральных выборок. Таким образом, предлагаемая методика устанавливает, как минимум, наличие перед Восточно-Японским землетрясением хорошо выраженных предвестников, использование которых позволило бы расширить интервал времени на предотвращение разрушительных последствий землетрясения с 10–30 мин (время достижения волной цунами побережья) до нескольких суток, что соответствует прогнозу по “сценарию форшоков” [Федотов и др., 2012; Федотов, Соломатин, 2015].

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

Тем не менее, для перехода от анализа прогнозируемости сейсмического потока к прогнозу собственно сильных землетрясений необходимо предварительно решить две проблемы: 1) отработать и верифицировать алгоритм автоматического выделения зон максимальной изменчивости сейсмического потока, что позволило бы уже сейчас перейти к оперативному сейсмическому мониторингу с отслеживанием гипоцентральных зон с максимальными тенденциями к активизации/затуханию сейсмического потока; 2) попытаться установить статистическую связь пространственных производных сейсмического потока с максимальной энергией землетрясений в их прогнозных зависимостях. Последнее позволяет определить допустимый предел возможных экстраполяций. Наиболее перспективны экстраполяции в будущее тенденций к увеличению выделяющейся сейсмической энергии E, позволяющие рассчитать в какой именно промежуток времени выделение сейсмической энергии может достичь потенциально опасных значений. Предельный диапазон выделения сейсмической энергии в определенный период времени можно оценить по вертикальной составляющей полосы ошибок на экстраполяционной части прогнозной зависимости. Аналогичная экстраполяция суммарного количества землетрясений позволяет рассчитать на этот же период времени ожидающееся число сейсмических толчков. Соотношение этих величин дает расчетное значение вероятной энергии одиночного сейсмического события.

При использовании данных Всемирного каталога Геологической службы США мониторинг на основе метода саморазвивающихся процессов с предварительным картированием нелинейности сейсмического потока для Южной Европы и Средиземноморья имеет смысл при радиусах гипоцентральных выборок от 7.5 до 60 км. При этом с повышением радиуса выборок увеличивается число сильных землетрясений, имеющих прогнозируемость по потоку сейсмической энергии, тогда как уменьшение радиуса приводит к возрастанию уровней прогнозируемости и прогнозной нелинейности потока энергии, а также к повышению точности определения пространственного положения экстремумов нелинейности потока сейсмической энергии, к которым приурочены сильные землетрясения.

ЗАКЛЮЧЕНИЕ

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