CC BY
31
Вычислительные технологии
Том 7, № 5, 2002
МЕТОДИКА ОЦЕНИВАНИЯ НЕКОТОРЫХ ПАРАМЕТРОВ ИМПУЛЬСНОГО ЭЛЕКТРОМАГНИТНОГО ИЗЛУЧЕНИЯ ЗЕМЛИ
Г. М. ВодинчАР Камчатский государственный педагогический университет Петропавловск-Камчатский, Россия e-mail: g_vodin@kgpu.kamchatka.ru
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 г.
