способ предсказания землетрясений

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

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



а магнитуду землетрясения М из соотношения

lg t [сут] = 0,54 М-3,37,

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

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

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

Известно множество косвенных признаков возможного землетрясения, регистрируемых различными физическими методами. Систематизированный перечень известных признаков см., например, Т. Рикитаке "Предсказание землетрясений", перевод с англ., М., Мир, 1979 г., стр.314, таблица 15.13.

При наличии надежных предвестников существует линейная зависимость между логарифмом времени предвестника t [сут] и магнитудой ожидаемого землетрясения М, формула Гутенберга-Рихтера:

lgt=0,79 М-1,88=aM-b.

(см. , например, Т. Рикитаке "Предсказание землетрясений", перевод с англ.. Мир, М., 1979 г., стр.242 - аналог).

В последние два десятилетия проводились обобщения данных об аномальных полях, предшествующих землетрясениям, в виде регрессионной зависимости Гутенберга-Рихтера. Для различных предвестников (более 1000 наблюдений) получены следующие зависимости:

lgt=0,77 М-4,4 (деформационное поле)

lgt=0,54 М-3,37 (сейсмические волны)

lgt=0,3 М-1,84 (электросопротивление)

lgt=0,18 М-2,5 (земные токи)

lgt=0,24 М-2,47 (дебит флюидов)

(см. , например, "Краткосрочный прогноз катастрофических землетрясений с помощью радиофизических наземно-космических методов", Доклады конференции, ОИФЗ им. О.Ю. Шмидта, РАН, М., 1998 г., стр.10-12 - аналог).

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

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

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

- большой разброс значений регрессионных коэффициентов а, b, снижающих точность результирующего прогноза.

Наиболее достоверным динамическим признаком-предвестником землетрясения является раскачка очага перед ударом, сопровождаемая распространением от него сверхнизких литосферных волн с периодами от 3,2 до 4,5 час (см., например, Давыдов В. Ф. "Землетрясения. Телеметрия предвестников". Монография, МГУЛ, М., 2001 г., стр.19, рис.9).

В качестве устойчивого признака-предвестника землетрясения при этом рассматривают период сверхнизких литосферных волн. Однако измерение медленно меняющихся процессов с периодами <10 сек является трудно выполнимой технической задачей (см., например, "Геофизические методы мониторинга природных сред" под редакцией В.Н. Сорокина, АН СССР, Институт общей физики, М., 1991 г., стр.267).

В известном способе измерений (см., например, патент РФ 2170446, G 01 V 9/00, 2001 г. - ближайший аналог) подобная задача решается методом выноса измерителя (наблюдателя) за пределы инерциальной системы и отслеживания колебательного процесса как бы со стороны путем использования эфемеридной информации космической навигационной системы типа Navstar.

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









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



и магнитуду



где Т - период отклонения координат (час), dekr - натуральный логарифм отношения амплитуд отклонения координат двух смежных периодов; (d=15), (l= 20) - коэффициенты регрессии.

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

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

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

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

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

а магнитуду землетрясения определяют из соотношения:

lgt[сут]=0,54 М-3,37.

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

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

фиг. 1 - функция электростатического поля, возникающая в атмосфере над очагом накануне удара;

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

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

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

- высокая достоверность, связанная с прямым измерением признака-предвестника в виде электростатического поля над очагом, появляющегося накануне удара;

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

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

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



где с - жесткость сжатой породы, m - колебательная масса очага, - сопротивление среды, пропорциональное первой степени скорости, Q(р) - вынуждающая сила, p - частота вынуждающей силы (см., например, Н.С. Пискунов "Дифференциальное и интегральное исчисление", учебник для втузов, 5-е изд., 1964. Стр. 526-532).

Решение этого неоднородного дифференциального уравнения состоит из двух слагаемых.

Общего решения однородного уравнения, соответствующего данному неоднородному



Частного решения неоднородного уравнения (для случая, когда частота вынуждающей силы совпадает с частотой свободных колебаний)



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



а амплитуда нарастает во времени. Касательные амплитуды определяются суммой функций



Разложив функцию e^-ht в ряд Маклорена для точки начала координат и ограничившись первыми двумя членами разложения, получено, что крутизна касательных представляется сложной зависимостью как периода (Т), так и магнитуды (М) ожидаемого удара. Следовательно, для достоверного прогнозирования характеристик землетрясения по параметрам колебательного процесса раскачки очага следует вычислять как период (Т), так и крутизну касательных (k).

Раскачка очага накануне удара сопровождается интенсивной эманацией радона в атмосферу. Непрерывное поступление радона приводит к накоплению нескомпенсированного заряда ионов в атмосфере и возникновению сильного электростатического поля над очагом. Вид функции электростатического поля иллюстрируется на фиг.1 (см., например, "Краткосрочный прогноз катастрофических землетрясений с помощью радиофизических наземно-космических методов", Доклады конференции, РАН, ОИФЗ им. О.Ю. Шмидта, М., 1998 г., стр.27). Вследствие раскачки имеет место параметрическая модуляция функции электростатического потенциала колебаниями очага. Поскольку малейшие подвижки земной коры сопровождаются эманацией газов (радона и др.) в атмосферу, то при наличии высокочувствительных измерителей электростатического поля представляется возможным зарегистрировать начало процесса без наличия "мертвой "зоны измерений, присущей способу-прототипу.

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

В соответствии с вышеизложенным решением дифференциального уравнения функцию, характеризующую раскачку очага землетрясения, представим в виде (фиг.2):

x=A(t)sint;

где A(t)=kt - нарастающая во времени амплитуда раскачки,

k - крутизна касательных.

Раскачка продолжается до тех пор, пока динамический напор колебательной массы очага не превысит предел прочности земной коры. Земная кора терпит разрыв при относительных напряжениях >10^-4 (см., например, "Теория предсказания землетрясений" в книге Т. Рикитаке "Предсказание землетрясений" перевод с англ.. Мир, М., 1979 г., стр.306-309). Поскольку земная кора очага уже находится в состоянии механических напряжений, то незначительная добавка динамического напора колебательной массы вызывает разрыв породы, проявляющийся в виде толчка землетрясения. Динамический напор максимален, когда максимальна скорость движения колебательной массы. Скорость - это первая производная модулирующей функции:

x"=(ktsint)";

Максимум скорости соответствует моменту времени, когда "синусоида" колебаний пересекает ось абсцисс или



Если x"mах по истечении времени t=t_y превысит предельно допустимую величину u_o, то наступит удар. Откуда время удара рассчитывают из соотношения



Таким образом, рассчитав по дискретной выборке измерений период Т и наклон касательных (k) модулирующей функции, вычисляют ожидаемое время удара t_у. Количественная оценка расчетного времени удара приводится в примере реализации. Отождествляя время удара t_y со временем существования признака-предвестника t= t_у магнитуду предстоящего землетрясения рассчитывают по регрессионной зависимости представленного выше способа аналога для сейсмических волн:

lgt, сут=0,54М-3,37

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

Заявляемый способ может быть реализован по схеме фиг.3. Функциональная схема устройства, фиг.3, содержит два параллельных канала измерений 1, разнесенных на измерительной базе 2 в составе последовательно соединенных функционального преобразователя-политрона 3 и интегратора 4. На входы измерительных каналов 1 подключены электростатические датчики 5, а выходы каналов нагружены на дифференциальную мостовую схему 6, выход которой подключен к последовательной цепочке элементов из пороговой схемы 7, аналогово-цифрового преобразователя (АЦП) 8, буфера-накопителя 9 и ПЭВМ 10 в стандартном наборе элементов: процессора 11, оперативного запоминающего устройства 12, винчестера 13, дисплея 14, принтера 15, клавиатуры 16. Процессор 11 ПЭВМ 10 подключен к программируемой схеме выборки измерений 17, синхронизующей работу пороговой схемы 7, АЦП 8 и буфера-накопителя 9. Измерительные каналы размещают в зоне подготавливаемого землетрясения. Осуществляют балансировку дифференциальной мостовой схемы 6 таким образом, чтобы ее выходное напряжение в отсутствии сигнала предвестника было близким к нулю и не превышало напряжения установленного порога пороговой схемы 7. При возникновении сильного электростатического поля в атмосфере над очагом накануне удара его величина в разнесенных по пространству измерителях будет существенно отличаться. В результате усиления в политроне 3 и интегрирования 4 сигналов различного уровня, снимаемых с электростатических датчиков 5, баланс дифференциальной мостовой схемы нарушается. Выходное напряжение схемы 6 превысит пороговое напряжение схемы 7, задаваемое программой от программируемой схемы выборки измерений 17. Программа, предварительно сформированная на ПЭВМ 10 и заложенная в программируемую схему выборки 17, запускает АЦП 8, осуществляющий квантование амплитуды поступающего разностного сигнала со схемы 6 и его дискретизацию во времени. Поток цифровых данных с выхода АЦП 8 заполняет буфер-накопитель 9 и по достижении заданного программой объема файла пересылается в ОЗУ 12. Шкала квантования сигнала по амплитуде, длительность одного импульса, длительность цикла измерений, объемы кадра и файла определяются программой. Имеется возможность адаптации к измерительному процессу путем изменения формируемой на ПЭВМ программы.

Атмосферное электрическое поле как часть глобальной электрической цепи Земля - ионосфера может меняться на поверхности Земли в пределах от 50 до 200 В/м в зависимости от географических координат. В присутствии источника ионизации, каковым является эманация радона в атмосферу в сейсмоактивной области, величина аномального электростатического поля возрастает до значений 1кВ/м. Изменение электростатического потенциала над очагом иллюстрируется фиг.1. Селектируемыми параметрами сигнала признака-предвестника являются: момент возникновения, период и скорость нарастания амплитуды сигнала во времени. Скрытая информация о параметрах предстоящего удара содержится в модулирующей функции сигнала предвестника, изменяющейся синфазно с раскачкой очага. Выделение модулирующей функции сигнала по дискретным выборкам измерений осуществляют программным расчетом на начальном участке измерений. Для чего в непрерывно получаемых выборках измерений по разработанному алгоритму находят смежные значения максимума и минимума. В окрестностях этих точек производная анализируемой функции должна равняться нулю. После достоверного вычисления точек "mах" и "min" регистрируемой функции определяют интервал времени между этими точками функции. При известном дискретном времени каждого отсчета, задаваемого программой работы АЦП, определяемый интервал времени между "mах" и "min" регистрируемой функции равен Т/2.

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

Подпрограмма расчета периода функции и наклона касательных:

Выборка данных находится в массиве A:mass=array[1..N]of real.

Выборка состоит из n записей (n<=N), сделанных через равные промежутки времени dt.

Процедура возвращает значение периода Т, а также тангенсы углов наклонов огибающих максимумы tgx и минимумы tgn. Если получено нулевое значение, то для реализации алгоритма определения периода и (или) наклона огибающих недостаточно данных.

Procedure FindExtr(A:mass; n:integer; dt:real; var T,tgx,tgn:real);

Function MaxA(x,y,z:real):Boolean;

Begin

MaxA:=false;

If(x<=y) and (z<y) or (x<y) and (z<=y) then MaxA:=true

End;

Function MinA(x,y,z:real):Boolean;

Begin

MinA:=false;

If(x>=y) and (z>y) or (x>y) and (z>=y) then MinA:=true

End;

Var I, imax, imin, iox, ion:integer; Amax, Amin, Aox, Aon: real;

Begin

T:=0; tgx:=0; tgn:=0;

imin:=l; imax:=1; Amin:=0; Amax:=0;

iox:=l; ion:=1; Aox:=0; Aon:=0;

ForI:=2to n-1 do

Begin

If minA(A[I-1],A[I],A[I+1]) then

begin

ion:=imin; Aon:=Amin; imin:=I; Amin:=A[I]

end;

If maxA(A[I-1],A[I],A[I+1]) then

begin

iox:=imax; Aox:=Amax; imax:=I; Amax:=A[I]

end;

End;

If (imin>1) and (imax>1) then

Begin

T:=abs(imax-imin)*dt;

If iox>1 then tgx:=(Amax-Aox)/((imax-iox)*dt);

If ion>1 then tgn:=(Amin-Aon)/((imin-ion)*dt);

End;

End.

Достоверность и статистическая устойчивость вычисления периода модулирующей функции могут быть получены другими программными методами (см., например, "Справочник по MathCad.", использование функции предсказания Predict (стр.242)).

Из способа ближайшего аналога (рис.3) следует, что момент удара характеризовался следующими параметрами:

- длина измерительной базы l=39.906.453 мм;

- максимальная амплитуда отклонения l=140 мм;

- период колебаний в среде с сопротивлением T₁=4 час=0,166 сут.

Приравняв колебательные скорости ближайшего аналога и рассматриваемого примера, получим



где



или

.

Поскольку зависимость времени существования предвестника-аналога имеет размерность [сут] , а вычисленная скорость [м/сек], то полученное выражение следует домножить на число секунд в одних сутках. Тогда



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



где =Аmax-Аmin.

Ожидаемая магнитуда удара рассматриваемого примера реализации:

lg 0,93=0,54 М-3,37; М=5,9 баллов.

Элементы измерителя выполнены по стандартным электронным схемам и могут быть реализованы на существующей элементной базе. В качестве ПЭВМ используется IBM PC/AT 486/487. Программируемая схема выборки, аналогово-цифровой преобразователь, буфер-накопитель выполнены на стандартных интегральных платах, совместимых с контроллерами IBM PC/AT. Программируемая схемы выборки выполнена на плате ЛA-TMS-31, АЦП и буфер-накопитель выполнены на плате ЛА-20 (см., например, Якубовский Б. и др. "Цифровые и аналоговые интегральные микросхемы", Справочник, Радио и связь, М., 1990).

Тактовая частота использованных интегральных микросхем 2,5 МГц, что при 8 разрядном, стандартном (0. ..256 уровней) шаге квантования обеспечивает время отсчета одного измерения не более 0,3 ms. Объем измерений, необходимый для восстановления модулирующей функции, составляет порядка 310⁵ дискретных отсчетов.

В измерителе использованы стандартные дифференциальная мостовая и пороговая (ключевая) схемы (см., например, "Справочник по радиоэлектронике", том 2 под редакцией А.А. Куликовского, Энергия, М., 1968 г., стр.484, рис.19-26 - дифференциальная схема, а также "Справочник по радиоэлектронным устройствам", том 1, справочник под редакцией А.А. Куликовского, Энергия, М., 1978 г., 4.3 Электронные ключи, стр.339-346).

Эффективность способа определяется такими показателями, как достоверность, точность, оперативность. Чтобы не пропустить ожидаемое событие, система должна работать круглосуточно в дежурном режиме. Достоверность определяется пороговым напряжением программы измерений, формируемой на ПЭВМ, а точность и оперативность - частотой дискретизации и шагом квантования амплитуды измеряемого сигнала. Характеристики используемых интегральных микросхем не накладывают ограничений на точность оценки измеряемых параметров регистрируемого процесса.