способ оперативного прогноза землетрясений

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

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

где В_х1(0), В_у1(0), В _х2(0), В_у2(0) - значения автокорреляционных функций сигналов в нуле в плоскостях х, у геофонов соответственно первого и второго пунктов.

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

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

На настоящий момент выявлено множество долгосрочных признаков - предвестников зоны подготавливаемого землетрясения [см., например, в кн. Т.Рикитаке «Предсказание землетрясений», перев. с англ. - М.: Мир, 1979 г., табл.15.13, стр.314-333]. Известные статические признаки имеют продолжительные интервалы (годы) существования, но не позволяют достаточно точно предсказать момент удара и его магнитуду. Другой класс - краткосрочные динамические признаки-предвестники. Они проявляются за несколько суток (часов) до удара, но в силу отсутствия технических средств для их измерения не могут быть своевременно обнаружены. Среди наиболее значимых краткосрочных динамических признаков-предвестников, измеренных средствами глобальной навигационной системы позиционирования (GPS, NAVSTAR,) является раскачка очага землетрясения перед ударом. Раскачка сопровождается распространением от гипоцентра очага сверхдлинных (с периодом ~10⁴ сек) литосферных волн [см., например, «Способ предсказания землетрясений», Патент RU № 2170446, G.01.V, 9/00, 2001 г.].

Известен и такой признак-предвестник землетрясений, как изменение спектра шумов сейсмического фона, так называемое «затишье» непосредственно перед ударом [см., например, Патент RU № 2181205, G.01.V, 9/00, 2002 г.]. Достоверное предсказание могут обеспечить те методы регистрации, которые основаны на измерении первопричин землетрясения. В теоретическом плане существует несколько геофизических моделей, претендующих на обоснование первопричин землетрясений [см., например, «Краткосрочный прогноз катастрофических землетрясений с помощью радиофизических наземно-космических методов», Доклады конференции, РАН, ОИФЗ им. О.Ю.Шмидта. - М., 1998 г., стр.14-16].

Одна из перспективных моделей сейсмического процесса учитывает восходящую диффузию легких газов из земной коры: водорода, гелия. Эстафетный механизм передачи упругой энергии легкими газами снизу вверх приводит к установлению периодических зон растяжения и сжатия. При этом отдельные блоки зоны расширения сливаются, образуя крупномасштабную зону неустойчивого состояния - очаг будущего землетрясения. Факт эманации легких газов из земной коры в атмосферу накануне удара см., например, Патенты RU № 2204852, 2003 г., № 2275659, 2006 г., № 2302020, 2007 г.

Известен «Способ предсказания землетрясений» на основе регистрации и обработки сигнала сейсмического фона. Патент RU № 2130195, 1999 г. - аналог. В способе-аналоге создают на выбранном профиле наблюдений измерительный полигон в виде прямоугольной решетки из N-безынерционных измерителей, размещенных в узлах решетки и отстоящих друг от друга на расстояние /4 при геометрических размерах сторон решетки, кратных длине волны сигнала- предвестника, преобразуют сейсмические волны в электрический сигнал, измеряют амплитуду А сигнала в каждом узле решетки со скважностью менее 1 сек, формируют матрицу цифровых отсчетов сигнала- предвестника размерностью m×n элементов в функции пространственных координат А(х,у), вычисляют параметры электрического сигнала матрицы: энергетический спектр сигналов S(F_x), S(F_y), пространственный период направление прихода волны: =arctg F_x/F_y, автокорреляционную функцию B(R) и по их значениям и времени существования судят о магнитуде и времени ожидаемого землетрясения.

К недостаткам способа-аналога следует отнести:

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

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

В технике при диагностировании объектов широко используют математические процедуры обработки функций электрических сигналов. Наибольшее распространение получили математические методы разложения функций сигналов в дискретные либо непрерывные ряды. Известны разложения функций в степенные либо ортогональные ряды Маклорена, Тейлора, Фурье [см., например, Н.С.Пискунов, «Дифференциальное и интегральное исчисления для ВТУЗов, т.2, учебник, 5-е изд. М.: Наука, 1964 г., Ряды, Ряды Фурье, стр.180-182, стр.218-221]. При этом коэффициенты разложения функций в ряд могут иметь какой-либо физический или наглядный смысл.

Ближайшим аналогом к заявляемому способу, где используют разложение сигнала сейсмического фона в ряд Фурье, является «Способ краткосрочного предсказания землетрясений», патент RU № 2181205, 2005 г. В способе ближайшего аналога регистрируют волны сейсмического фона в виде непрерывной последовательности дискретных отсчетов амплитуды сигнала A(t) в двух, разнесенных по координатам пунктах, рассчитывают спектр Фурье от зарегистрированной функции с объемом отсчетов в каждой выборке N 2F_max/ , вычисляют интервал корреляции , регистрируют начало изменения параметра и при его непрерывном отслеживании фиксируют время запаздывания изменения фазы данного признака между двумя пунктами, рассчитывают направляющий косинус траверзы прихода сверхнизких волн очага , определяют гипотетический центр очага как точку пересечения на траверзе радиус-векторов с косинусом угла при вершине

Вычисляют период Т₀ параметра и по его величине прогнозируют магнитуду М 110/Т₀ ² (час) и время удара t_x 2,3Т₀, где:

F_max - максимальная частота спектра сейсмического фона, Гц;

- среднеквадратическая ошибка вычисления спектра Фурье по дискретной выборке измерений;

а - длина базы между двумя пунктами, м;

V - скорость сейсмических волн в земной коре, м/с;

B₁(0), В ₂(0) - значения автокорреляционных функций сигнала в нуле для каждого пункта.

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

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

- неточность регрессионных зависимостей определения магнитуды и времени удара, в частности известно из соотношений Гутенберга-Рихтера, что чем больше время (Т₀), тем магнитуда ожидаемого удара должна быть больше.

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

Технический результат достигается тем, что способ оперативного прогноза землетрясений, включающий регистрацию сейсмических волн в виде дискретных отсчетов амплитуд (x_к, у _к) сигналов во взаимно ортогональных плоскостях в двух, разнесенных на измерительной базе пунктах, причем оси чувствительности датчиков пунктов по координате x ориентируют по направлению базы, обработку зарегистрированных данных путем разложения в функциональный ряд и вычисления автокорреляционных функций сигналов, дополнительно регистрируют естественный акустический шум Земли в виде акустограмм посредством геофонов, располагаемых в глубоких скважинах, визуализируют точечное множество отсчетов последовательности данных в виде диаграмм Пуанкаре в плоскостях (x_к, x_к+1 ) (у_к, у_к+1), отслеживают динамику изменения формы диаграмм Пуанкаре и определяют постоянную времени переходного процесса (Т), прогнозируют время ожидаемого удара t_у 4,7T, магнитуду из соотношения lg t_у[сут]=0,54M-3,37, а координаты гипоцентра очага отождествляют с точкой пересечения направляющих косинусов, отсчитываемых от оси базы:

где B_x1(0), B_y1(0), B_x2(0), B_y2(0) - значения автокорреляционных функций сигналов в нуле в плоскостях х, у геофонов соответственно первого и второго пунктов.

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

фиг.1 - функция сигнала сейсмического фона (распечатка с ПЭВМ);

фиг.2 - решение дифференциального уравнения движения земной коры при ее насыщении газовой компонентой для математической модели;

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

фиг.4 - плотность распределения вероятностей признака - предвестника изменения спектра сейсмического фона;

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

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

Техническая сущность способа состоит в следующем.

Насыщение земной коры газовой компонентой приводит к изменению ее вязкоупругих характеристик. Математическая модель вязкоупругих характеристик земной коры представлялась в виде системы дифференциальных уравнений:

i=1 n

Функции k_i(t) для всех i определяются как

k_i(0)=0

р - вероятность сброса упругих связей.

В начальный момент t=0 элементы неподвижны х _i=0, i=0 n+1.

Масса m_i совершает колебание x₀(t)=sin( t), масса m_i+1 неподвижна x_n+1(t)=0.

Блоки земной коры в модели контактируют через систему последовательных элементов с упругой связью. Частота собственных колебаний механической системы и размеры колебательных «зерен» зависят от коэффициента жесткости k_i(t+ ·dt), который изменяется от степени насыщения земной коры газовой компонентой. Решение дифференциальных уравнений модели осуществлялось по специально разработанной математической программе. Из решения дифференциальных уравнений модели следует, что размеры «зерен» увеличиваются от единиц метров до нескольких км, а частота сейсмического фона очаговой зоны изменяется от звукового диапазона (единицы кГц) до инфразвукового диапазона (долей Гц). Изменение спектра сейсмического фона очаговой зоны (как одна из реализации решения дифференциальных уравнений модели) иллюстрируется графиком фиг.2.

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

За счет насыщения газом коэффициенты упругости в связях неограниченно возрастают во времени, стремясь к обеспечению абсолютно жесткого сцепления, а вся область зоны подготавливаемого землетрясения превращается в монолитный блок. При этом частота собственных колебаний блока из-за его больших размеров (100-150 км в диаметре) составляет 10^-3 10^-4 Гц.

На фиг.2 - это участок сверхдлинных литосферных волн раскачки очага, измеренный средствами космической навигационной системы GPS [см. Патент RU № 2170446, 2001 г.], занимает интервал ~8 час.

Если осуществлять непрерывный контроль сейсмического фона в звуковом и инфразвуковом диапазоне, то можно расширить интервал упреждающего предсказания землетрясений на время прохождения динамического процесса раскачки очага перед ударом через этот интервал. На фиг.2 - это начальный и средний участки графика функции с временем существования несколько суток. Известно [см., например, Р.Дуда, П.Харт «Распознавание образов и анализ сцен», перев. с англ., М.: Мир, 1976 г., гл. Пространственное дифференцирование], что максимум информации о происходящих процессах содержится в образе объекта.

Скрытая информация о динамике процесса на его ранней стадии может быть извлечена из электрического сигнала сейсмического шума, если математическими процедурами обработки визуализировать и отслеживать изменение формы сигнала (образа объекта). Для этого программным методом последовательность данных отображают в виде диаграмм Пуанкаре [см., например, Пуанкаре А. Избранные труды в 3-х т., пер. с французского под ред. Н.Н.Боголюбова. - М.: Наука, 1974 г.]. Диаграмма Пуанкаре представляет собой точечное графическое отображение N-значений последовательности, например Х_к при к=1, 2, 3 N на двумерном поле, в котором ординатой точки является значение Х_к+1, а абсциссой - предшествующее значение Х_к. Нанося поочередно точки для к=1, 2, 3 N на график, получают точечное множество Х_к+1 , (Х_к), образующее фигуру, по которой можно судить о типе последовательности (образе объекта). Текст математической программы построения диаграммы Пуанкаре представлен в примере реализации. На фиг.3 иллюстрируются диаграммы Пуанкаре для ортогональных каналов (х, у) чистого акустического шума Земли (а) и после предвестия (б). Как следует из полученных графиков, форма сигнала (образ объекта) существенно изменяется. Известно [см., например, «Теоретические основы радиолокации» под ред. В.Е.Дулевича. - М., Сов. радио, 1984 г., стр.114], что плотность распределения вероятностей огибающей чистого шума подчиняется закону Релея. Появление в спектре шума регулярной составляющей изменяет плотность распределения вероятностей W(u) амплитуды огибающей, которая подчиняется обобщенному закону Релея

где I₀(au) - функция Бесселя нулевого порядка от мнимого аргумента;

u=V/ - относительная величина результирующего сигнала;

a=А/ - относительная величина регулярной составляющей.

Изменение плотности распределения вероятностей амплитуды результирующей огибающей иллюстрируется графиком фиг.4. При больших амплитудах регулярной составляющей (a=А/ ) обобщенный закон Релея стремится к нормальному. Этот же результат следует из визуализации диаграмм Пуанкаре. Изменяются как форма признака-предвестника, так и длина вектора Х_к+1 , Y_к+1, фиг.3б. Благодаря визуализации признака-предвестника представляется возможным отследить всю динамику переходного процесса. Переходные процессы описываются дифференциальным уравнением первой степени, общим решением которых является экспонента [см., например, Н.С.Пискунов. Дифференциальное и интегральное исчисления для ВТУЗов, т.1, учебник, 5-е издание. М., Наука, 1964 г., стр.458]. Экспоненциальная зависимость обладает тем свойством, что по ее отдельному участку можно восстановить всю траекторию, вычисляя постоянную времени (Т) переходного процесса из соотношения:

где u₁, u₂, u ₃ - длина вектора признака-предвестника из диаграмм Пуанкаре;

t=(t₂-t₁); (t₃-t₂ ) - интервал времени наблюдения между отсчетами u₁ , u₂, u₃;

u₀ - установившееся значение длины признака-предвестника

Расчетные соотношения для определения (Т) иллюстрируются графиком фиг.4. По постоянной времени переходного процесса прогнозируют характеристики ожидаемого землетрясения.

Время удара - это интервал времени, за который длина признака-предвестника с вероятностью, близкой к единице, достигает установившегося значения u₀ для экспоненты: t _уст 4,7T (с вероятностью 0,99).

Магнитуду удара находят из соотношения Гутенберга-Рихтера:

lg t_у[сут]=0.54M-3,37.

Гипоцентр землетрясения находят методом пеленгации источника, генерирующего регулярную составляющую в акустическом шуме Земли, с двух, разнесенных на измерительной базе пунктов. Известно, что направление переноса энергии волновым процессом совпадает с фазовым фронтом волны в данной точке. Направление вектора в пространстве задается его проекциями на осях координат. Проекции вектора направления переноса энергии пропорциональны мощности шумов, регистрируемых геофонами во взаимно ортогональных плоскостях (х, у). В соответствии с операциями ближайшего аналога значение автокорреляционных функций в нуле В₀ равно мощности процесса, т.е. сумме мощностей постоянной (средней, m²) и переменной (дисперсии, ²) составляющих, для чего определяют средние значения и дисперсию шума в каждом канале х, у. Поскольку направление оси х геофонов ориентируют по направлению базы, то косинус направляющие на гипоцентр очага соответственно составляют:

Гипоцентр очага отождествляют с точкой пересечения векторов, построенных на карте местности под углами ( ', ) относительно направления измерительной базы. Определение гипоцентра очага методом пеленгации иллюстрируется фиг.6.

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

Функциональная схема устройства, реализующего способ, фиг.7, содержит оборудованный на соответствующем геолого-структурном образовании полигон 1 с пробуренными обсадными скважинами 2, разнесенными на измерительной базе 3. В каждую из скважин на твердое основание коренной породы 4 опущены геофоны 5, 6 (типа МАГ-3С), осуществляющие преобразование акустического шума Земли в электрический сигнал во взаимно ортогональных плоскостях X, Y. Выходы геофонов через измерительные усилители 7, 8 (модель 2635 фирмы Bruel & Kjair) подключены к тракту обработки 9, из последовательно подключенных многофункционального блока 10 (модель 3560-L, Bruel & Kjair) в составе канального коммутатора 11, спектроанализатора 12, аналогово-цифрового преобразователя 13, устройства ввода данных 14, компьютера 15 в стандартном наборе элементов: процессора 16, оперативного запоминающего устройства 17, винчестера 18, дисплея 19, принтера 20, клавиатуры 21.

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

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

program Pua;

uses crt, graph;

const kp=200;

type

Prec=^rec;

rec=record dat: real; next: Prec end;

var f: text;

maxx, minx, m, d, S, x, у: real;

n, i, xx, yy, GrDr, GrMod: integer;

str: string;

p, q: Prec;

begin

clrscr;

repeat

writeln ('Введите имя файла с данными');

readln (str);

assign (f, str);

i: = IOresult;

if i<>0 then writeln ('Нет такого файла')

until i=0;

S:=0;

p:=nil;

n:=0;

while not eof (f) do begin

readln (f, x);

if n=1 then begin maxx:=x; minx:=x end;

if x>maxx then maxx:=x;

if x<minx then minx:=x;

S:=S+x;

n:=n+1;

new (q);

q^.dat:=x;

q^.next:=p;

p:=q

end;

m:=S/n;

d:=0;

q:=p;

while q<>nil do begin

x:=q^.dat;

d:=d+sqr(x-m)/n;

q:=q^.next

end;

writeln ('Всего чисел', n:5);

writeln ('maxx=', maxx:15:7, 'minx=', minx:15:7);

writeln ('Среднее значение=', m:15:7);

writeln ('Средний квадрат отклонения от среднего =',d:15:7);

writeln ('Для построения диаграммы Пуанкаре нажмите Enter');

readln;

GrDr:=detect;

initgraph(GrDr,GrMod, '');

line (10, 10, 10, 10+kp);

line (10, 10+kp, 10+kp, 10+kp);

q:=p;

if q<>nil then begin x:=q^.dat; q:=q^.next end;

while q<>nil do begin

y:=q^.dat;

xx:=10+round((x-minx)/(maxx-minx)*kp);

yy:=10+kp-round((y-minx)/(maxx-minx)*kp);

putpixel (xx, yy, 255);

x:=y;

q:=q^.next

end;

repeat until keypressed;

closegraph

end.

Из графиков фиг.3б установившееся значение u₀ в шкале квантования 0 256 уровней составило величину ~250. Значения u₁ и u₂ в интервале наблюдения t=8 час (t₁=8 час, t₂=16 час) соответствовали u₁=125, u₂=175. Откуда постоянная времени процесса:

Ожидаемое время удара t_у 4,7T=70 час=2,9 суток

Ожидаемая магнитуда удара М=7,1 балла.

Направляющие косинусы гипоцентра очага землетрясения (фиг.3а) для первого пункта

по направляющему косинусу второго пункта =48°.

Гипоцентр очага землетрясения находился в Охотском море, где оно и произошло 18.11.2002 г., магнитудой 7,3 балла.

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