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

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

Способ прогнозирования параметров землетрясения, включающий получение изображений очага по двум взаимно ортогональным по поляризации каналам съемки, формирование синтезированной матрицы изображения из попиксельных амплитуд сигналов поляризационных каналов, обработку матрицы синтезированного изображения, отличающийся тем, что съемку средствами орбитального комплекса осуществляют в режиме с равными промежутками времени t между смежными витками, вычисляют математическое ожидание сигналов изображений поляризационных каналов, синтезированную матрицу формируют из попиксельных отношений амплитуд сигнала с большим математическим ожиданием к меньшему, последовательно от начала синтезированной матрицы разбивают изображения на фрагменты по четыре смежных пиксела, делят каждый четырехугольник диагональю на два смежных треугольника и по формуле Герона вычисляют их площади, находят площадь рельефа (S_p) поверхности изображения как сумму площадей мозаики треугольников всех фрагментов, по серии последовательных во времени снимков очага отслеживают динамику изменения функции D=S_p/S₀ шероховатости рельефа, по зависимости вычисляют постоянную экспоненты Т наблюдаемого процесса, рассчитывают ожидаемое время удара t_у4,6T и его магнитуду (М) из соотношения lgt_y[сут]=0,54М-3,37, где t₂-t₁=t - интервал времени между двумя последовательными съемками, ч; D_1, D₂ - вычисленные функции шероховатости рельефа в моменты t₁, t₂; D₀ - предел шероховатости площади рельефа, при котором происходит «вспаривание» очага; S₀ - площадь синтезированной матрицы изображения из |m×n| пикселей.

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

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

За последние годы выявлено несколько новых, ранее недоступных для наблюдения признаков - предвестников грядущих землетрясений на основе использования методов дистанционного зондирования поверхности Земли из космоса. Среди них такие признаки как: раскачка очага землетрясения накануне удара, отслеживаемая средствами космической системы Navstar (Патент RU № 2.170.446, 2001 г. - аналог), образование над эпицентральной областью очага в атмосфере вертикального электростатического поля, напряженностью в несколько кВ/м (Патент RU №2.205.432, 2003 г. - аналог). Известен "Способ обнаружения очагов землетрясения", Патент RU № 2.209.452, 2003 г. - аналог.

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

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

Ближайшим аналогом к заявляемому техническому решению является "Способ предсказания землетрясений", Патент RU № 2.213.359, 2003 г., G.01.V, 9/00, 2003 г.

В способе ближайшего аналога получают изображения подстилающей поверхности А₁(х,у); А₂(х,у) по двум взаимно-ортогональным по поляризации каналам приема, формируют синтезированную матрицу путем попиксельного совмещения получаемых изображений с амплитудой сигнала в каждой точке методами пространственного дифференцирования выделяют контуры на синтезированном изображении, проводят фрактальный и линиаментный анализ сигналов фрагментов изображения внутри выделенных контуров и формируют градиентное поле направлений линиаментов, идентифицируют очаг землетрясения на изображении по геометрии контура, фрактальной размерности и узору рисунка поля линиаментов, вычисляют средневзвешенную сумму азимутов линиаментов и отслеживают изменение этой функции по серии последовательных изображений идентифицированного очага, фиксируют момент появления периодических отклонений функции средневзвешенной суммы азимутов линиаментов, рассчитывают период (Т) и коэффициент (k) наклона касательных отслеживаемой функции, прогнозируют Т время удара отсчитывамое от начала периодических отклонений функции и магнитуду удара из соотношения lgt_y[cут]=0,54M-3,37, где u ₀ - предельная величина сдвиговой деформации, при которой происходит разрыв земной коры.

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

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

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

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

Поставленная задача решается тем, что в способе прогнозирования параметров землетрясения, включающем получение изображений очага по двум взаимно-ортогональным по поляризации каналам съемки, формирование синтезированной матрицы изображения из попиксельных амплитуд сигналов поляризационных каналов, обработку матрицы синтезированного изображения, дополнительно, съемку средствами орбитального комплекса осуществляют в режиме с равными промежутками времени t между смежными витками, вычисляют математическое ожидание сигналов изображений поляризационных каналов, синтезированную матрицу формируют из попиксельных отношений амплитуд сигнала с большим математическим ожиданием к меньшему, последовательно, от начала синтезированной матрицы разбивают изображение на фрагменты по четыре смежных пикселя, делят каждый четырехугольник диагональю на два смежных треугольника и по формуле Герона вычисляют их площади, находят площадь рельефа S_p поверхности изображения как сумму площадей мозаики треугольников всех фрагментов, по серии последовательных во времени снимков очага отслеживают динамику изменения функции - D=S_p/S₀, шероховатости рельефа, по зависимости вычисляют постоянную экспоненты Т наблюдаемого процесса, рассчитывают ожидаемое время удара: (t_y)4,6T и его магнитуду (М) из соотношения lgt_y[сут]=0,54М-3,37, где:

t₂-t₁=t, интервал времени между двумя последовательными съемками, [час];

D₁,D₂ - вычисленные функции шероховатости рельефа в моменты t₁t₂;

D₀ - предел шероховатости рельефа, при котором происходит "вспарывание" очага;

S₀ - площадь синтезированной матрицы изображения из |m×n| пикселей.

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

фиг.1 - изменение альбедо (А) поверхности и поляризуемости () светового потока от напряженности электростатического поля в атмосфере над очагом;

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

фиг.3 - функции изменения шероховатости рельефа D=S_p/S₀ при различных T;

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

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

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

Наличие вертикального электростатического поля над поверхностью очага вызывает поляризацию молекул водяного пара, содержащегося в атмосфере. Исходя из классического принципа взаимодействия электромагнитного поля с веществом [см., например, "Преломление света". Физический энциклопедический словарь, под редакцией А.М. Прохорова, Из-во Сов.энциклопед., М., 1983 г., стр. 168] электроны и атомы вещества под действием световой волны совершают вынужденные колебания. Коэффициент преломления вещества (атмосферы) для падающего светового потока зависит как от концентрации вторичных излучателей в веществе, так и от соотношения длин волн падающего светового потока и собственного излучения вторичных вибраторов.

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

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

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

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

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

Поскольку площадь рельефа в каждом конкретном случае зависит от площади наблюдаемого очага (размеров матрицы изображения из |m×n| пикселей), вводят относительную величину, характеризующую только "шероховатость" рельефа: D=S _p/S₀, где S_p - расчетная площадь рельефа, S₀ - площадь прямоугольника, равная произведению числа строк (m) на число столбцов (n) и на площадь одного пикселя изображения. Интервал изменения функции D достаточно широк, что предопределяет высокую чувствительность способа к введенному показателю.

Раскачка очага землетрясения (а вместе с ней интенсивное выделение ионизированных газов в атмосферу) продолжается до тех пор, пока динамический напор колебательной массы не превысит предел прочности земной коры, при которой наступает ее разрыв, скол, сопровождаемый ударом. Поскольку "вспарывание" очага землетрясения происходит при одних и тех же значениях предельной прочности земной коры, то существует предел (D₀), к которому стремится связанная с данным процессом функция шероховатости рельефа изображения. Динамика изменения функции площади рельефа содержит информацию о магнитуде и времени удара. Связь между динамикой процесса (скоростью) и самой функцией представляется дифференциальным уравнением первой степени. Известно [см., например, Н.С.Пискунов, "Дифференциальное и интегральное исчисления" для Втузов, учебник, 5-е изд., М., Наука, 1964 г., стр 443-451], что общим решением дифференциального уравнения первой степени является экспонента. Начальные условия для численного решения дифференциального уравнения определяют по двум последовательным во времени, снимкам отслеживаемого очага. Из свойств экспонента следует, см. фиг.3):

где t₂-t₁=t, интервал времен

и между двумя последовательными съемками, час;

D ₁D₂ - вычесленные функции шероховатости рельефа в моменты t₁t₂;

Т - постоянная экспоненты (начальные условия).

Величину D₀ вычисляют апостериори, по серии состоявшихся землетрясений.

Множество решений исходного дифференциального уравнения при различных значениях постоянной экспоненты иллюстрируется фиг.3. Постоянная экспоненты характеризует переходный колебательный процесс от момента начала раскачки очага до удара. Чем больше колебательная масса очага, тем постоянная времени процесса T больше. Существует формула Гутенберга-Рихтера, связывающая время существования признака-предвестника (t) с магнитудой ожидаемого землетрясения М [см., например, «Краткосрочный прогноз катастрофических землетрясений с помощью радиофизических, наземно-космических методов», Доклады конференции, ОИФЗ им. О.Ю.Шмидта, РАН, М., 1988 г., стр.10].

Для сейсмических волн эта зависимость имеет вид:

lgt[cут]=0,54M-3,37.

Время существования признака предвестника определяют из свойств экспоненты. С вероятностью 0,99 экспонента достигает предельной величины при t=4,6Т. Таким образом, время удара, отсчитываемое от момента нарушения «гладкости» рельефа, составит величину t_у=4,6T.

Саму постоянную T находят из вышеприведенной зависимости вычислением значений D₁ и D₂ шероховатости рельефа на моменты съемки t₁,t₂.

Всю площадь рельефа представляют мозаикой треугольников. Площадь каждого треугольника вычисляют по формуле Герона:

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

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

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

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

Заявляемый способ может быть реализован по схеме фиг.5. Функциональная схема устройства фиг.5 содержит орбитальную станцию 1, типа МКС, с установленными на ее борту цифровыми соосными, ортогональными по поляризации телекамерами 2,3, типа BYP-30 (Япония), принимающих отраженный световой поток от подстилающей поверхности.

Полученные в каналах приема изображения записываются на бортовой магнитофон 5 типа «Нива». Съемку запланированных сейсмических регионов осуществляют по программам или разовым командам, закладываемым в бортовой комплекс управления 6, передаваемым из центра управления (ЦУП) 7.

Отснятые (по целеуказаниям из ЦУПа) кадры изображений поверхности отслеживаемых очагов передают по оперативному каналу связи (8) на наземный пункт ретрансляции сообщений 9, где осуществляют запись на видеомагнитофон 10 типа «Арктур». Скомпонованные файлы информации (вместе со служебными признаками, время съемки, координаты) перегоняются в Геофизический Центр МЧС 11, где создается долговременный архив 12 всех снятых очагов землетрясений, на базе стримеров, типа FT-120. Обработку полученных изображений очагов, по операциям заявляемого способа осуществляют на ПЭВМ 13 в стандартном наборе элементов: процессора 14, оперативного ЗУ 15, винчестера 16, дисплея 17, принтера 18, клавиатуры 19. По результатам обработки снимков, выделения и отслеживания динамики изменения признака-предвестника создают базу эталонных данных 20, которую выводят на сервер 21 сети «Интернет».

Полученные в поляризационных каналах съемки 2, 3 изображения вводят в ПЭВМ 13, из которых формируют синтезированную матрицу. Вычисление синтезированной матрицы по известному алгоритму попиксельного отношения амплитуд представляется стандартной математической операцией, входящей в комплект специализированного программного обеспечения MATH САД 7.0 PLUS [см., например, MATH САД, издание второе стереотипное, Информиздат, дом. Филинъ, М., 1997 г., Векторизация элементов матрицы, стр.211]. Исходная синтезированная матрица (фиг.2) содержала 512 строк и 480 столбцов, полученная в каналах видеосъемки с разрешением х=у=120 м на пиксель. При этом площадь поверхности очага на местности равна

S₀=[512×120×480×120]356 км².

Амплитуда пикселей квантовалась в стандартной шкале 0...256 уровней.

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

Текст программы расчета площади рельефа.

Program SqrGeron;

Const maxx=1000;

Var S:real;

Ss:string;

St1, st2:array [1...maxx] of byte;

F:text;

i,j kolx:integer;

dx, dy, dh:real;

procedure SqrSqr (var S:real; x1, x2, y2:real;

h11,h12,h21,h22:byte);

var dS, dS1, dS2:real;

function Geron (x1, y1, h1, x2, y2, h2, x3, y3, h3:real):real;

var 112, 123, 131, p:real;

begin

112:=sqrt (sqr(x1-x2)+sqr(y1-y2)+sqr(h1-h2));

123:=sqrt (sqr(x2-x3)+sqr(y2-y3)+sqr(h2-h3));

131:=sqrt (sqr(x3-x1)+sqr(y3-y1)+sqr(h3-h1));

p:=(112+123+131)/2;

Geron:=sqrt (p*(p-112)*(p-123)*(p-131));

end;

begin

dS1:=0;

if (h11<>255) and (h12<>255) and (h21<>255) then

dS1:=dS1+Geron (x1, y1, h11*dh, x2, y1, h12*dh, x1, y2, h21*dh);

if(h22<>255) and (h12<>255) and (h21<>255) then

dS1:=dS1+Geron (x2, y2, h22*dh, x2, y1, h12*dh, x1, y2, h21*dh);

dS2:=0;

if (h11<>255) and (h12<>255) and (h22<>255) then

dS2:=dS2+Geron (x1, y1, h11*dh, x2, y1, h12*dh, x2, y2, h22*dh);

if(h22<>255) and (h11<>255) and (h21<>255) then

dS2:=dS2+Geron (x2, y2, h22*dh, x1, y1, h11*dh, x1, y2, h21*dh);

if dS1>0 then dS:=dS1;

if(dS2>0) and (dS2<dS) then dS:=dS2;

S:=S+dS

end;

begin {main}

S:=0;

repeat

writeln (Enter file name, please');

readln (ss);

assign (f, ss);

{SI-} reset (f);{SI+}

i:=IOresult;

until i=0;

writeln ('Enter dx, dy,dh');

readln (dx, dy, dh);

kolx:=l;

if not eof (f) then

while not eoln (f) or (kolx>maxx) do begin

read (f, st1[kolx]); kolx:=kolx+1;

end;

if eoln(f) then begin;

readln (f);

j:=1;

while not eof (f) do begin

for i:=1 to kolx do read (f, st2[i]);

for i:=1 to kolx-1 do

SqrSqr (S, 0, dx, 0, dy, st1[i], st1[i+1], st2[i], st2[i+1]);

readln (f);

St1:=st2;

j:=j+1;

end;

end else writeln ('Very short array');

close (f);

assign (f, 'result.txt');

rewrite (f);

writeln (f, 'file name -', ss);

writeln (f, dx=', dx:10:5, 'dy:10:5, 'dh=',dh:10:5);

writeln (f, 'Strok -'j:5, 'Stolb -', kolx);

writeln (f, 'Area of region is', S:10:2);

writeln ('file name -', ss);

writeln ('dx=', dx:10:5, 'dy=', dy:10:5, 'dh=', dh:10:5);

writeln ('Strok -', j:5, 'Stolb -', kolx);

writeln ('Area of region is', S:10:2);

close (f)

end. {main}

Съемка осуществлялась в зоне видимости МКС с наземного пункта в начале зоны (на первом витке) и в конце зоны (на третьем витке). Интервал времени между двумя последовательными съемками составил: 1,5 часа на виток × 2=3 часа. Отношение в полученных снимках составило 2,2. Откуда постоянная экспоненты T=3/ln2,2=4,25 часов.

Ожидаемое время удара t_y =4,25×4,6=19,5=0,82 суток. Ожидаемая магнитуда удара lg0,82=0,54М-3,37. М=6,1 балла. Большой диапазон изменения признака-предвестника в виде площади рельефа поверхности синтезированной матрицы (порядка 12) и репрезентативность выборки измерений [(объем m×n) отсчетов] обеспечивают высокую чувствительность и статистическую устойчивость рассмотренного способа прогнозирования.