Концепция прогнозирования сейсмособытий на основе экспериментальных данных различной природы Текст научной статьи по специальности «Математика»

Известия Тульского государственного университета Естественные науки. 2009. Вып. 2. С. 147-158

= Информатика =

УДК 669.187

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

И.А. Афанасьев

Аннотация. Предложена концепция прогнозирования сейсмособытий с использованием данных приборной системы, а также метеорологических и геофизических данных. С помощью метода «Гусеница»-88А решается задача о прогнозировании времени будущих землетрясений. Проводится анализ геофизических и метеорологических данных для определения района предстоящих сейсмособытий.

Ключевые снова: временной ряд, метод «Гусеница»-88А, статистический анализ.

Введение

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

1) Мстеоданные. Влияние погоды на ссйсмоактивность и наоборот не вызывает сомнений. Эта взаимосвязь была замечена людьми тысячелетия назад. Например, во многих случаях перед началом толчков очевидцы замечали появление вдруг странного тумана, который, как правило, бывает только по утрам, да и то только при определенной влажности и температуре воздуха.

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

отклонение температуры воздуха от нормы, обычно в большую сторону — резкое наступление жары, засухи.

2) Данные о границах сейсмоактивных регионов.

3) Данные о разломах земной коры, особенно об активных. Как было сказано выше, землетрясения происходят с большей долей вероятности именно в зоне активных земных разломов.

4) Основные данные для исследования и прогнозирования — данные измерений геодинамических поляризационных потенциалов в лаборатории кафедры ОТ-СиЛП ТулГУ — экспериментальные сигналы.

Процесс прогнозирования сейсмособытий состоит из двух этапов: определение районов будущих значимых (с магнитудой более 6 баллов) сейсмособытий, и определение временной задержки, т.е. времени, через которое эти события наступят. Понятно, что точно определить эту задержку невозможно. Поэтому, она будет выражаться как I ± А1. Очевидно, что чем выше будет точность самого метода, чем больше мы будем учитывать значимых факторов, тем меньше будет в итоге А1, и, следовательно, выше точность прогноза.

Исходными данными для определения временной задержки будущих землетрясений являются данные измерений геодинамических поляризационных потенциалов в лаборатории кафедры ОТ-СиЛП ТулГУ. Они представляют собой временной ряд. Но, по сути, мы имеем не один, а два ряда: первый — экспериментальный сигнал, второй — ряд фактических землетрясений, точек на оси времени. Если рассматривать последний как временной ряд с тем же шагом по оси времени, что и первый, а также учесть тот факт, что экспериментальный сигнал каким-то образом отображает «предвестники землетрясений», виде «необычного поведения» временного ряда, то можно рассматривать не просто два параллельных временных ряда, а один двумерный [2]. Т.е. задача исследования сводится к анализу двумерного временного ряда. Для этого будем использовать метод «Гусеница»-88А . Именно он, в отличие от всех остальных, отлично справляется с задачей обработки двумерного временного ряда [1]. Рассмотрим подробно метод «Гусеница»-88А для анализа многомерных временных рядов.

1. Базовый алгоритм метода «Гусеница»-88Л

Пусть имеются два ряда ^ = (/0...../лг-О иС=(й.........gN-l) длины

N > 2. Их можно рассматривать как двумерный временной ряд (_Р, С) или как комплексный ряд С = (со,... , елт-1), где сд, = Д + ígk, 0 ^ к < N. Для анализа временных рядов выберем целый параметр Ь, 1 < Ь < N. Его называют «длина гусеницы» или «длина окна». Схему одномерного метода рассмотрим для ряда F. Далее будем предполагать, что ряды ненулевые.

Базовый алгоритм состоит из четырех шагов. Первые два шага образуют этап разложения, а последние два — этап восстановления.

1.1. Вложение. Построим последовательности Ь-мерных векторов для рядов ^ и (7:

К = (/¿-Ъ • • • , /¿+Ь-2)Т, Ог = . . . Jgi+L-2)rГ, 1^^К.

В каждой из последовательностей число векторов равно К = N — Ь + 1. Траекторией матрицей двумерного ряда, порожденной длиной окна Ь, назовем матрицу

X = : ... : Рк: : ... : Ск] = [Е, С].

Другими словами, траекторная матрица двумерного ряда представляет собой расположенные последовательно траекторные матрицы Е и С одномерных рядов ^ и С, полученные при одинаковой длине окна Ь. Таким образом, размерность получившейся матрицы равна Ь х 2К.

Ясно, как изменится первый шаг метода при произвольном числе рассматриваемых рядов произвольной (не обязательно одинаковой) длины. Будем

считать, что наблюдается система из 8 временных рядов ,

где к = Параметр А^, таким образом, есть длина к-го ряда. Вы-

бираем Ь такое, что 1 < Ь < Nk для любого к. Для каждого к вычислим Ки = N/1 — Ь + 1 векторов

Х',1" = (/м....../йи)''' 1

Тогда траекторная матрица многомерного ряда ... , будет иметь

вид

Х =

Л1 ... К\....Л1 ....лка

Х™:...: Х«

где Х^ — траекторная матрица ряда соответствующая длине окна Ь.

Х

к=1

дов и длине окна Ь траекторная матрица однозначно определяет ряды, по которым она построена.

1.2. Сингулярное разложение. Рассмотрим матрицу Э = ХХТ, где

Х

менного ряда, соответствующая длине окна Ь.

Поскольку Э положительно полуопределепа, ее собственные числа неотрицательны. Обозначим через Аь ... , Ль собственные числа матрицы Э, взятые в порядке убывания (Ах ^ ^ Аь ^ 0) и через и±,..., 1!ь ортонор-

Э

этим собственным числам. Пусть (1 = тах{г, таких, что А/ > 0}. Обозначив V* = =Хт?7г (г = 1,... , е£), получим разложение траекторной матрицы:

Х Х Х

где Х* = \/\1 и^.

Разложение (1) и является сингулярным разложением траекторией матрицы X. Отмстим, что ортонормированные векторы V* являются собственными векторами матрицы ХТХ, соответствующими тем же собственным числам А*. В стандартной терминологии \/Х7 называются сингулярными числами, £/*

Х

Х

элементарными матрицами. Набор (\/А7, V-,) мы будем называть г-й соб-

ственной тройкой сингулярного разложения.

На геометрическом языке система собственных векторов задаст ор-тонормированный базис в линейном пространстве, порождаемом столбцами Х

Х

В другой (статистической) терминологии векторы называются

собственными, векторы V-* — факторными, направление, задаваемое г-м собственным вектором иг — г-м главным направлением, вектор Zi = \/А7V-*, составленный из проекций векторов X* на г-е главное направление — вектором г-х главных компонент.

Прокомментируем результат этапа разложения для двумерного случая. При таком подходе Э = ХХТ = ЕЕТ + ООт,

14 = Тг, (ст)' и1 = (^)

и

<1 <1

Х _ тт ЦтТ-п ТтТс\_^ ГП цТ-п . л ттТс

Х — ик\ик п 1 ик с) — [ик^к п ■ ^к^к ° ■

к=1 к=1

Таким образом, для каждого из рядов получено разложение столбцов их траскторных матриц по общему базису (II 1,..., и а)'-

(1 (1

ЕПГ1 ———

ики№, С = ^,ики?С. (2)

к=1 к=1

Эти разложения, однако, не обязаны быть сингулярными разложениями матриц

1.3. Группировка. Процедура группировки формально одинакова для всех рассматриваемых разновидностей БвА. На основе разложения (1) процедура группировки делит все множество индексов {1,..., на га нспсрс-сскающихся подмножеств 1\,..., 1т.

Х

Х Х Х

вычисляются для I = /1,...,/т, тем самым разложение (1) может быть записано в сгруппированном виде:

Х Х Х

(3)

Процедура выбора множеств 1\,... , 1т и называется группировкой собственных троек.

1.4. Диагональное усреднение. На последнем шаге базового алгоритма каждая матрица сгруппированного разложения переводится в новый ряд длины N. Пусть элементами y-ij, 1 ^ L.

1 ^ j < К. Положим L* = min (L, К), К* = max (L, К) nN = L+ K-l. Пусть Zij = yij, если L < К и Zij = yji в остальных случаях.

Диагональное усреднение переводит матрицу Y в ряд (go,... по

формуле

1

k+1

gk =

L*

j,k-j+2

для 0 ^ к < L* — 1, для L* - 1 ^ к < К*

(4)

i=і 1

N - к

N-K*+l

Zjik^j+2 ДЛЯ К* ^ к < N.

j=k-K*+2

Это выражение соответствует усреднению элементов матрицы вдоль «диагоналей» I + ] = к + 2: выбор к = 0 даст ^о = уц, для к = 1 получаем ё! = (у 12 + г/21)/2 и т.д. Применив диагональное усреднение к матрицам Х]к. полученным на этапе группировки, приходим к разложению исходного ряда в сумму т рядов.

Рассмотрим случай для двух исходных рядов. Каждая матрица Х]к размера Ь х 2А' в сгруппированном разложении (3) разбивается на послсдова-

Х(Я : X(r2)

После

тельно расположенные матрицы размера Ь х К: Х/к =

(?)

Х

пенис по формуле (4). В результате каждое слагаемое в правой части (3) порождает двумерный ряд — аддитивную компоненту ряда (.Р, С).

Таким образом, результатом применения базового алгоритма ББА к (одномерному, многомерному или комплексному) ряду является его представление в виде суммы т рядов. Параметрами метода являются длина окна Ь и способ группировки элементарных матриц [3].

Итак, после диагонального усреднения и последующего продолжения, появляется возможность построения прогноза двумерного ряда. Хотя нас интересует прогноз только второй его составляющей — временного ряда фактических сейсмособытий [4]. При рекуррентном прогнозировании диагональное усреднение используется для получения восстановленного ряда, к последним точкам которого затем применяется линейная рекуррентная формула.

Таким образом, с помощью метода «Гуссница»-88А решается задача о прогнозировании времени будущих землетрясений.

2. Определение района сейсмособытий

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

Под эгидой «Глобальной программы по предупреждению землетрясений» (GSHAP) была разработана карта сейсмически опасных районов Земли [7]. Вся суша Земли на ней разделена на районы, классифицированные по степени опасности ссйсмособытий. Для последующего анализа и обработки данных мы введем разбиение сейсмически опасных районов на отдельные зоны. Таким образом, для последующего рассмотрения мы получим N сейсмоактивных зон.

Также в пашем распоряжении имеется карта поверхности Земли с обозначением тектонических плит и разломов, составленная компанией National Geographic [6].

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

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

Данные по температурам Т поступают постоянно в режиме реального времени и автоматически сравниваются с данными со средними показателями (по соответствующим точкам земли). Под точкой понимается целый район, область. При возникновении ситуации, при которой Т — Тср > Тпорог. где Тп0р0г — пороговое значение (определяемое экспериментально), данная точка запоминается — проецируется на карту, ту самую, на которой нанесены сейсмоопасныс районы и активные линии разломов. Будем называть такие точки — «красные точки» (далее — КТ). По сути, Тпорог — одна из тех характеристик, которая делает вклад в определение будущего района предполагаемого землетрясения. Очевидно, такие красные точки могут появляться теоретически сразу в нескольких районах земли, а могут не появляться вовсе.

Так как данные по температурам поступают к нам в режиме реального времени (раз в день),то отклонение температуры в отдельно взятой точке поверхности земли от Тп0р0г может быть эпизодическим, случайным, а может и продолжаться в течение нескольких дней (некоторого периода времени), поэтому целесообразно ввести еще одну характеристику — Е интенсивность

красной точки, т.с. сс длительность существования. Далее происходит процесс оценки положения, степени близости красных точек на карте относительно положения сейсмоактивных областей поверхности земли (лежит ли КТ в одной из таких областей или пет) и относительно линий активных разломов. Интуитивно понятно, что чем ближе красные точки к разломам, тем с большей степенью вероятности в указанных областях через определенный промежуток времени I произойдут землетрясения. Ясно, что значимые степени близости КТ к указанным областям на разных континентах могут оказаться своими.

Итак, если мы определили КТ, удовлетворяющую необходимому пороговому значению, определенной степени интенсивности Е, если положение данной КТ лежит на территории сейсмоактивной зоны, удовлетворяет степени близости Д к активной линии разлома (ближайшей), то район (границы) предполагаемого землетрясения определяется следующим образом: на карте строится область по внешним границам трех составляющих — КТ, участка лини разлома, сейсмоактивного района. В упрощенном варианте для простоты программной реализации можно границы предполагаемого землетрясения определять по границам сейсмоактивной зоны.

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

В дисперсионном анализе исследуется влияние одного или нескольких качественных показателей на количественный [5]. В однофакторном дисперсионном анализе на одну количественную переменную У оказывает влияние один фактор (один качественный показатель), наблюдаемый на к уровнях [5].

В нашем случае качественным показателем (фактором) является наличие/отсутствие в данной зоне тектонического разлома. Количественным показателем является число землетрясений в сейсмически активной зоне за определенный период времени. Наблюдаемые данные обозначим (число

сейсмоактивных зон), где г - индекс уровня (г = 1,... , к); — индекс наблю-

дения на г-м уровне, у = 1,..., п* (на каждом уровне может быть свое число наблюдений).

к

Общее число опытов (наблюдений) п = ^ п-1-

¿=1

По данным можно определить следующие характеристики:

. . _ 1 П{

У г — среднее значение переменной У на г-м уровне yi = — ^

Щ 3=1

У — среднее значение переменной У по всем значениям

, к гц 1 к

* = 1 j = 1 * = 1

5 — сумму квадратов отклонений наблюдений от общего среднего

к щ

¿=1 ¿=1

Бр — сумму квадратов отклонений средних групповых значений от общего среднего

Зр = '¡Г (у* - У)2'Пг'->

¿=1

Бост — сумму квадратов отклонений средних групповых значений от общего среднего

к гц

ОСТ

¿=1 ¿=1

ОСТ

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

Далее выполняются следующие действия: вычисление общей дисперсии

к тц

°2общ=—[=—[ ё ^ - у)2;

¿=1 ¿=1

вычисление факторной дисперсии

к

"Г = = ГЗТ Т, <-У‘ -

п — 1 к — 1

¿=1

вычисление остаточной дисперсии

-* к П{

2 ООСТ I -\2

а0СТ~—[ -:^ЪЪ(уц-у) ■

¿=1 ¿=1

В дисперсионном анализе проверяется гипотеза Но о равенстве средних групповых значений количественного показателя:

(Н0 :у1 = у2 = ... = Ук).

Критерием для проверки этой гипотезы является соотношение факторной дисперсии и остаточной дисперсии:

ОСТ

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

2.2. Определение влияния температуры. Вышеизложенный метод прогноза временных рядов (ББА) даст нам ряд прогнозных значений на определенный промежуток времени вперед. Т.с. фактически мы имеем «предсказание» сейсмособытий по всему Земному шару. Но для определения фактического места землетрясения нам необходима дополнительная информация. Как было сказано ранее, имеются многочисленные наблюдения о том, что в преддверии землетрясения наблюдались нестандартные погодные явления. Именно температурные отклонения могут служить этой дополнительной информацией, по которой мы сможем идентифицировать место (район) предстоящего землетрясения.

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

т "

А ср

Именно по данному показателю мы будем принимать или не принимать во внимание тс или иные сейсмоактивные зоны.

Для определения влияния аномальных отклонений температур на ссйсмособытия воспользуемся корреляционным анализом. Выше было введено понятие интенсивности (Е) КТ. Поэтому будем рассматривать только тс отклонения, которые продолжались не менее определенного времени (1кр).

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

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

т

'Р = ~, и

где и — общее число землетрясений в сейсмоактивной зоне; т — число землетрясений с магнитудой более 6 баллов при отклонении температуры в сейсмоактивной зоне.

Можно также определить доверительные границы для данной вероятности по формуле:

где ¿ф = а^ Ф(Рд/%) — обратная функция Лапласа; Р$ = 0,9.

Таким образом, мы получаем вероятности ссйсмособытия с магнитудой выше 6 баллов для каждой рассматриваемой сейсмоактивной зоны.

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

Здесь п — общее число землетрясений в зонах с разломом/без разлома, т — число землетрясений более 6 баллов при отклонении температуры в зонах с разломом / без разломов. Т.с. мы фактически получим две вероятности (вероятности ссйсмособытия выше 6 баллов при аномальном отклонении температуры через время Ьр ± Д£): Рр — для зон с разломами и Рбр — для зон без разломов.

Теперь окончательный метод прогнозирования ссйсмособытия выглядит следующим образом. Для прогнозирования ссйсмособытий у нас для обработки есть виды данных:

1. Прогнозные данные, полученные при анализе двумерного временного ряда методом ББА. Данный метод прогноза временных рядов даст нам

тр I бр

ряд прогнозных значений числа сейсмособытий на определенный промежуток времени вперед. Т.с. фактически мы имеем «предсказание» сейсмособытий по всему Земному шару.

2. Метеорологические данные (Красные точки). То есть данные по аномальным отклонениям температуры в рассматриваемых сейсмически активных зонах. Они будут использоваться для определения района будущих землетрясений.

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

Таблица 1. Возможные варианты прогнозов сейсмособытий

Прогноз методом «Гусеница^-БЭЛ Наличие КТ

да нет

да В1 В2

нет ВЗ В4

Рассмотрим каждый вариант подробно.

Вариант 1 (В1). В сейсмоактивной зоне есть красная точка и метод ,8,8А фиксирует предстоящее сейсмособытие.

В этом случае считается, что в данной зоне (или зонах, если красных

р

бр

Вариант И (В2). Красных точек нет. метод ,8,8А фиксирует предстоящее сейсмособытие.

В этом случае считается, что сейсмособытис произойдет с определенной

р бр

активных зонах.

Вариант 3 (ВЗ). В сейсмоактивной зоне есть красная точка, а метод ,8,8А не фиксирует предстоящее сейсмособытие.

В этом случае считается, что в данной зоне (или зонах, если красных точек несколько) сейсмособытия не произойдут, поскольку не каждое отклонение температуры может быть связано с землетрясением.

Вариант 4 (В4). В сейсмоактивной зоне нет КТ и метод ,8,8А не фиксирует предстоящее сейсмособытие.

В этом случае считается, что сейсмособытия нигде не произойдут.

Таким образом, была разработана концепция определения времени предстоящих сейсмособытий на основе анализа временных рядов методом «Гуссница»-88А, а также проведен статистический анализ геофизических и метеорологических данных для определения района (или районов) землетрясения.

Список литературы

1. Голяндина Н.Э. Метод «Гусеница^-SSA: прогноз временных рядов. Учеб. пособие. СПб.: СПбГУ, 2004. 78 с. ‘

2. Чистова Г.К. Модели и методы обработки сейсмических сигналов в системах распознавания. Пенза: ПГУ, 2004. 230 с.

3. Eisner J., Tsonis A. Singular Spectrum Analysis. A New Tool in Time Series Analysis. N.Y.: Plenum Press, 1996. 163 p.

4. Данилов Д.Л. Метод «Гусеница» для прогнозирования временных рядов. СПб.: Пресс ком, 1997. 213 с.

5. Гмурман В.Е. Теория вероятностей и математическая статистика. М.: Высшая школа, 1972. 355 с.

6. Официальный сайт National Geographic: www.nationalgeographic.com.

7. Сайт национальной геологической службы США: www.usgs.gov.

Поступило 05.06.2009

Афанасьев Игорь Александрович ([email protected]), аспирант, кафедра прикладной математики и информатики, Тульский государственный университет.

The earthquakes prediction concept based on experimental data of different nature

I.A. Afanasiev

Abstract. We perform the earthquakes prediction concept using the data from the instrument system, as well as meteorological and geophysical data. Using the method of «Catcrpillar»-SSA wc solve the problem of predicting the time of future earthquakes. Using the analysis of geophysical and meteorological data wc determine the area of upcoming earthquakes.

Keywords : time scries, the method of « Caterpillar »-SS A, a statistical analysis.

Afanasiev Igor ([email protected]), postgraduate student, department of applied mathematics and computer sciences, Tula State University.