способ обнаружения очагов землетрясений сетью сейсмостанций

Формула изобретения

Способ обнаружения очагов землетрясений сетью сейсмостанций, включающий формирование массивов дискретных отсчетов сигнала отношения скоростей продольных (V_P) и поперечных (V _S) волн от сейсмостанций при различных азимутах направлений трасс измерений, визуализацию массивов измерений в виде диаграмм Пуанкаре и определение по ним характера измеряемого процесса, отождествление момента появления в диаграмме регулярной составляющей с началом сейсмического процесса, пеленгацию обнаруженной аномалии несколькими сейсмостанциями путем формирования их диаграмм направленности двумя сейсмоприемниками из трехкомпонентных кондуктометрических датчиков давления, размещенных на базе размером, кратным длине волны измеряемого сигнала, перпендикулярной трассе измерений, определение гипоцентра очага землетрясения как точки пересечения радиус-векторов скоростей сейсмоволн от нескольких сейсмостанций, представление массивов измерений сигнала в виде феноменологического разностного уравнения отслеживаемого процесса, вычисление коэффициентов феноменологического уравнения и динамики изменения энергии сигнала по коэффициентам W_i(t)=a²+b² +c+d от выборки к выборке, определение постоянной времени Т _i динамического процесса и установившегося значения сигнала W_{уст i} W_0ie^-4,7, расчет результирующего затухания сигнала на трассе W_i W_0i-W_{уст i} и относительного затухания K_i= W_i/W_{уст i}, прогнозирование времени ожидаемого удара t_y «4,7T и магнитуды M=r·ln[K ],

где t - интервал времени между выборками измерений;

- среднегеометрическая постоянная времени измеряемого процесса по всем трассам;

- среднегеометрическое значение относительного затухания сигнала по всем трассам;

W_0i - энергия хаотического процесса на измеряемой трасе;

n - число станций в сети;

r - коэффициент регрессии;

а, b, с, d - расчетные коэффициенты феноменологического уравнения.

Описание изобретения к патенту

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

Достоверный прогноз землетрясений возможен при установлении их первопричины. Одной из первопричин, вызывающих землетрясения, является накачка земной коры дополнительной энергией восходящего потока газов (Н₂, Не, СН₄). Факт эманации газов из земной коры накануне сейсмического удара [см., например, патенты RU № 2204852, 2003 г., № 2275659, 2006 г., № 2302020, 2007 г.]. Накачка земной коры дополнительной энергией приводит к раскачке очага землетрясений, сопровождаемой распространением от гипоцентра сверхнизких литосферных волн [см. Научное открытие № 365, «Явление раскачки очага землетрясений перед сейсмическим ударом», РАЕН, 2008 г.].

Известен способ предсказания землетрясений, патент RU № 2170446, 2001 г., путем измерения сверхнизких литосферных волн раскачки очага землетрясений космическими средствами системы Navstar (GPS) - аналог.

В способе-аналоге размещают в сейсмическом районе приемные станции космической навигационной системы, разнесенные на протяженной измерительной базе, осуществляют непрерывное высокоточное измерение координат точек размещения приемных станций, регистрируют момент появления периодических отклонений x_i, y_i, z_i координат точек и отслеживают изменения этих отклонений во времени, вычисляют гипотетический центр очага как точку пересечения радиус-векторов в пространстве, длину и направляющие косинусов, которые определяют из соотношений: , рассчитывают время удара, отсчитываемое от момента появления периодических отклонений координат: и магнитуду , где Т - период отклонения координат, ч, dek r - натуральный логарифм отношения амплитуд отклонения координат двух смежных периодов; d, l - коэффициенты регрессии.

Недостатками способа-аналога являются:

- наличие скрытого участка нечувствительности (мертвой зоны), ограниченной среднеквадратичной ошибкой средств GPS, снижающей интервал прогнозирования;

- неточность регрессивных зависимостей расчета времени удара и магнитуды ввиду зависимости периода (Т) отклонения координат от времени.

Ближайшим аналогом к заявленному способу является способ оперативного прогноза землетрясений, патент RU № 2353957, G01V, 9/00, 2009 г.

Способ ближайшего аналога включает регистрацию сейсмических волн в виде дискретных отсчетов амплитуд (x_k, y_k) сигналов во взаимно ортогональных плоскостях в двух разнесенных на измерительной базе пунктах, причем оси чувствительных датчиков пунктов по координате x ориентируют по направлению базы обработку зарегистрированных данных путем разложения в функциональный ряд, вычисление автокорреляционных функций сигналов, дополнительно, регистрацию естественного акустического шума Земли в виде акустограмм посредством геофонов, располагаемых в глубоких скважинах, визуализацию точечного множества отсчетов последовательности данных в виде диаграмм Пуанкаре в плоскостях (x_k, x_k+1) (y_k, y_k+1 ), отслеживание динамики изменения формы диаграмм Пуанкаре, определение постоянной времени переходного процесса (Т), прогноз времени ожидаемого удара t_y 4,7T, магнитуды (М) из соотношения lgt_y[сут]=0,54M-3,37, а координаты гипоцентра очага отождествляют с точкой пересечения направляющих косинусов, отсчитываемых от оси базы: где B_x1(0), B_y1(0), B_x2 (0), B_y2(0) - значения автокорреляционных функций сигналов в нуле в плоскостях x, y геофонов соответственно первого и второго пунктов.

Недостатками ближайшего аналога являются:

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

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

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

Поставленная задача решается тем, что способ обнаружения очагов землетрясений сетью сейсмостанций включает формирование массивов дискретных отсчетов сигнала отношения скоростей продольных (V_P ) и поперечных (V_S) волн сейсмостанций при различных азимутах направлений трасс измерений, визуализацию массивов измерений в виде диаграмм Пуанкаре и определение по ним характера измеряемого процесса, отождествление момента появления в диаграмме регулярной составляющей с началом сейсмического процесса, пеленгацию обнаруженной аномалии несколькими сейсмостанциями путем формирования их диаграмм направленности двумя сейсмоприемниками из трехкомпонентных кондуктометрических датчиков давлений, размещенных на базе, с размерами кратной длине волны измеряемого сигнала, перпендикулярной трассе измерений, определение гипоцентра очага землетрясения как точки пересечения радиус-векторов скоростей сейсмоволн нескольких сейсмостанций, представление массивов измерений сигнала в виде феноменологического разностного уравнения отслеживаемого процесса, вычисление коэффициентов феноменологического уравнения и динамики изменения энергии сигнала по коэффициентам W_i(t)=a²+b² +c+d от выборки к выборке, определение постоянной времени T _i динамического процесса и установившегося значения сигнала W_ycтi W_{0 i}·l^-4,7, расчет результирующего затухания сигнала на трассе W_i W_{0 i}-W_{ycт i} и относительного затухания K_i= W_i/W_{ycт I}, прогнозирования времени ожидаемого удара t_y 4,7T и магнитуды M=r·ln[K ], где

t - интервал времени между выборками измерений;

- среднегеометрическая постоянная времени измеряемого процесса по всем трассам;

- среднегеометрическое значение относительного затухания сигнала по всем трассам;

W_0i - энергия хаотического процесса на измеряемой трассе;

n - число станций в сети;

r - коэффициент регрессии;

а, b, с, d - расчетные коэффициенты феноменологического уравнения.

Изобретение поясняется чертежами, где

фиг.1 - схема пеленгации очагов землетрясений сетью сейсмостанций Камчатского геофизического полигона;

фиг.2 - вид функции регистрируемого сигнала (V_p/V _s) в интервале наблюдений;

фиг.3 - визуализации диаграмм Пуанкаре функций сигналов: а) случайного сигнала, б), в), г) с регулярной составляющей;

фиг.4 - графическая интерпретация поля решений феноменологического уравнения: а) устойчивого состояния, б), в), г) накачки (сброса) энергии в среде на измеряемых трассах;

фиг.5 - динамика изменения функций энергии сигналов на измеряемых трассах;

фиг.6 - динамика изменения затухания сигнала от размеров зоны подготавливаемого землетрясения;

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

Техническая сущность изобретения заключается в следующем. Одной из первопричин землетрясений является дегазация земной коры и ее накачка дополнительной энергией восходящего потока газов высокого давления. Насыщение земной коры газовой компонентой изменяет ее вязкоупругие характеристики. Дисперсия плотности породы приводит к дисперсии скоростей распространения литосферных волн, т.е. изменению соотношения между скоростями распространения продольных V_p и поперечных V_s волн. Скорость распространения продольных волн в среде поперечных где Е - модуль упругости Юнга, G - модуль сдвига, - плотность среды [см. Курс физики. Г.А.Зисман, О.М.Тодес. Т.1, стр.301-302. М., Изд. Наука, 1964 г.]. В качестве сигнала предвестника изменения вязкоупругих характеристик среды выбрано отношение V_p/V_s, изменяющееся во времени и в пространстве от азимута трассы измерений. Одна из реализаций сигнала предвестника иллюстрируется графиком фиг.2. Известно [см., например, Р. Дуда, П.Харт. «Распознавание образов и анализ сцен». стр.319-331, перев. с англ., М., Мир, 1976 г.], что максимум информации о происходящих процессах содержится в образе объекта.

Скрытая информация о динамике процесса на его ранней стадии может быть извлечена из электрического сигнала сейсмического фона, если математическими процедурами обработки визуализировать и отслеживать изменение формы сигнала (образа объекта). Для этого программным методом последовательность данных отображают в виде диаграмм Пуанкаре [см., например, А.Пуанкаре. Избранные труды в 3-х т. Пер. с французского под ред. Н.Н.Боголюбова, Наука, 1974 г.]. Диаграмма Пуанкаре представляет собой точечное графическое отображение N-значений последовательности, например X_к при к=1, 2, 3 N на двумерном поле, в котором ординатой точки является значение Х_к+1, а абсциссой - предшествующее значение Х_к. Нанося поочередно точки для к=1, 2, 3 N на график, получают точечное множество Х_к+1 , (Х_к), образующее фигуру, по которой можно судить о типе последовательности (образе объекта).

Текст программы визуализации диаграмм Пуанкаре

На фиг.3 иллюстрируются диаграммы Пуанкаре для различных классов сигналов:

а - Белый гауссовский шум, отсутствие информации;

б - неустойчивое состояние, хаотическое колебание;

в - устойчивое состояние среды, сходящаяся последовательность измерений;

г - катастрофическое состояние среды, расходящаяся последовательность (сброс энергии).

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

Момент появления в сейсмическом фоне регулярной составляющей отождествляют с возникновением на измеряемой трассе аномальной зоны подготавливаемого землетрясения. После обнаружения аномалии на трассе измерений осуществляют ее пеленгацию. Известно [см., Г.Корн, Т.Корн. «Справочник по математике для научных работников и инженеров» разд. «Аналитическая геометрия», стр.73-74, перевод с англ., М, Наука, 1971 г.], что положение радиус-вектора в пространстве полностью определяется его косинус-направляющими. Из аналитической геометрии известно, что косинус-направляющая вектора равна отношению его проекции на данную ось к длине вектора. Длина вектора полной скорости (дл.V) находится, как корень квадратный из суммы квадратов его проекций: где V_x, V_y - скорости продольных волн для данной сейсмостанции по координатам x, y. Для каждой сейсмостанции определяют направление . Координаты аномалии находят как точку пересечения радиус-векторов скоростей от нескольких сейсмостанции. После обнаружения аномалии и ее пеленгации на измеряемых трассах осуществляют прогноз параметров ожидаемого сейсмического удара. Возникновение на измеряемой трассе газонасыщенного объема земной коры изменяет параметры регистрируемого сигнала. Динамика изменения параметров сигнала содержит всю информацию о предстоящем сейсмическом ударе. Для извлечения скрытой информации из последовательности выборок измерений используют их приближение феноменологическим уравнением в виде:

x_k =f(x_k-1, x_k-2, , x_k- , a)+ _k; f F.

Здесь f(.) - искомая функция, принадлежащая некоторому выбранному классу функций F; - параметр состояния, x_k=x(t_k) - наблюдаемый параметр системы, t_k, k=1, 2 N - дискретное время, N - объем выборки.

Метод разработан в классе степенных рекуррентных алгебраических полиномов [см., например, О.Н.Новоселов. «Идентификация и анализ динамических систем». Изд. МГУЛ, М, 2006 г., стр.34-44]

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

Коэффициенты a, b, c, d находились из условия минимума функции среднего квадрата отклонения который определяется путем приравнивания нулю частных производных:

Эта система представляет собой систему линейных алгебраических уравнений относительно искомых параметров a, b, c, d, где коэффициенты при неизвестных имеют вид:

Здесь N - количество элементов последовательности x_k; , =0..4; p, q=0..2; m[.] - оператор математического ожидания. Решение системы линейных уравнений осуществлялось специализированной программой на основе метода Гаусса.

Текст программы реализации алгоритма

Возникновение на измеряемой трассе газонасыщенного объема земной коры подготавливаемого землетрясения изменяет параметры регистрируемого сигнала. В исходном (устойчивом) состоянии соотношение скоростей V_p/V_s составляет (фиг.2) порядка 1,4. Насыщение земной коры газовой компонентой изменяет вязкоупругие характеристики среды, при этом отношение V_p/V_s стремится к единице, а сам сигнал как бы пропадает, его дисперсия стремится к нулю.

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

Коэффициенты a, b при первых степенях амплитуд сигналов (x_к, x_к+1 ) имеют физический смысл постоянной составляющей.

Коэффициенты с, d при квадратичных значениях характеризуют энергию переменной составляющей сигнала (дисперсию).

Полная энергия сигнала равна квадрату постоянной составляющей плюс дисперсия, т.е. W_i(t)=[a²+b²+c+d].

Одна из реализаций функции энергии сигнала для станции MKZ в интервале наблюдений (неустойчивое состояние) иллюстрируется графиком фиг.5

Известно, что сама физическая величина и скорость ее изменения связаны дифференциальным уравнением первой степени, общим решением которого является экспоненциальная зависимость [см., Н.С.Пискунов. «Дифференциальное и интегральное исчисление для ВТУЗов». М.: Наука, 1976 г., стр.458]. Экспоненциальная зависимость обладает тем свойством, что по ее дискретным отсчетам можно восстановить всю функцию. Экспоненциальная зависимость затухания сигнала на измеряемой трассе иллюстрируется графиком фиг.5. Энергия сигнала уменьшается от W_0i хаотического (Гауссова) процесса до установившегося минимума W_ycт.i , значение которого С вероятностью 0,99 экспонента достигает предельного значения при соотношении t/T_i=4,7. Признаком-предвестником землетрясения является увеличение размеров зоны подготавливаемого землетрясения (аномальной зоны фиг.1). Известна примерная зависимость ожидаемой магнитуды сейсмического удара (М) от размеров зоны подготавливаемого землетрясения [см., «Краткосрочный прогноз катастрофических землетрясений с помощью радиофизических наземно-космических методов». Доклады конференции ОИФЗ им. О.Ю.Шмидта. М.: РАН, 1998 г., стр.9].

M ln[R, км].

За счет поглощения в газонасыщенном объеме (в первом приближении) можно считать, что энергия сигнала W_i на трассе распространения убывает пропорционально размеру R_i аномальной зоны. Абсолютную величину затухания сигнала на трассе рассчитывают как W_i=W_0i-W_ycт.i= W(R_i).

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

K _i= W_i/W_ycт.i.

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

Количественные значения расчетных параметров приведены ниже в примере реализации.

Пример реализации способа

Реализация способа осуществлялась на Камчатском Геофизическом полигоне РАН. Схема полигона иллюстрируется фиг.1. Полигон включает 6 сейсмических станций: Русская (RUS), Петропавловск (PET), Шипунский (SPN), Мыс Козлова (MKZ), Крутоберегово (КВТ), Авача (AVH).

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

Фокальная зона Камчатки представляет сочленение континентальной и океанической плит, последняя из которых осуществляет движение под континентальную плиту. К северу, в зоне мониторинга станций КВТ, имеет место сочленение этих плит с более северной, т.е. имеет место тройное сочленение, сейсмическая активность которого очень высока. Фокальная зона землетрясений - от севера Камчатки на юг, вплоть до Индонезии. В мониторинге среды станциями RUS, PET, КВТ, MKZ, SPN, AVH использованы данные, полученные только от местных землетрясений.

Функциональная схема системы измерений, реализующей способ, иллюстрируется фиг.7. Система измерений (фиг.7) содержит размещенные в Камчатском регионе шесть сейсмостанций 1 6. Каждая из сейсмостанций, для азимутов трасс измерений в фокальной зоне, содержит группу из двух трехкомпонентных сейсмоприемников, 7, 8, оси чувствительности одноименных компонент которых взаимно параллельны. Сейсмоприемники размещены на измерительной базе 9, кратной размерам средней длине волны измеряемого спектра. Базы располагают перпендикулярно измеряемым трассам. Сейсмоприемники устанавливаются заглубленными, соосно друг другу. Через встроенные усилители 10, 11, 12 сигналы с выхода сейсмоприемников подают на соответствующие входы канального коммутатора 13. Канальный коммутатор реализует функцию отношений сумм сигналов и сигнал пеленга аномалии (V_px/дл.V_p ). Результирующие сигналы от канального коммутатора 13 оцифровывают в АЦП 14, 15, записывают в буферное ЗУ 16 для последующего ввода данных в компьютер 17 в стандартном наборе элементов: процессор 18, оперативное ЗУ 19, винчестер 20, дисплей 21, принтер 22, клавиатура 23. Компьютеры 17 объединены в локальную вычислительную сеть, и обеспечена синхронизация их работы во времени.

Все элементы устройства представляют собой существующие технические разработки. В устройстве использованы элементы виброизмерительной аппаратуры фирмы Bruel & Kjairi ENDEVCO, Дания, следующих моделей: встроенные усилители 10, 11, 12 модели 2626, 2628; канальный коммутатор, АЦП, буферное ЗУ (элементы 13, 14, 15, 16) - многофункциональный блок 3560-L.

Сейсмическая волна характеризуется избыточным давлением Р[н/м²]=pc| |, где р - плотность среды, с - скорость акустических волн, - скорость колебательного движения, u - смещение частиц от положения равновесия в процессе колебательного движения. Поскольку давление во фронте сейсмической волны пропорционально скорости, для ее измерения на сейсмостанциях используют датчики давления. Известен класс устройств мембранного типа, позволяющих измерять инфразвуковые колебания земной коры (см., например, «Кондуктометрический датчик колебаний». Патент RU № 2055352, G01N, 27/02, 1996 г.). Датчик имеет высокую чувствительность и линейность выходной характеристики в области инфрачастот 0,1 100 Гц, быстродействие 100ms. Промышленная разработка кондуктометрического датчика представляет собой трехкомпонентный сейсмоприемник (см., например, «Способ определения гипоцентра нефтегазового месторождения». Патент RU № 2150719, G01V, 1/100, 2000 г.). Он может использоваться в любой пространственной ориентации, т.е. осуществлять измерения во взаимно ортогональных плоскостях.

Анализ сейсмических данных проводился для периода с 01.01.2001 по 31.12.2004 гг. Интервалы времени, в которые рассматривались изменения состояний среды, составляли ~16, 8, 4, 2 суток до состоявшегося землетрясения и соответственно 2, 4, 8, 16 суток после. Обработка сигналов на каждой станции проводилась по зарегистрированным сейсмическим волнам с учетом направления их прихода, т.е. выбирались различные азимуты (азимут может иметь значение от 0 до 360).

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

Функция изменения энергии сигнала на трассе станции MKZ во всем интервале выборок измерений иллюстрируется графиком фиг.5.

Расчетная энергия гауссовского шума равна W_0i=0,42 мВт.

Постоянная времени сейсмического процесса

Установившееся значение энергии сигнала W_ycт=W_0i·e^-4,7=0,0038.

Абсолютное поглощение энергии сигнала аномальной зоной W=W_0i-W_ycт=0,416.

Относительная величина затухания сигнала .

Ожидаемая магнитуда сейсмического удара: M₁=r·ln|K₁|=1,2 ln|112|=5,7.

Ожидаемое время сейсмического удара t_y=4,7T₁ 3,9 сут.

Аналогичные расчеты для тех же временных интервалов станции RUS дали результаты: К₂ =89, М₂=5,4, t_y=3,4 сут.

Среднегеометрическая величина постоянной времени сейсмического процесса для 2-х сейсмостанций .

Среднестатистическое затухание сигнала в аномальной зоне , M =5,55 баллов.

Поскольку напрямую размеры аномальной зоны не измеряются, а используется зависимость затухания сигнала от размеров зоны, определяется поправочный коэффициент регрессии r при расчете М. В интервале времен наблюдения [2001-2003 гг.] по апостериорным данным состоявшихся землетрясений в диапазоне М [5,0 7,3] баллов получена величина регрессии r=1,2.

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