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

Вычислительные технологии

Том 7, № 5, 2002

МЕТОДИКА ОЦЕНИВАНИЯ НЕКОТОРЫХ ПАРАМЕТРОВ ИМПУЛЬСНОГО ЭЛЕКТРОМАГНИТНОГО ИЗЛУЧЕНИЯ ЗЕМЛИ

Г. М. ВодинчАР Камчатский государственный педагогический университет Петропавловск-Камчатский, Россия e-mail: [email protected]

The article analyses the strategy of point and interval estimation of mathematical expectation signal parameters, which characterize electromagnetic pulsed terrestrial emission. A non-stationary Poisson process is used as the signal model. The results of testing of the strategy are given.

Введение

В последнее время оформилось перспективное направление работ по прогнозированию землетрясений выделением из различных геофизических полей сигналов, в спектре которых проявляются частоты волн приливного воздействия. К таким сигналам относятся, например, высокочастотный сейсмический шум [1-3], импульсное электромагнитное излучение (ИЭМИ) из литосферы [4, 5]. В этих сигналах непосредственно перед землетрясением наблюдаются стабилизации фаз гармонических компонент с частотами приливных волн O1 и/или M2 [6] и их высших гармоник. Поэтому важной является задача надежно определять временные промежутки, когда над фоном шумов появляется сигнал с частотами волн O1 и/или M2 и определенными значениями начальных фаз. Эту задачу естественно решать в скользящем вдоль ряда данных временном окне. При этом речь идет о построении не только точечных, но и интервальных оценок, позволяющих определять допускаемые погрешности.

В настоящей работе рассматривается решение данной задачи для ИЭМИ Земли. В качестве характеристики ИЭМИ используется интенсивность — число превышающих шумовой порог импульсов в минуту. Частотный состав переменного во времени математического ожидания числа импульсов рассмотрен в [5], где установлено, что в спектре проявляются гармоники, модулированные приливными волнами суточной и полусуточной групп [6]. Построение интервальных оценок для случайных нестационарных сигналов рассмотрено, например, в [7-9], где оценивание, как правило, ведется по нескольким реализациям или предполагается, что сигнал является композицией детерминированного сигнала и стационарного шума. Для ИЭМИ оценки на основе подобной модели построены в [10]. В настоящей работе используется другая модель сигнала, более точно соответствующая его структуре.

© Г. М. Водинчар, 2002.

1. Модель сигнала

Будем считать, что целочисленное время n изменяется от — N до N и число импульсов, превышающих шумовой порог в промежуток [—N + n; — N + n +1], равно xn. Превышение ИЭМИ шумового порога образует пуассоновский поток событий, значит, xn распределено по закону Пуассона, и различные элементы последовательности xn попарно независимы. Характеризующий распределение Пуассона параметр An, являющийся математическим

ожиданием для xn, меняется со временем. Так как частотный состав An известен [5], счи-

p

таем, что An = A0 + (Ai cos Win + Bi sinWin) с известными частотами w^ Требуется по

i= 1

одной реализации сигнала xn оценить параметры Ai и Bi.

2. Построение точечных оценок

Представим xn в виде суммы детерминированного сигнала An и шума

p

xn = Ao + ^ (Ai cos Win + Bi sin Win) + Sn,

i=1

где sn xn An.

Ясно, что Msn = 0 и Dsn = An. Для получения точечных оценок применим метод наименьших квадратов. Этот метод дает хорошие статистические оценки тогда, когда шумовые компоненты являются одинаково распределенными независимыми величинами [11]. В нашем случае шумовые компоненты по-разному распределены, поэтому свойства МНК-оценок нуждаются в изучении.

N

Построим МНК-оценки. Ясно, что sin win = 0. Легко установить непосредственным

n=-N

£ sin (N + 0.5) Wi и

суммированием, что > cos Win = -. Для упрощения вычислений введем

n=-N sin0.5wi

б sin (N + 0.5) Wi в . Т N n N fía

обозначения: «in = cos Win — . —, Pin = sin Win. Тогда «in = 0 и Pin = 0.

(2N+1)sin0.5Wi n=-N n=-N

В новых обозначениях An = + Yl (Ai«in + Bi^in), где = Ao + ^ А+ 0.5) Wi. Обоз-

i=i i=i (2N +1)sin Wi

NN

начим далее Y] ainа^и втв?«. через (ai,aJ-)и (^i, в?) соответственно. Непосредст-

n=-N n=-N

венным вычислением можно установить следующие выражения для (ai,aJ-) и (в^в?):

( )= N 1 sin (2N +1) Wi — sin2 (N + 0.5) Wi («i,ai) N + 2+ 2 sin Wi (2N + 1) sin2 0.5Wi'

. . sin (N + 0.5) (Wi — Wj) sin (N + 0.5) (wí + w?) sin (N + 0.5) Wi sin (N + 0.5) w?

(a' a •) ==-+-—-

2 sin0.5 (wi — w?) 2 sin0.5 (wi + w?) (2N +1)sin0.5wi sin0.5w?

пРи i = j;

(e,,ei) = N +2 — üni^ülWi;

2 2 sin Wi

ол sin (N + 0.5) (Wi — Wj) sin (N + 0.5) (wi + w?) . . .

(ei,e?) = o • пк f -\--o • nKf \-Г1 при i = j.

2 sin 0.5 (wi — w?) 2 sin 0.5 (wi + w?)

Видно, что (аг, а) и (вг, вз) бесконечно велики при N ^ то, когда г = ^, и ограничены, когда г = .

Точечные МНК-оценки С *, А*, В* параметров С, соответственно являются ре-

шениями трех независимых систем уравнений:

N

^ +1) С* = Хп, (1)

n=-N

Р

(a^aj) A* = (a^,x) i = 1, 2,...,p, (2)

j=i

p

) B* = (ßi,x) i =1, 2,...,p. (3)

j=i

Из доказываемого ниже предложения 1 следует обратимость основных матриц систем (2) и (3) при достаточной длительности сигнала xn. Пусть матрицы U = Ци^ ||p и V =

||Vij ||pxp являются обратными для матриц || (aj, aj) ||pxp и || (ßj, ßj) ||pxp соответственно. Тогда

pp

A* = Uij (aj, x) и B* = Vj (ßj, x). Определим также отклонения полученных оценок j=i j=i

p

от истинных значений. Умножив почленно равенство xn = + ^2 (Ajajn + Bjßjn) + sn на

j=i

p

ain и просуммировав по n от — N до N, получим (aj, x) = ^2 (аг, aj) Aj + (аг, s). С учетом

j=i

выражения для (aj,x) из системы (2) величины Aj — A* являются решениями системы

p

^(aj,aj) (Aj — A*) = — (aj,s), i = 1, 2,...,p. j=i

pp

Тогда Aj — A* = — Uij (aj, s). Аналогично можно установить, что Вг — B* = — ^2 Vj (ßj, s).

j=i j=i 1 n

Ясно также, что C — C* =--У^ sn.

' 2N +1 n=-N n

3. Статистические свойства и распределения точечных оценок

Из выражений для Aj — A*, Bj — B*, C—C* очевидно, что полученные МНК-оценки являются несмещенными. Покажем их состоятельность, предварительно доказав вспомогательное Предложение 1. Если у симметрической матрицы T (N) = ||tjj (N)||pxp внедиаго-нальные элементы ограничены и lim tu (N) = то, то при достаточно больших N ма-

N ^те

трица обратима и все элементы обратной матрицы бесконечно малы при N ^ то.

Доказательство. Все собственные значения (N) матрицы T (N) в комплексной плоскости лежат в области Гершгорина [12], являющейся объединением кругов

|z — t„ (N )|<^|tij (N )|, i =1, 2,...

p.

j=i 3=i

При N ^ то центры этих кругов также стремятся к бесконечности, а их радиусы ограничены. Тогда расстояние от области Гершгорина до начала координат неограниченно возрастает и lim (N) = то. Это гарантирует, что все собственные значения отличны от 0 при

N ^те

достаточно больших N и матрица обратима. Матрица T (N) — симметрическая, значит, обладает ортонормированной системой собственных вектор-столбцов u. Тогда обратная матрица представима в виде

T-i (N)

1 1 -ui----7777 Up

(N) 1 0p (N)

[ ui- ■ -Up ]T.

Все координаты векторов u не превосходят по модулю 1, а lim —-—- = 0. Тогда все

N^те (N)

элементы T-i (N) при N ^ то бесконечно малы.

Предложение 2. Оценки Л*, B*, C* являются состоятельными. Доказательство. Найдем дисперсию Л*, учитывая, что полигармоническая последовательность Ап ограничена. Значит, существует положительное число А такое, что А > Ап

(р N \ N/P \2 N/P 4 2

Y, OjÄ = Y, ajn An < А Y,

j=i n=-N I n=-N \ j=i / n=-N\7=i

N p p

А Z) Ё) uijUifcöjnttfcn = А У) ujjMjfc (aj, afc) = Аий.

n=-Nj,k=i j, fc= i

По предложению 1 u^ ^ 0 при N ^ то, тогда lim DA* = 0 и A* — состоятельная

N ^те

оценка. Аналогично доказывается и состоятельность оценки B*. Из системы (1) получаем, что DC* = C/(2N + 1), значит, оценка C* также состоятельна.

Итак, установлено, что МНК-оценки являются несмещенными и состоятельными, это дает возможность дальнейшего построения на их основе интервальных оценок амплитуд и начальных фаз. С этой целью выясним характер распределения Л* и B* при N ^ то. Знание распределения C* не понадобится. Л* и B* являются линейными комбинациями случайных величин (aj, x) и (в, x) соответственно, для которых справедливо

Предложение 3. Распределения величин (aj ,x) и (в, x) при N ^ то асимптотически нормальны.

Доказательство. Непосредственным вычислением с учетом того, что третий центральный момент распределения Пуассона совпадает с математическим ожиданием [11], можно установить следующее:

N

У M (ajnXn - M (ajnXn))3 - - Aj N,

i * N—^nn IX

n=-N N

V M (впХп - M (jXn))3 - BjN.

*—' N ^те

n=-N

NN

> Б (о^ж га ) - АсЖ, > Б (в^гаХга) ~ АоЖ Тогда последовательности и (Д^Хп} удовлетворяют условиям Ляпунова [7] и по

N

центральной предельной теореме при N ^ то случайные величины (а, х) = ^ а?гахга и N

(в?, х) = ^ распределены асимптотически нормально.

При построении оценок используются сигналы длительностью 40 000 отсчетов и более [4, 5], поэтому можно считать величины (аз-, х) и (вз, х) распределенными нормально. Будем считать также нормально распределенными их линейные комбинации А* и В*.

4. Построение доверительных интервалов для амплитуд и начальных фаз

С вероятностью 7 нормально распределенные случайные величины А* и В* лежат соответственно в интервалах [А* - г7Л/ Б А*; А* + г7Л/БА*] и [В* - г7Л/ Бв*; Вг* + г7Л/ Б В*], где г7 определяется с помощью функции Лапласа Ф(г) из условия 2Ф(г7) = 7 [11]. Большое число используемых отсчетов позволяет приближенно использовать эти интервалы в качестве доверительных для Аг и Вг, заменив дисперсии А* и В* их несмещенными и состоятельными оценками.

N / р \2

Дисперсия А* имеет вид Е Е иу азп Ап. Заменив в этом выражении множитель

п=-N у=1

Р N / р \ 2

Ап на АП = С * + Е (А*+ В*в^п), получим для нее оценку = Е Е «уа^п А^

з=1 п=-N у=1 у

N /р \2

Аналогично в качестве оценки для Б В* возьмем = Е Е А^. Покажем,

1 п=-N у=1 у

что данные оценки обладают нужными свойствами.

Предложение 4. А^ является несмещенной и состоятельной оценкой Ап.

Доказательство. Несмещенность следует из линейных свойств математического ожи-

р

Е (азпесу (С *, А*) +

3=1

дания. Рассмотрим |БАП|

p p

Б ( C* + Е A*j + Е В* j j=i j=i

p

ej„cov (c*,b*)) + e a-näncov (a*,b*) j,fc=i

p _ _

< E(k-nlVDC* Б A* + DC* БВ*) + j=i

E Б A* БВ*. Так как и ^jn ограничены при N ^ то, а дисперсии оценок

j,fc=i

1j, В* бесконечно малы, то БАП

* i '

Предложение 5. является несмещенной и состоятельной оценкой Б А* Доказательство. Несмещенность непосредственно следует из предложения 4. Анало-

гично предыдущему доказательству . | < £ ( Е Uj a^n) (£ л/бапба*

2

N I p \ / p

■ra"'^'

ra,m=-N \ j=1 1 \j=1

Так как БАП ^ 0 при N ^ то, то она ограничена сверху некоторым числом D. Тогда n/P \2/p \2 N/P \2 N/P 4 2

|DSAi | < DE Е

uij ajn I I / v uij ajm = DEE E E Uij a jm

n,m=-N \j=1 у \j=1 у n=-N \j=1 у m=-N Vj'=1

Аналогично доказательству предложения 2 последнее выражение равно Du2^ Значит, дисперсия SA. бесконечно мала, и оценка состоятельна.

Аналогично предложению 5 можно показать несмещенность и состоятельность оценки .

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

на рисунке. В плоскости переменных A и —В выстраивается доверительный прямоугольник для точки Тогда угол, под которым этот прямоугольник виден из начала координат, будет доверительным интервалом для Наименьшее и наибольшее расстояния от точек прямоугольника до начала координат будут соответственно нижней и верхней границами доверительного интервала для М^. Следует отметить, что при покрытии истинной точки (А^ ;В) доверительным прямоугольником получаются верные интервальные оценки и для амплитуд, и для фаз. Но если доверительный прямоугольник не покрывает точку, возможно, что интервальная оценка амплитуды или фазы окажется верной.

5. Тестирование методики

Описанная выше методика интервального оценивания была программно реализована в системе Borland C+—Ъ Тестирование алгоритма и программы проводилось на основе последовательности из 1 млн независимых пуассоновских величин, параметр An которых изменялся по известному полигармоническому закону из четырех гармоник. В качестве периодов гармоник использовались периоды приливных волн O1,S1,M2,S2 [6] в минутах. Тест-последовательность была получена на основе последовательности псевдослучайных чисел, сгенерированной стандартной библиотечной функцией srand().

Оценка проводилась в скользящем временном окне длительностью 40 000 отсчетов (минут). После каждой оценки окно сдвигалось на 180 отсчетов. Интервалы выстраивались с доверительной вероятностью 7 = 0.95.

В таблице приводятся результаты тестирования. Два последних столбца содержат процент покрытия доверительными интервалами истинных значений параметров A и B.

Волна Период, мин Покрытие A, % Покрытие B, %

Oi 1549.160455 98.2 98.3

Si 1440.0 98.6 98.7

M2 745.236078 97.9 97.6

S2 720.0 98.2 98.2

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

Заключение

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

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

[1] Гордеев Е. И., Салтыков В. А., Синицин В. И., Чебров В. Н. К вопросу о связи высокочастотного сейсмического шума с лунно-солнечными приливами // Докл. РАН. 1995. Т. 340, №3. С. 386-388.

[2] Салтыков В. А., Синицин В. И., Чебров В.Н. Вариации приливной компоненты высокочастотного сейсмического шума в результате изменения напряженного состояния среды // Вулканология и сейсмология. 1997. №4. С. 73-83.

[3] Салтыков В. А. О воздействии земных приливов на сейсмические процессы // Проблемы сейсмичности Дальнего Востока / Под ред. А. В. Викулина. Петропавловск-Камчатский, 2000. С. 12-21.

[4] Кролевец А. Н., ПАвлюков В. К. Инициирование приливного отклика импульсного электромагнитного излучения из литосферы процессами в очагах землетрясений. Петропавловск-Камчатский, 1999. (Препр. / Камчатский пединститут; №1(01)).

[5] Кролевец А. Н., ПАвлюков В. К. Приливной отклик импульсного электромагнитного излучения и краткосрочный прогноз сильных землетрясений // Проблемы сейсмичности Дальнего Востока / Под ред. А. В. Викулина. Петропавловск-Камчатский, 2000. С. 171-181.

[6] Мельхиор П. Земные приливы. М., 1968.

[7] Бендат Дж., Пирсол А. Прикладной анализ случайных данных. М., 1989.

[8] МАрпл-мл. С. Л. Цифровой спектральный анализ и его приложения. М., 1990.

[9] Серебренников М. Г., Первозванский А. А. Выявление скрытых периодичностей, М., 1965.

[10] ВодинчАр Г.М. Погрешности оценивания параметров приливного отклика в импульсном электромагнитном излучении Земли // Вычисл. технологии. 2001. Т. 6, №3. С. 3-6.

[11] Справочник по теории вероятностей и математической статистике / В. С. Королюк и др. М., 1985.

[12] Гантмахер Ф. Р. Теория матриц. М., 1953.

Поступила в редакцию 11 марта 2002 г.